Commits June 2012

commits@lists.einsteintoolkit.org

4 participants
73 discussions

New changesets in Carpet Mercurial
by Carpet mercurial repositories 20 Jun '12

20 Jun '12

changeset: 3584:b3da6ac7b72a user: Roland Haas <roland.haas(a)physics.gatech.edu> date: Sat Apr 28 14:46:40 2012 -0700 summary: CarpetLibs: shift vertex centered directions in restrict_3d_vc_rf2 changeset: 3585:a5f32a94867b tag: tip user: Roland Haas <roland.haas(a)physics.gatech.edu> date: Sat May 05 15:15:09 2012 -0700 summary: CarpetLib: add index checking on source array in mixed restriction operator

1 0

[svn:einsteintoolkit] GRHydro/trunk/src/ (Rev. 361)
by rhaas＠tapir.caltech.edu 20 Jun '12

20 Jun '12

User: rhaas Date: 2012/06/20 11:27 AM Modified: /trunk/src/ GRHydro_Eigenproblem.F90 Log: really improve handling of cs2 < 0 in the "standard" (non-MHD, non-hot EOS) eigenvalue routine. the previous patch with this same commit message (r359) has the wrong commit message. It should instead read: --8<-- Introduce real-valued atmo mask that can be interpolated. Synchronize real-valued atmo mask and update integer-valued atmo mask based on synchronized value. This fix ensures that the results are independent of the number of processors used. It also gets rid of some extra primitive syncs that are otherwise necessary when sync_conserved_only = yes. Does not cause any new test failures (Roland). Patch by Christian Reisswig. --8<-- Patch by Christian Ott. File Changes: Directory: /trunk/src/ ====================== File [modified]: GRHydro_Eigenproblem.F90 Delta lines: +5 -2 =================================================================== --- trunk/src/GRHydro_Eigenproblem.F90 2012-06-20 16:19:10 UTC (rev 360) +++ trunk/src/GRHydro_Eigenproblem.F90 2012-06-20 16:27:09 UTC (rev 361) @@ -86,13 +86,16 @@ (1.0d0 + eps + press/rho) if(cs2.lt.0.0d0) then + !$OMP CRITICAL if (abs(cs2) .gt. 1.0d-4) then write(warnline,'(a50,6g16.7)') 'abs(cs2), rho, dpdrho, press*dpdeps/rho**2, eps, press/rho: ', abs(cs2), rho, dpdrho, press * dpdeps / (rho**2), eps, press/rho call CCTK_WARN(1,warnline) - call CCTK_WARN(0,"cs2 < 0! Check speed of sound calculation!") + call CCTK_WARN(1,"cs2 < 0! Check speed of sound calculation!") + cs2 = 0.0d0 else - cs2 = 0 + cs2 = 0.0d0 endif + !$OMP END CRITICAL endif

1 0

[svn:einsteintoolkit] GRHydro/trunk/src/ (Rev. 360)
by rhaas＠tapir.caltech.edu 20 Jun '12

20 Jun '12

User: rhaas Date: 2012/06/20 11:19 AM Modified: /trunk/src/ GRHydro_ParamCheck.F90 Log: GRHydro: do not access whichpsidcspeed in ParamCheck unless evolve_MHD File Changes: Directory: /trunk/src/ ====================== File [modified]: GRHydro_ParamCheck.F90 Delta lines: +5 -1 =================================================================== --- trunk/src/GRHydro_ParamCheck.F90 2012-06-20 16:18:04 UTC (rev 359) +++ trunk/src/GRHydro_ParamCheck.F90 2012-06-20 16:19:10 UTC (rev 360) @@ -126,7 +126,7 @@ end if - if(clean_divergence.ne.0) then + if(evolve_MHD.ne.0 .and. clean_divergence.ne.0) then if (CCTK_EQUALS(psidcspeed,"char speed")) then whichpsidcspeed = 0 else if (CCTK_EQUALS(psidcspeed,"light speed")) then @@ -138,6 +138,10 @@ end if end if + if(clean_divergence.eq.0 .and. evolve_MHD.ne.0) then + call CCTK_WARN(1, "You activated divergence cleaning but do not evolve a magnetic field. At best, nothing should happen, be preapered for a segfault otherwise.") + end if + if (Tmunu_damping_radius_min .gt. Tmunu_damping_radius_max) then call CCTK_WARN(0, "Minimum damping radius is greater than maximum damping radius!") end if

