From b96a7cd2c82bc5b30820651c0c0a5dbf8025c23e Mon Sep 17 00:00:00 2001 From: Matthias Maier Date: Thu, 5 Dec 2019 16:06:55 -0600 Subject: [PATCH] initial commit --- examples/step-69/CMakeLists.txt | 55 + examples/step-69/doc/builds-on | 0 examples/step-69/doc/intro.dox | 367 +++++ examples/step-69/doc/kind | 1 + examples/step-69/doc/results.dox | 0 examples/step-69/doc/tooltip | 0 examples/step-69/step-69.cc | 2159 ++++++++++++++++++++++++++++++ 7 files changed, 2582 insertions(+) create mode 100644 examples/step-69/CMakeLists.txt create mode 100644 examples/step-69/doc/builds-on create mode 100644 examples/step-69/doc/intro.dox create mode 100644 examples/step-69/doc/kind create mode 100644 examples/step-69/doc/results.dox create mode 100644 examples/step-69/doc/tooltip create mode 100644 examples/step-69/step-69.cc diff --git a/examples/step-69/CMakeLists.txt b/examples/step-69/CMakeLists.txt new file mode 100644 index 0000000000..ae880b8f4a --- /dev/null +++ b/examples/step-69/CMakeLists.txt @@ -0,0 +1,55 @@ +## +# CMake script +## + +# Set the name of the project and target: +SET(TARGET "step-69") + +# Declare all source files the target consists of. Here, this is only +# the one step-X.cc file, but as you expand your project you may wish +# to add other source files as well. If your project becomes much larger, +# you may want to either replace the following statement by something like +# FILE(GLOB_RECURSE TARGET_SRC "source/*.cc") +# FILE(GLOB_RECURSE TARGET_INC "include/*.h") +# SET(TARGET_SRC ${TARGET_SRC} ${TARGET_INC}) +# or switch altogether to the large project CMakeLists.txt file discussed +# in the "CMake in user projects" page accessible from the "User info" +# page of the documentation. +SET(TARGET_SRC + ${TARGET}.cc + ) + +# Usually, you will not need to modify anything beyond this point... + +CMAKE_MINIMUM_REQUIRED(VERSION 2.8.12) + +FIND_PACKAGE(deal.II 9.2.0 QUIET + HINTS ${deal.II_DIR} ${DEAL_II_DIR} ../ ../../ $ENV{DEAL_II_DIR} + ) +IF(NOT ${deal.II_FOUND}) + MESSAGE(FATAL_ERROR "\n" + "*** Could not locate a (sufficiently recent) version of deal.II. ***\n\n" + "You may want to either pass a flag -DDEAL_II_DIR=/path/to/deal.II to cmake\n" + "or set an environment variable \"DEAL_II_DIR\" that contains this path." + ) +ENDIF() + +# +# Are all dependencies fulfilled? +# +IF(NOT DEAL_II_WITH_MPI OR NOT DEAL_II_WITH_P4EST) # keep in one line + MESSAGE(FATAL_ERROR " +Error! This tutorial requires a deal.II library that was configured with the following options: + DEAL_II_WITH_MPI = ON + DEAL_II_WITH_P4EST = ON +However, the deal.II library found at ${DEAL_II_PATH} was configured with these options + DEAL_II_WITH_MPI = ${DEAL_II_WITH_MPI} + DEAL_II_WITH_P4EST = ${DEAL_II_WITH_P4EST} +which conflict with the requirements." + ) +ENDIF() + +DEAL_II_INITIALIZE_CACHED_VARIABLES() +SET(CLEAN_UP_FILES *.log *.vtu *.pvtu) +PROJECT(${TARGET}) +DEAL_II_INVOKE_AUTOPILOT() diff --git a/examples/step-69/doc/builds-on b/examples/step-69/doc/builds-on new file mode 100644 index 0000000000..e69de29bb2 diff --git a/examples/step-69/doc/intro.dox b/examples/step-69/doc/intro.dox new file mode 100644 index 0000000000..8d03e560ff --- /dev/null +++ b/examples/step-69/doc/intro.dox @@ -0,0 +1,367 @@ + + This program was contributed by Matthias Maier (Texas A&M University), + and Ignacio Tomas (Sandia National Laboratories, Albuquerque). + + +@dealiiTutorialDOI{10.5281/zenodo.TODO,https://zenodo.org/badge/DOI/10.5281/zenodo.TODO.svg} + +@note This tutorial step implements a first-order accurate guaranteed +maximum wavespeed method based on a first-order graph viscosity +for solving Euler's equations of gas dynamics @cite GuermondPopov2016. As +such it is presented primarily for educational purposes. For actual +research computations you might want to consider exploring a corresponding +high-performance implementation of a second-order accurate scheme that uses +convex limiting techniques, and strong stability-preserving (SSP) +time integration @cite GuermondEtAl2018. The repository can be found at +TODO. + + + + + + + +

Introduction

+ +This tutorial presents a first-order scheme for solving compressible +Euler's equations that is based on three ingredients: a +collocation-type discretization of Euler's equations in context of +finite elements; a graph-viscosity stabilization based on a +guaranteed upper bound of the local wave speed; and explicit +time-stepping. As such the ideas and techniques presented in this tutorial +step are drastically different from those used in Step-33, which focuses on +the use of automatic differentiation. From a programming perspective this +tutorial will focus on a number of techniques found in large-scale +computations: hybrid thread- and process- (MPI) parallelization; efficient +local numbering of degrees of freedom; concurrent post-processing and +write-out of results using worker threads; as well as checkpointing and +restart. + +It should be noted that first-order schemes in the context of hyperbolic +conservation laws require prohibitively many degrees of freedom to resolve +certain key features of the simulated fluid, and thus, typically only serve +as elementary building blocks in higher-order schemes +@cite GuermondEtAl2018. However, we hope that the reader still finds the +tutorial step to be a good starting point (in particular with respect to +the programming techniques) before jumping into full research codes such as +the second-order scheme @cite GuermondEtAl2018 maintained TODO. + + + + + + + + +

Euler's equations of gas dynamics

+ +The compressible Euler's equations of gas dynamics are written in +conservative form as follows: + +@f{align} +\mathbf{u}_t + \text{div} \, \mathbb{f}(\mathbf{u}) = \boldsymbol{0} , +@f} + +where $\mathbf{u}(\textbf{x},t):\mathbb{R}^{d} \times \mathbb{R} +\rightarrow \mathbb{R}^{d+2}$, and $\mathbb{f}(\mathbf{u}):\mathbb{R}^{d+2} +\rightarrow \mathbb{R}^{(d+2) \times d}$, and $d \geq 1$ is the space +dimension. We say that $\mathbf{u} \in \mathbb{R}^{d+2}$ is the state and +$\mathbb{f}(\mathbf{u}) \in \mathbb{R}^{(d+2) \times d}$ is the flux of +the system. In the case of Euler's equations the state is given by +$\textbf{u} = [\rho, \textbf{m},E]^{\top}$: where $\rho \in \mathbb{R}^+$ +denotes the density, $\textbf{m} \in \mathbb{R}^d$ is the momentum, and $E +\in \mathbb{R}^+$ is the total energy of the system. The flux of the system +$\mathbb{f}(\mathbf{u})$ is defined as + +@f{align*} +\mathbb{f}(\textbf{u}) += +\begin{bmatrix} + \textbf{m}^\top \\ + \rho^{-1} \textbf{m} \otimes \textbf{m} + \mathbb{I} p\\ + \tfrac{\textbf{m}^\top}{\rho} (E + p) +\end{bmatrix}, +@f} + +where $\mathbb{I} \in \mathbb{R}^{d \times d}$ is the identity matrix and +$\otimes$ denotes the tensor product. Here, we have introduced the pressure +$p$ that, in general, is defined by an closed-form equation of state. +For the tutorial we limit the discussion to the class of polytropic ideal gases +for which the pressure is given by +@f{align*} +p = p(\textbf{u}) := (\gamma -1) \Big(E - \frac{\|\textbf{m}\|^2}{2\,\rho} +\Big), +@f} + +where the factor $\gamma \in (1,5/3]$ denotes the +ratio of +specific heats, and $\|\,.\|$ denotes the Euclidian norm. + + + + + + +

Solution theory

+ +Hyperbolic conservation laws, such as +@f{align*} +\mathbf{u}_t + \text{div} \, \mathbb{f}(\mathbf{u}) = \boldsymbol{0}, +@f} +pose a significant challenge with respect to solution theory. An evident +observation is that rewriting the equation in variational form and testing with +the solution itself does not lead to an energy estimate because the pairing +$\langle \text{div} \, \mathbb{f}(\mathbf{u}), \mathbf{u}\rangle$ (understood as +the $L^2(\Omega)$ inner product or duality pairing) is not guaranteed to be +non-negative. Notions such as energy-stability or $L^2(\Omega)$-stability are +(in general) meaningles in this context. + +Historically, the most fruitful step taken in order to deepen the +understanding of hyperbolic conservation laws was to assume that the +solution is formally defined as $\mathbf{u} := \lim_{\epsilon \rightarrow +0^+} \mathbf{u}^{\epsilon}$ where $\mathbf{u}^{\epsilon}$ is the solution +of the parabolic regularization + +@f{align} +\mathbf{u}_t^{\epsilon} + \text{div} \, \mathbb{f}(\mathbf{u}^{\epsilon}) += {\epsilon} \Delta \mathbf{u}^{\epsilon}. +@f} + +Such solutions, which are understood as the solution recovered in the +zero-viscosity limit, are often refered to as viscosity solutions. +Global existence and uniqueness of such solutions is a widely open issue. +However, we know at least that if such viscosity solutions exists they have +to satisfy the constraint $\textbf{u}(\mathbf{x},t) \in \mathcal{B}$ for +all $\mathbf{x} \in \Omega$ and $t \geq 0$ where + +@f{align} + \mathcal{B} = \big\{ \textbf{u} = + [\rho, \textbf{m},E]^{\top} \in \mathbb{R}^{d+2} \, \big | + \quad + \rho > 0, + \quad + \ E - \tfrac{|\textbf{m}|_{\ell^2}^2}{2 \rho} > 0, + \quad + s(\mathbf{u}) \geq \min_{x \in \Omega} s(\mathbf{u}_0(\mathbf{x})) + \big\}. +@f} + +Here, $s(\mathbf{u})$ denotes the specific entropy + +@f{align} + s(\mathbf{u}) = \ln \Big(\frac{p(\mathbf{u})}{\rho^{\gamma}}\Big). +@f} + +In other words, a state $\mathbf{u}(\mathbf{x},t)\in\mathcal{B}$ obeys +positivity of the density, positivity of the internal energy, and a local +minimum principle on the specific entropy. This condition is a simplified +version of a class of pointwise stability constraints satisfied by the +exact (viscosity) solution. By pointwise we mean that the constraint has to +be satisfied at every point of the domain, not just in an averaged +(integral, or high order moments) sense. + +In context of a numerical approximation, a violation of such a constraint +has dire consequences: it almost surely leads to catrastrophic failure of +the numerical scheme; loss of hyperbolicity, and overall, loss of +well-posedness of the (discrete) problem. In the following we will +formulate a scheme that ensures that the discrete approximation of +$\mathbf{u}(\mathbf{x},t)$ remains in $\mathcal{B}$. + + + + + + +

