WIP: Sc/polytropic 2d wave speed clean

Project.toml.bkp

@@ -0,0 +1,95 @@
+name = "Trixi"
+uuid = "a7f1ee26-1774-49b1-8366-f1abc58fbfcb"
+authors = ["Michael Schlottke-Lakemper <[email protected]>", "Gregor Gassner <[email protected]>", "Hendrik Ranocha <[email protected]>", "Andrew R. Winters <[email protected]>", "Jesse Chan <[email protected]>"]
+version = "0.5.28-pre"
+
+[deps]
+AbbreviatedStackTraces = "ac637c84-cc71-43bf-9c33-c1b4316be3d4"
+CodeTracking = "da1fd8a2-8d9e-5ec2-8556-3022fb5608a2"
+ConstructionBase = "187b0558-2788-49d3-abe0-74a17ed4e7c9"
+DiffEqCallbacks = "459566f4-90b8-5000-8ac3-15dfb0a30def"
+EllipsisNotation = "da5c29d0-fa7d-589e-88eb-ea29b0a81949"
+FillArrays = "1a297f60-69ca-5386-bcde-b61e274b549b"
+ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210"
+HDF5 = "f67ccb44-e63f-5c2f-98bd-6dc0ccc4ba2f"
+IfElse = "615f187c-cbe4-4ef1-ba3b-2fcf58d6d173"
+Infiltrator = "5903a43b-9cc3-4c30-8d17-598619ec4e9b"
+LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
+LinearMaps = "7a12625a-238d-50fd-b39a-03d52299707e"
+LoopVectorization = "bdcacae8-1622-11e9-2a5c-532679323890"
+MPI = "da04e1cc-30fd-572f-bb4f-1f8673147195"
+MuladdMacro = "46d2c3a1-f734-5fdb-9937-b9b9aeba4221"
+Octavian = "6fd5a793-0b7e-452c-907f-f8bfe9c57db4"
+OffsetArrays = "6fe1bfb0-de20-5000-8ca7-80f57d26f881"
+OrdinaryDiffEq = "1dea7af3-3e70-54e6-95c3-0bf5283fa5ed"
+P4est = "7d669430-f675-4ae7-b43e-fab78ec5a902"
+Polyester = "f517fe37-dbe3-4b94-8317-1923a5111588"
+PrecompileTools = "aea7be01-6a6a-4083-8856-8a6e6704d82a"
+Printf = "de0858da-6303-5e67-8744-51eddeeeb8d7"
+RecipesBase = "3cdcf5f2-1ef4-517c-9805-6587b60abb01"
+Reexport = "189a3867-3050-52da-a836-e630ba90ab69"
+Requires = "ae029012-a4dd-5104-9daa-d747884805df"
+Revise = "295af30f-e4ad-537b-8983-00126c2a3abe"
+SciMLBase = "0bca4576-84f4-4d90-8ffe-ffa030f20462"
+Setfield = "efcf1570-3423-57d1-acb7-fd33fddbac46"
+SimpleUnPack = "ce78b400-467f-4804-87d8-8f486da07d0a"
+SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
+StartUpDG = "472ebc20-7c99-4d4b-9470-8fde4e9faa0f"
+Static = "aedffcd0-7271-4cad-89d0-dc628f76c6d3"
+StaticArrayInterface = "0d7ed370-da01-4f52-bd93-41d350b8b718"
+StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
+StrideArrays = "d1fa6d79-ef01-42a6-86c9-f7c551f8593b"
+StructArrays = "09ab397b-f2b6-538f-b94a-2f83cf4a842a"
+SummationByPartsOperators = "9f78cca6-572e-554e-b819-917d2f1cf240"
+TimerOutputs = "a759f4b9-e2f1-59dc-863e-4aeb61b1ea8f"
+Triangulate = "f7e6ffb2-c36d-4f8f-a77e-16e897189344"
+TriplotBase = "981d1d27-644d-49a2-9326-4793e63143c3"
+TriplotRecipes = "808ab39a-a642-4abf-81ff-4cb34ebbffa3"
+Trixi2Vtk = "bc1476a1-1ca6-4cc3-950b-c312b255ff95"
+
+[weakdeps]
+Makie = "ee78f7c6-11fb-53f2-987a-cfe4a2b5a57a"
+
+[extensions]
+TrixiMakieExt = "Makie"
+
+[compat]
+CodeTracking = "1.0.5"
+ConstructionBase = "1.3"
+DiffEqCallbacks = "2.25"
+EllipsisNotation = "1.0"
+FillArrays = "0.13.2, 1"
+ForwardDiff = "0.10.18"
+HDF5 = "0.14, 0.15, 0.16"
+IfElse = "0.1"
+LinearMaps = "2.7, 3.0"
+LoopVectorization = "0.12.118"
+MPI = "0.20"
+Makie = "0.19"
+MuladdMacro = "0.2.2"
+Octavian = "0.3.5"
+OffsetArrays = "1.3"
+P4est = "0.4"
+Polyester = "0.3.4, 0.5, 0.6, 0.7"
+PrecompileTools = "1.1"
+RecipesBase = "1.1"
+Reexport = "1.0"
+Requires = "1.1"
+SciMLBase = "1.90"
+Setfield = "0.8, 1"
+SimpleUnPack = "1.1"
+StartUpDG = "0.16"
+Static = "0.3, 0.4, 0.5, 0.6, 0.7, 0.8"
+StaticArrayInterface = "1.4"
+StaticArrays = "1"
+StrideArrays = "0.1.18"
+StructArrays = "0.6"
+SummationByPartsOperators = "0.5.25"
+TimerOutputs = "0.5"
+Triangulate = "2.0"
+TriplotBase = "0.1"
+TriplotRecipes = "0.1"
+julia = "1.8"
+
+[extras]
+Makie = "ee78f7c6-11fb-53f2-987a-cfe4a2b5a57a"

