[Developers] patch: CactusEinstein/IDAxiBrillBH cleanup
Jonathan Thornburg
jthorn at aei.mpg.de
Fri May 20 05:30:55 CDT 2005
[[I originall sent this patch to developers@ on 9.May.2005, but
the mailing list seems to have somehow lost it -- noone else seems
to have gotten it, and it's not in the web mailing list archives.
Hence this resend...]]
Hi,
Over the past month or so I've been working on a general cleanup of
CactusEinstein/IDAxiBrillBH, and this is now in a state where I'd
like to commit the changes described by the enclosed patch, to CVS.
The changes affect 5 files in this thorn:
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
src/interp2.F --> archive/interp2.F
* this file is unused; move it from src/ to new directory archive/
[I chose to keep it in archive/, rather than only in the
CVS attic, because the CVS attic isn't usable with normal
Unix tools ("cat, awk, and grep" et al).]
Finally, a metacomment about patch granularity:
As you can see, this patch contains a whole set of changes. In an
ideal world, each change would be a separate patch. In practice,
though, splitting things up this way is quite awkward -- each patch
has its own patch-review-by-Tom latency, and until patch N is committed
to CVS, it's cumbersome to generate the diffs for patch N+1. Committing
my changes immediately to a private CVS branch would solve that problem,
but as I understand it, even on a private CVS branch, commits would
still need Tom's approval.
Apart from "switch to darcs", is there an elegant solution to this
problem?
ciao,
--
-- Jonathan Thornburg <jthorn at aei.mpg.de>
Max-Planck-Institut fuer Gravitationsphysik (Albert-Einstein-Institut),
Golm, Germany, "Old Europe" http://www.aei.mpg.de/~jthorn/home.html
"Washing one's hands of the conflict between the powerful and the
powerless means to side with the powerful, not to be neutral."
-- quote by Freire / poster by Oxfam
-------------- next part --------------
Index: param.ccl
===================================================================
RCS file: /cactusdevcvs/CactusEinstein/IDAxiBrillBH/param.ccl,v
retrieving revision 1.4
diff -u -r1.4 param.ccl
--- param.ccl 3 May 2002 17:50:31 -0000 1.4
+++ param.ccl 20 May 2005 10:19:28 -0000
@@ -18,6 +18,78 @@
private:
+############################################################
+
+#
+# ***** debugging parameters *****
+#
+
+#
+# If this parameter is set to true, we output an ASCII data file
+# giving psi on the 2D grid. The file format should be directly usable
+# by a gnuplot 'splot' command.
+#
+Boolean output_psi2D \
+ "should we output the conformal factor psi on the 2D grid?"
+{
+} false
+
+string output_psi2D_file_name \
+ "if we output the conformal factor psi on the 2D grid, \
+ what file name should we use for the output file?"
+{
+".+" :: "any non-empty string that's a valid file name"
+} "psi2D.dat"
+
+#
+# this parameter controls the amount of (potentially very detailed)
+# debugging information this thorn prints
+#
+CCTK_INT debug \
+ "level of debugging information to print \
+ (0 = none, 2 = a little, 6 = a lot, 10 = huge amounts)"
+{
+0:* :: "any integer >= 0"
+} 0
+
+#
+# to keep the output size quasi-finite, some debug printing which is
+# logicially "per grid point" on the 2-D (eta,q) grid, is actually only
+# done for this single 2-D grid point
+#
+CCTK_INT debug_ii "i coordinate for per-2D-grid-point debug printing"
+{
+* :: "any integer"
+} 14
+CCTK_INT debug_jj "j coordinate for per-2D-grid-point debug printing"
+{
+* :: "any integer"
+} 15
+
+#
+# to keep the output size quasi-finite, some debug printing which is
+# logicially "per Cactus grid point" is actually only done for this single
+# Cactus grid point
+#
+CCTK_INT debug_i "i coordinate for per-grid-point debug printing"
+{
+* :: "any integer"
+} 14
+CCTK_INT debug_j "j coordinate for per-grid-point debug printing"
+{
+* :: "any integer"
+} 15
+CCTK_INT debug_k "k coordinate for per-grid-point debug printing"
+{
+* :: "any integer"
+} 10
+
+############################################################
+
+#
+# ***** parameters controlling the Brill wave itself *****
+#
+
REAL amp "Brill wave amplitude"
{
*:* :: "No restriction"
@@ -33,28 +105,65 @@
*:* :: "No restriction"
} 1.0
-REAL etamax "Eta value for outer edge of grid"
+REAL etamax "eta value for outer edge of grid"
{
*:* :: "No restriction"
} 5.0
-
-INT n "sin^n theta in brill wave"
+INT n "sin^n theta in Brill wave"
{
*:* :: "No restriction"
} 2
-INT ne "Eta resolution for solve"
+############################################################
+
+#
+# ***** parameters for the numerical solution of the *****
+# ***** Brill-wave equation on the 2-D (eta,q) grid *****
+#
+
+INT ne "eta resolution for solve"
{
*:* :: "No restriction"
} 300
-INT nq "Theta resolution for solve"
+INT nq "theta resolution for solve"
{
*:* :: "No restriction"
} 50
+REAL error_tolerance "tolerance parameter for elliptic solver"
+{
+(0:* :: "any positive real number"
+} 1.0e-12
+
+############################################################
+
+#
+# ***** interpolation parameters *****
+#
+
+#
+# This thorn first computes the Brill wave solution on an internal 2-D
+# grid, then interpolates this to the 3-D Cactus grid. The following
+# parameters control this interpolation.
+#
+
+STRING interpolator_name \
+ "name of CCTK_InterpLocalUniform() interpolation operator"
+{
+".*" :: "any string"
+} "uniform cartesian"
+
+STRING interpolator_pars "parameters for the interpolation operator"
+{
+".*" :: \
+ "any nonempty string acceptable to Util_TableSetFromString() \
+ and to the interpolator, or the empty string to use 'order=n', \
+ where n is specified by the interpolation_order parameter"
+} ""
+
INT interpolation_order "Order for interpolation" STEERABLE = ALWAYS
{
- 1:3 :: "Choose between first, second, and third-order"
+0:9 :: "any integer accepted by the interpolator"
} 1
-------------- next part --------------
Index: src/IDAxiBrillBH.F
===================================================================
RCS file: /cactusdevcvs/CactusEinstein/IDAxiBrillBH/src/IDAxiBrillBH.F,v
retrieving revision 1.16
diff -u -r1.16 IDAxiBrillBH.F
--- src/IDAxiBrillBH.F 5 Nov 2004 11:17:23 -0000 1.16
+++ src/IDAxiBrillBH.F 20 May 2005 10:19:47 -0000
@@ -45,8 +45,7 @@
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 @@
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 @@
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 @@
#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 @@
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 @@
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 @@
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 @@
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 @@
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 @@
$ 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 @@
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,
-------------- next part --------------
Index: doc/documentation.tex
===================================================================
RCS file: /cactusdevcvs/CactusEinstein/IDAxiBrillBH/doc/documentation.tex,v
retrieving revision 1.5
diff -u -r1.5 documentation.tex
--- doc/documentation.tex 30 May 2004 08:22:55 -0000 1.5
+++ doc/documentation.tex 20 May 2005 10:20:02 -0000
@@ -9,7 +9,7 @@
\begin{document}
\title{IDAxiBrillBH}
-\author{Paul Walker, Steve Brandt}
+\author{Paul Walker, Steve Brandt,\\some cleanups by Jonathan Thornburg}
\date{$ $Date: 2004/05/30 08:22:55 $ $}
\maketitle
@@ -17,6 +17,8 @@
% Do not delete next line
% START CACTUS THORNGUIDE
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
\begin{abstract}
Thorn IDAxiBrillBH provides analytic initial data for a vacuum
black hole spacetime: a single Schwarzschild black hole in
@@ -25,6 +27,8 @@
extrinsic curvature.
\end{abstract}
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
\section{Purpose}
The pioneer, Bernstein, studied a single black hole which is
@@ -40,10 +44,10 @@
The Hamiltonian constraint reduces to
\begin{equation}
\hat \Delta \psi = \frac{1}{8}\psi \hat R,
-\label{eqn:conformal_hamiltonian}
+\label{IDAxiBrillBH/eqn:conformal-hamiltonian}
\end{equation}
-where $\hat \Delta$ is the covariant Laplacian and $\hat R$ is Ricci
-tensor for the conformal three metric. This form allows
+where $\hat \Delta$ is the covariant Laplacian and $\hat R$ is the
+Ricci tensor for the conformal three metric. This form allows
us to choose an arbitrary conformal three metric, and then solve an
elliptic equation for the conformal factor, therefore satisfying the
constraint equations ($K_{ij} = 0$ trivially satisfies the momentum
@@ -52,13 +56,13 @@
Bernstein extended this to the black hole spacetime. Using
spherical-polar coordinates, one can write the 3-metric,
\begin{equation}
-\label{eqn:sph-cood}
+\label{IDAxiBrillBH/eqn:sph-coord}
ds^2 = \psi^4 (e^{2q} (dr^2 + r^2 d \theta^2) + r^2 \sin \theta d
\phi^2),
\end{equation}
where $q$ is the Brill ``packet'' which takes some functional form.
-Using this ansatz with (\ref{eqn:conformal_hamiltonian}) leads to
-an elliptic equation for $\psi$ which must be solved
+Using this ansatz with (\ref{IDAxiBrillBH/eqn:conformal-hamiltonian})
+leads to an elliptic equation for $\psi$ which must be solved
numerically. Applying the isometry condition on $\psi$ at a finite
radius, and applying $M/2r$ falloff conditions on $\psi$ at the
outer boundary (the ``Robin'' condition), along with a packet which
@@ -68,15 +72,16 @@
Schwarzschild solution. The typical $q$ function used in
axisymmetry, and considered here in the non-rotating case, is
\begin{equation}
-q = Q_0 \sin \theta^n \left [ \exp\left(\frac{\eta -
+\label{IDAxiBrillBH/eqn:Q}
+q = Q_0 \sin^n \theta \left [ \exp\left(\frac{\eta -
\eta_0^2}{\sigma^2}\right ) + \exp\left(\frac{\eta +
\eta_0^2}{\sigma^2}\right ) \right ].
\end{equation}
Note regularity along the axis requires that the exponent $n$ must be
even. Choose a logarithmic radial coordinate $\eta$, which is
related to the
-asymptoticlly flat coordinate $r$ by $\eta = ln (2r/m)$, where $m$ is a
-scale parameter. One can rewrite (\ref{eqn:sph-cood}) as
+asymptoticlly flat coordinate $r$ by $\eta = \ln (2r/m)$, where $m$ is a
+scale parameter. One can rewrite (\ref{IDAxiBrillBH/eqn:sph-coord}) as
\begin{equation}
ds^2 = \psi(\eta)^4 [ e^{2 q} (d \eta^2 + d\theta^2) + \sin^2
\theta d\phi^2].
@@ -86,20 +91,20 @@
logarithmic radial coordinate
\begin{equation}
-\label{eta_coord}
+\label{IDAxiBrillBH/eta-coord}
\eta = \ln{\frac{2r}{m}}.
\end{equation}
The scale parameter $m$ is equal to the mass of the Schwarzschild
black hole, if $q=0$. In this coordinate, the 3-metric is
\begin{equation}
-\label{eqn:metric_brill_eta}
+\label{IDAxiBrillBH/eqn:metric-brill-eta}
ds^2 = \tilde{\psi}^4 (e^{2q} (d\eta^2+d\theta^2)+\sin^2 \theta
d\phi^2),
\end{equation}
and the Schwarzschild solution is
\begin{equation}
-\label{eqn:psi}
+\label{IDAxiBrillBH/eqn:psi}
\tilde{\psi} = \sqrt{2M} \cosh (\frac{\eta}{2}).
\end{equation}
We also change the notation of $\psi$ for the conformal factor is same
@@ -108,7 +113,7 @@
$\psi$ differ by a factor of $\sqrt{r}$. The Hamiltonian
constraint is
\begin{equation}
-\label{eqn:ham}
+\label{IDAxiBrillBH/eqn:ham}
\frac{\partial^2 \tilde{\psi}}{\partial \eta^2} + \frac{\partial^2
\tilde{\psi}}{\partial \theta^2} + \cot \theta \frac{\partial
\tilde{\psi}}{\partial \theta} = - \frac{1}{4} \tilde{\psi}
@@ -122,14 +127,14 @@
\delta \tilde{\psi} & = & \tilde{\psi}+\tilde{\psi}_0 \\
& = & \tilde{\psi}-\sqrt{2m} \cosh(\frac{\eta}{2}).
\end{eqnarray}
-to the equation~(\ref{eqn:ham}), then we can linearize it as
+to the equation~(\ref{IDAxiBrillBH/eqn:ham}), then we can linearize it as
\begin{equation}
\frac{\partial^2 \delta\tilde{\psi}}{\partial \eta^2} + \frac{\partial^2
\delta\tilde{\psi}}{\partial \theta^2} + \cot \theta \frac{\partial
\delta\tilde{\psi}}{\partial \theta} = - \frac{1}{4}
(\delta\tilde{\psi} + \tilde{\psi}_0) (\frac{\partial^2 q}{\partial
\eta^2} + \frac{\partial^2 q}{\partial \theta^2} -1).
-\label{eqn:ham_linear}
+\label{IDAxiBrillBH/eqn:ham-linear}
\end{equation}
For the boundary conditions, we use for the inner boundary condition
an isometry condition:
@@ -142,19 +147,84 @@
\tilde{\psi})|_{\eta=\eta_{max}} = 0.
\end{equation}
-% [[ DPR: What is this: ?? ]]
-%This thorn provides
-% \begin{enumerate}
-% \item CactusEinstein
-% \end{enumerate}
-
-\section{Comments}
-
-We calculate equation~(\ref{eqn:ham_linear}) with spherical
-coordinates. However, Cactus needs Cartesian coordinates. Therefore,
-we interpolate $\psi$ to the Cartesian grid by using an interpolator.
-Note that the interpolator has linear, quadratic, and cubic
-interpolation.
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+\section{2-D Grid and Interpolation Parameters}
+
+This thorn solves equation~(\ref{IDAxiBrillBH/eqn:ham-linear}) on a 2-D
+$(\eta,\theta)$ grid. However, Cactus needs a 3-D grid, typically with
+Cartesian coordinates. Therefore, this thorn interpolate $\psi$ and its
+$(\eta,\theta)$ derivatives to the Cartesian grid.
+
+The parameters \verb|neta| and \verb|nq| specify the resolution of
+this thorn's 2-D grid in $\eta$ and $\theta$ respectively.%%%
+\footnote{%%%
+ Internally, this thorn uses ``$q$'' to refer
+ to $\theta$ in Fortran code, with the $q$ function
+ of~$(\protect\ref{IDAxiBrillBH/eqn:Q})$ being hidden
+ in the Mathematica files (and not present in the Fortran
+ code). Noone seems to know \emph{why} the code does
+ things this way\dots{} Unfortunately, this renaming
+ has leaked out into the parameter names\dots
+ }%%%
+{} The default values are a reasonable starting point, but you may
+need to increase them substantially if you need very high accuracy
+(very small constraint violations).
+
+To help judge what resolution may be needed, this thorn has an option
+to write out $\psi(\eta)$ and $\psi$ on the 2-D grid to an ASCII data file
+where it can be examined and/or plotted. To do this, set the Boolean
+parameter \verb|output_psi2D|, and possibly also the string parameter
+\verb|output_psi2D_file_name|.
+
+This thorn uses the standard Cactus \verb|CCTK_InterpLocalUniform()|
+local interpolation system for this interpolation. The interpolation
+operator is specified with the \verb|interpolator_name| parameter
+(this defaults to \verb|"uniform cartesian"|, the interpolation
+operator provided by thorn \textbf{CactusBase/LocalInterp}).
+
+The interpolation order and/or other parameters can be specified
+in either of two ways:%%%
+\footnote{%%%
+ Notice that, for historical reasons, the
+ interpolation parameter names are a bit
+ inconsistent: \texttt{interpolat\underline{ion}\_order}
+ versus \texttt{interpolat\underline{or}\_name}
+ and \texttt{interpolat\underline{or}\_pars}.
+ }%%%
+\begin{itemize}
+\item The integer parameter \verb|interpolation_order| may be
+ used directly to specify the interpolation order.
+\item More generally, the string parameter \verb|interpolator_pars|
+ may be set to any nonempty string (it defaults to the empty string).
+ If this is done, this parameter overrides \verb|interpolation_order|,
+ and explicitly specifies a parameter string for the interpolator.
+\end{itemize}
+Note that the default interpolator parameters specify \emph{linear}
+interpolation. This is rather inaccurate, and (due to aliasing effects
+between the 2-D and 3-D grids) will give a fair bit of noise in the
+metric components. You may want to specify a higher-order interpolator
+to reduce this noise.
+
+For example, for one test series where I (JT) needed very accurate
+initial data (I wanted the initial-data errors to be much less than
+the errors from 4th~order finite differencing on the 3-D Cactus grid),
+I had to use a resolution of $1000$ in $\eta$ and $2000$ in $\theta$,
+together with either 4th~order Lagrange or 3rd~order Hermite interpolation
+(provided by thorn \textbf{AEIThorns/AEILocalInterp}) to get sufficient
+accuracy.
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+\section{Debugging Parameters}
+
+This thorn has options to print very detailed debugging information
+about internal quantities at selected grid points. This is enabled
+by setting the integer parameter \verb|debug| to a positive value
+(the default is $0$, which means no debugging output). See
+\verb|param.ccl| and the source code \verb|src/IDAxiBrillBH.F| for details.
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\bibliographystyle{prsty}
\begin{thebibliography}{10}
@@ -167,6 +237,8 @@
K. Camarda, Ph.D thesis University of Illinois Urbana-Champaign, (1998)
\end{thebibliography}
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
% Do not delete next line
% END CACTUS THORNGUIDE
-------------- next part --------------
Things that Could be Improved in This Thorn
===========================================
This thorn works with non-unigrid Cactus, but it doesn't work very
efficiently: it does the 2D elliptic setup and solve all over again
for each local Cactus (3D) grid. It would be more efficient to move the
2D elliptic stuff out into a separately scheduled routine which would
be done once (for Carpet, it would be scheduled in level or global
mode), and have only the interpolation to the 3D grid done for each
local Cactus grid.
Part of this same change might also be to store the 2D elliptic results
in Cactus grid arrays rather than native Fortran 90 arrays. This would
let us use standard Cactus I/O methods to output the 2D elliptic results
instead of the hard-wired code we use right now.
More information about the Developers
mailing list