This routine uses now the new include for derivatives within the do-loops.

Added support for mixed centered/upwinded derivatives and for shift
upwinding.


git-svn-id: http://svn.einsteintoolkit.org/cactus/EinsteinAnalysis/EHFinder/trunk@53 2a26948c-0e4f-0410-aee8-f1d3e353619c

1 files changed, 52 insertions, 224 deletions

diff --git a/src/EHFinder_Sources.F90 b/src/EHFinder_Sources.F90
index 9cb287c..80d15dd 100644
--- a/src/EHFinder_Sources.F90
+++ b/src/EHFinder_Sources.F90
@@ -16,15 +16,17 @@ subroutine EHFinder_Sources(CCTK_ARGUMENTS)
   DECLARE_CCTK_FUNCTIONS
 
   CCTK_INT :: i, j, k
-  CCTK_REAL :: idx, idy, idz
+  CCTK_REAL :: idx, idy, idz, mdelta
   CCTK_REAL :: a, b, c
   CCTK_REAL :: g4tt, g4tx, g4ty, g4tz, g4xx, g4xy, g4xz, g4yy, g4yz, g4zz
   CCTK_REAL :: gxxc, gxyc, gxzc, gyyc, gyzc, gzzc, psito4
   CCTK_REAL :: idetg, alp2, tmp1, tmp2, tmp3
   CCTK_REAL :: ratio
   CCTK_REAL, dimension(3) :: maxpos
+  CCTK_REAL :: al, ar, bl, br, cl, cr
+  CCTK_REAL :: alminus, alplus, blminus, blplus, clminus, clplus
+  CCTK_REAL :: arminus, arplus, brminus, brplus, crminus, crplus
 
-!  print*,'Entered Sources'
 ! calculate 1/(2*delta) in each direction
 
 #include "include/physical_part.h"
@@ -33,223 +35,31 @@ subroutine EHFinder_Sources(CCTK_ARGUMENTS)
   idy = half / cctk_delta_space(2)
   idz = half / cctk_delta_space(3)
 
-!  if ( cctk_bbox(1) .eq. 0 ) eh_mask(1:ngx,:,:) = -2
-!  if ( cctk_bbox(2) .eq. 0 ) eh_mask(nx-ngx+1:nx,:,:) = -2
-!  if ( cctk_bbox(3) .eq. 0 ) eh_mask(:,1:ngy,:) = -2
-!  if ( cctk_bbox(4) .eq. 0 ) eh_mask(:,ny-ngy+1:ny,:) = -2
-!  if ( cctk_bbox(5) .eq. 0 ) eh_mask(:,:,1:ngz) = -2
-!  if ( cctk_bbox(6) .eq. 0 ) eh_mask(:,:,nz-ngz+1:nz) = -2
+  mdelta = maxval ( cctk_delta_space )
 
