diff options
| -rw-r--r-- | doc/documentation.tex | 79 | ||||
| -rw-r--r-- | param.ccl | 91 | ||||
| -rw-r--r-- | schedule.ccl | 2 | ||||
| -rw-r--r-- | src/ParamCheck.c | 31 | ||||
| -rw-r--r-- | src/Startup.c | 4 | ||||
| -rw-r--r-- | src/brilldata.F | 30 | ||||
| -rw-r--r-- | src/brillq.F | 69 | ||||
| -rw-r--r-- | src/finishbrilldata.F | 6 | ||||
| -rw-r--r-- | src/setupbrilldata2D.F | 58 | ||||
| -rw-r--r-- | src/setupbrilldata3D.F | 118 | ||||
| -rw-r--r-- | test/test_brilldata_1.par | 6 | ||||
| -rw-r--r-- | test/test_brilldata_2.par | 6 | ||||
| -rw-r--r-- | test/test_wsor.par | 9 |
13 files changed, 269 insertions, 240 deletions
diff --git a/doc/documentation.tex b/doc/documentation.tex index 49731f9..a17be7a 100644 --- a/doc/documentation.tex +++ b/doc/documentation.tex @@ -3,13 +3,13 @@ \begin{document} \title{IDBrillData} -\author{Carsten Gundlach} -\date{6 September 1999} +\author{Carsten Gundlach, Gabrielle Allen} +\date{$Date$} \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.} +\abstract{This thorn creates time symmetric 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} @@ -22,40 +22,38 @@ ds^2 = \Psi^4 \left[ e^{2q} \left( d\rho^2 + dz^2 \right) \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. +fall-off conditions, $\rho=\sqrt{x^2+y^2}$ 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 +Substituting this 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} +\hat{\nabla} \Psi - \frac{\Psi}{8} \hat{R} = 0 +\end{equation} +where the conformal Ricci scalar is found to be +\begin{eqnarray} +\hat{R} = -2 \left(e^{-2q} (\partial^2_z q + \partial^2_\rho q) + +\frac{1}{\rho^2} (3 (\partial_\phi q)^2 + 2 \partial_\phi q)\right) +\end{eqnarray} +Assuming the initial data to be time symmetric means that the momentum constraints are trivially satisfied. -% [[ DPR: The code does not use these two different initial_data -% keywords, but the axisym parameter instead. It should probably use -% the two different initial_data values. I have changed the doc below -% to correspond with the current code: ]] +The thorn considers several different forms of the function $q$ +depending on certain parameters that will be described below. -The thorn is activated by choosing the CactusEinstein/ADMBase parameter -``initial\_data'' to be: -%in one of the following two ways: +Brill initial data is activated by choosing the {\tt CactusEinstein/ADMBase} +parameter {\tt initial\_data} to be {\tt brilldata}. -\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} +In the case of axisymmetry, 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} +If the initial data is chosen to be {\tt +ADMBase::initial\_data = "brilldata2D"} then this elliptic equation +is solved rather than the equation above. -To choose axisymmetric data, set the parameter {\tt axisym} to true -(the default). Note that the data is computed on a cartesian grid in -any event. {\tt axisym = "no"} adds an angular dependence factor, -which is detailed below. \section{Parameters for the thorn} @@ -128,14 +126,23 @@ The elliptic solver is controlled by the additional parameters: \begin{itemize} -\item solver (KEYWORD): Elliptic solver used to solve the - hamiltonian constraint [sor/petsc/bam] (default ``sor''). +\item {\tt solver} (KEYWORD): Elliptic solver used to solve the + hamiltonian constraint [sor/petsc/bam] (default "sor"). -\item thresh (REAL): Threshold for elliptic solver (default +\item {\tt thresh} (REAL): Threshold for elliptic solver (default 0.00001). \end{itemize} +\section{Notes} + +Thorn {\tt IDBrillData} understands both the {\tt physical} and {\tt +static conformal} {\tt metric\_type}. In the case of a conformal +metric being chosen, the conformal factor is set to $\Psi$. Currently +the derivatives of the conformal factor are not calculated, so that +only {\tt staticconformal::conformal\_storage = "factor"} is +supported. + % Automatically created from the ccl files % Do not worry for now. @@ -6,6 +6,7 @@ shares: ADMBase EXTENDS KEYWORD initial_data { "brilldata" :: "Brill wave initial data" + "brilldata2D" :: "Brill wave initial data assuming axisymmetry" } USES KEYWORD metric_type @@ -55,78 +56,124 @@ REAL brill_thresh "How far (absolute norm) to go" # Brill wave parameters -BOOLEAN axisym "Axisymmetric initial data?" +KEYWORD q_function "Form of function q [0,1,2]" { -} "yes" + "exp" :: "" + "eppley" :: + "gundlach" :: +} "gundlach" -INT brill_q "Form of function q [0,1,2]" +# ------------------------------------------------------------- + +REAL exp_a "Exp Brill wave: Amplitude" +{ + : :: "Anything" +} 0.0 + +INT exp_b "Exp Brill wave: used in exponent in rho: rho^(2+b)" { - 0:2 :: "Only cases 0,1,2 defined at the moment" + : :: "Anything" } 2 -REAL brill_a "Brill wave: Amplitude" +REAL exp_rho0 "Exp Brill wave: radius of torus in rho" +{ + 0:* :: "Positive" +} 0.0 + +REAL exp_sigmaz "Exp Brill wave: sigma in z" +{ + (0:* :: "Positive" +} 1.0 + +# ------------------------------------------------------------- + +REAL eppley_a "Eppley Brill wave: Amplitude" { : :: "Anything" } 0.0 -INT brill_b "Brill wave: used in exponent in rho" +INT eppley_b "Eppley Brill wave: used in exponent in rho: rho^b" { : :: "Anything" } 2 -INT brill_c "Brill wave: (r^2 - r0^2)^(c/2)" +INT eppley_c "Eppley Brill wave: (r^2 - r0^2)^(c/2)" { : :: "Anything" } 2 -REAL brill_rho0 "Brill wave: radius of torus in rho" +REAL eppley_r0 "Eppley Brill wave: radius of torus in r" { - : :: "Anything" + (0:* :: "Positive" } 0.0 -REAL brill_r0 "Brill wave: radius of torus in r" +REAL eppley_sigmarho "Eppley Brill wave: sigma in rho" { - : :: "Anything" -} 0.0 + : :: "Anything" +} 1.0 -REAL brill_srho "Brill wave: sigma in rho" +REAL eppley_sigmar "Eppley Brill wave: sigma in r" { : :: "Anything" } 1.0 -REAL brill_sr "Brill wave: sigma in r" +# ------------------------------------------------------------- + +REAL gundlach_a "Gundlach Brill wave: Amplitude" +{ + : :: "Anything" +} 0.0 + +INT gundlach_b "Gundlach Brill wave: used in exponent in rho: rho^b" +{ + : :: "Anything" +} 2 + +REAL gundlach_sigmarho "Gundlach Brill wave: sigma in rho" { : :: "Anything" } 1.0 -REAL brill_sz "Brill wave: sigma in z" +REAL gundlach_r0 "Gundlach Brill wave: radius of torus in r" +{ + (0:* :: "Positive" +} 0.0 + +INT gundlach_c "Gundlach Brill wave: (r^2 - r0^2)^(c/2)" +{ + : :: "Anything" +} 2 + +REAL gundlach_sigmar "Gundlach Brill wave: sigma in r" { : :: "Anything" } 1.0 +# ------------------------------------------------------------- + # 3D Brill wave parameters -REAL brill_d "3D Brill wave: d rho^m cos^2(n (phi + phi0))" +REAL brill3d_d "3D Brill wave: d rho^m cos^2(n (phi + phi0))" { : :: "Anything" } 0.0 -REAL brill_e "3D Brill wave: d rho^m cos^2(n (phi + phi0))" +REAL brill3d_e "3D Brill wave: d rho^m cos^2(n (phi + phi0))" { : :: "Anything" } 1.0 -REAL brill_m "3D Brill wave: d rho^m cos^2(n (phi + phi0))" +REAL brill3d_m "3D Brill wave: d rho^m cos^2(n (phi + phi0))" { : :: "Anything" } 2.0 -REAL brill_n "3D Brill wave: d rho^m cos^2(n (phi + phi0))" +REAL brill3d_n "3D Brill wave: d rho^m cos^2(n (phi + phi0))" { : :: "Anything" } 2.0 -REAL brill_phi0 "3D Brill wave: d rho^m cos^2(n (phi + phi0))" +REAL brill3d_phi0 "3D Brill wave: d rho^m cos^2(n (phi + phi0))" { : :: "Anything" } 0.0 @@ -143,3 +190,7 @@ INT sor_maxit "Maximum number of iterations" { 0:* :: "Positive" } 100 + +BOOLEAN output_coeffs "output coefficients for elliptic solve" +{ +} "no"
\ No newline at end of file diff --git a/schedule.ccl b/schedule.ccl index 4d06086..0e8035e 100644 --- a/schedule.ccl +++ b/schedule.ccl @@ -1,7 +1,7 @@ # Schedule definitions for thorn IDBrillData # $Header$ -if (CCTK_Equals(initial_data,"brilldata")) +if (CCTK_Equals(initial_data,"brilldata") || CCTK_Equals(initial_data,"brilldata2D")) { STORAGE: brillconf diff --git a/src/ParamCheck.c b/src/ParamCheck.c index fe55173..7a04cec 100644 --- a/src/ParamCheck.c +++ b/src/ParamCheck.c @@ -52,7 +52,9 @@ void IDBrillData_ParamChecker(CCTK_ARGUMENTS) if( ! CCTK_Equals(metric_type, "physical") && ! CCTK_Equals(metric_type, "static conformal")) { - CCTK_PARAMWARN("Unknown ADMBase::metric_type - known types are \"physical\" and \"static conformal\""); + CCTK_VParamWarn("Unknown ADMBase::metric_type (%s): " + "known types are \"physical\" and \"static conformal\"", + metric_type); } if (CCTK_Equals(metric_type, "static conformal")) @@ -72,6 +74,33 @@ void IDBrillData_ParamChecker(CCTK_ARGUMENTS) { CCTK_INFO(" ... using physical metric"); } + + if (CCTK_Equals(initial_data,"brilldata")) + { + CCTK_INFO(" ... using full 3D solver"); + } + else if (CCTK_Equals(initial_data,"brilldata2D")) + { + CCTK_INFO(" ... using reduced 2D solver"); + if (brill3d_d != 0.0) + { + CCTK_WARN(0,"But 3D parameter brill3d_d set for the 2D solver !?!"); + } + } + + if (CCTK_Equals(q_function,"exp")) + { + CCTK_INFO(" ... Exponential brill function"); + } + else if (CCTK_Equals(q_function,"eppley")) + { + CCTK_INFO(" ... generalised Eppley Brill function"); + } + else if (CCTK_Equals(q_function,"exp")) + { + CCTK_INFO(" ... Gundlach-Holz Brill function"); + } + } /******************************************************************** diff --git a/src/Startup.c b/src/Startup.c index 8660cbe..a6d4c13 100644 --- a/src/Startup.c +++ b/src/Startup.c @@ -18,7 +18,9 @@ void BrillData_InitSymBound(CCTK_ARGUMENTS) sym[1] = 1; sym[2] = 1; - SetCartSymVN(cctkGH, sym,"IDBrillData::brillpsi"); + SetCartSymVN(cctkGH,sym,"IDBrillData::brillpsi"); + SetCartSymVN(cctkGH,sym,"IDBrillData::brillMlinear"); + SetCartSymVN(cctkGH,sym,"IDBrillData::brillNsource"); return; } diff --git a/src/brilldata.F b/src/brilldata.F index 7dc01e0..3c91b97 100644 --- a/src/brilldata.F +++ b/src/brilldata.F @@ -46,22 +46,33 @@ c Get indices for grid functions. call CCTK_VarIndex(iMcoeff,"idbrilldata::brillMlinear") if (iMcoeff.lt.0) then - call CCTK_WARN(0,"Grid variable index for Mcoeff not found") + call CCTK_WARN(0,"Grid variable index for brillMlinear not found") end if call CCTK_VarIndex(iNcoeff,"idbrilldata::brillNsource") if (iNcoeff.lt.0) then - call CCTK_WARN(0,"Grid variable index for Ncoeff not found") + call CCTK_WARN(0,"Grid variable index for brillNsource not found") end if c Set up background metric and coefficients for linear solve. - if (axisym.eq.1) then + if (CCTK_EQUALS(initial_data,"brilldata2D")) then call setupbrilldata2D(CCTK_ARGUMENTS) else call setupbrilldata3D(CCTK_ARGUMENTS) end if + if (output_coeffs .eq. 1) then + call CCTK_OutputVarAs(ierr,cctkGH,"IDBrillData::brillMlinear","Debug_brillMlinear") + call CCTK_OutputVarAs(ierr,cctkGH,"IDBrillData::brillNsource","Debug_brillNsource") + call CCTK_OutputVarAs(ierr,cctkGH,"ADMBase::gxx","Debug_gxx") + call CCTK_OutputVarAs(ierr,cctkGH,"ADMBase::gxy","Debug_gxy") + call CCTK_OutputVarAs(ierr,cctkGH,"ADMBase::gxz","Debug_gxz") + call CCTK_OutputVarAs(ierr,cctkGH,"ADMBase::gyy","Debug_gyy") + call CCTK_OutputVarAs(ierr,cctkGH,"ADMBase::gyz","Debug_gyz") + call CCTK_OutputVarAs(ierr,cctkGH,"ADMBase::gzz","Debug_gzz") + end if + c Tolerances for elliptic solve. AbsTol(1)= brill_thresh @@ -74,14 +85,17 @@ c Tolerances for elliptic solve. c Boundaries. - call Ell_SetStrKey(ierr,"yes","EllLinConfMetric::Bnd::Robin") - call Ell_SetRealKey(ierr,1.0D0,"EllLinConfMetric::Bnd::Robin::inf") - call Ell_SetIntKey(ierr,1,"EllLinConfMetric::Bnd::Robin::falloff") + call Ell_SetStrKey(ierr,"yes","EllLinMetric::Bnd::Robin") + call Ell_SetRealKey(ierr,1.0D0,"EllLinMetric::Bnd::Robin::inf") + call Ell_SetIntKey(ierr,1,"EllLinMetric::Bnd::Robin::falloff") + +c Elliptic solver -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 - call Ell_SetIntKey(ierr,sor_maxit,"Ell::SORmaxit"); + 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 diff --git a/src/brillq.F b/src/brillq.F index bbd9220..df6a65d 100644 --- a/src/brillq.F +++ b/src/brillq.F @@ -50,93 +50,86 @@ c q -> q [ 1 + d rho cos (n (phi)+phi0 ) / ( 1 + e rho ) ] implicit none DECLARE_CCTK_PARAMETERS + DECLARE_CCTK_FUNCTIONS logical firstcall integer qtype - integer b,c 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 - save a,b,c,r0,sr,rho0,srho,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 - sz = brill_sz - - if (axisym.eq.0) then - d = brill_d - e = brill_e - m = brill_m - n = brill_n - phi0 = brill_phi0 + if (CCTK_EQUALS(q_function,"exp")) then + qtype = 0 + else if (CCTK_EQUALS(q_function,"eppley")) then + qtype = 1 + else if (CCTK_EQUALS(q_function,"gundlach")) then + qtype = 2 + else + call CCTK_WARN(0,"Brill wave data type not recognised") end if firstcall = .false. end if + if (rho < 0) then + rho = -rho + end if + c Calculate q(rho,z) from a choice of forms. -c brill_q = 0. +c "exp" brill data if (qtype.eq.0) then - brillq = a*dabs(rho)**(2 + b) - . / dexp((rho - rho0)**2)/(rho**2 + z**2) + brillq = exp_a*dabs(rho)**(2 + exp_b) + . / dexp((rho - exp_rho0)**2)/(rho**2 + z**2) - if (sz.ne.0.0D0) then - brillq = brillq/dexp((z/sz)**2) + if (exp_sigmaz.ne.0.0D0) then + brillq = brillq/dexp((z/exp_sigmaz)**2) end if else if (qtype.eq.1) then -c brill_q = 1. This includes Eppleys choice of q. +c "Eppley" brill data. 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 = eppley_a*(rho/eppley_sigmarho)**eppley_b + & / ( 1.0D0 + ((rho**2 + z**2 - eppley_r0**2)/ + & eppley_sigmar**2)**(eppley_c/2) ) else if (qtype.eq.2) then -c brill_q = 2. This includes my (Carstens) notion of what a +c "Gundlach" brill data This includes my (Carstens) notion of what a c smooth "pure quadrupole" q should be, plus the choice of c Holz et al. - brillq = a*(rho/srho)**b - . / dexp(((rho**2 + z**2 - r0**2)/sr**2)**(c/2)) + brillq = gundlach_a*(rho/gundlach_sigmarho)**gundlach_b + & / dexp(((rho**2 + z**2 - gundlach_r0**2)/ + & gundlach_sigmar**2)**(gundlach_c/2)) else c Unknown type for q function. - call CCTK_WARN(0,"Unknown type of Brill function q") + call CCTK_WARN(0,"brillq: Unknown type of Brill function q") end if 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)) + if (CCTK_EQUALS(initial_data,"brilldata")) then + brillq = brillq*(1.0D0 + brill3d_d*rho**brill3d_m* + & cos(brill3D_n*(phi-brill3d_phi0))**2 + & / (1.0D0 + brill3d_e*rho**brill3d_m)) end if return diff --git a/src/finishbrilldata.F b/src/finishbrilldata.F index 9af9fbf..467c28d 100644 --- a/src/finishbrilldata.F +++ b/src/finishbrilldata.F @@ -55,7 +55,7 @@ c Brill metric calculated from q and Psi. y1 = y(i,j,k) z1 = z(i,j,k) - rho2 = x1**2 + y1**2 + rho2 = x1*x1 + y1*y1 rho1 = dsqrt(rho2) phi = phif(x1,y1) @@ -74,8 +74,8 @@ c The individual coefficients can be read off as 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) + gxx(i,j,k) = (e2q + (one - e2q)*y1*y1/rho2) + gyy(i,j,k) = (e2q + (one - e2q)*x1*x1/rho2) gzz(i,j,k) = e2q gxy(i,j,k) = - (one - e2q)*x1*y1/rho2 diff --git a/src/setupbrilldata2D.F b/src/setupbrilldata2D.F index f84e825..f2e5f13 100644 --- a/src/setupbrilldata2D.F +++ b/src/setupbrilldata2D.F @@ -33,11 +33,8 @@ c f DECLARE_CCTK_PARAMETERS integer i,j,k - integer nx,ny,nz integer ierr - integer, dimension(3) :: sym - CCTK_REAL x1,y1,z1,rho1 CCTK_REAL brillq,eps CCTK_REAL zp,zm,rhop,rhom @@ -50,53 +47,34 @@ 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 Epsilon for finite differencing. eps = cctk_delta_space(1) -c Initialize psi. - - brillpsi = one - -c Initialize metric. - - gxx = one - gyy = one - gzz = one - - gxy = zero - gxz = zero - gyz = zero - -c Define the symmetries for the brill GFs. - - sym(1) = 1 - sym(2) = 1 - sym(3) = 1 - - call SetCartSymVN(ierr,cctkGH,sym,'idbrilldata::brillMlinear') - call SetCartSymVN(ierr,cctkGH,sym,'idbrilldata::brillNsource') - c Set up coefficient arrays for elliptic solve. c Notice that the Cactus conventions are: c __ c \/ psi + Mlinear*psi + Nsource = 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) x1 = x(i,j,k) y1 = y(i,j,k) z1 = z(i,j,k) - rho1 = dsqrt(x1**2 + y1**2) + rho1 = dsqrt(x1*x1 + y1*y1) + + brillpsi(i,j,k) = one + + gxx(i,j,k) = one + gyy(i,j,k) = one + gzz(i,j,k) = one + + gxy(i,j,k) = zero + gxz(i,j,k) = zero + gyz(i,j,k) = zero 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 @@ -116,17 +94,15 @@ c will not be regular on the axis. . + brillq(rhom,z1,zero) . - 4.0D0*brillq(rho1,z1,zero))/eps**2 + brillNsource(i,j,k) = zero + end do end do end do -c Set coefficient Nsource = 0. - - brillNsource = zero - c Synchronization and boundaries. - call CCTK_SyncGroup(ierr,cctkGH,'idbrilldata::brillelliptic') +c call CCTK_SyncGroup(ierr,cctkGH,'idbrilldata::brillelliptic') call CartSymGN(ierr,cctkGH,'idbrilldata::brillelliptic') return diff --git a/src/setupbrilldata3D.F b/src/setupbrilldata3D.F index ecd0fd3..97cf28b 100644 --- a/src/setupbrilldata3D.F +++ b/src/setupbrilldata3D.F @@ -43,22 +43,13 @@ c c \ z rho phi phi / DECLARE_CCTK_FUNCTIONS integer i,j,k - integer nx,ny,nz integer ierr - integer, dimension(3) :: sym - 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 @@ -74,78 +65,44 @@ c Epsilon for finite differencing. eps = cctk_delta_space(1) -c Initialize psi. - - brillpsi = one - - 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) z1 = z(i,j,k) - rho2 = x1**2 + y1**2 + rho2 = x1*x1 + y1*y1 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 Define the symmetries for the brill GFs. - - sym(1) = 1 - sym(2) = 1 - sym(3) = 1 - - call SetCartSymVN(ierr,cctkGH,sym,'idbrilldata::brillMlinear') - call SetCartSymVN(ierr,cctkGH,sym,'idbrilldata::brillNsource') - -c Set up coefficient arrays for elliptic solve. -c Notice that the Cactus conventions are: -c __ -c \/ psi + Mlinear*psi + Nsource = 0 - - 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) +c Initialise psi - phi = phif(x1,y1) + brillpsi(i,j,k) = one - e2q = dexp(two*brillq(rho1,z1,phi)) +c Set up metric and coefficient arrays for elliptic solve. +c Notice that the Cactus conventions are: +c __ +c \/ psi + Mlinear*psi + Nsource = 0 c Find M using centered differences over a small c interval. +c Here we 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 + gxx(i,j,k) = e2q + (one - e2q)*y1*y1/rho2 + gyy(i,j,k) = e2q + (one - e2q)*x1*x1/rho2 + gzz(i,j,k) = e2q + gxy(i,j,k) = - (one - e2q)*x1*y1/rho2 + brillMlinear(i,j,k) = 0.25D0/e2q . *(brillq(rho1,z1+eps,phi) . + brillq(rho1,z1-eps,phi) @@ -154,8 +111,20 @@ c interval. . - 4.0D0*brillq(rho1,z1,phi)) . / eps**2 + brillMlinear(i,j,k) = brillMlinear(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 + else + gxx(i,j,k) = one + gyy(i,j,k) = one + gzz(i,j,k) = one + gxy(i,j,k) = zero + brillMlinear(i,j,k) = 0.25D0/e2q . *(brillq(rho1,z1+eps,phi) . + brillq(rho1,z1-eps,phi) @@ -165,32 +134,21 @@ c interval. 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. + gxz(i,j,k) = zero + gyz(i,j,k) = zero - if (rho1.gt.rhofudge) then - brillMlinear(i,j,k) = brillMlinear(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 +c Set coefficient Nsource = 0 + brillNsource(i,j,k) = zero end do end do end do -c Set coefficient Nsource = 0. - - brillNsource = zero c Synchronization and boundaries. - call CCTK_SyncGroup(ierr,cctkGH,'idbrilldata::brillelliptic') - call CartSymGN(ierr,cctkGH,'idbrilldata::brillelliptic') +c call CCTK_SyncGroup(ierr,cctkGH,'idbrilldata::brillelliptic') +c call CartSymGN(ierr,cctkGH,'idbrilldata::brillelliptic') return end diff --git a/test/test_brilldata_1.par b/test/test_brilldata_1.par index 6dc6543..5d673ba 100644 --- a/test/test_brilldata_1.par +++ b/test/test_brilldata_1.par @@ -17,13 +17,13 @@ cactus::cctk_itlast = 0 # Brill wave initial data -ADMBase::initial_data = "brilldata" +ADMBase::initial_data = "brilldata2D" ADMBase::metric_type = "physical" # Brill wave parameters -idbrilldata::brill_q = 2 -idbrilldata::brill_a = 5.0 +idbrilldata::q_function = "gundlach" +idbrilldata::gundlach_a = 5.0 # Elliptic solver. diff --git a/test/test_brilldata_2.par b/test/test_brilldata_2.par index 858a256..630cd64 100644 --- a/test/test_brilldata_2.par +++ b/test/test_brilldata_2.par @@ -17,14 +17,14 @@ cactus::cctk_itlast = 0 # Brill wave initial data -ADMBase::initial_data = "brilldata" +ADMBase::initial_data = "brilldata2D" ADMBase::metric_type = "static conformal" StaticConformal::conformal_storage = "factor" # Brill wave parameters -idbrilldata::brill_q = 2 -idbrilldata::brill_a = 5.0 +idbrilldata::q_function = "gundlach" +idbrilldata::gundlach_a = 5.0 # Elliptic solver. diff --git a/test/test_wsor.par b/test/test_wsor.par index f3fa563..a9b0337 100644 --- a/test/test_wsor.par +++ b/test/test_wsor.par @@ -1,4 +1,4 @@ -!DESC "Brill wave initial data (Holz type), amplitude 5, SOR, robin boundaries, physical metric" +!DESC "Brill wave initial data (Physical, SOR)" ActiveThorns = "admbase staticconformal ellbase ellsor boundary pugh pughslab pughreduce cartgrid3d idbrilldata ioascii iobasic ioutil time" @@ -16,13 +16,13 @@ cactus::cctk_itlast = 0 # Brill wave initial data -ADMBase::initial_data = "brilldata" +ADMBase::initial_data = "brilldata2D" ADMBase::metric_type = "physical" # Brill wave parameters -idbrilldata::brill_q = 2 -idbrilldata::brill_a = 5.0 +idbrilldata::q_function = "gundlach" +idbrilldata::gundlach_a = 5.0 # Elliptic solver. @@ -50,4 +50,3 @@ IOASCII::out1D_vars = "admbase::gxx admbase::gyy admbase::gzz admbase::gxy idbri IO::new_filename_scheme = "no" -pugh::enable_all_storage = "yes" # fixes a seg fault for now |