From 4adfcb611a12689096a604e2a510b7ff12c145f1 Mon Sep 17 00:00:00 2001 From: bmundim Date: Wed, 13 Apr 2011 03:29:13 +0000 Subject: RIT MHD development: 1) Fix a bug in the divergence cleaning implementation: a psidc term in the induction equation was implemented as a source term where it was supposed to be coded as part of the flux calculation. 2) Introduce divB as a diagnostic grid function. 3) Remove an old file: GRHydro_CalcUpdateM.F90 git-svn-id: http://svn.einsteintoolkit.org/cactus/EinsteinEvolve/GRHydro/trunk@228 c83d129a-5a75-4d5a-9c4d-ed3a5855bf45 --- src/GRHydro_HLLEM.F90 | 66 +++++++++++++++++++++++++++++++++++++++------------ 1 file changed, 51 insertions(+), 15 deletions(-) (limited to 'src/GRHydro_HLLEM.F90') diff --git a/src/GRHydro_HLLEM.F90 b/src/GRHydro_HLLEM.F90 index c5a15bb..56e19a1 100644 --- a/src/GRHydro_HLLEM.F90 +++ b/src/GRHydro_HLLEM.F90 @@ -172,6 +172,11 @@ subroutine GRHydro_HLLEM(CCTK_ARGUMENTS) avg_det = SPATIAL_DETERMINANT(gxxh,gxyh,gxzh,gyyh,gyzh,gzzh) + call UpperMetric(uxxh, uxyh, uxzh, uyyh, uyzh, uzzh, & + avg_det,gxxh, gxyh, gxzh, & + gyyh, gyzh, gzzh) + + vxtp = prim_p(2)-avg_betax/avg_alp vytp = prim_p(3)-avg_betay/avg_alp vztp = prim_p(4)-avg_betaz/avg_alp @@ -218,7 +223,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 - psidcf = cons_p(6) + 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) endif else if (flux_direction == 2) then @@ -228,7 +236,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 - psidcf = cons_p(7) + 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) endif else if (flux_direction == 3) then @@ -238,7 +249,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 - psidcf = cons_p(8) + 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) endif else @@ -247,10 +261,6 @@ subroutine GRHydro_HLLEM(CCTK_ARGUMENTS) else !!! The end of this branch is right at the bottom of the routine - call UpperMetric(uxxh, uxyh, uxzh, uyyh, uyzh, uzzh, & - avg_det,gxxh, gxyh, gxzh, & - gyyh, gyzh, gzzh) - if (flux_direction == 1) then usendh = uxxh else if (flux_direction == 2) then @@ -302,8 +312,14 @@ subroutine GRHydro_HLLEM(CCTK_ARGUMENTS) avg_det,avg_alp,avg_beta) if(clean_divergence.ne.0) then - psidcfp = cons_p(6) - psidcfm = cons_m(6) + 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) endif else if (flux_direction == 2) then @@ -329,8 +345,14 @@ 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 - psidcfp = cons_p(7) - psidcfm = cons_m(7) + 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) endif else if (flux_direction == 3) then @@ -356,8 +378,14 @@ 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 - psidcfp = cons_p(8) - psidcfm = cons_m(8) + 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) endif else @@ -389,8 +417,16 @@ subroutine GRHydro_HLLEM(CCTK_ARGUMENTS) if(clean_divergence.ne.0) then psidcdiff = psidcm - psidcp - psidcf = (charmax * psidcfp - charmin * psidcfm + & - charmax * charmin * psidcdiff) / charpm + +!!$ psidcf = (charmax * psidcfp - charmin * psidcfm + & +!!$ charmax * charmin * psidcdiff) / charpm + +!!$ Wavespeeds for psidc are +/-c, not Fast Magnetosonic? + + psidcf = (1.d0 * psidcfp - (-1.d0) * psidcfm + & + 1.d0 * (-1.d0) * psidcdiff) / 2.d0 + + endif end if !!! The end of the SpaceMask check for a trivial RP. -- cgit v1.2.3