WIP: Add FV method for T8codeMesh

.github/workflows/ci.yml

@@ -71,6 +71,7 @@ jobs:
           - p4est_part2
           - t8code_part1
           - t8code_part2
+          - t8code_fv
           - unstructured_dgmulti
           - parabolic
           - paper_self_gravitating_gas_dynamics

examples/t8code_2d_fv/elixir_advection_basic.jl

@@ -0,0 +1,43 @@
+using OrdinaryDiffEq
+using Trixi
+
+advection_velocity = (0.2, -0.7)
+equations = LinearScalarAdvectionEquation2D(advection_velocity)
+
+initial_condition = initial_condition_convergence_test
+
+solver = FV(surface_flux = flux_lax_friedrichs)
+
+mesh_function = Trixi.cmesh_new_periodic_hybrid
+# mesh_function = Trixi.cmesh_new_periodic_quad
+# mesh_function = Trixi.cmesh_new_periodic_tri
+cmesh = mesh_function(Trixi.mpi_comm())
+mesh = T8codeMesh(cmesh, solver,
+                  initial_refinement_level = 3)
+
+semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver)
+
+ode = semidiscretize(semi, (0.0, 1.0));
+
+summary_callback = SummaryCallback()
+
+analysis_interval = 100
+analysis_callback = AnalysisCallback(semi, interval = analysis_interval)
+
+alive_callback = AliveCallback(analysis_interval = analysis_interval)
+
+save_solution = SaveSolutionCallback(interval = 10,
+                                     solution_variables = cons2prim)
+
+stepsize_callback = StepsizeCallback(cfl = 0.5)
+
+callbacks = CallbackSet(summary_callback, analysis_callback, alive_callback,
+                        stepsize_callback, save_solution)
+
+###############################################################################
+# run the simulation
+
+sol = solve(ode, CarpenterKennedy2N54(williamson_condition = false),
+            dt = 5.0e-2, # solve needs some value here but it will be overwritten by the stepsize_callback
+            save_everystep = false, callback = callbacks);
+summary_callback()

examples/t8code_2d_fv/elixir_advection_gauss.jl

@@ -0,0 +1,42 @@
+using OrdinaryDiffEq
+using Trixi
+
+####################################################
+
+advection_velocity = (0.2, -0.7)
+equations = LinearScalarAdvectionEquation2D(advection_velocity)
+
+initial_condition = initial_condition_gauss
+
+solver = FV(surface_flux = flux_lax_friedrichs)
+
+initial_refinement_level = 4
+cmesh = Trixi.cmesh_new_periodic_hybrid(Trixi.mpi_comm())
+mesh = T8codeMesh(cmesh, solver, initial_refinement_level = initial_refinement_level)
+
+semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver)
+
+ode = semidiscretize(semi, (0.0, 20.0));
+
+summary_callback = SummaryCallback()
+
+analysis_interval = 100
+analysis_callback = AnalysisCallback(semi, interval = analysis_interval)
+
+alive_callback = AliveCallback(analysis_interval = analysis_interval)
+
+save_solution = SaveSolutionCallback(interval = 10,
+                                     solution_variables = cons2prim)
+
+stepsize_callback = StepsizeCallback(cfl = 0.5)
+
+callbacks = CallbackSet(summary_callback, save_solution, analysis_callback, alive_callback,
+                        stepsize_callback)
+
+###############################################################################
+# run the simulation
+
+sol = solve(ode, CarpenterKennedy2N54(williamson_condition = false),#Euler(),
+            dt = 5.0e-2, # solve needs some value here but it will be overwritten by the stepsize_callback
+            save_everystep = false, callback = callbacks);
+summary_callback()

examples/t8code_2d_fv/elixir_euler_blast_wave.jl

