[marek88]

I have a problem with unexpected results of the matrix*vector multiplication fvMatrix<scalar>::Amul().

The code below is supposed to calculate the laplacian of the field psi, where psi=y*y is the square of the y-coordinate.
Because d/dy d/dy (y*y) == 2 , I expect to calculate result==2, but result varies between -1e-5 and 5e-7 on a 20x20x1 mesh. The result is uniform except for one layer next to a wall.

I think there is a serious flaw about how I call the Amul() function. Could you please check for an obvious mistake?
(The code can be run for example on the mesh from the icoFoam/cavity tutorial. I have defined fixedValue boundary conditions for result.
I would like to solve this issue, because it is needed to verify parts of my CFD code.)

Thanks
Marek

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

    // define psi=y*y
    volScalarField psi( IOobject("psi",runTime.timeName(), mesh, IOobject::NO_READ, IOobject::AUTO_WRITE ),
    Foam::sqr(mesh.C().component(1)));
    // result=laplacian(psi) will be calculated below
    volScalarField result( IOobject("result",runTime.timeName(), mesh, IOobject::MUST_READ, IOobject::AUTO_WRITE ),
    mesh);

    //required to get the 3rd argument to Amul
    FieldField<Field, scalar> bouCoeffs(psi.boundaryField().interfaces().size());
    forAll(bouCoeffs,I){
      bouCoeffs.set(I,result.boundaryField()[I].valueBoundaryCoeffs(result));
    }

    fvMatrix<scalar> M=fvm::laplacian(result);
    M.Amul(result,psi,bouCoeffs,result.boundaryField().interfaces(),0); // should give result=laplacian(psi)

    runTime++;
    runTime.writeNow();
}