Skip to content

Commit

Permalink
TestRKAB: Updated initial data functions and paramete files
Browse files Browse the repository at this point in the history
  • Loading branch information
Lucas Sanches committed Oct 22, 2024
1 parent afab507 commit 8679164
Show file tree
Hide file tree
Showing 7 changed files with 287 additions and 378 deletions.
46 changes: 29 additions & 17 deletions TestRKAB/par/gaussian.par
View file
Original file line number Diff line number Diff line change
Expand Up @@ -5,28 +5,42 @@ 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

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
Expand All @@ -37,30 +51,28 @@ 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

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
"
56 changes: 36 additions & 20 deletions TestRKAB/par/standing.par
View file
Original file line number Diff line number Diff line change
Expand Up @@ -5,31 +5,48 @@ 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

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
Expand All @@ -40,29 +57,28 @@ 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

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
"
"
64 changes: 0 additions & 64 deletions TestRKAB/par/time.par
View file

This file was deleted.

1 change: 0 additions & 1 deletion TestRKAB/param.ccl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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"
Expand Down
125 changes: 125 additions & 0 deletions TestRKAB/src/gaussian.hxx
View file
Original file line number Diff line number Diff line change
@@ -0,0 +1,125 @@
#ifndef TEST_RKAB_GAUSSIAN_HXX
#define TEST_RKAB_GAUSSIAN_HXX

#include <cmath>

namespace TestRKAB::gauss {

template <typename T>
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<T>::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 <typename T>
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<T>::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 <typename T>
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<T>::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 <typename T>
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<T>::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 <typename T>
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<T>::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
Loading

0 comments on commit 8679164

Please sign in to comment.