WIP: Taal

.github/workflows/CompatHelper.yml

@@ -13,4 +13,4 @@ jobs:
         env:
           GITHUB_TOKEN: ${{ secrets.GITHUB_TOKEN }}
           COMPATHELPER_PRIV: ${{ secrets.DOCUMENTER_KEY }}  # optional
-        run: julia -e 'using CompatHelper; CompatHelper.main(; subdirs = ["", "test", "docs"])'
+        run: julia -e 'using CompatHelper; CompatHelper.main(; subdirs=["", "docs", "examples", "test"])'

.gitignore

@@ -8,7 +8,7 @@
 *.avi
 *.ogv
 **/Manifest.toml
-out/
+out*/
 docs/build
 public/
 coverage/

Project.toml

@@ -1,9 +1,11 @@
 name = "Trixi"
 uuid = "a7f1ee26-1774-49b1-8366-f1abc58fbfcb"
 authors = ["Michael Schlottke-Lakemper <[email protected]>", "Gregor Gassner <[email protected]>", "Hendrik Ranocha <[email protected]>", "Andrew R. Winters <[email protected]>"]
-version = "0.2.6-pre"
+version = "0.3.0-pre"
 
 [deps]
+DiffEqBase = "2b5f629d-d688-5b77-993f-72d75c75574e"
+DiffEqCallbacks = "459566f4-90b8-5000-8ac3-15dfb0a30def"
 EllipsisNotation = "da5c29d0-fa7d-589e-88eb-ea29b0a81949"
 HDF5 = "f67ccb44-e63f-5c2f-98bd-6dc0ccc4ba2f"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
@@ -19,6 +21,8 @@ Tullio = "bc48ee85-29a4-5162-ae0b-a64e1601d4bc"
 UnPack = "3a884ed6-31ef-47d7-9d2a-63182c4928ed"
 
 [compat]
+DiffEqBase = "6.47"
+DiffEqCallbacks = "2.14"
 EllipsisNotation = "0.4"
 HDF5 = "0.13"
 LinearMaps = "2.7"

README.md

@@ -29,7 +29,10 @@ installation and postprocessing procedures. Its features include:
   * Kinetic energy-preserving and entropy-stable split forms
   * Entropy-stable shock capturing
   * Positivity-preserving limiting
