CFD Online Logo CFD Online URL
www.cfd-online.com
[Sponsors]
Home > Forums > Software User Forums > OpenFOAM > OpenFOAM Programming & Development

lduMesh, lduAddressing and thisDb()

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

Like Tree5Likes
  • 1 Post By davide.c
  • 1 Post By mAlletto
  • 1 Post By olesen
  • 1 Post By dlahaye
  • 1 Post By dlahaye

Reply
 
LinkBack Thread Tools Search this Thread Display Modes
Old   April 27, 2022, 11:50
Default lduMesh, lduAddressing and thisDb()
  #1
New Member
 
Davide Cortellessa
Join Date: Jul 2019
Posts: 4
Rep Power: 6
davide.c is on a distinguished road
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 external-solver 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
        )
    )
or

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_;
}
Therefore, my first conclusion (which I hope to be wrong) is that if a solver require access to the member function thisDb() itís not possible to use the class lduPrimitiveMesh to define the addressing.

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_;
}
As such, my second conclusion is that itís not possible to have a fvMesh object with an lduAddressing different from the predefinite OpenFOAM one.

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
dlahaye likes this.
davide.c is offline   Reply With Quote

Old   April 29, 2022, 04:04
Default
  #2
Senior Member
 
Michael Alletto
Join Date: Jun 2018
Location: Bremen
Posts: 615
Rep Power: 15
mAlletto will become famous soon enough
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
dlahaye likes this.
mAlletto is offline   Reply With Quote

Old   April 29, 2022, 05:07
Default
  #3
Senior Member
 
Domenico Lahaye
Join Date: Dec 2013
Posts: 751
Blog Entries: 1
Rep Power: 17
dlahaye is on a distinguished road
@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 re-inventing 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/
dlahaye is offline   Reply With Quote

Old   April 30, 2022, 11:06
Default
  #4
Senior Member
 
Mark Olesen
Join Date: Mar 2009
Location: https://olesenm.github.io/
Posts: 1,708
Rep Power: 40
olesen has a spectacular aura aboutolesen has a spectacular aura about
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.
olesen is offline   Reply With Quote

Old   April 30, 2022, 11:09
Default
  #5
Senior Member
 
Mark Olesen
Join Date: Mar 2009
Location: https://olesenm.github.io/
Posts: 1,708
Rep Power: 40
olesen has a spectacular aura aboutolesen has a spectacular aura about
The thisDb() is just a reference to the associated objectRegistry - ie, a dumping ground (database) of various regIOobjects and a top-level Time - nothing in there particular to an faMesh or fvMesh.
olesen is offline   Reply With Quote

Old   May 2, 2022, 09:26
Default
  #6
Senior Member
 
Domenico Lahaye
Join Date: Dec 2013
Posts: 751
Blog Entries: 1
Rep Power: 17
dlahaye is on a distinguished road
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 fill-in) 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.
dlahaye is offline   Reply With Quote

Old   May 5, 2022, 04:01
Default
  #7
Senior Member
 
Mark Olesen
Join Date: Mar 2009
Location: https://olesenm.github.io/
Posts: 1,708
Rep Power: 40
olesen has a spectacular aura aboutolesen has a spectacular aura about
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 off-processor components. The local matrices correspond to a block of diagonals in the global view and the inter-processor connections are everything off-diagonal (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 off-processor components in the PETSc language) are handled by the ldu-interfaces, which correspond to processor-processor 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 ldu-interfaces to setup the off-processor 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 low-level and invasive change and you have to pre-calculate 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.
Attached Images
File Type: jpg petsc-matrix-layout.jpg (90.9 KB, 33 views)
dlahaye likes this.
olesen is offline   Reply With Quote

Old   May 9, 2022, 15:58
Default
  #8
Senior Member
 
Domenico Lahaye
Join Date: Dec 2013
Posts: 751
Blog Entries: 1
Rep Power: 17
dlahaye is on a distinguished road
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 pressure-solve in incompressible case, no AMD reordering and ILU with drop-tolerance for pressure-solve in transonic case). The profiling of the OpenFoam linear solvers remains difficult. This is currently hurting my simulations a lot.

PETSc provides support for finite-volume 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 drop-tolerance 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.
dlahaye is offline   Reply With Quote

Old   May 23, 2022, 02:48
Default
  #9
Senior Member
 
