diff --git a/mesh/array_indexer.py b/mesh/array_indexer.py index 200c04c94..8f4eb4559 100644 --- a/mesh/array_indexer.py +++ b/mesh/array_indexer.py @@ -26,9 +26,11 @@ def _buf_split(b): class ArrayIndexer(np.ndarray): - """ a class that wraps the data region of a single array (d) - and allows us to easily do array operations like d[i+1,j] - using the ip() method. """ + """a class that wraps the data region of a single cell-centered data + array (d) and allows us to easily do array operations like + d[i+1,j] using the ip() method. + + """ def __new__(self, d, grid=None): obj = np.asarray(d).view(self) @@ -328,3 +330,244 @@ def pretty_print(self, n=0, fmt=None, show_ghost=True): +---> x """ print(leg) + + +class ArrayIndexerFC(ArrayIndexer): + """a class that wraps the data region of a single face-centered data + array (d) and allows us to easily do array operations like + d[i+1,j] using the ip() method. + + """ + + def __new__(self, d, idir, grid=None): + obj = np.asarray(d).view(self) + obj.g = grid + obj.idir = idir + obj.c = len(d.shape) + + return obj + + def __array_finalize__(self, obj): + if obj is None: + return + self.g = getattr(obj, "g", None) + self.idir = getattr(obj, "idir", None) + self.c = getattr(obj, "c", None) + + def __array_wrap__(self, out_arr, context=None): + return np.ndarray.__array_wrap__(self, out_arr, context) + + def ip_jp(self, ishift, jshift, buf=0, n=0, s=1): + """return a view of the data shifted by ishift in the x direction and + jshift in the y direction. By default the view is the same + size as the valid region, but the buf can specify how many + ghost cells on each side to include. The component is n and s + is the stride + + """ + bxlo, bxhi, bylo, byhi = _buf_split(buf) + c = len(self.shape) + + if self.idir == 1: + # face-centered in x + if c == 2: + return np.asarray(self[self.g.ilo-bxlo+ishift:self.g.ihi+2+bxhi+ishift:s, + self.g.jlo-bylo+jshift:self.g.jhi+1+byhi+jshift:s]) + else: + return np.asarray(self[self.g.ilo-bxlo+ishift:self.g.ihi+2+bxhi+ishift:s, + self.g.jlo-bylo+jshift:self.g.jhi+1+byhi+jshift:s, n]) + elif self.idir == 2: + # face-centered in y + if c == 2: + return np.asarray(self[self.g.ilo-bxlo+ishift:self.g.ihi+1+bxhi+ishift:s, + self.g.jlo-bylo+jshift:self.g.jhi+2+byhi+jshift:s]) + else: + return np.asarray(self[self.g.ilo-bxlo+ishift:self.g.ihi+1+bxhi+ishift:s, + self.g.jlo-bylo+jshift:self.g.jhi+2+byhi+jshift:s, n]) + + def lap(self, n=0, buf=0): + raise NotImplementedError("lap not implemented for ArrayIndexerFC") + + def norm(self, n=0): + """ + find the norm of the quantity (index n) defined on the same grid, + in the domain's valid region + + """ + c = len(self.shape) + if self.idir == 1: + if c == 2: + return np.sqrt(self.g.dx * self.g.dy * + np.sum((self[self.g.ilo:self.g.ihi+2, self.g.jlo:self.g.jhi+1]**2).flat)) + + else: + _tmp = self[:, :, n] + return np.sqrt(self.g.dx * self.g.dy * + np.sum((_tmp[self.g.ilo:self.g.ihi+2, self.g.jlo:self.g.jhi+1]**2).flat)) + elif self.idir == 2: + if c == 2: + return np.sqrt(self.g.dx * self.g.dy * + np.sum((self[self.g.ilo:self.g.ihi+1, self.g.jlo:self.g.jhi+2]**2).flat)) + + else: + _tmp = self[:, :, n] + return np.sqrt(self.g.dx * self.g.dy * + np.sum((_tmp[self.g.ilo:self.g.ihi+1, self.g.jlo:self.g.jhi+2]**2).flat)) + + def copy(self): + """make a copy of the array, defined on the same grid""" + return ArrayIndexer(np.asarray(self).copy(), self.idir, grid=self.g) + + def is_symmetric(self, nodal=False, tol=1.e-14, asymmetric=False): + """return True is the data is left-right symmetric (to the tolerance + tol) For node-centered data, set nodal=True + + """ + raise NotImplementedError() + + def is_asymmetric(self, nodal=False, tol=1.e-14): + """return True is the data is left-right asymmetric (to the tolerance + tol)---e.g, the sign flips. For node-centered data, set nodal=True + + """ + raise NotImplementedError() + + def fill_ghost(self, n=0, bc=None): + """Fill the boundary conditions. This operates on a single component, + n. We do periodic, reflect-even, reflect-odd, and outflow + + We need a BC object to tell us what BC type on each boundary. + """ + + # there is only a single grid, so every boundary is on + # a physical boundary (except if we are periodic) + + # Note: we piggy-back on outflow and reflect-odd for + # Neumann and Dirichlet homogeneous BCs respectively, but + # this only works for a single ghost cell + + # -x boundary + if bc.xlb in ["outflow", "neumann", "reflect-even", "reflect-odd", "dirichlet"]: + raise NotImplementedError("boundary condition not implemented for -x") + elif bc.xlb == "periodic": + if self.idir == 1: + # face-centered in x + for i in range(self.g.ilo): + self[i, :, n] = self[self.g.ihi-self.g.ng+i+1, :, n] + elif self.idir == 2: + # face-centered in y + for i in range(self.g.ilo): + self[i, :, n] = self[self.g.ihi-self.g.ng+i+1, :, n] + + # +x boundary + if bc.xrb in ["outflow", "neumann", "reflect-even", "reflect-odd", "dirichlet"]: + raise NotImplementedError("boundary condition not implemented for +x") + elif bc.xrb == "periodic": + if self.idir == 1: + # face-centered in x + for i in range(self.g.ihi+2, 2*self.g.ng + self.g.nx + 1): + self[i, :, n] = self[i-self.g.ihi-1+self.g.ng, :, n] + elif self.idir == 2: + # face-centered in y + for i in range(self.g.ihi+1, 2*self.g.ng + self.g.nx): + self[i, :, n] = self[i-self.g.ihi-1+self.g.ng, :, n] + + # -y boundary + if bc.ylb in ["outflow", "neumann", "reflect-even", "reflect-odd", "dirichlet"]: + raise NotImplementedError("boundary condition not implemented for -y") + elif bc.ylb == "periodic": + if self.idir == 1: + # face-centered in x + for j in range(self.g.jlo): + self[:, j, n] = self[:, self.g.jhi-self.g.ng+j+1, n] + elif self.idir == 2: + # face-centered in y + for j in range(self.g.jlo): + self[:, j, n] = self[:, self.g.jhi-self.g.ng+j+1, n] + + # +y boundary + if bc.yrb in ["outflow", "neumann", "reflect-even", "reflect-odd", "dirichlet"]: + raise NotImplementedError("boundary condition not implemented for +y") + elif bc.yrb == "periodic": + if self.idir == 1: + # face-centered in x + for j in range(self.g.jhi+1, 2*self.g.ng + self.g.ny): + self[:, j, n] = self[:, j-self.g.jhi-1+self.g.ng, n] + elif self.idir == 2: + for j in range(self.g.jhi+2, 2*self.g.ng + self.g.ny + 1): + self[:, j, n] = self[:, j-self.g.jhi-1+self.g.ng, n] + + def pretty_print(self, n=0, fmt=None, show_ghost=True): + """ + Print out a small dataset to the screen with the ghost cells + a different color, to make things stand out + """ + + if fmt is None: + if self.dtype == np.int: + fmt = "%4d" + elif self.dtype == np.float64: + fmt = "%10.5g" + else: + raise ValueError("ERROR: dtype not supported") + + # print j descending, so it looks like a grid (y increasing + # with height) + if show_ghost: + if self.idir == 1: + ilo = 0 + ihi = self.g.qx + jlo = 0 + jhi = self.g.qy-1 + elif self.idir == 2: + ilo = 0 + ihi = self.g.qx-1 + jlo = 0 + jhi = self.g.qy + + else: + if self.idir == 1: + ilo = self.g.ilo + ihi = self.g.ihi+1 + jlo = self.g.jlo + jhi = self.g.jhi + elif self.idir == 2: + ilo = self.g.ilo + ihi = self.g.ihi + jlo = self.g.jlo + jhi = self.g.jhi+1 + + for j in reversed(range(jlo, jhi+1)): + for i in range(ilo, ihi+1): + + if self.idir == 1: + if (j < self.g.jlo or j > self.g.jhi or + i < self.g.ilo or i > self.g.ihi+1): + gc = 1 + else: + gc = 0 + elif self.idir == 2: + if (j < self.g.jlo or j > self.g.jhi+1 or + i < self.g.ilo or i > self.g.ihi): + gc = 1 + else: + gc = 0 + + if self.c == 2: + val = self[i, j] + else: + val = self[i, j, n] + + if gc: + print("\033[31m" + fmt % (val) + "\033[0m", end="") + else: + print(fmt % (val), end="") + + print(" ") + + leg = """ + ^ y + | + +---> x + """ + print(leg) diff --git a/mesh/face-centered-data-examples.ipynb b/mesh/face-centered-data-examples.ipynb new file mode 100644 index 000000000..0bf3dc134 --- /dev/null +++ b/mesh/face-centered-data-examples.ipynb @@ -0,0 +1,555 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# face-centered data examples" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "The `FaceCenterData2d` represents data that is centered on faces. Here we explore the methods provided" + ] + }, + { + "cell_type": "code", + "execution_count": 1, + "metadata": {}, + "outputs": [], + "source": [ + "from __future__ import print_function\n", + "import numpy as np\n", + "import mesh.boundary as bnd\n", + "import mesh.patch as patch\n", + "import matplotlib.pyplot as plt\n", + "%matplotlib inline\n", + "\n", + "# for unit testing, we want to ensure the same random numbers\n", + "np.random.seed(100)" + ] + }, + { + "cell_type": "code", + "execution_count": 2, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "2-d grid: nx = 4, ny = 6, ng = 2\n" + ] + } + ], + "source": [ + "g = patch.Grid2d(4, 6, ng=2)\n", + "print(g)" + ] + }, + { + "cell_type": "code", + "execution_count": 3, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "BCs: -x: periodic +x: periodic -y: periodic +y: periodic\n" + ] + } + ], + "source": [ + "bc = bnd.BC(xlb=\"periodic\", xrb=\"periodic\",\n", + " ylb=\"periodic\", yrb=\"periodic\")\n", + "print(bc)" + ] + }, + { + "cell_type": "code", + "execution_count": 4, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "fc data: idir = 1, nx = 4, ny = 6, ng = 2\n", + " nvars = 1\n", + " variables:\n", + " a: min: 0.0000000000 max: 0.0000000000\n", + " BCs: -x: periodic +x: periodic -y: periodic +y: periodic \n", + "\n" + ] + } + ], + "source": [ + "# create data that is centered on x-faces\n", + "d = patch.FaceCenterData2d(g, 1)\n", + "d.register_var(\"a\", bc)\n", + "d.create()\n", + "print(d)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Now we fill the grid with random data. `get_var()` returns an `ArrayIndexer` object that has methods for accessing views into the data. Here we use `a.v()` to get the \"valid\" region, i.e. excluding ghost cells." + ] + }, + { + "cell_type": "code", + "execution_count": 5, + "metadata": {}, + "outputs": [], + "source": [ + "a = d.get_var(\"a\")\n", + "a.v()[:,:] = np.random.rand(g.nx+1, g.ny)" + ] + }, + { + "cell_type": "code", + "execution_count": 6, + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "mesh.array_indexer.ArrayIndexerFC" + ] + }, + "execution_count": 6, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "type(a)" + ] + }, + { + "cell_type": "code", + "execution_count": 7, + "metadata": {}, + "outputs": [], + "source": [ + "# since we are nodal in x, it should be the case that ilo and ihi+1 are identical\n", + "a[a.g.ihi+1,:] = a[a.g.ilo,:]" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "when we pretty_print() the variable, we see the ghost cells colored red. Note that we just filled the interior above." + ] + }, + { + "cell_type": "code", + "execution_count": 8, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "\u001b[31m 0\u001b[0m\u001b[31m 0\u001b[0m\u001b[31m 0\u001b[0m\u001b[31m 0\u001b[0m\u001b[31m 0\u001b[0m\u001b[31m 0\u001b[0m\u001b[31m 0\u001b[0m\u001b[31m 0\u001b[0m\u001b[31m 0\u001b[0m \n", + "\u001b[31m 0\u001b[0m\u001b[31m 0\u001b[0m\u001b[31m 0\u001b[0m\u001b[31m 0\u001b[0m\u001b[31m 0\u001b[0m\u001b[31m 0\u001b[0m\u001b[31m 0\u001b[0m\u001b[31m 0\u001b[0m\u001b[31m 0\u001b[0m \n", + "\u001b[31m 0\u001b[0m\u001b[31m 0\u001b[0m 0.12157 0.2092 0.17194 0.33611 0.12157\u001b[31m 0\u001b[0m\u001b[31m 0\u001b[0m \n", + "\u001b[31m 0\u001b[0m\u001b[31m 0\u001b[0m 0.0047189 0.89132 0.81168 0.81765 0.0047189\u001b[31m 0\u001b[0m\u001b[31m 0\u001b[0m \n", + "\u001b[31m 0\u001b[0m\u001b[31m 0\u001b[0m 0.84478 0.57509 0.97862 0.94003 0.84478\u001b[31m 0\u001b[0m\u001b[31m 0\u001b[0m \n", + "\u001b[31m 0\u001b[0m\u001b[31m 0\u001b[0m 0.42452 0.13671 0.2197 0.4317 0.42452\u001b[31m 0\u001b[0m\u001b[31m 0\u001b[0m \n", + "\u001b[31m 0\u001b[0m\u001b[31m 0\u001b[0m 0.27837 0.82585 0.10838 0.27407 0.27837\u001b[31m 0\u001b[0m\u001b[31m 0\u001b[0m \n", + "\u001b[31m 0\u001b[0m\u001b[31m 0\u001b[0m 0.5434 0.67075 0.18533 0.81622 0.5434\u001b[31m 0\u001b[0m\u001b[31m 0\u001b[0m \n", + "\u001b[31m 0\u001b[0m\u001b[31m 0\u001b[0m\u001b[31m 0\u001b[0m\u001b[31m 0\u001b[0m\u001b[31m 0\u001b[0m\u001b[31m 0\u001b[0m\u001b[31m 0\u001b[0m\u001b[31m 0\u001b[0m\u001b[31m 0\u001b[0m \n", + "\u001b[31m 0\u001b[0m\u001b[31m 0\u001b[0m\u001b[31m 0\u001b[0m\u001b[31m 0\u001b[0m\u001b[31m 0\u001b[0m\u001b[31m 0\u001b[0m\u001b[31m 0\u001b[0m\u001b[31m 0\u001b[0m\u001b[31m 0\u001b[0m \n", + "\n", + " ^ y\n", + " |\n", + " +---> x\n", + " \n" + ] + } + ], + "source": [ + "a.pretty_print()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "`pretty_print()` can also take an argumet, specifying the format string to be used for the output." + ] + }, + { + "cell_type": "code", + "execution_count": 9, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "\u001b[31m 0\u001b[0m\u001b[31m 0\u001b[0m\u001b[31m 0\u001b[0m\u001b[31m 0\u001b[0m\u001b[31m 0\u001b[0m\u001b[31m 0\u001b[0m\u001b[31m 0\u001b[0m\u001b[31m 0\u001b[0m\u001b[31m 0\u001b[0m \n", + "\u001b[31m 0\u001b[0m\u001b[31m 0\u001b[0m\u001b[31m 0\u001b[0m\u001b[31m 0\u001b[0m\u001b[31m 0\u001b[0m\u001b[31m 0\u001b[0m\u001b[31m 0\u001b[0m\u001b[31m 0\u001b[0m\u001b[31m 0\u001b[0m \n", + "\u001b[31m 0\u001b[0m\u001b[31m 0\u001b[0m 0.122 0.209 0.172 0.336 0.122\u001b[31m 0\u001b[0m\u001b[31m 0\u001b[0m \n", + "\u001b[31m 0\u001b[0m\u001b[31m 0\u001b[0m0.00472 0.891 0.812 0.8180.00472\u001b[31m 0\u001b[0m\u001b[31m 0\u001b[0m \n", + "\u001b[31m 0\u001b[0m\u001b[31m 0\u001b[0m 0.845 0.575 0.979 0.94 0.845\u001b[31m 0\u001b[0m\u001b[31m 0\u001b[0m \n", + "\u001b[31m 0\u001b[0m\u001b[31m 0\u001b[0m 0.425 0.137 0.22 0.432 0.425\u001b[31m 0\u001b[0m\u001b[31m 0\u001b[0m \n", + "\u001b[31m 0\u001b[0m\u001b[31m 0\u001b[0m 0.278 0.826 0.108 0.274 0.278\u001b[31m 0\u001b[0m\u001b[31m 0\u001b[0m \n", + "\u001b[31m 0\u001b[0m\u001b[31m 0\u001b[0m 0.543 0.671 0.185 0.816 0.543\u001b[31m 0\u001b[0m\u001b[31m 0\u001b[0m \n", + "\u001b[31m 0\u001b[0m\u001b[31m 0\u001b[0m\u001b[31m 0\u001b[0m\u001b[31m 0\u001b[0m\u001b[31m 0\u001b[0m\u001b[31m 0\u001b[0m\u001b[31m 0\u001b[0m\u001b[31m 0\u001b[0m\u001b[31m 0\u001b[0m \n", + "\u001b[31m 0\u001b[0m\u001b[31m 0\u001b[0m\u001b[31m 0\u001b[0m\u001b[31m 0\u001b[0m\u001b[31m 0\u001b[0m\u001b[31m 0\u001b[0m\u001b[31m 0\u001b[0m\u001b[31m 0\u001b[0m\u001b[31m 0\u001b[0m \n", + "\n", + " ^ y\n", + " |\n", + " +---> x\n", + " \n" + ] + } + ], + "source": [ + "a.pretty_print(fmt=\"%7.3g\")" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "now fill the ghost cells -- notice that the left and right are periodic, the upper is outflow, and the lower is reflect, as specified when we registered the data above." + ] + }, + { + "cell_type": "code", + "execution_count": 10, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "\u001b[31m 0.10838\u001b[0m\u001b[31m 0.27407\u001b[0m\u001b[31m 0.27837\u001b[0m\u001b[31m 0.82585\u001b[0m\u001b[31m 0.10838\u001b[0m\u001b[31m 0.27407\u001b[0m\u001b[31m 0.27837\u001b[0m\u001b[31m 0.82585\u001b[0m\u001b[31m 0.10838\u001b[0m \n", + "\u001b[31m 0.18533\u001b[0m\u001b[31m 0.81622\u001b[0m\u001b[31m 0.5434\u001b[0m\u001b[31m 0.67075\u001b[0m\u001b[31m 0.18533\u001b[0m\u001b[31m 0.81622\u001b[0m\u001b[31m 0.5434\u001b[0m\u001b[31m 0.67075\u001b[0m\u001b[31m 0.18533\u001b[0m \n", + "\u001b[31m 0.17194\u001b[0m\u001b[31m 0.33611\u001b[0m 0.12157 0.2092 0.17194 0.33611 0.12157\u001b[31m 0.2092\u001b[0m\u001b[31m 0.17194\u001b[0m \n", + "\u001b[31m 0.81168\u001b[0m\u001b[31m 0.81765\u001b[0m 0.0047189 0.89132 0.81168 0.81765 0.0047189\u001b[31m 0.89132\u001b[0m\u001b[31m 0.81168\u001b[0m \n", + "\u001b[31m 0.97862\u001b[0m\u001b[31m 0.94003\u001b[0m 0.84478 0.57509 0.97862 0.94003 0.84478\u001b[31m 0.57509\u001b[0m\u001b[31m 0.97862\u001b[0m \n", + "\u001b[31m 0.2197\u001b[0m\u001b[31m 0.4317\u001b[0m 0.42452 0.13671 0.2197 0.4317 0.42452\u001b[31m 0.13671\u001b[0m\u001b[31m 0.2197\u001b[0m \n", + "\u001b[31m 0.10838\u001b[0m\u001b[31m 0.27407\u001b[0m 0.27837 0.82585 0.10838 0.27407 0.27837\u001b[31m 0.82585\u001b[0m\u001b[31m 0.10838\u001b[0m \n", + "\u001b[31m 0.18533\u001b[0m\u001b[31m 0.81622\u001b[0m 0.5434 0.67075 0.18533 0.81622 0.5434\u001b[31m 0.67075\u001b[0m\u001b[31m 0.18533\u001b[0m \n", + "\u001b[31m 0.17194\u001b[0m\u001b[31m 0.33611\u001b[0m\u001b[31m 0.12157\u001b[0m\u001b[31m 0.2092\u001b[0m\u001b[31m 0.17194\u001b[0m\u001b[31m 0.33611\u001b[0m\u001b[31m 0.12157\u001b[0m\u001b[31m 0.2092\u001b[0m\u001b[31m 0.17194\u001b[0m \n", + "\u001b[31m 0.81168\u001b[0m\u001b[31m 0.81765\u001b[0m\u001b[31m 0.0047189\u001b[0m\u001b[31m 0.89132\u001b[0m\u001b[31m 0.81168\u001b[0m\u001b[31m 0.81765\u001b[0m\u001b[31m 0.0047189\u001b[0m\u001b[31m 0.89132\u001b[0m\u001b[31m 0.81168\u001b[0m \n", + "\n", + " ^ y\n", + " |\n", + " +---> x\n", + " \n" + ] + } + ], + "source": [ + "d.fill_BC(\"a\")\n", + "a.pretty_print()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "We can find the L2 norm of the data easily" + ] + }, + { + "cell_type": "code", + "execution_count": 11, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + " 0.619671\n" + ] + } + ], + "source": [ + "print(\"{:12.6g}\".format(a.norm()))" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "and the min and max" + ] + }, + { + "cell_type": "code", + "execution_count": 12, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "0.004718856190972565 0.9786237847073697\n" + ] + } + ], + "source": [ + "print(a.min(), a.max())" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## `ArrayIndexer`" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "We we access the data, an `ArrayIndexer` object is returned. The `ArrayIndexer` sub-classes the NumPy `ndarray`, so it can do all of the methods that a NumPy array can, but in addition, we can use the `ip()`, `jp()`, or `ipjp()` methods to the `ArrayIndexer` object shift our view in the x, y, or x & y directions.\n", + "\n", + "To make this clearer, we'll change our data set to be nicely ordered numbers. We index the `ArrayIndex` the same way we would a NumPy array. The index space includes ghost cells, so the `ilo` and `ihi` attributes from the grid object are useful to index just the valid region. The `.v()` method is a shortcut that also gives a view into just the valid data.\n", + "\n", + "Note: when we use one of the `ip()`, `jp()`, `ipjp()`, or `v()` methods, the result is a regular NumPy `ndarray`, not an `ArrayIndexer` object. This is because it only spans part of the domain (e.g., no ghost cells), and therefore cannot be associated with the `Grid2d` object that the `ArrayIndexer` is built from." + ] + }, + { + "cell_type": "code", + "execution_count": 13, + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "mesh.array_indexer.ArrayIndexerFC" + ] + }, + "execution_count": 13, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "type(a)" + ] + }, + { + "cell_type": "code", + "execution_count": 14, + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "numpy.ndarray" + ] + }, + "execution_count": 14, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "type(a.v())" + ] + }, + { + "cell_type": "code", + "execution_count": 15, + "metadata": {}, + "outputs": [], + "source": [ + "a[:,:] = np.arange((g.qx+1)*g.qy).reshape(g.qx+1, g.qy)" + ] + }, + { + "cell_type": "code", + "execution_count": 16, + "metadata": { + "tags": [ + "nbval-ignore-output" + ] + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "\u001b[31m 9\u001b[0m\u001b[31m 19\u001b[0m\u001b[31m 29\u001b[0m\u001b[31m 39\u001b[0m\u001b[31m 49\u001b[0m\u001b[31m 59\u001b[0m\u001b[31m 69\u001b[0m\u001b[31m 79\u001b[0m\u001b[31m 89\u001b[0m \n", + "\u001b[31m 8\u001b[0m\u001b[31m 18\u001b[0m\u001b[31m 28\u001b[0m\u001b[31m 38\u001b[0m\u001b[31m 48\u001b[0m\u001b[31m 58\u001b[0m\u001b[31m 68\u001b[0m\u001b[31m 78\u001b[0m\u001b[31m 88\u001b[0m \n", + "\u001b[31m 7\u001b[0m\u001b[31m 17\u001b[0m 27 37 47 57 67\u001b[31m 77\u001b[0m\u001b[31m 87\u001b[0m \n", + "\u001b[31m 6\u001b[0m\u001b[31m 16\u001b[0m 26 36 46 56 66\u001b[31m 76\u001b[0m\u001b[31m 86\u001b[0m \n", + "\u001b[31m 5\u001b[0m\u001b[31m 15\u001b[0m 25 35 45 55 65\u001b[31m 75\u001b[0m\u001b[31m 85\u001b[0m \n", + "\u001b[31m 4\u001b[0m\u001b[31m 14\u001b[0m 24 34 44 54 64\u001b[31m 74\u001b[0m\u001b[31m 84\u001b[0m \n", + "\u001b[31m 3\u001b[0m\u001b[31m 13\u001b[0m 23 33 43 53 63\u001b[31m 73\u001b[0m\u001b[31m 83\u001b[0m \n", + "\u001b[31m 2\u001b[0m\u001b[31m 12\u001b[0m 22 32 42 52 62\u001b[31m 72\u001b[0m\u001b[31m 82\u001b[0m \n", + "\u001b[31m 1\u001b[0m\u001b[31m 11\u001b[0m\u001b[31m 21\u001b[0m\u001b[31m 31\u001b[0m\u001b[31m 41\u001b[0m\u001b[31m 51\u001b[0m\u001b[31m 61\u001b[0m\u001b[31m 71\u001b[0m\u001b[31m 81\u001b[0m \n", + "\u001b[31m 0\u001b[0m\u001b[31m 10\u001b[0m\u001b[31m 20\u001b[0m\u001b[31m 30\u001b[0m\u001b[31m 40\u001b[0m\u001b[31m 50\u001b[0m\u001b[31m 60\u001b[0m\u001b[31m 70\u001b[0m\u001b[31m 80\u001b[0m \n", + "\n", + " ^ y\n", + " |\n", + " +---> x\n", + " \n" + ] + } + ], + "source": [ + "a.pretty_print()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "We index our arrays as {i,j}, so x (indexed by i) is the row and y (indexed by j) is the column in the NumPy array. Note that python arrays are stored in row-major order, which means that all of the entries in the same row are adjacent in memory. This means that when we simply print out the `ndarray`, we see constant-x horizontally, which is the transpose of what we are used to." + ] + }, + { + "cell_type": "code", + "execution_count": 17, + "metadata": { + "tags": [ + "nbval-ignore-output" + ] + }, + "outputs": [ + { + "data": { + "text/plain": [ + "array([[22., 23., 24., 25., 26., 27.],\n", + " [32., 33., 34., 35., 36., 37.],\n", + " [42., 43., 44., 45., 46., 47.],\n", + " [52., 53., 54., 55., 56., 57.],\n", + " [62., 63., 64., 65., 66., 67.]])" + ] + }, + "execution_count": 17, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "a.v()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "We can offset our view into the array by one in x -- this would be like {i+1, j} when we loop over data. The `ip()` method is used here, and takes an argument which is the (positive) shift in the x (i) direction. So here's a shift by 1" + ] + }, + { + "cell_type": "code", + "execution_count": 18, + "metadata": { + "tags": [ + "nbval-ignore-output" + ] + }, + "outputs": [ + { + "data": { + "text/plain": [ + "array([[ 1., 2., 3., 4., 5., 6., 7., 8.],\n", + " [11., 12., 13., 14., 15., 16., 17., 18.],\n", + " [21., 22., 23., 24., 25., 26., 27., 28.],\n", + " [31., 32., 33., 34., 35., 36., 37., 38.],\n", + " [41., 42., 43., 44., 45., 46., 47., 48.],\n", + " [51., 52., 53., 54., 55., 56., 57., 58.],\n", + " [61., 62., 63., 64., 65., 66., 67., 68.]])" + ] + }, + "execution_count": 18, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "a.ip(-1, buf=1)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "A shifted view is necessarily smaller than the original array, and relies on ghost cells to bring new data into view. Because of this, the underlying data is no longer the same size as the original data, so we return it as an `ndarray` (which is actually just a view into the data in the `ArrayIndexer` object, so no copy is made.\n", + "\n", + "To see that it is simply a view, lets shift and edit the data" + ] + }, + { + "cell_type": "code", + "execution_count": 19, + "metadata": { + "scrolled": true + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "\u001b[31m 9\u001b[0m\u001b[31m 19\u001b[0m\u001b[31m 29\u001b[0m\u001b[31m 39\u001b[0m\u001b[31m 49\u001b[0m\u001b[31m 59\u001b[0m\u001b[31m 69\u001b[0m\u001b[31m 79\u001b[0m\u001b[31m 89\u001b[0m \n", + "\u001b[31m 8\u001b[0m\u001b[31m 18\u001b[0m\u001b[31m 28\u001b[0m\u001b[31m 38\u001b[0m\u001b[31m 48\u001b[0m\u001b[31m 58\u001b[0m\u001b[31m 68\u001b[0m\u001b[31m 78\u001b[0m\u001b[31m 88\u001b[0m \n", + "\u001b[31m 7\u001b[0m\u001b[31m 17\u001b[0m 27 37 47 57 67\u001b[31m 77\u001b[0m\u001b[31m 87\u001b[0m \n", + "\u001b[31m 6\u001b[0m\u001b[31m 16\u001b[0m 26 36 46 56 66\u001b[31m 76\u001b[0m\u001b[31m 86\u001b[0m \n", + "\u001b[31m 5\u001b[0m\u001b[31m 15\u001b[0m 25 35 45 55 65\u001b[31m 75\u001b[0m\u001b[31m 85\u001b[0m \n", + "\u001b[31m 4\u001b[0m\u001b[31m 14\u001b[0m 24 34 44 54 64\u001b[31m 74\u001b[0m\u001b[31m 84\u001b[0m \n", + "\u001b[31m 3\u001b[0m\u001b[31m 13\u001b[0m 23 33 0 53 63\u001b[31m 73\u001b[0m\u001b[31m 83\u001b[0m \n", + "\u001b[31m 2\u001b[0m\u001b[31m 12\u001b[0m 22 32 42 52 62\u001b[31m 72\u001b[0m\u001b[31m 82\u001b[0m \n", + "\u001b[31m 1\u001b[0m\u001b[31m 11\u001b[0m\u001b[31m 21\u001b[0m\u001b[31m 31\u001b[0m\u001b[31m 41\u001b[0m\u001b[31m 51\u001b[0m\u001b[31m 61\u001b[0m\u001b[31m 71\u001b[0m\u001b[31m 81\u001b[0m \n", + "\u001b[31m 0\u001b[0m\u001b[31m 10\u001b[0m\u001b[31m 20\u001b[0m\u001b[31m 30\u001b[0m\u001b[31m 40\u001b[0m\u001b[31m 50\u001b[0m\u001b[31m 60\u001b[0m\u001b[31m 70\u001b[0m\u001b[31m 80\u001b[0m \n", + "\n", + " ^ y\n", + " |\n", + " +---> x\n", + " \n" + ] + } + ], + "source": [ + "d = a.ip(1)\n", + "d[1,1] = 0.0\n", + "a.pretty_print()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Here, since d was really a view into $a_{i+1,j}$, and we accessed element (1,1) into that view (with 0,0 as the origin), we were really accessing the element (2,1) in the valid region" + ] + } + ], + "metadata": { + "celltoolbar": "Tags", + "kernelspec": { + "display_name": "Python 3", + "language": "python", + "name": "python3" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 3 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython3", + "version": "3.7.2" + } + }, + "nbformat": 4, + "nbformat_minor": 2 +} diff --git a/mesh/patch.py b/mesh/patch.py index fd2c4126d..86260b261 100644 --- a/mesh/patch.py +++ b/mesh/patch.py @@ -376,7 +376,7 @@ def get_var(self, name): return var raise KeyError("name {} is not valid".format(name)) else: - return ai.ArrayIndexer(d=self.data[:, :, n], grid=self.grid) + return self.get_var_by_index(n) def get_var_by_index(self, n): """ @@ -696,6 +696,62 @@ def create(self): self.initialized = 1 + def __str__(self): + """ print out some basic information about the FaceCenterData2d + object """ + + if self.initialized == 0: + my_str = "FaceCenterData2d object not yet initialized" + return my_str + + my_str = "fc data: idir = {}, nx = {}, ny = {}, ng = {}\n".format( + self.idir, self.grid.nx, self.grid.ny, self.grid.ng) + my_str += " nvars = {}\n".format(self.nvar) + my_str += " variables:\n" + + for n in range(self.nvar): + my_str += "%16s: min: %15.10f max: %15.10f\n" % \ + (self.names[n], self.min(self.names[n]), self.max(self.names[n])) + my_str += "%16s BCs: -x: %-12s +x: %-12s -y: %-12s +y: %-12s\n" %\ + (" ", self.BCs[self.names[n]].xlb, + self.BCs[self.names[n]].xrb, + self.BCs[self.names[n]].ylb, + self.BCs[self.names[n]].yrb) + + return my_str + + def get_var_by_index(self, n): + """ + Return a data array for the variable with index n in the + data array. Any changes made to this are automatically + reflected in the CellCenterData2d object. + + Parameters + ---------- + n : int + The index of the variable to access + + Returns + ------- + out : ndarray + The array of data corresponding to the index + + """ + return ai.ArrayIndexerFC(d=self.data[:, :, n], idir=self.idir, grid=self.grid) + + def get_vars(self): + """ + Return the entire data array. Any changes made to this + are automatically reflected in the CellCenterData2d object. + + Returns + ------- + out : ndarray + The array of data + + """ + return ai.ArrayIndexerFC(d=self.data, idir=self.idir, grid=self.grid) + def fill_BC(self, name): """ Fill the boundary conditions. This operates on a single state