-! calculate one sided derivatives where necessary for x, y and z
-
-!  where ( eh_mask .ge. 0 )
-!    where ( btest(eh_mask,0) )
-!      dfx = idx * ( -three * f + four * cshift(f,1,1) - cshift(f,2,1) )
-!    end where
-!    where ( btest(eh_mask,1) )
-!      dfx = idx * ( three * f - four * cshift(f,-1,1) + cshift(f,-2,1) )
-!    end where
-!    where ( btest(eh_mask,2) )
-!      dfy = idy * ( -three * f + four * cshift(f,1,2) - cshift(f,2,2) )
-!    end where
-!    where ( btest(eh_mask,3) )
-!      dfy = idy * ( three * f - four * cshift(f,-1,2) + cshift(f,-2,2) )
-!    end where
-!    where ( btest(eh_mask,4) )
-!      dfz = idz * ( -three * f + four * cshift(f,1,3) - cshift(f,2,3) )
-!    end where
-!    where ( btest(eh_mask,5) )
-!      dfz = idz * ( three * f - four * cshift(f,-1,3) + cshift(f,-2,3) )
-!    end where
-!  
-!    where ( ( .not. btest(eh_mask,0) ) .and. ( .not. btest(eh_mask,1) ) )
-!      dfx =  idx * ( cshift(f,1,1) - cshift(f,-1,1) )
-!    end where
-!    where ( ( .not. btest(eh_mask,2) ) .and. ( .not. btest(eh_mask,3) ) )
-!      dfy =  idy * ( cshift(f,1,2) - cshift(f,-1,2) )
-!    end where
-!    where ( ( .not. btest(eh_mask,4) ) .and. ( .not. btest(eh_mask,5) ) )
-!      dfz =  idz * ( cshift(f,1,3) - cshift(f,-1,3) )
-!    end where
-!  end where
-
-  do k = kzl, kzr
-    do j = jyl, jyr
-      do i = ixl, ixr
-        if ( eh_mask(i,j,k) .ge. 0 ) then
-          if ( ( .not. btest(eh_mask(i,j,k),0) ) .and. &
-               ( .not. btest(eh_mask(i,j,k),1) ) ) then
-            dfx(i,j,k) = idx * ( f(i+1,j,k) - f(i-1,j,k) )
-          else if ( btest(eh_mask(i,j,k),0) ) then
-            dfx(i,j,k) = idx * ( -three * f(i,j,k) + &
-                                  four * f(i+1,j,k) - f(i+2,j,k) )
-          else
-            dfx(i,j,k) = idx * ( three * f(i,j,k) - &
-                                 four * f(i-1,j,k) + f(i-2,j,k) )
-          end if
-          if ( ( .not. btest(eh_mask(i,j,k),2) ) .and. &
-               ( .not. btest(eh_mask(i,j,k),3) ) ) then
-            dfy(i,j,k) = idy * ( f(i,j+1,k) - f(i,j-1,k) )
-          else if ( btest(eh_mask(i,j,k),2) ) then
-             dfy(i,j,k) = idy * ( -three * f(i,j,k) + &
-                                   four * f(i,j+1,k) - f(i,j+2,k) )
-          else
-             dfy(i,j,k) = idy * ( three * f(i,j,k) - &
-                                  four * f(i,j-1,k) + f(i,j-2,k) )
-          end if
-          if ( ( .not. btest(eh_mask(i,j,k),4) ) .and. &
-               ( .not. btest(eh_mask(i,j,k),5) ) ) then
-            dfz(i,j,k) = idz * ( f(i,j,k+1) - f(i,j,k-1) )
-          else if ( btest(eh_mask(i,j,k),4) ) then
-             dfz(i,j,k) = idz * ( -three * f(i,j,k) + &
-                                   four * f(i,j,k+1) - f(i,j,k+2) )
+  if ( CCTK_EQUALS ( upwind_type, 'intrinsic' ) ) then
+    do k = kzl, kzr
+      do j = jyl, jyr
+        do i = ixl, ixr
+          if ( f(i,j,k) .gt. fmin_bound - three*mdelta ) then
+#include "include/centered_second2.h"
           else
-             dfz(i,j,k) = idz * ( three * f(i,j,k) - &
-                                  four * f(i,j,k-1) + f(i,j,k-2) )
+#include "include/upwind_second2.h"
           end if
-        end if
+        end do
       end do
     end do
-  end do
-
+  end if
     
