WIP: Sc/p4est view coupled

examples/p4est_2d_dgsem/elixir_advection_coupled.jl

@@ -0,0 +1,121 @@
+# The same setup as tree_2d_dgsem/elixir_advection_basic.jl
+# to verify the StructuredMesh implementation against TreeMesh
+
+using OrdinaryDiffEq
+using Trixi
+using Trixi2Vtk
+
+
+function initial_constant(x, t, equations)
+    return SVector(2.5,)
+end
+
+function initial_sinx(x, t, equations)
+    return SVector(sin(x[1]*pi),)
+end
+
+function initial_siny(x, t, equations)
+    return SVector(sin(x[2]*pi),)
+end
+
+function initial_linx(x, t, equations)
+    return SVector(x[1],)
+end
+
+function initial_liny(x, t, equations)
+    return SVector(x[2],)
+end
+
+function gaussian(x, t, equations)
+    if x[1] < 0
+        return SVector(exp((-(x[1]+0.5)^2 - x[2]^2)*15),)
+    else
+        return SVector(0.0,)
+    end
+end
+
+initial_condition = gaussian
+
+###############################################################################
+# semidiscretization of the linear advection equation
+
+advection_velocity = (1.0, 0.0)
+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)
+#                volume_integral=VolumeIntegralUpwind(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 P4estMesh with 8 x 8 trees and 16 x 16 elements
+parent_mesh = P4estMesh(trees_per_dimension, polydeg=3,
+                        coordinates_min=coordinates_min, coordinates_max=coordinates_max,
+                        initial_refinement_level=1, periodicity=(false, true))
+mesh1 = P4estMeshView(parent_mesh; indices_min = (1, 1), indices_max = (6, 8),
+                      coordinates_min = coordinates_min, coordinates_max = (0.5, 1.0),
+                      periodicity = (false, true))
+mesh2 = P4estMeshView(parent_mesh; indices_min = (7, 1), indices_max = (8, 8),
+                      coordinates_min = (0.5, -1.0), coordinates_max = coordinates_max,
+                      periodicity = (false, true))
+
+coupling_function1 = (x, u, equations_other, equations_own) -> u
+coupling_function2 = (x, u, equations_other, equations_own) -> u
+
+boundary_conditions1 = Dict(:x_neg => BoundaryConditionCoupled(2, (:end, :i_forward), Float64, coupling_function1),
+                            :x_pos => BoundaryConditionCoupled(2, (:begin, :i_forward), Float64, coupling_function1))
+boundary_conditions2 = Dict(:x_neg => BoundaryConditionCoupled(1, (:end, :i_forward), Float64, coupling_function2),
+                            :x_pos => BoundaryConditionCoupled(1, (:begin, :i_forward), Float64, coupling_function2))
+
+# A semidiscretization collects data structures and functions for the spatial discretization
+semi1 = SemidiscretizationHyperbolic(mesh1, equations, initial_condition, solver, boundary_conditions=boundary_conditions1)
+semi2 = SemidiscretizationHyperbolic(mesh2, equations, initial_condition, solver, boundary_conditions=boundary_conditions2)
+
+# Create a semidiscretization that bundles semi1 and semi2
+semi = SemidiscretizationCoupled(semi1, semi2)
+
+###############################################################################
+# ODE solvers, callbacks etc.
+
+# ode = semidiscretize(semi, (0.0, 20.3));
+ode = semidiscretize(semi, (0.0, 5.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()
+
+# Convert the snapshots into vtk format.
+# trixi2vtk("out/solution_*.h5")

examples/p4est_2d_dgsem/elixir_advection_p4estview.jl

@@ -0,0 +1,80 @@
+# The same setup as tree_2d_dgsem/elixir_advection_basic.jl
+# to verify the StructuredMesh implementation against TreeMesh
+
+using OrdinaryDiffEq
+using Trixi
+
+###############################################################################
+# semidiscretization of the linear advection equation
+
+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 P4estMesh with 8 x 8 trees and 16 x 16 elements and the mesh views
+parent_mesh = P4estMesh(trees_per_dimension, polydeg = 3,
+                        coordinates_min = coordinates_min,
+                        coordinates_max = coordinates_max,
+                        initial_refinement_level = 1, periodicity = (false, true))
+mesh1 = P4estMeshView(parent_mesh; indices_min = (1, 1), indices_max = (4, 8),
+                      coordinates_min = coordinates_min, coordinates_max = (0.0, 1.0),
+                      periodicity = (false, true))
+mesh2 = P4estMeshView(parent_mesh; indices_min = (5, 1), indices_max = (8, 8),
+                      coordinates_min = (0.0, -1.0), coordinates_max = coordinates_max,
+                      periodicity = (false, true))
+
+boundary_conditions = Dict(:x_neg => BoundaryConditionDirichlet(initial_condition_convergence_test),
+                           :x_pos => BoundaryConditionDirichlet(initial_condition_convergence_test))
+
+# A semidiscretization collects data structures and functions for the spatial discretization
+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 = SemidiscretizationCoupled(semi1, semi2)
+
+###############################################################################
+# 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_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()

src/Trixi.jl

@@ -227,7 +227,7 @@ export lake_at_rest_error
 export ncomponents, eachcomponent
 
 export TreeMesh, StructuredMesh, StructuredMeshView, UnstructuredMesh2D, P4estMesh,
-       T8codeMesh
+       P4estMeshView, T8codeMesh
 
 export DG,
        DGSEM, LobattoLegendreBasis,

src/callbacks_step/analysis_dg2d.jl

@@ -32,7 +32,8 @@ end
 function create_cache_analysis(analyzer,
                                mesh::Union{StructuredMesh{2}, StructuredMeshView{2},
                                            UnstructuredMesh2D,
-                                           P4estMesh{2}, T8codeMesh{2}},
+                                           P4estMesh{2}, P4estMeshView{2},
+                                           T8codeMesh{2}},
                                equations, dg::DG, cache,
                                RealT, uEltype)
 
@@ -109,7 +110,8 @@ end
 
 function calc_error_norms(func, u, t, analyzer,
                           mesh::Union{StructuredMesh{2}, StructuredMeshView{2},
-                                      UnstructuredMesh2D, P4estMesh{2}, T8codeMesh{2}},
+                                      UnstructuredMesh2D, P4estMesh{2},
+                                      P4estMeshView{2}, T8codeMesh{2}},
                           equations,
                           initial_condition, dg::DGSEM, cache, cache_analysis)
     @unpack vandermonde, weights = analyzer
