WIP: Enhanced p4est mesh view coupling

examples/Project.toml

@@ -0,0 +1,5 @@
+[deps]
+OrdinaryDiffEqLowStorageRK = "b0944070-b475-4768-8dec-fb6eb410534d"
+OrdinaryDiffEqSSPRK = "669c94d9-1f4b-4b64-b377-1aa079aa2388"
+Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80"
+Trixi = "a7f1ee26-1774-49b1-8366-f1abc58fbfcb"

examples/p4est_2d_dgsem/Project.toml

@@ -0,0 +1,9 @@
+[deps]
+AbbreviatedStackTraces = "ac637c84-cc71-43bf-9c33-c1b4316be3d4"
+BenchmarkTools = "6e4b80f9-dd63-53aa-95a3-0cdb28fa8baf"
+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,94 @@
+using OrdinaryDiffEqSSPRK, OrdinaryDiffEqLowStorageRK
+using Trixi
+
+###############################################################################
+# Simplest coupled setup consisting of two non-trivial mesh views.
+
+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)
+
+coordinates_min = (-1.0, -1.0) # minimum coordinates (min(x), min(y))
+coordinates_max = (1.0, 1.0) # maximum coordinates (max(x), max(y))
+
+trees_per_dimension = (8, 8)
+
+# Create parent P4estMesh with 8 x 8 trees and 8 x 8 elements
+# Since we couple through the boundaries, the periodicity does not matter here,
+# but it is to trigger parts of the code for the test.
+parent_mesh = P4estMesh(trees_per_dimension, polydeg = 3,
+                        coordinates_min = coordinates_min,
+                        coordinates_max = coordinates_max,
+                        initial_refinement_level = 0,
+                        periodicity = false)
+
+# Define the mesh views consisting of a small square in the center
+# and a square ring around it.
+cell_ids1 = vcat((1:18), (23:26), (31:34), (39:42), (47:64))
+mesh1 = P4estMeshView(parent_mesh, cell_ids1)
+cell_ids2 = vcat((19:22), (27:30), (35:38), (43:46))
+mesh2 = P4estMeshView(parent_mesh, cell_ids2)
+
+# Define a trivial coupling function.
+coupling_function = (x, u, equations_other, equations_own) -> u
+
+# Define a coupling function for each combination of interfaces.
+coupling_functions = Array{Function}(undef, 2, 2)
+coupling_functions[1, 1] = (x, u, equations_other, equations_own) -> u
+coupling_functions[1, 2] = (x, u, equations_other, equations_own) -> u
+coupling_functions[2, 1] = (x, u, equations_other, equations_own) -> u
+coupling_functions[2, 2] = (x, u, equations_other, equations_own) -> u
+
+# The mesh is coupled across the physical boundaries, which makes this setup
+# effectively double periodic.
+boundary_conditions = Dict(:x_neg => BoundaryConditionCoupledP4est(coupling_functions),
+                           :y_neg => BoundaryConditionCoupledP4est(coupling_functions),
+                           :y_pos => BoundaryConditionCoupledP4est(coupling_functions),
+                           :x_pos => BoundaryConditionCoupledP4est(coupling_functions))
+
+semi1 = SemidiscretizationHyperbolic(mesh1, equations, initial_condition_convergence_test,
+                                     solver,
+                                     boundary_conditions = boundary_conditions)
+semi2 = SemidiscretizationHyperbolic(mesh2, equations, initial_condition_convergence_test,
+                                     solver,
+                                     boundary_conditions = boundary_conditions)
+
+# Create a semidiscretization that bundles semi1 and semi2
+semi = SemidiscretizationCoupledP4est(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 = AnalysisCallbackCoupledP4est(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
+            ode_default_options()..., callback = callbacks);

examples/p4est_2d_dgsem/elixir_euler_mhd_coupled.jl

@@ -0,0 +1,126 @@
+using OrdinaryDiffEqSSPRK, OrdinaryDiffEqLowStorageRK
+using Trixi
+
+###############################################################################
+# One Euler system is coupled to one MHD system.
+
+# Pressure wave, same for the Euler system.
+function initial_condition_mhd(x, t, equations::IdealGlmMhdEquations2D)
+    rho = ((1.0 + 0.01 * sin(x[1] * 2 * pi)))
+    v1 = ((0.01 * sin((x[1] - 1 / 2) * 2 * pi)))
+    v2 = 0.0
+    v3 = 0.0
+    p = rho^equations.gamma
+    B1 = 0.0
+    B2 = 0.0
+    B3 = 0.0
+    psi = 0.0
+
+    return prim2cons(SVector(rho, v1, v2, v3, p, B1, B2, B3, psi), equations)
+end
+
+# Pressure wave, same as for the MHD system.
+function initial_condition_euler(x, t, equations::CompressibleEulerEquations2D)
+    rho = ((1.0 + 0.01 * sin(x[1] * 2 * pi)))
+    v1 = ((0.01 * sin((x[1] - 1 / 2) * 2 * pi)))
+    v2 = 0.0
+    p = rho .^ equations.gamma
+
+    return prim2cons(SVector(rho, v1, v2, p), equations)
+end
+
+# Define the parent mesh.
+coordinates_min = (-2.0, -2.0) # minimum coordinates (min(x), min(y))
+coordinates_max = (2.0, 2.0) # maximum coordinates (max(x), max(y))
+trees_per_dimension = (8, 8)
+# Here we set the priodicity to false for the coupling.
+# Since we couple through the physical boundaries the system is effectively periodic.
+parent_mesh = P4estMesh(trees_per_dimension, polydeg = 3,
+                        coordinates_min = coordinates_min,
+                        coordinates_max = coordinates_max,
+                        initial_refinement_level = 0, periodicity = (false, false))
+
+equations1 = IdealGlmMhdEquations2D(5 / 3)
+equations2 = CompressibleEulerEquations2D(5 / 3)
+
+# Define the coupling function between every possible pair of systems.
+coupling_functions = Array{Function}(undef, 2, 2)
+coupling_functions[1, 1] = (x, u, equations_other, equations_own) -> u
+coupling_functions[1, 2] = (x, u, equations_other, equations_own) -> SVector(u[1], u[2],
+                                                                             u[3], 0.0,
+                                                                             u[4], 0.0, 0.0,
+                                                                             0.0, 0.0)
+coupling_functions[2, 1] = (x, u, equations_other, equations_own) -> SVector(u[1], u[2],
+                                                                             u[3], u[5])
+coupling_functions[2, 2] = (x, u, equations_other, equations_own) -> u
+
+# semi 1 MHD.
+cell_ids1 = vcat(Vector(1:8), Vector(32:64))
+mesh1 = P4estMeshView(parent_mesh, cell_ids1)
+volume_flux = (flux_hindenlang_gassner, flux_nonconservative_powell)
+solver1 = DGSEM(polydeg = 3,
+                surface_flux = (flux_lax_friedrichs, flux_nonconservative_powell),
+                volume_integral = VolumeIntegralFluxDifferencing(volume_flux))
+boundary_conditions1 = Dict(:x_neg => BoundaryConditionCoupledP4est(coupling_functions),
+                            :y_neg => BoundaryConditionCoupledP4est(coupling_functions),
+                            :y_pos => BoundaryConditionCoupledP4est(coupling_functions),
+                            :x_pos => BoundaryConditionCoupledP4est(coupling_functions))
+semi1 = SemidiscretizationHyperbolic(mesh1, equations1, initial_condition_mhd, solver1,
+                                     boundary_conditions = boundary_conditions1)
+
+# semi 2 Euler
+cell_ids2 = Vector(9:31)
+mesh2 = P4estMeshView(parent_mesh, cell_ids2)
+solver2 = DGSEM(polydeg = 3, surface_flux = flux_hll,
+                volume_integral = VolumeIntegralWeakForm())
+boundary_conditions2 = Dict(:x_neg => BoundaryConditionCoupledP4est(coupling_functions),
+                            :y_neg => BoundaryConditionCoupledP4est(coupling_functions),
+                            :y_pos => BoundaryConditionCoupledP4est(coupling_functions),
+                            :x_pos => BoundaryConditionCoupledP4est(coupling_functions))
+semi2 = SemidiscretizationHyperbolic(mesh2, equations2, initial_condition_euler, solver2,
+                                     boundary_conditions = boundary_conditions2)
+
+# Create a semidiscretization that bundles semi1 and semi2
+semi = SemidiscretizationCoupledP4est(semi1, semi2)
+
+###############################################################################
+# ODE solvers, callbacks etc.
+
+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_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 = 0.8)
+
+glm_speed_callback = GlmSpeedCallback(glm_scale = 0.5, cfl = 0.8,
+                                      semi_indices = Vector([1]))
+
+# 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,
+                        glm_speed_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 = 0.0001, # 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/structured_2d_dgsem/Project.toml