-!  if ( cctk_bbox(1) .eq. 1 ) then
-!    dfx(1, 1+iy1:ny-iy2, 1+iz1:nz-iz2) = idx * ( &
-!                                  -3.0*f(1, 1+iy1:ny-iy2, 1+iz1:nz-iz2) + &
-!                                   4.0*f(2, 1+iy1:ny-iy2, 1+iz1:nz-iz2) - &
-!                                       f(3, 1+iy1:ny-iy2, 1+iz1:nz-iz2) )
-!  endif
-!  if ( cctk_bbox(2) .eq. 1 ) then
-!    dfx(nx, 1+iy1:ny-iy2, 1+iz1:nz-iz2) = idx * ( &
-!                                   3.0*f(nx, 1+iy1:ny-iy2, 1+iz1:nz-iz2) - &
-!                                   4.0*f(nx-1, 1+iy1:ny-iy2, 1+iz1:nz-iz2) + &
-!                                       f(nx-2, 1+iy1:ny-iy2, 1+iz1:nz-iz2) )
-!  endif
-!
-!! for y
-!
-!  if ( cctk_bbox(3) .eq. 1 ) then
-!    dfy(1+ix1:nx-ix2, 1, 1+iz1:nz-iz2) = idy * ( &
-!                                  -3.0*f(1+ix1:nx-ix2, 1, 1+iz1:nz-iz2) + &
-!                                   4.0*f(1+ix1:nx-ix2, 2, 1+iz1:nz-iz2) - &
-!                                       f(1+ix1:nx-ix2, 3, 1+iz1:nz-iz2) )
-!  endif
-!  if ( cctk_bbox(4) .eq. 1 ) then
-!    dfy(1+ix1:nx-ix2, ny, 1+iz1:nz-iz2) = idy * ( &
-!                                   3.0*f(1+ix1:nx-ix2, ny, 1+iz1:nz-iz2) - &
-!                                   4.0*f(1+ix1:nx-ix2, ny-1, 1+iz1:nz-iz2) + &
-!                                       f(1+ix1:nx-ix2, ny-2, 1+iz1:nz-iz2) )
-!  endif
-!
-!! for z
-!
-!  if ( cctk_bbox(5) .eq. 1 ) then
-!    dfz(1+ix1:nx-ix2, 1+iy1:ny-iy2, 1) = idz * ( &
-!                                  -3.0*f(1+ix1:nx-ix2, 1+iy1:ny-iy2, 1) + &
-!                                   4.0*f(1+ix1:nx-ix2, 1+iy1:ny-iy2, 2) - &
-!                                       f(1+ix1:nx-ix2, 1+iy1:ny-iy2, 3) )
-!  endif
-!  if ( cctk_bbox(6) .eq. 1 ) then
-!    dfz(1+ix1:nx-ix2, 1+iy1:ny-iy2, nz) = idz * ( &
-!                                   3.0*f(1+ix1:nx-ix2, 1+iy1:ny-iy2, nz) - &
-!                                   4.0*f(1+ix1:nx-ix2, 1+iy1:ny-iy2, nz-1) + &
-!                                       f(1+ix1:nx-ix2, 1+iy1:ny-iy2, nz-2) )
-!  endif
-!
-!! calculate normal centered derivatives for interior points
-!
-!  if ( one_sided ) then
-!    do k = 1+iz1, nz-iz2
-!      do j = 1+iy1, ny-iy2
-!        do i = 1+jx1, nx-jx2
-!          if ( f(i,j,k) - f(i-1,j,k) .ne. 0 ) then
-!            ratio = abs ( f(i+1,j,k) - f(i,j,k) ) / ( f(i,j,k) - f(i-1,j,k) )
-!            if ( ratio .lt. 1 ) then
-!              dfx(i,j,k) = idx * ( -3.0*f(i,j,k) + 4.0*f(i+1,j,k) - f(i+2,j,k) )
-!            else
-!              dfx(i,j,k) = idx * ( 3.0*f(i,j,k) - 4.0*f(i-1,j,k) + f(i-2,j,k) )
-!            end if
-!          else
-!            dfx(i,j,k) = idx * ( -3.0*f(i,j,k) + 4.0*f(i+1,j,k) - f(i+2,j,k) )
-!          end if
-!        end do
-!      end do
-!    end do
-!          
-!    do k = 1+iz1, nz-iz2
-!      do j = 1+jy1, ny-jy2
-!        do i = 1+ix1, nx-ix2
-!          if ( f(i,j,k) - f(i,j-1,k) .ne. 0 ) then
-!            ratio = abs ( f(i,j+1,k) - f(i,j,k) ) / ( f(i,j,k) - f(i,j-1,k) )
-!            if ( ratio .lt. 1 ) then
-!              dfy(i,j,k) = idy * ( -3.0*f(i,j,k) + 4.0*f(i,j+1,k) - f(i,j+2,k) )
-!            else
-!              dfy(i,j,k) = idy * ( 3.0*f(i,j,k) - 4.0*f(i,j-1,k) + f(i,j-2,k) )
-!            end if
-!          else
-!            dfy(i,j,k) = idy * ( -3.0*f(i,j,k) + 4.0*f(i,j+1,k) - f(i,j+2,k) )
-!          end if
-!        end do
-!      end do
-!    end do
-!            
-!    do k = 1+jz1, nz-jz2
-!      do j = 1+iy1, ny-iy2
-!        do i = 1+ix1, nx-ix2
-!          if ( f(i,j,k) - f(i,j,k-1) .ne. 0 ) then
-!            ratio = abs ( f(i,j,k+1) - f(i,j,k) ) / ( f(i,j,k) - f(i,j,k-1) )
-!            if ( ratio .lt. 1 ) then
-!              dfz(i,j,k) = idz * ( -3.0*f(i,j,k) + 4.0*f(i,j,k+1) - f(i,j,k+2) )
-!            else
-!              dfz(i,j,k) = idz * ( 3.0*f(i,j,k) - 4.0*f(i,j,k-1) + f(i,j,k-2) )
-!            end if
-!          else
-!            dfz(i,j,k) = idz * ( -3.0*f(i,j,k) + 4.0*f(i,j,k+1) - f(i,j,k+2) )
-!          end if
-!        end do
-!      end do
-!    end do
-!
-!  else
-!
-!    dfx(1+jx1:nx-jx2, 1+iy1:ny-iy2, 1+iz1:nz-iz2) = &
-!                  idx * ( f(2+jx1:nx+1-jx2, 1+iy1:ny-iy2, 1+iz1:nz-iz2) - &
-!                          f(jx1:nx-1-jx2, 1+iy1:ny-iy2, 1+iz1:nz-iz2) )
-!    dfy(1+ix1:nx-ix2, 1+jy1:ny-jy2, 1+iz1:nz-iz2) = &
-!                  idy * ( f(1+ix1:nx-ix2, 2+jy1:ny+1-jy2, 1+iz1:nz-iz2) - &
-!                          f(1+ix1:nx-ix2, jy1:ny-1-jy2, 1+iz1:nz-iz2) )
-!    dfz(1+ix1:nx-ix2, 1+iy1:ny-iy2, 1+jz1:nz-jz2) = &
-!                  idz * ( f(1+ix1:nx-ix2, 1+iy1:ny-iy2, 2+jz1:nz+1-jz2) - &
-!                          f(1+ix1:nx-ix2, 1+iy1:ny-iy2, jz1:nz-1-jz2) )
-!  end if
-
-!  if ( CCTK_MyProc(cctkGH) .eq. 0 ) then
-!    print*,dfx(1,1+iy1:ny-iy2, 1+iz1:nz-iz2)
-!    print*,dfx(nx,1+iy1:ny-iy2, 1+iz1:nz-iz2)
-!    print*
-!    print*,dfy(1+ix1:nx-ix2, 1, 1+iz1:nz-iz2)
-!    print*,dfy(1+ix1:nx-ix2, ny, 1+iz1:nz-iz2)
-!    print*
-!    print*,dfz(1+ix1:nx-ix2, 1+iy1:ny-iy2, 1)
-!    print*,dfz(1+ix1:nx-ix2, 1+iy1:ny-iy2, nz)
-!    print*
-!  end if
-! note I multiply by the square of the lapse in the components of
-! the 4 metric. This does not effect the solution to the quadratic
-! equation, but avoids potential problems of dividing with the square 
-! of the lapse in regions where it is close to zero.
-!
-!        
-!  if ( normalize ) then
-!    dfsq = 1.0 / sqrt(dfx**2+dfy**2+dfz**2)
-!    dfx = dfx*dfsq
-!    dfy = dfy*dfsq
-!    dfz = dfz*dfsq
-!  end if
-
-!  do k = 1+iz1, nz-iz2
-!    do j = 1+iy1, ny-iy2
-!      do i = 1+ix1, nx-ix2
+  if ( CCTK_EQUALS ( upwind_type, 'shift' ) ) then
+    do k = kzl, kzr
+      do j = jyl, jyr
+        do i = ixl, ixr
+#include "include/upwind_shift_second2.h"
+        end do
+      end do
+    end do
+  end if
 
   if ( CCTK_EQUALS ( metric_type, 'physical' ) ) then
     do k = kzl, kzr
