GEOS
GEOS copied to clipboard
Implement a structured mesh with an IJK interface.
In order to most effeciently run a wellbore simulation, we should have a structured mesh option with an IJK interface. Key components are:
- Structured mesh in cartesian and cylindrical coordinate frames.
- IJK interface for maps. For instance
elementToNode
map can just be a layer that accesses data using IJK backend.
@corbett5 I am thinking you would be most qualified to do this.
After the refactoring of #1418 we have a CellBlockManagerABC
with the following interface (this is a subset chosen for the discussion).
class CellBlockManagerABC
{
public:
...
virtual array2d< real64, nodes::REFERENCE_POSITION_PERM > getNodePositions() const = 0;
...
};
Here abstraction depend on the implementation of the mappings (that are done as LvArrays
).
For example, if I consider the CellBlockManagerABC::getNodePositions()
member function, the mapping is implemented under the form of an array2d
. But conceptually, it’s a N -> R^3
mapping. So we could imagine implementing some kind of
class NtoR3
{
public:
virtual inputSize() const = 0; // ~ to numNodes
virtual void operator( localIndex n, real64 & x, real64 & y, real64 & z ) const = 0 ;
}
and this class could be implemented, using an array2d< real64, nodes::REFERENCE_POSITION_PERM >
, but another implementation could be using the IJK
information not to store an array but recompute the positions on the fly on the GPU.
And this is for position, but it’s valid for more or less any mapping that we’re interested in.
So we could end up with some code like
class CellBlockManagerABC
{
public:
...
virtual NtoR3 getNodePositions() const = 0;
...
};
Obviously,
- It’s impossible to get
virtual
calls on GPU so some additional dispatch is necessary; putting pressure on compilation and/or JIT tooling. - We'll have to update some managers accordingly.
- We can start with solvers that are using a limited number of mappings (and not only the node mapping).
- We can also build convenience functions (like
NtoR3 to array2d
functions) to be able not to break all the code at once, or to deal with code which does not run on GPU and is not too demanding. - If we do not modify the solvers we'll still use some
for node in nodes
loop and not some triplefor i, for j, for k ...
loop. Which may not be what we want? We'll still use the internalIJK
structure, less memory, ... but maybe it will not be as efficient as we could? Maybe it's better than nothing anyway?