CanudaX_BSSNMoL: An Einstein Toolkit thorn for evolving Einstein’s equations in the BSSN formalism

Cheng-Hsin Cheng

April 15, 2026

Abstract

CanudaX_BSSNMoL is the CarpetX-compatible version of LeanBSSNMoL, which solves Einstein’s equations in the BSSN formalism. CanudaX_BSSNMoL uses finite-difference stencils which are accurate up to 8th order, compared with only 6th order in LeanBSSNMoL. It relies on the ODESolvers thorn for the time integration.

1 CanudaX_BSSNMoL

The Lean code was first introduced in [1], to evolve vacuum spacetimes with the BSSN formalism. It has since been modified into LeanBSSNMoL to run with the MoL thorn for the time integration and generalized to include matter terms through the TmunuBase thorn. LeanBSSNMoL was made publicly available through the Canuda numerical relativity library [2], and has since been distributed also as a part of the Einstein Toolkit.

CanudaX_BSSNMoL is the CarpetX-compatible version of LeanBSSNMoL which further introduces several new physics and computational capabilities:

The bulk of the code is written in C++ and should be simple to follow – emphasis has been given to readability. The main part of the code, where the right-hand side of the evolution equations is computed, can be found in the file calc_bssn_rhs.cxx. For the conformal factor the code uses the “W” version, i.e. \(W=\gamma ^{-1/6}\), and it employs the usual “1+log” and “Gamma-driver” gauge conditions by default.

2 Obtaining this thorn

CanudaX_BSSNMoL is included with the Einstein Toolkit and can also be obtained through the Canuda numerical relativity library [2].

References

[1]   U. Sperhake, “Binary black-hole evolutions of excision and puncture data,” Phys. Rev. D 76 (2007), 104015 doi:10.1103/PhysRevD.76.104015 [arXiv:gr-qc/0606079 [gr-qc]].

[2]   H. Witek, M. Zilhao, G. Bozzola, C.-H. Cheng, A. Dima, M. Elley, G. Ficarra, T. Ikeda, R. Luna, C. Richards, N. Sanchis-Gual, H. Okada da Silva. “Canuda: a public numerical relativity library to probe fundamental physics,” Zenodo (2023) doi: 10.5281/zenodo.3565474

[3]   T. W. Baumgarte and D. Hilditch, “Shock-avoiding slicing conditions: Tests and calibrations,” Phys. Rev. D 106, no.4, 044014 (2022) doi:10.1103/PhysRevD.106.044014 [arXiv:2207.06376 [gr-qc]].

[4]   Z. B. Etienne, “Improved moving-puncture techniques for compact binary simulations,” Phys. Rev. D 110, no.6, 064045 (2024) doi:10.1103/PhysRevD.110.064045 [arXiv:2404.01137 [gr-qc]].

[5]   M. Radia, U. Sperhake, A. Drew, K. Clough, P. Figueras, E. A. Lim, J. L. Ripley, J. C. Aurrekoetxea, T. França and T. Helfer, “Lessons for adaptive mesh refinement in numerical relativity,” Class. Quant. Grav. 39, no.13, 135006 (2022) doi:10.1088/1361-6382/ac6fa9 [arXiv:2112.10567 [gr-qc]].

3 Parameters

add_ko_dissipation
 Scope: private BOOLEAN
Description: Add Kreiss-Oliger dissipation?
Default: yes

beta_alp
 Scope: private REAL
Description: Exponent for lapse in front of Gammaî in shift condition
Range Default: 0.0
*:*
Anything possible, default is zero

beta_gamma
 Scope: private REAL
Description: Prefactor to the Gammaî term in shift condition
Range Default: 0.75
*:*
Anything possible

boundary_conditions
 Scope: private KEYWORD
Description: Outer boundary conditions for BSSN variables
Range Default: radiative
radiative
Use radiative boundary conditions from NewRadX
none
Do not impose boundary conditions here, but let CarpetX set it

