From 8448c19709c000e089723e276c871563f278b2dd Mon Sep 17 00:00:00 2001 From: allen Date: Sat, 25 Sep 1999 10:34:45 +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/IDBrillData/trunk@3 a678b1cf-93e1-4b43-a69d-d43939e66649 --- README | 9 +++ doc/documentation.tex | 142 +++++++++++++++++++++++++++++++++++++++++ interface.ccl | 17 +++++ par/brilldata.par | 50 +++++++++++++++ param.ccl | 119 +++++++++++++++++++++++++++++++++++ schedule.ccl | 12 ++++ src/brilldata.F | 67 ++++++++++++++++++++ src/brillq.F | 141 +++++++++++++++++++++++++++++++++++++++++ src/finishbrilldata.F | 112 +++++++++++++++++++++++++++++++++ src/make.code.defn | 8 +++ src/phif.F | 47 ++++++++++++++ src/setupbrilldata2D.F | 101 ++++++++++++++++++++++++++++++ src/setupbrilldata3D.F | 167 +++++++++++++++++++++++++++++++++++++++++++++++++ 13 files changed, 992 insertions(+) create mode 100644 README create mode 100644 doc/documentation.tex create mode 100644 interface.ccl create mode 100644 par/brilldata.par create mode 100644 param.ccl create mode 100644 schedule.ccl create mode 100644 src/brilldata.F create mode 100644 src/brillq.F create mode 100644 src/finishbrilldata.F create mode 100644 src/make.code.defn create mode 100644 src/phif.F create mode 100644 src/setupbrilldata2D.F create mode 100644 src/setupbrilldata3D.F diff --git a/README b/README new file mode 100644 index 0000000..a6b13a9 --- /dev/null +++ b/README @@ -0,0 +1,9 @@ +Cactus Code Thorn BrillData +Authors : Carsten Gundlach, Miguel Alcubierre +Managed by : C.G. + M.A. +-------------------------------------------------------------------------- + +Initializes Brill wave data. + + diff --git a/doc/documentation.tex b/doc/documentation.tex new file mode 100644 index 0000000..59a7721 --- /dev/null +++ b/doc/documentation.tex @@ -0,0 +1,142 @@ +% Thorn documentation template +\documentclass{article} +\begin{document} + +\title{IDBrillData} +\author{Carsten Gundlach} +\date{6 September 1999} +\maketitle + +\abstract{This thorn creates initial data for Brill wave spacetimes. +It can create both axisymmetric data (in a 3D cartesian grid), as +well as data with an angular dependency.} + +\section{Purpose} + +The purpose of this thorn is to create initial data for a Brill wave +spacetime. It does so by starting from a three--metric of the form +originally considered by Brill +\begin{equation} +ds^2 = \Psi^4 \left[ e^{2q} \left( d\rho^2 + dz^2 \right) ++ \rho^2 d\phi^2 \right] =\Psi^4 \hat{ds}^{2}, +\label{eqn:brillmetric} +\end{equation} +where $q$ is a free function subject to certain regularity and +fall-off conditions and $Psi$ is a conformal factor to be solved for. + +The thorn considers several different forms of the function $q$ +depending on certain parameters that will be described below. +Substituting the metric above into the Hamiltonian constraint results +in an elliptic equation for the conformal factor $Psi$ that can be +solved numerically once the function $q$ has been specified. The +initial data is also assumed to be time-symmetric, so the momentum +constraints are trivially satisfied. + +The thorn is activated by choosing the standard Cactus parameter +``initial\_data'' in one of the following two ways: + +\begin{itemize} + +\item initial\_data = ``brilldata'': Axisymmetric Brill wave initial data + (but calculated in a cartesian grid!). + +\item initial\_data = ``brilldata3D'': Brill wave initial data with an + angular dependency. + +\end{itemize} + + +\section{Parameters for the thorn} + +The thorn is controlled by the following parameters: + +\begin{itemize} + +\item brill\_q (INT): Form of the function $q$ [0,1,2] (default 2): + +\begin{itemize} + +\item brill\_q = 0: +\[ +q = a \; \frac{\rho^{2+b}}{r^2} \left( \frac{z}{\sigma_z} \right)^2 +e^{-(\rho - \rho_0^2)} +\] + +\item brill\_q = 1: +\[ +q = a \left( \frac{\rho}{\sigma_\rho} \right)^b \frac{1}{1 + \left[ +\left( r^2 - r_0^2 \right) / \sigma_r^2 \right]^{c/2}} +\] + +\item brill\_q = 2: +\[ +q = a \left( \frac{\rho}{\sigma_\rho} \right)^b e^{-\left[ +\left( r^2 - r_0^2 \right) / \sigma_r^2 \right]^{c/2}} +\] + +\item If one specifies 3D data (see above), the function $q$ is multiplied +by an additional factor with an angular dependency: +\[ +q \rightarrow q \left[ 1 + d \frac{\rho^m}{1 + e \rho^m} +\cos^2 \left( n \phi + \phi_0 \right) \right] +\] + +\end{itemize} + +\item brill\_a (REAL): Amplitude (default 0.0). + +\item brill\_b (REAL): $b$ in above expressions (default 2.0). + +\item brill\_c (REAL): $c$ in above expressions (default 2.0). + +\item brill\_d (REAL): $d$ in above expressions (default 0.0). + +\item brill\_e (REAL): $e$ in above expressions (default 1.0). + +\item brill\_m (REAL): $m$ in above expressions (default 2.0). + +\item brill\_n (REAL): $n$ in above expressions (default 2.0). + +\item brill\_r0 (REAL): $r_0$ in above expressions (default 0.0). + +\item brill\_rho0 (REAL): $\rho_0$ in above expressions + (default 0.0). + +\item brill\_phi0 (REAL): $\phi_0$ in above expressions + (default 0.0). + +\item brill\_sr (REAL): $\sigma_r$ in above expressions + (default 1.0). + +\item brill\_srho (REAL): $\sigma_\rho$ in above + expressions (default 1.0). + +\item savepsi (KEYWORD): Save conformal factor for output? + [``yes'',''no''] Normally, the conformal factor is calculated in the + grid function ``psi'', but this is set back to 1 at the end once the + physical metric has been constructed. Setting this parameter to + ``yes'' copies the conformal factor to the grid function + ``brillpsi0'' before resetting it to 1 (default ``no''). + +\end{itemize} + +The elliptic solver is controlled by additional the parameters: + +\begin{itemize} + +\item solver (KEYWORD): Elliptic solver used to solve the + hamiltonian constraint [sor/petsc/bam] (default ``sor''). + +\item thresh (REAL): Threshold for elliptic solver (default + 0.00001). + +\end{itemize} + +% Automatically created from the ccl files +% Do not worry for now. + +\include{interface} +\include{param} +\include{schedule} + +\end{document} diff --git a/interface.ccl b/interface.ccl new file mode 100644 index 0000000..f7b01c8 --- /dev/null +++ b/interface.ccl @@ -0,0 +1,17 @@ +implements: IDBrillData +inherits: einstein + +# Grid functions for linear elliptic equation. + +private: + +real linearelliptic TYPE=GF +{ + Mlinear, + Nsource +} "Coefficients for linear elliptic equation" + +real brillconf TYPE=GF +{ + brillpsi +} "Conformal factor for Brill data" diff --git a/par/brilldata.par b/par/brilldata.par new file mode 100644 index 0000000..52a4575 --- /dev/null +++ b/par/brilldata.par @@ -0,0 +1,50 @@ +################################################################# +# +# brilldata.par +# +# Simply Brill wave data of the Hoz et.al. type +# +################################################################# + +ActiveThorns = "iobasic ellbase ellsor IDBrillData pugh einstein cartgrid3d ioutil ioascii" + +# Brill wave initial data + +Einstein::initial_data = "brilldata" + +# Brill wave parameters [Holz et.al. type data] + +idbrilldata::brill_q = 2 # Form of function q [0,1,2] +idbrilldata::brill_a = 1.0 # Amplitude (default 0) + +idbrilldata::savepsi = "yes" + +# Elliptic solver. + +idbrilldata::solver = "sor" #[petsc,sor,bam] + +# General + +driver::global_nx = 32 +driver::global_ny = 32 +driver::global_nz = 32 + +grid::type = "byspacing" +grid::dxyz = 0.2 +grid::domain = "octant" + +cactus::cctk_itlast = 0.0 + +# Output. + +IO::outdir = "BrillData" + +IO::out_every = 1 + +IOASCII::out1D_vars = "einstein::gxx einstein::gyy einstein::gzz + einstein::gxy einstein::gxz einstein::gyz + IDBrillData: brillpsi0" + +IOBasic::outScalar_vars = "einstein::gxx einstein::gyy einstein::gzz + einstein::gxy einstein::gxz einstein::gyz + IDBrillData: brillpsi0" diff --git a/param.ccl b/param.ccl new file mode 100644 index 0000000..e624a2c --- /dev/null +++ b/param.ccl @@ -0,0 +1,119 @@ +# Parameter definitions for thorn IDBrillData + +shares:einstein + +EXTENDS KEYWORD initial_data "" +{ + "brilldata" :: "Brill wave initial data" +} "" + + +private: + +# Parameters for elliptic solve + +KEYWORD solver "Which elliptic solver to use" +{ + "sor" :: "Use SOR solver" + "petsc" :: "Use PETSc solver" + "bam" :: "Use bam solver" +} "sor" + +REAL thresh "How far (absolute norm) to go" +{ + 0.0: :: "Positive number please" +} 0.00001 + + +# Brill wave parameters + +BOOLEAN axisym "Axisymmetric initial data?" +{ +} "yes" + +INT brill_q "Form of function q [0,1,2]" +{ + 0:2 :: "Only cases 0,1,2 defined at the moment" +} 1 + +REAL brill_a "Brill wave: Amplitude" +{ + : :: +} 0.0 + +REAL brill_b "Brill wave: rho^b" +{ + : :: +} 2.0 + +REAL brill_c "Brill wave: (r^2 - r0^2)^(c/2)" +{ + : :: +} 2.0 + +REAL brill_rho0 "Brill wave: radius of torus in rho" +{ + : :: +} 0.0 + +REAL brill_r0 "Brill wave: radius of torus in r" +{ + : :: +} 0.0 + +REAL brill_srho "Brill wave: sigma in rho" +{ + : :: +} 1.0 + +REAL brill_sr "Brill wave: sigma in r" +{ + : :: +} 1.0 + +# 3D Brill wave parameters + +REAL brill_d "3D Brill wave: d rho^m cos^2(n (phi + phi0))" +{ + : :: +} 0.0 + +REAL brill_e "3D Brill wave: d rho^m cos^2(n (phi + phi0))" +{ + : :: +} 1.0 + +REAL brill_m "3D Brill wave: d rho^m cos^2(n (phi + phi0))" +{ + : :: +} 2.0 + +REAL brill_n "3D Brill wave: d rho^m cos^2(n (phi + phi0))" +{ + : :: +} 2.0 + +REAL brill_phi0 "3D Brill wave: d rho^m cos^2(n (phi + phi0))" +{ + : :: +} 0.0 + +KEYWORD outputpsi "Output conformal factor?" +{ + "yes" :: "Output brillpsi" + "no" :: "Don't output brillpsi" +}"yes" + + +# 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 + diff --git a/schedule.ccl b/schedule.ccl new file mode 100644 index 0000000..94bd666 --- /dev/null +++ b/schedule.ccl @@ -0,0 +1,12 @@ + +if (CCTK_Equals(initial_data,"brilldata")) +{ + schedule brilldata at CCTK_INITIAL + { + STORAGE: linearelliptic,brillconf + COMMUNICATION: linearelliptic,brillconf + LANG: Fortran + } "Construct Brill wave initial data" +} + + 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 + -- cgit v1.2.3