From 5dbceebfbc51ed8a947211bd448d546d7aad3b12 Mon Sep 17 00:00:00 2001 From: Erik Schnetter Date: Sun, 3 Apr 2011 15:52:06 -0400 Subject: Update auto-generated code as generated by the current version of Kranc --- ML_WaveToy/src/Boundaries.cc | 249 +++++++++++++++++++++++++++++++++++ ML_WaveToy/src/RegisterMoL.cc | 18 +++ ML_WaveToy/src/RegisterSymmetries.cc | 29 ++++ ML_WaveToy/src/Startup.cc | 10 ++ ML_WaveToy/src/WT_Gaussian.cc | 109 +++++++++++++++ ML_WaveToy/src/WT_RHS.cc | 130 ++++++++++++++++++ 6 files changed, 545 insertions(+) create mode 100644 ML_WaveToy/src/Boundaries.cc create mode 100644 ML_WaveToy/src/RegisterMoL.cc create mode 100644 ML_WaveToy/src/RegisterSymmetries.cc create mode 100644 ML_WaveToy/src/Startup.cc create mode 100644 ML_WaveToy/src/WT_Gaussian.cc create mode 100644 ML_WaveToy/src/WT_RHS.cc (limited to 'ML_WaveToy') diff --git a/ML_WaveToy/src/Boundaries.cc b/ML_WaveToy/src/Boundaries.cc new file mode 100644 index 0000000..1e5f99e --- /dev/null +++ b/ML_WaveToy/src/Boundaries.cc @@ -0,0 +1,249 @@ +/* File produced by Kranc */ + +#include "cctk.h" +#include "cctk_Arguments.h" +#include "cctk_Parameters.h" +#include "cctk_Faces.h" +#include "util_Table.h" +#include "Symmetry.h" + + +/* the boundary treatment is split into 3 steps: */ +/* 1. excision */ +/* 2. symmetries */ +/* 3. "other" boundary conditions, e.g. radiative */ + +/* to simplify scheduling and testing, the 3 steps */ +/* are currently applied in separate functions */ + + +extern "C" void ML_WaveToy_CheckBoundaries(CCTK_ARGUMENTS) +{ + DECLARE_CCTK_ARGUMENTS; + DECLARE_CCTK_PARAMETERS; + + return; +} + +extern "C" void ML_WaveToy_SelectBoundConds(CCTK_ARGUMENTS) +{ + DECLARE_CCTK_ARGUMENTS; + DECLARE_CCTK_PARAMETERS; + + CCTK_INT ierr = 0; + + if (CCTK_EQUALS(WT_rho_bound, "none" ) || + CCTK_EQUALS(WT_rho_bound, "static") || + CCTK_EQUALS(WT_rho_bound, "flat" ) || + CCTK_EQUALS(WT_rho_bound, "zero" ) ) + { + ierr = Boundary_SelectGroupForBC(cctkGH, CCTK_ALL_FACES, 1, -1, + "ML_WaveToy::WT_rho", WT_rho_bound); + if (ierr < 0) + CCTK_WARN(0, "Failed to register WT_rho_bound BC for ML_WaveToy::WT_rho!"); + } + + if (CCTK_EQUALS(WT_u_bound, "none" ) || + CCTK_EQUALS(WT_u_bound, "static") || + CCTK_EQUALS(WT_u_bound, "flat" ) || + CCTK_EQUALS(WT_u_bound, "zero" ) ) + { + ierr = Boundary_SelectGroupForBC(cctkGH, CCTK_ALL_FACES, 1, -1, + "ML_WaveToy::WT_u", WT_u_bound); + if (ierr < 0) + CCTK_WARN(0, "Failed to register WT_u_bound BC for ML_WaveToy::WT_u!"); + } + + if (CCTK_EQUALS(rho_bound, "none" ) || + CCTK_EQUALS(rho_bound, "static") || + CCTK_EQUALS(rho_bound, "flat" ) || + CCTK_EQUALS(rho_bound, "zero" ) ) + { + ierr = Boundary_SelectVarForBC(cctkGH, CCTK_ALL_FACES, 1, -1, + "ML_WaveToy::rho", rho_bound); + if (ierr < 0) + CCTK_WARN(0, "Failed to register rho_bound BC for ML_WaveToy::rho!"); + } + + if (CCTK_EQUALS(u_bound, "none" ) || + CCTK_EQUALS(u_bound, "static") || + CCTK_EQUALS(u_bound, "flat" ) || + CCTK_EQUALS(u_bound, "zero" ) ) + { + ierr = Boundary_SelectVarForBC(cctkGH, CCTK_ALL_FACES, 1, -1, + "ML_WaveToy::u", u_bound); + if (ierr < 0) + CCTK_WARN(0, "Failed to register u_bound BC for ML_WaveToy::u!"); + } + + if (CCTK_EQUALS(WT_rho_bound, "radiative")) + { + /* select radiation boundary condition */ + static CCTK_INT handle_WT_rho_bound = -1; + if (handle_WT_rho_bound < 0) handle_WT_rho_bound = Util_TableCreate(UTIL_TABLE_FLAGS_CASE_INSENSITIVE); + if (handle_WT_rho_bound < 0) CCTK_WARN(0, "could not create table!"); + if (Util_TableSetReal(handle_WT_rho_bound , WT_rho_bound_limit, "LIMIT") < 0) + CCTK_WARN(0, "could not set LIMIT value in table!"); + if (Util_TableSetReal(handle_WT_rho_bound ,WT_rho_bound_speed, "SPEED") < 0) + CCTK_WARN(0, "could not set SPEED value in table!"); + + ierr = Boundary_SelectGroupForBC(cctkGH, CCTK_ALL_FACES, 1, handle_WT_rho_bound, + "ML_WaveToy::WT_rho", "Radiation"); + + if (ierr < 0) + CCTK_WARN(0, "Failed to register Radiation BC for ML_WaveToy::WT_rho!"); + + } + + if (CCTK_EQUALS(WT_u_bound, "radiative")) + { + /* select radiation boundary condition */ + static CCTK_INT handle_WT_u_bound = -1; + if (handle_WT_u_bound < 0) handle_WT_u_bound = Util_TableCreate(UTIL_TABLE_FLAGS_CASE_INSENSITIVE); + if (handle_WT_u_bound < 0) CCTK_WARN(0, "could not create table!"); + if (Util_TableSetReal(handle_WT_u_bound , WT_u_bound_limit, "LIMIT") < 0) + CCTK_WARN(0, "could not set LIMIT value in table!"); + if (Util_TableSetReal(handle_WT_u_bound ,WT_u_bound_speed, "SPEED") < 0) + CCTK_WARN(0, "could not set SPEED value in table!"); + + ierr = Boundary_SelectGroupForBC(cctkGH, CCTK_ALL_FACES, 1, handle_WT_u_bound, + "ML_WaveToy::WT_u", "Radiation"); + + if (ierr < 0) + CCTK_WARN(0, "Failed to register Radiation BC for ML_WaveToy::WT_u!"); + + } + + if (CCTK_EQUALS(rho_bound, "radiative")) + { + /* select radiation boundary condition */ + static CCTK_INT handle_rho_bound = -1; + if (handle_rho_bound < 0) handle_rho_bound = Util_TableCreate(UTIL_TABLE_FLAGS_CASE_INSENSITIVE); + if (handle_rho_bound < 0) CCTK_WARN(0, "could not create table!"); + if (Util_TableSetReal(handle_rho_bound , rho_bound_limit, "LIMIT") < 0) + CCTK_WARN(0, "could not set LIMIT value in table!"); + if (Util_TableSetReal(handle_rho_bound ,rho_bound_speed, "SPEED") < 0) + CCTK_WARN(0, "could not set SPEED value in table!"); + + ierr = Boundary_SelectVarForBC(cctkGH, CCTK_ALL_FACES, 1, handle_rho_bound, + "ML_WaveToy::rho", "Radiation"); + + if (ierr < 0) + CCTK_WARN(0, "Failed to register Radiation BC for ML_WaveToy::rho!"); + + } + + if (CCTK_EQUALS(u_bound, "radiative")) + { + /* select radiation boundary condition */ + static CCTK_INT handle_u_bound = -1; + if (handle_u_bound < 0) handle_u_bound = Util_TableCreate(UTIL_TABLE_FLAGS_CASE_INSENSITIVE); + if (handle_u_bound < 0) CCTK_WARN(0, "could not create table!"); + if (Util_TableSetReal(handle_u_bound , u_bound_limit, "LIMIT") < 0) + CCTK_WARN(0, "could not set LIMIT value in table!"); + if (Util_TableSetReal(handle_u_bound ,u_bound_speed, "SPEED") < 0) + CCTK_WARN(0, "could not set SPEED value in table!"); + + ierr = Boundary_SelectVarForBC(cctkGH, CCTK_ALL_FACES, 1, handle_u_bound, + "ML_WaveToy::u", "Radiation"); + + if (ierr < 0) + CCTK_WARN(0, "Failed to register Radiation BC for ML_WaveToy::u!"); + + } + + if (CCTK_EQUALS(WT_rho_bound, "scalar")) + { + /* select scalar boundary condition */ + static CCTK_INT handle_WT_rho_bound = -1; + if (handle_WT_rho_bound < 0) handle_WT_rho_bound = Util_TableCreate(UTIL_TABLE_FLAGS_CASE_INSENSITIVE); + if (handle_WT_rho_bound < 0) CCTK_WARN(0, "could not create table!"); + if (Util_TableSetReal(handle_WT_rho_bound ,WT_rho_bound_scalar, "SCALAR") < 0) + CCTK_WARN(0, "could not set SCALAR value in table!"); + + ierr = Boundary_SelectGroupForBC(cctkGH, CCTK_ALL_FACES, 1, handle_WT_rho_bound, + "ML_WaveToy::WT_rho", "scalar"); + + if (ierr < 0) + CCTK_WARN(0, "Failed to register Scalar BC for ML_WaveToy::WT_rho!"); + + } + + if (CCTK_EQUALS(WT_u_bound, "scalar")) + { + /* select scalar boundary condition */ + static CCTK_INT handle_WT_u_bound = -1; + if (handle_WT_u_bound < 0) handle_WT_u_bound = Util_TableCreate(UTIL_TABLE_FLAGS_CASE_INSENSITIVE); + if (handle_WT_u_bound < 0) CCTK_WARN(0, "could not create table!"); + if (Util_TableSetReal(handle_WT_u_bound ,WT_u_bound_scalar, "SCALAR") < 0) + CCTK_WARN(0, "could not set SCALAR value in table!"); + + ierr = Boundary_SelectGroupForBC(cctkGH, CCTK_ALL_FACES, 1, handle_WT_u_bound, + "ML_WaveToy::WT_u", "scalar"); + + if (ierr < 0) + CCTK_WARN(0, "Failed to register Scalar BC for ML_WaveToy::WT_u!"); + + } + + if (CCTK_EQUALS(rho_bound, "scalar")) + { + /* select scalar boundary condition */ + static CCTK_INT handle_rho_bound = -1; + if (handle_rho_bound < 0) handle_rho_bound = Util_TableCreate(UTIL_TABLE_FLAGS_CASE_INSENSITIVE); + if (handle_rho_bound < 0) CCTK_WARN(0, "could not create table!"); + if (Util_TableSetReal(handle_rho_bound ,rho_bound_scalar, "SCALAR") < 0) + CCTK_WARN(0, "could not set SCALAR value in table!"); + + ierr = Boundary_SelectVarForBC(cctkGH, CCTK_ALL_FACES, 1, handle_rho_bound, + "ML_WaveToy::rho", "scalar"); + + if (ierr < 0) + CCTK_WARN(0, "Error in registering Scalar BC for ML_WaveToy::rho!"); + + } + + if (CCTK_EQUALS(u_bound, "scalar")) + { + /* select scalar boundary condition */ + static CCTK_INT handle_u_bound = -1; + if (handle_u_bound < 0) handle_u_bound = Util_TableCreate(UTIL_TABLE_FLAGS_CASE_INSENSITIVE); + if (handle_u_bound < 0) CCTK_WARN(0, "could not create table!"); + if (Util_TableSetReal(handle_u_bound ,u_bound_scalar, "SCALAR") < 0) + CCTK_WARN(0, "could not set SCALAR value in table!"); + + ierr = Boundary_SelectVarForBC(cctkGH, CCTK_ALL_FACES, 1, handle_u_bound, + "ML_WaveToy::u", "scalar"); + + if (ierr < 0) + CCTK_WARN(0, "Error in registering Scalar BC for ML_WaveToy::u!"); + + } + return; +} + + + +/* template for entries in parameter file: +#$bound$#ML_WaveToy::WT_rho_bound = "skip" +#$bound$#ML_WaveToy::WT_rho_bound_speed = 1.0 +#$bound$#ML_WaveToy::WT_rho_bound_limit = 0.0 +#$bound$#ML_WaveToy::WT_rho_bound_scalar = 0.0 + +#$bound$#ML_WaveToy::WT_u_bound = "skip" +#$bound$#ML_WaveToy::WT_u_bound_speed = 1.0 +#$bound$#ML_WaveToy::WT_u_bound_limit = 0.0 +#$bound$#ML_WaveToy::WT_u_bound_scalar = 0.0 + +#$bound$#ML_WaveToy::rho_bound = "skip" +#$bound$#ML_WaveToy::rho_bound_speed = 1.0 +#$bound$#ML_WaveToy::rho_bound_limit = 0.0 +#$bound$#ML_WaveToy::rho_bound_scalar = 0.0 + +#$bound$#ML_WaveToy::u_bound = "skip" +#$bound$#ML_WaveToy::u_bound_speed = 1.0 +#$bound$#ML_WaveToy::u_bound_limit = 0.0 +#$bound$#ML_WaveToy::u_bound_scalar = 0.0 + +*/ + diff --git a/ML_WaveToy/src/RegisterMoL.cc b/ML_WaveToy/src/RegisterMoL.cc new file mode 100644 index 0000000..21434f5 --- /dev/null +++ b/ML_WaveToy/src/RegisterMoL.cc @@ -0,0 +1,18 @@ +/* File produced by Kranc */ + +#include "cctk.h" +#include "cctk_Arguments.h" +#include "cctk_Parameters.h" + +extern "C" void ML_WaveToy_RegisterVars(CCTK_ARGUMENTS) +{ + DECLARE_CCTK_ARGUMENTS; + DECLARE_CCTK_PARAMETERS; + + CCTK_INT ierr = 0; + + /* Register all the evolved grid functions with MoL */ + ierr += MoLRegisterEvolved(CCTK_VarIndex("ML_WaveToy::rho"), CCTK_VarIndex("ML_WaveToy::rhorhs")); + ierr += MoLRegisterEvolved(CCTK_VarIndex("ML_WaveToy::u"), CCTK_VarIndex("ML_WaveToy::urhs")); + return; +} diff --git a/ML_WaveToy/src/RegisterSymmetries.cc b/ML_WaveToy/src/RegisterSymmetries.cc new file mode 100644 index 0000000..ae89b96 --- /dev/null +++ b/ML_WaveToy/src/RegisterSymmetries.cc @@ -0,0 +1,29 @@ +/* File produced by Kranc */ + +#include "cctk.h" +#include "cctk_Arguments.h" +#include "cctk_Parameters.h" +#include "Symmetry.h" + +extern "C" void ML_WaveToy_RegisterSymmetries(CCTK_ARGUMENTS) +{ + DECLARE_CCTK_ARGUMENTS; + DECLARE_CCTK_PARAMETERS; + + + /* array holding symmetry definitions */ + CCTK_INT sym[3]; + + + /* Register symmetries of grid functions */ + sym[0] = 1; + sym[1] = 1; + sym[2] = 1; + SetCartSymVN(cctkGH, sym, "ML_WaveToy::rho"); + + sym[0] = 1; + sym[1] = 1; + sym[2] = 1; + SetCartSymVN(cctkGH, sym, "ML_WaveToy::u"); + +} diff --git a/ML_WaveToy/src/Startup.cc b/ML_WaveToy/src/Startup.cc new file mode 100644 index 0000000..5f3d2bd --- /dev/null +++ b/ML_WaveToy/src/Startup.cc @@ -0,0 +1,10 @@ +/* File produced by Kranc */ + +#include "cctk.h" + +extern "C" int ML_WaveToy_Startup(void) +{ + const char * banner = "ML_WaveToy"; + CCTK_RegisterBanner(banner); + return 0; +} diff --git a/ML_WaveToy/src/WT_Gaussian.cc b/ML_WaveToy/src/WT_Gaussian.cc new file mode 100644 index 0000000..f4d0da6 --- /dev/null +++ b/ML_WaveToy/src/WT_Gaussian.cc @@ -0,0 +1,109 @@ +/* File produced by Kranc */ + +#define KRANC_C + +#include +#include +#include +#include +#include "cctk.h" +#include "cctk_Arguments.h" +#include "cctk_Parameters.h" +#include "GenericFD.h" +#include "Differencing.h" +#include "loopcontrol.h" + +/* Define macros used in calculations */ +#define INITVALUE (42) +#define QAD(x) (SQR(SQR(x))) +#define INV(x) ((1.0) / (x)) +#define SQR(x) ((x) * (x)) +#define CUB(x) ((x) * (x) * (x)) + +static void WT_Gaussian_Body(cGH const * restrict const cctkGH, int const dir, int const face, CCTK_REAL const normal[3], CCTK_REAL const tangentA[3], CCTK_REAL const tangentB[3], int const min[3], int const max[3], int const n_subblock_gfs, CCTK_REAL * restrict const subblock_gfs[]) +{ + DECLARE_CCTK_ARGUMENTS; + DECLARE_CCTK_PARAMETERS; + + + /* Declare finite differencing variables */ + + if (verbose > 1) + { + CCTK_VInfo(CCTK_THORNSTRING,"Entering WT_Gaussian_Body"); + } + + if (cctk_iteration % WT_Gaussian_calc_every != WT_Gaussian_calc_offset) + { + return; + } + + const char *groups[] = {"ML_WaveToy::WT_rho","ML_WaveToy::WT_u"}; + GenericFD_AssertGroupStorage(cctkGH, "WT_Gaussian", 2, groups); + + /* Include user-supplied include files */ + + /* Initialise finite differencing variables */ + ptrdiff_t const di = 1; + ptrdiff_t const dj = CCTK_GFINDEX3D(cctkGH,0,1,0) - CCTK_GFINDEX3D(cctkGH,0,0,0); + ptrdiff_t const dk = CCTK_GFINDEX3D(cctkGH,0,0,1) - CCTK_GFINDEX3D(cctkGH,0,0,0); + CCTK_REAL const dx = ToReal(CCTK_DELTA_SPACE(0)); + CCTK_REAL const dy = ToReal(CCTK_DELTA_SPACE(1)); + CCTK_REAL const dz = ToReal(CCTK_DELTA_SPACE(2)); + CCTK_REAL const dxi = INV(dx); + CCTK_REAL const dyi = INV(dy); + CCTK_REAL const dzi = INV(dz); + CCTK_REAL const khalf = 0.5; + CCTK_REAL const kthird = 1/3.0; + CCTK_REAL const ktwothird = 2.0/3.0; + CCTK_REAL const kfourthird = 4.0/3.0; + CCTK_REAL const keightthird = 8.0/3.0; + CCTK_REAL const hdxi = 0.5 * dxi; + CCTK_REAL const hdyi = 0.5 * dyi; + CCTK_REAL const hdzi = 0.5 * dzi; + + /* Initialize predefined quantities */ + CCTK_REAL const p1o12dx = 0.0833333333333333333333333333333*INV(dx); + CCTK_REAL const p1o12dy = 0.0833333333333333333333333333333*INV(dy); + CCTK_REAL const p1o12dz = 0.0833333333333333333333333333333*INV(dz); + CCTK_REAL const p1o144dxdy = 0.00694444444444444444444444444444*INV(dx)*INV(dy); + CCTK_REAL const p1o144dxdz = 0.00694444444444444444444444444444*INV(dx)*INV(dz); + CCTK_REAL const p1o144dydz = 0.00694444444444444444444444444444*INV(dy)*INV(dz); + CCTK_REAL const pm1o12dx2 = -0.0833333333333333333333333333333*INV(SQR(dx)); + CCTK_REAL const pm1o12dy2 = -0.0833333333333333333333333333333*INV(SQR(dy)); + CCTK_REAL const pm1o12dz2 = -0.0833333333333333333333333333333*INV(SQR(dz)); + + /* Loop over the grid points */ + #pragma omp parallel + LC_LOOP3 (WT_Gaussian, + i,j,k, min[0],min[1],min[2], max[0],max[1],max[2], + cctk_lsh[0],cctk_lsh[1],cctk_lsh[2]) + { + ptrdiff_t const index = di*i + dj*j + dk*k; + + /* Assign local copies of grid functions */ + + /* Include user supplied include files */ + + /* Precompute derivatives */ + + /* Calculate temporaries and grid functions */ + CCTK_REAL uL = 0; + + CCTK_REAL rhoL = 0; + + + /* Copy local copies back to grid functions */ + rho[index] = rhoL; + u[index] = uL; + } + LC_ENDLOOP3 (WT_Gaussian); +} + +extern "C" void WT_Gaussian(CCTK_ARGUMENTS) +{ + DECLARE_CCTK_ARGUMENTS; + DECLARE_CCTK_PARAMETERS; + + GenericFD_LoopOverEverything(cctkGH, &WT_Gaussian_Body); +} diff --git a/ML_WaveToy/src/WT_RHS.cc b/ML_WaveToy/src/WT_RHS.cc new file mode 100644 index 0000000..3f1729b --- /dev/null +++ b/ML_WaveToy/src/WT_RHS.cc @@ -0,0 +1,130 @@ +/* File produced by Kranc */ + +#define KRANC_C + +#include +#include +#include +#include +#include "cctk.h" +#include "cctk_Arguments.h" +#include "cctk_Parameters.h" +#include "GenericFD.h" +#include "Differencing.h" +#include "loopcontrol.h" + +/* Define macros used in calculations */ +#define INITVALUE (42) +#define QAD(x) (SQR(SQR(x))) +#define INV(x) ((1.0) / (x)) +#define SQR(x) ((x) * (x)) +#define CUB(x) ((x) * (x) * (x)) + +extern "C" void WT_RHS_SelectBCs(CCTK_ARGUMENTS) +{ + DECLARE_CCTK_ARGUMENTS; + DECLARE_CCTK_PARAMETERS; + + CCTK_INT ierr = 0; + ierr = Boundary_SelectGroupForBC(cctkGH, CCTK_ALL_FACES, GenericFD_GetBoundaryWidth(cctkGH), -1 /* no table */, "ML_WaveToy::WT_rhorhs","flat"); + if (ierr < 0) + CCTK_WARN(1, "Failed to register flat BC for ML_WaveToy::WT_rhorhs."); + ierr = Boundary_SelectGroupForBC(cctkGH, CCTK_ALL_FACES, GenericFD_GetBoundaryWidth(cctkGH), -1 /* no table */, "ML_WaveToy::WT_urhs","flat"); + if (ierr < 0) + CCTK_WARN(1, "Failed to register flat BC for ML_WaveToy::WT_urhs."); + return; +} + +static void WT_RHS_Body(cGH const * restrict const cctkGH, int const dir, int const face, CCTK_REAL const normal[3], CCTK_REAL const tangentA[3], CCTK_REAL const tangentB[3], int const min[3], int const max[3], int const n_subblock_gfs, CCTK_REAL * restrict const subblock_gfs[]) +{ + DECLARE_CCTK_ARGUMENTS; + DECLARE_CCTK_PARAMETERS; + + + /* Declare finite differencing variables */ + + if (verbose > 1) + { + CCTK_VInfo(CCTK_THORNSTRING,"Entering WT_RHS_Body"); + } + + if (cctk_iteration % WT_RHS_calc_every != WT_RHS_calc_offset) + { + return; + } + + const char *groups[] = {"ML_WaveToy::WT_rho","ML_WaveToy::WT_rhorhs","ML_WaveToy::WT_u","ML_WaveToy::WT_urhs"}; + GenericFD_AssertGroupStorage(cctkGH, "WT_RHS", 4, groups); + + /* Include user-supplied include files */ + + /* Initialise finite differencing variables */ + ptrdiff_t const di = 1; + ptrdiff_t const dj = CCTK_GFINDEX3D(cctkGH,0,1,0) - CCTK_GFINDEX3D(cctkGH,0,0,0); + ptrdiff_t const dk = CCTK_GFINDEX3D(cctkGH,0,0,1) - CCTK_GFINDEX3D(cctkGH,0,0,0); + CCTK_REAL const dx = ToReal(CCTK_DELTA_SPACE(0)); + CCTK_REAL const dy = ToReal(CCTK_DELTA_SPACE(1)); + CCTK_REAL const dz = ToReal(CCTK_DELTA_SPACE(2)); + CCTK_REAL const dxi = INV(dx); + CCTK_REAL const dyi = INV(dy); + CCTK_REAL const dzi = INV(dz); + CCTK_REAL const khalf = 0.5; + CCTK_REAL const kthird = 1/3.0; + CCTK_REAL const ktwothird = 2.0/3.0; + CCTK_REAL const kfourthird = 4.0/3.0; + CCTK_REAL const keightthird = 8.0/3.0; + CCTK_REAL const hdxi = 0.5 * dxi; + CCTK_REAL const hdyi = 0.5 * dyi; + CCTK_REAL const hdzi = 0.5 * dzi; + + /* Initialize predefined quantities */ + CCTK_REAL const p1o12dx = 0.0833333333333333333333333333333*INV(dx); + CCTK_REAL const p1o12dy = 0.0833333333333333333333333333333*INV(dy); + CCTK_REAL const p1o12dz = 0.0833333333333333333333333333333*INV(dz); + CCTK_REAL const p1o144dxdy = 0.00694444444444444444444444444444*INV(dx)*INV(dy); + CCTK_REAL const p1o144dxdz = 0.00694444444444444444444444444444*INV(dx)*INV(dz); + CCTK_REAL const p1o144dydz = 0.00694444444444444444444444444444*INV(dy)*INV(dz); + CCTK_REAL const pm1o12dx2 = -0.0833333333333333333333333333333*INV(SQR(dx)); + CCTK_REAL const pm1o12dy2 = -0.0833333333333333333333333333333*INV(SQR(dy)); + CCTK_REAL const pm1o12dz2 = -0.0833333333333333333333333333333*INV(SQR(dz)); + + /* Loop over the grid points */ + #pragma omp parallel + LC_LOOP3 (WT_RHS, + i,j,k, min[0],min[1],min[2], max[0],max[1],max[2], + cctk_lsh[0],cctk_lsh[1],cctk_lsh[2]) + { + ptrdiff_t const index = di*i + dj*j + dk*k; + + /* Assign local copies of grid functions */ + CCTK_REAL rhoL = rho[index]; + CCTK_REAL uL = u[index]; + + /* Include user supplied include files */ + + /* Precompute derivatives */ + CCTK_REAL const PDstandardNth11u = PDstandardNth11(&u[index]); + CCTK_REAL const PDstandardNth22u = PDstandardNth22(&u[index]); + CCTK_REAL const PDstandardNth33u = PDstandardNth33(&u[index]); + + /* Calculate temporaries and grid functions */ + CCTK_REAL urhsL = rhoL; + + CCTK_REAL rhorhsL = PDstandardNth11u + PDstandardNth22u + + PDstandardNth33u; + + + /* Copy local copies back to grid functions */ + rhorhs[index] = rhorhsL; + urhs[index] = urhsL; + } + LC_ENDLOOP3 (WT_RHS); +} + +extern "C" void WT_RHS(CCTK_ARGUMENTS) +{ + DECLARE_CCTK_ARGUMENTS; + DECLARE_CCTK_PARAMETERS; + + GenericFD_LoopOverInterior(cctkGH, &WT_RHS_Body); +} -- cgit v1.2.3