examples/p4est_2d_dgsem/Project.toml

@@ -0,0 +1,8 @@
+[deps]
+AbbreviatedStackTraces = "ac637c84-cc71-43bf-9c33-c1b4316be3d4"
+Infiltrator = "5903a43b-9cc3-4c30-8d17-598619ec4e9b"
+OrdinaryDiffEq = "1dea7af3-3e70-54e6-95c3-0bf5283fa5ed"
+Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80"
+Revise = "295af30f-e4ad-537b-8983-00126c2a3abe"
+Trixi = "a7f1ee26-1774-49b1-8366-f1abc58fbfcb"
+Trixi2Vtk = "bc1476a1-1ca6-4cc3-950b-c312b255ff95"

examples/p4est_2d_dgsem/elixir_advection_coupled.jl

@@ -0,0 +1,142 @@
+using OrdinaryDiffEq
+using Trixi
+
+###############################################################################
+# Coupled semidiscretization of two linear advection systems using converter functions such that
+# the upper half of the domain is coupled periodically, while the lower half is not coupled
+# and any incoming wave is completely absorbed.
+#
+# 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`
+trees_per_dimension = (8, 16)
+
+trees_per_dimension1 = trees_per_dimension
+
+# Create P4estMesh with 8 x 16 trees and 16 x 32 elements
+mesh1 = P4estMesh(trees_per_dimension1, polydeg = 3,
+                  coordinates_min=coordinates_min1, coordinates_max=coordinates_max1,
+                  initial_refinement_level = 1)
+
+
+# The user can define their own coupling functions.
+coupling_function1 = (x, u) -> (sign(x[2] - 0.0)*0.1 + 1.0)/1.1 * u
+
+boundary_conditions_x_neg1 = BoundaryConditionCoupled(2, (:end, :i_forward), Float64,
+                                                      coupling_function1)
+boundary_conditions_x_pos1 = BoundaryConditionCoupled(2, (:begin, :i_forward), Float64,
+                                                      coupling_function1)
+
+# 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 = boundary_conditions_x_neg1,
+                                                            # Connect right boundary with left boundary of right mesh
+                                                            x_pos = boundary_conditions_x_pos1,
+                                                            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))
+
+trees_per_dimension2 = trees_per_dimension
+
+# Create P4estMesh with 8 x 16 trees and 16 x 32 elements
+mesh2 = P4estMesh(trees_per_dimension2, polydeg = 3,
+                  coordinates_min=coordinates_min2, coordinates_max=coordinates_max2,
+                  initial_refinement_level = 1)
+
+
+coupling_function2 = (x, u) -> (sign(x[2] - 0.0)*0.1 + 1.0)/1.1 * u
+
+boundary_conditions_x_neg2 = BoundaryConditionCoupled(1, (:end, :i_forward), Float64,
+                                                      coupling_function2)
+boundary_conditions_x_pos2 = BoundaryConditionCoupled(1, (:begin, :i_forward), Float64,
+                                                      coupling_function2)
+
+semi2 = SemidiscretizationHyperbolic(mesh2, equations, initial_condition_convergence_test,
+                                     solver,
+                                     boundary_conditions = (
+                                                            # Connect left boundary with right boundary of left mesh
+                                                            x_neg = boundary_conditions_x_neg2,
+                                                            # Connect right boundary with left boundary of left mesh
+                                                            x_pos = boundary_conditions_x_pos2,
+                                                            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, 20.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 = 1,
+                                     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()

examples/p4est_3d_dgsem/Project.toml

@@ -0,0 +1,5 @@
+[deps]
+OrdinaryDiffEq = "1dea7af3-3e70-54e6-95c3-0bf5283fa5ed"
+Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80"
+Trixi = "a7f1ee26-1774-49b1-8366-f1abc58fbfcb"
+Trixi2Vtk = "bc1476a1-1ca6-4cc3-950b-c312b255ff95"

examples/structured_2d_dgsem/Project.toml

@@ -0,0 +1,5 @@
+[deps]
+OrdinaryDiffEq = "1dea7af3-3e70-54e6-95c3-0bf5283fa5ed"
+Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80"
+Trixi = "a7f1ee26-1774-49b1-8366-f1abc58fbfcb"
+Trixi2Vtk = "bc1476a1-1ca6-4cc3-950b-c312b255ff95"

examples/structured_3d_dgsem/Project.toml

@@ -0,0 +1,5 @@
+[deps]
+OrdinaryDiffEq = "1dea7af3-3e70-54e6-95c3-0bf5283fa5ed"
+Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80"
+Trixi = "a7f1ee26-1774-49b1-8366-f1abc58fbfcb"
+Trixi2Vtk = "bc1476a1-1ca6-4cc3-950b-c312b255ff95"

examples/tree_1d_dgsem/Project.toml

@@ -0,0 +1,5 @@
+[deps]
+OrdinaryDiffEq = "1dea7af3-3e70-54e6-95c3-0bf5283fa5ed"
+Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80"
+Trixi = "a7f1ee26-1774-49b1-8366-f1abc58fbfcb"
+Trixi2Vtk = "bc1476a1-1ca6-4cc3-950b-c312b255ff95"

examples/tree_2d_dgsem/Project.toml

@@ -0,0 +1,8 @@
+[deps]
+AbbreviatedStackTraces = "ac637c84-cc71-43bf-9c33-c1b4316be3d4"
+Infiltrator = "5903a43b-9cc3-4c30-8d17-598619ec4e9b"
+OrdinaryDiffEq = "1dea7af3-3e70-54e6-95c3-0bf5283fa5ed"
+Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80"
+Revise = "295af30f-e4ad-537b-8983-00126c2a3abe"
+Trixi = "a7f1ee26-1774-49b1-8366-f1abc58fbfcb"
+Trixi2Vtk = "bc1476a1-1ca6-4cc3-950b-c312b255ff95"

examples/tree_3d_dgsem/Project.toml

@@ -0,0 +1,5 @@
+[deps]
+OrdinaryDiffEq = "1dea7af3-3e70-54e6-95c3-0bf5283fa5ed"
+Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80"
+Trixi = "a7f1ee26-1774-49b1-8366-f1abc58fbfcb"
+Trixi2Vtk = "bc1476a1-1ca6-4cc3-950b-c312b255ff95"

examples/unstructured_2d_dgsem/Project.toml

@@ -0,0 +1,5 @@
+[deps]
+OrdinaryDiffEq = "1dea7af3-3e70-54e6-95c3-0bf5283fa5ed"
+Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80"
+Trixi = "a7f1ee26-1774-49b1-8366-f1abc58fbfcb"
+Trixi2Vtk = "bc1476a1-1ca6-4cc3-950b-c312b255ff95"

src/Trixi.jl

@@ -125,6 +125,7 @@ include("time_integration/time_integration.jl")
 include("callbacks_step/callbacks_step.jl")
 include("callbacks_stage/callbacks_stage.jl")
 include("semidiscretization/semidiscretization_euler_gravity.jl")
+include("time_integration/time_integration.jl")
 
 # `trixi_include` and special elixirs such as `convergence_test`
 include("auxiliary/special_elixirs.jl")
@@ -200,6 +201,8 @@ export boundary_condition_do_nothing,
        BoundaryConditionNavierStokesWall, NoSlip, Adiabatic, Isothermal,
        BoundaryConditionCoupled
 
+export coupling_converter_heaviside_2d
+
 export initial_condition_convergence_test, source_terms_convergence_test
 export source_terms_harmonic
 export initial_condition_poisson_nonperiodic, source_terms_poisson_nonperiodic,

src/coupling_converters/coupling_converters.jl

@@ -0,0 +1,13 @@
+# By default, Julia/LLVM does not use fused multiply-add operations (FMAs).
+# Since these FMAs can increase the performance of many numerical algorithms,
+# we need to opt-in explicitly.
+# See https://ranocha.de/blog/Optimizing_EC_Trixi for further details.
+@muladd begin
+#! format: noindent
+
+####################################################################################################
+# Include files with actual implementations for different systems of equations.
+
+include("coupling_converters_2d.jl")
+
+end # @muladd

src/coupling_converters/coupling_converters_2d.jl

@@ -0,0 +1,37 @@
+# By default, Julia/LLVM does not use fused multiply-add operations (FMAs).
+# Since these FMAs can increase the performance of many numerical algorithms,
+# we need to opt-in explicitly.
+# See https://ranocha.de/blog/Optimizing_EC_Trixi for further details.
+@muladd begin
+#! format: noindent
+
+
+@doc raw"""
+    Coupling converter function for a system of two LinearScalarAdvectionEquation2D.
+
+The coupling is given as a Heaviside step.
+```math
+c(x) = {c_0, for x \ge x_0 \times s
+        0, for x < x_0 \times s}
+```
+Here, `s` is the sign of the step function, x_0 the save_position
+of the step and c_0 the amplitude.
+"""
+function coupling_converter_heaviside_2d(x_0, c_0, s)
+    return (x, u) -> c_0 * (s*sign(x[2] - x_0) + 1.0)/2.0 * u
+end
+
+
+@doc raw"""
+    Coupling converter function for a system of two LinearScalarAdvectionEquation2D.
+
+The coupling is given as a linear function.
+```math
+c(x) = x * u(x)
+```
+"""
+function coupling_converter_linear_2d()
+    return (x, u) -> x[2]*u
+end
+
+end # @muladd