A Multi-Patch Wave Toy

Erik Schnetter <schnetter@aei.mpg.de>
Luis Lehner <lehner@lsu.edu>
Manuel Tiglio <tiglio@lsu.edu>

April 20, 2017

1 Physical System

The massless wave equation for a scalar field ϕ can be written as

μ(γμνd ν) = 0

where γμν = ggμν, and dμ μϕ. Using the expressions g = αh and

gμν = 1α2 βiα2 βjα2 γij βiβjα2 ,

where γij is the inverse of the three metric, the equation can be rewritten as

ϕ̇ = Π, (1) Π̇ = βi iΠ + α hi h α βiΠ + h α Hijd j + α hdit hβi α Πα ht h α , (2) di ̇ = iΠ (3)

with Hij α2γij βiβj. The non-shift speed modes with respect to a boundary with normal ni are

v± = λΠ + Hijn idj

and the shift speed modes are dA, with A transversal directions.

The physical energy is

E = 1 2 1 α Π2 + Hijd idj hdx3

and the way the equations above have been written this energy is not increasing in the stationary background case if homogeneous boundary conditions are given at outer boundaries, also at the discrete level (replacing above i by Di).

The implementation of the field equations in the code differs slightly from the above. Assuming a stationary background and expanding out the derivative operator, we obtain the equations in the final form

ϕ̇ = Π, (4) Π̇ = βi iΠ + α hi h α βiΠ + α hi h α Hij jϕ + α h h α Hij ijϕ. (5)

2 Non-linear addition to the multi-patch wave toy

Simple addition to the wave multi-patch toy to get started on this.

2.1 Physical System

This is just a non-linear wave equation obtained, a very minor modification to the wave equation thorn adding a different initial data and a slight modification to the right hand side. A reference to this is in a paper by Liebling to appear in Phys. Rev. D (2005). The non-linear wave equation is written as

u(γμνd ν) = ϕp

where γμν = ggμν, and dμ = μϕ and p must be an odd integer 3.

The initial data coded is given by

ϕ = Ae(r1R)2δ2 (6) Π = μϕ,r + Ω(yϕ,x xϕ,y) (7) di = ϕ,i (8)

with r̃2 = 𝜖xx2 + 𝜖yy2 + z2.

The parameters used for this initial data are given some distinct names to avoid conflicts with existing ones and are as follows:

Note, as the solution is not known, one must set for now the incoming fields to 0. CPBC might one day be put... though who knows :-)

3 Formulations

U = ρvxvyvz T = ρvi T (9) tU = Ai iU + (10) ||ni|| 1 (11)

3.1 dt

Setting ρ = tu.

RHS:

Hij = α2γij βiβj (12) tu = ρ (13) tρ = βi iρ + α 𝜖 i 𝜖 α βiρ + Hijv j (14) tvi = iρ (15)

Propagation matrix:

Ax = 2βx βxβx + α2γxx βxβy + γxyα2 βxβz + α2γxz 1 0 0 0 0 0 0 0 0 0 0 0 (16)

An = 2βin i βiβj + α2γij n i ni 0 (17)

Eigensystem:

λ1 = 0 , w1 = 0001 T (18) λ2 = 0 , w2 = 0010 T (19) λ3 = βx αγxx , w 3 = βx αγxxβxβx + α2γxxβxβy + α2γxyβxβz + α2γxz T (20) λ4 = βx + αγxx , w 4 = βx + αγxxβxβx + α2γxxβxβy + α2γxyβxβz + α2γxz T (21)

λt = 0 , wt = 0ti T (22) λ± = βin i ± αγij ni nj , w± = βin i ± αγij ni njHijnj T (23)

3.2 d0

Setting ρ = nu with na = Dat, leading to ρ = 0u = (1α)tu (1α)βiiu.

RHS:

Hij = γij (24) tu = αρ + βiv i (25) tρ = βi iρ + 1 𝜖i𝜖αHijv j + ρ 𝜖i𝜖βi (26) tvi = βj jvi + vjiβj + iαρ (27)

Propagation matrix:

Ax = βxαγxxαγxyαγxz α βx 0 0 0 0 βx 0 0 0 0 βx (28)

An = βini αγijni αni βknkδij (29)

Eigensystem:

