From b964c93c1969615fda7e7d68574b3ebd93634aea Mon Sep 17 00:00:00 2001 From: vincmarks Date: Thu, 17 Jul 2025 14:00:43 +0200 Subject: [PATCH 01/45] new paragraph: Boundary conditions --- docs/src/meshes/p4est_mesh.md | 3 +++ 1 file changed, 3 insertions(+) diff --git a/docs/src/meshes/p4est_mesh.md b/docs/src/meshes/p4est_mesh.md index 57bb06a7b4d..dc7c8082b26 100644 --- a/docs/src/meshes/p4est_mesh.md +++ b/docs/src/meshes/p4est_mesh.md @@ -939,3 +939,6 @@ rebalance_solver!(u0, mesh, equations, dg, cache, old_global_first_quadrant) reinitialize_boundaries!(semi.boundary_conditions, cache) # Needs to be called after `rebalance_solver!` ``` This code could then be placed in the [`resize!`](https://github.com/trixi-framework/Trixi.jl/blob/eaeb04113523500ed831e3ab459694f12f7a49ea/src/time_integration/methods_2N.jl#L251-L255) function of a corresponding multirate integrator to ensure load-balancing for simulations involving AMR. + +### Boundary conditions +For p4est meshes, boundary conditions are defined and stored in dictionaries ( as shown for example in `examples/p4est_2d_dgsem/elixir_advection_diffusion_nonperiodic_amr.jl`). If you’d like to apply the same condition to every face of the mesh, you can use the convenient functions `boundary_condition_default_p4est_2D` and `boundary_condition_default_p4est_3D`. \ No newline at end of file From d9790da8522bf1fb9fc0d28dfb46a4587766b6b9 Mon Sep 17 00:00:00 2001 From: vincmarks Date: Thu, 17 Jul 2025 14:01:03 +0200 Subject: [PATCH 02/45] new paragraph: Boundary conditions --- docs/src/meshes/structured_mesh.md | 3 +++ 1 file changed, 3 insertions(+) diff --git a/docs/src/meshes/structured_mesh.md b/docs/src/meshes/structured_mesh.md index 6db8bf66b13..f0ee805395d 100644 --- a/docs/src/meshes/structured_mesh.md +++ b/docs/src/meshes/structured_mesh.md @@ -11,3 +11,6 @@ such as the compressible Euler equations can use [`FluxRotated`](@ref) to wrap numerical fluxes implemented only for Cartesian meshes. This simplifies the re-use of existing functionality for the [`TreeMesh`](@ref) but is usually less efficient, cf. [PR #550](https://github.com/trixi-framework/Trixi.jl/pull/550). + +### Boundary conditions +For structured meshes, boundary conditions are defined and stored in named tuples ( as shown for example in `examples/structured_1d_dgsem/elixir_euler_source_terms_nonperiodic.jl`). If you’d like to apply the same condition to every face of the mesh, you can use the convenient functions `boundary_condition_default_structured_1D`, `boundary_condition_default_structured_2D` and `boundary_condition_default_structured_3D`. From c1ceea3cc60b708e8c45bb352f4fa7efa0cf7098 Mon Sep 17 00:00:00 2001 From: vincmarks Date: Thu, 17 Jul 2025 14:01:42 +0200 Subject: [PATCH 03/45] additional functions for default choice of boundary conditions --- .../sort_boundary_conditions.jl | 126 +++++++++++++++++- 1 file changed, 125 insertions(+), 1 deletion(-) diff --git a/src/solvers/dgsem_unstructured/sort_boundary_conditions.jl b/src/solvers/dgsem_unstructured/sort_boundary_conditions.jl index 0cb3bd7f409..971aaa7e331 100644 --- a/src/solvers/dgsem_unstructured/sort_boundary_conditions.jl +++ b/src/solvers/dgsem_unstructured/sort_boundary_conditions.jl @@ -119,4 +119,128 @@ function initialize!(boundary_types_container::UnstructuredSortedBoundaryTypes{N return boundary_types_container end -end # @muladd + + +""" + boundary_condition_default_2D + +Create a default boundary condition dictionary for p4est meshes in 2D +that use the standard boundary naming convention. + +This function applies the same boundary condition to all standard boundaries: +- `:x_neg`: negative x-direction boundary +- `:x_pos`: positive x-direction boundary +- `:y_neg`: negative y-direction boundary +- `:y_pos`: positive y-direction boundary + +# Arguments +- `boundary_condition`: The boundary condition function to apply to all boundaries + +# Returns +- `Dict{Symbol, Any}`: Dictionary mapping boundary names to the boundary condition + +""" +function boundary_condition_default_p4est_2D(boundarycondition) + + return Dict(:x_neg => boundarycondition, + :y_neg => boundarycondition, + :y_pos => boundarycondition, + :x_pos => boundarycondition) +end + +""" + boundary_condition_default_3D +Create a default boundary condition dictionary for p4est meshes in 3D +that use the standard boundary naming convention. +This function applies the same boundary condition to all standard boundaries: +- `:x_neg`: negative x-direction boundary +- `:x_pos`: positive x-direction boundary +- `:y_neg`: negative y-direction boundary +- `:y_pos`: positive y-direction boundary +- `:z_neg`: negative z-direction boundary +- `:z_pos`: positive z-direction boundary +# Arguments +- `boundary_condition`: The boundary condition function to apply to all boundaries +# Returns +- `Dict{Symbol, Any}`: Dictionary mapping boundary names to the boundary condition + +""" + +function boundary_condition_default_p4est_3D(boundarycondition) + + return Dict(:x_neg => boundarycondition, + :x_pos => boundarycondition, + :y_neg => boundarycondition, + :y_pos => boundarycondition, + :z_neg => boundarycondition, + :z_pos => boundarycondition) +end + + +""" + boundary_condition_default_structured_1D +Create a default boundary condition dictionary for structured meshes in 1D +that use the standard boundary naming convention. +This function applies the same boundary condition to all standard boundaries: +- `:x_neg`: negative x-direction boundary +- `:x_pos`: positive x-direction boundary +# Arguments +- `boundarycondition`: The boundary condition function to apply to all boundaries +# Returns +- Named tuple mapping boundary names to the boundary condition +""" + +function boundary_condition_default_structured_1D(boundarycondition) + + return (x_neg = boundarycondition, + x_pos = boundarycondition) +end +""" + boundary_condition_default_structured_2D +Create a default boundary condition dictionary for structured meshes in 2D +that use the standard boundary naming convention. +This function applies the same boundary condition to all standard boundaries: +- `:x_neg`: negative x-direction boundary +- `:x_pos`: positive x-direction boundary +- `:y_neg`: negative y-direction boundary +- `:y_pos`: positive y-direction boundary +# Arguments +- `boundarycondition`: The boundary condition function to apply to all boundaries +# Returns +- Named tuple mapping boundary names to the boundary condition +""" +function boundary_condition_default_structured_2D(boundarycondition) + + return (x_neg = boundarycondition, + x_pos = boundarycondition, + y_neg = boundarycondition, + y_pos = boundarycondition) +end +""" + boundary_condition_default_structured_3D +Create a default boundary condition dictionary for structured meshes in 3D +that use the standard boundary naming convention. +This function applies the same boundary condition to all standard boundaries: +- `:x_neg`: negative x-direction boundary +- `:x_pos`: positive x-direction boundary +- `:y_neg`: negative y-direction boundary +- `:y_pos`: positive y-direction boundary +- `:z_neg`: negative z-direction boundary +- `:z_pos`: positive z-direction boundary +# Arguments +- `boundarycondition`: The boundary condition function to apply to all boundaries +# Returns +- Named tuple mapping boundary names to the boundary condition +""" +function boundary_condition_default_structured_3D(boundarycondition) + + return (x_neg = boundarycondition, + x_pos = boundarycondition, + y_neg = boundarycondition, + y_pos = boundarycondition, + z_neg = boundarycondition, + z_pos = boundarycondition) + +end + +end#@muladd From ffc3e41aa2a8b07c02bc61a9493aee4ef4eab8af Mon Sep 17 00:00:00 2001 From: vincmarks Date: Thu, 17 Jul 2025 14:02:00 +0200 Subject: [PATCH 04/45] added the default condition functions in Trixi --- src/Trixi.jl | 5 +++++ 1 file changed, 5 insertions(+) diff --git a/src/Trixi.jl b/src/Trixi.jl index a707437655e..5aafef4db18 100644 --- a/src/Trixi.jl +++ b/src/Trixi.jl @@ -217,6 +217,11 @@ export initial_condition_constant, initial_condition_weak_blast_wave export boundary_condition_do_nothing, + boundary_condition_default_p4est_2D, + boundary_condition_default_p4est_3D, + boundary_condition_default_structured_1D, + boundary_condition_default_structured_2D, + boundary_condition_default_structured_3D, boundary_condition_periodic, BoundaryConditionDirichlet, BoundaryConditionNeumann, From 7737998bb5aab86789cedc0db148f2cfd83ce4b1 Mon Sep 17 00:00:00 2001 From: vincmarks Date: Thu, 7 Aug 2025 17:14:17 +0200 Subject: [PATCH 05/45] added the tree mesh default boundary condition --- .../sort_boundary_conditions.jl | 136 +++++++++++++----- 1 file changed, 102 insertions(+), 34 deletions(-) diff --git a/src/solvers/dgsem_unstructured/sort_boundary_conditions.jl b/src/solvers/dgsem_unstructured/sort_boundary_conditions.jl index 971aaa7e331..220114be604 100644 --- a/src/solvers/dgsem_unstructured/sort_boundary_conditions.jl +++ b/src/solvers/dgsem_unstructured/sort_boundary_conditions.jl @@ -122,13 +122,13 @@ end """ - boundary_condition_default_2D + boundary_condition_default_p4est_2D Create a default boundary condition dictionary for p4est meshes in 2D that use the standard boundary naming convention. This function applies the same boundary condition to all standard boundaries: -- `:x_neg`: negative x-direction boundary +- `:x_neg`: negative x-direction boundary - `:x_pos`: positive x-direction boundary - `:y_neg`: negative y-direction boundary - `:y_pos`: positive y-direction boundary @@ -140,12 +140,12 @@ This function applies the same boundary condition to all standard boundaries: - `Dict{Symbol, Any}`: Dictionary mapping boundary names to the boundary condition """ -function boundary_condition_default_p4est_2D(boundarycondition) - - return Dict(:x_neg => boundarycondition, - :y_neg => boundarycondition, - :y_pos => boundarycondition, - :x_pos => boundarycondition) +function boundary_condition_default_p4est_2D(boundary_condition) + + return Dict(:x_neg => boundary_condition, + :y_neg => boundary_condition, + :y_pos => boundary_condition, + :x_pos => boundary_condition) end """ @@ -166,14 +166,14 @@ This function applies the same boundary condition to all standard boundaries: """ -function boundary_condition_default_p4est_3D(boundarycondition) - - return Dict(:x_neg => boundarycondition, - :x_pos => boundarycondition, - :y_neg => boundarycondition, - :y_pos => boundarycondition, - :z_neg => boundarycondition, - :z_pos => boundarycondition) +function boundary_condition_default_p4est_3D(boundary_condition) + + return Dict(:x_neg => boundary_condition, + :x_pos => boundary_condition, + :y_neg => boundary_condition, + :y_pos => boundary_condition, + :z_neg => boundary_condition, + :z_pos => boundary_condition) end @@ -190,10 +190,10 @@ This function applies the same boundary condition to all standard boundaries: - Named tuple mapping boundary names to the boundary condition """ -function boundary_condition_default_structured_1D(boundarycondition) - - return (x_neg = boundarycondition, - x_pos = boundarycondition) +function boundary_condition_default_structured_1D(boundary_condition) + + return (x_neg = boundary_condition, + x_pos = boundary_condition) end """ boundary_condition_default_structured_2D @@ -209,12 +209,12 @@ This function applies the same boundary condition to all standard boundaries: # Returns - Named tuple mapping boundary names to the boundary condition """ -function boundary_condition_default_structured_2D(boundarycondition) - - return (x_neg = boundarycondition, - x_pos = boundarycondition, - y_neg = boundarycondition, - y_pos = boundarycondition) +function boundary_condition_default_structured_2D(boundary_condition) + + return (x_neg = boundary_condition, + x_pos = boundary_condition, + y_neg = boundary_condition, + y_pos = boundary_condition) end """ boundary_condition_default_structured_3D @@ -232,15 +232,83 @@ This function applies the same boundary condition to all standard boundaries: # Returns - Named tuple mapping boundary names to the boundary condition """ -function boundary_condition_default_structured_3D(boundarycondition) - - return (x_neg = boundarycondition, - x_pos = boundarycondition, - y_neg = boundarycondition, - y_pos = boundarycondition, - z_neg = boundarycondition, - z_pos = boundarycondition) +function boundary_condition_default_structured_3D(boundary_condition) + + return (x_neg = boundary_condition, + x_pos = boundary_condition, + y_neg = boundary_condition, + y_pos = boundary_condition, + z_neg = boundary_condition, + z_pos = boundary_condition) end + +""" + boundary_condition_default_tree_1D +Create a default boundary condition dictionary for tree meshes in 1D +that use the standard boundary naming convention. +This function applies the same boundary condition to all standard boundaries: +- `:x_neg`: negative x-direction boundary +- `:x_pos`: positive x-direction boundary +# Arguments +- `boundary_condition`: The boundary condition function to apply to all boundaries +# Returns +- Named tuple mapping boundary names to the boundary condition + +""" +function boundary_condition_default_tree_1D(boundary_condition) + + return (x_neg = boundary_condition, + x_pos = boundary_condition) +end + +""" + boundary_condition_default_tree_2D +Create a default boundary condition dictionary for tree meshes in 2D +that use the standard boundary naming convention. +This function applies the same boundary condition to all standard boundaries: +- `:x_neg`: negative x-direction boundary +- `:x_pos`: positive x-direction boundary +- `:y_neg`: negative y-direction boundary +- `:y_pos`: positive y-direction boundary +# Arguments +- `boundary_condition`: The boundary condition function to apply to all boundaries +# Returns +- Named tuple mapping boundary names to the boundary condition +""" +function boundary_condition_default_tree_2D(boundary_condition) + + return (x_neg = boundary_condition, + x_pos = boundary_condition, + y_neg = boundary_condition, + y_pos = boundary_condition) +end + +""" + boundary_condition_default_tree_3D +Create a default boundary condition dictionary for tree meshes in 3D +that use the standard boundary naming convention. +This function applies the same boundary condition to all standard boundaries: +- `:x_neg`: negative x-direction boundary +- `:x_pos`: positive x-direction boundary +- `:y_neg`: negative y-direction boundary +- `:y_pos`: positive y-direction boundary +- `:z_neg`: negative z-direction boundary +- `:z_pos`: positive z-direction boundary +# Arguments +- `boundary_condition`: The boundary condition function to apply to all boundaries +# Returns +- Named tuple mapping boundary names to the boundary condition +""" + +function boundary_condition_default_tree_3D(boundary_condition) + + return (x_neg = boundary_condition, + x_pos = boundary_condition, + y_neg = boundary_condition, + y_pos = boundary_condition, + z_neg = boundary_condition, + z_pos = boundary_condition) +end end#@muladd From dc011bf167528d1f35245ce8151139abb439cae1 Mon Sep 17 00:00:00 2001 From: vincmarks Date: Thu, 7 Aug 2025 17:15:23 +0200 Subject: [PATCH 06/45] add boundary condition paragraph --- docs/src/meshes/tree_mesh.md | 10 ++++++++++ 1 file changed, 10 insertions(+) diff --git a/docs/src/meshes/tree_mesh.md b/docs/src/meshes/tree_mesh.md index 9a88919a40c..26dfb4784f0 100644 --- a/docs/src/meshes/tree_mesh.md +++ b/docs/src/meshes/tree_mesh.md @@ -9,3 +9,13 @@ It is limited to hypercube domains (that is, lines in 1D, squares in 2D and cube Due to its Cartesian nature, (numerical) fluxes need to implement methods dispatching on the `orientation::Integer` as described in the [conventions](@ref conventions). + + +### Boundary conditions +For `Tree meshes`, boundary conditions are defined and stored in named tuples (as shown for example in `examples/tree_2d_dgsem/elixir_advection_diffusion_nonperiodic.jl`). If you’d like to apply the same condition to every face of the mesh, you can use the convenient `functions boundary_condition_default_tree_1D`, `boundary_condition_default_tree_2D` and `boundary_condition_default_tree_3D`. For example, in the two dimensional case: + +```julia +initial_condition = initial_condition_eriksson_johnson + +boundary_conditions = boundary_condition_default_tree_2D(BoundaryConditionDirichlet(initial_condition)) +``` \ No newline at end of file From 284972e480f11afa2e48e6c87d39102b3e33137a Mon Sep 17 00:00:00 2001 From: vincmarks Date: Thu, 7 Aug 2025 17:15:52 +0200 Subject: [PATCH 07/45] some reformulaition --- docs/src/meshes/p4est_mesh.md | 7 ++++++- docs/src/meshes/structured_mesh.md | 7 ++++++- 2 files changed, 12 insertions(+), 2 deletions(-) diff --git a/docs/src/meshes/p4est_mesh.md b/docs/src/meshes/p4est_mesh.md index dc7c8082b26..a420ca517ce 100644 --- a/docs/src/meshes/p4est_mesh.md +++ b/docs/src/meshes/p4est_mesh.md @@ -941,4 +941,9 @@ reinitialize_boundaries!(semi.boundary_conditions, cache) # Needs to be called a This code could then be placed in the [`resize!`](https://github.com/trixi-framework/Trixi.jl/blob/eaeb04113523500ed831e3ab459694f12f7a49ea/src/time_integration/methods_2N.jl#L251-L255) function of a corresponding multirate integrator to ensure load-balancing for simulations involving AMR. ### Boundary conditions -For p4est meshes, boundary conditions are defined and stored in dictionaries ( as shown for example in `examples/p4est_2d_dgsem/elixir_advection_diffusion_nonperiodic_amr.jl`). If you’d like to apply the same condition to every face of the mesh, you can use the convenient functions `boundary_condition_default_p4est_2D` and `boundary_condition_default_p4est_3D`. \ No newline at end of file +For `P4est meshes`, boundary conditions are defined and stored in dictionaries (as shown for example in `examples/p4est_2d_dgsem/elixir_advection_diffusion_nonperiodic_amr.jl`). If you’d like to apply the same condition to every face of the mesh, you can use the convenient functions `boundary_condition_default_p4est_2D` and `boundary_condition_default_p4est_3D`. For example, in the two dimensional case: + +```julia +boundary_condition = boundary_condition_slip_wall +boundary_conditions = boundary_condition_default_p4est_2D(boundary_condition) +``` diff --git a/docs/src/meshes/structured_mesh.md b/docs/src/meshes/structured_mesh.md index f0ee805395d..14f1ca04556 100644 --- a/docs/src/meshes/structured_mesh.md +++ b/docs/src/meshes/structured_mesh.md @@ -13,4 +13,9 @@ the re-use of existing functionality for the [`TreeMesh`](@ref) but is usually less efficient, cf. [PR #550](https://github.com/trixi-framework/Trixi.jl/pull/550). ### Boundary conditions -For structured meshes, boundary conditions are defined and stored in named tuples ( as shown for example in `examples/structured_1d_dgsem/elixir_euler_source_terms_nonperiodic.jl`). If you’d like to apply the same condition to every face of the mesh, you can use the convenient functions `boundary_condition_default_structured_1D`, `boundary_condition_default_structured_2D` and `boundary_condition_default_structured_3D`. +For `Structured meshes`, boundary conditions are defined and stored in named tuples (as shown for example in `examples/structured_1d_dgsem/elixir_euler_source_terms_nonperiodic.jl`). If you’d like to apply the same condition to every face of the mesh, you can use the convenient functions `boundary_condition_default_structured_1D`, `boundary_condition_default_structured_2D` and `boundary_condition_default_structured_3D`. For example, in the one dimensional case: + +```julia +boundary_condition = BoundaryConditionDirichlet(initial_condition) +boundary_conditions = boundary_condition_default_structured_1D(boundary_condition) +``` From 752aeedd187a16af55d95db7c26c522995366963 Mon Sep 17 00:00:00 2001 From: vincmarks Date: Thu, 7 Aug 2025 17:27:50 +0200 Subject: [PATCH 08/45] added the boundary_condition_default_tree funtions to trixi --- src/Trixi.jl | 3 +++ 1 file changed, 3 insertions(+) diff --git a/src/Trixi.jl b/src/Trixi.jl index 7593631681c..c62546e422d 100644 --- a/src/Trixi.jl +++ b/src/Trixi.jl @@ -225,6 +225,9 @@ export boundary_condition_do_nothing, boundary_condition_default_structured_1D, boundary_condition_default_structured_2D, boundary_condition_default_structured_3D, + boundary_condition_default_tree_1D, + boundary_condition_default_tree_2D, + boundary_condition_default_tree_3D, boundary_condition_periodic, BoundaryConditionDirichlet, BoundaryConditionNeumann, From a4378aa8292c08be824b4305bacc078e0c4e196d Mon Sep 17 00:00:00 2001 From: Vincent Marks Date: Thu, 25 Sep 2025 14:45:48 +0200 Subject: [PATCH 09/45] Update src/solvers/dgsem_unstructured/sort_boundary_conditions.jl Co-authored-by: Joshua Lampert <51029046+JoshuaLampert@users.noreply.github.com> --- src/solvers/dgsem_unstructured/sort_boundary_conditions.jl | 1 - 1 file changed, 1 deletion(-) diff --git a/src/solvers/dgsem_unstructured/sort_boundary_conditions.jl b/src/solvers/dgsem_unstructured/sort_boundary_conditions.jl index 7ac5be33101..d01c71be624 100644 --- a/src/solvers/dgsem_unstructured/sort_boundary_conditions.jl +++ b/src/solvers/dgsem_unstructured/sort_boundary_conditions.jl @@ -143,7 +143,6 @@ This function applies the same boundary condition to all standard boundaries: # Returns - `Dict{Symbol, Any}`: Dictionary mapping boundary names to the boundary condition - """ function boundary_condition_default_p4est_2D(boundary_condition) From 1331679f38f45a28f840bb3a8f5418b78d70a410 Mon Sep 17 00:00:00 2001 From: Vincent Marks Date: Thu, 25 Sep 2025 14:46:26 +0200 Subject: [PATCH 10/45] Update docs/src/meshes/tree_mesh.md Co-authored-by: Joshua Lampert <51029046+JoshuaLampert@users.noreply.github.com> --- docs/src/meshes/tree_mesh.md | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/docs/src/meshes/tree_mesh.md b/docs/src/meshes/tree_mesh.md index 26dfb4784f0..fe30e08332d 100644 --- a/docs/src/meshes/tree_mesh.md +++ b/docs/src/meshes/tree_mesh.md @@ -12,7 +12,7 @@ dispatching on the `orientation::Integer` as described in the ### Boundary conditions -For `Tree meshes`, boundary conditions are defined and stored in named tuples (as shown for example in `examples/tree_2d_dgsem/elixir_advection_diffusion_nonperiodic.jl`). If you’d like to apply the same condition to every face of the mesh, you can use the convenient `functions boundary_condition_default_tree_1D`, `boundary_condition_default_tree_2D` and `boundary_condition_default_tree_3D`. For example, in the two dimensional case: +For `TreeMesh`es, boundary conditions are defined and stored in named tuples (as shown for example in `examples/tree_2d_dgsem/elixir_advection_diffusion_nonperiodic.jl`). If you’d like to apply the same condition to every face of the mesh, you can use the convenient `functions boundary_condition_default_tree_1D`, `boundary_condition_default_tree_2D` and `boundary_condition_default_tree_3D`. For example, in the two dimensional case: ```julia initial_condition = initial_condition_eriksson_johnson From 170590df40379dc7b3d9314dfcf8599d9171f2f6 Mon Sep 17 00:00:00 2001 From: Vincent Marks Date: Thu, 25 Sep 2025 14:46:55 +0200 Subject: [PATCH 11/45] Update src/solvers/dgsem_unstructured/sort_boundary_conditions.jl Co-authored-by: Joshua Lampert <51029046+JoshuaLampert@users.noreply.github.com> --- src/solvers/dgsem_unstructured/sort_boundary_conditions.jl | 1 + 1 file changed, 1 insertion(+) diff --git a/src/solvers/dgsem_unstructured/sort_boundary_conditions.jl b/src/solvers/dgsem_unstructured/sort_boundary_conditions.jl index d01c71be624..f5c45c4c92e 100644 --- a/src/solvers/dgsem_unstructured/sort_boundary_conditions.jl +++ b/src/solvers/dgsem_unstructured/sort_boundary_conditions.jl @@ -154,6 +154,7 @@ end """ boundary_condition_default_3D + Create a default boundary condition dictionary for p4est meshes in 3D that use the standard boundary naming convention. This function applies the same boundary condition to all standard boundaries: From 4f0b06405b16eab625a50dd47675d3c0783c44df Mon Sep 17 00:00:00 2001 From: Vincent Marks Date: Thu, 25 Sep 2025 14:47:17 +0200 Subject: [PATCH 12/45] Update src/solvers/dgsem_unstructured/sort_boundary_conditions.jl Co-authored-by: Joshua Lampert <51029046+JoshuaLampert@users.noreply.github.com> --- src/solvers/dgsem_unstructured/sort_boundary_conditions.jl | 1 - 1 file changed, 1 deletion(-) diff --git a/src/solvers/dgsem_unstructured/sort_boundary_conditions.jl b/src/solvers/dgsem_unstructured/sort_boundary_conditions.jl index f5c45c4c92e..c01023fde82 100644 --- a/src/solvers/dgsem_unstructured/sort_boundary_conditions.jl +++ b/src/solvers/dgsem_unstructured/sort_boundary_conditions.jl @@ -168,7 +168,6 @@ This function applies the same boundary condition to all standard boundaries: - `boundary_condition`: The boundary condition function to apply to all boundaries # Returns - `Dict{Symbol, Any}`: Dictionary mapping boundary names to the boundary condition - """ function boundary_condition_default_p4est_3D(boundary_condition) From 446c35c7943f987527c855e6769eceb4c087ddc4 Mon Sep 17 00:00:00 2001 From: Vincent Marks Date: Thu, 25 Sep 2025 14:47:29 +0200 Subject: [PATCH 13/45] Update src/solvers/dgsem_unstructured/sort_boundary_conditions.jl Co-authored-by: Joshua Lampert <51029046+JoshuaLampert@users.noreply.github.com> --- src/solvers/dgsem_unstructured/sort_boundary_conditions.jl | 1 - 1 file changed, 1 deletion(-) diff --git a/src/solvers/dgsem_unstructured/sort_boundary_conditions.jl b/src/solvers/dgsem_unstructured/sort_boundary_conditions.jl index c01023fde82..b9d5431f7b7 100644 --- a/src/solvers/dgsem_unstructured/sort_boundary_conditions.jl +++ b/src/solvers/dgsem_unstructured/sort_boundary_conditions.jl @@ -193,7 +193,6 @@ This function applies the same boundary condition to all standard boundaries: # Returns - Named tuple mapping boundary names to the boundary condition """ - function boundary_condition_default_structured_1D(boundary_condition) return (x_neg = boundary_condition, From b2105989d421921e2e0cb2765140f123a928cb12 Mon Sep 17 00:00:00 2001 From: Vincent Marks Date: Thu, 25 Sep 2025 14:47:50 +0200 Subject: [PATCH 14/45] Update src/solvers/dgsem_unstructured/sort_boundary_conditions.jl Co-authored-by: Joshua Lampert <51029046+JoshuaLampert@users.noreply.github.com> --- src/solvers/dgsem_unstructured/sort_boundary_conditions.jl | 1 + 1 file changed, 1 insertion(+) diff --git a/src/solvers/dgsem_unstructured/sort_boundary_conditions.jl b/src/solvers/dgsem_unstructured/sort_boundary_conditions.jl index b9d5431f7b7..d0f9a12a7e3 100644 --- a/src/solvers/dgsem_unstructured/sort_boundary_conditions.jl +++ b/src/solvers/dgsem_unstructured/sort_boundary_conditions.jl @@ -219,6 +219,7 @@ function boundary_condition_default_structured_2D(boundary_condition) y_neg = boundary_condition, y_pos = boundary_condition) end + """ boundary_condition_default_structured_3D Create a default boundary condition dictionary for structured meshes in 3D From fdc03f477115449fd31faf3cb50b9052539a0cb6 Mon Sep 17 00:00:00 2001 From: Vincent Marks Date: Thu, 25 Sep 2025 14:48:19 +0200 Subject: [PATCH 15/45] Update src/solvers/dgsem_unstructured/sort_boundary_conditions.jl Co-authored-by: Joshua Lampert <51029046+JoshuaLampert@users.noreply.github.com> --- src/solvers/dgsem_unstructured/sort_boundary_conditions.jl | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/src/solvers/dgsem_unstructured/sort_boundary_conditions.jl b/src/solvers/dgsem_unstructured/sort_boundary_conditions.jl index d0f9a12a7e3..725cd755222 100644 --- a/src/solvers/dgsem_unstructured/sort_boundary_conditions.jl +++ b/src/solvers/dgsem_unstructured/sort_boundary_conditions.jl @@ -221,7 +221,8 @@ function boundary_condition_default_structured_2D(boundary_condition) end """ - boundary_condition_default_structured_3D + boundary_condition_default_structured_3D(boundary_condition) + Create a default boundary condition dictionary for structured meshes in 3D that use the standard boundary naming convention. This function applies the same boundary condition to all standard boundaries: From 1a13708aa88c339aa179cd9ea42ced3da7837cf4 Mon Sep 17 00:00:00 2001 From: Vincent Marks Date: Thu, 25 Sep 2025 14:48:45 +0200 Subject: [PATCH 16/45] Update src/solvers/dgsem_unstructured/sort_boundary_conditions.jl Co-authored-by: Joshua Lampert <51029046+JoshuaLampert@users.noreply.github.com> --- src/solvers/dgsem_unstructured/sort_boundary_conditions.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/solvers/dgsem_unstructured/sort_boundary_conditions.jl b/src/solvers/dgsem_unstructured/sort_boundary_conditions.jl index 725cd755222..047d46d896f 100644 --- a/src/solvers/dgsem_unstructured/sort_boundary_conditions.jl +++ b/src/solvers/dgsem_unstructured/sort_boundary_conditions.jl @@ -127,7 +127,7 @@ end """ - boundary_condition_default_p4est_2D + boundary_condition_default_p4est_2D(boundary_condition) Create a default boundary condition dictionary for p4est meshes in 2D that use the standard boundary naming convention. From b0c90fb7f2e5b96a833353781d5f6c15588a75e9 Mon Sep 17 00:00:00 2001 From: Vincent Marks Date: Thu, 25 Sep 2025 14:48:59 +0200 Subject: [PATCH 17/45] Update src/solvers/dgsem_unstructured/sort_boundary_conditions.jl Co-authored-by: Joshua Lampert <51029046+JoshuaLampert@users.noreply.github.com> --- src/solvers/dgsem_unstructured/sort_boundary_conditions.jl | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/src/solvers/dgsem_unstructured/sort_boundary_conditions.jl b/src/solvers/dgsem_unstructured/sort_boundary_conditions.jl index 047d46d896f..4cfb59734be 100644 --- a/src/solvers/dgsem_unstructured/sort_boundary_conditions.jl +++ b/src/solvers/dgsem_unstructured/sort_boundary_conditions.jl @@ -250,7 +250,8 @@ end """ - boundary_condition_default_tree_1D + boundary_condition_default_tree_1D(boundary_condition) + Create a default boundary condition dictionary for tree meshes in 1D that use the standard boundary naming convention. This function applies the same boundary condition to all standard boundaries: From a49e3254fb163d61e87fe026f84e043eb3314d81 Mon Sep 17 00:00:00 2001 From: Vincent Marks Date: Thu, 25 Sep 2025 14:49:12 +0200 Subject: [PATCH 18/45] Update src/solvers/dgsem_unstructured/sort_boundary_conditions.jl Co-authored-by: Joshua Lampert <51029046+JoshuaLampert@users.noreply.github.com> --- src/solvers/dgsem_unstructured/sort_boundary_conditions.jl | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/src/solvers/dgsem_unstructured/sort_boundary_conditions.jl b/src/solvers/dgsem_unstructured/sort_boundary_conditions.jl index 4cfb59734be..d43bc45eb79 100644 --- a/src/solvers/dgsem_unstructured/sort_boundary_conditions.jl +++ b/src/solvers/dgsem_unstructured/sort_boundary_conditions.jl @@ -270,7 +270,8 @@ function boundary_condition_default_tree_1D(boundary_condition) end """ - boundary_condition_default_tree_2D + boundary_condition_default_tree_2D(boundary_condition) + Create a default boundary condition dictionary for tree meshes in 2D that use the standard boundary naming convention. This function applies the same boundary condition to all standard boundaries: From 35dc6c9e5aa2d26a3ce0297a3760c57ad80fdef4 Mon Sep 17 00:00:00 2001 From: Vincent Marks Date: Thu, 25 Sep 2025 14:49:31 +0200 Subject: [PATCH 19/45] Update src/solvers/dgsem_unstructured/sort_boundary_conditions.jl Co-authored-by: Joshua Lampert <51029046+JoshuaLampert@users.noreply.github.com> --- src/solvers/dgsem_unstructured/sort_boundary_conditions.jl | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/src/solvers/dgsem_unstructured/sort_boundary_conditions.jl b/src/solvers/dgsem_unstructured/sort_boundary_conditions.jl index d43bc45eb79..b111ec641be 100644 --- a/src/solvers/dgsem_unstructured/sort_boundary_conditions.jl +++ b/src/solvers/dgsem_unstructured/sort_boundary_conditions.jl @@ -293,7 +293,8 @@ function boundary_condition_default_tree_2D(boundary_condition) end """ - boundary_condition_default_tree_3D + boundary_condition_default_tree_3D(boundary_condition) + Create a default boundary condition dictionary for tree meshes in 3D that use the standard boundary naming convention. This function applies the same boundary condition to all standard boundaries: From f668846d4faa9dd5ccc74f1da1b792f00cf84473 Mon Sep 17 00:00:00 2001 From: Vincent Marks Date: Thu, 25 Sep 2025 14:49:45 +0200 Subject: [PATCH 20/45] Update src/solvers/dgsem_unstructured/sort_boundary_conditions.jl Co-authored-by: Joshua Lampert <51029046+JoshuaLampert@users.noreply.github.com> --- src/solvers/dgsem_unstructured/sort_boundary_conditions.jl | 1 - 1 file changed, 1 deletion(-) diff --git a/src/solvers/dgsem_unstructured/sort_boundary_conditions.jl b/src/solvers/dgsem_unstructured/sort_boundary_conditions.jl index b111ec641be..d940dc5bc2b 100644 --- a/src/solvers/dgsem_unstructured/sort_boundary_conditions.jl +++ b/src/solvers/dgsem_unstructured/sort_boundary_conditions.jl @@ -309,7 +309,6 @@ This function applies the same boundary condition to all standard boundaries: # Returns - Named tuple mapping boundary names to the boundary condition """ - function boundary_condition_default_tree_3D(boundary_condition) return (x_neg = boundary_condition, From f08d444f6c4e76879bbc1f96d5a9c76b18318b85 Mon Sep 17 00:00:00 2001 From: Vincent Marks Date: Thu, 25 Sep 2025 14:50:04 +0200 Subject: [PATCH 21/45] Update src/solvers/dgsem_unstructured/sort_boundary_conditions.jl Co-authored-by: Joshua Lampert <51029046+JoshuaLampert@users.noreply.github.com> --- src/solvers/dgsem_unstructured/sort_boundary_conditions.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/solvers/dgsem_unstructured/sort_boundary_conditions.jl b/src/solvers/dgsem_unstructured/sort_boundary_conditions.jl index d940dc5bc2b..5306ee7977b 100644 --- a/src/solvers/dgsem_unstructured/sort_boundary_conditions.jl +++ b/src/solvers/dgsem_unstructured/sort_boundary_conditions.jl @@ -153,7 +153,7 @@ function boundary_condition_default_p4est_2D(boundary_condition) end """ - boundary_condition_default_3D + boundary_condition_default_3D(boundary_condition) Create a default boundary condition dictionary for p4est meshes in 3D that use the standard boundary naming convention. From 29a292397b7d8437ddeec8e55214ac4b3db26d6d Mon Sep 17 00:00:00 2001 From: Vincent Marks Date: Thu, 25 Sep 2025 14:50:55 +0200 Subject: [PATCH 22/45] Update src/solvers/dgsem_unstructured/sort_boundary_conditions.jl Co-authored-by: Joshua Lampert <51029046+JoshuaLampert@users.noreply.github.com> --- src/solvers/dgsem_unstructured/sort_boundary_conditions.jl | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/src/solvers/dgsem_unstructured/sort_boundary_conditions.jl b/src/solvers/dgsem_unstructured/sort_boundary_conditions.jl index 5306ee7977b..e652f5c497d 100644 --- a/src/solvers/dgsem_unstructured/sort_boundary_conditions.jl +++ b/src/solvers/dgsem_unstructured/sort_boundary_conditions.jl @@ -182,7 +182,8 @@ end """ - boundary_condition_default_structured_1D + boundary_condition_default_structured_1D(boundary_condition) + Create a default boundary condition dictionary for structured meshes in 1D that use the standard boundary naming convention. This function applies the same boundary condition to all standard boundaries: From ffbfe3916072a7cd0b4f5badeb8b30254065851a Mon Sep 17 00:00:00 2001 From: Vincent Marks Date: Thu, 25 Sep 2025 14:51:15 +0200 Subject: [PATCH 23/45] Update src/solvers/dgsem_unstructured/sort_boundary_conditions.jl Co-authored-by: Joshua Lampert <51029046+JoshuaLampert@users.noreply.github.com> --- src/solvers/dgsem_unstructured/sort_boundary_conditions.jl | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/src/solvers/dgsem_unstructured/sort_boundary_conditions.jl b/src/solvers/dgsem_unstructured/sort_boundary_conditions.jl index e652f5c497d..ca8cf267e32 100644 --- a/src/solvers/dgsem_unstructured/sort_boundary_conditions.jl +++ b/src/solvers/dgsem_unstructured/sort_boundary_conditions.jl @@ -200,7 +200,8 @@ function boundary_condition_default_structured_1D(boundary_condition) x_pos = boundary_condition) end """ - boundary_condition_default_structured_2D + boundary_condition_default_structured_2D(boundary_condition) + Create a default boundary condition dictionary for structured meshes in 2D that use the standard boundary naming convention. This function applies the same boundary condition to all standard boundaries: From 648596486a97f62688308058d6342eefedaa34cd Mon Sep 17 00:00:00 2001 From: Vincent Marks Date: Thu, 25 Sep 2025 18:26:18 +0200 Subject: [PATCH 24/45] Apply suggestions from code review Co-authored-by: Joshua Lampert <51029046+JoshuaLampert@users.noreply.github.com> --- docs/src/meshes/p4est_mesh.md | 2 +- docs/src/meshes/structured_mesh.md | 2 +- docs/src/meshes/tree_mesh.md | 2 +- src/solvers/dgsem_unstructured/sort_boundary_conditions.jl | 4 +--- 4 files changed, 4 insertions(+), 6 deletions(-) diff --git a/docs/src/meshes/p4est_mesh.md b/docs/src/meshes/p4est_mesh.md index a420ca517ce..4ebac036d32 100644 --- a/docs/src/meshes/p4est_mesh.md +++ b/docs/src/meshes/p4est_mesh.md @@ -941,7 +941,7 @@ reinitialize_boundaries!(semi.boundary_conditions, cache) # Needs to be called a This code could then be placed in the [`resize!`](https://github.com/trixi-framework/Trixi.jl/blob/eaeb04113523500ed831e3ab459694f12f7a49ea/src/time_integration/methods_2N.jl#L251-L255) function of a corresponding multirate integrator to ensure load-balancing for simulations involving AMR. ### Boundary conditions -For `P4est meshes`, boundary conditions are defined and stored in dictionaries (as shown for example in `examples/p4est_2d_dgsem/elixir_advection_diffusion_nonperiodic_amr.jl`). If you’d like to apply the same condition to every face of the mesh, you can use the convenient functions `boundary_condition_default_p4est_2D` and `boundary_condition_default_p4est_3D`. For example, in the two dimensional case: +For [`P4estMesh`](@ref)es, boundary conditions are defined and stored in dictionaries (as shown for example in `examples/p4est_2d_dgsem/elixir_advection_diffusion_nonperiodic_amr.jl`). If you’d like to apply the same condition to every face of the mesh, you can use the convenient functions `boundary_condition_default_p4est_2D` and `boundary_condition_default_p4est_3D`. For example, in the two dimensional case: ```julia boundary_condition = boundary_condition_slip_wall diff --git a/docs/src/meshes/structured_mesh.md b/docs/src/meshes/structured_mesh.md index 14f1ca04556..70f99e93785 100644 --- a/docs/src/meshes/structured_mesh.md +++ b/docs/src/meshes/structured_mesh.md @@ -13,7 +13,7 @@ the re-use of existing functionality for the [`TreeMesh`](@ref) but is usually less efficient, cf. [PR #550](https://github.com/trixi-framework/Trixi.jl/pull/550). ### Boundary conditions -For `Structured meshes`, boundary conditions are defined and stored in named tuples (as shown for example in `examples/structured_1d_dgsem/elixir_euler_source_terms_nonperiodic.jl`). If you’d like to apply the same condition to every face of the mesh, you can use the convenient functions `boundary_condition_default_structured_1D`, `boundary_condition_default_structured_2D` and `boundary_condition_default_structured_3D`. For example, in the one dimensional case: +For [`StructuredMesh`](@ref)es, boundary conditions are defined and stored in named tuples (as shown for example in `examples/structured_1d_dgsem/elixir_euler_source_terms_nonperiodic.jl`). If you’d like to apply the same condition to every face of the mesh, you can use the convenient functions `boundary_condition_default_structured_1D`, `boundary_condition_default_structured_2D` and `boundary_condition_default_structured_3D`. For example, in the one dimensional case: ```julia boundary_condition = BoundaryConditionDirichlet(initial_condition) diff --git a/docs/src/meshes/tree_mesh.md b/docs/src/meshes/tree_mesh.md index fe30e08332d..098d636e465 100644 --- a/docs/src/meshes/tree_mesh.md +++ b/docs/src/meshes/tree_mesh.md @@ -12,7 +12,7 @@ dispatching on the `orientation::Integer` as described in the ### Boundary conditions -For `TreeMesh`es, boundary conditions are defined and stored in named tuples (as shown for example in `examples/tree_2d_dgsem/elixir_advection_diffusion_nonperiodic.jl`). If you’d like to apply the same condition to every face of the mesh, you can use the convenient `functions boundary_condition_default_tree_1D`, `boundary_condition_default_tree_2D` and `boundary_condition_default_tree_3D`. For example, in the two dimensional case: +For [`TreeMesh`](@ref)es, boundary conditions are defined and stored in named tuples (as shown for example in `examples/tree_2d_dgsem/elixir_advection_diffusion_nonperiodic.jl`). If you’d like to apply the same condition to every face of the mesh, you can use the convenient `functions boundary_condition_default_tree_1D`, `boundary_condition_default_tree_2D` and `boundary_condition_default_tree_3D`. For example, in the two dimensional case: ```julia initial_condition = initial_condition_eriksson_johnson diff --git a/src/solvers/dgsem_unstructured/sort_boundary_conditions.jl b/src/solvers/dgsem_unstructured/sort_boundary_conditions.jl index ca8cf267e32..ccd20aa08b9 100644 --- a/src/solvers/dgsem_unstructured/sort_boundary_conditions.jl +++ b/src/solvers/dgsem_unstructured/sort_boundary_conditions.jl @@ -169,7 +169,6 @@ This function applies the same boundary condition to all standard boundaries: # Returns - `Dict{Symbol, Any}`: Dictionary mapping boundary names to the boundary condition """ - function boundary_condition_default_p4est_3D(boundary_condition) return Dict(:x_neg => boundary_condition, @@ -199,6 +198,7 @@ function boundary_condition_default_structured_1D(boundary_condition) return (x_neg = boundary_condition, x_pos = boundary_condition) end + """ boundary_condition_default_structured_2D(boundary_condition) @@ -250,7 +250,6 @@ function boundary_condition_default_structured_3D(boundary_condition) end - """ boundary_condition_default_tree_1D(boundary_condition) @@ -263,7 +262,6 @@ This function applies the same boundary condition to all standard boundaries: - `boundary_condition`: The boundary condition function to apply to all boundaries # Returns - Named tuple mapping boundary names to the boundary condition - """ function boundary_condition_default_tree_1D(boundary_condition) From a09495744bd879927cea10a71816f6da9b68bc4a Mon Sep 17 00:00:00 2001 From: vincmarks Date: Thu, 25 Sep 2025 19:23:46 +0200 Subject: [PATCH 25/45] added reference files to get an impression of the new function: boundary_condition_default(mesh, boundary_condition) --- docs/src/meshes/p4est_mesh.md | 8 ++------ docs/src/meshes/structured_mesh.md | 8 ++------ docs/src/meshes/tree_mesh.md | 9 ++------- 3 files changed, 6 insertions(+), 19 deletions(-) diff --git a/docs/src/meshes/p4est_mesh.md b/docs/src/meshes/p4est_mesh.md index a420ca517ce..2cc604060d7 100644 --- a/docs/src/meshes/p4est_mesh.md +++ b/docs/src/meshes/p4est_mesh.md @@ -941,9 +941,5 @@ reinitialize_boundaries!(semi.boundary_conditions, cache) # Needs to be called a This code could then be placed in the [`resize!`](https://github.com/trixi-framework/Trixi.jl/blob/eaeb04113523500ed831e3ab459694f12f7a49ea/src/time_integration/methods_2N.jl#L251-L255) function of a corresponding multirate integrator to ensure load-balancing for simulations involving AMR. ### Boundary conditions -For `P4est meshes`, boundary conditions are defined and stored in dictionaries (as shown for example in `examples/p4est_2d_dgsem/elixir_advection_diffusion_nonperiodic_amr.jl`). If you’d like to apply the same condition to every face of the mesh, you can use the convenient functions `boundary_condition_default_p4est_2D` and `boundary_condition_default_p4est_3D`. For example, in the two dimensional case: - -```julia -boundary_condition = boundary_condition_slip_wall -boundary_conditions = boundary_condition_default_p4est_2D(boundary_condition) -``` +For [`P4estMesh`](@ref)es, boundary conditions are defined and stored in dictionaries (see, for example, `examples/p4est_2d_dgsem/elixir_advection_diffusion_nonperiodic_amr.jl`). +If you want to apply the same boundary condition to all faces of the mesh, you can use the `boundary_condition_default(mesh, boundary_condition)` function, as demonstrated in `examples/p4est_2d_dgsem/elixir_advection_diffusion_nonperiodic_amr.jl` and `examples/p4est_3d_dgsem/elixir_euler_source_terms_nonperiodic.jl`. diff --git a/docs/src/meshes/structured_mesh.md b/docs/src/meshes/structured_mesh.md index 14f1ca04556..e915fefbf6d 100644 --- a/docs/src/meshes/structured_mesh.md +++ b/docs/src/meshes/structured_mesh.md @@ -13,9 +13,5 @@ the re-use of existing functionality for the [`TreeMesh`](@ref) but is usually less efficient, cf. [PR #550](https://github.com/trixi-framework/Trixi.jl/pull/550). ### Boundary conditions -For `Structured meshes`, boundary conditions are defined and stored in named tuples (as shown for example in `examples/structured_1d_dgsem/elixir_euler_source_terms_nonperiodic.jl`). If you’d like to apply the same condition to every face of the mesh, you can use the convenient functions `boundary_condition_default_structured_1D`, `boundary_condition_default_structured_2D` and `boundary_condition_default_structured_3D`. For example, in the one dimensional case: - -```julia -boundary_condition = BoundaryConditionDirichlet(initial_condition) -boundary_conditions = boundary_condition_default_structured_1D(boundary_condition) -``` +For [`StructuredMesh`](@ref)es, boundary conditions are defined and stored in dictionaries (see, for example, `examples/structured_1d_dgsem/elixir_euler_source_terms_nonperiodic.jl`). +If you want to apply the same boundary condition to all faces of the mesh, you can use the `boundary_condition_default(mesh, boundary_condition)` function, as demonstrated in `examples/structured_1d_dgsem/elixir_euler_source_terms_nonperiodic.jl`, `examples/structured_2d_dgsem/elixir_euler_rayleigh_taylor_instability.jl` and `examples/structured_3d_dgsem/elixir_euler_source_terms_nonperiodic_curved.jl`. diff --git a/docs/src/meshes/tree_mesh.md b/docs/src/meshes/tree_mesh.md index 26dfb4784f0..aef0692cc8e 100644 --- a/docs/src/meshes/tree_mesh.md +++ b/docs/src/meshes/tree_mesh.md @@ -12,10 +12,5 @@ dispatching on the `orientation::Integer` as described in the ### Boundary conditions -For `Tree meshes`, boundary conditions are defined and stored in named tuples (as shown for example in `examples/tree_2d_dgsem/elixir_advection_diffusion_nonperiodic.jl`). If you’d like to apply the same condition to every face of the mesh, you can use the convenient `functions boundary_condition_default_tree_1D`, `boundary_condition_default_tree_2D` and `boundary_condition_default_tree_3D`. For example, in the two dimensional case: - -```julia -initial_condition = initial_condition_eriksson_johnson - -boundary_conditions = boundary_condition_default_tree_2D(BoundaryConditionDirichlet(initial_condition)) -``` \ No newline at end of file +For [`TreeMesh`](@ref)es, boundary conditions are defined and stored in dictionaries (see, for example, `examples/tree_1d_dgsem/elixir_euler_source_terms_nonperiodic.jl`). +If you want to apply the same boundary condition to all faces of the mesh, you can use the `boundary_condition_default(mesh, boundary_condition)` function, as demonstrated in `examples/tree_1d_dgsem/elixir_euler_source_terms_nonperiodic.jl`, `examples/tree_2d_dgsem/elixir_euler_source_terms_nonperiodic.jl` and `examples/tree_3d_dgsem/elixir_hypdiff_nonperiodic.jl`. \ No newline at end of file From 231ac02c1f7e139f5002a8e41aa337f95cb410a9 Mon Sep 17 00:00:00 2001 From: vincmarks Date: Thu, 25 Sep 2025 19:28:17 +0200 Subject: [PATCH 26/45] adding examples for the boundary_condition_default(mesh, boundary_condition) function --- .../elixir_advection_diffusion_nonperiodic_amr.jl | 4 ++++ .../p4est_3d_dgsem/elixir_euler_source_terms_nonperiodic.jl | 3 +++ .../elixir_euler_source_terms_nonperiodic.jl | 3 +++ .../elixir_euler_rayleigh_taylor_instability.jl | 4 ++++ .../elixir_euler_source_terms_nonperiodic_curved.jl | 3 +++ .../tree_1d_dgsem/elixir_euler_source_terms_nonperiodic.jl | 3 +++ .../tree_2d_dgsem/elixir_euler_source_terms_nonperiodic.jl | 3 +++ examples/tree_3d_dgsem/elixir_hypdiff_nonperiodic.jl | 4 ++++ 8 files changed, 27 insertions(+) diff --git a/examples/p4est_2d_dgsem/elixir_advection_diffusion_nonperiodic_amr.jl b/examples/p4est_2d_dgsem/elixir_advection_diffusion_nonperiodic_amr.jl index 2ba5f08d4fc..f9de74d9af1 100644 --- a/examples/p4est_2d_dgsem/elixir_advection_diffusion_nonperiodic_amr.jl +++ b/examples/p4est_2d_dgsem/elixir_advection_diffusion_nonperiodic_amr.jl @@ -50,6 +50,10 @@ boundary_conditions_parabolic = Dict(:x_neg => BoundaryConditionDirichlet(initia :y_neg => BoundaryConditionDirichlet(initial_condition), :y_pos => BoundaryConditionDirichlet(initial_condition)) +# Alternative to the boundary_conditions_parabolic defined above: +# boundary_condition = BoundaryConditionDirichlet(initial_condition) +# boundary_conditions_parabolic = boundary_condition_default(mesh, boundary_condition) + # A semidiscretization collects data structures and functions for the spatial discretization semi = SemidiscretizationHyperbolicParabolic(mesh, (equations, equations_parabolic), diff --git a/examples/p4est_3d_dgsem/elixir_euler_source_terms_nonperiodic.jl b/examples/p4est_3d_dgsem/elixir_euler_source_terms_nonperiodic.jl index 6d57205a169..c81d179f1bb 100644 --- a/examples/p4est_3d_dgsem/elixir_euler_source_terms_nonperiodic.jl +++ b/examples/p4est_3d_dgsem/elixir_euler_source_terms_nonperiodic.jl @@ -35,6 +35,9 @@ mesh = P4estMesh(trees_per_dimension, polydeg = 1, coordinates_min = coordinates_min, coordinates_max = coordinates_max, periodicity = false, initial_refinement_level = 1) +# Alternative to the boundary_conditions defined above: +# boundary_conditions = boundary_condition_default(mesh, boundary_condition) + semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver, source_terms = source_terms_convergence_test, boundary_conditions = boundary_conditions) diff --git a/examples/structured_1d_dgsem/elixir_euler_source_terms_nonperiodic.jl b/examples/structured_1d_dgsem/elixir_euler_source_terms_nonperiodic.jl index 9284b79539e..1c83774dbfb 100644 --- a/examples/structured_1d_dgsem/elixir_euler_source_terms_nonperiodic.jl +++ b/examples/structured_1d_dgsem/elixir_euler_source_terms_nonperiodic.jl @@ -30,6 +30,9 @@ f1() = SVector(0.0) f2() = SVector(2.0) mesh = StructuredMesh((16,), (f1, f2), periodicity = false) +# Alternative to the boundary_conditions defined above: +# boundary_conditions = boundary_condition_default(mesh, boundary_condition) + semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver, source_terms = source_terms, boundary_conditions = boundary_conditions) diff --git a/examples/structured_2d_dgsem/elixir_euler_rayleigh_taylor_instability.jl b/examples/structured_2d_dgsem/elixir_euler_rayleigh_taylor_instability.jl index adf8b231019..805735b6ad4 100644 --- a/examples/structured_2d_dgsem/elixir_euler_rayleigh_taylor_instability.jl +++ b/examples/structured_2d_dgsem/elixir_euler_rayleigh_taylor_instability.jl @@ -76,6 +76,10 @@ boundary_conditions = (x_neg = boundary_condition_slip_wall, y_neg = boundary_condition_slip_wall, y_pos = boundary_condition_slip_wall) +# Alternative to the boundary_conditions defined above: +# boundary_condition = boundary_condition_slip_wall +# boundary_conditions = boundary_condition_default(mesh, boundary_condition) + # # Alternative setup: left/right periodic BCs and Dirichlet BCs on the top/bottom. # boundary_conditions = ( # x_neg=boundary_condition_periodic, diff --git a/examples/structured_3d_dgsem/elixir_euler_source_terms_nonperiodic_curved.jl b/examples/structured_3d_dgsem/elixir_euler_source_terms_nonperiodic_curved.jl index b0c5234e5b1..1ffc2bc2b80 100644 --- a/examples/structured_3d_dgsem/elixir_euler_source_terms_nonperiodic_curved.jl +++ b/examples/structured_3d_dgsem/elixir_euler_source_terms_nonperiodic_curved.jl @@ -57,6 +57,9 @@ end cells_per_dimension = (4, 4, 4) mesh = StructuredMesh(cells_per_dimension, mapping, periodicity = false) +# Alternative to the boundary_conditions defined above: +# boundary_conditions = boundary_condition_default(mesh, boundary_condition) + semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver, source_terms = source_terms_convergence_test, boundary_conditions = boundary_conditions) diff --git a/examples/tree_1d_dgsem/elixir_euler_source_terms_nonperiodic.jl b/examples/tree_1d_dgsem/elixir_euler_source_terms_nonperiodic.jl index e51ca955c91..715768d01a9 100644 --- a/examples/tree_1d_dgsem/elixir_euler_source_terms_nonperiodic.jl +++ b/examples/tree_1d_dgsem/elixir_euler_source_terms_nonperiodic.jl @@ -31,6 +31,9 @@ mesh = TreeMesh(coordinates_min, coordinates_max, n_cells_max = 10_000, periodicity = false) +# Alternative to the boundary_conditions defined above: +# boundary_conditions = boundary_condition_default(mesh, boundary_condition) + semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver, source_terms = source_terms_convergence_test, boundary_conditions = boundary_conditions) diff --git a/examples/tree_2d_dgsem/elixir_euler_source_terms_nonperiodic.jl b/examples/tree_2d_dgsem/elixir_euler_source_terms_nonperiodic.jl index 9f0bcbda2c2..e9780a3700a 100644 --- a/examples/tree_2d_dgsem/elixir_euler_source_terms_nonperiodic.jl +++ b/examples/tree_2d_dgsem/elixir_euler_source_terms_nonperiodic.jl @@ -32,6 +32,9 @@ mesh = TreeMesh(coordinates_min, coordinates_max, n_cells_max = 10_000, periodicity = false) +# Alternative to the boundary_conditions defined above: +# boundary_conditions = boundary_condition_default(mesh, boundary_condition) + semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver, source_terms = source_terms_convergence_test, boundary_conditions = boundary_conditions) diff --git a/examples/tree_3d_dgsem/elixir_hypdiff_nonperiodic.jl b/examples/tree_3d_dgsem/elixir_hypdiff_nonperiodic.jl index 8c1b9f76341..0dc23d736b9 100644 --- a/examples/tree_3d_dgsem/elixir_hypdiff_nonperiodic.jl +++ b/examples/tree_3d_dgsem/elixir_hypdiff_nonperiodic.jl @@ -22,6 +22,10 @@ mesh = TreeMesh(coordinates_min, coordinates_max, n_cells_max = 30_000, periodicity = (false, true, true)) +# Alternative to the boundary_conditions defined above if you want to apply the same boundary condition to all faces of the mesh: +# boundary_condition = boundary_condition_periodic +# boundary_conditions = boundary_condition_default(mesh, boundary_condition) + semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver, source_terms = source_terms_poisson_nonperiodic, boundary_conditions = boundary_conditions) From 75be46e6f229342dad30d710a9d2166c4a3c6f88 Mon Sep 17 00:00:00 2001 From: vincmarks Date: Thu, 25 Sep 2025 19:29:28 +0200 Subject: [PATCH 27/45] included the default_boundary_conditions.jl file and the boundary_condition_default function --- src/Trixi.jl | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/src/Trixi.jl b/src/Trixi.jl index c62546e422d..9e9e08ac6c5 100644 --- a/src/Trixi.jl +++ b/src/Trixi.jl @@ -141,6 +141,7 @@ include("auxiliary/t8code.jl") include("equations/equations.jl") include("meshes/meshes.jl") include("solvers/solvers.jl") +include("solvers/default_boundary_conditions.jl") include("equations/equations_parabolic.jl") # these depend on parabolic solver types include("semidiscretization/semidiscretization.jl") include("semidiscretization/semidiscretization_hyperbolic.jl") @@ -151,7 +152,6 @@ include("time_integration/time_integration.jl") include("callbacks_step/callbacks_step.jl") include("callbacks_stage/callbacks_stage.jl") include("semidiscretization/semidiscretization_euler_gravity.jl") - # Special elixirs such as `convergence_test` include("auxiliary/special_elixirs.jl") @@ -228,6 +228,7 @@ export boundary_condition_do_nothing, boundary_condition_default_tree_1D, boundary_condition_default_tree_2D, boundary_condition_default_tree_3D, + boundary_condition_default, boundary_condition_periodic, BoundaryConditionDirichlet, BoundaryConditionNeumann, From 0c00a4aee182c25f8e2ed163355c70f08e33cafe Mon Sep 17 00:00:00 2001 From: vincmarks Date: Thu, 25 Sep 2025 19:32:51 +0200 Subject: [PATCH 28/45] create a new file (default_boundary_conditions.jl)for the boundary_condition_default(mesh, boundary_condition) function instead of writing it in the sort_boundary_conditions.jl file --- src/solvers/default_boundary_conditions.jl | 241 ++++++++++++++++++ .../sort_boundary_conditions.jl | 194 +------------- 2 files changed, 242 insertions(+), 193 deletions(-) create mode 100644 src/solvers/default_boundary_conditions.jl diff --git a/src/solvers/default_boundary_conditions.jl b/src/solvers/default_boundary_conditions.jl new file mode 100644 index 00000000000..278addedcf1 --- /dev/null +++ b/src/solvers/default_boundary_conditions.jl @@ -0,0 +1,241 @@ +""" + boundary_condition_default_p4est_2D(boundary_condition) + +Create a default boundary condition dictionary for p4est meshes in 2D +that use the standard boundary naming convention. + +This function applies the same boundary condition to all standard boundaries: +- `:x_neg`: negative x-direction boundary +- `:x_pos`: positive x-direction boundary +- `:y_neg`: negative y-direction boundary +- `:y_pos`: positive y-direction boundary + +# Arguments +- `boundary_condition`: The boundary condition function to apply to all boundaries + +# Returns +- `Dict{Symbol, Any}`: Dictionary mapping boundary names to the boundary condition +""" +function boundary_condition_default_p4est_2D(boundary_condition) + + return Dict(:x_neg => boundary_condition, + :y_neg => boundary_condition, + :y_pos => boundary_condition, + :x_pos => boundary_condition) +end + +""" + boundary_condition_default_p4est_3D(boundary_condition) + +Create a default boundary condition dictionary for p4est meshes in 3D +that use the standard boundary naming convention. +This function applies the same boundary condition to all standard boundaries: +- `:x_neg`: negative x-direction boundary +- `:x_pos`: positive x-direction boundary +- `:y_neg`: negative y-direction boundary +- `:y_pos`: positive y-direction boundary +- `:z_neg`: negative z-direction boundary +- `:z_pos`: positive z-direction boundary +# Arguments +- `boundary_condition`: The boundary condition function to apply to all boundaries +# Returns +- `Dict{Symbol, Any}`: Dictionary mapping boundary names to the boundary condition +""" +function boundary_condition_default_p4est_3D(boundary_condition) + + return Dict(:x_neg => boundary_condition, + :x_pos => boundary_condition, + :y_neg => boundary_condition, + :y_pos => boundary_condition, + :z_neg => boundary_condition, + :z_pos => boundary_condition) +end + + +""" + boundary_condition_default_structured_1D(boundary_condition) + +Create a default boundary condition dictionary for structured meshes in 1D +that use the standard boundary naming convention. +This function applies the same boundary condition to all standard boundaries: +- `:x_neg`: negative x-direction boundary +- `:x_pos`: positive x-direction boundary +# Arguments +- `boundary_condition`: The boundary condition function to apply to all boundaries +# Returns +- Named tuple mapping boundary names to the boundary condition +""" +function boundary_condition_default_structured_1D(boundary_condition) + + return (x_neg = boundary_condition, + x_pos = boundary_condition) +end + +""" + boundary_condition_default_structured_2D(boundary_condition) + +Create a default boundary condition dictionary for structured meshes in 2D +that use the standard boundary naming convention. +This function applies the same boundary condition to all standard boundaries: +- `:x_neg`: negative x-direction boundary +- `:x_pos`: positive x-direction boundary +- `:y_neg`: negative y-direction boundary +- `:y_pos`: positive y-direction boundary +# Arguments +- `boundary_condition`: The boundary condition function to apply to all boundaries +# Returns +- Named tuple mapping boundary names to the boundary condition +""" +function boundary_condition_default_structured_2D(boundary_condition) + + return (x_neg = boundary_condition, + x_pos = boundary_condition, + y_neg = boundary_condition, + y_pos = boundary_condition) +end + +""" + boundary_condition_default_structured_3D(boundary_condition) + +Create a default boundary condition dictionary for structured meshes in 3D +that use the standard boundary naming convention. +This function applies the same boundary condition to all standard boundaries: +- `:x_neg`: negative x-direction boundary +- `:x_pos`: positive x-direction boundary +- `:y_neg`: negative y-direction boundary +- `:y_pos`: positive y-direction boundary +- `:z_neg`: negative z-direction boundary +- `:z_pos`: positive z-direction boundary +# Arguments +- `boundary_condition`: The boundary condition function to apply to all boundaries +# Returns +- Named tuple mapping boundary names to the boundary condition +""" +function boundary_condition_default_structured_3D(boundary_condition) + + return (x_neg = boundary_condition, + x_pos = boundary_condition, + y_neg = boundary_condition, + y_pos = boundary_condition, + z_neg = boundary_condition, + z_pos = boundary_condition) + +end + +""" + boundary_condition_default_tree_1D(boundary_condition) + +Create a default boundary condition dictionary for tree meshes in 1D +that use the standard boundary naming convention. +This function applies the same boundary condition to all standard boundaries: +- `:x_neg`: negative x-direction boundary +- `:x_pos`: positive x-direction boundary +# Arguments +- `boundary_condition`: The boundary condition function to apply to all boundaries +# Returns +- Named tuple mapping boundary names to the boundary condition +""" +function boundary_condition_default_tree_1D(boundary_condition) + + return (x_neg = boundary_condition, + x_pos = boundary_condition) +end + +""" + boundary_condition_default_tree_2D(boundary_condition) + +Create a default boundary condition dictionary for tree meshes in 2D +that use the standard boundary naming convention. +This function applies the same boundary condition to all standard boundaries: +- `:x_neg`: negative x-direction boundary +- `:x_pos`: positive x-direction boundary +- `:y_neg`: negative y-direction boundary +- `:y_pos`: positive y-direction boundary +# Arguments +- `boundary_condition`: The boundary condition function to apply to all boundaries +# Returns +- Named tuple mapping boundary names to the boundary condition +""" +function boundary_condition_default_tree_2D(boundary_condition) + + return (x_neg = boundary_condition, + x_pos = boundary_condition, + y_neg = boundary_condition, + y_pos = boundary_condition) +end + +""" + boundary_condition_default_tree_3D(boundary_condition) + +Create a default boundary condition dictionary for tree meshes in 3D +that use the standard boundary naming convention. +This function applies the same boundary condition to all standard boundaries: +- `:x_neg`: negative x-direction boundary +- `:x_pos`: positive x-direction boundary +- `:y_neg`: negative y-direction boundary +- `:y_pos`: positive y-direction boundary +- `:z_neg`: negative z-direction boundary +- `:z_pos`: positive z-direction boundary +# Arguments +- `boundary_condition`: The boundary condition function to apply to all boundaries +# Returns +- Named tuple mapping boundary names to the boundary condition +""" +function boundary_condition_default_tree_3D(boundary_condition) + + return (x_neg = boundary_condition, + x_pos = boundary_condition, + y_neg = boundary_condition, + y_pos = boundary_condition, + z_neg = boundary_condition, + z_pos = boundary_condition) +end + +""" + boundary_condition_default(mesh, boundary_condition) + +Create a default boundary condition dictionary based on the mesh type. +This function automatically determines the appropriate boundary condition +format based on the mesh type and applies the given boundary condition +to all boundaries. + +# Arguments +- `mesh`: The mesh object (TreeMesh, StructuredMesh, P4estMesh) +- `boundary_condition`: The boundary condition function to apply to all boundaries + +# Returns +- Dictionary or named tuple mapping boundary names to the boundary condition +""" +function boundary_condition_default(mesh, boundary_condition) + + if isa(mesh, TreeMesh) + ndims_mesh = ndims(mesh) + if ndims_mesh == 1 + return boundary_condition_default_tree_1D(boundary_condition) + elseif ndims_mesh == 2 + return boundary_condition_default_tree_2D(boundary_condition) + elseif ndims_mesh == 3 + return boundary_condition_default_tree_3D(boundary_condition) + end + elseif isa(mesh, StructuredMesh) + ndims_mesh = ndims(mesh) + if ndims_mesh == 1 + return boundary_condition_default_structured_1D(boundary_condition) + elseif ndims_mesh == 2 + return boundary_condition_default_structured_2D(boundary_condition) + elseif ndims_mesh == 3 + return boundary_condition_default_structured_3D(boundary_condition) + end + elseif isa(mesh, P4estMesh) + ndims_mesh = ndims(mesh) + if ndims_mesh == 2 + return boundary_condition_default_p4est_2D(boundary_condition) + elseif ndims_mesh == 3 + return boundary_condition_default_p4est_3D(boundary_condition) + end + + end + + error("No default boundary conditions available for mesh type $(typeof(mesh)). " * + "Please specify the boundary conditions manually.") +end \ No newline at end of file diff --git a/src/solvers/dgsem_unstructured/sort_boundary_conditions.jl b/src/solvers/dgsem_unstructured/sort_boundary_conditions.jl index 7ac5be33101..992501c08d6 100644 --- a/src/solvers/dgsem_unstructured/sort_boundary_conditions.jl +++ b/src/solvers/dgsem_unstructured/sort_boundary_conditions.jl @@ -124,196 +124,4 @@ end # @eval due to @muladd @eval Adapt.@adapt_structure(UnstructuredSortedBoundaryTypes) - - -""" - boundary_condition_default_p4est_2D - -Create a default boundary condition dictionary for p4est meshes in 2D -that use the standard boundary naming convention. - -This function applies the same boundary condition to all standard boundaries: -- `:x_neg`: negative x-direction boundary -- `:x_pos`: positive x-direction boundary -- `:y_neg`: negative y-direction boundary -- `:y_pos`: positive y-direction boundary - -# Arguments -- `boundary_condition`: The boundary condition function to apply to all boundaries - -# Returns -- `Dict{Symbol, Any}`: Dictionary mapping boundary names to the boundary condition - -""" -function boundary_condition_default_p4est_2D(boundary_condition) - - return Dict(:x_neg => boundary_condition, - :y_neg => boundary_condition, - :y_pos => boundary_condition, - :x_pos => boundary_condition) -end - -""" - boundary_condition_default_3D -Create a default boundary condition dictionary for p4est meshes in 3D -that use the standard boundary naming convention. -This function applies the same boundary condition to all standard boundaries: -- `:x_neg`: negative x-direction boundary -- `:x_pos`: positive x-direction boundary -- `:y_neg`: negative y-direction boundary -- `:y_pos`: positive y-direction boundary -- `:z_neg`: negative z-direction boundary -- `:z_pos`: positive z-direction boundary -# Arguments -- `boundary_condition`: The boundary condition function to apply to all boundaries -# Returns -- `Dict{Symbol, Any}`: Dictionary mapping boundary names to the boundary condition - -""" - -function boundary_condition_default_p4est_3D(boundary_condition) - - return Dict(:x_neg => boundary_condition, - :x_pos => boundary_condition, - :y_neg => boundary_condition, - :y_pos => boundary_condition, - :z_neg => boundary_condition, - :z_pos => boundary_condition) -end - - -""" - boundary_condition_default_structured_1D -Create a default boundary condition dictionary for structured meshes in 1D -that use the standard boundary naming convention. -This function applies the same boundary condition to all standard boundaries: -- `:x_neg`: negative x-direction boundary -- `:x_pos`: positive x-direction boundary -# Arguments -- `boundarycondition`: The boundary condition function to apply to all boundaries -# Returns -- Named tuple mapping boundary names to the boundary condition -""" - -function boundary_condition_default_structured_1D(boundary_condition) - - return (x_neg = boundary_condition, - x_pos = boundary_condition) -end -""" - boundary_condition_default_structured_2D -Create a default boundary condition dictionary for structured meshes in 2D -that use the standard boundary naming convention. -This function applies the same boundary condition to all standard boundaries: -- `:x_neg`: negative x-direction boundary -- `:x_pos`: positive x-direction boundary -- `:y_neg`: negative y-direction boundary -- `:y_pos`: positive y-direction boundary -# Arguments -- `boundarycondition`: The boundary condition function to apply to all boundaries -# Returns -- Named tuple mapping boundary names to the boundary condition -""" -function boundary_condition_default_structured_2D(boundary_condition) - - return (x_neg = boundary_condition, - x_pos = boundary_condition, - y_neg = boundary_condition, - y_pos = boundary_condition) -end -""" - boundary_condition_default_structured_3D -Create a default boundary condition dictionary for structured meshes in 3D -that use the standard boundary naming convention. -This function applies the same boundary condition to all standard boundaries: -- `:x_neg`: negative x-direction boundary -- `:x_pos`: positive x-direction boundary -- `:y_neg`: negative y-direction boundary -- `:y_pos`: positive y-direction boundary -- `:z_neg`: negative z-direction boundary -- `:z_pos`: positive z-direction boundary -# Arguments -- `boundarycondition`: The boundary condition function to apply to all boundaries -# Returns -- Named tuple mapping boundary names to the boundary condition -""" -function boundary_condition_default_structured_3D(boundary_condition) - - return (x_neg = boundary_condition, - x_pos = boundary_condition, - y_neg = boundary_condition, - y_pos = boundary_condition, - z_neg = boundary_condition, - z_pos = boundary_condition) - -end - - -""" - boundary_condition_default_tree_1D -Create a default boundary condition dictionary for tree meshes in 1D -that use the standard boundary naming convention. -This function applies the same boundary condition to all standard boundaries: -- `:x_neg`: negative x-direction boundary -- `:x_pos`: positive x-direction boundary -# Arguments -- `boundary_condition`: The boundary condition function to apply to all boundaries -# Returns -- Named tuple mapping boundary names to the boundary condition - -""" -function boundary_condition_default_tree_1D(boundary_condition) - - return (x_neg = boundary_condition, - x_pos = boundary_condition) -end - -""" - boundary_condition_default_tree_2D -Create a default boundary condition dictionary for tree meshes in 2D -that use the standard boundary naming convention. -This function applies the same boundary condition to all standard boundaries: -- `:x_neg`: negative x-direction boundary -- `:x_pos`: positive x-direction boundary -- `:y_neg`: negative y-direction boundary -- `:y_pos`: positive y-direction boundary -# Arguments -- `boundary_condition`: The boundary condition function to apply to all boundaries -# Returns -- Named tuple mapping boundary names to the boundary condition -""" -function boundary_condition_default_tree_2D(boundary_condition) - - return (x_neg = boundary_condition, - x_pos = boundary_condition, - y_neg = boundary_condition, - y_pos = boundary_condition) -end - -""" - boundary_condition_default_tree_3D -Create a default boundary condition dictionary for tree meshes in 3D -that use the standard boundary naming convention. -This function applies the same boundary condition to all standard boundaries: -- `:x_neg`: negative x-direction boundary -- `:x_pos`: positive x-direction boundary -- `:y_neg`: negative y-direction boundary -- `:y_pos`: positive y-direction boundary -- `:z_neg`: negative z-direction boundary -- `:z_pos`: positive z-direction boundary -# Arguments -- `boundary_condition`: The boundary condition function to apply to all boundaries -# Returns -- Named tuple mapping boundary names to the boundary condition -""" - -function boundary_condition_default_tree_3D(boundary_condition) - - return (x_neg = boundary_condition, - x_pos = boundary_condition, - y_neg = boundary_condition, - y_pos = boundary_condition, - z_neg = boundary_condition, - z_pos = boundary_condition) -end -end#@muladd +end#@muladd \ No newline at end of file From dbc823b1961195098201fa6088a0ca61a769aa2b Mon Sep 17 00:00:00 2001 From: vincmarks Date: Thu, 25 Sep 2025 19:56:56 +0200 Subject: [PATCH 29/45] fixed typo: named tuples instead of dicitonaries --- docs/src/meshes/structured_mesh.md | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/docs/src/meshes/structured_mesh.md b/docs/src/meshes/structured_mesh.md index e915fefbf6d..85a50f37ad3 100644 --- a/docs/src/meshes/structured_mesh.md +++ b/docs/src/meshes/structured_mesh.md @@ -13,5 +13,5 @@ the re-use of existing functionality for the [`TreeMesh`](@ref) but is usually less efficient, cf. [PR #550](https://github.com/trixi-framework/Trixi.jl/pull/550). ### Boundary conditions -For [`StructuredMesh`](@ref)es, boundary conditions are defined and stored in dictionaries (see, for example, `examples/structured_1d_dgsem/elixir_euler_source_terms_nonperiodic.jl`). +For [`StructuredMesh`](@ref)es, boundary conditions are defined and stored in named tuples (see, for example, `examples/structured_1d_dgsem/elixir_euler_source_terms_nonperiodic.jl`). If you want to apply the same boundary condition to all faces of the mesh, you can use the `boundary_condition_default(mesh, boundary_condition)` function, as demonstrated in `examples/structured_1d_dgsem/elixir_euler_source_terms_nonperiodic.jl`, `examples/structured_2d_dgsem/elixir_euler_rayleigh_taylor_instability.jl` and `examples/structured_3d_dgsem/elixir_euler_source_terms_nonperiodic_curved.jl`. From e86c2baacc423da369914f6c86545dc1d4661901 Mon Sep 17 00:00:00 2001 From: vincmarks Date: Sun, 19 Oct 2025 19:43:54 +0200 Subject: [PATCH 30/45] delete unnecessary functions --- src/Trixi.jl | 8 -------- 1 file changed, 8 deletions(-) diff --git a/src/Trixi.jl b/src/Trixi.jl index 3a027022a8c..7a67833ab87 100644 --- a/src/Trixi.jl +++ b/src/Trixi.jl @@ -220,14 +220,6 @@ export initial_condition_constant, initial_condition_weak_blast_wave export boundary_condition_do_nothing, - boundary_condition_default_p4est_2D, - boundary_condition_default_p4est_3D, - boundary_condition_default_structured_1D, - boundary_condition_default_structured_2D, - boundary_condition_default_structured_3D, - boundary_condition_default_tree_1D, - boundary_condition_default_tree_2D, - boundary_condition_default_tree_3D, boundary_condition_default, boundary_condition_periodic, BoundaryConditionDirichlet, From 2ead177e7ce8ad14583ca8ecc5663bb355ff58ce Mon Sep 17 00:00:00 2001 From: vincmarks Date: Sun, 19 Oct 2025 19:44:41 +0200 Subject: [PATCH 31/45] only use boundary_condition_default function --- src/solvers/default_boundary_conditions.jl | 83 +++++----------------- 1 file changed, 16 insertions(+), 67 deletions(-) diff --git a/src/solvers/default_boundary_conditions.jl b/src/solvers/default_boundary_conditions.jl index 278addedcf1..a4796b15208 100644 --- a/src/solvers/default_boundary_conditions.jl +++ b/src/solvers/default_boundary_conditions.jl @@ -1,5 +1,5 @@ """ - boundary_condition_default_p4est_2D(boundary_condition) + boundary_condition_default(mesh::P4estMesh{2,2}, boundary_condition) Create a default boundary condition dictionary for p4est meshes in 2D that use the standard boundary naming convention. @@ -16,7 +16,7 @@ This function applies the same boundary condition to all standard boundaries: # Returns - `Dict{Symbol, Any}`: Dictionary mapping boundary names to the boundary condition """ -function boundary_condition_default_p4est_2D(boundary_condition) +function boundary_condition_default(mesh::P4estMesh{2,2}, boundary_condition) return Dict(:x_neg => boundary_condition, :y_neg => boundary_condition, @@ -25,7 +25,7 @@ function boundary_condition_default_p4est_2D(boundary_condition) end """ - boundary_condition_default_p4est_3D(boundary_condition) + boundary_condition_default(mesh::P4estMesh{3,3}, boundary_condition) Create a default boundary condition dictionary for p4est meshes in 3D that use the standard boundary naming convention. @@ -41,7 +41,7 @@ This function applies the same boundary condition to all standard boundaries: # Returns - `Dict{Symbol, Any}`: Dictionary mapping boundary names to the boundary condition """ -function boundary_condition_default_p4est_3D(boundary_condition) +function boundary_condition_default(mesh::P4estMesh{3,3}, boundary_condition) return Dict(:x_neg => boundary_condition, :x_pos => boundary_condition, @@ -51,9 +51,8 @@ function boundary_condition_default_p4est_3D(boundary_condition) :z_pos => boundary_condition) end - """ - boundary_condition_default_structured_1D(boundary_condition) + boundary_condition_default(mesh::StructuredMesh1D, boundary_condition) Create a default boundary condition dictionary for structured meshes in 1D that use the standard boundary naming convention. @@ -65,14 +64,14 @@ This function applies the same boundary condition to all standard boundaries: # Returns - Named tuple mapping boundary names to the boundary condition """ -function boundary_condition_default_structured_1D(boundary_condition) +function boundary_condition_default(mesh::StructuredMesh{1}, boundary_condition) return (x_neg = boundary_condition, x_pos = boundary_condition) end """ - boundary_condition_default_structured_2D(boundary_condition) + boundary_condition_default(mesh::StructuredMesh2D, boundary_condition) Create a default boundary condition dictionary for structured meshes in 2D that use the standard boundary naming convention. @@ -86,7 +85,7 @@ This function applies the same boundary condition to all standard boundaries: # Returns - Named tuple mapping boundary names to the boundary condition """ -function boundary_condition_default_structured_2D(boundary_condition) +function boundary_condition_default(mesh::StructuredMesh{2}, boundary_condition) return (x_neg = boundary_condition, x_pos = boundary_condition, @@ -95,7 +94,7 @@ function boundary_condition_default_structured_2D(boundary_condition) end """ - boundary_condition_default_structured_3D(boundary_condition) + boundary_condition_default(mesh::StructuredMesh3D, boundary_condition) Create a default boundary condition dictionary for structured meshes in 3D that use the standard boundary naming convention. @@ -111,7 +110,7 @@ This function applies the same boundary condition to all standard boundaries: # Returns - Named tuple mapping boundary names to the boundary condition """ -function boundary_condition_default_structured_3D(boundary_condition) +function boundary_condition_default(mesh::StructuredMesh{3}, boundary_condition) return (x_neg = boundary_condition, x_pos = boundary_condition, @@ -123,7 +122,7 @@ function boundary_condition_default_structured_3D(boundary_condition) end """ - boundary_condition_default_tree_1D(boundary_condition) + boundary_condition_default(mesh::TreeMesh1D, boundary_condition) Create a default boundary condition dictionary for tree meshes in 1D that use the standard boundary naming convention. @@ -135,14 +134,13 @@ This function applies the same boundary condition to all standard boundaries: # Returns - Named tuple mapping boundary names to the boundary condition """ -function boundary_condition_default_tree_1D(boundary_condition) - +function boundary_condition_default(mesh::TreeMesh1D, boundary_condition) return (x_neg = boundary_condition, x_pos = boundary_condition) end """ - boundary_condition_default_tree_2D(boundary_condition) + boundary_condition_default(mesh::TreeMesh2D, boundary_condition) Create a default boundary condition dictionary for tree meshes in 2D that use the standard boundary naming convention. @@ -156,7 +154,7 @@ This function applies the same boundary condition to all standard boundaries: # Returns - Named tuple mapping boundary names to the boundary condition """ -function boundary_condition_default_tree_2D(boundary_condition) +function boundary_condition_default(mesh::TreeMesh2D, boundary_condition) return (x_neg = boundary_condition, x_pos = boundary_condition, @@ -165,7 +163,7 @@ function boundary_condition_default_tree_2D(boundary_condition) end """ - boundary_condition_default_tree_3D(boundary_condition) + boundary_condition_default(mesh::TreeMesh3D, boundary_condition) Create a default boundary condition dictionary for tree meshes in 3D that use the standard boundary naming convention. @@ -181,7 +179,7 @@ This function applies the same boundary condition to all standard boundaries: # Returns - Named tuple mapping boundary names to the boundary condition """ -function boundary_condition_default_tree_3D(boundary_condition) +function boundary_condition_default(mesh::TreeMesh3D, boundary_condition) return (x_neg = boundary_condition, x_pos = boundary_condition, @@ -189,53 +187,4 @@ function boundary_condition_default_tree_3D(boundary_condition) y_pos = boundary_condition, z_neg = boundary_condition, z_pos = boundary_condition) -end - -""" - boundary_condition_default(mesh, boundary_condition) - -Create a default boundary condition dictionary based on the mesh type. -This function automatically determines the appropriate boundary condition -format based on the mesh type and applies the given boundary condition -to all boundaries. - -# Arguments -- `mesh`: The mesh object (TreeMesh, StructuredMesh, P4estMesh) -- `boundary_condition`: The boundary condition function to apply to all boundaries - -# Returns -- Dictionary or named tuple mapping boundary names to the boundary condition -""" -function boundary_condition_default(mesh, boundary_condition) - - if isa(mesh, TreeMesh) - ndims_mesh = ndims(mesh) - if ndims_mesh == 1 - return boundary_condition_default_tree_1D(boundary_condition) - elseif ndims_mesh == 2 - return boundary_condition_default_tree_2D(boundary_condition) - elseif ndims_mesh == 3 - return boundary_condition_default_tree_3D(boundary_condition) - end - elseif isa(mesh, StructuredMesh) - ndims_mesh = ndims(mesh) - if ndims_mesh == 1 - return boundary_condition_default_structured_1D(boundary_condition) - elseif ndims_mesh == 2 - return boundary_condition_default_structured_2D(boundary_condition) - elseif ndims_mesh == 3 - return boundary_condition_default_structured_3D(boundary_condition) - end - elseif isa(mesh, P4estMesh) - ndims_mesh = ndims(mesh) - if ndims_mesh == 2 - return boundary_condition_default_p4est_2D(boundary_condition) - elseif ndims_mesh == 3 - return boundary_condition_default_p4est_3D(boundary_condition) - end - - end - - error("No default boundary conditions available for mesh type $(typeof(mesh)). " * - "Please specify the boundary conditions manually.") end \ No newline at end of file From 402c84cc48c788b43638cd34ac01ffc092642861 Mon Sep 17 00:00:00 2001 From: vincmarks Date: Sun, 19 Oct 2025 19:44:57 +0200 Subject: [PATCH 32/45] formatting --- src/solvers/dgsem_unstructured/sort_boundary_conditions.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/solvers/dgsem_unstructured/sort_boundary_conditions.jl b/src/solvers/dgsem_unstructured/sort_boundary_conditions.jl index 3998f2dc7a5..d6cf6e1ce6d 100644 --- a/src/solvers/dgsem_unstructured/sort_boundary_conditions.jl +++ b/src/solvers/dgsem_unstructured/sort_boundary_conditions.jl @@ -124,4 +124,4 @@ end # @eval due to @muladd @eval Adapt.@adapt_structure(UnstructuredSortedBoundaryTypes) -end#@muladd +end # @muladd From eb3e2fc5ec67fb60b3d465a5f12f834accb57df6 Mon Sep 17 00:00:00 2001 From: vincmarks Date: Sun, 19 Oct 2025 20:09:50 +0200 Subject: [PATCH 33/45] Syncen with Trixi.jl --- .github/workflows/ci.yml | 6 +- Project.toml | 2 +- ...er_convergence_implicit_sparse_jacobian.jl | 15 +- src/solvers/dgsem_p4est/dg_2d_parabolic.jl | 381 ++++++++++-------- src/solvers/dgsem_p4est/dg_3d_parabolic.jl | 343 ++++++++-------- src/solvers/dgsem_tree/dg_1d_parabolic.jl | 136 ++++--- src/solvers/dgsem_tree/dg_2d_parabolic.jl | 220 ++++++---- src/solvers/dgsem_tree/dg_3d_parabolic.jl | 312 ++++++++------ test/test_parabolic_3d.jl | 14 +- test/test_structured_2d.jl | 2 +- test/test_threaded.jl | 24 +- 11 files changed, 835 insertions(+), 620 deletions(-) diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml index 9e36e47c7f9..9c46176bc6f 100644 --- a/.github/workflows/ci.yml +++ b/.github/workflows/ci.yml @@ -91,8 +91,12 @@ jobs: os: ubuntu-latest arch: x64 trixi_test: threaded_legacy + - version: '1.12' + os: ubuntu-latest + arch: x64 + trixi_test: threaded_legacy - version: '1.10' - os: macos-13 + os: macos-15-intel arch: x64 trixi_test: mpi - version: '1.10' diff --git a/Project.toml b/Project.toml index a7a0b3fc48f..550b6bce797 100644 --- a/Project.toml +++ b/Project.toml @@ -1,7 +1,7 @@ name = "Trixi" uuid = "a7f1ee26-1774-49b1-8366-f1abc58fbfcb" authors = ["Michael Schlottke-Lakemper ", "Gregor Gassner ", "Hendrik Ranocha ", "Andrew R. Winters ", "Jesse Chan ", "Andrés Rueda-Ramírez "] -version = "0.13.9-DEV" +version = "0.13.11-DEV" [deps] Accessors = "7d9f7c33-5ae7-4f3b-8dc6-eff91059b697" diff --git a/examples/structured_2d_dgsem/elixir_euler_convergence_implicit_sparse_jacobian.jl b/examples/structured_2d_dgsem/elixir_euler_convergence_implicit_sparse_jacobian.jl index 909d9b3e6ec..fd821029f33 100644 --- a/examples/structured_2d_dgsem/elixir_euler_convergence_implicit_sparse_jacobian.jl +++ b/examples/structured_2d_dgsem/elixir_euler_convergence_implicit_sparse_jacobian.jl @@ -78,23 +78,24 @@ coloring_vec = column_colors(coloring_result) ### sparsity-aware semidiscretization and ode ### # Semidiscretization for actual simulation -semi_float_type = SemidiscretizationHyperbolic(mesh, equations, - initial_condition_convergence_test, - solver, - source_terms = source_terms_convergence_test) +# Must be called 'semi' in order for the convergence test to run successfully +semi = SemidiscretizationHyperbolic(mesh, equations, + initial_condition_convergence_test, + solver, + source_terms = source_terms_convergence_test) # Supply Jacobian prototype and coloring vector to the semidiscretization -ode_jac_sparse = semidiscretize(semi_float_type, tspan, +ode_jac_sparse = semidiscretize(semi, tspan, jac_prototype = jac_prototype, colorvec = coloring_vec) -# using "dense" `ode = semidiscretize(semi_float_type, tspan)` +# using "dense" `ode = semidiscretize(semi, tspan)` # is essentially infeasible, even single step takes ages! ############################################################################### ### callbacks & solve ### summary_callback = SummaryCallback() -analysis_callback = AnalysisCallback(semi_float_type, interval = 50) +analysis_callback = AnalysisCallback(semi, interval = 50) alive_callback = AliveCallback(alive_interval = 3) # Note: No `stepsize_callback` due to implicit solver diff --git a/src/solvers/dgsem_p4est/dg_2d_parabolic.jl b/src/solvers/dgsem_p4est/dg_2d_parabolic.jl index 7d263b5fa2e..c14bdea33c6 100644 --- a/src/solvers/dgsem_p4est/dg_2d_parabolic.jl +++ b/src/solvers/dgsem_p4est/dg_2d_parabolic.jl @@ -164,57 +164,8 @@ function calc_gradient!(gradients, u_transformed, t, # Calculate volume integral @trixi_timeit timer() "volume integral" begin - (; derivative_dhat) = dg.basis - (; contravariant_vectors) = cache.elements - - @threaded for element in eachelement(dg, cache) - - # Calculate gradients with respect to reference coordinates in one element - for j in eachnode(dg), i in eachnode(dg) - u_node = get_node_vars(u_transformed, equations_parabolic, dg, i, j, - element) - - for ii in eachnode(dg) - multiply_add_to_node_vars!(gradients_x, derivative_dhat[ii, i], - u_node, equations_parabolic, dg, - ii, j, element) - end - - for jj in eachnode(dg) - multiply_add_to_node_vars!(gradients_y, derivative_dhat[jj, j], - u_node, equations_parabolic, dg, - i, jj, element) - end - end - - # now that the reference coordinate gradients are computed, transform them node-by-node to physical gradients - # using the contravariant vectors - for j in eachnode(dg), i in eachnode(dg) - Ja11, Ja12 = get_contravariant_vector(1, contravariant_vectors, - i, j, element) - Ja21, Ja22 = get_contravariant_vector(2, contravariant_vectors, - i, j, element) - - gradients_reference_1 = get_node_vars(gradients_x, equations_parabolic, - dg, i, j, element) - gradients_reference_2 = get_node_vars(gradients_y, equations_parabolic, - dg, i, j, element) - - # note that the contravariant vectors are transposed compared with computations of flux - # divergences in `calc_volume_integral!`. See - # https://github.com/trixi-framework/Trixi.jl/pull/1490#discussion_r1213345190 - # for a more detailed discussion. - gradient_x_node = Ja11 * gradients_reference_1 + - Ja21 * gradients_reference_2 - gradient_y_node = Ja12 * gradients_reference_1 + - Ja22 * gradients_reference_2 - - set_node_vars!(gradients_x, gradient_x_node, equations_parabolic, dg, - i, j, element) - set_node_vars!(gradients_y, gradient_y_node, equations_parabolic, dg, - i, j, element) - end - end + calc_gradient_volume_integral!(gradients, u_transformed, + mesh, equations_parabolic, dg, cache) end # Prolong solution to interfaces. @@ -241,9 +192,9 @@ function calc_gradient!(gradients, u_transformed, t, # Calculate boundary fluxes @trixi_timeit timer() "boundary flux" begin - calc_boundary_flux_gradients!(cache_parabolic, t, boundary_conditions_parabolic, - mesh, equations_parabolic, dg.surface_integral, - dg) + calc_gradient_boundary_flux!(cache_parabolic, t, boundary_conditions_parabolic, + mesh, equations_parabolic, dg.surface_integral, + dg) end # Prolong solution to mortars. This reuses the hyperbolic version of `prolong2mortars` @@ -264,113 +215,8 @@ function calc_gradient!(gradients, u_transformed, t, # Calculate surface integrals @trixi_timeit timer() "surface integral" begin - (; boundary_interpolation) = dg.basis - (; surface_flux_values) = cache_parabolic.elements - (; contravariant_vectors) = cache.elements - - # Access the factors only once before beginning the loop to increase performance. - # We also use explicit assignments instead of `+=` to let `@muladd` turn these - # into FMAs (see comment at the top of the file). - factor_1 = boundary_interpolation[1, 1] - factor_2 = boundary_interpolation[nnodes(dg), 2] - @threaded for element in eachelement(dg, cache) - for l in eachnode(dg) - for v in eachvariable(equations_parabolic) - - # Compute x-component of gradients - - # surface at -x - normal_direction_x, _ = get_normal_direction(1, - contravariant_vectors, - 1, l, element) - gradients_x[v, 1, l, element] = (gradients_x[v, 1, l, element] + - surface_flux_values[v, l, 1, - element] * - factor_1 * normal_direction_x) - - # surface at +x - normal_direction_x, _ = get_normal_direction(2, - contravariant_vectors, - nnodes(dg), l, element) - gradients_x[v, nnodes(dg), l, element] = (gradients_x[v, nnodes(dg), - l, - element] + - surface_flux_values[v, l, - 2, - element] * - factor_2 * - normal_direction_x) - - # surface at -y - normal_direction_x, _ = get_normal_direction(3, - contravariant_vectors, - l, 1, element) - gradients_x[v, l, 1, element] = (gradients_x[v, l, 1, element] + - surface_flux_values[v, l, 3, - element] * - factor_1 * normal_direction_x) - - # surface at +y - normal_direction_x, _ = get_normal_direction(4, - contravariant_vectors, - l, nnodes(dg), element) - gradients_x[v, l, nnodes(dg), element] = (gradients_x[v, l, - nnodes(dg), - element] + - surface_flux_values[v, l, - 4, - element] * - factor_2 * - normal_direction_x) - - # Compute y-component of gradients - - # surface at -x - _, normal_direction_y = get_normal_direction(1, - contravariant_vectors, - 1, l, element) - gradients_y[v, 1, l, element] = (gradients_y[v, 1, l, element] + - surface_flux_values[v, l, 1, - element] * - factor_1 * normal_direction_y) - - # surface at +x - _, normal_direction_y = get_normal_direction(2, - contravariant_vectors, - nnodes(dg), l, element) - gradients_y[v, nnodes(dg), l, element] = (gradients_y[v, nnodes(dg), - l, - element] + - surface_flux_values[v, l, - 2, - element] * - factor_2 * - normal_direction_y) - - # surface at -y - _, normal_direction_y = get_normal_direction(3, - contravariant_vectors, - l, 1, element) - gradients_y[v, l, 1, element] = (gradients_y[v, l, 1, element] + - surface_flux_values[v, l, 3, - element] * - factor_1 * normal_direction_y) - - # surface at +y - _, normal_direction_y = get_normal_direction(4, - contravariant_vectors, - l, nnodes(dg), element) - gradients_y[v, l, nnodes(dg), element] = (gradients_y[v, l, - nnodes(dg), - element] + - surface_flux_values[v, l, - 4, - element] * - factor_2 * - normal_direction_y) - end - end - end + calc_gradient_surface_integral!(gradients, mesh, equations_parabolic, + dg, cache, cache_parabolic) end # Apply Jacobian from mapping to reference element @@ -831,7 +677,7 @@ end return nothing end -# TODO: parabolic, finish implementing `calc_boundary_flux_gradients!` and `calc_boundary_flux_divergence!` +# TODO: parabolic, finish implementing `calc_gradient_boundary_flux!` and `calc_boundary_flux_divergence!` function prolong2boundaries!(cache_parabolic, flux_viscous, mesh::P4estMesh{2}, equations_parabolic::AbstractEquationsParabolic, @@ -872,17 +718,77 @@ function prolong2boundaries!(cache_parabolic, flux_viscous, return nothing end -function calc_boundary_flux_gradients!(cache, t, - boundary_condition::Union{BoundaryConditionPeriodic, - BoundaryConditionDoNothing}, - mesh::P4estMesh, equations, surface_integral, - dg::DG) +function calc_gradient_volume_integral!(gradients, u_transformed, + mesh::P4estMesh{2}, # for dispatch only + equations_parabolic::AbstractEquationsParabolic, + dg::DG, cache) + @unpack derivative_dhat = dg.basis + @unpack contravariant_vectors = cache.elements + gradients_x, gradients_y = gradients + + @threaded for element in eachelement(dg, cache) + # Calculate volume terms in one element, + # corresponds to `kernel` functions for the hyperbolic part of the flux + for j in eachnode(dg), i in eachnode(dg) + u_node = get_node_vars(u_transformed, equations_parabolic, dg, + i, j, element) + + for ii in eachnode(dg) + multiply_add_to_node_vars!(gradients_x, derivative_dhat[ii, i], + u_node, equations_parabolic, dg, + ii, j, element) + end + + for jj in eachnode(dg) + multiply_add_to_node_vars!(gradients_y, derivative_dhat[jj, j], + u_node, equations_parabolic, dg, + i, jj, element) + end + end + + # now that the reference coordinate gradients are computed, transform them node-by-node to physical gradients + # using the contravariant vectors + for j in eachnode(dg), i in eachnode(dg) + Ja11, Ja12 = get_contravariant_vector(1, contravariant_vectors, + i, j, element) + Ja21, Ja22 = get_contravariant_vector(2, contravariant_vectors, + i, j, element) + + gradients_reference_1 = get_node_vars(gradients_x, equations_parabolic, dg, + i, j, element) + gradients_reference_2 = get_node_vars(gradients_y, equations_parabolic, dg, + i, j, element) + + # note that the contravariant vectors are transposed compared with computations of flux + # divergences in `calc_volume_integral!`. See + # https://github.com/trixi-framework/Trixi.jl/pull/1490#discussion_r1213345190 + # for a more detailed discussion. + gradient_x_node = Ja11 * gradients_reference_1 + + Ja21 * gradients_reference_2 + gradient_y_node = Ja12 * gradients_reference_1 + + Ja22 * gradients_reference_2 + + set_node_vars!(gradients_x, gradient_x_node, equations_parabolic, dg, + i, j, element) + set_node_vars!(gradients_y, gradient_y_node, equations_parabolic, dg, + i, j, element) + end + end + + return nothing +end + +function calc_gradient_boundary_flux!(cache, t, + boundary_condition::Union{BoundaryConditionPeriodic, + BoundaryConditionDoNothing}, + mesh::P4estMesh, equations, surface_integral, + dg::DG) @assert isempty(eachboundary(dg, cache)) end # Function barrier for type stability -function calc_boundary_flux_gradients!(cache, t, boundary_conditions, mesh::P4estMesh, - equations, surface_integral, dg::DG) +function calc_gradient_boundary_flux!(cache, t, boundary_conditions, mesh::P4estMesh, + equations, surface_integral, dg::DG) (; boundary_condition_types, boundary_indices) = boundary_conditions calc_boundary_flux_by_type!(cache, t, boundary_condition_types, boundary_indices, @@ -996,6 +902,143 @@ function calc_boundary_flux!(cache, t, return nothing end +function calc_gradient_surface_integral!(gradients, + mesh::P4estMesh{2}, # for dispatch only + equations_parabolic::AbstractEquationsParabolic, + dg::DGSEM, cache, cache_parabolic) + @unpack boundary_interpolation = dg.basis + @unpack surface_flux_values = cache_parabolic.elements + @unpack contravariant_vectors = cache.elements + + gradients_x, gradients_y = gradients + + # Access the factors only once before beginning the loop to increase performance. + # We also use explicit assignments instead of `+=` to let `@muladd` turn these + # into FMAs (see comment at the top of the file). + factor_1 = boundary_interpolation[1, 1] + factor_2 = boundary_interpolation[nnodes(dg), 2] + @threaded for element in eachelement(dg, cache) + for l in eachnode(dg) + for v in eachvariable(equations_parabolic) + + # Compute x-component of gradients + + # surface at -x + normal_direction_x, _ = get_normal_direction(1, + contravariant_vectors, + 1, l, + element) + gradients_x[v, 1, l, element] = (gradients_x[v, + 1, l, + element] + + surface_flux_values[v, + l, 1, + element] * + factor_1 * normal_direction_x) + + # surface at +x + normal_direction_x, _ = get_normal_direction(2, + contravariant_vectors, + nnodes(dg), l, + element) + gradients_x[v, nnodes(dg), l, element] = (gradients_x[v, + nnodes(dg), l, + element] + + surface_flux_values[v, + l, 2, + element] * + factor_2 * + normal_direction_x) + + # surface at -y + normal_direction_x, _ = get_normal_direction(3, + contravariant_vectors, + l, 1, + element) + gradients_x[v, l, 1, element] = (gradients_x[v, + l, 1, + element] + + surface_flux_values[v, + l, 3, + element] * + factor_1 * normal_direction_x) + + # surface at +y + normal_direction_x, _ = get_normal_direction(4, + contravariant_vectors, + l, nnodes(dg), + element) + gradients_x[v, l, nnodes(dg), element] = (gradients_x[v, + l, nnodes(dg), + element] + + surface_flux_values[v, + l, 4, + element] * + factor_2 * + normal_direction_x) + + # Compute y-component of gradients + + # surface at -x + _, normal_direction_y = get_normal_direction(1, + contravariant_vectors, + 1, l, + element) + gradients_y[v, 1, l, element] = (gradients_y[v, + 1, l, + element] + + surface_flux_values[v, + l, 1, + element] * + factor_1 * normal_direction_y) + + # surface at +x + _, normal_direction_y = get_normal_direction(2, + contravariant_vectors, + nnodes(dg), l, + element) + gradients_y[v, nnodes(dg), l, element] = (gradients_y[v, + nnodes(dg), l, + element] + + surface_flux_values[v, + l, 2, + element] * + factor_2 * + normal_direction_y) + + # surface at -y + _, normal_direction_y = get_normal_direction(3, + contravariant_vectors, + l, 1, + element) + gradients_y[v, l, 1, element] = (gradients_y[v, + l, 1, + element] + + surface_flux_values[v, + l, 3, + element] * + factor_1 * normal_direction_y) + + # surface at +y + _, normal_direction_y = get_normal_direction(4, + contravariant_vectors, + l, nnodes(dg), + element) + gradients_y[v, l, nnodes(dg), element] = (gradients_y[v, + l, nnodes(dg), + element] + + surface_flux_values[v, + l, 4, + element] * + factor_2 * + normal_direction_y) + end + end + end + + return nothing +end + function apply_jacobian_parabolic!(du, mesh::P4estMesh{2}, equations::AbstractEquationsParabolic, dg::DG, cache) diff --git a/src/solvers/dgsem_p4est/dg_3d_parabolic.jl b/src/solvers/dgsem_p4est/dg_3d_parabolic.jl index 6703d3014de..bd76bd2487d 100644 --- a/src/solvers/dgsem_p4est/dg_3d_parabolic.jl +++ b/src/solvers/dgsem_p4est/dg_3d_parabolic.jl @@ -42,74 +42,8 @@ function calc_gradient!(gradients, u_transformed, t, # Calculate volume integral @trixi_timeit timer() "volume integral" begin - (; derivative_dhat) = dg.basis - (; contravariant_vectors) = cache.elements - - @threaded for element in eachelement(dg, cache) - - # Calculate gradients with respect to reference coordinates in one element - for k in eachnode(dg), j in eachnode(dg), i in eachnode(dg) - u_node = get_node_vars(u_transformed, equations_parabolic, dg, - i, j, k, element) - - for ii in eachnode(dg) - multiply_add_to_node_vars!(gradients_x, derivative_dhat[ii, i], - u_node, equations_parabolic, dg, - ii, j, k, element) - end - - for jj in eachnode(dg) - multiply_add_to_node_vars!(gradients_y, derivative_dhat[jj, j], - u_node, equations_parabolic, dg, - i, jj, k, element) - end - - for kk in eachnode(dg) - multiply_add_to_node_vars!(gradients_z, derivative_dhat[kk, k], - u_node, equations_parabolic, dg, - i, j, kk, element) - end - end - - # now that the reference coordinate gradients are computed, transform them node-by-node to physical gradients - # using the contravariant vectors - for k in eachnode(dg), j in eachnode(dg), i in eachnode(dg) - Ja11, Ja12, Ja13 = get_contravariant_vector(1, contravariant_vectors, - i, j, k, element) - Ja21, Ja22, Ja23 = get_contravariant_vector(2, contravariant_vectors, - i, j, k, element) - Ja31, Ja32, Ja33 = get_contravariant_vector(3, contravariant_vectors, - i, j, k, element) - - gradients_reference_1 = get_node_vars(gradients_x, equations_parabolic, - dg, i, j, k, element) - gradients_reference_2 = get_node_vars(gradients_y, equations_parabolic, - dg, i, j, k, element) - gradients_reference_3 = get_node_vars(gradients_z, equations_parabolic, - dg, i, j, k, element) - - # note that the contravariant vectors are transposed compared with computations of flux - # divergences in `calc_volume_integral!`. See - # https://github.com/trixi-framework/Trixi.jl/pull/1490#discussion_r1213345190 - # for a more detailed discussion. - gradient_x_node = Ja11 * gradients_reference_1 + - Ja21 * gradients_reference_2 + - Ja31 * gradients_reference_3 - gradient_y_node = Ja12 * gradients_reference_1 + - Ja22 * gradients_reference_2 + - Ja32 * gradients_reference_3 - gradient_z_node = Ja13 * gradients_reference_1 + - Ja23 * gradients_reference_2 + - Ja33 * gradients_reference_3 - - set_node_vars!(gradients_x, gradient_x_node, equations_parabolic, dg, - i, j, k, element) - set_node_vars!(gradients_y, gradient_y_node, equations_parabolic, dg, - i, j, k, element) - set_node_vars!(gradients_z, gradient_z_node, equations_parabolic, dg, - i, j, k, element) - end - end + calc_gradient_volume_integral!(gradients, u_transformed, + mesh, equations_parabolic, dg, cache) end # Prolong solution to interfaces @@ -135,9 +69,9 @@ function calc_gradient!(gradients, u_transformed, t, # Calculate boundary fluxes @trixi_timeit timer() "boundary flux" begin - calc_boundary_flux_gradients!(cache_parabolic, t, boundary_conditions_parabolic, - mesh, equations_parabolic, dg.surface_integral, - dg) + calc_gradient_boundary_flux!(cache_parabolic, t, boundary_conditions_parabolic, + mesh, equations_parabolic, dg.surface_integral, + dg) end # Prolong solution to mortars. These should reuse the hyperbolic version of `prolong2mortars` @@ -160,92 +94,8 @@ function calc_gradient!(gradients, u_transformed, t, # Calculate surface integrals @trixi_timeit timer() "surface integral" begin - (; boundary_interpolation) = dg.basis - (; surface_flux_values) = cache_parabolic.elements - (; contravariant_vectors) = cache.elements - - # Access the factors only once before beginning the loop to increase performance. - # We also use explicit assignments instead of `+=` to let `@muladd` turn these - # into FMAs (see comment at the top of the file). - factor_1 = boundary_interpolation[1, 1] - factor_2 = boundary_interpolation[nnodes(dg), 2] - @threaded for element in eachelement(dg, cache) - for l in eachnode(dg), m in eachnode(dg) - for v in eachvariable(equations_parabolic) - for dim in 1:3 - grad = gradients[dim] - # surface at -x - normal_direction = get_normal_direction(1, - contravariant_vectors, - 1, l, m, element) - grad[v, 1, l, m, element] = (grad[v, 1, l, m, element] + - surface_flux_values[v, l, m, 1, - element] * - factor_1 * normal_direction[dim]) - - # surface at +x - normal_direction = get_normal_direction(2, - contravariant_vectors, - nnodes(dg), l, m, - element) - grad[v, nnodes(dg), l, m, element] = (grad[v, nnodes(dg), l, m, - element] + - surface_flux_values[v, l, - m, - 2, - element] * - factor_2 * - normal_direction[dim]) - - # surface at -y - normal_direction = get_normal_direction(3, - contravariant_vectors, - l, m, 1, element) - grad[v, l, 1, m, element] = (grad[v, l, 1, m, element] + - surface_flux_values[v, l, m, 3, - element] * - factor_1 * normal_direction[dim]) - - # surface at +y - normal_direction = get_normal_direction(4, - contravariant_vectors, - l, nnodes(dg), m, - element) - grad[v, l, nnodes(dg), m, element] = (grad[v, l, nnodes(dg), m, - element] + - surface_flux_values[v, l, - m, - 4, - element] * - factor_2 * - normal_direction[dim]) - - # surface at -z - normal_direction = get_normal_direction(5, - contravariant_vectors, - l, m, 1, element) - grad[v, l, m, 1, element] = (grad[v, l, m, 1, element] + - surface_flux_values[v, l, m, 5, - element] * - factor_1 * normal_direction[dim]) - - # surface at +z - normal_direction = get_normal_direction(6, - contravariant_vectors, - l, m, nnodes(dg), - element) - grad[v, l, m, nnodes(dg), element] = (grad[v, l, m, nnodes(dg), - element] + - surface_flux_values[v, l, - m, - 6, - element] * - factor_2 * - normal_direction[dim]) - end - end - end - end + calc_gradient_surface_integral!(gradients, mesh, equations_parabolic, + dg, cache, cache_parabolic) end # Apply Jacobian from mapping to reference element @@ -868,7 +718,84 @@ end return nothing end -# TODO: parabolic, finish implementing `calc_boundary_flux_gradients!` and `calc_boundary_flux_divergence!` +function calc_gradient_volume_integral!(gradients, u_transformed, + mesh::P4estMesh{3}, # for dispatch only + equations_parabolic::AbstractEquationsParabolic, + dg::DG, cache) + @unpack derivative_dhat = dg.basis + @unpack contravariant_vectors = cache.elements + gradients_x, gradients_y, gradients_z = gradients + + @threaded for element in eachelement(dg, cache) + # Calculate volume terms in one element, + # corresponds to `kernel` functions for the hyperbolic part of the flux + for k in eachnode(dg), j in eachnode(dg), i in eachnode(dg) + u_node = get_node_vars(u_transformed, equations_parabolic, dg, + i, j, k, element) + + for ii in eachnode(dg) + multiply_add_to_node_vars!(gradients_x, derivative_dhat[ii, i], + u_node, equations_parabolic, dg, + ii, j, k, element) + end + + for jj in eachnode(dg) + multiply_add_to_node_vars!(gradients_y, derivative_dhat[jj, j], + u_node, equations_parabolic, dg, + i, jj, k, element) + end + + for kk in eachnode(dg) + multiply_add_to_node_vars!(gradients_z, derivative_dhat[kk, k], + u_node, equations_parabolic, dg, + i, j, kk, element) + end + end + + # now that the reference coordinate gradients are computed, transform them node-by-node to physical gradients + # using the contravariant vectors + for k in eachnode(dg), j in eachnode(dg), i in eachnode(dg) + Ja11, Ja12, Ja13 = get_contravariant_vector(1, contravariant_vectors, + i, j, k, element) + Ja21, Ja22, Ja23 = get_contravariant_vector(2, contravariant_vectors, + i, j, k, element) + Ja31, Ja32, Ja33 = get_contravariant_vector(3, contravariant_vectors, + i, j, k, element) + + gradients_reference_1 = get_node_vars(gradients_x, equations_parabolic, dg, + i, j, k, element) + gradients_reference_2 = get_node_vars(gradients_y, equations_parabolic, dg, + i, j, k, element) + gradients_reference_3 = get_node_vars(gradients_z, equations_parabolic, dg, + i, j, k, element) + + # note that the contravariant vectors are transposed compared with computations of flux + # divergences in `calc_volume_integral!`. See + # https://github.com/trixi-framework/Trixi.jl/pull/1490#discussion_r1213345190 + # for a more detailed discussion. + gradient_x_node = Ja11 * gradients_reference_1 + + Ja21 * gradients_reference_2 + + Ja31 * gradients_reference_3 + gradient_y_node = Ja12 * gradients_reference_1 + + Ja22 * gradients_reference_2 + + Ja32 * gradients_reference_3 + gradient_z_node = Ja13 * gradients_reference_1 + + Ja23 * gradients_reference_2 + + Ja33 * gradients_reference_3 + + set_node_vars!(gradients_x, gradient_x_node, equations_parabolic, dg, + i, j, k, element) + set_node_vars!(gradients_y, gradient_y_node, equations_parabolic, dg, + i, j, k, element) + set_node_vars!(gradients_z, gradient_z_node, equations_parabolic, dg, + i, j, k, element) + end + end + + return nothing +end + +# TODO: parabolic, finish implementing `calc_gradient_boundary_flux!` and `calc_boundary_flux_divergence!` function prolong2boundaries!(cache_parabolic, flux_viscous, mesh::P4estMesh{3}, equations_parabolic::AbstractEquationsParabolic, @@ -1002,6 +929,106 @@ function calc_boundary_flux!(cache, t, return nothing end +function calc_gradient_surface_integral!(gradients, + mesh::P4estMesh{3}, # for dispatch only + equations_parabolic::AbstractEquationsParabolic, + dg::DGSEM, cache, cache_parabolic) + @unpack boundary_interpolation = dg.basis + @unpack surface_flux_values = cache_parabolic.elements + @unpack contravariant_vectors = cache.elements + + # Access the factors only once before beginning the loop to increase performance. + # We also use explicit assignments instead of `+=` to let `@muladd` turn these + # into FMAs (see comment at the top of the file). + factor_1 = boundary_interpolation[1, 1] + factor_2 = boundary_interpolation[nnodes(dg), 2] + @threaded for element in eachelement(dg, cache) + for l in eachnode(dg), m in eachnode(dg) + for v in eachvariable(equations_parabolic) + for dim in 1:3 + grad = gradients[dim] + # surface at -x + normal_direction = get_normal_direction(1, contravariant_vectors, + 1, l, m, + element) + grad[v, 1, l, m, element] = (grad[v, + 1, l, m, + element] + + surface_flux_values[v, + l, m, 1, + element] * + factor_1 * normal_direction[dim]) + + # surface at +x + normal_direction = get_normal_direction(2, contravariant_vectors, + nnodes(dg), l, m, + element) + grad[v, nnodes(dg), l, m, element] = (grad[v, + nnodes(dg), l, m, + element] + + surface_flux_values[v, + l, m, 2, + element] * + factor_2 * + normal_direction[dim]) + + # surface at -y + normal_direction = get_normal_direction(3, contravariant_vectors, + l, 1, m, + element) + grad[v, l, 1, m, element] = (grad[v, + l, 1, m, + element] + + surface_flux_values[v, + l, m, 3, + element] * + factor_1 * normal_direction[dim]) + + # surface at +y + normal_direction = get_normal_direction(4, contravariant_vectors, + l, nnodes(dg), m, + element) + grad[v, l, nnodes(dg), m, element] = (grad[v, + l, nnodes(dg), m, + element] + + surface_flux_values[v, + l, m, 4, + element] * + factor_2 * + normal_direction[dim]) + + # surface at -z + normal_direction = get_normal_direction(5, contravariant_vectors, + l, m, 1, + element) + grad[v, l, m, 1, element] = (grad[v, + l, m, 1, + element] + + surface_flux_values[v, + l, m, 5, + element] * + factor_1 * normal_direction[dim]) + + # surface at +z + normal_direction = get_normal_direction(6, contravariant_vectors, + l, m, nnodes(dg), + element) + grad[v, l, m, nnodes(dg), element] = (grad[v, + l, m, nnodes(dg), + element] + + surface_flux_values[v, + l, m, 6, + element] * + factor_2 * + normal_direction[dim]) + end + end + end + end + + return nothing +end + function apply_jacobian_parabolic!(du, mesh::P4estMesh{3}, equations::AbstractEquationsParabolic, dg::DG, cache) diff --git a/src/solvers/dgsem_tree/dg_1d_parabolic.jl b/src/solvers/dgsem_tree/dg_1d_parabolic.jl index 06a6a4488ec..4cbc1afe3a0 100644 --- a/src/solvers/dgsem_tree/dg_1d_parabolic.jl +++ b/src/solvers/dgsem_tree/dg_1d_parabolic.jl @@ -255,11 +255,11 @@ function calc_viscous_fluxes!(flux_viscous, gradients, u_transformed, mesh::Tree return nothing end -function calc_boundary_flux_gradients!(cache, t, - boundary_conditions_parabolic::BoundaryConditionPeriodic, - mesh::TreeMesh{1}, - equations_parabolic::AbstractEquationsParabolic, - surface_integral, dg::DG) +function calc_gradient_boundary_flux!(cache, t, + boundary_conditions_parabolic::BoundaryConditionPeriodic, + mesh::TreeMesh{1}, + equations_parabolic::AbstractEquationsParabolic, + surface_integral, dg::DG) return nothing end @@ -271,11 +271,11 @@ function calc_boundary_flux_divergence!(cache, t, return nothing end -function calc_boundary_flux_gradients!(cache, t, - boundary_conditions_parabolic::NamedTuple, - mesh::TreeMesh{1}, - equations_parabolic::AbstractEquationsParabolic, - surface_integral, dg::DG) +function calc_gradient_boundary_flux!(cache, t, + boundary_conditions_parabolic::NamedTuple, + mesh::TreeMesh{1}, # for dispatch only + equations_parabolic::AbstractEquationsParabolic, + surface_integral, dg::DG) @unpack surface_flux_values = cache.elements @unpack n_boundaries_per_direction = cache.boundaries @@ -284,12 +284,12 @@ function calc_boundary_flux_gradients!(cache, t, firsts = lasts - n_boundaries_per_direction .+ 1 # Calc boundary fluxes in each direction - calc_boundary_flux_by_direction_gradient!(surface_flux_values, t, + calc_gradient_boundary_flux_by_direction!(surface_flux_values, t, boundary_conditions_parabolic[1], equations_parabolic, surface_integral, dg, cache, 1, firsts[1], lasts[1]) - calc_boundary_flux_by_direction_gradient!(surface_flux_values, t, + calc_gradient_boundary_flux_by_direction!(surface_flux_values, t, boundary_conditions_parabolic[2], equations_parabolic, surface_integral, dg, cache, @@ -298,7 +298,7 @@ function calc_boundary_flux_gradients!(cache, t, return nothing end -function calc_boundary_flux_by_direction_gradient!(surface_flux_values::AbstractArray{<:Any, +function calc_gradient_boundary_flux_by_direction!(surface_flux_values::AbstractArray{<:Any, 3}, t, boundary_condition, @@ -413,6 +413,30 @@ function calc_boundary_flux_by_direction_divergence!(surface_flux_values::Abstra return nothing end +function calc_gradient_volume_integral!(gradients, u_transformed, + mesh::TreeMesh{1}, # for dispatch only + equations_parabolic::AbstractEquationsParabolic, + dg::DGSEM, cache) + @unpack derivative_dhat = dg.basis + + @threaded for element in eachelement(dg, cache) + # Calculate volume terms in one element, + # corresponds to `kernel` functions for the hyperbolic part of the flux + for i in eachnode(dg) + u_node = get_node_vars(u_transformed, equations_parabolic, dg, + i, element) + + for ii in eachnode(dg) + multiply_add_to_node_vars!(gradients, derivative_dhat[ii, i], + u_node, equations_parabolic, dg, + ii, element) + end + end + end + + return nothing +end + function calc_gradient_interface_flux!(surface_flux_values, mesh::TreeMesh{1}, equations_parabolic, dg::DG, parabolic_scheme, @@ -446,6 +470,36 @@ function calc_gradient_interface_flux!(surface_flux_values, return nothing end +function calc_gradient_surface_integral!(gradients, + mesh::TreeMesh{1}, # for dispatch only + equations_parabolic::AbstractEquationsParabolic, + dg::DGSEM, cache, cache_parabolic) + @unpack boundary_interpolation = dg.basis + @unpack surface_flux_values = cache_parabolic.elements + + # Note that all fluxes have been computed with outward-pointing normal vectors. + # Access the factors only once before beginning the loop to increase performance. + # We also use explicit assignments instead of `+=` to let `@muladd` turn these + # into FMAs (see comment at the top of the file). + factor_1 = boundary_interpolation[1, 1] + factor_2 = boundary_interpolation[nnodes(dg), 2] + @threaded for element in eachelement(dg, cache) + for v in eachvariable(equations_parabolic) + # surface at -x + gradients[v, 1, element] = (gradients[v, 1, element] - + surface_flux_values[v, 1, element] * + factor_1) + + # surface at +x + gradients[v, nnodes(dg), element] = (gradients[v, nnodes(dg), element] + + surface_flux_values[v, 2, element] * + factor_2) + end + end + + return nothing +end + # Calculate the gradient of the transformed variables function calc_gradient!(gradients, u_transformed, t, mesh::TreeMesh{1}, equations_parabolic, boundary_conditions_parabolic, @@ -458,21 +512,8 @@ function calc_gradient!(gradients, u_transformed, t, mesh::TreeMesh{1}, # Calculate volume integral @trixi_timeit timer() "volume integral" begin - @unpack derivative_dhat = dg.basis - @threaded for element in eachelement(dg, cache) - - # Calculate volume terms in one element - for i in eachnode(dg) - u_node = get_node_vars(u_transformed, equations_parabolic, dg, - i, element) - - for ii in eachnode(dg) - multiply_add_to_node_vars!(gradients, derivative_dhat[ii, i], - u_node, equations_parabolic, dg, - ii, element) - end - end - end + calc_gradient_volume_integral!(gradients, u_transformed, + mesh, equations_parabolic, dg, cache) end # Prolong solution to interfaces @@ -497,39 +538,18 @@ function calc_gradient!(gradients, u_transformed, t, mesh::TreeMesh{1}, dg) # Calculate boundary fluxes - @trixi_timeit timer() "boundary flux" calc_boundary_flux_gradients!(cache_parabolic, - t, - boundary_conditions_parabolic, - mesh, - equations_parabolic, - dg.surface_integral, - dg) + @trixi_timeit timer() "boundary flux" calc_gradient_boundary_flux!(cache_parabolic, + t, + boundary_conditions_parabolic, + mesh, + equations_parabolic, + dg.surface_integral, + dg) # Calculate surface integrals @trixi_timeit timer() "surface integral" begin - @unpack boundary_interpolation = dg.basis - @unpack surface_flux_values = cache_parabolic.elements - - # Note that all fluxes have been computed with outward-pointing normal vectors. - # Access the factors only once before beginning the loop to increase performance. - # We also use explicit assignments instead of `+=` to let `@muladd` turn these - # into FMAs (see comment at the top of the file). - factor_1 = boundary_interpolation[1, 1] - factor_2 = boundary_interpolation[nnodes(dg), 2] - @threaded for element in eachelement(dg, cache) - for v in eachvariable(equations_parabolic) - # surface at -x - gradients[v, 1, element] = (gradients[v, 1, element] - - surface_flux_values[v, 1, element] * - factor_1) - - # surface at +x - gradients[v, nnodes(dg), element] = (gradients[v, nnodes(dg), element] + - surface_flux_values[v, 2, - element] * - factor_2) - end - end + calc_gradient_surface_integral!(gradients, mesh, equations_parabolic, + dg, cache, cache_parabolic) end # Apply Jacobian from mapping to reference element diff --git a/src/solvers/dgsem_tree/dg_2d_parabolic.jl b/src/solvers/dgsem_tree/dg_2d_parabolic.jl index 232e13de88b..1488778d534 100644 --- a/src/solvers/dgsem_tree/dg_2d_parabolic.jl +++ b/src/solvers/dgsem_tree/dg_2d_parabolic.jl @@ -343,11 +343,11 @@ function get_unsigned_normal_vector_2d(direction) end end -function calc_boundary_flux_gradients!(cache, t, - boundary_conditions_parabolic::BoundaryConditionPeriodic, - mesh::Union{TreeMesh{2}, P4estMesh{2}}, - equations_parabolic::AbstractEquationsParabolic, - surface_integral, dg::DG) +function calc_gradient_boundary_flux!(cache, t, + boundary_conditions_parabolic::BoundaryConditionPeriodic, + mesh::Union{TreeMesh{2}, P4estMesh{2}}, + equations_parabolic::AbstractEquationsParabolic, + surface_integral, dg::DG) return nothing end @@ -359,11 +359,11 @@ function calc_boundary_flux_divergence!(cache, t, return nothing end -function calc_boundary_flux_gradients!(cache, t, - boundary_conditions_parabolic::NamedTuple, - mesh::TreeMesh{2}, - equations_parabolic::AbstractEquationsParabolic, - surface_integral, dg::DG) +function calc_gradient_boundary_flux!(cache, t, + boundary_conditions_parabolic::NamedTuple, + mesh::TreeMesh{2}, # for dispatch only + equations_parabolic::AbstractEquationsParabolic, + surface_integral, dg::DG) @unpack surface_flux_values = cache.elements @unpack n_boundaries_per_direction = cache.boundaries @@ -372,22 +372,22 @@ function calc_boundary_flux_gradients!(cache, t, firsts = lasts - n_boundaries_per_direction .+ 1 # Calc boundary fluxes in each direction - calc_boundary_flux_by_direction_gradient!(surface_flux_values, t, + calc_gradient_boundary_flux_by_direction!(surface_flux_values, t, boundary_conditions_parabolic[1], equations_parabolic, surface_integral, dg, cache, 1, firsts[1], lasts[1]) - calc_boundary_flux_by_direction_gradient!(surface_flux_values, t, + calc_gradient_boundary_flux_by_direction!(surface_flux_values, t, boundary_conditions_parabolic[2], equations_parabolic, surface_integral, dg, cache, 2, firsts[2], lasts[2]) - calc_boundary_flux_by_direction_gradient!(surface_flux_values, t, + calc_gradient_boundary_flux_by_direction!(surface_flux_values, t, boundary_conditions_parabolic[3], equations_parabolic, surface_integral, dg, cache, 3, firsts[3], lasts[3]) - calc_boundary_flux_by_direction_gradient!(surface_flux_values, t, + calc_gradient_boundary_flux_by_direction!(surface_flux_values, t, boundary_conditions_parabolic[4], equations_parabolic, surface_integral, dg, cache, @@ -396,7 +396,7 @@ function calc_boundary_flux_gradients!(cache, t, return nothing end -function calc_boundary_flux_by_direction_gradient!(surface_flux_values::AbstractArray{<:Any, +function calc_gradient_boundary_flux_by_direction!(surface_flux_values::AbstractArray{<:Any, 4}, t, boundary_condition, @@ -548,9 +548,13 @@ function prolong2mortars!(cache, flux_viscous::Vector{Array{uEltype, 4}}, # L2 mortars in x-direction for l in eachnode(dg) for v in eachvariable(equations_parabolic) - cache.mortars.u_upper[2, v, l, mortar] = flux_viscous_x[v, 1, l, + cache.mortars.u_upper[2, v, l, mortar] = flux_viscous_x[v, + 1, + l, upper_element] - cache.mortars.u_lower[2, v, l, mortar] = flux_viscous_x[v, 1, l, + cache.mortars.u_lower[2, v, l, mortar] = flux_viscous_x[v, + 1, + l, lower_element] end end @@ -558,9 +562,13 @@ function prolong2mortars!(cache, flux_viscous::Vector{Array{uEltype, 4}}, # L2 mortars in y-direction for l in eachnode(dg) for v in eachvariable(equations_parabolic) - cache.mortars.u_upper[2, v, l, mortar] = flux_viscous_y[v, l, 1, + cache.mortars.u_upper[2, v, l, mortar] = flux_viscous_y[v, + l, + 1, upper_element] - cache.mortars.u_lower[2, v, l, mortar] = flux_viscous_y[v, l, 1, + cache.mortars.u_lower[2, v, l, mortar] = flux_viscous_y[v, + l, + 1, lower_element] end end @@ -584,10 +592,12 @@ function prolong2mortars!(cache, flux_viscous::Vector{Array{uEltype, 4}}, # L2 mortars in y-direction for l in eachnode(dg) for v in eachvariable(equations_parabolic) - cache.mortars.u_upper[1, v, l, mortar] = flux_viscous_y[v, l, + cache.mortars.u_upper[1, v, l, mortar] = flux_viscous_y[v, + l, nnodes(dg), upper_element] - cache.mortars.u_lower[1, v, l, mortar] = flux_viscous_y[v, l, + cache.mortars.u_lower[1, v, l, mortar] = flux_viscous_y[v, + l, nnodes(dg), lower_element] end @@ -758,6 +768,37 @@ end return nothing end +function calc_gradient_volume_integral!(gradients, u_transformed, + mesh::TreeMesh{2}, # for dispatch only + equations_parabolic::AbstractEquationsParabolic, + dg::DGSEM, cache) + @unpack derivative_dhat = dg.basis + gradients_x, gradients_y = gradients + + @threaded for element in eachelement(dg, cache) + # Calculate volume terms in one element, + # corresponds to `kernel` functions for the hyperbolic part of the flux + for j in eachnode(dg), i in eachnode(dg) + u_node = get_node_vars(u_transformed, equations_parabolic, dg, + i, j, element) + + for ii in eachnode(dg) + multiply_add_to_node_vars!(gradients_x, derivative_dhat[ii, i], + u_node, equations_parabolic, dg, + ii, j, element) + end + + for jj in eachnode(dg) + multiply_add_to_node_vars!(gradients_y, derivative_dhat[jj, j], + u_node, equations_parabolic, dg, + i, jj, element) + end + end + end + + return nothing +end + function calc_gradient_interface_flux!(surface_flux_values, mesh::TreeMesh{2}, equations, dg::DG, parabolic_scheme, @@ -795,6 +836,67 @@ function calc_gradient_interface_flux!(surface_flux_values, return nothing end +function calc_gradient_surface_integral!(gradients, + mesh::TreeMesh{2}, # for dispatch only + equations_parabolic::AbstractEquationsParabolic, + dg::DGSEM, cache, cache_parabolic) + @unpack boundary_interpolation = dg.basis + @unpack surface_flux_values = cache_parabolic.elements + + gradients_x, gradients_y = gradients + + # Note that all fluxes have been computed with outward-pointing normal vectors. + # Access the factors only once before beginning the loop to increase performance. + # We also use explicit assignments instead of `+=` to let `@muladd` turn these + # into FMAs (see comment at the top of the file). + factor_1 = boundary_interpolation[1, 1] + factor_2 = boundary_interpolation[nnodes(dg), 2] + + @threaded for element in eachelement(dg, cache) + for l in eachnode(dg) + for v in eachvariable(equations_parabolic) + # surface at -x + gradients_x[v, 1, l, element] = (gradients_x[v, + 1, l, + element] - + surface_flux_values[v, + l, 1, + element] * + factor_1) + + # surface at +x + gradients_x[v, nnodes(dg), l, element] = (gradients_x[v, + nnodes(dg), l, + element] + + surface_flux_values[v, + l, 2, + element] * + factor_2) + + # surface at -y + gradients_y[v, l, 1, element] = (gradients_y[v, + l, 1, + element] - + surface_flux_values[v, + l, 3, + element] * + factor_1) + + # surface at +y + gradients_y[v, l, nnodes(dg), element] = (gradients_y[v, + l, nnodes(dg), + element] + + surface_flux_values[v, + l, 4, + element] * + factor_2) + end + end + end + + return nothing +end + # Calculate the gradient of the transformed variables function calc_gradient!(gradients, u_transformed, t, mesh::TreeMesh{2}, equations_parabolic, @@ -810,27 +912,8 @@ function calc_gradient!(gradients, u_transformed, t, # Calculate volume integral @trixi_timeit timer() "volume integral" begin - @unpack derivative_dhat = dg.basis - @threaded for element in eachelement(dg, cache) - - # Calculate volume terms in one element - for j in eachnode(dg), i in eachnode(dg) - u_node = get_node_vars(u_transformed, equations_parabolic, dg, - i, j, element) - - for ii in eachnode(dg) - multiply_add_to_node_vars!(gradients_x, derivative_dhat[ii, i], - u_node, equations_parabolic, dg, - ii, j, element) - end - - for jj in eachnode(dg) - multiply_add_to_node_vars!(gradients_y, derivative_dhat[jj, j], - u_node, equations_parabolic, dg, - i, jj, element) - end - end - end + calc_gradient_volume_integral!(gradients, u_transformed, + mesh, equations_parabolic, dg, cache) end # Prolong solution to interfaces @@ -855,10 +938,10 @@ function calc_gradient!(gradients, u_transformed, t, # Calculate boundary fluxes @trixi_timeit timer() "boundary flux" begin - calc_boundary_flux_gradients!(cache_parabolic, t, - boundary_conditions_parabolic, mesh, - equations_parabolic, - dg.surface_integral, dg) + calc_gradient_boundary_flux!(cache_parabolic, t, + boundary_conditions_parabolic, mesh, + equations_parabolic, + dg.surface_integral, dg) end # Prolong solution to mortars @@ -877,49 +960,8 @@ function calc_gradient!(gradients, u_transformed, t, # Calculate surface integrals @trixi_timeit timer() "surface integral" begin - @unpack boundary_interpolation = dg.basis - @unpack surface_flux_values = cache_parabolic.elements - - # Note that all fluxes have been computed with outward-pointing normal vectors. - # Access the factors only once before beginning the loop to increase performance. - # We also use explicit assignments instead of `+=` to let `@muladd` turn these - # into FMAs (see comment at the top of the file). - factor_1 = boundary_interpolation[1, 1] - factor_2 = boundary_interpolation[nnodes(dg), 2] - @threaded for element in eachelement(dg, cache) - for l in eachnode(dg) - for v in eachvariable(equations_parabolic) - # surface at -x - gradients_x[v, 1, l, element] = (gradients_x[v, 1, l, element] - - surface_flux_values[v, l, 1, - element] * - factor_1) - - # surface at +x - gradients_x[v, nnodes(dg), l, element] = (gradients_x[v, nnodes(dg), - l, element] + - surface_flux_values[v, l, - 2, - element] * - factor_2) - - # surface at -y - gradients_y[v, l, 1, element] = (gradients_y[v, l, 1, element] - - surface_flux_values[v, l, 3, - element] * - factor_1) - - # surface at +y - gradients_y[v, l, nnodes(dg), element] = (gradients_y[v, l, - nnodes(dg), - element] + - surface_flux_values[v, l, - 4, - element] * - factor_2) - end - end - end + calc_gradient_surface_integral!(gradients, mesh, equations_parabolic, + dg, cache, cache_parabolic) end # Apply Jacobian from mapping to reference element diff --git a/src/solvers/dgsem_tree/dg_3d_parabolic.jl b/src/solvers/dgsem_tree/dg_3d_parabolic.jl index a39d704199d..80fd41ebace 100644 --- a/src/solvers/dgsem_tree/dg_3d_parabolic.jl +++ b/src/solvers/dgsem_tree/dg_3d_parabolic.jl @@ -180,14 +180,20 @@ function prolong2boundaries!(cache_parabolic, flux_viscous, for k in eachnode(dg), j in eachnode(dg), v in eachvariable(equations_parabolic) # OBS! `boundaries_u` stores the interpolated *fluxes* and *not the solution*! - boundaries_u[1, v, j, k, boundary] = flux_viscous_x[v, nnodes(dg), - j, k, element] + boundaries_u[1, v, j, k, boundary] = flux_viscous_x[v, + nnodes(dg), + j, + k, + element] end else # Element in +x direction of boundary for k in eachnode(dg), j in eachnode(dg), v in eachvariable(equations_parabolic) # OBS! `boundaries_u` stores the interpolated *fluxes* and *not the solution*! - boundaries_u[2, v, j, k, boundary] = flux_viscous_x[v, 1, j, k, + boundaries_u[2, v, j, k, boundary] = flux_viscous_x[v, + 1, + j, + k, element] end end @@ -198,8 +204,10 @@ function prolong2boundaries!(cache_parabolic, flux_viscous, for k in eachnode(dg), i in eachnode(dg), v in eachvariable(equations_parabolic) # OBS! `boundaries_u` stores the interpolated *fluxes* and *not the solution*! - boundaries_u[1, v, i, k, boundary] = flux_viscous_y[v, i, - nnodes(dg), k, + boundaries_u[1, v, i, k, boundary] = flux_viscous_y[v, + i, + nnodes(dg), + k, element] end else @@ -207,7 +215,10 @@ function prolong2boundaries!(cache_parabolic, flux_viscous, for k in eachnode(dg), i in eachnode(dg), v in eachvariable(equations_parabolic) # OBS! `boundaries_u` stores the interpolated *fluxes* and *not the solution*! - boundaries_u[2, v, i, k, boundary] = flux_viscous_y[v, i, 1, k, + boundaries_u[2, v, i, k, boundary] = flux_viscous_y[v, + i, + 1, + k, element] end end @@ -218,7 +229,9 @@ function prolong2boundaries!(cache_parabolic, flux_viscous, for j in eachnode(dg), i in eachnode(dg), v in eachvariable(equations_parabolic) # OBS! `boundaries_u` stores the interpolated *fluxes* and *not the solution*! - boundaries_u[1, v, i, j, boundary] = flux_viscous_z[v, i, j, + boundaries_u[1, v, i, j, boundary] = flux_viscous_z[v, + i, + j, nnodes(dg), element] end @@ -227,7 +240,10 @@ function prolong2boundaries!(cache_parabolic, flux_viscous, for j in eachnode(dg), i in eachnode(dg), v in eachvariable(equations_parabolic) # OBS! `boundaries_u` stores the interpolated *fluxes* and *not the solution*! - boundaries_u[2, v, i, j, boundary] = flux_viscous_z[v, i, j, 1, + boundaries_u[2, v, i, j, boundary] = flux_viscous_z[v, + i, + j, + 1, element] end end @@ -293,11 +309,11 @@ function get_unsigned_normal_vector_3d(direction) end end -function calc_boundary_flux_gradients!(cache, t, - boundary_conditions_parabolic::BoundaryConditionPeriodic, - mesh::Union{TreeMesh{3}, P4estMesh{3}}, - equations_parabolic::AbstractEquationsParabolic, - surface_integral, dg::DG) +function calc_gradient_boundary_flux!(cache, t, + boundary_conditions_parabolic::BoundaryConditionPeriodic, + mesh::Union{TreeMesh{3}, P4estMesh{3}}, + equations_parabolic::AbstractEquationsParabolic, + surface_integral, dg::DG) return nothing end @@ -309,11 +325,11 @@ function calc_boundary_flux_divergence!(cache, t, return nothing end -function calc_boundary_flux_gradients!(cache, t, - boundary_conditions_parabolic::NamedTuple, - mesh::TreeMesh{3}, - equations_parabolic::AbstractEquationsParabolic, - surface_integral, dg::DG) +function calc_gradient_boundary_flux!(cache, t, + boundary_conditions_parabolic::NamedTuple, + mesh::TreeMesh{3}, # for dispatch only + equations_parabolic::AbstractEquationsParabolic, + surface_integral, dg::DG) @unpack surface_flux_values = cache.elements @unpack n_boundaries_per_direction = cache.boundaries @@ -322,32 +338,32 @@ function calc_boundary_flux_gradients!(cache, t, firsts = lasts - n_boundaries_per_direction .+ 1 # Calc boundary fluxes in each direction - calc_boundary_flux_by_direction_gradient!(surface_flux_values, t, + calc_gradient_boundary_flux_by_direction!(surface_flux_values, t, boundary_conditions_parabolic[1], equations_parabolic, surface_integral, dg, cache, 1, firsts[1], lasts[1]) - calc_boundary_flux_by_direction_gradient!(surface_flux_values, t, + calc_gradient_boundary_flux_by_direction!(surface_flux_values, t, boundary_conditions_parabolic[2], equations_parabolic, surface_integral, dg, cache, 2, firsts[2], lasts[2]) - calc_boundary_flux_by_direction_gradient!(surface_flux_values, t, + calc_gradient_boundary_flux_by_direction!(surface_flux_values, t, boundary_conditions_parabolic[3], equations_parabolic, surface_integral, dg, cache, 3, firsts[3], lasts[3]) - calc_boundary_flux_by_direction_gradient!(surface_flux_values, t, + calc_gradient_boundary_flux_by_direction!(surface_flux_values, t, boundary_conditions_parabolic[4], equations_parabolic, surface_integral, dg, cache, 4, firsts[4], lasts[4]) - calc_boundary_flux_by_direction_gradient!(surface_flux_values, t, + calc_gradient_boundary_flux_by_direction!(surface_flux_values, t, boundary_conditions_parabolic[5], equations_parabolic, surface_integral, dg, cache, 5, firsts[5], lasts[5]) - calc_boundary_flux_by_direction_gradient!(surface_flux_values, t, + calc_gradient_boundary_flux_by_direction!(surface_flux_values, t, boundary_conditions_parabolic[6], equations_parabolic, surface_integral, dg, cache, @@ -356,7 +372,7 @@ function calc_boundary_flux_gradients!(cache, t, return nothing end -function calc_boundary_flux_by_direction_gradient!(surface_flux_values::AbstractArray{<:Any, +function calc_gradient_boundary_flux_by_direction!(surface_flux_values::AbstractArray{<:Any, 5}, t, boundary_condition, @@ -889,6 +905,43 @@ end return nothing end +function calc_gradient_volume_integral!(gradients, u_transformed, + mesh::TreeMesh{3}, # for dispatch only + equations_parabolic::AbstractEquationsParabolic, + dg::DGSEM, cache) + @unpack derivative_dhat = dg.basis + gradients_x, gradients_y, gradients_z = gradients + + @threaded for element in eachelement(dg, cache) + # Calculate volume terms in one element, + # corresponds to `kernel` functions for the hyperbolic part of the flux + for k in eachnode(dg), j in eachnode(dg), i in eachnode(dg) + u_node = get_node_vars(u_transformed, equations_parabolic, dg, + i, j, k, element) + + for ii in eachnode(dg) + multiply_add_to_node_vars!(gradients_x, derivative_dhat[ii, i], + u_node, equations_parabolic, dg, + ii, j, k, element) + end + + for jj in eachnode(dg) + multiply_add_to_node_vars!(gradients_y, derivative_dhat[jj, j], + u_node, equations_parabolic, dg, + i, jj, k, element) + end + + for kk in eachnode(dg) + multiply_add_to_node_vars!(gradients_z, derivative_dhat[kk, k], + u_node, equations_parabolic, dg, + i, j, kk, element) + end + end + end + + return nothing +end + function calc_gradient_interface_flux!(surface_flux_values, mesh::TreeMesh{3}, equations, dg::DG, parabolic_scheme, @@ -927,6 +980,108 @@ function calc_gradient_interface_flux!(surface_flux_values, return nothing end +function calc_gradient_surface_integral!(gradients, + mesh::TreeMesh{3}, # for dispatch only + equations_parabolic::AbstractEquationsParabolic, + dg::DGSEM, cache, cache_parabolic) + @unpack boundary_interpolation = dg.basis + @unpack surface_flux_values = cache_parabolic.elements + + gradients_x, gradients_y, gradients_z = gradients + + # Note that all fluxes have been computed with outward-pointing normal vectors. + # Access the factors only once before beginning the loop to increase performance. + # We also use explicit assignments instead of `+=` to let `@muladd` turn these + # into FMAs (see comment at the top of the file). + factor_1 = boundary_interpolation[1, 1] + factor_2 = boundary_interpolation[nnodes(dg), 2] + @threaded for element in eachelement(dg, cache) + for m in eachnode(dg), l in eachnode(dg) + for v in eachvariable(equations_parabolic) + # surface at -x + gradients_x[v, 1, l, m, element] = (gradients_x[v, + 1, + l, + m, + element] - + surface_flux_values[v, + l, + m, + 1, + element] * + factor_1) + + # surface at +x + gradients_x[v, nnodes(dg), l, m, element] = (gradients_x[v, + nnodes(dg), + l, + m, + element] + + surface_flux_values[v, + l, + m, + 2, + element] * + factor_2) + + # surface at -y + gradients_y[v, l, 1, m, element] = (gradients_y[v, + l, + 1, + m, + element] - + surface_flux_values[v, + l, + m, + 3, + element] * + factor_1) + + # surface at +y + gradients_y[v, l, nnodes(dg), m, element] = (gradients_y[v, + l, + nnodes(dg), + m, + element] + + surface_flux_values[v, + l, + m, + 4, + element] * + factor_2) + + # surface at -z + gradients_z[v, l, m, 1, element] = (gradients_z[v, + l, + m, + 1, + element] - + surface_flux_values[v, + l, + m, + 5, + element] * + factor_1) + + # surface at +z + gradients_z[v, l, m, nnodes(dg), element] = (gradients_z[v, + l, + m, + nnodes(dg), + element] + + surface_flux_values[v, + l, + m, + 6, + element] * + factor_2) + end + end + end + + return nothing +end + # Calculate the gradient of the transformed variables function calc_gradient!(gradients, u_transformed, t, mesh::TreeMesh{3}, equations_parabolic, @@ -943,33 +1098,8 @@ function calc_gradient!(gradients, u_transformed, t, # Calculate volume integral @trixi_timeit timer() "volume integral" begin - @unpack derivative_dhat = dg.basis - @threaded for element in eachelement(dg, cache) - - # Calculate volume terms in one element - for k in eachnode(dg), j in eachnode(dg), i in eachnode(dg) - u_node = get_node_vars(u_transformed, equations_parabolic, dg, - i, j, k, element) - - for ii in eachnode(dg) - multiply_add_to_node_vars!(gradients_x, derivative_dhat[ii, i], - u_node, equations_parabolic, dg, - ii, j, k, element) - end - - for jj in eachnode(dg) - multiply_add_to_node_vars!(gradients_y, derivative_dhat[jj, j], - u_node, equations_parabolic, dg, - i, jj, k, element) - end - - for kk in eachnode(dg) - multiply_add_to_node_vars!(gradients_z, derivative_dhat[kk, k], - u_node, equations_parabolic, dg, - i, j, kk, element) - end - end - end + calc_gradient_volume_integral!(gradients, u_transformed, + mesh, equations_parabolic, dg, cache) end # Prolong solution to interfaces @@ -994,9 +1124,9 @@ function calc_gradient!(gradients, u_transformed, t, # Calculate boundary fluxes @trixi_timeit timer() "boundary flux" begin - calc_boundary_flux_gradients!(cache_parabolic, t, boundary_conditions_parabolic, - mesh, equations_parabolic, - dg.surface_integral, dg) + calc_gradient_boundary_flux!(cache_parabolic, t, boundary_conditions_parabolic, + mesh, equations_parabolic, + dg.surface_integral, dg) end # Prolong solution to mortars @@ -1017,76 +1147,8 @@ function calc_gradient!(gradients, u_transformed, t, # Calculate surface integrals @trixi_timeit timer() "surface integral" begin - @unpack boundary_interpolation = dg.basis - @unpack surface_flux_values = cache_parabolic.elements - - # Note that all fluxes have been computed with outward-pointing normal vectors. - # Access the factors only once before beginning the loop to increase performance. - # We also use explicit assignments instead of `+=` to let `@muladd` turn these - # into FMAs (see comment at the top of the file). - factor_1 = boundary_interpolation[1, 1] - factor_2 = boundary_interpolation[nnodes(dg), 2] - @threaded for element in eachelement(dg, cache) - for m in eachnode(dg), l in eachnode(dg) - for v in eachvariable(equations_parabolic) - # surface at -x - gradients_x[v, 1, l, m, element] = (gradients_x[v, 1, l, m, - element] - - surface_flux_values[v, l, m, 1, - element] * - factor_1) - - # surface at +x - gradients_x[v, nnodes(dg), l, m, element] = (gradients_x[v, - nnodes(dg), - l, m, - element] + - surface_flux_values[v, - l, - m, - 2, - element] * - factor_2) - - # surface at -y - gradients_y[v, l, 1, m, element] = (gradients_y[v, l, 1, m, - element] - - surface_flux_values[v, l, m, 3, - element] * - factor_1) - - # surface at +y - gradients_y[v, l, nnodes(dg), m, element] = (gradients_y[v, l, - nnodes(dg), - m, - element] + - surface_flux_values[v, - l, - m, - 4, - element] * - factor_2) - - # surface at -z - gradients_z[v, l, m, 1, element] = (gradients_z[v, l, m, 1, - element] - - surface_flux_values[v, l, m, 5, - element] * - factor_1) - - # surface at +z - gradients_z[v, l, m, nnodes(dg), element] = (gradients_z[v, l, m, - nnodes(dg), - element] + - surface_flux_values[v, - l, - m, - 6, - element] * - factor_2) - end - end - end + calc_gradient_surface_integral!(gradients, mesh, equations_parabolic, + dg, cache, cache_parabolic) end # Apply Jacobian from mapping to reference element diff --git a/test/test_parabolic_3d.jl b/test/test_parabolic_3d.jl index 1ee64ee7232..34a7cfd337e 100644 --- a/test/test_parabolic_3d.jl +++ b/test/test_parabolic_3d.jl @@ -501,17 +501,17 @@ end "elixir_navierstokes_crm.jl"), l2=[ 2.2353998537135728e-10, - 6.438639847871375e-8, - 4.8159768090963756e-8, - 6.287534321002294e-8, - 7.596958320246443e-5 + 6.438248564149087e-8, + 4.8158393932278836e-8, + 6.287149707414103e-8, + 7.591822518308327e-5 ], linf=[ 0.36753460935979454, 209.1460932734352, - 51.3036359192352, - 41.461135962058904, - 113004.3158457772 + 51.303610073892436, + 41.45867326162659, + 113004.34493505303 ], tspan=(0.0, 1e-10)) # Ensure that we do not have excessive memory allocations diff --git a/test/test_structured_2d.jl b/test/test_structured_2d.jl index 578401bba12..12dbc1be64d 100644 --- a/test/test_structured_2d.jl +++ b/test/test_structured_2d.jl @@ -255,7 +255,7 @@ end ]) # Ensure that we do not have excessive memory allocations # (e.g., from type instabilities) - @test_allocations(Trixi.rhs!, semi_float_type, sol, 1000) + @test_allocations(Trixi.rhs!, semi, sol, 1000) end @trixi_testset "elixir_eulermulti_convergence_ec.jl" begin diff --git a/test/test_threaded.jl b/test/test_threaded.jl index ad7c49c5923..85aa726da02 100644 --- a/test/test_threaded.jl +++ b/test/test_threaded.jl @@ -91,7 +91,11 @@ Trixi.MPI.Barrier(Trixi.mpi_comm()) # Ensure that we do not have excessive memory allocations # (e.g., from type instabilities) - @test_allocations(Trixi.rhs!, semi, sol, 5000) + if VERSION >= v"1.12" + @test_allocations(Trixi.rhs!, semi, sol, 7500) + else + @test_allocations(Trixi.rhs!, semi, sol, 5000) + end end @trixi_testset "elixir_euler_ec.jl" begin @@ -133,7 +137,11 @@ Trixi.MPI.Barrier(Trixi.mpi_comm()) tspan=(0.0, 1.0),) # Ensure that we do not have excessive memory allocations # (e.g., from type instabilities) - @test_allocations(Trixi.rhs!, semi, sol, 5000) + if VERSION >= v"1.12" + @test_allocations(Trixi.rhs!, semi, sol, 7500) + else + @test_allocations(Trixi.rhs!, semi, sol, 5000) + end end @trixi_testset "elixir_advection_diffusion.jl" begin @@ -180,7 +188,11 @@ Trixi.MPI.Barrier(Trixi.mpi_comm()) # Ensure that we do not have excessive memory allocations # (e.g., from type instabilities) - @test_allocations(Trixi.rhs!, semi, sol, 5000) + if VERSION >= v"1.12" + @test_allocations(Trixi.rhs!, semi, sol, 15000) + else + @test_allocations(Trixi.rhs!, semi, sol, 5000) + end end end @@ -259,7 +271,11 @@ end # Ensure that we do not have excessive memory allocations # (e.g., from type instabilities) - @test_allocations(Trixi.rhs!, semi, sol, 5000) + if VERSION >= v"1.12" + @test_allocations(Trixi.rhs!, semi, sol, 7500) + else + @test_allocations(Trixi.rhs!, semi, sol, 5000) + end end @trixi_testset "elixir_eulergravity_convergence.jl" begin From d97c162d3289e564fbcdbae5fc7b439b03268d66 Mon Sep 17 00:00:00 2001 From: Vincent Marks Date: Mon, 27 Oct 2025 17:45:39 +0100 Subject: [PATCH 34/45] Apply suggestions from code review Co-authored-by: Hendrik Ranocha --- src/solvers/default_boundary_conditions.jl | 23 ++++++++-------------- 1 file changed, 8 insertions(+), 15 deletions(-) diff --git a/src/solvers/default_boundary_conditions.jl b/src/solvers/default_boundary_conditions.jl index a4796b15208..d2bc5c7560a 100644 --- a/src/solvers/default_boundary_conditions.jl +++ b/src/solvers/default_boundary_conditions.jl @@ -1,7 +1,7 @@ """ boundary_condition_default(mesh::P4estMesh{2,2}, boundary_condition) -Create a default boundary condition dictionary for p4est meshes in 2D +Create a default boundary condition dictionary for [`P4estMesh`](@ref)es in 2D that use the standard boundary naming convention. This function applies the same boundary condition to all standard boundaries: @@ -17,7 +17,6 @@ This function applies the same boundary condition to all standard boundaries: - `Dict{Symbol, Any}`: Dictionary mapping boundary names to the boundary condition """ function boundary_condition_default(mesh::P4estMesh{2,2}, boundary_condition) - return Dict(:x_neg => boundary_condition, :y_neg => boundary_condition, :y_pos => boundary_condition, @@ -27,7 +26,7 @@ end """ boundary_condition_default(mesh::P4estMesh{3,3}, boundary_condition) -Create a default boundary condition dictionary for p4est meshes in 3D +Create a default boundary condition dictionary for [`P4estMesh`](@ref)es in 3D that use the standard boundary naming convention. This function applies the same boundary condition to all standard boundaries: - `:x_neg`: negative x-direction boundary @@ -42,7 +41,6 @@ This function applies the same boundary condition to all standard boundaries: - `Dict{Symbol, Any}`: Dictionary mapping boundary names to the boundary condition """ function boundary_condition_default(mesh::P4estMesh{3,3}, boundary_condition) - return Dict(:x_neg => boundary_condition, :x_pos => boundary_condition, :y_neg => boundary_condition, @@ -54,7 +52,7 @@ end """ boundary_condition_default(mesh::StructuredMesh1D, boundary_condition) -Create a default boundary condition dictionary for structured meshes in 1D +Create a default boundary condition dictionary for [`StructuredMesh`](@ref)es in 1D that use the standard boundary naming convention. This function applies the same boundary condition to all standard boundaries: - `:x_neg`: negative x-direction boundary @@ -65,7 +63,6 @@ This function applies the same boundary condition to all standard boundaries: - Named tuple mapping boundary names to the boundary condition """ function boundary_condition_default(mesh::StructuredMesh{1}, boundary_condition) - return (x_neg = boundary_condition, x_pos = boundary_condition) end @@ -73,7 +70,7 @@ end """ boundary_condition_default(mesh::StructuredMesh2D, boundary_condition) -Create a default boundary condition dictionary for structured meshes in 2D +Create a default boundary condition dictionary for [`StructuredMesh`](@ref)es in 2D that use the standard boundary naming convention. This function applies the same boundary condition to all standard boundaries: - `:x_neg`: negative x-direction boundary @@ -86,7 +83,6 @@ This function applies the same boundary condition to all standard boundaries: - Named tuple mapping boundary names to the boundary condition """ function boundary_condition_default(mesh::StructuredMesh{2}, boundary_condition) - return (x_neg = boundary_condition, x_pos = boundary_condition, y_neg = boundary_condition, @@ -96,7 +92,7 @@ end """ boundary_condition_default(mesh::StructuredMesh3D, boundary_condition) -Create a default boundary condition dictionary for structured meshes in 3D +Create a default boundary condition dictionary for [`StructuredMesh`](@ref)es in 3D that use the standard boundary naming convention. This function applies the same boundary condition to all standard boundaries: - `:x_neg`: negative x-direction boundary @@ -111,7 +107,6 @@ This function applies the same boundary condition to all standard boundaries: - Named tuple mapping boundary names to the boundary condition """ function boundary_condition_default(mesh::StructuredMesh{3}, boundary_condition) - return (x_neg = boundary_condition, x_pos = boundary_condition, y_neg = boundary_condition, @@ -124,7 +119,7 @@ end """ boundary_condition_default(mesh::TreeMesh1D, boundary_condition) -Create a default boundary condition dictionary for tree meshes in 1D +Create a default boundary condition dictionary for [`TreeMesh`](@ref)es in 1D that use the standard boundary naming convention. This function applies the same boundary condition to all standard boundaries: - `:x_neg`: negative x-direction boundary @@ -142,7 +137,7 @@ end """ boundary_condition_default(mesh::TreeMesh2D, boundary_condition) -Create a default boundary condition dictionary for tree meshes in 2D +Create a default boundary condition dictionary for [`TreeMesh`](@ref)es in 2D that use the standard boundary naming convention. This function applies the same boundary condition to all standard boundaries: - `:x_neg`: negative x-direction boundary @@ -155,7 +150,6 @@ This function applies the same boundary condition to all standard boundaries: - Named tuple mapping boundary names to the boundary condition """ function boundary_condition_default(mesh::TreeMesh2D, boundary_condition) - return (x_neg = boundary_condition, x_pos = boundary_condition, y_neg = boundary_condition, @@ -165,7 +159,7 @@ end """ boundary_condition_default(mesh::TreeMesh3D, boundary_condition) -Create a default boundary condition dictionary for tree meshes in 3D +Create a default boundary condition dictionary for [`TreeMesh`](@ref)es in 3D that use the standard boundary naming convention. This function applies the same boundary condition to all standard boundaries: - `:x_neg`: negative x-direction boundary @@ -180,7 +174,6 @@ This function applies the same boundary condition to all standard boundaries: - Named tuple mapping boundary names to the boundary condition """ function boundary_condition_default(mesh::TreeMesh3D, boundary_condition) - return (x_neg = boundary_condition, x_pos = boundary_condition, y_neg = boundary_condition, From 0583638f4c1790bea1e7b1437904bd094c4f0195 Mon Sep 17 00:00:00 2001 From: vincmarks Date: Mon, 27 Oct 2025 18:42:04 +0100 Subject: [PATCH 35/45] using the function ``boundary_condition_default `` --- ...xir_advection_diffusion_nonperiodic_amr.jl | 16 ++++++------- .../elixir_euler_source_terms_nonperiodic.jl | 20 ++++++++-------- .../elixir_euler_source_terms_nonperiodic.jl | 18 +++++++------- ...lixir_euler_rayleigh_taylor_instability.jl | 17 ++++++------- ...r_euler_source_terms_nonperiodic_curved.jl | 24 +++++++++---------- .../elixir_euler_source_terms_nonperiodic.jl | 18 +++++++------- .../elixir_euler_source_terms_nonperiodic.jl | 20 ++++++++-------- .../elixir_hypdiff_nonperiodic.jl | 14 +++++------ 8 files changed, 74 insertions(+), 73 deletions(-) diff --git a/examples/p4est_2d_dgsem/elixir_advection_diffusion_nonperiodic_amr.jl b/examples/p4est_2d_dgsem/elixir_advection_diffusion_nonperiodic_amr.jl index f9de74d9af1..35effcf9729 100644 --- a/examples/p4est_2d_dgsem/elixir_advection_diffusion_nonperiodic_amr.jl +++ b/examples/p4est_2d_dgsem/elixir_advection_diffusion_nonperiodic_amr.jl @@ -45,14 +45,14 @@ boundary_conditions = Dict(:x_neg => BoundaryConditionDirichlet(initial_conditio :y_pos => BoundaryConditionDirichlet(initial_condition), :x_pos => boundary_condition_do_nothing) -boundary_conditions_parabolic = Dict(:x_neg => BoundaryConditionDirichlet(initial_condition), - :x_pos => BoundaryConditionDirichlet(initial_condition), - :y_neg => BoundaryConditionDirichlet(initial_condition), - :y_pos => BoundaryConditionDirichlet(initial_condition)) - -# Alternative to the boundary_conditions_parabolic defined above: -# boundary_condition = BoundaryConditionDirichlet(initial_condition) -# boundary_conditions_parabolic = boundary_condition_default(mesh, boundary_condition) +# Assign a single boundary condition to all boundaries +boundary_condition = BoundaryConditionDirichlet(initial_condition) +boundary_conditions_parabolic = boundary_condition_default(mesh, boundary_condition) +# Alternatively, you can use +# boundary_conditions_parabolic = Dict(:x_neg => BoundaryConditionDirichlet(initial_condition), +# :x_pos => BoundaryConditionDirichlet(initial_condition), +# :y_neg => BoundaryConditionDirichlet(initial_condition), +# :y_pos => BoundaryConditionDirichlet(initial_condition)) # A semidiscretization collects data structures and functions for the spatial discretization semi = SemidiscretizationHyperbolicParabolic(mesh, diff --git a/examples/p4est_3d_dgsem/elixir_euler_source_terms_nonperiodic.jl b/examples/p4est_3d_dgsem/elixir_euler_source_terms_nonperiodic.jl index c81d179f1bb..2903df96c97 100644 --- a/examples/p4est_3d_dgsem/elixir_euler_source_terms_nonperiodic.jl +++ b/examples/p4est_3d_dgsem/elixir_euler_source_terms_nonperiodic.jl @@ -8,14 +8,6 @@ equations = CompressibleEulerEquations3D(1.4) initial_condition = initial_condition_convergence_test -boundary_condition = BoundaryConditionDirichlet(initial_condition) -boundary_conditions = Dict(:x_neg => boundary_condition, - :x_pos => boundary_condition, - :y_neg => boundary_condition, - :y_pos => boundary_condition, - :z_neg => boundary_condition, - :z_pos => boundary_condition) - # Up to version 0.13.0, `max_abs_speed_naive` was used as the default wave speed estimate of # `const flux_lax_friedrichs = FluxLaxFriedrichs(), i.e., `FluxLaxFriedrichs(max_abs_speed = max_abs_speed_naive)`. # In the `StepsizeCallback`, though, the less diffusive `max_abs_speeds` is employed which is consistent with `max_abs_speed`. @@ -35,8 +27,16 @@ mesh = P4estMesh(trees_per_dimension, polydeg = 1, coordinates_min = coordinates_min, coordinates_max = coordinates_max, periodicity = false, initial_refinement_level = 1) -# Alternative to the boundary_conditions defined above: -# boundary_conditions = boundary_condition_default(mesh, boundary_condition) +# Assign a single boundary condition to all boundaries +boundary_condition = BoundaryConditionDirichlet(initial_condition) +boundary_conditions = boundary_condition_default(mesh, boundary_condition) +# Alternatively, you can use +# boundary_conditions = Dict(:x_neg => boundary_condition, +# :x_pos => boundary_condition, +# :y_neg => boundary_condition, +# :y_pos => boundary_condition, +# :z_neg => boundary_condition, +# :z_pos => boundary_condition) semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver, source_terms = source_terms_convergence_test, diff --git a/examples/structured_1d_dgsem/elixir_euler_source_terms_nonperiodic.jl b/examples/structured_1d_dgsem/elixir_euler_source_terms_nonperiodic.jl index 1c83774dbfb..0d625b31086 100644 --- a/examples/structured_1d_dgsem/elixir_euler_source_terms_nonperiodic.jl +++ b/examples/structured_1d_dgsem/elixir_euler_source_terms_nonperiodic.jl @@ -10,13 +10,6 @@ initial_condition = initial_condition_convergence_test source_terms = source_terms_convergence_test -# you can either use a single function to impose the BCs weakly in all -# 1*ndims == 2 directions or you can pass a tuple containing BCs for -# each direction -boundary_condition = BoundaryConditionDirichlet(initial_condition) -boundary_conditions = (x_neg = boundary_condition, - x_pos = boundary_condition) - # Up to version 0.13.0, `max_abs_speed_naive` was used as the default wave speed estimate of # `const flux_lax_friedrichs = FluxLaxFriedrichs(), i.e., `FluxLaxFriedrichs(max_abs_speed = max_abs_speed_naive)`. # In the `StepsizeCallback`, though, the less diffusive `max_abs_speeds` is employed which is consistent with `max_abs_speed`. @@ -30,8 +23,15 @@ f1() = SVector(0.0) f2() = SVector(2.0) mesh = StructuredMesh((16,), (f1, f2), periodicity = false) -# Alternative to the boundary_conditions defined above: -# boundary_conditions = boundary_condition_default(mesh, boundary_condition) +# you can either use a single function to impose the BCs weakly in all +# 1*ndims == 2 directions or you can pass a tuple containing BCs for +# each direction +# Assign a single boundary condition to all boundaries +boundary_condition = BoundaryConditionDirichlet(initial_condition) +boundary_conditions = boundary_condition_default(mesh, boundary_condition) +# Alternatively, you can use +# boundary_conditions = (x_neg = boundary_condition, +# x_pos = boundary_condition) semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver, source_terms = source_terms, diff --git a/examples/structured_2d_dgsem/elixir_euler_rayleigh_taylor_instability.jl b/examples/structured_2d_dgsem/elixir_euler_rayleigh_taylor_instability.jl index 805735b6ad4..6579b9bb1ea 100644 --- a/examples/structured_2d_dgsem/elixir_euler_rayleigh_taylor_instability.jl +++ b/examples/structured_2d_dgsem/elixir_euler_rayleigh_taylor_instability.jl @@ -71,14 +71,15 @@ cells_per_dimension = (num_elements_per_dimension, num_elements_per_dimension * mesh = StructuredMesh(cells_per_dimension, mapping, periodicity = false) initial_condition = initial_condition_rayleigh_taylor_instability -boundary_conditions = (x_neg = boundary_condition_slip_wall, - x_pos = boundary_condition_slip_wall, - y_neg = boundary_condition_slip_wall, - y_pos = boundary_condition_slip_wall) - -# Alternative to the boundary_conditions defined above: -# boundary_condition = boundary_condition_slip_wall -# boundary_conditions = boundary_condition_default(mesh, boundary_condition) + +# Assign a single boundary condition to all boundaries +boundary_condition = boundary_condition_slip_wall +boundary_conditions = boundary_condition_default(mesh, boundary_condition) +# Alternatively, you can use +# boundary_conditions = (x_neg = boundary_condition_slip_wall, +# x_pos = boundary_condition_slip_wall, +# y_neg = boundary_condition_slip_wall, +# y_pos = boundary_condition_slip_wall) # # Alternative setup: left/right periodic BCs and Dirichlet BCs on the top/bottom. # boundary_conditions = ( diff --git a/examples/structured_3d_dgsem/elixir_euler_source_terms_nonperiodic_curved.jl b/examples/structured_3d_dgsem/elixir_euler_source_terms_nonperiodic_curved.jl index 1ffc2bc2b80..b55df45c77e 100644 --- a/examples/structured_3d_dgsem/elixir_euler_source_terms_nonperiodic_curved.jl +++ b/examples/structured_3d_dgsem/elixir_euler_source_terms_nonperiodic_curved.jl @@ -8,16 +8,6 @@ equations = CompressibleEulerEquations3D(1.4) initial_condition = initial_condition_convergence_test -# you can either use a single function to impose the BCs weakly in all -# 2*ndims == 4 directions or you can pass a tuple containing BCs for each direction -boundary_condition = BoundaryConditionDirichlet(initial_condition) -boundary_conditions = (x_neg = boundary_condition, - x_pos = boundary_condition, - y_neg = boundary_condition, - y_pos = boundary_condition, - z_neg = boundary_condition, - z_pos = boundary_condition) - # Up to version 0.13.0, `max_abs_speed_naive` was used as the default wave speed estimate of # `const flux_lax_friedrichs = FluxLaxFriedrichs(), i.e., `FluxLaxFriedrichs(max_abs_speed = max_abs_speed_naive)`. # In the `StepsizeCallback`, though, the less diffusive `max_abs_speeds` is employed which is consistent with `max_abs_speed`. @@ -57,8 +47,18 @@ end cells_per_dimension = (4, 4, 4) mesh = StructuredMesh(cells_per_dimension, mapping, periodicity = false) -# Alternative to the boundary_conditions defined above: -# boundary_conditions = boundary_condition_default(mesh, boundary_condition) +# you can either use a single function to impose the BCs weakly in all +# 2*ndims == 4 directions or you can pass a tuple containing BCs for each direction +# Assign a single boundary condition to all boundaries +boundary_condition = BoundaryConditionDirichlet(initial_condition) +boundary_conditions = boundary_condition_default(mesh, boundary_condition) +# Alternatively, you can use +# boundary_conditions = (x_neg = boundary_condition, +# x_pos = boundary_condition, +# y_neg = boundary_condition, +# y_pos = boundary_condition, +# z_neg = boundary_condition, +# z_pos = boundary_condition) semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver, source_terms = source_terms_convergence_test, diff --git a/examples/tree_1d_dgsem/elixir_euler_source_terms_nonperiodic.jl b/examples/tree_1d_dgsem/elixir_euler_source_terms_nonperiodic.jl index 715768d01a9..614aad34e41 100644 --- a/examples/tree_1d_dgsem/elixir_euler_source_terms_nonperiodic.jl +++ b/examples/tree_1d_dgsem/elixir_euler_source_terms_nonperiodic.jl @@ -8,13 +8,6 @@ equations = CompressibleEulerEquations1D(1.4) initial_condition = initial_condition_convergence_test -# you can either use a single function to impose the BCs weakly in all -# 1*ndims == 2 directions or you can pass a tuple containing BCs for -# each direction -boundary_condition = BoundaryConditionDirichlet(initial_condition) -boundary_conditions = (x_neg = boundary_condition, - x_pos = boundary_condition) - # Up to version 0.13.0, `max_abs_speed_naive` was used as the default wave speed estimate of # `const flux_lax_friedrichs = FluxLaxFriedrichs(), i.e., `FluxLaxFriedrichs(max_abs_speed = max_abs_speed_naive)`. # In the `StepsizeCallback`, though, the less diffusive `max_abs_speeds` is employed which is consistent with `max_abs_speed`. @@ -31,8 +24,15 @@ mesh = TreeMesh(coordinates_min, coordinates_max, n_cells_max = 10_000, periodicity = false) -# Alternative to the boundary_conditions defined above: -# boundary_conditions = boundary_condition_default(mesh, boundary_condition) +# you can either use a single function to impose the BCs weakly in all +# 1*ndims == 2 directions or you can pass a tuple containing BCs for +# each direction +# Assign a single boundary condition to all boundaries +boundary_condition = BoundaryConditionDirichlet(initial_condition) +boundary_conditions = boundary_condition_default(mesh, boundary_condition) +# Alternatively, you can use +# boundary_conditions = (x_neg = boundary_condition, +# x_pos = boundary_condition) semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver, source_terms = source_terms_convergence_test, diff --git a/examples/tree_2d_dgsem/elixir_euler_source_terms_nonperiodic.jl b/examples/tree_2d_dgsem/elixir_euler_source_terms_nonperiodic.jl index e9780a3700a..343b6789a2b 100644 --- a/examples/tree_2d_dgsem/elixir_euler_source_terms_nonperiodic.jl +++ b/examples/tree_2d_dgsem/elixir_euler_source_terms_nonperiodic.jl @@ -8,14 +8,6 @@ equations = CompressibleEulerEquations2D(1.4) initial_condition = initial_condition_convergence_test -# you can either use a single function to impose the BCs weakly in all -# 2*ndims == 4 directions or you can pass a tuple containing BCs for each direction -boundary_condition = BoundaryConditionDirichlet(initial_condition) -boundary_conditions = (x_neg = boundary_condition, - x_pos = boundary_condition, - y_neg = boundary_condition, - y_pos = boundary_condition) - # Up to version 0.13.0, `max_abs_speed_naive` was used as the default wave speed estimate of # `const flux_lax_friedrichs = FluxLaxFriedrichs(), i.e., `FluxLaxFriedrichs(max_abs_speed = max_abs_speed_naive)`. # In the `StepsizeCallback`, though, the less diffusive `max_abs_speeds` is employed which is consistent with `max_abs_speed`. @@ -32,8 +24,16 @@ mesh = TreeMesh(coordinates_min, coordinates_max, n_cells_max = 10_000, periodicity = false) -# Alternative to the boundary_conditions defined above: -# boundary_conditions = boundary_condition_default(mesh, boundary_condition) +# you can either use a single function to impose the BCs weakly in all +# 2*ndims == 4 directions or you can pass a tuple containing BCs for each direction +# Assign a single boundary condition to all boundaries +boundary_condition = BoundaryConditionDirichlet(initial_condition) +boundary_conditions = boundary_condition_default(mesh, boundary_condition) +# Alternatively, you can use +# boundary_conditions = (x_neg = boundary_condition, +# x_pos = boundary_condition, +# y_neg = boundary_condition, +# y_pos = boundary_condition) semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver, source_terms = source_terms_convergence_test, diff --git a/examples/tree_3d_dgsem/elixir_hypdiff_nonperiodic.jl b/examples/tree_3d_dgsem/elixir_hypdiff_nonperiodic.jl index 0dc23d736b9..9c2a30122b3 100644 --- a/examples/tree_3d_dgsem/elixir_hypdiff_nonperiodic.jl +++ b/examples/tree_3d_dgsem/elixir_hypdiff_nonperiodic.jl @@ -6,12 +6,6 @@ using Trixi equations = HyperbolicDiffusionEquations3D() initial_condition = initial_condition_poisson_nonperiodic -boundary_conditions = (x_neg = boundary_condition_poisson_nonperiodic, - x_pos = boundary_condition_poisson_nonperiodic, - y_neg = boundary_condition_periodic, - y_pos = boundary_condition_periodic, - z_neg = boundary_condition_periodic, - z_pos = boundary_condition_periodic) solver = DGSEM(polydeg = 4, surface_flux = flux_lax_friedrichs) @@ -22,7 +16,13 @@ mesh = TreeMesh(coordinates_min, coordinates_max, n_cells_max = 30_000, periodicity = (false, true, true)) -# Alternative to the boundary_conditions defined above if you want to apply the same boundary condition to all faces of the mesh: +boundary_conditions = (x_neg = boundary_condition_poisson_nonperiodic, + x_pos = boundary_condition_poisson_nonperiodic, + y_neg = boundary_condition_periodic, + y_pos = boundary_condition_periodic, + z_neg = boundary_condition_periodic, + z_pos = boundary_condition_periodic) +# Alternative assign a single boundary condition to all boundaries # boundary_condition = boundary_condition_periodic # boundary_conditions = boundary_condition_default(mesh, boundary_condition) From a6a9baea6ccaa7dd23a5d498d7bd0cfe014936f9 Mon Sep 17 00:00:00 2001 From: vincmarks Date: Mon, 27 Oct 2025 18:54:54 +0100 Subject: [PATCH 36/45] fixed typo --- examples/tree_3d_dgsem/elixir_hypdiff_nonperiodic.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/examples/tree_3d_dgsem/elixir_hypdiff_nonperiodic.jl b/examples/tree_3d_dgsem/elixir_hypdiff_nonperiodic.jl index 9c2a30122b3..b0c045a067b 100644 --- a/examples/tree_3d_dgsem/elixir_hypdiff_nonperiodic.jl +++ b/examples/tree_3d_dgsem/elixir_hypdiff_nonperiodic.jl @@ -22,7 +22,7 @@ boundary_conditions = (x_neg = boundary_condition_poisson_nonperiodic, y_pos = boundary_condition_periodic, z_neg = boundary_condition_periodic, z_pos = boundary_condition_periodic) -# Alternative assign a single boundary condition to all boundaries +# Alternatively assign a single boundary condition to all boundaries # boundary_condition = boundary_condition_periodic # boundary_conditions = boundary_condition_default(mesh, boundary_condition) From bfd6fdf592283e5fd451e9389a5fa63ebc12d983 Mon Sep 17 00:00:00 2001 From: vincmarks Date: Tue, 28 Oct 2025 14:31:44 +0100 Subject: [PATCH 37/45] using the function boundary_condition_default --- .../elixir_advection_extended.jl | 20 +++++++++++-------- 1 file changed, 12 insertions(+), 8 deletions(-) diff --git a/examples/tree_3d_dgsem/elixir_advection_extended.jl b/examples/tree_3d_dgsem/elixir_advection_extended.jl index 7c719d653d6..63ac204da84 100644 --- a/examples/tree_3d_dgsem/elixir_advection_extended.jl +++ b/examples/tree_3d_dgsem/elixir_advection_extended.jl @@ -9,13 +9,6 @@ equations = LinearScalarAdvectionEquation3D(advection_velocity) initial_condition = initial_condition_convergence_test -# you can either use a single function to impose the BCs weakly in all -# 1*ndims == 2 directions or you can pass a tuple containing BCs for -# each direction -# Note: "boundary_condition_periodic" indicates that it is a periodic boundary and can be omitted on -# fully periodic domains. Here, however, it is included to allow easy override during testing -boundary_conditions = boundary_condition_periodic - # 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) @@ -28,6 +21,17 @@ mesh = TreeMesh(coordinates_min, coordinates_max, n_cells_max = 30_000, # set maximum capacity of tree data structure periodicity = true) +# you can either use a single function to impose the BCs weakly in all +# 1*ndims == 2 directions or you can pass a tuple containing BCs for +# each direction +# Note: "boundary_condition_periodic" indicates that it is a periodic boundary and can be omitted on +# fully periodic domains. Here, however, it is included to allow easy override during testing +# Assign a single boundary condition to all boundaries +boundary_condition = boundary_condition_periodic +boundary_conditions = boundary_condition_default(mesh, boundary_condition) +# Alternatively, you can use +# boundary_conditions = boundary_condition_periodic + # A semidiscretization collects data structures and functions for the spatial discretization semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver, boundary_conditions = boundary_conditions) @@ -76,4 +80,4 @@ callbacks = CallbackSet(summary_callback, # 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); + ode_default_options()..., callback = callbacks) From 0d3588851a36461ef1d46a9fbc07df730243a1d4 Mon Sep 17 00:00:00 2001 From: vincmarks Date: Tue, 28 Oct 2025 14:32:53 +0100 Subject: [PATCH 38/45] not using the function boundary_condition_default here --- examples/tree_3d_dgsem/elixir_hypdiff_nonperiodic.jl | 3 --- 1 file changed, 3 deletions(-) diff --git a/examples/tree_3d_dgsem/elixir_hypdiff_nonperiodic.jl b/examples/tree_3d_dgsem/elixir_hypdiff_nonperiodic.jl index b0c045a067b..51ada9cc233 100644 --- a/examples/tree_3d_dgsem/elixir_hypdiff_nonperiodic.jl +++ b/examples/tree_3d_dgsem/elixir_hypdiff_nonperiodic.jl @@ -22,9 +22,6 @@ boundary_conditions = (x_neg = boundary_condition_poisson_nonperiodic, y_pos = boundary_condition_periodic, z_neg = boundary_condition_periodic, z_pos = boundary_condition_periodic) -# Alternatively assign a single boundary condition to all boundaries -# boundary_condition = boundary_condition_periodic -# boundary_conditions = boundary_condition_default(mesh, boundary_condition) semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver, source_terms = source_terms_poisson_nonperiodic, From f5538da7c7e5c6ec7d0eaed9c01eafa966f5c951 Mon Sep 17 00:00:00 2001 From: vincmarks Date: Tue, 28 Oct 2025 14:37:14 +0100 Subject: [PATCH 39/45] changed the file where the function `boundary_condition_default(mesh, boundary_condition)` is used in the examples/tree_3d_dgsem folder --- docs/src/meshes/tree_mesh.md | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/docs/src/meshes/tree_mesh.md b/docs/src/meshes/tree_mesh.md index f0d7d71e8b0..e1df92a64d0 100644 --- a/docs/src/meshes/tree_mesh.md +++ b/docs/src/meshes/tree_mesh.md @@ -13,4 +13,4 @@ dispatching on the `orientation::Integer` as described in the ### Boundary conditions For [`TreeMesh`](@ref)es, boundary conditions are defined and stored in named tuples (see, for example, `examples/tree_1d_dgsem/elixir_euler_source_terms_nonperiodic.jl`). -If you want to apply the same boundary condition to all faces of the mesh, you can use the `boundary_condition_default(mesh, boundary_condition)` function, as demonstrated in `examples/tree_1d_dgsem/elixir_euler_source_terms_nonperiodic.jl`, `examples/tree_2d_dgsem/elixir_euler_source_terms_nonperiodic.jl` and `examples/tree_3d_dgsem/elixir_hypdiff_nonperiodic.jl`. \ No newline at end of file +If you want to apply the same boundary condition to all faces of the mesh, you can use the `boundary_condition_default(mesh, boundary_condition)` function, as demonstrated in `examples/tree_1d_dgsem/elixir_euler_source_terms_nonperiodic.jl`, `examples/tree_2d_dgsem/elixir_euler_source_terms_nonperiodic.jl` and `examples/tree_3d_dgsem/elixir_advection_extended.jl`. \ No newline at end of file From 64152a158b82e3264c4d29b05c2a17f2bcf90a97 Mon Sep 17 00:00:00 2001 From: Vincent Marks Date: Wed, 29 Oct 2025 10:17:35 +0100 Subject: [PATCH 40/45] Apply suggestions from code review Co-authored-by: Hendrik Ranocha --- .../elixir_euler_rayleigh_taylor_instability.jl | 3 +-- .../tree_3d_dgsem/elixir_advection_extended.jl | 2 +- src/solvers/default_boundary_conditions.jl | 17 ++++++++--------- 3 files changed, 10 insertions(+), 12 deletions(-) diff --git a/examples/structured_2d_dgsem/elixir_euler_rayleigh_taylor_instability.jl b/examples/structured_2d_dgsem/elixir_euler_rayleigh_taylor_instability.jl index 6579b9bb1ea..1926a0f4c27 100644 --- a/examples/structured_2d_dgsem/elixir_euler_rayleigh_taylor_instability.jl +++ b/examples/structured_2d_dgsem/elixir_euler_rayleigh_taylor_instability.jl @@ -73,8 +73,7 @@ mesh = StructuredMesh(cells_per_dimension, mapping, periodicity = false) initial_condition = initial_condition_rayleigh_taylor_instability # Assign a single boundary condition to all boundaries -boundary_condition = boundary_condition_slip_wall -boundary_conditions = boundary_condition_default(mesh, boundary_condition) +boundary_conditions = boundary_condition_default(mesh, boundary_condition_slip_wall) # Alternatively, you can use # boundary_conditions = (x_neg = boundary_condition_slip_wall, # x_pos = boundary_condition_slip_wall, diff --git a/examples/tree_3d_dgsem/elixir_advection_extended.jl b/examples/tree_3d_dgsem/elixir_advection_extended.jl index 63ac204da84..d6655ca7a7f 100644 --- a/examples/tree_3d_dgsem/elixir_advection_extended.jl +++ b/examples/tree_3d_dgsem/elixir_advection_extended.jl @@ -80,4 +80,4 @@ callbacks = CallbackSet(summary_callback, # 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) + ode_default_options()..., callback = callbacks); diff --git a/src/solvers/default_boundary_conditions.jl b/src/solvers/default_boundary_conditions.jl index d2bc5c7560a..f7caf119fda 100644 --- a/src/solvers/default_boundary_conditions.jl +++ b/src/solvers/default_boundary_conditions.jl @@ -2,8 +2,7 @@ boundary_condition_default(mesh::P4estMesh{2,2}, boundary_condition) Create a default boundary condition dictionary for [`P4estMesh`](@ref)es in 2D -that use the standard boundary naming convention. - +that uses the standard boundary naming convention. This function applies the same boundary condition to all standard boundaries: - `:x_neg`: negative x-direction boundary - `:x_pos`: positive x-direction boundary @@ -27,7 +26,7 @@ end boundary_condition_default(mesh::P4estMesh{3,3}, boundary_condition) Create a default boundary condition dictionary for [`P4estMesh`](@ref)es in 3D -that use the standard boundary naming convention. +that uses the standard boundary naming convention. This function applies the same boundary condition to all standard boundaries: - `:x_neg`: negative x-direction boundary - `:x_pos`: positive x-direction boundary @@ -53,7 +52,7 @@ end boundary_condition_default(mesh::StructuredMesh1D, boundary_condition) Create a default boundary condition dictionary for [`StructuredMesh`](@ref)es in 1D -that use the standard boundary naming convention. +that uses the standard boundary naming convention. This function applies the same boundary condition to all standard boundaries: - `:x_neg`: negative x-direction boundary - `:x_pos`: positive x-direction boundary @@ -71,7 +70,7 @@ end boundary_condition_default(mesh::StructuredMesh2D, boundary_condition) Create a default boundary condition dictionary for [`StructuredMesh`](@ref)es in 2D -that use the standard boundary naming convention. +that uses the standard boundary naming convention. This function applies the same boundary condition to all standard boundaries: - `:x_neg`: negative x-direction boundary - `:x_pos`: positive x-direction boundary @@ -93,7 +92,7 @@ end boundary_condition_default(mesh::StructuredMesh3D, boundary_condition) Create a default boundary condition dictionary for [`StructuredMesh`](@ref)es in 3D -that use the standard boundary naming convention. +that uses the standard boundary naming convention. This function applies the same boundary condition to all standard boundaries: - `:x_neg`: negative x-direction boundary - `:x_pos`: positive x-direction boundary @@ -120,7 +119,7 @@ end boundary_condition_default(mesh::TreeMesh1D, boundary_condition) Create a default boundary condition dictionary for [`TreeMesh`](@ref)es in 1D -that use the standard boundary naming convention. +that uses the standard boundary naming convention. This function applies the same boundary condition to all standard boundaries: - `:x_neg`: negative x-direction boundary - `:x_pos`: positive x-direction boundary @@ -138,7 +137,7 @@ end boundary_condition_default(mesh::TreeMesh2D, boundary_condition) Create a default boundary condition dictionary for [`TreeMesh`](@ref)es in 2D -that use the standard boundary naming convention. +that uses the standard boundary naming convention. This function applies the same boundary condition to all standard boundaries: - `:x_neg`: negative x-direction boundary - `:x_pos`: positive x-direction boundary @@ -160,7 +159,7 @@ end boundary_condition_default(mesh::TreeMesh3D, boundary_condition) Create a default boundary condition dictionary for [`TreeMesh`](@ref)es in 3D -that use the standard boundary naming convention. +that uses the standard boundary naming convention. This function applies the same boundary condition to all standard boundaries: - `:x_neg`: negative x-direction boundary - `:x_pos`: positive x-direction boundary From 6e02a05a90b7f7fdc68862fe8430c2d10077ef78 Mon Sep 17 00:00:00 2001 From: Vincent Marks Date: Wed, 29 Oct 2025 20:42:28 +0100 Subject: [PATCH 41/45] Apply suggestions from code review Co-authored-by: Joshua Lampert <51029046+JoshuaLampert@users.noreply.github.com> --- .../elixir_euler_source_terms_nonperiodic.jl | 2 +- ...elixir_euler_source_terms_nonperiodic_curved.jl | 2 +- .../elixir_euler_source_terms_nonperiodic.jl | 2 +- .../tree_3d_dgsem/elixir_advection_extended.jl | 2 +- src/solvers/default_boundary_conditions.jl | 14 ++++++++++++++ 5 files changed, 18 insertions(+), 4 deletions(-) diff --git a/examples/structured_1d_dgsem/elixir_euler_source_terms_nonperiodic.jl b/examples/structured_1d_dgsem/elixir_euler_source_terms_nonperiodic.jl index 0d625b31086..11fc7ca6074 100644 --- a/examples/structured_1d_dgsem/elixir_euler_source_terms_nonperiodic.jl +++ b/examples/structured_1d_dgsem/elixir_euler_source_terms_nonperiodic.jl @@ -24,7 +24,7 @@ f2() = SVector(2.0) mesh = StructuredMesh((16,), (f1, f2), periodicity = false) # you can either use a single function to impose the BCs weakly in all -# 1*ndims == 2 directions or you can pass a tuple containing BCs for +# 2*ndims == 2 directions or you can pass a tuple containing BCs for # each direction # Assign a single boundary condition to all boundaries boundary_condition = BoundaryConditionDirichlet(initial_condition) diff --git a/examples/structured_3d_dgsem/elixir_euler_source_terms_nonperiodic_curved.jl b/examples/structured_3d_dgsem/elixir_euler_source_terms_nonperiodic_curved.jl index b55df45c77e..b4be67074eb 100644 --- a/examples/structured_3d_dgsem/elixir_euler_source_terms_nonperiodic_curved.jl +++ b/examples/structured_3d_dgsem/elixir_euler_source_terms_nonperiodic_curved.jl @@ -48,7 +48,7 @@ cells_per_dimension = (4, 4, 4) mesh = StructuredMesh(cells_per_dimension, mapping, periodicity = false) # you can either use a single function to impose the BCs weakly in all -# 2*ndims == 4 directions or you can pass a tuple containing BCs for each direction +# 2*ndims == 6 directions or you can pass a tuple containing BCs for each direction # Assign a single boundary condition to all boundaries boundary_condition = BoundaryConditionDirichlet(initial_condition) boundary_conditions = boundary_condition_default(mesh, boundary_condition) diff --git a/examples/tree_1d_dgsem/elixir_euler_source_terms_nonperiodic.jl b/examples/tree_1d_dgsem/elixir_euler_source_terms_nonperiodic.jl index 614aad34e41..667c9d76690 100644 --- a/examples/tree_1d_dgsem/elixir_euler_source_terms_nonperiodic.jl +++ b/examples/tree_1d_dgsem/elixir_euler_source_terms_nonperiodic.jl @@ -25,7 +25,7 @@ mesh = TreeMesh(coordinates_min, coordinates_max, periodicity = false) # you can either use a single function to impose the BCs weakly in all -# 1*ndims == 2 directions or you can pass a tuple containing BCs for +# 2*ndims == 2 directions or you can pass a tuple containing BCs for # each direction # Assign a single boundary condition to all boundaries boundary_condition = BoundaryConditionDirichlet(initial_condition) diff --git a/examples/tree_3d_dgsem/elixir_advection_extended.jl b/examples/tree_3d_dgsem/elixir_advection_extended.jl index d6655ca7a7f..4b3a86cee72 100644 --- a/examples/tree_3d_dgsem/elixir_advection_extended.jl +++ b/examples/tree_3d_dgsem/elixir_advection_extended.jl @@ -22,7 +22,7 @@ mesh = TreeMesh(coordinates_min, coordinates_max, periodicity = true) # you can either use a single function to impose the BCs weakly in all -# 1*ndims == 2 directions or you can pass a tuple containing BCs for +# 2*ndims == 6 directions or you can pass a tuple containing BCs for # each direction # Note: "boundary_condition_periodic" indicates that it is a periodic boundary and can be omitted on # fully periodic domains. Here, however, it is included to allow easy override during testing diff --git a/src/solvers/default_boundary_conditions.jl b/src/solvers/default_boundary_conditions.jl index f7caf119fda..07292147598 100644 --- a/src/solvers/default_boundary_conditions.jl +++ b/src/solvers/default_boundary_conditions.jl @@ -34,8 +34,10 @@ This function applies the same boundary condition to all standard boundaries: - `:y_pos`: positive y-direction boundary - `:z_neg`: negative z-direction boundary - `:z_pos`: positive z-direction boundary + # Arguments - `boundary_condition`: The boundary condition function to apply to all boundaries + # Returns - `Dict{Symbol, Any}`: Dictionary mapping boundary names to the boundary condition """ @@ -56,8 +58,10 @@ that uses the standard boundary naming convention. This function applies the same boundary condition to all standard boundaries: - `:x_neg`: negative x-direction boundary - `:x_pos`: positive x-direction boundary + # Arguments - `boundary_condition`: The boundary condition function to apply to all boundaries + # Returns - Named tuple mapping boundary names to the boundary condition """ @@ -76,8 +80,10 @@ This function applies the same boundary condition to all standard boundaries: - `:x_pos`: positive x-direction boundary - `:y_neg`: negative y-direction boundary - `:y_pos`: positive y-direction boundary + # Arguments - `boundary_condition`: The boundary condition function to apply to all boundaries + # Returns - Named tuple mapping boundary names to the boundary condition """ @@ -100,8 +106,10 @@ This function applies the same boundary condition to all standard boundaries: - `:y_pos`: positive y-direction boundary - `:z_neg`: negative z-direction boundary - `:z_pos`: positive z-direction boundary + # Arguments - `boundary_condition`: The boundary condition function to apply to all boundaries + # Returns - Named tuple mapping boundary names to the boundary condition """ @@ -123,8 +131,10 @@ that uses the standard boundary naming convention. This function applies the same boundary condition to all standard boundaries: - `:x_neg`: negative x-direction boundary - `:x_pos`: positive x-direction boundary + # Arguments - `boundary_condition`: The boundary condition function to apply to all boundaries + # Returns - Named tuple mapping boundary names to the boundary condition """ @@ -143,8 +153,10 @@ This function applies the same boundary condition to all standard boundaries: - `:x_pos`: positive x-direction boundary - `:y_neg`: negative y-direction boundary - `:y_pos`: positive y-direction boundary + # Arguments - `boundary_condition`: The boundary condition function to apply to all boundaries + # Returns - Named tuple mapping boundary names to the boundary condition """ @@ -167,8 +179,10 @@ This function applies the same boundary condition to all standard boundaries: - `:y_pos`: positive y-direction boundary - `:z_neg`: negative z-direction boundary - `:z_pos`: positive z-direction boundary + # Arguments - `boundary_condition`: The boundary condition function to apply to all boundaries + # Returns - Named tuple mapping boundary names to the boundary condition """ From 8537ca2c78f2bb10a64c6ab3164a2f9238b8302e Mon Sep 17 00:00:00 2001 From: vincmarks Date: Wed, 29 Oct 2025 21:39:35 +0100 Subject: [PATCH 42/45] using `mesh` as an argument in the docstrings --- src/solvers/default_boundary_conditions.jl | 47 +++++++++++++--------- 1 file changed, 27 insertions(+), 20 deletions(-) diff --git a/src/solvers/default_boundary_conditions.jl b/src/solvers/default_boundary_conditions.jl index 07292147598..7b26fbfc19a 100644 --- a/src/solvers/default_boundary_conditions.jl +++ b/src/solvers/default_boundary_conditions.jl @@ -11,16 +11,17 @@ This function applies the same boundary condition to all standard boundaries: # Arguments - `boundary_condition`: The boundary condition function to apply to all boundaries +- `mesh`: The [`P4estMesh`](@ref) in 2D for which boundary conditions are created # Returns - `Dict{Symbol, Any}`: Dictionary mapping boundary names to the boundary condition """ -function boundary_condition_default(mesh::P4estMesh{2,2}, boundary_condition) - return Dict(:x_neg => boundary_condition, - :y_neg => boundary_condition, - :y_pos => boundary_condition, - :x_pos => boundary_condition) -end +function boundary_condition_default(mesh::P4estMesh{2, 2}, boundary_condition) + return Dict(:x_neg => boundary_condition, + :y_neg => boundary_condition, + :y_pos => boundary_condition, + :x_pos => boundary_condition) +end """ boundary_condition_default(mesh::P4estMesh{3,3}, boundary_condition) @@ -37,17 +38,18 @@ This function applies the same boundary condition to all standard boundaries: # Arguments - `boundary_condition`: The boundary condition function to apply to all boundaries +- `mesh`: The [`P4estMesh`](@ref) in 3D for which boundary conditions are created # Returns - `Dict{Symbol, Any}`: Dictionary mapping boundary names to the boundary condition """ -function boundary_condition_default(mesh::P4estMesh{3,3}, boundary_condition) - return Dict(:x_neg => boundary_condition, - :x_pos => boundary_condition, - :y_neg => boundary_condition, - :y_pos => boundary_condition, - :z_neg => boundary_condition, - :z_pos => boundary_condition) +function boundary_condition_default(mesh::P4estMesh{3, 3}, boundary_condition) + return Dict(:x_neg => boundary_condition, + :x_pos => boundary_condition, + :y_neg => boundary_condition, + :y_pos => boundary_condition, + :z_neg => boundary_condition, + :z_pos => boundary_condition) end """ @@ -61,6 +63,7 @@ This function applies the same boundary condition to all standard boundaries: # Arguments - `boundary_condition`: The boundary condition function to apply to all boundaries +- `mesh`: The [`StructuredMesh`](@ref) in 1D for which boundary conditions are created # Returns - Named tuple mapping boundary names to the boundary condition @@ -83,11 +86,12 @@ This function applies the same boundary condition to all standard boundaries: # Arguments - `boundary_condition`: The boundary condition function to apply to all boundaries +- `mesh`: The [`StructuredMesh`](@ref) in 2D for which boundary conditions are created # Returns - Named tuple mapping boundary names to the boundary condition """ -function boundary_condition_default(mesh::StructuredMesh{2}, boundary_condition) +function boundary_condition_default(mesh::StructuredMesh{2}, boundary_condition) return (x_neg = boundary_condition, x_pos = boundary_condition, y_neg = boundary_condition, @@ -109,19 +113,19 @@ This function applies the same boundary condition to all standard boundaries: # Arguments - `boundary_condition`: The boundary condition function to apply to all boundaries +- `mesh`: The [`StructuredMesh`](@ref) in 3D for which boundary conditions are created # Returns - Named tuple mapping boundary names to the boundary condition """ -function boundary_condition_default(mesh::StructuredMesh{3}, boundary_condition) +function boundary_condition_default(mesh::StructuredMesh{3}, boundary_condition) return (x_neg = boundary_condition, x_pos = boundary_condition, y_neg = boundary_condition, y_pos = boundary_condition, z_neg = boundary_condition, z_pos = boundary_condition) - -end +end """ boundary_condition_default(mesh::TreeMesh1D, boundary_condition) @@ -134,6 +138,7 @@ This function applies the same boundary condition to all standard boundaries: # Arguments - `boundary_condition`: The boundary condition function to apply to all boundaries +- `mesh`: The [`TreeMesh`](@ref) in 1D for which boundary conditions are created # Returns - Named tuple mapping boundary names to the boundary condition @@ -156,11 +161,12 @@ This function applies the same boundary condition to all standard boundaries: # Arguments - `boundary_condition`: The boundary condition function to apply to all boundaries +- `mesh`: The [`TreeMesh`](@ref) in 2D for which boundary conditions are created # Returns - Named tuple mapping boundary names to the boundary condition """ -function boundary_condition_default(mesh::TreeMesh2D, boundary_condition) +function boundary_condition_default(mesh::TreeMesh2D, boundary_condition) return (x_neg = boundary_condition, x_pos = boundary_condition, y_neg = boundary_condition, @@ -182,15 +188,16 @@ This function applies the same boundary condition to all standard boundaries: # Arguments - `boundary_condition`: The boundary condition function to apply to all boundaries +- `mesh`: The [`TreeMesh`](@ref) in 3D for which boundary conditions are created # Returns - Named tuple mapping boundary names to the boundary condition """ -function boundary_condition_default(mesh::TreeMesh3D, boundary_condition) +function boundary_condition_default(mesh::TreeMesh3D, boundary_condition) return (x_neg = boundary_condition, x_pos = boundary_condition, y_neg = boundary_condition, y_pos = boundary_condition, z_neg = boundary_condition, z_pos = boundary_condition) -end \ No newline at end of file +end From 275c30a158464e6a649a141584076d9671de1a17 Mon Sep 17 00:00:00 2001 From: Joshua Lampert <51029046+JoshuaLampert@users.noreply.github.com> Date: Wed, 29 Oct 2025 21:55:27 +0100 Subject: [PATCH 43/45] format --- src/Trixi.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/Trixi.jl b/src/Trixi.jl index fe495a71da1..b86d9c9ff85 100644 --- a/src/Trixi.jl +++ b/src/Trixi.jl @@ -141,7 +141,7 @@ include("auxiliary/t8code.jl") include("equations/equations.jl") include("meshes/meshes.jl") include("solvers/solvers.jl") -include("solvers/default_boundary_conditions.jl") +include("solvers/default_boundary_conditions.jl") include("equations/equations_parabolic.jl") # these depend on parabolic solver types include("semidiscretization/semidiscretization.jl") include("semidiscretization/semidiscretization_hyperbolic.jl") From 4bd67435c0d2ec7c548b606859ae56b82c9f8fa7 Mon Sep 17 00:00:00 2001 From: vincmarks Date: Thu, 30 Oct 2025 19:08:45 +0100 Subject: [PATCH 44/45] changing the file name to boundary_condition_default.jl --- src/Trixi.jl | 2 +- ...ult_boundary_conditions.jl => boundary_condition_default.jl} | 0 2 files changed, 1 insertion(+), 1 deletion(-) rename src/solvers/{default_boundary_conditions.jl => boundary_condition_default.jl} (100%) diff --git a/src/Trixi.jl b/src/Trixi.jl index b86d9c9ff85..21a86f2a57f 100644 --- a/src/Trixi.jl +++ b/src/Trixi.jl @@ -141,7 +141,7 @@ include("auxiliary/t8code.jl") include("equations/equations.jl") include("meshes/meshes.jl") include("solvers/solvers.jl") -include("solvers/default_boundary_conditions.jl") +include("solvers/boundary_condition_default.jl") include("equations/equations_parabolic.jl") # these depend on parabolic solver types include("semidiscretization/semidiscretization.jl") include("semidiscretization/semidiscretization_hyperbolic.jl") diff --git a/src/solvers/default_boundary_conditions.jl b/src/solvers/boundary_condition_default.jl similarity index 100% rename from src/solvers/default_boundary_conditions.jl rename to src/solvers/boundary_condition_default.jl From dcba76a44b2d9041c32032da62c132a98f007a52 Mon Sep 17 00:00:00 2001 From: Vincent Marks Date: Thu, 30 Oct 2025 19:11:22 +0100 Subject: [PATCH 45/45] Apply suggestions from code review Co-authored-by: Daniel Doehring --- docs/src/meshes/p4est_mesh.md | 2 +- docs/src/meshes/structured_mesh.md | 2 +- docs/src/meshes/tree_mesh.md | 2 +- 3 files changed, 3 insertions(+), 3 deletions(-) diff --git a/docs/src/meshes/p4est_mesh.md b/docs/src/meshes/p4est_mesh.md index 2cc604060d7..69db1f190e2 100644 --- a/docs/src/meshes/p4est_mesh.md +++ b/docs/src/meshes/p4est_mesh.md @@ -941,5 +941,5 @@ reinitialize_boundaries!(semi.boundary_conditions, cache) # Needs to be called a This code could then be placed in the [`resize!`](https://github.com/trixi-framework/Trixi.jl/blob/eaeb04113523500ed831e3ab459694f12f7a49ea/src/time_integration/methods_2N.jl#L251-L255) function of a corresponding multirate integrator to ensure load-balancing for simulations involving AMR. ### Boundary conditions -For [`P4estMesh`](@ref)es, boundary conditions are defined and stored in dictionaries (see, for example, `examples/p4est_2d_dgsem/elixir_advection_diffusion_nonperiodic_amr.jl`). +For [`P4estMesh`](@ref)es, boundary conditions are defined and stored in [dictionaries](https://docs.julialang.org/en/v1/base/collections/#Base.Dict) (see, for example, `examples/p4est_2d_dgsem/elixir_advection_diffusion_nonperiodic_amr.jl`). If you want to apply the same boundary condition to all faces of the mesh, you can use the `boundary_condition_default(mesh, boundary_condition)` function, as demonstrated in `examples/p4est_2d_dgsem/elixir_advection_diffusion_nonperiodic_amr.jl` and `examples/p4est_3d_dgsem/elixir_euler_source_terms_nonperiodic.jl`. diff --git a/docs/src/meshes/structured_mesh.md b/docs/src/meshes/structured_mesh.md index 85a50f37ad3..7308a91956e 100644 --- a/docs/src/meshes/structured_mesh.md +++ b/docs/src/meshes/structured_mesh.md @@ -13,5 +13,5 @@ the re-use of existing functionality for the [`TreeMesh`](@ref) but is usually less efficient, cf. [PR #550](https://github.com/trixi-framework/Trixi.jl/pull/550). ### Boundary conditions -For [`StructuredMesh`](@ref)es, boundary conditions are defined and stored in named tuples (see, for example, `examples/structured_1d_dgsem/elixir_euler_source_terms_nonperiodic.jl`). +For [`StructuredMesh`](@ref)es, boundary conditions are defined and stored in [named tuples](https://docs.julialang.org/en/v1/manual/functions/#Named-Tuples) (see, for example, `examples/structured_1d_dgsem/elixir_euler_source_terms_nonperiodic.jl`). If you want to apply the same boundary condition to all faces of the mesh, you can use the `boundary_condition_default(mesh, boundary_condition)` function, as demonstrated in `examples/structured_1d_dgsem/elixir_euler_source_terms_nonperiodic.jl`, `examples/structured_2d_dgsem/elixir_euler_rayleigh_taylor_instability.jl` and `examples/structured_3d_dgsem/elixir_euler_source_terms_nonperiodic_curved.jl`. diff --git a/docs/src/meshes/tree_mesh.md b/docs/src/meshes/tree_mesh.md index e1df92a64d0..77dc673f7e9 100644 --- a/docs/src/meshes/tree_mesh.md +++ b/docs/src/meshes/tree_mesh.md @@ -12,5 +12,5 @@ dispatching on the `orientation::Integer` as described in the ### Boundary conditions -For [`TreeMesh`](@ref)es, boundary conditions are defined and stored in named tuples (see, for example, `examples/tree_1d_dgsem/elixir_euler_source_terms_nonperiodic.jl`). +For [`TreeMesh`](@ref)es, boundary conditions are defined and stored in [named tuples](https://docs.julialang.org/en/v1/manual/functions/#Named-Tuples) (see, for example, `examples/tree_1d_dgsem/elixir_euler_source_terms_nonperiodic.jl`). If you want to apply the same boundary condition to all faces of the mesh, you can use the `boundary_condition_default(mesh, boundary_condition)` function, as demonstrated in `examples/tree_1d_dgsem/elixir_euler_source_terms_nonperiodic.jl`, `examples/tree_2d_dgsem/elixir_euler_source_terms_nonperiodic.jl` and `examples/tree_3d_dgsem/elixir_advection_extended.jl`. \ No newline at end of file