@@ -0,0 +1,59 @@
+using OrdinaryDiffEq
+using Trixi
+
+####################################################
+
+equations = CompressibleEulerEquations2D(1.4)
+
+function initial_condition_blast_wave(x, t, equations::CompressibleEulerEquations2D)
+    # Modified From Hennemann & Gassner JCP paper 2020 (Sec. 6.3) -> "medium blast wave"
+    # 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)
+    phi = atan(y_norm, x_norm)
+    sin_phi, cos_phi = sincos(phi)
+
+    # Calculate primitive variables
+    rho = r > 0.5 ? 1.0 : 1.1691
+    v1 = r > 0.5 ? 0.0 : 0.1882 * cos_phi
+    v2 = r > 0.5 ? 0.0 : 0.1882 * sin_phi
+    p = r > 0.5 ? 1.0E-3 : 1.245
+
+    return prim2cons(SVector(rho, v1, v2, p), equations)
+end
+initial_condition = initial_condition_blast_wave
+
+solver = FV(surface_flux = flux_lax_friedrichs)
+
+initial_refinement_level = 4
+cmesh = Trixi.cmesh_new_periodic_hybrid2(Trixi.mpi_comm())
+mesh = T8codeMesh(cmesh, solver, initial_refinement_level = initial_refinement_level)
+
+semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver)
+
+ode = semidiscretize(semi, (0.0, 2.0));
+
+summary_callback = SummaryCallback()
+
+analysis_interval = 100
+analysis_callback = AnalysisCallback(semi, interval = analysis_interval)
+
+alive_callback = AliveCallback(analysis_interval = analysis_interval)
+
+save_solution = SaveSolutionCallback(interval = 10,
+                                     solution_variables = cons2prim)
+
+stepsize_callback = StepsizeCallback(cfl = 0.5)
+
+callbacks = CallbackSet(summary_callback, save_solution, analysis_callback, alive_callback,
+                        stepsize_callback)
+
+###############################################################################
+# run the simulation
+
+sol = solve(ode, CarpenterKennedy2N54(williamson_condition = false),
+            dt = 5.0e-2, # solve needs some value here but it will be overwritten by the stepsize_callback
+            save_everystep = false, callback = callbacks);
+summary_callback()

examples/t8code_2d_fv/elixir_euler_kelvin_helmholtz_instability.jl

@@ -0,0 +1,74 @@
+using OrdinaryDiffEq
+using Trixi
+
+###############################################################################
+# semidiscretization of the compressible Euler equations
+gamma = 1.4
+equations = CompressibleEulerEquations2D(gamma)
+
+"""
+    initial_condition_kelvin_helmholtz_instability(x, t, equations::CompressibleEulerEquations2D)
+
+A version of the classical Kelvin-Helmholtz instability based on
+- Andrés M. Rueda-Ramírez, Gregor J. Gassner (2021)
+  A Subcell Finite Volume Positivity-Preserving Limiter for DGSEM Discretizations
+  of the Euler Equations
+  [arXiv: 2102.06017](https://arxiv.org/abs/2102.06017)
+"""
+function initial_condition_kelvin_helmholtz_instability(x, t,
+                                                        equations::CompressibleEulerEquations2D)
+    # change discontinuity to tanh
+    # typical resolution 128^2, 256^2
+    # domain size is [-1,+1]^2
+    slope = 15
+    amplitude = 0.02
+    B = tanh(slope * x[2] + 7.5) - tanh(slope * x[2] - 7.5)
+    rho = 0.5 + 0.75 * B
+    v1 = 0.5 * (B - 1)
+    v2 = 0.1 * sin(2 * pi * x[1])
+    p = 1.0
+    return prim2cons(SVector(rho, v1, v2, p), equations)
+end
+initial_condition = initial_condition_kelvin_helmholtz_instability
+
+solver = FV(surface_flux = flux_lax_friedrichs)
+
+initial_refinement_level = 4
+# mesh_function = Trixi.cmesh_new_periodic_hybrid
+# mesh_function = Trixi.cmesh_new_periodic_quad
+# mesh_function = Trixi.cmesh_new_periodic_tri
+mesh_function = Trixi.cmesh_new_periodic_tri2
+cmesh = mesh_function(Trixi.mpi_comm())
+mesh = T8codeMesh(cmesh, solver, initial_refinement_level = initial_refinement_level)
+
+semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver)
+
+tspan = (0.0, 3.7)
+ode = semidiscretize(semi, tspan)
+
+summary_callback = SummaryCallback()
+
+analysis_interval = 100
+analysis_callback = AnalysisCallback(semi, interval = analysis_interval)
+
+alive_callback = AliveCallback(analysis_interval = analysis_interval)
+
+save_solution = SaveSolutionCallback(interval = 40,
+                                     save_initial_solution = true,
+                                     save_final_solution = true,
+                                     solution_variables = cons2prim)
+
+stepsize_callback = StepsizeCallback(cfl = 0.9)
+
+callbacks = CallbackSet(summary_callback,
+                        analysis_callback, alive_callback,
+                        save_solution,
+                        stepsize_callback)
+
+###############################################################################
+# run the simulation
+
+sol = solve(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);
+summary_callback() # print the timer summary

