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

Jump to bottom

Denote pseudo time for Euler-Gravity #2101

Open
wants to merge 11 commits into
base: main
Choose a base branch
from
Open

Denote pseudo time for Euler-Gravity #2101

Changes from 1 commit
Commits
Show all changes
11 commits
Select commit Hold shift + click to select a range
518bef5
Denote pseudo time for Euler-Gravity
DanielDoehring Oct 8, 2024
166d9e5
Update src/semidiscretization/semidiscretization_euler_gravity.jl
DanielDoehring Oct 8, 2024
b277d7b
Apply suggestions from code review
DanielDoehring Oct 8, 2024
ff3d09e
further comments
DanielDoehring Oct 9, 2024
9f7d482
Merge branch 'Denote_pseudo_time_for_Euler-Gravity' of github.com:Dan…
DanielDoehring Oct 9, 2024
6c83c68
comment
DanielDoehring Oct 10, 2024
11548ec
dtau
DanielDoehring Oct 10, 2024
7283633
Merge branch 'main' into Denote_pseudo_time_for_Euler-Gravity
DanielDoehring Oct 10, 2024
eea03eb
Merge branch 'main' into Denote_pseudo_time_for_Euler-Gravity
DanielDoehring Oct 14, 2024
e24cf90
Merge branch 'main' into Denote_pseudo_time_for_Euler-Gravity
ranocha Oct 16, 2024
964aca5
Apply suggestions from code review
DanielDoehring Oct 16, 2024
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
60 changes: 35 additions & 25 deletions src/semidiscretization/semidiscretization_euler_gravity.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -19,8 +19,8 @@ struct ParametersEulerGravity{RealT <: Real, TimestepGravity}
background_density :: RealT # aka rho0
gravitational_constant :: RealT # aka G
cfl :: RealT
resid_tol :: RealT
n_iterations_max :: Int
resid_tol :: RealT # Hyp.-Diff. Eq. steady state tolerance
n_iterations_max :: Int # Max. number of its of the pseudo-time gravity solver
DanielDoehring marked this conversation as resolved.
Show resolved Hide resolved
timestep_gravity :: TimestepGravity
end

Expand Down Expand Up @@ -217,6 +217,8 @@ end
calc_error_norms(func, u, t, analyzer, semi.semi_euler, cache_analysis)
end

# Coupled Euler and gravity solver at each Runge-Kutta stage,
# corresponding to Algorithm 2 in Schlottke-Lakemper et al. (2020).
DanielDoehring marked this conversation as resolved.
Show resolved Hide resolved
function rhs!(du_ode, u_ode, semi::SemidiscretizationEulerGravity, t)
@unpack semi_euler, semi_gravity, cache = semi

Expand Down Expand Up @@ -275,33 +277,33 @@ function update_gravity!(semi::SemidiscretizationEulerGravity, u_ode)
finalstep = false
@unpack n_iterations_max, cfl, resid_tol, timestep_gravity = parameters
iter = 0
t = zero(real(semi_gravity.solver))
tau = zero(real(semi_gravity.solver)) # Pseudo-time

# iterate gravity solver until convergence or maximum number of iterations are reached
@unpack equations = semi_gravity
while !finalstep
dt = @trixi_timeit timer() "calculate dt" begin
cfl * max_dt(u_gravity, t, semi_gravity.mesh,
dtau = @trixi_timeit timer() "calculate dtau" begin
cfl * max_dt(u_gravity, tau, semi_gravity.mesh,
have_constant_speed(equations), equations,
semi_gravity.solver, semi_gravity.cache)
end

# evolve solution by one pseudo-time step
time_start = time_ns()
timestep_gravity(cache, u_euler, t, dt, parameters, semi_gravity)
timestep_gravity(cache, u_euler, tau, dtau, parameters, semi_gravity)
runtime = time_ns() - time_start
put!(gravity_counter, runtime)

# update iteration counter
iter += 1
t += dt
tau += dtau

# check if we reached the maximum number of iterations
if n_iterations_max > 0 && iter >= n_iterations_max
@warn "Max iterations reached: Gravity solver failed to converge!" residual=maximum(abs,
@views du_gravity[1,
..,
:]) t=t dt=dt
:]) tau=tau dtau=dtau
finalstep = true
end

Expand All @@ -315,7 +317,8 @@ function update_gravity!(semi::SemidiscretizationEulerGravity, u_ode)
end

