Added support for a conformal metric. Added support for the new

interpolator parameters. Cleaned up the code a bit.


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

1 files changed, 150 insertions, 68 deletions

diff --git a/src/EHFinder_Generator_Sources.F90 b/src/EHFinder_Generator_Sources.F90
index 924588f..95cb22e 100644
--- a/src/EHFinder_Generator_Sources.F90
+++ b/src/EHFinder_Generator_Sources.F90
@@ -17,6 +17,11 @@ subroutine EHFinder_Generator_Sources(CCTK_ARGUMENTS)
 
   CCTK_INT :: i
   CCTK_INT :: interp_handle, table_handle, status, coord_system_handle
+
+  character(len=200) :: gen_interp
+  CCTK_INT :: gen_interp_len
+  character(len=7) :: gen_order
+
   CCTK_INT, dimension(1) :: lsh
   CCTK_POINTER, dimension(3) :: interp_coords
   CCTK_POINTER, dimension(14) :: out_arrays
@@ -28,78 +33,102 @@ subroutine EHFinder_Generator_Sources(CCTK_ARGUMENTS)
                                                       0, 0, 0, 0, 0, &
                                                       1, 2, 3, 0 /)
   CCTK_INT, dimension(14) :: out_types
-  CCTK_REAL :: alp2, dfux, dfuy, dfuz, factor
+  CCTK_REAL :: alp2, psi4, dfux, dfuy, dfuz, factor
   CCTK_REAL :: idetg, guxx, guxy, guxz, guyy, guyz, guzz
 
   out_types = CCTK_VARIABLE_REAL
 
-  if ( CCTK_EQUALS ( metric_type, 'physical' ) ) then
-    call CCTK_GrouplshGN ( status, cctkGH, 1, lsh, "ehfinder::generators" )
-    if ( status .lt. 0 ) then
-      call CCTK_WARN ( 0, "cannot get local size for surface arrays" )
-    end if
-!    print*,lsh 
+  ! Convert the generator_interpolator string parameter to a Fortran string.
+  call CCTK_FortranString ( gen_interp_len, generator_interpolator, &
+                                            gen_interp )
 
-    call CCTK_InterpHandle ( interp_handle, &
-                                     "Lagrange polynomial interpolation" )
-    if ( interp_handle .lt. 0 ) then
-      call CCTK_WARN( 0, "Cannot get handle for interpolation. Forgot to activate an implementation providing interpolation operators??" )
-    end if
-    ! For now order=2 is hard wired.
-    call Util_TableCreateFromString ( table_handle, "order=3" )
-    if ( table_handle .lt. 0 ) then
-      call CCTK_WARN( 0, "Cannot create parameter table for interpolator" )
-    end if
+  ! Get the corresponding interpolator handle.
+  call CCTK_InterpHandle ( interp_handle, gen_interp )
 
-    ! Get the 3D coordinate system handle.
-    call CCTK_CoordSystemHandle ( coord_system_handle, "cart3d" )
-    if ( coord_system_handle .lt. 0) then
-      call CCTK_WARN( 0, "Cannot get handle for cart3d coordinate system.  Forgot to activate an implementation providing coordinates ??" )
-    endif
+  if ( interp_handle .lt. 0 ) then
+    call CCTK_WARN( 0, "Cannot get handle for interpolation. Forgot to activate an implementation providing interpolation operators??" )
+  end if
+
+  ! Convert the interpolation order parameter to a Fortran string to be placed
+  ! in the interpolator table. Note that the order is assumed to contain only
+  ! 1 digit.
+  write(gen_order,'(a6,i1)') 'order=',generator_interpolation_order
+
+  ! Create the table directly from the string.
+  call Util_TableCreateFromString ( table_handle, gen_order )
+  if ( table_handle .lt. 0 ) then
+    call CCTK_WARN( 0, "Cannot create parameter table for interpolator" )
+  end if
 
