RIT MHD development:

 Revert previous commit changes to schedule.ccl

 Merge MHD and GRHydro routines: 
    Boundaries, CalcUpdate, PPM, Reconstruct(poly), RegisterGZ,
    RegisterVars

 Solve a few race conditions arising in the "magnetic" routines.
This is still being tested, since it passes tests with intel
compilers but fails for gnu ones.

 Alias function GRHydro_Con2PrimM_pt for C2P2CM test.



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

1 files changed, 76 insertions, 14 deletions

diff --git a/src/GRHydro_Boundaries.F90 b/src/GRHydro_Boundaries.F90
index 1c97c2c..1105eef 100644
--- a/src/GRHydro_Boundaries.F90
+++ b/src/GRHydro_Boundaries.F90
@@ -21,6 +21,9 @@
 #define sx(i,j,k) scon(i,j,k,1)
 #define sy(i,j,k) scon(i,j,k,2)
 #define sz(i,j,k) scon(i,j,k,3)
+#define Bvecx(i,j,k) Bvec(i,j,k,1)
+#define Bvecy(i,j,k) Bvec(i,j,k,2)
+#define Bvecz(i,j,k) Bvec(i,j,k,3)
 
  /*@@
    @routine    GRHydro_InitSymBound
@@ -64,6 +67,9 @@ subroutine GRHydro_InitSymBound(CCTK_ARGUMENTS)
   call SetCartSymVN(ierr, cctkGH, sym, "GRHydro::w_lorentz")
   call SetCartSymVN(ierr, cctkGH, sym, "HydroBase::eps")
   call SetCartSymVN(ierr, cctkGH, sym, "GRHydro::GRHydro_C2P_failed")
+  if(evolve_mhd.ne.0.and.clean_divergence.ne.0) then
+     call SetCartSymVN(ierr, cctkGH, sym, "GRHydro::psidc")
+  endif
 
 !!$ handle multiple tracer variables
   if(evolve_tracer.ne.0) then
@@ -87,6 +93,8 @@ subroutine GRHydro_InitSymBound(CCTK_ARGUMENTS)
   
   call SetCartSymVN(ierr, cctkGH, sym, "HydroBase::vel[0]")
   call SetCartSymVN(ierr, cctkGH, sym, "GRHydro::scon[0]")
+  if(evolve_mhd.ne.0)call SetCartSymVN(ierr, cctkGH, sym, "HydroBase::Bvec[0]")
+
   
   sym(1) = 1
   sym(2) = -1
@@ -94,6 +102,7 @@ subroutine GRHydro_InitSymBound(CCTK_ARGUMENTS)
   
   call SetCartSymVN(ierr, cctkGH, sym, "HydroBase::vel[1]")
   call SetCartSymVN(ierr, cctkGH, sym, "GRHydro::scon[1]")
+  if(evolve_mhd.ne.0)call SetCartSymVN(ierr, cctkGH, sym, "HydroBase::Bvec[1]")
   
   sym(1) = 1
   sym(2) = 1
@@ -101,6 +110,7 @@ subroutine GRHydro_InitSymBound(CCTK_ARGUMENTS)
   
   call SetCartSymVN(ierr, cctkGH, sym, "HydroBase::vel[2]")
   call SetCartSymVN(ierr, cctkGH, sym, "GRHydro::scon[2]")
+  if(evolve_mhd.ne.0)call SetCartSymVN(ierr, cctkGH, sym, "HydroBase::Bvec[2]")
   
 end subroutine GRHydro_InitSymBound
 
@@ -156,6 +166,14 @@ subroutine GRHydro_Boundaries(CCTK_ARGUMENTS)
          "HydroBase::eps", "Flat")
     ierr = Boundary_SelectGroupForBC(cctkGH, faces, GRHydro_stencil, -ione, &
          "HydroBase::vel", "Flat")
+    if(evolve_mhd.ne.0) then
+       ierr = Boundary_SelectGroupForBC(cctkGH, faces, GRHydro_stencil, -ione, &
+            "HydroBase::Bvec", "Flat")
+       if(clean_divergence.ne.0) then
+          ierr = Boundary_SelectGroupForBC(cctkGH, faces, GRHydro_stencil, -ione, &
+               "GRHydro::psidc", "Flat")
+       endif
+    endif
 
     if(evolve_tracer.ne.0) then 
       ierr = Boundary_SelectGroupForBC(cctkGH, faces, GRHydro_stencil, -ione, &
@@ -198,6 +216,14 @@ subroutine GRHydro_Boundaries(CCTK_ARGUMENTS)
          "HydroBase::eps", "None")
     ierr = Boundary_SelectGroupForBC(cctkGH, faces, GRHydro_stencil, -ione, &
          "HydroBase::vel", "None")
+    if(evolve_mhd.ne.0) then
+       ierr = Boundary_SelectGroupForBC(cctkGH, faces, GRHydro_stencil, -ione, &
+            "HydroBase::Bvec", "None")
+       if(clean_divergence.ne.0) then
+          ierr = Boundary_SelectGroupForBC(cctkGH, faces, GRHydro_stencil, -ione, &
+               "GRHydro::psidc", "None")
+       endif
+    endif
 
     if(evolve_tracer.ne.0) then 
       ierr = Boundary_SelectGroupForBC(cctkGH, faces, GRHydro_stencil, -ione, &
@@ -286,17 +312,35 @@ subroutine GRHydro_OutflowBoundaries(CCTK_ARGUMENTS)
               eps(i,j,k)        = eps(i,j,k-1)
               press(i,j,k)      = press(i,j,k-1)
               w_lorentz(i,j,k)  = w_lorentz(i,j,k-1)
-              
+              if(evolve_mhd.ne.0) then
+                 Bvecx(i,j,k)       = Bvecx(i,j,k-1)
+                 Bvecy(i,j,k)       = Bvecy(i,j,k-1)
+                 Bvecz(i,j,k)       = Bvecz(i,j,k-1)
+                 if(clean_divergence.ne.0) then
+                    psidc(i,j,k)    = psidc(i,j,k-1)
+                 endif
+              endif
+  
               psi4pt = 1.0d0
               det = SPATIAL_DETERMINANT(gxx(i,j,k),gxy(i,j,k),gxz(i,j,k),\
                    gyy(i,j,k),gyz(i,j,k),gzz(i,j,k))
-              call prim2con(GRHydro_eos_handle,psi4pt*gxx(i,j,k),&
-                   psi4pt*gxy(i,j,k),psi4pt*gxz(i,j,k),&
-                   psi4pt*gyy(i,j,k),psi4pt*gyz(i,j,k),psi4pt*gzz(i,j,k),&
-                   det, dens(i,j,k),sx(i,j,k),sy(i,j,k),sz(i,j,k),&
-                   tau(i,j,k),rho(i,j,k),velx(i,j,k),vely(i,j,k),velz(i,j,k),&
-                   eps(i,j,k),press(i,j,k),w_lorentz(i,j,k))
               
+              if(evolve_mhd.ne.0) then
+                 call prim2conM(GRHydro_eos_handle,psi4pt*gxx(i,j,k),&
+                      psi4pt*gxy(i,j,k),psi4pt*gxz(i,j,k),&
+                      psi4pt*gyy(i,j,k),psi4pt*gyz(i,j,k),psi4pt*gzz(i,j,k),&
+                      det, dens(i,j,k),sx(i,j,k),sy(i,j,k),sz(i,j,k),&
+                      tau(i,j,k),Bvecx(i,j,k),Bvecy(i,j,k),Bvecz(i,j,k),&
+                      rho(i,j,k),velx(i,j,k),vely(i,j,k),velz(i,j,k),&
+                      eps(i,j,k),press(i,j,k),w_lorentz(i,j,k))
+              else
+                 call prim2con(GRHydro_eos_handle,psi4pt*gxx(i,j,k),&
+                      psi4pt*gxy(i,j,k),psi4pt*gxz(i,j,k),&
+                      psi4pt*gyy(i,j,k),psi4pt*gyz(i,j,k),psi4pt*gzz(i,j,k),&
+                      det, dens(i,j,k),sx(i,j,k),sy(i,j,k),sz(i,j,k),&
+                      tau(i,j,k),rho(i,j,k),velx(i,j,k),vely(i,j,k),velz(i,j,k),&
+                      eps(i,j,k),press(i,j,k),w_lorentz(i,j,k))
+              endif
             end if
           
           end do
@@ -327,16 +371,34 @@ subroutine GRHydro_OutflowBoundaries(CCTK_ARGUMENTS)
               press(i,j,k)      = press(i,j,k+1)
               w_lorentz(i,j,k)  = w_lorentz(i,j,k+1)
               
+              if(evolve_mhd.ne.0) then
+                 Bvecx(i,j,k)       = Bvecx(i,j,k+1)
+                 Bvecy(i,j,k)       = Bvecy(i,j,k+1)
+                 Bvecz(i,j,k)       = Bvecz(i,j,k+1)
+                 if(clean_divergence.ne.0) then
+                    psidc(i,j,k)    = psidc(i,j,k+1)
+                 endif
+              endif
+
               psi4pt = 1.0d0
               det = SPATIAL_DETERMINANT(gxx(i,j,k),gxy(i,j,k),gxz(i,j,k),\
                    gyy(i,j,k),gyz(i,j,k),gzz(i,j,k))
-              call prim2con(GRHydro_eos_handle,psi4pt*gxx(i,j,k),&
-                   psi4pt*gxy(i,j,k),psi4pt*gxz(i,j,k),&
-                   psi4pt*gyy(i,j,k),psi4pt*gyz(i,j,k),psi4pt*gzz(i,j,k),&
-                   det, dens(i,j,k),sx(i,j,k),sy(i,j,k),sz(i,j,k),&
-                   tau(i,j,k),rho(i,j,k),velx(i,j,k),vely(i,j,k),velz(i,j,k),&
-                   eps(i,j,k),press(i,j,k),w_lorentz(i,j,k))
-              
+              if(evolve_mhd.ne.0) then
+                 call prim2conM(GRHydro_eos_handle,psi4pt*gxx(i,j,k),&
+                      psi4pt*gxy(i,j,k),psi4pt*gxz(i,j,k),&
+                      psi4pt*gyy(i,j,k),psi4pt*gyz(i,j,k),psi4pt*gzz(i,j,k),&
+                      det, dens(i,j,k),sx(i,j,k),sy(i,j,k),sz(i,j,k),&
+                      tau(i,j,k),Bvecx(i,j,k),Bvecy(i,j,k),Bvecz(i,j,k),&
+                      rho(i,j,k),velx(i,j,k),vely(i,j,k),velz(i,j,k),&
+                      eps(i,j,k),press(i,j,k),w_lorentz(i,j,k))
+              else
+                 call prim2con(GRHydro_eos_handle,psi4pt*gxx(i,j,k),&
+                      psi4pt*gxy(i,j,k),psi4pt*gxz(i,j,k),&
+                      psi4pt*gyy(i,j,k),psi4pt*gyz(i,j,k),psi4pt*gzz(i,j,k),&
+                      det, dens(i,j,k),sx(i,j,k),sy(i,j,k),sz(i,j,k),&
+                      tau(i,j,k),rho(i,j,k),velx(i,j,k),vely(i,j,k),velz(i,j,k),&
+                      eps(i,j,k),press(i,j,k),w_lorentz(i,j,k))
+              endif
             end if
           
           end do