From: Wolfgang Bangerth Date: Sat, 20 Apr 2024 03:48:15 +0000 (-0600) Subject: Some minor updates to step-89. X-Git-Tag: v9.6.0-rc1~343^2 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=6a6a1cac07eee325a7f91df644c1e2db402e957a;p=dealii.git Some minor updates to step-89. --- diff --git a/examples/step-89/doc/intro.dox b/examples/step-89/doc/intro.dox index 666eef93e2..c4dcd07747 100644 --- a/examples/step-89/doc/intro.dox +++ b/examples/step-89/doc/intro.dox @@ -8,7 +8,7 @@ and Magdalena Schreter-Fleischhacker. This tutorial is loosely based on the publication "High-order non-conforming discontinuous Galerkin methods for the acoustic conservation equations" -by Johannes Heinz, Peter Munch and Manfred Kaltenbacher @cite heinz2022high. +by Johannes Heinz, Peter Munch, and Manfred Kaltenbacher @cite heinz2022high. Johannes Heinz was supported by the European Union’s Framework Programme for Research and Innovation Horizon 2020 (2014-2020) under the Marie Skłodowská--Curie Grant @@ -34,8 +34,35 @@ The set of equations consists of the conservation of mass and momentum \frac{\partial \, \mathbf{u}}{\partial \, t} + \frac{1}{\rho}\nabla p = \mathbf{0}. @f] Here, $p$ is the acoustic pressure, $\mathbf{u}$ the acoustic particle velocity, $c$ the -speed of sound, and $\rho$ the mean density of the fluid in which waves are propagating. -For the discretization we make use of discontinuous Galerkin (DG) methods. DG methods are +speed of sound, and $\rho$ the mean density of the fluid in which waves are propagating. As stated +above, the two equations are simply a different way of writing the wave equation: If +you take the time derivative of the first equation, and the divergence of the second, +i.e., compute +@f[ + \frac{\partial^2 \, p}{\partial \, t^2} + \rho c^2 \nabla\cdot \frac{\partial \mathbf{u}}{\partial t} = 0,\\ + \frac{\partial \, \nabla \cdot \mathbf{u}}{\partial \, t} + \nabla \cdot \frac{1}{\rho}\nabla p = \mathbf{0}, +@f] +then you can substitute the second equation into the first one to obtain +@f[ + \frac{\partial^2 \, p}{\partial \, t^2} - \rho c^2 \nabla \cdot \frac{1}{\rho}\nabla p = \mathbf{0}, +@f] +which in the case of constant density $\rho$ results in the more familiar form of the wave +equation: +@f[ + \frac{\partial^2 \, p}{\partial \, t^2} - c^2 \Delta p = \mathbf{0}. +@f] + +The reason one may want to consider the original form is that it has the form of +a +[hyperbolic conservation law](https://en.wikipedia.org/wiki/Hyperbolic_partial_differential_equation#Hyperbolic_system_and_conservation_laws) +in which only first time and spatial derivatives +appear. Whereas both the more familiar, second order form of the wave equation +and the formulation as a first-order system conserve energy, it is often easier to +devise numerical schemes that have the right amount of dissipation (necessary +for numerical stability) using the well-developed machinery available for first-order +systems. + +For the discretization of this form, we make use of discontinuous Galerkin (DG) methods. DG methods are especially attractive for the acoustic conservation equations due to their low numerical dispersion errors. More importantly for this tutorial, DG methods natively extend to non-matching Nitsche-type methods @cite arnold2002unified. I.e., numerical fluxes are not only used on inner @@ -47,12 +74,13 @@ The discretized equations read \int_{\Omega} \mathbf{w}_h\cdot\frac{\partial \,\mathbf{u}_h}{\partial \, t} +\int_{\Omega} \mathbf{w}_h\cdot \frac{1}{\rho} \nabla p_h +\int_{\partial\Omega} \mathbf{w}_h \cdot\mathbf{n} \frac{1}{\rho}(p^*_h-p_h)=\mathbf{0}, @f] where $\mathbf{w}_h$ and $q_h$ are test functions. The numerical fluxes are -defined as follows @cite hochbruck2014efficient +defined as follows @cite hochbruck2014efficient : @f[ p_h^*=p_h-\frac{\tau^-}{\tau^-+\tau^+}[p_h]+\frac{\tau^-\tau^+}{\tau^-+\tau^+}\jump{\mathbf{u}_h},\\ \mathbf{u}_h^*=\mathbf{u}_h-\frac{\gamma^-}{\gamma^-+\gamma^+}[\mathbf{u}_h]+\frac{\gamma^-\gamma^+}{\gamma^-+\gamma^+}\jump{p_h}, @f] with the penalty parameters $\tau=\frac{\rho c}{2}$ and $\gamma=\frac{1}{2\rho c}$. +In these formulas, $[a] = a^- - a^+ $ denotes the jump of an arbitrary quantity $a$ over element faces (face between elements $K^-$ and $K^+$) and $\jump{a} = a^- \mathbf{n}^- + a^+ \mathbf{n}^+$. @@ -62,7 +90,8 @@ For homogeneous materials, the fluxes reduce to standard Lax--Friedrichs fluxes p_h^*=\average{p_h}+\tau\jump{\mathbf{u}_h},\\ \mathbf{u}_h^*=\average{\mathbf{u}_h}+\gamma\jump{p_h}. @f] -$\average{a}=\frac{a^- + a^+}{2}$ denots the averaging operator. +The expression $\average{a}=\frac{a^- + a^+}{2}$ denotes the averaging operator. +

%Non-matching discretizations

@@ -73,27 +102,37 @@ discretization than the surrounding air. In purely acoustic simulations, differe require different element sizes to resolve the same wave because the speed of sound is directly proportional to the wavelength (we will simulate this example later on). -Considering sliding rotating interfaces @cite duerrwaechter2021an also requires the ability to -handle non-matching discretizations: A cylindrical mesh is embedded in a surrounding mesh with a -cylindrical hole. Every time step, the cylinder rotates while the outer mesh remains at the -same position. In this situation it is hardly possible to construct a conforming mesh in every -time step. +A different example of the usefulness of non-matching grids is where one wants to move the mesh +in parts of the domain, but not others. A typical example is the simulation of windmills: +One might want to enclose the rotating wings into a co-rotating mesh (to avoid having to +remesh with every time step) but of course the mesh that describes the air above the +surrounding landscape and around the tower on which the windmill is located should not +rotate. In a case like this, one considers +sliding rotating interfaces @cite duerrwaechter2021an between the co-rotating part of the +mesh and the stationary part of the mesh, but this also requires the ability to +handle non-matching discretizations. Besides this, non-matching methods can be extended to Chimera methods (elements overlap). -Chimera methods can help to reduce the pressure on mesh generation tools since different regions -of a mesh (that may overlap) can be considered independently. +Chimera methods can help reduce the pressure on mesh generation tools since different regions +of a mesh (that may overlap) can be meshed independently. Different methods exist to treat non-matching interfaces. Within this tutorial, we concentrate on two methods: Point-to-point interpolation and Nitsche-type mortaring. +

%Point-to-point interpolation

@image html https://www.dealii.org/images/steps/developer/step_89_intro_point_to_point.svg "" width=25% -Point-to-point interpolation is a naive approach. The points in which values/gradients are -queried in the coupling terms are defined by the quadrature points on the element face of -element $K^-$. As it can be seen from the picture this approach might be subject to aliasing -in some cases. In the picture, information from element $K_1^+$ is completely neglected. +Point-to-point interpolation is a naive approach. Whenever you need to compute integrals over the boundary +of the cell at the left, for a coupled problem you then need to evaluate the solution or +shape functions on the right at quadrature points of the face on the left, i.e., of the face of +element $K^-$. You can just evaluate these be interpolating the information on the right +at these points, but this is in general expensive (read, for example, the documentation of +VectorTools::point_value(), which implements this kind of functionality). As it can be seen +from the picture this approach might be subject to aliasing +in some cases. +

%Nitsche-type mortaring

@@ -102,22 +141,61 @@ in some cases. In the picture, information from element $K_1^+$ is completely ne Mortaring is the process of computing intersections and is not related to the Mortar method which enforces the coupling via Lagrange multipliers. Obtained intersections are also referred to as mortars. On each mortar a new integration rule is defined. The integral of the face of element -$K^-$ is computed on the intersections. This way, the numerical integration is exact as long -as a sufficient number of integration points is used. -In this tutorial, the intersections are computed using `CGAL`. Therefore, `deal.II` has -to be configured with `DEAL_II_WITH_CGAL` for the Nitsche-type mortaring implementation. +$K^-$ is computed on the intersections. The idea is that if we want to integrate something +over a face $f\subset \partial K^-$, that we break that integral into pieces: +@f[ + \int_f \cdots dx = \sum_i \int_{f_i} \cdots dx +@f] +where each of the $f_i$ corresponds to the intersection of the original face $f$ +on the left with exactly one of the faces on the right; or, if we had multiple +faces on the left, then the $f_i$ would be the intersections of exactly one face +on the left and one face on the right. + +The point of this approach is first, that splitting the integral this way is +exact. Secondly, and maybe more importantly, the terms we are integrating +(the dots in the formula above) are now defined on one cell on each side, and +consequently are smooth (whereas a finite element solution considered across +multiple cells is, in general, *not* smooth). As a consequence, if we replace +the integrals by the numerical integration (quadrature), then the +result is exact as long as a sufficient number of integration points is used +(at least for affine element shapes; for general curved elements, the integrand +will contain rational expressions that are difficult to integrate exactly). + +In this tutorial, the intersections are computed using +[CGAL, the Computational Geometry Algorithms Library](https://www.cgal.org/). Therefore, `deal.II` has +to be configured with `DEAL_II_WITH_CGAL` for the Nitsche-type mortaring implementation. See the deal.II +[Readme file](https://dealii.org/current/readme.html) for information about installation. +

%FERemoteEvaluation

-FERemoteEvaluation is a wrapper class which provides a similar interface to, e.g., FEEvaluation to +In practice, for integrals as those mentioned above, we need to evaluate solutions +(and shape functions) from cells across the non-matching interface. This is awkward +enough if the other side is on the same processor, but outright difficult if the +cells on the other side of the interface are owned by a different process in a parallel +computation. + +On regular meshes (say, doing things as we do in step-40), this is addressed by making +sure that we are only computing integrals on +@ref GlossLocallyOwnedCell "locally owned cells" and keeping around one layer of +@ref GlossGhostCell "ghost cells" for which we can query information. +Ghost cells are the neighbors of locally owned cells, but in cases like +the picture above, where the meshes are not matching, the cells on the +other side of the interface are not *neighbors* in the logical sense -- though +they happen to be geometrically located adjacently. As a consequence, +we need to find a way to efficiently query information on cells that +are perhaps located on a different process. + +FERemoteEvaluation is a wrapper class which provides a similar interface to, e.g., the +FEEvaluation and FEFaceEvaluations classes to access values over non-matching interfaces in matrix-free loops. A detailed description on how to setup -the class and how to use it in actual code is given below on hands-on examples. Within this tutorial we only +the class and how to use it in actual code is given below using hands-on examples. Within this tutorial we only show the usage for non-matching discretizations. Note however, that FERemoteEvaluation can also be used in other settings such as volume coupling. Under the hood, Utilities::MPI::RemotePointEvaluation is used to query the solution or gradients at quadrature points. A detailed description how this is done can be found in step-87. The main difference in the usage of FERemoteEvaluation compared to FEEvaluation is that the interpolated values/gradients of the finite element solution at all the quadrature points are precomputed globally before -the loop that passes through the cells/faces of the mesh (i.e., near the place where the communication takes place) +the loop over the cells/faces of the mesh (i.e., near the place where the communication takes place) instead of performing the interpolation out of the vector on a cell-by-cell basis. The standard code to evaluate fluxes via FEEvaluation reads: @@ -142,7 +220,10 @@ const auto face_function = } }; -matrix_free.template loop({}, face_function, {}, dst, src); +matrix_free.template loop(/* cell_operation= */{}, + /* interior_face_operation= */ face_function, + /* boundary_face_operation= */{}, + dst, src); @endcode The code to evaluate fluxes via FERemoteEvaluation is shown below.