diff --git a/docs/coordinate_labels.md b/docs/coordinate_labels.md new file mode 100644 index 0000000..20796f4 --- /dev/null +++ b/docs/coordinate_labels.md @@ -0,0 +1,62 @@ +# Labeling of coordinates + + +## User (API) coordinate labels and internal labels + +*If you are a casual user of this library, you do not need to read this section. It details the inner workings of PyFVTool, but this knowledge is not required to simply use the PyFVTool as a library.* + +Internally, PyFVTool always uses an (x, y, z) convention for labeling coordinates, even for cylindrical and spherical grids. This is for historical reasons and efficient coding. Also, these three dimensions are always present, even in the case for 1D and 2D grids. + +To avoid confusing situations on the user (API) side, PyFVTool uses conventional coordinate labels towards the user, and then translates these to the appropriate internal labels. Thus, for 1D Cartesian grids, there is (x), for 1D cylindrical grids, there is \(r\), for 2D Cartesian grids, there is (x, y), for 2D cylindrical grids, there is (r, z), and so on. + +The 'internal' labeling uses a preceding underscore to distinguish it from the 'user' labeling. This follows the Python convention that indicates that these variables (properties) are internal to PyFVTool (private) and should not be touched by the external user. + +Below, the correspondence between the (conventional) user coordinate labels and the internal (underscored) variable names is given for the different meshes. + + +### Cell and mesh properties + +All mesh objects (subclasses of `MeshStructure`) have composite properties `cellsize`, `cellcenters` and `facecenters` which are defined along each of the coordinates defined by the specific grid of the mesh. + +For instance, the mesh `Grid2D` has `Grid2D.cellsize.x`, `Grid2D.cellsize.y` , `Grid2D.cellcenters.x` and so on. These are the coordinates labeled according to the 'user' convention. They correspond, in this case, to internal variables `._x` and `._y` . + +The correspondence between conventional user coordinate labels and the internal variable names is as given in the table. + + +| |`_x`|`_y` |`_z` | +|-------------------|----|-------|-----| +|`Grid1D` |`x` | | | +|`CylindricalGrid1D`|`r` | | | +|`SphericalGrid1D` |`r` | | | +|`Grid2D` |`x` |`y` | | +|`CylindricalGrid2D`|`r` |`z` | | +|`PolarGrid2D` |`r` |`theta`| | +|`Grid3D` |`x` |`y` |`z` | +|`CylindricalGrid3D`|`r` |`theta`|`z` | +|`SphericalGrid3D` |`r` |`theta`|`phi`| + +`SphericalGrid1D` has not yet been implemented. + + + + + +### FaceVariable + +`FaceVariable` objects handle vectorial quantities, defined with respect to the specific mesh coordinate system. Each of the components of the vector is in a separate variable (property) of the object, referred to as xvalue, rvalue and so on. The relation between the conventional user labels of the vector components of the vector and the internal variable names is listed in the table. + +| |`_xvalue`|`_yvalue` |`_zvalue` | +|-------------------|---------|------------|----------| +|`Grid1D` |`xvalue` | | | +|`CylindricalGrid1D`|`rvalue` | | | +|`SphericalGrid1D` |`rvalue` | | | +|`Grid2D` |`xvalue` |`yvalue` | | +|`CylindricalGrid2D`|`rvalue` |`zvalue` | | +|`PolarGrid2D` |`rvalue` |`thetavalue`| | +|`Grid3D` |`xvalue` |`yvalue` |`zvalue` | +|`CylindricalGrid3D`|`rvalue` |`thetavalue`|`zvalue` | +|`SphericalGrid3D` |`rvalue` |`thetavalue`|`phivalue`| + +`SphericalGrid1D` has not yet been implemented. + + diff --git a/examples-notebooks/PyFVTool-introduction-demo.ipynb b/examples-notebooks/PyFVTool-introduction-demo.ipynb index d27f8c1..3a1da71 100644 --- a/examples-notebooks/PyFVTool-introduction-demo.ipynb +++ b/examples-notebooks/PyFVTool-introduction-demo.ipynb @@ -59,21 +59,24 @@ "name": "stdout", "output_type": "stream", "text": [ - "dims : [10]\n", - "cellsize : x : [0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1]\n", - "y : [0.]\n", - "z : [0.]\n", + "dims: [10]\n", + "cellsize: _x: [0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1]\n", + "_y: [0.]\n", + "_z: [0.]\n", + "coordlabels: {'x': '_x'}\n", "\n", - "cellcenters : x : [0.05 0.15 0.25 0.35 0.45 0.55 0.65 0.75 0.85 0.95]\n", - "y : [0.]\n", - "z : [0.]\n", + "cellcenters: _x: [0.05 0.15 0.25 0.35 0.45 0.55 0.65 0.75 0.85 0.95]\n", + "_y: [0.]\n", + "_z: [0.]\n", + "coordlabels: {'x': '_x'}\n", "\n", - "facecenters : x : [0. 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1. ]\n", - "y : [0.]\n", - "z : [0.]\n", + "facecenters: _x: [0. 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1. ]\n", + "_y: [0.]\n", + "_z: [0.]\n", + "coordlabels: {'x': '_x'}\n", "\n", - "corners : [1]\n", - "edges : [1]\n", + "corners: [1]\n", + "edges: [1]\n", "\n" ] } @@ -102,38 +105,45 @@ " | Grid1D(*args)\n", " |\n", " | Mesh based on a 1D Cartesian grid (x)\n", + " | =====================================\n", " |\n", - " | Method resolution order:\n", - " | Grid1D\n", - " | MeshStructure\n", - " | builtins.object\n", + " | This class can be instantiated in different ways: from a list of cell face\n", + " | locations or from the number of cells and domain length.\n", " |\n", - " | Methods defined here:\n", + " | Instantiation Options:\n", + " | ----------------------\n", + " | - Grid1D(Nx, Lx)\n", + " | - Grid1D(face_locationsX)\n", " |\n", - " | __init__(self, *args)\n", - " | Create a Grid1D mesh object from a list of cell face locations or from\n", - " | number of cells and domain length.\n", " |\n", - " | Parameters\n", - " | ----------\n", + " | Parameters\n", + " | ----------\n", + " | Grid1D(Nx, Lx)\n", " | Nx : int\n", " | Number of cells in the x direction.\n", " | Lx : float\n", " | Length of the domain in the x direction.\n", - " | face_locations : ndarray\n", + " |\n", + " | Grid1D(face_locationsX)\n", + " | face_locationsX : ndarray\n", " | Locations of the cell faces in the x direction.\n", " |\n", - " | Returns\n", - " | -------\n", + " | Examples\n", + " | --------\n", + " | >>> import numpy as np\n", + " | >>> from pyfvtool import Grid1D\n", + " | >>> mesh = Grid1D(10, 10.0)\n", + " | >>> print(mesh)\n", + " |\n", + " | Method resolution order:\n", " | Grid1D\n", - " | A 1D mesh object.\n", + " | MeshStructure\n", + " | builtins.object\n", " |\n", - " | Examples\n", - " | --------\n", - " | >>> import numpy as np\n", - " | >>> from pyfvtool import Grid1D\n", - " | >>> mesh = Grid1D(10, 10.0)\n", - " | >>> print(mesh)\n", + " | Methods defined here:\n", + " |\n", + " | __init__(self, *args)\n", + " | Initialize self. See help(type(self)) for accurate signature.\n", " |\n", " | __repr__(self)\n", " | Return repr(self).\n", @@ -146,7 +156,18 @@ " | __str__(self)\n", " | Return str(self).\n", " |\n", - " | shift_origin(self, x=0.0, y=0.0, z=0.0)\n", + " | cellVolumes(self)\n", + " | Get the volumes of all finite volume cells in the mesh\n", + " |\n", + " | Returns\n", + " | -------\n", + " | np.ndarray\n", + " | containing all cell volumes, arranged according to gridcells\n", + " |\n", + " | TODO: these could perhaps be calculated statically, when initializing\n", + " | the mesh.\n", + " |\n", + " | shift_origin(self, _x=0.0, _y=0.0, _z=0.0)\n", " |\n", " | visualize(self)\n", " |\n", @@ -187,13 +208,21 @@ " | __init__(self, dims, cellsize, cellcenters, facecenters, corners, edges)\n", " | Initialize self. See help(type(self)) for accurate signature.\n", " |\n", - " | __repr__(self)\n", - " | Return repr(self).\n", - " |\n", " | __str__(self)\n", " | Return str(self).\n", " |\n", - " | shift_origin(self, x=0.0, y=0.0, z=0.0)\n", + " | cellVolumes(self)\n", + " | Get the volumes of all finite volume cells in the mesh\n", + " |\n", + " | Returns\n", + " | -------\n", + " | np.ndarray\n", + " | containing all cell volumes, arranged according to gridcells\n", + " |\n", + " | TODO: these could perhaps be calculated statically, when initializing\n", + " | the mesh.\n", + " |\n", + " | shift_origin(self, _x=0.0, _y=0.0, _z=0.0)\n", " |\n", " | visualize(self)\n", " |\n", @@ -293,8 +322,14 @@ "Lx, Ly, Lz = float(1.0), float(2.0), float(3.0)\n", "m = pf.Grid3D(Nx, Ny, Nz, Lx, Ly, Lz)\n", "\n", - "X, Y, Z = np.meshgrid(m.cellcenters.x, m.cellcenters.y, m.cellcenters.z, indexing='ij')\n", - "Xf, Yf, Zf = np.meshgrid(m.facecenters.x, m.facecenters.y, m.facecenters.z, indexing='ij')" + "X, Y, Z = np.meshgrid(m.cellcenters.x,\n", + " m.cellcenters.y,\n", + " m.cellcenters.z,\n", + " indexing='ij')\n", + "Xf, Yf, Zf = np.meshgrid(m.facecenters.x,\n", + " m.facecenters.y,\n", + " m.facecenters.z, \n", + " indexing='ij')" ] }, { @@ -303,18 +338,14 @@ "id": "f4c714e3-22ce-4b6d-9056-3b1ab09a5c66", "metadata": {}, "outputs": [ - { - "name": "stdout", - "output_type": "stream", - "text": [ - "x : [0.25 0.75]\n", - "y : [0.33333333 1. 1.66666667]\n", - "z : [0.375 1.125 1.875 2.625]\n" - ] - }, { "data": { - "text/plain": [] + "text/plain": [ + "_x: [0.25 0.75]\n", + "_y: [0.33333333 1. 1.66666667]\n", + "_z: [0.375 1.125 1.875 2.625]\n", + "coordlabels: {'x': '_x', 'y': '_y', 'z': '_z'}" + ] }, "execution_count": 10, "metadata": {}, @@ -391,21 +422,24 @@ "name": "stdout", "output_type": "stream", "text": [ - "domain : dims : [10]\n", - "cellsize : x : [0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1]\n", - "y : [0.]\n", - "z : [0.]\n", + "domain : dims: [10]\n", + "cellsize: _x: [0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1]\n", + "_y: [0.]\n", + "_z: [0.]\n", + "coordlabels: {'x': '_x'}\n", "\n", - "cellcenters : x : [0.05 0.15 0.25 0.35 0.45 0.55 0.65 0.75 0.85 0.95]\n", - "y : [0.]\n", - "z : [0.]\n", + "cellcenters: _x: [0.05 0.15 0.25 0.35 0.45 0.55 0.65 0.75 0.85 0.95]\n", + "_y: [0.]\n", + "_z: [0.]\n", + "coordlabels: {'x': '_x'}\n", "\n", - "facecenters : x : [0. 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1. ]\n", - "y : [0.]\n", - "z : [0.]\n", + "facecenters: _x: [0. 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1. ]\n", + "_y: [0.]\n", + "_z: [0.]\n", + "coordlabels: {'x': '_x'}\n", "\n", - "corners : [1]\n", - "edges : [1]\n", + "corners: [1]\n", + "edges: [1]\n", "\n", "left : a : [1.]\n", "b : [0.]\n", @@ -508,21 +542,24 @@ "name": "stdout", "output_type": "stream", "text": [ - "domain : dims : [3 4]\n", - "cellsize : x : [0.33333333 0.33333333 0.33333333 0.33333333 0.33333333]\n", - "y : [0.5 0.5 0.5 0.5 0.5 0.5]\n", - "z : [0.]\n", + "domain : dims: [3 4]\n", + "cellsize: _x: [0.33333333 0.33333333 0.33333333 0.33333333 0.33333333]\n", + "_y: [0.5 0.5 0.5 0.5 0.5 0.5]\n", + "_z: [0.]\n", + "coordlabels: {'x': '_x', 'y': '_y'}\n", "\n", - "cellcenters : x : [0.16666667 0.5 0.83333333]\n", - "y : [0.25 0.75 1.25 1.75]\n", - "z : [0.]\n", + "cellcenters: _x: [0.16666667 0.5 0.83333333]\n", + "_y: [0.25 0.75 1.25 1.75]\n", + "_z: [0.]\n", + "coordlabels: {'x': '_x', 'y': '_y'}\n", "\n", - "facecenters : x : [0. 0.33333333 0.66666667 1. ]\n", - "y : [0. 0.5 1. 1.5 2. ]\n", - "z : [0.]\n", + "facecenters: _x: [0. 0.33333333 0.66666667 1. ]\n", + "_y: [0. 0.5 1. 1.5 2. ]\n", + "_z: [0.]\n", + "coordlabels: {'x': '_x', 'y': '_y'}\n", "\n", - "corners : [ 0 24 5 29]\n", - "edges : [1]\n", + "corners: [ 0 24 5 29]\n", + "edges: [1]\n", "\n", "left : a : [1. 1. 1. 1.]\n", "b : [0. 0. 0. 0.]\n", @@ -1294,14 +1331,9 @@ "\n", "hfig1, ax1 = plt.subplots()\n", "\n", - "# xx, uu = pyfvm.utilities.get_CellVariable_profile1D(ci[0])\n", - "# ax1.plot(xx, uu, 'b-', label='initial value')\n", "ax1.plot(x, ca[0], 'bo', label='initial analytic value')\n", "ax1.plot(x, ci[0].internalCellValues, 'b-', label='initial numerical value')\n", "\n", - "\n", - "# xx, uu = pyfvm.utilities.get_CellVariable_profile1D(ci[-1])\n", - "# ax1.plot(xx, uu, 'r-', label='final value')\n", "ax1.plot(x, ca[-1], 'ro', label='final analytic value')\n", "ax1.plot(x, ci[-1].internalCellValues, 'r-', label='final numerical value')\n", "\n", @@ -1346,7 +1378,6 @@ " global ci\n", " global ca\n", " global Nx\n", - " # line.set_data(*pyfvm.utilities.get_CellVariable_profile1D(ci[ii]))\n", " line.set_data(x, ci[ii].internalCellValues\n", " circ.set_data(x, ca[ii]) \n", " return line, circ\n", @@ -1816,15 +1847,15 @@ "source": [ "# Solving a 1D diffusion equation with a fixed concentration \n", "# at the left boundary and a closed boundary on the right side\n", - "Nx = 20 # number of finite volume cells\n", - "Lx = 1.0 # [m] length of the domain \n", + "Nr = 20 # number of finite volume cells\n", + "Lr = 1.0 # [m] length of the domain \n", "c_left = 1.0 # left boundary concentration\n", "c_init = 0.0 # initial concentration\n", "D_val = 1e-5 # diffusion coefficient (gas phase)\n", "t_simulation = 3600.0 # [s] simulation time\n", "dt = 60.0 # [s] time step\n", "\n", - "m1 = pf.Grid1D(Nx, Lx) # mesh object\n", + "m1 = pf.Grid1D(Nr, Lr) # mesh object\n", "bc = pf.BoundaryConditions(m1) # Neumann boundary condition by default\n", "\n", "# switch the left boundary to Dirichlet: fixed concentration\n", @@ -1861,6 +1892,14 @@ "metadata": {}, "outputs": [], "source": [] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "aa3638aa-3ec4-4302-92f5-fd91ae031834", + "metadata": {}, + "outputs": [], + "source": [] } ], "metadata": { diff --git a/examples-notebooks/cylindrical1D_diffusion.ipynb b/examples-notebooks/cylindrical1D_diffusion.ipynb index a03a6bc..921d72a 100644 --- a/examples-notebooks/cylindrical1D_diffusion.ipynb +++ b/examples-notebooks/cylindrical1D_diffusion.ipynb @@ -37,8 +37,7 @@ "from pyfvtool import CellVariable\n", "from pyfvtool import transientTerm, diffusionTerm\n", "from pyfvtool import harmonicMean, boundaryConditionsTerm\n", - "from pyfvtool import solveMatrixPDE\n", - "from pyfvtool import get_CellVariable_profile1D" + "from pyfvtool import solveMatrixPDE" ] }, { @@ -274,7 +273,7 @@ "\n", "for i in range(0,1001):\n", " if i in sample_i:\n", - " r, phi = get_CellVariable_profile1D(c)\n", + " r, phi = c.plotprofile()\n", " plt.plot(r, phi)\n", " c_an = np.array([crank522(rr,t,a_val,D_val) for rr in r_an])\n", " plt.plot(r_an,c_an, 'k:')\n", @@ -318,7 +317,7 @@ } ], "source": [ - "r_num, phi_num = get_CellVariable_profile1D(c)\n", + "r_num, phi_num = c.plotprofile()\n", "phi_an_r_num = np.array([crank522(rr, t, a_val, D_val) for rr in r_num])\n", "norm_err = (phi_num - phi_an_r_num)/phi_an_r_num.max()\n", "plt.plot(r_num, norm_err, 'o')\n", @@ -336,6 +335,13 @@ "# It is also a basic benchmark test for correct function of PyFVTool\n", "assert np.alltrue(np.abs(norm_err) < 0.0005)" ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [] } ], "metadata": { diff --git a/examples-notebooks/cylindrical1D_diffusion_source_photothermal.ipynb b/examples-notebooks/cylindrical1D_diffusion_source_photothermal.ipynb index df88350..92d8eef 100644 --- a/examples-notebooks/cylindrical1D_diffusion_source_photothermal.ipynb +++ b/examples-notebooks/cylindrical1D_diffusion_source_photothermal.ipynb @@ -470,8 +470,7 @@ "from pyfvtool import transientTerm, diffusionTerm\n", "from pyfvtool import harmonicMean, boundaryConditionsTerm\n", "from pyfvtool import constantSourceTerm\n", - "from pyfvtool import solveMatrixPDE\n", - "from pyfvtool import get_CellVariable_profile1D" + "from pyfvtool import solveMatrixPDE" ] }, { @@ -541,7 +540,7 @@ }, "outputs": [], "source": [ - "dotqval = alpha*S(mesh.cellcenters.x, P, w)\n", + "dotqval = alpha*S(mesh.cellcenters.r, P, w)\n", "fv_dotq = CellVariable(mesh, dotqval)" ] }, @@ -610,7 +609,7 @@ "OKlegend = False\n", "for i in range(0,101):\n", " if i in sample_i:\n", - " r, phi = get_CellVariable_profile1D(fv_T)\n", + " r, phi = fv_T.plotprofile()\n", " ln0, = plt.plot(r, phi)\n", " ln1, = plt.plot(rr0, DeltaT_Sheldon(rr0, i*deltat, P, w, k, cp, rho, t_c), 'k:')\n", " if not OKlegend:\n", @@ -657,7 +656,7 @@ ], "source": [ "# take last calculated time `t`\n", - "r_num, phi_num = get_CellVariable_profile1D(fv_T)\n", + "r_num, phi_num = fv_T.plotprofile()\n", "phi_an_r_num = DeltaT_Sheldon(r_num, t, P, w, k, cp, rho, t_c)\n", "norm_err = (phi_num - phi_an_r_num)/phi_an_r_num.max()\n", "plt.plot(r_num, norm_err, 'o')\n", @@ -673,6 +672,13 @@ "source": [ "assert np.alltrue(np.abs(norm_err) < 0.0025)" ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [] } ], "metadata": { diff --git a/examples-notebooks/cylindrical1D_diffusion_steady.ipynb b/examples-notebooks/cylindrical1D_diffusion_steady.ipynb index 96c1111..9a07eb1 100644 --- a/examples-notebooks/cylindrical1D_diffusion_steady.ipynb +++ b/examples-notebooks/cylindrical1D_diffusion_steady.ipynb @@ -34,8 +34,7 @@ "from pyfvtool import transientTerm, diffusionTerm\n", "from pyfvtool import harmonicMean, boundaryConditionsTerm\n", "from pyfvtool import constantSourceTerm\n", - "from pyfvtool import solveMatrixPDE\n", - "from pyfvtool import get_CellVariable_profile1D" + "from pyfvtool import solveMatrixPDE" ] }, { @@ -253,7 +252,7 @@ }, "outputs": [], "source": [ - "rnum, Tnum = get_CellVariable_profile1D(T)" + "rnum, Tnum = T.plotprofile()" ] }, { @@ -381,6 +380,13 @@ "source": [ "assert np.alltrue(np.abs(norm_err) < 0.001)" ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [] } ], "metadata": { diff --git a/examples-notebooks/cylindrical2D_convection_Taylor.ipynb b/examples-notebooks/cylindrical2D_convection_Taylor.ipynb index dd2c898..a9145a0 100644 --- a/examples-notebooks/cylindrical2D_convection_Taylor.ipynb +++ b/examples-notebooks/cylindrical2D_convection_Taylor.ipynb @@ -98,7 +98,7 @@ "source": [ "# calculate simple finite-volume integral over r\n", "def integral_dr(phi0):\n", - " v = pf.cellVolume(phi0.domain).internalCellValues\n", + " v = phi0.domain.cellVolumes()\n", " c = phi0.internalCellValues\n", " return (v*c).sum(axis=0)" ] @@ -257,8 +257,8 @@ "metadata": {}, "outputs": [], "source": [ - "rr = msh.cellcenters.x\n", - "zz = msh.facecenters.y" + "rr = msh.cellcenters.r\n", + "zz = msh.facecenters.z" ] }, { @@ -288,8 +288,8 @@ "metadata": {}, "outputs": [], "source": [ - "u.xvalue[:] = 0\n", - "u.yvalue[:] = uu[:, np.newaxis]" + "u.rvalue[:] = 0\n", + "u.zvalue[:] = uu[:, np.newaxis]" ] }, { @@ -311,7 +311,7 @@ ], "source": [ "for i in [1, 10, -1]:\n", - " plt.plot(rr*1e6, u.yvalue[:, i])\n", + " plt.plot(rr*1e6, u.zvalue[:, i])\n", "plt.xlabel('$r$ / µm') \n", "plt.ylabel('$u_z(r)$ / m s$^{-1}$');" ] @@ -417,7 +417,7 @@ } ], "source": [ - "initInt = pf.domainInt(phi)\n", + "initInt = phi.domainIntegral()\n", "print(initInt)" ] }, @@ -502,7 +502,7 @@ } ], "source": [ - "print(t, initInt, pf.domainInt(phi), pf.domainInt(phi_new))" + "print(t, initInt, phi.domainIntegral(), phi_new.domainIntegral())" ] }, { @@ -584,7 +584,7 @@ } ], "source": [ - "print(t, initInt, pf.domainInt(phi), pf.domainInt(phi_new))" + "print(t, initInt, phi.domainIntegral(), phi_new.domainIntegral())" ] }, { @@ -632,8 +632,8 @@ "metadata": {}, "outputs": [], "source": [ - "DX = phi.domain.facecenters.y[loadix0]\n", - "X = phi.domain.facecenters.y[loadix1] - DX\n", + "DX = phi.domain.facecenters.z[loadix0]\n", + "X = phi.domain.facecenters.z[loadix1] - DX\n", "C_0 = phiprofs[0][1][(loadix0+loadix1)//2] # slot#0 contains initial condition" ] }, @@ -672,7 +672,7 @@ " else:\n", " lbl1 = None\n", " lbl2 = None\n", - " plt.plot(phi.domain.cellcenters.y, phiprof, \n", + " plt.plot(phi.domain.cellcenters.z, phiprof, \n", " label=lbl1)\n", " if tprof >= X/umax:\n", " plt.plot(zzz, TaylorA3_vec(zzz-DX, tprof, X, C_0, umax),\n", @@ -705,7 +705,7 @@ "outputs": [], "source": [ "(tprof, phiprof) = phiprofs[-1]\n", - "z_num, c_num = phi.domain.cellcenters.y, phiprof\n", + "z_num, c_num = phi.domain.cellcenters.z, phiprof\n", "c_an_z_num = TaylorA3_vec(z_num-DX, tprof, X, C_0, umax)\n", "norm_err = (c_an_z_num - c_num)/c_an_z_num.max()" ] diff --git a/examples-notebooks/how-to-use-pypardiso-as-external-solver.ipynb b/examples-notebooks/how-to-use-pypardiso-as-external-solver.ipynb index bbd2f5b..ea1ae64 100644 --- a/examples-notebooks/how-to-use-pypardiso-as-external-solver.ipynb +++ b/examples-notebooks/how-to-use-pypardiso-as-external-solver.ipynb @@ -174,7 +174,7 @@ "source": [ "# calculate simple finite-volume integral over r\n", "def integral_dr(phi0):\n", - " v = pf.cellVolume(phi0.domain).internalCellValues\n", + " v = phi0.domain.cellVolumes()\n", " c = phi0.internalCellValues\n", " return (v*c).sum(axis=0)" ] @@ -303,8 +303,8 @@ "metadata": {}, "outputs": [], "source": [ - "rr = msh.cellcenters.x\n", - "zz = msh.facecenters.y" + "rr = msh.cellcenters.r\n", + "zz = msh.facecenters.z" ] }, { @@ -334,8 +334,8 @@ "metadata": {}, "outputs": [], "source": [ - "u.xvalue[:] = 0\n", - "u.yvalue[:] = uu[:, np.newaxis]" + "u.rvalue[:] = 0\n", + "u.zvalue[:] = uu[:, np.newaxis]" ] }, { @@ -431,7 +431,7 @@ "metadata": {}, "outputs": [], "source": [ - "initInt = pf.domainInt(phi)\n", + "initInt = phi.domainIntegral()\n", "# print(initInt)" ] }, @@ -533,7 +533,7 @@ "source": [ "# print(t, initInt, pf.domainInt(phi))\n", "# test conservation of mass\n", - "assert np.isclose(initInt, pf.domainInt(phi))" + "assert np.isclose(initInt, phi.domainIntegral())" ] }, { @@ -601,7 +601,7 @@ } ], "source": [ - "print(t, initInt, pf.domainInt(phi))" + "print(t, initInt, phi.domainIntegral())" ] }, { @@ -655,7 +655,7 @@ " lbl1 = 'FVM'\n", " else:\n", " lbl1 = None\n", - " plt.plot(phi.domain.cellcenters.y, phiprof, \n", + " plt.plot(phi.domain.cellcenters.z, phiprof, \n", " label=lbl1)\n", "plt.xlabel('z / m')\n", "plt.legend();" diff --git a/examples-notebooks/tutorial-wave_equation_1D.ipynb b/examples-notebooks/tutorial-wave_equation_1D.ipynb index ab81ea4..5196a89 100644 --- a/examples-notebooks/tutorial-wave_equation_1D.ipynb +++ b/examples-notebooks/tutorial-wave_equation_1D.ipynb @@ -139,10 +139,10 @@ "\n", "hfig1, ax1 = plt.subplots()\n", "\n", - "xx, uu = pf.get_CellVariable_profile1D(ui[0])\n", + "xx, uu = ui[0].plotprofile()\n", "ax1.plot(xx, uu, 'b-', label='initial value')\n", "\n", - "xx, uu = pf.get_CellVariable_profile1D(ui[-1])\n", + "xx, uu = ui[-1].plotprofile()\n", "ax1.plot(xx, uu, 'r-', label='final value')\n", "\n", "ax1.set_xlim((0, Lx))\n", diff --git a/examples/diffusion_explicit_2d.py b/examples/diffusion_explicit_2d.py index 330a10e..1d21359 100644 --- a/examples/diffusion_explicit_2d.py +++ b/examples/diffusion_explicit_2d.py @@ -9,8 +9,8 @@ # define the domain L = 5.0 # domain length -Nx = 100 # number of cells -meshstruct = pf.CylindricalGrid2D(Nx, Nx, L, L) +N = 100 # number of cells +meshstruct = pf.CylindricalGrid2D(N, N, L, L) BC = pf.BoundaryConditions(meshstruct) # all Neumann boundary condition structure BC.bottom.a[:] = 0.0 BC.bottom.b[:] = 1.0 @@ -18,7 +18,7 @@ BC.top.a[:] = 0.0 BC.top.b[:] = 1.0 BC.top.c[:] = 0.0 # top boundary -x = meshstruct.cellcenters.x +r = meshstruct.cellcenters.r ## define the transfer coeffs D_val = 1.0 alfa = pf.CellVariable(meshstruct, 1.0) @@ -37,11 +37,11 @@ c_old.update_value(c) # analytical solution -c_analytical = 1-erf(x/(2*np.sqrt(D_val*t))) +c_analytical = 1-erf(r/(2*np.sqrt(D_val*t))) plt.figure(1) plt.clf() -plt.plot(x, c.internalCellValues[2,:], 'k', label='PyFVTool') -plt.plot(x, c_analytical, 'r--', label='analytic') +plt.plot(r, c.internalCellValues[2,:], 'k', label='PyFVTool') +plt.plot(r, c_analytical, 'r--', label='analytic') plt.legend() plt.show() diff --git a/examples/spherical1D_diffusion.py b/examples/spherical1D_diffusion.py index 74d5494..5fb503a 100644 --- a/examples/spherical1D_diffusion.py +++ b/examples/spherical1D_diffusion.py @@ -40,13 +40,13 @@ ## Solution variable and define initial condition c_init = 0.0 c_old = pf.CellVariable(m, c_init, BC) -r = c_old.domain.cellcenters.x +r = c_old.domain.cellcenters.r c_old.internalCellValues[r<1.0] = 1.0 # TO DO: find a better syntax for this # calculate volumes of FV cellslices # We use this for demonstrating mass conservation -cellA = m.facecenters.x[0:-1] -cellB = m.facecenters.x[1:] +cellA = m.facecenters.r[0:-1] +cellB = m.facecenters.r[1:] cellvol = 4/3 * pi * (cellB**3 - cellA**3) cellsum = np.sum(cellvol) diff --git a/pyproject.toml b/pyproject.toml index 6ffa0f6..13e7dbe 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -1,7 +1,7 @@ # minimalist pyproject.toml [project] name = "pyfvtool" -version = "0.2.1" +version = "0.3.0" authors = [ { name="Ali A. Eftekhari", email="e.eftekhari@gmail.com" }, { name="Gavin M. Weir"}, diff --git a/src/pyfvtool/__init__.py b/src/pyfvtool/__init__.py index e5f041a..2aa3680 100644 --- a/src/pyfvtool/__init__.py +++ b/src/pyfvtool/__init__.py @@ -1,4 +1,4 @@ -__version__ = "0.2.1" +__version__ = "0.3.0" __author__ = [ "Ali A. Eftekhari", @@ -29,11 +29,8 @@ from .averaging import linearMean, arithmeticMean, upwindMean,\ harmonicMean, geometricMean, tvdMean from .pdesolver import solveMatrixPDE, solvePDE, solveExplicitPDE -from .cell import CellVariable, cellVolume, BC2GhostCells -from .cell import copyCellVariable +from .cell import CellVariable from .cell import funceval, celleval -from .cell import domainInt, domainIntegrate -from .cell import get_CellVariable_profile1D from .cell import cellLocations from .face import FaceVariable, faceeval from .face import faceLocations @@ -54,4 +51,11 @@ from .legacy import createMesh3D from .legacy import createMeshCylindrical3D from .legacy import createMeshSpherical3D + from .legacy import get_CellVariable_profile1D + from .legacy import get_CellVariable_profile2D + from .legacy import get_CellVariable_profile3D + from .legacy import domainInt + from .legacy import cellVolume + + \ No newline at end of file diff --git a/src/pyfvtool/advection.py b/src/pyfvtool/advection.py index a66cbbc..fbb1341 100644 --- a/src/pyfvtool/advection.py +++ b/src/pyfvtool/advection.py @@ -14,34 +14,34 @@ def _upwind_min_max(u: FaceVariable, u_upwind: FaceVariable): if issubclass(type(u.domain), Grid1D): - ux_min = np.copy(u.xvalue) - ux_max = np.copy(u.xvalue) - ux_min[u_upwind.xvalue > 0.0] = 0.0 - ux_max[u_upwind.xvalue < 0.0] = 0.0 + ux_min = np.copy(u._xvalue) + ux_max = np.copy(u._xvalue) + ux_min[u_upwind._xvalue > 0.0] = 0.0 + ux_max[u_upwind._xvalue < 0.0] = 0.0 return ux_min, ux_max elif issubclass(type(u.domain), Grid2D): - ux_min = np.copy(u.xvalue) - ux_max = np.copy(u.xvalue) - uy_min = np.copy(u.yvalue) - uy_max = np.copy(u.yvalue) - ux_min[u_upwind.xvalue > 0.0] = 0.0 - ux_max[u_upwind.xvalue < 0.0] = 0.0 - uy_min[u_upwind.yvalue > 0.0] = 0.0 - uy_max[u_upwind.yvalue < 0.0] = 0.0 + ux_min = np.copy(u._xvalue) + ux_max = np.copy(u._xvalue) + uy_min = np.copy(u._yvalue) + uy_max = np.copy(u._yvalue) + ux_min[u_upwind._xvalue > 0.0] = 0.0 + ux_max[u_upwind._xvalue < 0.0] = 0.0 + uy_min[u_upwind._yvalue > 0.0] = 0.0 + uy_max[u_upwind._yvalue < 0.0] = 0.0 return ux_min, ux_max, uy_min, uy_max elif issubclass(type(u.domain), Grid3D): - ux_min = np.copy(u.xvalue) - ux_max = np.copy(u.xvalue) - uy_min = np.copy(u.yvalue) - uy_max = np.copy(u.yvalue) - uz_min = np.copy(u.zvalue) - uz_max = np.copy(u.zvalue) - ux_min[u_upwind.xvalue > 0.0] = 0.0 - ux_max[u_upwind.xvalue < 0.0] = 0.0 - uy_min[u_upwind.yvalue > 0.0] = 0.0 - uy_max[u_upwind.yvalue < 0.0] = 0.0 - uz_min[u_upwind.zvalue > 0.0] = 0.0 - uz_max[u_upwind.zvalue < 0.0] = 0.0 + ux_min = np.copy(u._xvalue) + ux_max = np.copy(u._xvalue) + uy_min = np.copy(u._yvalue) + uy_max = np.copy(u._yvalue) + uz_min = np.copy(u._zvalue) + uz_max = np.copy(u._zvalue) + ux_min[u_upwind._xvalue > 0.0] = 0.0 + ux_max[u_upwind._xvalue < 0.0] = 0.0 + uy_min[u_upwind._yvalue > 0.0] = 0.0 + uy_max[u_upwind._yvalue < 0.0] = 0.0 + uz_min[u_upwind._zvalue > 0.0] = 0.0 + uz_max[u_upwind._zvalue < 0.0] = 0.0 return ux_min, ux_max, uy_min, uy_max, uz_min, uz_max @@ -58,13 +58,13 @@ def convectionTerm1D(u: FaceVariable): # extract data from the mesh structure Nx = u.domain.dims[0] G = u.domain.cell_numbers() - #DX = u.domain.cellsize.x - DXe = u.domain.cellsize.x[2:] - DXw = u.domain.cellsize.x[0:-2] - DXp = u.domain.cellsize.x[1:-1] + #DX = u.domain.cellsize._x + DXe = u.domain.cellsize._x[2:] + DXw = u.domain.cellsize._x[0:-2] + DXp = u.domain.cellsize._x[1:-1] # reassign the east, west for code readability - ue = u.xvalue[1:Nx+1]/(DXp+DXe) - uw = u.xvalue[0:Nx]/(DXp+DXw) + ue = u._xvalue[1:Nx+1]/(DXp+DXe) + uw = u._xvalue[0:Nx]/(DXp+DXw) # build the sparse matrix based on the numbering system iix = np.tile(G[1:Nx+1].ravel(), 3) jjx = np.hstack([G[0:Nx], G[1:Nx+1], G[2:Nx+2]]) @@ -84,7 +84,7 @@ def convectionUpwindTerm1D(u: FaceVariable, *args): # extract data from the mesh structure Nx = u.domain.dims[0] G = u.domain.cell_numbers() - DXp = u.domain.cellsize.x[1:-1] + DXp = u.domain.cellsize._x[1:-1] # find the velocity direction for the upwind scheme ux_min, ux_max = _upwind_min_max(u, u_upwind) ue_min = ux_min[1:Nx+1] @@ -123,8 +123,8 @@ def convectionTvdRHS1D(u: FaceVariable, phi: CellVariable, u_upwind = u # extract data from the mesh structure Nx = u.domain.dims[0] - DXp = u.domain.cellsize.x[1:-1] - dx = 0.5*(u.domain.cellsize.x[0:-1]+u.domain.cellsize.x[1:]) + DXp = u.domain.cellsize._x[1:-1] + dx = 0.5*(u.domain.cellsize._x[0:-1]+u.domain.cellsize._x[1:]) RHS = np.zeros(Nx+2) psi_p = np.zeros(Nx+1) psi_m = np.zeros(Nx+1) @@ -158,15 +158,15 @@ def convectionTermCylindrical1D(u: FaceVariable): # extract data from the mesh structure Nx = u.domain.dims[0] G = u.domain.cell_numbers() - #DX = u.domain.cellsize.x - DXe = u.domain.cellsize.x[2:] - DXw = u.domain.cellsize.x[0:-2] - DXp = u.domain.cellsize.x[1:-1] - rp = u.domain.cellcenters.x - rf = u.domain.facecenters.x + #DX = u.domain.cellsize._x + DXe = u.domain.cellsize._x[2:] + DXw = u.domain.cellsize._x[0:-2] + DXp = u.domain.cellsize._x[1:-1] + rp = u.domain.cellcenters._x + rf = u.domain.facecenters._x # reassign the east, west for code readability - ue = rf[1:Nx+1]*u.xvalue[1:Nx+1]/(rp*(DXp+DXe)) - uw = rf[0:Nx]*u.xvalue[0:Nx]/(rp*(DXp+DXw)) + ue = rf[1:Nx+1]*u._xvalue[1:Nx+1]/(rp*(DXp+DXe)) + uw = rf[0:Nx]*u._xvalue[0:Nx]/(rp*(DXp+DXw)) # calculate the coefficients for the internal cells AE = ue AW = -uw @@ -190,9 +190,9 @@ def convectionUpwindTermCylindrical1D(u: FaceVariable, *args): # extract data from the mesh structure Nx = u.domain.dims[0] G = u.domain.cell_numbers() - DXp = u.domain.cellsize.x[1:-1] - rp = u.domain.cellcenters.x - rf = u.domain.facecenters.x + DXp = u.domain.cellsize._x[1:-1] + rp = u.domain.cellcenters._x + rf = u.domain.facecenters._x # find the velocity direction for the upwind scheme ux_min, ux_max = _upwind_min_max(u, u_upwind) ue_min = ux_min[1:Nx+1] @@ -231,10 +231,10 @@ def convectionTvdRHSCylindrical1D(u: FaceVariable, phi: CellVariable, u_upwind = u # extract data from the mesh structure Nx = u.domain.dims[0] - DXp = u.domain.cellsize.x[1:-1] - r = u.domain.cellcenters.x - rf = u.domain.facecenters.x - dx = 0.5*(u.domain.cellsize.x[0:-1]+u.domain.cellsize.x[1:]) + DXp = u.domain.cellsize._x[1:-1] + r = u.domain.cellcenters._x + rf = u.domain.facecenters._x + dx = 0.5*(u.domain.cellsize._x[0:-1]+u.domain.cellsize._x[1:]) RHS = np.zeros(Nx+2) psi_p = np.zeros(Nx+1) psi_m = np.zeros(Nx+1) @@ -268,19 +268,19 @@ def convectionTerm2D(u: FaceVariable): # extract data from the mesh structure Nx, Ny = u.domain.dims G = u.domain.cell_numbers() - DXe = u.domain.cellsize.x[2:][:, np.newaxis] - DXw = u.domain.cellsize.x[0:-2][:, np.newaxis] - DXp = u.domain.cellsize.x[1:-1][:, np.newaxis] - DYn = u.domain.cellsize.y[2:] - DYs = u.domain.cellsize.y[0:-2] - DYp = u.domain.cellsize.y[1:-1] + DXe = u.domain.cellsize._x[2:][:, np.newaxis] + DXw = u.domain.cellsize._x[0:-2][:, np.newaxis] + DXp = u.domain.cellsize._x[1:-1][:, np.newaxis] + DYn = u.domain.cellsize._y[2:] + DYs = u.domain.cellsize._y[0:-2] + DYp = u.domain.cellsize._y[1:-1] # define the vectors to store the sparse matrix data mn = Nx*Ny # reassign the east, west for code readability - ue = u.xvalue[1:Nx+1, :]/(DXp+DXe) - uw = u.xvalue[0:Nx, :]/(DXp+DXw) - vn = u.yvalue[:, 1:Ny+1]/(DYp+DYn) - vs = u.yvalue[:, 0:Ny]/(DYp+DYs) + ue = u._xvalue[1:Nx+1, :]/(DXp+DXe) + uw = u._xvalue[0:Nx, :]/(DXp+DXw) + vn = u._yvalue[:, 1:Ny+1]/(DYp+DYn) + vs = u._yvalue[:, 0:Ny]/(DYp+DYs) # calculate the coefficients for the internal cells AE = ue.ravel() AW = -uw.ravel() @@ -319,8 +319,8 @@ def convectionUpwindTerm2D(u: FaceVariable, *args): u_upwind = u Nx, Ny = u.domain.dims G = u.domain.cell_numbers() - DXp = u.domain.cellsize.x[1:-1][:, np.newaxis] - DYp = u.domain.cellsize.y[1:-1] + DXp = u.domain.cellsize._x[1:-1][:, np.newaxis] + DYp = u.domain.cellsize._y[1:-1] # define the vectors to store the sparse matrix data mn = Nx*Ny # find the velocity direction for the upwind scheme @@ -380,10 +380,10 @@ def convectionTvdRHS2D(u: FaceVariable, phi: CellVariable, FL, *args): # extract data from the mesh structure Nx, Ny = u.domain.dims G = u.domain.cell_numbers() - DXp = u.domain.cellsize.x[1:-1][:, np.newaxis] - DYp = u.domain.cellsize.y[1:-1][np.newaxis, :] - dx = 0.5*(u.domain.cellsize.x[0:-1]+u.domain.cellsize.x[1:])[:, np.newaxis] - dy = 0.5*(u.domain.cellsize.y[0:-1]+u.domain.cellsize.y[1:])[np.newaxis, :] + DXp = u.domain.cellsize._x[1:-1][:, np.newaxis] + DYp = u.domain.cellsize._y[1:-1][np.newaxis, :] + dx = 0.5*(u.domain.cellsize._x[0:-1]+u.domain.cellsize._x[1:])[:, np.newaxis] + dy = 0.5*(u.domain.cellsize._y[0:-1]+u.domain.cellsize._y[1:])[np.newaxis, :] psiX_p = np.zeros((Nx+1, Ny)) psiX_m = np.zeros((Nx+1, Ny)) psiY_p = np.zeros((Nx, Ny+1)) @@ -448,21 +448,21 @@ def convectionTermCylindrical2D(u: FaceVariable): # extract data from the mesh structure Nx, Ny = u.domain.dims G = u.domain.cell_numbers() - DXe = u.domain.cellsize.x[2:][:, np.newaxis] - DXw = u.domain.cellsize.x[0:-2][:, np.newaxis] - DXp = u.domain.cellsize.x[1:-1][:, np.newaxis] - DYn = u.domain.cellsize.y[2:] - DYs = u.domain.cellsize.y[0:-2] - DYp = u.domain.cellsize.y[1:-1] - rp = u.domain.cellcenters.x[:, np.newaxis] - rf = u.domain.facecenters.x[:, np.newaxis] + DXe = u.domain.cellsize._x[2:][:, np.newaxis] + DXw = u.domain.cellsize._x[0:-2][:, np.newaxis] + DXp = u.domain.cellsize._x[1:-1][:, np.newaxis] + DYn = u.domain.cellsize._y[2:] + DYs = u.domain.cellsize._y[0:-2] + DYp = u.domain.cellsize._y[1:-1] + rp = u.domain.cellcenters._x[:, np.newaxis] + rf = u.domain.facecenters._x[:, np.newaxis] # define the vectors to store the sparse matrix data mn = Nx*Ny # reassign the east, west for code readability - ue = rf[1:Nx+1, :]*u.xvalue[1:Nx+1, :]/(rp*(DXp+DXe)) - uw = rf[0:Nx, :]*u.xvalue[0:Nx, :]/(rp*(DXp+DXe)) - vn = u.yvalue[:, 1:Ny+1]/(DYp+DYn) - vs = u.yvalue[:, 0:Ny]/(DYp+DYs) + ue = rf[1:Nx+1, :]*u._xvalue[1:Nx+1, :]/(rp*(DXp+DXe)) + uw = rf[0:Nx, :]*u._xvalue[0:Nx, :]/(rp*(DXp+DXe)) + vn = u._yvalue[:, 1:Ny+1]/(DYp+DYn) + vs = u._yvalue[:, 0:Ny]/(DYp+DYs) # calculate the coefficients for the internal cells AE = ue.ravel() AW = -uw.ravel() @@ -501,10 +501,10 @@ def convectionUpwindTermCylindrical2D(u: FaceVariable, *args): u_upwind = u Nx, Ny = u.domain.dims G = u.domain.cell_numbers() - DXp = u.domain.cellsize.x[1:-1][:, np.newaxis] - DYp = u.domain.cellsize.y[1:-1] - rp = u.domain.cellcenters.x[:, np.newaxis] - rf = u.domain.facecenters.x[:, np.newaxis] + DXp = u.domain.cellsize._x[1:-1][:, np.newaxis] + DYp = u.domain.cellsize._y[1:-1] + rp = u.domain.cellcenters._x[:, np.newaxis] + rf = u.domain.facecenters._x[:, np.newaxis] re = rf[1:Nx+1, :] rw = rf[0:Nx, :] # define the vectors to store the sparse matrix data @@ -566,14 +566,14 @@ def convectionTvdRHSCylindrical2D(u: FaceVariable, phi: CellVariable, FL, *args) # extract data from the mesh structure Nx, Ny = u.domain.dims G = u.domain.cell_numbers() - DXp = u.domain.cellsize.x[1:-1][:, np.newaxis] - DYp = u.domain.cellsize.y[1:-1][np.newaxis, :] - rp = u.domain.cellcenters.x[:, np.newaxis] - rf = u.domain.facecenters.x[:, np.newaxis] + DXp = u.domain.cellsize._x[1:-1][:, np.newaxis] + DYp = u.domain.cellsize._y[1:-1][np.newaxis, :] + rp = u.domain.cellcenters._x[:, np.newaxis] + rf = u.domain.facecenters._x[:, np.newaxis] re = rf[1:Nx+1, :] rw = rf[0:Nx, :] - dx = 0.5*(u.domain.cellsize.x[0:-1]+u.domain.cellsize.x[1:])[:, np.newaxis] - dy = 0.5*(u.domain.cellsize.y[0:-1]+u.domain.cellsize.y[1:])[np.newaxis, :] + dx = 0.5*(u.domain.cellsize._x[0:-1]+u.domain.cellsize._x[1:])[:, np.newaxis] + dy = 0.5*(u.domain.cellsize._y[0:-1]+u.domain.cellsize._y[1:])[np.newaxis, :] psiX_p = np.zeros((Nx+1, Ny)) psiX_m = np.zeros((Nx+1, Ny)) psiY_p = np.zeros((Nx, Ny+1)) @@ -630,21 +630,21 @@ def convectionTermPolar2D(u: FaceVariable): # extract data from the mesh structure Nx, Ny = u.domain.dims G = u.domain.cell_numbers() - DXe = u.domain.cellsize.x[2:][:, np.newaxis] - DXw = u.domain.cellsize.x[0:-2][:, np.newaxis] - DXp = u.domain.cellsize.x[1:-1][:, np.newaxis] - DYn = u.domain.cellsize.y[2:] - DYs = u.domain.cellsize.y[0:-2] - DYp = u.domain.cellsize.y[1:-1] - rp = u.domain.cellcenters.x[:, np.newaxis] - rf = u.domain.facecenters.x[:, np.newaxis] + DXe = u.domain.cellsize._x[2:][:, np.newaxis] + DXw = u.domain.cellsize._x[0:-2][:, np.newaxis] + DXp = u.domain.cellsize._x[1:-1][:, np.newaxis] + DYn = u.domain.cellsize._y[2:] + DYs = u.domain.cellsize._y[0:-2] + DYp = u.domain.cellsize._y[1:-1] + rp = u.domain.cellcenters._x[:, np.newaxis] + rf = u.domain.facecenters._x[:, np.newaxis] # define the vectors to store the sparse matrix data mn = Nx*Ny # reassign the east, west for code readability - ue = rf[1:Nx+1, :]*u.xvalue[1:Nx+1, :]/(rp*(DXp+DXe)) - uw = rf[0:Nx, :]*u.xvalue[0:Nx, :]/(rp*(DXp+DXe)) - vn = u.yvalue[:, 1:Ny+1]/(rp*(DYp+DYn)) - vs = u.yvalue[:, 0:Ny]/(rp*(DYp+DYs)) + ue = rf[1:Nx+1, :]*u._xvalue[1:Nx+1, :]/(rp*(DXp+DXe)) + uw = rf[0:Nx, :]*u._xvalue[0:Nx, :]/(rp*(DXp+DXe)) + vn = u._yvalue[:, 1:Ny+1]/(rp*(DYp+DYn)) + vs = u._yvalue[:, 0:Ny]/(rp*(DYp+DYs)) # calculate the coefficients for the internal cells AE = ue.ravel() AW = -uw.ravel() @@ -683,10 +683,10 @@ def convectionUpwindTermPolar2D(u: FaceVariable, *args): u_upwind = u Nx, Ny = u.domain.dims G = u.domain.cell_numbers() - DXp = u.domain.cellsize.x[1:-1][:, np.newaxis] - DYp = u.domain.cellsize.y[1:-1] - rp = u.domain.cellcenters.x[:, np.newaxis] - rf = u.domain.facecenters.x[:, np.newaxis] + DXp = u.domain.cellsize._x[1:-1][:, np.newaxis] + DYp = u.domain.cellsize._y[1:-1] + rp = u.domain.cellcenters._x[:, np.newaxis] + rf = u.domain.facecenters._x[:, np.newaxis] re = rf[1:Nx+1, :] rw = rf[0:Nx, :] # define the vectors to store the sparse matrix data @@ -748,14 +748,14 @@ def convectionTvdRHSPolar2D(u: FaceVariable, phi: CellVariable, FL, *args): # extract data from the mesh structure Nx, Ny = u.domain.dims G = u.domain.cell_numbers() - DXp = u.domain.cellsize.x[1:-1][:, np.newaxis] - DYp = u.domain.cellsize.y[1:-1][np.newaxis, :] - rp = u.domain.cellcenters.x[:, np.newaxis] - rf = u.domain.facecenters.x[:, np.newaxis] + DXp = u.domain.cellsize._x[1:-1][:, np.newaxis] + DYp = u.domain.cellsize._y[1:-1][np.newaxis, :] + rp = u.domain.cellcenters._x[:, np.newaxis] + rf = u.domain.facecenters._x[:, np.newaxis] re = rf[1:Nx+1, :] rw = rf[0:Nx, :] - dx = 0.5*(u.domain.cellsize.x[0:-1]+u.domain.cellsize.x[1:])[:, np.newaxis] - dy = 0.5*(u.domain.cellsize.y[0:-1]+u.domain.cellsize.y[1:])[np.newaxis, :] + dx = 0.5*(u.domain.cellsize._x[0:-1]+u.domain.cellsize._x[1:])[:, np.newaxis] + dy = 0.5*(u.domain.cellsize._y[0:-1]+u.domain.cellsize._y[1:])[np.newaxis, :] psiX_p = np.zeros((Nx+1, Ny)) psiX_m = np.zeros((Nx+1, Ny)) psiY_p = np.zeros((Nx, Ny+1)) @@ -812,25 +812,25 @@ def convectionTerm3D(u: FaceVariable): # extract data from the mesh structure Nx, Ny, Nz = u.domain.dims G = u.domain.cell_numbers() - DXe = u.domain.cellsize.x[2:][:, np.newaxis, np.newaxis] - DXw = u.domain.cellsize.x[0:-2][:, np.newaxis, np.newaxis] - DXp = u.domain.cellsize.x[1:-1][:, np.newaxis, np.newaxis] - DYn = u.domain.cellsize.y[2:][np.newaxis, :, np.newaxis] - DYs = u.domain.cellsize.y[0:-2][np.newaxis, :, np.newaxis] - DYp = u.domain.cellsize.y[1:-1][np.newaxis, :, np.newaxis] - DZf = u.domain.cellsize.z[2:][np.newaxis, np.newaxis, :] - DZb = u.domain.cellsize.z[0:-2][np.newaxis, np.newaxis, :] - DZp = u.domain.cellsize.z[1:-1][np.newaxis, np.newaxis, :] + DXe = u.domain.cellsize._x[2:][:, np.newaxis, np.newaxis] + DXw = u.domain.cellsize._x[0:-2][:, np.newaxis, np.newaxis] + DXp = u.domain.cellsize._x[1:-1][:, np.newaxis, np.newaxis] + DYn = u.domain.cellsize._y[2:][np.newaxis, :, np.newaxis] + DYs = u.domain.cellsize._y[0:-2][np.newaxis, :, np.newaxis] + DYp = u.domain.cellsize._y[1:-1][np.newaxis, :, np.newaxis] + DZf = u.domain.cellsize._z[2:][np.newaxis, np.newaxis, :] + DZb = u.domain.cellsize._z[0:-2][np.newaxis, np.newaxis, :] + DZp = u.domain.cellsize._z[1:-1][np.newaxis, np.newaxis, :] # define the vectors to stores the sparse matrix data mn = Nx*Ny*Nz # reassign the east, west, north, and south velocity vectors for the # code readability - ue = u.xvalue[1:Nx+1, :, :]/(DXp+DXe) - uw = u.xvalue[0:Nx, :, :]/(DXp+DXw) - vn = u.yvalue[:, 1:Ny+1, :]/(DYp+DYn) - vs = u.yvalue[:, 0:Ny, :]/(DYp+DYs) - wf = u.zvalue[:, :, 1:Nz+1]/(DZp+DZf) - wb = u.zvalue[:, :, 0:Nz]/(DZp+DZb) + ue = u._xvalue[1:Nx+1, :, :]/(DXp+DXe) + uw = u._xvalue[0:Nx, :, :]/(DXp+DXw) + vn = u._yvalue[:, 1:Ny+1, :]/(DYp+DYn) + vs = u._yvalue[:, 0:Ny, :]/(DYp+DYs) + wf = u._zvalue[:, :, 1:Nz+1]/(DZp+DZf) + wb = u._zvalue[:, :, 0:Nz]/(DZp+DZb) # calculate the coefficients for the internal cells AE = ue.ravel() AW = -uw.ravel() @@ -878,9 +878,9 @@ def convectionUpwindTerm3D(u: FaceVariable, *args): u_upwind = u Nx, Ny, Nz = u.domain.dims G = u.domain.cell_numbers() - DXp = u.domain.cellsize.x[1:-1][:, np.newaxis, np.newaxis] - DYp = u.domain.cellsize.y[1:-1][np.newaxis, :, np.newaxis] - DZp = u.domain.cellsize.z[1:-1][np.newaxis, np.newaxis, :] + DXp = u.domain.cellsize._x[1:-1][:, np.newaxis, np.newaxis] + DYp = u.domain.cellsize._y[1:-1][np.newaxis, :, np.newaxis] + DZp = u.domain.cellsize._z[1:-1][np.newaxis, np.newaxis, :] # define the vectors to stores the sparse matrix data mn = Nx*Ny*Nz @@ -973,16 +973,16 @@ def convectionTvdRHS3D(u: FaceVariable, phi: CellVariable, FL, *args): # extract data from the mesh structure Nx, Ny, Nz = u.domain.dims G = u.domain.cell_numbers() - DXp = u.domain.cellsize.x[1:-1][:, np.newaxis, np.newaxis] - DYp = u.domain.cellsize.y[1:-1][np.newaxis, :, np.newaxis] - DZp = u.domain.cellsize.z[1:-1][np.newaxis, np.newaxis, :] + DXp = u.domain.cellsize._x[1:-1][:, np.newaxis, np.newaxis] + DYp = u.domain.cellsize._y[1:-1][np.newaxis, :, np.newaxis] + DZp = u.domain.cellsize._z[1:-1][np.newaxis, np.newaxis, :] # define the vectors to stores the sparse matrix data - dx = 0.5*(u.domain.cellsize.x[0:-1] + - u.domain.cellsize.x[1:])[:, np.newaxis, np.newaxis] - dy = 0.5*(u.domain.cellsize.y[0:-1] + - u.domain.cellsize.y[1:])[np.newaxis, :, np.newaxis] - dz = 0.5*(u.domain.cellsize.z[0:-1] + - u.domain.cellsize.z[1:])[np.newaxis, np.newaxis, :] + dx = 0.5*(u.domain.cellsize._x[0:-1] + + u.domain.cellsize._x[1:])[:, np.newaxis, np.newaxis] + dy = 0.5*(u.domain.cellsize._y[0:-1] + + u.domain.cellsize._y[1:])[np.newaxis, :, np.newaxis] + dz = 0.5*(u.domain.cellsize._z[0:-1] + + u.domain.cellsize._z[1:])[np.newaxis, np.newaxis, :] psiX_p = np.zeros((Nx+1, Ny, Nz)) psiX_m = np.zeros((Nx+1, Ny, Nz)) psiY_p = np.zeros((Nx, Ny+1, Nz)) @@ -1065,27 +1065,27 @@ def convectionTermCylindrical3D(u: FaceVariable): # extract data from the mesh structure Nr, Ntheta, Nz = u.domain.dims G = u.domain.cell_numbers() - DRe = u.domain.cellsize.x[2:][:, np.newaxis, np.newaxis] - DRw = u.domain.cellsize.x[0:-2][:, np.newaxis, np.newaxis] - DRp = u.domain.cellsize.x[1:-1][:, np.newaxis, np.newaxis] - DTHETAn = u.domain.cellsize.y[2:][np.newaxis, :, np.newaxis] - DTHETAs = u.domain.cellsize.y[0:-2][np.newaxis, :, np.newaxis] - DTHETAp = u.domain.cellsize.y[1:-1][np.newaxis, :, np.newaxis] - DZf = u.domain.cellsize.z[2:][np.newaxis, np.newaxis, :] - DZb = u.domain.cellsize.z[0:-2][np.newaxis, np.newaxis, :] - DZp = u.domain.cellsize.z[1:-1][np.newaxis, np.newaxis, :] - rp = u.domain.cellcenters.x[:, np.newaxis, np.newaxis] - rf = u.domain.facecenters.x[:, np.newaxis, np.newaxis] + DRe = u.domain.cellsize._x[2:][:, np.newaxis, np.newaxis] + DRw = u.domain.cellsize._x[0:-2][:, np.newaxis, np.newaxis] + DRp = u.domain.cellsize._x[1:-1][:, np.newaxis, np.newaxis] + DTHETAn = u.domain.cellsize._y[2:][np.newaxis, :, np.newaxis] + DTHETAs = u.domain.cellsize._y[0:-2][np.newaxis, :, np.newaxis] + DTHETAp = u.domain.cellsize._y[1:-1][np.newaxis, :, np.newaxis] + DZf = u.domain.cellsize._z[2:][np.newaxis, np.newaxis, :] + DZb = u.domain.cellsize._z[0:-2][np.newaxis, np.newaxis, :] + DZp = u.domain.cellsize._z[1:-1][np.newaxis, np.newaxis, :] + rp = u.domain.cellcenters._x[:, np.newaxis, np.newaxis] + rf = u.domain.facecenters._x[:, np.newaxis, np.newaxis] # define the vectors to stores the sparse matrix data mn = Nr*Ntheta*Nz # reassign the east, west, north, and south velocity vectors for the # code readability - ue = rf[1:Nr+1]*u.xvalue[1:Nr+1, :, :]/(rp*(DRp+DRe)) - uw = rf[0:Nr]*u.xvalue[0:Nr, :, :]/(rp*(DRp+DRw)) - vn = u.yvalue[:, 1:Ntheta+1, :]/(rp*(DTHETAp+DTHETAn)) - vs = u.yvalue[:, 0:Ntheta, :]/(rp*(DTHETAp+DTHETAs)) - wf = u.zvalue[:, :, 1:Nz+1]/(DZp+DZf) - wb = u.zvalue[:, :, 0:Nz]/(DZp+DZb) + ue = rf[1:Nr+1]*u._xvalue[1:Nr+1, :, :]/(rp*(DRp+DRe)) + uw = rf[0:Nr]*u._xvalue[0:Nr, :, :]/(rp*(DRp+DRw)) + vn = u._yvalue[:, 1:Ntheta+1, :]/(rp*(DTHETAp+DTHETAn)) + vs = u._yvalue[:, 0:Ntheta, :]/(rp*(DTHETAp+DTHETAs)) + wf = u._zvalue[:, :, 1:Nz+1]/(DZp+DZf) + wb = u._zvalue[:, :, 0:Nz]/(DZp+DZb) # calculate the coefficients for the internal cells AE = ue.ravel() @@ -1134,17 +1134,17 @@ def convectionUpwindTermCylindrical3D(u: FaceVariable, *args): u_upwind = u Nr, Ntheta, Nz = u.domain.dims G = u.domain.cell_numbers() - DRe = u.domain.cellsize.x[2:][:, np.newaxis, np.newaxis] - DRw = u.domain.cellsize.x[0:-2][:, np.newaxis, np.newaxis] - DRp = u.domain.cellsize.x[1:-1][:, np.newaxis, np.newaxis] - DTHETAn = u.domain.cellsize.y[2:][np.newaxis, :, np.newaxis] - DTHETAs = u.domain.cellsize.y[0:-2][np.newaxis, :, np.newaxis] - DTHETAp = u.domain.cellsize.y[1:-1][np.newaxis, :, np.newaxis] - DZf = u.domain.cellsize.z[2:][np.newaxis, np.newaxis, :] - DZb = u.domain.cellsize.z[0:-2][np.newaxis, np.newaxis, :] - DZp = u.domain.cellsize.z[1:-1][np.newaxis, np.newaxis, :] - rp = u.domain.cellcenters.x[:, np.newaxis, np.newaxis] - rf = u.domain.facecenters.x[:, np.newaxis, np.newaxis] + DRe = u.domain.cellsize._x[2:][:, np.newaxis, np.newaxis] + DRw = u.domain.cellsize._x[0:-2][:, np.newaxis, np.newaxis] + DRp = u.domain.cellsize._x[1:-1][:, np.newaxis, np.newaxis] + DTHETAn = u.domain.cellsize._y[2:][np.newaxis, :, np.newaxis] + DTHETAs = u.domain.cellsize._y[0:-2][np.newaxis, :, np.newaxis] + DTHETAp = u.domain.cellsize._y[1:-1][np.newaxis, :, np.newaxis] + DZf = u.domain.cellsize._z[2:][np.newaxis, np.newaxis, :] + DZb = u.domain.cellsize._z[0:-2][np.newaxis, np.newaxis, :] + DZp = u.domain.cellsize._z[1:-1][np.newaxis, np.newaxis, :] + rp = u.domain.cellcenters._x[:, np.newaxis, np.newaxis] + rf = u.domain.facecenters._x[:, np.newaxis, np.newaxis] mn = Nr*Ntheta*Nz re = rf[1:Nr+1, :, :] rw = rf[0:Nr, :, :] @@ -1227,30 +1227,30 @@ def convectionTvdRHSCylindrical3D(u: FaceVariable, phi: CellVariable, FL, *args) u_upwind = u Nr, Ntheta, Nz = u.domain.dims G = u.domain.cell_numbers() - DRe = u.domain.cellsize.x[2:][:, np.newaxis, np.newaxis] - DRw = u.domain.cellsize.x[0:-2][:, np.newaxis, np.newaxis] - DRp = u.domain.cellsize.x[1:-1][:, np.newaxis, np.newaxis] - DTHETAn = u.domain.cellsize.y[2:][np.newaxis, :, np.newaxis] - DTHETAs = u.domain.cellsize.y[0:-2][np.newaxis, :, np.newaxis] - DTHETAp = u.domain.cellsize.y[1:-1][np.newaxis, :, np.newaxis] - DZf = u.domain.cellsize.z[2:][np.newaxis, np.newaxis, :] - DZb = u.domain.cellsize.z[0:-2][np.newaxis, np.newaxis, :] - DZp = u.domain.cellsize.z[1:-1][np.newaxis, np.newaxis, :] - dr = 0.5*(u.domain.cellsize.x[0:-1] + - u.domain.cellsize.x[1:])[:, np.newaxis, np.newaxis] + DRe = u.domain.cellsize._x[2:][:, np.newaxis, np.newaxis] + DRw = u.domain.cellsize._x[0:-2][:, np.newaxis, np.newaxis] + DRp = u.domain.cellsize._x[1:-1][:, np.newaxis, np.newaxis] + DTHETAn = u.domain.cellsize._y[2:][np.newaxis, :, np.newaxis] + DTHETAs = u.domain.cellsize._y[0:-2][np.newaxis, :, np.newaxis] + DTHETAp = u.domain.cellsize._y[1:-1][np.newaxis, :, np.newaxis] + DZf = u.domain.cellsize._z[2:][np.newaxis, np.newaxis, :] + DZb = u.domain.cellsize._z[0:-2][np.newaxis, np.newaxis, :] + DZp = u.domain.cellsize._z[1:-1][np.newaxis, np.newaxis, :] + dr = 0.5*(u.domain.cellsize._x[0:-1] + + u.domain.cellsize._x[1:])[:, np.newaxis, np.newaxis] dtheta = 0.5 * \ - (u.domain.cellsize.y[0:-1]+u.domain.cellsize.y[1:] + (u.domain.cellsize._y[0:-1]+u.domain.cellsize._y[1:] )[np.newaxis, :, np.newaxis] - dz = 0.5*(u.domain.cellsize.z[0:-1] + - u.domain.cellsize.z[1:])[np.newaxis, np.newaxis, :] + dz = 0.5*(u.domain.cellsize._z[0:-1] + + u.domain.cellsize._z[1:])[np.newaxis, np.newaxis, :] psiX_p = np.zeros((Nr+1, Ntheta, Nz)) psiX_m = np.zeros((Nr+1, Ntheta, Nz)) psiY_p = np.zeros((Nr, Ntheta+1, Nz)) psiY_m = np.zeros((Nr, Ntheta+1, Nz)) psiZ_p = np.zeros((Nr, Ntheta, Nz+1)) psiZ_m = np.zeros((Nr, Ntheta, Nz+1)) - rp = u.domain.cellcenters.x[:, np.newaxis, np.newaxis] - rf = u.domain.facecenters.x[:, np.newaxis, np.newaxis] + rp = u.domain.cellcenters._x[:, np.newaxis, np.newaxis] + rf = u.domain.facecenters._x[:, np.newaxis, np.newaxis] # calculate the upstream to downstream gradient ratios for u>0 (+ ratio) # x direction diff --git a/src/pyfvtool/averaging.py b/src/pyfvtool/averaging.py index b5f3e46..cc474fe 100644 --- a/src/pyfvtool/averaging.py +++ b/src/pyfvtool/averaging.py @@ -25,16 +25,16 @@ def cell_size_array(m: MeshStructure): if issubclass(type(m), Grid1D): - dx = m.cellsize.x + dx = m.cellsize._x return dx elif issubclass(type(m), Grid2D): - dx = m.cellsize.x[:,np.newaxis] - dy = m.cellsize.y[np.newaxis,:] + dx = m.cellsize._x[:,np.newaxis] + dy = m.cellsize._y[np.newaxis,:] return dx, dy elif issubclass(type(m), Grid3D): - dx = m.cellsize.x[:,np.newaxis,np.newaxis] - dy = m.cellsize.y[np.newaxis,:,np.newaxis] - dz = m.cellsize.z[np.newaxis,np.newaxis,:] + dx = m.cellsize._x[:,np.newaxis,np.newaxis] + dy = m.cellsize._y[np.newaxis,:,np.newaxis] + dz = m.cellsize._z[np.newaxis,np.newaxis,:] return dx, dy, dz @@ -335,7 +335,7 @@ def upwindMean(phi: CellVariable, u: FaceVariable): # currently, it assumes a uniform mesh. phi_tmp = np.copy(phi.value) if issubclass(type(phi.domain), Grid1D): - ux = u.xvalue + ux = u._xvalue # assign the value of the left boundary to the left ghost cell phi_tmp[0] = 0.5*(phi.value[0]+phi.value[1]) # assign the value of the right boundary to the right ghost cell @@ -346,8 +346,8 @@ def upwindMean(phi: CellVariable, u: FaceVariable): np.array([]), np.array([])) elif issubclass(type(phi.domain), Grid2D): - ux = u.xvalue - uy = u.yvalue + ux = u._xvalue + uy = u._yvalue # assign the value of the left boundary to the left ghost cells phi_tmp[0,:] = 0.5*(phi.value[0,:]+phi.value[1,:]) # assign the value of the right boundary to the right ghost cells @@ -365,9 +365,9 @@ def upwindMean(phi: CellVariable, u: FaceVariable): 0.5*(uy==0.0)*(phi.value[1:-1,0:-1]+phi.value[1:-1,1:]), np.array([])) elif issubclass(type(phi.domain), Grid3D): - ux = u.xvalue - uy = u.yvalue - uz = u.zvalue + ux = u._xvalue + uy = u._yvalue + uz = u._zvalue # assign the value of the left boundary to the left ghost cells phi_tmp[0,:,:] = 0.5*(phi.value[0,:,:]+phi.value[1,:,:]) # assign the value of the right boundary to the right ghost cells @@ -438,12 +438,12 @@ def tvdMean(phi: CellVariable, u: FaceVariable, FL): # if issubclass(type(phi.domain), Mesh1D): # # extract data from the mesh structure # Nx = u.domain.dims[0] -# dx = 0.5*(u.domain.cellsize.x[0:-1]+u.domain.cellsize.x[1:]) +# dx = 0.5*(u.domain.cellsize._x[0:-1]+u.domain.cellsize._x[1:]) # phi_p = np.zeros(Float64, Nx+1) # phi_m = np.zeros(Float64, Nx+1) # # extract the velocity data -# ux = u.xvalue +# ux = u._xvalue # # calculate the upstream to downstream gradient ratios for u>0 (+ ratio) # dphi_p = (phi.value[1:Nx+2]-phi.value[0:Nx+1])/dx @@ -465,9 +465,9 @@ def tvdMean(phi: CellVariable, u: FaceVariable, FL): # # extract data from the mesh structure # Nx = u.domain.dims[0] # Ny = u.domain.dims[2] -# dx=0.5*(u.domain.cellsize.x[0:-1]+u.domain.cellsize.x[1:]) +# dx=0.5*(u.domain.cellsize._x[0:-1]+u.domain.cellsize._x[1:]) # dy=np.zeros( 1, Ny+1) -# dy[:]=0.5*(u.domain.cellsize.y[0:-1]+u.domain.cellsize.y[1:]) +# dy[:]=0.5*(u.domain.cellsize._y[0:-1]+u.domain.cellsize._y[1:]) # phi_p = np.zeros(Float64, Nx+1) # phi_m = np.zeros(Float64, Nx+1) # phiX_p = np.zeros(Float64, Nx+1, Ny) @@ -476,8 +476,8 @@ def tvdMean(phi: CellVariable, u: FaceVariable, FL): # phiY_m = np.zeros(Float64, Nx,Ny+1) # # extract the velocity data -# ux = u.xvalue -# uy = u.yvalue +# ux = u._xvalue +# uy = u._yvalue # # calculate the upstream to downstream gradient ratios for u>0 (+ ratio) # # x direction @@ -517,15 +517,15 @@ def tvdMean(phi: CellVariable, u: FaceVariable, FL): # Nx = u.domain.dims[0] # Ny = u.domain.dims[2] # Nz = u.domain.dims[3] -# dx=0.5*(u.domain.cellsize.x[0:-1]+u.domain.cellsize.x[1:]) +# dx=0.5*(u.domain.cellsize._x[0:-1]+u.domain.cellsize._x[1:]) # dy=np.zeros( 1, Ny+1) -# dy[:]=0.5*(u.domain.cellsize.y[0:-1]+u.domain.cellsize.y[1:]) +# dy[:]=0.5*(u.domain.cellsize._y[0:-1]+u.domain.cellsize._y[1:]) # dz=np.zeros( 1, 1, Nz+1) -# dz[:]=0.5*(u.domain.cellsize.z[0:-1]+u.domain.cellsize.z[1:]) +# dz[:]=0.5*(u.domain.cellsize._z[0:-1]+u.domain.cellsize._z[1:]) # # extract the velocity data -# ux = u.xvalue -# uy = u.yvalue -# uz = u.zvalue +# ux = u._xvalue +# uy = u._yvalue +# uz = u._zvalue # # define the tvd face vectors # phiX_p = np.zeros(Float64, Nx+1,Ny,Nz) diff --git a/src/pyfvtool/boundary.py b/src/pyfvtool/boundary.py index b3a4b1d..fe6881a 100644 --- a/src/pyfvtool/boundary.py +++ b/src/pyfvtool/boundary.py @@ -178,8 +178,8 @@ def cellValuesWithBoundaries1D(phi, BC): # extract data from the mesh structure # Nx = MeshStructure.numberofcells - dx_1 = BC.domain.cellsize.x[1] - dx_end = BC.domain.cellsize.x[-1] + dx_1 = BC.domain.cellsize._x[1] + dx_end = BC.domain.cellsize._x[-1] # boundary condition (a d\\phi/dx + b \\phi = c, a column vector of [d a]) # a (phi(i)-phi(i-1))/dx + b (phi(i)+phi(i-1))/2 = c @@ -205,10 +205,10 @@ def cellValuesWithBoundaries2D(phi, BC): # extract data from the mesh structure Nx, Ny = BC.domain.dims - dx_1 = BC.domain.cellsize.x[0] - dx_end = BC.domain.cellsize.x[-1] - dy_1 = BC.domain.cellsize.y[0] - dy_end = BC.domain.cellsize.y[-1] + dx_1 = BC.domain.cellsize._x[0] + dx_end = BC.domain.cellsize._x[-1] + dy_1 = BC.domain.cellsize._y[0] + dy_end = BC.domain.cellsize._y[-1] # define the output matrix phiBC = np.zeros((Nx+2, Ny+2)) @@ -261,12 +261,12 @@ def cellValuesWithBoundaries3D(phi, BC): """ Nx, Ny, Nz = BC.domain.dims - dx_1 = BC.domain.cellsize.x[0] - dx_end = BC.domain.cellsize.x[-1] - dy_1 = BC.domain.cellsize.y[0] - dy_end = BC.domain.cellsize.y[-1] - dz_1 = BC.domain.cellsize.z[0] - dz_end = BC.domain.cellsize.z[-1] + dx_1 = BC.domain.cellsize._x[0] + dx_end = BC.domain.cellsize._x[-1] + dy_1 = BC.domain.cellsize._y[0] + dy_end = BC.domain.cellsize._y[-1] + dz_1 = BC.domain.cellsize._z[0] + dz_end = BC.domain.cellsize._z[-1] i_ind = int_range(1,Nx)[:, np.newaxis, np.newaxis] j_ind = int_range(1,Ny)[np.newaxis, :, np.newaxis] @@ -360,13 +360,13 @@ def cellValuesWithBoundariesCylindrical3D(phi, BC): """ Nx, Ny, Nz = BC.domain.dims - dx_1 = BC.domain.cellsize.x[0] - dx_end = BC.domain.cellsize.x[-1] - dy_1 = BC.domain.cellsize.y[0] - dy_end = BC.domain.cellsize.y[-1] - dz_1 = BC.domain.cellsize.z[0] - dz_end = BC.domain.cellsize.z[-1] - rp = BC.domain.cellcenters.x[:, np.newaxis] + dx_1 = BC.domain.cellsize._x[0] + dx_end = BC.domain.cellsize._x[-1] + dy_1 = BC.domain.cellsize._y[0] + dy_end = BC.domain.cellsize._y[-1] + dz_1 = BC.domain.cellsize._z[0] + dz_end = BC.domain.cellsize._z[-1] + rp = BC.domain.cellcenters._x[:, np.newaxis] i_ind = int_range(1,Nx)[:, np.newaxis, np.newaxis] j_ind = int_range(1,Ny)[np.newaxis, :, np.newaxis] @@ -461,11 +461,11 @@ def cellValuesWithBoundariesPolar2D(phi, BC): # extract data from the mesh structure Nx, Ny = BC.domain.dims - dx_1 = BC.domain.cellsize.x[0] - dx_end = BC.domain.cellsize.x[-1] - dy_1 = BC.domain.cellsize.y[0] - dy_end = BC.domain.cellsize.y[-1] - rp = BC.domain.cellcenters.x + dx_1 = BC.domain.cellsize._x[0] + dx_end = BC.domain.cellsize._x[-1] + dy_1 = BC.domain.cellsize._y[0] + dy_end = BC.domain.cellsize._y[-1] + rp = BC.domain.cellcenters._x # define the output matrix phiBC = np.zeros((Nx+2, Ny+2)) @@ -566,8 +566,8 @@ def cellValuesWithBoundaries(phi, BC) -> np.ndarray: def boundaryConditionsTerm1D(BC: BoundaryConditions1D): Nx = BC.domain.dims[0] - dx_1 = BC.domain.cellsize.x[0] - dx_end = BC.domain.cellsize.x[-1] + dx_1 = BC.domain.cellsize._x[0] + dx_end = BC.domain.cellsize._x[-1] G = int_range(1, Nx+2)-1 nb = 8 # number of boundary nodes @@ -648,10 +648,10 @@ def boundaryConditionsTerm1D(BC: BoundaryConditions1D): def boundaryConditionsTerm2D(BC: BoundaryConditions2D): Nx, Ny = BC.domain.dims - dx_1 = BC.domain.cellsize.x[0] - dx_end = BC.domain.cellsize.x[-1] - dy_1 = BC.domain.cellsize.y[0] - dy_end = BC.domain.cellsize.y[-1] + dx_1 = BC.domain.cellsize._x[0] + dx_end = BC.domain.cellsize._x[-1] + dy_1 = BC.domain.cellsize._y[0] + dy_end = BC.domain.cellsize._y[-1] G = BC.domain.cell_numbers() nb = 8*(Nx+Ny+2) # number of boundary nodes @@ -814,12 +814,12 @@ def boundaryConditionsTerm3D(BC: BoundaryConditions3D): # extract data from the mesh structure Nx, Ny, Nz = BC.domain.dims G=BC.domain.cell_numbers() - dx_1 = BC.domain.cellsize.x[0] - dx_end = BC.domain.cellsize.x[-1] - dy_1 = BC.domain.cellsize.y[0] - dy_end = BC.domain.cellsize.y[-1] - dz_1 = BC.domain.cellsize.z[0] - dz_end = BC.domain.cellsize.z[-1] + dx_1 = BC.domain.cellsize._x[0] + dx_end = BC.domain.cellsize._x[-1] + dy_1 = BC.domain.cellsize._y[0] + dy_end = BC.domain.cellsize._y[-1] + dz_1 = BC.domain.cellsize._z[0] + dz_end = BC.domain.cellsize._z[-1] i_ind = int_range(1,Nx)[:, np.newaxis, np.newaxis] j_ind = int_range(1,Ny)[np.newaxis, :, np.newaxis] @@ -1077,11 +1077,11 @@ def boundaryConditionsTerm3D(BC: BoundaryConditions3D): def boundaryConditionsTermPolar2D(BC: BoundaryConditions2D): Nx, Ny = BC.domain.dims - dx_1 = BC.domain.cellsize.x[0] - dx_end = BC.domain.cellsize.x[-1] - dy_1 = BC.domain.cellsize.y[0] - dy_end = BC.domain.cellsize.y[-1] - rp = BC.domain.cellcenters.x + dx_1 = BC.domain.cellsize._x[0] + dx_end = BC.domain.cellsize._x[-1] + dy_1 = BC.domain.cellsize._y[0] + dy_end = BC.domain.cellsize._y[-1] + rp = BC.domain.cellcenters._x G = BC.domain.cell_numbers() nb = 8*(Nx+Ny+2) # number of boundary nodes @@ -1244,13 +1244,13 @@ def boundaryConditionsTermCylindrical3D(BC: BoundaryConditions3D): # extract data from the mesh structure Nx, Ny, Nz = BC.domain.dims G=BC.domain.cell_numbers() - dx_1 = BC.domain.cellsize.x[0] - dx_end = BC.domain.cellsize.x[-1] - dy_1 = BC.domain.cellsize.y[0] - dy_end = BC.domain.cellsize.y[-1] - dz_1 = BC.domain.cellsize.z[0] - dz_end = BC.domain.cellsize.z[-1] - rp = BC.domain.cellcenters.x[:, np.newaxis] + dx_1 = BC.domain.cellsize._x[0] + dx_end = BC.domain.cellsize._x[-1] + dy_1 = BC.domain.cellsize._y[0] + dy_end = BC.domain.cellsize._y[-1] + dz_1 = BC.domain.cellsize._z[0] + dz_end = BC.domain.cellsize._z[-1] + rp = BC.domain.cellcenters._x[:, np.newaxis] i_ind = int_range(1,Nx)[:, np.newaxis, np.newaxis] j_ind = int_range(1,Ny)[np.newaxis, :, np.newaxis] diff --git a/src/pyfvtool/calculus.py b/src/pyfvtool/calculus.py index ee5d5e8..521a5e9 100644 --- a/src/pyfvtool/calculus.py +++ b/src/pyfvtool/calculus.py @@ -28,35 +28,35 @@ def gradientTerm(phi: CellVariable): >>> m = pf.Grid1D(10, 1.0) >>> phi = pf.CellVariable(m, 1.0) >>> gradPhi = pf.gradientTerm(phi) - >>> gradPhi.xvalue + >>> gradPhi._xvalue """ # calculates the gradient of a variable # the output is a face variable if issubclass(type(phi.domain), Grid1D): - dx = 0.5*(phi.domain.cellsize.x[0:-1]+phi.domain.cellsize.x[1:]) + dx = 0.5*(phi.domain.cellsize._x[0:-1]+phi.domain.cellsize._x[1:]) return FaceVariable(phi.domain, (phi.value[1:]-phi.value[0:-1])/dx, np.array([]), np.array([])) elif (type(phi.domain) is Grid2D) or (type(phi.domain) is CylindricalGrid2D): - dx = 0.5*(phi.domain.cellsize.x[0:-1]+phi.domain.cellsize.x[1:]) - dy = 0.5*(phi.domain.cellsize.y[0:-1]+phi.domain.cellsize.y[1:]) + dx = 0.5*(phi.domain.cellsize._x[0:-1]+phi.domain.cellsize._x[1:]) + dy = 0.5*(phi.domain.cellsize._y[0:-1]+phi.domain.cellsize._y[1:]) return FaceVariable(phi.domain, (phi.value[1:, 1:-1]-phi.value[0:-1, 1:-1])/dx[:,np.newaxis], (phi.value[1:-1, 1:]-phi.value[1:-1, 0:-1])/dy, np.array([])) elif (type(phi.domain) is PolarGrid2D): - dx = 0.5*(phi.domain.cellsize.x[0:-1]+phi.domain.cellsize.x[1:]) - dtheta = 0.5*(phi.domain.cellsize.y[0:-1]+phi.domain.cellsize.y[1:]) - rp = phi.domain.cellcenters.x + dx = 0.5*(phi.domain.cellsize._x[0:-1]+phi.domain.cellsize._x[1:]) + dtheta = 0.5*(phi.domain.cellsize._y[0:-1]+phi.domain.cellsize._y[1:]) + rp = phi.domain.cellcenters._x return FaceVariable(phi.domain, (phi.value[1:, 1:-1]-phi.value[0:-1, 1:-1])/dx[:,np.newaxis], (phi.value[1:-1, 1:]-phi.value[1:-1, 0:-1])/(dtheta[np.newaxis,:]*rp[:,np.newaxis]), np.array([])) elif (type(phi.domain) is Grid3D): - dx = 0.5*(phi.domain.cellsize.x[0:-1]+phi.domain.cellsize.x[1:]) - dy = 0.5*(phi.domain.cellsize.y[0:-1]+phi.domain.cellsize.y[1:]) - dz = 0.5*(phi.domain.cellsize.z[0:-1]+phi.domain.cellsize.z[1:]) + dx = 0.5*(phi.domain.cellsize._x[0:-1]+phi.domain.cellsize._x[1:]) + dy = 0.5*(phi.domain.cellsize._y[0:-1]+phi.domain.cellsize._y[1:]) + dz = 0.5*(phi.domain.cellsize._z[0:-1]+phi.domain.cellsize._z[1:]) return FaceVariable(phi.domain, (phi.value[1:, 1:-1, 1:-1] - phi.value[0:-1, 1:-1, 1:-1])/dx[:,np.newaxis,np.newaxis], @@ -65,10 +65,10 @@ def gradientTerm(phi: CellVariable): (phi.value[1:-1, 1:-1, 1:] - phi.value[1:-1, 1:-1, 0:-1])/dz[np.newaxis,np.newaxis,:]) elif (type(phi.domain) is CylindricalGrid3D): - dx = 0.5*(phi.domain.cellsize.x[0:-1]+phi.domain.cellsize.x[1:]) - dy = 0.5*(phi.domain.cellsize.y[0:-1]+phi.domain.cellsize.y[1:]) - dz = 0.5*(phi.domain.cellsize.z[0:-1]+phi.domain.cellsize.z[1:]) - rp = phi.domain.cellcenters.x + dx = 0.5*(phi.domain.cellsize._x[0:-1]+phi.domain.cellsize._x[1:]) + dy = 0.5*(phi.domain.cellsize._y[0:-1]+phi.domain.cellsize._y[1:]) + dz = 0.5*(phi.domain.cellsize._z[0:-1]+phi.domain.cellsize._z[1:]) + rp = phi.domain.cellcenters._x return FaceVariable(phi.domain, (phi.value[1:, 1:-1, 1:-1] - phi.value[0:-1, 1:-1, 1:-1])/dx[:,np.newaxis,np.newaxis], @@ -84,11 +84,11 @@ def divergenceTerm1D(F: FaceVariable): # extract data from the mesh structure Nx = F.domain.dims[0] G = F.domain.cell_numbers() - DX = F.domain.cellsize.x[1:-1] + DX = F.domain.cellsize._x[1:-1] # define the vector of cell index row_index = G[1:Nx+1] # main diagonal # compute the divergence - div_x = (F.xvalue[1:Nx+1]-F.xvalue[0:Nx])/DX + div_x = (F._xvalue[1:Nx+1]-F._xvalue[0:Nx])/DX # define the RHS Vector RHSdiv = np.zeros(Nx+2) # assign the values of the RHS vector @@ -103,15 +103,15 @@ def divergenceTermCylindrical1D(F: FaceVariable): # extract data from the mesh structure Nx = F.domain.dims[0] G = F.domain.cell_numbers() - DX = F.domain.cellsize.x[1:-1] - rp = F.domain.cellcenters.x - rf = F.domain.facecenters.x + DX = F.domain.cellsize._x[1:-1] + rp = F.domain.cellcenters._x + rf = F.domain.facecenters._x # define the vector of cell index row_index = G[1:Nx+1] # main diagonal # reassign the east, west, north, and south flux vectors for the # code readability - Fe = F.xvalue[1:Nx+1] - Fw = F.xvalue[0:Nx] + Fe = F._xvalue[1:Nx+1] + Fw = F._xvalue[0:Nx] re = rf[1:Nx+1] rw = rf[0:Nx] # compute the divergence @@ -129,16 +129,16 @@ def divergenceTerm2D(F: FaceVariable): # extract data from the mesh structure Nx, Ny = F.domain.dims G= F.domain.cell_numbers() - DX = F.domain.cellsize.x[1:-1][:, np.newaxis] - DY = F.domain.cellsize.y[1:-1][np.newaxis, :] + DX = F.domain.cellsize._x[1:-1][:, np.newaxis] + DY = F.domain.cellsize._y[1:-1][np.newaxis, :] # define the vector of cell index row_index = G[1:Nx+1,1:Ny+1].ravel() # main diagonal # reassign the east, west, north, and south flux vectors for the # code readability - Fe = F.xvalue[1:Nx+1,:] - Fw = F.xvalue[0:Nx,:] - Fn = F.yvalue[:,1:Ny+1] - Fs = F.yvalue[:,0:Ny] + Fe = F._xvalue[1:Nx+1,:] + Fw = F._xvalue[0:Nx,:] + Fn = F._yvalue[:,1:Ny+1] + Fs = F._yvalue[:,0:Ny] # compute the divergence div_x = (Fe - Fw)/DX div_y = (Fn - Fs)/DY @@ -161,18 +161,18 @@ def divergenceTermCylindrical2D(F:FaceVariable): # extract data from the mesh structure Nr, Nz = F.domain.dims G= F.domain.cell_numbers() - dr = F.domain.cellsize.x[1:-1][:, np.newaxis] - dz = F.domain.cellsize.y[1:-1][np.newaxis, :] - rp = F.domain.cellcenters.x[:, np.newaxis] - rf = F.domain.facecenters.x[:, np.newaxis] + dr = F.domain.cellsize._x[1:-1][:, np.newaxis] + dz = F.domain.cellsize._y[1:-1][np.newaxis, :] + rp = F.domain.cellcenters._x[:, np.newaxis] + rf = F.domain.facecenters._x[:, np.newaxis] # define the vector of cell index row_index = G[1:Nr+1,1:Nz+1].ravel() # main diagonal # reassign the east, west, north, and south flux vectors for the # code readability - Fe = F.xvalue[1:Nr+1,:] - Fw = F.xvalue[0:Nr,:] - Fn = F.yvalue[:,1:Nz+1] - Fs = F.yvalue[:,0:Nz] + Fe = F._xvalue[1:Nr+1,:] + Fw = F._xvalue[0:Nr,:] + Fn = F._yvalue[:,1:Nz+1] + Fs = F._yvalue[:,0:Nz] re = rf[1:Nr+1] rw = rf[0:Nr] # compute the divergence @@ -196,18 +196,18 @@ def divergenceTermPolar2D(F:FaceVariable): # extract data from the mesh structure Nr, Ntheta = F.domain.dims G=F.domain.cell_numbers() - dr = F.domain.cellsize.x[1:-1][:, np.newaxis] - dtheta= F.domain.cellsize.y[1:-1][np.newaxis, :] - rp = F.domain.cellcenters.x[:, np.newaxis] - rf = F.domain.facecenters.x[:, np.newaxis] + dr = F.domain.cellsize._x[1:-1][:, np.newaxis] + dtheta= F.domain.cellsize._y[1:-1][np.newaxis, :] + rp = F.domain.cellcenters._x[:, np.newaxis] + rf = F.domain.facecenters._x[:, np.newaxis] # define the vector of cell index row_index = G[1:Nr+1,1:Ntheta+1].ravel() # main diagonal # reassign the east, west, north, and south flux vectors for the # code readability - Fe = F.xvalue[1:Nr+1,:] - Fw = F.xvalue[0:Nr,:] - Fn = F.yvalue[:,1:Ntheta+1] - Fs = F.yvalue[:,0:Ntheta] + Fe = F._xvalue[1:Nr+1,:] + Fw = F._xvalue[0:Nr,:] + Fn = F._yvalue[:,1:Ntheta+1] + Fs = F._yvalue[:,0:Ntheta] re = rf[1:Nr+1] rw = rf[0:Nr] # compute the divergence @@ -230,19 +230,19 @@ def divergenceTerm3D(F:FaceVariable): # extract data from the mesh structure Nx, Ny, Nz = F.domain.dims G=F.domain.cell_numbers() - dx = F.domain.cellsize.x[1:-1][:,np.newaxis,np.newaxis] - dy = F.domain.cellsize.y[1:-1][np.newaxis,:,np.newaxis] - dz = F.domain.cellsize.z[1:-1][np.newaxis,np.newaxis,:] + dx = F.domain.cellsize._x[1:-1][:,np.newaxis,np.newaxis] + dy = F.domain.cellsize._y[1:-1][np.newaxis,:,np.newaxis] + dz = F.domain.cellsize._z[1:-1][np.newaxis,np.newaxis,:] # define the vector of cell index row_index = G[1:Nx+1,1:Ny+1,1:Nz+1].ravel() # main diagonal # reassign the east, west, north, and south flux vectors for the # code readability - Fe = F.xvalue[1:Nx+1,:,:] - Fw = F.xvalue[0:Nx,:,:] - Fn = F.yvalue[:,1:Ny+1,:] - Fs = F.yvalue[:,0:Ny,:] - Ff = F.zvalue[:,:,1:Nz+1] - Fb = F.zvalue[:,:,0:Nz] + Fe = F._xvalue[1:Nx+1,:,:] + Fw = F._xvalue[0:Nx,:,:] + Fn = F._yvalue[:,1:Ny+1,:] + Fs = F._yvalue[:,0:Ny,:] + Ff = F._zvalue[:,:,1:Nz+1] + Fb = F._zvalue[:,:,0:Nz] # compute the divergence div_x = (Fe - Fw)/dx div_y = (Fn - Fs)/dy @@ -266,20 +266,20 @@ def divergenceTermCylindrical3D(F:FaceVariable): # extract data from the mesh structure Nx, Ny, Nz = F.domain.dims G=F.domain.cell_numbers() - dx = F.domain.cellsize.x[1:-1][:,np.newaxis,np.newaxis] - dy = F.domain.cellsize.y[1:-1][np.newaxis,:,np.newaxis] - dz = F.domain.cellsize.z[1:-1][np.newaxis,np.newaxis,:] - rp = F.domain.cellcenters.x[:,np.newaxis,np.newaxis] + dx = F.domain.cellsize._x[1:-1][:,np.newaxis,np.newaxis] + dy = F.domain.cellsize._y[1:-1][np.newaxis,:,np.newaxis] + dz = F.domain.cellsize._z[1:-1][np.newaxis,np.newaxis,:] + rp = F.domain.cellcenters._x[:,np.newaxis,np.newaxis] # define the vector of cell index row_index = G[1:Nx+1,1:Ny+1,1:Nz+1].ravel() # main diagonal # reassign the east, west, north, and south flux vectors for the # code readability - Fe = F.xvalue[1:Nx+1,:,:] - Fw = F.xvalue[0:Nx,:,:] - Fn = F.yvalue[:,1:Ny+1,:] - Fs = F.yvalue[:,0:Ny,:] - Ff = F.zvalue[:,:,1:Nz+1] - Fb = F.zvalue[:,:,0:Nz] + Fe = F._xvalue[1:Nx+1,:,:] + Fw = F._xvalue[0:Nx,:,:] + Fn = F._yvalue[:,1:Ny+1,:] + Fs = F._yvalue[:,0:Ny,:] + Ff = F._zvalue[:,:,1:Nz+1] + Fb = F._zvalue[:,:,0:Nz] # compute the divergence div_x = (Fe - Fw)/dx div_y = (Fn - Fs)/(dy*rp) @@ -336,16 +336,26 @@ def divergenceTerm(F: FaceVariable): def gradientTermFixedBC(phi): """ - Warning: unless you know for sure that you need this function, do not use it! - This function calculates the gradient of parameter phi in x,y, and z directions. It takes care of the often nonphysical - values of the ghost cells. Note that phi is not a variable but a parameter calculated with a function over a domain. - Make sure that phi is calculated by BC2GhostCells (usually but not necessarily in combination with celleval); - otherwise, do not use this function as it leads to wrong values at the boundaries. - It checks for the availability of the ghost variables and use them, otherwise estimate them, assuming a zero gradient - on the boundaries. - Note: I'm not happy with this implementation but it was the fastest solution that came into my mind while onboard the Geilo-Oslo train. - I have to find a better way to do this. The problem is that it is almost always used for a cell variable calculated as f(phi) so having a boundary condition - does not really help. I have to think about it. + Warning: + + Unless you know for sure that you need this function, do not use it! + This function calculates the gradient of parameter phi in x,y, and z + directions. It takes care of the often nonphysical values of the ghost + cells. Note that phi is not a variable but a parameter calculated with a + function over a domain. + + Make sure that phi is calculated by BC2GhostCells (usually but not + necessarily in combination with celleval); otherwise, do not use this + function as it leads to wrong values at the boundaries. + + It checks for the availability of the ghost variables and use them, + otherwise estimate them, assuming a zero gradient on the boundaries. + + Note: I'm not happy with this implementation but it was the fastest + solution that came into my mind while onboard the Geilo-Oslo train. + I have to find a better way to do this. The problem is that it is almost + always used for a cell variable calculated as f(phi) so having a boundary + condition does not really help. I have to think about it. parameters ---------- @@ -363,23 +373,23 @@ def gradientTermFixedBC(phi): >>> import numpy as np >>> m = pf.Grid1D(10, 1.0) >>> phi = pf.CellVariable(m, 1.0) - >>> sin_phi = pf.celleval(np.sin, BC2GhostCells(sw)) + >>> sin_phi = pf.celleval(np.sin, sw.BC2GhostCells()) >>> gradPhi = pf.gradientTermFixedBC(sin_phi) """ faceGrad = gradientTerm(phi) if issubclass(type(phi.domain), Grid1D): - faceGrad.xvalue[0] = 2*faceGrad.xvalue[0] - faceGrad.xvalue[-1] = 2*faceGrad.xvalue[-1] + faceGrad._xvalue[0] = 2*faceGrad._xvalue[0] + faceGrad._xvalue[-1] = 2*faceGrad._xvalue[-1] elif issubclass(type(phi.domain), Grid2D): - faceGrad.xvalue[0, :] = 2*faceGrad.xvalue[0, :] - faceGrad.xvalue[-1, :] = 2*faceGrad.xvalue[-1, :] - faceGrad.yvalue[:, 0] = 2*faceGrad.yvalue[:, 0] - faceGrad.yvalue[:, -1] = 2*faceGrad.yvalue[:, -1] + faceGrad._xvalue[0, :] = 2*faceGrad._xvalue[0, :] + faceGrad._xvalue[-1, :] = 2*faceGrad._xvalue[-1, :] + faceGrad._yvalue[:, 0] = 2*faceGrad._yvalue[:, 0] + faceGrad._yvalue[:, -1] = 2*faceGrad._yvalue[:, -1] elif issubclass(type(phi.domain), Grid3D): - faceGrad.xvalue[0, :, :] = 2*faceGrad.xvalue[0, :, :] - faceGrad.xvalue[-1, :, :] = 2*faceGrad.xvalue[-1, :, :] - faceGrad.yvalue[:, 0, :] = 2*faceGrad.yvalue[:, 0, :] - faceGrad.yvalue[:, -1, :] = 2*faceGrad.yvalue[:, -1, :] - faceGrad.zvalue[:, :, 0] = 2*faceGrad.zvalue[:, :, 0] - faceGrad.zvalue[:, :, -1] = 2*faceGrad.zvalue[:, :, -1] + faceGrad._xvalue[0, :, :] = 2*faceGrad._xvalue[0, :, :] + faceGrad._xvalue[-1, :, :] = 2*faceGrad._xvalue[-1, :, :] + faceGrad._yvalue[:, 0, :] = 2*faceGrad._yvalue[:, 0, :] + faceGrad._yvalue[:, -1, :] = 2*faceGrad._yvalue[:, -1, :] + faceGrad._zvalue[:, :, 0] = 2*faceGrad._zvalue[:, :, 0] + faceGrad._zvalue[:, :, -1] = 2*faceGrad._zvalue[:, :, -1] return faceGrad \ No newline at end of file diff --git a/src/pyfvtool/cell.py b/src/pyfvtool/cell.py index 6f1996e..b23b707 100644 --- a/src/pyfvtool/cell.py +++ b/src/pyfvtool/cell.py @@ -10,6 +10,8 @@ from .boundary import BoundaryConditionsBase, BoundaryConditions from .boundary import cellValuesWithBoundaries + + class CellVariable: @overload @@ -103,39 +105,6 @@ def internalCellValues(self, values): self.value[1:-1, 1:-1] = values elif issubclass(type(self.domain), Grid3D): self.value[1:-1, 1:-1, 1:-1] = values - - def update_bc_cells(self, BC: BoundaryConditionsBase): - phi_temp = CellVariable(self.domain, self.internalCellValues, BC) - self.update_value(phi_temp) - - def update_value(self, new_cell): - np.copyto(self.value, new_cell.value) - - def bc_to_ghost(self): - """ - assign the boundary values to the ghost cells - """ - if issubclass(type(self.domain), Grid1D): - self.value[0] = 0.5*(self.value[1]+self.value[0]) - self.value[-1] = 0.5*(self.value[-2]+self.value[-1]) - elif issubclass(type(self.domain), Grid2D): - self.value[0, 1:-1] = 0.5*(self.value[1, 1:-1]+self.value[0, 1:-1]) - self.value[-1, 1:-1] = 0.5*(self.value[-2, 1:-1]+self.value[-1, 1:-1]) - self.value[1:-1, 0] = 0.5*(self.value[1:-1, 1]+self.value[1:-1, 0]) - self.value[1:-1, -1] = 0.5*(self.value[1:-1, -2]+self.value[1:-1, -1]) - elif issubclass(type(self.domain), Grid3D): - self.value[0, 1:-1, 1:-1] = 0.5*(self.value[1, 1:-1, 1:-1]+self.value[0, 1:-1, 1:-1]) - self.value[-1, 1:-1, 1:-1] = 0.5*(self.value[-2, 1:-1, 1:-1]+self.value[-1, 1:-1, 1:-1]) - self.value[1:-1, 0, 1:-1] = 0.5*(self.value[1:-1, 1, 1:-1]+self.value[1:-1, 0, 1:-1]) - self.value[1:-1, -1, 1:-1] = 0.5*(self.value[1:-1, -2, 1:-1]+self.value[1:-1, -1, 1:-1]) - self.value[1:-1, 1:-1, 0] = 0.5*(self.value[1:-1, 1:-1, 1]+self.value[1:-1, 1:-1, 0]) - self.value[1:-1, 1:-1, -1] = 0.5*(self.value[1:-1, 1:-1, -2]+self.value[1:-1, 1:-1, -1]) - - def copy(self): - """ - Create a copy of the CellVariable - """ - return CellVariable(self.domain, np.copy(self.value)) def __add__(self, other): if type(other) is CellVariable: @@ -240,43 +209,224 @@ def __abs__(self): return CellVariable(self.domain, np.abs(self.value)) + def update_bc_cells(self, BC: BoundaryConditionsBase): + phi_temp = CellVariable(self.domain, self.internalCellValues, BC) + self.update_value(phi_temp) + def update_value(self, new_cell): + np.copyto(self.value, new_cell.value) -def copyCellVariable(phi: CellVariable) -> CellVariable: - """ - Create a copy of a CellVariable + def bc_to_ghost(self): + """ + assign the boundary values to the ghost cells + """ + if issubclass(type(self.domain), Grid1D): + self.value[0] = 0.5*(self.value[1]+self.value[0]) + self.value[-1] = 0.5*(self.value[-2]+self.value[-1]) + elif issubclass(type(self.domain), Grid2D): + self.value[0, 1:-1] = 0.5*(self.value[1, 1:-1]+self.value[0, 1:-1]) + self.value[-1, 1:-1] = 0.5*(self.value[-2, 1:-1]+self.value[-1, 1:-1]) + self.value[1:-1, 0] = 0.5*(self.value[1:-1, 1]+self.value[1:-1, 0]) + self.value[1:-1, -1] = 0.5*(self.value[1:-1, -2]+self.value[1:-1, -1]) + elif issubclass(type(self.domain), Grid3D): + self.value[0, 1:-1, 1:-1] = 0.5*(self.value[1, 1:-1, 1:-1]+self.value[0, 1:-1, 1:-1]) + self.value[-1, 1:-1, 1:-1] = 0.5*(self.value[-2, 1:-1, 1:-1]+self.value[-1, 1:-1, 1:-1]) + self.value[1:-1, 0, 1:-1] = 0.5*(self.value[1:-1, 1, 1:-1]+self.value[1:-1, 0, 1:-1]) + self.value[1:-1, -1, 1:-1] = 0.5*(self.value[1:-1, -2, 1:-1]+self.value[1:-1, -1, 1:-1]) + self.value[1:-1, 1:-1, 0] = 0.5*(self.value[1:-1, 1:-1, 1]+self.value[1:-1, 1:-1, 0]) + self.value[1:-1, 1:-1, -1] = 0.5*(self.value[1:-1, 1:-1, -2]+self.value[1:-1, 1:-1, -1]) + + def copy(self): + """ + Create a copy of the CellVariable + + + Returns + ------- + CellVariable + Copy of the CellVariable. + """ + return CellVariable(self.domain, np.copy(self.value)) + + def plotprofile(self): + """ + Create a 'profile' of a cell variable for plotting, export, etc. + + 'Plot profiles' are sets of arrays that contain the axes' coordinates, + cell values, including the values at the boundaries. + + For 2D and 3D visualization, it is perhaps better to use only the + internalCellValues (e.g. for a false color map à la plt.pcolormesh), + and not include the values at the boundaries. + + + 1D meshes + ========= + This generates a pair of vectors containing the abscissa and ordinates + for plotting the values of the cell variable over the entire calculation + domain. It includes the values at the outer faces of the domain, by + taking into account the values of the ghost cells. + + Returns + ------- + x : np.ndarray + x (or r) coordinates. + phi0 : np.ndarray + Value of the CellVariables at those points. + + + 2D meshes + ========= + This generates a set of vectors containing the (x, y) or (r, z) + coordinates and the values of the cell variable at those coordinates + for plotting the values of the cell variable over the entire calculation + domain. It includes the values at the outer faces of the domain, by + taking into account the values of the ghost cells. + + Returns + ------- + x : np.ndarray + x (or r) coordinates. + y : np.ndarray + y (or z) coordinates. + phi0 : np.ndarray + Value of the CellVariables at those coordinates. + + + 3D meshes + ========= + This generates a set of vectors containing the (x, y, z) or (r, theta, z) + coordinates and the values of the cell variable at those coordinates + for plotting the values of the cell variable over the entire calculation + domain. It includes the values at the outer faces of the domain, by + taking into account the values of the ghost cells. + + Returns + ------- + x : np.ndarray + x (or r) coordinates. + y : np.ndarray + y (or theta) coordinates. + z : np.ndarray + z coordinates. + phi0 : np.ndarray + Value of the CellVariables at those coordinates. + """ + # + # TODO: + # Add a keyword that specifies how the outer FaceValues will be estimated + # Currently, this is just the average of the last inner cell and the boundary + # (ghost) cell. + # In certain cases it may be visually desirable to use an extrapolation of + # the last inner cell values. + # + if isinstance(self.domain, Grid1D): + x = np.hstack([self.domain.facecenters._x[0], + self.domain.cellcenters._x, + self.domain.facecenters._x[-1]]) + phi0 = np.hstack([0.5*(self.value[0]+self.value[1]), + self.value[1:-1], + 0.5*(self.value[-2]+self.value[-1])]) + # The size of the ghost cell is always equal to the size of the + # first (or last) cell within the domain. The value at the + # boundary can therefore be obtained by direct averaging with a + # weight factor of 0.5. + return (x, phi0) + elif isinstance(self.domain, Grid2D): + x = np.hstack([self.domain.facecenters._x[0], + self.domain.cellcenters._x, + self.domain.facecenters._x[-1]]) + y = np.hstack([self.domain.facecenters._y[0], + self.domain.cellcenters._y, + self.domain.facecenters._y[-1]]) + phi0 = np.copy(self.value) + phi0[:, 0] = 0.5*(phi0[:, 0]+phi0[:, 1]) + phi0[0, :] = 0.5*(phi0[0, :]+phi0[1, :]) + phi0[:, -1] = 0.5*(phi0[:, -1]+phi0[:, -2]) + phi0[-1, :] = 0.5*(phi0[-1, :]+phi0[-2, :]) + phi0[0, 0] = phi0[0, 1] + phi0[0, -1] = phi0[0, -2] + phi0[-1, 0] = phi0[-1, 1] + phi0[-1, -1] = phi0[-1, -2] + return (x, y, phi0) + elif isinstance(self.domain, Grid3D): + x = np.hstack([self.domain.facecenters._x[0], + self.domain.cellcenters._x, + self.domain.facecenters._x[-1]])[:, np.newaxis, np.newaxis] + y = np.hstack([self.domain.facecenters._y[0], + self.domain.cellcenters._y, + self.domain.facecenters._y[-1]])[np.newaxis, :, np.newaxis] + z = np.hstack([self.domain.facecenters._z[0], + self.domain.cellcenters._z, + self.domain.facecenters._z[-1]])[np.newaxis, np.newaxis, :] + phi0 = np.copy(self.value) + phi0[:,0,:]=0.5*(phi0[:,0,:]+phi0[:,1,:]) + phi0[:,-1,:]=0.5*(phi0[:,-2,:]+phi0[:,-1,:]) + phi0[:,:,0]=0.5*(phi0[:,:,0]+phi0[:,:,0]) + phi0[:,:,-1]=0.5*(phi0[:,:,-2]+phi0[:,:,-1]) + phi0[0,:,:]=0.5*(phi0[1,:,:]+phi0[2,:,:]) + phi0[-1,:,:]=0.5*(phi0[-2,:,:]+phi0[-1,:,:]) + return (x, y, z, phi0) + else: + raise NotImplementedError("plotprofile() not implemented for mesh type '{0:s}'".\ + format(self.domain.__class__.__name__)) + + + def domainIntegral(self) -> float: + """ + Calculate the finite-volume integral of a CellVariable over entire domain + + The finite-volume integral over the entire mesh domain gives the total + amount of `CellVariable` present in the system. Calculation of this + integral is useful for checking conservation of the quantity concerned + (in case of 'no-flux' BCs), or for monitoring its evolution due to + exchanges via the boundaries (other BCs) or to the presence of source + terms. + + May later become a built-in method of the CellVariable class, but for now + this implementation as a function is chosen for consistency with FVTool. + + + Returns + ------- + float + Total finite-volume integral over entire domain. + + """ + v = self.domain.cellVolumes() + c = self.internalCellValues + return (v*c).flatten().sum() + + + def BC2GhostCells(self): + """ + assign the boundary values to the ghost cells + + Returns + ------- + CellVariable + the new cell variable + """ + phi = self.copy() + if issubclass(type(phi.domain), Grid1D): + phi.value[0] = 0.5*(phi.value[1]+phi.value[0]) + phi.value[-1] = 0.5*(phi.value[-2]+phi.value[-1]) + elif issubclass(type(phi.domain), Grid2D): + phi.value[0, 1:-1] = 0.5*(phi.value[1, 1:-1]+phi.value[0, 1:-1]) + phi.value[-1, 1:-1] = 0.5*(phi.value[-2, 1:-1]+phi.value[-1, 1:-1]) + phi.value[1:-1, 0] = 0.5*(phi.value[1:-1, 1]+phi.value[1:-1, 0]) + phi.value[1:-1, -1] = 0.5*(phi.value[1:-1, -2]+phi.value[1:-1, -1]) + elif issubclass(type(phi.domain), Grid3D): + phi.value[0, 1:-1, 1:-1] = 0.5*(phi.value[1, 1:-1, 1:-1]+phi.value[0, 1:-1, 1:-1]) + phi.value[-1, 1:-1, 1:-1] = 0.5*(phi.value[-2, 1:-1, 1:-1]+phi.value[-1, 1:-1, 1:-1]) + phi.value[1:-1, 0, 1:-1] = 0.5*(phi.value[1:-1, 1, 1:-1]+phi.value[1:-1, 0, 1:-1]) + phi.value[1:-1, -1, 1:-1] = 0.5*(phi.value[1:-1, -2, 1:-1]+phi.value[1:-1, -1, 1:-1]) + phi.value[1:-1, 1:-1, 0] = 0.5*(phi.value[1:-1, 1:-1, 1]+phi.value[1:-1, 1:-1, 0]) + phi.value[1:-1, 1:-1, -1] = 0.5*(phi.value[1:-1, 1:-1, -2]+phi.value[1:-1, 1:-1, -1]) + return phi - Parameters - ---------- - phi : CellVariable - CellVariable to be copied. - Returns - ------- - CellVariable - Copy of the CellVariable. - """ - return CellVariable(phi.domain, np.copy(phi.value)) - - -def cellVolume(m: MeshStructure): - BC = BoundaryConditions(m) - if (type(m) is Grid1D): - c=m.cellsize.x[1:-1] - elif (type(m) is CylindricalGrid1D): - c=2.0*np.pi*m.cellsize.x[1:-1]*m.cellcenters.x - elif (type(m) is Grid2D): - c=m.cellsize.x[1:-1][:, np.newaxis]*m.cellsize.y[1:-1][np.newaxis, :] - elif (type(m) is CylindricalGrid2D): - c=2.0*np.pi*m.cellcenters.x[:, np.newaxis]*m.cellsize.x[1:-1][:, np.newaxis]*m.cellsize.y[1:-1][np.newaxis, :] - elif (type(m) is PolarGrid2D): - c=m.cellcenters.x*m.cellsize.x[1:-1][:, np.newaxis]*m.cellsize.y[1:-1][np.newaxis, :] - elif (type(m) is Grid3D): - c=m.cellsize.x[1:-1][:,np.newaxis,np.newaxis]*m.cellsize.y[1:-1][np.newaxis,:,np.newaxis]*m.cellsize.z[1:-1][np.newaxis,np.newaxis,:] - elif (type(m) is CylindricalGrid3D): - c=m.cellcenters.x*m.cellsize.x[1:-1][:,np.newaxis,np.newaxis]*m.cellsize.y[1:-1][np.newaxis,:,np.newaxis]*m.cellsize.z[np.newaxis,np.newaxis,:] - return CellVariable(m, c, BC) def cellLocations(m: MeshStructure): @@ -285,7 +435,8 @@ def cellLocations(m: MeshStructure): It can later be used in defining properties that are variable in space. - Incompletely tested + Incompletely tested, and there may be other, more direct, ways to + calculate properties that vary in space, e.g. using cellcenters directly? Parameters ---------- @@ -318,26 +469,27 @@ def cellLocations(m: MeshStructure): if (type(m) is Grid1D)\ or (type(m) is CylindricalGrid1D): - X = CellVariable(m, m.cellcenters.x) + X = CellVariable(m, m.cellcenters._x) return X elif (type(m) is Grid2D)\ or (type(m) is CylindricalGrid2D)\ or (type(m) is PolarGrid2D): - X = CellVariable(m, np.tile(m.cellcenters.x[:, np.newaxis], (1, N[1]))) - Y = CellVariable(m, np.tile(m.cellcenters.y[:, np.newaxis].T, (N[0], 1))) + X = CellVariable(m, np.tile(m.cellcenters._x[:, np.newaxis], (1, N[1]))) + Y = CellVariable(m, np.tile(m.cellcenters._y[:, np.newaxis].T, (N[0], 1))) return X, Y elif (type(m) is Grid3D)\ or (type(m) is CylindricalGrid3D): - X = CellVariable(m, np.tile(m.cellcenters.x[:, np.newaxis, np.newaxis], (1, N[1], N[2]))) - Y = CellVariable(m, np.tile((m.cellcenters.y[:, np.newaxis].T)[:,:,np.newaxis], (N[0], 1, N[2]))) + X = CellVariable(m, np.tile(m.cellcenters._x[:, np.newaxis, np.newaxis], (1, N[1], N[2]))) + Y = CellVariable(m, np.tile((m.cellcenters._y[:, np.newaxis].T)[:,:,np.newaxis], (N[0], 1, N[2]))) z = np.zeros((1,1,N[2])) - z[0, 0, :] = m.cellcenters.z + z[0, 0, :] = m.cellcenters._z Z = CellVariable(m, np.tile(z, (N[0], N[1], 1))) return X, Y, Z raise TypeError('mesh type not implemented') return None + def funceval(f, *args): if len(args)==1: return CellVariable(args[0].domain, @@ -371,236 +523,3 @@ def celleval(f, *args): -def domainInt(phi: CellVariable) -> float: - """ - Calculate the finite-volume integral of a CellVariable over entire domain - - The finite-volume integral over the entire mesh domain gives the total - amount of `CellVariable` present in the system. Calculation of this - integral is useful for checking conservation of the quantity concerned - (in case of 'no-flux' BCs), or for monitoring its evolution due to - exchanges via the boundaries (other BCs) or to the presence of source - terms. - - May later become a built-in method of the CellVariable class, but for now - this implementation as a function is chosen for consistency with FVTool. - An alias `domainIntegrate`has been added that sounds better than - `domainInt`. - - Parameters - ---------- - phi : CellVariable - Variable whose finite-volume integral will calculated. - - Returns - ------- - float - Total finite-volume integral over entire domain. - - """ - v = cellVolume(phi.domain).internalCellValues - c = phi.internalCellValues - return (v*c).flatten().sum() - -def domainIntegrate(phi: CellVariable) -> float: - """ - Alias for `domainInt()` - - See `domainInt()` - """ - return domainInt(phi) - - -# TODO: -# get_CellVariable_profile1D can become a method of CellVariable -# (shared with 2D and 3D versions) -# TODO: -# Add a keyword that specifies how the outer FaceValues will be estimated -# Currently, this is just the average of the last inner cell and the boundary -# (ghost) cell. -# In certain cases it may be visually desirable to use an extrapolation of -# the last inner cell values. - -def get_CellVariable_profile1D(phi: CellVariable): - """ - Create a profile of a cell variable for plotting, export, etc. (1D). - - This generates a pair of vectors containing the abscissa and ordinates - for plotting the values of the cell variable over the entire calculation - domain. It includes the values at the outer faces of the domain, by - taking into account the values of the ghost cells. - - This function may later become a method of the CellVariable class, but - is a function now for simplicity and consistency with other CellVariable - utility functions (e.g. `domainIntegrate`). - - - Parameters - ---------- - phi : CellVariable - - - Returns - ------- - x : np.ndarray - x (or r) coordinates. - phi0 : np.ndarray - Value of the CellVariables at those points. - - """ - - x = np.hstack([phi.domain.facecenters.x[0], - phi.domain.cellcenters.x, - phi.domain.facecenters.x[-1]]) - phi0 = np.hstack([0.5*(phi.value[0]+phi.value[1]), - phi.value[1:-1], - 0.5*(phi.value[-2]+phi.value[-1])]) - # The size of the ghost cell is always equal to the size of the - # first (or last) cell within the domain. The value at the - # boundary can therefore be obtained by direct averaging with a - # weight factor of 0.5. - return (x, phi0) - - -# TODO: -# get_CellVariable_profile2D can become a method of CellVariable -# (shared with 1D and 3D versions) -# TODO: -# Add a keyword that specifies how the outer FaceValues will be estimated -# Currently, this is just the average of the last inner cell and the boundary -# (ghost) cell. -# In certain cases it may be visually desirable to use an extrapolation of -# the last inner cell values. -# Perhaps for 2D and 3D visualization it is perhaps best to only use the -# innervalues (e.g. for a false color map à la plt.pcolormesh) - -def get_CellVariable_profile2D(phi: CellVariable): - """ - Create a profile of a cell variable for plotting, export, etc. (2D). - - This generates a set of vectors containing the (x, y) or (r, z) - coordinates and the values of the cell variable at those coordinates - for plotting the values of the cell variable over the entire calculation - domain. It includes the values at the outer faces of the domain, by - taking into account the values of the ghost cells. - - This function may later become a method of the CellVariable class, but - is a function now for simplicity and consistency with other CellVariable - utility functions (e.g. `domainIntegrate`). - - Parameters - ---------- - phi : CellVariable - - - Returns - ------- - x : np.ndarray - x (or r) coordinates. - y : np.ndarray - y (or z) coordinates. - phi0 : np.ndarray - Value of the CellVariables at those coordinates. - - """ - x = np.hstack([phi.domain.facecenters.x[0], - phi.domain.cellcenters.x, - phi.domain.facecenters.x[-1]]) - y = np.hstack([phi.domain.facecenters.y[0], - phi.domain.cellcenters.y, - phi.domain.facecenters.y[-1]]) - phi0 = np.copy(phi.value) - phi0[:, 0] = 0.5*(phi0[:, 0]+phi0[:, 1]) - phi0[0, :] = 0.5*(phi0[0, :]+phi0[1, :]) - phi0[:, -1] = 0.5*(phi0[:, -1]+phi0[:, -2]) - phi0[-1, :] = 0.5*(phi0[-1, :]+phi0[-2, :]) - phi0[0, 0] = phi0[0, 1] - phi0[0, -1] = phi0[0, -2] - phi0[-1, 0] = phi0[-1, 1] - phi0[-1, -1] = phi0[-1, -2] - return (x, y, phi0) - - - -# TODO: -# get_CellVariable_profile3D can become a method of CellVariable -# (shared with 1D and 2D versions) -# TODO: -# Add a keyword that specifies how the outer FaceValues will be estimated -# Currently, this is just the average of the last inner cell and the boundary -# (ghost) cell. -# In certain cases it may be visually desirable to use an extrapolation of -# the last inner cell values. -# Perhaps for 2D and 3D visualization it is perhaps best to only use the -# innervalues (e.g. for a false color map à la plt.pcolormesh) - -def get_CellVariable_profile3D(phi: CellVariable): - """ - Create a profile of a cell variable for plotting, export, etc. (3D). - - This generates a set of vectors containing the (x, y, z) or (r, theta, z) - coordinates and the values of the cell variable at those coordinates - for plotting the values of the cell variable over the entire calculation - domain. It includes the values at the outer faces of the domain, by - taking into account the values of the ghost cells. - - This function may later become a method of the CellVariable class, but - is a function now for simplicity and consistency with other CellVariable - utility functions (e.g. `domainIntegrate`). - - Parameters - ---------- - phi : CellVariable - - - Returns - ------- - x : np.ndarray - x (or r) coordinates. - y : np.ndarray - y (or theta) coordinates. - z : np.ndarray - z coordinates. - phi0 : np.ndarray - Value of the CellVariables at those coordinates. - """ - - x = np.hstack([phi.domain.facecenters.x[0], - phi.domain.cellcenters.x, - phi.domain.facecenters.x[-1]])[:, np.newaxis, np.newaxis] - y = np.hstack([phi.domain.facecenters.y[0], - phi.domain.cellcenters.y, - phi.domain.facecenters.y[-1]])[np.newaxis, :, np.newaxis] - z = np.hstack([phi.domain.facecenters.z[0], - phi.domain.cellcenters.z, - phi.domain.facecenters.z[-1]])[np.newaxis, np.newaxis, :] - phi0 = np.copy(phi.value) - phi0[:,0,:]=0.5*(phi0[:,0,:]+phi0[:,1,:]) - phi0[:,-1,:]=0.5*(phi0[:,-2,:]+phi0[:,-1,:]) - phi0[:,:,0]=0.5*(phi0[:,:,0]+phi0[:,:,0]) - phi0[:,:,-1]=0.5*(phi0[:,:,-2]+phi0[:,:,-1]) - phi0[0,:,:]=0.5*(phi0[1,:,:]+phi0[2,:,:]) - phi0[-1,:,:]=0.5*(phi0[-2,:,:]+phi0[-1,:,:]) - return (x, y, z, phi0) - -def BC2GhostCells(phi0): - """ - assign the boundary values to the ghost cells and returns the new cell variable - """ - phi = copyCellVariable(phi0) - if issubclass(type(phi.domain), Grid1D): - phi.value[0] = 0.5*(phi.value[1]+phi.value[0]) - phi.value[-1] = 0.5*(phi.value[-2]+phi.value[-1]) - elif issubclass(type(phi.domain), Grid2D): - phi.value[0, 1:-1] = 0.5*(phi.value[1, 1:-1]+phi.value[0, 1:-1]) - phi.value[-1, 1:-1] = 0.5*(phi.value[-2, 1:-1]+phi.value[-1, 1:-1]) - phi.value[1:-1, 0] = 0.5*(phi.value[1:-1, 1]+phi.value[1:-1, 0]) - phi.value[1:-1, -1] = 0.5*(phi.value[1:-1, -2]+phi.value[1:-1, -1]) - elif issubclass(type(phi.domain), Grid3D): - phi.value[0, 1:-1, 1:-1] = 0.5*(phi.value[1, 1:-1, 1:-1]+phi.value[0, 1:-1, 1:-1]) - phi.value[-1, 1:-1, 1:-1] = 0.5*(phi.value[-2, 1:-1, 1:-1]+phi.value[-1, 1:-1, 1:-1]) - phi.value[1:-1, 0, 1:-1] = 0.5*(phi.value[1:-1, 1, 1:-1]+phi.value[1:-1, 0, 1:-1]) - phi.value[1:-1, -1, 1:-1] = 0.5*(phi.value[1:-1, -2, 1:-1]+phi.value[1:-1, -1, 1:-1]) - phi.value[1:-1, 1:-1, 0] = 0.5*(phi.value[1:-1, 1:-1, 1]+phi.value[1:-1, 1:-1, 0]) - phi.value[1:-1, 1:-1, -1] = 0.5*(phi.value[1:-1, 1:-1, -2]+phi.value[1:-1, 1:-1, -1]) - return phi \ No newline at end of file diff --git a/src/pyfvtool/diffusion.py b/src/pyfvtool/diffusion.py index 0e64c9f..928704d 100644 --- a/src/pyfvtool/diffusion.py +++ b/src/pyfvtool/diffusion.py @@ -14,12 +14,12 @@ def diffusionTerm1D(D: FaceVariable) -> csr_array: # extract data from the mesh structure Nx = D.domain.dims[0] G = D.domain.cell_numbers() - DX = D.domain.cellsize.x + DX = D.domain.cellsize._x dx = 0.5*(DX[0:-1]+DX[1:]) # extract the velocity data # note: size(Dx) = [1:m+1, 1:n] and size(Dy) = [1:m, 1:n+1] - Dx = D.xvalue + Dx = D._xvalue # reassign the east, west, north, and south velocity vectors for the # code readability @@ -45,14 +45,14 @@ def diffusionTermCylindrical1D(D: FaceVariable) -> csr_array: # extract data from the mesh structure Nx = D.domain.dims[0] G = D.domain.cell_numbers() - DX = D.domain.cellsize.x + DX = D.domain.cellsize._x dx = 0.5*(DX[0:-1]+DX[1:]) - rp = D.domain.cellcenters.x - rf = D.domain.facecenters.x + rp = D.domain.cellcenters._x + rf = D.domain.facecenters._x # extract the velocity data # note: size(Dx) = [1:m+1, 1:n] and size(Dy) = [1:m, 1:n+1] - Dx = D.xvalue + Dx = D._xvalue # reassign the east, west, north, and south velocity vectors for the # code readability @@ -79,17 +79,17 @@ def diffusionTerm2D(D: FaceVariable) -> csr_array: # extract data from the mesh structure Nx, Ny = D.domain.dims G = D.domain.cell_numbers() - DX = D.domain.cellsize.x - DY = D.domain.cellsize.y + DX = D.domain.cellsize._x + DY = D.domain.cellsize._y dx = 0.5*(DX[0:-1]+DX[1:]) dy = 0.5*(DY[0:-1]+DY[1:]) mn = Nx*Ny # reassign the east, west for code readability (use broadcasting in Julia) - De = D.xvalue[1:Nx+1, :]/(dx[1:Nx+1]*DX[1:Nx+1])[:, np.newaxis] - Dw = D.xvalue[0:Nx, :]/(dx[0:Nx]*DX[1:Nx+1])[:, np.newaxis] - Dn = D.yvalue[:, 1:Ny+1]/(dy[1:Ny+1]*DY[1:Ny+1]) - Ds = D.yvalue[:, 0:Ny]/(dy[0:Ny]*DY[1:Ny+1]) + De = D._xvalue[1:Nx+1, :]/(dx[1:Nx+1]*DX[1:Nx+1])[:, np.newaxis] + Dw = D._xvalue[0:Nx, :]/(dx[0:Nx]*DX[1:Nx+1])[:, np.newaxis] + Dn = D._yvalue[:, 1:Ny+1]/(dy[1:Ny+1]*DY[1:Ny+1]) + Ds = D._yvalue[:, 0:Ny]/(dy[0:Ny]*DY[1:Ny+1]) # calculate the coefficients for the internal cells AE = De.ravel() @@ -126,20 +126,20 @@ def diffusionTermCylindrical2D(D: FaceVariable) -> csr_array: # extract data from the mesh structure Nx, Ny = D.domain.dims G = D.domain.cell_numbers() - DX = D.domain.cellsize.x - DY = D.domain.cellsize.y + DX = D.domain.cellsize._x + DY = D.domain.cellsize._y dx = 0.5*(DX[0:-1]+DX[1:]) dy = 0.5*(DY[0:-1]+DY[1:]) - rp = D.domain.cellcenters.x - rf = D.domain.facecenters.x[:, np.newaxis] + rp = D.domain.cellcenters._x + rf = D.domain.facecenters._x[:, np.newaxis] mn = Nx*Ny # reassign the east, west for code readability (use broadcasting in Julia) - De = rf[1:Nx+1]*D.xvalue[1:Nx+1, :] / \ + De = rf[1:Nx+1]*D._xvalue[1:Nx+1, :] / \ (rp*dx[1:Nx+1]*DX[1:Nx+1])[:, np.newaxis] - Dw = rf[0:Nx]*D.xvalue[0:Nx, :]/(rp*dx[0:Nx]*DX[1:Nx+1])[:, np.newaxis] - Dn = D.yvalue[:, 1:Ny+1]/(dy[1:Ny+1]*DY[1:Ny+1]) - Ds = D.yvalue[:, 0:Ny]/(dy[0:Ny]*DY[1:Ny+1]) + Dw = rf[0:Nx]*D._xvalue[0:Nx, :]/(rp*dx[0:Nx]*DX[1:Nx+1])[:, np.newaxis] + Dn = D._yvalue[:, 1:Ny+1]/(dy[1:Ny+1]*DY[1:Ny+1]) + Ds = D._yvalue[:, 0:Ny]/(dy[0:Ny]*DY[1:Ny+1]) # calculate the coefficients for the internal cells AE = De.ravel() @@ -176,22 +176,22 @@ def diffusionTermPolar2D(D: FaceVariable) -> csr_array: # extract data from the mesh structure Nx, Ny = D.domain.dims G = D.domain.cell_numbers() - DX = D.domain.cellsize.x - DY = D.domain.cellsize.y + DX = D.domain.cellsize._x + DY = D.domain.cellsize._y dx = 0.5*(DX[0:-1]+DX[1:]) dy = 0.5*(DY[0:-1]+DY[1:]) - rp = D.domain.cellcenters.x[:, np.newaxis] - rf = D.domain.facecenters.x[:, np.newaxis] + rp = D.domain.cellcenters._x[:, np.newaxis] + rf = D.domain.facecenters._x[:, np.newaxis] mn = Nx*Ny # reassign the east, west for code readability (use broadcasting in Julia) - De = rf[1:Nx+1]*D.xvalue[1:Nx+1, :] / \ + De = rf[1:Nx+1]*D._xvalue[1:Nx+1, :] / \ (rp*dx[1:Nx+1][:, np.newaxis]*DX[1:Nx+1][:, np.newaxis]) - Dw = rf[0:Nx]*D.xvalue[0:Nx, :] / \ + Dw = rf[0:Nx]*D._xvalue[0:Nx, :] / \ (rp*dx[0:Nx][:, np.newaxis]*DX[1:Nx+1][:, np.newaxis]) - Dn = D.yvalue[:, 1:Ny+1] / \ + Dn = D._yvalue[:, 1:Ny+1] / \ (rp*rp*dy[1:Ny+1][np.newaxis, :]*DY[1:Ny+1][np.newaxis, :]) - Ds = D.yvalue[:, 0:Ny]/(rp*rp*dy[0:Ny][np.newaxis, :] + Ds = D._yvalue[:, 0:Ny]/(rp*rp*dy[0:Ny][np.newaxis, :] * DY[1:Ny+1][np.newaxis, :]) # calculate the coefficients for the internal cells @@ -229,9 +229,9 @@ def diffusionTerm3D(D: FaceVariable) -> csr_array: # extract data from the mesh structure Nx, Ny, Nz = D.domain.dims G = D.domain.cell_numbers() - DX = D.domain.cellsize.x - DY = D.domain.cellsize.y - DZ = D.domain.cellsize.z + DX = D.domain.cellsize._x + DY = D.domain.cellsize._y + DZ = D.domain.cellsize._z dx = 0.5*(DX[0:-1]+DX[1:]) dy = 0.5*(DY[0:-1]+DY[1:]) dz = 0.5*(DZ[0:-1]+DZ[1:]) @@ -241,15 +241,15 @@ def diffusionTerm3D(D: FaceVariable) -> csr_array: # reassign the east, west, north, and south velocity vectors for the # code readability (use broadcasting) - De = D.xvalue[1:Nx+1, :, :] / \ + De = D._xvalue[1:Nx+1, :, :] / \ (dx[1:Nx+1]*DX[1:Nx+1])[:, np.newaxis, np.newaxis] - Dw = D.xvalue[0:Nx, :, :]/(dx[0:Nx]*DX[1:Nx+1])[:, np.newaxis, np.newaxis] - Dn = D.yvalue[:, 1:Ny+1, :] / \ + Dw = D._xvalue[0:Nx, :, :]/(dx[0:Nx]*DX[1:Nx+1])[:, np.newaxis, np.newaxis] + Dn = D._yvalue[:, 1:Ny+1, :] / \ (dy[1:Ny+1]*DY[1:Ny+1])[np.newaxis, :, np.newaxis] - Ds = D.yvalue[:, 0:Ny, :]/(dy[0:Ny]*DY[1:Ny+1])[np.newaxis, :, np.newaxis] - Df = D.zvalue[:, :, 1:Nz+1] / \ + Ds = D._yvalue[:, 0:Ny, :]/(dy[0:Ny]*DY[1:Ny+1])[np.newaxis, :, np.newaxis] + Df = D._zvalue[:, :, 1:Nz+1] / \ (dz[1:Nz+1]*DZ[1:Nz+1])[np.newaxis, np.newaxis, :] - Db = D.zvalue[:, :, 0:Nz]/(dz[0:Nz]*DZ[1:Nz+1])[np.newaxis, np.newaxis, :] + Db = D._zvalue[:, :, 0:Nz]/(dz[0:Nz]*DZ[1:Nz+1])[np.newaxis, np.newaxis, :] # calculate the coefficients for the internal cells AE = De.ravel() @@ -292,33 +292,33 @@ def diffusionTermCylindrical3D(D: FaceVariable) -> csr_array: # extract data from the mesh structure Nx, Ny, Nz = D.domain.dims G = D.domain.cell_numbers() - DX = D.domain.cellsize.x - DY = D.domain.cellsize.y - DZ = D.domain.cellsize.z + DX = D.domain.cellsize._x + DY = D.domain.cellsize._y + DZ = D.domain.cellsize._z dx = 0.5*(DX[0:-1]+DX[1:]) dy = 0.5*(DY[0:-1]+DY[1:]) dz = 0.5*(DZ[0:-1]+DZ[1:]) - rp = D.domain.cellcenters.x[:, np.newaxis, np.newaxis] - rf = D.domain.facecenters.x[:, np.newaxis, np.newaxis] + rp = D.domain.cellcenters._x[:, np.newaxis, np.newaxis] + rf = D.domain.facecenters._x[:, np.newaxis, np.newaxis] # define the vectors to store the sparse matrix data mn = Nx*Ny*Nz # reassign the east, west, north, and south velocity vectors for the # code readability (use broadcasting) - De = rf[1:Nx+1]*D.xvalue[1:Nx+1, :, :] / \ + De = rf[1:Nx+1]*D._xvalue[1:Nx+1, :, :] / \ (rp*dx[1:Nx+1][:, np.newaxis, np.newaxis] * DX[1:Nx+1][:, np.newaxis, np.newaxis]) - Dw = rf[0:Nx]*D.xvalue[0:Nx, :, :] / \ + Dw = rf[0:Nx]*D._xvalue[0:Nx, :, :] / \ (rp*dx[0:Nx][:, np.newaxis, np.newaxis] * DX[1:Nx+1][:, np.newaxis, np.newaxis]) - Dn = D.yvalue[:, 1:Ny+1, :]/(rp*rp*dy[1:Ny+1][np.newaxis, + Dn = D._yvalue[:, 1:Ny+1, :]/(rp*rp*dy[1:Ny+1][np.newaxis, :, np.newaxis]*DY[1:Ny+1][np.newaxis, :, np.newaxis]) - Ds = D.yvalue[:, 0:Ny, :]/(rp*rp*dy[0:Ny][np.newaxis, + Ds = D._yvalue[:, 0:Ny, :]/(rp*rp*dy[0:Ny][np.newaxis, :, np.newaxis]*DY[1:Ny+1][np.newaxis, :, np.newaxis]) - Df = D.zvalue[:, :, 1:Nz+1] / \ + Df = D._zvalue[:, :, 1:Nz+1] / \ (dz[1:Nz+1]*DZ[1:Nz+1])[np.newaxis, np.newaxis, :] - Db = D.zvalue[:, :, 0:Nz]/(dz[0:Nz]*DZ[1:Nz+1])[np.newaxis, np.newaxis, :] + Db = D._zvalue[:, :, 0:Nz]/(dz[0:Nz]*DZ[1:Nz+1])[np.newaxis, np.newaxis, :] # calculate the coefficients for the internal cells AE = De.ravel() diff --git a/src/pyfvtool/face.py b/src/pyfvtool/face.py index f574550..d6c42bb 100644 --- a/src/pyfvtool/face.py +++ b/src/pyfvtool/face.py @@ -5,6 +5,7 @@ from .mesh import Grid1D, Grid2D, Grid3D from .mesh import CylindricalGrid1D, CylindricalGrid2D from .mesh import PolarGrid2D, CylindricalGrid3D +from .mesh import SphericalGrid3D @@ -31,9 +32,9 @@ def __init__(self, mesh: MeshStructure, faceval : np.ndarray): @overload def __init__(self, mesh: MeshStructure, - xvalue: np.ndarray, - yvalue: np.ndarray, - zvalue: np.ndarray): + _xvalue: np.ndarray, + _yvalue: np.ndarray, + _zvalue: np.ndarray): ... @@ -41,224 +42,438 @@ def __init__(self, mesh: MeshStructure, *args): if len(args)==3: - xvalue = args[0] - yvalue = args[1] - zvalue = args[2] + _xvalue = args[0] + _yvalue = args[1] + _zvalue = args[2] elif len(args)==1: faceval = args[0] if issubclass(type(mesh), Grid1D): Nx = mesh.dims if np.isscalar(faceval): - xvalue = faceval*np.ones(Nx+1) - yvalue = np.array([]) - zvalue = np.array([]) + _xvalue = faceval*np.ones(Nx+1) + _yvalue = np.array([]) + _zvalue = np.array([]) else: - xvalue = faceval[0]*np.ones(Nx+1) - yvalue = np.array([]) - zvalue = np.array([]) + _xvalue = faceval[0]*np.ones(Nx+1) + _yvalue = np.array([]) + _zvalue = np.array([]) elif issubclass(type(mesh), Grid2D): Nx, Ny = mesh.dims if np.isscalar(faceval): - xvalue = faceval*np.ones((Nx+1, Ny)) - yvalue = faceval*np.ones((Nx, Ny+1)) - zvalue = np.array([]) + _xvalue = faceval*np.ones((Nx+1, Ny)) + _yvalue = faceval*np.ones((Nx, Ny+1)) + _zvalue = np.array([]) else: - xvalue = faceval[0]*np.ones((Nx+1, Ny)) - yvalue = faceval[1]*np.ones((Nx, Ny+1)) - zvalue = np.array([]) + _xvalue = faceval[0]*np.ones((Nx+1, Ny)) + _yvalue = faceval[1]*np.ones((Nx, Ny+1)) + _zvalue = np.array([]) elif issubclass(type(mesh), Grid3D): Nx, Ny, Nz = mesh.dims if np.isscalar(faceval): - xvalue = faceval*np.ones((Nx+1, Ny, Nz)) - yvalue = faceval*np.ones((Nx, Ny+1, Nz)) - zvalue = faceval*np.ones((Nx, Ny, Nz+1)) + _xvalue = faceval*np.ones((Nx+1, Ny, Nz)) + _yvalue = faceval*np.ones((Nx, Ny+1, Nz)) + _zvalue = faceval*np.ones((Nx, Ny, Nz+1)) else: - xvalue = faceval[0]*np.ones((Nx+1, Ny, Nz)) - yvalue = faceval[1]*np.ones((Nx, Ny+1, Nz)) - zvalue = faceval[2]*np.ones((Nx, Ny, Nz+1)) + _xvalue = faceval[0]*np.ones((Nx+1, Ny, Nz)) + _yvalue = faceval[1]*np.ones((Nx, Ny+1, Nz)) + _zvalue = faceval[2]*np.ones((Nx, Ny, Nz+1)) else: raise TypeError('Unexpected number of arguments') - self.domain = mesh - self.xvalue = xvalue - self.yvalue = yvalue - self.zvalue = zvalue - + self._xvalue = _xvalue + self._yvalue = _yvalue + self._zvalue = _zvalue + + + @property + def xvalue(self): + if (type(self.domain) is Grid1D)\ + or (type(self.domain) is Grid2D)\ + or (type(self.domain) is Grid3D): + return self._xvalue + elif (type(self.domain) is CylindricalGrid1D)\ + or (type(self.domain) is CylindricalGrid2D)\ + or (type(self.domain) is CylindricalGrid3D): + raise AttributeError('xvalue does not exist for cylindrical grids') + elif (type(self.domain) is PolarGrid2D)\ + or (type(self.domain) is SphericalGrid3D): + raise AttributeError('xvalue does not exist for polar or spherical grids') + else: + raise NotImplementedError("FaceVariable not implemented for mesh type '{0:s}'".\ + format(self.domain.__class__.__name__)) + + @xvalue.setter + def xvalue(self, value): + if (type(self.domain) is Grid1D)\ + or (type(self.domain) is Grid2D)\ + or (type(self.domain) is Grid3D): + self._xvalue = value + elif (type(self.domain) is CylindricalGrid1D)\ + or (type(self.domain) is CylindricalGrid2D)\ + or (type(self.domain) is CylindricalGrid3D): + raise AttributeError('xvalue does not exist for cylindrical grids') + elif (type(self.domain) is PolarGrid2D)\ + or (type(self.domain) is SphericalGrid3D): + raise AttributeError('xvalue does not exist for polar or spherical grids') + else: + raise NotImplementedError("FaceVariable not implemented for mesh type '{0:s}'".\ + format(self.domain.__class__.__name__)) + + + @property + def yvalue(self): + if (type(self.domain) is Grid2D)\ + or (type(self.domain) is Grid3D): + return self._yvalue + elif (type(self.domain) is Grid1D): + raise AttributeError('yvalue does not exist for 1D Cartesian grids') + elif (type(self.domain) is CylindricalGrid1D)\ + or (type(self.domain) is CylindricalGrid2D)\ + or (type(self.domain) is CylindricalGrid3D): + raise AttributeError('yvalue does not exist for cylindrical grids') + elif (type(self.domain) is PolarGrid2D)\ + or (type(self.domain) is SphericalGrid3D): + raise AttributeError('yvalue does not exist for polar or spherical grids') + else: + raise NotImplementedError("FaceVariable not implemented for mesh type '{0:s}'".\ + format(self.domain.__class__.__name__)) + + @yvalue.setter + def yvalue(self, value): + if (type(self.domain) is Grid2D)\ + or (type(self.domain) is Grid3D): + self._yvalue = value + elif (type(self.domain) is Grid1D): + raise AttributeError('yvalue does not exist for 1D Cartesian grids') + elif (type(self.domain) is CylindricalGrid1D)\ + or (type(self.domain) is CylindricalGrid2D)\ + or (type(self.domain) is CylindricalGrid3D): + raise AttributeError('yvalue does not exist for cylindrical grids') + elif (type(self.domain) is PolarGrid2D)\ + or (type(self.domain) is SphericalGrid3D): + raise AttributeError('yvalue does not exist for polar or spherical grids') + else: + raise NotImplementedError("FaceVariable not implemented for mesh type '{0:s}'".\ + format(self.domain.__class__.__name__)) + + @property + def zvalue(self): + if (type(self.domain) is Grid3D)\ + or (type(self.domain) is CylindricalGrid3D): + return self._zvalue + elif (type(self.domain) is Grid2D)\ + or (type(self.domain) is Grid1D): + raise AttributeError('zvalue does not exist for 1D and 2D Cartesian grids') + elif (type(self.domain) is CylindricalGrid1D): + raise AttributeError('zvalue does not exist for 1D cylindrical grids') + elif (type(self.domain) is CylindricalGrid2D): + return self._yvalue + elif (type(self.domain) is PolarGrid2D)\ + or (type(self.domain) is SphericalGrid3D): + raise AttributeError('zvalue does not exist for polar or spherical grids') + else: + raise NotImplementedError("FaceVariable not implemented for mesh type '{0:s}'".\ + format(self.domain.__class__.__name__)) + + @zvalue.setter + def zvalue(self, value): + if (type(self.domain) is Grid3D)\ + or (type(self.domain) is CylindricalGrid3D): + self._zvalue = value + elif (type(self.domain) is Grid2D)\ + or (type(self.domain) is Grid1D): + raise AttributeError('zvalue does not exist for 1D and 2D Cartesian grids') + elif (type(self.domain) is CylindricalGrid1D): + raise AttributeError('zvalue does not exist for 1D cylindrical grids') + elif (type(self.domain) is CylindricalGrid2D): + self._yvalue = value + elif (type(self.domain) is PolarGrid2D)\ + or (type(self.domain) is SphericalGrid3D): + raise AttributeError('zvalue does not exist for polar or spherical grids') + else: + raise NotImplementedError("FaceVariable not implemented for mesh type '{0:s}'".\ + format(self.domain.__class__.__name__)) + + + @property + def rvalue(self): + if (type(self.domain) is Grid1D)\ + or (type(self.domain) is Grid2D)\ + or (type(self.domain) is Grid3D): + raise AttributeError('rvalue does not exist for Cartesian grids') + elif (type(self.domain) is CylindricalGrid1D)\ + or (type(self.domain) is CylindricalGrid2D)\ + or (type(self.domain) is CylindricalGrid3D)\ + or (type(self.domain) is PolarGrid2D)\ + or (type(self.domain) is SphericalGrid3D): + return self._xvalue + else: + raise NotImplementedError("FaceVariable not implemented for mesh type '{0:s}'".\ + format(self.domain.__class__.__name__)) + + @rvalue.setter + def rvalue(self, value): + if (type(self.domain) is Grid1D)\ + or (type(self.domain) is Grid2D)\ + or (type(self.domain) is Grid3D): + self._xvalue = value + elif (type(self.domain) is CylindricalGrid1D)\ + or (type(self.domain) is CylindricalGrid2D)\ + or (type(self.domain) is CylindricalGrid3D)\ + or (type(self.domain) is PolarGrid2D)\ + or (type(self.domain) is SphericalGrid3D): + self._xvalue = value + else: + raise NotImplementedError("FaceVariable not implemented for mesh type '{0:s}'".\ + format(self.domain.__class__.__name__)) + + + @property + def thetavalue(self): + if (type(self.domain) is PolarGrid2D)\ + or (type(self.domain) is CylindricalGrid3D)\ + or (type(self.domain) is SphericalGrid3D): + return self._yvalue + elif (type(self.domain) is Grid1D)\ + or (type(self.domain) is Grid2D)\ + or (type(self.domain) is Grid3D): + raise AttributeError('thetavalue does not exist for Cartesian grids') + elif (type(self.domain) is CylindricalGrid1D)\ + or (type(self.domain) is CylindricalGrid2D): + raise AttributeError('thetavalue does not exist for 1D and 2D cylindrical grids') + else: + raise NotImplementedError("FaceVariable not implemented for mesh type '{0:s}'".\ + format(self.domain.__class__.__name__)) + + @thetavalue.setter + def thetavalue(self, value): + if (type(self.domain) is PolarGrid2D)\ + or (type(self.domain) is CylindricalGrid3D)\ + or (type(self.domain) is SphericalGrid3D): + self._yvalue = value + elif (type(self.domain) is Grid1D)\ + or (type(self.domain) is Grid2D)\ + or (type(self.domain) is Grid3D): + raise AttributeError('thetavalue does not exist for Cartesian grids') + elif (type(self.domain) is CylindricalGrid1D)\ + or (type(self.domain) is CylindricalGrid2D): + raise AttributeError('thetavalue does not exist for 1D and 2D cylindrical grids') + else: + raise NotImplementedError("FaceVariable not implemented for mesh type '{0:s}'".\ + format(self.domain.__class__.__name__)) + + + @property + def phivalue(self): + if (type(self.domain) is SphericalGrid3D): + return self._zvalue + elif (type(self.domain) is Grid1D)\ + or (type(self.domain) is Grid2D)\ + or (type(self.domain) is Grid3D): + raise AttributeError('phivalue does not exist for Cartesian grids') + elif (type(self.domain) is CylindricalGrid1D)\ + or (type(self.domain) is CylindricalGrid2D)\ + or (type(self.domain) is PolarGrid2D)\ + or (type(self.domain) is CylindricalGrid3D): + raise AttributeError('phivalue does not exist for cylindrical and polar grids') + else: + raise NotImplementedError("FaceVariable not implemented for mesh type '{0:s}'".\ + format(self.domain.__class__.__name__)) + + @phivalue.setter + def phivalue(self, value): + if (type(self.domain) is SphericalGrid3D): + self._zvalue = value + elif (type(self.domain) is Grid1D)\ + or (type(self.domain) is Grid2D)\ + or (type(self.domain) is Grid3D): + raise AttributeError('phivalue does not exist for Cartesian grids') + elif (type(self.domain) is CylindricalGrid1D)\ + or (type(self.domain) is CylindricalGrid2D)\ + or (type(self.domain) is PolarGrid2D)\ + or (type(self.domain) is CylindricalGrid3D): + raise AttributeError('phivalue does not exist for cylindrical and polar grids') + else: + raise NotImplementedError("FaceVariable not implemented for mesh type '{0:s}'".\ + format(self.domain.__class__.__name__)) + + + def __add__(self, other): if type(other) is FaceVariable: return FaceVariable(self.domain, - self.xvalue+other.xvalue, - self.yvalue+other.yvalue, - self.zvalue+other.zvalue) + self._xvalue+other._xvalue, + self._yvalue+other._yvalue, + self._zvalue+other._zvalue) else: return FaceVariable(self.domain, - self.xvalue+other, - self.yvalue+other, - self.zvalue+other) + self._xvalue+other, + self._yvalue+other, + self._zvalue+other) def __radd__(self, other): if type(other) is FaceVariable: return FaceVariable(self.domain, - self.xvalue+other.xvalue, - self.yvalue+other.yvalue, - self.zvalue+other.zvalue) + self._xvalue+other._xvalue, + self._yvalue+other._yvalue, + self._zvalue+other._zvalue) else: return FaceVariable(self.domain, - self.xvalue+other, - self.yvalue+other, - self.zvalue+other) + self._xvalue+other, + self._yvalue+other, + self._zvalue+other) def __rsub__(self, other): if type(other) is FaceVariable: - return FaceVariable(self.domain, other.xvalue-self.xvalue, - other.yvalue-self.yvalue, - other.zvalue-self.zvalue) + return FaceVariable(self.domain, other._xvalue-self._xvalue, + other._yvalue-self._yvalue, + other._zvalue-self._zvalue) else: return FaceVariable(self.domain, - other-self.xvalue, - other-self.yvalue, - other-self.zvalue) + other-self._xvalue, + other-self._yvalue, + other-self._zvalue) def __sub__(self, other): if type(other) is FaceVariable: - return FaceVariable(self.domain, self.xvalue-other.xvalue, - self.yvalue-other.yvalue, - self.zvalue-other.zvalue) + return FaceVariable(self.domain, self._xvalue-other._xvalue, + self._yvalue-other._yvalue, + self._zvalue-other._zvalue) else: - return FaceVariable(self.domain, self.xvalue-other, - self.yvalue-other, - self.zvalue-other) + return FaceVariable(self.domain, self._xvalue-other, + self._yvalue-other, + self._zvalue-other) def __mul__(self, other): if type(other) is FaceVariable: - return FaceVariable(self.domain, self.xvalue*other.xvalue, - self.yvalue*other.yvalue, - self.zvalue*other.zvalue) + return FaceVariable(self.domain, self._xvalue*other._xvalue, + self._yvalue*other._yvalue, + self._zvalue*other._zvalue) else: - return FaceVariable(self.domain, self.xvalue*other, - self.yvalue*other, - self.zvalue*other) + return FaceVariable(self.domain, self._xvalue*other, + self._yvalue*other, + self._zvalue*other) def __rmul__(self, other): if type(other) is FaceVariable: - return FaceVariable(self.domain, self.xvalue*other.xvalue, - self.yvalue*other.yvalue, - self.zvalue*other.zvalue) + return FaceVariable(self.domain, self._xvalue*other._xvalue, + self._yvalue*other._yvalue, + self._zvalue*other._zvalue) else: - return FaceVariable(self.domain, self.xvalue*other, - self.yvalue*other, - self.zvalue*other) + return FaceVariable(self.domain, self._xvalue*other, + self._yvalue*other, + self._zvalue*other) def __truediv__(self, other): if type(other) is FaceVariable: - return FaceVariable(self.domain, self.xvalue/other.xvalue, - self.yvalue/other.yvalue, - self.zvalue/other.zvalue) + return FaceVariable(self.domain, self._xvalue/other._xvalue, + self._yvalue/other._yvalue, + self._zvalue/other._zvalue) else: - return FaceVariable(self.domain, self.xvalue/other, - self.yvalue/other, - self.zvalue/other) + return FaceVariable(self.domain, self._xvalue/other, + self._yvalue/other, + self._zvalue/other) def __rtruediv__(self, other): if type(other) is FaceVariable: - return FaceVariable(self.domain, other.xvalue/self.xvalue, - other.yvalue/self.yvalue, - other.zvalue/self.zvalue) + return FaceVariable(self.domain, other._xvalue/self._xvalue, + other._yvalue/self._yvalue, + other._zvalue/self._zvalue) else: - return FaceVariable(self.domain, other/self.xvalue, - other/self.yvalue, - other/self.zvalue) + return FaceVariable(self.domain, other/self._xvalue, + other/self._yvalue, + other/self._zvalue) def __neg__(self): - return FaceVariable(self.domain, -self.xvalue, - -self.yvalue, - -self.zvalue) + return FaceVariable(self.domain, -self._xvalue, + -self._yvalue, + -self._zvalue) def __pow__(self, other): if type(other) is FaceVariable: - return FaceVariable(self.domain, self.xvalue**other.xvalue, - self.yvalue**other.yvalue, - self.zvalue**other.zvalue) + return FaceVariable(self.domain, self._xvalue**other._xvalue, + self._yvalue**other._yvalue, + self._zvalue**other._zvalue) else: - return FaceVariable(self.domain, self.xvalue**other, - self.yvalue**other, - self.zvalue**other) + return FaceVariable(self.domain, self._xvalue**other, + self._yvalue**other, + self._zvalue**other) def __rpow__(self, other): if type(other) is FaceVariable: - return FaceVariable(self.domain, other.xvalue**self.xvalue, - other.yvalue**self.yvalue, - other.zvalue**self.zvalue) + return FaceVariable(self.domain, other._xvalue**self._xvalue, + other._yvalue**self._yvalue, + other._zvalue**self._zvalue) else: - return FaceVariable(self.domain, other**self.xvalue, - other**self.yvalue, - other**self.zvalue) + return FaceVariable(self.domain, other**self._xvalue, + other**self._yvalue, + other**self._zvalue) def __gt__(self, other): if type(other) is FaceVariable: - return FaceVariable(self.domain, self.xvalue > other.xvalue, - self.yvalue > other.yvalue, - self.zvalue > other.zvalue) + return FaceVariable(self.domain, self._xvalue > other._xvalue, + self._yvalue > other._yvalue, + self._zvalue > other._zvalue) else: - return FaceVariable(self.domain, self.xvalue > other, - self.yvalue > other, - self.zvalue > other) + return FaceVariable(self.domain, self._xvalue > other, + self._yvalue > other, + self._zvalue > other) def __ge__(self, other): if type(other) is FaceVariable: - return FaceVariable(self.domain, self.xvalue >= other.xvalue, - self.yvalue >= other.yvalue, - self.zvalue >= other.zvalue) + return FaceVariable(self.domain, self._xvalue >= other._xvalue, + self._yvalue >= other._yvalue, + self._zvalue >= other._zvalue) else: - return FaceVariable(self.domain, self.xvalue >= other, - self.yvalue >= other, - self.zvalue >= other) + return FaceVariable(self.domain, self._xvalue >= other, + self._yvalue >= other, + self._zvalue >= other) def __lt__(self, other): if type(other) is FaceVariable: - return FaceVariable(self.domain, self.xvalue < other.xvalue, - self.yvalue < other.yvalue, - self.zvalue < other.zvalue) + return FaceVariable(self.domain, self._xvalue < other._xvalue, + self._yvalue < other._yvalue, + self._zvalue < other._zvalue) else: - return FaceVariable(self.domain, self.xvalue < other, - self.yvalue < other, - self.zvalue < other) + return FaceVariable(self.domain, self._xvalue < other, + self._yvalue < other, + self._zvalue < other) def __le__(self, other): if type(other) is FaceVariable: - return FaceVariable(self.domain, self.xvalue <= other.xvalue, - self.yvalue <= other.yvalue, - self.zvalue <= other.zvalue) + return FaceVariable(self.domain, self._xvalue <= other._xvalue, + self._yvalue <= other._yvalue, + self._zvalue <= other._zvalue) else: - return FaceVariable(self.domain, self.xvalue <= other, - self.yvalue <= other, - self.zvalue <= other) + return FaceVariable(self.domain, self._xvalue <= other, + self._yvalue <= other, + self._zvalue <= other) def __and__(self, other): if type(other) is FaceVariable: - return FaceVariable(self.domain, np.logical_and(self.xvalue, other.xvalue), - np.logical_and(self.yvalue, other.yvalue), - np.logical_and(self.zvalue, other.zvalue)) + return FaceVariable(self.domain, np.logical_and(self._xvalue, other._xvalue), + np.logical_and(self._yvalue, other._yvalue), + np.logical_and(self._zvalue, other._zvalue)) else: - return FaceVariable(self.domain, np.logical_and(self.xvalue, other), - np.logical_and(self.yvalue, other), - np.logical_and(self.zvalue, other)) + return FaceVariable(self.domain, np.logical_and(self._xvalue, other), + np.logical_and(self._yvalue, other), + np.logical_and(self._zvalue, other)) def __or__(self, other): if type(other) is FaceVariable: - return FaceVariable(self.domain, np.logical_or(self.xvalue, other.xvalue), - np.logical_or(self.yvalue, other.yvalue), - np.logical_or(self.zvalue, other.zvalue)) + return FaceVariable(self.domain, np.logical_or(self._xvalue, other._xvalue), + np.logical_or(self._yvalue, other._yvalue), + np.logical_or(self._zvalue, other._zvalue)) else: - return FaceVariable(self.domain, np.logical_or(self.xvalue, other), - np.logical_or(self.yvalue, other), - np.logical_or(self.zvalue, other)) + return FaceVariable(self.domain, np.logical_or(self._xvalue, other), + np.logical_or(self._yvalue, other), + np.logical_or(self._zvalue, other)) def __abs__(self): - return FaceVariable(self.domain, np.abs(self.xvalue), - np.abs(self.yvalue), - np.abs(self.zvalue)) + return FaceVariable(self.domain, np.abs(self._xvalue), + np.abs(self._yvalue), + np.abs(self._zvalue)) @@ -302,7 +517,7 @@ def faceLocations(m: MeshStructure): if (type(m) is Grid1D)\ or (type(m) is CylindricalGrid1D): X = FaceVariable(m, 0) - X.xvalue = m.facecenters.x + X._xvalue = m.facecenters._x return X elif (type(m) is Grid2D)\ @@ -310,10 +525,10 @@ def faceLocations(m: MeshStructure): or (type(m) is PolarGrid2D): X = FaceVariable(m, 0) Y = FaceVariable(m, 0) - X.xvalue = np.tile(m.facecenters.x[:, np.newaxis], (1, N[1])) - X.yvalue = np.tile(m.cellcenters.y[:, np.newaxis].T, (N[0]+1, 1)) - Y.xvalue = np.tile(m.cellcenters.x[:, np.newaxis], (1, N[1]+1)) - Y.yvalue = np.tile(m.facecenters.y[:, np.newaxis].T, (N[0], 1)) + X._xvalue = np.tile(m.facecenters._x[:, np.newaxis], (1, N[1])) + X._yvalue = np.tile(m.cellcenters._y[:, np.newaxis].T, (N[0]+1, 1)) + Y._xvalue = np.tile(m.cellcenters._x[:, np.newaxis], (1, N[1]+1)) + Y._yvalue = np.tile(m.facecenters._y[:, np.newaxis].T, (N[0], 1)) return X, Y elif (type(m) is Grid3D)\ @@ -322,21 +537,21 @@ def faceLocations(m: MeshStructure): Y = FaceVariable(m, 0) Z = FaceVariable(m, 0) z = np.zeros((1,1,N[2])) - z[0, 0, :] = m.cellcenters.z + z[0, 0, :] = m.cellcenters._z - X.xvalue = np.tile(m.facecenters.x[:, np.newaxis, np.newaxis], (1, N[1], N[2])) - X.yvalue = np.tile((m.cellcenters.y[:, np.newaxis].T)[:, :, np.newaxis], (N[0]+1, 1, N[2])) - X.zvalue = np.tile(z, (N[0]+1, N[1], 1)) + X._xvalue = np.tile(m.facecenters._x[:, np.newaxis, np.newaxis], (1, N[1], N[2])) + X._yvalue = np.tile((m.cellcenters._y[:, np.newaxis].T)[:, :, np.newaxis], (N[0]+1, 1, N[2])) + X._zvalue = np.tile(z, (N[0]+1, N[1], 1)) - Y.xvalue = np.tile(m.cellcenters.x[:, np.newaxis, np.newaxis], (1, N[1]+1, N[2])) - Y.yvalue = np.tile((m.facecenters.y[:, np.newaxis].T)[:, :, np.newaxis], (N[0], 1, N[2])) - Y.zvalue = np.tile(z, (N[0], N[1]+1, 1)) + Y._xvalue = np.tile(m.cellcenters._x[:, np.newaxis, np.newaxis], (1, N[1]+1, N[2])) + Y._yvalue = np.tile((m.facecenters._y[:, np.newaxis].T)[:, :, np.newaxis], (N[0], 1, N[2])) + Y._zvalue = np.tile(z, (N[0], N[1]+1, 1)) z = np.zeros((1,1,N[2]+1)) - z[0, 0, :] = m.cellcenters.z - Z.xvalue = np.tile(m.cellcenters.x[:, np.newaxis, np.newaxis], (1, N[1], N[2]+1)) - Z.yvalue = np.tile((m.facecenters.y[:, np.newaxis].T)[:, :, np.newaxis], (N[0], 1, N[2]+1)) - Z.zvalue = np.tile(z, (N[0], N[1], 1)) + z[0, 0, :] = m.cellcenters._z + Z._xvalue = np.tile(m.cellcenters._x[:, np.newaxis, np.newaxis], (1, N[1], N[2]+1)) + Z._yvalue = np.tile((m.facecenters._y[:, np.newaxis].T)[:, :, np.newaxis], (N[0], 1, N[2]+1)) + Z._zvalue = np.tile(z, (N[0], N[1], 1)) return X, Y, Z raise TypeError('mesh type not implemented') return None @@ -353,41 +568,41 @@ def faceeval(f, *args): """ if len(args)==1: return FaceVariable(args[0].domain, - f(args[0].xvalue), - f(args[0].yvalue), - f(args[0].zvalue)) + f(args[0]._xvalue), + f(args[0]._yvalue), + f(args[0]._zvalue)) elif len(args)==2: return FaceVariable(args[0].domain, - f(args[0].xvalue, args[1].xvalue), - f(args[0].yvalue, args[1].yvalue), - f(args[0].zvalue, args[1].zvalue)) + f(args[0]._xvalue, args[1]._xvalue), + f(args[0]._yvalue, args[1]._yvalue), + f(args[0]._zvalue, args[1]._zvalue)) elif len(args)==3: return FaceVariable(args[0].domain, - f(args[0].xvalue, args[1].xvalue, args[2].xvalue), - f(args[0].yvalue, args[1].yvalue, args[2].yvalue), - f(args[0].zvalue, args[1].zvalue, args[2].zvalue)) + f(args[0]._xvalue, args[1]._xvalue, args[2]._xvalue), + f(args[0]._yvalue, args[1]._yvalue, args[2]._yvalue), + f(args[0]._zvalue, args[1]._zvalue, args[2]._zvalue)) elif len(args)==4: return FaceVariable(args[0].domain, - f(args[0].xvalue, args[1].xvalue, args[2].xvalue, args[3].xvalue), - f(args[0].yvalue, args[1].yvalue, args[2].yvalue, args[3].yvalue), - f(args[0].zvalue, args[1].zvalue, args[2].zvalue, args[3].zvalue)) + f(args[0]._xvalue, args[1]._xvalue, args[2]._xvalue, args[3]._xvalue), + f(args[0]._yvalue, args[1]._yvalue, args[2]._yvalue, args[3]._yvalue), + f(args[0]._zvalue, args[1]._zvalue, args[2]._zvalue, args[3]._zvalue)) elif len(args)==5: return FaceVariable(args[0].domain, - f(args[0].xvalue, args[1].xvalue, args[2].xvalue, args[3].xvalue, args[4].xvalue), - f(args[0].yvalue, args[1].yvalue, args[2].yvalue, args[3].yvalue, args[4].yvalue), - f(args[0].zvalue, args[1].zvalue, args[2].zvalue, args[3].zvalue, args[4].zvalue)) + f(args[0]._xvalue, args[1]._xvalue, args[2]._xvalue, args[3]._xvalue, args[4]._xvalue), + f(args[0]._yvalue, args[1]._yvalue, args[2]._yvalue, args[3]._yvalue, args[4]._yvalue), + f(args[0]._zvalue, args[1]._zvalue, args[2]._zvalue, args[3]._zvalue, args[4]._zvalue)) elif len(args)==6: return FaceVariable(args[0].domain, - f(args[0].xvalue, args[1].xvalue, args[2].xvalue, args[3].xvalue, args[4].xvalue, args[5].xvalue), - f(args[0].yvalue, args[1].yvalue, args[2].yvalue, args[3].yvalue, args[4].yvalue, args[5].yvalue), - f(args[0].zvalue, args[1].zvalue, args[2].zvalue, args[3].zvalue, args[4].zvalue, args[5].zvalue)) + f(args[0]._xvalue, args[1]._xvalue, args[2]._xvalue, args[3]._xvalue, args[4]._xvalue, args[5]._xvalue), + f(args[0]._yvalue, args[1]._yvalue, args[2]._yvalue, args[3]._yvalue, args[4]._yvalue, args[5]._yvalue), + f(args[0]._zvalue, args[1]._zvalue, args[2]._zvalue, args[3]._zvalue, args[4]._zvalue, args[5]._zvalue)) elif len(args)==7: return FaceVariable(args[0].domain, - f(args[0].xvalue, args[1].xvalue, args[2].xvalue, args[3].xvalue, args[4].xvalue, args[5].xvalue, args[6].xvalue), - f(args[0].yvalue, args[1].yvalue, args[2].yvalue, args[3].yvalue, args[4].yvalue, args[5].yvalue, args[6].yvalue), - f(args[0].zvalue, args[1].zvalue, args[2].zvalue, args[3].zvalue, args[4].zvalue, args[5].zvalue, args[6].zvalue)) + f(args[0]._xvalue, args[1]._xvalue, args[2]._xvalue, args[3]._xvalue, args[4]._xvalue, args[5]._xvalue, args[6]._xvalue), + f(args[0]._yvalue, args[1]._yvalue, args[2]._yvalue, args[3]._yvalue, args[4]._yvalue, args[5]._yvalue, args[6]._yvalue), + f(args[0]._zvalue, args[1]._zvalue, args[2]._zvalue, args[3]._zvalue, args[4]._zvalue, args[5]._zvalue, args[6]._zvalue)) elif len(args)==8: return FaceVariable(args[0].domain, - f(args[0].xvalue, args[1].xvalue, args[2].xvalue, args[3].xvalue, args[4].xvalue, args[5].xvalue, args[6].xvalue, args[7].xvalue), - f(args[0].yvalue, args[1].yvalue, args[2].yvalue, args[3].yvalue, args[4].yvalue, args[5].yvalue, args[6].yvalue, args[7].yvalue), - f(args[0].zvalue, args[1].zvalue, args[2].zvalue, args[3].zvalue, args[4].zvalue, args[5].zvalue, args[6].zvalue, args[7].zvalue)) + f(args[0]._xvalue, args[1]._xvalue, args[2]._xvalue, args[3]._xvalue, args[4]._xvalue, args[5]._xvalue, args[6]._xvalue, args[7]._xvalue), + f(args[0]._yvalue, args[1]._yvalue, args[2]._yvalue, args[3]._yvalue, args[4]._yvalue, args[5]._yvalue, args[6]._yvalue, args[7]._yvalue), + f(args[0]._zvalue, args[1]._zvalue, args[2]._zvalue, args[3]._zvalue, args[4]._zvalue, args[5]._zvalue, args[6]._zvalue, args[7]._zvalue)) diff --git a/src/pyfvtool/legacy.py b/src/pyfvtool/legacy.py index 97b9c2e..6a1dda6 100644 --- a/src/pyfvtool/legacy.py +++ b/src/pyfvtool/legacy.py @@ -10,6 +10,7 @@ from .mesh import Grid1D, CylindricalGrid1D, SphericalGrid1D from .mesh import Grid2D, CylindricalGrid2D, PolarGrid2D from .mesh import Grid3D, CylindricalGrid3D, SphericalGrid3D +from .mesh import MeshStructure def boundaryConditionTerm(BC): @@ -28,6 +29,7 @@ def createCellVariable(*args, **kwargs): def createFaceVariable(mesh, faceval): + """Legacy factory function for FaceVariable""" return FaceVariable(mesh, faceval) @@ -74,3 +76,25 @@ def createMeshCylindrical3D(*args) -> CylindricalGrid3D: def createMeshSpherical3D(*args) -> SphericalGrid3D: """Legacy factory function for SphericalGrid3D""" return SphericalGrid3D(*args) + + +def get_CellVariable_profile1D(phi: CellVariable): + """Legacy function""" + return phi.plotprofile() + +def get_CellVariable_profile2D(phi: CellVariable): + return phi.plotprofile() + +def get_CellVariable_profile3D(phi: CellVariable): + return phi.plotprofile() + +def domainInt(phi: CellVariable) -> float: + return phi.domainIntegral() + + +def cellVolume(m: MeshStructure): + BC = BoundaryConditions(m) + c = m.cellVolumes() + return CellVariable(m, c, BC) + + diff --git a/src/pyfvtool/mesh.py b/src/pyfvtool/mesh.py index 23418dc..c5f7ffa 100644 --- a/src/pyfvtool/mesh.py +++ b/src/pyfvtool/mesh.py @@ -9,61 +9,168 @@ #%% # General data structures for describing meshes -class CellSize: - def __init__(self, x: np.ndarray, y: np.ndarray, z: np.ndarray): - self.x = x - self.y = y - self.z = z + +class CellProp: + def __init__(self, _x: np.ndarray, _y: np.ndarray, _z: np.ndarray, + coordlabels: dict): + self._x = _x + self._y = _y + self._z = _z + self.coordlabels = coordlabels def __str__(self): temp = vars(self) + result = "" for item in temp: - print(item, ':', temp[item]) - return "" - + result += f"{item}: {temp[item]}\n" + return result + def __repr__(self): - temp = vars(self) - for item in temp: - print(item, ':', temp[item]) - return "" + return str(self) + + + # The following coordinate-labeling properties should probably be read-only. + # For this reason, the 'setters' have been commented out, but are kept for + # future reference. + + @property + def x(self): + if 'x' in self.coordlabels: + if self.coordlabels['x']=='_x': + return self._x + else: + raise AttributeError("Unexpectedly, user label 'x' does not correspond to '_x'") + else: + raise AttributeError("This mesh has no coordinate labeled 'x'.") + + # @x.setter + # def x(self, value): + # if 'x' in self.coordlabels: + # if self.coordlabels['x']=='_x': + # self._x = value + # else: + # raise AttributeError("Unexpectedly, user label 'x' does not correspond to '_x'") + # else: + # raise AttributeError("This mesh has no coordinate labeled 'x'.") + + @property + def r(self): + if 'r' in self.coordlabels: + if self.coordlabels['r']=='_x': + return self._x + else: + raise AttributeError("Unexpectedly, user label 'r' does not correspond to '_x'") + else: + raise AttributeError("This mesh has no coordinate labeled 'r'.") + + # @r.setter + # def r(self, value): + # if 'r' in self.coordlabels: + # if self.coordlabels['r']=='_x': + # self._x = value + # else: + # raise AttributeError("Unexpectedly, user label 'r' does not correspond to '_x'") + # else: + # raise AttributeError("This mesh has no coordinate labeled 'r'.") + + @property + def z(self): + if 'z' in self.coordlabels: + if self.coordlabels['z']=='_y': + return self._y + elif self.coordlabels['z']=='_z': + return self._z + else: + raise AttributeError(f"Unexpected label correspondence: 'z' -> '{self.coordlabels['z']}'") + else: + raise AttributeError("This mesh has no coordinate labeled 'z'.") + + # @z.setter + # def z(self, value): + # if 'z' in self.coordlabels: + # if self.coordlabels['z']=='_y': + # self._y = value + # if self.coordlabels['z']=='_z': + # self._z = value + # else: + # raise AttributeError(f"Unexpected label correspondence: 'z' -> '{self.coordlabels['z']}'") + # else: + # raise AttributeError("This mesh has no coordinate labeled 'z'.") + + @property + def y(self): + if 'y' in self.coordlabels: + if self.coordlabels['y']=='_y': + return self._y + else: + raise AttributeError(f"Unexpected label correspondence: 'y' -> '{self.coordlabels['y']}'") + else: + raise AttributeError("This mesh has no coordinate labeled 'y'.") + + # @y.setter + # def y(self, value): + # if 'y' in self.coordlabels: + # if self.coordlabels['y']=='_y': + # self._y = value + # else: + # raise AttributeError(f"Unexpected label correspondence: 'y' -> '{self.coordlabels['y']}'") + # else: + # raise AttributeError("This mesh has no coordinate labeled 'y'.") + + @property + def theta(self): + if 'theta' in self.coordlabels: + if self.coordlabels['theta']=='_y': + return self._y + else: + raise AttributeError(f"Unexpected label correspondence: 'theta' -> '{self.coordlabels['y']}'") + else: + raise AttributeError("This mesh has no coordinate labeled 'theta'.") + + # @theta.setter + # def theta(self, value): + # if 'theta' in self.coordlabels: + # if self.coordlabels['theta']=='_y': + # self._y = value + # else: + # raise AttributeError(f"Unexpected label correspondence: 'theta' -> '{self.coordlabels['y']}'") + # else: + # raise AttributeError("This mesh has no coordinate labeled 'theta'.") + + @property + def phi(self): + if 'phi' in self.coordlabels: + if self.coordlabels['phi']=='_z': + return self._z + else: + raise AttributeError(f"Unexpected label correspondence: 'phi' -> '{self.coordlabels['y']}'") + else: + raise AttributeError("This mesh has no coordinate labeled 'phi'.") + + # @phi.setter + # def phi(self, value): + # if 'phi' in self.coordlabels: + # if self.coordlabels['phi']=='_z': + # self._z = value + # else: + # raise AttributeError(f"Unexpected label correspondence: 'phi' -> '{self.coordlabels['y']}'") + # else: + # raise AttributeError("This mesh has no coordinate labeled 'phi'.") -class CellLocation: - def __init__(self, x: np.ndarray, y: np.ndarray, z: np.ndarray): - self.x = x - self.y = y - self.z = z - def __str__(self): - temp = vars(self) - for item in temp: - print(item, ':', temp[item]) - return "" - def __repr__(self): - temp = vars(self) - for item in temp: - print(item, ':', temp[item]) - return "" +class CellSize(CellProp): + pass -class FaceLocation: - def __init__(self, x: np.ndarray, y: np.ndarray, z: np.ndarray): - self.x = x - self.y = y - self.z = z +class CellLocation(CellProp): + pass - def __str__(self): - temp = vars(self) - for item in temp: - print(item, ':', temp[item]) - return "" - def __repr__(self): - temp = vars(self) - for item in temp: - print(item, ':', temp[item]) - return "" +class FaceLocation(CellProp): + pass + class MeshStructure: @@ -78,15 +185,10 @@ def __init__(self, dims, cellsize, def __str__(self): temp = vars(self) + result = "" for item in temp: - print(item, ':', temp[item]) - return "" - - def __repr__(self): - temp = vars(self) - for item in temp: - print(item, ':', temp[item]) - return "" + result += f"{item}: {temp[item]}\n" + return result def _facelocation_to_cellsize(self, facelocation): return np.hstack([facelocation[1]-facelocation[0], @@ -96,13 +198,53 @@ def _facelocation_to_cellsize(self, facelocation): def visualize(self): pass - def shift_origin(self, x=0.0, y=0.0, z=0.0): - self.cellcenters.x += x - self.cellcenters.y += y - self.cellcenters.z += z - self.facecenters.x += x - self.facecenters.y += y - self.facecenters.z += z + + def shift_origin(self, _x=0.0, _y=0.0, _z=0.0): + self.cellcenters._x += _x + self.cellcenters._y += _y + self.cellcenters._z += _z + self.facecenters._x += _x + self.facecenters._y += _y + self.facecenters._z += _z + + + def cellVolumes(self): + """Get the volumes of all finite volume cells in the mesh + + Returns + ------- + np.ndarray + containing all cell volumes, arranged according to gridcells + + TODO: these could perhaps be calculated statically, when initializing + the mesh. + + """ + if (type(self) is Grid1D): + c = self.cellsize._x[1:-1] + elif (type(self) is CylindricalGrid1D): + c = 2.0*np.pi*self.cellsize._x[1:-1]*self.cellcenters._x + elif (type(self) is Grid2D): + c = self.cellsize._x[1:-1][:, np.newaxis]\ + *self.cellsize._y[1:-1][np.newaxis, :] + elif (type(self) is CylindricalGrid2D): + c = 2.0*np.pi*self.cellcenters._x[:, np.newaxis]\ + *self.cellsize._x[1:-1][:, np.newaxis]\ + *self.cellsize._y[1:-1][np.newaxis, :] + elif (type(self) is PolarGrid2D): + c = self.cellcenters._x\ + *self.cellsize._x[1:-1][:, np.newaxis]\ + *self.cellsize._y[1:-1][np.newaxis, :] + elif (type(self) is Grid3D): + c = self.cellsize._x[1:-1][:,np.newaxis,np.newaxis]\ + *self.cellsize._y[1:-1][np.newaxis,:,np.newaxis]\ + *self.cellsize._z[1:-1][np.newaxis,np.newaxis,:] + elif (type(self) is CylindricalGrid3D): + c = self.cellcenters._x\ + *self.cellsize._x[1:-1][:,np.newaxis,np.newaxis]\ + *self.cellsize._y[1:-1][np.newaxis,:,np.newaxis]\ + *self.cellsize._z[np.newaxis,np.newaxis,:] + return c @@ -112,106 +254,51 @@ def shift_origin(self, x=0.0, y=0.0, z=0.0): class Grid1D(MeshStructure): """Mesh based on a 1D Cartesian grid (x) + ===================================== - """ + This class can be instantiated in different ways: from a list of cell face + locations or from the number of cells and domain length. - @overload - def __init__(self, Nx: int, Lx: float): - """ - TODO: docstring (multiple dispatch) - These @overload docstrings do NOT show up in help(pf.Grid1D). Therefore, - put all versions in the main __init__ docstring - - Parameters - ---------- + Instantiation Options: + ---------------------- + - Grid1D(Nx, Lx) + - Grid1D(face_locationsX) + + + Parameters + ---------- + Grid1D(Nx, Lx) Nx : int - DESCRIPTION. + Number of cells in the x direction. Lx : float - DESCRIPTION. - - Returns - ------- - None. + Length of the domain in the x direction. + + Grid1D(face_locationsX) + face_locationsX : ndarray + Locations of the cell faces in the x direction. + + Examples + -------- + >>> import numpy as np + >>> from pyfvtool import Grid1D + >>> mesh = Grid1D(10, 10.0) + >>> print(mesh) + """ - """ + @overload + def __init__(self, Nx: int, Lx: float): ... @overload def __init__(self, face_locations: np.ndarray): - """ - TODO: docstring (multiple dispatch) - These @overload docstrings do NOT show up in help(pf.Grid1D). Therefore, - put all versions in the main __init__ docstring - - - Parameters - ---------- - face_locations : np.ndarray - DESCRIPTION. - - Returns - ------- - None. - - """ ... @overload def __init__(self, dims, cellsize, cellcenters, facecenters, corners, edges): - """ - TODO: docstring (multiple dispatch) - These @overload docstrings do NOT show up in help(pf.Grid1D). Therefore, - put all versions in the main __init__ docstring - - - Parameters - ---------- - dims : TYPE - DESCRIPTION. - cellsize : TYPE - DESCRIPTION. - cellcenters : TYPE - DESCRIPTION. - facecenters : TYPE - DESCRIPTION. - corners : TYPE - DESCRIPTION. - edges : TYPE - DESCRIPTION. - - Returns - ------- - None. - - """ ... def __init__(self, *args): - """Create a Grid1D mesh object from a list of cell face locations or from - number of cells and domain length. - - Parameters - ---------- - Nx : int - Number of cells in the x direction. - Lx : float - Length of the domain in the x direction. - face_locations : ndarray - Locations of the cell faces in the x direction. - - Returns - ------- - Grid1D - A 1D mesh object. - - Examples - -------- - >>> import numpy as np - >>> from pyfvtool import Grid1D - >>> mesh = Grid1D(10, 10.0) - >>> print(mesh) - """ if (len(args)==6): dims, cell_size, cell_location, face_location, corners, edges\ = args @@ -221,15 +308,7 @@ def __init__(self, *args): super().__init__(dims, cell_size, cell_location, face_location, corners, edges) - def _mesh_1d_param(self, *args): - # In the future, when implementing specific coordinate labels (e.g. (r,z) for 2D - # cylindrical), we may create subclasses for CellSize, CellLocation, - # FaceLocation, and pass the suitable subclasses as the first three positional - # arguments to this _mesh_2d_param method, before *args. This will - # allow this method to apply the suitable subclass handling the - # coordinate labels. Of course, we should start with renaming the existing - # 'internal' x,y,z to _x,_y,_z (same for xvalues, yvalues, zvalues) - + def _mesh_1d_param(self, *args, coordlabels={'x':'_x'}): if len(args) == 1: # Use face locations facelocationX = args[0] @@ -237,26 +316,34 @@ def _mesh_1d_param(self, *args): cell_size_x = np.hstack([facelocationX[1]-facelocationX[0], facelocationX[1:]-facelocationX[0:-1], facelocationX[-1]-facelocationX[-2]]) - cell_size = CellSize(cell_size_x, np.array([0.0]), np.array([0.0])) + cell_size = CellSize(cell_size_x, np.array([0.0]), np.array([0.0]), + coordlabels) cell_location = CellLocation( - 0.5*(facelocationX[1:]+facelocationX[0:-1]), np.array([0.0]), np.array([0.0])) + 0.5*(facelocationX[1:]+facelocationX[0:-1]), + np.array([0.0]), + np.array([0.0]), + coordlabels) face_location = FaceLocation( - facelocationX, np.array([0.0]), np.array([0.0])) + facelocationX, np.array([0.0]), np.array([0.0]), + coordlabels) elif len(args) == 2: # Use number of cells and domain length Nx = args[0] Width = args[1] dx = Width/Nx cell_size = CellSize( - dx*np.ones(Nx+2), np.array([0.0]), np.array([0.0])) + dx*np.ones(Nx+2), np.array([0.0]), np.array([0.0]), + coordlabels) cell_location = CellLocation( int_range(1, Nx)*dx-dx/2, np.array([0.0]), - np.array([0.0])) + np.array([0.0]), + coordlabels) face_location = FaceLocation( int_range(0, Nx)*dx, np.array([0.0]), - np.array([0.0])) + np.array([0.0]), + coordlabels) dims = np.array([Nx], dtype=int) cellsize = cell_size cellcenters = cell_location @@ -277,53 +364,57 @@ def cell_numbers(self): class CylindricalGrid1D(Grid1D): """Mesh based on a 1D cylindrical grid (r) + ======================================= + + This class can be instantiated in different ways: from a list of cell face + locations or from the number of cells and domain length. - The volume elements are cylindrical shells around a central axis + Instantiation Options: + ---------------------- + - CylindricalGrid1D(Nr, Lr) + - CylindricalGrid1D(face_locationsR) + + + Parameters + ---------- + CylindricalGrid1D(Nr, Lr) + Nr : int + Number of cells in the r direction. + Lr : float + Length of the domain in the r direction. + + CylindricalGrid1D(face_locationsR) + face_locationsR : ndarray + Locations of the cell faces in the r direction. + + Examples + -------- + >>> import numpy as np + >>> from pyfvtool import CylindricalGrid1D + >>> mesh = CylindricalGrid1D(10, 10.0) + >>> print(mesh) """ + @overload - def __init__(self, Nx: int, Lx: float): + def __init__(self, Nr: int, Lr: float): ... @overload - def __init__(self, face_locations: np.ndarray): + def __init__(self, face_locationsR: np.ndarray): ... @overload def __init__(self, dims, cellsize, - cellcenters, facecenters, corners, edges): + cellcenters, facecenters, corners, edges): ... def __init__(self, *args): - """Create a CylindricalGrid1D object from a list of cell face locations or from - number of cells and domain length. - - Parameters - ---------- - Nx : int - Number of cells in the x direction. - Lx : float - Length of the domain in the x direction. - face_locations : ndarray - Locations of the cell faces in the x direction. - - Returns - ------- - CylindricalGrid1D - A 1D cylindrical mesh object. - - Examples - -------- - >>> import numpy as np - >>> from pyfvtool import CylindricalGrid1D - >>> mesh = CylindricalGrid1D(10, 10.0) - >>> print(mesh) - """ if (len(args)==6): dims, cell_size, cell_location, face_location, corners, edges\ = args else: dims, cell_size, cell_location, face_location, corners, edges\ - = self._mesh_1d_param(*args) + = self._mesh_1d_param(*args, coordlabels={'r':'_x'}) super().__init__(dims, cell_size, cell_location, face_location, corners, edges) @@ -335,15 +426,43 @@ def __repr__(self): class SphericalGrid1D(Grid1D): """Mesh based on a 1D spherical grid (r) + ===================================== + + This class can be instantiated in different ways: from a list of cell face + locations or from the number of cells and domain length. + + Instantiation Options: + ---------------------- + - SphericalGrid1D(Nr, Lr) + - SphericalGrid1D(face_locationsR) + + + Parameters + ---------- + SphericalGrid1D(Nr, Lr) + Nr : int + Number of cells in the r direction. + Lr : float + Length of the domain in the r direction. + + SphericalGrid1D(face_locationsR) + face_locationsR : ndarray + Locations of the cell faces in the r direction. - The volume elements are concentric spherical shells + Examples + -------- + >>> import numpy as np + >>> from pyfvtool import SphericalGrid1D + >>> mesh = SphericalGrid1D(10, 10.0) + >>> print(mesh) """ + @overload - def __init__(self, Nx: int, Lx: float): + def __init__(self, Nr: int, Lr: float): ... @overload - def __init__(self, face_locations: np.ndarray): + def __init__(self, face_locationsR: np.ndarray): ... @overload @@ -352,47 +471,17 @@ def __init__(self, dims, cellsize, ... def __init__(self, *args): - """ - Create a SphericalGrid1D object from a list of cell face locations or from - number of cells and domain length. - - Parameters - ---------- - Nx : int - Number of cells in the x direction. - Lx : float - Length of the domain in the x direction. - face_locations : ndarray - Locations of the cell faces in the x direction. - - Returns - ------- - SphericalGrid1D - A 1D spherical mesh object. - - Examples - -------- - >>> import numpy as np - >>> from pyfvtool import SphericalGrid1D - >>> mesh = SphericalGrid1D(10, 10.0) - >>> print(mesh) - - Notes - ----- - The mesh is created in spherical coordinates. - """ if (len(args)==6): dims, cell_size, cell_location, face_location, corners, edges\ = args else: dims, cell_size, cell_location, face_location, corners, edges\ - = self._mesh_1d_param(*args) + = self._mesh_1d_param(*args, coordlabels={'r':'_x'}) super().__init__(dims, cell_size, cell_location, face_location, corners, edges) def __repr__(self): - print(f"1D Spherical mesh with Nr={self.dims[0]} cells") - return "" + return f"1D Spherical mesh with Nr={self.dims[0]} cells" @@ -401,8 +490,45 @@ def __repr__(self): class Grid2D(MeshStructure): """Mesh based on a 2D Cartesian grid (x, y) + ======================================== + + This class can be instantiated in different ways: from a list of cell face + locations or from the number of cells and domain length. + + Instantiation Options: + ---------------------- + - Grid2D(Nx, Ny, Lx, Ly) + - Grid2D(face_locationsX, face_locationsY) + + Parameters + ---------- + Grid2D(Nx, Ny, Lx, Ly) + Nx : int + Number of cells in the x direction. + Ny : int + Number of cells in the y direction. + Lx : float + Length of the domain in the x direction. + Ly : float + Length of the domain in the y direction. + + Grid2D(face_locationsX, face_locationsY) + face_locationsX : ndarray + Locations of the cell faces in the x direction. + face_locationsY : ndarray + Locations of the cell faces in the y direction. + + + Examples + -------- + >>> import numpy as np + >>> from pyfvtool import Grid2D + >>> mesh = Grid2D(10, 10, 10.0, 10.0) + >>> print(mesh) """ + + @overload def __init__(self, Nx: int, Ny: int, Lx: float, Ly: float): ... @@ -419,36 +545,7 @@ def __init__(self, dims, cellsize, ... def __init__(self, *args): - """Create a Grid2D object from a list of cell face locations or from - number of cells and domain length. - - Parameters - ---------- - Nx : int - Number of cells in the x direction. - Ny : int - Number of cells in the y direction. - Lx : float - Length of the domain in the x direction. - Ly : float - Length of the domain in the y direction. - face_locationsX : ndarray - Locations of the cell faces in the x direction. - face_locationsY : ndarray - Locations of the cell faces in the y direction. - - Returns - ------- - Grid2D - A 2D mesh object. - - Examples - -------- - >>> import numpy as np - >>> from pyfvtool import Grid2D - >>> mesh = Grid2D(10, 10, 10.0, 10.0) - >>> print(mesh) - """ + if (len(args)==6): dims, cell_size, cell_location, face_location, corners, edges\ = args @@ -458,15 +555,8 @@ def __init__(self, *args): super().__init__(dims, cell_size, cell_location, face_location, corners, edges) - def _mesh_2d_param(self, *args): - # In the future, when implementing specific coordinate labels (e.g. (r,z) for 2D - # cylindrical), we may create subclasses for CellSize, CellLocation, - # FaceLocation, and pass the suitable subclasses as the first three positional - # arguments to this _mesh_2d_param method, before *args. This will - # allow this method to apply the suitable subclass handling the - # coordinate labels. Of course, we should start with renaming the existing - # 'internal' x,y,z to _x,_y,_z (same for xvalues, yvalues, zvalues) - + def _mesh_2d_param(self, *args, coordlabels={'x':'_x', + 'y':'_y'}): if len(args) == 2: # Use face locations facelocationX = args[0] @@ -475,15 +565,18 @@ def _mesh_2d_param(self, *args): Ny = facelocationY.size-1 cell_size = CellSize(self._facelocation_to_cellsize(facelocationX), self._facelocation_to_cellsize(facelocationY), - np.array([0.0])) + np.array([0.0]), + coordlabels) cell_location = CellLocation( 0.5*(facelocationX[1:]+facelocationX[0:-1]), 0.5*(facelocationY[1:]+facelocationY[0:-1]), - np.array([0.0])) + np.array([0.0]), + coordlabels) face_location = FaceLocation( facelocationX, facelocationY, - np.array([0.0])) + np.array([0.0]), + coordlabels) elif len(args) == 4: # Use number of cells and domain length Nx = args[0] @@ -495,15 +588,18 @@ def _mesh_2d_param(self, *args): cell_size = CellSize( dx*np.ones(Nx+2), dy*np.ones(Ny+2), - np.array([0.0])) + np.array([0.0]), + coordlabels) cell_location = CellLocation( int_range(1, Nx)*dx-dx/2, int_range(1, Ny)*dy-dy/2, - np.array([0.0])) + np.array([0.0]), + coordlabels) face_location = FaceLocation( int_range(0, Nx)*dx, int_range(0, Ny)*dy, - np.array([0.0])) + np.array([0.0]), + coordlabels) dims = np.array([Nx, Ny], dtype=int) cellsize = cell_size @@ -515,8 +611,7 @@ def _mesh_2d_param(self, *args): return dims, cellsize, cellcenters, facecenters, corners, edges def __repr__(self): - print(f"2D Cartesian mesh with {self.dims[0]}x{self.dims[1]} cells") - return "" + return f"2D Cartesian mesh with {self.dims[0]}x{self.dims[1]} cells" def cell_numbers(self): Nx, Ny = self.dims @@ -524,18 +619,58 @@ def cell_numbers(self): return G.reshape(Nx+2, Ny+2) + class CylindricalGrid2D(Grid2D): """Mesh based on a 2D cylindrical grid (r, z) - + ========================================== + + This class can be instantiated in different ways: from a list of cell face + locations or from the number of cells and domain length. + + Instantiation Options: + ---------------------- + - CylindricalGrid2D(Nr, Nz, Lr, Lz) + - CylindricalGrid2D(face_locationsR, face_locationsZ) + + Parameters + ---------- + CylindricalGrid2D(Nr, Nz, Lr, Lz) + Nr : int + Number of cells in the r direction. + Nz : int + Number of cells in the z direction. + Lr : float + Length of the domain in the r direction. + Lz : float + Length of the domain in the z direction. + + CylindricalGrid2D(face_locationsR, face_locationsZ) + face_locationsR : ndarray + Locations of the cell faces in the r direction. + face_locationsZ : ndarray + Locations of the cell faces in the z direction. + + Returns + ------- + CylindricalGrid2D + A 2D cylindrical mesh object. + + Examples + -------- + >>> import numpy as np + >>> from pyfvtool import CylindricalGrid2D + >>> mesh = CylindricalGrid2D(10, 10, 10.0, 10.0) + >>> print(mesh) """ + @overload - def __init__(self, Nx: int, Ny: int, - Lx: float, Ly: float): + def __init__(self, Nr: int, Nz: int, + Lr: float, Lz: float): ... @overload - def __init__(self, face_locationsX: np.ndarray, - face_locationsY: np.ndarray): + def __init__(self, face_locationsR: np.ndarray, + face_locationsZ: np.ndarray): ... @overload @@ -544,46 +679,17 @@ def __init__(self, dims, cellsize, ... def __init__(self, *args): - """Create a CylindricalGrid2D object from a list of cell face locations or from - number of cells and domain length. - - Parameters - ---------- - Nx : int - Number of cells in the x direction. - Ny : int - Number of cells in the y direction. - Lx : float - Length of the domain in the x direction. - Ly : float - Length of the domain in the y direction. - face_locationsX : ndarray - Locations of the cell faces in the x direction. - face_locationsY : ndarray - Locations of the cell faces in the y direction. - - Returns - ------- - CylindricalGrid2D - A 2D cylindrical mesh object. - - Examples - -------- - >>> import numpy as np - >>> from pyfvtool import CylindricalGrid2D - >>> mesh = CylindricalGrid2D(10, 10, 10.0, 10.0) - >>> print(mesh) - """ + if (len(args)==6): dims, cell_size, cell_location, face_location, corners, edges\ = args else: dims, cell_size, cell_location, face_location, corners, edges\ - = self._mesh_2d_param(*args) + = self._mesh_2d_param(*args, coordlabels={'r':'_x', + 'z':'_y'}) super().__init__(dims, cell_size, cell_location, face_location, corners, edges) - def __repr__(self): print( f"2D Cylindrical mesh with Nr={self.dims[0]}xNz={self.dims[1]} cells") @@ -593,60 +699,58 @@ def __repr__(self): class PolarGrid2D(Grid2D): """Mesh based on a 2D polar grid (r, theta) - + ======================================== + + This class can be instantiated in different ways: from a list of cell face + locations or from the number of cells and domain length. + + Instantiation Options: + ---------------------- + - PolarGrid2D(Nr, Ntheta, Lr, Ltheta) + - PolarGrid2D(face_locationsR, face_locationsTheta) + + Parameters: + ----------- + PolarGrid2D(Nr, Ntheta, Lr, Ltheta): + Nr : int + Number of cells in the r direction. + Ntheta : int + Number of cells in the theta direction. + Lr : float + Length of the domain in the r direction. + Ltheta : float + Length of the domain in the theta direction. + + PolarGrid2D(face_locationsR, face_locationsTheta): + face_locationsR : ndarray + Locations of the cell faces in the r direction. + face_locationsTheta : ndarray + Locations of the cell faces in the theta direction. + + Examples: + --------- + >>> import numpy as np + >>> from pyfvtool import PolarGrid2D + >>> mesh = PolarGrid2D(10, 10, 10.0, 10.0) + >>> print(mesh) """ @overload - def __init__(self, Nx: int, Ny: int, Lx: float, Ly: float): + def __init__(self, Nr: int, Ntheta: int, Lr: float, Ltheta: float): ... @overload - def __init__(self, face_locationsX: np.ndarray, - face_locationsY: np.ndarray): + def __init__(self, face_locationsR: np.ndarray, + face_locationsTheta: np.ndarray): ... @overload def __init__(self, dims, cellsize, - cellcenters, facecenters, corners, edges): + cellcenters, facecenters, corners, edges): ... def __init__(self, *args): - """ - Create a PolarGrid2D object from a list of cell face locations or from - number of cells and domain length. - - Parameters - ---------- - Nx : int - Number of cells in the x direction. - Ny : int - Number of cells in the y direction. - Lx : float - Length of the domain in the x direction. - Ly : float - Length of the domain in the y direction. - face_locationsX : ndarray - Locations of the cell faces in the x direction. - face_locationsY : ndarray - Locations of the cell faces in the y direction. - - Returns - ------- - PolarGrid2D - A 2D radial mesh object. - - Examples - -------- - >>> import numpy as np - >>> from pyfvtool import PolarGrid2D - >>> mesh = PolarGrid2D(10, 10, 10.0, 10.0) - >>> print(mesh) - - Notes - ----- - The mesh is created in radial (cylindrical) coordinates. - """ if (len(args)==6): dims, cell_size, cell_location,\ face_location, corners, edges= args @@ -658,14 +762,13 @@ def __init__(self, *args): if (theta_max > 2*np.pi): warn("Recreate the mesh with an upper bound of 2*pi for \\theta or there will be unknown consequences!") dims, cell_size, cell_location, face_location, corners, edges\ - = self._mesh_2d_param(*args) + = self._mesh_2d_param(*args, coordlabels={'r' :'_x', + 'theta':'_y'}) super().__init__(dims, cell_size, cell_location, face_location, corners, edges) def __repr__(self): - print( - f"2D Polar mesh with Nr={self.dims[0]}xN_theta={self.dims[1]} cells") - return "" + return f"2D Polar mesh with N_r={self.dims[0]}xN_theta={self.dims[1]} cells" @@ -674,30 +777,24 @@ def __repr__(self): class Grid3D(MeshStructure): - """Mesh based on a 3D Cartesian grid (x, y, z)""" - @overload - def __init__(self, Nx: int, Ny: int, Nz: int, - Lx: float, Ly: float, Lz: float): - ... - - @overload - def _init__(self, face_locationsX: np.ndarray, - face_locationsY: np.ndarray, - face_locationsZ: np.ndarray): - ... - - @overload - def __init__(self, dims, cellsize, - cellcenters, facecenters, corners, edges): - ... - - def __init__(self, *args): - """ - Create a Grid3D object from a list of cell face locations or from - number of cells and domain length. + """Mesh based on a 3D Cartesian grid (x, y, z) + =========================================== - Parameters - ---------- + This class can be instantiated in different ways: from a list of cell face + locations or from the number of cells and domain length. There are multiple + overloaded '__init__' methods available to provide flexibility in + instantiation. + + + Instantiation Options: + ---------------------- + - Grid3D(Nx, Ny, Nz, Lx, Ly, Lz) + - Grid3D(face_locationsX, face_locationsY, face_locationsZ) + + + Parameters + ---------- + Grid3D(Nx, Ny, Nz, Lx, Ly, Lz) Nx : int Number of cells in the x direction. Ny : int @@ -710,6 +807,8 @@ def __init__(self, *args): Length of the domain in the y direction. Lz : float Length of the domain in the z direction. + + Grid3D(face_locationsX, face_locationsY, face_locationsZ) face_locationsX : ndarray Locations of the cell faces in the x direction. face_locationsY : ndarray @@ -717,18 +816,32 @@ def __init__(self, *args): face_locationsZ : ndarray Locations of the cell faces in the z direction. - Returns - ------- - Grid3D - A 3D mesh object. - - Examples - -------- - >>> import numpy as np - >>> from pyfvtool import Grid3D - >>> mesh = Grid3D(10, 10, 10, 10.0, 10.0, 10.0) - >>> print(mesh) - """ + + Examples + -------- + >>> import numpy as np + >>> from pyfvtool import Grid3D + >>> mesh = Grid3D(10, 10, 10, 10.0, 10.0, 10.0) + >>> print(mesh) + + """ + @overload + def __init__(self, Nx: int, Ny: int, Nz: int, + Lx: float, Ly: float, Lz: float): + ... + + @overload + def _init__(self, face_locationsX: np.ndarray, + face_locationsY: np.ndarray, + face_locationsZ: np.ndarray): + ... + + @overload + def __init__(self, dims, cellsize, + cellcenters, facecenters, corners, edges): + ... + + def __init__(self, *args): direct_init = False # Flag to indicate if this is a 'direct' __init__ # not requiring any parsing of arguments. # These 'direct' instantiantions are used @@ -748,15 +861,9 @@ def __init__(self, *args): face_location, corners, edges) - def _mesh_3d_param(self, *args): - # In the future, when implementing specific coordinate labels (e.g. (r,z) for 2D - # cylindrical), we may create subclasses for CellSize, CellLocation, - # FaceLocation, and pass the suitable subclasses as the first three positional - # arguments to this _mesh_3d_param method, before *args. This will - # allow this method to apply the suitable subclass handling the - # coordinate labels. Of course, we should start with renaming the existing - # 'internal' x,y,z to _x,_y,_z (same for xvalues, yvalues, zvalues) - + def _mesh_3d_param(self, *args, coordlabels={'x':'_x', + 'y':'_y', + 'z':'_z'}): if len(args) == 3: # Use face locations facelocationX = args[0] @@ -767,15 +874,18 @@ def _mesh_3d_param(self, *args): Nz = facelocationZ.size-1 cell_size = CellSize(self._facelocation_to_cellsize(facelocationX), self._facelocation_to_cellsize(facelocationY), - self._facelocation_to_cellsize(facelocationZ)) + self._facelocation_to_cellsize(facelocationZ), + coordlabels) cell_location = CellLocation( 0.5*(facelocationX[1:]+facelocationX[0:-1]), 0.5*(facelocationY[1:]+facelocationY[0:-1]), - 0.5*(facelocationZ[1:]+facelocationZ[0:-1])) + 0.5*(facelocationZ[1:]+facelocationZ[0:-1]), + coordlabels) face_location = FaceLocation( facelocationX, facelocationY, - facelocationZ) + facelocationZ, + coordlabels) elif len(args) == 6: # Use number of cells and domain length Nx = args[0] @@ -790,15 +900,18 @@ def _mesh_3d_param(self, *args): cell_size = CellSize( dx*np.ones(Nx+2), dy*np.ones(Ny+2), - dz*np.ones(Nz+2)) + dz*np.ones(Nz+2), + coordlabels) cell_location = CellLocation( int_range(1, Nx)*dx-dx/2, int_range(1, Ny)*dy-dy/2, - int_range(1, Nz)*dz-dz/2) + int_range(1, Nz)*dz-dz/2, + coordlabels) face_location = FaceLocation( int_range(0, Nx)*dx, int_range(0, Ny)*dy, - int_range(0, Nz)*dz) + int_range(0, Nz)*dz, + coordlabels) G = int_range(1, (Nx+2)*(Ny+2)*(Nz+2))-1 G = G.reshape(Nx+2, Ny+2, Nz+2) dims = np.array([Nx, Ny, Nz], dtype=int) @@ -820,21 +933,67 @@ def cell_numbers(self): return G.reshape(Nx+2, Ny+2, Nz+2) def __repr__(self): - print( - f"3D Cartesian mesh with Nx={self.dims[0]}xNy={self.dims[1]}xNz={self.dims[1]} cells") - return "" + return f"3D Cartesian mesh with "\ + f"Nx={self.dims[0]}xNy={self.dims[1]}xNz={self.dims[1]} cells" + + class CylindricalGrid3D(Grid3D): - """Mesh based on a 3D cylindrical grid (r, theta, z)""" + """Mesh based on a 3D cylindrical grid (r, theta, z) + ================================================= + + This class can be instantiated in different ways: from a list of cell face + locations or from the number of cells and domain length. + + + Instantiation Options: + ---------------------- + - CylindricalGrid3D(Nr, Ntheta, Nz, Lr, Ltheta, Lz) + - CylindricalGrid3D(face_locationsR, face_locationsTheta, face_locationsZ) + + Parameters + ---------- + CylindricalGrid3D(Nr, Ntheta, Nz, Lr, Ltheta, Lz) + Nr : int + Number of cells in the r direction. + Ntheta : int + Number of cells in the theta direction. + Nz : int + Number of cells in the z direction. + Lr : float + Length of the domain in the r direction. + Ltheta : float + Length of the domain in the theta direction. + Lz : float + Length of the domain in the z direction. + + + CylindricalGrid3D(face_locationsR, face_locationsTheta, face_locationsZ) + face_locationsR : ndarray + Locations of the cell faces in the r direction. + face_locationsTheta : ndarray + Locations of the cell faces in the theta direction. + face_locationsZ : ndarray + Locations of the cell faces in the z direction. + + + Examples + -------- + >>> import numpy as np + >>> from pyfvtool import CylindricalGrid3D + >>> mesh = CylindricalGrid3D(10, 10, 10, 10.0, 10.0, 10.0) + >>> print(mesh) + """ + @overload - def __init__(self, Nx: int, Ny: int, Nz: int, - Lx: float, Ly: float, Lz: float): + def __init__(self, Nr: int, Ntheta: int, Nz: int, + Lr: float, Ltheta: float, Lz: float): ... @overload - def __init__(self, face_locationsX: np.ndarray, - face_locationsY: np.ndarray, + def __init__(self, face_locationsR: np.ndarray, + face_locationsTheta: np.ndarray, face_locationsZ: np.ndarray): ... @@ -846,47 +1005,7 @@ def __init__(self, dims, cellsize, def __init__(self, *args): """ - Create a CylindricalGrid3D object from a list of cell face locations or from - number of cells and domain length. - TO DO: docstring (and coordinate labels) -> r, theta, z - - Parameters - ---------- - Nx : int - Number of cells in the x direction. - Ny : int - Number of cells in the y direction. - Nz : int - Number of cells in the z direction. - Lx : float - Length of the domain in the x direction. - Ly : float - Length of the domain in the y direction. - Lz : float - Length of the domain in the z direction. - face_locationsX : ndarray - Locations of the cell faces in the x direction. - face_locationsY : ndarray - Locations of the cell faces in the y direction. - face_locationsZ : ndarray - Locations of the cell faces in the z direction. - - Returns - ------- - CylindricalGrid3D - A 3D cylindrical mesh object. - - Examples - -------- - >>> import numpy as np - >>> from pyfvtool import CylindricalGrid3D - >>> mesh = CylindricalGrid3D(10, 10, 10, 10.0, 10.0, 10.0) - >>> print(mesh) - - Notes - ----- - The mesh is created in cylindrical coordinates. """ direct_init = False # Flag to indicate if this is a 'direct' __init__ # not requiring any parsing of arguments. @@ -910,29 +1029,76 @@ def __init__(self, *args): warn("Recreate the mesh with an upper bound of 2*pi for theta or there will be unknown consequences!") dims, cell_size, cell_location, face_location, corners, edges\ - = self._mesh_3d_param(*args) + = self._mesh_3d_param(*args, coordlabels={'r' :'_x', + 'theta':'_y', + 'z' :'_z'}) super().__init__(dims, cell_size, cell_location, face_location, corners, edges) def __repr__(self): - print( - f"3D Cylindrical mesh with Nr={self.dims[0]}xN_theta={self.dims[1]}xNz={self.dims[1]} cells") - return "" + return f"3D Cylindrical mesh with Nr={self.dims[0]}x"\ + f"N_theta={self.dims[1]}xNz={self.dims[1]} cells" + class SphericalGrid3D(Grid3D): - """Mesh based on a 3D spherical grid (r, theta, phi)""" + """Mesh based on a 3D spherical grid (r, theta, phi) + ================================================= + + Create a SphericalGrid3D object from a list of cell face locations or from + the number of cells and domain length. + + + Instantiation Options: + ---------------------- + - SphericalGrid3D(Nr, Ntheta, Nphi, Lr, Ltheta, Lphi) + - SphericalGrid3D(face_locationsR, face_locationsTheta, face_locationsPhi) + + Parameters + ---------- + SphericalGrid3D(Nr, Ntheta, Nphi, Lr, Ltheta, Lphi) + Nr : int + Number of cells in the r direction. + Ntheta : int + Number of cells in the theta direction. + Nphi : int + Number of cells in the phi direction. + Lr : float + Length of the domain in the r direction. + Ltheta : float + Length of the domain in the theta direction. + Lphi : float + Length of the domain in the phi direction. + + + SphericalGrid3D(face_locationsR, face_locationsTheta, face_locationsPhi) + face_locationsR : ndarray + Locations of the cell faces in the r direction. + face_locationsTheta : ndarray + Locations of the cell faces in the theta direction. + face_locationsPhi : ndarray + Locations of the cell faces in the phi direction. + + + Examples + -------- + >>> import numpy as np + >>> from pyfvtool import SphericalGrid3D + >>> mesh = SphericalGrid3D(10, 10, 10, 10.0, 10.0, 10.0) + >>> print(mesh) + """ + @overload - def __init__(self, Nx: int, Ny: int, Nz: int, - Lx: float, Ly: float, Lz: float): + def __init__(self, Nr: int, Ntheta: int, Nphi: int, + Lr: float, Ltheta: float, Lphi: float): ... @overload - def __init__(self, face_locationsX: np.ndarray, - face_locationsY: np.ndarray, - face_locationsZ: np.ndarray): + def __init__(self, face_locationsR: np.ndarray, + face_locationsTheta: np.ndarray, + face_locationsPhi: np.ndarray): ... @overload @@ -942,47 +1108,6 @@ def __init__(self, dims, cellsize, def __init__(self, *args): - """ - Create a SphericalGrid3D object from a list of cell face locations or from - number of cells and domain length. - - Parameters - ---------- - Nx : int - Number of cells in the x direction. - Ny : int - Number of cells in the y direction. - Nz : int - Number of cells in the z direction. - Lx : float - Length of the domain in the x direction. - Ly : float - Length of the domain in the y direction. - Lz : float - Length of the domain in the z direction. - face_locationsX : ndarray - Locations of the cell faces in the x direction. - face_locationsY : ndarray - Locations of the cell faces in the y direction. - face_locationsZ : ndarray - Locations of the cell faces in the z direction. - - Returns - ------- - SphericalGrid3D - A 3D spherical mesh object. - - Examples - -------- - >>> import numpy as np - >>> from pyfvtool import SphericalGrid3D - >>> mesh = SphericalGrid3D(10, 10, 10, 10.0, 10.0, 10.0) - >>> print(mesh) - - Notes - ----- - The mesh is created in spherical coordinates. - """ direct_init = False # Flag to indicate if this is a 'direct' __init__ # not requiring any parsing of arguments. # These 'direct' instantiantions are used @@ -1004,13 +1129,13 @@ def __init__(self, *args): warn("Recreate the mesh with an upper bound of 2*pi for \\phi"\ " or there will be unknown consequences!") dims, cell_size, cell_location, face_location, corners, edges\ - = self._mesh_3d_param(*args) + = self._mesh_3d_param(*args, coordlabels={'r' :'_x', + 'theta':'_y', + 'phi' :'_z'}) super().__init__(dims, cell_size, cell_location, face_location, corners, edges) def __repr__(self): - print( - f"3D Shperical mesh with Nr={self.dims[0]}xN_theta={self.dims[1]}xN_phi={self.dims[1]} cells") - return "" - + return f"3D Shperical mesh with Nr={self.dims[0]}x"\ + "N_theta={self.dims[1]}xN_phi={self.dims[1]} cells" diff --git a/src/pyfvtool/pdesolver.py b/src/pyfvtool/pdesolver.py index 4082977..68697cf 100644 --- a/src/pyfvtool/pdesolver.py +++ b/src/pyfvtool/pdesolver.py @@ -98,15 +98,12 @@ def solvePDE(phi: CellVariable, bcterm: tuple, eqnterms: list, The updated CellVariable. """ - # TODO: Presently the bcterm and the eqnterms are simply the objects - # returned by the respective xxxTerm routines. This is done for - # simplicity. Later, we could consider - # a specific `Term` class, with properties M, RHS and arithmetic - # magic methods defined suitably. Certain terms will have `None` - # for either M or RHS etc. - # A first version may simply be a 'signed tuple' class. (tuple with - # additional sign property such that we can write `-Term` in all - # cases) + # Presently the bcterm and the eqnterms are simply the objects + # returned by the respective xxxTerm routines. This is done for + # simplicity. Later, we might perhaps consider + # a specific `Term` class, with properties M, RHS and arithmetic + # magic methods defined suitably. Certain terms will have `None` + # for either M or RHS etc. if externalsolver is None: solver = spsolve else: diff --git a/src/pyfvtool/utilities.py b/src/pyfvtool/utilities.py index 4a1149a..4e3fcfe 100644 --- a/src/pyfvtool/utilities.py +++ b/src/pyfvtool/utilities.py @@ -10,6 +10,8 @@ def int_range(a:int, b:int) -> np.ndarray: """ return np.linspace(a, b, b-a+1, dtype=int) + + def fluxLimiter(flName: str, eps =2e-16): """ returns a flux limiter function @@ -86,3 +88,42 @@ def FL(r): return FL + +class SignedTuple(tuple): + """Just like a tuple, but supporting negation and unary positive operations + + Can only contain objects that support negation. It does nothing when + the unary positive operator is applied. + + This class is intended for returning matrix equation terms (M and RHS) + such that they can be used with solvePDE without reserve. At present, + no such signed tuples are needed, however. In the cases where tuples of + M and RHS are returned, a negation operation makes no sense (e.g. on + the boundaryConditionsTerm). All other cases, only use pure RHS or pure + M, which already support negation. + """ + def __init__(self, ii): + # When this tuple-subclass object has been created (__new__), + # it is already initialized with the tuple values. + # + # for xx in self: + # print(xx) + # + # So, here, in this __init__ we only need to check if the SignedTuple + # object only contains 'negatable' objects. No need to call super() + # + for i in ii: + if not hasattr(i,'__neg__'): + raise ValueError('Object incompatible with SignedTuple') + + def __pos__(self): + return self + + def __neg__(self): + l = [] + for x in self: + l.append(-x) + return type(self)(l) + + + diff --git a/src/pyfvtool/visualization.py b/src/pyfvtool/visualization.py index 682efeb..f603fbd 100644 --- a/src/pyfvtool/visualization.py +++ b/src/pyfvtool/visualization.py @@ -4,8 +4,6 @@ from .mesh import CylindricalGrid2D from .mesh import PolarGrid2D, CylindricalGrid3D from .cell import CellVariable -from .cell import get_CellVariable_profile1D, get_CellVariable_profile2D -from .cell import get_CellVariable_profile3D import matplotlib.pyplot as plt @@ -39,36 +37,13 @@ def visualizeCells(phi: CellVariable, >>> phi = pf.CellVariable(m, 1.0) >>> pf.visualizeCells(phi) """ - if issubclass(type(phi.domain), Grid1D): - x, phi0 = get_CellVariable_profile1D(phi) - # TODO: - # get_CellVariable_profile1D can become a method of CellVariable - # (shared with 2D and 3D versions) + if isinstance(phi.domain, Grid1D): + x, phi0 = phi.plotprofile() plt.plot(x, phi0) # plt.show() elif (type(phi.domain) is Grid2D) or (type(phi.domain) is CylindricalGrid2D): - x, y, phi0 = get_CellVariable_profile2D(phi) - # TODO: - # get_CellVariable_profile2D can become a method of CellVariable - # (shared with 1D and 3D versions) - ## Kept old code below for reference. Can be removed. - # x = np.hstack([phi.domain.facecenters.x[0], - # phi.domain.cellcenters.x, - # phi.domain.facecenters.x[-1]]) - # y = np.hstack([phi.domain.facecenters.y[0], - # phi.domain.cellcenters.y, - # phi.domain.facecenters.y[-1]]) - # phi0 = np.copy(phi.value) - # phi0[:, 0] = 0.5*(phi0[:, 0]+phi0[:, 1]) - # phi0[0, :] = 0.5*(phi0[0, :]+phi0[1, :]) - # phi0[:, -1] = 0.5*(phi0[:, -1]+phi0[:, -2]) - # phi0[-1, :] = 0.5*(phi0[-1, :]+phi0[-2, :]) - # phi0[0, 0] = phi0[0, 1] - # phi0[0, -1] = phi0[0, -2] - # phi0[-1, 0] = phi0[-1, 1] - # phi0[-1, -1] = phi0[-1, -2] - ## + x, y, phi0 = phi.plotprofile() if vmin is None: vmin = phi0.min() if vmax is None: @@ -79,55 +54,13 @@ def visualizeCells(phi: CellVariable, # plt.show() elif (type(phi.domain) is PolarGrid2D): - x, y, phi0 = get_CellVariable_profile2D(phi) - # TODO: - # get_CellVariable_profile2D can become a method of CellVariable - # (shared with 1D and 3D versions) - ## Kept old code below for reference. Can be removed. - # x = np.hstack([phi.domain.facecenters.x[0], - # phi.domain.cellcenters.x, - # phi.domain.facecenters.x[-1]]) - # y = np.hstack([phi.domain.facecenters.y[0], - # phi.domain.cellcenters.y, - # phi.domain.facecenters.y[-1]]) - # phi0 = np.copy(phi.value) - # phi0[:, 0] = 0.5*(phi0[:, 0]+phi0[:, 1]) - # phi0[0, :] = 0.5*(phi0[0, :]+phi0[1, :]) - # phi0[:, -1] = 0.5*(phi0[:, -1]+phi0[:, -2]) - # phi0[-1, :] = 0.5*(phi0[-1, :]+phi0[-2, :]) - # phi0[0, 0] = phi0[0, 1] - # phi0[0, -1] = phi0[0, -2] - # phi0[-1, 0] = phi0[-1, 1] - # phi0[-1, -1] = phi0[-1, -2] - ## + x, y, phi0 = phi.plotprofile() plt.subplot(111, polar="true") plt.pcolor(y, x, phi0) # plt.show() elif (type(phi.domain) is Grid3D): - x, y, z, phi0 = get_CellVariable_profile3D(phi) - # TODO: - # get_CellVariable_profile3D can become a method of CellVariable - # (shared with 1D and 2D versions) - ## Kept old code below for reference. Can be removed. - # x = np.hstack([phi.domain.facecenters.x[0], - # phi.domain.cellcenters.x, - # phi.domain.facecenters.x[-1]])[:, np.newaxis, np.newaxis] - # y = np.hstack([phi.domain.facecenters.y[0], - # phi.domain.cellcenters.y, - # phi.domain.facecenters.y[-1]])[np.newaxis, :, np.newaxis] - # z = np.hstack([phi.domain.facecenters.z[0], - # phi.domain.cellcenters.z, - # phi.domain.facecenters.z[-1]])[np.newaxis, np.newaxis, :] - # phi0 = np.copy(phi.value) - # phi0[:,0,:]=0.5*(phi0[:,0,:]+phi0[:,1,:]) - # phi0[:,-1,:]=0.5*(phi0[:,-2,:]+phi0[:,-1,:]) - # phi0[:,:,0]=0.5*(phi0[:,:,0]+phi0[:,:,0]) - # phi0[:,:,-1]=0.5*(phi0[:,:,-2]+phi0[:,:,-1]) - # phi0[0,:,:]=0.5*(phi0[1,:,:]+phi0[2,:,:]) - # phi0[-1,:,:]=0.5*(phi0[-2,:,:]+phi0[-1,:,:]) - ## - + x, y, z, phi0 = phi.plotprofile() vmin = np.min(phi0) vmax = np.max(phi0) mynormalize = lambda a:((a - vmin)/(vmax-vmin)) @@ -139,13 +72,7 @@ def visualizeCells(phi: CellVariable, fig = plt.figure() ax = fig.add_subplot(111, projection = "3d") - # r = linspace(1.25, 1.25, 50) - # p = linspace(0, 2π, 50) - # R = repmat(r, 1, 50) - # P = repmat(p', 50, 1) - # Zc = rand(50, 50) # (P.^2-1).^2 - # Z = repmat(linspace(0, 2, 50), 1, 50) - # X, Y = R.*cos.(P), R.*sin.(P) + ax.plot_surface(X[0,:,:], Y[0,:,:], Z[0,:,:], facecolors=plt.cm.viridis(mynormalize(phi0[0,:,:])), alpha=0.8) @@ -167,26 +94,8 @@ def visualizeCells(phi: CellVariable, # plt.show() elif (type(phi.domain) is CylindricalGrid3D): - r, theta, z, phi0 = get_CellVariable_profile3D(phi) - ## Kept old code below for reference. Can be removed. - # r = np.hstack([phi.domain.facecenters.x[0], - # phi.domain.cellcenters.x, - # phi.domain.facecenters.x[-1]])[:, np.newaxis, np.newaxis] - # theta = np.hstack([phi.domain.facecenters.y[0], - # phi.domain.cellcenters.y, - # phi.domain.facecenters.y[-1]])[np.newaxis, :, np.newaxis] - # z = np.hstack([phi.domain.facecenters.z[0], - # phi.domain.cellcenters.z, - # phi.domain.facecenters.z[-1]])[np.newaxis, np.newaxis, :] - # phi0 = np.copy(phi.value) - # phi0[:, 0, :] = 0.5*(phi0[:, 0, :]+phi0[:, 1, :]) - # phi0[:, -1, :] = 0.5*(phi0[:, -2, :]+phi0[:, -1, :]) - # phi0[:, :, 0] = 0.5*(phi0[:, :, 0]+phi0[:, :, 0]) - # phi0[:, :, -1] = 0.5*(phi0[:, :, -2]+phi0[:, :, -1]) - # phi0[0, :, :] = 0.5*(phi0[1, :, :]+phi0[2, :, :]) - # phi0[-1, :, :] = 0.5*(phi0[-2, :, :]+phi0[-1, :, :]) - ## - + r, theta, z, phi0 = phi.plotprofile() + x = r*np.cos(theta) y = r*np.sin(theta) vmin = np.min(phi0) diff --git a/tests/test_coordinate_labels.py b/tests/test_coordinate_labels.py new file mode 100644 index 0000000..3a8d356 --- /dev/null +++ b/tests/test_coordinate_labels.py @@ -0,0 +1,360 @@ +# -*- coding: utf-8 -*- +""" +Created on Fri Mar 29 16:41:50 2024 + +@author: werts-moltech +""" + +# Test coordinate labeling in different mesh geometries + + +import numpy as np + +from pyfvtool import Grid1D, CylindricalGrid1D, SphericalGrid1D +from pyfvtool import Grid2D, CylindricalGrid2D, PolarGrid2D +from pyfvtool import Grid3D, CylindricalGrid3D, SphericalGrid3D + + +errors_expected = 0 +errors_caught = 0 + +msh = Grid1D(10, 1.) +xx = msh.cellcenters.x +assert np.all(xx == msh.cellcenters._x) +errors_expected+=1 +try: + rr = msh.cellcenters.r +except AttributeError: + errors_caught+=1 +xx = msh.cellsize.x +assert np.all(xx == msh.cellsize._x) +errors_expected+=1 +try: + rr = msh.cellsize.r +except AttributeError: + errors_caught+=1 +xx = msh.facecenters.x +assert np.all(xx == msh.facecenters._x) +errors_expected+=1 +try: + rr = msh.facecenters.r +except AttributeError: + errors_caught+=1 + + + +msh = CylindricalGrid1D(10, 1.) +rr = msh.cellcenters.r +assert np.all(rr == msh.cellcenters._x) +errors_expected+=1 +try: + xx = msh.cellcenters.x +except AttributeError: + errors_caught+=1 +rr = msh.cellsize.r +assert np.all(rr == msh.cellsize._x) +errors_expected+=1 +try: + xx = msh.cellsize.x +except AttributeError: + errors_caught+=1 +rr = msh.facecenters.r +assert np.all(rr == msh.facecenters._x) +errors_expected+=1 +try: + xx = msh.facecenters.x +except AttributeError: + errors_caught+=1 + + + +msh = CylindricalGrid2D(10, 10, 1., 1.) +rr = msh.cellcenters.r +zz = msh.cellcenters.z +assert np.all(rr == msh.cellcenters._x) +assert np.all(zz == msh.cellcenters._y) +errors_expected+=1 +try: + xx = msh.cellcenters.x +except AttributeError: + errors_caught+=1 +rr = msh.cellsize.r +zz = msh.cellsize.z +assert np.all(rr == msh.cellsize._x) +assert np.all(zz == msh.cellsize._y) +errors_expected+=1 +try: + xx = msh.cellsize.x +except AttributeError: + errors_caught+=1 +rr = msh.facecenters.r +zz = msh.facecenters.z +assert np.all(rr == msh.facecenters._x) +assert np.all(zz == msh.facecenters._y) +errors_expected+=1 +try: + xx = msh.facecenters.x +except AttributeError: + errors_caught+=1 +errors_expected+=1 +try: + yy = msh.cellcenters.y +except AttributeError: + errors_caught+=1 +errors_expected+=1 +try: + yy = msh.cellsize.y +except AttributeError: + errors_caught+=1 +errors_expected+=1 +try: + yy = msh.facecenters.y +except AttributeError: + errors_caught+=1 + + + +msh = Grid2D(10, 10, 1., 1.) +xx = msh.cellcenters.x +yy = msh.cellcenters.y +assert np.all(xx == msh.cellcenters._x) +assert np.all(yy == msh.cellcenters._y) +errors_expected+=1 +try: + rr = msh.cellcenters.r +except AttributeError: + errors_caught+=1 +xx = msh.cellsize.x +yy = msh.cellsize.y +assert np.all(xx == msh.cellsize._x) +assert np.all(yy == msh.cellsize._y) +errors_expected+=1 +try: + rr = msh.cellsize.r +except AttributeError: + errors_caught+=1 +xx = msh.facecenters.x +yy = msh.facecenters.y +assert np.all(xx == msh.facecenters._x) +assert np.all(yy == msh.facecenters._y) +errors_expected+=1 +try: + rr = msh.facecenters.r +except AttributeError: + errors_caught+=1 +errors_expected+=1 +try: + zz = msh.cellcenters.z +except AttributeError: + errors_caught+=1 +errors_expected+=1 +try: + zz = msh.cellsize.z +except AttributeError: + errors_caught+=1 +errors_expected+=1 +try: + zz = msh.facecenters.z +except AttributeError: + errors_caught+=1 + + + + +msh = PolarGrid2D(10, 10, 1., 1.) +rr = msh.cellcenters.r +theta = msh.cellcenters.theta +assert np.all(rr == msh.cellcenters._x) +assert np.all(theta == msh.cellcenters._y) +errors_expected+=1 +try: + xx = msh.cellcenters.x +except AttributeError: + errors_caught+=1 +rr = msh.cellsize.r +theta = msh.cellsize.theta +assert np.all(rr == msh.cellsize._x) +assert np.all(theta == msh.cellsize._y) +errors_expected+=1 +try: + xx = msh.cellsize.x +except AttributeError: + errors_caught+=1 +rr = msh.facecenters.r +theta = msh.facecenters.theta +assert np.all(rr == msh.facecenters._x) +assert np.all(theta == msh.facecenters._y) +errors_expected+=1 +try: + xx = msh.facecenters.x +except AttributeError: + errors_caught+=1 +errors_expected+=1 +try: + yy = msh.cellcenters.y +except AttributeError: + errors_caught+=1 +errors_expected+=1 +try: + yy = msh.cellsize.y +except AttributeError: + errors_caught+=1 +errors_expected+=1 +try: + yy = msh.facecenters.y +except AttributeError: + errors_caught+=1 + + + + +msh = Grid3D(10, 10, 10, 1., 1., 1.) +xx = msh.cellcenters.x +yy = msh.cellcenters.y +zz = msh.cellcenters.z +assert np.all(xx == msh.cellcenters._x) +assert np.all(yy == msh.cellcenters._y) +assert np.all(zz == msh.cellcenters._z) +errors_expected+=1 +try: + rr = msh.cellcenters.r +except AttributeError: + errors_caught+=1 +xx = msh.cellsize.x +yy = msh.cellsize.y +zz = msh.cellsize.z +assert np.all(xx == msh.cellsize._x) +assert np.all(yy == msh.cellsize._y) +assert np.all(zz == msh.cellsize._z) +errors_expected+=1 +try: + rr = msh.cellsize.r +except AttributeError: + errors_caught+=1 +xx = msh.facecenters.x +yy = msh.facecenters.y +zz = msh.facecenters.z +assert np.all(xx == msh.facecenters._x) +assert np.all(yy == msh.facecenters._y) +assert np.all(zz == msh.facecenters._z) +errors_expected+=1 +try: + rr = msh.facecenters.r +except AttributeError: + errors_caught+=1 + + + +msh = CylindricalGrid3D(10, 10, 10, 1., 2*np.pi, 1.) +rr = msh.cellcenters.r +theta = msh.cellcenters.theta +zz = msh.cellcenters.z +assert np.all(rr == msh.cellcenters._x) +assert np.all(theta == msh.cellcenters._y) +assert np.all(zz == msh.cellcenters._z) +errors_expected+=1 +try: + xx = msh.cellcenters.x +except AttributeError: + errors_caught+=1 +rr = msh.cellsize.r +zz = msh.cellsize.z +assert np.all(rr == msh.cellsize._x) +assert np.all(zz == msh.cellsize._z) +errors_expected+=1 +try: + xx = msh.cellsize.x +except AttributeError: + errors_caught+=1 +rr = msh.facecenters.r +zz = msh.facecenters.z +assert np.all(rr == msh.facecenters._x) +assert np.all(zz == msh.facecenters._z) +errors_expected+=1 +try: + xx = msh.facecenters.x +except AttributeError: + errors_caught+=1 +errors_expected+=1 +try: + yy = msh.cellcenters.y +except AttributeError: + errors_caught+=1 +errors_expected+=1 +try: + yy = msh.cellsize.y +except AttributeError: + errors_caught+=1 +errors_expected+=1 +try: + yy = msh.facecenters.y +except AttributeError: + errors_caught+=1 + + + + + +msh = SphericalGrid3D(10, 15, 20, 1., 2*np.pi, 2*np.pi) +rr = msh.cellcenters.r +theta = msh.cellcenters.theta +phi = msh.cellcenters.phi +assert np.all(rr == msh.cellcenters._x) +assert np.all(theta == msh.cellcenters._y) +assert np.all(phi == msh.cellcenters._z) +errors_expected+=1 +try: + xx = msh.cellcenters.x +except AttributeError: + errors_caught+=1 +rr = msh.cellsize.r +theta = msh.cellsize.theta +phi = msh.cellsize.phi +assert np.all(rr == msh.cellsize._x) +assert np.all(theta == msh.cellsize._y) +assert np.all(phi == msh.cellsize._z) +errors_expected+=1 +try: + xx = msh.cellsize.x +except AttributeError: + errors_caught+=1 +rr = msh.facecenters.r +theta = msh.facecenters.theta +phi = msh.facecenters.phi +assert np.all(rr == msh.facecenters._x) +assert np.all(theta == msh.facecenters._y) +assert np.all(phi == msh.facecenters._z) +errors_expected+=1 +try: + xx = msh.facecenters.x +except AttributeError: + errors_caught+=1 +errors_expected+=1 +try: + yy = msh.cellcenters.y +except AttributeError: + errors_caught+=1 +errors_expected+=1 +try: + yy = msh.cellsize.y +except AttributeError: + errors_caught+=1 +errors_expected+=1 +try: + yy = msh.facecenters.z +except AttributeError: + errors_caught+=1 + + + + + + + + +print('Coordinate label errors expected: ', errors_expected, + ' caught: ', errors_caught) + + +def test_success(): + assert errors_caught == errors_expected \ No newline at end of file diff --git a/tests/test_facevariable.py b/tests/test_facevariable.py new file mode 100644 index 0000000..80d4edf --- /dev/null +++ b/tests/test_facevariable.py @@ -0,0 +1,131 @@ +# -*- coding: utf-8 -*- +""" +Created on Wed Mar 27 12:35:16 2024 + +@author: werts-moltech +""" + +# Test some aspects of FaceVariable handling, +# in particular xvalue, rvalue etc. for different mesh types + +import numpy as np + +from pyfvtool import Grid1D, CylindricalGrid1D, SphericalGrid1D +from pyfvtool import Grid2D, CylindricalGrid2D, PolarGrid2D +from pyfvtool import Grid3D, CylindricalGrid3D, SphericalGrid3D + +from pyfvtool import FaceVariable + + + + + + +errors_expected = 0 +errors_caught = 0 + + + +g1d = Grid1D(10, 1.0) +fv = FaceVariable(g1d, 1.0) +fv.xvalue[:] = 3.0 +errors_expected += 1 +try: + fv.yvalue[:] = 3.0 +except AttributeError: + errors_caught +=1 + + + +g2d = Grid2D(10, 10, 1.0, 1.0) +fv = FaceVariable(g2d, 1.0) +fv.xvalue[:] = 3.0 +fv.yvalue[:] = 4.0 +errors_expected += 1 +try: + fv.zvalue[:] = 3.0 +except AttributeError: + errors_caught +=1 +errors_expected += 1 +try: + fv.rvalue[:] = 3.0 +except AttributeError: + errors_caught +=1 +errors_expected += 1 +try: + fv.thetavalue[:] = 3.0 +except AttributeError: + errors_caught +=1 + + + +g3d = Grid3D(10, 10, 10, 1.0, 1.0, 1.0) +fv = FaceVariable(g3d, 1.0) +fv.xvalue[:] = 3.0 +fv.yvalue[:] = 4.0 +fv.zvalue[:] = 5.0 +errors_expected += 1 +try: + fv.rvalue[:] = 7.0 +except AttributeError: + errors_caught +=1 + + +c1d = CylindricalGrid1D(10, 1.0) +fv = FaceVariable(c1d, 1.0) +fv.rvalue[:] = 3.0 +print(fv.rvalue) +errors_expected += 1 +try: + fv.xvalue[:] = 3.0 +except AttributeError: + errors_caught +=1 +errors_expected += 1 +try: + fv.thetavalue[:] = 3.0 +except AttributeError: + errors_caught +=1 + + +s1d = SphericalGrid1D(10, 1.0) +fv = FaceVariable(s1d, 1.0) +errors_expected += 1 +try: + fv.xvalue[:] = 3.0 +except NotImplementedError: + errors_caught +=1 + + + +c2d = CylindricalGrid2D(10, 10, 1.0, 1.0) +fv = FaceVariable(c2d, 1.0) +errors_expected += 1 +fv.zvalue[:] = 3.0 +# peek inside (testing only) +assert np.all(fv._yvalue == fv.zvalue) +try: + fv.yvalue[:] = 3.0 +except AttributeError: + errors_caught +=1 + + + +s3d = SphericalGrid3D(10, 10, 10, 1.0, 2*np.pi, 2*np.pi) +fv = FaceVariable(s3d, 1.0) +fv.rvalue[:] = 3.0 +fv.thetavalue[:] = 4.0 +print(fv.thetavalue) +fv.phivalue[:] = 5.0 +assert np.all(fv._xvalue == fv.rvalue) +assert np.all(fv._yvalue == fv.thetavalue) +assert np.all(fv._zvalue == fv.phivalue) + + + +print('FaceVariable errors expected: ', errors_expected, + ' caught: ', errors_caught) + + +def test_success(): + assert errors_caught == errors_expected + diff --git a/tests/test_pdesolver_mason_weaver.py b/tests/test_pdesolver_mason_weaver.py index e33005d..b578561 100644 --- a/tests/test_pdesolver_mason_weaver.py +++ b/tests/test_pdesolver_mason_weaver.py @@ -7,7 +7,7 @@ Part. Part. Syst. Charact. 2017, 34, 1700095. doi:10.1002/ppsc.201700095. -using solvePDE v0.2.1 +using solvePDE v0.3.0 """ import numpy as np import matplotlib.pyplot as plt @@ -31,7 +31,7 @@ BC_c = pf.BoundaryConditions(msh) c = pf.CellVariable(msh, 1.0, BC_c) -total_c = [pf.domainInt(c)] +total_c = [c.domainIntegral()] # advection field u = pf.FaceVariable(msh, (sg,)) @@ -70,7 +70,7 @@ bcterm, eqn) it+=1 - total_c.append(pf.domainInt(c)) + total_c.append(c.domainIntegral()) if (it % Nskip == 0): pf.visualizeCells(c)