1 0

[svn:einsteintoolkit] GRHydro/trunk/ (Rev. 359)
by rhaas＠tapir.caltech.edu 20 Jun '12

20 Jun '12

User: rhaas Date: 2012/06/20 11:18 AM Modified: /trunk/ interface.ccl, schedule.ccl /trunk/src/ GRHydro_Boundaries.F90, GRHydro_UpdateMask.F90 Log: * improve handling of cs2 < 0 in the "standard" (non-MHD, non-hot EOS) eigenvalue routine. Patch by Christian Ott. File Changes: Directory: /trunk/src/ ====================== File [modified]: GRHydro_Boundaries.F90 Delta lines: +6 -1 =================================================================== --- trunk/src/GRHydro_Boundaries.F90 2012-06-20 16:17:32 UTC (rev 358) +++ trunk/src/GRHydro_Boundaries.F90 2012-06-20 16:18:04 UTC (rev 359) @@ -99,6 +99,8 @@ !endif endif + call SetCartSymVN(ierr, cctkGH, sym, "GRHydro::atmosphere_mask_real") + sym(1) = -1 sym(2) = 1 sym(3) = 1 @@ -210,6 +212,8 @@ "GRHydro::tau", "Flat") ierr = ierr + Boundary_SelectGroupForBC(cctkGH, faces, GRHydro_stencil, -ione, & "GRHydro::scon", "Flat") + ierr = ierr + Boundary_SelectGroupForBC(cctkGH, faces, GRHydro_stencil, -ione, & + "GRHydro::GRHydro_atmosphere_mask_real", "Flat") if (sync_conserved_only .eq. 0) then ierr = ierr + Boundary_SelectGroupForBC(cctkGH, faces, GRHydro_stencil, -ione, & "HydroBase::w_lorentz", "Flat") @@ -282,7 +286,8 @@ "GRHydro::tau", "None") ierr = ierr + Boundary_SelectGroupForBC(cctkGH, faces, GRHydro_stencil, -ione, & "GRHydro::scon", "None") - + ierr = ierr + Boundary_SelectGroupForBC(cctkGH, faces, GRHydro_stencil, -ione, & + "GRHydro::GRHydro_atmosphere_mask_real", "None") if (sync_conserved_only .eq. 0) then ierr = ierr + Boundary_SelectGroupForBC(cctkGH, faces, GRHydro_stencil, -ione, & "HydroBase::w_lorentz", "None") File [modified]: GRHydro_UpdateMask.F90 Delta lines: +35 -2 =================================================================== --- trunk/src/GRHydro_UpdateMask.F90 2012-06-20 16:17:32 UTC (rev 358) +++ trunk/src/GRHydro_UpdateMask.F90 2012-06-20 16:18:04 UTC (rev 359) @@ -48,7 +48,9 @@ densrhs(i,j,k) = 0.0d0 srhs(i,j,k,:) = 0.0d0 taurhs(i,j,k) = 0.0d0 - atmosphere_mask(i,j,k) = 1 + ! Set real-valued mask! This will be sync'ed and right after syncing translated to + ! our standard integer based mask (so that atmosphere_mask is still valid!). + atmosphere_mask_real(i,j,k) = 1 end if end do end do @@ -67,7 +69,9 @@ densrhs(i,j,k) = 0.0d0 srhs(i,j,k,:) = 0.0d0 taurhs(i,j,k) = 0.0d0 - atmosphere_mask(i,j,k) = 1 + ! Set real-valued mask! This will be sync'ed and right after syncing translated to + ! our standard integer based mask (so that atmosphere_mask is still valid!). + atmosphere_mask_real(i,j,k) = 1 end if end do end do @@ -77,6 +81,33 @@ end subroutine GRHydroUpdateAtmosphereMask + +subroutine GRHydroPostSyncAtmosphereMask(CCTK_ARGUMENTS) + + implicit none + + DECLARE_CCTK_ARGUMENTS + DECLARE_CCTK_PARAMETERS + + CCTK_INT :: i,j,k + +!! This sets the integer atmo mask based on the real-valued (and sync'ed) atmo mask + + !$OMP PARALLEL DO PRIVATE(k,j,i) + do k = 1, cctk_lsh(3) + do j = 1, cctk_lsh(2) + do i = 1, cctk_lsh(1) + if ( atmosphere_mask_real(i,j,k) .ne. 0) then + atmosphere_mask(i,j,k) = 1 + end if + end do + end do + end do + !$OMP END PARALLEL DO + +end subroutine GRHydroPostSyncAtmosphereMask + + /*@@ @routine GRHydro_SetupMask @date Thu Jun 20 13:27:28 2002 @@ -116,6 +147,7 @@ endif atmosphere_mask = 0 + atmosphere_mask_real = 0 call CCTK_INFO("Setting up the atmosphere mask: all points are not_atmosphere") @@ -263,6 +295,7 @@ velz(i,j,k), eps(i,j,k), press(i,j,k), w_lorentz(i,j,k)) if (wk_atmosphere .eq. 0) then atmosphere_mask(i, j, k) = 0 + atmosphere_mask_real(i, j, k) = 0 end if endif Directory: /trunk/ ================== File [modified]: interface.ccl Delta lines: +9 -0 =================================================================== --- trunk/interface.ccl 2012-06-20 16:17:32 UTC (rev 358) +++ trunk/interface.ccl 2012-06-20 16:18:04 UTC (rev 359) @@ -424,6 +424,15 @@ atmosphere_mask } "Flags to say whether a point needs to be reset to the atmosphere" +# This real mask is set during UpdateAtmosphereMask and sync'ed afterwards (including possible interpatch interpolation) +# After syn'ing and before any con2prim fun, we set the integer mask above based on the real-valued mask. +# This ensures that any routine using the int mask is still correctly working. +real GRHydro_atmosphere_mask_real type = GF Timelevels = 1 tags='Prolongation="copy"' +{ + atmosphere_mask_real +} "Flags to say whether a point needs to be reset to the atmosphere. This is sync'ed (and possibly interpolated)!" + + int GRHydro_atmosphere_descriptors type=SCALAR { atmosphere_field_descriptor File [modified]: schedule.ccl Delta lines: +11 -7 =================================================================== --- trunk/schedule.ccl 2012-06-20 16:17:32 UTC (rev 358) +++ trunk/schedule.ccl 2012-06-20 16:18:04 UTC (rev 359) @@ -477,6 +477,7 @@ if (CCTK_EQUALS(evolution_method, "GRHydro")) { STORAGE:GRHydro_atmosphere_mask +STORAGE:GRHydro_atmosphere_mask_real STORAGE:GRHydro_atmosphere_descriptors schedule GRHydro_SetupMask AT CCTK_Initial BEFORE HydroBase_Initial @@ -1024,6 +1025,7 @@ SYNC: Y_e_con SYNC: lvel SYNC: lBvec + SYNC: GRHydro_atmosphere_mask_real # Must be synced otherwise we get inconsistent results when using multiple procs } "Select GRHydro boundary conditions" } else @@ -1039,6 +1041,7 @@ SYNC: psidc SYNC: GRHydro_cons_tracers SYNC: Y_e_con + SYNC: GRHydro_atmosphere_mask_real # Must be synced otherwise we get inconsistent results when using multiple procs } "Select GRHydro boundary conditions" # after a regrid Cactus relies on the boundary thorns to fill in the outer @@ -1050,11 +1053,6 @@ { } "Apply boundary conditions to primitives" - # This is necessary to provide initial guesses for Con2Prim after atmo reset (which happens only at last MoL post-step)! - schedule group GRHydro_PrimitiveBoundaries IN HydroBase_PostStep BEFORE HydroBase_Boundaries IF GRHydro::InLastMoLPostStep - { - } "Apply boundary conditions to primitives" - if(CCTK_IsImplementationActive("Coordinates")) { @@ -1160,14 +1158,20 @@ ### This is executed in the last MoL_PosStep ### ################################################################ +# First we need to convert real (and synchronized mask) to an integer mask! +schedule GRHydroPostSyncAtmosphereMask IN HydroBase_PostStep AFTER HydroBase_Boundaries +{ + LANG: Fortran +} "Set integer atmosphere mask from synchronized real atmosphere mask" + if (CCTK_Equals(Bvec_evolution_method,"GRHydro")) { - schedule GRHydro_AtmosphereResetM IN MoL_PostStep AFTER ADMBase_SetADMVars AFTER GRHydro_SetLastMoLPostStep BEFORE HydroBase_PostStep IF GRHydro::InLastMoLPostStep + schedule GRHydro_AtmosphereResetM IN HydroBase_PostStep BEFORE HydroBase_Con2Prim IF GRHydro::InLastMoLPostStep { LANG: Fortran } "Reset the atmosphere - MHD version" } else { - schedule GRHydro_AtmosphereReset IN MoL_PostStep AFTER ADMBase_SetADMVars BEFORE HydroBase_PostStep IF GRHydro::InLastMoLPostStep + schedule GRHydro_AtmosphereReset IN HydroBase_PostStep BEFORE HydroBase_Con2Prim IF GRHydro::InLastMoLPostStep { LANG: Fortran } "Reset the atmosphere"

1 0

[svn:einsteintoolkit] GRHydro/trunk/src/ (Rev. 358)
by rhaas＠tapir.caltech.edu 20 Jun '12

20 Jun '12

User: rhaas Date: 2012/06/20 11:17 AM Modified: /trunk/src/ GRHydro_HLLEM.F90 Log: Change the characteristic speed for psidc. It was set to +-c when it should be the maximum coming from the constraint that the normal vector to the characteristic hypersurface of a wave propagating along an arbitrary x direction be spacelike. This results into the following bound for the wavespeed: +- \alpha sqrt(\gamma^{xx}) - \beta^x) Patch by Bruno Mundim. File Changes: Directory: /trunk/src/ ====================== File [modified]: GRHydro_HLLEM.F90 Delta lines: +20 -5 =================================================================== --- trunk/src/GRHydro_HLLEM.F90 2012-06-20 16:16:52 UTC (rev 357) +++ trunk/src/GRHydro_HLLEM.F90 2012-06-20 16:17:32 UTC (rev 358) @@ -55,6 +55,7 @@ CCTK_REAL :: pressstarp,pressstarm,velxlowp,velxlowm,velylowp,velylowm,velzlowp,velzlowm CCTK_REAL :: psidcp, psidcm, psidcf, psidcdiff, psidcfp, psidcfm + CCTK_REAL :: charmax_dc, charmin_dc, charpm_dc CCTK_INT :: type_bits, trivial, not_trivial @@ -440,13 +441,27 @@ end if case(1) !!$ Wavespeeds for psidc are +/-c, not Fast Magnetosonic? - psidcf = 0.5d0 * (1.d0 * psidcfp - (-1.d0) * psidcfm + & - 1.d0 * (-1.d0) * psidcdiff) + !!$ psidcf = 0.5d0 * (1.d0 * psidcfp - (-1.d0) * psidcfm + & + !!$ 1.d0 * (-1.d0) * psidcdiff) + + !!$ The fastest speed for psidc comes from the condition + !!$ that the normal vector to the characteristic hypersurface + !!$ be spacelike (Eq. 60 of Anton et al.) + + charmax_dc = sqrt(usendh) - avg_beta/avg_alp + charmin_dc = -1.d0*sqrt(usendh) - avg_beta/avg_alp + charpm_dc = charmax_dc - charmin_dc + + psidcf = (charmax_dc * psidcfp - charmin_dc * psidcfm + & + charmax_dc * charmin_dc * psidcdiff) / charpm_dc + if(decouple_normal_Bfield .ne. 0) then ! same expression for HLLE and LLF !!$ B^i field decouples from the others and has same propagation - !!$ speed as divergence - f1(5+flux_direction) = 0.5d0 * (1.d0 * fplus(5+flux_direction) - (-1.d0) * fminus(5+flux_direction) + & - 1.d0 * (-1.d0) * qdiff(5+flux_direction)) + !!$ speed as divergence -null direction, + !!$ \pm sqrt(g^{xx}} - beta^x/alpha + f1(5+flux_direction) = (charmax_dc * fplus(5+flux_direction) & + - charmin_dc * fminus(5+flux_direction) + & + charmax_dc * charmin_dc * qdiff(5+flux_direction)) / charpm_dc end if case(2) charmax = setcharmax

