eschnett/cactusamrex - Change LZG2DA6SVFAC6E2EBLGPZLV423AME7FEBJ7SZP5OYFUYTO6HCAVQC

DGWaveToy: Implement WaveToy with DG derivatives

Created by Erik Schnetter on June 27, 2020

LZG2DA6SVFAC6E2EBLGPZLV423AME7FEBJ7SZP5OYFUYTO6HCAVQC

Dependencies

In channels

main

Change contents

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

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

[0.16]

Cactus Code Thorn DGCoordinates
Author(s)    : Erik Schnetter <schnetter@gmail.com>
Maintainer(s): Erik Schnetter <schnetter@gmail.com>
Licence      : LGPL
--------------------------------------------------------------------------
1. Purpose
Provide coordinate grid functions

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

[0.16]

# Interface definition for thorn DGCoordinates
IMPLEMENTS: DGCoordinates
INCLUDES HEADER: dg.hxx IN dg.hxx
USES INCLUDE HEADER: dg.hxx
USES INCLUDE HEADER: div.hxx
USES INCLUDE HEADER: loop.hxx
USES INCLUDE HEADER: vect.hxx
PUBLIC:
CCTK_REAL coords TYPE=gf { coordx coordy coordz } "Coordinates"
CCTK_REAL volume TYPE=gf { cvol } "Cell volume"

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

[0.16]

# Parameter definitions for thorn DGCoordinates
RESTRICTED:
CCTK_INT dg_npoints "Discontinuous Galerkin number of collocation points"
{
  2 :: ""
  3 :: ""
  4 :: ""
  8 :: ""
} 2

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

[0.16]

# Schedule definitions for thorn DGCoordinates
SCHEDULE DGCoordinates_Setup AT basegrid
{
  LANG: C
  WRITES: coords(everywhere)
  WRITES: volume(everywhere)
} "Set coordinate grid functions"

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

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

[0.1086]

#include "dg.hxx"
#include <div.hxx>
#include <loop.hxx>
#include <vect.hxx>
#include <cctk.h>
#include <cctk_Arguments_Checked.h>
#include <cctk_Parameters.h>
#include <cassert>
#include <cstdlib>
namespace DGCoordinates {
using namespace Arith;
using namespace DG;
using namespace std;
using Loop::dim;
extern "C" void DGCoordinates_Setup(CCTK_ARGUMENTS) {
  DECLARE_CCTK_ARGUMENTS_DGCoordinates_Setup;
  DECLARE_CCTK_PARAMETERS;
  const Loop::GF3D<CCTK_REAL, 1, 1, 1> coordx_(cctkGH, coordx);
  const Loop::GF3D<CCTK_REAL, 1, 1, 1> coordy_(cctkGH, coordy);
  const Loop::GF3D<CCTK_REAL, 1, 1, 1> coordz_(cctkGH, coordz);
  const Loop::GF3D<CCTK_REAL, 1, 1, 1> cvol_(cctkGH, cvol);
  const Loop::GridDescBase grid(cctkGH);
  for (int d = 0; d < dim; ++d)
    assert(grid.nghostzones[d] == 1);
  for (int d = 0; d < dim; ++d)
    assert(mod_floor(grid.gsh[d] - 2 * grid.nghostzones[d], dg_npoints) == 0);
  for (int d = 0; d < dim; ++d)
    assert(mod_floor(grid.lbnd[d], dg_npoints) == 0);
  for (int d = 0; d < dim; ++d)
    assert(mod_floor(grid.lsh[d] - 2 * grid.nghostzones[d], dg_npoints) == 0);
  switch (dg_npoints) {
  case 2: {
    constexpr int N = 2;
    grid.loop_all<1, 1, 1>(grid.nghostzones, [&](const Loop::PointDesc &p) {
      const get_coords<CCTK_REAL, N> coords(grid, p);
      coordx_(p.I) = coords.coord[0];
      coordy_(p.I) = coords.coord[1];
      coordz_(p.I) = coords.coord[2];
      cvol_(p.I) = coords.vol;
    });
    break;
  }
  case 3: {
    constexpr int N = 3;
    grid.loop_all<1, 1, 1>(grid.nghostzones, [&](const Loop::PointDesc &p) {
      const get_coords<CCTK_REAL, N> coords(grid, p);
      coordx_(p.I) = coords.coord[0];
      coordy_(p.I) = coords.coord[1];
      coordz_(p.I) = coords.coord[2];
      cvol_(p.I) = coords.vol;
    });
    break;
  }
  case 4: {
    constexpr int N = 4;
    grid.loop_all<1, 1, 1>(grid.nghostzones, [&](const Loop::PointDesc &p) {
      const get_coords<CCTK_REAL, N> coords(grid, p);
      coordx_(p.I) = coords.coord[0];
      coordy_(p.I) = coords.coord[1];
      coordz_(p.I) = coords.coord[2];
      cvol_(p.I) = coords.vol;
    });
    break;
  }
  case 8: {
    constexpr int N = 8;
    grid.loop_all<1, 1, 1>(grid.nghostzones, [&](const Loop::PointDesc &p) {
      const get_coords<CCTK_REAL, N> coords(grid, p);
      coordx_(p.I) = coords.coord[0];
      coordy_(p.I) = coords.coord[1];
      coordz_(p.I) = coords.coord[2];
      cvol_(p.I) = coords.vol;
    });
    break;
  }
  default:
    assert(0);
  }
}
} // namespace DGCoordinates

