CFD Online Discussion Forums (https://www.cfd-online.com/Forums/)
-   OpenFOAM Programming & Development (https://www.cfd-online.com/Forums/openfoam-programming-development/)
-   -   Matrix manipulation (https://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/bashrclydia@lydia:~/OpenFOAM\$ cd run/tutorials/basic/scalarTransportFoam/pitzDaily/lydia@lydia:~/OpenFOAM/run/tutorials/basic/scalarTransportFoam/pitzDaily\$ gdb scalarTransportFoamGNU gdb (Ubuntu 7.7-0ubuntu3) 7.7Copyright (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) startTemporary 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:55        Foam::argList args(argc, argv);(gdb) startThe program being debugged has been started already.Start it from the beginning? (y or n) yTemporary 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:55        Foam::argList args(argc, argv);(gdb) b fvMatrixSolve.C:140Breakpoint 3 at 0x44f6f1: fvMatrixSolve.C:140. (3 locations)(gdb) continueContinuing./*---------------------------------------------------------------------------*\| =========                 |                                                 || \\      /  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-9739c53ec43fExec   : /home/lydia/OpenFOAM/OpenFOAM-2.2.2/platforms/linux64GccDPDebug/bin/scalarTransportFoamDate   : May 20 2014Time   : 18:53:06Host   : "lydia-Latitude-D630"PID    : 3591Case   : /home/lydia/OpenFOAM/run/tutorials/basic/scalarTransportFoam/pitzDailynProcs : 1sigFpe : Enabling floating point exception trapping (FOAM_SIGFPE).fileModificationChecking : Monitoring run-time modified files using timeStampMasterallowSystemOperations : Disallowing user-supplied system call operations// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //Create timeCreate mesh for time = 0Reading field TReading field UReading transportPropertiesReading diffusivity DTReading/calculating face flux field phiNo finite volume options presentSIMPLE: no convergence criteria found. Calculations will run for 0.0002 steps.Calculating scalar transportCourant Number mean: 0.447747 max: 1.94843Time = 0.0001DILUPBiCG:  Solving for T, Initial residual = 1, Final residual = 6.02093e-07, No Iterations 12Time = 0.0002DILUPBiCG:  Solving for T, Initial residual = 0.152376, Final residual = 6.30498e-07, No Iterations 11End[Inferior 1 (process 3591) exited normally](gdb) pfvmatrixfull this x.dat(gdb) shell cat x.datcat: x.dat: No such file or directory(gdb) quitlydia@lydia:~/OpenFOAM/run/tutorials/basic/scalarTransportFoam/pitzDaily\$ l0/  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 (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) startTemporary 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:55           Foam::argList args(argc, argv);(gdb) b fvScalarMatrix.C:158Breakpoint 2 at 0x7fffecd961a2: file fvMatrices/fvScalarMatrix/fvScalarMatrix.C, line 158.(gdb) continueContinuing./*---------------------------------------------------------------------------*\| =========                 |                                                 || \\      /  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-9739c53ec43fExec   : /panfs/home/schulzel/OpenFOAM/OpenFOAM-2.2.2_debug/platforms/linux64GccDPDebug/bin/scalarTransportFoamDate   : May 21 2014Time   : 12:55:52Host   : "service0"PID    : 2631Case   : /panfs/w3/schulzel/NUMAP/06_testcase-advectionFoam/advectionFoam-gdbnProcs : 1sigFpe : Enabling floating point exception trapping (FOAM_SIGFPE).fileModificationChecking : Monitoring run-time modified files using timeStampMasterallowSystemOperations : Disallowing user-supplied system call operations// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //Create timeCreate mesh for time = 0Reading field TReading field UReading transportPropertiesReading diffusivity DTReading/calculating face flux field phiNo finite volume options presentSIMPLE: no convergence criteria found. Calculations will run for 0.02 steps.Calculating scalar transportCourant Number mean: 0.01 max: 0.01Time = 0.01Breakpoint 2, Foam::fvMatrix<double>::solveSegregated (this=0x8366f0, solverControls=...) at fvMatrices/fvScalarMatrix/fvScalarMatrix.C:161warning: Source file is more recent than executable.161         ((gdb) l156         scalarField totalSource(source_);157         addBoundarySource(totalSource, false);158159         // Solver call160         solverPerformance solverPerf = lduMatrix::solver::New161         (162             psi.name(),163             *this,164             boundaryCoeffs_,165             internalCoeffs_,(gdb) pfvmatrixsparse this e.datGetö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 08:02.