eschnett/cactusamrex - Change ZJVGO4A2JCKZGLPXJWBBUO7TJZCT7AI5RGPS24LQD35T3W5C7FGAC

Punctures: New thorn

Created by Erik Schnetter on June 9, 2020

ZJVGO4A2JCKZGLPXJWBBUO7TJZCT7AI5RGPS24LQD35T3W5C7FGAC

Dependencies

In channels

master

main

Change contents

File addition: Punctures (drwxr-xr-x)
[1.0]

File addition: README (-rw-r--r--)

[0.12]

Cactus Code Thorn Punctures
Author(s)    : Erik Schnetter <schnetter@gmail.com>
Maintainer(s): Erik Schnetter <schnetter@gmail.com>
Licence      : LGPL
--------------------------------------------------------------------------
1. Purpose
Provide multi-black-hole initial conditions based on the "puncture" idea:
Steven Brandt, Bernd Brügmann, "A simple construction of initial data
for multiple black holes", Phys. Rev. Lett. 78:3606-3609, 1997, DOI:
10.1103/PhysRevLett.78.3606, arXiv:gr-qc/9703066.

File addition: configuration.ccl (-rw-r--r--)

[0.12]

# Configuration definitions for thorn Punctures
REQUIRES Arith

File addition: interface.ccl (-rw-r--r--)

[0.12]

# Interface definition for thorn Punctures
IMPLEMENTS: Punctures
INHERITS: ADMBase
USES INCLUDE HEADER: loop.hxx
USES INCLUDE HEADER: mat.hxx
USES INCLUDE HEADER: sum.hxx
USES INCLUDE HEADER: vec.hxx
USES INCLUDE HEADER: vect.hxx
void FUNCTION CallScheduleGroup(
  CCTK_POINTER IN cctkGH,
  CCTK_STRING IN groupname)
REQUIRES FUNCTION CallScheduleGroup
void FUNCTION SolvePoisson(
  CCTK_INT IN gi_sol,
  CCTK_INT IN gi_rhs,
  CCTK_INT IN gi_res,
  CCTK_REAL IN reltol,
  CCTK_REAL IN abstol,
  CCTK_REAL OUT res_initial,
  CCTK_REAL OUT res_final)
REQUIRES FUNCTION SolvePoisson
CCTK_REAL urhs TYPE=gf TAGS='index={0 0 0} checkpoint="no"' "Right hand side"
CCTK_REAL usol TYPE=gf TAGS='index={0 0 0} checkpoint="no"' "Conformal factor u"
CCTK_REAL ures TYPE=gf TAGS='index={0 0 0} checkpoint="no"' "Residual"

File addition: par (drwxr-xr-x)
[0.12]

File addition: brill-lindquist.par (-rw-r--r--)

[0.1461]

# run.me:
# run.cores: 4
# run.memory: 1.0e9
# run.time: 60.0
ActiveThorns = "
    ADMBase
    CarpetX
    Formaline
    IOUtil
    Punctures
"
Cactus::cctk_show_schedule = no
Cactus::presync_mode = "mixed-error"
CarpetX::verbose = no
CarpetX::poison_undefined_values = yes
Cactus::cctk_itlast = 0
CarpetX::xmin = -1.0
CarpetX::ymin = -1.0
CarpetX::zmin = -1.0
CarpetX::xmax = +1.0
CarpetX::ymax = +1.0
CarpetX::zmax = +1.0
CarpetX::ncells_x = 8
CarpetX::ncells_y = 8
CarpetX::ncells_z = 8
CarpetX::periodic_x = no
CarpetX::periodic_y = no
CarpetX::periodic_z = no
CarpetX::periodic = no
# CarpetX::reflection_x = yes
# CarpetX::reflection_y = yes
# CarpetX::reflection_z = yes
# CarpetX::reflection_upper_x = yes
# CarpetX::reflection_upper_y = yes
# CarpetX::reflection_upper_z = yes
CarpetX::ghost_size = 1
ADMBase::initial_data = "Punctures"
ADMBase::initial_lapse = "Punctures"
ADMBase::initial_shift = "Punctures"
ADMBase::initial_dtlapse = "Punctures"
ADMBase::initial_dtshift = "Punctures"
Punctures::npunctures = 1
IO::out_dir = $parfile
IO::out_every = 1
IO::out_mode = "np"
IO::out_proc_every = 1
CarpetX::out_silo_vars = "all"

