From 0cc8f9d9bcace88c21f9385baa5fe789bab3f485 Mon Sep 17 00:00:00 2001 From: Lucas Sanches Date: Wed, 23 Oct 2024 15:00:38 -0500 Subject: [PATCH] TestRKAB: Code refactoring. Added noise initial data --- TestRKAB/par/gaussian.par | 1 + TestRKAB/par/noise.par | 79 ++++++++++++++++ TestRKAB/par/standing.par | 1 + TestRKAB/param.ccl | 17 +++- TestRKAB/schedule.ccl | 15 +-- TestRKAB/src/energy.cxx | 27 ++++++ TestRKAB/src/error.cxx | 85 +++++++++++++++++ TestRKAB/src/initial.cxx | 71 ++++++++++++++ TestRKAB/src/make.code.defn | 2 +- TestRKAB/src/rhs.cxx | 44 +++++++++ TestRKAB/src/testrkab.cxx | 181 ------------------------------------ 11 files changed, 334 insertions(+), 189 deletions(-) create mode 100644 TestRKAB/par/noise.par create mode 100644 TestRKAB/src/energy.cxx create mode 100644 TestRKAB/src/error.cxx create mode 100644 TestRKAB/src/initial.cxx create mode 100644 TestRKAB/src/rhs.cxx delete mode 100644 TestRKAB/src/testrkab.cxx diff --git a/TestRKAB/par/gaussian.par b/TestRKAB/par/gaussian.par index 1e590be6d..18233a6a2 100644 --- a/TestRKAB/par/gaussian.par +++ b/TestRKAB/par/gaussian.par @@ -25,6 +25,7 @@ $dt = $cfl * $dx $itlast = ceil($final_time / $dt) $out_every_it = ceil($out_every_time / $dt) +TestRKAB::compute_error = yes TestRKAB::initial_condition = "Gaussian" TestRKAB::gaussian_width = $cfl / sqrt(2) diff --git a/TestRKAB/par/noise.par b/TestRKAB/par/noise.par new file mode 100644 index 000000000..cc5636e83 --- /dev/null +++ b/TestRKAB/par/noise.par @@ -0,0 +1,79 @@ +ActiveThorns = " + CarpetX + IOUtil + ODESolvers + TestRKAB +" + +###### Settings ###### +$final_time = 1.0 +$out_every_time = 0.25 + +$start = -1.0 +$end = 1.0 + +$ncells = 160 +$cfl = 0.25 +###################### + +# We are using a 5 point centered stencil +$ghost_size = 2 + +$dx = ($end - $start)/$ncells +$dt = $cfl * $dx + +$itlast = ceil($final_time / $dt) +$out_every_it = ceil($out_every_time / $dt) + +TestRKAB::compute_error = no +TestRKAB::initial_condition = "noise" +TestRKAB::noise_seed = 100.0 +TestRKAB::noise_boundary = 1.0e-7 + +CarpetX::poison_undefined_values = no +CarpetX::verbose = no + +Cactus::cctk_show_schedule = yes +Cactus::presync_mode = "mixed-error" + +CarpetX::xmin = $start +CarpetX::ymin = $start +CarpetX::zmin = $start + +CarpetX::xmax = $end +CarpetX::ymax = $end +CarpetX::zmax = $end + +CarpetX::ncells_x = $ncells +CarpetX::ncells_y = $ncells +CarpetX::ncells_z = $ncells + +CarpetX::periodic = yes +CarpetX::periodic_x = yes +CarpetX::periodic_y = yes +CarpetX::periodic_z = yes + +CarpetX::ghost_size = $ghost_size + +CarpetX::dtfac = $cfl + +CarpetX::blocking_factor_x = 1 +CarpetX::blocking_factor_y = 1 +CarpetX::blocking_factor_z = 1 + +Cactus::terminate = "iteration" +Cactus::cctk_itlast = $itlast + +ODESolvers::method = "RKAB" +ODESolvers::verbose = yes + +IO::out_dir = $parfile +IO::out_every = $out_every_it + +CarpetX::out_norm_vars = "" +CarpetX::out_tsv_vars = "" + +CarpetX::out_silo_vars = " + TestRKAB::state + TestRKAB::energy +" diff --git a/TestRKAB/par/standing.par b/TestRKAB/par/standing.par index 99f60527b..73c2626ed 100644 --- a/TestRKAB/par/standing.par +++ b/TestRKAB/par/standing.par @@ -28,6 +28,7 @@ $dt = $cfl * $dx $itlast = ceil($final_time / $dt) $out_every_it = ceil($out_every_time / $dt) +TestRKAB::compute_error = yes TestRKAB::initial_condition = "standing wave" TestRKAB::amplitude = 1.0 TestRKAB::standing_wave_kx = $wave_k diff --git a/TestRKAB/param.ccl b/TestRKAB/param.ccl index b9344baf8..e2094f5f4 100644 --- a/TestRKAB/param.ccl +++ b/TestRKAB/param.ccl @@ -4,8 +4,13 @@ KEYWORD initial_condition "Initial condition" { "standing wave" :: "Standing wave" "Gaussian" :: "Gaussian" + "noise" :: "noise" } "standing wave" +BOOLEAN compute_error "Turn on error computation for Gaussian and standing wave initial data" +{ +} no + CCTK_REAL amplitude "Initial amplitude" { *:* :: "" @@ -29,4 +34,14 @@ CCTK_REAL standing_wave_kz "ky for standing wave" CCTK_REAL gaussian_width "width of Gaussian" { (0:* :: "" -} 1.0 \ No newline at end of file +} 1.0 + +CCTK_REAL noise_seed "Seed of the RNG for noise initial data" +{ + *:* :: "" +} 100.0 + +CCTK_REAL noise_boundary "Nois will be generated in the [-noise_boundary, noise_boundary] range" +{ + *:* :: "" +} 1.0e-7 \ No newline at end of file diff --git a/TestRKAB/schedule.ccl b/TestRKAB/schedule.ccl index 7e6e6cc33..a2f35fdd5 100644 --- a/TestRKAB/schedule.ccl +++ b/TestRKAB/schedule.ccl @@ -37,13 +37,16 @@ SCHEDULE TestRKAB_Sync IN ODESolvers_PostStep # Analysis -SCHEDULE TestRKAB_Error IN ODESolvers_PostStep AFTER TestRKAB_Sync +if (CCTK_EQUALS(initial_condition, "Gaussian") || CCTK_EQUALS(initial_condition, "standing wave")) { - LANG: C - READS: state(interior) - WRITES: error(interior) - SYNC: error -} "Calculate error in scalar wave state" + SCHEDULE TestRKAB_Error IN ODESolvers_PostStep AFTER TestRKAB_Sync + { + LANG: C + READS: state(interior) + WRITES: error(interior) + SYNC: error + } "Calculate error in scalar wave state" +} SCHEDULE TestRKAB_Energy IN ODESolvers_PostStep AFTER TestRKAB_Sync { diff --git a/TestRKAB/src/energy.cxx b/TestRKAB/src/energy.cxx new file mode 100644 index 000000000..f5c58a796 --- /dev/null +++ b/TestRKAB/src/energy.cxx @@ -0,0 +1,27 @@ +#include +#include +#include + +#include + +namespace TestRKAB { +using namespace Arith; + +extern "C" void TestRKAB_Energy(CCTK_ARGUMENTS) { + DECLARE_CCTK_ARGUMENTSX_TestRKAB_Energy; + DECLARE_CCTK_PARAMETERS; + + grid.loop_int_device<0, 0, 0>( + grid.nghostzones, + [=] CCTK_DEVICE(const Loop::PointDesc &p) CCTK_ATTRIBUTE_ALWAYS_INLINE { + const auto local_Pi{Pi(p.I)}; + const auto local_Dx{Dx(p.I)}; + const auto local_Dy{Dy(p.I)}; + const auto local_Dz{Dz(p.I)}; + + energy_density(p.I) = 0.5 * (local_Pi * local_Pi + local_Dx * local_Dx + + local_Dy * local_Dy + local_Dz * local_Dz); + }); +} + +} // namespace TestRKAB \ No newline at end of file diff --git a/TestRKAB/src/error.cxx b/TestRKAB/src/error.cxx new file mode 100644 index 000000000..5fbe37c41 --- /dev/null +++ b/TestRKAB/src/error.cxx @@ -0,0 +1,85 @@ +#include +#include +#include + +#include + +#include "standing_wave.hxx" +#include "gaussian.hxx" + +namespace TestRKAB { +using namespace Arith; + +extern "C" void TestRKAB_Error(CCTK_ARGUMENTS) { + DECLARE_CCTK_ARGUMENTSX_TestRKAB_Error; + DECLARE_CCTK_PARAMETERS; + + if (CCTK_EQUALS(initial_condition, "standing wave")) { + grid.loop_int_device<0, 0, 0>( + grid.nghostzones, + [=] CCTK_DEVICE(const Loop::PointDesc &p) CCTK_ATTRIBUTE_ALWAYS_INLINE { + using std::fabs; + + const auto t{cctk_time}; + + const auto A{amplitude}; + const auto kx{standing_wave_kx}; + const auto ky{standing_wave_ky}; + const auto kz{standing_wave_kz}; + + const auto expected_phi{sw::phi(A, kx, ky, kz, t, p.x, p.y, p.z)}; + const auto expected_Pi{sw::Pi(A, kx, ky, kz, t, p.x, p.y, p.z)}; + const auto expected_Dx{sw::Dx(A, kx, ky, kz, t, p.x, p.y, p.z)}; + const auto expected_Dy{sw::Dy(A, kx, ky, kz, t, p.x, p.y, p.z)}; + const auto expected_Dz{sw::Dz(A, kx, ky, kz, t, p.x, p.y, p.z)}; + + const auto actual_phi{phi(p.I)}; + const auto actual_Pi{Pi(p.I)}; + const auto actual_Dx{Dx(p.I)}; + const auto actual_Dy{Dy(p.I)}; + const auto actual_Dz{Dz(p.I)}; + + phi_err(p.I) = fabs(expected_phi - actual_phi); + Pi_err(p.I) = fabs(expected_Pi - actual_Pi); + Dx_err(p.I) = fabs(expected_Dx - actual_Dx); + Dy_err(p.I) = fabs(expected_Dy - actual_Dy); + Dz_err(p.I) = fabs(expected_Dz - actual_Dz); + }); + + } else if (CCTK_EQUALS(initial_condition, "Gaussian")) { + grid.loop_int_device<0, 0, 0>( + grid.nghostzones, + [=] CCTK_DEVICE(const Loop::PointDesc &p) CCTK_ATTRIBUTE_ALWAYS_INLINE { + using std::fabs; + + const auto t{cctk_time}; + + const auto expected_phi{ + gauss::phi(amplitude, gaussian_width, t, p.x, p.y, p.z)}; + const auto expected_Pi{ + gauss::Pi(amplitude, gaussian_width, t, p.x, p.y, p.z)}; + const auto expected_Dx{ + gauss::Dx(amplitude, gaussian_width, t, p.x, p.y, p.z)}; + const auto expected_Dy{ + gauss::Dy(amplitude, gaussian_width, t, p.x, p.y, p.z)}; + const auto expected_Dz{ + gauss::Dz(amplitude, gaussian_width, t, p.x, p.y, p.z)}; + + const auto actual_phi{phi(p.I)}; + const auto actual_Pi{Pi(p.I)}; + const auto actual_Dx{Dx(p.I)}; + const auto actual_Dy{Dy(p.I)}; + const auto actual_Dz{Dz(p.I)}; + + phi_err(p.I) = fabs(expected_phi - actual_phi); + Pi_err(p.I) = fabs(expected_Pi - actual_Pi); + Dx_err(p.I) = fabs(expected_Dx - actual_Dx); + Dy_err(p.I) = fabs(expected_Dy - actual_Dy); + Dz_err(p.I) = fabs(expected_Dz - actual_Dz); + }); + } else { + CCTK_VERROR("Unknown initial condition \"%s\"", initial_condition); + } +} + +} // namespace TestRKAB \ No newline at end of file diff --git a/TestRKAB/src/initial.cxx b/TestRKAB/src/initial.cxx new file mode 100644 index 000000000..a7f334f0e --- /dev/null +++ b/TestRKAB/src/initial.cxx @@ -0,0 +1,71 @@ +#include +#include +#include + +#include + +#include "standing_wave.hxx" +#include "gaussian.hxx" + +#include + +namespace TestRKAB { +using namespace Arith; + +extern "C" void TestRKAB_Initial(CCTK_ARGUMENTS) { + DECLARE_CCTK_ARGUMENTSX_TestRKAB_Initial; + DECLARE_CCTK_PARAMETERS; + + if (CCTK_EQUALS(initial_condition, "standing wave")) { + grid.loop_int_device<0, 0, 0>( + grid.nghostzones, + [=] CCTK_DEVICE(const Loop::PointDesc &p) CCTK_ATTRIBUTE_ALWAYS_INLINE { + const auto t{cctk_time}; + + const auto A{amplitude}; + const auto kx{standing_wave_kx}; + const auto ky{standing_wave_ky}; + const auto kz{standing_wave_kz}; + + phi(p.I) = sw::phi(A, kx, ky, kz, t, p.x, p.y, p.z); + Pi(p.I) = sw::Pi(A, kx, ky, kz, t, p.x, p.y, p.z); + Dx(p.I) = sw::Dx(A, kx, ky, kz, t, p.x, p.y, p.z); + Dy(p.I) = sw::Dy(A, kx, ky, kz, t, p.x, p.y, p.z); + Dz(p.I) = sw::Dz(A, kx, ky, kz, t, p.x, p.y, p.z); + }); + + } else if (CCTK_EQUALS(initial_condition, "Gaussian")) { + grid.loop_int_device<0, 0, 0>( + grid.nghostzones, + [=] CCTK_DEVICE(const Loop::PointDesc &p) CCTK_ATTRIBUTE_ALWAYS_INLINE { + const auto t{cctk_time}; + + phi(p.I) = gauss::phi(amplitude, gaussian_width, t, p.x, p.y, p.z); + Pi(p.I) = gauss::Pi(amplitude, gaussian_width, t, p.x, p.y, p.z); + Dx(p.I) = gauss::Dx(amplitude, gaussian_width, t, p.x, p.y, p.z); + Dy(p.I) = gauss::Dy(amplitude, gaussian_width, t, p.x, p.y, p.z); + Dz(p.I) = gauss::Dz(amplitude, gaussian_width, t, p.x, p.y, p.z); + }); + + } else if (CCTK_EQUALS(initial_condition, "noise")) { + std::uniform_real_distribution noise_distrib{-noise_boundary, + noise_boundary}; + std::mt19937 noise_engine{noise_seed}; + + grid.loop_int_device<0, 0, 0>(grid.nghostzones, + [&] CCTK_DEVICE(const Loop::PointDesc &p) + CCTK_ATTRIBUTE_ALWAYS_INLINE { + const auto t{cctk_time}; + + phi(p.I) = noise_distrib(noise_engine); + Pi(p.I) = noise_distrib(noise_engine); + Dx(p.I) = noise_distrib(noise_engine); + Dy(p.I) = noise_distrib(noise_engine); + Dz(p.I) = noise_distrib(noise_engine); + }); + } else { + CCTK_VERROR("Unknown initial condition \"%s\"", initial_condition); + } +} + +} // namespace TestRKAB \ No newline at end of file diff --git a/TestRKAB/src/make.code.defn b/TestRKAB/src/make.code.defn index ae2f57373..e40696f8c 100644 --- a/TestRKAB/src/make.code.defn +++ b/TestRKAB/src/make.code.defn @@ -1,7 +1,7 @@ # Main make.code.defn file for thorn TestSubcyclingMC # Source files in this directory -SRCS = testrkab.cxx +SRCS = energy.cxx error.cxx initial.cxx rhs.cxx # Subdirectories containing source files SUBDIRS = diff --git a/TestRKAB/src/rhs.cxx b/TestRKAB/src/rhs.cxx new file mode 100644 index 000000000..b3c6e04d9 --- /dev/null +++ b/TestRKAB/src/rhs.cxx @@ -0,0 +1,44 @@ + +#include +#include +#include + +#include + +namespace TestRKAB { +using namespace Arith; + +// Finite difference helpers +enum class fd_dir : std::size_t { x = 0, y = 1, z = 2 }; + +template +static inline auto fd_c_1_4(const Loop::PointDesc &p, + const Loop::GF3D2 &gf) noexcept -> T { + constexpr auto d{static_cast(direction)}; + const auto num{gf(p.I - 2 * p.DI[d]) - 8.0 * gf(p.I - 1 * p.DI[d]) + + 8.0 * gf(p.I + 1 * p.DI[d]) - 1.0 * gf(p.I + 2 * p.DI[d])}; + const auto den{1.0 / (12.0 * p.DX[d])}; + return num * den; +} + +extern "C" void TestRKAB_RHS(CCTK_ARGUMENTS) { + DECLARE_CCTK_ARGUMENTSX_TestRKAB_RHS; + DECLARE_CCTK_PARAMETERS; + + grid.loop_int_device<0, 0, 0>( + grid.nghostzones, + [=] CCTK_DEVICE(const Loop::PointDesc &p) CCTK_ATTRIBUTE_ALWAYS_INLINE { + phi_rhs(p.I) = Pi(p.I); + Pi_rhs(p.I) = fd_c_1_4(p, Dx) + fd_c_1_4(p, Dy) + + fd_c_1_4(p, Dz); + Dx_rhs(p.I) = fd_c_1_4(p, Pi); + Dy_rhs(p.I) = fd_c_1_4(p, Pi); + Dz_rhs(p.I) = fd_c_1_4(p, Pi); + }); +} + +extern "C" void TestRKAB_Sync(CCTK_ARGUMENTS) { + // do nothing +} + +} // namespace TestRKAB diff --git a/TestRKAB/src/testrkab.cxx b/TestRKAB/src/testrkab.cxx deleted file mode 100644 index 8f5abd2ef..000000000 --- a/TestRKAB/src/testrkab.cxx +++ /dev/null @@ -1,181 +0,0 @@ -#include - -#include -#include - -#include -#include -#include - -#include - -#include "standing_wave.hxx" -#include "gaussian.hxx" - -namespace TestRKAB { -using namespace Arith; - -// Finite difference helpers -enum class fd_dir : std::size_t { x = 0, y = 1, z = 2 }; - -template -static inline auto fd_c_1_4(const Loop::PointDesc &p, - const Loop::GF3D2 &gf) noexcept -> T { - constexpr auto d{static_cast(direction)}; - const auto num{gf(p.I - 2 * p.DI[d]) - 8.0 * gf(p.I - 1 * p.DI[d]) + - 8.0 * gf(p.I + 1 * p.DI[d]) - 1.0 * gf(p.I + 2 * p.DI[d])}; - const auto den{1.0 / (12.0 * p.DX[d])}; - return num * den; -} - -// Scheduled functions -extern "C" void TestRKAB_Initial(CCTK_ARGUMENTS) { - DECLARE_CCTK_ARGUMENTSX_TestRKAB_Initial; - DECLARE_CCTK_PARAMETERS; - - if (CCTK_EQUALS(initial_condition, "standing wave")) { - grid.loop_int_device<0, 0, 0>( - grid.nghostzones, - [=] CCTK_DEVICE(const Loop::PointDesc &p) CCTK_ATTRIBUTE_ALWAYS_INLINE { - const auto t{cctk_time}; - - const auto A{amplitude}; - const auto kx{standing_wave_kx}; - const auto ky{standing_wave_ky}; - const auto kz{standing_wave_kz}; - - phi(p.I) = sw::phi(A, kx, ky, kz, t, p.x, p.y, p.z); - Pi(p.I) = sw::Pi(A, kx, ky, kz, t, p.x, p.y, p.z); - Dx(p.I) = sw::Dx(A, kx, ky, kz, t, p.x, p.y, p.z); - Dy(p.I) = sw::Dy(A, kx, ky, kz, t, p.x, p.y, p.z); - Dz(p.I) = sw::Dz(A, kx, ky, kz, t, p.x, p.y, p.z); - }); - - } else if (CCTK_EQUALS(initial_condition, "Gaussian")) { - grid.loop_int_device<0, 0, 0>( - grid.nghostzones, - [=] CCTK_DEVICE(const Loop::PointDesc &p) CCTK_ATTRIBUTE_ALWAYS_INLINE { - const auto t{cctk_time}; - - phi(p.I) = gauss::phi(amplitude, gaussian_width, t, p.x, p.y, p.z); - Pi(p.I) = gauss::Pi(amplitude, gaussian_width, t, p.x, p.y, p.z); - Dx(p.I) = gauss::Dx(amplitude, gaussian_width, t, p.x, p.y, p.z); - Dy(p.I) = gauss::Dy(amplitude, gaussian_width, t, p.x, p.y, p.z); - Dz(p.I) = gauss::Dz(amplitude, gaussian_width, t, p.x, p.y, p.z); - }); - - } else { - CCTK_VERROR("Unknown initial condition \"%s\"", initial_condition); - } -} - -extern "C" void TestRKAB_RHS(CCTK_ARGUMENTS) { - DECLARE_CCTK_ARGUMENTSX_TestRKAB_RHS; - DECLARE_CCTK_PARAMETERS; - - grid.loop_int_device<0, 0, 0>( - grid.nghostzones, - [=] CCTK_DEVICE(const Loop::PointDesc &p) CCTK_ATTRIBUTE_ALWAYS_INLINE { - phi_rhs(p.I) = Pi(p.I); - Pi_rhs(p.I) = fd_c_1_4(p, Dx) + fd_c_1_4(p, Dy) + - fd_c_1_4(p, Dz); - Dx_rhs(p.I) = fd_c_1_4(p, Pi); - Dy_rhs(p.I) = fd_c_1_4(p, Pi); - Dz_rhs(p.I) = fd_c_1_4(p, Pi); - }); -} - -extern "C" void TestRKAB_Sync(CCTK_ARGUMENTS) { - // do nothing -} - -extern "C" void TestRKAB_Error(CCTK_ARGUMENTS) { - DECLARE_CCTK_ARGUMENTSX_TestRKAB_Error; - DECLARE_CCTK_PARAMETERS; - - if (CCTK_EQUALS(initial_condition, "standing wave")) { - grid.loop_int_device<0, 0, 0>( - grid.nghostzones, - [=] CCTK_DEVICE(const Loop::PointDesc &p) CCTK_ATTRIBUTE_ALWAYS_INLINE { - using std::fabs; - - const auto t{cctk_time}; - - const auto A{amplitude}; - const auto kx{standing_wave_kx}; - const auto ky{standing_wave_ky}; - const auto kz{standing_wave_kz}; - - const auto expected_phi{sw::phi(A, kx, ky, kz, t, p.x, p.y, p.z)}; - const auto expected_Pi{sw::Pi(A, kx, ky, kz, t, p.x, p.y, p.z)}; - const auto expected_Dx{sw::Dx(A, kx, ky, kz, t, p.x, p.y, p.z)}; - const auto expected_Dy{sw::Dy(A, kx, ky, kz, t, p.x, p.y, p.z)}; - const auto expected_Dz{sw::Dz(A, kx, ky, kz, t, p.x, p.y, p.z)}; - - const auto actual_phi{phi(p.I)}; - const auto actual_Pi{Pi(p.I)}; - const auto actual_Dx{Dx(p.I)}; - const auto actual_Dy{Dy(p.I)}; - const auto actual_Dz{Dz(p.I)}; - - phi_err(p.I) = fabs(expected_phi - actual_phi); - Pi_err(p.I) = fabs(expected_Pi - actual_Pi); - Dx_err(p.I) = fabs(expected_Dx - actual_Dx); - Dy_err(p.I) = fabs(expected_Dy - actual_Dy); - Dz_err(p.I) = fabs(expected_Dz - actual_Dz); - }); - - } else if (CCTK_EQUALS(initial_condition, "Gaussian")) { - grid.loop_int_device<0, 0, 0>( - grid.nghostzones, - [=] CCTK_DEVICE(const Loop::PointDesc &p) CCTK_ATTRIBUTE_ALWAYS_INLINE { - using std::fabs; - - const auto t{cctk_time}; - - const auto expected_phi{ - gauss::phi(amplitude, gaussian_width, t, p.x, p.y, p.z)}; - const auto expected_Pi{ - gauss::Pi(amplitude, gaussian_width, t, p.x, p.y, p.z)}; - const auto expected_Dx{ - gauss::Dx(amplitude, gaussian_width, t, p.x, p.y, p.z)}; - const auto expected_Dy{ - gauss::Dy(amplitude, gaussian_width, t, p.x, p.y, p.z)}; - const auto expected_Dz{ - gauss::Dz(amplitude, gaussian_width, t, p.x, p.y, p.z)}; - - const auto actual_phi{phi(p.I)}; - const auto actual_Pi{Pi(p.I)}; - const auto actual_Dx{Dx(p.I)}; - const auto actual_Dy{Dy(p.I)}; - const auto actual_Dz{Dz(p.I)}; - - phi_err(p.I) = fabs(expected_phi - actual_phi); - Pi_err(p.I) = fabs(expected_Pi - actual_Pi); - Dx_err(p.I) = fabs(expected_Dx - actual_Dx); - Dy_err(p.I) = fabs(expected_Dy - actual_Dy); - Dz_err(p.I) = fabs(expected_Dz - actual_Dz); - }); - } else { - CCTK_VERROR("Unknown initial condition \"%s\"", initial_condition); - } -} - -extern "C" void TestRKAB_Energy(CCTK_ARGUMENTS) { - DECLARE_CCTK_ARGUMENTSX_TestRKAB_Energy; - DECLARE_CCTK_PARAMETERS; - - grid.loop_int_device<0, 0, 0>( - grid.nghostzones, - [=] CCTK_DEVICE(const Loop::PointDesc &p) CCTK_ATTRIBUTE_ALWAYS_INLINE { - const auto local_Pi{Pi(p.I)}; - const auto local_Dx{Dx(p.I)}; - const auto local_Dy{Dy(p.I)}; - const auto local_Dz{Dz(p.I)}; - - energy_density(p.I) = 0.5 * (local_Pi * local_Pi + local_Dx * local_Dx + - local_Dy * local_Dy + local_Dz * local_Dz); - }); -} - -} // namespace TestRKAB