[Linse]

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!

[andyru]

Hi,

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

volScalarField rUA = 1.0/UEqn.A();

Here you can acces them an manipulate them.

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)

What is your physical application?


Greetings

Andy

[santiagomarquezd]

Originally Posted by Linse

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]

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

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]

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]

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

Regards

[j-avdeev]

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"

$ (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

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 file:

$n (next) will step to the next line in the current file, but will not go into functions and
included files.
$s (step) will go to the next line, also in functions and included files.

- 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]

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]

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

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]

Well,

Originally Posted by j-avdeev

it would be cool for me to know this moment, where equation like

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)

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
$ pfvmatrixfull this AdvDiff.dat
$ shell cat AdvDiff.dat

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]

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?
Thanks a lot in advance.

Regards,
Lydia

[nisi]

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]

Dear Juan,
thanks for your reply. I downloaded the latest version of gdbof - unfortunately I get the same result as described above."

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:~/OpenFOAM$ cd run/tutorials/basic/scalarTransportFoam/pitzDaily/
lydia@lydia:~/OpenFOAM/run/tutorials/basic/scalarTransportFoam/pitzDaily$ gdb scalarTransportFoam
GNU gdb (Ubuntu 7.7-0ubuntu3) 7.7
Copyright (C) 2014 Free Software Foundation, Inc.
License GPLv3+: GNU GPL version 3 or later <http://gnu.org/licenses/gpl.html>
This is free software: you are free to change and redistribute it.
There is NO WARRANTY, to 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 bug reporting instructions, please see:
<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 help, type "help".
Type "apropos word" to search for commands related to "word"...
Reading symbols from scalarTransportFoam...done.
(gdb) start
Temporary breakpoint 1 at 0x4384a6: file /home/lydia/OpenFOAM/OpenFOAM-2.2.2/src/OpenFOAM/lnInclude/setRootCase.H, line 5.
Starting program: /home/lydia/OpenFOAM/OpenFOAM-2.2.2/platforms/linux64GccDPDebug/bin/scalarTransportFoam 

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

Temporary breakpoint 2, main (argc=1, argv=0x7fffffffcfc8) at /home/lydia/OpenFOAM/OpenFOAM-2.2.2/src/OpenFOAM/lnInclude/setRootCase.H:5
5        Foam::argList args(argc, argv);
(gdb) b 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

Reading field T

Reading field U

Reading transportProperties

Reading diffusivity DT

Reading/calculating face flux field phi

No finite volume options present


SIMPLE: no convergence criteria found. Calculations will run for 0.0002 steps.


Calculating scalar transport

Courant Number mean: 0.447747 max: 1.94843
Time = 0.0001

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

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

[Inferior 1 (process 3591) exited normally]
(gdb) pfvmatrixfull this x.dat
(gdb) shell cat x.dat
cat: x.dat: No such file or directory
(gdb) quit
lydia@lydia:~/OpenFOAM/run/tutorials/basic/scalarTransportFoam/pitzDaily$ l
0/  0.0001/  0.0002/  constant/  fileAux  fileAux.gdbof  system/ 


the fileAux.gdbof looks like that:

No symbol "this" in current context.

As before - any hints about wrong usage are very appreciated!
Thank you in advance.
Regards,
Lydia

[nisi]

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

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]

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?

PHP Code:

gdb scalarTransportFoam 
GNU gdb (GDB) SUSE (7.0-0.4.16)
Copyright (C) 2009 Free Software Foundation, Inc.
License GPLv3+: GNU GPL version 3 or later <http://gnu.org/licenses/gpl.html>
This is free software: you are free to change and redistribute it.
There is NO WARRANTY, to the extent permitted by law.  Type "show copying"
and "show warranty" for details.
This GDB was configured as "x86_64-suse-linux".
For bug reporting instructions, please see:
<http://www.gnu.org/software/gdb/bugs/>...
Reading symbols from /panfs/home/schulzel/OpenFOAM/OpenFOAM-2.2.2_debug/platforms/linux64GccDPDebug/bin/scalarTransportFoam...done.
(gdb) start
Temporary breakpoint 1 at 0x43c9ed: file /home/schulzel/OpenFOAM/OpenFOAM-2.2.2_debug/src/OpenFOAM/lnInclude/setRootCase.H, line 5.
Starting program: /panfs/home/schulzel/OpenFOAM/OpenFOAM-2.2.2_debug/platforms/linux64GccDPDebug/bin/scalarTransportFoam 
[Thread debugging using libthread_db enabled]

Temporary breakpoint 1, main (argc=1, argv=0x7fffffffc728) at /home/schulzel/OpenFOAM/OpenFOAM-2.2.2_debug/src/OpenFOAM/lnInclude/setRootCase.H:5
5           Foam::argList args(argc, argv);
(gdb) b fvScalarMatrix.C:158
Breakpoint 2 at 0x7fffecd961a2: file fvMatrices/fvScalarMatrix/fvScalarMatrix.C, line 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
Case   : /panfs/w3/schulzel/NUMAP/06_testcase-advectionFoam/advectionFoam-gdb
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

Reading field T

Reading field U

Reading transportProperties

Reading diffusivity DT

Reading/calculating face flux field phi

No finite volume options present


SIMPLE: no convergence criteria found. Calculations will run for 0.02 steps.


Calculating scalar transport

Courant Number mean: 0.01 max: 0.01
Time = 0.01


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


Any help is much appreciated!!
Regards,
Lydia

[Lydia]

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