1 0

[svn:einsteintoolkit] GRHydro/trunk/src/ (Rev. 357)
by rhaas＠tapir.caltech.edu 20 Jun '12

20 Jun '12

User: rhaas Date: 2012/06/20 11:16 AM Modified: /trunk/src/ GRHydro_Con2Prim.F90 Log: GRHydro: do not use integers as booleans gfortran (and the Fortran standard I guess) does not allow this. File Changes: Directory: /trunk/src/ ====================== File [modified]: GRHydro_Con2Prim.F90 Delta lines: +2 -2 =================================================================== --- trunk/src/GRHydro_Con2Prim.F90 2012-06-20 16:15:59 UTC (rev 356) +++ trunk/src/GRHydro_Con2Prim.F90 2012-06-20 16:16:52 UTC (rev 357) @@ -486,7 +486,7 @@ ! Check if Newton-Raphson resulted in something reasonable! - if (c2p_resort_to_bisection) then + if (c2p_resort_to_bisection.ne.0) then tmp = (utau + pnew + udens)**2 - s2 plow = max(pmin, sqrt(s2) - utau - udens) @@ -872,7 +872,7 @@ ! Check if Newton-Raphson resulted in something reasonable! - if (c2p_resort_to_bisection) then + if (c2p_resort_to_bisection.ne.0) then tmp = (utau + pnew + udens)**2 - s2 plow = max(pminl, sqrt(s2) - utau - udens)

