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

1 files changed, 51 insertions, 15 deletions

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.