File addition: dg.hxx (-rw-r--r--)

[0.1086]

#ifndef DG_HXX
#define DG_HXX
#include <div.hxx>
#include <loop.hxx>
#include <vect.hxx>
#include <array>
namespace DG {
using namespace Arith;
using namespace Loop;
using namespace std;
template <typename T, int N> struct dg_basis;
// Collocation points are in the range [-1; 1].
// Derivative coefficients include 1 ghost point.
// Note that the derivative operator is non-zero for ghosts only for the
// boundary points, and is zero on the diagonal.
template <typename T> struct dg_basis<T, 2> {
  static constexpr int N = 2;
  static constexpr array<T, N> coords{-1, 1};
  static constexpr array<T, N> vols{0.5, 1, 0.5};
  typedef array<T, N + 2> ATN2;
  static constexpr array<ATN2, N> deriv{ATN2{-0.5, 0, 0.5, 0},
                                        ATN2{0, -0.5, 0, 0.5}};
};
template <typename T> struct dg_basis<T, 3> {
  static constexpr int N = 3;
  static constexpr array<T, N> coords{-1, 0, 0};
  static constexpr array<T, N> vols{0.5, 1, 0.5};
  typedef array<T, N + 2> ATN2;
  static constexpr array<ATN2, N> deriv{ATN2{-1.5, 0, 2, -0.5, 0}, //
                                        ATN2{0, -0.5, 0, 0.5, 0},  //
                                        ATN2{0, 0.5, -2, 0, 1.5}};
};
template <typename T> struct dg_basis<T, 4> {
  static constexpr int N = 4;
  static constexpr array<T, N> coords{
      -1, -0.4472135954999579392818347337462552470881,
      0.4472135954999579392818347337462552470881, 1};
  static constexpr array<T, N> vols{
      0.166666666666666666666666666667, 0.83333333333333333333333333333,
      0.83333333333333333333333333333, 0.166666666666666666666666666667};
  typedef array<T, N + 2> ATN2;
  static constexpr array<ATN2, N> deriv{
      ATN2{-3, 0, 4.0450849718747371205114670859,
           -1.54508497187473712051146708591, 0.500000000000000000000000000000,
           0},
      ATN2{0, -0.80901699437494742410229341718, 0,
           1.1180339887498948482045868344, -0.309016994374947424102293417183,
           0},
      ATN2{0, 0.309016994374947424102293417183, -1.1180339887498948482045868344,
           0, 0.80901699437494742410229341718, 0},
      ATN2{0, -0.500000000000000000000000000000,
           1.54508497187473712051146708591, -4.0450849718747371205114670859, 0,
           3}};
};
template <typename T> struct dg_basis<T, 8> {
  static constexpr int N = 8;
  static constexpr array<T, N> coords{
      -1,
      -0.8717401485096066153374457612206634381038,
      -0.5917001814331423021445107313979531899457,
      -0.2092992179024788687686572603453512552955,
      0.2092992179024788687686572603453512552955,
      0.5917001814331423021445107313979531899457,
      0.8717401485096066153374457612206634381038,
      1};
  static constexpr array<T, N> vols{
      0.035714285714285714285714285714, 0.210704227143506039382992065776,
      0.34112269248350436476424067711,  0.41245879465870388156705297140,
      0.41245879465870388156705297140,  0.34112269248350436476424067711,
      0.210704227143506039382992065776, 0.035714285714285714285714285714};
  typedef array<T, N + 2> ATN2;
  static constexpr array<ATN2, N> deriv{
      ATN2{-14, 0, 18.9375986071173705131546424941,
           -7.5692898193484870087981787234, 4.2979081642651752142490066598,
           -2.810188989257949038537763687, 1.94165942554412230232935244532,
           -1.29768738832023198239705918927, 0.50000000000000000000000000000,
           0},
      ATN2{0, -3.209915703002990336815351572, 0, 4.5435850645665642173963830608,
           -2.11206121431454222839296010401, 1.2942320509135015076932056372,
           -0.8694480983314929341612144455, 0.5735654149402641373099571270,
           -0.21995751477130436303001970391, 0},
      ATN2{0, 0.7924766813205145268013597361, -2.806475794736433436470922316, 0,
           2.875517405972505217473636209, -1.37278583180602848091037564048,
           0.84502255650651048982986561020, -0.5370395861576610609376193431,
           0.203284568900592744214055743889, 0},
      ATN2{0, -0.3721504357285948658470918627, 1.0789446887904527025084635997,
           -2.378187233515505824569788251, 0, 2.388924359158239214879541510,
           -1.1353580168811114404647212743, 0.6611573509003112232014379102,
           -0.2433307127237910097078416316, 0},
      ATN2{0, 0.2433307127237910097078416316, -0.6611573509003112232014379102,
           1.1353580168811114404647212743, -2.388924359158239214879541510, 0,
           2.378187233515505824569788251, -1.0789446887904527025084635997,
           0.3721504357285948658470918627, 0},
      ATN2{0, -0.203284568900592744214055743889, 0.5370395861576610609376193431,
           -0.84502255650651048982986561020, 1.37278583180602848091037564048,
           -2.875517405972505217473636209, 0, 2.806475794736433436470922316,
           -0.7924766813205145268013597361, 0},
      ATN2{0, 0.21995751477130436303001970391, -0.5735654149402641373099571270,
           0.8694480983314929341612144455, -1.2942320509135015076932056372,
           2.11206121431454222839296010401, -4.5435850645665642173963830608, 0,
           3.209915703002990336815351572, 0},
      ATN2{0, -0.50000000000000000000000000000, 1.29768738832023198239705918927,
           -1.94165942554412230232935244532, 2.810188989257949038537763687,
           -4.2979081642651752142490066598, 7.5692898193484870087981787234,
           -18.9375986071173705131546424941, 0, 14}};
};
////////////////////////////////////////////////////////////////////////////////
template <typename T, int N> struct get_coords {
  vect<T, dim> coord;
  T vol;
  get_coords(const Loop::GridDescBase &grid, const Loop::PointDesc &p) {
    vect<int, dim> i0, i1;
    vect<T, dim> x0, dx;
    for (int d = 0; d < dim; ++d) {
      i1[d] = mod_floor(p.I[d] - grid.nghostzones[d], N);
      i0[d] = p.I[d] - i1[d];
      x0[d] = grid.x0[d] + (grid.lbnd[d] + i0[d] - T(0.5)) * grid.dx[d];
      dx[d] = N * grid.dx[d];
      coord[d] = x0[d] + dx[d] / 2 + dx[d] / 2 * dg_basis<T, N>().coords[i1[d]];
    }
    vol = 1;
    for (int d = 0; d < dim; ++d)
      vol *= dx[d] / (2 * N) * dg_basis<T, N>().coords[i1[d]];
  }
};
template <int N, typename T>
vect<T, dim> get_derivs(const Loop::GridDescBase &grid,
                        const Loop::GF3D<const T, 1, 1, 1> &u,
                        const Loop::PointDesc &p) {
  vect<T, dim> deriv;
  vect<int, dim> i0, i1;
  vect<T, dim> x0, dx;
  for (int d = 0; d < dim; ++d) {
    i1[d] = mod_floor(p.I[d] - grid.nghostzones[d], int(N));
    i0[d] = p.I[d] - i1[d];
    x0[d] = grid.x0[d] + (grid.lbnd[d] + i0[d]) * grid.dx[d];
    dx[d] = N * grid.dx[d];
    deriv[d] = 0;
    for (int j = 0; j < int(N) + 2; ++j) {
      vect<int, dim> J = p.I;
      J[d] = i0[d] + j - 1;
      deriv[d] += dg_basis<T, N>().deriv[i1[d]][j] * u(J);
    }
    deriv[d] /= dx[d] / 2;
  }
  return deriv;
}
} // namespace DG
#endif // #ifndef DG_HXX

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

