Viscous CFL: TreeMesh{1}

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 = 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
+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()

examples/tree_1d_dgsem/elixir_navierstokes_convergence_periodic_cfl.jl

@@ -0,0 +1,143 @@
+
+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 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 = StepsizeCallback(cfl = 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

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_acoustics),
-                                                                                        StepsizeCallback(cfl_euler),
+                                                              typeof(integrator_euler)}(StepsizeCallback(cfl_acoustics,
+                                                                                                         0.0),
+                                                                                        StepsizeCallback(cfl_euler,
+                                                                                                         0.0),
                                                                                         mean_values,
                                                                                         integrator_euler)
     condition = (u, t, integrator) -> true