Preparing for release. Added comments. Changed the way the inverse three metric

is calculated, so I don't have to calculate it both for the characteristic
upwinding and for the rhs evaluation and to make it easier to convert between
conformal metric and physical metric. This requires a set of grid functions
for the inverse 3 metric.


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

3 files changed, 108 insertions, 113 deletions

diff --git a/src/include/metric.h b/src/include/metric.h
index d4890df..97e7a18 100644
--- a/src/include/metric.h
+++ b/src/include/metric.h
@@ -1,19 +1,18 @@
-! Calculation of inverse three metric, lapse squared and shift
+! Calculation of inverse three metric. At this point there is no distinction
+! between the case of a physical or static conformal metric.
 ! $Header$
 
 if ( eh_mask(i,j,k,l) .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
+  g3xx(i,j,k) = tmp1 * idetg
+  g3xy(i,j,k) = tmp2 * idetg
+  g3xz(i,j,k) = tmp3 * idetg
+  g3yy(i,j,k) = ( gxx(i,j,k)*gzz(i,j,k) - gxz(i,j,k)**2 ) * idetg
+  g3yz(i,j,k) = ( gxy(i,j,k)*gxz(i,j,k) - gxx(i,j,k)*gyz(i,j,k) ) * idetg
+  g3zz(i,j,k) = ( 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
index 04b1d06..deec132 100644
--- a/src/include/upwind_characteristic_second2.h
+++ b/src/include/upwind_characteristic_second2.h
@@ -40,9 +40,17 @@ if ( eh_mask(i,j,k,l) .ge. 0 ) then
 ! Calculate the characteristic direction using these preliminary
 ! derivatives and multiply it with the timestep.
 
-  dfup(1) = g3xx * dfx(i,j,k,l) + g3xy * dfy(i,j,k,l) + g3xz * dfz(i,j,k,l)
-  dfup(2) = g3xy * dfx(i,j,k,l) + g3yy * dfy(i,j,k,l) + g3yz * dfz(i,j,k,l)
-  dfup(3) = g3xz * dfx(i,j,k,l) + g3yz * dfy(i,j,k,l) + g3zz * dfz(i,j,k,l)
+  dfup(1) = g3xx(i,j,k) * dfx(i,j,k,l) + &
+            g3xy(i,j,k) * dfy(i,j,k,l) + &
+            g3xz(i,j,k) * dfz(i,j,k,l)
+  dfup(2) = g3xy(i,j,k) * dfx(i,j,k,l) + &
+            g3yy(i,j,k) * dfy(i,j,k,l) + &
+            g3yz(i,j,k) * dfz(i,j,k,l)
+  dfup(3) = g3xz(i,j,k) * dfx(i,j,k,l) + &
+            g3yz(i,j,k) * dfy(i,j,k,l) + &
+            g3zz(i,j,k) * dfz(i,j,k,l)
+
+  alp2 = alp(i,j,k)**2
 
   cfactor = one / sqrt ( alp2 * ( dfup(1) * dfx(i,j,k,l) + &
                                   dfup(2) * dfy(i,j,k,l) + &
diff --git a/src/include/upwind_first.h b/src/include/upwind_first.h
index 1f8c309..366da99 100644
--- a/src/include/upwind_first.h
+++ b/src/include/upwind_first.h
@@ -1,105 +1,93 @@
 ! Calculate first order upwinded differences.
 ! $Header$
 
-idx = one / cctk_delta_space(1)
-idy = one / cctk_delta_space(2)
-idz = one / cctk_delta_space(3)
+if ( eh_mask(i,j,k,l) .ge. 0 ) then
+  if ( ( .not. btest(eh_mask(i,j,k,l),0) ) .and. &
+       ( .not. btest(eh_mask(i,j,k,l),1) ) ) then
+    al = idx * ( f(i,j,k,l) - f(i-1,j,k,l) )
+    ar = idx * ( f(i+1,j,k,l) - f(i,j,k,l) )
+  else if ( btest(eh_mask(i,j,k,l),0) ) then
+    al = zero
+    ar = idx * ( f(i+1,j,k,l) - f(i,j,k,l) )
+  else
+    al = idx * ( f(i,j,k,l) - f(i-1,j,k,l) )
+    ar = zero
+  end if
+  if ( ( .not. btest(eh_mask(i,j,k,l),2) ) .and. &
+       ( .not. btest(eh_mask(i,j,k,l),3) ) ) then
+    bl = idy * ( f(i,j,k,l) - f(i,j-1,k,l) )
+    br = idy * ( f(i,j+1,k,l) - f(i,j,k,l) )
+  else if ( btest(eh_mask(i,j,k,l),2) ) then
+    bl = zero
+    br = idy * ( f(i,j+1,k,l) - f(i,j,k,l) )
+  else
+    bl = idy * ( f(i,j,k,l) - f(i,j-1,k,l) )
+    br = zero
+  end if
+  if ( ( .not. btest(eh_mask(i,j,k,l),4) ) .and. &
+       ( .not. btest(eh_mask(i,j,k,l),5) ) ) then
+    cl = idz * ( f(i,j,k,l) - f(i,j,k-1,l) )
+    cr = idz * ( f(i,j,k+1,l) - f(i,j,k,l) )
+  else if ( btest(eh_mask(i,j,k,l),4) ) then
+    cl = zero
+    cr = idz * ( f(i,j,k+1,l) - f(i,j,k,l) )
+  else
+    cl = idz * ( f(i,j,k,l) - f(i,j,k-1,l) )
+    cr = zero
+  end if
 
-do l = 1, eh_number_level_sets
-  do k = 1, nz
-    do j = 1, ny
-      do i = 1, nx
-        if ( eh_mask(i,j,k,l) .ge. 0 ) then
-          if ( ( .not. btest(eh_mask(i,j,k,l),0) ) .and. &
-               ( .not. btest(eh_mask(i,j,k,l),1) ) ) then
-            al = idx * ( f(i,j,k,l) - f(i-1,j,k,l) )
-            ar = idx * ( f(i+1,j,k,l) - f(i,j,k,l) )
-          else if ( btest(eh_mask(i,j,k,l),0) ) then
-            al = zero
-            ar = idx * ( f(i+1,j,k,l) - f(i,j,k,l) )
-          else
-            al = idx * ( f(i,j,k,l) - f(i-1,j,k,l) )
-            ar = zero
-          end if
-          if ( ( .not. btest(eh_mask(i,j,k,l),2) ) .and. &
-               ( .not. btest(eh_mask(i,j,k,l),3) ) ) then
-            bl = idy * ( f(i,j,k,l) - f(i,j-1,k,l) )
-            br = idy * ( f(i,j+1,k,l) - f(i,j,k,l) )
-          else if ( btest(eh_mask(i,j,k,l),2) ) then
-            bl = zero
-            br = idy * ( f(i,j+1,k,l) - f(i,j,k,l) )
-          else
-            bl = idy * ( f(i,j,k,l) - f(i,j-1,k,l) )
-            br = zero
-          end if
-          if ( ( .not. btest(eh_mask(i,j,k,l),4) ) .and. &
-               ( .not. btest(eh_mask(i,j,k,l),5) ) ) then
-            cl = idz * ( f(i,j,k,l) - f(i,j,k-1,l) )
-            cr = idz * ( f(i,j,k+1,l) - f(i,j,k,l) )
-          else if ( btest(eh_mask(i,j,k,l),4) ) then
-            cl = zero
-            cr = idz * ( f(i,j,k+1,l) - f(i,j,k,l) )
-          else
-            cl = idz * ( f(i,j,k,l) - f(i,j,k-1,l) )
-            cr = zero
-          end if
-  
-          alminus = half*(al-abs(al))
-          alplus = half*(al+abs(al))
-          arminus = half*(ar-abs(ar))
-          arplus = half*(ar+abs(ar))
-  
-          blminus = half*(bl-abs(bl))
-          blplus = half*(bl+abs(bl))
-          brminus = half*(br-abs(br))
-          brplus = half*(br+abs(br))
-  
-          clminus = half*(cl-abs(cl))
-          clplus = half*(cl+abs(cl))
-          crminus = half*(cr-abs(cr))
-          crplus = half*(cr+abs(cr))
-  
-          if ( f(i,j,k,l) .gt. 0 ) then
-            if ( abs(alplus) .gt. abs(arminus) ) then
-              dfx(i,j,k,l) = alplus
-            else
-              dfx(i,j,k,l) = arminus
-            end if
-            if ( abs(blplus) .gt. abs(brminus) ) then
-              dfy(i,j,k,l) = blplus
-            else
-              dfy(i,j,k,l) = brminus
-            end if
-            if ( abs(clplus) .gt. abs(crminus) ) then
-              dfz(i,j,k,l) = clplus
-            else
-              dfz(i,j,k,l) = crminus
-            end if
-          endif
-  
-          if ( f(i,j,k,l) .lt. 0 ) then
-            if ( abs(alminus) .gt. abs(arplus) ) then
-              dfx(i,j,k,l) = alminus
-            else
-              dfx(i,j,k,l) = arplus
-            end if
-            if ( abs(blminus) .gt. abs(brplus) ) then
-              dfy(i,j,k,l) = blminus
-            else
-              dfy(i,j,k,l) = brplus
-            end if
-            if ( abs(clminus) .gt. abs(crplus) ) then
-              dfz(i,j,k,l) = clminus
-            else
-              dfz(i,j,k,l) = crplus
-            end if
-          end if
-        else
-          dfx(i,j,k,l) = zero
-          dfy(i,j,k,l) = zero
-          dfz(i,j,k,l) = zero
-        end if
-      end do
-    end do
-  end do
-end do
+  alminus = half*(al-abs(al))
+  alplus = half*(al+abs(al))
+  arminus = half*(ar-abs(ar))
+  arplus = half*(ar+abs(ar))
+
+  blminus = half*(bl-abs(bl))
+  blplus = half*(bl+abs(bl))
+  brminus = half*(br-abs(br))
+  brplus = half*(br+abs(br))
+
+  clminus = half*(cl-abs(cl))
+  clplus = half*(cl+abs(cl))
+  crminus = half*(cr-abs(cr))
+  crplus = half*(cr+abs(cr))
+
+  if ( f(i,j,k,l) .gt. 0 ) then
+    if ( abs(alplus) .gt. abs(arminus) ) then
+      dfx(i,j,k,l) = alplus
+    else
+      dfx(i,j,k,l) = arminus
+    end if
+    if ( abs(blplus) .gt. abs(brminus) ) then
+      dfy(i,j,k,l) = blplus
+    else
+      dfy(i,j,k,l) = brminus
+    end if
+    if ( abs(clplus) .gt. abs(crminus) ) then
+      dfz(i,j,k,l) = clplus
+    else
+      dfz(i,j,k,l) = crminus
+    end if
+  endif
+
+  if ( f(i,j,k,l) .lt. 0 ) then
+    if ( abs(alminus) .gt. abs(arplus) ) then
+      dfx(i,j,k,l) = alminus
+    else
+      dfx(i,j,k,l) = arplus
+    end if
+    if ( abs(blminus) .gt. abs(brplus) ) then
+      dfy(i,j,k,l) = blminus
+    else
+      dfy(i,j,k,l) = brplus
+    end if
+    if ( abs(clminus) .gt. abs(crplus) ) then
+      dfz(i,j,k,l) = clminus
+    else
+      dfz(i,j,k,l) = crplus
+    end if
+  end if
+else
+  dfx(i,j,k,l) = zero
+  dfy(i,j,k,l) = zero
+  dfz(i,j,k,l) = zero
+end if