@@ -0,0 +1,11 @@
+[deps]
+AbbreviatedStackTraces = "ac637c84-cc71-43bf-9c33-c1b4316be3d4"
+BenchmarkTools = "6e4b80f9-dd63-53aa-95a3-0cdb28fa8baf"
+Infiltrator = "5903a43b-9cc3-4c30-8d17-598619ec4e9b"
+OrdinaryDiffEq = "1dea7af3-3e70-54e6-95c3-0bf5283fa5ed"
+OrdinaryDiffEqLowStorageRK = "b0944070-b475-4768-8dec-fb6eb410534d"
+OrdinaryDiffEqSSPRK = "669c94d9-1f4b-4b64-b377-1aa079aa2388"
+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/structured_2d_dgsem/test.jl

@@ -0,0 +1,80 @@
+# #Ignore this
+# using Pkg
+# Pkg.activate("..")
+
+using Trixi
+using OrdinaryDiffEqSSPRK, OrdinaryDiffEqLowStorageRK
+using Plots
+
+equations = ShallowWaterEquations2D(gravity_constant = 9.81)
+# Define the initial condition
+
+function initial_condition_test(x, t, equations::ShallowWaterEquations2D)
+    x1, x2 = x
+    v1 = 0.0
+    v2 = 0.0
+    b = 0.0
+    H = 1.0
+    return prim2cons(SVector(H, v1, v2, b), equations)
+end
+
+initial_condition = initial_condition_test
+# Define the dirichlet boundary condition
+boundary_conditions = BoundaryConditionDirichlet(initial_condition_test)
+
+# Define the problem solver
+
+volume_flux = (flux_wintermeyer_etal, flux_nonconservative_wintermeyer_etal)
+solver = DGSEM(polydeg = 3,
+               surface_flux = (flux_lax_friedrichs, flux_nonconservative_fjordholm_etal),
+               volume_integral = VolumeIntegralFluxDifferencing(volume_flux))
+
+xl = 0.0;
+yl = 0.0;
+xr = 3.0;
+yr = 3.0;
+coordinates_min = (xl, yl)
+coordinates_max = (xr, yr)
+cells_per_dimension = (30, 30)
+
+# Define the mesh
+
+# If we use tree mesh, it is good.
+#mesh = TreeMesh(coordinates_min, coordinates_max, initial_refinement_level = 5, n_cells_max = 10_000, periodicity = false)         
+
+# If we use a structured mesh is wrong.
+mesh = StructuredMesh(cells_per_dimension, coordinates_min, coordinates_max,
+                      periodicity = false)
+
+# Create the semi discretization object
+semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver,
+                                    boundary_conditions = boundary_conditions)
+
+###############################################################################
+# ODE solver
+
+tspan = (0.0, 2.0)
+ode = semidiscretize(semi, tspan)
+
+###############################################################################
+# Callbacks
+
+summary_callback = SummaryCallback()
+
+analysis_interval = 1000
+analysis_callback = AnalysisCallback(semi, interval = analysis_interval,
+                                     extra_analysis_integrals = (lake_at_rest_error,))
+
+alive_callback = AliveCallback(analysis_interval = analysis_interval)
+
+stepsize_callback = StepsizeCallback(cfl = 0.5)
+
+callbacks = CallbackSet(summary_callback, analysis_callback, alive_callback,
+                        stepsize_callback)
+
+###############################################################################
+# run the simulation
+
+sol = solve(ode, CarpenterKennedy2N54(williamson_condition = false);
+            dt = 1.0, # solve needs some value here but it will be overwritten by the stepsize_callback
+            ode_default_options()..., callback = callbacks);

