TPFA_ResSim.grid

  1"""Tools for working with model grid coordinates.
  2
  3Most functions here are barely in use, mostly serving as reference.
  4After all, it is surprisingly hard to remember
  5which direction and index is for x and which is for y.
  6
  7The index ordering is "C-style" (numpy default).
  8This choice means that `x` is the 1st coord., `y` is 2nd,
  9and is hardcoded in the reservoir simulator model code
 10(in what takes place **between** `np.ravel` and `np.reshape`,
 11both of which are configured to use row-major index ordering.
 12"F-style" (column-major) indexing implementation is perfectly possible,
 13but would imply an undue amount hassle).
 14Conveniently, it also means that `x` and `y` tend to occur in alphabetic order.
 15Thus, in printing a matrix of a field, the `x` coordinate corresponds to the row index.
 16By contrast, the plotting module depicts `x` from left to right, `y` from bottom to top.
 17"""
 18
 19from dataclasses import dataclass
 20
 21import numpy as np
 22
 23
 24@dataclass
 25class Grid2D:
 26    """Defines a 2D rectangular grid.
 27
 28    Example (2 x-nodes, 5 y-nodes):
 29    >>> grid = Grid2D(Lx=6, Ly=10, Nx=3, Ny=5)
 30
 31    The nodes are centered in the cells:
 32    >>> X, Y = grid.mesh
 33    >>> X
 34    array([[1., 1., 1., 1., 1.],
 35           [3., 3., 3., 3., 3.],
 36           [5., 5., 5., 5., 5.]])
 37
 38    You can compute cell boundaries (i.e. non-central nodes) by adding or subtracting
 39    `hx`/2 and `hy`/2 (i.e. you will miss either boundary at 0 or `Lx` or `Ly`).
 40
 41    Test of round-trip capability of grid mapping computations:
 42    >>> ij = (0, 4)
 43    >>> grid.xy2sub(X[ij], Y[ij]) == ij
 44    array([ True,  True])
 45
 46    >>> grid.sub2xy(*ij) == (X[ij], Y[ij])
 47    array([ True,  True])
 48    """
 49
 50    Lx: float = 1.0
 51    """Physical x-length of domain."""
 52    Ly: float = 1.0
 53    """Physical y-length of domain."""
 54    Nx: int = 32
 55    """Number of grid cells (and their centres) in x dir."""
 56    Ny: int = 32
 57    """Number of grid cells (and their centres) in y dir."""
 58
 59    @property
 60    def shape(self):
 61        """`(Nx, Ny)`"""
 62        return self.Nx, self.Ny
 63
 64    @property
 65    def size(self):
 66        """Total number of elements."""
 67        return np.prod(self.shape)
 68
 69    @property
 70    def domain(self):
 71        """`((0, 0), (Lx, Ly))`"""
 72        return ((0, 0), (self.Lx, self.Ly))
 73
 74    @property
 75    def Nxy(self):
 76        """`Nx` * `Ny`"""
 77        return np.prod(self.shape)
 78
 79    @property
 80    def hx(self):
 81        """x-length of cells"""
 82        return self.Lx/self.Nx
 83
 84    @property
 85    def hy(self):
 86        """y-length of cells"""
 87        return self.Ly/self.Ny
 88
 89    @property
 90    def h2(self):
 91        """`hx` * `hy`"""
 92        return self.hx*self.hy
 93
 94    @property
 95    def mesh(self):
 96        """Generate 2D coordinate grid of cell centres."""
 97        xx = np.linspace(0, self.Lx, self.Nx, endpoint=False) + self.hx/2
 98        yy = np.linspace(0, self.Ly, self.Ny, endpoint=False) + self.hy/2
 99        return np.meshgrid(xx, yy, indexing="ij")
100
101    def sub2ind(self, ix, iy):
102        """Convert index `(ix, iy)` to index in flattened array."""
103        idx = np.ravel_multi_index((ix, iy), self.shape)
104        return idx
105
106    def ind2sub(self, ind):
107        """Inv. of `self.sub2ind`."""
108        ix, iy = np.unravel_index(ind, self.shape)
109        return np.asarray([ix, iy])
110
111    def xy2sub(self, x, y):
112        """Convert physical coordinate tuple to `(ix, iy)`, ix ∈ {0, ..., Nx-1}.
113
114        .. warning:: `xy2sub` and `xy2ind` *round* to nearest cell center.
115            I.e. they are not injective.
116            The alternative would be to return some kind
117            of interpolation weights distributing `(x, y)` over multiple nodes.
118        """
119        x = np.asarray(x)
120        y = np.asarray(y)
121        # Don't silence errors! Validation is useful in optimisation (e.g.)
122        assert np.all(x <= self.Lx)
123        assert np.all(y <= self.Ly)
124        # Set upper border values to slightly interior
125        x = x.clip(max=self.Lx-1e-8)
126        y = y.clip(max=self.Ly-1e-8)
127        ix = np.floor(x / self.Lx * self.Nx).astype(int)
128        iy = np.floor(y / self.Ly * self.Ny).astype(int)
129        return np.asarray([ix, iy])
130
131    def xy2ind(self, x, y):
132        """Convert physical coord to flat indx."""
133        i, j = self.xy2sub(x, y)
134        return self.sub2ind(i, j)
135
136    def sub2xy(self, ix, iy):
137        """Inverse of `self.xy2sub`."""
138        x = (np.asarray(ix) + .5) * self.hx
139        y = (np.asarray(iy) + .5) * self.hy
140        return np.asarray([x, y])
141
142    def ind2xy(self, ind):
143        """Inverse of `self.xy2ind`."""
144        i, j = self.ind2sub(ind)
145        return self.sub2xy(i, j)