Book a Demo!
CoCalc Logo Icon
StoreFeaturesDocsShareSupportNewsAboutPoliciesSign UpSign In
trixi-framework
GitHub Repository: trixi-framework/Trixi.jl
Path: blob/main/examples/tree_1d_dgsem/elixir_advection_diffusion_ldg.jl
2055 views
1
using OrdinaryDiffEqLowStorageRK
2
using Trixi
3
4
###############################################################################
5
# semidiscretization of the linear advection diffusion equation
6
7
advection_velocity = 0.1
8
equations = LinearScalarAdvectionEquation1D(advection_velocity)
9
diffusivity() = 0.5
10
equations_parabolic = LaplaceDiffusion1D(diffusivity(), equations)
11
12
# Create DG solver with polynomial degree = 3 and (local) Lax-Friedrichs/Rusanov flux as surface flux
13
solver = DGSEM(polydeg = 3, surface_flux = flux_lax_friedrichs)
14
15
coordinates_min = -convert(Float64, pi) # minimum coordinate
16
coordinates_max = convert(Float64, pi) # maximum coordinate
17
18
# Create a uniformly refined mesh with periodic boundaries
19
mesh = TreeMesh(coordinates_min, coordinates_max,
20
initial_refinement_level = 4,
21
n_cells_max = 30_000, # set maximum capacity of tree data structure
22
periodicity = true)
23
24
function x_trans_periodic(x, domain_length = SVector(oftype(x[1], 2 * pi)),
25
center = SVector(oftype(x[1], 0)))
26
x_normalized = x .- center
27
x_shifted = x_normalized .% domain_length
28
x_offset = ((x_shifted .< -0.5f0 * domain_length) -
29
(x_shifted .> 0.5f0 * domain_length)) .*
30
domain_length
31
return center + x_shifted + x_offset
32
end
33
34
# Define initial condition
35
function initial_condition_diffusive_convergence_test(x, t,
36
equation::LinearScalarAdvectionEquation1D)
37
# Store translated coordinate for easy use of exact solution
38
x_trans = x_trans_periodic(x - equation.advection_velocity * t)
39
40
nu = diffusivity()
41
c = 0
42
A = 1
43
omega = 1
44
scalar = c + A * sin(omega * sum(x_trans)) * exp(-nu * omega^2 * t)
45
return SVector(scalar)
46
end
47
initial_condition = initial_condition_diffusive_convergence_test
48
49
# define periodic boundary conditions everywhere
50
boundary_conditions = boundary_condition_periodic
51
boundary_conditions_parabolic = boundary_condition_periodic
52
53
# A semidiscretization collects data structures and functions for the spatial discretization
54
semi = SemidiscretizationHyperbolicParabolic(mesh, (equations, equations_parabolic),
55
initial_condition,
56
solver;
57
solver_parabolic = ViscousFormulationLocalDG(),
58
boundary_conditions = (boundary_conditions,
59
boundary_conditions_parabolic))
60
61
###############################################################################
62
# ODE solvers, callbacks etc.
63
64
# Create ODE problem with time span from 0.0 to 1.0
65
tspan = (0.0, 1.0)
66
ode = semidiscretize(semi, tspan)
67
68
# At the beginning of the main loop, the SummaryCallback prints a summary of the simulation setup
69
# and resets the timers
70
summary_callback = SummaryCallback()
71
72
# The AnalysisCallback allows to analyse the solution in regular intervals and prints the results
73
analysis_callback = AnalysisCallback(semi, interval = 100)
74
75
# The AliveCallback prints short status information in regular intervals
76
alive_callback = AliveCallback(analysis_interval = 100)
77
78
# The SaveRestartCallback allows to save a file from which a Trixi.jl simulation can be restarted
79
save_restart = SaveRestartCallback(interval = 100,
80
save_final_restart = true)
81
82
# Create a CallbackSet to collect all callbacks such that they can be passed to the ODE solver
83
callbacks = CallbackSet(summary_callback, analysis_callback, alive_callback, save_restart)
84
85
###############################################################################
86
# run the simulation
87
88
# OrdinaryDiffEq's `solve` method evolves the solution in time and executes the passed callbacks
89
# For CI purposes, we use fixed time-stepping for this elixir.
90
sol = solve(ode, RDPK3SpFSAL35(); dt = 1.0e-3, adaptive = false,
91
ode_default_options()..., callback = callbacks)
92
93