From: Sam Scheuerman Date: Thu, 12 Sep 2024 16:57:52 +0000 (-0600) Subject: Add the step-93 tutorial program. X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=799ef37c398eb8b057e7cd3087d0022b789ba737;p=dealii.git Add the step-93 tutorial program. --- diff --git a/doc/doxygen/tutorial/tutorial.h.in b/doc/doxygen/tutorial/tutorial.h.in index ee6a84be73..9acd8df457 100644 --- a/doc/doxygen/tutorial/tutorial.h.in +++ b/doc/doxygen/tutorial/tutorial.h.in @@ -700,6 +700,12 @@ * Solving the Laplace-Beltrami equation on a surface using the trace finite element method. *
Keywords: MeshWorker::mesh_loop(), NonMatching::FEImmersedSurfaceValues * + * + * + * step-93 + * Using nonlocal dofs to solve a simple optimization problem with the Poisson equation + *
Keywords: nonlocal dofs, optimization, Poisson + * * * * diff --git a/doc/news/changes/major/20241106Scheuerman b/doc/news/changes/major/20241106Scheuerman new file mode 100644 index 0000000000..a89ab2880d --- /dev/null +++ b/doc/news/changes/major/20241106Scheuerman @@ -0,0 +1,4 @@ +New: The step-93 tutorial program shows how to use nonlocal dofs in the +deal.II framework, in the context of a simple optimization problem. +
+(Sam Scheuerman, Wolfgang Bangerth, 2024/11/06) diff --git a/examples/step-93/CMakeLists.txt b/examples/step-93/CMakeLists.txt new file mode 100644 index 0000000000..af81e345d2 --- /dev/null +++ b/examples/step-93/CMakeLists.txt @@ -0,0 +1,39 @@ +## +# CMake script for the step-93 tutorial program: +## + +# Set the name of the project and target: +set(TARGET "step-93") + +# 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 3.13.4) + +find_package(deal.II 9.7.0 + HINTS ${deal.II_DIR} ${DEAL_II_DIR} ../ ../../ ../../../ ../../build-dealii $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() + +deal_ii_initialize_cached_variables() +project(${TARGET}) +deal_ii_invoke_autopilot() diff --git a/examples/step-93/doc/builds-on b/examples/step-93/doc/builds-on new file mode 100644 index 0000000000..11f3e73b22 --- /dev/null +++ b/examples/step-93/doc/builds-on @@ -0,0 +1 @@ +step-4 step-8 diff --git a/examples/step-93/doc/intro.dox b/examples/step-93/doc/intro.dox new file mode 100644 index 0000000000..03a121ced6 --- /dev/null +++ b/examples/step-93/doc/intro.dox @@ -0,0 +1,370 @@ +This program was contributed by Sam Scheuerman and Wolfgang +Bangerth. +
+ + +

Introduction

