[marek88]

I need some help to understand what is going on in fvMatrix::solve and lduMatrix::Amul.
I attach a complete case directory which solves

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:

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

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:

-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:

-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

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:

#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]

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]

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]

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]

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?

#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;

}

[lengweee]

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