28769960495

Committed 06 Jul 2026 05:29AM UTC coverage: 96.853%. Remained the same

Build # 28769960495

Build Type

push

github

Committed by

ranocha

Commit Message

set development version to v0.16.25-DEV

Coverage Stats

49584 of 51195 relevant lines covered (96.85%)

104745925.55 hits per line

Source File
Press 'n' to go to next uncovered line, 'b' for previous

100.0

/examples/tree_2d_dgsem/elixir_diffusion_steady_state_linear_map.jl

using Trixi

###############################################################################

# Build pure diffusion (Laplace) operator
diffusivity() = 1
equations = LinearDiffusionEquation2D(diffusivity())
solver = DGSEM(polydeg = 5)

coordinates_min = (0.0, 0.0)
coordinates_max = (1.0, 1.0)
mesh = TreeMesh(coordinates_min, coordinates_max,
                initial_refinement_level = 2,
                n_cells_max = 80_000,
                periodicity = false)

# Analytical/continuous steady-state solution
function continuous_solution(x, t, equations)
    a = (1 - cosh(2 * pi)) / (sinh(2 * pi))

    u_sol = (cosh(2 * pi * x[2]) + a * sinh(2 * pi * x[2])) * sinpi(2 * x[1])
    return SVector(u_sol)
end
initial_condition = continuous_solution

function bc_homogeneous(x, t, equations)
    return SVector(0)
end
bc_homogeneous_dirichlet = BoundaryConditionDirichlet(bc_homogeneous)

function bc_sin(x, t, equations)
    return SVector(sinpi(2 * x[1]))
end
bc_sin_dirichlet = BoundaryConditionDirichlet(bc_sin)

boundary_conditions = (; x_neg = bc_homogeneous_dirichlet,
                       y_neg = bc_sin_dirichlet,
                       y_pos = bc_sin_dirichlet,
                       x_pos = bc_homogeneous_dirichlet)

# Build the parabolic semidiscretization with a local DG formulation for diffusion
semi = SemidiscretizationParabolic(mesh, equations, initial_condition, solver;
                                   solver_parabolic = ParabolicFormulationLocalDG(),
                                   boundary_conditions = boundary_conditions)

# Note that `linear_structure` does not access the `initial_condition`/steady-state solution
A_map, b = linear_structure(semi)

# Direct solve, with explicit sparse matrix construction.
# Has trouble due to poor conditioning, visible in top right corner
#=
using SparseArrays
A_matrix = sparse(A_map)
u_ls = A_matrix \ b
=#

# Iterative solve, works directly on the linear map, no explicit matrix construction required!
using Krylov

# This solves the Laplace equation (i.e., steady-state diffusion/heat equation).
# Note that we do not use CG since the operator `A_map` itself is only
# symmetric and positive definite with respect to the inner product induced by
# the DGSEM quadrature, not the standard Euclidean inner product. Standard CG
# would require multiplying the result of `A_map * u` (and `b`) by the mass matrix.
u_ls, stats = gmres(A_map, b, atol = 1.0e-11, rtol = 1.0e-10)

###############################################################################

# Construct the ODE problem for easy plotting and comparison to analytical solution
# Note that we only solve a steady-state problem, i.e., `tspan` has a range of zero.
tspan = (0.0, 0.0)
ode = semidiscretize(semi, tspan)

summary_callback = SummaryCallback()
# Analysis callback quantifies discretization/interpolation error of the exact solution
analysis_callback = AnalysisCallback(semi)
callbacks = CallbackSet(summary_callback,
                        analysis_callback)

# Choice of ODE Solver does not matter here
using OrdinaryDiffEqLowStorageRK

sol = solve(ode, CarpenterKennedy2N54(williamson_condition = false);
            dt = 1e-4,
            ode_default_options()..., callback = callbacks);

# Check interpolation errors due to choice of number of cells & polynomial degree
interpolation_errors = analysis_callback(sol)

using Plots
# Plot analytical solution (interpolated initial condition)
plot(sol)

# Inject linear system solution for plotting & error computation
sol.u[end] = u_ls
plot(sol)

# Check linear system solution errors
linear_solution_errors = analysis_callback(sol)