1 0

[svn:einsteintoolkit] GRHydro/trunk/ (Rev. 356)
by rhaas＠tapir.caltech.edu 20 Jun '12

20 Jun '12

User: rhaas Date: 2012/06/20 11:16 AM Modified: /trunk/ param.ccl Log: * make parameter c2p_resort_to_bisection steerable Patch by Christian Ott. File Changes: Directory: /trunk/ ================== File [modified]: param.ccl Delta lines: +1 -1 =================================================================== --- trunk/param.ccl 2012-06-20 16:15:15 UTC (rev 355) +++ trunk/param.ccl 2012-06-20 16:15:59 UTC (rev 356) @@ -314,7 +314,7 @@ 0: :: "greater than zero" } 1.e-20 -boolean c2p_resort_to_bisection "If the pressure guess is unphysical, should we try with bisection (slower, but more robust)" +boolean c2p_resort_to_bisection "If the pressure guess is unphysical, should we try with bisection (slower, but more robust)" STEERABLE=ALWAYS { } no

1 0

[svn:einsteintoolkit] GRHydro/trunk/src/ (Rev. 355)
by rhaas＠tapir.caltech.edu 20 Jun '12

20 Jun '12

User: rhaas Date: 2012/06/20 11:15 AM Modified: /trunk/src/ GRHydro_Con2Prim.F90 Log: C2P_hot: Resort to bisection as well (if user requested so). Patch by Christian Reisswig. File Changes: Directory: /trunk/src/ ====================== File [modified]: GRHydro_Con2Prim.F90 Delta lines: +44 -3 =================================================================== --- trunk/src/GRHydro_Con2Prim.F90 2012-06-20 16:14:28 UTC (rev 354) +++ trunk/src/GRHydro_Con2Prim.F90 2012-06-20 16:15:15 UTC (rev 355) @@ -491,7 +491,7 @@ tmp = (utau + pnew + udens)**2 - s2 plow = max(pmin, sqrt(s2) - utau - udens) - if ((pnew .lt. plow .or. tmp .le. 0.0d0) .and. c2p_resort_to_bisection) then + if (pnew .lt. plow .or. tmp .le. 0.0d0) 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) @@ -624,7 +624,7 @@ CCTK_REAL GRHydro_C2P_failed CCTK_REAL udens, usx, usy, usz, utau, pold, pnew, epsold, epsnew, w2, & w2mhalf, temp1, drhobydpress, depsbydpress, dpressbydeps, dpressbydrho, pmin, epsmin - CCTK_REAL epsminl,pminl + CCTK_REAL epsminl,pminl,plow,tmp CCTK_REAL local_perc_ptol character(len=256) warnline @@ -672,6 +672,10 @@ call EOS_Omni_press(handle,keytemp,GRHydro_eos_rf_prec,n,& rho,epsilon,temp,ye,xpress,keyerr,anyerr) 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) ! error handling if(anyerr.ne.0) then !$OMP CRITICAL @@ -861,8 +865,45 @@ dpressbydeps*depsbydpress pold = pnew - pnew = max(pold - f/df, pminl) + ! Try to obtain new pressure via Newton-Raphson. + + pnew = pold - f/df + + ! Check if Newton-Raphson resulted in something reasonable! + + if (c2p_resort_to_bisection) then + + tmp = (utau + pnew + udens)**2 - s2 + plow = max(pminl, sqrt(s2) - utau - udens) + + if (pnew .lt. plow .or. tmp .le. 0.0d0) 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 + + end if + + else + + ! This is the standard algorithm without resorting to bisection. + + pnew = max(pminl, pnew) + tmp = (utau + pnew + udens)**2 - s2 + + endif + + ! Check if we are still in the physical range now. + ! If not, we set the C2P failed mask, and set pnew = pold (which will exit the loop). + + if ((tmp .le. 0.0d0) .and. GRHydro_C2P_failed .eq. 0) then + GRHydro_C2P_failed = 1 + pnew = pold + endif + !!$ Recalculate primitive variables and function rho = udens * sqrt( (utau + pnew + udens)**2 - s2)/(utau + pnew + udens)

