From aa3f2938e6edb094392582654a3771ed22492a77 Mon Sep 17 00:00:00 2001 From: allen Date: Tue, 21 May 2002 21:09:47 +0000 Subject: Homogenise parameters, changed documentation git-svn-id: http://svn.einsteintoolkit.org/cactus/EinsteinInitialData/IDBrillData/trunk@82 a678b1cf-93e1-4b43-a69d-d43939e66649 --- doc/documentation.tex | 154 ++++++++++++++++++++++++++-------------------- param.ccl | 22 +++---- src/brilldata.F | 8 +-- src/finishbrilldata.F | 63 ++++++++----------- src/setupbrilldata2D.F | 2 +- src/setupbrilldata3D.F | 6 +- test/test_brilldata_1.par | 2 +- test/test_brilldata_2.par | 3 +- test/test_wsor.par | 4 +- 9 files changed, 131 insertions(+), 133 deletions(-) diff --git a/doc/documentation.tex b/doc/documentation.tex index a17be7a..2703ebb 100644 --- a/doc/documentation.tex +++ b/doc/documentation.tex @@ -13,9 +13,9 @@ 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 +The purpose of this thorn is to create (time symmetric) 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}, @@ -25,7 +25,29 @@ where $q$ is a free function subject to certain regularity and fall-off conditions, $\rho=\sqrt{x^2+y^2}$ and $\Psi$ is a conformal factor to be solved for. -Substituting this metric into the Hamiltonian constraint gives an +Thorn {\tt IDBrillData} provides three choices for the $q$ function: +an exponential form, ({\tt IDBrillData::q\_function = "exp"}) +\begin{equation} +q = a \; \frac{\rho^{2+b}}{r^2} \left( \frac{z}{\sigma_z} \right)^2 +e^{-(\rho - \rho_0^2)} \left[ 1 + d \frac{\rho^m}{1 + e \rho^m} +\cos^2 \left( n \phi + \phi_0 \right) \right] +\end{equation} +a generalized form of the $q$ function first written down by Eppley +({\tt IDBrillData::q\_function = "eppley"}) +\begin{equation} +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}}\left[ 1 + d \frac{\rho^m}{1 + e \rho^m} +\cos^2 \left( n \phi + \phi_0 \right) \right] +\end{equation} +and the (default) Gundlach $q$ function which includes the Holz form +({\tt IDBrillData::q\_function = "gundlach"}) +\begin{equation} +q = a \left( \frac{\rho}{\sigma_\rho} \right)^b e^{-\left[ +\left( r^2 - r_0^2 \right) / \sigma_r^2 \right]^{c/2}} \left[ 1 + d \frac{\rho^m}{1 + e \rho^m} +\cos^2 \left( n \phi + \phi_0 \right) \right] +\end{equation} + +Substituting the metric into the Hamiltonian constraint gives an elliptic equation for the conformal factor $\Psi$ which is then numerically solved for a given function $q$: \begin{equation} @@ -39,14 +61,9 @@ where the conformal Ricci scalar is found to be Assuming the initial data to be time symmetric means that the momentum constraints are trivially satisfied. -The thorn considers several different forms of the function $q$ -depending on certain parameters that will be described below. - -Brill initial data is activated by choosing the {\tt CactusEinstein/ADMBase} -parameter {\tt initial\_data} to be {\tt brilldata}. - -In the case of axisymmetry, the Hamiltonian constraint can be written -as an elliptic equation for $\Psi$ with just the flat space Laplacian, +In the case of axisymmetry (that is $d=0$ in the above expressions for +$q$), the Hamiltonian constraint can be written as an elliptic +equation for $\Psi$ with just the flat space Laplacian, \begin{equation} \nabla_{flat} \Psi + \frac{\Psi}{4} (\partial_z^2 q + \partial_\rho^2 q) = 0 \end{equation} @@ -55,82 +72,60 @@ ADMBase::initial\_data = "brilldata2D"} then this elliptic equation is solved rather than the equation above. -\section{Parameters for the thorn} +\section{Generating Initial Data with IDBrillData} -The thorn is controlled by the following parameters: +Brill initial data is activated by choosing the {\tt CactusEinstein/ADMBase} +parameter {\tt initial\_data} to be {\tt brilldata}, or for the case of +axisymmetry {\tt brilldata2D} can also be used. -\begin{itemize} +The parameter {\tt IDBrillData::q\_function} chooses the form of the +$q$ function to be used, defaulting to the Gundlach expression. -\item brill\_q (INT): Form of the function $q$ [0,1,2] (default 2): +Additional {\tt IDBrillData} parameters for each form of $q$ fix the +remaining freedom: \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 Exponential $q$: {\tt IDBrillData::q\_function = "exp"} -\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}} -\] +$(a, b,\sigma_z,\rho_0)=$ ({\tt exp\_a, exp\_b, exp\_sigmaz,exp\_rho0}) -\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 Eppley $q$: {\tt IDBrillData::q\_function = "eppley"} -\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] -\] +$(a, b,\sigma_\rho, r_0,\sigma_r,c)=$ ({\tt eppley\_a, eppley\_b, eppley\_sigmarho, eppley\_r0, eppley\_sigmar, eppley\_c}) -\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 Gundlach $q$: {\tt IDBrillData::q\_function = "gundlach"} -\item brill\_e (REAL): $e$ in above expressions (default 1.0). +$(a, b,\sigma_\rho, r_0,\sigma_r,c)=$ ({\tt gundlach\_a, gundlach\_b, gundlach\_sigmarho, gundlach\_r0, gundlach\_sigmar, gundlach\_c}) -\item brill\_m (REAL): $m$ in above expressions (default 2.0). +\item Non-axisymmetric part for each choice of $q$ -\item brill\_n (REAL): $n$ in above expressions (default 2.0). +$(d, m, e, n, \phi0)=$ ({\tt brill3d\_d, brill3d\_m, brill3d\_e, brill3d\_n, brill3d\_phi0}) -\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). - \end{itemize} -The elliptic solver is controlled by the additional parameters: +Note that the default $q$ expression is +$$ +q = {\tt gundlach\_a} \quad \rho^2 e^{-r^2} +$$ + +{\tt IDBrillData} can use the elliptic solvers (type LinMetric) +provided by {\tt CactusEinstein/EllSOR}, {\tt AEIThorns/BAM\_Elliptic} +or {\tt CactusElliptic/EllPETSc} to solve the equation resulting from +the Hamiltonian constraint. +In all cases the parameter {\tt thresh} sets the threshold for the elliptic +solve. The choice of elliptic solver is made +through the parameter {\tt brill\_solver}: \begin{itemize} -\item {\tt solver} (KEYWORD): Elliptic solver used to solve the - hamiltonian constraint [sor/petsc/bam] (default "sor"). +\item {\tt sor}: Understands the Robin boundary condition, additional +parameters control the maximum number of iterations ({\tt sor\_maxit}). -\item {\tt thresh} (REAL): Threshold for elliptic solver (default - 0.00001). +\item {\tt bam}: {\tt BAM\_Elliptic} does not properly implement the +elliptic infrastructure of {\tt EllBase}, and the {\tt BAM\_Elliptic} +parameter to use the Robin boundary condition must be set independently +of {tt IDBrillWave::brill\_bound}. \end{itemize} @@ -143,6 +138,29 @@ the derivatives of the conformal factor are not calculated, so that only {\tt staticconformal::conformal\_storage = "factor"} is supported. +\section{References} + +\subsection{Specification of Brill Waves} + +\begin{enumerate} + +\item Dieter Brill, {\bf Ann. Phys.}, 7, 466, 1959. + +\item Ken Eppley, {\bf Sources of Gravitational Radiation}, edited by L. Smarr (Cambridge University Press, +Cambridge, England, 1979), p. 275. + +\end{enumerate} + +\subsection{Numerical Evolutions of Brill Waves} + +\begin{enumerate} + +\item {\it Gravitational Collapse of Gravitational Waves in 3D Numerical Relativity}, + Miguel Alcubierre, Gabrielle Allen, Bernd Bruegmann, Gerd Lanfermann, Edward Seidel, Wai-Mo Suen, Malcolm Tobias, +{\bf Phys. Rev. D61}, 041501, 2000. + +\end{enumerate} + % Automatically created from the ccl files % Do not worry for now. diff --git a/param.ccl b/param.ccl index 1313ac5..4586bab 100644 --- a/param.ccl +++ b/param.ccl @@ -20,35 +20,35 @@ private: # Parameters for elliptic solve -KEYWORD brill_solver "Which elliptic solver to use" +KEYWORD solver "Which elliptic solver to use" { "sor" :: "Use SOR solver" "petsc" :: "Use PETSc solver" "bam" :: "Use bam solver" } "sor" -KEYWORD brill_bound "Which boundary condition to use" +KEYWORD bound "Which boundary condition to use" { - "const" :: "constant boundary: set brill_const_v0" - "robin" :: "Robin boundary: set brill_robin_falloff, brill_robin_inf" -} "const" + "const" :: "constant boundary: set const_v0" + "robin" :: "Robin boundary: set robin_falloff, robin_inf" +} "robin" -INT brill_robin_falloff "Fall-off of Robin BC" +INT robin_falloff "Fall-off of Robin BC" { 0: :: "any positive integer value" } 1 -REAL brill_const_v0 "Value of constant BC" +REAL const_v0 "Value of constant BC" { : :: "anything goes" } 1.0 -REAL brill_robin_inf "Value at infinity of Robin BC" +REAL robin_inf "Value at infinity of Robin BC" { : :: "anything goes" } 1.0 -REAL brill_thresh "How far (absolute norm) to go" +REAL thresh "How far (absolute norm) to go" { 0.0: :: "Positive number please" } 0.00001 @@ -181,7 +181,7 @@ REAL brill3d_phi0 "3D Brill wave: d rho^m cos^2(n (phi + phi0))" # Additional parameters -REAL brill_rhofudge "delta rho for axis fudge" +REAL rhofudge "delta rho for axis fudge" { 0: :: "Positive please" } 0.00001 @@ -189,7 +189,7 @@ REAL brill_rhofudge "delta rho for axis fudge" INT sor_maxit "Maximum number of iterations" { 0:* :: "Positive" -} 100 +} 10000 BOOLEAN output_coeffs "output coefficients for elliptic solve" { diff --git a/src/brilldata.F b/src/brilldata.F index 3c91b97..e38afaa 100644 --- a/src/brilldata.F +++ b/src/brilldata.F @@ -75,7 +75,7 @@ c Set up background metric and coefficients for linear solve. c Tolerances for elliptic solve. - AbsTol(1)= brill_thresh + AbsTol(1)= thresh AbsTol(2)= -1.0D0 AbsTol(3)= -1.0D0 @@ -94,14 +94,14 @@ c Elliptic solver c Just in case some solvers do not use the standard interface conformal_state = 0 - if (CCTK_EQUALS(brill_solver,"sor")) then + if (CCTK_EQUALS(solver,"sor")) then call Ell_SetIntKey(ierr,sor_maxit,"Ell::SORmaxit") call Ell_LinMetricSolver(ierr,cctkGH, . metric_index,ipsi,iMcoeff,iNcoeff,AbsTol,RelTol,"sor") - else if (CCTK_EQUALS(brill_solver,"petsc")) then + else if (CCTK_EQUALS(solver,"petsc")) then call Ell_LinMetricSolver(ierr,cctkGH, . metric_index,ipsi,iMcoeff,iNcoeff,AbsTol,RelTol,"petsc") - else if (CCTK_EQUALS(brill_solver,"bam")) then + else if (CCTK_EQUALS(solver,"bam")) then call Ell_LinMetricSolver(ierr,cctkGH, . metric_index,ipsi,iMcoeff,iNcoeff,AbsTol,RelTol,"bam") end if diff --git a/src/finishbrilldata.F b/src/finishbrilldata.F index 467c28d..bd13c79 100644 --- a/src/finishbrilldata.F +++ b/src/finishbrilldata.F @@ -21,10 +21,9 @@ DECLARE_CCTK_FUNCTIONS integer i,j,k - integer nx,ny,nz CCTK_REAL x1,y1,z1,rho1,rho2 - CCTK_REAL phi,psi4,e2q,rhofudge + CCTK_REAL phi,psi4,e2q CCTK_REAL zero,one CCTK_REAL brillq,phif @@ -34,22 +33,12 @@ 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. - - 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 + do k=1,cctk_lsh(3) + do j=1,cctk_lsh(2) + do i=1,cctk_lsh(1) x1 = x(i,j,k) y1 = y(i,j,k) @@ -84,13 +73,23 @@ c The individual coefficients can be read off as c This fudge assumes that q = O(rho^2) near the axis. Which c it should be, or the data will be singular. - gxx(i,j,k) = 0.0d0 - gyy(i,j,k) = 0.0d0 - gzz(i,j,k) = 0.0d0 - gxy(i,j,k) = 0.0d0 + gxx(i,j,k) = zero + gyy(i,j,k) = zero + gzz(i,j,k) = zero + gxy(i,j,k) = zero end if + gxz(i,j,k) = zero + gyz(i,j,k) = zero + + kxx(i,j,k) = zero + kyy(i,j,k) = zero + kzz(i,j,k) = zero + kxy(i,j,k) = zero + kxz(i,j,k) = zero + kyz(i,j,k) = zero + end do end do end do @@ -99,9 +98,9 @@ c it should be, or the data will be singular. conformal_state = 1 - do k=1,nz - do j=1,ny - do i=1,nx + do k=1,cctk_lsh(3) + do j=1,cctk_lsh(2) + do i=1,cctk_lsh(1) psi(i,j,k) = brillpsi(i,j,k) @@ -113,9 +112,9 @@ c it should be, or the data will be singular. conformal_state = 0 - do k=1,nz - do j=1,ny - do i=1,nx + do k=1,cctk_lsh(3) + do j=1,cctk_lsh(2) + do i=1,cctk_lsh(1) psi4 = brillpsi(i,j,k)**4 @@ -130,19 +129,5 @@ c it should be, or the data will be singular. end if -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/setupbrilldata2D.F b/src/setupbrilldata2D.F index f2e5f13..185bef5 100644 --- a/src/setupbrilldata2D.F +++ b/src/setupbrilldata2D.F @@ -103,7 +103,7 @@ c will not be regular on the axis. c Synchronization and boundaries. c call CCTK_SyncGroup(ierr,cctkGH,'idbrilldata::brillelliptic') - call CartSymGN(ierr,cctkGH,'idbrilldata::brillelliptic') +c call CartSymGN(ierr,cctkGH,'idbrilldata::brillelliptic') return end diff --git a/src/setupbrilldata3D.F b/src/setupbrilldata3D.F index 97cf28b..b77d886 100644 --- a/src/setupbrilldata3D.F +++ b/src/setupbrilldata3D.F @@ -47,7 +47,7 @@ c c \ z rho phi phi / CCTK_REAL x1,y1,z1,rho1,rho2 CCTK_REAL phi,phif,e2q - CCTK_REAL brillq,rhofudge,eps + CCTK_REAL brillq,eps CCTK_REAL zero,one,two,three c Define numbers @@ -57,10 +57,6 @@ c Define numbers two = 2.0D0 three = 3.0D0 -c Parameters. - - rhofudge = brill_rhofudge - c Epsilon for finite differencing. eps = cctk_delta_space(1) diff --git a/test/test_brilldata_1.par b/test/test_brilldata_1.par index 5d673ba..0d15a58 100644 --- a/test/test_brilldata_1.par +++ b/test/test_brilldata_1.par @@ -27,7 +27,7 @@ idbrilldata::gundlach_a = 5.0 # Elliptic solver. -idbrilldata::brill_solver = "bam" +idbrilldata::solver = "bam" bam_elliptic::bam_bound = "bamrobin" # Output. diff --git a/test/test_brilldata_2.par b/test/test_brilldata_2.par index 630cd64..4e0b101 100644 --- a/test/test_brilldata_2.par +++ b/test/test_brilldata_2.par @@ -28,8 +28,7 @@ idbrilldata::gundlach_a = 5.0 # Elliptic solver. -idbrilldata::brill_solver = "bam" - +idbrilldata::solver = "bam" bam_elliptic::bam_bound = "bamrobin" # Output. diff --git a/test/test_wsor.par b/test/test_wsor.par index a9b0337..4c04382 100644 --- a/test/test_wsor.par +++ b/test/test_wsor.par @@ -26,8 +26,8 @@ idbrilldata::gundlach_a = 5.0 # Elliptic solver. -idbrilldata::brill_solver = "sor" -idbrilldata::brill_bound = "robin" +idbrilldata::solver = "sor" +idbrilldata::bound = "robin" idbrilldata::sor_maxit = 100 # Warning -- since this doesn't converge for so few iterations, the # result will depend on number of cpus, even beyond the sor tolerance -- cgit v1.2.3