+ +In all previous tutorial programs, we have only ever dealt with +discretizations of partial differential equations posed on a +Triangulation where every degree of freedom was associated with +a local geometric element (vertices, lines, quadrangles, etc.) which +is clearly part of one cell or another (or, in the worst case, some +small number of adjacent cells). This was true +whether we were solving a single differential equation (say, in step-4 +or step-6), or a system of coupled equations (step-8, step-22, and +many others). In other words, the *unknowns* we were considering were +always functions of space (and sometimes time) that we discretized on +a Triangulation. + +There are also situations where we have unknowns in the equations +we consider that are *not* functions but are *scalars* instead. The +example we will consider here to illustrate this kind of situation is +an optimization problem: Let's assume we have a body to which we have +attached four heaters. How should we set the power levels $C_1,\ldots, +C_4$ of these four heaters so that the (steady state) temperature $u(\mathbf x)$ of +the body is as close as possible to a desired state? The issue here is +that we have one unknown that is a function (namely, the temperature +$u(\mathbf x)$) for which we can apply the usual finite element +approximation, but we also have four scalars $C_i$ that are *not* +functions and that, after discretization of the whole problem, are not +tied to one of the cells of the mesh. As a consequence, we call these +unknowns *non-local degrees of freedom* to represent that they are not +tied to a specific locality: They represent something that does not +depend to the spatial variable $\mathbf x$ at all, and is not the +result of a discretization process. + +Before going into how we solve this issue, let us first state the concrete set of equations. + +We want to find a solution $u$ to the Poisson equation + +@f{align*}{ + -\Delta u =& f \text{ in }\Omega\\ + u =&0 \text{ on } \partial\Omega, +@f} + +which is as close as +possible to a desired solution $\overline{u}$. We can formulate this as an optimization problem + +@f{align*}{ + \min_{u, f} & \frac 12 ||u-\overline{u}||^2\\ + \text{s.t.} & -\Delta u = f. +@f} + +Instead of allowing *any* source function $f$, we will assume that the heat source is described by four +heating pads $f_k(\mathbf x)$ that are characteristic functions of a subset of the domain, times +a heater setting $C^k$ (perhaps the power level you set each heating pad to). +In other words, we have that $f(\mathbf x) = \sum_{k=0}^m C^k f_k(\mathbf x)$. +The optimization problem then reads as follows: +@f{align*}{ + \min_{u, C^k} & \frac 12 ||u-\overline{u}||^2\\ + \text{s.t.} & -\Delta u = \sum_k C^k f_k. +@f} + +To write this optimization problem in a more compact form, we introduce a Lagrangian + +@f{align*}{ + L = \frac{1}{2}||u - \overline{u}||^2 - \left(\lambda, \left(-\Delta u - \sum_k C^k f_k\right)\right), +@f} + +where $\lambda = \lambda(\mathbf{x})$ is a Lagrange multiplier and $(\cdot, \cdot)$ denotes the usual $L^2$ scalar +product (note that $\lambda$ is a function that lives in the space of test functions, not just a scalar). We do this because finding a critical point of $L$ +is equivalent to solving the original optimization problem, and we can find a critical point +of $L$ by finding roots of its derivatives. At this point, we have a choice: we can either +derive the optimality conditions first, then discretize the problem (in the literature this +is called the optimize then discretize approach), or we can first discretize $L$ and then +determine the optimality conditions for the discrete system (discretize then optimize). We +will use the second approach here, because the formulation is a more straightforward. +The discretized Lagrangian is + +@f{align*}{ + L =& \frac{1}{2}||u_h - \overline{u}||^2 - (\lambda_h, (-\Delta u_h - f))\\ + =& \frac{1}{2}\left[\sum_{i=0}^n\sum_{j=0}^n U^iU^j\int_{\Omega}\varphi_i\varphi_j + - 2\sum_{i=0}^n U^i\int_\Omega \varphi_i \overline{u} + \int_\Omega\overline{u}^2\right] + - \sum_{i=0}^n\sum_{j=0}^n\left[\Lambda^iU^j\int_{\Omega}\nabla\varphi_i\cdot\nabla\varphi_j\right] + + \sum_{i=0}^n \sum_{k=0}^m \left[\Lambda^i C^k \int_{\Omega}\varphi_i f_k\right]. +@f} + +where + +@f{align*}{ + u_h(x) =& \sum_i^n U^i\varphi_i(x),\\ + \lambda_h =& \sum_j^n \Lambda^j\varphi_j,\text{ and}\\ + f =& \sum_{k=0}^m C^k f_k. +@f} + +Note the application of integration by parts above, which makes use of the Dirichlet boundary conditions to eliminate the boundary integral. In practice, we assume that each $f_k$ is something simple. In this example step, each $f_k$ is a characteristic function as mentioned, and $m=4$. These functions, plotted together, look like this: + +
Visualization of four characteristic functions together on the same plot
+ + +_Note:_ At this point, we can clearly see where the nonlocal degrees of freedom are coming from. The term $\int_\Omega \varphi_i f_k$ is an integral which cannot a priori be confined to a single cell or small collection of cells, but rather must be evaluated on all cells that lie in the domain of $f$. + +To find a critical point, we take a partial derivative with respect to each $U^i$, + $\Lambda^j$, and $C^k$. The derivatives have the following form: + +@f{align*}{ + \frac{\partial L}{\partial U^i} =& \sum_{j=0}^nU^j\int_{\Omega}\varphi_i\varphi_j + - \int_\Omega \varphi_i\overline{u} + - \sum_{j=0}^n\left[\Lambda^j\int_{\Omega}\nabla\varphi_i\cdot\nabla\varphi_j\right]\\ + \frac{\partial L}{\partial \Lambda^j} =& -\sum_{i=0}^n\left[U^i\int_{\Omega}\nabla\varphi_i\cdot\nabla\varphi_j\right] + + \sum_{k=0}^m C^k\int_{\Omega}\varphi_j f_k\\ + \frac{\partial L}{\partial C^k} =& \sum_{j=0}^n \Lambda^j\int_{\Omega}\varphi_j f_k. +@f} + +Note the paired indices on the left- and right-hand sides. As with many other minimization problems, we can observe that for a critical point to occur, +these derivatives must all be 0. This gives us $n + n + m$ linear equations to solve +for the coefficient vectors $U$, $\Lambda$, and $C$: + +@f{align*}{ + \frac{\partial L}{\partial U^i} =& \sum_{j=0}^nU^j\int_{\Omega}\varphi_i\varphi_j + - \int_\Omega\varphi_i\overline{u} + - \sum_{j=0}^n\left[\Lambda^j\int_{\Omega}\nabla\varphi_i\cdot\nabla\varphi_j\right] = 0\\ + \frac{\partial L}{\partial \Lambda^j} =& -\sum_{i=0}^n\left[U^i\int_{\Omega}\nabla\varphi_i\cdot\nabla\varphi_j\right] + + \sum_{k=0}^m C^k\int_{\Omega}\varphi_j f_k = 0\\ + \frac{\partial L}{\partial C^k} =& \sum_{j=0}^n \Lambda^j\int_{\Omega}\varphi_j f_k = 0.\\ +@f} + +The summations are suggestive of matrix-vector multiplication, and in fact we can write the equations +in such a form. Take $(\mathcal{M}_{i, j}) = \int_\Omega\varphi_i\varphi_j$, $(\mathcal{N}_{i, j}) + = \int_\Omega \nabla\varphi_i\cdot\nabla\varphi_j$, $(\overline{\mathcal{U}}_i) = \int_\Omega\varphi_i\overline{u}$, + and $(\mathcal{F}_{j, k}) = \int_\Omega\varphi_j + f_k$, and note that $\mathcal{F}$ is an $n\times m$ matrix. Then we get the matrix system + + @f{align*}{ + \mathcal{M}U - \mathcal{U} - \mathcal{N}^T\Lambda =& 0\\ + -\mathcal{N}U + \mathcal{F}^TC =& 0\\ + \mathcal{F}\Lambda =& 0. + @f} + +Moving known quantities to the right-hand-side and combinging the matrix equations gives us an +idealized system to solve: + +@f{align*}{ + \left(\begin{array}{c c c} + \mathcal{M} & -\mathcal{N}^T & 0\\ + -\mathcal{N} & 0 & \mathcal{F}^T\\ + 0 & \mathcal{F} & 0 + \end{array}\right) + \left(\begin{array}{c} + U\\\Lambda\\ C + \end{array}\right) = + \left(\begin{array}{c} + \overline{\mathcal{U}}\\0\\0 + \end{array}\right). +@f} + +Unfortunately, setting up this type of block matrix, where some of the dofs couple in a +nonlocal way, does not work well with deal.II's internal structure. So, we will need to +modify how this problem is constructed to work well with deal.II, which is discussed below. + +

What is the problem?

+ +The straightforward way to deal with the equations above is to define +an FESystem object with two components, one for the temperature +$u(\mathbf x)$ and one for the Lagrange multiplier $\lambda(\mathbf +x)$. This would work in much the same way as you define system +elements in all other @ref vector_valued "vector valued problems". +You will then end up with $N_U$ variables for the finite element +approximation of the temperature, and $N_\Lambda$ variables for the +finite element approximation of the Lagrange multiplier, and all of +these $N_U+N_\Lambda$ unknowns live on the triangulation in the form +of a DoFHandler object. While in principle $N_U$ and $N_\Lambda$ can +be different (if you chose different element types, or different polynomial +degrees to discretize $u$ and $\lambda$), in practice it is best to keep them the same, so we have +that $N_U = N_\Lambda = n$. + +Now the problem is that what you really want is a linear system of +$N=n + n + m$ unknowns that includes the $m = 4$ degrees of +freedom for the source strengths we are trying to optimize for: We +need a matrix of size $N\times N$, and solution and right hand side +vectors of size $N$. Moreover, if you try to build a block linear +system because that's how your linear solver is structured, you +probably want a $3\times 3$ block structure for the matrix, and three +blocks for the vectors, even though the finite element and the +DoFHandler only know of two solution components. + +Conceptually, of course, there is nothing wrong with creating matrices +and vectors that have sizes that don't match the number of unknowns in +a DoFHandler. You just need a way to translate between the index of a +degree of freedom in a DoFHandler and the corresponding row and column +in a linear system. Here, perhaps one would make the choice that the +mesh-based degrees of freedom (which the DoFHandler numbers from zero +to $n+n-1$) come first, and the four $C_k$ additional scalar +unknowns come last. However, the fact that we need to spell out such +a convention should give you a hint that something with this kind of +design is wrong. Most codes which would use a convention like this +will, almost certainly, leave this convention implicit (something +for the programer to know about) rather than enforcing it +explicitly (say, in the form of a function `mesh_dof_to_linear_system_index()`). +It is then very easy for a number of bugs to emerge. +For example, we can not renumber all degrees of freedom with the +functions in DoFRenumbering in arbitrary ways if we want the $m$ +additional variables always at the end. + +Precisely because it is so easy to make mistakes and because it is so +restricting to have implicit conventions like "Yes, your matrix can be +larger than the number of DoFs in the DoFHandler; but it can't be +smaller", deal.II a long time ago made the decision that a lot of +functions test that the sizes of vectors and matrices are equal to the +number of unknowns in a DoFHandler. For example, +DoFTools::make_sparsity_pattern() takes a DoFHandler and expects that +the sparsity pattern has as many rows and columns as there are DoFs in +the DoFHandler. So does +DoFTools::make_hanging_node_constraints(). VectorTools::interpolate_boundary_values() +uses the same convention. Many many other functions do too, and they +all in fact check that this is true and will terminate the program if +it is not. This is because one can find many many bugs in codes this +way, and it also allows for an assumption that is so obvious that we +have never stated it so far: The degree of freedom with index $i$ +corresponds to row $i$ of the linear system and column $i$ of the +matrix. If we allowed sizes to differ, we would need to have a way to +make this connection. Further, that way would need to be be explicit (by way of a +function that translates from DoF index to row and column, as already mentioned above) or there +would be no end to the possibility of introducing bugs. + +In other words, the fact that we want to have $2n$ degrees +of freedom in the DoFHandler, but matrices and vectors with +$2n + m$ rows and columns does not work with deal.II's +underlying design decisions. Creating matrices and vectors that happen +to be larger is going to trigger lots of assertions, and for good +reasons. + + + +

