Commits July 2012

commits@lists.einsteintoolkit.org

7 participants
70 discussions

[svn:einsteintoolkit] GRHydro/trunk/src/ (Rev. 402)
by rhaas＠tapir.caltech.edu 17 Jul '12

17 Jul '12

User: rhaas Date: 2012/07/17 12:08 PM Modified: /trunk/src/ GRHydro_Con2Prim.F90 Log: GRHydro: remove comparison of logical to .true. I appreciate the effort to make the code compile with gcc by no longer using integers in place of logicals. I really do. However it seems that F0RTRAN in its infinite wisdom requires you to write ... .and. mustbisect .eqv. .true. (note the all important 'v') when testing logicals. I have instead opted to remove the whole .eqv. term in favor of ... .and. mustbisect From: Roland Haas <roland.haas(a)physics.gatech.edu> File Changes: Directory: /trunk/src/ ====================== File [modified]: GRHydro_Con2Prim.F90 Delta lines: +2 -2 =================================================================== --- trunk/src/GRHydro_Con2Prim.F90 2012-07-17 17:08:33 UTC (rev 401) +++ trunk/src/GRHydro_Con2Prim.F90 2012-07-17 17:08:35 UTC (rev 402) @@ -556,7 +556,7 @@ tmp = (utau + pnew + udens)**2 - s2 plow = max(pmin, sqrt(s2) - utau - udens) - if (pnew .lt. plow .or. tmp .le. 0.0d0 .or. mustbisect .eq. .true.) then + if (pnew .lt. plow .or. tmp .le. 0.0d0 .or. mustbisect) then ! Ok, Newton-Raphson ended up finding something unphysical. ! Let's try to find our root via bisection (which converges slower but is more robust) @@ -1039,7 +1039,7 @@ tmp = (utau + pnew + udens)**2 - s2 plow = max(pminl, sqrt(s2) - utau - udens) - if (pnew .lt. plow .or. tmp .le. 0.0d0 .or. mustbisect .eq. .true.) then + if (pnew .lt. plow .or. tmp .le. 0.0d0 .or. mustbisect) then ! Ok, Newton-Raphson ended up finding something unphysical. ! Let's try to find our root via bisection (which converges slower but is more robust)

1 0

[svn:einsteintoolkit] GRHydro/trunk/src/ (Rev. 401)
by rhaas＠tapir.caltech.edu 17 Jul '12

17 Jul '12

