From 95518c5670774dfccfd66db2fc0df15ec91b251f Mon Sep 17 00:00:00 2001 From: Simon Candelaresi <10759273+SimonCan@users.noreply.github.com> Date: Fri, 16 Jun 2023 20:35:24 +0100 Subject: [PATCH 01/15] Initial support for surface coupling of two systems (#1452) * Corrected bugs from coupled to main merger. * Added polytropic equation. * Added further polytropic equations and examples. * Added coupling equations. * Added coupling equation between Euler. * Commented debugging bits, like infiltrator. * Add missing `using` * Fix `destats` deprecation warning * Added coupled elixir. * Removed commented testing code. * Added other_list to the coupled semi discretisation elixir. * Removed flux coupling equation. * Removed surface coupling equation for polytropic Euler. * Removed polytropic Euler equation. * Removed any code related to BoundaryConditionCoupledAB. * Removed flux coupling code. * Removed numerical fluxes for BoundaryConditionCoupledAB. * Removed surface fluxes for BoundaryConditionCoupledAB. * Removed coupled elixir. * Remove Coupled StructuredMesh from visualization test. * Remove duplicate function definitions * make advection elixir go further * Removed initial_condition_peak. * Removed src/equations/hyperbolic_diffusion_2d.jl. * Removed 3d coupling. * Remopved 3d capability. * Removed 3d capability. * Removed 3d plotting of coupled data. * Remove extra dependencies * Remove whitespace changes * Fix type instability * Some temporary fixes. * Fix type in semidiscretization * Removed analysis_callback for simple coupled elixir. * Removed analysis callbacks for the coupled case. * Removed AnalysysCallback for coupled elixir. * Removed polytropic coupling elixir. * Update summary output * Update src/solvers/dgsem_structured/dg_2d.jl * Format summary * Fix save solution callback * Remove unused code * Move timeit call before dispatch on semi * Avoid copy on calculcate_dt * Avoid copy on save_solution_file * Remove unnnecessary override of wrap_array * Undo changes to analysis callback * Remove equations_list * Remove unused functions * nmeshes -> nsystems * Further cleanup * Move BoundaryConditionCoupled to the correct location * Visualize setup * Formatting improvmenets * Change 1:nsystems(semi) to eachsystem(semi) * Remove redundant ndofs(...) function * copy_to_coupled_boundary --> copy_to_coupled_boundary! * Move all SemidiscretizationCoupled-specific code to semi/semi_coupled.jl * Use uEltype for BCCoupled * Add comment * I --> Indices * Add comment * Remove Base.summary for SemiCoupled since it appears to be unused * Add parens * Int64 -> Int * Add xref for Documenter.jl * Fixup comment * Remove unused `total_volume` * Remove obsolete comment * summary_semi --> print_summary_semi for clarity * Make SemiCoupled ctor more convenient * Fix docstring * Add description to elixir * Rename elixir * Remove unused kwarg * Fix argument order and simplify interface for IO functions * Explicitly return nothing in functions that should do - nothing * Update comment * Add AnalysisCallback to coupled semidiscretization (#1505) * Add AnalysisCallback to coupled semidiscretization * First non-crashing version of the AnalysisCallbackCoupled * Added comment to offending line in the analysis callback. * Fix stupid bug * Rename variable for easier testing * Clean up code * Remove type instability * Prettify output * Add test * Enable `convergence_test` for SemidiscretizationCoupled * Increased the frequency of the solution write out for a more usable animation. * Reverted analysis intervals. * Updated the l2 and linf errors for the elixir_advection_basic_coupled.jl test to reflect the increased simulation time. * Corrected bracket typo in structured_2d test. * Renamed plural variable names to lists. * Further renaiming plural variable names. * Added convergence_test for elixir_advection_basic_coupled. * Fix coverage for convergence_test * Add test for analysis_callback(sol) * Split timers between systems * fully specialize on rhs and parameters when constructing an ODEProblem * switch example back to OrdinaryDiffEq.jl * Reverted coupled example to use Trixi.solve instead of OrdinaryDiffEq solve. This should fix issues with LoadError: Failed to precompile OrdinaryDiffEq in the thread_legacy test. * Changed Julia version in project toml to 1.9 to fix OrdinaryDiffEq issues on the github test. * Change 1:nsystems(semi) to eachsystem(semi) * Use `get_system_u_ode` * Move all SemidiscretizationCoupled-specific code to semi/semi_coupled.jl * Changed file name name of elixir_advection_basic_coupled.jl to elixir_advection_coupled.jl. * Reverted Julia version to 1.8 in Project toml file. * Apply suggestions from code review * -use type SciMLBase.FullSpecialize instead of instance * Use get_system_u_ode instead of manual view * Reorder elixir ingredients * Make comment reflect code again * Use solve from OrdinaryDiffEq * Use more precise type for array * Test EOCs for each system separately * Allow test to run for the full duration --------- Co-authored-by: SimonCan Co-authored-by: Hendrik Ranocha * Remove unused `total_volume(...)` * Make `save_solution_file` work for SemiEulerGravity again (and make it multi-system aware) * Update src/semidiscretization/semidiscretization_euler_gravity.jl Co-authored-by: Hendrik Ranocha * Apply formatting --------- Co-authored-by: Michael Schlottke-Lakemper Co-authored-by: Hendrik Ranocha Co-authored-by: Hendrik Ranocha Co-authored-by: Benjamin Bolm <74359358+bennibolm@users.noreply.github.com> Co-authored-by: Erik Faulhaber <44124897+efaulhaber@users.noreply.github.com> Co-authored-by: Lucas Gemein <74359570+NichtLucas@users.noreply.github.com> --- .../elixir_advection_coupled.jl | 117 ++++ src/Trixi.jl | 8 +- src/auxiliary/special_elixirs.jl | 16 +- src/callbacks_step/analysis.jl | 30 +- src/callbacks_step/save_solution.jl | 91 +-- src/callbacks_step/stepsize.jl | 19 +- src/callbacks_step/summary.jl | 22 +- src/meshes/mesh_io.jl | 8 +- src/semidiscretization/semidiscretization.jl | 6 +- .../semidiscretization_coupled.jl | 610 ++++++++++++++++++ .../semidiscretization_euler_gravity.jl | 17 +- test/test_special_elixirs.jl | 7 + test/test_structured_2d.jl | 13 + 13 files changed, 890 insertions(+), 74 deletions(-) create mode 100644 examples/structured_2d_dgsem/elixir_advection_coupled.jl create mode 100644 src/semidiscretization/semidiscretization_coupled.jl diff --git a/examples/structured_2d_dgsem/elixir_advection_coupled.jl b/examples/structured_2d_dgsem/elixir_advection_coupled.jl new file mode 100644 index 00000000000..1e54e411db6 --- /dev/null +++ b/examples/structured_2d_dgsem/elixir_advection_coupled.jl @@ -0,0 +1,117 @@ +using OrdinaryDiffEq +using Trixi + + +############################################################################### +# Coupled semidiscretization of two linear advection systems, which are connected periodically +# +# In this elixir, we have a square domain that is divided into a left half and a right half. On each +# half of the domain, a completely independent SemidiscretizationHyperbolic is created for the +# linear advection equations. The two systems are coupled in the x-direction and have periodic +# boundaries in the y-direction. For a high-level overview, see also the figure below: +# +# (-1, 1) ( 1, 1) +# ┌────────────────────┬────────────────────┐ +# │ ↑ periodic ↑ │ ↑ periodic ↑ │ +# │ │ │ +# │ │ │ +# │ ========= │ ========= │ +# │ system #1 │ system #2 │ +# │ ========= │ ========= │ +# │ │ │ +# │ │ │ +# │ │ │ +# │ │ │ +# │ coupled -->│<-- coupled │ +# │ │ │ +# │<-- coupled │ coupled -->│ +# │ │ │ +# │ │ │ +# │ ↓ periodic ↓ │ ↓ periodic ↓ │ +# └────────────────────┴────────────────────┘ +# (-1, -1) ( 1, -1) + +advection_velocity = (0.2, -0.7) +equations = LinearScalarAdvectionEquation2D(advection_velocity) + +# Create DG solver with polynomial degree = 3 and (local) Lax-Friedrichs/Rusanov flux as surface flux +solver = DGSEM(polydeg=3, surface_flux=flux_lax_friedrichs) + +# First mesh is the left half of a [-1,1]^2 square +coordinates_min1 = (-1.0, -1.0) # minimum coordinates (min(x), min(y)) +coordinates_max1 = ( 0.0, 1.0) # maximum coordinates (max(x), max(y)) + +# Define identical resolution as a variable such that it is easier to change from `trixi_include` +cells_per_dimension = (8, 16) + +cells_per_dimension1 = cells_per_dimension + +mesh1 = StructuredMesh(cells_per_dimension1, coordinates_min1, coordinates_max1) + +# A semidiscretization collects data structures and functions for the spatial discretization +semi1 = SemidiscretizationHyperbolic(mesh1, equations, initial_condition_convergence_test, solver, + boundary_conditions=( + # Connect left boundary with right boundary of right mesh + x_neg=BoundaryConditionCoupled(2, (:end, :i_forward), Float64), + # Connect right boundary with left boundary of right mesh + x_pos=BoundaryConditionCoupled(2, (:begin, :i_forward), Float64), + y_neg=boundary_condition_periodic, + y_pos=boundary_condition_periodic)) + + +# Second mesh is the right half of a [-1,1]^2 square +coordinates_min2 = (0.0, -1.0) # minimum coordinates (min(x), min(y)) +coordinates_max2 = (1.0, 1.0) # maximum coordinates (max(x), max(y)) + +cells_per_dimension2 = cells_per_dimension + +mesh2 = StructuredMesh(cells_per_dimension2, coordinates_min2, coordinates_max2) + +semi2 = SemidiscretizationHyperbolic(mesh2, equations, initial_condition_convergence_test, solver, + boundary_conditions=( + # Connect left boundary with right boundary of left mesh + x_neg=BoundaryConditionCoupled(1, (:end, :i_forward), Float64), + # Connect right boundary with left boundary of left mesh + x_pos=BoundaryConditionCoupled(1, (:begin, :i_forward), Float64), + y_neg=boundary_condition_periodic, + y_pos=boundary_condition_periodic)) + +# Create a semidiscretization that bundles semi1 and semi2 +semi = SemidiscretizationCoupled(semi1, semi2) + +############################################################################### +# ODE solvers, callbacks etc. + +# Create ODE problem with time span from 0.0 to 2.0 +ode = semidiscretize(semi, (0.0, 2.0)); + +# At the beginning of the main loop, the SummaryCallback prints a summary of the simulation setup +# and resets the timers +summary_callback = SummaryCallback() + +# The AnalysisCallback allows to analyse the solution in regular intervals and prints the results +analysis_callback1 = AnalysisCallback(semi1, interval=100) +analysis_callback2 = AnalysisCallback(semi2, interval=100) +analysis_callback = AnalysisCallbackCoupled(semi, analysis_callback1, analysis_callback2) + +# The SaveSolutionCallback allows to save the solution to a file in regular intervals +save_solution = SaveSolutionCallback(interval=100, + solution_variables=cons2prim) + +# The StepsizeCallback handles the re-calculation of the maximum Δt after each time step +stepsize_callback = StepsizeCallback(cfl=1.6) + +# Create a CallbackSet to collect all callbacks such that they can be passed to the ODE solver +callbacks = CallbackSet(summary_callback, analysis_callback, save_solution, stepsize_callback) + + +############################################################################### +# run the simulation + +# OrdinaryDiffEq's `solve` method evolves the solution in time and executes the passed callbacks +sol = solve(ode, CarpenterKennedy2N54(williamson_condition=false), + dt=1.0, # solve needs some value here but it will be overwritten by the stepsize_callback + save_everystep=false, callback=callbacks); + +# Print the timer summary +summary_callback() diff --git a/src/Trixi.jl b/src/Trixi.jl index d5579aeea33..86e349c7dad 100644 --- a/src/Trixi.jl +++ b/src/Trixi.jl @@ -117,6 +117,7 @@ include("semidiscretization/semidiscretization.jl") include("semidiscretization/semidiscretization_hyperbolic.jl") include("semidiscretization/semidiscretization_hyperbolic_parabolic.jl") include("semidiscretization/semidiscretization_euler_acoustics.jl") +include("semidiscretization/semidiscretization_coupled.jl") include("callbacks_step/callbacks_step.jl") include("callbacks_stage/callbacks_stage.jl") include("semidiscretization/semidiscretization_euler_gravity.jl") @@ -184,7 +185,8 @@ export boundary_condition_do_nothing, boundary_condition_noslip_wall, boundary_condition_slip_wall, boundary_condition_wall, - BoundaryConditionNavierStokesWall, NoSlip, Adiabatic, Isothermal + BoundaryConditionNavierStokesWall, NoSlip, Adiabatic, Isothermal, + BoundaryConditionCoupled export initial_condition_convergence_test, source_terms_convergence_test export source_terms_harmonic @@ -229,12 +231,14 @@ export SemidiscretizationEulerAcoustics export SemidiscretizationEulerGravity, ParametersEulerGravity, timestep_gravity_erk52_3Sstar!, timestep_gravity_carpenter_kennedy_erk54_2N! +export SemidiscretizationCoupled + export SummaryCallback, SteadyStateCallback, AnalysisCallback, AliveCallback, SaveRestartCallback, SaveSolutionCallback, TimeSeriesCallback, VisualizationCallback, AveragingCallback, AMRCallback, StepsizeCallback, GlmSpeedCallback, LBMCollisionCallback, EulerAcousticsCouplingCallback, - TrivialCallback + TrivialCallback, AnalysisCallbackCoupled export load_mesh, load_time diff --git a/src/auxiliary/special_elixirs.jl b/src/auxiliary/special_elixirs.jl index da73b42e572..25bca8939ce 100644 --- a/src/auxiliary/special_elixirs.jl +++ b/src/auxiliary/special_elixirs.jl @@ -85,9 +85,21 @@ function convergence_test(mod::Module, elixir::AbstractString, iterations; kwarg println("#"^100) end - # number of variables - _, equations, _, _ = mesh_equations_solver_cache(mod.semi) + # Use raw error values to compute EOC + analyze_convergence(errors, iterations, mod.semi) +end + +# Analyze convergence for any semidiscretization +# Note: this intermediate method is to allow dispatching on the semidiscretization +function analyze_convergence(errors, iterations, semi::AbstractSemidiscretization) + _, equations, _, _ = mesh_equations_solver_cache(semi) variablenames = varnames(cons2cons, equations) + analyze_convergence(errors, iterations, variablenames) +end + +# This method is called with the collected error values to actually compute and print the EOC +function analyze_convergence(errors, iterations, + variablenames::Union{Tuple, AbstractArray}) nvariables = length(variablenames) # Reshape errors to get a matrix where the i-th row represents the i-th iteration diff --git a/src/callbacks_step/analysis.jl b/src/callbacks_step/analysis.jl index 2e038401df7..7fa2e21a244 100644 --- a/src/callbacks_step/analysis.jl +++ b/src/callbacks_step/analysis.jl @@ -84,11 +84,13 @@ function Base.show(io::IO, ::MIME"text/plain", end end +# This is the convenience constructor that gets called from the elixirs function AnalysisCallback(semi::AbstractSemidiscretization; kwargs...) mesh, equations, solver, cache = mesh_equations_solver_cache(semi) AnalysisCallback(mesh, equations, solver, cache; kwargs...) end +# This is the actual constructor function AnalysisCallback(mesh, equations::AbstractEquations, solver, cache; interval = 0, save_analysis = false, @@ -132,9 +134,18 @@ function AnalysisCallback(mesh, equations::AbstractEquations, solver, cache; initialize = initialize!) end +# This method gets called from OrdinaryDiffEq's `solve(...)` function initialize!(cb::DiscreteCallback{Condition, Affect!}, u_ode, t, integrator) where {Condition, Affect! <: AnalysisCallback} semi = integrator.p + du_ode = first(get_tmp_cache(integrator)) + initialize!(cb, u_ode, du_ode, t, integrator, semi) +end + +# This is the actual initialization method +# Note: we have this indirection to allow initializing a callback from the AnalysisCallbackCoupled +function initialize!(cb::DiscreteCallback{Condition, Affect!}, u_ode, du_ode, t, + integrator, semi) where {Condition, Affect! <: AnalysisCallback} initial_state_integrals = integrate(u_ode, semi) _, equations, _, _ = mesh_equations_solver_cache(semi) @@ -202,13 +213,21 @@ function initialize!(cb::DiscreteCallback{Condition, Affect!}, u_ode, t, # Note: For details see the actual callback function below analysis_callback.start_gc_time = Base.gc_time_ns() - analysis_callback(integrator) + analysis_callback(u_ode, du_ode, integrator, semi) return nothing end -# TODO: Taal refactor, allow passing an IO object (which could be devnull to avoid cluttering the console) +# This method gets called from OrdinaryDiffEq's `solve(...)` function (analysis_callback::AnalysisCallback)(integrator) semi = integrator.p + du_ode = first(get_tmp_cache(integrator)) + u_ode = integrator.u + analysis_callback(u_ode, du_ode, integrator, semi) +end + +# This method gets called internally as the main entry point to the AnalysiCallback +# TODO: Taal refactor, allow passing an IO object (which could be devnull to avoid cluttering the console) +function (analysis_callback::AnalysisCallback)(u_ode, du_ode, integrator, semi) mesh, equations, solver, cache = mesh_equations_solver_cache(semi) @unpack dt, t = integrator iter = integrator.stats.naccept @@ -300,15 +319,14 @@ function (analysis_callback::AnalysisCallback)(integrator) end # Calculate current time derivative (needed for semidiscrete entropy time derivative, residual, etc.) - du_ode = first(get_tmp_cache(integrator)) # `integrator.f` is usually just a call to `rhs!` # However, we want to allow users to modify the ODE RHS outside of Trixi.jl # and allow us to pass a combined ODE RHS to OrdinaryDiffEq, e.g., for # hyperbolic-parabolic systems. - @notimeit timer() integrator.f(du_ode, integrator.u, semi, t) - u = wrap_array(integrator.u, mesh, equations, solver, cache) + @notimeit timer() integrator.f(du_ode, u_ode, semi, t) + u = wrap_array(u_ode, mesh, equations, solver, cache) du = wrap_array(du_ode, mesh, equations, solver, cache) - l2_error, linf_error = analysis_callback(io, du, u, integrator.u, t, semi) + l2_error, linf_error = analysis_callback(io, du, u, u_ode, t, semi) mpi_println("─"^100) mpi_println() diff --git a/src/callbacks_step/save_solution.jl b/src/callbacks_step/save_solution.jl index 55f17bbc1c7..1fe0d6b1e15 100644 --- a/src/callbacks_step/save_solution.jl +++ b/src/callbacks_step/save_solution.jl @@ -141,14 +141,7 @@ function initialize_save_cb!(solution_callback::SaveSolutionCallback, u, t, inte mpi_isroot() && mkpath(solution_callback.output_directory) semi = integrator.p - mesh, _, _, _ = mesh_equations_solver_cache(semi) - @trixi_timeit timer() "I/O" begin - if mesh.unsaved_changes - mesh.current_filename = save_mesh_file(mesh, - solution_callback.output_directory) - mesh.unsaved_changes = false - end - end + @trixi_timeit timer() "I/O" save_mesh(semi, solution_callback.output_directory) if solution_callback.save_initial_solution solution_callback(integrator) @@ -157,6 +150,16 @@ function initialize_save_cb!(solution_callback::SaveSolutionCallback, u, t, inte return nothing end +# Save mesh for a general semidiscretization (default) +function save_mesh(semi::AbstractSemidiscretization, output_directory, timestep = 0) + mesh, _, _, _ = mesh_equations_solver_cache(semi) + + if mesh.unsaved_changes + mesh.current_filename = save_mesh_file(mesh, output_directory) + mesh.unsaved_changes = false + end +end + # this method is called to determine whether the callback should be activated function (solution_callback::SaveSolutionCallback)(u, t, integrator) @unpack interval_or_dt, save_final_solution = solution_callback @@ -174,41 +177,15 @@ end # this method is called when the callback is activated function (solution_callback::SaveSolutionCallback)(integrator) u_ode = integrator.u - @unpack t, dt = integrator - iter = integrator.stats.naccept semi = integrator.p - mesh, _, _, _ = mesh_equations_solver_cache(semi) + iter = integrator.stats.naccept @trixi_timeit timer() "I/O" begin - @trixi_timeit timer() "save mesh" if mesh.unsaved_changes - mesh.current_filename = save_mesh_file(mesh, - solution_callback.output_directory, - iter) - mesh.unsaved_changes = false - end - - element_variables = Dict{Symbol, Any}() - @trixi_timeit timer() "get element variables" begin - get_element_variables!(element_variables, u_ode, semi) - callbacks = integrator.opts.callback - if callbacks isa CallbackSet - for cb in callbacks.continuous_callbacks - get_element_variables!(element_variables, u_ode, semi, cb; - t = integrator.t, - iter = integrator.stats.naccept) - end - for cb in callbacks.discrete_callbacks - get_element_variables!(element_variables, u_ode, semi, cb; - t = integrator.t, - iter = integrator.stats.naccept) - end - end - end - - @trixi_timeit timer() "save solution" save_solution_file(u_ode, t, dt, iter, - semi, - solution_callback, - element_variables) + # Call high-level functions that dispatch on semidiscretization type + @trixi_timeit timer() "save mesh" save_mesh(semi, + solution_callback.output_directory, + iter) + save_solution_file(semi, u_ode, solution_callback, integrator) end # avoid re-evaluating possible FSAL stages @@ -216,13 +193,43 @@ function (solution_callback::SaveSolutionCallback)(integrator) return nothing end +@inline function save_solution_file(semi::AbstractSemidiscretization, u_ode, + solution_callback, + integrator; system = "") + @unpack t, dt = integrator + iter = integrator.stats.naccept + + element_variables = Dict{Symbol, Any}() + @trixi_timeit timer() "get element variables" begin + get_element_variables!(element_variables, u_ode, semi) + callbacks = integrator.opts.callback + if callbacks isa CallbackSet + for cb in callbacks.continuous_callbacks + get_element_variables!(element_variables, u_ode, semi, cb; + t = integrator.t, iter = iter) + end + for cb in callbacks.discrete_callbacks + get_element_variables!(element_variables, u_ode, semi, cb; + t = integrator.t, iter = iter) + end + end + end + + @trixi_timeit timer() "save solution" save_solution_file(u_ode, t, dt, iter, semi, + solution_callback, + element_variables, + system = system) +end + @inline function save_solution_file(u_ode, t, dt, iter, semi::AbstractSemidiscretization, solution_callback, - element_variables = Dict{Symbol, Any}()) + element_variables = Dict{Symbol, Any}(); + system = "") mesh, equations, solver, cache = mesh_equations_solver_cache(semi) u = wrap_array_native(u_ode, mesh, equations, solver, cache) save_solution_file(u, t, dt, iter, mesh, equations, solver, cache, - solution_callback, element_variables) + solution_callback, + element_variables; system = system) end # TODO: Taal refactor, move save_mesh_file? diff --git a/src/callbacks_step/stepsize.jl b/src/callbacks_step/stepsize.jl index 9e9f2d4885b..8b5cb958318 100644 --- a/src/callbacks_step/stepsize.jl +++ b/src/callbacks_step/stepsize.jl @@ -64,14 +64,11 @@ end t = integrator.t u_ode = integrator.u semi = integrator.p - mesh, equations, solver, cache = mesh_equations_solver_cache(semi) @unpack cfl_number = stepsize_callback - u = wrap_array(u_ode, mesh, equations, solver, cache) - dt = @trixi_timeit timer() "calculate dt" begin - cfl_number * max_dt(u, t, mesh, have_constant_speed(equations), equations, - solver, cache) - end + # Dispatch based on semidiscretization + dt = @trixi_timeit timer() "calculate dt" calculate_dt(u_ode, t, cfl_number, + semi) set_proposed_dt!(integrator, dt) integrator.opts.dtmax = dt @@ -83,6 +80,16 @@ end return nothing end +# General case for a single semidiscretization +function calculate_dt(u_ode, t, cfl_number, semi::AbstractSemidiscretization) + mesh, equations, solver, cache = mesh_equations_solver_cache(semi) + u = wrap_array(u_ode, mesh, equations, solver, cache) + + dt = cfl_number * max_dt(u, t, mesh, + have_constant_speed(equations), equations, + solver, cache) +end + # Time integration methods from the DiffEq ecosystem without adaptive time stepping on their own # such as `CarpenterKennedy2N54` require passing `dt=...` in `solve(ode, ...)`. Since we don't have # an integrator at this stage but only the ODE, this method will be used there. It's called in diff --git a/src/callbacks_step/summary.jl b/src/callbacks_step/summary.jl index a73b2a1913b..08e13d0b98d 100644 --- a/src/callbacks_step/summary.jl +++ b/src/callbacks_step/summary.jl @@ -152,15 +152,7 @@ function initialize_summary_callback(cb::DiscreteCallback, u, t, integrator) :indentation_level => 0) semi = integrator.p - show(io_context, MIME"text/plain"(), semi) - println(io, "\n") - mesh, equations, solver, _ = mesh_equations_solver_cache(semi) - show(io_context, MIME"text/plain"(), mesh) - println(io, "\n") - show(io_context, MIME"text/plain"(), equations) - println(io, "\n") - show(io_context, MIME"text/plain"(), solver) - println(io, "\n") + print_summary_semidiscretization(io_context, semi) callbacks = integrator.opts.callback if callbacks isa CallbackSet @@ -208,6 +200,18 @@ function initialize_summary_callback(cb::DiscreteCallback, u, t, integrator) return nothing end +function print_summary_semidiscretization(io::IO, semi::AbstractSemidiscretization) + show(io, MIME"text/plain"(), semi) + println(io, "\n") + mesh, equations, solver, _ = mesh_equations_solver_cache(semi) + show(io, MIME"text/plain"(), mesh) + println(io, "\n") + show(io, MIME"text/plain"(), equations) + println(io, "\n") + show(io, MIME"text/plain"(), solver) + println(io, "\n") +end + function (cb::DiscreteCallback{Condition, Affect!})(io::IO = stdout) where {Condition, Affect! <: typeof(summary_callback) diff --git a/src/meshes/mesh_io.jl b/src/meshes/mesh_io.jl index b9c462fa15a..ede85d80106 100644 --- a/src/meshes/mesh_io.jl +++ b/src/meshes/mesh_io.jl @@ -95,11 +95,15 @@ end # of the mesh, like its size and the type of boundary mapping function. # Then, within Trixi2Vtk, the StructuredMesh and its node coordinates are reconstructured from # these attributes for plotting purposes -function save_mesh_file(mesh::StructuredMesh, output_directory) +function save_mesh_file(mesh::StructuredMesh, output_directory; system = "") # Create output directory (if it does not exist) mkpath(output_directory) - filename = joinpath(output_directory, "mesh.h5") + if isempty(system) + filename = joinpath(output_directory, "mesh.h5") + else + filename = joinpath(output_directory, @sprintf("mesh_%s.h5", system)) + end # Open file (clobber existing content) h5open(filename, "w") do file diff --git a/src/semidiscretization/semidiscretization.jl b/src/semidiscretization/semidiscretization.jl index 8fef66d261e..ac312c57c89 100644 --- a/src/semidiscretization/semidiscretization.jl +++ b/src/semidiscretization/semidiscretization.jl @@ -76,7 +76,8 @@ function semidiscretize(semi::AbstractSemidiscretization, tspan) # mpi_isparallel() && MPI.Barrier(mpi_comm()) # See https://github.com/trixi-framework/Trixi.jl/issues/328 iip = true # is-inplace, i.e., we modify a vector when calling rhs! - return ODEProblem{iip}(rhs!, u0_ode, tspan, semi) + specialize = SciMLBase.FullSpecialize # specialize on rhs! and parameters (semi) + return ODEProblem{iip, specialize}(rhs!, u0_ode, tspan, semi) end """ @@ -93,7 +94,8 @@ function semidiscretize(semi::AbstractSemidiscretization, tspan, # mpi_isparallel() && MPI.Barrier(mpi_comm()) # See https://github.com/trixi-framework/Trixi.jl/issues/328 iip = true # is-inplace, i.e., we modify a vector when calling rhs! - return ODEProblem{iip}(rhs!, u0_ode, tspan, semi) + specialize = SciMLBase.FullSpecialize # specialize on rhs! and parameters (semi) + return ODEProblem{iip, specialize}(rhs!, u0_ode, tspan, semi) end """ diff --git a/src/semidiscretization/semidiscretization_coupled.jl b/src/semidiscretization/semidiscretization_coupled.jl new file mode 100644 index 00000000000..b7adff78425 --- /dev/null +++ b/src/semidiscretization/semidiscretization_coupled.jl @@ -0,0 +1,610 @@ +""" + SemidiscretizationCoupled + +A struct used to bundle multiple semidiscretizations. +[`semidiscretize`](@ref) will return an `ODEProblem` that synchronizes time steps between the semidiscretizations. +Each call of `rhs!` will call `rhs!` for each semidiscretization individually. +The semidiscretizations can be coupled by gluing meshes together using [`BoundaryConditionCoupled`](@ref). + +!!! warning "Experimental code" + This is an experimental feature and can change any time. +""" +struct SemidiscretizationCoupled{S, Indices, EquationList} <: AbstractSemidiscretization + semis::S + u_indices::Indices # u_ode[u_indices[i]] is the part of u_ode corresponding to semis[i] + performance_counter::PerformanceCounter +end + +""" + SemidiscretizationCoupled(semis...) + +Create a coupled semidiscretization that consists of the semidiscretizations passed as arguments. +""" +function SemidiscretizationCoupled(semis...) + @assert all(semi -> ndims(semi) == ndims(semis[1]), semis) "All semidiscretizations must have the same dimension!" + + # Number of coefficients for each semidiscretization + n_coefficients = zeros(Int, length(semis)) + for i in 1:length(semis) + _, equations, _, _ = mesh_equations_solver_cache(semis[i]) + n_coefficients[i] = ndofs(semis[i]) * nvariables(equations) + end + + # Compute range of coefficients associated with each semidiscretization and allocate coupled BCs + u_indices = Vector{UnitRange{Int}}(undef, length(semis)) + for i in 1:length(semis) + offset = sum(n_coefficients[1:(i - 1)]) + 1 + u_indices[i] = range(offset, length = n_coefficients[i]) + + allocate_coupled_boundary_conditions(semis[i]) + end + + performance_counter = PerformanceCounter() + + SemidiscretizationCoupled{typeof(semis), typeof(u_indices), typeof(performance_counter) + }(semis, u_indices, performance_counter) +end + +function Base.show(io::IO, semi::SemidiscretizationCoupled) + @nospecialize semi # reduce precompilation time + + print(io, "SemidiscretizationCoupled($(semi.semis))") +end + +function Base.show(io::IO, ::MIME"text/plain", semi::SemidiscretizationCoupled) + @nospecialize semi # reduce precompilation time + + if get(io, :compact, false) + show(io, semi) + else + summary_header(io, "SemidiscretizationCoupled") + summary_line(io, "#spatial dimensions", ndims(semi.semis[1])) + summary_line(io, "#systems", nsystems(semi)) + for i in eachsystem(semi) + summary_line(io, "system", i) + mesh, equations, solver, _ = mesh_equations_solver_cache(semi.semis[i]) + summary_line(increment_indent(io), "mesh", mesh |> typeof |> nameof) + summary_line(increment_indent(io), "equations", equations |> typeof |> nameof) + summary_line(increment_indent(io), "initial condition", + semi.semis[i].initial_condition) + # no boundary conditions since that could be too much + summary_line(increment_indent(io), "source terms", semi.semis[i].source_terms) + summary_line(increment_indent(io), "solver", solver |> typeof |> nameof) + end + summary_line(io, "total #DOFs", ndofs(semi)) + summary_footer(io) + end +end + +function print_summary_semidiscretization(io::IO, semi::SemidiscretizationCoupled) + show(io, MIME"text/plain"(), semi) + println(io, "\n") + for i in eachsystem(semi) + mesh, equations, solver, _ = mesh_equations_solver_cache(semi.semis[i]) + summary_header(io, "System #$i") + + summary_line(io, "mesh", mesh |> typeof |> nameof) + show(increment_indent(io), MIME"text/plain"(), mesh) + + summary_line(io, "equations", equations |> typeof |> nameof) + show(increment_indent(io), MIME"text/plain"(), equations) + + summary_line(io, "solver", solver |> typeof |> nameof) + show(increment_indent(io), MIME"text/plain"(), solver) + + summary_footer(io) + println(io, "\n") + end +end + +@inline Base.ndims(semi::SemidiscretizationCoupled) = ndims(semi.semis[1]) + +@inline nsystems(semi::SemidiscretizationCoupled) = length(semi.semis) + +@inline eachsystem(semi::SemidiscretizationCoupled) = Base.OneTo(nsystems(semi)) + +@inline Base.real(semi::SemidiscretizationCoupled) = promote_type(real.(semi.semis)...) + +@inline Base.eltype(semi::SemidiscretizationCoupled) = promote_type(eltype.(semi.semis)...) + +@inline function ndofs(semi::SemidiscretizationCoupled) + sum(ndofs, semi.semis) +end + +@inline function nelements(semi::SemidiscretizationCoupled) + return sum(semi.semis) do semi_ + mesh, equations, solver, cache = mesh_equations_solver_cache(semi_) + + nelements(mesh, solver, cache) + end +end + +function compute_coefficients(t, semi::SemidiscretizationCoupled) + @unpack u_indices = semi + + u_ode = Vector{real(semi)}(undef, u_indices[end][end]) + + for i in eachsystem(semi) + # Call `compute_coefficients` in `src/semidiscretization/semidiscretization.jl` + u_ode[u_indices[i]] .= compute_coefficients(t, semi.semis[i]) + end + + return u_ode +end + +@inline function get_system_u_ode(u_ode, index, semi::SemidiscretizationCoupled) + @view u_ode[semi.u_indices[index]] +end + +function rhs!(du_ode, u_ode, semi::SemidiscretizationCoupled, t) + @unpack u_indices = semi + + time_start = time_ns() + + @trixi_timeit timer() "copy to coupled boundaries" begin + for semi_ in semi.semis + copy_to_coupled_boundary!(semi_.boundary_conditions, u_ode, semi) + end + end + + # Call rhs! for each semidiscretization + for i in eachsystem(semi) + u_loc = get_system_u_ode(u_ode, i, semi) + du_loc = get_system_u_ode(du_ode, i, semi) + + @trixi_timeit timer() "system #$i" rhs!(du_loc, u_loc, semi.semis[i], t) + end + + runtime = time_ns() - time_start + put!(semi.performance_counter, runtime) + + return nothing +end + +################################################################################ +### AnalysisCallback +################################################################################ + +""" + AnalysisCallbackCoupled(semi, callbacks...) + +Combine multiple analysis callbacks for coupled simulations with a +[`SemidiscretizationCoupled`](@ref). For each coupled system, an indididual +[`AnalysisCallback`](@ref) **must** be created and passed to the `AnalysisCallbackCoupled` **in +order**, i.e., in the same sequence as the indidvidual semidiscretizations are stored in the +`SemidiscretizationCoupled`. + +!!! warning "Experimental code" + This is an experimental feature and can change any time. +""" +struct AnalysisCallbackCoupled{CB} + callbacks::CB +end + +function Base.show(io::IO, ::MIME"text/plain", + cb_coupled::DiscreteCallback{<:Any, <:AnalysisCallbackCoupled}) + @nospecialize cb_coupled # reduce precompilation time + + if get(io, :compact, false) + show(io, cb_coupled) + else + analysis_callback_coupled = cb_coupled.affect! + + summary_header(io, "AnalysisCallbackCoupled") + for (i, cb) in enumerate(analysis_callback_coupled.callbacks) + summary_line(io, "Callback #$i", "") + show(increment_indent(io), MIME"text/plain"(), cb) + end + summary_footer(io) + end +end + +# Convenience constructor for the coupled callback that gets called directly from the elixirs +function AnalysisCallbackCoupled(semi_coupled, callbacks...) + if length(callbacks) != nsystems(semi_coupled) + error("an AnalysisCallbackCoupled requires one AnalysisCallback for each semidiscretization") + end + + analysis_callback_coupled = AnalysisCallbackCoupled{typeof(callbacks)}(callbacks) + + # This callback is triggered if any of its subsidiary callbacks' condition is triggered + condition = (u, t, integrator) -> any(callbacks) do callback + callback.condition(u, t, integrator) + end + + DiscreteCallback(condition, analysis_callback_coupled, + save_positions = (false, false), + initialize = initialize!) +end + +# This method gets called during initialization from OrdinaryDiffEq's `solve(...)` +function initialize!(cb_coupled::DiscreteCallback{Condition, Affect!}, u_ode_coupled, t, + integrator) where {Condition, Affect! <: AnalysisCallbackCoupled} + analysis_callback_coupled = cb_coupled.affect! + semi_coupled = integrator.p + du_ode_coupled = first(get_tmp_cache(integrator)) + + # Loop over coupled systems' callbacks and initialize them individually + for i in eachsystem(semi_coupled) + cb = analysis_callback_coupled.callbacks[i] + semi = semi_coupled.semis[i] + u_ode = get_system_u_ode(u_ode_coupled, i, semi_coupled) + du_ode = get_system_u_ode(du_ode_coupled, i, semi_coupled) + initialize!(cb, u_ode, du_ode, t, integrator, semi) + end +end + +# This method gets called from OrdinaryDiffEq's `solve(...)` +function (analysis_callback_coupled::AnalysisCallbackCoupled)(integrator) + semi_coupled = integrator.p + u_ode_coupled = integrator.u + du_ode_coupled = first(get_tmp_cache(integrator)) + + # Loop over coupled systems' callbacks and call them individually + for i in eachsystem(semi_coupled) + @unpack condition = analysis_callback_coupled.callbacks[i] + analysis_callback = analysis_callback_coupled.callbacks[i].affect! + u_ode = get_system_u_ode(u_ode_coupled, i, semi_coupled) + + # Check condition and skip callback if it is not yet its turn + if !condition(u_ode, integrator.t, integrator) + continue + end + + semi = semi_coupled.semis[i] + du_ode = get_system_u_ode(du_ode_coupled, i, semi_coupled) + analysis_callback(u_ode, du_ode, integrator, semi) + end +end + +# used for error checks and EOC analysis +function (cb::DiscreteCallback{Condition, Affect!})(sol) where {Condition, + Affect! <: + AnalysisCallbackCoupled} + semi_coupled = sol.prob.p + u_ode_coupled = sol.u[end] + @unpack callbacks = cb.affect! + + uEltype = real(semi_coupled) + l2_error_collection = uEltype[] + linf_error_collection = uEltype[] + for i in eachsystem(semi_coupled) + analysis_callback = callbacks[i].affect! + @unpack analyzer = analysis_callback + cache_analysis = analysis_callback.cache + + semi = semi_coupled.semis[i] + u_ode = get_system_u_ode(u_ode_coupled, i, semi_coupled) + + l2_error, linf_error = calc_error_norms(u_ode, sol.t[end], analyzer, semi, + cache_analysis) + append!(l2_error_collection, l2_error) + append!(linf_error_collection, linf_error) + end + + (; l2 = l2_error_collection, linf = linf_error_collection) +end + +################################################################################ +### SaveSolutionCallback +################################################################################ + +# Save mesh for a coupled semidiscretization, which contains multiple meshes internally +function save_mesh(semi::SemidiscretizationCoupled, output_directory, timestep = 0) + for i in eachsystem(semi) + mesh, _, _, _ = mesh_equations_solver_cache(semi.semis[i]) + + if mesh.unsaved_changes + mesh.current_filename = save_mesh_file(mesh, output_directory, system = i) + mesh.unsaved_changes = false + end + end +end + +@inline function save_solution_file(semi::SemidiscretizationCoupled, u_ode, + solution_callback, + integrator) + @unpack semis = semi + + for i in eachsystem(semi) + u_ode_slice = get_system_u_ode(u_ode, i, semi) + save_solution_file(semis[i], u_ode_slice, solution_callback, integrator, system = i) + end +end + +################################################################################ +### StepsizeCallback +################################################################################ + +# In case of coupled system, use minimum timestep over all systems +function calculate_dt(u_ode, t, cfl_number, 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]) + end + + return dt +end + +################################################################################ +### Equations +################################################################################ + +""" + BoundaryConditionCoupled(other_semi_index, indices, uEltype) + +Boundary condition to glue two meshes together. Solution values at the boundary +of another mesh will be used as boundary values. This requires the use +of [`SemidiscretizationCoupled`](@ref). The other mesh is specified by `other_semi_index`, +which is the index of the mesh in the tuple of semidiscretizations. + +Note that the elements and nodes of the two meshes at the coupled boundary must coincide. +This is currently only implemented for [`StructuredMesh`](@ref). + +# Arguments +- `other_semi_index`: the index in `SemidiscretizationCoupled` of the semidiscretization + from which the values are copied +- `indices::Tuple`: node/cell indices at the boundary of the mesh in the other + semidiscretization. See examples below. +- `uEltype::Type`: element type of solution + +# Examples +```julia +# Connect the left boundary of mesh 2 to our boundary such that our positive +# boundary direction will match the positive y direction of the other boundary +BoundaryConditionCoupled(2, (:begin, :i), Float64) + +# Connect the same two boundaries oppositely oriented +BoundaryConditionCoupled(2, (:begin, :i_backwards), Float64) + +# Using this as y_neg boundary will connect `our_cells[i, 1, j]` to `other_cells[j, end-i, end]` +BoundaryConditionCoupled(2, (:j, :i_backwards, :end), Float64) +``` + +!!! warning "Experimental code" + This is an experimental feature and can change any time. +""" +mutable struct BoundaryConditionCoupled{NDIMS, NDIMST2M1, uEltype <: Real, Indices} + # NDIMST2M1 == NDIMS * 2 - 1 + # Buffer for boundary values: [variable, nodes_i, nodes_j, cell_i, cell_j] + u_boundary :: Array{uEltype, NDIMST2M1} # NDIMS * 2 - 1 + other_semi_index :: Int + other_orientation :: Int + indices :: Indices + + function BoundaryConditionCoupled(other_semi_index, indices, uEltype) + NDIMS = length(indices) + u_boundary = Array{uEltype, NDIMS * 2 - 1}(undef, ntuple(_ -> 0, NDIMS * 2 - 1)) + + if indices[1] in (:begin, :end) + other_orientation = 1 + elseif indices[2] in (:begin, :end) + other_orientation = 2 + else # indices[3] in (:begin, :end) + other_orientation = 3 + end + + new{NDIMS, NDIMS * 2 - 1, uEltype, typeof(indices)}(u_boundary, other_semi_index, + other_orientation, indices) + end +end + +function Base.eltype(boundary_condition::BoundaryConditionCoupled) + eltype(boundary_condition.u_boundary) +end + +function (boundary_condition::BoundaryConditionCoupled)(u_inner, orientation, direction, + cell_indices, surface_node_indices, + surface_flux_function, equations) + # get_node_vars(boundary_condition.u_boundary, equations, solver, surface_node_indices..., cell_indices...), + # but we don't have a solver here + u_boundary = SVector(ntuple(v -> boundary_condition.u_boundary[v, + surface_node_indices..., + cell_indices...], + 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) + end + + return flux +end + +function allocate_coupled_boundary_conditions(semi::AbstractSemidiscretization) + n_boundaries = 2 * ndims(semi) + mesh, equations, solver, _ = mesh_equations_solver_cache(semi) + + for direction in 1:n_boundaries + boundary_condition = semi.boundary_conditions[direction] + + allocate_coupled_boundary_condition(boundary_condition, direction, mesh, equations, + solver) + end +end + +# Don't do anything for other BCs than BoundaryConditionCoupled +function allocate_coupled_boundary_condition(boundary_condition, direction, mesh, equations, + solver) + return nothing +end + +# In 2D +function allocate_coupled_boundary_condition(boundary_condition::BoundaryConditionCoupled{2 + }, + direction, mesh, equations, dg::DGSEM) + if direction in (1, 2) + cell_size = size(mesh, 2) + else + cell_size = size(mesh, 1) + end + + uEltype = eltype(boundary_condition) + boundary_condition.u_boundary = Array{uEltype, 3}(undef, nvariables(equations), + nnodes(dg), + cell_size) +end + +# Don't do anything for other BCs than BoundaryConditionCoupled +function copy_to_coupled_boundary!(boundary_condition, u_ode, semi) + return nothing +end + +function copy_to_coupled_boundary!(boundary_conditions::Union{Tuple, NamedTuple}, u_ode, + semi) + for boundary_condition in boundary_conditions + copy_to_coupled_boundary!(boundary_condition, u_ode, semi) + end +end + +# In 2D +function copy_to_coupled_boundary!(boundary_condition::BoundaryConditionCoupled{2}, u_ode, + semi) + @unpack u_indices = semi + @unpack other_semi_index, other_orientation, indices = boundary_condition + + mesh, equations, solver, cache = mesh_equations_solver_cache(semi.semis[other_semi_index]) + u = wrap_array(get_system_u_ode(u_ode, other_semi_index, semi), mesh, equations, solver, + cache) + + linear_indices = LinearIndices(size(mesh)) + + if other_orientation == 1 + cells = axes(mesh, 2) + else # other_orientation == 2 + cells = axes(mesh, 1) + end + + # Copy solution data to the coupled boundary using "delayed indexing" with + # a start value and a step size to get the correct face and orientation. + node_index_range = eachnode(solver) + i_node_start, i_node_step = index_to_start_step_2d(indices[1], node_index_range) + j_node_start, j_node_step = index_to_start_step_2d(indices[2], node_index_range) + + i_cell_start, i_cell_step = index_to_start_step_2d(indices[1], axes(mesh, 1)) + j_cell_start, j_cell_step = index_to_start_step_2d(indices[2], axes(mesh, 2)) + + i_cell = i_cell_start + j_cell = j_cell_start + + for cell in cells + i_node = i_node_start + j_node = j_node_start + + for i in eachnode(solver) + for v in 1:size(u, 1) + boundary_condition.u_boundary[v, i, cell] = u[v, i_node, j_node, + linear_indices[i_cell, + j_cell]] + end + i_node += i_node_step + j_node += j_node_step + end + i_cell += i_cell_step + j_cell += j_cell_step + end +end + +################################################################################ +### DGSEM/structured +################################################################################ + +@inline function calc_boundary_flux_by_direction!(surface_flux_values, u, t, orientation, + boundary_condition::BoundaryConditionCoupled, + mesh::StructuredMesh, equations, + surface_integral, dg::DG, cache, + direction, node_indices, + surface_node_indices, element) + @unpack node_coordinates, contravariant_vectors, inverse_jacobian = cache.elements + @unpack surface_flux = surface_integral + + cell_indices = get_boundary_indices(element, orientation, mesh) + + u_inner = get_node_vars(u, equations, dg, node_indices..., element) + + # If the mapping is orientation-reversing, the contravariant vectors' orientation + # is reversed as well. The normal vector must be oriented in the direction + # from `left_element` to `right_element`, or the numerical flux will be computed + # incorrectly (downwind direction). + sign_jacobian = sign(inverse_jacobian[node_indices..., element]) + + # Contravariant vector Ja^i is the normal vector + normal = sign_jacobian * get_contravariant_vector(orientation, contravariant_vectors, + node_indices..., element) + + # If the mapping is orientation-reversing, the normal vector will be reversed (see above). + # However, the flux now has the wrong sign, since we need the physical flux in normal direction. + flux = sign_jacobian * boundary_condition(u_inner, normal, direction, cell_indices, + surface_node_indices, surface_flux, equations) + + for v in eachvariable(equations) + surface_flux_values[v, surface_node_indices..., direction, element] = flux[v] + end +end + +function get_boundary_indices(element, orientation, mesh::StructuredMesh{2}) + cartesian_indices = CartesianIndices(size(mesh)) + if orientation == 1 + # Get index of element in y-direction + cell_indices = (cartesian_indices[element][2],) + else # orientation == 2 + # Get index of element in x-direction + cell_indices = (cartesian_indices[element][1],) + end + + return cell_indices +end + +################################################################################ +### Special elixirs +################################################################################ + +# Analyze convergence for SemidiscretizationCoupled +function analyze_convergence(errors_coupled, iterations, + semi_coupled::SemidiscretizationCoupled) + # Extract errors: the errors are currently stored as + # | iter 1 sys 1 var 1...n | iter 1 sys 2 var 1...n | ... | iter 2 sys 1 var 1...n | ... + # but for calling `analyze_convergence` below, we need the following layout + # sys n: | iter 1 var 1...n | iter 1 var 1...n | ... | iter 2 var 1...n | ... + # That is, we need to extract and join the data for a single system + errors = Dict{Symbol, Vector{Float64}}[] + for i in eachsystem(semi_coupled) + push!(errors, Dict(:l2 => Float64[], :linf => Float64[])) + end + offset = 0 + for iter in 1:iterations, i in eachsystem(semi_coupled) + # Extract information on current semi + semi = semi_coupled.semis[i] + _, equations, _, _ = mesh_equations_solver_cache(semi) + variablenames = varnames(cons2cons, equations) + + # Compute offset + first = offset + 1 + last = offset + length(variablenames) + offset += length(variablenames) + + # Append errors to appropriate storage + append!(errors[i][:l2], errors_coupled[:l2][first:last]) + append!(errors[i][:linf], errors_coupled[:linf][first:last]) + end + + eoc_mean_values = Vector{Dict{Symbol, Any}}(undef, nsystems(semi_coupled)) + for i in eachsystem(semi_coupled) + # Use visual cues to separate output from multiple systems + println() + println("="^100) + println("# System $i") + println("="^100) + + # Extract information on current semi + semi = semi_coupled.semis[i] + _, equations, _, _ = mesh_equations_solver_cache(semi) + variablenames = varnames(cons2cons, equations) + + eoc_mean_values[i] = analyze_convergence(errors[i], iterations, variablenames) + end + + return eoc_mean_values +end diff --git a/src/semidiscretization/semidiscretization_euler_gravity.jl b/src/semidiscretization/semidiscretization_euler_gravity.jl index 665f2be9bfa..8fe9de1d2b2 100644 --- a/src/semidiscretization/semidiscretization_euler_gravity.jl +++ b/src/semidiscretization/semidiscretization_euler_gravity.jl @@ -477,18 +477,29 @@ end @inline function save_solution_file(u_ode, t, dt, iter, semi::SemidiscretizationEulerGravity, solution_callback, - element_variables = Dict{Symbol, Any}()) + element_variables = Dict{Symbol, Any}(); + system = "") + # If this is called already as part of a multi-system setup (i.e., system is non-empty), + # we build a combined system name + if !isempty(system) + system_euler = system * "_euler" + system_gravity = system * "_gravity" + else + system_euler = "euler" + system_gravity = "gravity" + end + u_euler = wrap_array_native(u_ode, semi.semi_euler) filename_euler = save_solution_file(u_euler, t, dt, iter, mesh_equations_solver_cache(semi.semi_euler)..., solution_callback, element_variables, - system = "euler") + system = system_euler) u_gravity = wrap_array_native(semi.cache.u_ode, semi.semi_gravity) filename_gravity = save_solution_file(u_gravity, t, dt, iter, mesh_equations_solver_cache(semi.semi_gravity)..., solution_callback, element_variables, - system = "gravity") + system = system_gravity) return filename_euler, filename_gravity end diff --git a/test/test_special_elixirs.jl b/test/test_special_elixirs.jl index 742a3abc376..23017059eaa 100644 --- a/test/test_special_elixirs.jl +++ b/test/test_special_elixirs.jl @@ -30,6 +30,12 @@ coverage = occursin("--code-coverage", cmd) && !occursin("--code-coverage=none", @test isapprox(mean_convergence[:l2], [4.0], rtol=0.05) end + @timed_testset "structured_2d_dgsem coupled" begin + mean_convergence = convergence_test(@__MODULE__, joinpath(EXAMPLES_DIR, "structured_2d_dgsem", "elixir_advection_coupled.jl"), 3) + @test isapprox(mean_convergence[1][:l2], [4.0], rtol=0.05) + @test isapprox(mean_convergence[2][:l2], [4.0], rtol=0.05) + end + @timed_testset "p4est_2d_dgsem" begin # Run convergence test on unrefined mesh no_refine = @cfunction((p4est, which_tree, quadrant) -> Cint(0), Cint, (Ptr{Trixi.p4est_t}, Ptr{Trixi.p4est_topidx_t}, Ptr{Trixi.p4est_quadrant_t})) @@ -57,6 +63,7 @@ coverage = occursin("--code-coverage", cmd) && !occursin("--code-coverage=none", @test_nowarn_mod convergence_test(@__MODULE__, joinpath(EXAMPLES_DIR, "tree_2d_dgsem", "elixir_advection_basic.jl"), 2, tspan=(0.0, 0.01)) @test_nowarn_mod convergence_test(@__MODULE__, joinpath(EXAMPLES_DIR, "tree_2d_dgsem", "elixir_advection_extended.jl"), 2, initial_refinement_level=0, tspan=(0.0, 0.1)) @test_nowarn_mod convergence_test(@__MODULE__, joinpath(EXAMPLES_DIR, "structured_2d_dgsem", "elixir_advection_basic.jl"), 2, tspan=(0.0, 0.01)) + @test_nowarn_mod convergence_test(@__MODULE__, joinpath(EXAMPLES_DIR, "structured_2d_dgsem", "elixir_advection_coupled.jl"), 2, tspan=(0.0, 0.01)) @test_nowarn_mod convergence_test(@__MODULE__, joinpath(EXAMPLES_DIR, "structured_2d_dgsem", "elixir_advection_extended.jl"), 2, cells_per_dimension=(1, 1), tspan=(0.0, 0.1)) end end diff --git a/test/test_structured_2d.jl b/test/test_structured_2d.jl index feaf66c4a7f..16fc72f0a46 100644 --- a/test/test_structured_2d.jl +++ b/test/test_structured_2d.jl @@ -19,6 +19,19 @@ isdir(outdir) && rm(outdir, recursive=true) linf = [6.627000273229378e-5]) end + @trixi_testset "elixir_advection_coupled.jl" begin + @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_advection_coupled.jl"), + l2 = [7.816742843181738e-6, 7.816742843196112e-6], + linf = [6.314906965543265e-5, 6.314906965410039e-5], + coverage_override = (maxiters=10^5,)) + + @testset "analysis_callback(sol) for AnalysisCallbackCoupled" begin + errors = analysis_callback(sol) + @test errors.l2 ≈ [7.816742843181738e-6, 7.816742843196112e-6] rtol=1.0e-4 + @test errors.linf ≈ [6.314906965543265e-5, 6.314906965410039e-5] rtol=1.0e-4 + end + end + @trixi_testset "elixir_advection_extended.jl" begin @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_advection_extended.jl"), l2 = [4.220397559713772e-6], From deb027adefdd88fccf6cec5ce4ca5c76106a0439 Mon Sep 17 00:00:00 2001 From: "dependabot[bot]" <49699333+dependabot[bot]@users.noreply.github.com> Date: Fri, 16 Jun 2023 22:53:10 +0200 Subject: [PATCH 02/15] Bump crate-ci/typos from 1.14.12 to 1.15.0 (#1524) * Bump crate-ci/typos from 1.14.12 to 1.15.0 Bumps [crate-ci/typos](https://github.com/crate-ci/typos) from 1.14.12 to 1.15.0. - [Release notes](https://github.com/crate-ci/typos/releases) - [Changelog](https://github.com/crate-ci/typos/blob/master/CHANGELOG.md) - [Commits](https://github.com/crate-ci/typos/compare/v1.14.12...v1.15.0) --- updated-dependencies: - dependency-name: crate-ci/typos dependency-type: direct:production update-type: version-update:semver-minor ... Signed-off-by: dependabot[bot] * Fix typos --------- Signed-off-by: dependabot[bot] Co-authored-by: dependabot[bot] <49699333+dependabot[bot]@users.noreply.github.com> Co-authored-by: Michael Schlottke-Lakemper --- .github/workflows/SpellCheck.yml | 2 +- docs/src/visualization.md | 2 +- examples/p4est_2d_dgsem/elixir_euler_supersonic_cylinder.jl | 2 +- 3 files changed, 3 insertions(+), 3 deletions(-) diff --git a/.github/workflows/SpellCheck.yml b/.github/workflows/SpellCheck.yml index c4ab3a98557..bc324c689bc 100644 --- a/.github/workflows/SpellCheck.yml +++ b/.github/workflows/SpellCheck.yml @@ -10,4 +10,4 @@ jobs: - name: Checkout Actions Repository uses: actions/checkout@v3 - name: Check spelling - uses: crate-ci/typos@v1.14.12 + uses: crate-ci/typos@v1.15.0 diff --git a/docs/src/visualization.md b/docs/src/visualization.md index e29313cc080..8f72bb4b1c6 100644 --- a/docs/src/visualization.md +++ b/docs/src/visualization.md @@ -339,7 +339,7 @@ create a [`PlotData1D`](@ref) with the keyword argument `curve` set to your list Let's give an example of this with the basic advection equation from above by creating a plot along the circle marked in green: -![2d-plot-along-cirlce](https://user-images.githubusercontent.com/72009492/130951042-e1849447-8e55-4798-9361-c8badb9f3a49.png) +![2d-plot-along-circle](https://user-images.githubusercontent.com/72009492/130951042-e1849447-8e55-4798-9361-c8badb9f3a49.png) We can write a function like this, that outputs a list of points on a circle: ```julia diff --git a/examples/p4est_2d_dgsem/elixir_euler_supersonic_cylinder.jl b/examples/p4est_2d_dgsem/elixir_euler_supersonic_cylinder.jl index 42370e861ce..366be700f9f 100644 --- a/examples/p4est_2d_dgsem/elixir_euler_supersonic_cylinder.jl +++ b/examples/p4est_2d_dgsem/elixir_euler_supersonic_cylinder.jl @@ -3,7 +3,7 @@ # Boundary conditions are supersonic Mach 3 inflow at the left portion of the domain # and supersonic outflow at the right portion of the domain. The top and bottom of the # channel as well as the cylinder are treated as Euler slip wall boundaries. -# This flow results in strong shock refletions / interactions as well as Kelvin-Helmholtz +# This flow results in strong shock reflections / interactions as well as Kelvin-Helmholtz # instabilities at later times as two Mach stems form above and below the cylinder. # # For complete details on the problem setup see Section 5.7 of the paper: From 642da1af9f9cc390f1d3d2a47a5fd07628a632f0 Mon Sep 17 00:00:00 2001 From: Joshua Lampert <51029046+JoshuaLampert@users.noreply.github.com> Date: Mon, 19 Jun 2023 09:07:15 +0200 Subject: [PATCH 03/15] Use `Pointerwrapper`s in `P4estMesh` (#1434) * use PointerWrappers * fix typo * fix typo * fixes * fixes * fixes * unsafe_load_sc returns PointerWrapper * bug fixes * bug fixes * bug fixes * add comment about * rename unsafe_load_* to load_pointerwrapper_* * format * fix merge conflicts * fix bad format again * fix * add unsafe_wrap_sc again * fix * fix * Update src/auxiliary/p4est.jl Co-authored-by: Michael Schlottke-Lakemper * introduce generic type PointerOrWrapper{T} * format --------- Co-authored-by: Michael Schlottke-Lakemper Co-authored-by: Hendrik Ranocha --- src/auxiliary/p4est.jl | 94 +++++++----- src/callbacks_step/amr.jl | 30 ++-- src/callbacks_step/amr_dg.jl | 3 +- src/callbacks_step/analysis.jl | 2 +- src/meshes/p4est_mesh.jl | 117 +++++++------- src/solvers/dgsem_p4est/containers.jl | 122 +++++++-------- src/solvers/dgsem_p4est/containers_2d.jl | 2 +- src/solvers/dgsem_p4est/containers_3d.jl | 2 +- .../dgsem_p4est/containers_parallel.jl | 145 +++++++++--------- src/solvers/dgsem_p4est/dg_parallel.jl | 80 +++++----- 10 files changed, 311 insertions(+), 286 deletions(-) diff --git a/src/auxiliary/p4est.jl b/src/auxiliary/p4est.jl index 93b5166cd81..968af339cbd 100644 --- a/src/auxiliary/p4est.jl +++ b/src/auxiliary/p4est.jl @@ -24,35 +24,38 @@ function init_p4est() return nothing end +# for convenience to either pass a Ptr or a PointerWrapper +const PointerOrWrapper = Union{Ptr{T}, PointerWrapper{T}} where {T} + # Convert sc_array of type T to Julia array -function unsafe_wrap_sc(::Type{T}, sc_array::Ptr{sc_array}) where {T} - sc_array_obj = unsafe_load(sc_array) +function unsafe_wrap_sc(::Type{T}, sc_array_ptr::Ptr{sc_array}) where {T} + sc_array_obj = unsafe_load(sc_array_ptr) return unsafe_wrap_sc(T, sc_array_obj) end function unsafe_wrap_sc(::Type{T}, sc_array_obj::sc_array) where {T} elem_count = sc_array_obj.elem_count array = sc_array_obj.array - return unsafe_wrap(Array, Ptr{T}(array), elem_count) end -# Load the ith element (1-indexed) of an sc array of type T -function unsafe_load_sc(::Type{T}, sc_array::Ptr{sc_array}, i = 1) where {T} - sc_array_obj = unsafe_load(sc_array) - return unsafe_load_sc(T, sc_array_obj, i) -end +function unsafe_wrap_sc(::Type{T}, sc_array_pw::PointerWrapper{sc_array}) where {T} + elem_count = sc_array_pw.elem_count[] + array = sc_array_pw.array -function unsafe_load_sc(::Type{T}, sc_array_obj::sc_array, i = 1) where {T} - element_size = sc_array_obj.elem_size - @assert element_size == sizeof(T) + return unsafe_wrap(Array, Ptr{T}(pointer(array)), elem_count) +end - return unsafe_load(Ptr{T}(sc_array_obj.array), i) +# Load the ith element (1-indexed) of an sc array of type T as PointerWrapper +function load_pointerwrapper_sc(::Type{T}, sc_array::PointerWrapper{sc_array}, + i::Integer = 1) where {T} + return PointerWrapper(T, pointer(sc_array.array) + (i - 1) * sizeof(T)) end # Create new `p4est` from a p4est_connectivity # 2D -function new_p4est(connectivity::Ptr{p4est_connectivity_t}, initial_refinement_level) +function new_p4est(connectivity::PointerOrWrapper{p4est_connectivity_t}, + initial_refinement_level) comm = P4est.uses_mpi() ? mpi_comm() : 0 # Use Trixi.jl's MPI communicator if p4est supports MPI p4est_new_ext(comm, connectivity, @@ -65,7 +68,8 @@ function new_p4est(connectivity::Ptr{p4est_connectivity_t}, initial_refinement_l end # 3D -function new_p4est(connectivity::Ptr{p8est_connectivity_t}, initial_refinement_level) +function new_p4est(connectivity::PointerOrWrapper{p8est_connectivity_t}, + initial_refinement_level) comm = P4est.uses_mpi() ? mpi_comm() : 0 # Use Trixi.jl's MPI communicator if p4est supports MPI p8est_new_ext(comm, connectivity, 0, initial_refinement_level, true, 2 * sizeof(Int), C_NULL, C_NULL) @@ -73,13 +77,13 @@ end # Save `p4est` data to file # 2D -function save_p4est!(file, p4est::Ptr{p4est_t}) +function save_p4est!(file, p4est::PointerOrWrapper{p4est_t}) # Don't save user data of the quads p4est_save(file, p4est, false) end # 3D -function save_p4est!(file, p8est::Ptr{p8est_t}) +function save_p4est!(file, p8est::PointerOrWrapper{p8est_t}) # Don't save user data of the quads p8est_save(file, p8est, false) end @@ -107,27 +111,33 @@ read_inp_p4est(meshfile, ::Val{3}) = p8est_connectivity_read_inp(meshfile) # Refine `p4est` if refine_fn_c returns 1 # 2D -function refine_p4est!(p4est::Ptr{p4est_t}, recursive, refine_fn_c, init_fn_c) +function refine_p4est!(p4est::PointerOrWrapper{p4est_t}, recursive, refine_fn_c, + init_fn_c) p4est_refine(p4est, recursive, refine_fn_c, init_fn_c) end # 3D -function refine_p4est!(p8est::Ptr{p8est_t}, recursive, refine_fn_c, init_fn_c) +function refine_p4est!(p8est::PointerOrWrapper{p8est_t}, recursive, refine_fn_c, + init_fn_c) p8est_refine(p8est, recursive, refine_fn_c, init_fn_c) end # Refine `p4est` if coarsen_fn_c returns 1 # 2D -function coarsen_p4est!(p4est::Ptr{p4est_t}, recursive, coarsen_fn_c, init_fn_c) +function coarsen_p4est!(p4est::PointerOrWrapper{p4est_t}, recursive, coarsen_fn_c, + init_fn_c) p4est_coarsen(p4est, recursive, coarsen_fn_c, init_fn_c) end # 3D -function coarsen_p4est!(p8est::Ptr{p8est_t}, recursive, coarsen_fn_c, init_fn_c) +function coarsen_p4est!(p8est::PointerOrWrapper{p8est_t}, recursive, coarsen_fn_c, + init_fn_c) p8est_coarsen(p8est, recursive, coarsen_fn_c, init_fn_c) end # Create new ghost layer from p4est, only connections via faces are relevant # 2D -ghost_new_p4est(p4est::Ptr{p4est_t}) = p4est_ghost_new(p4est, P4est.P4EST_CONNECT_FACE) +function ghost_new_p4est(p4est::PointerOrWrapper{p4est_t}) + p4est_ghost_new(p4est, P4est.P4EST_CONNECT_FACE) +end # 3D # In 3D it is not sufficient to use `P8EST_CONNECT_FACE`. Consider the neighbor elements of a mortar # in 3D. We have to determine which MPI ranks are involved in this mortar. @@ -147,28 +157,37 @@ ghost_new_p4est(p4est::Ptr{p4est_t}) = p4est_ghost_new(p4est, P4est.P4EST_CONNEC # `P8EST_CONNECT_FACE`. But if it is not in the ghost layer, it will not be available in # `iterate_p4est` and thus we cannot determine its MPI rank # (see https://github.com/cburstedde/p4est/blob/439bc9aae849555256ddfe4b03d1f9fe8d18ff0e/src/p8est_iterate.h#L66-L72). -ghost_new_p4est(p8est::Ptr{p8est_t}) = p8est_ghost_new(p8est, P4est.P8EST_CONNECT_FULL) +function ghost_new_p4est(p8est::PointerOrWrapper{p8est_t}) + p8est_ghost_new(p8est, P4est.P8EST_CONNECT_FULL) +end # Check if ghost layer is valid # 2D -function ghost_is_valid_p4est(p4est::Ptr{p4est_t}, ghost_layer::Ptr{p4est_ghost_t}) +function ghost_is_valid_p4est(p4est::PointerOrWrapper{p4est_t}, + ghost_layer::Ptr{p4est_ghost_t}) return p4est_ghost_is_valid(p4est, ghost_layer) end # 3D -function ghost_is_valid_p4est(p4est::Ptr{p8est_t}, ghost_layer::Ptr{p8est_ghost_t}) +function ghost_is_valid_p4est(p4est::PointerOrWrapper{p8est_t}, + ghost_layer::Ptr{p8est_ghost_t}) return p8est_ghost_is_valid(p4est, ghost_layer) end # Destroy ghost layer # 2D -ghost_destroy_p4est(ghost_layer::Ptr{p4est_ghost_t}) = p4est_ghost_destroy(ghost_layer) +function ghost_destroy_p4est(ghost_layer::PointerOrWrapper{p4est_ghost_t}) + p4est_ghost_destroy(ghost_layer) +end # 3D -ghost_destroy_p4est(ghost_layer::Ptr{p8est_ghost_t}) = p8est_ghost_destroy(ghost_layer) +function ghost_destroy_p4est(ghost_layer::PointerOrWrapper{p8est_ghost_t}) + p8est_ghost_destroy(ghost_layer) +end # Let `p4est` iterate over each cell volume and cell face. # Call iter_volume_c for each cell and iter_face_c for each face. # 2D -function iterate_p4est(p4est::Ptr{p4est_t}, user_data; ghost_layer = C_NULL, +function iterate_p4est(p4est::PointerOrWrapper{p4est_t}, user_data; + ghost_layer = C_NULL, iter_volume_c = C_NULL, iter_face_c = C_NULL) if user_data === C_NULL user_data_ptr = user_data @@ -191,7 +210,8 @@ function iterate_p4est(p4est::Ptr{p4est_t}, user_data; ghost_layer = C_NULL, end # 3D -function iterate_p4est(p8est::Ptr{p8est_t}, user_data; ghost_layer = C_NULL, +function iterate_p4est(p8est::PointerOrWrapper{p8est_t}, user_data; + ghost_layer = C_NULL, iter_volume_c = C_NULL, iter_face_c = C_NULL) if user_data === C_NULL user_data_ptr = user_data @@ -216,23 +236,25 @@ end # Load i-th element of the sc_array info.sides of the type p[48]est_iter_face_side_t # 2D version -function unsafe_load_side(info::Ptr{p4est_iter_face_info_t}, i = 1) - return unsafe_load_sc(p4est_iter_face_side_t, unsafe_load(info).sides, i) +function load_pointerwrapper_side(info::PointerWrapper{p4est_iter_face_info_t}, + i::Integer = 1) + return load_pointerwrapper_sc(p4est_iter_face_side_t, info.sides, i) end # 3D version -function unsafe_load_side(info::Ptr{p8est_iter_face_info_t}, i = 1) - return unsafe_load_sc(p8est_iter_face_side_t, unsafe_load(info).sides, i) +function load_pointerwrapper_side(info::PointerWrapper{p8est_iter_face_info_t}, + i::Integer = 1) + return load_pointerwrapper_sc(p8est_iter_face_side_t, info.sides, i) end # Load i-th element of the sc_array p4est.trees of the type p[48]est_tree_t # 2D version -function unsafe_load_tree(p4est::Ptr{p4est_t}, i = 1) - return unsafe_load_sc(p4est_tree_t, unsafe_load(p4est).trees, i) +function load_pointerwrapper_tree(p4est::PointerWrapper{p4est_t}, i::Integer = 1) + return load_pointerwrapper_sc(p4est_tree_t, p4est.trees, i) end # 3D version -function unsafe_load_tree(p8est::Ptr{p8est_t}, i = 1) - return unsafe_load_sc(p8est_tree_t, unsafe_load(p8est).trees, i) +function load_pointerwrapper_tree(p8est::PointerWrapper{p8est_t}, i::Integer = 1) + return load_pointerwrapper_sc(p8est_tree_t, p8est.trees, i) end end # @muladd diff --git a/src/callbacks_step/amr.jl b/src/callbacks_step/amr.jl index d6e19b79886..bef49b4c482 100644 --- a/src/callbacks_step/amr.jl +++ b/src/callbacks_step/amr.jl @@ -348,24 +348,24 @@ end # Copy controller values to quad user data storage, will be called below function copy_to_quad_iter_volume(info, user_data) - info_obj = unsafe_load(info) + info_pw = PointerWrapper(info) # Load tree from global trees array, one-based indexing - tree = unsafe_load_tree(info_obj.p4est, info_obj.treeid + 1) + tree_pw = load_pointerwrapper_tree(info_pw.p4est, info_pw.treeid[] + 1) # Quadrant numbering offset of this quadrant - offset = tree.quadrants_offset + offset = tree_pw.quadrants_offset[] # Global quad ID - quad_id = offset + info_obj.quadid + quad_id = offset + info_pw.quadid[] # Access user_data = lambda - user_data_ptr = Ptr{Int}(user_data) + user_data_pw = PointerWrapper(Int, user_data) # Load controller_value = lambda[quad_id + 1] - controller_value = unsafe_load(user_data_ptr, quad_id + 1) + controller_value = user_data_pw[quad_id + 1] # Access quadrant's user data ([global quad ID, controller_value]) - quad_data_ptr = Ptr{Int}(unsafe_load(info_obj.quad.p.user_data)) + quad_data_pw = PointerWrapper(Int, info_pw.quad.p.user_data[]) # Save controller value to quadrant's user data. - unsafe_store!(quad_data_ptr, controller_value, 2) + quad_data_pw[2] = controller_value return nothing end @@ -599,22 +599,22 @@ function current_element_levels(mesh::TreeMesh, solver, cache) end function extract_levels_iter_volume(info, user_data) - info_obj = unsafe_load(info) + info_pw = PointerWrapper(info) # Load tree from global trees array, one-based indexing - tree = unsafe_load_tree(info_obj.p4est, info_obj.treeid + 1) + tree_pw = load_pointerwrapper_tree(info_pw.p4est, info_pw.treeid[] + 1) # Quadrant numbering offset of this quadrant - offset = tree.quadrants_offset + offset = tree_pw.quadrants_offset[] # Global quad ID - quad_id = offset + info_obj.quadid + quad_id = offset + info_pw.quadid[] # Julia element ID element_id = quad_id + 1 - current_level = unsafe_load(info_obj.quad.level) + current_level = info_pw.quad.level[] # Unpack user_data = current_levels and save current element level - ptr = Ptr{Int}(user_data) - unsafe_store!(ptr, current_level, element_id) + pw = PointerWrapper(Int, user_data) + pw[element_id] = current_level return nothing end diff --git a/src/callbacks_step/amr_dg.jl b/src/callbacks_step/amr_dg.jl index 19bbebd9254..1dcfdccdea8 100644 --- a/src/callbacks_step/amr_dg.jl +++ b/src/callbacks_step/amr_dg.jl @@ -9,8 +9,7 @@ function rebalance_solver!(u_ode::AbstractVector, mesh::ParallelP4estMesh, equations, dg::DGSEM, cache, old_global_first_quadrant) # mpi ranks are 0-based, this array uses 1-based indices - global_first_quadrant = unsafe_wrap(Array, - unsafe_load(mesh.p4est).global_first_quadrant, + global_first_quadrant = unsafe_wrap(Array, mesh.p4est.global_first_quadrant, mpi_nranks() + 1) if global_first_quadrant[mpi_rank() + 1] == old_global_first_quadrant[mpi_rank() + 1] && diff --git a/src/callbacks_step/analysis.jl b/src/callbacks_step/analysis.jl index 7fa2e21a244..8cf43a1d15e 100644 --- a/src/callbacks_step/analysis.jl +++ b/src/callbacks_step/analysis.jl @@ -508,7 +508,7 @@ function print_amr_information(callbacks, mesh::P4estMesh, solver, cache) elements_per_level = zeros(P4EST_MAXLEVEL + 1) - for tree in unsafe_wrap_sc(p4est_tree_t, unsafe_load(mesh.p4est).trees) + for tree in unsafe_wrap_sc(p4est_tree_t, mesh.p4est.trees) elements_per_level .+= tree.quadrants_per_level end diff --git a/src/meshes/p4est_mesh.jl b/src/meshes/p4est_mesh.jl index ddd6cf473e4..60db285e04f 100644 --- a/src/meshes/p4est_mesh.jl +++ b/src/meshes/p4est_mesh.jl @@ -13,9 +13,9 @@ to manage trees and mesh refinement. """ mutable struct P4estMesh{NDIMS, RealT <: Real, IsParallel, P, Ghost, NDIMSP2, NNODES} <: AbstractMesh{NDIMS} - p4est::P # Either Ptr{p4est_t} or Ptr{p8est_t} - is_parallel::IsParallel - ghost::Ghost # Either Ptr{p4est_ghost_t} or Ptr{p8est_ghost_t} + p4est :: P # Either PointerWrapper{p4est_t} or PointerWrapper{p8est_t} + is_parallel :: IsParallel + ghost :: Ghost # Either PointerWrapper{p4est_ghost_t} or PointerWrapper{p8est_ghost_t} # Coordinates at the nodes specified by the tensor product of `nodes` (NDIMS times). # This specifies the geometry interpolation for each tree. tree_node_coordinates::Array{RealT, NDIMSP2} # [dimension, i, j, k, tree] @@ -43,18 +43,21 @@ mutable struct P4estMesh{NDIMS, RealT <: Real, IsParallel, P, Ghost, NDIMSP2, NN is_parallel = False() end + p4est_pw = PointerWrapper(p4est) + ghost = ghost_new_p4est(p4est) + ghost_pw = PointerWrapper(ghost) mesh = new{NDIMS, eltype(tree_node_coordinates), typeof(is_parallel), - typeof(p4est), typeof(ghost), NDIMS + 2, length(nodes)}(p4est, - is_parallel, - ghost, - tree_node_coordinates, - nodes, - boundary_names, - current_filename, - unsaved_changes, - p4est_partition_allow_for_coarsening) + typeof(p4est_pw), typeof(ghost_pw), NDIMS + 2, length(nodes)}(p4est_pw, + is_parallel, + ghost_pw, + tree_node_coordinates, + nodes, + boundary_names, + current_filename, + unsaved_changes, + p4est_partition_allow_for_coarsening) # Destroy `p4est` structs when the mesh is garbage collected finalizer(destroy_mesh, mesh) @@ -70,14 +73,14 @@ const ParallelP4estMesh{NDIMS} = P4estMesh{NDIMS, <:Real, <:True} @inline mpi_parallel(mesh::ParallelP4estMesh) = True() function destroy_mesh(mesh::P4estMesh{2}) - connectivity = unsafe_load(mesh.p4est).connectivity + connectivity = mesh.p4est.connectivity p4est_ghost_destroy(mesh.ghost) p4est_destroy(mesh.p4est) p4est_connectivity_destroy(connectivity) end function destroy_mesh(mesh::P4estMesh{3}) - connectivity = unsafe_load(mesh.p4est).connectivity + connectivity = mesh.p4est.connectivity p8est_ghost_destroy(mesh.ghost) p8est_destroy(mesh.p4est) p8est_connectivity_destroy(connectivity) @@ -87,11 +90,10 @@ end @inline Base.real(::P4estMesh{NDIMS, RealT}) where {NDIMS, RealT} = RealT @inline function ntrees(mesh::P4estMesh) - trees = unsafe_load(mesh.p4est).trees - return unsafe_load(trees).elem_count + return mesh.p4est.trees.elem_count[] end # returns Int32 by default which causes a weird method error when creating the cache -@inline ncells(mesh::P4estMesh) = Int(unsafe_load(mesh.p4est).local_num_quadrants) +@inline ncells(mesh::P4estMesh) = Int(mesh.p4est.local_num_quadrants[]) function Base.show(io::IO, mesh::P4estMesh) print(io, "P4estMesh{", ndims(mesh), ", ", real(mesh), "}") @@ -387,14 +389,14 @@ function p4est_mesh_from_hohqmesh_abaqus(meshfile, initial_refinement_level, n_dimensions, RealT) # Create the mesh connectivity using `p4est` connectivity = read_inp_p4est(meshfile, Val(n_dimensions)) - connectivity_obj = unsafe_load(connectivity) + connectivity_pw = PointerWrapper(connectivity) # These need to be of the type Int for unsafe_wrap below to work - n_trees::Int = connectivity_obj.num_trees - n_vertices::Int = connectivity_obj.num_vertices + n_trees::Int = connectivity_pw.num_trees[] + n_vertices::Int = connectivity_pw.num_vertices[] # Extract a copy of the element vertices to compute the tree node coordinates - vertices = unsafe_wrap(Array, connectivity_obj.vertices, (3, n_vertices)) + vertices = unsafe_wrap(Array, connectivity_pw.vertices, (3, n_vertices)) # Readin all the information from the mesh file into a string array file_lines = readlines(open(meshfile)) @@ -445,14 +447,14 @@ function p4est_mesh_from_standard_abaqus(meshfile, mapping, polydeg, initial_refinement_level, n_dimensions, RealT) # Create the mesh connectivity using `p4est` connectivity = read_inp_p4est(meshfile, Val(n_dimensions)) - connectivity_obj = unsafe_load(connectivity) + connectivity_pw = PointerWrapper(connectivity) # These need to be of the type Int for unsafe_wrap below to work - n_trees::Int = connectivity_obj.num_trees - n_vertices::Int = connectivity_obj.num_vertices + n_trees::Int = connectivity_pw.num_trees[] + n_vertices::Int = connectivity_pw.num_vertices[] - vertices = unsafe_wrap(Array, connectivity_obj.vertices, (3, n_vertices)) - tree_to_vertex = unsafe_wrap(Array, connectivity_obj.tree_to_vertex, + vertices = unsafe_wrap(Array, connectivity_pw.vertices, (3, n_vertices)) + tree_to_vertex = unsafe_wrap(Array, connectivity_pw.tree_to_vertex, (2^n_dimensions, n_trees)) basis = LobattoLegendreBasis(RealT, polydeg) @@ -1511,17 +1513,18 @@ end function update_ghost_layer!(mesh::P4estMesh) ghost_destroy_p4est(mesh.ghost) - mesh.ghost = ghost_new_p4est(mesh.p4est) + mesh.ghost = PointerWrapper(ghost_new_p4est(mesh.p4est)) end function init_fn(p4est, which_tree, quadrant) # Unpack quadrant's user data ([global quad ID, controller_value]) - ptr = Ptr{Int}(unsafe_load(quadrant.p.user_data)) + # Use `unsafe_load` here since `quadrant.p.user_data isa Ptr{Ptr{Nothing}}` + # and we only need the first (only!) entry + pw = PointerWrapper(Int, unsafe_load(quadrant.p.user_data)) # Initialize quad ID as -1 and controller_value as 0 (don't refine or coarsen) - unsafe_store!(ptr, -1, 1) - unsafe_store!(ptr, 0, 2) - + pw[1] = -1 + pw[2] = 0 return nothing end @@ -1539,8 +1542,10 @@ end function refine_fn(p4est, which_tree, quadrant) # Controller value has been copied to the quadrant's user data storage before. # Unpack quadrant's user data ([global quad ID, controller_value]). - ptr = Ptr{Int}(unsafe_load(quadrant.p.user_data)) - controller_value = unsafe_load(ptr, 2) + # Use `unsafe_load` here since `quadrant.p.user_data isa Ptr{Ptr{Nothing}}` + # and we only need the first (only!) entry + pw = PointerWrapper(Int, unsafe_load(quadrant.p.user_data)) + controller_value = pw[2] if controller_value > 0 # return true (refine) @@ -1586,9 +1591,9 @@ function coarsen_fn(p4est, which_tree, quadrants_ptr) # Controller value has been copied to the quadrant's user data storage before. # Load controller value from quadrant's user data ([global quad ID, controller_value]). - function controller_value(i) - unsafe_load(Ptr{Int}(unsafe_load(quadrants[i].p.user_data)), 2) - end + # Use `unsafe_load` here since `quadrant.p.user_data isa Ptr{Ptr{Nothing}}` + # and we only need the first (only!) entry + controller_value(i) = PointerWrapper(Int, unsafe_load(quadrants[i].p.user_data))[2] # `p4est` calls this function for each 2^ndims quads that could be coarsened to a single one. # Only coarsen if all these 2^ndims quads have been marked for coarsening. @@ -1671,20 +1676,19 @@ end # Copy global quad ID to quad's user data storage, will be called below function save_original_id_iter_volume(info, user_data) - info_obj = unsafe_load(info) + info_pw = PointerWrapper(info) # Load tree from global trees array, one-based indexing - tree = unsafe_load_tree(info_obj.p4est, info_obj.treeid + 1) + tree_pw = load_pointerwrapper_tree(info_pw.p4est, info_pw.treeid[] + 1) # Quadrant numbering offset of this quadrant - offset = tree.quadrants_offset + offset = tree_pw.quadrants_offset[] # Global quad ID - quad_id = offset + info_obj.quadid + quad_id = offset + info_pw.quadid[] # Unpack quadrant's user data ([global quad ID, controller_value]) - ptr = Ptr{Int}(unsafe_load(info_obj.quad.p.user_data)) + pw = PointerWrapper(Int, info_pw.quad.p.user_data[]) # Save global quad ID - unsafe_store!(ptr, quad_id, 1) - + pw[1] = quad_id return nothing end @@ -1708,24 +1712,23 @@ end # Extract information about which cells have been changed function collect_changed_iter_volume(info, user_data) - info_obj = unsafe_load(info) + info_pw = PointerWrapper(info) # The original element ID has been saved to user_data before. # Load original quad ID from quad's user data ([global quad ID, controller_value]). - quad_data_ptr = Ptr{Int}(unsafe_load(info_obj.quad.p.user_data)) - original_id = unsafe_load(quad_data_ptr, 1) + quad_data_pw = PointerWrapper(Int, info_pw.quad.p.user_data[]) + original_id = quad_data_pw[1] # original_id of cells that have been newly created is -1 if original_id >= 0 # Unpack user_data = original_cells - user_data_ptr = Ptr{Int}(user_data) + user_data_pw = PointerWrapper(Int, user_data) # If quad has an original_id, it existed before refinement/coarsening, # and therefore wasn't changed. # Mark original_id as "not changed during refinement/coarsening" in original_cells - unsafe_store!(user_data_ptr, 0, original_id + 1) + user_data_pw[original_id + 1] = 0 end - return nothing end @@ -1756,29 +1759,27 @@ end # Extract newly created cells function collect_new_iter_volume(info, user_data) - info_obj = unsafe_load(info) + info_pw = PointerWrapper(info) # The original element ID has been saved to user_data before. # Unpack quadrant's user data ([global quad ID, controller_value]). - quad_data_ptr = Ptr{Int}(unsafe_load(info_obj.quad.p.user_data)) - original_id = unsafe_load(quad_data_ptr, 1) + original_id = PointerWrapper(Int, info_pw.quad.p.user_data[])[1] # original_id of cells that have been newly created is -1 if original_id < 0 # Load tree from global trees array, one-based indexing - tree = unsafe_load_tree(info_obj.p4est, info_obj.treeid + 1) + tree_pw = load_pointerwrapper_tree(info_pw.p4est, info_pw.treeid[] + 1) # Quadrant numbering offset of this quadrant - offset = tree.quadrants_offset + offset = tree_pw.quadrants_offset[] # Global quad ID - quad_id = offset + info_obj.quadid + quad_id = offset + info_pw.quadid[] # Unpack user_data = original_cells - user_data_ptr = Ptr{Int}(user_data) + user_data_pw = PointerWrapper(Int, user_data) # Mark cell as "newly created during refinement/coarsening/balancing" - unsafe_store!(user_data_ptr, 1, quad_id + 1) + user_data_pw[quad_id + 1] = 1 end - return nothing end diff --git a/src/solvers/dgsem_p4est/containers.jl b/src/solvers/dgsem_p4est/containers.jl index 9b87de777a6..2b9c6987d24 100644 --- a/src/solvers/dgsem_p4est/containers.jl +++ b/src/solvers/dgsem_p4est/containers.jl @@ -276,18 +276,18 @@ function init_boundaries!(boundaries, mesh::P4estMesh) end # Function barrier for type stability -function init_boundaries_iter_face_inner(info, boundaries, boundary_id, mesh) +function init_boundaries_iter_face_inner(info_pw, boundaries, boundary_id, mesh) # Extract boundary data - side = unsafe_load_side(info) + side_pw = load_pointerwrapper_side(info_pw) # Get local tree, one-based indexing - tree = unsafe_load_tree(mesh.p4est, side.treeid + 1) + tree_pw = load_pointerwrapper_tree(mesh.p4est, side_pw.treeid[] + 1) # Quadrant numbering offset of this quadrant - offset = tree.quadrants_offset + offset = tree_pw.quadrants_offset[] # Verify before accessing is.full, but this should never happen - @assert side.is_hanging == false + @assert side_pw.is_hanging[] == false - local_quad_id = side.is.full.quadid + local_quad_id = side_pw.is.full.quadid[] # Global ID of this quad quad_id = offset + local_quad_id @@ -296,13 +296,13 @@ function init_boundaries_iter_face_inner(info, boundaries, boundary_id, mesh) boundaries.neighbor_ids[boundary_id] = quad_id + 1 # Face at which the boundary lies - face = side.face + face = side_pw.face[] # Save boundaries.node_indices dimension specific in containers_[23]d.jl init_boundary_node_indices!(boundaries, face, boundary_id) # One-based indexing - boundaries.name[boundary_id] = mesh.boundary_names[face + 1, side.treeid + 1] + boundaries.name[boundary_id] = mesh.boundary_names[face + 1, side_pw.treeid[] + 1] return nothing end @@ -479,32 +479,33 @@ end # Function barrier for type stability function init_surfaces_iter_face_inner(info, user_data) @unpack interfaces, mortars, boundaries = user_data - elem_count = unsafe_load(info).sides.elem_count + info_pw = PointerWrapper(info) + elem_count = info_pw.sides.elem_count[] if elem_count == 2 # Two neighboring elements => Interface or mortar # Extract surface data - sides = (unsafe_load_side(info, 1), unsafe_load_side(info, 2)) + sides_pw = (load_pointerwrapper_side(info_pw, 1), + load_pointerwrapper_side(info_pw, 2)) - if sides[1].is_hanging == false && sides[2].is_hanging == false + if sides_pw[1].is_hanging[] == false && sides_pw[2].is_hanging[] == false # No hanging nodes => normal interface if interfaces !== nothing - init_interfaces_iter_face_inner(info, sides, user_data) + init_interfaces_iter_face_inner(info_pw, sides_pw, user_data) end else # Hanging nodes => mortar if mortars !== nothing - init_mortars_iter_face_inner(info, sides, user_data) + init_mortars_iter_face_inner(info_pw, sides_pw, user_data) end end elseif elem_count == 1 # One neighboring elements => boundary if boundaries !== nothing - init_boundaries_iter_face_inner(info, user_data) + init_boundaries_iter_face_inner(info_pw, user_data) end end - return nothing end @@ -519,18 +520,18 @@ function init_surfaces!(interfaces, mortars, boundaries, mesh::P4estMesh) end # Initialization of interfaces after the function barrier -function init_interfaces_iter_face_inner(info, sides, user_data) +function init_interfaces_iter_face_inner(info_pw, sides_pw, user_data) @unpack interfaces, interface_id, mesh = user_data user_data.interface_id += 1 # Get Tuple of local trees, one-based indexing - trees = (unsafe_load_tree(mesh.p4est, sides[1].treeid + 1), - unsafe_load_tree(mesh.p4est, sides[2].treeid + 1)) + trees_pw = (load_pointerwrapper_tree(mesh.p4est, sides_pw[1].treeid[] + 1), + load_pointerwrapper_tree(mesh.p4est, sides_pw[2].treeid[] + 1)) # Quadrant numbering offsets of the quadrants at this interface - offsets = SVector(trees[1].quadrants_offset, - trees[2].quadrants_offset) + offsets = SVector(trees_pw[1].quadrants_offset[], + trees_pw[2].quadrants_offset[]) - local_quad_ids = SVector(sides[1].is.full.quadid, sides[2].is.full.quadid) + local_quad_ids = SVector(sides_pw[1].is.full.quadid[], sides_pw[2].is.full.quadid[]) # Global IDs of the neighboring quads quad_ids = offsets + local_quad_ids @@ -540,31 +541,30 @@ function init_interfaces_iter_face_inner(info, sides, user_data) interfaces.neighbor_ids[2, interface_id] = quad_ids[2] + 1 # Face at which the interface lies - faces = (sides[1].face, sides[2].face) + faces = (sides_pw[1].face[], sides_pw[2].face[]) # Save interfaces.node_indices dimension specific in containers_[23]d.jl - init_interface_node_indices!(interfaces, faces, - unsafe_load(info).orientation, interface_id) + init_interface_node_indices!(interfaces, faces, info_pw.orientation[], interface_id) return nothing end # Initialization of boundaries after the function barrier -function init_boundaries_iter_face_inner(info, user_data) +function init_boundaries_iter_face_inner(info_pw, user_data) @unpack boundaries, boundary_id, mesh = user_data user_data.boundary_id += 1 # Extract boundary data - side = unsafe_load_side(info) + side_pw = load_pointerwrapper_side(info_pw) # Get local tree, one-based indexing - tree = unsafe_load_tree(mesh.p4est, side.treeid + 1) + tree_pw = load_pointerwrapper_tree(mesh.p4est, side_pw.treeid[] + 1) # Quadrant numbering offset of this quadrant - offset = tree.quadrants_offset + offset = tree_pw.quadrants_offset[] # Verify before accessing is.full, but this should never happen - @assert side.is_hanging == false + @assert side_pw.is_hanging[] == false - local_quad_id = side.is.full.quadid + local_quad_id = side_pw.is.full.quadid[] # Global ID of this quad quad_id = offset + local_quad_id @@ -573,52 +573,52 @@ function init_boundaries_iter_face_inner(info, user_data) boundaries.neighbor_ids[boundary_id] = quad_id + 1 # Face at which the boundary lies - face = side.face + face = side_pw.face[] # Save boundaries.node_indices dimension specific in containers_[23]d.jl init_boundary_node_indices!(boundaries, face, boundary_id) # One-based indexing - boundaries.name[boundary_id] = mesh.boundary_names[face + 1, side.treeid + 1] + boundaries.name[boundary_id] = mesh.boundary_names[face + 1, side_pw.treeid[] + 1] return nothing end # Initialization of mortars after the function barrier -function init_mortars_iter_face_inner(info, sides, user_data) +function init_mortars_iter_face_inner(info_pw, sides_pw, user_data) @unpack mortars, mortar_id, mesh = user_data user_data.mortar_id += 1 # Get Tuple of local trees, one-based indexing - trees = (unsafe_load_tree(mesh.p4est, sides[1].treeid + 1), - unsafe_load_tree(mesh.p4est, sides[2].treeid + 1)) + trees_pw = (load_pointerwrapper_tree(mesh.p4est, sides_pw[1].treeid[] + 1), + load_pointerwrapper_tree(mesh.p4est, sides_pw[2].treeid[] + 1)) # Quadrant numbering offsets of the quadrants at this interface - offsets = SVector(trees[1].quadrants_offset, - trees[2].quadrants_offset) + offsets = SVector(trees_pw[1].quadrants_offset[], + trees_pw[2].quadrants_offset[]) - if sides[1].is_hanging == true + if sides_pw[1].is_hanging[] == true # Left is small, right is large - faces = (sides[1].face, sides[2].face) + faces = (sides_pw[1].face[], sides_pw[2].face[]) - local_small_quad_ids = sides[1].is.hanging.quadid + local_small_quad_ids = sides_pw[1].is.hanging.quadid[] # Global IDs of the two small quads small_quad_ids = offsets[1] .+ local_small_quad_ids # Just be sure before accessing is.full - @assert sides[2].is_hanging == false - large_quad_id = offsets[2] + sides[2].is.full.quadid - else # sides[2].is_hanging == true + @assert sides_pw[2].is_hanging[] == false + large_quad_id = offsets[2] + sides_pw[2].is.full.quadid[] + else # sides_pw[2].is_hanging[] == true # Right is small, left is large. # init_mortar_node_indices! below expects side 1 to contain the small elements. - faces = (sides[2].face, sides[1].face) + faces = (sides_pw[2].face[], sides_pw[1].face[]) - local_small_quad_ids = sides[2].is.hanging.quadid + local_small_quad_ids = sides_pw[2].is.hanging.quadid[] # Global IDs of the two small quads small_quad_ids = offsets[2] .+ local_small_quad_ids # Just be sure before accessing is.full - @assert sides[1].is_hanging == false - large_quad_id = offsets[1] + sides[1].is.full.quadid + @assert sides_pw[1].is_hanging[] == false + large_quad_id = offsets[1] + sides_pw[1].is.full.quadid[] end # Write data to mortar container, 1 and 2 are the small elements @@ -627,7 +627,7 @@ function init_mortars_iter_face_inner(info, sides, user_data) # Last entry is the large element mortars.neighbor_ids[end, mortar_id] = large_quad_id + 1 - init_mortar_node_indices!(mortars, faces, unsafe_load(info).orientation, mortar_id) + init_mortar_node_indices!(mortars, faces, info_pw.orientation[], mortar_id) return nothing end @@ -638,34 +638,36 @@ end # - boundaries # and collect the numbers in `user_data` in this order. function count_surfaces_iter_face(info, user_data) - elem_count = unsafe_load(info).sides.elem_count + info_pw = PointerWrapper(info) + elem_count = info_pw.sides.elem_count[] if elem_count == 2 # Two neighboring elements => Interface or mortar # Extract surface data - sides = (unsafe_load_side(info, 1), unsafe_load_side(info, 2)) + sides_pw = (load_pointerwrapper_side(info_pw, 1), + load_pointerwrapper_side(info_pw, 2)) - if sides[1].is_hanging == false && sides[2].is_hanging == false + if sides_pw[1].is_hanging[] == false && sides_pw[2].is_hanging[] == false # No hanging nodes => normal interface # Unpack user_data = [interface_count] and increment interface_count - ptr = Ptr{Int}(user_data) - id = unsafe_load(ptr, 1) - unsafe_store!(ptr, id + 1, 1) + pw = PointerWrapper(Int, user_data) + id = pw[1] + pw[1] = id + 1 else # Hanging nodes => mortar # Unpack user_data = [mortar_count] and increment mortar_count - ptr = Ptr{Int}(user_data) - id = unsafe_load(ptr, 2) - unsafe_store!(ptr, id + 1, 2) + pw = PointerWrapper(Int, user_data) + id = pw[2] + pw[2] = id + 1 end elseif elem_count == 1 # One neighboring elements => boundary # Unpack user_data = [boundary_count] and increment boundary_count - ptr = Ptr{Int}(user_data) - id = unsafe_load(ptr, 3) - unsafe_store!(ptr, id + 1, 3) + pw = PointerWrapper(Int, user_data) + id = pw[3] + pw[3] = id + 1 end return nothing diff --git a/src/solvers/dgsem_p4est/containers_2d.jl b/src/solvers/dgsem_p4est/containers_2d.jl index 4f7d903897a..11747f1f175 100644 --- a/src/solvers/dgsem_p4est/containers_2d.jl +++ b/src/solvers/dgsem_p4est/containers_2d.jl @@ -52,7 +52,7 @@ function calc_node_coordinates!(node_coordinates, p4est_root_len = 1 << P4EST_MAXLEVEL p4est_quadrant_len(l) = 1 << (P4EST_MAXLEVEL - l) - trees = unsafe_wrap_sc(p4est_tree_t, unsafe_load(mesh.p4est).trees) + trees = unsafe_wrap_sc(p4est_tree_t, mesh.p4est.trees) for tree in eachindex(trees) offset = trees[tree].quadrants_offset diff --git a/src/solvers/dgsem_p4est/containers_3d.jl b/src/solvers/dgsem_p4est/containers_3d.jl index 6cdc2cf9611..e9994fe4569 100644 --- a/src/solvers/dgsem_p4est/containers_3d.jl +++ b/src/solvers/dgsem_p4est/containers_3d.jl @@ -43,7 +43,7 @@ function calc_node_coordinates!(node_coordinates, p4est_root_len = 1 << P4EST_MAXLEVEL p4est_quadrant_len(l) = 1 << (P4EST_MAXLEVEL - l) - trees = unsafe_wrap_sc(p8est_tree_t, unsafe_load(mesh.p4est).trees) + trees = unsafe_wrap_sc(p8est_tree_t, mesh.p4est.trees) for tree in eachindex(trees) offset = trees[tree].quadrants_offset diff --git a/src/solvers/dgsem_p4est/containers_parallel.jl b/src/solvers/dgsem_p4est/containers_parallel.jl index 42d6ea44c5e..e7ee1f81478 100644 --- a/src/solvers/dgsem_p4est/containers_parallel.jl +++ b/src/solvers/dgsem_p4est/containers_parallel.jl @@ -311,21 +311,24 @@ function init_surfaces_iter_face_inner(info, # surfaces at once or any subset of them, some of the unpacked values above may be `nothing` if # they're not supposed to be initialized during this call. That is why we need additional # `!== nothing` checks below before initializing individual faces. - if unsafe_load(info).sides.elem_count == 2 + info_pw = PointerWrapper(info) + if info_pw.sides.elem_count[] == 2 # Two neighboring elements => Interface or mortar # Extract surface data - sides = (unsafe_load_side(info, 1), unsafe_load_side(info, 2)) + sides_pw = (load_pointerwrapper_side(info_pw, 1), + load_pointerwrapper_side(info_pw, 2)) - if sides[1].is_hanging == false && sides[2].is_hanging == false + if sides_pw[1].is_hanging[] == false && sides_pw[2].is_hanging[] == false # No hanging nodes => normal interface or MPI interface - if sides[1].is.full.is_ghost == true || sides[2].is.full.is_ghost == true # remote side => MPI interface + if sides_pw[1].is.full.is_ghost[] == true || + sides_pw[2].is.full.is_ghost[] == true # remote side => MPI interface if mpi_interfaces !== nothing - init_mpi_interfaces_iter_face_inner(info, sides, user_data) + init_mpi_interfaces_iter_face_inner(info_pw, sides_pw, user_data) end else if interfaces !== nothing - init_interfaces_iter_face_inner(info, sides, user_data) + init_interfaces_iter_face_inner(info_pw, sides_pw, user_data) end end else @@ -333,18 +336,18 @@ function init_surfaces_iter_face_inner(info, # First, we check which side is hanging, i.e., on which side we have the refined cells. # Then we check if any of the refined cells or the coarse cell are "ghost" cells, i.e., they # belong to another rank. That way we can determine if this is a regular mortar or MPI mortar - if sides[1].is_hanging == true - @assert sides[2].is_hanging == false - if any(sides[1].is.hanging.is_ghost .== true) || - sides[2].is.full.is_ghost == true + if sides_pw[1].is_hanging[] == true + @assert sides_pw[2].is_hanging[] == false + if any(sides_pw[1].is.hanging.is_ghost[] .== true) || + sides_pw[2].is.full.is_ghost[] == true face_has_ghost_side = true else face_has_ghost_side = false end - else # sides[2].is_hanging == true - @assert sides[1].is_hanging == false - if sides[1].is.full.is_ghost == true || - any(sides[2].is.hanging.is_ghost .== true) + else # sides_pw[2].is_hanging[] == true + @assert sides_pw[1].is_hanging[] == false + if sides_pw[1].is.full.is_ghost[] == true || + any(sides_pw[2].is.hanging.is_ghost[] .== true) face_has_ghost_side = true else face_has_ghost_side = false @@ -352,15 +355,15 @@ function init_surfaces_iter_face_inner(info, end # Initialize mortar or MPI mortar if face_has_ghost_side && mpi_mortars !== nothing - init_mpi_mortars_iter_face_inner(info, sides, user_data) + init_mpi_mortars_iter_face_inner(info_pw, sides_pw, user_data) elseif !face_has_ghost_side && mortars !== nothing - init_mortars_iter_face_inner(info, sides, user_data) + init_mortars_iter_face_inner(info_pw, sides_pw, user_data) end end - elseif unsafe_load(info).sides.elem_count == 1 + elseif info_pw.sides.elem_count[] == 1 # One neighboring elements => boundary if boundaries !== nothing - init_boundaries_iter_face_inner(info, user_data) + init_boundaries_iter_face_inner(info_pw, user_data) end end @@ -381,23 +384,23 @@ function init_surfaces!(interfaces, mortars, boundaries, mpi_interfaces, mpi_mor end # Initialization of MPI interfaces after the function barrier -function init_mpi_interfaces_iter_face_inner(info, sides, user_data) +function init_mpi_interfaces_iter_face_inner(info_pw, sides_pw, user_data) @unpack mpi_interfaces, mpi_interface_id, mesh = user_data user_data.mpi_interface_id += 1 - if sides[1].is.full.is_ghost == true + if sides_pw[1].is.full.is_ghost[] == true local_side = 2 - elseif sides[2].is.full.is_ghost == true + elseif sides_pw[2].is.full.is_ghost[] == true local_side = 1 else error("should not happen") end # Get local tree, one-based indexing - tree = unsafe_load_tree(mesh.p4est, sides[local_side].treeid + 1) + tree_pw = load_pointerwrapper_tree(mesh.p4est, sides_pw[local_side].treeid[] + 1) # Quadrant numbering offset of the local quadrant at this interface - offset = tree.quadrants_offset - tree_quad_id = sides[local_side].is.full.quadid # quadid in the local tree + offset = tree_pw.quadrants_offset[] + tree_quad_id = sides_pw[local_side].is.full.quadid[] # quadid in the local tree # ID of the local neighboring quad, cumulative over local trees local_quad_id = offset + tree_quad_id @@ -406,52 +409,52 @@ function init_mpi_interfaces_iter_face_inner(info, sides, user_data) mpi_interfaces.local_sides[mpi_interface_id] = local_side # Face at which the interface lies - faces = (sides[1].face, sides[2].face) + faces = (sides_pw[1].face[], sides_pw[2].face[]) # Save mpi_interfaces.node_indices dimension specific in containers_[23]d_parallel.jl init_mpi_interface_node_indices!(mpi_interfaces, faces, local_side, - unsafe_load(info).orientation, + info_pw.orientation[], mpi_interface_id) return nothing end # Initialization of MPI mortars after the function barrier -function init_mpi_mortars_iter_face_inner(info, sides, user_data) +function init_mpi_mortars_iter_face_inner(info_pw, sides_pw, user_data) @unpack mpi_mortars, mpi_mortar_id, mesh = user_data user_data.mpi_mortar_id += 1 # Get Tuple of adjacent trees, one-based indexing - trees = (unsafe_load_tree(mesh.p4est, sides[1].treeid + 1), - unsafe_load_tree(mesh.p4est, sides[2].treeid + 1)) + trees_pw = (load_pointerwrapper_tree(mesh.p4est, sides_pw[1].treeid[] + 1), + load_pointerwrapper_tree(mesh.p4est, sides_pw[2].treeid[] + 1)) # Quadrant numbering offsets of the quadrants at this mortar - offsets = SVector(trees[1].quadrants_offset, - trees[2].quadrants_offset) + offsets = SVector(trees_pw[1].quadrants_offset[], + trees_pw[2].quadrants_offset[]) - if sides[1].is_hanging == true + if sides_pw[1].is_hanging[] == true hanging_side = 1 full_side = 2 - else # sides[2].is_hanging == true + else # sides_pw[2].is_hanging[] == true hanging_side = 2 full_side = 1 end # Just be sure before accessing is.full or is.hanging later - @assert sides[full_side].is_hanging == false - @assert sides[hanging_side].is_hanging == true + @assert sides_pw[full_side].is_hanging[] == false + @assert sides_pw[hanging_side].is_hanging[] == true # Find small quads that are locally available - local_small_quad_positions = findall(sides[hanging_side].is.hanging.is_ghost .== + local_small_quad_positions = findall(sides_pw[hanging_side].is.hanging.is_ghost[] .== false) # Get id of local small quadrants within their tree # Indexing CBinding.Caccessor via a Vector does not work here -> use map instead - tree_small_quad_ids = map(p -> sides[hanging_side].is.hanging.quadid[p], + tree_small_quad_ids = map(p -> sides_pw[hanging_side].is.hanging.quadid[][p], local_small_quad_positions) local_small_quad_ids = offsets[hanging_side] .+ tree_small_quad_ids # ids cumulative over local trees # Determine if large quadrant is available and if yes, determine its id - if sides[full_side].is.full.is_ghost == false - local_large_quad_id = offsets[full_side] + sides[full_side].is.full.quadid + if sides_pw[full_side].is.full.is_ghost[] == false + local_large_quad_id = offsets[full_side] + sides_pw[full_side].is.full.quadid[] else local_large_quad_id = -1 # large quad is ghost end @@ -470,9 +473,8 @@ function init_mpi_mortars_iter_face_inner(info, sides, user_data) mpi_mortars.local_neighbor_positions[mpi_mortar_id] = local_neighbor_positions # init_mortar_node_indices! expects side 1 to contain small elements - faces = (sides[hanging_side].face, sides[full_side].face) - init_mortar_node_indices!(mpi_mortars, faces, unsafe_load(info).orientation, - mpi_mortar_id) + faces = (sides_pw[hanging_side].face[], sides_pw[full_side].face[]) + init_mortar_node_indices!(mpi_mortars, faces, info_pw.orientation[], mpi_mortar_id) return nothing end @@ -485,42 +487,45 @@ end # - (MPI) mortars at subdomain boundaries # and collect the numbers in `user_data` in this order. function count_surfaces_iter_face_parallel(info, user_data) - if unsafe_load(info).sides.elem_count == 2 + info_pw = PointerWrapper(info) + if info_pw.sides.elem_count[] == 2 # Two neighboring elements => Interface or mortar # Extract surface data - sides = (unsafe_load_side(info, 1), unsafe_load_side(info, 2)) + sides_pw = (load_pointerwrapper_side(info_pw, 1), + load_pointerwrapper_side(info_pw, 2)) - if sides[1].is_hanging == false && sides[2].is_hanging == false + if sides_pw[1].is_hanging[] == false && sides_pw[2].is_hanging[] == false # No hanging nodes => normal interface or MPI interface - if sides[1].is.full.is_ghost == true || sides[2].is.full.is_ghost == true # remote side => MPI interface + if sides_pw[1].is.full.is_ghost[] == true || + sides_pw[2].is.full.is_ghost[] == true # remote side => MPI interface # Unpack user_data = [mpi_interface_count] and increment mpi_interface_count - ptr = Ptr{Int}(user_data) - id = unsafe_load(ptr, 4) - unsafe_store!(ptr, id + 1, 4) + pw = PointerWrapper(Int, user_data) + id = pw[4] + pw[4] = id + 1 else # Unpack user_data = [interface_count] and increment interface_count - ptr = Ptr{Int}(user_data) - id = unsafe_load(ptr, 1) - unsafe_store!(ptr, id + 1, 1) + pw = PointerWrapper(Int, user_data) + id = pw[1] + pw[1] = id + 1 end else # Hanging nodes => mortar or MPI mortar # First, we check which side is hanging, i.e., on which side we have the refined cells. # Then we check if any of the refined cells or the coarse cell are "ghost" cells, i.e., they # belong to another rank. That way we can determine if this is a regular mortar or MPI mortar - if sides[1].is_hanging == true - @assert sides[2].is_hanging == false - if any(sides[1].is.hanging.is_ghost .== true) || - sides[2].is.full.is_ghost == true + if sides_pw[1].is_hanging[] == true + @assert sides_pw[2].is_hanging[] == false + if any(sides_pw[1].is.hanging.is_ghost[] .== true) || + sides_pw[2].is.full.is_ghost[] == true face_has_ghost_side = true else face_has_ghost_side = false end - else # sides[2].is_hanging == true - @assert sides[1].is_hanging == false - if sides[1].is.full.is_ghost == true || - any(sides[2].is.hanging.is_ghost .== true) + else # sides_pw[2].is_hanging[] == true + @assert sides_pw[1].is_hanging[] == false + if sides_pw[1].is.full.is_ghost[] == true || + any(sides_pw[2].is.hanging.is_ghost[] .== true) face_has_ghost_side = true else face_has_ghost_side = false @@ -528,23 +533,23 @@ function count_surfaces_iter_face_parallel(info, user_data) end if face_has_ghost_side # Unpack user_data = [mpi_mortar_count] and increment mpi_mortar_count - ptr = Ptr{Int}(user_data) - id = unsafe_load(ptr, 5) - unsafe_store!(ptr, id + 1, 5) + pw = PointerWrapper(Int, user_data) + id = pw[5] + pw[5] = id + 1 else # Unpack user_data = [mortar_count] and increment mortar_count - ptr = Ptr{Int}(user_data) - id = unsafe_load(ptr, 2) - unsafe_store!(ptr, id + 1, 2) + pw = PointerWrapper(Int, user_data) + id = pw[2] + pw[2] = id + 1 end end - elseif unsafe_load(info).sides.elem_count == 1 + elseif info_pw.sides.elem_count[] == 1 # One neighboring elements => boundary # Unpack user_data = [boundary_count] and increment boundary_count - ptr = Ptr{Int}(user_data) - id = unsafe_load(ptr, 3) - unsafe_store!(ptr, id + 1, 3) + pw = PointerWrapper(Int, user_data) + id = pw[3] + pw[3] = id + 1 end return nothing diff --git a/src/solvers/dgsem_p4est/dg_parallel.jl b/src/solvers/dgsem_p4est/dg_parallel.jl index ac122d048c1..324bc7f3cd6 100644 --- a/src/solvers/dgsem_p4est/dg_parallel.jl +++ b/src/solvers/dgsem_p4est/dg_parallel.jl @@ -263,15 +263,13 @@ function init_mpi_cache!(mpi_cache::P4estMPICache, mesh::ParallelP4estMesh, uEltype) # Determine local and total number of elements - n_elements_global = Int(unsafe_load(mesh.p4est).global_num_quadrants) - n_elements_by_rank = vcat(Int.(unsafe_wrap(Array, - unsafe_load(mesh.p4est).global_first_quadrant, + n_elements_global = Int(mesh.p4est.global_num_quadrants[]) + n_elements_by_rank = vcat(Int.(unsafe_wrap(Array, mesh.p4est.global_first_quadrant, mpi_nranks())), n_elements_global) |> diff # diff sufficient due to 0-based quad indices n_elements_by_rank = OffsetArray(n_elements_by_rank, 0:(mpi_nranks() - 1)) # Account for 1-based indexing in Julia - first_element_global_id = Int(unsafe_load(unsafe_load(mesh.p4est).global_first_quadrant, - mpi_rank() + 1)) + 1 + first_element_global_id = Int(mesh.p4est.global_first_quadrant[mpi_rank() + 1]) + 1 @assert n_elements_global==sum(n_elements_by_rank) "error in total number of elements" # TODO reuse existing structures @@ -379,17 +377,19 @@ function init_neighbor_rank_connectivity_iter_face_inner(info, user_data) @unpack interfaces, interface_id, global_interface_ids, neighbor_ranks_interface, mortars, mortar_id, global_mortar_ids, neighbor_ranks_mortar, mesh = user_data + info_pw = PointerWrapper(info) # Get the global interface/mortar ids and neighbor rank if current face belongs to an MPI # interface/mortar - if unsafe_load(info).sides.elem_count == 2 # MPI interfaces/mortars have two neighboring elements + if info_pw.sides.elem_count[] == 2 # MPI interfaces/mortars have two neighboring elements # Extract surface data - sides = (unsafe_load_side(info, 1), unsafe_load_side(info, 2)) + sides_pw = (load_pointerwrapper_side(info_pw, 1), + load_pointerwrapper_side(info_pw, 2)) - if sides[1].is_hanging == false && sides[2].is_hanging == false # No hanging nodes for MPI interfaces - if sides[1].is.full.is_ghost == true + if sides_pw[1].is_hanging[] == false && sides_pw[2].is_hanging[] == false # No hanging nodes for MPI interfaces + if sides_pw[1].is.full.is_ghost[] == true remote_side = 1 local_side = 2 - elseif sides[2].is.full.is_ghost == true + elseif sides_pw[2].is.full.is_ghost[] == true remote_side = 2 local_side = 1 else # both sides are on this rank -> skip since it's a regular interface @@ -397,16 +397,17 @@ function init_neighbor_rank_connectivity_iter_face_inner(info, user_data) end # Sanity check, current face should belong to current MPI interface - local_tree = unsafe_load_tree(mesh.p4est, sides[local_side].treeid + 1) # one-based indexing - local_quad_id = local_tree.quadrants_offset + - sides[local_side].is.full.quadid + local_tree_pw = load_pointerwrapper_tree(mesh.p4est, + sides_pw[local_side].treeid[] + 1) # one-based indexing + local_quad_id = local_tree_pw.quadrants_offset[] + + sides_pw[local_side].is.full.quadid[] @assert interfaces.local_neighbor_ids[interface_id] == local_quad_id + 1 # one-based indexing # Get neighbor ID from ghost layer proc_offsets = unsafe_wrap(Array, - unsafe_load(unsafe_load(info).ghost_layer).proc_offsets, + info_pw.ghost_layer.proc_offsets, mpi_nranks() + 1) - ghost_id = sides[remote_side].is.full.quadid # indexes the ghost layer, 0-based + ghost_id = sides_pw[remote_side].is.full.quadid[] # indexes the ghost layer, 0-based neighbor_rank = findfirst(r -> proc_offsets[r] <= ghost_id < proc_offsets[r + 1], 1:mpi_nranks()) - 1 # MPI ranks are 0-based @@ -415,21 +416,18 @@ function init_neighbor_rank_connectivity_iter_face_inner(info, user_data) # Global interface id is the globally unique quadrant id of the quadrant on the primary # side (1) multiplied by the number of faces per quadrant plus face if local_side == 1 - offset = unsafe_load(unsafe_load(mesh.p4est).global_first_quadrant, - mpi_rank() + 1) # one-based indexing + offset = mesh.p4est.global_first_quadrant[mpi_rank() + 1] # one-based indexing primary_quad_id = offset + local_quad_id else - offset = unsafe_load(unsafe_load(mesh.p4est).global_first_quadrant, - neighbor_rank + 1) # one-based indexing - primary_quad_id = offset + - unsafe_load(sides[1].is.full.quad.p.piggy3.local_num) + offset = mesh.p4est.global_first_quadrant[neighbor_rank + 1] # one-based indexing + primary_quad_id = offset + sides_pw[1].is.full.quad.p.piggy3.local_num[] end - global_interface_id = 2 * ndims(mesh) * primary_quad_id + sides[1].face + global_interface_id = 2 * ndims(mesh) * primary_quad_id + sides_pw[1].face[] global_interface_ids[interface_id] = global_interface_id user_data.interface_id += 1 else # hanging node - if sides[1].is_hanging == true + if sides_pw[1].is_hanging[] == true hanging_side = 1 full_side = 2 else @@ -437,26 +435,26 @@ function init_neighbor_rank_connectivity_iter_face_inner(info, user_data) full_side = 1 end # Verify before accessing is.full / is.hanging - @assert sides[hanging_side].is_hanging == true && - sides[full_side].is_hanging == false + @assert sides_pw[hanging_side].is_hanging[] == true && + sides_pw[full_side].is_hanging[] == false # If all quadrants are locally available, this is a regular mortar -> skip - if sides[full_side].is.full.is_ghost == false && - all(sides[hanging_side].is.hanging.is_ghost .== false) + if sides_pw[full_side].is.full.is_ghost[] == false && + all(sides_pw[hanging_side].is.hanging.is_ghost[] .== false) return nothing end - trees = (unsafe_load_tree(mesh.p4est, sides[1].treeid + 1), - unsafe_load_tree(mesh.p4est, sides[2].treeid + 1)) + trees_pw = (load_pointerwrapper_tree(mesh.p4est, sides_pw[1].treeid[] + 1), + load_pointerwrapper_tree(mesh.p4est, sides_pw[2].treeid[] + 1)) # Find small quads that are remote and determine which rank owns them - remote_small_quad_positions = findall(sides[hanging_side].is.hanging.is_ghost .== + remote_small_quad_positions = findall(sides_pw[hanging_side].is.hanging.is_ghost[] .== true) proc_offsets = unsafe_wrap(Array, - unsafe_load(unsafe_load(info).ghost_layer).proc_offsets, + info_pw.ghost_layer.proc_offsets, mpi_nranks() + 1) # indices of small remote quads inside the ghost layer, 0-based - ghost_ids = map(pos -> sides[hanging_side].is.hanging.quadid[pos], + ghost_ids = map(pos -> sides_pw[hanging_side].is.hanging.quadid[][pos], remote_small_quad_positions) neighbor_ranks = map(ghost_ids) do ghost_id return findfirst(r -> proc_offsets[r] <= ghost_id < proc_offsets[r + 1], @@ -464,28 +462,26 @@ function init_neighbor_rank_connectivity_iter_face_inner(info, user_data) end # Determine global quad id of large element to determine global MPI mortar id # Furthermore, if large element is ghost, add its owner rank to neighbor_ranks - if sides[full_side].is.full.is_ghost == true - ghost_id = sides[full_side].is.full.quadid + if sides_pw[full_side].is.full.is_ghost[] == true + ghost_id = sides_pw[full_side].is.full.quadid[] large_quad_owner_rank = findfirst(r -> proc_offsets[r] <= ghost_id < proc_offsets[r + 1], 1:mpi_nranks()) - 1 # MPI ranks are 0-based push!(neighbor_ranks, large_quad_owner_rank) - offset = unsafe_load(unsafe_load(mesh.p4est).global_first_quadrant, - large_quad_owner_rank + 1) # one-based indexing + offset = mesh.p4est.global_first_quadrant[large_quad_owner_rank + 1] # one-based indexing large_quad_id = offset + - unsafe_load(sides[full_side].is.full.quad.p.piggy3.local_num) + sides_pw[full_side].is.full.quad.p.piggy3.local_num[] else - offset = unsafe_load(unsafe_load(mesh.p4est).global_first_quadrant, - mpi_rank() + 1) # one-based indexing - large_quad_id = offset + trees[full_side].quadrants_offset + - sides[full_side].is.full.quadid + offset = mesh.p4est.global_first_quadrant[mpi_rank() + 1] # one-based indexing + large_quad_id = offset + trees_pw[full_side].quadrants_offset[] + + sides_pw[full_side].is.full.quadid[] end neighbor_ranks_mortar[mortar_id] = neighbor_ranks # Global mortar id is the globally unique quadrant id of the large quadrant multiplied by the # number of faces per quadrant plus face global_mortar_ids[mortar_id] = 2 * ndims(mesh) * large_quad_id + - sides[full_side].face + sides_pw[full_side].face[] user_data.mortar_id += 1 end From eb4c91a6ba3c1fe1c1f0bbb9c332ab0068171113 Mon Sep 17 00:00:00 2001 From: David Knapp Date: Mon, 19 Jun 2023 10:48:56 +0200 Subject: [PATCH 04/15] Added a dispatchable constructur for DG with TensorProductWedges (#1389) * Added a dispatchable constructur for DG with TensorProductWedges Adapted the timestepping for DGMulti to use the minimal polynomial degree of the tensorproduct * Use the maximal poly-degree * Update src/solvers/dgmulti/types.jl Co-authored-by: Hendrik Ranocha * Change DGMulti-call for Tensorwedges * Adding comment about the polydeg-choice * Adapt constructor * Add an example for TensorWedges * Update example * Update examples/dgmulti_3d/elixir_euler_tensorWedge.jl Co-authored-by: Hendrik Ranocha * Update examples/dgmulti_3d/elixir_euler_tensorWedge.jl Co-authored-by: Hendrik Ranocha * Update examples/dgmulti_3d/elixir_euler_tensorWedge.jl Co-authored-by: Hendrik Ranocha * Readd Gauss Accidentaly deleted during Merge-Conflict resolvement online * Remove explicit load of StartUpDG in elixir * Update max_dt * Apply suggestions from code review Co-authored-by: Hendrik Ranocha * Update comments * Add tensorWedge-test * Update src/solvers/dgmulti/dg.jl Co-authored-by: Jesse Chan <1156048+jlchan@users.noreply.github.com> * Renamed tensorwedge elixir * Update testset title * Remove typos * temporal push * fix max_dt dispatch * Update TensorProductWedges-test * Update test * Address CI-Warning * Update dt_polydeg_scaling * Correct dt_polydeg_scaling A previous commit changed it from (polydeg + 1) to polydeg * Update DGMulti-call * Rename file * Change Path in test-file * Apply suggestions from code review Co-authored-by: Hendrik Ranocha * Update examples/dgmulti_3d/elixir_advection_tensor_wedge.jl Co-authored-by: Hendrik Ranocha * Reformat src/Trixi.jl * Format solvers/dgmulti/dg.jl * Formatting --------- Co-authored-by: Knapp Co-authored-by: Hendrik Ranocha Co-authored-by: Jesse Chan <1156048+jlchan@users.noreply.github.com> Co-authored-by: Jesse Chan --- .../elixir_advection_tensor_wedge.jl | 56 +++++++++++++++++++ src/Trixi.jl | 3 +- src/solvers/dgmulti/dg.jl | 11 ++-- src/solvers/dgmulti/types.jl | 15 +++++ test/test_dgmulti_3d.jl | 6 ++ 5 files changed, 86 insertions(+), 5 deletions(-) create mode 100644 examples/dgmulti_3d/elixir_advection_tensor_wedge.jl diff --git a/examples/dgmulti_3d/elixir_advection_tensor_wedge.jl b/examples/dgmulti_3d/elixir_advection_tensor_wedge.jl new file mode 100644 index 00000000000..4f43f2571a3 --- /dev/null +++ b/examples/dgmulti_3d/elixir_advection_tensor_wedge.jl @@ -0,0 +1,56 @@ +using OrdinaryDiffEq +using Trixi +using LinearAlgebra + +############################################################################### +equations = LinearScalarAdvectionEquation3D(1.0, 1.0, 1.0) + +initial_condition = initial_condition_convergence_test + +# Define the polynomial degrees for the polynoms of the triangular base and the line +# of the tensor-prism +tensor_polydeg = (3, 4) + +dg = DGMulti(element_type = Wedge(), + approximation_type = Polynomial(), + surface_flux = flux_lax_friedrichs, + polydeg = tensor_polydeg) + + +cells_per_dimension = (8, 8, 8) +mesh = DGMultiMesh(dg, + cells_per_dimension, + coordinates_min = (-1.0, -1.0, -1.0), + coordinates_max = (1.0, 1.0, 1.0), + periodicity = true) + + +semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, dg, + boundary_conditions=boundary_condition_periodic) + +############################################################################### +# ODE solvers, callbacks etc. + +tspan = (0.0, 5.0) +ode = semidiscretize(semi, tspan) + +summary_callback = SummaryCallback() + +analysis_interval = 100 +analysis_callback = AnalysisCallback(semi, interval=analysis_interval, uEltype=real(dg)) + +alive_callback = AliveCallback(analysis_interval=analysis_interval) + +# The StepsizeCallback handles the re-calculation of the maximum Δt after each time step +stepsize_callback = StepsizeCallback(cfl=1.0) + +callbacks = CallbackSet(summary_callback, analysis_callback, alive_callback, stepsize_callback) + + +############################################################################### +# run the simulation + +sol = solve(ode, CarpenterKennedy2N54(williamson_condition = false), dt = 1.0, + save_everystep=false, callback=callbacks); + +summary_callback() # print the timer summary \ No newline at end of file diff --git a/src/Trixi.jl b/src/Trixi.jl index 86e349c7dad..66878f4b459 100644 --- a/src/Trixi.jl +++ b/src/Trixi.jl @@ -83,7 +83,8 @@ import SummationByPartsOperators: integrate, semidiscretize, upwind_operators # DGMulti solvers -@reexport using StartUpDG: StartUpDG, Polynomial, Gauss, SBP, Line, Tri, Quad, Hex, Tet +@reexport using StartUpDG: StartUpDG, Polynomial, Gauss, TensorProductWedge, SBP, Line, Tri, + Quad, Hex, Tet, Wedge using StartUpDG: RefElemData, MeshData, AbstractElemShape # TODO: include_optimized diff --git a/src/solvers/dgmulti/dg.jl b/src/solvers/dgmulti/dg.jl index d51c7cabf9d..bc76aa1a9d2 100644 --- a/src/solvers/dgmulti/dg.jl +++ b/src/solvers/dgmulti/dg.jl @@ -226,6 +226,11 @@ function estimate_dt(mesh::DGMultiMesh, dg::DGMulti) return StartUpDG.estimate_h(rd, mesh.md) / StartUpDG.inverse_trace_constant(rd) end +dt_polydeg_scaling(dg::DGMulti) = inv(dg.basis.N + 1) +function dt_polydeg_scaling(dg::DGMulti{3, <:Wedge, <:TensorProductWedge}) + inv(maximum(dg.basis.N) + 1) +end + # for the stepsize callback function max_dt(u, t, mesh::DGMultiMesh, constant_speed::False, equations, dg::DGMulti{NDIMS}, @@ -247,8 +252,7 @@ function max_dt(u, t, mesh::DGMultiMesh, # `polydeg+1`. This is because `nnodes(dg)` returns the total number of # multi-dimensional nodes for DGMulti solver types, while `nnodes(dg)` returns # the number of 1D nodes for `DGSEM` solvers. - polydeg = rd.N - return 2 * dt_min / (polydeg + 1) + return 2 * dt_min * dt_polydeg_scaling(dg) end function max_dt(u, t, mesh::DGMultiMesh, @@ -270,8 +274,7 @@ function max_dt(u, t, mesh::DGMultiMesh, # `polydeg+1`. This is because `nnodes(dg)` returns the total number of # multi-dimensional nodes for DGMulti solver types, while `nnodes(dg)` returns # the number of 1D nodes for `DGSEM` solvers. - polydeg = rd.N - return 2 * dt_min / (polydeg + 1) + return 2 * dt_min * dt_polydeg_scaling(dg) end # interpolates from solution coefficients to face quadrature points diff --git a/src/solvers/dgmulti/types.jl b/src/solvers/dgmulti/types.jl index c225e334e8e..f1f7b158dec 100644 --- a/src/solvers/dgmulti/types.jl +++ b/src/solvers/dgmulti/types.jl @@ -96,6 +96,21 @@ function DGMulti(; polydeg = nothing, polydeg = polydeg, kwargs...) end +# dispatchable constructor for DGMulti using a TensorProductWedge +function DGMulti(element_type::Wedge, + approximation_type, + volume_integral, + surface_integral; + polydeg::Tuple, + kwargs...) + factor_a = RefElemData(Tri(), approximation_type, polydeg[1]; kwargs...) + factor_b = RefElemData(Line(), approximation_type, polydeg[2]; kwargs...) + + tensor = TensorProductWedge(factor_a, factor_b) + rd = RefElemData(element_type, tensor; kwargs...) + return DG(rd, nothing, surface_integral, volume_integral) +end + # dispatchable constructor for DGMulti to allow for specialization function DGMulti(element_type::AbstractElemShape, approximation_type, diff --git a/test/test_dgmulti_3d.jl b/test/test_dgmulti_3d.jl index 22c0a0fd3ba..68fa1d13304 100644 --- a/test/test_dgmulti_3d.jl +++ b/test/test_dgmulti_3d.jl @@ -135,6 +135,12 @@ isdir(outdir) && rm(outdir, recursive=true) ) end + @trixi_testset "elixir_advection_tensor_wedge.jl" begin + @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_advection_tensor_wedge.jl"), + l2 = [2.30487910e-04] , + linf = [6.31795281e-04] ) + end + end # Clean up afterwards: delete Trixi.jl output directory From a837dd7726f25a88c7bab8bed2dc6955f20534bf Mon Sep 17 00:00:00 2001 From: Andrew Winters Date: Mon, 19 Jun 2023 14:50:27 +0200 Subject: [PATCH 05/15] Update test values for 1D discontinuous examples (#1511) * set true discontinuous initial conditions for 1D * update Burgers tests * update MHD tests in 1D * update structured_1d_dgsem/elixir_advection_shockcapturing elixir and test * remove workaround in elixir_shallowwater_ec.jl and update test * accidentally removed adjusted tolerance in structured mesh test * remove workaround in elixir_shallowwater_shockcapturing.jl and update test * remove workaround in elixir_shallowwater_well_balanced.jl and update tests * bug fix in energy_total computation for ShallowWaterTwoLayerEquations1D * remove workaround in elixir_shallowwater_twolayer_dam_break.jl and update test * update docs test for adding a nonconservative equation * update Nonconservative terms in 1D (linear advection) test * update all necessary IdealGlmMhdMulticomponentEquations1D tests * update all necessary CompressibleEulerMulticomponentEquations1D tests * update necessary CompressibleEuler1D tests except for neural network shock capturing * Update test values for NN-based SC * add short description to the NEWS.md * add extra test to avoid coverage drop * Update src/solvers/dg.jl Co-authored-by: Hendrik Ranocha --------- Co-authored-by: Michael Schlottke-Lakemper Co-authored-by: Hendrik Ranocha --- NEWS.md | 1 + .../files/adding_nonconservative_equation.jl | 2 +- .../elixir_advection_shockcapturing.jl | 6 +- .../tree_1d_dgsem/elixir_burgers_shock.jl | 6 +- ..._eulermulti_two_interacting_blast_waves.jl | 2 +- .../elixir_mhd_briowu_shock_tube.jl | 2 +- examples/tree_1d_dgsem/elixir_mhdmulti_ec.jl | 4 +- examples/tree_1d_dgsem/elixir_mhdmulti_es.jl | 4 +- .../tree_1d_dgsem/elixir_shallowwater_ec.jl | 71 ++++++--------- .../elixir_shallowwater_shock_capturing.jl | 81 +++++++---------- .../elixir_shallowwater_twolayer_dam_break.jl | 77 +++++++---------- .../elixir_shallowwater_well_balanced.jl | 56 +++--------- src/equations/shallow_water_two_layer_1d.jl | 86 +++++++++---------- src/solvers/dg.jl | 9 ++ test/test_structured_1d.jl | 4 +- test/test_tree_1d.jl | 2 +- test/test_tree_1d_burgers.jl | 8 +- test/test_tree_1d_euler.jl | 44 +++++----- test/test_tree_1d_eulermulti.jl | 32 ++++--- test/test_tree_1d_mhd.jl | 12 +-- test/test_tree_1d_mhdmulti.jl | 48 +++++------ test/test_tree_1d_shallowwater.jl | 28 +++--- test/test_tree_1d_shallowwater_twolayer.jl | 6 +- 23 files changed, 267 insertions(+), 324 deletions(-) diff --git a/NEWS.md b/NEWS.md index 9b46ba565fe..35c7039b2ef 100644 --- a/NEWS.md +++ b/NEWS.md @@ -9,6 +9,7 @@ for human readability. #### Added - Experimental support for 3D parabolic diffusion terms has been added. +- Capability to set truly discontinuous initial conditions in 1D. #### Changed diff --git a/docs/literate/src/files/adding_nonconservative_equation.jl b/docs/literate/src/files/adding_nonconservative_equation.jl index 08dd631058e..110fa486070 100644 --- a/docs/literate/src/files/adding_nonconservative_equation.jl +++ b/docs/literate/src/files/adding_nonconservative_equation.jl @@ -147,7 +147,7 @@ plot(sol) # above. error_1 = analysis_callback(sol).l2 |> first -@test isapprox(error_1, 0.0002961027497) #src +@test isapprox(error_1, 0.00029609575838969394) #src # Next, we increase the grid resolution by one refinement level and run the # simulation again. diff --git a/examples/structured_1d_dgsem/elixir_advection_shockcapturing.jl b/examples/structured_1d_dgsem/elixir_advection_shockcapturing.jl index 9a81acfe51c..313812fe08d 100644 --- a/examples/structured_1d_dgsem/elixir_advection_shockcapturing.jl +++ b/examples/structured_1d_dgsem/elixir_advection_shockcapturing.jl @@ -9,7 +9,9 @@ advection_velocity = 1.0 """ initial_condition_composite(x, t, equations::LinearScalarAdvectionEquation1D) -Slight simplification of +Wave form that is a combination of a Gaussian pulse, a square wave, a triangle wave, +and half an ellipse with periodic boundary conditions. +Slight simplification from - Jiang, Shu (1996) Efficient Implementation of Weighted ENO Schemes [DOI: 10.1006/jcph.1996.0130](https://doi.org/10.1006/jcph.1996.0130) @@ -60,7 +62,7 @@ volume_integral = VolumeIntegralShockCapturingHG(indicator_sc; solver = DGSEM(basis, surface_flux, volume_integral) # Create curved mesh -cells_per_dimension = (125,) +cells_per_dimension = (120,) coordinates_min = (-1.0,) # minimum coordinate coordinates_max = (1.0,) # maximum coordinate mesh = StructuredMesh(cells_per_dimension, coordinates_min, coordinates_max, diff --git a/examples/tree_1d_dgsem/elixir_burgers_shock.jl b/examples/tree_1d_dgsem/elixir_burgers_shock.jl index 987fb320ad6..00b5314e19f 100644 --- a/examples/tree_1d_dgsem/elixir_burgers_shock.jl +++ b/examples/tree_1d_dgsem/elixir_burgers_shock.jl @@ -21,7 +21,7 @@ surface_flux = flux_lax_friedrichs volume_integral = VolumeIntegralShockCapturingHG(indicator_sc; volume_flux_dg=surface_flux, volume_flux_fv=surface_flux) - + solver = DGSEM(basis, surface_flux, volume_integral) coordinate_min = 0.0 @@ -59,7 +59,7 @@ end boundary_conditions = (x_neg=boundary_condition_inflow, x_pos=boundary_condition_outflow) - + initial_condition = initial_condition_shock semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver, @@ -79,7 +79,7 @@ analysis_callback = AnalysisCallback(semi, interval=analysis_interval) alive_callback = AliveCallback(analysis_interval=analysis_interval) -stepsize_callback = StepsizeCallback(cfl=0.9) +stepsize_callback = StepsizeCallback(cfl=0.85) callbacks = CallbackSet(summary_callback, diff --git a/examples/tree_1d_dgsem/elixir_eulermulti_two_interacting_blast_waves.jl b/examples/tree_1d_dgsem/elixir_eulermulti_two_interacting_blast_waves.jl index 353093e5f70..81966194180 100644 --- a/examples/tree_1d_dgsem/elixir_eulermulti_two_interacting_blast_waves.jl +++ b/examples/tree_1d_dgsem/elixir_eulermulti_two_interacting_blast_waves.jl @@ -88,7 +88,7 @@ ode = semidiscretize(semi, tspan) summary_callback = SummaryCallback() -analysis_interval = 100 +analysis_interval = 1000 analysis_callback = AnalysisCallback(semi, interval=analysis_interval) diff --git a/examples/tree_1d_dgsem/elixir_mhd_briowu_shock_tube.jl b/examples/tree_1d_dgsem/elixir_mhd_briowu_shock_tube.jl index 1c07fc4fdde..c5727109d92 100644 --- a/examples/tree_1d_dgsem/elixir_mhd_briowu_shock_tube.jl +++ b/examples/tree_1d_dgsem/elixir_mhd_briowu_shock_tube.jl @@ -94,7 +94,7 @@ amr_callback = AMRCallback(semi, amr_controller, adapt_initial_condition=true, adapt_initial_condition_only_refine=true) -stepsize_callback = StepsizeCallback(cfl=0.8) +stepsize_callback = StepsizeCallback(cfl=0.65) callbacks = CallbackSet(summary_callback, analysis_callback, alive_callback, diff --git a/examples/tree_1d_dgsem/elixir_mhdmulti_ec.jl b/examples/tree_1d_dgsem/elixir_mhdmulti_ec.jl index 34fdce6634e..69ea0551bed 100644 --- a/examples/tree_1d_dgsem/elixir_mhdmulti_ec.jl +++ b/examples/tree_1d_dgsem/elixir_mhdmulti_ec.jl @@ -4,8 +4,8 @@ using Trixi ############################################################################### # semidiscretization of the ideal MHD equations -equations = IdealGlmMhdMulticomponentEquations1D(gammas = (2.0, 2.0, 2.0), - gas_constants = (2.0, 2.0, 2.0)) +equations = IdealGlmMhdMulticomponentEquations1D(gammas = (2.0, 2.0, 2.0), + gas_constants = (2.0, 2.0, 2.0)) initial_condition = initial_condition_weak_blast_wave diff --git a/examples/tree_1d_dgsem/elixir_mhdmulti_es.jl b/examples/tree_1d_dgsem/elixir_mhdmulti_es.jl index 8ca32194b9e..93cf3e0fdb2 100644 --- a/examples/tree_1d_dgsem/elixir_mhdmulti_es.jl +++ b/examples/tree_1d_dgsem/elixir_mhdmulti_es.jl @@ -4,8 +4,8 @@ using Trixi ############################################################################### # semidiscretization of the ideal MHD equations -equations = IdealGlmMhdMulticomponentEquations1D(gammas = (2.0, 2.0, 2.0), - gas_constants = (2.0, 2.0, 2.0)) +equations = IdealGlmMhdMulticomponentEquations1D(gammas = (2.0, 2.0, 2.0), + gas_constants = (2.0, 2.0, 2.0)) initial_condition = initial_condition_weak_blast_wave diff --git a/examples/tree_1d_dgsem/elixir_shallowwater_ec.jl b/examples/tree_1d_dgsem/elixir_shallowwater_ec.jl index be6a2cb166c..1469afec1ca 100644 --- a/examples/tree_1d_dgsem/elixir_shallowwater_ec.jl +++ b/examples/tree_1d_dgsem/elixir_shallowwater_ec.jl @@ -8,9 +8,34 @@ using Trixi equations = ShallowWaterEquations1D(gravity_constant=9.81) -# Note, this initial condition is used to compute errors in the analysis callback but the initialization is -# overwritten by `initial_condition_ec_discontinuous_bottom` below. -initial_condition = initial_condition_weak_blast_wave +# Initial condition with a truly discontinuous water height, velocity, and bottom +# topography function as an academic testcase for entropy conservation. +# The errors from the analysis callback are not important but `∑∂S/∂U ⋅ Uₜ` should +# be around machine roundoff. +# Works as intended for TreeMesh1D with `initial_refinement_level=4`. If the mesh +# refinement level is changed the initial condition below may need changed as well to +# ensure that the discontinuities lie on an element interface. +function initial_condition_ec_discontinuous_bottom(x, t, equations::ShallowWaterEquations1D) + # Set the background values + H = 4.25 + v = 0.0 + b = sin(x[1]) # arbitrary continuous function + + # Setup the discontinuous water height and velocity + if x[1] >= 0.125 && x[1] <= 0.25 + H = 5.0 + v = 0.1882 + end + + # Setup a discontinuous bottom topography + if x[1] >= -0.25 && x[1] <= -0.125 + b = 2.0 + 0.5 * sin(2.0 * pi * x[1]) + end + + return prim2cons(SVector(H, v, b), equations) +end + +initial_condition = initial_condition_ec_discontinuous_bottom ############################################################################### # Get the DG approximation space @@ -37,46 +62,6 @@ semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver) tspan = (0.0, 2.0) ode = semidiscretize(semi, tspan) -############################################################################### -# Workaround to set a discontinuous bottom topography and initial condition for debugging and testing. - -# alternative version of the initial conditinon used to setup a truly discontinuous -# bottom topography function and initial condition for this academic testcase of entropy conservation. -# The errors from the analysis callback are not important but `∑∂S/∂U ⋅ Uₜ` should be around machine roundoff -# In contrast to the usual signature of initial conditions, this one get passed the -# `element_id` explicitly. In particular, this initial conditions works as intended -# only for the TreeMesh1D with `initial_refinement_level=4`. -function initial_condition_ec_discontinuous_bottom(x, t, element_id, equations::ShallowWaterEquations1D) - # Set the background values - H = 4.25 - v = 0.0 - b = sin(x[1]) # arbitrary continuous function - - # setup the discontinuous water height and velocity - if element_id == 10 - H = 5.0 - v = 0.1882 - end - - # Setup a discontinuous bottom topography using the element id number - if element_id == 7 - b = 2.0 + 0.5 * sin(2.0 * pi * x[1]) - end - - return prim2cons(SVector(H, v, b), equations) -end - -# point to the data we want to augment -u = Trixi.wrap_array(ode.u0, semi) -# reset the initial condition -for element in eachelement(semi.solver, semi.cache) - for i in eachnode(semi.solver) - x_node = Trixi.get_node_coords(semi.cache.elements.node_coordinates, equations, semi.solver, i, element) - u_node = initial_condition_ec_discontinuous_bottom(x_node, first(tspan), element, equations) - Trixi.set_node_vars!(u, u_node, equations, semi.solver, i, element) - end -end - ############################################################################### # Callbacks diff --git a/examples/tree_1d_dgsem/elixir_shallowwater_shock_capturing.jl b/examples/tree_1d_dgsem/elixir_shallowwater_shock_capturing.jl index 50241126a28..62346d7b5ab 100644 --- a/examples/tree_1d_dgsem/elixir_shallowwater_shock_capturing.jl +++ b/examples/tree_1d_dgsem/elixir_shallowwater_shock_capturing.jl @@ -7,24 +7,37 @@ using Trixi equations = ShallowWaterEquations1D(gravity_constant=9.812, H0=1.75) -function initial_condition_stone_throw(x, t, equations::ShallowWaterEquations1D) - # Set up polar coordinates - inicenter = 0.15 - x_norm = x[1] - inicenter[1] - r = abs(x_norm) +# Initial condition with a truly discontinuous velocity and bottom topography. +# Works as intended for TreeMesh1D with `initial_refinement_level=3`. If the mesh +# refinement level is changed the initial condition below may need changed as well to +# ensure that the discontinuities lie on an element interface. +function initial_condition_stone_throw_discontinuous_bottom(x, t, equations::ShallowWaterEquations1D) + + # Calculate primitive variables + + # flat lake + H = equations.H0 + + # Discontinuous velocity + v = 0.0 + if x[1] >= -0.75 && x[1] <= 0.0 + v = -1.0 + elseif x[1] >= 0.0 && x[1] <= 0.75 + v = 1.0 + end - # Calculate primitive variables - H = equations.H0 - # v = 0.0 # for well-balanced test - v = r < 0.6 ? 1.75 : 0.0 # for stone throw + b = ( 1.5 / exp( 0.5 * ((x[1] - 1.0)^2 ) ) + + 0.75 / exp( 0.5 * ((x[1] + 1.0)^2 ) ) ) - b = ( 1.5 / exp( 0.5 * ((x[1] - 1.0)^2 ) ) - + 0.75 / exp( 0.5 * ((x[1] + 1.0)^2 ) ) ) + # Force a discontinuous bottom topography + if x[1] >= -1.5 && x[1] <= 0.0 + b = 0.5 + end - return prim2cons(SVector(H, v, b), equations) + return prim2cons(SVector(H, v, b), equations) end -initial_condition = initial_condition_stone_throw +initial_condition = initial_condition_stone_throw_discontinuous_bottom boundary_condition = boundary_condition_slip_wall @@ -62,49 +75,13 @@ semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver, boundary_conditions = boundary_condition) ############################################################################### -# ODE solvers, callbacks etc. +# ODE solver tspan = (0.0, 3.0) ode = semidiscretize(semi, tspan) -# Hack in a discontinuous bottom for a more interesting test -function initial_condition_stone_throw_discontinuous_bottom(x, t, element_id, equations::ShallowWaterEquations1D) - - inicenter = 0.15 - x_norm = x[1] - inicenter[1] - r = abs(x_norm) - - # Calculate primitive variables - H = equations.H0 # flat lake - # Discontinuous velocity set via element id number - v = 0.0 - if element_id == 4 - v = -1.0 - elseif element_id == 5 - v = 1.0 - end - - b = ( 1.5 / exp( 0.5 * ((x[1] - 1.0)^2 ) ) - + 0.75 / exp( 0.5 * ((x[1] + 1.0)^2 ) ) ) - - # Setup a discontinuous bottom topography using the element id number - if element_id == 3 || element_id == 4 - b = 0.5 - end - - return prim2cons(SVector(H, v, b), equations) -end - -# point to the data we want to augment -u = Trixi.wrap_array(ode.u0, semi) -# reset the initial condition -for element in eachelement(semi.solver, semi.cache) - for i in eachnode(semi.solver) - x_node = Trixi.get_node_coords(semi.cache.elements.node_coordinates, equations, semi.solver, i, element) - u_node = initial_condition_stone_throw_discontinuous_bottom(x_node, first(tspan), element, equations) - Trixi.set_node_vars!(u, u_node, equations, semi.solver, i, element) - end -end +############################################################################### +# Callbacks summary_callback = SummaryCallback() diff --git a/examples/tree_1d_dgsem/elixir_shallowwater_twolayer_dam_break.jl b/examples/tree_1d_dgsem/elixir_shallowwater_twolayer_dam_break.jl index 67c1098b178..60770d158fa 100644 --- a/examples/tree_1d_dgsem/elixir_shallowwater_twolayer_dam_break.jl +++ b/examples/tree_1d_dgsem/elixir_shallowwater_twolayer_dam_break.jl @@ -3,20 +3,34 @@ using OrdinaryDiffEq using Trixi ############################################################################### -# Semidiscretization of the two-layer shallow water equations for a dam break test with a -# discontinuous bottom topography function to test entropy conservation +# Semidiscretization of the two-layer shallow water equations for a dam break +# test with a discontinuous bottom topography function to test entropy conservation equations = ShallowWaterTwoLayerEquations1D(gravity_constant=9.81, H0=2.0, rho_upper=0.9, rho_lower=1.0) -############################################################################### -# Workaround to set a discontinuous bottom topography and initial condition. +# Initial condition of a dam break with a discontinuous water heights and bottom topography. +# Works as intended for TreeMesh1D with `initial_refinement_level=5`. If the mesh +# refinement level is changed the initial condition below may need changed as well to +# ensure that the discontinuities lie on an element interface. +function initial_condition_dam_break(x, t, equations::ShallowWaterTwoLayerEquations1D) + v1_upper = 0.0 + v1_lower = 0.0 + + # Set the discontinuity + if x[1] <= 10.0 + H_lower = 2.0 + H_upper = 4.0 + b = 0.0 + else + H_lower = 1.5 + H_upper = 3.0 + b = 0.5 + end + + return prim2cons(SVector(H_upper, v1_upper, H_lower, v1_lower, b), equations) +end -# This test case uses a special work around to setup a truly discontinuous bottom topography -# function and initial condition for this academic testcase of entropy conservation. First, a -# dummy initial condition is introduced to create the semidiscretization. Then the initial condition -# is reset with the true discontinuous values from initial_condition_dam_break. Note, that this -# initial condition only works for TreeMesh1D with `initial_refinement_level=5`. -initial_condition = initial_condition_convergence_test +initial_condition = initial_condition_dam_break ############################################################################### # Get the DG approximation space @@ -25,7 +39,6 @@ volume_flux = (flux_wintermeyer_etal, flux_nonconservative_wintermeyer_etal) solver = DGSEM(polydeg=3, surface_flux=(flux_fjordholm_etal, flux_nonconservative_fjordholm_etal), volume_integral=VolumeIntegralFluxDifferencing(volume_flux)) - ############################################################################### # Get the TreeMesh and setup a non-periodic mesh @@ -39,54 +52,22 @@ mesh = TreeMesh(coordinates_min, coordinates_max, boundary_condition = boundary_condition_slip_wall # create the semidiscretization object -semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver, +semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver, boundary_conditions=boundary_condition) ############################################################################### -# ODE solvers, callbacks etc. +# ODE solvers tspan = (0.0,0.4) ode = semidiscretize(semi, tspan) -# Initial conditions dam break test case -function initial_condition_dam_break(x, t, element_id, equations::ShallowWaterTwoLayerEquations1D) - v1_upper = 0.0 - v1_lower = 0.0 - - # Set the discontinuity - if element_id <= 16 - H_lower = 2.0 - H_upper = 4.0 - b = 0.0 - else - H_lower = 1.5 - H_upper = 3.0 - b = 0.5 - end - - return prim2cons(SVector(H_upper, v1_upper, H_lower, v1_lower, b), equations) -end - - -# point to the data we want to augment -u = Trixi.wrap_array(ode.u0, semi) -# reset the initial condition -for element in eachelement(semi.solver, semi.cache) - for i in eachnode(semi.solver) - x_node = Trixi.get_node_coords(semi.cache.elements.node_coordinates, - equations, semi.solver, i, element) - u_node = initial_condition_dam_break(x_node, first(tspan), element, equations) - Trixi.set_node_vars!(u, u_node, equations, semi.solver, i, element) - end -end - - - +############################################################################### +# Callbacks summary_callback = SummaryCallback() analysis_interval = 500 -analysis_callback = AnalysisCallback(semi, interval=analysis_interval, save_analysis=false, +analysis_callback = AnalysisCallback(semi, interval=analysis_interval, save_analysis=false, extra_analysis_integrals=(energy_total, energy_kinetic, energy_internal,)) stepsize_callback = StepsizeCallback(cfl=1.0) diff --git a/examples/tree_1d_dgsem/elixir_shallowwater_well_balanced.jl b/examples/tree_1d_dgsem/elixir_shallowwater_well_balanced.jl index 83835656839..e07bc04d76a 100644 --- a/examples/tree_1d_dgsem/elixir_shallowwater_well_balanced.jl +++ b/examples/tree_1d_dgsem/elixir_shallowwater_well_balanced.jl @@ -8,21 +8,28 @@ using Trixi equations = ShallowWaterEquations1D(gravity_constant=9.81, H0=3.25) -# An initial condition with constant total water height and zero velocities to test well-balancedness. -# Note, this routine is used to compute errors in the analysis callback but the initialization is -# overwritten by `initial_condition_discontinuous_well_balancedness` below. -function initial_condition_well_balancedness(x, t, equations::ShallowWaterEquations1D) +# Setup a truly discontinuous bottom topography function for this academic +# testcase of well-balancedness. The errors from the analysis callback are +# not important but the error for this lake-at-rest test case +# `∑|H0-(h+b)|` should be around machine roundoff. +# Works as intended for TreeMesh1D with `initial_refinement_level=3`. If the mesh +# refinement level is changed the initial condition below may need changed as well to +# ensure that the discontinuities lie on an element interface. +function initial_condition_discontinuous_well_balancedness(x, t, equations::ShallowWaterEquations1D) # Set the background values H = equations.H0 v = 0.0 + b = 0.0 - # bottom topography inspired by from Pond.control in [HOHQMesh](https://github.com/trixi-framework/HOHQMesh) - b = (1.5 / exp( 0.5 * ((x[1] - 1.0)^2) )+ 0.75 / exp(0.5 * ((x[1] + 1.0)^2))) + # Setup a discontinuous bottom topography + if x[1] >= 0.5 && x[1] <= 0.75 + b = 2.0 + 0.5 * sin(2.0 * pi * x[1]) + end return prim2cons(SVector(H, v, b), equations) end -initial_condition = initial_condition_well_balancedness +initial_condition = initial_condition_discontinuous_well_balancedness ############################################################################### # Get the DG approximation space @@ -50,41 +57,6 @@ semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver) tspan = (0.0, 100.0) ode = semidiscretize(semi, tspan) -############################################################################### -# Workaround to set a discontinuous bottom topography and initial condition for debugging and testing. - -# alternative version of the initial conditinon used to setup a truly discontinuous -# bottom topography function for this academic testcase of well-balancedness. -# The errors from the analysis callback are not important but the error for this lake at rest test case -# `∑|H0-(h+b)|` should be around machine roundoff. -# In contrast to the usual signature of initial conditions, this one get passed the -# `element_id` explicitly. In particular, this initial conditions works as intended -# only for the TreeMesh1D with `initial_refinement_level=3`. -function initial_condition_discontinuous_well_balancedness(x, t, element_id, equations::ShallowWaterEquations1D) - # Set the background values - H = equations.H0 - v = 0.0 - b = 0.0 - - # Setup a discontinuous bottom topography using the element id number - if element_id == 7 - b = 2.0 + 0.5 * sin(2.0 * pi * x[1]) - end - - return prim2cons(SVector(H, v, b), equations) -end - -# point to the data we want to augment -u = Trixi.wrap_array(ode.u0, semi) -# reset the initial condition -for element in eachelement(semi.solver, semi.cache) - for i in eachnode(semi.solver) - x_node = Trixi.get_node_coords(semi.cache.elements.node_coordinates, equations, semi.solver, i, element) - u_node = initial_condition_discontinuous_well_balancedness(x_node, first(tspan), element, equations) - Trixi.set_node_vars!(u, u_node, equations, semi.solver, i, element) - end -end - ############################################################################### # Callbacks diff --git a/src/equations/shallow_water_two_layer_1d.jl b/src/equations/shallow_water_two_layer_1d.jl index edf7d5e32ff..02899171509 100644 --- a/src/equations/shallow_water_two_layer_1d.jl +++ b/src/equations/shallow_water_two_layer_1d.jl @@ -11,28 +11,28 @@ Two-Layer Shallow Water equations (2LSWE) in one space dimension. The equations are given by ```math \begin{alignat*}{4} -&\frac{\partial}{\partial t}h_{upper} -&&+ \frac{\partial}{\partial x}\left(h_{upper} v_{1,upper}\right) +&\frac{\partial}{\partial t}h_{upper} +&&+ \frac{\partial}{\partial x}\left(h_{upper} v_{1,upper}\right) &&= 0 \\ -&\frac{\partial}{\partial t}\left(h_{upper}v_{1,upper}\right) -&&+ \frac{\partial}{\partial x}\left(h_{upper}v_{1,upper}^2 + \dfrac{gh_{upper}^2}{2}\right) +&\frac{\partial}{\partial t}\left(h_{upper}v_{1,upper}\right) +&&+ \frac{\partial}{\partial x}\left(h_{upper}v_{1,upper}^2 + \dfrac{gh_{upper}^2}{2}\right) &&= -gh_{upper}\frac{\partial}{\partial x}\left(b+h_{lower}\right)\\ -&\frac{\partial}{\partial t}h_{lower} -&&+ \frac{\partial}{\partial x}\left(h_{lower}v_{1,lower}\right) +&\frac{\partial}{\partial t}h_{lower} +&&+ \frac{\partial}{\partial x}\left(h_{lower}v_{1,lower}\right) &&= 0 \\ -&\frac{\partial}{\partial t}\left(h_{lower}v_{1,lower}\right) -&&+ \frac{\partial}{\partial x}\left(h_{lower}v_{1,lower}^2 + \dfrac{gh_{lower}^2}{2}\right) +&\frac{\partial}{\partial t}\left(h_{lower}v_{1,lower}\right) +&&+ \frac{\partial}{\partial x}\left(h_{lower}v_{1,lower}^2 + \dfrac{gh_{lower}^2}{2}\right) &&= -gh_{lower}\frac{\partial}{\partial x}\left(b+\dfrac{\rho_{upper}}{\rho_{lower}}h_{upper}\right). \end{alignat*} ``` -The unknown quantities of the 2LSWE are the water heights of the {lower} layer ``h_{lower}`` and the -{upper} layer ``h_{upper}`` with respective velocities ``v_{1,upper}`` and ``v_{1,lower}``. The gravitational constant is -denoted by `g`, the layer densitites by ``\rho_{upper}``and ``\rho_{lower}`` and the (possibly) variable -bottom topography function ``b(x)``. The conservative variable water height ``h_{lower}`` is measured -from the bottom topography ``b`` and ``h_{upper}`` relative to ``h_{lower}``, therefore one also defines the +The unknown quantities of the 2LSWE are the water heights of the {lower} layer ``h_{lower}`` and the +{upper} layer ``h_{upper}`` with respective velocities ``v_{1,upper}`` and ``v_{1,lower}``. The gravitational constant is +denoted by `g`, the layer densitites by ``\rho_{upper}``and ``\rho_{lower}`` and the (possibly) variable +bottom topography function ``b(x)``. The conservative variable water height ``h_{lower}`` is measured +from the bottom topography ``b`` and ``h_{upper}`` relative to ``h_{lower}``, therefore one also defines the total water heights as ``H_{upper} = h_{upper} + h_{upper} + b`` and ``H_{lower} = h_{lower} + b``. -The densities must be chosen such that ``\rho_{upper} < \rho_{lower}``, to make sure that the heavier fluid +The densities must be chosen such that ``\rho_{upper} < \rho_{lower}``, to make sure that the heavier fluid ``\rho_{lower}`` is in the bottom layer and the lighter fluid ``\rho_{upper}`` in the {upper} layer. The additional quantity ``H_0`` is also available to store a reference value for the total water @@ -41,13 +41,13 @@ height that is useful to set initial conditions or test the "lake-at-rest" well- The bottom topography function ``b(x)`` is set inside the initial condition routine for a particular problem setup. -In addition to the unknowns, Trixi currently stores the bottom topography values at the -approximation points despite being fixed in time. This is done for convenience of computing the -bottom topography gradients on the fly during the approximation as well as computing auxiliary +In addition to the unknowns, Trixi currently stores the bottom topography values at the +approximation points despite being fixed in time. This is done for convenience of computing the +bottom topography gradients on the fly during the approximation as well as computing auxiliary quantities like the total water height ``H`` or the entropy variables. This affects the implementation and use of these equations in various ways: * The flux values corresponding to the bottom topography must be zero. -* The bottom topography values must be included when defining initial conditions, boundary +* The bottom topography values must be included when defining initial conditions, boundary conditions or source terms. * [`AnalysisCallback`](@ref) analyzes this variable. * Trixi's visualization tools will visualize the bottom topography by default. @@ -101,7 +101,7 @@ end initial_condition_convergence_test(x, t, equations::ShallowWaterTwoLayerEquations1D) A smooth initial condition used for convergence tests in combination with -[`source_terms_convergence_test`](@ref) (and +[`source_terms_convergence_test`](@ref) (and [`BoundaryConditionDirichlet(initial_condition_convergence_test)`](@ref) in non-periodic domains). """ function initial_condition_convergence_test(x, t, @@ -121,9 +121,9 @@ end """ source_terms_convergence_test(u, x, t, equations::ShallowWaterTwoLayerEquations1D) -Source terms used for convergence tests in combination with -[`initial_condition_convergence_test`](@ref) -(and [`BoundaryConditionDirichlet(initial_condition_convergence_test)`](@ref) +Source terms used for convergence tests in combination with +[`initial_condition_convergence_test`](@ref) +(and [`BoundaryConditionDirichlet(initial_condition_convergence_test)`](@ref) in non-periodic domains). """ @inline function source_terms_convergence_test(u, x, t, @@ -167,8 +167,8 @@ the internal value. For details see Section 9.2.5 of the book: - Eleuterio F. Toro (2001) - Shock-Capturing Methods for Free-Surface Shallow Flows - 1st edition + Shock-Capturing Methods for Free-Surface Shallow Flows + 1st edition ISBN 0471987662 """ @inline function boundary_condition_slip_wall(u_inner, orientation_or_normal, direction, @@ -219,7 +219,7 @@ end Non-symmetric two-point volume flux discretizing the nonconservative (source) term that contains the gradient of the bottom topography [`ShallowWaterTwoLayerEquations2D`](@ref) and an -additional term that couples the momentum of both layers. This is a slightly modified version +additional term that couples the momentum of both layers. This is a slightly modified version to account for the additional source term compared to the standard SWE described in the paper. Further details are available in the paper: @@ -238,7 +238,7 @@ Further details are available in the paper: z = zero(eltype(u_ll)) - # Bottom gradient nonconservative term: (0, g*h_upper*(b+h_lower)_x, + # Bottom gradient nonconservative term: (0, g*h_upper*(b+h_lower)_x, # 0, g*h_lower*(b+r*h_upper)_x, 0) f = SVector(z, equations.gravity * h_upper_ll * (b_rr + h_lower_rr), @@ -254,9 +254,9 @@ end !!! warning "Experimental code" This numerical flux is experimental and may change in any future release. - -Non-symmetric two-point surface flux discretizing the nonconservative (source) term that contains -the gradients of the bottom topography and an additional term that couples the momentum of both + +Non-symmetric two-point surface flux discretizing the nonconservative (source) term that contains +the gradients of the bottom topography and an additional term that couples the momentum of both layers [`ShallowWaterTwoLayerEquations2D`](@ref). Further details are available in the paper: @@ -286,7 +286,7 @@ formulation. z = zero(eltype(u_ll)) - # Bottom gradient nonconservative term: (0, g*h_upper*(b+h_lower)_x, + # Bottom gradient nonconservative term: (0, g*h_upper*(b+h_lower)_x, # 0, g*h_lower*(b+r*h_upper)_x, 0) f = SVector(z, g * h_upper_ll * (b_ll + h_lower_ll) + @@ -303,13 +303,13 @@ end flux_fjordholm_etal(u_ll, u_rr, orientation, equations::ShallowWaterTwoLayerEquations1D) -Total energy conservative (mathematical entropy for shallow water equations). When the bottom -topography is nonzero this should only be used as a surface flux otherwise the scheme will not be +Total energy conservative (mathematical entropy for shallow water equations). When the bottom +topography is nonzero this should only be used as a surface flux otherwise the scheme will not be well-balanced. For well-balancedness in the volume flux use [`flux_wintermeyer_etal`](@ref). Details are available in Eq. (4.1) in the paper: - Ulrik S. Fjordholm, Siddhartha Mishra and Eitan Tadmor (2011) - Well-balanced and energy stable schemes for the shallow water equations with discontinuous + Well-balanced and energy stable schemes for the shallow water equations with discontinuous topography [DOI: 10.1016/j.jcp.2011.03.042](https://doi.org/10.1016/j.jcp.2011.03.042) and the application to two layers is shown in the paper: - Ulrik Skre Fjordholm (2012) @@ -348,11 +348,11 @@ end """ flux_wintermeyer_etal(u_ll, u_rr, orientation, equations::ShallowWaterTwoLayerEquations1D) - + Total energy conservative (mathematical entropy for two-layer shallow water equations) split form. When the bottom topography is nonzero this scheme will be well-balanced when used as a `volume_flux`. The `surface_flux` should still use, e.g., [`flux_fjordholm_etal`](@ref). To obtain the flux for the -two-layer shallow water equations the flux that is described in the paper for the normal shallow +two-layer shallow water equations the flux that is described in the paper for the normal shallow water equations is used within each layer. Further details are available in Theorem 1 of the paper: @@ -391,8 +391,8 @@ end flux_es_fjordholm_etal(u_ll, u_rr, orientation, equations::ShallowWaterTwoLayerEquations1D) -Entropy stable surface flux for the two-layer shallow water equations. Uses the entropy -conservative flux_fjordholm_etal and adds a Lax-Friedrichs type dissipation dependent on the jump +Entropy stable surface flux for the two-layer shallow water equations. Uses the entropy +conservative flux_fjordholm_etal and adds a Lax-Friedrichs type dissipation dependent on the jump of entropy variables. Further details are available in the paper: @@ -460,10 +460,10 @@ formulation. end # Calculate approximation for maximum wave speed for local Lax-Friedrichs-type dissipation as the -# maximum velocity magnitude plus the maximum speed of sound. This function uses approximate -# eigenvalues using the speed of the barotropic mode as there is no simple way to calculate them +# maximum velocity magnitude plus the maximum speed of sound. This function uses approximate +# eigenvalues using the speed of the barotropic mode as there is no simple way to calculate them # analytically. -# +# # A good overview of the derivation is given in: # - Jonas Nycander, Andrew McC. Hogg, Leela M. Frankcombe (2008) # Open boundary conditions for nonlinear channel Flows @@ -488,7 +488,7 @@ end return (max(abs(v_m_ll) + c_ll, abs(v_m_rr) + c_rr)) end -# Specialized `DissipationLocalLaxFriedrichs` to avoid spurious dissipation in the bottom +# Specialized `DissipationLocalLaxFriedrichs` to avoid spurious dissipation in the bottom # topography @inline function (dissipation::DissipationLocalLaxFriedrichs)(u_ll, u_rr, orientation_or_normal_direction, @@ -530,7 +530,7 @@ end end # Convert conservative variables to entropy variables -# Note, only the first four are the entropy variables, the fifth entry still just carries the +# Note, only the first four are the entropy variables, the fifth entry still just carries the # bottom topography values for convenience @inline function cons2entropy(u, equations::ShallowWaterTwoLayerEquations1D) h_upper, _, h_lower, _, b = u @@ -567,7 +567,7 @@ end # Calculate total energy for a conservative state `cons` @inline function energy_total(cons, equations::ShallowWaterTwoLayerEquations1D) - h_upper, h_lower, h_v1_upper, h_v2_lower, b = cons + h_upper, h_v1_upper, h_lower, h_v2_lower, b = cons # Set new variables for better readability g = equations.gravity rho_upper = equations.rho_upper diff --git a/src/solvers/dg.jl b/src/solvers/dg.jl index 838fa2d5819..2536cfe0bf2 100644 --- a/src/solvers/dg.jl +++ b/src/solvers/dg.jl @@ -628,6 +628,15 @@ function compute_coefficients!(u, func, t, mesh::AbstractMesh{1}, equations, dg: for i in eachnode(dg) x_node = get_node_coords(cache.elements.node_coordinates, equations, dg, i, element) + # Changing the node positions passed to the initial condition by the minimum + # amount possible with the current type of floating point numbers allows setting + # discontinuous initial data in a simple way. In particular, a check like `if x < x_jump` + # works if the jump location `x_jump` is at the position of an interface. + if i == 1 + x_node = SVector(nextfloat(x_node[1])) + elseif i == nnodes(dg) + x_node = SVector(prevfloat(x_node[1])) + end u_node = func(x_node, t, equations) set_node_vars!(u, u_node, equations, dg, i, element) end diff --git a/test/test_structured_1d.jl b/test/test_structured_1d.jl index a27d3c219e1..ec8c7a138d5 100644 --- a/test/test_structured_1d.jl +++ b/test/test_structured_1d.jl @@ -27,8 +27,8 @@ isdir(outdir) && rm(outdir, recursive=true) @trixi_testset "elixir_advection_shockcapturing.jl" begin @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_advection_shockcapturing.jl"), - l2 = [0.08004076716881656], - linf = [0.6342577638501385], + l2 = [0.08015029105233593], + linf = [0.610709468736576], atol = 1.0e-5) end diff --git a/test/test_tree_1d.jl b/test/test_tree_1d.jl index e37e1efc3e6..7737a93a15a 100644 --- a/test/test_tree_1d.jl +++ b/test/test_tree_1d.jl @@ -271,7 +271,7 @@ end sol = solve(ode, Tsit5(), abstol=1.0e-6, reltol=1.0e-6, save_everystep=false, callback=callbacks); - @test analysis_callback(sol).l2 ≈ [0.00029610274971929974, 5.573684084938363e-6] + @test analysis_callback(sol).l2 ≈ [0.00029609575838969394, 5.5681704039507985e-6] end diff --git a/test/test_tree_1d_burgers.jl b/test/test_tree_1d_burgers.jl index 788c7ab4199..8c4cfaa406d 100644 --- a/test/test_tree_1d_burgers.jl +++ b/test/test_tree_1d_burgers.jl @@ -22,14 +22,14 @@ EXAMPLES_DIR = pkgdir(Trixi, "examples", "tree_1d_dgsem") @trixi_testset "elixir_burgers_shock.jl" begin @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_burgers_shock.jl"), - l2 = [0.4407585104869119], - linf = [1.000000000000001]) + l2 = [0.4422505602587537], + linf = [1.0000000000000009]) end @trixi_testset "elixir_burgers_rarefaction.jl" begin @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_burgers_rarefaction.jl"), - l2 = [0.40287062735307044], - linf = [1.0042992585765542]) + l2 = [0.4038224690923722], + linf = [1.0049201454652736]) end end diff --git a/test/test_tree_1d_euler.jl b/test/test_tree_1d_euler.jl index 40f2a38b0e1..5fb74b80bce 100644 --- a/test/test_tree_1d_euler.jl +++ b/test/test_tree_1d_euler.jl @@ -22,8 +22,8 @@ EXAMPLES_DIR = pkgdir(Trixi, "examples", "tree_1d_dgsem") @trixi_testset "elixir_euler_density_wave.jl" begin @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_density_wave.jl"), - l2 = [0.0011482554820185795, 0.00011482554830363504, 5.741277417754598e-6], - linf = [0.004090978306820037, 0.00040909783134346345, 2.0454891732413216e-5]) + l2 = [0.0011482554820217855, 0.00011482554830323462, 5.741277429325267e-6], + linf = [0.004090978306812376, 0.0004090978313582294, 2.045489210189544e-5]) end @trixi_testset "elixir_euler_density_wave.jl with initial_condition_constant" begin @@ -41,14 +41,14 @@ EXAMPLES_DIR = pkgdir(Trixi, "examples", "tree_1d_dgsem") @trixi_testset "elixir_euler_ec.jl" begin @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_ec.jl"), - l2 = [0.11915540925414216, 0.15489191247295198, 0.44543052524765375], - linf = [0.2751485868543495, 0.2712764982000735, 0.9951407418216425]) + l2 = [0.11821957357197649, 0.15330089521538678, 0.4417674632047301], + linf = [0.24280567569982958, 0.29130548795961936, 0.8847009003152442]) end @trixi_testset "elixir_euler_ec.jl with flux_kennedy_gruber" begin @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_ec.jl"), - l2 = [0.07905582221868049, 0.10180958900546237, 0.29596551476711125], - linf = [0.23515297345769826, 0.2958208108392532, 0.8694224308790321], + l2 = [0.07803455838661963, 0.10032577312032283, 0.29228156303827935], + linf = [0.2549869853794955, 0.3376472164661263, 0.9650477546553962], maxiters = 10, surface_flux = flux_kennedy_gruber, volume_flux = flux_kennedy_gruber) @@ -56,8 +56,8 @@ EXAMPLES_DIR = pkgdir(Trixi, "examples", "tree_1d_dgsem") @trixi_testset "elixir_euler_ec.jl with flux_shima_etal" begin @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_ec.jl"), - l2 = [0.07909267609417114, 0.1018246500951966, 0.2959649187481973], - linf = [0.23631829743146504, 0.2977756307879202, 0.8642794698697331], + l2 = [0.07800654460172655, 0.10030365573277883, 0.2921481199111959], + linf = [0.25408579350400395, 0.3388657679031271, 0.9776486386921928], maxiters = 10, surface_flux = flux_shima_etal, volume_flux = flux_shima_etal) @@ -65,8 +65,8 @@ EXAMPLES_DIR = pkgdir(Trixi, "examples", "tree_1d_dgsem") @trixi_testset "elixir_euler_ec.jl with flux_chandrashekar" begin @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_ec.jl"), - l2 = [0.07905306555214126, 0.10181180378499956, 0.2959171937479504], - linf = [0.24057642004451651, 0.29691454643616433, 0.886425723870524], + l2 = [0.07801923089205756, 0.10039557434912669, 0.2922210399923278], + linf = [0.2576521982607225, 0.3409717926625057, 0.9772961936567048], maxiters = 10, surface_flux = flux_chandrashekar, volume_flux = flux_chandrashekar) @@ -74,8 +74,8 @@ EXAMPLES_DIR = pkgdir(Trixi, "examples", "tree_1d_dgsem") @trixi_testset "elixir_euler_ec.jl with flux_hll" begin @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_ec.jl"), - l2 = [0.07959780803600519, 0.10342491934977621, 0.2978851659149904], - linf = [0.19228754121840885, 0.2524152253292552, 0.725604944702432], + l2 = [0.07852272782240548, 0.10209790867523805, 0.293873048809011], + linf = [0.19244768908604093, 0.2515941686151897, 0.7258000837553769], maxiters = 10, surface_flux = flux_hll, volume_flux = flux_ranocha) @@ -83,8 +83,8 @@ EXAMPLES_DIR = pkgdir(Trixi, "examples", "tree_1d_dgsem") @trixi_testset "elixir_euler_shockcapturing.jl" begin @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_shockcapturing.jl"), - l2 = [0.11665968950973675, 0.15105507394693413, 0.43503082674771115], - linf = [0.1867400345208743, 0.24621854448555328, 0.703826406555577]) + l2 = [0.11606096465319675, 0.15028768943458806, 0.4328230323046703], + linf = [0.18031710091067965, 0.2351582421501841, 0.6776805692092567]) end @trixi_testset "elixir_euler_sedov_blast_wave.jl" begin @@ -96,8 +96,8 @@ EXAMPLES_DIR = pkgdir(Trixi, "examples", "tree_1d_dgsem") @trixi_testset "elixir_euler_sedov_blast_wave_pure_fv.jl" begin @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_sedov_blast_wave_pure_fv.jl"), - l2 = [1.075075094036344, 0.06766902169711514, 0.9221426570128292], - linf = [3.3941512671408542, 0.16862631133303882, 2.6572394126490315], + l2 = [1.0735456065491455, 0.07131078703089379, 0.9205739468590453], + linf = [3.4296365168219216, 0.17635583964559245, 2.6574584326179505], # Let this test run longer to cover some lines in flux_hllc coverage_override = (maxiters=10^5, tspan=(0.0, 0.1))) end @@ -129,22 +129,22 @@ EXAMPLES_DIR = pkgdir(Trixi, "examples", "tree_1d_dgsem") @trixi_testset "elixir_euler_blast_wave.jl" begin @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_blast_wave.jl"), - l2 = [0.21651329948737183, 0.28091709900008616, 0.5580778880050432], - linf = [1.513525457073142, 1.5328754303137992, 2.0467706106669556], + l2 = [0.21934822867340323, 0.28131919126002686, 0.554361702716662], + linf = [1.5180897390290355, 1.3967085956620369, 2.0663825294019595], maxiters = 30) end @trixi_testset "elixir_euler_blast_wave_neuralnetwork_perssonperaire.jl" begin @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_blast_wave_neuralnetwork_perssonperaire.jl"), - l2 = [2.13605618e-01, 2.79953055e-01, 5.54424459e-01], - linf = [1.55151701e+00, 1.55696782e+00, 2.05525953e+00], + l2 = [0.21814833203212694, 0.2818328665444332, 0.5528379124720818], + linf = [1.5548653877320868, 1.4474018998129738, 2.071919577393772], maxiters = 30) end @trixi_testset "elixir_euler_blast_wave_neuralnetwork_rayhesthaven.jl" begin @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_blast_wave_neuralnetwork_rayhesthaven.jl"), - l2 = [2.18148857e-01, 2.83182959e-01, 5.59096194e-01], - linf = [1.62706876e+00, 1.61680275e+00, 2.05876517e+00], + l2 = [0.22054468879127423, 0.2828269190680846, 0.5542369885642424], + linf = [1.5623359741479623, 1.4290121654488288, 2.1040405133123072], maxiters = 30) end end diff --git a/test/test_tree_1d_eulermulti.jl b/test/test_tree_1d_eulermulti.jl index eac54a9372c..e880f98e2d0 100644 --- a/test/test_tree_1d_eulermulti.jl +++ b/test/test_tree_1d_eulermulti.jl @@ -11,39 +11,47 @@ EXAMPLES_DIR = pkgdir(Trixi, "examples", "tree_1d_dgsem") @trixi_testset "elixir_eulermulti_ec.jl" begin @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_eulermulti_ec.jl"), - l2 = [1.54891912e-01, 4.45430525e-01, 1.70222013e-02, 3.40444026e-02, 6.80888053e-02], - linf = [2.71276498e-01, 9.95140742e-01, 3.93069410e-02, 7.86138820e-02, 1.57227764e-01]) + l2 = [0.15330089521538684, 0.4417674632047301, 0.016888510510282385, 0.03377702102056477, + 0.06755404204112954], + linf = [0.29130548795961864, 0.8847009003152357, 0.034686525099975274, 0.06937305019995055, + 0.1387461003999011]) end @trixi_testset "elixir_eulermulti_es.jl" begin @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_eulermulti_es.jl"), - l2 = [1.53387916e-01, 4.41585576e-01, 3.93605635e-02, 7.87211270e-02], - linf = [2.49632117e-01, 7.21088064e-01, 6.38328770e-02, 1.27665754e-01]) + l2 = [0.1522380497572071, 0.43830846465313206, 0.03907262116499431, 0.07814524232998862], + linf = [0.24939193075537294, 0.7139395740052739, 0.06324208768391237, 0.12648417536782475]) end @trixi_testset "elixir_eulermulti_convergence_ec.jl" begin @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_eulermulti_convergence_ec.jl"), - l2 = [8.57523604e-05, 1.63878043e-04, 1.94126993e-05, 3.88253986e-05], - linf = [3.05932773e-04, 6.24480393e-04, 7.25312144e-05, 1.45062429e-04]) + l2 = [8.575236038539227e-5, 0.00016387804318585358, 1.9412699303977585e-5, 3.882539860795517e-5], + linf = [0.00030593277277124464, 0.0006244803933350696, 7.253121435135679e-5, 0.00014506242870271358]) end @trixi_testset "elixir_eulermulti_convergence_es.jl" begin @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_eulermulti_convergence_es.jl"), - l2 = [1.8983933794407234e-5, 6.207744299844731e-5, 1.5466205761868047e-6, 3.0932411523736094e-6, 6.186482304747219e-6, 1.2372964609494437e-5], - linf = [0.00012014372605895218, 0.0003313207215800418, 6.50836791016296e-6, 1.301673582032592e-5, 2.603347164065184e-5, 5.206694328130368e-5]) + l2 = [1.8983933794407234e-5, 6.207744299844731e-5, 1.5466205761868047e-6, 3.0932411523736094e-6, + 6.186482304747219e-6, 1.2372964609494437e-5], + linf = [0.00012014372605895218, 0.0003313207215800418, 6.50836791016296e-6, 1.301673582032592e-5, + 2.603347164065184e-5, 5.206694328130368e-5]) end @trixi_testset "elixir_eulermulti_convergence_es.jl with flux_chandrashekar" begin @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_eulermulti_convergence_es.jl"), - l2 = [1.88845048e-05, 5.49106005e-05, 9.42673716e-07, 1.88534743e-06, 3.77069486e-06, 7.54138973e-06], - linf = [1.16223512e-04, 3.07922197e-04, 3.21774233e-06, 6.43548465e-06, 1.28709693e-05, 2.57419386e-05], + l2 = [1.888450477353845e-5, 5.4910600482795386e-5, 9.426737161533622e-7, 1.8853474323067245e-6, + 3.770694864613449e-6, 7.541389729226898e-6], + linf = [0.00011622351152063004, 0.0003079221967086099, 3.2177423254231563e-6, 6.435484650846313e-6, + 1.2870969301692625e-5, 2.574193860338525e-5], volume_flux = flux_chandrashekar) end @trixi_testset "elixir_eulermulti_two_interacting_blast_waves.jl" begin @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_eulermulti_two_interacting_blast_waves.jl"), - l2 = [1.28886761e+00, 8.27133526e+01, 3.50680272e-03, 1.36987844e-02, 1.91795185e-02], - linf = [2.96413045e+01, 1.32258448e+03, 9.19191937e-02, 3.10929710e-01, 4.41798976e-01], + l2 = [1.288867611915533, 82.71335258388848, 0.00350680272313187, 0.013698784353152794, + 0.019179518517518084], + linf = [29.6413044707026, 1322.5844802186496, 0.09191919374782143, 0.31092970966717925, + 0.4417989757182038], tspan = (0.0, 0.0001)) end diff --git a/test/test_tree_1d_mhd.jl b/test/test_tree_1d_mhd.jl index 938959831c1..e3a0cda3250 100644 --- a/test/test_tree_1d_mhd.jl +++ b/test/test_tree_1d_mhd.jl @@ -32,14 +32,14 @@ EXAMPLES_DIR = pkgdir(Trixi, "examples", "tree_1d_dgsem") @trixi_testset "elixir_mhd_ec.jl" begin @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_mhd_ec.jl"), - l2 = [5.86009540e-02, 8.16048158e-02, 5.46791194e-02, 5.46791194e-02, 1.54509265e-01, 4.13046273e-17, 5.47637521e-02, 5.47637521e-02], - linf = [1.10014999e-01, 1.81982581e-01, 9.13611439e-02, 9.13611439e-02, 4.23831370e-01, 1.11022302e-16, 9.93731761e-02, 9.93731761e-02]) + l2 = [0.05815183849746399, 0.08166807325621023, 0.054659228513541165, 0.054659228513541165, 0.15578125987042743, 4.130462730494e-17, 0.05465258887150046, 0.05465258887150046], + linf = [0.12165312668363826, 0.1901920742264952, 0.10059813883022554, 0.10059813883022554, 0.44079257431070706, 1.1102230246251565e-16, 0.10528911365809579, 0.10528911365809579]) end @trixi_testset "elixir_mhd_briowu_shock_tube.jl" begin @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_mhd_briowu_shock_tube.jl"), - l2 = [0.17764301067932906, 0.19693621875378622, 0.3635136528288642, 0.0, 0.3757321708837591, 8.593007507325741e-16, 0.36473438378159656, 0.0], - linf = [0.5601530250396535, 0.43867368105486537, 1.0960903616351099, 0.0, 1.0551794137886303, 4.107825191113079e-15, 1.5374410890043144, 0.0], + l2 = [0.17477712356961989, 0.19489623595086944, 0.3596546157640463, 0.0, 0.3723215736814466, 1.2060075775846403e-15, 0.36276754492568164, 0.0], + linf = [0.5797109945880677, 0.4372991899547103, 1.0906536287185835, 0.0, 1.0526758874956808, 5.995204332975845e-15, 1.5122922036932964, 0.0], coverage_override = (maxiters=6,)) end @@ -51,8 +51,8 @@ EXAMPLES_DIR = pkgdir(Trixi, "examples", "tree_1d_dgsem") @trixi_testset "elixir_mhd_ryujones_shock_tube.jl" begin @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_mhd_ryujones_shock_tube.jl"), - l2 = [2.34809441e-01, 3.92255943e-01, 8.23575546e-02, 1.75599624e-01, 9.61613519e-01, 6.60825891e-17, 2.15346454e-01, 1.07006529e-01], - linf = [6.40732148e-01, 9.44889516e-01, 3.54932707e-01, 8.54060243e-01, 2.07757711e+00, 1.11022302e-16, 4.92584725e-01, 2.49526561e-01], + l2 = [0.23469781891518154, 0.3916675299696121, 0.08245195301016353, 0.1745346945706147, 0.9606363432904367, 6.608258910237605e-17, 0.21542929107153735, 0.10705457908737925], + linf = [0.6447951791685409, 0.9461857095377463, 0.35074627554617605, 0.8515177411529542, 2.0770652030507053, 1.1102230246251565e-16, 0.49670855513788204, 0.24830199967863564], tspan = (0.0, 0.1)) end diff --git a/test/test_tree_1d_mhdmulti.jl b/test/test_tree_1d_mhdmulti.jl index 2985e6d5663..5214ed26d38 100644 --- a/test/test_tree_1d_mhdmulti.jl +++ b/test/test_tree_1d_mhdmulti.jl @@ -11,33 +11,33 @@ EXAMPLES_DIR = pkgdir(Trixi, "examples", "tree_1d_dgsem") @trixi_testset "elixir_mhdmulti_ec.jl" begin @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_mhdmulti_ec.jl"), - l2 = [0.08160481582829862, 0.05467911944326103, 0.05467911944326103, 0.15450926504459692, - 4.130462730494e-17, 0.054763752050210085, 0.054763752050210085, 0.008371564857135208, - 0.016743129714270416, 0.03348625942854083], - linf = [0.18198258075330706, 0.09136114386311774, 0.09136114386311774, 0.423831369951313, - 1.1102230246251565e-16, 0.09937317613143604, 0.09937317613143604, 0.0157164284712992, - 0.0314328569425984, 0.0628657138851968]) + l2 = [0.08166807325620999, 0.054659228513541616, 0.054659228513541616, 0.15578125987042812, + 4.130462730494e-17, 0.054652588871500665, 0.054652588871500665, 0.008307405499637766, + 0.01661481099927553, 0.03322962199855106], + linf = [0.19019207422649645, 0.10059813883022888, 0.10059813883022888, 0.4407925743107146, + 1.1102230246251565e-16, 0.10528911365809623, 0.10528911365809623, 0.01737901809766182, + 0.03475803619532364, 0.06951607239064728]) end @trixi_testset "elixir_mhdmulti_ec.jl with flux_derigs_etal" begin @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_mhdmulti_ec.jl"), - l2 = [0.08153372259925547, 0.05464109003345891, 0.05464109003345891, 0.1540576724164453, - 4.130462730494e-17, 0.054734930802131036, 0.054734930802131036, 0.008391254781284321, - 0.016782509562568642, 0.033565019125137284], - linf = [0.17492544007323832, 0.09029632168248182, 0.09029632168248182, 0.40798609353896564, - 1.1102230246251565e-16, 0.09872923637833075, 0.09872923637833075, 0.01609818847160674, - 0.03219637694321348, 0.06439275388642696], + l2 = [0.08151404166186461, 0.054640238302693274, 0.054640238302693274, 0.15536125426328573, + 4.130462730494e-17, 0.054665489963920275, 0.054665489963920275, 0.008308349501359825, + 0.01661669900271965, 0.0332333980054393], + linf = [0.1824424257860952, 0.09734687137001484, 0.09734687137001484, 0.4243089502087325, + 1.1102230246251565e-16, 0.09558639591092555, 0.09558639591092555, 0.017364773041550624, + 0.03472954608310125, 0.0694590921662025], volume_flux = flux_derigs_etal) end @trixi_testset "elixir_mhdmulti_es.jl" begin @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_mhdmulti_es.jl"), - l2 = [0.07968782477167513, 0.05398115008116676, 0.05398115008116676, 0.15015281822439228, - 4.130462730494e-17, 0.053629890024921495, 0.053629890024921495, 0.008279068245579706, - 0.016558136491159413, 0.033116272982318826], - linf = [0.14118014632124837, 0.07820697032983395, 0.07820697032983395, 0.3390558674728652, - 1.1102230246251565e-16, 0.06998787893467828, 0.06998787893467828, 0.014943825414763745, - 0.02988765082952749, 0.05977530165905498]) + l2 = [0.07994082660130175, 0.053940174914031976, 0.053940174914031976, 0.15165513559250643, + 4.130462730494e-17, 0.05363207135290325, 0.05363207135290325, 0.008258265884659555, + 0.01651653176931911, 0.03303306353863822], + linf = [0.14101014428198477, 0.07762441749521025, 0.07762441749521025, 0.3381334453289866, + 1.1102230246251565e-16, 0.07003646400675223, 0.07003646400675223, 0.014962483760600165, + 0.02992496752120033, 0.05984993504240066]) end @trixi_testset "elixir_mhdmulti_convergence.jl" begin @@ -52,12 +52,12 @@ EXAMPLES_DIR = pkgdir(Trixi, "examples", "tree_1d_dgsem") @trixi_testset "elixir_mhdmulti_briowu_shock_tube.jl" begin @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_mhdmulti_briowu_shock_tube.jl"), - l2 = [0.1946577804333822, 0.3591196215528672, 0.0, 0.36875476066849383, - 4.7644020131827105e-16, 0.36668249926193885, 0.0, 0.05775369214541893, - 0.11550738429083786], - linf = [0.4345551123140612, 1.0874941615375844, 0.0, 1.0493729052116585, - 3.219646771412954e-15, 1.5160434573973656, 0.0, 0.18616213071936066, - 0.3723242614387213], + l2 = [0.1877830835572639, 0.3455841730726793, 0.0, 0.35413123388836687, + 8.745556626531982e-16, 0.3629920109231055, 0.0, 0.05329005553971236, + 0.10658011107942472], + linf = [0.4288187627971754, 1.0386547815614993, 0.0, 0.9541678878162702, + 5.773159728050814e-15, 1.4595119339458051, 0.0, 0.18201910908829552, + 0.36403821817659104], coverage_override = (maxiters=6,)) end diff --git a/test/test_tree_1d_shallowwater.jl b/test/test_tree_1d_shallowwater.jl index f8901a3dcb6..1c3bac1fab6 100644 --- a/test/test_tree_1d_shallowwater.jl +++ b/test/test_tree_1d_shallowwater.jl @@ -10,22 +10,30 @@ EXAMPLES_DIR = pkgdir(Trixi, "examples", "tree_1d_dgsem") @testset "Shallow Water" begin @trixi_testset "elixir_shallowwater_ec.jl" begin @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_shallowwater_ec.jl"), - l2 = [0.8122354510732459, 1.01586214815876, 0.43404255061704217], - linf = [1.4883285368551107, 3.8717508164234276, 1.7711213427919539], + l2 = [0.244729018751225, 0.8583565222389505, 0.07330427577586297], + linf = [2.1635021283528504, 3.8717508164234453, 1.7711213427919539], + tspan = (0.0, 0.25)) + end + + @trixi_testset "elixir_shallowwater_ec.jl with initial_condition_weak_blast_wave" begin + @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_shallowwater_ec.jl"), + l2 = [0.39464782107209717, 2.03880864210846, 4.1623084150546725e-10], + linf = [0.778905801278281, 3.2409883402608273, 7.419800190922032e-10], + initial_condition=initial_condition_weak_blast_wave, tspan = (0.0, 0.25)) end @trixi_testset "elixir_shallowwater_well_balanced.jl" begin @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_shallowwater_well_balanced.jl"), - l2 = [1.2427984842961743, 1.0332499675061871e-14, 1.2427984842961741], - linf = [1.619041478244762, 1.266865149831811e-14, 1.6190414782447629], + l2 = [0.10416666834254829, 1.4352935256803184e-14, 0.10416666834254838], + linf = [1.9999999999999996, 3.248036646353028e-14, 2.0], tspan = (0.0, 0.25)) end @trixi_testset "elixir_shallowwater_well_balanced.jl with FluxHydrostaticReconstruction" begin @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_shallowwater_well_balanced.jl"), - l2 = [1.2427984842961743, 1.2663646513352053e-14, 1.2427984842961741], - linf = [1.619041478244762, 2.4566658711604395e-14, 1.6190414782447629], + l2 = [0.10416666834254835, 1.1891029971551825e-14, 0.10416666834254838], + linf = [2.0000000000000018, 2.4019608337954543e-14, 2.0], surface_flux=(FluxHydrostaticReconstruction(flux_lax_friedrichs, hydrostatic_reconstruction_audusse_etal), flux_nonconservative_audusse_etal), tspan = (0.0, 0.25)) end @@ -59,14 +67,14 @@ EXAMPLES_DIR = pkgdir(Trixi, "examples", "tree_1d_dgsem") tspan = (0.0, 0.025)) end - @trixi_testset "elixir_shallowwater_well_balanced_nonperiodic.jl with dirichlet boundary" begin + @trixi_testset "elixir_shallowwater_well_balanced_nonperiodic.jl with Dirichlet boundary" begin @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_shallowwater_well_balanced_nonperiodic.jl"), l2 = [1.725964362045055e-8, 5.0427180314307505e-16, 1.7259643530442137e-8], linf = [3.844551077492042e-8, 3.469453422316143e-15, 3.844551077492042e-8], tspan = (0.0, 0.25)) end - @trixi_testset "elixir_shallowwater_well_nonperiodic.jl with wall boundary" begin + @trixi_testset "elixir_shallowwater_well_balanced_nonperiodic.jl with wall boundary" begin @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_shallowwater_well_balanced_nonperiodic.jl"), l2 = [1.7259643614361866e-8, 3.5519018243195145e-16, 1.7259643530442137e-8], linf = [3.844551010878661e-8, 9.846474508971374e-16, 3.844551077492042e-8], @@ -76,8 +84,8 @@ EXAMPLES_DIR = pkgdir(Trixi, "examples", "tree_1d_dgsem") @trixi_testset "elixir_shallowwater_shock_capturing.jl" begin @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_shallowwater_shock_capturing.jl"), - l2 = [0.2884024818919076, 0.5252262013521178, 0.2890348477852955], - linf = [0.7565706154863958, 2.076621603471687, 0.8646939843534258], + l2 = [0.07424140641160326, 0.2148642632748155, 0.0372579849000542], + linf = [1.1209754279344226, 1.3230788645853582, 0.8646939843534251], tspan = (0.0, 0.05)) end end diff --git a/test/test_tree_1d_shallowwater_twolayer.jl b/test/test_tree_1d_shallowwater_twolayer.jl index 6c0ad2941cc..0d8a83806f9 100644 --- a/test/test_tree_1d_shallowwater_twolayer.jl +++ b/test/test_tree_1d_shallowwater_twolayer.jl @@ -38,9 +38,9 @@ EXAMPLES_DIR = pkgdir(Trixi, "examples", "tree_1d_dgsem") @trixi_testset "elixir_shallowwater_twolayer_dam_break.jl with flux_lax_friedrichs" begin @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_shallowwater_twolayer_dam_break.jl"), - l2 = [0.35490827242437256, 1.6715402155795918, 0.6960264969949427, - 0.9351481433409805, 0.7938172946965545], - linf = [0.6417127471419837, 1.9742107034120873, 1.135774587483082, 1.236125279347084, 1.1], + l2 = [0.10010269243463918, 0.5668733957648654, 0.08759617327649398, + 0.4538443183566172, 0.013638618139749523], + linf = [0.5854202777756559, 2.1278930820498934, 0.5193686074348809, 1.8071213168086229, 0.5], surface_flux = (flux_lax_friedrichs, flux_nonconservative_fjordholm_etal), tspan = (0.0, 0.25)) end From 9ea37cce8f510dcc135ad568280fd770a693ba10 Mon Sep 17 00:00:00 2001 From: David Knapp Date: Mon, 19 Jun 2023 16:16:48 +0200 Subject: [PATCH 06/15] Add a pre-commit script and an adapted formatting file (#1534) * Add a pre-commit script and an adapted formatting file * Formatting I am aware of the irony * Fix typo * Added documentation * Fix Typo --- docs/src/styleguide.md | 11 ++++++++++ utils/pre-commit | 40 ++++++++++++++++++++++++++++++++++++ utils/trixi-format-file.jl | 42 ++++++++++++++++++++++++++++++++++++++ 3 files changed, 93 insertions(+) create mode 100755 utils/pre-commit create mode 100755 utils/trixi-format-file.jl diff --git a/docs/src/styleguide.md b/docs/src/styleguide.md index de367c086cc..9e6b6a8c265 100644 --- a/docs/src/styleguide.md +++ b/docs/src/styleguide.md @@ -65,3 +65,14 @@ utils/trixi-format.jl ``` You can get more information about using the convenience script by running it with the `--help`/`-h` flag. + +### Checking formatting before committing +It can be convenient to check the formatting of source code automatically before each commit. +We use git-hooks for it and provide a `pre-commit` script in the `utils` folder. The script uses +[JuliaFormatter.jl](https://github.com/domluna/JuliaFormatter.jl) just like formatting script that +runs over the whole Trixi.jl directory. +You can copy the `pre-commit`-script into `.git/hooks/pre-commit` and it will check your formatting +before each commit. If errors are found the commit is aborted and you can add the corrections via +```shell +git add -p +``` \ No newline at end of file diff --git a/utils/pre-commit b/utils/pre-commit new file mode 100755 index 00000000000..2977b9a200b --- /dev/null +++ b/utils/pre-commit @@ -0,0 +1,40 @@ +#!/bin/sh + +# Copy this file into .git/hooks/pre-commit to execute before each commit. +# It checks and corrects the format for each file. +# If incorrect formatting is found you can add the correction via git add -p + +echo "Checking format before committing" + +if git ref-parse --verify HEAD >/dev/null 2>&1 +then + against=HEAD +else + # Initial commit: diff against an empty tree object + against=280fc57fade28e35046c3e884e587ffef05d3867 +fi + +# Redirect output to stderr. +exec 1>&2 + +# Create a list of files to format. +files=() + +for file in `git diff --cached --name-only` +do + # only indent existing files, this is necessary since if we rename or delete + # a file it is added to the committed files and we thus would try to indent a + # nonexisting file. + if [ ! -e $file ] + then + continue + fi + # We only indent .jl files + FILE_ENDING="${file##*.}" + if [ $FILE_ENDING = "jl" ] + then + files+=($file) + fi +done + +julia utils/trixi-format-file.jl "${files[@]}" diff --git a/utils/trixi-format-file.jl b/utils/trixi-format-file.jl new file mode 100755 index 00000000000..c4d8e7c9032 --- /dev/null +++ b/utils/trixi-format-file.jl @@ -0,0 +1,42 @@ +#!/usr/bin/env julia + +using Pkg +Pkg.activate(; temp = true, io = devnull) +Pkg.add("JuliaFormatter"; preserve = PRESERVE_ALL, io = devnull) + +using JuliaFormatter: format_file + +function main() + # Show help + if "-h" in ARGS || "--help" in ARGS + println("usage: trixi-format.jl PATH [PATH...]") + println() + println("positional arguments:") + println() + println(" PATH One or more paths (directories or files) to format. Default: '.'") + return nothing + end + + file_list = ARGS + if isempty(ARGS) + exit(0) + end + non_formatted_files = Vector{String}() + for file in file_list + println("Checking file " * file) + if !format_file(file) + push!(non_formatted_files, file) + end + end + if isempty(non_formatted_files) + exit(0) + else + @error "Some files have not been formatted! Formatting has been applied, run 'git add -p' to update changes." + for file in non_formatted_files + println(file) + end + exit(1) + end +end + +main() From 11f6fa786340534c10f6356886e19d6291cf618a Mon Sep 17 00:00:00 2001 From: Hendrik Ranocha Date: Mon, 19 Jun 2023 16:57:26 +0200 Subject: [PATCH 07/15] set version to v0.5.29 --- Project.toml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Project.toml b/Project.toml index 9d51e4dcffc..303e97d6324 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.5.29-pre" +version = "0.5.29" [deps] CodeTracking = "da1fd8a2-8d9e-5ec2-8556-3022fb5608a2" From e84946146c022588ef2d4260e4fa4c34b1ab2d9d Mon Sep 17 00:00:00 2001 From: Hendrik Ranocha Date: Mon, 19 Jun 2023 16:57:48 +0200 Subject: [PATCH 08/15] set development version to v0.5.30-pre --- Project.toml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Project.toml b/Project.toml index 303e97d6324..d3983262591 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.5.29" +version = "0.5.30-pre" [deps] CodeTracking = "da1fd8a2-8d9e-5ec2-8556-3022fb5608a2" From 9519c26b2fb88398ab12ddc4f3c6acc62c02a426 Mon Sep 17 00:00:00 2001 From: "dependabot[bot]" <49699333+dependabot[bot]@users.noreply.github.com> Date: Mon, 19 Jun 2023 21:24:37 +0200 Subject: [PATCH 09/15] Bump crate-ci/typos from 1.15.0 to 1.15.1 (#1537) Bumps [crate-ci/typos](https://github.com/crate-ci/typos) from 1.15.0 to 1.15.1. - [Release notes](https://github.com/crate-ci/typos/releases) - [Changelog](https://github.com/crate-ci/typos/blob/master/CHANGELOG.md) - [Commits](https://github.com/crate-ci/typos/compare/v1.15.0...v1.15.1) --- updated-dependencies: - dependency-name: crate-ci/typos dependency-type: direct:production update-type: version-update:semver-patch ... Signed-off-by: dependabot[bot] Co-authored-by: dependabot[bot] <49699333+dependabot[bot]@users.noreply.github.com> --- .github/workflows/SpellCheck.yml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/.github/workflows/SpellCheck.yml b/.github/workflows/SpellCheck.yml index bc324c689bc..75886465f85 100644 --- a/.github/workflows/SpellCheck.yml +++ b/.github/workflows/SpellCheck.yml @@ -10,4 +10,4 @@ jobs: - name: Checkout Actions Repository uses: actions/checkout@v3 - name: Check spelling - uses: crate-ci/typos@v1.15.0 + uses: crate-ci/typos@v1.15.1 From 830d1d7a5fcc76d9ee205c72abff15f328372d02 Mon Sep 17 00:00:00 2001 From: Joshua Lampert <51029046+JoshuaLampert@users.noreply.github.com> Date: Tue, 20 Jun 2023 09:02:05 +0200 Subject: [PATCH 10/15] Fix pre-commit (#1536) * fix pre-commit * consistent indent --- utils/pre-commit | 34 +++++++++++++++++----------------- 1 file changed, 17 insertions(+), 17 deletions(-) diff --git a/utils/pre-commit b/utils/pre-commit index 2977b9a200b..73ad061baef 100755 --- a/utils/pre-commit +++ b/utils/pre-commit @@ -1,4 +1,4 @@ -#!/bin/sh +#!/bin/bash # Copy this file into .git/hooks/pre-commit to execute before each commit. # It checks and corrects the format for each file. @@ -8,10 +8,10 @@ echo "Checking format before committing" if git ref-parse --verify HEAD >/dev/null 2>&1 then - against=HEAD + against=HEAD else - # Initial commit: diff against an empty tree object - against=280fc57fade28e35046c3e884e587ffef05d3867 + # Initial commit: diff against an empty tree object + against=280fc57fade28e35046c3e884e587ffef05d3867 fi # Redirect output to stderr. @@ -22,19 +22,19 @@ files=() for file in `git diff --cached --name-only` do - # only indent existing files, this is necessary since if we rename or delete - # a file it is added to the committed files and we thus would try to indent a - # nonexisting file. - if [ ! -e $file ] - then - continue - fi - # We only indent .jl files - FILE_ENDING="${file##*.}" - if [ $FILE_ENDING = "jl" ] - then - files+=($file) - fi + # only indent existing files, this is necessary since if we rename or delete + # a file it is added to the committed files and we thus would try to indent a + # nonexisting file. + if [ ! -e $file ] + then + continue + fi + # We only indent .jl files + FILE_ENDING="${file##*.}" + if [ $FILE_ENDING = "jl" ] + then + files+=($file) + fi done julia utils/trixi-format-file.jl "${files[@]}" From 196f139069370bd1d450b1d3865170c447c2b74f Mon Sep 17 00:00:00 2001 From: Johannes Markert <10619309+jmark@users.noreply.github.com> Date: Tue, 20 Jun 2023 09:12:02 +0200 Subject: [PATCH 11/15] Update styleguide.md (#1538) Added missing apostrophe. --- docs/src/styleguide.md | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/docs/src/styleguide.md b/docs/src/styleguide.md index 9e6b6a8c265..c2562a0c651 100644 --- a/docs/src/styleguide.md +++ b/docs/src/styleguide.md @@ -55,7 +55,7 @@ julia -e 'using Pkg; Pkg.add("JuliaFormatter")' ``` You can then recursively format all Julia files in the Trixi.jl repo by executing ```shell -julia -e 'using JuliaFormatter; format(".") +julia -e 'using JuliaFormatter; format(".")' ``` from inside the Trixi.jl repository. For convenience, there is also a script you can directly run from your terminal shell, which will automatically install JuliaFormatter in a @@ -75,4 +75,4 @@ You can copy the `pre-commit`-script into `.git/hooks/pre-commit` and it will ch before each commit. If errors are found the commit is aborted and you can add the corrections via ```shell git add -p -``` \ No newline at end of file +``` From 3303ed8c0b8af262a16e96f2ad6dcd84034bbe7b Mon Sep 17 00:00:00 2001 From: Johannes Markert <10619309+jmark@users.noreply.github.com> Date: Tue, 20 Jun 2023 11:45:41 +0200 Subject: [PATCH 12/15] Update styleguide.md (#1540) * Update styleguide.md Updated the formatting command line fixing the issue https://github.com/trixi-framework/Trixi.jl/issues/1539 * Update styleguide.md Removed superfluous whitespace. --- docs/src/styleguide.md | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/docs/src/styleguide.md b/docs/src/styleguide.md index c2562a0c651..60e227204ca 100644 --- a/docs/src/styleguide.md +++ b/docs/src/styleguide.md @@ -53,9 +53,9 @@ of your PR), you need to install JuliaFormatter.jl first by running ```shell julia -e 'using Pkg; Pkg.add("JuliaFormatter")' ``` -You can then recursively format all Julia files in the Trixi.jl repo by executing +You can then recursively format the core Julia files in the Trixi.jl repo by executing ```shell -julia -e 'using JuliaFormatter; format(".")' +julia -e 'using JuliaFormatter; format(["benchmark", "ext", "src", "utils"])' ``` from inside the Trixi.jl repository. For convenience, there is also a script you can directly run from your terminal shell, which will automatically install JuliaFormatter in a From 414dafa310738414d43d16d0393ca13588b2c101 Mon Sep 17 00:00:00 2001 From: Joshua Lampert <51029046+JoshuaLampert@users.noreply.github.com> Date: Wed, 21 Jun 2023 11:31:35 +0200 Subject: [PATCH 13/15] Fix some typos in docs (#1541) * fix some typos in docs * fix copy mistake * update link --- docs/src/github-git.md | 4 ++-- src/auxiliary/auxiliary.jl | 10 +++++----- src/auxiliary/math.jl | 6 +++++- src/auxiliary/mpi.jl | 4 ++-- src/equations/acoustic_perturbation_2d.jl | 2 +- src/equations/compressible_euler_2d.jl | 2 +- src/equations/compressible_euler_3d.jl | 2 +- src/equations/compressible_euler_multicomponent_1d.jl | 8 ++++---- src/equations/compressible_euler_multicomponent_2d.jl | 8 ++++---- src/equations/compressible_navier_stokes_2d.jl | 8 ++++---- src/equations/compressible_navier_stokes_3d.jl | 4 ++-- src/equations/hyperbolic_diffusion_2d.jl | 2 +- src/equations/hyperbolic_diffusion_3d.jl | 2 +- src/equations/ideal_glm_mhd_3d.jl | 2 +- src/equations/shallow_water_two_layer_1d.jl | 2 +- src/equations/shallow_water_two_layer_2d.jl | 3 ++- src/meshes/structured_mesh.jl | 4 ++-- src/solvers/dgmulti/types.jl | 6 +++--- src/solvers/dgsem_tree/indicators.jl | 2 +- 19 files changed, 43 insertions(+), 38 deletions(-) diff --git a/docs/src/github-git.md b/docs/src/github-git.md index ad5991d87af..57b63073e79 100644 --- a/docs/src/github-git.md +++ b/docs/src/github-git.md @@ -112,7 +112,7 @@ branch, and the corresponding pull request will be updated automatically. Please note that a review has nothing to do with the lack of experience of the person developing changes: We try to review all code before it gets added to `main`, even from the most experienced developers. This is good practice and -helps to keep the error rate low while ensuring the the code is developed in a +helps to keep the error rate low while ensuring that the code is developed in a consistent fashion. Furthermore, do not take criticism of your code personally - we just try to keep Trixi.jl as accessible and easy to use for everyone. @@ -121,7 +121,7 @@ Once your branch is reviewed and declared ready for merging by the reviewer, make sure that all the latest changes have been pushed. Then, one of the developers will merge your PR. If you are one of the developers, you can also go to the pull request page on GitHub and and click on **Merge pull request**. -Voilá, you are done! Your branch will have been merged to +Voilà, you are done! Your branch will have been merged to `main` and the source branch will have been deleted in the GitHub repository (if you are not working in your own fork). diff --git a/src/auxiliary/auxiliary.jl b/src/auxiliary/auxiliary.jl index 115d055c0ca..1f7d30d6aa8 100644 --- a/src/auxiliary/auxiliary.jl +++ b/src/auxiliary/auxiliary.jl @@ -132,7 +132,7 @@ end default_example() Return the path to an example elixir that can be used to quickly see Trixi.jl in action on a -[`TreeMesh`]@(ref). See also [`examples_dir`](@ref) and [`get_examples`](@ref). +[`TreeMesh`](@ref). See also [`examples_dir`](@ref) and [`get_examples`](@ref). """ function default_example() joinpath(examples_dir(), "tree_2d_dgsem", "elixir_advection_basic.jl") @@ -142,7 +142,7 @@ end default_example_unstructured() Return the path to an example elixir that can be used to quickly see Trixi.jl in action on an -[`UnstructuredMesh2D`]@(ref). This simulation is run on the example curved, unstructured mesh +[`UnstructuredMesh2D`](@ref). This simulation is run on the example curved, unstructured mesh given in the Trixi.jl documentation regarding unstructured meshes. """ function default_example_unstructured() @@ -155,7 +155,7 @@ end Return the default options for OrdinaryDiffEq's `solve`. Pass `ode_default_options()...` to `solve` to only return the solution at the final time and enable **MPI aware** error-based step size control, whenever MPI is used. -For example, use `solve(ode, alg; ode_default_options()...)` +For example, use `solve(ode, alg; ode_default_options()...)`. """ function ode_default_options() if mpi_isparallel() @@ -213,8 +213,8 @@ might be provided by other packages such as [Polyester.jl](https://github.com/Ju This macro does not necessarily work for general `for` loops. For example, it does not necessarily support general iterables such as `eachline(filename)`. -Some discussion can be found at https://discourse.julialang.org/t/overhead-of-threads-threads/53964 -and https://discourse.julialang.org/t/threads-threads-with-one-thread-how-to-remove-the-overhead/58435. +Some discussion can be found at [https://discourse.julialang.org/t/overhead-of-threads-threads/53964](https://discourse.julialang.org/t/overhead-of-threads-threads/53964) +and [https://discourse.julialang.org/t/threads-threads-with-one-thread-how-to-remove-the-overhead/58435](https://discourse.julialang.org/t/threads-threads-with-one-thread-how-to-remove-the-overhead/58435). """ macro threaded(expr) # Use `esc(quote ... end)` for nested macro calls as suggested in diff --git a/src/auxiliary/math.jl b/src/auxiliary/math.jl index 27c1bed5ca4..4ecf7dd3fcc 100644 --- a/src/auxiliary/math.jl +++ b/src/auxiliary/math.jl @@ -51,7 +51,7 @@ Given ε = 1.0e-4, we use the following algorithm. - Agner Fog. Lists of instruction latencies, throughputs and micro-operation breakdowns for Intel, AMD, and VIA CPUs. - https://www.agner.org/optimize/instruction_tables.pdf + [https://www.agner.org/optimize/instruction_tables.pdf](https://www.agner.org/optimize/instruction_tables.pdf) """ @inline function ln_mean(x, y) epsilon_f2 = 1.0e-4 @@ -166,8 +166,10 @@ checks necessary in the presence of `NaN`s (or signed zeros). # Examples +```jldoctest julia> max(2, 5, 1) 5 +``` """ @inline max(args...) = @fastmath max(args...) @@ -183,8 +185,10 @@ checks necessary in the presence of `NaN`s (or signed zeros). # Examples +```jldoctest julia> min(2, 5, 1) 1 +``` """ @inline min(args...) = @fastmath min(args...) diff --git a/src/auxiliary/mpi.jl b/src/auxiliary/mpi.jl index 2c485b4832c..c85c23670b0 100644 --- a/src/auxiliary/mpi.jl +++ b/src/auxiliary/mpi.jl @@ -72,7 +72,7 @@ You must pass this function as a keyword argument to OrdinaryDiffEq.jl's `solve` when using error-based step size control with MPI parallel execution of Trixi.jl. -See the "Advanced Adaptive Stepsize Control" section of the [documentation](https://docs.sciml.ai/DiffEqDocs/stable/basics/common_solver_opts/) +See the "Advanced Adaptive Stepsize Control" section of the [documentation](https://docs.sciml.ai/DiffEqDocs/stable/basics/common_solver_opts/). """ ode_norm(u::Number, t) = @fastmath abs(u) function ode_norm(u::AbstractArray, t) @@ -125,6 +125,6 @@ You should pass this function as a keyword argument to OrdinaryDiffEq.jl's `solve` when using error-based step size control with MPI parallel execution of Trixi.jl. -See the "Miscellaneous" section of the [documentation](https://docs.sciml.ai/DiffEqDocs/stable/basics/common_solver_opts/) +See the "Miscellaneous" section of the [documentation](https://docs.sciml.ai/DiffEqDocs/stable/basics/common_solver_opts/). """ ode_unstable_check(dt, u, semi, t) = isnan(dt) diff --git a/src/equations/acoustic_perturbation_2d.jl b/src/equations/acoustic_perturbation_2d.jl index 786630a14c7..f4ce770e1e9 100644 --- a/src/equations/acoustic_perturbation_2d.jl +++ b/src/equations/acoustic_perturbation_2d.jl @@ -145,7 +145,7 @@ function initial_condition_convergence_test(x, t, end """ - source_terms_convergence_test(u, x, t, equations::AcousticPerturbationEquations2D) + source_terms_convergence_test(u, x, t, equations::AcousticPerturbationEquations2D) Source terms used for convergence tests in combination with [`initial_condition_convergence_test`](@ref). diff --git a/src/equations/compressible_euler_2d.jl b/src/equations/compressible_euler_2d.jl index 66e3c7bff84..89f04ef1e05 100644 --- a/src/equations/compressible_euler_2d.jl +++ b/src/equations/compressible_euler_2d.jl @@ -31,7 +31,7 @@ The compressible Euler equations ``` for an ideal gas with ratio of specific heats `gamma` in two space dimensions. -Here, ``\rho`` is the density, ``v_1``,`v_2` the velocities, ``e`` the specific total energy **rather than** specific internal energy, and +Here, ``\rho`` is the density, ``v_1``, ``v_2`` the velocities, ``e`` the specific total energy **rather than** specific internal energy, and ```math p = (\gamma - 1) \left( \rho e - \frac{1}{2} \rho (v_1^2+v_2^2) \right) ``` diff --git a/src/equations/compressible_euler_3d.jl b/src/equations/compressible_euler_3d.jl index c16a454b176..cd081cfc42a 100644 --- a/src/equations/compressible_euler_3d.jl +++ b/src/equations/compressible_euler_3d.jl @@ -36,7 +36,7 @@ The compressible Euler equations ``` for an ideal gas with ratio of specific heats `gamma` in three space dimensions. -Here, ``\rho`` is the density, ``v_1``,`v_2`, `v_3` the velocities, ``e`` the specific total energy **rather than** specific internal energy, and +Here, ``\rho`` is the density, ``v_1``, ``v_2``, ``v_3`` the velocities, ``e`` the specific total energy **rather than** specific internal energy, and ```math p = (\gamma - 1) \left( \rho e - \frac{1}{2} \rho (v_1^2+v_2^2+v_3^2) \right) ``` diff --git a/src/equations/compressible_euler_multicomponent_1d.jl b/src/equations/compressible_euler_multicomponent_1d.jl index 4a50d60471a..23ac222b976 100644 --- a/src/equations/compressible_euler_multicomponent_1d.jl +++ b/src/equations/compressible_euler_multicomponent_1d.jl @@ -44,8 +44,8 @@ specific heat capacity at constant volume of component ``i``. In case of more than one component, the specific heat ratios `gammas` and the gas constants `gas_constants` should be passed as tuples, e.g., `gammas=(1.4, 1.667)`. -The remaining variables like the specific heats at constant volume 'cv' or the specific heats at -constant pressure 'cp' are then calculated considering a calorically perfect gas. +The remaining variables like the specific heats at constant volume `cv` or the specific heats at +constant pressure `cp` are then calculated considering a calorically perfect gas. """ struct CompressibleEulerMulticomponentEquations1D{NVARS, NCOMP, RealT <: Real} <: AbstractCompressibleEulerMulticomponentEquations{1, NVARS, NCOMP} @@ -247,8 +247,8 @@ end Entropy conserving two-point flux by - Ayoub Gouasmi, Karthik Duraisamy (2020) - "Formulation of Entropy-Stable schemes for the multicomponent compressible Euler equations"" - arXiv:1904.00972v3 [math.NA] 4 Feb 2020 + "Formulation of Entropy-Stable schemes for the multicomponent compressible Euler equations" + [arXiv:1904.00972v3](https://arxiv.org/abs/1904.00972) [math.NA] 4 Feb 2020 """ @inline function flux_chandrashekar(u_ll, u_rr, orientation::Integer, equations::CompressibleEulerMulticomponentEquations1D) diff --git a/src/equations/compressible_euler_multicomponent_2d.jl b/src/equations/compressible_euler_multicomponent_2d.jl index 5a015777cb1..7b437f4a1b4 100644 --- a/src/equations/compressible_euler_multicomponent_2d.jl +++ b/src/equations/compressible_euler_multicomponent_2d.jl @@ -48,8 +48,8 @@ specific heat capacity at constant volume of component ``i``. In case of more than one component, the specific heat ratios `gammas` and the gas constants `gas_constants` in [kJ/(kg*K)] should be passed as tuples, e.g., `gammas=(1.4, 1.667)`. -The remaining variables like the specific heats at constant volume 'cv' or the specific heats at -constant pressure 'cp' are then calculated considering a calorically perfect gas. +The remaining variables like the specific heats at constant volume `cv` or the specific heats at +constant pressure `cp` are then calculated considering a calorically perfect gas. """ struct CompressibleEulerMulticomponentEquations2D{NVARS, NCOMP, RealT <: Real} <: AbstractCompressibleEulerMulticomponentEquations{2, NVARS, NCOMP} @@ -275,8 +275,8 @@ end Adaption of the entropy conserving two-point flux by - Ayoub Gouasmi, Karthik Duraisamy (2020) - "Formulation of Entropy-Stable schemes for the multicomponent compressible Euler equations"" - arXiv:1904.00972v3 [math.NA] 4 Feb 2020 + "Formulation of Entropy-Stable schemes for the multicomponent compressible Euler equations" + [arXiv:1904.00972v3](https://arxiv.org/abs/1904.00972) [math.NA] 4 Feb 2020 """ @inline function flux_chandrashekar(u_ll, u_rr, orientation::Integer, equations::CompressibleEulerMulticomponentEquations2D) diff --git a/src/equations/compressible_navier_stokes_2d.jl b/src/equations/compressible_navier_stokes_2d.jl index 33badba15d9..9b06e0b5abf 100644 --- a/src/equations/compressible_navier_stokes_2d.jl +++ b/src/equations/compressible_navier_stokes_2d.jl @@ -73,8 +73,8 @@ where w_2 = \frac{\rho v_1}{p},\, w_3 = \frac{\rho v_2}{p},\, w_4 = -\frac{\rho}{p} ``` -#!!! warning "Experimental code" -# This code is experimental and may be changed or removed in any future release. +!!! warning "Experimental code" + This code is experimental and may be changed or removed in any future release. """ struct CompressibleNavierStokesDiffusion2D{GradientVariables, RealT <: Real, E <: AbstractCompressibleEulerEquations{2}} <: @@ -94,8 +94,8 @@ struct CompressibleNavierStokesDiffusion2D{GradientVariables, RealT <: Real, end """ -#!!! warning "Experimental code" -# This code is experimental and may be changed or removed in any future release. +!!! warning "Experimental code" + This code is experimental and may be changed or removed in any future release. `GradientVariablesPrimitive` and `GradientVariablesEntropy` are gradient variable type parameters for `CompressibleNavierStokesDiffusion2D`. By default, the gradient variables are set to be diff --git a/src/equations/compressible_navier_stokes_3d.jl b/src/equations/compressible_navier_stokes_3d.jl index 8930489295d..0b770dff1ca 100644 --- a/src/equations/compressible_navier_stokes_3d.jl +++ b/src/equations/compressible_navier_stokes_3d.jl @@ -73,8 +73,8 @@ where w_2 = \frac{\rho v_1}{p},\, w_3 = \frac{\rho v_2}{p},\, w_4 = \frac{\rho v_3}{p},\, w_5 = -\frac{\rho}{p} ``` -#!!! warning "Experimental code" -# This code is experimental and may be changed or removed in any future release. +!!! warning "Experimental code" + This code is experimental and may be changed or removed in any future release. """ struct CompressibleNavierStokesDiffusion3D{GradientVariables, RealT <: Real, E <: AbstractCompressibleEulerEquations{3}} <: diff --git a/src/equations/hyperbolic_diffusion_2d.jl b/src/equations/hyperbolic_diffusion_2d.jl index 25536a060f8..511d1b8935d 100644 --- a/src/equations/hyperbolic_diffusion_2d.jl +++ b/src/equations/hyperbolic_diffusion_2d.jl @@ -10,7 +10,7 @@ The linear hyperbolic diffusion equations in two space dimensions. A description of this system can be found in Sec. 2.5 of the book "I Do Like CFD, Too: Vol 1". -The book is freely available at http://www.cfdbooks.com/ and further analysis can be found in +The book is freely available at [http://www.cfdbooks.com/](http://www.cfdbooks.com/) and further analysis can be found in the paper by Nishikawa [DOI: 10.1016/j.jcp.2007.07.029](https://doi.org/10.1016/j.jcp.2007.07.029) """ struct HyperbolicDiffusionEquations2D{RealT <: Real} <: diff --git a/src/equations/hyperbolic_diffusion_3d.jl b/src/equations/hyperbolic_diffusion_3d.jl index bf6a00140d4..ed807511b67 100644 --- a/src/equations/hyperbolic_diffusion_3d.jl +++ b/src/equations/hyperbolic_diffusion_3d.jl @@ -10,7 +10,7 @@ The linear hyperbolic diffusion equations in three space dimensions. A description of this system can be found in Sec. 2.5 of the book "I Do Like CFD, Too: Vol 1". -The book is freely available at http://www.cfdbooks.com/ and further analysis can be found in +The book is freely available at [http://www.cfdbooks.com/](http://www.cfdbooks.com/) and further analysis can be found in the paper by Nishikawa [DOI: 10.1016/j.jcp.2007.07.029](https://doi.org/10.1016/j.jcp.2007.07.029) """ struct HyperbolicDiffusionEquations3D{RealT <: Real} <: diff --git a/src/equations/ideal_glm_mhd_3d.jl b/src/equations/ideal_glm_mhd_3d.jl index 401fcd2daf1..2e149d2849f 100644 --- a/src/equations/ideal_glm_mhd_3d.jl +++ b/src/equations/ideal_glm_mhd_3d.jl @@ -41,7 +41,7 @@ end # Set initial conditions at physical location `x` for time `t` """ -initial_condition_constant(x, t, equations::IdealGlmMhdEquations3D) + initial_condition_constant(x, t, equations::IdealGlmMhdEquations3D) A constant initial condition to test free-stream preservation. """ diff --git a/src/equations/shallow_water_two_layer_1d.jl b/src/equations/shallow_water_two_layer_1d.jl index 02899171509..e126eec7c25 100644 --- a/src/equations/shallow_water_two_layer_1d.jl +++ b/src/equations/shallow_water_two_layer_1d.jl @@ -392,7 +392,7 @@ end equations::ShallowWaterTwoLayerEquations1D) Entropy stable surface flux for the two-layer shallow water equations. Uses the entropy -conservative flux_fjordholm_etal and adds a Lax-Friedrichs type dissipation dependent on the jump +conservative [`flux_fjordholm_etal`](@ref) and adds a Lax-Friedrichs type dissipation dependent on the jump of entropy variables. Further details are available in the paper: diff --git a/src/equations/shallow_water_two_layer_2d.jl b/src/equations/shallow_water_two_layer_2d.jl index b5e52d636e4..a54831c711f 100644 --- a/src/equations/shallow_water_two_layer_2d.jl +++ b/src/equations/shallow_water_two_layer_2d.jl @@ -695,8 +695,9 @@ end """ flux_es_fjordholm_etal(u_ll, u_rr, orientation_or_normal_direction, equations::ShallowWaterTwoLayerEquations1D) + Entropy stable surface flux for the two-layer shallow water equations. Uses the entropy conservative -flux_fjordholm_etal and adds a Lax-Friedrichs type dissipation dependent on the jump of entropy +[`flux_fjordholm_etal`](@ref) and adds a Lax-Friedrichs type dissipation dependent on the jump of entropy variables. Further details are available in the paper: diff --git a/src/meshes/structured_mesh.jl b/src/meshes/structured_mesh.jl index 5872681933a..df067db833d 100644 --- a/src/meshes/structured_mesh.jl +++ b/src/meshes/structured_mesh.jl @@ -33,7 +33,7 @@ Create a StructuredMesh of the given size and shape that uses `RealT` as coordin the reference mesh to the physical domain. If no `mapping_as_string` is defined, this function must be defined with the name `mapping` to allow for restarts. - This will be changed in the future, see https://github.com/trixi-framework/Trixi.jl/issues/541. + This will be changed in the future, see [https://github.com/trixi-framework/Trixi.jl/issues/541](https://github.com/trixi-framework/Trixi.jl/issues/541). - `RealT::Type`: the type that should be used for coordinates. - `periodicity`: either a `Bool` deciding if all of the boundaries are periodic or an `NTuple{NDIMS, Bool}` deciding for each dimension if the boundaries in this dimension are periodic. @@ -41,7 +41,7 @@ Create a StructuredMesh of the given size and shape that uses `RealT` as coordin - `mapping_as_string::String`: the code that defines the `mapping`. If `CodeTracking` can't find the function definition, it can be passed directly here. The code string must define the mapping function with the name `mapping`. - This will be changed in the future, see https://github.com/trixi-framework/Trixi.jl/issues/541. + This will be changed in the future, see [https://github.com/trixi-framework/Trixi.jl/issues/541](https://github.com/trixi-framework/Trixi.jl/issues/541). """ function StructuredMesh(cells_per_dimension, mapping; RealT = Float64, periodicity = true, unsaved_changes = true, diff --git a/src/solvers/dgmulti/types.jl b/src/solvers/dgmulti/types.jl index f1f7b158dec..fe6510856b0 100644 --- a/src/solvers/dgmulti/types.jl +++ b/src/solvers/dgmulti/types.jl @@ -180,9 +180,9 @@ GeometricTermsType(mesh_type::Curved, element_type::AbstractElemShape) = NonAffi # other potential mesh types to add later: Polynomial{polydeg_geo}? """ - DGMultiMesh(dg::DGMulti{NDIMS}, vertex_coordinates, EToV; - is_on_boundary=nothing, - periodicity=ntuple(_->false, NDIMS)) where {NDIMS} + DGMultiMesh(dg::DGMulti{NDIMS}, vertex_coordinates, EToV; + is_on_boundary=nothing, + periodicity=ntuple(_->false, NDIMS)) where {NDIMS} - `dg::DGMulti` contains information associated with to the reference element (e.g., quadrature, basis evaluation, differentiation, etc). diff --git a/src/solvers/dgsem_tree/indicators.jl b/src/solvers/dgsem_tree/indicators.jl index 2eb0af87148..b8f8a796f2b 100644 --- a/src/solvers/dgsem_tree/indicators.jl +++ b/src/solvers/dgsem_tree/indicators.jl @@ -159,7 +159,7 @@ and `basis` if this indicator should be used for shock capturing. - Löhner (1987) "An adaptive finite element scheme for transient problems in CFD" [doi: 10.1016/0045-7825(87)90098-3](https://doi.org/10.1016/0045-7825(87)90098-3) -- http://flash.uchicago.edu/site/flashcode/user_support/flash4_ug_4p62/node59.html#SECTION05163100000000000000 +- [https://flash.rochester.edu/site/flashcode/user_support/flash4_ug_4p62/node59.html#SECTION05163100000000000000](https://flash.rochester.edu/site/flashcode/user_support/flash4_ug_4p62/node59.html#SECTION05163100000000000000) """ struct IndicatorLöhner{RealT <: Real, Variable, Cache} <: AbstractIndicator f_wave::RealT # TODO: Taal documentation From 054a917a09127570dabc458f1350550f2ddb6a09 Mon Sep 17 00:00:00 2001 From: Andrew Winters Date: Wed, 21 Jun 2023 17:56:07 +0200 Subject: [PATCH 14/15] update links to Flash manual (#1544) Co-authored-by: Hendrik Ranocha --- examples/p4est_2d_dgsem/elixir_euler_sedov.jl | 4 ++-- examples/p4est_3d_dgsem/elixir_euler_sedov.jl | 4 ++-- .../elixir_eulergravity_jeans_instability.jl | 6 +++--- .../elixir_eulergravity_sedov_blast_wave.jl | 10 +++++----- examples/structured_1d_dgsem/elixir_euler_sedov.jl | 8 ++++---- examples/structured_2d_dgsem/elixir_euler_sedov.jl | 10 +++++----- examples/structured_3d_dgsem/elixir_euler_sedov.jl | 10 +++++----- examples/tree_1d_dgsem/elixir_euler_positivity.jl | 4 ++-- .../tree_1d_dgsem/elixir_euler_sedov_blast_wave.jl | 4 ++-- .../elixir_euler_sedov_blast_wave_pure_fv.jl | 4 ++-- examples/tree_2d_dgsem/elixir_euler_positivity.jl | 4 ++-- .../tree_2d_dgsem/elixir_euler_sedov_blast_wave.jl | 4 ++-- ...er_sedov_blast_wave_neuralnetwork_perssonperaire.jl | 4 ++-- .../tree_3d_dgsem/elixir_euler_sedov_blast_wave.jl | 6 +++--- examples/unstructured_2d_dgsem/elixir_euler_sedov.jl | 4 ++-- 15 files changed, 43 insertions(+), 43 deletions(-) diff --git a/examples/p4est_2d_dgsem/elixir_euler_sedov.jl b/examples/p4est_2d_dgsem/elixir_euler_sedov.jl index 9f5247e8c4d..d5d8e0c78bf 100644 --- a/examples/p4est_2d_dgsem/elixir_euler_sedov.jl +++ b/examples/p4est_2d_dgsem/elixir_euler_sedov.jl @@ -11,7 +11,7 @@ equations = CompressibleEulerEquations2D(1.4) initial_condition_sedov_blast_wave(x, t, equations::CompressibleEulerEquations2D) The Sedov blast wave setup based on Flash -- http://flash.uchicago.edu/site/flashcode/user_support/flash_ug_devel/node184.html#SECTION010114000000000000000 +- https://flash.rochester.edu/site/flashcode/user_support/flash_ug_devel/node187.html#SECTION010114000000000000000 """ function initial_condition_sedov_blast_wave(x, t, equations::CompressibleEulerEquations2D) # Set up polar coordinates @@ -20,7 +20,7 @@ function initial_condition_sedov_blast_wave(x, t, equations::CompressibleEulerEq y_norm = x[2] - inicenter[2] r = sqrt(x_norm^2 + y_norm^2) - # Setup based on http://flash.uchicago.edu/site/flashcode/user_support/flash_ug_devel/node184.html#SECTION010114000000000000000 + # Setup based on https://flash.rochester.edu/site/flashcode/user_support/flash_ug_devel/node187.html#SECTION010114000000000000000 r0 = 0.21875 # = 3.5 * smallest dx (for domain length=4 and max-ref=6) E = 1.0 p0_inner = 3 * (equations.gamma - 1) * E / (3 * pi * r0^2) diff --git a/examples/p4est_3d_dgsem/elixir_euler_sedov.jl b/examples/p4est_3d_dgsem/elixir_euler_sedov.jl index 00da4132851..6fa285b5565 100644 --- a/examples/p4est_3d_dgsem/elixir_euler_sedov.jl +++ b/examples/p4est_3d_dgsem/elixir_euler_sedov.jl @@ -11,7 +11,7 @@ equations = CompressibleEulerEquations3D(1.4) initial_condition_medium_sedov_blast_wave(x, t, equations::CompressibleEulerEquations3D) The Sedov blast wave setup based on Flash -- http://flash.uchicago.edu/site/flashcode/user_support/flash_ug_devel/node184.html#SECTION010114000000000000000 +- https://flash.rochester.edu/site/flashcode/user_support/flash_ug_devel/node187.html#SECTION010114000000000000000 with smaller strength of the initial discontinuity. """ function initial_condition_medium_sedov_blast_wave(x, t, equations::CompressibleEulerEquations3D) @@ -22,7 +22,7 @@ function initial_condition_medium_sedov_blast_wave(x, t, equations::Compressible z_norm = x[3] - inicenter[3] r = sqrt(x_norm^2 + y_norm^2 + z_norm^2) - # Setup based on http://flash.uchicago.edu/site/flashcode/user_support/flash_ug_devel/node184.html#SECTION010114000000000000000 + # Setup based on https://flash.rochester.edu/site/flashcode/user_support/flash_ug_devel/node187.html#SECTION010114000000000000000 r0 = 0.21875 # = 3.5 * smallest dx (for domain length=4 and max-ref=6) E = 1.0 p0_inner = 3 * (equations.gamma - 1) * E / (4 * pi * r0^2) diff --git a/examples/paper_self_gravitating_gas_dynamics/elixir_eulergravity_jeans_instability.jl b/examples/paper_self_gravitating_gas_dynamics/elixir_eulergravity_jeans_instability.jl index 1774e39513d..fb445616cd4 100644 --- a/examples/paper_self_gravitating_gas_dynamics/elixir_eulergravity_jeans_instability.jl +++ b/examples/paper_self_gravitating_gas_dynamics/elixir_eulergravity_jeans_instability.jl @@ -15,7 +15,7 @@ The classical Jeans instability taken from - Dominik Derigs, Andrew R. Winters, Gregor J. Gassner, Stefanie Walch (2016) A Novel High-Order, Entropy Stable, 3D AMR MHD Solver with Guaranteed Positive Pressure [arXiv: 1605.03572](https://arxiv.org/abs/1605.03572) -- Flash manual https://flash.uchicago.edu/site/flashcode/user_support/flash_ug_devel.pdf +- Flash manual https://flash.rochester.edu/site/flashcode/user_support/flash_ug_devel/node189.html#SECTION010131000000000000000 in CGS (centimeter, gram, second) units. """ function initial_condition_jeans_instability(x, t, @@ -32,7 +32,7 @@ function initial_condition_jeans_instability(x, t, pres0 = 1.5e7 # dyn/cm^2 delta0 = 1e-3 # set wave vector values for perturbation (units 1/cm) - # see FLASH manual: https://flash.uchicago.edu/site/flashcode/user_support/flash_ug_devel.pdf + # see FLASH manual: https://flash.rochester.edu/site/flashcode/user_support/flash_ug_devel/node189.html#SECTION010131000000000000000 kx = 2.0*pi/0.5 # 2π/λ_x, λ_x = 0.5 ky = 0.0 # 2π/λ_y, λ_y = 1e10 k_dot_x = kx*x[1] + ky*x[2] @@ -49,7 +49,7 @@ function initial_condition_jeans_instability(x, t, equations::HyperbolicDiffusionEquations2D) # gravity equation: -Δϕ = -4πGρ # Constants taken from the FLASH manual - # https://flash.uchicago.edu/site/flashcode/user_support/flash_ug_devel.pdf + # https://flash.rochester.edu/site/flashcode/user_support/flash_ug_devel/node189.html#SECTION010131000000000000000 rho0 = 1.5e7 delta0 = 1e-3 diff --git a/examples/paper_self_gravitating_gas_dynamics/elixir_eulergravity_sedov_blast_wave.jl b/examples/paper_self_gravitating_gas_dynamics/elixir_eulergravity_sedov_blast_wave.jl index f7bb5bbb01c..8933224a2c7 100644 --- a/examples/paper_self_gravitating_gas_dynamics/elixir_eulergravity_sedov_blast_wave.jl +++ b/examples/paper_self_gravitating_gas_dynamics/elixir_eulergravity_sedov_blast_wave.jl @@ -15,14 +15,14 @@ Adaptation of the Sedov blast wave with self-gravity taken from A purely hyperbolic discontinuous Galerkin approach for self-gravitating gas dynamics [arXiv: 2008.10593](https://arxiv.org/abs/2008.10593) based on -- http://flash.uchicago.edu/site/flashcode/user_support/flash4_ug_4p62/node184.html#SECTION010114000000000000000 +- https://flash.rochester.edu/site/flashcode/user_support/flash_ug_devel/node187.html#SECTION010114100000000000000 Should be used together with [`boundary_condition_sedov_self_gravity`](@ref). """ function initial_condition_sedov_self_gravity(x, t, equations::CompressibleEulerEquations2D) # Set up polar coordinates r = sqrt(x[1]^2 + x[2]^2) - # Setup based on http://flash.uchicago.edu/site/flashcode/user_support/flash4_ug_4p62/node184.html#SECTION010114000000000000000 + # Setup based on https://flash.rochester.edu/site/flashcode/user_support/flash_ug_devel/node187.html#SECTION010114100000000000000 r0 = 0.125 # = 4.0 * smallest dx (for domain length=8 and max-ref=8) E = 1.0 p_inner = (equations.gamma - 1) * E / (pi * r0^2) @@ -59,7 +59,7 @@ Adaptation of the Sedov blast wave with self-gravity taken from A purely hyperbolic discontinuous Galerkin approach for self-gravitating gas dynamics [arXiv: 2008.10593](https://arxiv.org/abs/2008.10593) based on -- http://flash.uchicago.edu/site/flashcode/user_support/flash4_ug_4p62/node184.html#SECTION010114000000000000000 +- https://flash.rochester.edu/site/flashcode/user_support/flash_ug_devel/node187.html#SECTION010114100000000000000 Should be used together with [`initial_condition_sedov_self_gravity`](@ref). """ function boundary_condition_sedov_self_gravity(u_inner, orientation, direction, x, t, @@ -122,7 +122,7 @@ Adaptation of the Sedov blast wave with self-gravity taken from A purely hyperbolic discontinuous Galerkin approach for self-gravitating gas dynamics [arXiv: 2008.10593](https://arxiv.org/abs/2008.10593) based on -- http://flash.uchicago.edu/site/flashcode/user_support/flash4_ug_4p62/node184.html#SECTION010114000000000000000 +- https://flash.rochester.edu/site/flashcode/user_support/flash_ug_devel/node187.html#SECTION010114100000000000000 Should be used together with [`boundary_condition_sedov_self_gravity`](@ref). """ function initial_condition_sedov_self_gravity(x, t, equations::HyperbolicDiffusionEquations2D) @@ -143,7 +143,7 @@ Adaptation of the Sedov blast wave with self-gravity taken from A purely hyperbolic discontinuous Galerkin approach for self-gravitating gas dynamics [arXiv: 2008.10593](https://arxiv.org/abs/2008.10593) based on -- http://flash.uchicago.edu/site/flashcode/user_support/flash4_ug_4p62/node184.html#SECTION010114000000000000000 +- https://flash.rochester.edu/site/flashcode/user_support/flash_ug_devel/node187.html#SECTION010114100000000000000 Should be used together with [`initial_condition_sedov_self_gravity`](@ref). """ function boundary_condition_sedov_self_gravity(u_inner, orientation, direction, x, t, diff --git a/examples/structured_1d_dgsem/elixir_euler_sedov.jl b/examples/structured_1d_dgsem/elixir_euler_sedov.jl index ee950b3aaaa..9d7be21a5c1 100644 --- a/examples/structured_1d_dgsem/elixir_euler_sedov.jl +++ b/examples/structured_1d_dgsem/elixir_euler_sedov.jl @@ -11,7 +11,7 @@ equations = CompressibleEulerEquations1D(1.4) initial_condition_sedov_blast_wave(x, t, equations::CompressibleEulerEquations1D) The Sedov blast wave setup based on Flash -- http://flash.uchicago.edu/site/flashcode/user_support/flash_ug_devel/node184.html#SECTION010114000000000000000 +- https://flash.rochester.edu/site/flashcode/user_support/flash_ug_devel/node187.html#SECTION010114000000000000000 """ function initial_condition_sedov_blast_wave(x, t, equations::CompressibleEulerEquations1D) # Set up polar coordinates @@ -19,7 +19,7 @@ function initial_condition_sedov_blast_wave(x, t, equations::CompressibleEulerEq x_norm = x[1] - inicenter[1] r = abs(x_norm) - # Setup based on http://flash.uchicago.edu/site/flashcode/user_support/flash_ug_devel/node184.html#SECTION010114000000000000000 + # Setup based on https://flash.rochester.edu/site/flashcode/user_support/flash_ug_devel/node187.html#SECTION010114000000000000000 r0 = 0.21875 # = 3.5 * smallest dx (for domain length=4 and max-ref=6) # r0 = 0.5 # = more reasonable setup E = 1.0 @@ -78,8 +78,8 @@ save_solution = SaveSolutionCallback(interval=100, stepsize_callback = StepsizeCallback(cfl=0.5) -callbacks = CallbackSet(summary_callback, - analysis_callback, +callbacks = CallbackSet(summary_callback, + analysis_callback, alive_callback, save_solution, stepsize_callback) diff --git a/examples/structured_2d_dgsem/elixir_euler_sedov.jl b/examples/structured_2d_dgsem/elixir_euler_sedov.jl index ed1bfab3be2..efc3b6627c0 100644 --- a/examples/structured_2d_dgsem/elixir_euler_sedov.jl +++ b/examples/structured_2d_dgsem/elixir_euler_sedov.jl @@ -11,7 +11,7 @@ equations = CompressibleEulerEquations2D(1.4) initial_condition_sedov_blast_wave(x, t, equations::CompressibleEulerEquations2D) The Sedov blast wave setup based on Flash -- http://flash.uchicago.edu/site/flashcode/user_support/flash_ug_devel/node184.html#SECTION010114000000000000000 +- https://flash.rochester.edu/site/flashcode/user_support/flash_ug_devel/node187.html#SECTION010114000000000000000 """ function initial_condition_sedov_blast_wave(x, t, equations::CompressibleEulerEquations2D) # Set up polar coordinates @@ -20,7 +20,7 @@ function initial_condition_sedov_blast_wave(x, t, equations::CompressibleEulerEq y_norm = x[2] - inicenter[2] r = sqrt(x_norm^2 + y_norm^2) - # Setup based on http://flash.uchicago.edu/site/flashcode/user_support/flash_ug_devel/node184.html#SECTION010114000000000000000 + # Setup based on https://flash.rochester.edu/site/flashcode/user_support/flash_ug_devel/node187.html#SECTION010114000000000000000 r0 = 0.21875 # = 3.5 * smallest dx (for domain length=4 and max-ref=6) E = 1.0 p0_inner = 3 * (equations.gamma - 1) * E / (3 * pi * r0^2) @@ -59,12 +59,12 @@ function mapping(xi, eta) y = eta + 0.125 * (cos(1.5 * pi * xi) * cos(0.5 * pi * eta)) x = xi + 0.125 * (cos(0.5 * pi * xi) * cos(2 * pi * y)) - + return SVector(x, y) end - + cells_per_dimension = (16, 16) - + mesh = StructuredMesh(cells_per_dimension, mapping, periodicity=true) # create the semidiscretization diff --git a/examples/structured_3d_dgsem/elixir_euler_sedov.jl b/examples/structured_3d_dgsem/elixir_euler_sedov.jl index 8f428495b4f..e0595437c99 100644 --- a/examples/structured_3d_dgsem/elixir_euler_sedov.jl +++ b/examples/structured_3d_dgsem/elixir_euler_sedov.jl @@ -11,7 +11,7 @@ equations = CompressibleEulerEquations3D(1.4) initial_condition_medium_sedov_blast_wave(x, t, equations::CompressibleEulerEquations3D) The Sedov blast wave setup based on Flash -- http://flash.uchicago.edu/site/flashcode/user_support/flash_ug_devel/node184.html#SECTION010114000000000000000 +- https://flash.rochester.edu/site/flashcode/user_support/flash_ug_devel/node187.html#SECTION010114000000000000000 with smaller strength of the initial discontinuity. """ function initial_condition_medium_sedov_blast_wave(x, t, equations::CompressibleEulerEquations3D) @@ -22,11 +22,11 @@ function initial_condition_medium_sedov_blast_wave(x, t, equations::Compressible z_norm = x[3] - inicenter[3] r = sqrt(x_norm^2 + y_norm^2 + z_norm^2) - # Setup based on http://flash.uchicago.edu/site/flashcode/user_support/flash_ug_devel/node184.html#SECTION010114000000000000000 + # Setup based on https://flash.rochester.edu/site/flashcode/user_support/flash_ug_devel/node187.html#SECTION010114000000000000000 r0 = 0.21875 # = 3.5 * smallest dx (for domain length=4 and max-ref=6) E = 1.0 p0_inner = 3 * (equations.gamma - 1) * E / (4 * pi * r0^2) - p0_outer = 1.0e-3 + p0_outer = 1.0e-3 # Calculate primitive variables rho = 1.0 @@ -52,8 +52,8 @@ indicator_sc = IndicatorHennemannGassner(equations, basis, volume_integral = VolumeIntegralShockCapturingHG(indicator_sc; volume_flux_dg=volume_flux, volume_flux_fv=surface_flux) - -solver = DGSEM(polydeg=polydeg, surface_flux=surface_flux, volume_integral=volume_integral) + +solver = DGSEM(polydeg=polydeg, surface_flux=surface_flux, volume_integral=volume_integral) # Mapping as described in https://arxiv.org/abs/2012.12040 function mapping(xi, eta, zeta) diff --git a/examples/tree_1d_dgsem/elixir_euler_positivity.jl b/examples/tree_1d_dgsem/elixir_euler_positivity.jl index 7942937151a..966661e8894 100644 --- a/examples/tree_1d_dgsem/elixir_euler_positivity.jl +++ b/examples/tree_1d_dgsem/elixir_euler_positivity.jl @@ -11,7 +11,7 @@ equations = CompressibleEulerEquations1D(1.4) initial_condition_sedov_blast_wave(x, t, equations::CompressibleEulerEquations1D) The Sedov blast wave setup based on Flash -- http://flash.uchicago.edu/site/flashcode/user_support/flash_ug_devel/node184.html#SECTION010114000000000000000 +- https://flash.rochester.edu/site/flashcode/user_support/flash_ug_devel/node187.html#SECTION010114000000000000000 """ function initial_condition_sedov_blast_wave(x, t, equations::CompressibleEulerEquations1D) # Set up polar coordinates @@ -19,7 +19,7 @@ function initial_condition_sedov_blast_wave(x, t, equations::CompressibleEulerEq x_norm = x[1] - inicenter[1] r = abs(x_norm) - # Setup based on http://flash.uchicago.edu/site/flashcode/user_support/flash_ug_devel/node184.html#SECTION010114000000000000000 + # Setup based on https://flash.rochester.edu/site/flashcode/user_support/flash_ug_devel/node187.html#SECTION010114000000000000000 r0 = 0.21875 # = 3.5 * smallest dx (for domain length=4 and max-ref=6) # r0 = 0.5 # = more reasonable setup E = 1.0 diff --git a/examples/tree_1d_dgsem/elixir_euler_sedov_blast_wave.jl b/examples/tree_1d_dgsem/elixir_euler_sedov_blast_wave.jl index 746a7cf1bac..106ccacf4f5 100644 --- a/examples/tree_1d_dgsem/elixir_euler_sedov_blast_wave.jl +++ b/examples/tree_1d_dgsem/elixir_euler_sedov_blast_wave.jl @@ -11,7 +11,7 @@ equations = CompressibleEulerEquations1D(1.4) initial_condition_sedov_blast_wave(x, t, equations::CompressibleEulerEquations1D) The Sedov blast wave setup based on Flash -- http://flash.uchicago.edu/site/flashcode/user_support/flash_ug_devel/node184.html#SECTION010114000000000000000 +- https://flash.rochester.edu/site/flashcode/user_support/flash_ug_devel/node187.html#SECTION010114000000000000000 """ function initial_condition_sedov_blast_wave(x, t, equations::CompressibleEulerEquations1D) # Set up polar coordinates @@ -19,7 +19,7 @@ function initial_condition_sedov_blast_wave(x, t, equations::CompressibleEulerEq x_norm = x[1] - inicenter[1] r = abs(x_norm) - # Setup based on http://flash.uchicago.edu/site/flashcode/user_support/flash_ug_devel/node184.html#SECTION010114000000000000000 + # Setup based on https://flash.rochester.edu/site/flashcode/user_support/flash_ug_devel/node187.html#SECTION010114000000000000000 r0 = 0.21875 # = 3.5 * smallest dx (for domain length=4 and max-ref=6) # r0 = 0.5 # = more reasonable setup E = 1.0 diff --git a/examples/tree_1d_dgsem/elixir_euler_sedov_blast_wave_pure_fv.jl b/examples/tree_1d_dgsem/elixir_euler_sedov_blast_wave_pure_fv.jl index 00b80dbae92..ebe8fa7cebf 100644 --- a/examples/tree_1d_dgsem/elixir_euler_sedov_blast_wave_pure_fv.jl +++ b/examples/tree_1d_dgsem/elixir_euler_sedov_blast_wave_pure_fv.jl @@ -11,7 +11,7 @@ equations = CompressibleEulerEquations1D(1.4) initial_condition_sedov_blast_wave(x, t, equations::CompressibleEulerEquations1D) The Sedov blast wave setup based on Flash -- http://flash.uchicago.edu/site/flashcode/user_support/flash_ug_devel/node184.html#SECTION010114000000000000000 +- https://flash.rochester.edu/site/flashcode/user_support/flash_ug_devel/node187.html#SECTION010114000000000000000 """ function initial_condition_sedov_blast_wave(x, t, equations::CompressibleEulerEquations1D) # Set up polar coordinates @@ -19,7 +19,7 @@ function initial_condition_sedov_blast_wave(x, t, equations::CompressibleEulerEq x_norm = x[1] - inicenter[1] r = abs(x_norm) - # Setup based on http://flash.uchicago.edu/site/flashcode/user_support/flash_ug_devel/node184.html#SECTION010114000000000000000 + # Setup based on https://flash.rochester.edu/site/flashcode/user_support/flash_ug_devel/node187.html#SECTION010114000000000000000 r0 = 0.21875 # = 3.5 * smallest dx (for domain length=4 and max-ref=6) # r0 = 0.5 # = more reasonable setup E = 1.0 diff --git a/examples/tree_2d_dgsem/elixir_euler_positivity.jl b/examples/tree_2d_dgsem/elixir_euler_positivity.jl index e40dc3b47af..4c7dd7eb6cf 100644 --- a/examples/tree_2d_dgsem/elixir_euler_positivity.jl +++ b/examples/tree_2d_dgsem/elixir_euler_positivity.jl @@ -11,7 +11,7 @@ equations = CompressibleEulerEquations2D(gamma) initial_condition_sedov_blast_wave(x, t, equations::CompressibleEulerEquations2D) The Sedov blast wave setup based on Flash -- http://flash.uchicago.edu/site/flashcode/user_support/flash_ug_devel/node184.html#SECTION010114000000000000000 +- https://flash.rochester.edu/site/flashcode/user_support/flash_ug_devel/node187.html#SECTION010114000000000000000 """ function initial_condition_sedov_blast_wave(x, t, equations::CompressibleEulerEquations2D) # Set up polar coordinates @@ -20,7 +20,7 @@ function initial_condition_sedov_blast_wave(x, t, equations::CompressibleEulerEq y_norm = x[2] - inicenter[2] r = sqrt(x_norm^2 + y_norm^2) - # Setup based on http://flash.uchicago.edu/site/flashcode/user_support/flash_ug_devel/node184.html#SECTION010114000000000000000 + # Setup based on https://flash.rochester.edu/site/flashcode/user_support/flash_ug_devel/node187.html#SECTION010114000000000000000 r0 = 0.21875 # = 3.5 * smallest dx (for domain length=4 and max-ref=6) # r0 = 0.5 # = more reasonable setup E = 1.0 diff --git a/examples/tree_2d_dgsem/elixir_euler_sedov_blast_wave.jl b/examples/tree_2d_dgsem/elixir_euler_sedov_blast_wave.jl index da7e1d55c91..512e5822374 100644 --- a/examples/tree_2d_dgsem/elixir_euler_sedov_blast_wave.jl +++ b/examples/tree_2d_dgsem/elixir_euler_sedov_blast_wave.jl @@ -11,7 +11,7 @@ equations = CompressibleEulerEquations2D(gamma) initial_condition_sedov_blast_wave(x, t, equations::CompressibleEulerEquations2D) The Sedov blast wave setup based on Flash -- http://flash.uchicago.edu/site/flashcode/user_support/flash_ug_devel/node184.html#SECTION010114000000000000000 +- https://flash.rochester.edu/site/flashcode/user_support/flash_ug_devel/node187.html#SECTION010114000000000000000 """ function initial_condition_sedov_blast_wave(x, t, equations::CompressibleEulerEquations2D) # Set up polar coordinates @@ -20,7 +20,7 @@ function initial_condition_sedov_blast_wave(x, t, equations::CompressibleEulerEq y_norm = x[2] - inicenter[2] r = sqrt(x_norm^2 + y_norm^2) - # Setup based on http://flash.uchicago.edu/site/flashcode/user_support/flash_ug_devel/node184.html#SECTION010114000000000000000 + # Setup based on https://flash.rochester.edu/site/flashcode/user_support/flash_ug_devel/node187.html#SECTION010114000000000000000 r0 = 0.21875 # = 3.5 * smallest dx (for domain length=4 and max-ref=6) # r0 = 0.5 # = more reasonable setup E = 1.0 diff --git a/examples/tree_2d_dgsem/elixir_euler_sedov_blast_wave_neuralnetwork_perssonperaire.jl b/examples/tree_2d_dgsem/elixir_euler_sedov_blast_wave_neuralnetwork_perssonperaire.jl index 56715789377..5fd32da2e5c 100644 --- a/examples/tree_2d_dgsem/elixir_euler_sedov_blast_wave_neuralnetwork_perssonperaire.jl +++ b/examples/tree_2d_dgsem/elixir_euler_sedov_blast_wave_neuralnetwork_perssonperaire.jl @@ -23,7 +23,7 @@ equations = CompressibleEulerEquations2D(gamma) initial_condition_sedov_blast_wave(x, t, equations::CompressibleEulerEquations2D) The Sedov blast wave setup based on Flash -- http://flash.uchicago.edu/site/flashcode/user_support/flash_ug_devel/node184.html#SECTION010114000000000000000 +- https://flash.rochester.edu/site/flashcode/user_support/flash_ug_devel/node187.html#SECTION010114000000000000000 """ function initial_condition_sedov_blast_wave(x, t, equations::CompressibleEulerEquations2D) # Set up polar coordinates @@ -32,7 +32,7 @@ function initial_condition_sedov_blast_wave(x, t, equations::CompressibleEulerEq y_norm = x[2] - inicenter[2] r = sqrt(x_norm^2 + y_norm^2) - # Setup based on http://flash.uchicago.edu/site/flashcode/user_support/flash_ug_devel/node184.html#SECTION010114000000000000000 + # Setup based on https://flash.rochester.edu/site/flashcode/user_support/flash_ug_devel/node187.html#SECTION010114000000000000000 r0 = 0.21875 # = 3.5 * smallest dx (for domain length=4 and max-ref=6) # r0 = 0.5 # = more reasonable setup E = 1.0 diff --git a/examples/tree_3d_dgsem/elixir_euler_sedov_blast_wave.jl b/examples/tree_3d_dgsem/elixir_euler_sedov_blast_wave.jl index 336c09e9212..3641878149a 100644 --- a/examples/tree_3d_dgsem/elixir_euler_sedov_blast_wave.jl +++ b/examples/tree_3d_dgsem/elixir_euler_sedov_blast_wave.jl @@ -15,14 +15,14 @@ Adaptation of the Sedov blast wave with self-gravity taken from A purely hyperbolic discontinuous Galerkin approach for self-gravitating gas dynamics [arXiv: 2008.10593](https://arxiv.org/abs/2008.10593) based on -- http://flash.uchicago.edu/site/flashcode/user_support/flash4_ug_4p62/node184.html#SECTION010114000000000000000 +- https://flash.rochester.edu/site/flashcode/user_support/flash_ug_devel/node187.html#SECTION010114000000000000000 Should be used together with [`boundary_condition_sedov_self_gravity`](@ref). """ function initial_condition_sedov_self_gravity(x, t, equations::CompressibleEulerEquations3D) # Calculate radius as distance from origin r = sqrt(x[1]^2 + x[2]^2 + x[3]^2) - # Setup based on http://flash.uchicago.edu/site/flashcode/user_support/flash4_ug_4p62/node184.html#SECTION010114000000000000000 + # Setup based on https://flash.rochester.edu/site/flashcode/user_support/flash_ug_devel/node187.html#SECTION010114000000000000000 r0 = 0.25 # = 4.0 * smallest dx (for domain length=8 and max-ref=7) E = 1.0 p_inner = (equations.gamma - 1) * E / (4/3 * pi * r0^3) @@ -60,7 +60,7 @@ Adaptation of the Sedov blast wave with self-gravity taken from A purely hyperbolic discontinuous Galerkin approach for self-gravitating gas dynamics [arXiv: 2008.10593](https://arxiv.org/abs/2008.10593) based on -- http://flash.uchicago.edu/site/flashcode/user_support/flash4_ug_4p62/node184.html#SECTION010114000000000000000 +- https://flash.rochester.edu/site/flashcode/user_support/flash_ug_devel/node187.html#SECTION010114000000000000000 Should be used together with [`initial_condition_sedov_self_gravity`](@ref). """ function boundary_condition_sedov_self_gravity(u_inner, orientation, direction, x, t, diff --git a/examples/unstructured_2d_dgsem/elixir_euler_sedov.jl b/examples/unstructured_2d_dgsem/elixir_euler_sedov.jl index 3d5a391bd90..570a2084691 100644 --- a/examples/unstructured_2d_dgsem/elixir_euler_sedov.jl +++ b/examples/unstructured_2d_dgsem/elixir_euler_sedov.jl @@ -11,7 +11,7 @@ equations = CompressibleEulerEquations2D(1.4) initial_condition_sedov_blast_wave(x, t, equations::CompressibleEulerEquations2D) The Sedov blast wave setup based on Flash -- http://flash.uchicago.edu/site/flashcode/user_support/flash_ug_devel/node184.html#SECTION010114000000000000000 +- https://flash.rochester.edu/site/flashcode/user_support/flash_ug_devel/node187.html#SECTION010114000000000000000 """ function initial_condition_sedov_blast_wave(x, t, equations::CompressibleEulerEquations2D) # Set up polar coordinates @@ -20,7 +20,7 @@ function initial_condition_sedov_blast_wave(x, t, equations::CompressibleEulerEq y_norm = x[2] - inicenter[2] r = sqrt(x_norm^2 + y_norm^2) - # Setup based on http://flash.uchicago.edu/site/flashcode/user_support/flash_ug_devel/node184.html#SECTION010114000000000000000 + # Setup based on https://flash.rochester.edu/site/flashcode/user_support/flash_ug_devel/node187.html#SECTION010114000000000000000 r0 = 0.21875 # = 3.5 * smallest dx (for domain length=4 and max-ref=6) E = 1.0 p0_inner = 3 * (equations.gamma - 1) * E / (3 * pi * r0^2) From bea4bfed27ce896edcfa3dbe79861a4faef4dec7 Mon Sep 17 00:00:00 2001 From: Hendrik Ranocha Date: Thu, 22 Jun 2023 17:06:19 +0200 Subject: [PATCH 15/15] `splitting_lax_friedrichs` for `LinearScalarAdvection1D` (#1546) * splitting_lax_friedrichs for LinearScalarAdvection1D * Update src/equations/linear_scalar_advection_1d.jl Co-authored-by: Andrew Winters --------- Co-authored-by: Andrew Winters --- .../tree_1d_fdsbp/elixir_advection_upwind.jl | 57 +++++++++++++++++++ src/equations/inviscid_burgers_1d.jl | 2 +- src/equations/linear_scalar_advection_1d.jl | 38 +++++++++++++ test/test_tree_1d_fdsbp.jl | 18 ++++++ 4 files changed, 114 insertions(+), 1 deletion(-) create mode 100644 examples/tree_1d_fdsbp/elixir_advection_upwind.jl diff --git a/examples/tree_1d_fdsbp/elixir_advection_upwind.jl b/examples/tree_1d_fdsbp/elixir_advection_upwind.jl new file mode 100644 index 00000000000..5c50e1a6c64 --- /dev/null +++ b/examples/tree_1d_fdsbp/elixir_advection_upwind.jl @@ -0,0 +1,57 @@ +# !!! warning "Experimental implementation (upwind SBP)" +# This is an experimental feature and may change in future releases. + +using OrdinaryDiffEq +using Trixi + +############################################################################### +# semidiscretization of the linear scalar advection equation equation + +equations = LinearScalarAdvectionEquation1D(1.0) + +function initial_condition_sin(x, t, equation::LinearScalarAdvectionEquation1D) + return SVector(sinpi(x[1] - equations.advection_velocity[1] * t)) +end + +D_upw = upwind_operators(SummationByPartsOperators.Mattsson2017, + derivative_order = 1, + accuracy_order = 4, + xmin = -1.0, xmax = 1.0, + N = 16) +flux_splitting = splitting_lax_friedrichs +solver = FDSBP(D_upw, + surface_integral = SurfaceIntegralUpwind(flux_splitting), + volume_integral = VolumeIntegralUpwind(flux_splitting)) + +coordinates_min = -1.0 +coordinates_max = 1.0 +mesh = TreeMesh(coordinates_min, coordinates_max, + initial_refinement_level = 4, + n_cells_max = 10_000) + +semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition_sin, solver) + + +############################################################################### +# ODE solvers, callbacks etc. + +tspan = (0.0, 2.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) + +callbacks = CallbackSet(summary_callback, + analysis_callback, alive_callback) + + +############################################################################### +# run the simulation + +sol = solve(ode, RDPK3SpFSAL49(); abstol=1.0e-6, reltol=1.0e-6, + ode_default_options()..., callback=callbacks); +summary_callback() # print the timer summary diff --git a/src/equations/inviscid_burgers_1d.jl b/src/equations/inviscid_burgers_1d.jl index 8d4410b6ffe..6a2cfb6aa8e 100644 --- a/src/equations/inviscid_burgers_1d.jl +++ b/src/equations/inviscid_burgers_1d.jl @@ -132,7 +132,7 @@ end equations::InviscidBurgersEquation1D) Naive local Lax-Friedrichs style flux splitting of the form `f⁺ = 0.5 (f + λ u)` -and `f⁻ = 0.5 (f - λ u)` where λ = abs(u). +and `f⁻ = 0.5 (f - λ u)` where `λ = abs(u)`. Returns a tuple of the fluxes "minus" (associated with waves going into the negative axis direction) and "plus" (associated with waves going into the diff --git a/src/equations/linear_scalar_advection_1d.jl b/src/equations/linear_scalar_advection_1d.jl index 7769cb61fbf..6c6b9dd3721 100644 --- a/src/equations/linear_scalar_advection_1d.jl +++ b/src/equations/linear_scalar_advection_1d.jl @@ -172,6 +172,44 @@ end return abs.(equation.advection_velocity) end +""" + splitting_lax_friedrichs(u, orientation::Integer, + equations::LinearScalarAdvectionEquation1D) + splitting_lax_friedrichs(u, which::Union{Val{:minus}, Val{:plus}} + orientation::Integer, + equations::LinearScalarAdvectionEquation1D) + +Naive local Lax-Friedrichs style flux splitting of the form `f⁺ = 0.5 (f + λ u)` +and `f⁻ = 0.5 (f - λ u)` where `λ` is the absolute value of the advection +velocity. + +Returns a tuple of the fluxes "minus" (associated with waves going into the +negative axis direction) and "plus" (associated with waves going into the +positive axis direction). If only one of the fluxes is required, use the +function signature with argument `which` set to `Val{:minus}()` or `Val{:plus}()`. + +!!! warning "Experimental implementation (upwind SBP)" + This is an experimental feature and may change in future releases. +""" +@inline function splitting_lax_friedrichs(u, orientation::Integer, + equations::LinearScalarAdvectionEquation1D) + fm = splitting_lax_friedrichs(u, Val{:minus}(), orientation, equations) + fp = splitting_lax_friedrichs(u, Val{:plus}(), orientation, equations) + return fm, fp +end + +@inline function splitting_lax_friedrichs(u, ::Val{:plus}, orientation::Integer, + equations::LinearScalarAdvectionEquation1D) + a = equations.advection_velocity[1] + return a > 0 ? flux(u, orientation, equations) : zero(u) +end + +@inline function splitting_lax_friedrichs(u, ::Val{:minus}, orientation::Integer, + equations::LinearScalarAdvectionEquation1D) + a = equations.advection_velocity[1] + return a < 0 ? flux(u, orientation, equations) : zero(u) +end + # Convert conservative variables to primitive @inline cons2prim(u, equation::LinearScalarAdvectionEquation1D) = u diff --git a/test/test_tree_1d_fdsbp.jl b/test/test_tree_1d_fdsbp.jl index a966b3836f3..118385c34b3 100644 --- a/test/test_tree_1d_fdsbp.jl +++ b/test/test_tree_1d_fdsbp.jl @@ -7,6 +7,24 @@ include("test_trixi.jl") EXAMPLES_DIR = pkgdir(Trixi, "examples", "tree_1d_fdsbp") +@testset "Linear scalar advection" begin + @trixi_testset "elixir_advection_upwind.jl" begin + @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_advection_upwind.jl"), + l2 = [1.7735637157305526e-6], + linf = [1.0418854521951328e-5], + tspan = (0.0, 0.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 +end + @testset "Inviscid Burgers" begin @trixi_testset "elixir_burgers_basic.jl" begin @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_burgers_basic.jl"),