How do we address this?

+ +So then how are we going to address this? At its core, what needs to +happen is that the DoFHandler knows about the additional degrees of +freedom. The question is how we are going to teach it about these +"nonlocal" degrees of freedom -- they are, after all, not ones that +are associated with vertices, edges, faces, or cells of the mesh as +the usual finite element degrees of freedom are. The approach we are +going to use is not what one would generally call "obvious" or +"natural", but it does fit into the infrastructure deal.II provides +(even though this infrastructure was originally built for other +purposes). In fact, at the time this program is written in 2024, there +is substantial debate on whether we can implement *better* ways to achieve +the same result that do not require this kind of "abuse" of features intended +for other purposes; if consensus is ever reached on these alternative +approaches, this program may also be converted. + +Whatever the future, let us turn to how our goal can be achieved today. +The key to the following is that deal.II (i) allows using different +finite elements on different cells, originally intended to support +$hp$ finite elements; and (ii) has elements that have no degrees of +freedom. We have used these features in other tutorial programs: +step-27 and step-75 make use of feature (i) to implement the $hp$ +finite element method, and step-46 makes use of feature (ii) for +multiphysics problems where some +variables live only on some cells. + +We will (ab)use the general approach used in step-46 for our purposes +here. Specifically, on the first four cells of the mesh, we will +define a finite element that has three components: a temperature component and a +Lagrange multiplier component (both of which are the usual FE_Q +element), plus a third vector component that is piecewise constant +(FE_DGQ(0)). On all other +cells, the finite element used has the same temperature and Lagrange +multiplier components, but the third component is of type FE_Nothing +and so has no degrees of freedom at all. Note that with this setup, +the first four cells will have one additional degree of freedom. +The end result of this +approach is that nearly every function you have ever used will work as +one would expect. In particular: + +- We have exactly four additional degrees of freedom, located on the + first four cells. +- We can use DoFRenumbering::component_wise() to sort degrees of freedom + explicitly into the order described above (rather than using an + implicit convention). In particular, if we really wanted to, we could + also ask for the nonlocal degrees of freedom to be numbered *first* + rather than last -- because they are on equal footing to all of the other + field-based variables. +- We can use functions such as + DoFTools::count_dofs_per_fe_block() or DoFTools::count_dofs_per_fe_component() + to create a $3\times 3$ block structure + of the linear system, because the finite elements used on all cells have + three vector components. +- The size of matrices and vectors corresponds to the number of degrees + of freedom in the DoFHandler. +- The additional degrees of freedom end up in the output files created by + DataOut because they are finite element fields -- just fields of a very + specific kind. + +What we have with this setup is a good way to not break things in deal.II. +However, it is not so natural that we can use the full extent of features +in deal.II without issue. In particular, to the DoFHandler the +nonlocal degrees of freedom only live on one cell each. On the other +hand, from the perspective of the PDE, we need these nonlocal DoFs on +all cells that intersect the areas in which the heating is applied. So, +we cannot use `cell->get_dof_indices()` to find DoF +indices for the third vector component on most cells, because the DoF handler +does not see that these dofs reach beyond the cell they are defined on. As a consequence, +we will have to do some extra work when building sparsity patterns, and later when building the +system matrix/system right hand side. + + +

What concrete steps do we need to take?

+ +There are four places in the standard finite element code where we will +need to take special care: assigning FE_DGQ elements, constructing the sparsity pattern, constructing +the system matrix/system right hand side, and outputting results. + +1. First, we need to determine on which cells the FE_Nothing element becomes +an FE_DGQ element. We do this in `setup_system()`, before dofs are distributed. +Since we are using FE_DGQ elements, we don't need to worry about boundary +conditions affecting the nonlocal dofs. We simply enumerate the cells +until we have collected enough dofs, which in this program depends on the +dimension of the problem. Each time we encounter a new cell, we +change the hp finite element from an FESystem with a 3rd FE_Nothing +component to an FESystem with a third FE_DGQ component. Calling +DoFHandler::distribute_dofs() after this is all we need to do to have the system set +up correctly. Next, we need to grab the indices of the FE_DGQ dofs. We +do this by calling DoFTools::extract_dofs(), using a component mask to +pick just the third component of the FESystem. We store these indices +in a standard vector, and use them for the three parts below. + +2. When we construct the sparsity pattern, deal.II implicitly assumes that +only dofs that share a cell could have a non-zero entry in the matrix. +Consequently, when working with non-local dofs, the sparsity pattern that +DoFTools::make_sparsity_pattern() creates will omit entries, on all +but the first four cells, that in fact +we want to fill. To fix this, after creating the sparsity pattern we +need to loop over all cells in the mesh and check whether the dofs on +that cell overlap with regions where each $f_k$ is nonzero. If they do, +then we need to add an entry to the sparsity pattern at location +(`dof_index`, `nonlocal_dof_index`) and (`nonlocal_dof_index`, `dof_index`). +Note that we need to add two entries because it doesn't matter whether +the nonlocal dof is the row or column index, the entry could be +nonzero. In this example, we further refine this by only adding entries +where the $\Lambda$ component of the FESystem is nonzero, because we only +deal with the nonlocal functions in the equations stemming from the +$\Lambda^j$ and $C^k$ derivatives. + +3. We take a similar approach when constructing the system matrix. We do +the standard loop over cells, quadrature points, and dof indices. Normally, +we would simply integrate the finite elements on the cell and add these +values to a local matrix. However, because we are dealing with nonlocal dofs +this approach won't work. Instead, we use the standard approach for +constructing the $\mathcal{M}$ and $\mathcal{N}$ portions of the block system +matrix, but for the $\mathcal{F}$ portions we need to numerically integrate +each finite element against the nonlocal $f_k$ functions, and enter the +resulting value in the correct location in the system matrix. To do the integration, +we simply instantiate a Point object for the quadrature point and pass this +to the $f_k$ functions, keeping all other aspects of integration the same. We +then use the extracted dof indices to add the integrated value to the system +matrix, using the loop's current dof index as the other index. The system right hand +side is constructed as usual, integrating the right hand side function on each cell. + +4. Finally, to properly visualize the solution and constituent parts of the program, +we need to correctly interpret the nonlocal dofs. We do this by interpolating the +$f_k$ functions onto a smaller DoFHandler (one which has only one FE_Q element), +then multiplying by the corresponding coefficient $c_k$. For this program, we also +add together all the individual $C^kf_k$ to visualize $f$, and interpolate the +target function as well. These steps (for the $f_k$) all occur in a loop, which makes +scaling the number of nonlocal dofs relatively straightforward. + +To see how we actually do these things, let's look at the code. diff --git a/examples/step-93/doc/kind b/examples/step-93/doc/kind new file mode 100644 index 0000000000..c1d9154931 --- /dev/null +++ b/examples/step-93/doc/kind @@ -0,0 +1 @@ +techniques diff --git a/examples/step-93/doc/results.dox b/examples/step-93/doc/results.dox new file mode 100644 index 0000000000..0f3af8f0ee --- /dev/null +++ b/examples/step-93/doc/results.dox @@ -0,0 +1,188 @@ +