@@ -356,12 +166,40 @@ subroutine EHFinder_Sources(CCTK_ARGUMENTS)
                         g4xz * dfx(i,j,k) * dfz(i,j,k) + &
                         g4yz * dfy(i,j,k) * dfz(i,j,k) )
 
-            if ( b .lt. zero ) then
+            if ( b .le. zero ) then
               sf(i,j,k) = ( -b + sqrt(b**2 - four*a*c) ) / ( two*a )
             else
               sf(i,j,k) = two*c / ( -b - sqrt(b**2 - 4*a*c) )
             end if
             sf(i,j,k) = sf(i,j,k)*sign(one,alp(i,j,k))
+
+!            if ( i .eq. 41 .and. j .eq. 20 .and. k .eq. 25 ) then
+!              print*,i,j,k
+!              print*,f(i,j,k),eh_mask(i,j,k)
+!              print*,f(i-1,j,k),f(i+1,j,k),eh_mask(i-1,j,k),eh_mask(i+1,j,k)
+!              print*,f(i,j-1,k),f(i,j+1,k),eh_mask(i,j-1,k),eh_mask(i,j+1,k)
+!              print*,f(i,j,k-1),f(i,j,k+1),eh_mask(i,j,k-1),eh_mask(i,j,k+1)
+!              print*,dfx(i,j,k),dfy(i,j,k),dfz(i,j,k)
+!              print*
+!              print*,x(i,j,k),y(i,j,k),z(i,j,k)
+!              print*,alp(i,j,k)
+!              print*,betax(i,j,k),betay(i,j,k),betaz(i,j,k)
+!              print*,gxx(i,j,k),gxy(i,j,k),gxz(i,j,k)
+!              print*,gyy(i,j,k),gyz(i,j,k),gzz(i,j,k)
+!              print*,psi(i,j,k)
+!              print*
+!              print*,tmp1,tmp2,tmp3
+!              print*,gxxc*tmp1 + gxyc*tmp2 + gxzc*tmp3
+!              print*,idetg
+!              print*
+!              print*,g4xx,g4xy,g4xz
+!              print*,g4yy,g4yz,g4zz
+!              print*
+!              print*,a,b,c
+!              print*,sf(i,j,k)
+!              print*
+!              print*
+!            end if
           else
             sf(i,j,k) = zero
           end if
@@ -370,15 +208,5 @@ subroutine EHFinder_Sources(CCTK_ARGUMENTS)
     end do 
   end if
  
-!  print*,maxval(sf)
-!  print*
-!  sf = zero
-
-!  if ( CCTK_MyProc(cctkGH) .eq. 0 ) then
-!    call CCTK_WARN(0,'stopping')
-!  end if
-
-!  print*,'Exited Sources'
-  
   return
 end subroutine EHFinder_Sources