
[Sponsors] 
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: 6 
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 The output of the code is: Quote:
However, for A*x I expect Quote:
If I display the matrix elements of A, I get: Quote:
Quote:
Code:
A.Amul(b,x,interfaceBouCoeffs,interfaces,0); 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; } 

September 1, 2011, 05:00 

#2 
Senior Member
Laurence R. McGlashan
Join Date: Mar 2009
Posts: 370
Rep Power: 14 
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: 6 
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 facetocell 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: 14 
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: 6 
If anyone is interested, I have added a nested forloop 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<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; } 

Thread Tools  
Display Modes  


Similar Threads  
Thread  Thread Starter  Forum  Replies  Last Post 
Implementation of boundary conditions for FVM  Tom  Main CFD Forum  7  August 26, 2014 05:58 
Boundary Conditions  Thomas P. Abraham  Main CFD Forum  20  July 7, 2013 05:05 
OpenFOAM Variable Velocity Boundary Conditions  NickolasPl  OpenFOAM Programming & Development  2  May 19, 2011 05:37 
Installation OF1.5dev  ttdtud  OpenFOAM Installation  46  May 5, 2009 02:32 
Boundary conditions?  Tom  Main CFD Forum  0  November 5, 2002 02:54 