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.



git-svn-id: http://svn.einsteintoolkit.org/cactus/EinsteinEvolve/GRHydro/trunk@196 c83d129a-5a75-4d5a-9c4d-ed3a5855bf45

9 files changed, 102 insertions, 87 deletions

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 \