CFD Online Discussion Forums (http://www.cfd-online.com/Forums/)
-   OpenFOAM (http://www.cfd-online.com/Forums/openfoam/)
-   -   boundary conditions in fvMatrix multiplication make me despair (http://www.cfd-online.com/Forums/openfoam/92046-boundary-conditions-fvmatrix-multiplication-make-me-despair.html)

 marek88 August 31, 2011 11:16

boundary conditions in fvMatrix multiplication make me despair

1 Attachment(s)
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; }```

 l_r_mcglashan September 1, 2011 05:00

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.

 marek88 September 1, 2011 05:57

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

 l_r_mcglashan September 1, 2011 06:05

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.

 marek88 September 1, 2011 07:40

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<scalar> 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="<<nl<<x.internalField()<<nl<<endl;     Info<<"A*x="<<nl<<b.internalField()<<nl<<endl; }```

 All times are GMT -4. The time now is 19:46.