# Matrix manipulation

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

 February 13, 2012, 11:23 Matrix manipulation #1 Senior Member   Bernhard Linseisen Join Date: May 2010 Location: Magdeburg/Geneva Posts: 176 Blog Entries: 1 Rep Power: 7 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!

 February 13, 2012, 13:42 #2 New Member   Andreas Ruopp Join Date: Aug 2009 Location: Stuttgart / Germany Posts: 29 Rep Power: 7 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) (Manipulating source fvc::grad(vsf1+vsf2)) What is your physical application? Greetings Andy

February 21, 2012, 22:28
#3
Senior Member

Santiago Marquez Damian
Join Date: Aug 2009
Location: Santa Fe, Santa Fe, Argentina
Posts: 418
Rep Power: 14
Quote:
 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.
__________________
Santiago MÁRQUEZ DAMIÁN, Ph.D.
Post-doctoral Fellow
Research Center for Computational Mechanics (CIMEC) - CONICET/FICH-UNL
T.E.: 54-342-4511594 Ext. 1005
Güemes 3450 - (3000) Santa Fe
Santa Fe - Argentina
http://www.cimec.org.ar

April 4, 2012, 05:58
#4
Member

Avdeev Evgeniy
Join Date: Jan 2011
Location: Togliatty, Russia
Posts: 52
Blog Entries: 1
Rep Power: 12
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?

 April 4, 2012, 07:08 #5 Senior Member   Elvis Join Date: Mar 2009 Location: Sindelfingen, Germany Posts: 577 Blog Entries: 5 Rep Power: 13 I guess they want to do something similar that Santiago Marquez Damian described in his post Can someone tell me more about vulashaka

 April 4, 2012, 08:52 #6 Senior Member     Santiago Marquez Damian Join Date: Aug 2009 Location: Santa Fe, Santa Fe, Argentina Posts: 418 Rep Power: 14 Maybe you need http://openfoamwiki.net/index.php/Contrib_gdbOF, it has an specific tool to do that. Regards __________________ Santiago MÁRQUEZ DAMIÁN, Ph.D. Post-doctoral Fellow Research Center for Computational Mechanics (CIMEC) - CONICET/FICH-UNL T.E.: 54-342-4511594 Ext. 1005 Güemes 3450 - (3000) Santa Fe Santa Fe - Argentina http://www.cimec.org.ar

 April 23, 2012, 17:27 #7 Member   Avdeev Evgeniy Join Date: Jan 2011 Location: Togliatty, Russia Posts: 52 Blog Entries: 1 Rep Power: 12 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.

 April 23, 2012, 17:33 #8 Senior Member     Santiago Marquez Damian Join Date: Aug 2009 Location: Santa Fe, Santa Fe, Argentina Posts: 418 Rep Power: 14 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 __________________ Santiago MÁRQUEZ DAMIÁN, Ph.D. Post-doctoral Fellow Research Center for Computational Mechanics (CIMEC) - CONICET/FICH-UNL T.E.: 54-342-4511594 Ext. 1005 Güemes 3450 - (3000) Santa Fe Santa Fe - Argentina http://www.cimec.org.ar

 April 24, 2012, 16:11 #9 Member   Avdeev Evgeniy Join Date: Jan 2011 Location: Togliatty, Russia Posts: 52 Blog Entries: 1 Rep Power: 12 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.

April 24, 2012, 16:28
#10
Senior Member

Santiago Marquez Damian
Join Date: Aug 2009
Location: Santa Fe, Santa Fe, Argentina
Posts: 418
Rep Power: 14
Well,

Quote:
 Originally Posted by j-avdeev 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.
__________________
Santiago MÁRQUEZ DAMIÁN, Ph.D.
Post-doctoral Fellow
Research Center for Computational Mechanics (CIMEC) - CONICET/FICH-UNL
T.E.: 54-342-4511594 Ext. 1005
Güemes 3450 - (3000) Santa Fe
Santa Fe - Argentina
http://www.cimec.org.ar

 May 2, 2014, 08:59 cannot write pvfmatrixfull with gdbof #11 New Member   Lydia Schulze Join Date: Jan 2012 Location: Karlsruhe, Germany Posts: 18 Rep Power: 5 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

 May 15, 2014, 16:35 #12 New Member   Juan Marcelo Gimenez Join Date: Dec 2009 Location: Santa Fe, Argentina Posts: 11 Rep Power: 7 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

 May 20, 2014, 13:05 problem not yet solved - cannot write pvfmatrixfull with gdbof_v1.03a #13 New Member   Lydia Schulze Join Date: Jan 2012 Location: Karlsruhe, Germany Posts: 18 Rep Power: 5 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/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 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:.Find the GDB manual and other documentation resources online at:.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! Thank you in advance. Regards, Lydia

 May 20, 2014, 13:49 #14 New Member   Juan Marcelo Gimenez Join Date: Dec 2009 Location: Santa Fe, Argentina Posts: 11 Rep Power: 7 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::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.

 May 21, 2014, 10:12 #15 New Member   Lydia Schulze Join Date: Jan 2012 Location: Karlsruhe, Germany Posts: 18 Rep Power: 5 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 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:...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::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

 May 25, 2014, 08:05 I/O error #16 New Member   Lydia Schulze Join Date: Jan 2012 Location: Karlsruhe, Germany Posts: 18 Rep Power: 5 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 Last edited by Lydia; May 25, 2014 at 11:31.

 Thread Tools Display Modes Linear Mode

 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 OffTrackbacks are On Pingbacks are On Refbacks are On Forum Rules

 Similar Threads Thread Thread Starter Forum Replies Last Post colopolo CFX 13 October 4, 2011 22:03 arjun shanker Main CFD Forum 0 November 12, 2010 23:12 lakeat OpenFOAM Running, Solving & CFD 42 August 26, 2009 21:47 Demidov Main CFD Forum 0 November 14, 2004 05:51 xueying Main CFD Forum 2 September 24, 2002 09:44

All times are GMT -4. The time now is 23:20.