calculate_constraints
 Scope: private BOOLEAN
Description: Calculate the BSSN constraints?
Default: no

calculate_constraints_every
 Scope: private INT
Description: Calculate the BSSN constraints every N iterations
Range Default: 1
*:*
0 or a negative value means never compute them

chi_gamma
 Scope: private REAL
Description: Coefficient to Yo-term in the ∖tilde{Gamma}equation. See gr-qc/0209066
Range Default: 0.666666666666667
*:*
2/3 is a good value; the sign must be the same as betak,k

conf_fac_floor
 Scope: private REAL
Description: Minimal value of conformal factor
Range Default: 1.0d-04
*:*
Any value possible

derivs_order
 Scope: private INT
Description: Accuracy of spatial finite difference stencils in RHS calculation
Range Default: 4
4
4th order accurate stencils
6
6th order accurate stencils
8
8th order accurate stencils

diss_eps
 Scope: private REAL
Description: Strength coefficient for Kreiss-Oliger dissipation
Range Default: 0.15
0.0:*
0 to turn off dissipation, most commonly use 0.15

diss_order
 Scope: private INT
Description: Accuracy of Kreiss-Oliger dissipation stencils; should be 1 highehr than derivs_order
Range Default: 5
5
5th order accurate stencils
7
7th order accurate stencils
9
9th order accurate stencils

eta_beta
 Scope: private REAL
Description: Damping parameter in live shift
Range Default: 1
0:*
non-negative

eta_transition
 Scope: private BOOLEAN
Description: Use an r-dependent damping parameter eta?
Default: no

eta_transition_eps
 Scope: private REAL
Description: Minimal value of radius squared for eta_transition
Range Default: 1.0d-06
0:*
Any value possible

eta_transition_r
 Scope: private REAL
Description: Transition radius for r-dependent eta
Range Default: 1
0:*
non-negative

impose_conf_fac_floor_at_initial
 Scope: private BOOLEAN
Description: Use floor value on initial data?
Default: no

kappa_alpha
 Scope: private REAL
Description: Coefficient in slicing condition - default is 2 for 1+log
Range Default: 2.0
*:*
Any value possible

make_aa_tracefree
 Scope: private BOOLEAN
Description: Remove trace of ∖tilde{A}_{ij}after each timestep?
Default: yes

n_aij
 Scope: private INT
Description: n power of outgoing boundary rˆn fall off rate for conformal tracefree extrinsic curvature
Range Default: 2
0:2
2 is reasonable

n_alpha
 Scope: private INT
Description: n power of outgoing boundary rˆn fall off rate for lapse
Range Default: 1
0:2
1 is my guess

n_beta
 Scope: private INT
Description: n power of outgoing boundary rˆn fall off rate for shift
Range Default: 1
0:2
1 is my guess

n_conf_fac
 Scope: private INT
Description: n power of outgoing boundary rˆn fall off rate for conformal factor
Range Default: 1
0:2
1 is reasonable

n_gammat
 Scope: private INT
Description: n power of outgoing boundary rˆn fall off rate for conformal connection
Range Default: 1
0:2
Maybe 1?

n_hij
 Scope: private INT
Description: n power of outgoing boundary rˆn fall off rate for conformal metric
Range Default: 1
0:2
1 is reasonable

n_trk
 Scope: private INT
Description: n power of outgoing boundary rˆn fall off rate for trace of extrinsic curvature
Range Default: 2
0:2
2 is reasonable

precollapsed_lapse
 Scope: private BOOLEAN
Description: Initialize lapse as alp*psiˆ{ -2}?
Default: no

refinement_criterion
 Scope: private KEYWORD
Description: How to compute regrid error for AMR?
Range Default: fixed
fixed
Distance to origin
see [1] below
Conformal factor minus 1
conf_fac_laplacian
Laplacian of conformal factor
see [1] below
Sum absolute values of each term in Hamiltonian constraint

