Cactus Code Thorn TestODESolvers
Author(s)    : Erik Schnetter <schnetter@gmail.com>
Maintainer(s): Erik Schnetter <schnetter@gmail.com>
Licence      : LGPL
--------------------------------------------------------------------------
1. Purpose
Test the ODESolvers thorn

# Configuration definitions for thorn TestODESolvers

# Interface definition for thorn TestODESolvers
IMPLEMENTS: TestODESolvers
CCTK_REAL state TYPE=gf TAGS='rhs="TestODESolvers::rhs"' "State vector"
CCTK_REAL rhs TYPE=gf "RHS of state vector"
CCTK_REAL error TYPE=gf "Error in state vector"

ActiveThorns = "
    CarpetX
    IOUtil
    ODESolvers
    TestODESolvers
"
CarpetX::ncells_x = 1
CarpetX::ncells_y = 1
CarpetX::ncells_z = 1
CarpetX::blocking_factor_x = 1
CarpetX::blocking_factor_y = 1
CarpetX::blocking_factor_z = 1
CarpetX::ghost_size = 0
CarpetX::dtfac = 0.1
ODESolvers::method = "Euler"
TestODESolvers::order = 1
CarpetX::out_plotfile_groups = ""
CarpetX::out_tsv = yes
IO::out_dir = $parfile
IO::out_every = 1

# Parameter definitions for thorn TestODESolvers
CCTK_INT order "order p, where u(t) = (1+t)^p" STEERABLE=always
{
  0:* :: ""
} 1

# Schedule definitions for thorn TestODESolvers
STORAGE: state rhs
SCHEDULE TestODESolvers_initial AT initial
{
  LANG: C
  WRITES: state
} "Initialize state vector"
SCHEDULE TestODESolvers_RHS IN ODESolvers_RHS
{
  LANG: C
  READS: state
  WRITES: rhs
} "RHS function"
SCHEDULE TestODESolvers_boundary AT poststep
{
  LANG: C
  SYNC: state
} "Apply boundary conditions to state vector"
SCHEDULE TestODESolvers_error AT analysis
{
  LANG: C
  READS: state
  WRITES: error
} "Calculate error in state vector"

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

#include <cctk.h>
#include <cctk_Parameters.h>
#include <cctk_Arguments_Checked.h>
#include <cmath>
namespace TestODESolvers {
using namespace std;
extern "C" void TestODESolvers_initial(CCTK_ARGUMENTS) {
  DECLARE_CCTK_ARGUMENTS_TestODESolvers_initial;
  DECLARE_CCTK_PARAMETERS;
  const CCTK_REAL u0 = pow(1 + cctk_time, order);
  for (int k = 0; k < cctk_lsh[2]; ++k) {
    for (int j = 0; j < cctk_lsh[1]; ++j) {
#pragma omp simd
      for (int i = 0; i < cctk_lsh[0]; ++i) {
        int ind = CCTK_GFINDEX3D(cctkGH, i, j, k);
        state[ind] = u0;
      }
    }
  }
}
extern "C" void TestODESolvers_boundary(CCTK_ARGUMENTS) {
  DECLARE_CCTK_ARGUMENTS_TestODESolvers_boundary;
  DECLARE_CCTK_PARAMETERS;
  // do nothing
}
extern "C" void TestODESolvers_RHS(CCTK_ARGUMENTS) {
  DECLARE_CCTK_ARGUMENTS_TestODESolvers_RHS;
  DECLARE_CCTK_PARAMETERS;
  // u(t) = (1+t)^p
  // d/dt u = p (1+t)^(p-1)
  for (int k = 0; k < cctk_lsh[2]; ++k) {
    for (int j = 0; j < cctk_lsh[1]; ++j) {
#pragma omp simd
      for (int i = 0; i < cctk_lsh[0]; ++i) {
        int ind = CCTK_GFINDEX3D(cctkGH, i, j, k);
        rhs[ind] = order * pow(1 + cctk_time, order - 1);
      }
    }
  }
}
extern "C" void TestODESolvers_error(CCTK_ARGUMENTS) {
  DECLARE_CCTK_ARGUMENTS_TestODESolvers_error;
  DECLARE_CCTK_PARAMETERS;
  const CCTK_REAL u0 = pow(1 + cctk_time, order);
  for (int k = 0; k < cctk_lsh[2]; ++k) {
    for (int j = 0; j < cctk_lsh[1]; ++j) {
#pragma omp simd
      for (int i = 0; i < cctk_lsh[0]; ++i) {
        int ind = CCTK_GFINDEX3D(cctkGH, i, j, k);
        error[ind] = state[ind] - u0;
      }
    }
  }
}
} // namespace TestODESolvers

ActiveThorns = "
    CarpetX
    IOUtil
    ODESolvers
    TestODESolvers
"
CarpetX::ncells_x = 1
CarpetX::ncells_y = 1
CarpetX::ncells_z = 1
CarpetX::blocking_factor_x = 1
CarpetX::blocking_factor_y = 1
CarpetX::blocking_factor_z = 1
CarpetX::ghost_size = 0
CarpetX::dtfac = 0.1
ODESolvers::method = "Euler"
TestODESolvers::order = 1
CarpetX::out_plotfile_groups = ""
CarpetX::out_tsv = yes
IO::out_dir = $parfile
IO::out_every = 1

ActiveThorns = "
    CarpetX
    IOUtil
    ODESolvers
    TestODESolvers
"
CarpetX::ncells_x = 1
CarpetX::ncells_y = 1
CarpetX::ncells_z = 1
CarpetX::blocking_factor_x = 1
CarpetX::blocking_factor_y = 1
CarpetX::blocking_factor_z = 1
CarpetX::ghost_size = 0
CarpetX::dtfac = 0.1
ODESolvers::method = "RK2"
TestODESolvers::order = 2
CarpetX::out_plotfile_groups = ""
CarpetX::out_tsv = yes
IO::out_dir = $parfile
IO::out_every = 1

CarpetX/TestODESolvers