SSP IMEX and Linearized SSP IMEX #34
Open
Conversation
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
|
Your PR requires formatting changes to meet the project's style guidelines. Click here to view the suggested changes.diff --git a/examples/trixi_adv_diff_imex.jl b/examples/trixi_adv_diff_imex.jl
index 1fa5e95..595707e 100644
--- a/examples/trixi_adv_diff_imex.jl
+++ b/examples/trixi_adv_diff_imex.jl
@@ -34,14 +34,18 @@ coordinates_min = (-1.0, -1.0) # minimum coordinates (min(x), min(y))
coordinates_max = (1.0, 1.0) # maximum coordinates (max(x), max(y))
# Create a uniformly refined mesh with periodic boundaries
-mesh = TreeMesh(coordinates_min, coordinates_max,
- initial_refinement_level = 4,
- periodicity = true,
- n_cells_max = 30_000) # set maximum capacity of tree data structure
+mesh = TreeMesh(
+ coordinates_min, coordinates_max,
+ initial_refinement_level = 4,
+ periodicity = true,
+ n_cells_max = 30_000
+) # set maximum capacity of tree data structure
# Define initial condition
-function initial_condition_diffusive_convergence_test(x, t,
- equation::LinearScalarAdvectionEquation2D)
+function initial_condition_diffusive_convergence_test(
+ x, t,
+ equation::LinearScalarAdvectionEquation2D
+ )
# Store translated coordinate for easy use of exact solution
RealT = eltype(x)
x_trans = x - equation.advection_velocity * t
@@ -62,12 +66,16 @@ boundary_conditions = boundary_condition_periodic
boundary_conditions_parabolic = boundary_condition_periodic
# A semidiscretization collects data structures and functions for the spatial discretization
-semi = SemidiscretizationHyperbolicParabolic(mesh,
- (equations, equations_parabolic),
- initial_condition, solver;
- solver_parabolic = ViscousFormulationBassiRebay1(),
- boundary_conditions = (boundary_conditions,
- boundary_conditions_parabolic))
+semi = SemidiscretizationHyperbolicParabolic(
+ mesh,
+ (equations, equations_parabolic),
+ initial_condition, solver;
+ solver_parabolic = ViscousFormulationBassiRebay1(),
+ boundary_conditions = (
+ boundary_conditions,
+ boundary_conditions_parabolic,
+ )
+)
###############################################################################
# ODE solvers, callbacks etc.
@@ -94,10 +102,10 @@ callbacks = CallbackSet(summary_callback, analysis_callback, alive_callback)
# run the simulation
sol = solve(
- ode,
- Implicit.RKImplicitExplicitEuler();
+ ode,
+ Implicit.RKImplicitExplicitEuler();
dt = 0.001, # solve needs some value here but it will be overwritten by the stepsize_callback
- ode_default_options()..., callback = callbacks,
- # verbose=1,
- krylov_algo = :gmres,
+ ode_default_options()..., callback = callbacks,
+ # verbose=1,
+ krylov_algo = :gmres,
);
diff --git a/examples/trixi_fvm_lid_driven.jl b/examples/trixi_fvm_lid_driven.jl
index 2f9ba74..9840c90 100644
--- a/examples/trixi_fvm_lid_driven.jl
+++ b/examples/trixi_fvm_lid_driven.jl
@@ -8,8 +8,10 @@ prandtl_number() = 0.72
mu = 0.001
equations = CompressibleEulerEquations2D(1.4)
-equations_parabolic = CompressibleNavierStokesDiffusion2D(equations, mu = mu,
- Prandtl = prandtl_number())
+equations_parabolic = CompressibleNavierStokesDiffusion2D(
+ equations, mu = mu,
+ Prandtl = prandtl_number()
+)
# Create DG solver with polynomial degree = 3 and (local) Lax-Friedrichs/Rusanov flux as surface flux
solver = DGSEM(polydeg = 0, surface_flux = flux_lax_friedrichs)
@@ -19,10 +21,12 @@ coordinates_max = (1.0, 1.0) # maximum coordinates (max(x), max(y))
# Create a uniformly refined mesh
trees_per_dimension = (16, 16)
-mesh = P4estMesh(trees_per_dimension,
- polydeg = 0, initial_refinement_level = 0,
- coordinates_min = coordinates_min, coordinates_max = coordinates_max,
- periodicity = (false, false))
+mesh = P4estMesh(
+ trees_per_dimension,
+ polydeg = 0, initial_refinement_level = 0,
+ coordinates_min = coordinates_min, coordinates_max = coordinates_max,
+ periodicity = (false, false)
+)
function initial_condition_cavity(x, t, equations::CompressibleEulerEquations2D)
Ma = 0.1
@@ -40,21 +44,29 @@ heat_bc = Adiabatic((x, t, equations_parabolic) -> 0.0)
boundary_condition_lid = BoundaryConditionNavierStokesWall(velocity_bc_lid, heat_bc)
boundary_condition_cavity = BoundaryConditionNavierStokesWall(velocity_bc_cavity, heat_bc)
-boundary_conditions = Dict(:x_neg => boundary_condition_slip_wall,
- :y_neg => boundary_condition_slip_wall,
- :y_pos => boundary_condition_slip_wall,
- :x_pos => boundary_condition_slip_wall)
+boundary_conditions = Dict(
+ :x_neg => boundary_condition_slip_wall,
+ :y_neg => boundary_condition_slip_wall,
+ :y_pos => boundary_condition_slip_wall,
+ :x_pos => boundary_condition_slip_wall
+)
-boundary_conditions_parabolic = Dict(:x_neg => boundary_condition_cavity,
- :y_neg => boundary_condition_cavity,
- :y_pos => boundary_condition_lid,
- :x_pos => boundary_condition_cavity)
+boundary_conditions_parabolic = Dict(
+ :x_neg => boundary_condition_cavity,
+ :y_neg => boundary_condition_cavity,
+ :y_pos => boundary_condition_lid,
+ :x_pos => boundary_condition_cavity
+)
# A semidiscretization collects data structures and functions for the spatial discretization
-semi = SemidiscretizationHyperbolicParabolic(mesh, (equations, equations_parabolic),
- initial_condition, solver;
- boundary_conditions = (boundary_conditions,
- boundary_conditions_parabolic))
+semi = SemidiscretizationHyperbolicParabolic(
+ mesh, (equations, equations_parabolic),
+ initial_condition, solver;
+ boundary_conditions = (
+ boundary_conditions,
+ boundary_conditions_parabolic,
+ )
+)
###############################################################################
# ODE solvers, callbacks etc.
@@ -73,11 +85,11 @@ callbacks = CallbackSet(summary_callback, alive_callback)
# run the simulation
sol = solve(
- ode,
- Implicit.RKLinearImplicitExplicitEuler();
- #Implicit.KS22();
+ ode,
+ Implicit.RKLinearImplicitExplicitEuler();
+ #Implicit.KS22();
dt = 0.00001, # solve needs some value here but it will be overwritten by the stepsize_callback
- ode_default_options()..., callback = callbacks,
- # verbose=1,
- krylov_algo = :gmres,
+ ode_default_options()..., callback = callbacks,
+ # verbose=1,
+ krylov_algo = :gmres,
);
diff --git a/examples/trixi_imex.jl b/examples/trixi_imex.jl
index f39e5f1..e8e29f4 100644
--- a/examples/trixi_imex.jl
+++ b/examples/trixi_imex.jl
@@ -25,10 +25,10 @@ trixi_include(joinpath(examples_dir(), "tree_2d_dgsem", "elixir_advection_diffus
# run the simulation
sol = solve(
- ode,
- Implicit.RKImplicitExplicitEuler();
+ ode,
+ Implicit.RKImplicitExplicitEuler();
dt = 0.001, # solve needs some value here but it will be overwritten by the stepsize_callback
- ode_default_options()..., callback = callbacks,
- # verbose=1,
- krylov_algo = :gmres,
+ ode_default_options()..., callback = callbacks,
+ # verbose=1,
+ krylov_algo = :gmres,
);
diff --git a/examples/trixi_lid_driven_linear.jl b/examples/trixi_lid_driven_linear.jl
index 5327b5c..093f530 100644
--- a/examples/trixi_lid_driven_linear.jl
+++ b/examples/trixi_lid_driven_linear.jl
@@ -27,13 +27,13 @@ ode = semidiscretize(semi, (0.0, 10.0))
# run the simulation
sol = solve(
- ode,
- Implicit.RKLSSPIMEX332();
- #Implicit.KS22();
+ ode,
+ Implicit.RKLSSPIMEX332();
+ #Implicit.KS22();
dt = 0.001, # solve needs some value here but it will be overwritten by the stepsize_callback
- ode_default_options()..., callback = callbacks,
- # verbose=1,
- krylov_algo = :gmres,
+ ode_default_options()..., callback = callbacks,
+ # verbose=1,
+ krylov_algo = :gmres,
);
## dt_vec = [0.001, 0.001/2, 0.0001]
diff --git a/examples/trixi_linear_imex.jl b/examples/trixi_linear_imex.jl
index 9816efc..5e3f9c6 100644
--- a/examples/trixi_linear_imex.jl
+++ b/examples/trixi_linear_imex.jl
@@ -25,10 +25,10 @@ trixi_include(joinpath(examples_dir(), "tree_2d_dgsem", "elixir_advection_diffus
# run the simulation
sol = solve(
- ode,
- Implicit.RKLinearImplicitExplicitEuler();
+ ode,
+ Implicit.RKLinearImplicitExplicitEuler();
dt = 0.001, # solve needs some value here but it will be overwritten by the stepsize_callback
- ode_default_options()..., callback = callbacks,
- # verbose=1,
- krylov_algo = :gmres,
+ ode_default_options()..., callback = callbacks,
+ # verbose=1,
+ krylov_algo = :gmres,
);
diff --git a/examples/trixi_navier_stokes_imex.jl b/examples/trixi_navier_stokes_imex.jl
index f1c641b..a30cc8b 100644
--- a/examples/trixi_navier_stokes_imex.jl
+++ b/examples/trixi_navier_stokes_imex.jl
@@ -26,12 +26,12 @@ trixi_include(joinpath(examples_dir(), "tree_2d_dgsem", "elixir_navierstokes_con
# run the simulation
sol = solve(
- ode,
- Implicit.RKImplicitExplicitEuler();
- dt = 0.001/2, # solve needs some value here but it will be overwritten by the stepsize_callback
- ode_default_options()..., callback = callbacks,
- # verbose=1,
- krylov_algo = :gmres,
+ ode,
+ Implicit.RKImplicitExplicitEuler();
+ dt = 0.001 / 2, # solve needs some value here but it will be overwritten by the stepsize_callback
+ ode_default_options()..., callback = callbacks,
+ # verbose=1,
+ krylov_algo = :gmres,
);
## dt_vec = [0.001, 0.001/2, 0.0001]
diff --git a/examples/trixi_navier_stokes_linear_imex.jl b/examples/trixi_navier_stokes_linear_imex.jl
index c115991..0aa790a 100644
--- a/examples/trixi_navier_stokes_linear_imex.jl
+++ b/examples/trixi_navier_stokes_linear_imex.jl
@@ -25,11 +25,11 @@ trixi_include(joinpath(examples_dir(), "tree_2d_dgsem", "elixir_navierstokes_con
# run the simulation
sol = solve(
- ode,
- Implicit.RKLinearImplicitExplicitEuler();
- dt = 0.001/2, # solve needs some value here but it will be overwritten by the stepsize_callback
- ode_default_options()..., callback = callbacks,
- # verbose=1,
- krylov_algo = :gmres,
+ ode,
+ Implicit.RKLinearImplicitExplicitEuler();
+ dt = 0.001 / 2, # solve needs some value here but it will be overwritten by the stepsize_callback
+ ode_default_options()..., callback = callbacks,
+ # verbose=1,
+ krylov_algo = :gmres,
);
diff --git a/libs/Implicit/src/Implicit.jl b/libs/Implicit/src/Implicit.jl
index d2f000f..ab8fa34 100644
--- a/libs/Implicit/src/Implicit.jl
+++ b/libs/Implicit/src/Implicit.jl
@@ -13,8 +13,8 @@ struct MOperator{JOp}
end
struct LMOperator{JOp}
- J::JOp
-dt::Float64
+ J::JOp
+ dt::Float64
end
Base.size(M::LMOperator) = size(M.J)
@@ -625,10 +625,10 @@ end
function stage!(integrator, alg::Direct)
F!(du, u, p) = integrator.f(du, u, p, integrator.t)
- J = JacobianOperator(F!, integrator.du, integrator.u, integrator.p)
- M = MOperator(J, integrator.dt)
+ J = JacobianOperator(F!, integrator.du, integrator.u, integrator.p)
+ M = MOperator(J, integrator.dt)
kc = KrylovConstructor(integrator.res)
- workspace = krylov_workspace(:gmres, kc)
+ workspace = krylov_workspace(:gmres, kc)
for stage in 1:stages(alg)
alg(integrator.res, integrator.u, integrator.dt, integrator.f, integrator.du, integrator.u_tmp, integrator.p, integrator.t, integrator.stages, stage, workspace, M, integrator.RK)
diff --git a/libs/Implicit/src/imex.jl b/libs/Implicit/src/imex.jl
index 41f1f77..c260175 100644
--- a/libs/Implicit/src/imex.jl
+++ b/libs/Implicit/src/imex.jl
@@ -8,11 +8,11 @@ struct RKImplicitExplicitEuler <: RKIMEX{1} end
function (::RKImplicitExplicitEuler)(res, uₙ, Δt, f1!, f2!, du, du_tmp, u, p, t, stages, stage, RK)
if stage == 1
# Stage 1:
- ## f2 is the conservative part
- ## f1 is the parabolic part
- f2!(du, uₙ, p, t + RK.c[stage] * Δt )
- f1!(du_tmp, u, p, t + RK.c[stage] * Δt )
- return res .= u .- uₙ .- RK.a[stage, stage] * Δt .* du - RK.a[stage,stage] * Δt * du_tmp
+ ## f2 is the conservative part
+ ## f1 is the parabolic part
+ f2!(du, uₙ, p, t + RK.c[stage] * Δt)
+ f1!(du_tmp, u, p, t + RK.c[stage] * Δt)
+ return res .= u .- uₙ .- RK.a[stage, stage] * Δt .* du - RK.a[stage, stage] * Δt * du_tmp
else
@. u = uₙ + RK.b[1] * Δt * stages[1]
end
@@ -55,7 +55,7 @@ function ImplicitExplicitEulerTableau()
end
-function SimpleImplicitExplicitOptions(callback, tspan; maxiters=typemax(Int), verbose=0, krylov_algo=:gmres, krylov_kwargs=(;), kwargs...)
+function SimpleImplicitExplicitOptions(callback, tspan; maxiters = typemax(Int), verbose = 0, krylov_algo = :gmres, krylov_kwargs = (;), kwargs...)
return SimpleImplicitExplicitOptions{typeof(callback)}(
callback, false, Inf, maxiters,
[last(tspan)],
@@ -66,14 +66,14 @@ function SimpleImplicitExplicitOptions(callback, tspan; maxiters=typemax(Int), v
end
mutable struct SimpleImplicitExplicit{
- RealT<:Real,uType,Params,Sol,F,F1,F2,M,Alg<:SimpleImplicitExplicitAlgorithm,
- SimpleImplicitExplicitOptions,RKTableau,
-} <: AbstractTimeIntegrator
+ RealT <: Real, uType, Params, Sol, F, F1, F2, M, Alg <: SimpleImplicitExplicitAlgorithm,
+ SimpleImplicitExplicitOptions, RKTableau,
+ } <: AbstractTimeIntegrator
u::uType
du::uType
du_tmp::uType
u_tmp::uType
- stages::NTuple{M,uType}
+ stages::NTuple{M, uType}
res::uType
t::RealT
dt::RealT # current time step
@@ -93,16 +93,16 @@ end
function Base.getproperty(integrator::SimpleImplicitExplicit, field::Symbol)
if field === :stats
- return (naccept=getfield(integrator, :iter),)
+ return (naccept = getfield(integrator, :iter),)
end
# general fallback
return getfield(integrator, field)
end
function init(
- ode::ODEProblem, alg::SimpleImplicitExplicitAlgorithm{N};
- dt, callback::Union{CallbackSet,Nothing}=nothing, kwargs...,
-) where {N}
+ ode::ODEProblem, alg::SimpleImplicitExplicitAlgorithm{N};
+ dt, callback::Union{CallbackSet, Nothing} = nothing, kwargs...,
+ ) where {N}
u = copy(ode.u0)
du = zero(u)
res = zero(u)
@@ -111,12 +111,13 @@ function init(
t = first(ode.tspan)
iter = 0
integrator = SimpleImplicitExplicit(
- u, du, copy(du), u_tmp, stages, res, t, dt, zero(dt), iter, ode.p,
- (prob=ode,), ode.f.f1, ode.f.f1, ode.f.f2, alg,
+ u, du, copy(du), u_tmp, stages, res, t, dt, zero(dt), iter, ode.p,
+ (prob = ode,), ode.f.f1, ode.f.f1, ode.f.f2, alg,
SimpleImplicitExplicitOptions(
callback, ode.tspan;
kwargs...,
- ), false, RKTableau(alg))
+ ), false, RKTableau(alg)
+ )
# initialize callbacks
if callback isa CallbackSet
@@ -133,10 +134,10 @@ end
# Fakes `solve`: https://diffeq.sciml.ai/v6.8/basics/overview/#Solving-the-Problems-1
function solve(
- ode::ODEProblem, alg::SimpleImplicitExplicitAlgorithm;
- dt, callback=nothing, kwargs...,
-)
- integrator = init(ode, alg, dt=dt, callback=callback; kwargs...)
+ ode::ODEProblem, alg::SimpleImplicitExplicitAlgorithm;
+ dt, callback = nothing, kwargs...,
+ )
+ integrator = init(ode, alg, dt = dt, callback = callback; kwargs...)
# Start actual solve
return solve!(integrator)
@@ -174,7 +175,7 @@ function step!(integrator::SimpleImplicitExplicit)
# if the next iteration would push the simulation beyond the end time, set dt accordingly
if integrator.t + integrator.dt > t_end ||
- isapprox(integrator.t + integrator.dt, t_end)
+ isapprox(integrator.t + integrator.dt, t_end)
integrator.dt = t_end - integrator.t
terminate!(integrator)
end
@@ -214,22 +215,23 @@ function stage!(integrator, alg::RKIMEX)
F! = nonlinear_problem(alg, integrator.f2)
# TODO: Pass in `stages[1:(stage-1)]` or full tuple?
_, stats = Ariadne.newton_krylov!(
- F!, integrator.u_tmp, (integrator.u, integrator.dt, integrator.f1, integrator.du, integrator.du_tmp, integrator.p, integrator.t, integrator.stages, stage, integrator.RK), integrator.res;
- verbose=integrator.opts.verbose, krylov_kwargs=integrator.opts.krylov_kwargs,
- algo=integrator.opts.algo, tol_abs=6.0e-6,
+ F!, integrator.u_tmp, (integrator.u, integrator.dt, integrator.f1, integrator.du, integrator.du_tmp, integrator.p, integrator.t, integrator.stages, stage, integrator.RK), integrator.res;
+ verbose = integrator.opts.verbose, krylov_kwargs = integrator.opts.krylov_kwargs,
+ algo = integrator.opts.algo, tol_abs = 6.0e-6,
)
@assert stats.solved
# Store the solution for each stage in stages
- ## For a split Problem we need to compute rhs_conservative and rhs_parabolic
+ ## For a split Problem we need to compute rhs_conservative and rhs_parabolic
integrator.f2(integrator.du, integrator.u_tmp, integrator.p, integrator.t + integrator.RK.c[stage] * integrator.dt)
- integrator.stages[stage] .= integrator.du
+ integrator.stages[stage] .= integrator.du
integrator.f1(integrator.du, integrator.u_tmp, integrator.p, integrator.t + integrator.RK.c[stage] * integrator.dt)
integrator.stages[stage] .+= integrator.du
- if stage == stages(alg)
+ if stage == stages(alg)
alg(integrator.res, integrator.u, integrator.dt, integrator.f1, integrator.f2, integrator.du, integrator.du_tmp, integrator.u_tmp, integrator.p, integrator.t, integrator.stages, stage + 1, integrator.RK)
end
end
+ return
end
# get a cache where the RHS can be stored
diff --git a/libs/Implicit/src/linear_imex.jl b/libs/Implicit/src/linear_imex.jl
index e2aa144..e002dce 100644
--- a/libs/Implicit/src/linear_imex.jl
+++ b/libs/Implicit/src/linear_imex.jl
@@ -3,13 +3,13 @@ abstract type SimpleLinearImplicitExplicitAlgorithm{N} end
abstract type RKLIMEX{N} <: SimpleLinearImplicitExplicitAlgorithm{N} end
-struct IMEXRKButcher{T1<:AbstractArray,T2<:AbstractArray} <: RKTableau
+struct IMEXRKButcher{T1 <: AbstractArray, T2 <: AbstractArray} <: RKTableau
a::T1
b::T2
c::T2
ah::T1
bh::T2
- ch::T2
+ ch::T2
end
struct RKLinearImplicitExplicitEuler <: RKLIMEX{1} end
@@ -24,79 +24,79 @@ function mul!(out::AbstractVector, M::LMOperator, v::AbstractVector)
end
function (::RKLinearImplicitExplicitEuler)(res, uₙ, Δt, f1!, f2!, du, du_tmp, u, p, t, stages, ustages, jstages, stage, RK, M, lin_du_tmp, lin_du_tmp1, workspace)
- if stage == 1
+ return if stage == 1
# Stage 1:
- ## f2 is the conservative part
- ## f1 is the parabolic part
- mul!(lin_du_tmp, M.J, uₙ)
-# mul!(lin_du_tmp1, J, u)
- f2!(du, uₙ, p, t + RK.c[stage] * Δt)
+ ## f2 is the conservative part
+ ## f1 is the parabolic part
+ mul!(lin_du_tmp, M.J, uₙ)
+ # mul!(lin_du_tmp1, J, u)
+ f2!(du, uₙ, p, t + RK.c[stage] * Δt)
f1!(du_tmp, uₙ, p, t + RK.c[stage] * Δt)
-
- res .= uₙ .+ RK.a[stage, stage] * Δt .* (du .+ du_tmp .- lin_du_tmp)
- krylov_solve!(workspace, M, res, atol = 1e-6, rtol = 1e-6)
- @. u = workspace.x
-# f2!(du, workspace.x, p, t + RK.c[stage] * Δt)
-# f1!(du_tmp, workspace.x, p, t + RK.c[stage] * Δt)
-# stages[stage] .= du .+ du_tmp
-# @. u = uₙ + RK.b[1] * Δt * stages[1]
+
+ res .= uₙ .+ RK.a[stage, stage] * Δt .* (du .+ du_tmp .- lin_du_tmp)
+ krylov_solve!(workspace, M, res, atol = 1.0e-6, rtol = 1.0e-6)
+ @. u = workspace.x
+ # f2!(du, workspace.x, p, t + RK.c[stage] * Δt)
+ # f1!(du_tmp, workspace.x, p, t + RK.c[stage] * Δt)
+ # stages[stage] .= du .+ du_tmp
+ # @. u = uₙ + RK.b[1] * Δt * stages[1]
end
end
function (::RKLIMEX{3})(res, uₙ, Δt, f1!, f2!, du, du_tmp, u, p, t, stages, ustages, jstages, stage, RK, M, lin_du_tmp, lin_du_tmp1, workspace)
- F!(du, u, p) = f1!(du, u, p, t) ## parabolic
- if stage == 1
+ F!(du, u, p) = f1!(du, u, p, t) ## parabolic
+ return if stage == 1
# Stage 1:
- ## f2 is the conservative part
- ## f1 is the parabolic part
- J = JacobianOperator(F!, du, uₙ, p)
- M = LMOperator(J, RK.ah[stage,stage] * Δt)
- krylov_solve!(workspace, M, uₙ, atol = 1e-6, rtol = 1e-6)
- @. u = workspace.x
- J = JacobianOperator(F!, du, u, p)
- mul!(jstages[stage], J, u)
-
- f2!(du, workspace.x, p, t + RK.c[stage] * Δt)
+ ## f2 is the conservative part
+ ## f1 is the parabolic part
+ J = JacobianOperator(F!, du, uₙ, p)
+ M = LMOperator(J, RK.ah[stage, stage] * Δt)
+ krylov_solve!(workspace, M, uₙ, atol = 1.0e-6, rtol = 1.0e-6)
+ @. u = workspace.x
+ J = JacobianOperator(F!, du, u, p)
+ mul!(jstages[stage], J, u)
+
+ f2!(du, workspace.x, p, t + RK.c[stage] * Δt)
f1!(du_tmp, workspace.x, p, t + RK.c[stage] * Δt)
- stages[stage] .= du .+ du_tmp - jstages[stage]
- ustages[stage] .= u
-# @. u = uₙ + RK.b[1] * Δt * stages[1]
- elseif stage == 2
-
- J = JacobianOperator(F!, du, ustages[1], p)
- M = LMOperator(J, RK.ah[stage,stage] * Δt)
- mul!(lin_du_tmp, M.J, uₙ)
-# mul!(lin_du_tmp1, J, u)
- res .= uₙ + RK.a[stage,1] * Δt * stages[1] - Δt * RK.ah[stage,1] * jstages[1]
-
- krylov_solve!(workspace, M, res, atol = 1e-6, rtol = 1e-6)
- @. u = workspace.x
- J = JacobianOperator(F!, du, u, p)
- mul!(jstages[stage], J, u)
- f2!(du, workspace.x, p, t + RK.c[stage] * Δt)
+ stages[stage] .= du .+ du_tmp - jstages[stage]
+ ustages[stage] .= u
+ # @. u = uₙ + RK.b[1] * Δt * stages[1]
+ elseif stage == 2
+
+ J = JacobianOperator(F!, du, ustages[1], p)
+ M = LMOperator(J, RK.ah[stage, stage] * Δt)
+ mul!(lin_du_tmp, M.J, uₙ)
+ # mul!(lin_du_tmp1, J, u)
+ res .= uₙ + RK.a[stage, 1] * Δt * stages[1] - Δt * RK.ah[stage, 1] * jstages[1]
+
+ krylov_solve!(workspace, M, res, atol = 1.0e-6, rtol = 1.0e-6)
+ @. u = workspace.x
+ J = JacobianOperator(F!, du, u, p)
+ mul!(jstages[stage], J, u)
+ f2!(du, workspace.x, p, t + RK.c[stage] * Δt)
f1!(du_tmp, workspace.x, p, t + RK.c[stage] * Δt)
- stages[stage] .= du .+ du_tmp - jstages[stage]
- ustages[stage] .= u
-
- elseif stage == 3
-
- J = JacobianOperator(F!, du, ustages[2], p)
- M = LMOperator(J, RK.ah[stage,stage] * Δt)
- mul!(lin_du_tmp, M.J, uₙ)
-# mul!(lin_du_tmp1, J, u)
- res .= uₙ + RK.a[stage,1] * Δt * stages[1] + RK.a[stage,2] * Δt * stages[2] - Δt * RK.ah[stage,1] * jstages[1] - Δt * RK.ah[stage,2] * jstages[2]
- krylov_solve!(workspace, M, res, atol = 1e-6, rtol = 1e-6)
- @. u = workspace.x
- J = JacobianOperator(F!, du, u, p)
- mul!(jstages[stage], J, u)
- f2!(du, workspace.x, p, t + RK.c[stage] * Δt)
+ stages[stage] .= du .+ du_tmp - jstages[stage]
+ ustages[stage] .= u
+
+ elseif stage == 3
+
+ J = JacobianOperator(F!, du, ustages[2], p)
+ M = LMOperator(J, RK.ah[stage, stage] * Δt)
+ mul!(lin_du_tmp, M.J, uₙ)
+ # mul!(lin_du_tmp1, J, u)
+ res .= uₙ + RK.a[stage, 1] * Δt * stages[1] + RK.a[stage, 2] * Δt * stages[2] - Δt * RK.ah[stage, 1] * jstages[1] - Δt * RK.ah[stage, 2] * jstages[2]
+ krylov_solve!(workspace, M, res, atol = 1.0e-6, rtol = 1.0e-6)
+ @. u = workspace.x
+ J = JacobianOperator(F!, du, u, p)
+ mul!(jstages[stage], J, u)
+ f2!(du, workspace.x, p, t + RK.c[stage] * Δt)
f1!(du_tmp, workspace.x, p, t + RK.c[stage] * Δt)
- stages[stage] .= du .+ du_tmp - jstages[stage]
- ustages[stage] .= u
+ stages[stage] .= du .+ du_tmp - jstages[stage]
+ ustages[stage] .= u
- @. u = uₙ + RK.b[1] * Δt * stages[1] + RK.b[2] * Δt * stages[2] + RK.b[3] * Δt * stages[3] - RK.bh[1] * Δt * jstages[1] - RK.bh[2] * Δt * jstages[2] - RK.bh[3] * Δt * jstages[3]
- end
+ @. u = uₙ + RK.b[1] * Δt * stages[1] + RK.b[2] * Δt * stages[2] + RK.b[3] * Δt * stages[3] - RK.bh[1] * Δt * jstages[1] - RK.bh[2] * Δt * jstages[2] - RK.bh[3] * Δt * jstages[3]
+ end
end
stages(::SimpleLinearImplicitExplicitAlgorithm{N}) where {N} = N
@@ -117,7 +117,7 @@ mutable struct SimpleLinearImplicitExplicitOptions{Callback}
end
function RKTableau(alg::RKLSSPIMEX332)
-return RKLSSPIMEX332Tableau()
+ return RKLSSPIMEX332Tableau()
end
function RKLSSPIMEX332Tableau()
@@ -127,30 +127,30 @@ function RKLSSPIMEX332Tableau()
a[2, 1] = 0.5
a[3, 1] = 0.5
a[3, 2] = 0.5
-
+
b = zeros(Float64, nstage)
- b[1] = 1/3
- b[2] = 1/3
- b[3] = 1/3
+ b[1] = 1 / 3
+ b[2] = 1 / 3
+ b[3] = 1 / 3
c = zeros(Float64, nstage)
c[2] = 0.5
c[3] = 1.0
ah = zeros(Float64, nstage, nstage)
- ah[1, 1] = 1/4
- ah[2, 2] = 1/4
- ah[3, 1] = 1/3
- ah[3, 2] = 1/3
- ah[3, 3] = 1/3
-
+ ah[1, 1] = 1 / 4
+ ah[2, 2] = 1 / 4
+ ah[3, 1] = 1 / 3
+ ah[3, 2] = 1 / 3
+ ah[3, 3] = 1 / 3
+
bh = zeros(Float64, nstage)
- bh[1] = 1/3
- bh[2] = 1/3
- bh[3] = 1/3
+ bh[1] = 1 / 3
+ bh[2] = 1 / 3
+ bh[3] = 1 / 3
ch = zeros(Float64, nstage)
- ch[1] = 1/4
- ch[2] = 1/4
+ ch[1] = 1 / 4
+ ch[2] = 1 / 4
ch[3] = 1.0
return IMEXRKButcher(a, b, c, ah, bh, ch)
end
@@ -174,7 +174,7 @@ function LinearImplicitExplicitEulerTableau()
end
-function SimpleLinearImplicitExplicitOptions(callback, tspan; maxiters=typemax(Int), verbose=0, krylov_algo=:gmres, krylov_kwargs=(;), kwargs...)
+function SimpleLinearImplicitExplicitOptions(callback, tspan; maxiters = typemax(Int), verbose = 0, krylov_algo = :gmres, krylov_kwargs = (;), kwargs...)
return SimpleLinearImplicitExplicitOptions{typeof(callback)}(
callback, false, Inf, maxiters,
[last(tspan)],
@@ -185,18 +185,18 @@ function SimpleLinearImplicitExplicitOptions(callback, tspan; maxiters=typemax(I
end
mutable struct SimpleLinearImplicitExplicit{
- RealT<:Real,uType,Params,Sol,F,F1,F2,M,Alg<:SimpleLinearImplicitExplicitAlgorithm,
- SimpleLinearImplicitExplicitOptions,RKTableau,
-} <: AbstractTimeIntegrator
+ RealT <: Real, uType, Params, Sol, F, F1, F2, M, Alg <: SimpleLinearImplicitExplicitAlgorithm,
+ SimpleLinearImplicitExplicitOptions, RKTableau,
+ } <: AbstractTimeIntegrator
u::uType
du::uType
du_tmp::uType
lin_du_tmp::uType
lin_du_tmp1::uType
u_tmp::uType
- stages::NTuple{M,uType}
- ustages::NTuple{M,uType}
- jstages::NTuple{M,uType}
+ stages::NTuple{M, uType}
+ ustages::NTuple{M, uType}
+ jstages::NTuple{M, uType}
res::uType
t::RealT
dt::RealT # current time step
@@ -216,16 +216,16 @@ end
function Base.getproperty(integrator::SimpleLinearImplicitExplicit, field::Symbol)
if field === :stats
- return (naccept=getfield(integrator, :iter),)
+ return (naccept = getfield(integrator, :iter),)
end
# general fallback
return getfield(integrator, field)
end
function init(
- ode::ODEProblem, alg::SimpleLinearImplicitExplicitAlgorithm{N};
- dt, callback::Union{CallbackSet,Nothing}=nothing, kwargs...,
-) where {N}
+ ode::ODEProblem, alg::SimpleLinearImplicitExplicitAlgorithm{N};
+ dt, callback::Union{CallbackSet, Nothing} = nothing, kwargs...,
+ ) where {N}
u = copy(ode.u0)
du = zero(u)
res = zero(u)
@@ -236,12 +236,13 @@ function init(
t = first(ode.tspan)
iter = 0
integrator = SimpleLinearImplicitExplicit(
- u, du, copy(du),copy(du), copy(du), u_tmp,stages, ustages, jstages, res, t, dt, zero(dt), iter, ode.p,
- (prob=ode,), ode.f.f1, ode.f.f1, ode.f.f2, alg,
+ u, du, copy(du), copy(du), copy(du), u_tmp, stages, ustages, jstages, res, t, dt, zero(dt), iter, ode.p,
+ (prob = ode,), ode.f.f1, ode.f.f1, ode.f.f2, alg,
SimpleLinearImplicitExplicitOptions(
callback, ode.tspan;
kwargs...,
- ), false, RKTableau(alg))
+ ), false, RKTableau(alg)
+ )
# initialize callbacks
if callback isa CallbackSet
@@ -258,10 +259,10 @@ end
# Fakes `solve`: https://diffeq.sciml.ai/v6.8/basics/overview/#Solving-the-Problems-1
function solve(
- ode::ODEProblem, alg::SimpleLinearImplicitExplicitAlgorithm;
- dt, callback=nothing, kwargs...,
-)
- integrator = init(ode, alg, dt=dt, callback=callback; kwargs...)
+ ode::ODEProblem, alg::SimpleLinearImplicitExplicitAlgorithm;
+ dt, callback = nothing, kwargs...,
+ )
+ integrator = init(ode, alg, dt = dt, callback = callback; kwargs...)
# Start actual solve
return solve!(integrator)
@@ -299,7 +300,7 @@ function step!(integrator::SimpleLinearImplicitExplicit)
# if the next iteration would push the simulation beyond the end time, set dt accordingly
if integrator.t + integrator.dt > t_end ||
- isapprox(integrator.t + integrator.dt, t_end)
+ isapprox(integrator.t + integrator.dt, t_end)
integrator.dt = t_end - integrator.t
terminate!(integrator)
end
@@ -334,17 +335,18 @@ function step!(integrator::SimpleLinearImplicitExplicit)
end
function stage!(integrator, alg::RKLIMEX)
- F!(du, u, p) = integrator.f1(du, u, p, integrator.t) ## parabolic
- J = JacobianOperator(F!, integrator.du, integrator.u, integrator.p)
- M = LMOperator(J, integrator.dt)
+ F!(du, u, p) = integrator.f1(du, u, p, integrator.t) ## parabolic
+ J = JacobianOperator(F!, integrator.du, integrator.u, integrator.p)
+ M = LMOperator(J, integrator.dt)
kc = KrylovConstructor(integrator.res)
workspace = krylov_workspace(:gmres, kc)
- for stage in 1:stages(alg)
+ for stage in 1:stages(alg)
# Store the solution for each stage in stages
- ## For a split Problem we need to compute rhs_conservative and rhs_parabolic
- alg(integrator.res, integrator.u, integrator.dt, integrator.f1, integrator.f2, integrator.du, integrator.du_tmp, integrator.u_tmp, integrator.p, integrator.t, integrator.stages, integrator.ustages, integrator.jstages, stage, integrator.RK, M, integrator.lin_du_tmp, integrator.lin_du_tmp1, workspace)
+ ## For a split Problem we need to compute rhs_conservative and rhs_parabolic
+ alg(integrator.res, integrator.u, integrator.dt, integrator.f1, integrator.f2, integrator.du, integrator.du_tmp, integrator.u_tmp, integrator.p, integrator.t, integrator.stages, integrator.ustages, integrator.jstages, stage, integrator.RK, M, integrator.lin_du_tmp, integrator.lin_du_tmp1, workspace)
end
+ return
end
# get a cache where the RHS can be stored
diff --git a/libs/Implicit/src/skeleton.jl b/libs/Implicit/src/skeleton.jl
index 1b7bb22..6a614d7 100644
--- a/libs/Implicit/src/skeleton.jl
+++ b/libs/Implicit/src/skeleton.jl
@@ -7,11 +7,11 @@ struct RKImplicitExplicitEuler <: RKIMEX{1} end
function (::RKImplicitExplicitEuler)(res, uₙ, Δt, f1!, f2!, du, du_tmp, u, p, t, stages, stage, RK)
if stage == 1
# Stage 1:
- ## f2 is the conservative part
- ## f1 is the parabolic part
- f2!(du, uₙ, p, t + RK.c[stage] * Δt )
- f1!(du_tmp, u, p, t + RK.c[stage] * Δt )
- return res .= u .- uₙ .- RK.a[stage, stage] * Δt .* du - RK.a[stage,stage] * Δt * du_tmp
+ ## f2 is the conservative part
+ ## f1 is the parabolic part
+ f2!(du, uₙ, p, t + RK.c[stage] * Δt)
+ f1!(du_tmp, u, p, t + RK.c[stage] * Δt)
+ return res .= u .- uₙ .- RK.a[stage, stage] * Δt .* du - RK.a[stage, stage] * Δt * du_tmp
else
@. u = uₙ + RK.b[1] * Δt * stages[1]
end
@@ -54,7 +54,7 @@ function ImplicitExplicitEulerTableau()
end
-function SimpleImplicitExplicitOptions(callback, tspan; maxiters=typemax(Int), verbose=0, krylov_algo=:gmres, krylov_kwargs=(;), kwargs...)
+function SimpleImplicitExplicitOptions(callback, tspan; maxiters = typemax(Int), verbose = 0, krylov_algo = :gmres, krylov_kwargs = (;), kwargs...)
return SimpleImplicitExplicitOptions{typeof(callback)}(
callback, false, Inf, maxiters,
[last(tspan)],
@@ -65,14 +65,14 @@ function SimpleImplicitExplicitOptions(callback, tspan; maxiters=typemax(Int), v
end
mutable struct SimpleImplicitExplicit{
- RealT<:Real,uType,Params,Sol,F,F1,F2,M,Alg<:SimpleImplicitExplicitAlgorithm,
- SimpleImplicitExplicitOptions,RKTableau,
-} <: AbstractTimeIntegrator
+ RealT <: Real, uType, Params, Sol, F, F1, F2, M, Alg <: SimpleImplicitExplicitAlgorithm,
+ SimpleImplicitExplicitOptions, RKTableau,
+ } <: AbstractTimeIntegrator
u::uType
du::uType
du_tmp::uType
u_tmp::uType
- stages::NTuple{M,uType}
+ stages::NTuple{M, uType}
res::uType
t::RealT
dt::RealT # current time step
@@ -92,16 +92,16 @@ end
function Base.getproperty(integrator::SimpleImplicitExplicit, field::Symbol)
if field === :stats
- return (naccept=getfield(integrator, :iter),)
+ return (naccept = getfield(integrator, :iter),)
end
# general fallback
return getfield(integrator, field)
end
function init(
- ode::ODEProblem, alg::SimpleImplicitExplicitAlgorithm{N};
- dt, callback::Union{CallbackSet,Nothing}=nothing, kwargs...,
-) where {N}
+ ode::ODEProblem, alg::SimpleImplicitExplicitAlgorithm{N};
+ dt, callback::Union{CallbackSet, Nothing} = nothing, kwargs...,
+ ) where {N}
u = copy(ode.u0)
du = zero(u)
res = zero(u)
@@ -110,12 +110,13 @@ function init(
t = first(ode.tspan)
iter = 0
integrator = SimpleImplicitExplicit(
- u, du, copy(du), u_tmp, stages, res, t, dt, zero(dt), iter, ode.p,
- (prob=ode,), ode.f.f1, ode.f.f1, ode.f.f2, alg,
+ u, du, copy(du), u_tmp, stages, res, t, dt, zero(dt), iter, ode.p,
+ (prob = ode,), ode.f.f1, ode.f.f1, ode.f.f2, alg,
SimpleImplicitExplicitOptions(
callback, ode.tspan;
kwargs...,
- ), false, RKTableau(alg))
+ ), false, RKTableau(alg)
+ )
# initialize callbacks
if callback isa CallbackSet
@@ -132,10 +133,10 @@ end
# Fakes `solve`: https://diffeq.sciml.ai/v6.8/basics/overview/#Solving-the-Problems-1
function solve(
- ode::ODEProblem, alg::SimpleImplicitExplicitAlgorithm;
- dt, callback=nothing, kwargs...,
-)
- integrator = init(ode, alg, dt=dt, callback=callback; kwargs...)
+ ode::ODEProblem, alg::SimpleImplicitExplicitAlgorithm;
+ dt, callback = nothing, kwargs...,
+ )
+ integrator = init(ode, alg, dt = dt, callback = callback; kwargs...)
# Start actual solve
return solve!(integrator)
@@ -173,7 +174,7 @@ function step!(integrator::SimpleImplicitExplicit)
# if the next iteration would push the simulation beyond the end time, set dt accordingly
if integrator.t + integrator.dt > t_end ||
- isapprox(integrator.t + integrator.dt, t_end)
+ isapprox(integrator.t + integrator.dt, t_end)
integrator.dt = t_end - integrator.t
terminate!(integrator)
end
@@ -213,22 +214,23 @@ function stage!(integrator, alg::RKIMEX)
F! = nonlinear_problem(alg, integrator.f2)
# TODO: Pass in `stages[1:(stage-1)]` or full tuple?
_, stats = Ariadne.newton_krylov!(
- F!, integrator.u_tmp, (integrator.u, integrator.dt, integrator.f1, integrator.du, integrator.du_tmp, integrator.p, integrator.t, integrator.stages, stage, integrator.RK), integrator.res;
- verbose=integrator.opts.verbose, krylov_kwargs=integrator.opts.krylov_kwargs,
- algo=integrator.opts.algo, tol_abs=6.0e-6,
+ F!, integrator.u_tmp, (integrator.u, integrator.dt, integrator.f1, integrator.du, integrator.du_tmp, integrator.p, integrator.t, integrator.stages, stage, integrator.RK), integrator.res;
+ verbose = integrator.opts.verbose, krylov_kwargs = integrator.opts.krylov_kwargs,
+ algo = integrator.opts.algo, tol_abs = 6.0e-6,
)
@assert stats.solved
# Store the solution for each stage in stages
- ## For a split Problem we need to compute rhs_conservative and rhs_parabolic
+ ## For a split Problem we need to compute rhs_conservative and rhs_parabolic
integrator.f2(integrator.du, integrator.u_tmp, integrator.p, integrator.t + integrator.RK.c[stage] * integrator.dt)
- integrator.stages[stage] .= integrator.du
+ integrator.stages[stage] .= integrator.du
integrator.f1(integrator.du, integrator.u_tmp, integrator.p, integrator.t + integrator.RK.c[stage] * integrator.dt)
integrator.stages[stage] .+= integrator.du
- if stage == stages(alg)
+ if stage == stages(alg)
alg(integrator.res, integrator.u, integrator.dt, integrator.f1, integrator.f2, integrator.du, integrator.du_tmp, integrator.u_tmp, integrator.p, integrator.t, integrator.stages, stage + 1, integrator.RK)
end
end
+ return
end
# get a cache where the RHS can be stored |
vchuravy
reviewed
Jul 17, 2025
vchuravy
reviewed
Jul 17, 2025
libs/Implicit/src/linear_imex.jl
Outdated
| f1!(du_tmp, uₙ, p, t + RK.c[stage] * Δt) | ||
|
|
||
| res .= uₙ .+ RK.a[stage, stage] * Δt .* (du .+ du_tmp .- lin_du_tmp) | ||
| krylov_solve!(workspace, M, copy(res)) |
MarcoArtiano
changed the title
Collaborator
Author
|
This has been written already in a general form, that it already works with user defined splitting using the SemiDiscretizationSplit struct defined in trixi-framework/Trixi.jl#2058 |
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
Add this suggestion to a batch that can be applied as a single commit.
This suggestion is invalid because no changes were made to the code.
Suggestions cannot be applied while the pull request is closed.
Suggestions cannot be applied while viewing a subset of changes.
Only one suggestion per line can be applied in a batch.
Add this suggestion to a batch that can be applied as a single commit.
Applying suggestions on deleted lines is not supported.
You must change the existing code in this line in order to create a valid suggestion.
Outdated suggestions cannot be applied.
This suggestion has been applied or marked resolved.
Suggestions cannot be applied from pending reviews.
Suggestions cannot be applied on multi-line comments.
Suggestions cannot be applied while the pull request is queued to merge.
Suggestion cannot be applied right now. Please check back later.
No description provided.