diff --git a/NEWS.md b/NEWS.md index 66dd1e58dc6..88afec79987 100644 --- a/NEWS.md +++ b/NEWS.md @@ -11,6 +11,7 @@ for human readability. to 1D and 3D on `TreeMesh`. - Implementation of 1D Linearized Euler Equations. - New analysis callback for 2D `P4estMesh` to compute integrated quantities along a boundary surface, e.g., pressure lift and drag coefficients. +- Optional tuple parameter for `GlmSpeedCallback` called `semi_indices` to specify for which semidiscretization of a `SemidiscretizationCoupled` we need to update the GLM speed. ## Changes when updating to v0.7 from v0.6.x diff --git a/Project.toml b/Project.toml index dc5fe949f2d..f4f060a0932 100644 --- a/Project.toml +++ b/Project.toml @@ -1,7 +1,7 @@ name = "Trixi" uuid = "a7f1ee26-1774-49b1-8366-f1abc58fbfcb" authors = ["Michael Schlottke-Lakemper ", "Gregor Gassner ", "Hendrik Ranocha ", "Andrew R. Winters ", "Jesse Chan "] -version = "0.7.10-pre" +version = "0.7.11-pre" [deps] CodeTracking = "da1fd8a2-8d9e-5ec2-8556-3022fb5608a2" diff --git a/docs/src/multi-physics_coupling.md b/docs/src/multi-physics_coupling.md index eec92bc21de..2c82c1b19e5 100644 --- a/docs/src/multi-physics_coupling.md +++ b/docs/src/multi-physics_coupling.md @@ -36,6 +36,17 @@ By passing the equations we can make use of their parameters, if they are requir Examples can be seen in `examples/structured_2d_dgsem/elixir_advection_coupled.jl`. +## GlmSpeedCallback for coupled MHD simulations + +When simulating an MHD system and the [`GlmSpeedCallback`](@ref) is required, +we need to specify for which semidiscretization we need the GLM speed updated. +This can be done with an additional parameter called `semi_indices`, which +is a tuple containing the semidiscretization indices for all systems +that require the GLM speed updated. + +An example elixir can be found at [`examples/structured_2d_dgsem/elixir_mhd_coupled.jl`](https://github.com/trixi-framework/Trixi.jl/blob/main/examples/structured_2d_dgsem/elixir_mhd_coupled.jl). + + ## Warning about binary compatibility Currently the coordinate values on the nodes can differ by machine precision when simulating the mesh and when splitting the mesh in multiple domains. diff --git a/examples/structured_2d_dgsem/elixir_mhd_coupled.jl b/examples/structured_2d_dgsem/elixir_mhd_coupled.jl new file mode 100644 index 00000000000..d3aa4ecf582 --- /dev/null +++ b/examples/structured_2d_dgsem/elixir_mhd_coupled.jl @@ -0,0 +1,136 @@ +using OrdinaryDiffEq +using Trixi + +############################################################################### +# Two semidiscretizations of the ideal GLM-MHD systems using converter functions such that +# they are coupled across the domain boundaries to generate a periodic system. +# +# In this elixir, we have a square domain that is divided into a left and right half. +# On each half of the domain, an independent SemidiscretizationHyperbolic is created for +# each set of ideal GLM-MHD equations. The two systems are coupled in the x-direction +# and are periodic in the y-direction. +# For a high-level overview, see also the figure below: +# +# (-2, 2) ( 2, 2) +# ┌────────────────────┬────────────────────┐ +# │ ↑ periodic ↑ │ ↑ periodic ↑ │ +# │ │ │ +# │ ========= │ ========= │ +# │ system #1 │ system #2 │ +# │ ========= │ ========= │ +# │ │ │ +# │<-- coupled │<-- coupled │ +# │ coupled -->│ coupled -->│ +# │ │ │ +# │ ↓ periodic ↓ │ ↓ periodic ↓ │ +# └────────────────────┴────────────────────┘ +# (-2, -2) ( 2, -2) + +gamma = 5 / 3 +equations = IdealGlmMhdEquations2D(gamma) + +cells_per_dimension = (32, 64) + +# Extend the definition of the non-conservative Powell flux functions. +import Trixi.flux_nonconservative_powell +function flux_nonconservative_powell(u_ll, u_rr, + normal_direction_ll::AbstractVector, + equations::IdealGlmMhdEquations2D) + flux_nonconservative_powell(u_ll, u_rr, normal_direction_ll, normal_direction_ll, + equations) +end +volume_flux = (flux_hindenlang_gassner, flux_nonconservative_powell) +solver = DGSEM(polydeg = 3, + surface_flux = (flux_lax_friedrichs, flux_nonconservative_powell), + volume_integral = VolumeIntegralFluxDifferencing(volume_flux)) + +########### +# system #1 +########### + +initial_condition1 = initial_condition_convergence_test +coordinates_min1 = (-1 / sin(pi / 4), -1 / sin(pi / 4)) +coordinates_max1 = (0.0, 1 / sin(pi / 4)) +mesh1 = StructuredMesh(cells_per_dimension, + coordinates_min1, + coordinates_max1, + periodicity = (false, true)) + +coupling_function1 = (x, u, equations_other, equations_own) -> u +boundary_conditions1 = (x_neg = BoundaryConditionCoupled(2, (:end, :i_forward), Float64, + coupling_function1), + x_pos = BoundaryConditionCoupled(2, (:begin, :i_forward), Float64, + coupling_function1), + y_neg = boundary_condition_periodic, + y_pos = boundary_condition_periodic) + +semi1 = SemidiscretizationHyperbolic(mesh1, equations, initial_condition1, solver, + boundary_conditions = boundary_conditions1) + +########### +# system #2 +########### + +initial_condition2 = initial_condition_convergence_test +coordinates_min2 = (0.0, -1 / sin(pi / 4)) +coordinates_max2 = (1 / sin(pi / 4), 1 / sin(pi / 4)) +mesh2 = StructuredMesh(cells_per_dimension, + coordinates_min2, + coordinates_max2, + periodicity = (false, true)) + +coupling_function2 = (x, u, equations_other, equations_own) -> u +boundary_conditions2 = (x_neg = BoundaryConditionCoupled(1, (:end, :i_forward), Float64, + coupling_function2), + x_pos = BoundaryConditionCoupled(1, (:begin, :i_forward), Float64, + coupling_function2), + y_neg = boundary_condition_periodic, + y_pos = boundary_condition_periodic) + +semi2 = SemidiscretizationHyperbolic(mesh2, equations, initial_condition2, solver, + boundary_conditions = boundary_conditions2) + +# Create a semidiscretization that bundles all the semidiscretizations. +semi = SemidiscretizationCoupled(semi1, semi2) + +############################################################################### +# ODE solvers, callbacks etc. + +tspan = (0.0, 0.1) +ode = semidiscretize(semi, tspan) + +summary_callback = SummaryCallback() + +analysis_interval = 100 + +analysis_callback1 = AnalysisCallback(semi1, interval = 100) +analysis_callback2 = AnalysisCallback(semi2, interval = 100) +analysis_callback = AnalysisCallbackCoupled(semi, analysis_callback1, analysis_callback2) + +alive_callback = AliveCallback(analysis_interval = analysis_interval) + +save_solution = SaveSolutionCallback(interval = 50, + save_initial_solution = true, + save_final_solution = true, + solution_variables = cons2prim) + +cfl = 1.0 + +stepsize_callback = StepsizeCallback(cfl = cfl) + +glm_speed_callback = GlmSpeedCallback(glm_scale = 0.5, cfl = cfl, + semi_indices = [1, 2]) + +callbacks = CallbackSet(summary_callback, + analysis_callback, alive_callback, + save_solution, + stepsize_callback, + glm_speed_callback) + +############################################################################### +# run the simulation + +sol = solve(ode, CarpenterKennedy2N54(williamson_condition = false), + dt = 0.01, # 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 diff --git a/src/callbacks_step/glm_speed.jl b/src/callbacks_step/glm_speed.jl index 036f61a522b..8ee406af5f9 100644 --- a/src/callbacks_step/glm_speed.jl +++ b/src/callbacks_step/glm_speed.jl @@ -6,7 +6,7 @@ #! format: noindent """ - GlmSpeedCallback(; glm_scale=0.5, cfl) + GlmSpeedCallback(; glm_scale=0.5, cfl, semi_indices=()) Update the divergence cleaning wave speed `c_h` according to the time step computed in [`StepsizeCallback`](@ref) for the ideal GLM-MHD equations. @@ -14,18 +14,26 @@ The `cfl` number should be set to the same value as for the time step size calcu `glm_scale` ensures that the GLM wave speed is lower than the fastest physical waves in the MHD solution and should thus be set to a value within the interval [0,1]. Note that `glm_scale = 0` deactivates the divergence cleaning. + +In case of coupled semidiscretizations, specify for which `semi_index`, i.e. index of the +semidiscretization, the divergence cleaning should be applied. See also +[`SemidiscretizationCoupled`](@ref). +Note: `SemidiscretizationCoupled` and all related features are considered experimental and +may change at any time. """ struct GlmSpeedCallback{RealT <: Real} glm_scale::RealT cfl::RealT + semi_indices::Vector{Int} end function Base.show(io::IO, cb::DiscreteCallback{<:Any, <:GlmSpeedCallback}) @nospecialize cb # reduce precompilation time glm_speed_callback = cb.affect! - @unpack glm_scale, cfl = glm_speed_callback - print(io, "GlmSpeedCallback(glm_scale=", glm_scale, ", cfl=", cfl, ")") + @unpack glm_scale, cfl, semi_indices = glm_speed_callback + print(io, "GlmSpeedCallback(glm_scale=", glm_scale, ", cfl=", cfl, "semi_indices=", + semi_indices, ")") end function Base.show(io::IO, ::MIME"text/plain", @@ -40,15 +48,16 @@ function Base.show(io::IO, ::MIME"text/plain", setup = [ "GLM wave speed scaling" => glm_speed_callback.glm_scale, "Expected CFL number" => glm_speed_callback.cfl, + "Selected semidiscretizations" => glm_speed_callback.semi_indices, ] summary_box(io, "GlmSpeedCallback", setup) end end -function GlmSpeedCallback(; glm_scale = 0.5, cfl) +function GlmSpeedCallback(; glm_scale = 0.5, cfl, semi_indices = Int[]) @assert 0<=glm_scale<=1 "glm_scale must be between 0 and 1" - glm_speed_callback = GlmSpeedCallback(glm_scale, cfl) + glm_speed_callback = GlmSpeedCallback(glm_scale, cfl, semi_indices) DiscreteCallback(glm_speed_callback, glm_speed_callback, # the first one is the condition, the second the affect! save_positions = (false, false), @@ -65,19 +74,29 @@ function (glm_speed_callback::GlmSpeedCallback)(u, t, integrator) return true end -# This method is called as callback after the StepsizeCallback during the time integration. -@inline function (glm_speed_callback::GlmSpeedCallback)(integrator) - dt = get_proposed_dt(integrator) - semi = integrator.p - mesh, equations, solver, cache = mesh_equations_solver_cache(semi) +function update_cleaning_speed!(semi, glm_speed_callback, dt) @unpack glm_scale, cfl = glm_speed_callback + mesh, equations, solver, cache = mesh_equations_solver_cache(semi) + # compute time step for GLM linear advection equation with c_h=1 (redone due to the possible AMR) c_h_deltat = calc_dt_for_cleaning_speed(cfl, mesh, equations, solver, cache) # c_h is proportional to its own time step divided by the complete MHD time step equations.c_h = glm_scale * c_h_deltat / dt + return semi +end + +# This method is called as callback after the StepsizeCallback during the time integration. +@inline function (glm_speed_callback::GlmSpeedCallback)(integrator) + dt = get_proposed_dt(integrator) + semi = integrator.p + + # Call the appropriate update function (this indirection allows to specialize on, + # e.g., the semidiscretization type) + update_cleaning_speed!(semi, glm_speed_callback, dt) + # avoid re-evaluating possible FSAL stages u_modified!(integrator, false) diff --git a/src/callbacks_step/time_series_dg.jl b/src/callbacks_step/time_series_dg.jl index 3781a10662d..5ba072bf560 100644 --- a/src/callbacks_step/time_series_dg.jl +++ b/src/callbacks_step/time_series_dg.jl @@ -9,9 +9,9 @@ function save_time_series_file(time_series_callback, mesh::Union{TreeMesh, UnstructuredMesh2D}, equations, dg::DG) - @unpack (interval, solution_variables, variable_names, + @unpack (interval, variable_names, output_directory, filename, point_coordinates, - point_data, time, step, time_series_cache) = time_series_callback + point_data, time, step) = time_series_callback n_points = length(point_data) h5open(joinpath(output_directory, filename), "w") do file diff --git a/src/semidiscretization/semidiscretization_coupled.jl b/src/semidiscretization/semidiscretization_coupled.jl index b8116d3cbd3..cc629d1a674 100644 --- a/src/semidiscretization/semidiscretization_coupled.jl +++ b/src/semidiscretization/semidiscretization_coupled.jl @@ -349,6 +349,36 @@ function calculate_dt(u_ode, t, cfl_number, semi::SemidiscretizationCoupled) return dt end +function update_cleaning_speed!(semi_coupled::SemidiscretizationCoupled, + glm_speed_callback, dt) + @unpack glm_scale, cfl, semi_indices = glm_speed_callback + + if length(semi_indices) == 0 + throw("Since you have more than one semidiscretization you need to specify the 'semi_indices' for which the GLM speed needs to be calculated.") + end + + # Check that all MHD semidiscretizations received a GLM cleaning speed update. + for (semi_index, semi) in enumerate(semi_coupled.semis) + if (typeof(semi.equations) <: AbstractIdealGlmMhdEquations && + !(semi_index in semi_indices)) + error("Equation of semidiscretization $semi_index needs to be included in 'semi_indices' of 'GlmSpeedCallback'.") + end + end + + for semi_index in semi_indices + semi = semi_coupled.semis[semi_index] + mesh, equations, solver, cache = mesh_equations_solver_cache(semi) + + # compute time step for GLM linear advection equation with c_h=1 (redone due to the possible AMR) + c_h_deltat = calc_dt_for_cleaning_speed(cfl, mesh, equations, solver, cache) + + # c_h is proportional to its own time step divided by the complete MHD time step + equations.c_h = glm_scale * c_h_deltat / dt + end + + return semi_coupled +end + ################################################################################ ### Equations ################################################################################ @@ -436,10 +466,28 @@ function (boundary_condition::BoundaryConditionCoupled)(u_inner, orientation, di Val(nvariables(equations)))) # Calculate boundary flux - if iseven(direction) # u_inner is "left" of boundary, u_boundary is "right" of boundary - flux = surface_flux_function(u_inner, u_boundary, orientation, equations) - else # u_boundary is "left" of boundary, u_inner is "right" of boundary - flux = surface_flux_function(u_boundary, u_inner, orientation, equations) + if surface_flux_function isa Tuple + # In case of conservative (index 1) and non-conservative (index 2) fluxes, + # add the non-conservative one with a factor of 1/2. + if iseven(direction) # u_inner is "left" of boundary, u_boundary is "right" of boundary + flux = (surface_flux_function[1](u_inner, u_boundary, orientation, + equations) + + 0.5 * + surface_flux_function[2](u_inner, u_boundary, orientation, + equations)) + else # u_boundary is "left" of boundary, u_inner is "right" of boundary + flux = (surface_flux_function[1](u_boundary, u_inner, orientation, + equations) + + 0.5 * + surface_flux_function[2](u_boundary, u_inner, orientation, + equations)) + end + else + if iseven(direction) # u_inner is "left" of boundary, u_boundary is "right" of boundary + flux = surface_flux_function(u_inner, u_boundary, orientation, equations) + else # u_boundary is "left" of boundary, u_inner is "right" of boundary + flux = surface_flux_function(u_boundary, u_inner, orientation, equations) + end end return flux diff --git a/src/time_integration/methods_SSP.jl b/src/time_integration/methods_SSP.jl index 9d1e06488b4..18f9c613249 100644 --- a/src/time_integration/methods_SSP.jl +++ b/src/time_integration/methods_SSP.jl @@ -88,7 +88,7 @@ mutable struct SimpleIntegratorSSP{RealT <: Real, uType, Params, Sol, F, Alg, t::RealT tdir::RealT dt::RealT # current time step - dtcache::RealT # ignored + dtcache::RealT # manually set time step iter::Int # current number of time steps (iteration) p::Params # will be the semidiscretization from Trixi sol::Sol # faked @@ -102,7 +102,7 @@ end """ add_tstop!(integrator::SimpleIntegratorSSP, t) -Add a time stop during the time integration process. +Add a time stop during the time integration process. This function is called after the periodic SaveSolutionCallback to specify the next stop to save the solution. """ function add_tstop!(integrator::SimpleIntegratorSSP, t) @@ -145,7 +145,7 @@ function solve(ode::ODEProblem, alg = SimpleSSPRK33()::SimpleAlgorithmSSP; t = first(ode.tspan) tdir = sign(ode.tspan[end] - ode.tspan[1]) iter = 0 - integrator = SimpleIntegratorSSP(u, du, r0, t, tdir, dt, zero(dt), iter, ode.p, + integrator = SimpleIntegratorSSP(u, du, r0, t, tdir, dt, dt, iter, ode.p, (prob = ode,), ode.f, alg, SimpleIntegratorSSPOptions(callback, ode.tspan; kwargs...), @@ -185,6 +185,8 @@ function solve!(integrator::SimpleIntegratorSSP) error("time step size `dt` is NaN") end + modify_dt_for_tstops!(integrator) + # if the next iteration would push the simulation beyond the end time, set dt accordingly if integrator.t + integrator.dt > t_end || isapprox(integrator.t + integrator.dt, t_end) @@ -192,8 +194,6 @@ function solve!(integrator::SimpleIntegratorSSP) terminate!(integrator) end - modify_dt_for_tstops!(integrator) - @. integrator.r0 = integrator.u for stage in eachindex(alg.c) t_stage = integrator.t + integrator.dt * alg.c[stage] @@ -232,7 +232,7 @@ function solve!(integrator::SimpleIntegratorSSP) end end - # Empty the tstops array. + # Empty the tstops array. # This cannot be done in terminate!(integrator::SimpleIntegratorSSP) because DiffEqCallbacks.PeriodicCallbackAffect would return at error. extract_all!(integrator.opts.tstops) @@ -253,12 +253,12 @@ u_modified!(integrator::SimpleIntegratorSSP, ::Bool) = false # used by adaptive timestepping algorithms in DiffEq function set_proposed_dt!(integrator::SimpleIntegratorSSP, dt) - integrator.dt = dt + (integrator.dt = dt; integrator.dtcache = dt) end # used by adaptive timestepping algorithms in DiffEq function get_proposed_dt(integrator::SimpleIntegratorSSP) - return integrator.dt + return ifelse(integrator.opts.adaptive, integrator.dt, integrator.dtcache) end # stop the time integration diff --git a/test/test_structured_2d.jl b/test/test_structured_2d.jl index d124b8f8bc3..f095c97b19e 100644 --- a/test/test_structured_2d.jl +++ b/test/test_structured_2d.jl @@ -941,6 +941,61 @@ end @test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 1000 end end + +@trixi_testset "elixir_mhd_coupled.jl" begin + @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_mhd_coupled.jl"), + l2=[ + 1.0743426980507015e-7, 0.030901698521864966, + 0.030901698662039206, 0.04370160129981656, + 8.259193827852516e-8, 0.03090169908364623, + 0.030901699039770684, 0.04370160128147447, + 8.735923402748945e-9, 1.0743426996067106e-7, + 0.03090169852186498, 0.030901698662039206, + 0.04370160129981657, 8.259193829690747e-8, + 0.03090169908364624, 0.030901699039770726, + 0.04370160128147445, 8.73592340076897e-9, + ], + linf=[ + 9.021023431587949e-7, 0.043701454182710486, + 0.043701458294527366, 0.061803146322536154, + 9.487023335807976e-7, 0.043701561010342616, + 0.04370147392153734, 0.06180318786081025, + 3.430673132525334e-8, 9.02102342825728e-7, + 0.043701454182710764, 0.043701458294525895, + 0.06180314632253597, 9.487023254761695e-7, + 0.04370156101034084, 0.04370147392153745, + 0.06180318786081015, 3.430672973680963e-8, + ], + coverage_override=(maxiters = 10^5,)) + + @testset "analysis_callback(sol) for AnalysisCallbackCoupled" begin + errors = analysis_callback(sol) + @test errors.l2≈[ + 1.0743426980507015e-7, 0.030901698521864966, 0.030901698662039206, + 0.04370160129981656, 8.259193827852516e-8, 0.03090169908364623, + 0.030901699039770684, 0.04370160128147447, 8.735923402748945e-9, + 1.0743426996067106e-7, 0.03090169852186498, 0.030901698662039206, + 0.04370160129981657, 8.259193829690747e-8, 0.03090169908364624, + 0.030901699039770726, 0.04370160128147445, 8.73592340076897e-9, + ] rtol=1.0e-4 + @test errors.linf≈[ + 9.021023431587949e-7, 0.043701454182710486, 0.043701458294527366, + 0.061803146322536154, 9.487023335807976e-7, 0.043701561010342616, + 0.04370147392153734, 0.06180318786081025, 3.430673132525334e-8, + 9.02102342825728e-7, 0.043701454182710764, 0.043701458294525895, + 0.06180314632253597, 9.487023254761695e-7, 0.04370156101034084, + 0.04370147392153745, 0.06180318786081015, 3.430672973680963e-8, + ] rtol=1.0e-4 + # Ensure that we do not have excessive memory allocations + # (e.g., from type instabilities) + let + t = sol.t[end] + u_ode = sol.u[end] + du_ode = similar(u_ode) + @test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 1000 + end + end +end end # Clean up afterwards: delete Trixi.jl output directory diff --git a/test/test_tree_2d_euler.jl b/test/test_tree_2d_euler.jl index 97a549ca4b8..d593dc540f7 100644 --- a/test/test_tree_2d_euler.jl +++ b/test/test_tree_2d_euler.jl @@ -235,6 +235,38 @@ end end end +@trixi_testset "elixir_euler_shockcapturing_subcell.jl (fixed time step)" begin + # Testing local SSP method without stepsize callback + # Additionally, tests combination with SaveSolutionCallback using time interval + @test_trixi_include(joinpath(EXAMPLES_DIR, + "elixir_euler_shockcapturing_subcell.jl"), + dt=2.0e-3, + tspan=(0.0, 0.25), + save_solution=SaveSolutionCallback(dt = 0.1 + 1.0e-8), + callbacks=CallbackSet(summary_callback, save_solution, + analysis_callback, alive_callback), + l2=[ + 0.05624855363458103, + 0.06931288786158463, + 0.06931283188960778, + 0.6200535829842072, + ], + linf=[ + 0.29029967648805566, + 0.6494728865862608, + 0.6494729363533714, + 3.0949621505674787, + ]) + # Ensure that we do not have excessive memory allocations + # (e.g., from type instabilities) + let + t = sol.t[end] + u_ode = sol.u[end] + du_ode = similar(u_ode) + @test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 15000 + end +end + @trixi_testset "elixir_euler_blast_wave.jl" begin @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_blast_wave.jl"), l2=[