-
Notifications
You must be signed in to change notification settings - Fork 112
Commit
This commit does not belong to any branch on this repository, and may belong to a fork outside of the repository.
- Loading branch information
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.