Skip to content

Commit

Permalink
TestRKAB: Fixed exact gaussian functions
Browse files Browse the repository at this point in the history
  • Loading branch information
@lucass-carneiro
lucass-carneiro committed Oct 8, 2024
1 parent 26a5930 commit e154cf6
Show file tree
Hide file tree
Showing 2 changed files with 42 additions and 43 deletions.
10 changes: 5 additions & 5 deletions TestRKAB/par/gaussian.par
View file
Original file line number Diff line number Diff line change
Expand Up @@ -8,7 +8,7 @@ ActiveThorns = "
$ncells = 80
$cfl = 0.25

$itlast = 1
$itlast = 50
$out_every = 1

TestRKAB::initial_condition = "Gaussian"
Expand Down Expand Up @@ -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
Expand Down
75 changes: 37 additions & 38 deletions TestRKAB/src/testrkab.cxx
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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));
}
}
Expand All @@ -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)));
}
}

Expand All @@ -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));
}
}

Expand All @@ -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));
}
}

Expand All @@ -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));
}
}

Expand Down

0 comments on commit e154cf6

Please sign in to comment.