@@ -179,6 +181,7 @@ end
 function integrate_via_indices(func::Func, u,
                                mesh::Union{StructuredMesh{2}, StructuredMeshView{2},
                                            UnstructuredMesh2D, P4estMesh{2},
+                                           P4estMeshView{2},
                                            T8codeMesh{2}},
                                equations,
                                dg::DGSEM, cache, args...; normalize = true) where {Func}
@@ -208,7 +211,8 @@ end
 
 function integrate(func::Func, u,
                    mesh::Union{TreeMesh{2}, StructuredMesh{2}, StructuredMeshView{2},
-                               UnstructuredMesh2D, P4estMesh{2}, T8codeMesh{2}},
+                               UnstructuredMesh2D, P4estMesh{2}, P4estMeshView{2},
+                               T8codeMesh{2}},
                    equations, dg::DG, cache; normalize = true) where {Func}
     integrate_via_indices(u, mesh, equations, dg, cache;
                           normalize = normalize) do u, i, j, element, equations, dg
@@ -218,7 +222,7 @@ function integrate(func::Func, u,
 end
 
 function integrate(func::Func, u,
-                   mesh::Union{TreeMesh{2}, P4estMesh{2}},
+                   mesh::Union{TreeMesh{2}, P4estMesh{2}, P4estMeshView{2}},
                    equations, equations_parabolic,
                    dg::DGSEM,
                    cache, cache_parabolic; normalize = true) where {Func}
@@ -237,7 +241,8 @@ end
 
 function analyze(::typeof(entropy_timederivative), du, u, t,
                  mesh::Union{TreeMesh{2}, StructuredMesh{2}, StructuredMeshView{2},
-                             UnstructuredMesh2D, P4estMesh{2}, T8codeMesh{2}},
+                             UnstructuredMesh2D, P4estMesh{2}, P4estMeshView{2},
+                             T8codeMesh{2}},
                  equations, dg::DG, cache)
     # Calculate ∫(∂S/∂u ⋅ ∂u/∂t)dΩ
     integrate_via_indices(u, mesh, equations, dg, cache,
@@ -282,7 +287,7 @@ end
 
 function analyze(::Val{:l2_divb}, du, u, t,
                  mesh::Union{StructuredMesh{2}, UnstructuredMesh2D, P4estMesh{2},
-                             T8codeMesh{2}},
+                             P4estMeshView{2}, T8codeMesh{2}},
                  equations::IdealGlmMhdEquations2D, dg::DGSEM, cache)
     @unpack contravariant_vectors = cache.elements
     integrate_via_indices(u, mesh, equations, dg, cache, cache,
@@ -350,7 +355,7 @@ end
 
 function analyze(::Val{:linf_divb}, du, u, t,
                  mesh::Union{StructuredMesh{2}, UnstructuredMesh2D, P4estMesh{2},
-                             T8codeMesh{2}},
+                             P4estMeshView{2}, T8codeMesh{2}},
                  equations::IdealGlmMhdEquations2D, dg::DGSEM, cache)
     @unpack derivative_matrix, weights = dg.basis
     @unpack contravariant_vectors = cache.elements

