CFD Online Discussion Forums

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 17:01.