Results

+ +When you run the program with the step target function (in 2D), the output looks something like this: + +@code +Number of active cells: 16384 +Number of degrees of freedom: 66049+66049+4 = 132102 +Number of nonlocal dofs: 4 +Beginning solve... +Wall time: 63.9265s +Solved in 39973 MINRES iterations. +The c coefficients are + c1: 28.7408 + c2: -6.51604 + c3: -6.51604 + c4: -1.62044 +@endcode + +When you run the program with the Gaussian target function (in 2D), the output should look like this: + +@code +Number of active cells: 16384 +Number of degrees of freedom: 66049+66049+4 = 132102 +Number of nonlocal dofs: 4 +Beginning solve... +Wall time: 98.4858s +Solved in 62131 MINRES iterations. +The c coefficients are + c1: 23.553 + c2: -4.86562 + c3: -4.86562 + c4: -1.42344 +@endcode + +The goal of this program is to determine which temperature settings best match the target function, so first +let's see what these targets look like: + + + + + + + + + + +
+
Step Target %Function
+
+
Gaussian Target %Function
+
+ Step target function + + Gaussian target function +
+ +After solving the Lagrangian system, we arrive at solutions $U_\text{step}$ and $U_\text{gauss}$ that +look like this: + + + + + + + + + + +
+
$U_\text{step}$
+
+
$U_\text{gauss}$
+
+ Solution for step shaped target function + + Solution for Gaussian target function +
+ +Notice that $U_\text{gauss}$ matches the target much better than +$U_\text{step}$. Intuitively, this makes sense: in general, solutions +to the heat equation look something like Gaussians, so the +Gaussian target function is a much more "natural" thing to match than +a sharp step function. We can also see this in the optimal heat +profiles. + + + + + + + + + + +
+
Heat plate settings for matching step function
+
+
Heat plate settings for matching Gaussian
+
+ Heat plate settings for matching step function + + Heat plate settings for matching Gaussian +
+ +Notice that for the Gaussian target, the 4 plates are set to less extreme values. In contrast, to try to match the step function, higher and lower temperatures must be applied. + +While it does not contain much useful information, we can also plot the Lagrange multiplier $\Lambda$, which has an interesting shape: + + + + + + + + + + +
+
$\Lambda_\text{step}$
+
+
$\Lambda_\text{gauss}$
+
+ Lagrange multiplier for step target function + + Lagrange multiplier for Gaussian target function +
+ + + +

Possibilities for extensions