λt = βin i , wt = 0 γijnjnktk + γjknjnkti T (30) λ± = βin i ± αγij ni nj , w± = ±γij ni njγijnj T (31)

3.3 dk

(This section very probably contains errors, say Erik on 2005-04-13.)

Setting ρ = ku with a “Killing” vector ka = txa, leading to ρ = (1α)tu.

RHS:

Hij = γij βiβj α2 (32) tu = αρ (33) tρ = βi α iαρ + 1 𝜖i𝜖 βiρ + αHijv j (34) tvi = iαρ (35)

Propagation matrix:

Ax = 2βxα βxβx α2 + γxxα βxβy α2 + γxyα βxβz α2 + γxz α 0 0 0 0 0 0 0 0 0 0 0 (36)

Eigensystem:

λ1 = 0 , w1 = 0 βxβz + α2γxz0βxβx α2γxx T (37) λ2 = 0 , w2 = 0 βxβy + α2γxyβxβx α2γxx0 T (38) λ3 = βx αγxx , w 3 = βx γxxα00 T (39) λ4 = βx + αγxx , w 4 = βx + γxxα00 T (40)

References

[1]   Gioel Calabrese, Luis Lehner, Dave Neilsen, Jorge Pullin, Oscar Reula, Olivier Sarbach, Manuel Tiglio, Novel finite-differencing techniques for numerical relativity: application to black-hole excision, Class. Quantum Grav. 20, L245 (2003), gr-qc/0302072.

4 Parameters




amplitude
Scope: private  REAL



Description: Amplitude



Range   Default: 1.0
*:*






anl
Scope: private  REAL



Description: Amplitude of the non-linear Gaussian



Range   Default: 0.0
*:*






bound
Scope: private  STRING



Description: Boundary condition



Range   Default: static
.*
any registered boundary condition






compute_second_derivative_from_first_derivative
Scope: private  BOOLEAN



Description: Take first derivative twice to compute second derivate



  Default: no






deltanl
Scope: private  REAL



Description: sigma of the non-linear Gaussian



Range   Default: 0.0
*:*






eps
Scope: private  REAL



Description: A small number



Range   Default: 1.0e-10
0:*






epsx
Scope: private  REAL



Description: eps in x-direction of the non-linear Gaussian



Range   Default: 0.0
*:*






epsy
Scope: private  REAL



Description: eps in y-direction of the non-linear Gaussian



Range   Default: 0.0
*:*






initial_data
Scope: private  KEYWORD



Description: Type of initial data



Range   Default: plane
linear
x and y coordinates
plane
Plane wave
Gaussian
Gaussian wave packet
GaussianNonLinear
Gaussian wave packet for the non-linear RHS
GeneralMultipole
Multipole with arbitrary l and m
multipole
L=1 initial data, Gaussian in r
multipole l=1, m=0
L=1 m=0 initial data, Gaussian in r, u=0
multipole l=1, m=1
L=1 m=1 initial data, Gaussian in r, u=0
multipole l=1, m=-1
L=1 m=-1 initial data, Gaussian in r, u=0
multipole l=2
L=2 initial data, Gaussian in r
multipole l=2, u=0
L=2 initial data, Gaussian in r, u=0
multipole l=2, m=1
L=2 m=1 initial data, Gaussian in r, u=0
multipole l=2, m=-1
L=2 m=-1 initial data, Gaussian in r, u=0
multipole l=2, m=2
L=2 m=2 initial data, Gaussian in r, u=0
multipole l=2, m=-2
L=2 m=-2 initial data, Gaussian in r, u=0
multipole l=2, m=-2
L=2 m=-2 initial data, Gaussian in r, u=0
multipole l=4, m=0
L=4 m=0 initial data, Gaussian in r, u=0
multipole l=4, m=1
L=4 m=1 initial data, Gaussian in r, u=0
multipole l=4, m=-1
L=4 m=-1 initial data, Gaussian in r, u=0
multipole l=4, m=2
L=4 m=-2 initial data, Gaussian in r, u=0
multipole l=4, m=-2
L=4 m=-2 initial data, Gaussian in r, u=0
multipole l=4, m=3
L=4 m=3 initial data, Gaussian in r, u=0
multipole l=4, m=-3
L=4 m=-3 initial data, Gaussian in r, u=0
multipole l=4, m=4
L=4 m=4 initial data, Gaussian in r, u=0
multipole l=4, m=-4
L=4 m=-4 initial data, Gaussian in r, u=0
noise
Random noise
debug
number of current patch and grid point index






