Initial support for surface coupling of two systems

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

@@ -115,6 +115,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")
@@ -177,7 +178,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
@@ -218,12 +220,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

@@ -83,9 +83,20 @@ 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", cb::DiscreteCallback{<:Any, <:Ana
 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,
@@ -127,8 +129,17 @@ function AnalysisCallback(mesh, equations::AbstractEquations, solver, cache;
 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)
 
@@ -197,14 +208,21 @@ function initialize!(cb::DiscreteCallback{Condition,Affect!}, u_ode, t, integrat
   # 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
@@ -289,15 +307,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)
-    du = wrap_array(du_ode,       mesh, equations, solver, cache)
-    l2_error, linf_error = analysis_callback(io, du, u, integrator.u, t, semi)
+    @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, u_ode, t, semi)
 
     mpi_println("─"^100)
     mpi_println()

src/callbacks_step/save_solution.jl

@@ -135,13 +135,7 @@ function initialize_save_cb!(solution_callback::SaveSolutionCallback, u, t, inte
   mpi_isroot() && mkpath(solution_callback.output_directory)
 
   semi = integrator.p
-  mesh, _, _, _ = mesh_equations_solver_cache(semi)
-  @trixi_timeit timer() "I/O" begin
-    if mesh.unsaved_changes
-      mesh.current_filename = save_mesh_file(mesh, solution_callback.output_directory)
-      mesh.unsaved_changes = false
-    end
-  end
+  @trixi_timeit timer() "I/O" save_mesh(semi, solution_callback.output_directory)
 
   if solution_callback.save_initial_solution
     solution_callback(integrator)
@@ -151,6 +145,17 @@ function initialize_save_cb!(solution_callback::SaveSolutionCallback, u, t, inte
 end
 
 
+# Save mesh for a general semidiscretization (default)
+function save_mesh(semi::AbstractSemidiscretization, output_directory, timestep=0)
+  mesh, _, _, _ = mesh_equations_solver_cache(semi)
+
+  if mesh.unsaved_changes
+    mesh.current_filename = save_mesh_file(mesh, output_directory)
+    mesh.unsaved_changes = false
+  end
+end
+
+
 # this method is called to determine whether the callback should be activated
 function (solution_callback::SaveSolutionCallback)(u, t, integrator)
   @unpack interval_or_dt, save_final_solution = solution_callback
@@ -169,32 +174,13 @@ end
 # this method is called when the callback is activated
 function (solution_callback::SaveSolutionCallback)(integrator)
   u_ode = integrator.u
-  @unpack t, dt = integrator
-  iter = integrator.stats.naccept
   semi = integrator.p
-  mesh, _, _, _ = mesh_equations_solver_cache(semi)
+  iter = integrator.stats.naccept
 
   @trixi_timeit timer() "I/O" begin
-    @trixi_timeit timer() "save mesh" if mesh.unsaved_changes
-      mesh.current_filename = save_mesh_file(mesh, solution_callback.output_directory, iter)
-      mesh.unsaved_changes = false
-    end
-
-    element_variables = Dict{Symbol, Any}()
-    @trixi_timeit timer() "get element variables" begin
-      get_element_variables!(element_variables, u_ode, semi)
-      callbacks = integrator.opts.callback
-      if callbacks isa CallbackSet
-        for cb in callbacks.continuous_callbacks
-          get_element_variables!(element_variables, u_ode, semi, cb; t=integrator.t, iter=integrator.stats.naccept)
-        end
-        for cb in callbacks.discrete_callbacks
-          get_element_variables!(element_variables, u_ode, semi, cb; t=integrator.t, iter=integrator.stats.naccept)
-        end
-      end
-    end
-
-    @trixi_timeit timer() "save solution" save_solution_file(u_ode, t, dt, iter, semi, solution_callback, element_variables)
+    # Call high-level functions that dispatch on semidiscretization type
+    @trixi_timeit timer() "save mesh" save_mesh(semi, solution_callback.output_directory, iter)
+    save_solution_file(semi, u_ode, solution_callback, integrator)
   end
 
   # avoid re-evaluating possible FSAL stages
@@ -203,12 +189,38 @@ function (solution_callback::SaveSolutionCallback)(integrator)
 end
 
 
+@inline function save_solution_file(semi::AbstractSemidiscretization, u_ode, solution_callback,
+                                    integrator; system="")
+  @unpack t, dt = integrator
+  iter = integrator.stats.naccept
+
+  element_variables = Dict{Symbol, Any}()
+  @trixi_timeit timer() "get element variables" begin
+    get_element_variables!(element_variables, u_ode, semi)
+    callbacks = integrator.opts.callback
+    if callbacks isa CallbackSet
+      for cb in callbacks.continuous_callbacks
+        get_element_variables!(element_variables, u_ode, semi, cb; t=integrator.t, iter=iter)
+      end
+      for cb in callbacks.discrete_callbacks
+        get_element_variables!(element_variables, u_ode, semi, cb; t=integrator.t, iter=iter)
+      end
+    end
+  end
+
+  @trixi_timeit timer() "save solution" save_solution_file(u_ode, t, dt, iter, semi,
+                                                           solution_callback, element_variables,
+                                                           system=system)
+end
+
+
 @inline function save_solution_file(u_ode, t, dt, iter,
                                     semi::AbstractSemidiscretization, solution_callback,
-                                    element_variables=Dict{Symbol,Any}())
+                                    element_variables=Dict{Symbol,Any}(); system="")
   mesh, equations, solver, cache = mesh_equations_solver_cache(semi)
   u = wrap_array_native(u_ode, mesh, equations, solver, cache)
-  save_solution_file(u, t, dt, iter, mesh, equations, solver, cache, solution_callback, element_variables)
+  save_solution_file(u, t, dt, iter, mesh, equations, solver, cache, solution_callback,
+                     element_variables; system=system)
 end
 
 

src/callbacks_step/stepsize.jl

@@ -68,13 +68,11 @@ end
     t = integrator.t
     u_ode = integrator.u
     semi = integrator.p
-    mesh, equations, solver, cache = mesh_equations_solver_cache(semi)
     @unpack cfl_number = stepsize_callback
-    u = wrap_array(u_ode, mesh, equations, solver, cache)
 
-    dt = @trixi_timeit timer() "calculate dt" cfl_number * max_dt(u, t, mesh,
-                                                                  have_constant_speed(equations), equations,
-                                                                  solver, cache)
+    # Dispatch based on semidiscretization
+    dt = @trixi_timeit timer() "calculate dt" calculate_dt(u_ode, t, cfl_number, semi)
+
     set_proposed_dt!(integrator, dt)
     integrator.opts.dtmax = dt
     integrator.dtcache = dt
@@ -86,6 +84,17 @@ end
 end
 
 
+# General case for a single semidiscretization
+function calculate_dt(u_ode, t, cfl_number, semi::AbstractSemidiscretization)
+  mesh, equations, solver, cache = mesh_equations_solver_cache(semi)
+  u = wrap_array(u_ode, mesh, equations, solver, cache)
+
+  dt = cfl_number * max_dt(u, t, mesh,
+                           have_constant_speed(equations), equations,
+                           solver, cache)
+end
+
+
 # Time integration methods from the DiffEq ecosystem without adaptive time stepping on their own
 # such as `CarpenterKennedy2N54` require passing `dt=...` in `solve(ode, ...)`. Since we don't have
 # an integrator at this stage but only the ODE, this method will be used there. It's called in