GRHydro: Add basic vector potential support

   Basic cell-centered, algebraic gauge vector potential
   method with place-holders for lorenz gauge.
   Initial Avec constrained to poloidal at the moment.

From: Tanja Bode <tanja.bode@physics.gatech.edu>

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

1 files changed, 147 insertions, 0 deletions

diff --git a/src/GRHydro_BvecfromAvec.F90 b/src/GRHydro_BvecfromAvec.F90
new file mode 100644
index 0000000..ea0a58d
--- /dev/null
+++ b/src/GRHydro_BvecfromAvec.F90
@@ -0,0 +1,147 @@
+  /*@@
+   @file      GRHydro_BvecfromAvec 
+   @date      Aug 31, 2010
+   @author    Tanja Bode
+   @desc 
+   Calculate B^i (at cell center) from Avec (at edges)
+   @enddesc 
+ @@*/
+
+#include "cctk.h"
+#include "cctk_Arguments.h"
+#include "cctk_Parameters.h"
+#include "cctk_Functions.h"
+#include "GRHydro_Macros.h"
+
+#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 Avecx(i,j,k) Avec(i,j,k,1)
+#define Avecy(i,j,k) Avec(i,j,k,2)
+#define Avecz(i,j,k) Avec(i,j,k,3)
+
+! Second order f.d. for volume-centered quantity from volume-centered quantity
+#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. for volume-centered quantity from volume-centered quantity
+#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)
+
+ /*@@
+   @routine   GRHydro_BvecfromAvec
+   @date      Oct 24
+   @author    Tanja Bode
+   @desc 
+   Calculate B^i (at cell center) from Avec (at edges)
+   @enddesc 
+   @calls     
+   @calledby   
+   @history 
+ 
+   @endhistory 
+
+@@*/
+
+subroutine GRHydro_BvecfromAvec(CCTK_ARGUMENTS)
+
+  implicit none
+
+  CCTK_REAL :: Az_y, Ay_z, Ax_z, Az_x, Ay_x, Ax_y
+  CCTK_REAL :: det, isdet
+  CCTK_REAL :: dx, dy, dz, idx, idy, idz
+  CCTK_INT  :: i, j, k, nx, ny, nz, GRHydro_UseGeneralCoordinates
+
+  logical, allocatable, dimension (:,:,:) :: force_spatial_second_order
+
+  DECLARE_CCTK_ARGUMENTS
+  DECLARE_CCTK_PARAMETERS
+  !!DECLARE_CCTK_FUNCTIONS
+
+  if (GRHydro_UseGeneralCoordinates(cctkGH).ne.0) then
+     call CCTK_WARN(0,"Bvec from Avec only in Cartesian at the moment."); 
+  end if
+
+  nx = cctk_lsh(1)
+  ny = cctk_lsh(2)
+  nz = cctk_lsh(3)
+  dx = CCTK_DELTA_SPACE(1)
+  dy = CCTK_DELTA_SPACE(2)
+  dz = CCTK_DELTA_SPACE(3)
+  idx = 1.d0/dx
+  idy = 1.d0/dy
+  idz = 1.d0/dz
+ 
+  allocate (force_spatial_second_order(nx,ny,nz))
+  force_spatial_second_order = .FALSE.
+  
+  if (spatial_order > 2) then
+    !$OMP PARALLEL DO PRIVATE(i, j, k)
+    do k = 1 + GRHydro_stencil, nz - GRHydro_stencil
+      do j = 1 + GRHydro_stencil, ny - GRHydro_stencil
+        do i = 1 + GRHydro_stencil, nx - GRHydro_stencil
+          if ((i < 3).or.(i > cctk_lsh(1) - 2).or. &
+               (j < 3).or.(j > cctk_lsh(2) - 2).or. &
+               (k < 3).or.(k > cctk_lsh(3) - 2) ) then
+            force_spatial_second_order(i,j,k) = .TRUE.
+          else if ( use_mask > 0 ) then
+            if (minval(emask(i-2:i+2,j-2:j+2,k-2:k+2)) < 0.75d0) then
+              force_spatial_second_order(i,j,k) = .TRUE.
+            end if
+          end if
+        end do
+      end do
+    end do
+    !$OMP END PARALLEL DO
+  end if
+  
+ 
+  call CCTK_WARN(1,"Bvec from Avec start."); 
+  !$OMP PARALLEL DO PRIVATE( det, isdet, local_spatial_order, &
+  !$OMP                      Az_y, Ay_z, Ax_z, Az_x, Ay_x, Ax_y)
+  do k = GRHydro_stencil, cctk_lsh(3)-GRHydro_stencil+1
+    do j = GRHydro_stencil, cctk_lsh(2)-GRHydro_stencil+1
+      do i = GRHydro_stencil, cctk_lsh(1)-GRHydro_stencil+1
+
+          local_spatial_order = spatial_order
+          if (force_spatial_second_order(i,j,k)) then
+            local_spatial_order = 2
+          end if
+
+          if (local_spatial_order .eq. 2) then
+             Az_y = DIFF_Y_2(Avecz)
+             Ay_z = DIFF_Z_2(Avecy)
+             Ax_z = DIFF_Z_2(Avecx)
+             Az_x = DIFF_X_2(Avecz)
+             Ay_x = DIFF_X_2(Avecy)
+             Ax_y = DIFF_Y_2(Avecx)
+          else
+             Az_y = DIFF_Y_4(Avecz)
+             Ay_z = DIFF_Z_4(Avecy)
+             Ax_z = DIFF_Z_4(Avecx)
+             Az_x = DIFF_X_4(Avecz)
+             Ay_x = DIFF_X_4(Avecy)
+             Ax_y = DIFF_Y_4(Avecx)
+          end if
+
+          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))
+          isdet = 1.d0/sqrt(det)
+
+          Bvecx(i,j,k) = isdet*( Az_y - Ay_z )
+          Bvecy(i,j,k) = isdet*( Ax_z - Az_x )
+          Bvecz(i,j,k) = isdet*( Ay_x - Ax_y )
+
+      end do
+    end do
+  end do
+  !$OMP END PARALLEL DO
+  call CCTK_WARN(1,"Bvec from Avec end."); 
+
+
+end subroutine GRHydro_BvecfromAvec
+