Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

port 2N gravity time step to Taal #263

Merged
merged 1 commit into from
Oct 29, 2020
Merged

port 2N gravity time step to Taal #263

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
6 changes: 3 additions & 3 deletions examples/2d/elixir_euler_gravity_jeans_instability.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -39,9 +39,9 @@ semi_gravity = SemidiscretizationHyperbolic(mesh, equations_gravity, initial_con
# combining both semidiscretizations for Euler + self-gravity
parameters = ParametersEulerGravity(background_density=1.5e7, # aka rho0
gravitational_constant=6.674e-8, # aka G
cfl=2.4,
cfl=1.6,
n_iterations_max=1000,
timestep_gravity=timestep_gravity_erk52_3Sstar!)
timestep_gravity=timestep_gravity_carpenter_kennedy_erk54_2N!)

semi = SemidiscretizationEulerGravity(semi_euler, semi_gravity, parameters)

Expand All @@ -60,7 +60,7 @@ save_solution = SaveSolutionCallback(interval=10,
save_final_solution=true,
solution_variables=:primitive)
# TODO: Taal, IO
# restart_interval = 10
# restart_interval = 100

analysis_interval = 100
alive_callback = AliveCallback(analysis_interval=analysis_interval)
Expand Down
3 changes: 2 additions & 1 deletion src/Trixi.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -97,7 +97,8 @@ export DG,

export SemidiscretizationHyperbolic, semidiscretize, compute_coefficients, integrate

export SemidiscretizationEulerGravity, ParametersEulerGravity, timestep_gravity_erk52_3Sstar!
export SemidiscretizationEulerGravity, ParametersEulerGravity,
timestep_gravity_erk52_3Sstar!, timestep_gravity_carpenter_kennedy_erk54_2N!

