Include files for a new experimental upwinding scheme based on the

characteristics of the evolution equation.


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

2 files changed, 165 insertions, 0 deletions

diff --git a/src/include/metric.h b/src/include/metric.h
new file mode 100644
index 0000000..f3a22ff
--- /dev/null
+++ b/src/include/metric.h
@@ -0,0 +1,19 @@
+! Calculation of inverse three metric, lapse squared and shift
+! $Header$
+
+if ( eh_mask(i,j,k) .ge. 0 ) then
+  alp2 = alp(i,j,k)**2
+
+  tmp1 = gyy(i,j,k)*gzz(i,j,k) - gyz(i,j,k)**2
+  tmp2 = gxz(i,j,k)*gyz(i,j,k) - gxy(i,j,k)*gzz(i,j,k)
+  tmp3 = gxy(i,j,k)*gyz(i,j,k) - gxz(i,j,k)*gyy(i,j,k)
+
+  idetg = one / ( gxx(i,j,k)*tmp1 + gxy(i,j,k)*tmp2 + gxz(i,j,k)*tmp3 )
+
+  g3xx = tmp1 * idetg
+  g3xy = tmp2 * idetg
+  g3xz = tmp3 * idetg
+  g3yy = ( gxx(i,j,k)*gzz(i,j,k) - gxz(i,j,k)**2 ) * idetg
+  g3yz = ( gxy(i,j,k)*gxz(i,j,k) - gxx(i,j,k)*gyz(i,j,k) ) * idetg
+  g3zz = ( gxx(i,j,k)*gyy(i,j,k) - gxy(i,j,k)**2 ) * idetg
+end if
diff --git a/src/include/upwind_characteristic_second2.h b/src/include/upwind_characteristic_second2.h
new file mode 100644
index 0000000..e851859
--- /dev/null
+++ b/src/include/upwind_characteristic_second2.h
@@ -0,0 +1,146 @@
+! Calculation of characteristic upwinded differences except at boundaries.
+! $Header$
+
+! If this is an active cell we have to compute its derivative.
+if ( eh_mask(i,j,k) .ge. 0 ) then
+
+  ! If it is not a boundary cell in the x-direction then calculate as a
+  ! starting point the centered derivative in the x-direction.
+  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) )
+
+!   Otherwise we will have to calculate the appropriate one sided derivative.
+  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
+
+  ! Do the same for the y direction.
+  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
+
+  ! And for the z direction.
+  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) )
+  else
+    dfz(i,j,k) = idz * ( three * f(i,j,k) - four * f(i,j,k-1) + f(i,j,k-2) )
+  end if
+
+! Calculate the characteristic direction using these preliminary
+! derivatives and multiply it with the timestep.
+
+  dfup(1) = g3xx * dfx(i,j,k) + g3xy * dfy(i,j,k) + g3xz * dfz(i,j,k)
+  dfup(2) = g3xy * dfx(i,j,k) + g3yy * dfy(i,j,k) + g3yz * dfz(i,j,k)
+  dfup(3) = g3xz * dfx(i,j,k) + g3yz * dfy(i,j,k) + g3zz * dfz(i,j,k)
+
+  cfactor = one / sqrt ( alp2 * ( dfup(1) * dfx(i,j,k) + &
+                                  dfup(2) * dfy(i,j,k) + &
+                                  dfup(3) * dfz(i,j,k) ) )
+
+  cdx(1) = ( -betax(i,j,k) + alp2 * dfup(1) * cfactor ) * cctk_delta_time
+  cdx(2) = ( -betay(i,j,k) + alp2 * dfup(2) * cfactor ) * cctk_delta_time
+  cdx(3) = ( -betaz(i,j,k) + alp2 * dfup(3) * cfactor ) * cctk_delta_time
+
+  if ( ( .not. btest(eh_mask(i,j,k),0) ) .and. &
+       ( .not. btest(eh_mask(i,j,k),1) ) ) then
+    ! Look at the characteristic direction to figure out in which direction
+    ! the stencil should be chosen.
+    if ( cdx(1) .gt. zero ) then           ! Choose the stencil to the left.
+
+      ! If the neighbour to the left is also an active non-boundary cell its
+      ! safe to assume that i-2, i-1,i and i+1 have good values.
+      ! Otherwise the neigbour must be a boundary cell, so we can only use a
+      ! centered derivative, which we have already computed.
+      if ( eh_mask(i-1,j,k) .eq. 0 ) then
+
+        ! Choose the stencil, that gives the smalles second derivative.
+        ! If it happens to be the centered, do not do anything since we have
+        ! already calculated the centered derivative.
+        if ( f(i-2,j,k) - two * f(i-1,j,k) + f(i,j,k) .le. &
+             f(i-1,j,k) - two * f(i,j,k) + f(i+1,j,k) ) then
+          dfx(i,j,k) = idx * ( three * f(i,j,k) - &
+                               four * f(i-1,j,k) + f(i-2,j,k) )
+        end if
+      end if
+
+    else if ( cdx(1) .lt. zero ) then      ! Choose the stencil to the right.
+
+      ! If the neighbour to the right is also an active non-boundary cell its
+      ! safe to assume that i-1, i,i+1 and i+2 have good values.
+      ! Otherwise the neigbour must be a boundary cell, so we can only use a
+      ! centered derivative, which we have already computed.
+      if ( eh_mask(i+1,j,k) .eq. 0 ) then
+
+        ! Choose the stencil, that gives the smalles second derivative.
+        ! If it happens to be the centered, do not do anything since we have
+        ! already calculated the centered derivative.
+        if ( f(i-1,j,k) - two * f(i,j,k) + f(i+1,j,k) .ge. &
+             f(i,j,k) - two * f(i+1,j,k) + f(i+2,j,k) ) then
+          dfx(i,j,k) = idx * ( -three * f(i,j,k) + &
+                               four * f(i+1,j,k) - f(i+2,j,k) )
+        end if
+      end if
+    end if
+  end if
+  
+  ! Do the same for the y direction.
+  if ( ( .not. btest(eh_mask(i,j,k),2) ) .and. &
+       ( .not. btest(eh_mask(i,j,k),3) ) ) then
+
+    if ( cdx(2) .gt. zero ) then           ! Choose the stencil to the left.
+      if ( eh_mask(i,j-1,k) .eq. 0 ) then
+        if ( f(i,j-2,k) - two * f(i,j-1,k) + f(i,j,k) .le. &
+             f(i,j-1,k) - two * f(i,j,k) + f(i,j+1,k) ) then
+          dfy(i,j,k) = idy * ( three * f(i,j,k) - &
+                               four * f(i,j-1,k) + f(i,j-2,k) )
+        end if
+      end if
+    else if ( cdx(2) .lt. zero ) then      ! Choose the stencil to the left.
+      if ( eh_mask(i,j+1,k) .eq. 0 ) then
+        if ( f(i,j-1,k) - two * f(i,j,k) + f(i,j+1,k) .ge. &
+             f(i,j,k) - two * f(i,j+1,k) + f(i,j+2,k) ) then
+          dfy(i,j,k) = idy * ( -three * f(i,j,k) + &
+                               four * f(i,j+1,k) - f(i,j+2,k) )
+        end if
+      end if
+    end if
+  end if
+
+  ! And for the z direction.
+  if ( ( .not. btest(eh_mask(i,j,k),4) ) .and. &
+       ( .not. btest(eh_mask(i,j,k),5) ) ) then
+
+    if ( cdx(3) .gt. zero ) then
+      if ( eh_mask(i,j,k-1) .eq. 0 ) then
+        if ( f(i,j,k-2) - two * f(i,j,k-1) + f(i,j,k) .le. &
+             f(i,j,k-1) - two * f(i,j,k) + f(i,j,k+1) ) then
+          dfz(i,j,k) = idz * ( three * f(i,j,k) - &
+                               four * f(i,j,k-1) + f(i,j,k-2) )
+        end if
+      end if
+    else if ( cdx(3) .lt. zero ) then
+      if ( eh_mask(i,j,k+1) .eq. 0 ) then
+        if ( f(i,j,k-1) - two * f(i,j,k) + f(i,j,k+1) .ge. &
+             f(i,j,k) - two * f(i,j,k+1) + f(i,j,k+2) ) then
+          dfz(i,j,k) = idz * ( -three * f(i,j,k) + &
+                               four * f(i,j,k+1) - f(i,j,k+2) )
+        end if
+      end if
+    end if
+  end if
+else
+  dfx(i,j,k) = zero
+  dfy(i,j,k) = zero
+  dfz(i,j,k) = zero
+end if