Central SBP finite difference solver for UnstructuredMesh2D (trixi-…

…framework#1773) * containers and kernels for FDSBP solver on UnstructuredMesh2D * add elixirs and corresponding tests * apply formatter to new and edited files * add advection equation test to up coverage * Apply suggestions from code review Co-authored-by: Hendrik Ranocha <[email protected]> * update variable name to N --------- Co-authored-by: Hendrik Ranocha <[email protected]>
amrueda · Dec 13, 2023 · 14796ea · 14796ea
1 parent 766c798
commit 14796ea
Show file tree

Hide file tree

Showing 10 changed files with 643 additions and 11 deletions.
diff --git a/examples/unstructured_2d_fdsbp/elixir_advection_basic.jl b/examples/unstructured_2d_fdsbp/elixir_advection_basic.jl
@@ -0,0 +1,69 @@
+
+using Downloads: download
+using OrdinaryDiffEq
+using Trixi
+
+###############################################################################
+# semidiscretization of the linear advection equation
+
+advection_velocity = (0.2, -0.7)
+equations = LinearScalarAdvectionEquation2D(advection_velocity)
+
+###############################################################################
+# Get the FDSBP approximation operator
+
+D_SBP = derivative_operator(SummationByPartsOperators.MattssonAlmquistVanDerWeide2018Accurate(),
+                            derivative_order = 1, accuracy_order = 4,
+                            xmin = -1.0, xmax = 1.0, N = 15)
+solver = FDSBP(D_SBP,
+               surface_integral = SurfaceIntegralStrongForm(flux_lax_friedrichs),
+               volume_integral = VolumeIntegralStrongForm())
+
+###############################################################################
+# Get the curved quad mesh from a file (downloads the file if not available locally)
+
+default_mesh_file = joinpath(@__DIR__, "mesh_periodic_square_with_twist.mesh")
+isfile(default_mesh_file) ||
+    download("https://gist.githubusercontent.com/andrewwinters5000/12ce661d7c354c3d94c74b964b0f1c96/raw/8275b9a60c6e7ebbdea5fc4b4f091c47af3d5273/mesh_periodic_square_with_twist.mesh",
+             default_mesh_file)
+mesh_file = default_mesh_file
+
+mesh = UnstructuredMesh2D(mesh_file, periodicity = true)
+
+###############################################################################
+# create the semidiscretization object
+
+semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition_convergence_test,
+                                    solver)
+
+###############################################################################
+# ODE solvers, callbacks etc.
+
+# Create ODE problem with time span from 0.0 to 1.0
+ode = semidiscretize(semi, (0.0, 1.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_callback = AnalysisCallback(semi, interval = 100)
+
+# 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
+
+sol = solve(ode, CarpenterKennedy2N54(williamson_condition = false),
+            dt = 1.0, # solve needs some value here but it will be overwritten by the stepsize_callback
+            save_everystep = false, callback = callbacks);
+summary_callback() # print the timer summary
diff --git a/examples/unstructured_2d_fdsbp/elixir_euler_free_stream.jl b/examples/unstructured_2d_fdsbp/elixir_euler_free_stream.jl
@@ -0,0 +1,77 @@
+
+using Downloads: download
+using OrdinaryDiffEq
+using Trixi
+
+###############################################################################
+# semidiscretization of the compressible Euler equations
+
+equations = CompressibleEulerEquations2D(1.4)
+
+# Free-stream initial condition
+initial_condition = initial_condition_constant
+
+# Boundary conditions for free-stream testing
+boundary_condition_free_stream = BoundaryConditionDirichlet(initial_condition)
+boundary_conditions = Dict(:Body => boundary_condition_free_stream,
+                           :Button1 => boundary_condition_free_stream,
+                           :Button2 => boundary_condition_free_stream,
+                           :Eye1 => boundary_condition_free_stream,
+                           :Eye2 => boundary_condition_free_stream,
+                           :Smile => boundary_condition_free_stream,
+                           :Bowtie => boundary_condition_free_stream)
+
+###############################################################################
+# Get the FDSBP approximation space
+
+D_SBP = derivative_operator(SummationByPartsOperators.MattssonAlmquistVanDerWeide2018Accurate(),
+                            derivative_order = 1, accuracy_order = 4,
+                            xmin = -1.0, xmax = 1.0, N = 12)
+solver = FDSBP(D_SBP,
+               surface_integral = SurfaceIntegralStrongForm(flux_hll),
+               volume_integral = VolumeIntegralStrongForm())
+
+###############################################################################
+# Get the curved quad mesh from a file (downloads the file if not available locally)
+
+default_mesh_file = joinpath(@__DIR__, "mesh_gingerbread_man.mesh")
+isfile(default_mesh_file) ||
+    download("https://gist.githubusercontent.com/andrewwinters5000/2c6440b5f8a57db131061ad7aa78ee2b/raw/1f89fdf2c874ff678c78afb6fe8dc784bdfd421f/mesh_gingerbread_man.mesh",
+             default_mesh_file)
+mesh_file = default_mesh_file
+
+mesh = UnstructuredMesh2D(mesh_file)
+
+###############################################################################
+# create the semi discretization object
+
+semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver,
+                                    boundary_conditions = boundary_conditions)
+
+###############################################################################
+# 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)
+
+alive_callback = AliveCallback(analysis_interval = analysis_interval)
+
+save_solution = SaveSolutionCallback(interval = 100,
+                                     save_initial_solution = true,
+                                     save_final_solution = true)
+
+callbacks = CallbackSet(summary_callback, analysis_callback,
+                        alive_callback, save_solution)
+
+###############################################################################
+# run the simulation
+
+# set small tolerances for the free-stream preservation test
+sol = solve(ode, SSPRK43(), abstol = 1.0e-12, reltol = 1.0e-12,
+            save_everystep = false, callback = callbacks)
+summary_callback() # print the timer summary
diff --git a/examples/unstructured_2d_fdsbp/elixir_euler_source_terms.jl b/examples/unstructured_2d_fdsbp/elixir_euler_source_terms.jl
@@ -0,0 +1,65 @@
+
+using Downloads: download
+using OrdinaryDiffEq
+using Trixi
+
+###############################################################################
+# semidiscretization of the compressible Euler equations
+
+equations = CompressibleEulerEquations2D(1.4)
+
+initial_condition = initial_condition_convergence_test
+
+###############################################################################
+# Get the FDSBP approximation operator
+
+D_SBP = derivative_operator(SummationByPartsOperators.MattssonNordström2004(),
+                            derivative_order = 1, accuracy_order = 4,
+                            xmin = -1.0, xmax = 1.0, N = 10)
+solver = FDSBP(D_SBP,
+               surface_integral = SurfaceIntegralStrongForm(flux_lax_friedrichs),
+               volume_integral = VolumeIntegralStrongForm())
+
+###############################################################################
+# Get the curved quad mesh from a file (downloads the file if not available locally)
+
+default_mesh_file = joinpath(@__DIR__, "mesh_periodic_square_with_twist.mesh")
+isfile(default_mesh_file) ||
+    download("https://gist.githubusercontent.com/andrewwinters5000/12ce661d7c354c3d94c74b964b0f1c96/raw/8275b9a60c6e7ebbdea5fc4b4f091c47af3d5273/mesh_periodic_square_with_twist.mesh",
+             default_mesh_file)
+mesh_file = default_mesh_file
+
+mesh = UnstructuredMesh2D(mesh_file, periodicity = true)
+
+###############################################################################
+# create the semi discretization object
+
+semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver,
+                                    source_terms = source_terms_convergence_test)
+
+###############################################################################
+# ODE solvers, callbacks etc.
+
+tspan = (0.0, 1.0)
+ode = semidiscretize(semi, tspan)
+
+summary_callback = SummaryCallback()
+
+analysis_interval = 100
+analysis_callback = AnalysisCallback(semi, interval = analysis_interval)
+
+alive_callback = AliveCallback(analysis_interval = analysis_interval)
+
+save_solution = SaveSolutionCallback(interval = 100,
+                                     save_initial_solution = true,
+                                     save_final_solution = true)
+
+callbacks = CallbackSet(summary_callback, analysis_callback,
+                        alive_callback, save_solution)
+
+###############################################################################
+# run the simulation
+
+sol = solve(ode, SSPRK43(), abstol = 1.0e-9, reltol = 1.0e-9,
+            save_everystep = false, callback = callbacks)
+summary_callback() # print the timer summary
diff --git a/src/solvers/dg.jl b/src/solvers/dg.jl
@@ -41,8 +41,8 @@ standard textbooks.
   Applications
   [doi: 10.1007/978-0-387-72067-8](https://doi.org/10.1007/978-0-387-72067-8)
 
-`VolumeIntegralWeakForm()` is only implemented for conserved terms as 
-non-conservative terms should always be discretized in conjunction with a flux-splitting scheme, 
+`VolumeIntegralWeakForm()` is only implemented for conserved terms as
+non-conservative terms should always be discretized in conjunction with a flux-splitting scheme,
 see [`VolumeIntegralFluxDifferencing`](@ref).
 This treatment is required to achieve, e.g., entropy-stability or well-balancedness.
 """
@@ -415,7 +415,8 @@ function Base.show(io::IO, mime::MIME"text/plain", dg::DG)
         summary_line(io, "surface integral", dg.surface_integral |> typeof |> nameof)
         show(increment_indent(io), mime, dg.surface_integral)
         summary_line(io, "volume integral", dg.volume_integral |> typeof |> nameof)
-        if !(dg.volume_integral isa VolumeIntegralWeakForm)
+        if !(dg.volume_integral isa VolumeIntegralWeakForm) &&
+           !(dg.volume_integral isa VolumeIntegralStrongForm)
             show(increment_indent(io), mime, dg.volume_integral)
         end
         summary_footer(io)
@@ -598,6 +599,7 @@ include("dgsem/dgsem.jl")
 # and boundary conditions weakly. Thus, these methods can re-use a lot of
 # functionality implemented for DGSEM.
 include("fdsbp_tree/fdsbp.jl")
+include("fdsbp_unstructured/fdsbp.jl")
 
 function allocate_coefficients(mesh::AbstractMesh, equations, dg::DG, cache)
     # We must allocate a `Vector` in order to be able to `resize!` it (AMR).

diff --git a/src/solvers/dgsem_unstructured/containers_2d.jl b/src/solvers/dgsem_unstructured/containers_2d.jl
@@ -45,7 +45,7 @@ end
     eachelement(elements::UnstructuredElementContainer2D)
 
 Return an iterator over the indices that specify the location in relevant data structures
-for the elements in `elements`. 
+for the elements in `elements`.
 In particular, not the elements themselves are returned.
 """
 @inline function eachelement(elements::UnstructuredElementContainer2D)
@@ -84,32 +84,33 @@ function init_elements!(elements::UnstructuredElementContainer2D, mesh, basis)
     # loop through elements and call the correct constructor based on whether the element is curved
     for element in eachelement(elements)
         if mesh.element_is_curved[element]
-            init_element!(elements, element, basis.nodes,
+            init_element!(elements, element, basis,
                           view(mesh.surface_curves, :, element))
         else # straight sided element
             for i in 1:4, j in 1:2
                 # pull the (x,y) values of these corners out of the global corners array
                 four_corners[i, j] = mesh.corners[j, mesh.element_node_ids[i, element]]
             end
-            init_element!(elements, element, basis.nodes, four_corners)
+            init_element!(elements, element, basis, four_corners)
         end
     end
 end
 
 # initialize all the values in the container of a general element (either straight sided or curved)
-function init_element!(elements, element, nodes, corners_or_surface_curves)
-    calc_node_coordinates!(elements.node_coordinates, element, nodes,
+function init_element!(elements, element, basis::LobattoLegendreBasis,
+                       corners_or_surface_curves)
+    calc_node_coordinates!(elements.node_coordinates, element, get_nodes(basis),
                            corners_or_surface_curves)
 
-    calc_metric_terms!(elements.jacobian_matrix, element, nodes,
+    calc_metric_terms!(elements.jacobian_matrix, element, get_nodes(basis),
                        corners_or_surface_curves)
 
     calc_inverse_jacobian!(elements.inverse_jacobian, element, elements.jacobian_matrix)
 
     calc_contravariant_vectors!(elements.contravariant_vectors, element,
                                 elements.jacobian_matrix)
 
-    calc_normal_directions!(elements.normal_directions, element, nodes,
+    calc_normal_directions!(elements.normal_directions, element, get_nodes(basis),
                             corners_or_surface_curves)
 
     return elements

diff --git a/src/solvers/fdsbp_tree/fdsbp_2d.jl b/src/solvers/fdsbp_tree/fdsbp_2d.jl
@@ -9,7 +9,7 @@
 #! format: noindent
 
 # 2D caches
-function create_cache(mesh::TreeMesh{2}, equations,
+function create_cache(mesh::Union{TreeMesh{2}, UnstructuredMesh2D}, equations,
                       volume_integral::VolumeIntegralStrongForm, dg, uEltype)
     prototype = Array{SVector{nvariables(equations), uEltype}, ndims(mesh)}(undef,
                                                                             ntuple(_ -> nnodes(dg),