Variational versus collocation-type discretizations

+ +Following Step-9, Step-12, and Step-33, at this point it might look tempting +to base a discretization of Euler's equations on a (semi-discrete) variational +formulation: + +@f{align*} + (\partial_t\mathbf{u}_{h},\textbf{v}_h)_{L^2(\Omega)} + - ( \mathbb{f}(\mathbf{u}_{h}) ,\text{grad} \, \textbf{v}_{h})_{L^2(\Omega)} + + s_h(\mathbf{u}_{h},\textbf{v}_h)_{L^2(\Omega)} = \boldsymbol{0} + \quad\forall \textbf{v}_h \in \mathbb{V}_h. +@f} + +Here, $\mathbb{V}_h$ is an appropriate finite element space, and +$s_h(\cdot,\cdot)_{L^2(\Omega)}$ is some linear stabilization method +(possibly complemented with some ad-hoc shock-capturing technique, see for +instance @cite GuermondErn2004 Chapter 5 and references therein). Most +time-dependent discretization approaches described in the deal.II tutorials +are based on such a (semi-discrete) variational approach. Fundamentally, +from an analysis perspective, variational discretizations are conceived in +order to provide some notion of global (integral) stabiliy, meaning an estimate +of the form + +@f{align*} + |\!|\!| \mathbf{u}_{h}(t) |\!|\!| \leq |\!|\!| \mathbf{u}_{h}(0) |\!|\!| +@f} + +holds true, where $|\!|\!| \cdot |\!|\!| $ could represent the +$L^2(\Omega)$-norm or, more generally, some discrete (possibly mesh +dependent) energy-norm. Variational discretizations of hyperbolic +conservation laws have been very popular since the mid eighties, in +particular combined with SUPG-type stabilization and/or upwinding +techniques (see the early work of @cite Brooks1982 and @cite Johnson1986). They +have proven to be some of the best approaches for simulations in the subsonic +shockless regime and similarly benign regimes. + +However, in the transonic and supersonic regime, and shock-hydrodynamics +applications the use of variational schemes might be questionable. In fact, at +the time of this writing, most shock-hydrodynamics codes are still firmly +grounded on finite volumes methods. The main reason for failure of variational +schemes in such extreme regimes is the lack of pointwise stability. This stems +from the fact that a priori bounds on integrated quantities (e.g. +integrals of moments) have in general no implications on pointwise properties +of the solution. While some of these problems might be alleviated by the +(perpetual) chase of the right shock capturing scheme, finite difference-like +and finite volume schemes still have an edge in many regards. + +In this tutorial step we therefore depart from variational schemes. We will +present a completely algebraic formulation (with the flavor of a +collocation-type scheme) that preserves constraints pointwise, i.e., + +@f{align*} + \textbf{u}_h(\mathbf{x}_i,t) \in \mathcal{B} + \;\text{at every node}\;\mathbf{x}_i\;\text{of the mesh}. +@f} + +However, contrary to finite difference/volume schemes, the +scheme implemented in this step maximizes the use of finite element +software infrastructure, works in any mesh, in any space dimension, and is +theoretically guaranteed to always work, all the time, no exception. This +illustrates that deal.ii can be used far beyond the context of variational +schemes in Hilbert spaces and that a large number of classes, modules and +namespaces from deal.ii can be adapted for such purpose. + + + + + +

Description of the scheme

+ +Let $\mathbb{V}_h$ be scalar-valued finite dimensional space spanned by a +basis $\{\phi_i\}_{i \in \mathcal{V}}$ where: $\phi_i:\Omega \rightarrow +\mathbb{R}$ and $\mathcal{V}$ is the set of all indices (nonnegative +integers) identifying each scalar Degree of Freedom (DOF) in the mesh. +Therefore a scalar finite element functiona $u_h \in \mathbb{V}_h$ it can +be written as $u_h = \sum_{i \in \mathcal{V}} U_i \phi_i$ with $U_i \in +\mathbb{R}$. We introduce the notation for vector-valued approximation +spaces $\pmb{\mathbb{V}}_h := \{\mathbb{V}_h\}^{d+2}$. Let $\mathbf{u}_h +\in \pmb{\mathbb{V}}_h$, then it can be written as $\mathbf{u}_h = \sum_{i +\in \mathcal{V}} \mathbf{U}_i \phi_i$ where $\mathbf{U}_i \in +\mathbb{R}^{d+2}$ and $\phi_i$ is a scalar-valued shape function. + +Note. +For simplicity we will consider the usual Lagrange finite elements. In such +context $\{\mathbf{x}_i\}_{i \in \mathcal{V}}$ be the set of all "support +points" (see @ref GlossSupport "this glossary entry") where $\mathbf{x}_i \in +\mathbb{R}^d$. Then each integer index $i \in \mathcal{V}$ +uniquely identifies a support point $\mathbf{x}_i$ and/or scalar-valued shape +function $\phi_i$. + +With this notation we can define the scheme as + +@f{align*} + m_i \frac{\mathbf{U}_i^{n+1} - \mathbf{U}_i^{n}}{\tau} + + \sum_{j \in \mathcal{I}(i)} \mathbb{f}(\mathbf{U}_j^{n})\cdot + \mathbf{c}_{ij} - d_{ij} \mathbf{U}_j^{n} = \boldsymbol{0} +@f} + +Where + - $m_i := \int_{\Omega} \phi_i \, \mathrm{d}\mathbf{x}$ + - $\tau$ is the time step size + - $\mathbf{c}_{ij} := \int_{\Omega} \nabla\phi_j\phi_i \, + \mathrm{d}\mathbf{x}$ (note that $\mathbf{c}_{ij}\in \mathbb{R}^d$) + - $\mathcal{I}(i) := \{j \in \mathcal{V} \ | \ \mathbf{c}_{ij} \not \equiv + \boldsymbol{0}\} \cup \{i\}$. We will refer to $\mathcal{I}(i)$ as the + "stencil" (or adjacency list) at the support point $i$. + - $\mathbb{f}(\mathbf{U}_j^{n})$ is just the flux $\mathbb{f}$ of the + hyperbolic system evaluated at the state $\mathbf{U}_j^{n}$ stored at the + support point $j$. + - $d_{ij} := \max \{ \lambda_{\text{max}} + (\mathbf{U}_i^{n},\mathbf{U}_j^{n}, \textbf{n}_{ij}), + \lambda_{\text{max}} (\mathbf{U}_j^{n}, \mathbf{U}_i^{n}, + \textbf{n}_{ji}) \} \|\mathbf{c}_{ij}\|_{\ell^2} $ + - $\textbf{n}_{ij} = \frac{\mathbf{c}_{ij}}{ \|\mathbf{c}_{ij}\|_{\ell^2} }$ + +The definition of $\lambda_{\text{max}} (\mathbf{U},\mathbf{V}, +\textbf{n})$ is far from trivial and we will postpone their definition in +order to focus on the computational/coding issues of this tutorial Step. +For the time being let's note that + - $m_i$ and $\mathbf{c}_{ij}$ do not evolve in time. It makes sense to + compute and store them once, and later recall them at very time step. + They are part of what we are going to call off-line data. + - At every time step we have to evaluate $\mathbb{f}(\mathbf{U}_j^{n})$ and + $d_{ij} := \max \{ \lambda_{\text{max}} + (\mathbf{U}_i^{n},\mathbf{U}_j^{n}, \textbf{n}_{ij}), + \lambda_{\text{max}} (\mathbf{U}_j^{n}, \mathbf{U}_i^{n}, + \textbf{n}_{ji}) \} \|\mathbf{c}_{ij}\|_{\ell^2} $ + +Before we start with the description of the implementation of this scheme, it +is worth saying a thing or two about the "assembly" of this system. Consider +for instance a hypothetical pseudo-code, illustrating +a possible strategy to compute the solution $\textbf{U}^{n+1}$: + +@f{align*} +&\textbf{For } i \in \mathcal{V} \\ +&\ \ \ \ \{\mathbf{c}_{ij}\}_{j \in \mathcal{I}(i)} := +\texttt{gather_cij_vectors}(\textbf{c}, \mathcal{I}(i)) \\ +&\ \ \ \ \{\textbf{U}_j^n\}_{j \in \mathcal{I}(i)} := +\texttt{gather_state_vectors}(\textbf{U}^n, \mathcal{I}(i)) \\ +&\ \ \ \ \ \textbf{U}_i^{n+1} := \mathbf{U}_i^{n} \\ +&\ \ \ \ \textbf{For } j \in \mathcal{I}(i) \\ +&\ \ \ \ \ \ \ \ \textbf{U}_i^{n+1} := \textbf{U}_i^{n+1} - \frac{\tau}{m_i} + \mathbb{f}(\mathbf{U}_j^{n})\cdot + \mathbf{c}_{ij} + d_{ij} \mathbf{U}_j^{n} \\ +&\ \ \ \ \textbf{EndFor} \\ +&\textbf{EndFor} +@f} +We note here that: +- This "assembly" does not require any form of quadrature or cell-loops. +- Here $\textbf{c}$ and $\textbf{U}^n$ are a global matrix and a global vector +containing all the vectors $\mathbf{c}_{ij}$ and all the states +$\mathbf{U}_j^n$ respectively. +- $\texttt{gather_cij_vectors}$ and $\texttt{gather_state_vectors}$ are +hypothetical implementations that collect (from global matrices and vectors) +only the quantities required to compute the update at the node $i$. +- Note that: if we assume a cartesian mesh in two space +dimensions, first-order polynomial space $\mathbb{Q}^1$, and that +$\mathbf{x}_i$ is an interior node (i.e. $\mathbf{x}_i$ is not on the boundary +of the domain ) then: $\{\textbf{U}_j^n\}_{j \in \mathcal{I}(i)}$ should contain +nine state-vectors (i.e. all the states in the patch/macro element associated to +the shape function $\phi_i$). This is one of the major differences with the +usual cell-based loop where the gather functionality (encoded in +FEValuesBase.get_function_values() ) only collects values for the +local cell (just a subset of the patch). + +It is worth noting that, from a practitioner's point of view +fully-algebraic schemes (i.e. no bilinear forms, no cell-loops, and no +quadratures) are not unusual at all in the CFD community. There is rich history +of application of this kind of schemes, also called "edge-based" or +"graph-based" finite element schemes (see for instance @cite Rainald2008 for +more historical references). + +This pseudo-code was introduced only to prepare the mindset of the reader for +what is going to be presented in the in the next section. The +actual implementation described in the next section is somewhat different from +what is described in the pseudo-code but shares the same core mentality: we do +not loop on cells but rather we loop on the edges of the sparsity graph (hence +the name "edge-based" code) in order to assemble the system. + + + + + +

