From 6b266c308ed6ed21b401e3a670cc017286d79036 Mon Sep 17 00:00:00 2001 From: ryoji Date: Wed, 1 Dec 1999 13:36:25 +0000 Subject: 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 --- README | 20 ++ interface.ccl | 9 + param.ccl | 65 +++++ schedule.ccl | 14 + src/DistortedBHIVP.F | 665 +++++++++++++++++++++++++++++++++++++++++++++ src/DistortedBHIVP.F.sav | 681 +++++++++++++++++++++++++++++++++++++++++++++++ src/Stab3d.F | 280 +++++++++++++++++++ src/bhbrill3d.m | 132 +++++++++ src/bhbrill3d.x | 126 +++++++++ src/gij.x | 117 ++++++++ src/interp3.F | 264 ++++++++++++++++++ src/make.code.defn | 9 + src/psi_1st_deriv.x | 20 ++ src/psi_2nd_deriv.x | 71 +++++ src/qfunc.x | 29 ++ test/test_dbh.par | 55 ++++ 16 files changed, 2557 insertions(+) create mode 100644 README create mode 100644 interface.ccl create mode 100644 param.ccl create mode 100644 schedule.ccl create mode 100644 src/DistortedBHIVP.F create mode 100644 src/DistortedBHIVP.F.sav create mode 100644 src/Stab3d.F create mode 100644 src/bhbrill3d.m create mode 100644 src/bhbrill3d.x create mode 100644 src/gij.x create mode 100644 src/interp3.F create mode 100644 src/make.code.defn create mode 100644 src/psi_1st_deriv.x create mode 100644 src/psi_2nd_deriv.x create mode 100644 src/qfunc.x create mode 100755 test/test_dbh.par 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

+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" + +################################################################## + -- cgit v1.2.3