CFD Online Logo CFD Online URL
www.cfd-online.com
[Sponsors]
Home > Forums > OpenFOAM

boundary conditions in fvMatrix multiplication make me despair

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

Like Tree1Likes
  • 1 Post By marek88

Reply
 
LinkBack Thread Tools Display Modes
Old   August 31, 2011, 11:16
Default boundary conditions in fvMatrix multiplication make me despair
  #1
New Member
 
Marek
Join Date: Aug 2011
Posts: 6
Rep Power: 5
marek88 is on a distinguished road
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
File Type: gz testFoam.tar.gz (1.8 KB, 29 views)
marek88 is offline   Reply With Quote

Old   September 1, 2011, 05:00
Default
  #2
Senior Member
 
Laurence R. McGlashan
Join Date: Mar 2009
Posts: 370
Rep Power: 14
l_r_mcglashan will become famous soon enough
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
l_r_mcglashan is offline   Reply With Quote

Old   September 1, 2011, 05:57
Default
  #3
New Member
 
Marek
Join Date: Aug 2011
Posts: 6
Rep Power: 5
marek88 is on a distinguished road
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
marek88 is offline   Reply With Quote

Old   September 1, 2011, 06:05
Default
  #4
Senior Member
 
Laurence R. McGlashan
Join Date: Mar 2009
Posts: 370
Rep Power: 14
l_r_mcglashan will become famous soon enough
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
l_r_mcglashan is offline   Reply With Quote

Old   September 1, 2011, 07:40
Default
  #5
New Member
 
Marek
Join Date: Aug 2011
Posts: 6
Rep Power: 5
marek88 is on a distinguished road
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<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;

}
fumiya likes this.
marek88 is offline   Reply With Quote

Reply

Thread Tools
Display Modes

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 Off
Trackbacks are On
Pingbacks are On
Refbacks are On


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.5-dev ttdtud OpenFOAM Installation 46 May 5, 2009 02:32
Boundary conditions? Tom Main CFD Forum 0 November 5, 2002 02:54


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