From 3e7cb8e3d8e802d9af520114e0685576c5fd253e Mon Sep 17 00:00:00 2001 From: Daniel_Doehring Date: Mon, 12 Aug 2024 15:19:47 +0200 Subject: [PATCH 01/17] 1st version viscous cfl --- .../elixir_advection_diffusion_cfl.jl | 95 ++++++++++++ src/Trixi.jl | 2 +- src/callbacks_step/callbacks_step.jl | 1 + src/callbacks_step/stepsize_dg1d.jl | 24 ++- .../stepsize_hyperbolic_parabolic.jl | 140 ++++++++++++++++++ src/equations/laplace_diffusion_1d.jl | 6 + 6 files changed, 263 insertions(+), 5 deletions(-) create mode 100644 examples/tree_1d_dgsem/elixir_advection_diffusion_cfl.jl create mode 100644 src/callbacks_step/stepsize_hyperbolic_parabolic.jl diff --git a/examples/tree_1d_dgsem/elixir_advection_diffusion_cfl.jl b/examples/tree_1d_dgsem/elixir_advection_diffusion_cfl.jl new file mode 100644 index 00000000000..006e87c2fca --- /dev/null +++ b/examples/tree_1d_dgsem/elixir_advection_diffusion_cfl.jl @@ -0,0 +1,95 @@ + +using OrdinaryDiffEq +using Trixi + +############################################################################### +# semidiscretization of the linear advection diffusion equation + +advection_velocity = 0.1 +equations = LinearScalarAdvectionEquation1D(advection_velocity) +diffusivity() = 0.1 +equations_parabolic = LaplaceDiffusion1D(diffusivity(), equations) + +# 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) + +coordinates_min = -pi # minimum coordinate +coordinates_max = pi # maximum coordinate + +# Create a uniformly refined mesh with periodic boundaries +mesh = TreeMesh(coordinates_min, coordinates_max, + initial_refinement_level = 4, + n_cells_max = 30_000, # set maximum capacity of tree data structure + periodicity = true) + +function x_trans_periodic(x, domain_length = SVector(2 * pi), center = SVector(0.0)) + x_normalized = x .- center + x_shifted = x_normalized .% domain_length + x_offset = ((x_shifted .< -0.5 * domain_length) - (x_shifted .> 0.5 * domain_length)) .* + domain_length + return center + x_shifted + x_offset +end + +# Define initial condition +function initial_condition_diffusive_convergence_test(x, t, + equation::LinearScalarAdvectionEquation1D) + # Store translated coordinate for easy use of exact solution + x_trans = x_trans_periodic(x - equation.advection_velocity * t) + + nu = diffusivity() + c = 0.0 + A = 1.0 + L = 2 + f = 1 / L + omega = 1.0 + scalar = c + A * sin(omega * sum(x_trans)) * exp(-nu * omega^2 * t) + return SVector(scalar) +end +initial_condition = initial_condition_diffusive_convergence_test + +# define periodic boundary conditions everywhere +boundary_conditions = boundary_condition_periodic +boundary_conditions_parabolic = boundary_condition_periodic + +# A semidiscretization collects data structures and functions for the spatial discretization +semi = SemidiscretizationHyperbolicParabolic(mesh, (equations, equations_parabolic), + initial_condition, + solver; + boundary_conditions = (boundary_conditions, + boundary_conditions_parabolic)) + +############################################################################### +# 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_callback = AnalysisCallback(semi, interval = 100) + +# The AliveCallback prints short status information in regular intervals +alive_callback = AliveCallback(analysis_interval = 100) + +# Stepsize callback which selects the timestep to the most restrictive influence +stepsize_callback = StepsizeCallbackHyperbolicParabolic(cfl_convective = 1.6, + cfl_diffusive = 0.4) + +# 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, + 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() diff --git a/src/Trixi.jl b/src/Trixi.jl index 23a8cfe1d0a..623df916f8e 100644 --- a/src/Trixi.jl +++ b/src/Trixi.jl @@ -262,7 +262,7 @@ export SemidiscretizationCoupled export SummaryCallback, SteadyStateCallback, AnalysisCallback, AliveCallback, SaveRestartCallback, SaveSolutionCallback, TimeSeriesCallback, VisualizationCallback, AveragingCallback, - AMRCallback, StepsizeCallback, + AMRCallback, StepsizeCallback, StepsizeCallbackHyperbolicParabolic, GlmSpeedCallback, LBMCollisionCallback, EulerAcousticsCouplingCallback, TrivialCallback, AnalysisCallbackCoupled, AnalysisSurfaceIntegral, DragCoefficientPressure, LiftCoefficientPressure, diff --git a/src/callbacks_step/callbacks_step.jl b/src/callbacks_step/callbacks_step.jl index 09d197bf225..2b050acca5c 100644 --- a/src/callbacks_step/callbacks_step.jl +++ b/src/callbacks_step/callbacks_step.jl @@ -61,6 +61,7 @@ include("averaging.jl") include("amr.jl") include("stepsize.jl") +include("stepsize_hyperbolic_parabolic.jl") include("glm_speed.jl") include("lbm_collision.jl") include("euler_acoustics_coupling.jl") diff --git a/src/callbacks_step/stepsize_dg1d.jl b/src/callbacks_step/stepsize_dg1d.jl index edc25ec78f6..d0341b4fc6e 100644 --- a/src/callbacks_step/stepsize_dg1d.jl +++ b/src/callbacks_step/stepsize_dg1d.jl @@ -18,7 +18,7 @@ function max_dt(u, t, mesh::TreeMesh{1}, lambda1, = max_abs_speeds(u_node, equations) max_lambda1 = max(max_lambda1, lambda1) end - inv_jacobian = cache.elements.inverse_jacobian[element] + inv_jacobian = cache.elements.inverse_jacobian[element] # Δx max_scaled_speed = max(max_scaled_speed, inv_jacobian * max_lambda1) end @@ -33,13 +33,29 @@ function max_dt(u, t, mesh::TreeMesh{1}, for element in eachelement(dg, cache) max_lambda1, = max_abs_speeds(equations) - inv_jacobian = cache.elements.inverse_jacobian[element] + inv_jacobian = cache.elements.inverse_jacobian[element] # Δx max_scaled_speed = max(max_scaled_speed, inv_jacobian * max_lambda1) end return 2 / (nnodes(dg) * max_scaled_speed) end +function max_dt(u, t, mesh::TreeMesh{1}, + constant_diffusivity::True, + equations_parabolic::AbstractEquationsParabolic, dg::DG, cache) + # to avoid a division by zero if the speed vanishes everywhere, + # e.g. for steady-state linear advection + max_scaled_speed = nextfloat(zero(t)) + + for element in eachelement(dg, cache) + max_lambda1, = max_diffusivity(equations_parabolic) + inv_jacobian = cache.elements.inverse_jacobian[element] # Δx + max_scaled_speed = max(max_scaled_speed, inv_jacobian^2 * max_lambda1) + end + + return 2 / (nnodes(dg) * max_scaled_speed) +end + function max_dt(u, t, mesh::StructuredMesh{1}, constant_speed::False, equations, dg::DG, cache) # to avoid a division by zero if the speed vanishes everywhere, @@ -53,7 +69,7 @@ function max_dt(u, t, mesh::StructuredMesh{1}, u_node = get_node_vars(u, equations, dg, i, element) lambda1, = max_abs_speeds(u_node, equations) - inv_jacobian = cache.elements.inverse_jacobian[i, element] + inv_jacobian = cache.elements.inverse_jacobian[i, element] # Δx max_lambda1 = max(max_lambda1, inv_jacobian * lambda1) end @@ -74,7 +90,7 @@ function max_dt(u, t, mesh::StructuredMesh{1}, max_lambda1, = max_abs_speeds(equations) for i in eachnode(dg) - inv_jacobian = cache.elements.inverse_jacobian[i, element] + inv_jacobian = cache.elements.inverse_jacobian[i, element] # Δx max_scaled_speed = max(max_scaled_speed, inv_jacobian * max_lambda1) end end diff --git a/src/callbacks_step/stepsize_hyperbolic_parabolic.jl b/src/callbacks_step/stepsize_hyperbolic_parabolic.jl new file mode 100644 index 00000000000..998fc03c2b0 --- /dev/null +++ b/src/callbacks_step/stepsize_hyperbolic_parabolic.jl @@ -0,0 +1,140 @@ +# By default, Julia/LLVM does not use fused multiply-add operations (FMAs). +# Since these FMAs can increase the performance of many numerical algorithms, +# we need to opt-in explicitly. +# See https://ranocha.de/blog/Optimizing_EC_Trixi for further details. +@muladd begin +#! format: noindent + +""" + StepsizeCallbackHyperbolicParabolic(; cfl_convective=1.0, cfl_diffusive=1.0) + +Set the time step size according to a CFL condition with CFL number `cfl` +if the time integration method isn't adaptive itself. +""" +mutable struct StepsizeCallbackHyperbolicParabolic{RealT} + cfl_convective::RealT + cfl_diffusive::RealT +end + +function Base.show(io::IO, + cb::DiscreteCallback{<:Any, <:StepsizeCallbackHyperbolicParabolic}) + @nospecialize cb # reduce precompilation time + + stepsize_callback = cb.affect! + @unpack cfl_convective = stepsize_callback + @unpack cfl_diffusive = stepsize_callback + print(io, "StepsizeCallbackHyperbolicParabolic(cfl_convective=", cfl_convective, + "cfl_diffusive=", cfl_diffusive, ")") +end + +function Base.show(io::IO, ::MIME"text/plain", + cb::DiscreteCallback{<:Any, <:StepsizeCallbackHyperbolicParabolic}) + @nospecialize cb # reduce precompilation time + + if get(io, :compact, false) + show(io, cb) + else + stepsize_callback = cb.affect! + + setup = [ + "Convective CFL number" => stepsize_callback.cfl_convective, + "Diffusive CFL number" => stepsize_callback.cfl_diffusive, + ] + summary_box(io, "StepsizeCallbackHyperbolicParabolic", setup) + end +end + +function StepsizeCallbackHyperbolicParabolic(; + cfl_convective::Real = 1.0, + cfl_diffusive::Real = 1.0) + stepsize_callback = StepsizeCallbackHyperbolicParabolic(cfl_convective, + cfl_diffusive) + + DiscreteCallback(stepsize_callback, stepsize_callback, # the first one is the condition, the second the affect! + save_positions = (false, false), + initialize = initialize!) +end + +function initialize!(cb::DiscreteCallback{Condition, Affect!}, u, t, + integrator) where {Condition, + Affect! <: StepsizeCallbackHyperbolicParabolic} + cb.affect!(integrator) +end + +# this method is called to determine whether the callback should be activated +function (stepsize_callback::StepsizeCallbackHyperbolicParabolic)(u, t, integrator) + return true +end + +# This method is called as callback during the time integration. +@inline function (stepsize_callback::StepsizeCallbackHyperbolicParabolic)(integrator) + # TODO: Taal decide, shall we set the time step even if the integrator is adaptive? + if !integrator.opts.adaptive + t = integrator.t + u_ode = integrator.u + semi = integrator.p + @unpack cfl_convective, cfl_diffusive = stepsize_callback + + # Dispatch based on semidiscretization + dt = @trixi_timeit timer() "calculate dt" calculate_dt(u_ode, t, + cfl_convective, + cfl_diffusive, + semi) + + set_proposed_dt!(integrator, dt) + integrator.opts.dtmax = dt + integrator.dtcache = dt + end + + # avoid re-evaluating possible FSAL stages + u_modified!(integrator, false) + return nothing +end + +# Case for a hyperbolic-parabolic semidiscretization +function calculate_dt(u_ode, t, cfl_convective, cfl_diffusive, + semi::SemidiscretizationHyperbolicParabolic) + mesh, equations, solver, cache = mesh_equations_solver_cache(semi) + equations_parabolic = semi.equations_parabolic + + u = wrap_array(u_ode, mesh, equations, solver, cache) + + dt_convective = cfl_convective * max_dt(u, t, mesh, + have_constant_speed(equations), + equations, solver, cache) + + dt_diffusive = cfl_diffusive * max_dt(u, t, mesh, + have_constant_diffusivity(equations_parabolic), + equations_parabolic, solver, cache) + + return min(dt_convective, dt_diffusive) +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 +# many examples in `solve(ode, ..., dt=stepsize_callback(ode), ...)`. +function (cb::DiscreteCallback{Condition, Affect!})(ode::ODEProblem) where {Condition, + Affect! <: + StepsizeCallbackHyperbolicParabolic + } + stepsize_callback = cb.affect! + @unpack cfl_convective, cfl_diffusive = stepsize_callback + u_ode = ode.u0 + t = first(ode.tspan) + semi = ode.p + mesh, equations, solver, cache = mesh_equations_solver_cache(semi) + equations_parabolic = semi.equations_parabolic + u = wrap_array(u_ode, mesh, equations, solver, cache) + + dt_convective = dt_convective * + max_dt(u, t, mesh, have_constant_speed(equations), + equations, solver, cache) + + dt_diffusive = dt_diffusive * + max_dt(u, t, mesh, have_constant_diffusivity(equations_parabolic), + equations_parabolic, solver, cache) + + return min(dt_convective, dt_diffusive) +end +end # @muladd diff --git a/src/equations/laplace_diffusion_1d.jl b/src/equations/laplace_diffusion_1d.jl index 64a72ef3a13..5603ff9b5e6 100644 --- a/src/equations/laplace_diffusion_1d.jl +++ b/src/equations/laplace_diffusion_1d.jl @@ -24,6 +24,12 @@ function flux(u, gradients, orientation::Integer, equations_parabolic::LaplaceDi return equations_parabolic.diffusivity * dudx end +@inline have_constant_diffusivity(::LaplaceDiffusion1D) = True() + +@inline function max_diffusivity(equations_parabolic::LaplaceDiffusion1D) + return equations_parabolic.diffusivity +end + # Dirichlet and Neumann boundary conditions for use with parabolic solvers in weak form. # Note that these are general, so they apply to LaplaceDiffusion in any spatial dimension. @inline function (boundary_condition::BoundaryConditionDirichlet)(flux_inner, u_inner, From e8f6651583aad56d78e1c73f2c9baf36a607822f Mon Sep 17 00:00:00 2001 From: Daniel_Doehring Date: Mon, 12 Aug 2024 16:03:41 +0200 Subject: [PATCH 02/17] non-constant diffusivities --- .../elixir_advection_diffusion_cfl.jl | 95 ------------ ...r_navierstokes_convergence_periodic_cfl.jl | 141 ++++++++++++++++++ src/callbacks_step/stepsize_dg1d.jl | 44 +++++- .../stepsize_hyperbolic_parabolic.jl | 21 ++- .../compressible_navier_stokes_1d.jl | 9 ++ 5 files changed, 201 insertions(+), 109 deletions(-) delete mode 100644 examples/tree_1d_dgsem/elixir_advection_diffusion_cfl.jl create mode 100644 examples/tree_1d_dgsem/elixir_navierstokes_convergence_periodic_cfl.jl diff --git a/examples/tree_1d_dgsem/elixir_advection_diffusion_cfl.jl b/examples/tree_1d_dgsem/elixir_advection_diffusion_cfl.jl deleted file mode 100644 index 006e87c2fca..00000000000 --- a/examples/tree_1d_dgsem/elixir_advection_diffusion_cfl.jl +++ /dev/null @@ -1,95 +0,0 @@ - -using OrdinaryDiffEq -using Trixi - -############################################################################### -# semidiscretization of the linear advection diffusion equation - -advection_velocity = 0.1 -equations = LinearScalarAdvectionEquation1D(advection_velocity) -diffusivity() = 0.1 -equations_parabolic = LaplaceDiffusion1D(diffusivity(), equations) - -# 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) - -coordinates_min = -pi # minimum coordinate -coordinates_max = pi # maximum coordinate - -# Create a uniformly refined mesh with periodic boundaries -mesh = TreeMesh(coordinates_min, coordinates_max, - initial_refinement_level = 4, - n_cells_max = 30_000, # set maximum capacity of tree data structure - periodicity = true) - -function x_trans_periodic(x, domain_length = SVector(2 * pi), center = SVector(0.0)) - x_normalized = x .- center - x_shifted = x_normalized .% domain_length - x_offset = ((x_shifted .< -0.5 * domain_length) - (x_shifted .> 0.5 * domain_length)) .* - domain_length - return center + x_shifted + x_offset -end - -# Define initial condition -function initial_condition_diffusive_convergence_test(x, t, - equation::LinearScalarAdvectionEquation1D) - # Store translated coordinate for easy use of exact solution - x_trans = x_trans_periodic(x - equation.advection_velocity * t) - - nu = diffusivity() - c = 0.0 - A = 1.0 - L = 2 - f = 1 / L - omega = 1.0 - scalar = c + A * sin(omega * sum(x_trans)) * exp(-nu * omega^2 * t) - return SVector(scalar) -end -initial_condition = initial_condition_diffusive_convergence_test - -# define periodic boundary conditions everywhere -boundary_conditions = boundary_condition_periodic -boundary_conditions_parabolic = boundary_condition_periodic - -# A semidiscretization collects data structures and functions for the spatial discretization -semi = SemidiscretizationHyperbolicParabolic(mesh, (equations, equations_parabolic), - initial_condition, - solver; - boundary_conditions = (boundary_conditions, - boundary_conditions_parabolic)) - -############################################################################### -# 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_callback = AnalysisCallback(semi, interval = 100) - -# The AliveCallback prints short status information in regular intervals -alive_callback = AliveCallback(analysis_interval = 100) - -# Stepsize callback which selects the timestep to the most restrictive influence -stepsize_callback = StepsizeCallbackHyperbolicParabolic(cfl_convective = 1.6, - cfl_diffusive = 0.4) - -# 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, - 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() diff --git a/examples/tree_1d_dgsem/elixir_navierstokes_convergence_periodic_cfl.jl b/examples/tree_1d_dgsem/elixir_navierstokes_convergence_periodic_cfl.jl new file mode 100644 index 00000000000..b3b15ed2ef0 --- /dev/null +++ b/examples/tree_1d_dgsem/elixir_navierstokes_convergence_periodic_cfl.jl @@ -0,0 +1,141 @@ + +using OrdinaryDiffEq +using Trixi + +############################################################################### +# semidiscretization of the compressible Navier-Stokes equations + +prandtl_number() = 0.72 +mu() = 6.25e-4 # equivalent to Re = 1600 + +equations = CompressibleEulerEquations1D(1.4) +equations_parabolic = CompressibleNavierStokesDiffusion1D(equations, mu = mu(), + Prandtl = prandtl_number()) + +# This convergence test setup was originally derived by Andrew Winters (@andrewwinters5000) +# (Simplified version of the 2D) +function initial_condition_navier_stokes_convergence_test(x, t, equations) + # Amplitude and shift + A = 0.5 + c = 2.0 + + # convenience values for trig. functions + pi_x = pi * x[1] + pi_t = pi * t + + rho = c + A * sin(pi_x) * cos(pi_t) + v1 = sin(pi_x) * cos(pi_t) + p = rho^2 + + return prim2cons(SVector(rho, v1, p), equations) +end +initial_condition = initial_condition_navier_stokes_convergence_test + +@inline function source_terms_navier_stokes_convergence_test(u, x, t, equations) + # we currently need to hardcode these parameters until we fix the "combined equation" issue + # see also https://github.com/trixi-framework/Trixi.jl/pull/1160 + inv_gamma_minus_one = inv(equations.gamma - 1) + Pr = prandtl_number() + mu_ = mu() + + # Same settings as in `initial_condition` + # Amplitude and shift + A = 0.5 + c = 2.0 + + # convenience values for trig. functions + pi_x = pi * x[1] + pi_t = pi * t + + # compute the manufactured solution and all necessary derivatives + rho = c + A * sin(pi_x) * cos(pi_t) + rho_t = -pi * A * sin(pi_x) * sin(pi_t) + rho_x = pi * A * cos(pi_x) * cos(pi_t) + rho_xx = -pi * pi * A * sin(pi_x) * cos(pi_t) + + v1 = sin(pi_x) * cos(pi_t) + v1_t = -pi * sin(pi_x) * sin(pi_t) + v1_x = pi * cos(pi_x) * cos(pi_t) + v1_xx = -pi * pi * sin(pi_x) * cos(pi_t) + + p = rho * rho + p_t = 2.0 * rho * rho_t + p_x = 2.0 * rho * rho_x + p_xx = 2.0 * rho * rho_xx + 2.0 * rho_x * rho_x + + E = p * inv_gamma_minus_one + 0.5 * rho * v1^2 + E_t = p_t * inv_gamma_minus_one + 0.5 * rho_t * v1^2 + rho * v1 * v1_t + E_x = p_x * inv_gamma_minus_one + 0.5 * rho_x * v1^2 + rho * v1 * v1_x + + # Some convenience constants + T_const = equations.gamma * inv_gamma_minus_one / Pr + inv_rho_cubed = 1.0 / (rho^3) + + # compute the source terms + # density equation + du1 = rho_t + rho_x * v1 + rho * v1_x + + # x-momentum equation + du2 = (rho_t * v1 + rho * v1_t + + p_x + rho_x * v1^2 + 2.0 * rho * v1 * v1_x - + # stress tensor from x-direction + v1_xx * mu_) + + # total energy equation + du3 = (E_t + v1_x * (E + p) + v1 * (E_x + p_x) - + # stress tensor and temperature gradient terms from x-direction + v1_xx * v1 * mu_ - + v1_x * v1_x * mu_ - + T_const * inv_rho_cubed * + (p_xx * rho * rho - + 2.0 * p_x * rho * rho_x + + 2.0 * p * rho_x * rho_x - + p * rho * rho_xx) * mu_) + + return SVector(du1, du2, du3) +end + +volume_flux = flux_ranocha +solver = DGSEM(polydeg = 3, surface_flux = flux_hllc, + volume_integral = VolumeIntegralFluxDifferencing(volume_flux)) + +coordinates_min = -1.0 +coordinates_max = 1.0 +mesh = TreeMesh(coordinates_min, coordinates_max, + initial_refinement_level = 4, + n_cells_max = 100_000) + +semi = SemidiscretizationHyperbolicParabolic(mesh, (equations, equations_parabolic), + initial_condition, solver, + source_terms = source_terms_navier_stokes_convergence_test) + +############################################################################### +# ODE solvers, callbacks etc. + +tspan = (0.0, 10.0) +ode = semidiscretize(semi, tspan) + +summary_callback = SummaryCallback() + +analysis_interval = 1000 +analysis_callback = AnalysisCallback(semi, interval = analysis_interval) + +alive_callback = AliveCallback(analysis_interval = analysis_interval) + +# Stepsize callback which selects the timestep to the most restrictive influence +stepsize_callback = StepsizeCallbackHyperbolicParabolic(cfl_convective = 1.8, + cfl_diffusive = 0.3) + +callbacks = CallbackSet(summary_callback, + analysis_callback, + alive_callback, + stepsize_callback) + +############################################################################### +# run the simulation + +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 diff --git a/src/callbacks_step/stepsize_dg1d.jl b/src/callbacks_step/stepsize_dg1d.jl index d0341b4fc6e..ec3de6efb00 100644 --- a/src/callbacks_step/stepsize_dg1d.jl +++ b/src/callbacks_step/stepsize_dg1d.jl @@ -6,7 +6,8 @@ #! format: noindent function max_dt(u, t, mesh::TreeMesh{1}, - constant_speed::False, equations, dg::DG, cache) + constant_speed::False, equations, + dg::DG, cache) # to avoid a division by zero if the speed vanishes everywhere, # e.g. for steady-state linear advection max_scaled_speed = nextfloat(zero(t)) @@ -26,7 +27,30 @@ function max_dt(u, t, mesh::TreeMesh{1}, end function max_dt(u, t, mesh::TreeMesh{1}, - constant_speed::True, equations, dg::DG, cache) + constant_diffusivity::False, equations, + equations_parabolic::AbstractEquationsParabolic, + dg::DG, cache) + # to avoid a division by zero if the speed vanishes everywhere, + # e.g. for steady-state linear advection + max_scaled_speed = nextfloat(zero(t)) + + for element in eachelement(dg, cache) + max_lambda1 = zero(max_scaled_speed) + for i in eachnode(dg) + u_node = get_node_vars(u, equations, dg, i, element) + lambda1, = max_diffusivity(u_node, equations_parabolic) + max_lambda1 = max(max_lambda1, lambda1) + end + inv_jacobian = cache.elements.inverse_jacobian[element] # Δx + max_scaled_speed = max(max_scaled_speed, inv_jacobian^2 * max_lambda1) + end + + return 2 / (nnodes(dg) * max_scaled_speed) +end + +function max_dt(u, t, mesh::TreeMesh{1}, + constant_speed::True, equations, + dg::DG, cache) # to avoid a division by zero if the speed vanishes everywhere, # e.g. for steady-state linear advection max_scaled_speed = nextfloat(zero(t)) @@ -41,8 +65,9 @@ function max_dt(u, t, mesh::TreeMesh{1}, end function max_dt(u, t, mesh::TreeMesh{1}, - constant_diffusivity::True, - equations_parabolic::AbstractEquationsParabolic, dg::DG, cache) + constant_diffusivity::True, equations, + equations_parabolic::AbstractEquationsParabolic, + dg::DG, cache) # to avoid a division by zero if the speed vanishes everywhere, # e.g. for steady-state linear advection max_scaled_speed = nextfloat(zero(t)) @@ -57,7 +82,8 @@ function max_dt(u, t, mesh::TreeMesh{1}, end function max_dt(u, t, mesh::StructuredMesh{1}, - constant_speed::False, equations, dg::DG, cache) + constant_speed::False, equations, + dg::DG, cache) # to avoid a division by zero if the speed vanishes everywhere, # e.g. for steady-state linear advection max_scaled_speed = nextfloat(zero(t)) @@ -81,7 +107,8 @@ function max_dt(u, t, mesh::StructuredMesh{1}, end function max_dt(u, t, mesh::StructuredMesh{1}, - constant_speed::True, equations, dg::DG, cache) + constant_speed::True, equations, + dg::DG, cache) # to avoid a division by zero if the speed vanishes everywhere, # e.g. for steady-state linear advection max_scaled_speed = nextfloat(zero(t)) @@ -97,4 +124,9 @@ function max_dt(u, t, mesh::StructuredMesh{1}, return 2 / (nnodes(dg) * max_scaled_speed) end + +# Note: `max_dt` is not implemented for `StructuredMesh{1}` since +# for the `StructuredMesh` there is no support of parabolic terms (yet), see the overview in the docs: +# https://trixi-framework.github.io/Trixi.jl/stable/overview/#overview-semidiscretizations + end # @muladd diff --git a/src/callbacks_step/stepsize_hyperbolic_parabolic.jl b/src/callbacks_step/stepsize_hyperbolic_parabolic.jl index 998fc03c2b0..320eca4fe86 100644 --- a/src/callbacks_step/stepsize_hyperbolic_parabolic.jl +++ b/src/callbacks_step/stepsize_hyperbolic_parabolic.jl @@ -8,8 +8,12 @@ """ StepsizeCallbackHyperbolicParabolic(; cfl_convective=1.0, cfl_diffusive=1.0) -Set the time step size according to a CFL condition with CFL number `cfl` -if the time integration method isn't adaptive itself. +Set the time step size according to a CFL condition with +CFL numbers `cfl_convective` for the convective stability estimate and +`cfl_diffusive` for the diffusive stability estimate. +This callback should be used if convection-dominance cannot be assured over the +course of the simulation. +This callback is only relevant for non-adaptive time integrators. """ mutable struct StepsizeCallbackHyperbolicParabolic{RealT} cfl_convective::RealT @@ -100,11 +104,11 @@ function calculate_dt(u_ode, t, cfl_convective, cfl_diffusive, u = wrap_array(u_ode, mesh, equations, solver, cache) dt_convective = cfl_convective * max_dt(u, t, mesh, - have_constant_speed(equations), - equations, solver, cache) + have_constant_speed(equations), equations, + solver, cache) dt_diffusive = cfl_diffusive * max_dt(u, t, mesh, - have_constant_diffusivity(equations_parabolic), + have_constant_diffusivity(equations_parabolic), equations, equations_parabolic, solver, cache) return min(dt_convective, dt_diffusive) @@ -128,11 +132,12 @@ function (cb::DiscreteCallback{Condition, Affect!})(ode::ODEProblem) where {Cond u = wrap_array(u_ode, mesh, equations, solver, cache) dt_convective = dt_convective * - max_dt(u, t, mesh, have_constant_speed(equations), - equations, solver, cache) + max_dt(u, t, mesh, + have_constant_speed(equations), equations, solver, cache) dt_diffusive = dt_diffusive * - max_dt(u, t, mesh, have_constant_diffusivity(equations_parabolic), + max_dt(u, t, mesh, + have_constant_diffusivity(equations_parabolic), equations, equations_parabolic, solver, cache) return min(dt_convective, dt_diffusive) diff --git a/src/equations/compressible_navier_stokes_1d.jl b/src/equations/compressible_navier_stokes_1d.jl index 024ce2d7e87..c0aebf5d2be 100644 --- a/src/equations/compressible_navier_stokes_1d.jl +++ b/src/equations/compressible_navier_stokes_1d.jl @@ -178,6 +178,15 @@ function flux(u, gradients, orientation::Integer, return SVector(f1, f2, f3) end +# The dynamic viscosity `mu` may be a function of temperature +@inline have_constant_diffusivity(::CompressibleNavierStokesDiffusion1D) = False() + +@inline function max_diffusivity(u, + equations_parabolic::CompressibleNavierStokesDiffusion1D) + # Eigenvalues of diffusivity matrix: mu, mu * kappa + return dynamic_viscosity(u, equations_parabolic) * max(1, equations_parabolic.kappa) +end + # Convert conservative variables to primitive @inline function cons2prim(u, equations::CompressibleNavierStokesDiffusion1D) rho, rho_v1, _ = u From 0ce6f12e38abbd8f912ad3702f27c87934babc36 Mon Sep 17 00:00:00 2001 From: Daniel_Doehring Date: Mon, 12 Aug 2024 16:11:49 +0200 Subject: [PATCH 03/17] Comment --- src/equations/compressible_navier_stokes_1d.jl | 15 ++++++++++++++- 1 file changed, 14 insertions(+), 1 deletion(-) diff --git a/src/equations/compressible_navier_stokes_1d.jl b/src/equations/compressible_navier_stokes_1d.jl index c0aebf5d2be..1b3e4a35c9b 100644 --- a/src/equations/compressible_navier_stokes_1d.jl +++ b/src/equations/compressible_navier_stokes_1d.jl @@ -183,7 +183,20 @@ end @inline function max_diffusivity(u, equations_parabolic::CompressibleNavierStokesDiffusion1D) - # Eigenvalues of diffusivity matrix: mu, mu * kappa + # For the diffusive estimate we use the eigenvalues of the diffusivity matrix, as suggested in + + # FLEXI: A high order discontinuous Galerkin framework for hyperbolic–parabolic conservation laws + # https://doi.org/10.1016/j.camwa.2020.05.004 + + # Here, the diffusive flux is given by + # [0, [0, + # mu * τ, = mu * dv/dx, + # mu * (v * τ + q1)] mu * (v * dv/dx + kappa * dT/dx)] + # which can be rewritten as + # [0, 0, 0 + # 0, mu, 0, + # 0, mu * v, mu * q1] * grad(u) . + # Thus, the eigenvalues of the diffusivity matrix are {0, mu, mu * kappa}. return dynamic_viscosity(u, equations_parabolic) * max(1, equations_parabolic.kappa) end From 131b5abd2189a2d7c4697252fab5caac3ba00935 Mon Sep 17 00:00:00 2001 From: Daniel_Doehring Date: Mon, 12 Aug 2024 16:17:20 +0200 Subject: [PATCH 04/17] tests --- ...lixir_navierstokes_convergence_periodic_cfl.jl | 4 +++- test/test_parabolic_1d.jl | 15 +++++++++++++++ 2 files changed, 18 insertions(+), 1 deletion(-) diff --git a/examples/tree_1d_dgsem/elixir_navierstokes_convergence_periodic_cfl.jl b/examples/tree_1d_dgsem/elixir_navierstokes_convergence_periodic_cfl.jl index b3b15ed2ef0..d4ee18608d2 100644 --- a/examples/tree_1d_dgsem/elixir_navierstokes_convergence_periodic_cfl.jl +++ b/examples/tree_1d_dgsem/elixir_navierstokes_convergence_periodic_cfl.jl @@ -122,7 +122,9 @@ analysis_callback = AnalysisCallback(semi, interval = analysis_interval) alive_callback = AliveCallback(analysis_interval = analysis_interval) -# Stepsize callback which selects the timestep to the most restrictive influence +# Stepsize callback which selects the timestep according to the most restrictive CFL condition. +# For coarser grids, linear stability is governed by the convective CFL condition, +# while for high refinements (e.g. initial_refinement_level = 8) the flow becomes diffusion-dominated. stepsize_callback = StepsizeCallbackHyperbolicParabolic(cfl_convective = 1.8, cfl_diffusive = 0.3) diff --git a/test/test_parabolic_1d.jl b/test/test_parabolic_1d.jl index 38bebdcce1d..09ce10a6d7c 100644 --- a/test/test_parabolic_1d.jl +++ b/test/test_parabolic_1d.jl @@ -84,6 +84,21 @@ end end end +@trixi_testset "TreeMesh1D: elixir_navierstokes_convergence_periodic_cfl.jl" begin + @test_trixi_include(joinpath(examples_dir(), "tree_1d_dgsem", + "elixir_navierstokes_convergence_periodic_cfl.jl"), + l2=[0.00011338560756751962, 6.240158271610694e-5, 0.0002848510206540238], + linf=[0.0006233189520368221, 0.0003592942992138859, 0.0016105764529221744]) + # 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 + @trixi_testset "TreeMesh1D: elixir_navierstokes_convergence_periodic.jl: GradientVariablesEntropy" begin @test_trixi_include(joinpath(examples_dir(), "tree_1d_dgsem", "elixir_navierstokes_convergence_periodic.jl"), From 8b538fb9a012e3623c67ca2ca65847c7f4734834 Mon Sep 17 00:00:00 2001 From: Daniel_Doehring Date: Mon, 12 Aug 2024 16:18:03 +0200 Subject: [PATCH 05/17] fmt --- test/test_parabolic_1d.jl | 12 ++++++++++-- 1 file changed, 10 insertions(+), 2 deletions(-) diff --git a/test/test_parabolic_1d.jl b/test/test_parabolic_1d.jl index 09ce10a6d7c..08d704f4b66 100644 --- a/test/test_parabolic_1d.jl +++ b/test/test_parabolic_1d.jl @@ -87,8 +87,16 @@ end @trixi_testset "TreeMesh1D: elixir_navierstokes_convergence_periodic_cfl.jl" begin @test_trixi_include(joinpath(examples_dir(), "tree_1d_dgsem", "elixir_navierstokes_convergence_periodic_cfl.jl"), - l2=[0.00011338560756751962, 6.240158271610694e-5, 0.0002848510206540238], - linf=[0.0006233189520368221, 0.0003592942992138859, 0.0016105764529221744]) + l2=[ + 0.00011338560756751962, + 6.240158271610694e-5, + 0.0002848510206540238, + ], + linf=[ + 0.0006233189520368221, + 0.0003592942992138859, + 0.0016105764529221744, + ]) # Ensure that we do not have excessive memory allocations # (e.g., from type instabilities) let From 55b52b11e179fad8e278008273c2cdfab351b12a Mon Sep 17 00:00:00 2001 From: Daniel_Doehring Date: Mon, 12 Aug 2024 16:21:31 +0200 Subject: [PATCH 06/17] Comment --- src/callbacks_step/stepsize_dg1d.jl | 1 + 1 file changed, 1 insertion(+) diff --git a/src/callbacks_step/stepsize_dg1d.jl b/src/callbacks_step/stepsize_dg1d.jl index ec3de6efb00..8e372fc83b1 100644 --- a/src/callbacks_step/stepsize_dg1d.jl +++ b/src/callbacks_step/stepsize_dg1d.jl @@ -128,5 +128,6 @@ end # Note: `max_dt` is not implemented for `StructuredMesh{1}` since # for the `StructuredMesh` there is no support of parabolic terms (yet), see the overview in the docs: # https://trixi-framework.github.io/Trixi.jl/stable/overview/#overview-semidiscretizations +# Thus, there is also no need to implement the `max_dt` function for the `StructuredMesh{1}` case (yet). end # @muladd From 5c008162022c03b964e75d68b665ab31aa7cfca9 Mon Sep 17 00:00:00 2001 From: Daniel_Doehring Date: Mon, 12 Aug 2024 16:35:04 +0200 Subject: [PATCH 07/17] Bring back adv diff since this has const diff --- .../elixir_advection_diffusion_cfl.jl | 97 +++++++++++++++++++ test/test_parabolic_1d.jl | 16 +++ 2 files changed, 113 insertions(+) create mode 100644 examples/tree_1d_dgsem/elixir_advection_diffusion_cfl.jl diff --git a/examples/tree_1d_dgsem/elixir_advection_diffusion_cfl.jl b/examples/tree_1d_dgsem/elixir_advection_diffusion_cfl.jl new file mode 100644 index 00000000000..0d7335b4218 --- /dev/null +++ b/examples/tree_1d_dgsem/elixir_advection_diffusion_cfl.jl @@ -0,0 +1,97 @@ + +using OrdinaryDiffEq +using Trixi + +############################################################################### +# semidiscretization of the linear advection diffusion equation + +advection_velocity = 0.1 +equations = LinearScalarAdvectionEquation1D(advection_velocity) +diffusivity() = 0.1 +equations_parabolic = LaplaceDiffusion1D(diffusivity(), equations) + +# 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) + +coordinates_min = -pi # minimum coordinate +coordinates_max = pi # maximum coordinate + +# Create a uniformly refined mesh with periodic boundaries +mesh = TreeMesh(coordinates_min, coordinates_max, + initial_refinement_level = 4, + n_cells_max = 30_000, # set maximum capacity of tree data structure + periodicity = true) + +function x_trans_periodic(x, domain_length = SVector(2 * pi), center = SVector(0.0)) + x_normalized = x .- center + x_shifted = x_normalized .% domain_length + x_offset = ((x_shifted .< -0.5 * domain_length) - (x_shifted .> 0.5 * domain_length)) .* + domain_length + return center + x_shifted + x_offset +end + +# Define initial condition +function initial_condition_diffusive_convergence_test(x, t, + equation::LinearScalarAdvectionEquation1D) + # Store translated coordinate for easy use of exact solution + x_trans = x_trans_periodic(x - equation.advection_velocity * t) + + nu = diffusivity() + c = 0.0 + A = 1.0 + L = 2 + f = 1 / L + omega = 1.0 + scalar = c + A * sin(omega * sum(x_trans)) * exp(-nu * omega^2 * t) + return SVector(scalar) +end +initial_condition = initial_condition_diffusive_convergence_test + +# define periodic boundary conditions everywhere +boundary_conditions = boundary_condition_periodic +boundary_conditions_parabolic = boundary_condition_periodic + +# A semidiscretization collects data structures and functions for the spatial discretization +semi = SemidiscretizationHyperbolicParabolic(mesh, (equations, equations_parabolic), + initial_condition, + solver; + boundary_conditions = (boundary_conditions, + boundary_conditions_parabolic)) + +############################################################################### +# 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_callback = AnalysisCallback(semi, interval = 100) + +# The AliveCallback prints short status information in regular intervals +alive_callback = AliveCallback(analysis_interval = 100) + +# Stepsize callback which selects the timestep according to the most restrictive CFL condition. +# For coarser grids, linear stability is governed by the convective CFL condition, +# while for high refinements the flow becomes diffusion-dominated. +stepsize_callback = StepsizeCallbackHyperbolicParabolic(cfl_convective = 1.6, + cfl_diffusive = 0.4) + +# 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, + 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() \ No newline at end of file diff --git a/test/test_parabolic_1d.jl b/test/test_parabolic_1d.jl index 08d704f4b66..d568f712def 100644 --- a/test/test_parabolic_1d.jl +++ b/test/test_parabolic_1d.jl @@ -28,6 +28,22 @@ isdir(outdir) && rm(outdir, recursive = true) end end +@trixi_testset "TreeMesh1D: elixir_advection_diffusion_cfl.jl" begin + @test_trixi_include(joinpath(examples_dir(), "tree_1d_dgsem", + "elixir_advection_diffusion_cfl.jl"), + initial_refinement_level=4, tspan=(0.0, 0.4), polydeg=3, + l2=[1.0840828945107538e-5], + linf=[3.911865515410229e-5]) + # 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 + @trixi_testset "TreeMesh1D: elixir_advection_diffusion.jl (AMR)" begin @test_trixi_include(joinpath(examples_dir(), "tree_1d_dgsem", "elixir_advection_diffusion.jl"), From b67b7c610b3a439fb35ee79c25aab8d886ecabfe Mon Sep 17 00:00:00 2001 From: Daniel_Doehring Date: Mon, 12 Aug 2024 16:41:15 +0200 Subject: [PATCH 08/17] fmt --- examples/tree_1d_dgsem/elixir_advection_diffusion_cfl.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/examples/tree_1d_dgsem/elixir_advection_diffusion_cfl.jl b/examples/tree_1d_dgsem/elixir_advection_diffusion_cfl.jl index 0d7335b4218..0c58b589eb4 100644 --- a/examples/tree_1d_dgsem/elixir_advection_diffusion_cfl.jl +++ b/examples/tree_1d_dgsem/elixir_advection_diffusion_cfl.jl @@ -94,4 +94,4 @@ sol = solve(ode, CarpenterKennedy2N54(williamson_condition = false), save_everystep = false, callback = callbacks); # Print the timer summary -summary_callback() \ No newline at end of file +summary_callback() From 4fb09f0e3c937326d212e3772a1432a38bc69e19 Mon Sep 17 00:00:00 2001 From: Daniel_Doehring Date: Tue, 13 Aug 2024 09:12:04 +0200 Subject: [PATCH 09/17] test values for t = 0.4 --- test/test_parabolic_1d.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/test/test_parabolic_1d.jl b/test/test_parabolic_1d.jl index d568f712def..b48ddae7ce1 100644 --- a/test/test_parabolic_1d.jl +++ b/test/test_parabolic_1d.jl @@ -32,8 +32,8 @@ end @test_trixi_include(joinpath(examples_dir(), "tree_1d_dgsem", "elixir_advection_diffusion_cfl.jl"), initial_refinement_level=4, tspan=(0.0, 0.4), polydeg=3, - l2=[1.0840828945107538e-5], - linf=[3.911865515410229e-5]) + l2=[8.389469356681668e-6], + linf=[2.8474476202633436e-5]) # Ensure that we do not have excessive memory allocations # (e.g., from type instabilities) let From dd78bc7d26a846b021413e9f78ce01a6f401b266 Mon Sep 17 00:00:00 2001 From: Daniel_Doehring Date: Tue, 13 Aug 2024 17:01:56 +0200 Subject: [PATCH 10/17] use only one stepsize method --- .../elixir_advection_diffusion_cfl.jl | 2 +- src/Trixi.jl | 2 +- src/callbacks_step/callbacks_step.jl | 1 - src/callbacks_step/stepsize.jl | 84 ++++++++-- .../stepsize_hyperbolic_parabolic.jl | 145 ------------------ 5 files changed, 71 insertions(+), 163 deletions(-) delete mode 100644 src/callbacks_step/stepsize_hyperbolic_parabolic.jl diff --git a/examples/tree_1d_dgsem/elixir_advection_diffusion_cfl.jl b/examples/tree_1d_dgsem/elixir_advection_diffusion_cfl.jl index 0c58b589eb4..3c51b72f6ad 100644 --- a/examples/tree_1d_dgsem/elixir_advection_diffusion_cfl.jl +++ b/examples/tree_1d_dgsem/elixir_advection_diffusion_cfl.jl @@ -78,7 +78,7 @@ alive_callback = AliveCallback(analysis_interval = 100) # Stepsize callback which selects the timestep according to the most restrictive CFL condition. # For coarser grids, linear stability is governed by the convective CFL condition, # while for high refinements the flow becomes diffusion-dominated. -stepsize_callback = StepsizeCallbackHyperbolicParabolic(cfl_convective = 1.6, +stepsize_callback = StepsizeCallback(cfl = 1.6, cfl_diffusive = 0.4) # Create a CallbackSet to collect all callbacks such that they can be passed to the ODE solver diff --git a/src/Trixi.jl b/src/Trixi.jl index 623df916f8e..23a8cfe1d0a 100644 --- a/src/Trixi.jl +++ b/src/Trixi.jl @@ -262,7 +262,7 @@ export SemidiscretizationCoupled export SummaryCallback, SteadyStateCallback, AnalysisCallback, AliveCallback, SaveRestartCallback, SaveSolutionCallback, TimeSeriesCallback, VisualizationCallback, AveragingCallback, - AMRCallback, StepsizeCallback, StepsizeCallbackHyperbolicParabolic, + AMRCallback, StepsizeCallback, GlmSpeedCallback, LBMCollisionCallback, EulerAcousticsCouplingCallback, TrivialCallback, AnalysisCallbackCoupled, AnalysisSurfaceIntegral, DragCoefficientPressure, LiftCoefficientPressure, diff --git a/src/callbacks_step/callbacks_step.jl b/src/callbacks_step/callbacks_step.jl index 2b050acca5c..09d197bf225 100644 --- a/src/callbacks_step/callbacks_step.jl +++ b/src/callbacks_step/callbacks_step.jl @@ -61,7 +61,6 @@ include("averaging.jl") include("amr.jl") include("stepsize.jl") -include("stepsize_hyperbolic_parabolic.jl") include("glm_speed.jl") include("lbm_collision.jl") include("euler_acoustics_coupling.jl") diff --git a/src/callbacks_step/stepsize.jl b/src/callbacks_step/stepsize.jl index 8b5cb958318..352b57c443c 100644 --- a/src/callbacks_step/stepsize.jl +++ b/src/callbacks_step/stepsize.jl @@ -6,21 +6,32 @@ #! format: noindent """ - StepsizeCallback(; cfl=1.0) + StepsizeCallback(; cfl=1.0, cfl_diffusive = 0.0) Set the time step size according to a CFL condition with CFL number `cfl` if the time integration method isn't adaptive itself. +One can additionally supply a diffusive CFL number `cfl_diffusive` to +limit the admissible timestep also respecting diffusive restrictions. +In that case, a number larger than zero needs to be supplied. +By default, `cfl_diffusive` is set to zero which means that only the convective +CFL number is considered. """ mutable struct StepsizeCallback{RealT} - cfl_number::RealT + cfl_convective::RealT + cfl_diffusive::RealT end function Base.show(io::IO, cb::DiscreteCallback{<:Any, <:StepsizeCallback}) @nospecialize cb # reduce precompilation time stepsize_callback = cb.affect! - @unpack cfl_number = stepsize_callback - print(io, "StepsizeCallback(cfl_number=", cfl_number, ")") + @unpack cfl_convective, cfl_diffusive = stepsize_callback + if cfl_diffusive == 0.0 + print(io, "StepsizeCallback(cfl_convective=", cfl_convective, ")") + else + print(io, "StepsizeCallback(cfl_convective=", cfl_convective, + "cfl_diffusive=", cfl_diffusive, ")") + end end function Base.show(io::IO, ::MIME"text/plain", @@ -32,15 +43,22 @@ function Base.show(io::IO, ::MIME"text/plain", else stepsize_callback = cb.affect! - setup = [ - "CFL number" => stepsize_callback.cfl_number, - ] + if stepsize_callback.cfl_diffusive == 0.0 + setup = [ + "CFL number" => stepsize_callback.cfl_convective, + ] + else + setup = [ + "CFL number" => stepsize_callback.cfl_convective, + "Diffusive CFL number" => stepsize_callback.cfl_diffusive, + ] + end summary_box(io, "StepsizeCallback", setup) end end -function StepsizeCallback(; cfl::Real = 1.0) - stepsize_callback = StepsizeCallback(cfl) +function StepsizeCallback(; cfl::Real = 1.0, cfl_diffusive::Real = 0.0) + stepsize_callback = StepsizeCallback(cfl, cfl_diffusive) DiscreteCallback(stepsize_callback, stepsize_callback, # the first one is the condition, the second the affect! save_positions = (false, false), @@ -64,10 +82,11 @@ end t = integrator.t u_ode = integrator.u semi = integrator.p - @unpack cfl_number = stepsize_callback + @unpack cfl_convective, cfl_diffusive = stepsize_callback # Dispatch based on semidiscretization - dt = @trixi_timeit timer() "calculate dt" calculate_dt(u_ode, t, cfl_number, + dt = @trixi_timeit timer() "calculate dt" calculate_dt(u_ode, t, + cfl_convective, cfl_diffusive, semi) set_proposed_dt!(integrator, dt) @@ -81,15 +100,38 @@ end end # General case for a single semidiscretization -function calculate_dt(u_ode, t, cfl_number, semi::AbstractSemidiscretization) +function calculate_dt(u_ode, t, cfl_convective, cfl_diffusive, 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, + dt = cfl_convective * max_dt(u, t, mesh, have_constant_speed(equations), equations, solver, cache) end +# Case for a hyperbolic-parabolic semidiscretization +function calculate_dt(u_ode, t, cfl_convective, cfl_diffusive, + semi::SemidiscretizationHyperbolicParabolic) + mesh, equations, solver, cache = mesh_equations_solver_cache(semi) + equations_parabolic = semi.equations_parabolic + + u = wrap_array(u_ode, mesh, equations, solver, cache) + + dt_convective = cfl_convective * max_dt(u, t, mesh, + have_constant_speed(equations), equations, + solver, cache) + + if cfl_diffusive > 0.0 # Check if diffusive CFL should be considered + dt_diffusive = cfl_diffusive * max_dt(u, t, mesh, + have_constant_diffusivity(equations_parabolic), equations, + equations_parabolic, solver, cache) + + return min(dt_convective, dt_diffusive) + else + return dt_convective + end +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 @@ -99,15 +141,27 @@ function (cb::DiscreteCallback{Condition, Affect!})(ode::ODEProblem) where {Cond StepsizeCallback } stepsize_callback = cb.affect! - @unpack cfl_number = stepsize_callback + @unpack cfl_convective, cfl_diffusive = stepsize_callback u_ode = ode.u0 t = first(ode.tspan) semi = ode.p mesh, equations, solver, cache = mesh_equations_solver_cache(semi) u = wrap_array(u_ode, mesh, equations, solver, cache) - return cfl_number * + dt_convective = cfl_convective * max_dt(u, t, mesh, have_constant_speed(equations), equations, solver, cache) + + if hasfield(semi, :equations_parabolic) # Check if we have a hyperbolic-parabolic semidiscretization + if cfl_diffusive > 0.0 # Check if diffusive CFL should be considered + dt_diffusive = cfl_diffusive * + max_dt(u, t, mesh, have_constant_diffusivity(semi.equations_parabolic), + equations, semi.equations_parabolic, solver, cache) + + return min(dt_convective, dt_diffusive) + else + return dt_convective + end + end end include("stepsize_dg1d.jl") diff --git a/src/callbacks_step/stepsize_hyperbolic_parabolic.jl b/src/callbacks_step/stepsize_hyperbolic_parabolic.jl deleted file mode 100644 index 320eca4fe86..00000000000 --- a/src/callbacks_step/stepsize_hyperbolic_parabolic.jl +++ /dev/null @@ -1,145 +0,0 @@ -# By default, Julia/LLVM does not use fused multiply-add operations (FMAs). -# Since these FMAs can increase the performance of many numerical algorithms, -# we need to opt-in explicitly. -# See https://ranocha.de/blog/Optimizing_EC_Trixi for further details. -@muladd begin -#! format: noindent - -""" - StepsizeCallbackHyperbolicParabolic(; cfl_convective=1.0, cfl_diffusive=1.0) - -Set the time step size according to a CFL condition with -CFL numbers `cfl_convective` for the convective stability estimate and -`cfl_diffusive` for the diffusive stability estimate. -This callback should be used if convection-dominance cannot be assured over the -course of the simulation. -This callback is only relevant for non-adaptive time integrators. -""" -mutable struct StepsizeCallbackHyperbolicParabolic{RealT} - cfl_convective::RealT - cfl_diffusive::RealT -end - -function Base.show(io::IO, - cb::DiscreteCallback{<:Any, <:StepsizeCallbackHyperbolicParabolic}) - @nospecialize cb # reduce precompilation time - - stepsize_callback = cb.affect! - @unpack cfl_convective = stepsize_callback - @unpack cfl_diffusive = stepsize_callback - print(io, "StepsizeCallbackHyperbolicParabolic(cfl_convective=", cfl_convective, - "cfl_diffusive=", cfl_diffusive, ")") -end - -function Base.show(io::IO, ::MIME"text/plain", - cb::DiscreteCallback{<:Any, <:StepsizeCallbackHyperbolicParabolic}) - @nospecialize cb # reduce precompilation time - - if get(io, :compact, false) - show(io, cb) - else - stepsize_callback = cb.affect! - - setup = [ - "Convective CFL number" => stepsize_callback.cfl_convective, - "Diffusive CFL number" => stepsize_callback.cfl_diffusive, - ] - summary_box(io, "StepsizeCallbackHyperbolicParabolic", setup) - end -end - -function StepsizeCallbackHyperbolicParabolic(; - cfl_convective::Real = 1.0, - cfl_diffusive::Real = 1.0) - stepsize_callback = StepsizeCallbackHyperbolicParabolic(cfl_convective, - cfl_diffusive) - - DiscreteCallback(stepsize_callback, stepsize_callback, # the first one is the condition, the second the affect! - save_positions = (false, false), - initialize = initialize!) -end - -function initialize!(cb::DiscreteCallback{Condition, Affect!}, u, t, - integrator) where {Condition, - Affect! <: StepsizeCallbackHyperbolicParabolic} - cb.affect!(integrator) -end - -# this method is called to determine whether the callback should be activated -function (stepsize_callback::StepsizeCallbackHyperbolicParabolic)(u, t, integrator) - return true -end - -# This method is called as callback during the time integration. -@inline function (stepsize_callback::StepsizeCallbackHyperbolicParabolic)(integrator) - # TODO: Taal decide, shall we set the time step even if the integrator is adaptive? - if !integrator.opts.adaptive - t = integrator.t - u_ode = integrator.u - semi = integrator.p - @unpack cfl_convective, cfl_diffusive = stepsize_callback - - # Dispatch based on semidiscretization - dt = @trixi_timeit timer() "calculate dt" calculate_dt(u_ode, t, - cfl_convective, - cfl_diffusive, - semi) - - set_proposed_dt!(integrator, dt) - integrator.opts.dtmax = dt - integrator.dtcache = dt - end - - # avoid re-evaluating possible FSAL stages - u_modified!(integrator, false) - return nothing -end - -# Case for a hyperbolic-parabolic semidiscretization -function calculate_dt(u_ode, t, cfl_convective, cfl_diffusive, - semi::SemidiscretizationHyperbolicParabolic) - mesh, equations, solver, cache = mesh_equations_solver_cache(semi) - equations_parabolic = semi.equations_parabolic - - u = wrap_array(u_ode, mesh, equations, solver, cache) - - dt_convective = cfl_convective * max_dt(u, t, mesh, - have_constant_speed(equations), equations, - solver, cache) - - dt_diffusive = cfl_diffusive * max_dt(u, t, mesh, - have_constant_diffusivity(equations_parabolic), equations, - equations_parabolic, solver, cache) - - return min(dt_convective, dt_diffusive) -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 -# many examples in `solve(ode, ..., dt=stepsize_callback(ode), ...)`. -function (cb::DiscreteCallback{Condition, Affect!})(ode::ODEProblem) where {Condition, - Affect! <: - StepsizeCallbackHyperbolicParabolic - } - stepsize_callback = cb.affect! - @unpack cfl_convective, cfl_diffusive = stepsize_callback - u_ode = ode.u0 - t = first(ode.tspan) - semi = ode.p - mesh, equations, solver, cache = mesh_equations_solver_cache(semi) - equations_parabolic = semi.equations_parabolic - u = wrap_array(u_ode, mesh, equations, solver, cache) - - dt_convective = dt_convective * - max_dt(u, t, mesh, - have_constant_speed(equations), equations, solver, cache) - - dt_diffusive = dt_diffusive * - max_dt(u, t, mesh, - have_constant_diffusivity(equations_parabolic), equations, - equations_parabolic, solver, cache) - - return min(dt_convective, dt_diffusive) -end -end # @muladd From bd76a7e595436915f2c7b13b04e21af4a573a782 Mon Sep 17 00:00:00 2001 From: Daniel_Doehring Date: Tue, 13 Aug 2024 17:02:49 +0200 Subject: [PATCH 11/17] fmt --- src/callbacks_step/stepsize.jl | 20 ++++++++++++-------- 1 file changed, 12 insertions(+), 8 deletions(-) diff --git a/src/callbacks_step/stepsize.jl b/src/callbacks_step/stepsize.jl index 352b57c443c..21347d7e458 100644 --- a/src/callbacks_step/stepsize.jl +++ b/src/callbacks_step/stepsize.jl @@ -85,8 +85,9 @@ end @unpack cfl_convective, cfl_diffusive = stepsize_callback # Dispatch based on semidiscretization - dt = @trixi_timeit timer() "calculate dt" calculate_dt(u_ode, t, - cfl_convective, cfl_diffusive, + dt = @trixi_timeit timer() "calculate dt" calculate_dt(u_ode, t, + cfl_convective, + cfl_diffusive, semi) set_proposed_dt!(integrator, dt) @@ -100,7 +101,8 @@ end end # General case for a single semidiscretization -function calculate_dt(u_ode, t, cfl_convective, cfl_diffusive, semi::AbstractSemidiscretization) +function calculate_dt(u_ode, t, cfl_convective, cfl_diffusive, + semi::AbstractSemidiscretization) mesh, equations, solver, cache = mesh_equations_solver_cache(semi) u = wrap_array(u_ode, mesh, equations, solver, cache) @@ -123,8 +125,8 @@ function calculate_dt(u_ode, t, cfl_convective, cfl_diffusive, if cfl_diffusive > 0.0 # Check if diffusive CFL should be considered dt_diffusive = cfl_diffusive * max_dt(u, t, mesh, - have_constant_diffusivity(equations_parabolic), equations, - equations_parabolic, solver, cache) + have_constant_diffusivity(equations_parabolic), equations, + equations_parabolic, solver, cache) return min(dt_convective, dt_diffusive) else @@ -149,13 +151,15 @@ function (cb::DiscreteCallback{Condition, Affect!})(ode::ODEProblem) where {Cond u = wrap_array(u_ode, mesh, equations, solver, cache) dt_convective = cfl_convective * - max_dt(u, t, mesh, have_constant_speed(equations), equations, solver, cache) + max_dt(u, t, mesh, have_constant_speed(equations), equations, + solver, cache) if hasfield(semi, :equations_parabolic) # Check if we have a hyperbolic-parabolic semidiscretization if cfl_diffusive > 0.0 # Check if diffusive CFL should be considered dt_diffusive = cfl_diffusive * - max_dt(u, t, mesh, have_constant_diffusivity(semi.equations_parabolic), - equations, semi.equations_parabolic, solver, cache) + max_dt(u, t, mesh, + have_constant_diffusivity(semi.equations_parabolic), + equations, semi.equations_parabolic, solver, cache) return min(dt_convective, dt_diffusive) else From a677363f83d90eea73d63d3054283494fbdb496f Mon Sep 17 00:00:00 2001 From: Daniel_Doehring Date: Tue, 13 Aug 2024 17:05:17 +0200 Subject: [PATCH 12/17] else case --- src/callbacks_step/stepsize.jl | 2 ++ 1 file changed, 2 insertions(+) diff --git a/src/callbacks_step/stepsize.jl b/src/callbacks_step/stepsize.jl index 21347d7e458..4f22631d89f 100644 --- a/src/callbacks_step/stepsize.jl +++ b/src/callbacks_step/stepsize.jl @@ -165,6 +165,8 @@ function (cb::DiscreteCallback{Condition, Affect!})(ode::ODEProblem) where {Cond else return dt_convective end + else + return dt_convective end end From acdbd20799f4ab8abd40a3115d44d65fd6ef3dc8 Mon Sep 17 00:00:00 2001 From: Daniel_Doehring Date: Tue, 13 Aug 2024 17:10:01 +0200 Subject: [PATCH 13/17] fmt --- examples/tree_1d_dgsem/elixir_advection_diffusion_cfl.jl | 2 +- .../elixir_navierstokes_convergence_periodic_cfl.jl | 4 ++-- 2 files changed, 3 insertions(+), 3 deletions(-) diff --git a/examples/tree_1d_dgsem/elixir_advection_diffusion_cfl.jl b/examples/tree_1d_dgsem/elixir_advection_diffusion_cfl.jl index 3c51b72f6ad..c96d20e7b03 100644 --- a/examples/tree_1d_dgsem/elixir_advection_diffusion_cfl.jl +++ b/examples/tree_1d_dgsem/elixir_advection_diffusion_cfl.jl @@ -79,7 +79,7 @@ alive_callback = AliveCallback(analysis_interval = 100) # For coarser grids, linear stability is governed by the convective CFL condition, # while for high refinements the flow becomes diffusion-dominated. stepsize_callback = StepsizeCallback(cfl = 1.6, - cfl_diffusive = 0.4) + cfl_diffusive = 0.4) # 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, diff --git a/examples/tree_1d_dgsem/elixir_navierstokes_convergence_periodic_cfl.jl b/examples/tree_1d_dgsem/elixir_navierstokes_convergence_periodic_cfl.jl index d4ee18608d2..5fccffb688c 100644 --- a/examples/tree_1d_dgsem/elixir_navierstokes_convergence_periodic_cfl.jl +++ b/examples/tree_1d_dgsem/elixir_navierstokes_convergence_periodic_cfl.jl @@ -125,8 +125,8 @@ alive_callback = AliveCallback(analysis_interval = analysis_interval) # Stepsize callback which selects the timestep according to the most restrictive CFL condition. # For coarser grids, linear stability is governed by the convective CFL condition, # while for high refinements (e.g. initial_refinement_level = 8) the flow becomes diffusion-dominated. -stepsize_callback = StepsizeCallbackHyperbolicParabolic(cfl_convective = 1.8, - cfl_diffusive = 0.3) +stepsize_callback = StepsizeCallback(cfl = 1.8, + cfl_diffusive = 0.3) callbacks = CallbackSet(summary_callback, analysis_callback, From a3bd884a03e87b4a1152566aa44bf73a00ae2ccc Mon Sep 17 00:00:00 2001 From: Daniel_Doehring Date: Tue, 13 Aug 2024 17:36:35 +0200 Subject: [PATCH 14/17] remove hasfield query --- src/callbacks_step/stepsize.jl | 22 +++++++++++----------- 1 file changed, 11 insertions(+), 11 deletions(-) diff --git a/src/callbacks_step/stepsize.jl b/src/callbacks_step/stepsize.jl index 4f22631d89f..0783e56586e 100644 --- a/src/callbacks_step/stepsize.jl +++ b/src/callbacks_step/stepsize.jl @@ -154,20 +154,20 @@ function (cb::DiscreteCallback{Condition, Affect!})(ode::ODEProblem) where {Cond max_dt(u, t, mesh, have_constant_speed(equations), equations, solver, cache) - if hasfield(semi, :equations_parabolic) # Check if we have a hyperbolic-parabolic semidiscretization - if cfl_diffusive > 0.0 # Check if diffusive CFL should be considered - dt_diffusive = cfl_diffusive * - max_dt(u, t, mesh, - have_constant_diffusivity(semi.equations_parabolic), - equations, semi.equations_parabolic, solver, cache) - - return min(dt_convective, dt_diffusive) - else - return dt_convective - end + #if hasfield(semi, :equations_parabolic) # Check if we have a hyperbolic-parabolic semidiscretization + if cfl_diffusive > 0.0 # Check if diffusive CFL should be considered + dt_diffusive = cfl_diffusive * + max_dt(u, t, mesh, + have_constant_diffusivity(semi.equations_parabolic), + equations, semi.equations_parabolic, solver, cache) + + return min(dt_convective, dt_diffusive) else return dt_convective end + #else + # return dt_convective + #end end include("stepsize_dg1d.jl") From 8a8e6d6480bfc027bda9063be025afc2b9d05c3d Mon Sep 17 00:00:00 2001 From: Daniel_Doehring Date: Tue, 13 Aug 2024 17:40:13 +0200 Subject: [PATCH 15/17] try kw --- src/callbacks_step/euler_acoustics_coupling.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/callbacks_step/euler_acoustics_coupling.jl b/src/callbacks_step/euler_acoustics_coupling.jl index 52dc55befdc..7c7f8c07e3e 100644 --- a/src/callbacks_step/euler_acoustics_coupling.jl +++ b/src/callbacks_step/euler_acoustics_coupling.jl @@ -130,8 +130,8 @@ function EulerAcousticsCouplingCallback(ode_euler, mean_values, alg, cfl_acousti euler_acoustics_coupling = EulerAcousticsCouplingCallback{typeof(cfl_acoustics), typeof(mean_values), - typeof(integrator_euler)}(StepsizeCallback(cfl_acoustics), - StepsizeCallback(cfl_euler), + typeof(integrator_euler)}(StepsizeCallback(cfl = cfl_acoustics), + StepsizeCallback(cfl = cfl_euler), mean_values, integrator_euler) condition = (u, t, integrator) -> true From 9bd9a700160e62c6b42ca28fb6ea29b9b053e3f4 Mon Sep 17 00:00:00 2001 From: Daniel_Doehring Date: Tue, 13 Aug 2024 17:45:12 +0200 Subject: [PATCH 16/17] no kw --- src/callbacks_step/euler_acoustics_coupling.jl | 6 ++++-- 1 file changed, 4 insertions(+), 2 deletions(-) diff --git a/src/callbacks_step/euler_acoustics_coupling.jl b/src/callbacks_step/euler_acoustics_coupling.jl index 7c7f8c07e3e..c3c6f2e681a 100644 --- a/src/callbacks_step/euler_acoustics_coupling.jl +++ b/src/callbacks_step/euler_acoustics_coupling.jl @@ -130,8 +130,10 @@ function EulerAcousticsCouplingCallback(ode_euler, mean_values, alg, cfl_acousti euler_acoustics_coupling = EulerAcousticsCouplingCallback{typeof(cfl_acoustics), typeof(mean_values), - typeof(integrator_euler)}(StepsizeCallback(cfl = cfl_acoustics), - StepsizeCallback(cfl = cfl_euler), + typeof(integrator_euler)}(StepsizeCallback(cfl_acoustics, + 0.0), + StepsizeCallback(cfl_euler, + 0.0), mean_values, integrator_euler) condition = (u, t, integrator) -> true From 4ac4970557916c073f43e509f7a574779772b3d0 Mon Sep 17 00:00:00 2001 From: Daniel_Doehring Date: Tue, 13 Aug 2024 18:23:51 +0200 Subject: [PATCH 17/17] cfl for coupled systems --- src/callbacks_step/stepsize.jl | 9 ++++----- src/semidiscretization/semidiscretization_coupled.jl | 5 +++-- 2 files changed, 7 insertions(+), 7 deletions(-) diff --git a/src/callbacks_step/stepsize.jl b/src/callbacks_step/stepsize.jl index 0783e56586e..1ee03110225 100644 --- a/src/callbacks_step/stepsize.jl +++ b/src/callbacks_step/stepsize.jl @@ -154,8 +154,10 @@ function (cb::DiscreteCallback{Condition, Affect!})(ode::ODEProblem) where {Cond max_dt(u, t, mesh, have_constant_speed(equations), equations, solver, cache) - #if hasfield(semi, :equations_parabolic) # Check if we have a hyperbolic-parabolic semidiscretization - if cfl_diffusive > 0.0 # Check if diffusive CFL should be considered + # Check if diffusive CFL should be considered. + # NOTE: + # For non-zero `cfl_diffusive`, `semi` is expected to be a `SemidiscretizationHyperbolicParabolic`. + if cfl_diffusive > 0.0 dt_diffusive = cfl_diffusive * max_dt(u, t, mesh, have_constant_diffusivity(semi.equations_parabolic), @@ -165,9 +167,6 @@ function (cb::DiscreteCallback{Condition, Affect!})(ode::ODEProblem) where {Cond else return dt_convective end - #else - # return dt_convective - #end end include("stepsize_dg1d.jl") diff --git a/src/semidiscretization/semidiscretization_coupled.jl b/src/semidiscretization/semidiscretization_coupled.jl index 6b009cfad20..a02f545edee 100644 --- a/src/semidiscretization/semidiscretization_coupled.jl +++ b/src/semidiscretization/semidiscretization_coupled.jl @@ -350,10 +350,11 @@ end ################################################################################ # In case of coupled system, use minimum timestep over all systems -function calculate_dt(u_ode, t, cfl_number, semi::SemidiscretizationCoupled) +function calculate_dt(u_ode, t, cfl_convective, cfl_diffusive, + semi::SemidiscretizationCoupled) dt = minimum(eachsystem(semi)) do i u_ode_slice = get_system_u_ode(u_ode, i, semi) - calculate_dt(u_ode_slice, t, cfl_number, semi.semis[i]) + calculate_dt(u_ode_slice, t, cfl_convective, cfl_diffusive, semi.semis[i]) end return dt