1 0

[svn:einsteintoolkit] GRHydro/trunk/ (Rev. 354)
by rhaas＠tapir.caltech.edu 20 Jun '12

20 Jun '12

User: rhaas Date: 2012/06/20 11:14 AM Modified: /trunk/ param.ccl /trunk/doc/ documentation.tex /trunk/src/ GRHydro_Con2Prim.F90 Log: GRHydro: C2P: resort to bisection if Newton-Raphson fails. This is only activated if user requested this. Patch by Christian Reisswig. File Changes: Directory: /trunk/doc/ ====================== File [modified]: documentation.tex Delta lines: +1 -1 =================================================================== --- trunk/doc/documentation.tex 2012-06-06 14:37:40 UTC (rev 353) +++ trunk/doc/documentation.tex 2012-06-20 16:14:28 UTC (rev 354) @@ -1991,7 +1991,7 @@ encouragement that it was possible! Thirdly, the support of the Cactus team, especially Tom Goodale, -Gabrielle Allen and Thomas Radke made (and, especially Thomas, continue to make) life a +Gabrielle Allen and Thomas Radke made (and continues to make) life a lot easier. Finally, for their work in coding, ideas and suggestions, or just Directory: /trunk/src/ ====================== File [modified]: GRHydro_Con2Prim.F90 Delta lines: +50 -10 =================================================================== --- trunk/src/GRHydro_Con2Prim.F90 2012-06-06 14:37:40 UTC (rev 353) +++ trunk/src/GRHydro_Con2Prim.F90 2012-06-20 16:14:28 UTC (rev 354) @@ -319,8 +319,6 @@ endif #endif - - end do end do end do @@ -368,7 +366,7 @@ ! begin EOS Omni vars integer :: n,keytemp,anyerr,keyerr(1) - real*8 :: xpress,xtemp,xye + real*8 :: xpress,xtemp,xye,tmp,plow n=1;keytemp=0;anyerr=0;keyerr(1)=0 xpress=0.0d0;xtemp=0.0d0;xye=0.0d0 ! end EOS Omni vars @@ -387,6 +385,10 @@ call EOS_Omni_press(handle,keytemp,GRHydro_eos_rf_prec,n,& rho,epsilon,xtemp,xye,xpress,keyerr,anyerr) 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) !!$ Check that the variables have a chance of being physical @@ -398,10 +400,10 @@ else !$OMP CRITICAL !!!! call CCTK_WARN(GRHydro_NaN_verbose, "c2p failed and being told not to reset the pressure") + !call CCTK_WARN(0, "c2p failed and being told not to reset the pressure") GRHydro_C2P_failed = 1 !$OMP END CRITICAL endif - endif !!$ Calculate rho and epsilon @@ -477,14 +479,52 @@ pold = pnew - pnew = max(pold - f/df, pmin) + ! Try to obtain new pressure via Newton-Raphson. + + pnew = pold - f/df + + ! Check if Newton-Raphson resulted in something reasonable! + + if (c2p_resort_to_bisection) then + + tmp = (utau + pnew + udens)**2 - s2 + plow = max(pmin, sqrt(s2) - utau - udens) + + if ((pnew .lt. plow .or. tmp .le. 0.0d0) .and. c2p_resort_to_bisection) 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 + + end if + + else + + ! This is the standard algorithm without resorting to bisection. + + pnew = max(pmin, pnew) + tmp = (utau + pnew + udens)**2 - s2 + + endif + + ! Check if we are still in the physical range now. + ! If not, we set the C2P failed mask, and set pnew = pold (which will exit the loop). + + if ((tmp .le. 0.0d0) .and. GRHydro_C2P_failed .eq. 0) then + GRHydro_C2P_failed = 1 + pnew = pold + endif + + !!$ Recalculate primitive variables and function - - rho = udens * sqrt( (utau + pnew + udens)**2 - s2)/(utau + pnew + udens) - w_lorentz = (utau + pnew + udens) / sqrt( (utau + pnew + udens)**2 - & - s2) - epsilon = (sqrt( (utau + pnew + udens)**2 - s2) - pnew * w_lorentz - & + + rho = udens * sqrt( tmp)/(utau + pnew + udens) + + w_lorentz = (utau + pnew + udens) / sqrt( tmp) + epsilon = (sqrt( tmp) - pnew * w_lorentz - & udens)/udens call EOS_Omni_press(handle,keytemp,GRHydro_eos_rf_prec,n,& Directory: /trunk/ ================== File [modified]: param.ccl Delta lines: +4 -0 =================================================================== --- trunk/param.ccl 2012-06-06 14:37:40 UTC (rev 353) +++ trunk/param.ccl 2012-06-20 16:14:28 UTC (rev 354) @@ -314,6 +314,10 @@ 0: :: "greater than zero" } 1.e-20 +boolean c2p_resort_to_bisection "If the pressure guess is unphysical, should we try with bisection (slower, but more robust)" +{ +} no + real GRHydro_eps_min "Minimum value of specific internal energy - this is now only used in GRHydro_InitData's GRHydro_Only_Atmo routine" { 0: :: "Positive"