examples/structured_3d_dgsem/Project.toml

@@ -0,0 +1,10 @@
+[deps]
+AbbreviatedStackTraces = "ac637c84-cc71-43bf-9c33-c1b4316be3d4"
+BenchmarkTools = "6e4b80f9-dd63-53aa-95a3-0cdb28fa8baf"
+Infiltrator = "5903a43b-9cc3-4c30-8d17-598619ec4e9b"
+OrdinaryDiffEqLowStorageRK = "b0944070-b475-4768-8dec-fb6eb410534d"
+OrdinaryDiffEqSSPRK = "669c94d9-1f4b-4b64-b377-1aa079aa2388"
+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/t8code_2d_dgsem/Project.toml

@@ -0,0 +1,9 @@
+[deps]
+AbbreviatedStackTraces = "ac637c84-cc71-43bf-9c33-c1b4316be3d4"
+BenchmarkTools = "6e4b80f9-dd63-53aa-95a3-0cdb28fa8baf"
+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"

src/Trixi.jl

@@ -147,6 +147,7 @@ include("semidiscretization/semidiscretization_hyperbolic.jl")
 include("semidiscretization/semidiscretization_hyperbolic_parabolic.jl")
 include("semidiscretization/semidiscretization_euler_acoustics.jl")
 include("semidiscretization/semidiscretization_coupled.jl")
