Operators like div/grad

[HenningScheufler]

See proofofconcept branch

Setup the testcase to compare the implemented operators to openfoam

explicit operators

• [x] interpolation scheme
• [x] grad operator
  • [x] scalar
  • [x] include boundary conditions
  • [x] setup test cases
• [x] div operator
• [ ] ddt operator
• [ ] laplacian operator
• [ ] source terms

implicit operators

• [ ] source terms
• [ ] div operator
• [ ] ddt operator
• [ ] laplacian operator

[HenningScheufler]

The current implementation of the gradient scheme looks at followed:

void GaussGreenKernel::operator()(const GPUExecutor& exec)
{
    using executor = typename GPUExecutor::exec;
    const NeoFOAM::labelField& owner = mesh_.faceOwner();
    const NeoFOAM::labelField& neighbour = mesh_.faceNeighbour();
    const NeoFOAM::vectorField& Sf = mesh_.faceAreas();
    const NeoFOAM::scalarField& V = mesh_.cellVolumes();
    auto s_gradPhi = gradPhi_.field();
    auto s_phi = phi_.field();
    auto s_owner = owner.field();
    auto s_neighbour = neighbour.field();
    auto s_Sf = Sf.field();
    auto s_V = V.field();

    Kokkos::parallel_for(
        "gaussGreenGrad", Kokkos::RangePolicy<executor>(0, mesh_.nInternalFaces()), KOKKOS_LAMBDA(const int i) {
            int32_t own = s_owner[i];
            int32_t nei = s_neighbour[i];
            NeoFOAM::scalar phif = 0.5 * (s_phi[nei] + s_phi[own]);
            NeoFOAM::Vector value_own = s_Sf[i] * (phif / s_V[own]);
            NeoFOAM::Vector value_nei = s_Sf[i] * (phif / s_V[nei]);
            Kokkos::atomic_add(&s_gradPhi[own], value_own);
            Kokkos::atomic_sub(&s_gradPhi[nei], value_nei);
        }
    );
}

With this implementation, the mesh's job is to provide the connectivity information necessary to access neighboring cells. For simple computations like the gradient, this is not an issue, but more complicated methods require different connectivity data. Consequently, the mesh class would become very large, resulting in a massive interface. To avoid this, OpenFOAM is able to register stencils for the current mesh. See the least square methods documentation:

https://www.openfoam.com/documentation/guides/latest/api/classFoam_1_1extendedCentredCellToFaceStencil.html https://www.openfoam.com/documentation/guides/latest/api/classFoam_1_1fv_1_1LeastSquaresGrad.html

The gradient operator constructs and registers its connectivity information in the mesh database (registry) and looks it up with every construction of the operator. This greatly improves flexibility and makes the design more modular. It also defines the role for unstructuredMesh:

* To provide the data needed to compute/store the connectivity.

And for operators:

* To compute the numerical operator.
* To store/lookup the connectivity information required for this operation.

The unstructuredMesh would have a member class:

StencilDataBase

that stores the stencil with the following data structure:

std::unordered_map<std::string, std::any> stencilDB_;

The lookup in an operator would look something like this:

const auto& cellFaceStencil = cellFaceStencil::New(mesh) // construct if not available look up if possible

[HenningScheufler]

What are your thoughts about this? @bevanwsjones @greole

[greole]

I guess, that means we need an approach for object registries, before we can finalize this? Also, for my understanding we should assume that different operators might need different stencils, that way we could implement low and high order schemes next to each other if the stencil implementation is flexible enough. That also raises the question about interpolation schemes since they would define the requirement for a different stencil implementation.

[bevanwsjones]

So just some thoughts,

With this implementation, the mesh's job is to provide the connectivity information necessary to access neighbouring cells.

This is correct, although it is also the mesh's job to provide the geometric data for the mesh including points, volumes, etc. Like the stencils, however, some of this data is similar between spacial discretisation methods but some data might be unique, being a bit naive on the others I can't say for sure. However, the approach taken here with the stencils should work for the geometric data as well.

Consequently, the mesh class would become very large, resulting in a massive interface.