# Integrate gravity solver for 2N-type low-storage schemes
function timestep_gravity_2N!(cache, u_euler, t, dt, gravity_parameters, semi_gravity,
function timestep_gravity_2N!(cache, u_euler, tau, dtau, gravity_parameters,
semi_gravity,
a, b, c)
G = gravity_parameters.gravitational_constant
rho0 = gravity_parameters.background_density
Expand All @@ -325,20 +328,20 @@ function timestep_gravity_2N!(cache, u_euler, t, dt, gravity_parameters, semi_gr
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]
tau_stage = tau + dtau * c[stage]

# rhs! has the source term for the harmonic problem
# We don't need a `@trixi_timeit timer() "rhs!"` here since that's already
# included in the `rhs!` call.
rhs!(du_ode, u_ode, semi_gravity, t_stage)
rhs!(du_ode, u_ode, semi_gravity, tau_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
b_stage_dt = b[stage] * dtau
@trixi_timeit timer() "Runge-Kutta step" begin
@threaded for idx in eachindex(u_ode)
u_tmp1_ode[idx] = du_ode[idx] - u_tmp1_ode[idx] * a_stage
Expand All @@ -350,7 +353,7 @@ function timestep_gravity_2N!(cache, u_euler, t, dt, gravity_parameters, semi_gr
return nothing
end

function timestep_gravity_carpenter_kennedy_erk54_2N!(cache, u_euler, t, dt,
function timestep_gravity_carpenter_kennedy_erk54_2N!(cache, u_euler, tau, dtau,
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,
Expand All @@ -363,12 +366,12 @@ function timestep_gravity_carpenter_kennedy_erk54_2N!(cache, u_euler, t, dt,
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)
timestep_gravity_2N!(cache, u_euler, tau, dtau, 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,
function timestep_gravity_3Sstar!(cache, u_euler, tau, dtau, gravity_parameters,
semi_gravity,
gamma1, gamma2, gamma3, beta, delta, c)
G = gravity_parameters.gravitational_constant
Expand All @@ -380,12 +383,12 @@ function timestep_gravity_3Sstar!(cache, u_euler, t, dt, gravity_parameters,
u_tmp2_ode .= u_ode
du_gravity = wrap_array(du_ode, semi_gravity)
for stage in eachindex(c)
t_stage = t + dt * c[stage]
tau_stage = tau + dtau * c[stage]

# rhs! has the source term for the harmonic problem
# We don't need a `@trixi_timeit timer() "rhs!"` here since that's already
# included in the `rhs!` call.
rhs!(du_ode, u_ode, semi_gravity, t_stage)
rhs!(du_ode, u_ode, semi_gravity, tau_stage)

# Source term: Jeans instability OR coupling convergence test OR blast wave
# put in gravity source term proportional to Euler density
Expand All @@ -396,9 +399,10 @@ function timestep_gravity_3Sstar!(cache, u_euler, t, dt, gravity_parameters,
gamma1_stage = gamma1[stage]
gamma2_stage = gamma2[stage]
gamma3_stage = gamma3[stage]
beta_stage_dt = beta[stage] * dt
beta_stage_dt = beta[stage] * dtau
@trixi_timeit timer() "Runge-Kutta step" begin
@threaded for idx in eachindex(u_ode)
# See Algorithm 1 (3S* method) from Schlottke-Lakemper et al. (2020)
DanielDoehring marked this conversation as resolved.
Show resolved Hide resolved
u_tmp1_ode[idx] += delta_stage * u_ode[idx]
u_ode[idx] = (gamma1_stage * u_ode[idx] +
gamma2_stage * u_tmp1_ode[idx] +
Expand All @@ -411,7 +415,8 @@ function timestep_gravity_3Sstar!(cache, u_euler, t, dt, gravity_parameters,
return nothing
end

function timestep_gravity_erk51_3Sstar!(cache, u_euler, t, dt, gravity_parameters,
# First-order 5-stage, 3S*-storage optimized method
DanielDoehring marked this conversation as resolved.
Show resolved Hide resolved
function timestep_gravity_erk51_3Sstar!(cache, u_euler, tau, dtau, gravity_parameters,
semi_gravity)
# New 3Sstar coefficients optimized for polynomials of degree polydeg=3
# and examples/parameters_hypdiff_lax_friedrichs.toml
Expand All @@ -434,11 +439,13 @@ function timestep_gravity_erk51_3Sstar!(cache, u_euler, t, dt, gravity_parameter
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,
timestep_gravity_3Sstar!(cache, u_euler, tau, dtau, gravity_parameters,
semi_gravity,
gamma1, gamma2, gamma3, beta, delta, c)
end

function timestep_gravity_erk52_3Sstar!(cache, u_euler, t, dt, gravity_parameters,
# Second-order 5-stage, 3S*-storage optimized method
DanielDoehring marked this conversation as resolved.
Show resolved Hide resolved
function timestep_gravity_erk52_3Sstar!(cache, u_euler, tau, dtau, gravity_parameters,
semi_gravity)
# New 3Sstar coefficients optimized for polynomials of degree polydeg=3
# and examples/parameters_hypdiff_lax_friedrichs.toml
Expand All @@ -461,11 +468,13 @@ function timestep_gravity_erk52_3Sstar!(cache, u_euler, t, dt, gravity_parameter
c = SVector(0.0000000000000000E+00, 4.5158640252832094E-01, 1.0221535725056414E+00,
1.4280257701954349E+00, 7.1581334196229851E-01)

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

function timestep_gravity_erk53_3Sstar!(cache, u_euler, t, dt, gravity_parameters,
# Third-order 5-stage, 3S*-storage optimized method
DanielDoehring marked this conversation as resolved.
Show resolved Hide resolved
function timestep_gravity_erk53_3Sstar!(cache, u_euler, tau, dtau, gravity_parameters,
semi_gravity)
# New 3Sstar coefficients optimized for polynomials of degree polydeg=3
# and examples/parameters_hypdiff_lax_friedrichs.toml
Expand All @@ -488,7 +497,8 @@ function timestep_gravity_erk53_3Sstar!(cache, u_euler, t, dt, gravity_parameter
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,
timestep_gravity_3Sstar!(cache, u_euler, tau, dtau, gravity_parameters,
semi_gravity,
gamma1, gamma2, gamma3, beta, delta, c)
end

Expand Down
Loading