aboutsummaryrefslogtreecommitdiff
diff options
context:
space:
mode:
authorryoji <ryoji@971fb155-194f-0410-9daf-e2eca44e59f5>1999-12-01 13:36:25 +0000
committerryoji <ryoji@971fb155-194f-0410-9daf-e2eca44e59f5>1999-12-01 13:36:25 +0000
commit6b266c308ed6ed21b401e3a670cc017286d79036 (patch)
tree759c532bf74cfc87e989c4ff2d2ab01e60b42274
parentf50dd221f56adcee042c811725c2dae63edf23c0 (diff)
This commit was generated by cvs2svn to compensate for changes in r2,
which included commits to RCS files with non-trunk default branches. git-svn-id: http://svn.einsteintoolkit.org/cactus/EinsteinInitialData/DistortedBHIVP/trunk@3 971fb155-194f-0410-9daf-e2eca44e59f5
Diffstat
-rw-r--r--README20
-rw-r--r--interface.ccl9
-rw-r--r--param.ccl65
-rw-r--r--schedule.ccl14
-rw-r--r--src/DistortedBHIVP.F665
-rw-r--r--src/DistortedBHIVP.F.sav681
-rw-r--r--src/Stab3d.F280
-rw-r--r--src/bhbrill3d.m132
-rw-r--r--src/bhbrill3d.x126
-rw-r--r--src/gij.x117
-rw-r--r--src/interp3.F264
-rw-r--r--src/make.code.defn9
-rw-r--r--src/psi_1st_deriv.x20
-rw-r--r--src/psi_2nd_deriv.x71
-rw-r--r--src/qfunc.x29
-rwxr-xr-xtest/test_dbh.par55
16 files changed, 2557 insertions, 0 deletions
diff --git a/README b/README
new file mode 100644
index 0000000..183b23b
--- /dev/null
+++ b/README
@@ -0,0 +1,20 @@
+Cactus Code Thorn DistortedBHIVP
+Authors : ...
+Managed by : ... <...@...........>
+Version : ...
+CVS info : $Header$
+--------------------------------------------------------------------------
+
+1. Purpose of the thorn
+
+This thorn does ...
+
+2. Dependencies of the thorn
+
+This thorn additionally requires implementations and thorns ...
+
+3. Thorn distribution
+
+This thorn is available to ...
+
+4. Additional information
diff --git a/interface.ccl b/interface.ccl
new file mode 100644
index 0000000..dd29b48
--- /dev/null
+++ b/interface.ccl
@@ -0,0 +1,9 @@
+# Interface definition for thorn DistortedBHIVP
+# $Header$
+
+implements: distortedbhivp
+inherits: einstein
+
+public:
+
+
diff --git a/param.ccl b/param.ccl
new file mode 100644
index 0000000..5e581ea
--- /dev/null
+++ b/param.ccl
@@ -0,0 +1,65 @@
+# Parameter definitions for thorn DistortedBHIVP
+# $Header$
+
+shares:einstein
+
+EXTENDS KEYWORD initial_lapse ""
+{
+ "schwarz" :: "Set lapse to Schwarzschild"
+}
+
+EXTENDS KEYWORD initial_data ""
+{
+"distortedbh" :: "Non-Aix BH + Brill Wave IVP"
+}
+
+private:
+
+REAL amp "Brill wave amplitude"
+{
+ *:* :: "No restriction"
+} 0.1
+
+REAL eta0 "Brill wave center (in eta coords)"
+{
+ *:* :: "No restriction"
+} 0.0
+
+REAL c "Azimuthal dependence of Brill wave"
+{
+ *:* :: "No restriction"
+} 0.0
+
+REAL sigma "Brill wave width (in eta)"
+{
+ *:* :: "No restriction"
+} 1.0
+
+REAL etamax "Eta value for outer edge of grid"
+{
+ *:* :: "No restriction"
+} 5.0
+
+
+INT n "sin^n theta in brill wave"
+{
+ *:* :: "No restriction"
+} 2
+
+INT ne "Eta resolution for solve"
+{
+ *:* :: "No restriction"
+} 202
+
+INT nq "Theta resolution for solve"
+{
+ *:* :: "No restriction"
+} 54
+
+INT np "Eta resolution for solve"
+{
+ *:* :: "No restriction"
+} 5
+
+
+
diff --git a/schedule.ccl b/schedule.ccl
new file mode 100644
index 0000000..c679a46
--- /dev/null
+++ b/schedule.ccl
@@ -0,0 +1,14 @@
+# Schedule definitions for thorn DistortedBHIVP
+# $Header$
+
+if (CCTK_Equals(initial_data,"distortedbh"))
+{
+ schedule DistortedBHIVP at CCTK_INITIAL
+ {
+ LANG: Fortran
+ } "Construct DistortedBHIVP"
+}
+
+
+
+
diff --git a/src/DistortedBHIVP.F b/src/DistortedBHIVP.F
new file mode 100644
index 0000000..9682cf4
--- /dev/null
+++ b/src/DistortedBHIVP.F
@@ -0,0 +1,665 @@
+c /*@@
+c @file DistortedBHIVP.F
+c @date
+c @author
+c @desc
+c
+c @enddesc
+c@@ */
+
+#include "cctk.h"
+#include "cctk_parameters.h"
+#include "cctk_arguments.h"
+
+c /*@@
+c @routine DistortedBHIVP
+c @date
+c @author
+c @desc
+c
+c @enddesc
+c @calls
+c @calledby
+c @history
+c
+c @endhistory
+c@@ */
+
+c Need include file from Einstein
+#include "CactusEinstein/Einstein/src/Einstein.h"
+
+ subroutine DistortedBHIVP(CCTK_FARGUMENTS)
+
+ implicit none
+
+ DECLARE_CCTK_FARGUMENTS
+ DECLARE_CCTK_PARAMETERS
+
+c Perhaps this and others should go into cctk.h
+ integer CCTK_Equals
+
+ real*8 :: deta,dq,dphi
+ real*8, allocatable :: ac(:,:,:),ae(:,:,:),aw(:,:,:),an(:,:,:),
+ $ as(:,:,:),aq(:,:,:),ab(:,:,:),rhs(:,:,:),
+ $ qf(:,:,:),qfetaeta(:,:,:),qfqq(:,:,:),qfphi(:,:,:),
+ $ qfphiphi(:,:,:),
+ $ psisph(:,:,:),psiprim(:,:,:),detapsisph(:,:,:),
+ $ dqpsisph(:,:,:),dphipsisph(:,:,:),detaetapsisph(:,:,:),
+ $ detaqpsisph(:,:,:),detaphipsisph(:,:,:),dqqpsisph(:,:,:),
+ $ dqphipsisph(:,:,:),dphiphipsisph(:,:,:),
+ $ psi3d(:,:,:),detapsi3d(:,:,:),dqpsi3d(:,:,:),
+ $ dphipsi3d(:,:,:),detaetapsi3d(:,:,:),detaqpsi3d(:,:,:),
+ $ detaphipsi3d(:,:,:),dqqpsi3d(:,:,:),dqphipsi3d(:,:,:),
+ $ dphiphipsi3d(:,:,:)
+ real*8, allocatable :: etagrd(:),qgrd(:),phigrd(:)
+ real*8, allocatable :: eta(:,:,:),abseta(:,:,:),sign_eta(:,:,:),
+ $ q(:,:,:),phi(:,:,:)
+ real*8 o1,o2,o3,o4,o5,o6,o7,o8,o9,o10,o11,o12,o13,o14,o15,o16,o17,
+ $ o18,o19,o20,o21,o22,o23,o24,o25,o26,o27,o28,o29,o30,o31,o32,
+ $ o33,o34,o35,o36,o37,o38,o39,o40,o41,o42,o43,o44,o45,o46,o47,
+ $ o48,o49,o50,o51,o52,o53,o54,o55,o56,o57,o58,o59,o60,o61,o62,
+ $ o63,o64,o65,o66,o67
+ real*8 rmax,adm
+ real*8,parameter :: dbh_eps = 1.0d-9
+ real*8 r,pi
+ integer :: nx,ny,nz
+ integer i,j,k,ier,nquads,nocts,order
+ integer npoints,handle,ierror
+
+ conformal_state = CONFORMAL_METRIC
+
+ 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 Distorted Schwarzchild BH parameters
+
+ print *,"Brill wave + Distorted BH solve"
+ write(*,123)amp,eta0,c,sigma,n
+ print*,'etamax=',etamax
+ 123 format(1x, 'Pars: amp',f8.5,' eta0',f8.5,' c',f8.5,' sigma',f8.5,' n',i3)
+
+c Sovle on this sized cartesian grid
+c 3D grid size NE x NT x NP
+c Add 2 zones for eta coordinate and 4 for theta
+c and phi coordenate.
+ ne = ne + 2
+ nq = nq + 2
+ np = np + 2
+c
+ allocate(ac(ne,nq,np),ae(ne,nq,np),aw(ne,nq,np),an(ne,nq,np),
+ $ as(ne,nq,np),aq(ne,nq,np),ab(ne,nq,np),rhs(ne,nq,np),
+ $ qf(ne,nq,np),qfetaeta(ne,nq,np),qfqq(ne,nq,np),
+ $ qfphi(ne,nq,np),qfphiphi(ne,nq,np),
+ $ psisph(ne,nq,np),psiprim(ne,nq,np),detapsisph(ne,nq,np),
+ $ dqpsisph(ne,nq,np),dphipsisph(ne,nq,np),
+ $ detaetapsisph(ne,nq,np),detaqpsisph(ne,nq,np),
+ $ detaphipsisph(ne,nq,np),dqqpsisph(ne,nq,np),
+ $ dqphipsisph(ne,nq,np),dphiphipsisph(ne,nq,np),
+ $ psi3d(nx,ny,nz),detapsi3d(nx,ny,nz),dqpsi3d(nx,ny,nz),
+ $ dphipsi3d(nx,ny,nz),detaetapsi3d(nx,ny,nz),
+ $ detaqpsi3d(nx,ny,nz),detaphipsi3d(nx,ny,nz),
+ $ dqqpsi3d(nx,ny,nz),dqphipsi3d(nx,ny,nz),
+ $ dphiphipsi3d(nx,ny,nz))
+ allocate(etagrd(ne),qgrd(nq),phigrd(np))
+ allocate(eta(nx,ny,nz),abseta(nx,ny,nz),sign_eta(nx,ny,nz),
+ $ q(nx,ny,nz),phi(nx,ny,nz))
+c
+c Initialize some array
+c
+ nocts = 4
+ nquads = 2
+ dphi = nocts*0.5*pi/(np-2)
+ dq = nquads*0.5*pi/(nq-2)
+ deta = etamax/(ne-3)
+
+ do k = 1,np
+ phigrd(k) = (k-1.5)*dphi
+ enddo
+ do j = 1,nq
+ qgrd(j) = (j-1.5)*dq
+ enddo
+ do i=1,ne
+ etagrd(i) = (i-2)*deta
+ enddo
+c
+c Initialize q-function and its derivatives: should be generalized
+c
+ do k = 1,np
+ do j = 1,nq
+ do i = 1,ne
+#include "CactusEinstein/DistortedBHIVP/src/qfunc.x"
+ enddo
+ enddo
+ enddo
+
+c
+c Initialize psi to the Schwarzschild solution:
+c
+ psiprim = 0.
+
+ do k = 1,np
+ do j = 1,nq
+ do i = 1,ne
+ psisph(i,j,k) = 2.*cosh(0.5*etagrd(i))
+ enddo
+ enddo
+ enddo
+c
+c Initialize stencil coefficients:
+c
+ ac = 0.
+ ae = 0.
+ aw = 0.
+ an = 0.
+ as = 0.
+ aq = 0.
+ ab = 0.
+ rhs = 0.
+c
+
+ do k = 2,np-1
+ do j = 2,nq-1
+ do i = 2,ne-1
+ ac(i,j,k) = -2./deta**2-2./dq**2-2.*exp(2.*
+ & qf(i,j,k))/(dphi**2*
+ & sin(qgrd(j))**2)+0.25*
+ & (qfetaeta(i,j,k)+qfqq(i,j,k)+2.*
+ & exp(2.*qf(i,j,k))*qfphiphi(i,j,k)/
+ & sin(qgrd(j))**2+3.*exp(2.*
+ & qf(i,j,k))*qfphi(i,j,k)**2/
+ & sin(qgrd(j))**2-1.)
+ ae(i,j,k) = 1./deta**2
+ aw(i,j,k) = 1./deta**2
+ an(i,j,k) = 1./dq**2+0.5/(dq*tan(qgrd(j)))
+
+ as(i,j,k) = 1./dq**2-0.5/(dq*tan(qgrd(j)))
+
+ aq(i,j,k) = exp(2.*qf(i,j,k))/(dphi**2*
+ & sin(qgrd(j))**2)+exp(2.*
+ & qf(i,j,k))*qfphi(i,j,k)/(dphi*
+ & sin(qgrd(j))**2)
+
+ ab(i,j,k) = exp(2.*qf(i,j,k))/(dphi**2*
+ & sin(qgrd(j))**2)-exp(2.*
+ & qf(i,j,k))*qfphi(i,j,k)/(dphi*
+ & sin(qgrd(j))**2)
+
+ rhs(i,j,k) = -0.25*(qfetaeta(i,j,k)+
+ & qfqq(i,j,k)+2.*exp(2.*qf(i,j,k))*
+ & qfphiphi(i,j,k)/sin(qgrd(j))**2+3.*
+ & exp(2.*qf(i,j,k))*qfphi(i,j,k)**2/
+ & sin(qgrd(j))**2)*psisph(i,j,k)
+ enddo
+ enddo
+ enddo
+c
+c Apply boundary conditions to the faces of the cube:
+c
+c i=2:
+ do k = 3,np-2
+ do j = 3,nq-2
+ ae(2,j,k) = ae(2,j,k) + aw(2,j,k)
+ aw(2,j,k) = 0.
+c
+c i=ne-1:
+ ac(ne-1,j,k) = ac(ne-1,j,k)+4.*ae(ne-1,j,k)/(3.+deta)
+ aw(ne-1,j,k) = aw(ne-1,j,k) - ae(ne-1,j,k)/(3.+deta)
+ ae(ne-1,j,k) = 0.
+ enddo
+ enddo
+c
+c j=2:
+ do k = 3,np-2
+ do i = 3,ne-2
+ ac(i,2,k) = ac(i,2,k) + as(i,2,k)
+ as(i,2,k) = 0.
+c
+c j=nq-1:
+ ac(i,nq-1,k) = ac(i,nq-1,k) + an(i,nq-1,k)
+ an(i,nq-1,k) = 0.
+ enddo
+ enddo
+c
+c k=2:
+ do j = 3,nq-2
+ do i = 3,ne-2
+ ac(i,j,2) = ac(i,j,2) + ab(i,j,2)
+ ab(i,j,2) = 0.
+c
+c k=np-1:
+ ac(i,j,np-1) = ac(i,j,np-1) + aq(i,j,np-1)
+ aq(i,j,np-1) = 0.
+ enddo
+ enddo
+c
+c Apply boundary conditions to the edges of the cube:
+c
+c i=2, j=2:
+ do k = 3,np-2
+ ae(2,2,k) = ae(2,2,k) + aw(2,2,k)
+ ac(2,2,k) = ac(2,2,k) + as(2,2,k)
+ aw(2,2,k) = 0.
+ as(2,2,k) = 0.
+c
+c i=ne-1, j=2:
+ aw(ne-1,2,k) = aw(ne-1,2,k) - ae(ne-1,2,k)/(3.+deta)
+ ac(ne-1,2,k) = ac(ne-1,2,k) + as(ne-1,2,k) +
+ & 4.*ae(ne-1,2,k)/(3.+deta)
+ ae(ne-1,2,k) = 0.
+ as(ne-1,2,k) = 0.
+c
+c i=2, j=nq-1:
+ ae(2,nq-1,k) = ae(2,nq-1,k) + aw(2,nq-1,k)
+ ac(2,nq-1,k) = ac(2,nq-1,k) + an(2,nq-1,k)
+ aw(2,nq-1,k) = 0.
+ an(2,nq-1,k) = 0.
+c
+c i=ne-1, j=nq-1:
+ aw(ne-1,nq-1,k) = aw(ne-1,nq-1,k) - ae(ne-1,nq-1,k)/
+ & (3.+deta)
+ ac(ne-1,nq-1,k) = ac(ne-1,nq-1,k) + an(ne-1,nq-1,k) +
+ & 4.*ae(ne-1,nq-1,k)/(3.+deta)
+ ae(ne-1,nq-1,k) = 0.
+ an(ne-1,nq-1,k) = 0.
+ enddo
+c
+c i=2, k=2:
+ do j = 3,nq-2
+ ae(2,j,2) = ae(2,j,2) + aw(2,j,2)
+ ac(2,j,2) = ac(2,j,2) + ab(2,j,2)
+ aw(2,j,2) = 0.
+ ab(2,j,2) = 0.
+c
+c i=ne-1, k=2:
+ aw(ne-1,j,2) = aw(ne-1,j,2) - ae(ne-1,j,2)/(3.+deta)
+ ac(ne-1,j,2) = ac(ne-1,j,2) + ab(ne-1,j,2) +
+ & 4.*ae(ne-1,j,2)/(3.+deta)
+ ae(ne-1,j,2) = 0.
+ ab(ne-1,j,2) = 0.
+c
+c i=2, k=np-1:
+ ae(2,j,np-1) = ae(2,j,np-1) + aw(2,j,np-1)
+ ac(2,j,np-1) = ac(2,j,np-1) + aq(2,j,np-1)
+ aw(2,j,np-1) = 0.
+ aq(2,j,np-1) = 0.
+c
+c i=ne-1, k=np-1:
+ aw(ne-1,j,np-1) = aw(ne-1,j,np-1) - ae(ne-1,j,np-1)/
+ & (3.+deta)
+ ac(ne-1,j,np-1) = ac(ne-1,j,np-1) + aq(ne-1,j,np-1) +
+ & 4.*ae(ne-1,j,np-1)/(3.+deta)
+ ae(ne-1,j,np-1) = 0.
+ aq(ne-1,j,np-1) = 0.
+ enddo
+c
+c j=2, k=2:
+ do i = 3,ne-2
+ ac(i,2,2) = ac(i,2,2) + as(i,2,2) + ab(i,2,2)
+ as(i,2,2) = 0.
+ ab(i,2,2) = 0.
+c
+c j=nq-1, k=2:
+ ac(i,nq-1,2) = ac(i,nq-1,2) + an(i,nq-1,2) +
+ & ab(i,nq-1,2)
+ an(i,nq-1,2) = 0.
+ ab(i,nq-1,2) = 0.
+c
+c j=2, k=np-1:
+ ac(i,2,np-1) = ac(i,2,np-1) + as(i,2,np-1) +
+ & aq(i,2,np-1)
+ as(i,2,np-1) = 0.
+ aq(i,2,np-1) = 0.
+c
+c j=nq-1, k=np-1:
+ ac(i,nq-1,np-1) = ac(i,nq-1,np-1) + an(i,nq-1,np-1) +
+ & aq(i,nq-1,np-1)
+ an(i,nq-1,np-1) = 0.
+ aq(i,nq-1,np-1) = 0.
+ enddo
+c
+c Apply boundary conditions to the corners of the cube:
+c
+c i=2, j=2, k=2:
+ ae(2,2,2) = ae(2,2,2) + aw(2,2,2)
+ ac(2,2,2) = ac(2,2,2) + as(2,2,2) + ab(2,2,2)
+ aw(2,2,2) = 0.
+ as(2,2,2) = 0.
+ ab(2,2,2) = 0.
+c
+c i=ne-1, j=2, k=2:
+ aw(ne-1,2,2) = aw(ne-1,2,2) - ae(ne-1,2,2)/(3.+deta)
+ ac(ne-1,2,2) = ac(ne-1,2,2) + as(ne-1,2,2) + ab(ne-1,2,2) +
+ & 4.*ae(ne-1,2,2)/(3.+deta)
+ ae(ne-1,2,2) = 0.
+ as(ne-1,2,2) = 0.
+ ab(ne-1,2,2) = 0.
+c
+c i=2, j=nq-1, k=2:
+ ae(2,nq-1,2) = ae(2,nq-1,2) + aw(2,nq-1,2)
+ ac(2,nq-1,2) = ac(2,nq-1,2) + an(2,nq-1,2) + ab(2,nq-1,2)
+ aw(2,nq-1,2) = 0.
+ an(2,nq-1,2) = 0.
+ ab(2,nq-1,2) = 0.
+c
+c i=2, j=2, k=np-1:
+ ae(2,2,np-1) = ae(2,2,np-1) + aw(2,2,np-1)
+ ac(2,2,np-1) = ac(2,2,np-1) + as(2,2,np-1) + aq(2,2,np-1)
+ aw(2,2,np-1) = 0.
+ as(2,2,np-1) = 0.
+ aq(2,2,np-1) = 0.
+c
+c i=ne-1, j=nq-1, k=2:
+ aw(ne-1,nq-1,2) = aw(ne-1,nq-1,2) - ae(ne-1,nq-1,2)/(3.+deta)
+ ac(ne-1,nq-1,2) = ac(ne-1,nq-1,2) + an(ne-1,nq-1,2) + ab(ne-1,nq-1,2) +
+ & 4.*ae(ne-1,nq-1,2)/(3.+deta)
+ ae(ne-1,nq-1,2) = 0.
+ an(ne-1,nq-1,2) = 0.
+ ab(ne-1,nq-1,2) = 0.
+c
+c i=ne-1, j=2, k=np-1:
+ aw(ne-1,2,np-1) = aw(ne-1,2,np-1) - ae(ne-1,2,np-1)/(3.+deta)
+ ac(ne-1,2,np-1) = ac(ne-1,2,np-1) + as(ne-1,2,np-1) + aq(ne-1,2,np-1) +
+ & 4.*ae(ne-1,2,np-1)/(3.+deta)
+ ae(ne-1,2,np-1) = 0.
+ as(ne-1,2,np-1) = 0.
+ aq(ne-1,2,np-1) = 0.
+c
+c i=2, j=nq-1, k=np-1:
+ ae(2,nq-1,np-1) = ae(2,nq-1,np-1) + aw(2,nq-1,np-1)
+ ac(2,nq-1,np-1) = ac(2,nq-1,np-1) + an(2,nq-1,np-1) + aq(2,nq-1,np-1)
+ aw(2,nq-1,np-1) = 0.
+ an(2,nq-1,np-1) = 0.
+ aq(2,nq-1,np-1) = 0.
+c
+c i=ne-1, j=nq-1, k=np-1:
+ aw(ne-1,nq-1,np-1) = aw(ne-1,nq-1,np-1) - ae(ne-1,nq-1,np-1)/(3.+deta)
+ ac(ne-1,nq-1,np-1) = ac(ne-1,nq-1,np-1) + an(ne-1,nq-1,np-1) +
+ & aq(ne-1,nq-1,np-1) + 4.*ae(ne-1,nq-1,np-1)/(3.+deta)
+ ae(ne-1,nq-1,np-1) = 0.
+ an(ne-1,nq-1,np-1) = 0.
+ aq(ne-1,nq-1,np-1) = 0.
+c
+c Solve for psi:
+c
+ call bicgst3d(ac,ae,aw,an,as,aq,ab,psiprim,rhs,dbh_eps,rmax,ier,
+ $ ne,nq,np)
+c
+ if (rmax.gt.dbh_eps) then
+ write(*,*) '***WARNING: bicgst3d did not converge.'
+ endif
+ if (ier.eq.-1) then
+ write(*,*) '***WARNING: ier=-1'
+ endif
+ print *,'psiprim = ',maxval(psiprim),' ',minval(psiprim)
+c
+c Now, apply boundary conditions to psiprim:
+c
+ do k = 1,np
+ do j = 1,nq
+ psiprim(1,j,k) = psiprim(3,j,k)
+ psiprim(ne,j,k) = (4.*psiprim(ne-1,j,k)-psiprim(ne-2,j,k))/
+ $ (3.+deta)
+ enddo
+ enddo
+ do k = 1,np
+ do i = 1,ne
+ psiprim(i,1,k) = psiprim(i,2,k)
+ psiprim(i,nq,k) = psiprim(i,nq-1,k)
+ enddo
+ enddo
+ do j = 1,nq
+ do i = 1,ne
+ psiprim(i,j,1) = psiprim(i,j,2)
+ psiprim(i,j,np) = psiprim(i,j,np-1)
+ enddo
+ enddo
+c
+c Here, compute the derivatives of the spherical conformal factor
+c
+c goto 110
+
+ do k = 1,np
+ do j = 1,nq
+ do i = 2,ne-1
+ detapsisph(i,j,k)=0.5*(psiprim(i+1,j,k)-psiprim(i-1,j,k))
+ $ /deta + sinh(0.5*etagrd(i))
+ enddo
+ detapsisph(1,j,k) = -detapsisph(3,j,k)
+ enddo
+ enddo
+c
+ do k = 1,np
+ do j = 2,nq-1
+ do i = 1,ne
+ dqpsisph(i,j,k)=0.5*(psiprim(i,j+1,k)-psiprim(i,j-1,k))/
+ $ dq
+ enddo
+ enddo
+ do i = 1,ne
+ dqpsisph(i,1,k) = -dqpsisph(i,2,k)
+ dqpsisph(i,nq,k) = -dqpsisph(i,nq-1,k)
+ enddo
+ enddo
+c
+ do k = 2,np-1
+ do j = 1,nq
+ do i = 1,ne
+ dphipsisph(i,j,k)=0.5*(psiprim(i,j,k+1)-psiprim(i,j,k-1))
+ $ /dphi
+ enddo
+ enddo
+ enddo
+ do j = 1,nq
+ do i = 1,ne
+ dphipsisph(i,j,1) = -dphipsisph(i,j,2)
+ dphipsisph(i,j,np) = -dphipsisph(i,j,np-1)
+ enddo
+ enddo
+c
+ do k = 1,np
+ do j = 1,nq
+ do i = 2,ne-1
+ detaetapsisph(i,j,k)=(psiprim(i+1,j,k)-2.*psiprim(i,j,k)+
+ & psiprim(i-1,j,k))/deta**2+sqrt(0.25)*
+ & cosh(0.5*etagrd(i))
+ enddo
+ detaetapsisph(1,j,k) = detaetapsisph(3,j,k)
+ enddo
+ enddo
+c
+ do k = 1,np
+ do j = 2,nq-1
+ do i = 1,ne
+ detaqpsisph(i,j,k)=0.5*(detapsisph(i,j+1,k)-detapsisph(i,
+ $ j-1,k))/dq
+ enddo
+ enddo
+ do i = 1,ne
+ detaqpsisph(i,1,k) = -detaqpsisph(i,2,k)
+ detaqpsisph(i,nq,k) = -detaqpsisph(i,nq-1,k)
+ enddo
+ enddo
+c
+ do k = 2,np-1
+ do j = 1,nq
+ do i = 1,ne
+ detaphipsisph(i,j,k)=0.5*(detapsisph(i,j,k+1)-detapsisph(
+ $ i,j,k-1))/dphi
+ enddo
+ enddo
+ enddo
+ do j = 1,nq
+ do i = 1,ne
+ detaphipsisph(i,j,1) = -detaphipsisph(i,j,2)
+ detaphipsisph(i,j,np) = -detaphipsisph(i,j,np-1)
+ enddo
+ enddo
+c
+ do k = 1,np
+ do j = 2,nq-1
+ do i = 1,ne
+ dqqpsisph(i,j,k)=0.5*(dqpsisph(i,j+1,k)-dqpsisph(i,j-1,k))/
+ $ dq
+ enddo
+ enddo
+ do i = 1,ne
+ dqqpsisph(i,1,k) = dqqpsisph(i,2,k)
+ dqqpsisph(i,nq,k) = dqqpsisph(i,nq-1,k)
+ enddo
+ enddo
+c
+ do k = 2,np-1
+ do j = 1,nq
+ do i = 1,ne
+ dqphipsisph(i,j,k)=0.5*(dqpsisph(i,j,k+1)-dqpsisph(i,j,k-1))/
+ $ dphi
+ enddo
+ enddo
+ enddo
+ do j = 1,nq
+ do i = 1,ne
+ dqphipsisph(i,j,1) = -dqphipsisph(i,j,2)
+ dqphipsisph(i,j,np) = -dqphipsisph(i,j,np-1)
+ enddo
+ enddo
+c
+ do k = 2,np-1
+ do j = 1,nq
+ do i = 1,ne
+ dphiphipsisph(i,j,k)=0.5*(dphipsisph(i,j,k+1)-
+ $ dphipsisph(i,j,k-1))/dphi
+ enddo
+ enddo
+ enddo
+ do j = 1,nq
+ do i = 1,ne
+ dphiphipsisph(i,j,1) = dphiphipsisph(i,j,2)
+ dphiphipsisph(i,j,np) = dphiphipsisph(i,j,np-1)
+ enddo
+ enddo
+c
+ do k = 1,np
+ do j = 1,nq
+ psisph(:,j,k)=psiprim(:,j,k)+2.0*cosh(0.5*etagrd)
+ enddo
+ enddo
+
+c
+c Now compute on the Cartesian coordinate.
+c
+c Compute eta,q,phi at the each points of cartesian grid
+
+ eta = 0.5d0*dlog(x**2+y**2+z**2)
+ abseta = abs(eta)
+ q = datan2(sqrt(x**2+y**2),z)
+ phi = datan2(y,x)
+
+ do k = 1,nz
+ do j = 1,ny
+ do i = 1,nx
+
+ if(eta(i,j,k) .lt. 0)then
+ sign_eta(i,j,k) = -1.0d0
+ else
+ sign_eta(i,j,k) = 1.0d0
+ endif
+
+ enddo
+ enddo
+ enddo
+
+ call CCTK_GetInterpHandle (handle, "simple_local")
+
+ npoints = nx*ny*nz
+
+ call CCTK_Interp (ierror,cctkGH,handle,npoints,3,10,10,
+ $ ne,nq,np,abseta,q,phi,
+ $ CCTK_VARIABLE_REAL,CCTK_VARIABLE_REAL,CCTK_VARIABLE_REAL,
+ $ etagrd(1),qgrd(1),phigrd(1)-pi,deta,dq,dphi,
+ $ psisph,detapsisph,dqpsisph,dphipsisph,detaetapsisph,
+ $ detaqpsisph,detaphipsisph,dqqpsisph,dqphipsisph,dphiphipsisph,
+ $ 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,
+ $ CCTK_VARIABLE_REAL,
+ $ psi3d,detapsi3d,dqpsi3d,dphipsi3d,detaetapsi3d,detaqpsi3d,
+ $ detaphipsi3d,dqqpsi3d,dqphipsi3d,dphiphipsi3d,
+ $ 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,
+ $ CCTK_VARIABLE_REAL)
+
+ psi = psi3d*exp(-0.5*eta)
+ detapsi3d = sign_eta*detapsi3d
+ detaqpsi3d = sign_eta*detaqpsi3d
+ detaphipsi3d = sign_eta*detaphipsi3d
+
+ do k=1,nz
+ do j=1,ny
+ do i=1,nx
+c psix = \partial psi / \partial x / psi
+#include "CactusEinstein/DistortedBHIVP/src/psi_1st_deriv.x"
+c psixx = \partial^2\psi / \partial x^2 / psi
+#include "CactusEinstein/DistortedBHIVP/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
+c
+ do k=1,nz
+ do j=1,ny
+ do i=1,nx
+#include "CactusEinstein/DistortedBHIVP/src/gij.x"
+ enddo
+ enddo
+ enddo
+
+c Courvature
+ kxx = 0.0d0
+ kxy = 0.0d0
+ kxz = 0.0d0
+ kyy = 0.0d0
+ kyz = 0.0d0
+ kzz = 0.0d0
+
+ 110 continue
+c Set ADM mass
+ i = ne-15
+ adm = 0.0
+ do k=2,np-1
+ do j=2,nq-1
+ adm=adm+(psisph(i,j,k)-(psisph(i+1,j,k)-psisph(i-1,j,k))/
+ $ deta)*exp(0.5*etagrd(i))
+ enddo
+ enddo
+ adm=adm/(nq-2)/(np-2)
+ print *,'ADM mass: ',adm
+ if (CCTK_Equals(initial_lapse,"schwarz")==1) then
+ write (*,*)"Initial with schwarzschild-like lapse"
+ write (*,*)"using alp = (2.*r - adm)/(2.*r+adm)."
+ alp = (2.*r - adm)/(2.*r+adm)
+ endif
+
+ deallocate(ac,ae,aw,an,as,aq,ab,rhs,
+ $ qf,qfetaeta,qfqq,qfphi,qfphiphi,
+ $ psisph,psiprim,detapsisph,dqpsisph,dphipsisph,detaetapsisph,
+ $ detaqpsisph,detaphipsisph,dqqpsisph,dqphipsisph,
+ $ dphiphipsisph,
+ $ psi3d,detapsi3d,dqpsi3d,
+ $ dphipsi3d,detaetapsi3d,detaqpsi3d,
+ $ detaphipsi3d,dqqpsi3d,dqphipsi3d,
+ $ dphiphipsi3d)
+ deallocate(etagrd,qgrd,phigrd)
+ deallocate(eta,abseta,sign_eta,q,phi)
+
+ return
+ end
diff --git a/src/DistortedBHIVP.F.sav b/src/DistortedBHIVP.F.sav
new file mode 100644
index 0000000..f02c56b
--- /dev/null
+++ b/src/DistortedBHIVP.F.sav
@@ -0,0 +1,681 @@
+c /*@@
+c @file DistortedBHIVP.F
+c @date
+c @author
+c @desc
+c
+c @enddesc
+c@@ */
+
+#include "cctk.h"
+#include "cctk_parameters.h"
+#include "cctk_arguments.h"
+
+c /*@@
+c @routine DistortedBHIVP
+c @date
+c @author
+c @desc
+c
+c @enddesc
+c @calls
+c @calledby
+c @history
+c
+c @endhistory
+c@@ */
+
+c Need include file from Einstein
+#include "CactusEinstein/Einstein/src/Einstein.h"
+
+ subroutine DistortedBHIVP(CCTK_FARGUMENTS)
+
+ implicit none
+
+ DECLARE_CCTK_FARGUMENTS
+ DECLARE_CCTK_PARAMETERS
+
+c Perhaps this and others should go into cctk.h
+ integer CCTK_Equals
+
+ real*8 :: amp,eta0,sigma,c,etamax,deta,dq,dphi
+ integer :: ne,nq,np,n
+ real*8, allocatable :: ac(:,:,:),ae(:,:,:),aw(:,:,:),an(:,:,:),
+ $ as(:,:,:),aq(:,:,:),ab(:,:,:),rhs(:,:,:),
+ $ qf(:,:,:),qfetaeta(:,:,:),qfqq(:,:,:),qfphi(:,:,:),
+ $ qfphiphi(:,:,:),
+ $ psisph(:,:,:),psiprim(:,:,:),detapsisph(:,:,:),
+ $ dqpsisph(:,:,:),dphipsisph(:,:,:),detaetapsisph(:,:,:),
+ $ detaqpsisph(:,:,:),detaphipsisph(:,:,:),dqqpsisph(:,:,:),
+ $ dqphipsisph(:,:,:),dphiphipsisph(:,:,:),
+ $ psi3d(:,:,:),detapsi3d(:,:,:),dqpsi3d(:,:,:),
+ $ dphipsi3d(:,:,:),detaetapsi3d(:,:,:),detaqpsi3d(:,:,:),
+ $ detaphipsi3d(:,:,:),dqqpsi3d(:,:,:),dqphipsi3d(:,:,:),
+ $ dphiphipsi3d(:,:,:)
+ real*8, allocatable :: etagrd(:),qgrd(:),phigrd(:)
+ real*8, allocatable :: eta(:,:,:),abseta(:,:,:),sign_eta(:,:,:),
+ $ q(:,:,:),phi(:,:,:)
+ real*8 o1,o2,o3,o4,o5,o6,o7,o8,o9,o10,o11,o12,o13,o14,o15,o16,o17,
+ $ o18,o19,o20,o21,o22,o23,o24,o25,o26,o27,o28,o29,o30,o31,o32,
+ $ o33,o34,o35,o36,o37,o38,o39,o40,o41,o42,o43,o44,o45,o46,o47,
+ $ o48,o49,o50,o51,o52,o53,o54,o55,o56,o57,o58,o59,o60,o61,o62,
+ $ o63,o64,o65,o66,o67
+ real*8 rmax,adm
+ real*8,parameter :: dbh_eps = 1.0d-9
+ real*8 r,pi
+ integer :: nx,ny,nz
+ integer i,j,k,ier,nquads,nocts,order
+ integer npoints,handle,ierror
+
+ conformal_state = CONFORMAL_METRIC
+
+ 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 Distorted Schwarzchild BH parameters
+ amp = distortedbh_amp
+ eta0 = distortedbh_eta0
+ c = distortedbh_c
+ sigma = distortedbh_sigma
+ n = distortedbh_n
+ etamax = distortedbh_etamax
+
+ print *,"Brill wave + Distorted BH solve"
+ write(*,123)amp,eta0,c,sigma,n
+ print*,'etamax=',etamax
+ 123 format(1x, 'Pars: amp',f8.5,' eta0',f8.5,' c',f8.5,' sigma',f8.5,' n',i3)
+
+c Sovle on this sized cartesian grid
+c 3D grid size NE x NT x NP
+c Add 2 zones for eta coordinate and 4 for theta
+c and phi coordenate.
+ ne = distortedbh_ne + 2
+ nq = distortedbh_nq + 2
+ np = distortedbh_np + 2
+c
+ allocate(ac(ne,nq,np),ae(ne,nq,np),aw(ne,nq,np),an(ne,nq,np),
+ $ as(ne,nq,np),aq(ne,nq,np),ab(ne,nq,np),rhs(ne,nq,np),
+ $ qf(ne,nq,np),qfetaeta(ne,nq,np),qfqq(ne,nq,np),
+ $ qfphi(ne,nq,np),qfphiphi(ne,nq,np),
+ $ psisph(ne,nq,np),psiprim(ne,nq,np),detapsisph(ne,nq,np),
+ $ dqpsisph(ne,nq,np),dphipsisph(ne,nq,np),
+ $ detaetapsisph(ne,nq,np),detaqpsisph(ne,nq,np),
+ $ detaphipsisph(ne,nq,np),dqqpsisph(ne,nq,np),
+ $ dqphipsisph(ne,nq,np),dphiphipsisph(ne,nq,np),
+ $ psi3d(nx,ny,nz),detapsi3d(nx,ny,nz),dqpsi3d(nx,ny,nz),
+ $ dphipsi3d(nx,ny,nz),detaetapsi3d(nx,ny,nz),
+ $ detaqpsi3d(nx,ny,nz),detaphipsi3d(nx,ny,nz),
+ $ dqqpsi3d(nx,ny,nz),dqphipsi3d(nx,ny,nz),
+ $ dphiphipsi3d(nx,ny,nz))
+ allocate(etagrd(ne),qgrd(nq),phigrd(np))
+ allocate(eta(nx,ny,nz),abseta(nx,ny,nz),sign_eta(nx,ny,nz),
+ $ q(nx,ny,nz),phi(nx,ny,nz))
+c
+c Initialize some array
+c
+ nocts = 4
+ nquads = 2
+ dphi = nocts*0.5*pi/(np-2)
+ dq = nquads*0.5*pi/(nq-2)
+ deta = etamax/(ne-3)
+
+ do k = 1,np
+ phigrd(k) = (k-1.5)*dphi
+ do j = 1,nq
+ qgrd(j) = (j-1.5)*dq
+ do i=1,ne
+ etagrd(i) = (i-2)*deta
+ enddo
+ enddo
+ enddo
+c
+c Initialize q-function and its derivatives: should be generalized
+c
+ do k = 1,np
+ do j = 1,nq
+ do i = 1,ne
+#include "CactusEinstein/DistortedBHIVP/src/qfunc.x"
+ enddo
+ enddo
+ enddo
+
+c
+c Initialize psi to the Schwarzschild solution:
+c
+ do k = 1,np
+ do j = 1,nq
+ do i = 1,ne
+ psisph(i,j,k) = 2.*cosh(0.5*etagrd(i))
+ psiprim(i,j,k) = 0.
+ enddo
+ enddo
+ enddo
+c
+c Initialize stencil coefficients:
+c
+ do k = 1,np
+ do j = 1,nq
+ do i = 1,ne
+ ac(i,j,k) = 0.
+ ae(i,j,k) = 0.
+ aw(i,j,k) = 0.
+ an(i,j,k) = 0.
+ as(i,j,k) = 0.
+ aq(i,j,k) = 0.
+ ab(i,j,k) = 0.
+ rhs(i,j,k) = 0.
+ enddo
+ enddo
+ enddo
+c
+ do k = 2,np-1
+ do j = 2,nq-1
+ do i = 2,ne-1
+ ac(i,j,k) = -2./deta**2-2./dq**2-2.*exp(2.*
+ & qf(i,j,k))/(dphi**2*
+ & sin(qgrd(j))**2)+0.25*
+ & (qfetaeta(i,j,k)+qfqq(i,j,k)+2.*
+ & exp(2.*qf(i,j,k))*qfphiphi(i,j,k)/
+ & sin(qgrd(j))**2+3.*exp(2.*
+ & qf(i,j,k))*qfphi(i,j,k)**2/
+ & sin(qgrd(j))**2-1.)
+
+ ae(i,j,k) = 1./deta**2
+
+ aw(i,j,k) = 1./deta**2
+
+ an(i,j,k) = 1./dq**2+0.5/(dq*tan(qgrd(j)))
+
+ as(i,j,k) = 1./dq**2-0.5/(dq*tan(qgrd(j)))
+
+ aq(i,j,k) = exp(2.*qf(i,j,k))/(dphi**2*
+ & sin(qgrd(j))**2)+exp(2.*
+ & qf(i,j,k))*qfphi(i,j,k)/(dphi*
+ & sin(qgrd(j))**2)
+
+ ab(i,j,k) = exp(2.*qf(i,j,k))/(dphi**2*
+ & sin(qgrd(j))**2)-exp(2.*
+ & qf(i,j,k))*qfphi(i,j,k)/(dphi*
+ & sin(qgrd(j))**2)
+
+ rhs(i,j,k) = -0.25*(qfetaeta(i,j,k)+
+ & qfqq(i,j,k)+2.*exp(2.*qf(i,j,k))*
+ & qfphiphi(i,j,k)/sin(qgrd(j))**2+3.*
+ & exp(2.*qf(i,j,k))*qfphi(i,j,k)**2/
+ & sin(qgrd(j))**2)*psisph(i,j,k)
+ enddo
+ enddo
+ enddo
+c
+c Apply boundary conditions to the faces of the cube:
+c
+c i=2:
+ do k = 3,np-2
+ do j = 3,nq-2
+ ae(2,j,k) = ae(2,j,k) + aw(2,j,k)
+ aw(2,j,k) = 0.
+c
+c i=ne-1:
+ ac(ne-1,j,k) = ac(ne-1,j,k)+4.*ae(ne-1,j,k)/(3.+deta)
+ aw(ne-1,j,k) = aw(ne-1,j,k) - ae(ne-1,j,k)/(3.+deta)
+ ae(ne-1,j,k) = 0.
+ enddo
+ enddo
+c
+c j=2:
+ do k = 3,np-2
+ do i = 3,ne-2
+ ac(i,2,k) = ac(i,2,k) + as(i,2,k)
+ as(i,2,k) = 0.
+c
+c j=nq-1:
+ ac(i,nq-1,k) = ac(i,nq-1,k) + an(i,nq-1,k)
+ an(i,nq-1,k) = 0.
+ enddo
+ enddo
+c
+c k=2:
+ do j = 3,nq-2
+ do i = 3,ne-2
+ ac(i,j,2) = ac(i,j,2) + ab(i,j,2)
+ ab(i,j,2) = 0.
+c
+c k=np-1:
+ ac(i,j,np-1) = ac(i,j,np-1) + aq(i,j,np-1)
+ aq(i,j,np-1) = 0.
+ enddo
+ enddo
+c
+c Apply boundary conditions to the edges of the cube:
+c
+c i=2, j=2:
+ do k = 3,np-2
+ ae(2,2,k) = ae(2,2,k) + aw(2,2,k)
+ ac(2,2,k) = ac(2,2,k) + as(2,2,k)
+ aw(2,2,k) = 0.
+ as(2,2,k) = 0.
+c
+c i=ne-1, j=2:
+ aw(ne-1,2,k) = aw(ne-1,2,k) - ae(ne-1,2,k)/(3.+deta)
+ ac(ne-1,2,k) = ac(ne-1,2,k) + as(ne-1,2,k) +
+ & 4.*ae(ne-1,2,k)/(3.+deta)
+ ae(ne-1,2,k) = 0.
+ as(ne-1,2,k) = 0.
+c
+c i=2, j=nq-1:
+ ae(2,nq-1,k) = ae(2,nq-1,k) + aw(2,nq-1,k)
+ ac(2,nq-1,k) = ac(2,nq-1,k) + an(2,nq-1,k)
+ aw(2,nq-1,k) = 0.
+ an(2,nq-1,k) = 0.
+c
+c i=ne-1, j=nq-1:
+ aw(ne-1,nq-1,k) = aw(ne-1,nq-1,k) - ae(ne-1,nq-1,k)/
+ & (3.+deta)
+ ac(ne-1,nq-1,k) = ac(ne-1,nq-1,k) + an(ne-1,nq-1,k) +
+ & 4.*ae(ne-1,nq-1,k)/(3.+deta)
+ ae(ne-1,nq-1,k) = 0.
+ an(ne-1,nq-1,k) = 0.
+ enddo
+c
+c i=2, k=2:
+ do j = 3,nq-2
+ ae(2,j,2) = ae(2,j,2) + aw(2,j,2)
+ ac(2,j,2) = ac(2,j,2) + ab(2,j,2)
+ aw(2,j,2) = 0.
+ ab(2,j,2) = 0.
+c
+c i=ne-1, k=2:
+ aw(ne-1,j,2) = aw(ne-1,j,2) - ae(ne-1,j,2)/(3.+deta)
+ ac(ne-1,j,2) = ac(ne-1,j,2) + ab(ne-1,j,2) +
+ & 4.*ae(ne-1,j,2)/(3.+deta)
+ ae(ne-1,j,2) = 0.
+ ab(ne-1,j,2) = 0.
+c
+c i=2, k=np-1:
+ ae(2,j,np-1) = ae(2,j,np-1) + aw(2,j,np-1)
+ ac(2,j,np-1) = ac(2,j,np-1) + aq(2,j,np-1)
+ aw(2,j,np-1) = 0.
+ aq(2,j,np-1) = 0.
+c
+c i=ne-1, k=np-1:
+ aw(ne-1,j,np-1) = aw(ne-1,j,np-1) - ae(ne-1,j,np-1)/
+ & (3.+deta)
+ ac(ne-1,j,np-1) = ac(ne-1,j,np-1) + aq(ne-1,j,np-1) +
+ & 4.*ae(ne-1,j,np-1)/(3.+deta)
+ ae(ne-1,j,np-1) = 0.
+ aq(ne-1,j,np-1) = 0.
+ enddo
+c
+c j=2, k=2:
+ do i = 3,ne-2
+ ac(i,2,2) = ac(i,2,2) + as(i,2,2) + ab(i,2,2)
+ as(i,2,2) = 0.
+ ab(i,2,2) = 0.
+c
+c j=nq-1, k=2:
+ ac(i,nq-1,2) = ac(i,nq-1,2) + an(i,nq-1,2) +
+ & ab(i,nq-1,2)
+ an(i,nq-1,2) = 0.
+ ab(i,nq-1,2) = 0.
+c
+c j=2, k=np-1:
+ ac(i,2,np-1) = ac(i,2,np-1) + as(i,2,np-1) +
+ & aq(i,2,np-1)
+ as(i,2,np-1) = 0.
+ aq(i,2,np-1) = 0.
+c
+c j=nq-1, k=np-1:
+ ac(i,nq-1,np-1) = ac(i,nq-1,np-1) + an(i,nq-1,np-1) +
+ & aq(i,nq-1,np-1)
+ an(i,nq-1,np-1) = 0.
+ aq(i,nq-1,np-1) = 0.
+ enddo
+c
+c Apply boundary conditions to the corners of the cube:
+c
+c i=2, j=2, k=2:
+ ae(2,2,2) = ae(2,2,2) + aw(2,2,2)
+ ac(2,2,2) = ac(2,2,2) + as(2,2,2) + ab(2,2,2)
+ aw(2,2,2) = 0.
+ as(2,2,2) = 0.
+ ab(2,2,2) = 0.
+c
+c i=ne-1, j=2, k=2:
+ aw(ne-1,2,2) = aw(ne-1,2,2) - ae(ne-1,2,2)/(3.+deta)
+ ac(ne-1,2,2) = ac(ne-1,2,2) + as(ne-1,2,2) + ab(ne-1,2,2) +
+ & 4.*ae(ne-1,2,2)/(3.+deta)
+ ae(ne-1,2,2) = 0.
+ as(ne-1,2,2) = 0.
+ ab(ne-1,2,2) = 0.
+c
+c i=2, j=nq-1, k=2:
+ ae(2,nq-1,2) = ae(2,nq-1,2) + aw(2,nq-1,2)
+ ac(2,nq-1,2) = ac(2,nq-1,2) + an(2,nq-1,2) + ab(2,nq-1,2)
+ aw(2,nq-1,2) = 0.
+ an(2,nq-1,2) = 0.
+ ab(2,nq-1,2) = 0.
+c
+c i=2, j=2, k=np-1:
+ ae(2,2,np-1) = ae(2,2,np-1) + aw(2,2,np-1)
+ ac(2,2,np-1) = ac(2,2,np-1) + as(2,2,np-1) + aq(2,2,np-1)
+ aw(2,2,np-1) = 0.
+ as(2,2,np-1) = 0.
+ aq(2,2,np-1) = 0.
+c
+c i=ne-1, j=nq-1, k=2:
+ aw(ne-1,nq-1,2) = aw(ne-1,nq-1,2) - ae(ne-1,nq-1,2)/(3.+deta)
+ ac(ne-1,nq-1,2) = ac(ne-1,nq-1,2) + an(ne-1,nq-1,2) + ab(ne-1,nq-1,2) +
+ & 4.*ae(ne-1,nq-1,2)/(3.+deta)
+ ae(ne-1,nq-1,2) = 0.
+ an(ne-1,nq-1,2) = 0.
+ ab(ne-1,nq-1,2) = 0.
+c
+c i=ne-1, j=2, k=np-1:
+ aw(ne-1,2,np-1) = aw(ne-1,2,np-1) - ae(ne-1,2,np-1)/(3.+deta)
+ ac(ne-1,2,np-1) = ac(ne-1,2,np-1) + as(ne-1,2,np-1) + aq(ne-1,2,np-1) +
+ & 4.*ae(ne-1,2,np-1)/(3.+deta)
+ ae(ne-1,2,np-1) = 0.
+ as(ne-1,2,np-1) = 0.
+ aq(ne-1,2,np-1) = 0.
+c
+c i=2, j=nq-1, k=np-1:
+ ae(2,nq-1,np-1) = ae(2,nq-1,np-1) + aw(2,nq-1,np-1)
+ ac(2,nq-1,np-1) = ac(2,nq-1,np-1) + an(2,nq-1,np-1) + aq(2,nq-1,np-1)
+ aw(2,nq-1,np-1) = 0.
+ an(2,nq-1,np-1) = 0.
+ aq(2,nq-1,np-1) = 0.
+c
+c i=ne-1, j=nq-1, k=np-1:
+ aw(ne-1,nq-1,np-1) = aw(ne-1,nq-1,np-1) - ae(ne-1,nq-1,np-1)/(3.+deta)
+ ac(ne-1,nq-1,np-1) = ac(ne-1,nq-1,np-1) + an(ne-1,nq-1,np-1) +
+ & aq(ne-1,nq-1,np-1) + 4.*ae(ne-1,nq-1,np-1)/(3.+deta)
+ ae(ne-1,nq-1,np-1) = 0.
+ an(ne-1,nq-1,np-1) = 0.
+ aq(ne-1,nq-1,np-1) = 0.
+c
+c Solve for psi:
+c
+ call bicgst3d(ac,ae,aw,an,as,aq,ab,psiprim,rhs,dbh_eps,rmax,ier,
+ $ ne,nq,np)
+c
+ if (rmax.gt.dbh_eps) then
+ write(*,*) '***WARNING: bicgst3d did not converge.'
+ endif
+ if (ier.eq.-1) then
+ write(*,*) '***WARNING: ier=-1'
+ endif
+ print *,'psiprim = ',maxval(psiprim),' ',minval(psiprim)
+c
+c Now, apply boundary conditions to psiprim:
+c
+ do k = 1,np
+ do j = 1,nq
+ psiprim(1,j,k) = psiprim(3,j,k)
+ psiprim(ne,j,k) = (4.*psiprim(ne-1,j,k)-psiprim(ne-2,j,k))/
+ $ (3.+deta)
+ enddo
+ enddo
+ do k = 1,np
+ do i = 1,ne
+ psiprim(i,1,k) = psiprim(i,2,k)
+ psiprim(i,nq,k) = psiprim(i,nq-1,k)
+ enddo
+ enddo
+ do j = 1,nq
+ do i = 1,ne
+ psiprim(i,j,1) = psiprim(i,j,2)
+ psiprim(i,j,np) = psiprim(i,j,np-1)
+ enddo
+ enddo
+c
+c Here, compute the derivatives of the spherical conformal factor
+c
+c goto 110
+
+ do k = 1,np
+ do j = 1,nq
+ do i = 2,ne-1
+ detapsisph(i,j,k)=0.5*(psiprim(i+1,j,k)-psiprim(i-1,j,k))
+ $ /deta + sinh(0.5*etagrd(i))
+ enddo
+ detapsisph(1,j,k) = -detapsisph(3,j,k)
+ enddo
+ enddo
+c
+ do k = 1,np
+ do j = 2,nq-1
+ do i = 1,ne
+ dqpsisph(i,j,k)=0.5*(psiprim(i,j+1,k)-psiprim(i,j-1,k))/
+ $ dq
+ enddo
+ enddo
+ do i = 1,ne
+ dqpsisph(i,1,k) = -dqpsisph(i,2,k)
+ dqpsisph(i,nq,k) = -dqpsisph(i,nq-1,k)
+ enddo
+ enddo
+c
+ do k = 2,np-1
+ do j = 1,nq
+ do i = 1,ne
+ dphipsisph(i,j,k)=0.5*(psiprim(i,j,k+1)-psiprim(i,j,k-1))
+ $ /dphi
+ enddo
+ enddo
+ enddo
+ do j = 1,nq
+ do i = 1,ne
+ dphipsisph(i,j,1) = -dphipsisph(i,j,2)
+ dphipsisph(i,j,np) = -dphipsisph(i,j,np-1)
+ enddo
+ enddo
+c
+ do k = 1,np
+ do j = 1,nq
+ do i = 2,ne-1
+ detaetapsisph(i,j,k)=(psiprim(i+1,j,k)-2.*psiprim(i,j,k)+
+ & psiprim(i-1,j,k))/deta**2+sqrt(0.25)*
+ & cosh(0.5*etagrd(i))
+ enddo
+ detaetapsisph(1,j,k) = detaetapsisph(3,j,k)
+ enddo
+ enddo
+c
+ do k = 1,np
+ do j = 2,nq-1
+ do i = 1,ne
+ detaqpsisph(i,j,k)=0.5*(detapsisph(i,j+1,k)-detapsisph(i,
+ $ j-1,k))/dq
+ enddo
+ enddo
+ do i = 1,ne
+ detaqpsisph(i,1,k) = -detaqpsisph(i,2,k)
+ detaqpsisph(i,nq,k) = -detaqpsisph(i,nq-1,k)
+ enddo
+ enddo
+c
+ do k = 2,np-1
+ do j = 1,nq
+ do i = 1,ne
+ detaphipsisph(i,j,k)=0.5*(detapsisph(i,j,k+1)-detapsisph(
+ $ i,j,k-1))/dphi
+ enddo
+ enddo
+ enddo
+ do j = 1,nq
+ do i = 1,ne
+ detaphipsisph(i,j,1) = -detaphipsisph(i,j,2)
+ detaphipsisph(i,j,np) = -detaphipsisph(i,j,np-1)
+ enddo
+ enddo
+c
+ do k = 1,np
+ do j = 2,nq-1
+ do i = 1,ne
+ dqqpsisph(i,j,k)=0.5*(dqpsisph(i,j+1,k)-dqpsisph(i,j-1,k))/
+ $ dq
+ enddo
+ enddo
+ do i = 1,ne
+ dqqpsisph(i,1,k) = dqqpsisph(i,2,k)
+ dqqpsisph(i,nq,k) = dqqpsisph(i,nq-1,k)
+ enddo
+ enddo
+c
+ do k = 2,np-1
+ do j = 1,nq
+ do i = 1,ne
+ dqphipsisph(i,j,k)=0.5*(dqpsisph(i,j,k+1)-dqpsisph(i,j,k-1))/
+ $ dphi
+ enddo
+ enddo
+ enddo
+ do j = 1,nq
+ do i = 1,ne
+ dqphipsisph(i,j,1) = -dqphipsisph(i,j,2)
+ dqphipsisph(i,j,np) = -dqphipsisph(i,j,np-1)
+ enddo
+ enddo
+c
+ do k = 2,np-1
+ do j = 1,nq
+ do i = 1,ne
+ dphiphipsisph(i,j,k)=0.5*(dphipsisph(i,j,k+1)-
+ $ dphipsisph(i,j,k-1))/dphi
+ enddo
+ enddo
+ enddo
+ do j = 1,nq
+ do i = 1,ne
+ dphiphipsisph(i,j,1) = dphiphipsisph(i,j,2)
+ dphiphipsisph(i,j,np) = dphiphipsisph(i,j,np-1)
+ enddo
+ enddo
+c
+ do k = 1,np
+ do j = 1,nq
+ psisph(:,j,k)=psiprim(:,j,k)+2.0*cosh(0.5*etagrd)
+ enddo
+ enddo
+
+c
+c Now compute on the Cartesian coordinate.
+c
+c Compute eta,q,phi at the each points of cartesian grid
+ do k = 1,nz
+ do j = 1,ny
+ do i = 1,nx
+ eta(i,j,k) = 0.5d0*dlog(x(i,j,k)**2+y(i,j,k)**2+z(i,j,k)**2)
+ abseta(i,j,k) = abs(eta(i,j,k))
+ if(eta(i,j,k) .lt. 0)then
+ sign_eta(i,j,k) = -1.0d0
+ else
+ sign_eta(i,j,k) = 1.0d0
+ endif
+ q(i,j,k) = datan2(sqrt(x(i,j,k)**2+y(i,j,k)**2),z(i,j,k))
+ phi(i,j,k) = datan2(y(i,j,k),x(i,j,k))
+ enddo
+ enddo
+ enddo
+
+ call CCTK_GetInterpHandle (handle, "simple_local")
+
+ npoints = nx*ny*nz
+
+ call CCTK_Interp (cctkGH,ierror,handle,npoints,3,10,10,
+ $ ne,nq,np,abseta,q,phi,
+ $ CCTK_VARIABLE_REAL,CCTK_VARIABLE_REAL,CCTK_VARIABLE_REAL,
+ $ etagrd(1),qgrd(1),phigrd(1)-pi,deta,dq,dphi,
+ $ psisph,detapsisph,dqpsisph,dphipsisph,detaetapsisph,
+ $ detaqpsisph,detaphipsisph,dqqpsisph,dqphipsisph,dphiphipsisph,
+ $ 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,
+ $ CCTK_VARIABLE_REAL,
+ $ psi3d,detapsi3d,dqpsi3d,dphipsi3d,detaetapsi3d,detaqpsi3d,
+ $ detaphipsi3d,dqqpsi3d,dqphipsi3d,dphiphipsi3d,
+ $ 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,
+ $ CCTK_VARIABLE_REAL)
+
+ do k=1,nz
+ do j=1,ny
+ do i=1,nx
+ psi(i,j,k) = psi3d(i,j,k)*exp(-0.5*eta(i,j,k))
+
+ detapsi3d(i,j,k) = sign_eta(i,j,k)*detapsi3d(i,j,k)
+
+c psix = \partial psi / \partial x / psi
+#include "CactusEinstein/DistortedBHIVP/src/psi_1st_deriv.x"
+
+ detaqpsi3d(i,j,k) = sign_eta(i,j,k)*detaqpsi3d(i,j,k)
+ detaphipsi3d(i,j,k) = sign_eta(i,j,k)*detaphipsi3d(i,j,k)
+
+c psixx = \partial^2\psi / \partial x^2 / psi
+#include "CactusEinstein/DistortedBHIVP/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
+c
+ do k=1,nz
+ do j=1,ny
+ do i=1,nx
+ eta(i,j,k) = 0.5d0*dlog(x(i,j,k)**2+y(i,j,k)**2+z(i,j,k)**2)
+ q(i,j,k) = datan2(sqrt(x(i,j,k)**2+y(i,j,k)**2),z(i,j,k))
+ phi(i,j,k) = datan2(y(i,j,k),x(i,j,k))
+#include "CactusEinstein/DistortedBHIVP/src/gij.x"
+ enddo
+ enddo
+ enddo
+
+c Courvature
+ kxx = 0.0d0
+ kxy = 0.0d0
+ kxz = 0.0d0
+ kyy = 0.0d0
+ kyz = 0.0d0
+ kzz = 0.0d0
+
+ 110 continue
+c Set ADM mass
+ i = ne-15
+ adm = 0.0
+ do k=2,np-1
+ do j=2,nq-1
+ adm=adm+(psisph(i,j,k)-(psisph(i+1,j,k)-psisph(i-1,j,k))/
+ $ deta)*exp(0.5*etagrd(i))
+ enddo
+ enddo
+ adm=adm/(nq-2)/(np-2)
+ print *,'ADM mass: ',adm
+ if (CCTK_Equals(initial_lapse,"schwarz")==1) then
+ write (*,*)"Initial with schwarzschild-like lapse"
+ write (*,*)"using alp = (2.*r - adm)/(2.*r+adm)."
+ alp = (2.*r - adm)/(2.*r+adm)
+ endif
+
+ deallocate(ac,ae,aw,an,as,aq,ab,rhs,
+ $ qf,qfetaeta,qfqq,qfphi,qfphiphi,
+ $ psisph,psiprim,detapsisph,dqpsisph,dphipsisph,detaetapsisph,
+ $ detaqpsisph,detaphipsisph,dqqpsisph,dqphipsisph,
+ $ dphiphipsisph,
+ $ psi3d,detapsi3d,dqpsi3d,
+ $ dphipsi3d,detaetapsi3d,detaqpsi3d,
+ $ detaphipsi3d,dqqpsi3d,dqphipsi3d,
+ $ dphiphipsi3d)
+ deallocate(etagrd,qgrd,phigrd)
+ deallocate(eta,abseta,sign_eta,q,phi)
+
+ return
+ end
diff --git a/src/Stab3d.F b/src/Stab3d.F
new file mode 100644
index 0000000..1aaa7da
--- /dev/null
+++ b/src/Stab3d.F
@@ -0,0 +1,280 @@
+ subroutine bicgst3d(cc,ce,cw,cn,cs,ct,cb,u,rhs,
+ & eps,rmax,ier,im,jm,km)
+c
+c This routine was lifted from stab.f. Minor modifications have
+c been made.
+c
+ implicit none
+c
+ integer,intent(in) :: im,jm,km
+ real*8,intent(inout) :: cc(im,jm,km),cn(im,jm,km),cs(im,jm,km),
+ $ ce(im,jm,km),cw(im,jm,km),ct(im,jm,km),cb(im,jm,km)
+ real*8,intent(out) :: eps
+ real*8,intent(out) :: rmax
+ real*8 :: u(im,jm,km),rhs(im,jm,km)
+c Local variable
+ integer ncyc
+ integer iscale,i,j,k,ier
+c
+c******************************************************************************
+c
+
+ ncyc = (im-2)*(jm-2)*(km-2)
+
+ ier=0
+
+*
+* Determine whether we can diagonally scale the problem to speed
+* convergence. Can only be done if there are no zeros on the main
+* diagonal (ie. central difference coefficient).
+*
+ iscale=1
+ do k = 2,km-1
+ do j = 2,jm-1
+ do i = 2,im-1
+ if (cc(i,j,k).eq.0.) iscale = 0
+ enddo
+ enddo
+ enddo
+
+*
+* Do the diagonal scaling if we can.
+*
+ if (iscale.eq.1) then
+
+
+ do k = 2,km-1
+ do j = 2,jm-1
+ do i = 2,im-1
+ rhs(i,j,k)=rhs(i,j,k)/cc(i,j,k)
+ cb(i,j,k)=cb(i,j,k)/cc(i,j,k)
+ ct(i,j,k)=ct(i,j,k)/cc(i,j,k)
+ cw(i,j,k)=cw(i,j,k)/cc(i,j,k)
+ ce(i,j,k)=ce(i,j,k)/cc(i,j,k)
+ cs(i,j,k)=cs(i,j,k)/cc(i,j,k)
+ cn(i,j,k)=cn(i,j,k)/cc(i,j,k)
+ cc(i,j,k)=1.
+ enddo
+ enddo
+ enddo
+
+ endif
+
+*
+* Now call the bicgstab routine
+*
+
+ call bicgstab (cc,cn,cs,ce,cw,ct,cb,u,rhs,eps,ncyc,rmax,ier,
+ & im,jm,km)
+ if (rmax.gt.eps) then
+ ier = -1
+ print*,'Did not converge'
+ print*,' maximum residual = ',rmax
+ print*,' tolerance = ',eps
+ endif
+
+ return
+ end
+c
+c******************************************************************
+c
+ subroutine bicgstab (cc,cn,cs,ce,cw,ct,cb,x,r,tol,ncyc,rnorm,
+ $ ier,im,jm,km)
+c
+c This routine was lifted from stab.f. Minor modifications have
+c been made.
+c
+ implicit none
+c
+ integer,intent(in) :: ncyc,im,jm,km
+ real*8,intent(out) :: cc(im*jm*km),cn(im*jm*km)
+ real*8,intent(out) :: cs(im*jm*km),ce(im*jm*km)
+ real*8,intent(out) :: cw(im*jm*km),ct(im*jm*km)
+ real*8,intent(out) :: cb(im*jm*km)
+ real*8,intent(out) :: tol,rnorm
+ integer,intent(out) :: ier
+ real*8 x(im*jm*km),r(im*jm*km)
+c Local variables
+ integer :: i,j,k,kk
+ real*8 :: p(im*jm*km),Ap(im*jm*km),w(im*jm*km),As(im*jm*km)
+ real*8 :: omega, chi,chi1,chi2, delta, deltap, pp
+*
+***********************************************************************
+*
+
+ do i = 1,im*jm*km
+ p(i) = 0.
+ Ap(i) = 0.
+ w(i) = 0.
+ As(i) = 0.
+ enddo
+
+ kk = 0
+ 10 call usermv(cc,cn,cs,ce,cw,ct,cb,x,Ap,im,jm,km)
+
+
+ do i = 1,im*jm*km
+ r(i) = r(i)-Ap(i)
+ p(i) = r(i)
+ enddo
+
+c delta = sum(r)
+ delta = 0.
+
+ do i = 1,im*jm*km
+ delta = delta+r(i)
+ enddo
+
+ if (delta.eq.0.) then
+ ier=-1
+ return
+ endif
+
+ call usermv(cc,cn,cs,ce,cw,ct,cb,p,Ap,im,jm,km)
+
+c phi = sum(Ap)
+ pp = 0.
+
+ do i = 1,im*jm*km
+ pp = pp+Ap(i)
+ enddo
+ pp = pp/delta
+
+ if (pp.eq.0.) then
+ ier=-1
+ return
+ endif
+
+c rnorm = sum(r**2)
+ rnorm = 0.
+
+ do i = 1,im*jm*km
+ rnorm = rnorm + r(i)**2
+ enddo
+ rnorm=sqrt(rnorm)
+
+c Test if initial guess is great (residual less than tolerance)
+ if (rnorm .lt. tol) return
+
+ 1 continue
+ kk = kk + 1
+
+ omega = 1./pp
+
+
+ do i = 1,im*jm*km
+ w(i) = r(i) - omega*Ap(i)
+ enddo
+
+ call usermv(cc,cn,cs,ce,cw,ct,cb,w,As,im,jm,km)
+
+c chi1 = sum(As*w)
+ chi1 = 0.
+
+ do i = 1,im*jm*km
+ chi1 = chi1+As(i)*w(i)
+ enddo
+c chi2 = sum(As**2)
+ chi2 = 0.
+
+ do i = 1,im*jm*km
+ chi2 = chi2+As(i)**2
+ enddo
+
+ chi = chi1/chi2
+
+
+ do i = 1,im*jm*km
+ r(i) = w(i) - chi*As(i)
+ x(i) = x(i) + omega*p(i) + chi*w(i)
+ enddo
+
+ deltap = delta
+c delta = sum(r)
+ delta = 0.
+
+ do i = 1,im*jm*km
+ delta = delta+r(i)
+ enddo
+
+ if (delta.eq.0.) then
+ goto 10
+ endif
+
+
+ do i = 1,im*jm*km
+ p(i) = r(i) + (p(i)-chi*Ap(i))*omega*
+ & delta/(deltap*chi)
+ enddo
+
+ call usermv(cc,cn,cs,ce,cw,ct,cb,p,Ap,im,jm,km)
+
+c phi = sum(Ap)
+ pp = 0.
+
+ do i = 1,im*jm*km
+ pp = pp+Ap(i)
+ enddo
+ pp=pp/delta
+
+ if (pp.eq.0.) then
+ goto 10
+ endif
+
+ if (kk.gt.ncyc) then
+ print*,' BI-CGStab solver reached maximum nuber of iterations.'
+ ier=-1
+ return
+ endif
+
+c rnorm = sum(r**2)
+ rnorm = 0.
+
+ do i = 1,im*jm*km
+ rnorm = rnorm+r(i)**2
+ enddo
+ rnorm=sqrt(rnorm)
+
+ if (rnorm .gt. tol) goto 1
+ print*,'-----------'
+ write(*,*) "ncyc=",kk
+ print*,'-----------------------------------------------------------'
+ return
+ end
+c
+c******************************************************************
+c
+ subroutine usermv(cc,cn,cs,ce,cw,ct,cb,x,y,im,jm,km)
+c
+c This routine was lifted from stab.f. Minor modifications have
+c been made.
+c
+c Be careful that the c's are zero on their outer boundary!!
+c
+ implicit none
+c
+ integer,intent(in) :: im,jm,km
+ real*8,intent(out) :: cc(im*jm*km),cn(im*jm*km)
+ real*8,intent(out) :: cs(im*jm*km),ce(im*jm*km)
+ real*8,intent(out) :: cw(im*jm*km),ct(im*jm*km)
+ real*8,intent(out) :: cb(im*jm*km)
+ real*8 :: x(im*jm*km), y(im*jm*km)
+c Local variables
+ integer :: i, j, k
+*
+***********************************************************************
+*
+*
+ do i = im*jm+im+1,im*(jm*km-jm-1)
+ y(i) = cw(i)*x(i-1)
+ & +cc(i)*x(i)
+ & +ce(i)*x(i+1)
+ & +cn(i)*x(i+im)
+ & +cs(i)*x(i-im)
+ & +ct(i)*x(i+im*jm)
+ & +cb(i)*x(i-im*jm)
+ enddo
+c
+ return
+ end
+
diff --git a/src/bhbrill3d.m b/src/bhbrill3d.m
new file mode 100644
index 0000000..ef7bccc
--- /dev/null
+++ b/src/bhbrill3d.m
@@ -0,0 +1,132 @@
+$Path = Union[$Path,{"~/SetTensor"}];
+Needs["SetTensor`"];
+
+Dimension = 3;
+x[1] = eta; x[2] = q; x[3] = phi
+qf[eta_,q_,phi_] := amp (Exp[-(eta-eta0)^2/sigma^2]+Exp[-(eta+eta0)^2/sigma^2]) Sin[q]^n (1+c Cos[phi]^2)
+
+md = {
+{Exp[2 qf[eta,q,phi]],0,0},
+{0,Exp[2 qf[eta,q,phi]],0},
+{0,0,Sin[q]^2}} psisph[eta,q,phi]^4;
+InitializeMetric[md];
+
+Clear[exc];
+DefineTensor[exc];
+SetTensor[exc[la,lb],{{0,0,0},{0,0,0},{0,0,0}}];
+
+tmp = RicciR[la,lb] Metricg[ua,ub]+exc[la,lb] Metricg[ua,ub] exc[lc,ld] Metricg[uc,ud]-
+ exc[la,lb] exc[lc,ld] Metricg[ua,uc] Metricg[ub,ud];
+tmp = RicciToAffine[tmp];
+tmp = EvalMT[tmp];
+tmp = ExpandAll[-Exp[2 qf[eta,q,phi]] psisph[eta,q,phi]^5/8 tmp]
+sav=tmp
+tmp = SubFun[sav,psisph[eta,q,phi],2 Cosh[eta/2]+psisph[eta,q,phi]]
+
+(* Make the stencil... *)
+
+stencil = ExpandAll[tmp /. {
+ D[psisph[eta,q,phi],eta]->(psisph[i+1,j,k]-psisph[i-1,j,k])/(2 deta),
+ D[psisph[eta,q,phi],eta,eta]->(psisph[i+1,j,k]+psisph[i-1,j,k]-2 psisph[i,j,k])/(deta deta),
+ D[psisph[eta,q,phi],q]->(psisph[i,j+1,k]-psisph[i,j-1,k])/(2 dq),
+ D[psisph[eta,q,phi],q,q]->(psisph[i,j+1,k]+psisph[i,j-1,k]-2 psisph[i,j,k])/(dq dq),
+ D[psisph[eta,q,phi],phi]->(psisph[i,j,k+1]-psisph[i,j,k-1])/(2 dphi),
+ D[psisph[eta,q,phi],phi,phi]->(psisph[i,j,k+1]+psisph[i,j,k-1]-2 psisph[i,j,k])/(dphi dphi),
+ psisph[eta,q,phi]->psisph[i,j,k]
+ }];
+
+ac = Coefficient[stencil,psisph[i,j,k]]
+an = Coefficient[stencil,psisph[i+1,j,k]]
+as = Coefficient[stencil,psisph[i-1,j,k]]
+ae = Coefficient[stencil,psisph[i,j,k+1]]
+aw = Coefficient[stencil,psisph[i,j,k-1]]
+aq = Coefficient[stencil,psisph[i,j+1,k]]
+ab = Coefficient[stencil,psisph[i,j-1,k]]
+rhs = -SubFun[tmp,psisph[eta,q,phi],0]
+
+FortranOutputOfDepList = "(i,j,k)";
+$FortranReplace = Union[{
+ "UND"->"_",
+ "(eta,q,phi)"->"(i,j,k)"
+}];
+fd = FortranOpen["bhbrill3d.x"];
+FortranWrite[fd,An[i,j,k],an ];
+FortranWrite[fd,As[i,j,k],as ];
+FortranWrite[fd,Ae[i,j,k],ae ];
+FortranWrite[fd,Aw[i,j,k],aw ];
+FortranWrite[fd,Aq[i,j,k],aq ];
+FortranWrite[fd,Ab[i,j,k],ab ];
+FortranWrite[fd,Ac[i,j,k],ac ];
+FortranWrite[fd,Rhs[i,j,k],rhs ];
+FortranClose[fd];
+
+(* Next part, write out conformal g's and d's *)
+
+
+xv = Exp[eta] Sin[q] Cos[phi];
+yv = Exp[eta] Sin[q] Sin[phi];
+zv = Exp[eta] Cos[q];
+
+mc = Table[ D[ {xv,yv,zv}[[i]], {eta,q,phi}[[j]] ],{i,1,3},{j,1,3}];
+mci = Simplify[Inverse[mc]];
+
+Clear[mct];
+DefineTensor[mct,{{1,2},1}];
+Iter[mct[ua,lb],
+ mct[ua,lb]=mc[[ua,-lb]];
+];
+
+Clear[mcti];
+DefineTensor[mcti,{{1,2},1}];
+Iter[mcti[ua,lb],
+ mcti[ua,lb]=mci[[ua,-lb]];
+];
+
+gijtmp = Exp[2 eta]/psisph[eta,q,phi]^4 Metricg[lc,ld] mcti[uc,la] mcti[ud,lb]
+
+Clear[i2];
+DefineTensor[i2,{{2,1},1}];
+
+fd = FortranOpen["gij.x"];
+Iter[i2[ua,ub],
+ v1 = {x,y,z}[[ua]];
+ v2 = {x,y,z}[[ub]];
+ metv = ToExpression["g"<>ToString[v1]<>ToString[v2]<>"[i,j,k]"];
+ gg[v1,v2]=Simplify[EvalMT[gijtmp,{la->-ua,lb->-ub}]];
+ FortranWrite[fd,metv,gg[v1,v2]];
+ For[ii=1,ii<=3,ii++,
+ v3 = {x,y,z}[[ii]];
+ dmetv = ToExpression["d"<>ToString[v3]<>ToString[metv]];
+ res = OD[gg[v1,v2],lc] mcti[uc,ld]/2;
+ res = EvalMT[res,ld-> -ii];
+ res = Simplify[res];
+ FortranWrite[fd,dmetv,res];
+ ];
+];
+FortranClose[fd];
+
+$FortranReplace = {
+ "UND"->"_",
+ "(eta,q,phi)"->""
+};
+
+fd = FortranOpen["psi_1st_deriv.x"];
+For[ii=1,ii<=3,ii++,
+ v1 = {x,y,z}[[ii]];
+ psv =ToExpression["psi"<>ToString[v1]<>"[i,j,k]"];
+ rhs = CD[Exp[-eta/2] psi3d[eta,q,phi],lc] mcti[uc,la];
+ rhs = EvalMT[rhs,{la->-ii}]/(Exp[-eta/2] psi3d[eta,q,phi]);
+ FortranWrite[fd,psv,rhs];
+];
+FortranClose[fd];
+
+fd = FortranOpen["psi_2nd_deriv.x"];
+Iter[i2[ua,ub],
+ v1 = {x,y,z}[[ua]];
+ v2 = {x,y,z}[[ub]];
+ psv = ToExpression["psi"<>ToString[v1]<>ToString[v2]<>"[i,j,k]"];
+ rhs = OD[OD[Exp[-eta/2] psi3d[eta,q,phi],lc] mcti[uc,la],ld] mcti[ud,lb];
+ rhs = EvalMT[rhs,{la->-ua,lb->-ub}]/(Exp[-eta/2] psi3d[eta,q,phi]);
+ FortranWrite[fd,psv,rhs];
+];
+FortranClose[fd];
diff --git a/src/bhbrill3d.x b/src/bhbrill3d.x
new file mode 100644
index 0000000..01a731b
--- /dev/null
+++ b/src/bhbrill3d.x
@@ -0,0 +1,126 @@
+ o1 = deta**2
+ o2 = 1/o1
+ o3 = dq**2
+ o4 = 1/o3
+ o5 = 1/dq
+ o6 = tan(q)
+ o7 = 1/o6
+ o8 = dphi**2
+ o9 = 1/o8
+ o10 = eta**2
+ o11 = sigma**2
+ o12 = 1/o11
+ o13 = -(o10*o12)
+ o14 = -2.00000000000000d0*eta*eta0*o12
+ o15 = eta0**2
+ o16 = -(o12*o15)
+ o17 = o13 + o14 + o16
+ o18 = exp(o17)
+ o19 = sin(q)
+ o20 = o19**n
+ o21 = 2.00000000000000d0*amp*o18*o20
+ o22 = cos(phi)
+ o23 = o22**2
+ o24 = 2.00000000000000d0*amp*c*o18*o20*o23
+ o25 = 2.00000000000000d0*eta*eta0*o12
+ o26 = o13 + o16 + o25
+ o27 = exp(o26)
+ o28 = 2.00000000000000d0*amp*o20*o27
+ o29 = 2.00000000000000d0*amp*c*o20*o23*o27
+ o30 = o21 + o24 + o28 + o29
+ o31 = exp(o30)
+ o32 = o19**2
+ o33 = 1/o32
+ o34 = 1.00000000000000d0*o31*o33*o9
+ o35 = 1/dphi
+ o36 = o13 + o14 + o16 + o21 + o24 + o28 + o29
+ o37 = exp(o36)
+ o38 = sin(phi)
+ o39 = -2.00000000000000d0 + n
+ o40 = o19**o39
+ o41 = o13 + o16 + o21 + o24 + o25 + o28 + o29
+ o42 = exp(o41)
+ o43 = cos(q)
+ o44 = o43**2
+ o45 = n**2
+ o46 = o38**2
+ o47 = o11**2
+ o48 = 1/o47
+ o49 = amp**2
+ o50 = c**2
+ o51 = -2.00000000000000d0*o10*o12
+ o52 = -2.00000000000000d0*o12*o15
+ o53 = o21 + o24 + o28 + o29 + o51 + o52
+ o54 = exp(o53)
+ o55 = 2.00000000000000d0*n
+ o56 = -2.00000000000000d0 + o55
+ o57 = o19**o56
+ o58 = -4.0000000000000d0*eta*eta0*o12
+ o59 = o21 + o24 + o28 + o29 + o51 + o52 + o58
+ o60 = exp(o59)
+ o61 = 4.0000000000000d0*eta*eta0*o12
+ o62 = o21 + o24 + o28 + o29 + o51 + o52 + o61
+ o63 = exp(o62)
+ o64 = o6**2
+ o65 = 1/o64
+ o66 = 5.0000000000000d-1*eta
+ o67 = cosh(o66)
+ Ae(i,j,k) = o2
+ Aw(i,j,k) = o2
+ An(i,j,k) = o4 + 5.0000000000000d-1*o5*o7
+ As(i,j,k) = o4 - 5.0000000000000d-1*o5*o7
+ Ab(i,j,k) = o34 - 2.00000000000000d0*amp*c*o22*o35*o37*o38*o40 -
+ & 2.00000000000000d0*amp*c*o22*o35*o38*o40*o42
+ Aq(i,j,k) = o34 + 2.00000000000000d0*amp*c*o22*o35*o37*o38*o40 +
+ & 2.00000000000000d0*amp*c*o22*o35*o38*o40*o42
+ Ac(i,j,k) = -1.25000000000000d-1 - 2.00000000000000d0*o2 - 2.500
+ & 00000000000d-1*amp*n*o18*o20 - 5.0000000000000d-1*amp*o12*o18*o2
+ & 0 - 2.50000000000000d-1*amp*c*n*o18*o20*o23 - 5.0000000000000d-1
+ & *amp*c*o12*o18*o20*o23 - 2.50000000000000d-1*amp*n*o20*o27 - 5.0
+ & 000000000000d-1*amp*o12*o20*o27 - 2.50000000000000d-1*amp*c*n*o2
+ & 0*o23*o27 - 5.0000000000000d-1*amp*c*o12*o20*o23*o27 - 1.2500000
+ & 0000000d-1*o33 - 2.00000000000000d0*o4 - amp*c*o23*o37*o40 - amp
+ & *c*o23*o40*o42 - 2.50000000000000d-1*amp*n*o18*o40*o44 - 2.50000
+ & 000000000d-1*amp*c*n*o18*o23*o40*o44 - 2.50000000000000d-1*amp*n
+ & *o27*o40*o44 - 2.50000000000000d-1*amp*c*n*o23*o27*o40*o44 + 2.5
+ & 0000000000000d-1*amp*o18*o40*o44*o45 + 2.50000000000000d-1*amp*c
+ & *o18*o23*o40*o44*o45 + 2.50000000000000d-1*amp*o27*o40*o44*o45 +
+ & 2.50000000000000d-1*amp*c*o23*o27*o40*o44*o45 + amp*c*o37*o40*o
+ & 46 + amp*c*o40*o42*o46 + 2.00000000000000d0*amp*eta*eta0*o18*o20
+ & *o48 + amp*o10*o18*o20*o48 + amp*o15*o18*o20*o48 + 2.00000000000
+ & 000d0*amp*c*eta*eta0*o18*o20*o23*o48 + amp*c*o10*o18*o20*o23*o48
+ & + amp*c*o15*o18*o20*o23*o48 - 2.00000000000000d0*amp*eta*eta0*o
+ & 20*o27*o48 + amp*o10*o20*o27*o48 + amp*o15*o20*o27*o48 - 2.00000
+ & 000000000d0*amp*c*eta*eta0*o20*o23*o27*o48 + amp*c*o10*o20*o23*o
+ & 27*o48 + amp*c*o15*o20*o23*o27*o48 + 6.0000000000000d0*o23*o46*o
+ & 49*o50*o54*o57 + 3.00000000000000d0*o23*o46*o49*o50*o57*o60 + 3.
+ & 00000000000000d0*o23*o46*o49*o50*o57*o63 + 1.25000000000000d-1*o
+ & 65 - 2.00000000000000d0*o31*o33*o9
+ t1 = -2.50000000000000d-1*o67 + 5.0000000000000d-1*amp*n*o18*o20
+ & *o67 + amp*o12*o18*o20*o67 + 5.0000000000000d-1*amp*c*n*o18*o20*
+ & o23*o67 + amp*c*o12*o18*o20*o23*o67 + 5.0000000000000d-1*amp*n*o
+ & 20*o27*o67 + amp*o12*o20*o27*o67 + 5.0000000000000d-1*amp*c*n*o2
+ & 0*o23*o27*o67 + amp*c*o12*o20*o23*o27*o67 + 2.50000000000000d-1*
+ & o33*o67 + 2.00000000000000d0*amp*c*o23*o37*o40*o67 + 2.000000000
+ & 00000d0*amp*c*o23*o40*o42*o67 + 5.0000000000000d-1*amp*n*o18*o40
+ & *o44*o67 + 5.0000000000000d-1*amp*c*n*o18*o23*o40*o44*o67 + 5.00
+ & 00000000000d-1*amp*n*o27*o40*o44*o67 + 5.0000000000000d-1*amp*c*
+ & n*o23*o27*o40*o44*o67 - 5.0000000000000d-1*amp*o18*o40*o44*o45*o
+ & 67 - 5.0000000000000d-1*amp*c*o18*o23*o40*o44*o45*o67 - 5.000000
+ & 0000000d-1*amp*o27*o40*o44*o45*o67
+ t2 = -5.0000000000000d-1*amp*c*o23*o27*o40*o44*o45*o67 - 2.00000
+ & 000000000d0*amp*c*o37*o40*o46*o67 - 2.00000000000000d0*amp*c*o40
+ & *o42*o46*o67 - 4.0000000000000d0*amp*eta*eta0*o18*o20*o48*o67 -
+ & 2.00000000000000d0*amp*o10*o18*o20*o48*o67 - 2.00000000000000d0*
+ & amp*o15*o18*o20*o48*o67 - 4.0000000000000d0*amp*c*eta*eta0*o18*o
+ & 20*o23*o48*o67 - 2.00000000000000d0*amp*c*o10*o18*o20*o23*o48*o6
+ & 7 - 2.00000000000000d0*amp*c*o15*o18*o20*o23*o48*o67 + 4.0000000
+ & 000000d0*amp*eta*eta0*o20*o27*o48*o67 - 2.00000000000000d0*amp*o
+ & 10*o20*o27*o48*o67 - 2.00000000000000d0*amp*o15*o20*o27*o48*o67
+ & + 4.0000000000000d0*amp*c*eta*eta0*o20*o23*o27*o48*o67 - 2.00000
+ & 000000000d0*amp*c*o10*o20*o23*o27*o48*o67 - 2.00000000000000d0*a
+ & mp*c*o15*o20*o23*o27*o48*o67 - 1.20000000000000d1*o23*o46*o49*o5
+ & 0*o54*o57*o67 - 6.0000000000000d0*o23*o46*o49*o50*o57*o60*o67 -
+ & 6.0000000000000d0*o23*o46*o49*o50*o57*o63*o67 - 2.50000000000000
+ & d-1*o65*o67
+ Rhs(i,j,k) = t1 + t2
diff --git a/src/gij.x b/src/gij.x
new file mode 100644
index 0000000..5a9d688
--- /dev/null
+++ b/src/gij.x
@@ -0,0 +1,117 @@
+ o1 = 2.00000000000000d0*phi(i,j,k)
+ o2 = cos(o1)
+ o3 = cos(phi(i,j,k))
+ o4 = o3**2
+ o5 = c*o4
+ o6 = 1.00000000000000d0 + o5
+ o7 = -eta0
+ o8 = eta(i,j,k) + o7
+ o9 = o8**2
+ o10 = sigma**2
+ o11 = 1/o10
+ o12 = -(o11*o9)
+ o13 = exp(o12)
+ o14 = eta(i,j,k) + eta0
+ o15 = o14**2
+ o16 = -(o11*o15)
+ o17 = exp(o16)
+ o18 = o13 + o17
+ o19 = sin(q(i,j,k))
+ o20 = o19**n
+ o21 = 2.00000000000000d0*amp*o18*o20*o6
+ o22 = exp(o21)
+ o23 = -1.00000000000000d0 + o22
+ o24 = -eta(i,j,k)
+ o25 = exp(o24)
+ o26 = sin(phi(i,j,k))
+ o27 = o26**2
+ o28 = 1/o19
+ o29 = sin(o1)
+ o30 = cos(q(i,j,k))
+ o31 = o30**2
+ o32 = -1.00000000000000d0 + n
+ o33 = o19**o32
+ o34 = amp*n*o18*o22*o3*o31*o33*o6
+ o35 = 2.00000000000000d0*amp*c*o18*o22*o27*o3*o33
+ o36 = 4.0000000000000d0*eta(i,j,k)*eta0*o11
+ o37 = exp(o36)
+ o38 = eta(i,j,k)*o37
+ o39 = -(eta0*o37)
+ o40 = eta(i,j,k) + eta0 + o38 + o39
+ o41 = o16 + o21
+ o42 = exp(o41)
+ o43 = 1.00000000000000d0 + n
+ o44 = o19**o43
+ o45 = -2.00000000000000d0*amp*o11*o3*o40*o42*o44*o6
+ o46 = -2.00000000000000d0*amp*c*o18*o22*o26*o33*o4
+ o47 = amp*n*o18*o22*o26*o31*o33*o6
+ o48 = -2.00000000000000d0*amp*o11*o26*o40*o42*o44*o6
+ o49 = c*o2
+ o50 = 2.00000000000000d0 + c + o49
+ o51 = -1.00000000000000d0 + o37
+ o52 = -2.00000000000000d0*eta0*o51
+ o53 = 1.00000000000000d0 + o37
+ o54 = 2.00000000000000d0*eta(i,j,k)*o53
+ o55 = n*o10*o53
+ o56 = o52 + o54 + o55
+ o57 = o16 + o21 + o24
+ o58 = exp(o57)
+ o59 = o11*o9
+ o60 = exp(o59)
+ o61 = o11*o15
+ o62 = exp(o61)
+ o63 = o60 + o62
+ o64 = o19**2
+ gxx(i,j,k) = 5.0000000000000d-1*(1.00000000000000d0 + o22 + o2*o
+ & 23)
+c dxgxx(i,j,k) = 5.0000000000000d-1*o25*(1.00000000000000d0*o22*o2
+c & 6*o28*o29 - 2.00000000000000d0*o27*o28*o3 + 2.00000000000000d0*a
+c & mp*c*o18*o2*o22*o27*o3*o33 + o34 + o35 + o45 + amp*n*o18*o2*o22*
+c & o3*o31*o33*o6 - 2.00000000000000d0*amp*o11*o2*o3*o40*o42*o44*o6)
+c dygxx(i,j,k) = 5.0000000000000d-1*o25*(1.00000000000000d0*o28*o2
+c & 9*o3 - 2.00000000000000d0*o22*o26*o28*o4 - 2.00000000000000d0*am
+c & p*c*o18*o2*o22*o26*o33*o4 + o46 + o47 + o48 + amp*n*o18*o2*o22*o
+c & 26*o31*o33*o6 - 2.00000000000000d0*amp*o11*o2*o26*o40*o42*o44*o6
+c & )
+c dzgxx(i,j,k) = -5.0000000000000d-1*amp*o11*o20*o30*o4*o50*o56*o5
+c & 8
+ gxy(i,j,k) = 5.0000000000000d-1*o23*o29
+c dxgxy(i,j,k) = 5.0000000000000d-1*o25*(1.00000000000000d0*o2*o26
+c & *o28 - o2*o22*o26*o28 + 2.00000000000000d0*amp*c*o18*o22*o27*o29
+c & *o3*o33 + amp*n*o18*o22*o29*o3*o31*o33*o6 - 2.00000000000000d0*a
+c & mp*o11*o29*o3*o40*o42*o44*o6)
+c dygxy(i,j,k) = 5.0000000000000d-1*o25*(-(o2*o28*o3) + 1.00000000
+c & 000000d0*o2*o22*o28*o3 - 2.00000000000000d0*amp*c*o18*o22*o26*o2
+c & 9*o33*o4 + amp*n*o18*o22*o26*o29*o31*o33*o6 - 2.00000000000000d0
+c & *amp*o11*o26*o29*o40*o42*o44*o6)
+c dzgxy(i,j,k) = -5.0000000000000d-1*amp*o11*o20*o29*o30*o56*o58*o
+c & 6
+ gxz(i,j,k) = 0
+c dxgxz(i,j,k) = 0
+c dygxz(i,j,k) = 0
+c dzgxz(i,j,k) = 0
+ gyy(i,j,k) = 5.0000000000000d-1*(1.00000000000000d0 + o22 - o2*o
+ & 23)
+c dxgyy(i,j,k) = 5.0000000000000d-1*o25*(1.00000000000000d0*o26*o2
+c & 8*o29 - 2.00000000000000d0*o22*o27*o28*o3 - 2.00000000000000d0*a
+c & mp*c*o18*o2*o22*o27*o3*o33 + o34 + o35 + o45 - amp*n*o18*o2*o22*
+c & o3*o31*o33*o6 + 2.00000000000000d0*amp*o11*o2*o3*o40*o42*o44*o6)
+c dygyy(i,j,k) = 5.0000000000000d-1*o25*(1.00000000000000d0*o22*o2
+c & 8*o29*o3 - 2.00000000000000d0*o26*o28*o4 + 2.00000000000000d0*am
+c & p*c*o18*o2*o22*o26*o33*o4 + o46 + o47 + o48 - amp*n*o18*o2*o22*o
+c & 26*o31*o33*o6 + 2.00000000000000d0*amp*o11*o2*o26*o40*o42*o44*o6
+c & )
+c dzgyy(i,j,k) = -5.0000000000000d-1*amp*o11*o20*o27*o30*o50*o56*o
+c & 58
+ gyz(i,j,k) = 0
+c dxgyz(i,j,k) = 0
+c dygyz(i,j,k) = 0
+c dzgyz(i,j,k) = 0
+ gzz(i,j,k) = o22
+c dxgzz(i,j,k) = amp*o11*o3*o33*(n*o10*o31*o6*o63 + 2.000000000000
+c & 00d0*(c*o10*o27*o63 - o40*o6*o60*o64))*exp(o12 + o16 + o21 + o24
+c & )
+c dygzz(i,j,k) = amp*o26*o33*(-2.00000000000000d0*c*o18*o4 + n*o18
+c & *o31*o6 - 2.00000000000000d0*o11*o17*o40*o6*o64)*exp(o21 + o24)
+c dzgzz(i,j,k) = -(amp*o11*o20*o30*o56*o58*o6)
+
diff --git a/src/interp3.F b/src/interp3.F
new file mode 100644
index 0000000..0bbdf90
--- /dev/null
+++ b/src/interp3.F
@@ -0,0 +1,264 @@
+c /*@@
+c @routine interp3
+c @date Fri Feb 14 08:46:53 1997
+c @author Ryoji Takahashi
+c @desc
+c Interpolates from 3D data var with coordinates x, y, and z and
+c sizes nx , ny, and nz onto 1D data out with position outx, outy
+c ,and outz nout points.
+c <p>
+c This has linear, quadratic and cubic interpolators in it.
+c Or will one day very soon.
+c @enddesc
+c @calls
+c @calledby numerical_nonaxi
+c@@*/
+
+ subroutine interp3(dat,x,y,z,nx,ny,nz,out,outx,outy,outz,nout
+ $ ,order)
+ implicit none
+ integer nx,ny,nz,nout
+ real*8 dat(nx,ny,nz), x(nx), y(ny), z(nz)
+ real*8 out(nout),outx(nout),outy(nout),outz(nout)
+ integer order
+c Interpolation goes from ibelow to ibelow+[1,2,3] depending on order
+ integer i,j,k,l,m,ibelow,jbelow,kbelow,pt
+ real*8 xsym,ysym,zsym,findx,findy,findz,frac
+ real*8 ydir(order+1)
+ real*8 zdir(order+1)
+ real*8 f(4), fp(4,4), fl(4)
+ real*8 dbh_linear, dbh_quad, dbh_cubic
+ real*8 dx, dy, dz, PI
+ integer twobhjsad
+
+ pi = 4.0d0*atan(1.0d0)
+
+ dx = x(2) - x(1)
+ dy = y(2) - y(1)
+ dz = z(2) - z(1)
+
+c Loop over all out points
+ do pt=1,nout
+ zsym = 1.0D0
+ ysym = 1.0D0
+ xsym = 1.0D0
+c Check bounds
+ findx = outx(pt)
+ if (findx .lt. x(1)) then
+ write (*,*) "Below inner bound at ",pt,outx(pt)
+ STOP
+ endif
+ if (findx .gt. x(nx)) then
+ write (*,*) "Above x bounds at ",pt,outx(pt),x(nx)
+ STOP
+ endif
+ findy = outy(pt)
+ if (findy .lt. y(1)) then
+ write (*,*) "Below y inner bound at ",pt,outy(pt)
+ STOP
+ endif
+ if (findy .gt. y(ny)) then
+ write (*,*) "Below y inner bound at ",pt,outy(pt)
+ STOP
+ endif
+ findz = outz(pt)+pi
+ if (findz .lt. z(1)) then
+ write (*,*) "Below z inner bound at ",pt,outz(pt)
+ STOP
+ endif
+ if (findz .gt. z(nz)) then
+ write (*,*) "Below z inner bound at ",pt,outz(pt)
+ STOP
+ endif
+
+c Locate ourselves in i,j space
+c do i=1,nx
+c if (x(i) .lt. findx) then
+c ibelow = i
+c endif
+c enddo
+c Assume a regular grid
+ ibelow = (findx-x(1))/dx+1
+
+ if (order .eq. 3 .and. ibelow .gt. 1) then
+ ibelow = ibelow - 1
+ endif
+ if (ibelow + order .gt. nx) then
+ ibelow = nx - order
+ endif
+
+c do i=1,ny
+c if (y(i) .lt. findy) then
+c jbelow = i
+c endif
+c enddo
+c Assume a regular grid
+ jbelow = (findy-y(1))/dy+1
+ if (order .eq. 3 .and. jbelow .gt. 1) then
+ jbelow = jbelow - 1
+ endif
+ if (jbelow + order .gt. ny) then
+ jbelow = ny - order
+ endif
+
+c do i=1,nz
+c if (z(i) .lt. findz) then
+c kbelow = i
+c endif
+c enddo
+c Assume a regular grid
+ kbelow = (findz-z(1))/dz+1
+ if (order .eq. 3 .and. kbelow .gt. 1) then
+ kbelow = kbelow - 1
+ endif
+ if (kbelow + order .gt. nz) then
+ kbelow = nz - order
+ endif
+
+c write (*,*) "PT :",findx,findy
+c write (*,*) "SYM:",sym
+c write (*,*) "BOUND X ",ibelow,x(ibelow),x(ibelow+1)
+c write (*,*) "BOUND Y ",jbelow,y(jbelow),y(jbelow+1)
+c write (*,*) "BOUND z ",kbelow,z(kbelow),z(kbelow+1)
+
+c So do the interpolation
+ if (order .eq. 1) then
+c Interp in the x direction
+! frac = (findx-x(ibelow))/(x(ibelow+1)-x(ibelow))
+ do l = 1,2
+ do m = 1,2
+ f(1) = dat(ibelow,jbelow+l-1,kbelow+m-1)
+ f(2) = dat(ibelow+1,jbelow+l-1,kbelow+m-1)
+ fp(l,m) = dbh_linear(f,x(ibelow),dx,findx)
+ enddo
+ enddo
+c Now take our 2x2 plane and interp to the center of both
+c in the y direction
+! frac = (findy-y(jbelow))/(y(jbelow+1)-y(jbelow))
+ do m = 1,2
+ f(1) = fp(1,m)
+ f(2) = fp(2,m)
+ fl(m) = dbh_linear(f,y(jbelow),dy,findy)
+ enddo
+c And finally, interp in the z direction
+ out(pt) = xsym * ysym * zsym *
+ $ dbh_linear(fl,z(kbelow),dz,findz)
+ else if (order .eq. 2) then
+c Load up for calls to poly2inter
+ do l=1,3
+ do m=1,3
+ f(1) = dat(ibelow,jbelow+l-1,kbelow+m-1)
+ f(2) = dat(ibelow+1,jbelow+l-1,kbelow+m-1)
+ f(3) = dat(ibelow+2,jbelow+l-1,kbelow+m-1)
+ fp(l,m) = dbh_quad(f,x(ibelow),dx,findx)
+ enddo
+ enddo
+c Now take our 2x2 plane and interp to the center of both
+c in the y direction
+ do m=1,3
+ f(1) = fp(1,m)
+ f(2) = fp(2,m)
+ f(3) = fp(3,m)
+ fl(m) = dbh_quad(f,y(jbelow),dy,findy)
+ enddo
+c And finally, interp in the z direction
+ out(pt) = xsym * ysym * zsym *
+ $ dbh_quad(fl,z(kbelow),dz,findz)
+ else if (order .eq. 3) then
+c Load up for calls to cubic
+ do l=1,4
+ do m=1,4
+ f(1) = dat(ibelow,jbelow+l-1,kbelow+m-1)
+ f(2) = dat(ibelow+1,jbelow+l-1,kbelow+m-1)
+ f(3) = dat(ibelow+2,jbelow+l-1,kbelow+m-1)
+ f(4) = dat(ibelow+3,jbelow+l-1,kbelow+m-1)
+ fp(l,m) = dbh_cubic(f,x(ibelow),dx,findx)
+ enddo
+ enddo
+c Now take our 2x2 plane and interp to the center of both
+c in the y direction
+ do m=1,4
+ f(1) = fp(1,m)
+ f(2) = fp(2,m)
+ f(3) = fp(3,m)
+ f(4) = fp(4,m)
+ fl(m) = dbh_cubic(f,y(jbelow),dy,findz)
+ enddo
+c And finally, interp in the z direction
+ out(pt) = xsym * ysym * zsym *
+ $ dbh_cubic(fl,z(kbelow),dz,findz)
+ else
+ write (*,*) "ORDER set wrong in interp3d",order
+ stop
+ endif
+
+ enddo
+ return
+ end
+
+ real*8 function dbh_linear(f, x0, dx, findx)
+ implicit none
+ real*8 f(2),x0,dx,findx
+ real*8 frac
+
+ frac = (findx-x0)/dx
+ dbh_linear = (frac)*f(2) + (1.0-frac)*f(1)
+
+ return
+ end
+
+ real*8 function dbh_quad(f, x0, dx, findx)
+ implicit none
+ real*8 f(3),x0, dx, findx, dbh_quad
+ real*8 f0,f1,f2
+ real*8 a,b,c, dx2, x02, o2dx2
+c Mathematica tells us
+c - List(List(Rule(c,(2*dx**2*f0 + 3*dx*f0*x0 - 4*dx*f1*x0
+c - + dx*f2*x0 + f0*x0**2 - 2*f1*x0**2 + f2*x0**2)
+c - /(2*dx**2)), Rule(b,(-3*dx*f0 + 4*dx*f1 - dx*f2 - 2*f0*x0
+c - + 4*f1*x0 - 2*f2*x0)/(2*dx**2)),Rule(a,(f0 - 2*f1 +
+c - f2)/(2*dx**2))))
+
+ f0 = f(1)
+ f1 = f(2)
+ f2 = f(3)
+ dx2 = dx**2
+ x02 = x0**2
+ o2dx2 = 1.0D0/(2.0D0*dx2)
+
+ c = (2.0D0*dx2*f0 + dx*x0*(3.0D0*f0 - 4.0D0*f1 + f2) +
+ $ x02*(f0 - 2.0D0*f1 + f2))*o2dx2
+ b = (dx * (-3.0D0*f0 + 4.0D0*f1 - f2) + x0 * (- 2.0D0*f0 +
+ $ 4.0D0*f1 - 2.0D0*f2))*o2dx2
+
+ a = (f0 - 2.0D0*f1 + f2)*o2dx2
+
+ dbh_quad = a*findx**2 + b*findx + c
+
+ end
+
+
+ real*8 function dbh_cubic(f, x0, dx, findx)
+ implicit none
+ real*8 a,b,c,d,dbh_cubic
+ real*8 f(4),x0,dx,findx
+
+ a = -(f(1)-3.0*f(2)+3.0*f(3)-f(4)) / (6.0*(dx**3))
+
+ b = (f(1)-2.0*f(2)+f(3))/(2.0*(dx**2)) +
+ $ (f(1)-3.0*f(2)+3.0*f(3)-f(4))*(dx+x0)/(2.0*(dx**3))
+
+ c = ((dx**2)*(-11.0*f(1) + 18.0*f(2) - 9.0*f(3) + 2.0*f(4)) +
+ $ dx*x0* (-12.0*f(1) + 30.0*f(2) - 24.0*f(3) + 6.0*f(4)) +
+ $ (x0**2)*( -3.0*f(1) + 9.0*f(2) - 9.0*f(3) + 3.0*f(4))) /
+ $ (6.0*(dx**3))
+
+ d = ((dx**3)* ( 6.0*f(1) ) +
+ $ (dx**2)*x0*( 11.0*f(1) - 18.0*f(2) + 9.0*f(3) - 2.0*f(4)) +
+ $ (x0**2)*dx*( 6.0*f(1) - 15.0*f(2) + 12.0*f(3) - 3.0*f(4)) +
+ $ (x0**3)* ( 1.0*f(1) - 3.0*f(2) + 3.0*f(3) - 1.0*f(4)))/
+ $ (6.0*(dx**3))
+
+ dbh_cubic = ((a*findx + b)*findx + c)*findx + d
+ return
+ end
diff --git a/src/make.code.defn b/src/make.code.defn
new file mode 100644
index 0000000..e34babf
--- /dev/null
+++ b/src/make.code.defn
@@ -0,0 +1,9 @@
+# Main make.code.defn file for thorn DistortedBH
+# $Header$
+
+# Source files in this directory
+SRCS = DistortedBHIVP.F Stab3d.F
+
+# Subdirectories containing source files
+SUBDIRS =
+
diff --git a/src/psi_1st_deriv.x b/src/psi_1st_deriv.x
new file mode 100644
index 0000000..d0a0be2
--- /dev/null
+++ b/src/psi_1st_deriv.x
@@ -0,0 +1,20 @@
+ o1 = 5.0000000000000d-1*eta(i,j,k)
+ o2 = exp(o1)
+ o3 = psi3d(i,j,k)
+ o4 = 1/o3
+ o5 = cos(phi(i,j,k))
+ o6 = cos(q(i,j,k))
+ o7 = dqpsi3d(i,j,k)
+ o8 = -1.50000000000000d0*eta(i,j,k)
+ o9 = exp(o8)
+ o10 = dphipsi3d(i,j,k)
+ o11 = sin(phi(i,j,k))
+ o12 = sin(q(i,j,k))
+ o13 = 1/o12
+ o14 = detapsi3d(i,j,k)
+ psix(i,j,k) = o2*o4*(-(o10*o11*o13*o9) + o12*o14*o5*o9 - 5.00000
+ & 00000000d-1*o12*o3*o5*o9 + o5*o6*o7*o9)
+ psiy(i,j,k) = o2*o4*(o11*o12*o14*o9 - 5.0000000000000d-1*o11*o12
+ & *o3*o9 + 1.00000000000000d0*o10*o13*o5*o9 + o11*o6*o7*o9)
+ psiz(i,j,k) = o2*o4*(o14*o6*o9 - 5.0000000000000d-1*o3*o6*o9 - o
+ & 12*o7*o9)
diff --git a/src/psi_2nd_deriv.x b/src/psi_2nd_deriv.x
new file mode 100644
index 0000000..5890b02
--- /dev/null
+++ b/src/psi_2nd_deriv.x
@@ -0,0 +1,71 @@
+ o1 = 5.0000000000000d-1*eta(i,j,k)
+ o2 = exp(o1)
+ o3 = psi3d(i,j,k)
+ o4 = 1/o3
+ o5 = cos(phi(i,j,k))
+ o6 = o5**2
+ o7 = cos(q(i,j,k))
+ o8 = o7**2
+ o9 = detapsi3d(i,j,k)
+ o10 = -2.50000000000000d0*eta(i,j,k)
+ o11 = exp(o10)
+ o12 = dqqpsi3d(i,j,k)
+ o13 = detaphipsi3d(i,j,k)
+ o14 = sin(phi(i,j,k))
+ o15 = dphipsi3d(i,j,k)
+ o16 = o14**2
+ o17 = sin(q(i,j,k))
+ o18 = o17**2
+ o19 = 1/o18
+ o20 = dphiphipsi3d(i,j,k)
+ o21 = detaqpsi3d(i,j,k)
+ o22 = dqpsi3d(i,j,k)
+ o23 = detaetapsi3d(i,j,k)
+ o24 = tan(q(i,j,k))
+ o25 = o24**2
+ o26 = 1/o25
+ o27 = dqphipsi3d(i,j,k)
+ o28 = 1/o24
+ psixx(i,j,k) = o2*o4*(1.00000000000000d0*o11*o16*o19*o20 + 1.000
+ & 00000000000d0*o11*o16*o22*o28 - 5.0000000000000d-1*o11*o16*o3 -
+ & 2.00000000000000d0*o11*o13*o14*o5 + 2.00000000000000d0*o11*o14*o
+ & 15*o5 + 1.00000000000000d0*o11*o14*o15*o19*o5 + 1.00000000000000
+ & d0*o11*o14*o15*o26*o5 - 2.00000000000000d0*o11*o14*o27*o28*o5 +
+ & o11*o18*o23*o6 + 7.5000000000000d-1*o11*o18*o3*o6 + 2.0000000000
+ & 0000d0*o11*o17*o21*o6*o7 - 3.00000000000000d0*o11*o17*o22*o6*o7
+ & + o11*o12*o6*o8 - 5.0000000000000d-1*o11*o3*o6*o8 + o11*o16*o9 -
+ & 2.00000000000000d0*o11*o18*o6*o9 + o11*o6*o8*o9)
+ psixy(i,j,k) = o2*o4*(-(o11*o13*o16) + 1.50000000000000d0*o11*o1
+ & 5*o16 + 1.00000000000000d0*o11*o15*o16*o26 - o11*o16*o27*o28 - o
+ & 11*o14*o19*o20*o5 + o11*o14*o18*o23*o5 - o11*o14*o22*o28*o5 + 5.
+ & 0000000000000d-1*o11*o14*o3*o5 + 7.5000000000000d-1*o11*o14*o18*
+ & o3*o5 + o11*o13*o6 - 5.0000000000000d-1*o11*o15*o6 - o11*o15*o19
+ & *o6 + 1.00000000000000d0*o11*o27*o28*o6 + 2.00000000000000d0*o11
+ & *o14*o17*o21*o5*o7 - 3.00000000000000d0*o11*o14*o17*o22*o5*o7 +
+ & o11*o12*o14*o5*o8 - 5.0000000000000d-1*o11*o14*o3*o5*o8 - o11*o1
+ & 4*o5*o9 - 2.00000000000000d0*o11*o14*o18*o5*o9 + o11*o14*o5*o8*o
+ & 9)
+ psixz(i,j,k) = o2*o4*(o11*o14*o27 - o11*o13*o14*o28 + 5.00000000
+ & 00000d-1*o11*o14*o15*o28 - o11*o18*o21*o5 + 1.50000000000000d0*o
+ & 11*o18*o22*o5 - o11*o12*o17*o5*o7 + o11*o17*o23*o5*o7 + 1.250000
+ & 00000000d0*o11*o17*o3*o5*o7 + o11*o21*o5*o8 - 1.50000000000000d0
+ & *o11*o22*o5*o8 - 3.00000000000000d0*o11*o17*o5*o7*o9)
+ psiyy(i,j,k) = o2*o4*(o11*o16*o18*o23 + 7.5000000000000d-1*o11*o
+ & 16*o18*o3 + 2.00000000000000d0*o11*o13*o14*o5 - 2.00000000000000
+ & d0*o11*o14*o15*o5 - o11*o14*o15*o19*o5 - o11*o14*o15*o26*o5 + 2.
+ & 00000000000000d0*o11*o14*o27*o28*o5 + 1.00000000000000d0*o11*o19
+ & *o20*o6 + 1.00000000000000d0*o11*o22*o28*o6 - 5.0000000000000d-1
+ & *o11*o3*o6 + 2.00000000000000d0*o11*o16*o17*o21*o7 - 3.000000000
+ & 00000d0*o11*o16*o17*o22*o7 + o11*o12*o16*o8 - 5.0000000000000d-1
+ & *o11*o16*o3*o8 - 2.00000000000000d0*o11*o16*o18*o9 + o11*o6*o9 +
+ & o11*o16*o8*o9)
+ psiyz(i,j,k) = o2*o4*(-(o11*o14*o18*o21) + 1.50000000000000d0*o1
+ & 1*o14*o18*o22 - o11*o27*o5 + 1.00000000000000d0*o11*o13*o28*o5 -
+ & 5.0000000000000d-1*o11*o15*o28*o5 - o11*o12*o14*o17*o7 + o11*o1
+ & 4*o17*o23*o7 + 1.25000000000000d0*o11*o14*o17*o3*o7 + o11*o14*o2
+ & 1*o8 - 1.50000000000000d0*o11*o14*o22*o8 - 3.00000000000000d0*o1
+ & 1*o14*o17*o7*o9)
+ psizz(i,j,k) = o2*o4*(o11*o12*o18 - 5.0000000000000d-1*o11*o18*o
+ & 3 - 2.00000000000000d0*o11*o17*o21*o7 + 3.00000000000000d0*o11*o
+ & 17*o22*o7 + o11*o23*o8 + 7.5000000000000d-1*o11*o3*o8 + o11*o18*
+ & o9 - 2.00000000000000d0*o11*o8*o9)
diff --git a/src/qfunc.x b/src/qfunc.x
new file mode 100644
index 0000000..d26919d
--- /dev/null
+++ b/src/qfunc.x
@@ -0,0 +1,29 @@
+ qf(i,j,k) = amp*(1.00000000000000d0+c*cos(phigrd(k))**2)*
+ $ (exp(-((etagrd(i)-eta0)**2/sigma**2))+exp(-((etagrd(i)+eta0)**2
+ $ /sigma**2)))*sin(qgrd(j))**n
+
+ qfetaeta(i,j,k) = amp*(1.00000000000000d0+c*cos(phigrd(k))
+ $ **2)*((4.0000000000000d0*(etagrd(i)-eta0)**2*exp(-((etagrd(i)-eta0)
+ $ **2/sigma**2)))/sigma**4-(2.00000000000000d0*exp(-((etagrd(i)
+ $ -eta0)**2/sigma**2)))/sigma**2+(4.0000000000000d0*(etagrd(i)+
+ $ eta0)**2*exp(-((etagrd(i)+eta0)**2/sigma**2)))/sigma**4-(2.00
+ $ 000000000000d0*exp(-((etagrd(i)+eta0)**2/sigma**2)))/sigma**2
+ $ )*sin(qgrd(j))**n
+
+ qfqq(i,j,k) = amp*(-1.00000000000000d0+n)*n*(1.00000000
+ $ 000000d0+c*cos(phigrd(k))**2)*cos(qgrd(j))**2*(exp(-((etagrd(i)-eta0)**2/
+ $ sigma**2))+exp(-((etagrd(i)+eta0)**2/sigma**2)))*sin(qgrd(j))**(-2.
+ $ 00000000000000d0+n)-amp*n*(1.00000000000000d0+c*cos(phigrd(k))
+ $ **2)*(exp(-((etagrd(i)-eta0)**2/sigma**2))+exp(-((etagrd(i)+eta0)**2
+ $ /sigma**2)))*sin(qgrd(j))**n
+
+ qfphi(i,j,k) = -2.00000000000000d0*amp*c*cos(phigrd(k))*(exp(-
+ $ ((etagrd(i)-eta0)**2/sigma**2))+exp(-((etagrd(i)+eta0)**2/sigma**2))
+ $ )*sin(phigrd(k))*sin(qgrd(j))**n
+
+ qfphiphi(i,j,k) = -2.00000000000000d0*amp*c*cos(phigrd(k))**2*
+ $ (exp(-((etagrd(i)-eta0)**2/sigma**2))+exp(-((etagrd(i)+eta0)**2/
+ $ sigma**2)))*sin(qgrd(j))**n+2.00000000000000d0*amp*c*(exp(-((
+ $ etagrd(i)-eta0)**2/sigma**2))+exp(-((etagrd(i)+eta0)**2/sigma**2)))*
+ $ sin(phigrd(k))**2*sin(qgrd(j))**n
+
diff --git a/test/test_dbh.par b/test/test_dbh.par
new file mode 100755
index 0000000..5e7c431
--- /dev/null
+++ b/test/test_dbh.par
@@ -0,0 +1,55 @@
+######################################################################
+#!DESC "Initial data for axisymmetric distored black hole"
+######################################################################
+
+ActiveThorns = "time iobasic ADMconstraints pugh interp cartgrid3d einstein ADM DistortedBHIVP ioascii ioutil"
+
+# GENERAL
+
+einstein::evolution_system = "ADM"
+
+driver::global_nx = 32
+driver::global_ny = 32
+driver::global_nz = 32
+
+grid::type = "byspacing"
+grid::dxyz = 0.4
+grid::domain = "full"
+
+time::dtfac = 0.25
+
+cactus::cctk_initial_time = 0.
+cactus::cctk_itlast = 1
+
+adm::method = "stagleap"
+adm::bound = "static"
+admconstraints::bound = "static"
+
+# MODEL
+
+einstein::initial_data = "distortedbh"
+distortedbhivp::amp = -0.1
+distortedbhivp::eta0 = 0.0
+distortedbhivp::c = 0.5
+distortedbhivp::sigma = 1.0
+distortedbhivp::etamax = 6.0
+distortedbhivp::n = 4
+distortedbhivp::ne = 102
+distortedbhivp::nq = 27
+distortedbhivp::np = 27
+interp::order = 1
+
+
+# GAUGE
+einstein::slicing = "geodesic"
+
+# OUTPUT ########################################################
+IO::outdir = "test_dbh"
+IOASCII::out1d_every = 1
+IOBasic::outScalar_vars = "einstein::metric einstein::curv admconstraints::ADMconstraints"
+IOASCII::out1D_vars = "einstein::metric einstein::curv einstein::lapse admconstraints::ADMconstraints einstein::grr einstein::psi einstein::psix einstein::psiy einstein::psiz einstein::psixx einstein::psixy einstein::psixz einstein::psiyy einstein::psiyz einstein::psizz"
+IOBasic::outInfo_every = 1
+IOBasic::outInfo_vars = "admconstraints::ham einstein::gxx"
+
+##################################################################
+
generated by cgit v1.2.3 (git 2.39.1) at 2024-05-05 17:04:44 +0000