[1] 

conf\_fac\_minus\_one

[1] 

ham\_constraint\_absolute

rescale_shift_initial
 Scope: private BOOLEAN
Description: Initialize shift as psiˆ { -2}beta ?
Default: no

reset_dethh
 Scope: private BOOLEAN
Description: Reset determinant of conformal metric?
Default: yes

slicing_condition
 Scope: private KEYWORD
Description: Lapse slicing condition
Range Default: 1+log
geodesic
geodesic slicing
harmonic
harmonic slicing
1+log
1+log slicing
shock-avoiding
shock-avoiding slicing from arXiv:2207.06376 [gr-qc]
SSL
slow-start lapse from arXiv:2404.01137 [gr-qc]

ssl_amp
 Scope: private REAL
Description: Amplitude of slow-start lapse Gaussian
Range Default: 0.6
0:*
Anything positive; 0.6 was used in arXiv:2404.01137 [gr-qc]

ssl_sigma
 Scope: private REAL
Description: Width of slow-start lapse Gaussian
Range Default: 20.0
0:*
Anything positive; 20 was used in arXiv:2404.01137 [gr-qc]

use_advection_stencils
 Scope: private BOOLEAN
Description: Use lopsided stencils for advection derivatives?
Default: yes

use_local_error_estimate
 Scope: private BOOLEAN
Description: Compute local regrid error to refine/derefine cells?
Default: no

zeta_alpha
 Scope: private REAL
Description: Prefactor to the advection term in slicing condition
Range Default: 1
*:*
Anything possible

zeta_beta
 Scope: private REAL
Description: Prefactor to the advection term in shift condition
Range Default: 1
0:*
non-negative

stress_energy_state
 Scope: restricted BOOLEAN
Description: Add matter contributions from TmunuBaseX to RHS and constraints?
Default: no

z_is_radial
 Scope: restricted BOOLEAN
Description: Use with multipatch?
Default: no

4 Interfaces

General

Implements:

canudax_bssnmol

Inherits:

admbasex

carpetxregrid

tmunubasex

Grid Variables

4.0.1 PUBLIC GROUPS

  Group Names     Variable Names   Details    
alpha centering centering={0 0 0}
alpha compact 0
description Lapse function
dimensions 3
distribution DEFAULT
group type GF
tags parities={+1 +1 +1} rhs=”alpha_rhs”
timelevels 1
variable type REAL
alpha_rhs centering centering={0 0 0}
alpha_rhs compact 0
description Time derivative of lapse function
dimensions 3
distribution DEFAULT
group type GF
tags parities={+1 +1 +1} checkpoint=”no”
timelevels 1
variable type REAL
beta centering centering={0 0 0}
betx compact 0
bety description Shift vector betaî
betz dimensions 3
distribution DEFAULT
group type GF
tags parities={-1 +1 +1 +1 -1 +1 +1 +1 -1} rhs=”beta_rhs”
timelevels 1
variable type REAL
beta_rhs centering centering={0 0 0}
betx_rhs compact 0
bety_rhs description Time derivative of shift vector
betz_rhs dimensions 3
distribution DEFAULT
group type GF
tags parities={-1 +1 +1 +1 -1 +1 +1 +1 -1} checkpoint=”no”
timelevels 1
variable type REAL
conf_fac centering centering={0 0 0}
conf_fac compact 0
description Conformal factor W=det(gamma_{ij}){ˆ -1/6}
dimensions 3
distribution DEFAULT
group type GF
tags parities={+1 +1 +1} rhs=”conf_fac_rhs”
timelevels 1
variable type REAL
conf_fac_rhs centering centering={0 0 0}
conf_fac_rhs compact 0
description Time derivative of conformal factor
dimensions 3
distribution DEFAULT
group type GF
tags parities={+1 +1 +1} checkpoint=”no”
timelevels 1
variable type REAL

  Group Names     Variable Names   Details    
