
[Sponsors] 
April 27, 2022, 11:50 
lduMesh, lduAddressing and thisDb()

#1 
New Member
Davide Cortellessa
Join Date: Jul 2019
Posts: 4
Rep Power: 6 
Dear foamers,
Within the context of my PhD thesis, I need to have a matrix with a sparsity pattern different from the OpenFOAM usual one. Therefore, I’m using the class lduMatrix, which takes the addressing from an object of the class lduPrimitiveMesh, which can be constructed providing size, lower and upper addressing. The problem comes when I try to use GAMG solver or preconditioner, or a solver from the petsc4Foam interface from the externalsolver module. This is because in these codes, at some point, they use the function thisDb(), for example: Code:
if ( !mesh.thisDb().foundObject<GAMGAgglomeration> ( GAMGAgglomeration::typeName ) ) Code:
const petscControls& controls = petscControls::New ( matrix_.mesh().thisDb().time() ); The issue is that lduPrimitiveMesh inherits from lduMesh the member function: Code:
const Foam::objectRegistry& Foam::lduMesh::thisDb() const { NotImplemented; const objectRegistry* orPtr_ = nullptr; return *orPtr_; } As a consequence I started looking at the class fvMesh. Unfortunately, if I’m not mistaken, it’s not possible to change the lduAddrssing of an existing fvMesh object, nor it’s possible to provide one at its creation, as the access function is also responsible for constructing the addressing the first time the access is attempted: Code:
const Foam::lduAddressing& Foam::fvMesh::lduAddr() const { if (!lduPtr_) { DebugInFunction << "Calculating fvMeshLduAddressing from nFaces:" << nFaces() << endl; lduPtr_ = new fvMeshLduAddressing(*this); return *lduPtr_; } return *lduPtr_; } Of course, I realize that it’s impossible that I’m the first one in the world who attempts to solve in OpenFOAM a linear system with a different sparsity pattern, but, as of now, I’m not able to see how to do it. Thanks in advance for any inputs! Davide 

April 29, 2022, 04:04 

#2 
Senior Member
Michael Alletto
Join Date: Jun 2018
Location: Bremen
Posts: 615
Rep Power: 15 
Hello Davide,
actually changing the addressing is possible and it is done in the overset method in the com version of OpenFOAM. I started with a description of overPimpleDymFOAM here: https://openfoamwiki.net/index.php/OverPimpleDyMFoam. In the file fvMeshPrimitiveLduAddressing.C and function addAddressing the lowerAddr() and upperAddr() arrays are extended to account for the additional coupling due to the overset interpolation. This function is called within the dynamicOversetFvMesh class which is a child of fvMesh. It provided also the functionality to access the database. So you can have a look at this file to inspire you how to extend your addressing. Best Michael 

April 29, 2022, 05:07 

#3 
Senior Member

@davide: sincere thanks for raising this important issue!
Two ideas: 1/ you mention petsc4Foam. Would it be feasible to take this a step further and employ the DMPLEX infrastructure that PETSc provides [1] to avoid reinventing the (admittedly sophisticated) wheel; 2/ does it make sense to leverage the linear algebra to Julia using FvCFD.jl (or similar) hoping to accelerate overall development? Domenico. [1]: https://petsc.org/release/docs/manual/dmplex/ 

April 30, 2022, 11:06 

#4 
Senior Member
Mark Olesen
Join Date: Mar 2009
Location: https://olesenm.github.io/
Posts: 1,708
Rep Power: 40 
The newer PETSc versions have a fast CCO to CSR that is used by the petsc4Foam interface, so probably still safer to go with that instead of the dmplex. However, if you find that the dmplex approach wraps better, I'd be happy to take a look.


April 30, 2022, 11:09 

#5 
Senior Member
Mark Olesen
Join Date: Mar 2009
Location: https://olesenm.github.io/
Posts: 1,708
Rep Power: 40 
The thisDb() is just a reference to the associated objectRegistry  ie, a dumping ground (database) of various regIOobjects and a toplevel Time  nothing in there particular to an faMesh or fvMesh.