File addition: headon.par (-rw-r--r--)

[0.1461]

# run.me:
# run.cores: 4
# run.memory: 1.0e9
# run.time: 60.0
ActiveThorns = "
    ADMBase
    CarpetX
    Formaline
    IOUtil
    Punctures
"
Cactus::cctk_show_schedule = no
Cactus::presync_mode = "mixed-error"
CarpetX::verbose = no
CarpetX::poison_undefined_values = yes
Cactus::cctk_itlast = 0
CarpetX::xmin = -1.0
CarpetX::ymin = -1.0
CarpetX::zmin = -1.0
CarpetX::xmax = +1.0
CarpetX::ymax = +1.0
CarpetX::zmax = +1.0
CarpetX::ncells_x = 8
CarpetX::ncells_y = 8
CarpetX::ncells_z = 8
CarpetX::periodic_x = no
CarpetX::periodic_y = no
CarpetX::periodic_z = no
CarpetX::periodic = no
# CarpetX::reflection_x = yes
# CarpetX::reflection_y = yes
# CarpetX::reflection_z = yes
# CarpetX::reflection_upper_x = yes
# CarpetX::reflection_upper_y = yes
# CarpetX::reflection_upper_z = yes
CarpetX::ghost_size = 1
ADMBase::initial_data = "Punctures"
ADMBase::initial_lapse = "Punctures"
ADMBase::initial_shift = "Punctures"
ADMBase::initial_dtlapse = "Punctures"
ADMBase::initial_dtshift = "Punctures"
Punctures::npunctures = 2
Punctures::mass[0] =  0.5
Punctures::posz[0] = +0.5
Punctures::mass[1] =  0.5
Punctures::posz[1] = -0.5
IO::out_dir = $parfile
IO::out_every = 1
IO::out_mode = "np"
IO::out_proc_every = 1
CarpetX::out_silo_vars = "all"

File addition: minkowski.par (-rw-r--r--)

[0.1461]

# run.me:
# run.cores: 4
# run.memory: 1.0e9
# run.time: 60.0
ActiveThorns = "
    ADMBase
    CarpetX
    Formaline
    IOUtil
    Punctures
"
Cactus::cctk_show_schedule = no
Cactus::presync_mode = "mixed-error"
CarpetX::verbose = no
CarpetX::poison_undefined_values = yes
Cactus::cctk_itlast = 0
CarpetX::xmin = -1.0
CarpetX::ymin = -1.0
CarpetX::zmin = -1.0
CarpetX::xmax = +1.0
CarpetX::ymax = +1.0
CarpetX::zmax = +1.0
CarpetX::ncells_x = 8
CarpetX::ncells_y = 8
CarpetX::ncells_z = 8
CarpetX::periodic_x = no
CarpetX::periodic_y = no
CarpetX::periodic_z = no
CarpetX::periodic = no
# CarpetX::reflection_x = yes
# CarpetX::reflection_y = yes
# CarpetX::reflection_z = yes
# CarpetX::reflection_upper_x = yes
# CarpetX::reflection_upper_y = yes
# CarpetX::reflection_upper_z = yes
CarpetX::ghost_size = 1
ADMBase::initial_data = "Punctures"
ADMBase::initial_lapse = "Punctures"
ADMBase::initial_shift = "Punctures"
ADMBase::initial_dtlapse = "Punctures"
ADMBase::initial_dtshift = "Punctures"
IO::out_dir = $parfile
IO::out_every = 1
IO::out_mode = "np"
IO::out_proc_every = 1
CarpetX::out_silo_vars = "all"

File addition: qc0-unigrid.par (-rw-r--r--)

[0.1461]

# run.me:
# run.cores: 4
# run.memory: 1.0e9
# run.time: 60.0
ActiveThorns = "
    ADMBase
    CarpetX
    Formaline
    IOUtil
    Punctures
