aboutsummaryrefslogtreecommitdiff
path: root/src/AHFinder_int.F
diff options
context:
space:
mode:
Diffstat (limited to 'src/AHFinder_int.F')
-rw-r--r--src/AHFinder_int.F177
1 files changed, 91 insertions, 86 deletions
diff --git a/src/AHFinder_int.F b/src/AHFinder_int.F
index b3feabe..060834a 100644
--- a/src/AHFinder_int.F
+++ b/src/AHFinder_int.F
@@ -3,7 +3,7 @@
@file AHFinder_int.F
@date April 1998
@author Miguel Alcubierre
- @desc
+ @desc
Find surface integrals. The integrals are done
in parallel, but the number of processors is
assumed to be either the square of an integer
@@ -34,10 +34,16 @@
integer auxi
integer ierror
- integer interp_handle,coord_system_handle,reduction_handle
+ integer interp_handle,coord_system_handle
+ integer reduction_handle,param_table_handle
+ character(30) options_string
CCTK_INT red_tmp
- integer, dimension(9) :: in_array_indices
+ CCTK_POINTER, dimension(3) :: interp_coords
+ CCTK_INT, dimension(9) :: in_array_indices
+ CCTK_POINTER, dimension(9) :: out_arrays
+ CCTK_INT, dimension(9) :: out_array_type_codes
+
CCTK_REAL LEGEN
CCTK_REAL xp,yp,zp,rp
@@ -50,7 +56,6 @@
CCTK_REAL sina,cosa,intw,grad,sigma,h0
CCTK_REAL theta,phi
CCTK_REAL aux,aux1,aux2
- INTEGER num_arrays
CCTK_REAL tempv(0:lmax),tempm(lmax,lmax)
@@ -166,7 +171,7 @@
call CCTK_ReductionArrayHandle(reduction_handle,"sum")
- if (reduction_handle.lt.0) then
+ if (reduction_handle.lt.0) then
call CCTK_WARN(1,"Cannot get handle for sum reduction ! Forgot to activate an implementation providing reduction operators ??")
end if
@@ -577,7 +582,7 @@
interror1 = red_tmp
call CCTK_ReduceLocalScalar(ierror,cctkGH,-1,reduction_handle,
- . interror2,red_tmp,CCTK_VARIABLE_INT)
+ . interror2,red_tmp,CCTK_VARIABLE_INT)
if (ierror.ne.0) then
call CCTK_WARN(1,"Reduction of norm failed!")
end if
@@ -608,29 +613,32 @@
npoints = (l_ntheta+1)*(l_nphi+1)
-! Interpolator and coordinate system handle.
-
+! parameter, local interpolator, and coordinate system handle.
+ param_table_handle = -1
interp_handle = -1
coord_system_handle = -1
- if (interpolation_order .eq. 1) then
- call CCTK_InterpHandle (interp_handle, "first-order uniform cartesian")
- else if (interpolation_order .eq. 2) then
- call CCTK_InterpHandle (interp_handle, "second-order uniform cartesian")
- else if (interpolation_order .eq. 3) then
- call CCTK_InterpHandle (interp_handle, "third-order uniform cartesian")
+ options_string = "order = " // char(ichar('0') + interpolation_order)
+ call Util_TableCreateFromString (param_table_handle, options_string)
+ if (param_table_handle .lt. 0) then
+ call CCTK_WARN(0,"Cannot create parameter table for interpolator")
endif
+ call CCTK_InterpHandle (interp_handle,"Lagrange polynomial interpolation")
if (interp_handle .lt. 0) then
- call CCTK_WARN(0,"Cannot get handle for cart3d coordinate system ! Forgot to activate an implementation providing coordinates ??")
+ call CCTK_WARN(0,"Cannot get handle for interpolation ! Forgot to activate an implementation providing interpolation operators ??")
endif
call CCTK_CoordSystemHandle (coord_system_handle, "cart3d")
-
if (coord_system_handle .lt. 0) then
- call CCTK_WARN(0,"Cannot get handle for interpolation ! Forgot to activate an implementation providing interpolation operators ??")
+ call CCTK_WARN(0,"Cannot get handle for cart3d coordinate system ! Forgot to activate an implementation providing coordinates ??")
endif
+! fill in the input/output arrays for the interpolator
+ interp_coords(1) = CCTK_PointerTo(xa)
+ interp_coords(2) = CCTK_PointerTo(ya)
+ interp_coords(3) = CCTK_PointerTo(za)
+
call CCTK_VarIndex (in_array_indices(1), "admbase::gxx")
call CCTK_VarIndex (in_array_indices(2), "admbase::gyy")
call CCTK_VarIndex (in_array_indices(3), "admbase::gzz")
@@ -641,69 +649,66 @@
call CCTK_VarIndex (in_array_indices(8), "ahfinder::ahmask")
call CCTK_VarIndex (in_array_indices(9), "ahfinder::ahfgradn")
+ out_arrays(1) = CCTK_PointerTo(txx)
+ out_arrays(2) = CCTK_PointerTo(tyy)
+ out_arrays(3) = CCTK_PointerTo(tzz)
+ out_arrays(4) = CCTK_PointerTo(txy)
+ out_arrays(5) = CCTK_PointerTo(txz)
+ out_arrays(6) = CCTK_PointerTo(tyz)
+ out_arrays(7) = CCTK_PointerTo(exp)
+ out_arrays(8) = CCTK_PointerTo(intmask)
+ out_arrays(9) = CCTK_PointerTo(gradn)
+
+ out_array_type_codes = CCTK_VARIABLE_REAL
+
+
! For minimization we need the interpolated metric,
! the expansion and the mask.
if (.not.flow) then
- num_arrays = 8
- call CCTK_InterpGV (ierror, cctkGH, interp_handle, coord_system_handle,
- . npoints, num_arrays, num_arrays,
- . xa, ya, za,
- . CCTK_VARIABLE_REAL, CCTK_VARIABLE_REAL,
- . CCTK_VARIABLE_REAL,
- . in_array_indices(1), in_array_indices(2),
- . in_array_indices(3), in_array_indices(4),
- . in_array_indices(5), in_array_indices(6),
- . in_array_indices(7), in_array_indices(8),
- . txx, tyy, tzz, txy, txz, tyz, exp, intmask,
- . CCTK_VARIABLE_REAL, CCTK_VARIABLE_REAL,
- . CCTK_VARIABLE_REAL, CCTK_VARIABLE_REAL,
- . CCTK_VARIABLE_REAL, CCTK_VARIABLE_REAL,
- . CCTK_VARIABLE_REAL, CCTK_VARIABLE_REAL)
+ call CCTK_InterpGridArrays (ierror, cctkGH, 3, interp_handle,
+ . param_table_handle, coord_system_handle,
+ . npoints, CCTK_VARIABLE_REAL, interp_coords,
+ . 8, in_array_indices,
+ . 8, out_array_type_codes, out_arrays)
+ if (ierror < 0) then
+ call CCTK_WARN (1, "interpolator call returned an error code");
+ endif
+
! For N flow, we need the interpolated metric, the expansion,
! the norm of the gradient of the horizon function, and the mask.
else if (nw.ne.zero) then
- num_arrays = 9
- call CCTK_InterpGV (ierror, cctkGH, interp_handle, coord_system_handle,
- . npoints, num_arrays, num_arrays,
- . xa, ya, za,
- . CCTK_VARIABLE_REAL, CCTK_VARIABLE_REAL,
- . CCTK_VARIABLE_REAL,
- . in_array_indices(1), in_array_indices(2),
- . in_array_indices(3), in_array_indices(4),
- . in_array_indices(5), in_array_indices(6),
- . in_array_indices(7), in_array_indices(8),
- . in_array_indices(9),
- . txx, tyy, tzz, txy, txz, tyz, exp, intmask, gradn,
- . CCTK_VARIABLE_REAL, CCTK_VARIABLE_REAL,
- . CCTK_VARIABLE_REAL, CCTK_VARIABLE_REAL,
- . CCTK_VARIABLE_REAL, CCTK_VARIABLE_REAL,
- . CCTK_VARIABLE_REAL, CCTK_VARIABLE_REAL,
- . CCTK_VARIABLE_REAL)
+ call CCTK_InterpGridArrays (ierror, cctkGH, 3, interp_handle,
+ . param_table_handle, coord_system_handle,
+ . npoints, CCTK_VARIABLE_REAL, interp_coords,
+ . 9, in_array_indices,
+ . 9, out_array_type_codes, out_arrays)
+ if (ierror < 0) then
+ call CCTK_WARN (1, "interpolator call returned an error code");
+ endif
! For H or C flows, we need the interpolated expansion, the
! norm of the gradient of the horizon function, and the mask.
else
- num_arrays = 3
- call CCTK_InterpGV (ierror, cctkGH, interp_handle, coord_system_handle,
- . npoints, num_arrays, num_arrays,
- . xa, ya, za,
- . CCTK_VARIABLE_REAL, CCTK_VARIABLE_REAL,
- . CCTK_VARIABLE_REAL,
- . in_array_indices(7), in_array_indices(8),
- . in_array_indices(9),
- . exp, intmask, gradn,
- . CCTK_VARIABLE_REAL, CCTK_VARIABLE_REAL,
- . CCTK_VARIABLE_REAL)
+ call CCTK_InterpGridArrays (ierror, cctkGH, 3, interp_handle,
+ . param_table_handle, coord_system_handle,
+ . npoints, CCTK_VARIABLE_REAL, interp_coords,
+ . 3, in_array_indices(7),
+ . 3, out_array_type_codes(7), out_arrays(7))
+ if (ierror < 0) then
+ call CCTK_WARN (1, "interpolator call returned an error code");
+ endif
end if
+! release parameter table
+ call Util_TableDestroy (ierror, param_table_handle)
! ***************************************
! *** CHECK IF WE ARE INSIDE MASK ***
@@ -937,27 +942,27 @@
! Reduce integrals.
call CCTK_ReduceLocalScalar(ierror,cctkGH,-1,reduction_handle,
- . intarea,aux,CCTK_VARIABLE_REAL)
+ . intarea,aux,CCTK_VARIABLE_REAL)
intarea=aux
call CCTK_ReduceLocalScalar(ierror,cctkGH,-1,reduction_handle,
- . intexp,aux,CCTK_VARIABLE_REAL)
+ . intexp,aux,CCTK_VARIABLE_REAL)
intexp=aux
call CCTK_ReduceLocalScalar(ierror,cctkGH,-1,reduction_handle,
- . intexp2,aux,CCTK_VARIABLE_REAL)
+ . intexp2,aux,CCTK_VARIABLE_REAL)
intexp2=aux
call CCTK_ReduceLocalScalar(ierror,cctkGH,-1,reduction_handle,
- . intexpdel2,aux,CCTK_VARIABLE_REAL)
+ . intexpdel2,aux,CCTK_VARIABLE_REAL)
intexpdel2=aux
call CCTK_ReduceLocalScalar(ierror,cctkGH,-1,reduction_handle,
- . inside_min_neg_count,red_tmp,CCTK_VARIABLE_INT)
+ . inside_min_neg_count,red_tmp,CCTK_VARIABLE_INT)
inside_min_neg_count=red_tmp
call CCTK_ReduceLocalScalar(ierror,cctkGH,-1,reduction_handle,
- . inside_min_count,red_tmp,CCTK_VARIABLE_INT)
+ . inside_min_count,red_tmp,CCTK_VARIABLE_INT)
inside_min_count=red_tmp
end if
@@ -1033,9 +1038,9 @@
zw = za(i,j) - zc
! Calculate metric on surface.
-
- det = txx(i,j)*tyy(i,j)*tzz(i,j)
- . + two*txy(i,j)*txz(i,j)*tyz(i,j)
+
+ det = txx(i,j)*tyy(i,j)*tzz(i,j)
+ . + two*txy(i,j)*txz(i,j)*tyz(i,j)
. - txx(i,j)*tyz(i,j)**2 - tyy(i,j)*txz(i,j)**2
. - tzz(i,j)*txy(i,j)**2
@@ -1044,7 +1049,7 @@
gupij(1,1) = idet*(tyy(i,j)*tzz(i,j)-tyz(i,j)**2)
gupij(2,2) = idet*(txx(i,j)*tzz(i,j)-txz(i,j)**2)
gupij(3,3) = idet*(txx(i,j)*tyy(i,j)-txy(i,j)**2)
-
+
gupij(1,2) = idet*(txz(i,j)*tyz(i,j)-tzz(i,j)*txy(i,j))
gupij(2,1) = gupij(1,2)
gupij(1,3) = idet*(txy(i,j)*tyz(i,j)-tyy(i,j)*txz(i,j))
@@ -1122,21 +1127,21 @@
auxi = lmax+1
- if (hw.ne.zero) then
+ if (hw.ne.zero) then
call CCTK_ReduceLocalArray1D(ierror,cctkGH,-1,reduction_handle,
- . hflow0,tempv,auxi,CCTK_VARIABLE_REAL)
+ . hflow0,tempv,auxi,CCTK_VARIABLE_REAL)
hflow0 = tempv
end if
- if (cw.ne.zero) then
+ if (cw.ne.zero) then
call CCTK_ReduceLocalArray1D(ierror,cctkGH,-1,reduction_handle,
- . cflow0,tempv,auxi,CCTK_VARIABLE_REAL)
+ . cflow0,tempv,auxi,CCTK_VARIABLE_REAL)
cflow0 = tempv
end if
if (nw.ne.zero) then
call CCTK_ReduceLocalArray1D(ierror,cctkGH,-1,reduction_handle,
- . nflow0,tempv,auxi,CCTK_VARIABLE_REAL)
+ . nflow0,tempv,auxi,CCTK_VARIABLE_REAL)
nflow0 = tempv
end if
@@ -1146,41 +1151,41 @@
auxi = lmax*lmax
- if (hw.ne.zero) then
+ if (hw.ne.zero) then
call CCTK_ReduceLocalArray1D(ierror,cctkGH,-1,
- . reduction_handle,hflowc,tempm,auxi,CCTK_VARIABLE_REAL)
+ . reduction_handle,hflowc,tempm,auxi,CCTK_VARIABLE_REAL)
hflowc = tempm
- end if
+ end if
if (cw.ne.zero) then
call CCTK_ReduceLocalArray1D(ierror,cctkGH,-1,
- . reduction_handle,cflowc,tempm,auxi,CCTK_VARIABLE_REAL)
+ . reduction_handle,cflowc,tempm,auxi,CCTK_VARIABLE_REAL)
cflowc = tempm
end if
if (nw.ne.zero) then
call CCTK_ReduceLocalArray1D(ierror,cctkGH,-1,
- . reduction_handle,nflowc,tempm,auxi,CCTK_VARIABLE_REAL)
+ . reduction_handle,nflowc,tempm,auxi,CCTK_VARIABLE_REAL)
nflowc = tempm
end if
if (.not.refy) then
- if (hw.ne.zero) then
+ if (hw.ne.zero) then
call CCTK_ReduceLocalArray1D(ierror,cctkGH,-1,
- . reduction_handle,hflows,tempm,auxi,CCTK_VARIABLE_REAL)
+ . reduction_handle,hflows,tempm,auxi,CCTK_VARIABLE_REAL)
hflows = tempm
- end if
+ end if
if (cw.ne.zero) then
call CCTK_ReduceLocalArray1D(ierror,cctkGH,-1,
- . reduction_handle,cflows,tempm,auxi,CCTK_VARIABLE_REAL)
+ . reduction_handle,cflows,tempm,auxi,CCTK_VARIABLE_REAL)
cflows = tempm
end if
if (nw.ne.zero) then
call CCTK_ReduceLocalArray1D(ierror,cctkGH,-1,
- . reduction_handle,nflows,tempm,auxi,CCTK_VARIABLE_REAL)
+ . reduction_handle,nflows,tempm,auxi,CCTK_VARIABLE_REAL)
nflows = tempm
end if
@@ -1236,7 +1241,7 @@
if (flow) then
! Axisymmetric.
-
+
hflow0 = intw*hflow0
cflow0 = intw*cflow0
nflow0 = intw*nflow0
generated by cgit v1.2.3 (git 2.39.1) at 2024-06-17 18:02:02 +0000