Michael Alletto
Join Date: Jun 2018
Location: Bremen
Posts: 615
Rep Power: 15
mAlletto will become famous soon enough
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
mAlletto is offline   Reply With Quote

Old   May 23, 2022, 05:48
Default
  #10
Senior Member
 
Domenico Lahaye
Join Date: Dec 2013
Posts: 751
Blog Entries: 1
Rep Power: 17
dlahaye is on a distinguished road
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://www-users.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 non-overlapping 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.prace-ri.eu : this reference provides examples and listings that illustrate the above discussion;

[2] https://prace-ri.eu/wp-content/uploa...-Framework.pdf

[3]: https://www.linkedin.com/pulse/openf...menico-lahaye/
mAlletto likes this.
dlahaye is offline   Reply With Quote

Old   August 5, 2022, 09:15
Default
  #11
Senior Member
 
Domenico Lahaye
Join Date: Dec 2013
Posts: 751
Blog Entries: 1
Rep Power: 17
dlahaye is on a distinguished road
Question:

what measures does OpenFOAM have in place to overlay computations and computations in vector-vector and matrix-vector operations in its precondition Krylov subspace solvers?

Thanks.
dlahaye is offline   Reply With Quote

Old   August 6, 2022, 03:59
Default
  #12
Senior Member
 
Michael Alletto
Join Date: Jun 2018
Location: Bremen
Posts: 615
Rep Power: 15
mAlletto will become famous soon enough
Hello, what do you mean exactly by overlay computations and computations?
mAlletto is offline   Reply With Quote

Old   August 7, 2022, 04:53
Default
  #13
Senior Member
 
Santiago Lopez Castano
Join Date: Nov 2012
Posts: 354
Rep Power: 15
Santiago is on a distinguished road
Quote:
Originally Posted by dlahaye View Post
@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 re-inventing 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/
Why the trouble of going through all these retro-fittings and code-obliteration/obsfucation trying to adapt PETSc into FOAM? Wouldn't it be just easier to implement your solver directly within the TAO framework? I mean, they offer everything, and more, of what FOAM offers. If you are fluent in the TAO library, I really don't see the reason why you're going through all this effort. I, perhaps, am too stupid to see it though.

This feels to me like a convoluted form of spaghetti programming, but not in FORTRAN anymore.
Santiago is offline   Reply With Quote

Old   August 8, 2022, 04:37
Default
  #14
Senior Member
 
Domenico Lahaye
Join Date: Dec 2013
Posts: 751
Blog Entries: 1
Rep Power: 17
dlahaye is on a distinguished road
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
dlahaye is offline   Reply With Quote

Old   August 8, 2022, 04:51
Default
  #15
Senior Member
 
Domenico Lahaye
Join Date: Dec 2013
Posts: 751
Blog Entries: 1
Rep Power: 17
dlahaye is on a distinguished road
@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.
Santiago likes this.
dlahaye is offline   Reply With Quote

Old   August 12, 2022, 15:14
Default
  #16
Senior Member
 
Domenico Lahaye
Join Date: Dec 2013
Posts: 751
Blog Entries: 1
Rep Power: 17
dlahaye is on a distinguished road
Overlapping communication and computation is well explained at
https://netlib.org/linalg/html_templ...igrearrangedcg.

Is this implemented in OpenFoam?
dlahaye is offline   Reply With Quote

Old   August 13, 2022, 04:40
Default
  #17
Senior Member
 
Michael Alletto
Join Date: Jun 2018
Location: Bremen
Posts: 615
Rep Power: 15
mAlletto will become famous soon enough
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
mAlletto is offline   Reply With Quote

Old   August 14, 2022, 13:27
Default
  #18
Senior Member
 
Domenico Lahaye
Join Date: Dec 2013
Posts: 751
Blog Entries: 1
Rep Power: 17
dlahaye is on a distinguished road
Dear Michael,

Thank you for sharing this valuable reference. It is references like these that increase my wish to understand the underlying technology.

Kind wishes, Domenico.
dlahaye is offline   Reply With Quote

Reply

Tags
lduaddressing, ldumatrix, ldumesh, lduprimitivemesh, thisdb

Thread Tools Search this Thread
Search this Thread:

Advanced Search
Display Modes

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 Off
Trackbacks are Off
Pingbacks are On
Refbacks are On



All times are GMT -4. The time now is 12:27.