[0.1086]

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

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

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

[0.10723]

Cactus Code Thorn DGWaveToy
Author(s)    : Erik Schnetter <schnetter@gmail.com>
Maintainer(s): Erik Schnetter <schnetter@gmail.com>
Licence      : LGPL
--------------------------------------------------------------------------
1. Purpose
Solve the scalar wave equation
2. Equations
dt^2 u = dx^2 u
ft = dt u
fx = dx u
dt u = ft
dt ft = dx fx
dt fx = dx ft

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

[0.10723]

# Interface definition for thorn DGWaveToy
IMPLEMENTS: DGWaveToy
INHERITS: DGCoordinates
USES INCLUDE HEADER: dg.hxx
USES INCLUDE HEADER: loop.hxx
USES INCLUDE HEADER: vect.hxx
PUBLIC:
CCTK_REAL u TYPE=gf TAGS='parities={1 1 1} rhs="u_rhs"' { u } "Scalar"
CCTK_REAL f TYPE=gf TAGS='parities={1 1 1   -1 1 1   1 -1 1   1 1 -1} rhs="f_rhs"' { ft fx fy fz } "Flux"
CCTK_REAL u_rhs TYPE=gf TAGS='parities={1 1 1} checkpoint="no"' { u_rhs } "RHS for scalar"
CCTK_REAL f_rhs TYPE=gf TAGS='parities={1 1 1   -1 1 1   1 -1 1   1 1 -1} checkpoint="no"' { ft_rhs fx_rhs fy_rhs fz_rhs } "RHS for flux"
CCTK_REAL eps TYPE=gf TAGS='parities={1 1 1} checkpoint="no"' { eps } "Energy density"
CCTK_REAL u_error TYPE=gf TAGS='parities={1 1 1} checkpoint="no"' { u_error } "Error of scalar"
CCTK_REAL f_error TYPE=gf TAGS='parities={1 1 1   -1 1 1   1 -1 1   1 1 -1} checkpoint="no"' { ft_error fx_error fy_error fz_error } "Error of flux"

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