hmetric centering centering={0 0 0}
hxx compact 0
hxy description Conformal metric ∖tilde{gamma}_{ij}
hxz dimensions 3
hyy distribution DEFAULT
hyz group type GF
hzz tags parities={+1 +1 +1 -1 -1 +1 -1 +1 -1 +1 +1 +1 +1 -1 -1 +1 +1 +1} rhs=”hmetric_rhs”
timelevels 1
variable type REAL
hmetric_rhs centering centering={0 0 0}
hxx_rhs compact 0
hxy_rhs description Time derivative of conformal metric
hxz_rhs dimensions 3
hyy_rhs distribution DEFAULT
hyz_rhs group type GF
hzz_rhs tags parities={+1 +1 +1 -1 -1 +1 -1 +1 -1 +1 +1 +1 +1 -1 -1 +1 +1 +1} checkpoint=”no”
timelevels 1
variable type REAL
hcurv centering centering={0 0 0}
axx compact 0
axy description Conformal tracefree part of extrinsic curvature ∖tilde{A}_{ij}
axz dimensions 3
ayy distribution DEFAULT
ayz group type GF
azz tags parities={+1 +1 +1 -1 -1 +1 -1 +1 -1 +1 +1 +1 +1 -1 -1 +1 +1 +1} rhs=”hcurv_rhs”
timelevels 1
variable type REAL
hcurv_rhs centering centering={0 0 0}
axx_rhs compact 0
axy_rhs description Time derivative of conformal tracefree extrinsic curvature
axz_rhs dimensions 3
ayy_rhs distribution DEFAULT
ayz_rhs group type GF
azz_rhs tags parities={+1 +1 +1 -1 -1 +1 -1 +1 -1 +1 +1 +1 +1 -1 -1 +1 +1 +1} checkpoint=”no”
timelevels 1
variable type REAL
trk centering centering={0 0 0}
tracek compact 0
description Trace of extrinsic curvature
dimensions 3
distribution DEFAULT
group type GF
tags parities={+1 +1 +1} rhs=”trk_rhs”
timelevels 1
variable type REAL
trk_rhs centering centering={0 0 0}
tracek_rhs compact 0
description Time derivative of Trace of extrinsic curvature
dimensions 3
distribution DEFAULT
group type GF
tags parities={+1 +1 +1} checkpoint=”no”
timelevels 1
variable type REAL

  Group Names     Variable Names   Details    
gammat centering centering={0 0 0}
gammatx compact 0
gammaty description Conformal connection function ∖tilde{Gamma}î
gammatz dimensions 3
distribution DEFAULT
group type GF
tags parities={-1 +1 +1 +1 -1 +1 +1 +1 -1} rhs=”gammat_rhs”
timelevels 1
variable type REAL
gammat_rhs centering centering={0 0 0}
gammatx_rhs compact 0
gammaty_rhs description Time derivative of conformal connection function
gammatz_rhs dimensions 3
distribution DEFAULT
group type GF
tags parities={-1 +1 +1 +1 -1 +1 +1 +1 -1} checkpoint=”no”
timelevels 1
variable type REAL
ham centering centering={0 0 0}
hc compact 0
description Hamiltonian constraint
dimensions 3
distribution DEFAULT
group type GF
tags checkpoint=”no”
timelevels 1
variable type REAL
mom centering centering={0 0 0}
mcx compact 0
mcy description Momentum constraints
mcz dimensions 3
distribution DEFAULT
group type GF
tags parities={-1 +1 +1 +1 -1 +1 +1 +1 -1} checkpoint=”no”
timelevels 1
variable type REAL

Adds header:

evolve_utils.hxx

fd_stencils.hxx

Uses header:

loop.hxx

loop_device.hxx

newradx.hxx

5 Schedule

