# boundary conditions in fvMatrix multiplication make me despair

 Register Blogs Members List Search Today's Posts Mark Forums Read

August 31, 2011, 11:16
boundary conditions in fvMatrix multiplication make me despair
#1
New Member

Marek
Join Date: Aug 2011
Posts: 6
Rep Power: 14
I need some help to understand what is going on in fvMatrix::solve and lduMatrix::Amul.
I attach a complete case directory which solves
Code:
`fvm::laplacian(x)==b`
for b=(1 1 1 1 1), x unknown, on a mesh of 5 cells making up the unit cube.

The output of the code is:
Quote:
 x= 5(-0.05 -0.11 -0.13 -0.11 -0.05) A*x= 5(-0.3 0.2 0.2 0.2 -0.3)
The result x is what I expect from the laplace equation.
However, for A*x I expect
Quote:
 A*x = (0.2 0.2 0.2 0.2 0.2) instead of A*x = (-0.3 0.2 0.2 0.2 -0.3)
Obviously the boundary conditions are not included in A.
If I display the matrix elements of A, I get:
Quote:
 -5.0 5.0 0.0 0.0 0.0 5.0 -10. 5.0 0.0 0.0 0.0 5.0 -10. 5.0 0.0 0.0 0.0 5.0 -10. 5.0 0.0 0.0 0.0 5.0 -5.0
but (given Dirichlet boundary conditions) A should look like:
Quote:
 -15. 5.0 0.0 0.0 0.0 5.0 -10. 5.0 0.0 0.0 0.0 5.0 -10. 5.0 0.0 0.0 0.0 5.0 -10. 5.0 0.0 0.0 0.0 5.0 -15.
May I ask you to compile and run the attached code (wmake && blockMesh && testFoam) and if possible suggest how I have to change the line
Code:
`A.Amul(b,x,interfaceBouCoeffs,interfaces,0);`
such that the expected b=(0.2 0.2 0.2 0.2 0.2) is reproduced. This simple step is really important for my work, because I cannot continue if a trivial matrix multiplication results in unexplained values.

Thank you very much
Marek

This is the solver code, also included in the attachment:
Code:
```#include "fvCFD.H"
#   include "FieldField.H"

int main(int argc,char *argv[]){
#   include "setRootCase.H"
#   include "createTime.H"
#   include "createMesh.H"

volScalarField x( IOobject("x",runTime.timeName(), mesh, IOobject::MUST_READ, IOobject::AUTO_WRITE ),
mesh);
volScalarField b( IOobject("b",runTime.timeName(), mesh, IOobject::NO_READ, IOobject::AUTO_WRITE ),
mesh,1.);

fvMatrix<scalar> A=fvm::laplacian(x);
solve(A==b);

const FieldField< Field, scalar > interfaceBouCoeffs(0);
const lduInterfaceFieldPtrsList interfaces(0);
A.Amul(b,x,interfaceBouCoeffs,interfaces,0);

Info<<"x="<<nl<<x.internalField()<<nl<<endl;
Info<<"A*x="<<nl<<b.internalField()<<nl<<endl;
}```
Attached Files
 testFoam.tar.gz (1.8 KB, 40 views)

 September 1, 2011, 05:00 #2 Senior Member   Laurence R. McGlashan Join Date: Mar 2009 Posts: 370 Rep Power: 23 Hi (I like the thread title by the way!), Matrix A does have the 'effect' of the boundary on the internal cells -> try printing out A.internalCoeffs(). You could add the internalCoeffs to the relevant cells before using Amul, but I'm not sure what you're trying to do to be honest, there's probably a 'better' way of going about what you're actually trying to achieve. __________________ Laurence R. McGlashan :: Website

 September 1, 2011, 05:57 #3 New Member   Marek Join Date: Aug 2011 Posts: 6 Rep Power: 14 Hello Laurence, I am implementing a new discretisation scheme, and I need to verify where things go wrong. I narrowed the problem down to the boundary and internal coefficients of boundary patches. As you suggested, I am now adding the contributions from A.boundaryCoeffs() and A.internalCoeffs() to the matrix multiplication. It is somewhat involved, because I loop over each patch, then over each face, and finally have to get the face-to-cell addressing right. Do you know an OpenFOAM function that adds the boundary contributions to my solution b? I tried to find the corresponding code in the iterative solvers (PCG.C, PBiCG.C), but even there I was not able to find out how the boundary contributions are taken care of. Thanks Marek

 September 1, 2011, 06:05 #4 Senior Member   Laurence R. McGlashan Join Date: Mar 2009 Posts: 370 Rep Power: 23 Hmmm, well if you look at the protected functions in fvMatrix, there are ones called "addToInternalField()" etc. I assume they do what you want. I'm afraid that's about my limit! Maybe grep those functions to see where they are called. __________________ Laurence R. McGlashan :: Website

 September 1, 2011, 07:40 #5 New Member   Marek Join Date: Aug 2011 Posts: 6 Rep Power: 14 If anyone is interested, I have added a nested for-loop after the A.Amul() multiplication which solves the problem. Is it possible to replace the inner loop with some existing function instead of multiplying element by element? Code: ```#include "fvCFD.H" # include "FieldField.H" int main(int argc,char *argv[]){ # include "setRootCase.H" # include "createTime.H" # include "createMesh.H" volScalarField x( IOobject("x",runTime.timeName(), mesh, IOobject::MUST_READ, IOobject::AUTO_WRITE ), mesh); volScalarField b( IOobject("b",runTime.timeName(), mesh, IOobject::NO_READ, IOobject::AUTO_WRITE ), mesh,1.); fvMatrix A=fvm::laplacian(x); solve(A==b); A.Amul(b,x,A.boundaryCoeffs(),x.boundaryField().interfaces(),0); forAll(b.boundaryField(),I){ const fvPatch &p=b.boundaryField()[I].patch(); forAll(p,J){ b[p.faceCells()[J]]+= x.boundaryField().boundaryInternalField()[I][J]*A.internalCoeffs()[I][J] +x.boundaryField()[I][J]*A.boundaryCoeffs()[I][J]; } } Info<<"x="<

 August 25, 2018, 03:55 #6 New Member   wei leng Join Date: Mar 2014 Posts: 3 Rep Power: 12 Please find the code of exporting openfoam matrix to MATLAB and checking that A*x = b in the following thread. Extracting the matrix data from an fvVectorMatrix