export SummaryCallback, SteadyStateCallback, AMRCallback, StepsizeCallback,
AnalysisCallback, SaveSolutionCallback, AliveCallback
Expand Down
84 changes: 84 additions & 0 deletions src/semidiscretization_euler_gravity.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -242,6 +242,55 @@ function update_gravity!(semi::SemidiscretizationEulerGravity, u_ode::AbstractVe
end


# Integrate gravity solver for 2N-type low-storage schemes
function timestep_gravity_2N!(cache, u_euler, t, dt, gravity_parameters, semi_gravity,
a, b, c)
G = gravity_parameters.gravitational_constant
rho0 = gravity_parameters.background_density
grav_scale = -4.0*pi*G

@unpack u_ode, du_ode, u_tmp1_ode = cache
u_tmp1_ode .= zero(eltype(u_tmp1_ode))
du_gravity = wrap_array(du_ode, semi_gravity)
for stage in eachindex(c)
t_stage = t + dt * c[stage]

# rhs! has the source term for the harmonic problem
@timeit_debug timer() "rhs!" rhs!(du_ode, u_ode, semi_gravity, t_stage)

# Source term: Jeans instability OR coupling convergence test OR blast wave
# put in gravity source term proportional to Euler density
# OBS! subtract off the background density ρ_0 (spatial mean value)
@views @. du_gravity[1, .., :] += grav_scale * (u_euler[1, .., :] - rho0)

a_stage = a[stage]
b_stage_dt = b[stage] * dt
@timeit_debug timer() "Runge-Kutta step" begin
Threads.@threads for idx in eachindex(u_ode)
u_tmp1_ode[idx] = du_ode[idx] - u_tmp1_ode[idx] * a_stage
u_ode[idx] += u_tmp1_ode[idx] * b_stage_dt
end
end
end

return nothing
end

function timestep_gravity_carpenter_kennedy_erk54_2N!(cache, u_euler, t, dt, gravity_parameters, semi_gravity)
# Coefficients for Carpenter's 5-stage 4th-order low-storage Runge-Kutta method
a = @SVector [0.0, 567301805773.0 / 1357537059087.0,2404267990393.0 / 2016746695238.0,
3550918686646.0 / 2091501179385.0, 1275806237668.0 / 842570457699.0]
b = @SVector [1432997174477.0 / 9575080441755.0, 5161836677717.0 / 13612068292357.0,
1720146321549.0 / 2090206949498.0, 3134564353537.0 / 4481467310338.0,
2277821191437.0 / 14882151754819.0]
c = @SVector [0.0, 1432997174477.0 / 9575080441755.0, 2526269341429.0 / 6820363962896.0,
2006345519317.0 / 3224310063776.0, 2802321613138.0 / 2924317926251.0]

timestep_gravity_2N!(cache, u_euler, t, dt, gravity_parameters, semi_gravity, a, b, c)
end


# Integrate gravity solver for 3S*-type low-storage schemes
function timestep_gravity_3Sstar!(cache, u_euler, t, dt, gravity_parameters, semi_gravity,
gamma1, gamma2, gamma3, beta, delta, c)
G = gravity_parameters.gravitational_constant
Expand All @@ -254,6 +303,8 @@ function timestep_gravity_3Sstar!(cache, u_euler, t, dt, gravity_parameters, sem
du_gravity = wrap_array(du_ode, semi_gravity)
for stage in eachindex(c)
t_stage = t + dt * c[stage]

# rhs! has the source term for the harmonic problem
@timeit_debug timer() "rhs!" rhs!(du_ode, u_ode, semi_gravity, t_stage)

# Source term: Jeans instability OR coupling convergence test OR blast wave
Expand All @@ -280,6 +331,22 @@ function timestep_gravity_3Sstar!(cache, u_euler, t, dt, gravity_parameters, sem
return nothing
end

# This function is commented out since it's less performant than timestep_gravity_erk52_3Sstar!
# Nevertheless, we want to keep the coefficients as reference values.
# function timestep_gravity_erk51_3Sstar!(cache, u_euler, t, dt, gravity_parameters, semi_gravity)
# # New 3Sstar coefficients optimized for polynomials of degree polydeg=3
# # and examples/parameters_hyp_diff_llf.toml
# # 5 stages, order 1
# gamma1 = @SVector [0.0000000000000000E+00, 5.2910412316555866E-01, 2.8433964362349406E-01, -1.4467571130907027E+00, 7.5592215948661057E-02]
# gamma2 = @SVector [1.0000000000000000E+00, 2.6366970460864109E-01, 3.7423646095836322E-01, 7.8786901832431289E-01, 3.7754129043053775E-01]
# gamma3 = @SVector [0.0000000000000000E+00, 0.0000000000000000E+00, 0.0000000000000000E+00, 8.0043329115077388E-01, 1.3550099149374278E-01]
# beta = @SVector [1.9189497208340553E-01, 5.4506406707700059E-02, 1.2103893164085415E-01, 6.8582252490550921E-01, 8.7914657211972225E-01]
# delta = @SVector [1.0000000000000000E+00, 7.8593091509463076E-01, 1.2639038717454840E-01, 1.7726945920209813E-01, 0.0000000000000000E+00]
# c = @SVector [0.0000000000000000E+00, 1.9189497208340553E-01, 1.9580448818599061E-01, 2.4241635859769023E-01, 5.0728347557552977E-01]

# timestep_gravity_3Sstar!(cache, u_euler, t, dt, gravity_parameters, semi_gravity,
# gamma1, gamma2, gamma3, beta, delta, c)
# end

function timestep_gravity_erk52_3Sstar!(cache, u_euler, t, dt, gravity_parameters, semi_gravity)
# New 3Sstar coefficients optimized for polynomials of degree polydeg=3
Expand All @@ -296,6 +363,23 @@ function timestep_gravity_erk52_3Sstar!(cache, u_euler, t, dt, gravity_parameter
gamma1, gamma2, gamma3, beta, delta, c)
end

# This function is commented out since it's less performant than timestep_gravity_erk52_3Sstar!
# Nevertheless, we want to keep the coefficients as reference values.
# function timestep_gravity_erk53_3Sstar!(cache, u_euler, t, dt, gravity_parameters, semi_gravity)
# # New 3Sstar coefficients optimized for polynomials of degree polydeg=3
# # and examples/parameters_hyp_diff_llf.toml
# # 5 stages, order 3
# gamma1 = @SVector [0.0000000000000000E+00, 6.9362208054011210E-01, 9.1364483229179472E-01, 1.3129305757628569E+00, -1.4615811339132949E+00]
# gamma2 = @SVector [1.0000000000000000E+00, 1.3224582239681788E+00, 2.4213162353103135E-01, -3.8532017293685838E-01, 1.5603355704723714E+00]
# gamma3 = @SVector [0.0000000000000000E+00, 0.0000000000000000E+00, 0.0000000000000000E+00, 3.8306787039991996E-01, -3.5683121201711010E-01]
# beta = @SVector [8.4476964977404881E-02, 3.0834660698015803E-01, 3.2131664733089232E-01, 2.8783574345390539E-01, 8.2199204703236073E-01]
# delta = @SVector [1.0000000000000000E+00, -7.6832695815481578E-01, 1.2497251501714818E-01, 1.4496404749796306E+00, 0.0000000000000000E+00]
# c = @SVector [0.0000000000000000E+00, 8.4476964977404881E-02, 2.8110631488732202E-01, 5.7093842145029405E-01, 7.2999896418559662E-01]

# timestep_gravity_3Sstar!(cache, u_euler, t, dt, gravity_parameters, semi_gravity,
# gamma1, gamma2, gamma3, beta, delta, c)
# end


# TODO: Taal decide, where should specific parts like these be?
@inline function save_solution_file(u_ode::AbstractVector, t, dt, iter,
Expand Down
6 changes: 3 additions & 3 deletions test/test_examples_2d.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -237,9 +237,9 @@ const EXAMPLES_DIR = joinpath(pathof(Trixi) |> dirname |> dirname, "examples", "

@testset "elixir_euler_gravity_jeans_instability.jl" begin
test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_gravity_jeans_instability.jl"),
l2 = [21174.129216837846, 978.8980277332366, 4.723889542064964e-6, 52935.30182880739],
linf = [29951.765533944592, 1388.65770524552, 2.0555889477917298e-5, 74879.49312048778],
tspan = (0.0, 0.6))
l2 = [10733.633835334986, 13356.780418347276, 1.6387671498796526e-6, 26834.07694603585],
linf = [15194.296497141942, 18881.481420845837, 8.91191153038996e-6, 37972.99718572572],
tspan = (0.0, 0.1))
end

@testset "elixir_euler_gravity_sedov_blast_wave_shockcapturing_amr.jl" begin
Expand Down