 |
G+Smo
23.12.0
Geometry + Simulation Modules
|
|
|
G+Smo
23.12.0
Geometry + Simulation Modules
|
|
We consider the Poisson boundary value problem (for simplicity here only with Dirichlet conditions):
\begin{eqnarray*} -\Delta u &=& f \quad \text{in} \quad \Omega, \\ u &=& g \quad \text{on} \quad \partial\Omega. \end{eqnarray*}
Its variational form is: Find \( u\in V_g \) such that
\[ \int_\Omega \nabla u \nabla v d x = \int_\Omega f v dx \]
for all \( u\in V_0 \), which corresponds to the linear system
\[ A \underline u = \underline f. \]
The computational domain is a non-overlapping multi-patch domain:
\[ \overline{ \Omega } = \bigcup_{k=1}^K \overline{ \Omega_k } \]
\[ \Omega_k \cap \Omega_l = \emptyset \quad \text{ for all } k\not=l. \]
On each patch, we have a isogeometric function space \( V_k \). We assume that the function spaces are fully matching. The global function space is
\[ V_h = \{ v\in V \;:\; v|_{\Omega_k} \in V_k \} \]
where \( V \) is the proper continuous space.
After reading the command line arguments, we read in the geometry file and modify the geometry as desired:
Then, we set up the right-hand-side f and a gsBoundaryConditions object based on the command line arguments.
Then, we extract the basis from the multi patch object to obtain a multi basis.
We set the spline degree and the grid size as desired.
Now, we can discuss the IETI system.
We are considering patch-wise formulations of the variational formulation:
\[ \int_{\Omega_k} \nabla u \nabla v d x = \int_{\Omega_k} f v dx \]
Each of these bilinear forms and linear forms lead to matrices and vectors
\[ A_k \qquad\text{and}\qquad \underline f_k. \]
Using the jump matrices \( B_k \), which are provided by the class gsIetiMapper, we obtain the following problem:
\[ \begin{pmatrix} A_1 & & & & B_1^\top \\ & A_2 & & & B_2^\top \\ & &\ddots& & \vdots \\ & & & A_K & B_K^\top \\ B_1 & B_2 &\cdots& B_K & 0 \\ \end{pmatrix} \begin{pmatrix} \underline u_1 \\ \underline u_2 \\ \vdots \\ \underline u_K \\ \underline \lambda \end{pmatrix} = \begin{pmatrix} \underline f_1 \\ \underline f_2 \\ \vdots \\ \underline f_K \\ 0 \end{pmatrix} \]
This is the kind of systems, the class gsIetiSystem is working with. The problem is that – in general – the matrices \( A_k \) are singular. To circumvent this problem, in FETI-dp / IETI-dp methods, we use primal degrees of freedom, like vertex values or edge averages.
These primal dofs are handled by the classes gsPrimalSystem and again gsIetiMapper. The IETI mapper provides vectors \( c_k^{(i)} \) that realize the constraint that the primal dof is set to zero. For each patch, these primal constraints build a matrix \( C_k \). The IETI mapper also provides a vector of indices that associates each of these primal constraints to a (global) index referring to the primal dof.
Finally, we obtain one subspace (the K+1st) that refers to the primal degrees of freedom.
We have the following IETI system:
\[ \begin{pmatrix} \tilde A_1 & & & & \tilde B_1^\top \\ & \tilde A_2 & & & \tilde B_2^\top \\ & & \ddots & & \vdots \\ & & & \tilde A_{K+1} & \tilde B_{K+1}^\top \\ \tilde B_1 & \tilde B_2 & \cdots & \tilde B_{K+1} & 0 \\ \end{pmatrix} \begin{pmatrix} \underline{ \tilde u}_1 \\ \underline{ \tilde u}_2 \\ \vdots \\ \underline{ \tilde u}_{K+1} \\ \underline \lambda \end{pmatrix} = \begin{pmatrix} \underline{ \tilde f}_1 \\ \underline{ \tilde f}_2 \\ \vdots \\ \underline{ \tilde f}_{K+1} \\ 0 \end{pmatrix}, \qquad (*) \]
which is passed to the class gsIetiSystem. For \( k=1,\ldots,K \), we have the local saddle point system
\[ \tilde{A}_k = \begin{pmatrix} A_k & C_k^\top \\ C_k & 0 \end{pmatrix}, \]
the extended jump matrix
\[ \tilde{B}_k = B_k (P_0^{(k)})^\top = \begin{pmatrix} B_k & 0_k \end{pmatrix} \]
and the extended right-hand side
\[ \underline{\tilde f}_k = P_0^{(k)} \underline f_k = \begin{pmatrix} \underline f_k \\ 0 \end{pmatrix}. \]
These modifications are provided by gsPrimalSystem. The class gsPrimalSystem moreover computes the energy minimizing basis functions and represents them as matrices \( \Psi_k = (P_c^{(k)})^\top \tilde{A}_k^{-1} R_c \), where \( R_c=(0,I)^\top \) and \( P_c=(I,0) \). It then computes
\[ \tilde A_{K+1} = \sum_{k=1}^K \Psi_k A_k \Psi_k^\top \quad\text{and}\quad \tilde B_{K+1} = \sum_{k=1}^K B_k \Psi_k^\top \quad\text{and}\quad \underline{\tilde f}_{K+1} = \sum_{k=1}^K \Psi_k \underline f_k. \]
By setting a corresponding flag, it is possible to eliminate the pointwise constraints. A constraint is pointwise the correspoding row of \( C_k \) has only one non-zero entries. Typically, the corner values are pointwise constraints.
Instead of directly solving the IETI-saddle point system (*), we solve the corresponding Schur complement system
\[ F \underline \lambda = \underline g, \]
where
\[ F = \sum_{k=1}^{K+1} \tilde B_k \tilde A_k^{-1} \tilde B_k^\top \quad\text{and}\quad \underline g = \sum_{k=1}^{K+1} \tilde B_k \tilde A_k^{-1} \underline{\tilde f}_k. \]
The Schur complement system is preconditioned with the scaled Dirichlet preconditioner
\[ M_{sD} = \sum_{k=1}^{K} \hat B_k D_k^{-1} S_k D_k^{-1} \hat B_k^\top, \]
where \( S_k \) is the Schur complement of \( A_k \) (not: \(\tilde A_k\)) with respect to the dofs on the interfaces. \( \hat B_k \) is the restriction of \( B_k \) to the dofs on the interfaces. \( D_k \) is a scaling matrix, e.g., based on multiplicity scaling. The class gsScaledDirichletPrec not only collects the contributions of the preconditioner; it also helps with restricting to the skeleton and computing the Schur complement.
In the example file ieti2_example.cpp, we solve the saddle point system (*) with a MINRES solver, preconditioned with a block-diagonal preconditioner, where the last block is \( M_{sD} \).
First, we set up the gsIetiMapper object, which helps us with assembling and keeping track of the dof indices.
For initializing the IETI mapper, it is necessary to provide a gsDofMapper for the global problem and the function values for the eliminated dofs (typically the Dirichlet conditions). If systems are considered, the method it is expected to have one IETI mapper per variable.
For obtaining a gsDofMapper, we use the assembler:
Now we set the primal degrees of freedom. The following commands compute the (vectors that make up the) constraint matrices \( C_k \) and store them in the gsIetiMapper object.
The following command computes the jump matrices \( B_k \) and stores them in the gsIetiMapper object.
Now, we set up the classes for the IETI system, the primal system, and the scaled Dirichlet preconditioner.
Now, we assemble the linear system and the right-hand side on each patch using a for loop. We use local versions of the geometry (gsMultiPatch), the basis (gsMultiBasis) and the boundary conditions (gsBoundaryConditions).
The function initFeSpace provides a new dof mapper and the Dirichlet data. This is necessary since it might happen that a 2d-patch touches the Dirichlet boundary just with a corner or that a 3d-patch touches the Dirichlet boundary with a corner or an edge. These cases are not covered by bc.getConditionsForPatch
The example file ieti2_example.cpp shows how to assemble with a classical assembler.
Now, we tell the preconditioner about the local matrix \( A_k \), the local vector \( \underline f_k \) and the matrix \( B_k \). This is done before the class gsPrimalSystem modifies these objects. The following command does all necessary steps:
This can also be done in a step-by-step way, for example if one does not want to use a sparse Cholesky solver (which requires symmetric positive definite problems and would be done by default). This is shown in ieti2_example.cpp.
As next step, we treat the primal dofs. We first compute energy minimizing primal basis functions. Then the following commands amends the variables jumpMatrix, localMatrix, and localRhs, which hold \( B_k \), \( A_k \) and \( \underline f_k \) before calling the function and \( \tilde B_k \), \( \tilde A_k \), and \( \underline{\tilde f}_k \) thereafter. Finally, this function collects the contributions for \( \tilde B_{K+1} \), \( \tilde A_{K+1} \), and \( \underline{\tilde f}_{K+1} \) in the object.
This can also be done in a step-by-step way, for example if one does not want to use a sparse LU solver (which would be done by default). This is shown in the file ieti2_example.cpp.
Now, we add \( \tilde B_k \), \( \tilde A_k \), and \( \underline{\tilde f}_k \) to the IETI-system:
Additionally, one can also provide a solver object that realizes \( \tilde A_k^{-1} \), possibly the one which has already been used for the setup of the bases for the primal doefs. By default, again, a LU solver would be used.
Finally, we add \( \tilde B_{K+1} \), \( \tilde A_{K+1} \), and \( \underline{\tilde f}_{K+1} \) to the IETI-system:
The gsScaledDirichletPrec can automatically compute the multiplicity scaling. Other scaling matrices can be manually defined.
Now, we compute the right-hand-side
\[ \underline g = \sum_{k=1}^{K+1} \tilde B_k \tilde A_k^{-1} \underline{\tilde f}_k. \]
We choose a random initial guess.
As a next step, we solve the problem
\[ F \underline \lambda = \underline g \]
using the scaled Dirichlet preconditioner with a conjugate gradient solver.
So far, we have only computed the Lagrange multipliers \( \underline \lambda \), from which we want to recover the solution.
Here, we first compute the solutions \( \underline{\tilde u}_k \) for \( k=1,\ldots,K+1 \). Here, the first K solutions refer to the patches and have homogenous values for the primal dofs. So, we distribute the primal dof values from the primal problem (= last subdomain) to the patches (=first K subdomains) and obtain the solutions \( \underline u_k \) for \( k=1,\ldots,K \). Then, finally, the IETI mapper is able to combine everything into one solution vector \( \underline u \).
Here is the full file examples/ieti_example.cpp. Clicking on a function or class name will lead you to its reference documentation.