This section lists all the variables which are assigned storage by thorn CanudaX/CanudaX_BSSNMoL. 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: Conditional:
ADMBaseX::metric ADMBaseX::lapse
beta ham
conf_fac_rhs mom
hmetric_rhs  
hcurv_rhs  
trk_rhs  
gammat_rhs  
alpha_rhs  
beta_rhs  
ADMBaseX::shift  
ADMBaseX::curv  
conf_fac  
hmetric  
hcurv  
trk  
gammat  
alpha  
   

Scheduled Functions

CCTK_INITIAL (conditional)

  canudax_bssnmol_initialgroup

  canudax_bssnmol: initial bin

 

  After: tmunubasex_settmunuvars
  Type: group

ODESolvers_PostStep (conditional)

  canudax_bssnmol_poststepgroup

  canudax_bssnmol: post-step operations

 

  Before: admbasex_setadmvars
  Type: group

CanudaX_BSSNMoL_RHSGroup (conditional)

  canudax_bssnmol_apply_boundary_conditions

  apply radiative boundary conditions

 

  After: canudax_bssnmol_calc_bssn_rhs
  Language: c
  Reads: alpha(interior)
    beta(interior)
    conf_fac(interior)
    hmetric(interior)
    trk(interior)
    hcurv(interior)
    gammat(interior)
    alpha_rhs(interior)
    beta_rhs(interior)
    conf_fac_rhs(interior)
    trk_rhs(interior)
    hmetric_rhs(interior)
    hcurv_rhs(interior)
    gammat_rhs(interior)
  Type: function
  Writes: alpha_rhs(interior)
    beta_rhs(interior)
    conf_fac_rhs(interior)
    trk_rhs(interior)
    hmetric_rhs(interior)
    hcurv_rhs(interior)
    gammat_rhs(interior)

CanudaX_BSSNMoL_PostStepGroup (conditional)

  canudax_bssnmol_sync

  sync state vector

 

  Language: c
  Options: global
  Sync: alpha
    beta
    conf_fac
    hmetric
    trk
    hcurv
    gammat
  Type: function

CanudaX_BSSNMoL_PostStepGroup (conditional)

  canudax_bssnmol_impose_conf_fac_floor

  make sure conformal factor does not drop below specified value

 

  Before: canudax_bssnmol_sync
  Language: c
  Reads: conf_fac(interior)
  Type: function
  Writes: conf_fac(interior)

CanudaX_BSSNMoL_PostStepGroup (conditional)

  canudax_bssnmol_reset_detmetric

  reset dethh = 1

 

  After: canudax_bssnmol_impose_conf_fac_floor
  Before: canudax_bssnmol_sync
  Language: c
  Reads: hmetric(interior)
  Type: function
  Writes: hmetric(interior)

CanudaX_BSSNMoL_PostStepGroup (conditional)

  canudax_bssnmol_remove_tra

  remove trace of a

 

  After: canudax_bssnmol_reset_detmetric
  Before: canudax_bssnmol_sync
  Language: c
  Reads: hcurv(interior)
    hmetric(interior)
  Type: function
  Writes: hcurv(interior)
    hmetric(interior)

CanudaX_BSSNMoL_PostStepGroup (conditional)

  canudax_bssnmol_bssn2adm

  convert variables back to the adm ones

 

  After: canudax_bssnmol_sync
  Before: admbasex_setadmvars
  Language: c
  Options: local
  Reads: alpha(everywhere)
    beta(everywhere)
    conf_fac(everywhere)
    hmetric(everywhere)
    trk(everywhere)
    hcurv(everywhere)
    gammat(everywhere)
  Type: function
  Writes: admbasex::metric(everywhere)
    admbasex::curv(everywhere)
    admbasex::lapse(everywhere)
    admbasex::shift(everywhere)