src/callbacks_step/save_solution_dg.jl

@@ -9,7 +9,7 @@ function save_solution_file(u, time, dt, timestep,
                             mesh::Union{SerialTreeMesh, StructuredMesh,
                                         StructuredMeshView,
                                         UnstructuredMesh2D, SerialP4estMesh,
-                                        SerialT8codeMesh},
+                                        P4estMeshView, SerialT8codeMesh},
                             equations, dg::DG, cache,
                             solution_callback,
                             element_variables = Dict{Symbol, Any}(),

src/callbacks_step/stepsize_dg2d.jl

@@ -77,7 +77,7 @@ end
 
 function max_dt(u, t,
                 mesh::Union{StructuredMesh{2}, UnstructuredMesh2D, P4estMesh{2},
-                            T8codeMesh{2}, StructuredMeshView{2}},
+                            T8codeMesh{2}, StructuredMeshView{2}, P4estMeshView{2}},
                 constant_speed::False, equations, dg::DG, cache)
     # to avoid a division by zero if the speed vanishes everywhere,
     # e.g. for steady-state linear advection
@@ -113,7 +113,7 @@ end
 
 function max_dt(u, t,
                 mesh::Union{StructuredMesh{2}, UnstructuredMesh2D, P4estMesh{2},
-                            T8codeMesh{2}, StructuredMeshView{2}},
+                            T8codeMesh{2}, StructuredMeshView{2}, P4estMeshView{2}},
                 constant_speed::True, equations, dg::DG, cache)
     @unpack contravariant_vectors, inverse_jacobian = cache.elements
 

src/meshes/meshes.jl

@@ -13,6 +13,7 @@ include("unstructured_mesh.jl")
 include("face_interpolant.jl")
 include("transfinite_mappings_3d.jl")
 include("p4est_mesh.jl")
+include("p4est_mesh_view.jl")
 include("t8code_mesh.jl")
 include("mesh_io.jl")
 include("dgmulti_meshes.jl")

src/meshes/p4est_mesh.jl

@@ -88,6 +88,7 @@ end
 
 @inline Base.ndims(::P4estMesh{NDIMS}) where {NDIMS} = NDIMS
 @inline Base.real(::P4estMesh{NDIMS, RealT}) where {NDIMS, RealT} = RealT
+@inline Base.size(mesh::P4estMesh) = size(mesh.tree_node_coordinates)[end]
 
 @inline function ntrees(mesh::P4estMesh)
     return mesh.p4est.trees.elem_count[]
@@ -2051,4 +2052,9 @@ function collect_new_cells(mesh::P4estMesh)
 
     return new_cells
 end
+
+Base.size(mesh::P4estMesh) = return (16, 16)
+Base.size(mesh::P4estMesh, i::Int) = return 16
+Base.axes(mesh::P4estMesh) = map(Base.OneTo, size(mesh))
+Base.axes(mesh::P4estMesh, i) = Base.OneTo(size(mesh, i))
 end # @muladd