File addition: standing-128.par (-rw-r--r--)

[0.12054]

ActiveThorns = "
    CarpetX
    DGCoordinates
    DGWaveToy
    Formaline
    IOUtil
    ODESolvers
"
$ncells = 128
Cactus::cctk_show_schedule = yes
Cactus::presync_mode = "mixed-error"
Cactus::terminate = "time"
Cactus::cctk_final_time = 1.0
CarpetX::verbose = yes
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 = $ncells
CarpetX::ncells_y = $ncells
CarpetX::ncells_z = $ncells
CarpetX::periodic_x = no
CarpetX::periodic_y = no
CarpetX::periodic_z = no
CarpetX::periodic = yes
CarpetX::ghost_size = 1
CarpetX::max_num_levels = 1
DGCoordinates::dg_npoints = 2
ODESolvers::method = "RK2"
CarpetX::dtfac = 0.25
IO::out_dir = $parfile
IO::out_every = $ncells
CarpetX::out_silo_vars = ""   # "all"

File addition: standing-64.par (-rw-r--r--)

[0.12054]

ActiveThorns = "
    CarpetX
    DGCoordinates
    DGWaveToy
    Formaline
    IOUtil
    ODESolvers
"
$ncells = 64
Cactus::cctk_show_schedule = yes
Cactus::presync_mode = "mixed-error"
Cactus::terminate = "time"
Cactus::cctk_final_time = 1.0
CarpetX::verbose = yes
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 = $ncells
CarpetX::ncells_y = $ncells
CarpetX::ncells_z = $ncells
CarpetX::periodic_x = no
CarpetX::periodic_y = no
CarpetX::periodic_z = no
CarpetX::periodic = yes
CarpetX::ghost_size = 1
CarpetX::max_num_levels = 1
DGCoordinates::dg_npoints = 2
ODESolvers::method = "RK2"
CarpetX::dtfac = 0.25
IO::out_dir = $parfile
IO::out_every = $ncells
CarpetX::out_silo_vars = ""   # "all"

File addition: standing-n4-64.par (-rw-r--r--)

[0.12054]

ActiveThorns = "
    CarpetX
    DGCoordinates
    DGWaveToy
    Formaline
    IOUtil
    ODESolvers
"
$ncells = 64
Cactus::cctk_show_schedule = yes
Cactus::presync_mode = "mixed-error"
Cactus::terminate = "time"
Cactus::cctk_final_time = 1.0
CarpetX::verbose = yes
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 = $ncells
CarpetX::ncells_y = $ncells
CarpetX::ncells_z = $ncells
CarpetX::periodic_x = no
CarpetX::periodic_y = no
CarpetX::periodic_z = no
CarpetX::periodic = yes
CarpetX::ghost_size = 1
CarpetX::max_num_levels = 1
DGCoordinates::dg_npoints = 4
ODESolvers::method = "RK4"
CarpetX::dtfac = 0.25
IO::out_dir = $parfile
IO::out_every = $ncells
CarpetX::out_silo_vars = ""   # "all"

File addition: standing-n4.par (-rw-r--r--)

[0.12054]

ActiveThorns = "
    CarpetX
    DGCoordinates
    DGWaveToy
    Formaline
    IOUtil
    ODESolvers
"
$ncells = 32
Cactus::cctk_show_schedule = yes
Cactus::presync_mode = "mixed-error"
Cactus::terminate = "time"
Cactus::cctk_final_time = 1.0
CarpetX::verbose = yes
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 = $ncells
CarpetX::ncells_y = $ncells
CarpetX::ncells_z = $ncells
CarpetX::periodic_x = no
CarpetX::periodic_y = no
CarpetX::periodic_z = no
CarpetX::periodic = yes
CarpetX::ghost_size = 1
CarpetX::max_num_levels = 1
DGCoordinates::dg_npoints = 4
ODESolvers::method = "RK4"
CarpetX::dtfac = 0.25
IO::out_dir = $parfile
IO::out_every = $ncells
CarpetX::out_silo_vars = ""   # "all"

File addition: standing-n8-64.par (-rw-r--r--)

[0.12054]

ActiveThorns = "
    CarpetX
    DGCoordinates
    DGWaveToy
    Formaline
    IOUtil
    ODESolvers
"
$ncells = 64
Cactus::cctk_show_schedule = yes
Cactus::presync_mode = "mixed-error"
Cactus::terminate = "time"
Cactus::cctk_final_time = 1.0
CarpetX::verbose = yes
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 = $ncells
CarpetX::ncells_y = $ncells
CarpetX::ncells_z = $ncells
CarpetX::periodic_x = no
CarpetX::periodic_y = no
CarpetX::periodic_z = no
CarpetX::periodic = yes
CarpetX::ghost_size = 1
CarpetX::max_num_levels = 1
DGCoordinates::dg_npoints = 8
ODESolvers::method = "DP87"
CarpetX::dtfac = 0.25
IO::out_dir = $parfile
IO::out_every = $ncells
CarpetX::out_silo_vars = ""   # "all"

File addition: standing-n8.par (-rw-r--r--)

[0.12054]

