Skip to content

Commit 1c29507

Browse files
DanielDoehringJoshuaLampertranocha
authored
Use Jacobian-free Newton-Krylov for SDIRK time integrators (#2505)
* example linear newton krylov * Nonlinear example * comment * Update test/Project.toml * Update test/Project.toml * check different sciml * v * v * v * v * v * v * v * v * v * v * Update test/Project.toml * Update Project.toml * ci * compat * v * v * v * v * v * v * Update test/Project.toml Co-authored-by: Joshua Lampert <51029046+JoshuaLampert@users.noreply.github.com> * Update test/Project.toml Co-authored-by: Joshua Lampert <51029046+JoshuaLampert@users.noreply.github.com> * Update test/Project.toml * Update test/Project.toml * Update test/test_parabolic_2d.jl * up * Update test/test_parabolic_2d.jl * Update test/test_parabolic_2d.jl * tolerances * ci vals * ode solver tols * Update test/test_parabolic_2d.jl --------- Co-authored-by: Joshua Lampert <51029046+JoshuaLampert@users.noreply.github.com> Co-authored-by: Hendrik Ranocha <ranocha@users.noreply.github.com>
1 parent aa192cc commit 1c29507

6 files changed

Lines changed: 298 additions & 3 deletions

File tree

‎examples/p4est_2d_dgsem/elixir_navierstokes_viscous_shock_newton_krylov.jl‎

Lines changed: 183 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,183 @@
1+
using Trixi
2+
3+
using OrdinaryDiffEqSDIRK
4+
using LinearSolve # For Jacobian-free Newton-Krylov (GMRES) solver
5+
using ADTypes # For automatic differentiation via finite differences
6+
7+
# This is the classic 1D viscous shock wave problem with analytical solution
8+
# for a special value of the Prandtl number.
9+
# The original references are:
10+
#
11+
# - R. Becker (1922)
12+
# Stoßwelle und Detonation.
13+
# [DOI: 10.1007/BF01329605](https://doi.org/10.1007/BF01329605)
14+
#
15+
# English translations:
16+
# Impact waves and detonation. Part I.
17+
# https://ntrs.nasa.gov/api/citations/19930090862/downloads/19930090862.pdf
18+
# Impact waves and detonation. Part II.
19+
# https://ntrs.nasa.gov/api/citations/19930090863/downloads/19930090863.pdf
20+
#
21+
# - M. Morduchow, P. A. Libby (1949)
22+
# On a Complete Solution of the One-Dimensional Flow Equations
23+
# of a Viscous, Head-Conducting, Compressible Gas
24+
# [DOI: 10.2514/8.11882](https://doi.org/10.2514/8.11882)
25+
#
26+
#
27+
# The particular problem considered here is described in
28+
# - L. G. Margolin, J. M. Reisner, P. M. Jordan (2017)
29+
# Entropy in self-similar shock profiles
30+
# [DOI: 10.1016/j.ijnonlinmec.2017.07.003](https://doi.org/10.1016/j.ijnonlinmec.2017.07.003)
31+
32+
### Fixed parameters ###
33+
34+
# Special value for which nonlinear solver can be omitted
35+
# Corresponds essentially to fixing the Mach number
36+
alpha = 0.5
37+
# We want kappa = cp * mu = mu_bar to ensure constant enthalpy
38+
prandtl_number() = 3 / 4
39+
40+
### Free choices: ###
41+
gamma() = 5 / 3
42+
43+
# In Margolin et al., the Navier-Stokes equations are given for an
44+
# isotropic stress tensor τ, i.e., ∇ ⋅ τ = μ Δu
45+
mu_isotropic() = 0.15
46+
mu_bar() = mu_isotropic() / (gamma() - 1) # Re-scaled viscosity
47+
48+
rho_0() = 1
49+
v() = 1 # Shock speed
50+
51+
domain_length = 4.0
52+
53+
### Derived quantities ###
54+
55+
Ma() = 2 / sqrt(3 - gamma()) # Mach number for alpha = 0.5
56+
c_0() = v() / Ma() # Speed of sound ahead of the shock
57+
58+
# From constant enthalpy condition
59+
p_0() = c_0()^2 * rho_0() / gamma()
60+
61+
l() = mu_bar() / (rho_0() * v()) * 2 * gamma() / (gamma() + 1) # Appropriate length scale
62+
63+
"""
64+
initial_condition_viscous_shock(x, t, equations)
65+
66+
Classic 1D viscous shock wave problem with analytical solution
67+
for a special value of the Prandtl number.
68+
The version implemented here is described in
69+
- L. G. Margolin, J. M. Reisner, P. M. Jordan (2017)
70+
Entropy in self-similar shock profiles
71+
[DOI: 10.1016/j.ijnonlinmec.2017.07.003](https://doi.org/10.1016/j.ijnonlinmec.2017.07.003)
72+
"""
73+
function initial_condition_viscous_shock(x, t, equations)
74+
y = x[1] - v() * t # Translated coordinate
75+
76+
# Coordinate transformation. See eq. (33) in Margolin et al. (2017)
77+
chi = 2 * exp(y / (2 * l()))
78+
79+
w = 1 + 1 / (2 * chi^2) * (1 - sqrt(1 + 2 * chi^2))
80+
81+
rho = rho_0() / w
82+
u = v() * (1 - w)
83+
p = p_0() * 1 / w * (1 + (gamma() - 1) / 2 * Ma()^2 * (1 - w^2))
84+
85+
return prim2cons(SVector(rho, u, 0, p), equations)
86+
end
87+
initial_condition = initial_condition_viscous_shock
88+
89+
###############################################################################
90+
# semidiscretization of the ideal compressible Navier-Stokes equations
91+
92+
equations = CompressibleEulerEquations2D(gamma())
93+
94+
# Trixi implements the stress tensor in deviatoric form, thus we need to
95+
# convert the "isotropic viscosity" to the "deviatoric viscosity"
96+
mu_deviatoric() = mu_bar() * 3 / 4
97+
equations_parabolic = CompressibleNavierStokesDiffusion2D(equations, mu = mu_deviatoric(),
98+
Prandtl = prandtl_number(),
99+
gradient_variables = GradientVariablesEntropy())
100+
101+
solver = DGSEM(polydeg = 3, surface_flux = flux_hlle)
102+
103+
coordinates_min = (-domain_length / 2, -domain_length / 2)
104+
coordinates_max = (domain_length / 2, domain_length / 2)
105+
106+
trees_per_dimension = (12, 3)
107+
mesh = P4estMesh(trees_per_dimension,
108+
polydeg = 3, initial_refinement_level = 0,
109+
coordinates_min = coordinates_min, coordinates_max = coordinates_max,
110+
periodicity = (false, true))
111+
112+
### Inviscid boundary conditions ###
113+
114+
# Prescribe pure influx based on initial conditions
115+
function boundary_condition_inflow(u_inner, normal_direction::AbstractVector, x, t,
116+
surface_flux_function,
117+
equations::CompressibleEulerEquations2D)
118+
u_cons = initial_condition_viscous_shock(x, t, equations)
119+
flux = Trixi.flux(u_cons, normal_direction, equations)
120+
121+
return flux
122+
end
123+
124+
boundary_conditions = Dict(:x_neg => boundary_condition_inflow,
125+
:x_pos => boundary_condition_do_nothing)
126+
127+
### Viscous boundary conditions ###
128+
# For the viscous BCs, we use the known analytical solution
129+
velocity_bc = NoSlip() do x, t, equations_parabolic
130+
Trixi.velocity(initial_condition_viscous_shock(x,
131+
t,
132+
equations_parabolic),
133+
equations_parabolic)
134+
end
135+
136+
heat_bc = Isothermal() do x, t, equations_parabolic
137+
Trixi.temperature(initial_condition_viscous_shock(x,
138+
t,
139+
equations_parabolic),
140+
equations_parabolic)
141+
end
142+
143+
boundary_condition_parabolic = BoundaryConditionNavierStokesWall(velocity_bc, heat_bc)
144+
145+
boundary_conditions_parabolic = Dict(:x_neg => boundary_condition_parabolic,
146+
:x_pos => boundary_condition_parabolic)
147+
148+
semi = SemidiscretizationHyperbolicParabolic(mesh, (equations, equations_parabolic),
149+
initial_condition, solver;
150+
boundary_conditions = (boundary_conditions,
151+
boundary_conditions_parabolic))
152+
153+
###############################################################################
154+
# ODE solvers, callbacks etc.
155+
156+
tspan = (0.0, 1.0)
157+
ode = semidiscretize(semi, tspan)
158+
159+
summary_callback = SummaryCallback()
160+
alive_callback = AliveCallback(alive_interval = 1)
161+
analysis_callback = AnalysisCallback(semi, interval = 100)
162+
163+
callbacks = CallbackSet(summary_callback, alive_callback, analysis_callback)
164+
165+
###############################################################################
166+
# run the simulation
167+
168+
# Tolerances for GMRES residual, see https://jso.dev/Krylov.jl/stable/solvers/unsymmetric/#Krylov.gmres
169+
atol_lin_solve = 1e-5
170+
rtol_lin_solve = 1e-5
171+
172+
# Jacobian-free Newton-Krylov (GMRES) solver
173+
linsolve = KrylovJL_GMRES(atol = atol_lin_solve, rtol = rtol_lin_solve)
174+
175+
# Use (diagonally) implicit Runge-Kutta, see
176+
# https://docs.sciml.ai/DiffEqDocs/stable/tutorials/advanced_ode_example/#Using-Jacobian-Free-Newton-Krylov
177+
ode_alg = Kvaerno4(autodiff = AutoFiniteDiff(), linsolve = linsolve)
178+
179+
atol_ode_solve = 1e-4
180+
rtol_ode_solve = 1e-4
181+
sol = solve(ode, ode_alg;
182+
abstol = atol_ode_solve, reltol = rtol_ode_solve,
183+
ode_default_options()..., callback = callbacks);

‎examples/tree_1d_dgsem/elixir_diffusion_ldg.jl‎

Lines changed: 2 additions & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -44,15 +44,14 @@ boundary_conditions_parabolic = boundary_condition_periodic
4444
solver_parabolic = ViscousFormulationLocalDG()
4545
semi = SemidiscretizationHyperbolicParabolic(mesh, (equations, equations_parabolic),
4646
initial_condition,
47-
solver;
48-
solver_parabolic,
47+
solver; solver_parabolic,
4948
boundary_conditions = (boundary_conditions,
5049
boundary_conditions_parabolic))
5150

5251
###############################################################################
5352
# ODE solvers, callbacks etc.
5453

55-
# Create ODE problem with time span from 0.0 to 1.0
54+
# Create ODE problem with time span from 0.0 to 0.1
5655
tspan = (0.0, 0.1)
5756
ode = semidiscretize(semi, tspan)
5857

‎examples/tree_1d_dgsem/elixir_diffusion_ldg_newton_krylov.jl‎

Lines changed: 71 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,71 @@
1+
using Trixi
2+
3+
using OrdinaryDiffEqSDIRK
4+
using LinearSolve # For Jacobian-free Newton-Krylov (GMRES) solver
5+
using ADTypes # For automatic differentiation via finite differences
6+
7+
###############################################################################
8+
# semidiscretization of the linear (advection) diffusion equation
9+
10+
advection_velocity = 0.0 # Note: This renders the equation mathematically purely parabolic
11+
equations = LinearScalarAdvectionEquation1D(advection_velocity)
12+
diffusivity() = 0.5
13+
equations_parabolic = LaplaceDiffusion1D(diffusivity(), equations)
14+
15+
# surface flux does not matter for pure diffusion problem
16+
solver = DGSEM(polydeg = 3, surface_flux = flux_central)
17+
18+
coordinates_min = -convert(Float64, pi)
19+
coordinates_max = convert(Float64, pi)
20+
21+
mesh = TreeMesh(coordinates_min, coordinates_max,
22+
initial_refinement_level = 4,
23+
n_cells_max = 30_000)
24+
25+
function initial_condition_pure_diffusion_1d_convergence_test(x, t,
26+
equation)
27+
nu = diffusivity()
28+
c = 0
29+
A = 1
30+
omega = 1
31+
scalar = c + A * sin(omega * sum(x)) * exp(-nu * omega^2 * t)
32+
return SVector(scalar)
33+
end
34+
initial_condition = initial_condition_pure_diffusion_1d_convergence_test
35+
36+
solver_parabolic = ViscousFormulationLocalDG()
37+
semi = SemidiscretizationHyperbolicParabolic(mesh, (equations, equations_parabolic),
38+
initial_condition,
39+
solver; solver_parabolic)
40+
41+
###############################################################################
42+
# ODE solvers, callbacks etc.
43+
44+
tspan = (0.0, 2.0)
45+
ode = semidiscretize(semi, tspan)
46+
47+
summary_callback = SummaryCallback()
48+
analysis_callback = AnalysisCallback(semi, interval = 10)
49+
alive_callback = AliveCallback(alive_interval = 1)
50+
51+
callbacks = CallbackSet(summary_callback, analysis_callback, alive_callback)
52+
53+
###############################################################################
54+
# run the simulation
55+
56+
# Tolerances for GMRES residual, see https://jso.dev/Krylov.jl/stable/solvers/unsymmetric/#Krylov.gmres
57+
atol_lin_solve = 1e-6
58+
rtol_lin_solve = 1e-5
59+
60+
# Jacobian-free Newton-Krylov (GMRES) solver
61+
linsolve = KrylovJL_GMRES(atol = atol_lin_solve, rtol = rtol_lin_solve)
62+
63+
# Use (diagonally) implicit Runge-Kutta, see
64+
# https://docs.sciml.ai/DiffEqDocs/stable/tutorials/advanced_ode_example/#Using-Jacobian-Free-Newton-Krylov
65+
ode_alg = KenCarp47(autodiff = AutoFiniteDiff(), linsolve = linsolve)
66+
67+
atol_ode_solve = 1e-5
68+
rtol_ode_solve = 1e-4
69+
sol = solve(ode, ode_alg;
70+
abstol = atol_ode_solve, reltol = rtol_ode_solve,
71+
ode_default_options()..., callback = callbacks);

‎test/Project.toml‎

Lines changed: 2 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -12,6 +12,7 @@ ECOS = "e2685f51-7e38-5353-a97d-a921fd2c8199"
1212
ExplicitImports = "7d51a73a-1435-4ff3-83d9-f097790105c7"
1313
ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210"
1414
LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
15+
LinearSolve = "7ed4a6bd-45f5-4d41-b270-4a48e9bafcae"
1516
MPI = "da04e1cc-30fd-572f-bb4f-1f8673147195"
1617
NLsolve = "2774e3e8-f4cf-5e23-947b-6d7e65073b56"
1718
OrdinaryDiffEqFeagin = "101fe9f7-ebb6-4678-b671-3a81e7194747"
@@ -43,6 +44,7 @@ ECOS = "1.1.2"
4344
ExplicitImports = "1.0.1"
4445
ForwardDiff = "0.10.36, 1"
4546
LinearAlgebra = "1"
47+
LinearSolve = "2.36.1, 3"
4648
MPI = "0.20.6"
4749
NLsolve = "4.5.1"
4850
OrdinaryDiffEqFeagin = "1"

‎test/test_parabolic_1d.jl‎

Lines changed: 14 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -58,6 +58,20 @@ end
5858
end
5959
end
6060

61+
@trixi_testset "TreeMesh1D: elixir_diffusion_ldg_newton_krylov.jl" begin
62+
@test_trixi_include(joinpath(examples_dir(), "tree_1d_dgsem",
63+
"elixir_diffusion_ldg_newton_krylov.jl"),
64+
l2=[4.2710445174631516e-6], linf=[2.28491835256861e-5])
65+
# Ensure that we do not have excessive memory allocations
66+
# (e.g., from type instabilities)
67+
let
68+
t = sol.t[end]
69+
u_ode = sol.u[end]
70+
du_ode = similar(u_ode)
71+
@test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 1000
72+
end
73+
end
74+
6175
@trixi_testset "TreeMesh1D: elixir_advection_diffusion_restart.jl" begin
6276
@test_trixi_include(joinpath(examples_dir(), "tree_1d_dgsem",
6377
"elixir_advection_diffusion_restart.jl"),

‎test/test_parabolic_2d.jl‎

Lines changed: 26 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -866,6 +866,32 @@ end
866866
end
867867
end
868868

869+
@trixi_testset "P4estMesh2D: elixir_navierstokes_viscous_shock_newton_krylov.jl" begin
870+
@test_trixi_include(joinpath(examples_dir(), "p4est_2d_dgsem",
871+
"elixir_navierstokes_viscous_shock_newton_krylov.jl"),
872+
tspan=(0.0, 0.1),
873+
l2=[
874+
3.468233560427797e-5,
875+
2.64864594855224e-5,
876+
7.879490760481979e-10,
877+
2.8748482665365446e-5
878+
],
879+
linf=[
880+
0.00018754529350140103,
881+
0.00014045634087878067,
882+
9.043610782328732e-9,
883+
0.00014499382160382268
884+
])
885+
# Ensure that we do not have excessive memory allocations
886+
# (e.g., from type instabilities)
887+
let
888+
t = sol.t[end]
889+
u_ode = sol.u[end]
890+
du_ode = similar(u_ode)
891+
@test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 1000
892+
end
893+
end
894+
869895
@trixi_testset "elixir_navierstokes_SD7003airfoil.jl" begin
870896
@test_trixi_include(joinpath(examples_dir(), "p4est_2d_dgsem",
871897
"elixir_navierstokes_SD7003airfoil.jl"),

0 commit comments

Comments
 (0)