User: bmundim Date: 2010/12/21 11:35 PM

Added: /trunk/src/ GRHydro_DivergenceClean.F90

Removed: /trunk/src/ GRHydro_ParamCheckM.F90

Modified: /trunk/ interface.ccl, param.ccl, schedule.ccl /trunk/src/ GRHydro_CalcUpdateM.F90, GRHydro_HLLEM.F90, GRHydro_HLLE_EOSGM.F90, GRHydro_ParamCheck.F90, GRHydro_RiemannSolveM.F90, GRHydro_SourceM.F90, make.code.defn

Log: RIT MHD development: an update.

Several bugfixes, including properly isolating divergence cleaning code and initialization of the "psidc" gridfunction in GRHydro_DivergenceClean.F90, as well as properly initializing rhoenth in GRHydro_HLLEM (thanks to Roland for catching that). Minor changes to scheduler to properly call MHD and non-MHD routines in parallel for Riemann solving. Remove GRHydro_ParamCheckM.F90 and add GRHydro_DivergenceClean.F90.

File Changes:

Directory: /trunk/src/ ======================

File [modified]: GRHydro_CalcUpdateM.F90 Delta lines: +10 -7 =================================================================== --- trunk/src/GRHydro_CalcUpdateM.F90 2010-12-21 10:58:39 UTC (rev 195) +++ trunk/src/GRHydro_CalcUpdateM.F90 2010-12-22 05:35:43 UTC (rev 196) @@ -89,10 +89,11 @@ 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 @@ 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

File [added]: GRHydro_DivergenceClean.F90 Delta lines: +40 -0 =================================================================== --- trunk/src/GRHydro_DivergenceClean.F90 (rev 0) +++ trunk/src/GRHydro_DivergenceClean.F90 2010-12-22 05:35:43 UTC (rev 196) @@ -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 + + +

File [modified]: GRHydro_HLLEM.F90 Delta lines: +6 -1 =================================================================== --- trunk/src/GRHydro_HLLEM.F90 2010-12-21 10:58:39 UTC (rev 195) +++ trunk/src/GRHydro_HLLEM.F90 2010-12-22 05:35:43 UTC (rev 196) @@ -187,6 +187,9 @@ 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 @@ 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

File [modified]: GRHydro_HLLE_EOSGM.F90 Delta lines: +4 -1 =================================================================== --- trunk/src/GRHydro_HLLE_EOSGM.F90 2010-12-21 10:58:39 UTC (rev 195) +++ trunk/src/GRHydro_HLLE_EOSGM.F90 2010-12-22 05:35:43 UTC (rev 196) @@ -438,7 +438,10 @@ 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

File [modified]: GRHydro_ParamCheck.F90 Delta lines: +9 -2 =================================================================== --- trunk/src/GRHydro_ParamCheck.F90 2010-12-21 10:58:39 UTC (rev 195) +++ trunk/src/GRHydro_ParamCheck.F90 2010-12-22 05:35:43 UTC (rev 196) @@ -87,9 +87,9 @@ 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 @@ 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

File [removed]: GRHydro_ParamCheckM.F90 Delta lines: +0 -50 =================================================================== --- trunk/src/GRHydro_ParamCheckM.F90 2010-12-21 10:58:39 UTC (rev 195) +++ trunk/src/GRHydro_ParamCheckM.F90 2010-12-22 05:35:43 UTC (rev 196) @@ -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 -

File [modified]: GRHydro_RiemannSolveM.F90 Delta lines: +0 -1 =================================================================== --- trunk/src/GRHydro_RiemannSolveM.F90 2010-12-21 10:58:39 UTC (rev 195) +++ trunk/src/GRHydro_RiemannSolveM.F90 2010-12-22 05:35:43 UTC (rev 196) @@ -362,7 +362,6 @@ 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)

File [modified]: GRHydro_SourceM.F90 Delta lines: +32 -24 =================================================================== --- trunk/src/GRHydro_SourceM.F90 2010-12-21 10:58:39 UTC (rev 195) +++ trunk/src/GRHydro_SourceM.F90 2010-12-22 05:35:43 UTC (rev 196) @@ -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 @@ 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 @@ endif

if (clean_divergence .ne. 0) then - Bvecrhs=0.d0 psidcrhs=0.d0 endif

@@ -294,19 +294,23 @@ 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 @@ 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

File [modified]: make.code.defn Delta lines: +1 -1 =================================================================== --- trunk/src/make.code.defn 2010-12-21 10:58:39 UTC (rev 195) +++ trunk/src/make.code.defn 2010-12-22 05:35:43 UTC (rev 196) @@ -7,6 +7,7 @@ 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 @@ GRHydro_HLLE_EOSGM.F90 \ GRHydro_HLLEM.F90 \ GRHydro_PPMM.F90 \ - GRHydro_ParamCheckM.F90 \ GRHydro_Prim2ConM.F90 \ GRHydro_RiemannSolveM.F90 \ GRHydro_ReconstructM.F90 \

Directory: /trunk/ ==================

File [modified]: interface.ccl Delta lines: +1 -1 =================================================================== --- trunk/interface.ccl 2010-12-21 10:58:39 UTC (rev 195) +++ trunk/interface.ccl 2010-12-22 05:35:43 UTC (rev 196) @@ -392,7 +392,7 @@

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"

File [modified]: param.ccl Delta lines: +1 -1 =================================================================== --- trunk/param.ccl 2010-12-21 10:58:39 UTC (rev 195) +++ trunk/param.ccl 2010-12-22 05:35:43 UTC (rev 196) @@ -76,7 +76,7 @@

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

File [modified]: schedule.ccl Delta lines: +18 -18 =================================================================== --- trunk/schedule.ccl 2010-12-21 10:58:39 UTC (rev 195) +++ trunk/schedule.ccl 2010-12-22 05:35:43 UTC (rev 196) @@ -31,11 +31,11 @@ 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 @@ 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 @@ { } "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 @@ 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 @@ 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