libtrixi/LibTrixi.jl/examples/libelixir_p4est2d_dgsem_euler_sedov.jl at 223866de090dddfcf2f328008b1cf891d4fc1b4a · trixi-framework/libtrixi · GitHub

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
using LibTrixi
using Trixi
using OrdinaryDiffEq

# The function to create the simulation state needs to be named `init_simstate`
function init_simstate()

    ###############################################################################
    # semidiscretization of the compressible Euler equations

    equations = CompressibleEulerEquations2D(1.4)

    """
        initial_condition_sedov_blast_wave(x, t, equations::CompressibleEulerEquations2D)

    The Sedov blast wave setup based on Flash
    - http://flash.uchicago.edu/site/flashcode/user_support/flash_ug_devel/node184.html#SECTION010114000000000000000
    """
    function initial_condition_sedov_blast_wave(x, t, equations::CompressibleEulerEquations2D)
    # Set up polar coordinates
    inicenter = SVector(0.0, 0.0)
    x_norm = x[1] - inicenter[1]
    y_norm = x[2] - inicenter[2]
    r = sqrt(x_norm^2 + y_norm^2)

    # Setup based on http://flash.uchicago.edu/site/flashcode/user_support/flash_ug_devel/node184.html#SECTION010114000000000000000
    r0 = 0.21875 # = 3.5 * smallest dx (for domain length=4 and max-ref=6)
    E = 1.0
    p0_inner = 3 * (equations.gamma - 1) * E / (3 * pi * r0^2)
    p0_outer = 1.0e-5 # = true Sedov setup

    # Calculate primitive variables
    rho = 1.0
    v1  = 0.0
    v2  = 0.0
    p   = r > r0 ? p0_outer : p0_inner

    return prim2cons(SVector(rho, v1, v2, p), equations)
    end

    initial_condition = initial_condition_sedov_blast_wave

    # Get the DG approximation space
    surface_flux = flux_lax_friedrichs
    volume_flux = flux_ranocha
    polydeg = 4
    basis = LobattoLegendreBasis(polydeg)
    indicator_sc = IndicatorHennemannGassner(equations, basis,
                                            alpha_max=1.0,
                                            alpha_min=0.001,
                                            alpha_smooth=true,
                                            variable=density_pressure)
    volume_integral = VolumeIntegralShockCapturingHG(indicator_sc;
                                                    volume_flux_dg=volume_flux,
                                                    volume_flux_fv=surface_flux)

    solver = DGSEM(polydeg=polydeg, surface_flux=surface_flux, volume_integral=volume_integral)

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

    coordinates_min = (-1.0, -1.0)
    coordinates_max = ( 1.0,  1.0)

    trees_per_dimension = (4, 4)
    mesh = P4estMesh(trees_per_dimension,
                    polydeg=4, initial_refinement_level=2,
                    coordinates_min=coordinates_min, coordinates_max=coordinates_max,
                    periodicity=true)

    semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver)

    ###############################################################################
    # ODE solvers, callbacks etc.

    # tspan = (0.0, 12.5) # original timespan
    tspan = (0.0, 2.5)
    ode = semidiscretize(semi, tspan)

    summary_callback = SummaryCallback()

    analysis_interval = 300
    analysis_callback = AnalysisCallback(semi, interval=analysis_interval)

    alive_callback = AliveCallback(analysis_interval=analysis_interval)

    save_solution = SaveSolutionCallback(interval=300,
                                        save_initial_solution=true,
                                        save_final_solution=true)

    stepsize_callback = StepsizeCallback(cfl=0.5)

    callbacks = CallbackSet(summary_callback,
                            analysis_callback,
                            alive_callback,
                            save_solution,
                            stepsize_callback)

    ###############################################################################
    # create the time integrator

    # OrdinaryDiffEq's `integrator`
    integrator = init(ode, CarpenterKennedy2N54(williamson_condition=false),
                      dt=1.0, # solve needs some value here but it will be overwritten by the stepsize_callback ?!
                      save_everystep=false, callback=callbacks);

    ###############################################################################
    # Create simulation state

    # registry only used for tests
    registry = LibTrixiDataRegistry(undef, 1)
    simstate = SimulationState(semi, integrator, registry)

    return simstate
end