ActiveThorns = "
    CarpetX
    DGCoordinates
    DGWaveToy
    Formaline
    IOUtil
    ODESolvers
"
$ncells = 32
Cactus::cctk_show_schedule = yes
Cactus::presync_mode = "mixed-error"
Cactus::terminate = "time"
Cactus::cctk_final_time = 1.0
CarpetX::verbose = yes
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 = $ncells
CarpetX::ncells_y = $ncells
CarpetX::ncells_z = $ncells
CarpetX::periodic_x = no
CarpetX::periodic_y = no
CarpetX::periodic_z = no
CarpetX::periodic = yes
CarpetX::ghost_size = 1
CarpetX::max_num_levels = 1
DGCoordinates::dg_npoints = 8
ODESolvers::method = "DP87"
CarpetX::dtfac = 0.25
IO::out_dir = $parfile
IO::out_every = $ncells
CarpetX::out_silo_vars = ""   # "all"

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

[0.12054]

ActiveThorns = "
    CarpetX
    DGCoordinates
    DGWaveToy
    Formaline
    IOUtil
    ODESolvers
"
$ncells = 32
Cactus::cctk_show_schedule = yes
Cactus::presync_mode = "mixed-error"
Cactus::terminate = "time"
Cactus::cctk_final_time = 1.0
CarpetX::verbose = yes
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 = $ncells
CarpetX::ncells_y = $ncells
CarpetX::ncells_z = $ncells
CarpetX::periodic_x = no
CarpetX::periodic_y = no
CarpetX::periodic_z = no
CarpetX::periodic = yes
CarpetX::ghost_size = 1
CarpetX::max_num_levels = 1
DGCoordinates::dg_npoints = 2
ODESolvers::method = "RK2"
CarpetX::dtfac = 0.25
IO::out_dir = $parfile
IO::out_every = $ncells
CarpetX::out_silo_vars = ""   # "all"

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

[0.10723]

# Parameter definitions for thorn DGWaveToy
SHARES: DGCoordinates
USES CCTK_INT dg_npoints

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

[0.10723]

# Schedule definitions for thorn DGWaveToy
SCHEDULE DGWaveToy_Init AT initial
{
  LANG: C
  READS: DGCoordinates::coords(everywhere)
  WRITES: u(everywhere)
  WRITES: f(everywhere)
} "Initialise wave equation"
SCHEDULE DGWaveToy_RHS IN ODESolvers_RHS
{
  LANG: C
  READS: u(everywhere)
  READS: f(everywhere)
  WRITES: u_rhs(everywhere)
  WRITES: f_rhs(interior)
  SYNC: f_rhs
} "Calculate RHS of wave equation"
SCHEDULE DGWaveToy_RHSBoundaries IN ODESolvers_RHS AFTER DGWaveToy_RHS
{
  LANG: C
  WRITES: f_rhs(boundary)
} "Set RHS boundaries of wave equation"
SCHEDULE DGWaveToy_Energy AT analysis
{
  LANG: C
  READS: f(everywhere)
  READS: f_rhs(interior)
  WRITES: eps(interior)
  SYNC: eps
} "Calculate energy of wave equation"
SCHEDULE DGWaveToy_EnergyBoundaries AT analysis AFTER DGWaveToy_Energy
{
  LANG: C
  WRITES: eps(boundary)
} "Set energy boundary of wave equation"
SCHEDULE DGWaveToy_Error AT analysis
{
  LANG: C
  READS: DGCoordinates::coords(everywhere)
  READS: u(everywhere)
  READS: f(everywhere)
  WRITES: u_error(everywhere)
  WRITES: f_error(everywhere)
} "Calculate error of wave equation"

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

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

[0.19096]

