diff options
| author | bmundim <bmundim@c83d129a-5a75-4d5a-9c4d-ed3a5855bf45> | 2010-12-31 18:54:31 +0000 |
|---|---|---|
| committer | bmundim <bmundim@c83d129a-5a75-4d5a-9c4d-ed3a5855bf45> | 2010-12-31 18:54:31 +0000 |
| commit | 73d4c00c79692336f19aeb09fc948681aac421f7 (patch) | |
| tree | 58ed66c3f936468517df79cf65749e0c7f4bc254 | |
| parent | 4d03a2ddb4348a2e9caf153a37ac109956dba7f1 (diff) | |
RIT MHD development:
Revert previous commit changes to schedule.ccl
Merge MHD and GRHydro routines:
Boundaries, CalcUpdate, PPM, Reconstruct(poly), RegisterGZ,
RegisterVars
Solve a few race conditions arising in the "magnetic" routines.
This is still being tested, since it passes tests with intel
compilers but fails for gnu ones.
Alias function GRHydro_Con2PrimM_pt for C2P2CM test.
git-svn-id: http://svn.einsteintoolkit.org/cactus/EinsteinEvolve/GRHydro/trunk@202 c83d129a-5a75-4d5a-9c4d-ed3a5855bf45
| -rw-r--r-- | interface.ccl | 19 | ||||
| -rw-r--r-- | schedule.ccl | 325 | ||||
| -rw-r--r-- | src/GRHydro_Boundaries.F90 | 90 | ||||
| -rw-r--r-- | src/GRHydro_CalcUpdate.F90 | 35 | ||||
| -rw-r--r-- | src/GRHydro_Con2PrimM.F90 | 7 | ||||
| -rw-r--r-- | src/GRHydro_Con2PrimM_pt.c | 2 | ||||
| -rw-r--r-- | src/GRHydro_PPMM.F90 | 314 | ||||
| -rw-r--r-- | src/GRHydro_Prim2ConM.F90 | 12 | ||||
| -rw-r--r-- | src/GRHydro_Reconstruct.F90 | 435 | ||||
| -rw-r--r-- | src/GRHydro_ReconstructPoly.F90 | 435 | ||||
| -rw-r--r-- | src/GRHydro_RegisterGZ.cc | 176 | ||||
| -rw-r--r-- | src/GRHydro_RegisterVars.cc | 14 | ||||
| -rw-r--r-- | src/GRHydro_TmunuM.F90 | 15 | ||||
| -rw-r--r-- | src/make.code.defn | 10 |
14 files changed, 1171 insertions, 718 deletions
diff --git a/interface.ccl b/interface.ccl index ce31620..6efa839 100644 --- a/interface.ccl +++ b/interface.ccl @@ -54,6 +54,22 @@ void FUNCTION UpperMet(CCTK_REAL OUT uxx, CCTK_REAL OUT uxy, \ # CCTK_REAL IN GRHydro_rho_min, CCTK_REAL IN pmin, \ # CCTK_INT IN GRHydro_reflevel, CCTK_REAL IN GRHydro_C2P_failed) +void FUNCTION Con2PrimGenM(CCTK_INT INOUT handle, CCTK_REAL INOUT gamma_eos, CCTK_REAL INOUT dens, \ + CCTK_REAL INOUT sx, CCTK_REAL INOUT sy, CCTK_REAL INOUT sz, \ + CCTK_REAL INOUT tau, CCTK_REAL INOUT rho, \ + CCTK_REAL INOUT velx, CCTK_REAL INOUT vely, CCTK_REAL INOUT velz, \ + CCTK_REAL INOUT epsilon, CCTK_REAL INOUT pressure, \ + CCTK_REAL INOUT w_lorentz, \ + CCTK_REAL INOUT gxx, CCTK_REAL INOUT gxy, CCTK_REAL INOUT gxz, \ + CCTK_REAL INOUT gyy, CCTK_REAL INOUT gyz, CCTK_REAL INOUT gzz, \ + CCTK_REAL INOUT uxx, CCTK_REAL INOUT uxy, CCTK_REAL INOUT uxz, \ + CCTK_REAL INOUT uyy, CCTK_REAL INOUT uyz, CCTK_REAL INOUT uzz, \ + CCTK_REAL INOUT det, \ + CCTK_REAL INOUT Bvecx, CCTK_REAL INOUT Bvecy, CCTK_REAL INOUT Bvecz, \ + CCTK_REAL INOUT Bvecsq, \ + CCTK_INT OUT epsnegative, \ + CCTK_REAL OUT retval) + void FUNCTION Con2PrimPoly(CCTK_INT IN handle, CCTK_REAL OUT dens, \ CCTK_REAL OUT sx, CCTK_REAL OUT sy, \ CCTK_REAL OUT sz, CCTK_REAL OUT tau, \ @@ -124,6 +140,7 @@ PROVIDES FUNCTION SpatialDet WITH SpatialDeterminant LANGUAGE Fortran PROVIDES FUNCTION UpperMet WITH UpperMetric LANGUAGE Fortran #PROVIDES FUNCTION Con2Prim WITH Con2Prim_pt LANGUAGE Fortran PROVIDES FUNCTION Con2PrimPoly WITH Con2Prim_ptPolytype LANGUAGE Fortran +PROVIDES FUNCTION Con2PrimGenM WITH GRHydro_Con2PrimM_pt LANGUAGE Fortran PROVIDES FUNCTION Prim2ConGen WITH prim2con LANGUAGE Fortran PROVIDES FUNCTION Prim2ConPoly WITH prim2conpolytype LANGUAGE Fortran PROVIDES FUNCTION Prim2ConGenM WITH prim2conM LANGUAGE Fortran @@ -287,7 +304,7 @@ int GRHydro_eos_scalars type = SCALAR GRHydro_polytrope_handle } "Handle number for EOS" -real GRHydro_minima type = SCALAR +CCTK_REAL GRHydro_minima type = SCALAR { GRHydro_rho_min # GRHydro_dens_min diff --git a/schedule.ccl b/schedule.ccl index 5fec310..4b67d7a 100644 --- a/schedule.ccl +++ b/schedule.ccl @@ -200,18 +200,10 @@ schedule GRHydro_ParamCheck AT PARAMCHECK ### Standard symmetry registration ### ###################################### -if(CCTK_Equals(Bvec_evolution_method,"GRHydro")) +schedule GRHydro_InitSymBound AT BASEGRID { - schedule GRHydro_InitSymBoundM AT BASEGRID - { - LANG: Fortran - } "Schedule symmetries - MHD version" -} else { - schedule GRHydro_InitSymBound AT BASEGRID - { - LANG: Fortran - } "Schedule symmetries" -} + LANG: Fortran +} "Schedule symmetries" ########################################################## ### Schedule the flux weighting calculation at initial ### @@ -230,23 +222,13 @@ SCHEDULE GROUP GZPatchSystem_register { } "Tell Cactus that this group exists, but is not scheduled from here" -if(CCTK_Equals(Bvec_evolution_method,"GRHydro")) +SCHEDULE GRHydro_register_GZPatchSystem IN GZPatchSystem_register { - SCHEDULE GRHydro_register_GZPatchSystemM IN GZPatchSystem_register - { - LANG: C - OPTIONS: meta - } "register to-be-interpatch-synchronized variables with GZPatchSystem - MHD version" - -} else { - SCHEDULE GRHydro_register_GZPatchSystem IN GZPatchSystem_register - { - LANG: C - OPTIONS: meta - } "register to-be-interpatch-synchronized variables with GZPatchSystem" -} + LANG: C + OPTIONS: meta +} "register to-be-interpatch-synchronized variables with GZPatchSystem" -################################## +################################# ### Set the handle for the EOS ### ################################## @@ -415,7 +397,6 @@ if (EoS_Change) SYNC: hydrobase::eps SYNC: hydrobase::vel SYNC: hydrobase::temperature - SYNC: Y_e_con } "Reset the hydro variables if the EoS Gamma and K change between ID and evolution" } } @@ -434,18 +415,10 @@ if (EoS_Change) ### Standard registration of variables to MoL ### ################################################# -if(CCTK_Equals(Bvec_evolution_method,"GRHydro")) +schedule GRHydro_Register IN MoL_Register { - schedule GRHydro_RegisterM IN MoL_Register - { - LANG: C - } "Register variables for MoL - MHD version" -} else { - schedule GRHydro_Register IN MoL_Register - { - LANG: C - } "Register variables for MoL" -} + LANG: C +} "Register variables for MoL" #################################################### ### Setup of any scalars for efficiency purposes ### @@ -685,34 +658,18 @@ if (CCTK_Equals(method_type, "RSA FV")) if (CCTK_Equals(GRHydro_eos_type,"General")) { - if(CCTK_Equals(Bvec_evolution_method,"GRHydro")) + schedule Reconstruction IN FluxTerms AS Reconstruct { - schedule ReconstructionM IN FluxTerms AS Reconstruct - { - LANG: Fortran - } "Reconstruct the functions at the cell boundaries - MHD version" - } else { - schedule Reconstruction IN FluxTerms AS Reconstruct - { - LANG: Fortran - } "Reconstruct the functions at the cell boundaries" - } + LANG: Fortran + } "Reconstruct the functions at the cell boundaries" } else if (CCTK_Equals(GRHydro_eos_type,"Polytype")) { - if(CCTK_Equals(Bvec_evolution_method,"GRHydro")) + schedule ReconstructionPolytype IN FluxTerms AS Reconstruct { - schedule ReconstructionPolytypeM IN FluxTerms AS Reconstruct - { - LANG: Fortran - } "Reconstruct the functions at the cell boundaries - MHD version" - } else { - schedule ReconstructionPolytype IN FluxTerms AS Reconstruct - { - LANG: Fortran - } "Reconstruct the functions at the cell boundaries" - } + LANG: Fortran + } "Reconstruct the functions at the cell boundaries" } if (set_trivial_rp_grid_function) @@ -834,18 +791,10 @@ else if (CCTK_Equals(method_type, "Flux split FD")) ### After calculating the fluxes, calculate the update. ### ########################################################### -if(CCTK_Equals(Bvec_evolution_method,"GRHydro")) +schedule UpdateCalculation IN FluxTerms AFTER Riemann AS UpdateCalcul { - schedule UpdateCalculationM IN FluxTerms AFTER Riemann AS UpdateCalcul - { - LANG: Fortran - } "Calculate the update term from the fluxes - MHD version" -} else { - schedule UpdateCalculation IN FluxTerms AFTER Riemann AS UpdateCalcul - { - LANG: Fortran - } "Calculate the update term from the fluxes" -} + LANG: Fortran +} "Calculate the update term from the fluxes" ################################# ### Advance the loop counter. ### @@ -1071,18 +1020,10 @@ if (evolve_tracer) if (outflow_boundaries) { - if(CCTK_Equals(Bvec_evolution_method,"GRHydro")) + schedule GRHydro_OutflowBoundaries IN HydroBase_PostStep BEFORE GRHydro_Bound { - schedule GRHydro_OutflowBoundariesM IN HydroBase_PostStep BEFORE GRHydro_Bound - { - LANG: Fortran - } "Outflow boundaries over only some of the domain - MHD version" - } else { - schedule GRHydro_OutflowBoundaries IN HydroBase_PostStep BEFORE GRHydro_Bound - { - LANG: Fortran - } "Outflow boundaries over only some of the domain" - } + LANG: Fortran + } "Outflow boundaries over only some of the domain" } # This should not be used anymore and should be removed after some time @@ -1090,167 +1031,87 @@ schedule group Do_GRHydro_Boundaries IN HydroBase_Boundaries { } "GRHydro Boundary conditions group" -if(CCTK_Equals(Y_e_evolution_method,"GRHydro")) { - if(CCTK_Equals(Bvec_evolution_method,"GRHydro")) +if(CCTK_Equals(Y_e_evolution_method,"GRHydro")) +{ + if (evolve_tracer) { - if (evolve_tracer) + schedule GRHydro_Boundaries IN HydroBase_Select_Boundaries AS GRHydro_Bound { - schedule GRHydro_BoundariesM IN HydroBase_Select_Boundaries AS GRHydro_Bound - { - LANG: Fortran - OPTIONS: LEVEL - SYNC: dens - SYNC: tau - SYNC: scon - SYNC: w_lorentz - SYNC: HydroBase::rho - SYNC: HydroBase::press - SYNC: HydroBase::eps - SYNC: HydroBase::vel - SYNC: HydroBase::Bvec - SYNC: GRHydro_cons_tracers - SYNC: GRHydro_tracers - SYNC: hydrobase::temperature - SYNC: hydrobase::Y_e - SYNC: Y_e_con - } "Select GRHydro boundary conditions - MHD version" - } else { - schedule GRHydro_BoundariesM IN HydroBase_Select_Boundaries AS GRHydro_Bound - { - LANG: Fortran - OPTIONS: LEVEL - SYNC: dens - SYNC: tau - SYNC: scon - SYNC: w_lorentz - SYNC: HydroBase::rho - SYNC: HydroBase::press - SYNC: HydroBase::eps - SYNC: HydroBase::vel - SYNC: HydroBase::Bvec - SYNC: hydrobase::temperature - SYNC: hydrobase::Y_e - SYNC: Y_e_con - } "Select GRHydro boundary conditions - MHD version" - } + LANG: Fortran + OPTIONS: LEVEL + SYNC: dens + SYNC: tau + SYNC: scon + SYNC: w_lorentz + SYNC: HydroBase::rho + SYNC: HydroBase::press + SYNC: HydroBase::eps + SYNC: HydroBase::vel + SYNC: HydroBase::Bvec + SYNC: GRHydro_cons_tracers + SYNC: GRHydro_tracers + SYNC: hydrobase::temperature + SYNC: hydrobase::Y_e + SYNC: Y_e_con + } "Select GRHydro boundary conditions" } else { - if (evolve_tracer) + schedule GRHydro_Boundaries IN HydroBase_Select_Boundaries AS GRHydro_Bound { - schedule GRHydro_Boundaries IN HydroBase_Select_Boundaries AS GRHydro_Bound - { - LANG: Fortran - OPTIONS: LEVEL - SYNC: dens - SYNC: tau - SYNC: scon - SYNC: w_lorentz - SYNC: HydroBase::rho - SYNC: HydroBase::press - SYNC: HydroBase::eps - SYNC: HydroBase::vel - SYNC: GRHydro_cons_tracers - SYNC: GRHydro_tracers - SYNC: hydrobase::temperature - SYNC: hydrobase::Y_e - SYNC: Y_e_con - } "Select GRHydro boundary conditions" - } else - { - schedule GRHydro_Boundaries IN HydroBase_Select_Boundaries AS GRHydro_Bound - { - LANG: Fortran - OPTIONS: LEVEL - SYNC: dens - SYNC: tau - SYNC: scon - SYNC: w_lorentz - SYNC: HydroBase::rho - SYNC: HydroBase::press - SYNC: HydroBase::eps - SYNC: HydroBase::vel - SYNC: hydrobase::temperature - SYNC: hydrobase::Y_e - SYNC: Y_e_con - } "Select GRHydro boundary conditions" - } + LANG: Fortran + OPTIONS: LEVEL + SYNC: dens + SYNC: tau + SYNC: scon + SYNC: w_lorentz + SYNC: HydroBase::rho + SYNC: HydroBase::press + SYNC: HydroBase::eps + SYNC: HydroBase::vel + SYNC: HydroBase::Bvec + SYNC: hydrobase::temperature + SYNC: hydrobase::Y_e + SYNC: Y_e_con + } "Select GRHydro boundary conditions" } } else { - if(CCTK_Equals(Bvec_evolution_method,"GRHydro")) + if (evolve_tracer) { - if (evolve_tracer) + schedule GRHydro_Boundaries IN HydroBase_Select_Boundaries AS GRHydro_Bound { - schedule GRHydro_BoundariesM IN HydroBase_Select_Boundaries AS GRHydro_Bound - { - LANG: Fortran - OPTIONS: LEVEL - SYNC: dens - SYNC: tau - SYNC: scon - SYNC: w_lorentz - SYNC: HydroBase::rho - SYNC: HydroBase::press - SYNC: HydroBase::eps - SYNC: HydroBase::vel - SYNC: HydroBase::Bvec - SYNC: GRHydro_cons_tracers - SYNC: GRHydro_tracers - SYNC: hydrobase::temperature - } "Select GRHydro boundary conditions - MHD version" - } else { - schedule GRHydro_BoundariesM IN HydroBase_Select_Boundaries AS GRHydro_Bound - { - LANG: Fortran - OPTIONS: LEVEL - SYNC: dens - SYNC: tau - SYNC: scon - SYNC: w_lorentz - SYNC: HydroBase::rho - SYNC: HydroBase::press - SYNC: HydroBase::eps - SYNC: HydroBase::vel - SYNC: HydroBase::Bvec - SYNC: hydrobase::temperature - } "Select GRHydro boundary conditions - MHD version" - } + LANG: Fortran + OPTIONS: LEVEL + SYNC: dens + SYNC: tau + SYNC: scon + SYNC: w_lorentz + SYNC: HydroBase::rho + SYNC: HydroBase::press + SYNC: HydroBase::eps + SYNC: HydroBase::vel + SYNC: HydroBase::Bvec + SYNC: GRHydro_cons_tracers + SYNC: GRHydro_tracers + SYNC: hydrobase::temperature + } "Select GRHydro boundary conditions - MHD version" } else { - if (evolve_tracer) - { - schedule GRHydro_Boundaries IN HydroBase_Select_Boundaries AS GRHydro_Bound - { - LANG: Fortran - OPTIONS: LEVEL - SYNC: dens - SYNC: tau - SYNC: scon - SYNC: w_lorentz - SYNC: HydroBase::rho - SYNC: HydroBase::press - SYNC: HydroBase::eps - SYNC: HydroBase::vel - SYNC: GRHydro_cons_tracers - SYNC: GRHydro_tracers - SYNC: hydrobase::temperature - } "Select GRHydro boundary conditions" - } else + schedule GRHydro_Boundaries IN HydroBase_Select_Boundaries AS GRHydro_Bound { - schedule GRHydro_Boundaries IN HydroBase_Select_Boundaries AS GRHydro_Bound - { - LANG: Fortran - OPTIONS: LEVEL - SYNC: dens - SYNC: tau - SYNC: scon - SYNC: w_lorentz - SYNC: HydroBase::rho - SYNC: HydroBase::press - SYNC: HydroBase::eps - SYNC: HydroBase::vel - SYNC: hydrobase::temperature - } "Select GRHydro boundary conditions" - } + LANG: Fortran + OPTIONS: LEVEL + SYNC: dens + SYNC: tau + SYNC: scon + SYNC: w_lorentz + SYNC: HydroBase::rho + SYNC: HydroBase::press + SYNC: HydroBase::eps + SYNC: HydroBase::vel + SYNC: HydroBase::Bvec + SYNC: hydrobase::temperature + } "Select GRHydro boundary conditions - MHD version" } -} +} + ############################################################ ### Compute first differences of rho for mesh refinement ### ############################################################ diff --git a/src/GRHydro_Boundaries.F90 b/src/GRHydro_Boundaries.F90 index 1c97c2c..1105eef 100644 --- a/src/GRHydro_Boundaries.F90 +++ b/src/GRHydro_Boundaries.F90 @@ -21,6 +21,9 @@ #define sx(i,j,k) scon(i,j,k,1) #define sy(i,j,k) scon(i,j,k,2) #define sz(i,j,k) scon(i,j,k,3) +#define Bvecx(i,j,k) Bvec(i,j,k,1) +#define Bvecy(i,j,k) Bvec(i,j,k,2) +#define Bvecz(i,j,k) Bvec(i,j,k,3) /*@@ @routine GRHydro_InitSymBound @@ -64,6 +67,9 @@ subroutine GRHydro_InitSymBound(CCTK_ARGUMENTS) call SetCartSymVN(ierr, cctkGH, sym, "GRHydro::w_lorentz") call SetCartSymVN(ierr, cctkGH, sym, "HydroBase::eps") call SetCartSymVN(ierr, cctkGH, sym, "GRHydro::GRHydro_C2P_failed") + if(evolve_mhd.ne.0.and.clean_divergence.ne.0) then + call SetCartSymVN(ierr, cctkGH, sym, "GRHydro::psidc") + endif !!$ handle multiple tracer variables if(evolve_tracer.ne.0) then @@ -87,6 +93,8 @@ subroutine GRHydro_InitSymBound(CCTK_ARGUMENTS) call SetCartSymVN(ierr, cctkGH, sym, "HydroBase::vel[0]") call SetCartSymVN(ierr, cctkGH, sym, "GRHydro::scon[0]") + if(evolve_mhd.ne.0)call SetCartSymVN(ierr, cctkGH, sym, "HydroBase::Bvec[0]") + sym(1) = 1 sym(2) = -1 @@ -94,6 +102,7 @@ subroutine GRHydro_InitSymBound(CCTK_ARGUMENTS) call SetCartSymVN(ierr, cctkGH, sym, "HydroBase::vel[1]") call SetCartSymVN(ierr, cctkGH, sym, "GRHydro::scon[1]") + if(evolve_mhd.ne.0)call SetCartSymVN(ierr, cctkGH, sym, "HydroBase::Bvec[1]") sym(1) = 1 sym(2) = 1 @@ -101,6 +110,7 @@ subroutine GRHydro_InitSymBound(CCTK_ARGUMENTS) call SetCartSymVN(ierr, cctkGH, sym, "HydroBase::vel[2]") call SetCartSymVN(ierr, cctkGH, sym, "GRHydro::scon[2]") + if(evolve_mhd.ne.0)call SetCartSymVN(ierr, cctkGH, sym, "HydroBase::Bvec[2]") end subroutine GRHydro_InitSymBound @@ -156,6 +166,14 @@ subroutine GRHydro_Boundaries(CCTK_ARGUMENTS) "HydroBase::eps", "Flat") ierr = Boundary_SelectGroupForBC(cctkGH, faces, GRHydro_stencil, -ione, & "HydroBase::vel", "Flat") + if(evolve_mhd.ne.0) then + ierr = Boundary_SelectGroupForBC(cctkGH, faces, GRHydro_stencil, -ione, & + "HydroBase::Bvec", "Flat") + if(clean_divergence.ne.0) then + ierr = Boundary_SelectGroupForBC(cctkGH, faces, GRHydro_stencil, -ione, & + "GRHydro::psidc", "Flat") + endif + endif if(evolve_tracer.ne.0) then ierr = Boundary_SelectGroupForBC(cctkGH, faces, GRHydro_stencil, -ione, & @@ -198,6 +216,14 @@ subroutine GRHydro_Boundaries(CCTK_ARGUMENTS) "HydroBase::eps", "None") ierr = Boundary_SelectGroupForBC(cctkGH, faces, GRHydro_stencil, -ione, & "HydroBase::vel", "None") + if(evolve_mhd.ne.0) then + ierr = Boundary_SelectGroupForBC(cctkGH, faces, GRHydro_stencil, -ione, & + "HydroBase::Bvec", "None") + if(clean_divergence.ne.0) then + ierr = Boundary_SelectGroupForBC(cctkGH, faces, GRHydro_stencil, -ione, & + "GRHydro::psidc", "None") + endif + endif if(evolve_tracer.ne.0) then ierr = Boundary_SelectGroupForBC(cctkGH, faces, GRHydro_stencil, -ione, & @@ -286,17 +312,35 @@ subroutine GRHydro_OutflowBoundaries(CCTK_ARGUMENTS) eps(i,j,k) = eps(i,j,k-1) press(i,j,k) = press(i,j,k-1) w_lorentz(i,j,k) = w_lorentz(i,j,k-1) - + if(evolve_mhd.ne.0) then + Bvecx(i,j,k) = Bvecx(i,j,k-1) + Bvecy(i,j,k) = Bvecy(i,j,k-1) + Bvecz(i,j,k) = Bvecz(i,j,k-1) + if(clean_divergence.ne.0) then + psidc(i,j,k) = psidc(i,j,k-1) + endif + endif + psi4pt = 1.0d0 det = SPATIAL_DETERMINANT(gxx(i,j,k),gxy(i,j,k),gxz(i,j,k),\ gyy(i,j,k),gyz(i,j,k),gzz(i,j,k)) - call prim2con(GRHydro_eos_handle,psi4pt*gxx(i,j,k),& - psi4pt*gxy(i,j,k),psi4pt*gxz(i,j,k),& - psi4pt*gyy(i,j,k),psi4pt*gyz(i,j,k),psi4pt*gzz(i,j,k),& - det, dens(i,j,k),sx(i,j,k),sy(i,j,k),sz(i,j,k),& - tau(i,j,k),rho(i,j,k),velx(i,j,k),vely(i,j,k),velz(i,j,k),& - eps(i,j,k),press(i,j,k),w_lorentz(i,j,k)) + if(evolve_mhd.ne.0) then + call prim2conM(GRHydro_eos_handle,psi4pt*gxx(i,j,k),& + psi4pt*gxy(i,j,k),psi4pt*gxz(i,j,k),& + psi4pt*gyy(i,j,k),psi4pt*gyz(i,j,k),psi4pt*gzz(i,j,k),& + det, dens(i,j,k),sx(i,j,k),sy(i,j,k),sz(i,j,k),& + tau(i,j,k),Bvecx(i,j,k),Bvecy(i,j,k),Bvecz(i,j,k),& + rho(i,j,k),velx(i,j,k),vely(i,j,k),velz(i,j,k),& + eps(i,j,k),press(i,j,k),w_lorentz(i,j,k)) + else + call prim2con(GRHydro_eos_handle,psi4pt*gxx(i,j,k),& + psi4pt*gxy(i,j,k),psi4pt*gxz(i,j,k),& + psi4pt*gyy(i,j,k),psi4pt*gyz(i,j,k),psi4pt*gzz(i,j,k),& + det, dens(i,j,k),sx(i,j,k),sy(i,j,k),sz(i,j,k),& + tau(i,j,k),rho(i,j,k),velx(i,j,k),vely(i,j,k),velz(i,j,k),& + eps(i,j,k),press(i,j,k),w_lorentz(i,j,k)) + endif end if end do @@ -327,16 +371,34 @@ subroutine GRHydro_OutflowBoundaries(CCTK_ARGUMENTS) press(i,j,k) = press(i,j,k+1) w_lorentz(i,j,k) = w_lorentz(i,j,k+1) + if(evolve_mhd.ne.0) then + Bvecx(i,j,k) = Bvecx(i,j,k+1) + Bvecy(i,j,k) = Bvecy(i,j,k+1) + Bvecz(i,j,k) = Bvecz(i,j,k+1) + if(clean_divergence.ne.0) then + psidc(i,j,k) = psidc(i,j,k+1) + endif + endif + psi4pt = 1.0d0 det = SPATIAL_DETERMINANT(gxx(i,j,k),gxy(i,j,k),gxz(i,j,k),\ gyy(i,j,k),gyz(i,j,k),gzz(i,j,k)) - call prim2con(GRHydro_eos_handle,psi4pt*gxx(i,j,k),& - psi4pt*gxy(i,j,k),psi4pt*gxz(i,j,k),& - psi4pt*gyy(i,j,k),psi4pt*gyz(i,j,k),psi4pt*gzz(i,j,k),& - det, dens(i,j,k),sx(i,j,k),sy(i,j,k),sz(i,j,k),& - tau(i,j,k),rho(i,j,k),velx(i,j,k),vely(i,j,k),velz(i,j,k),& - eps(i,j,k),press(i,j,k),w_lorentz(i,j,k)) - + if(evolve_mhd.ne.0) then + call prim2conM(GRHydro_eos_handle,psi4pt*gxx(i,j,k),& + psi4pt*gxy(i,j,k),psi4pt*gxz(i,j,k),& + psi4pt*gyy(i,j,k),psi4pt*gyz(i,j,k),psi4pt*gzz(i,j,k),& + det, dens(i,j,k),sx(i,j,k),sy(i,j,k),sz(i,j,k),& + tau(i,j,k),Bvecx(i,j,k),Bvecy(i,j,k),Bvecz(i,j,k),& + rho(i,j,k),velx(i,j,k),vely(i,j,k),velz(i,j,k),& + eps(i,j,k),press(i,j,k),w_lorentz(i,j,k)) + else + call prim2con(GRHydro_eos_handle,psi4pt*gxx(i,j,k),& + psi4pt*gxy(i,j,k),psi4pt*gxz(i,j,k),& + psi4pt*gyy(i,j,k),psi4pt*gyz(i,j,k),psi4pt*gzz(i,j,k),& + det, dens(i,j,k),sx(i,j,k),sy(i,j,k),sz(i,j,k),& + tau(i,j,k),rho(i,j,k),velx(i,j,k),vely(i,j,k),velz(i,j,k),& + eps(i,j,k),press(i,j,k),w_lorentz(i,j,k)) + endif end if end do diff --git a/src/GRHydro_CalcUpdate.F90 b/src/GRHydro_CalcUpdate.F90 index 56f5222..81de7db 100644 --- a/src/GRHydro_CalcUpdate.F90 +++ b/src/GRHydro_CalcUpdate.F90 @@ -80,7 +80,22 @@ subroutine UpdateCalculation(CCTK_ARGUMENTS) taurhs(i,j,k) = taurhs(i,j,k) + & (alp_l * tauflux(i-xoffset,j-yoffset,k-zoffset) - & alp_r * tauflux(i,j,k)) * idx - + if(evolve_mhd.ne.0) then + Bvecrhs(i,j,k,1) = Bvecrhs(i,j,k,1) + & + (alp_l * Bvecxflux(i-xoffset,j-yoffset,k-zoffset) - & + alp_r * Bvecxflux(i,j,k)) * idx + Bvecrhs(i,j,k,2) = Bvecrhs(i,j,k,2) + & + (alp_l * Bvecyflux(i-xoffset,j-yoffset,k-zoffset) - & + alp_r * Bvecyflux(i,j,k)) * idx + Bvecrhs(i,j,k,3) = Bvecrhs(i,j,k,3) + & + (alp_l * Bveczflux(i-xoffset,j-yoffset,k-zoffset) - & + alp_r * Bveczflux(i,j,k)) * idx + if(clean_divergence.ne.0) then + psidcrhs(i,j,k) = psidcrhs(i,j,k) + & + (alp_l * psidcflux(i-xoffset,j-yoffset,k-zoffset) - & + alp_r * psidcflux(i,j,k)) * idx + endif + endif if (evolve_tracer .ne. 0) then do itracer=1,number_of_tracers @@ -217,7 +232,23 @@ subroutine UpdateCalculation(CCTK_ARGUMENTS) szflux(i,j,k)) * idx taurhs(i,j,k) = taurhs(i,j,k) + & (tauflux(i-xoffset,j-yoffset,k-zoffset) - & - tauflux(i,j,k)) * idx + tauflux(i,j,k)) * idx + if(evolve_mhd.ne.0) then + Bvecrhs(i,j,k,1) = Bvecrhs(i,j,k,1) + & + (Bvecxflux(i-xoffset,j-yoffset,k-zoffset) - & + Bvecxflux(i,j,k)) * idx + Bvecrhs(i,j,k,2) = Bvecrhs(i,j,k,2) + & + (Bvecyflux(i-xoffset,j-yoffset,k-zoffset) - & + Bvecyflux(i,j,k)) * idx + Bvecrhs(i,j,k,3) = Bvecrhs(i,j,k,3) + & + (Bveczflux(i-xoffset,j-yoffset,k-zoffset) - & + Bveczflux(i,j,k)) * idx + if(clean_divergence.ne.0) then + psidcrhs(i,j,k) = psidcrhs(i,j,k) + & + (psidcflux(i-xoffset,j-yoffset,k-zoffset) - & + psidcflux(i,j,k)) * idx + endif + endif enddo enddo diff --git a/src/GRHydro_Con2PrimM.F90 b/src/GRHydro_Con2PrimM.F90 index 984569a..7367e45 100644 --- a/src/GRHydro_Con2PrimM.F90 +++ b/src/GRHydro_Con2PrimM.F90 @@ -93,7 +93,8 @@ subroutine Conservative2PrimitiveM(CCTK_ARGUMENTS) !$OMP PARALLEL DO PRIVATE(i,j,itracer,& - !$OMP uxx, uxy, uxz, uyy, uyz, uzz, det, epsnegative) + !$OMP uxx, uxy, uxz, uyy, uyz, uzz, det, epsnegative, & + !$OMP b2,local_gam,xrho,xpress,local_K,local_pgam,sc) do k = 1, nz do j = 1, ny do i = 1, nx @@ -560,8 +561,10 @@ subroutine Conservative2PrimitivePolytypeM(CCTK_ARGUMENTS) !!$ do k = GRHydro_stencil + 1, nz - GRHydro_stencil !!$ do j = GRHydro_stencil + 1, ny - GRHydro_stencil !!$ do i = GRHydro_stencil + 1, nx - GRHydro_stencil + !$OMP PARALLEL DO PRIVATE(i,j,itracer,& - !$OMP uxx, uxy, uxz, uyy, uyz, uzz, det) + !$OMP uxx, uxy, uxz, uyy, uyz, uzz, det, epsnegative, & + !$OMP b2, xrho, xpress, local_K, local_pgam, sc) do k = 1, nz do j = 1, ny do i = 1, nx diff --git a/src/GRHydro_Con2PrimM_pt.c b/src/GRHydro_Con2PrimM_pt.c index bf888ff..f5ffd56 100644 --- a/src/GRHydro_Con2PrimM_pt.c +++ b/src/GRHydro_Con2PrimM_pt.c @@ -393,7 +393,7 @@ void CCTK_FCALL CCTK_FNAME(GRHydro_Con2PrimM_pt) ( fprintf(stdout,"velx = %26.16e \n",*velx ); fprintf(stdout,"vely = %26.16e \n",*vely ); fprintf(stdout,"velz = %26.16e \n",*velz ); - fprintf(stdout,"gam = %26.16e \n",gam ); + fprintf(stdout,"gammaeos = %26.16e \n",gammaeos ); fflush(stdout); #endif diff --git a/src/GRHydro_PPMM.F90 b/src/GRHydro_PPMM.F90 index 0ff4d72..713c24e 100644 --- a/src/GRHydro_PPMM.F90 +++ b/src/GRHydro_PPMM.F90 @@ -12,23 +12,6 @@ #include "cctk_Parameters.h" #include "GRHydro_Macros.h" -!!subroutine PPM_TVD(origm, orig, origp, bextm, bextp) -!! CCTK_REAL :: origm, orig, origp, bextm, bextp -!! CCTK_REAL :: dloc, dupw, delta -!! -!! dupw = orig - origm -!! dloc = origp - orig -!! if (dupw*dloc < 0.d0) then -!! delta=0.d0 -!! else if (abs(dupw) < abs(dloc)) then -!! delta=dupw -!! else -!! delta=dloc -!! end if -!! bextm = orig - 0.5d0 * delta -!! bextp = orig + 0.5d0 * delta -!!end subroutine PPM_TVD - /*@@ @routine SimplePPM_1dM @date Thu Feb 14 19:08:52 2002 @@ -52,9 +35,9 @@ subroutine SimplePPM_1dM(handle,poly,nx,dx,& - rho,velx,vely,velz,Bvcx,Bvcy,Bvcz,eps,press,& - rhominus,velxminus,velyminus,velzminus,Bvcxminus,Bvcyminus,Bvczminus,epsminus,& - rhoplus,velxplus,velyplus,velzplus,Bvcxplus,Bvcyplus,Bvczplus,epsplus,& + rho,velx,vely,velz,Bvcx,Bvcy,Bvcz,psidc,eps,press,& + rhominus,velxminus,velyminus,velzminus,Bvcxminus,Bvcyminus,Bvczminus,psidcminus,epsminus,& + rhoplus,velxplus,velyplus,velzplus,Bvcxplus,Bvcyplus,Bvczplus,psidcplus,epsplus,dc_flag,& trivial_rp, hydro_excision_mask,& gxx, gxy, gxz, gyy, gyz, gzz, psi4, beta, alp, w_lorentz, & dir, ni, nj, nrx, nry, nrz, ev_l, ev_r, xw) @@ -69,27 +52,17 @@ subroutine SimplePPM_1dM(handle,poly,nx,dx,& CCTK_INT :: handle,poly,nx CCTK_REAL :: dx CCTK_REAL, dimension(nx) :: rho,velx,vely,velz,eps - CCTK_REAL, dimension(nx) :: Bvcx,Bvcy,Bvcz + CCTK_REAL, dimension(nx) :: Bvcx,Bvcy,Bvcz,psidc CCTK_REAL, dimension(nx) :: rhominus,velxminus,velyminus,velzminus,epsminus - CCTK_REAL, dimension(nx) :: Bvcxminus,Bvcyminus,Bvczminus + CCTK_REAL, dimension(nx) :: Bvcxminus,Bvcyminus,Bvczminus,psidcminus CCTK_REAL, dimension(nx) :: rhoplus,velxplus,velyplus,velzplus,epsplus - CCTK_REAL, dimension(nx) :: Bvcxplus,Bvcyplus,Bvczplus - CCTK_REAL, dimension(nx) :: rhominusl,velxminusl,velyminusl,velzminusl - CCTK_REAL, dimension(nx) :: Bvcxminusl,Bvcyminusl,Bvczminusl - CCTK_REAL, dimension(nx) :: epsminusl - CCTK_REAL, dimension(nx) :: rhoplusl,velxplusl,velyplusl,velzplusl,epsplusl - CCTK_REAL, dimension(nx) :: Bvcxplusl,Bvcyplusl,Bvczplusl - CCTK_REAL, dimension(nx) :: rhominusr,velxminusr,velyminusr,velzminusr - CCTK_REAL, dimension(nx) :: Bvcxminusr,Bvcyminusr,Bvczminusr - CCTK_REAL, dimension(nx) :: epsminusr - CCTK_REAL, dimension(nx) :: rhoplusr,velxplusr,velyplusr,velzplusr,epsplusr - CCTK_REAL, dimension(nx) :: Bvcxplusr,Bvcyplusr,Bvczplusr - - CCTK_INT :: i,s + CCTK_REAL, dimension(nx) :: Bvcxplus,Bvcyplus,Bvczplus,psidcplus + + CCTK_INT :: i,s,dc_flag CCTK_REAL, dimension(nx) :: drho,dvelx,dvely,dvelz,deps - CCTK_REAL, dimension(nx) :: dBvcx,dBvcy,dBvcz + CCTK_REAL, dimension(nx) :: dBvcx,dBvcy,dBvcz,dpsidc CCTK_REAL, dimension(nx) :: dmrho,dmvelx,dmvely,dmvelz,dmeps - CCTK_REAL, dimension(nx) :: dmBvcx,dmBvcy,dmBvcz + CCTK_REAL, dimension(nx) :: dmBvcx,dmBvcy,dmBvcz,dmpsidc CCTK_REAL, dimension(nx) :: press,dpress,d2rho,tilde_flatten CCTK_REAL :: dpress2,dvel,w,flatten,eta,etatilde @@ -126,6 +99,8 @@ trivial_rp = .true. dBvcx(i) = 0.5d0 * (Bvcx(i+1) - Bvcx(i-1)) dBvcy(i) = 0.5d0 * (Bvcy(i+1) - Bvcy(i-1)) dBvcz(i) = 0.5d0 * (Bvcz(i+1) - Bvcz(i-1)) + if(dc_flag.ne.0)dpsidc(i) = 0.5d0 * (psidc(i+1) - psidc(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 @@ -155,6 +130,10 @@ trivial_rp = .true. STEEP(Bvcx, dBvcx, dmBvcx) STEEP(Bvcy, dBvcy, dmBvcy) STEEP(Bvcz, dBvcz, dmBvcz) + if(dc_flag.ne.0) then + STEEP(psidc,dpsidc,dmpsidc) + endif + end do if (poly .eq. 0) then do i = 2, nx - 1 @@ -187,6 +166,11 @@ trivial_rp = .true. Bvczplus(i) = 0.5d0 * (Bvcz(i) + Bvcz(i+1)) + & (dmBvcz(i) - dmBvcz(i+1)) / 6.d0 Bvczminus(i+1) = Bvczplus(i) + if(dc_flag.ne.0) then + psidcplus(i) = 0.5d0 * (psidc(i) + psidc(i+1)) + & + (dmpsidc(i) - dmpsidc(i+1)) / 6.d0 + psidcminus(i+1) = psidcplus(i) + endif end do if (poly .eq. 0) then do i = 2, nx-2 @@ -283,6 +267,10 @@ trivial_rp = .true. LEFTMINUS(Bvcy, Bvcyminus) LEFTPLUS(Bvcz, Bvczplus) LEFTMINUS(Bvcz, Bvczminus) + if(dc_flag.ne.0) then + LEFTPLUS(psidc, psidcplus) + LEFTMINUS(psidc, psidcminus) + endif if (poly .eq. 0) then LEFTPLUS(eps, epsplus) LEFTMINUS(eps, epsminus) @@ -302,6 +290,10 @@ trivial_rp = .true. RIGHTMINUS(Bvcy, Bvcyminus) RIGHTPLUS(Bvcz, Bvczplus) RIGHTMINUS(Bvcz, Bvczminus) + if(dc_flag.ne.0) then + RIGHTPLUS(psidc, psidcplus) + RIGHTMINUS(psidc, psidcminus) + endif if (poly .eq. 0) then RIGHTPLUS(eps, epsplus) RIGHTMINUS(eps, epsminus) @@ -321,6 +313,10 @@ trivial_rp = .true. CHECK(Bvcy, Bvcyminus(i+1)) CHECK(Bvcz, Bvczplus(i)) CHECK(Bvcz, Bvczminus(i+1)) + if(dc_flag.ne.0) then + CHECK(psidc, psidcplus(i)) + CHECK(psidc, psidcminus(i+1)) + endif if (poly .eq. 0) then CHECK(eps, epsplus(i)) CHECK(eps, epsminus(i+1)) @@ -391,6 +387,10 @@ trivial_rp = .true. Bvcyminus(i) = flatten * Bvcyminus(i) + (1.d0 - flatten) * Bvcy(i) Bvczplus(i) = flatten * Bvczplus(i) + (1.d0 - flatten) * Bvcz(i) Bvczminus(i) = flatten * Bvczminus(i) + (1.d0 - flatten) * Bvcz(i) + if(dc_flag.ne.0) then + psidcplus(i) = flatten * psidcplus(i) + (1.d0 - flatten) * psidc(i) + psidcminus(i) = flatten * psidcminus(i) + (1.d0 - flatten) * psidc(i) + endif 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) @@ -418,6 +418,10 @@ trivial_rp = .true. Bvcyminus(i) = flatten * Bvcyminus(i) + (1.d0 - flatten) * Bvcy(i) Bvczplus(i) = flatten * Bvczplus(i) + (1.d0 - flatten) * Bvcz(i) Bvczminus(i) = flatten * Bvczminus(i) + (1.d0 - flatten) * Bvcz(i) + if(dc_flag.ne.0) then + psidcplus(i) = flatten * psidcplus(i) + (1.d0 - flatten) * psidc(i) + psidcminus(i) = flatten * psidcminus(i) + (1.d0 - flatten) * psidc(i) + endif 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) @@ -461,6 +465,9 @@ do i = GRHydro_stencil, nx - GRHydro_stencil + 1 MON(Bvcxminus,Bvcx,Bvcxplus) MON(Bvcyminus,Bvcy,Bvcyplus) MON(Bvczminus,Bvcz,Bvczplus) + if(dc_flag.ne.0) then + MON(psidcminus,psidc,psidcplus) + endif end do if (poly .eq. 0) then do i = GRHydro_stencil, nx - GRHydro_stencil + 1 @@ -517,7 +524,13 @@ do i = GRHydro_stencil, nx - GRHydro_stencil + 1 Bvcyplus(i-1)=Bvcy(i) Bvczminus(i-1)=Bvcz(i) Bvczplus(i-1)=Bvcz(i) - if (poly .eq. 0) then + if (dc_flag.ne.0) then + psidcminus(i)=psidc(i) + psidcplus(i)=psidc(i) + psidcminus(i-1)=psidc(i) + psidcplus(i-1)=psidc(i) + end if + if (poly .eq. 0) then epsminus(i)=eps(i) epsplus(i)=eps(i) epsminus(i-1)=eps(i) @@ -536,6 +549,9 @@ do i = GRHydro_stencil, nx - GRHydro_stencil + 1 call PPM_TVD(Bvcx(i-1), Bvcx(i), Bvcx(i+1), Bvcxminus(i), Bvcxplus(i)) call PPM_TVD(Bvcy(i-1), Bvcy(i), Bvcy(i+1), Bvcyminus(i), Bvcyplus(i)) call PPM_TVD(Bvcz(i-1), Bvcz(i), Bvcz(i+1), Bvczminus(i), Bvczplus(i)) + if(dc_flag.ne.0) then + call PPM_TVD(psidc(i-1), psidc(i), psidc(i+1), psidcminus(i), psidcplus(i)) + endif if (poly .eq. 0) then call PPM_TVD(eps(i-1), eps(i), eps(i+1), epsminus(i), epsplus(i)) end if @@ -574,6 +590,12 @@ do i = GRHydro_stencil, nx - GRHydro_stencil + 1 Bvcyplus(i+1)=Bvcy(i) Bvczminus(i+1)=Bvcz(i) Bvczplus(i+1)=Bvcz(i) + if (dc_flag.ne.0) then + psidcminus(i)=psidc(i) + psidcplus(i)=psidc(i) + psidcminus(i+1)=psidc(i) + psidcplus(i+1)=psidc(i) + endif if (poly .eq. 0) then epsminus(i)=eps(i) epsplus(i)=eps(i) @@ -593,6 +615,9 @@ do i = GRHydro_stencil, nx - GRHydro_stencil + 1 call PPM_TVD(Bvcx(i-1), Bvcx(i), Bvcx(i+1), Bvcxminus(i), Bvcxplus(i)) call PPM_TVD(Bvcy(i-1), Bvcy(i), Bvcy(i+1), Bvcyminus(i), Bvcyplus(i)) call PPM_TVD(Bvcz(i-1), Bvcz(i), Bvcz(i+1), Bvczminus(i), Bvczplus(i)) + if(dc_flag.ne.0) then + call PPM_TVD(Bvcz(i-1), Bvcz(i), Bvcz(i+1), Bvczminus(i), Bvczplus(i)) + endif if (poly .eq. 0) then call PPM_TVD(eps(i-1), eps(i), eps(i+1), epsminus(i), epsplus(i)) end if @@ -605,216 +630,3 @@ return end subroutine SimplePPM_1dM - -!!subroutine SimplePPM_tracer_1d(nx,dx,rho,velx,vely,velz, & -!! tracer,tracerminus,tracerplus,press) -!! -!! USE GRHydro_Scalars -!! -!! implicit none -!! -!! DECLARE_CCTK_PARAMETERS -!! -!! CCTK_INT :: nx -!! CCTK_REAL :: dx -!! CCTK_REAL, dimension(nx) :: rho,velx,vely,velz -!! CCTK_REAL, dimension(nx,number_of_tracers) :: tracer,tracerminus,tracerplus -!! CCTK_REAL :: tracerflatomega -!! -!! -!! CCTK_INT :: i,s,itracer -!! CCTK_REAL, dimension(nx) :: press,dpress,tilde_flatten -!! CCTK_REAL, dimension(nx,number_of_tracers) :: dmtracer, dtracer, tracerflat!, d2tracer -!! CCTK_REAL :: dpress2,w,flatten,dvel -!! CCTK_REAL :: eta, etatilde -!! -!!!!$ Average slopes delta_m(a). See (1.7) of Colella and Woodward, p.178 -!!!!$ This is the expression for an even grid. -!! -!! -!! do i = 2, nx - 1 -!! dpress(i) = press(i+1) - press(i-1) -!! end do -!! -!! do itracer=1,number_of_tracers -!! do i = 2, nx - 1 -!! dtracer(i,itracer) = 0.5d0 * (tracer(i+1,itracer) - tracer(i-1,itracer)) -!!! d2tracer(i,itracer) = (tracer(i+1) - 2.d0 * tracer(i) + tracer(i-1))! / 6.d0 / dx / dx -!!! ! since we use d2tracer only for the condition d2tracer(i+1)*d2tracer(i-1)<0 -!!! ! the denominator is not necessary -!! enddo -!! enddo -!! -!!!!$ Steepened slope. See (1.8) of Colella and Woodward, p.178 -!! -!! do itracer=1,number_of_tracers -!! do i = 2, nx - 1 -!! if( (tracer(i+1,itracer) - tracer(i,itracer)) * & -!! (tracer(i,itracer) - tracer(i-1,itracer)) > 0.0d0 ) then -!! dmtracer(i,itracer) = sign(1.0d0,dtracer(i,itracer)) * & -!! min(abs(dtracer(i,itracer)), 2.0d0 * & -!! abs(tracer(i,itracer) - tracer(i-1,itracer)), & -!! 2.0d0 * abs(tracer(i+1,itracer) - tracer(i,itracer))) -!! else -!! dmtracer(i,itracer) = 0.0d0 -!! endif -!! end do -!! enddo -!! -!!!!$ Initial boundary states. See (1.9) of Colella and Woodward, p.178 -!! -!! do itracer=1,number_of_tracers -!! do i = 2, nx - 2 -!! tracerplus(i,itracer) = 0.5d0 * (tracer(i,itracer) + tracer(i+1,itracer)) + & -!! (dmtracer(i,itracer) - dmtracer(i+1,itracer)) / 6.d0 -!! tracerminus(i+1,itracer) = tracerplus(i,itracer) -!! enddo -!! enddo -!! -!! -!!!!$Discontinuity steepening. See (1.14-17) of C&W. -!!!!$This is the detect routine which mat be activated with the ppm_detect parameter -!!!!$Note that this part really also depends on the grid being even. -!!!!$Note also that we do not have access to the gas constant gamma. -!!!!$So this is just dropped from eq. (3.2) of C&W. -!!!!$We can get around this by just rescaling the constant k0 (ppm_k0 here). -!! -!!!!! We might play around with this for the tracer. CURRENTLY TURNED OFF -!! -!!#if 0 -!! if (ppm_detect .eq. 1000) then -!! do itracer=1,number_of_tracers -!! -!! do i = 3, nx - 2 -!! if ( (dtracer(i+1,itracer)*dtracer(i-1,itracer) > 0.d0) & !make sure this is not an extremum -!! .and.(abs(tracer(i+1,itracer)-tracer(i-1,itracer)) - & !this is to prevent steepening -!! !of relatively small composition jumps -!! ppm_epsilon_shock * min(tracer(i+1,itracer), tracer(i-1,itracer)) > 0.d0 ) & -!! .and. & ! the actual criterion from Plewa & Mueller -!! ((tracer(i+1,itracer)-tracer(i-1,itracer)) / & -!! (tracer(i+2,itracer)-tracer(i-2,itracer)) > ppm_omega1 ) ) then -!! -!! etatilde = (tracer(i-2,itracer) - tracer(i+2,itracer) + & -!! 4.d0 * dtracer(i,itracer)) / (dtracer(i,itracer) * 12.d0) -!! -!! write(*,*) "Additional Steepening in Zone: ",i -!! -!! else -!! etatilde = 0.d0 -!! end if -!! eta = max(0.d0, min(1.d0, ppm_eta1 * (etatilde - ppm_eta2))) -!! if (ppm_k0 * abs(dtracer(i,itracer)) * min(press(i-1),press(i+1)) < & -!! abs(dpress(i)) * min(tracer(i-1,itracer), tracer(i+1,itracer))) then -!! eta = 0.d0 -!! end if -!! tracerminus(i,itracer) = tracerminus(i,itracer) * (1.d0 - eta) + & -!! (tracer(i-1,itracer) + 0.5d0 * dmtracer(i-1,itracer)) * eta -!! tracerplus(i,itracer) = tracerplus(i,itracer) * (1.d0 - eta) + & -!! (tracer(i+1,itracer) - 0.5d0 * dmtracer(i+1,itracer)) * eta -!! end do -!! -!! enddo -!! -!! end if -!!#endif -!! -!!!!$ Zone flattening. See appendix of C&W, p. 197-8. -!! -!! do i = 3, nx - 2 -!! dpress2 = press(i+2) - press(i-2) -!! dvel = velx(i+1) - velx(i-1) -!! if ( (abs(dpress(i)) > ppm_epsilon * min(press(i-1),press(i+1))) .and. & -!! (dvel < 0.d0) ) then -!! w = 1.d0 -!! else -!! w = 0.d0 -!! end if -!! if (abs(dpress2) < ppm_small) then -!! tilde_flatten(i) = 1.d0 -!! else -!! tilde_flatten(i) = max(0.d0, 1.d0 - w * max(0.d0, ppm_omega2 * & -!! (dpress(i) / dpress2 - ppm_omega1))) -!! end if -!! end do -!! -!! if (PPM3) then -!! do itracer=1,number_of_tracers -!! do i = 3, nx - 2 -!! flatten = tilde_flatten(i) -!! tracerplus(i,itracer) = flatten * tracerplus(i,itracer) & -!! + (1.d0 - flatten) * tracer(i,itracer) -!! tracerminus(i,itracer) = flatten * tracerminus(i,itracer) & -!! + (1.d0 - flatten) * tracer(i,itracer) -!! end do -!! enddo -!! else !!$ Really implement C&W, page 197; which requires stencil 4. -!! do itracer=1,number_of_tracers -!! do i = 4, nx - 3 -!! s=sign(1.d0, -dpress(i)) -!! flatten = max(tilde_flatten(i), tilde_flatten(i+s)) -!! tracerplus(i,itracer) = flatten * tracerplus(i,itracer) + & -!! (1.d0 - flatten) * tracer(i,itracer) -!! tracerminus(i,itracer) = flatten * tracerminus(i,itracer) & -!! + (1.d0 - flatten) * tracer(i,itracer) -!! end do -!! enddo -!! end if -!! -!! -!!!! Additional flattening a la Plewa & Mueller -!! -!!#if 1 -!! do itracer=1,number_of_tracers -!! do i = 2, nx - 1 -!! if ( ( tracer(i+1,itracer) - tracer(i,itracer) ) * & -!! ( tracer(i,itracer) - tracer(i-1,itracer) ) < 0.0d0 ) then -!! tracerflat(i,itracer) = 1.0d0 -!! else -!! tracerflat(i,itracer) = 0.0d0 -!! endif -!! enddo -!! enddo -!! -!! do itracer=1,number_of_tracers -!! do i = 3, nx -2 -!! -!! tracerflatomega = 0.5d0 * max(tracerflat(i-1,itracer),2.0d0*tracerflat(i,itracer), & -!! tracerflat(i+1,itracer)) * ppm_omega_tracer -!! -!! tracerplus(i,itracer) = tracerflatomega*tracer(i,itracer) + & -!! (1.0d0 - tracerflatomega)*tracerplus(i,itracer) -!! -!! tracerminus(i,itracer) = tracerflatomega*tracer(i,itracer) + & -!! (1.0d0 - tracerflatomega)*tracerminus(i,itracer) -!! -!! enddo -!! enddo -!! -!! -!!#endif -!! -!!!!$ Monotonicity. See (1.10) of C&W. -!! -!! -!! do itracer=1,number_of_tracers -!! do i = GRHydro_stencil, nx - GRHydro_stencil + 1 -!! if (((tracerplus(i,itracer)-tracer(i,itracer))* & -!! (tracer(i,itracer)-tracerminus(i,itracer)) .le. 0.d0)) then -!! tracerminus(i,itracer) = tracer(i,itracer) -!! tracerplus(i,itracer) = tracer(i,itracer) -!! else if ((tracerplus(i,itracer) - tracerminus(i,itracer)) * (tracer(i,itracer) - 0.5d0 * & -!! (tracerplus(i,itracer) + tracerminus(i,itracer))) > & -!! (tracerplus(i,itracer) - tracerminus(i,itracer))**2 / 6.d0) then -!! tracerminus(i,itracer) = 3.d0 * tracer(i,itracer) - 2.d0 * tracerplus(i,itracer) -!! else if ((tracerplus(i,itracer) - tracerminus(i,itracer)) * (tracer(i,itracer) - 0.5d0 * & -!! (tracerplus(i,itracer) + tracerminus(i,itracer))) < & -!! -(tracerplus(i,itracer) - tracerminus(i,itracer))**2 / 6.d0 ) then -!! tracerplus(i,itracer) = 3.d0 * tracer(i,itracer) - 2.d0 * tracerminus(i,itracer) -!! end if -!! end do -!! enddo -!! -!! -!! -!!end subroutine SimplePPM_tracer_1d -!! diff --git a/src/GRHydro_Prim2ConM.F90 b/src/GRHydro_Prim2ConM.F90 index 840a181..90c4644 100644 --- a/src/GRHydro_Prim2ConM.F90 +++ b/src/GRHydro_Prim2ConM.F90 @@ -116,23 +116,25 @@ subroutine prim2conM(handle, gxx, gxy, gxz, gyy, gyz, gzz, det, ddens, & implicit none DECLARE_CCTK_PARAMETERS + DECLARE_CCTK_FUNCTIONS CCTK_REAL :: gxx, gxy, gxz, gyy, gyz, gzz, det - CCTK_REAL :: ddens, dsx, dsy, dsz, dtau, dBvcx, dBvcy, dBvcz, & + CCTK_REAL, dimension(1) :: ddens, dsx, dsy, dsz, dtau, dBvcx, dBvcy, dBvcz, & drho, dvelx, dvely, dvelz,& - deps, dpress, w, vlowx, vlowy, vlowz - CCTK_REAL :: Bdotv,ab0,b2,blowx,blowy,blowz + deps, dpress, vlowx, vlowy, vlowz + CCTK_REAL :: w + CCTK_REAL, dimension(1) :: Bdotv,ab0,b2,blowx,blowy,blowz CCTK_INT :: handle character(len=256) NaN_WarnLine ! begin EOS Omni vars integer :: n, keytemp, anyerr, keyerr(1) - real*8 :: xpress,xeps,xtemp,xye + real*8 :: xpress(1),xeps(1),xtemp(1),xye(1) n = 1; keytemp = 0; anyerr = 0; keyerr(1) = 0 xpress = 0.0d0; xeps = 0.0d0; xtemp = 0.0d0; xye = 0.0d0 ! end EOS Omni vars - w = 1.d0 / sqrt(1.d0 - DOT2(dvelx,dvely,dvelz)) + w = 1.d0 / sqrt(1.d0 - DOT2(dvelx(1),dvely(1),dvelz(1))) !!$ BEGIN: Check for NaN value if (w .ne. w) then diff --git a/src/GRHydro_Reconstruct.F90 b/src/GRHydro_Reconstruct.F90 index 9168147..23ec0a2 100644 --- a/src/GRHydro_Reconstruct.F90 +++ b/src/GRHydro_Reconstruct.F90 @@ -20,6 +20,9 @@ #define sx(i,j,k) scon(i,j,k,1) #define sy(i,j,k) scon(i,j,k,2) #define sz(i,j,k) scon(i,j,k,3) +#define Bvecx(i,j,k) Bvec(i,j,k,1) +#define Bvecy(i,j,k) Bvec(i,j,k,2) +#define Bvecz(i,j,k) Bvec(i,j,k,3) /*@@ @@ -56,7 +59,7 @@ subroutine Reconstruction(CCTK_ARGUMENTS) &type_bitsz, trivialz, not_trivialz CCTK_REAL, dimension(:,:,:),allocatable :: & - &psi4, lbetax, lbetay, lbetaz + &psi4, lbetax, lbetay, lbetaz, dum, dump, dumm !!$ CCTK_REAL, dimension(cctk_lsh(1),cctk_lsh(2),cctk_lsh(3)) :: & !!$ &psi4, lbetax, lbetay, lbetaz @@ -69,10 +72,15 @@ subroutine Reconstruction(CCTK_ARGUMENTS) call CCTK_WARN(0, "Allocation problems with trivial_rp") end if +!!$ The dum variables are used as dummies if MHD on but divergence cleaning off allocate(psi4(cctk_lsh(1),cctk_lsh(2),cctk_lsh(3)),& lbetax(cctk_lsh(1),cctk_lsh(2),cctk_lsh(3)),& lbetay(cctk_lsh(1),cctk_lsh(2),cctk_lsh(3)),& - lbetaz(cctk_lsh(1),cctk_lsh(2),cctk_lsh(3)),STAT=ierr) + lbetaz(cctk_lsh(1),cctk_lsh(2),cctk_lsh(3)),& + dum(cctk_lsh(1),cctk_lsh(2),cctk_lsh(3)),& + dump(cctk_lsh(1),cctk_lsh(2),cctk_lsh(3)),& + dumm(cctk_lsh(1),cctk_lsh(2),cctk_lsh(3)),STAT=ierr) + if (ierr .ne. 0) then call CCTK_WARN(0, "Allocation problems with lbeta") end if @@ -126,6 +134,19 @@ subroutine Reconstruction(CCTK_ARGUMENTS) velzplus = 0.0d0 velzminus = 0.0d0 + if(evolve_mhd.ne.0) then + Bvecxplus = 0.0d0 + Bvecxminus = 0.0d0 + Bvecyplus = 0.0d0 + Bvecyminus = 0.0d0 + Bveczplus = 0.0d0 + Bveczminus = 0.0d0 + if(clean_divergence.ne.0) then + psidcplus = 0.0d0 + psidcminus = 0.0d0 + endif + endif + if (evolve_tracer .ne. 0) then tracerplus = 0.0d0 tracerminus = 0.0d0 @@ -182,6 +203,20 @@ subroutine Reconstruction(CCTK_ARGUMENTS) call CCTK_WARN(0, "Variable type to reconstruct not recognized.") end if +!!$ B-field is both prim and con + if(evolve_mhd.ne.0) then + call tvdreconstruct(nx, ny, nz, xoffset, yoffset, zoffset, & + Bvec(:,:,:,1), Bvecxplus, Bvecxminus, trivial_rp, hydro_excision_mask) + call tvdreconstruct(nx, ny, nz, xoffset, yoffset, zoffset, & + Bvec(:,:,:,2), Bvecyplus, Bvecyminus, trivial_rp, hydro_excision_mask) + call tvdreconstruct(nx, ny, nz, xoffset, yoffset, zoffset, & + Bvec(:,:,:,3), Bveczplus, Bveczminus, trivial_rp, hydro_excision_mask) + if(clean_divergence.ne.0) then + call tvdreconstruct(nx, ny, nz, xoffset, yoffset, zoffset, & + psidc, psidcplus, psidcminus, trivial_rp, hydro_excision_mask) + endif + endif + !$OMP PARALLEL DO PRIVATE(i, j) do k = 1, nz do j = 1, ny @@ -211,31 +246,91 @@ subroutine Reconstruction(CCTK_ARGUMENTS) else if (CCTK_EQUALS(recon_method,"ppm")) then if (flux_direction == 1) then - !$OMP PARALLEL DO PRIVATE(i, j) - do k = GRHydro_stencil, nz - GRHydro_stencil + 1 - do j = GRHydro_stencil, ny - GRHydro_stencil + 1 - call SimplePPM_1d(GRHydro_eos_handle,0,nx,CCTK_DELTA_SPACE(1),& - rho(:,j,k),velx(:,j,k),vely(:,j,k),velz(:,j,k),eps(:,j,k),& - press(:,j,k),rhominus(:,j,k),velxminus(:,j,k),velyminus(:,j,k),& - velzminus(:,j,k),epsminus(:,j,k),rhoplus(:,j,k),& - velxplus(:,j,k),velyplus(:,j,k),velzplus(:,j,k),epsplus(:,j,k),& - trivial_rp(:,j,k), hydro_excision_mask(:,j,k),& - gxx(:,j,k), gxy(:,j,k), gxz(:,j,k), gyy(:,j,k), gyz(:,j,k), & - gzz(:,j,k), psi4(:,j,k), lbetax(:,j,k), alp(:,j,k),& - w_lorentz(:,j,k), & - 1, j, k, nx, ny, nz, GRHydro_mppm_eigenvalue_x_left, & - GRHydro_mppm_eigenvalue_x_right, & - GRHydro_mppm_xwind) - do i = 1, nx - if (trivial_rp(i,j,k)) then - SpaceMask_SetStateBitsF90(space_mask, i, j, k, type_bitsx, trivialx) - else - SpaceMask_SetStateBitsF90(space_mask, i, j, k, type_bitsx, not_trivialx) - end if - end do - end do - end do - !$OMP END PARALLEL DO + if(evolve_mhd.ne.0) then + if(clean_divergence.ne.0) then + !$OMP PARALLEL DO PRIVATE(i, j) + do k = GRHydro_stencil, nz - GRHydro_stencil + 1 + do j = GRHydro_stencil, ny - GRHydro_stencil + 1 + call SimplePPM_1dM(GRHydro_eos_handle,0,nx,CCTK_DELTA_SPACE(1),& + rho(:,j,k),velx(:,j,k),vely(:,j,k),velz(:,j,k),& + Bvecx(:,j,k),Bvecy(:,j,k),Bvecz(:,j,k),psidc(:,j,k),eps(:,j,k),press(:,j,k),& + rhominus(:,j,k),velxminus(:,j,k),velyminus(:,j,k),velzminus(:,j,k),& + Bvecxminus(:,j,k),Bvecyminus(:,j,k),Bveczminus(:,j,k),psidcminus(:,j,k),epsminus(:,j,k),& + rhoplus(:,j,k),velxplus(:,j,k),velyplus(:,j,k),velzplus(:,j,k),& + Bvecxplus(:,j,k),Bvecyplus(:,j,k),Bveczplus(:,j,k),psidcplus(:,j,k),epsplus(:,j,k),& + 1,trivial_rp(:,j,k), hydro_excision_mask(:,j,k),& + gxx(:,j,k), gxy(:,j,k), gxz(:,j,k), gyy(:,j,k), gyz(:,j,k), & + gzz(:,j,k), psi4(:,j,k), lbetax(:,j,k), alp(:,j,k),& + w_lorentz(:,j,k), & + 1, j, k, nx, ny, nz, GRHydro_mppm_eigenvalue_x_left, & + GRHydro_mppm_eigenvalue_x_right, & + GRHydro_mppm_xwind) + do i = 1, nx + if (trivial_rp(i,j,k)) then + SpaceMask_SetStateBitsF90(space_mask, i, j, k, type_bitsx, trivialx) + else + SpaceMask_SetStateBitsF90(space_mask, i, j, k, type_bitsx, not_trivialx) + end if + end do + end do + end do + !$OMP END PARALLEL DO + else !clean_divergence + !$OMP PARALLEL DO PRIVATE(i, j) + do k = GRHydro_stencil, nz - GRHydro_stencil + 1 + do j = GRHydro_stencil, ny - GRHydro_stencil + 1 + call SimplePPM_1dM(GRHydro_eos_handle,0,nx,CCTK_DELTA_SPACE(1),& + rho(:,j,k),velx(:,j,k),vely(:,j,k),velz(:,j,k),& + Bvecx(:,j,k),Bvecy(:,j,k),Bvecz(:,j,k),dum(:,j,k),eps(:,j,k),press(:,j,k),& + rhominus(:,j,k),velxminus(:,j,k),velyminus(:,j,k),velzminus(:,j,k),& + Bvecxminus(:,j,k),Bvecyminus(:,j,k),Bveczminus(:,j,k),dumm(:,j,k),epsminus(:,j,k),& + rhoplus(:,j,k),velxplus(:,j,k),velyplus(:,j,k),velzplus(:,j,k),& + Bvecxplus(:,j,k),Bvecyplus(:,j,k),Bveczplus(:,j,k),dump(:,j,k),epsplus(:,j,k),& + 0,trivial_rp(:,j,k), hydro_excision_mask(:,j,k),& + gxx(:,j,k), gxy(:,j,k), gxz(:,j,k), gyy(:,j,k), gyz(:,j,k), & + gzz(:,j,k), psi4(:,j,k), lbetax(:,j,k), alp(:,j,k),& + w_lorentz(:,j,k), & + 1, j, k, nx, ny, nz, GRHydro_mppm_eigenvalue_x_left, & + GRHydro_mppm_eigenvalue_x_right, & + GRHydro_mppm_xwind) + do i = 1, nx + if (trivial_rp(i,j,k)) then + SpaceMask_SetStateBitsF90(space_mask, i, j, k, type_bitsx, trivialx) + else + SpaceMask_SetStateBitsF90(space_mask, i, j, k, type_bitsx, not_trivialx) + end if + end do + end do + end do + !$OMP END PARALLEL DO + endif !clean_divergence + else !evolve_mhd + !$OMP PARALLEL DO PRIVATE(i, j) + do k = GRHydro_stencil, nz - GRHydro_stencil + 1 + do j = GRHydro_stencil, ny - GRHydro_stencil + 1 + call SimplePPM_1d(GRHydro_eos_handle,0,nx,CCTK_DELTA_SPACE(1),& + rho(:,j,k),velx(:,j,k),vely(:,j,k),velz(:,j,k),eps(:,j,k),& + press(:,j,k),rhominus(:,j,k),velxminus(:,j,k),velyminus(:,j,k),& + velzminus(:,j,k),epsminus(:,j,k),rhoplus(:,j,k),& + velxplus(:,j,k),velyplus(:,j,k),velzplus(:,j,k),epsplus(:,j,k),& + trivial_rp(:,j,k), hydro_excision_mask(:,j,k),& + gxx(:,j,k), gxy(:,j,k), gxz(:,j,k), gyy(:,j,k), gyz(:,j,k), & + gzz(:,j,k), psi4(:,j,k), lbetax(:,j,k), alp(:,j,k),& + w_lorentz(:,j,k), & + 1, j, k, nx, ny, nz, GRHydro_mppm_eigenvalue_x_left, & + GRHydro_mppm_eigenvalue_x_right, & + GRHydro_mppm_xwind) + do i = 1, nx + if (trivial_rp(i,j,k)) then + SpaceMask_SetStateBitsF90(space_mask, i, j, k, type_bitsx, trivialx) + else + SpaceMask_SetStateBitsF90(space_mask, i, j, k, type_bitsx, not_trivialx) + end if + end do + end do + end do + !$OMP END PARALLEL DO + endif !evolve_mhd if(evolve_tracer.ne.0) then !$OMP PARALLEL DO PRIVATE(j) do k = GRHydro_stencil, nz - GRHydro_stencil + 1 @@ -262,31 +357,91 @@ subroutine Reconstruction(CCTK_ARGUMENTS) end if else if (flux_direction == 2) then - !$OMP PARALLEL DO PRIVATE(i, j) - do k = GRHydro_stencil, nz - GRHydro_stencil + 1 - do j = GRHydro_stencil, nx - GRHydro_stencil + 1 - call SimplePPM_1d(GRHydro_eos_handle,0,ny,CCTK_DELTA_SPACE(2),& - rho(j,:,k),vely(j,:,k),velz(j,:,k),velx(j,:,k),eps(j,:,k),& - press(j,:,k),rhominus(j,:,k),velyminus(j,:,k),velzminus(j,:,k),& - velxminus(j,:,k),epsminus(j,:,k),rhoplus(j,:,k),& - velyplus(j,:,k),velzplus(j,:,k),velxplus(j,:,k),epsplus(j,:,k),& - trivial_rp(j,:,k), hydro_excision_mask(j,:,k),& - gyy(j,:,k), gyz(j,:,k), gxy(j,:,k), gzz(j,:,k), gxz(j,:,k), & - gxx(j,:,k), psi4(j,:,k), lbetay(j,:,k), alp(j,:,k),& - w_lorentz(j,:,k), & - 2, j, k, nx, ny, nz, GRHydro_mppm_eigenvalue_y_left, & - GRHydro_mppm_eigenvalue_y_right, & - GRHydro_mppm_xwind) - do i = 1, ny - if (trivial_rp(j,i,k)) then - SpaceMask_SetStateBitsF90(space_mask, j, i, k, type_bitsy, trivialy) - else - SpaceMask_SetStateBitsF90(space_mask, j, i, k, type_bitsy, not_trivialy) - end if + if(evolve_mhd.ne.0) then + if(clean_divergence.ne.0) then + !$OMP PARALLEL DO PRIVATE(i, j) + do k = GRHydro_stencil, nz - GRHydro_stencil + 1 + do j = GRHydro_stencil, ny - GRHydro_stencil + 1 + call SimplePPM_1dM(GRHydro_eos_handle,0,ny,CCTK_DELTA_SPACE(2),& + rho(j,:,k),vely(j,:,k),velz(j,:,k),velx(j,:,k),& + Bvecy(j,:,k),Bvecz(j,:,k),Bvecx(j,:,k),psidc(j,:,k),eps(j,:,k),press(j,:,k),& + rhominus(j,:,k),velyminus(j,:,k),velzminus(j,:,k),velxminus(j,:,k),& + Bvecyminus(j,:,k),Bveczminus(j,:,k),Bvecxminus(j,:,k),psidcminus(j,:,k),epsminus(j,:,k),& + rhoplus(j,:,k),velyplus(j,:,k),velzplus(j,:,k),velxplus(j,:,k),& + Bvecyplus(j,:,k),Bveczplus(j,:,k),Bvecxplus(j,:,k),psidcplus(j,:,k),epsplus(j,:,k),& + 1,trivial_rp(j,:,k), hydro_excision_mask(j,:,k),& + gyy(j,:,k), gyz(j,:,k), gxy(j,:,k), gzz(j,:,k), gxz(j,:,k), & + gxx(j,:,k), psi4(j,:,k), lbetay(j,:,k), alp(j,:,k),& + w_lorentz(j,:,k), & + 2, j, k, nx, ny, nz, GRHydro_mppm_eigenvalue_y_left, & + GRHydro_mppm_eigenvalue_y_right, & + GRHydro_mppm_xwind) + do i = 1, ny + if (trivial_rp(j,i,k)) then + SpaceMask_SetStateBitsF90(space_mask, j, i, k, type_bitsy, trivialy) + else + SpaceMask_SetStateBitsF90(space_mask, j, i, k, type_bitsy, not_trivialy) + end if + end do + end do + end do + !$OMP END PARALLEL DO + else !clean_divergence + !$OMP PARALLEL DO PRIVATE(i, j) + do k = GRHydro_stencil, nz - GRHydro_stencil + 1 + do j = GRHydro_stencil, ny - GRHydro_stencil + 1 + call SimplePPM_1dM(GRHydro_eos_handle,0,ny,CCTK_DELTA_SPACE(2),& + rho(j,:,k),vely(j,:,k),velz(j,:,k),velx(j,:,k),& + Bvecy(j,:,k),Bvecz(j,:,k),Bvecx(j,:,k),dum(j,:,k),eps(j,:,k),press(j,:,k),& + rhominus(j,:,k),velyminus(j,:,k),velzminus(j,:,k),velxminus(j,:,k),& + Bvecyminus(j,:,k),Bveczminus(j,:,k),Bvecxminus(j,:,k),dumm(j,:,k),epsminus(j,:,k),& + rhoplus(j,:,k),velyplus(j,:,k),velzplus(j,:,k),velxplus(j,:,k),& + Bvecyplus(j,:,k),Bveczplus(j,:,k),Bvecxplus(j,:,k),dump(j,:,k),epsplus(j,:,k),& + 0,trivial_rp(j,:,k), hydro_excision_mask(j,:,k),& + gyy(j,:,k), gyz(j,:,k), gxy(j,:,k), gzz(j,:,k), gxz(j,:,k), & + gxx(j,:,k), psi4(j,:,k), lbetay(j,:,k), alp(j,:,k),& + w_lorentz(j,:,k), & + 2, j, k, nx, ny, nz, GRHydro_mppm_eigenvalue_y_left, & + GRHydro_mppm_eigenvalue_y_right, & + GRHydro_mppm_xwind) + do i = 1, ny + if (trivial_rp(j,i,k)) then + SpaceMask_SetStateBitsF90(space_mask, j, i, k, type_bitsy, trivialy) + else + SpaceMask_SetStateBitsF90(space_mask, j, i, k, type_bitsy, not_trivialy) + end if + end do + end do + end do + !$OMP END PARALLEL DO + endif !clean_divergence + else !evolve_mhd + !$OMP PARALLEL DO PRIVATE(i, j) + do k = GRHydro_stencil, nz - GRHydro_stencil + 1 + do j = GRHydro_stencil, ny - GRHydro_stencil + 1 + call SimplePPM_1d(GRHydro_eos_handle,0,ny,CCTK_DELTA_SPACE(2),& + rho(j,:,k),vely(j,:,k),velz(j,:,k),velx(j,:,k),eps(j,:,k),& + press(j,:,k),rhominus(j,:,k),velyminus(j,:,k),velzminus(j,:,k),& + velxminus(j,:,k),epsminus(j,:,k),rhoplus(j,:,k),& + velyplus(j,:,k),velzplus(j,:,k),velxplus(j,:,k),epsplus(j,:,k),& + trivial_rp(j,:,k), hydro_excision_mask(j,:,k),& + gyy(j,:,k), gyz(j,:,k), gxy(j,:,k), gzz(j,:,k), gxz(j,:,k), & + gxx(j,:,k), psi4(j,:,k), lbetay(j,:,k), alp(j,:,k),& + w_lorentz(j,:,k), & + 2, j, k, nx, ny, nz, GRHydro_mppm_eigenvalue_y_left, & + GRHydro_mppm_eigenvalue_y_right, & + GRHydro_mppm_xwind) + do i = 1, ny + if (trivial_rp(j,i,k)) then + SpaceMask_SetStateBitsF90(space_mask, j, i, k, type_bitsy, trivialy) + else + SpaceMask_SetStateBitsF90(space_mask, j, i, k, type_bitsy, not_trivialy) + end if + end do end do end do - end do - !$OMP END PARALLEL DO + !$OMP END PARALLEL DO + endif !evolve_mhd if(evolve_tracer.ne.0) then !$OMP PARALLEL DO PRIVATE(j) do k = GRHydro_stencil, nz - GRHydro_stencil + 1 @@ -313,31 +468,91 @@ subroutine Reconstruction(CCTK_ARGUMENTS) end if else if (flux_direction == 3) then - !$OMP PARALLEL DO PRIVATE(i, j) - do k = GRHydro_stencil, ny - GRHydro_stencil + 1 - do j = GRHydro_stencil, nx - GRHydro_stencil + 1 - call SimplePPM_1d(GRHydro_eos_handle,0,nz,CCTK_DELTA_SPACE(3),& - rho(j,k,:),velz(j,k,:),velx(j,k,:),vely(j,k,:),eps(j,k,:),& - press(j,k,:),rhominus(j,k,:),velzminus(j,k,:),velxminus(j,k,:),& - velyminus(j,k,:),epsminus(j,k,:),rhoplus(j,k,:),& - velzplus(j,k,:),velxplus(j,k,:),velyplus(j,k,:),epsplus(j,k,:),& - trivial_rp(j,k,:), hydro_excision_mask(j,k,:),& - gzz(j,k,:), gxz(j,k,:), gyz(j,k,:), gxx(j,k,:), gxy(j,k,:), & - gyy(j,k,:), psi4(j,k,:), lbetaz(j,k,:), alp(j,k,:),& - w_lorentz(j,k,:), & - 3, j, k, nx, ny, nz, GRHydro_mppm_eigenvalue_z_left, & - GRHydro_mppm_eigenvalue_z_right, & - GRHydro_mppm_xwind) - do i = 1, nz - if (trivial_rp(j,k,i)) then - SpaceMask_SetStateBitsF90(space_mask, j, k, i, type_bitsz, trivialz) - else - SpaceMask_SetStateBitsF90(space_mask, j, k, i, type_bitsz, not_trivialz) - end if + if(evolve_mhd.ne.0) then + if(clean_divergence.ne.0) then + !$OMP PARALLEL DO PRIVATE(i, j) + do k = GRHydro_stencil, ny - GRHydro_stencil + 1 + do j = GRHydro_stencil, nx - GRHydro_stencil + 1 + call SimplePPM_1dM(GRHydro_eos_handle,0,nz,CCTK_DELTA_SPACE(3),& + rho(j,k,:),velz(j,k,:),velx(j,k,:),vely(j,k,:),& + Bvecz(j,k,:),Bvecx(j,k,:),Bvecy(j,k,:),psidc(j,k,:),eps(j,k,:),press(j,k,:),& + rhominus(j,k,:),velzminus(j,k,:),velxminus(j,k,:),velyminus(j,k,:),& + Bveczminus(j,k,:),Bvecxminus(j,k,:),Bvecyminus(j,k,:),psidcminus(j,k,:),epsminus(j,k,:),& + rhoplus(j,k,:),velzplus(j,k,:),velxplus(j,k,:),velyplus(j,k,:),& + Bveczplus(j,k,:),Bvecxplus(j,k,:),Bvecyplus(j,k,:),psidcplus(j,k,:),epsplus(j,k,:),& + 1,trivial_rp(j,k,:), hydro_excision_mask(j,k,:),& + gzz(j,k,:), gxz(j,k,:), gyz(j,k,:), gxx(j,k,:), gxy(j,k,:), & + gyy(j,k,:), psi4(j,k,:), lbetaz(j,k,:), alp(j,k,:),& + w_lorentz(j,k,:), & + 3, j, k, nx, ny, nz, GRHydro_mppm_eigenvalue_z_left, & + GRHydro_mppm_eigenvalue_z_right, & + GRHydro_mppm_xwind) + do i = 1, nz + if (trivial_rp(j,k,i)) then + SpaceMask_SetStateBitsF90(space_mask, j, k, i, type_bitsz, trivialz) + else + SpaceMask_SetStateBitsF90(space_mask, j, k, i, type_bitsz, not_trivialz) + end if + end do + end do + end do + !$OMP END PARALLEL DO + else !clean_divergence + !$OMP PARALLEL DO PRIVATE(i, j) + do k = GRHydro_stencil, ny - GRHydro_stencil + 1 + do j = GRHydro_stencil, nx - GRHydro_stencil + 1 + call SimplePPM_1dM(GRHydro_eos_handle,0,nz,CCTK_DELTA_SPACE(3),& + rho(j,k,:),velz(j,k,:),velx(j,k,:),vely(j,k,:),& + Bvecz(j,k,:),Bvecx(j,k,:),Bvecy(j,k,:),dum(j,k,:),eps(j,k,:),press(j,k,:),& + rhominus(j,k,:),velzminus(j,k,:),velxminus(j,k,:),velyminus(j,k,:),& + Bveczminus(j,k,:),Bvecxminus(j,k,:),Bvecyminus(j,k,:),dumm(j,k,:),epsminus(j,k,:),& + rhoplus(j,k,:),velzplus(j,k,:),velxplus(j,k,:),velyplus(j,k,:),& + Bveczplus(j,k,:),Bvecxplus(j,k,:),Bvecyplus(j,k,:),dump(j,k,:),epsplus(j,k,:),& + 0,trivial_rp(j,k,:), hydro_excision_mask(j,k,:),& + gzz(j,k,:), gxz(j,k,:), gyz(j,k,:), gxx(j,k,:), gxy(j,k,:), & + gyy(j,k,:), psi4(j,k,:), lbetaz(j,k,:), alp(j,k,:),& + w_lorentz(j,k,:), & + 3, j, k, nx, ny, nz, GRHydro_mppm_eigenvalue_z_left, & + GRHydro_mppm_eigenvalue_z_right, & + GRHydro_mppm_xwind) + do i = 1, nz + if (trivial_rp(j,k,i)) then + SpaceMask_SetStateBitsF90(space_mask, j, k, i, type_bitsz, trivialz) + else + SpaceMask_SetStateBitsF90(space_mask, j, k, i, type_bitsz, not_trivialz) + end if + end do + end do + end do + !$OMP END PARALLEL DO + endif !clean_divergence + else !evolve_mhd + !$OMP PARALLEL DO PRIVATE(i, j) + do k = GRHydro_stencil, ny - GRHydro_stencil + 1 + do j = GRHydro_stencil, nx - GRHydro_stencil + 1 + call SimplePPM_1d(GRHydro_eos_handle,0,nz,CCTK_DELTA_SPACE(3),& + rho(j,k,:),velz(j,k,:),velx(j,k,:),vely(j,k,:),eps(j,k,:),& + press(j,k,:),rhominus(j,k,:),velzminus(j,k,:),velxminus(j,k,:),& + velyminus(j,k,:),epsminus(j,k,:),rhoplus(j,k,:),& + velzplus(j,k,:),velxplus(j,k,:),velyplus(j,k,:),epsplus(j,k,:),& + trivial_rp(j,k,:), hydro_excision_mask(j,k,:),& + gzz(j,k,:), gxz(j,k,:), gyz(j,k,:), gxx(j,k,:), gxy(j,k,:), & + gyy(j,k,:), psi4(j,k,:), lbetaz(j,k,:), alp(j,k,:),& + w_lorentz(j,k,:), & + 3, j, k, nx, ny, nz, GRHydro_mppm_eigenvalue_z_left, & + GRHydro_mppm_eigenvalue_z_right, & + GRHydro_mppm_xwind) + do i = 1, nz + if (trivial_rp(j,k,i)) then + SpaceMask_SetStateBitsF90(space_mask, j, k, i, type_bitsz, trivialz) + else + SpaceMask_SetStateBitsF90(space_mask, j, k, i, type_bitsz, not_trivialz) + end if + end do end do end do - end do - !$OMP END PARALLEL DO + !$OMP END PARALLEL DO + endif !evolve_mhd if(evolve_tracer.ne.0) then !$OMP PARALLEL DO PRIVATE(j) do k = GRHydro_stencil, ny - GRHydro_stencil + 1 @@ -426,6 +641,25 @@ subroutine Reconstruction(CCTK_ARGUMENTS) call CCTK_WARN(0, "Variable type to reconstruct not recognized.") !$OMP END CRITICAL end if + +!!$ B-fields are both prim and con + if(evolve_mhd.ne.0) then + call GRHydro_ENOReconstruct1d(eno_order,cctk_lsh(1),& + Bvecx(:,j,k),Bvecxminus(:,j,k),Bvecxplus(:,j,k),& + trivial_rp(:,j,k), hydro_excision_mask(:,j,k)) + call GRHydro_ENOReconstruct1d(eno_order,cctk_lsh(1),& + Bvecy(:,j,k),Bvecyminus(:,j,k),Bvecyplus(:,j,k),& + trivial_rp(:,j,k), hydro_excision_mask(:,j,k)) + call GRHydro_ENOReconstruct1d(eno_order,cctk_lsh(1),& + Bvecz(:,j,k),Bveczminus(:,j,k),Bveczplus(:,j,k),& + trivial_rp(:,j,k), hydro_excision_mask(:,j,k)) + if(clean_divergence.ne.0) then + call GRHydro_ENOReconstruct1d(eno_order,cctk_lsh(1),& + psidc(:,j,k),psidcminus(:,j,k),psidcplus(:,j,k),& + trivial_rp(:,j,k), hydro_excision_mask(:,j,k)) + endif + endif + do i = 1, cctk_lsh(1) if (trivial_rp(i,j,k)) then SpaceMask_SetStateBitsF90(space_mask, i, j, k, type_bitsx, trivialx) @@ -477,6 +711,24 @@ subroutine Reconstruction(CCTK_ARGUMENTS) call CCTK_WARN(0, "Variable type to reconstruct not recognized.") !$OMP END CRITICAL end if + + if(evolve_mhd.ne.0) then + call GRHydro_ENOReconstruct1d(eno_order,cctk_lsh(2),& + Bvecx(j,:,k),Bvecxminus(j,:,k),Bvecxplus(j,:,k),& + trivial_rp(j,:,k), hydro_excision_mask(j,:,k)) + call GRHydro_ENOReconstruct1d(eno_order,cctk_lsh(2),& + Bvecy(j,:,k),Bvecyminus(j,:,k),Bvecyplus(j,:,k),& + trivial_rp(j,:,k), hydro_excision_mask(j,:,k)) + call GRHydro_ENOReconstruct1d(eno_order,cctk_lsh(2),& + Bvecz(j,:,k),Bveczminus(j,:,k),Bveczplus(j,:,k),& + trivial_rp(j,:,k), hydro_excision_mask(j,:,k)) + if(clean_divergence.ne.0) then + call GRHydro_ENOReconstruct1d(eno_order,cctk_lsh(2),& + psidc(j,:,k),psidcminus(j,:,k),psidcplus(j,:,k),& + trivial_rp(j,:,k), hydro_excision_mask(j,:,k)) + endif + endif + do i = 1, cctk_lsh(2) if (trivial_rp(j,i,k)) then SpaceMask_SetStateBitsF90(space_mask, j, i, k, type_bitsy, trivialy) @@ -528,6 +780,24 @@ subroutine Reconstruction(CCTK_ARGUMENTS) call CCTK_WARN(0, "Variable type to reconstruct not recognized.") !$OMP END CRITICAL end if + + if(evolve_mhd.ne.0) then + call GRHydro_ENOReconstruct1d(eno_order,cctk_lsh(3),& + Bvecx(j,k,:),Bvecxminus(j,k,:),Bvecxplus(j,k,:),& + trivial_rp(j,k,:), hydro_excision_mask(j,k,:)) + call GRHydro_ENOReconstruct1d(eno_order,cctk_lsh(3),& + Bvecy(j,k,:),Bvecyminus(j,k,:),Bvecyplus(j,k,:),& + trivial_rp(j,k,:), hydro_excision_mask(j,k,:)) + call GRHydro_ENOReconstruct1d(eno_order,cctk_lsh(3),& + Bvecz(j,k,:),Bveczminus(j,k,:),Bveczplus(j,k,:),& + trivial_rp(j,k,:), hydro_excision_mask(j,k,:)) + if(clean_divergence.ne.0) then + call GRHydro_ENOReconstruct1d(eno_order,cctk_lsh(3),& + psidc(j,k,:),psidcminus(j,k,:),psidcplus(j,k,:),& + trivial_rp(j,k,:), hydro_excision_mask(j,k,:)) + endif + endif + do i = 1, cctk_lsh(3) if (trivial_rp(j,k,i)) then SpaceMask_SetStateBitsF90(space_mask, j, k, i, type_bitsz, trivialz) @@ -547,6 +817,7 @@ subroutine Reconstruction(CCTK_ARGUMENTS) deallocate(trivial_rp) deallocate(psi4, lbetax, lbetay, lbetaz) + deallocate(dum,dump,dumm) !$OMP WORKSHARE where ( (rhoplus < GRHydro_rho_min).or.(rhominus < GRHydro_rho_min) ) @@ -585,12 +856,24 @@ subroutine Reconstruction(CCTK_ARGUMENTS) if (CCTK_EQUALS(recon_vars,"primitive").or.& CCTK_EQUALS(recon_method,"ppm")) then if (use_eosgeneral == 0) then - call primitive2conservative(CCTK_PASS_FTOF) + if(evolve_mhd.ne.0) then + call primitive2conservativeM(CCTK_PASS_FTOF) + else + call primitive2conservative(CCTK_PASS_FTOF) + endif else - call primitive2conservativegeneral(CCTK_PASS_FTOF) + if(evolve_mhd.ne.0) then + call primitive2conservativegeneralM(CCTK_PASS_FTOF) + else + call primitive2conservativegeneral(CCTK_PASS_FTOF) + endif end if else if (CCTK_EQUALS(recon_vars,"conservative")) then - call Conservative2PrimitiveBounds(CCTK_PASS_FTOF) + if(evolve_mhd.ne.0) then + call Conservative2PrimitiveBoundsM(CCTK_PASS_FTOF) + else + call Conservative2PrimitiveBounds(CCTK_PASS_FTOF) + endif else call CCTK_WARN(0,"Variable type to reconstruct not recognized.") end if diff --git a/src/GRHydro_ReconstructPoly.F90 b/src/GRHydro_ReconstructPoly.F90 index 06c969b..cd735e1 100644 --- a/src/GRHydro_ReconstructPoly.F90 +++ b/src/GRHydro_ReconstructPoly.F90 @@ -20,6 +20,9 @@ #define sx(i,j,k) scon(i,j,k,1) #define sy(i,j,k) scon(i,j,k,2) #define sz(i,j,k) scon(i,j,k,3) +#define Bvecx(i,j,k) Bvec(i,j,k,1) +#define Bvecy(i,j,k) Bvec(i,j,k,2) +#define Bvecz(i,j,k) Bvec(i,j,k,3) /*@@ @@ -56,7 +59,7 @@ subroutine ReconstructionPolytype(CCTK_ARGUMENTS) &type_bitsz, trivialz, not_trivialz CCTK_REAL, dimension(:,:,:),allocatable :: & - &psi4, lbetax, lbetay, lbetaz + &psi4, lbetax, lbetay, lbetaz, dum, dump, dumm !!$ CCTK_REAL, dimension(cctk_lsh(1),cctk_lsh(2),cctk_lsh(3)) :: & !!$ &psi4, lbetax, lbetay, lbetaz @@ -72,7 +75,11 @@ subroutine ReconstructionPolytype(CCTK_ARGUMENTS) allocate(psi4(cctk_lsh(1),cctk_lsh(2),cctk_lsh(3)),& lbetax(cctk_lsh(1),cctk_lsh(2),cctk_lsh(3)),& lbetay(cctk_lsh(1),cctk_lsh(2),cctk_lsh(3)),& - lbetaz(cctk_lsh(1),cctk_lsh(2),cctk_lsh(3)),STAT=ierr) + lbetaz(cctk_lsh(1),cctk_lsh(2),cctk_lsh(3)),& + dum(cctk_lsh(1),cctk_lsh(2),cctk_lsh(3)),& + dump(cctk_lsh(1),cctk_lsh(2),cctk_lsh(3)),& + dumm(cctk_lsh(1),cctk_lsh(2),cctk_lsh(3)),STAT=ierr) + if (ierr .ne. 0) then call CCTK_WARN(0, "Allocation problems with lbeta") end if @@ -125,6 +132,18 @@ subroutine ReconstructionPolytype(CCTK_ARGUMENTS) velyminus = 0.0d0 velzplus = 0.0d0 velzminus = 0.0d0 + if(evolve_mhd.ne.0) then + Bvecxplus = 0.0d0 + Bvecxminus = 0.0d0 + Bvecyplus = 0.0d0 + Bvecyminus = 0.0d0 + Bveczplus = 0.0d0 + Bveczminus = 0.0d0 + if(clean_divergence.ne.0) then + psidcplus = 0.0d0 + psidcminus=0.0d0 + endif + endif if (evolve_tracer .ne. 0) then tracerplus = 0.0d0 @@ -178,6 +197,20 @@ subroutine ReconstructionPolytype(CCTK_ARGUMENTS) call CCTK_WARN(0, "Variable type to reconstruct not recognized.") end if +!!$ B-field always needs reconstruction, as it is both prim and con + if(evolve_mhd.ne.0) then + call tvdreconstruct(nx, ny, nz, xoffset, yoffset, zoffset, & + Bvec(:,:,:,1), Bvecxplus, Bvecxminus, trivial_rp, hydro_excision_mask) + call tvdreconstruct(nx, ny, nz, xoffset, yoffset, zoffset, & + Bvec(:,:,:,2), Bvecyplus, Bvecyminus, trivial_rp, hydro_excision_mask) + call tvdreconstruct(nx, ny, nz, xoffset, yoffset, zoffset, & + Bvec(:,:,:,3), Bveczplus, Bveczminus, trivial_rp, hydro_excision_mask) + if(clean_divergence.ne.0) then + call tvdreconstruct(nx, ny, nz, xoffset, yoffset, zoffset, & + psidc, psidcplus, psidcminus, trivial_rp, hydro_excision_mask) + endif + endif + !$OMP PARALLEL DO PRIVATE(i, j) do k = 1, nz do j = 1, ny @@ -207,31 +240,91 @@ subroutine ReconstructionPolytype(CCTK_ARGUMENTS) else if (CCTK_EQUALS(recon_method,"ppm")) then if (flux_direction == 1) then - !$OMP PARALLEL DO PRIVATE(i, j) - do k = GRHydro_stencil, nz - GRHydro_stencil + 1 - do j = GRHydro_stencil, ny - GRHydro_stencil + 1 - call SimplePPM_1d(GRHydro_eos_handle,1,nx,CCTK_DELTA_SPACE(1),& - rho(:,j,k),velx(:,j,k),vely(:,j,k),velz(:,j,k),eps(:,j,k),& - press(:,j,k),rhominus(:,j,k),velxminus(:,j,k),velyminus(:,j,k),& - velzminus(:,j,k),epsminus(:,j,k),rhoplus(:,j,k),& - velxplus(:,j,k),velyplus(:,j,k),velzplus(:,j,k),epsplus(:,j,k),& - trivial_rp(:,j,k), hydro_excision_mask(:,j,k),& - gxx(:,j,k), gxy(:,j,k), gxz(:,j,k), gyy(:,j,k), gyz(:,j,k), & - gzz(:,j,k), psi4(:,j,k), lbetax(:,j,k), alp(:,j,k),& - w_lorentz(:,j,k), & - 1, j, k, nx, ny, nz, GRHydro_mppm_eigenvalue_x_left, & - GRHydro_mppm_eigenvalue_x_right, & - GRHydro_mppm_xwind) - do i = 1, nx - if (trivial_rp(i,j,k)) then - SpaceMask_SetStateBitsF90(space_mask, i, j, k, type_bitsx, trivialx) - else - SpaceMask_SetStateBitsF90(space_mask, i, j, k, type_bitsx, not_trivialx) - end if + if(evolve_mhd.ne.0) then + if(clean_divergence.ne.0) then + !$OMP PARALLEL DO PRIVATE(i, j) + do k = GRHydro_stencil, nz - GRHydro_stencil + 1 + do j = GRHydro_stencil, ny - GRHydro_stencil + 1 + call SimplePPM_1dM(GRHydro_eos_handle,1,nx,CCTK_DELTA_SPACE(1),& + rho(:,j,k),velx(:,j,k),vely(:,j,k),velz(:,j,k),& + Bvecx(:,j,k),Bvecy(:,j,k),Bvecz(:,j,k),psidc(:,j,k),eps(:,j,k),press(:,j,k),& + rhominus(:,j,k),velxminus(:,j,k),velyminus(:,j,k),velzminus(:,j,k),& + Bvecxminus(:,j,k),Bvecyminus(:,j,k),Bveczminus(:,j,k),psidcminus(:,j,k),epsminus(:,j,k),& + rhoplus(:,j,k),velxplus(:,j,k),velyplus(:,j,k),velzplus(:,j,k),& + Bvecxplus(:,j,k),Bvecyplus(:,j,k),Bveczplus(:,j,k),psidcplus(:,j,k),epsplus(:,j,k),& + 1,trivial_rp(:,j,k), hydro_excision_mask(:,j,k),& + gxx(:,j,k), gxy(:,j,k), gxz(:,j,k), gyy(:,j,k), gyz(:,j,k), & + gzz(:,j,k), psi4(:,j,k), lbetax(:,j,k), alp(:,j,k),& + w_lorentz(:,j,k), & + 1, j, k, nx, ny, nz, GRHydro_mppm_eigenvalue_x_left, & + GRHydro_mppm_eigenvalue_x_right, & + GRHydro_mppm_xwind) + do i = 1, nx + if (trivial_rp(i,j,k)) then + SpaceMask_SetStateBitsF90(space_mask, i, j, k, type_bitsx, trivialx) + else + SpaceMask_SetStateBitsF90(space_mask, i, j, k, type_bitsx, not_trivialx) + end if + end do + end do + end do + !$OMP END PARALLEL DO + else !clean_divergence + !$OMP PARALLEL DO PRIVATE(i, j) + do k = GRHydro_stencil, nz - GRHydro_stencil + 1 + do j = GRHydro_stencil, ny - GRHydro_stencil + 1 + call SimplePPM_1dM(GRHydro_eos_handle,1,nx,CCTK_DELTA_SPACE(1),& + rho(:,j,k),velx(:,j,k),vely(:,j,k),velz(:,j,k),& + Bvecx(:,j,k),Bvecy(:,j,k),Bvecz(:,j,k),dum(:,j,k),eps(:,j,k),press(:,j,k),& + rhominus(:,j,k),velxminus(:,j,k),velyminus(:,j,k),velzminus(:,j,k),& + Bvecxminus(:,j,k),Bvecyminus(:,j,k),Bveczminus(:,j,k),dumm(:,j,k),epsminus(:,j,k),& + rhoplus(:,j,k),velxplus(:,j,k),velyplus(:,j,k),velzplus(:,j,k),& + Bvecxplus(:,j,k),Bvecyplus(:,j,k),Bveczplus(:,j,k),dump(:,j,k),epsplus(:,j,k),& + 0,trivial_rp(:,j,k), hydro_excision_mask(:,j,k),& + gxx(:,j,k), gxy(:,j,k), gxz(:,j,k), gyy(:,j,k), gyz(:,j,k), & + gzz(:,j,k), psi4(:,j,k), lbetax(:,j,k), alp(:,j,k),& + w_lorentz(:,j,k), & + 1, j, k, nx, ny, nz, GRHydro_mppm_eigenvalue_x_left, & + GRHydro_mppm_eigenvalue_x_right, & + GRHydro_mppm_xwind) + do i = 1, nx + if (trivial_rp(i,j,k)) then + SpaceMask_SetStateBitsF90(space_mask, i, j, k, type_bitsx, trivialx) + else + SpaceMask_SetStateBitsF90(space_mask, i, j, k, type_bitsx, not_trivialx) + end if + end do + end do + end do + !$OMP END PARALLEL DO + endif !clean_divergence + else !evolve_mhd + !$OMP PARALLEL DO PRIVATE(i, j) + do k = GRHydro_stencil, nz - GRHydro_stencil + 1 + do j = GRHydro_stencil, ny - GRHydro_stencil + 1 + call SimplePPM_1d(GRHydro_eos_handle,1,nx,CCTK_DELTA_SPACE(1),& + rho(:,j,k),velx(:,j,k),vely(:,j,k),velz(:,j,k),eps(:,j,k),& + press(:,j,k),rhominus(:,j,k),velxminus(:,j,k),velyminus(:,j,k),& + velzminus(:,j,k),epsminus(:,j,k),rhoplus(:,j,k),& + velxplus(:,j,k),velyplus(:,j,k),velzplus(:,j,k),epsplus(:,j,k),& + trivial_rp(:,j,k), hydro_excision_mask(:,j,k),& + gxx(:,j,k), gxy(:,j,k), gxz(:,j,k), gyy(:,j,k), gyz(:,j,k), & + gzz(:,j,k), psi4(:,j,k), lbetax(:,j,k), alp(:,j,k),& + w_lorentz(:,j,k), & + 1, j, k, nx, ny, nz, GRHydro_mppm_eigenvalue_x_left, & + GRHydro_mppm_eigenvalue_x_right, & + GRHydro_mppm_xwind) + do i = 1, nx + if (trivial_rp(i,j,k)) then + SpaceMask_SetStateBitsF90(space_mask, i, j, k, type_bitsx, trivialx) + else + SpaceMask_SetStateBitsF90(space_mask, i, j, k, type_bitsx, not_trivialx) + end if + end do end do end do - end do - !$OMP END PARALLEL DO + !$OMP END PARALLEL DO + endif !evolve_mhd if(evolve_tracer.ne.0) then !$OMP PARALLEL DO PRIVATE(j) do k = GRHydro_stencil, nz - GRHydro_stencil + 1 @@ -258,31 +351,91 @@ subroutine ReconstructionPolytype(CCTK_ARGUMENTS) end if else if (flux_direction == 2) then - !$OMP PARALLEL DO PRIVATE(i, j) - do k = GRHydro_stencil, nz - GRHydro_stencil + 1 - do j = GRHydro_stencil, nx - GRHydro_stencil + 1 - call SimplePPM_1d(GRHydro_eos_handle,1,ny,CCTK_DELTA_SPACE(2),& - rho(j,:,k),vely(j,:,k),velz(j,:,k),velx(j,:,k),eps(j,:,k),& - press(j,:,k),rhominus(j,:,k),velyminus(j,:,k),velzminus(j,:,k),& - velxminus(j,:,k),epsminus(j,:,k),rhoplus(j,:,k),& - velyplus(j,:,k),velzplus(j,:,k),velxplus(j,:,k),epsplus(j,:,k),& - trivial_rp(j,:,k), hydro_excision_mask(j,:,k),& - gyy(j,:,k), gyz(j,:,k), gxy(j,:,k), gzz(j,:,k), gxz(j,:,k), & - gxx(j,:,k), psi4(j,:,k), lbetay(j,:,k), alp(j,:,k),& - w_lorentz(j,:,k), & - 2, j, k, nx, ny, nz, GRHydro_mppm_eigenvalue_y_left, & - GRHydro_mppm_eigenvalue_y_right, & - GRHydro_mppm_xwind) - do i = 1, ny - if (trivial_rp(j,i,k)) then - SpaceMask_SetStateBitsF90(space_mask, j, i, k, type_bitsy, trivialy) - else - SpaceMask_SetStateBitsF90(space_mask, j, i, k, type_bitsy, not_trivialy) - end if + if(evolve_mhd.ne.0) then + if(clean_divergence.ne.0) then + !$OMP PARALLEL DO PRIVATE(i, j) + do k = GRHydro_stencil, nz - GRHydro_stencil + 1 + do j = GRHydro_stencil, nx - GRHydro_stencil + 1 + call SimplePPM_1dM(GRHydro_eos_handle,1,ny,CCTK_DELTA_SPACE(2),& + rho(j,:,k),vely(j,:,k),velz(j,:,k),velx(j,:,k),& + Bvecy(j,:,k),Bvecz(j,:,k),Bvecx(j,:,k),psidc(j,:,k),eps(j,:,k),press(j,:,k),& + rhominus(j,:,k),velyminus(j,:,k),velzminus(j,:,k),velxminus(j,:,k),& + Bvecyminus(j,:,k),Bveczminus(j,:,k),Bvecxminus(j,:,k),psidcminus(j,:,k),epsminus(j,:,k),& + rhoplus(j,:,k),velyplus(j,:,k),velzplus(j,:,k),velxplus(j,:,k),& + Bvecyplus(j,:,k),Bveczplus(j,:,k),Bvecxplus(j,:,k),psidcplus(j,:,k),epsplus(j,:,k),& + 1,trivial_rp(j,:,k), hydro_excision_mask(j,:,k),& + gyy(j,:,k), gyz(j,:,k), gxy(j,:,k), gzz(j,:,k), gxz(j,:,k), & + gxx(j,:,k), psi4(j,:,k), lbetay(j,:,k), alp(j,:,k),& + w_lorentz(j,:,k), & + 2, j, k, nx, ny, nz, GRHydro_mppm_eigenvalue_y_left, & + GRHydro_mppm_eigenvalue_y_right, & + GRHydro_mppm_xwind) + do i = 1, ny + if (trivial_rp(j,i,k)) then + SpaceMask_SetStateBitsF90(space_mask, j, i, k, type_bitsy, trivialy) + else + SpaceMask_SetStateBitsF90(space_mask, j, i, k, type_bitsy, not_trivialy) + end if + end do + end do + end do + !$OMP END PARALLEL DO + else !clean_divergence + !$OMP PARALLEL DO PRIVATE(i, j) + do k = GRHydro_stencil, nz - GRHydro_stencil + 1 + do j = GRHydro_stencil, nx - GRHydro_stencil + 1 + call SimplePPM_1dM(GRHydro_eos_handle,1,ny,CCTK_DELTA_SPACE(2),& + rho(j,:,k),vely(j,:,k),velz(j,:,k),velx(j,:,k),& + Bvecy(j,:,k),Bvecz(j,:,k),Bvecx(j,:,k),dum(j,:,k),eps(j,:,k),press(j,:,k),& + rhominus(j,:,k),velyminus(j,:,k),velzminus(j,:,k),velxminus(j,:,k),& + Bvecyminus(j,:,k),Bveczminus(j,:,k),Bvecxminus(j,:,k),dumm(j,:,k),epsminus(j,:,k),& + rhoplus(j,:,k),velyplus(j,:,k),velzplus(j,:,k),velxplus(j,:,k),& + Bvecyplus(j,:,k),Bveczplus(j,:,k),Bvecxplus(j,:,k),dump(j,:,k),epsplus(j,:,k),& + 0,trivial_rp(j,:,k), hydro_excision_mask(j,:,k),& + gyy(j,:,k), gyz(j,:,k), gxy(j,:,k), gzz(j,:,k), gxz(j,:,k), & + gxx(j,:,k), psi4(j,:,k), lbetay(j,:,k), alp(j,:,k),& + w_lorentz(j,:,k), & + 2, j, k, nx, ny, nz, GRHydro_mppm_eigenvalue_y_left, & + GRHydro_mppm_eigenvalue_y_right, & + GRHydro_mppm_xwind) + do i = 1, ny + if (trivial_rp(j,i,k)) then + SpaceMask_SetStateBitsF90(space_mask, j, i, k, type_bitsy, trivialy) + else + SpaceMask_SetStateBitsF90(space_mask, j, i, k, type_bitsy, not_trivialy) + end if + end do + end do + end do + !$OMP END PARALLEL DO + endif !clean_divergence + else !evolve_mhd + !$OMP PARALLEL DO PRIVATE(i, j) + do k = GRHydro_stencil, nz - GRHydro_stencil + 1 + do j = GRHydro_stencil, nx - GRHydro_stencil + 1 + call SimplePPM_1d(GRHydro_eos_handle,1,ny,CCTK_DELTA_SPACE(2),& + rho(j,:,k),vely(j,:,k),velz(j,:,k),velx(j,:,k),eps(j,:,k),& + press(j,:,k),rhominus(j,:,k),velyminus(j,:,k),velzminus(j,:,k),& + velxminus(j,:,k),epsminus(j,:,k),rhoplus(j,:,k),& + velyplus(j,:,k),velzplus(j,:,k),velxplus(j,:,k),epsplus(j,:,k),& + trivial_rp(j,:,k), hydro_excision_mask(j,:,k),& + gyy(j,:,k), gyz(j,:,k), gxy(j,:,k), gzz(j,:,k), gxz(j,:,k), & + gxx(j,:,k), psi4(j,:,k), lbetay(j,:,k), alp(j,:,k),& + w_lorentz(j,:,k), & + 2, j, k, nx, ny, nz, GRHydro_mppm_eigenvalue_y_left, & + GRHydro_mppm_eigenvalue_y_right, & + GRHydro_mppm_xwind) + do i = 1, ny + if (trivial_rp(j,i,k)) then + SpaceMask_SetStateBitsF90(space_mask, j, i, k, type_bitsy, trivialy) + else + SpaceMask_SetStateBitsF90(space_mask, j, i, k, type_bitsy, not_trivialy) + end if + end do end do end do - end do - !$OMP END PARALLEL DO + !$OMP END PARALLEL DO + endif !evolve_mhd if(evolve_tracer.ne.0) then !$OMP PARALLEL DO PRIVATE(j) do k = GRHydro_stencil, nz - GRHydro_stencil + 1 @@ -309,31 +462,91 @@ subroutine ReconstructionPolytype(CCTK_ARGUMENTS) end if else if (flux_direction == 3) then - !$OMP PARALLEL DO PRIVATE(i, j) - do k = GRHydro_stencil, ny - GRHydro_stencil + 1 - do j = GRHydro_stencil, nx - GRHydro_stencil + 1 - call SimplePPM_1d(GRHydro_eos_handle,1,nz,CCTK_DELTA_SPACE(3),& - rho(j,k,:),velz(j,k,:),velx(j,k,:),vely(j,k,:),eps(j,k,:),& - press(j,k,:),rhominus(j,k,:),velzminus(j,k,:),velxminus(j,k,:),& - velyminus(j,k,:),epsminus(j,k,:),rhoplus(j,k,:),& - velzplus(j,k,:),velxplus(j,k,:),velyplus(j,k,:),epsplus(j,k,:),& - trivial_rp(j,k,:), hydro_excision_mask(j,k,:),& - gzz(j,k,:), gxz(j,k,:), gyz(j,k,:), gxx(j,k,:), gxy(j,k,:), & - gyy(j,k,:), psi4(j,k,:), lbetaz(j,k,:), alp(j,k,:),& - w_lorentz(j,k,:), & - 3, j, k, nx, ny, nz, GRHydro_mppm_eigenvalue_z_left, & - GRHydro_mppm_eigenvalue_z_right, & - GRHydro_mppm_xwind) - do i = 1, nz - if (trivial_rp(j,k,i)) then - SpaceMask_SetStateBitsF90(space_mask, j, k, i, type_bitsz, trivialz) - else - SpaceMask_SetStateBitsF90(space_mask, j, k, i, type_bitsz, not_trivialz) - end if + if(evolve_mhd.ne.0) then + if(clean_divergence.ne.0) then + !$OMP PARALLEL DO PRIVATE(i, j) + do k = GRHydro_stencil, ny - GRHydro_stencil + 1 + do j = GRHydro_stencil, nx - GRHydro_stencil + 1 + call SimplePPM_1dM(GRHydro_eos_handle,1,nz,CCTK_DELTA_SPACE(3),& + rho(j,k,:),velz(j,k,:),velx(j,k,:),vely(j,k,:),& + Bvecz(j,k,:),Bvecx(j,k,:),Bvecy(j,k,:),psidc(:,j,k),eps(j,k,:),press(j,k,:),& + rhominus(j,k,:),velzminus(j,k,:),velxminus(j,k,:),velyminus(j,k,:),& + Bveczminus(j,k,:),Bvecxminus(j,k,:),Bvecyminus(j,k,:),psidcminus(:,j,k),epsminus(j,k,:),& + rhoplus(j,k,:),velzplus(j,k,:),velxplus(j,k,:),velyplus(j,k,:),& + Bveczplus(j,k,:),Bvecxplus(j,k,:),Bvecyplus(j,k,:),psidcplus(:,j,k),epsplus(j,k,:),& + 1,trivial_rp(j,k,:), hydro_excision_mask(j,k,:),& + gzz(j,k,:), gxz(j,k,:), gyz(j,k,:), gxx(j,k,:), gxy(j,k,:), & + gyy(j,k,:), psi4(j,k,:), lbetaz(j,k,:), alp(j,k,:),& + w_lorentz(j,k,:), & + 3, j, k, nx, ny, nz, GRHydro_mppm_eigenvalue_z_left, & + GRHydro_mppm_eigenvalue_z_right, & + GRHydro_mppm_xwind) + do i = 1, nz + if (trivial_rp(j,k,i)) then + SpaceMask_SetStateBitsF90(space_mask, j, k, i, type_bitsz, trivialz) + else + SpaceMask_SetStateBitsF90(space_mask, j, k, i, type_bitsz, not_trivialz) + end if + end do + end do + end do + !$OMP END PARALLEL DO + else !clean_divergence + !$OMP PARALLEL DO PRIVATE(i, j) + do k = GRHydro_stencil, ny - GRHydro_stencil + 1 + do j = GRHydro_stencil, nx - GRHydro_stencil + 1 + call SimplePPM_1dM(GRHydro_eos_handle,1,nz,CCTK_DELTA_SPACE(3),& + rho(j,k,:),velz(j,k,:),velx(j,k,:),vely(j,k,:),& + Bvecz(j,k,:),Bvecx(j,k,:),Bvecy(j,k,:),dum(:,j,k),eps(j,k,:),press(j,k,:),& + rhominus(j,k,:),velzminus(j,k,:),velxminus(j,k,:),velyminus(j,k,:),& + Bveczminus(j,k,:),Bvecxminus(j,k,:),Bvecyminus(j,k,:),dumm(:,j,k),epsminus(j,k,:),& + rhoplus(j,k,:),velzplus(j,k,:),velxplus(j,k,:),velyplus(j,k,:),& + Bveczplus(j,k,:),Bvecxplus(j,k,:),Bvecyplus(j,k,:),dump(:,j,k),epsplus(j,k,:),& + 0,trivial_rp(j,k,:), hydro_excision_mask(j,k,:),& + gzz(j,k,:), gxz(j,k,:), gyz(j,k,:), gxx(j,k,:), gxy(j,k,:), & + gyy(j,k,:), psi4(j,k,:), lbetaz(j,k,:), alp(j,k,:),& + w_lorentz(j,k,:), & + 3, j, k, nx, ny, nz, GRHydro_mppm_eigenvalue_z_left, & + GRHydro_mppm_eigenvalue_z_right, & + GRHydro_mppm_xwind) + do i = 1, nz + if (trivial_rp(j,k,i)) then + SpaceMask_SetStateBitsF90(space_mask, j, k, i, type_bitsz, trivialz) + else + SpaceMask_SetStateBitsF90(space_mask, j, k, i, type_bitsz, not_trivialz) + end if + end do + end do + end do + !$OMP END PARALLEL DO + endif !clean_divergence + else !evolve_mhd + !$OMP PARALLEL DO PRIVATE(i, j) + do k = GRHydro_stencil, ny - GRHydro_stencil + 1 + do j = GRHydro_stencil, nx - GRHydro_stencil + 1 + call SimplePPM_1d(GRHydro_eos_handle,1,nz,CCTK_DELTA_SPACE(3),& + rho(j,k,:),velz(j,k,:),velx(j,k,:),vely(j,k,:),eps(j,k,:),& + press(j,k,:),rhominus(j,k,:),velzminus(j,k,:),velxminus(j,k,:),& + velyminus(j,k,:),epsminus(j,k,:),rhoplus(j,k,:),& + velzplus(j,k,:),velxplus(j,k,:),velyplus(j,k,:),epsplus(j,k,:),& + trivial_rp(j,k,:), hydro_excision_mask(j,k,:),& + gzz(j,k,:), gxz(j,k,:), gyz(j,k,:), gxx(j,k,:), gxy(j,k,:), & + gyy(j,k,:), psi4(j,k,:), lbetaz(j,k,:), alp(j,k,:),& + w_lorentz(j,k,:), & + 3, j, k, nx, ny, nz, GRHydro_mppm_eigenvalue_z_left, & + GRHydro_mppm_eigenvalue_z_right, & + GRHydro_mppm_xwind) + do i = 1, nz + if (trivial_rp(j,k,i)) then + SpaceMask_SetStateBitsF90(space_mask, j, k, i, type_bitsz, trivialz) + else + SpaceMask_SetStateBitsF90(space_mask, j, k, i, type_bitsz, not_trivialz) + end if + end do end do end do - end do - !$OMP END PARALLEL DO + !$OMP END PARALLEL DO + endif !evolve_mhd if(evolve_tracer.ne.0) then !$OMP PARALLEL DO PRIVATE(j) do k = GRHydro_stencil, ny - GRHydro_stencil + 1 @@ -415,6 +628,25 @@ subroutine ReconstructionPolytype(CCTK_ARGUMENTS) call CCTK_WARN(0, "Variable type to reconstruct not recognized.") !$OMP END CRITICAL end if + + if(evolve_mhd.ne.0) then +!!$ Always do B-field + call GRHydro_ENOReconstruct1d(eno_order,cctk_lsh(1),& + Bvecx(:,j,k),Bvecxminus(:,j,k),Bvecxplus(:,j,k),& + trivial_rp(:,j,k), hydro_excision_mask(:,j,k)) + call GRHydro_ENOReconstruct1d(eno_order,cctk_lsh(1),& + Bvecy(:,j,k),Bvecyminus(:,j,k),Bvecyplus(:,j,k),& + trivial_rp(:,j,k), hydro_excision_mask(:,j,k)) + call GRHydro_ENOReconstruct1d(eno_order,cctk_lsh(1),& + Bvecz(:,j,k),Bveczminus(:,j,k),Bveczplus(:,j,k),& + trivial_rp(:,j,k), hydro_excision_mask(:,j,k)) + if(clean_divergence.ne.0) then + call GRHydro_ENOReconstruct1d(eno_order,cctk_lsh(1),& + psidc(:,j,k),psidcminus(:,j,k),psidcplus(:,j,k),& + trivial_rp(:,j,k), hydro_excision_mask(:,j,k)) + endif + endif + do i = 1, cctk_lsh(1) if (trivial_rp(i,j,k)) then SpaceMask_SetStateBitsF90(space_mask, i, j, k, type_bitsx, trivialx) @@ -460,6 +692,24 @@ subroutine ReconstructionPolytype(CCTK_ARGUMENTS) call CCTK_WARN(0, "Variable type to reconstruct not recognized.") !$OMP END CRITICAL end if + + if(evolve_mhd.ne.0) then + call GRHydro_ENOReconstruct1d(eno_order,cctk_lsh(2),& + Bvecx(j,:,k),Bvecxminus(j,:,k),Bvecxplus(j,:,k),& + trivial_rp(j,:,k), hydro_excision_mask(j,:,k)) + call GRHydro_ENOReconstruct1d(eno_order,cctk_lsh(2),& + Bvecy(j,:,k),Bvecyminus(j,:,k),Bvecyplus(j,:,k),& + trivial_rp(j,:,k), hydro_excision_mask(j,:,k)) + call GRHydro_ENOReconstruct1d(eno_order,cctk_lsh(2),& + Bvecz(j,:,k),Bveczminus(j,:,k),Bveczplus(j,:,k),& + trivial_rp(j,:,k), hydro_excision_mask(j,:,k)) + if(clean_divergence.ne.0) then + call GRHydro_ENOReconstruct1d(eno_order,cctk_lsh(2),& + psidc(j,:,k),psidcminus(j,:,k),psidcplus(j,:,k),& + trivial_rp(j,:,k), hydro_excision_mask(j,:,k)) + endif + endif + do i = 1, cctk_lsh(2) if (trivial_rp(j,i,k)) then SpaceMask_SetStateBitsF90(space_mask, j, i, k, type_bitsy, trivialy) @@ -505,6 +755,24 @@ subroutine ReconstructionPolytype(CCTK_ARGUMENTS) call CCTK_WARN(0, "Variable type to reconstruct not recognized.") !$OMP END CRITICAL end if + + if(evolve_mhd.ne.0) then + call GRHydro_ENOReconstruct1d(eno_order,cctk_lsh(3),& + Bvecx(j,k,:),Bvecxminus(j,k,:),Bvecxplus(j,k,:),& + trivial_rp(j,k,:), hydro_excision_mask(j,k,:)) + call GRHydro_ENOReconstruct1d(eno_order,cctk_lsh(3),& + Bvecy(j,k,:),Bvecyminus(j,k,:),Bvecyplus(j,k,:),& + trivial_rp(j,k,:), hydro_excision_mask(j,k,:)) + call GRHydro_ENOReconstruct1d(eno_order,cctk_lsh(3),& + Bvecz(j,k,:),Bveczminus(j,k,:),Bveczplus(j,k,:),& + trivial_rp(j,k,:), hydro_excision_mask(j,k,:)) + if(clean_divergence.ne.0) then + call GRHydro_ENOReconstruct1d(eno_order,cctk_lsh(3),& + psidc(j,k,:),psidcminus(j,k,:),psidcplus(j,k,:),& + trivial_rp(j,k,:), hydro_excision_mask(j,k,:)) + endif + endif + do i = 1, cctk_lsh(3) if (trivial_rp(j,k,i)) then SpaceMask_SetStateBitsF90(space_mask, j, k, i, type_bitsz, trivialz) @@ -524,6 +792,7 @@ subroutine ReconstructionPolytype(CCTK_ARGUMENTS) deallocate(trivial_rp) deallocate(psi4, lbetax, lbetay, lbetaz) + deallocate(dum, dump, dumm) !$OMP WORKSHARE where ( (rhoplus < GRHydro_rho_min).or.(rhominus < GRHydro_rho_min).or.& @@ -563,17 +832,29 @@ subroutine ReconstructionPolytype(CCTK_ARGUMENTS) if (CCTK_EQUALS(recon_vars,"primitive").or.& CCTK_EQUALS(recon_method,"ppm")) then if (use_eosgeneral == 0) then - call Prim2ConservativePolytype(CCTK_PASS_FTOF) + if(evolve_mhd.ne.0) then + call Prim2ConservativePolytypeM(CCTK_PASS_FTOF) + else + call Prim2ConservativePolytype(CCTK_PASS_FTOF) + endif else - call primitive2conservativegeneral(CCTK_PASS_FTOF) + if(evolve_mhd.ne.0) then + call primitive2conservativegeneral(CCTK_PASS_FTOF) + else + call primitive2conservativegeneral(CCTK_PASS_FTOF) + endif end if else if (CCTK_EQUALS(recon_vars,"conservative")) then - call Con2PrimBoundsPolytype(CCTK_PASS_FTOF) + if(evolve_mhd.ne.0) then + call Con2PrimBoundsPolytype(CCTK_PASS_FTOF) + else + call Con2PrimBoundsPolytype(CCTK_PASS_FTOF) + endif else - call CCTK_WARN(0,"Variable type to reconstruct not recognized.") + call CCTK_WARN(0,"Variable type to reconstruct not recognized.") end if - + return - + end subroutine ReconstructionPolytype diff --git a/src/GRHydro_RegisterGZ.cc b/src/GRHydro_RegisterGZ.cc index 30fe5a4..da6f76e 100644 --- a/src/GRHydro_RegisterGZ.cc +++ b/src/GRHydro_RegisterGZ.cc @@ -10,6 +10,7 @@ #include "cctk.h" #include "cctk_Arguments.h" +#include "cctk_Parameters.h" using namespace std; @@ -26,33 +27,75 @@ using namespace std; // extern "C"void GRHydro_register_GZPatchSystem(CCTK_ARGUMENTS) { + DECLARE_CCTK_PARAMETERS; if (CCTK_IsFunctionAliased("GZPatchSystem_register_sync")) { CCTK_VInfo(CCTK_THORNSTRING, "registering to-be-interpatch-synchronized variables " "with GZPatchSystem"); - - string var[8] = {"HydroBase::rho", "HydroBase::press", "HydroBase::eps", - "HydroBase::vel", - "GRHydro::dens", "GRHydro::tau", "GRHydro::w_lorentz", - "GRHydro::scon"}; - for (int i = 0; i < 8; i++) - { - int status = - GZPatchSystem_register_sync(var[i].c_str()); - if (status < 0) - { - CCTK_VWarn(0, __LINE__, __FILE__, CCTK_THORNSTRING, - "***** GRHydro_register_GZPatchSystem():\n" - " error registering var group %s to be " - "interpatch-synchronized!\n" - " (GZPatchSystem_register_sync() error code %d)\n", - var[i].c_str(), status); - } - } - } + if(CCTK_EQUALS(Bvec_evolution_method,"GRHydro")) { + if(clean_divergence) { + string var[10] = {"HydroBase::rho", "HydroBase::press", "HydroBase::eps", + "HydroBase::vel", + "GRHydro::dens", "GRHydro::tau", "GRHydro::w_lorentz", + "GRHydro::scon", "HydroBase::Bvec", "GRHydro::psidc"}; + for (int i = 0; i < 10; i++) + { + int status = + GZPatchSystem_register_sync(var[i].c_str()); + if (status < 0) + { + CCTK_VWarn(0, __LINE__, __FILE__, CCTK_THORNSTRING, + "***** GRHydro_register_GZPatchSystem():\n" + " error registering var group %s to be " + "interpatch-synchronized!\n" + " (GZPatchSystem_register_sync() error code %d)\n", + var[i].c_str(), status); + } + } + } else { + string var[9] = {"HydroBase::rho", "HydroBase::press", "HydroBase::eps", + "HydroBase::vel", + "GRHydro::dens", "GRHydro::tau", "GRHydro::w_lorentz", + "GRHydro::scon", "HydroBase::Bvec"}; + for (int i = 0; i < 9; i++) + { + int status = + GZPatchSystem_register_sync(var[i].c_str()); + if (status < 0) + { + CCTK_VWarn(0, __LINE__, __FILE__, CCTK_THORNSTRING, + "***** GRHydro_register_GZPatchSystem():\n" + " error registering var group %s to be " + "interpatch-synchronized!\n" + " (GZPatchSystem_register_sync() error code %d)\n", + var[i].c_str(), status); + } + } + } + } else { + string var[8] = {"HydroBase::rho", "HydroBase::press", "HydroBase::eps", + "HydroBase::vel", + "GRHydro::dens", "GRHydro::tau", "GRHydro::w_lorentz", + "GRHydro::scon"}; + for (int i = 0; i < 8; i++) + { + int status = + GZPatchSystem_register_sync(var[i].c_str()); + if (status < 0) + { + CCTK_VWarn(0, __LINE__, __FILE__, CCTK_THORNSTRING, + "***** GRHydro_register_GZPatchSystem():\n" + " error registering var group %s to be " + "interpatch-synchronized!\n" + " (GZPatchSystem_register_sync() error code %d)\n", + var[i].c_str(), status); + } + } + } + } else { CCTK_WARN(1, "Function GZPatchSystem_register_sync not registered!"); @@ -62,31 +105,78 @@ extern "C"void GRHydro_register_GZPatchSystem(CCTK_ARGUMENTS) { CCTK_VInfo(CCTK_THORNSTRING, "registering to-be-cxformed variables with GZPatchSystem"); - - string var[11] = {"HydroBase::rho", "HydroBase::press", "HydroBase::eps", - "HydroBase::vel", - "GRHydro::dens", "GRHydro::tau", "GRHydro::w_lorentz", - "GRHydro::scon", - "ADMBase::metric", "ADMBase::curv", "ADMBase::shift"}; - for (int i = 0; i < 11; i++) - for (int j = 0; j < 3; j++) - { - int ps_status = - GZPatchSystem_register_cxform(j, var[i].c_str()); - if (ps_status < 0) - { - CCTK_VWarn(0, __LINE__, __FILE__, CCTK_THORNSTRING, - "***** GRHydro_register_GZPatchSystem():\n" - " error registering var group %s to be " - "interpatch-cxformhronized!\n" - " (GZPatchSystem_register_cxform() error code %d)\n", - var[i].c_str(), ps_status); + + if(CCTK_EQUALS(Bvec_evolution_method,"GRHydro")) { + if(clean_divergence) { + string var[13] = {"HydroBase::rho", "HydroBase::press", "HydroBase::eps", + "HydroBase::vel", + "GRHydro::dens", "GRHydro::tau", "GRHydro::w_lorentz", + "GRHydro::scon", "HydroBase::Bvec", "GRHydro::psidc", + "ADMBase::metric", "ADMBase::curv", "ADMBase::shift"}; + for (int i = 0; i < 13; i++) + for (int j = 0; j < 3; j++) + { + int ps_status = + GZPatchSystem_register_cxform(j, var[i].c_str()); + if (ps_status < 0) + { + CCTK_VWarn(0, __LINE__, __FILE__, CCTK_THORNSTRING, + "***** GRHydro_register_GZPatchSystem():\n" + " error registering var group %s to be " + "interpatch-cxformhronized!\n" + " (GZPatchSystem_register_cxform() error code %d)\n", + var[i].c_str(), ps_status); + } + } + } else { + string var[12] = {"HydroBase::rho", "HydroBase::press", "HydroBase::eps", + "HydroBase::vel", + "GRHydro::dens", "GRHydro::tau", "GRHydro::w_lorentz", + "GRHydro::scon", "HydroBase::Bvec", + "ADMBase::metric", "ADMBase::curv", "ADMBase::shift"}; + for (int i = 0; i < 12; i++) + for (int j = 0; j < 3; j++) + { + int ps_status = + GZPatchSystem_register_cxform(j, var[i].c_str()); + if (ps_status < 0) + { + CCTK_VWarn(0, __LINE__, __FILE__, CCTK_THORNSTRING, + "***** GRHydro_register_GZPatchSystem():\n" + " error registering var group %s to be " + "interpatch-cxformhronized!\n" + " (GZPatchSystem_register_cxform() error code %d)\n", + var[i].c_str(), ps_status); + } + } } - } + + } else { + string var[11] = {"HydroBase::rho", "HydroBase::press", "HydroBase::eps", + "HydroBase::vel", + "GRHydro::dens", "GRHydro::tau", "GRHydro::w_lorentz", + "GRHydro::scon", + "ADMBase::metric", "ADMBase::curv", "ADMBase::shift"}; + for (int i = 0; i < 11; i++) + for (int j = 0; j < 3; j++) + { + int ps_status = + GZPatchSystem_register_cxform(j, var[i].c_str()); + if (ps_status < 0) + { + CCTK_VWarn(0, __LINE__, __FILE__, CCTK_THORNSTRING, + "***** GRHydro_register_GZPatchSystem():\n" + " error registering var group %s to be " + "interpatch-cxformhronized!\n" + " (GZPatchSystem_register_cxform() error code %d)\n", + var[i].c_str(), ps_status); + } + } + } } else - { - CCTK_WARN(1, "Function GZPatchSystem_register_cxform not registered!"); - } + { + CCTK_WARN(1, "Function GZPatchSystem_register_cxform not registered!"); + } } diff --git a/src/GRHydro_RegisterVars.cc b/src/GRHydro_RegisterVars.cc index 5600d2b..1e35036 100644 --- a/src/GRHydro_RegisterVars.cc +++ b/src/GRHydro_RegisterVars.cc @@ -62,6 +62,13 @@ extern "C"void GRHydro_Register(CCTK_ARGUMENTS) register_evolved("GRHydro::dens", "GRHydro::densrhs"); register_evolved("GRHydro::scon", "GRHydro::srhs"); + if (CCTK_EQUALS(Bvec_evolution_method, "GRHydro")) { + register_evolved("HydroBase::Bvec", "GRHydro::Bvecrhs"); + if(clean_divergence) { + register_evolved("GRHydro::psidc" , "GRHydro::psidcrhs"); + } + } + // tau if (CCTK_EQUALS(GRHydro_eos_type, "General")) register_evolved("GRHydro::tau" , "GRHydro::taurhs"); @@ -115,6 +122,13 @@ extern "C"void GRHydro_Register(CCTK_ARGUMENTS) register_constrained("GRHydro::dens"); register_constrained("GRHydro::scon"); register_constrained("GRHydro::tau"); + + if (CCTK_EQUALS(Bvec_evolution_method, "GRHydro")) { + register_constrained("HydroBase::Bvec"); + if(clean_divergence) { + register_constrained("GRHydro::psidc"); + } + } } } diff --git a/src/GRHydro_TmunuM.F90 b/src/GRHydro_TmunuM.F90 index 5464e45..6315a44 100644 --- a/src/GRHydro_TmunuM.F90 +++ b/src/GRHydro_TmunuM.F90 @@ -74,16 +74,17 @@ DECLARE_CCTK_ARGUMENTS - CCTK_REAL velxlow, velylow, velzlow - CCTK_REAL betaxlow, betaylow, betazlow, beta2 - CCTK_REAL Bvecxlow,Bvecylow,Bveczlow - CCTK_REAL bdotv,b2,bxlow,bylow,bzlow,btlow,dum - CCTK_REAL utlow,rhohstarw2,pstar - CCTK_REAL bdotbeta,vdotbeta - CCTK_INT i,j,k + CCTK_REAL :: velxlow, velylow, velzlow + CCTK_REAL :: betaxlow, betaylow, betazlow, beta2 + CCTK_REAL :: Bvecxlow,Bvecylow,Bveczlow + CCTK_REAL :: bdotv,b2,bxlow,bylow,bzlow,btlow,dum + CCTK_REAL :: utlow,rhohstarw2,pstar + CCTK_REAL :: bdotbeta,vdotbeta + CCTK_INT :: i,j,k !$OMP PARALLEL DO PRIVATE(i,j,velxlow, velylow, velzlow,& + !$OMP Bvecxlow,Bvecylow,Bveczlow, bdotv,dum,b2,bxlow,bylow,bzlow,& !$OMP betaxlow, betaylow, betazlow, beta2, bdotbeta,vdotbeta,utlow, btlow,& !$OMP rhohstarw2,pstar) diff --git a/src/make.code.defn b/src/make.code.defn index 30e5805..872b058 100644 --- a/src/make.code.defn +++ b/src/make.code.defn @@ -7,7 +7,7 @@ SRCS = Utils.F90 \ GRHydro_Boundaries.F90 \ GRHydro_CalcUpdate.F90 \ GRHydro_Con2Prim.F90 \ - GRHydro_DivergenceClean.F90 \ + GRHydro_DivergenceClean.F90 \ GRHydro_Eigenproblem.F90 \ GRHydro_Eigenproblem_Marquina.F90 \ GRHydro_ENOReconstruct.F90 \ @@ -45,10 +45,6 @@ SRCS = Utils.F90 \ GRHydro_Differences.F90 \ GRHydro_EoSChangeGamma.F90 \ GRHydro_Tmunu.F90 \ - GRHydro_RegisterGZM.cc \ - GRHydro_RegisterVarsM.cc \ - GRHydro_BoundariesM.F90 \ - GRHydro_CalcUpdateM.F90 \ GRHydro_Con2PrimM.F90 \ GRHydro_Con2PrimM_pt.c \ GRHydro_Con2PrimM_polytype_pt.c \ @@ -59,13 +55,13 @@ SRCS = Utils.F90 \ GRHydro_PPMM.F90 \ GRHydro_Prim2ConM.F90 \ GRHydro_RiemannSolveM.F90 \ - GRHydro_ReconstructM.F90 \ - GRHydro_ReconstructPolyM.F90 \ GRHydro_SourceM.F90 \ GRHydro_TmunuM.F90 \ GRHydro_UpdateMaskM.F90 \ GRHydro_UtilsM.F90 +### GRHydro_RegisterGZM.cc \ +### GRHydro_RegisterVarsM.cc \ ### GRHydro_Weights.c \ |