User: rhaas Date: 2012/07/17 12:08 PM Modified: /trunk/src/ GRHydro_Con2Prim.F90 Log: GRHydro: Con2Prim improvement to bisection algorithm. This fixes some NANs associated with NS collapse test case when full GR is used. From: Christian Reisswig <reisswig(a)scriwalker.(none)> File Changes: Directory: /trunk/src/ ====================== File [modified]: GRHydro_Con2Prim.F90 Delta lines: +175 -8 =================================================================== --- trunk/src/GRHydro_Con2Prim.F90 2012-07-17 17:08:31 UTC (rev 400) +++ trunk/src/GRHydro_Con2Prim.F90 2012-07-17 17:08:33 UTC (rev 401) @@ -118,6 +118,21 @@ do j = 1, ny do i = 1, nx +#if 0 + !$OMP CRITICAL + if (rho(i,j,k).le.0.0d0) then + call CCTK_WARN(1,"rho less than zero at entry of Con2Prim") + write(warnline,"(A10,i5)") "reflevel: ",GRHydro_reflevel + call CCTK_WARN(1,warnline) + write(warnline,"(3i5,1P10E15.6)") i,j,k,x(i,j,k),y(i,j,k),z(i,j,k) + call CCTK_WARN(1,warnline) + write(warnline,"(1P10E15.6)") rho(i,j,k),dens(i,j,k),eps(i,j,k) + call CCTK_WARN(1,warnline) + call CCTK_WARN(0,"Aborting!!!") + endif + !$OMP END CRITICAL +#endif + !do not compute if in atmosphere or in excised region if ((atmosphere_mask(i,j,k) .ne. 0) .or. & (hydro_excision_mask(i,j,k) .ne. 0)) cycle @@ -320,6 +335,46 @@ endif #endif +!! Some debugging convenience output +# if 0 + !$OMP CRITICAL + if (dens(i,j,k).ne.dens(i,j,k)) then + call CCTK_WARN(1,"NAN at exit of Con2Prim") + write(warnline,"(A10,i5)") "reflevel: ",GRHydro_reflevel + call CCTK_WARN(1,warnline) + write(warnline,"(3i5,1P10E15.6)") i,j,k,x(i,j,k),y(i,j,k),z(i,j,k) + call CCTK_WARN(1,warnline) + write(warnline,"(1P10E15.6)") rho(i,j,k),dens(i,j,k),eps(i,j,k) + call CCTK_WARN(1,warnline) + call CCTK_WARN(0,"Aborting!!!") + endif + if (rho(i,j,k).ne.rho(i,j,k)) then + call CCTK_WARN(1,"NAN in rho at exit of Con2Prim") + write(warnline,"(A10,i5)") "reflevel: ",GRHydro_reflevel + call CCTK_WARN(1,warnline) + write(warnline,"(3i5,1P10E15.6)") i,j,k,x(i,j,k),y(i,j,k),z(i,j,k) + call CCTK_WARN(1,warnline) + write(warnline,"(1P10E15.6)") rho(i,j,k),dens(i,j,k),eps(i,j,k) + call CCTK_WARN(1,warnline) + call CCTK_WARN(0,"Aborting!!!") + endif + !$OMP END CRITICAL + + !$OMP CRITICAL + if (rho(i,j,k).le.0.0d0) then + call CCTK_WARN(1,"rho less than zero at exit of Con2Prim") + write(warnline,"(A10,i5)") "reflevel: ",GRHydro_reflevel + call CCTK_WARN(1,warnline) + write(warnline,"(3i5,1P10E15.6)") i,j,k,x(i,j,k),y(i,j,k),z(i,j,k) + call CCTK_WARN(1,warnline) + write(warnline,"(1P10E15.6)") rho(i,j,k),dens(i,j,k),eps(i,j,k) + call CCTK_WARN(1,warnline) + call CCTK_WARN(0,"Aborting!!!") + endif + !$OMP END CRITICAL +#endif + + end do end do end do @@ -363,7 +418,7 @@ CCTK_REAL udens, usx, usy, usz, utau, pold, pnew, epsold, epsnew, w2, & w2mhalf, temp1, drhobydpress, depsbydpress, dpressbydeps, dpressbydrho, pmin, epsmin character(len=200) warnline - logical epsnegative + logical epsnegative, mustbisect ! begin EOS Omni vars integer :: n,keytemp,anyerr,keyerr(1) @@ -385,12 +440,16 @@ !!$ Set initial guess for pressure: call EOS_Omni_press(handle,keytemp,GRHydro_eos_rf_prec,n,& rho,epsilon,xtemp,xye,xpress,keyerr,anyerr) - pold = max(1.d-10,xpress) + !pold = max(1.d-10,xpress) ! This is the lowest admissible pressure, otherwise we are off the physical regime ! Sometimes, we may end up with bad initial guesses. In that case we need to start off with something ! reasonable at least plow = max(pmin, sqrt(s2) - utau - udens) + + ! Start out with some reasonable initial guess + pold = max(plow+1.d-10,xpress) + mustbisect = .false. !!$ Check that the variables have a chance of being physical @@ -423,6 +482,11 @@ f = pold - xpress + if (f .ne. f) then + ! Ok, this yielded nonsense, let's enforce bisection! + mustbisect = .true. + endif + !!$Find the root count = 0 @@ -492,14 +556,15 @@ tmp = (utau + pnew + udens)**2 - s2 plow = max(pmin, sqrt(s2) - utau - udens) - if (pnew .lt. plow .or. tmp .le. 0.0d0) then + if (pnew .lt. plow .or. tmp .le. 0.0d0 .or. mustbisect .eq. .true.) then ! Ok, Newton-Raphson ended up finding something unphysical. ! Let's try to find our root via bisection (which converges slower but is more robust) pnew = (plow + pold) / 2 tmp = (utau + pnew + udens)**2 - s2 - + + mustbisect = .false. end if else @@ -524,6 +589,48 @@ rho = udens * sqrt( tmp)/(utau + pnew + udens) +!! Some debugging convenience output +#if 0 + if (rho .ne. rho) then + !$OMP CRITICAL + write(warnline,'(a38, i16)') 'NAN in rho! Root finding iteration: ', count + call CCTK_WARN(1,warnline) + write(warnline,'(a28,i2)') 'on carpet reflevel: ',GRHydro_reflevel + call CCTK_WARN(1,warnline) + write(warnline,'(a20,3g16.7)') 'xyz location: ',x,y,z + call CCTK_WARN(1,warnline) + write(warnline,'(a20,g16.7)') 'radius: ',r + call CCTK_WARN(1,warnline) + write(warnline,'(a20,g16.7)') 'udens: ', udens + call CCTK_WARN(1,warnline) + write(warnline,'(a20,g16.7)') 'dens: ', dens + call CCTK_WARN(1,warnline) + !write(warnline,'(a20,g16.7)') 'rho-old: ', rhoold + !call CCTK_WARN(1,warnline) + write(warnline,'(a20,g16.7)') 'pnew: ', pnew + call CCTK_WARN(1,warnline) + write(warnline,'(a20,g16.7)') 'pold: ', pold + call CCTK_WARN(1,warnline) + write(warnline,'(a20,g16.7)') 'xpress: ', xpress + call CCTK_WARN(1,warnline) + write(warnline,'(a20,g16.7)') 'f: ', f + call CCTK_WARN(1,warnline) + write(warnline,'(a20,g16.7)') 'df: ', df + call CCTK_WARN(1,warnline) + write(warnline,'(a20,g16.7)') 'f/df: ', f/df + call CCTK_WARN(1,warnline) + write(warnline,'(a30,g16.7)') '(utau + pnew + udens)**2 - s2: ', (utau + pnew + udens)**2 - s2 + call CCTK_WARN(1,warnline) + write(warnline,'(a30,g16.7)') '(utau + pnew + udens): ', (utau + pnew + udens) + call CCTK_WARN(1,warnline) + write(warnline,'(a20,g16.7)') 'drhobydpress: ', drhobydpress + call CCTK_WARN(1,warnline) + write(warnline,'(a20,g16.7)') 'depsbydpress: ', depsbydpress + call CCTK_WARN(0,warnline) + !$OMP END CRITICAL + endif +#endif + w_lorentz = (utau + pnew + udens) / sqrt( tmp) epsilon = (sqrt( tmp) - pnew * w_lorentz - & udens)/udens @@ -533,6 +640,11 @@ f = pnew - xpress + if (f .ne. f) then + ! Ok, this yielded nonsense, let's enforce bisection! + mustbisect = .true. + endif + enddo !!$ Polish the root @@ -604,6 +716,47 @@ epsnegative = .true. endif +#if 0 + if (rho .le. 0.0d0) then + !$OMP CRITICAL + write(warnline,'(a38, i16)') 'Epsilon negative!' + call CCTK_WARN(1,warnline) + write(warnline,'(a28,i2)') 'on carpet reflevel: ',GRHydro_reflevel + call CCTK_WARN(1,warnline) + write(warnline,'(a20,3g16.7)') 'xyz location: ',x,y,z + call CCTK_WARN(1,warnline) + write(warnline,'(a20,g16.7)') 'radius: ',r + call CCTK_WARN(1,warnline) + write(warnline,'(a20,g16.7)') 'udens: ', udens + call CCTK_WARN(1,warnline) + write(warnline,'(a20,g16.7)') 'dens: ', dens + call CCTK_WARN(1,warnline) + !write(warnline,'(a20,g16.7)') 'rho-old: ', rhoold + !call CCTK_WARN(1,warnline) + write(warnline,'(a20,g16.7)') 'pnew: ', pnew + call CCTK_WARN(1,warnline) + write(warnline,'(a20,g16.7)') 'pold: ', pold + call CCTK_WARN(1,warnline) + write(warnline,'(a20,g16.7)') 'xpress: ', xpress + call CCTK_WARN(1,warnline) + write(warnline,'(a20,g16.7)') 'f: ', f + call CCTK_WARN(1,warnline) + write(warnline,'(a20,g16.7)') 'df: ', df + call CCTK_WARN(1,warnline) + write(warnline,'(a20,g16.7)') 'f/df: ', f/df + call CCTK_WARN(1,warnline) + write(warnline,'(a30,g16.7)') '(utau + pnew + udens)**2 - s2: ', (utau + pnew + udens)**2 - s2 + call CCTK_WARN(1,warnline) + write(warnline,'(a30,g16.7)') '(utau + pnew + udens): ', (utau + pnew + udens) + call CCTK_WARN(1,warnline) + write(warnline,'(a20,g16.7)') 'drhobydpress: ', drhobydpress + call CCTK_WARN(1,warnline) + write(warnline,'(a20,g16.7)') 'depsbydpress: ', depsbydpress + call CCTK_WARN(0,warnline) + !$OMP END CRITICAL + endif +#endif + end subroutine Con2Prim_pt subroutine Con2Prim_pt_hot(cctk_iteration, ii,jj,kk,handle, dens, & @@ -629,7 +782,7 @@ CCTK_REAL local_perc_ptol character(len=256) warnline - logical epsnegative + logical epsnegative, mustbisect integer :: failwarnmode integer :: failinfomode @@ -643,6 +796,8 @@ nf=0;nfudgemax=30 ! end EOS Omni vars + mustbisect = .false. + failinfomode = 1 failwarnmode = 0 @@ -672,11 +827,13 @@ temp0 = temp call EOS_Omni_press(handle,keytemp,GRHydro_eos_rf_prec,n,& rho,epsilon,temp,ye,xpress,keyerr,anyerr) - pold = max(1.d-15,xpress) + !pold = max(1.d-15,xpress) ! This is the lowest admissible pressure, otherwise we are off the physical regime ! Sometimes, we may end up with bad initial guesses. In that case we need to start off with something ! reasonable at least plow = max(pminl, sqrt(s2) - utau - udens) + ! Start out with some reasonable initial guess + pold = max(plow+1.d-10,xpress) ! error handling if(anyerr.ne.0) then !$OMP CRITICAL @@ -836,6 +993,10 @@ endif f = pold - xpress + if (f .ne. f) then + ! Ok, this yielded nonsense, let's enforce bisection! + mustbisect = .true. + endif !!$Find the root @@ -878,14 +1039,15 @@ tmp = (utau + pnew + udens)**2 - s2 plow = max(pminl, sqrt(s2) - utau - udens) - if (pnew .lt. plow .or. tmp .le. 0.0d0) then + if (pnew .lt. plow .or. tmp .le. 0.0d0 .or. mustbisect .eq. .true.) then ! Ok, Newton-Raphson ended up finding something unphysical. ! Let's try to find our root via bisection (which converges slower but is more robust) pnew = (plow + pold) / 2 tmp = (utau + pnew + udens)**2 - s2 - + + mustbisect = .false. end if else @@ -985,6 +1147,11 @@ f = pnew - xpress + if (f .ne. f) then + ! Ok, this yielded nonsense, let's enforce bisection! + mustbisect = .true. + endif + enddo

