Book a Demo!
CoCalc Logo Icon
StoreFeaturesDocsShareSupportNewsAboutPoliciesSign UpSign In
trixi-framework
GitHub Repository: trixi-framework/Trixi.jl
 Path: blob/main/examples/dgmulti_1d/elixir_advection_diffusion_sbp.jl
5586 views
1
using OrdinaryDiffEqLowStorageRK

2
using Trixi

3


4
dg = DGMulti(polydeg = 3, element_type = Line(), approximation_type = GaussSBP(),

5
             surface_integral = SurfaceIntegralWeakForm(flux_lax_friedrichs))

6


7
advection_velocity = 0.1

8
equations = LinearScalarAdvectionEquation1D(advection_velocity)

9
diffusivity() = 5.0e-2

10
equations_parabolic = LaplaceDiffusion1D(diffusivity(), equations)

11


12
function x_trans_periodic(x, domain_length = SVector(oftype(x[1], 2 * pi)),

13
                          center = SVector(oftype(x[1], 0)))

14
    x_normalized = x .- center

15
    x_shifted = x_normalized .% domain_length

16
    x_offset = ((x_shifted .< -0.5f0 * domain_length) -

17
                (x_shifted .> 0.5f0 * domain_length)) .*

18
               domain_length

19
    return center + x_shifted + x_offset

20
end

21


22
# Define initial condition

23
function initial_condition_diffusive_convergence_test(x, t,

24
                                                      equation::LinearScalarAdvectionEquation1D)

25
    # Store translated coordinate for easy use of exact solution

26
    x_trans = x_trans_periodic(x - equation.advection_velocity * t)

27


28
    nu = diffusivity()

29
    c = 0

30
    A = 1

31
    omega = 1

32
    scalar = c + A * sin(omega * sum(x_trans)) * exp(-nu * omega^2 * t)

33
    return SVector(scalar)

34
end

35
initial_condition = initial_condition_diffusive_convergence_test

36


37
cells_per_dimension = (12,)

38
mesh = DGMultiMesh(dg, cells_per_dimension,

39
                   coordinates_min = (-pi,), coordinates_max = (pi,),

40
                   periodicity = true)

41


42
# define periodic boundary conditions everywhere

43
boundary_conditions = boundary_condition_periodic

44
boundary_conditions_parabolic = boundary_condition_periodic

45


46
semi = SemidiscretizationHyperbolicParabolic(mesh, (equations, equations_parabolic),

47
                                             initial_condition, dg;

48
                                             boundary_conditions = (boundary_conditions,

49
                                                                    boundary_conditions_parabolic))

50


51
tspan = (0.0, 1.0)

52
ode = semidiscretize(semi, tspan)

53


54
summary_callback = SummaryCallback()

55


56
alive_callback = AliveCallback(alive_interval = 10)

57


58
analysis_interval = 100

59
analysis_callback = AnalysisCallback(semi, interval = analysis_interval, uEltype = real(dg))

60


61
callbacks = CallbackSet(summary_callback, alive_callback,

62
                        analysis_callback)

63


64
###############################################################################

65
# run the simulation

66


67
time_int_tol = 1e-6

68
sol = solve(ode, RDPK3SpFSAL49(); abstol = time_int_tol, reltol = time_int_tol,

69
            ode_default_options()..., callback = callbacks)

70


71