initial_data_analytic_derivatives
Scope: private  BOOLEAN



Description: Calculate spatial derivatives of the initial data analytically?



  Default: no






lapse
Scope: private  REAL



Description: Lapse function multiplier



Range   Default: 1.0
*:*
must not be zero






mass
Scope: private  REAL



Description: Mass M



Range   Default: 1.0
*:*






metric
Scope: private  KEYWORD



Description: Global metric



Range   Default: Minkowski
Minkowski
Minkowski
Kerr-Schild
Kerr-Schild
Kerr
Kerr Black Hole






multipole_l
Scope: private  INT



Description: For GeneralMultipole initial data: degree of spherical harmonic function



Range   Default: 2
0:*
A positive integer






multipole_m
Scope: private  INT



Description: For GeneralMultipole initial data: order of spherical harmonic function



Range   Default: 2
*:*
An integer -l<=m<=l






multipole_s
Scope: private  INT



Description: For GeneralMultipole initial data: spin weight spherical harmonic function



Range   Default: (none)
*:*
A positive integer






munl
Scope: private  REAL



Description: Speed of the non-linear Gaussian



Range   Default: 0.0
*:*






nonlinearrhs
Scope: private  BOOLEAN



Description: Add a non-linear term to the RHS?



  Default: no






numevolvedvars
Scope: private  INT



Description: The number of evolved variables in this thorn



Range   Default: 5
5:5
five






omenl
Scope: private  REAL



Description: Omega of the non-linear Gaussian



Range   Default: 0.0
*:*






outer_bound
Scope: private  STRING



Description: outer boundary



Range   Default: solution
zero
set all characteristics to zero
solution
use the same analytic solution as for the initial data
dirichlet
set all fields to zero
radiative
radiative boundary
none
no boundary condition






outer_penalty_bound
Scope: private  STRING



Description: outer penalty boundary



Range   Default: zero
zero
set all characteristics to zero
solution
use the same analytic solution as for the initial data






powerrhs
Scope: private  REAL



Description: Exponent in the non-linear RHS term



Range   Default: 5.0
3.0:13.0






radius
Scope: private  REAL



Description: Radius of the Gaussian



Range   Default: 0.0
0:*






recalculate_rhs
Scope: private  BOOLEAN



Description: Recalculate the RHSs in the ANALYSIS timebin



  Default: yes






rhsbound
Scope: private  STRING



Description: Boundary condition during RHS evaluation



Range   Default: none
.*
any registered boundary condition






rnl
Scope: private  REAL



Description: How fat the non-linear Gaussian is (r-R)



Range   Default: 0.0
*:*






shift
Scope: private  REAL



Description: Shift vector addition



Range   Default: 0.0
*:*






shift_interpolation_type
Scope: private  KEYWORD



Description: Setting for the interpolating vector field bî (only used for the db formulation)



Range   Default: shift
shift
Set bî = betaî (corresponds to d0 formulation)
zero
Set bî = 0 (corresponds to dk formulation)






shift_omega
Scope: private  REAL



Description: Rotational shift vector addition about z axis



Range   Default: 0.0
*:*






space_offset
Scope: private  REAL



Description: Space offset



Range   Default: 0.0
*:*






spin
Scope: private  REAL



Description: Spin a=J/M2̂



Range   Default: 0.0
-1:+1






time_offset
Scope: private  REAL



Description: Time offset



Range   Default: 0.0
*:*






wave_number
Scope: private  REAL



Description: Wave number



Range   Default: 0.0
*:*






width
Scope: private  REAL



Description: Width of the Gaussian