+include("semidiscretization/semidiscretization_coupled_p4est.jl")
 include("time_integration/time_integration.jl")
 include("callbacks_step/callbacks_step.jl")
 include("callbacks_stage/callbacks_stage.jl")
@@ -229,7 +230,7 @@ export boundary_condition_do_nothing,
        BoundaryConditionNavierStokesWall,
        NoSlip, Slip,
        Adiabatic, Isothermal,
-       BoundaryConditionCoupled
+       BoundaryConditionCoupled, BoundaryConditionCoupledP4est
 
 export initial_condition_convergence_test, source_terms_convergence_test,
        source_terms_lorentz, source_terms_collision_ion_electron,
@@ -291,14 +292,14 @@ export SemidiscretizationEulerGravity, ParametersEulerGravity,
        timestep_gravity_erk53_3Sstar!,
        timestep_gravity_carpenter_kennedy_erk54_2N!
 
-export SemidiscretizationCoupled
+export SemidiscretizationCoupled, SemidiscretizationCoupledP4est
 
 export SummaryCallback, SteadyStateCallback, AnalysisCallback, AliveCallback,
        SaveRestartCallback, SaveSolutionCallback, TimeSeriesCallback, VisualizationCallback,
        AveragingCallback,
        AMRCallback, StepsizeCallback,
        GlmSpeedCallback, LBMCollisionCallback, EulerAcousticsCouplingCallback,
-       TrivialCallback, AnalysisCallbackCoupled,
+       TrivialCallback, AnalysisCallbackCoupled, AnalysisCallbackCoupledP4est,
        AnalysisSurfaceIntegral, DragCoefficientPressure2D, LiftCoefficientPressure2D,
        DragCoefficientShearStress2D, LiftCoefficientShearStress2D,
        DragCoefficientPressure3D, LiftCoefficientPressure3D