Unless you DOP the design, in which case the mesh is a struct and all data is public. However, then you would have a laundry list of data members, only some of whom would be 'active' at any given time. This also sounds like a management nightmare.

• To provide the data needed to compute/store the connectivity.

We would need also to store geometric data, as above.

Overall, this is a very good approach for flexibility, and I cannot see any performance issues either. We need to take care with boundaries as well, but I can't see it not working.

[bevanwsjones]

disclaimer, it was getting late so I have not completely gone through all of this, but perhaps it makes sense. Will try edit it more tomorrow

TLDR: In addition to Gregor's comment, I would also suggest to more concretely 'define' some terms, and the underlying data structures and naming that would be needed to enable this.

From here on in it's a bit of a ramble, I apologise, also if there are people who have solved the below please feel free to correct me.

Ok to start, and just so we speak the same language, I want to define some terms more concretely. So connectivity refers to the relationship data between various mesh entities. For example, it specifies which faces or cells a particular cell is connected to; for example, cell 10 might be connected to faces 33 and 134. On the other hand, a stencil represents local connectivity information, rebased to 0. It provides the local connectivity details for a given "subject" entity, akin to the way we notate discrete mathematics for a general cell indexed as 'i'. This stencil facilitates our equations by translating the local index to a global one (in a single-rank sense). Although, in many cases, this is implicit in the way FVM layers out its connectivity data, so most don't realise they are using these stencils. Below is a first-order cell-cell example:

std::vector<std::vector<std::size_t> > cell_cell_connectivity;  // for a row-major 3x3 [[1, 3], [0, 2, 4], [1, 5] ... [5, 7]]

// The 1st order cell stencil is then obtained by extracting the relevant row from the connectivity data structure.
const auto& cell_1_stencil = cell_cell_connectivity[1];  // for the stencil of cell 1.

// To get the nth stencil connection we get
cont auto connecting_0th_cell = cell_1_stencil[0];  // 0
cont auto connecting_1st_cell = cell_1_stencil[1];  // 2

And now, to complicate matters, let's extend this for second order. First, we need to expand the cell_cell_connectivity because it's unlikely that we would perform searches for every access. But where and how should we store this additional data? Do we keep it in a separate cell_cell_connectivity_2 vector or append it to the existing data? Just to illustrate this a bit more concretely:

// option 1:
std::vector<std::vector<std::size_t> > cell_cell_connectivity_1;  // for a row-major 3x3 [[1, 3], [0, 2, 4], [1, 5] ... [5, 7]]
std::vector<std::vector<std::size_t> > cell_cell_connectivity_2;  // for a row-major 3x3 [[2, 4, 6], [3, 7, 5], [0, 4, 8] ... [2, 4, 6]]

// The 1st order cell stencil is then obtained by extracting the relevant row from the connectivity data structure.
const auto& cell_1_stencil = {cell_cell_connectivity_1[1], cell_cell_connectivity_2[1]};  // for the stencil of cell 1. 

// To get the nth stencil connection we get
cont auto 0th_connecting_cell = cell_1_stencil[0];  // 0
cont auto 1th_connecting_cell = cell_1_stencil[3];  // 3  - the stencil is then automatically switches to  cell_cell_connectivity_2 lookup.

// option 2:
std::vector<std::vector<std::size_t> > cell_cell_connectivity;  // for a row-major 3x3 [[1, 3, 2, 4, 6], [0, 2, 4, 3, 7, 5], [1, 5, 0, 4, 8] ... [5, 7, 2, 4, 6]]; // what is the ordering conversion here? do we just append or sort? 

// The 1st order cell stencil is then obtained by extracting the relevant row from the connectivity data structure.
const auto& cell_1_stencil = cell_cell_connectivity;  // for the stencil of cell 1. 

// To get the nth stencil connection we get
cont auto 0th_connecting_cell = cell_1_stencil[0];  // 0
cont auto 1th_connecting_cell = cell_1_stencil[3];  // 3  - the stencil just walks through

