CFD Online Discussion Forums

CFD Online Discussion Forums (http://www.cfd-online.com/Forums/)
-   OpenFOAM Programming & Development (http://www.cfd-online.com/Forums/openfoam-programming-development/)
-   -   fvMatrix::Amul() usage (http://www.cfd-online.com/Forums/openfoam-programming-development/91699-fvmatrix-amul-usage.html)

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.)

Thanks
Marek

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 ),
    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();
}



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