diff options
-rw-r--r-- | param.ccl | 44 | ||||
-rw-r--r-- | src/brillq.F | 129 | ||||
-rw-r--r-- | src/finishbrilldata.F | 16 |
3 files changed, 99 insertions, 90 deletions
@@ -1,4 +1,6 @@ # Parameter definitions for thorn IDBrillData +# $Header$ + shares:einstein @@ -25,17 +27,17 @@ KEYWORD brill_bound "Which boundary condition to use" "robin" :: "Robin boundary: set brill_robin_falloff, brill_robin_inf" } "const" -CCTK_REAL brill_const_v0 "Value of constant BC" -{ - : :: "anything goes" -} 1.0 - -CCTK_INT brill_robin_falloff "Fall-off of Robin BC" +INT brill_robin_falloff "Fall-off of Robin BC" { 0: :: "any positive integer value" } 1 -CCTK_REAL brill_robin_inf "Value at infinity of Robin BC" +REAL brill_const_v0 "Value of constant BC" +{ + : :: "anything goes" +} 1.0 + +REAL brill_robin_inf "Value at infinity of Robin BC" { : :: "anything goes" } 1.0 @@ -55,22 +57,22 @@ BOOLEAN axisym "Axisymmetric initial data?" INT brill_q "Form of function q [0,1,2]" { 0:2 :: "Only cases 0,1,2 defined at the moment" -} 1 +} 2 REAL brill_a "Brill wave: Amplitude" { : :: "Anything" } 0.0 -REAL brill_b "Brill wave: rho^b" +INT brill_b "Brill wave: rho^b" { : :: "Anything" -} 2.0 +} 2 -REAL brill_c "Brill wave: (r^2 - r0^2)^(c/2)" +INT brill_c "Brill wave: (r^2 - r0^2)^(c/2)" { : :: "Anything" -} 2.0 +} 2 REAL brill_rho0 "Brill wave: radius of torus in rho" { @@ -92,6 +94,11 @@ REAL brill_sr "Brill wave: sigma in r" : :: "Anything" } 1.0 +REAL brill_sz "Brill wave: sigma in z" +{ + : :: "Anything" +} 1.0 + # 3D Brill wave parameters REAL brill_d "3D Brill wave: d rho^m cos^2(n (phi + phi0))" @@ -128,13 +135,14 @@ KEYWORD outputpsi "Output conformal factor?" # Additional parameters -REAL brill_eps "epsilon for finite differencing" -{ - 0: :: "Positive please" -} 1.e-6 - REAL brill_rhofudge "delta rho for axis fudge" { 0: :: "Positive please" -} 1.e-6 +} 0.00001 + + + + + + diff --git a/src/brillq.F b/src/brillq.F index 8d62064..ab299ec 100644 --- a/src/brillq.F +++ b/src/brillq.F @@ -10,44 +10,42 @@ #include "cctk.h" #include "cctk_Parameters.h" -#include "cctk_Arguments.h" function brillq(rho,z,phi) -C Calculates the function q that appear in the conformal -C metric for Brill waves: -C -C ds^2 = psi^4 ( e^(2q) (drho^2 + dz^2) + rho^2 dphi^2 ) -C -C There are three different choices for the form of q depending -C on the value of the parameter "brill_q": -C -C brill_q = 0: -C 2 -C 2+b - (rho-rho0) 2 2 -C q = a rho e (z/sz) / r -C -C -C brill_q = 1: (includes Eppleys form as special case) -C -C -C b 2 2 2 c/2 -C q = a (rho/srho) / { 1 + [ ( r - r0 ) / sr ] } -C -C -C brill_q = 2: (includes Holz et al form as special case) -C -C 2 2 2 c/2 -C b - [ ( r - r0 ) / sr ] -C q = a (rho/srho) e -C -C -C If we want 3D initial data, the function q is multiplied by an -C additional factor: -C -C m 2 m -C q -> q [ 1 + d rho cos (n (phi)+phi0 ) / ( 1 + e rho ) ] - +c Calculates the function q that appear in the conformal +c metric for Brill waves: +c +c ds^2 = psi^4 ( e^(2q) (drho^2 + dz^2) + rho^2 dphi^2 ) +c +c There are three different choices for the form of q depending +c on the value of the parameter "brill_q": +c +c brill_q = 0: +c 2 +c 2+b - (rho-rho0) 2 2 +c q = a rho e (z/sz) / r +c +c +c brill_q = 1: (includes Eppleys form as special case) +c +c +c b 2 2 2 c/2 +c q = a (rho/srho) / { 1 + [ ( r - r0 ) / sr ] } +c +c +c brill_q = 2: (includes Holz et al form as special case) +c +c 2 2 2 c/2 +c b - [ ( r - r0 ) / sr ] +c q = a (rho/srho) e +c +c +c If we want 3D initial data, the function q is multiplied by an +c additional factor: +c +c m 2 m +c q -> q [ 1 + d rho cos (n (phi)+phi0 ) / ( 1 + e rho ) ] implicit none @@ -56,17 +54,19 @@ C q -> q [ 1 + d rho cos (n (phi)+phi0 ) / ( 1 + e rho ) ] logical firstcall integer qtype + integer b,c - real*8 brillq,rho,z,phi - real*8 a,b,c,r0,sr,rho0,srho,brill_sz - real*8 d,e,m,n,phi0 + CCTK_REAL brillq,rho,z,phi + CCTK_REAL a,r0,sr,rho0,srho,sz + CCTK_REAL d,e,m,n,phi0 data firstcall /.true./ - save firstcall,qtype,a,b,c,r0,sr,rho0,srho,brill_sz + save firstcall,qtype + save a,b,c,r0,sr,rho0,srho,sz save d,e,m,n,phi0 -C Get parameters at first call. +c Get parameters at first call. if (firstcall) then @@ -80,7 +80,7 @@ C Get parameters at first call. sr = brill_sr rho0 = brill_rho0 srho = brill_srho - brill_sz = brill_sz + sz = brill_sz if (axisym.eq.0) then d = brill_d @@ -94,50 +94,51 @@ C Get parameters at first call. end if - if (rho.lt.0) call CCTK_WARN(1,"Warning: negative rho in brillq:") +c Calculate q(rho,z) from a choice of forms. -C Calculate q(rho,z) from a choice of forms. Type 0, 1 and 2 are those -C of Bruegmann. +c brill_q = 0. -C brill_q = 0. + if (qtype.eq.0) then - if (qtype .eq. 0) then + brillq = a*dabs(rho)**(2 + b) + . / dexp((rho - rho0)**2)/(rho**2 + z**2) - brillq = a * rho**(2.d0+b) - $ / dexp((rho - rho0)**2) / (rho**2 + z**2) - if (brill_sz .ne. 0.d0) then - brillq = brillq / exp((z/brill_sz)**2) + if (sz.ne.0.0D0) then + brillq = brillq/dexp((z/sz)**2) end if - else if (qtype .eq. 1) then + else if (qtype.eq.1) then + +c brill_q = 1. This includes Eppleys choice of q. +c But note that q(Eppley) = 2q(Cactus). -C brill_q = 1. This includes Eppleys choice of q. -C But note that q(Eppley) = 2q(Cactus). + brillq = a*(rho/srho)**b + . / ( 1.0D0 + ((rho**2 + z**2 - r0**2)/sr**2)**(c/2) ) - brillq = a * (rho/srho)**b - $ / ( 1.d0 + ((rho**2 + z**2 - r0**2) / sr**2)**(c/2) ) + else if (qtype.eq.2) then - else if (qtype .eq. 2) then +c brill_q = 2. This includes my (Carstens) notion of what a +c smooth "pure quadrupole" q should be, plus the choice of +c Holz et al. -C brill_q = 2. This includes my (Carstens) notion of what a -C smooth "pure quadrupole" q should be. + brillq = a*(rho/srho)**b + . / dexp(((rho**2 + z**2 - r0**2)/sr**2)**(c/2)) - brillq = a * (rho/srho)**b - $ / dexp( ((rho**2 + z**2 - r0**2) / sr**2)**(c/2) ) + print *,a,b,c,r0,sr,srho,z,rho,brillq else -C Unknown type for q function. +c Unknown type for q function. - call CCTK_WARN(0,"Unknown type of Brill function q.") + call CCTK_WARN(0,"Unknown type of Brill function q") end if -C If desired, multiply with phi dependent factor. +c If desired, multiply with a phi dependent factor. if (axisym.eq.0) then brillq = brillq*(1.0D0 + d*rho**m*cos(n*(phi-phi0))**2 - . /(1.0D0 + e*rho**m)) + . / (1.0D0 + e*rho**m)) end if return diff --git a/src/finishbrilldata.F b/src/finishbrilldata.F index 36060e9..597c3ae 100644 --- a/src/finishbrilldata.F +++ b/src/finishbrilldata.F @@ -25,15 +25,15 @@ integer nx,ny,nz CCTK_REAL x1,y1,z1,rho1,rho2 - CCTK_REAL brillq,psi4,e2q,rhofudge - CCTK_REAL phi,phif - CCTK_REAL zero + CCTK_REAL phi,psi4,e2q,rhofudge + CCTK_REAL zero,one - external brillq + CCTK_REAL brillq,phif c Numbers. zero = 0.0D0 + one = 1.0D0 c Set up grid size. @@ -74,12 +74,12 @@ c e^2q (dx^2 + dy^2 + dz^2) + (1-e^2q) (xdy-ydx)^2/rho^2 c c The individual coefficients can be read off as - if (rho1 .gt. rhofudge) then + if (rho1.gt.rhofudge) then - gxx(i,j,k) = psi4*(e2q + (1.d0 - e2q)*y1**2/rho2) - gyy(i,j,k) = psi4*(e2q + (1.d0 - e2q)*x1**2/rho2) + gxx(i,j,k) = psi4*(e2q + (one - e2q)*y1**2/rho2) + gyy(i,j,k) = psi4*(e2q + (one - e2q)*x1**2/rho2) gzz(i,j,k) = psi4*e2q - gxy(i,j,k) = - psi4*(1.d0 - e2q)*x1*y1/rho2 + gxy(i,j,k) = - psi4*(one - e2q)*x1*y1/rho2 else |