src/Trixi.jl

@@ -242,6 +242,8 @@ export DG,
 export VolumeIntegralSubcellLimiting, BoundsCheckCallback,
        SubcellLimiterIDP, SubcellLimiterIDPCorrection
 
+export FV
+
 export nelements, nnodes, nvariables,
        eachelement, eachnode, eachvariable
 

src/callbacks_step/analysis.jl

@@ -26,8 +26,8 @@ or `extra_analysis_errors = (:conservation_error,)`.
 If you want to omit the computation (to safe compute-time) of the [`default_analysis_errors`](@ref), specify
 `analysis_errors = Symbol[]`.
 Note: `default_analysis_errors` are `:l2_error` and `:linf_error` for all equations.
-If you want to compute `extra_analysis_errors` such as `:conservation_error` solely, i.e., 
-without `:l2_error, :linf_error` you need to specify 
+If you want to compute `extra_analysis_errors` such as `:conservation_error` solely, i.e.,
+without `:l2_error, :linf_error` you need to specify
 `analysis_errors = [:conservation_error]` instead of `extra_analysis_errors = [:conservation_error]`.
 
 Further scalar functions `func` in `extra_analysis_integrals` are applied to the numerical
@@ -141,6 +141,49 @@ function AnalysisCallback(mesh, equations::AbstractEquations, solver, cache;
                      initialize = initialize!)
 end
 
+function AnalysisCallback(mesh, equations::AbstractEquations, solver::FV, cache;
+                          interval = 0,
+                          save_analysis = false,
+                          output_directory = "out",
+                          analysis_filename = "analysis.dat",
+                          extra_analysis_errors = Symbol[],
+                          analysis_errors = union(default_analysis_errors(equations),
+                                                  extra_analysis_errors),
+                          extra_analysis_integrals = (),
+                          analysis_integrals = union(default_analysis_integrals(equations),
+                                                     extra_analysis_integrals),
+                          RealT = real(solver),
+                          uEltype = eltype(cache.elements[1]), # TODO: this routine is equal to the existing one except this default type here.
+                          kwargs...)
+    # Decide when the callback is activated.
+    # With error-based step size control, some steps can be rejected. Thus,
+    #   `integrator.iter >= integrator.stats.naccept`
+    #    (total #steps)       (#accepted steps)
+    # We need to check the number of accepted steps since callbacks are not
+    # activated after a rejected step.
+    condition = (u, t, integrator) -> interval > 0 &&
+        ((integrator.stats.naccept % interval == 0 &&
+          !(integrator.stats.naccept == 0 && integrator.iter > 0)) ||
+         isfinished(integrator))
+
+    analyzer = SolutionAnalyzer(solver; kwargs...)
+    cache_analysis = create_cache_analysis(analyzer, mesh, equations, solver, cache,
+                                           RealT, uEltype)
+
+    analysis_callback = AnalysisCallback(0.0, 0.0, 0, 0.0,
+                                         interval, save_analysis, output_directory,
+                                         analysis_filename,
+                                         analyzer,
+                                         analysis_errors, Tuple(analysis_integrals),
+                                         SVector(ntuple(_ -> zero(uEltype),
+                                                        Val(nvariables(equations)))),
+                                         cache_analysis)
+
+    DiscreteCallback(condition, analysis_callback,
+                     save_positions = (false, false),
+                     initialize = initialize!)
+end
+
 # This method gets called from OrdinaryDiffEq's `solve(...)`
 function initialize!(cb::DiscreteCallback{Condition, Affect!}, u_ode, t,
                      integrator) where {Condition, Affect! <: AnalysisCallback}