Implementation of the scheme

+ + + + + + + + + diff --git a/examples/step-69/doc/kind b/examples/step-69/doc/kind new file mode 100644 index 0000000000..86a44aa1ef --- /dev/null +++ b/examples/step-69/doc/kind @@ -0,0 +1 @@ +time dependent diff --git a/examples/step-69/doc/results.dox b/examples/step-69/doc/results.dox new file mode 100644 index 0000000000..e69de29bb2 diff --git a/examples/step-69/doc/tooltip b/examples/step-69/doc/tooltip new file mode 100644 index 0000000000..e69de29bb2 diff --git a/examples/step-69/step-69.cc b/examples/step-69/step-69.cc new file mode 100644 index 0000000000..1a95bb76e5 --- /dev/null +++ b/examples/step-69/step-69.cc @@ -0,0 +1,2159 @@ +/* --------------------------------------------------------------------- + * + * Copyright (C) 2012 - 2019 by the deal.II authors + * + * This file is part of the deal.II library. + * + * The deal.II library is free software; you can use it, redistribute + * it, and/or modify it under the terms of the GNU Lesser General + * Public License as published by the Free Software Foundation; either + * version 2.1 of the License, or (at your option) any later version. + * The full text of the license can be found in the file LICENSE.md at + * the top level directory of deal.II. + * + * --------------------------------------------------------------------- + + * + * Authors: Matthias Maier, Texas A&M University; + * Ignacio Tomas, Texas A&M University, Sandia National Laboratories + */ + +// @sect3{Include files} +// The set of include files is quite standard. The most intriguing part at this +// point in time is that: either though this code is a "thread and mpi parallel" +// we are using neither Trilinos nor PETSC vectors. Actually we are using dealii +// distributed vectors la_parallel_vector.h and the regular dealii +// sparse matrices sparse_matrix.h +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include +#include +#include +#include +#include + +#include + + +// @sect3{Declaration/s of the namespace Step69} + +namespace Step69 +{ + enum Boundary : dealii::types::boundary_id + { + do_nothing = 0, + slip = 2, + dirichlet = 3, + }; + + + // @sect4{Declaration of Discretization class template} + + // The main goal of this class is to digest the input file and act as a + // "container" of members that may be changed/decided at run time (through the + // input file). It was natural to derive this class from + // dealii::ParameterAcceptor. This class is in charge of + // initializing mpi comunicator, geometry dimensions, triangulation, mapping, + // finite element, mapping, and quadratures. If we think of the class + // Discretization as a "container": it doesn't contain any + // memmory demanding class member such a dof_handlers, vectors or matrices. + // The most memmory thirsty class member is the + // dealii::parallel::distributed::Triangulation. + + template + class Discretization : public dealii::ParameterAcceptor + { + public: + Discretization(const MPI_Comm & mpi_communicator, + dealii::TimerOutput &computing_timer, + const std::string & subsection = "Discretization"); + + void setup(); + + const MPI_Comm &mpi_communicator; + + dealii::parallel::distributed::Triangulation triangulation; + + const dealii::MappingQ mapping; + const dealii::FE_Q finite_element; + const dealii::QGauss quadrature; + const dealii::QGauss face_quadrature; + + private: + dealii::TimerOutput &computing_timer; + + double immersed_disc_length; + double immersed_disc_height; + double immersed_disc_object_position; + double immersed_disc_object_diameter; + + unsigned int refinement; + }; + + // @sect4{Declaration of OfflineData class template} + + // The class OfflineData is initializes (and "owns") + // pretty much all the components of the discretization that + // do not evolve in time. In particular: dof_handler, sparsity + // patterns, boundary maps, lumped mass matrix, and other matrices + // and vectors that do not change in time are members of this class. + // The term "offline" here refers to the idea that all the class members + // of OfflineData are initialized and assigned values + // a "time step zero" and are not meant to be modified at any other later + // time step. For instance, the sparsity pattern should not + // change as we advance in time (we are not doing any form of adaptivity in + // space). Similarly, the entries of the vector + // lumped_mass_matrix should not be modified as we advance in + // time either. + // + // Placeholder: Say something about BoundaryNormalMap. + + template + class OfflineData : public dealii::ParameterAcceptor + { + public: + using BoundaryNormalMap = + std::map /* normal */, + dealii::types::boundary_id /* id */, + dealii::Point> /* position */>; + + OfflineData(const MPI_Comm & mpi_communicator, + dealii::TimerOutput & computing_timer, + const Discretization &discretization, + const std::string & subsection = "OfflineData"); + + void setup(); + void assemble(); + + dealii::DoFHandler dof_handler; + + std::shared_ptr partitioner; + + unsigned int n_locally_owned; + unsigned int n_locally_relevant; + + dealii::SparsityPattern sparsity_pattern; + + BoundaryNormalMap boundary_normal_map; + + dealii::SparseMatrix lumped_mass_matrix; + std::array, dim> cij_matrix; + std::array, dim> nij_matrix; + dealii::SparseMatrix norm_matrix; + + private: + const MPI_Comm & mpi_communicator; + dealii::TimerOutput &computing_timer; + + dealii::SmartPointer> discretization; + }; + + // @sect4{Declaration of ProblemDescription class template} + + // Most of the implementations of the members of this class will be utility + // classes/functions specific for Euler's equations: + // - The type alias rank1_type will be used for the states + // $\mathbf{U}_i^n$ + // - The type alias rank2_type will be used for the fluxes + // $\mathbb{f}(\mathbf{U}_j^n)$. + // - The implementation of momentum will extract $\textbf{m}$ + // (out of the state vector $[\rho,\textbf{m},E]$) and store it in a + // Tensor<1, dim> for our convenience. + // - The implementation of internal_energy will compute + // $E - \frac{|\textbf{m}|^2}{2\rho}$ from the state vector + // $[\rho,\textbf{m},E]$. + // + // The purpose of the remaining class members component_names, + // pressure, and speed_of_sound, + // is evident from their names. Most notably, the last + // one compute_lambda_max is in charge of computing + // $\lambda_{max}(\mathbf{U},\mathbf{V},\mathbf{n})$ which is required + // to compute the first order viscosity $d_{ij}$ as detailed in the section + // Description of the scheme. + + template + class ProblemDescription + { + public: + static constexpr unsigned int problem_dimension = 2 + dim; + + using rank1_type = dealii::Tensor<1, problem_dimension>; + using rank2_type = + dealii::Tensor<1, problem_dimension, dealii::Tensor<1, dim>>; + + const static std::array component_names; + + static constexpr double gamma = 7. / 5.; + + static DEAL_II_ALWAYS_INLINE inline dealii::Tensor<1, dim> + momentum(const rank1_type U); + + static DEAL_II_ALWAYS_INLINE inline double + internal_energy(const rank1_type U); + + static DEAL_II_ALWAYS_INLINE inline double pressure(const rank1_type U); + + static DEAL_II_ALWAYS_INLINE inline double + speed_of_sound(const rank1_type U); + + static DEAL_II_ALWAYS_INLINE inline rank2_type f(const rank1_type U); + + static DEAL_II_ALWAYS_INLINE inline double + compute_lambda_max(const rank1_type U_i, + const rank1_type U_j, + const dealii::Tensor<1, dim> &n_ij); + }; + + // @sect4{Declaration of InitialValues class template} + + // Placeholder here + + template + class InitialValues : public dealii::ParameterAcceptor + { + public: + using rank1_type = typename ProblemDescription::rank1_type; + + InitialValues(const std::string &subsection = "InitialValues"); + + std::function &point, double t)> + initial_state; + + private: + void parse_parameters_callback(); + + dealii::Tensor<1, dim> initial_direction; + dealii::Tensor<1, 3> initial_1d_state; + }; + + // @sect4{Declaration of TimeStep class template} + + // Placeholder here + + template + class TimeStep : public dealii::ParameterAcceptor + { + public: + static constexpr unsigned int problem_dimension = + ProblemDescription::problem_dimension; + + using rank1_type = typename ProblemDescription::rank1_type; + using rank2_type = typename ProblemDescription::rank2_type; + + typedef std::array, + problem_dimension> + vector_type; + + TimeStep(const MPI_Comm & mpi_communicator, + dealii::TimerOutput & computing_timer, + const OfflineData & offline_data, + const InitialValues &initial_values, + const std::string & subsection = "TimeStep"); + + void prepare(); + + double step(vector_type &U, double t); + + private: + const MPI_Comm & mpi_communicator; + dealii::TimerOutput &computing_timer; + + dealii::SmartPointer> offline_data; + dealii::SmartPointer> initial_values; + + dealii::SparseMatrix dij_matrix; + + vector_type temp; + + double cfl_update; + }; + + // @sect4{Declaration of SchlierenPostprocessor class template} + + // At its core, the Schilieren class implements the class member + // compute_schlieren. The main purpose of this class member + // is to compute auxiliary finite element field schlieren + // at each node, defined as + // \f[ \text{schlieren}[i] = e^{\beta \frac{ |\nabla r_i| + // - \min_j |\nabla r_j| }{\max_j |\nabla r_j| - \min_j |\nabla r_j| } } \f] + // where $r$ in principle could be any scalar finite element field. + // The natural candidate is choosing $r := \rho$. Schlieren postprocessing + // is a standard methodology to enhance the contrast of the visualization + // inspired in actual X-ray and shadowgraphy experimental techniques of + // visualization. + + template + class SchlierenPostprocessor : public dealii::ParameterAcceptor + { + public: + static constexpr unsigned int problem_dimension = + ProblemDescription::problem_dimension; + + using rank1_type = typename ProblemDescription::rank1_type; + + using vector_type = + std::array, + problem_dimension>; + + SchlierenPostprocessor( + const MPI_Comm & mpi_communicator, + dealii::TimerOutput & computing_timer, + const OfflineData &offline_data, + const std::string & subsection = "SchlierenPostprocessor"); + + void prepare(); + + void compute_schlieren(const vector_type &U); + + dealii::LinearAlgebra::distributed::Vector schlieren; + + private: + const MPI_Comm & mpi_communicator; + dealii::TimerOutput &computing_timer; + + dealii::SmartPointer> offline_data; + + dealii::Vector r; + + unsigned int schlieren_index; + double schlieren_beta; + }; + + // @sect4{Declaration of TimeLoop class template} + + // Placeholder here + + template + class TimeLoop : public dealii::ParameterAcceptor + { + public: + using vector_type = typename TimeStep::vector_type; + + TimeLoop(const MPI_Comm &mpi_comm); + + void run(); + + private: + vector_type interpolate_initial_values(double t = 0); + + void output(const vector_type &U, + const std::string &name, + double t, + unsigned int cycle, + bool checkpoint = false); + + const MPI_Comm & mpi_communicator; + std::ostringstream timer_output; + dealii::TimerOutput computing_timer; + + dealii::ConditionalOStream pcout; + + std::string base_name; + double t_final; + double output_granularity; + bool enable_compute_error; + + bool resume; + + Discretization discretization; + OfflineData offline_data; + InitialValues initial_values; + TimeStep time_step; + SchlierenPostprocessor schlieren_postprocessor; + + std::unique_ptr filestream; + + std::thread output_thread; + vector_type output_vector; + }; + +} // namespace Step69 + + +// @sect3{Implementation of the classes in namespace Step69} + +namespace Step69 +{ + using namespace dealii; + + // @sect4{Implementation of the members of the class Discretization} + + // Not much is done here other that initializing the corresponding + // class members in the initialization list. + + template + Discretization::Discretization(const MPI_Comm & mpi_communicator, + TimerOutput & computing_timer, + const std::string &subsection) + : ParameterAcceptor(subsection) + , mpi_communicator(mpi_communicator) + , triangulation(mpi_communicator) + , mapping(1) + , finite_element(1) + , quadrature(3) + , face_quadrature(3) + , computing_timer(computing_timer) + { + immersed_disc_length = 4.; + add_parameter("immersed disc - length", + immersed_disc_length, + "Immersed disc: length of computational domain"); + + immersed_disc_height = 2.; + add_parameter("immersed disc - height", + immersed_disc_height, + "Immersed disc: height of computational domain"); + + immersed_disc_object_position = 0.6; + add_parameter("immersed disc - object position", + immersed_disc_object_position, + "Immersed disc: x position of immersed disc center point"); + + immersed_disc_object_diameter = 0.5; + add_parameter("immersed disc - object diameter", + immersed_disc_object_diameter, + "Immersed disc: diameter of immersed disc"); + + refinement = 5; + add_parameter("initial refinement", + refinement, + "Initial refinement of the geometry"); + } + + // Note that in the previous constructor we only passed the MPI + // communicator to the triangulationbut we still have not + // initialized the underlying geometry/mesh. In order to define the geometry + // we will use the class create_immersed_disc_geometry + // that uses the tools in GridGenerator in order to create a + // rectangular domain with a whole. + + // The following is just a dummy implementation/placeholder that does + // nothing other than throwing an exception if we want to run this program + // with a space dimension that is not 2. + + template + void + create_immersed_disc_geometry(parallel::distributed::Triangulation &, + const double /*length*/, + const double /*height*/, + const double /*step_position*/, + const double /*step_height*/) + { + AssertThrow(false, ExcNotImplemented()); + } + + // For the two-dimensional case we have the following template + // specialization that creates the geometry. + + template <> + void create_immersed_disc_geometry<2>( + parallel::distributed::Triangulation<2> &triangulation, + const double length, + const double height, + const double disc_position, + const double disc_diameter) + { + constexpr int dim = 2; + + Triangulation tria1, tria2, tria3, tria4; + + GridGenerator::hyper_cube_with_cylindrical_hole( + tria1, disc_diameter / 2., disc_diameter, 0.5, 1, false); + + GridGenerator::subdivided_hyper_rectangle( + tria2, + {2, 1}, + Point<2>(-disc_diameter, disc_diameter), + Point<2>(disc_diameter, height / 2.)); + + GridGenerator::subdivided_hyper_rectangle( + tria3, + {2, 1}, + Point<2>(-disc_diameter, -disc_diameter), + Point<2>(disc_diameter, -height / 2.)); + + GridGenerator::subdivided_hyper_rectangle( + tria4, + {6, 4}, + Point<2>(disc_diameter, -height / 2.), + Point<2>(length - disc_position, height / 2.)); + + GridGenerator::merge_triangulations({&tria1, &tria2, &tria3, &tria4}, + triangulation, + 1.e-12, + true); + + triangulation.set_manifold(0, PolarManifold<2>(Point<2>())); + + for (auto cell : triangulation.active_cell_iterators()) + for (unsigned int v = 0; v < GeometryInfo::vertices_per_cell; ++v) + { + auto &vertex = cell->vertex(v); + if (vertex[0] <= -disc_diameter + 1.e-6) + vertex[0] = -disc_position; + } + + for (auto cell : triangulation.active_cell_iterators()) + { + for (unsigned int f = 0; f < GeometryInfo<2>::faces_per_cell; ++f) + { + const auto face = cell->face(f); + + if (!face->at_boundary()) + continue; + + const auto center = face->center(); + + if (center[0] > length - disc_position - 1.e-6) + face->set_boundary_id(Boundary::do_nothing); + else if (center[0] < -disc_position + 1.e-6) + face->set_boundary_id(Boundary::dirichlet); + else + face->set_boundary_id(Boundary::slip); + } + } + } + + // This the the last class member to be implemented of the class. + // Discretization: it initializes the actual mesh + // of the triangulation by calling create_immersed_disc_geometry. + + template + void Discretization::setup() + { + TimerOutput::Scope t(computing_timer, "discretization - setup"); + + triangulation.clear(); + + create_immersed_disc_geometry(triangulation, + immersed_disc_length, + immersed_disc_height, + immersed_disc_object_position, + immersed_disc_object_diameter); + + triangulation.refine_global(refinement); + } + + + // @sect4{Implementation of the members of the class OfflineData} + + // Not much is done here other that initializing the corresponding + // class members in the initialization list. + // Constructor of the class OfflineData. + + template + OfflineData::OfflineData(const MPI_Comm & mpi_communicator, + TimerOutput & computing_timer, + const Discretization &discretization, + const std::string & subsection) + : ParameterAcceptor(subsection) + , mpi_communicator(mpi_communicator) + , computing_timer(computing_timer) + , discretization(&discretization) + {} + + // Now the class member OfflineData::setup() will take care + // of initializating + // - The dof_handler. + // - The IndexSets corresponding to locally owned and locally relevant DOFs. + // - The partitioner. + + template + void OfflineData::setup() + { + IndexSet locally_owned; + IndexSet locally_relevant; + + { + TimerOutput::Scope t(computing_timer, "offline_data - distribute dofs"); + + dof_handler.initialize(discretization->triangulation, + discretization->finite_element); + DoFRenumbering::Cuthill_McKee(dof_handler); + + locally_owned = dof_handler.locally_owned_dofs(); + n_locally_owned = locally_owned.n_elements(); + } + + { + TimerOutput::Scope t( + computing_timer, + "offline_data - create partitioner and affine constraints"); + + locally_relevant.clear(); + DoFTools::extract_locally_relevant_dofs(dof_handler, locally_relevant); + n_locally_relevant = locally_relevant.n_elements(); + + partitioner.reset(new Utilities::MPI::Partitioner(locally_owned, + locally_relevant, + mpi_communicator)); + } + + const auto dofs_per_cell = discretization->finite_element.dofs_per_cell; + + // Here we create the sparsity patterns for the off-line data. There are + // quite a few peculiarities that deserve our attention: + // - Our "local" dynamic sparsity pattern (dsp) + // will be of dimensions + // n_locally_relevant $\times$ + // n_locally_relevant (this choice is definitely unusual). + // The goal behind such choice is to reduce communication: we do + // not want to request (to another mpi process) ghosted offline data (such + // as the vectors $\mathbf{c}_{ij}$ when $j$ is not locally owned) for + // every time step. It is more efficient to simply take more memory by + // storing (locally) all relevant off-line data. + // - We loop on all locally owned and ghosted cells (see @ref + // GlossArtificialCell "this glossary entry") order to: + //
    + //
  • Extract the dof_indices associated to the cell DOFs + // (having global numbering) and renumber them using + // partitioner->global_to_local(index). For the case + // of locally owned DOFs: such renumbering consist in applying a + // shift (i.e. we subtract a number) such that now they will + // become a number in the integer interval + // $[0,$n_locally_owned$)$. However, for the case of + // "ghosted DOFs" (i.e. not locally owned) the situation is quite + // different, since the global indices associated to ghosted DOFs + // will not be (in general) a contiguous set of integers. + //
    • + //
  • Once, we are done with that, we add the corresponding entries to + // the rows of the dynamic sparsity pattern with + // dsp.add_entries
    • + //