-    interp_coords(1) = CCTK_PointerTo(xg)
-    interp_coords(2) = CCTK_PointerTo(yg)
-    interp_coords(3) = CCTK_PointerTo(zg)
-
-    out_arrays(1) = CCTK_PointerTo(alpg)
-    out_arrays(2) = CCTK_PointerTo(betaxg)
-    out_arrays(3) = CCTK_PointerTo(betayg)
-    out_arrays(4) = CCTK_PointerTo(betazg)
-    out_arrays(5) = CCTK_PointerTo(gxxg)
-    out_arrays(6) = CCTK_PointerTo(gxyg)
-    out_arrays(7) = CCTK_PointerTo(gxzg)
-    out_arrays(8) = CCTK_PointerTo(gyyg)
-    out_arrays(9) = CCTK_PointerTo(gyzg)
-    out_arrays(10) = CCTK_PointerTo(gzzg)
-    out_arrays(11) = CCTK_PointerTo(dfxg)
-    out_arrays(12) = CCTK_PointerTo(dfyg)
-    out_arrays(13) = CCTK_PointerTo(dfzg)
-
-    call CCTK_VarIndex ( in_arrays(1), "admbase::alp" )
-    call CCTK_VarIndex ( in_arrays(2), "admbase::betax" )
-    call CCTK_VarIndex ( in_arrays(3), "admbase::betay" )
-    call CCTK_VarIndex ( in_arrays(4), "admbase::betaz" )
-    call CCTK_VarIndex ( in_arrays(5), "admbase::gxx" )
-    call CCTK_VarIndex ( in_arrays(6), "admbase::gxy" )
-    call CCTK_VarIndex ( in_arrays(7), "admbase::gxz" )
-    call CCTK_VarIndex ( in_arrays(8), "admbase::gyy" )
-    call CCTK_VarIndex ( in_arrays(9), "admbase::gyz" )
-    call CCTK_VarIndex ( in_arrays(10), "admbase::gzz" )
-    call CCTK_VarIndex ( in_arrays(11), "ehfinder::f" )
+  ! Get the 3D coordinate system handle.
+  call CCTK_CoordSystemHandle ( coord_system_handle, "cart3d" )
+  if ( coord_system_handle .lt. 0) then
+    call CCTK_WARN( 0, "Cannot get handle for cart3d coordinate system.  Forgot to activate an implementation providing coordinates ??" )
+  endif
+
+  ! Find out how many interpolation points are located on this processor.
+  call CCTK_GrouplshGN ( status, cctkGH, 1, lsh, "ehfinder::generators" )
+  if ( status .lt. 0 ) then
+    call CCTK_WARN ( 0, "cannot get local size for surface arrays" )
+  end if
+
+  ! Set the pointers to the points to be interpolated to.
+  interp_coords(1) = CCTK_PointerTo(xg)
+  interp_coords(2) = CCTK_PointerTo(yg)
+  interp_coords(3) = CCTK_PointerTo(zg)
+
+  ! Set the pointers to the output arrays.
+  out_arrays(1) = CCTK_PointerTo(alpg)
+  out_arrays(2) = CCTK_PointerTo(betaxg)
+  out_arrays(3) = CCTK_PointerTo(betayg)
+  out_arrays(4) = CCTK_PointerTo(betazg)
+  out_arrays(5) = CCTK_PointerTo(gxxg)
+  out_arrays(6) = CCTK_PointerTo(gxyg)
+  out_arrays(7) = CCTK_PointerTo(gxzg)
+  out_arrays(8) = CCTK_PointerTo(gyyg)
+  out_arrays(9) = CCTK_PointerTo(gyzg)
+  out_arrays(10) = CCTK_PointerTo(gzzg)
+  out_arrays(11) = CCTK_PointerTo(dfxg)
+  out_arrays(12) = CCTK_PointerTo(dfyg)
+  out_arrays(13) = CCTK_PointerTo(dfzg)
+  out_arrays(14) = CCTK_PointerTo(psig)
+
+  ! Set the indices to the input grid functions.
+  call CCTK_VarIndex ( in_arrays(1), "admbase::alp" )
+  call CCTK_VarIndex ( in_arrays(2), "admbase::betax" )
+  call CCTK_VarIndex ( in_arrays(3), "admbase::betay" )
+  call CCTK_VarIndex ( in_arrays(4), "admbase::betaz" )
+  call CCTK_VarIndex ( in_arrays(5), "admbase::gxx" )
+  call CCTK_VarIndex ( in_arrays(6), "admbase::gxy" )
+  call CCTK_VarIndex ( in_arrays(7), "admbase::gxz" )
+  call CCTK_VarIndex ( in_arrays(8), "admbase::gyy" )
+  call CCTK_VarIndex ( in_arrays(9), "admbase::gyz" )
+  call CCTK_VarIndex ( in_arrays(10), "admbase::gzz" )
+  call CCTK_VarIndex ( in_arrays(11), "ehfinder::f" )
+  call CCTK_VarIndex ( in_arrays(12), "staticconformal::psi" )
+
+  ! Check the metric type. At present physical and static_conformal are
+  ! supported.
+  if ( CCTK_EQUALS ( metric_type, 'physical' ) ) then
 
