diff options
author | jthorn <jthorn@0a4070d5-58f5-498f-b6c0-2693e757fa0f> | 2005-05-20 13:25:38 +0000 |
---|---|---|
committer | jthorn <jthorn@0a4070d5-58f5-498f-b6c0-2693e757fa0f> | 2005-05-20 13:25:38 +0000 |
commit | e20c593a8ea763e2d6c1b24ca56e9cfe4fc07e07 (patch) | |
tree | 75ca7ecfca0d0b61dedc2dce188628d9cdd8ff4d /src | |
parent | 2f9378fe1f03a73153bbc82d5efb2f5d282d15c4 (diff) |
This is a general code cleanup, and also adds some new features to
this thorn.
[Unforunately, for the reasons explained in today's E-mail discussion
on the developers@cactuscode.org mailing list, it's not practical for
me to break this up into multiple smaller single-purposed commits.]
param.ccl
* add comments
* add additional interpolator parameters (described in detail below)
* add "output psi on 2D grid" parameters (described in detail below)
* add debugging parameters
doc/documentation.tex
* add "IDAxiBrillBH/" prefix to make LaTeX labels unique across thorns
* correct a few typos
* clarify explanation of how we solve on a 2-D grid and then interpolate
to the Cactus 3-D grid
* explain the new interpolator and output-psi-on-2D-grid parameters
* add some comments on choosing the 2-D grid resolution parameters
* mention debug parameters
doc/TODO
* new file noting some problems with this thorn which I found, but
didn't fix
src/IDAxiBrillBH.F
* use new CCTK_WARN(CCTK_WARN_ABORT, ...) instead of old CCTK_WARN(0, ...)
* more flexible interpolator parameters:
This thorn first solves the Brill-wave equation on an internal
(axisymmetric) 2D grid, then uses uses the standard Cactus
CCTK_InterpLocalUniform() local-interpolator API to interpolate
this to the 3D grid points. Before this patch, this thorn
hard-wired the interpolation operator to "uniform cartesian",
only allowed orders 1, 2, and 3, and didn't allow any other
interpolator parameters to be set. This patch adds two new
parameters to allow any interpolation operator and parameters
to be specified. Either the old or the new parameters can be
used (see doc/documentation.tex for details of how this works);
the defaults leave the behavior of this thorn unchanged.
* add option to output psi on 2D grid
When using this thorn, I found it hard to choose the parameters
defining the resolution and extent of the internal 2D grid.
To help with this, I added an option to output the solution
on that grid to an ASCII data file, so it can be examined
and plotted.
* add debugging code
I've added debug options which, if enabled, print the values
of various internal quantities at specified 2D and/or 3D grid
points. I've left these in the final "production" code as
possibly being useful for future debugging.
* add many comments
* correct or remove some old comments which were out of date
* adjust whitespace to make the code a bit more readable:
I've changed code like
allocate(cc(neb,nqb),ce(neb,nqb),cw(neb,nqb),cn(neb,nqb),cs(neb,nqb),
$ rhs(neb,nqb),psi2d(neb,nqb),detapsi2d(neb,nqb),dqpsi2d(neb,nqb),
$ detaetapsi2d(neb,nqb),detaqpsi2d(neb,nqb),dqqpsi2d(neb,nqb),
$ etagrd(neb),qgrd(nqb))
to make it clearer which arrays have the same size:
allocate( cc(neb,nqb))
allocate( ce(neb,nqb))
allocate( cw(neb,nqb))
allocate( cn(neb,nqb))
allocate( cs(neb,nqb))
allocate( rhs(neb,nqb))
allocate( psi2d(neb,nqb))
allocate( detapsi2d(neb,nqb))
allocate( dqpsi2d(neb,nqb))
allocate(detaetapsi2d(neb,nqb))
allocate( detaqpsi2d(neb,nqb))
allocate( dqqpsi2d(neb,nqb))
allocate(etagrd(neb))
allocate( qgrd(nqb))
* change some Fortran write statements to CCTK_INFO
[so they're properly synchronized with output from C routines,
even when stdout+stderr are redirected to a log file]
* include more information about what went wrong, in various error messages
[eg if the interpolator returns an error code, the code
now includes that error code in the error message]
* rename a few variables to make their purpose clearer:
ep1 --> error_at_this_grid_point
ep2 --> max_error_in_grid
git-svn-id: http://svn.einsteintoolkit.org/cactus/EinsteinInitialData/IDAxiBrillBH/trunk@69 0a4070d5-58f5-498f-b6c0-2693e757fa0f
Diffstat (limited to 'src')
-rw-r--r-- | src/IDAxiBrillBH.F | 394 |
1 files changed, 291 insertions, 103 deletions
diff --git a/src/IDAxiBrillBH.F b/src/IDAxiBrillBH.F index 52dfe3d..7637776 100644 --- a/src/IDAxiBrillBH.F +++ b/src/IDAxiBrillBH.F @@ -45,8 +45,7 @@ c@@*/ real*8, allocatable :: psi2dv(:,:,:),detapsi2dv(:,:,:), $ dqpsi2dv(:,:,:),detaetapsi2dv(:,:,:),detaqpsi2dv(:,:,:), $ dqqpsi2dv(:,:,:) - parameter(axibheps = 1.0d-12) - real*8 ep1,ep2 + real*8 error_at_this_grid_point,max_error_in_grid real*8 o1,o2,o3,o4,o5,o6,o7,o8,o9 real*8 o10,o11,o12,o13,o14,o15,o16,o17,o18,o19 real*8 o20,o21,o22,o23,o24,o25,o26,o27,o28,o29 @@ -60,13 +59,24 @@ c@@*/ integer i22 real*8 pi real*8 adm + real*8 exp_mhalf_eta, psi3d integer :: nx,ny,nz integer i,j,k,nquads integer npoints,ierror integer neb, nqb + integer posn + integer io_status + + integer, parameter :: max_string_length = 500 + + integer param_table_handle, interp_handle + character*max_string_length :: message_buffer + + integer fstring_length + character*max_string_length :: interpolator_name_fstring + character*max_string_length :: interpolator_pars_fstring + character*max_string_length :: output_psi2D_file_name_fstring - 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 @@ -81,31 +91,104 @@ c Set up the grid spacings 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... +c +c ***** set up interpolator handle and parameter table ***** +c + +c +c ... we do this now, rather than waiting till we need them, +c so any errors (eg user forgot to activate a suitable interpolator thorn) +c get printed right away, rather than making the user wait for them +c ... n.b. we must first convert our C-style string parameters +c to Fortran strings +c + call CCTK_FortranString(fstring_length, + $ interpolator_name, + $ interpolator_name_fstring) + call CCTK_InterpHandle(interp_handle, + $ interpolator_name_fstring(1:fstring_length)) + if (interp_handle .lt. 0) then + call CCTK_WARN(CCTK_WARN_ABORT, "Cannot get interpolator handle! Did you forgot to activate a suitable local-interpolator thorn?") + endif + +c +c build the interpolator parameter table from a suitable string: +c if interpolator_pars is a nonempty string, use that, otherwise +c use "order=n" where n is given by the interpolation_order +c parameter +c + call CCTK_FortranString(fstring_length, + $ interpolator_pars, + $ interpolator_pars_fstring) + if (fstring_length .eq. 0) then + interpolator_pars_fstring + $ = "order=" // char(ichar('0') + interpolation_order) + endif + call Util_TableCreateFromString + $ (param_table_handle, + $ interpolator_pars_fstring(1:fstring_length)) + if (param_table_handle .lt. 0) then + write(message_buffer, '(A,I)') + $ 'failed to create interpolator param table: error code ', + $ param_table_handle + call CCTK_WARN(CCTK_WARN_ABORT, message_buffer) + endif + +c +c ***** solve Brill-wave equation on 2D (eta,theta) grid ***** +c +c The code uses the abbreviation 'q' for theta. This is a bit confusing, +c because this is *not* the same quantity as $q(\eta,\theta)$ described +c in the thorn guide (which is is stored in the 3-D array q(nx,ny,nz)). +c +c The (eta,theta) grid spans the range +c eta in [0,etamax] + ghost zones (neb points, spacing deta) +c theta in [0,pi ] + ghost zones (nqb points, spacing dq ) +c +c 2D grid size NE x NQ, plus 2 zones for boundaries +c c 21/11/00 TR: dont change parameters in place c but keep a copy in local variables c Otherwise the changed parameters cause trouble c after recovery. +c neb = ne+2 nqb = nq+2 - ! do I need to call free? - allocate(cc(neb,nqb),ce(neb,nqb),cw(neb,nqb),cn(neb,nqb),cs(neb,nqb), - $ rhs(neb,nqb),psi2d(neb,nqb),detapsi2d(neb,nqb),dqpsi2d(neb,nqb), - $ detaetapsi2d(neb,nqb),detaqpsi2d(neb,nqb),dqqpsi2d(neb,nqb), - $ etagrd(neb),qgrd(nqb)) - allocate(eta(nx,ny,nz),abseta(nx,ny,nz),sign_eta(nx,ny,nz), - $ q(nx,ny,nz),phi(nx,ny,nz), - $ psi2dv(nx,ny,nz),detapsi2dv(nx,ny,nz),dqpsi2dv(nx,ny,nz), - $ detaetapsi2dv(nx,ny,nz),detaqpsi2dv(nx,ny,nz), - $ dqqpsi2dv(nx,ny,nz)) + + allocate( cc(neb,nqb)) + allocate( ce(neb,nqb)) + allocate( cw(neb,nqb)) + allocate( cn(neb,nqb)) + allocate( cs(neb,nqb)) + allocate( rhs(neb,nqb)) + allocate( psi2d(neb,nqb)) + allocate( detapsi2d(neb,nqb)) + allocate( dqpsi2d(neb,nqb)) + allocate(detaetapsi2d(neb,nqb)) + allocate( detaqpsi2d(neb,nqb)) + allocate( dqqpsi2d(neb,nqb)) + + allocate(etagrd(neb)) + allocate( qgrd(nqb)) + + allocate( eta(nx,ny,nz)) + allocate( abseta(nx,ny,nz)) + allocate( sign_eta(nx,ny,nz)) + allocate( q(nx,ny,nz)) + allocate( phi(nx,ny,nz)) + allocate( psi2dv(nx,ny,nz)) + allocate( detapsi2dv(nx,ny,nz)) + allocate( dqpsi2dv(nx,ny,nz)) + allocate(detaetapsi2dv(nx,ny,nz)) + allocate( detaqpsi2dv(nx,ny,nz)) + allocate( dqqpsi2dv(nx,ny,nz)) + c Initialize some arrays psi2d = 1.0d0 detapsi2d = 0.0d0 nquads = 2 - dq = nquads*0.5*pi/(nqb-2) + dq = nquads*0.5d0*pi/(nqb-2) deta = etamax/(neb-3) do j=1,nqb @@ -115,7 +198,7 @@ c Initialize some arrays #include "CactusEinstein/IDAxiBrillBH/src/bhbrill.x" enddo enddo -c Boundary conditions +c Boundary conditions do j=1,nqb ce(2,j)=ce(2,j)+cw(2,j) cw(2,j)=0.0 @@ -132,7 +215,8 @@ c Boundary conditions cn(i,nqb-1)=0.0 enddo -c Do the solve +c Do the solve + axibheps = error_tolerance call CCTK_INFO("Calling axisymmetric solver") call mgparm (levels,5,id5,id9,idi,idg,neb,nqb) @@ -142,45 +226,55 @@ c Do the solve call CCTK_INFO("Solve complete") -c The solution is now available. -c Debugging is needed, a stop statement should -c be called if the IVP solve is not successful - +c +c The solution is (hopefully) now available. +c if(ier .ne. 0) then - call CCTK_WARN(0,"Solution to BH+Brill Wave not found") + write(message_buffer, '(A,I)') + $ 'failed to solve elliptic equation: ier=', ier + call CCTK_WARN(CCTK_WARN_ABORT, message_buffer) end if print *,'rmax = ',rmax print *,'axibheps = ',axibheps print *,'psi2d = ',maxval(psi2d),' ',minval(psi2d) - ep2 = 0.0 + max_error_in_grid = 0.0d0 do j=2,nqb-1 do i=2,neb-1 - ep1 = rhs(i,j)-psi2d(i,j)*cc(i,j)-psi2d(i,j+1)*cn(i,j)-psi2d(i,j-1)*cs(i,j)- - & psi2d(i+1,j)*ce(i,j)-psi2d(i-1,j)*cw(i,j) - ep2 = max(abs(ep1),ep2) + error_at_this_grid_point = rhs(i,j) + $ - psi2d(i,j )*cc(i,j) + $ - psi2d(i,j+1)*cn(i,j) + $ - psi2d(i,j-1)*cs(i,j) + $ - psi2d(i+1,j)*ce(i,j) + $ - psi2d(i-1,j)*cw(i,j) + max_error_in_grid = max(max_error_in_grid, + $ abs(error_at_this_grid_point)) enddo enddo - print *,'Resulting eps =',ep2 + print *,'Resulting eps =',max_error_in_grid +c Boundary conditions do j=1,nqb - psi2d(1,j)=psi2d(3,j) - psi2d(neb,j)=-deta*psi2d(neb-1,j)+psi2d(neb-2,j) + psi2d(1 ,j) = psi2d(3,j) + psi2d(neb,j) = - deta*psi2d(neb-1,j) + psi2d(neb-2,j) enddo do i=1,neb - psi2d(i,1)=psi2d(i,2) - psi2d(i,nqb)=psi2d(i,nqb-1) + psi2d(i,1 ) = psi2d(i,2) + psi2d(i,nqb) = psi2d(i,nqb-1) enddo +c derivatives of psi do j=2,nqb-1 do i=2,neb-1 - dqpsi2d(i,j)=0.5*(psi2d(i,j+1)-psi2d(i,j-1))/dq - dqqpsi2d(i,j)=(psi2d(i,j+1)+psi2d(i,j-1)-2.*psi2d(i,j))/dq**2 - detapsi2d(i,j)=sinh(0.5*etagrd(i))+0.5*(psi2d(i+1,j)-psi2d(i-1,j))/deta - detaetapsi2d(i,j)=0.5*cosh(0.5*etagrd(i))+ - $ (psi2d(i+1,j)+psi2d(i-1,j)-2.*psi2d(i,j))/deta**2 + dqpsi2d (i,j) = 0.5d0*(psi2d(i,j+1)-psi2d(i,j-1))/dq + dqqpsi2d (i,j) = (psi2d(i,j+1)+psi2d(i,j-1)-2.0d0*psi2d(i,j))/dq**2 + detapsi2d(i,j) = sinh(0.5d0*etagrd(i)) + $ + 0.5d0*(psi2d(i+1,j)-psi2d(i-1,j))/deta + detaetapsi2d(i,j) + $ = 0.5d0*cosh(0.5d0*etagrd(i)) + $ + (psi2d(i+1,j)+psi2d(i-1,j)-2.0d0*psi2d(i,j))/deta**2 enddo enddo @@ -214,7 +308,7 @@ c be called if the IVP solve is not successful do j=2,nqb-1 do i=2,neb-1 - detaqpsi2d(i,j)=0.5*(detapsi2d(i,j+1)-detapsi2d(i,j-1))/dq + detaqpsi2d(i,j)=0.5d0*(detapsi2d(i,j+1)-detapsi2d(i,j-1))/dq enddo enddo @@ -229,13 +323,84 @@ c be called if the IVP solve is not successful enddo do j=1,nqb - psi2d(:,j)=psi2d(:,j)+2.*cosh(0.5*etagrd) + psi2d(:,j)=psi2d(:,j)+2.0d0*cosh(0.5d0*etagrd) enddo -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 + if (debug .ge. 6) then + print *, '### 2-D grid results (= inputs to interpolation) ###' + print *, 'effective 2-D grid size: neb,nqb =', neb,nqb + print *, 'debug_{ii,jj} =', debug_ii,debug_jj + print *, 'at this 2-D grid point...' + print *, ' eta =', etagrd(debug_ii) + print *, ' theta [q] =', qgrd(debug_jj) + print *, ' psi2d =', psi2d(debug_ii,debug_jj) + print *, ' detapsi2d =', detapsi2d(debug_ii,debug_jj) + print *, ' dqpsi2d =', dqpsi2d(debug_ii,debug_jj) + print *, ' detaetapsi2d =', detaetapsi2d(debug_ii,debug_jj) + print *, ' detaqpsi2d =', detaqpsi2d(debug_ii,debug_jj) + print *, ' dqqpsi2d =', dqqpsi2d(debug_ii,debug_jj) + endif + +c +c write conformal factor psi on 2D grid to output file if requested +c + if (output_psi2D .ne. 0) then + write (message_buffer, '(A,A,A)') + $ 'writing 2D psi to "', + $ output_psi2D_file_name_fstring(1:fstring_length), + $ '"' + call CCTK_INFO(message_buffer) + + call CCTK_FortranString(fstring_length, + $ output_psi2D_file_name, + $ output_psi2D_file_name_fstring) + open (9, iostat=io_status, status='replace', + $ file=output_psi2D_file_name_fstring) + if (io_status .ne. 0) then + write (message_buffer, '(A,A,A,I)') + $ 'error opening psi2D output file "', + $ output_psi2D_file_name_fstring(1:fstring_length), + $ '": io_status=', io_status + call CCTK_WARN(CCTK_WARN_ABORT, message_buffer) + endif + + write (9, '(a)') + $ '# eta theta (=q) psi (2D) psi (3D)' + do i = 1,neb + exp_mhalf_eta = exp(-0.5d0 * etagrd(i)) + do j = 1,nqb + psi3d = psi2d(i,j) * exp_mhalf_eta + write (9, '(f12.8,a, f12.8,a, g20.10e3,a, g20.10e3)') + $ etagrd(i), ' ', qgrd(j), ' ', + $ psi2d(i,j), ' ', psi3d + end do + write (9, '(a)') '' + end do + + close (9, iostat=io_status) + if (io_status .ne. 0) then + write(message_buffer, '(A,I)') + $ 'error closing psi2D output file: io_status=', io_status + call CCTK_WARN(CCTK_WARN_ABORT, message_buffer) + endif + endif -c Conformal factor + + +c +c ***** interpolate psi and its derivatives to the xyz grid positions ***** +c ***** and compute the ADM variables there ***** +c +c More precisely, at this point we have q, psi, and psi's derivatives +c on the (eta,theta) grid. We want to interpolate these values to the +c (eta,theta) locations of each of the (x,y,z) grid points. +c +c The (eta,theta) grid only spans the range eta >= 0, so we actually +c interpolate using the (x,y,z) grid points' |eta| values, and fix +c things up afterwords. +c + + call CCTK_INFO("interpolating solution to xyz grid points") eta = 0.5d0 * dlog (x**2 + y**2 + z**2) abseta = abs (eta) @@ -245,60 +410,66 @@ c Conformal factor do k=1,nz do j=1,ny do i=1,nx -c eta(i,j,k) = 0.5d0*dlog(x(i,j,k)**2+y(i,j,k)**2+z(i,j,k)**2) -c abseta(i,j,k) = abs(eta(i,j,k)) if(eta(i,j,k) .lt. 0)then sign_eta(i,j,k) = -1 else sign_eta(i,j,k) = 1 endif -c q(i,j,k) = datan2(sqrt(x(i,j,k)**2+y(i,j,k)**2),z(i,j,k)) -c TYPO HERE ??????????? -c | -c | -c phi(i,j,k)= datan2(y(i,j,k),x(i,j,k)) enddo enddo enddo - npoints = nx*ny*nz - -! 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") + if (debug .ge. 6) then +c posn = (0-origin) 1-D positionn of this point in interpolator arrays +c (this is useful because interpolator error messages refer to it) + posn = (debug_i-1) + nx*(debug_j-1) + nx*ny*(debug_k-1) + print *, '### 3-D interpolation coordinates and inputs ###' + print *, '3-D grid size: n[xyz] =', nx,ny,nz + print *, 'debug_[ijk] =', debug_i,debug_j,debug_k, '==> posn =', posn + print *, 'at this 3-D grid point, ...' + print *, ' x =', x(debug_i,debug_j,debug_k) + print *, ' y =', y(debug_i,debug_j,debug_k) + print *, ' z =', z(debug_i,debug_j,debug_k) + print *, ' eta =', eta(debug_i,debug_j,debug_k) + print *, ' abseta =', abseta(debug_i,debug_j,debug_k) + print *, 'theta [q] =', q(debug_i,debug_j,debug_k) + print *, ' phi =', phi(debug_i,debug_j,debug_k) endif - 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 +c set up the interpolator array pointers + npoints = nx*ny*nz -! fill in the input/output arrays for the interpolator coord_origin(1) = etagrd(1) - coord_origin(2) = qgrd(1) + coord_origin(2) = qgrd(1) coord_delta(1) = etagrd(2) - etagrd(1) - coord_delta(2) = qgrd(2) - qgrd(1) + coord_delta(2) = qgrd(2) - qgrd(1) + + if (debug .ge. 6) then + print *, '### 2-D grid origin: eta=', coord_origin(1) + print *, ' theta=', coord_origin(2) + print *, ' delta: eta=', coord_delta(1) + print *, ' theta=', coord_delta(2) + end if interp_coords(1) = CCTK_PointerTo(abseta) interp_coords(2) = CCTK_PointerTo(q) - in_array_dims(1) = neb; in_array_dims(2) = nqb + 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(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) + 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(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) + out_arrays(5) = CCTK_PointerTo( detaqpsi2dv) + out_arrays(6) = CCTK_PointerTo( dqqpsi2dv) call CCTK_InterpLocalUniform (ierror, 2, $ interp_handle, param_table_handle, @@ -306,50 +477,51 @@ c phi(i,j,k)= datan2(y(i,j,k),x(i,j,k)) $ npoints, type_codes(1), interp_coords, $ 6, in_array_dims, type_codes, in_arrays, $ 6, type_codes, out_arrays) + if (ierror < 0) then + write(message_buffer, '(A,I)') + $ 'error in interpolator: ierror=', ierror + call CCTK_WARN(CCTK_WARN_ABORT, message_buffer) + endif call Util_TableDestroy (ierror, param_table_handle) - psi = psi2dv * exp (-0.5 * eta) + if (debug .ge. 6) then + print *, '### interpolation results (at this 3-D grid point) ###' + print *, ' psi2dv =', psi2dv(debug_i,debug_j,debug_k) + end if +c +c ***** compute the ADMBase conformal factor, its derivatives, ***** +c ***** metric, and extrinsic curvature from the interpolation output ***** +c + psi = psi2dv * exp (-0.5d0 * eta) detapsi2dv = sign_eta * detapsi2dv detaqpsi2dv = sign_eta * detaqpsi2dv 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 - do k=1,nz 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 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)) +c Conformal metric +c gxx = ... +c Derivatives of the metric (currently commented-out) +c dxgxx = 1/2 \partial gxx / \partial x #include "CactusEinstein/IDAxiBrillBH/src/gij.x" enddo enddo enddo -c Curvature +c Extrinsic Curvature is identically zero kxx = 0.0D0 kxy = 0.0D0 kxz = 0.0D0 @@ -357,25 +529,41 @@ c Curvature kyz = 0.0D0 kzz = 0.0D0 - 111 continue -c Set ADM mass + if (debug .ge. 6) then + print *, '### final results (again at this 3-D grid point) ###' + print *, ' psi =', psi(debug_i,debug_j,debug_k) + print *, ' gxx =', gxx(debug_i,debug_j,debug_k) + print *, ' gxy =', gxy(debug_i,debug_j,debug_k) + print *, ' gxz =', gxz(debug_i,debug_j,debug_k) + print *, ' gyy =', gyy(debug_i,debug_j,debug_k) + print *, ' gyz =', gyz(debug_i,debug_j,debug_k) + print *, ' gzz =', gzz(debug_i,debug_j,debug_k) + endif + +c +c ***** diagnostics and final cleanup +c + +c ADM mass + call CCTK_INFO("computing ADM mass") i = neb-15 - adm = 0.0 + adm = 0.0d0 do j=2,nqb-1 - adm=adm+(psi2d(i,j)-(psi2d(i+1,j)-psi2d(i-1,j))/deta)*exp(0.5* - $ etagrd(i)) + adm = adm + $ + (psi2d(i,j)-(psi2d(i+1,j)-psi2d(i-1,j))/deta) + $ *exp(0.5d0*etagrd(i)) enddo adm=adm/(nqb-2) print *,'ADM mass: ',adm + 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) + write (*,*)"using alp = (2r - adm)/(2r+adm)" + alp = (2.0d0*r - adm)/(2.0d0*r+adm) endif -c conformal_state = CONFORMAL_METRIC (What is CONFORMAL_METRIC?) +c conformal_state = 'all' derivatives were calculated conformal_state = 3 -c 3 ==> 'all' derivatives were calculated deallocate(cc,ce,cw,cn,cs,rhs,psi2d,detapsi2d,dqpsi2d, $ detaetapsi2d,detaqpsi2d,dqqpsi2d, |