1 0

[svn:einsteintoolkit] GRHydro/trunk/ (Rev. 400)
by rhaas＠tapir.caltech.edu 17 Jul '12

17 Jul '12

User: rhaas Date: 2012/07/17 12:08 PM Modified: /trunk/ schedule.ccl Log: GRHydro: schedule AtmosphereReset before PrimitiveInitialGuessesBoundaries which used to be called PrimitiveBoundaries From: Roland Haas <roland.haas(a)physics.gatech.edu> File Changes: Directory: /trunk/ ================== File [modified]: schedule.ccl Delta lines: +1 -1 =================================================================== --- trunk/schedule.ccl 2012-07-17 17:08:30 UTC (rev 399) +++ trunk/schedule.ccl 2012-07-17 17:08:31 UTC (rev 400) @@ -1009,7 +1009,7 @@ # Synchronize the atmosphere mask before doing other boundaries! -schedule group GRHydro_AtmosphereMaskBoundaries IN HydroBase_PostStep BEFORE (HydroBase_Boundaries,GRHydro_PrimitiveBoundaries) +schedule group GRHydro_AtmosphereMaskBoundaries IN HydroBase_PostStep BEFORE (HydroBase_Boundaries,GRHydro_PrimitiveInitialGuessesBoundaries) { } "Apply boundary conditions to primitives"