#include <dg.hxx>
#include <loop.hxx>
#include <vect.hxx>
#include <cctk.h>
#include <cctk_Arguments_Checked.h>
#include <cctk_Parameters.h>
#include <cassert>
#include <cmath>
namespace DGWaveToy {
using namespace DG;
using namespace Arith;
using namespace std;
using Loop::dim;
extern "C" void DGWaveToy_Energy(CCTK_ARGUMENTS) {
  DECLARE_CCTK_ARGUMENTS_DGWaveToy_Energy;
  DECLARE_CCTK_PARAMETERS;
  const Loop::GF3D<const CCTK_REAL, 1, 1, 1> ft_rhs_(cctkGH, ft_rhs);
  const Loop::GF3D<const CCTK_REAL, 1, 1, 1> fx_(cctkGH, fx);
  const Loop::GF3D<const CCTK_REAL, 1, 1, 1> fy_(cctkGH, fy);
  const Loop::GF3D<const CCTK_REAL, 1, 1, 1> fz_(cctkGH, fz);
  const Loop::GF3D<CCTK_REAL, 1, 1, 1> eps_(cctkGH, eps);
  const Loop::GridDescBase grid(cctkGH);
  switch (dg_npoints) {
  case 2: {
    constexpr int N = 2;
    grid.loop_int<1, 1, 1>(grid.nghostzones, [&](const Loop::PointDesc &p) {
      const vect<CCTK_REAL, dim> dfx = get_derivs<N>(grid, fx_, p);
      const vect<CCTK_REAL, dim> dfy = get_derivs<N>(grid, fy_, p);
      const vect<CCTK_REAL, dim> dfz = get_derivs<N>(grid, fz_, p);
      eps_(p.I) = (pow(ft_rhs_(p.I), 2) + pow(dfx[0], 2) + pow(dfy[1], 2) +
                   pow(dfz[2], 2));
    });
    break;
  }
  case 4: {
    constexpr int N = 4;
    grid.loop_int<1, 1, 1>(grid.nghostzones, [&](const Loop::PointDesc &p) {
      const vect<CCTK_REAL, dim> dfx = get_derivs<N>(grid, fx_, p);
      const vect<CCTK_REAL, dim> dfy = get_derivs<N>(grid, fy_, p);
      const vect<CCTK_REAL, dim> dfz = get_derivs<N>(grid, fz_, p);
      eps_(p.I) = (pow(ft_rhs_(p.I), 2) + pow(dfx[0], 2) + pow(dfy[1], 2) +
                   pow(dfz[2], 2));
    });
    break;
  }
  case 8: {
    constexpr int N = 8;
    grid.loop_int<1, 1, 1>(grid.nghostzones, [&](const Loop::PointDesc &p) {
      const vect<CCTK_REAL, dim> dfx = get_derivs<N>(grid, fx_, p);
      const vect<CCTK_REAL, dim> dfy = get_derivs<N>(grid, fy_, p);
      const vect<CCTK_REAL, dim> dfz = get_derivs<N>(grid, fz_, p);
      eps_(p.I) = (pow(ft_rhs_(p.I), 2) + pow(dfx[0], 2) + pow(dfy[1], 2) +
                   pow(dfz[2], 2));
    });
    break;
  }
  default:
    assert(0);
  }
}
extern "C" void DGWaveToy_EnergyBoundaries(CCTK_ARGUMENTS) {
  DECLARE_CCTK_ARGUMENTS_DGWaveToy_EnergyBoundaries;
  DECLARE_CCTK_PARAMETERS;
  const Loop::GF3D<CCTK_REAL, 1, 1, 1> eps_(cctkGH, eps);
  const Loop::GridDescBase grid(cctkGH);
  grid.loop_bnd<1, 1, 1>(grid.nghostzones,
                         [&](const Loop::PointDesc &p) { eps_(p.I) = 0; });
}
} // namespace DGWaveToy

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

[0.19096]

#include <dg.hxx>
#include <loop.hxx>
#include <vect.hxx>
#include <cctk.h>
#include <cctk_Arguments_Checked.h>
#include <cctk_Parameters.h>
#include <cmath>
namespace DGWaveToy {
using namespace DG;
using namespace std;
using Loop::dim;
extern "C" void DGWaveToy_Error(CCTK_ARGUMENTS) {
  DECLARE_CCTK_ARGUMENTS_DGWaveToy_Error;
  DECLARE_CCTK_PARAMETERS;
  const Loop::GF3D<const CCTK_REAL, 1, 1, 1> coordx_(cctkGH, coordx);
  const Loop::GF3D<const CCTK_REAL, 1, 1, 1> coordy_(cctkGH, coordy);
  const Loop::GF3D<const CCTK_REAL, 1, 1, 1> coordz_(cctkGH, coordz);
  const Loop::GF3D<const CCTK_REAL, 1, 1, 1> u_(cctkGH, u);
  const Loop::GF3D<const CCTK_REAL, 1, 1, 1> ft_(cctkGH, ft);
  const Loop::GF3D<const CCTK_REAL, 1, 1, 1> fx_(cctkGH, fx);
  const Loop::GF3D<const CCTK_REAL, 1, 1, 1> fy_(cctkGH, fy);
  const Loop::GF3D<const CCTK_REAL, 1, 1, 1> fz_(cctkGH, fz);
  const Loop::GF3D<CCTK_REAL, 1, 1, 1> u_error_(cctkGH, u_error);
  const Loop::GF3D<CCTK_REAL, 1, 1, 1> ft_error_(cctkGH, ft_error);
  const Loop::GF3D<CCTK_REAL, 1, 1, 1> fx_error_(cctkGH, fx_error);
  const Loop::GF3D<CCTK_REAL, 1, 1, 1> fy_error_(cctkGH, fy_error);
  const Loop::GF3D<CCTK_REAL, 1, 1, 1> fz_error_(cctkGH, fz_error);
  const CCTK_REAL kx = 2 * M_PI;
  const CCTK_REAL ky = 2 * M_PI;
  const CCTK_REAL kz = 2 * M_PI;
  const CCTK_REAL omega = sqrt(pow(kx, 2) + pow(ky, 2) + pow(kz, 2));
  Loop::loop_all<1, 1, 1>(cctkGH, [&](const Loop::PointDesc &p) {
    const CCTK_REAL t = cctk_time;
    const CCTK_REAL x = coordx_(p.I);
    const CCTK_REAL y = coordy_(p.I);
    const CCTK_REAL z = coordz_(p.I);
    u_error_(p.I) =
        u_(p.I) - cos(omega * t) * cos(kx * x) * cos(ky * y) * cos(kz * z);
    ft_error_(p.I) = ft_(p.I) - -omega * sin(omega * t) * cos(kx * x) *
                                    cos(ky * y) * cos(kz * z);
    fx_error_(p.I) = fx_(p.I) - -kx * cos(omega * t) * sin(kx * x) *
                                    cos(ky * y) * cos(kz * z);
    fy_error_(p.I) = fy_(p.I) - -ky * cos(omega * t) * cos(kx * x) *
                                    sin(ky * y) * cos(kz * z);
    fz_error_(p.I) = fz_(p.I) - -kz * cos(omega * t) * cos(kx * x) *
                                    cos(ky * y) * sin(kz * z);
  });
}
} // namespace DGWaveToy

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

