- Generated by
1.9.8
|
Saddle-Point
|
In this project, we show how to solve a system of partial differential equations (PDE) with two components on a domain \(\Omega\) in 2D. The final version of the program implementation is on the branch master. One can use the following command
to switch to the branch that contains the final version of the program. The manufactured solution that we present later is solved on the circle with radius \(\pi\). We focus on the system of PDE with the following strong and weak forms, and this problem comes from a saddle point formulation when computing the harmonic Riesz representation of a given functional in the dual space of \(H^1\).
Find \((p,u) \in H^1(\Omega) \times H^1(\Omega)\) such that
\begin{align*} u - \Delta_{\Gamma} u &= g, \\ \Delta u &= 0, \\ -\Delta p &= f_1, \end{align*}
where \(\Omega\) is the domain, \(\Gamma=\partial\Omega\) is its boundary, and \(f_1, g \in H^{-1}(\Omega)\).
The full problem we are trying to solve has the corresponding weak form
\begin{align*} \int_{\Gamma} u\,v + \int_{\Gamma} \nabla_{\Gamma} u \cdot \nabla_{\Gamma} v + \int_{\Omega} A \,\nabla p \cdot \nabla v &= \int_{\Omega} f_1\, v + \int_{\Gamma} g\, v \qquad \forall v, \\ \int_{\Omega} A \,\nabla u \cdot \nabla q &= 0 \qquad \forall q, \end{align*}
for all \(v, q \in C^{\infty}(\bar{\Omega})\), the smooth functions on the closure of the domain \(\Omega\). We have added the coefficient matrix \(A\) to the weak form. We have defined coefficient in the header file step-6.h, and we simply let it be \(1\), unless the user decides to change it to some other value.
This project starts with code from tutorial step-6, and its implementation has borrowed from step-20 for the vector-valued FESystem and from step-38 for the tangential gradient and the surface-Laplacian term.
param.prm file as input at run time. We used Tutorial step-28 as a reference.param.prm file. The FunctionParser class implements a function object that gets its value by parsing a string describing this function. It is a wrapper class for the muparser library. We referenced tutorial step-36 for using FunctionParser.For convenience, we note that the source directory has the following files:
And the include directory has the following files:
Our manufactured solution has two components. It is implemented in the solution.cc file in the source directory. The output shown in the Results section corresponds to this manufactured solution.
The first component of the manufactured solution is \(p(x,y) = e^x \cos y\), and the second component of the manufactured solution is \(u(x,y) = -2xy\). The corresponding right-hand sides are specified in the param.prm file.
Going back to the weak formulation, the first component \(p\) is not involved in the integral over the boundary, so computing the corresponding right-hand side is not too involved — we just need to use \(f = -\Delta p\) from the strong form. On the other hand, the right-hand-side function \(g = u - \Delta_{\Gamma} u\) requires a bit more work, and its computation can be found in tutorial step-38. We use the following formula for the surface Laplacian:
\begin{align*} - \Delta_{\Gamma} u &= \Delta u - n^T D^2 u\, n - (n \cdot \nabla u)\, \kappa,\\ - \Delta_{\Gamma} u &= - 8 x y, \end{align*}
where \(\kappa\) is the total curvature of the boundary \(\Gamma\). Since we take our domain to be the circle with radius \(\pi\) in this manufactured solution (and of radius \(1\) in the additional manufactured solution), the outward normal \(n\) equals the vector \(\langle x,y\rangle\), and \(\kappa = 1\). Hence, we obtain the desired right-hand-side function \(g\).
There is an additional manufactured solution that we have implemented in the program. This manufactured solution is commented out, and to run it one has to uncomment some corresponding lines. Since the output is similar to the test case we presented, we do not discuss this manufactured solution in the results section. It has the following first component:
\[ p(x,y) = \sin(\pi x)\,\sin(\pi y), \]
and the following second component:
\[ u(x,y) = -2xy. \]
To run the program with this manufactured solution instead, we uncomment the corresponding equations in the solution.cc file and also change the right-hand-side function in the param.prm file to the correct right-hand side for this exact solution. The correct right-hand side is also implemented in the rhs_function.h file. To use it, we specify TYPE_SIN_SIN when we create the object RHSFunction<dim> and use the following lines
instead of the lines
and we would be able to run the program with this manufactured solution instead.
The constructor of this class specifies the parameters, the finite elements, and associates the DoFHandler with the triangulation. (Note that the triangulation isn't set up with a mesh at this point.)
This function sets up all the variables that describe the linear system, such as the DoFHandler, matrices, and vectors. It also generates the sparsity pattern, distributes the degrees of freedom, builds hanging-node and Dirichlet boundary constraints (using the manufactured solution as the boundary data), and initializes the solution and right-hand side vectors.
We have to assemble the local matrix and map the local matrix and vector on each cell into the global system. Here is where we did the most modification, to include additional for loops that handle the integral over the boundary \(\Gamma\) and to implement the tangential gradient. Hence we have additional code that not only allocates quadrature and FEValues objects for the cell terms, but also for face terms.
We also use extractors here to get the components of the vector-valued shape functions. We basically use them as subscripts on the FEValues objects below, and after subscription they refer only to the scalar component at position 0 or the scalar component at position 1:
Finally, in this step we transfer the local contributions from cell_matrix and cell_rhs into the global matrix and right-hand-side vector, and similarly do it again for the loop corresponding to the integral over the boundary.
To solve the system that we assembled, we rely on the function SaddlePointProblem<dim>::solve(). Here we use an iterative solver, namely the Conjugate Gradient (CG) method. To tell the CG algorithm when to stop, we need a SolverControl object; the stopping criterion is
After that, we have several more functions, such as
SaddlePointProblem<dim>::output_results,SaddlePointProblem<dim>::output_exact_results,SaddlePointProblem<dim>::compute_error,SaddlePointProblem<dim>::write_convergence_table,SaddlePointProblem<dim>::run,which more or less do something straightforward and are documented in detail in their respective source listings; see step-6.cc.
1.9.8