Ease the interleave data process by creating empty arrays on all sibling objects

[jgostick]

At the moment if you have two physics objects, phys1 and phys2`, and you create phys1['pore.blah'] = 1.0, then you do ``phys2['pore.blah']`` it will do a BUNCHA stuff behind the scenes to eventually return nans. The process is complicated and take about 200 lines.

What if, instead, upon the creation of phys1['pore.blah'] = 1.0 we automatically create phys2['pore.blah'] = np.nan. This would literally reduce the size of the interleave_data to just a few lines.

Downside: creating new arrays requires checking all the other objects...but the setitem method already does a tonne of that stuff anyway.

[jgostick]

An alternative is to only store data on 'parent' or 'full domain' objects (i.e. network and phases), and the subdomain objects just 'look-up' data from the parent as needed, knowing what indices its associated with.

In this scenario, phys and geom would store NO data on themselves at all. Their parent objects would have all the data, including the label indicating where the subdomains apply (which is already how things work). This label would literally be used as a mask to fetch the data for the subdomain, like phys['pore.blah'] would actually call phase['pore.blah']['pore.' + phys.name'].

When calling phys.num_pores() it would actually do phase['pore.' + phys.name].sum(). And phys.pores(target=phase) would do np.where(phase['pore.' + phys.name])[0].

The main downside to this is that you can no longer add labels to the subdomain objects. I am 99.9% sure nobody has ever used this feature. I certainly haven't.

I think this is a radical change, but is worth seriously considering.

[jgostick]

BTW, I should add that interleave_data is now NON-EXISTENT. It's just phase['pore.blah'] and you literally fetch the whole actual array.

[jgostick]

Here is a mock-up of this approach:

import numpy as np


class Domain(dict):

    def __init__(self, N):
        self.N = N

    def __setitem__(self, key, value):
        if np.isscalar(value):
            value = np.ones(self.N, type(value))*value
        val = np.atleast_1d(value)
        assert val.shape[0] == self.N, "Received array is wrong length"
        super().__setitem__(key, val)

    def create_array(self, key, value, dtype=None):
        val = np.array(value)
        if dtype is None:
            dtype = bool if val.dtype == bool else float
        if val.ndim == 1:
            arr = np.zeros(self.N, dtype=dtype)
        else:
            arr = np.zeros([self.N, val.shape[1]], dtype=dtype)
        self[key] = arr


class Sub:

    def __init__(self, domain, name, locs):
        self.name = name
        self.domain = domain
        mask = np.zeros(domain.N, dtype=bool)
        mask[locs] = True
        self.domain[name] = mask

    def __setitem__(self, key, value):
        mask = self.domain[self.name]
        try:
            self.domain[key][mask] = value
        except KeyError:
            self.domain.create_array(key, value)
            self.__setitem__(key, value)

    def __getitem__(self, key):
        mask = self.domain[self.name]
        arr = self.domain[key][mask]
        return arr


pn = Domain(N=4)
geo = Sub(domain=pn, name='foo', locs=[0, 2])
geo['test'] = [True, False]
geo['bar'] = [1, 1]
pn['blah'] = [3, 2, 2, 3]
print(pn)

[jgostick]

We have decided to go with the latter option, so that all data is stored on the network/phase (i.e. FullDomain), and the geometry/physics (i.e. SubDomain) just reads/writes to the respective FullDomain. This means that when a value is written to the SubDomain, the FULL array is created on the FullDomain but only partially filled with values. Interleaving is then no longer necessary as the FULL (i.e. interleaved) array is now already in existence.