From 9c1798b974d3ed1f9c48b5d29457e6edf01345c4 Mon Sep 17 00:00:00 2001 From: tradke Date: Mon, 17 May 2004 09:53:11 +0000 Subject: Use CCTK_InterpLocalUniform() instead of the deprecated CCTK_InterpLocal() API. git-svn-id: http://svn.einsteintoolkit.org/cactus/EinsteinInitialData/IDAxiBrillBH/trunk@62 0a4070d5-58f5-498f-b6c0-2693e757fa0f --- src/IDAxiBrillBH.F | 174 +++++++++++++++++++++++++++++------------------------ 1 file changed, 97 insertions(+), 77 deletions(-) diff --git a/src/IDAxiBrillBH.F b/src/IDAxiBrillBH.F index 58b4477..a80b16e 100644 --- a/src/IDAxiBrillBH.F +++ b/src/IDAxiBrillBH.F @@ -1,40 +1,39 @@ c/*@@ c @file IDAxiBrillBH.F -c @date -c @author -c @desc -c -c @enddesc +c @date +c @author +c @desc +c +c @enddesc c@@*/ -#include "cctk.h" +#include "cctk.h" #include "cctk_Parameters.h" #include "cctk_Arguments.h" +#include "cctk_Functions.h" c/*@@ c @routine IDAxiBrillBH -c @date -c @author -c @desc -c -c @enddesc -c @calls -c @calledby -c @history +c @date +c @author +c @desc +c +c @enddesc +c @calls +c @calledby +c @history c -c @endhistory +c @endhistory c@@*/ - + subroutine IDAxiBrillBH(CCTK_ARGUMENTS) - + implicit none - + DECLARE_CCTK_ARGUMENTS DECLARE_CCTK_PARAMETERS + DECLARE_CCTK_FUNCTIONS -c Perhaps this and others should go into cctk.h - integer CCTK_Equals - real*8 axibheps, rmax, dq, deta integer levels,id5,id9,idi,idg,ier real*8, allocatable :: cc(:,:),ce(:,:),cw(:,:),cn(:,:),cs(:,:), @@ -63,16 +62,25 @@ c Perhaps this and others should go into cctk.h real*8 adm integer :: nx,ny,nz integer i,j,k,nquads - integer npoints,handle,ierror + integer npoints,ierror integer neb, nqb - + + integer param_table_handle, interp_handle + character(30) options_string + CCTK_REAL, dimension(2) :: coord_origin, coord_delta + CCTK_POINTER, dimension(2) :: interp_coords + CCTK_POINTER, dimension(6) :: in_arrays, out_arrays + CCTK_INT, dimension(2) :: in_array_dims + CCTK_INT, dimension(6), parameter :: type_codes = CCTK_VARIABLE_REAL + + pi = 4.0d0*atan(1.0d0) - + c Set up the grid spacings nx = cctk_lsh(1) ny = cctk_lsh(2) nz = cctk_lsh(3) - + c Solve on this sized cartesian grid c 2D grid size NE x NQ c Add 2 zones for boundaries... @@ -101,27 +109,27 @@ c Initialize some arrays deta = etamax/(neb-3) do j=1,nqb - qgrd(j) = (j-1.5)*dq - do i=1,neb - etagrd(i) = (i-2)*deta + qgrd(j) = (j-1.5)*dq + do i=1,neb + etagrd(i) = (i-2)*deta #include "CactusEinstein/IDAxiBrillBH/src/bhbrill.x" - enddo + enddo enddo c Boundary conditions do j=1,nqb - ce(2,j)=ce(2,j)+cw(2,j) - cw(2,j)=0.0 + ce(2,j)=ce(2,j)+cw(2,j) + cw(2,j)=0.0 - cw(neb-1,j)=cw(neb-1,j)+ce(neb-1,j) - cc(neb-1,j)=cc(neb-1,j)-deta*ce(neb-1,j) - ce(neb-1,j)=0.0 + cw(neb-1,j)=cw(neb-1,j)+ce(neb-1,j) + cc(neb-1,j)=cc(neb-1,j)-deta*ce(neb-1,j) + ce(neb-1,j)=0.0 enddo do i=1,neb - cc(i,2)=cc(i,2)+cs(i,2) - cs(i,2)=0.0 - cc(i,nqb-1)=cc(i,nqb-1)+cn(i,nqb-1) - cn(i,nqb-1)=0.0 + cc(i,2)=cc(i,2)+cs(i,2) + cs(i,2)=0.0 + cc(i,nqb-1)=cc(i,nqb-1)+cn(i,nqb-1) + cn(i,nqb-1)=0.0 enddo c Do the solve @@ -224,7 +232,7 @@ c be called if the IVP solve is not successful psi2d(:,j)=psi2d(:,j)+2.*cosh(0.5*etagrd) enddo -c Now evaluate each of the following at x(i,j,k), y(i,j,k) and +c Now evaluate each of the following at x(i,j,k), y(i,j,k) and c z(i,j,k) where i,j,k go from 1 to nx,ny,nz c Conformal factor @@ -255,37 +263,49 @@ c phi(i,j,k)= datan2(y(i,j,k),x(i,j,k)) npoints = nx*ny*nz -! Interpolator handle. - - handle = -1 - - if (interpolation_order .eq. 1) then - call CCTK_InterpHandle (handle, "first-order uniform cartesian") - else if (interpolation_order .eq. 2) then - call CCTK_InterpHandle (handle, "second-order uniform cartesian") - else if (interpolation_order .eq. 3) then - call CCTK_InterpHandle (handle, "third-order uniform cartesian") +! Parameter table and interpolator handles. + 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 - if (handle .lt. 0) then - call CCTK_WARN (0, "Couldn't get handle for interpolation operator") + call CCTK_InterpHandle (interp_handle, "uniform cartesian") + if (interp_handle .lt. 0) then + call CCTK_WARN(0,"Cannot get handle for interpolation ! Forgot to activate an implementation providing interpolation operators ??") endif - call CCTK_InterpLocal (ierror, cctkGH, handle, npoints, 2, 6, 6, - $ neb, nqb, etagrd, qgrd, - $ CCTK_VARIABLE_REAL, CCTK_VARIABLE_REAL, - $ abseta, q, - $ CCTK_VARIABLE_REAL,CCTK_VARIABLE_REAL, - $ psi2d, detapsi2d, dqpsi2d, detaetapsi2d, - $ detaqpsi2d, dqqpsi2d, - $ CCTK_VARIABLE_REAL,CCTK_VARIABLE_REAL, - $ CCTK_VARIABLE_REAL,CCTK_VARIABLE_REAL, - $ CCTK_VARIABLE_REAL,CCTK_VARIABLE_REAL, - $ psi2dv, detapsi2dv, dqpsi2dv, detaetapsi2dv, - $ detaqpsi2dv, dqqpsi2dv, - $ CCTK_VARIABLE_REAL,CCTK_VARIABLE_REAL, - $ CCTK_VARIABLE_REAL,CCTK_VARIABLE_REAL, - $ CCTK_VARIABLE_REAL,CCTK_VARIABLE_REAL) +! fill in the input/output arrays for the interpolator + coord_origin(1) = etagrd(1) + coord_origin(2) = qgrd(1) + coord_delta(1) = etagrd(2) - etagrd(1) + coord_delta(2) = qgrd(2) - qgrd(1) + + interp_coords(1) = CCTK_PointerTo(abseta) + interp_coords(2) = CCTK_PointerTo(q) + + in_array_dims(1) = neb; in_array_dims(2) = nqb + + in_arrays(1) = CCTK_PointerTo(psi2d) + in_arrays(2) = CCTK_PointerTo(detapsi2d) + in_arrays(3) = CCTK_PointerTo(dqpsi2d) + in_arrays(4) = CCTK_PointerTo(detaetapsi2d) + in_arrays(5) = CCTK_PointerTo(detaqpsi2d) + in_arrays(6) = CCTK_PointerTo(dqqpsi2d) + + out_arrays(1) = CCTK_PointerTo(psi2dv) + out_arrays(2) = CCTK_PointerTo(detapsi2dv) + out_arrays(3) = CCTK_PointerTo(dqpsi2dv) + out_arrays(4) = CCTK_PointerTo(detaetapsi2dv) + out_arrays(5) = CCTK_PointerTo(detaqpsi2dv) + out_arrays(6) = CCTK_PointerTo(dqqpsi2dv) + + call CCTK_InterpLocalUniform (ierror, 2, + $ interp_handle, param_table_handle, + $ coord_origin, coord_delta, + $ npoints, type_codes(1), interp_coords, + $ 6, in_array_dims, type_codes, in_arrays, + $ 6, type_codes, out_arrays) psi = psi2dv * exp (-0.5 * eta) detapsi2dv = sign_eta * detapsi2dv @@ -294,24 +314,24 @@ c phi(i,j,k)= datan2(y(i,j,k),x(i,j,k)) do k=1,nz do j=1,ny do i=1,nx - + c psi(i,j,k) = psi2dv(i,j,k)*exp(-0.5*eta(i,j,k)) c detapsi2dv(i,j,k) = sign_eta(i,j,k)*detapsi2dv(i,j,k) - + c psix = \partial psi / \partial x / psi #include "CactusEinstein/IDAxiBrillBH/src/psi_1st_deriv.x" - + c detaqpsi2dv(i,j,k) = sign_eta(i,j,k)*detaqpsi2dv(i,j,k) - + c psixx = \partial^2\psi / \partial x^2 / psi #include "CactusEinstein/IDAxiBrillBH/src/psi_2nd_deriv.x" enddo enddo enddo - + c Conformal metric c gxx = ... - + c Derivatives of the metric c dxgxx = 1/2 \partial gxx / \partial x @@ -319,7 +339,7 @@ c dxgxx = 1/2 \partial gxx / \partial x do j=1,ny do i=1,nx c THESE WERE ALREADY CALCULATED ABOVE !!! -c eta(i,j,k) = 0.5d0*dlog(x(i,j,k)**2+y(i,j,k)**2+z(i,j,k)**2) +c eta(i,j,k) = 0.5d0*dlog(x(i,j,k)**2+y(i,j,k)**2+z(i,j,k)**2) c q(i,j,k) = datan2(sqrt(x(i,j,k)**2+y(i,j,k)**2),z(i,j,k)) c phi(i,j,k) = datan2(y(i,j,k),x(i,j,k)) #include "CactusEinstein/IDAxiBrillBH/src/gij.x" @@ -334,7 +354,7 @@ c Curvature kyy = 0.0D0 kyz = 0.0D0 kzz = 0.0D0 - + 111 continue c Set ADM mass i = neb-15 @@ -345,16 +365,16 @@ c Set ADM mass enddo adm=adm/(nqb-2) print *,'ADM mass: ',adm - if (CCTK_Equals(initial_lapse,"schwarz")==1) then + if (CCTK_EQUALS(initial_lapse,"schwarz")) then write (*,*)"Initial with schwarzschild-like lapse" write (*,*)"using alp = (2.*r - adm)/(2.*r+adm)." alp = (2.*r - adm)/(2.*r+adm) endif - + c conformal_state = CONFORMAL_METRIC (What is CONFORMAL_METRIC?) conformal_state = 3 c 3 ==> 'all' derivatives were calculated - + deallocate(cc,ce,cw,cn,cs,rhs,psi2d,detapsi2d,dqpsi2d, $ detaetapsi2d,detaqpsi2d,dqqpsi2d, $ etagrd,qgrd, -- cgit v1.2.3