Upwind SBP on curved meshes (#1857)

* baseline implementation of the curvilinear USBP for testing * working version of curvilinear upwind solver. Needs significant cleanup of debugging statements and different variants of the rotated flux vector splittings * cleanup of the fdsbp_2d file * clean-up FVS routines in the compressible Euler file * cleanup and remove unnecessary containers * add tests for the new solver * remove extra space * run formatter * Apply suggestions from code review Co-authored-by: Hendrik Ranocha <[email protected]> * add specialized calc_metric_terms function for upwind type * revert change to the surface integral function * add reference for curvilinear van Leer splitting * new splitting_drikakis_tsangaris in Cartesian and generalized coordinates * added test for Cartesian splitting_drikakis_tsangaris * run formatter * Update src/equations/compressible_euler_2d.jl * remove orientation_or_normal from Steger-Warming --------- Co-authored-by: Hendrik Ranocha <[email protected]>
trixi-framework · Mar 6, 2024 · 8cdb938 · 8cdb938
1 parent e205b06
commit 8cdb938
Show file tree

Hide file tree

Showing 11 changed files with 659 additions and 28 deletions.
diff --git a/examples/unstructured_2d_fdsbp/elixir_euler_free_stream_upwind.jl b/examples/unstructured_2d_fdsbp/elixir_euler_free_stream_upwind.jl
@@ -0,0 +1,86 @@
+# !!! warning "Experimental implementation (upwind SBP)"
+#     This is an experimental feature and may change in future releases.
+
+using OrdinaryDiffEq
+using Trixi
+
+###############################################################################
+# semidiscretization of the compressible Euler equations
+
+equations = CompressibleEulerEquations2D(1.4)
+
+initial_condition = initial_condition_constant
+
+# Boundary conditions for free-stream preservation test
+boundary_condition_free_stream = BoundaryConditionDirichlet(initial_condition)
+
+boundary_conditions = Dict(:outerCircle => boundary_condition_free_stream,
+                           :cone1 => boundary_condition_free_stream,
+                           :cone2 => boundary_condition_free_stream,
+                           :iceCream => boundary_condition_free_stream)
+
+###############################################################################
+# Get the Upwind FDSBP approximation space
+
+# TODO: FDSBP
+# Note, one must set `xmin=-1` and `xmax=1` due to the reuse
+# of interpolation routines from `calc_node_coordinates!` to create
+# the physical coordinates in the mappings.
+D_upw = upwind_operators(SummationByPartsOperators.Mattsson2017,
+                         derivative_order = 1,
+                         accuracy_order = 8,
+                         xmin = -1.0, xmax = 1.0,
+                         N = 17)
+
+flux_splitting = splitting_vanleer_haenel
+solver = FDSBP(D_upw,
+               surface_integral = SurfaceIntegralStrongForm(FluxUpwind(flux_splitting)),
+               volume_integral = VolumeIntegralUpwind(flux_splitting))
+
+###############################################################################
+# Get the curved quad mesh from a file (downloads the file if not available locally)
+
+# Mesh with second-order boundary polynomials requires an upwind SBP operator
+# with (at least) 4th order boundary closure to guarantee the approximation is
+# free-stream preserving
+mesh_file = Trixi.download("https://gist.githubusercontent.com/andrewwinters5000/ec9a345f09199ebe471d35d5c1e4e08f/raw/15975943d8642e42f8292235314b6f1b30aa860d/mesh_inner_outer_boundaries.mesh",
+                           joinpath(@__DIR__, "mesh_inner_outer_boundaries.mesh"))
+
+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 = 1000
+analysis_callback = AnalysisCallback(semi, interval = analysis_interval)
+
+alive_callback = AliveCallback(analysis_interval = analysis_interval)
+
+save_solution = SaveSolutionCallback(interval = 1000,
+                                     save_initial_solution = true,
+                                     save_final_solution = true)
+
+callbacks = CallbackSet(summary_callback,
+                        analysis_callback,
+                        save_solution,
+                        alive_callback)
+
+###############################################################################
+# 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_upwind.jl b/examples/unstructured_2d_fdsbp/elixir_euler_source_terms_upwind.jl
@@ -0,0 +1,87 @@
+# !!! warning "Experimental implementation (upwind SBP)"
+#     This is an experimental feature and may change in future releases.
+
+using OrdinaryDiffEq
+using Trixi
+
+###############################################################################
+# semidiscretization of the compressible Euler equations
+
+equations = CompressibleEulerEquations2D(1.4)
+
+initial_condition = initial_condition_convergence_test
+
+source_term = source_terms_convergence_test
+
+boundary_condition_eoc = BoundaryConditionDirichlet(initial_condition)
+
+boundary_conditions = Dict(:Top => boundary_condition_eoc,
+                           :Bottom => boundary_condition_eoc,
+                           :Right => boundary_condition_eoc,
+                           :Left => boundary_condition_eoc)
+
+###############################################################################
+# Get the Upwind FDSBP approximation space
+
+# TODO: FDSBP
+# Note, one must set `xmin=-1` and `xmax=1` due to the reuse
+# of interpolation routines from `calc_node_coordinates!` to create
+# the physical coordinates in the mappings.
+D_upw = upwind_operators(SummationByPartsOperators.Mattsson2017,
+                         derivative_order = 1,
+                         accuracy_order = 4,
+                         xmin = -1.0, xmax = 1.0,
+                         N = 9)
+
+flux_splitting = splitting_drikakis_tsangaris
+solver = FDSBP(D_upw,
+               surface_integral = SurfaceIntegralStrongForm(FluxUpwind(flux_splitting)),
+               volume_integral = VolumeIntegralUpwind(flux_splitting))
+
+###############################################################################
+# Get the curved quad mesh from a file (downloads the file if not available locally)
+
+# Mesh with first-order boundary polynomials requires an upwind SBP operator
+# with (at least) 2nd order boundary closure to guarantee the approximation is
+# free-stream preserving
+mesh_file = Trixi.download("https://gist.githubusercontent.com/andrewwinters5000/a4f4743008bf3233957a9ea6ac7a62e0/raw/8b36cc6649153fe0a5723b200368a210a1d74eaf/mesh_refined_box.mesh",
+                           joinpath(@__DIR__, "mesh_refined_box.mesh"))
+
+mesh = UnstructuredMesh2D(mesh_file)
+
+###############################################################################
+# create the semidiscretization object
+
+semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver,
+                                    source_terms = source_term,
+                                    boundary_conditions = boundary_conditions)
+
+###############################################################################
+# ODE solvers, callbacks etc.
+
+tspan = (0.0, 1.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)
+
+save_solution = SaveSolutionCallback(interval = 1000,
+                                     save_initial_solution = true,
+                                     save_final_solution = true)
+
+callbacks = CallbackSet(summary_callback,
+                        analysis_callback,
+                        save_solution,
+                        alive_callback)
+
+###############################################################################
+# run the simulation
+
+sol = solve(ode, SSPRK43(), abstol = 1.0e-6, reltol = 1.0e-6,
+            save_everystep = false, callback = callbacks)
+
+summary_callback() # print the timer summary
diff --git a/src/Trixi.jl b/src/Trixi.jl
@@ -191,7 +191,8 @@ export flux, flux_central, flux_lax_friedrichs, flux_hll, flux_hllc, flux_hlle,
        FluxUpwind
 
 export splitting_steger_warming, splitting_vanleer_haenel,
-       splitting_coirier_vanleer, splitting_lax_friedrichs
+       splitting_coirier_vanleer, splitting_lax_friedrichs,
+       splitting_drikakis_tsangaris
 
 export initial_condition_constant,
        initial_condition_gauss,