# boundary conditions in fvMatrix multiplication make me despair

 Register Blogs Members List Search Today's Posts Mark Forums Read

August 31, 2011, 11:16
boundary conditions in fvMatrix multiplication make me despair
#1
New Member

Marek
Join Date: Aug 2011
Posts: 6
Rep Power: 7
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;
}
Attached Files
 testFoam.tar.gz (1.8 KB, 32 views)

 September 1, 2011, 05:00 #2 Senior Member   Laurence R. McGlashan Join Date: Mar 2009 Posts: 370 Rep Power: 15 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. __________________ Laurence R. McGlashan :: Website

 September 1, 2011, 05:57 #3 New Member   Marek Join Date: Aug 2011 Posts: 6 Rep Power: 7 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

 September 1, 2011, 06:05 #4 Senior Member   Laurence R. McGlashan Join Date: Mar 2009 Posts: 370 Rep Power: 15 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. __________________ Laurence R. McGlashan :: Website

 September 1, 2011, 07:40 #5 New Member   Marek Join Date: Aug 2011 Posts: 6 Rep Power: 7 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 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="<

 Thread Tools Display Modes Linear Mode

 Posting Rules You may not post new threads You may not post replies You may not post attachments You may not edit your posts BB code is On Smilies are On [IMG] code is On HTML code is OffTrackbacks are On Pingbacks are On Refbacks are On Forum Rules

 Similar Threads Thread Thread Starter Forum Replies Last Post Tom Main CFD Forum 7 August 26, 2014 05:58 Thomas P. Abraham Main CFD Forum 20 July 7, 2013 05:05 NickolasPl OpenFOAM Programming & Development 2 May 19, 2011 05:37 ttdtud OpenFOAM Installation 46 May 5, 2009 02:32 Tom Main CFD Forum 0 November 5, 2002 02:54

All times are GMT -4. The time now is 12:52.