diff options
| -rw-r--r-- | interface.ccl | 2 | ||||
| -rw-r--r-- | param.ccl | 2 | ||||
| -rw-r--r-- | schedule.ccl | 36 | ||||
| -rw-r--r-- | src/GRHydro_CalcUpdateM.F90 | 17 | ||||
| -rw-r--r-- | src/GRHydro_DivergenceClean.F90 | 40 | ||||
| -rw-r--r-- | src/GRHydro_HLLEM.F90 | 7 | ||||
| -rw-r--r-- | src/GRHydro_HLLE_EOSGM.F90 | 5 | ||||
| -rw-r--r-- | src/GRHydro_ParamCheck.F90 | 11 | ||||
| -rw-r--r-- | src/GRHydro_ParamCheckM.F90 | 50 | ||||
| -rw-r--r-- | src/GRHydro_RiemannSolveM.F90 | 1 | ||||
| -rw-r--r-- | src/GRHydro_SourceM.F90 | 56 | ||||
| -rw-r--r-- | src/make.code.defn | 2 |
12 files changed, 122 insertions, 107 deletions
diff --git a/interface.ccl b/interface.ccl index fc9e930..ce31620 100644 --- a/interface.ccl +++ b/interface.ccl @@ -392,7 +392,7 @@ int evolve_Y_e type = SCALAR tags='checkpoint="no"' "Are we evolving Y_e? Set in int evolve_temper type = SCALAR tags='checkpoint="no"' "Are we evolving temperature? Set in Paramcheck" -int MHD type = SCALAR tags='checkpoint="no"' "Are we doing MHD? Set in ParamCheck" +int evolve_MHD type = SCALAR tags='checkpoint="no"' "Are we doing MHD? Set in ParamCheck" int GRHydro_reflevel type = SCALAR tags='checkpoint="no"' "Refinement level GRHydro is working on right now" @@ -76,7 +76,7 @@ CCTK_INT GRHydro_hydro_excision "Turns excision automatically on in HydroBase" A CCTK_INT GRHydro_MaxNumEvolvedVars "The maximum number of evolved variables used by GRHydro" ACCUMULATOR-BASE=MethodofLines::MoL_Num_Evolved_Vars { - 5:9 :: "dens scon[3] tau Bvec[3] ye" + 5:10 :: "dens scon[3] tau Bvec[3] psidc ye" } 5 CCTK_INT GRHydro_MaxNumConstrainedVars "The maximum number of constrained variables used by GRHydro" ACCUMULATOR-BASE=MethodofLines::MoL_Num_Constrained_Vars diff --git a/schedule.ccl b/schedule.ccl index 42295fe..1f79f13 100644 --- a/schedule.ccl +++ b/schedule.ccl @@ -31,10 +31,10 @@ if (timelevels == 3) if(CCTK_Equals(Bvec_evolution_method,"GRHydro")) { STORAGE: HydroBase::Bvec[3] - } - if (clean_divergence) - { - STORAGE:psidc[3] + if (clean_divergence) + { + STORAGE:psidc[3] + } } } else @@ -57,13 +57,13 @@ else if(CCTK_Equals(Bvec_evolution_method,"GRHydro")) { STORAGE: HydroBase::Bvec[2] - } - if (clean_divergence) - { - STORAGE:psidc[2] + if (clean_divergence) + { + STORAGE:psidc[2] + } } } -STORAGE:MHD +STORAGE:evolve_MHD STORAGE:evolve_Y_e STORAGE:evolve_temper STORAGE:GRHydro_reflevel @@ -99,6 +99,14 @@ schedule group GRHydro_Initial IN HydroBase_Initial BEFORE SetTmunu { } "GRHydro initial data group" +if(CCTK_Equals(Bvec_evolution_method,"GRHydro") && clean_divergence) +{ + schedule GRHydro_InitDivergenceClean IN HydroBase_Initial + { + LANG: Fortran + } "Set psi for divergence cleaning initially to zero" +} + ################################################# ### Storage for the extra timelevels for the ### ### use of MoL with Einstein ### @@ -188,14 +196,6 @@ schedule GRHydro_ParamCheck AT PARAMCHECK LANG: Fortran } "Check parameters" -if(CCTK_Equals(Bvec_evolution_method,"GRHydro")) -{ - schedule GRHydro_ParamCheckM AT PARAMCHECK AFTER GRHydro_ParamCheck - { - LANG: Fortran - } "Check parameters - MHD version" -} - ###################################### ### Standard symmetry registration ### ###################################### @@ -733,7 +733,7 @@ if (CCTK_Equals(method_type, "RSA FV")) else if (CCTK_Equals(GRHydro_eos_type,"General")) { if(CCTK_Equals(Bvec_evolution_method,"GRHydro")) { - schedule RiemannSolve IN FluxTerms AFTER Reconstruct AS Riemann + schedule RiemannSolveM IN FluxTerms AFTER Reconstruct AS Riemann { LANG: Fortran STORAGE: EOS_temps diff --git a/src/GRHydro_CalcUpdateM.F90 b/src/GRHydro_CalcUpdateM.F90 index e5296ce..9bfdb75 100644 --- a/src/GRHydro_CalcUpdateM.F90 +++ b/src/GRHydro_CalcUpdateM.F90 @@ -89,10 +89,11 @@ subroutine UpdateCalculationM(CCTK_ARGUMENTS) 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 - psidcrhs(i,j,k) = psidcrhs(i,j,k) + & - (alp_l * psidcflux(i-xoffset,j-yoffset,k-zoffset) - & - alp_r * psidcflux(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 if (evolve_tracer .ne. 0) then do itracer=1,number_of_tracers @@ -236,9 +237,11 @@ subroutine UpdateCalculationM(CCTK_ARGUMENTS) Bvecrhs(i,j,k,3) = Bvecrhs(i,j,k,3) + & (Bveczflux(i-xoffset,j-yoffset,k-zoffset) - & Bveczflux(i,j,k)) * idx - psidcrhs(i,j,k) = psidcrhs(i,j,k) + & - (psidcflux(i-xoffset,j-yoffset,k-zoffset) - & - psidcflux(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 enddo enddo diff --git a/src/GRHydro_DivergenceClean.F90 b/src/GRHydro_DivergenceClean.F90 new file mode 100644 index 0000000..23c2918 --- /dev/null +++ b/src/GRHydro_DivergenceClean.F90 @@ -0,0 +1,40 @@ + /*@@ + @file GRHydro_DivergenceClean.F90 + @date Nov 24, 2010 + @author + @desc + Routines controlling divergence cleaning + @enddesc + @@*/ + +#include "cctk.h" +#include "cctk_Parameters.h" +#include "cctk_Arguments.h" + + /*@@ + @routine GRHydro_InitDivergenceClean + @date Nov 24, 2010 + @author Joshua Faber + @desc + Set Psi=0 initially for use with divergence cleaning + @enddesc + @calls + @calledby + @history + + @endhistory + +@@*/ + +subroutine GRHydro_InitDivergenceClean(CCTK_ARGUMENTS) + + implicit none + + DECLARE_CCTK_ARGUMENTS + + psidc=0.0 + +end subroutine GRHydro_InitDivergenceClean + + + diff --git a/src/GRHydro_HLLEM.F90 b/src/GRHydro_HLLEM.F90 index aec2899..c5a15bb 100644 --- a/src/GRHydro_HLLEM.F90 +++ b/src/GRHydro_HLLEM.F90 @@ -187,6 +187,9 @@ subroutine GRHydro_HLLEM(CCTK_ARGUMENTS) prim_m(2),prim_m(3),prim_m(4),cons_m(6),cons_m(7),cons_m(8), & velxlowm,velylowm,velzlowm,Bvecxlowm,Bvecylowm,Bveczlowm, & bdotvm,b2m,v2m,wm,bxlowm,bylowm,bzlowm) + + rhoenth_p = prim_p(1)*(1.d0+prim_p(5))+prim_p(6) + rhoenth_m = prim_m(1)*(1.d0+prim_m(5))+prim_m(6) ab0p = wp*bdotvp ab0m = wm*bdotvm @@ -400,7 +403,9 @@ subroutine GRHydro_HLLEM(CCTK_ARGUMENTS) Bvecxflux(i, j, k) = f1(6) Bvecyflux(i, j, k) = f1(7) Bveczflux(i, j, k) = f1(8) - psidcflux(i,j,k) = psidcf + if(clean_divergence.ne.0) then + psidcflux(i,j,k) = psidcf + endif if(evolve_Y_e.ne.0) then if (densflux(i, j, k) > 0.d0) then diff --git a/src/GRHydro_HLLE_EOSGM.F90 b/src/GRHydro_HLLE_EOSGM.F90 index 8fe2ef0..c89d047 100644 --- a/src/GRHydro_HLLE_EOSGM.F90 +++ b/src/GRHydro_HLLE_EOSGM.F90 @@ -438,7 +438,10 @@ subroutine GRHydro_HLLEGeneralM(CCTK_ARGUMENTS) Bvecxflux(i, j, k) = f1(6) Bvecyflux(i, j, k) = f1(7) Bveczflux(i, j, k) = f1(8) - psidcflux(i,j,k) = psidcf + if(clean_divergence.ne.0) then + psidcflux(i,j,k) = psidcf + endif + if(evolve_Y_e.ne.0) then if (densflux(i, j, k) > 0.d0) then diff --git a/src/GRHydro_ParamCheck.F90 b/src/GRHydro_ParamCheck.F90 index 93dcaa3..f6177f7 100644 --- a/src/GRHydro_ParamCheck.F90 +++ b/src/GRHydro_ParamCheck.F90 @@ -87,9 +87,9 @@ subroutine GRHydro_ParamCheck(CCTK_ARGUMENTS) end if if (CCTK_EQUALS(Bvec_evolution_method,"GRHydro")) then - MHD = 1 + evolve_MHD = 1 else - MHD = 0 + evolve_MHD = 0 endif if (CCTK_EQUALS(Y_e_evolution_method,"GRHydro")) then @@ -104,6 +104,13 @@ subroutine GRHydro_ParamCheck(CCTK_ARGUMENTS) evolve_temper = 0 endif + if (CCTK_EQUALS(riemann_solver,"Roe").and.CCTK_EQUALS(Bvec_evolution_method,"GRHydro")) then + call CCTK_PARAMWARN("Roe solver is not implemented yet for MHD") + end if + + if (CCTK_EQUALS(riemann_solver,"Marquina").and.CCTK_EQUALS(Bvec_evolution_method,"GRHydro")) then + call CCTK_PARAMWARN("Marquina solver is not implemented yet for MHD") + end if end subroutine GRHydro_ParamCheck diff --git a/src/GRHydro_ParamCheckM.F90 b/src/GRHydro_ParamCheckM.F90 deleted file mode 100644 index 84615eb..0000000 --- a/src/GRHydro_ParamCheckM.F90 +++ /dev/null @@ -1,50 +0,0 @@ - /*@@ - @file GRHydro_ParamCheckM.F90 - @date Sep 23, 2010 - @author - @desc - Parameter checking routine. - @enddesc - @@*/ - -#include "cctk.h" -#include "cctk_Parameters.h" -#include "cctk_Arguments.h" -#include "cctk_Functions.h" - - /*@@ - @routine GRHydro_ParamCheckM - @date Sep 23, 2010 - @author Joshua Faber, Scott Noble, Bruno Mundim - @desc - Checks the parameters. - @enddesc - @calls - @calledby - @history - - @endhistory - -@@*/ - -subroutine GRHydro_ParamCheckM(CCTK_ARGUMENTS) - - implicit none - - DECLARE_CCTK_ARGUMENTS - DECLARE_CCTK_PARAMETERS - DECLARE_CCTK_FUNCTIONS - - - - if (CCTK_EQUALS(riemann_solver,"Roe").and.CCTK_EQUALS(Bvec_evolution_method,"GRHydro")) then - call CCTK_PARAMWARN("Roe solver is not implemented yet for MHD") - end if - - if (CCTK_EQUALS(riemann_solver,"Marquina").and.CCTK_EQUALS(Bvec_evolution_method,"GRHydro")) then - call CCTK_PARAMWARN("Marquina solver is not implemented yet for MHD") - end if - - -end subroutine GRHydro_ParamCheckM - diff --git a/src/GRHydro_RiemannSolveM.F90 b/src/GRHydro_RiemannSolveM.F90 index d0eb991..8944b26 100644 --- a/src/GRHydro_RiemannSolveM.F90 +++ b/src/GRHydro_RiemannSolveM.F90 @@ -362,7 +362,6 @@ subroutine RiemannSolveGeneralM(CCTK_ARGUMENTS) Bvecxflux(i,j,k)= Bvecxflux(i,j,k)+ 0.5d0 * tmp_flux(6) Bvecyflux(i,j,k)= Bvecyflux(i,j,k)+ 0.5d0 * tmp_flux(7) Bveczflux(i,j,k)= Bveczflux(i,j,k)+ 0.5d0 * tmp_flux(8) - psidcflux(i,j,k)= psidcflux(i,j,k)+ 0.5d0 * psidcf if(clean_divergence.ne.0) then psidcf = cons_m(6) diff --git a/src/GRHydro_SourceM.F90 b/src/GRHydro_SourceM.F90 index 15c7aa0..2a5de76 100644 --- a/src/GRHydro_SourceM.F90 +++ b/src/GRHydro_SourceM.F90 @@ -9,18 +9,18 @@ ! Second order f.d. -#define DIFF_X_2(q) 0.5d0 * (q(i+1,j,k) - q(i-1,j,k)) * idx -#define DIFF_Y_2(q) 0.5d0 * (q(i,j+1,k) - q(i,j-1,k)) * idy -#define DIFF_Z_2(q) 0.5d0 * (q(i,j,k+1) - q(i,j,k-1)) * idz +#define DIFF_X_2(q) (0.5d0 * (q(i+1,j,k) - q(i-1,j,k)) * idx) +#define DIFF_Y_2(q) (0.5d0 * (q(i,j+1,k) - q(i,j-1,k)) * idy) +#define DIFF_Z_2(q) (0.5d0 * (q(i,j,k+1) - q(i,j,k-1)) * idz) ! Fourth order f.d. -#define DIFF_X_4(q) (-q(i+2,j,k) + 8.d0 * q(i+1,j,k) - 8.d0 * q(i-1,j,k) + \ - q(i-2,j,k)) / 12.d0 * idx -#define DIFF_Y_4(q) (-q(i,j+2,k) + 8.d0 * q(i,j+1,k) - 8.d0 * q(i,j-1,k) + \ - q(i,j-2,k)) / 12.d0 * idy -#define DIFF_Z_4(q) (-q(i,j,k+2) + 8.d0 * q(i,j,k+1) - 8.d0 * q(i,j,k-1) + \ - q(i,j,k-2)) / 12.d0 * idz +#define DIFF_X_4(q) ((-q(i+2,j,k) + 8.d0 * q(i+1,j,k) - 8.d0 * q(i-1,j,k) + \ + q(i-2,j,k)) / 12.d0 * idx) +#define DIFF_Y_4(q) ((-q(i,j+2,k) + 8.d0 * q(i,j+1,k) - 8.d0 * q(i,j-1,k) + \ + q(i,j-2,k)) / 12.d0 * idy) +#define DIFF_Z_4(q) ((-q(i,j,k+2) + 8.d0 * q(i,j,k+1) - 8.d0 * q(i,j,k-1) + \ + q(i,j,k-2)) / 12.d0 * idz) #include "cctk.h" #include "cctk_Parameters.h" @@ -103,6 +103,7 @@ subroutine SourceTermsM(CCTK_ARGUMENTS) densrhs = 0.d0 srhs = 0.d0 taurhs = 0.d0 + Bvecrhs=0.d0 if (evolve_tracer .ne. 0) then cons_tracerrhs = 0.d0 @@ -113,7 +114,6 @@ subroutine SourceTermsM(CCTK_ARGUMENTS) endif if (clean_divergence .ne. 0) then - Bvecrhs=0.d0 psidcrhs=0.d0 endif @@ -294,19 +294,23 @@ subroutine SourceTermsM(CCTK_ARGUMENTS) dx_alp = DIFF_X_2(alp) dy_alp = DIFF_Y_2(alp) dz_alp = DIFF_Z_2(alp) - dx_psidc = DIFF_X_2(psidc) - dy_psidc = DIFF_Y_2(psidc) - dz_psidc = DIFF_Z_2(psidc) + + if(clean_divergence.ne.0) then + dx_psidc = DIFF_X_2(psidc) + dy_psidc = DIFF_Y_2(psidc) + dz_psidc = DIFF_Z_2(psidc) + endif else dx_alp = DIFF_X_4(alp) dy_alp = DIFF_Y_4(alp) dz_alp = DIFF_Z_4(alp) - dx_psidc = DIFF_X_4(psidc) - dy_psidc = DIFF_Y_4(psidc) - dz_psidc = DIFF_Z_4(psidc) - + if(clean_divergence.ne.0) then + dx_psidc = DIFF_X_4(psidc) + dy_psidc = DIFF_Y_4(psidc) + dz_psidc = DIFF_Z_4(psidc) + endif end if if (local_spatial_order .eq. 2) then @@ -432,19 +436,23 @@ subroutine SourceTermsM(CCTK_ARGUMENTS) rhohstarW2*invalp*(vlowx*dz_betax + vlowy*dz_betay + vlowz*dz_betaz) -& bt*(bxlow*dz_betax + bylow*dz_betay + bzlow*dz_betaz) - bvcx_source = uxx*dx_psidc+uxy*dy_psidc+uxz*dz_psidc - bvcy_source = uxy*dx_psidc+uyy*dy_psidc+uyz*dz_psidc - bvcz_source = uxz*dx_psidc+uyz*dy_psidc+uzz*dz_psidc + if(clean_divergence.ne.0) then + bvcx_source = uxx*dx_psidc+uxy*dy_psidc+uxz*dz_psidc + bvcy_source = uxy*dx_psidc+uyy*dy_psidc+uyz*dz_psidc + bvcz_source = uxz*dx_psidc+uyz*dy_psidc+uzz*dz_psidc + endif densrhs(i,j,k) = 0.d0 srhs(i,j,k,1) = alp(i,j,k)*sqrtdet*sx_source 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 - Bvecrhs(i,j,k,1) = alp(i,j,k)*sqrtdet*bvcx_source - Bvecrhs(i,j,k,2) = alp(i,j,k)*sqrtdet*bvcy_source - Bvecrhs(i,j,k,3) = alp(i,j,k)*sqrtdet*bvcz_source - psidcrhs(i,j,k) = -1.d0*ch_dc*ch_dc/cp_dc/cp_dc*psidc(i,j,k) + if(clean_divergence.ne.0) then + Bvecrhs(i,j,k,1) = alp(i,j,k)*sqrtdet*bvcx_source + Bvecrhs(i,j,k,2) = alp(i,j,k)*sqrtdet*bvcy_source + Bvecrhs(i,j,k,3) = alp(i,j,k)*sqrtdet*bvcz_source + psidcrhs(i,j,k) = -1.d0*ch_dc*ch_dc/cp_dc/cp_dc*psidc(i,j,k) + endif enddo enddo diff --git a/src/make.code.defn b/src/make.code.defn index 2a00a1b..b28c974 100644 --- a/src/make.code.defn +++ b/src/make.code.defn @@ -7,6 +7,7 @@ SRCS = Utils.F90 \ GRHydro_Boundaries.F90 \ GRHydro_CalcUpdate.F90 \ GRHydro_Con2Prim.F90 \ + GRHydro_DivergenceClean.F90 \ GRHydro_Eigenproblem.F90 \ GRHydro_Eigenproblem_Marquina.F90 \ GRHydro_ENOReconstruct.F90 \ @@ -55,7 +56,6 @@ SRCS = Utils.F90 \ GRHydro_HLLE_EOSGM.F90 \ GRHydro_HLLEM.F90 \ GRHydro_PPMM.F90 \ - GRHydro_ParamCheckM.F90 \ GRHydro_Prim2ConM.F90 \ GRHydro_RiemannSolveM.F90 \ GRHydro_ReconstructM.F90 \ |