1 0

[svn:einsteintoolkit] GRHydro/trunk/ (Rev. 399)
by rhaas＠tapir.caltech.edu 17 Jul '12

17 Jul '12

User: rhaas Date: 2012/07/17 12:08 PM Modified: /trunk/ schedule.ccl Log: * schedule GRHydroCopyIntegerMask in POST_RECOVER_VARIABLES before MoL_PostStep rather than HydroBase_PostStep, since the latter is scheduled IN MoL_PostStep -> rely on correct scheduling rather than on luck to get the correct atmosphere map after recovery. From: Christian Ott <cott(a)bethe.tapir.caltech.edu> File Changes: Directory: /trunk/ ================== File [modified]: schedule.ccl Delta lines: +1 -1 =================================================================== --- trunk/schedule.ccl 2012-07-17 17:08:28 UTC (rev 398) +++ trunk/schedule.ccl 2012-07-17 17:08:30 UTC (rev 399) @@ -487,7 +487,7 @@ # alternatively drop support for wk_atmosphere and schedule GRHydro_SetupMask # only in CCTK_BaseGrid -schedule GRHydroCopyIntegerMask AT CCTK_POST_RECOVER_VARIABLES BEFORE HydroBase_PostStep +schedule GRHydroCopyIntegerMask AT CCTK_POST_RECOVER_VARIABLES BEFORE MoL_PostStep { LANG: Fortran } "Initialize the real valued atmosphere mask after checkpoint recovery"

1 0

[svn:einsteintoolkit] GRHydro/trunk/ (Rev. 398)
by rhaas＠tapir.caltech.edu 17 Jul '12

17 Jul '12

User: rhaas Date: 2012/07/17 12:08 PM Modified: /trunk/ param.ccl /trunk/src/ GRHydro_PPM.F90 Log: GRHydro: ePPM: Allow for tolerance in vel^2 limiting. From: Christian Reisswig <reisswig(a)scriwalker.(none)> File Changes: Directory: /trunk/src/ ====================== File [modified]: GRHydro_PPM.F90 Delta lines: +1 -1 =================================================================== --- trunk/src/GRHydro_PPM.F90 2012-07-17 17:08:26 UTC (rev 397) +++ trunk/src/GRHydro_PPM.F90 2012-07-17 17:08:28 UTC (rev 398) @@ -791,7 +791,7 @@ !! Check if velocity is below speed of light. If not, reduce reconstruction to first order! #define VELCHECK(vxminus,vx,vxplus,vyminus,vy,vyplus,vzminus,vz,vzplus) \ - if (VEL2M(vxminus(i),vyminus(i),vzminus(i)) .ge. 1.0d0 .or. VEL2P(vxplus(i),vyplus(i),vzplus(i)) .ge. 1.0d0) then &&\ + if (VEL2M(vxminus(i),vyminus(i),vzminus(i)) .ge. (1.0d0-enhanced_ppm_C4) .or. VEL2P(vxplus(i),vyplus(i),vzplus(i)) .ge. (1.0d0-enhanced_ppm_C4)) then &&\ vxminus(i) = vx(i) &&\ vxplus(i) = vx(i) &&\ vyminus(i) = vy(i) &&\ Directory: /trunk/ ================== File [modified]: param.ccl Delta lines: +5 -0 =================================================================== --- trunk/param.ccl 2012-07-17 17:08:26 UTC (rev 397) +++ trunk/param.ccl 2012-07-17 17:08:28 UTC (rev 398) @@ -218,7 +218,12 @@ 0:* :: "must be greater than 0." } 0.1 +real enhanced_ppm_C4 "Parameter for enhancecd ppm limiter. This is for the additional velocity limiter" STEERABLE=ALWAYS +{ + 0:1 :: "must be greater than or equal to 0, and not larger than 1" +} 0.0 + int eno_order "The order of accuracy of the ENO reconstruction" { 1:* :: "Anything above 1, but above 5 is pointless"