May 2, 2022, 09:26 

#6 
Senior Member

Dear Mark,
Thank you for your input! My objective is to take advantage of the reordering (approximate minimal degree (AMD)) schemes and advanced preconditioners (incomplete LU factorization with fillin) in parallel that PETSc provides. I see two ways to approach this. 1/ The first approach is to intervene on the linear solver level and replace OpenFoam native solvers by PETSc equivalent. This is straightforward for a sequential implementation. I have no clear understanding of what is required for a parallel implementation. Does petsc4Foam leverage the parallel linear algebra from OpenFoam to PETSc? Does it include the mesh decomposition? 2/ The second approach is to intervene at the discretization level and replace finite volume discretization in OpenFoam by PETSc equivalent using DMPLEX and PETScFV. This is more work, but components already exist. The advantage of the second approach is to take full advantage of PETSc. I suggest to try this approach for laplacianFoam and scalarTransportFoam first. What is your point on the above? Thanks, Domenico. 

May 5, 2022, 04:01 

#7 
Senior Member
Mark Olesen
Join Date: Mar 2009
Location: https://olesenm.github.io/
Posts: 1,708
Rep Power: 40 
Hi Domenico,
Before launching into coding anything, you need to take stock of the components and plan where this going to go. Lets work backwards from the goal (new linear solver implementation in native petsc/mkl whatever). The matrix that you want to solve will be spread across multiple processor ranks. In the PETSc world they use a global view of that matrix, but still distinguish between the local matrix and the offprocessor components. The local matrices correspond to a block of diagonals in the global view and the interprocessor connections are everything offdiagonal (eg, attached slide is from an overview by Simone Bna a few years back). You will have the same thing in OpenFOAM, but using a slightly different representation. As you noted, OpenFOAM stores as an LDU addressing (ie, a special sorted CCO addressing format) for the local matrices. A global matrix view is not used (or needed). The connections between local matrices (ie, the offprocessor components in the PETSc language) are handled by the lduinterfaces, which correspond to processorprocessor boundary patches (or any coupled patch). The petsc4Foam layer is what is responsible for transcribing the ldu representation into the PETSc global view, including walking over the lduinterfaces to setup the offprocessor matrix components for PETSc. Sure we have a memory and cpu overhead of copying the components across, but we have a real gain since the LDU format is much easier to use than the same thing with CSR. If you want to avoid this transcription, you would need to assemble the matrix contributions directly into the PETSc matrix locations. This means replacing the assembly provided by fvMatrix with a different internal layout (CSR for example). This is a fairly lowlevel and invasive change and you have to precalculate the CSR layout before you can start inserting. In this case you might instead decide to collect where the contributions go before actually positioning them into the CSR (since we would want to avoid unnecessary shuffling of the CSR contents during assembly). However this sounds very much like fvMatrix + petsc4Foam. If we go up one step even further and use PETSc for the discretization schemes etc (relaxation? constraints? source terms?)  this is not only radical, but there isn't much left that is OpenFOAM beyond the geometry representation at that stage. 

May 9, 2022, 15:58 

#8 
Senior Member

Dear Marc,
Thank you for sharing your insight. I agree that OpenFoam and PETSc implement the same parallel matrix and vector layout. I understand that the OpenFoam LDU format blends well with the finite volume discretization. The fact is that the linear solvers in OpenFoam are naive (no genuine algebraic multigrid for pressuresolve in incompressible case, no AMD reordering and ILU with droptolerance for pressuresolve in transonic case). The profiling of the OpenFoam linear solvers remains difficult. This is currently hurting my simulations a lot. PETSc provides support for finitevolume discretization. Support for boundary layers, thermodynamics, combustion and radiative heat support, however, is missing in PETSc. Some blending of OpenFoam and PETSc is thus required. I agree that I need to prepare my homework more carefully. Your input provides a strong basis for further developments. petsc4Foam should allow to 1/ in the incompressible case compare GAMG in OpenFoam with genuine AMG (no initial G) in PETSc for various mesh sizes and number of cores; 2/ in the transonic case compare solver in OpenFoam with ILU with droptolerance and AMD reordering in PETSc for various mesh sizes and number of cores; Once the argument is sufficiently well developed, one could discuss a tighter integration of OpenFoam and PETSc in a next step. Does the above make sense? Should I provide more details? Can this be discussed in a wider circle at the forthcoming OpenFoam workshop? Cheers, Domenico. 

