From f1cbd3e0e16cdfb29bdcf18ab92a01e263acd4c9 Mon Sep 17 00:00:00 2001 From: rideout Date: Fri, 3 May 2002 17:50:32 +0000 Subject: IDAxiBrillBH converted to CactusEinstein2 git-svn-id: http://svn.einsteintoolkit.org/cactus/EinsteinInitialData/IDAxiBrillBH/trunk@39 0a4070d5-58f5-498f-b6c0-2693e757fa0f --- README | 5 +- interface.ccl | 6 +- param.ccl | 7 +- schedule.ccl | 9 +- src/AxiBrillBHIVP.F | 339 ------------------------------------------------- src/IDAxiBrillBH.F | 5 +- src/ParamCheck.c | 64 ++++++++++ src/make.code.defn | 5 +- test/test_axibrill.par | 19 ++- 9 files changed, 94 insertions(+), 365 deletions(-) delete mode 100644 src/AxiBrillBHIVP.F create mode 100644 src/ParamCheck.c diff --git a/README b/README index f2c0a9d..c0f2baf 100644 --- a/README +++ b/README @@ -3,6 +3,5 @@ Authors : Steve Brandt, Paul Walker, Ryoji Takahashi CVS info : $Header$ -------------------------------------------------------------------------- -This thorn does calculates initial data for a black hole distorted -by a even parity perturbation - +This thorn calculates initial data for a black hole distorted by an +even parity perturbation. diff --git a/interface.ccl b/interface.ccl index 7954581..1ed1806 100644 --- a/interface.ccl +++ b/interface.ccl @@ -1,9 +1,7 @@ -# Interface definition for thorn AxiBrillBHIVP +# Interface definition for thorn IDAxiBrillBH # $Header$ implements: idaxibrillbh -inherits: einstein +inherits: ADMBase grid StaticConformal public: - - diff --git a/param.ccl b/param.ccl index 47f98d2..e458cf7 100644 --- a/param.ccl +++ b/param.ccl @@ -1,7 +1,7 @@ -# Parameter definitions for thorn AxiBrillBHIVP +# Parameter definitions for thorn IDAxiBrillBH # $Header$ -shares:einstein +shares: ADMBase EXTENDS KEYWORD initial_lapse { @@ -13,6 +13,9 @@ EXTENDS KEYWORD initial_data "axibrillbh" :: "Axisymmetry Initial data for Black hole + Brill wave" } +USES KEYWORD metric_type + + private: REAL amp "Brill wave amplitude" diff --git a/schedule.ccl b/schedule.ccl index ac43b1c..e04c602 100644 --- a/schedule.ccl +++ b/schedule.ccl @@ -1,9 +1,14 @@ -# Schedule definitions for thorn AxiBrillBHIVP +# Schedule definitions for thorn IDAxiBrillBH # $Header$ if (CCTK_Equals(initial_data,"axibrillbh")) { - schedule IDAxiBrillBH at CCTK_INITIAL + schedule IDAxiBrillBH_ParamChecker at CCTK_PARAMCHECK + { + LANG: C + } "Check that the metric_type is recognised" + + schedule IDAxiBrillBH in ADMBase_InitialData { LANG: Fortran } "Construct IDAxiBrillBH" diff --git a/src/AxiBrillBHIVP.F b/src/AxiBrillBHIVP.F deleted file mode 100644 index 6e35d1b..0000000 --- a/src/AxiBrillBHIVP.F +++ /dev/null @@ -1,339 +0,0 @@ -c/*@@ -c @file AxiBrillBHIVP.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 AxiBrillBHIVP -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 AxiBrillBHIVP(CCTK_ARGUMENTS) - - implicit none - - DECLARE_CCTK_ARGUMENTS - DECLARE_CCTK_PARAMETERS - -c Perhaps this and others should go into cctk.h - integer CCTK_Equals - - - real*8 axibheps, rmax, dq, deta - integer levels,id5,id9,idi,idg,ier - real*8, allocatable :: cc(:,:),ce(:,:),cw(:,:),cn(:,:),cs(:,:), - $ rhs(:,:),psi2d(:,:),detapsi2d(:,:),dqpsi2d(:,:), - $ detaetapsi2d(:,:),detaqpsi2d(:,:),dqqpsi2d(:,:) - real*8, allocatable :: etagrd(:),qgrd(:) - real*8, allocatable :: eta(:,:,:),abseta(:,:,:),sign_eta(:,:,:), - $ q(:,:,:),phi(:,:,:) - real*8, allocatable :: psi2dv(:,:,:),detapsi2dv(:,:,:), - $ dqpsi2dv(:,:,:),detaetapsi2dv(:,:,:),detaqpsi2dv(:,:,:), - $ dqqpsi2dv(:,:,:) - parameter(axibheps = 1.0d-12) - real*8 ep1,ep2 - real*8 o1,o2,o3,o4,o5,o6,o7,o8,o9 - real*8 o10,o11,o12,o13,o14,o15,o16,o17,o18,o19 - real*8 o20,o21,o22,o23,o24,o25,o26,o27,o28,o29 - real*8 o30,o31,o32,o33,o34,o35,o36,o37,o38,o39 - real*8 o40,o41,o42,o43,o44,o45,o46,o47,o48,o49 - real*8 o50,o51,o52,o53,o54,o55,o56,o57,o58,o59 - real*8 o60,o61,o62,o63,o64,o65,o66,o67,o68,o69 - real*8 o70,o71,o72,o73,o74,o75,o76,o77,o78,o79 - real*8 o80,o81,o82,o83,o84,o85,o86,o87,o88,o89 - real*8 o90,o91,o92,o93,o94,o95,o96,o97,o98,o99 - real*8 pi - real*8 adm - integer :: nx,ny,nz - integer i,j,k,nquads - integer npoints,handle,ierror - - 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 Brill wave parameters - - print *,"Brill wave + BH Axisymmetric solve" - write (*,123)amp,eta0,sigma,n - print *,'etamax=',etamax - 123 format(1x,'Pars : Amp',f8.5,' eta0',f8.5,' sigma',f8.5,' n ',i3) - -c Solve on this sized cartesian grid -c 2D grid size NE x NQ -c Add 2 zones for boundaries... - ne = ne+2 - nq = nq+2 - ! do I need to call free? - allocate(cc(ne,nq),ce(ne,nq),cw(ne,nq),cn(ne,nq),cs(ne,nq), - $ rhs(ne,nq),psi2d(ne,nq),detapsi2d(ne,nq),dqpsi2d(ne,nq), - $ detaetapsi2d(ne,nq),detaqpsi2d(ne,nq),dqqpsi2d(ne,nq), - $ etagrd(ne),qgrd(nq)) - allocate(eta(nx,ny,nz),abseta(nx,ny,nz),sign_eta(nx,ny,nz), - $ q(nx,ny,nz),phi(nx,ny,nz), - $ psi2dv(nx,ny,nz),detapsi2dv(nx,ny,nz),dqpsi2dv(nx,ny,nz), - $ detaetapsi2dv(nx,ny,nz),detaqpsi2dv(nx,ny,nz), - $ dqqpsi2dv(nx,ny,nz)) -c Initialize some arrays - psi2d = 1.0d0 - detapsi2d = 0.0d0 - - nquads = 2 - dq = nquads*0.5*pi/(nq-2) - deta = etamax/(ne-3) - - do j=1,nq - qgrd(j) = (j-1.5)*dq - do i=1,ne - etagrd(i) = (i-2)*deta -#include "Development/IDAxiBrillBH/src/bhbrill.x" - enddo - enddo -c Boundary conditions - do j=1,nq - ce(2,j)=ce(2,j)+cw(2,j) - cw(2,j)=0.0 - - cw(ne-1,j)=cw(ne-1,j)+ce(ne-1,j) - cc(ne-1,j)=cc(ne-1,j)-deta*ce(ne-1,j) - ce(ne-1,j)=0.0 - - enddo - do i=1,ne - cc(i,2)=cc(i,2)+cs(i,2) - cs(i,2)=0.0 - cc(i,nq-1)=cc(i,nq-1)+cn(i,nq-1) - cn(i,nq-1)=0.0 - enddo - -c Do the solve - print *, " Calling axisymmetric solver" - call mgparm (levels,5,id5,id9,idi,idg,ne,nq) - call mg5 (ne,2,ne-1,nq,2,nq-1, - $ cc,cn,cs,cw,ce,psi2d,rhs, - $ id5,id9,idi,idg,1,axibheps,rmax,ier) - print *, " Solve complete" -c The solution is now available. -c Debugging is needed, a stop statement should -c be called if the IVP solve is not successful - - if(ier .ne. 0) stop 'bad solution to brill wave problem' - print *,'rmax = ',rmax - print *,'axibheps = ',axibheps - print *,'psi2d = ',maxval(psi2d),' ',minval(psi2d) - - ep2 = 0.0 - do j=2,nq-1 - do i=2,ne-1 - ep1 = rhs(i,j)-psi2d(i,j)*cc(i,j)-psi2d(i,j+1)*cn(i,j)-psi2d(i,j-1)*cs(i,j)- - & psi2d(i+1,j)*ce(i,j)-psi2d(i-1,j)*cw(i,j) - ep2 = max(abs(ep1),ep2) - enddo - enddo - print *,'Resulting eps =',ep2 - - ! what a pain all this is.... - do j=1,nq - psi2d(1,j)=psi2d(3,j) - psi2d(ne,j)=-deta*psi2d(ne-1,j)+psi2d(ne-2,j) - enddo - do i=1,ne - psi2d(i,1)=psi2d(i,2) - psi2d(i,nq)=psi2d(i,nq-1) - enddo -c goto 111 - do j=2,nq-1 - do i=2,ne-1 - dqpsi2d(i,j)=0.5*(psi2d(i,j+1)-psi2d(i,j-1))/dq - dqqpsi2d(i,j)=(psi2d(i,j+1)+psi2d(i,j-1)-2.*psi2d(i,j))/dq**2 - detapsi2d(i,j)=sinh(0.5*etagrd(i))+0.5*(psi2d(i+1,j)-psi2d(i-1,j))/deta - detaetapsi2d(i,j)=0.5*cosh(0.5*etagrd(i))+ - $ (psi2d(i+1,j)+psi2d(i-1,j)-2.*psi2d(i,j))/deta**2 - enddo - enddo - do j=1,nq - detapsi2d(1,j)=-detapsi2d(3,j) - detapsi2d(ne,j)=detapsi2d(ne-2,j) ! simplified - - detaetapsi2d(1,j)=detaetapsi2d(3,j) - detaetapsi2d(ne,j)=detaetapsi2d(ne-2,j) ! simplified... - - dqqpsi2d(1,j)=dqqpsi2d(3,j) - dqqpsi2d(ne,j)=dqqpsi2d(ne-2,j) ! simplified - - dqpsi2d(1,j)=dqpsi2d(3,j) - dqpsi2d(ne,j)=-dq*dqpsi2d(ne-1,j)+dqpsi2d(ne-2,j) - enddo - do i=1,ne - detapsi2d(i,1)=detapsi2d(i,2) - detapsi2d(i,nq)=detapsi2d(i,nq-1) - - detaetapsi2d(i,1)=detaetapsi2d(i,2) - detaetapsi2d(i,nq)=detaetapsi2d(i,nq-1) - - dqqpsi2d(i,1)=dqqpsi2d(i,2) - dqqpsi2d(i,nq)=dqqpsi2d(i,nq-1) - - dqpsi2d(i,1)=-dqpsi2d(i,2) - dqpsi2d(i,nq)=-dqpsi2d(i,nq-1) - enddo - do j=2,nq-1 - do i=2,ne-1 - detaqpsi2d(i,j)=0.5*(detapsi2d(i,j+1)-detapsi2d(i,j-1))/dq - enddo - enddo - do j=1,nq - detaqpsi2d(1,j)=-detaqpsi2d(3,j) - detaqpsi2d(ne,j)=detaqpsi2d(ne-2,j) ! simplified - enddo - do i=1,ne - detaqpsi2d(i,1)=-detaqpsi2d(i,2) - detaqpsi2d(i,nq)=-detaqpsi2d(i,nq-1) - enddo - do j=1,nq - psi2d(:,j)=psi2d(:,j)+2.*cosh(0.5*etagrd) - enddo - -c Now evaluate each of the following at x(i,j,k), y(i,j,k) and -c z(i,j,k) where i,j,k go from 1 to nx,ny,nz - -c Conformal factor - - 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 -c eta(i,j,k) = 0.5d0*dlog(x(i,j,k)**2+y(i,j,k)**2+z(i,j,k)**2) -c abseta(i,j,k) = abs(eta(i,j,k)) - if(eta(i,j,k) .lt. 0)then - sign_eta(i,j,k) = -1 - else - sign_eta(i,j,k) = 1 - endif -c q(i,j,k) = datan2(sqrt(x(i,j,k)**2+y(i,j,k)**2),z(i,j,k)) -c TYPO HERE ??????????? -c | -c | -c phi(i,j,k)= datan2(y(i,j,k),x(i,j,k)) - enddo - enddo - enddo - - call CCTK_InterpHandle (handle, "simple_local") - - npoints = nx*ny*nz - - call CCTK_Interp (ierror,cctkGH,handle,npoints,2,6,6, - $ ne,nq,abseta,q, - $ CCTK_VARIABLE_REAL,CCTK_VARIABLE_REAL, - $ etagrd(1),qgrd(1),deta,dq, - $ psi2d,detapsi2d,dqpsi2d,detaetapsi2d,detaqpsi2d,dqqpsi2d, - $ CCTK_VARIABLE_REAL,CCTK_VARIABLE_REAL,CCTK_VARIABLE_REAL, - $ CCTK_VARIABLE_REAL,CCTK_VARIABLE_REAL,CCTK_VARIABLE_REAL, - $ psi2dv,detapsi2dv,dqpsi2dv,detaetapsi2dv,detaqpsi2dv, - $ dqqpsi2dv, - $ CCTK_VARIABLE_REAL,CCTK_VARIABLE_REAL,CCTK_VARIABLE_REAL, - $ CCTK_VARIABLE_REAL,CCTK_VARIABLE_REAL,CCTK_VARIABLE_REAL) - - psi = psi2dv * exp (-0.5 * eta) - detapsi2dv = sign_eta * detapsi2dv - detaqpsi2dv = sign_eta * detaqpsi2dv - - do k=1,nz - do j=1,ny - do i=1,nx - -c psi(i,j,k) = psi2dv(i,j,k)*exp(-0.5*eta(i,j,k)) -c detapsi2dv(i,j,k) = sign_eta(i,j,k)*detapsi2dv(i,j,k) - -c psix = \partial psi / \partial x / psi -#include "Development/IDAxiBrillBH/src/psi_1st_deriv.x" - -c detaqpsi2dv(i,j,k) = sign_eta(i,j,k)*detaqpsi2dv(i,j,k) - -c psixx = \partial^2\psi / \partial x^2 / psi -#include "Development/IDAxiBrillBH/src/psi_2nd_deriv.x" - enddo - enddo - enddo - -c Conformal metric -c gxx = ... - -c Derivatives of the metric -c dxgxx = 1/2 \partial gxx / \partial x - - do k=1,nz - do j=1,ny - do i=1,nx -c THESE WERE ALREADY CALCULATED ABOVE !!! -c eta(i,j,k) = 0.5d0*dlog(x(i,j,k)**2+y(i,j,k)**2+z(i,j,k)**2) -c q(i,j,k) = datan2(sqrt(x(i,j,k)**2+y(i,j,k)**2),z(i,j,k)) -c phi(i,j,k) = datan2(y(i,j,k),x(i,j,k)) -#include "Development/IDAxiBrillBH/src/gij.x" - enddo - enddo - enddo - -c Curvature - kxx = 0.0D0 - kxy = 0.0D0 - kxz = 0.0D0 - kyy = 0.0D0 - kyz = 0.0D0 - kzz = 0.0D0 - - 111 continue -c Set ADM mass - i = ne-15 - adm = 0.0 - do j=2,nq-1 - adm=adm+(psi2d(i,j)-(psi2d(i+1,j)-psi2d(i-1,j))/deta)*exp(0.5* - $ etagrd(i)) - enddo - adm=adm/(nq-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 - - conformal_state = CONFORMAL_METRIC - - deallocate(cc,ce,cw,cn,cs,rhs,psi2d,detapsi2d,dqpsi2d, - $ detaetapsi2d,detaqpsi2d,dqqpsi2d, - $ etagrd,qgrd, - $ eta,abseta,sign_eta,q,phi,psi2dv,detapsi2dv,dqpsi2dv, - $ detaetapsi2dv,detaqpsi2dv,dqqpsi2dv) - - return - end - diff --git a/src/IDAxiBrillBH.F b/src/IDAxiBrillBH.F index 84f510d..2a72bd7 100644 --- a/src/IDAxiBrillBH.F +++ b/src/IDAxiBrillBH.F @@ -350,7 +350,9 @@ c Set ADM mass alp = (2.*r - adm)/(2.*r+adm) endif - conformal_state = CONFORMAL_METRIC +c conformal_state = CONFORMAL_METRIC (What is CONFORMAL_METRIC?) + conformal_state = 3 +c 3 ==> 'all' derivatives were calculated deallocate(cc,ce,cw,cn,cs,rhs,psi2d,detapsi2d,dqpsi2d, $ detaetapsi2d,detaqpsi2d,dqqpsi2d, @@ -360,4 +362,3 @@ c Set ADM mass return end - diff --git a/src/ParamCheck.c b/src/ParamCheck.c new file mode 100644 index 0000000..f2259c0 --- /dev/null +++ b/src/ParamCheck.c @@ -0,0 +1,64 @@ + /*@@ + @file ParamCheck.c + @date Thu May 2 19:23:47 CEST 2002 + @author David Rideout + @desc + Check the parameters for IDAxiBrillBH + @enddesc + @version $Header$ + @@*/ + +#include "cctk.h" + +#include "cctk_Arguments.h" +#include "cctk_Parameters.h" + +static const char *rcsid = "$Header$"; + +CCTK_FILEVERSION(CactusEinstein_IDAxiBrillBH_ParamCheck_c) + +/******************************************************************** + ********************* Local Data Types *********************** + ********************************************************************/ + +/******************************************************************** + ********************* Local Routine Prototypes ********************* + ********************************************************************/ + +/******************************************************************** + ***************** Scheduled Routine Prototypes ********************* + ********************************************************************/ + +void IDAxiBrillBH_ParamChecker(CCTK_ARGUMENTS); + +/******************************************************************** + ********************* Other Routine Prototypes ********************* + ********************************************************************/ + +/******************************************************************** + ********************* Local Data ***************************** + ********************************************************************/ + +/******************************************************************** + ********************* External Routines ********************** + ********************************************************************/ + +void IDAxiBrillBH_ParamChecker(CCTK_ARGUMENTS) +{ + DECLARE_CCTK_ARGUMENTS + DECLARE_CCTK_PARAMETERS + + /* Do we know how to deal with this type of metric ? */ + if( ! CCTK_EQUALS(metric_type, "static conformal")) + if (CCTK_EQUALS(metric_type, "physical")) + { + CCTK_PARAMWARN("\n\tPlease add code into IDAxiBrillBH.F to compute\n\tthe physical metric from the conformal metric."); + } else + { + CCTK_PARAMWARN("Unknown ADMBase::metric_type - known types are \"physical\" and \"static conformal\""); + } +} + +/******************************************************************** + ********************* Local Routines ************************* + ********************************************************************/ diff --git a/src/make.code.defn b/src/make.code.defn index 7bd9c1e..6671b79 100644 --- a/src/make.code.defn +++ b/src/make.code.defn @@ -1,9 +1,8 @@ -# Main make.code.defn file for thorn AxiBrillBHIVP +# Main make.code.defn file for thorn IDAxiBrillBH # $Header$ # Source files in this directory -SRCS = IDAxiBrillBH.F mg59p.F shmgp.F77 +SRCS = IDAxiBrillBH.F mg59p.F shmgp.F77 ParamCheck.c # Subdirectories containing source files SUBDIRS = - diff --git a/test/test_axibrill.par b/test/test_axibrill.par index 6b5029c..94b161f 100644 --- a/test/test_axibrill.par +++ b/test/test_axibrill.par @@ -2,14 +2,15 @@ !DESC "Initial data for axisymmetric distored black hole" ###################################################################### -ActiveThorns = "Boundary time iobasic ADMconstraints pugh pughslab pughreduce pughinterp cartgrid3d einstein ADM IDAxiBrillBH ioascii ioutil" +ActiveThorns = "Boundary time iobasic ADMconstraints pugh pughslab pughreduce pughinterp cartgrid3d ADMBase StaticConformal ADM admmacros coordgauge ADMAnalysis IDAxiBrillBH ioascii ioutil" #pugh::enable_all_storage = "yes" #cactus::cctk_show_banners = 0 -# GENERAL +# GENERAL -einstein::evolution_system = "ADM" +ADMBase::evolution_method = "ADM" +ADMBase::metric_type = "static conformal" driver::global_nx = 32 driver::global_ny = 32 @@ -27,9 +28,8 @@ cactus::cctk_itlast = 1 adm::method = "stagleap" adm::bound = "static" -# MODEL - -einstein::initial_data = "axibrillbh" +# MODEL +ADMBase::initial_data = "axibrillbh" idaxibrillbh::amp = 0.5 idaxibrillbh::eta0 = 0.0 idaxibrillbh::sigma = 1.0 @@ -41,16 +41,16 @@ idaxibrillbh::nq = 27 idaxibrillbh::interpolation_order = 1 # GAUGE -einstein::slicing = "geodesic" +ADMBase::lapse_evolution_method = "geodesic" # OUTPUT ######################################################## IO::outdir = "test_axibrill" IO::out_fileinfo = "none" IOASCII::out1d_every = 1 -IOASCII::out1D_vars = "einstein::curv einstein::metric einstein::grr admconstraints::ham" +IOASCII::out1D_vars = "ADMBase::curv admbase::metric ADMAnalysis::grr admconstraints::ham" IOBasic::outInfo_every = 1 -IOBasic::outInfo_vars = "einstein::gxx" +IOBasic::outInfo_vars = "admbase::gxx" IOASCII::out1D_xline_y = 0.2 IOASCII::out1D_xline_z = 0.2 @@ -62,4 +62,3 @@ IOASCII::out1D_zline_x = 0.2 IOASCII::out1D_zline_y = 0.2 ################################################################## - -- cgit v1.2.3