[0.19096]

#include <loop.hxx>
#include <cctk.h>
#include <cctk_Arguments_Checked.h>
#include <cctk_Parameters.h>
#include <cmath>
namespace DGWaveToy {
using namespace std;
extern "C" void DGWaveToy_Init(CCTK_ARGUMENTS) {
  DECLARE_CCTK_ARGUMENTS_DGWaveToy_Init;
  DECLARE_CCTK_PARAMETERS;
  const Loop::GF3D<const CCTK_REAL, 1, 1, 1> coordx_(cctkGH, coordx);
  const Loop::GF3D<const CCTK_REAL, 1, 1, 1> coordy_(cctkGH, coordy);
  const Loop::GF3D<const CCTK_REAL, 1, 1, 1> coordz_(cctkGH, coordz);
  const Loop::GF3D<CCTK_REAL, 1, 1, 1> u_(cctkGH, u);
  const Loop::GF3D<CCTK_REAL, 1, 1, 1> ft_(cctkGH, ft);
  const Loop::GF3D<CCTK_REAL, 1, 1, 1> fx_(cctkGH, fx);
  const Loop::GF3D<CCTK_REAL, 1, 1, 1> fy_(cctkGH, fy);
  const Loop::GF3D<CCTK_REAL, 1, 1, 1> fz_(cctkGH, fz);
 
  const CCTK_REAL kx = 2 * M_PI;
  const CCTK_REAL ky = 2 * M_PI;
  const CCTK_REAL kz = 2 * M_PI;
  const CCTK_REAL omega = sqrt(pow(kx, 2) + pow(ky, 2) + pow(kz, 2));
  Loop::loop_all<1, 1, 1>(cctkGH, [&](const Loop::PointDesc &p) {
    const CCTK_REAL t = cctk_time;
    const CCTK_REAL x = coordx_(p.I);
    const CCTK_REAL y = coordy_(p.I);
    const CCTK_REAL z = coordz_(p.I);
    u_(p.I) = cos(omega * t) * cos(kx * x) * cos(ky * y) * cos(kz * z);
    ft_(p.I) =
        -omega * sin(omega * t) * cos(kx * x) * cos(ky * y) * cos(kz * z);
    fx_(p.I) = -kx * cos(omega * t) * sin(kx * x) * cos(ky * y) * cos(kz * z);
    fy_(p.I) = -ky * cos(omega * t) * cos(kx * x) * sin(ky * y) * cos(kz * z);
    fz_(p.I) = -kz * cos(omega * t) * cos(kx * x) * cos(ky * y) * sin(kz * z);
  });
}
} // namespace DGWaveToy

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

[0.19096]

# Main make.code.defn file for thorn DGCoordinates
# Source files in this directory
SRCS = energy.cxx error.cxx init.cxx rhs.cxx
# Subdirectories containing source files
SUBDIRS =

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

[0.19096]

