Merge branch 'main' into jc/p4est_parabolic_BCs

.github/workflows/SpellCheck.yml

@@ -10,4 +10,4 @@ jobs:
       - name: Checkout Actions Repository
         uses: actions/checkout@v3
       - name: Check spelling
-        uses: crate-ci/typos@v1.14.12
+        uses: crate-ci/typos@v1.15.0

docs/src/visualization.md

@@ -339,7 +339,7 @@ create a [`PlotData1D`](@ref) with the keyword argument `curve` set to your list
 Let's give an example of this with the basic advection equation from above by creating
 a plot along the circle marked in green:
 
-![2d-plot-along-cirlce](https://user-images.githubusercontent.com/72009492/130951042-e1849447-8e55-4798-9361-c8badb9f3a49.png)
+![2d-plot-along-circle](https://user-images.githubusercontent.com/72009492/130951042-e1849447-8e55-4798-9361-c8badb9f3a49.png)
 
 We can write a function like this, that outputs a list of points on a circle:
 ```julia

examples/p4est_2d_dgsem/elixir_euler_supersonic_cylinder.jl

@@ -3,7 +3,7 @@
 # Boundary conditions are supersonic Mach 3 inflow at the left portion of the domain
 # and supersonic outflow at the right portion of the domain. The top and bottom of the
 # channel as well as the cylinder are treated as Euler slip wall boundaries.
-# This flow results in strong shock refletions / interactions as well as Kelvin-Helmholtz
+# This flow results in strong shock reflections / interactions as well as Kelvin-Helmholtz
 # instabilities at later times as two Mach stems form above and below the cylinder.
 #
 # For complete details on the problem setup see Section 5.7 of the paper:

examples/structured_2d_dgsem/elixir_advection_coupled.jl

@@ -0,0 +1,117 @@
+using OrdinaryDiffEq
+using Trixi
+
+
+###############################################################################
+# Coupled semidiscretization of two linear advection systems, which are connected periodically
+#
+# In this elixir, we have a square domain that is divided into a left half and a right half. On each
+# half of the domain, a completely independent SemidiscretizationHyperbolic is created for the
+# linear advection equations. The two systems are coupled in the x-direction and have periodic
+# boundaries in the y-direction. For a high-level overview, see also the figure below:
+#
+# (-1,  1)                                   ( 1,  1)
+#     ┌────────────────────┬────────────────────┐
+#     │    ↑ periodic ↑    │    ↑ periodic ↑    │
+#     │                    │                    │
+#     │                    │                    │
+#     │     =========      │     =========      │
+#     │     system #1      │     system #2      │
+#     │     =========      │     =========      │
+#     │                    │                    │
+#     │                    │                    │
+#     │                    │                    │
+#     │                    │                    │
+#     │         coupled -->│<-- coupled         │
+#     │                    │                    │
+#     │<-- coupled         │         coupled -->│
+#     │                    │                    │
+#     │                    │                    │
+#     │    ↓ periodic ↓    │    ↓ periodic ↓    │
+#     └────────────────────┴────────────────────┘
+# (-1, -1)                                   ( 1, -1)
+
+advection_velocity = (0.2, -0.7)
+equations = LinearScalarAdvectionEquation2D(advection_velocity)
+
+# 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)
+
+# First mesh is the left half of a [-1,1]^2 square
+coordinates_min1 = (-1.0, -1.0) # minimum coordinates (min(x), min(y))
+coordinates_max1 = ( 0.0,  1.0) # maximum coordinates (max(x), max(y))
+
+# Define identical resolution as a variable such that it is easier to change from `trixi_include`
+cells_per_dimension = (8, 16)
+
+cells_per_dimension1 = cells_per_dimension
+
+mesh1 = StructuredMesh(cells_per_dimension1, coordinates_min1, coordinates_max1)
+
+# A semidiscretization collects data structures and functions for the spatial discretization
+semi1 = SemidiscretizationHyperbolic(mesh1, equations, initial_condition_convergence_test, solver,
+                                     boundary_conditions=(
+                                       # Connect left boundary with right boundary of right mesh
+                                       x_neg=BoundaryConditionCoupled(2, (:end, :i_forward), Float64),
+                                       # Connect right boundary with left boundary of right mesh
+                                       x_pos=BoundaryConditionCoupled(2, (:begin, :i_forward),  Float64),
+                                       y_neg=boundary_condition_periodic,
+                                       y_pos=boundary_condition_periodic))
+
+
+# Second mesh is the right half of a [-1,1]^2 square
+coordinates_min2 = (0.0, -1.0) # minimum coordinates (min(x), min(y))
+coordinates_max2 = (1.0,  1.0) # maximum coordinates (max(x), max(y))
+
+cells_per_dimension2 = cells_per_dimension
+
+mesh2 = StructuredMesh(cells_per_dimension2, coordinates_min2, coordinates_max2)
+
+semi2 = SemidiscretizationHyperbolic(mesh2, equations, initial_condition_convergence_test, solver,
+                                     boundary_conditions=(
+                                       # Connect left boundary with right boundary of left mesh
+                                       x_neg=BoundaryConditionCoupled(1, (:end, :i_forward), Float64),
+                                       # Connect right boundary with left boundary of left mesh
+                                       x_pos=BoundaryConditionCoupled(1, (:begin, :i_forward),  Float64),
+                                       y_neg=boundary_condition_periodic,
+                                       y_pos=boundary_condition_periodic))
+
+# Create a semidiscretization that bundles semi1 and semi2
+semi = SemidiscretizationCoupled(semi1, semi2)
+
+###############################################################################
+# ODE solvers, callbacks etc.
+
+# Create ODE problem with time span from 0.0 to 2.0
+ode = semidiscretize(semi, (0.0, 2.0));
+
+# 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_callback1 = AnalysisCallback(semi1, interval=100)
+analysis_callback2 = AnalysisCallback(semi2, interval=100)
+analysis_callback = AnalysisCallbackCoupled(semi, analysis_callback1, analysis_callback2)
+
+# 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, 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()

src/Trixi.jl

@@ -117,6 +117,7 @@ include("semidiscretization/semidiscretization.jl")
 include("semidiscretization/semidiscretization_hyperbolic.jl")
 include("semidiscretization/semidiscretization_hyperbolic_parabolic.jl")
 include("semidiscretization/semidiscretization_euler_acoustics.jl")
+include("semidiscretization/semidiscretization_coupled.jl")
 include("callbacks_step/callbacks_step.jl")
 include("callbacks_stage/callbacks_stage.jl")
 include("semidiscretization/semidiscretization_euler_gravity.jl")
@@ -184,7 +185,8 @@ export boundary_condition_do_nothing,
        boundary_condition_noslip_wall,
        boundary_condition_slip_wall,
        boundary_condition_wall,
-       BoundaryConditionNavierStokesWall, NoSlip, Adiabatic, Isothermal
+       BoundaryConditionNavierStokesWall, NoSlip, Adiabatic, Isothermal,
+       BoundaryConditionCoupled
 
 export initial_condition_convergence_test, source_terms_convergence_test
 export source_terms_harmonic
@@ -229,12 +231,14 @@ export SemidiscretizationEulerAcoustics
 export SemidiscretizationEulerGravity, ParametersEulerGravity,
        timestep_gravity_erk52_3Sstar!, timestep_gravity_carpenter_kennedy_erk54_2N!
 
+export SemidiscretizationCoupled
+
 export SummaryCallback, SteadyStateCallback, AnalysisCallback, AliveCallback,
        SaveRestartCallback, SaveSolutionCallback, TimeSeriesCallback, VisualizationCallback,
        AveragingCallback,
        AMRCallback, StepsizeCallback,
        GlmSpeedCallback, LBMCollisionCallback, EulerAcousticsCouplingCallback,
-       TrivialCallback
+       TrivialCallback, AnalysisCallbackCoupled
 
 export load_mesh, load_time
 

src/auxiliary/special_elixirs.jl

@@ -85,9 +85,21 @@ function convergence_test(mod::Module, elixir::AbstractString, iterations; kwarg
         println("#"^100)
     end
 
-    # number of variables
-    _, equations, _, _ = mesh_equations_solver_cache(mod.semi)
+    # Use raw error values to compute EOC
+    analyze_convergence(errors, iterations, mod.semi)
+end
+
+# Analyze convergence for any semidiscretization
+# Note: this intermediate method is to allow dispatching on the semidiscretization
+function analyze_convergence(errors, iterations, semi::AbstractSemidiscretization)
+    _, equations, _, _ = mesh_equations_solver_cache(semi)
     variablenames = varnames(cons2cons, equations)
+    analyze_convergence(errors, iterations, variablenames)
+end
+
+# This method is called with the collected error values to actually compute and print the EOC
+function analyze_convergence(errors, iterations,
+                             variablenames::Union{Tuple, AbstractArray})
     nvariables = length(variablenames)
 
     # Reshape errors to get a matrix where the i-th row represents the i-th iteration

src/callbacks_step/analysis.jl

@@ -84,11 +84,13 @@ function Base.show(io::IO, ::MIME"text/plain",
     end
 end
 
+# This is the convenience constructor that gets called from the elixirs
 function AnalysisCallback(semi::AbstractSemidiscretization; kwargs...)
     mesh, equations, solver, cache = mesh_equations_solver_cache(semi)
     AnalysisCallback(mesh, equations, solver, cache; kwargs...)
 end
 
+# This is the actual constructor
 function AnalysisCallback(mesh, equations::AbstractEquations, solver, cache;
                           interval = 0,
                           save_analysis = false,
@@ -132,9 +134,18 @@ function AnalysisCallback(mesh, equations::AbstractEquations, solver, cache;
                      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}
     semi = integrator.p
+    du_ode = first(get_tmp_cache(integrator))
+    initialize!(cb, u_ode, du_ode, t, integrator, semi)
+end
+
+# This is the actual initialization method
+# Note: we have this indirection to allow initializing a callback from the AnalysisCallbackCoupled
+function initialize!(cb::DiscreteCallback{Condition, Affect!}, u_ode, du_ode, t,
+                     integrator, semi) where {Condition, Affect! <: AnalysisCallback}
     initial_state_integrals = integrate(u_ode, semi)
     _, equations, _, _ = mesh_equations_solver_cache(semi)
 
@@ -202,13 +213,21 @@ function initialize!(cb::DiscreteCallback{Condition, Affect!}, u_ode, t,
     # Note: For details see the actual callback function below
     analysis_callback.start_gc_time = Base.gc_time_ns()
 
-    analysis_callback(integrator)
+    analysis_callback(u_ode, du_ode, integrator, semi)
     return nothing
 end
 
-# TODO: Taal refactor, allow passing an IO object (which could be devnull to avoid cluttering the console)
+# This method gets called from OrdinaryDiffEq's `solve(...)`
 function (analysis_callback::AnalysisCallback)(integrator)
     semi = integrator.p
+    du_ode = first(get_tmp_cache(integrator))
+    u_ode = integrator.u
+    analysis_callback(u_ode, du_ode, integrator, semi)
+end
+
+# This method gets called internally as the main entry point to the AnalysiCallback
+# TODO: Taal refactor, allow passing an IO object (which could be devnull to avoid cluttering the console)
+function (analysis_callback::AnalysisCallback)(u_ode, du_ode, integrator, semi)
     mesh, equations, solver, cache = mesh_equations_solver_cache(semi)
     @unpack dt, t = integrator
     iter = integrator.stats.naccept
@@ -300,15 +319,14 @@ function (analysis_callback::AnalysisCallback)(integrator)
         end
 
         # Calculate current time derivative (needed for semidiscrete entropy time derivative, residual, etc.)
-        du_ode = first(get_tmp_cache(integrator))
         # `integrator.f` is usually just a call to `rhs!`
         # However, we want to allow users to modify the ODE RHS outside of Trixi.jl
         # and allow us to pass a combined ODE RHS to OrdinaryDiffEq, e.g., for
         # hyperbolic-parabolic systems.
-        @notimeit timer() integrator.f(du_ode, integrator.u, semi, t)
-        u = wrap_array(integrator.u, mesh, equations, solver, cache)
+        @notimeit timer() integrator.f(du_ode, u_ode, semi, t)
+        u = wrap_array(u_ode, mesh, equations, solver, cache)
         du = wrap_array(du_ode, mesh, equations, solver, cache)
-        l2_error, linf_error = analysis_callback(io, du, u, integrator.u, t, semi)
+        l2_error, linf_error = analysis_callback(io, du, u, u_ode, t, semi)
 
         mpi_println("─"^100)
         mpi_println()