May 23, 2022, 02:48 

#9 
Senior Member
Michael Alletto
Join Date: Jun 2018
Location: Bremen
Posts: 615
Rep Power: 15 
Hello all,
do you know if there are some text books or other documents where algorithms to to solve in parallel linear systems of equations are described. Are there documents where the parallel algorithms implemented in OF are described? Best Michael 

May 23, 2022, 05:48 

#10 
Senior Member

Yes and possibly no.
Yes, the literature on parallel linear solvers is vast. Some pointers are http://ddm.org and the book by Youssef Saad https://wwwusers.cse.umn.edu/~saad/books.html (among many others). It really depends on what you are looking for. I am happy to provide more detailed info in case requested. Yes, OpenFoam implements a classical nonoverlapping domain decomposition algorithm. Documents on the implementation in OpenFoam include [1] and [2]. Possibly no, I have not able to work through the class hierarchy of primitive mesh, ldumesh and fvMatrix to grasp a full understanding of the implementation of the subdomain coupling through interfaces in OpenFoam. A first attempt to describe the implementation is [3]. A better understanding would be valuable to obtain. [1] M. Culpo, “Current Bottlenecks in the Scalability of OpenFOAM on Massively Parallel Computers”, CINECA technical report available at www.praceri.eu : this reference provides examples and listings that illustrate the above discussion; [2] https://praceri.eu/wpcontent/uploa...Framework.pdf [3]: https://www.linkedin.com/pulse/openf...menicolahaye/ 

August 5, 2022, 09:15 

#11 
Senior Member

Question:
what measures does OpenFOAM have in place to overlay computations and computations in vectorvector and matrixvector operations in its precondition Krylov subspace solvers? Thanks. 

August 6, 2022, 03:59 

#12 
Senior Member
Michael Alletto
Join Date: Jun 2018
Location: Bremen
Posts: 615
Rep Power: 15 
Hello, what do you mean exactly by overlay computations and computations?


August 7, 2022, 04:53 

#13  
Senior Member
Santiago Lopez Castano
Join Date: Nov 2012
Posts: 354
Rep Power: 15 
Quote:
This feels to me like a convoluted form of spaghetti programming, but not in FORTRAN anymore. 

August 8, 2022, 04:37 

#14 
Senior Member

I mean to allow the following two processes to occur simultaneously:
a/ local computations on a given processor; b/ data transfer to/from the given processor; See e.g. Section 2.2 of https://dl.acm.org/doi/pdf/10.1145/3...90a_6cZiK0284w 

August 8, 2022, 04:51 

#15 
Senior Member

@Santiago: thanks you for raising your valid concerns.
My interest is exploring ways to accelerate computations in OpenFOAM. The physics of turbulent reactive flow e.g. is far less developed in PETSc/TAO. Extending the physics of PETSc/TAO is indeed an alternative route to explore. 

August 12, 2022, 15:14 

#16 
Senior Member

Overlapping communication and computation is well explained at
https://netlib.org/linalg/html_templ...igrearrangedcg. Is this implemented in OpenFoam? 

August 13, 2022, 04:40 

#17 
Senior Member
Michael Alletto
Join Date: Jun 2018
Location: Bremen
Posts: 615
Rep Power: 15 
I cannot remember anything similar in OpenFOAM. However I saw a presentation https://www.linkedin.com/posts/charl...um=android_app
Where they could achieve super linear scaling on CPU with a cache of 128 M. Linear scaling was achieved until 50K cells per core. Afterwards communication started to impact the performance. Best Michael 

Tags 
lduaddressing, ldumatrix, ldumesh, lduprimitivemesh, thisdb 
Thread Tools  Search this Thread 
Display Modes  

