diff --git a/TestRKAB/par/gaussian.par b/TestRKAB/par/gaussian.par index a8e9f4f6d..1e590be6d 100644 --- a/TestRKAB/par/gaussian.par +++ b/TestRKAB/par/gaussian.par @@ -5,14 +5,28 @@ ActiveThorns = " TestRKAB " +###### Settings ###### +$final_time = 10.0 +$out_every_time = 0.25 + +$start = -1.0 +$end = 1.0 + $ncells = 80 $cfl = 0.25 +###################### + +# We are using a 5 point centered stencil +$ghost_size = 2 -$itlast = 10 -$out_every = 1 +$dx = ($end - $start)/$ncells +$dt = $cfl * $dx + +$itlast = ceil($final_time / $dt) +$out_every_it = ceil($out_every_time / $dt) TestRKAB::initial_condition = "Gaussian" -TestRKAB::gaussian_width = 0.17677669529 # sqrt(2)*W = 0.25 +TestRKAB::gaussian_width = $cfl / sqrt(2) CarpetX::poison_undefined_values = no CarpetX::verbose = no @@ -20,13 +34,13 @@ CarpetX::verbose = no Cactus::cctk_show_schedule = yes Cactus::presync_mode = "mixed-error" -CarpetX::xmin = -1.0 -CarpetX::ymin = -1.0 -CarpetX::zmin = -1.0 +CarpetX::xmin = $start +CarpetX::ymin = $start +CarpetX::zmin = $start -CarpetX::xmax = +1.0 -CarpetX::ymax = +1.0 -CarpetX::zmax = +1.0 +CarpetX::xmax = $end +CarpetX::ymax = $end +CarpetX::zmax = $end CarpetX::ncells_x = $ncells CarpetX::ncells_y = $ncells @@ -37,15 +51,14 @@ CarpetX::periodic_x = yes CarpetX::periodic_y = yes CarpetX::periodic_z = yes -#CarpetX::ghost_size = 3 +CarpetX::ghost_size = $ghost_size + CarpetX::dtfac = $cfl -CarpetX::blocking_factor_x = 2 -CarpetX::blocking_factor_y = 2 -CarpetX::blocking_factor_z = 2 +CarpetX::blocking_factor_x = 1 +CarpetX::blocking_factor_y = 1 +CarpetX::blocking_factor_z = 1 -#Cactus::terminate = "time" -#Cactus::cctk_final_time = 1.0 Cactus::terminate = "iteration" Cactus::cctk_itlast = $itlast @@ -53,14 +66,13 @@ ODESolvers::method = "RKAB" ODESolvers::verbose = yes IO::out_dir = $parfile -IO::out_every = $out_every +IO::out_every = $out_every_it CarpetX::out_norm_vars = "" CarpetX::out_tsv_vars = "" CarpetX::out_silo_vars = " TestRKAB::state - TestRKAB::pre TestRKAB::error TestRKAB::energy " diff --git a/TestRKAB/par/standing.par b/TestRKAB/par/standing.par index 6ee087d81..99f60527b 100644 --- a/TestRKAB/par/standing.par +++ b/TestRKAB/par/standing.par @@ -5,17 +5,34 @@ ActiveThorns = " TestRKAB " +###### Settings ###### +$final_time = 35.0 +$out_every_time = 1.0 + $ncells = 80 $cfl = 0.25 +###################### + +# We are using a 5 point centered stencil +$ghost_size = 2 + +# Make sure that the distance between the first left ghost and the right boundary is 2 pi +$boundary = ($ncells * $pi) / (1.0 + $ncells) + +# Make sure that we fit one period in between the boundaries +$wave_k = ($ncells + 1) / ($ncells * $pi) -$itlast = 10 -$out_every = 1 +$dx = 2.0 * $boundary / $ncells +$dt = $cfl * $dx + +$itlast = ceil($final_time / $dt) +$out_every_it = ceil($out_every_time / $dt) TestRKAB::initial_condition = "standing wave" TestRKAB::amplitude = 1.0 -TestRKAB::standing_wave_kx = 1.0 -TestRKAB::standing_wave_ky = 1.0 -TestRKAB::standing_wave_kz = 1.0 +TestRKAB::standing_wave_kx = $wave_k +TestRKAB::standing_wave_ky = $wave_k +TestRKAB::standing_wave_kz = $wave_k CarpetX::poison_undefined_values = no CarpetX::verbose = no @@ -23,13 +40,13 @@ CarpetX::verbose = no Cactus::cctk_show_schedule = yes Cactus::presync_mode = "mixed-error" -CarpetX::xmin = -1.0 -CarpetX::ymin = -1.0 -CarpetX::zmin = -1.0 +CarpetX::xmin = -$boundary +CarpetX::ymin = -$boundary +CarpetX::zmin = -$boundary -CarpetX::xmax = +1.0 -CarpetX::ymax = +1.0 -CarpetX::zmax = +1.0 +CarpetX::xmax = $boundary +CarpetX::ymax = $boundary +CarpetX::zmax = $boundary CarpetX::ncells_x = $ncells CarpetX::ncells_y = $ncells @@ -40,15 +57,14 @@ CarpetX::periodic_x = yes CarpetX::periodic_y = yes CarpetX::periodic_z = yes -CarpetX::ghost_size = 3 +CarpetX::ghost_size = $ghost_size + CarpetX::dtfac = $cfl -CarpetX::blocking_factor_x = 2 -CarpetX::blocking_factor_y = 2 -CarpetX::blocking_factor_z = 2 +CarpetX::blocking_factor_x = 1 +CarpetX::blocking_factor_y = 1 +CarpetX::blocking_factor_z = 1 -#Cactus::terminate = "time" -#Cactus::cctk_final_time = 1.0 Cactus::terminate = "iteration" Cactus::cctk_itlast = $itlast @@ -56,13 +72,13 @@ ODESolvers::method = "RKAB" ODESolvers::verbose = yes IO::out_dir = $parfile -IO::out_every = $out_every +IO::out_every = $out_every_it +CarpetX::out_norm_vars = "" CarpetX::out_tsv_vars = "" CarpetX::out_silo_vars = " TestRKAB::state - TestRKAB::pre TestRKAB::error TestRKAB::energy -" +" \ No newline at end of file diff --git a/TestRKAB/par/time.par b/TestRKAB/par/time.par deleted file mode 100644 index c5427269f..000000000 --- a/TestRKAB/par/time.par +++ /dev/null @@ -1,64 +0,0 @@ -ActiveThorns = " - CarpetX - IOUtil - ODESolvers - TestRKAB -" - -$ncells = 80 -$cfl = 0.25 - -$itlast = 10 -$out_every = 1 - -TestRKAB::initial_condition = "time" - -CarpetX::poison_undefined_values = no -CarpetX::verbose = no - -Cactus::cctk_show_schedule = yes -Cactus::presync_mode = "mixed-error" - -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 = yes -CarpetX::periodic_x = yes -CarpetX::periodic_y = yes -CarpetX::periodic_z = yes - -CarpetX::ghost_size = 3 -CarpetX::dtfac = $cfl - -CarpetX::blocking_factor_x = 2 -CarpetX::blocking_factor_y = 2 -CarpetX::blocking_factor_z = 2 - -#Cactus::terminate = "time" -#Cactus::cctk_final_time = 1.0 -Cactus::terminate = "iteration" -Cactus::cctk_itlast = $itlast - -ODESolvers::method = "RKAB" -ODESolvers::verbose = yes - -IO::out_dir = $parfile -IO::out_every = $out_every - -CarpetX::out_tsv_vars = "" - -CarpetX::out_silo_vars = " - TestRKAB::state - TestRKAB::pre - TestRKAB::error - TestRKAB::energy -" diff --git a/TestRKAB/param.ccl b/TestRKAB/param.ccl index 0bc16ed31..b9344baf8 100644 --- a/TestRKAB/param.ccl +++ b/TestRKAB/param.ccl @@ -4,7 +4,6 @@ KEYWORD initial_condition "Initial condition" { "standing wave" :: "Standing wave" "Gaussian" :: "Gaussian" - "time" :: "Simulation time" } "standing wave" CCTK_REAL amplitude "Initial amplitude" diff --git a/TestRKAB/src/gaussian.hxx b/TestRKAB/src/gaussian.hxx new file mode 100644 index 000000000..e330a5501 --- /dev/null +++ b/TestRKAB/src/gaussian.hxx @@ -0,0 +1,125 @@ +#ifndef TEST_RKAB_GAUSSIAN_HXX +#define TEST_RKAB_GAUSSIAN_HXX + +#include + +namespace TestRKAB::gauss { + +template +static inline auto phi(T W, T A, T t, T x, T y, T z) noexcept -> T { + using std::sqrt, std::cosh, std::sinh; + + const auto r{sqrt(x * x + y * y + z * z)}; + + if (r < sqrt(std::numeric_limits::epsilon())) { + return (2 * A * t) / (exp(pow(t, 2) / (2. * pow(W, 2))) * pow(W, 2)); + } else { + return (A * + (exp(-0.5 * pow(t - sqrt(pow(x, 2) + pow(y, 2) + pow(z, 2)), 2) / + pow(W, 2)) - + exp(-0.5 * pow(t + sqrt(pow(x, 2) + pow(y, 2) + pow(z, 2)), 2) / + pow(W, 2)))) / + sqrt(pow(x, 2) + pow(y, 2) + pow(z, 2)); + } +} + +template +static inline auto Pi(T W, T A, T t, T x, T y, T z) noexcept -> T { + using std::sqrt, std::cosh, std::sinh; + + const auto r{sqrt(x * x + y * y + z * z)}; + + if (r < sqrt(std::numeric_limits::epsilon())) { + return (2 * A * (-t + W) * (t + W)) / + (exp(pow(t, 2) / (2. * pow(W, 2))) * pow(W, 4)); + } else { + return (2 * A * + (sqrt(pow(x, 2) + pow(y, 2) + pow(z, 2)) * + cosh((t * sqrt(pow(x, 2) + pow(y, 2) + pow(z, 2))) / + pow(W, 2)) - + t * sinh((t * sqrt(pow(x, 2) + pow(y, 2) + pow(z, 2))) / + pow(W, 2)))) / + (exp((pow(t, 2) + pow(x, 2) + pow(y, 2) + pow(z, 2)) / + (2. * pow(W, 2))) * + pow(W, 2) * sqrt(pow(x, 2) + pow(y, 2) + pow(z, 2))); + } +} + +template +static inline auto Dx(T W, T A, T t, T x, T y, T z) noexcept -> T { + using std::sqrt, std::cosh, std::sinh; + + const auto r{sqrt(x * x + y * y + z * z)}; + + if (r < sqrt(std::numeric_limits::epsilon())) { + return 0; + } else { + return (A * x * + (pow(W, 2) + pow(x, 2) + pow(y, 2) + pow(z, 2) + + t * sqrt(pow(x, 2) + pow(y, 2) + pow(z, 2)) + + exp((2 * t * sqrt(pow(x, 2) + pow(y, 2) + pow(z, 2))) / + pow(W, 2)) * + t * sqrt(pow(x, 2) + pow(y, 2) + pow(z, 2)) - + exp((2 * t * sqrt(pow(x, 2) + pow(y, 2) + pow(z, 2))) / + pow(W, 2)) * + (pow(W, 2) + pow(x, 2) + pow(y, 2) + pow(z, 2)))) / + (exp((pow(t, 2) + pow(x, 2) + pow(y, 2) + pow(z, 2) + + 2 * t * sqrt(pow(x, 2) + pow(y, 2) + pow(z, 2))) / + (2. * pow(W, 2))) * + pow(W, 2) * pow(pow(x, 2) + pow(y, 2) + pow(z, 2), 1.5)); + } +} + +template +static inline auto Dy(T W, T A, T t, T x, T y, T z) noexcept -> T { + using std::sqrt, std::cosh, std::sinh; + + const auto r{sqrt(x * x + y * y + z * z)}; + + if (r < sqrt(std::numeric_limits::epsilon())) { + return 0; + } else { + return (A * y * + (pow(W, 2) + pow(x, 2) + pow(y, 2) + pow(z, 2) + + t * sqrt(pow(x, 2) + pow(y, 2) + pow(z, 2)) + + exp((2 * t * sqrt(pow(x, 2) + pow(y, 2) + pow(z, 2))) / + pow(W, 2)) * + t * sqrt(pow(x, 2) + pow(y, 2) + pow(z, 2)) - + exp((2 * t * sqrt(pow(x, 2) + pow(y, 2) + pow(z, 2))) / + pow(W, 2)) * + (pow(W, 2) + pow(x, 2) + pow(y, 2) + pow(z, 2)))) / + (exp((pow(t, 2) + pow(x, 2) + pow(y, 2) + pow(z, 2) + + 2 * t * sqrt(pow(x, 2) + pow(y, 2) + pow(z, 2))) / + (2. * pow(W, 2))) * + pow(W, 2) * pow(pow(x, 2) + pow(y, 2) + pow(z, 2), 1.5)); + } +} + +template +static inline auto Dz(T W, T A, T t, T x, T y, T z) noexcept -> T { + using std::sqrt, std::cosh, std::sinh; + + const auto r{sqrt(x * x + y * y + z * z)}; + + if (r < sqrt(std::numeric_limits::epsilon())) { + return 0; + } else { + return (A * z * + (pow(W, 2) + pow(x, 2) + pow(y, 2) + pow(z, 2) + + t * sqrt(pow(x, 2) + pow(y, 2) + pow(z, 2)) + + exp((2 * t * sqrt(pow(x, 2) + pow(y, 2) + pow(z, 2))) / + pow(W, 2)) * + t * sqrt(pow(x, 2) + pow(y, 2) + pow(z, 2)) - + exp((2 * t * sqrt(pow(x, 2) + pow(y, 2) + pow(z, 2))) / + pow(W, 2)) * + (pow(W, 2) + pow(x, 2) + pow(y, 2) + pow(z, 2)))) / + (exp((pow(t, 2) + pow(x, 2) + pow(y, 2) + pow(z, 2) + + 2 * t * sqrt(pow(x, 2) + pow(y, 2) + pow(z, 2))) / + (2. * pow(W, 2))) * + pow(W, 2) * pow(pow(x, 2) + pow(y, 2) + pow(z, 2), 1.5)); + } +} + +} // namespace TestRKAB::gauss + +#endif // TEST_RKAB_GAUSSIAN_HXX diff --git a/TestRKAB/src/standing_wave.hxx b/TestRKAB/src/standing_wave.hxx new file mode 100644 index 000000000..24bf4e19c --- /dev/null +++ b/TestRKAB/src/standing_wave.hxx @@ -0,0 +1,56 @@ +#ifndef TEST_RKAB_STANDING_WAVE_HXX +#define TEST_RKAB_STANDING_WAVE_HXX + +#include + +namespace TestRKAB::sw { + +template +static inline auto phi(T A, T kx, T ky, T kz, T t, T x, T y, + T z) noexcept -> T { + using std::sqrt, std::sin, std::cos, std::acos; + const auto pi{acos(T{-1})}; + const auto omega{sqrt(kx * kx + ky * ky + kz * kz)}; + return A * cos(2 * omega * pi * t) * cos(2 * kx * pi * x) * + cos(2 * ky * pi * y) * cos(2 * kz * pi * z); +} + +template +static inline auto Pi(T A, T kx, T ky, T kz, T t, T x, T y, T z) noexcept -> T { + using std::sqrt, std::sin, std::cos, std::acos; + const auto pi{acos(T{-1})}; + const auto omega{sqrt(kx * kx + ky * ky + kz * kz)}; + return -2 * A * omega * pi * cos(2 * kx * pi * x) * cos(2 * ky * pi * y) * + cos(2 * kz * pi * z) * sin(2 * omega * pi * t); +} + +template +static inline auto Dx(T A, T kx, T ky, T kz, T t, T x, T y, T z) noexcept -> T { + using std::sqrt, std::sin, std::cos, std::acos; + const auto pi{acos(T{-1})}; + const auto omega{sqrt(kx * kx + ky * ky + kz * kz)}; + return -2 * A * kx * pi * cos(2 * omega * pi * t) * cos(2 * ky * pi * y) * + cos(2 * kz * pi * z) * sin(2 * kx * pi * x); +} + +template +static inline auto Dy(T A, T kx, T ky, T kz, T t, T x, T y, T z) noexcept -> T { + using std::sqrt, std::sin, std::cos, std::acos; + const auto pi{acos(T{-1})}; + const auto omega{sqrt(kx * kx + ky * ky + kz * kz)}; + return -2 * A * ky * pi * cos(2 * omega * pi * t) * cos(2 * kx * pi * x) * + cos(2 * kz * pi * z) * sin(2 * ky * pi * y); +} + +template +static inline auto Dz(T A, T kx, T ky, T kz, T t, T x, T y, T z) noexcept -> T { + using std::sqrt, std::sin, std::cos, std::acos; + const auto pi{acos(T{-1})}; + const auto omega{sqrt(kx * kx + ky * ky + kz * kz)}; + return -2 * A * kz * pi * cos(2 * omega * pi * t) * cos(2 * kx * pi * x) * + cos(2 * ky * pi * y) * sin(2 * kz * pi * z); +} + +} // namespace TestRKAB::sw + +#endif // TEST_RKAB_STANDING_WAVE_HXX \ No newline at end of file diff --git a/TestRKAB/src/testrkab.cxx b/TestRKAB/src/testrkab.cxx index 3615dbca6..8f5abd2ef 100644 --- a/TestRKAB/src/testrkab.cxx +++ b/TestRKAB/src/testrkab.cxx @@ -1,5 +1,3 @@ -#include -#include #include #include @@ -9,179 +7,14 @@ #include #include -#include -#include -#include #include +#include "standing_wave.hxx" +#include "gaussian.hxx" + namespace TestRKAB { using namespace Arith; -constexpr int dim = 3; - -// Standing wave functions -template -static inline auto sw_phi(T A, T kx, T ky, T kz, T t, T x, T y, - T z) noexcept -> T { - using std::sqrt, std::sin, std::cos; - const auto pi{acos(T{-1})}; - const auto omega{sqrt(kx * kx + ky * ky + kz * kz)}; - return A * cos(2 * omega * pi * t) * sin(2 * kx * pi * x) * - sin(2 * ky * pi * y) * sin(2 * kz * pi * z); -} - -template -static inline auto sw_Pi(T A, T kx, T ky, T kz, T t, T x, T y, - T z) noexcept -> T { - using std::sqrt, std::sin, std::cos; - const auto pi{acos(T{-1})}; - const auto omega{sqrt(kx * kx + ky * ky + kz * kz)}; - return -2 * A * omega * pi * sin(2 * omega * pi * t) * sin(2 * kx * pi * x) * - sin(2 * ky * pi * y) * sin(2 * kz * pi * z); -} - -template -static inline auto sw_Dx(T A, T kx, T ky, T kz, T t, T x, T y, - T z) noexcept -> T { - using std::sqrt, std::sin, std::cos; - const auto pi{acos(T{-1})}; - const auto omega{sqrt(kx * kx + ky * ky + kz * kz)}; - return 2 * A * kx * pi * cos(2 * omega * pi * t) * cos(2 * kx * pi * x) * - sin(2 * ky * pi * y) * sin(2 * kz * pi * z); -} - -template -static inline auto sw_Dy(T A, T kx, T ky, T kz, T t, T x, T y, - T z) noexcept -> T { - using std::sqrt, std::sin, std::cos; - const auto pi{acos(T{-1})}; - const auto omega{sqrt(kx * kx + ky * ky + kz * kz)}; - return 2 * A * ky * pi * cos(2 * omega * pi * t) * sin(2 * kx * pi * x) * - cos(2 * ky * pi * y) * sin(2 * kz * pi * z); -} - -template -static inline auto sw_Dz(T A, T kx, T ky, T kz, T t, T x, T y, - T z) noexcept -> T { - using std::sqrt, std::sin, std::cos; - const auto pi{acos(T{-1})}; - const auto omega{sqrt(kx * kx + ky * ky + kz * kz)}; - return 2 * A * kz * pi * cos(2 * omega * pi * t) * sin(2 * kx * pi * x) * - sin(2 * ky * pi * y) * cos(2 * kz * pi * z); -} - -// Gaussian functions -template -static inline auto gaussian_phi(T A, T W, T t, T x, T y, T z) noexcept -> T { - using std::sqrt, std::exp, std::pow, std::sinh, std::cosh; - - const auto tol{sqrt(std::numeric_limits::epsilon())}; - const auto r{sqrt(x * x + y * y + z * z)}; - - if (r < tol) { - return (2 * A * t) / (exp(pow(t, 2) / (2. * pow(W, 2))) * pow(W, 2)); - } else { - return (A * - (exp(-0.5 * pow(t - sqrt(pow(x, 2) + pow(y, 2) + pow(z, 2)), 2) / - pow(W, 2)) - - exp(-0.5 * pow(t + sqrt(pow(x, 2) + pow(y, 2) + pow(z, 2)), 2) / - pow(W, 2)))) / - sqrt(pow(x, 2) + pow(y, 2) + pow(z, 2)); - } -} - -template -static inline auto gaussian_Pi(T A, T W, T t, T x, T y, T z) noexcept -> T { - using std::sqrt, std::exp, std::pow, std::sinh, std::cosh; - - const auto tol{sqrt(std::numeric_limits::epsilon())}; - const auto r{sqrt(x * x + y * y + z * z)}; - - if (r < tol) { - return (2 * A * (-t + W) * (t + W)) / - (exp(pow(t, 2) / (2. * pow(W, 2))) * pow(W, 4)); - } else { - return (2 * A * - (sqrt(pow(x, 2) + pow(y, 2) + pow(z, 2)) * - cosh((t * sqrt(pow(x, 2) + pow(y, 2) + pow(z, 2))) / - pow(W, 2)) - - t * sinh((t * sqrt(pow(x, 2) + pow(y, 2) + pow(z, 2))) / - pow(W, 2)))) / - (exp((pow(t, 2) + pow(x, 2) + pow(y, 2) + pow(z, 2)) / - (2. * pow(W, 2))) * - pow(W, 2) * sqrt(pow(x, 2) + pow(y, 2) + pow(z, 2))); - } -} - -template -static inline auto gaussian_Dx(T A, T W, T t, T x, T y, T z) noexcept -> T { - using std::sqrt, std::exp, std::pow, std::sinh, std::cosh; - - const auto tol{sqrt(std::numeric_limits::epsilon())}; - const auto r{sqrt(x * x + y * y + z * z)}; - - if (r < tol) { - return 0; - } else { - return (A * x * - (2 * t * sqrt(pow(x, 2) + pow(y, 2) + pow(z, 2)) * - cosh((t * sqrt(pow(x, 2) + pow(y, 2) + pow(z, 2))) / - pow(W, 2)) - - 2 * (pow(W, 2) + pow(x, 2) + pow(y, 2) + pow(z, 2)) * - sinh((t * sqrt(pow(x, 2) + pow(y, 2) + pow(z, 2))) / - pow(W, 2)))) / - (exp((pow(t, 2) + pow(x, 2) + pow(y, 2) + pow(z, 2)) / - (2. * pow(W, 2))) * - pow(W, 2) * pow(pow(x, 2) + pow(y, 2) + pow(z, 2), 1.5)); - } -} - -template -static inline auto gaussian_Dy(T A, T W, T t, T x, T y, T z) noexcept -> T { - using std::sqrt, std::exp, std::pow, std::sinh, std::cosh; - - const auto tol{sqrt(std::numeric_limits::epsilon())}; - const auto r{sqrt(x * x + y * y + z * z)}; - - if (r < tol) { - return A; - } else { - return (A * y * - (2 * t * sqrt(pow(x, 2) + pow(y, 2) + pow(z, 2)) * - cosh((t * sqrt(pow(x, 2) + pow(y, 2) + pow(z, 2))) / - pow(W, 2)) - - 2 * (pow(W, 2) + pow(x, 2) + pow(y, 2) + pow(z, 2)) * - sinh((t * sqrt(pow(x, 2) + pow(y, 2) + pow(z, 2))) / - pow(W, 2)))) / - (exp((pow(t, 2) + pow(x, 2) + pow(y, 2) + pow(z, 2)) / - (2. * pow(W, 2))) * - pow(W, 2) * pow(pow(x, 2) + pow(y, 2) + pow(z, 2), 1.5)); - } -} - -template -static inline auto gaussian_Dz(T A, T W, T t, T x, T y, T z) noexcept -> T { - using std::sqrt, std::exp, std::pow, std::sinh, std::cosh; - - const auto tol{sqrt(std::numeric_limits::epsilon())}; - const auto r{sqrt(x * x + y * y + z * z)}; - - if (r < tol) { - return A; - } else { - return (A * z * - (2 * t * sqrt(pow(x, 2) + pow(y, 2) + pow(z, 2)) * - cosh((t * sqrt(pow(x, 2) + pow(y, 2) + pow(z, 2))) / - pow(W, 2)) - - 2 * (pow(W, 2) + pow(x, 2) + pow(y, 2) + pow(z, 2)) * - sinh((t * sqrt(pow(x, 2) + pow(y, 2) + pow(z, 2))) / - pow(W, 2)))) / - (exp((pow(t, 2) + pow(x, 2) + pow(y, 2) + pow(z, 2)) / - (2. * pow(W, 2))) * - pow(W, 2) * pow(pow(x, 2) + pow(y, 2) + pow(z, 2), 1.5)); - } -} - // Finite difference helpers enum class fd_dir : std::size_t { x = 0, y = 1, z = 2 }; @@ -192,7 +25,7 @@ static inline auto fd_c_1_4(const Loop::PointDesc &p, 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 den * num; + return num * den; } // Scheduled functions @@ -205,24 +38,17 @@ extern "C" void TestRKAB_Initial(CCTK_ARGUMENTS) { grid.nghostzones, [=] CCTK_DEVICE(const Loop::PointDesc &p) CCTK_ATTRIBUTE_ALWAYS_INLINE { const auto t{cctk_time}; - const auto t_pre{t - cctk_delta_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); - - phi_pre(p.I) = sw_phi(A, kx, ky, kz, t_pre, p.x, p.y, p.z); - Pi_pre(p.I) = sw_Pi(A, kx, ky, kz, t_pre, p.x, p.y, p.z); - Dx_pre(p.I) = sw_Dx(A, kx, ky, kz, t_pre, p.x, p.y, p.z); - Dy_pre(p.I) = sw_Dy(A, kx, ky, kz, t_pre, p.x, p.y, p.z); - Dz_pre(p.I) = sw_Dz(A, kx, ky, kz, t_pre, p.x, p.y, p.z); + 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")) { @@ -230,45 +56,16 @@ extern "C" void TestRKAB_Initial(CCTK_ARGUMENTS) { grid.nghostzones, [=] CCTK_DEVICE(const Loop::PointDesc &p) CCTK_ATTRIBUTE_ALWAYS_INLINE { const auto t{cctk_time}; - const auto t_pre{t - cctk_delta_time}; - const auto A{amplitude}; - const auto W{gaussian_width}; - - phi(p.I) = gaussian_phi(A, W, t, p.x, p.y, p.z); - Pi(p.I) = gaussian_Pi(A, W, t, p.x, p.y, p.z); - Dx(p.I) = gaussian_Dx(A, W, t, p.x, p.y, p.z); - Dy(p.I) = gaussian_Dy(A, W, t, p.x, p.y, p.z); - Dz(p.I) = gaussian_Dz(A, W, t, p.x, p.y, p.z); - - phi_pre(p.I) = gaussian_phi(A, W, t_pre, p.x, p.y, p.z); - Pi_pre(p.I) = gaussian_Pi(A, W, t_pre, p.x, p.y, p.z); - Dx_pre(p.I) = gaussian_Dx(A, W, t_pre, p.x, p.y, p.z); - Dy_pre(p.I) = gaussian_Dy(A, W, t_pre, p.x, p.y, p.z); - Dz_pre(p.I) = gaussian_Dz(A, W, t_pre, p.x, p.y, p.z); + 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, "time")) { - 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 t_pre{t - cctk_delta_time}; - - phi(p.I) = t; - Pi(p.I) = t; - Dx(p.I) = t; - Dy(p.I) = t; - Dz(p.I) = t; - - phi_pre(p.I) = t_pre; - Pi_pre(p.I) = t_pre; - Dx_pre(p.I) = t_pre; - Dy_pre(p.I) = t_pre; - Dz_pre(p.I) = t_pre; - }); } else { - CCTK_ERROR("Unknown initial condition"); + CCTK_VERROR("Unknown initial condition \"%s\"", initial_condition); } } @@ -276,28 +73,16 @@ extern "C" void TestRKAB_RHS(CCTK_ARGUMENTS) { DECLARE_CCTK_ARGUMENTSX_TestRKAB_RHS; DECLARE_CCTK_PARAMETERS; - if (CCTK_EQUALS(initial_condition, "time")) { - grid.loop_int_device<0, 0, 0>(grid.nghostzones, - [=] CCTK_DEVICE(const Loop::PointDesc &p) - CCTK_ATTRIBUTE_ALWAYS_INLINE { - phi_rhs(p.I) = 1.0; - Pi_rhs(p.I) = 1.0; - Dx_rhs(p.I) = 1.0; - Dy_rhs(p.I) = 1.0; - Dz_rhs(p.I) = 1.0; - }); - } else { - 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); - }); - } + 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) { @@ -321,11 +106,11 @@ extern "C" void TestRKAB_Error(CCTK_ARGUMENTS) { 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 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)}; @@ -348,14 +133,16 @@ extern "C" void TestRKAB_Error(CCTK_ARGUMENTS) { const auto t{cctk_time}; - const auto A{amplitude}; - const auto W{gaussian_width}; - - const auto expected_phi{gaussian_phi(A, W, t, p.x, p.y, p.z)}; - const auto expected_Pi{gaussian_Pi(A, W, t, p.x, p.y, p.z)}; - const auto expected_Dx{gaussian_Dx(A, W, t, p.x, p.y, p.z)}; - const auto expected_Dy{gaussian_Dy(A, W, t, p.x, p.y, p.z)}; - const auto expected_Dz{gaussian_Dz(A, W, t, p.x, p.y, p.z)}; + 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)}; @@ -369,30 +156,8 @@ extern "C" void TestRKAB_Error(CCTK_ARGUMENTS) { Dy_err(p.I) = fabs(expected_Dy - actual_Dy); Dz_err(p.I) = fabs(expected_Dz - actual_Dz); }); - - } else if (CCTK_EQUALS(initial_condition, "time")) { - 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 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(t - actual_phi); - Pi_err(p.I) = fabs(t - actual_Pi); - Dx_err(p.I) = fabs(t - actual_Dx); - Dy_err(p.I) = fabs(t - actual_Dy); - Dz_err(p.I) = fabs(t - actual_Dz); - }); - } else { - CCTK_ERROR("Unknown initial condition"); + CCTK_VERROR("Unknown initial condition \"%s\"", initial_condition); } }