CCTK_ANALYSIS (conditional)

  canudax_bssnmol_bssn_constraints

  compute constraints

 

  Language: c
  Reads: alpha(everywhere)
    beta(everywhere)
    conf_fac(everywhere)
    hmetric(everywhere)
    trk(everywhere)
    hcurv(everywhere)
    gammat(everywhere)
    tmunubasex::ettt(interior)
    tmunubasex::etti(interior)
    tmunubasex::etij(interior)
  Type: function
  Writes: ham(interior)
    mom(interior)

ODESolvers_RHS (conditional)

  canudax_bssnmol_rhsgroup

  canudax_bssnmol: calculate bssn rhs

 

  Type: group

CCTK_PARAMCHECK (conditional)

  canudax_bssnmol_paramcheck

  check canudax_bssnmol parameters for consistency

 

  Language: c
  Type: function

CanudaX_BSSNMoL_InitialGroup (conditional)

  canudax_bssnmol_adm2bssn

  convert initial data into bssn variables

 

  Language: c
  Options: local
  Reads: admbasex::metric(everywhere)
    admbasex::curv(everywhere)
    admbasex::lapse(everywhere)
    admbasex::shift(everywhere)
  Type: function
  Writes: alpha(everywhere)
    beta(everywhere)
    conf_fac(everywhere)
    hmetric(everywhere)
    hcurv(everywhere)
    trk(everywhere)

CanudaX_BSSNMoL_InitialGroup (conditional)

  canudax_bssnmol_gammat

  set gammat initial data

 

  After: canudax_bssnmol_adm2bssn
  Language: c
  Options: local
  Reads: hmetric(everywhere)
  Sync: gammat
  Type: function
  Writes: gammat(interior)

CanudaX_BSSNMoL_InitialGroup (conditional)

  canudax_bssnmol_estimate_error

  estimate error for regridding

 

  After: canudax_bssnmol_sync
  Language: c
  Reads: conf_fac(everywhere)
  Type: function
  Writes: carpetxregrid::regrid_error(interior)

ODESolvers_EstimateError (conditional)

  canudax_bssnmol_estimate_error

  estimate error for regridding

 

  Language: c
  Reads: conf_fac(everywhere)
  Type: function
  Writes: carpetxregrid::regrid_error(interior)

CanudaX_BSSNMoL_RHSGroup (conditional)

  canudax_bssnmol_calc_bssn_rhs

  mol rhs calculation

 

  Language: c
  Reads: alpha(everywhere)
    beta(everywhere)
    conf_fac(everywhere)
    hmetric(everywhere)
    trk(everywhere)
    hcurv(everywhere)
    gammat(everywhere)
    tmunubasex::ettt(interior)
    tmunubasex::etti(interior)
    tmunubasex::etij(interior)
  Type: function
  Writes: alpha_rhs(interior)
    beta_rhs(interior)
    conf_fac_rhs(interior)
    trk_rhs(interior)
    hmetric_rhs(interior)
    hcurv_rhs(interior)
    gammat_rhs(interior)

CanudaX_BSSNMoL_RHSGroup (conditional)

  canudax_bssnmol_add_dissipation

  add ko dissipation

 

  After: canudax_bssnmol_calc_bssn_rhs
  Before: canudax_bssnmol_apply_boundary_conditions
  Language: c
  Reads: alpha(everywhere)
    beta(everywhere)
    conf_fac(everywhere)
    hmetric(everywhere)
    trk(everywhere)
    hcurv(everywhere)
    gammat(everywhere)
    alpha_rhs(interior)
    beta_rhs(interior)
    conf_fac_rhs(interior)
    trk_rhs(interior)
    hmetric_rhs(interior)
    hcurv_rhs(interior)
    gammat_rhs(interior)
  Type: function
  Writes: alpha_rhs(interior)
    beta_rhs(interior)
    conf_fac_rhs(interior)
    trk_rhs(interior)
    hmetric_rhs(interior)
    hcurv_rhs(interior)
    gammat_rhs(interior)