From 26a5930abb43a3e2ed6178b07e8f542cdd399308 Mon Sep 17 00:00:00 2001 From: Lucas Date: Tue, 8 Oct 2024 16:28:10 -0500 Subject: [PATCH] TestRKAB: Fixed bug in RHS of evolution code --- TestRKAB/interface.ccl | 2 - TestRKAB/par/gaussian.par | 62 ++++------ TestRKAB/par/gaussian_single_mesh.par | 63 ---------- TestRKAB/par/standing.par | 11 +- TestRKAB/schedule.ccl | 10 +- TestRKAB/src/testrkab.cxx | 170 +++++++++++++++++++++++++- 6 files changed, 200 insertions(+), 118 deletions(-) delete mode 100644 TestRKAB/par/gaussian_single_mesh.par diff --git a/TestRKAB/interface.ccl b/TestRKAB/interface.ccl index ad192ed66..e0def2407 100644 --- a/TestRKAB/interface.ccl +++ b/TestRKAB/interface.ccl @@ -2,8 +2,6 @@ IMPLEMENTS: TestRKAB -INHERITS: CarpetX - USES INCLUDE HEADER: loop_device.hxx USES INCLUDE HEADER: vect.hxx diff --git a/TestRKAB/par/gaussian.par b/TestRKAB/par/gaussian.par index 8886f0dd4..8045d2663 100644 --- a/TestRKAB/par/gaussian.par +++ b/TestRKAB/par/gaussian.par @@ -1,14 +1,15 @@ ActiveThorns = " - BoxInBox CarpetX IOUtil ODESolvers TestRKAB " -$out_every = 16 -$nlevels = 3 -$ncells = 64 +$ncells = 80 +$cfl = 0.25 + +$itlast = 1 +$out_every = 1 TestRKAB::initial_condition = "Gaussian" TestRKAB::gaussian_width = 0.17677669529 # sqrt(2)*W = 0.25 @@ -19,66 +20,45 @@ CarpetX::verbose = no Cactus::cctk_show_schedule = yes Cactus::presync_mode = "mixed-error" -CarpetX::xmin = -4.0 -CarpetX::ymin = -4.0 -CarpetX::zmin = -4.0 +CarpetX::xmin = -1.0 +CarpetX::ymin = -1.0 +CarpetX::zmin = -1.0 -CarpetX::xmax = +4.0 -CarpetX::ymax = +4.0 -CarpetX::zmax = +4.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::max_num_levels = $nlevels -CarpetX::regrid_every = 128 -CarpetX::regrid_error_threshold = 0.9 - -CarpetX::prolongation_type = "ddf" -CarpetX::prolongation_order = 5 -CarpetX::ghost_size = 3 -CarpetX::dtfac = 1.0 # 0.25 * 2^2 +#CarpetX::ghost_size = 3 +CarpetX::dtfac = $cfl CarpetX::blocking_factor_x = 2 CarpetX::blocking_factor_y = 2 CarpetX::blocking_factor_z = 2 -#CarpetX::max_grid_size_x = 10000000 -#CarpetX::max_grid_size_y = 10000000 -#CarpetX::max_grid_size_z = 10000000 -#CarpetX::max_tile_size_x = 10000000 -#CarpetX::max_tile_size_y = 10000000 -#CarpetX::max_tile_size_z = 10000000 - -BoxInBox::num_regions = 1 -# Region 1 -BoxInBox::shape_1 = "cube" -BoxInBox::num_levels_1 = $nlevels -BoxInBox::radius_1 = [ -1.0, 1.0, 0.5 ] - Cactus::terminate = "time" -Cactus::cctk_final_time = 2.0 +Cactus::cctk_final_time = 1.0 +#Cactus::terminate = "iteration" #Cactus::cctk_itlast = $itlast -ODESolvers::method = "RK4" -#ODESolvers::verbose = yes +ODESolvers::method = "RKAB" +ODESolvers::verbose = yes IO::out_dir = $parfile IO::out_every = $out_every -#CarpetX::out_norm_vars = "all" +CarpetX::out_norm_vars = "" +CarpetX::out_tsv_vars = "" -CarpetX::out_tsv_vars = " +CarpetX::out_silo_vars = " TestRKAB::state TestRKAB::error " -# -#CarpetX::out_silo_vars = " -# TestRKAB::state -# TestRKAB::error -#" diff --git a/TestRKAB/par/gaussian_single_mesh.par b/TestRKAB/par/gaussian_single_mesh.par deleted file mode 100644 index 641bd687b..000000000 --- a/TestRKAB/par/gaussian_single_mesh.par +++ /dev/null @@ -1,63 +0,0 @@ -ActiveThorns = " - CarpetX - IOUtil - ODESolvers - TestRKAB -" - -$ncells = 80 -$cfl = 0.25 - -$itlast = 1 -$out_every = 1 - -TestRKAB::initial_condition = "Gaussian" -TestRKAB::gaussian_width = 0.17677669529 # sqrt(2)*W = 0.25 - -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_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 = 4.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_norm_vars = "" -CarpetX::out_tsv_vars = "" - -CarpetX::out_silo_vars = " - TestRKAB::state - TestRKAB::error -" diff --git a/TestRKAB/par/standing.par b/TestRKAB/par/standing.par index 4ca82e800..92d66c298 100644 --- a/TestRKAB/par/standing.par +++ b/TestRKAB/par/standing.par @@ -8,7 +8,7 @@ ActiveThorns = " $ncells = 80 $cfl = 0.25 -$itlast = 1 +$itlast = 50 $out_every = 1 TestRKAB::initial_condition = "standing wave" @@ -35,6 +35,7 @@ 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 @@ -46,10 +47,10 @@ CarpetX::blocking_factor_x = 2 CarpetX::blocking_factor_y = 2 CarpetX::blocking_factor_z = 2 -Cactus::terminate = "time" -Cactus::cctk_final_time = 4.0 -#Cactus::terminate = "iteration" -#Cactus::cctk_itlast = $itlast +#Cactus::terminate = "time" +#Cactus::cctk_final_time = 1.0 +Cactus::terminate = "iteration" +Cactus::cctk_itlast = $itlast ODESolvers::method = "RKAB" ODESolvers::verbose = yes diff --git a/TestRKAB/schedule.ccl b/TestRKAB/schedule.ccl index fc6fb3f35..7e6e6cc33 100644 --- a/TestRKAB/schedule.ccl +++ b/TestRKAB/schedule.ccl @@ -1,6 +1,6 @@ # Schedule definitions for thorn TestRKAB -STORAGE: state pre rhs +STORAGE: state pre rhs energy error # Init SCHEDULE TestRKAB_Initial AT initial @@ -10,12 +10,14 @@ SCHEDULE TestRKAB_Initial AT initial SYNC: state pre } "Initialize scalar wave state and previous state" + + # Step SCHEDULE TestRKAB_Sync AT postregrid { LANG: C OPTIONS: global - SYNC: state + SYNC: state pre } "Synchronize state" SCHEDULE TestRKAB_RHS IN ODESolvers_RHS @@ -29,9 +31,11 @@ SCHEDULE TestRKAB_Sync IN ODESolvers_PostStep { LANG: C OPTIONS: global - SYNC: state + SYNC: state pre } "Synchronize state" + + # Analysis SCHEDULE TestRKAB_Error IN ODESolvers_PostStep AFTER TestRKAB_Sync { diff --git a/TestRKAB/src/testrkab.cxx b/TestRKAB/src/testrkab.cxx index 791d8b75a..4100e6b74 100644 --- a/TestRKAB/src/testrkab.cxx +++ b/TestRKAB/src/testrkab.cxx @@ -70,6 +70,119 @@ static inline auto sw_Dz(T A, T kx, T ky, T kz, T t, T x, T y, 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 * exp(pow(t, 2) / (2. * pow(W, 2))) * t) / pow(W, 2); + } else { + return (A * (exp(pow(t - sqrt(pow(x, 2) + pow(y, 2) + pow(z, 2)), 2) / + (2. * pow(W, 2))) - + exp(pow(t + sqrt(pow(x, 2) + pow(y, 2) + pow(z, 2)), 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 * exp(pow(t, 2) / (2. * pow(W, 2))) * + (pow(t, 2) + pow(W, 2))) / + pow(W, 4); + } else { + return (-2 * A * + exp((pow(t, 2) + pow(x, 2) + pow(y, 2) + pow(z, 2)) / + (2. * pow(W, 2))) * + (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)))) / + (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 (-2 * A * + exp((pow(t, 2) + pow(x, 2) + pow(y, 2) + pow(z, 2)) / + (2. * pow(W, 2))) * + x * + (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)) + + (-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)))) / + (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 (-2 * A * + exp((pow(t, 2) + pow(x, 2) + pow(y, 2) + pow(z, 2)) / + (2. + pow(y, 2) + pow(z, 2)) * + cosh((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)) * + sinh((t * sqrt(pow(x, 2) + pow(y, 2) + pow(z, 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 (-2 * A * + exp((pow(t, 2) + pow(x, 2) + pow(y, 2) + pow(z, 2)) / + (2. * pow(W, 2))) * + z * + (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)) + + (-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)))) / + (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 }; @@ -83,6 +196,7 @@ static inline auto fd_c_1_4(const Loop::PointDesc &p, return den * num; } +// Scheduled functions extern "C" void TestRKAB_Initial(CCTK_ARGUMENTS) { DECLARE_CCTK_ARGUMENTSX_TestRKAB_Initial; DECLARE_CCTK_PARAMETERS; @@ -113,7 +227,28 @@ extern "C" void TestRKAB_Initial(CCTK_ARGUMENTS) { }); } else if (CCTK_EQUALS(initial_condition, "Gaussian")) { - CCTK_ERROR("Gaussian initial data not implemented yet"); + 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}; + + 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); + }); + } else { CCTK_ERROR("Unknown initial condition"); } @@ -127,7 +262,7 @@ extern "C" void TestRKAB_RHS(CCTK_ARGUMENTS) { 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, Dx) + + 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); @@ -144,7 +279,6 @@ extern "C" void TestRKAB_Error(CCTK_ARGUMENTS) { 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 { @@ -177,7 +311,35 @@ extern "C" void TestRKAB_Error(CCTK_ARGUMENTS) { }); } else if (CCTK_EQUALS(initial_condition, "Gaussian")) { - CCTK_VERROR("Gaussian error not implemented yet"); + 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 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 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_ERROR("Unknown initial condition"); }