Skip to content

Multidomain Grid Function Spaces

Santiago Ospina De Los Ríos requested to merge feature/multidomain-ordering into master

During the implementation of !555 (merged), I realised that most of the (heavy-)machinery I used to generate subdomain spaces in dune-copasi could be unnecessary if we exposed some of the already implemented features of the ordering library from PDELab.

Why is the current ordering implementation suitable for this?

Let me explain: the ordering library is very good at separating the concerns (offsets) required for every DOF in the GFS tree. In particular, every GFS node tagged with entity blocking (or leaf nodes otherwise) will create a so-called grid view ordering. This grid view ordering generates a vector of offsets that identify each degree of freedom index for a given entity (Remember that an ordering is nothing more than a map from DOFIndex to ContainerIndex). Now, the trick here is that when we compose several GFS nodes that have entity blocking into a GFS node that has no entity blocking (see example below), the grid view ordering for each child node ordering is independent of its sibling nodes. This implies that a complete set of offsets is built for each of those child nodes.

using GFS = GridFunctionSpace<....>; // entity-blocked ordering node
using VBE = ISTL::VectorBackend<>;
using RootGFS = PowerGridFunctionSpace<GFS, 2, VBE, LexicographicOrderingTag>; // non entity-blocked ordering node!

Therefore, multiple domains are native to the ordering library! We just need a different entity set per child node and this is exactly what dune-multidomaingrid does. Eureka!!

What is missing?

BUT not all PDELab is ready for this:

  • The assembly algorithm and vtk writers fetch the entity set and the grid view from root GFS. Currently, this is obtained from the first leaf node. If we do not change this, the assembly would only happen in one part of a bigger domain. The root node shall have another entity set that covers all child nodes.
  • Currently, the local function space exects that every child has a local finite element. This may not be true since some domains (child GFSs) may have no entities nor DOFs in certain regions.
  • The interpolation algorithm tries to interpolate even when there are no degrees of freedom to fill in the coefficient vector. This of course would fail in the case of being out of domain because there is not a finite element to interpolate from.

What does this MR do?

This MR basically solves all the issues listed above and presents a working test for a multidomain grid function space. To have a test, I added dune-multidomaingrid as a suggestion. I am not sure if this is appropriated or not.

How would this work?

Please check out the testpoisson-multidomain.cc file. There, a grid is split into two parts, each part is given to a leaf (entity-blocked) GFS to create a power GFS of degree 2. This root node is assigned with an entity set that holds the whole grid so that assembly passes through all entities. Then, on the assembly, the local operator forwards the correct child space to the base ConvectionDiffusionFEM local operator and everything there works as before.

Is there something that needs to be double checked?

Yes! This MR is still a Draft!

  • Parallel runs with multiple domains has not been checked and is not likely to be working now.

Is this Backwards compatible?

Yes and No:

  • Entity sets are now checked to be the same on nodes within the same entity blocking. Why? Because that's the main assumption of local orderings. However, we only check that they are using the same underlying (shared) index set. We do this, to avoid having different index sets in memory that represent the exact same partition.

  • On the other hand, the typical creation of GFS with grid view creates different (exclusive) entity sets for each node. Doing so and composing them with an entity blocked parent results in an exception. This due to the different underlying index sets. One could argue that this is a backwards compatibility issue but it is also an anti-pattern which very likely reduces a lot of bandwidth for memory bound problems. This is exactly what happens in the fixedsize test: !557 (diffs). The fix is as simple as using the entity set constructor of the GFS.

This could be made fully backwards compatible if it feels too invasive...

Can this MR be accepted?

  • Implemented
    • New ownership of entity set for root GFS
    • Local function spaces accept leaf nodes with no entities attached to them
    • Interpolation is not performed on empty leaf local function spaces
    • Entity indices used in parallel data handler are from leaf nodes instead of from root nodes
  • Added/Updated tests
  • Added/Updated documentation
  • Pipelines passing
  • Added entry to CHANGELOG.md

Related issues

Closes #

Edited by Santiago Ospina De Los Ríos

Merge request reports