+ // Finally we use dsp to initialize the actual sparsity + // pattern sparsity_pattern. + + { + TimerOutput::Scope t(computing_timer, + "offline_data - create sparsity pattern"); + + DynamicSparsityPattern dsp(n_locally_relevant, n_locally_relevant); + + std::vector dof_indices(dofs_per_cell); + + for (auto cell : dof_handler.active_cell_iterators()) + { + if (cell->is_artificial()) + continue; + + cell->get_dof_indices(dof_indices); + std::transform(dof_indices.begin(), + dof_indices.end(), + dof_indices.begin(), + [&](auto index) { + return partitioner->global_to_local(index); + }); + + for (const auto dof : dof_indices) + dsp.add_entries(dof, dof_indices.begin(), dof_indices.end()); + } + + sparsity_pattern.copy_from(dsp); + } + + // We initialize the off-line data matrices. Note that these matrices do + // not require an mpi communicator (that's the idea). + + { + TimerOutput::Scope t(computing_timer, "offline_data - set up matrices"); + + lumped_mass_matrix.reinit(sparsity_pattern); + norm_matrix.reinit(sparsity_pattern); + for (auto &matrix : cij_matrix) + matrix.reinit(sparsity_pattern); + for (auto &matrix : nij_matrix) + matrix.reinit(sparsity_pattern); + } + } + + // In this last brace we finished with the implementation of the + // OfflineData::setup(). + // + // Now we define a collection of assembly utilities: + // - CopyData: This will only be used to compute the off-line + // data using WorkStream. It acts as a container: it is just a + // struct where WorkStream stores the local cell contributions. + // - get_entry: it reads the value stored at the entry + // pointed by the iterator it of matrix + // - set_entry: it sets value at the entry + // pointed by the iterator it of matrix. + // - gather_get_entry: we note that + // $\mathbf{c}_{ij} \in \mathbb{R}^d$. If $d=2$ then + // $\mathbf{c}_{ij} = [\mathbf{c}_{ij}^1,\mathbf{c}_{ij}^2]^\top$. + // Which basically implies + // that we need one matrix per space dimension to store the + // $\mathbf{c}_{ij}$ vectors. Similar observation follows for the matrix + // $\mathbf{n}_{ij}$. The purpose of gather_get_entry + // is to retrieve those entries a store them into a + // Tensor<1, dim> for our convenience. + // - gather (first interface): + // - gather (second interface): + // - scatter: + + namespace + { + template + struct CopyData + { + bool is_artificial; + std::vector local_dof_indices; + typename OfflineData::BoundaryNormalMap local_boundary_normal_map; + dealii::FullMatrix cell_lumped_mass_matrix; + std::array, dim> cell_cij_matrix; + }; + + + template + DEAL_II_ALWAYS_INLINE inline typename Matrix::value_type + get_entry(const Matrix &matrix, const Iterator &it) + { + const auto global_index = it->global_index(); + const typename Matrix::const_iterator matrix_iterator(&matrix, + global_index); + return matrix_iterator->value(); + } + + + template + inline DEAL_II_ALWAYS_INLINE void + set_entry(Matrix & matrix, + const Iterator & it, + typename Matrix::value_type value) + { + const auto global_index = it->global_index(); + typename Matrix::iterator matrix_iterator(&matrix, global_index); + matrix_iterator->value() = value; + } + + + template + DEAL_II_ALWAYS_INLINE inline dealii::Tensor<1, k> + gather_get_entry(const std::array &U, const T2 it) + { + dealii::Tensor<1, k> result; + for (unsigned int j = 0; j < k; ++j) + result[j] = get_entry(U[j], it); + return result; + } + + + template + DEAL_II_ALWAYS_INLINE inline dealii::Tensor<1, k> + gather(const std::array &U, const T2 i, const T3 l) + { + dealii::Tensor<1, k> result; + for (unsigned int j = 0; j < k; ++j) + result[j] = U[j](i, l); + return result; + } + + + template + DEAL_II_ALWAYS_INLINE inline dealii::Tensor<1, k> + gather(const std::array &U, const T2 i) + { + dealii::Tensor<1, k> result; + for (unsigned int j = 0; j < k; ++j) + result[j] = U[j].local_element(i); + return result; + } + + + template + DEAL_II_ALWAYS_INLINE inline void + scatter(std::array &U, const T2 &result, const T3 i) + { + for (unsigned int j = 0; j < k1; ++j) + U[j].local_element(i) = result[j]; + } + } // end of namespace. + + // The following piece of code implements the class member + // OfflineData::assemble() which (in short) + // computes the lumped mass entries $m_i$, the vectors $\mathbf{c}_{ij}$, + // the vector $\mathbf{n}_{ij} = \frac{\mathbf{c}_{ij}}{|\mathbf{c}_{ij}|}$, + // and the boundary normals. The information about boundary normals is + // collected into the map BoundaryNormalMap: which maps the + // global index of the DOF/node into the tuple + // $\{\text{normal}, \text{boundary id},\text{position} \}$. + // + // In order to exploit thread parallelization we use WorkStream approach + // detailed in the @ref threads "Parallel computing with multiple processors + // accessing shared memory". As customary this requires + // definition of + // - Scratch data: in this case it is scratch_data. + // - The worker: in the case it is local_assemble_system that + // actually computes the local (i.e. current cell) contributions. + // - A copy data: a struct that contains all the local assembly + // contributions, in this case called CopyData(). + // - A copy data routine: in this case it is + // copy_local_to_global in charge of actually coping these + // local contributions into the global objects (matrices and/or vectors) + // + // Most the following lines are spent in the definition of the worker + // local_assemble_system and the copy routine + // copy_local_to_global. There is not much to say about the + // WorkStream framework since the vast majority of ideas are reasonably + // well-documented in Step-9, Step-13 and Step-32 among others. + + template + void OfflineData::assemble() + { + lumped_mass_matrix = 0.; + norm_matrix = 0.; + for (auto &matrix : cij_matrix) + matrix = 0.; + for (auto &matrix : nij_matrix) + matrix = 0.; + + const unsigned int dofs_per_cell = + discretization->finite_element.dofs_per_cell; + const unsigned int n_q_points = discretization->quadrature.size(); + + /* This is the implementation of the scratch data required by WorkStream */ + MeshWorker::ScratchData scratch_data( + discretization->mapping, + discretization->finite_element, + discretization->quadrature, + update_values | update_gradients | update_quadrature_points | + update_JxW_values, + discretization->face_quadrature, + update_normal_vectors | update_values | update_JxW_values); + + { + TimerOutput::Scope t( + computing_timer, + "offline_data - assemble lumped mass matrix, and c_ij"); + + /* This is the implementation of the "worker" required by WorkStream */ + const auto local_assemble_system = [&](const auto &cell, + auto & scratch, + auto & copy) { + auto &is_artificial = copy.is_artificial; + auto &local_dof_indices = copy.local_dof_indices; + auto &local_boundary_normal_map = copy.local_boundary_normal_map; + auto &cell_lumped_mass_matrix = copy.cell_lumped_mass_matrix; + auto &cell_cij_matrix = copy.cell_cij_matrix; + + is_artificial = cell->is_artificial(); + if (is_artificial) + return; + + local_boundary_normal_map.clear(); + cell_lumped_mass_matrix.reinit(dofs_per_cell, dofs_per_cell); + for (auto &matrix : cell_cij_matrix) + matrix.reinit(dofs_per_cell, dofs_per_cell); + + const auto &fe_values = scratch.reinit(cell); + + local_dof_indices.resize(dofs_per_cell); + cell->get_dof_indices(local_dof_indices); + + std::transform(local_dof_indices.begin(), + local_dof_indices.end(), + local_dof_indices.begin(), + [&](auto index) { + return partitioner->global_to_local(index); + }); + + for (unsigned int q_point = 0; q_point < n_q_points; ++q_point) + { + const auto JxW = fe_values.JxW(q_point); + + for (unsigned int j = 0; j < dofs_per_cell; ++j) + { + const auto value_JxW = fe_values.shape_value(j, q_point) * JxW; + const auto grad_JxW = fe_values.shape_grad(j, q_point) * JxW; + + cell_lumped_mass_matrix(j, j) += value_JxW; + + for (unsigned int i = 0; i < dofs_per_cell; ++i) + { + const auto value = fe_values.shape_value(i, q_point); + for (unsigned int d = 0; d < dim; ++d) + cell_cij_matrix[d](i, j) += (value * grad_JxW)[d]; + + } /* for i */ + } /* for j */ + } /* for q */ + + for (unsigned int f = 0; f < GeometryInfo::faces_per_cell; ++f) + { + const auto face = cell->face(f); + const auto id = face->boundary_id(); + + if (!face->at_boundary()) + continue; + + const auto &fe_face_values = scratch.reinit(cell, f); + + const unsigned int n_face_q_points = + fe_face_values.get_quadrature().size(); + + for (unsigned int j = 0; j < dofs_per_cell; ++j) + { + if (!discretization->finite_element.has_support_on_face(j, f)) + continue; + + Tensor<1, dim> normal; + if (id == Boundary::slip) + { + for (unsigned int q = 0; q < n_face_q_points; ++q) + normal += fe_face_values.normal_vector(q) * + fe_face_values.shape_value(j, q); + } + + const auto index = local_dof_indices[j]; + + Point position; + const auto global_index = partitioner->local_to_global(index); + for (unsigned int v = 0; + v < GeometryInfo::vertices_per_cell; + ++v) + if (cell->vertex_dof_index(v, 0) == global_index) + position = cell->vertex(v); + + const auto old_id = + std::get<1>(local_boundary_normal_map[index]); + local_boundary_normal_map[index] = + std::make_tuple(normal, std::max(old_id, id), position); + } /* j */ + } /* f */ + }; /* done with the definition of the worker */ + + /* This is the copy data routine for WorkStream */ + const auto copy_local_to_global = [&](const auto ©) { + const auto &is_artificial = copy.is_artificial; + const auto &local_dof_indices = copy.local_dof_indices; + const auto &local_boundary_normal_map = copy.local_boundary_normal_map; + const auto &cell_lumped_mass_matrix = copy.cell_lumped_mass_matrix; + const auto &cell_cij_matrix = copy.cell_cij_matrix; + + if (is_artificial) + return; + + for (const auto &it : local_boundary_normal_map) + { + auto &[normal, id, position] = boundary_normal_map[it.first]; + auto &[new_normal, new_id, new_position] = it.second; + + normal += new_normal; + id = std::max(id, new_id); + position = new_position; + } + + lumped_mass_matrix.add(local_dof_indices, cell_lumped_mass_matrix); + + for (int k = 0; k < dim; ++k) + { + cij_matrix[k].add(local_dof_indices, cell_cij_matrix[k]); + nij_matrix[k].add(local_dof_indices, cell_cij_matrix[k]); + } + }; /* end of the copy data routine */ + + WorkStream::run(dof_handler.begin_active(), + dof_handler.end(), + local_assemble_system, + copy_local_to_global, + scratch_data, + CopyData()); + } /* We are done with m_i and c_{ij} */ + + // At this point in time we are done with the computation of $m_i$ and + // $\mathbf{c}_{ij}$, but so far the matrix nij_matrix + // contains a just copy of the matrix cij_matrix. + // That's not what we really + // want: we have to normalize its entries. In addition, we have not even + // touched the entries of the matrix norm_matrix yet. We would + // like to exploit thread paralellization in order to carry out such + // operations, but WorkStream executes parallel cell-loops, so it might not + // the right tool. We want to execute node-loops: we + // want to visit every node $i$ in the mesh/sparsity graph, and for every + // such node we want to visit to every $j$ such that + // $\mathbf{c}_{ij} \not \equiv 0$. From an algebraic point of view, this is + // equivalent to: visiting every row in the matrix (equivalently sparsity + // pattern) and for each one of these rows execute a loop on the columns. + // Node-loops is a core theme of this tutorial step (see the pseudo-code + // in the introduction). + // + // We have the thread paralellization capability + // parallel::apply_to_subranges that is somehow more general than the + // WorkStream framework, an in particular it can be used for our node-loops. + // This functionality requires four input arguments: + // - A begin iterator: indices.begin() + // - A end iterator: indices.end() + // - A function f(i1,i2), where i1 and i2 define a + // sub-range with the range spanned by the the end and begin iterators + // of the previous two bullets. The function f(i1,i2) is + // called on_subranges in this example. It applies an + // operation for every "abstract element" in the subrange. In this case + // each "element" is a row rows of the sparsity pattern. + // - Grainsize: minimum number of "elements" (in this case rows) processed + // by + // each thread. We decided for a minimum of 4096 rows. + // + // We start by defining the operation on_subranges to be + // applied at each row in the sub-range. Given a fixed + // row_index we want to visit every entry in such row. In order + // to execute such columns-loops we use + // std::for_each + // from the standard library, where: + // sparsity_pattern.begin(row_index) + // gives us an iterator starting at the first column, + // sparsity_pattern.end(row_index) is an iterator pointing at + // the last column of the row. The last + // argument required by std::for_each is the operation applied at each + // column (a lambda expression in this case) of such row. We note that + // because of the nature of the data that we want to modify (we want to + // modify entries of a entire row at a time) threads cannot collide + // attempting to write the same entry (we do not need a scheduler). This + // advantage appears to be a particular characteristic of edge-based finite + // element schemes when they are properly implemented. + + // boost::irange + + { + TimerOutput::Scope t(computing_timer, + "offline_data - compute |c_ij|, and n_ij"); + + const auto on_subranges = [&](auto i1, const auto i2) { + for (; i1 < i2; ++i1) + { + const auto row_index = *i1; + + /* First column-loop: we compute/store the entries of the matrix + norm_matrix */ + std::for_each(sparsity_pattern.begin(row_index), + sparsity_pattern.end(row_index), + [&](const auto &jt) { + const auto value = + gather_get_entry(cij_matrix, &jt); + const double norm = value.norm(); + set_entry(norm_matrix, &jt, norm); + }); + + /* Second column-loop: we normalize the entries of the matrix + nij_matrix */ + for (auto &matrix : nij_matrix) + { + auto nij_entry = matrix.begin(row_index); + std::for_each(norm_matrix.begin(row_index), + norm_matrix.end(row_index), + [&](const auto &it) { + const auto norm = it.value(); + nij_entry->value() /= norm; + ++nij_entry; + }); + } + + } /* row_index */ + }; /* done with the definition of "on_subranges" */ + + const auto indices = boost::irange(0, n_locally_relevant); + parallel::apply_to_subranges(indices.begin(), + indices.end(), + on_subranges, + 4096); + + for (auto &it : boundary_normal_map) + { + auto &[normal, id, _] = it.second; + normal /= (normal.norm() + std::numeric_limits::epsilon()); + } + } + + // Placeholder here. + + { + TimerOutput::Scope t(computing_timer, + "offline_data - fix slip boundary c_ij"); + + const auto local_assemble_system = [&](const auto &cell, + auto & scratch, + auto & copy) { + auto &is_artificial = copy.is_artificial; + auto &local_dof_indices = copy.local_dof_indices; + + auto &cell_cij_matrix = copy.cell_cij_matrix; + + is_artificial = cell->is_artificial(); + if (is_artificial) + return; + + for (auto &matrix : cell_cij_matrix) + matrix.reinit(dofs_per_cell, dofs_per_cell); + + local_dof_indices.resize(dofs_per_cell); + cell->get_dof_indices(local_dof_indices); + std::transform(local_dof_indices.begin(), + local_dof_indices.end(), + local_dof_indices.begin(), + [&](auto index) { + return partitioner->global_to_local(index); + }); + + for (auto &matrix : cell_cij_matrix) + matrix = 0.; + + for (unsigned int f = 0; f < GeometryInfo::faces_per_cell; ++f) + { + const auto face = cell->face(f); + const auto id = face->boundary_id(); + + if (!face->at_boundary()) + continue; + + if (id != Boundary::slip) + continue; + + const auto &fe_face_values = scratch.reinit(cell, f); + + const unsigned int n_face_q_points = + fe_face_values.get_quadrature().size(); + + for (unsigned int q = 0; q < n_face_q_points; ++q) + { + const auto JxW = fe_face_values.JxW(q); + const auto normal_q = fe_face_values.normal_vector(q); + + for (unsigned int j = 0; j < dofs_per_cell; ++j) + { + if (!discretization->finite_element.has_support_on_face(j, + f)) + continue; + + const auto &[normal_j, _1, _2] = + boundary_normal_map[local_dof_indices[j]]; + + const auto value_JxW = + fe_face_values.shape_value(j, q) * JxW; + + for (unsigned int i = 0; i < dofs_per_cell; ++i) + { + const auto value = fe_face_values.shape_value(i, q); + + for (unsigned int d = 0; d < dim; ++d) + cell_cij_matrix[d](i, j) += + (normal_j[d] - normal_q[d]) * (value * value_JxW); + } /* i */ + } /* j */ + } /* q */ + } /* f */ + }; + + const auto copy_local_to_global = [&](const auto ©) { + const auto &is_artificial = copy.is_artificial; + const auto &local_dof_indices = copy.local_dof_indices; + const auto &cell_cij_matrix = copy.cell_cij_matrix; + + if (is_artificial) + return; + + for (int k = 0; k < dim; ++k) + cij_matrix[k].add(local_dof_indices, cell_cij_matrix[k]); + }; + + WorkStream::run(dof_handler.begin_active(), + dof_handler.end(), + local_assemble_system, + copy_local_to_global, + scratch_data, + CopyData()); + } + } /* assemble() */ + + // Placeholder here. + + template + DEAL_II_ALWAYS_INLINE inline dealii::Tensor<1, dim> + ProblemDescription::momentum(const rank1_type U) + { + dealii::Tensor<1, dim> result; + std::copy(&U[1], &U[1 + dim], &result[0]); + return result; + } + + // Placeholder here. + + template + DEAL_II_ALWAYS_INLINE inline double + ProblemDescription::internal_energy(const rank1_type U) + { + const double &rho = U[0]; + const auto m = momentum(U); + const double &E = U[dim + 1]; + return E - 0.5 * m.norm_square() / rho; + } + + // Placeholder here. + + template + DEAL_II_ALWAYS_INLINE inline double + ProblemDescription::pressure(const rank1_type U) + { + return (gamma - 1.) * internal_energy(U); + } + + // Placeholder here. + + + template + DEAL_II_ALWAYS_INLINE inline double + ProblemDescription::speed_of_sound(const rank1_type U) + { + const double &rho = U[0]; + const double p = pressure(U); + + return std::sqrt(gamma * p / rho); + } + + // Placeholder here. + + template + DEAL_II_ALWAYS_INLINE inline typename ProblemDescription::rank2_type + ProblemDescription::f(const rank1_type U) + { + const double &rho = U[0]; + const auto m = momentum(U); + const auto p = pressure(U); + const double &E = U[dim + 1]; + + rank2_type result; + + result[0] = m; + for (unsigned int i = 0; i < dim; ++i) + { + result[1 + i] = m * m[i] / rho; + result[1 + i][i] += p; + } + result[dim + 1] = m / rho * (E + p); + + return result; + } + + // Placeholder here. + + namespace + { + template + DEAL_II_ALWAYS_INLINE inline std::array riemann_data_from_state( + const typename ProblemDescription::rank1_type U, + const dealii::Tensor<1, dim> & n_ij) + { + dealii::Tensor<1, 3> projected_U; + projected_U[0] = U[0]; + + const auto m = ProblemDescription::momentum(U); + projected_U[1] = n_ij * m; + + const auto perpendicular_m = m - projected_U[1] * n_ij; + projected_U[2] = U[1 + dim] - 0.5 * perpendicular_m.norm_square() / U[0]; + + std::array result; + result[0] = projected_U[0]; + result[1] = projected_U[1] / projected_U[0]; + result[2] = ProblemDescription<1>::pressure(projected_U); + result[3] = ProblemDescription<1>::speed_of_sound(projected_U); + + return result; + } + + + DEAL_II_ALWAYS_INLINE inline double positive_part(const double number) + { + return (std::abs(number) + number) / 2.0; + } + + + DEAL_II_ALWAYS_INLINE inline double negative_part(const double number) + { + return (std::fabs(number) - number) / 2.0; + } + + + DEAL_II_ALWAYS_INLINE inline double + lambda1_minus(const std::array &riemann_data, + const double p_star) + { + constexpr double gamma = ProblemDescription<1>::gamma; + const auto &[rho_Z, u_Z, p_Z, a_Z] = riemann_data; + + const double factor = (gamma + 1.0) / 2.0 / gamma; + const double tmp = positive_part((p_star - p_Z) / p_Z); + return u_Z - a_Z * std::sqrt(1.0 + factor * tmp); + } + + + DEAL_II_ALWAYS_INLINE inline double + lambda3_plus(const std::array &riemann_data, const double p_star) + { + constexpr double gamma = ProblemDescription<1>::gamma; + const auto &[rho_Z, u_Z, p_Z, a_Z] = riemann_data; + + const double factor = (gamma + 1.0) / 2.0 / gamma; + const double tmp = positive_part((p_star - p_Z) / p_Z); + return u_Z + a_Z * std::sqrt(1.0 + factor * tmp); + } + + + DEAL_II_ALWAYS_INLINE inline double + lambda_max_two_rarefaction(const std::array &riemann_data_i, + const std::array &riemann_data_j) + { + constexpr double gamma = ProblemDescription<1>::gamma; + const auto &[rho_i, u_i, p_i, a_i] = riemann_data_i; + const auto &[rho_j, u_j, p_j, a_j] = riemann_data_j; + + const double numerator = a_i + a_j - (gamma - 1.) / 2. * (u_j - u_i); + + const double denominator = + a_i * std::pow(p_i / p_j, -1. * (gamma - 1.) / 2. / gamma) + a_j * 1.; + + const double p_star = + p_j * std::pow(numerator / denominator, 2. * gamma / (gamma - 1)); + + const double lambda1 = lambda1_minus(riemann_data_i, p_star); + const double lambda3 = lambda3_plus(riemann_data_j, p_star); + + return std::max(positive_part(lambda3), negative_part(lambda1)); + }; + + + DEAL_II_ALWAYS_INLINE inline double + lambda_max_expansion(const std::array &riemann_data_i, + const std::array &riemann_data_j) + { + const auto &[rho_i, u_i, p_i, a_i] = riemann_data_i; + const auto &[rho_j, u_j, p_j, a_j] = riemann_data_j; + + return std::max(std::abs(u_i), std::abs(u_j)) + 5. * std::max(a_i, a_j); + } + } // namespace + + // Placeholder here. + + template + DEAL_II_ALWAYS_INLINE inline double + ProblemDescription::compute_lambda_max( + const rank1_type U_i, + const rank1_type U_j, + const dealii::Tensor<1, dim> &n_ij) + { + const auto riemann_data_i = riemann_data_from_state(U_i, n_ij); + const auto riemann_data_j = riemann_data_from_state(U_j, n_ij); + + const double lambda_1 = + lambda_max_two_rarefaction(riemann_data_i, riemann_data_j); + + const double lambda_2 = + lambda_max_expansion(riemann_data_i, riemann_data_j); + + return std::min(lambda_1, lambda_2); + } + + // Placeholder here. + + template <> + const std::array // + ProblemDescription<1>::component_names{"rho", "m", "E"}; + + template <> + const std::array // + ProblemDescription<2>::component_names{"rho", "m_1", "m_2", "E"}; + + template <> + const std::array // + ProblemDescription<3>::component_names{"rho", "m_1", "m_2", "m_3", "E"}; + + // Placeholder here. + + template + InitialValues::InitialValues(const std::string &subsection) + : ParameterAcceptor(subsection) + { + ParameterAcceptor::parse_parameters_call_back.connect( + std::bind(&InitialValues::parse_parameters_callback, this)); + + initial_direction[0] = 1.; + add_parameter("initial direction", + initial_direction, + "Initial direction of the uniform flow field"); + + static constexpr auto gamma = ProblemDescription::gamma; + initial_1d_state[0] = gamma; + initial_1d_state[1] = 3.; + initial_1d_state[2] = 1.; + add_parameter("initial 1d state", + initial_1d_state, + "Initial 1d state (rho, u, p) of the uniform flow field"); + } + + // Placeholder here. + + template + void InitialValues::parse_parameters_callback() + { + AssertThrow(initial_direction.norm() != 0., + ExcMessage( + "Initial shock front direction is set to the zero vector.")); + initial_direction /= initial_direction.norm(); + + static constexpr auto gamma = ProblemDescription::gamma; + + const auto from_1d_state = + [=](const dealii::Tensor<1, 3, double> &state_1d) -> rank1_type { + const auto &rho = state_1d[0]; + const auto &u = state_1d[1]; + const auto &p = state_1d[2]; + + rank1_type state; + + state[0] = rho; + for (unsigned int i = 0; i < dim; ++i) + state[1 + i] = rho * u * initial_direction[i]; + state[dim + 1] = p / (gamma - 1.) + 0.5 * rho * u * u; + + return state; + }; + + initial_state = [=](const dealii::Point & /*point*/, double /*t*/) { + return from_1d_state(initial_1d_state); + }; + } + + // Placeholder here. + + template + TimeStep::TimeStep(const MPI_Comm & mpi_communicator, + dealii::TimerOutput & computing_timer, + const OfflineData & offline_data, + const InitialValues &initial_values, + const std::string & subsection /*= "TimeStep"*/) + : ParameterAcceptor(subsection) + , mpi_communicator(mpi_communicator) + , computing_timer(computing_timer) + , offline_data(&offline_data) + , initial_values(&initial_values) + { + cfl_update = 0.80; + add_parameter("cfl update", + cfl_update, + "relative CFL constant used for update"); + } + + // Placeholder here. + + template + void TimeStep::prepare() + { + TimerOutput::Scope time(computing_timer, + "time_step - prepare scratch space"); + + const auto &partitioner = offline_data->partitioner; + for (auto &it : temp) + it.reinit(partitioner); + + const auto &sparsity = offline_data->sparsity_pattern; + dij_matrix.reinit(sparsity); + } + + // Placeholder here. + + template + double TimeStep::step(vector_type &U, double t) + { + const auto &n_locally_owned = offline_data->n_locally_owned; + const auto &n_locally_relevant = offline_data->n_locally_relevant; + + const auto indices_owned = boost::irange(0, n_locally_owned); + const auto indices_relevant = + boost::irange(0, n_locally_relevant); + + const auto &sparsity = offline_data->sparsity_pattern; + + const auto &lumped_mass_matrix = offline_data->lumped_mass_matrix; + const auto &norm_matrix = offline_data->norm_matrix; + const auto &nij_matrix = offline_data->nij_matrix; + const auto &cij_matrix = offline_data->cij_matrix; + + const auto &boundary_normal_map = offline_data->boundary_normal_map; + + { + TimerOutput::Scope time(computing_timer, "time_step - 1 compute d_ij"); + + const auto on_subranges = [&](auto i1, const auto i2) { + for (const auto i : boost::make_iterator_range(i1, i2)) + { + const auto U_i = gather(U, i); + + /* Column-loop */ + for (auto jt = sparsity.begin(i); jt != sparsity.end(i); ++jt) + { + const auto j = jt->column(); + + if (j >= i) + continue; + + const auto U_j = gather(U, j); + + const auto n_ij = gather_get_entry(nij_matrix, jt); + const double norm = get_entry(norm_matrix, jt); + + const auto lambda_max = + ProblemDescription::compute_lambda_max(U_i, U_j, n_ij); + + double d = norm * lambda_max; + + if (boundary_normal_map.count(i) != 0 && + boundary_normal_map.count(j) != 0) + { + const auto n_ji = gather(nij_matrix, j, i); + const auto lambda_max_2 = + ProblemDescription::compute_lambda_max(U_j, + U_i, + n_ji); + const double norm_2 = norm_matrix(j, i); + + d = std::max(d, norm_2 * lambda_max_2); + } + + set_entry(dij_matrix, jt, d); + dij_matrix(j, i) = d; + } /* End of column-loop */ + } /* End of row-loop */ + }; /* End of definition of on_subranges */ + + parallel::apply_to_subranges(indices_relevant.begin(), + indices_relevant.end(), + on_subranges, + 4096); + } /* End of the computation of the off-diagonal entries of dij_matrix */ + + std::atomic tau_max{std::numeric_limits::infinity()}; + + { + TimerOutput::Scope time(computing_timer, + "time_step - 2 compute d_ii, and tau_max"); + + const auto on_subranges = [&](auto i1, const auto i2) { + double tau_max_on_subrange = std::numeric_limits::infinity(); + + for (const auto i : boost::make_iterator_range(i1, i2)) + { + double d_sum = 0.; + + for (auto jt = sparsity.begin(i); jt != sparsity.end(i); ++jt) + { + const auto j = jt->column(); + + if (j == i) + continue; + + d_sum -= get_entry(dij_matrix, jt); + } + + dij_matrix.diag_element(i) = d_sum; + + const double mass = lumped_mass_matrix.diag_element(i); + const double tau = cfl_update * mass / (-2. * d_sum); + tau_max_on_subrange = std::min(tau_max_on_subrange, tau); + } + + double current_tau_max = tau_max.load(); + while ( + current_tau_max > tau_max_on_subrange && + !tau_max.compare_exchange_weak(current_tau_max, tau_max_on_subrange)) + ; + }; + + parallel::apply_to_subranges(indices_relevant.begin(), + indices_relevant.end(), + on_subranges, + 4096); + + tau_max.store(Utilities::MPI::min(tau_max.load(), mpi_communicator)); + + AssertThrow(!std::isnan(tau_max) && !std::isinf(tau_max) && tau_max > 0., + ExcMessage("I'm sorry, Dave. I'm afraid I can't " + "do that. - We crashed.")); + } /* End of the computation of the diagonal entries of dij_matrix */ + + { + TimerOutput::Scope time(computing_timer, "time_step - 3 perform update"); + + const auto on_subranges = [&](auto i1, const auto i2) { + for (const auto i : boost::make_iterator_range(i1, i2)) + { + Assert(i < n_locally_owned, ExcInternalError()); + + const auto U_i = gather(U, i); + + const auto f_i = ProblemDescription::f(U_i); + const double m_i = lumped_mass_matrix.diag_element(i); + + auto U_i_new = U_i; + + for (auto jt = sparsity.begin(i); jt != sparsity.end(i); ++jt) + { + const auto j = jt->column(); + + const auto U_j = gather(U, j); + const auto f_j = ProblemDescription::f(U_j); + + const auto c_ij = gather_get_entry(cij_matrix, jt); + const auto d_ij = get_entry(dij_matrix, jt); + + for (unsigned int k = 0; k < problem_dimension; ++k) + { + U_i_new[k] += + tau_max / m_i * + (-(f_j[k] - f_i[k]) * c_ij + d_ij * (U_j[k] - U_i[k])); + } + } + + scatter(temp, U_i_new, i); + } + }; + + /* Only iterate over locally owned subset! */ + parallel::apply_to_subranges(indices_owned.begin(), + indices_owned.end(), + on_subranges, + 4096); + } /* End of the computation of the new solution */ + + { + TimerOutput::Scope time(computing_timer, + "time_step - 4 fix boundary states"); + + const auto on_subranges = [&](const auto it1, const auto it2) { + for (auto it = it1; it != it2; ++it) + { + const auto i = it->first; + + /* Only iterate over locally owned subset */ + if (i >= n_locally_owned) + continue; + + const auto &[normal, id, position] = it->second; + + /* Skip constrained degrees of freedom */ + if (++sparsity.begin(i) == sparsity.end(i)) + continue; + + auto U_i = gather(temp, i); + + /* On boundary 1 remove the normal component of the momentum: */ + + if (id == Boundary::slip) + { + auto m = ProblemDescription::momentum(U_i); + m -= 1. * (m * normal) * normal; + for (unsigned int k = 0; k < dim; ++k) + U_i[k + 1] = m[k]; + } + + /* On boundary 2 enforce initial conditions: */ + + if (id == Boundary::dirichlet) + { + U_i = initial_values->initial_state(position, t + tau_max); + } + + scatter(temp, U_i, i); + } + }; + + on_subranges(boundary_normal_map.begin(), boundary_normal_map.end()); + } + + for (auto &it : temp) + it.update_ghost_values(); + + U.swap(temp); + + return tau_max; + } /* End of TimeStep::step */ + + // Placeholder here. + + template + SchlierenPostprocessor::SchlierenPostprocessor( + const MPI_Comm & mpi_communicator, + dealii::TimerOutput & computing_timer, + const OfflineData &offline_data, + const std::string & subsection /*= "SchlierenPostprocessor"*/) + : ParameterAcceptor(subsection) + , mpi_communicator(mpi_communicator) + , computing_timer(computing_timer) + , offline_data(&offline_data) + { + schlieren_beta = 10.; + add_parameter("schlieren beta", + schlieren_beta, + "Beta factor used in Schlieren-type postprocessor"); + + schlieren_index = 0; + add_parameter("schlieren index", + schlieren_index, + "Use the corresponding component of the state vector for the " + "schlieren plot"); + } + + template + void SchlierenPostprocessor::prepare() + { + TimerOutput::Scope t(computing_timer, + "schlieren_postprocessor - prepare scratch space"); + + const auto &n_locally_relevant = offline_data->n_locally_relevant; + const auto &partitioner = offline_data->partitioner; + + r.reinit(n_locally_relevant); + schlieren.reinit(partitioner); + } + + template + void SchlierenPostprocessor::compute_schlieren(const vector_type &U) + { + TimerOutput::Scope t(computing_timer, + "schlieren_postprocessor - compute schlieren plot"); + + const auto &sparsity = offline_data->sparsity_pattern; + const auto &lumped_mass_matrix = offline_data->lumped_mass_matrix; + const auto &cij_matrix = offline_data->cij_matrix; + const auto &boundary_normal_map = offline_data->boundary_normal_map; + + const auto &n_locally_owned = offline_data->n_locally_owned; + const auto indices = boost::irange(0, n_locally_owned); + + std::atomic r_i_max{0.}; + std::atomic r_i_min{std::numeric_limits::infinity()}; + + { + const auto on_subranges = [&](auto i1, const auto i2) { + double r_i_max_on_subrange = 0.; + double r_i_min_on_subrange = std::numeric_limits::infinity(); + + for (; i1 < i2; ++i1) + { + const auto i = *i1; + + Assert(i < n_locally_owned, ExcInternalError()); + + Tensor<1, dim> r_i; + + for (auto jt = sparsity.begin(i); jt != sparsity.end(i); ++jt) + { + const auto j = jt->column(); + + if (i == j) + continue; + + const auto U_js = U[schlieren_index].local_element(j); + const auto c_ij = gather_get_entry(cij_matrix, jt); + + r_i += c_ij * U_js; + } + + const auto bnm_it = boundary_normal_map.find(i); + if (bnm_it != boundary_normal_map.end()) + { + const auto [normal, id, _] = bnm_it->second; + if (id == Boundary::slip) + { + r_i -= 1. * (r_i * normal) * normal; + } + else + { + r_i = 0.; + } + } + + const double m_i = lumped_mass_matrix.diag_element(i); + r[i] = r_i.norm() / m_i; + + r_i_max_on_subrange = std::max(r_i_max_on_subrange, r[i]); + r_i_min_on_subrange = std::min(r_i_min_on_subrange, r[i]); + } + + double current_r_i_max = r_i_max.load(); + while ( + current_r_i_max < r_i_max_on_subrange && + !r_i_max.compare_exchange_weak(current_r_i_max, r_i_max_on_subrange)) + ; + + double current_r_i_min = r_i_min.load(); + while ( + current_r_i_min > r_i_min_on_subrange && + !r_i_min.compare_exchange_weak(current_r_i_min, r_i_min_on_subrange)) + ; + }; + + parallel::apply_to_subranges(indices.begin(), + indices.end(), + on_subranges, + 4096); + } + + r_i_max.store(Utilities::MPI::max(r_i_max.load(), mpi_communicator)); + r_i_min.store(Utilities::MPI::min(r_i_min.load(), mpi_communicator)); + + { + const auto on_subranges = [&](auto i1, const auto i2) { + for (; i1 < i2; ++i1) + { + const auto i = *i1; + + Assert(i < n_locally_owned, ExcInternalError()); + + schlieren.local_element(i) = + 1. - std::exp(-schlieren_beta * (r[i] - r_i_min) / + (r_i_max - r_i_min)); + } + }; + + parallel::apply_to_subranges(indices.begin(), + indices.end(), + on_subranges, + 4096); + } + + schlieren.update_ghost_values(); + } + + // Placeholder here. + + template + TimeLoop::TimeLoop(const MPI_Comm &mpi_comm) + : ParameterAcceptor("A - TimeLoop") + , mpi_communicator(mpi_comm) + , computing_timer(mpi_communicator, + timer_output, + TimerOutput::never, + TimerOutput::cpu_and_wall_times) + , pcout(std::cout, Utilities::MPI::this_mpi_process(mpi_communicator) == 0) + , discretization(mpi_communicator, computing_timer, "B - Discretization") + , offline_data(mpi_communicator, + computing_timer, + discretization, + "C - OfflineData") + , initial_values("D - InitialValues") + , time_step(mpi_communicator, + computing_timer, + offline_data, + initial_values, + "E - TimeStep") + , schlieren_postprocessor(mpi_communicator, + computing_timer, + offline_data, + "F - SchlierenPostprocessor") + { + base_name = "test"; + add_parameter("basename", base_name, "Base name for all output files"); + + t_final = 4.; + add_parameter("final time", t_final, "Final time"); + + output_granularity = 0.02; + add_parameter("output granularity", + output_granularity, + "time interval for output"); + + resume = false; + add_parameter("resume", resume, "Resume an interrupted computation."); + } + + // Placeholder here. + + namespace + { + void print_head(dealii::ConditionalOStream &pcout, + std::string header, + std::string secondary = "") + { + const auto header_size = header.size(); + const auto padded_header = std::string((34 - header_size) / 2, ' ') + + header + + std::string((35 - header_size) / 2, ' '); + + const auto secondary_size = secondary.size(); + const auto padded_secondary = + std::string((34 - secondary_size) / 2, ' ') + secondary + + std::string((35 - secondary_size) / 2, ' '); + + /* clang-format off */ + pcout << std::endl; + pcout << " ####################################################" << std::endl; + pcout << " ######### #########" << std::endl; + pcout << " #########" << padded_header << "#########" << std::endl; + pcout << " #########" << padded_secondary << "#########" << std::endl; + pcout << " ######### #########" << std::endl; + pcout << " ####################################################" << std::endl; + pcout << std::endl; + /* clang-format on */ + } + } // namespace + + + // Implementation of the class member interpolate_initial_values + // . + + template + void TimeLoop::run() + { + if (Utilities::MPI::this_mpi_process(mpi_communicator) == 0) + { + pcout << "Reading parameters and allocating objects... " << std::flush; + ParameterAcceptor::initialize("step-69.prm"); + pcout << "done" << std::endl; + } + else + { + ParameterAcceptor::initialize("step-69.prm"); + } + + print_head(pcout, "create triangulation"); + discretization.setup(); + + print_head(pcout, "compute offline data"); + offline_data.setup(); + offline_data.assemble(); + + print_head(pcout, "set up time step"); + time_step.prepare(); + schlieren_postprocessor.prepare(); + + double t = 0.; + unsigned int output_cycle = 0; + + print_head(pcout, "interpolate initial values"); + auto U = interpolate_initial_values(); + + if (resume) + { + print_head(pcout, "restore interrupted computation"); + + const auto & triangulation = discretization.triangulation; + const unsigned int i = triangulation.locally_owned_subdomain(); + std::string name = base_name + "-checkpoint-" + + dealii::Utilities::int_to_string(i, 4) + ".archive"; + std::ifstream file(name, std::ios::binary); + + boost::archive::binary_iarchive ia(file); + ia >> t >> output_cycle; + + for (auto &it1 : U) + { + for (auto &it2 : it1) + ia >> it2; + it1.update_ghost_values(); + } + } + + output(U, base_name + "-solution", t, output_cycle++); + + print_head(pcout, "enter main loop"); + + for (unsigned int cycle = 1; t < t_final; ++cycle) + { + std::ostringstream head; + head << "Cycle " << Utilities::int_to_string(cycle, 6) << " (" + << std::fixed << std::setprecision(1) << t / t_final * 100 << "%)"; + std::ostringstream secondary; + secondary << "at time t = " << std::setprecision(8) << std::fixed << t; + print_head(pcout, head.str(), secondary.str()); + + t += time_step.step(U, t); + + if (t > output_cycle * output_granularity) + output(U, base_name + "-solution", t, output_cycle++, true); + + } /* end of loop */ + + if (output_thread.joinable()) + output_thread.join(); + + computing_timer.print_summary(); + pcout << timer_output.str() << std::endl; + } + + // Implementation of the class member interpolate_initial_values + // . + + template + typename TimeLoop::vector_type + TimeLoop::interpolate_initial_values(double t) + { + pcout << "TimeLoop::interpolate_initial_values(t = " << t << ")" + << std::endl; + TimerOutput::Scope timer(computing_timer, + "time_loop - setup scratch space"); + + vector_type U; + + const auto &partitioner = offline_data.partitioner; + for (auto &it : U) + it.reinit(partitioner); + + constexpr auto problem_dimension = + ProblemDescription::problem_dimension; + + for (unsigned int i = 0; i < problem_dimension; ++i) + VectorTools::interpolate(offline_data.dof_handler, + ScalarFunctionFromFunctionObject( + [&](const auto &p) { + return initial_values.initial_state(p, t)[i]; + }), + U[i]); + + for (auto &it : U) + it.update_ghost_values(); + + return U; + } + + // Implementation of the class member TimeLoop . + + template + void TimeLoop::output(const typename TimeLoop::vector_type &U, + const std::string & name, + double t, + unsigned int cycle, + bool checkpoint) + { + pcout << "TimeLoop::output(t = " << t + << ", checkpoint = " << checkpoint << ")" << std::endl; + + if (output_thread.joinable()) + { + TimerOutput::Scope timer(computing_timer, "time_loop - stalled output"); + output_thread.join(); + } + + constexpr auto problem_dimension = + ProblemDescription::problem_dimension; + const auto &component_names = ProblemDescription::component_names; + + for (unsigned int i = 0; i < problem_dimension; ++i) + { + output_vector[i] = U[i]; + output_vector[i].update_ghost_values(); + } + + schlieren_postprocessor.compute_schlieren(output_vector); + + const auto output_worker = [this, name, t, cycle, checkpoint]() { + constexpr auto problem_dimension = + ProblemDescription::problem_dimension; + const auto &dof_handler = offline_data.dof_handler; + const auto &triangulation = discretization.triangulation; + const auto &mapping = discretization.mapping; + + if (checkpoint) + { + const unsigned int i = triangulation.locally_owned_subdomain(); + std::string name = base_name + "-checkpoint-" + + dealii::Utilities::int_to_string(i, 4) + + ".archive"; + + if (std::filesystem::exists(name)) + std::filesystem::rename(name, name + "~"); + + std::ofstream file(name, std::ios::binary | std::ios::trunc); + + boost::archive::binary_oarchive oa(file); + oa << t << cycle; + for (const auto &it1 : output_vector) + for (const auto &it2 : it1) + oa << it2; + } + + dealii::DataOut data_out; + data_out.attach_dof_handler(dof_handler); + + for (unsigned int i = 0; i < problem_dimension; ++i) + data_out.add_data_vector(output_vector[i], component_names[i]); + + data_out.add_data_vector(schlieren_postprocessor.schlieren, + "schlieren_plot"); + + data_out.build_patches(mapping, discretization.finite_element.degree - 1); + + DataOutBase::VtkFlags flags(t, + cycle, + true, + DataOutBase::VtkFlags::best_speed); + data_out.set_flags(flags); + + const auto filename = [&](const unsigned int i) -> std::string { + const auto seq = dealii::Utilities::int_to_string(i, 4); + return name + "-" + Utilities::int_to_string(cycle, 6) + "-" + seq + + ".vtu"; + }; + + const unsigned int i = triangulation.locally_owned_subdomain(); + std::ofstream output(filename(i)); + data_out.write_vtu(output); + + if (dealii::Utilities::MPI::this_mpi_process(mpi_communicator) == 0) + { + const unsigned int n_mpi_processes = + dealii::Utilities::MPI::n_mpi_processes(mpi_communicator); + std::vector filenames; + for (unsigned int i = 0; i < n_mpi_processes; ++i) + filenames.push_back(filename(i)); + + std::ofstream output(name + "-" + Utilities::int_to_string(cycle, 6) + + ".pvtu"); + data_out.write_pvtu_record(output, filenames); + } + }; + + output_thread = std::move(std::thread(output_worker)); + } + +} // namespace Step69 + + +int main(int argc, char *argv[]) +{ + constexpr int dim = 2; + + using namespace dealii; + using namespace Step69; + + Utilities::MPI::MPI_InitFinalize mpi_initialization(argc, argv); + + MPI_Comm mpi_communicator(MPI_COMM_WORLD); + TimeLoop time_loop(mpi_communicator); + + time_loop.run(); +} -- 2.39.5