Dirichlet BC in TwoScale method#
Symbolic work to verify how to deal with Dirichlet BC when passing from fine scale to coarse scale.
Attention
Approach (I) to (IV): Works only for shifted enriched function
initialize symbols#
The following symbolic variables represent the different matrices involved in the computation:
Fine scale system#
Applying the fine-scale Dirichlet operator to \(A\) and \(B\) gives:
TS approach (I)#
Applying the TS operator to \(A\),\(B\) gives:
Applying the coarse Dirichlet operator to \(AP_C\) and \(BP_C\) gives:
TS approach (II)#
Applying TS operator to \(AD\),\(BD\) gives:
Check same matrix#
Normally, if this approach is correct, matrices \(A_C^2\) and \(A_C\) could be the same. The difference gives
Let’s split \(DC\) into standard and enriched parts as follows:
This leads to:
As \(U=\mathbb{I}-D\), we have:
What can be said about \(D.P_S.D_{CS}\):
\(D_{CS}\) is filtering out columns (i.e. set them to zero) of \(P_S\) related to eliminated DOFs by coarse Dirichlet BC.
\(D\) is filtering out rows (i.e. set them to zero) of \(P_S.D_{CS}\) related to eliminated DOFs by fine Dirichlet BC.
De facto, as the imposed Dirichlet boundary conditions at both scales are expected to be the same, the following assertion will be true in most cases:
The eliminated DOFs at fine scale are included in the set of rows of the columns eliminated by \(D_{CS}\)
This is wrong if for example coarse and fine BC do not stop at the same location due to discretisation. In this case we expect to adapt coarse mesh so that it corresponds to this stopping location.
If this assertion is true, the rows eliminated by \(D\) are already null, thus:
\(D.P_S.D_{CS}=P_S.D_{CS}\)
and
\(U.P_S.D_{CS}=(\mathbb{I}-D).P_S.D_{CS}=P_S.D_{CS}-D.P_S.D_{CS}=P_S.D_{CS}-P_S.D_{CS}=0\)
The expression of \(A_C^2-A_C\) then becomes:
But it is clearly not the case for \(D.P_E.D_{CE}\) as coarse enriched eliminated DOFs if any are not a priori related to equivalent DOFs of \(D_{CS}\) so \(D.P_E.D_{CE}\neq P_E.D_{CE}\)
But what can be said about \(D.A.D-A\) considering that \(D=\mathbb{I}-U\)
The expression of \(A_C^2-A_C\) becomes:
And as already mentioned \(U.P_S.D_{CS}=0\) thus expression simplify further:
Which is not null as there is no reason for \(U.P_E\) to be null. So it cannot be simplified further except if we consider that enriched DOFs are not eliminated and in this case \(D_{CE}=\mathbb{I}\) and \(A_C^2-A_C\) becomes:
In conclusion, \(A_C^2\) and \(A_C\) are not the same.
check same system#
Matrices are not the same but maybe systems are giving the same solutions. From last expression of \(dA_C^2=A_C^2-A_C\) (the one with \(D_{CE}\) not forcefully \(\mathbb{I}\) we can write:
And if both systems give the same solution \(A_C^2.X_C\) should be equal to \(B_C^2\):
\(A_C^2.X_C=(A_C+dA_C^2).X_C=B_C+dA_C^2.A_C^{-1}.B_C=B_C^2\)
or \(B_C^2-B_C=dA_C^2.A_C^{-1}.B_C\)
we have, splitting \(X_{DC}\) into its two contributions \(X_{DCS}\) \(X_{DCE}\):
using \(D.P_S.D_{CS}=P_S.D_{CS}\) and \(U.P_S.D_{CS}=0\)
\(B_C^2-B_C\) is
But what is \(X_{DCE}\) ? It is the vector of potential imposed value applied to some enriched DOFs. If any enriched DOFs are imposed they will be used to eliminate rank deficiency in problem and will thus certainly be set to zero. So we can expect that if \(X_{DCE}\) exists it will be a null vector. So we will set this hypothesis:
\(X_{DCE}=0\)
With this hypothesis expression becomes:
On the other hand we have \(dA_C^2.A_C^{-1}.B_C\) which is not easy to manipulate. But if we could have expressed \(B_C^2-B_C\) as \(H.B_C\) then relation to pouve would be:
\(H.B_C=dA_C^2.A_C^{-1}.B_C\)
\((H-dA_C^2.A_C^{-1}).B_C=0\)
Then
if \(B_C=0\), then \(H\) does not exist
if \(B_C\neq0\) either \(H-dA_C^2.A_C^{-1}\) is orthogonal to \(B_C\) or \(H-dA_C^2.A_C^{-1}=0\)
In this last case
\(H=dA_C^2.A_C^{-1}\)
\(H.A_C=dA_C^2\)
not very simple to handle.
It is hard to conclude
TS approch (III)#
This approach is a simple reorganization of the approach (I) to optimize implementation. First we rewrite final system using block \(D_C\)
Then we set :
System becomes
We see that for \(A_C^3\) construction we don’t need to keep \(P_S\), \(P_E\) during the TS loop. Only \(Q_S\),\(Q_E\) are required. But for \(B_C^3\) we need \(P_E\), \(P_S\) and \(A\) !!!
But one can observe that \(P_S.X_{DCS}\) is constant during TS loop and can be set as a \(W_S\) vector. Thus
Is also constant during TS loop. Thus we can write \(B_C^3\) as:
Now, concerning \(A.P_E.X_{DCE}\) if we look only at this formal expression it has to be updated at all TS iteration as \(P_E\) changes. But if we follow the same hypotheses as above (\(X_{DCE}=0\)) the final system is:
From an implementation point of view, it is perfect as \(P_E\),\(P_S\) are not needed anymore and applying \(Q_S\),\(Q_E\) does the job of applying the TS operator and the coarse Dirichlet BC at the same time.
But in the TS Loop we project the solution \(X_C\) of the problem at coarse scale onto the fine-scale field using \(P_E\),\(P_S\) as follows:
Here we are going to use the fact that \(P.\mathbb{I}=P.(D_C+U_C)=Q+P.U_{C}\) so the complement of \(Q\) for \(P\) is:
So the projected solution is:
And as \(X_C\) filtered by \(U_C\) in fact gives \(X_{DC}\) as they are related to the same DOFs, we have \(L.X_C=P.U_C.X_C=P.X_{DC}=[W_S+P_E.X_{DCE}]\)
and as we consider above that \(X_{DCE}=0\)
\(L.X_C=[W_S+0]\)
From an implementation point of view it is again perfect as \(P_E\),\(P_S\) are not needed anymore and only \(Q_S\),\(Q_E\) and \(W_S\) are required
TS approach (IV)#
This approach is motivated by the use of PETSc nested matrix format. In this case sub-blocks are treated independently. The idea is to consider that no enriched DOFs are eliminated and to leverage matrix manipulation by really eliminating Dirichlet boundary condition from the system. Also, the enriched space is now really limited to only the enriched DOFs of enriched nodes (i.e. no Dirichlet). This approach is thus a reorganization of the approach (I) like approach (III). First we add the following definition:
The standard set is split into Dirichlet-imposed DOFs and free DOFs called hereafter reduced DOFs
The enriched set only has DOFs related to enriched nodes (no Dirichlet) also called reduced DOFs hereafter
\(R\) an operator to passe from the reduced set to the full set
And similarly to approach (III) we are going to mix boundary condition and operator. This is done here by splitting \(P_S\) in its Dirichlet column block \(P_{SD}\) and its reduced column block \(P_{SR}\). For enriched operator we use again the name \(P_E\) but now columns are in this approach, restricted only to enriched nodes (no Dirichlet):
The solution at coarse level is:
Using the first step of approach (I) we have:
And now applying the reduction operator:
Like in approach (III) we can set
which are constant during TS loop and can be computed once during the initialization of the loop. The system is then expressed as:
And regarding projection on fine computation we have:
So like in approach (III) it is possible to mix boundary condition and operator application and get the following minimal computations :
create \(P_E\)
create \(P_{SR}\) and temporary \(P_{SD}\)
compute \(W\)
remove \(P_{SD}\)
compute once at the initialization of the loop:
\(Z\)
\(P_{SR}^t.(B-Z)\) block of the VecNest vector \(B_C^4\)
\(P_{SR}^t.A.P_{SR}\) block of the MatNest matrix \(A_C^4\)
In the loop at each TS iteration:
compute \(P_E\)
compute each block of the MatNest matrix \(A_C^4\) related to \(P_E\)
compute \(P_{E}^t.(B-Z)\) block of the VecNest vector \(B_C^4\)
Solve the system
compute \(STS\)
TS approch (V)#
This approach is just the generalization of approach (IV) to a general enriched function. Imposing Dirichlet is now done on the combination of the standard and enriched DOF as enriched function is not anymore null at enriched nodes. This corresponds to the addition of kinematic equation of the form:
\(X_{D}-G_E.X_{RE}=X_{RD}\)
where
\(X_{D}\) are the standard DOFs involved in a kinematic equation (i.e. DOFs on enriched node where we impose Dirichlet BC)
\(G_E\) is an operator constructed from enriched function values:
rows correspond to the set of DOFs related to \(X_D\) (standard DOFs eliminated)
columns correspond to the set of enriched DOFs related to \(X_D\)
\(X_{RD}\) remains the same: the imposed Dirichlet values
Note
These extra equations exist only if the Dirichlet boundary condition is applied on an enriched node. Otherwise it is a usual BC elimination. This can be imposed by setting a zero row in \(G_E\) for these non enriched nodes. If all Dirichlet boundary conditions are imposed on non enriched nodes \(G_E\) is null and consequently there is no difference with approach (IV).
These extra equations can be added naturally to the restriction operator \(R\) by eliminating \(X_{D}\) which gives
The solution at coarse level is given by:
where we recognize the expression for \(X_D\) on the first line.
The system and the projection are:
Let add the following:
Then the system and the projection are:
From a formal point of view this approach is the same as approach (IV) with \(T_E\) taking the role of \(P_E\).
Let look at \(T_E\).It is \(P_E\) plus \(P_{SD}G_E\). \(P_{SD}\) are the column of \(P_S\) related to Dirichlet DOFs. Multiplied by \(G_E\) it modify column of \(P_E\) related to enriched DOFs on node where Dirichlet BC are applied. And more precisely this operation correspond in fact to use shifted enriched function for those particular DOFs/column of \(P_E\)
Thus the computational path is almost the same as the one used in approach (IV):
create \(P_E\) and identify column of \(P_E\) related to node with Dirichlet BC
create \(P_{SR}\) and temporary \(P_{SD}\)
compute \(W\)
remove \(P_{SD}\)
compute once at the initialization of the loop:
\(Z\)
\(P_{SR}^t.(B-Z)\) block of the VecNest vector \(B_C^5\)
\(P_{SR}^t.A.P_{SR}\) block of the MatNest matrix \(A_C^5\)
In the loop at each TS iteration:
compute \(T_E\) (like \(P_E\) but with shifted enriched function for Diriclet column)
compute each block of the MatNest matrix \(A_C^5\) related to \(T_E\)
compute \(T_{E}^t.(B-Z)\) block of the VecNest vector \(B_C^5\)
Solve the system
compute \(STS\)
TS approch (VI)#
In approach (IV) and (V) elimination of Dirichlet boundary condition by restriction operator leads to a system with sizes that are no longer a multiple of the block size (bs=nb components per node). To keep this property, as in conventional Dirichlet treatment in FEniCSx/PETSc, using a null row/column operator plus associated identity addition can do the job. This reactivates approach (II) as now we can imagine to use fine system with Dirichlet elimination to generate blocks related to standard DOFs.
Let’s start by renaming everything considering general enriched function and potentially imposed Dirichlet boundary condition to enriched nodes. This last point adds the following equations:
\(X_{SD}-G_E.X_{E}=C_{SD}\)
or
\(X_{SD}=G_E.X_{E}+C_{SD}\)
with:
So conventionally, we define a \(T\) kinematic operator as follows:
And the coarse solution \(X_C\) vector can be expressed as:
where we recognize the expression of \(X_{SD}\) on the first line.
Conventionally the system
\(A_{CC}.X_C=B_C\)
becomes
\(A_{CC}.(T.X_I+X_H)=B_C\)
or
\(A_{CC}.T.X_I=B_C-A_{CC}.X_H\)
And by pre multiplying by \(T^t\) to re obtain a square matrix:
\(T^t.A_{CC}.T.X_I=T^t.(B_C-A_{CC}.X_H)\)
But as mentioned above in this approach we don’t want to have a reduced size system. Thus we embedded this system into a larger one by simply adding dummy equations corresponding to constrained set with the addition of dummy unknowns \(X_{Dumy}\).
\(T\) becomes
And the solution is then given by:
The system is then \(A_{CC}.(T.X_{CSDdumy}+X_H)=B_C\)
or
\(A_{CC}.T.X_{CSDdumy}=B_C-A_{CC}.X_H\)
And by premultiplying by \(T^t\) to obtain a square matrix again:
\(T^t.A_{CC}.T.X_{CSDdumy}=T^t.(B_C-A_{CC}.X_H)\)
But now the sub block \(SD\times SD\) is null so we need to add 1 to the diagonal of this sub block to make it invertible. We call \(U_C\) like in previous approach the matrix containing this modified sub block :
The final form of the solved system is then:
\((T^t.A_{CC}.T+U_C).X_{CSDdumy}=T^t.(B_C-A_{CC}.X_H)\)
And the coarse solution is:
\(U_C= T.(T^t.A_{CC}.T+U_C)^{-1}.T^t.(B_C-A_{CC}.X_H)+X_H\)
When \(G_E\) is null (with shift enrichment function it can be the case) then these two steps can be grouped in one with the solution of the system being directly the searched coarse solution. \(T\) is then simply \(D_C\) defined in previous approach:
\((D_C^t.A_{CC}.D_C+U_C).X_{C}=D_C^t.(B_C-A_{CC}.X_H)+X_H\)
Let’s now split TS operator considering constraints at the coarse level and Dirichlet at fine level.
We have
Note that imposed constraint \(C_{SD}\) correspond exactly to Dirichlet imposed value at fine scale on common coarse location. And thus by construction \(PS_{DI}\) is null and \(PS_{DD}\) is the identity as only \(D\) set coarse DOFs are related by a coefficient of 1 with the \(D\) set fine DOFs. Here we also consider that fine Dirichlet boundary condition are embedded in coarse Dirichlet BC so d set DOFs are necessarily related only to D set coarse DOFs and then \(PS_{dI}\) is null. So \(P\) simplifies to:
To simplify for now we group TS operator in standard and enriched part:
And split the same way the fine system:
The Dirichlet boundary condition operators at fine scale are:
Thus the fine system with eliminated Dirichlet is:
Let’s now apply the TS operator to both systems:
and apply the coarse constraints \(T\) to form the system
dA is not null as seen in approach (II) but the (SD+SI)x(SD+SI) block is.
dB is not null but SDx1 is. And if we further impose that d set fine boundary condition is a linear interpolation of coarse D then \(C_d=PS_{dD}.C_{SD}\). Then it is in this case the(SD+SI)x1 which is null.
The question is then: can we mix use of \(A\),\(B\) and \(AD\),\(BD\) per block to get directly \(A_C^6\) correctly without adding \(U_C\) which impose extra matrix operation with PETSc.
Not clear for now ……
Let us \(Q\),\(W\),\(Z\) like in previous approaches:
And lets imagine the modified \(Q\):
Then we can obtain a corect sub block (stdxstd) with this operator:
But it cost the price of having \(Q\) and \(Qb\) in memory, which is too much, because \(Qb\) is ok for sub block (stdxstd) but not for the other sub block.
It is cheaper to use \(U_C\) even if adding it will cost some reallocation but as it is done only once at the beginning of the loop it can be ok.
Conlusion this approch is not giving anything.