Range   Default: 1.0
(0:*






mol_num_evolved_vars
Scope: shared from METHODOFLINES INT



5 Interfaces

General

Implements:

llamawavetoy

Inherits:

grid

coordinates

globalderivative

summationbyparts

interpolate

Grid Variables

5.0.1 PRIVATE GROUPS




  Group Names    Variable Names    Details   




scalar   compact0
u   descriptionThe scalar of the scalar wave equation fame
  dimensions3
  distributionDEFAULT
  group typeGF
  tagstensortypealias=”scalar”
  timelevels2
 variable typeREAL




density   compact0
rho   descriptionTime derivative of u
  dimensions3
  distributionDEFAULT
  group typeGF
  tagstensortypealias=”scalar”
  timelevels2
 variable typeREAL




dx_scalar   compact0
dx_u   descriptionSpatial derivatives of u
dy_u   dimensions3
dz_u   distributionDEFAULT
  group typeGF
  tagstensortypealias=”scalar”
  timelevels2
 variable typeREAL




dxx_scalar   compact0
dxx_u   descriptionSpatial derivatives of u
dyy_u   dimensions3
dzz_u   distributionDEFAULT
dxy_u   group typeGF
dxz_u   tagstensortypealias=”scalar”
dyz_u   timelevels2
 variable typeREAL




dx_density   compact0
dx_rho   descriptionSpatial derivatives of rho
dy_rho   dimensions3
dz_rho   distributionDEFAULT
  group typeGF
  tagstensortypealias=”scalar”
  timelevels2
 variable typeREAL




velocity   compact0
vx   descriptionSpatial derivative of u
vy   dimensions3
vz   distributionDEFAULT
  group typeGF
  tagstensortypealias=”scalar”
  timelevels2
 variable typeREAL








  Group Names    Variable Names    Details   




debug   compact0
vxdebug   descriptionSpatial derivative of u
vydebug   dimensions3
vzdebug   distributionDEFAULT
  group typeGF
  tagstensortypealias=”scalar”
  timelevels1
 variable typeREAL




scalardot   compact0
udot   descriptionRHS of u
  dimensions3
  distributionDEFAULT
  group typeGF
  tagstensortypealias=”scalar” Prolongation=”None”
  timelevels1
 variable typeREAL




densitydot   compact0
rhodot   descriptionRHS of rho
  dimensions3
  distributionDEFAULT
  group typeGF
  tagstensortypealias=”scalar” Prolongation=”None”
  timelevels1
 variable typeREAL




velocitydot   compact0
vxdot   descriptionRHS ov v
vydot   dimensions3
vzdot   distributionDEFAULT
  group typeGF
  tagstensortypealias=”scalar” Prolongation=”None”
  timelevels1
 variable typeREAL




constraints   compact0
wx   descriptionIntegrability condition
wy   dimensions3
wz   distributionDEFAULT
  group typeGF
  tagstensortypealias=”scalar” tensorparity=-1 Prolongation=”None”
  timelevels1
 variable typeREAL




difference_v   compact0
diff_vx   descriptionDifference between v_i and d/dxî u
diff_vy   dimensions3
diff_vz   distributionDEFAULT
  group typeGF
  tagstensortypealias=”scalar” Prolongation=”None”
  timelevels1
 variable typeREAL








  Group Names    Variable Names    Details   




velocity_squared   compact0
v2   descriptionVelocity squared
  dimensions3
  distributionDEFAULT
  group typeGF
  tagstensortypealias=”scalar” Prolongation=”None”
  timelevels1
 variable typeREAL




scalarenergy   compact0
energy   descriptionEnergy of the scalar field
  dimensions3
  distributionDEFAULT
  group typeGF
  tagstensortypealias=”scalar” Prolongation=”None”
  timelevels1
 variable typeREAL




errors   compact0
error   descriptionError of the solution
error_rho   dimensions3
error_vx   distributionDEFAULT
error_vy   group typeGF
error_vz   tagstensortypealias=”scalar” Prolongation=”None”
exact   timelevels1
exact_rho  variable typeREAL




errorsperiodic   compact0
errorperiodic   descriptionError for a solution which is known to be periodic
errorperiodic_rho   dimensions3
  distributionDEFAULT
  group typeGF
  tagstensortypealias=”scalar” Prolongation=”None”
  timelevels1
 variable typeREAL




metric   compact0
gxx   descriptionSpatial background metric
gxy   dimensions3
gxz   distributionDEFAULT
gyy   group typeGF
gyz   tagstensortypealias=”scalar” Prolongation=”None”
gzz   timelevels1
 variable typeREAL




inverse_metric   compact0
guxx   descriptionInverse of the spatial background metric
guxy   dimensions3
guxz   distributionDEFAULT
guyy   group typeGF
guyz   tagstensortypealias=”scalar” Prolongation=”None”
guzz   timelevels1
 variable typeREAL








  Group Names    Variable Names    Details   




lapse   compact0
alpha   descriptionSpatial background metric
  dimensions3
  distributionDEFAULT
  group typeGF
  tagstensortypealias=”scalar” Prolongation=”None”
  timelevels1
 variable typeREAL




shift   compact0
betax   descriptionSpatial background metric
betay   dimensions3
betaz   distributionDEFAULT
  group typeGF
  tagstensortypealias=”scalar” Prolongation=”None”
  timelevels1
 variable typeREAL




volume_element   compact0
epsilon   descriptionVolume element due to the spatial background metric
  dimensions3
  distributionDEFAULT
  group typeGF
  tagstensortypealias=”scalar” Prolongation=”None”
  timelevels1
 variable typeREAL




min_spacing   compact0
min_spacing   descriptionMinimum grid spacing
  dimensions3
  distributionDEFAULT
  group typeGF
  tagstensortypealias=”scalar” Prolongation=”None”
  timelevels1
 variable typeREAL




6 Schedule

This section lists all the variables which are assigned storage by thorn Llama/LlamaWaveToy. Storage can either last for the duration of the run (Always means that if this thorn is activated storage will be assigned, Conditional means that if this thorn is activated storage will be assigned for the duration of the run if some condition is met), or can be turned on for the duration of a schedule function.

Storage

 

Always:  
scalar[2] density[2] velocity[2]  
dx_scalar[2]  
dxx_scalar[2]  
dx_density[2]  
errors  
metric inverse_metric lapse shift volume_element 
scalardot densitydot velocitydot  
   

Scheduled Functions

CCTK_STARTUP

  lwt_startup

  register banner with cactus

 

 Language:c
 Options: meta
 Type: function

MoL_Register (conditional)

  lwt_register_mol

  register variables with mol

 

 Language:c
 Options: meta
 Type: function

CCTK_ANALYSIS (conditional)

  lwt_calc_rhs

  calculate the rhs

 

 Language:fortran
 Storage: scalardot
   densitydot
   velocitydot
 Sync: scalardot
   densitydot
   velocitydot
 Triggers: scalardot
   densitydot
   velocitydot
 Type: function

CCTK_ANALYSIS

  lwt_calcenergy

  calculate the energy of the scalar field

 

 Language:fortran
 Storage: scalarenergy
 Sync: scalarenergy
 Triggers: scalarenergy
 Type: function

CCTK_ANALYSIS

  lwt_error

  calculate errors of the solution

 

 Language:fortran
 Storage: errorsperiodic
 Triggers: errors
   errorsperiodic
 Type: function

CCTK_ANALYSIS

  lwt_min_spacing

  calculate the smallest grid spacing

 

 Language:fortran
 Storage: min_spacing
 Sync: min_spacing
 Triggers: min_spacing
 Type: function

CCTK_INITIAL

  lwt_init_metric

  initialise the metric

 

 Language:fortran
 Type: function

CCTK_INITIAL

  lwt_calc_inverse_metric

  transform the metric

 

 After: lwt_init_metric
  Language:fortran
 Type: function

CCTK_INITIAL

  lwt_init

  initialise the system

 

 After: lwt_init_metric
  Language:fortran
 Type: function

MoL_CalcRHS

  lwt_calc_rhs

  calculate the rhs

 

 Language:fortran
 Type: function

MoL_PostStep

  lwt_outerboundary

  apply outer boundaries

 

 Language:fortran
 Type: function

MoL_RHSBoundaries

  lwt_rhs_outerboundary

  apply mol rhs outer boundaries (eg. radiative boundary condition)

 

 Language:fortran
 Type: function

MoL_PostStep

  lwt_boundaries

  select the boundary condition

 

 After: lwt_outerbound
  Language:fortran
 Options: level
 Sync: scalar
   density
   velocity
 Type: function

MoL_PostStep

  applybcs

  apply boundary conditions

 

 After:lwt_boundaries
 Type:group

Aliased Functions

 

Alias Name:        Function Name:
ApplyBCs LWT_ApplyBCs