CFD Online Discussion Forums - Confusion about how the implicit coupled interfaces are implemented?

CFD Online Discussion Forums (https://www.cfd-online.com/Forums/)

- OpenFOAM Programming & Development (https://www.cfd-online.com/Forums/openfoam-programming-development/)

- - Confusion about how the implicit coupled interfaces are implemented? (https://www.cfd-online.com/Forums/openfoam-programming-development/196965-confusion-about-how-implicit-coupled-interfaces-implemented.html)

Confusion about how the implicit coupled interfaces are implemented?

Thinking about the matrix

is a discretization of linearized operator acted on state vector

, there result is

, which must be converged to source term

.

Considering parallism using domain decomposition method. x and b are divided as

and

and

are located and modifiable in processor j, and the access of

and

on another processor i needs communications.

Now, assume the matrix A can be decomposed as
$A = \begin{bmatrix} A_{11} & A_{12} & A_{13}\\ A_{21} & A_{22} & A_{23}\\ A_{31} & A_{32} & A_{33} \end{bmatrix}$
So every time function `Amul()` is invoked in Krylov type solver. The

on processor i must be calculated as:

$b_i=A_{ii}\cdot x_i+\sum_{j\ne i}A_{ij}\cdot x_j=(L_{ii}+D_{ii}+U_{ii})\cdot x_i + C_{ii}\cdot x_i + \sum_{j\ne i}A_{ij}\cdot x_j \\ =b_i'+b_i'' + b_i'''$

In Foam::lduMatrix, $L_{ii},D_{ii},U_{ii}$ are stored as lower_, diag_ and upper_ scalarField. The $C_{ii}$ is the cyclic operator (though I do not know how it is implemented). And the off-diagonal block operators $A_{ij}, i\ne j$ must be implemented in processor interfaces.

So my question is: where is off-diagonal block $A_{ij}$ stored or implemented. It must be related to the action of processor j to processor i, which is related to a processor interface.

Now, considering a face F. There are several possibilities:
1. F is an internal face, which connect cell O numbered o and cell N numbered n on processor i, so there will be matrix coefficients $a_{on}, a_{no}$ in $A_{ii}$ (there might be some source term and diagonal term depending on the scheme used). This is the simplest case;
2. F is a boundary face, which only have an owner cell O numbered as o in processor i. It can only form source term

and diagonal term on $D_{ii}$ . This is the most normal boundary condition;
3. If face F is an processor interface, it will have two owner cells. One is owner cell

on processor i, another one is owner cell

on processor j. And there must be an off-diagonal term. I noticed that there is `initMatrixInterfaces`before LDU matrix*vector operation in `Amul()` function. And after that `updateMatrixInterfaces` is called.
I see the parameter list is :

Code:

initMatrixInterfaces|updateMatrixInterfaces(

add,

interfaceBouCoeffs|interfaceIntCoeffs,

interfaces,

psi,

result,

cmpt

)

I think `psi` is the state vector to be solved and `cmpt` is the direction to be solved if the state vector is not a scalar field. The `interfaces` is the list of
processor or coupled patches. The `interfaceBouCoeffs/interfaceIntCoeffs` is the most confusing part. I do not know what it is. I do not know why it is not stored in lduMatrix but in fvMatrix<Type>.
4. If face F is an cyclic-like interface, it will have two `owner` in one processor. So in the so called implicit coupling implementation of cyclicFvPatchField, the function `initInterfaceMatrixUpdate` is empty. However, I still cannot understand why it is implicit.
5. When face F is on a cyclic processor interface... I am still unable to dig into it.:confused:

I think the coupling implementation in openfoam is not 'zero halo'. In the terms of domain decomposition, pure 'zero halo' is like Dirichlet-Dirichlet coupling. It is not implicit.

I just figured out what is going on in coupled BC.
In openfoam, face addressing (owner, neighbour ) is different from lduAddressing (lower, upper). There are two implementation in OFv1706. One is `fvMeshLduAddressing`, another one is `fvMeshPrimitiveLduAddressing`.
1. in `fvMeshLduAddressing`, lower==owner[0:nInternalFaces], upper==neighbour;
2. in `fvMeshPrimitiveLduAddressing` which is used in overset, there can be more elements in lower and upper.

`fvMeshLduAddressing` is most commonly used. Because of LDU addressing, the function of lduMatrix is limited and it is designed to be like that. lduMatrix is a scalar matrix which only stores none zero coefficients in `upper_`, `lower_`, `diag_`. So the coupled BC cannot be represented in lduMatrix because of this limitation.

So the coupled BC is implementated like this: the coupled BC related coefficients are stored in fvMatrix<T> as `internalCoeffs_` and `boundaryCoeffs_`. And the type is FieldField<Field,T>.

Now all the coefficients can be accessed in fvMatrix<T>. And the coupled BC
(both cyclic BC and processor BC) can be handled in an unified manner.

However, I am still confusing about how fvVectorMatrix are solved in segrated or coupled manner. The code is very confusing.:confused:

Quote:

Originally Posted by chengdi (Post 675969)

I also got to this point on my own and reading the OpenFOAM book, however I want to know how exactly are boundaryCoeffs_ and internalCoeffs_ are handled for processor boundary.

I don't get the weight business with gradientInternalCoeffs, valueInternalCoeffs etc.

It is irrelevant for me, I exactly know the fvMatrix coefficients in terms of upper, diagnol and source, I just want to run the darn thing in parallel