From 59c381790dc39733e20107641dca2b772ca49080 Mon Sep 17 00:00:00 2001 From: ryoji Date: Wed, 9 Mar 2005 06:57:33 +0000 Subject: fixed of bug in schedule.ccl. git-svn-id: http://svn.einsteintoolkit.org/cactus/EinsteinInitialData/IDAxiOddBrillBH/trunk@51 b6f3ac56-194f-0410-8878-cdf6079d7f1b --- interface.ccl | 2 +- param.ccl | 33 +++++++++++++++++++++++++++++---- schedule.ccl | 2 +- src/IDAxiOddBrillBH.F | 43 ++++++++++++++++++++++++++++++++----------- src/gij.x | 49 +++++++++++++++++++++++++++---------------------- 5 files changed, 90 insertions(+), 39 deletions(-) diff --git a/interface.ccl b/interface.ccl index d8c3ed9..2f29f81 100644 --- a/interface.ccl +++ b/interface.ccl @@ -2,7 +2,7 @@ # $Header$ implements: axiodd -inherits: admbase staticconformal grid admanalysis driver +inherits: ADMBase grid StaticConformal private: diff --git a/param.ccl b/param.ccl index 7db98ea..e321d87 100644 --- a/param.ccl +++ b/param.ccl @@ -20,17 +20,37 @@ EXTENDS KEYWORD initial_data private: -REAL amp "Brill wave amplitude" +REAL amp "Brill wave amplitude for extrinsic curveture" { *:* :: "No restriction" } 0.1 -REAL eta0 "Brill wave center (in eta coords)" +REAL amp_me "Brill wave amplitude for metric" { *:* :: "No restriction" } 0.0 -REAL sigma "Brill wave width (in eta)" +REAL xi "for recoil" +{ + *:* :: "No restriction" +} 0.0 + +REAL eta0 "Brill wave center (in eta coords) for extrisic curvature" +{ + *:* :: "No restriction" +} 0.0 + +REAL eta0_me "Brill wave center (in eta coords) for metric" +{ + *:* :: "No restriction" +} 0.0 + +REAL sigma "Brill wave width (in eta) for extrinsic curveture" +{ + *:* :: "No restriction" +} 1.0 + +REAL sigma_me "Brill wave width (in eta) for metric" { *:* :: "No restriction" } 1.0 @@ -46,7 +66,12 @@ REAL etamax "Eta value for outer edge of grid" } 5.0 -INT n "sin**n theta in brill wave" +INT n "sin**n theta in brill wave for extrinsic curveture" +{ + *:* :: "No restriction" +} 3 + +INT n_me "sin**n theta in brill wave for metric" { *:* :: "No restriction" } 3 diff --git a/schedule.ccl b/schedule.ccl index 79254f1..c9c43df 100644 --- a/schedule.ccl +++ b/schedule.ccl @@ -10,7 +10,7 @@ if (CCTK_Equals(initial_data,"axioddbh")) LANG: C } "Check Parameters" - schedule IDAxiOddBrillBH at CCTK_INITIAL + schedule IDAxiOddBrillBH IN ADMBase_InitialData { LANG: Fortran } "Construct IDAxiOddBrillBH" diff --git a/src/IDAxiOddBrillBH.F b/src/IDAxiOddBrillBH.F index 838ac2f..781fd36 100644 --- a/src/IDAxiOddBrillBH.F +++ b/src/IDAxiOddBrillBH.F @@ -36,7 +36,8 @@ c@@*/ real*8 deta,dq real*8, allocatable :: ac(:,:),ae(:,:),aw(:,:),an(:,:),as(:,:), - $ rhs(:,:),ksq(:,:), psi2dv(:,:),dpsi2dv(:,:),ddpsi2dv(:,:), + $ rhs(:,:),qfetaeta(:,:),qfqq(:,:),ksq(:,:), + $ psi2dv(:,:),dpsi2dv(:,:),ddpsi2dv(:,:), $ detapsisph(:,:),dqpsisph(:,:),detaetapsisph(:,:), $ detaqpsisph(:,:),dqqpsisph(:,:), $ etagrd(:),qgrd(:) @@ -85,11 +86,17 @@ c -------------- write(*,*) '"Solve Bradt-Seidel Initial Data Sets"' endif - write(*,*) 'Brill Wave: amplitude :',amp - write(*,*) 'Brill Wave: center :',eta0 - write(*,*) 'Brill Wave: sigma :',sigma + write(*,*) 'Brill Wave: amplitude for K :',amp + write(*,*) 'Brill Wave: center for K:',eta0 + write(*,*) 'Brill Wave: sigma for K:',sigma write(*,*) 'Bowen & York Momenta: :',byJ - write(*,*) 'Brill Wave: sin^n theta :',n + write(*,*) 'Brill Wave: sin^n theta for K:',n + + write(*,*) 'Brill Wave: amplitude for g :',amp_me + write(*,*) 'Brill Wave: center for g:',eta0_me + write(*,*) 'Brill Wave: sigma for g:',sigma_me + write(*,*) 'Brill Wave: sin^n theta for g:',n_me + c Check if we should create and store conformal factor stuff @@ -124,7 +131,8 @@ c Add 2 zones for eta coordinate and 4 for theta nq = nq+2 c allocate(ac(ne,nq),ae(ne,nq),aw(ne,nq),an(ne,nq),as(ne,nq), - $ rhs(ne,nq),ksq(ne,nq),psi2dv(ne,nq),dpsi2dv(ne,nq), + $ rhs(ne,nq),qfetaeta(ne,nq),qfqq(ne,nq),ksq(ne,nq), + $ psi2dv(ne,nq),dpsi2dv(ne,nq), $ ddpsi2dv(ne,nq),detapsisph(ne,nq),dqpsisph(ne,nq), $ detaetapsisph(ne,nq),detaqpsisph(ne,nq),dqqpsisph(ne,nq), $ etagrd(ne),qgrd(nq)) @@ -141,7 +149,16 @@ c do i=1,ne etagrd(i) = (i-2)*deta enddo - +c +c Initialize q-function and its derivatives: should be generalized +c + do j = 1,nq + do i = 1,ne +#include "qfunc_even.x" + enddo + enddo + + c c Specify the initial guess for the conformal factor: c @@ -190,7 +207,9 @@ c $ (dpsi2dv(i+1,j)-2.*dpsi2dv(i,j)+dpsi2dv(i-1,j))/ $ deta**2+(dpsi2dv(i,j+1)-2.*dpsi2dv(i,j)+ $ dpsi2dv(i,j-1))/dq**2+0.5*(dpsi2dv(i,j+1)- - $ dpsi2dv(i,j-1))/(dq*tan(qgrd(j)))-0.25*dpsi2dv(i,j)+ + $ dpsi2dv(i,j-1))/(dq*tan(qgrd(j))) + $ +0.25*dpsi2dv(i,j)* + $ (qfetaeta(i,j)+qfqq(i,j)-1.0d0)+ $ 0.125*ksq(i,j)/(psi2dv(i,j)+dpsi2dv(i,j))**7 enddo enddo @@ -220,8 +239,9 @@ c Compute stencil coefficients: c do j = 2,nq-1 do i = 2,ne-1 - ac(i,j) = -2./deta**2-2./dq**2-sin(qgrd(j))**2-0.25- - $ 0.875*ksq(i,j)/(psi2dv(i,j)+ + ac(i,j) = -2./deta**2-2./dq**2-sin(qgrd(j))**2+0.25 + $ *(qfetaeta(i,j)+qfqq(i,j) - 1.0d0) + $ -0.875*ksq(i,j)/(psi2dv(i,j)+ $ dpsi2dv(i,j))**8 c ae(i,j) = 1./deta**2 @@ -549,7 +569,8 @@ c print *,'Angular momentum: ',Jmom alp = (2.*r - adm)/(2.*r+adm) endif - deallocate(ac,ae,aw,an,as,rhs,ksq,psi2dv,dpsi2dv,ddpsi2dv, + deallocate(ac,ae,aw,an,as,rhs,qfetaeta,qfqq, + $ ksq,psi2dv,dpsi2dv,ddpsi2dv, $ detapsisph,dqpsisph,detaetapsisph,detaqpsisph,dqqpsisph, $ etagrd,qgrd) diff --git a/src/gij.x b/src/gij.x index 693b57d..1f57c56 100644 --- a/src/gij.x +++ b/src/gij.x @@ -1,24 +1,29 @@ - gxx(i,j,k) = 1.00000000000000d0 -c dxgxx(i,j,k) = 0 -c dygxx(i,j,k) = 0 -c dzgxx(i,j,k) = 0 - gxy(i,j,k) = 0 -c dxgxy(i,j,k) = 0 -c dygxy(i,j,k) = 0 -c dzgxy(i,j,k) = 0 + o1 = 2.d0*phi(i,j,k) + o2 = cos(o1) + o3 = -xi + o4 = cos(q(i,j,k)) + o5 = o4*xi + o6 = 1.d0 + o3 + o5 + o7 = -eta0_me + o8 = eta(i,j,k) + o7 + o9 = o8**2 + o10 = sigma_me**2 + o11 = 1/o10 + o12 = -(o11*o9) + o13 = exp(o12) + o14 = eta(i,j,k) + eta0_me + o15 = o14**2 + o16 = -(o11*o15) + o17 = exp(o16) + o18 = o13 + o17 + o19 = sin(q(i,j,k)) + o20 = o19**n_me + o21 = 2.d0*amp_me*o18*o20*o6 + o22 = exp(o21) + o23 = -1.d0 + o22 + gxx(i,j,k) = 5.d-1*(1.d0 + o22 + o2*o23) + gxy(i,j,k) = o23*cos(phi(i,j,k))*sin(phi(i,j,k)) 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) = 1.00000000000000d0 -c dxgyy(i,j,k) = 0 -c dygyy(i,j,k) = 0 -c dzgyy(i,j,k) = 0 + gyy(i,j,k) = 5.d-1*(1.d0 + o22 - o2*o23) 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) = 1.00000000000000d0 -c dxgzz(i,j,k) = 0 -c dygzz(i,j,k) = 0 -c dzgzz(i,j,k) = 0 + gzz(i,j,k) = o22 -- cgit v1.2.3