CFD Online Discussion Forums (http://www.cfd-online.com/Forums/)
-   OpenFOAM Programming & Development (http://www.cfd-online.com/Forums/openfoam-programming-development/)
-   -   Matrix manipulation (http://www.cfd-online.com/Forums/openfoam-programming-development/97251-matrix-manipulation.html)

 Linse February 13, 2012 11:23

Matrix manipulation

Hello all!

According to some standard books (Patankar, Versteeg, etc.), the common way of solving CFD-problems (abbreviatedly) consists of getting the equations into a form of aX + bX + cX = D, which is put into a matrix A of coefficients a, b and c and two vectors of X and D, basically giving a system [A] X = D.
Several steps are needed for solving this equation system, among others
- getting a, b and c into a form that makes [A] consisting of one diagonal of values instead of three (use of recurrence relations)
- calculating the actual value for X based on the results of the recurrence calculations

My question now is:
In which parts of the code are these steps done?
And how could these steps be manipulated?
Ideally I would like to access the equations/matrices on a cell-level and change/control them cell-wise...

I guess I will have to manipulate the matrices themselves, but how is that done? Any information/pointing to information would be highly welcome!

 andyru February 13, 2012 13:42

Hi,

if you look into the solver ico-Foam for example, the main coefficients of the matrix are calculated with function A():

Code:

volScalarField rUA = 1.0/UEqn.A();
Here you can acces them an manipulate them.

Code:

forAll(rUA,celli)
{
...
}

So, this depends, what you want to do?

I'm facing almost the same problem and I need to manipulate the source term grad(vsf) (http://www.cfd-online.com/Forums/ope...tml#post343643)

Greetings

Andy

 santiagomarquezd February 21, 2012 22:28

Quote:
 Originally Posted by Linse (Post 344151) Hello all! According to some standard books (Patankar, Versteeg, etc.), the common way of solving CFD-problems (abbreviatedly) consists of getting the equations into a form of aX + bX + cX = D, which is put into a matrix A of coefficients a, b and c and two vectors of X and D, basically giving a system [A] X = D. Several steps are needed for solving this equation system, among others - getting a, b and c into a form that makes [A] consisting of one diagonal of values instead of three (use of recurrence relations) - calculating the actual value for X based on the results of the recurrence calculations My question now is: In which parts of the code are these steps done? And how could these steps be manipulated? Ideally I would like to access the equations/matrices on a cell-level and change/control them cell-wise... I guess I will have to manipulate the matrices themselves, but how is that done? Any information/pointing to information would be highly welcome! Thank you all in advance!
The questions are not so clear. I could help but I can't understand the problem.

Regards.

 j-avdeev April 4, 2012 05:58

Propably, I have same problem.
I want to save matrices A to separate file (A from AX=D).

http://openfoamwiki.net/index.php/Op...x_coefficients
Quote:
 OpenFOAM stores: diagonal coefficients as a scalar field, indexed by cell volume; and off-diagonal coefficients as two scalar fields, indexed by face centres. lduAddressing is used to related the two indexing methods.
but how to write it in separate file?

 elvis April 4, 2012 07:08

I guess they want to do something similar that
Santiago Marquez Damian described in his post http://www.cfd-online.com/Forums/ope...tml#post302921

 santiagomarquezd April 4, 2012 08:52

Maybe you need http://openfoamwiki.net/index.php/Contrib_gdbOF, it has an specific tool to do that.

Regards

 j-avdeev April 23, 2012 17:27

Santiago, yes, GdbOF looks like what I want.
I've read GdbOF-Manual.pdf - it has a lot of, but very general examples like

Example 3 "View Boundary Field values with gdbOF"
Code:

\$ (gdb) ppatchvalues vSF 0
and not written exactly when, from which directory I have to type this command.
Maybe it's very simple for describing - but I need at least one step-by-step example.

I like example from (gdb)+OF
Code:

go to cavity tutorial directory
\$cd cavity
run icoFoam in debug
\$gdb icoFoam
create breakpoint
\$b icoFoam.C:77

\$run
There are two ways of stepping in the ﬁle:

\$n (next) will step to the next line in the current ﬁle, but will not go into functions and
included ﬁles.
\$s (step) will go to the next line, also in functions and included ﬁles.

- I want similar, but for GdbOF (with pfvmatrixfull or pfvmatrixsparse commands).

Can someone describe sequence of commands for describing? Or correct me, if I do something wrong.

 santiagomarquezd April 23, 2012 17:33

Well, it suppose some basic knowledge about gdb as a pre-requisite, tell us which solver do you want to debug and what system within it.

Regards

 j-avdeev April 24, 2012 16:11

Solver - let it be simpleFoam.
Which solver to use - not matter for me (at least I think so now).

What system within solver... == what algoritm of solver?

About simpleFoam exist article here The_SIMPLE_algorithm_in_OpenFOAM
but equations of solver far from Ax=b (it would be cool for me to know this moment, where equation like
Code:

fvm::div(phi, U) - laplacian(nu, U)
become single matrix like A or I wrong?)

Also if needs case for example of debugging - it could be
\$FOAM_TUTORIALS/incompressible/simpleFoam/pitzDaily

Thanks.

 santiagomarquezd April 24, 2012 16:28

Well,

Quote:
 Originally Posted by j-avdeev (Post 356727) it would be cool for me to know this moment, where equation like Code: fvm::div(phi, U) - laplacian(nu, U) become single matrix like A or I wrong?)
it's a long story about FVM discretization and C++ implementation, but in order to be helpful, suppose you're solving a tutorial for scalarTransportFoam and you want to see the equation for following lines in scalarTransportFoam.C (the scalar advection diffusion equation without sources)

Code:

00058            solve
00059            (
00060                fvm::ddt(T)
00061              + fvm::div(phi, T)
00062              - fvm::laplacian(DT, T)
00063            );

then go the case directory, for eample \$FOAM_TUTORIALS/basic/scalarTransportFoam/pitzDaily and start gdb

\$ gdb scalarTransportFoam

once you you're within gdb:

\$ start
\$ b fvMatrixSolve.C:140
\$ continue

these lines will show you the matrix for this system. You also can read the AdvDiff.dat in octave or Matlab as a matrix.

Regards.

 Lydia May 2, 2014 08:59

cannot write pvfmatrixfull with gdbof

Hello Santiago,
thanks for your tool. Sounds very promising!
However, while making first tests I came across the following problem:

I followed your description above. Using OpenFOAM_2.2.2_Debug and running the scalarTransportFoam pitzDaily tutorial with gdbof. When I want to write the pfvmatrixfull with:
\$(gdb) pfvmatrixfull this test.txt
\$(gdb) shell cat test.txt

The programme does not write the matrix but a file with the message:
No symbol "this" in current context.

Can anyone give me a hint, where my mistake is?

Regards,
Lydia

 nisi May 15, 2014 16:35

Dear Lydia,

Thanks for reporting.

We have solved that bug and also we have included new commands which improve our tool. You can download the updated version of gdbOF from the wiki page: http://openfoamwiki.net/index.php/Contrib_gdbOF

Regards, Juan

 Lydia May 20, 2014 13:05

problem not yet solved - cannot write pvfmatrixfull with gdbof_v1.03a

Dear Juan,

Maybe I'm doing something wrong...for making sure that is not the case I'm posting my complete shell protocol below.
PHP Code:

lydia@lydia:~/OpenFOAM\$ . ./OpenFOAM-2.2.2/etc/bashrc
lydia
@lydia:~/OpenFOAMcd run/tutorials/basic/scalarTransportFoam/pitzDaily/
lydia@lydia:~/OpenFOAM/run/tutorials/basic/scalarTransportFoam/pitzDailygdb scalarTransportFoam
GNU gdb
(Ubuntu 7.7-0ubuntu37.7
(C2014 Free Software FoundationInc.
This is free softwareyou are free to change and redistribute it.
There is NO WARRANTYto the extent permitted by law.  Type "show copying"
and "show warranty" for details.
This GDB was configured as "x86_64-linux-gnu".
Type "show configuration" for configuration details.
For
<
http://www.gnu.org/software/gdb/bugs/>.
Find the GDB manual and other documentation resources online at:
<
http://www.gnu.org/software/gdb/documentation/>.
For helptype "help".
Type "apropos word" to search for commands related to "word"...
(
gdbstart
Temporary breakpoint 1 at 0x4384a6
file /home/lydia/OpenFOAM/OpenFOAM-2.2.2/src/OpenFOAM/lnInclude/setRootCase.Hline 5.
Starting program
: /home/lydia/OpenFOAM/OpenFOAM-2.2.2/platforms/linux64GccDPDebug/bin/scalarTransportFoam

Temporary breakpoint 1
main (argc=1argv=0x7fffffffcfc8at /home/lydia/OpenFOAM/OpenFOAM-2.2.2/src/OpenFOAM/lnInclude/setRootCase.H:5
5        Foam
::argList args(argcargv);
(
gdbstart
The program being debugged has been started already
.
Start it from the beginning? (or ny
Temporary breakpoint 2 at 0x4384a6
file /home/lydia/OpenFOAM/OpenFOAM-2.2.2/src/OpenFOAM/lnInclude/setRootCase.Hline 5.
Starting program
: /home/lydia/OpenFOAM/OpenFOAM-2.2.2/platforms/linux64GccDPDebug/bin/scalarTransportFoam

Temporary breakpoint 2
main (argc=1argv=0x7fffffffcfc8at /home/lydia/OpenFOAM/OpenFOAM-2.2.2/src/OpenFOAM/lnInclude/setRootCase.H:5
5        Foam
::argList args(argcargv);
(
gdbb fvMatrixSolve.C:140
Breakpoint 3 at 0x44f6f1
fvMatrixSolve.C:140. (3 locations)
(
gdb) continue
Continuing.
/*---------------------------------------------------------------------------*\
| =========                 |                                                 |
| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
|  \\    /   O peration     | Version:  2.2.2                                 |
|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
|    \\/     M anipulation  |                                                 |
\*---------------------------------------------------------------------------*/
Build  2.2.2-9739c53ec43f
Exec
: /home/lydia/OpenFOAM/OpenFOAM-2.2.2/platforms/linux64GccDPDebug/bin/scalarTransportFoam
Date
May 20 2014
Time
18:53:06
Host
"lydia-Latitude-D630"
PID    3591
Case   : /home/lydia/OpenFOAM/run/tutorials/basic/scalarTransportFoam/pitzDaily
nProcs
1
sigFpe
Enabling floating point exception trapping (FOAM_SIGFPE).
fileModificationChecking Monitoring run-time modified files using timeStampMaster
allowSystemOperations
Disallowing user-supplied system call operations

// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
Create time

Create mesh
for time 0

/calculating face flux field phi

No finite volume options present

SIMPLE
no convergence criteria foundCalculations will run for 0.0002 steps.

Calculating scalar transport

Courant Number mean
0.447747 max1.94843
Time
0.0001

DILUPBiCG
:  Solving for TInitial residual 1, Final residual 6.02093e-07No Iterations 12
Time
0.0002

DILUPBiCG
:  Solving for TInitial residual 0.152376, Final residual 6.30498e-07No Iterations 11
End

[Inferior 1 (process 3591exited normally]
(
gdbpfvmatrixfull this x.dat
(gdbshell cat x.dat
cat
x.datNo such file or directory
(gdbquit
lydia
@lydia:~/OpenFOAM/run/tutorials/basic/scalarTransportFoam/pitzDailyl
0
/  0.0001/  0.0002/  constant/  fileAux  fileAux.gdbof  system

the fileAux.gdbof looks like that:
HTML Code:

No symbol "this" in current context.
As before - any hints about wrong usage are very appreciated!
Regards,
Lydia

 nisi May 20, 2014 13:49

Lydia,

In this case the problem is that the execution is not stopping in the specified breakpoint (it is not a gdbOF problem). If you check, you can see "End [Inferior 1 (process 3591) exited normally]", this says that the execution has finished with the last time-step, then the pointer "this" is out of scope (in fact there is nothing on memory at this point)

The main idea is to put a breakpoint exactly before the call
Code:

solverPerf = lduMatrix::solver::New
(
psi.name() + pTraits<Type>::componentNames[cmpt],
*this,
bouCoeffsCmpt,
intCoeffsCmpt,
interfaces,
solverControls
)->solve(psiCmpt, sourceCmpt, cmpt);

so, check in your own OF 2.2.2 version where this line is.
Then you can export the assembled matrix using our gdbOF commands.

 Lydia May 21, 2014 10:12

Thank you for the fast answer, Juan!

I'm starting to understand...

Unfortunately I didn't manage to solve the complete problem. As proposed I looked up where the breakpoint has to be set in OpenFoam_2.2.2. Then I started it and tried using pfvmatrixfull for writing the matrix.

I had to kill gdb since it used 100% of the CPU for more than 2 hours...so I still don't have the written matrix.

I used a 4 x 1 x 1 cell mesh and the scalarTransportSolver:

Maybe someone has time to have a look at the following protocol:
...Where is my mistake?:confused:
PHP Code:

gdb scalarTransportFoam
GNU gdb
(GDBSUSE (7.0-0.4.16)
This is free softwareyou are free to change and redistribute it.
There is NO WARRANTYto the extent permitted by law.  Type "show copying"
and "show warranty" for details.
This GDB was configured as "x86_64-suse-linux".
For
<
http://www.gnu.org/software/gdb/bugs/>...
(
gdbstart
Temporary breakpoint 1 at 0x43c9ed
file /home/schulzel/OpenFOAM/OpenFOAM-2.2.2_debug/src/OpenFOAM/lnInclude/setRootCase.Hline 5.
Starting program
: /panfs/home/schulzel/OpenFOAM/OpenFOAM-2.2.2_debug/platforms/linux64GccDPDebug/bin/scalarTransportFoam

Temporary breakpoint 1main (argc=1argv=0x7fffffffc728at /home/schulzel/OpenFOAM/OpenFOAM-2.2.2_debug/src/OpenFOAM/lnInclude/setRootCase.H:5
5           Foam
::argList args(argcargv);
(
gdbb fvScalarMatrix.C:158
Breakpoint 2 at 0x7fffecd961a2
file fvMatrices/fvScalarMatrix/fvScalarMatrix.Cline 158.
(gdb) continue
Continuing.
/*---------------------------------------------------------------------------*\
| =========                 |                                                 |
| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
|  \\    /   O peration     | Version:  2.2.2_debug                           |
|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
|    \\/     M anipulation  |                                                 |
\*---------------------------------------------------------------------------*/
Build  2.2.2-9739c53ec43f
Exec
: /panfs/home/schulzel/OpenFOAM/OpenFOAM-2.2.2_debug/platforms/linux64GccDPDebug/bin/scalarTransportFoam
Date
May 21 2014
Time
12:55:52
Host
"service0"
PID    2631
nProcs
1
sigFpe
Enabling floating point exception trapping (FOAM_SIGFPE).
fileModificationChecking Monitoring run-time modified files using timeStampMaster
allowSystemOperations
Disallowing user-supplied system call operations

// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
Create time

Create mesh
for time 0

/calculating face flux field phi

No finite volume options present

SIMPLE
no convergence criteria foundCalculations will run for 0.02 steps.

Calculating scalar transport

Courant Number mean
0.01 max0.01
Time
0.01

Breakpoint 2
Foam::fvMatrix<double>::solveSegregated (this=0x8366f0solverControls=...) at fvMatrices/fvScalarMatrix/fvScalarMatrix.C:161
warning
Source file is more recent than executable.
161         (
(
gdbl
156         scalarField totalSource
(source_);
158
159
// Solver call
160         solverPerformance solverPerf lduMatrix::solver::New
161         (
162             psi.name(),
163             *this,
164             boundaryCoeffs_,
165             internalCoeffs_,
(
gdbpfvmatrixsparse this e.dat
Getötet

Any help is much appreciated!!
Regards,
Lydia

 Lydia May 25, 2014 08:05

I/O error

Hello,

I managed to solve my problem now! Thank you for your help.
I tried it with an other solver and an other tutorial example(interFoam damBreak) and there it worked. I'm not yet sure, why it didn't work with the pitzDaily scalarTransportFoam example, but was probably some silly mistake.

Best regards,
Lydia

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