+    ! Set the operand indices table entry, corresponding
+    ! to interpolation of admbase::alp (1), admbase::shift (3),
+    ! admbase::metric(6) and the first derivatives of ehfinder::f (3) for
+    ! a total of 13 output arrays.
     call Util_TableSetIntArray ( status, table_handle, 13, &
                                op_indices(1:13), "operand_indices" )
     if ( status .lt. 0 ) then
       call CCTK_WARN ( 0, "Cannot set operand indices array in parameter table" )
     endif
 
+    ! Set the corresponding table entry for the operation codes.
     call Util_TableSetIntArray ( status, table_handle, 13, &
                                op_codes(1:13), "operation_codes" )
     if ( status .lt. 0 ) then
       call CCTK_WARN ( 0, "Cannot set operation codes array in parameter table" )
     endif
 
-
+    ! Call the interpolator.
     call CCTK_InterpGridArrays ( status, cctkGH, 3, interp_handle, &
                                  table_handle, coord_system_handle, &
                                  lsh(1), CCTK_VARIABLE_REAL, &
@@ -109,9 +138,15 @@ subroutine EHFinder_Generator_Sources(CCTK_ARGUMENTS)
     if ( status .lt. 0 ) then
       call CCTK_INFO ( 'Interpolation failed.' )
     end if
+
+    ! For each point on this processor calculate the right hand side of the
+    ! characteristic evolution equation.
     do i = 1, lsh(1)
+
+      ! calculate the square of the lapse.
       alp2 = alpg(i)**2
 
+      ! Calculate the inverse of the 3-metric.
       guxx = gyyg(i) * gzzg(i) - gyzg(i)**2
       guxy = gxzg(i) * gyzg(i) - gxyg(i) * gzzg(i)
       guxz = gxyg(i) * gyzg(i) - gxzg(i) * gyyg(i)
@@ -126,6 +161,69 @@ subroutine EHFinder_Generator_Sources(CCTK_ARGUMENTS)
       guyz = ( gxyg(i) * gxzg(i) - gxxg(i) * gyzg(i) ) * idetg
       guzz = ( gxxg(i) * gyyg(i) - gxyg(i)**2 ) * idetg
       
+      ! Raise the index of the partial derivatives of f.
+      dfux = guxx * dfxg(i) + guxy * dfyg(i) + guxz * dfzg(i)
+      dfuy = guxy * dfxg(i) + guyy * dfyg(i) + guyz * dfzg(i)
+      dfuz = guxz * dfxg(i) + guyz * dfyg(i) + guzz * dfzg(i)
+
+      ! Calculate the overall multiplication factor.
+      factor = alp2 / sqrt ( alp2 * ( dfux * dfxg(i) + &
+                                      dfuy * dfyg(i) + &
+                                      dfuz * dfzg(i) ) )
+
+      ! Finally obtain dx^i/dt.
+      dxg(i) = - betaxg(i) + factor * dfux
+      dyg(i) = - betayg(i) + factor * dfuy
+      dzg(i) = - betazg(i) + factor * dfuz
+    end do
+  else if ( CCTK_EQUALS ( metric_type, 'static conformal' ) ) then
+    call Util_TableSetIntArray ( status, table_handle, 14, &
+                               op_indices, "operand_indices" )
+    if ( status .lt. 0 ) then
+      call CCTK_WARN ( 0, "Cannot set operand indices array in parameter table" )
+    endif
+
+    call Util_TableSetIntArray ( status, table_handle, 14, &
+                               op_codes, "operation_codes" )
+    if ( status .lt. 0 ) then
+      call CCTK_WARN ( 0, "Cannot set operation codes array in parameter table" )
+    endif
+
+    call CCTK_InterpGridArrays ( status, cctkGH, 3, interp_handle, &
+                                 table_handle, coord_system_handle, &
+                                 lsh(1), CCTK_VARIABLE_REAL, &
+                                 interp_coords, 12, in_arrays, &
+                                 14, out_types, out_arrays )
+
+    if ( status .lt. 0 ) then
+      call CCTK_INFO ( 'Interpolation failed.' )
+    end if
+
+    do i = 1, lsh(1)
+      alp2 = alpg(i)**2
+
+!     The inverse of psi^4
+      psi4 = one / psig(i)**4
+
+      guxx = gyyg(i) * gzzg(i) - gyzg(i)**2
+      guxy = gxzg(i) * gyzg(i) - gxyg(i) * gzzg(i)
+      guxz = gxyg(i) * gyzg(i) - gxzg(i) * gyyg(i)
+
+!     The determinant divided by psi^4.
+      idetg = psi4 / ( gxxg(i) * guxx + &
+                       gxyg(i) * guxy + &
+                       gxzg(i) * guxz )
+
+!     The inverse metric. Since the determinant is already divided
+!     by psi^4, this gives the inverse of the physical metric.
+      guxx = idetg * guxx
+      guxy = idetg * guxy
+      guxz = idetg * guxz
+
+      guyy = ( gxxg(i) * gzzg(i) - gxzg(i)**2 ) * idetg
+      guyz = ( gxyg(i) * gxzg(i) - gxxg(i) * gyzg(i) ) * idetg
+      guzz = ( gxxg(i) * gyyg(i) - gxyg(i)**2 ) * idetg
+      
       dfux = guxx * dfxg(i) + guxy * dfyg(i) + guxz * dfzg(i)
       dfuy = guxy * dfxg(i) + guyy * dfyg(i) + guyz * dfzg(i)
       dfuz = guxz * dfxg(i) + guyz * dfyg(i) + guzz * dfzg(i)
@@ -135,22 +233,6 @@ subroutine EHFinder_Generator_Sources(CCTK_ARGUMENTS)
       dxg(i) = - betaxg(i) + factor * dfux
       dyg(i) = - betayg(i) + factor * dfuy
       dzg(i) = - betazg(i) + factor * dfuz
-!      print*, alpg
-!      print*
-!      print*, betaxg, betayg, betazg
-!      print*
-!      print*, gxxg, gxyg, gxzg, gyyg, gyzg, gzzg
-!      print*
-!      print*, dfxg, dfyg, dfzg
-!      print*
-!      print*, alp2
-!      print*
-!      print*, dfux, dfuy, dfuz
-!      print*
-!      print*, factor
-!      print*
-!      print*, dxg(i), dyg(i), dzg(i)
-!      stop
     end do
   end if
   return