Add a routine to calculate mag pressure/comoving field

    
    * When Tmunu is not calculated, ie when no storage is selected,
    the comoving magnetic field/magnetic pressure wasn't calculated
    at all. This commit allows that now.



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

3 files changed, 77 insertions, 0 deletions

diff --git a/schedule.ccl b/schedule.ccl
index fde07cb..392086a 100644
--- a/schedule.ccl
+++ b/schedule.ccl
@@ -1139,6 +1139,16 @@ if (CCTK_Equals(Bvec_evolution_method,"GRHydro"))
   {
     LANG: Fortran
   } "Compute the energy-momentum tensor - MHD version"
+# bcom is already calculated in GRHydro_TmunuM. 
+# When Tmunu is not calculated, then the following 
+# routine is used to evaluate bcom, if requested.
+  if (calculate_bcom){
+    schedule GRHydro_CalcBcom IN HydroBase_PostStep
+    {
+      LANG: Fortran
+    } "Compute comoving magnetic field, pressure, etc..."
+  }
+
 } else {
   schedule GRHydro_Tmunu IN AddToTmunu 
   {
diff --git a/src/GRHydro_CalcBcom.F90 b/src/GRHydro_CalcBcom.F90
new file mode 100644
index 0000000..b1acd27
--- /dev/null
+++ b/src/GRHydro_CalcBcom.F90
@@ -0,0 +1,66 @@
+ /*@@
+   @file      GRHydro_CalcBcom.F90
+   @date      Jan 01, 2012
+   @author    Bruno Mundim
+   @histpry
+    Based on GRHydro_TmunuM.F90 file.
+   @desc 
+     The calculation of the magnetic pressure and the comoving magnetic
+     field.
+   @enddesc 
+ @@*/
+#include "cctk.h"
+#include "cctk_Arguments.h"
+
+#define velx(i,j,k) vel(i,j,k,1)
+#define vely(i,j,k) vel(i,j,k,2)
+#define velz(i,j,k) vel(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)
+#define bcomx(i,j,k) bcom(i,j,k,1)
+#define bcomy(i,j,k) bcom(i,j,k,2)
+#define bcomz(i,j,k) bcom(i,j,k,3)
+
+ subroutine GRHydro_CalcBcom(CCTK_ARGUMENTS)
+
+   implicit none
+
+   DECLARE_CCTK_ARGUMENTS
+
+   CCTK_REAL :: velxlow, velylow, velzlow
+   CCTK_REAL :: Bvecxlow,Bvecylow,Bveczlow
+   CCTK_REAL :: bdotv,b2,bxlow,bylow,bzlow,btlow,dum
+   CCTK_INT :: i,j,k
+
+   if(stress_energy_state.ne.0) then
+     return
+   end if
+
+   !$OMP PARALLEL DO PRIVATE(i,j,k,velxlow, velylow, velzlow,&
+   !$OMP Bvecxlow,Bvecylow,Bveczlow, bdotv,dum,b2,bxlow,bylow,bzlow)
+
+   do k = 1, cctk_lsh(3)
+     do j = 1, cctk_lsh(2)
+       do i = 1, cctk_lsh(1)
+
+          call calc_vlow_blow(gxx(i,j,k),gxy(i,j,k),gxz(i,j,k),&
+               gyy(i,j,k),gyz(i,j,k),gzz(i,j,k), &
+               velx(i,j,k),vely(i,j,k),velz(i,j,k),Bvecx(i,j,k),Bvecy(i,j,k),Bvecz(i,j,k), &
+               velxlow,velylow,velzlow,Bvecxlow,Bvecylow,Bveczlow, &
+               bdotv,b2,dum,dum,bxlow,bylow,bzlow)
+
+           bcom_sq(i,j,k) = b2 
+           bcom0(i,j,k) = w_lorentz(i,j,k)*bdotv/alp(i,j,k)
+           bcomx(i,j,k) = Bvecx(i,j,k)/w_lorentz(i,j,k) + bcom0(i,j,k)*(alp(i,j,k)*velx(i,j,k)-betax(i,j,k))
+           bcomy(i,j,k) = Bvecy(i,j,k)/w_lorentz(i,j,k) + bcom0(i,j,k)*(alp(i,j,k)*vely(i,j,k)-betay(i,j,k))
+           bcomz(i,j,k) = Bvecz(i,j,k)/w_lorentz(i,j,k) + bcom0(i,j,k)*(alp(i,j,k)*velz(i,j,k)-betaz(i,j,k))
+
+       end do
+     end do
+   end do
+   !$OMP END PARALLEL DO
+
+   return
+ 
+ end subroutine GRHydro_CalcBcom
diff --git a/src/make.code.defn b/src/make.code.defn
index 414313b..50cd1bb 100644
--- a/src/make.code.defn
+++ b/src/make.code.defn
@@ -5,6 +5,7 @@
 
 SRCS = 	Utils.F90 \
 	GRHydro_Boundaries.F90 \
+	GRHydro_CalcBcom.F90 \
 	GRHydro_CalcUpdate.F90 \
 	GRHydro_Con2Prim.F90 \
    GRHydro_DivergenceClean.F90 \