+ +There are a few ways that this program could be extended, which we list below. + +1. As mentioned in the code documentation, this program does not make +use of any preconditioners before solving. This is because, for a 2D +problem, the code runs fast enough that optimization is not +necessary. However, as shown in the screen output above, the number of +iterations required to solve the linear system is quite large. Thus, +for larger problems, it would be good if the solver ran more +quickly. See the "Possibilities for extensions" section of step-6 for +a more detailed discussion on how to change preconditioners. We should +note that since the block matrix we use has many zeros on the +diagonal, preconditioners like PreconditionJacobi will not work +because they divide by diagonal entries. Instead, block +preconditioners such as those discussed in step-20 or step-22 (among +many others) will likely be useful. + +2. To validate the optimization problem is working correctly, we could +try to match a target function which is itself a solution to the +Poisson equation with prescribed heat profile. If the optimization +problem is being solved correctly, it should be able to perfectly +match this solution. To create such a function, we would need to first +solve the Poisson problem on a scalar field, with a RHS described by +the chosen heat profile. See step-7 for more information on the method +of manufactured solutions. + +3. The program at the moment has the number of nonlocal degrees of freedom +hardcoded as $2^d$ (see the constructor). We then assign each of these +degrees of freedom to one of the first $2^d$ cells. This is not going to +be much of a problem because there are always enough cells for this as +long as you start with a mesh that is at least once refined. What would +we do if we had a number of nonlocal DoFs that is not easily predictable, +and that may be larger than the number of cells? Perhaps a better approach +would be to come up with a way to assign *all* of these to the first cell, +because there is *always* at least one cell. The way to achieve this would +be to replace the use of FE_DGQ(0) (an element with exactly one degree of +freedom) by an element that has more than one -- and in particular exactly +the right number of degrees of freedom -- and that can be used on the first +cell; all other cells would then use FE_Nothing. At the time of writing +this program, there is no element class that can easily be given a *specific* +number of degrees of freedom, but it would not be very difficult to write +such a class. diff --git a/examples/step-93/doc/tooltip b/examples/step-93/doc/tooltip new file mode 100644 index 0000000000..9a059f3a4d --- /dev/null +++ b/examples/step-93/doc/tooltip @@ -0,0 +1 @@ +Working with nonlocal DoFs. Solving an optimal control problem with constrained input. diff --git a/examples/step-93/step-93.cc b/examples/step-93/step-93.cc new file mode 100644 index 0000000000..5057f4fd55 --- /dev/null +++ b/examples/step-93/step-93.cc @@ -0,0 +1,821 @@ +/* ------------------------------------------------------------------------ + * + * SPDX-License-Identifier: LGPL-2.1-or-later + * Copyright (C) 2024 by the deal.II authors + * + * This file is part of the deal.II library. + * + * Part of the source code is dual licensed under Apache-2.0 WITH + * LLVM-exception OR LGPL-2.1-or-later. Detailed license information + * governing the source code and code contributions can be found in + * LICENSE.md and CONTRIBUTING.md at the top level directory of deal.II. + * + * ------------------------------------------------------------------------ + * + * Authors : Sam Scheuerman, Wolfgang Bangerth, Colorado State University, 2024. + */ + + +// @sect3{Include files and other top matter} + +#include +#include +#include + +#include +#include + +#include +#include +#include +#include + +#include +#include + +#include +#include + +#include +#include +#include +#include +#include +#include +#include + +#include +#include + +#include +#include + + +namespace Step93 +{ + using namespace dealii; + + + // @sect3{Target and indicator functions} + + // We start by defining a function class for $\bar u$; it represents the heat + // profile we want to match. + // + // This class has two member variables: + // + // - `center`: A Point object representing the center of the function + // - `radius`: A double representing the radius of the circular step, or the + // standard deviation of a Gaussian + template + class TargetFunction : public Function + { + public: + TargetFunction() + : Function(1) + , center() + , radius(0.1) + {} + + // Next, we define an overloaded constructor for TargetFunction, so we can + // choose to set the center and radius + // + // Parameters: + // + // - `center`: A constant Point pointer used to set the member variable + // center + // - `radius`: A constant double used to set the member variable radius + TargetFunction(const unsigned int n_components, + const Point ¢er, + const double radius = .3) + : Function(n_components) + , center(center) + , radius(radius) + {} + + virtual double value(const Point &p, + const unsigned int component = 0) const override; + + private: + const Point center; + const double radius; + }; + + // The value() function returns the value of the function at point `p` and + // component index `component`. In this case, if the component corresponds to + // the solution $u$, then the function returns a value based on either a step + // function or a Gaussian. If the component corresponds to any of the other + // variables, the function always returns 0. + // + // Note: In the documentation we discuss using a Gaussian target function, but + // this is not included in the code for readability. To run the code with a + // Gaussian target function, simply replace the code in the `component == 0` + // braces by + // + // `return std::exp(-((p-center).norm()*(p-center).norm())/(radius*radius));` + template + double TargetFunction::value(const Point &p, + const unsigned int component) const + { + if (component == 0) + { + if ((p - center).norm() <= radius) + return 1; + else + return 0; + } + else + return 0; + } + + // The next class we define is one for a circular indicator function. Unlike + // the target function, this does not need a component argument because we + // have to manually address where it gets used in the code. Objects of this + // function type correspond to the nonlocal dofs. + // + // Parameters: + // + // - `center`: A constant Point object giving the center of the indicator + // region. + // - `radius`: The radius of the region. + template + class CircularIndicatorFunction : public Function + { + public: + CircularIndicatorFunction(); + + CircularIndicatorFunction(const Point ¢er, const double radius); + + virtual double value(const Point &p, + const unsigned int component = 0) const override; + + private: + const Point center; + const double radius; + }; + + template + CircularIndicatorFunction::CircularIndicatorFunction() + : center(Point()) + , radius(20) + {} + + + template + CircularIndicatorFunction::CircularIndicatorFunction( + const Point ¢er, + const double radius) + : center(center) + , radius(radius) + {} + + template + double CircularIndicatorFunction::value( + const Point &p, + const unsigned int /* component */) const + { + if ((center - p).norm() <= radius) + return 1; + else + return 0; + } + + + // @sect3{The principal class} + + // The main class is very similar to step-4 in structure, given that this is + // a relatively simple program that does not use adaptive mesh refinement. + // However, there are three new member variables: + // + // - `nonlocal_dofs`: A `std::vector` of dof indices that stores the + // dof index for the nonlocal dofs. + // + // - `heat_functions`: A `std::vector` of CircularIndicatorFunction objects; + // these are the heat sources. + // + // - `target_function`: This is the function we want to match. We store it + // as a class variable because it is used both in assemble_system() and + // output_results(). + template + class Step93 + { + public: + Step93(); + + void run(); + + private: + void make_grid(); + void setup_system(); + void assemble_system(); + void solve(); + void output_results() const; + + Triangulation triangulation; + DoFHandler dof_handler; + + hp::FECollection fe_collection; + hp::QCollection quadrature_collection; + + AffineConstraints constraints; + + SparsityPattern sparsity_pattern; + SparseMatrix system_matrix; + + Vector solution; + Vector system_rhs; + + std::vector nonlocal_dofs; + + std::vector> heat_functions; + + const TargetFunction target_function; + }; + + + // @sect4{The Step93 constructor} + // In this constructor, we set up several of the fundamental data + // structures this program needs. Specifically, we generate the + // finite element collection, which is basically a list of all the + // possible finite elements we could use on each cell. This + // collection has two elements: one FESystem that has two degree 2 + // FE_Q elements and one FE_Nothing element (to be used on all cells + // that are not "special"), and one FESystem that has two degree 2 + // FE_Q elements and one degree 0 FE_DGQ element (to be used on + // those "special" cells we use to anchor the non-local degrees of + // freedom, as discussed in the introduction). + // + // Where we have a collection of elements, we then also need a + // collection of quadratures -- which here has only a single element + // because we can use the same quadrature on all cells. + + // The default constructor also initializes + // dof_handler and target_function, + // and it generates a vector of CircularIndicatorFunctions objects which + // represent the heat functions. + template + Step93::Step93() + : dof_handler(triangulation) + , target_function(3, + (dim == 1 ? Point(0.5) : + dim == 2 ? Point(0.5, 0.5) : + Point(0.5, 0.5, 0.5))) + { + // Here, we generate the finite element collection, which is basically a + // list of all the possible finite elements we could use on each cell. This + // collection has two elements: one FESystem that has two degree 2 FE_Q + // elements and one FE_Nothing element, and one FESystem that has two + // degree 2 FE_Q elements and one degree 0 FE_DGQ element. + fe_collection.push_back( + FESystem(FE_Q(2), 2, FE_Nothing(), 1)); + fe_collection.push_back(FESystem(FE_Q(2), 2, FE_DGQ(0), 1)); + + // The quadrature collection is just one degree 3 QGauss element. + quadrature_collection.push_back(QGauss(3)); + + // Here, we choose the center points for the heat functions to be the + // vertices of a lattice, in this case the corners of a hypercube + // centered at 0, with side length 1. The block of code below + // then creates the CircularIndicatorFunction + // objects for the main class. To make these, we check the dimension of the + // problem and create 2, 4, or 8 function objects. + switch (dim) + { + case (1): + heat_functions.emplace_back(Point({0.5}), 0.2); + heat_functions.emplace_back(Point({-0.5}), 0.2); + break; + case (2): + heat_functions.emplace_back(Point({0.5, 0.5}), 0.2); + heat_functions.emplace_back(Point({0.5, -0.5}), 0.2); + heat_functions.emplace_back(Point({-0.5, 0.5}), 0.2); + heat_functions.emplace_back(Point({-0.5, -0.5}), 0.2); + break; + case (3): + heat_functions.emplace_back(Point({0.5, 0.5, 0.5}), 0.2); + heat_functions.emplace_back(Point({0.5, 0.5, -0.5}), 0.2); + heat_functions.emplace_back(Point({0.5, -0.5, 0.5}), 0.2); + heat_functions.emplace_back(Point({0.5, -0.5, -0.5}), 0.2); + heat_functions.emplace_back(Point({-0.5, 0.5, 0.5}), 0.2); + heat_functions.emplace_back(Point({-0.5, 0.5, -0.5}), 0.2); + heat_functions.emplace_back(Point({-0.5, -0.5, 0.5}), 0.2); + heat_functions.emplace_back(Point({-0.5, -0.5, -0.5}), 0.2); + break; + default: + DEAL_II_ASSERT_UNREACHABLE(); + break; + } + } + + // @sect4{Step93::make_grid()} + + // The make_grid() function makes a hypercube grid, see step-4: + template + void Step93::make_grid() + { + GridGenerator::hyper_cube(triangulation, -1, 1); + triangulation.refine_global(7); + + std::cout << "Number of active cells: " << triangulation.n_active_cells() + << std::endl; + } + + + // @sect4{Step93::setup_system()} + + // The `setup_system()` function is similar to step-4, except we have to add a + // few steps to prepare for the nonlocal dofs. + template + void Step93::setup_system() + { + // Here, we loop over the cells and set the FESystem index to 1, which + // corresponds to the system with 2 FE_Q elements and one FE_DGQ element. We + // do this until we have enough dofs for each heat function. Note that we + // use the global active cell index to measure when to stop the + // loop. This allows the loop to be run in parallel with no alteration. + // Then, we call DoFHandler::distribute_dofs() to actually enumerate all + // degrees of freedom. + for (const auto &cell : dof_handler.active_cell_iterators()) + { + if (cell->global_active_cell_index() < heat_functions.size()) + { + cell->set_active_fe_index(1); + } + else + { + break; + } + } + dof_handler.distribute_dofs(fe_collection); + + // Once we've assigned dofs, the code block below counts the + // number of dofs in the system, and outputs to the console. In + // other contexts, we might want to use *block* matrices (see, for + // example, step-20 or step-22) to build more efficient linear + // solvers; here, we will just put everything into one big matrix + // and so knowing the number of unknowns for each of the variables + // is purely for informational purposes. + const std::vector dofs_per_component = + DoFTools::count_dofs_per_fe_component(dof_handler); + const unsigned int n_dofs_u = dofs_per_component[0], + n_dofs_l = dofs_per_component[1], + n_dofs_c = dofs_per_component[2]; + std::cout << "Number of degrees of freedom: " << n_dofs_u << "+" << n_dofs_l + << "+" << n_dofs_c << " = " << n_dofs_u + n_dofs_l + n_dofs_c + << std::endl; + + // Finally, we need to extract the indices of the finite elements which + // correspond to the non-local dofs. + // + // First, we make a component mask, which is `false` except for + // the third component. This will extract only the dofs from the + // third component of the FE system. Next, we actually extract the + // dofs, and store them in an IndexSet variable. Finally, we add + // each extracted index to the member array `nonlocal_dofs`. + const ComponentMask component_mask_c({false, false, true}); + const IndexSet indices_c = + DoFTools::extract_dofs(dof_handler, component_mask_c); + + for (const types::global_dof_index non_local_index : indices_c) + nonlocal_dofs.push_back(non_local_index); + + std::cout << "Number of nonlocal dofs: " << nonlocal_dofs.size() + << std::endl; + + + // The mesh we created above is not locally refined, and so there + // are no hanging node constraints to keep track of. But it does + // not hurt to just use the same setup we have used starting in + // step-6 of building a constraints object that contains hanging + // node constraints (anticipating that perhaps we'd want to do + // adaptive mesh refinement in a later step) into which we then + // also put the constraints for boundary values on both the $u$ + // and $\lambda$ variables. + // + // Because the nonlocal degrees of freedom use discontinuous + // elements, they do not contribute to boundary values + // (discontinuous elements do not have degrees of freedom + // logically located on the boundary that could be interpolated) + // and we do not need to exclude these solution components + // explicitly when calling + // VectorTools::interpolate_boundary_values() + constraints.clear(); + DoFTools::make_hanging_node_constraints(dof_handler, constraints); + VectorTools::interpolate_boundary_values(dof_handler, + 0, + Functions::ZeroFunction(3), + constraints); + constraints.close(); + + // The remainder of the function deals with building the sparsity + // pattern. It consists of two parts: The entries that result from + // the usual integration of products of shape functions, and then + // the entries that result from integrals that contain nonlocal + // degrees of freedom (which one can think of as associated with + // shape functions that are constant across the entire + // domain). The first part is easily built using standard tools: + DynamicSparsityPattern dsp(dof_handler.n_dofs(), dof_handler.n_dofs()); + DoFTools::make_sparsity_pattern(dof_handler, dsp, constraints, false); + + // The other part is more awkward. We have matrix entries that + // result from terms such as $\int_{\Omega}\varphi_j f_k$ where + // each $\varphi_j$ is a shape function associated to $\lambda$ + // and each $f_k$ is a + // characteristic function of a part of the domain. As a + // consequence, we end up with a nonzero matrix entry $A_{jk}$ + // (along with its transpose $A_{kj}$) if there is overlap between + // $\varphi_j$ and $f_k$. In practice, because we will use + // quadrature, this means that we end up with a quadrature point + // on a cell on which $\varphi_j$ lives and at which $f_k$ is not + // zero. (We will implicitly assume that a shape function that + // lives on the current cell is nonzero at all quadrature points + // -- an assumption that is generally true unless one chooses + // specific quadrature formulas.) Determining which sparsity + // pattern entries we need to add then essentially comes down to + // "simulating" what would happen if we actually computed + // integrals and which matrix entries would end up being non-zero + // in the process. As a consequence, the following code's general + // structure looks very similar to what we will do for the + // nonlocal contributions in the `assemble_system()` function + // below. + // + // To get this started, we create an `hp_fe_values` that we + // will only use to query the quadrature point locations. The + // non-local dofs will need to interact with the second component + // of the fe system (namely, $\lambda$), so we also declare a + // variable that will help us extract this scalar field for use + // below. + hp::FEValues hp_fe_values(fe_collection, + quadrature_collection, + update_quadrature_points); + + const FEValuesExtractors::Scalar lambda(1); + + // Then, we loop over the cells, then over the quadrature points, and + // finally over the indices, as if we were constructing a mass matrix. + // However, what we instead do here is check two things. First, we check if + // the quadrature point is within the radius of a circular indicator + // function that represents our non-local dof. If so + // then we add an entry to the sparse matrix at the + // (nonlocal dof index, lambda dof index) entry and the (lambda dof index, + // nonlocal dof index) entry for all lambda degrees of freedom. (Because the + // matrix we solve with has both the lambda-nonlocal interacting block and + // its transpose, we need to add two entries each time.) + std::vector local_dof_indices; + for (const auto &cell : dof_handler.active_cell_iterators()) + { + hp_fe_values.reinit(cell); + + const FEValues &fe_values = hp_fe_values.get_present_fe_values(); + + local_dof_indices.resize(fe_values.dofs_per_cell); + cell->get_dof_indices(local_dof_indices); + + for (const unsigned int q_index : fe_values.quadrature_point_indices()) + { + const Point q_point = fe_values.quadrature_point(q_index); + for (const unsigned int i : fe_values.dof_indices()) + { + if (fe_values.get_fe().system_to_component_index(i).first == + 1) // 'i' is a lambda shape function + { + for (unsigned int j = 0; j < heat_functions.size(); ++j) + if (heat_functions[j].value(q_point) != 0) + { + dsp.add(local_dof_indices[i], nonlocal_dofs[j]); + dsp.add(nonlocal_dofs[j], local_dof_indices[i]); + } + } + } + } + } + + // The rest (below) is standard setup code, see step-4: + sparsity_pattern.copy_from(dsp); + + system_matrix.reinit(sparsity_pattern); + + solution.reinit(dof_handler.n_dofs()); + system_rhs.reinit(dof_handler.n_dofs()); + } + + + // @sect4{Step93::assemble_system()} + + // The `assemble_system()` function works very similar to how is + // does in other tutorial programs (cf. step-4, step-6, step-8, and + // for the vector-valued case see step-22). However, there is an + // additional component to constructing the system matrix, because + // we need to handle the nonlocal dofs manually. + template + void Step93::assemble_system() + { + // First, we do a standard loop setup for constructing the system matrix. + hp::FEValues hp_fe_values(fe_collection, + quadrature_collection, + update_values | update_gradients | + update_quadrature_points | + update_JxW_values); + + FullMatrix cell_matrix; + Vector cell_rhs; + + std::vector local_dof_indices; + + const FEValuesExtractors::Scalar u(0); + const FEValuesExtractors::Scalar lambda(1); + const FEValuesExtractors::Scalar c(2); + + for (const auto &cell : dof_handler.active_cell_iterators()) + { + const unsigned int dofs_per_cell = cell->get_fe().n_dofs_per_cell(); + + cell_matrix.reinit(dofs_per_cell, dofs_per_cell); + cell_rhs.reinit(dofs_per_cell); + hp_fe_values.reinit(cell); + + cell_matrix = 0; + cell_rhs = 0; + + const FEValues &fe_values = hp_fe_values.get_present_fe_values(); + + local_dof_indices.resize(fe_values.dofs_per_cell); + + cell->get_dof_indices(local_dof_indices); + + + // In the loop over quadrature points, we start by building + // all of the usual terms that are bilinear in shape functions + // corresponding to the $u$ and $\lambda$ variables: + for (const unsigned int q_index : fe_values.quadrature_point_indices()) + { + const double JxW = fe_values.JxW(q_index); + for (const unsigned int i : fe_values.dof_indices()) + { + const double phi_i_u = fe_values[u].value(i, q_index), + phi_i_l = fe_values[lambda].value(i, q_index); + + const Tensor<1, dim> grad_i_u = + fe_values[u].gradient(i, q_index), + grad_i_l = + fe_values[lambda].gradient(i, q_index); + + for (const unsigned int j : fe_values.dof_indices()) + { + const double phi_j_u = fe_values[u].value(j, q_index); + + const Tensor<1, dim> grad_j_u = + fe_values[u].gradient(j, q_index), + grad_j_l = + fe_values[lambda].gradient(j, + q_index); + + cell_matrix(i, j) += phi_i_u * phi_j_u * JxW; + cell_matrix(i, j) += -grad_i_u * grad_j_l * JxW; + cell_matrix(i, j) += -grad_i_l * grad_j_u * JxW; + } + + const Point q_point = fe_values.quadrature_point(q_index); + cell_rhs(i) += (phi_i_u * target_function.value(q_point) * JxW); + + + // For the integrals that involve the nonlocal dofs, + // we make use of the quadrature point again. To + // compute the integrals, we loop over the + // heat functions, adding the numeric integral of each + // heat equation with each $\lambda$ shape function, + // at the appropriate indices (which we found in + // `setup_system()`). Note that if we try to add 0 to + // a matrix entry we have not previously indicated should + // be nonzero, there will not be a problem; but if we + // try to add a nonzero value to an entry not + // previously added to the sparsity pattern, we will + // get an error. In other words, the following lines + // of the code check that we adjusted the sparsity + // pattern correctly in the previous function. + for (unsigned int j = 0; j < heat_functions.size(); ++j) + { + system_matrix.add(local_dof_indices[i], + nonlocal_dofs[j], + heat_functions[j].value(q_point) * + phi_i_l * JxW); + system_matrix.add(nonlocal_dofs[j], + local_dof_indices[i], + heat_functions[j].value(q_point) * + phi_i_l * JxW); + } + } + } + + // Finally, we copy the local contributions to the linear + // system into the global matrix and right hand side vector, + // taking into account hanging node and boundary values + // constraints: + constraints.distribute_local_to_global( + cell_matrix, cell_rhs, local_dof_indices, system_matrix, system_rhs); + } + } + + + // @sect4{Step93::solve()} + + // The solve() function works similar to how it is done in step-6 + // and step-8, except we need to use a different solver because the + // linear problem we are trying to solve is a saddle point problem + // for which the Conjugate Gradient algorithm is not + // applicable. But, because the matrix is symmetric, we can use + // SolverMinRes, an iterative solver specialized for symmetric + // indefinite problems. This solver could be improved with the use + // of preconditioners, but we don't do that here for simplicity (see + // the Possibilities for Extensions section below). + // + // As you will see in the output, given that we are not using a + // preconditioner, we need a *lot* of iterations to solve this + // linear system. We set the maximum to one million, more than we + // need of course, but an indication that this is not an efficient + // solver. For smaller problems, one can also use a direct solver + // (see step-29) for which you would just replace the main part of + // this function by the following three lines of code: + // + // `SparseDirectUMFPACK direct_solver;` + // + // `direct_solver.initialize(system_matrix);` + // + // `direct_solver.vmult(solution, system_rhs);` + template + void Step93::solve() + { + std::cout << "Beginning solve..." << std::endl; + Timer timer; + + SolverControl solver_control(1'000'000, 1e-6 * system_rhs.l2_norm()); + SolverMinRes> solver(solver_control); + + solver.solve(system_matrix, solution, system_rhs, PreconditionIdentity()); + + timer.stop(); + std::cout << "Wall time: " << timer.wall_time() << "s" << std::endl; + std::cout << "Solved in " << solver_control.last_step() + << " MINRES iterations." << std::endl; + } + + + // @sect4{Step93::output_results()} + + // The `output_results()` function is a bit more robust for this program than + // is typical. This is because, in order to visualize the heat sources we have + // optimized, we need to do extra work and interpolate them onto a mesh. We do + // this by instantiating a new DoFHandler object and then using the helper + // function VectorTools::interpolate(). + // + // The top of the function is as always when using vector-valued + // elements (see, for example, step-22) and simply outputs all of + // the solution variables on the mesh cells they are defined on: + template + void Step93::output_results() const + { + const std::vector solution_names = {"u", "lambda", "c"}; + + const std::vector + interpretation = {DataComponentInterpretation::component_is_scalar, + DataComponentInterpretation::component_is_scalar, + DataComponentInterpretation::component_is_scalar}; + + DataOut data_out; + data_out.add_data_vector(dof_handler, + solution, + solution_names, + interpretation); + + // The non-local degrees of freedom are of course defined on the + // first several cells, but at least logically are considered to + // live everywhere (or, if you prefer, nowhere at all since they + // do not represent spatial functions). But, conceptually, we use + // them as multipliers for the heat sources, and so while the + // coefficients $C^k$ are non-local, the heat source $\sum_k C^k + // f_k(\mathbf x)$ *is* a spatially variable function. It would be + // nice if we could visualize that as well. The same is true for + // the target heat distribution $\bar u$ we are trying to match. + // + // To do so, we + // create a new dof handler to output the target function + // and heat plate values, and associate it with + // a finite element with a degree that matches what we used to solve for $u$ + // and $\lambda$, although in reality this is an arbitrary choice: + DoFHandler new_dof_handler(triangulation); + + const FE_Q new_fe(2); + new_dof_handler.distribute_dofs(new_fe); + + // To get started with the visualization, we need a vector which + // stores the interpolated target function. We create the vector, + // interpolate the target function $\bar u$ onto the mesh, then + // add the data to our data_out object. + Vector target(new_dof_handler.n_dofs()); + VectorTools::interpolate(new_dof_handler, target_function, target); + data_out.add_data_vector(new_dof_handler, target, "u_bar"); + + // In order to visualize the sum of the heat sources $\sum_k C^k + // f_k(\mathbf x)$, we create a vector which will store the + // interpolated values of this function. Then, we loop through + // the heat functions, create a vector to store the interpolated + // data, call the VectorTools::interpolate() function to fill the + // vector, multiply the interpolated data by the nonlocal dof + // value $C^k$ (so that the heat plate is set to the correct + // temperature), and then add this data to the sum of heat + // sources. Because we can, we also add the vector for each source + // individually to the DataOut object, so that they can be + // visualized individually. + Vector full_heat_profile(new_dof_handler.n_dofs()); + + for (unsigned int i = 0; i < heat_functions.size(); ++i) + { + Vector hot_plate_i(new_dof_handler.n_dofs()); + + VectorTools::interpolate(new_dof_handler, + heat_functions[i], + hot_plate_i); + + hot_plate_i *= solution[nonlocal_dofs[i]]; + full_heat_profile += hot_plate_i; + + const std::string data_name = + "Heat_Source_" + Utilities::int_to_string(i); + data_out.add_data_vector(new_dof_handler, hot_plate_i, data_name); + } + + // Once all the heat functions have been combined, we add them to the + // data_out object, and output everything into a file: + data_out.add_data_vector(new_dof_handler, + full_heat_profile, + "Full_Heat_Profile"); + + data_out.build_patches(); + + std::ofstream output("solution.vtu"); + data_out.write_vtu(output); + + // Finally, we output the nonlocal coefficient values to the console: + std::cout << "The c coefficients are " << std::endl; + for (long unsigned int i = 0; i < nonlocal_dofs.size(); ++i) + { + std::cout << "\tc" << i + 1 << ": " << solution[nonlocal_dofs[i]] + << std::endl; + } + } + + + // @sect4{Step93::run()} + + // The run() function runs through each step of the program, nothing new here: + template + void Step93::run() + { + make_grid(); + setup_system(); + assemble_system(); + solve(); + output_results(); + } +} // namespace Step93 + + +// @sect3{The main() function} + +// The `main()` function looks essentially like that of most other tutorial +// programs. +int main() +{ + try + { + Step93::Step93<2> heat_optimization_problem; + heat_optimization_problem.run(); + } + catch (std::exception &exc) + { + std::cerr << std::endl + << std::endl + << "----------------------------------------------------" + << std::endl; + std::cerr << "Exception on processing: " << std::endl + << exc.what() << std::endl + << "Aborting!" << std::endl + << "----------------------------------------------------" + << std::endl; + + return 1; + } + catch (...) + { + std::cerr << std::endl + << std::endl + << "----------------------------------------------------" + << std::endl; + std::cerr << "Unknown exception!" << std::endl + << "Aborting!" << std::endl + << "----------------------------------------------------" + << std::endl; + return 1; + } + + return 0; +} diff --git a/tests/examples/step-93.diff b/tests/examples/step-93.diff new file mode 100644 index 0000000000..a2d4bc2106 --- /dev/null +++ b/tests/examples/step-93.diff @@ -0,0 +1,33 @@ +diff --git a/examples/step-93/step-93.cc b/examples/step-93/step-93.cc +index 2e7944c2b7..a0aac6263b 100644 +--- a/examples/step-93/step-93.cc ++++ b/examples/step-93/step-93.cc +@@ -374,7 +374,7 @@ namespace Step93 + void Step93::make_grid() + { + GridGenerator::hyper_cube(triangulation, -1, 1); +- triangulation.refine_global(7); ++ triangulation.refine_global(3); + + std::cout << "Number of active cells: " << triangulation.n_active_cells() + << std::endl; +@@ -655,7 +655,7 @@ namespace Step93 + { + // Notice that we also time how long this process takes. + std::cout << "Beginning solve" << std::endl; +- Timer timer; ++ //Timer timer; + + SolverControl solver_control(5000000, 1e-6 * system_rhs.l2_norm()); + SolverCG> solver(solver_control); +@@ -665,8 +665,8 @@ namespace Step93 + /* SparseDirectUMFPACK direct_solver; + direct_solver.initialize(system_matrix); + direct_solver.vmult(solution, system_rhs); */ +- timer.stop(); +- std::cout << "Wall time: " << timer.wall_time() << "s" << std::endl; ++ //timer.stop(); ++ //std::cout << "Wall time: " << timer.wall_time() << "s" << std::endl; + } + + // The output_results() function is a bit more robust for this program than