#include <dg.hxx>
#include <loop.hxx>
#include <vect.hxx>
#include <cctk.h>
#include <cctk_Arguments_Checked.h>
#include <cctk_Parameters.h>
#include <cassert>
#include <cmath>
namespace DGWaveToy {
using namespace DG;
using namespace std;
using Loop::dim;
extern "C" void DGWaveToy_RHS(CCTK_ARGUMENTS) {
  DECLARE_CCTK_ARGUMENTS_DGWaveToy_RHS;
  DECLARE_CCTK_PARAMETERS;
  const Loop::GF3D<const CCTK_REAL, 1, 1, 1> u_(cctkGH, u);
  const Loop::GF3D<const CCTK_REAL, 1, 1, 1> ft_(cctkGH, ft);
  const Loop::GF3D<const CCTK_REAL, 1, 1, 1> fx_(cctkGH, fx);
  const Loop::GF3D<const CCTK_REAL, 1, 1, 1> fy_(cctkGH, fy);
  const Loop::GF3D<const CCTK_REAL, 1, 1, 1> fz_(cctkGH, fz);
  const Loop::GF3D<CCTK_REAL, 1, 1, 1> u_rhs_(cctkGH, u_rhs);
  const Loop::GF3D<CCTK_REAL, 1, 1, 1> ft_rhs_(cctkGH, ft_rhs);
  const Loop::GF3D<CCTK_REAL, 1, 1, 1> fx_rhs_(cctkGH, fx_rhs);
  const Loop::GF3D<CCTK_REAL, 1, 1, 1> fy_rhs_(cctkGH, fy_rhs);
  const Loop::GF3D<CCTK_REAL, 1, 1, 1> fz_rhs_(cctkGH, fz_rhs);
  const Loop::GridDescBase grid(cctkGH);
  grid.loop_all<1, 1, 1>(grid.nghostzones, [&](const Loop::PointDesc &p) {
    u_rhs_(p.I) = ft_(p.I);
  });
  switch (dg_npoints) {
  case 2: {
    constexpr int N = 2;
    grid.loop_int<1, 1, 1>(grid.nghostzones, [&](const Loop::PointDesc &p) {
      const vect<CCTK_REAL, dim> dft = get_derivs<N>(grid, ft_, p);
      const vect<CCTK_REAL, dim> dfx = get_derivs<N>(grid, fx_, p);
      const vect<CCTK_REAL, dim> dfy = get_derivs<N>(grid, fy_, p);
      const vect<CCTK_REAL, dim> dfz = get_derivs<N>(grid, fz_, p);
      ft_rhs_(p.I) = dfx[0] + dfy[1] + dfz[2];
      fx_rhs_(p.I) = dft[0];
      fy_rhs_(p.I) = dft[1];
      fz_rhs_(p.I) = dft[2];
    });
    break;
  }
  case 4: {
    constexpr int N = 4;
    grid.loop_int<1, 1, 1>(grid.nghostzones, [&](const Loop::PointDesc &p) {
      const vect<CCTK_REAL, dim> dft = get_derivs<N>(grid, ft_, p);
      const vect<CCTK_REAL, dim> dfx = get_derivs<N>(grid, fx_, p);
      const vect<CCTK_REAL, dim> dfy = get_derivs<N>(grid, fy_, p);
      const vect<CCTK_REAL, dim> dfz = get_derivs<N>(grid, fz_, p);
      ft_rhs_(p.I) = dfx[0] + dfy[1] + dfz[2];
      fx_rhs_(p.I) = dft[0];
      fy_rhs_(p.I) = dft[1];
      fz_rhs_(p.I) = dft[2];
    });
    break;
  }
  case 8: {
    constexpr int N = 8;
    grid.loop_int<1, 1, 1>(grid.nghostzones, [&](const Loop::PointDesc &p) {
      const vect<CCTK_REAL, dim> dft = get_derivs<N>(grid, ft_, p);
      const vect<CCTK_REAL, dim> dfx = get_derivs<N>(grid, fx_, p);
      const vect<CCTK_REAL, dim> dfy = get_derivs<N>(grid, fy_, p);
      const vect<CCTK_REAL, dim> dfz = get_derivs<N>(grid, fz_, p);
      ft_rhs_(p.I) = dfx[0] + dfy[1] + dfz[2];
      fx_rhs_(p.I) = dft[0];
      fy_rhs_(p.I) = dft[1];
      fz_rhs_(p.I) = dft[2];
    });
    break;
  }
  default:
    assert(0);
  }
}
extern "C" void DGWaveToy_RHSBoundaries(CCTK_ARGUMENTS) {
  DECLARE_CCTK_ARGUMENTS_DGWaveToy_RHSBoundaries;
  DECLARE_CCTK_PARAMETERS;
  const Loop::GF3D<CCTK_REAL, 1, 1, 1> ft_rhs_(cctkGH, ft_rhs);
  const Loop::GF3D<CCTK_REAL, 1, 1, 1> fx_rhs_(cctkGH, fx_rhs);
  const Loop::GF3D<CCTK_REAL, 1, 1, 1> fy_rhs_(cctkGH, fy_rhs);
  const Loop::GF3D<CCTK_REAL, 1, 1, 1> fz_rhs_(cctkGH, fz_rhs);
  const Loop::GridDescBase grid(cctkGH);
  grid.loop_bnd<1, 1, 1>(grid.nghostzones, [&](const Loop::PointDesc &p) {
    ft_rhs_(p.I) = 0;
    fx_rhs_(p.I) = 0;
    fy_rhs_(p.I) = 0;
    fz_rhs_(p.I) = 0;
  });
}
} // namespace DGWaveToy

Insertion in azure-pipelines/carpetx.th at line 142 [4.11]
[3.38]
[2.2354]
```
CarpetX/DGCoordinates
CarpetX/DGWaveToy
```

DGWaveToy: Implement WaveToy with DG derivatives

Dependencies

In channels

Change contents

File addition: DGCoordinates (drwxr-xr-x)

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

File addition: interface.ccl (-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: coordinates.cxx (-rw-r--r--)

File addition: dg.hxx (-rw-r--r--)

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

File addition: DGWaveToy (drwxr-xr-x)

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

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

File addition: par (drwxr-xr-x)

File addition: standing-128.par (-rw-r--r--)

File addition: standing-64.par (-rw-r--r--)

File addition: standing-n4-64.par (-rw-r--r--)

File addition: standing-n4.par (-rw-r--r--)

File addition: standing-n8-64.par (-rw-r--r--)

File addition: standing-n8.par (-rw-r--r--)

File addition: standing.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: energy.cxx (-rw-r--r--)

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

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

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

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

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