From a291620d05a0bc6b985065e209ca712f2cd98c6c Mon Sep 17 00:00:00 2001 From: Michael Schlottke-Lakemper Date: Wed, 6 Sep 2023 12:04:08 +0200 Subject: [PATCH 01/56] Add StructuredMeshView as proxy between solver and actual StructuredMesh --- .../elixir_advection_smview.jl | 64 +++++++++++++++++++ src/Trixi.jl | 3 +- src/callbacks_step/analysis_dg2d.jl | 20 +++--- src/meshes/meshes.jl | 1 + src/meshes/structured_mesh_view.jl | 19 ++++++ src/solvers/dg.jl | 4 +- src/solvers/dgsem_structured/containers.jl | 2 +- src/solvers/dgsem_structured/containers_2d.jl | 8 ++- src/solvers/dgsem_structured/dg.jl | 4 +- src/solvers/dgsem_structured/dg_2d.jl | 19 +++--- src/solvers/dgsem_tree/dg_2d.jl | 7 +- 11 files changed, 121 insertions(+), 30 deletions(-) create mode 100644 examples/structured_2d_dgsem/elixir_advection_smview.jl create mode 100644 src/meshes/structured_mesh_view.jl diff --git a/examples/structured_2d_dgsem/elixir_advection_smview.jl b/examples/structured_2d_dgsem/elixir_advection_smview.jl new file mode 100644 index 00000000000..b12a4c4bf75 --- /dev/null +++ b/examples/structured_2d_dgsem/elixir_advection_smview.jl @@ -0,0 +1,64 @@ +# The same setup as tree_2d_dgsem/elixir_advection_basic.jl +# to verify the StructuredMesh implementation against TreeMesh + +using OrdinaryDiffEq +using Trixi + +############################################################################### +# semidiscretization of the linear advection equation + +advection_velocity = (0.2, -0.7) +equations = LinearScalarAdvectionEquation2D(advection_velocity) + +# Create DG solver with polynomial degree = 3 and (local) Lax-Friedrichs/Rusanov flux as surface flux +solver = DGSEM(polydeg=3, surface_flux=flux_lax_friedrichs) + +coordinates_min = (-1.0, -1.0) # minimum coordinates (min(x), min(y)) +coordinates_max = ( 1.0, 1.0) # maximum coordinates (max(x), max(y)) + +cells_per_dimension = (16, 16) + +# Create curved mesh with 16 x 16 elements +# mesh = StructuredMesh(cells_per_dimension, coordinates_min, coordinates_max) +parent = StructuredMesh(cells_per_dimension, coordinates_min, coordinates_max) +mesh = StructuredMeshView(parent) + +# A semidiscretization collects data structures and functions for the spatial discretization +semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition_convergence_test, solver) + + +############################################################################### +# ODE solvers, callbacks etc. + +# Create ODE problem with time span from 0.0 to 1.0 +ode = semidiscretize(semi, (0.0, 1.0)); + +# At the beginning of the main loop, the SummaryCallback prints a summary of the simulation setup +# and resets the timers +summary_callback = SummaryCallback() + +# The AnalysisCallback allows to analyse the solution in regular intervals and prints the results +analysis_callback = AnalysisCallback(semi, interval=10) + +# The SaveSolutionCallback allows to save the solution to a file in regular intervals +# save_solution = SaveSolutionCallback(interval=100, +# solution_variables=cons2prim) + +# The StepsizeCallback handles the re-calculation of the maximum Δt after each time step +# stepsize_callback = StepsizeCallback(cfl=1.6) + +# Create a CallbackSet to collect all callbacks such that they can be passed to the ODE solver +# callbacks = CallbackSet(summary_callback, analysis_callback, save_solution, stepsize_callback) +callbacks = CallbackSet(summary_callback, analysis_callback) + + +############################################################################### +# run the simulation + +# OrdinaryDiffEq's `solve` method evolves the solution in time and executes the passed callbacks +sol = solve(ode, CarpenterKennedy2N54(williamson_condition=false), + dt=5.0e-2, # solve needs some value here but it will be overwritten by the stepsize_callback + save_everystep=false, callback=callbacks); + +# Print the timer summary +summary_callback() diff --git a/src/Trixi.jl b/src/Trixi.jl index ec4d20558e5..dcf33c1ca03 100644 --- a/src/Trixi.jl +++ b/src/Trixi.jl @@ -213,7 +213,8 @@ export entropy, energy_total, energy_kinetic, energy_internal, energy_magnetic, export lake_at_rest_error export ncomponents, eachcomponent -export TreeMesh, StructuredMesh, UnstructuredMesh2D, P4estMesh, T8codeMesh +export TreeMesh, StructuredMesh, UnstructuredMesh2D, P4estMesh, T8codeMesh, + StructuredMeshView export DG, DGSEM, LobattoLegendreBasis, diff --git a/src/callbacks_step/analysis_dg2d.jl b/src/callbacks_step/analysis_dg2d.jl index aecabf0e4b7..54503a51da6 100644 --- a/src/callbacks_step/analysis_dg2d.jl +++ b/src/callbacks_step/analysis_dg2d.jl @@ -30,7 +30,7 @@ function create_cache_analysis(analyzer, mesh::TreeMesh{2}, end function create_cache_analysis(analyzer, - mesh::Union{StructuredMesh{2}, UnstructuredMesh2D, + mesh::Union{StructuredMesh{2}, StructuredMeshView{2}, UnstructuredMesh2D, P4estMesh{2}, T8codeMesh{2}}, equations, dg::DG, cache, RealT, uEltype) @@ -107,8 +107,9 @@ function calc_error_norms(func, u, t, analyzer, end function calc_error_norms(func, u, t, analyzer, - mesh::Union{StructuredMesh{2}, UnstructuredMesh2D, - P4estMesh{2}, T8codeMesh{2}}, equations, + mesh::Union{StructuredMesh{2}, StructuredMeshView{2}, + UnstructuredMesh2D, P4estMesh{2}, T8codeMesh{2}}, + equations, initial_condition, dg::DGSEM, cache, cache_analysis) @unpack vandermonde, weights = analyzer @unpack node_coordinates, inverse_jacobian = cache.elements @@ -175,8 +176,9 @@ function integrate_via_indices(func::Func, u, end function integrate_via_indices(func::Func, u, - mesh::Union{StructuredMesh{2}, UnstructuredMesh2D, - P4estMesh{2}, T8codeMesh{2}}, equations, + mesh::Union{StructuredMesh{2}, StructuredMeshView{2}, + UnstructuredMesh2D, P4estMesh{2}, T8codeMesh{2}}, + equations, dg::DGSEM, cache, args...; normalize = true) where {Func} @unpack weights = dg.basis @@ -203,8 +205,8 @@ function integrate_via_indices(func::Func, u, end function integrate(func::Func, u, - mesh::Union{TreeMesh{2}, StructuredMesh{2}, UnstructuredMesh2D, - P4estMesh{2}, T8codeMesh{2}}, + mesh::Union{TreeMesh{2}, StructuredMesh{2}, StructuredMeshView{2}, + UnstructuredMesh2D, P4estMesh{2}, T8codeMesh{2}}, equations, dg::DG, cache; normalize = true) where {Func} integrate_via_indices(u, mesh, equations, dg, cache; normalize = normalize) do u, i, j, element, equations, dg @@ -232,8 +234,8 @@ function integrate(func::Func, u, end function analyze(::typeof(entropy_timederivative), du, u, t, - mesh::Union{TreeMesh{2}, StructuredMesh{2}, UnstructuredMesh2D, - P4estMesh{2}, T8codeMesh{2}}, + mesh::Union{TreeMesh{2}, StructuredMesh{2}, StructuredMeshView{2}, + UnstructuredMesh2D, P4estMesh{2}, T8codeMesh{2}}, equations, dg::DG, cache) # Calculate ∫(∂S/∂u ⋅ ∂u/∂t)dΩ integrate_via_indices(u, mesh, equations, dg, cache, diff --git a/src/meshes/meshes.jl b/src/meshes/meshes.jl index ed2158b169a..4d6016e5564 100644 --- a/src/meshes/meshes.jl +++ b/src/meshes/meshes.jl @@ -7,6 +7,7 @@ include("tree_mesh.jl") include("structured_mesh.jl") +include("structured_mesh_view.jl") include("surface_interpolant.jl") include("unstructured_mesh.jl") include("face_interpolant.jl") diff --git a/src/meshes/structured_mesh_view.jl b/src/meshes/structured_mesh_view.jl new file mode 100644 index 00000000000..53e0a7da96f --- /dev/null +++ b/src/meshes/structured_mesh_view.jl @@ -0,0 +1,19 @@ +mutable struct StructuredMeshView{NDIMS, RealT <: Real} <: AbstractMesh{NDIMS} + parent::StructuredMesh{NDIMS, RealT} + mapping::Any # Not relevant for performance +end + +function StructuredMeshView(parent::StructuredMesh{NDIMS, RealT}) where {NDIMS, RealT} + return StructuredMeshView{NDIMS, RealT}(parent, parent.mapping) +end + +# Check if mesh is periodic +isperiodic(mesh::StructuredMeshView) = all(mesh.parent.periodicity) +isperiodic(mesh::StructuredMeshView, dimension) = mesh.parent.periodicity[dimension] + +@inline Base.ndims(::StructuredMeshView{NDIMS}) where {NDIMS} = NDIMS +@inline Base.real(::StructuredMeshView{NDIMS, RealT}) where {NDIMS, RealT} = RealT +Base.size(mesh::StructuredMeshView) = mesh.parent.cells_per_dimension +Base.size(mesh::StructuredMeshView, i) = mesh.parent.cells_per_dimension[i] +Base.axes(mesh::StructuredMeshView) = map(Base.OneTo, mesh.parent.cells_per_dimension) +Base.axes(mesh::StructuredMeshView, i) = Base.OneTo(mesh.parent.cells_per_dimension[i]) diff --git a/src/solvers/dg.jl b/src/solvers/dg.jl index 36bbc6de361..4883f34ed16 100644 --- a/src/solvers/dg.jl +++ b/src/solvers/dg.jl @@ -403,8 +403,8 @@ function get_element_variables!(element_variables, u, mesh, equations, dg::DG, c dg, cache) end -const MeshesDGSEM = Union{TreeMesh, StructuredMesh, UnstructuredMesh2D, P4estMesh, - T8codeMesh} +const MeshesDGSEM = Union{TreeMesh, StructuredMesh, StructuredMeshView, UnstructuredMesh2D, + P4estMesh, T8codeMesh} @inline function ndofs(mesh::MeshesDGSEM, dg::DG, cache) nelements(cache.elements) * nnodes(dg)^ndims(mesh) diff --git a/src/solvers/dgsem_structured/containers.jl b/src/solvers/dgsem_structured/containers.jl index 41eabf7c6bf..14f49641b55 100644 --- a/src/solvers/dgsem_structured/containers.jl +++ b/src/solvers/dgsem_structured/containers.jl @@ -23,7 +23,7 @@ struct ElementContainer{NDIMS, RealT <: Real, uEltype <: Real, NDIMSP1, NDIMSP2, end # Create element container and initialize element data -function init_elements(mesh::StructuredMesh{NDIMS, RealT}, +function init_elements(mesh::Union{StructuredMesh{NDIMS, RealT},StructuredMeshView{NDIMS, RealT}}, equations::AbstractEquations, basis, ::Type{uEltype}) where {NDIMS, RealT <: Real, uEltype <: Real} diff --git a/src/solvers/dgsem_structured/containers_2d.jl b/src/solvers/dgsem_structured/containers_2d.jl index fb6db48e0a5..fafaea7ff8e 100644 --- a/src/solvers/dgsem_structured/containers_2d.jl +++ b/src/solvers/dgsem_structured/containers_2d.jl @@ -6,7 +6,8 @@ #! format: noindent # Initialize data structures in element container -function init_elements!(elements, mesh::StructuredMesh{2}, basis::LobattoLegendreBasis) +function init_elements!(elements, mesh::Union{StructuredMesh{2},StructuredMeshView{2}}, + basis::LobattoLegendreBasis) @unpack node_coordinates, left_neighbors, jacobian_matrix, contravariant_vectors, inverse_jacobian = elements @@ -35,7 +36,7 @@ end # `mesh.mapping` is passed as an additional argument for type stability (function barrier) function calc_node_coordinates!(node_coordinates, element, cell_x, cell_y, mapping, - mesh::StructuredMesh{2}, + mesh::Union{StructuredMesh{2},StructuredMeshView{2}}, basis::LobattoLegendreBasis) @unpack nodes = basis @@ -148,7 +149,8 @@ function calc_inverse_jacobian!(inverse_jacobian::AbstractArray{<:Any, 3}, eleme end # Save id of left neighbor of every element -function initialize_left_neighbor_connectivity!(left_neighbors, mesh::StructuredMesh{2}, +function initialize_left_neighbor_connectivity!(left_neighbors, + mesh::Union{StructuredMesh{2},StructuredMeshView{2}}, linear_indices) # Neighbors in x-direction for cell_y in 1:size(mesh, 2) diff --git a/src/solvers/dgsem_structured/dg.jl b/src/solvers/dgsem_structured/dg.jl index 5cf4c4ef78c..571add15887 100644 --- a/src/solvers/dgsem_structured/dg.jl +++ b/src/solvers/dgsem_structured/dg.jl @@ -8,8 +8,8 @@ # This method is called when a SemidiscretizationHyperbolic is constructed. # It constructs the basic `cache` used throughout the simulation to compute # the RHS etc. -function create_cache(mesh::StructuredMesh, equations::AbstractEquations, dg::DG, ::Any, - ::Type{uEltype}) where {uEltype <: Real} +function create_cache(mesh::Union{StructuredMesh,StructuredMeshView}, + equations::AbstractEquations, dg::DG, ::Any, ::Type{uEltype}) where {uEltype <: Real} elements = init_elements(mesh, equations, dg.basis, uEltype) cache = (; elements) diff --git a/src/solvers/dgsem_structured/dg_2d.jl b/src/solvers/dgsem_structured/dg_2d.jl index 3e8ce759b30..3feb17f874d 100644 --- a/src/solvers/dgsem_structured/dg_2d.jl +++ b/src/solvers/dgsem_structured/dg_2d.jl @@ -6,7 +6,7 @@ #! format: noindent function rhs!(du, u, t, - mesh::StructuredMesh{2}, equations, + mesh::Union{StructuredMesh{2},StructuredMeshView{2}}, equations, initial_condition, boundary_conditions, source_terms::Source, dg::DG, cache) where {Source} # Reset du @@ -51,8 +51,9 @@ end @inline function weak_form_kernel!(du, u, element, - mesh::Union{StructuredMesh{2}, UnstructuredMesh2D, - P4estMesh{2}, T8codeMesh{2}}, + mesh::Union{StructuredMesh{2}, StructuredMeshView{2}, + UnstructuredMesh2D, P4estMesh{2}, + T8codeMesh{2}}, nonconservative_terms::False, equations, dg::DGSEM, cache, alpha = true) # true * [some floating point value] == [exactly the same floating point value] @@ -380,7 +381,7 @@ end end function calc_interface_flux!(cache, u, - mesh::StructuredMesh{2}, + mesh::Union{StructuredMesh{2},StructuredMeshView{2}}, nonconservative_terms, # can be True/False equations, surface_integral, dg::DG) @unpack elements = cache @@ -409,7 +410,7 @@ end @inline function calc_interface_flux!(surface_flux_values, left_element, right_element, orientation, u, - mesh::StructuredMesh{2}, + mesh::Union{StructuredMesh{2},StructuredMeshView{2}}, nonconservative_terms::False, equations, surface_integral, dg::DG, cache) # This is slow for LSA, but for some reason faster for Euler (see #519) @@ -544,8 +545,8 @@ end # TODO: Taal dimension agnostic function calc_boundary_flux!(cache, u, t, boundary_condition::BoundaryConditionPeriodic, - mesh::StructuredMesh{2}, equations, surface_integral, - dg::DG) + mesh::Union{StructuredMesh{2},StructuredMeshView{2}}, + equations, surface_integral, dg::DG) @assert isperiodic(mesh) end @@ -609,8 +610,8 @@ function calc_boundary_flux!(cache, u, t, boundary_conditions::NamedTuple, end function apply_jacobian!(du, - mesh::Union{StructuredMesh{2}, UnstructuredMesh2D, - P4estMesh{2}, T8codeMesh{2}}, + mesh::Union{StructuredMesh{2}, StructuredMeshView{2}, + UnstructuredMesh2D, P4estMesh{2}, T8codeMesh{2}}, equations, dg::DG, cache) @unpack inverse_jacobian = cache.elements diff --git a/src/solvers/dgsem_tree/dg_2d.jl b/src/solvers/dgsem_tree/dg_2d.jl index c30d0a8e01a..7c7293f6107 100644 --- a/src/solvers/dgsem_tree/dg_2d.jl +++ b/src/solvers/dgsem_tree/dg_2d.jl @@ -180,8 +180,8 @@ end function calc_volume_integral!(du, u, mesh::Union{TreeMesh{2}, StructuredMesh{2}, - UnstructuredMesh2D, P4estMesh{2}, - T8codeMesh{2}}, + StructuredMeshView{2}, UnstructuredMesh2D, + P4estMesh{2}, T8codeMesh{2}}, nonconservative_terms, equations, volume_integral::VolumeIntegralWeakForm, dg::DGSEM, cache) @@ -1078,7 +1078,8 @@ end return nothing end -function calc_surface_integral!(du, u, mesh::Union{TreeMesh{2}, StructuredMesh{2}}, +function calc_surface_integral!(du, u, mesh::Union{TreeMesh{2}, StructuredMesh{2}, + StructuredMeshView{2}}, equations, surface_integral::SurfaceIntegralWeakForm, dg::DG, cache) @unpack boundary_interpolation = dg.basis From f3d20567dc2e6b76655ee417be5ada5d684e9efd Mon Sep 17 00:00:00 2001 From: Michael Schlottke-Lakemper Date: Wed, 6 Sep 2023 13:13:22 +0200 Subject: [PATCH 02/56] Enable StructuredMeshView on submesh --- .../elixir_advection_smview.jl | 12 ++-- src/meshes/structured_mesh_view.jl | 66 +++++++++++++++++-- src/solvers/dgsem_structured/containers_2d.jl | 2 +- 3 files changed, 68 insertions(+), 12 deletions(-) diff --git a/examples/structured_2d_dgsem/elixir_advection_smview.jl b/examples/structured_2d_dgsem/elixir_advection_smview.jl index b12a4c4bf75..23b65b18cc7 100644 --- a/examples/structured_2d_dgsem/elixir_advection_smview.jl +++ b/examples/structured_2d_dgsem/elixir_advection_smview.jl @@ -11,17 +11,19 @@ advection_velocity = (0.2, -0.7) equations = LinearScalarAdvectionEquation2D(advection_velocity) # Create DG solver with polynomial degree = 3 and (local) Lax-Friedrichs/Rusanov flux as surface flux -solver = DGSEM(polydeg=3, surface_flux=flux_lax_friedrichs) +solver = DGSEM(polydeg=2, surface_flux=flux_lax_friedrichs) -coordinates_min = (-1.0, -1.0) # minimum coordinates (min(x), min(y)) +coordinates_min = (-3.0, -1.0) # minimum coordinates (min(x), min(y)) +# coordinates_max = ( 1.0, 1.0) # maximum coordinates (max(x), max(y)) coordinates_max = ( 1.0, 1.0) # maximum coordinates (max(x), max(y)) -cells_per_dimension = (16, 16) +# cells_per_dimension = (16, 16) +cells_per_dimension = (32, 16) # Create curved mesh with 16 x 16 elements # mesh = StructuredMesh(cells_per_dimension, coordinates_min, coordinates_max) parent = StructuredMesh(cells_per_dimension, coordinates_min, coordinates_max) -mesh = StructuredMeshView(parent) +mesh = StructuredMeshView(parent; index_min = (17, 1), index_max = (32, 16)) # A semidiscretization collects data structures and functions for the spatial discretization semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition_convergence_test, solver) @@ -60,5 +62,7 @@ sol = solve(ode, CarpenterKennedy2N54(williamson_condition=false), dt=5.0e-2, # solve needs some value here but it will be overwritten by the stepsize_callback save_everystep=false, callback=callbacks); +errors_ref = (l2 = [8.312427642603623e-6], linf = [6.626865824577166e-5]) + # Print the timer summary summary_callback() diff --git a/src/meshes/structured_mesh_view.jl b/src/meshes/structured_mesh_view.jl index 53e0a7da96f..eb2fb1d2419 100644 --- a/src/meshes/structured_mesh_view.jl +++ b/src/meshes/structured_mesh_view.jl @@ -1,19 +1,71 @@ mutable struct StructuredMeshView{NDIMS, RealT <: Real} <: AbstractMesh{NDIMS} parent::StructuredMesh{NDIMS, RealT} mapping::Any # Not relevant for performance + index_min::NTuple{NDIMS, Int} + index_max::NTuple{NDIMS, Int} +end + +function StructuredMeshView(parent::StructuredMesh{NDIMS, RealT}; + index_min, index_max) where {NDIMS, RealT} + @assert index_min <= index_max + @assert all(index_min .> 0) + @assert index_max <= size(parent) + + return StructuredMeshView{NDIMS, RealT}(parent, parent.mapping, index_min, index_max) end function StructuredMeshView(parent::StructuredMesh{NDIMS, RealT}) where {NDIMS, RealT} - return StructuredMeshView{NDIMS, RealT}(parent, parent.mapping) + return StructuredMeshView(parent, ntuple(_ -> 1, Val(NDIMS)), size(parent)) end # Check if mesh is periodic -isperiodic(mesh::StructuredMeshView) = all(mesh.parent.periodicity) -isperiodic(mesh::StructuredMeshView, dimension) = mesh.parent.periodicity[dimension] +function isperiodic(mesh::StructuredMeshView) + @unpack parent = mesh + return isperiodic(parent) && size(parent) == size(mesh) +end +function isperiodic(mesh::StructuredMeshView, dimension) + @unpack parent, index_min, index_max = mesh + return (isperiodic(parent, dimension) && + index_min[dimension] == 1 && + index_max[dimension] == size(parent, dimension)) +end @inline Base.ndims(::StructuredMeshView{NDIMS}) where {NDIMS} = NDIMS @inline Base.real(::StructuredMeshView{NDIMS, RealT}) where {NDIMS, RealT} = RealT -Base.size(mesh::StructuredMeshView) = mesh.parent.cells_per_dimension -Base.size(mesh::StructuredMeshView, i) = mesh.parent.cells_per_dimension[i] -Base.axes(mesh::StructuredMeshView) = map(Base.OneTo, mesh.parent.cells_per_dimension) -Base.axes(mesh::StructuredMeshView, i) = Base.OneTo(mesh.parent.cells_per_dimension[i]) +function Base.size(mesh::StructuredMeshView) + @unpack index_min, index_max = mesh + return index_max .- index_min .+ 1 +end +function Base.size(mesh::StructuredMeshView, i) + @unpack index_min, index_max = mesh + return index_max[i] - index_min[i] + 1 +end +Base.axes(mesh::StructuredMeshView) = map(Base.OneTo, size(mesh)) +Base.axes(mesh::StructuredMeshView, i) = Base.OneTo(size(mesh, i)) + + +function calc_node_coordinates!(node_coordinates, element, + cell_x, cell_y, mapping, + mesh::StructuredMeshView{2}, + basis::LobattoLegendreBasis) + @unpack nodes = basis + @unpack parent, index_min, index_max = mesh + + # Get cell length in reference mesh + dx = 2 / size(parent, 1) + dy = 2 / size(parent, 2) + + # Calculate index offsets with respect to parent + parent_offset_x = index_min[1] - 1 + parent_offset_y = index_min[2] - 1 + + # Calculate node coordinates of reference mesh + cell_x_offset = -1 + (cell_x - 1 + parent_offset_x) * dx + dx / 2 + cell_y_offset = -1 + (cell_y - 1 + parent_offset_y) * dy + dy / 2 + + for j in eachnode(basis), i in eachnode(basis) + # node_coordinates are the mapped reference node_coordinates + node_coordinates[:, i, j, element] .= mapping(cell_x_offset + dx / 2 * nodes[i], + cell_y_offset + dy / 2 * nodes[j]) + end +end diff --git a/src/solvers/dgsem_structured/containers_2d.jl b/src/solvers/dgsem_structured/containers_2d.jl index fafaea7ff8e..9e77a9ce680 100644 --- a/src/solvers/dgsem_structured/containers_2d.jl +++ b/src/solvers/dgsem_structured/containers_2d.jl @@ -36,7 +36,7 @@ end # `mesh.mapping` is passed as an additional argument for type stability (function barrier) function calc_node_coordinates!(node_coordinates, element, cell_x, cell_y, mapping, - mesh::Union{StructuredMesh{2},StructuredMeshView{2}}, + mesh::StructuredMesh{2}, basis::LobattoLegendreBasis) @unpack nodes = basis From 37b9242ff7776c3c4784a32c52eadc79926d41bd Mon Sep 17 00:00:00 2001 From: SimonCan Date: Thu, 7 Sep 2023 09:18:34 +0100 Subject: [PATCH 03/56] Attempt to using coupled with meshviews. --- .../elixir_advection_smview.jl | 44 ++++++++++++++----- src/meshes/structured_mesh_view.jl | 11 ++++- .../semidiscretization_coupled.jl | 4 +- src/solvers/dgsem_structured/dg.jl | 4 +- src/solvers/dgsem_structured/dg_2d.jl | 2 +- 5 files changed, 48 insertions(+), 17 deletions(-) diff --git a/examples/structured_2d_dgsem/elixir_advection_smview.jl b/examples/structured_2d_dgsem/elixir_advection_smview.jl index 23b65b18cc7..7d28299ef79 100644 --- a/examples/structured_2d_dgsem/elixir_advection_smview.jl +++ b/examples/structured_2d_dgsem/elixir_advection_smview.jl @@ -10,23 +10,47 @@ using Trixi advection_velocity = (0.2, -0.7) equations = LinearScalarAdvectionEquation2D(advection_velocity) -# Create DG solver with polynomial degree = 3 and (local) Lax-Friedrichs/Rusanov flux as surface flux +# Create DG solver with polynomial degree = 2 and (local) Lax-Friedrichs/Rusanov flux as surface flux solver = DGSEM(polydeg=2, surface_flux=flux_lax_friedrichs) -coordinates_min = (-3.0, -1.0) # minimum coordinates (min(x), min(y)) -# coordinates_max = ( 1.0, 1.0) # maximum coordinates (max(x), max(y)) +coordinates_min = ( -1.0, -1.0) # maximum coordinates (max(x), max(y)) coordinates_max = ( 1.0, 1.0) # maximum coordinates (max(x), max(y)) -# cells_per_dimension = (16, 16) -cells_per_dimension = (32, 16) +cells_per_dimension = (16, 16) +# cells_per_dimension = (32, 16) # Create curved mesh with 16 x 16 elements -# mesh = StructuredMesh(cells_per_dimension, coordinates_min, coordinates_max) parent = StructuredMesh(cells_per_dimension, coordinates_min, coordinates_max) -mesh = StructuredMeshView(parent; index_min = (17, 1), index_max = (32, 16)) + +#mesh = StructuredMeshView(parent; index_min = (1, 1), index_max = (16, 16)) +#semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition_convergence_test, solver) + +mesh1 = StructuredMeshView(parent; index_min = (1, 1), index_max = (8, 16)) +mesh2 = StructuredMeshView(parent; index_min = (9, 1), index_max = (16, 16)) + +# Define the coupled boundary conditions +boundary_conditions1=( + # Connect left boundary with right boundary of left mesh + x_neg=boundary_condition_periodic, + x_pos=boundary_condition_periodic, +# x_pos=BoundaryConditionCoupled(2, (:begin, :i_forward), Float64), + y_neg=boundary_condition_periodic, + y_pos=boundary_condition_periodic) +boundary_conditions2=( + # Connect left boundary with right boundary of left mesh +# x_neg=BoundaryConditionCoupled(1, (:end, :i_forward), Float64), + x_neg=boundary_condition_periodic, + x_pos=boundary_condition_periodic, + y_neg=boundary_condition_periodic, + y_pos=boundary_condition_periodic) # A semidiscretization collects data structures and functions for the spatial discretization -semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition_convergence_test, solver) +# semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition_convergence_test, solver) +semi1 = SemidiscretizationHyperbolic(mesh1, equations, initial_condition_convergence_test, solver, + boundary_conditions=boundary_conditions1) +semi2 = SemidiscretizationHyperbolic(mesh2, equations, initial_condition_convergence_test, solver, + boundary_conditions=boundary_conditions2) +semi = SemidiscretizationCoupled(semi1, semi2) ############################################################################### @@ -40,7 +64,7 @@ ode = semidiscretize(semi, (0.0, 1.0)); summary_callback = SummaryCallback() # The AnalysisCallback allows to analyse the solution in regular intervals and prints the results -analysis_callback = AnalysisCallback(semi, interval=10) +#analysis_callback = AnalysisCallback(semi, interval=10) # The SaveSolutionCallback allows to save the solution to a file in regular intervals # save_solution = SaveSolutionCallback(interval=100, @@ -51,7 +75,7 @@ analysis_callback = AnalysisCallback(semi, interval=10) # Create a CallbackSet to collect all callbacks such that they can be passed to the ODE solver # callbacks = CallbackSet(summary_callback, analysis_callback, save_solution, stepsize_callback) -callbacks = CallbackSet(summary_callback, analysis_callback) +callbacks = CallbackSet(summary_callback) ############################################################################### diff --git a/src/meshes/structured_mesh_view.jl b/src/meshes/structured_mesh_view.jl index eb2fb1d2419..26d9ccd05bf 100644 --- a/src/meshes/structured_mesh_view.jl +++ b/src/meshes/structured_mesh_view.jl @@ -1,8 +1,12 @@ +@muladd begin +#! format: noindent + mutable struct StructuredMeshView{NDIMS, RealT <: Real} <: AbstractMesh{NDIMS} parent::StructuredMesh{NDIMS, RealT} mapping::Any # Not relevant for performance index_min::NTuple{NDIMS, Int} index_max::NTuple{NDIMS, Int} + periodicity::NTuple{NDIMS, Bool} end function StructuredMeshView(parent::StructuredMesh{NDIMS, RealT}; @@ -11,7 +15,8 @@ function StructuredMeshView(parent::StructuredMesh{NDIMS, RealT}; @assert all(index_min .> 0) @assert index_max <= size(parent) - return StructuredMeshView{NDIMS, RealT}(parent, parent.mapping, index_min, index_max) + + return StructuredMeshView{NDIMS, RealT}(parent, parent.mapping, index_min, index_max, periodicity) end function StructuredMeshView(parent::StructuredMesh{NDIMS, RealT}) where {NDIMS, RealT} @@ -47,7 +52,8 @@ Base.axes(mesh::StructuredMeshView, i) = Base.OneTo(size(mesh, i)) function calc_node_coordinates!(node_coordinates, element, cell_x, cell_y, mapping, mesh::StructuredMeshView{2}, - basis::LobattoLegendreBasis) +# basis::LobattoLegendreBasis) + basis) @unpack nodes = basis @unpack parent, index_min, index_max = mesh @@ -69,3 +75,4 @@ function calc_node_coordinates!(node_coordinates, element, cell_y_offset + dy / 2 * nodes[j]) end end +end # @muladd diff --git a/src/semidiscretization/semidiscretization_coupled.jl b/src/semidiscretization/semidiscretization_coupled.jl index b7adff78425..a74bed59eb8 100644 --- a/src/semidiscretization/semidiscretization_coupled.jl +++ b/src/semidiscretization/semidiscretization_coupled.jl @@ -513,7 +513,7 @@ end @inline function calc_boundary_flux_by_direction!(surface_flux_values, u, t, orientation, boundary_condition::BoundaryConditionCoupled, - mesh::StructuredMesh, equations, + mesh::Union{StructuredMesh,StructuredMeshView}, equations, surface_integral, dg::DG, cache, direction, node_indices, surface_node_indices, element) @@ -544,7 +544,7 @@ end end end -function get_boundary_indices(element, orientation, mesh::StructuredMesh{2}) +function get_boundary_indices(element, orientation, mesh::Union{StructuredMesh{2},StructuredMeshView{2}}) cartesian_indices = CartesianIndices(size(mesh)) if orientation == 1 # Get index of element in y-direction diff --git a/src/solvers/dgsem_structured/dg.jl b/src/solvers/dgsem_structured/dg.jl index 571add15887..1255a4ac0d5 100644 --- a/src/solvers/dgsem_structured/dg.jl +++ b/src/solvers/dgsem_structured/dg.jl @@ -30,7 +30,7 @@ end @inline function calc_boundary_flux_by_direction!(surface_flux_values, u, t, orientation, boundary_condition::BoundaryConditionPeriodic, - mesh::StructuredMesh, equations, + mesh::Union{StructuredMesh,StructuredMeshView}, equations, surface_integral, dg::DG, cache, direction, node_indices, surface_node_indices, element) @@ -40,7 +40,7 @@ end @inline function calc_boundary_flux_by_direction!(surface_flux_values, u, t, orientation, boundary_condition, - mesh::StructuredMesh, equations, + mesh::Union{StructuredMesh,StructureMeshView}, equations, surface_integral, dg::DG, cache, direction, node_indices, surface_node_indices, element) diff --git a/src/solvers/dgsem_structured/dg_2d.jl b/src/solvers/dgsem_structured/dg_2d.jl index 3feb17f874d..f2fad34fc66 100644 --- a/src/solvers/dgsem_structured/dg_2d.jl +++ b/src/solvers/dgsem_structured/dg_2d.jl @@ -551,7 +551,7 @@ function calc_boundary_flux!(cache, u, t, boundary_condition::BoundaryConditionP end function calc_boundary_flux!(cache, u, t, boundary_conditions::NamedTuple, - mesh::StructuredMesh{2}, equations, surface_integral, + mesh::Union{StructuredMesh{2},StructuredMeshView{2}}, equations, surface_integral, dg::DG) @unpack surface_flux_values = cache.elements linear_indices = LinearIndices(size(mesh)) From f0684f920445a98ec2bf5f49b9004962687b142d Mon Sep 17 00:00:00 2001 From: Simon Candelaresi <10759273+SimonCan@users.noreply.github.com> Date: Fri, 20 Oct 2023 17:50:29 +0100 Subject: [PATCH 04/56] Update structured_mesh_view.jl Format. --- src/meshes/structured_mesh_view.jl | 1 - 1 file changed, 1 deletion(-) diff --git a/src/meshes/structured_mesh_view.jl b/src/meshes/structured_mesh_view.jl index 26d9ccd05bf..3d4788793a5 100644 --- a/src/meshes/structured_mesh_view.jl +++ b/src/meshes/structured_mesh_view.jl @@ -15,7 +15,6 @@ function StructuredMeshView(parent::StructuredMesh{NDIMS, RealT}; @assert all(index_min .> 0) @assert index_max <= size(parent) - return StructuredMeshView{NDIMS, RealT}(parent, parent.mapping, index_min, index_max, periodicity) end From c555b186c12560cc3e603efd2241461361134393 Mon Sep 17 00:00:00 2001 From: SimonCan Date: Thu, 9 Nov 2023 16:34:32 +0000 Subject: [PATCH 05/56] Applied autoformatter on smview elixir. --- .../elixir_advection_smview.jl | 52 +++++++++---------- 1 file changed, 26 insertions(+), 26 deletions(-) diff --git a/examples/structured_2d_dgsem/elixir_advection_smview.jl b/examples/structured_2d_dgsem/elixir_advection_smview.jl index 7d28299ef79..69d7ce2cd22 100644 --- a/examples/structured_2d_dgsem/elixir_advection_smview.jl +++ b/examples/structured_2d_dgsem/elixir_advection_smview.jl @@ -11,10 +11,10 @@ advection_velocity = (0.2, -0.7) equations = LinearScalarAdvectionEquation2D(advection_velocity) # Create DG solver with polynomial degree = 2 and (local) Lax-Friedrichs/Rusanov flux as surface flux -solver = DGSEM(polydeg=2, surface_flux=flux_lax_friedrichs) +solver = DGSEM(polydeg = 2, surface_flux = flux_lax_friedrichs) -coordinates_min = ( -1.0, -1.0) # maximum coordinates (max(x), max(y)) -coordinates_max = ( 1.0, 1.0) # maximum coordinates (max(x), max(y)) +coordinates_min = (-1.0, -1.0) # maximum coordinates (max(x), max(y)) +coordinates_max = (1.0, 1.0) # maximum coordinates (max(x), max(y)) cells_per_dimension = (16, 16) # cells_per_dimension = (32, 16) @@ -29,30 +29,31 @@ mesh1 = StructuredMeshView(parent; index_min = (1, 1), index_max = (8, 16)) mesh2 = StructuredMeshView(parent; index_min = (9, 1), index_max = (16, 16)) # Define the coupled boundary conditions -boundary_conditions1=( - # Connect left boundary with right boundary of left mesh - x_neg=boundary_condition_periodic, - x_pos=boundary_condition_periodic, -# x_pos=BoundaryConditionCoupled(2, (:begin, :i_forward), Float64), - y_neg=boundary_condition_periodic, - y_pos=boundary_condition_periodic) -boundary_conditions2=( - # Connect left boundary with right boundary of left mesh -# x_neg=BoundaryConditionCoupled(1, (:end, :i_forward), Float64), - x_neg=boundary_condition_periodic, - x_pos=boundary_condition_periodic, - y_neg=boundary_condition_periodic, - y_pos=boundary_condition_periodic) +boundary_conditions1 = ( + # Connect left boundary with right boundary of left mesh + x_neg = boundary_condition_periodic, + x_pos = boundary_condition_periodic, + # x_pos=BoundaryConditionCoupled(2, (:begin, :i_forward), Float64), + y_neg = boundary_condition_periodic, + y_pos = boundary_condition_periodic) +boundary_conditions2 = ( + # Connect left boundary with right boundary of left mesh + # x_neg=BoundaryConditionCoupled(1, (:end, :i_forward), Float64), + x_neg = boundary_condition_periodic, + x_pos = boundary_condition_periodic, + y_neg = boundary_condition_periodic, + y_pos = boundary_condition_periodic) # A semidiscretization collects data structures and functions for the spatial discretization # semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition_convergence_test, solver) -semi1 = SemidiscretizationHyperbolic(mesh1, equations, initial_condition_convergence_test, solver, - boundary_conditions=boundary_conditions1) -semi2 = SemidiscretizationHyperbolic(mesh2, equations, initial_condition_convergence_test, solver, - boundary_conditions=boundary_conditions2) +semi1 = SemidiscretizationHyperbolic(mesh1, equations, initial_condition_convergence_test, + solver, + boundary_conditions = boundary_conditions1) +semi2 = SemidiscretizationHyperbolic(mesh2, equations, initial_condition_convergence_test, + solver, + boundary_conditions = boundary_conditions2) semi = SemidiscretizationCoupled(semi1, semi2) - ############################################################################### # ODE solvers, callbacks etc. @@ -77,14 +78,13 @@ summary_callback = SummaryCallback() # callbacks = CallbackSet(summary_callback, analysis_callback, save_solution, stepsize_callback) callbacks = CallbackSet(summary_callback) - ############################################################################### # run the simulation # OrdinaryDiffEq's `solve` method evolves the solution in time and executes the passed callbacks -sol = solve(ode, CarpenterKennedy2N54(williamson_condition=false), - dt=5.0e-2, # solve needs some value here but it will be overwritten by the stepsize_callback - save_everystep=false, callback=callbacks); +sol = solve(ode, CarpenterKennedy2N54(williamson_condition = false), + dt = 5.0e-2, # solve needs some value here but it will be overwritten by the stepsize_callback + save_everystep = false, callback = callbacks); errors_ref = (l2 = [8.312427642603623e-6], linf = [6.626865824577166e-5]) From 11645a937b397b4d68dda7d7e754eb28e2e70568 Mon Sep 17 00:00:00 2001 From: SimonCan Date: Thu, 9 Nov 2023 16:35:15 +0000 Subject: [PATCH 06/56] Applied autoformatter on meshview. --- src/meshes/structured_mesh_view.jl | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/src/meshes/structured_mesh_view.jl b/src/meshes/structured_mesh_view.jl index 3d4788793a5..fae86af7113 100644 --- a/src/meshes/structured_mesh_view.jl +++ b/src/meshes/structured_mesh_view.jl @@ -15,7 +15,8 @@ function StructuredMeshView(parent::StructuredMesh{NDIMS, RealT}; @assert all(index_min .> 0) @assert index_max <= size(parent) - return StructuredMeshView{NDIMS, RealT}(parent, parent.mapping, index_min, index_max, periodicity) + return StructuredMeshView{NDIMS, RealT}(parent, parent.mapping, index_min, + index_max, periodicity) end function StructuredMeshView(parent::StructuredMesh{NDIMS, RealT}) where {NDIMS, RealT} @@ -47,11 +48,10 @@ end Base.axes(mesh::StructuredMeshView) = map(Base.OneTo, size(mesh)) Base.axes(mesh::StructuredMeshView, i) = Base.OneTo(size(mesh, i)) - function calc_node_coordinates!(node_coordinates, element, cell_x, cell_y, mapping, mesh::StructuredMeshView{2}, -# basis::LobattoLegendreBasis) + # basis::LobattoLegendreBasis) basis) @unpack nodes = basis @unpack parent, index_min, index_max = mesh From a9caff896f15417c6ed04d53beac446700df2181 Mon Sep 17 00:00:00 2001 From: SimonCan Date: Wed, 22 Nov 2023 13:41:51 +0000 Subject: [PATCH 07/56] Corected minor typo with StructuredMeshView. --- src/solvers/dgsem_structured/dg.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/solvers/dgsem_structured/dg.jl b/src/solvers/dgsem_structured/dg.jl index 1255a4ac0d5..fa19982c52c 100644 --- a/src/solvers/dgsem_structured/dg.jl +++ b/src/solvers/dgsem_structured/dg.jl @@ -40,7 +40,7 @@ end @inline function calc_boundary_flux_by_direction!(surface_flux_values, u, t, orientation, boundary_condition, - mesh::Union{StructuredMesh,StructureMeshView}, equations, + mesh::Union{StructuredMesh,StructuredMeshView}, equations, surface_integral, dg::DG, cache, direction, node_indices, surface_node_indices, element) From 1c7a833c4a48028e9cc390ddd7a87322987d9f44 Mon Sep 17 00:00:00 2001 From: SimonCan Date: Wed, 22 Nov 2023 14:39:42 +0000 Subject: [PATCH 08/56] Corrected x-boundaries for smview example elixir. --- .../structured_2d_dgsem/elixir_advection_smview.jl | 14 ++++++++------ 1 file changed, 8 insertions(+), 6 deletions(-) diff --git a/examples/structured_2d_dgsem/elixir_advection_smview.jl b/examples/structured_2d_dgsem/elixir_advection_smview.jl index 69d7ce2cd22..6fc57036a39 100644 --- a/examples/structured_2d_dgsem/elixir_advection_smview.jl +++ b/examples/structured_2d_dgsem/elixir_advection_smview.jl @@ -31,16 +31,18 @@ mesh2 = StructuredMeshView(parent; index_min = (9, 1), index_max = (16, 16)) # Define the coupled boundary conditions boundary_conditions1 = ( # Connect left boundary with right boundary of left mesh - x_neg = boundary_condition_periodic, - x_pos = boundary_condition_periodic, - # x_pos=BoundaryConditionCoupled(2, (:begin, :i_forward), Float64), + x_neg=BoundaryConditionCoupled(2, (:end, :i_forward), Float64), + x_pos=BoundaryConditionCoupled(2, (:begin, :i_forward), Float64), +# x_neg = boundary_condition_periodic, +# x_pos = boundary_condition_periodic, y_neg = boundary_condition_periodic, y_pos = boundary_condition_periodic) boundary_conditions2 = ( # Connect left boundary with right boundary of left mesh - # x_neg=BoundaryConditionCoupled(1, (:end, :i_forward), Float64), - x_neg = boundary_condition_periodic, - x_pos = boundary_condition_periodic, + x_neg=BoundaryConditionCoupled(1, (:end, :i_forward), Float64), + x_pos=BoundaryConditionCoupled(1, (:begin, :i_forward), Float64), +# x_neg = boundary_condition_periodic, +# x_pos = boundary_condition_periodic, y_neg = boundary_condition_periodic, y_pos = boundary_condition_periodic) From 607b623d98312f2eb01e98c58cb4f08271022a53 Mon Sep 17 00:00:00 2001 From: SimonCan Date: Wed, 22 Nov 2023 14:40:06 +0000 Subject: [PATCH 09/56] Removed parent's periodicity for smview. --- src/meshes/structured_mesh_view.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/meshes/structured_mesh_view.jl b/src/meshes/structured_mesh_view.jl index fae86af7113..5e4d9d84c32 100644 --- a/src/meshes/structured_mesh_view.jl +++ b/src/meshes/structured_mesh_view.jl @@ -6,7 +6,6 @@ mutable struct StructuredMeshView{NDIMS, RealT <: Real} <: AbstractMesh{NDIMS} mapping::Any # Not relevant for performance index_min::NTuple{NDIMS, Int} index_max::NTuple{NDIMS, Int} - periodicity::NTuple{NDIMS, Bool} end function StructuredMeshView(parent::StructuredMesh{NDIMS, RealT}; @@ -16,7 +15,7 @@ function StructuredMeshView(parent::StructuredMesh{NDIMS, RealT}; @assert index_max <= size(parent) return StructuredMeshView{NDIMS, RealT}(parent, parent.mapping, index_min, - index_max, periodicity) + index_max) end function StructuredMeshView(parent::StructuredMesh{NDIMS, RealT}) where {NDIMS, RealT} @@ -28,6 +27,7 @@ function isperiodic(mesh::StructuredMeshView) @unpack parent = mesh return isperiodic(parent) && size(parent) == size(mesh) end + function isperiodic(mesh::StructuredMeshView, dimension) @unpack parent, index_min, index_max = mesh return (isperiodic(parent, dimension) && From bc2303a0cbb39f154231d9f962929c5a66c81ea7 Mon Sep 17 00:00:00 2001 From: SimonCan Date: Wed, 22 Nov 2023 16:37:29 +0000 Subject: [PATCH 10/56] Removed redundant and problematic redifintion of StructuredMeshView function. --- src/meshes/structured_mesh_view.jl | 7 ++----- 1 file changed, 2 insertions(+), 5 deletions(-) diff --git a/src/meshes/structured_mesh_view.jl b/src/meshes/structured_mesh_view.jl index 5e4d9d84c32..451f5baebcb 100644 --- a/src/meshes/structured_mesh_view.jl +++ b/src/meshes/structured_mesh_view.jl @@ -9,7 +9,8 @@ mutable struct StructuredMeshView{NDIMS, RealT <: Real} <: AbstractMesh{NDIMS} end function StructuredMeshView(parent::StructuredMesh{NDIMS, RealT}; - index_min, index_max) where {NDIMS, RealT} + index_min=ntuple(_ -> 1, Val(NDIMS)), + index_max=size(parent)) where {NDIMS, RealT} @assert index_min <= index_max @assert all(index_min .> 0) @assert index_max <= size(parent) @@ -18,10 +19,6 @@ function StructuredMeshView(parent::StructuredMesh{NDIMS, RealT}; index_max) end -function StructuredMeshView(parent::StructuredMesh{NDIMS, RealT}) where {NDIMS, RealT} - return StructuredMeshView(parent, ntuple(_ -> 1, Val(NDIMS)), size(parent)) -end - # Check if mesh is periodic function isperiodic(mesh::StructuredMeshView) @unpack parent = mesh From d4624f76a51c8311769ab5144407a43969683cb8 Mon Sep 17 00:00:00 2001 From: SimonCan Date: Thu, 23 Nov 2023 13:53:30 +0000 Subject: [PATCH 11/56] Applied auto formatter on files affected by the structured mesh view. --- .../elixir_advection_smview.jl | 16 ++++++++-------- src/callbacks_step/analysis_dg2d.jl | 6 ++++-- src/meshes/structured_mesh_view.jl | 4 ++-- .../semidiscretization_coupled.jl | 7 +++++-- src/solvers/dgsem_structured/containers.jl | 3 ++- src/solvers/dgsem_structured/containers_2d.jl | 5 +++-- src/solvers/dgsem_structured/dg.jl | 13 +++++++++---- src/solvers/dgsem_structured/dg_2d.jl | 12 +++++++----- src/solvers/dgsem_tree/dg_2d.jl | 5 +++-- 9 files changed, 43 insertions(+), 28 deletions(-) diff --git a/examples/structured_2d_dgsem/elixir_advection_smview.jl b/examples/structured_2d_dgsem/elixir_advection_smview.jl index 6fc57036a39..2ec97cc9831 100644 --- a/examples/structured_2d_dgsem/elixir_advection_smview.jl +++ b/examples/structured_2d_dgsem/elixir_advection_smview.jl @@ -31,18 +31,18 @@ mesh2 = StructuredMeshView(parent; index_min = (9, 1), index_max = (16, 16)) # Define the coupled boundary conditions boundary_conditions1 = ( # Connect left boundary with right boundary of left mesh - x_neg=BoundaryConditionCoupled(2, (:end, :i_forward), Float64), - x_pos=BoundaryConditionCoupled(2, (:begin, :i_forward), Float64), -# x_neg = boundary_condition_periodic, -# x_pos = boundary_condition_periodic, + x_neg = BoundaryConditionCoupled(2, (:end, :i_forward), Float64), + x_pos = BoundaryConditionCoupled(2, (:begin, :i_forward), Float64), + # x_neg = boundary_condition_periodic, + # x_pos = boundary_condition_periodic, y_neg = boundary_condition_periodic, y_pos = boundary_condition_periodic) boundary_conditions2 = ( # Connect left boundary with right boundary of left mesh - x_neg=BoundaryConditionCoupled(1, (:end, :i_forward), Float64), - x_pos=BoundaryConditionCoupled(1, (:begin, :i_forward), Float64), -# x_neg = boundary_condition_periodic, -# x_pos = boundary_condition_periodic, + x_neg = BoundaryConditionCoupled(1, (:end, :i_forward), Float64), + x_pos = BoundaryConditionCoupled(1, (:begin, :i_forward), Float64), + # x_neg = boundary_condition_periodic, + # x_pos = boundary_condition_periodic, y_neg = boundary_condition_periodic, y_pos = boundary_condition_periodic) diff --git a/src/callbacks_step/analysis_dg2d.jl b/src/callbacks_step/analysis_dg2d.jl index 8f0b5fcde5a..de6b9a2a4a6 100644 --- a/src/callbacks_step/analysis_dg2d.jl +++ b/src/callbacks_step/analysis_dg2d.jl @@ -30,7 +30,8 @@ function create_cache_analysis(analyzer, mesh::TreeMesh{2}, end function create_cache_analysis(analyzer, - mesh::Union{StructuredMesh{2}, StructuredMeshView{2}, UnstructuredMesh2D, + mesh::Union{StructuredMesh{2}, StructuredMeshView{2}, + UnstructuredMesh2D, P4estMesh{2}, T8codeMesh{2}}, equations, dg::DG, cache, RealT, uEltype) @@ -177,7 +178,8 @@ end function integrate_via_indices(func::Func, u, mesh::Union{StructuredMesh{2}, StructuredMeshView{2}, - UnstructuredMesh2D, P4estMesh{2}, T8codeMesh{2}}, + UnstructuredMesh2D, P4estMesh{2}, + T8codeMesh{2}}, equations, dg::DGSEM, cache, args...; normalize = true) where {Func} @unpack weights = dg.basis diff --git a/src/meshes/structured_mesh_view.jl b/src/meshes/structured_mesh_view.jl index 451f5baebcb..a440a7d3357 100644 --- a/src/meshes/structured_mesh_view.jl +++ b/src/meshes/structured_mesh_view.jl @@ -9,8 +9,8 @@ mutable struct StructuredMeshView{NDIMS, RealT <: Real} <: AbstractMesh{NDIMS} end function StructuredMeshView(parent::StructuredMesh{NDIMS, RealT}; - index_min=ntuple(_ -> 1, Val(NDIMS)), - index_max=size(parent)) where {NDIMS, RealT} + index_min = ntuple(_ -> 1, Val(NDIMS)), + index_max = size(parent)) where {NDIMS, RealT} @assert index_min <= index_max @assert all(index_min .> 0) @assert index_max <= size(parent) diff --git a/src/semidiscretization/semidiscretization_coupled.jl b/src/semidiscretization/semidiscretization_coupled.jl index 1cd42debd3a..872be63b744 100644 --- a/src/semidiscretization/semidiscretization_coupled.jl +++ b/src/semidiscretization/semidiscretization_coupled.jl @@ -513,7 +513,9 @@ end @inline function calc_boundary_flux_by_direction!(surface_flux_values, u, t, orientation, boundary_condition::BoundaryConditionCoupled, - mesh::Union{StructuredMesh,StructuredMeshView}, equations, + mesh::Union{StructuredMesh, + StructuredMeshView}, + equations, surface_integral, dg::DG, cache, direction, node_indices, surface_node_indices, element) @@ -544,7 +546,8 @@ end end end -function get_boundary_indices(element, orientation, mesh::Union{StructuredMesh{2},StructuredMeshView{2}}) +function get_boundary_indices(element, orientation, + mesh::Union{StructuredMesh{2}, StructuredMeshView{2}}) cartesian_indices = CartesianIndices(size(mesh)) if orientation == 1 # Get index of element in y-direction diff --git a/src/solvers/dgsem_structured/containers.jl b/src/solvers/dgsem_structured/containers.jl index 14f49641b55..5e8fb60c22f 100644 --- a/src/solvers/dgsem_structured/containers.jl +++ b/src/solvers/dgsem_structured/containers.jl @@ -23,7 +23,8 @@ struct ElementContainer{NDIMS, RealT <: Real, uEltype <: Real, NDIMSP1, NDIMSP2, end # Create element container and initialize element data -function init_elements(mesh::Union{StructuredMesh{NDIMS, RealT},StructuredMeshView{NDIMS, RealT}}, +function init_elements(mesh::Union{StructuredMesh{NDIMS, RealT}, + StructuredMeshView{NDIMS, RealT}}, equations::AbstractEquations, basis, ::Type{uEltype}) where {NDIMS, RealT <: Real, uEltype <: Real} diff --git a/src/solvers/dgsem_structured/containers_2d.jl b/src/solvers/dgsem_structured/containers_2d.jl index 9e77a9ce680..8a0722fc5d5 100644 --- a/src/solvers/dgsem_structured/containers_2d.jl +++ b/src/solvers/dgsem_structured/containers_2d.jl @@ -6,7 +6,7 @@ #! format: noindent # Initialize data structures in element container -function init_elements!(elements, mesh::Union{StructuredMesh{2},StructuredMeshView{2}}, +function init_elements!(elements, mesh::Union{StructuredMesh{2}, StructuredMeshView{2}}, basis::LobattoLegendreBasis) @unpack node_coordinates, left_neighbors, jacobian_matrix, contravariant_vectors, inverse_jacobian = elements @@ -150,7 +150,8 @@ end # Save id of left neighbor of every element function initialize_left_neighbor_connectivity!(left_neighbors, - mesh::Union{StructuredMesh{2},StructuredMeshView{2}}, + mesh::Union{StructuredMesh{2}, + StructuredMeshView{2}}, linear_indices) # Neighbors in x-direction for cell_y in 1:size(mesh, 2) diff --git a/src/solvers/dgsem_structured/dg.jl b/src/solvers/dgsem_structured/dg.jl index fa19982c52c..00e321fba65 100644 --- a/src/solvers/dgsem_structured/dg.jl +++ b/src/solvers/dgsem_structured/dg.jl @@ -8,8 +8,9 @@ # This method is called when a SemidiscretizationHyperbolic is constructed. # It constructs the basic `cache` used throughout the simulation to compute # the RHS etc. -function create_cache(mesh::Union{StructuredMesh,StructuredMeshView}, - equations::AbstractEquations, dg::DG, ::Any, ::Type{uEltype}) where {uEltype <: Real} +function create_cache(mesh::Union{StructuredMesh, StructuredMeshView}, + equations::AbstractEquations, dg::DG, ::Any, + ::Type{uEltype}) where {uEltype <: Real} elements = init_elements(mesh, equations, dg.basis, uEltype) cache = (; elements) @@ -30,7 +31,9 @@ end @inline function calc_boundary_flux_by_direction!(surface_flux_values, u, t, orientation, boundary_condition::BoundaryConditionPeriodic, - mesh::Union{StructuredMesh,StructuredMeshView}, equations, + mesh::Union{StructuredMesh, + StructuredMeshView}, + equations, surface_integral, dg::DG, cache, direction, node_indices, surface_node_indices, element) @@ -40,7 +43,9 @@ end @inline function calc_boundary_flux_by_direction!(surface_flux_values, u, t, orientation, boundary_condition, - mesh::Union{StructuredMesh,StructuredMeshView}, equations, + mesh::Union{StructuredMesh, + StructuredMeshView}, + equations, surface_integral, dg::DG, cache, direction, node_indices, surface_node_indices, element) diff --git a/src/solvers/dgsem_structured/dg_2d.jl b/src/solvers/dgsem_structured/dg_2d.jl index f8437ca337d..1eb244e06aa 100644 --- a/src/solvers/dgsem_structured/dg_2d.jl +++ b/src/solvers/dgsem_structured/dg_2d.jl @@ -6,7 +6,7 @@ #! format: noindent function rhs!(du, u, t, - mesh::Union{StructuredMesh{2},StructuredMeshView{2}}, equations, + mesh::Union{StructuredMesh{2}, StructuredMeshView{2}}, equations, initial_condition, boundary_conditions, source_terms::Source, dg::DG, cache) where {Source} # Reset du @@ -388,7 +388,7 @@ end end function calc_interface_flux!(cache, u, - mesh::Union{StructuredMesh{2},StructuredMeshView{2}}, + mesh::Union{StructuredMesh{2}, StructuredMeshView{2}}, nonconservative_terms, # can be True/False equations, surface_integral, dg::DG) @unpack elements = cache @@ -417,7 +417,8 @@ end @inline function calc_interface_flux!(surface_flux_values, left_element, right_element, orientation, u, - mesh::Union{StructuredMesh{2},StructuredMeshView{2}}, + mesh::Union{StructuredMesh{2}, + StructuredMeshView{2}}, nonconservative_terms::False, equations, surface_integral, dg::DG, cache) # This is slow for LSA, but for some reason faster for Euler (see #519) @@ -552,13 +553,14 @@ end # TODO: Taal dimension agnostic function calc_boundary_flux!(cache, u, t, boundary_condition::BoundaryConditionPeriodic, - mesh::Union{StructuredMesh{2},StructuredMeshView{2}}, + mesh::Union{StructuredMesh{2}, StructuredMeshView{2}}, equations, surface_integral, dg::DG) @assert isperiodic(mesh) end function calc_boundary_flux!(cache, u, t, boundary_conditions::NamedTuple, - mesh::Union{StructuredMesh{2},StructuredMeshView{2}}, equations, surface_integral, + mesh::Union{StructuredMesh{2}, StructuredMeshView{2}}, + equations, surface_integral, dg::DG) @unpack surface_flux_values = cache.elements linear_indices = LinearIndices(size(mesh)) diff --git a/src/solvers/dgsem_tree/dg_2d.jl b/src/solvers/dgsem_tree/dg_2d.jl index 5e2acdffb30..a9c7b961ba3 100644 --- a/src/solvers/dgsem_tree/dg_2d.jl +++ b/src/solvers/dgsem_tree/dg_2d.jl @@ -1085,8 +1085,9 @@ end return nothing end -function calc_surface_integral!(du, u, mesh::Union{TreeMesh{2}, StructuredMesh{2}, - StructuredMeshView{2}}, +function calc_surface_integral!(du, u, + mesh::Union{TreeMesh{2}, StructuredMesh{2}, + StructuredMeshView{2}}, equations, surface_integral::SurfaceIntegralWeakForm, dg::DG, cache) @unpack boundary_interpolation = dg.basis From 6f4d21e36777f1bdcf2a6f22d10ebe1b98dd2f8f Mon Sep 17 00:00:00 2001 From: SimonCan Date: Fri, 24 Nov 2023 14:58:03 +0000 Subject: [PATCH 12/56] Applied autoforatter. --- src/solvers/dg.jl | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/src/solvers/dg.jl b/src/solvers/dg.jl index dd7e4df96de..1fa8bd77fe2 100644 --- a/src/solvers/dg.jl +++ b/src/solvers/dg.jl @@ -435,7 +435,8 @@ function get_node_variables!(node_variables, mesh, equations, dg::DG, cache) get_node_variables!(node_variables, mesh, equations, dg.volume_integral, dg, cache) end -const MeshesDGSEM = Union{TreeMesh, StructuredMesh, StructuredMeshView, UnstructuredMesh2D, +const MeshesDGSEM = Union{TreeMesh, StructuredMesh, StructuredMeshView, + UnstructuredMesh2D, P4estMesh, T8codeMesh} @inline function ndofs(mesh::MeshesDGSEM, dg::DG, cache) From d6d0f1195376a6b7e891b765c4f56ae9b41e1462 Mon Sep 17 00:00:00 2001 From: SimonCan Date: Wed, 6 Dec 2023 13:03:05 +0000 Subject: [PATCH 13/56] Added entries for meshview IO. --- examples/structured_2d_dgsem/elixir_advection_smview.jl | 3 ++- src/meshes/mesh_io.jl | 2 +- src/meshes/structured_mesh_view.jl | 3 ++- 3 files changed, 5 insertions(+), 3 deletions(-) diff --git a/examples/structured_2d_dgsem/elixir_advection_smview.jl b/examples/structured_2d_dgsem/elixir_advection_smview.jl index 2ec97cc9831..c285e0adcb3 100644 --- a/examples/structured_2d_dgsem/elixir_advection_smview.jl +++ b/examples/structured_2d_dgsem/elixir_advection_smview.jl @@ -78,7 +78,8 @@ summary_callback = SummaryCallback() # Create a CallbackSet to collect all callbacks such that they can be passed to the ODE solver # callbacks = CallbackSet(summary_callback, analysis_callback, save_solution, stepsize_callback) -callbacks = CallbackSet(summary_callback) +callbacks = CallbackSet(summary_callback, + smview_callback) ############################################################################### # run the simulation diff --git a/src/meshes/mesh_io.jl b/src/meshes/mesh_io.jl index 92e38ce1bf3..66b9991f0b4 100644 --- a/src/meshes/mesh_io.jl +++ b/src/meshes/mesh_io.jl @@ -95,7 +95,7 @@ end # of the mesh, like its size and the type of boundary mapping function. # Then, within Trixi2Vtk, the StructuredMesh and its node coordinates are reconstructured from # these attributes for plotting purposes -function save_mesh_file(mesh::StructuredMesh, output_directory; system = "") +function save_mesh_file(mesh::Union{StructuredMesh, StructuredMeshView}, output_directory; system = "") # Create output directory (if it does not exist) mkpath(output_directory) diff --git a/src/meshes/structured_mesh_view.jl b/src/meshes/structured_mesh_view.jl index a440a7d3357..a81e1855991 100644 --- a/src/meshes/structured_mesh_view.jl +++ b/src/meshes/structured_mesh_view.jl @@ -6,6 +6,7 @@ mutable struct StructuredMeshView{NDIMS, RealT <: Real} <: AbstractMesh{NDIMS} mapping::Any # Not relevant for performance index_min::NTuple{NDIMS, Int} index_max::NTuple{NDIMS, Int} + unsaved_changes::Bool end function StructuredMeshView(parent::StructuredMesh{NDIMS, RealT}; @@ -16,7 +17,7 @@ function StructuredMeshView(parent::StructuredMesh{NDIMS, RealT}; @assert index_max <= size(parent) return StructuredMeshView{NDIMS, RealT}(parent, parent.mapping, index_min, - index_max) + index_max, parent.unsaved_changes) end # Check if mesh is periodic From a64bbc2ce7a39d4dfb1ea0cd0ad9a2cfd0b5a8d5 Mon Sep 17 00:00:00 2001 From: SimonCan Date: Thu, 7 Dec 2023 22:45:02 +0000 Subject: [PATCH 14/56] Added structuresd mesh view to parametric types. --- src/callbacks_step/save_solution_dg.jl | 2 +- src/meshes/mesh_io.jl | 2 +- src/meshes/structured_mesh_view.jl | 7 +++++-- 3 files changed, 7 insertions(+), 4 deletions(-) diff --git a/src/callbacks_step/save_solution_dg.jl b/src/callbacks_step/save_solution_dg.jl index 350aee7336a..26effed85a4 100644 --- a/src/callbacks_step/save_solution_dg.jl +++ b/src/callbacks_step/save_solution_dg.jl @@ -6,7 +6,7 @@ #! format: noindent function save_solution_file(u, time, dt, timestep, - mesh::Union{SerialTreeMesh, StructuredMesh, + mesh::Union{SerialTreeMesh, StructuredMesh, StructuredMeshView, UnstructuredMesh2D, SerialP4estMesh, SerialT8codeMesh}, equations, dg::DG, cache, diff --git a/src/meshes/mesh_io.jl b/src/meshes/mesh_io.jl index 1a6ae472689..9bfb759039c 100644 --- a/src/meshes/mesh_io.jl +++ b/src/meshes/mesh_io.jl @@ -255,7 +255,7 @@ function load_mesh_serial(mesh_file::AbstractString; n_cells_max, RealT) end mesh = TreeMesh(SerialTree{ndims}, max(n_cells_max, capacity)) load_mesh!(mesh, mesh_file) - elseif mesh_type == "StructuredMesh" + elseif (mesh_type == "StructuredMesh") || (mesh_type == "StructuredMeshView") size_, mapping_as_string = h5open(mesh_file, "r") do file return read(attributes(file)["size"]), read(attributes(file)["mapping"]) diff --git a/src/meshes/structured_mesh_view.jl b/src/meshes/structured_mesh_view.jl index a81e1855991..da8991892c6 100644 --- a/src/meshes/structured_mesh_view.jl +++ b/src/meshes/structured_mesh_view.jl @@ -4,6 +4,8 @@ mutable struct StructuredMeshView{NDIMS, RealT <: Real} <: AbstractMesh{NDIMS} parent::StructuredMesh{NDIMS, RealT} mapping::Any # Not relevant for performance + mapping_as_string::String + current_filename::String index_min::NTuple{NDIMS, Int} index_max::NTuple{NDIMS, Int} unsaved_changes::Bool @@ -16,8 +18,9 @@ function StructuredMeshView(parent::StructuredMesh{NDIMS, RealT}; @assert all(index_min .> 0) @assert index_max <= size(parent) - return StructuredMeshView{NDIMS, RealT}(parent, parent.mapping, index_min, - index_max, parent.unsaved_changes) + return StructuredMeshView{NDIMS, RealT}(parent, parent.mapping, parent.mapping_as_string, + parent.current_filename, + index_min, index_max, parent.unsaved_changes) end # Check if mesh is periodic From 076fd24a1d51750c8df36f93c7156bd97a132f7c Mon Sep 17 00:00:00 2001 From: SimonCan Date: Thu, 14 Dec 2023 14:42:13 +0000 Subject: [PATCH 15/56] Added unctionality of writing out mesh files for StructuredMeshView for different time steps. --- src/meshes/mesh_io.jl | 25 ++++++++++++++++++- .../semidiscretization_coupled.jl | 2 +- 2 files changed, 25 insertions(+), 2 deletions(-) diff --git a/src/meshes/mesh_io.jl b/src/meshes/mesh_io.jl index 9bfb759039c..3ed419ef433 100644 --- a/src/meshes/mesh_io.jl +++ b/src/meshes/mesh_io.jl @@ -96,7 +96,7 @@ end # of the mesh, like its size and the type of boundary mapping function. # Then, within Trixi2Vtk, the StructuredMesh and its node coordinates are reconstructured from # these attributes for plotting purposes -function save_mesh_file(mesh::Union{StructuredMesh, StructuredMeshView}, output_directory; system = "") +function save_mesh_file(mesh::StructuredMesh, output_directory; system = "", timestep = 0) # Create output directory (if it does not exist) mkpath(output_directory) @@ -118,6 +118,29 @@ function save_mesh_file(mesh::Union{StructuredMesh, StructuredMeshView}, output_ return filename end +# Does not save the mesh itself to an HDF5 file. Instead saves important attributes +# of the mesh, like its size and the type of boundary mapping function. +# Then, within Trixi2Vtk, the StructuredMesh and its node coordinates are reconstructured from +# these attributes for plotting purposes +function save_mesh_file(mesh::StructuredMeshView, output_directory; system = "", timestep = 0) + # Create output directory (if it does not exist) + mkpath(output_directory) + + # filename = joinpath(output_directory, @sprintf("mesh_%s.h5", system)) + filename = joinpath(output_directory, @sprintf("mesh_%s_%06d.h5", system, timestep)) + + # Open file (clobber existing content) + h5open(filename, "w") do file + # Add context information as attributes + attributes(file)["mesh_type"] = get_name(mesh) + attributes(file)["ndims"] = ndims(mesh) + attributes(file)["size"] = collect(size(mesh)) + attributes(file)["mapping"] = mesh.mapping_as_string + end + + return filename +end + # Does not save the mesh itself to an HDF5 file. Instead saves important attributes # of the mesh, like its size and the corresponding `.mesh` file used to construct the mesh. # Then, within Trixi2Vtk, the UnstructuredMesh2D and its node coordinates are reconstructured diff --git a/src/semidiscretization/semidiscretization_coupled.jl b/src/semidiscretization/semidiscretization_coupled.jl index 3b85c63cdc7..a2390f9bad8 100644 --- a/src/semidiscretization/semidiscretization_coupled.jl +++ b/src/semidiscretization/semidiscretization_coupled.jl @@ -296,7 +296,7 @@ function save_mesh(semi::SemidiscretizationCoupled, output_directory, timestep = mesh, _, _, _ = mesh_equations_solver_cache(semi.semis[i]) if mesh.unsaved_changes - mesh.current_filename = save_mesh_file(mesh, output_directory, system = i) + mesh.current_filename = save_mesh_file(mesh, output_directory; system = i, timestep = timestep) mesh.unsaved_changes = false end end From 65d0413df0978fa6bf478991316bf08647204563 Mon Sep 17 00:00:00 2001 From: SimonCan Date: Fri, 5 Jan 2024 13:23:16 +0000 Subject: [PATCH 16/56] Added temptative code for dynamically changing mesh view sizes. --- src/meshes/structured_mesh_view.jl | 14 ++++++++++++++ 1 file changed, 14 insertions(+) diff --git a/src/meshes/structured_mesh_view.jl b/src/meshes/structured_mesh_view.jl index da8991892c6..4d1fbd85ce0 100644 --- a/src/meshes/structured_mesh_view.jl +++ b/src/meshes/structured_mesh_view.jl @@ -18,6 +18,20 @@ function StructuredMeshView(parent::StructuredMesh{NDIMS, RealT}; @assert all(index_min .> 0) @assert index_max <= size(parent) +# # Calculate the domain boundaries. +# coordinates_min = +# coordinates_max = +# mapping = coordinates2mapping(coordinates_min, coordinates_max) +# mapping_as_string = """ +# coordinates_min = $coordinates_min +# coordinates_max = $coordinates_max +# mapping = coordinates2mapping(coordinates_min, coordinates_max) +# """ +# +# return StructuredMeshView{NDIMS, RealT}(parent, mapping, mapping_as_string, +# parent.current_filename, +# index_min, index_max, parent.unsaved_changes) +# return StructuredMeshView{NDIMS, RealT}(parent, parent.mapping, parent.mapping_as_string, parent.current_filename, index_min, index_max, parent.unsaved_changes) From 798ea2247d740fc08e0be324467b2cd8e32fdcbd Mon Sep 17 00:00:00 2001 From: SimonCan Date: Thu, 11 Jan 2024 15:41:00 +0000 Subject: [PATCH 17/56] Applied autoformatter. --- src/callbacks_step/save_solution_dg.jl | 9 +++-- src/meshes/mesh_io.jl | 8 +++-- src/meshes/structured_mesh_view.jl | 34 ++++++++++--------- .../semidiscretization_coupled.jl | 6 ++-- 4 files changed, 33 insertions(+), 24 deletions(-) diff --git a/src/callbacks_step/save_solution_dg.jl b/src/callbacks_step/save_solution_dg.jl index 26effed85a4..424821e6b33 100644 --- a/src/callbacks_step/save_solution_dg.jl +++ b/src/callbacks_step/save_solution_dg.jl @@ -6,7 +6,8 @@ #! format: noindent function save_solution_file(u, time, dt, timestep, - mesh::Union{SerialTreeMesh, StructuredMesh, StructuredMeshView, + mesh::Union{SerialTreeMesh, StructuredMesh, + StructuredMeshView, UnstructuredMesh2D, SerialP4estMesh, SerialT8codeMesh}, equations, dg::DG, cache, @@ -33,7 +34,8 @@ function save_solution_file(u, time, dt, timestep, # compute the solution variables via broadcasting, and reinterpret the # result as a plain array of floating point numbers data = Array(reinterpret(eltype(u), - solution_variables.(reinterpret(SVector{nvariables(equations), + solution_variables.(reinterpret(SVector{ + nvariables(equations), eltype(u)}, u), Ref(equations)))) @@ -115,7 +117,8 @@ function save_solution_file(u, time, dt, timestep, # compute the solution variables via broadcasting, and reinterpret the # result as a plain array of floating point numbers data = Array(reinterpret(eltype(u), - solution_variables.(reinterpret(SVector{nvariables(equations), + solution_variables.(reinterpret(SVector{ + nvariables(equations), eltype(u)}, u), Ref(equations)))) diff --git a/src/meshes/mesh_io.jl b/src/meshes/mesh_io.jl index 3ed419ef433..f10a666e289 100644 --- a/src/meshes/mesh_io.jl +++ b/src/meshes/mesh_io.jl @@ -96,7 +96,8 @@ end # of the mesh, like its size and the type of boundary mapping function. # Then, within Trixi2Vtk, the StructuredMesh and its node coordinates are reconstructured from # these attributes for plotting purposes -function save_mesh_file(mesh::StructuredMesh, output_directory; system = "", timestep = 0) +function save_mesh_file(mesh::StructuredMesh, output_directory; system = "", + timestep = 0) # Create output directory (if it does not exist) mkpath(output_directory) @@ -122,12 +123,13 @@ end # of the mesh, like its size and the type of boundary mapping function. # Then, within Trixi2Vtk, the StructuredMesh and its node coordinates are reconstructured from # these attributes for plotting purposes -function save_mesh_file(mesh::StructuredMeshView, output_directory; system = "", timestep = 0) +function save_mesh_file(mesh::StructuredMeshView, output_directory; system = "", + timestep = 0) # Create output directory (if it does not exist) mkpath(output_directory) # filename = joinpath(output_directory, @sprintf("mesh_%s.h5", system)) - filename = joinpath(output_directory, @sprintf("mesh_%s_%06d.h5", system, timestep)) + filename = joinpath(output_directory, @sprintf("mesh_%s_%06d.h5", system, timestep)) # Open file (clobber existing content) h5open(filename, "w") do file diff --git a/src/meshes/structured_mesh_view.jl b/src/meshes/structured_mesh_view.jl index 4d1fbd85ce0..37ba13531d1 100644 --- a/src/meshes/structured_mesh_view.jl +++ b/src/meshes/structured_mesh_view.jl @@ -18,23 +18,25 @@ function StructuredMeshView(parent::StructuredMesh{NDIMS, RealT}; @assert all(index_min .> 0) @assert index_max <= size(parent) -# # Calculate the domain boundaries. -# coordinates_min = -# coordinates_max = -# mapping = coordinates2mapping(coordinates_min, coordinates_max) -# mapping_as_string = """ -# coordinates_min = $coordinates_min -# coordinates_max = $coordinates_max -# mapping = coordinates2mapping(coordinates_min, coordinates_max) -# """ -# -# return StructuredMeshView{NDIMS, RealT}(parent, mapping, mapping_as_string, -# parent.current_filename, -# index_min, index_max, parent.unsaved_changes) -# - return StructuredMeshView{NDIMS, RealT}(parent, parent.mapping, parent.mapping_as_string, + # # Calculate the domain boundaries. + # coordinates_min = + # coordinates_max = + # mapping = coordinates2mapping(coordinates_min, coordinates_max) + # mapping_as_string = """ + # coordinates_min = $coordinates_min + # coordinates_max = $coordinates_max + # mapping = coordinates2mapping(coordinates_min, coordinates_max) + # """ + # + # return StructuredMeshView{NDIMS, RealT}(parent, mapping, mapping_as_string, + # parent.current_filename, + # index_min, index_max, parent.unsaved_changes) + # + return StructuredMeshView{NDIMS, RealT}(parent, parent.mapping, + parent.mapping_as_string, parent.current_filename, - index_min, index_max, parent.unsaved_changes) + index_min, index_max, + parent.unsaved_changes) end # Check if mesh is periodic diff --git a/src/semidiscretization/semidiscretization_coupled.jl b/src/semidiscretization/semidiscretization_coupled.jl index a2390f9bad8..386eea4ca75 100644 --- a/src/semidiscretization/semidiscretization_coupled.jl +++ b/src/semidiscretization/semidiscretization_coupled.jl @@ -296,7 +296,8 @@ function save_mesh(semi::SemidiscretizationCoupled, output_directory, timestep = mesh, _, _, _ = mesh_equations_solver_cache(semi.semis[i]) if mesh.unsaved_changes - mesh.current_filename = save_mesh_file(mesh, output_directory; system = i, timestep = timestep) + mesh.current_filename = save_mesh_file(mesh, output_directory; system = i, + timestep = timestep) mesh.unsaved_changes = false end end @@ -433,7 +434,8 @@ function allocate_coupled_boundary_condition(boundary_condition, direction, mesh end # In 2D -function allocate_coupled_boundary_condition(boundary_condition::BoundaryConditionCoupled{2}, +function allocate_coupled_boundary_condition(boundary_condition::BoundaryConditionCoupled{2 + }, direction, mesh, equations, dg::DGSEM) if direction in (1, 2) cell_size = size(mesh, 2) From 5b9bcc307698ff3181bf869021c4941e7eba674e Mon Sep 17 00:00:00 2001 From: iomsn Date: Thu, 29 Feb 2024 17:52:11 +0000 Subject: [PATCH 18/56] Corrected the calculation of the coordinate mapping for mesh views. --- src/meshes/structured_mesh_view.jl | 37 +++++++++++++++--------------- 1 file changed, 19 insertions(+), 18 deletions(-) diff --git a/src/meshes/structured_mesh_view.jl b/src/meshes/structured_mesh_view.jl index 37ba13531d1..4d998a00d12 100644 --- a/src/meshes/structured_mesh_view.jl +++ b/src/meshes/structured_mesh_view.jl @@ -18,25 +18,26 @@ function StructuredMeshView(parent::StructuredMesh{NDIMS, RealT}; @assert all(index_min .> 0) @assert index_max <= size(parent) - # # Calculate the domain boundaries. - # coordinates_min = - # coordinates_max = - # mapping = coordinates2mapping(coordinates_min, coordinates_max) - # mapping_as_string = """ - # coordinates_min = $coordinates_min - # coordinates_max = $coordinates_max - # mapping = coordinates2mapping(coordinates_min, coordinates_max) - # """ - # - # return StructuredMeshView{NDIMS, RealT}(parent, mapping, mapping_as_string, - # parent.current_filename, - # index_min, index_max, parent.unsaved_changes) - # - return StructuredMeshView{NDIMS, RealT}(parent, parent.mapping, - parent.mapping_as_string, + # Calculate the domain boundaries. + deltas = (parent.mapping.coordinates_max .- parent.mapping.coordinates_min)./parent.cells_per_dimension + coordinates_min = parent.mapping.coordinates_min .+ deltas .* (index_min .- 1) + coordinates_max = parent.mapping.coordinates_min .+ deltas .* index_max + mapping = coordinates2mapping(coordinates_min, coordinates_max) + mapping_as_string = """ + coordinates_min = $coordinates_min + coordinates_max = $coordinates_max + mapping = coordinates2mapping(coordinates_min, coordinates_max) + """ + + return StructuredMeshView{NDIMS, RealT}(parent, mapping, mapping_as_string, parent.current_filename, - index_min, index_max, - parent.unsaved_changes) + index_min, index_max, parent.unsaved_changes) + + # return StructuredMeshView{NDIMS, RealT}(parent, parent.mapping, + # parent.mapping_as_string, + # parent.current_filename, + # index_min, index_max, + # parent.unsaved_changes) end # Check if mesh is periodic From f0ab06d0e40c7ffcb858dc14db1b3af72fb5f1ff Mon Sep 17 00:00:00 2001 From: iomsn Date: Fri, 1 Mar 2024 14:33:44 +0000 Subject: [PATCH 19/56] Added cells_per_dimension to the StructuredMeshView. --- src/meshes/structured_mesh_view.jl | 23 ++++++++++------------- 1 file changed, 10 insertions(+), 13 deletions(-) diff --git a/src/meshes/structured_mesh_view.jl b/src/meshes/structured_mesh_view.jl index 4d998a00d12..fd51b6a9dc8 100644 --- a/src/meshes/structured_mesh_view.jl +++ b/src/meshes/structured_mesh_view.jl @@ -3,6 +3,7 @@ mutable struct StructuredMeshView{NDIMS, RealT <: Real} <: AbstractMesh{NDIMS} parent::StructuredMesh{NDIMS, RealT} + cells_per_dimension::NTuple{NDIMS, Int} mapping::Any # Not relevant for performance mapping_as_string::String current_filename::String @@ -19,6 +20,7 @@ function StructuredMeshView(parent::StructuredMesh{NDIMS, RealT}; @assert index_max <= size(parent) # Calculate the domain boundaries. + cells_per_dimension = index_max .- index_min .+ 1 deltas = (parent.mapping.coordinates_max .- parent.mapping.coordinates_min)./parent.cells_per_dimension coordinates_min = parent.mapping.coordinates_min .+ deltas .* (index_min .- 1) coordinates_max = parent.mapping.coordinates_min .+ deltas .* index_max @@ -29,7 +31,7 @@ function StructuredMeshView(parent::StructuredMesh{NDIMS, RealT}; mapping = coordinates2mapping(coordinates_min, coordinates_max) """ - return StructuredMeshView{NDIMS, RealT}(parent, mapping, mapping_as_string, + return StructuredMeshView{NDIMS, RealT}(parent, cells_per_dimension, mapping, mapping_as_string, parent.current_filename, index_min, index_max, parent.unsaved_changes) @@ -69,27 +71,22 @@ Base.axes(mesh::StructuredMeshView, i) = Base.OneTo(size(mesh, i)) function calc_node_coordinates!(node_coordinates, element, cell_x, cell_y, mapping, mesh::StructuredMeshView{2}, - # basis::LobattoLegendreBasis) + # basis::LobattoLegendreBasis) basis) @unpack nodes = basis - @unpack parent, index_min, index_max = mesh # Get cell length in reference mesh - dx = 2 / size(parent, 1) - dy = 2 / size(parent, 2) - - # Calculate index offsets with respect to parent - parent_offset_x = index_min[1] - 1 - parent_offset_y = index_min[2] - 1 + dx = 2 / size(mesh, 1) + dy = 2 / size(mesh, 2) # Calculate node coordinates of reference mesh - cell_x_offset = -1 + (cell_x - 1 + parent_offset_x) * dx + dx / 2 - cell_y_offset = -1 + (cell_y - 1 + parent_offset_y) * dy + dy / 2 + cell_x_offset = -1 + (cell_x - 1) * dx + dx / 2 + cell_y_offset = -1 + (cell_y - 1) * dy + dy / 2 for j in eachnode(basis), i in eachnode(basis) # node_coordinates are the mapped reference node_coordinates node_coordinates[:, i, j, element] .= mapping(cell_x_offset + dx / 2 * nodes[i], - cell_y_offset + dy / 2 * nodes[j]) - end + cell_y_offset + dy / 2 * nodes[j]) + end end end # @muladd From 82b2014bff5bc4a740d24ed39569833a15bd1418 Mon Sep 17 00:00:00 2001 From: iomsn Date: Thu, 7 Mar 2024 13:41:44 +0000 Subject: [PATCH 20/56] Added StructuredMeshView to max_dt clculations. --- src/callbacks_step/stepsize_dg2d.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/callbacks_step/stepsize_dg2d.jl b/src/callbacks_step/stepsize_dg2d.jl index c6d32c0f6dc..41251506a0d 100644 --- a/src/callbacks_step/stepsize_dg2d.jl +++ b/src/callbacks_step/stepsize_dg2d.jl @@ -77,7 +77,7 @@ end function max_dt(u, t, mesh::Union{StructuredMesh{2}, UnstructuredMesh2D, P4estMesh{2}, - T8codeMesh{2}}, + T8codeMesh{2}, StructuredMeshView{2}}, constant_speed::False, equations, dg::DG, cache) # to avoid a division by zero if the speed vanishes everywhere, # e.g. for steady-state linear advection @@ -113,7 +113,7 @@ end function max_dt(u, t, mesh::Union{StructuredMesh{2}, UnstructuredMesh2D, P4estMesh{2}, - T8codeMesh{2}}, + T8codeMesh{2}, StructuredMeshView{2}}, constant_speed::True, equations, dg::DG, cache) @unpack contravariant_vectors, inverse_jacobian = cache.elements From 51eb64abd5cabc7918737c0e3d49b5de86b70d9e Mon Sep 17 00:00:00 2001 From: iomsn Date: Mon, 11 Mar 2024 14:33:23 +0000 Subject: [PATCH 21/56] Applied autoformatter. --- src/meshes/structured_mesh_view.jl | 13 ++++++++----- 1 file changed, 8 insertions(+), 5 deletions(-) diff --git a/src/meshes/structured_mesh_view.jl b/src/meshes/structured_mesh_view.jl index fd51b6a9dc8..bf143a2b169 100644 --- a/src/meshes/structured_mesh_view.jl +++ b/src/meshes/structured_mesh_view.jl @@ -21,7 +21,8 @@ function StructuredMeshView(parent::StructuredMesh{NDIMS, RealT}; # Calculate the domain boundaries. cells_per_dimension = index_max .- index_min .+ 1 - deltas = (parent.mapping.coordinates_max .- parent.mapping.coordinates_min)./parent.cells_per_dimension + deltas = (parent.mapping.coordinates_max .- parent.mapping.coordinates_min) ./ + parent.cells_per_dimension coordinates_min = parent.mapping.coordinates_min .+ deltas .* (index_min .- 1) coordinates_max = parent.mapping.coordinates_min .+ deltas .* index_max mapping = coordinates2mapping(coordinates_min, coordinates_max) @@ -31,9 +32,11 @@ function StructuredMeshView(parent::StructuredMesh{NDIMS, RealT}; mapping = coordinates2mapping(coordinates_min, coordinates_max) """ - return StructuredMeshView{NDIMS, RealT}(parent, cells_per_dimension, mapping, mapping_as_string, + return StructuredMeshView{NDIMS, RealT}(parent, cells_per_dimension, mapping, + mapping_as_string, parent.current_filename, - index_min, index_max, parent.unsaved_changes) + index_min, index_max, + parent.unsaved_changes) # return StructuredMeshView{NDIMS, RealT}(parent, parent.mapping, # parent.mapping_as_string, @@ -86,7 +89,7 @@ function calc_node_coordinates!(node_coordinates, element, for j in eachnode(basis), i in eachnode(basis) # node_coordinates are the mapped reference node_coordinates node_coordinates[:, i, j, element] .= mapping(cell_x_offset + dx / 2 * nodes[i], - cell_y_offset + dy / 2 * nodes[j]) - end + cell_y_offset + dy / 2 * nodes[j]) + end end end # @muladd From 619382a06f6a48f62a4b64898bb88ed449dac0a9 Mon Sep 17 00:00:00 2001 From: iomsn Date: Mon, 11 Mar 2024 14:33:46 +0000 Subject: [PATCH 22/56] Applied autoformatter. --- src/semidiscretization/semidiscretization_coupled.jl | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/src/semidiscretization/semidiscretization_coupled.jl b/src/semidiscretization/semidiscretization_coupled.jl index ba9528138aa..1bd05a2dd19 100644 --- a/src/semidiscretization/semidiscretization_coupled.jl +++ b/src/semidiscretization/semidiscretization_coupled.jl @@ -466,7 +466,8 @@ function allocate_coupled_boundary_condition(boundary_condition, direction, mesh end # In 2D -function allocate_coupled_boundary_condition(boundary_condition::BoundaryConditionCoupled{2 +function allocate_coupled_boundary_condition(boundary_condition::BoundaryConditionCoupled{ + 2 }, direction, mesh, equations, dg::DGSEM) if direction in (1, 2) From 5ecf1889274b5901b707901db4e091397509e009 Mon Sep 17 00:00:00 2001 From: iomsn Date: Mon, 11 Mar 2024 15:05:39 +0000 Subject: [PATCH 23/56] Cleaned up meshview coupled example. --- .../elixir_advection_smview.jl | 147 ++++++++++-------- 1 file changed, 82 insertions(+), 65 deletions(-) diff --git a/examples/structured_2d_dgsem/elixir_advection_smview.jl b/examples/structured_2d_dgsem/elixir_advection_smview.jl index c285e0adcb3..52715ce5212 100644 --- a/examples/structured_2d_dgsem/elixir_advection_smview.jl +++ b/examples/structured_2d_dgsem/elixir_advection_smview.jl @@ -1,95 +1,112 @@ -# The same setup as tree_2d_dgsem/elixir_advection_basic.jl -# to verify the StructuredMesh implementation against TreeMesh - using OrdinaryDiffEq using Trixi ############################################################################### -# semidiscretization of the linear advection equation +# Coupled semidiscretization of two linear advection systems using converter functions +# and mesh views for the semidiscretizations. First we define a parent mesh +# for the entire physical domain, then we define the two mesh views on the parent. +# +# In this elixir, we have a square domain that is divided into a left and right. +# On each half of the domain, a completely independent SemidiscretizationHyperbolic +# is created for the linear advection equations. The two systems are coupled in the +# x-direction. +# For a high-level overview, see also the figure below: +# +# (-1, 1) ( 1, 1) +# ┌────────────────────┬────────────────────┐ +# │ ↑ periodic ↑ │ ↑ periodic ↑ │ +# │ │ │ +# │ ========= │ ========= │ +# │ system #1 │ system #2 │ +# │ ========= │ ========= │ +# │ │ │ +# │<-- coupled │<-- coupled │ +# │ coupled -->│ coupled -->│ +# │ │ │ +# │ ↓ periodic ↓ │ ↓ periodic ↓ │ +# └────────────────────┴────────────────────┘ +# (-1, -1) ( 1, -1) advection_velocity = (0.2, -0.7) equations = LinearScalarAdvectionEquation2D(advection_velocity) # Create DG solver with polynomial degree = 2 and (local) Lax-Friedrichs/Rusanov flux as surface flux -solver = DGSEM(polydeg = 2, surface_flux = flux_lax_friedrichs) +solver = DGSEM(polydeg = 3)#, surface_flux = flux_lax_friedrichs) -coordinates_min = (-1.0, -1.0) # maximum coordinates (max(x), max(y)) -coordinates_max = (1.0, 1.0) # maximum coordinates (max(x), max(y)) +# Domain size of the parent mesh. +coordinates_min = (-1.0, -1.0) +coordinates_max = (1.0, 1.0) +# Cell dimensions of the parent mesh. cells_per_dimension = (16, 16) -# cells_per_dimension = (32, 16) -# Create curved mesh with 16 x 16 elements -parent = StructuredMesh(cells_per_dimension, coordinates_min, coordinates_max) +# Create drentmesh with 16 x 16 elements +parent_mesh = StructuredMesh(cells_per_dimension, coordinates_min, coordinates_max) -#mesh = StructuredMeshView(parent; index_min = (1, 1), index_max = (16, 16)) -#semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition_convergence_test, solver) +# Create the two mesh views, each of which takes half of the parent mesh. +mesh1 = StructuredMeshView(parent_mesh; index_min = (1, 1), index_max = (8, 16)) +mesh2 = StructuredMeshView(parent_mesh; index_min = (9, 1), index_max = (16, 16)) -mesh1 = StructuredMeshView(parent; index_min = (1, 1), index_max = (8, 16)) -mesh2 = StructuredMeshView(parent; index_min = (9, 1), index_max = (16, 16)) +# The coupling function is simply the identity, as we are dealing with two identical systems. +coupling_function = (x, u, equations_other, equations_own) -> u # Define the coupled boundary conditions boundary_conditions1 = ( # Connect left boundary with right boundary of left mesh - x_neg = BoundaryConditionCoupled(2, (:end, :i_forward), Float64), - x_pos = BoundaryConditionCoupled(2, (:begin, :i_forward), Float64), - # x_neg = boundary_condition_periodic, - # x_pos = boundary_condition_periodic, + x_neg = BoundaryConditionCoupled(2, (:end, :i_forward), Float64, coupling_function), + x_pos = BoundaryConditionCoupled(2, (:begin, :i_forward), Float64, coupling_function), y_neg = boundary_condition_periodic, y_pos = boundary_condition_periodic) boundary_conditions2 = ( # Connect left boundary with right boundary of left mesh - x_neg = BoundaryConditionCoupled(1, (:end, :i_forward), Float64), - x_pos = BoundaryConditionCoupled(1, (:begin, :i_forward), Float64), - # x_neg = boundary_condition_periodic, - # x_pos = boundary_condition_periodic, + x_neg = BoundaryConditionCoupled(1, (:end, :i_forward), Float64, coupling_function), + x_pos = BoundaryConditionCoupled(1, (:begin, :i_forward), Float64, coupling_function), y_neg = boundary_condition_periodic, y_pos = boundary_condition_periodic) # A semidiscretization collects data structures and functions for the spatial discretization -# semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition_convergence_test, solver) semi1 = SemidiscretizationHyperbolic(mesh1, equations, initial_condition_convergence_test, solver, boundary_conditions = boundary_conditions1) -semi2 = SemidiscretizationHyperbolic(mesh2, equations, initial_condition_convergence_test, - solver, - boundary_conditions = boundary_conditions2) -semi = SemidiscretizationCoupled(semi1, semi2) - -############################################################################### -# ODE solvers, callbacks etc. - -# Create ODE problem with time span from 0.0 to 1.0 -ode = semidiscretize(semi, (0.0, 1.0)); - -# At the beginning of the main loop, the SummaryCallback prints a summary of the simulation setup -# and resets the timers -summary_callback = SummaryCallback() - -# The AnalysisCallback allows to analyse the solution in regular intervals and prints the results -#analysis_callback = AnalysisCallback(semi, interval=10) - -# The SaveSolutionCallback allows to save the solution to a file in regular intervals -# save_solution = SaveSolutionCallback(interval=100, -# solution_variables=cons2prim) - -# The StepsizeCallback handles the re-calculation of the maximum Δt after each time step -# stepsize_callback = StepsizeCallback(cfl=1.6) - -# Create a CallbackSet to collect all callbacks such that they can be passed to the ODE solver -# callbacks = CallbackSet(summary_callback, analysis_callback, save_solution, stepsize_callback) -callbacks = CallbackSet(summary_callback, - smview_callback) - -############################################################################### -# run the simulation - -# OrdinaryDiffEq's `solve` method evolves the solution in time and executes the passed callbacks -sol = solve(ode, CarpenterKennedy2N54(williamson_condition = false), - dt = 5.0e-2, # solve needs some value here but it will be overwritten by the stepsize_callback - save_everystep = false, callback = callbacks); - -errors_ref = (l2 = [8.312427642603623e-6], linf = [6.626865824577166e-5]) - -# Print the timer summary -summary_callback() +#semi2 = SemidiscretizationHyperbolic(mesh2, equations, initial_condition_convergence_test, +# solver, +# boundary_conditions = boundary_conditions2) +#semi = SemidiscretizationCoupled(semi1, semi2) + +################################################################################ +## ODE solvers, callbacks etc. +# +## Create ODE problem with time span from 0.0 to 1.0 +#ode = semidiscretize(semi, (0.0, 1.0)); +# +## At the beginning of the main loop, the SummaryCallback prints a summary of the simulation setup +## and resets the timers +#summary_callback = SummaryCallback() +# +## The AnalysisCallback allows to analyse the solution in regular intervals and prints the results +##analysis_callback = AnalysisCallback(semi, interval=10) +# +## The SaveSolutionCallback allows to save the solution to a file in regular intervals +## save_solution = SaveSolutionCallback(interval=100, +## solution_variables=cons2prim) +# +## The StepsizeCallback handles the re-calculation of the maximum Δt after each time step +## stepsize_callback = StepsizeCallback(cfl=1.6) +# +## Create a CallbackSet to collect all callbacks such that they can be passed to the ODE solver +## callbacks = CallbackSet(summary_callback, analysis_callback, save_solution, stepsize_callback) +#callbacks = CallbackSet(summary_callback, +# smview_callback) +# +################################################################################ +## run the simulation +# +## OrdinaryDiffEq's `solve` method evolves the solution in time and executes the passed callbacks +#sol = solve(ode, CarpenterKennedy2N54(williamson_condition = false), +# dt = 5.0e-2, # solve needs some value here but it will be overwritten by the stepsize_callback +# save_everystep = false, callback = callbacks); +# +#errors_ref = (l2 = [8.312427642603623e-6], linf = [6.626865824577166e-5]) +# +## Print the timer summary +#summary_callback() From 781a164954a321a312490f313dc1e8160aa088da Mon Sep 17 00:00:00 2001 From: iomsn Date: Mon, 11 Mar 2024 15:19:46 +0000 Subject: [PATCH 24/56] Added 2d structured mesh views to periodicity checks. --- src/semidiscretization/semidiscretization_hyperbolic.jl | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/src/semidiscretization/semidiscretization_hyperbolic.jl b/src/semidiscretization/semidiscretization_hyperbolic.jl index f61378a7dca..883f86c4b6e 100644 --- a/src/semidiscretization/semidiscretization_hyperbolic.jl +++ b/src/semidiscretization/semidiscretization_hyperbolic.jl @@ -259,7 +259,8 @@ function check_periodicity_mesh_boundary_conditions(mesh::Union{TreeMesh{1}, end function check_periodicity_mesh_boundary_conditions(mesh::Union{TreeMesh{2}, - StructuredMesh{2}}, + StructuredMesh{2}, + StructuredMeshView{2}}, boundary_conditions::Union{NamedTuple, Tuple}) check_periodicity_mesh_boundary_conditions_x(mesh, boundary_conditions[1], From 4bc13bf9aa75c382a22b47a3441b27db25206f28 Mon Sep 17 00:00:00 2001 From: iomsn Date: Mon, 11 Mar 2024 15:32:13 +0000 Subject: [PATCH 25/56] Cleaned up coupled strucutred mesh view example so that it can be used as test. --- .../elixir_advection_smview.jl | 84 ++++++++++--------- 1 file changed, 43 insertions(+), 41 deletions(-) diff --git a/examples/structured_2d_dgsem/elixir_advection_smview.jl b/examples/structured_2d_dgsem/elixir_advection_smview.jl index 52715ce5212..3c45e166929 100644 --- a/examples/structured_2d_dgsem/elixir_advection_smview.jl +++ b/examples/structured_2d_dgsem/elixir_advection_smview.jl @@ -68,45 +68,47 @@ boundary_conditions2 = ( semi1 = SemidiscretizationHyperbolic(mesh1, equations, initial_condition_convergence_test, solver, boundary_conditions = boundary_conditions1) -#semi2 = SemidiscretizationHyperbolic(mesh2, equations, initial_condition_convergence_test, -# solver, -# boundary_conditions = boundary_conditions2) -#semi = SemidiscretizationCoupled(semi1, semi2) +semi2 = SemidiscretizationHyperbolic(mesh2, equations, initial_condition_convergence_test, + solver, + boundary_conditions = boundary_conditions2) +semi = SemidiscretizationCoupled(semi1, semi2) -################################################################################ -## ODE solvers, callbacks etc. -# -## Create ODE problem with time span from 0.0 to 1.0 -#ode = semidiscretize(semi, (0.0, 1.0)); -# -## At the beginning of the main loop, the SummaryCallback prints a summary of the simulation setup -## and resets the timers -#summary_callback = SummaryCallback() -# -## The AnalysisCallback allows to analyse the solution in regular intervals and prints the results -##analysis_callback = AnalysisCallback(semi, interval=10) -# -## The SaveSolutionCallback allows to save the solution to a file in regular intervals -## save_solution = SaveSolutionCallback(interval=100, -## solution_variables=cons2prim) -# -## The StepsizeCallback handles the re-calculation of the maximum Δt after each time step -## stepsize_callback = StepsizeCallback(cfl=1.6) -# -## Create a CallbackSet to collect all callbacks such that they can be passed to the ODE solver -## callbacks = CallbackSet(summary_callback, analysis_callback, save_solution, stepsize_callback) -#callbacks = CallbackSet(summary_callback, -# smview_callback) -# -################################################################################ -## run the simulation -# -## OrdinaryDiffEq's `solve` method evolves the solution in time and executes the passed callbacks -#sol = solve(ode, CarpenterKennedy2N54(williamson_condition = false), -# dt = 5.0e-2, # solve needs some value here but it will be overwritten by the stepsize_callback -# save_everystep = false, callback = callbacks); -# -#errors_ref = (l2 = [8.312427642603623e-6], linf = [6.626865824577166e-5]) -# -## Print the timer summary -#summary_callback() +############################################################################### +# ODE solvers, callbacks etc. + +# Create ODE problem with time span from 0.0 to 1.0 +ode = semidiscretize(semi, (0.0, 1.0)); + +# At the beginning of the main loop, the SummaryCallback prints a summary of the simulation setup +# and resets the timers +summary_callback = SummaryCallback() + +# The AnalysisCallback allows to analyse the solution in regular intervals and prints the results +analysis_callback1 = AnalysisCallback(semi1, interval=100) +analysis_callback2 = AnalysisCallback(semi2, interval=100) +analysis_callback = AnalysisCallbackCoupled(semi, analysis_callback1, analysis_callback2) + +alive_callback = AliveCallback(analysis_interval=analysis_interval) + +# The SaveSolutionCallback allows to save the solution to a file in regular intervals +save_solution = SaveSolutionCallback(interval=100, + save_initial_solution=true, + save_final_solution=true, + solution_variables=cons2prim) + +# The StepsizeCallback handles the re-calculation of the maximum Δt after each time step +stepsize_callback = StepsizeCallback(cfl=1.6) + +# Create a CallbackSet to collect all callbacks such that they can be passed to the ODE solver +callbacks = CallbackSet(summary_callback, analysis_callback, save_solution, stepsize_callback) + +############################################################################### +# run the simulation + +# OrdinaryDiffEq's `solve` method evolves the solution in time and executes the passed callbacks +sol = solve(ode, CarpenterKennedy2N54(williamson_condition = false), + dt = 5.0e-2, # solve needs some value here but it will be overwritten by the stepsize_callback + save_everystep = false, callback = callbacks); + +# Print the timer summary +summary_callback() From 2d7461574589e46790a69cf17fa2ddf856660a86 Mon Sep 17 00:00:00 2001 From: iomsn Date: Mon, 11 Mar 2024 15:42:55 +0000 Subject: [PATCH 26/56] Added 2d structured mesh view to the tests. --- test/test_structured_2d.jl | 33 +++++++++++++++++++++++++++++++++ 1 file changed, 33 insertions(+) diff --git a/test/test_structured_2d.jl b/test/test_structured_2d.jl index 64a1faf05b8..1b60b2a4b74 100644 --- a/test/test_structured_2d.jl +++ b/test/test_structured_2d.jl @@ -70,6 +70,39 @@ end end end +@trixi_testset "elixir_advection_smview.jl" begin + @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_advection_smview.jl"), + l2=[ + 4.5131319539071844e-5, + 4.5131319538970356e-5, + ], + linf=[ + 0.00022262992334731724, + 0.00022262994922361834, + ], + coverage_override=(maxiters = 10^5,)) + + @testset "analysis_callback(sol) for AnalysisCallbackCoupled" begin + errors = analysis_callback(sol) + @test errors.l2≈[ + 4.5131319539071844e-5, + 4.5131319538970356e-5, + ] rtol=1.0e-4 + @test errors.linf≈[ + 0.00022262992334731724, + 0.00022262994922361834, + ] rtol=1.0e-4 + # Ensure that we do not have excessive memory allocations + # (e.g., from type instabilities) + let + t = sol.t[end] + u_ode = sol.u[end] + du_ode = similar(u_ode) + @test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 1000 + end + end +end + @trixi_testset "elixir_advection_extended.jl" begin @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_advection_extended.jl"), l2=[4.220397559713772e-6], From 7e32f7a4dd443f5746af6ff203a0f8c8be2d49f4 Mon Sep 17 00:00:00 2001 From: iomsn Date: Mon, 11 Mar 2024 15:43:12 +0000 Subject: [PATCH 27/56] Added analysis_interval variable. --- examples/structured_2d_dgsem/elixir_advection_smview.jl | 1 + 1 file changed, 1 insertion(+) diff --git a/examples/structured_2d_dgsem/elixir_advection_smview.jl b/examples/structured_2d_dgsem/elixir_advection_smview.jl index 3c45e166929..6ae569e0550 100644 --- a/examples/structured_2d_dgsem/elixir_advection_smview.jl +++ b/examples/structured_2d_dgsem/elixir_advection_smview.jl @@ -88,6 +88,7 @@ analysis_callback1 = AnalysisCallback(semi1, interval=100) analysis_callback2 = AnalysisCallback(semi2, interval=100) analysis_callback = AnalysisCallbackCoupled(semi, analysis_callback1, analysis_callback2) +analysis_interval = 100 alive_callback = AliveCallback(analysis_interval=analysis_interval) # The SaveSolutionCallback allows to save the solution to a file in regular intervals From ad94584c32c551ab9990fceebad3d1303c84d5ad Mon Sep 17 00:00:00 2001 From: iomsn Date: Tue, 12 Mar 2024 11:37:47 +0000 Subject: [PATCH 28/56] Applied updated autoformatter. --- .../elixir_advection_smview.jl | 31 +++++++++++-------- src/callbacks_step/save_solution_dg.jl | 6 ++-- .../semidiscretization_coupled.jl | 4 +-- .../semidiscretization_hyperbolic.jl | 2 +- 4 files changed, 22 insertions(+), 21 deletions(-) diff --git a/examples/structured_2d_dgsem/elixir_advection_smview.jl b/examples/structured_2d_dgsem/elixir_advection_smview.jl index 6ae569e0550..121dc9946cf 100644 --- a/examples/structured_2d_dgsem/elixir_advection_smview.jl +++ b/examples/structured_2d_dgsem/elixir_advection_smview.jl @@ -53,14 +53,18 @@ coupling_function = (x, u, equations_other, equations_own) -> u # Define the coupled boundary conditions boundary_conditions1 = ( # Connect left boundary with right boundary of left mesh - x_neg = BoundaryConditionCoupled(2, (:end, :i_forward), Float64, coupling_function), - x_pos = BoundaryConditionCoupled(2, (:begin, :i_forward), Float64, coupling_function), + x_neg = BoundaryConditionCoupled(2, (:end, :i_forward), Float64, + coupling_function), + x_pos = BoundaryConditionCoupled(2, (:begin, :i_forward), Float64, + coupling_function), y_neg = boundary_condition_periodic, y_pos = boundary_condition_periodic) boundary_conditions2 = ( # Connect left boundary with right boundary of left mesh - x_neg = BoundaryConditionCoupled(1, (:end, :i_forward), Float64, coupling_function), - x_pos = BoundaryConditionCoupled(1, (:begin, :i_forward), Float64, coupling_function), + x_neg = BoundaryConditionCoupled(1, (:end, :i_forward), Float64, + coupling_function), + x_pos = BoundaryConditionCoupled(1, (:begin, :i_forward), Float64, + coupling_function), y_neg = boundary_condition_periodic, y_pos = boundary_condition_periodic) @@ -84,24 +88,25 @@ ode = semidiscretize(semi, (0.0, 1.0)); summary_callback = SummaryCallback() # The AnalysisCallback allows to analyse the solution in regular intervals and prints the results -analysis_callback1 = AnalysisCallback(semi1, interval=100) -analysis_callback2 = AnalysisCallback(semi2, interval=100) +analysis_callback1 = AnalysisCallback(semi1, interval = 100) +analysis_callback2 = AnalysisCallback(semi2, interval = 100) analysis_callback = AnalysisCallbackCoupled(semi, analysis_callback1, analysis_callback2) analysis_interval = 100 -alive_callback = AliveCallback(analysis_interval=analysis_interval) +alive_callback = AliveCallback(analysis_interval = analysis_interval) # The SaveSolutionCallback allows to save the solution to a file in regular intervals -save_solution = SaveSolutionCallback(interval=100, - save_initial_solution=true, - save_final_solution=true, - solution_variables=cons2prim) +save_solution = SaveSolutionCallback(interval = 100, + save_initial_solution = true, + save_final_solution = true, + solution_variables = cons2prim) # The StepsizeCallback handles the re-calculation of the maximum Δt after each time step -stepsize_callback = StepsizeCallback(cfl=1.6) +stepsize_callback = StepsizeCallback(cfl = 1.6) # Create a CallbackSet to collect all callbacks such that they can be passed to the ODE solver -callbacks = CallbackSet(summary_callback, analysis_callback, save_solution, stepsize_callback) +callbacks = CallbackSet(summary_callback, analysis_callback, save_solution, + stepsize_callback) ############################################################################### # run the simulation diff --git a/src/callbacks_step/save_solution_dg.jl b/src/callbacks_step/save_solution_dg.jl index 424821e6b33..7367886ca94 100644 --- a/src/callbacks_step/save_solution_dg.jl +++ b/src/callbacks_step/save_solution_dg.jl @@ -34,8 +34,7 @@ function save_solution_file(u, time, dt, timestep, # compute the solution variables via broadcasting, and reinterpret the # result as a plain array of floating point numbers data = Array(reinterpret(eltype(u), - solution_variables.(reinterpret(SVector{ - nvariables(equations), + solution_variables.(reinterpret(SVector{nvariables(equations), eltype(u)}, u), Ref(equations)))) @@ -117,8 +116,7 @@ function save_solution_file(u, time, dt, timestep, # compute the solution variables via broadcasting, and reinterpret the # result as a plain array of floating point numbers data = Array(reinterpret(eltype(u), - solution_variables.(reinterpret(SVector{ - nvariables(equations), + solution_variables.(reinterpret(SVector{nvariables(equations), eltype(u)}, u), Ref(equations)))) diff --git a/src/semidiscretization/semidiscretization_coupled.jl b/src/semidiscretization/semidiscretization_coupled.jl index 1bd05a2dd19..b8116d3cbd3 100644 --- a/src/semidiscretization/semidiscretization_coupled.jl +++ b/src/semidiscretization/semidiscretization_coupled.jl @@ -466,9 +466,7 @@ function allocate_coupled_boundary_condition(boundary_condition, direction, mesh end # In 2D -function allocate_coupled_boundary_condition(boundary_condition::BoundaryConditionCoupled{ - 2 - }, +function allocate_coupled_boundary_condition(boundary_condition::BoundaryConditionCoupled{2}, direction, mesh, equations, dg::DGSEM) if direction in (1, 2) cell_size = size(mesh, 2) diff --git a/src/semidiscretization/semidiscretization_hyperbolic.jl b/src/semidiscretization/semidiscretization_hyperbolic.jl index 883f86c4b6e..dcd211671c8 100644 --- a/src/semidiscretization/semidiscretization_hyperbolic.jl +++ b/src/semidiscretization/semidiscretization_hyperbolic.jl @@ -260,7 +260,7 @@ end function check_periodicity_mesh_boundary_conditions(mesh::Union{TreeMesh{2}, StructuredMesh{2}, - StructuredMeshView{2}}, + StructuredMeshView{2}}, boundary_conditions::Union{NamedTuple, Tuple}) check_periodicity_mesh_boundary_conditions_x(mesh, boundary_conditions[1], From 3c7265e9408d24210088a59ccfdd422dc73d4ee5 Mon Sep 17 00:00:00 2001 From: iomsn Date: Wed, 13 Mar 2024 11:15:50 +0000 Subject: [PATCH 29/56] Removed unused lines of code. Added comment about muladd. --- src/meshes/structured_mesh_view.jl | 10 ++++------ 1 file changed, 4 insertions(+), 6 deletions(-) diff --git a/src/meshes/structured_mesh_view.jl b/src/meshes/structured_mesh_view.jl index bf143a2b169..b82ac5e44bc 100644 --- a/src/meshes/structured_mesh_view.jl +++ b/src/meshes/structured_mesh_view.jl @@ -1,3 +1,7 @@ +# By default, Julia/LLVM does not use fused multiply-add operations (FMAs). +# Since these FMAs can increase the performance of many numerical algorithms, +# we need to opt-in explicitly. +# See https://ranocha.de/blog/Optimizing_EC_Trixi for further details. @muladd begin #! format: noindent @@ -37,12 +41,6 @@ function StructuredMeshView(parent::StructuredMesh{NDIMS, RealT}; parent.current_filename, index_min, index_max, parent.unsaved_changes) - - # return StructuredMeshView{NDIMS, RealT}(parent, parent.mapping, - # parent.mapping_as_string, - # parent.current_filename, - # index_min, index_max, - # parent.unsaved_changes) end # Check if mesh is periodic From e977aea1be239a81fee69c2a987de8a0356b9d78 Mon Sep 17 00:00:00 2001 From: iomsn Date: Wed, 13 Mar 2024 11:21:41 +0000 Subject: [PATCH 30/56] Added explanatory comments. --- src/meshes/structured_mesh_view.jl | 15 +++++++++++++++ 1 file changed, 15 insertions(+) diff --git a/src/meshes/structured_mesh_view.jl b/src/meshes/structured_mesh_view.jl index b82ac5e44bc..00b6d23089b 100644 --- a/src/meshes/structured_mesh_view.jl +++ b/src/meshes/structured_mesh_view.jl @@ -5,6 +5,11 @@ @muladd begin #! format: noindent +""" + StructuredMeshView{NDIMS, RealT <: Real} <: AbstractMesh{NDIMS} + +A view on a structured curved mesh. +""" mutable struct StructuredMeshView{NDIMS, RealT <: Real} <: AbstractMesh{NDIMS} parent::StructuredMesh{NDIMS, RealT} cells_per_dimension::NTuple{NDIMS, Int} @@ -16,6 +21,16 @@ mutable struct StructuredMeshView{NDIMS, RealT <: Real} <: AbstractMesh{NDIMS} unsaved_changes::Bool end +""" + StructuredMeshView(parent; index_min, index_max) + +Create a StructuredMeshView on a StructuredMesh parent. + +# Arguments +- `parent`: the parent StructuredMesh. +- `index_min`: starting indices of the parent mesh. +- `index_max`: ending indices of the parent mesh. +""" function StructuredMeshView(parent::StructuredMesh{NDIMS, RealT}; index_min = ntuple(_ -> 1, Val(NDIMS)), index_max = size(parent)) where {NDIMS, RealT} From 40e1acb56e6dd5aba3acc58a8340707a9e21d2fb Mon Sep 17 00:00:00 2001 From: iomsn Date: Wed, 13 Mar 2024 11:23:10 +0000 Subject: [PATCH 31/56] Removed unused code. --- src/meshes/mesh_io.jl | 1 - 1 file changed, 1 deletion(-) diff --git a/src/meshes/mesh_io.jl b/src/meshes/mesh_io.jl index 82caccfaea0..ce2ef410aa8 100644 --- a/src/meshes/mesh_io.jl +++ b/src/meshes/mesh_io.jl @@ -129,7 +129,6 @@ function save_mesh_file(mesh::StructuredMeshView, output_directory; system = "", # Create output directory (if it does not exist) mkpath(output_directory) - # filename = joinpath(output_directory, @sprintf("mesh_%s.h5", system)) filename = joinpath(output_directory, @sprintf("mesh_%s_%06d.h5", system, timestep)) # Open file (clobber existing content) From 0aac5dff993ed2012331b25e7f627ae75bc6d6dc Mon Sep 17 00:00:00 2001 From: Simon Candelaresi <10759273+SimonCan@users.noreply.github.com> Date: Mon, 18 Mar 2024 15:45:55 +0000 Subject: [PATCH 32/56] Update examples/structured_2d_dgsem/elixir_advection_smview.jl Co-authored-by: Daniel Doehring --- examples/structured_2d_dgsem/elixir_advection_smview.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/examples/structured_2d_dgsem/elixir_advection_smview.jl b/examples/structured_2d_dgsem/elixir_advection_smview.jl index 121dc9946cf..5f54304652e 100644 --- a/examples/structured_2d_dgsem/elixir_advection_smview.jl +++ b/examples/structured_2d_dgsem/elixir_advection_smview.jl @@ -4,7 +4,7 @@ using Trixi ############################################################################### # Coupled semidiscretization of two linear advection systems using converter functions # and mesh views for the semidiscretizations. First we define a parent mesh -# for the entire physical domain, then we define the two mesh views on the parent. +# for the entire physical domain, then we define the two mesh views on this parent. # # In this elixir, we have a square domain that is divided into a left and right. # On each half of the domain, a completely independent SemidiscretizationHyperbolic From cb22c76172c83fc77cf994d5578651bc791249ef Mon Sep 17 00:00:00 2001 From: Simon Candelaresi <10759273+SimonCan@users.noreply.github.com> Date: Mon, 18 Mar 2024 15:53:46 +0000 Subject: [PATCH 33/56] Update examples/structured_2d_dgsem/elixir_advection_smview.jl Co-authored-by: Daniel Doehring --- examples/structured_2d_dgsem/elixir_advection_smview.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/examples/structured_2d_dgsem/elixir_advection_smview.jl b/examples/structured_2d_dgsem/elixir_advection_smview.jl index 5f54304652e..c7bafa8ec3f 100644 --- a/examples/structured_2d_dgsem/elixir_advection_smview.jl +++ b/examples/structured_2d_dgsem/elixir_advection_smview.jl @@ -6,7 +6,7 @@ using Trixi # and mesh views for the semidiscretizations. First we define a parent mesh # for the entire physical domain, then we define the two mesh views on this parent. # -# In this elixir, we have a square domain that is divided into a left and right. +# In this elixir, we have a square domain that is divided into left and right subdomains. # On each half of the domain, a completely independent SemidiscretizationHyperbolic # is created for the linear advection equations. The two systems are coupled in the # x-direction. From 48716e306b9dd9500994392fc8ccbe85ae3c781c Mon Sep 17 00:00:00 2001 From: Simon Candelaresi <10759273+SimonCan@users.noreply.github.com> Date: Mon, 18 Mar 2024 15:54:10 +0000 Subject: [PATCH 34/56] Update examples/structured_2d_dgsem/elixir_advection_smview.jl Co-authored-by: Daniel Doehring --- examples/structured_2d_dgsem/elixir_advection_smview.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/examples/structured_2d_dgsem/elixir_advection_smview.jl b/examples/structured_2d_dgsem/elixir_advection_smview.jl index c7bafa8ec3f..3d5af3b0e39 100644 --- a/examples/structured_2d_dgsem/elixir_advection_smview.jl +++ b/examples/structured_2d_dgsem/elixir_advection_smview.jl @@ -7,7 +7,7 @@ using Trixi # for the entire physical domain, then we define the two mesh views on this parent. # # In this elixir, we have a square domain that is divided into left and right subdomains. -# On each half of the domain, a completely independent SemidiscretizationHyperbolic +# On each half of the domain, a completely independent `SemidiscretizationHyperbolic` # is created for the linear advection equations. The two systems are coupled in the # x-direction. # For a high-level overview, see also the figure below: From e34e01d1a15a4a81db2215c5d257811e3cf273ae Mon Sep 17 00:00:00 2001 From: Simon Candelaresi <10759273+SimonCan@users.noreply.github.com> Date: Mon, 18 Mar 2024 15:54:23 +0000 Subject: [PATCH 35/56] Update src/meshes/mesh_io.jl Co-authored-by: Daniel Doehring --- src/meshes/mesh_io.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/meshes/mesh_io.jl b/src/meshes/mesh_io.jl index ce2ef410aa8..5757a232ef3 100644 --- a/src/meshes/mesh_io.jl +++ b/src/meshes/mesh_io.jl @@ -123,7 +123,7 @@ end # Does not save the mesh itself to an HDF5 file. Instead saves important attributes # of the mesh, like its size and the type of boundary mapping function. # Then, within Trixi2Vtk, the StructuredMesh and its node coordinates are reconstructured from -# these attributes for plotting purposes +# these attributes for plotting purposes. function save_mesh_file(mesh::StructuredMeshView, output_directory; system = "", timestep = 0) # Create output directory (if it does not exist) From 70435c093f7a980c55f32e4094e36b1c613e3b91 Mon Sep 17 00:00:00 2001 From: Simon Candelaresi <10759273+SimonCan@users.noreply.github.com> Date: Mon, 18 Mar 2024 15:55:54 +0000 Subject: [PATCH 36/56] Update src/meshes/structured_mesh_view.jl Co-authored-by: Daniel Doehring --- src/meshes/structured_mesh_view.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/meshes/structured_mesh_view.jl b/src/meshes/structured_mesh_view.jl index 00b6d23089b..c79634bb8a9 100644 --- a/src/meshes/structured_mesh_view.jl +++ b/src/meshes/structured_mesh_view.jl @@ -37,9 +37,9 @@ function StructuredMeshView(parent::StructuredMesh{NDIMS, RealT}; @assert index_min <= index_max @assert all(index_min .> 0) @assert index_max <= size(parent) - - # Calculate the domain boundaries. cells_per_dimension = index_max .- index_min .+ 1 + + # Compute cells sizes `deltas` deltas = (parent.mapping.coordinates_max .- parent.mapping.coordinates_min) ./ parent.cells_per_dimension coordinates_min = parent.mapping.coordinates_min .+ deltas .* (index_min .- 1) From 4f6a9d2e60a25af15d26fa459476fd07fa518210 Mon Sep 17 00:00:00 2001 From: iomsn Date: Mon, 18 Mar 2024 15:58:57 +0000 Subject: [PATCH 37/56] Corrected comment on solver. --- examples/structured_2d_dgsem/elixir_advection_smview.jl | 2 +- src/meshes/structured_mesh_view.jl | 1 + 2 files changed, 2 insertions(+), 1 deletion(-) diff --git a/examples/structured_2d_dgsem/elixir_advection_smview.jl b/examples/structured_2d_dgsem/elixir_advection_smview.jl index 3d5af3b0e39..dded36769c9 100644 --- a/examples/structured_2d_dgsem/elixir_advection_smview.jl +++ b/examples/structured_2d_dgsem/elixir_advection_smview.jl @@ -30,7 +30,7 @@ using Trixi advection_velocity = (0.2, -0.7) equations = LinearScalarAdvectionEquation2D(advection_velocity) -# Create DG solver with polynomial degree = 2 and (local) Lax-Friedrichs/Rusanov flux as surface flux +# Create DG solver with polynomial degree = 3 solver = DGSEM(polydeg = 3)#, surface_flux = flux_lax_friedrichs) # Domain size of the parent mesh. diff --git a/src/meshes/structured_mesh_view.jl b/src/meshes/structured_mesh_view.jl index c79634bb8a9..66fe1d6c1df 100644 --- a/src/meshes/structured_mesh_view.jl +++ b/src/meshes/structured_mesh_view.jl @@ -37,6 +37,7 @@ function StructuredMeshView(parent::StructuredMesh{NDIMS, RealT}; @assert index_min <= index_max @assert all(index_min .> 0) @assert index_max <= size(parent) + cells_per_dimension = index_max .- index_min .+ 1 # Compute cells sizes `deltas` From d313117ab7ba1feee8918475eb824584f968495d Mon Sep 17 00:00:00 2001 From: iomsn Date: Mon, 18 Mar 2024 16:03:25 +0000 Subject: [PATCH 38/56] Changed order of imported meshes. --- src/Trixi.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/Trixi.jl b/src/Trixi.jl index ca87bf3ec7f..018c1d14ffc 100644 --- a/src/Trixi.jl +++ b/src/Trixi.jl @@ -226,8 +226,8 @@ export entropy, energy_total, energy_kinetic, energy_internal, energy_magnetic, export lake_at_rest_error export ncomponents, eachcomponent -export TreeMesh, StructuredMesh, UnstructuredMesh2D, P4estMesh, T8codeMesh, - StructuredMeshView +export TreeMesh, StructuredMesh, StructuredMeshView, UnstructuredMesh2D, P4estMesh, + T8codeMesh export DG, DGSEM, LobattoLegendreBasis, From 8ebecea9a9e380298538e613a3c7b504c2a4b0c4 Mon Sep 17 00:00:00 2001 From: Simon Candelaresi <10759273+SimonCan@users.noreply.github.com> Date: Mon, 18 Mar 2024 16:04:00 +0000 Subject: [PATCH 39/56] Update src/meshes/structured_mesh_view.jl Co-authored-by: Daniel Doehring --- src/meshes/structured_mesh_view.jl | 1 + 1 file changed, 1 insertion(+) diff --git a/src/meshes/structured_mesh_view.jl b/src/meshes/structured_mesh_view.jl index 66fe1d6c1df..ccefa2b7afb 100644 --- a/src/meshes/structured_mesh_view.jl +++ b/src/meshes/structured_mesh_view.jl @@ -43,6 +43,7 @@ function StructuredMeshView(parent::StructuredMesh{NDIMS, RealT}; # Compute cells sizes `deltas` deltas = (parent.mapping.coordinates_max .- parent.mapping.coordinates_min) ./ parent.cells_per_dimension + # Calculate the domain boundaries. coordinates_min = parent.mapping.coordinates_min .+ deltas .* (index_min .- 1) coordinates_max = parent.mapping.coordinates_min .+ deltas .* index_max mapping = coordinates2mapping(coordinates_min, coordinates_max) From 7b2f0a30eec8de948241897d995891815b6b2620 Mon Sep 17 00:00:00 2001 From: iomsn Date: Mon, 18 Mar 2024 16:34:10 +0000 Subject: [PATCH 40/56] Added comment about coupling interface indexing. --- examples/structured_2d_dgsem/elixir_advection_smview.jl | 2 ++ 1 file changed, 2 insertions(+) diff --git a/examples/structured_2d_dgsem/elixir_advection_smview.jl b/examples/structured_2d_dgsem/elixir_advection_smview.jl index dded36769c9..d9ce64d1ade 100644 --- a/examples/structured_2d_dgsem/elixir_advection_smview.jl +++ b/examples/structured_2d_dgsem/elixir_advection_smview.jl @@ -51,6 +51,8 @@ mesh2 = StructuredMeshView(parent_mesh; index_min = (9, 1), index_max = (16, 16) coupling_function = (x, u, equations_other, equations_own) -> u # Define the coupled boundary conditions +# The indices (:end, :i_forward) and (:begin, :i_forward) denote the interface indexing. +# For a system with coupling in x and y see examples/structured_2d_dgsem/elixir_advection_coupled.jl. boundary_conditions1 = ( # Connect left boundary with right boundary of left mesh x_neg = BoundaryConditionCoupled(2, (:end, :i_forward), Float64, From b68b6c5d1ccbe548fab2af2893e02581b7a57db1 Mon Sep 17 00:00:00 2001 From: iomsn Date: Mon, 18 Mar 2024 17:37:50 +0000 Subject: [PATCH 41/56] Applied autoformatter. --- src/meshes/structured_mesh_view.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/meshes/structured_mesh_view.jl b/src/meshes/structured_mesh_view.jl index ccefa2b7afb..3441ca1a19d 100644 --- a/src/meshes/structured_mesh_view.jl +++ b/src/meshes/structured_mesh_view.jl @@ -39,7 +39,7 @@ function StructuredMeshView(parent::StructuredMesh{NDIMS, RealT}; @assert index_max <= size(parent) cells_per_dimension = index_max .- index_min .+ 1 - + # Compute cells sizes `deltas` deltas = (parent.mapping.coordinates_max .- parent.mapping.coordinates_min) ./ parent.cells_per_dimension From 8d19521dc1bae9d67fa973dec47765c9857f2ef0 Mon Sep 17 00:00:00 2001 From: Simon Candelaresi <10759273+SimonCan@users.noreply.github.com> Date: Tue, 19 Mar 2024 11:12:34 +0000 Subject: [PATCH 42/56] Update src/meshes/structured_mesh_view.jl Co-authored-by: Daniel Doehring --- src/meshes/structured_mesh_view.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/meshes/structured_mesh_view.jl b/src/meshes/structured_mesh_view.jl index 3441ca1a19d..1bbcb363513 100644 --- a/src/meshes/structured_mesh_view.jl +++ b/src/meshes/structured_mesh_view.jl @@ -40,7 +40,7 @@ function StructuredMeshView(parent::StructuredMesh{NDIMS, RealT}; cells_per_dimension = index_max .- index_min .+ 1 - # Compute cells sizes `deltas` + # Compute cell sizes `deltas` deltas = (parent.mapping.coordinates_max .- parent.mapping.coordinates_min) ./ parent.cells_per_dimension # Calculate the domain boundaries. From 7e0ec93f3b0bb3bfe90fab2ed1c3db274ec8da1d Mon Sep 17 00:00:00 2001 From: iomsn Date: Tue, 19 Mar 2024 11:28:35 +0000 Subject: [PATCH 43/56] Renamed index_min and index_max to indices_min and indices_max. --- .../elixir_advection_smview.jl | 4 +- src/meshes/structured_mesh_view.jl | 42 +++++++++---------- 2 files changed, 23 insertions(+), 23 deletions(-) diff --git a/examples/structured_2d_dgsem/elixir_advection_smview.jl b/examples/structured_2d_dgsem/elixir_advection_smview.jl index d9ce64d1ade..b8f8ee953e0 100644 --- a/examples/structured_2d_dgsem/elixir_advection_smview.jl +++ b/examples/structured_2d_dgsem/elixir_advection_smview.jl @@ -44,8 +44,8 @@ cells_per_dimension = (16, 16) parent_mesh = StructuredMesh(cells_per_dimension, coordinates_min, coordinates_max) # Create the two mesh views, each of which takes half of the parent mesh. -mesh1 = StructuredMeshView(parent_mesh; index_min = (1, 1), index_max = (8, 16)) -mesh2 = StructuredMeshView(parent_mesh; index_min = (9, 1), index_max = (16, 16)) +mesh1 = StructuredMeshView(parent_mesh; indices_min = (1, 1), indices_max = (8, 16)) +mesh2 = StructuredMeshView(parent_mesh; indices_min = (9, 1), indices_max = (16, 16)) # The coupling function is simply the identity, as we are dealing with two identical systems. coupling_function = (x, u, equations_other, equations_own) -> u diff --git a/src/meshes/structured_mesh_view.jl b/src/meshes/structured_mesh_view.jl index 1bbcb363513..9a09adefd08 100644 --- a/src/meshes/structured_mesh_view.jl +++ b/src/meshes/structured_mesh_view.jl @@ -16,36 +16,36 @@ mutable struct StructuredMeshView{NDIMS, RealT <: Real} <: AbstractMesh{NDIMS} mapping::Any # Not relevant for performance mapping_as_string::String current_filename::String - index_min::NTuple{NDIMS, Int} - index_max::NTuple{NDIMS, Int} + indices_min::NTuple{NDIMS, Int} + indices_max::NTuple{NDIMS, Int} unsaved_changes::Bool end """ - StructuredMeshView(parent; index_min, index_max) + StructuredMeshView(parent; indices_min, indices_max) Create a StructuredMeshView on a StructuredMesh parent. # Arguments - `parent`: the parent StructuredMesh. -- `index_min`: starting indices of the parent mesh. -- `index_max`: ending indices of the parent mesh. +- `indices_min`: starting indices of the parent mesh. +- `indices_max`: ending indices of the parent mesh. """ function StructuredMeshView(parent::StructuredMesh{NDIMS, RealT}; - index_min = ntuple(_ -> 1, Val(NDIMS)), - index_max = size(parent)) where {NDIMS, RealT} - @assert index_min <= index_max - @assert all(index_min .> 0) - @assert index_max <= size(parent) + indices_min = ntuple(_ -> 1, Val(NDIMS)), + indices_max = size(parent)) where {NDIMS, RealT} + @assert indices_min <= indices_max + @assert all(indices_min .> 0) + @assert indices_max <= size(parent) - cells_per_dimension = index_max .- index_min .+ 1 + cells_per_dimension = indices_max .- indices_min .+ 1 # Compute cell sizes `deltas` deltas = (parent.mapping.coordinates_max .- parent.mapping.coordinates_min) ./ parent.cells_per_dimension # Calculate the domain boundaries. - coordinates_min = parent.mapping.coordinates_min .+ deltas .* (index_min .- 1) - coordinates_max = parent.mapping.coordinates_min .+ deltas .* index_max + coordinates_min = parent.mapping.coordinates_min .+ deltas .* (indices_min .- 1) + coordinates_max = parent.mapping.coordinates_min .+ deltas .* indices_max mapping = coordinates2mapping(coordinates_min, coordinates_max) mapping_as_string = """ coordinates_min = $coordinates_min @@ -56,7 +56,7 @@ function StructuredMeshView(parent::StructuredMesh{NDIMS, RealT}; return StructuredMeshView{NDIMS, RealT}(parent, cells_per_dimension, mapping, mapping_as_string, parent.current_filename, - index_min, index_max, + indices_min, indices_max, parent.unsaved_changes) end @@ -67,21 +67,21 @@ function isperiodic(mesh::StructuredMeshView) end function isperiodic(mesh::StructuredMeshView, dimension) - @unpack parent, index_min, index_max = mesh + @unpack parent, indices_min, indices_max = mesh return (isperiodic(parent, dimension) && - index_min[dimension] == 1 && - index_max[dimension] == size(parent, dimension)) + indices_min[dimension] == 1 && + indices_max[dimension] == size(parent, dimension)) end @inline Base.ndims(::StructuredMeshView{NDIMS}) where {NDIMS} = NDIMS @inline Base.real(::StructuredMeshView{NDIMS, RealT}) where {NDIMS, RealT} = RealT function Base.size(mesh::StructuredMeshView) - @unpack index_min, index_max = mesh - return index_max .- index_min .+ 1 + @unpack indices_min, indices_max = mesh + return indices_max .- indices_min .+ 1 end function Base.size(mesh::StructuredMeshView, i) - @unpack index_min, index_max = mesh - return index_max[i] - index_min[i] + 1 + @unpack indices_min, indices_max = mesh + return indices_max[i] - indices_min[i] + 1 end Base.axes(mesh::StructuredMeshView) = map(Base.OneTo, size(mesh)) Base.axes(mesh::StructuredMeshView, i) = Base.OneTo(size(mesh, i)) From 92c682db043ab3611e69b780ec0ce35fe81767f6 Mon Sep 17 00:00:00 2001 From: iomsn Date: Tue, 19 Mar 2024 11:56:21 +0000 Subject: [PATCH 44/56] Removed relative tolerance from meshview test. --- test/test_structured_2d.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/test/test_structured_2d.jl b/test/test_structured_2d.jl index 1b60b2a4b74..6ca31c9576f 100644 --- a/test/test_structured_2d.jl +++ b/test/test_structured_2d.jl @@ -87,7 +87,7 @@ end @test errors.l2≈[ 4.5131319539071844e-5, 4.5131319538970356e-5, - ] rtol=1.0e-4 + ] @test errors.linf≈[ 0.00022262992334731724, 0.00022262994922361834, From dc93b0eac74c960f96cdcc240e1d8f9b8c5131f8 Mon Sep 17 00:00:00 2001 From: iomsn Date: Tue, 19 Mar 2024 12:02:55 +0000 Subject: [PATCH 45/56] Removed further relative tolerance as it is not needed. --- test/test_structured_2d.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/test/test_structured_2d.jl b/test/test_structured_2d.jl index 6ca31c9576f..ea4988b90f7 100644 --- a/test/test_structured_2d.jl +++ b/test/test_structured_2d.jl @@ -91,7 +91,7 @@ end @test errors.linf≈[ 0.00022262992334731724, 0.00022262994922361834, - ] rtol=1.0e-4 + ] # Ensure that we do not have excessive memory allocations # (e.g., from type instabilities) let From e141342a260ed8b1a0344a22ef2cd7c286591e1e Mon Sep 17 00:00:00 2001 From: Simon Candelaresi <10759273+SimonCan@users.noreply.github.com> Date: Tue, 19 Mar 2024 14:25:29 +0000 Subject: [PATCH 46/56] Update examples/structured_2d_dgsem/elixir_advection_smview.jl Co-authored-by: Daniel Doehring --- examples/structured_2d_dgsem/elixir_advection_smview.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/examples/structured_2d_dgsem/elixir_advection_smview.jl b/examples/structured_2d_dgsem/elixir_advection_smview.jl index b8f8ee953e0..fb7eb762f5b 100644 --- a/examples/structured_2d_dgsem/elixir_advection_smview.jl +++ b/examples/structured_2d_dgsem/elixir_advection_smview.jl @@ -31,7 +31,7 @@ advection_velocity = (0.2, -0.7) equations = LinearScalarAdvectionEquation2D(advection_velocity) # Create DG solver with polynomial degree = 3 -solver = DGSEM(polydeg = 3)#, surface_flux = flux_lax_friedrichs) +solver = DGSEM(polydeg = 3, surface_flux = flux_lax_friedrichs) # Domain size of the parent mesh. coordinates_min = (-1.0, -1.0) From 062d5980ae3705a7d21a8b3492c652559e6f145f Mon Sep 17 00:00:00 2001 From: iomsn Date: Tue, 19 Mar 2024 14:26:42 +0000 Subject: [PATCH 47/56] Fixed typo on parent mesh generation. --- examples/structured_2d_dgsem/elixir_advection_smview.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/examples/structured_2d_dgsem/elixir_advection_smview.jl b/examples/structured_2d_dgsem/elixir_advection_smview.jl index fb7eb762f5b..d8d27031090 100644 --- a/examples/structured_2d_dgsem/elixir_advection_smview.jl +++ b/examples/structured_2d_dgsem/elixir_advection_smview.jl @@ -40,7 +40,7 @@ coordinates_max = (1.0, 1.0) # Cell dimensions of the parent mesh. cells_per_dimension = (16, 16) -# Create drentmesh with 16 x 16 elements +# Create parent mesh with 16 x 16 elements parent_mesh = StructuredMesh(cells_per_dimension, coordinates_min, coordinates_max) # Create the two mesh views, each of which takes half of the parent mesh. From 824fcc6ff3d6ab0bfef19dd33d0ad949f5391163 Mon Sep 17 00:00:00 2001 From: SimonCan Date: Thu, 21 Mar 2024 08:52:32 +0000 Subject: [PATCH 48/56] Applied autoformatter. --- test/test_structured_2d.jl | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/test/test_structured_2d.jl b/test/test_structured_2d.jl index ea4988b90f7..e09e6632352 100644 --- a/test/test_structured_2d.jl +++ b/test/test_structured_2d.jl @@ -84,11 +84,11 @@ end @testset "analysis_callback(sol) for AnalysisCallbackCoupled" begin errors = analysis_callback(sol) - @test errors.l2≈[ + @test errors.l2 ≈ [ 4.5131319539071844e-5, 4.5131319538970356e-5, - ] - @test errors.linf≈[ + ] + @test errors.linf ≈ [ 0.00022262992334731724, 0.00022262994922361834, ] From 6afd4b92b1c2f0786d264a9efcc03f3ed6064285 Mon Sep 17 00:00:00 2001 From: SimonCan Date: Tue, 16 Apr 2024 14:30:58 +0100 Subject: [PATCH 49/56] Updated test results for elixir_advection_smview.jl. --- test/test_structured_2d.jl | 16 ++++++++-------- 1 file changed, 8 insertions(+), 8 deletions(-) diff --git a/test/test_structured_2d.jl b/test/test_structured_2d.jl index 89faecb36c2..b0d21ef9cba 100644 --- a/test/test_structured_2d.jl +++ b/test/test_structured_2d.jl @@ -73,24 +73,24 @@ end @trixi_testset "elixir_advection_smview.jl" begin @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_advection_smview.jl"), l2=[ - 4.5131319539071844e-5, - 4.5131319538970356e-5, + 8.311947673083206e-6, + 8.311947673068427e-6, ], linf=[ - 0.00022262992334731724, - 0.00022262994922361834, + 6.627000273318195e-5, + 6.62700027264096e-5, ], coverage_override=(maxiters = 10^5,)) @testset "analysis_callback(sol) for AnalysisCallbackCoupled" begin errors = analysis_callback(sol) @test errors.l2 ≈ [ - 4.5131319539071844e-5, - 4.5131319538970356e-5, + 8.311947673083206e-6, + 8.311947673068427e-6, ] @test errors.linf ≈ [ - 0.00022262992334731724, - 0.00022262994922361834, + 6.627000273318195e-5, + 6.62700027264096e-5, ] # Ensure that we do not have excessive memory allocations # (e.g., from type instabilities) From 96fbb2b3787e051d18617077056bd2f5d0a6a3d8 Mon Sep 17 00:00:00 2001 From: Simon Candelaresi <10759273+SimonCan@users.noreply.github.com> Date: Mon, 29 Apr 2024 09:38:11 +0100 Subject: [PATCH 50/56] Update src/meshes/mesh_io.jl Co-authored-by: Michael Schlottke-Lakemper --- src/meshes/mesh_io.jl | 2 ++ 1 file changed, 2 insertions(+) diff --git a/src/meshes/mesh_io.jl b/src/meshes/mesh_io.jl index 5757a232ef3..d1746b73335 100644 --- a/src/meshes/mesh_io.jl +++ b/src/meshes/mesh_io.jl @@ -97,6 +97,8 @@ end # of the mesh, like its size and the type of boundary mapping function. # Then, within Trixi2Vtk, the StructuredMesh and its node coordinates are reconstructured from # these attributes for plotting purposes +# Note: the `timestep` argument is needed for compatibility with the method for +# `StructuredMeshView` function save_mesh_file(mesh::StructuredMesh, output_directory; system = "", timestep = 0) # Create output directory (if it does not exist) From c529e856ddb481ff2d488bf3cd5d437c8e8bdc79 Mon Sep 17 00:00:00 2001 From: Simon Candelaresi <10759273+SimonCan@users.noreply.github.com> Date: Mon, 29 Apr 2024 09:39:57 +0100 Subject: [PATCH 51/56] Update src/meshes/mesh_io.jl Co-authored-by: Michael Schlottke-Lakemper --- src/meshes/mesh_io.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/meshes/mesh_io.jl b/src/meshes/mesh_io.jl index d1746b73335..0a4e4f8e0a9 100644 --- a/src/meshes/mesh_io.jl +++ b/src/meshes/mesh_io.jl @@ -282,7 +282,7 @@ function load_mesh_serial(mesh_file::AbstractString; n_cells_max, RealT) end mesh = TreeMesh(SerialTree{ndims}, max(n_cells_max, capacity)) load_mesh!(mesh, mesh_file) - elseif (mesh_type == "StructuredMesh") || (mesh_type == "StructuredMeshView") + elseif mesh_type in ("StructuredMesh", "StructuredMeshView") size_, mapping_as_string = h5open(mesh_file, "r") do file return read(attributes(file)["size"]), read(attributes(file)["mapping"]) From cbdf68590081a0386e492fe3375b611f99d455e9 Mon Sep 17 00:00:00 2001 From: SimonCan Date: Mon, 29 Apr 2024 09:41:57 +0100 Subject: [PATCH 52/56] Renamed elixir_advection_smview.jl to elixir_advection_meshview.jl since 'structured' is redundant for a file in this directory. --- .../{elixir_advection_smview.jl => elixir_advection_meshview.jl} | 0 1 file changed, 0 insertions(+), 0 deletions(-) rename examples/structured_2d_dgsem/{elixir_advection_smview.jl => elixir_advection_meshview.jl} (100%) diff --git a/examples/structured_2d_dgsem/elixir_advection_smview.jl b/examples/structured_2d_dgsem/elixir_advection_meshview.jl similarity index 100% rename from examples/structured_2d_dgsem/elixir_advection_smview.jl rename to examples/structured_2d_dgsem/elixir_advection_meshview.jl From 56394d20bbde882ed06ccf2f0734c6d587bde1d0 Mon Sep 17 00:00:00 2001 From: SimonCan Date: Mon, 29 Apr 2024 09:42:55 +0100 Subject: [PATCH 53/56] Renamed elixir_advection_smview.jl to elixir_advection_meshview.jl. --- test/test_structured_2d.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/test/test_structured_2d.jl b/test/test_structured_2d.jl index b0d21ef9cba..b61053d6011 100644 --- a/test/test_structured_2d.jl +++ b/test/test_structured_2d.jl @@ -70,8 +70,8 @@ end end end -@trixi_testset "elixir_advection_smview.jl" begin - @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_advection_smview.jl"), +@trixi_testset "elixir_advection_meshview.jl" begin + @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_advection_meshview.jl"), l2=[ 8.311947673083206e-6, 8.311947673068427e-6, From 282dc9b68a06d88ca0157b984ed25bcdb910347e Mon Sep 17 00:00:00 2001 From: SimonCan Date: Mon, 29 Apr 2024 10:25:27 +0100 Subject: [PATCH 54/56] Moved save_mesh_file function for StructuredMeshView. --- src/meshes/mesh_io.jl | 23 ----------------------- src/meshes/structured_mesh_view.jl | 23 +++++++++++++++++++++++ 2 files changed, 23 insertions(+), 23 deletions(-) diff --git a/src/meshes/mesh_io.jl b/src/meshes/mesh_io.jl index 0a4e4f8e0a9..9173ba18564 100644 --- a/src/meshes/mesh_io.jl +++ b/src/meshes/mesh_io.jl @@ -122,29 +122,6 @@ function save_mesh_file(mesh::StructuredMesh, output_directory; system = "", return filename end -# Does not save the mesh itself to an HDF5 file. Instead saves important attributes -# of the mesh, like its size and the type of boundary mapping function. -# Then, within Trixi2Vtk, the StructuredMesh and its node coordinates are reconstructured from -# these attributes for plotting purposes. -function save_mesh_file(mesh::StructuredMeshView, output_directory; system = "", - timestep = 0) - # Create output directory (if it does not exist) - mkpath(output_directory) - - filename = joinpath(output_directory, @sprintf("mesh_%s_%06d.h5", system, timestep)) - - # Open file (clobber existing content) - h5open(filename, "w") do file - # Add context information as attributes - attributes(file)["mesh_type"] = get_name(mesh) - attributes(file)["ndims"] = ndims(mesh) - attributes(file)["size"] = collect(size(mesh)) - attributes(file)["mapping"] = mesh.mapping_as_string - end - - return filename -end - # Does not save the mesh itself to an HDF5 file. Instead saves important attributes # of the mesh, like its size and the corresponding `.mesh` file used to construct the mesh. # Then, within Trixi2Vtk, the UnstructuredMesh2D and its node coordinates are reconstructured diff --git a/src/meshes/structured_mesh_view.jl b/src/meshes/structured_mesh_view.jl index 9a09adefd08..925fc10ab2a 100644 --- a/src/meshes/structured_mesh_view.jl +++ b/src/meshes/structured_mesh_view.jl @@ -107,4 +107,27 @@ function calc_node_coordinates!(node_coordinates, element, cell_y_offset + dy / 2 * nodes[j]) end end + +# Does not save the mesh itself to an HDF5 file. Instead saves important attributes +# of the mesh, like its size and the type of boundary mapping function. +# Then, within Trixi2Vtk, the StructuredMesh and its node coordinates are reconstructured from +# these attributes for plotting purposes. +function save_mesh_file(mesh::StructuredMeshView, output_directory; system = "", + timestep = 0) + # Create output directory (if it does not exist) + mkpath(output_directory) + + filename = joinpath(output_directory, @sprintf("mesh_%s_%06d.h5", system, timestep)) + + # Open file (clobber existing content) + h5open(filename, "w") do file + # Add context information as attributes + attributes(file)["mesh_type"] = get_name(mesh) + attributes(file)["ndims"] = ndims(mesh) + attributes(file)["size"] = collect(size(mesh)) + attributes(file)["mapping"] = mesh.mapping_as_string + end + + return filename +end end # @muladd From 29a422b3333ad51fba8077ff0db91325eeddd674 Mon Sep 17 00:00:00 2001 From: Simon Candelaresi <10759273+SimonCan@users.noreply.github.com> Date: Mon, 29 Apr 2024 10:27:40 +0100 Subject: [PATCH 55/56] Update src/meshes/structured_mesh_view.jl Co-authored-by: Michael Schlottke-Lakemper --- src/meshes/structured_mesh_view.jl | 1 - 1 file changed, 1 deletion(-) diff --git a/src/meshes/structured_mesh_view.jl b/src/meshes/structured_mesh_view.jl index 925fc10ab2a..bd55115cc90 100644 --- a/src/meshes/structured_mesh_view.jl +++ b/src/meshes/structured_mesh_view.jl @@ -89,7 +89,6 @@ Base.axes(mesh::StructuredMeshView, i) = Base.OneTo(size(mesh, i)) function calc_node_coordinates!(node_coordinates, element, cell_x, cell_y, mapping, mesh::StructuredMeshView{2}, - # basis::LobattoLegendreBasis) basis) @unpack nodes = basis From 30bb10de31ef941c2e5c5700fa3eaf4b4a8cd8ea Mon Sep 17 00:00:00 2001 From: SimonCan Date: Mon, 29 Apr 2024 12:07:40 +0100 Subject: [PATCH 56/56] Removed redundant check for structured mesh views. --- test/test_structured_2d.jl | 9 --------- 1 file changed, 9 deletions(-) diff --git a/test/test_structured_2d.jl b/test/test_structured_2d.jl index b61053d6011..d124b8f8bc3 100644 --- a/test/test_structured_2d.jl +++ b/test/test_structured_2d.jl @@ -83,15 +83,6 @@ end coverage_override=(maxiters = 10^5,)) @testset "analysis_callback(sol) for AnalysisCallbackCoupled" begin - errors = analysis_callback(sol) - @test errors.l2 ≈ [ - 8.311947673083206e-6, - 8.311947673068427e-6, - ] - @test errors.linf ≈ [ - 6.627000273318195e-5, - 6.62700027264096e-5, - ] # Ensure that we do not have excessive memory allocations # (e.g., from type instabilities) let