diff options
| author | knarf <knarf@c83d129a-5a75-4d5a-9c4d-ed3a5855bf45> | 2013-03-04 17:01:33 +0000 |
|---|---|---|
| committer | knarf <knarf@c83d129a-5a75-4d5a-9c4d-ed3a5855bf45> | 2013-03-04 17:01:33 +0000 |
| commit | 27532e3585845c715469a88714fe6ab0d0ec4ebc (patch) | |
| tree | 708c6a121385d491016c80f2f57fccfca8092abb | |
| parent | 718f43b6350665df6df9c58ab2c278f9919e7014 (diff) | |
convert tabs to (8) spaces. No other changes.
git-svn-id: http://svn.einsteintoolkit.org/cactus/EinsteinEvolve/GRHydro/trunk@485 c83d129a-5a75-4d5a-9c4d-ed3a5855bf45
| -rw-r--r-- | src/GRHydro_Boundaries.F90 | 170 | ||||
| -rw-r--r-- | src/GRHydro_Con2Prim.F90 | 112 | ||||
| -rw-r--r-- | src/GRHydro_Con2PrimM_polytype_pt.c | 124 | ||||
| -rw-r--r-- | src/GRHydro_Con2PrimM_pt.c | 100 | ||||
| -rw-r--r-- | src/GRHydro_Con2PrimM_pt_EOSOmni.c | 146 | ||||
| -rw-r--r-- | src/GRHydro_Con2PrimM_pt_EOSOmni.c.hopefullyfixed.c | 156 | ||||
| -rw-r--r-- | src/GRHydro_EoSChangeGamma.F90 | 2 | ||||
| -rw-r--r-- | src/GRHydro_InterfacesAM.h | 10 | ||||
| -rw-r--r-- | src/GRHydro_InterfacesM.h | 10 | ||||
| -rw-r--r-- | src/GRHydro_Macros.h | 2 | ||||
| -rw-r--r-- | src/GRHydro_PPM.F90 | 1646 | ||||
| -rw-r--r-- | src/GRHydro_Tmunu.F90 | 8 | ||||
| -rw-r--r-- | src/GRHydro_TransformTensorBasis.c | 26 | ||||
| -rw-r--r-- | src/GRHydro_UpdateMask.F90 | 16 | ||||
| -rw-r--r-- | src/GRHydro_WENOReconstruct.F90 | 12 |
15 files changed, 1270 insertions, 1270 deletions
diff --git a/src/GRHydro_Boundaries.F90 b/src/GRHydro_Boundaries.F90 index 99a8e89..e3be957 100644 --- a/src/GRHydro_Boundaries.F90 +++ b/src/GRHydro_Boundaries.F90 @@ -90,9 +90,9 @@ subroutine GRHydro_InitSymBound(CCTK_ARGUMENTS) !if (sync_conserved_only .eq. 0) then if (general_coordinates .ne. 0) then - call SetCartSymVN(ierr, cctkGH, sym, "GRHydro::lvel[0]") + call SetCartSymVN(ierr, cctkGH, sym, "GRHydro::lvel[0]") else - call SetCartSymVN(ierr, cctkGH, sym, "HydroBase::vel[0]") + call SetCartSymVN(ierr, cctkGH, sym, "HydroBase::vel[0]") endif !endif call SetCartSymVN(ierr, cctkGH, sym, "GRHydro::scon[0]") @@ -109,15 +109,15 @@ subroutine GRHydro_InitSymBound(CCTK_ARGUMENTS) !if (sync_conserved_only .eq. 0) then if (general_coordinates .ne. 0) then - call SetCartSymVN(ierr, cctkGH, sym, "GRHydro::lvel[1]") + call SetCartSymVN(ierr, cctkGH, sym, "GRHydro::lvel[1]") else - call SetCartSymVN(ierr, cctkGH, sym, "HydroBase::vel[1]") + call SetCartSymVN(ierr, cctkGH, sym, "HydroBase::vel[1]") endif !endif call SetCartSymVN(ierr, cctkGH, sym, "GRHydro::scon[1]") if(evolve_mhd.ne.0) then !if (sync_conserved_only .eq. 0) then - call SetCartSymVN(ierr, cctkGH, sym, "HydroBase::Bvec[1]") + call SetCartSymVN(ierr, cctkGH, sym, "HydroBase::Bvec[1]") !endif call SetCartSymVN(ierr, cctkGH, sym, "GRHydro::Bcons[1]") endif @@ -128,9 +128,9 @@ subroutine GRHydro_InitSymBound(CCTK_ARGUMENTS) !if (sync_conserved_only .eq. 0) then if (general_coordinates .ne. 0) then - call SetCartSymVN(ierr, cctkGH, sym, "GRHydro::lvel[2]") + call SetCartSymVN(ierr, cctkGH, sym, "GRHydro::lvel[2]") else - call SetCartSymVN(ierr, cctkGH, sym, "HydroBase::vel[2]") + call SetCartSymVN(ierr, cctkGH, sym, "HydroBase::vel[2]") endif !endif call SetCartSymVN(ierr, cctkGH, sym, "GRHydro::scon[2]") @@ -196,30 +196,30 @@ subroutine GRHydro_Boundaries(CCTK_ARGUMENTS) "GRHydro::scon", "Flat") if (sync_conserved_only .eq. 0) then ierr = ierr + Boundary_SelectGroupForBC(cctkGH, faces, GRHydro_stencil, -ione, & - "HydroBase::w_lorentz", "Flat") + "HydroBase::w_lorentz", "Flat") ierr = ierr + Boundary_SelectGroupForBC(cctkGH, faces, GRHydro_stencil, -ione, & - "HydroBase::rho", "Flat") + "HydroBase::rho", "Flat") ierr = ierr + Boundary_SelectGroupForBC(cctkGH, faces, GRHydro_stencil, -ione, & - "HydroBase::press", "Flat") + "HydroBase::press", "Flat") ierr = ierr + Boundary_SelectGroupForBC(cctkGH, faces, GRHydro_stencil, -ione, & - "HydroBase::eps", "Flat") + "HydroBase::eps", "Flat") if (general_coordinates .ne. 0) then - ierr = ierr + Boundary_SelectGroupForBC(cctkGH, faces, GRHydro_stencil, -ione, & - "GRHydro::lvel", "Flat") + ierr = ierr + Boundary_SelectGroupForBC(cctkGH, faces, GRHydro_stencil, -ione, & + "GRHydro::lvel", "Flat") else - ierr = ierr + Boundary_SelectGroupForBC(cctkGH, faces, GRHydro_stencil, -ione, & - "HydroBase::vel", "Flat") + ierr = ierr + Boundary_SelectGroupForBC(cctkGH, faces, GRHydro_stencil, -ione, & + "HydroBase::vel", "Flat") endif endif if(evolve_mhd.ne.0) then if (sync_conserved_only .eq. 0) then - if (general_coordinates .ne. 0) then - ierr = ierr + Boundary_SelectGroupForBC(cctkGH, faces, GRHydro_stencil, -ione, & - "GRHydro::lBvec", "Flat") - else - ierr = ierr + Boundary_SelectGroupForBC(cctkGH, faces, GRHydro_stencil, -ione, & - "HydroBase::Bvec", "Flat") - endif + if (general_coordinates .ne. 0) then + ierr = ierr + Boundary_SelectGroupForBC(cctkGH, faces, GRHydro_stencil, -ione, & + "GRHydro::lBvec", "Flat") + else + ierr = ierr + Boundary_SelectGroupForBC(cctkGH, faces, GRHydro_stencil, -ione, & + "HydroBase::Bvec", "Flat") + endif endif ierr = ierr + Boundary_SelectGroupForBC(cctkGH, faces, GRHydro_stencil, -ione, & "GRHydro::Bcons", "Flat") @@ -243,8 +243,8 @@ subroutine GRHydro_Boundaries(CCTK_ARGUMENTS) if(evolve_tracer.ne.0) then if (sync_conserved_only .eq. 0) then - ierr = ierr + Boundary_SelectGroupForBC(cctkGH, faces, GRHydro_stencil, -ione, & - "GRHydro::GRHydro_tracers", "Flat") + ierr = ierr + Boundary_SelectGroupForBC(cctkGH, faces, GRHydro_stencil, -ione, & + "GRHydro::GRHydro_tracers", "Flat") endif ierr = ierr + Boundary_SelectGroupForBC(cctkGH, faces, GRHydro_stencil, -ione, & "GRHydro::GRHydro_cons_tracers", "Flat") @@ -252,8 +252,8 @@ subroutine GRHydro_Boundaries(CCTK_ARGUMENTS) if(evolve_y_e.ne.0) then if (sync_conserved_only .eq. 0) then - ierr = ierr + Boundary_SelectGroupForBC(cctkGH, faces, GRHydro_stencil, -ione, & - "HydroBase::Y_e", "Flat") + ierr = ierr + Boundary_SelectGroupForBC(cctkGH, faces, GRHydro_stencil, -ione, & + "HydroBase::Y_e", "Flat") endif ierr = ierr + Boundary_SelectGroupForBC(cctkGH, faces, GRHydro_stencil, -ione, & "GRHydro::Y_e_con", "Flat") @@ -261,10 +261,10 @@ subroutine GRHydro_Boundaries(CCTK_ARGUMENTS) if(evolve_temper.ne.0) then if (sync_conserved_only .eq. 0) then - ierr = ierr + Boundary_SelectGroupForBC(cctkGH, faces, GRHydro_stencil, -ione, & - "HydroBase::temperature", "Flat") - ierr = ierr + Boundary_SelectGroupForBC(cctkGH, faces, GRHydro_stencil, -ione, & - "HydroBase::entropy", "Flat") + ierr = ierr + Boundary_SelectGroupForBC(cctkGH, faces, GRHydro_stencil, -ione, & + "HydroBase::temperature", "Flat") + ierr = ierr + Boundary_SelectGroupForBC(cctkGH, faces, GRHydro_stencil, -ione, & + "HydroBase::entropy", "Flat") endif endif @@ -281,31 +281,31 @@ subroutine GRHydro_Boundaries(CCTK_ARGUMENTS) if (sync_conserved_only .eq. 0) then ierr = ierr + Boundary_SelectGroupForBC(cctkGH, faces, GRHydro_stencil, -ione, & - "HydroBase::w_lorentz", "None") + "HydroBase::w_lorentz", "None") ierr = ierr + Boundary_SelectGroupForBC(cctkGH, faces, GRHydro_stencil, -ione, & - "HydroBase::rho", "None") + "HydroBase::rho", "None") ierr = ierr + Boundary_SelectGroupForBC(cctkGH, faces, GRHydro_stencil, -ione, & - "HydroBase::press", "None") + "HydroBase::press", "None") ierr = ierr + Boundary_SelectGroupForBC(cctkGH, faces, GRHydro_stencil, -ione, & - "HydroBase::eps", "None") + "HydroBase::eps", "None") if (general_coordinates .ne. 0) then - ierr = ierr + Boundary_SelectGroupForBC(cctkGH, faces, GRHydro_stencil, -ione, & - "GRHydro::lvel", "None") + ierr = ierr + Boundary_SelectGroupForBC(cctkGH, faces, GRHydro_stencil, -ione, & + "GRHydro::lvel", "None") else - ierr = ierr + Boundary_SelectGroupForBC(cctkGH, faces, GRHydro_stencil, -ione, & - "HydroBase::vel", "None") + ierr = ierr + Boundary_SelectGroupForBC(cctkGH, faces, GRHydro_stencil, -ione, & + "HydroBase::vel", "None") endif endif if(evolve_mhd.ne.0) then if (sync_conserved_only .eq. 0) then - if (general_coordinates .ne. 0) then - ierr = ierr + Boundary_SelectGroupForBC(cctkGH, faces, GRHydro_stencil, -ione, & - "GRHydro::lBvec", "None") - else - ierr = ierr + Boundary_SelectGroupForBC(cctkGH, faces, GRHydro_stencil, -ione, & - "HydroBase::Bvec", "None") - endif + if (general_coordinates .ne. 0) then + ierr = ierr + Boundary_SelectGroupForBC(cctkGH, faces, GRHydro_stencil, -ione, & + "GRHydro::lBvec", "None") + else + ierr = ierr + Boundary_SelectGroupForBC(cctkGH, faces, GRHydro_stencil, -ione, & + "HydroBase::Bvec", "None") + endif endif ierr = ierr + Boundary_SelectGroupForBC(cctkGH, faces, GRHydro_stencil, -ione, & "GRHydro::Bcons", "None") @@ -329,8 +329,8 @@ subroutine GRHydro_Boundaries(CCTK_ARGUMENTS) if(evolve_tracer.ne.0) then if (sync_conserved_only .eq. 0) then - ierr = ierr + Boundary_SelectGroupForBC(cctkGH, faces, GRHydro_stencil, -ione, & - "GRHydro::GRHydro_tracers", "None") + ierr = ierr + Boundary_SelectGroupForBC(cctkGH, faces, GRHydro_stencil, -ione, & + "GRHydro::GRHydro_tracers", "None") endif ierr = ierr + Boundary_SelectGroupForBC(cctkGH, faces, GRHydro_stencil, -ione, & "GRHydro::GRHydro_cons_tracers", "None") @@ -338,8 +338,8 @@ subroutine GRHydro_Boundaries(CCTK_ARGUMENTS) if(evolve_y_e.ne.0) then if (sync_conserved_only .eq. 0) then - ierr = ierr + Boundary_SelectGroupForBC(cctkGH, faces, GRHydro_stencil, -ione, & - "HydroBase::Y_e", "None") + ierr = ierr + Boundary_SelectGroupForBC(cctkGH, faces, GRHydro_stencil, -ione, & + "HydroBase::Y_e", "None") endif ierr = ierr + Boundary_SelectGroupForBC(cctkGH, faces, GRHydro_stencil, -ione, & "GRHydro::Y_e_con", "None") @@ -347,10 +347,10 @@ subroutine GRHydro_Boundaries(CCTK_ARGUMENTS) if(evolve_temper.ne.0) then if (sync_conserved_only .eq. 0) then - ierr = ierr + Boundary_SelectGroupForBC(cctkGH, faces, GRHydro_stencil, -ione, & - "HydroBase::temperature", "None") - ierr = ierr + Boundary_SelectGroupForBC(cctkGH, faces, GRHydro_stencil, -ione, & - "HydroBase::entropy", "None") + ierr = ierr + Boundary_SelectGroupForBC(cctkGH, faces, GRHydro_stencil, -ione, & + "HydroBase::temperature", "None") + ierr = ierr + Boundary_SelectGroupForBC(cctkGH, faces, GRHydro_stencil, -ione, & + "HydroBase::entropy", "None") endif endif @@ -412,13 +412,13 @@ subroutine GRHydro_SelectPrimitiveInitialGuessesBoundaries(CCTK_ARGUMENTS) "HydroBase::vel", "Flat") endif if(evolve_mhd.ne.0) then - if (general_coordinates .ne. 0) then - ierr = ierr + Boundary_SelectGroupForBC(cctkGH, faces, GRHydro_stencil, -ione, & - "GRHydro::lBvec", "Flat") - else - ierr = ierr + Boundary_SelectGroupForBC(cctkGH, faces, GRHydro_stencil, -ione, & - "HydroBase::Bvec", "Flat") - endif + if (general_coordinates .ne. 0) then + ierr = ierr + Boundary_SelectGroupForBC(cctkGH, faces, GRHydro_stencil, -ione, & + "GRHydro::lBvec", "Flat") + else + ierr = ierr + Boundary_SelectGroupForBC(cctkGH, faces, GRHydro_stencil, -ione, & + "HydroBase::Bvec", "Flat") + endif endif !if(evolve_tracer.ne.0) then @@ -459,13 +459,13 @@ subroutine GRHydro_SelectPrimitiveInitialGuessesBoundaries(CCTK_ARGUMENTS) endif if(evolve_mhd.ne.0) then - if (general_coordinates .ne. 0) then - ierr = ierr + Boundary_SelectGroupForBC(cctkGH, faces, GRHydro_stencil, -ione, & - "GRHydro::lBvec", "None") - else - ierr = ierr + Boundary_SelectGroupForBC(cctkGH, faces, GRHydro_stencil, -ione, & - "HydroBase::Bvec", "None") - endif + if (general_coordinates .ne. 0) then + ierr = ierr + Boundary_SelectGroupForBC(cctkGH, faces, GRHydro_stencil, -ione, & + "GRHydro::lBvec", "None") + else + ierr = ierr + Boundary_SelectGroupForBC(cctkGH, faces, GRHydro_stencil, -ione, & + "HydroBase::Bvec", "None") + endif endif @@ -548,23 +548,23 @@ subroutine GRHydro_SelectPrimitiveBoundaries(CCTK_ARGUMENTS) "HydroBase::vel", "Flat") endif if(evolve_mhd.ne.0) then - if (general_coordinates .ne. 0) then - ierr = ierr + Boundary_SelectGroupForBC(cctkGH, faces, GRHydro_stencil, -ione, & - "GRHydro::lBvec", "Flat") - else - ierr = ierr + Boundary_SelectGroupForBC(cctkGH, faces, GRHydro_stencil, -ione, & - "HydroBase::Bvec", "Flat") - endif + if (general_coordinates .ne. 0) then + ierr = ierr + Boundary_SelectGroupForBC(cctkGH, faces, GRHydro_stencil, -ione, & + "GRHydro::lBvec", "Flat") + else + ierr = ierr + Boundary_SelectGroupForBC(cctkGH, faces, GRHydro_stencil, -ione, & + "HydroBase::Bvec", "Flat") + endif endif if(evolve_tracer.ne.0) then - ierr = ierr + Boundary_SelectGroupForBC(cctkGH, faces, GRHydro_stencil, -ione, & - "GRHydro::GRHydro_tracers", "Flat") + ierr = ierr + Boundary_SelectGroupForBC(cctkGH, faces, GRHydro_stencil, -ione, & + "GRHydro::GRHydro_tracers", "Flat") endif if(evolve_y_e.ne.0) then - ierr = ierr + Boundary_SelectGroupForBC(cctkGH, faces, GRHydro_stencil, -ione, & - "HydroBase::Y_e", "Flat") + ierr = ierr + Boundary_SelectGroupForBC(cctkGH, faces, GRHydro_stencil, -ione, & + "HydroBase::Y_e", "Flat") endif if(evolve_temper.ne.0) then @@ -595,13 +595,13 @@ subroutine GRHydro_SelectPrimitiveBoundaries(CCTK_ARGUMENTS) endif if(evolve_mhd.ne.0) then - if (general_coordinates .ne. 0) then - ierr = ierr + Boundary_SelectGroupForBC(cctkGH, faces, GRHydro_stencil, -ione, & - "GRHydro::lBvec", "None") - else - ierr = ierr + Boundary_SelectGroupForBC(cctkGH, faces, GRHydro_stencil, -ione, & - "HydroBase::Bvec", "None") - endif + if (general_coordinates .ne. 0) then + ierr = ierr + Boundary_SelectGroupForBC(cctkGH, faces, GRHydro_stencil, -ione, & + "GRHydro::lBvec", "None") + else + ierr = ierr + Boundary_SelectGroupForBC(cctkGH, faces, GRHydro_stencil, -ione, & + "HydroBase::Bvec", "None") + endif endif diff --git a/src/GRHydro_Con2Prim.F90 b/src/GRHydro_Con2Prim.F90 index 5af9518..1387c65 100644 --- a/src/GRHydro_Con2Prim.F90 +++ b/src/GRHydro_Con2Prim.F90 @@ -120,16 +120,16 @@ subroutine Conservative2Primitive(CCTK_ARGUMENTS) #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 + 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 @@ -341,41 +341,41 @@ subroutine Conservative2Primitive(CCTK_ARGUMENTS) !! 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 + !$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 + 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 @@ -562,13 +562,13 @@ subroutine Con2Prim_pt(handle, dens, sx, sy, sz, tau, rho, velx, vely, & 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) + ! 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. + pnew = (plow + pold) / 2 + tmp = (utau + pnew + udens)**2 - s2 + + mustbisect = .false. end if else @@ -1054,13 +1054,13 @@ subroutine Con2Prim_pt_hot(cctk_iteration, ii,jj,kk,handle, dens, & 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) + ! 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. + pnew = (plow + pold) / 2 + tmp = (utau + pnew + udens)**2 - s2 + + mustbisect = .false. end if else diff --git a/src/GRHydro_Con2PrimM_polytype_pt.c b/src/GRHydro_Con2PrimM_polytype_pt.c index 97d7648..4931415 100644 --- a/src/GRHydro_Con2PrimM_polytype_pt.c +++ b/src/GRHydro_Con2PrimM_polytype_pt.c @@ -60,9 +60,9 @@ #define NEWT_TOL2 (1.0e-15) /* TOL of new 1D^*_{v^2} gnr2 method */ #define MIN_NEWT_TOL2 (1.0e-10) /* TOL of new 1D^*_{v^2} gnr2 method */ -#define W_TOO_BIG (1.e20) /* \gamma^2 (\rho_0 + u + p) is assumed +#define W_TOO_BIG (1.e20) /* \gamma^2 (\rho_0 + u + p) is assumed to always be smaller than this. This - is used to detect solver failures */ + is used to detect solver failures */ #define FAIL_VAL (1.e30) /* Generic value to which we set variables when a problem arises */ @@ -80,23 +80,23 @@ struct LocGlob { // Declarations: static CCTK_INT oned_newton_raphson( CCTK_REAL x[], CCTK_INT n, CCTK_REAL gammaeos, struct LocGlob *lgp, - void (*funcd) (CCTK_REAL [], CCTK_REAL [], - CCTK_REAL [], CCTK_REAL [][1], - CCTK_REAL *, CCTK_REAL *, CCTK_INT, CCTK_REAL, - struct LocGlob *lgp) ); + void (*funcd) (CCTK_REAL [], CCTK_REAL [], + CCTK_REAL [], CCTK_REAL [][1], + CCTK_REAL *, CCTK_REAL *, CCTK_INT, CCTK_REAL, + struct LocGlob *lgp) ); static void func_W( CCTK_REAL [], CCTK_REAL [], CCTK_REAL [], CCTK_REAL [][1], - CCTK_REAL *f, CCTK_REAL *df, CCTK_INT n, CCTK_REAL gammaeos, - struct LocGlob * lgp); + CCTK_REAL *f, CCTK_REAL *df, CCTK_INT n, CCTK_REAL gammaeos, + struct LocGlob * lgp); static CCTK_INT gnr2( CCTK_REAL x[], CCTK_INT n, CCTK_REAL gammaeos, struct LocGlob *lgp, - void (*funcd) (CCTK_REAL [], CCTK_REAL [], CCTK_REAL [], - CCTK_REAL [][1],CCTK_REAL *,CCTK_REAL *,CCTK_INT, CCTK_REAL, - struct LocGlob *lgp) ); + void (*funcd) (CCTK_REAL [], CCTK_REAL [], CCTK_REAL [], + CCTK_REAL [][1],CCTK_REAL *,CCTK_REAL *,CCTK_INT, CCTK_REAL, + struct LocGlob *lgp) ); static void func_rho(CCTK_REAL x[], CCTK_REAL dx[], CCTK_REAL resid[], - CCTK_REAL jac[][1], CCTK_REAL *f, CCTK_REAL *df, CCTK_INT n, CCTK_REAL gammaeos, - struct LocGlob *lgp); + CCTK_REAL jac[][1], CCTK_REAL *f, CCTK_REAL *df, CCTK_INT n, CCTK_REAL gammaeos, + struct LocGlob *lgp); void CCTK_FCALL CCTK_FNAME(GRHydro_Con2PrimM_Polytype_pt) ( CCTK_INT *handle, CCTK_REAL *gamma_eos, @@ -171,11 +171,11 @@ RETURN VALUE: of retval = (i*100 + j) where given tolerances; 2 -> Newton-Raphson procedure encountered a numerical divergence (occurrence of "nan" or "+/-inf" ; - + j = 0 -> success 1 -> failure: some sort of failure in Newton-Raphson; 2 -> failure: unphysical vsq = v^2 value at initial guess; - 3 -> failure: W<0 or W>W_TOO_BIG + 3 -> failure: W<0 or W>W_TOO_BIG 4 -> failure: v^2 > 1 ( used to be 5 -> failure: rho,uu <= 0 but now sets epsnegative to non-zero ) @@ -226,37 +226,37 @@ void CCTK_FCALL CCTK_FNAME(GRHydro_Con2PrimM_Polytype_pt) ( #if(DEBUG_CON2PRIMM) fprintf(stdout," *dens = %26.16e \n", *dens_in ); - fprintf(stdout," *sx = %26.16e \n", *sx_in ); - fprintf(stdout," *sy = %26.16e \n", *sy_in ); - fprintf(stdout," *sz = %26.16e \n", *sz_in ); - fprintf(stdout," *Sc = %26.16e \n", *sc_in ); - fprintf(stdout," *Bconsx = %26.16e \n", *Bconsx_in ); - fprintf(stdout," *Bconsy = %26.16e \n", *Bconsy_in ); - fprintf(stdout," *Bconsz = %26.16e \n", *Bconsz_in ); - fprintf(stdout," *rho = %26.16e \n", *rho ); - fprintf(stdout," *velx = %26.16e \n", *velx ); - fprintf(stdout," *vely = %26.16e \n", *vely ); - fprintf(stdout," *velz = %26.16e \n", *velz ); + fprintf(stdout," *sx = %26.16e \n", *sx_in ); + fprintf(stdout," *sy = %26.16e \n", *sy_in ); + fprintf(stdout," *sz = %26.16e \n", *sz_in ); + fprintf(stdout," *Sc = %26.16e \n", *sc_in ); + fprintf(stdout," *Bconsx = %26.16e \n", *Bconsx_in ); + fprintf(stdout," *Bconsy = %26.16e \n", *Bconsy_in ); + fprintf(stdout," *Bconsz = %26.16e \n", *Bconsz_in ); + fprintf(stdout," *rho = %26.16e \n", *rho ); + fprintf(stdout," *velx = %26.16e \n", *velx ); + fprintf(stdout," *vely = %26.16e \n", *vely ); + fprintf(stdout," *velz = %26.16e \n", *velz ); fprintf(stdout," *epsilon = %26.16e \n", *epsilon ); fprintf(stdout," *pressure = %26.16e \n", *pressure ); - fprintf(stdout," *Bx = %26.16e \n", *Bx ); - fprintf(stdout," *By = %26.16e \n", *By ); - fprintf(stdout," *Bz = %26.16e \n", *Bz ); - fprintf(stdout," *bsq = %26.16e \n", *bsq ); + fprintf(stdout," *Bx = %26.16e \n", *Bx ); + fprintf(stdout," *By = %26.16e \n", *By ); + fprintf(stdout," *Bz = %26.16e \n", *Bz ); + fprintf(stdout," *bsq = %26.16e \n", *bsq ); fprintf(stdout," *w_lorentz = %26.16e \n", *w_lorentz ); - fprintf(stdout," *gxx = %26.16e \n", *gxx ); - fprintf(stdout," *gxy = %26.16e \n", *gxy ); - fprintf(stdout," *gxz = %26.16e \n", *gxz ); - fprintf(stdout," *gyy = %26.16e \n", *gyy ); - fprintf(stdout," *gyz = %26.16e \n", *gyz ); - fprintf(stdout," *gzz = %26.16e \n", *gzz ); - fprintf(stdout," *uxx = %26.16e \n", *uxx ); - fprintf(stdout," *uxy = %26.16e \n", *uxy ); - fprintf(stdout," *uxz = %26.16e \n", *uxz ); - fprintf(stdout," *uyy = %26.16e \n", *uyy ); - fprintf(stdout," *uyz = %26.16e \n", *uyz ); - fprintf(stdout," *uzz = %26.16e \n", *uzz ); - fprintf(stdout," *det = %26.16e \n", *det ); + fprintf(stdout," *gxx = %26.16e \n", *gxx ); + fprintf(stdout," *gxy = %26.16e \n", *gxy ); + fprintf(stdout," *gxz = %26.16e \n", *gxz ); + fprintf(stdout," *gyy = %26.16e \n", *gyy ); + fprintf(stdout," *gyz = %26.16e \n", *gyz ); + fprintf(stdout," *gzz = %26.16e \n", *gzz ); + fprintf(stdout," *uxx = %26.16e \n", *uxx ); + fprintf(stdout," *uxy = %26.16e \n", *uxy ); + fprintf(stdout," *uxz = %26.16e \n", *uxz ); + fprintf(stdout," *uyy = %26.16e \n", *uyy ); + fprintf(stdout," *uyz = %26.16e \n", *uyz ); + fprintf(stdout," *uzz = %26.16e \n", *uzz ); + fprintf(stdout," *det = %26.16e \n", *det ); fprintf(stdout," *epsnegative = %10d \n", *epsnegative ); fprintf(stdout," *retval = %26.16e \n", *retval ); fflush(stdout); @@ -326,7 +326,7 @@ void CCTK_FCALL CCTK_FNAME(GRHydro_Con2PrimM_Polytype_pt) ( gammasq = 1. / (1. - vsq); gamma = sqrt(gammasq); - + // Always calculate rho from D and gamma so that using D in EOS remains consistent // i.e. you don't get positive values for dP/d(vsq) . rho0 = lg.D / gamma ; @@ -342,8 +342,8 @@ void CCTK_FCALL CCTK_FNAME(GRHydro_Con2PrimM_Polytype_pt) ( // Make sure that W is large enough so that v^2 < 1 : i_increase = 0; while( (( W_last*W_last*W_last * ( W_last + 2.*lg.Bsq ) - - lg.QdotBsq*(2.*W_last + lg.Bsq) ) <= W_last*W_last*(lg.Qtsq-lg.Bsq*lg.Bsq)) - && (i_increase < 10) ) { + - lg.QdotBsq*(2.*W_last + lg.Bsq) ) <= W_last*W_last*(lg.Qtsq-lg.Bsq*lg.Bsq)) + && (i_increase < 10) ) { W_last *= 10.; i_increase++; } @@ -357,7 +357,7 @@ void CCTK_FCALL CCTK_FNAME(GRHydro_Con2PrimM_Polytype_pt) ( *retval = 1.0*oned_newton_raphson( x_1d, 1, gammaeos, &lg, func_W ) ; W = x_1d[0]; - + /* Problem with solver, so return denoting error before doing anything further */ if( ((*retval) != 0.) || (W == FAIL_VAL) ) { *retval = *retval*100.+1.; @@ -487,9 +487,9 @@ static void validate_x(CCTK_REAL x[2], CCTK_REAL x0[2] ) *****************************************************************/ static CCTK_INT oned_newton_raphson( CCTK_REAL x[], CCTK_INT n, CCTK_REAL gammaeos, struct LocGlob *lgp, - void (*funcd) (CCTK_REAL [], CCTK_REAL [], CCTK_REAL [], - CCTK_REAL [][1], CCTK_REAL *, - CCTK_REAL *, CCTK_INT, CCTK_REAL, struct LocGlob *) ) + void (*funcd) (CCTK_REAL [], CCTK_REAL [], CCTK_REAL [], + CCTK_REAL [][1], CCTK_REAL *, + CCTK_REAL *, CCTK_INT, CCTK_REAL, struct LocGlob *) ) { CCTK_REAL f, df, dx[1], x_old[1], resid[1], jac[1][1]; @@ -528,8 +528,8 @@ static CCTK_INT oned_newton_raphson( CCTK_REAL x[], CCTK_INT n, CCTK_REAL gammae // //METHOD specific: // i_increase = 0; // while( (( x[0]*x[0]*x[0] * ( x[0] + 2.*lgp->Bsq ) - -// lgp->QdotBsq*(2.*x[0] + lgp->Bsq) ) <= x[0]*x[0]*(lgp->Qtsq-lgp->Bsq*lgp->Bsq)) -// && (i_increase < 10) ) { +// lgp->QdotBsq*(2.*x[0] + lgp->Bsq) ) <= x[0]*x[0]*(lgp->Qtsq-lgp->Bsq*lgp->Bsq)) +// && (i_increase < 10) ) { // x[0] -= (1.*i_increase) * dx[0] / 10. ; // i_increase++; // } @@ -561,7 +561,7 @@ static CCTK_INT oned_newton_raphson( CCTK_REAL x[], CCTK_INT n, CCTK_REAL gammae if( doing_extra == 1 ) i_extra++ ; if( ((fabs(errx) <= NEWT_TOL)&&(doing_extra == 0)) || - (i_extra > EXTRA_NEWT_ITER) || (n_iter >= (MAX_NEWT_ITER-1)) ) { + (i_extra > EXTRA_NEWT_ITER) || (n_iter >= (MAX_NEWT_ITER-1)) ) { keep_iterating = 0; } @@ -574,7 +574,7 @@ static CCTK_INT oned_newton_raphson( CCTK_REAL x[], CCTK_INT n, CCTK_REAL gammae if( (!finite(f)) || (!finite(df)) || (!finite(x[0])) ) { #if(DEBUG_CON2PRIMM) fprintf(stderr,"\ngnr not finite, f,df,x_o,x,W_o,W,rho_o,rho = %26.20e %26.20e %26.20e %26.20e %26.20e %26.20e %26.20e %26.20e \n", - f,df,x[0],x_old[0],lgp->W_for_gnr2_old,lgp->W_for_gnr2,lgp->rho_for_gnr2_old,lgp->rho_for_gnr2); fflush(stderr); + f,df,x[0],x_old[0],lgp->W_for_gnr2_old,lgp->W_for_gnr2,lgp->rho_for_gnr2_old,lgp->rho_for_gnr2); fflush(stderr); #endif return(2); } @@ -611,8 +611,8 @@ static CCTK_INT oned_newton_raphson( CCTK_REAL x[], CCTK_INT n, CCTK_REAL gammae n = dimension of x[]; *********************************************************************************/ static void func_W(CCTK_REAL x[], CCTK_REAL dx[], CCTK_REAL resid[], - CCTK_REAL jac[][1], CCTK_REAL *f, CCTK_REAL *df, CCTK_INT n, CCTK_REAL gammaeos, - struct LocGlob *lgp) + CCTK_REAL jac[][1], CCTK_REAL *f, CCTK_REAL *df, CCTK_INT n, CCTK_REAL gammaeos, + struct LocGlob *lgp) { CCTK_INT retval, ntries; CCTK_REAL W, x_rho[1], rho, rho_g ; @@ -689,8 +689,8 @@ static void func_W(CCTK_REAL x[], CCTK_REAL dx[], CCTK_REAL resid[], *****************************************************************/ static CCTK_INT gnr2( CCTK_REAL x[], CCTK_INT n, CCTK_REAL gammaeos, struct LocGlob *lgp, - void (*funcd) (CCTK_REAL [], CCTK_REAL [], CCTK_REAL [], - CCTK_REAL [][1],CCTK_REAL *,CCTK_REAL *,CCTK_INT,CCTK_REAL, struct LocGlob *lgp) ) + void (*funcd) (CCTK_REAL [], CCTK_REAL [], CCTK_REAL [], + CCTK_REAL [][1],CCTK_REAL *,CCTK_REAL *,CCTK_INT,CCTK_REAL, struct LocGlob *lgp) ) { CCTK_REAL f, df, dx[1], x_old[1], resid[1], jac[1][1]; @@ -746,7 +746,7 @@ static CCTK_INT gnr2( CCTK_REAL x[], CCTK_INT n, CCTK_REAL gammaeos, struct LocG // See if we've done the extra iterations, or have done too many iterations: if( ((fabs(errx) <= NEWT_TOL2)&&(doing_extra == 0)) || - (i_extra > EXTRA_NEWT_ITER) || (n_iter >= (MAX_NEWT_ITER-1)) ) { + (i_extra > EXTRA_NEWT_ITER) || (n_iter >= (MAX_NEWT_ITER-1)) ) { keep_iterating = 0; } @@ -759,7 +759,7 @@ static CCTK_INT gnr2( CCTK_REAL x[], CCTK_INT n, CCTK_REAL gammaeos, struct LocG if( (!finite(f)) || (!finite(df)) || (!finite(x[0])) ) { #if(DEBUG_CON2PRIMM) fprintf(stderr,"\ngnr2 not finite, f,df,x_o,x,W_o,W,rho_o,rho = %26.20e %26.20e %26.20e %26.20e %26.20e %26.20e %26.20e %26.20e \n", - f,df,x[0],x_old[0],lgp->W_for_gnr2_old,lgp->W_for_gnr2,lgp->rho_for_gnr2_old,lgp->rho_for_gnr2); fflush(stderr); + f,df,x[0],x_old[0],lgp->W_for_gnr2_old,lgp->W_for_gnr2,lgp->rho_for_gnr2_old,lgp->rho_for_gnr2); fflush(stderr); #endif return(2); } @@ -797,8 +797,8 @@ static CCTK_INT gnr2( CCTK_REAL x[], CCTK_INT n, CCTK_REAL gammaeos, struct LocG *********************************************************************************/ // for the isentropic version: eq. (27) static void func_rho(CCTK_REAL x[], CCTK_REAL dx[], CCTK_REAL resid[], - CCTK_REAL jac[][1], CCTK_REAL *f, CCTK_REAL *df, CCTK_INT n, CCTK_REAL gammaeos, - struct LocGlob *lgp) + CCTK_REAL jac[][1], CCTK_REAL *f, CCTK_REAL *df, CCTK_INT n, CCTK_REAL gammaeos, + struct LocGlob *lgp) { CCTK_REAL rho, W; diff --git a/src/GRHydro_Con2PrimM_pt.c b/src/GRHydro_Con2PrimM_pt.c index 865defa..b980887 100644 --- a/src/GRHydro_Con2PrimM_pt.c +++ b/src/GRHydro_Con2PrimM_pt.c @@ -61,9 +61,9 @@ #define NEWT_TOL2 (1.0e-15) /* TOL of new 1D^*_{v^2} gnr2 method */ #define MIN_NEWT_TOL2 (1.0e-10) /* TOL of new 1D^*_{v^2} gnr2 method */ -#define W_TOO_BIG (1.e20) /* \gamma^2 (\rho_0 + u + p) is assumed +#define W_TOO_BIG (1.e20) /* \gamma^2 (\rho_0 + u + p) is assumed to always be smaller than this. This - is used to detect solver failures */ + is used to detect solver failures */ #define FAIL_VAL (1.e30) /* Generic value to which we set variables when a problem arises */ @@ -80,12 +80,12 @@ struct LocGlob { static CCTK_REAL vsq_calc(CCTK_REAL W, struct LocGlob *lgp); static CCTK_INT twod_newton_raphson( CCTK_REAL x[], CCTK_REAL gammaeos, struct LocGlob *lgp, - void (*funcd) (CCTK_REAL [], CCTK_REAL [], - CCTK_REAL [], CCTK_REAL [][2], - CCTK_REAL *, CCTK_REAL *, CCTK_REAL, struct LocGlob *) ); + void (*funcd) (CCTK_REAL [], CCTK_REAL [], + CCTK_REAL [], CCTK_REAL [][2], + CCTK_REAL *, CCTK_REAL *, CCTK_REAL, struct LocGlob *) ); static void func_vsq( CCTK_REAL [], CCTK_REAL [], CCTK_REAL [], CCTK_REAL [][2], - CCTK_REAL *f, CCTK_REAL *df, CCTK_REAL gammaeos, struct LocGlob *lgp); + CCTK_REAL *f, CCTK_REAL *df, CCTK_REAL gammaeos, struct LocGlob *lgp); static CCTK_REAL x1_of_x0(CCTK_REAL x0, struct LocGlob *lgp ) ; @@ -177,11 +177,11 @@ RETURN VALUE: of retval = (i*100 + j) where given tolerances; 2 -> Newton-Raphson procedure encountered a numerical divergence (occurrence of "nan" or "+/-inf" ; - + j = 0 -> success 1 -> failure: some sort of failure in Newton-Raphson; 2 -> failure: unphysical vsq = v^2 value at initial guess; - 3 -> failure: W<0 or W>W_TOO_BIG + 3 -> failure: W<0 or W>W_TOO_BIG 4 -> failure: v^2 > 1 ( used to be 5 -> failure: rho,uu <= 0 but now sets epsnegative to non-zero ) @@ -235,37 +235,37 @@ void CCTK_FCALL CCTK_FNAME(GRHydro_Con2PrimM_pt) ( #if(DEBUG_CON2PRIMM) fprintf(stdout," *dens = %26.16e \n", *dens_in ); - fprintf(stdout," *sx = %26.16e \n", *sx_in ); - fprintf(stdout," *sy = %26.16e \n", *sy_in ); - fprintf(stdout," *sz = %26.16e \n", *sz_in ); - fprintf(stdout," *tau = %26.16e \n", *tau_in ); - fprintf(stdout," *Bconsx = %26.16e \n", *Bconsx_in ); - fprintf(stdout," *Bconsy = %26.16e \n", *Bconsy_in ); + fprintf(stdout," *sx = %26.16e \n", *sx_in ); + fprintf(stdout," *sy = %26.16e \n", *sy_in ); + fprintf(stdout," *sz = %26.16e \n", *sz_in ); + fprintf(stdout," *tau = %26.16e \n", *tau_in ); + fprintf(stdout," *Bconsx = %26.16e \n", *Bconsx_in ); + fprintf(stdout," *Bconsy = %26.16e \n", *Bconsy_in ); fprintf(stdout," *Bconsz = %26.16e \n", *Bconsz_in ); - fprintf(stdout," *rho = %26.16e \n", *rho ); - fprintf(stdout," *velx = %26.16e \n", *velx ); - fprintf(stdout," *vely = %26.16e \n", *vely ); - fprintf(stdout," *velz = %26.16e \n", *velz ); + fprintf(stdout," *rho = %26.16e \n", *rho ); + fprintf(stdout," *velx = %26.16e \n", *velx ); + fprintf(stdout," *vely = %26.16e \n", *vely ); + fprintf(stdout," *velz = %26.16e \n", *velz ); fprintf(stdout," *epsilon = %26.16e \n", *epsilon ); fprintf(stdout," *pressure = %26.16e \n", *pressure ); - fprintf(stdout," *Bx = %26.16e \n", *Bx ); - fprintf(stdout," *By = %26.16e \n", *By ); - fprintf(stdout," *Bz = %26.16e \n", *Bz ); - fprintf(stdout," *bsq = %26.16e \n", *bsq ); + fprintf(stdout," *Bx = %26.16e \n", *Bx ); + fprintf(stdout," *By = %26.16e \n", *By ); + fprintf(stdout," *Bz = %26.16e \n", *Bz ); + fprintf(stdout," *bsq = %26.16e \n", *bsq ); fprintf(stdout," *w_lorentz = %26.16e \n", *w_lorentz ); - fprintf(stdout," *gxx = %26.16e \n", *gxx ); - fprintf(stdout," *gxy = %26.16e \n", *gxy ); - fprintf(stdout," *gxz = %26.16e \n", *gxz ); - fprintf(stdout," *gyy = %26.16e \n", *gyy ); - fprintf(stdout," *gyz = %26.16e \n", *gyz ); - fprintf(stdout," *gzz = %26.16e \n", *gzz ); - fprintf(stdout," *uxx = %26.16e \n", *uxx ); - fprintf(stdout," *uxy = %26.16e \n", *uxy ); - fprintf(stdout," *uxz = %26.16e \n", *uxz ); - fprintf(stdout," *uyy = %26.16e \n", *uyy ); - fprintf(stdout," *uyz = %26.16e \n", *uyz ); - fprintf(stdout," *uzz = %26.16e \n", *uzz ); - fprintf(stdout," *det = %26.16e \n", *det ); + fprintf(stdout," *gxx = %26.16e \n", *gxx ); + fprintf(stdout," *gxy = %26.16e \n", *gxy ); + fprintf(stdout," *gxz = %26.16e \n", *gxz ); + fprintf(stdout," *gyy = %26.16e \n", *gyy ); + fprintf(stdout," *gyz = %26.16e \n", *gyz ); + fprintf(stdout," *gzz = %26.16e \n", *gzz ); + fprintf(stdout," *uxx = %26.16e \n", *uxx ); + fprintf(stdout," *uxy = %26.16e \n", *uxy ); + fprintf(stdout," *uxz = %26.16e \n", *uxz ); + fprintf(stdout," *uyy = %26.16e \n", *uyy ); + fprintf(stdout," *uyz = %26.16e \n", *uyz ); + fprintf(stdout," *uzz = %26.16e \n", *uzz ); + fprintf(stdout," *det = %26.16e \n", *det ); fprintf(stdout," *epsnegative = %10d \n", *epsnegative ); fprintf(stdout," *retval = %26.16e \n", *retval ); fflush(stdout); @@ -328,7 +328,7 @@ void CCTK_FCALL CCTK_FNAME(GRHydro_Con2PrimM_pt) ( gammasq = 1. / (1. - vsq); gamma = sqrt(gammasq); - + // Always calculate rho from D and gamma so that using D in EOS remains consistent // i.e. you don't get positive values for dP/d(vsq) . rho0 = lg.D / gamma ; @@ -342,8 +342,8 @@ void CCTK_FCALL CCTK_FNAME(GRHydro_Con2PrimM_pt) ( // Make sure that W is large enough so that v^2 < 1 : i_increase = 0; while( (( W_last*W_last*W_last * ( W_last + 2.*lg.Bsq ) - - lg.QdotBsq*(2.*W_last + lg.Bsq) ) <= W_last*W_last*(lg.Qtsq-lg.Bsq*lg.Bsq)) - && (i_increase < 10) ) { + - lg.QdotBsq*(2.*W_last + lg.Bsq) ) <= W_last*W_last*(lg.Qtsq-lg.Bsq*lg.Bsq)) + && (i_increase < 10) ) { W_last *= 10.; i_increase++; } @@ -355,7 +355,7 @@ void CCTK_FCALL CCTK_FNAME(GRHydro_Con2PrimM_pt) ( W = x_2d[0]; vsq = x_2d[1]; - + /* Problem with solver, so return denoting error before doing anything further */ if( ((*retval) != 0.) || (W == FAIL_VAL) ) { *retval = *retval*100.+1.; @@ -441,13 +441,13 @@ void CCTK_FCALL CCTK_FNAME(GRHydro_Con2PrimM_pt) ( ****************************************************************************/ static CCTK_REAL vsq_calc(CCTK_REAL W, struct LocGlob *lgp) { - CCTK_REAL Wsq,Xsq,Bsq_W; - - Wsq = W*W ; - Bsq_W = (lgp->Bsq + W); - Xsq = Bsq_W * Bsq_W; + CCTK_REAL Wsq,Xsq,Bsq_W; + + Wsq = W*W ; + Bsq_W = (lgp->Bsq + W); + Xsq = Bsq_W * Bsq_W; - return( ( Wsq * lgp->Qtsq + lgp->QdotBsq * (Bsq_W + W)) / (Wsq*Xsq) ); + return( ( Wsq * lgp->Qtsq + lgp->QdotBsq * (Bsq_W + W)) / (Wsq*Xsq) ); } @@ -509,9 +509,9 @@ static void validate_x(CCTK_REAL x[2], CCTK_REAL x0[2] ) *****************************************************************/ static CCTK_INT twod_newton_raphson( CCTK_REAL x[], CCTK_REAL gammaeos, struct LocGlob *lgp, - void (*funcd) (CCTK_REAL [], CCTK_REAL [], CCTK_REAL [], - CCTK_REAL [][2], CCTK_REAL *, - CCTK_REAL *, CCTK_REAL, struct LocGlob *) ) + void (*funcd) (CCTK_REAL [], CCTK_REAL [], CCTK_REAL [], + CCTK_REAL [][2], CCTK_REAL *, + CCTK_REAL *, CCTK_REAL, struct LocGlob *) ) { CCTK_REAL f, df, dx[2], x_old[2]; CCTK_REAL resid[2], jac[2][2]; @@ -580,7 +580,7 @@ static CCTK_INT twod_newton_raphson( CCTK_REAL x[], CCTK_REAL gammaeos, struct L if( doing_extra == 1 ) i_extra++ ; if( ((fabs(errx) <= NEWT_TOL)&&(doing_extra == 0)) - || (i_extra > EXTRA_NEWT_ITER) || (n_iter >= (MAX_NEWT_ITER-1)) ) { + || (i_extra > EXTRA_NEWT_ITER) || (n_iter >= (MAX_NEWT_ITER-1)) ) { keep_iterating = 0; } @@ -629,7 +629,7 @@ static CCTK_INT twod_newton_raphson( CCTK_REAL x[], CCTK_REAL gammaeos, struct L *********************************************************************************/ static void func_vsq(CCTK_REAL x[], CCTK_REAL dx[], CCTK_REAL resid[], - CCTK_REAL jac[][2], CCTK_REAL *f, CCTK_REAL *df, CCTK_REAL gammaeos, struct LocGlob *lgp) + CCTK_REAL jac[][2], CCTK_REAL *f, CCTK_REAL *df, CCTK_REAL gammaeos, struct LocGlob *lgp) { diff --git a/src/GRHydro_Con2PrimM_pt_EOSOmni.c b/src/GRHydro_Con2PrimM_pt_EOSOmni.c index 4fd6d1e..1d7e15c 100644 --- a/src/GRHydro_Con2PrimM_pt_EOSOmni.c +++ b/src/GRHydro_Con2PrimM_pt_EOSOmni.c @@ -61,9 +61,9 @@ #define NEWT_TOL2 (1.0e-15) /* TOL of new 1D^*_{v^2} gnr2 method */ #define MIN_NEWT_TOL2 (1.0e-10) /* TOL of new 1D^*_{v^2} gnr2 method */ -#define W_TOO_BIG (1.e20) /* \gamma^2 (\rho_0 + u + p) is assumed +#define W_TOO_BIG (1.e20) /* \gamma^2 (\rho_0 + u + p) is assumed to always be smaller than this. This - is used to detect solver failures */ + is used to detect solver failures */ #define FAIL_VAL (1.e30) /* Generic value to which we set variables when a problem arises */ @@ -86,21 +86,21 @@ struct eosomnivars { static CCTK_REAL vsq_calc(CCTK_REAL W, struct LocGlob *lgp); static CCTK_INT twod_newton_raphson( CCTK_REAL x[], CCTK_REAL gammaeos, struct LocGlob *lgp, - void (*funcd) (CCTK_REAL [], CCTK_REAL [], - CCTK_REAL [], CCTK_REAL [][2], - CCTK_REAL *, CCTK_REAL *, CCTK_REAL, struct LocGlob *) ); + void (*funcd) (CCTK_REAL [], CCTK_REAL [], + CCTK_REAL [], CCTK_REAL [][2], + CCTK_REAL *, CCTK_REAL *, CCTK_REAL, struct LocGlob *) ); static CCTK_INT threed_newton_raphson_omni( CCTK_REAL x[], struct eosomnivars *eosvars, struct LocGlob *lgp, - void (*funcd) (CCTK_REAL [], CCTK_REAL [], - CCTK_REAL [], CCTK_REAL [][3], - CCTK_REAL *, CCTK_REAL *, - struct eosomnivars *eosvars, struct LocGlob *) ); + void (*funcd) (CCTK_REAL [], CCTK_REAL [], + CCTK_REAL [], CCTK_REAL [][3], + CCTK_REAL *, CCTK_REAL *, + struct eosomnivars *eosvars, struct LocGlob *) ); static void func_vsq( CCTK_REAL [], CCTK_REAL [], CCTK_REAL [], CCTK_REAL [][2], - CCTK_REAL *f, CCTK_REAL *df, CCTK_REAL gammaeos, struct LocGlob *lgp); + CCTK_REAL *f, CCTK_REAL *df, CCTK_REAL gammaeos, struct LocGlob *lgp); static void func_vsq_eosomni( CCTK_REAL [], CCTK_REAL [], CCTK_REAL [], CCTK_REAL [][3], - CCTK_REAL *f, CCTK_REAL *df, struct eosomnivars *eosvars, struct LocGlob *lgp); + CCTK_REAL *f, CCTK_REAL *df, struct eosomnivars *eosvars, struct LocGlob *lgp); static CCTK_REAL x1_of_x0(CCTK_REAL x0, struct LocGlob *lgp ) ; @@ -196,11 +196,11 @@ RETURN VALUE: of retval = (i*100 + j) where given tolerances; 2 -> Newton-Raphson procedure encountered a numerical divergence (occurrence of "nan" or "+/-inf" ; - + j = 0 -> success 1 -> failure: some sort of failure in Newton-Raphson; 2 -> failure: unphysical vsq = v^2 value at initial guess; - 3 -> failure: W<0 or W>W_TOO_BIG + 3 -> failure: W<0 or W>W_TOO_BIG 4 -> failure: v^2 > 1 ( used to be 5 -> failure: rho,uu <= 0 but now sets epsnegative to non-zero ) @@ -264,39 +264,39 @@ void CCTK_FCALL CCTK_FNAME(GRHydro_Con2PrimM_pt) ( #if(DEBUG_CON2PRIMM) fprintf(stdout," *dens = %26.16e \n", *dens_in ); - fprintf(stdout," *sx = %26.16e \n", *sx_in ); - fprintf(stdout," *sy = %26.16e \n", *sy_in ); - fprintf(stdout," *sz = %26.16e \n", *sz_in ); - fprintf(stdout," *tau = %26.16e \n", *tau_in ); - fprintf(stdout," *Bconsx = %26.16e \n", *Bconsx_in ); - fprintf(stdout," *Bconsy = %26.16e \n", *Bconsy_in ); + fprintf(stdout," *sx = %26.16e \n", *sx_in ); + fprintf(stdout," *sy = %26.16e \n", *sy_in ); + fprintf(stdout," *sz = %26.16e \n", *sz_in ); + fprintf(stdout," *tau = %26.16e \n", *tau_in ); + fprintf(stdout," *Bconsx = %26.16e \n", *Bconsx_in ); + fprintf(stdout," *Bconsy = %26.16e \n", *Bconsy_in ); fprintf(stdout," *Bconsz = %26.16e \n", *Bconsz_in ); - fprintf(stdout," *rho = %26.16e \n", *rho ); - fprintf(stdout," *velx = %26.16e \n", *velx ); - fprintf(stdout," *vely = %26.16e \n", *vely ); - fprintf(stdout," *velz = %26.16e \n", *velz ); + fprintf(stdout," *rho = %26.16e \n", *rho ); + fprintf(stdout," *velx = %26.16e \n", *velx ); + fprintf(stdout," *vely = %26.16e \n", *vely ); + fprintf(stdout," *velz = %26.16e \n", *velz ); fprintf(stdout," *epsilon = %26.16e \n", *epsilon ); fprintf(stdout," *temp_in = %26.16e \n", *temp_in ); fprintf(stdout," *y_e_in = %26.16e \n", *y_e_in ); fprintf(stdout," *pressure = %26.16e \n", *pressure ); - fprintf(stdout," *Bx = %26.16e \n", *Bx ); - fprintf(stdout," *By = %26.16e \n", *By ); - fprintf(stdout," *Bz = %26.16e \n", *Bz ); - fprintf(stdout," *bsq = %26.16e \n", *bsq ); + fprintf(stdout," *Bx = %26.16e \n", *Bx ); + fprintf(stdout," *By = %26.16e \n", *By ); + fprintf(stdout," *Bz = %26.16e \n", *Bz ); + fprintf(stdout," *bsq = %26.16e \n", *bsq ); fprintf(stdout," *w_lorentz = %26.16e \n", *w_lorentz ); - fprintf(stdout," *gxx = %26.16e \n", *gxx ); - fprintf(stdout," *gxy = %26.16e \n", *gxy ); - fprintf(stdout," *gxz = %26.16e \n", *gxz ); - fprintf(stdout," *gyy = %26.16e \n", *gyy ); - fprintf(stdout," *gyz = %26.16e \n", *gyz ); - fprintf(stdout," *gzz = %26.16e \n", *gzz ); - fprintf(stdout," *uxx = %26.16e \n", *uxx ); - fprintf(stdout," *uxy = %26.16e \n", *uxy ); - fprintf(stdout," *uxz = %26.16e \n", *uxz ); - fprintf(stdout," *uyy = %26.16e \n", *uyy ); - fprintf(stdout," *uyz = %26.16e \n", *uyz ); - fprintf(stdout," *uzz = %26.16e \n", *uzz ); - fprintf(stdout," *det = %26.16e \n", *det ); + fprintf(stdout," *gxx = %26.16e \n", *gxx ); + fprintf(stdout," *gxy = %26.16e \n", *gxy ); + fprintf(stdout," *gxz = %26.16e \n", *gxz ); + fprintf(stdout," *gyy = %26.16e \n", *gyy ); + fprintf(stdout," *gyz = %26.16e \n", *gyz ); + fprintf(stdout," *gzz = %26.16e \n", *gzz ); + fprintf(stdout," *uxx = %26.16e \n", *uxx ); + fprintf(stdout," *uxy = %26.16e \n", *uxy ); + fprintf(stdout," *uxz = %26.16e \n", *uxz ); + fprintf(stdout," *uyy = %26.16e \n", *uyy ); + fprintf(stdout," *uyz = %26.16e \n", *uyz ); + fprintf(stdout," *uzz = %26.16e \n", *uzz ); + fprintf(stdout," *det = %26.16e \n", *det ); fprintf(stdout," *epsnegative = %10d \n", *epsnegative ); fprintf(stdout," *retval = %26.16e \n", *retval ); fflush(stdout); @@ -360,7 +360,7 @@ void CCTK_FCALL CCTK_FNAME(GRHydro_Con2PrimM_pt) ( gammasq = 1. / (1. - vsq); gamma = sqrt(gammasq); - + // Always calculate rho from D and gamma so that using D in EOS remains consistent // i.e. you don't get positive values for dP/d(vsq) . rho0 = lg.D / gamma ; @@ -379,13 +379,13 @@ void CCTK_FCALL CCTK_FNAME(GRHydro_Con2PrimM_pt) ( W_last = w*gammasq ; - //fprintf(stdout," p = %26.16e \n", p ); + //fprintf(stdout," p = %26.16e \n", p ); // Make sure that W is large enough so that v^2 < 1 : i_increase = 0; while( (( W_last*W_last*W_last * ( W_last + 2.*lg.Bsq ) - - lg.QdotBsq*(2.*W_last + lg.Bsq) ) <= W_last*W_last*(lg.Qtsq-lg.Bsq*lg.Bsq)) - && (i_increase < 10) ) { + - lg.QdotBsq*(2.*W_last + lg.Bsq) ) <= W_last*W_last*(lg.Qtsq-lg.Bsq*lg.Bsq)) + && (i_increase < 10) ) { W_last *= 10.; i_increase++; } @@ -406,7 +406,7 @@ void CCTK_FCALL CCTK_FNAME(GRHydro_Con2PrimM_pt) ( W = x_3d[0]; vsq = x_3d[1]; - + /* Problem with solver, so return denoting error before doing anything further */ if( ((*retval) != 0.) || (W == FAIL_VAL) ) { *retval = *retval*100.+1.; @@ -455,11 +455,11 @@ void CCTK_FCALL CCTK_FNAME(GRHydro_Con2PrimM_pt) ( if(rho0 <= 0.) { // *epsnegative = 1; fprintf(stdout," *epsnegative = %10d \n", *epsnegative ); - fprintf(stdout," rho0 = %26.16e \n", rho0 ); - fprintf(stdout," u = %26.16e \n", u ); - fprintf(stdout," W = %26.16e \n", W ); - fprintf(stdout," vsq = %26.16e \n", vsq ); - fprintf(stdout," uold = %26.16e \n", uold ); + fprintf(stdout," rho0 = %26.16e \n", rho0 ); + fprintf(stdout," u = %26.16e \n", u ); + fprintf(stdout," W = %26.16e \n", W ); + fprintf(stdout," vsq = %26.16e \n", vsq ); + fprintf(stdout," uold = %26.16e \n", uold ); return; } @@ -512,13 +512,13 @@ void CCTK_FCALL CCTK_FNAME(GRHydro_Con2PrimM_pt) ( ****************************************************************************/ static CCTK_REAL vsq_calc(CCTK_REAL W, struct LocGlob *lgp) { - CCTK_REAL Wsq,Xsq,Bsq_W; - - Wsq = W*W ; - Bsq_W = (lgp->Bsq + W); - Xsq = Bsq_W * Bsq_W; + CCTK_REAL Wsq,Xsq,Bsq_W; + + Wsq = W*W ; + Bsq_W = (lgp->Bsq + W); + Xsq = Bsq_W * Bsq_W; - return( ( Wsq * lgp->Qtsq + lgp->QdotBsq * (Bsq_W + W)) / (Wsq*Xsq) ); + return( ( Wsq * lgp->Qtsq + lgp->QdotBsq * (Bsq_W + W)) / (Wsq*Xsq) ); } @@ -578,9 +578,9 @@ static void validate_x(CCTK_REAL x[2], CCTK_REAL x0[2] ) *****************************************************************/ static CCTK_INT twod_newton_raphson( CCTK_REAL x[], CCTK_REAL gammaeos, struct LocGlob *lgp, - void (*funcd) (CCTK_REAL [], CCTK_REAL [], CCTK_REAL [], - CCTK_REAL [][2], CCTK_REAL *, - CCTK_REAL *, CCTK_REAL, struct LocGlob *) ) + void (*funcd) (CCTK_REAL [], CCTK_REAL [], CCTK_REAL [], + CCTK_REAL [][2], CCTK_REAL *, + CCTK_REAL *, CCTK_REAL, struct LocGlob *) ) { CCTK_REAL f, df, dx[2], x_old[2]; CCTK_REAL resid[2], jac[2][2]; @@ -649,7 +649,7 @@ static CCTK_INT twod_newton_raphson( CCTK_REAL x[], CCTK_REAL gammaeos, struct L if( doing_extra == 1 ) i_extra++ ; if( ((fabs(errx) <= NEWT_TOL)&&(doing_extra == 0)) - || (i_extra > EXTRA_NEWT_ITER) || (n_iter >= (MAX_NEWT_ITER-1)) ) { + || (i_extra > EXTRA_NEWT_ITER) || (n_iter >= (MAX_NEWT_ITER-1)) ) { keep_iterating = 0; } @@ -689,9 +689,9 @@ static CCTK_INT twod_newton_raphson( CCTK_REAL x[], CCTK_REAL gammaeos, struct L *****************************************************************/ static CCTK_INT threed_newton_raphson_omni( CCTK_REAL x[], struct eosomnivars *eosvars, struct LocGlob *lgp, - void (*funcd) (CCTK_REAL [], CCTK_REAL [], CCTK_REAL [], - CCTK_REAL [][3], CCTK_REAL *, - CCTK_REAL *, struct eosomnivars *, struct LocGlob *) ) + void (*funcd) (CCTK_REAL [], CCTK_REAL [], CCTK_REAL [], + CCTK_REAL [][3], CCTK_REAL *, + CCTK_REAL *, struct eosomnivars *, struct LocGlob *) ) { CCTK_REAL f, df, dx[3], x_old[3]; CCTK_REAL resid[3], jac[3][3]; @@ -766,7 +766,7 @@ static CCTK_INT threed_newton_raphson_omni( CCTK_REAL x[], struct eosomnivars *e if( doing_extra == 1 ) i_extra++ ; if( ((fabs(errx) <= NEWT_TOL)&&(doing_extra == 0)) - || (i_extra > EXTRA_NEWT_ITER) || (n_iter >= (MAX_NEWT_ITER-1)) ) { + || (i_extra > EXTRA_NEWT_ITER) || (n_iter >= (MAX_NEWT_ITER-1)) ) { keep_iterating = 0; } @@ -813,7 +813,7 @@ static CCTK_INT threed_newton_raphson_omni( CCTK_REAL x[], struct eosomnivars *e *********************************************************************************/ static void func_vsq(CCTK_REAL x[], CCTK_REAL dx[], CCTK_REAL resid[], - CCTK_REAL jac[][2], CCTK_REAL *f, CCTK_REAL *df, CCTK_REAL gammaeos, struct LocGlob *lgp) + CCTK_REAL jac[][2], CCTK_REAL *f, CCTK_REAL *df, CCTK_REAL gammaeos, struct LocGlob *lgp) { @@ -900,8 +900,8 @@ static void func_vsq(CCTK_REAL x[], CCTK_REAL dx[], CCTK_REAL resid[], *********************************************************************************/ static void func_vsq_eosomni(CCTK_REAL x[], CCTK_REAL dx[], CCTK_REAL resid[], - CCTK_REAL jac[][3], CCTK_REAL *f, CCTK_REAL *df, - struct eosomnivars *eosvars, struct LocGlob *lgp) + CCTK_REAL jac[][3], CCTK_REAL *f, CCTK_REAL *df, + struct eosomnivars *eosvars, struct LocGlob *lgp) { CCTK_REAL W, vsq, u, p_tmp,epsilon,Wsq, dPdvsq, dPdW; @@ -1034,16 +1034,16 @@ static CCTK_REAL pressure_rho0_eps_eosomni(CCTK_REAL rho,CCTK_REAL epsilon, CCTK epspt[0]=epsilon; EOS_Omni_press(eosvars->eoshandle,eosvars->eoskeytemp,eosvars->eosprec,1, - &rho,&epsilon,eosvars->eos_temp, - eosvars->eos_y_e,press,eosvars->eoskeyerr,eosvars->eosanyerr); + &rho,&epsilon,eosvars->eos_temp, + eosvars->eos_y_e,press,eosvars->eoskeyerr,eosvars->eosanyerr); EOS_Omni_DPressByDRho(eosvars->eoshandle,eosvars->eoskeytemp,eosvars->eosprec,1, - &rho,&epsilon,eosvars->eos_temp, - eosvars->eos_y_e,dpdrho,eosvars->eoskeyerr,eosvars->eosanyerr); + &rho,&epsilon,eosvars->eos_temp, + eosvars->eos_y_e,dpdrho,eosvars->eoskeyerr,eosvars->eosanyerr); EOS_Omni_DPressByDEps(eosvars->eoshandle,eosvars->eoskeytemp,eosvars->eosprec,1, - &rho,&epsilon,eosvars->eos_temp, - eosvars->eos_y_e,dpdeps,eosvars->eoskeyerr,eosvars->eosanyerr); + &rho,&epsilon,eosvars->eos_temp, + eosvars->eos_y_e,dpdeps,eosvars->eoskeyerr,eosvars->eosanyerr); return press[0]; diff --git a/src/GRHydro_Con2PrimM_pt_EOSOmni.c.hopefullyfixed.c b/src/GRHydro_Con2PrimM_pt_EOSOmni.c.hopefullyfixed.c index 8f2beb9..a08e933 100644 --- a/src/GRHydro_Con2PrimM_pt_EOSOmni.c.hopefullyfixed.c +++ b/src/GRHydro_Con2PrimM_pt_EOSOmni.c.hopefullyfixed.c @@ -61,9 +61,9 @@ #define NEWT_TOL2 (1.0e-15) /* TOL of new 1D^*_{v^2} gnr2 method */ #define MIN_NEWT_TOL2 (1.0e-10) /* TOL of new 1D^*_{v^2} gnr2 method */ -#define W_TOO_BIG (1.e20) /* \gamma^2 (\rho_0 + u + p) is assumed +#define W_TOO_BIG (1.e20) /* \gamma^2 (\rho_0 + u + p) is assumed to always be smaller than this. This - is used to detect solver failures */ + is used to detect solver failures */ #define FAIL_VAL (1.e30) /* Generic value to which we set variables when a problem arises */ @@ -86,21 +86,21 @@ struct eosomnivars { static CCTK_REAL vsq_calc(CCTK_REAL W, struct LocGlob *lgp); static CCTK_INT twod_newton_raphson( CCTK_REAL x[], CCTK_REAL gammaeos, struct LocGlob *lgp, - void (*funcd) (CCTK_REAL [], CCTK_REAL [], - CCTK_REAL [], CCTK_REAL [][2], - CCTK_REAL *, CCTK_REAL *, CCTK_REAL, struct LocGlob *) ); + void (*funcd) (CCTK_REAL [], CCTK_REAL [], + CCTK_REAL [], CCTK_REAL [][2], + CCTK_REAL *, CCTK_REAL *, CCTK_REAL, struct LocGlob *) ); static CCTK_INT threed_newton_raphson_omni( CCTK_REAL x[], struct eosomnivars *eosvars, struct LocGlob *lgp, - void (*funcd) (CCTK_REAL [], CCTK_REAL [], - CCTK_REAL [], CCTK_REAL [][3], - CCTK_REAL *, CCTK_REAL *, - struct eosomnivars *eosvars, struct LocGlob *) ); + void (*funcd) (CCTK_REAL [], CCTK_REAL [], + CCTK_REAL [], CCTK_REAL [][3], + CCTK_REAL *, CCTK_REAL *, + struct eosomnivars *eosvars, struct LocGlob *) ); static void func_vsq( CCTK_REAL [], CCTK_REAL [], CCTK_REAL [], CCTK_REAL [][2], - CCTK_REAL *f, CCTK_REAL *df, CCTK_REAL gammaeos, struct LocGlob *lgp); + CCTK_REAL *f, CCTK_REAL *df, CCTK_REAL gammaeos, struct LocGlob *lgp); static void func_vsq_eosomni( CCTK_REAL [], CCTK_REAL [], CCTK_REAL [], CCTK_REAL [][3], - CCTK_REAL *f, CCTK_REAL *df, struct eosomnivars *eosvars, struct LocGlob *lgp); + CCTK_REAL *f, CCTK_REAL *df, struct eosomnivars *eosvars, struct LocGlob *lgp); static CCTK_REAL x1_of_x0(CCTK_REAL x0, struct LocGlob *lgp ) ; @@ -196,11 +196,11 @@ RETURN VALUE: of retval = (i*100 + j) where given tolerances; 2 -> Newton-Raphson procedure encountered a numerical divergence (occurrence of "nan" or "+/-inf" ; - + j = 0 -> success 1 -> failure: some sort of failure in Newton-Raphson; 2 -> failure: unphysical vsq = v^2 value at initial guess; - 3 -> failure: W<0 or W>W_TOO_BIG + 3 -> failure: W<0 or W>W_TOO_BIG 4 -> failure: v^2 > 1 ( used to be 5 -> failure: rho,uu <= 0 but now sets epsnegative to non-zero ) @@ -264,39 +264,39 @@ void CCTK_FCALL CCTK_FNAME(GRHydro_Con2PrimM_pt) ( #if(DEBUG_CON2PRIMM) fprintf(stdout," *dens = %26.16e \n", *dens_in ); - fprintf(stdout," *sx = %26.16e \n", *sx_in ); - fprintf(stdout," *sy = %26.16e \n", *sy_in ); - fprintf(stdout," *sz = %26.16e \n", *sz_in ); - fprintf(stdout," *tau = %26.16e \n", *tau_in ); - fprintf(stdout," *Bconsx = %26.16e \n", *Bconsx_in ); - fprintf(stdout," *Bconsy = %26.16e \n", *Bconsy_in ); + fprintf(stdout," *sx = %26.16e \n", *sx_in ); + fprintf(stdout," *sy = %26.16e \n", *sy_in ); + fprintf(stdout," *sz = %26.16e \n", *sz_in ); + fprintf(stdout," *tau = %26.16e \n", *tau_in ); + fprintf(stdout," *Bconsx = %26.16e \n", *Bconsx_in ); + fprintf(stdout," *Bconsy = %26.16e \n", *Bconsy_in ); fprintf(stdout," *Bconsz = %26.16e \n", *Bconsz_in ); - fprintf(stdout," *rho = %26.16e \n", *rho ); - fprintf(stdout," *velx = %26.16e \n", *velx ); - fprintf(stdout," *vely = %26.16e \n", *vely ); - fprintf(stdout," *velz = %26.16e \n", *velz ); + fprintf(stdout," *rho = %26.16e \n", *rho ); + fprintf(stdout," *velx = %26.16e \n", *velx ); + fprintf(stdout," *vely = %26.16e \n", *vely ); + fprintf(stdout," *velz = %26.16e \n", *velz ); fprintf(stdout," *epsilon = %26.16e \n", *epsilon ); fprintf(stdout," *temp_in = %26.16e \n", *temp_in ); fprintf(stdout," *y_e_in = %26.16e \n", *y_e_in ); fprintf(stdout," *pressure = %26.16e \n", *pressure ); - fprintf(stdout," *Bx = %26.16e \n", *Bx ); - fprintf(stdout," *By = %26.16e \n", *By ); - fprintf(stdout," *Bz = %26.16e \n", *Bz ); - fprintf(stdout," *bsq = %26.16e \n", *bsq ); + fprintf(stdout," *Bx = %26.16e \n", *Bx ); + fprintf(stdout," *By = %26.16e \n", *By ); + fprintf(stdout," *Bz = %26.16e \n", *Bz ); + fprintf(stdout," *bsq = %26.16e \n", *bsq ); fprintf(stdout," *w_lorentz = %26.16e \n", *w_lorentz ); - fprintf(stdout," *gxx = %26.16e \n", *gxx ); - fprintf(stdout," *gxy = %26.16e \n", *gxy ); - fprintf(stdout," *gxz = %26.16e \n", *gxz ); - fprintf(stdout," *gyy = %26.16e \n", *gyy ); - fprintf(stdout," *gyz = %26.16e \n", *gyz ); - fprintf(stdout," *gzz = %26.16e \n", *gzz ); - fprintf(stdout," *uxx = %26.16e \n", *uxx ); - fprintf(stdout," *uxy = %26.16e \n", *uxy ); - fprintf(stdout," *uxz = %26.16e \n", *uxz ); - fprintf(stdout," *uyy = %26.16e \n", *uyy ); - fprintf(stdout," *uyz = %26.16e \n", *uyz ); - fprintf(stdout," *uzz = %26.16e \n", *uzz ); - fprintf(stdout," *det = %26.16e \n", *det ); + fprintf(stdout," *gxx = %26.16e \n", *gxx ); + fprintf(stdout," *gxy = %26.16e \n", *gxy ); + fprintf(stdout," *gxz = %26.16e \n", *gxz ); + fprintf(stdout," *gyy = %26.16e \n", *gyy ); + fprintf(stdout," *gyz = %26.16e \n", *gyz ); + fprintf(stdout," *gzz = %26.16e \n", *gzz ); + fprintf(stdout," *uxx = %26.16e \n", *uxx ); + fprintf(stdout," *uxy = %26.16e \n", *uxy ); + fprintf(stdout," *uxz = %26.16e \n", *uxz ); + fprintf(stdout," *uyy = %26.16e \n", *uyy ); + fprintf(stdout," *uyz = %26.16e \n", *uyz ); + fprintf(stdout," *uzz = %26.16e \n", *uzz ); + fprintf(stdout," *det = %26.16e \n", *det ); fprintf(stdout," *epsnegative = %10d \n", *epsnegative ); fprintf(stdout," *retval = %26.16e \n", *retval ); fflush(stdout); @@ -360,7 +360,7 @@ void CCTK_FCALL CCTK_FNAME(GRHydro_Con2PrimM_pt) ( gammasq = 1. / (1. - vsq); gamma = sqrt(gammasq); - + // Always calculate rho from D and gamma so that using D in EOS remains consistent // i.e. you don't get positive values for dP/d(vsq) . rho0 = lg.D / gamma ; @@ -379,13 +379,13 @@ void CCTK_FCALL CCTK_FNAME(GRHydro_Con2PrimM_pt) ( W_last = w*gammasq ; - //fprintf(stdout," p = %26.16e \n", p ); + //fprintf(stdout," p = %26.16e \n", p ); // Make sure that W is large enough so that v^2 < 1 : i_increase = 0; while( (( W_last*W_last*W_last * ( W_last + 2.*lg.Bsq ) - - lg.QdotBsq*(2.*W_last + lg.Bsq) ) <= W_last*W_last*(lg.Qtsq-lg.Bsq*lg.Bsq)) - && (i_increase < 10) ) { + - lg.QdotBsq*(2.*W_last + lg.Bsq) ) <= W_last*W_last*(lg.Qtsq-lg.Bsq*lg.Bsq)) + && (i_increase < 10) ) { W_last *= 10.; i_increase++; } @@ -406,7 +406,7 @@ void CCTK_FCALL CCTK_FNAME(GRHydro_Con2PrimM_pt) ( W = x_3d[0]; vsq = x_3d[1]; - + /* Problem with solver, so return denoting error before doing anything further */ if( ((*retval) != 0.) || (W == FAIL_VAL) ) { *retval = *retval*100.+1.; @@ -452,20 +452,20 @@ void CCTK_FCALL CCTK_FNAME(GRHydro_Con2PrimM_pt) ( if((*handle==1 || *handle==2) && (rho0 <= 0.) || (u <= 0.)) { *epsnegative = 1; fprintf(stdout," *epsnegative = %10d \n", *epsnegative ); - fprintf(stdout," rho0 = %26.16e \n", rho0 ); - fprintf(stdout," u = %26.16e \n", u ); - fprintf(stdout," W = %26.16e \n", W ); - fprintf(stdout," vsq = %26.16e \n", vsq ); - fprintf(stdout," uold = %26.16e \n", uold ); + fprintf(stdout," rho0 = %26.16e \n", rho0 ); + fprintf(stdout," u = %26.16e \n", u ); + fprintf(stdout," W = %26.16e \n", W ); + fprintf(stdout," vsq = %26.16e \n", vsq ); + fprintf(stdout," uold = %26.16e \n", uold ); return; } else if(rho0 <= 0.) { // *epsnegative = 1; fprintf(stdout," *epsnegative = %10d \n", *epsnegative ); - fprintf(stdout," rho0 = %26.16e \n", rho0 ); - fprintf(stdout," u = %26.16e \n", u ); - fprintf(stdout," W = %26.16e \n", W ); - fprintf(stdout," vsq = %26.16e \n", vsq ); - fprintf(stdout," uold = %26.16e \n", uold ); + fprintf(stdout," rho0 = %26.16e \n", rho0 ); + fprintf(stdout," u = %26.16e \n", u ); + fprintf(stdout," W = %26.16e \n", W ); + fprintf(stdout," vsq = %26.16e \n", vsq ); + fprintf(stdout," uold = %26.16e \n", uold ); return; } @@ -523,13 +523,13 @@ void CCTK_FCALL CCTK_FNAME(GRHydro_Con2PrimM_pt) ( ****************************************************************************/ static CCTK_REAL vsq_calc(CCTK_REAL W, struct LocGlob *lgp) { - CCTK_REAL Wsq,Xsq,Bsq_W; - - Wsq = W*W ; - Bsq_W = (lgp->Bsq + W); - Xsq = Bsq_W * Bsq_W; + CCTK_REAL Wsq,Xsq,Bsq_W; + + Wsq = W*W ; + Bsq_W = (lgp->Bsq + W); + Xsq = Bsq_W * Bsq_W; - return( ( Wsq * lgp->Qtsq + lgp->QdotBsq * (Bsq_W + W)) / (Wsq*Xsq) ); + return( ( Wsq * lgp->Qtsq + lgp->QdotBsq * (Bsq_W + W)) / (Wsq*Xsq) ); } @@ -589,9 +589,9 @@ static void validate_x(CCTK_REAL x[2], CCTK_REAL x0[2] ) *****************************************************************/ static CCTK_INT twod_newton_raphson( CCTK_REAL x[], CCTK_REAL gammaeos, struct LocGlob *lgp, - void (*funcd) (CCTK_REAL [], CCTK_REAL [], CCTK_REAL [], - CCTK_REAL [][2], CCTK_REAL *, - CCTK_REAL *, CCTK_REAL, struct LocGlob *) ) + void (*funcd) (CCTK_REAL [], CCTK_REAL [], CCTK_REAL [], + CCTK_REAL [][2], CCTK_REAL *, + CCTK_REAL *, CCTK_REAL, struct LocGlob *) ) { CCTK_REAL f, df, dx[2], x_old[2]; CCTK_REAL resid[2], jac[2][2]; @@ -660,7 +660,7 @@ static CCTK_INT twod_newton_raphson( CCTK_REAL x[], CCTK_REAL gammaeos, struct L if( doing_extra == 1 ) i_extra++ ; if( ((fabs(errx) <= NEWT_TOL)&&(doing_extra == 0)) - || (i_extra > EXTRA_NEWT_ITER) || (n_iter >= (MAX_NEWT_ITER-1)) ) { + || (i_extra > EXTRA_NEWT_ITER) || (n_iter >= (MAX_NEWT_ITER-1)) ) { keep_iterating = 0; } @@ -700,9 +700,9 @@ static CCTK_INT twod_newton_raphson( CCTK_REAL x[], CCTK_REAL gammaeos, struct L *****************************************************************/ static CCTK_INT threed_newton_raphson_omni( CCTK_REAL x[], struct eosomnivars *eosvars, struct LocGlob *lgp, - void (*funcd) (CCTK_REAL [], CCTK_REAL [], CCTK_REAL [], - CCTK_REAL [][3], CCTK_REAL *, - CCTK_REAL *, struct eosomnivars *, struct LocGlob *) ) + void (*funcd) (CCTK_REAL [], CCTK_REAL [], CCTK_REAL [], + CCTK_REAL [][3], CCTK_REAL *, + CCTK_REAL *, struct eosomnivars *, struct LocGlob *) ) { CCTK_REAL f, df, dx[3], x_old[3]; CCTK_REAL resid[3], jac[3][3]; @@ -777,7 +777,7 @@ static CCTK_INT threed_newton_raphson_omni( CCTK_REAL x[], struct eosomnivars *e if( doing_extra == 1 ) i_extra++ ; if( ((fabs(errx) <= NEWT_TOL)&&(doing_extra == 0)) - || (i_extra > EXTRA_NEWT_ITER) || (n_iter >= (MAX_NEWT_ITER-1)) ) { + || (i_extra > EXTRA_NEWT_ITER) || (n_iter >= (MAX_NEWT_ITER-1)) ) { keep_iterating = 0; } @@ -824,7 +824,7 @@ static CCTK_INT threed_newton_raphson_omni( CCTK_REAL x[], struct eosomnivars *e *********************************************************************************/ static void func_vsq(CCTK_REAL x[], CCTK_REAL dx[], CCTK_REAL resid[], - CCTK_REAL jac[][2], CCTK_REAL *f, CCTK_REAL *df, CCTK_REAL gammaeos, struct LocGlob *lgp) + CCTK_REAL jac[][2], CCTK_REAL *f, CCTK_REAL *df, CCTK_REAL gammaeos, struct LocGlob *lgp) { @@ -911,8 +911,8 @@ static void func_vsq(CCTK_REAL x[], CCTK_REAL dx[], CCTK_REAL resid[], *********************************************************************************/ static void func_vsq_eosomni(CCTK_REAL x[], CCTK_REAL dx[], CCTK_REAL resid[], - CCTK_REAL jac[][3], CCTK_REAL *f, CCTK_REAL *df, - struct eosomnivars *eosvars, struct LocGlob *lgp) + CCTK_REAL jac[][3], CCTK_REAL *f, CCTK_REAL *df, + struct eosomnivars *eosvars, struct LocGlob *lgp) { CCTK_REAL W, vsq, u, p_tmp,epsilon,Wsq, dPdvsq, dPdW; @@ -1045,16 +1045,16 @@ static CCTK_REAL pressure_rho0_eps_eosomni(CCTK_REAL rho,CCTK_REAL epsilon, CCTK epspt[0]=epsilon; EOS_Omni_press(eosvars->eoshandle,eosvars->eoskeytemp,eosvars->eosprec,1, - &rho,&epsilon,eosvars->eos_temp, - eosvars->eos_y_e,press,eosvars->eoskeyerr,eosvars->eosanyerr); + &rho,&epsilon,eosvars->eos_temp, + eosvars->eos_y_e,press,eosvars->eoskeyerr,eosvars->eosanyerr); EOS_Omni_DPressByDRho(eosvars->eoshandle,eosvars->eoskeytemp,eosvars->eosprec,1, - &rho,&epsilon,eosvars->eos_temp, - eosvars->eos_y_e,dpdrho,eosvars->eoskeyerr,eosvars->eosanyerr); + &rho,&epsilon,eosvars->eos_temp, + eosvars->eos_y_e,dpdrho,eosvars->eoskeyerr,eosvars->eosanyerr); EOS_Omni_DPressByDEps(eosvars->eoshandle,eosvars->eoskeytemp,eosvars->eosprec,1, - &rho,&epsilon,eosvars->eos_temp, - eosvars->eos_y_e,dpdeps,eosvars->eoskeyerr,eosvars->eosanyerr); + &rho,&epsilon,eosvars->eos_temp, + eosvars->eos_y_e,dpdeps,eosvars->eoskeyerr,eosvars->eosanyerr); return press[0]; diff --git a/src/GRHydro_EoSChangeGamma.F90 b/src/GRHydro_EoSChangeGamma.F90 index 52c2d3d..ce8440f 100644 --- a/src/GRHydro_EoSChangeGamma.F90 +++ b/src/GRHydro_EoSChangeGamma.F90 @@ -339,7 +339,7 @@ subroutine GRHydro_EoSChangeGammaK_Shibata(CCTK_ARGUMENTS) eos_k_initial = initial_k press = (local_Gamma - 1.d0) / (initial_Gamma - 1.0d0 ) * eos_k_initial * & - rho ** initial_Gamma + rho ** initial_Gamma eps = eos_k_initial * & rho ** initial_Gamma / & diff --git a/src/GRHydro_InterfacesAM.h b/src/GRHydro_InterfacesAM.h index c2565d8..154b8fc 100644 --- a/src/GRHydro_InterfacesAM.h +++ b/src/GRHydro_InterfacesAM.h @@ -4,7 +4,7 @@ module Con2PrimAM_fortran_interfaces interface subroutine GRHydro_Con2PrimM_pt( handle, keytemp, prec, & - local_gam, dens, & + local_gam, dens, & sx, sy, sz, & tau, & Bconsx, Bconsy, Bconsz, & @@ -20,8 +20,8 @@ module Con2PrimAM_fortran_interfaces uxx, uxy, uxz,& uyy, uyz, uzz,& det,& - epsnegative, & - retval) + epsnegative, & + retval) implicit none CCTK_INT handle @@ -65,8 +65,8 @@ module Con2PrimAM_fortran_interfaces uxx, uxy, uxz,& uyy, uyz, uzz,& det,& - epsnegative, & - retval) + epsnegative, & + retval) implicit none CCTK_INT handle diff --git a/src/GRHydro_InterfacesM.h b/src/GRHydro_InterfacesM.h index 5463cc7..758667f 100644 --- a/src/GRHydro_InterfacesM.h +++ b/src/GRHydro_InterfacesM.h @@ -4,7 +4,7 @@ module Con2PrimM_fortran_interfaces interface subroutine GRHydro_Con2PrimM_pt( handle, keytemp, prec, & - local_gam, dens, & + local_gam, dens, & sx, sy, sz, & tau, & Bconsx, Bconsy, Bconsz, & @@ -20,8 +20,8 @@ module Con2PrimM_fortran_interfaces uxx, uxy, uxz,& uyy, uyz, uzz,& det,& - epsnegative, & - retval) + epsnegative, & + retval) implicit none CCTK_INT handle @@ -65,8 +65,8 @@ module Con2PrimM_fortran_interfaces uxx, uxy, uxz,& uyy, uyz, uzz,& det,& - epsnegative, & - retval) + epsnegative, & + retval) implicit none CCTK_INT handle diff --git a/src/GRHydro_Macros.h b/src/GRHydro_Macros.h index 925cfb4..2f1f532 100644 --- a/src/GRHydro_Macros.h +++ b/src/GRHydro_Macros.h @@ -7,7 +7,7 @@ (gxy_)*( (x1_)*(y2_)+(y1_)*(x2_) )+(gxz_)*( (x1_)*(z2_)+(z1_)*(x2_) )+\ (gyz_)*( (y1_)*(z2_)+(z1_)*(y2_) ) ) -#define DOTP2(gxx_,gxy_,gxz_,gyy_,gyz_,gzz_,x_,y_,z_) \ +#define DOTP2(gxx_,gxy_,gxz_,gyy_,gyz_,gzz_,x_,y_,z_) \ ( (gxx_)*(x_)**2+(gyy_)*(y_)**2+(gzz_)*(z_)**2+ \ 2.0*( (gxy_)*(x_)*(y_)+(gxz_)*(x_)*(z_)+(gyz_)*(y_)*(z_) ) ) diff --git a/src/GRHydro_PPM.F90 b/src/GRHydro_PPM.F90 index a6c46dc..f71a0e1 100644 --- a/src/GRHydro_PPM.F90 +++ b/src/GRHydro_PPM.F90 @@ -130,13 +130,13 @@ trivial_rp = .true. #define STEEP(x,dx,dmx) \ - if ( (x(i+1) - x(i)) * (x(i) - x(i-1)) > 0.d0 ) then &&\ - dmx(i) = sign(one, dx(i)) * \ - min(abs(dx(i)), 2.d0 * abs(x(i) - x(i-1)), \ - 2.d0 * abs(x(i+1) - x(i))) &&\ - else &&\ - dmx(i) = 0.d0 &&\ - end if + if ( (x(i+1) - x(i)) * (x(i) - x(i-1)) > 0.d0 ) then &&\ + dmx(i) = sign(one, dx(i)) * \ + min(abs(dx(i)), 2.d0 * abs(x(i) - x(i-1)), \ + 2.d0 * abs(x(i+1) - x(i))) &&\ + else &&\ + dmx(i) = 0.d0 &&\ + end if if (use_enhanced_ppm .eq. 0) then @@ -146,57 +146,57 @@ trivial_rp = .true. !!$ This is the expression for an even grid. do i = 2, nx - 1 - drho(i) = 0.5d0 * (rho(i+1) - rho(i-1)) - dvelx(i) = 0.5d0 * (vx(i+1) - vx(i-1)) - dvely(i) = 0.5d0 * (vy(i+1) - vy(i-1)) - dvelz(i) = 0.5d0 * (vz(i+1) - vz(i-1)) - dpress(i) = press(i+1) - press(i-1) - d2rho(i) = (rho(i+1) - 2.d0 * rho(i) + rho(i-1))! / 6.d0 / dx / dx - ! since we use d2rho only for the condition d2rho(i+1)*d2rhoi(i-1)<0 - ! the denominator is not necessary + drho(i) = 0.5d0 * (rho(i+1) - rho(i-1)) + dvelx(i) = 0.5d0 * (vx(i+1) - vx(i-1)) + dvely(i) = 0.5d0 * (vy(i+1) - vy(i-1)) + dvelz(i) = 0.5d0 * (vz(i+1) - vz(i-1)) + dpress(i) = press(i+1) - press(i-1) + d2rho(i) = (rho(i+1) - 2.d0 * rho(i) + rho(i-1))! / 6.d0 / dx / dx + ! since we use d2rho only for the condition d2rho(i+1)*d2rhoi(i-1)<0 + ! the denominator is not necessary end do if (poly .eq. 0) then - do i = 2, nx - 1 - deps(i) = 0.5d0 * (eps(i+1) - eps(i-1)) - end do + do i = 2, nx - 1 + deps(i) = 0.5d0 * (eps(i+1) - eps(i-1)) + end do end if !!$ Steepened slope. See (1.8) of Colella and Woodward, p.178 do i = 2, nx - 1 - STEEP(rho, drho, dmrho) - STEEP(vx, dvelx, dmvelx) - STEEP(vy, dvely, dmvely) - STEEP(vz, dvelz, dmvelz) + STEEP(rho, drho, dmrho) + STEEP(vx, dvelx, dmvelx) + STEEP(vy, dvely, dmvely) + STEEP(vz, dvelz, dmvelz) end do if (poly .eq. 0) then - do i = 2, nx - 1 - STEEP(eps, deps, dmeps) - end do + do i = 2, nx - 1 + STEEP(eps, deps, dmeps) + end do end if !!$ Initial boundary states. See (1.9) of Colella and Woodward, p.178 do i = 2, nx-2 - rhoplus(i) = 0.5d0 * (rho(i) + rho(i+1)) + & - (dmrho(i) - dmrho(i+1)) / 6.d0 - rhominus(i+1) = rhoplus(i) - velxplus(i) = 0.5d0 * (vx(i) + vx(i+1)) + & - (dmvelx(i) - dmvelx(i+1)) / 6.d0 - velxminus(i+1) = velxplus(i) - velyplus(i) = 0.5d0 * (vy(i) + vy(i+1)) + & - (dmvely(i) - dmvely(i+1)) / 6.d0 - velyminus(i+1) = velyplus(i) - velzplus(i) = 0.5d0 * (vz(i) + vz(i+1)) + & - (dmvelz(i) - dmvelz(i+1)) / 6.d0 - velzminus(i+1) = velzplus(i) + rhoplus(i) = 0.5d0 * (rho(i) + rho(i+1)) + & + (dmrho(i) - dmrho(i+1)) / 6.d0 + rhominus(i+1) = rhoplus(i) + velxplus(i) = 0.5d0 * (vx(i) + vx(i+1)) + & + (dmvelx(i) - dmvelx(i+1)) / 6.d0 + velxminus(i+1) = velxplus(i) + velyplus(i) = 0.5d0 * (vy(i) + vy(i+1)) + & + (dmvely(i) - dmvely(i+1)) / 6.d0 + velyminus(i+1) = velyplus(i) + velzplus(i) = 0.5d0 * (vz(i) + vz(i+1)) + & + (dmvelz(i) - dmvelz(i+1)) / 6.d0 + velzminus(i+1) = velzplus(i) end do if (poly .eq. 0) then - do i = 2, nx-2 - epsplus(i) = 0.5d0 * (eps(i) + eps(i+1)) + & - (dmeps(i) - dmeps(i+1)) / 6.d0 - epsminus(i+1) = epsplus(i) - end do + do i = 2, nx-2 + epsplus(i) = 0.5d0 * (eps(i) + eps(i+1)) + & + (dmeps(i) - dmeps(i+1)) / 6.d0 + epsminus(i+1) = epsplus(i) + end do end if else !! This is the modified PPM algorithm by Colella & Sekora 2008 and McCorquodale & Colella 2011. @@ -213,32 +213,32 @@ trivial_rp = .true. #define LIMIT(a,ah,C, alim) \ if ((min(a(i),a(i+1)) .le. ah) .and. (ah .le. max(a(i),a(i+1)))) then &&\ - alim = ah &&\ + alim = ah &&\ else &&\ - D2a = 3.0d0 * (a(i) - 2.0d0*ah + a(i+1)) &&\ - D2aL = (a(i-1) - 2.0d0*a(i) + a(i+1)) &&\ - D2aR = (a(i) - 2.0d0*a(i+1) + a(i+2)) &&\ - D2aLim = sign(one, D2a)*min(C*abs(D2aL), C*abs(D2aR), abs(D2a))/3.0d0 &&\ - if (D2a*D2aR .ge. 0 .and. D2a*D2aL .ge. 0) then &&\ - alim = 0.5d0*(a(i)+a(i+1)) - D2aLim &&\ - else &&\ - alim = 0.5d0*(a(i)+a(i+1)) &&\ - end if &&\ + D2a = 3.0d0 * (a(i) - 2.0d0*ah + a(i+1)) &&\ + D2aL = (a(i-1) - 2.0d0*a(i) + a(i+1)) &&\ + D2aR = (a(i) - 2.0d0*a(i+1) + a(i+2)) &&\ + D2aLim = sign(one, D2a)*min(C*abs(D2aL), C*abs(D2aR), abs(D2a))/3.0d0 &&\ + if (D2a*D2aR .ge. 0 .and. D2a*D2aL .ge. 0) then &&\ + alim = 0.5d0*(a(i)+a(i+1)) - D2aLim &&\ + else &&\ + alim = 0.5d0*(a(i)+a(i+1)) &&\ + end if &&\ endif #define LIMIT_EPS(a,ah,C, alim) \ if ((min(a(i),a(i+1)) .le. ah) .and. (ah .le. max(a(i),a(i+1)))) then &&\ - alim = ah &&\ + alim = ah &&\ else &&\ - D2a = 3.0d0 * (a(i) - 2.0d0*ah + a(i+1)) &&\ - D2aL = (a(i-1) - 2.0d0*a(i) + a(i+1)) &&\ - D2aR = (a(i) - 2.0d0*a(i+1) + a(i+2)) &&\ - D2aLim = sign(one, D2a)*min(C*abs(D2aL), C*abs(D2aR), abs(D2a))/3.0d0 &&\ - if (D2a*D2aR .ge. 0 .and. D2a*D2aL .ge. 0 .and. abs(D2aLim) .lt. abs(0.5d0*(a(i)+a(i+1)))) then &&\ - alim = 0.5d0*(a(i)+a(i+1)) - D2aLim &&\ - else &&\ - alim = 0.5d0*(a(i)+a(i+1)) &&\ - end if &&\ + D2a = 3.0d0 * (a(i) - 2.0d0*ah + a(i+1)) &&\ + D2aL = (a(i-1) - 2.0d0*a(i) + a(i+1)) &&\ + D2aR = (a(i) - 2.0d0*a(i+1) + a(i+2)) &&\ + D2aLim = sign(one, D2a)*min(C*abs(D2aL), C*abs(D2aR), abs(D2a))/3.0d0 &&\ + if (D2a*D2aR .ge. 0 .and. D2a*D2aL .ge. 0 .and. abs(D2aLim) .lt. abs(0.5d0*(a(i)+a(i+1)))) then &&\ + alim = 0.5d0*(a(i)+a(i+1)) - D2aLim &&\ + else &&\ + alim = 0.5d0*(a(i)+a(i+1)) &&\ + end if &&\ endif @@ -266,11 +266,11 @@ trivial_rp = .true. velzminus(i+1) = velzplus(i) end do if (poly .eq. 0) then - do i = 2, nx-2 - APPROX_AT_CELL_INTERFACE(eps, epsplus(i)) - LIMIT_EPS(eps, epsplus(i), enhanced_ppm_C2, epsplus(i)) - epsminus(i+1) = epsplus(i) - end do + do i = 2, nx-2 + APPROX_AT_CELL_INTERFACE(eps, epsplus(i)) + LIMIT_EPS(eps, epsplus(i), enhanced_ppm_C2, epsplus(i)) + epsminus(i+1) = epsplus(i) + end do end if else @@ -295,11 +295,11 @@ trivial_rp = .true. velzminus(i+1) = velzplus(i) end do if (poly .eq. 0) then - do i = 3, nx-3 - APPROX_AT_CELL_INTERFACE_STENCIL4(eps, epsplus(i)) - LIMIT(eps, epsplus(i), enhanced_ppm_C2, epsplus(i)) - epsminus(i+1) = epsplus(i) - end do + do i = 3, nx-3 + APPROX_AT_CELL_INTERFACE_STENCIL4(eps, epsplus(i)) + LIMIT(eps, epsplus(i), enhanced_ppm_C2, epsplus(i)) + epsminus(i+1) = epsplus(i) + end do end if endif @@ -323,36 +323,36 @@ trivial_rp = .true. if (ppm_detect .ne. 0) then if (use_enhanced_ppm .eq. 1) then - ! make sure drho, d2rho and dmrho are computed everywhere! - do i = 2, nx-1 - drho(i) = 0.5d0 * (rho(i+1) - rho(i-1)) - d2rho(i) = (rho(i+1) - 2.d0 * rho(i) + rho(i-1)) - end do - do i = 2, nx-1 - STEEP(rho, drho, dmrho) - end do + ! make sure drho, d2rho and dmrho are computed everywhere! + do i = 2, nx-1 + drho(i) = 0.5d0 * (rho(i+1) - rho(i-1)) + d2rho(i) = (rho(i+1) - 2.d0 * rho(i) + rho(i-1)) + end do + do i = 2, nx-1 + STEEP(rho, drho, dmrho) + end do endif do i = 3, nx - 2 - if ( (d2rho(i+1)*d2rho(i-1) < 0.d0).and.(abs(rho(i+1)-rho(i-1)) - & - ppm_epsilon_shock * min(abs(rho(i+1)), abs(rho(i-1))) > 0.d0) ) then - etatilde = (rho(i-2) - rho(i+2) + 4.d0 * drho(i)) / (drho(i) * 12.d0) - else - etatilde = 0.d0 - end if - eta = max(0.d0, min(1.d0, ppm_eta1 * (etatilde - ppm_eta2))) - if (ppm_k0 * abs(drho(i)) * min(press(i-1),press(i+1)) < & - abs(dpress(i)) * min(rho(i-1), rho(i+1))) then - eta = 0.d0 - end if - if (eta > 0.d0) then - trivial_rp(i-1) = .false. - trivial_rp(i) = .false. - end if - rhominus(i) = rhominus(i) * (1.d0 - eta) + & - (rho(i-1) + 0.5d0 * dmrho(i-1)) * eta - rhoplus(i) = rhoplus(i) * (1.d0 - eta) + & - (rho(i+1) - 0.5d0 * dmrho(i+1)) * eta + if ( (d2rho(i+1)*d2rho(i-1) < 0.d0).and.(abs(rho(i+1)-rho(i-1)) - & + ppm_epsilon_shock * min(abs(rho(i+1)), abs(rho(i-1))) > 0.d0) ) then + etatilde = (rho(i-2) - rho(i+2) + 4.d0 * drho(i)) / (drho(i) * 12.d0) + else + etatilde = 0.d0 + end if + eta = max(0.d0, min(1.d0, ppm_eta1 * (etatilde - ppm_eta2))) + if (ppm_k0 * abs(drho(i)) * min(press(i-1),press(i+1)) < & + abs(dpress(i)) * min(rho(i-1), rho(i+1))) then + eta = 0.d0 + end if + if (eta > 0.d0) then + trivial_rp(i-1) = .false. + trivial_rp(i) = .false. + end if + rhominus(i) = rhominus(i) * (1.d0 - eta) + & + (rho(i-1) + 0.5d0 * dmrho(i-1)) * eta + rhoplus(i) = rhoplus(i) * (1.d0 - eta) + & + (rho(i+1) - 0.5d0 * dmrho(i+1)) * eta end do end if @@ -366,95 +366,95 @@ trivial_rp = .true. l_ev_r=0.d0 xwind=0.d0 do i=3, nx - 3 - agxx = 0.5d0*( gxx(i) + gxx(i+1) ) - agxy = 0.5d0*( gxy(i) + gxy(i+1) ) - agxz = 0.5d0*( gxz(i) + gxz(i+1) ) - agyy = 0.5d0*( gyy(i) + gyy(i+1) ) - agyz = 0.5d0*( gyz(i) + gyz(i+1) ) - agzz = 0.5d0*( gzz(i) + gzz(i+1) ) - det = SPATIAL_DETERMINANT(agxx, agxy, agxz, \ - agyy, agyz, agzz) - call UpperMetric (uxx, uxy, uxz, uyy, uyz, uzz, & - det, agxx, agxy, agxz, agyy, agyz, agzz) - call eigenvalues(handle,& - D_UPW(rho), D_UPW(velx), D_UPW(vely), D_UPW(velz), & - D_UPW(eps), D_UPW(w_lorentz), lam, & - agxx, agxy, agxz, agyy, agyz, agzz, & - uxx, D_UPW(alp), D_UPW(beta)) - l_ev_l(i)=lam(1) - l_ev_r(i)=lam(5) - xwind(i) = (lam(1) + lam(5)) / (abs(lam(1)) + abs(lam(5))) - xwind(i) = min(1.d0, max(-1.d0, xwind(i))) + agxx = 0.5d0*( gxx(i) + gxx(i+1) ) + agxy = 0.5d0*( gxy(i) + gxy(i+1) ) + agxz = 0.5d0*( gxz(i) + gxz(i+1) ) + agyy = 0.5d0*( gyy(i) + gyy(i+1) ) + agyz = 0.5d0*( gyz(i) + gyz(i+1) ) + agzz = 0.5d0*( gzz(i) + gzz(i+1) ) + det = SPATIAL_DETERMINANT(agxx, agxy, agxz, \ + agyy, agyz, agzz) + call UpperMetric (uxx, uxy, uxz, uyy, uyz, uzz, & + det, agxx, agxy, agxz, agyy, agyz, agzz) + call eigenvalues(handle,& + D_UPW(rho), D_UPW(velx), D_UPW(vely), D_UPW(velz), & + D_UPW(eps), D_UPW(w_lorentz), lam, & + agxx, agxy, agxz, agyy, agyz, agzz, & + uxx, D_UPW(alp), D_UPW(beta)) + l_ev_l(i)=lam(1) + l_ev_r(i)=lam(5) + xwind(i) = (lam(1) + lam(5)) / (abs(lam(1)) + abs(lam(5))) + xwind(i) = min(1.d0, max(-1.d0, xwind(i))) #define LEFTPLUS(x,xplus) xplus(i) = abs(xwind(i)) * LEFT1(x) + \ - (1.d0-abs(xwind(i))) * xplus(i) + (1.d0-abs(xwind(i))) * xplus(i) #define LEFTMINUS(x,xminus) xminus(i+1)= abs(xwind(i)) * LEFT1(x) + \ - (1.d0-abs(xwind(i))) * xminus(i+1) + (1.d0-abs(xwind(i))) * xminus(i+1) #define RIGHTPLUS(x,xplus) xplus(i) = abs(xwind(i)) * RIGHT1(x) + \ - (1.d0-abs(xwind(i))) * xplus(i) + (1.d0-abs(xwind(i))) * xplus(i) #define RIGHTMINUS(x,xminus) xminus(i+1)= abs(xwind(i)) * RIGHT1(x) + \ - (1.d0-abs(xwind(i))) * xminus(i+1) + (1.d0-abs(xwind(i))) * xminus(i+1) #define CHECK(x,xc) if (x(i+1) .gt. x(i)) then && xc=min(x(i+1),max(x(i),xc)) && else && xc=min(x(i),max(x(i+1),xc)) && endif !!$ xwind(i)=0.d0 - if (xwind(i) .lt. 0.0d0) then - LEFTPLUS(rho, rhoplus) - LEFTMINUS(rho, rhominus) - LEFTPLUS(velx, velxplus) - LEFTMINUS(velx, velxminus) - LEFTPLUS(vely, velyplus) - LEFTMINUS(vely, velyminus) - LEFTPLUS(velz, velzplus) - LEFTMINUS(velz, velzminus) - if (poly .eq. 0) then - LEFTPLUS(eps, epsplus) - LEFTMINUS(eps, epsminus) - end if - else - RIGHTPLUS(rho, rhoplus) - RIGHTMINUS(rho, rhominus) - RIGHTPLUS(velx, velxplus) - RIGHTMINUS(velx, velxminus) - RIGHTPLUS(vely, velyplus) - RIGHTMINUS(vely, velyminus) - RIGHTPLUS(velz, velzplus) - RIGHTMINUS(velz, velzminus) - if (poly .eq. 0) then - RIGHTPLUS(eps, epsplus) - RIGHTMINUS(eps, epsminus) - end if - end if - CHECK(rho, rhoplus(i)) - CHECK(rho, rhominus(i+1)) - CHECK(velx, velxplus(i)) - CHECK(velx, velxminus(i+1)) - CHECK(vely, velyplus(i)) - CHECK(vely, velyminus(i+1)) - CHECK(velz, velzplus(i)) - CHECK(velz, velzminus(i+1)) - if (poly .eq. 0) then - CHECK(eps, epsplus(i)) - CHECK(eps, epsminus(i+1)) - end if + if (xwind(i) .lt. 0.0d0) then + LEFTPLUS(rho, rhoplus) + LEFTMINUS(rho, rhominus) + LEFTPLUS(velx, velxplus) + LEFTMINUS(velx, velxminus) + LEFTPLUS(vely, velyplus) + LEFTMINUS(vely, velyminus) + LEFTPLUS(velz, velzplus) + LEFTMINUS(velz, velzminus) + if (poly .eq. 0) then + LEFTPLUS(eps, epsplus) + LEFTMINUS(eps, epsminus) + end if + else + RIGHTPLUS(rho, rhoplus) + RIGHTMINUS(rho, rhominus) + RIGHTPLUS(velx, velxplus) + RIGHTMINUS(velx, velxminus) + RIGHTPLUS(vely, velyplus) + RIGHTMINUS(vely, velyminus) + RIGHTPLUS(velz, velzplus) + RIGHTMINUS(velz, velzminus) + if (poly .eq. 0) then + RIGHTPLUS(eps, epsplus) + RIGHTMINUS(eps, epsminus) + end if + end if + CHECK(rho, rhoplus(i)) + CHECK(rho, rhominus(i+1)) + CHECK(velx, velxplus(i)) + CHECK(velx, velxminus(i+1)) + CHECK(vely, velyplus(i)) + CHECK(vely, velyminus(i+1)) + CHECK(velz, velzplus(i)) + CHECK(velz, velzminus(i+1)) + if (poly .eq. 0) then + CHECK(eps, epsplus(i)) + CHECK(eps, epsminus(i+1)) + end if !!$ if ((dir .eq. 1) .and. (ni .eq. 4) .and. (nj .eq. 4)) then !!$ write (*,*) rhoplus(i), rhominus(i+1) !!$ end if end do !!$ mppm debug output if (ppm_mppm_debug_eigenvalues .gt. 0) then - if (dir .eq. 1) then - ev_l(:,ni,nj) = l_ev_l - ev_r(:,ni,nj) = l_ev_r - xw(:,ni,nj) = xwind - else if (dir .eq. 2) then - ev_l(ni,:,nj) = l_ev_l - ev_r(ni,:,nj) = l_ev_r - xw(ni,:,nj) = xwind - else if (dir .eq. 3) then - ev_l(ni,nj,:) = l_ev_l - ev_r(ni,nj,:) = l_ev_r - xw(ni,nj,:) = xwind - else - write (*,*) "flux direction not 1 to 3 ?" - end if + if (dir .eq. 1) then + ev_l(:,ni,nj) = l_ev_l + ev_r(:,ni,nj) = l_ev_r + xw(:,ni,nj) = xwind + else if (dir .eq. 2) then + ev_l(ni,:,nj) = l_ev_l + ev_r(ni,:,nj) = l_ev_r + xw(ni,:,nj) = xwind + else if (dir .eq. 3) then + ev_l(ni,nj,:) = l_ev_l + ev_r(ni,nj,:) = l_ev_r + xw(ni,nj,:) = xwind + else + write (*,*) "flux direction not 1 to 3 ?" + end if end if end if end if @@ -481,46 +481,46 @@ trivial_rp = .true. if (use_enhanced_ppm .eq. 0) then ! In 1984 PPM, flattening is applied before constraining parabolic profiles. if (PPM3) then !!$ Implement C&W, page 197, but with a workaround which allows to use stencil=3. - do i = 3, nx - 2 - flatten = tilde_flatten(i) - if (abs(1.d0 - flatten) > 0.d0) then - trivial_rp(i-1) = .false. - trivial_rp(i) = .false. - end if - rhoplus(i) = flatten * rhoplus(i) + (1.d0 - flatten) * rho(i) - rhominus(i) = flatten * rhominus(i) + (1.d0 - flatten) * rho(i) - velxplus(i) = flatten * velxplus(i) + (1.d0 - flatten) * vx(i) - velxminus(i) = flatten * velxminus(i) + (1.d0 - flatten) * vx(i) - velyplus(i) = flatten * velyplus(i) + (1.d0 - flatten) * vy(i) - velyminus(i) = flatten * velyminus(i) + (1.d0 - flatten) * vy(i) - velzplus(i) = flatten * velzplus(i) + (1.d0 - flatten) * vz(i) - velzminus(i) = flatten * velzminus(i) + (1.d0 - flatten) * vz(i) - if (poly .eq. 0) then - epsplus(i) = flatten * epsplus(i) + (1.d0 - flatten) * eps(i) - epsminus(i) = flatten * epsminus(i) + (1.d0 - flatten) * eps(i) - end if - end do + do i = 3, nx - 2 + flatten = tilde_flatten(i) + if (abs(1.d0 - flatten) > 0.d0) then + trivial_rp(i-1) = .false. + trivial_rp(i) = .false. + end if + rhoplus(i) = flatten * rhoplus(i) + (1.d0 - flatten) * rho(i) + rhominus(i) = flatten * rhominus(i) + (1.d0 - flatten) * rho(i) + velxplus(i) = flatten * velxplus(i) + (1.d0 - flatten) * vx(i) + velxminus(i) = flatten * velxminus(i) + (1.d0 - flatten) * vx(i) + velyplus(i) = flatten * velyplus(i) + (1.d0 - flatten) * vy(i) + velyminus(i) = flatten * velyminus(i) + (1.d0 - flatten) * vy(i) + velzplus(i) = flatten * velzplus(i) + (1.d0 - flatten) * vz(i) + velzminus(i) = flatten * velzminus(i) + (1.d0 - flatten) * vz(i) + if (poly .eq. 0) then + epsplus(i) = flatten * epsplus(i) + (1.d0 - flatten) * eps(i) + epsminus(i) = flatten * epsminus(i) + (1.d0 - flatten) * eps(i) + end if + end do else !!$ Really implement C&W, page 197; which requires stencil 4. - do i = 4, nx - 3 - s=sign(one, -dpress(i)) - flatten = max(tilde_flatten(i), tilde_flatten(i+s)) - if (abs(1.d0 - flatten) > 0.d0) then - trivial_rp(i-1) = .false. - trivial_rp(i) = .false. - end if - rhoplus(i) = flatten * rhoplus(i) + (1.d0 - flatten) * rho(i) - rhominus(i) = flatten * rhominus(i) + (1.d0 - flatten) * rho(i) - velxplus(i) = flatten * velxplus(i) + (1.d0 - flatten) * vx(i) - velxminus(i) = flatten * velxminus(i) + (1.d0 - flatten) * vx(i) - velyplus(i) = flatten * velyplus(i) + (1.d0 - flatten) * vy(i) - velyminus(i) = flatten * velyminus(i) + (1.d0 - flatten) * vy(i) - velzplus(i) = flatten * velzplus(i) + (1.d0 - flatten) * vz(i) - velzminus(i) = flatten * velzminus(i) + (1.d0 - flatten) * vz(i) - if (poly .eq. 0) then - epsplus(i) = flatten * epsplus(i) + (1.d0 - flatten) * eps(i) - epsminus(i) = flatten * epsminus(i) + (1.d0 - flatten) * eps(i) - end if - end do + do i = 4, nx - 3 + s=sign(one, -dpress(i)) + flatten = max(tilde_flatten(i), tilde_flatten(i+s)) + if (abs(1.d0 - flatten) > 0.d0) then + trivial_rp(i-1) = .false. + trivial_rp(i) = .false. + end if + rhoplus(i) = flatten * rhoplus(i) + (1.d0 - flatten) * rho(i) + rhominus(i) = flatten * rhominus(i) + (1.d0 - flatten) * rho(i) + velxplus(i) = flatten * velxplus(i) + (1.d0 - flatten) * vx(i) + velxminus(i) = flatten * velxminus(i) + (1.d0 - flatten) * vx(i) + velyplus(i) = flatten * velyplus(i) + (1.d0 - flatten) * vy(i) + velyminus(i) = flatten * velyminus(i) + (1.d0 - flatten) * vy(i) + velzplus(i) = flatten * velzplus(i) + (1.d0 - flatten) * vz(i) + velzminus(i) = flatten * velzminus(i) + (1.d0 - flatten) * vz(i) + if (poly .eq. 0) then + epsplus(i) = flatten * epsplus(i) + (1.d0 - flatten) * eps(i) + epsminus(i) = flatten * epsminus(i) + (1.d0 - flatten) * eps(i) + end if + end do end if endif @@ -557,50 +557,50 @@ trivial_rp = .true. daplus = aplus(i)-a(i) &&\ daminus = a(i)-aminus(i) &&\ if (daplus*daminus .le. 0 .or. (a(i-2)-a(i))*(a(i)-a(i+2)) .le. 0) then &&\ - D2a = - (12.0d0*a(i) - 6.0d0*(aminus(i)+aplus(i))) && \ - D2aC = a(i-1) - 2.0d0*a(i) + a(i+1) &&\ - D2aL = a(i-2) - 2.0d0*a(i-1) + a(i) &&\ - D2aLL = a(i-3) - 2.0d0*a(i-2) + a(i-1) &&\ - D2aR = a(i) - 2.0d0*a(i+1) + a(i+2) &&\ - D2aRR = a(i+1) - 2.0d0*a(i+2) + a(i+3) &&\ - D3a = D2aR - D2aC &&\ - D3aL = D2aC - D2aL &&\ - D3aR = D2aRR - D2aR &&\ - D3aLL = D2aL - D2aLL &&\ - if (sign(one, D2a) .eq. sign(one, D2aC) .and. sign(one, D2a) .eq. sign(one, D2aL) .and. sign(one, D2a) .eq. sign(one, D2aR)) then &&\ - D2aLim = sign(one, D2a) * min(C*abs(D2aL), C*abs(D2aR), C*abs(D2aC), abs(D2a)) &&\ - else &&\ - D2aLim = 0 &&\ - end if &&\ - if (abs(D2a) .le. 1.d-12*max(abs(a(i-2)), abs(a(i-1)), abs(a(i)), abs(a(i+1)), abs(a(i+2)))) then &&\ - rhi = 0 &&\ - else &&\ - rhi = D2aLim / D2a &&\ - endif &&\ - if (.not. (rhi .ge. 1.0d0 - 1.d-12)) then &&\ - D3aMin = min(D3aLL, D3aL, D3a, D3aR) &&\ - D3aMax = max(D3aLL, D3aL, D3a, D3aR) &&\ - if (C3 * max(abs(D3aMin), abs(D3aMax)) .le. D3aMax-D3aMin) then &&\ - if (daplus*daminus .lt. 0) then &&\ - aplus(i) = a(i) + daplus * rhi &&\ - aminus(i) = a(i) - daminus * rhi &&\ - else if (abs(daminus) .ge. 2.0d0*abs(daplus)) then &&\ - aminus(i) = a(i) - 2.0d0*(1.0d0-rhi)*daplus - rhi*daminus &&\ - else if (abs(daplus) .ge. 2.0d0*abs(daminus)) then &&\ - aplus(i) = a(i) + 2.0d0*(1.0d0-rhi)*daminus + rhi*daplus &&\ - endif &&\ - endif &&\ - endif &&\ - trivial_rp(i-1) = .false. &&\ - trivial_rp(i) = .false. &&\ + D2a = - (12.0d0*a(i) - 6.0d0*(aminus(i)+aplus(i))) && \ + D2aC = a(i-1) - 2.0d0*a(i) + a(i+1) &&\ + D2aL = a(i-2) - 2.0d0*a(i-1) + a(i) &&\ + D2aLL = a(i-3) - 2.0d0*a(i-2) + a(i-1) &&\ + D2aR = a(i) - 2.0d0*a(i+1) + a(i+2) &&\ + D2aRR = a(i+1) - 2.0d0*a(i+2) + a(i+3) &&\ + D3a = D2aR - D2aC &&\ + D3aL = D2aC - D2aL &&\ + D3aR = D2aRR - D2aR &&\ + D3aLL = D2aL - D2aLL &&\ + if (sign(one, D2a) .eq. sign(one, D2aC) .and. sign(one, D2a) .eq. sign(one, D2aL) .and. sign(one, D2a) .eq. sign(one, D2aR)) then &&\ + D2aLim = sign(one, D2a) * min(C*abs(D2aL), C*abs(D2aR), C*abs(D2aC), abs(D2a)) &&\ + else &&\ + D2aLim = 0 &&\ + end if &&\ + if (abs(D2a) .le. 1.d-12*max(abs(a(i-2)), abs(a(i-1)), abs(a(i)), abs(a(i+1)), abs(a(i+2)))) then &&\ + rhi = 0 &&\ + else &&\ + rhi = D2aLim / D2a &&\ + endif &&\ + if (.not. (rhi .ge. 1.0d0 - 1.d-12)) then &&\ + D3aMin = min(D3aLL, D3aL, D3a, D3aR) &&\ + D3aMax = max(D3aLL, D3aL, D3a, D3aR) &&\ + if (C3 * max(abs(D3aMin), abs(D3aMax)) .le. D3aMax-D3aMin) then &&\ + if (daplus*daminus .lt. 0) then &&\ + aplus(i) = a(i) + daplus * rhi &&\ + aminus(i) = a(i) - daminus * rhi &&\ + else if (abs(daminus) .ge. 2.0d0*abs(daplus)) then &&\ + aminus(i) = a(i) - 2.0d0*(1.0d0-rhi)*daplus - rhi*daminus &&\ + else if (abs(daplus) .ge. 2.0d0*abs(daminus)) then &&\ + aplus(i) = a(i) + 2.0d0*(1.0d0-rhi)*daminus + rhi*daplus &&\ + endif &&\ + endif &&\ + endif &&\ + trivial_rp(i-1) = .false. &&\ + trivial_rp(i) = .false. &&\ else &&\ - if (abs(daplus) .ge. 2.0d0*abs(daminus)) then &&\ - aplus(i) = a(i) + 2.0d0*daminus &&\ - else if (abs(daminus) .ge. 2.0d0*abs(daplus)) then &&\ - aminus(i) = a(i) - 2.0d0*daplus &&\ - end if &&\ - trivial_rp(i-1) = .false. &&\ - trivial_rp(i) = .false. &&\ + if (abs(daplus) .ge. 2.0d0*abs(daminus)) then &&\ + aplus(i) = a(i) + 2.0d0*daminus &&\ + else if (abs(daminus) .ge. 2.0d0*abs(daplus)) then &&\ + aminus(i) = a(i) - 2.0d0*daplus &&\ + end if &&\ + trivial_rp(i-1) = .false. &&\ + trivial_rp(i) = .false. &&\ endif !! Monotonicity of PPM 2011 (McCorquodale & Colella 2011), Sec. 2.4.1, Eq. 23-34. @@ -611,69 +611,69 @@ trivial_rp = .true. daplus = aplus(i)-a(i) &&\ daminus = a(i)-aminus(i) &&\ if (daplus*daminus .le. 0 .or. (a(i-2)-a(i))*(a(i)-a(i+2)) .le. 0) then &&\ - D2a = - (12.0d0*a(i) - 6.0d0*(aminus(i)+aplus(i))) && \ - D2aC = a(i-1) - 2.0d0*a(i) + a(i+1) &&\ - D2aL = a(i-2) - 2.0d0*a(i-1) + a(i) &&\ - D2aLL = a(i-3) - 2.0d0*a(i-2) + a(i-1) &&\ - D2aR = a(i) - 2.0d0*a(i+1) + a(i+2) &&\ - D2aRR = a(i+1) - 2.0d0*a(i+2) + a(i+3) &&\ - D3a = D2aR - D2aC &&\ - D3aL = D2aC - D2aL &&\ - D3aR = D2aRR - D2aR &&\ - D3aLL = D2aL - D2aLL &&\ - if (sign(one, D2a) .eq. sign(one, D2aC) .and. sign(one, D2a) .eq. sign(one, D2aL) .and. sign(one, D2a) .eq. sign(one, D2aR)) then &&\ - D2aLim = sign(one, D2a) * min(C*abs(D2aL), C*abs(D2aR), C*abs(D2aC), abs(D2a)) &&\ - else &&\ - D2aLim = 0 &&\ - end if &&\ - if (abs(D2a) .le. 1.d-12*max(abs(a(i-2)), abs(a(i-1)), abs(a(i)), abs(a(i+1)), abs(a(i+2)))) then &&\ - rhi = 0 &&\ - else &&\ - rhi = D2aLim / D2a &&\ - endif &&\ - if (.not. (rhi .ge. 1.0d0 - 1.d-12)) then &&\ - D3aMin = min(D3aLL, D3aL, D3a, D3aR) &&\ - D3aMax = max(D3aLL, D3aL, D3a, D3aR) &&\ - if (C3 * max(abs(D3aMin), abs(D3aMax)) .le. D3aMax-D3aMin) then &&\ - if (abs(daplus) .le. abs(a(i)) .and. abs(daminus) .le. abs(a(i))) then &&\ - if (daplus*daminus .lt. 0) then &&\ - aplus(i) = a(i) + daplus * rhi &&\ - aminus(i) = a(i) - daminus * rhi &&\ - else if (abs(daminus) .ge. 2.0d0*abs(daplus)) then &&\ - aminus(i) = a(i) - 2.0d0*(1.0d0-rhi)*daplus - rhi*daminus &&\ - else if (abs(daplus) .ge. 2.0d0*abs(daminus)) then &&\ - aplus(i) = a(i) + 2.0d0*(1.0d0-rhi)*daminus + rhi*daplus &&\ - endif &&\ - else &&\ - if (daplus*daminus .lt. 0) then &&\ - aplus(i) = a(i) &&\ - aminus(i) = a(i) &&\ - else if (abs(daminus) .ge. 2.0d0*abs(daplus)) then &&\ - aminus(i) = a(i) &&\ - else if (abs(daplus) .ge. 2.0d0*abs(daminus)) then &&\ - aplus(i) = a(i) &&\ - endif &&\ - endif &&\ - endif &&\ - endif &&\ - trivial_rp(i-1) = .false. &&\ - trivial_rp(i) = .false. &&\ + D2a = - (12.0d0*a(i) - 6.0d0*(aminus(i)+aplus(i))) && \ + D2aC = a(i-1) - 2.0d0*a(i) + a(i+1) &&\ + D2aL = a(i-2) - 2.0d0*a(i-1) + a(i) &&\ + D2aLL = a(i-3) - 2.0d0*a(i-2) + a(i-1) &&\ + D2aR = a(i) - 2.0d0*a(i+1) + a(i+2) &&\ + D2aRR = a(i+1) - 2.0d0*a(i+2) + a(i+3) &&\ + D3a = D2aR - D2aC &&\ + D3aL = D2aC - D2aL &&\ + D3aR = D2aRR - D2aR &&\ + D3aLL = D2aL - D2aLL &&\ + if (sign(one, D2a) .eq. sign(one, D2aC) .and. sign(one, D2a) .eq. sign(one, D2aL) .and. sign(one, D2a) .eq. sign(one, D2aR)) then &&\ + D2aLim = sign(one, D2a) * min(C*abs(D2aL), C*abs(D2aR), C*abs(D2aC), abs(D2a)) &&\ + else &&\ + D2aLim = 0 &&\ + end if &&\ + if (abs(D2a) .le. 1.d-12*max(abs(a(i-2)), abs(a(i-1)), abs(a(i)), abs(a(i+1)), abs(a(i+2)))) then &&\ + rhi = 0 &&\ + else &&\ + rhi = D2aLim / D2a &&\ + endif &&\ + if (.not. (rhi .ge. 1.0d0 - 1.d-12)) then &&\ + D3aMin = min(D3aLL, D3aL, D3a, D3aR) &&\ + D3aMax = max(D3aLL, D3aL, D3a, D3aR) &&\ + if (C3 * max(abs(D3aMin), abs(D3aMax)) .le. D3aMax-D3aMin) then &&\ + if (abs(daplus) .le. abs(a(i)) .and. abs(daminus) .le. abs(a(i))) then &&\ + if (daplus*daminus .lt. 0) then &&\ + aplus(i) = a(i) + daplus * rhi &&\ + aminus(i) = a(i) - daminus * rhi &&\ + else if (abs(daminus) .ge. 2.0d0*abs(daplus)) then &&\ + aminus(i) = a(i) - 2.0d0*(1.0d0-rhi)*daplus - rhi*daminus &&\ + else if (abs(daplus) .ge. 2.0d0*abs(daminus)) then &&\ + aplus(i) = a(i) + 2.0d0*(1.0d0-rhi)*daminus + rhi*daplus &&\ + endif &&\ + else &&\ + if (daplus*daminus .lt. 0) then &&\ + aplus(i) = a(i) &&\ + aminus(i) = a(i) &&\ + else if (abs(daminus) .ge. 2.0d0*abs(daplus)) then &&\ + aminus(i) = a(i) &&\ + else if (abs(daplus) .ge. 2.0d0*abs(daminus)) then &&\ + aplus(i) = a(i) &&\ + endif &&\ + endif &&\ + endif &&\ + endif &&\ + trivial_rp(i-1) = .false. &&\ + trivial_rp(i) = .false. &&\ else &&\ - if (abs(daplus) .le. abs(a(i)) .and. abs(daminus) .le. abs(a(i))) then &&\ - if (abs(daplus) .ge. 2.0d0*abs(daminus)) then &&\ - aplus(i) = a(i) + 2.0d0*daminus &&\ - else if (abs(daminus) .ge. 2.0d0*abs(daplus)) then &&\ - aminus(i) = a(i) - 2.0d0*daplus &&\ - end if &&\ - else &&\ - if (abs(daplus) .ge. 2.0d0*abs(daminus)) then &&\ - aplus(i) = a(i) &&\ - else if (abs(daminus) .ge. 2.0d0*abs(daplus)) then &&\ - aminus(i) = a(i) &&\ - end if &&\ - endif &&\ - trivial_rp(i-1) = .false. &&\ - trivial_rp(i) = .false. &&\ + if (abs(daplus) .le. abs(a(i)) .and. abs(daminus) .le. abs(a(i))) then &&\ + if (abs(daplus) .ge. 2.0d0*abs(daminus)) then &&\ + aplus(i) = a(i) + 2.0d0*daminus &&\ + else if (abs(daminus) .ge. 2.0d0*abs(daplus)) then &&\ + aminus(i) = a(i) - 2.0d0*daplus &&\ + end if &&\ + else &&\ + if (abs(daplus) .ge. 2.0d0*abs(daminus)) then &&\ + aplus(i) = a(i) &&\ + else if (abs(daminus) .ge. 2.0d0*abs(daplus)) then &&\ + aminus(i) = a(i) &&\ + end if &&\ + endif &&\ + trivial_rp(i-1) = .false. &&\ + trivial_rp(i) = .false. &&\ endif @@ -684,40 +684,40 @@ trivial_rp = .true. daplus = aplus(i)-a(i) &&\ daminus = a(i)-aminus(i) &&\ if (daplus*daminus .le. 0 .or. (a(i-2)-a(i))*(a(i)-a(i+2)) .le. 0) then &&\ - D2a = - (12.0d0*a(i) - 6.0d0*(aminus(i)+aplus(i))) && \ - D2aC = a(i-1) - 2.0d0*a(i) + a(i+1) &&\ - D2aL = a(i-2) - 2.0d0*a(i-1) + a(i) &&\ - D2aR = a(i) - 2.0d0*a(i+1) + a(i+2) &&\ - if (sign(one, D2a) .eq. sign(one, D2aC) .and. sign(one, D2a) .eq. sign(one, D2aL) .and. sign(one, D2a) .eq. sign(one, D2aR)) then &&\ - D2aLim = sign(one, D2a) * min(C*abs(D2aL), C*abs(D2aR), C*abs(D2aC), abs(D2a)) &&\ - else &&\ - D2aLim = 0 &&\ - end if &&\ - if (abs(D2a) .le. 1.d-12*max(abs(a(i-2)), abs(a(i-1)), abs(a(i)), abs(a(i+1)), abs(a(i+2)))) then &&\ - rhi = 0 &&\ - else &&\ - rhi = D2aLim / D2a &&\ - endif &&\ - if (.not. (rhi .ge. 1.0d0 - 1.d-12)) then &&\ - if (daplus*daminus .lt. 0) then &&\ - aplus(i) = a(i) + daplus * rhi &&\ - aminus(i) = a(i) - daminus * rhi &&\ - else if (abs(daminus) .ge. 2.0d0*abs(daplus)) then &&\ - aminus(i) = a(i) - 2.0d0*(1.0d0-rhi)*daplus - rhi*daminus &&\ - else if (abs(daplus) .ge. 2.0d0*abs(daminus)) then &&\ - aplus(i) = a(i) + 2.0d0*(1.0d0-rhi)*daminus + rhi*daplus &&\ - endif &&\ - endif &&\ - trivial_rp(i-1) = .false. &&\ - trivial_rp(i) = .false. &&\ + D2a = - (12.0d0*a(i) - 6.0d0*(aminus(i)+aplus(i))) && \ + D2aC = a(i-1) - 2.0d0*a(i) + a(i+1) &&\ + D2aL = a(i-2) - 2.0d0*a(i-1) + a(i) &&\ + D2aR = a(i) - 2.0d0*a(i+1) + a(i+2) &&\ + if (sign(one, D2a) .eq. sign(one, D2aC) .and. sign(one, D2a) .eq. sign(one, D2aL) .and. sign(one, D2a) .eq. sign(one, D2aR)) then &&\ + D2aLim = sign(one, D2a) * min(C*abs(D2aL), C*abs(D2aR), C*abs(D2aC), abs(D2a)) &&\ + else &&\ + D2aLim = 0 &&\ + end if &&\ + if (abs(D2a) .le. 1.d-12*max(abs(a(i-2)), abs(a(i-1)), abs(a(i)), abs(a(i+1)), abs(a(i+2)))) then &&\ + rhi = 0 &&\ + else &&\ + rhi = D2aLim / D2a &&\ + endif &&\ + if (.not. (rhi .ge. 1.0d0 - 1.d-12)) then &&\ + if (daplus*daminus .lt. 0) then &&\ + aplus(i) = a(i) + daplus * rhi &&\ + aminus(i) = a(i) - daminus * rhi &&\ + else if (abs(daminus) .ge. 2.0d0*abs(daplus)) then &&\ + aminus(i) = a(i) - 2.0d0*(1.0d0-rhi)*daplus - rhi*daminus &&\ + else if (abs(daplus) .ge. 2.0d0*abs(daminus)) then &&\ + aplus(i) = a(i) + 2.0d0*(1.0d0-rhi)*daminus + rhi*daplus &&\ + endif &&\ + endif &&\ + trivial_rp(i-1) = .false. &&\ + trivial_rp(i) = .false. &&\ else &&\ - if (abs(daplus) .ge. 2.0d0*abs(daminus)) then &&\ - aplus(i) = a(i) + 2.0d0*daminus &&\ - else if (abs(daminus) .ge. 2.0d0*abs(daplus)) then &&\ - aminus(i) = a(i) - 2.0d0*daplus &&\ - end if &&\ - trivial_rp(i-1) = .false. &&\ - trivial_rp(i) = .false. &&\ + if (abs(daplus) .ge. 2.0d0*abs(daminus)) then &&\ + aplus(i) = a(i) + 2.0d0*daminus &&\ + else if (abs(daminus) .ge. 2.0d0*abs(daplus)) then &&\ + aminus(i) = a(i) - 2.0d0*daplus &&\ + end if &&\ + trivial_rp(i-1) = .false. &&\ + trivial_rp(i) = .false. &&\ endif @@ -729,133 +729,133 @@ trivial_rp = .true. daplus = aplus(i)-a(i) &&\ daminus = a(i)-aminus(i) &&\ if (daplus*daminus .le. 0 .or. (a(i-2)-a(i))*(a(i)-a(i+2)) .le. 0) then &&\ - D2a = - (12.0d0*a(i) - 6.0d0*(aminus(i)+aplus(i))) && \ - D2aC = a(i-1) - 2.0d0*a(i) + a(i+1) &&\ - D2aL = a(i-2) - 2.0d0*a(i-1) + a(i) &&\ - D2aR = a(i) - 2.0d0*a(i+1) + a(i+2) &&\ - if (sign(one, D2a) .eq. sign(one, D2aC) .and. sign(one, D2a) .eq. sign(one, D2aL) .and. sign(one, D2a) .eq. sign(one, D2aR)) then &&\ - D2aLim = sign(one, D2a) * min(C*abs(D2aL), C*abs(D2aR), C*abs(D2aC), abs(D2a)) &&\ - else &&\ - D2aLim = 0 &&\ - end if &&\ - if (abs(D2a) .le. 1.d-12*max(abs(a(i-2)), abs(a(i-1)), abs(a(i)), abs(a(i+1)), abs(a(i+2)))) then &&\ - rhi = 0 &&\ - else &&\ - rhi = D2aLim / D2a &&\ - endif &&\ - if (.not. (rhi .ge. 1.0d0 - 1.d-12)) then &&\ - if (abs(daplus) .le. abs(a(i)) .and. abs(daminus) .le. abs(a(i))) then &&\ - if (daplus*daminus .lt. 0) then &&\ - aplus(i) = a(i) + daplus * rhi &&\ - aminus(i) = a(i) - daminus * rhi &&\ - else if (abs(daminus) .ge. 2.0d0*abs(daplus)) then &&\ - aminus(i) = a(i) - 2.0d0*(1.0d0-rhi)*daplus - rhi*daminus &&\ - else if (abs(daplus) .ge. 2.0d0*abs(daminus)) then &&\ - aplus(i) = a(i) + 2.0d0*(1.0d0-rhi)*daminus + rhi*daplus &&\ - endif &&\ - else &&\ - if (daplus*daminus .lt. 0) then &&\ - aplus(i) = a(i) &&\ - aminus(i) = a(i) &&\ - else if (abs(daminus) .ge. 2.0d0*abs(daplus)) then &&\ - aminus(i) = a(i) &&\ - else if (abs(daplus) .ge. 2.0d0*abs(daminus)) then &&\ - aplus(i) = a(i) &&\ - endif &&\ - endif &&\ - endif &&\ - trivial_rp(i-1) = .false. &&\ - trivial_rp(i) = .false. &&\ + D2a = - (12.0d0*a(i) - 6.0d0*(aminus(i)+aplus(i))) && \ + D2aC = a(i-1) - 2.0d0*a(i) + a(i+1) &&\ + D2aL = a(i-2) - 2.0d0*a(i-1) + a(i) &&\ + D2aR = a(i) - 2.0d0*a(i+1) + a(i+2) &&\ + if (sign(one, D2a) .eq. sign(one, D2aC) .and. sign(one, D2a) .eq. sign(one, D2aL) .and. sign(one, D2a) .eq. sign(one, D2aR)) then &&\ + D2aLim = sign(one, D2a) * min(C*abs(D2aL), C*abs(D2aR), C*abs(D2aC), abs(D2a)) &&\ + else &&\ + D2aLim = 0 &&\ + end if &&\ + if (abs(D2a) .le. 1.d-12*max(abs(a(i-2)), abs(a(i-1)), abs(a(i)), abs(a(i+1)), abs(a(i+2)))) then &&\ + rhi = 0 &&\ + else &&\ + rhi = D2aLim / D2a &&\ + endif &&\ + if (.not. (rhi .ge. 1.0d0 - 1.d-12)) then &&\ + if (abs(daplus) .le. abs(a(i)) .and. abs(daminus) .le. abs(a(i))) then &&\ + if (daplus*daminus .lt. 0) then &&\ + aplus(i) = a(i) + daplus * rhi &&\ + aminus(i) = a(i) - daminus * rhi &&\ + else if (abs(daminus) .ge. 2.0d0*abs(daplus)) then &&\ + aminus(i) = a(i) - 2.0d0*(1.0d0-rhi)*daplus - rhi*daminus &&\ + else if (abs(daplus) .ge. 2.0d0*abs(daminus)) then &&\ + aplus(i) = a(i) + 2.0d0*(1.0d0-rhi)*daminus + rhi*daplus &&\ + endif &&\ + else &&\ + if (daplus*daminus .lt. 0) then &&\ + aplus(i) = a(i) &&\ + aminus(i) = a(i) &&\ + else if (abs(daminus) .ge. 2.0d0*abs(daplus)) then &&\ + aminus(i) = a(i) &&\ + else if (abs(daplus) .ge. 2.0d0*abs(daminus)) then &&\ + aplus(i) = a(i) &&\ + endif &&\ + endif &&\ + endif &&\ + trivial_rp(i-1) = .false. &&\ + trivial_rp(i) = .false. &&\ else &&\ - if (abs(daplus) .le. abs(a(i)) .and. abs(daminus) .le. abs(a(i))) then &&\ - if (abs(daplus) .ge. 2.0d0*abs(daminus)) then &&\ - aplus(i) = a(i) + 2.0d0*daminus &&\ - else if (abs(daminus) .ge. 2.0d0*abs(daplus)) then &&\ - aminus(i) = a(i) - 2.0d0*daplus &&\ - end if &&\ - else &&\ - if (abs(daplus) .ge. 2.0d0*abs(daminus)) then &&\ - aplus(i) = a(i) &&\ - else if (abs(daminus) .ge. 2.0d0*abs(daplus)) then &&\ - aminus(i) = a(i) &&\ - end if &&\ - endif &&\ - trivial_rp(i-1) = .false. &&\ - trivial_rp(i) = .false. &&\ + if (abs(daplus) .le. abs(a(i)) .and. abs(daminus) .le. abs(a(i))) then &&\ + if (abs(daplus) .ge. 2.0d0*abs(daminus)) then &&\ + aplus(i) = a(i) + 2.0d0*daminus &&\ + else if (abs(daminus) .ge. 2.0d0*abs(daplus)) then &&\ + aminus(i) = a(i) - 2.0d0*daplus &&\ + end if &&\ + else &&\ + if (abs(daplus) .ge. 2.0d0*abs(daminus)) then &&\ + aplus(i) = a(i) &&\ + else if (abs(daminus) .ge. 2.0d0*abs(daplus)) then &&\ + aminus(i) = a(i) &&\ + end if &&\ + endif &&\ + trivial_rp(i-1) = .false. &&\ + trivial_rp(i) = .false. &&\ endif if (use_enhanced_ppm .eq. 1) then !! Constrain parabolic profiles, PPM 2011/2008 if (PPM3) then - do i = 3, nx - 2 - MON_WITH_LOCAL_EXTREMUM(rhominus,rho,rhoplus, enhanced_ppm_C2) - MON_WITH_LOCAL_EXTREMUM(velxminus,vx,velxplus, enhanced_ppm_C2) - MON_WITH_LOCAL_EXTREMUM(velyminus,vy,velyplus, enhanced_ppm_C2) - MON_WITH_LOCAL_EXTREMUM(velzminus,vz,velzplus, enhanced_ppm_C2) - end do - if (poly .eq. 0) then - do i = 3, nx - 2 - MON_WITH_LOCAL_EXTREMUM_EPS(epsminus,eps,epsplus, enhanced_ppm_C2) - end do - end if - + do i = 3, nx - 2 + MON_WITH_LOCAL_EXTREMUM(rhominus,rho,rhoplus, enhanced_ppm_C2) + MON_WITH_LOCAL_EXTREMUM(velxminus,vx,velxplus, enhanced_ppm_C2) + MON_WITH_LOCAL_EXTREMUM(velyminus,vy,velyplus, enhanced_ppm_C2) + MON_WITH_LOCAL_EXTREMUM(velzminus,vz,velzplus, enhanced_ppm_C2) + end do + if (poly .eq. 0) then + do i = 3, nx - 2 + MON_WITH_LOCAL_EXTREMUM_EPS(epsminus,eps,epsplus, enhanced_ppm_C2) + end do + end if + else - do i = 4, nx - 3 - MON_WITH_LOCAL_EXTREMUM_STENCIL4(rhominus,rho,rhoplus, enhanced_ppm_C2, enhanced_ppm_C3) - MON_WITH_LOCAL_EXTREMUM_STENCIL4(velxminus,vx,velxplus, enhanced_ppm_C2, enhanced_ppm_C3) - MON_WITH_LOCAL_EXTREMUM_STENCIL4(velyminus,vy,velyplus, enhanced_ppm_C2, enhanced_ppm_C3) - MON_WITH_LOCAL_EXTREMUM_STENCIL4(velzminus,vz,velzplus, enhanced_ppm_C2, enhanced_ppm_C3) - end do - if (poly .eq. 0) then - do i = 4, nx - 3 - MON_WITH_LOCAL_EXTREMUM_STENCIL4_EPS(epsminus,eps,epsplus, enhanced_ppm_C2, enhanced_ppm_C3) - end do - endif + do i = 4, nx - 3 + MON_WITH_LOCAL_EXTREMUM_STENCIL4(rhominus,rho,rhoplus, enhanced_ppm_C2, enhanced_ppm_C3) + MON_WITH_LOCAL_EXTREMUM_STENCIL4(velxminus,vx,velxplus, enhanced_ppm_C2, enhanced_ppm_C3) + MON_WITH_LOCAL_EXTREMUM_STENCIL4(velyminus,vy,velyplus, enhanced_ppm_C2, enhanced_ppm_C3) + MON_WITH_LOCAL_EXTREMUM_STENCIL4(velzminus,vz,velzplus, enhanced_ppm_C2, enhanced_ppm_C3) + end do + if (poly .eq. 0) then + do i = 4, nx - 3 + MON_WITH_LOCAL_EXTREMUM_STENCIL4_EPS(epsminus,eps,epsplus, enhanced_ppm_C2, enhanced_ppm_C3) + end do + endif end if ! Apply flattening after constraining the parabolic profiles if (PPM3) then !!$ Implement C&W, page 197, but with a workaround which allows to use stencil=3. - do i = 3, nx - 2 - flatten = tilde_flatten(i) - if (abs(1.d0 - flatten) > 0.d0) then - trivial_rp(i-1) = .false. - trivial_rp(i) = .false. - end if - rhoplus(i) = flatten * rhoplus(i) + (1.d0 - flatten) * rho(i) - rhominus(i) = flatten * rhominus(i) + (1.d0 - flatten) * rho(i) - velxplus(i) = flatten * velxplus(i) + (1.d0 - flatten) * vx(i) - velxminus(i) = flatten * velxminus(i) + (1.d0 - flatten) * vx(i) - velyplus(i) = flatten * velyplus(i) + (1.d0 - flatten) * vy(i) - velyminus(i) = flatten * velyminus(i) + (1.d0 - flatten) * vy(i) - velzplus(i) = flatten * velzplus(i) + (1.d0 - flatten) * vz(i) - velzminus(i) = flatten * velzminus(i) + (1.d0 - flatten) * vz(i) - if (poly .eq. 0) then - epsplus(i) = flatten * epsplus(i) + (1.d0 - flatten) * eps(i) - epsminus(i) = flatten * epsminus(i) + (1.d0 - flatten) * eps(i) - end if - end do + do i = 3, nx - 2 + flatten = tilde_flatten(i) + if (abs(1.d0 - flatten) > 0.d0) then + trivial_rp(i-1) = .false. + trivial_rp(i) = .false. + end if + rhoplus(i) = flatten * rhoplus(i) + (1.d0 - flatten) * rho(i) + rhominus(i) = flatten * rhominus(i) + (1.d0 - flatten) * rho(i) + velxplus(i) = flatten * velxplus(i) + (1.d0 - flatten) * vx(i) + velxminus(i) = flatten * velxminus(i) + (1.d0 - flatten) * vx(i) + velyplus(i) = flatten * velyplus(i) + (1.d0 - flatten) * vy(i) + velyminus(i) = flatten * velyminus(i) + (1.d0 - flatten) * vy(i) + velzplus(i) = flatten * velzplus(i) + (1.d0 - flatten) * vz(i) + velzminus(i) = flatten * velzminus(i) + (1.d0 - flatten) * vz(i) + if (poly .eq. 0) then + epsplus(i) = flatten * epsplus(i) + (1.d0 - flatten) * eps(i) + epsminus(i) = flatten * epsminus(i) + (1.d0 - flatten) * eps(i) + end if + end do else !!$ Really implement C&W, page 197; which requires stencil 4. - do i = 4, nx - 3 - s=sign(one, -dpress(i)) - flatten = max(tilde_flatten(i), tilde_flatten(i+s)) - if (abs(1.d0 - flatten) > 0.d0) then - trivial_rp(i-1) = .false. - trivial_rp(i) = .false. - end if - rhoplus(i) = flatten * rhoplus(i) + (1.d0 - flatten) * rho(i) - rhominus(i) = flatten * rhominus(i) + (1.d0 - flatten) * rho(i) - velxplus(i) = flatten * velxplus(i) + (1.d0 - flatten) * vx(i) - velxminus(i) = flatten * velxminus(i) + (1.d0 - flatten) * vx(i) - velyplus(i) = flatten * velyplus(i) + (1.d0 - flatten) * vy(i) - velyminus(i) = flatten * velyminus(i) + (1.d0 - flatten) * vy(i) - velzplus(i) = flatten * velzplus(i) + (1.d0 - flatten) * vz(i) - velzminus(i) = flatten * velzminus(i) + (1.d0 - flatten) * vz(i) - if (poly .eq. 0) then - epsplus(i) = flatten * epsplus(i) + (1.d0 - flatten) * eps(i) - epsminus(i) = flatten * epsminus(i) + (1.d0 - flatten) * eps(i) - end if - end do + do i = 4, nx - 3 + s=sign(one, -dpress(i)) + flatten = max(tilde_flatten(i), tilde_flatten(i+s)) + if (abs(1.d0 - flatten) > 0.d0) then + trivial_rp(i-1) = .false. + trivial_rp(i) = .false. + end if + rhoplus(i) = flatten * rhoplus(i) + (1.d0 - flatten) * rho(i) + rhominus(i) = flatten * rhominus(i) + (1.d0 - flatten) * rho(i) + velxplus(i) = flatten * velxplus(i) + (1.d0 - flatten) * vx(i) + velxminus(i) = flatten * velxminus(i) + (1.d0 - flatten) * vx(i) + velyplus(i) = flatten * velyplus(i) + (1.d0 - flatten) * vy(i) + velyminus(i) = flatten * velyminus(i) + (1.d0 - flatten) * vy(i) + velzplus(i) = flatten * velzplus(i) + (1.d0 - flatten) * vz(i) + velzminus(i) = flatten * velzminus(i) + (1.d0 - flatten) * vz(i) + if (poly .eq. 0) then + epsplus(i) = flatten * epsplus(i) + (1.d0 - flatten) * eps(i) + epsminus(i) = flatten * epsminus(i) + (1.d0 - flatten) * eps(i) + end if + end do end if endif @@ -863,16 +863,16 @@ trivial_rp = .true. if (use_enhanced_ppm .eq. 0) then !! Constrain parabolic profiles, PPM 1984 do i = GRHydro_stencil, nx - GRHydro_stencil + 1 - ! original Colella&Woodward monotonicity preservation - MON(rhominus,rho,rhoplus) - MON(velxminus,vx,velxplus) - MON(velyminus,vy,velyplus) - MON(velzminus,vz,velzplus) + ! original Colella&Woodward monotonicity preservation + MON(rhominus,rho,rhoplus) + MON(velxminus,vx,velxplus) + MON(velyminus,vy,velyplus) + MON(velzminus,vz,velzplus) end do if (poly .eq. 0) then - do i = GRHydro_stencil, nx - GRHydro_stencil + 1 - MON(epsminus,eps,epsplus) - end do + do i = GRHydro_stencil, nx - GRHydro_stencil + 1 + MON(epsminus,eps,epsplus) + end do end if endif @@ -1067,13 +1067,13 @@ subroutine SimplePPM_temperature_1d(& #define STEEP(x,dx,dmx) \ - if ( (x(i+1) - x(i)) * (x(i) - x(i-1)) > 0.d0 ) then &&\ - dmx(i) = sign(one, dx(i)) * \ - min(abs(dx(i)), 2.d0 * abs(x(i) - x(i-1)), \ - 2.d0 * abs(x(i+1) - x(i))) &&\ - else &&\ - dmx(i) = 0.d0 &&\ - end if + if ( (x(i+1) - x(i)) * (x(i) - x(i-1)) > 0.d0 ) then &&\ + dmx(i) = sign(one, dx(i)) * \ + min(abs(dx(i)), 2.d0 * abs(x(i) - x(i-1)), \ + 2.d0 * abs(x(i+1) - x(i))) &&\ + else &&\ + dmx(i) = 0.d0 &&\ + end if if (use_enhanced_ppm .eq. 0) then @@ -1117,32 +1117,32 @@ subroutine SimplePPM_temperature_1d(& #define LIMITT(a,ah,C, alim) \ if ((min(a(i),a(i+1)) .le. ah) .and. (ah .le. max(a(i),a(i+1)))) then &&\ - alim = ah &&\ + alim = ah &&\ else &&\ - D2a = 3.0d0 * (a(i) - 2.0d0*ah + a(i+1)) &&\ - D2aL = (a(i-1) - 2.0d0*a(i) + a(i+1)) &&\ - D2aR = (a(i) - 2.0d0*a(i+1) + a(i+2)) &&\ - D2aLim = sign(one, D2a)*min(C*abs(D2aL), C*abs(D2aR), abs(D2a))/3.0d0 &&\ - if (D2a*D2aR .ge. 0 .and. D2a*D2aL .ge. 0) then &&\ - alim = 0.5d0*(a(i)+a(i+1)) - D2aLim &&\ - else &&\ - alim = 0.5d0*(a(i)+a(i+1)) &&\ - end if &&\ + D2a = 3.0d0 * (a(i) - 2.0d0*ah + a(i+1)) &&\ + D2aL = (a(i-1) - 2.0d0*a(i) + a(i+1)) &&\ + D2aR = (a(i) - 2.0d0*a(i+1) + a(i+2)) &&\ + D2aLim = sign(one, D2a)*min(C*abs(D2aL), C*abs(D2aR), abs(D2a))/3.0d0 &&\ + if (D2a*D2aR .ge. 0 .and. D2a*D2aL .ge. 0) then &&\ + alim = 0.5d0*(a(i)+a(i+1)) - D2aLim &&\ + else &&\ + alim = 0.5d0*(a(i)+a(i+1)) &&\ + end if &&\ endif #define LIMITT_EPS(a,ah,C, alim) \ if ((min(a(i),a(i+1)) .le. ah) .and. (ah .le. max(a(i),a(i+1)))) then &&\ - alim = ah &&\ + alim = ah &&\ else &&\ - D2a = 3.0d0 * (a(i) - 2.0d0*ah + a(i+1)) &&\ - D2aL = (a(i-1) - 2.0d0*a(i) + a(i+1)) &&\ - D2aR = (a(i) - 2.0d0*a(i+1) + a(i+2)) &&\ - D2aLim = sign(one, D2a)*min(C*abs(D2aL), C*abs(D2aR), abs(D2a))/3.0d0 &&\ - if (D2a*D2aR .ge. 0 .and. D2a*D2aL .ge. 0 .and. abs(D2aLim) .lt. abs(0.5d0*(a(i)+a(i+1)))) then &&\ - alim = 0.5d0*(a(i)+a(i+1)) - D2aLim &&\ - else &&\ - alim = 0.5d0*(a(i)+a(i+1)) &&\ - end if&&\ + D2a = 3.0d0 * (a(i) - 2.0d0*ah + a(i+1)) &&\ + D2aL = (a(i-1) - 2.0d0*a(i) + a(i+1)) &&\ + D2aR = (a(i) - 2.0d0*a(i+1) + a(i+2)) &&\ + D2aLim = sign(one, D2a)*min(C*abs(D2aL), C*abs(D2aR), abs(D2a))/3.0d0 &&\ + if (D2a*D2aR .ge. 0 .and. D2a*D2aL .ge. 0 .and. abs(D2aLim) .lt. abs(0.5d0*(a(i)+a(i+1)))) then &&\ + alim = 0.5d0*(a(i)+a(i+1)) - D2aLim &&\ + else &&\ + alim = 0.5d0*(a(i)+a(i+1)) &&\ + end if&&\ endif @@ -1202,18 +1202,18 @@ subroutine SimplePPM_temperature_1d(& if (use_enhanced_ppm .eq. 0) then ! In 1984 PPM, flattening is applied before constraining parabolic profiles. if (PPM3) then !!$ Implement C&W, page 197, but with a workaround which allows to use stencil=3. - do i = 3, nx - 2 - flatten = tilde_flatten(i) + do i = 3, nx - 2 + flatten = tilde_flatten(i) tempplus(i) = flatten * tempplus(i) + (1.d0 - flatten) * temperature(i) tempminus(i) = flatten * tempminus(i) + (1.d0 - flatten) * temperature(i) - end do + end do else !!$ Really implement C&W, page 197; which requires stencil 4. - do i = 4, nx - 3 - s=sign(one, -dpress(i)) - flatten = max(tilde_flatten(i), tilde_flatten(i+s)) + do i = 4, nx - 3 + s=sign(one, -dpress(i)) + flatten = max(tilde_flatten(i), tilde_flatten(i+s)) tempplus(i) = flatten * tempplus(i) + (1.d0 - flatten) * temperature(i) tempminus(i) = flatten * tempminus(i) + (1.d0 - flatten) * temperature(i) - end do + end do end if endif @@ -1240,46 +1240,46 @@ subroutine SimplePPM_temperature_1d(& daplus = aplus(i)-a(i) &&\ daminus = a(i)-aminus(i) &&\ if (daplus*daminus .le. 0 .or. (a(i-2)-a(i))*(a(i)-a(i+2)) .le. 0) then &&\ - D2a = - (12.0d0*a(i) - 6.0d0*(aminus(i)+aplus(i))) && \ - D2aC = a(i-1) - 2.0d0*a(i) + a(i+1) &&\ - D2aL = a(i-2) - 2.0d0*a(i-1) + a(i) &&\ - D2aLL = a(i-3) - 2.0d0*a(i-2) + a(i-1) &&\ - D2aR = a(i) - 2.0d0*a(i+1) + a(i+2) &&\ - D2aRR = a(i+1) - 2.0d0*a(i+2) + a(i+3) &&\ - D3a = D2aR - D2aC &&\ - D3aL = D2aC - D2aL &&\ - D3aR = D2aRR - D2aR &&\ - D3aLL = D2aL - D2aLL &&\ - if (sign(one, D2a) .eq. sign(one, D2aC) .and. sign(one, D2a) .eq. sign(one, D2aL) .and. sign(one, D2a) .eq. sign(one, D2aR)) then &&\ - D2aLim = sign(one, D2a) * min(C*abs(D2aL), C*abs(D2aR), C*abs(D2aC), abs(D2a)) &&\ - else &&\ - D2aLim = 0 &&\ - end if &&\ - if (abs(D2a) .le. 1.d-12*max(abs(a(i-2)), abs(a(i-1)), abs(a(i)), abs(a(i+1)), abs(a(i+2)))) then &&\ - rhi = 0 &&\ - else &&\ - rhi = D2aLim / D2a &&\ - endif &&\ - if (.not. (rhi .ge. 1.0d0 - 1.d-12)) then &&\ - D3aMin = min(D3aLL, D3aL, D3a, D3aR) &&\ - D3aMax = max(D3aLL, D3aL, D3a, D3aR) &&\ - if (C3 * max(abs(D3aMin), abs(D3aMax)) .le. D3aMax-D3aMin) then &&\ - if (daplus*daminus .lt. 0) then &&\ - aplus(i) = a(i) + daplus * rhi &&\ - aminus(i) = a(i) - daminus * rhi &&\ - else if (abs(daminus) .ge. 2.0d0*abs(daplus)) then &&\ - aminus(i) = a(i) - 2.0d0*(1.0d0-rhi)*daplus - rhi*daminus &&\ - else if (abs(daplus) .ge. 2.0d0*abs(daminus)) then &&\ - aplus(i) = a(i) + 2.0d0*(1.0d0-rhi)*daminus + rhi*daplus &&\ - endif &&\ - endif &&\ - endif &&\ + D2a = - (12.0d0*a(i) - 6.0d0*(aminus(i)+aplus(i))) && \ + D2aC = a(i-1) - 2.0d0*a(i) + a(i+1) &&\ + D2aL = a(i-2) - 2.0d0*a(i-1) + a(i) &&\ + D2aLL = a(i-3) - 2.0d0*a(i-2) + a(i-1) &&\ + D2aR = a(i) - 2.0d0*a(i+1) + a(i+2) &&\ + D2aRR = a(i+1) - 2.0d0*a(i+2) + a(i+3) &&\ + D3a = D2aR - D2aC &&\ + D3aL = D2aC - D2aL &&\ + D3aR = D2aRR - D2aR &&\ + D3aLL = D2aL - D2aLL &&\ + if (sign(one, D2a) .eq. sign(one, D2aC) .and. sign(one, D2a) .eq. sign(one, D2aL) .and. sign(one, D2a) .eq. sign(one, D2aR)) then &&\ + D2aLim = sign(one, D2a) * min(C*abs(D2aL), C*abs(D2aR), C*abs(D2aC), abs(D2a)) &&\ + else &&\ + D2aLim = 0 &&\ + end if &&\ + if (abs(D2a) .le. 1.d-12*max(abs(a(i-2)), abs(a(i-1)), abs(a(i)), abs(a(i+1)), abs(a(i+2)))) then &&\ + rhi = 0 &&\ + else &&\ + rhi = D2aLim / D2a &&\ + endif &&\ + if (.not. (rhi .ge. 1.0d0 - 1.d-12)) then &&\ + D3aMin = min(D3aLL, D3aL, D3a, D3aR) &&\ + D3aMax = max(D3aLL, D3aL, D3a, D3aR) &&\ + if (C3 * max(abs(D3aMin), abs(D3aMax)) .le. D3aMax-D3aMin) then &&\ + if (daplus*daminus .lt. 0) then &&\ + aplus(i) = a(i) + daplus * rhi &&\ + aminus(i) = a(i) - daminus * rhi &&\ + else if (abs(daminus) .ge. 2.0d0*abs(daplus)) then &&\ + aminus(i) = a(i) - 2.0d0*(1.0d0-rhi)*daplus - rhi*daminus &&\ + else if (abs(daplus) .ge. 2.0d0*abs(daminus)) then &&\ + aplus(i) = a(i) + 2.0d0*(1.0d0-rhi)*daminus + rhi*daplus &&\ + endif &&\ + endif &&\ + endif &&\ else &&\ - if (abs(daplus) .ge. 2.0d0*abs(daminus)) then &&\ - aplus(i) = a(i) + 2.0d0*daminus &&\ - else if (abs(daminus) .ge. 2.0d0*abs(daplus)) then &&\ - aminus(i) = a(i) - 2.0d0*daplus &&\ - end if &&\ + if (abs(daplus) .ge. 2.0d0*abs(daminus)) then &&\ + aplus(i) = a(i) + 2.0d0*daminus &&\ + else if (abs(daminus) .ge. 2.0d0*abs(daplus)) then &&\ + aminus(i) = a(i) - 2.0d0*daplus &&\ + end if &&\ endif !! Monotonicity of PPM 2011 (McCorquodale & Colella 2011), Sec. 2.4.1, Eq. 23-34. @@ -1290,65 +1290,65 @@ subroutine SimplePPM_temperature_1d(& daplus = aplus(i)-a(i) &&\ daminus = a(i)-aminus(i) &&\ if (daplus*daminus .le. 0 .or. (a(i-2)-a(i))*(a(i)-a(i+2)) .le. 0) then &&\ - D2a = - (12.0d0*a(i) - 6.0d0*(aminus(i)+aplus(i))) && \ - D2aC = a(i-1) - 2.0d0*a(i) + a(i+1) &&\ - D2aL = a(i-2) - 2.0d0*a(i-1) + a(i) &&\ - D2aLL = a(i-3) - 2.0d0*a(i-2) + a(i-1) &&\ - D2aR = a(i) - 2.0d0*a(i+1) + a(i+2) &&\ - D2aRR = a(i+1) - 2.0d0*a(i+2) + a(i+3) &&\ - D3a = D2aR - D2aC &&\ - D3aL = D2aC - D2aL &&\ - D3aR = D2aRR - D2aR &&\ - D3aLL = D2aL - D2aLL &&\ - if (sign(one, D2a) .eq. sign(one, D2aC) .and. sign(one, D2a) .eq. sign(one, D2aL) .and. sign(one, D2a) .eq. sign(one, D2aR)) then &&\ - D2aLim = sign(one, D2a) * min(C*abs(D2aL), C*abs(D2aR), C*abs(D2aC), abs(D2a)) &&\ - else &&\ - D2aLim = 0 &&\ - end if &&\ - if (abs(D2a) .le. 1.d-12*max(abs(a(i-2)), abs(a(i-1)), abs(a(i)), abs(a(i+1)), abs(a(i+2)))) then &&\ - rhi = 0 &&\ - else &&\ - rhi = D2aLim / D2a &&\ - endif &&\ - if (.not. (rhi .ge. 1.0d0 - 1.d-12)) then &&\ - D3aMin = min(D3aLL, D3aL, D3a, D3aR) &&\ - D3aMax = max(D3aLL, D3aL, D3a, D3aR) &&\ - if (C3 * max(abs(D3aMin), abs(D3aMax)) .le. D3aMax-D3aMin) then &&\ - if (abs(daplus) .le. abs(a(i)) .and. abs(daminus) .le. abs(a(i))) then &&\ - if (daplus*daminus .lt. 0) then &&\ - aplus(i) = a(i) + daplus * rhi &&\ - aminus(i) = a(i) - daminus * rhi &&\ - else if (abs(daminus) .ge. 2.0d0*abs(daplus)) then &&\ - aminus(i) = a(i) - 2.0d0*(1.0d0-rhi)*daplus - rhi*daminus &&\ - else if (abs(daplus) .ge. 2.0d0*abs(daminus)) then &&\ - aplus(i) = a(i) + 2.0d0*(1.0d0-rhi)*daminus + rhi*daplus &&\ - endif &&\ - else &&\ - if (daplus*daminus .lt. 0) then &&\ - aplus(i) = a(i) &&\ - aminus(i) = a(i) &&\ - else if (abs(daminus) .ge. 2.0d0*abs(daplus)) then &&\ - aminus(i) = a(i) &&\ - else if (abs(daplus) .ge. 2.0d0*abs(daminus)) then &&\ - aplus(i) = a(i) &&\ - endif &&\ - endif &&\ - endif &&\ - endif &&\ + D2a = - (12.0d0*a(i) - 6.0d0*(aminus(i)+aplus(i))) && \ + D2aC = a(i-1) - 2.0d0*a(i) + a(i+1) &&\ + D2aL = a(i-2) - 2.0d0*a(i-1) + a(i) &&\ + D2aLL = a(i-3) - 2.0d0*a(i-2) + a(i-1) &&\ + D2aR = a(i) - 2.0d0*a(i+1) + a(i+2) &&\ + D2aRR = a(i+1) - 2.0d0*a(i+2) + a(i+3) &&\ + D3a = D2aR - D2aC &&\ + D3aL = D2aC - D2aL &&\ + D3aR = D2aRR - D2aR &&\ + D3aLL = D2aL - D2aLL &&\ + if (sign(one, D2a) .eq. sign(one, D2aC) .and. sign(one, D2a) .eq. sign(one, D2aL) .and. sign(one, D2a) .eq. sign(one, D2aR)) then &&\ + D2aLim = sign(one, D2a) * min(C*abs(D2aL), C*abs(D2aR), C*abs(D2aC), abs(D2a)) &&\ + else &&\ + D2aLim = 0 &&\ + end if &&\ + if (abs(D2a) .le. 1.d-12*max(abs(a(i-2)), abs(a(i-1)), abs(a(i)), abs(a(i+1)), abs(a(i+2)))) then &&\ + rhi = 0 &&\ + else &&\ + rhi = D2aLim / D2a &&\ + endif &&\ + if (.not. (rhi .ge. 1.0d0 - 1.d-12)) then &&\ + D3aMin = min(D3aLL, D3aL, D3a, D3aR) &&\ + D3aMax = max(D3aLL, D3aL, D3a, D3aR) &&\ + if (C3 * max(abs(D3aMin), abs(D3aMax)) .le. D3aMax-D3aMin) then &&\ + if (abs(daplus) .le. abs(a(i)) .and. abs(daminus) .le. abs(a(i))) then &&\ + if (daplus*daminus .lt. 0) then &&\ + aplus(i) = a(i) + daplus * rhi &&\ + aminus(i) = a(i) - daminus * rhi &&\ + else if (abs(daminus) .ge. 2.0d0*abs(daplus)) then &&\ + aminus(i) = a(i) - 2.0d0*(1.0d0-rhi)*daplus - rhi*daminus &&\ + else if (abs(daplus) .ge. 2.0d0*abs(daminus)) then &&\ + aplus(i) = a(i) + 2.0d0*(1.0d0-rhi)*daminus + rhi*daplus &&\ + endif &&\ + else &&\ + if (daplus*daminus .lt. 0) then &&\ + aplus(i) = a(i) &&\ + aminus(i) = a(i) &&\ + else if (abs(daminus) .ge. 2.0d0*abs(daplus)) then &&\ + aminus(i) = a(i) &&\ + else if (abs(daplus) .ge. 2.0d0*abs(daminus)) then &&\ + aplus(i) = a(i) &&\ + endif &&\ + endif &&\ + endif &&\ + endif &&\ else &&\ - if (abs(daplus) .le. abs(a(i)) .and. abs(daminus) .le. abs(a(i))) then &&\ - if (abs(daplus) .ge. 2.0d0*abs(daminus)) then &&\ - aplus(i) = a(i) + 2.0d0*daminus &&\ - else if (abs(daminus) .ge. 2.0d0*abs(daplus)) then &&\ - aminus(i) = a(i) - 2.0d0*daplus &&\ - end if &&\ - else &&\ - if (abs(daplus) .ge. 2.0d0*abs(daminus)) then &&\ - aplus(i) = a(i) &&\ - else if (abs(daminus) .ge. 2.0d0*abs(daplus)) then &&\ - aminus(i) = a(i) &&\ - end if &&\ - endif &&\ + if (abs(daplus) .le. abs(a(i)) .and. abs(daminus) .le. abs(a(i))) then &&\ + if (abs(daplus) .ge. 2.0d0*abs(daminus)) then &&\ + aplus(i) = a(i) + 2.0d0*daminus &&\ + else if (abs(daminus) .ge. 2.0d0*abs(daplus)) then &&\ + aminus(i) = a(i) - 2.0d0*daplus &&\ + end if &&\ + else &&\ + if (abs(daplus) .ge. 2.0d0*abs(daminus)) then &&\ + aplus(i) = a(i) &&\ + else if (abs(daminus) .ge. 2.0d0*abs(daplus)) then &&\ + aminus(i) = a(i) &&\ + end if &&\ + endif &&\ endif @@ -1358,36 +1358,36 @@ subroutine SimplePPM_temperature_1d(& daplus = aplus(i)-a(i) &&\ daminus = a(i)-aminus(i) &&\ if (daplus*daminus .le. 0 .or. (a(i-2)-a(i))*(a(i)-a(i+2)) .le. 0) then &&\ - D2a = - (12.0d0*a(i) - 6.0d0*(aminus(i)+aplus(i))) && \ - D2aC = a(i-1) - 2.0d0*a(i) + a(i+1) &&\ - D2aL = a(i-2) - 2.0d0*a(i-1) + a(i) &&\ - D2aR = a(i) - 2.0d0*a(i+1) + a(i+2) &&\ - if (sign(one, D2a) .eq. sign(one, D2aC) .and. sign(one, D2a) .eq. sign(one, D2aL) .and. sign(one, D2a) .eq. sign(one, D2aR)) then &&\ - D2aLim = sign(one, D2a) * min(C*abs(D2aL), C*abs(D2aR), C*abs(D2aC), abs(D2a)) &&\ - else &&\ - D2aLim = 0 &&\ - end if &&\ - if (abs(D2a) .le. 1.d-12*max(abs(a(i-2)), abs(a(i-1)), abs(a(i)), abs(a(i+1)), abs(a(i+2)))) then &&\ - rhi = 0 &&\ - else &&\ - rhi = D2aLim / D2a &&\ - endif &&\ - if (.not. (rhi .ge. 1.0d0 - 1.d-12)) then &&\ - if (daplus*daminus .lt. 0) then &&\ - aplus(i) = a(i) + daplus * rhi &&\ - aminus(i) = a(i) - daminus * rhi &&\ - else if (abs(daminus) .ge. 2.0d0*abs(daplus)) then &&\ - aminus(i) = a(i) - 2.0d0*(1.0d0-rhi)*daplus - rhi*daminus &&\ - else if (abs(daplus) .ge. 2.0d0*abs(daminus)) then &&\ - aplus(i) = a(i) + 2.0d0*(1.0d0-rhi)*daminus + rhi*daplus &&\ - endif &&\ - endif &&\ + D2a = - (12.0d0*a(i) - 6.0d0*(aminus(i)+aplus(i))) && \ + D2aC = a(i-1) - 2.0d0*a(i) + a(i+1) &&\ + D2aL = a(i-2) - 2.0d0*a(i-1) + a(i) &&\ + D2aR = a(i) - 2.0d0*a(i+1) + a(i+2) &&\ + if (sign(one, D2a) .eq. sign(one, D2aC) .and. sign(one, D2a) .eq. sign(one, D2aL) .and. sign(one, D2a) .eq. sign(one, D2aR)) then &&\ + D2aLim = sign(one, D2a) * min(C*abs(D2aL), C*abs(D2aR), C*abs(D2aC), abs(D2a)) &&\ + else &&\ + D2aLim = 0 &&\ + end if &&\ + if (abs(D2a) .le. 1.d-12*max(abs(a(i-2)), abs(a(i-1)), abs(a(i)), abs(a(i+1)), abs(a(i+2)))) then &&\ + rhi = 0 &&\ + else &&\ + rhi = D2aLim / D2a &&\ + endif &&\ + if (.not. (rhi .ge. 1.0d0 - 1.d-12)) then &&\ + if (daplus*daminus .lt. 0) then &&\ + aplus(i) = a(i) + daplus * rhi &&\ + aminus(i) = a(i) - daminus * rhi &&\ + else if (abs(daminus) .ge. 2.0d0*abs(daplus)) then &&\ + aminus(i) = a(i) - 2.0d0*(1.0d0-rhi)*daplus - rhi*daminus &&\ + else if (abs(daplus) .ge. 2.0d0*abs(daminus)) then &&\ + aplus(i) = a(i) + 2.0d0*(1.0d0-rhi)*daminus + rhi*daplus &&\ + endif &&\ + endif &&\ else &&\ - if (abs(daplus) .ge. 2.0d0*abs(daminus)) then &&\ - aplus(i) = a(i) + 2.0d0*daminus &&\ - else if (abs(daminus) .ge. 2.0d0*abs(daplus)) then &&\ - aminus(i) = a(i) - 2.0d0*daplus &&\ - end if &&\ + if (abs(daplus) .ge. 2.0d0*abs(daminus)) then &&\ + aplus(i) = a(i) + 2.0d0*daminus &&\ + else if (abs(daminus) .ge. 2.0d0*abs(daplus)) then &&\ + aminus(i) = a(i) - 2.0d0*daplus &&\ + end if &&\ endif @@ -1399,84 +1399,84 @@ subroutine SimplePPM_temperature_1d(& daplus = aplus(i)-a(i) &&\ daminus = a(i)-aminus(i) &&\ if (daplus*daminus .le. 0 .or. (a(i-2)-a(i))*(a(i)-a(i+2)) .le. 0) then &&\ - D2a = - (12.0d0*a(i) - 6.0d0*(aminus(i)+aplus(i))) && \ - D2aC = a(i-1) - 2.0d0*a(i) + a(i+1) &&\ - D2aL = a(i-2) - 2.0d0*a(i-1) + a(i) &&\ - D2aR = a(i) - 2.0d0*a(i+1) + a(i+2) &&\ - if (sign(one, D2a) .eq. sign(one, D2aC) .and. sign(one, D2a) .eq. sign(one, D2aL) .and. sign(one, D2a) .eq. sign(one, D2aR)) then &&\ - D2aLim = sign(one, D2a) * min(C*abs(D2aL), C*abs(D2aR), C*abs(D2aC), abs(D2a)) &&\ - else &&\ - D2aLim = 0 &&\ - end if &&\ - if (abs(D2a) .le. 1.d-12*max(abs(a(i-2)), abs(a(i-1)), abs(a(i)), abs(a(i+1)), abs(a(i+2)))) then &&\ - rhi = 0 &&\ - else &&\ - rhi = D2aLim / D2a &&\ - endif &&\ - if (.not. (rhi .ge. 1.0d0 - 1.d-12)) then &&\ - if (abs(daplus) .le. abs(a(i)) .and. abs(daminus) .le. abs(a(i))) then &&\ - if (daplus*daminus .lt. 0) then &&\ - aplus(i) = a(i) + daplus * rhi &&\ - aminus(i) = a(i) - daminus * rhi &&\ - else if (abs(daminus) .ge. 2.0d0*abs(daplus)) then &&\ - aminus(i) = a(i) - 2.0d0*(1.0d0-rhi)*daplus - rhi*daminus &&\ - else if (abs(daplus) .ge. 2.0d0*abs(daminus)) then &&\ - aplus(i) = a(i) + 2.0d0*(1.0d0-rhi)*daminus + rhi*daplus &&\ - endif &&\ - else &&\ - if (daplus*daminus .lt. 0) then &&\ - aplus(i) = a(i) &&\ - aminus(i) = a(i) &&\ - else if (abs(daminus) .ge. 2.0d0*abs(daplus)) then &&\ - aminus(i) = a(i) &&\ - else if (abs(daplus) .ge. 2.0d0*abs(daminus)) then &&\ - aplus(i) = a(i) &&\ - endif &&\ - endif &&\ - endif &&\ + D2a = - (12.0d0*a(i) - 6.0d0*(aminus(i)+aplus(i))) && \ + D2aC = a(i-1) - 2.0d0*a(i) + a(i+1) &&\ + D2aL = a(i-2) - 2.0d0*a(i-1) + a(i) &&\ + D2aR = a(i) - 2.0d0*a(i+1) + a(i+2) &&\ + if (sign(one, D2a) .eq. sign(one, D2aC) .and. sign(one, D2a) .eq. sign(one, D2aL) .and. sign(one, D2a) .eq. sign(one, D2aR)) then &&\ + D2aLim = sign(one, D2a) * min(C*abs(D2aL), C*abs(D2aR), C*abs(D2aC), abs(D2a)) &&\ + else &&\ + D2aLim = 0 &&\ + end if &&\ + if (abs(D2a) .le. 1.d-12*max(abs(a(i-2)), abs(a(i-1)), abs(a(i)), abs(a(i+1)), abs(a(i+2)))) then &&\ + rhi = 0 &&\ + else &&\ + rhi = D2aLim / D2a &&\ + endif &&\ + if (.not. (rhi .ge. 1.0d0 - 1.d-12)) then &&\ + if (abs(daplus) .le. abs(a(i)) .and. abs(daminus) .le. abs(a(i))) then &&\ + if (daplus*daminus .lt. 0) then &&\ + aplus(i) = a(i) + daplus * rhi &&\ + aminus(i) = a(i) - daminus * rhi &&\ + else if (abs(daminus) .ge. 2.0d0*abs(daplus)) then &&\ + aminus(i) = a(i) - 2.0d0*(1.0d0-rhi)*daplus - rhi*daminus &&\ + else if (abs(daplus) .ge. 2.0d0*abs(daminus)) then &&\ + aplus(i) = a(i) + 2.0d0*(1.0d0-rhi)*daminus + rhi*daplus &&\ + endif &&\ + else &&\ + if (daplus*daminus .lt. 0) then &&\ + aplus(i) = a(i) &&\ + aminus(i) = a(i) &&\ + else if (abs(daminus) .ge. 2.0d0*abs(daplus)) then &&\ + aminus(i) = a(i) &&\ + else if (abs(daplus) .ge. 2.0d0*abs(daminus)) then &&\ + aplus(i) = a(i) &&\ + endif &&\ + endif &&\ + endif &&\ else &&\ - if (abs(daplus) .le. abs(a(i)) .and. abs(daminus) .le. abs(a(i))) then &&\ - if (abs(daplus) .ge. 2.0d0*abs(daminus)) then &&\ - aplus(i) = a(i) + 2.0d0*daminus &&\ - else if (abs(daminus) .ge. 2.0d0*abs(daplus)) then &&\ - aminus(i) = a(i) - 2.0d0*daplus &&\ - end if &&\ - else &&\ - if (abs(daplus) .ge. 2.0d0*abs(daminus)) then &&\ - aplus(i) = a(i) &&\ - else if (abs(daminus) .ge. 2.0d0*abs(daplus)) then &&\ - aminus(i) = a(i) &&\ - end if &&\ - endif &&\ + if (abs(daplus) .le. abs(a(i)) .and. abs(daminus) .le. abs(a(i))) then &&\ + if (abs(daplus) .ge. 2.0d0*abs(daminus)) then &&\ + aplus(i) = a(i) + 2.0d0*daminus &&\ + else if (abs(daminus) .ge. 2.0d0*abs(daplus)) then &&\ + aminus(i) = a(i) - 2.0d0*daplus &&\ + end if &&\ + else &&\ + if (abs(daplus) .ge. 2.0d0*abs(daminus)) then &&\ + aplus(i) = a(i) &&\ + else if (abs(daminus) .ge. 2.0d0*abs(daplus)) then &&\ + aminus(i) = a(i) &&\ + end if &&\ + endif &&\ endif if (use_enhanced_ppm .eq. 1) then !! Constrain parabolic profiles, PPM 2011/2008 if (PPM3) then - do i = 3, nx - 2 + do i = 3, nx - 2 MONT_WITH_LOCAL_EXTREMUM_EPS(tempminus,temperature,tempplus, enhanced_ppm_C2) enddo else - do i = 4, nx - 3 + do i = 4, nx - 3 MONT_WITH_LOCAL_EXTREMUM_STENCIL4_EPS(tempminus,temperature,tempplus, enhanced_ppm_C2, enhanced_ppm_C3) enddo end if ! Apply flattening after constraining the parabolic profiles if (PPM3) then !!$ Implement C&W, page 197, but with a workaround which allows to use stencil=3. - do i = 3, nx - 2 - flatten = tilde_flatten(i) + do i = 3, nx - 2 + flatten = tilde_flatten(i) tempplus(i) = flatten * tempplus(i) + (1.d0 - flatten) * temperature(i) tempminus(i) = flatten * tempminus(i) + (1.d0 - flatten) * temperature(i) - end do + end do else !!$ Really implement C&W, page 197; which requires stencil 4. - do i = 4, nx - 3 - s=sign(one, -dpress(i)) - flatten = max(tilde_flatten(i), tilde_flatten(i+s)) + do i = 4, nx - 3 + s=sign(one, -dpress(i)) + flatten = max(tilde_flatten(i), tilde_flatten(i+s)) tempplus(i) = flatten * tempplus(i) + (1.d0 - flatten) * temperature(i) tempminus(i) = flatten * tempminus(i) + (1.d0 - flatten) * temperature(i) - end do + end do end if endif @@ -1484,7 +1484,7 @@ subroutine SimplePPM_temperature_1d(& if (use_enhanced_ppm .eq. 0) then !! Constrain parabolic profiles, PPM 1984 do i = GRHydro_stencil, nx - GRHydro_stencil + 1 - ! original Colella&Woodward monotonicity preservation + ! original Colella&Woodward monotonicity preservation MONT(tempminus,temperature,tempplus) enddo endif @@ -1850,38 +1850,38 @@ subroutine SimplePPM_ye_1d(nx,dx,rho,velx,vely,velz, & #undef LIMIT #define LIMIT(a,ah,C, alim) \ if ((min(a(i),a(i+1)) .le. ah) .and. (ah .le. max(a(i),a(i+1)))) then &&\ - alim = ah &&\ + alim = ah &&\ else &&\ - D2a = 3.0d0 * (a(i) - 2.0d0*ah + a(i+1)) &&\ - D2aL = (a(i-1) - 2.0d0*a(i) + a(i+1)) &&\ - D2aR = (a(i) - 2.0d0*a(i+1) + a(i+2)) &&\ - if (D2a*D2aR .ge. 0 .and. D2a*D2aL .ge. 0) then &&\ - alim = 0.5d0*(a(i)+a(i+1)) - sign(one, D2a)*min(C*abs(D2aL), C*abs(D2aR), abs(D2a))/3.0d0 &&\ - else &&\ - alim = 0.5d0*(a(i)+a(i+1)) &&\ - end if &&\ + D2a = 3.0d0 * (a(i) - 2.0d0*ah + a(i+1)) &&\ + D2aL = (a(i-1) - 2.0d0*a(i) + a(i+1)) &&\ + D2aR = (a(i) - 2.0d0*a(i+1) + a(i+2)) &&\ + if (D2a*D2aR .ge. 0 .and. D2a*D2aL .ge. 0) then &&\ + alim = 0.5d0*(a(i)+a(i+1)) - sign(one, D2a)*min(C*abs(D2aL), C*abs(D2aR), abs(D2a))/3.0d0 &&\ + else &&\ + alim = 0.5d0*(a(i)+a(i+1)) &&\ + end if &&\ endif if (PPM3) then - !! We initialize "plus" \equiv a_j+1/2 with (16) via APPROX_AT_CELL_INTERFACE, - !! then checking for (13) of Colella & Sekora 2008 and applying - !! (18) and (19) if (13) is not satisfied. This is done with LIMIT. - !! There, we also switch to lower order if we are near atmosphere values. - do i = 2, nx-2 - APPROX_AT_CELL_INTERFACE(Y_e, Y_e_plus(i)) - LIMIT(Y_e, Y_e_plus(i), enhanced_ppm_C2, Y_e_plus(i)) - Y_e_minus(i+1) = Y_e_plus(i) - end do - + !! We initialize "plus" \equiv a_j+1/2 with (16) via APPROX_AT_CELL_INTERFACE, + !! then checking for (13) of Colella & Sekora 2008 and applying + !! (18) and (19) if (13) is not satisfied. This is done with LIMIT. + !! There, we also switch to lower order if we are near atmosphere values. + do i = 2, nx-2 + APPROX_AT_CELL_INTERFACE(Y_e, Y_e_plus(i)) + LIMIT(Y_e, Y_e_plus(i), enhanced_ppm_C2, Y_e_plus(i)) + Y_e_minus(i+1) = Y_e_plus(i) + end do + else - !! Same as above but for 4 stencil points using (17) of Colella & Sekora 2008 as - !! initial states. - do i = 3, nx-3 - APPROX_AT_CELL_INTERFACE_STENCIL4(Y_e, Y_e_plus(i)) - LIMIT(Y_e, Y_e_plus(i), enhanced_ppm_C2, Y_e_plus(i)) - Y_e_minus(i+1) = Y_e_plus(i) - end do + !! Same as above but for 4 stencil points using (17) of Colella & Sekora 2008 as + !! initial states. + do i = 3, nx-3 + APPROX_AT_CELL_INTERFACE_STENCIL4(Y_e, Y_e_plus(i)) + LIMIT(Y_e, Y_e_plus(i), enhanced_ppm_C2, Y_e_plus(i)) + Y_e_minus(i+1) = Y_e_plus(i) + end do endif end if @@ -1917,22 +1917,22 @@ subroutine SimplePPM_ye_1d(nx,dx,rho,velx,vely,velz, & if (use_enhanced_ppm .eq. 0) then if (PPM3) then - do i = 3, nx - 2 - flatten = tilde_flatten(i) - Y_e_plus(i) = flatten * Y_e_plus(i) & - + (1.d0 - flatten) * Y_e(i) - Y_e_minus(i) = flatten * Y_e_minus(i) & - + (1.d0 - flatten) * Y_e(i) - end do + do i = 3, nx - 2 + flatten = tilde_flatten(i) + Y_e_plus(i) = flatten * Y_e_plus(i) & + + (1.d0 - flatten) * Y_e(i) + Y_e_minus(i) = flatten * Y_e_minus(i) & + + (1.d0 - flatten) * Y_e(i) + end do else !!$ Really implement C&W, page 197; which requires stencil 4. - do i = 4, nx - 3 - s=sign(one, -dpress(i)) - flatten = max(tilde_flatten(i), tilde_flatten(i+s)) - Y_e_plus(i) = flatten * Y_e_plus(i) + & - (1.d0 - flatten) * Y_e(i) - Y_e_minus(i) = flatten * Y_e_minus(i) & - + (1.d0 - flatten) * Y_e(i) - end do + do i = 4, nx - 3 + s=sign(one, -dpress(i)) + flatten = max(tilde_flatten(i), tilde_flatten(i+s)) + Y_e_plus(i) = flatten * Y_e_plus(i) + & + (1.d0 - flatten) * Y_e(i) + Y_e_minus(i) = flatten * Y_e_minus(i) & + + (1.d0 - flatten) * Y_e(i) + end do end if endif @@ -1944,46 +1944,46 @@ subroutine SimplePPM_ye_1d(nx,dx,rho,velx,vely,velz, & daplus = aplus(i)-a(i) &&\ daminus = a(i)-aminus(i) &&\ if (daplus*daminus .le. 0 .or. (a(i-2)-a(i))*(a(i)-a(i+2)) .le. 0) then &&\ - D2a = - (12.0d0*a(i) - 6.0d0*(aminus(i)+aplus(i))) && \ - D2aC = a(i-1) - 2.0d0*a(i) + a(i+1) &&\ - D2aL = a(i-2) - 2.0d0*a(i-1) + a(i) &&\ - D2aLL = a(i-3) - 2.0d0*a(i-2) + a(i-1) &&\ - D2aR = a(i) - 2.0d0*a(i+1) + a(i+2) &&\ - D2aRR = a(i+1) - 2.0d0*a(i+2) + a(i+3) &&\ - D3a = D2aR - D2aC &&\ - D3aL = D2aC - D2aL &&\ - D3aR = D2aRR - D2aR &&\ - D3aLL = D2aL - D2aLL &&\ - if (sign(one, D2a) .eq. sign(one, D2aC) .and. sign(one, D2a) .eq. sign(one, D2aL) .and. sign(one, D2a) .eq. sign(one, D2aR)) then &&\ - D2aLim = sign(one, D2a) * min(C*abs(D2aL), C*abs(D2aR), C*abs(D2aC), abs(D2a)) &&\ - else &&\ - D2aLim = 0 &&\ - end if &&\ - if (abs(D2a) .le. 1.d-12*max(abs(a(i-2)), abs(a(i-1)), abs(a(i)), abs(a(i+1)), abs(a(i+2)))) then &&\ - rhi = 0 &&\ - else &&\ - rhi = D2aLim / D2a &&\ - endif &&\ - if (.not. (rhi .ge. 1.0d0 - 1.d-12)) then &&\ - D3aMin = min(D3aLL, D3aL, D3a, D3aR) &&\ - D3aMax = max(D3aLL, D3aL, D3a, D3aR) &&\ - if (C3 * max(abs(D3aMin), abs(D3aMax)) .le. D3aMax-D3aMin) then &&\ - if (daplus*daminus .lt. 0) then &&\ - aplus(i) = a(i) + daplus * rhi &&\ - aminus(i) = a(i) - daminus * rhi &&\ - else if (abs(daminus) .ge. 2.0d0*abs(daplus)) then &&\ - aminus(i) = a(i) - 2.0d0*(1.0d0-rhi)*daplus - rhi*daminus &&\ - else if (abs(daplus) .ge. 2.0d0*abs(daminus)) then &&\ - aplus(i) = a(i) + 2.0d0*(1.0d0-rhi)*daminus + rhi*daplus &&\ - endif &&\ - endif &&\ - endif &&\ + D2a = - (12.0d0*a(i) - 6.0d0*(aminus(i)+aplus(i))) && \ + D2aC = a(i-1) - 2.0d0*a(i) + a(i+1) &&\ + D2aL = a(i-2) - 2.0d0*a(i-1) + a(i) &&\ + D2aLL = a(i-3) - 2.0d0*a(i-2) + a(i-1) &&\ + D2aR = a(i) - 2.0d0*a(i+1) + a(i+2) &&\ + D2aRR = a(i+1) - 2.0d0*a(i+2) + a(i+3) &&\ + D3a = D2aR - D2aC &&\ + D3aL = D2aC - D2aL &&\ + D3aR = D2aRR - D2aR &&\ + D3aLL = D2aL - D2aLL &&\ + if (sign(one, D2a) .eq. sign(one, D2aC) .and. sign(one, D2a) .eq. sign(one, D2aL) .and. sign(one, D2a) .eq. sign(one, D2aR)) then &&\ + D2aLim = sign(one, D2a) * min(C*abs(D2aL), C*abs(D2aR), C*abs(D2aC), abs(D2a)) &&\ + else &&\ + D2aLim = 0 &&\ + end if &&\ + if (abs(D2a) .le. 1.d-12*max(abs(a(i-2)), abs(a(i-1)), abs(a(i)), abs(a(i+1)), abs(a(i+2)))) then &&\ + rhi = 0 &&\ + else &&\ + rhi = D2aLim / D2a &&\ + endif &&\ + if (.not. (rhi .ge. 1.0d0 - 1.d-12)) then &&\ + D3aMin = min(D3aLL, D3aL, D3a, D3aR) &&\ + D3aMax = max(D3aLL, D3aL, D3a, D3aR) &&\ + if (C3 * max(abs(D3aMin), abs(D3aMax)) .le. D3aMax-D3aMin) then &&\ + if (daplus*daminus .lt. 0) then &&\ + aplus(i) = a(i) + daplus * rhi &&\ + aminus(i) = a(i) - daminus * rhi &&\ + else if (abs(daminus) .ge. 2.0d0*abs(daplus)) then &&\ + aminus(i) = a(i) - 2.0d0*(1.0d0-rhi)*daplus - rhi*daminus &&\ + else if (abs(daplus) .ge. 2.0d0*abs(daminus)) then &&\ + aplus(i) = a(i) + 2.0d0*(1.0d0-rhi)*daminus + rhi*daplus &&\ + endif &&\ + endif &&\ + endif &&\ else &&\ - if (abs(daplus) .ge. 2.0d0*abs(daminus)) then &&\ - aplus(i) = a(i) + 2.0d0*daminus &&\ - else if (abs(daminus) .ge. 2.0d0*abs(daplus)) then &&\ - aminus(i) = a(i) - 2.0d0*daplus &&\ - end if &&\ + if (abs(daplus) .ge. 2.0d0*abs(daminus)) then &&\ + aplus(i) = a(i) + 2.0d0*daminus &&\ + else if (abs(daminus) .ge. 2.0d0*abs(daplus)) then &&\ + aminus(i) = a(i) - 2.0d0*daplus &&\ + end if &&\ endif @@ -1997,36 +1997,36 @@ subroutine SimplePPM_ye_1d(nx,dx,rho,velx,vely,velz, & daplus = aplus(i)-a(i) &&\ daminus = a(i)-aminus(i) &&\ if (daplus*daminus .le. 0 .or. (a(i-2)-a(i))*(a(i)-a(i+2)) .le. 0) then &&\ - D2a = - (12.0d0*a(i) - 6.0d0*(aminus(i)+aplus(i))) && \ - D2aC = a(i-1) - 2.0d0*a(i) + a(i+1) &&\ - D2aL = a(i-2) - 2.0d0*a(i-1) + a(i) &&\ - D2aR = a(i) - 2.0d0*a(i+1) + a(i+2) &&\ - if (sign(one, D2a) .eq. sign(one, D2aC) .and. sign(one, D2a) .eq. sign(one, D2aL) .and. sign(one, D2a) .eq. sign(one, D2aR)) then &&\ - D2aLim = sign(one, D2a) * min(C*abs(D2aL), C*abs(D2aR), C*abs(D2aC), abs(D2a)) &&\ - else &&\ - D2aLim = 0 &&\ - end if &&\ - if (abs(D2a) .le. 1.d-12*max(abs(a(i-2)), abs(a(i-1)), abs(a(i)), abs(a(i+1)), abs(a(i+2)))) then &&\ - rhi = 0 &&\ - else &&\ - rhi = D2aLim / D2a &&\ - endif &&\ - if (.not. (rhi .ge. 1.0d0 - 1.d-12)) then &&\ - if (daplus*daminus .lt. 0) then &&\ - aplus(i) = a(i) + daplus * rhi &&\ - aminus(i) = a(i) - daminus * rhi &&\ - else if (abs(daminus) .ge. 2.0d0*abs(daplus)) then &&\ - aminus(i) = a(i) - 2.0d0*(1.0d0-rhi)*daplus - rhi*daminus &&\ - else if (abs(daplus) .ge. 2.0d0*abs(daminus)) then &&\ - aplus(i) = a(i) + 2.0d0*(1.0d0-rhi)*daminus + rhi*daplus &&\ - endif &&\ - endif &&\ + D2a = - (12.0d0*a(i) - 6.0d0*(aminus(i)+aplus(i))) && \ + D2aC = a(i-1) - 2.0d0*a(i) + a(i+1) &&\ + D2aL = a(i-2) - 2.0d0*a(i-1) + a(i) &&\ + D2aR = a(i) - 2.0d0*a(i+1) + a(i+2) &&\ + if (sign(one, D2a) .eq. sign(one, D2aC) .and. sign(one, D2a) .eq. sign(one, D2aL) .and. sign(one, D2a) .eq. sign(one, D2aR)) then &&\ + D2aLim = sign(one, D2a) * min(C*abs(D2aL), C*abs(D2aR), C*abs(D2aC), abs(D2a)) &&\ + else &&\ + D2aLim = 0 &&\ + end if &&\ + if (abs(D2a) .le. 1.d-12*max(abs(a(i-2)), abs(a(i-1)), abs(a(i)), abs(a(i+1)), abs(a(i+2)))) then &&\ + rhi = 0 &&\ + else &&\ + rhi = D2aLim / D2a &&\ + endif &&\ + if (.not. (rhi .ge. 1.0d0 - 1.d-12)) then &&\ + if (daplus*daminus .lt. 0) then &&\ + aplus(i) = a(i) + daplus * rhi &&\ + aminus(i) = a(i) - daminus * rhi &&\ + else if (abs(daminus) .ge. 2.0d0*abs(daplus)) then &&\ + aminus(i) = a(i) - 2.0d0*(1.0d0-rhi)*daplus - rhi*daminus &&\ + else if (abs(daplus) .ge. 2.0d0*abs(daminus)) then &&\ + aplus(i) = a(i) + 2.0d0*(1.0d0-rhi)*daminus + rhi*daplus &&\ + endif &&\ + endif &&\ else &&\ - if (abs(daplus) .ge. 2.0d0*abs(daminus)) then &&\ - aplus(i) = a(i) - 2.0d0*(aminus(i) - a(i)) &&\ - else if (abs(daminus) .ge. 2.0d0*abs(daplus)) then &&\ - aminus(i) = a(i) - 2.0d0*(aplus(i) - a(i)) &&\ - end if &&\ + if (abs(daplus) .ge. 2.0d0*abs(daminus)) then &&\ + aplus(i) = a(i) - 2.0d0*(aminus(i) - a(i)) &&\ + else if (abs(daminus) .ge. 2.0d0*abs(daplus)) then &&\ + aminus(i) = a(i) - 2.0d0*(aplus(i) - a(i)) &&\ + end if &&\ endif @@ -2036,33 +2036,33 @@ subroutine SimplePPM_ye_1d(nx,dx,rho,velx,vely,velz, & if (use_enhanced_ppm .eq. 1) then ! Constrain parabolic profiles if (PPM3) then - do i = GRHydro_stencil, nx - GRHydro_stencil + 1 - MON_WITH_LOCAL_EXTREMUM(Y_e_minus,Y_e,Y_e_plus, enhanced_ppm_C2) - end do + do i = GRHydro_stencil, nx - GRHydro_stencil + 1 + MON_WITH_LOCAL_EXTREMUM(Y_e_minus,Y_e,Y_e_plus, enhanced_ppm_C2) + end do else - do i = GRHydro_stencil, nx - GRHydro_stencil + 1 - MON_WITH_LOCAL_EXTREMUM_STENCIL4(Y_e_minus,Y_e,Y_e_plus, enhanced_ppm_C2, enhanced_ppm_C3) - end do + do i = GRHydro_stencil, nx - GRHydro_stencil + 1 + MON_WITH_LOCAL_EXTREMUM_STENCIL4(Y_e_minus,Y_e,Y_e_plus, enhanced_ppm_C2, enhanced_ppm_C3) + end do endif ! In PPM 2011/2008, flattening is applied after constraining parabolic profiles. if (PPM3) then - do i = 3, nx - 2 - flatten = tilde_flatten(i) - Y_e_plus(i) = flatten * Y_e_plus(i) & - + (1.d0 - flatten) * Y_e(i) - Y_e_minus(i) = flatten * Y_e_minus(i) & - + (1.d0 - flatten) * Y_e(i) - end do + do i = 3, nx - 2 + flatten = tilde_flatten(i) + Y_e_plus(i) = flatten * Y_e_plus(i) & + + (1.d0 - flatten) * Y_e(i) + Y_e_minus(i) = flatten * Y_e_minus(i) & + + (1.d0 - flatten) * Y_e(i) + end do else !!$ Really implement C&W, page 197; which requires stencil 4. - do i = 4, nx - 3 - s=sign(one, -dpress(i)) - flatten = max(tilde_flatten(i), tilde_flatten(i+s)) - Y_e_plus(i) = flatten * Y_e_plus(i) + & - (1.d0 - flatten) * Y_e(i) - Y_e_minus(i) = flatten * Y_e_minus(i) & - + (1.d0 - flatten) * Y_e(i) - end do + do i = 4, nx - 3 + s=sign(one, -dpress(i)) + flatten = max(tilde_flatten(i), tilde_flatten(i+s)) + Y_e_plus(i) = flatten * Y_e_plus(i) + & + (1.d0 - flatten) * Y_e(i) + Y_e_minus(i) = flatten * Y_e_minus(i) & + + (1.d0 - flatten) * Y_e(i) + end do end if else diff --git a/src/GRHydro_Tmunu.F90 b/src/GRHydro_Tmunu.F90 index 6a2eee3..d0b5ea6 100644 --- a/src/GRHydro_Tmunu.F90 +++ b/src/GRHydro_Tmunu.F90 @@ -127,10 +127,10 @@ dampfac = 0.0 continue ! no need to add anything to Tmunu at the current point (it's zero anyway!) endif - else - dampfac = 1.0 - endif - + else + dampfac = 1.0 + endif + !!$ Calculate Tmunu (the lower components!). eTtt(i,j,k) = eTtt(i,j,k) + dampfac * (rhoenthalpy*utlow**2 + press(i,j,k)*(beta2 - alp(i,j,k)**2)) diff --git a/src/GRHydro_TransformTensorBasis.c b/src/GRHydro_TransformTensorBasis.c index 6e86a9a..8d71167 100644 --- a/src/GRHydro_TransformTensorBasis.c +++ b/src/GRHydro_TransformTensorBasis.c @@ -311,12 +311,12 @@ void GRHydroTransformPrimToLocalBasis(CCTK_ARGUMENTS) CCTK_INT idx2 = CCTK_VECTGFINDEX3D(cctkGH, i,j,k,1); CCTK_INT idx3 = CCTK_VECTGFINDEX3D(cctkGH, i,j,k,2); - /* Transform primitive velocity from global to local basis. - * Since velocity has contravariant index, use Jacobian */ - lvel[idx1] = vel[idx1]*J11[idx] + vel[idx2]*J12[idx] + vel[idx3]*J13[idx]; - lvel[idx2] = vel[idx1]*J21[idx] + vel[idx2]*J22[idx] + vel[idx3]*J23[idx]; - lvel[idx3] = vel[idx1]*J31[idx] + vel[idx2]*J32[idx] + vel[idx3]*J33[idx]; - + /* Transform primitive velocity from global to local basis. + * Since velocity has contravariant index, use Jacobian */ + lvel[idx1] = vel[idx1]*J11[idx] + vel[idx2]*J12[idx] + vel[idx3]*J13[idx]; + lvel[idx2] = vel[idx1]*J21[idx] + vel[idx2]*J22[idx] + vel[idx3]*J23[idx]; + lvel[idx3] = vel[idx1]*J31[idx] + vel[idx2]*J32[idx] + vel[idx3]*J33[idx]; + if(*evolve_MHD) { /* Transform primitive B-field from global to local basis. @@ -423,21 +423,21 @@ void GRHydroTransformPrimToGlobalBasis(CCTK_ARGUMENTS) CCTK_INT idx2 = CCTK_VECTGFINDEX3D(cctkGH, i,j,k,1); CCTK_INT idx3 = CCTK_VECTGFINDEX3D(cctkGH, i,j,k,2); - /* Transform primitive velocity from local to global basis. - * Since velocity has contravariant index, use inverse Jacobian */ - vel[idx1] = lvel[idx1]*iJ11[idx] + lvel[idx2]*iJ12[idx] + lvel[idx3]*iJ13[idx]; - vel[idx2] = lvel[idx1]*iJ21[idx] + lvel[idx2]*iJ22[idx] + lvel[idx3]*iJ23[idx]; - vel[idx3] = lvel[idx1]*iJ31[idx] + lvel[idx2]*iJ32[idx] + lvel[idx3]*iJ33[idx]; + /* Transform primitive velocity from local to global basis. + * Since velocity has contravariant index, use inverse Jacobian */ + vel[idx1] = lvel[idx1]*iJ11[idx] + lvel[idx2]*iJ12[idx] + lvel[idx3]*iJ13[idx]; + vel[idx2] = lvel[idx1]*iJ21[idx] + lvel[idx2]*iJ22[idx] + lvel[idx3]*iJ23[idx]; + vel[idx3] = lvel[idx1]*iJ31[idx] + lvel[idx2]*iJ32[idx] + lvel[idx3]*iJ33[idx]; if(*evolve_MHD) { /* Transform primitive B-field from local to global basis. - * Since B-field has contravariant index, use inverse Jacobian */ + * Since B-field has contravariant index, use inverse Jacobian */ Bvec[idx1] = lBvec[idx1]*iJ11[idx] + lBvec[idx2]*iJ12[idx] + lBvec[idx3]*iJ13[idx]; Bvec[idx2] = lBvec[idx1]*iJ21[idx] + lBvec[idx2]*iJ22[idx] + lBvec[idx3]*iJ23[idx]; Bvec[idx3] = lBvec[idx1]*iJ31[idx] + lBvec[idx2]*iJ32[idx] + lBvec[idx3]*iJ33[idx]; } - + } } } diff --git a/src/GRHydro_UpdateMask.F90 b/src/GRHydro_UpdateMask.F90 index a679dfb..2fd47df 100644 --- a/src/GRHydro_UpdateMask.F90 +++ b/src/GRHydro_UpdateMask.F90 @@ -98,11 +98,11 @@ subroutine GRHydroPostSyncAtmosphereMask(CCTK_ARGUMENTS) !$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 ( abs(atmosphere_mask_real(i,j,k)) .gt. 0.5) then - atmosphere_mask(i,j,k) = 1 - end if - end do + do i = 1, cctk_lsh(1) + if ( abs(atmosphere_mask_real(i,j,k)) .gt. 0.5) then + atmosphere_mask(i,j,k) = 1 + end if + end do end do end do !$OMP END PARALLEL DO @@ -139,9 +139,9 @@ subroutine GRHydroCopyIntegerMask(CCTK_ARGUMENTS) !$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) - atmosphere_mask_real(i,j,k) = atmosphere_mask(i,j,k) - end do + do i = 1, cctk_lsh(1) + atmosphere_mask_real(i,j,k) = atmosphere_mask(i,j,k) + end do end do end do !$OMP END PARALLEL DO diff --git a/src/GRHydro_WENOReconstruct.F90 b/src/GRHydro_WENOReconstruct.F90 index fb03f66..9c40d95 100644 --- a/src/GRHydro_WENOReconstruct.F90 +++ b/src/GRHydro_WENOReconstruct.F90 @@ -259,12 +259,12 @@ subroutine GRHydro_WENOReconstruct1d(order, nx, v, vminus, vplus, trivial_rp, & !!$ Calculate the reconstruction do j = 1, 5 - vplus(i) = vplus(i) + wplus1 * weno_coeffs(1,j)*v(i-3+j) & - + wplus2 * weno_coeffs(2,j)*v(i-3+j) & - + wplus3 * weno_coeffs(3,j)*v(i-3+j) - vminus(i) = vminus(i) + wminus1 * weno_coeffs(3,6-j)*v(i-3+j) & - + wminus2 * weno_coeffs(2,6-j)*v(i-3+j) & - + wminus3 * weno_coeffs(1,6-j)*v(i-3+j) + vplus(i) = vplus(i) + wplus1 * weno_coeffs(1,j)*v(i-3+j) & + + wplus2 * weno_coeffs(2,j)*v(i-3+j) & + + wplus3 * weno_coeffs(3,j)*v(i-3+j) + vminus(i) = vminus(i) + wminus1 * weno_coeffs(3,6-j)*v(i-3+j) & + + wminus2 * weno_coeffs(2,6-j)*v(i-3+j) & + + wminus3 * weno_coeffs(1,6-j)*v(i-3+j) end do !vminus(i) = v(i) !vplus(i) = v(i) |