diff options
| -rw-r--r-- | interface.ccl | 54 | ||||
| -rw-r--r-- | schedule.ccl | 9 | ||||
| -rw-r--r-- | src/GRHydro_Boundaries.F90 | 25 | ||||
| -rw-r--r-- | src/GRHydro_CalcUpdate.F90 | 43 | ||||
| -rw-r--r-- | src/GRHydro_Con2PrimM.F90 | 82 | ||||
| -rw-r--r-- | src/GRHydro_Con2PrimM_polytype_pt.c | 25 | ||||
| -rw-r--r-- | src/GRHydro_Con2PrimM_pt.c | 27 | ||||
| -rw-r--r-- | src/GRHydro_FluxM.F90 | 17 | ||||
| -rw-r--r-- | src/GRHydro_HLLEM.F90 | 157 | ||||
| -rw-r--r-- | src/GRHydro_InterfacesM.h | 20 | ||||
| -rw-r--r-- | src/GRHydro_Prim2ConM.F90 | 56 | ||||
| -rw-r--r-- | src/GRHydro_Reconstruct.F90 | 157 | ||||
| -rw-r--r-- | src/GRHydro_ReconstructPoly.F90 | 241 | ||||
| -rw-r--r-- | src/GRHydro_RegisterVars.cc | 5 | ||||
| -rw-r--r-- | src/GRHydro_SourceM.F90 | 81 | ||||
| -rw-r--r-- | src/GRHydro_UpdateMaskM.F90 | 20 |
16 files changed, 646 insertions, 373 deletions
diff --git a/interface.ccl b/interface.ccl index de24a28..50d791a 100644 --- a/interface.ccl +++ b/interface.ccl @@ -73,17 +73,18 @@ void FUNCTION Con2PrimGen(CCTK_INT INOUT handle, CCTK_REAL INOUT dens, \ 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 tau, CCTK_REAL INOUT Bconsx, CCTK_REAL INOUT Bconsy, CCTK_REAL INOUT Bconsz, \ + 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 Bvecx, CCTK_REAL INOUT Bvecy, CCTK_REAL INOUT Bvecz, \ + CCTK_REAL INOUT Bvecsq, \ 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) @@ -101,6 +102,23 @@ void FUNCTION Con2PrimPoly(CCTK_INT IN handle, CCTK_REAL OUT dens, \ CCTK_REAL IN r, CCTK_REAL IN rho_min, \ CCTK_INT IN GRHydro_reflevel, CCTK_REAL IN GRHydro_C2P_failed) +void FUNCTION Con2PrimPolyM(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 sc, CCTK_REAL INOUT Bconsx, CCTK_REAL INOUT Bconsy, CCTK_REAL INOUT Bconsz, \ + 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 Bvecx, CCTK_REAL INOUT Bvecy, CCTK_REAL INOUT Bvecz, \ + CCTK_REAL INOUT Bvecsq, \ + 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_INT OUT epsnegative, \ + CCTK_REAL OUT retval) + #void FUNCTION Prim2Con CCTK_INT handle, CCTK_REAL gxx, CCTK_REAL gxy, CCTK_REAL gxz, CCTK_REAL gyy, CCTK_REAL gyz, CCTK_REAL gzz, CCTK_REAL det, CCTK_REAL dens, CCTK_REAL sx, CCTK_REAL sy, CCTK_REAL sz, CCTK_REAL tau, CCTK_REAL rho, CCTK_REAL velx, CCTK_REAL vely, CCTK_REAL velz, CCTK_REAL epsilon, CCTK_REAL press, CCTK_REAL w_lorentz void FUNCTION Prim2ConGen(CCTK_INT IN handle, \ @@ -134,11 +152,12 @@ void FUNCTION Prim2ConGenM(CCTK_INT IN handle, \ CCTK_REAL IN det, CCTK_REAL OUT dens, \ CCTK_REAL OUT sx, CCTK_REAL OUT sy, \ CCTK_REAL OUT sz, CCTK_REAL OUT tau, \ - CCTK_REAL IN Bvecx, CCTK_REAL IN Bvecy, \ - CCTK_REAL IN Bvecz, CCTK_REAL IN rho, CCTK_REAL IN velx, \ + CCTK_REAL OUT Bconsx, CCTK_REAL OUT Bconsy, \ + CCTK_REAL OUT Bconsz, CCTK_REAL IN rho, CCTK_REAL IN velx, \ CCTK_REAL IN vely, \ CCTK_REAL IN velz, CCTK_REAL IN epsilon, \ - CCTK_REAL OUT press, CCTK_REAL OUT w_lorentz) + CCTK_REAL OUT press, CCTK_REAL IN Bvecx, CCTK_REAL IN Bvecy, \ + CCTK_REAL IN Bvecz, CCTK_REAL OUT w_lorentz) void FUNCTION Prim2ConPolyM(CCTK_INT IN handle, \ CCTK_REAL IN gxx, CCTK_REAL IN gxy, \ @@ -147,11 +166,12 @@ void FUNCTION Prim2ConPolyM(CCTK_INT IN handle, \ CCTK_REAL IN det, CCTK_REAL OUT dens, \ CCTK_REAL OUT sx, CCTK_REAL OUT sy, \ CCTK_REAL OUT sz, CCTK_REAL OUT tau, \ - CCTK_REAL IN Bvecx, CCTK_REAL IN Bvecy, \ - CCTK_REAL IN Bvecz, CCTK_REAL IN rho, CCTK_REAL IN velx, \ + CCTK_REAL OUT Bconsx, CCTK_REAL OUT Bconsy, \ + CCTK_REAL OUT Bconsz, CCTK_REAL IN rho, CCTK_REAL IN velx, \ CCTK_REAL IN vely, \ CCTK_REAL IN velz, CCTK_REAL OUT epsilon, \ - CCTK_REAL OUT press, CCTK_REAL OUT w_lorentz) + CCTK_REAL OUT press, CCTK_REAL IN Bvecx, CCTK_REAL IN Bvecy, \ + CCTK_REAL IN Bvecz, CCTK_REAL OUT w_lorentz) PROVIDES FUNCTION SpatialDet WITH SpatialDeterminant LANGUAGE Fortran PROVIDES FUNCTION UpperMet WITH UpperMetric LANGUAGE Fortran @@ -159,6 +179,7 @@ PROVIDES FUNCTION UpperMet WITH UpperMetric LANGUAGE Fortran PROVIDES FUNCTION Con2PrimPoly WITH Con2Prim_ptPolytype LANGUAGE Fortran PROVIDES FUNCTION Con2PrimGenM WITH GRHydro_Con2PrimM_pt LANGUAGE Fortran PROVIDES FUNCTION Con2PrimGen WITH Con2Prim_pt LANGUAGE Fortran +PROVIDES FUNCTION Con2PrimPolyM WITH GRHydro_Con2PrimM_Polytype_pt LANGUAGE Fortran PROVIDES FUNCTION Prim2ConGen WITH prim2con LANGUAGE Fortran PROVIDES FUNCTION Prim2ConPoly WITH prim2conpolytype LANGUAGE Fortran PROVIDES FUNCTION Prim2ConGenM WITH prim2conM LANGUAGE Fortran @@ -310,6 +331,8 @@ real tau type = GF Timelevels = 3 tags='ProlongationParameter="HydroBase::prolo real scon[3] type = GF Timelevels = 3 tags='ProlongationParameter="HydroBase::prolongation_type" tensortypealias="D" tensorweight=+1.0 interpolator="matter"' "generalized momenta" +real Bcons[3] type = GF Timelevels = 3 tags='ProlongationParameter="HydroBase::prolongation_type" tensortypealias="D" tensorweight=+1.0 interpolator="matter"' "B-field conservative variable" + real Y_e_con type = GF Timelevels = 3 tags='ProlongationParameter="HydroBase::prolongation_type" tensortypealias="Scalar" tensorweight=+1.0 interpolator="matter"' "Conserved electron fraction" #real bcom[3] type = GF Timelevels = 3 tags='Prolongation="none" tensortypealias="D" tensorweight=+1.0 interpolator="matter"' "comoving magnetic field components" @@ -328,7 +351,7 @@ real psidc type = GF Timelevels = 3 tags='ProlongationParameter="HydroBase::prol real densrhs type = GF Timelevels = 1 tags='Prolongation="None" checkpoint="no"' "Update term for dens" real taurhs type = GF Timelevels = 1 tags='Prolongation="None" checkpoint="no"' "Update term for tau" real srhs[3] type = GF Timelevels = 1 tags='Prolongation="None" checkpoint="no"' "Update term for s" -real Bvecrhs[3] type = GF Timelevels = 1 tags='Prolongation="None" checkpoint="no"' "Update term for Bvec" +real Bconsrhs[3] type = GF Timelevels = 1 tags='Prolongation="None" checkpoint="no"' "Update term for Bcons" real psidcrhs type = GF Timelevels = 1 tags='Prolongation="None" checkpoint="no"' "Update term for psidc" @@ -391,7 +414,7 @@ real GRHydro_fluxes type = GF Timelevels = 1 tags='Prolongation="None" checkpoin real GRHydro_Bfluxes type = GF Timelevels = 1 tags='Prolongation="None" checkpoint="no"' { - Bvecxflux, Bvecyflux, Bveczflux + Bconsxflux, Bconsyflux, Bconszflux } "Fluxes for each B-field variable" real GRHydro_psifluxes type = GF Timelevels = 1 tags='Prolongation="None" checkpoint="no"' @@ -417,13 +440,18 @@ real GRHydro_con_bext type = GF Timelevels = 1 tags='Prolongation="None" checkpo real GRHydro_MHD_con_bext type = GF Timelevels = 1 tags='Prolongation="None" checkpoint="no"' { - Bvecxplus,Bvecyplus,Bveczplus,Bvecxminus,Bvecyminus,Bveczminus + Bconsxplus,Bconsyplus,Bconszplus,Bconsxminus,Bconsyminus,Bconszminus } "Conservative variables extended to the cell boundaries" +real GRHydro_MHD_prim_bext type = GF Timelevels = 1 tags='Prolongation="None" checkpoint="no"' +{ + Bvecxplus,Bvecyplus,Bveczplus,Bvecxminus,Bvecyminus,Bveczminus +} "Primitive mhd variables extended to the cell boundaries" + real GRHydro_MHD_psidc_bext type = GF Timelevels = 1 tags='Prolongation="None" checkpoint="no"' { psidcplus,psidcminus -} "Conservative variables extended to the cell boundaries for diverence cleaning" +} "Divergence cleaning variable extended to the cell boundaries for diverence cleaning" # real fluxweightvolume type = GF Timelevels = 1 # { diff --git a/schedule.ccl b/schedule.ccl index d9e8552..163818f 100644 --- a/schedule.ccl +++ b/schedule.ccl @@ -30,6 +30,7 @@ if (timelevels == 3) if(CCTK_Equals(Bvec_evolution_method,"GRHydro")) { STORAGE: HydroBase::Bvec[3] + STORAGE: GRHydro::Bcons[3] if (clean_divergence) { STORAGE:psidc[3] @@ -55,6 +56,7 @@ else if(CCTK_Equals(Bvec_evolution_method,"GRHydro")) { STORAGE: HydroBase::Bvec[2] + STORAGE: GRHydro::Bcons[2] if (clean_divergence) { STORAGE:psidc[2] @@ -70,7 +72,7 @@ STORAGE:taurhs STORAGE:srhs if(CCTK_Equals(Bvec_evolution_method,"GRHydro")) { - STORAGE:Bvecrhs + STORAGE:Bconsrhs if (clean_divergence) { STORAGE:psidcrhs @@ -596,6 +598,7 @@ if (CCTK_Equals(method_type, "RSA FV")) STORAGE:GRHydro_con_bext STORAGE:GRHydro_fluxes STORAGE:GRHydro_MHD_con_bext + STORAGE:GRHydro_MHD_prim_bext STORAGE:GRHydro_MHD_psidc_bext STORAGE:GRHydro_Bfluxes STORAGE:GRHydro_psifluxes @@ -626,6 +629,7 @@ if (CCTK_Equals(method_type, "RSA FV")) STORAGE:GRHydro_con_bext STORAGE:GRHydro_fluxes STORAGE:GRHydro_MHD_con_bext + STORAGE:GRHydro_MHD_prim_bext STORAGE:GRHydro_MHD_psidc_bext STORAGE:GRHydro_Bfluxes STORAGE:GRHydro_psifluxes @@ -928,14 +932,15 @@ schedule GRHydro_Boundaries IN HydroBase_Select_Boundaries AS GRHydro_Bound SYNC: HydroBase::press SYNC: HydroBase::eps SYNC: HydroBase::vel + SYNC: Bcons SYNC: HydroBase::Bvec + SYNC: psidc SYNC: GRHydro_cons_tracers SYNC: GRHydro_tracers SYNC: hydrobase::temperature SYNC: hydrobase::entropy SYNC: hydrobase::Y_e SYNC: Y_e_con - SYNC: HydroBase::Bvec } "Select GRHydro boundary conditions" ############################################################ diff --git a/src/GRHydro_Boundaries.F90 b/src/GRHydro_Boundaries.F90 index 9e95adc..765755b 100644 --- a/src/GRHydro_Boundaries.F90 +++ b/src/GRHydro_Boundaries.F90 @@ -24,6 +24,9 @@ #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) +#define Bconsx(i,j,k) Bcons(i,j,k,1) +#define Bconsy(i,j,k) Bcons(i,j,k,2) +#define Bconsz(i,j,k) Bcons(i,j,k,3) /*@@ @routine GRHydro_InitSymBound @@ -93,8 +96,10 @@ 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]") - + if(evolve_mhd.ne.0) then + call SetCartSymVN(ierr, cctkGH, sym, "HydroBase::Bvec[0]") + call SetCartSymVN(ierr, cctkGH, sym, "GRHydro::Bcons[0]") + endif sym(1) = 1 sym(2) = -1 @@ -102,15 +107,21 @@ 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]") - + if(evolve_mhd.ne.0) then + call SetCartSymVN(ierr, cctkGH, sym, "HydroBase::Bvec[1]") + call SetCartSymVN(ierr, cctkGH, sym, "GRHydro::Bcons[1]") + endif + sym(1) = 1 sym(2) = 1 sym(3) = -1 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]") + if(evolve_mhd.ne.0) then + call SetCartSymVN(ierr, cctkGH, sym, "HydroBase::Bvec[2]") + call SetCartSymVN(ierr, cctkGH, sym, "GRHydro::Bcons[2]") + endif end subroutine GRHydro_InitSymBound @@ -169,6 +180,8 @@ subroutine GRHydro_Boundaries(CCTK_ARGUMENTS) if(evolve_mhd.ne.0) then ierr = Boundary_SelectGroupForBC(cctkGH, faces, GRHydro_stencil, -ione, & "HydroBase::Bvec", "Flat") + ierr = Boundary_SelectGroupForBC(cctkGH, faces, GRHydro_stencil, -ione, & + "GRHydro::Bcons", "Flat") if(clean_divergence.ne.0) then ierr = Boundary_SelectGroupForBC(cctkGH, faces, GRHydro_stencil, -ione, & "GRHydro::psidc", "Flat") @@ -219,6 +232,8 @@ subroutine GRHydro_Boundaries(CCTK_ARGUMENTS) if(evolve_mhd.ne.0) then ierr = Boundary_SelectGroupForBC(cctkGH, faces, GRHydro_stencil, -ione, & "HydroBase::Bvec", "None") + ierr = Boundary_SelectGroupForBC(cctkGH, faces, GRHydro_stencil, -ione, & + "GRHydro::Bcons", "None") if(clean_divergence.ne.0) then ierr = Boundary_SelectGroupForBC(cctkGH, faces, GRHydro_stencil, -ione, & "GRHydro::psidc", "None") diff --git a/src/GRHydro_CalcUpdate.F90 b/src/GRHydro_CalcUpdate.F90 index 452b070..1b3b96d 100644 --- a/src/GRHydro_CalcUpdate.F90 +++ b/src/GRHydro_CalcUpdate.F90 @@ -81,15 +81,15 @@ subroutine UpdateCalculation(CCTK_ARGUMENTS) (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 + Bconsrhs(i,j,k,1) = Bconsrhs(i,j,k,1) + & + (alp_l * Bconsxflux(i-xoffset,j-yoffset,k-zoffset) - & + alp_r * Bconsxflux(i,j,k)) * idx + Bconsrhs(i,j,k,2) = Bconsrhs(i,j,k,2) + & + (alp_l * Bconsyflux(i-xoffset,j-yoffset,k-zoffset) - & + alp_r * Bconsyflux(i,j,k)) * idx + Bconsrhs(i,j,k,3) = Bconsrhs(i,j,k,3) + & + (alp_l * Bconszflux(i-xoffset,j-yoffset,k-zoffset) - & + alp_r * Bconszflux(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) - & @@ -102,6 +102,13 @@ subroutine UpdateCalculation(CCTK_ARGUMENTS) Bvec(i+xoffset,j+yoffset,k+zoffset,flux_direction)) divB(i,j,k) = divB(i,j,k) + ( alp_l * Bvec_l - alp_r * Bvec_r ) * idx endif + +!!$ if(i.eq.75.and.j.eq.5.and.k.eq.5)write(6,*)'update:', & +!!$ Bconsrhs(i,j,k,1),Bconsxflux(i,j,k), & +!!$ Bconsxflux(i-xoffset,j-yoffset,k-zoffset), & +!!$ Bconsrhs(i,j,k,3),Bconszflux(i,j,k), & +!!$ Bconszflux(i-xoffset,j-yoffset,k-zoffset) + endif if (evolve_tracer .ne. 0) then @@ -241,15 +248,15 @@ subroutine UpdateCalculation(CCTK_ARGUMENTS) (tauflux(i-xoffset,j-yoffset,k-zoffset) - & 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 + Bconsrhs(i,j,k,1) = Bconsrhs(i,j,k,1) + & + (Bconsxflux(i-xoffset,j-yoffset,k-zoffset) - & + Bconsxflux(i,j,k)) * idx + Bconsrhs(i,j,k,2) = Bconsrhs(i,j,k,2) + & + (Bconsyflux(i-xoffset,j-yoffset,k-zoffset) - & + Bconsyflux(i,j,k)) * idx + Bconsrhs(i,j,k,3) = Bconsrhs(i,j,k,3) + & + (Bconszflux(i-xoffset,j-yoffset,k-zoffset) - & + Bconszflux(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) - & diff --git a/src/GRHydro_Con2PrimM.F90 b/src/GRHydro_Con2PrimM.F90 index 582eec5..e062976 100644 --- a/src/GRHydro_Con2PrimM.F90 +++ b/src/GRHydro_Con2PrimM.F90 @@ -129,8 +129,7 @@ subroutine Conservative2PrimitiveM(CCTK_ARGUMENTS) if(evolve_Y_e.ne.0) then Y_e(i,j,k) = Y_e_con(i,j,k) / dens(i,j,k) endif - - + if ( dens(i,j,k) .le. sqrt(det)*GRHydro_rho_min*(1.d0+GRHydro_atmo_tolerance) ) then b2=gxx(i,j,k)*Bvec(i,j,k,1)**2+gyy(i,j,k)*Bvec(i,j,k,2)**2+gzz(i,j,k)*Bvec(i,j,k,3)**2+ & @@ -163,11 +162,11 @@ subroutine Conservative2PrimitiveM(CCTK_ARGUMENTS) if(evolve_temper.eq.0) then call GRHydro_Con2PrimM_pt(GRHydro_eos_handle, local_gam, dens(i,j,k), & scon(i,j,k,1),scon(i,j,k,2),scon(i,j,k,3), tau(i,j,k), & - rho(i,j,k),& - vel(i,j,k,1),vel(i,j,k,2),vel(i,j,k,3),eps(i,j,k),press(i,j,k),w_lorentz(i,j,k),& + Bcons(i,j,k,1), Bcons(i,j,k,2), Bcons(i,j,k,3), rho(i,j,k),& + vel(i,j,k,1),vel(i,j,k,2),vel(i,j,k,3),eps(i,j,k),press(i,j,k),& + Bvec(i,j,k,1), Bvec(i,j,k,2), Bvec(i,j,k,3),b2, w_lorentz(i,j,k),& gxx(i,j,k),gxy(i,j,k),gxz(i,j,k),gyy(i,j,k),gyz(i,j,k),gzz(i,j,k), & uxx,uxy,uxz,uyy,uyz,uzz,det, & - Bvec(i,j,k,1), Bvec(i,j,k,2),Bvec(i,j,k,3),b2,& epsnegative,GRHydro_C2P_failed(i,j,k)) else stop "Please implement finite T con2prim routine in MHD part!" @@ -239,13 +238,13 @@ subroutine Conservative2PrimitiveM(CCTK_ARGUMENTS) call GRHydro_Con2PrimM_Polytype_pt(GRHydro_eos_handle, local_pgam, dens(i,j,k), & scon(i,j,k,1),scon(i,j,k,2),scon(i,j,k,3), sc, & - rho(i,j,k),& - vel(i,j,k,1),vel(i,j,k,2),vel(i,j,k,3),eps(i,j,k),press(i,j,k),w_lorentz(i,j,k),& + Bcons(i,j,k,1), Bcons(i,j,k,2), Bcons(i,j,k,3), rho(i,j,k),& + vel(i,j,k,1),vel(i,j,k,2),vel(i,j,k,3),eps(i,j,k),press(i,j,k),& + Bvec(i,j,k,1), Bvec(i,j,k,2), Bvec(i,j,k,3),b2, w_lorentz(i,j,k),& gxx(i,j,k),gxy(i,j,k),gxz(i,j,k),gyy(i,j,k),gyz(i,j,k),gzz(i,j,k), & uxx,uxy,uxz,uyy,uyz,uzz,det, & - Bvec(i,j,k,1), Bvec(i,j,k,2),Bvec(i,j,k,3),b2,& epsnegative,GRHydro_C2P_failed(i,j,k)) - + #endif end if @@ -294,7 +293,7 @@ subroutine Conservative2PrimitiveBoundsM(CCTK_ARGUMENTS) uxxr, uxyr, uxzr, uyyr, uyzr, uzzr, pmin, epsmin CCTK_REAL :: gxxl,gxyl,gxzl,gyyl,gyzl,gzzl,avg_detl,& gxxr,gxyr,gxzr,gyyr,gyzr,gzzr,avg_detr - CCTK_REAL :: b2minus, b2plus, local_gam, local_pgam,local_K,sc + CCTK_REAL :: b2minus, b2plus, local_gam, local_pgam,local_K,scminus,scplus CCTK_INT :: epsnegative character(len=100) warnline @@ -387,11 +386,11 @@ subroutine Conservative2PrimitiveBoundsM(CCTK_ARGUMENTS) call GRHydro_Con2PrimM_pt(GRHydro_eos_handle, local_gam,densminus(i,j,k), & sxminus(i,j,k),syminus(i,j,k),szminus(i,j,k), tauminus(i,j,k), & - rhominus(i,j,k),& - velxminus(i,j,k),velyminus(i,j,k),velzminus(i,j,k),epsminus(i,j,k),pressminus(i,j,k),w_lorentzminus(i,j,k),& + Bconsxminus(i,j,k), Bconsyminus(i,j,k), Bconszminus(i,j,k), rhominus(i,j,k),& + velxminus(i,j,k),velyminus(i,j,k),velzminus(i,j,k),epsminus(i,j,k),pressminus(i,j,k),& + Bvecxminus(i,j,k), Bvecyminus(i,j,k), Bveczminus(i,j,k),b2minus, w_lorentzminus(i,j,k),& gxxl,gxyl,gxzl,gyyl,gyzl,gzzl, & uxxl,uxyl,uxzl,uyyl,uyzl,uzzl,avg_detl, & - Bvecxminus(i,j,k), Bvecyminus(i,j,k),Bveczminus(i,j,k),b2minus,& epsnegative,GRHydro_C2P_failed(i,j,k)) if (epsnegative .ne. 0) then @@ -407,16 +406,17 @@ subroutine Conservative2PrimitiveBoundsM(CCTK_ARGUMENTS) call EOS_Omni_press(GRHydro_polytrope_handle,keytemp,GRHydro_eos_rf_prec,n,& xrho,1.0d0,xtemp,xye,xpress,keyerr,anyerr) local_pgam=log(xpress/local_K)/log(xrho) - sc = local_K*densminus(i,j,k) + scminus = local_K*densminus(i,j,k) call GRHydro_Con2PrimM_Polytype_pt(GRHydro_eos_handle, local_pgam, densminus(i,j,k), & - sxminus(i,j,k),syminus(i,j,k),szminus(i,j,k), sc, & - rhominus(i,j,k),& - velxminus(i,j,k),velyminus(i,j,k),velzminus(i,j,k),epsminus(i,j,k),pressminus(i,j,k),w_lorentzminus(i,j,k),& + sxminus(i,j,k),syminus(i,j,k),szminus(i,j,k), scminus, & + Bconsxminus(i,j,k), Bconsyminus(i,j,k), Bconszminus(i,j,k), rhominus(i,j,k),& + velxminus(i,j,k),velyminus(i,j,k),velzminus(i,j,k),epsminus(i,j,k),pressminus(i,j,k),& + Bvecxminus(i,j,k), Bvecyminus(i,j,k), Bveczminus(i,j,k),b2minus, w_lorentzminus(i,j,k),& gxxl,gxyl,gxzl,gyyl,gyzl,gzzl, & uxxl,uxyl,uxzl,uyyl,uyzl,uzzl,avg_detl, & - Bvecxminus(i,j,k), Bvecyminus(i,j,k),Bveczminus(i,j,k),b2minus,& epsnegative,GRHydro_C2P_failed(i,j,k)) + end if if (epsminus(i,j,k) .lt. 0.0d0) then @@ -443,13 +443,13 @@ subroutine Conservative2PrimitiveBoundsM(CCTK_ARGUMENTS) epsnegative = 0 call GRHydro_Con2PrimM_pt(GRHydro_eos_handle, local_gam, densplus(i,j,k), & sxplus(i,j,k),syplus(i,j,k),szplus(i,j,k), tauplus(i,j,k),& - rhoplus(i,j,k),& - velxplus(i,j,k),velyplus(i,j,k),velzplus(i,j,k),epsplus(i,j,k),pressplus(i,j,k),w_lorentzplus(i,j,k),& + Bconsxplus(i,j,k), Bconsyplus(i,j,k), Bconszplus(i,j,k), rhoplus(i,j,k),& + velxplus(i,j,k),velyplus(i,j,k),velzplus(i,j,k),epsplus(i,j,k),pressplus(i,j,k),& + Bvecxplus(i,j,k), Bvecyplus(i,j,k), Bveczplus(i,j,k),b2plus, w_lorentzplus(i,j,k),& gxxr,gxyr,gxzr,gyyr,gyzr,gzzr, & uxxr,uxyr,uxzr,uyyr,uyzr,uzzr,avg_detr, & - Bvecxplus(i,j,k), Bvecyplus(i,j,k),Bveczplus(i,j,k),b2plus,& epsnegative,GRHydro_C2P_failed(i,j,k)) - + if (epsnegative .ne. 0) then !$OMP CRITICAL call CCTK_WARN(GRHydro_NaN_verbose+2, 'Specific internal energy just went below 0, trying polytype!!') @@ -463,15 +463,15 @@ subroutine Conservative2PrimitiveBoundsM(CCTK_ARGUMENTS) call EOS_Omni_press(GRHydro_polytrope_handle,keytemp,GRHydro_eos_rf_prec,n,& xrho,1.0d0,xtemp,xye,xpress,keyerr,anyerr) local_pgam=log(xpress/local_K)/log(xrho) - sc = local_K*densplus(i,j,k) + scplus = local_K*densplus(i,j,k) call GRHydro_Con2PrimM_Polytype_pt(GRHydro_eos_handle, local_pgam, densplus(i,j,k), & - sxplus(i,j,k),syplus(i,j,k),szplus(i,j,k), sc,& - rhoplus(i,j,k),& - velxplus(i,j,k),velyplus(i,j,k),velzplus(i,j,k),epsplus(i,j,k),pressplus(i,j,k),w_lorentzplus(i,j,k),& + sxplus(i,j,k),syplus(i,j,k),szplus(i,j,k), scplus,& + Bconsxplus(i,j,k), Bconsyplus(i,j,k), Bconszplus(i,j,k), rhoplus(i,j,k),& + velxplus(i,j,k),velyplus(i,j,k),velzplus(i,j,k),epsplus(i,j,k),pressplus(i,j,k),& + Bvecxplus(i,j,k), Bvecyplus(i,j,k), Bveczplus(i,j,k),b2plus, w_lorentzplus(i,j,k),& gxxr,gxyr,gxzr,gyyr,gyzr,gzzr, & uxxr,uxyr,uxzr,uyyr,uyzr,uzzr,avg_detr, & - Bvecxplus(i,j,k), Bvecyplus(i,j,k),Bveczplus(i,j,k),b2plus,& epsnegative,GRHydro_C2P_failed(i,j,k)) end if @@ -610,14 +610,14 @@ subroutine Conservative2PrimitivePolytypeM(CCTK_ARGUMENTS) sc = local_K*dens(i,j,k) call GRHydro_Con2PrimM_Polytype_pt(GRHydro_eos_handle, local_pgam, dens(i,j,k), & - scon(i,j,k,1),scon(i,j,k,2),scon(i,j,k,3), sc,& - rho(i,j,k),& - vel(i,j,k,1),vel(i,j,k,2),vel(i,j,k,3),eps(i,j,k),press(i,j,k),w_lorentz(i,j,k),& - gxx(i,j,k),gxy(i,j,k),gxz(i,j,k),gyy(i,j,k),gyz(i,j,k),gzz(i,j,k), & - uxx,uxy,uxz,uyy,uyz,uzz,det, & - Bvec(i,j,k,1), Bvec(i,j,k,2),Bvec(i,j,k,3),b2,& - epsnegative,GRHydro_C2P_failed(i,j,k)) - + scon(i,j,k,1),scon(i,j,k,2),scon(i,j,k,3), sc, & + Bcons(i,j,k,1), Bcons(i,j,k,2), Bcons(i,j,k,3), rho(i,j,k),& + vel(i,j,k,1),vel(i,j,k,2),vel(i,j,k,3),eps(i,j,k),press(i,j,k),& + Bvec(i,j,k,1), Bvec(i,j,k,2), Bvec(i,j,k,3),b2, w_lorentz(i,j,k),& + gxx(i,j,k),gxy(i,j,k),gxz(i,j,k),gyy(i,j,k),gyz(i,j,k),gzz(i,j,k), & + uxx,uxy,uxz,uyy,uyz,uzz,det, & + epsnegative,GRHydro_C2P_failed(i,j,k)) + end do end do end do @@ -723,20 +723,20 @@ subroutine Con2PrimBoundsPolytypeM(CCTK_ARGUMENTS) call GRHydro_Con2PrimM_Polytype_pt(GRHydro_eos_handle, local_pgam,densminus(i,j,k), & sxminus(i,j,k),syminus(i,j,k),szminus(i,j,k), scminus,& - rhominus(i,j,k),& - velxminus(i,j,k),velyminus(i,j,k),velzminus(i,j,k),epsminus(i,j,k),pressminus(i,j,k),w_lorentzminus(i,j,k),& + Bconsxminus(i,j,k), Bconsyminus(i,j,k), Bconszminus(i,j,k), rhominus(i,j,k),& + velxminus(i,j,k),velyminus(i,j,k),velzminus(i,j,k),epsminus(i,j,k),pressminus(i,j,k),& + Bvecxminus(i,j,k), Bvecyminus(i,j,k), Bveczminus(i,j,k),b2minus, w_lorentzminus(i,j,k),& gxxl,gxyl,gxzl,gyyl,gyzl,gzzl, & uxxl,uxyl,uxzl,uyyl,uyzl,uzzl,avg_detl, & - Bvecxminus(i,j,k), Bvecyminus(i,j,k),Bveczminus(i,j,k),b2minus,& epsnegative,GRHydro_C2P_failed(i,j,k)) call GRHydro_Con2PrimM_Polytype_pt(GRHydro_eos_handle, local_pgam,densplus(i,j,k), & sxplus(i,j,k),syplus(i,j,k),szplus(i,j,k), scplus,& - rhoplus(i,j,k),& - velxplus(i,j,k),velyplus(i,j,k),velzplus(i,j,k),epsplus(i,j,k),pressplus(i,j,k),w_lorentzplus(i,j,k),& + Bconsxplus(i,j,k), Bconsyplus(i,j,k), Bconszplus(i,j,k), rhoplus(i,j,k),& + velxplus(i,j,k),velyplus(i,j,k),velzplus(i,j,k),epsplus(i,j,k),pressplus(i,j,k),& + Bvecxplus(i,j,k), Bvecyplus(i,j,k),Bveczplus(i,j,k),b2plus,w_lorentzplus(i,j,k),& gxxr,gxyr,gxzr,gyyr,gyzr,gzzr, & uxxr,uxyr,uxzr,uyyr,uyzr,uzzr,avg_detr, & - Bvecxplus(i,j,k), Bvecyplus(i,j,k),Bveczplus(i,j,k),b2plus,& epsnegative,GRHydro_C2P_failed(i,j,k)) end do end do diff --git a/src/GRHydro_Con2PrimM_polytype_pt.c b/src/GRHydro_Con2PrimM_polytype_pt.c index 78445cd..2f3deb2 100644 --- a/src/GRHydro_Con2PrimM_polytype_pt.c +++ b/src/GRHydro_Con2PrimM_polytype_pt.c @@ -103,17 +103,18 @@ void CCTK_FCALL CCTK_FNAME(GRHydro_Con2PrimM_Polytype_pt) ( CCTK_REAL *dens_in, CCTK_REAL *sx_in, CCTK_REAL *sy_in, CCTK_REAL *sz_in, CCTK_REAL *sc_in, + CCTK_REAL *Bconsx_in, CCTK_REAL *Bconsy_in, CCTK_REAL *Bconsz_in, CCTK_REAL *rho, CCTK_REAL *velx, CCTK_REAL *vely, CCTK_REAL *velz, CCTK_REAL *epsilon, CCTK_REAL *pressure, + CCTK_REAL *Bx, CCTK_REAL *By, CCTK_REAL *Bz, + CCTK_REAL *bsq, CCTK_REAL *w_lorentz, CCTK_REAL *gxx, CCTK_REAL *gxy, CCTK_REAL *gxz, CCTK_REAL *gyy, CCTK_REAL *gyz, CCTK_REAL *gzz, CCTK_REAL *uxx, CCTK_REAL *uxy, CCTK_REAL *uxz, CCTK_REAL *uyy, CCTK_REAL *uyz, CCTK_REAL *uzz, CCTK_REAL *det, - CCTK_REAL *Bx, CCTK_REAL *By, CCTK_REAL *Bz, - CCTK_REAL *bsq, CCTK_INT *epsnegative, CCTK_REAL *retval); @@ -184,17 +185,18 @@ void CCTK_FCALL CCTK_FNAME(GRHydro_Con2PrimM_Polytype_pt) ( CCTK_REAL *dens_in, CCTK_REAL *sx_in, CCTK_REAL *sy_in, CCTK_REAL *sz_in, CCTK_REAL *sc_in, + CCTK_REAL *Bconsx_in, CCTK_REAL *Bconsy_in, CCTK_REAL *Bconsz_in, CCTK_REAL *rho, CCTK_REAL *velx, CCTK_REAL *vely, CCTK_REAL *velz, CCTK_REAL *epsilon, CCTK_REAL *pressure, + CCTK_REAL *Bx, CCTK_REAL *By, CCTK_REAL *Bz, + CCTK_REAL *bsq, CCTK_REAL *w_lorentz, CCTK_REAL *gxx, CCTK_REAL *gxy, CCTK_REAL *gxz, CCTK_REAL *gyy, CCTK_REAL *gyz, CCTK_REAL *gzz, CCTK_REAL *uxx, CCTK_REAL *uxy, CCTK_REAL *uxz, CCTK_REAL *uyy, CCTK_REAL *uyz, CCTK_REAL *uzz, CCTK_REAL *det, - CCTK_REAL *Bx, CCTK_REAL *By, CCTK_REAL *Bz, - CCTK_REAL *bsq, CCTK_INT *epsnegative, CCTK_REAL *retval) @@ -228,12 +230,19 @@ void CCTK_FCALL CCTK_FNAME(GRHydro_Con2PrimM_Polytype_pt) ( 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," *w_lorentz = %26.16e \n", *w_lorentz ); fprintf(stdout," *gxx = %26.16e \n", *gxx ); fprintf(stdout," *gxy = %26.16e \n", *gxy ); @@ -248,10 +257,6 @@ void CCTK_FCALL CCTK_FNAME(GRHydro_Con2PrimM_Polytype_pt) ( fprintf(stdout," *uyz = %26.16e \n", *uyz ); fprintf(stdout," *uzz = %26.16e \n", *uzz ); fprintf(stdout," *det = %26.16e \n", *det ); - 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," *epsnegative = %10d \n", *epsnegative ); fprintf(stdout," *retval = %26.16e \n", *retval ); fflush(stdout); @@ -267,6 +272,10 @@ void CCTK_FCALL CCTK_FNAME(GRHydro_Con2PrimM_Polytype_pt) ( usy = (*uxy)*sx + (*uyy)*sy + (*uyz)*sz; usz = (*uxz)*sx + (*uyz)*sy + (*uzz)*sz; + *Bx = (*Bconsx_in) * inv_sqrt_detg; + *By = (*Bconsy_in) * inv_sqrt_detg; + *Bz = (*Bconsz_in) * inv_sqrt_detg; + // Calculate various scalars (Q.B, Q^2, etc) from the conserved variables: lg.Bsq = diff --git a/src/GRHydro_Con2PrimM_pt.c b/src/GRHydro_Con2PrimM_pt.c index cfcf4c0..a2d9853 100644 --- a/src/GRHydro_Con2PrimM_pt.c +++ b/src/GRHydro_Con2PrimM_pt.c @@ -108,18 +108,17 @@ void CCTK_FCALL CCTK_FNAME(GRHydro_Con2PrimM_pt) ( CCTK_INT *handle, CCTK_REAL *gamma_eos, CCTK_REAL *dens_in, CCTK_REAL *sx_in, CCTK_REAL *sy_in, CCTK_REAL *sz_in, - CCTK_REAL *tau_in, + CCTK_REAL *tau_in, CCTK_REAL *Bconsx_in, CCTK_REAL *Bconsy_in, CCTK_REAL *Bconsz_in, CCTK_REAL *rho, CCTK_REAL *velx, CCTK_REAL *vely, CCTK_REAL *velz, CCTK_REAL *epsilon, CCTK_REAL *pressure, + CCTK_REAL *Bx, CCTK_REAL *By, CCTK_REAL *Bz, CCTK_REAL *bsq, CCTK_REAL *w_lorentz, CCTK_REAL *gxx, CCTK_REAL *gxy, CCTK_REAL *gxz, CCTK_REAL *gyy, CCTK_REAL *gyz, CCTK_REAL *gzz, CCTK_REAL *uxx, CCTK_REAL *uxy, CCTK_REAL *uxz, CCTK_REAL *uyy, CCTK_REAL *uyz, CCTK_REAL *uzz, CCTK_REAL *det, - CCTK_REAL *Bx, CCTK_REAL *By, CCTK_REAL *Bz, - CCTK_REAL *bsq, CCTK_INT *epsnegative, CCTK_REAL *retval); @@ -189,18 +188,19 @@ void CCTK_FCALL CCTK_FNAME(GRHydro_Con2PrimM_pt) ( CCTK_INT *handle, CCTK_REAL *gamma_eos, CCTK_REAL *dens_in, CCTK_REAL *sx_in, CCTK_REAL *sy_in, CCTK_REAL *sz_in, - CCTK_REAL *tau_in, + CCTK_REAL *tau_in, + CCTK_REAL *Bconsx_in, CCTK_REAL *Bconsy_in, CCTK_REAL *Bconsz_in, CCTK_REAL *rho, CCTK_REAL *velx, CCTK_REAL *vely, CCTK_REAL *velz, CCTK_REAL *epsilon, CCTK_REAL *pressure, + CCTK_REAL *Bx, CCTK_REAL *By, CCTK_REAL *Bz, + CCTK_REAL *bsq, CCTK_REAL *w_lorentz, CCTK_REAL *gxx, CCTK_REAL *gxy, CCTK_REAL *gxz, CCTK_REAL *gyy, CCTK_REAL *gyz, CCTK_REAL *gzz, CCTK_REAL *uxx, CCTK_REAL *uxy, CCTK_REAL *uxz, CCTK_REAL *uyy, CCTK_REAL *uyz, CCTK_REAL *uzz, CCTK_REAL *det, - CCTK_REAL *Bx, CCTK_REAL *By, CCTK_REAL *Bz, - CCTK_REAL *bsq, CCTK_INT *epsnegative, CCTK_REAL *retval) @@ -231,12 +231,19 @@ void CCTK_FCALL CCTK_FNAME(GRHydro_Con2PrimM_pt) ( 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," *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," *w_lorentz = %26.16e \n", *w_lorentz ); fprintf(stdout," *gxx = %26.16e \n", *gxx ); fprintf(stdout," *gxy = %26.16e \n", *gxy ); @@ -251,10 +258,6 @@ void CCTK_FCALL CCTK_FNAME(GRHydro_Con2PrimM_pt) ( fprintf(stdout," *uyz = %26.16e \n", *uyz ); fprintf(stdout," *uzz = %26.16e \n", *uzz ); fprintf(stdout," *det = %26.16e \n", *det ); - 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," *epsnegative = %10d \n", *epsnegative ); fprintf(stdout," *retval = %26.16e \n", *retval ); fflush(stdout); @@ -271,6 +274,10 @@ void CCTK_FCALL CCTK_FNAME(GRHydro_Con2PrimM_pt) ( usy = (*uxy)*sx + (*uyy)*sy + (*uyz)*sz; usz = (*uxz)*sx + (*uyz)*sy + (*uzz)*sz; + *Bx = (*Bconsx_in) * inv_sqrt_detg; + *By = (*Bconsy_in) * inv_sqrt_detg; + *Bz = (*Bconsz_in) * inv_sqrt_detg; + // Calculate various scalars (Q.B, Q^2, etc) from the conserved variables: lg.Bsq = diff --git a/src/GRHydro_FluxM.F90 b/src/GRHydro_FluxM.F90 index fadbfe6..fa1f753 100644 --- a/src/GRHydro_FluxM.F90 +++ b/src/GRHydro_FluxM.F90 @@ -23,13 +23,10 @@ subroutine num_x_fluxM(dens,sx,sy,sz,tau,Bx,By,Bz,& CCTK_REAL :: vxt,vyt,vzt,bsubx,bsuby,bsubz,ab0,w CCTK_REAL :: det,alp,beta,pressstar CCTK_REAL :: velm - CCTK_REAL :: sdet,psipstar,psiBx,psiBy,psiBz + CCTK_REAL :: sdet,psipstar sdet=sqrt(det) psipstar=pressstar*sdet - psiBx=Bx*sdet - psiBy=By*sdet - psiBz=Bz*sdet !!$ We actually need all three values of vtilde = v^i - beta^i/alp, as well as velm=vxt+beta/alp @@ -38,19 +35,19 @@ subroutine num_x_fluxM(dens,sx,sy,sz,tau,Bx,By,Bz,& !!$ In the notation of Anton et al.: [psi^6 D] vtilde^i densf = dens * vxt - sxf = sx*vxt+psipstar-bsubx*psiBx/w + sxf = sx*vxt+psipstar-bsubx*Bx/w - syf = sy*vxt-bsuby*psiBx/w + syf = sy*vxt-bsuby*Bx/w - szf = sz*vxt-bsubz*psiBx/w + szf = sz*vxt-bsubz*Bx/w !!$ [psi^6 tau] vtilde^i +p* v^i - alp b^0 B^i/w - tauf = tau*vxt + psipstar*velm - ab0*psiBx/w + tauf = tau*vxt + psipstar*velm - ab0*Bx/w !!$ [psi^6 (B^k vtilde^i - B^i vtilde^k)] bxf = 0.0 - byf = psiBy * vxt - psiBx*vyt - bzf = psiBz * vxt - psiBx*vzt + byf = By * vxt - Bx*vyt + bzf = Bz * vxt - Bx*vzt end subroutine num_x_fluxM diff --git a/src/GRHydro_HLLEM.F90 b/src/GRHydro_HLLEM.F90 index 56e19a1..8cb2e49 100644 --- a/src/GRHydro_HLLEM.F90 +++ b/src/GRHydro_HLLEM.F90 @@ -41,9 +41,9 @@ subroutine GRHydro_HLLEM(CCTK_ARGUMENTS) integer :: i, j, k, m CCTK_REAL, dimension(8) :: cons_p,cons_m,fplus,fminus,f1,qdiff - CCTK_REAL, dimension(7) :: prim_p, prim_m + CCTK_REAL, dimension(10) :: prim_p, prim_m CCTK_REAL, dimension(5) :: lamminus,lamplus - CCTK_REAL :: charmin, charmax, charpm,avg_alp,avg_det + CCTK_REAL :: charmin, charmax, charpm,avg_alp,avg_det, sdet CCTK_REAL :: gxxh, gxyh, gxzh, gyyh, gyzh, gzzh, uxxh, uxyh, & uxzh, uyyh, uyzh, uzzh, avg_beta, usendh, alp_l, alp_r, & cs2_p, cs2_m, dpdeps_p, dpdeps_m @@ -98,18 +98,20 @@ subroutine GRHydro_HLLEM(CCTK_ARGUMENTS) cons_p(3) = syplus(i,j,k) cons_p(4) = szplus(i,j,k) cons_p(5) = tauplus(i,j,k) - cons_p(6) = Bvecxplus(i,j,k) - cons_p(7) = Bvecyplus(i,j,k) - cons_p(8) = Bveczplus(i,j,k) + cons_p(6) = Bconsxplus(i,j,k) + cons_p(7) = Bconsyplus(i,j,k) + cons_p(8) = Bconszplus(i,j,k) cons_m(1) = densminus(i+xoffset,j+yoffset,k+zoffset) cons_m(2) = sxminus(i+xoffset,j+yoffset,k+zoffset) cons_m(3) = syminus(i+xoffset,j+yoffset,k+zoffset) cons_m(4) = szminus(i+xoffset,j+yoffset,k+zoffset) cons_m(5) = tauminus(i+xoffset,j+yoffset,k+zoffset) - cons_m(6) = Bvecxminus(i+xoffset,j+yoffset,k+zoffset) - cons_m(7) = Bvecyminus(i+xoffset,j+yoffset,k+zoffset) - cons_m(8) = Bveczminus(i+xoffset,j+yoffset,k+zoffset) + cons_m(6) = Bconsxminus(i+xoffset,j+yoffset,k+zoffset) + cons_m(7) = Bconsyminus(i+xoffset,j+yoffset,k+zoffset) + cons_m(8) = Bconszminus(i+xoffset,j+yoffset,k+zoffset) + +!!$ if(i.eq.75.and.j.eq.5.and.k.eq.5)write(6,*)'HLLEM0:',cons_m(6),cons_p(6),cons_m(7),cons_p(7),cons_m(8),cons_p(8) prim_p(1) = rhoplus(i,j,k) prim_p(2) = velxplus(i,j,k) @@ -118,6 +120,9 @@ subroutine GRHydro_HLLEM(CCTK_ARGUMENTS) prim_p(5) = epsplus(i,j,k) prim_p(6) = pressplus(i,j,k) prim_p(7) = w_lorentzplus(i,j,k) + prim_p(8) = Bvecxplus(i,j,k) + prim_p(9) = Bvecyplus(i,j,k) + prim_p(10) = Bveczplus(i,j,k) prim_m(1) = rhominus(i+xoffset,j+yoffset,k+zoffset) prim_m(2) = velxminus(i+xoffset,j+yoffset,k+zoffset) @@ -126,6 +131,9 @@ subroutine GRHydro_HLLEM(CCTK_ARGUMENTS) prim_m(5) = epsminus(i+xoffset,j+yoffset,k+zoffset) prim_m(6) = pressminus(i+xoffset,j+yoffset,k+zoffset) prim_m(7) = w_lorentzminus(i+xoffset,j+yoffset,k+zoffset) + prim_m(8) = Bvecxminus(i+xoffset,j+yoffset,k+zoffset) + prim_m(9) = Bvecyminus(i+xoffset,j+yoffset,k+zoffset) + prim_m(10)= Bveczminus(i+xoffset,j+yoffset,k+zoffset) if(clean_divergence.ne.0) then psidcp = psidcplus(i,j,k) @@ -171,6 +179,7 @@ subroutine GRHydro_HLLEM(CCTK_ARGUMENTS) gzzh = 0.5d0 * (gzz(i+xoffset,j+yoffset,k+zoffset) + gzz(i,j,k)) avg_det = SPATIAL_DETERMINANT(gxxh,gxyh,gxzh,gyyh,gyzh,gzzh) + sdet = sqrt(avg_det) call UpperMetric(uxxh, uxyh, uxzh, uyyh, uyzh, uzzh, & avg_det,gxxh, gxyh, gxzh, & @@ -185,11 +194,11 @@ subroutine GRHydro_HLLEM(CCTK_ARGUMENTS) vztm = prim_m(4)-avg_betaz/avg_alp call calc_vlow_blow(gxxh,gxyh,gxzh,gyyh,gyzh,gzzh, & - prim_p(2),prim_p(3),prim_p(4),cons_p(6),cons_p(7),cons_p(8), & + prim_p(2),prim_p(3),prim_p(4),prim_p(8),prim_p(9),prim_p(10), & velxlowp,velylowp,velzlowp,Bvecxlowp,Bvecylowp,Bveczlowp, & bdotvp,b2p,v2p,wp,bxlowp,bylowp,bzlowp) call calc_vlow_blow(gxxh,gxyh,gxzh,gyyh,gyzh,gzzh, & - prim_m(2),prim_m(3),prim_m(4),cons_m(6),cons_m(7),cons_m(8), & + prim_m(2),prim_m(3),prim_m(4),prim_m(8),prim_m(9),prim_m(10), & velxlowm,velylowm,velzlowm,Bvecxlowm,Bvecylowm,Bveczlowm, & bdotvm,b2m,v2m,wm,bxlowm,bylowm,bzlowm) @@ -223,10 +232,10 @@ subroutine GRHydro_HLLEM(CCTK_ARGUMENTS) vxtp,vytp,vztp,pressstarp,bxlowp,bylowp,bzlowp,ab0p,wp, & avg_det,avg_alp,avg_beta) if(clean_divergence.ne.0) then - f1(6)=f1(6)+uxxh*psidcp - f1(7)=f1(7)+uxyh*psidcp - f1(8)=f1(8)+uxzh*psidcp - psidcf = ch_dc*ch_dc*cons_p(6) + f1(6)=f1(6) + avg_alp*sdet*uxxh*psidcp - cons_p(6)*avg_betax + f1(7)=f1(7) + avg_alp*sdet*uxyh*psidcp - cons_p(6)*avg_betay + f1(8)=f1(8) + avg_alp*sdet*uxzh*psidcp - cons_p(6)*avg_betaz + psidcf = avg_alp*cons_p(6)/sdet-psidcp*avg_betax endif else if (flux_direction == 2) then @@ -236,10 +245,10 @@ subroutine GRHydro_HLLEM(CCTK_ARGUMENTS) vytp,vztp,vxtp,pressstarp,bylowp,bzlowp,bxlowp,ab0p,wp, & avg_det,avg_alp,avg_beta) if(clean_divergence.ne.0) then - f1(6)=f1(6)+uxyh*psidcp - f1(7)=f1(7)+uyyh*psidcp - f1(8)=f1(8)+uyzh*psidcp - psidcf = ch_dc*ch_dc*cons_p(7) + f1(6)=f1(6) + avg_alp*sdet*uxyh*psidcp - cons_p(7)*avg_betax + f1(7)=f1(7) + avg_alp*sdet*uyyh*psidcp - cons_p(7)*avg_betay + f1(8)=f1(8) + avg_alp*sdet*uyzh*psidcp - cons_p(7)*avg_betaz + psidcf = avg_alp*cons_p(7)/sdet-psidcp*avg_betay endif else if (flux_direction == 3) then @@ -249,10 +258,10 @@ subroutine GRHydro_HLLEM(CCTK_ARGUMENTS) vztp,vxtp,vytp,pressstarp,bzlowp,bxlowp,bylowp,ab0p,wp, & avg_det,avg_alp,avg_beta) if(clean_divergence.ne.0) then - f1(6)=f1(6)+uxzh*psidcp - f1(7)=f1(7)+uyzh*psidcp - f1(8)=f1(8)+uzzh*psidcp - psidcf = ch_dc*ch_dc*cons_p(8) + f1(6)=f1(6) + avg_alp*sdet*uxzh*psidcp - cons_p(8)*avg_betax + f1(7)=f1(7) + avg_alp*sdet*uyzh*psidcp - cons_p(8)*avg_betay + f1(8)=f1(8) + avg_alp*sdet*uzzh*psidcp - cons_p(8)*avg_betaz + psidcf = avg_alp*cons_p(8)/sdet-psidcp*avg_betaz endif else @@ -291,12 +300,12 @@ subroutine GRHydro_HLLEM(CCTK_ARGUMENTS) if (flux_direction == 1) then call eigenvaluesM(GRHydro_eos_handle,& prim_m(1),prim_m(2),prim_m(3),prim_m(4),prim_m(5),prim_m(7), & - cons_m(6),cons_m(7),cons_m(8),& + prim_m(8),prim_m(9),prim_m(10),& lamminus,gxxh,gxyh,gxzh,gyyh,gyzh,gzzh,& usendh,avg_alp,avg_beta) call eigenvaluesM(GRHydro_eos_handle, & prim_p(1),prim_p(2),prim_p(3),prim_p(4),prim_p(5),prim_p(7), & - cons_p(6),cons_p(7),cons_p(8),& + prim_p(8),prim_p(9),prim_p(10),& lamplus,gxxh,gxyh,gxzh,gyyh,gyzh,gzzh,& usendh,avg_alp,avg_beta) call num_x_fluxM(cons_p(1),cons_p(2),cons_p(3),cons_p(4),cons_p(5),& @@ -312,25 +321,25 @@ subroutine GRHydro_HLLEM(CCTK_ARGUMENTS) avg_det,avg_alp,avg_beta) if(clean_divergence.ne.0) then - fminus(6)=fminus(6)+uxxh*psidcm - fminus(7)=fminus(7)+uxyh*psidcm - fminus(8)=fminus(8)+uxzh*psidcm - fplus(6)=fplus(6)+uxxh*psidcp - fplus(7)=fplus(7)+uxyh*psidcp - fplus(8)=fplus(8)+uxzh*psidcp - psidcfp = ch_dc*ch_dc*cons_p(6) - psidcfm = ch_dc*ch_dc*cons_m(6) + fminus(6)=fminus(6) + avg_alp*sdet*uxxh*psidcm - cons_m(6)*avg_betax + fminus(7)=fminus(7) + avg_alp*sdet*uxyh*psidcm - cons_m(6)*avg_betay + fminus(8)=fminus(8) + avg_alp*sdet*uxzh*psidcm - cons_m(6)*avg_betaz + fplus(6)=fplus(6) + avg_alp*sdet*uxxh*psidcp - cons_p(6)*avg_betax + fplus(7)=fplus(7) + avg_alp*sdet*uxyh*psidcp - cons_p(6)*avg_betay + fplus(8)=fplus(8) + avg_alp*sdet*uxzh*psidcp - cons_p(6)*avg_betaz + psidcfp = avg_alp*cons_p(6)/sdet-avg_betax*psidcp + psidcfm = avg_alp*cons_m(6)/sdet-avg_betax*psidcm endif else if (flux_direction == 2) then call eigenvaluesM(GRHydro_eos_handle,& prim_m(1),prim_m(3),prim_m(4),prim_m(2),prim_m(5),prim_m(7), & - cons_m(7),cons_m(8),cons_m(6),& + prim_m(9),prim_m(10),prim_m(8),& lamminus,gyyh,gyzh,gxyh,gzzh,gxzh,gxxh,& usendh,avg_alp,avg_beta) call eigenvaluesM(GRHydro_eos_handle, & prim_p(1),prim_p(3),prim_p(4),prim_p(2),prim_p(5),prim_p(7), & - cons_p(7),cons_p(8),cons_p(6),& + prim_p(9),prim_p(10),prim_p(8),& lamplus,gyyh,gyzh,gxyh,gzzh,gxzh,gxxh,& usendh,avg_alp,avg_beta) call num_x_fluxM(cons_p(1),cons_p(3),cons_p(4),cons_p(2),cons_p(5),& @@ -345,25 +354,25 @@ subroutine GRHydro_HLLEM(CCTK_ARGUMENTS) vytm,vztm,vxtm,pressstarm,bylowm,bzlowm,bxlowm,ab0m,wm, & avg_det,avg_alp,avg_beta) if(clean_divergence.ne.0) then - fminus(6)=fminus(6)+uxyh*psidcm - fminus(7)=fminus(7)+uyyh*psidcm - fminus(8)=fminus(8)+uyzh*psidcm - fplus(6)=fplus(6)+uxyh*psidcp - fplus(7)=fplus(7)+uyyh*psidcp - fplus(8)=fplus(8)+uyzh*psidcp - psidcfp = ch_dc*ch_dc*cons_p(7) - psidcfm = ch_dc*ch_dc*cons_m(7) + fminus(6)=fminus(6) + avg_alp*sdet*uxyh*psidcm - cons_m(7)*avg_betax + fminus(7)=fminus(7) + avg_alp*sdet*uyyh*psidcm - cons_m(7)*avg_betay + fminus(8)=fminus(8) + avg_alp*sdet*uyzh*psidcm - cons_m(7)*avg_betaz + fplus(6)=fplus(6) + avg_alp*sdet*uxyh*psidcp - cons_p(7)*avg_betax + fplus(7)=fplus(7) + avg_alp*sdet*uyyh*psidcp - cons_p(7)*avg_betay + fplus(8)=fplus(8) + avg_alp*sdet*uyzh*psidcp - cons_p(7)*avg_betaz + psidcfp = avg_alp*cons_p(7)/sdet-avg_betay*psidcp + psidcfm = avg_alp*cons_m(7)/sdet-avg_betay*psidcm endif else if (flux_direction == 3) then call eigenvaluesM(GRHydro_eos_handle,& prim_m(1),prim_m(4),prim_m(2),prim_m(3),prim_m(5),prim_m(7), & - cons_m(8),cons_m(6),cons_m(7),& + prim_m(10),prim_m(8),prim_m(9),& lamminus,gzzh,gxzh,gyzh,gxxh,gxyh,gyyh,& usendh,avg_alp,avg_beta) call eigenvaluesM(GRHydro_eos_handle,& prim_p(1),prim_p(4),prim_p(2),prim_p(3),prim_p(5),prim_p(7), & - cons_p(8),cons_p(6),cons_p(7),& + prim_p(10),prim_p(8),prim_p(9),& lamplus,gzzh,gxzh,gyzh,gxxh,gxyh,gyyh,& usendh,avg_alp,avg_beta) call num_x_fluxM(cons_p(1),cons_p(4),cons_p(2),cons_p(3),cons_p(5),& @@ -378,20 +387,21 @@ subroutine GRHydro_HLLEM(CCTK_ARGUMENTS) vztm,vxtm,vytm,pressstarm,bzlowm,bxlowm,bylowm,ab0m,wm, & avg_det,avg_alp,avg_beta) if(clean_divergence.ne.0) then - fminus(6)=fminus(6)+uxzh*psidcm - fminus(7)=fminus(7)+uyzh*psidcm - fminus(8)=fminus(8)+uzzh*psidcm - fplus(6)=fplus(6)+uxzh*psidcp - fplus(7)=fplus(7)+uyzh*psidcp - fplus(8)=fplus(8)+uzzh*psidcp - psidcfp = ch_dc*ch_dc*cons_p(8) - psidcfm = ch_dc*ch_dc*cons_m(8) + fminus(6)=fminus(6) + avg_alp*sdet*uxzh*psidcm - cons_m(8)*avg_betax + fminus(7)=fminus(7) + avg_alp*sdet*uyzh*psidcm - cons_m(8)*avg_betay + fminus(8)=fminus(8) + avg_alp*sdet*uzzh*psidcm - cons_m(8)*avg_betaz + fplus(6)=fplus(6) + avg_alp*sdet*uxzh*psidcp - cons_p(8)*avg_betax + fplus(7)=fplus(7) + avg_alp*sdet*uyzh*psidcp - cons_p(8)*avg_betay + fplus(8)=fplus(8) + avg_alp*sdet*uzzh*psidcp - cons_p(8)*avg_betaz + psidcfp = avg_alp*cons_p(8)/sdet-avg_betaz*psidcp + psidcfm = avg_alp*cons_m(8)/sdet-avg_betaz*psidcm endif else call CCTK_WARN(0, "Flux direction not x,y,z") end if - + + !!$ Find minimum and maximum wavespeeds charmin = min(0.d0, lamplus(1), lamplus(2), lamplus(3), & @@ -403,18 +413,23 @@ subroutine GRHydro_HLLEM(CCTK_ARGUMENTS) lamminus(4),lamminus(5)) charpm = charmax - charmin + +!!$ if(i.eq.75.and.j.eq.5.and.k.eq.5)write(6,*)'HLLEM:',fplus(6),fminus(6),fplus(8),fminus(8),charmax,charmin,cons_p(6),cons_m(6) !!$ Calculate flux by standard formula - do m = 1,8 - - qdiff(m) = cons_m(m) - cons_p(m) - - f1(m) = (charmax * fplus(m) - charmin * fminus(m) + & - charmax * charmin * qdiff(m)) / charpm - + do m = 1,8 + + qdiff(m) = cons_m(m) - cons_p(m) + + f1(m) = (charmax * fplus(m) - charmin * fminus(m) + & + charmax * charmin * qdiff(m)) / charpm + end do +!!$ if(i.eq.75.and.j.eq.5.and.k.eq.5)write(6,*)'HLLEM2:',fplus(6),fminus(6),fplus(8),fminus(8),charmax,charmin,qdiff(6),qdiff(8),f1(6),f1(8) + + if(clean_divergence.ne.0) then psidcdiff = psidcm - psidcp @@ -436,9 +451,13 @@ subroutine GRHydro_HLLEM(CCTK_ARGUMENTS) syflux(i, j, k) = f1(3) szflux(i, j, k) = f1(4) tauflux(i, j, k) = f1(5) - Bvecxflux(i, j, k) = f1(6) - Bvecyflux(i, j, k) = f1(7) - Bveczflux(i, j, k) = f1(8) + + Bconsxflux(i, j, k) = f1(6) + Bconsyflux(i, j, k) = f1(7) + Bconszflux(i, j, k) = f1(8) + +!!$ if(i.eq.75.and.j.eq.5.and.k.eq.5)write(6,*)'HLLEM2:',fplus(6),fminus(6),fplus(8),fminus(8),f1(6),f1(8) + if(clean_divergence.ne.0) then psidcflux(i,j,k) = psidcf endif @@ -623,12 +642,12 @@ subroutine GRHydro_HLLE_TracerM(CCTK_ARGUMENTS) if (flux_direction == 1) then call eigenvaluesM(GRHydro_eos_handle,& prim_m(1),prim_m(2),prim_m(3),prim_m(4),prim_m(5),prim_m(7), & - cons_m(6),cons_m(7),cons_m(8),& + mag_m(1),mag_m(2),mag_m(3),& lamminus,gxxh,gxyh,gxzh,gyyh,gyzh,gzzh,& usendh,avg_alp,avg_beta) call eigenvaluesM(GRHydro_eos_handle, & prim_p(1),prim_p(2),prim_p(3),prim_p(4),prim_p(5),prim_p(7), & - cons_p(6),cons_p(7),cons_p(8),& + mag_p(1),mag_p(2),mag_p(3),& lamplus,gxxh,gxyh,gxzh,gyyh,gyzh,gzzh,& usendh,avg_alp,avg_beta) fplus(:) = (velxplus(i,j,k) - avg_beta / avg_alp) * & @@ -638,12 +657,12 @@ subroutine GRHydro_HLLE_TracerM(CCTK_ARGUMENTS) else if (flux_direction == 2) then call eigenvaluesM(GRHydro_eos_handle,& prim_m(1),prim_m(3),prim_m(4),prim_m(2),prim_m(5),prim_m(7), & - cons_m(7),cons_m(8),cons_m(6),& + mag_m(2),mag_m(3),mag_m(1),& lamminus,gyyh,gyzh,gxyh,gzzh,gxzh,gxxh,& usendh,avg_alp,avg_beta) call eigenvaluesM(GRHydro_eos_handle, & prim_p(1),prim_p(3),prim_p(4),prim_p(2),prim_p(5),prim_p(7), & - cons_p(7),cons_p(8),cons_p(6),& + mag_p(2),mag_p(3),mag_p(1),& lamplus,gyyh,gyzh,gxyh,gzzh,gxzh,gxxh,& usendh,avg_alp,avg_beta) fplus(:) = (velyplus(i,j,k) - avg_beta / avg_alp) * & @@ -653,12 +672,12 @@ subroutine GRHydro_HLLE_TracerM(CCTK_ARGUMENTS) else if (flux_direction == 3) then call eigenvaluesM(GRHydro_eos_handle,& prim_m(1),prim_m(4),prim_m(2),prim_m(3),prim_m(5),prim_m(7), & - cons_m(8),cons_m(6),cons_m(7),& + mag_m(3),mag_m(1),mag_m(2),& lamminus,gzzh,gxzh,gyzh,gxxh,gxyh,gyyh,& usendh,avg_alp,avg_beta) call eigenvaluesM(GRHydro_eos_handle,& prim_p(1),prim_p(4),prim_p(2),prim_p(3),prim_p(5),prim_p(7), & - cons_p(8),cons_p(6),cons_p(7),& + mag_p(3),mag_p(1),mag_p(2),& lamplus,gzzh,gxzh,gyzh,gxxh,gxyh,gyyh,& usendh,avg_alp,avg_beta) fplus(:) = (velzplus(i,j,k) - avg_beta / avg_alp) * & diff --git a/src/GRHydro_InterfacesM.h b/src/GRHydro_InterfacesM.h index c6b49af..8cd85fa 100644 --- a/src/GRHydro_InterfacesM.h +++ b/src/GRHydro_InterfacesM.h @@ -7,17 +7,18 @@ module Con2PrimM_fortran_interfaces local_gam, dens, & sx, sy, sz, & tau, & + Bconsx, Bconsy, Bconsz, & rho, & velx, vely, velz,& epsilon, pressure,& + Bx, By, Bz, & + bsq,& w_lorentz, & gxx, gxy, gxz, & gyy, gyz, gzz, & uxx, uxy, uxz,& uyy, uyz, uzz,& det,& - Bx, By, Bz, & - bsq,& epsnegative, & retval) @@ -27,17 +28,18 @@ module Con2PrimM_fortran_interfaces CCTK_REAL dens CCTK_REAL sx, sy, sz CCTK_REAL tau + CCTK_REAL Bconsx, Bconsy, Bconsz CCTK_REAL rho CCTK_REAL velx, vely, velz CCTK_REAL epsilon, pressure + CCTK_REAL Bx, By, Bz + CCTK_REAL bsq CCTK_REAL w_lorentz CCTK_REAL gxx, gxy, gxz CCTK_REAL gyy, gyz, gzz CCTK_REAL uxx, uxy, uxz CCTK_REAL uyy, uyz, uzz CCTK_REAL det - CCTK_REAL Bx, By, Bz - CCTK_REAL bsq CCTK_INT epsnegative CCTK_REAL retval end subroutine GRHydro_Con2PrimM_pt @@ -46,17 +48,18 @@ module Con2PrimM_fortran_interfaces dens, & sx, sy, sz, & sc, & + Bconsx, Bconsy, Bconsz, & rho, & velx, vely, velz,& epsilon, pressure,& + Bx, By, Bz, & + bsq,& w_lorentz, & gxx, gxy, gxz, & gyy, gyz, gzz, & uxx, uxy, uxz,& uyy, uyz, uzz,& det,& - Bx, By, Bz, & - bsq,& epsnegative, & retval) @@ -66,17 +69,18 @@ module Con2PrimM_fortran_interfaces CCTK_REAL dens CCTK_REAL sx, sy, sz CCTK_REAL sc + CCTK_REAL Bconsx, Bconsy, Bconsz CCTK_REAL rho CCTK_REAL velx, vely, velz CCTK_REAL epsilon, pressure + CCTK_REAL Bx, By, Bz + CCTK_REAL bsq CCTK_REAL w_lorentz CCTK_REAL gxx, gxy, gxz CCTK_REAL gyy, gyz, gzz CCTK_REAL uxx, uxy, uxz CCTK_REAL uyy, uyz, uzz CCTK_REAL det - CCTK_REAL Bx, By, Bz - CCTK_REAL bsq CCTK_INT epsnegative CCTK_REAL retval end subroutine GRHydro_Con2PrimM_Polytype_pt diff --git a/src/GRHydro_Prim2ConM.F90 b/src/GRHydro_Prim2ConM.F90 index e50c426..b8d01cf 100644 --- a/src/GRHydro_Prim2ConM.F90 +++ b/src/GRHydro_Prim2ConM.F90 @@ -22,6 +22,9 @@ #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) +#define Bconsx(i,j,k) Bcons(i,j,k,1) +#define Bconsy(i,j,k) Bcons(i,j,k,2) +#define Bconsz(i,j,k) Bcons(i,j,k,3) #define DOT(x1,y1,z1,x2,y2,z2) ( DOTP(gxx,gxy,gxz,gyy,gyz,gzz,x1,y1,z1,x2,y2,z2) ) #define DOT2(x1,y1,z1) ( DOTP2(gxx,gxy,gxz,gyy,gyz,gzz,x1,y1,z1) ) @@ -77,16 +80,18 @@ subroutine primitive2conservativeM(CCTK_ARGUMENTS) call prim2conM(GRHydro_eos_handle, gxxl,gxyl,gxzl,gyyl,gyzl,gzzl, & avg_detl,densminus(i,j,k),sxminus(i,j,k),& syminus(i,j,k),szminus(i,j,k),tauminus(i,j,k),& - Bvecxminus(i,j,k),Bvecyminus(i,j,k),Bveczminus(i,j,k), rhominus(i,j,k), & + Bconsxminus(i,j,k),Bconsyminus(i,j,k),Bconszminus(i,j,k), rhominus(i,j,k), & velxminus(i,j,k),velyminus(i,j,k),velzminus(i,j,k),& - epsminus(i,j,k),pressminus(i,j,k),w_lorentzminus(i, j, k)) + epsminus(i,j,k),pressminus(i,j,k),Bvecxminus(i,j,k), & + Bvecyminus(i,j,k), Bveczminus(i,j,k), w_lorentzminus(i, j, k)) call prim2conM(GRHydro_eos_handle, gxxr,gxyr,gxzr,gyyr,gyzr,gzzr, & avg_detr, densplus(i,j,k),sxplus(i,j,k),& syplus(i,j,k),szplus(i,j ,k),tauplus(i,j,k),& - Bvecxplus(i,j,k),Bvecyplus(i,j,k),Bveczplus(i,j,k), & + Bconsxplus(i,j,k),Bconsyplus(i,j,k),Bconszplus(i,j,k), & rhoplus(i,j,k),velxplus(i,j,k),velyplus(i,j,k),& velzplus(i,j,k),epsplus(i,j,k),pressplus(i,j,k),& + Bvecxminus(i,j,k), Bvecyminus(i,j,k), Bveczminus(i,j,k), & w_lorentzplus(i,j,k)) end do @@ -111,7 +116,7 @@ end subroutine primitive2conservativeM @@*/ subroutine prim2conM(handle, gxx, gxy, gxz, gyy, gyz, gzz, det, ddens, & - dsx, dsy, dsz, dtau , dBvcx, dBvcy, dBvcz, drho, dvelx, dvely, dvelz, deps, dpress, w) + dsx, dsy, dsz, dtau , dBconsx, dBconsy, dBconsz, drho, dvelx, dvely, dvelz, deps, dpress, dBvcx, dBvcy, dBvcz, w) implicit none @@ -119,9 +124,10 @@ subroutine prim2conM(handle, gxx, gxy, gxz, gyy, gyz, gzz, det, ddens, & DECLARE_CCTK_FUNCTIONS CCTK_REAL :: gxx, gxy, gxz, gyy, gyz, gzz, det - CCTK_REAL, dimension(1) :: ddens, dsx, dsy, dsz, dtau, dBvcx, dBvcy, dBvcz, & + CCTK_REAL, dimension(1) :: ddens, dsx, dsy, dsz, dtau, & + dBconsx, dBconsy, dBconsz, & drho, dvelx, dvely, dvelz,& - deps, dpress, vlowx, vlowy, vlowz + deps, dpress, dBvcx, dBvcy, dBvcz, vlowx, vlowy, vlowz CCTK_REAL :: w CCTK_REAL, dimension(1) :: Bdotv,ab0,b2,blowx,blowy,blowz CCTK_INT :: handle @@ -171,6 +177,10 @@ subroutine prim2conM(handle, gxx, gxy, gxz, gyy, gyz, gzz, det, ddens, & ab0*blowz) dtau = sqrt(det) * ((drho*(1+deps)+dpress+b2)*w*w - dpress-b2/2.0-ab0**2) - ddens + dBconsx = sqrt(det)*dBvcx + dBconsy = sqrt(det)*dBvcy + dBconsz = sqrt(det)*dBvcz + end subroutine prim2conM @@ -212,9 +222,10 @@ subroutine Primitive2ConservativeCellsM(CCTK_ARGUMENTS) gxy(i,j,k),gxz(i,j,k),& gyy(i,j,k),gyz(i,j,k),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),& + tau(i,j,k),Bconsx(i,j,k),Bconsy(i,j,k),Bconsz(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)) + eps(i,j,k),press(i,j,k),Bvecx(i,j,k), & + Bvecy(i,j,k), Bvecz(i,j,k), w_lorentz(i,j,k)) end do end do @@ -290,19 +301,21 @@ subroutine Prim2ConservativePolytypeM(CCTK_ARGUMENTS) gyyl,gyzl,gzzl, & avg_detl,densminus(i,j,k),sxminus(i,j,k),& syminus(i,j,k),szminus(i,j,k),tauminus(i,j,k),& - Bvecxminus(i,j,k),Bvecyminus(i,j,k),Bveczminus(i,j,k),& + Bconsxminus(i,j,k),Bconsyminus(i,j,k),Bconszminus(i,j,k),& rhominus(i,j,k),& velxminus(i,j,k),velyminus(i,j,k),velzminus(i,j,k),& - epsminus(i,j,k),pressminus(i,j,k),w_lorentzminus(i, j, k)) + epsminus(i,j,k),pressminus(i,j,k),Bvecxminus(i,j,k), & + Bvecyminus(i,j,k),Bveczminus(i,j,k),w_lorentzminus(i, j, k)) call prim2conpolytypeM(GRHydro_eos_handle, gxxr,gxyr,gxzr,& gyyr,gyzr,gzzr, & avg_detr,densplus(i,j,k),sxplus(i,j,k),& syplus(i,j,k),szplus(i,j,k),tauplus(i,j,k),& - Bvecxplus(i,j,k),Bvecyplus(i,j,k),Bveczplus(i,j,k),& + Bconsxplus(i,j,k),Bconsyplus(i,j,k),Bconszplus(i,j,k),& rhoplus(i,j,k),& velxplus(i,j,k),velyplus(i,j,k),velzplus(i,j,k),& - epsplus(i,j,k),pressplus(i,j,k),w_lorentzplus(i, j, k)) + epsplus(i,j,k),pressplus(i,j,k),Bvecxplus(i,j,k), & + Bvecyplus(i,j,k),Bveczplus(i,j,k),w_lorentzplus(i, j, k)) if (densminus(i,j,k) .ne. densminus(i,j,k)) then !$OMP CRITICAL @@ -331,17 +344,17 @@ end subroutine Prim2ConservativePolytypeM @@*/ subroutine prim2conpolytypeM(handle, gxx, gxy, gxz, gyy, gyz, & - gzz, det, ddens, dsx, dsy, dsz, dtau, dBvcx, dBvcy, dBvcz, & - drho, dvelx, dvely, dvelz, deps, dpress, w) + gzz, det, ddens, dsx, dsy, dsz, dtau, dBconsx, dBconsy, dBconsz, & + drho, dvelx, dvely, dvelz, deps, dpress, dBvcx, dBvcy, dBvcz, w) implicit none DECLARE_CCTK_PARAMETERS CCTK_REAL :: gxx, gxy, gxz, gyy, gyz, gzz, det - CCTK_REAL :: ddens, dsx, dsy, dsz, dtau, dBvcx, dBvcy, dBvcz, & - drho, dvelx, dvely, dvelz,& - deps, dpress, w_tmp, w, vlowx, vlowy, vlowz, sqrtdet + CCTK_REAL :: ddens, dsx, dsy, dsz, dtau, dBconsx, dBconsy, dBconsz, & + drho, dvelx, dvely, dvelz, deps, dpress, dBvcx, dBvcy, dBvcz, & + w_tmp, w, vlowx, vlowy, vlowz, sqrtdet CCTK_INT :: handle CCTK_REAL :: Bdotv,ab0,b2,blowx,blowy,blowz character(len=256) NaN_WarnLine @@ -399,6 +412,10 @@ subroutine prim2conpolytypeM(handle, gxx, gxy, gxz, gyy, gyz, & dsz = sqrt(det) * ((drho*(1+deps)+dpress+b2)*w*w * vlowz - & ab0*blowz) dtau = sqrt(det) * ((drho*(1+deps)+dpress+b2)*w*w - dpress-b2/2.0-ab0**2) - ddens + + dBconsx = sqrt(det)*dBvcx + dBconsy = sqrt(det)*dBvcy + dBconsz = sqrt(det)*dBvcz end subroutine prim2conpolytypeM @@ -440,9 +457,10 @@ subroutine Primitive2ConservativePolyCellsM(CCTK_ARGUMENTS) call prim2conpolytypeM(GRHydro_eos_handle,gxx(i,j,k),gxy(i,j,k),& gxz(i,j,k),gyy(i,j,k),gyz(i,j,k),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),& + tau(i,j,k),Bconsx(i,j,k),Bconsy(i,j,k),Bconsz(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)) + eps(i,j,k),press(i,j,k),Bvecx(i,j,k), Bvecy(i,j,k), & + Bvecz(i,j,k), w_lorentz(i,j,k)) end do end do diff --git a/src/GRHydro_Reconstruct.F90 b/src/GRHydro_Reconstruct.F90 index 065d2bd..cf84024 100644 --- a/src/GRHydro_Reconstruct.F90 +++ b/src/GRHydro_Reconstruct.F90 @@ -23,6 +23,9 @@ #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) +#define Bconsx(i,j,k) Bcons(i,j,k,1) +#define Bconsy(i,j,k) Bcons(i,j,k,2) +#define Bconsz(i,j,k) Bcons(i,j,k,3) /*@@ @@ -208,6 +211,17 @@ subroutine Reconstruction(CCTK_ARGUMENTS) vel(:,:,:,3), velzplus, velzminus, trivial_rp, hydro_excision_mask) call tvdreconstruct(nx, ny, nz, xoffset, yoffset, zoffset, & eps, epsplus, epsminus, trivial_rp, hydro_excision_mask) + 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) + endif + +!!$ write(6,*)'recon0:',Bvecxplus(75,5,5),Bvecxminus(75,5,5) + else if (CCTK_EQUALS(recon_vars,"conservative")) then call tvdreconstruct(nx, ny, nz, xoffset, yoffset, zoffset, & dens, densplus, densminus, trivial_rp, hydro_excision_mask) @@ -219,22 +233,22 @@ subroutine Reconstruction(CCTK_ARGUMENTS) scon(:,:,:,3), szplus, szminus, trivial_rp, hydro_excision_mask) call tvdreconstruct(nx, ny, nz, xoffset, yoffset, zoffset, & tau, tauplus, tauminus, trivial_rp, hydro_excision_mask) + if(evolve_mhd.ne.0) then + call tvdreconstruct(nx, ny, nz, xoffset, yoffset, zoffset, & + Bcons(:,:,:,1), Bconsxplus, Bconsxminus, trivial_rp, hydro_excision_mask) + call tvdreconstruct(nx, ny, nz, xoffset, yoffset, zoffset, & + Bcons(:,:,:,2), Bconsyplus, Bconsyminus, trivial_rp, hydro_excision_mask) + call tvdreconstruct(nx, ny, nz, xoffset, yoffset, zoffset, & + Bcons(:,:,:,3), Bconszplus, Bconszminus, trivial_rp, hydro_excision_mask) + endif + else 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 + if(evolve_mhd.ne.0.and.clean_divergence.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 + psidc, psidcplus, psidcminus, trivial_rp, hydro_excision_mask) endif !$OMP PARALLEL DO PRIVATE(i, j) @@ -640,6 +654,17 @@ subroutine Reconstruction(CCTK_ARGUMENTS) call GRHydro_ENOReconstruct1d(eno_order,cctk_lsh(1),& eps(:,j,k),epsminus(:,j,k),epsplus(:,j,k),& trivial_rp(:,j,k), hydro_excision_mask(:,j,k)) + 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)) + endif else if (CCTK_EQUALS(recon_vars,"conservative")) then call GRHydro_ENOReconstruct1d(eno_order,cctk_lsh(1),& dens(:,j,k),densminus(:,j,k),densplus(:,j,k),& @@ -656,28 +681,27 @@ subroutine Reconstruction(CCTK_ARGUMENTS) call GRHydro_ENOReconstruct1d(eno_order,cctk_lsh(1),& tau(:,j,k),tauminus(:,j,k),tauplus(:,j,k),& trivial_rp(:,j,k), hydro_excision_mask(:,j,k)) + if(evolve_mhd.ne.0) then + call GRHydro_ENOReconstruct1d(eno_order,cctk_lsh(1),& + Bconsx(:,j,k),Bconsxminus(:,j,k),Bconsxplus(:,j,k),& + trivial_rp(:,j,k), hydro_excision_mask(:,j,k)) + call GRHydro_ENOReconstruct1d(eno_order,cctk_lsh(1),& + Bconsy(:,j,k),Bconsyminus(:,j,k),Bconsyplus(:,j,k),& + trivial_rp(:,j,k), hydro_excision_mask(:,j,k)) + call GRHydro_ENOReconstruct1d(eno_order,cctk_lsh(1),& + Bconsz(:,j,k),Bconszminus(:,j,k),Bconszplus(:,j,k),& + trivial_rp(:,j,k), hydro_excision_mask(:,j,k)) + endif else !$OMP CRITICAL 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)) + if(evolve_mhd.ne.0.and.clean_divergence.ne.0) then call GRHydro_ENOReconstruct1d(eno_order,cctk_lsh(1),& - Bvecz(:,j,k),Bveczminus(:,j,k),Bveczplus(:,j,k),& + psidc(:,j,k),psidcminus(:,j,k),psidcplus(:,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) @@ -710,7 +734,18 @@ subroutine Reconstruction(CCTK_ARGUMENTS) call GRHydro_ENOReconstruct1d(eno_order,cctk_lsh(2),& eps(j,:,k),epsminus(j,:,k),epsplus(j,:,k),& trivial_rp(j,:,k), hydro_excision_mask(j,:,k)) - else if (CCTK_EQUALS(recon_vars,"conservative")) then + 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)) + endif + else if (CCTK_EQUALS(recon_vars,"conservative")) then call GRHydro_ENOReconstruct1d(eno_order,cctk_lsh(2),& dens(j,:,k),densminus(j,:,k),densplus(j,:,k),& trivial_rp(j,:,k), hydro_excision_mask(j,:,k)) @@ -726,27 +761,27 @@ subroutine Reconstruction(CCTK_ARGUMENTS) call GRHydro_ENOReconstruct1d(eno_order,cctk_lsh(2),& tau(j,:,k),tauminus(j,:,k),tauplus(j,:,k),& trivial_rp(j,:,k), hydro_excision_mask(j,:,k)) + if(evolve_mhd.ne.0) then + call GRHydro_ENOReconstruct1d(eno_order,cctk_lsh(2),& + Bconsx(j,:,k),Bconsxminus(j,:,k),Bconsxplus(j,:,k),& + trivial_rp(j,:,k), hydro_excision_mask(j,:,k)) + call GRHydro_ENOReconstruct1d(eno_order,cctk_lsh(2),& + Bconsy(j,:,k),Bconsyminus(j,:,k),Bconsyplus(j,:,k),& + trivial_rp(j,:,k), hydro_excision_mask(j,:,k)) + call GRHydro_ENOReconstruct1d(eno_order,cctk_lsh(2),& + Bconsz(j,:,k),Bconszminus(j,:,k),Bconszplus(j,:,k),& + trivial_rp(j,:,k), hydro_excision_mask(j,:,k)) + endif else !$OMP CRITICAL 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 + if(evolve_mhd.ne.0.and.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 do i = 1, cctk_lsh(2) @@ -779,6 +814,17 @@ subroutine Reconstruction(CCTK_ARGUMENTS) call GRHydro_ENOReconstruct1d(eno_order,cctk_lsh(3),& eps(j,k,:),epsminus(j,k,:),epsplus(j,k,:),& trivial_rp(j,k,:), hydro_excision_mask(j,k,:)) + 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,:)) + endif else if (CCTK_EQUALS(recon_vars,"conservative")) then call GRHydro_ENOReconstruct1d(eno_order,cctk_lsh(3),& dens(j,k,:),densminus(j,k,:),densplus(j,k,:),& @@ -795,27 +841,27 @@ subroutine Reconstruction(CCTK_ARGUMENTS) call GRHydro_ENOReconstruct1d(eno_order,cctk_lsh(3),& tau(j,k,:),tauminus(j,k,:),tauplus(j,k,:),& trivial_rp(j,k,:), hydro_excision_mask(j,k,:)) + if(evolve_mhd.ne.0) then + call GRHydro_ENOReconstruct1d(eno_order,cctk_lsh(3),& + Bconsx(j,k,:),Bconsxminus(j,k,:),Bconsxplus(j,k,:),& + trivial_rp(j,k,:), hydro_excision_mask(j,k,:)) + call GRHydro_ENOReconstruct1d(eno_order,cctk_lsh(3),& + Bconsy(j,k,:),Bconsyminus(j,k,:),Bconsyplus(j,k,:),& + trivial_rp(j,k,:), hydro_excision_mask(j,k,:)) + call GRHydro_ENOReconstruct1d(eno_order,cctk_lsh(3),& + Bconsz(j,k,:),Bconszminus(j,k,:),Bconszplus(j,k,:),& + trivial_rp(j,k,:), hydro_excision_mask(j,k,:)) + endif else !$OMP CRITICAL 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,:)) + if(evolve_mhd.ne.0.and.clean_divergence.ne.0) then call GRHydro_ENOReconstruct1d(eno_order,cctk_lsh(3),& - Bvecz(j,k,:),Bveczminus(j,k,:),Bveczplus(j,k,:),& + psidc(j,k,:),psidcminus(j,k,:),psidcplus(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) @@ -877,6 +923,9 @@ subroutine Reconstruction(CCTK_ARGUMENTS) CCTK_EQUALS(recon_method,"ppm")) then if(evolve_mhd.ne.0) then call primitive2conservativeM(CCTK_PASS_FTOF) + +!!$ write(6,*)'recon1:',Bvecxplus(75,5,5),Bvecxminus(75,5,5),Bconsxplus(75,5,5),Bconsxminus(75,5,5) + else call primitive2conservative(CCTK_PASS_FTOF) endif diff --git a/src/GRHydro_ReconstructPoly.F90 b/src/GRHydro_ReconstructPoly.F90 index db78bd9..826ef1c 100644 --- a/src/GRHydro_ReconstructPoly.F90 +++ b/src/GRHydro_ReconstructPoly.F90 @@ -23,6 +23,9 @@ #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) +#define Bconsx(i,j,k) Bcons(i,j,k,1) +#define Bconsy(i,j,k) Bcons(i,j,k,2) +#define Bconsz(i,j,k) Bcons(i,j,k,3) /*@@ @@ -202,7 +205,14 @@ subroutine ReconstructionPolytype(CCTK_ARGUMENTS) vel(:,:,:,2), velyplus, velyminus, trivial_rp, hydro_excision_mask) call tvdreconstruct(nx, ny, nz, xoffset, yoffset, zoffset, & vel(:,:,:,3), velzplus, velzminus, trivial_rp, hydro_excision_mask) - + 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) + endif else if (CCTK_EQUALS(recon_vars,"conservative")) then call tvdreconstruct(nx, ny, nz, xoffset, yoffset, zoffset, & dens, densplus, densminus, trivial_rp, hydro_excision_mask) @@ -212,23 +222,21 @@ subroutine ReconstructionPolytype(CCTK_ARGUMENTS) scon(:,:,:,2), syplus, syminus, trivial_rp, hydro_excision_mask) call tvdreconstruct(nx, ny, nz, xoffset, yoffset, zoffset, & scon(:,:,:,3), szplus, szminus, trivial_rp, hydro_excision_mask) - + if(evolve_mhd.ne.0) then + call tvdreconstruct(nx, ny, nz, xoffset, yoffset, zoffset, & + Bcons(:,:,:,1), Bconsxplus, Bconsxminus, trivial_rp, hydro_excision_mask) + call tvdreconstruct(nx, ny, nz, xoffset, yoffset, zoffset, & + Bcons(:,:,:,2), Bconsyplus, Bconsyminus, trivial_rp, hydro_excision_mask) + call tvdreconstruct(nx, ny, nz, xoffset, yoffset, zoffset, & + Bcons(:,:,:,3), Bconszplus, Bconszminus, trivial_rp, hydro_excision_mask) + endif else 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) + if(evolve_mhd.ne.0.and.clean_divergence.ne.0) then 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 + psidc, psidcplus, psidcminus, trivial_rp, hydro_excision_mask) endif !$OMP PARALLEL DO PRIVATE(i, j) @@ -630,6 +638,17 @@ subroutine ReconstructionPolytype(CCTK_ARGUMENTS) call GRHydro_ENOReconstruct1d(eno_order,cctk_lsh(1),& velz(:,j,k),velzminus(:,j,k),velzplus(:,j,k),& trivial_rp(:,j,k), hydro_excision_mask(:,j,k)) + 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)) + endif else if (CCTK_EQUALS(recon_vars,"conservative")) then call GRHydro_ENOReconstruct1d(eno_order,cctk_lsh(1),& dens(:,j,k),densminus(:,j,k),densplus(:,j,k),& @@ -643,28 +662,27 @@ subroutine ReconstructionPolytype(CCTK_ARGUMENTS) call GRHydro_ENOReconstruct1d(eno_order,cctk_lsh(1),& sz(:,j,k),szminus(:,j,k),szplus(:,j,k),& trivial_rp(:,j,k), hydro_excision_mask(:,j,k)) + if(evolve_mhd.ne.0) then + call GRHydro_ENOReconstruct1d(eno_order,cctk_lsh(1),& + Bconsx(:,j,k),Bconsxminus(:,j,k),Bconsxplus(:,j,k),& + trivial_rp(:,j,k), hydro_excision_mask(:,j,k)) + call GRHydro_ENOReconstruct1d(eno_order,cctk_lsh(1),& + Bconsy(:,j,k),Bconsyminus(:,j,k),Bconsyplus(:,j,k),& + trivial_rp(:,j,k), hydro_excision_mask(:,j,k)) + call GRHydro_ENOReconstruct1d(eno_order,cctk_lsh(1),& + Bconsz(:,j,k),Bconszminus(:,j,k),Bconszplus(:,j,k),& + trivial_rp(:,j,k), hydro_excision_mask(:,j,k)) + endif else !$OMP CRITICAL 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)) + if(evolve_mhd.ne.0.and.clean_divergence.ne.0) then call GRHydro_ENOReconstruct1d(eno_order,cctk_lsh(1),& - Bvecz(:,j,k),Bveczminus(:,j,k),Bveczplus(:,j,k),& + psidc(:,j,k),psidcminus(:,j,k),psidcplus(:,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) @@ -675,70 +693,78 @@ subroutine ReconstructionPolytype(CCTK_ARGUMENTS) end if end do end do - end do - !$OMP END PARALLEL DO - else if (flux_direction == 2) then - !$OMP PARALLEL DO PRIVATE(i, j) - do k = GRHydro_stencil, cctk_lsh(3) - GRHydro_stencil + 1 - do j = GRHydro_stencil, cctk_lsh(1) - GRHydro_stencil + 1 - if (CCTK_EQUALS(recon_vars,"primitive")) then - call GRHydro_ENOReconstruct1d(eno_order,cctk_lsh(2),& - rho(j,:,k),rhominus(j,:,k),rhoplus(j,:,k),& - trivial_rp(j,:,k), hydro_excision_mask(j,:,k)) - call GRHydro_ENOReconstruct1d(eno_order,cctk_lsh(2),& - velx(j,:,k),velxminus(j,:,k),velxplus(j,:,k),& - trivial_rp(j,:,k), hydro_excision_mask(j,:,k)) - call GRHydro_ENOReconstruct1d(eno_order,cctk_lsh(2),& - vely(j,:,k),velyminus(j,:,k),velyplus(j,:,k),& - trivial_rp(j,:,k), hydro_excision_mask(j,:,k)) - call GRHydro_ENOReconstruct1d(eno_order,cctk_lsh(2),& - velz(j,:,k),velzminus(j,:,k),velzplus(j,:,k),& - trivial_rp(j,:,k), hydro_excision_mask(j,:,k)) - else if (CCTK_EQUALS(recon_vars,"conservative")) then - call GRHydro_ENOReconstruct1d(eno_order,cctk_lsh(2),& - dens(j,:,k),densminus(j,:,k),densplus(j,:,k),& - trivial_rp(j,:,k), hydro_excision_mask(j,:,k)) - call GRHydro_ENOReconstruct1d(eno_order,cctk_lsh(2),& - sx(j,:,k),sxminus(j,:,k),sxplus(j,:,k),& - trivial_rp(j,:,k), hydro_excision_mask(j,:,k)) - call GRHydro_ENOReconstruct1d(eno_order,cctk_lsh(2),& - sy(j,:,k),syminus(j,:,k),syplus(j,:,k),& - trivial_rp(j,:,k), hydro_excision_mask(j,:,k)) - call GRHydro_ENOReconstruct1d(eno_order,cctk_lsh(2),& - sz(j,:,k),szminus(j,:,k),szplus(j,:,k),& - trivial_rp(j,:,k), hydro_excision_mask(j,:,k)) - else - !$OMP CRITICAL - call CCTK_WARN(0, "Variable type to reconstruct not recognized.") - !$OMP END CRITICAL - end if + end do + !$OMP END PARALLEL DO + else if (flux_direction == 2) then + !$OMP PARALLEL DO PRIVATE(i, j) + do k = GRHydro_stencil, cctk_lsh(3) - GRHydro_stencil + 1 + do j = GRHydro_stencil, cctk_lsh(1) - GRHydro_stencil + 1 + if (CCTK_EQUALS(recon_vars,"primitive")) then + call GRHydro_ENOReconstruct1d(eno_order,cctk_lsh(2),& + rho(j,:,k),rhominus(j,:,k),rhoplus(j,:,k),& + trivial_rp(j,:,k), hydro_excision_mask(j,:,k)) + call GRHydro_ENOReconstruct1d(eno_order,cctk_lsh(2),& + velx(j,:,k),velxminus(j,:,k),velxplus(j,:,k),& + trivial_rp(j,:,k), hydro_excision_mask(j,:,k)) + call GRHydro_ENOReconstruct1d(eno_order,cctk_lsh(2),& + vely(j,:,k),velyminus(j,:,k),velyplus(j,:,k),& + trivial_rp(j,:,k), hydro_excision_mask(j,:,k)) + call GRHydro_ENOReconstruct1d(eno_order,cctk_lsh(2),& + velz(j,:,k),velzminus(j,:,k),velzplus(j,:,k),& + trivial_rp(j,:,k), hydro_excision_mask(j,:,k)) + 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)) + + else if (CCTK_EQUALS(recon_vars,"conservative")) then + call GRHydro_ENOReconstruct1d(eno_order,cctk_lsh(2),& + dens(j,:,k),densminus(j,:,k),densplus(j,:,k),& + trivial_rp(j,:,k), hydro_excision_mask(j,:,k)) + call GRHydro_ENOReconstruct1d(eno_order,cctk_lsh(2),& + sx(j,:,k),sxminus(j,:,k),sxplus(j,:,k),& + trivial_rp(j,:,k), hydro_excision_mask(j,:,k)) + call GRHydro_ENOReconstruct1d(eno_order,cctk_lsh(2),& + sy(j,:,k),syminus(j,:,k),syplus(j,:,k),& + trivial_rp(j,:,k), hydro_excision_mask(j,:,k)) + call GRHydro_ENOReconstruct1d(eno_order,cctk_lsh(2),& + sz(j,:,k),szminus(j,:,k),szplus(j,:,k),& + trivial_rp(j,:,k), hydro_excision_mask(j,:,k)) + call GRHydro_ENOReconstruct1d(eno_order,cctk_lsh(2),& + Bconsx(j,:,k),Bconsxminus(j,:,k),Bconsxplus(j,:,k),& + trivial_rp(j,:,k), hydro_excision_mask(j,:,k)) + call GRHydro_ENOReconstruct1d(eno_order,cctk_lsh(2),& + Bconsy(j,:,k),Bconsyminus(j,:,k),Bconsyplus(j,:,k),& + trivial_rp(j,:,k), hydro_excision_mask(j,:,k)) + call GRHydro_ENOReconstruct1d(eno_order,cctk_lsh(2),& + Bconsz(j,:,k),Bconszminus(j,:,k),Bconszplus(j,:,k),& + trivial_rp(j,:,k), hydro_excision_mask(j,:,k)) + else + !$OMP CRITICAL + 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) - else - SpaceMask_SetStateBitsF90(space_mask, j, i, k, type_bitsy, not_trivialy) - end if - end do + if(evolve_mhd.ne.0.and.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 + + do i = 1, cctk_lsh(2) + 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 else if (flux_direction == 3) then !$OMP PARALLEL DO PRIVATE(i, j) @@ -757,6 +783,17 @@ subroutine ReconstructionPolytype(CCTK_ARGUMENTS) call GRHydro_ENOReconstruct1d(eno_order,cctk_lsh(3),& velz(j,k,:),velzminus(j,k,:),velzplus(j,k,:),& trivial_rp(j,k,:), hydro_excision_mask(j,k,:)) + 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,:)) + endif else if (CCTK_EQUALS(recon_vars,"conservative")) then call GRHydro_ENOReconstruct1d(eno_order,cctk_lsh(3),& dens(j,k,:),densminus(j,k,:),densplus(j,k,:),& @@ -770,27 +807,27 @@ subroutine ReconstructionPolytype(CCTK_ARGUMENTS) call GRHydro_ENOReconstruct1d(eno_order,cctk_lsh(3),& sz(j,k,:),szminus(j,k,:),szplus(j,k,:),& trivial_rp(j,k,:), hydro_excision_mask(j,k,:)) + if(evolve_mhd.ne.0) then + call GRHydro_ENOReconstruct1d(eno_order,cctk_lsh(3),& + Bconsx(j,k,:),Bconsxminus(j,k,:),Bconsxplus(j,k,:),& + trivial_rp(j,k,:), hydro_excision_mask(j,k,:)) + call GRHydro_ENOReconstruct1d(eno_order,cctk_lsh(3),& + Bconsy(j,k,:),Bconsyminus(j,k,:),Bconsyplus(j,k,:),& + trivial_rp(j,k,:), hydro_excision_mask(j,k,:)) + call GRHydro_ENOReconstruct1d(eno_order,cctk_lsh(3),& + Bconsz(j,k,:),Bconszminus(j,k,:),Bconszplus(j,k,:),& + trivial_rp(j,k,:), hydro_excision_mask(j,k,:)) + endif else !$OMP CRITICAL 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,:)) + if(evolve_mhd.ne.0.and.clean_divergence.ne.0) then call GRHydro_ENOReconstruct1d(eno_order,cctk_lsh(3),& - Bvecz(j,k,:),Bveczminus(j,k,:),Bveczplus(j,k,:),& + psidc(j,k,:),psidcminus(j,k,:),psidcplus(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) diff --git a/src/GRHydro_RegisterVars.cc b/src/GRHydro_RegisterVars.cc index 3ac815f..0d99955 100644 --- a/src/GRHydro_RegisterVars.cc +++ b/src/GRHydro_RegisterVars.cc @@ -63,7 +63,8 @@ extern "C"void GRHydro_Register(CCTK_ARGUMENTS) register_evolved("GRHydro::scon", "GRHydro::srhs"); if (CCTK_EQUALS(Bvec_evolution_method, "GRHydro")) { - register_evolved("HydroBase::Bvec", "GRHydro::Bvecrhs"); + register_constrained("HydroBase::Bvec"); + register_evolved("GRHydro::Bcons", "GRHydro::Bconsrhs"); if(clean_divergence) { register_evolved("GRHydro::psidc" , "GRHydro::psidcrhs"); } @@ -124,7 +125,7 @@ extern "C"void GRHydro_Register(CCTK_ARGUMENTS) register_constrained("GRHydro::tau"); if (CCTK_EQUALS(Bvec_evolution_method, "GRHydro")) { - register_constrained("HydroBase::Bvec"); + register_constrained("HydroBase::Bcons"); if(clean_divergence) { register_constrained("GRHydro::psidc"); } diff --git a/src/GRHydro_SourceM.F90 b/src/GRHydro_SourceM.F90 index 257a24c..ce3bcdf 100644 --- a/src/GRHydro_SourceM.F90 +++ b/src/GRHydro_SourceM.F90 @@ -33,6 +33,12 @@ #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) +#define Bconsx(i,j,k) Bcons(i,j,k,1) +#define Bconsy(i,j,k) Bcons(i,j,k,2) +#define Bconsz(i,j,k) Bcons(i,j,k,3) +#define Bconsrhsx(i,j,k) Bconsrhs(i,j,k,1) +#define Bconsrhsy(i,j,k) Bconsrhs(i,j,k,2) +#define Bconsrhsz(i,j,k) Bconsrhs(i,j,k,3) /*@@ @routine SourceTermsM @@ -76,8 +82,9 @@ subroutine SourceTermsM(CCTK_ARGUMENTS) CCTK_REAL :: halfshiftdgx, halfshiftdgy, halfshiftdgz CCTK_REAL :: halfTdgx, halfTdgy, halfTdgz CCTK_REAL :: invalp, invalp2 - CCTK_REAL :: dx_psidc, dy_psidc, dz_psidc CCTK_REAL :: bvcx_source, bvcy_source, bvcz_source + CCTK_REAL :: dx_det_bydet, dy_det_bydet, dz_det_bydet + CCTK_REAL :: gdg_x, gdg_y, gdg_z !! g^{ik} d_k g_{ij} CCTK_REAL :: Bvecxlow,Bvecylow,Bveczlow,bdotv,b2,dum,bxlow,bylow,bzlow CCTK_REAL :: bt,bx,by,bz,rhohstarW2,pstar @@ -103,7 +110,7 @@ subroutine SourceTermsM(CCTK_ARGUMENTS) densrhs = 0.d0 srhs = 0.d0 taurhs = 0.d0 - Bvecrhs=0.d0 + Bconsrhs=0.d0 if (evolve_tracer .ne. 0) then cons_tracerrhs = 0.d0 @@ -157,7 +164,7 @@ subroutine SourceTermsM(CCTK_ARGUMENTS) !$OMP vlowx,vlowy,vlowz,Bvecxlow,Bvecylow,Bveczlow, & !$OMP bdotv,b2,dum,bxlow,bylow,bzlow,bt,bx,by,bz,rhohstarW2,pstar,& !$OMP t00,t0x,t0y,t0z,txx,txy,txz,tyy,tyz,tzz,& - !$OMP dx_alp,dy_alp,dz_alp,dx_psidc,dy_psidc,dz_psidc,& + !$OMP dx_alp,dy_alp,dz_alp,& !$OMP tau_source,sx_source,sy_source,sz_source,& !$OMP bvcx_source,bvcy_source,bvcz_source,& !$OMP uxx, uxy, uxz, uyy, uyz, uzz,& @@ -435,8 +442,74 @@ subroutine SourceTermsM(CCTK_ARGUMENTS) srhs(i,j,k,2) = alp(i,j,k)*sqrtdet*sy_source srhs(i,j,k,3) = alp(i,j,k)*sqrtdet*sz_source taurhs(i,j,k) = alp(i,j,k)*sqrtdet*tau_source + +!!$ if(i.eq.75.and.j.eq.5.and.k.eq.5)write(6,*)'source:',i,j,k, & +!!$ Bconsrhs(i,j,k,1),Bconsrhs(i,j,k,3) + if(clean_divergence.ne.0) then - psidcrhs(i,j,k) = -1.d0*ch_dc*ch_dc/cp_dc/cp_dc*psidc(i,j,k) + psidcrhs(i,j,k) = -1.d0 * ( ch_dc*ch_dc/cp_dc/cp_dc*alp(i,j,k) + & + dx_betax + dy_betay + dz_betaz ) * psidc(i,j,k) + & + Bconsx(i,j,k) * ( dx_alp - half*alp(i,j,k) * & + ( uxx*dx_gxx + uyy*dx_gyy + uzz*dx_gzz + 2.d0*uxy*dx_gxy + & + 2.d0*uxz*dx_gxz + 2.d0*uyz*dx_gyz ) )/ sqrtdet + & + Bconsy(i,j,k) * (dy_alp - half*alp(i,j,k) * & + ( uxx*dy_gxx + uyy*dy_gyy + uzz*dy_gzz + 2.d0*uxy*dy_gxy + & + 2.d0*uxz*dy_gxz + 2.d0*uyz*dy_gyz ) )/ sqrtdet + & + Bconsz(i,j,k) * (dz_alp - half*alp(i,j,k) * & + ( uxx*dz_gxx + uyy*dy_gyy + uzz*dy_gzz + 2.d0*uxy*dy_gxy + & + 2.d0*uxz*dz_gxz + 2.d0*uyz*dy_gyz ) )/ sqrtdet + + !!$ g^{jk} d_i g_{kj} = d_i (g) / det + dx_det_bydet = uxx*dx_gxx + uyy*dx_gyy + uzz*dx_gzz + & + 2.d0*(uxy*dx_gxy+uxz*dx_gxz+uyz*dx_gyz) + dy_det_bydet = uxx*dy_gxx + uyy*dy_gyy + uzz*dy_gzz + & + 2.d0*(uxy*dy_gxy+uxz*dy_gxz+uyz*dy_gyz) + dz_det_bydet = uxx*dz_gxx + uyy*dz_gyy + uzz*dz_gzz + & + 2.d0*(uxy*dz_gxy+uxz*dz_gxz+uyz*dz_gyz) + + !!$ g^{ik} d_k g_{li} + gdg_x = uxx*dx_gxx + uxy*dy_gxx + uxz*dz_gxx + & + uxy*dx_gxy + uyy*dy_gxy + uyz*dz_gxy + & + uxz*dx_gxz + uyz*dy_gxz + uzz*dz_gxz + + gdg_y = uxx*dx_gxy + uxy*dy_gxy + uxz*dz_gxy + & + uxy*dx_gyy + uyy*dy_gyy + uyz*dz_gyy + & + uxz*dx_gyz + uyz*dy_gyz + uzz*dz_gyz + + gdg_z = uxx*dx_gxz + uxy*dy_gxz + uxz*dz_gxz + & + uxy*dx_gyz + uyy*dy_gyz + uyz*dz_gyz + & + uxz*dx_gzz + uyz*dy_gzz + uzz*dz_gzz + + bvcx_source = -1.d0 * ( Bconsx(i,j,k)*dx_betax + & + Bconsy(i,j,k)*dy_betax + Bconsz(i,j,k)*dz_betax ) + & + psidc(i,j,k)*sqrtdet*( uxx*dx_alp+uxy*dy_alp+uxz*dz_alp ) + & + psidc(i,j,k)*alp(i,j,k)*sqrtdet*( uxx*dx_det_bydet + & + uxy*dy_det_bydet + uxz*dz_det_bydet ) - & + psidc(i,j,k)*alp(i,j,k)*sqrtdet*( uxx*gdg_x + uxy*gdg_y + & + uxz*gdg_z ) + + bvcy_source = -1.d0 * ( Bconsx(i,j,k)*dx_betay + & + Bconsy(i,j,k)*dy_betay + Bconsz(i,j,k)*dz_betay ) + & + psidc(i,j,k)*sqrtdet*( uxy*dx_alp+uyy*dy_alp+uyz*dz_alp ) + & + psidc(i,j,k)*alp(i,j,k)*sqrtdet*( uxy*dx_det_bydet + & + uyy*dy_det_bydet + uyz*dz_det_bydet ) - & + psidc(i,j,k)*alp(i,j,k)*sqrtdet*( uxy*gdg_x + uyy*gdg_y + & + uyz*gdg_z ) + + bvcz_source = -1.d0 * ( Bconsx(i,j,k)*dx_betaz + & + Bconsy(i,j,k)*dy_betaz + Bconsz(i,j,k)*dz_betaz ) + & + psidc(i,j,k)*sqrtdet*( uxz*dx_alp+uyz*dy_alp+uzz*dz_alp ) + & + psidc(i,j,k)*alp(i,j,k)*sqrtdet*( uxz*dx_det_bydet + & + uyz*dy_det_bydet + uzz*dz_det_bydet ) - & + psidc(i,j,k)*alp(i,j,k)*sqrtdet*( uxz*gdg_x + uyz*gdg_y + & + uzz*gdg_z ) + +!!$ if(i.eq.5.and.j.eq.5.and.k.eq.5)write(6,*)'source:',i,j,k, & +!!$ bvcx_source,bvcz_source + + Bconsrhsx(i,j,k) = bvcx_source + Bconsrhsy(i,j,k) = bvcy_source + Bconsrhsz(i,j,k) = bvcz_source endif enddo diff --git a/src/GRHydro_UpdateMaskM.F90 b/src/GRHydro_UpdateMaskM.F90 index 25a0b9a..9b6f65a 100644 --- a/src/GRHydro_UpdateMaskM.F90 +++ b/src/GRHydro_UpdateMaskM.F90 @@ -75,9 +75,10 @@ subroutine GRHydro_AtmosphereResetM(CCTK_ARGUMENTS) gxx(i,j,k), gxy(i,j,k), gxz(i,j,k), & gyy(i,j,k), gyz(i,j,k), gzz(i,j,k), det, & dens(i,j,k), scon(i,j,k,1), scon(i,j,k,2), scon(i,j,k,3), & - tau(i,j,k), Bvec(i,j,k,1),Bvec(i,j,k,2),Bvec(i,j,k,3),& + tau(i,j,k), Bcons(i,j,k,1),Bcons(i,j,k,2),Bcons(i,j,k,3),& rho(i,j,k), vel(i,j,k,1), vel(i,j,k,2), & - vel(i,j,k,3), eps(i,j,k), press(i,j,k), w_lorentz(i,j,k)) + vel(i,j,k,3), eps(i,j,k), press(i,j,k), & + Bvec(i,j,k,1),Bvec(i,j,k,2),Bvec(i,j,k,3),w_lorentz(i,j,k)) if (wk_atmosphere .eq. 0) then atmosphere_mask(i, j, k) = 0 SpaceMask_SetStateBitsF90(space_mask, i, j, k, type_bits,\ @@ -143,9 +144,10 @@ subroutine GRHydro_InitialAtmosphereResetM(CCTK_ARGUMENTS) gxx(i,j,k), gxy(i,j,k), gxz(i,j,k), & gyy(i,j,k), gyz(i,j,k), gzz(i,j,k), det, & dens(i,j,k), scon(i,j,k,1), scon(i,j,k,2), scon(i,j,k,3), & - tau(i,j,k), Bvec(i,j,k,1),Bvec(i,j,k,2),Bvec(i,j,k,3),& + tau(i,j,k), Bcons(i,j,k,1),Bcons(i,j,k,2),Bcons(i,j,k,3),& rho(i,j,k), vel(i,j,k,1), vel(i,j,k,2), & - vel(i,j,k,3), eps(i,j,k), press(i,j,k), w_lorentz(i,j,k)) + vel(i,j,k,3), eps(i,j,k), press(i,j,k), & + Bvec(i,j,k,1),Bvec(i,j,k,2),Bvec(i,j,k,3),w_lorentz(i,j,k)) end if if (timelevels .gt. 1) then if (rho_p(i,j,k) .le. GRHydro_rho_min) then @@ -165,9 +167,10 @@ subroutine GRHydro_InitialAtmosphereResetM(CCTK_ARGUMENTS) gxx_p(i,j,k), gxy_p(i,j,k), gxz_p(i,j,k), & gyy_p(i,j,k), gyz_p(i,j,k), gzz_p(i,j,k), det, & dens_p(i,j,k), scon_p(i,j,k,1), scon_p(i,j,k,2), scon_p(i,j,k,3), & - tau_p(i,j,k), Bvec_p(i,j,k,1),Bvec_p(i,j,k,2),Bvec_p(i,j,k,3),& + tau_p(i,j,k), Bcons_p(i,j,k,1),Bcons_p(i,j,k,2),Bcons_p(i,j,k,3),& rho_p(i,j,k), vel_p(i,j,k,1), vel_p(i,j,k,2), & - vel_p(i,j,k,3), eps_p(i,j,k), press_p(i,j,k), w_lorentz_p(i,j,k)) + vel_p(i,j,k,3), eps_p(i,j,k), press_p(i,j,k), & + Bvec_p(i,j,k,1),Bvec_p(i,j,k,2),Bvec_p(i,j,k,3),w_lorentz_p(i,j,k)) endif end if if (timelevels .gt. 2) then @@ -187,9 +190,10 @@ subroutine GRHydro_InitialAtmosphereResetM(CCTK_ARGUMENTS) gxx_p_p(i,j,k), gxy_p_p(i,j,k), gxz_p_p(i,j,k), & gyy_p_p(i,j,k), gyz_p_p(i,j,k), gzz_p_p(i,j,k), det, & dens_p_p(i,j,k), scon_p_p(i,j,k,1), scon_p_p(i,j,k,2), scon_p_p(i,j,k,3), & - tau_p_p(i,j,k), Bvec_p_p(i,j,k,1),Bvec_p_p(i,j,k,2),Bvec_p_p(i,j,k,3),& + tau_p_p(i,j,k), Bcons_p_p(i,j,k,1),Bcons_p_p(i,j,k,2),Bcons_p_p(i,j,k,3),& rho_p_p(i,j,k), vel_p_p(i,j,k,1), vel_p_p(i,j,k,2), & - vel_p_p(i,j,k,3), eps_p_p(i,j,k), press_p_p(i,j,k), w_lorentz_p_p(i,j,k)) + vel_p_p(i,j,k,3), eps_p_p(i,j,k), press_p_p(i,j,k), & + Bvec_p_p(i,j,k,1),Bvec_p_p(i,j,k,2),Bvec_p_p(i,j,k,3),w_lorentz_p_p(i,j,k)) endif endif |