"
Cactus::cctk_show_schedule = no
Cactus::presync_mode = "mixed-error"
CarpetX::verbose = no
CarpetX::poison_undefined_values = yes
Cactus::cctk_itlast = 0
CarpetX::xmin = -10.0
CarpetX::ymin = -10.0
CarpetX::zmin = -10.0
CarpetX::xmax = +10.0
CarpetX::ymax = +10.0
CarpetX::zmax = +10.0
CarpetX::ncells_x = 64
CarpetX::ncells_y = 64
CarpetX::ncells_z = 64
CarpetX::periodic_x = no
CarpetX::periodic_y = no
CarpetX::periodic_z = no
CarpetX::periodic = no
# CarpetX::reflection_x = yes
# CarpetX::reflection_y = yes
# CarpetX::reflection_z = yes
# CarpetX::reflection_upper_x = yes
# CarpetX::reflection_upper_y = yes
# CarpetX::reflection_upper_z = yes
CarpetX::ghost_size = 1
ADMBase::initial_data = "Punctures"
ADMBase::initial_lapse = "Punctures"
ADMBase::initial_shift = "Punctures"
ADMBase::initial_dtlapse = "Punctures"
ADMBase::initial_dtshift = "Punctures"
# QC-0 setup
Punctures::npunctures = 2
Punctures::mass[0] =  0.453
Punctures::posx[0] = +1.168642873
Punctures::momy[0] = +0.3331917498
Punctures::mass[1] =  0.453
Punctures::posx[1] = -1.168642873
Punctures::momy[1] = -0.3331917498
IO::out_dir = $parfile
IO::out_every = 1
IO::out_mode = "np"
IO::out_proc_every = 1
CarpetX::out_silo_vars = "all"

File addition: qc0.par (-rw-r--r--)

[0.1461]

# run.me:
# run.cores: 40
# run.memory: 1.0e9
# run.time: 600.0
ActiveThorns = "
    ADMBase
    CarpetX
    ErrorEstimator
    Formaline
    IOUtil
    Punctures
"
Cactus::cctk_show_schedule = no
Cactus::presync_mode = "mixed-error"
CarpetX::verbose = no
CarpetX::poison_undefined_values = yes
Cactus::cctk_itlast = 0
CarpetX::xmin = -16.0
CarpetX::ymin = -16.0
CarpetX::zmin = -16.0
CarpetX::xmax = +16.0
CarpetX::ymax = +16.0
CarpetX::zmax = +16.0
CarpetX::ncells_x = 64
CarpetX::ncells_y = 64
CarpetX::ncells_z = 64
CarpetX::periodic_x = no
CarpetX::periodic_y = no
CarpetX::periodic_z = no
CarpetX::periodic = no
# CarpetX::reflection_x = yes
# CarpetX::reflection_y = yes
# CarpetX::reflection_z = yes
# CarpetX::reflection_upper_x = yes
# CarpetX::reflection_upper_y = yes
# CarpetX::reflection_upper_z = yes
CarpetX::ghost_size = 1
CarpetX::max_num_levels = 5
CarpetX::regrid_every = 0
CarpetX::regrid_error_threshold = 8.0
ErrorEstimator::region_shape = "cube"
ErrorEstimator::scale_by_resolution = yes
ADMBase::initial_data = "Punctures"
ADMBase::initial_lapse = "Punctures"
ADMBase::initial_shift = "Punctures"
ADMBase::initial_dtlapse = "Punctures"
ADMBase::initial_dtshift = "Punctures"
# QC-0 setup
Punctures::npunctures = 2
Punctures::mass[0] =  0.453
Punctures::posx[0] = +1.168642873
Punctures::momy[0] = +0.3331917498
Punctures::mass[1] =  0.453
Punctures::posx[1] = -1.168642873
Punctures::momy[1] = -0.3331917498
IO::out_dir = $parfile
IO::out_every = 1
IO::out_mode = "np"
IO::out_proc_every = 1
CarpetX::out_silo_vars = "all"

File addition: param.ccl (-rw-r--r--)

[0.12]

