[davide.c]

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:

if
    (
        !mesh.thisDb().foundObject<GAMGAgglomeration>
        (
            GAMGAgglomeration::typeName
        )
    )

or

const petscControls& controls = petscControls::New
    (
        matrix_.mesh().thisDb().time()
    );

The issue is that lduPrimitiveMesh inherits from lduMesh the member function:

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:

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

[mAlletto]

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]

@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/

[olesen]

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]

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.

[dlahaye]

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.

[olesen]

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.

[dlahaye]

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.

[mAlletto]

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

[dlahaye]

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/

[dlahaye]

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.

[mAlletto]

Hello, what do you mean exactly by overlay computations and computations?

[Santiago]

Originally Posted by dlahaye

@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.

[dlahaye]

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]

@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.

[dlahaye]

Overlapping communication and computation is well explained at
https://netlib.org/linalg/html_templ...igrearrangedcg.

Is this implemented in OpenFoam?

[mAlletto]

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

[dlahaye]

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.