-
Notifications
You must be signed in to change notification settings - Fork 113
Commit
This commit does not belong to any branch on this repository, and may belong to a fork outside of the repository.
Add StructuredMeshView as proxy between solver and actual StructuredM…
…esh (#1624) * Add StructuredMeshView as proxy between solver and actual StructuredMesh * Enable StructuredMeshView on submesh * Attempt to using coupled with meshviews. * Update structured_mesh_view.jl Format. * Applied autoformatter on smview elixir. * Applied autoformatter on meshview. * Corected minor typo with StructuredMeshView. * Corrected x-boundaries for smview example elixir. * Removed parent's periodicity for smview. * Removed redundant and problematic redifintion of StructuredMeshView function. * Applied auto formatter on files affected by the structured mesh view. * Applied autoforatter. * Added entries for meshview IO. * Added structuresd mesh view to parametric types. * Added unctionality of writing out mesh files for StructuredMeshView for different time steps. * Added temptative code for dynamically changing mesh view sizes. * Applied autoformatter. * Corrected the calculation of the coordinate mapping for mesh views. * Added cells_per_dimension to the StructuredMeshView. * Added StructuredMeshView to max_dt clculations. * Applied autoformatter. * Applied autoformatter. * Cleaned up meshview coupled example. * Added 2d structured mesh views to periodicity checks. * Cleaned up coupled strucutred mesh view example so that it can be used as test. * Added 2d structured mesh view to the tests. * Added analysis_interval variable. * Applied updated autoformatter. * Removed unused lines of code. Added comment about muladd. * Added explanatory comments. * Removed unused code. * Update examples/structured_2d_dgsem/elixir_advection_smview.jl Co-authored-by: Daniel Doehring <[email protected]> * Update examples/structured_2d_dgsem/elixir_advection_smview.jl Co-authored-by: Daniel Doehring <[email protected]> * Update examples/structured_2d_dgsem/elixir_advection_smview.jl Co-authored-by: Daniel Doehring <[email protected]> * Update src/meshes/mesh_io.jl Co-authored-by: Daniel Doehring <[email protected]> * Update src/meshes/structured_mesh_view.jl Co-authored-by: Daniel Doehring <[email protected]> * Corrected comment on solver. * Changed order of imported meshes. * Update src/meshes/structured_mesh_view.jl Co-authored-by: Daniel Doehring <[email protected]> * Added comment about coupling interface indexing. * Applied autoformatter. * Update src/meshes/structured_mesh_view.jl Co-authored-by: Daniel Doehring <[email protected]> * Renamed index_min and index_max to indices_min and indices_max. * Removed relative tolerance from meshview test. * Removed further relative tolerance as it is not needed. * Update examples/structured_2d_dgsem/elixir_advection_smview.jl Co-authored-by: Daniel Doehring <[email protected]> * Fixed typo on parent mesh generation. * Applied autoformatter. * Updated test results for elixir_advection_smview.jl. * Update src/meshes/mesh_io.jl Co-authored-by: Michael Schlottke-Lakemper <[email protected]> * Update src/meshes/mesh_io.jl Co-authored-by: Michael Schlottke-Lakemper <[email protected]> * Renamed elixir_advection_smview.jl to elixir_advection_meshview.jl since 'structured' is redundant for a file in this directory. * Renamed elixir_advection_smview.jl to elixir_advection_meshview.jl. * Moved save_mesh_file function for StructuredMeshView. * Update src/meshes/structured_mesh_view.jl Co-authored-by: Michael Schlottke-Lakemper <[email protected]> * Removed redundant check for structured mesh views. --------- Co-authored-by: SimonCan <[email protected]> Co-authored-by: Simon Candelaresi <[email protected]> Co-authored-by: iomsn <[email protected]> Co-authored-by: Daniel Doehring <[email protected]>
- Loading branch information
5 people
authored
May 10, 2024
1 parent
ff3106c
commit 8a9fc7b
Showing
17 changed files
with
347 additions
and
39 deletions.
There are no files selected for viewing
122 changes: 122 additions & 0 deletions
122
examples/structured_2d_dgsem/elixir_advection_meshview.jl
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| Original file line number | Diff line number | Diff line change |
|---|---|---|
| @@ -0,0 +1,122 @@ | ||
| using OrdinaryDiffEq | ||
| 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 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` | ||
| # 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 = 3 | ||
| solver = DGSEM(polydeg = 3, surface_flux = flux_lax_friedrichs) | ||
|
|
||
| # 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) | ||
|
|
||
| # 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. | ||
| 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 | ||
|
|
||
| # 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, | ||
| 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), | ||
| y_neg = boundary_condition_periodic, | ||
| y_pos = boundary_condition_periodic) | ||
|
|
||
| # A semidiscretization collects data structures and functions for the spatial discretization | ||
| 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_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 | ||
| 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() |
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| Original file line number | Diff line number | Diff line change |
|---|---|---|
| @@ -0,0 +1,132 @@ | ||
| # 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 | ||
|
|
||
| """ | ||
| 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} | ||
| mapping::Any # Not relevant for performance | ||
| mapping_as_string::String | ||
| current_filename::String | ||
| indices_min::NTuple{NDIMS, Int} | ||
| indices_max::NTuple{NDIMS, Int} | ||
| unsaved_changes::Bool | ||
| end | ||
|
|
||
| """ | ||
| StructuredMeshView(parent; indices_min, indices_max) | ||
| Create a StructuredMeshView on a StructuredMesh parent. | ||
| # Arguments | ||
| - `parent`: the parent StructuredMesh. | ||
| - `indices_min`: starting indices of the parent mesh. | ||
| - `indices_max`: ending indices of the parent mesh. | ||
| """ | ||
| function StructuredMeshView(parent::StructuredMesh{NDIMS, RealT}; | ||
| 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 = 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 .* (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 | ||
| coordinates_max = $coordinates_max | ||
| mapping = coordinates2mapping(coordinates_min, coordinates_max) | ||
| """ | ||
|
|
||
| return StructuredMeshView{NDIMS, RealT}(parent, cells_per_dimension, mapping, | ||
| mapping_as_string, | ||
| parent.current_filename, | ||
| indices_min, indices_max, | ||
| parent.unsaved_changes) | ||
| end | ||
|
|
||
| # Check if mesh is periodic | ||
| function isperiodic(mesh::StructuredMeshView) | ||
| @unpack parent = mesh | ||
| return isperiodic(parent) && size(parent) == size(mesh) | ||
| end | ||
|
|
||
| function isperiodic(mesh::StructuredMeshView, dimension) | ||
| @unpack parent, indices_min, indices_max = mesh | ||
| return (isperiodic(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 indices_min, indices_max = mesh | ||
| return indices_max .- indices_min .+ 1 | ||
| end | ||
| function Base.size(mesh::StructuredMeshView, i) | ||
| @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)) | ||
|
|
||
| function calc_node_coordinates!(node_coordinates, element, | ||
| cell_x, cell_y, mapping, | ||
| mesh::StructuredMeshView{2}, | ||
| basis) | ||
| @unpack nodes = basis | ||
|
|
||
| # Get cell length in reference mesh | ||
| dx = 2 / size(mesh, 1) | ||
| dy = 2 / size(mesh, 2) | ||
|
|
||
| # Calculate node coordinates of reference mesh | ||
| 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 | ||
| 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 |
Oops, something went wrong.