CFD Online Discussion Forums

CFD Online Discussion Forums (
-   OpenFOAM Programming & Development (
-   -   fvMatrix::Amul() usage (

marek88 August 19, 2011 10:16

fvMatrix::Amul() usage
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.)



#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 ),
    // result=laplacian(psi) will be calculated below
    volScalarField result( IOobject("result",runTime.timeName(), mesh, IOobject::MUST_READ, IOobject::AUTO_WRITE ),

    //required to get the 3rd argument to Amul
    FieldField<Field, scalar> bouCoeffs(psi.boundaryField().interfaces().size());

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


All times are GMT -4. The time now is 22:56.