1	using Trixi
2
3	###############################################################################
4
5	# Build pure diffusion (Laplace) operator
6	diffusivity() = 1	1✔
7	equations = LinearDiffusionEquation2D(diffusivity())
8	solver = DGSEM(polydeg = 5)
9
10	coordinates_min = (0.0, 0.0)
11	coordinates_max = (1.0, 1.0)
12	mesh = TreeMesh(coordinates_min, coordinates_max,
13	initial_refinement_level = 2,
14	n_cells_max = 80_000,
15	periodicity = false)
16
17	# Analytical/continuous steady-state solution
18	function continuous_solution(x, t, equations)	14,128✔
19	a = (1 - cosh(2 * pi)) / (sinh(2 * pi))	14,128✔
20
21	u_sol = (cosh(2 * pi * x[2]) + a * sinh(2 * pi * x[2])) * sinpi(2 * x[1])	14,128✔
22	return SVector(u_sol)	14,128✔
23	end
24	initial_condition = continuous_solution
25
26	function bc_homogeneous(x, t, equations)	6,096✔
27	return SVector(0)	6,096✔
28	end
29	bc_homogeneous_dirichlet = BoundaryConditionDirichlet(bc_homogeneous)
30
31	function bc_sin(x, t, equations)	6,096✔
32	return SVector(sinpi(2 * x[1]))	6,096✔
33	end
34	bc_sin_dirichlet = BoundaryConditionDirichlet(bc_sin)
35
36	boundary_conditions = (; x_neg = bc_homogeneous_dirichlet,
37	y_neg = bc_sin_dirichlet,
38	y_pos = bc_sin_dirichlet,
39	x_pos = bc_homogeneous_dirichlet)
40
41	# Build the parabolic semidiscretization with a local DG formulation for diffusion
42	semi = SemidiscretizationParabolic(mesh, equations, initial_condition, solver;
43	solver_parabolic = ParabolicFormulationLocalDG(),
44	boundary_conditions = boundary_conditions)
45
46	# Note that `linear_structure` does not access the `initial_condition`/steady-state solution
47	A_map, b = linear_structure(semi)
48
49	# Direct solve, with explicit sparse matrix construction.
50	# Has trouble due to poor conditioning, visible in top right corner
51	#=
52	using SparseArrays
53	A_matrix = sparse(A_map)
54	u_ls = A_matrix \ b
55	=#
56
57	# Iterative solve, works directly on the linear map, no explicit matrix construction required!
58	using Krylov
59
60	# This solves the Laplace equation (i.e., steady-state diffusion/heat equation).
61	# Note that we do not use CG since the operator `A_map` itself is only
62	# symmetric and positive definite with respect to the inner product induced by
63	# the DGSEM quadrature, not the standard Euclidean inner product. Standard CG
64	# would require multiplying the result of `A_map * u` (and `b`) by the mass matrix.
65	u_ls, stats = gmres(A_map, b, atol = 1.0e-11, rtol = 1.0e-10)
66
67	###############################################################################
68
69	# Construct the ODE problem for easy plotting and comparison to analytical solution
70	# Note that we only solve a steady-state problem, i.e., `tspan` has a range of zero.
71	tspan = (0.0, 0.0)
72	ode = semidiscretize(semi, tspan)
73
74	summary_callback = SummaryCallback()
75	# Analysis callback quantifies discretization/interpolation error of the exact solution
76	analysis_callback = AnalysisCallback(semi)
77	callbacks = CallbackSet(summary_callback,
78	analysis_callback)
79
80	# Choice of ODE Solver does not matter here
81	using OrdinaryDiffEqLowStorageRK
82
83	sol = solve(ode, CarpenterKennedy2N54(williamson_condition = false);
84	dt = 1e-4,
85	ode_default_options()..., callback = callbacks);
86
87	# Check interpolation errors due to choice of number of cells & polynomial degree
88	interpolation_errors = analysis_callback(sol)
89
90	using Plots
91	# Plot analytical solution (interpolated initial condition)
92	plot(sol)
93
94	# Inject linear system solution for plotting & error computation
95	sol.u[end] = u_ls
96	plot(sol)
97
98	# Check linear system solution errors
99	linear_solution_errors = analysis_callback(sol)

trixi-framework / Trixi.jl / 28769960495

Source File Press 'n' to go to next uncovered line, 'b' for previous

Source File
Press 'n' to go to next uncovered line, 'b' for previous