[WIP] Three equations model for sloshing simulations

Project.toml

@@ -99,7 +99,7 @@ StaticArrays = "1.5"
 StrideArrays = "0.1.26"
 StructArrays = "0.6.11"
 SummationByPartsOperators = "0.5.41"
-T8code = "0.5"
+T8code = "0.6"
 TimerOutputs = "0.5.7"
 Triangulate = "2.2"
 TriplotBase = "0.1"

examples/t8code_2d_dgsem/elixir_advection_amr_solution_independent.jl

@@ -115,6 +115,11 @@ analysis_callback = AnalysisCallback(semi, interval = analysis_interval,
 
 alive_callback = AliveCallback(analysis_interval = analysis_interval)
 
+save_solution = SaveSolutionCallback(interval = 100,
+                                     save_initial_solution = true,
+                                     save_final_solution = true,
+                                     solution_variables = cons2prim)
+
 amr_controller = ControllerThreeLevel(semi,
                                       TrixiExtension.IndicatorSolutionIndependent(semi),
                                       base_level = 4,
@@ -124,12 +129,20 @@ amr_controller = ControllerThreeLevel(semi,
 amr_callback = AMRCallback(semi, amr_controller,
                            interval = 5,
                            adapt_initial_condition = true,
-                           adapt_initial_condition_only_refine = true)
+                           adapt_initial_condition_only_refine = true,
+                           dynamic_load_balancing = false)
+# We disable `dynamic_load_balancing` for now, since t8code does not support
+# partitioning for coarsening yet. That is, a complete family of elements always
+# stays on rank and is not split up due to partitioning. Without this feature
+# dynamic AMR simulations are not pefectly deterministic regarding to
+# convergent tests. Once this feature is available in t8code load balancing is
+# enabled again.
 
 stepsize_callback = StepsizeCallback(cfl = 1.6)
 
 callbacks = CallbackSet(summary_callback,
                         analysis_callback, alive_callback,
+                        save_solution,
                         amr_callback, stepsize_callback);
 
 ###############################################################################
@@ -139,3 +152,7 @@ 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
+
+# Finalize `T8codeMesh` to make sure MPI related objects in t8code are
+# released before `MPI` finalizes.
+!isinteractive() && finalize(mesh)

examples/t8code_2d_dgsem/elixir_advection_amr_unstructured_flag.jl

@@ -54,19 +54,35 @@ analysis_callback = AnalysisCallback(semi, interval = analysis_interval,
 
 alive_callback = AliveCallback(analysis_interval = analysis_interval)
 
+save_restart = SaveRestartCallback(interval = 100,
+                                   save_final_restart = true)
+
+save_solution = SaveSolutionCallback(interval = 100,
+                                     save_initial_solution = true,
+                                     save_final_solution = true,
+                                     solution_variables = cons2prim)
+
 amr_controller = ControllerThreeLevel(semi, IndicatorMax(semi, variable = first),
                                       base_level = 1,
                                       med_level = 2, med_threshold = 0.1,
                                       max_level = 3, max_threshold = 0.6)
 amr_callback = AMRCallback(semi, amr_controller,
                            interval = 5,
                            adapt_initial_condition = true,
-                           adapt_initial_condition_only_refine = true)
+                           adapt_initial_condition_only_refine = true,
+                           dynamic_load_balancing = false)
+# We disable `dynamic_load_balancing` for now, since t8code does not support
+# partitioning for coarsening yet. That is, a complete family of elements always
+# stays on rank and is not split up due to partitioning. Without this feature
+# dynamic AMR simulations are not pefectly deterministic regarding to
+# convergent tests. Once this feature is available in t8code load balancing is
+# enabled again.
 
 stepsize_callback = StepsizeCallback(cfl = 0.7)
 
 callbacks = CallbackSet(summary_callback,
                         analysis_callback, alive_callback,
+                        save_restart, save_solution,
                         amr_callback, stepsize_callback)
 
 ###############################################################################
@@ -77,3 +93,7 @@ sol = solve(ode, CarpenterKennedy2N54(williamson_condition = false),
             save_everystep = false, callback = callbacks);
 
 summary_callback() # print the timer summary
+
+# Finalize `T8codeMesh` to make sure MPI related objects in t8code are
+# released before `MPI` finalizes.
+!isinteractive() && finalize(mesh)

examples/t8code_2d_dgsem/elixir_advection_basic.jl

@@ -55,3 +55,7 @@ sol = solve(ode, CarpenterKennedy2N54(williamson_condition = false),
 
 # Print the timer summary
 summary_callback()
+
+# Finalize `T8codeMesh` to make sure MPI related objects in t8code are
+# released before `MPI` finalizes.
+!isinteractive() && finalize(mesh)

examples/t8code_2d_dgsem/elixir_advection_extended.jl

@@ -0,0 +1,89 @@
+
+using OrdinaryDiffEq
+using Trixi
+
+###############################################################################
+# semidiscretization of the linear advection equation
+
+advection_velocity = (0.2, -0.7)
+equations = LinearScalarAdvectionEquation2D(advection_velocity)
+
+initial_condition = initial_condition_convergence_test
+
+# BCs must be passed as Dict
+boundary_condition = BoundaryConditionDirichlet(initial_condition)
+boundary_conditions = Dict(:x_neg => boundary_condition,
+                           :x_pos => boundary_condition,
+                           :y_neg => boundary_condition,
+                           :y_pos => boundary_condition)
+
+# Create DG solver with polynomial degree = 3 and (local) Lax-Friedrichs/Rusanov flux as surface flux
+solver = DGSEM(polydeg = 3, surface_flux = flux_lax_friedrichs)
+
+# The initial condition is 2-periodic
+coordinates_min = (-1.5, 1.3) # minimum coordinates (min(x), min(y))
+coordinates_max = (0.5, 5.3) # maximum coordinates (max(x), max(y))
+
+trees_per_dimension = (19, 37)
+
+# Create curved mesh with 19 x 37 elements
+mesh = T8codeMesh(trees_per_dimension, polydeg = 3,
+                  coordinates_min = coordinates_min, coordinates_max = coordinates_max,
+                  periodicity = false)
+
+# A semidiscretization collects data structures and functions for the spatial discretization
+semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver,
+                                    boundary_conditions = boundary_conditions)
+
+###############################################################################
+# ODE solvers, callbacks etc.
+
+# Create ODE problem with time span from 0.0 to 1.0
+tspan = (0.0, 1.0)
+ode = semidiscretize(semi, tspan);
+
+# At the beginning of the main loop, the SummaryCallback prints a summary of the simulation setup
+# and resets the timers
+summary_callback = SummaryCallback()
+
+# The AnalysisCallback allows to analyse the solution in regular intervals and prints the results
+analysis_interval = 100
+analysis_callback = AnalysisCallback(semi, interval = analysis_interval,
+                                     extra_analysis_integrals = (entropy, energy_total))
+
+# The AliveCallback prints short status information in regular intervals
+alive_callback = AliveCallback(analysis_interval = analysis_interval)
+
+# The SaveRestartCallback allows to save a file from which a Trixi.jl simulation can be restarted
+save_restart = SaveRestartCallback(interval = 100,
+                                   save_final_restart = true)
+
+# The SaveSolutionCallback allows to save the solution to a file in regular intervals
+save_solution = SaveSolutionCallback(interval = 100,
+                                     save_initial_solution = true,
+                                     save_final_solution = true,
+                                     solution_variables = cons2prim)
+
+# The StepsizeCallback handles the re-calculation of the maximum Δt after each time step
+stepsize_callback = StepsizeCallback(cfl = 1.6)
+
+# Create a CallbackSet to collect all callbacks such that they can be passed to the ODE solver
+callbacks = CallbackSet(summary_callback,
+                        analysis_callback, alive_callback,
+                        save_restart, save_solution,
+                        stepsize_callback)
+
+###############################################################################
+# run the simulation
+
+# OrdinaryDiffEq's `solve` method evolves the solution in time and executes the passed callbacks
+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);
+
+# Print the timer summary
+summary_callback()
+
+# Finalize `T8codeMesh` to make sure MPI related objects in t8code are
+# released before `MPI` finalizes.
+!isinteractive() && finalize(mesh)

examples/t8code_2d_dgsem/elixir_advection_nonconforming_flag.jl

@@ -69,11 +69,16 @@ summary_callback = SummaryCallback()
 # The AnalysisCallback allows to analyse the solution in regular intervals and prints the results
 analysis_callback = AnalysisCallback(semi, interval = 100)
 
+# The SaveSolutionCallback allows to save the solution to a file in regular intervals
+save_solution = SaveSolutionCallback(interval = 100,
+                                     solution_variables = cons2prim)
+
 # The StepsizeCallback handles the re-calculation of the maximum Δt after each time step
 stepsize_callback = StepsizeCallback(cfl = 1.6)
 
 # Create a CallbackSet to collect all callbacks such that they can be passed to the ODE solver
-callbacks = CallbackSet(summary_callback, analysis_callback, stepsize_callback)
+callbacks = CallbackSet(summary_callback, analysis_callback, save_solution,
+                        stepsize_callback)
 
 ###############################################################################
 # run the simulation
@@ -85,3 +90,7 @@ sol = solve(ode, CarpenterKennedy2N54(williamson_condition = false),
 
 # Print the timer summary
 summary_callback()
+
+# Finalize `T8codeMesh` to make sure MPI related objects in t8code are
+# released before `MPI` finalizes.
+!isinteractive() && finalize(mesh)

examples/t8code_2d_dgsem/elixir_advection_restart.jl

@@ -0,0 +1,47 @@
+
+using OrdinaryDiffEq
+using Trixi
+
+###############################################################################
+# create a restart file
+
+elixir_file = "elixir_advection_extended.jl"
+restart_file = "restart_000000021.h5"
+
+trixi_include(@__MODULE__, joinpath(@__DIR__, elixir_file))
+
+###############################################################################
+# adapt the parameters that have changed compared to "elixir_advection_extended.jl"
+
+# Note: If you get a restart file from somewhere else, you need to provide
+# appropriate setups in the elixir loading a restart file
+
+restart_filename = joinpath("out", restart_file)
+mesh = load_mesh(restart_filename)
+
+semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver,
+                                    boundary_conditions = boundary_conditions)
+
+tspan = (load_time(restart_filename), 2.0)
+dt = load_dt(restart_filename)
+ode = semidiscretize(semi, tspan, restart_filename);
+
+# Do not overwrite the initial snapshot written by elixir_advection_extended.jl.
+save_solution.condition.save_initial_solution = false
+
+integrator = init(ode, CarpenterKennedy2N54(williamson_condition = false),
+                  dt = dt, # solve needs some value here but it will be overwritten by the stepsize_callback
+                  save_everystep = false, callback = callbacks, maxiters = 100_000);
+
+# Get the last time index and work with that.
+load_timestep!(integrator, restart_filename)
+
+###############################################################################
+# run the simulation
+
+sol = solve!(integrator)
+summary_callback() # print the timer summary
+
+# Finalize `T8codeMesh` to make sure MPI related objects in t8code are
+# released before `MPI` finalizes.
+!isinteractive() && finalize(mesh)

examples/t8code_2d_dgsem/elixir_advection_restart_amr.jl

@@ -0,0 +1,66 @@
+
+using OrdinaryDiffEq
+using Trixi
+
+###############################################################################
+# create a restart file
+
+elixir_file = "elixir_advection_extended.jl"
+restart_file = "restart_000000021.h5"
+
+trixi_include(@__MODULE__, joinpath(@__DIR__, elixir_file))
+
+###############################################################################
+# adapt the parameters that have changed compared to "elixir_advection_extended.jl"
+
+# Note: If you get a restart file from somewhere else, you need to provide
+# appropriate setups in the elixir loading a restart file
+
+restart_filename = joinpath("out", restart_file)
+mesh = load_mesh(restart_filename)
+
+semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver,
+                                    boundary_conditions = boundary_conditions)
+
+tspan = (load_time(restart_filename), 2.0)
+dt = load_dt(restart_filename)
+ode = semidiscretize(semi, tspan, restart_filename);
+
+# Do not overwrite the initial snapshot written by elixir_advection_extended.jl.
+save_solution.condition.save_initial_solution = false
+
+# Add AMR callback
+amr_controller = ControllerThreeLevel(semi, IndicatorMax(semi, variable = first),
+                                      base_level = 0,
+                                      med_level = 0, med_threshold = 0.8,
+                                      max_level = 1, max_threshold = 1.2)
+amr_callback = AMRCallback(semi, amr_controller,
+                           interval = 5,
+                           adapt_initial_condition = true,
+                           adapt_initial_condition_only_refine = true,
+                           dynamic_load_balancing = false)
+# We disable `dynamic_load_balancing` for now, since t8code does not support
+# partitioning for coarsening yet. That is, a complete family of elements always
+# stays on rank and is not split up due to partitioning. Without this feature
+# dynamic AMR simulations are not pefectly deterministic regarding to
+# convergent tests. Once this feature is available in t8code load balancing is
+# enabled again.
+
+callbacks_ext = CallbackSet(amr_callback, callbacks.discrete_callbacks...)
+
+integrator = init(ode, CarpenterKennedy2N54(williamson_condition = false),
+                  dt = dt, # solve needs some value here but it will be overwritten by the stepsize_callback
+                  save_everystep = false, callback = callbacks_ext, maxiters = 100_000);
+
+# Get the last time index and work with that.
+load_timestep!(integrator, restart_filename)
+
+###############################################################################
+# run the simulation
+
+sol = solve!(integrator)
+summary_callback() # print the timer summary
+
+# Finalize `T8codeMesh` to make sure MPI related objects in t8code are
+# released before `MPI` finalizes.
+!isinteractive() && finalize(mesh)