There are various pros and cons to the above methods. I would probably favor option 2. With option 1, you are effectively performing range-checking every time you access to determine when you need to switch. However, option 1 is advantageous in terms of reusing the first cell_cell_connectivity_1 connectivity data for the remaining 1st-order stencils (although, is this still necessary in high order?). Conversely, if you need 1st order stencils, you would require an additional vector to indicate where each order's data ends. Option 1's caching capability is also not optimal.

However, all of this discussion pertains only to the Finite Volume Method (FVM) approach, where there are single nodes at cell centres, faces, etc. With high-order methods (although I'm not deeply familiar with them), each cell will have multiple nodes on its edges, faces, volumes, etc.

I'm not suggesting that any of these challenges are insurmountable, but they do require careful consideration. Of course, if this flexible approach has already been implemented elsewhere, that would be cool and I look forward to being educated :).

[greole]

@bevanwsjones This is a very valuable discussion. We should copy the definitions to our documentation.

[greole]

Here are some thoughts after re-reading @bevanwsjones comment:

1. Do we have any literature on this?
2. Can we do a quick survey how other codes are handling stencils?
3. Regarding std::vector<std::vector<std::size_t> > cell_cell_connectivity_1 i guess the vector is for illustration, otherwise wouldn't it better to a consecutive array for the indices and a second array to mark the end of the stencils?
4. For my understanding, in your example the second order comes from the fact that we can now address the fourth element [3] the element has only three direct neighbours?

[greole]

Talking to a former colleague, he recommended me some papers on higher order fvm:

1. https://arc.aiaa.org/doi/10.2514/6.2009-954
2. https://doi.org/10.1016/j.jcp.2013.09.045
3. https://doi.org/10.1016/j.camwa.2016.06.024

Maybe something worth looking at.

[bevanwsjones]

I will clean up my original post later this evening, just to make it more readable for those coming later.

Regarding std::vector<std::vector > cell_cell_connectivity_1 i guess the vector is for illustration, otherwise wouldn't it better to a consecutive array for the indices and a second array to mark the end of the stencils?

Yes, it's more for illustration. I should probably get back to the Adjacency Graph stuff at some point. Although, Henning and I were discussing it; perhaps it's better to generalize the data container and call it a 'flat table' or something similar. Either way, I think we should use something like it for connectivity when we get around to it.

For my understanding, in your example the second order comes from the fact that we can now address the fourth element [3] the element has only three direct nieghbours?

Yes, exactly. So, we have a stencil where you can access all 2nd order (in a stencil sense) stencil cells.

I just quickly skimmed through the third link. I don't have access to the others. Here are the stencils they use:

[https://doi.org/10.1016/j.camwa.2016.06.024] Here, they used the two highlighted stencils, but you can see the other options. For 3rd order, they used the bottom right, neighbour's neighbours. This is where the extension of connectivity comes in.

Note: On each face, they have two integration points (Gauss integral points).

[https://doi.org/10.1016/j.camwa.2016.06.024]

This is why high order requires a different description of cell points because not all points are vertices."

[HenningScheufler]

Regarding the design of the operators, every operator can register or read additional connectivity information in the stencilDB:

linear::linear(const executor& exec, const unstructuredMesh& mesh)
    : surfaceInterpolationKernel(exec, mesh),
      mesh_(mesh),
      geometryScheme_(FvccGeometryScheme::readOrCreate(mesh))
{

};

The connectivity information (here FvccGeometryScheme) register a shared_ptr inside of the stencilDB or lookup the class if it is already available.

const std::shared_ptr<FvccGeometryScheme> FvccGeometryScheme::readOrCreate(const unstructuredMesh& uMesh)
{
    StencilDataBase& stencil_db = uMesh.stencilDB();
    if (!stencil_db.contains("FvccGeometryScheme"))
    {
        stencil_db.insert(std::string("FvccGeometryScheme"), std::make_shared<FvccGeometryScheme>(uMesh));
    }
    return stencil_db.get<std::shared_ptr<FvccGeometryScheme>>("FvccGeometryScheme");
}

So the operator would contain all the data necessary to perform its calculation. The stencilDB just gives the option to avoid the reconstruction of the connectivity with construction of the object