# Parameter definitions for thorn Punctures
SHARES: ADMBase
EXTENDS KEYWORD initial_data "Initial metric and extrinsic curvature datasets"
{
  "Punctures" :: "Puncture black holes"
}
EXTENDS KEYWORD initial_lapse "Initial lapse value"
{
  "Punctures" :: "Lapse suitable for puncture black holes"
}
EXTENDS KEYWORD initial_shift "Initial shift value"
{
  "Punctures" :: "Shift suitable for puncture black holes"
}
EXTENDS KEYWORD initial_dtlapse "Initial dtlapse value"
{
  "Punctures" :: "Lapse time derivative suitable for puncture black holes"
}
EXTENDS KEYWORD initial_dtshift "Initial dtshift value"
{
  "Punctures" :: "Shift time derivative suitable for puncture black holes"
}
PRIVATE:
CCTK_INT npunctures "Number of punctures"
{
  0:11 :: ""
} 0
CCTK_REAL mass[11] "Puncture mass"
{
  0.0:* :: ""
} 1.0
CCTK_REAL posx[11] "Puncture position"
{
  *:* :: ""
} 0.0
CCTK_REAL posy[11] "Puncture position"
{
  *:* :: ""
} 0.0
CCTK_REAL posz[11] "Puncture position"
{
  *:* :: ""
} 0.0
CCTK_REAL momx[11] "Puncture momentum"
{
  *:* :: ""
} 0.0
CCTK_REAL momy[11] "Puncture momentum"
{
  *:* :: ""
} 0.0
CCTK_REAL momz[11] "Puncture momentum"
{
  *:* :: ""
} 0.0
CCTK_REAL amomx[11] "Puncture angular momentum"
{
  *:* :: ""
} 0.0
CCTK_REAL amomy[11] "Puncture angular momentum"
{
  *:* :: ""
} 0.0
CCTK_REAL amomz[11] "Puncture angular momentum"
{
  *:* :: ""
} 0.0
CCTK_REAL rmin "Minimum radius (to avoid singularities)"
{
  0.0:* :: ""
} 1.0e-8
CCTK_REAL abstol "Absolute tolerance for non-linear solver"
{
  0.0:* :: ""
} 1.0e-6

File addition: schedule.ccl (-rw-r--r--)

[0.12]

# Schedule definitions for thorn Punctures
SCHEDULE Punctures_init AT initial
{
  LANG: C
  WRITES: usol(everywhere)
} "Set up initial guess"
SCHEDULE Punctures_solve AT initial AFTER  Punctures_init
{
  LANG: C
  OPTIONS: level
  READS: usol(everywhere)
  WRITES: usol(everywhere)
  WRITES: urhs(everywhere)
  WRITES: ures(interior)
} "Solve Hamiltonian constraint"
  SCHEDULE GROUP Punctures_solve1
  {
  } "Perform one linear solve"
  
    SCHEDULE Punctures_rhs IN Punctures_solve1
    {
      LANG: C
      READS: usol(everywhere)
      WRITES: urhs(everywhere)
    } "Set up right hand side"
  
  SCHEDULE GROUP Punctures_solve2
  {
  } "Perform one linear solve"
  
    SCHEDULE Punctures_boundary IN Punctures_solve2
    {
      LANG: C
      READS: usol(everywhere)
      WRITES: usol(everywhere)
    } "Apply boundary conditions to puncture solution"
SCHEDULE Punctures_ADMBase AT initial AFTER Punctures_solve
{
  LANG: C
  READS: usol(everywhere)
  WRITES: ADMBase::metric(everywhere)
  WRITES: ADMBase::curv(everywhere)
  WRITES: ADMBase::lapse(everywhere)
  WRITES: ADMBase::shift(everywhere)
  WRITES: ADMBase::dtlapse(everywhere)
  WRITES: ADMBase::dtshift(everywhere)
} "Set ADMBase variables"

File addition: src (drwxr-xr-x)
[0.12]

File addition: make.code.defn (-rw-r--r--)

[0.10831]

# Main make.code.defn file for thorn Punctures
# Source files in this directory
SRCS = punctures.cxx
# Subdirectories containing source files
SUBDIRS =

File addition: punctures.cxx (-rw-r--r--)

[0.10831]