1 0

[svn:einsteintoolkit] GRHydro/trunk/src/ (Rev. 397)
by rhaas＠tapir.caltech.edu 17 Jul '12

17 Jul '12

User: rhaas Date: 2012/07/17 12:08 PM Modified: /trunk/src/ GRHydro_PPM.F90 Log: GRHydro: ePPM: Corrected a logic bug in vel2 constraint. Also, use averaged metric. From: Christian Reisswig <reisswig(a)scriwalker.(none)> File Changes: Directory: /trunk/src/ ====================== File [modified]: GRHydro_PPM.F90 Delta lines: +20 -5 =================================================================== --- trunk/src/GRHydro_PPM.F90 2012-07-17 17:08:25 UTC (rev 396) +++ trunk/src/GRHydro_PPM.F90 2012-07-17 17:08:26 UTC (rev 397) @@ -767,16 +767,31 @@ endif +! Reconstruct quantity on plus face +#define AVGP(a) \ + 0.5*(a(i)+a(i+1)) -!! Compute scalar product v^i v_i -#define VEL2(dvelx,dvely,dvelz) \ - (gxx(i)*dvelx*dvelx + gyy(i)*dvely*dvely + gzz(i) \ - *dvelz*dvelz + 2*gxy(i)*dvelx*dvely + 2*gxz(i)*dvelx *dvelz + 2*gyz(i) \ +! Reconstruct quantity on minus face +#define AVGM(a) \ + 0.5*(a(i)+a(i-1)) + + +!! Compute scalar product v^i v_i on plus face +#define VEL2P(dvelx,dvely,dvelz) \ + (AVGP(gxx)*dvelx*dvelx + AVGP(gyy)*dvely*dvely + AVGP(gzz) \ + *dvelz*dvelz + 2*AVGP(gxy)*dvelx*dvely + 2*AVGP(gxz)*dvelx *dvelz + 2*AVGP(gyz) \ *dvely*dvelz) +!! Compute scalar product v^i v_i on minus face +#define VEL2M(dvelx,dvely,dvelz) \ + (AVGM(gxx)*dvelx*dvelx + AVGM(gyy)*dvely*dvely + AVGM(gzz) \ + *dvelz*dvelz + 2*AVGM(gxy)*dvelx*dvely + 2*AVGM(gxz)*dvelx *dvelz + 2*AVGM(gyz) \ + *dvely*dvelz) + + !! Check if velocity is below speed of light. If not, reduce reconstruction to first order! #define VELCHECK(vxminus,vx,vxplus,vyminus,vy,vyplus,vzminus,vz,vzplus) \ - if (1.0d0 .ge. VEL2(vxminus(i),vyminus(i),vzminus(i)) .or. 1.0d0 .ge. VEL2(vxplus(i),vyplus(i),vzplus(i))) then &&\ + if (VEL2M(vxminus(i),vyminus(i),vzminus(i)) .ge. 1.0d0 .or. VEL2P(vxplus(i),vyplus(i),vzplus(i)) .ge. 1.0d0) then &&\ vxminus(i) = vx(i) &&\ vxplus(i) = vx(i) &&\ vyminus(i) = vy(i) &&\

1 0

[svn:einsteintoolkit] GRHydro/trunk/src/ (Rev. 396)
by rhaas＠tapir.caltech.edu 17 Jul '12

17 Jul '12

User: rhaas Date: 2012/07/17 12:08 PM Modified: /trunk/src/ GRHydro_PPM.F90 Log: GRHydro: ePPM: Constrain reconstructed vel to below speed of light! From: Christian Reisswig <reisswig(a)scriwalker.(none)> File Changes: Directory: /trunk/src/ ====================== File [modified]: GRHydro_PPM.F90 Delta lines: +23 -0 =================================================================== --- trunk/src/GRHydro_PPM.F90 2012-07-17 17:08:22 UTC (rev 395) +++ trunk/src/GRHydro_PPM.F90 2012-07-17 17:08:25 UTC (rev 396) @@ -224,6 +224,7 @@ endif + if (PPM3) then !! We initialize "plus" \equiv a_j+1/2 with (16) via APPROX_AT_CELL_INTERFACE, @@ -767,6 +768,24 @@ +!! Compute scalar product v^i v_i +#define VEL2(dvelx,dvely,dvelz) \ + (gxx(i)*dvelx*dvelx + gyy(i)*dvely*dvely + gzz(i) \ + *dvelz*dvelz + 2*gxy(i)*dvelx*dvely + 2*gxz(i)*dvelx *dvelz + 2*gyz(i) \ + *dvely*dvelz) + +!! Check if velocity is below speed of light. If not, reduce reconstruction to first order! +#define VELCHECK(vxminus,vx,vxplus,vyminus,vy,vyplus,vzminus,vz,vzplus) \ + if (1.0d0 .ge. VEL2(vxminus(i),vyminus(i),vzminus(i)) .or. 1.0d0 .ge. VEL2(vxplus(i),vyplus(i),vzplus(i))) then &&\ + vxminus(i) = vx(i) &&\ + vxplus(i) = vx(i) &&\ + vyminus(i) = vy(i) &&\ + vyplus(i) = vy(i) &&\ + vzminus(i) = vz(i) &&\ + vzplus(i) = vz(i) &&\ + endif + + if (use_enhanced_ppm .eq. 1) then !! Constrain parabolic profiles, PPM 2011/2008 if (PPM3) then @@ -775,6 +794,8 @@ MON_WITH_LOCAL_EXTREMUM(velxminus,velx,velxplus, enhanced_ppm_C2) MON_WITH_LOCAL_EXTREMUM(velyminus,vely,velyplus, enhanced_ppm_C2) MON_WITH_LOCAL_EXTREMUM(velzminus,velz,velzplus, enhanced_ppm_C2) + ! Ensure that velocity is below speed of light! + VELCHECK(velxminus,velx,velxplus,velyminus,vely,velyplus,velzminus,velz,velzplus) end do if (poly .eq. 0) then do i = 3, nx - 2 @@ -788,6 +809,8 @@ MON_WITH_LOCAL_EXTREMUM_STENCIL4(velxminus,velx,velxplus, enhanced_ppm_C2, enhanced_ppm_C3) MON_WITH_LOCAL_EXTREMUM_STENCIL4(velyminus,vely,velyplus, enhanced_ppm_C2, enhanced_ppm_C3) MON_WITH_LOCAL_EXTREMUM_STENCIL4(velzminus,velz,velzplus, enhanced_ppm_C2, enhanced_ppm_C3) + ! Ensure that velocity is below speed of light! + VELCHECK(velxminus,velx,velxplus,velyminus,vely,velyplus,velzminus,velz,velzplus) end do if (poly .eq. 0) then do i = 4, nx - 3

