From e154cf67d92afc93ee3057dde52fcd25b74e37aa Mon Sep 17 00:00:00 2001 From: Lucas Date: Tue, 8 Oct 2024 17:14:46 -0500 Subject: [PATCH] TestRKAB: Fixed exact gaussian functions --- TestRKAB/par/gaussian.par | 10 +++--- TestRKAB/src/testrkab.cxx | 75 +++++++++++++++++++-------------------- 2 files changed, 42 insertions(+), 43 deletions(-) diff --git a/TestRKAB/par/gaussian.par b/TestRKAB/par/gaussian.par index 8045d2663..cf41229f7 100644 --- a/TestRKAB/par/gaussian.par +++ b/TestRKAB/par/gaussian.par @@ -8,7 +8,7 @@ ActiveThorns = " $ncells = 80 $cfl = 0.25 -$itlast = 1 +$itlast = 50 $out_every = 1 TestRKAB::initial_condition = "Gaussian" @@ -44,10 +44,10 @@ 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 +#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/src/testrkab.cxx b/TestRKAB/src/testrkab.cxx index 4100e6b74..b85f35486 100644 --- a/TestRKAB/src/testrkab.cxx +++ b/TestRKAB/src/testrkab.cxx @@ -79,12 +79,13 @@ static inline auto gaussian_phi(T A, T W, T t, T x, T y, T z) noexcept -> T { 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); + return (2 * A * t) / (exp(pow(t, 2) / (2. * pow(W, 2))) * 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))))) / + 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)); } } @@ -97,19 +98,18 @@ static inline auto gaussian_Pi(T A, T W, T t, T x, T y, T z) noexcept -> T { 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); + return (2 * A * (-t + W) * (t + W)) / + (exp(pow(t, 2) / (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))) * + 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)) + + 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))); + (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))); } } @@ -123,17 +123,16 @@ static inline auto gaussian_Dx(T A, T W, T t, T x, T y, T z) noexcept -> T { 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)) * + 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)) + - (-pow(W, 2) + 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)))) / - (pow(W, 2) * pow(pow(x, 2) + pow(y, 2) + pow(z, 2), 1.5)); + (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)); } } @@ -147,15 +146,16 @@ static inline auto gaussian_Dy(T A, T W, T t, T x, T y, T z) noexcept -> T { 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)); + 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)); } } @@ -169,17 +169,16 @@ static inline auto gaussian_Dz(T A, T W, T t, T x, T y, T z) noexcept -> T { 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)) * + 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)) + - (-pow(W, 2) + 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)))) / - (pow(W, 2) * pow(pow(x, 2) + pow(y, 2) + pow(z, 2), 1.5)); + (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)); } }