1 0

New changesets in Carpet Mercurial
by Carpet mercurial repositories 18 Jun '12

18 Jun '12

changeset: 3579:b806e7562361 user: Roland Haas <roland.haas(a)physics.gatech.edu> date: Sat Jun 02 09:23:05 2012 -0700 summary: CarpetEvolutionMask: make parameters steerable changeset: 3580:d2d5d16d4e1c user: Roland Haas <roland.haas(a)physics.gatech.edu> date: Sat Jun 02 09:33:37 2012 -0700 summary: CarpetEvolutionMask: use dd.ghost_widths rather than cctk_nghostzones changeset: 3581:f690c246e8d5 user: Roland Haas <roland.haas(a)physics.gatech.edu> date: Sat Jun 02 09:34:57 2012 -0700 summary: CarpetEvolutionMask: add parameter enlarge_evolved_region to fine tune changeset: 3582:e29b32263030 user: Roland Haas <roland.haas(a)physics.gatech.edu> date: Sat Jun 02 09:38:54 2012 -0700 summary: CarpetEvolutionMask: enlarge evolved region by one in vertex centering changeset: 3583:0616c103e7a0 tag: tip user: Roland Haas <roland.haas(a)physics.gatech.edu> date: Sat Jun 02 09:41:06 2012 -0700 summary: CarpetEvolutionMask: add TODOs

1 0

2026

2025

2024

2023

2022

2021

2020

2019

2018

2017

2016

2015

2014

2013

2012

2011

2010

Commits June 2012