First-order Finite Volume Method organized with t8code

examples/t8code_2d_fv/elixir_advection_basic.jl

@@ -0,0 +1,52 @@
+#!/usr/bin/env julia
+# Has to be at least Julia v1.9.0.
+
+# Running in parallel:
+#   ${JULIA_DEPOT_PATH}/.julia/bin/mpiexecjl --project=. -n 3 julia hybrid-t8code-mesh.jl
+#
+# More information: https://juliaparallel.org/MPI.jl/stable/usage/
+using OrdinaryDiffEq
+using Trixi
+
+advection_velocity = (0.2, -0.7)
+equations = LinearScalarAdvectionEquation2D(advection_velocity)
+
+initial_condition = initial_condition_convergence_test
+
+solver = FV(order = 1, surface_flux = flux_lax_friedrichs)
+
+coordinates_min = (-1.0, -1.0) # minimum coordinates (min(x), min(y))
+coordinates_max = (1.0, 1.0) # maximum coordinates (max(x), max(y))
+
+trees_per_dimension = (8, 8)
+
+mesh = T8codeMesh(trees_per_dimension, polydeg = 3,
+                  coordinates_min = coordinates_min, coordinates_max = coordinates_max,
+                  initial_refinement_level = 1)
+
+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 = 1,
+                                     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_advection_convergence.jl

@@ -0,0 +1,49 @@
+#!/usr/bin/env julia
+# Has to be at least Julia v1.9.0.
+
+# Running in parallel:
+#   ${JULIA_DEPOT_PATH}/.julia/bin/mpiexecjl --project=. -n 3 julia hybrid-t8code-mesh.jl
+#
+# More information: https://juliaparallel.org/MPI.jl/stable/usage/
+using OrdinaryDiffEq
+using Trixi
+
+advection_velocity = (0.2, -0.7)
+equations = LinearScalarAdvectionEquation2D(advection_velocity)
+
+initial_condition = initial_condition_convergence_test
+
+solver = FV(order = 1, 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 = T8codeFVMesh{2}(mesh_function, initial_refinement_level)
+
+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 = 1,
+                                     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_advection_gauss.jl

@@ -0,0 +1,48 @@
+#!/usr/bin/env julia
+# Has to be at least Julia v1.9.0.
+
+# Running in parallel:
+#   ${JULIA_DEPOT_PATH}/.julia/bin/mpiexecjl --project=. -n 3 julia hybrid-t8code-mesh.jl
+#
+# More information: https://juliaparallel.org/MPI.jl/stable/usage/
+using OrdinaryDiffEq
+using Trixi
+
+####################################################
+
+advection_velocity = (0.2, -0.7)
+equations = LinearScalarAdvectionEquation2D(advection_velocity)
+
+initial_condition = initial_condition_gauss
+
+solver = FV(order = 1, surface_flux = flux_lax_friedrichs)
+
+initial_refinement_level = 4
+mesh = T8codeFVMesh{2}(Trixi.cmesh_new_periodic_hybrid, 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 = 1,
+                                     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_hybrid_mesh.jl

@@ -0,0 +1,65 @@
+#!/usr/bin/env julia
+# Has to be at least Julia v1.9.0.
+
+# Running in parallel:
+#   ${JULIA_DEPOT_PATH}/.julia/bin/mpiexecjl --project=. -n 3 julia hybrid-t8code-mesh.jl
+#
+# More information: https://juliaparallel.org/MPI.jl/stable/usage/
+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(order = 1, surface_flux = flux_lax_friedrichs)
+
+initial_refinement_level = 4
+mesh = T8codeFVMesh{2}(Trixi.cmesh_new_periodic_hybrid2, 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,80 @@
+#!/usr/bin/env julia
+# Has to be at least Julia v1.9.0.
+
+# Running in parallel:
+#   ${JULIA_DEPOT_PATH}/.julia/bin/mpiexecjl --project=. -n 3 julia hybrid-t8code-mesh.jl
+#
+# More information: https://juliaparallel.org/MPI.jl/stable/usage/
+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(order = 1, 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
+mesh = T8codeFVMesh{2}(mesh_function, 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

@@ -221,7 +221,7 @@ export entropy, energy_total, energy_kinetic, energy_internal, energy_magnetic,
 export lake_at_rest_error
 export ncomponents, eachcomponent
 
-export TreeMesh, StructuredMesh, UnstructuredMesh2D, P4estMesh, T8codeMesh
+export TreeMesh, StructuredMesh, UnstructuredMesh2D, P4estMesh, T8codeMesh, T8codeFVMesh
 
 export DG,
        DGSEM, LobattoLegendreBasis,
@@ -240,6 +240,8 @@ export DG,
 export VolumeIntegralSubcellLimiting, BoundsCheckCallback,
        SubcellLimiterIDP, SubcellLimiterIDPCorrection
 
+export FV
+
 export nelements, nnodes, nvariables,
        eachelement, eachnode, eachvariable
 

src/callbacks_step/analysis.jl

@@ -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 = RealT,#TODO:eltype(cache.elements),
+                          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}