diff options
author | allen <allen@a678b1cf-93e1-4b43-a69d-d43939e66649> | 1999-09-25 10:34:45 +0000 |
---|---|---|
committer | allen <allen@a678b1cf-93e1-4b43-a69d-d43939e66649> | 1999-09-25 10:34:45 +0000 |
commit | 8448c19709c000e089723e276c871563f278b2dd (patch) | |
tree | 9ad226574a910f24367351a24226afb99937a9b8 /src | |
parent | 46976958af59db7a29d02472721dc2d281181e34 (diff) |
This commit was generated by cvs2svn to compensate for changes in r2, which
included commits to RCS files with non-trunk default branches.
git-svn-id: http://svn.einsteintoolkit.org/cactus/EinsteinInitialData/IDBrillData/trunk@3 a678b1cf-93e1-4b43-a69d-d43939e66649
Diffstat (limited to 'src')
-rw-r--r-- | src/brilldata.F | 67 | ||||
-rw-r--r-- | src/brillq.F | 141 | ||||
-rw-r--r-- | src/finishbrilldata.F | 112 | ||||
-rw-r--r-- | src/make.code.defn | 8 | ||||
-rw-r--r-- | src/phif.F | 47 | ||||
-rw-r--r-- | src/setupbrilldata2D.F | 101 | ||||
-rw-r--r-- | src/setupbrilldata3D.F | 167 |
7 files changed, 643 insertions, 0 deletions
diff --git a/src/brilldata.F b/src/brilldata.F new file mode 100644 index 0000000..d7f8348 --- /dev/null +++ b/src/brilldata.F @@ -0,0 +1,67 @@ +#include "cctk.h" +#include "cctk_parameters.h" +#include "cctk_arguments.h" + + subroutine brilldata(CCTK_FARGUMENTS) + +c Author: Carsten Gundlach. +c +c Driver routine for calculating Brill wave initial data. + + implicit none + + DECLARE_CCTK_FARGUMENTS + DECLARE_CCTK_PARAMETERS + DECLARE_CCTK_FUNCTIONS + +c Declare arrays for the linear elliptic solver call: + + integer Mlin_index,Nsrc_index,field_index,ierr + integer metpsi_index(7) + + CCTK_REAL AbsTol(3), RelTol(3) + +c Set up background metric and coefficients for linear solve. + + if (axisym.eq.1) then + call setupbrilldata2D(CCTK_FARGUMENTS) + else + call setupbrilldata3D(CCTK_FARGUMENTS) + end if + +c Call the Linear Elliptic solver interface +c to find conformal factor. + + call CCTK_VarIndex (metpsi_index(1), "einstein::gxx") + call CCTK_VarIndex (metpsi_index(2), "einstein::gxy") + call CCTK_VarIndex (metpsi_index(3), "einstein::gxz") + call CCTK_VarIndex (metpsi_index(4), "einstein::gyy") + call CCTK_VarIndex (metpsi_index(5), "einstein::gyz") + call CCTK_VarIndex (metpsi_index(6), "einstein::gzz") + call CCTK_VarIndex (metpsi_index(7), "einstein::psi") + + call CCTK_VarIndex (field_index, "IDBrillData::brillpsi") + call CCTK_VarIndex (Mlin_index, "IDBrillData::Mlinear") + call CCTK_VarIndex (Nsrc_index, "IDBrillData::Nsource") + + AbsTol(1)= 0.0001 + AbsTol(2)= -1 + AbsTol(3)= -1 + + RelTol(1)= -1 + RelTol(2)= -1 + RelTol(3)= -1 + + call Ell_LinConfMetricSolver(ierr,cctkGH,metpsi_index, + .field_index,Mlin_index,Nsrc_index,AbsTol,RelTol,"sor") + +c Synchronize conformal factor. + + call CCTK_SyncGroup(cctkGH,"einstein::confac") + +c Reconstruct physical metric. + + call finishbrilldata(CCTK_FARGUMENTS) + + return + end diff --git a/src/brillq.F b/src/brillq.F new file mode 100644 index 0000000..fb28639 --- /dev/null +++ b/src/brillq.F @@ -0,0 +1,141 @@ +#include "cctk.h" +#include "cctk_parameters.h" +#include "cctk_arguments.h" + + function brillq(rho,z,phi) + +C Authors: Carsten Gundlach, Miguel Alcubierre. +C +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 + + DECLARE_CCTK_PARAMETERS + + logical firstcall + + integer CCTK_Equals + integer qtype + + real*8 brillq,rho,z,phi + real*8 a,b,c,r0,sr,rho0,srho,brill_sz + real*8 d,e,m,n,phi0 + + data firstcall /.true./ + + save firstcall,qtype,a,b,c,r0,sr,rho0,srho,brill_sz + save d,e,m,n,phi0 + +C Get parameters at first call. + + if (firstcall) then + + qtype = brill_q + + a = brill_a + b = brill_b + c = brill_c + + r0 = brill_r0 + sr = brill_sr + rho0 = brill_rho0 + srho = brill_srho + brill_sz = brill_sz + + if (axisym.eq.0) then + d = brill_d + e = brill_e + m = brill_m + n = brill_n + phi0 = brill_phi0 + end if + + firstcall = .false. + + end if + + if (rho.lt.0) write(*,*) "Warning: negative rho in brillq:", rho + +C Calculate q(rho,z) from a choice of forms. Type 0, 1 and 2 are those +C of Bruegmann. + +C brill_q = 0. + + if (qtype .eq. 0) then + + 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) + end if + + else if (qtype .eq. 1) then + +C brill_q = 1. This includes Eppleys choice of q. +C But note that q(Eppley) = 2q(Cactus). + + brillq = a * (rho/srho)**b + $ / ( 1.d0 + ((rho**2 + z**2 - r0**2) / sr**2)**(c/2) ) + + 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. + + brillq = a * (rho/srho)**b + $ / dexp( ((rho**2 + z**2 - r0**2) / sr**2)**(c/2) ) + + else + +C Unknown type for q function. + + write(*,*) "Unknown type of Brill function q." + STOP + + end if + +C If desired, multiply with 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)) + end if + + return + end + + + diff --git a/src/finishbrilldata.F b/src/finishbrilldata.F new file mode 100644 index 0000000..605abb7 --- /dev/null +++ b/src/finishbrilldata.F @@ -0,0 +1,112 @@ +#include "cctk.h" +#include "cctk_parameters.h" +#include "cctk_arguments.h" +#include "../../../CactusEinstein/Einstein/src/Einstein.h" + + subroutine finishbrilldata(CCTK_FARGUMENTS) + +C Author: Carsten Gundlach. +C +C Set up 3-metric, extrinsic curvature and BM variables +C once the conformal factor has been found. + + implicit none + + DECLARE_CCTK_FARGUMENTS + DECLARE_CCTK_PARAMETERS + + integer i,j,k + integer nx,ny,nz + integer CCTK_Equals + + CCTK_REAL x1,y1,z1,rho1,rho2 + CCTK_REAL brillq,psi4,e2q,rhofudge + CCTK_REAL phi,phif + CCTK_REAL zero + + external brillq + +C Numbers. + + zero = 0.0D0 + +C Set up grid size. + + nx = cctk_lsh(1) + ny = cctk_lsh(2) + nz = cctk_lsh(3) + +C Parameters. + + rhofudge = brill_rhofudge + +C Replace flat metric left over from elliptic solve by +C Brill metric calculated from q and Psi. + + do k=1,nz + do j=1,ny + do i=1,nx + + x1 = x(i,j,k) + y1 = y(i,j,k) + z1 = z(i,j,k) + + rho2 = x1**2 + y1**2 + rho1 = dsqrt(rho2) + + phi = phif(x1,y1) + + psi4 = brillpsi(i,j,k)**4 + e2q = dexp(2.d0*brillq(rho1,z1,phi)) + +C Fudge division by rho^2 on axis. (Physically, y^/rho^2, +C x^2/rho^2 and xy/rho^2 are of course regular.) +C Transform Brills form of the physical metric to Cartesian +C coordinates via +C +C e^2q (drho^2 + dz^2) + rho^2 dphi^2 = +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 + + gzz(i,j,k) = psi4*e2q + gxx(i,j,k) = psi4*(e2q + (1.d0 - e2q)*y1**2/rho2) + gyy(i,j,k) = psi4*(e2q + (1.d0 - e2q)*x1**2/rho2) + gxy(i,j,k) = - psi4*(1.d0 - e2q)*x1*y1/rho2 + + else + +C This fudge assumes that q = O(rho^2) near the axis. Which +C it should be, or the data will be singular. + + gzz(i,j,k) = psi4 + gxx(i,j,k) = psi4 + gyy(i,j,k) = psi4 + gxy(i,j,k) = zero + + end if + + end do + end do + end do + +C In any case, + + gxz = zero + gyz = zero + +C Vanishing extrinsic curvature completes the Cauchy data. + + kxx = zero + kyy = zero + kzz = zero + kxy = zero + kxz = zero + kyz = zero + + return + end + + diff --git a/src/make.code.defn b/src/make.code.defn new file mode 100644 index 0000000..ad8953d --- /dev/null +++ b/src/make.code.defn @@ -0,0 +1,8 @@ +# -*-Makefile-*- +# Sub-makefile for thorn IDBrillData + +SRCS = brilldata.F setupbrilldata2D.F setupbrilldata3D.F \ + brillq.F phif.F finishbrilldata.F + + + diff --git a/src/phif.F b/src/phif.F new file mode 100644 index 0000000..c24f2ca --- /dev/null +++ b/src/phif.F @@ -0,0 +1,47 @@ +#include "cctk.h" + + CCTK_REAL function phif(x,y) + +c Author: Miguel Alcubierre, October 1998. +c +c Function to find angle phi from coordinates {x,y}. + + implicit none + + CCTK_REAL x,y + real*8 zero,half,one,two,pi + +c Define numbers. + + zero = 0.0D0 + half = 0.5D0 + one = 1.0D0 + two = 2.0D0 + + pi = acos(-one) + +c Find angle between 0 and pi/2 such that tan(phi) = |y/x|. + + if (dabs(x).gt.dabs(y)) then + phif = atan(dabs(y/x)) + else if (dabs(x).lt.dabs(y)) then + phif = half*pi - atan(dabs(x/y)) + else + phif = 0.25D0*pi + end if + +c Use signs of {x,y} to move to correct quadrant. + + if ((x.eq.zero).and.(y.eq.zero)) then + phif = zero + else if ((x.le.zero).and.(y.ge.zero)) then + phif = pi - phif + else if ((x.le.zero).and.(y.le.zero)) then + phif = pi + phif + else if ((x.ge.zero).and.(y.le.zero)) then + phif = two*pi - phif + end if + + return + end + diff --git a/src/setupbrilldata2D.F b/src/setupbrilldata2D.F new file mode 100644 index 0000000..2c4b351 --- /dev/null +++ b/src/setupbrilldata2D.F @@ -0,0 +1,101 @@ +#include "cctk.h" +#include "cctk_parameters.h" +#include "cctk_arguments.h" +#include "../../../CactusEinstein/Einstein/src/Einstein.h" + + subroutine setupbrilldata2D(CCTK_FARGUMENTS) + +C Author: Carsten Gundlach. +C +C Set up axisymmetric Brill data for elliptic solve. The elliptic +C equation we need to solve is: +C +C __ 2 2 +C \/ psi = psi ( d q + d q ) / 4 +C f z rho +C +C where: +C +C __ +C \/ = Flat space Laplacian. +C f + + implicit none + + DECLARE_CCTK_FARGUMENTS + DECLARE_CCTK_PARAMETERS + + integer i,j,k + integer nx,ny,nz + + CCTK_REAL x1,y1,z1,rho1 + CCTK_REAL brillq,eps + CCTK_REAL zero,one + + external brillq + +C Numbers. + + zero = 0.0D0 + one = 1.0D0 + +C Set up grid size. + + nx = cctk_lsh(1) + ny = cctk_lsh(2) + nz = cctk_lsh(3) + +C Parameters. + + eps = brill_eps + +C Initialize psi. + + brillpsi = one + +C Set up flat metric for the elliptic solve. +C [Delta(g) + Mlinear] psi = Nsource. + + psi = one + + gxx = one + gyy = one + gzz = one + gxy = zero + gxz = zero + gyz = zero + +C Set up coefficient Mlinear = - (1/4) (q,zz + q,rhorho). + + do k=1,nz + do j=1,ny + do i=1,nx + + x1 = x(i,j,k) + y1 = y(i,j,k) + z1 = z(i,j,k) + + rho1 = dsqrt(x1**2 + y1**2) + +C Centered derivatives. Note that here we may be calling brillq +C with a small negative rho, but that should be ok as long as +C brillq is even in rho - physically it must be, or the data +C will not be regular on the axis. + + Mlinear(i,j,k) = - 0.25d0 + $ *(brillq(rho1,z1+eps,zero) + $ + brillq(rho1,z1-eps,zero) + $ + brillq(rho1+eps,z1,zero) + $ + brillq(rho1-eps,z1,zero) + $ - 4.d0*brillq(rho1,z1,zero))/eps**2 + + end do + end do + end do + +C Set coefficient Nsource = 0. + + Nsource = zero + + return + end diff --git a/src/setupbrilldata3D.F b/src/setupbrilldata3D.F new file mode 100644 index 0000000..719528b --- /dev/null +++ b/src/setupbrilldata3D.F @@ -0,0 +1,167 @@ +#include "cctk.h" +#include "cctk_parameters.h" +#include "cctk_arguments.h" +#include "../../../CactusEinstein/Einstein/src/Einstein.h" + + subroutine setupbrilldata3D(CCTK_FARGUMENTS) + +c Author: Miguel Alcubierre, October 1998. +c +c Set up 3D Brill data for elliptic solve. The elliptic +c equation we need to solve is: +c +c __ +c \/ psi = psi R / 8 +c c c +c +c where: +c __ +c \/ = Laplacian operator for conformal metric. +c c +c +c R = Ricci scalar for conformal metric. +c c +c +c The Ricci scalar for the conformal metric turns out to be: +c +c / -2q 2 2 -2 2 2 \ +c R = 2 | e ( d q + d q ) + rho ( 3 (d q) + 2 d q ) | +c c \ z rho phi phi / + + + implicit none + + DECLARE_CCTK_FARGUMENTS + DECLARE_CCTK_PARAMETERS + + integer CCTK_Equals + integer i,j,k + integer nx,ny,nz + + CCTK_REAL x1,y1,z1,rho1,rho2 + CCTK_REAL phi,phif,e2q + CCTK_REAL brillq,rhofudge,eps + CCTK_REAL zero,one,two,three + +c Set up grid size. + + nx = cctk_lsh(1) + ny = cctk_lsh(2) + nz = cctk_lsh(3) + +c Define numbers + + zero = 0.0D0 + one = 1.0D0 + two = 2.0D0 + three = 3.0D0 + +c Parameters. + + rhofudge = brill_rhofudge + eps = brill_eps + +c Initialize psi. + + brillpsi = one + +c Set up conformal metric. + + psi = one + + do k=1,nz + do j=1,ny + do i=1,nx + + x1 = x(i,j,k) + y1 = y(i,j,k) + z1 = z(i,j,k) + + rho2 = x1**2 + y1**2 + rho1 = dsqrt(rho2) + + phi = phif(x1,y1) + + e2q = dexp(two*brillq(rho1,z1,phi)) + + if (rho1.gt.rhofudge) then + gxx(i,j,k) = e2q + (one - e2q)*y1**2/rho2 + gyy(i,j,k) = e2q + (one - e2q)*x1**2/rho2 + gzz(i,j,k) = e2q + gxy(i,j,k) = - (one - e2q)*x1*y1/rho2 + else + gxx(i,j,k) = one + gyy(i,j,k) = one + gzz(i,j,k) = one + gxy(i,j,k) = zero + end if + + end do + end do + end do + + gxz = zero + gyz = zero + +c Set up coefficient Mlinear = - (1/8) Rc. + + do k=1,nz + do j=1,ny + do i=1,nx + + x1 = x(i,j,k) + y1 = y(i,j,k) + z1 = z(i,j,k) + + rho2 = x1**2 + y1**2 + rho1 = dsqrt(rho2) + + phi = phif(x1,y1) + + e2q = dexp(two*brillq(rho1,z1,phi)) + +c Find M using centered differences over a small +c interval. + + if (rho1.gt.rhofudge) then + Mlinear(i,j,k) = - 0.25D0/e2q + . *(brillq(rho1,z1+eps,phi) + . + brillq(rho1,z1-eps,phi) + . + brillq(rho1+eps,z1,phi) + . + brillq(rho1-eps,z1,phi) + . - 4.0D0*brillq(rho1,z1,phi)) + . / eps**2 + else + Mlinear(i,j,k) = - 0.25D0/e2q + . *(brillq(rho1,z1+eps,phi) + . + brillq(rho1,z1-eps,phi) + . + two*brillq(eps,z1,phi) + . - two*brillq(rho1,z1,phi)) + . / eps**2 + end if + +c Here I assume that for very small rho, the +c phi derivatives are essentially zero. This +c must always be true otherwise the function +c is not regular on the axis. + + if (rho1.gt.rhofudge) then + Mlinear(i,j,k) = Mlinear(i,j,k) - 0.25D0/rho2 + . *(three*0.25D0*(brillq(rho1,z1,phi+eps) + . - brillq(rho1,z1,phi-eps))**2 + . + two*(brillq(rho1,z1,phi+eps) + . - two*brillq(rho1,z1,phi) + . + brillq(rho1,z1,phi-eps)))/eps**2 + end if + + end do + end do + end do + +c Set coefficient Nsource = 0. + + Nsource = zero + + return + end + |