-* Explicit low-storage Runge-Kutta time integration
+* Compatible with the [SciML ecosystem for ordinary differential equations](https://diffeq.sciml.ai/latest/)
+  * [Explicit low-storage Runge-Kutta time integration](https://diffeq.sciml.ai/latest/solvers/ode_solve/#Explicit-Strong-Stability-Preserving-Runge-Kutta-Methods-for-Hyperbolic-PDEs-(Conservation-Laws))
+  * [Strong stability preserving methods](https://diffeq.sciml.ai/latest/solvers/ode_solve/#Low-Storage-Methods)
+  * CFL-based and error-based time step control
 * Square/cubic domains with periodic and weakly-enforced boundary conditions
 * Multiple governing equations:
   * Compressible Euler equations

examples/2d/jeans_instability.jl

@@ -0,0 +1,99 @@
+
+using OrdinaryDiffEq
+using Trixi
+
+###############################################################################
+# semidiscretization of the compressible Euler equations
+gamma = 5/3
+equations_euler = CompressibleEulerEquations2D(gamma)
+
+# TODO: Taal, define initial_conditions_jeans_instability here for Euler
+initial_conditions = Trixi.initial_conditions_jeans_instability
+
+polydeg = 3
+solver_euler = DGSEM(polydeg, flux_hll)
+
+coordinates_min = (0, 0)
+coordinates_max = (1, 1)
+mesh = TreeMesh(coordinates_min, coordinates_max,
+                initial_refinement_level=4,
+                n_cells_max=10_000)
+
+semi_euler = SemidiscretizationHyperbolic(mesh, equations_euler, initial_conditions, solver_euler)
+
+
+###############################################################################
+# semidiscretization of the hyperbolic diffusion equations
+resid_tol = 1.0e-4
+equations_gravity = HyperbolicDiffusionEquations2D(resid_tol)
+
+# TODO: Taal, define initial_conditions_jeans_instability here for gravity
+
+solver_gravity = DGSEM(polydeg, flux_lax_friedrichs)
+
+semi_gravity = SemidiscretizationHyperbolic(mesh, equations_gravity, initial_conditions, solver_gravity,
+                                            source_terms=source_terms_harmonic)
+
+
+###############################################################################
+# combining both semidiscretizations for Euler + self-gravity
+parameters = ParametersEulerGravity(background_density=1.5e7, # aka rho0
+                                    gravitational_constant=6.674e-8, # aka G
+                                    cfl=2.4,
+                                    n_iterations_max=1000,
+                                    timestep_gravity=timestep_gravity_erk52_3Sstar!)
+
+semi = SemidiscretizationEulerGravity(semi_euler, semi_gravity, parameters)
+
+
+###############################################################################
+# ODE solvers, callbacks etc.
+tspan = (0.0, 5.0)
+ode = semidiscretize(semi, tspan);
+
+summary_callback = SummaryCallback()
+
+stepsize_callback = StepsizeCallback(cfl=1.0)
+
+save_solution = SaveSolutionCallback(interval=10,
+                                     save_initial_solution=true,
+                                     save_final_solution=true,
+                                     solution_variables=:primitive)
+# TODO: Taal, IO
+# restart_interval = 10
+
+analysis_interval = 100
+alive_callback = AliveCallback(analysis_interval=analysis_interval)
+
+Trixi.pretty_form_repl(::Val{:energy_potential}) = "∑e_potential"
+Trixi.pretty_form_file(::Val{:energy_potential}) = "e_potential"
+
+function Trixi.analyze(::Val{:energy_potential}, du, u_euler, t, semi::SemidiscretizationEulerGravity)
+
+  u_gravity = Trixi.wrap_array(semi.cache.u_ode, semi.semi_gravity)
+
+  mesh, equations_euler, dg, cache = Trixi.mesh_equations_solver_cache(semi.semi_euler)
+  _, equations_gravity, _, _ = Trixi.mesh_equations_solver_cache(semi.semi_gravity)
+
+  e_potential = Trixi.integrate(mesh, equations_euler, dg, cache, u_euler, equations_gravity, u_gravity) do u, i, j, element, equations_euler, dg, equations_gravity, u_gravity
+    u_euler_local   = Trixi.get_node_vars(u_euler,   equations_euler,   dg, i, j, element)
+    u_gravity_local = Trixi.get_node_vars(u_gravity, equations_gravity, dg, i, j, element)
+    # OBS! subtraction is specific to Jeans instability test where rho0 = 1.5e7
+    return (u_euler_local[1] - 1.5e7) * u_gravity_local[1]
+  end
+  return e_potential
+end
+
+analysis_callback = AnalysisCallback(semi_euler, interval=analysis_interval,
+                                     save_analysis=true,
+                                     extra_analysis_integrals=(entropy, energy_total, energy_kinetic, energy_internal, Val(:energy_potential)))
+
+callbacks = CallbackSet(summary_callback, stepsize_callback, save_solution, analysis_callback, alive_callback)
+
+
+###############################################################################
+# run the simulation
+sol = solve(ode, CarpenterKennedy2N54(williamson_condition=false), dt=stepsize_callback(ode),
+            save_everystep=false, callback=callbacks);
+summary_callback() # print the timer summary
+println("Number of gravity subcycles: ", semi.gravity_counter.ncalls_since_readout)

examples/2d/parameters.jl

@@ -0,0 +1,61 @@
+# TODO: Taal refactor, rename to
+# - linear_advection.jl
+# - advection_basic.jl
+# or something similar? parameters.jl isn't really helpful...
+
+using OrdinaryDiffEq
+using Trixi
+
+###############################################################################
+# semidiscretization of the linear advection equation
+
+advectionvelocity = (1.0, 1.0)
+# advectionvelocity = (0.2, -0.3)
+equations = LinearScalarAdvectionEquation2D(advectionvelocity)
+
+initial_conditions = initial_conditions_convergence_test
+
+surface_flux = flux_lax_friedrichs
+solver = DGSEM(3, surface_flux)
+
+coordinates_min = (-1, -1)
+coordinates_max = ( 1,  1)
+mesh = TreeMesh(coordinates_min, coordinates_max,
+                initial_refinement_level=4,
+                n_cells_max=30_000)
+
+
+semi = SemidiscretizationHyperbolic(mesh, equations, initial_conditions, solver)
+
+
+###############################################################################
+# ODE solvers, callbacks etc.
+
+tspan = (0.0, 1.0)
+ode = semidiscretize(semi, tspan);
+
+summary_callback = SummaryCallback()
+
+stepsize_callback = StepsizeCallback(cfl=1.6)
+
+save_solution = SaveSolutionCallback(interval=100,
+                                     save_initial_solution=true,
+                                     save_final_solution=true,
+                                     solution_variables=:conservative)
+# TODO: Taal, IO
+# restart_interval = 10
+
+analysis_interval = 100
+alive_callback = AliveCallback(analysis_interval=analysis_interval)
+analysis_callback = AnalysisCallback(semi, interval=analysis_interval,
+                                     extra_analysis_integrals=(entropy, energy_total))
+
+callbacks = CallbackSet(summary_callback, stepsize_callback, save_solution, analysis_callback, alive_callback)
+
+
+###############################################################################
+# run the simulation
+
+sol = solve(ode, CarpenterKennedy2N54(williamson_condition=false), dt=stepsize_callback(ode),
+            save_everystep=false, callback=callbacks);
+summary_callback() # print the timer summary

examples/2d/parameters_amr.jl

@@ -0,0 +1,71 @@
+# TODO: Taal refactor, rename to
+# - linear_advection_amr.jl
+# - advection_amr.jl
+# or something similar?
+
+using OrdinaryDiffEq
+using Trixi
+
+###############################################################################
+# semidiscretization of the linear advection equation
+
+advectionvelocity = (1.0, 1.0)
+# advectionvelocity = (0.2, -0.3)
+equations = LinearScalarAdvectionEquation2D(advectionvelocity)
+
+initial_conditions = initial_conditions_gauss
+
+surface_flux = flux_lax_friedrichs
+solver = DGSEM(3, surface_flux)
+
+coordinates_min = (-5, -5)
+coordinates_max = ( 5,  5)
+mesh = TreeMesh(coordinates_min, coordinates_max,
+                initial_refinement_level=4,
+                n_cells_max=30_000)
+
+
+semi = SemidiscretizationHyperbolic(mesh, equations, initial_conditions, solver)
+
+
+###############################################################################
+# ODE solvers, callbacks etc.
+
+tspan = (0.0, 10.0)
+ode = semidiscretize(semi, tspan)
+
+summary_callback = SummaryCallback()
+
+amr_indicator = IndicatorThreeLevel(semi, IndicatorMax(semi),
+                                    base_level=4,
+                                    med_level=5, med_threshold=0.1,
+                                    max_level=6, max_threshold=0.6)
+amr_callback = AMRCallback(semi, amr_indicator,
+                           interval=5,
+                           adapt_initial_conditions=true,
+                           adapt_initial_conditions_only_refine=true)
+
+stepsize_callback = StepsizeCallback(cfl=1.6)
+
+save_solution = SaveSolutionCallback(interval=100,
+                                     save_initial_solution=true,
+                                     save_final_solution=true,
+                                     solution_variables=:primitive)
+# TODO: Taal, IO
+# restart_interval = 10
+
+analysis_interval = 100
+alive_callback = AliveCallback(analysis_interval=analysis_interval)
+analysis_callback = AnalysisCallback(semi, interval=analysis_interval,
+                                     extra_analysis_integrals=(entropy,))
+
+# TODO: Taal decide, first AMR or save solution etc.
+callbacks = CallbackSet(summary_callback, amr_callback, stepsize_callback, save_solution, analysis_callback, alive_callback);
+
+
+###############################################################################
+# run the simulation
+
+sol = solve(ode, CarpenterKennedy2N54(williamson_condition=false), dt=stepsize_callback(ode),
+            save_everystep=false, callback=callbacks);
+summary_callback() # print the timer summary

examples/2d/parameters_blast_wave_shockcapturing.jl

@@ -0,0 +1,65 @@
+# TODO: Taal refactor, rename to
+# - euler_blast_wave_shockcapturing.jl
+# or something similar?
+
+using OrdinaryDiffEq
+using Trixi
+
+###############################################################################
+# semidiscretization of the compressible Euler equations
+
+equations = CompressibleEulerEquations2D(1.4)
+
+initial_conditions = initial_conditions_blast_wave
+
+surface_flux = flux_lax_friedrichs
+volume_flux  = flux_chandrashekar
+basis = LobattoLegendreBasis(3)
+indicator_sc = IndicatorHennemannGassner(equations, basis,
+                                         alpha_max=0.5,
+                                         alpha_min=0.001,
+                                         alpha_smooth=true,
+                                         variable=density_pressure)
+volume_integral = VolumeIntegralShockCapturingHG(indicator_sc;
+                                                 volume_flux_dg=volume_flux,
+                                                 volume_flux_fv=surface_flux)
+solver = DGSEM(basis, surface_flux, volume_integral)
+
+coordinates_min = (-2, -2)
+coordinates_max = ( 2,  2)
+mesh = TreeMesh(coordinates_min, coordinates_max,
+                initial_refinement_level=6,
+                n_cells_max=10_000)
+
+
+semi = SemidiscretizationHyperbolic(mesh, equations, initial_conditions, solver)
+
+
+###############################################################################
+# ODE solvers, callbacks etc.
+
+tspan = (0.0, 12.5)
+ode = semidiscretize(semi, tspan)
+
+summary_callback = SummaryCallback()
+
+stepsize_callback = StepsizeCallback(cfl=1.0)
+
+save_solution = SaveSolutionCallback(interval=100,
+                                     save_initial_solution=true,
+                                     save_final_solution=true,
+                                     solution_variables=:primitive)
+
+analysis_interval = 100
+alive_callback = AliveCallback(analysis_interval=analysis_interval)
+analysis_callback = AnalysisCallback(semi, interval=analysis_interval)
+
+callbacks = CallbackSet(summary_callback, stepsize_callback, save_solution, analysis_callback, alive_callback)
+
+
+###############################################################################
+# run the simulation
+
+sol = solve(ode, CarpenterKennedy2N54(williamson_condition=false), dt=stepsize_callback(ode),
+            save_everystep=false, callback=callbacks);
+summary_callback() # print the timer summary

examples/2d/parameters_blast_wave_shockcapturing.toml

@@ -35,6 +35,6 @@ initial_refinement_level = 6
 ####################################################################################################
 # I/O
 # save_initial_solution = false
-solution_interval = 10
+solution_interval = 100
 solution_variables = "primitive"
 restart_interval = 0