|
| 1 | + |
| 2 | +using OrdinaryDiffEqSSPRK, OrdinaryDiffEqLowStorageRK |
| 3 | +using Trixi |
| 4 | +using TrixiShallowWater |
| 5 | + |
| 6 | +############################################################################### |
| 7 | +# semidiscretization of the shallow water equations with a continuous |
| 8 | +# bottom topography function and a perturbation in the water height |
| 9 | + |
| 10 | +equations = ShallowWaterEquationsWetDry2D(gravity = 9.812, H0 = 2.1) |
| 11 | + |
| 12 | +function initial_condition_perturbation(x, t, equations::ShallowWaterEquationsWetDry2D) |
| 13 | + # Calculate primitive variables |
| 14 | + H = equations.H0 |
| 15 | + v1 = 0.0 |
| 16 | + v2 = 0.0 |
| 17 | + |
| 18 | + x1, x2 = x |
| 19 | + b = (1.75 / exp(0.5 * ((x1 - 1.0)^2 + (x2 + 1.0)^2)) |
| 20 | + + |
| 21 | + 0.8 / exp(0.5 * ((x1 + 1.0)^2 + (x2 - 1.0)^2)) |
| 22 | + - |
| 23 | + 0.5 / exp(3.5 * ((x1 - 0.4)^2 + (x2 - 0.325)^2))) |
| 24 | + |
| 25 | + # Waterheight perturbation |
| 26 | + H = H + 0.5 * exp(-40.0 * ((x[1])^2 + (x[2])^2)) |
| 27 | + |
| 28 | + return prim2cons(SVector(H, v1, v2, b), equations) |
| 29 | +end |
| 30 | + |
| 31 | +initial_condition = initial_condition_perturbation |
| 32 | + |
| 33 | +boundary_condition = Dict(:x_neg => boundary_condition_slip_wall, |
| 34 | + :y_neg => boundary_condition_slip_wall, |
| 35 | + :x_pos => boundary_condition_slip_wall, |
| 36 | + :y_pos => boundary_condition_slip_wall) |
| 37 | + |
| 38 | +############################################################################### |
| 39 | +# Get the DG approximation space |
| 40 | + |
| 41 | +volume_flux = (flux_wintermeyer_etal, flux_nonconservative_wintermeyer_etal) |
| 42 | + |
| 43 | +surface_flux = (FluxHydrostaticReconstruction(flux_hll_chen_noelle, |
| 44 | + hydrostatic_reconstruction_chen_noelle), |
| 45 | + flux_nonconservative_chen_noelle) |
| 46 | + |
| 47 | +# Create the solver |
| 48 | +solver = DGSEM(polydeg = 3, surface_flux = surface_flux, |
| 49 | + volume_integral = VolumeIntegralFluxDifferencing(volume_flux)) |
| 50 | + |
| 51 | +############################################################################### |
| 52 | +# This setup is for the curved, split form well-balancedness testing |
| 53 | + |
| 54 | +# Affine type mapping to take the [-1,1]^2 domain from the mesh file |
| 55 | +# and warp it as described in https://arxiv.org/abs/2012.12040 |
| 56 | +function mapping_twist(xi, eta) |
| 57 | + y = eta + 0.2 * cos(1.5 * pi * xi) * cos(0.5 * pi * eta) |
| 58 | + x = xi + 0.2 * cos(0.5 * pi * xi) * cos(2.0 * pi * y) |
| 59 | + return SVector(x, y) |
| 60 | +end |
| 61 | + |
| 62 | +# The mesh below can be made periodic if desired |
| 63 | + |
| 64 | +# Create P4estMesh with 4 x 4 trees |
| 65 | +trees_per_dimension = (4, 4) |
| 66 | +mesh = P4estMesh(trees_per_dimension, polydeg = 3, |
| 67 | + mapping = mapping_twist, |
| 68 | + initial_refinement_level = 0, |
| 69 | + periodicity = false) |
| 70 | + |
| 71 | +# Refine bottom left quadrant of each tree to level 3 |
| 72 | +function refine_fn(p4est, which_tree, quadrant) |
| 73 | + quadrant_obj = unsafe_load(quadrant) |
| 74 | + if quadrant_obj.x == 0 && quadrant_obj.y == 0 && quadrant_obj.level < 3 |
| 75 | + # return true (refine) |
| 76 | + return Cint(1) |
| 77 | + else |
| 78 | + # return false (don't refine) |
| 79 | + return Cint(0) |
| 80 | + end |
| 81 | +end |
| 82 | + |
| 83 | +# Refine recursively until each bottom left quadrant of a tree has level 3 |
| 84 | +# The mesh will be rebalanced before the simulation starts |
| 85 | +refine_fn_c = @cfunction(refine_fn, Cint, |
| 86 | + (Ptr{Trixi.p4est_t}, Ptr{Trixi.p4est_topidx_t}, |
| 87 | + Ptr{Trixi.p4est_quadrant_t})) |
| 88 | +Trixi.refine_p4est!(mesh.p4est, true, refine_fn_c, C_NULL) |
| 89 | + |
| 90 | +# Create the semi discretization object |
| 91 | +semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver, |
| 92 | + boundary_conditions = boundary_condition) |
| 93 | +############################################################################### |
| 94 | +# ODE solvers, callbacks, etc. |
| 95 | + |
| 96 | +tspan = (0.0, 3.0) |
| 97 | +ode = semidiscretize(semi, tspan) |
| 98 | + |
| 99 | +############################################################################### |
| 100 | + |
| 101 | +summary_callback = SummaryCallback() |
| 102 | + |
| 103 | +analysis_interval = 1000 |
| 104 | +analysis_callback = AnalysisCallback(semi, interval = analysis_interval, |
| 105 | + extra_analysis_errors = (:conservation_error,)) |
| 106 | + |
| 107 | +alive_callback = AliveCallback(analysis_interval = analysis_interval) |
| 108 | + |
| 109 | +save_solution = SaveSolutionCallback(dt = 0.04, |
| 110 | + save_initial_solution = true, |
| 111 | + save_final_solution = true) |
| 112 | + |
| 113 | +# Define the perturbation of water height as a variable to use in the AMR indicator |
| 114 | +@inline function waterheight_total(u, equations::ShallowWaterEquationsWetDry2D) |
| 115 | + return u[1] + u[4] |
| 116 | +end |
| 117 | + |
| 118 | +amr_controller = ControllerThreeLevel(semi, |
| 119 | + IndicatorMax(semi, variable = waterheight_total), |
| 120 | + base_level = 1, |
| 121 | + med_level = 2, med_threshold = 2.01, |
| 122 | + max_level = 4, max_threshold = 2.15) |
| 123 | + |
| 124 | +amr_callback = AMRCallback(semi, amr_controller, |
| 125 | + interval = 1, |
| 126 | + adapt_initial_condition = false, |
| 127 | + adapt_initial_condition_only_refine = false) |
| 128 | + |
| 129 | +stepsize_callback = StepsizeCallback(cfl = 1.0) |
| 130 | + |
| 131 | +callbacks = CallbackSet(summary_callback, |
| 132 | + analysis_callback, |
| 133 | + alive_callback, |
| 134 | + save_solution, |
| 135 | + amr_callback, |
| 136 | + stepsize_callback) |
| 137 | + |
| 138 | +############################################################################### |
| 139 | +# run the simulation |
| 140 | +sol = solve(ode, CarpenterKennedy2N54(williamson_condition = false), |
| 141 | + dt = 1.0, # solve needs some value here but it will be overwritten by the stepsize_callback |
| 142 | + adaptive = false, |
| 143 | + save_everystep = false, callback = callbacks); |
0 commit comments