examples/t8code_2d_dgsem/elixir_advection_unstructured_flag.jl

@@ -53,13 +53,18 @@ summary_callback = SummaryCallback()
 # prints the results.
 analysis_callback = AnalysisCallback(semi, interval = 100)
 
+# The SaveSolutionCallback allows to save the solution to a file in regular intervals
+save_solution = SaveSolutionCallback(interval = 100,
+                                     solution_variables = cons2prim)
+
 # The StepsizeCallback handles the re-calculation of the maximum Δt after each
 # time step.
 stepsize_callback = StepsizeCallback(cfl = 1.4)
 
 # Create a CallbackSet to collect all callbacks such that they can be passed to
 # the ODE solver.
-callbacks = CallbackSet(summary_callback, analysis_callback, stepsize_callback)
+callbacks = CallbackSet(summary_callback, analysis_callback, save_solution,
+                        stepsize_callback)
 
 ###############################################################################
 # Run the simulation.
@@ -71,3 +76,7 @@ sol = solve(ode, CarpenterKennedy2N54(williamson_condition = false),
 
 # Print the timer summary.
 summary_callback()
+
+# Finalize `T8codeMesh` to make sure MPI related objects in t8code are
+# released before `MPI` finalizes.
+!isinteractive() && finalize(mesh)

examples/t8code_2d_dgsem/elixir_euler_free_stream.jl

@@ -72,10 +72,16 @@ analysis_callback = AnalysisCallback(semi, interval = analysis_interval)
 
 alive_callback = AliveCallback(analysis_interval = analysis_interval)
 
+save_solution = SaveSolutionCallback(interval = 100,
+                                     save_initial_solution = true,
+                                     save_final_solution = true,
+                                     solution_variables = cons2prim)
+
 stepsize_callback = StepsizeCallback(cfl = 2.0)
 
 callbacks = CallbackSet(summary_callback,
                         analysis_callback, alive_callback,
+                        save_solution,
                         stepsize_callback)
 
 ###############################################################################
@@ -85,3 +91,7 @@ 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
+
+# Finalize `T8codeMesh` to make sure MPI related objects in t8code are
+# released before `MPI` finalizes.
+!isinteractive() && finalize(mesh)