#include <loop.hxx>
#include <mat.hxx>
#include <sum.hxx>
#include <vec.hxx>
#include <vect.hxx>
#include <cctk.h>
#include <cctk_Arguments_Checked.h>
#include <cctk_Parameters.h>
#include <cassert>
#include <cmath>
namespace Punctures {
using namespace Arith;
using namespace Loop;
using namespace std;
////////////////////////////////////////////////////////////////////////////////
template <typename T> T pown(T x, int n) {
  if (n < 0) {
    x = 1 / x;
    n = -n;
  }
  T r = 1;
  while (n > 0) {
    if (n % 2)
      r *= x;
    x *= x;
    n /= 2;
  }
  return r;
}
const mat<CCTK_REAL, 3, DN, DN> g([](int a, int b) { return a == b; });
const auto epsilon = [](int a, int b,
                        int c) -> CCTK_REAL CCTK_ATTRIBUTE_ALWAYS_INLINE {
  if (a == 0 && b == 1 && c == 2)
    return 1;
  if (a == 0 && b == 2 && c == 1)
    return -1;
  if (a == 1 && b == 0 && c == 2)
    return -1;
  if (a == 1 && b == 2 && c == 0)
    return 1;
  if (a == 2 && b == 0 && c == 1)
    return 1;
  if (a == 2 && b == 1 && c == 0)
    return -1;
  return 0;
};
////////////////////////////////////////////////////////////////////////////////
template <typename T>
CCTK_ATTRIBUTE_ALWAYS_INLINE void fcalc(const PointDesc &p, const T &u,
                                        T &alpha, mat<T, 3, DN, DN> &K, T &rhs,
                                        T &psi) {
  DECLARE_CCTK_PARAMETERS;
  T alpha1 = 0;
  K = mat<T, 3, DN, DN>();
  for (int i = 0; i < npunctures; ++i) {
    const vec<T, 3, UP> x{posx[i], posy[i], posz[i]};
    const vec<T, 3, UP> P{momx[i], momy[i], momz[i]};
    const vec<T, 3, UP> S{amomx[i], amomy[i], amomz[i]};
    const T r = sqrt(sum<3>([&](int a) { return pow(p.X[a] - x(a), 2); }));
    const vec<T, 3, UP> n([&](int a) { return r < rmin ? 0 : x(a) / r; });
    alpha1 += mass[i] / (2 * fmax(rmin, r));
    K += mat<T, 3, DN, DN>([&](int a, int b) {
      return 3 / (2 * pow(fmax(rmin, r), 2)) *
                 (P(a) * n(b) + P(b) * n(a) -
                  (g(a, b) - n(a) * n(b)) *
                      sum<3>([&](int c) { return P(c) * n(c); })) +
             3 / pow(fmax(rmin, r), 3) * sum<3>([&](int c, int d) {
               return epsilon(a, c, d) * S(c) * n(d) * n(b) +
                      epsilon(b, c, d) * S(c) * n(d) * n(a);
             });
    });
  }
  alpha = 1 / alpha1;
  if (isinf1(alpha)) {
    // Infinitely far away, or there are no black holes
    rhs = 0;
    psi = 1;
    return;
  }
  if (alpha == 0) {
    assert(0); // handled by rmin
    // At a puncture
    K = mat<T, 3, DN, DN>();
    rhs = 0;
    psi = INFINITY;
    return;
  }
  const T beta = 1 / T(8) * pow(alpha, 7) *
                 sum<3>([&](int a, int b) { return K(a, b) * K(a, b); });
  rhs = -beta * pow(1 + alpha * u, -7);
  psi = 1 / alpha + u;
}
template <typename T> T frhs(const PointDesc &p, const T &u) {
  T alpha;
  mat<T, 3, DN, DN> K;
  T rhs;
  T psi;
  fcalc(p, u, alpha, K, rhs, psi);
  return rhs;
}
// u = 1 + C / r + ...
template <typename T> T fbnd(const PointDesc &p) { return 1; }
template <typename T> T fpsi(const PointDesc &p, const T &u) {
  T alpha;
  mat<T, 3, DN, DN> K;
  T rhs;
  T psi;
  fcalc(p, u, alpha, K, rhs, psi);
  return psi;
}
template <typename T> mat<T, 3, DN, DN> fK(const PointDesc &p, const T &u) {
  T alpha;
  mat<T, 3, DN, DN> K;
  T rhs;
  T psi;
  fcalc(p, u, alpha, K, rhs, psi);
  return K;
}
////////////////////////////////////////////////////////////////////////////////
extern "C" void Punctures_init(CCTK_ARGUMENTS) {
  DECLARE_CCTK_ARGUMENTS_Punctures_init;
  DECLARE_CCTK_PARAMETERS;
  const GF3D<CCTK_REAL, 0, 0, 0> usol_(cctkGH, usol);
  // Initialize all points (including ghosts)
  loop_all<0, 0, 0>(cctkGH, [&](const PointDesc &p) { usol_(p.I) = 1.0; });
  // Note: The boundary conditions are applied on the outermost layer
  // of the interior points, not on the boundary points.
  // Set boundary conditions
  loop_all<0, 0, 0>(cctkGH, [&](const PointDesc &p) {
    if (cctk_lbnd[0] + p.i <= 1 || cctk_lbnd[0] + p.i >= cctk_gsh[0] - 1 ||
        cctk_lbnd[1] + p.j <= 1 || cctk_lbnd[1] + p.j >= cctk_gsh[1] - 1 ||
        cctk_lbnd[2] + p.k <= 1 || cctk_lbnd[2] + p.k >= cctk_gsh[2] - 1)
      usol_(p.I) = fbnd<CCTK_REAL>(p);
  });
}
extern "C" void Punctures_solve(CCTK_ARGUMENTS) {
  DECLARE_CCTK_ARGUMENTS_Punctures_solve;
  DECLARE_CCTK_PARAMETERS;
  const int max_iters = 100;
  for (int iter = 1;; ++iter) {
    CCTK_VINFO("Nonlinear iteration #%d", iter);
    // Set up RHS
    CallScheduleGroup(cctkGH, "Punctures_solve1");
    const int gi_sol = CCTK_GroupIndex("Punctures::usol");
    assert(gi_sol >= 0);
    const int gi_rhs = CCTK_GroupIndex("Punctures::urhs");
    assert(gi_rhs >= 0);
    const int gi_res = CCTK_GroupIndex("Punctures::ures");
    assert(gi_res >= 0);
    // linear solver accuracy
    const CCTK_REAL reltol = 0.0;
    const CCTK_REAL abstol = 1.0e-10;
    CCTK_REAL res_initial, res_final;
    SolvePoisson(gi_sol, gi_rhs, gi_res, reltol, abstol, &res_initial,
                 &res_final);
    CCTK_VINFO("Linear residual before solve: %g", double(res_initial));
    CCTK_VINFO("Linear residual after solve:  %g", double(res_final));
    // Correct boundaries
    CallScheduleGroup(cctkGH, "Punctures_solve2");
    if (res_initial <= abstol) {
      CCTK_VINFO("Nonlinear iterations finished after %d iterations; accuracy "
                 "goal reached",
                 iter);
      break;
    }
    if (iter == max_iters) {
      CCTK_VWARN(CCTK_WARN_ALERT,
                 "Nonlinear iterations finished after %d iterations; accuracy "
                 "goal NOT reached",
                 iter);
      break;
    }
  }
}
extern "C" void Punctures_rhs(CCTK_ARGUMENTS) {
  DECLARE_CCTK_ARGUMENTS_Punctures_rhs;
  DECLARE_CCTK_PARAMETERS;
  const GF3D<const CCTK_REAL, 0, 0, 0> usol_(cctkGH, usol);
  const GF3D<CCTK_REAL, 0, 0, 0> urhs_(cctkGH, urhs);
  // Set RHS
  loop_all<0, 0, 0>(
      cctkGH, [&](const PointDesc &p) { urhs_(p.I) = frhs(p, usol_(p.I)); });
}
extern "C" void Punctures_boundary(CCTK_ARGUMENTS) {
  DECLARE_CCTK_ARGUMENTS_Punctures_boundary;
  DECLARE_CCTK_PARAMETERS;
  const GF3D<CCTK_REAL, 0, 0, 0> usol_(cctkGH, usol);
  // Set boundary conditions
  loop_all<0, 0, 0>(cctkGH, [&](const PointDesc &p) {
    if (cctk_lbnd[0] + p.i <= 1 || cctk_lbnd[0] + p.i >= cctk_gsh[0] - 1 ||
        cctk_lbnd[1] + p.j <= 1 || cctk_lbnd[1] + p.j >= cctk_gsh[1] - 1 ||
        cctk_lbnd[2] + p.k <= 1 || cctk_lbnd[2] + p.k >= cctk_gsh[2] - 1)
      usol_(p.I) = fbnd<CCTK_REAL>(p);
  });
}
extern "C" void Punctures_ADMBase(CCTK_ARGUMENTS) {
  DECLARE_CCTK_ARGUMENTS_Punctures_ADMBase;
  DECLARE_CCTK_PARAMETERS;
  const GF3D<const CCTK_REAL, 0, 0, 0> usol_(cctkGH, usol);
  const GF3D<CCTK_REAL, 0, 0, 0> gxx_(cctkGH, gxx);
  const GF3D<CCTK_REAL, 0, 0, 0> gxy_(cctkGH, gxy);
  const GF3D<CCTK_REAL, 0, 0, 0> gxz_(cctkGH, gxz);
  const GF3D<CCTK_REAL, 0, 0, 0> gyy_(cctkGH, gyy);
  const GF3D<CCTK_REAL, 0, 0, 0> gyz_(cctkGH, gyz);
  const GF3D<CCTK_REAL, 0, 0, 0> gzz_(cctkGH, gzz);
  const GF3D<CCTK_REAL, 0, 0, 0> Kxx_(cctkGH, kxx);
  const GF3D<CCTK_REAL, 0, 0, 0> Kxy_(cctkGH, kxy);
  const GF3D<CCTK_REAL, 0, 0, 0> Kxz_(cctkGH, kxz);
  const GF3D<CCTK_REAL, 0, 0, 0> Kyy_(cctkGH, kyy);
  const GF3D<CCTK_REAL, 0, 0, 0> Kyz_(cctkGH, kyz);
  const GF3D<CCTK_REAL, 0, 0, 0> Kzz_(cctkGH, kzz);
  const GF3D<CCTK_REAL, 0, 0, 0> alp_(cctkGH, alp);
  const GF3D<CCTK_REAL, 0, 0, 0> betax_(cctkGH, betax);
  const GF3D<CCTK_REAL, 0, 0, 0> betay_(cctkGH, betay);
  const GF3D<CCTK_REAL, 0, 0, 0> betaz_(cctkGH, betaz);
  const GF3D<CCTK_REAL, 0, 0, 0> dtalp_(cctkGH, dtalp);
  const GF3D<CCTK_REAL, 0, 0, 0> dtbetax_(cctkGH, dtbetax);
  const GF3D<CCTK_REAL, 0, 0, 0> dtbetay_(cctkGH, dtbetay);
  const GF3D<CCTK_REAL, 0, 0, 0> dtbetaz_(cctkGH, dtbetaz);
  loop_all<0, 0, 0>(cctkGH, [&](const PointDesc &p) {
    const CCTK_REAL psi = fpsi(p, usol_(p.I));
    const mat<CCTK_REAL, 3, DN, DN> K = fK(p, usol_(p.I));
    const mat<CCTK_REAL, 3, DN, DN> gph = pow(psi, 4) * g;
    const mat<CCTK_REAL, 3, DN, DN> Kph = pow(psi, -2) * K;
    gxx_(p.I) = gph(0, 0);
    gxy_(p.I) = gph(0, 1);
    gxz_(p.I) = gph(0, 2);
    gyy_(p.I) = gph(1, 1);
    gyz_(p.I) = gph(1, 2);
    gzz_(p.I) = gph(2, 2);
    Kxx_(p.I) = Kph(0, 0);
    Kxy_(p.I) = Kph(0, 1);
    Kxz_(p.I) = Kph(0, 2);
    Kyy_(p.I) = Kph(1, 1);
    Kyz_(p.I) = Kph(1, 2);
    Kzz_(p.I) = Kph(2, 2);
    alp_(p.I) = 1; // TODO
    betax_(p.I) = 0;
    betay_(p.I) = 0;
    betaz_(p.I) = 0;
    dtalp_(p.I) = 0;
    dtbetax_(p.I) = 0;
    dtbetay_(p.I) = 0;
    dtbetaz_(p.I) = 0;
  });
}
} // namespace Punctures

Insertion in azure-pipelines/carpetx.th at line 149 [4.11]
[2.8081]
[3.519]
```
CarpetX/Punctures
```

Punctures: New thorn

Dependencies

In channels

Change contents

File addition: Punctures (drwxr-xr-x)

File addition: README (-rw-r--r--)

File addition: configuration.ccl (-rw-r--r--)

File addition: interface.ccl (-rw-r--r--)

File addition: par (drwxr-xr-x)

File addition: brill-lindquist.par (-rw-r--r--)

File addition: headon.par (-rw-r--r--)

File addition: minkowski.par (-rw-r--r--)

File addition: qc0-unigrid.par (-rw-r--r--)

File addition: qc0.par (-rw-r--r--)

File addition: param.ccl (-rw-r--r--)

File addition: schedule.ccl (-rw-r--r--)

File addition: src (drwxr-xr-x)

File addition: make.code.defn (-rw-r--r--)

File addition: punctures.cxx (-rw-r--r--)

Insertion in azure-pipelines/carpetx.th at line 149 [4.11]