libtrixi/LibTrixi.jl/examples/libelixir_t8code3d_euler_tracer.jl at 41415f14150203cebba280dbb203658c704daba0 · 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
115
116
117
118
119
120
121
122
123
124
# A manufactured solution of a circular wind with constant angular velocity
# on a planetary-sized cubed sphere mesh with a blob detected by AMR
#
# Note that this libelixir is based on an elixir by Erik Faulhaber for Trixi.jl
# Source: https://github.com/trixi-framework/Trixi.jl/blob/main/examples/p4est_3d_dgsem/elixir_euler_circular_wind_nonconforming.jl

using OrdinaryDiffEqLowStorageRK
using Trixi
using LinearAlgebra
using LibTrixi


function initial_condition_circular_wind(x, t, equations::CompressibleEulerEquations3D)
    radius_earth = 6.371229e6
    lambda, phi, r = cart_to_sphere(x)

    p = 1e5
    v1 = -10 * x[2] / radius_earth
    v2 = 10 * x[1] / radius_earth
    v3 = 0.0
    rho = 1.0 + 0.1  * exp(-50 * ((lambda-1.0)^2/2.0 + (phi-0.4)^2)) +
                0.08 * exp(-100 * ((lambda-0.8)^2/4.0 + (phi-0.5)^2))

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

@inline function source_terms_circular_wind(u, x, t,
                                            equations::CompressibleEulerEquations3D)
    radius_earth = 6.371229e6
    rho = u[1]

    du1 = 0.0
    du2 = -rho * (10 / radius_earth) * (10 * x[1] / radius_earth)
    du3 = -rho * (10 / radius_earth) * (10 * x[2] / radius_earth)
    du4 = 0.0
    du5 = 0.0

    return SVector(du1, du2, du3, du4, du5)
end

function cart_to_sphere(x)
    r = norm(x)
    lambda = atan(x[2], x[1])
    if lambda < 0
        lambda += 2 * pi
    end
    phi = asin(x[3] / r)

    return lambda, phi, r
end


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

    # compressible Euler equations
    gamma = 1.4
    equations = CompressibleEulerEquations3D(gamma)

    # setup of the problem
    initial_condition = initial_condition_circular_wind

    boundary_conditions = (; inside = boundary_condition_slip_wall,
                           outside = boundary_condition_slip_wall)

    # estimate for speed of sound
    surface_flux = FluxLMARS(374)
    solver = DGSEM(polydeg = 3, surface_flux = surface_flux)

    # increase trees_per_cube_face to 4 to get nicer results
    lat_lon_levels = 2
    layers = 1
    mesh = Trixi.T8codeMeshCubedSphere(lat_lon_levels, layers, 6.371229e6, 30000.0,
                                       polydeg = 3, initial_refinement_level = 0)

    semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver;
                                        source_terms = source_terms_circular_wind,
                                        boundary_conditions = boundary_conditions)

    # increase number of days
    days = 0.1
    tspan = (0.0, days * 24 * 60 * 60.0)

    ode = semidiscretize(semi, tspan)

    summary_callback = SummaryCallback()

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

    alive_callback = AliveCallback(analysis_interval = analysis_interval)

    save_solution = SaveSolutionCallback(interval = 2000,
                                         save_initial_solution = true,
                                         save_final_solution = true,
                                         solution_variables = cons2prim,
                                         output_directory = "out_tracer")

    amr_controller = ControllerThreeLevel(semi, IndicatorMax(semi, variable = first),
                                          base_level = 0,
                                          med_level = 1, med_threshold = 1.004,
                                          max_level = 3, max_threshold = 1.11)

    amr_callback = AMRCallback(semi, amr_controller,
                               interval = 2000,
                               adapt_initial_condition = true,
                               adapt_initial_condition_only_refine = true)

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

    # use a Runge-Kutta method with automatic (error based) time step size control
    integrator = init(ode, RDPK3SpFSAL49(thread = Trixi.False());
                      abstol = 1.0e-6, reltol = 1.0e-6,
                      ode_default_options()..., callback = callbacks, maxiters=1e7);

    # create simulation state
    simstate = SimulationState(semi, integrator)

    return simstate
end