1 0

[svn:einsteintoolkit] GRHydro/trunk/src/ (Rev. 395)
by rhaas＠tapir.caltech.edu 17 Jul '12

17 Jul '12

User: rhaas Date: 2012/07/17 12:08 PM Modified: /trunk/src/ GRHydro_Prim2Con.F90, GRHydro_UpdateMask.F90 Log: * ensure consistency in calls to prim2con_hot (fixes issue with my previous patch) From: Christian Ott <cott(a)bethe.tapir.caltech.edu> File Changes: Directory: /trunk/src/ ====================== File [modified]: GRHydro_Prim2Con.F90 Delta lines: +5 -4 =================================================================== --- trunk/src/GRHydro_Prim2Con.F90 2012-07-17 17:08:20 UTC (rev 394) +++ trunk/src/GRHydro_Prim2Con.F90 2012-07-17 17:08:22 UTC (rev 395) @@ -305,8 +305,8 @@ ! This is a way of recovering even on finer refinement levels: ! Use the average temperature at the interface instead of the ! reconstructed specific internal energy. + if(GRHydro_eos_hot_prim2con_warn.ne.0) then !$OMP CRITICAL - if(GRHydro_eos_hot_prim2con_warn.ne.0) then call CCTK_WARN(1,"EOS error in prim2con_hot: NOW using averaged temp!") write(warnline,"(i8,3i5,1P10E15.6)") cctk_iteration,ii,jj,kk,x,y,z,r call CCTK_WARN(1,warnline) @@ -316,8 +316,8 @@ call CCTK_WARN(1,warnline) write(warnline,"(A10,i5)") "reflevel: ", GRHydro_reflevel call CCTK_WARN(1,warnline) + !$OMP END CRITICAL endif - !$OMP END CRITICAL keytemp=1 temp = temp0 call EOS_Omni_press(handle,keytemp,GRHydro_eos_rf_prec,n,& @@ -441,8 +441,9 @@ det = SPATIAL_DETERMINANT(g11(i,j,k),g12(i,j,k),g13(i,j,k), \ g22(i,j,k),g23(i,j,k),g33(i,j,k)) - call prim2con_hot(GRHydro_eos_handle,GRHydro_reflevel,i,j,k,& - x(i,j,k),y(i,j,k),z(i,j,k),g11(i,j,k),& + call prim2con_hot(GRHydro_eos_handle,GRHydro_reflevel,cctk_iteration,& + i,j,k,& + x(i,j,k),y(i,j,k),z(i,j,k),r(i,j,k),g11(i,j,k),& g12(i,j,k),g13(i,j,k),& g22(i,j,k),g23(i,j,k),g33(i,j,k),& det, dens(i,j,k),sx(i,j,k),sy(i,j,k),sz(i,j,k),& File [modified]: GRHydro_UpdateMask.F90 Delta lines: +8 -4 =================================================================== --- trunk/src/GRHydro_UpdateMask.F90 2012-07-17 17:08:20 UTC (rev 394) +++ trunk/src/GRHydro_UpdateMask.F90 2012-07-17 17:08:22 UTC (rev 395) @@ -303,7 +303,8 @@ press(i,j,k),keyerr,anyerr) call prim2con_hot(GRHydro_eos_handle, GRHydro_reflevel,& - i,j,k,x(i,j,k),y(i,j,k),z(i,j,k),g11(i,j,k),g12(i,j,k),& + cctk_iteration,i,j,k,x(i,j,k),y(i,j,k),z(i,j,k),r(i,j,k),& + g11(i,j,k),g12(i,j,k),& g13(i,j,k),g22(i,j,k),g23(i,j,k),g33(i,j,k), & det,dens(i,j,k),scon(i,j,k,1), scon(i,j,k,2), scon(i,j,k,3), & tau(i,j,k), rho(i,j,k), velx(i,j,k), vely(i,j,k), & @@ -465,7 +466,8 @@ press(i,j,k),keyerr,anyerr) call prim2con_hot(GRHydro_eos_handle, GRHydro_reflevel,& - i,j,k,x(i,j,k),y(i,j,k),z(i,j,k),g11(i,j,k),g12(i,j,k),& + cctk_iteration,i,j,k,x(i,j,k),y(i,j,k),z(i,j,k),r(i,j,k),& + g11(i,j,k),g12(i,j,k),& g13(i,j,k),g22(i,j,k),g23(i,j,k),g33(i,j,k), & det,dens(i,j,k),scon(i,j,k,1), scon(i,j,k,2), scon(i,j,k,3), & tau(i,j,k), rho(i,j,k), velx(i,j,k), vely(i,j,k), & @@ -505,7 +507,8 @@ press_p(i,j,k),keyerr,anyerr) call prim2con_hot(GRHydro_eos_handle, GRHydro_reflevel,& - i,j,k,x(i,j,k),y(i,j,k),z(i,j,k),g11_p(i,j,k),g12_p(i,j,k),& + cctk_iteration,i,j,k,x(i,j,k),y(i,j,k),z(i,j,k),r(i,j,k),& + g11_p(i,j,k),g12_p(i,j,k),& g13_p(i,j,k),g22_p(i,j,k),g23_p(i,j,k),g33_p(i,j,k), & det,dens_p(i,j,k),scon_p(i,j,k,1), scon_p(i,j,k,2), scon_p(i,j,k,3), & tau_p(i,j,k), rho_p(i,j,k), velx_p(i,j,k), vely_p(i,j,k), & @@ -545,7 +548,8 @@ rho_p_p(i,j,k),eps_p_p(i,j,k),temperature_p_p(i,j,k),y_e_p_p(i,j,k),& press_p_p(i,j,k),keyerr,anyerr) call prim2con_hot(GRHydro_eos_handle, GRHydro_reflevel,& - i,j,k,x(i,j,k),y(i,j,k),z(i,j,k),g11_p_p(i,j,k),g12_p_p(i,j,k),& + cctk_iteration,i,j,k,x(i,j,k),y(i,j,k),z(i,j,k),r(i,j,k),& + g11_p_p(i,j,k),g12_p_p(i,j,k),& g13_p_p(i,j,k),g22_p_p(i,j,k),g23_p_p(i,j,k),g33_p_p(i,j,k), & det,dens_p_p(i,j,k),scon_p_p(i,j,k,1), scon_p_p(i,j,k,2), scon_p_p(i,j,k,3), & tau_p_p(i,j,k), rho_p_p(i,j,k), velx_p_p(i,j,k), vely_p_p(i,j,k), &

1 0

[svn:einsteintoolkit] GRHydro/trunk/ (Rev. 394)
by rhaas＠tapir.caltech.edu 17 Jul '12

17 Jul '12

User: rhaas Date: 2012/07/17 12:08 PM Modified: /trunk/ param.ccl Log: * add parameter used in Prim2Con_hot From: Christian Ott <cott(a)bethe.tapir.caltech.edu> File Changes: Directory: /trunk/ ================== File [modified]: param.ccl Delta lines: +4 -0 =================================================================== --- trunk/param.ccl 2012-07-17 17:08:18 UTC (rev 393) +++ trunk/param.ccl 2012-07-17 17:08:20 UTC (rev 394) @@ -295,6 +295,10 @@ { } "no" +boolean GRHydro_eos_hot_prim2con_warn "Warn about temperature workaround in prim2con" STEERABLE=ALWAYS +{ +} "yes" + boolean GRHydro_c2p_reset_eps_tau_hot_eos "As a last resort, reset tau" STEERABLE=ALWAYS { } "no"

1 0

[svn:einsteintoolkit] GRHydro/trunk/src/ (Rev. 393)
by rhaas＠tapir.caltech.edu 17 Jul '12

17 Jul '12

User: rhaas Date: 2012/07/17 12:08 PM Modified: /trunk/src/ GRHydro_Prim2Con.F90 Log: * improve treatment of hot EOS warning From: Christian Ott <cott(a)bethe.tapir.caltech.edu> File Changes: Directory: /trunk/src/ ====================== File [modified]: GRHydro_Prim2Con.F90 Delta lines: +24 -18 =================================================================== --- trunk/src/GRHydro_Prim2Con.F90 2012-07-11 15:26:05 UTC (rev 392) +++ trunk/src/GRHydro_Prim2Con.F90 2012-07-17 17:08:18 UTC (rev 393) @@ -153,7 +153,9 @@ xtemp(1) = 0.5d0*(temperature(i,j,k) + & temperature(i-xoffset,j-yoffset,k-zoffset)) call prim2con_hot(GRHydro_eos_handle, GRHydro_reflevel,& - i,j,k,x(i,j,k),y(i,j,k),z(i,j,k),g11l,g12l,g13l,g22l,& + cctk_iteration,i,j,k,x(i,j,k),y(i,j,k),z(i,j,k),& + r(i,j,k),& + g11l,g12l,g13l,g22l,& g23l,g33l, & avg_detl,densminus(i,j,k),sxminus(i,j,k),& syminus(i,j,k),szminus(i,j,k),tauminus(i,j,k),& @@ -170,8 +172,10 @@ xtemp(1) = 0.5d0*(temperature(i,j,k) + & temperature(i+xoffset,j+yoffset,k+zoffset)) call prim2con_hot(GRHydro_eos_handle, GRHydro_reflevel, & - i,j,k,x(i,j,k),y(i,j,k),z(i,j,k),g11r,g12r,g13r,& - g22r,g23r,g33r, & + cctk_iteration,i,j,k,x(i,j,k),y(i,j,k),z(i,j,k),& + r(i,j,k),& + g11r,g12r,g13r,& + g22r,g23r,g33r,& avg_detr, densplus(i,j,k),sxplus(i,j,k),& syplus(i,j,k),szplus(i,j ,k),tauplus(i,j,k),& rhoplus(i,j,k),velxplus(i,j,k),velyplus(i,j,k),& @@ -241,8 +245,8 @@ end subroutine prim2con -subroutine prim2con_hot(handle, GRHydro_reflevel, ii, jj, kk, & - x, y, z, gxx, gxy, gxz, gyy, gyz, gzz, det, ddens, & +subroutine prim2con_hot(handle, GRHydro_reflevel, cctk_iteration, ii, jj, kk, & + x, y, z, r, gxx, gxy, gxz, gyy, gyz, gzz, det, ddens, & dsx, dsy, dsz, dtau , drho, dvelx, dvely, dvelz, deps, dpress, w, & temp,ye) @@ -254,8 +258,8 @@ CCTK_REAL :: gxx, gxy, gxz, gyy, gyz, gzz, det CCTK_REAL :: ddens, dsx, dsy, dsz, dtau, drho(1), dvelx, dvely, dvelz,& deps(1), dpress(1), w, vlowx, vlowy, vlowz - CCTK_REAL :: temp(1),ye(1), x, y, z - CCTK_INT :: handle, GRHydro_reflevel, ii, jj, kk + CCTK_REAL :: temp(1),ye(1), x, y, z, r + CCTK_INT :: handle, GRHydro_reflevel, cctk_iteration, ii, jj, kk CCTK_REAL :: sdet, h character(len=512) warnline @@ -287,7 +291,7 @@ if(anyerr.ne.0) then !$OMP CRITICAL call CCTK_WARN(1,"EOS error in prim2con_hot: lev 2") - write(warnline,"(3i5,1P10E15.6)") ii,jj,kk,x,y,z + write(warnline,"(i8,3i5,1P10E15.6)") cctk_iteration,ii,jj,kk,x,y,z,r call CCTK_WARN(1,warnline) write(warnline,"(1P10E15.6)") drho,deps,temp,ye call CCTK_WARN(1,warnline) @@ -302,15 +306,17 @@ ! Use the average temperature at the interface instead of the ! reconstructed specific internal energy. !$OMP CRITICAL - call CCTK_WARN(1,"EOS error in prim2con_hot: NOW using averaged temp!") - write(warnline,"(3i5,1P10E15.6)") ii,jj,kk,x,y,z - call CCTK_WARN(1,warnline) - write(warnline,"(1P10E15.6)") drho,deps,temp0,ye - call CCTK_WARN(1,warnline) - write(warnline,"(A7,i8)") "code: ",keyerr(1) - call CCTK_WARN(1,warnline) - write(warnline,"(A10,i5)") "reflevel: ", GRHydro_reflevel - call CCTK_WARN(1,warnline) + if(GRHydro_eos_hot_prim2con_warn.ne.0) then + call CCTK_WARN(1,"EOS error in prim2con_hot: NOW using averaged temp!") + write(warnline,"(i8,3i5,1P10E15.6)") cctk_iteration,ii,jj,kk,x,y,z,r + call CCTK_WARN(1,warnline) + write(warnline,"(1P10E15.6)") drho,deps,temp0,ye + call CCTK_WARN(1,warnline) + write(warnline,"(A7,i8)") "code: ",keyerr(1) + call CCTK_WARN(1,warnline) + write(warnline,"(A10,i5)") "reflevel: ", GRHydro_reflevel + call CCTK_WARN(1,warnline) + endif !$OMP END CRITICAL keytemp=1 temp = temp0 @@ -320,7 +326,7 @@ if(anyerr.ne.0) then !$OMP CRITICAL call CCTK_WARN(1,"EOS error in prim2con_hot") - write(warnline,"(3i5,1P10E15.6)") ii,jj,kk,x,y,z + write(warnline,"(i8,3i5,1P10E15.6)") cctk_iteration,ii,jj,kk,x,y,z,r call CCTK_WARN(1,warnline) write(warnline,"(1P10E15.6)") drho,deps,temp,ye call CCTK_WARN(1,warnline)

1 0

2026

2025

2024

2023

2022

2021

2020

2019

2018

2017

2016

2015

2014

2013

2012

2011

2010

Commits July 2012