diff --git a/ci/docs.yml b/ci/docs.yml index 63b32c7c4..3f1eac96d 100644 --- a/ci/docs.yml +++ b/ci/docs.yml @@ -22,4 +22,6 @@ dependencies: - myst-nb - sphinx-design - nbsphinx + - gmpy2 + - mpmath - matplotlib diff --git a/ci/environment.yml b/ci/environment.yml index 1d53cd2bc..c3df18390 100644 --- a/ci/environment.yml +++ b/ci/environment.yml @@ -12,3 +12,5 @@ dependencies: - pytest-cov - scipy - xarray + - gmpy2 + - mpmath diff --git a/ci/upstream-dev-environment.yml b/ci/upstream-dev-environment.yml index e03435753..81db543d1 100644 --- a/ci/upstream-dev-environment.yml +++ b/ci/upstream-dev-environment.yml @@ -10,3 +10,5 @@ dependencies: - pytest - pytest-cov - scipy + - gmpy2 + - mpmath diff --git a/docs/examples.rst b/docs/examples.rst index 7121a2be4..a9bd7fc85 100644 --- a/docs/examples.rst +++ b/docs/examples.rst @@ -20,3 +20,4 @@ your own examples or suggesting one for us to create, please see the :doc:`contr examples/002-grid-topology.ipynb examples/003-area-calc.ipynb examples/004-working-with-mpas-grids.ipynb + examples/005-multiprecision-usage.ipynb diff --git a/docs/examples/000-template.ipynb b/docs/examples/000-template.ipynb index 2c2d42a08..d0fb939cb 100644 --- a/docs/examples/000-template.ipynb +++ b/docs/examples/000-template.ipynb @@ -2,7 +2,11 @@ "cells": [ { "cell_type": "markdown", - "metadata": {}, + "metadata": { + "pycharm": { + "name": "#%% md\n" + } + }, "source": [ "# Blank Template Notebook" ] @@ -10,7 +14,11 @@ { "cell_type": "code", "execution_count": 1, - "metadata": {}, + "metadata": { + "pycharm": { + "name": "#%%\n" + } + }, "outputs": [], "source": [ "import numpy as np\n", @@ -20,7 +28,11 @@ { "cell_type": "code", "execution_count": null, - "metadata": {}, + "metadata": { + "pycharm": { + "name": "#%%\n" + } + }, "outputs": [], "source": [] } diff --git a/docs/examples/001-read-grid-data.ipynb b/docs/examples/001-read-grid-data.ipynb new file mode 100644 index 000000000..343698ab3 --- /dev/null +++ b/docs/examples/001-read-grid-data.ipynb @@ -0,0 +1,1671 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "metadata": { + "pycharm": { + "name": "#%% md\n" + } + }, + "source": [ + "# Reading In Grid Data\n", + "\n", + "UXarray offers support for loading and representing unstructured grids\n", + "by utilizing existing Xarray functionality paired with new routines that\n", + "are specifically written for operating on unstructured grids.\n" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "collapsed": false, + "pycharm": { + "name": "#%% md\n" + } + }, + "source": [ + "## Overview\n", + "\n", + "When working with unstructured grids, the grid definition and data variables\n", + "are often stored as separate files (most of the time as netCDF). This means\n", + "that there are multiple separate files that need to be read in at once.\n", + "\n", + "For example, the NOAA Geoflow project consists of 4 files (1 grid file and 3\n", + "data files). This project follows the UGRID conventions. Special thanks to\n", + "John Clyne, Shilpi Gupta, and the VAPOR team for providing these data!\n", + "\n", + "```\n", + "geoflow-small\n", + "│ grid.nc\n", + "│ v1.nc\n", + "│ v2.nc\n", + "│ v3.nc\n", + "```" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "collapsed": false, + "pycharm": { + "name": "#%% md\n" + } + }, + "source": [ + "## Opening The Grid and Data Files with Xarray\n", + "\n", + "First, read in the grid definition and data variable netCCDF files by using\n", + "`xarray.open_dataset`. In addition, `xr.open_mf_dataset` can be used for\n", + "combining multiple data variables into a single `xarray.Dataset`." + ] + }, + { + "cell_type": "code", + "execution_count": 4, + "metadata": { + "pycharm": { + "name": "#%%\n" + } + }, + "outputs": [], + "source": [ + "import xarray as xr" + ] + }, + { + "cell_type": "code", + "execution_count": 5, + "metadata": { + "pycharm": { + "name": "#%%\n" + } + }, + "outputs": [], + "source": [ + "# Base data path\n", + "base_path = \"../../test/meshfiles/ugrid/geoflow-small/\"\n", + "\n", + "# Path to Grid file\n", + "grid_path = base_path + \"grid.nc\"\n", + "\n", + "# Paths to Data Variable files\n", + "var_names = ['v1.nc', 'v2.nc', 'v3.nc']\n", + "data_paths = [base_path + name for name in var_names]" + ] + }, + { + "cell_type": "code", + "execution_count": 6, + "metadata": { + "pycharm": { + "name": "#%%\n" + } + }, + "outputs": [ + { + "name": "stderr", + "output_type": "stream", + "text": [ + "C:\\Users\\Philip\\projects\\uxarray-redesign\\venv\\lib\\site-packages\\xarray\\backends\\plugins.py:71: RuntimeWarning: Engine 'cfgrib' loading failed:\n", + "Cannot find the ecCodes library\n", + " warnings.warn(f\"Engine {name!r} loading failed:\\n{ex}\", RuntimeWarning)\n" + ] + } + ], + "source": [ + "# Open the Grid file\n", + "grid_ds = xr.open_dataset(grid_path)\n", + "\n", + "# Open a single file or multiple files\n", + "v1_ds = xr.open_dataset(data_paths[0])\n", + "multi_data_ds = xr.open_mfdataset(data_paths)" + ] + }, + { + "cell_type": "code", + "execution_count": 7, + "metadata": { + "collapsed": false, + "pycharm": { + "name": "#%%\n" + } + }, + "outputs": [ + { + "data": { + "text/html": [ + "
\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "
<xarray.Dataset>\n",
+       "Dimensions:          (nMeshFaces: 3840, nFaceNodes: 4, nMeshNodes: 6000,\n",
+       "                      meshLayers: 20)\n",
+       "Coordinates:\n",
+       "    mesh_node_x      (nMeshNodes) float64 ...\n",
+       "    mesh_node_y      (nMeshNodes) float64 ...\n",
+       "Dimensions without coordinates: nMeshFaces, nFaceNodes, nMeshNodes, meshLayers\n",
+       "Data variables:\n",
+       "    mesh             int32 ...\n",
+       "    mesh_face_nodes  (nMeshFaces, nFaceNodes) float64 ...\n",
+       "    mesh_depth       (meshLayers, nMeshNodes) float64 ...
" + ], + "text/plain": [ + "\n", + "Dimensions: (nMeshFaces: 3840, nFaceNodes: 4, nMeshNodes: 6000,\n", + " meshLayers: 20)\n", + "Coordinates:\n", + " mesh_node_x (nMeshNodes) float64 ...\n", + " mesh_node_y (nMeshNodes) float64 ...\n", + "Dimensions without coordinates: nMeshFaces, nFaceNodes, nMeshNodes, meshLayers\n", + "Data variables:\n", + " mesh int32 ...\n", + " mesh_face_nodes (nMeshFaces, nFaceNodes) float64 ...\n", + " mesh_depth (meshLayers, nMeshNodes) float64 ..." + ] + }, + "execution_count": 7, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "grid_ds" + ] + }, + { + "cell_type": "code", + "execution_count": 8, + "metadata": { + "collapsed": false, + "pycharm": { + "name": "#%%\n" + } + }, + "outputs": [ + { + "data": { + "text/html": [ + "
\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "
<xarray.Dataset>\n",
+       "Dimensions:  (time: 1, meshLayers: 20, nMeshNodes: 6000)\n",
+       "Coordinates:\n",
+       "  * time     (time) float64 13.0\n",
+       "Dimensions without coordinates: meshLayers, nMeshNodes\n",
+       "Data variables:\n",
+       "    v1       (time, meshLayers, nMeshNodes) float64 dask.array<chunksize=(1, 20, 6000), meta=np.ndarray>\n",
+       "    v2       (time, meshLayers, nMeshNodes) float64 dask.array<chunksize=(1, 20, 6000), meta=np.ndarray>\n",
+       "    v3       (time, meshLayers, nMeshNodes) float64 dask.array<chunksize=(1, 20, 6000), meta=np.ndarray>
" + ], + "text/plain": [ + "\n", + "Dimensions: (time: 1, meshLayers: 20, nMeshNodes: 6000)\n", + "Coordinates:\n", + " * time (time) float64 13.0\n", + "Dimensions without coordinates: meshLayers, nMeshNodes\n", + "Data variables:\n", + " v1 (time, meshLayers, nMeshNodes) float64 dask.array\n", + " v2 (time, meshLayers, nMeshNodes) float64 dask.array\n", + " v3 (time, meshLayers, nMeshNodes) float64 dask.array" + ] + }, + "execution_count": 8, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "multi_data_ds" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "collapsed": false, + "pycharm": { + "name": "#%% md\n" + } + }, + "source": [ + "## Representing The Unstructured Grid with UXarray" + ] + }, + { + "cell_type": "code", + "execution_count": 9, + "metadata": { + "collapsed": false, + "pycharm": { + "name": "#%%\n" + } + }, + "outputs": [], + "source": [ + "import uxarray as ux" + ] + }, + { + "cell_type": "code", + "execution_count": 10, + "metadata": { + "collapsed": false, + "pycharm": { + "name": "#%%\n" + } + }, + "outputs": [], + "source": [ + "# Construct a UXarray Grid object from `grid_ds`\n", + "grid = ux.Grid(grid_ds)" + ] + }, + { + "cell_type": "code", + "execution_count": 11, + "metadata": { + "collapsed": false, + "pycharm": { + "name": "#%%\n" + } + }, + "outputs": [ + { + "data": { + "text/html": [ + "
\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "
<xarray.Dataset>\n",
+       "Dimensions:          (nMeshFaces: 3840, nFaceNodes: 4, nMeshNodes: 6000,\n",
+       "                      meshLayers: 20)\n",
+       "Coordinates:\n",
+       "    mesh_node_x      (nMeshNodes) float64 ...\n",
+       "    mesh_node_y      (nMeshNodes) float64 ...\n",
+       "Dimensions without coordinates: nMeshFaces, nFaceNodes, nMeshNodes, meshLayers\n",
+       "Data variables:\n",
+       "    mesh             int32 ...\n",
+       "    mesh_face_nodes  (nMeshFaces, nFaceNodes) float64 ...\n",
+       "    mesh_depth       (meshLayers, nMeshNodes) float64 ...
" + ], + "text/plain": [ + "\n", + "Dimensions: (nMeshFaces: 3840, nFaceNodes: 4, nMeshNodes: 6000,\n", + " meshLayers: 20)\n", + "Coordinates:\n", + " mesh_node_x (nMeshNodes) float64 ...\n", + " mesh_node_y (nMeshNodes) float64 ...\n", + "Dimensions without coordinates: nMeshFaces, nFaceNodes, nMeshNodes, meshLayers\n", + "Data variables:\n", + " mesh int32 ...\n", + " mesh_face_nodes (nMeshFaces, nFaceNodes) float64 ...\n", + " mesh_depth (meshLayers, nMeshNodes) float64 ..." + ] + }, + "execution_count": 11, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "grid._ds" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "collapsed": false, + "pycharm": { + "name": "#%% md\n" + } + }, + "source": [ + "As can be seen, the Grid object has its own `grid.ds` of type `xarray.Dataset`\n", + "to define the grid topology. However, the Grid object has further attributes,\n", + "variables, and functions to specify the unstructured grid and be executed on\n", + "it, which can be explored in the next notebooks." + ] + } + ], + "metadata": { + "kernelspec": { + "display_name": "Python 3.10.5 ('uxarray_build')", + "language": "python", + "name": "python3" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 3 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython3", + "version": "3.10.5" + }, + "orig_nbformat": 4, + "vscode": { + "interpreter": { + "hash": "8e8ae2f388051fced6c30f82a529eeca8cf1e059ab06a64326e2a2ad0ec3c36c" + } + } + }, + "nbformat": 4, + "nbformat_minor": 2 +} diff --git a/docs/examples/002-access-grid-info.ipynb b/docs/examples/002-access-grid-info.ipynb new file mode 100644 index 000000000..66277e5fc --- /dev/null +++ b/docs/examples/002-access-grid-info.ipynb @@ -0,0 +1,329 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "metadata": { + "collapsed": false, + "pycharm": { + "name": "#%% md\n" + } + }, + "source": [ + "# Accessing Grid Information\n", + "\n", + "Unstructured grids can be represented in one of many different conventions\n", + "(UGRID, SCRIP, EXODUS, etc). These conventions have different definitions\n", + "and representations of the attributes and variables used to describe\n", + "the unstructured grid topology. Even more, the [UGRID conventions](\n", + "https://ugrid-conventions.github.io/ugrid-conventions/) does not\n", + "enforce standard variable namings for most of the attributes and variables\n", + "(other than just a few required ones).\n", + "\n", + "UXarray unifies all of these conventions at the data loading step by\n", + "representing grids in the UGRID convention regardless of the original grid\n", + "type that is read in from the file. Furthermore, it uses a set of\n", + "standardized names for topology attributes and variables, while still\n", + "providing the user with the original attribute names and variables that\n", + "came from the grid definition file.\n", + "\n", + "## Overview\n", + "\n", + "This notebook will showcase the different methods available for accessing\n", + "the grid topology attributes and variables stored in the `UXarray.Grid`\n", + "object.\n", + "\n", + "For more details on how to load in data, check out our [previous usage\n", + "example](https://uxarray.readthedocs.io/en/latest/examples/read-grid-data.html)\n", + "\n", + "**Methods**\n", + "1. Indexing with Original Variable Names\n", + "2. Indexing with UXarray Variable Dictionary\n", + "3. UXarray's Standardized UGRID Names (Most convenient)" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "collapsed": false, + "pycharm": { + "name": "#%% md\n" + } + }, + "source": [ + "## Data\n", + "\n", + "We will be using two grid files, both of which are in the UGRID convention.\n", + "However, the key difference between them is the names used to describe the\n", + "attributes and variables.\n", + "\n", + "Let us first read in the data:" + ] + }, + { + "cell_type": "code", + "execution_count": 4, + "metadata": { + "collapsed": false, + "pycharm": { + "name": "#%%\n" + } + }, + "outputs": [], + "source": [ + "import uxarray as ux\n", + "import xarray as xr" + ] + }, + { + "cell_type": "code", + "execution_count": 5, + "metadata": { + "collapsed": false, + "pycharm": { + "name": "#%%\n" + } + }, + "outputs": [], + "source": [ + "# Base data path\n", + "base_path = \"../../test/meshfiles/\"\n", + "\n", + "# Path to Grid files\n", + "ugrid_01_path = base_path + \"/ugrid/outCSne30/outCSne30.ug\"\n", + "ugrid_02_path = base_path + \"/ugrid/geoflow-small/grid.nc\"\n", + "\n", + "# Load grid files and create UXarray Grid objects\n", + "ugrid_01_ds = xr.open_dataset(ugrid_01_path)\n", + "ugrid_02_ds = xr.open_dataset(ugrid_02_path)\n", + "\n", + "ugrid_01 = ux.Grid(ugrid_01_ds)\n", + "ugrid_02 = ux.Grid(ugrid_02_ds)" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "collapsed": false, + "pycharm": { + "name": "#%% md\n" + } + }, + "source": [ + "The output of the bottom cell showcases the slight differences\n", + "in variable names:" + ] + }, + { + "cell_type": "code", + "execution_count": 16, + "metadata": { + "collapsed": false, + "pycharm": { + "name": "#%%\n" + } + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "\n", + "Variable Names\n", + "ugrid_01 variable names: ['Mesh2', 'Mesh2_face_nodes', 'Mesh2_node_x', 'Mesh2_node_y', 'nMesh2_face', 'nMaxMesh2_face_nodes', 'nMesh2_node']\n", + "ugrid_02 variable names: ['mesh', 'mesh_face_nodes', 'mesh_depth', 'mesh_node_x', 'mesh_node_y', 'nMeshFaces', 'nFaceNodes', 'nMeshNodes', 'meshLayers']\n" + ] + } + ], + "source": [ + "# Extract variable names\n", + "ugrid_01_names = list(ugrid_01._ds.keys()) + \\\n", + " list(ugrid_01._ds.coords) + \\\n", + " list(ugrid_01._ds.dims)\n", + "ugrid_02_names = list(ugrid_02._ds.keys()) + \\\n", + " list(ugrid_02._ds.coords) + \\\n", + " list(ugrid_02._ds.dims)\n", + "\n", + "print(\"\\nAttribute and variable names for each grid:\")\n", + "print(\"ugrid_01 variable names:\", ugrid_01_names)\n", + "print(\"ugrid_02 variable names:\", ugrid_02_names)" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "collapsed": false, + "pycharm": { + "name": "#%% md\n" + } + }, + "source": [ + "## 1. Indexing with Original Variable Names\n", + "\n", + "The simplest approach is to use the original variable name as an index\n", + "into the grid dataset, `Grid._ds`. Since `ugrid_01` and `ugrid_02` have\n", + "different names for most of their topology attributes and variables, the\n", + "index for accessing them will be different for both grids." + ] + }, + { + "cell_type": "code", + "execution_count": 20, + "metadata": { + "collapsed": false, + "pycharm": { + "name": "#%%\n" + } + }, + "outputs": [], + "source": [ + "x_1 = ugrid_01._ds['Mesh2_node_x']\n", + "y_1 = ugrid_01._ds['Mesh2_node_y']\n", + "face_nodes_1 = ugrid_01._ds['Mesh2_face_nodes']\n", + "n_face_nodes_1 = ugrid_01._ds['nMaxMesh2_face_nodes']\n", + "\n", + "x_2 = ugrid_02._ds['mesh_node_x']\n", + "y_2 = ugrid_02._ds['mesh_node_y']\n", + "face_nodes_2 = ugrid_02._ds['mesh_face_nodes']\n", + "n_face_nodes_2 = ugrid_02._ds['nFaceNodes']" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "collapsed": false, + "pycharm": { + "name": "#%% md\n" + } + }, + "source": [ + "## 2. Indexing with UXarray Variable Dictionary\n", + "\n", + "UXarray provides a dictionary, `Grid._ds_var_names`, to map the original\n", + "topology attribute and variable names that come from the grid file into\n", + "a standardized set of names. In other words, the dictionary uses a\n", + "standardized set of UGRID attribute and variable names as keys, and the\n", + "original variable names that come from the grid file as values.\n", + "\n", + "This allows us to use the same index into either dataset. However, this\n", + "makes the indexing code much longer than the previous method." + ] + }, + { + "cell_type": "code", + "execution_count": 24, + "metadata": { + "collapsed": false, + "pycharm": { + "name": "#%%\n" + } + }, + "outputs": [], + "source": [ + "var_names_dict = ugrid_01.grid_var_names\n", + "x_1 = ugrid_01._ds[var_names_dict['Mesh2_node_x']]\n", + "y_1 = ugrid_01._ds[var_names_dict['Mesh2_node_y']]\n", + "face_nodes_1 = ugrid_01._ds[var_names_dict['Mesh2_face_nodes']]\n", + "n_face_nodes_1 = ugrid_01._ds[var_names_dict['nMesh2_node']]\n", + "\n", + "var_names_dict = ugrid_02.grid_var_names\n", + "x_2 = ugrid_02._ds[var_names_dict['Mesh2_node_x']]\n", + "y_2 = ugrid_02._ds[var_names_dict['Mesh2_node_y']]\n", + "face_nodes_2 = ugrid_02._ds[var_names_dict['Mesh2_face_nodes']]\n", + "n_face_nodes_2 = ugrid_02._ds[var_names_dict['nMesh2_node']]" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "collapsed": false, + "pycharm": { + "name": "#%% md\n" + } + }, + "source": [ + "Please note, for instance, we accessed the actual variable `mesh_node_x`\n", + "of `ugrid_02` via indexing the dictionary with the standardized name\n", + "`Mesh2_node_x`, likewise in `ugrid_01`." + ] + }, + { + "cell_type": "markdown", + "metadata": { + "collapsed": false, + "pycharm": { + "name": "#%% md\n" + } + }, + "source": [ + "## 3. UXarray's Standardized UGRID Names\n", + "The last way of accessing grid topology attributes and variables is to use\n", + "the standardized UGRID namings provided by UXarray. This method still\n", + "utilizes the dictionary, `ds_var_names`, under the hood to return a\n", + "reference to the variable or attribute that is stored withing\n", + "`UXarray.Grid._ds`.\n", + "\n", + "This eliminates the need to remember the original variable names and\n", + "needing to index through the `ds_var_names` dictionary. Because of this,\n", + "we find this as the most convenient approach." + ] + }, + { + "cell_type": "code", + "execution_count": 22, + "metadata": { + "collapsed": false, + "pycharm": { + "name": "#%%\n" + } + }, + "outputs": [], + "source": [ + "x_1 = ugrid_01.Mesh2_node_x\n", + "y_1 = ugrid_01.Mesh2_node_y\n", + "face_nodes_1 = ugrid_01.Mesh2_face_nodes\n", + "n_face_nodes_1 = ugrid_01.nMesh2_node\n", + "\n", + "x_2 = ugrid_02.Mesh2_node_x\n", + "y_2 = ugrid_02.Mesh2_node_y\n", + "face_nodes_2 = ugrid_02.Mesh2_face_nodes\n", + "n_face_nodes_2 = ugrid_02.nMesh2_node" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "collapsed": false, + "pycharm": { + "name": "#%% md\n" + } + }, + "source": [ + "In conclusion, there are three ways of accessing the grid attributes and\n", + "variables. Even though the UXarray developers recommend using the\n", + "standardized UGRID names method, users can find each various pros/cons with\n", + "each of these access ways." + ] + } + ], + "metadata": { + "kernelspec": { + "display_name": "Python 3", + "language": "python", + "name": "python3" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 2 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython2", + "version": "2.7.6" + } + }, + "nbformat": 4, + "nbformat_minor": 0 +} diff --git a/docs/examples/003-area-calc.ipynb b/docs/examples/003-area-calc.ipynb index c60cebc8d..275043bc3 100644 --- a/docs/examples/003-area-calc.ipynb +++ b/docs/examples/003-area-calc.ipynb @@ -2,7 +2,12 @@ "cells": [ { "cell_type": "markdown", - "metadata": {}, + "metadata": { + "collapsed": false, + "pycharm": { + "name": "#%% md\n" + } + }, "source": [ "# Face Area Calculations\n", "\n", @@ -23,7 +28,12 @@ }, { "cell_type": "markdown", - "metadata": {}, + "metadata": { + "collapsed": false, + "pycharm": { + "name": "#%% md\n" + } + }, "source": [ "We will be using the `outCSne30.ug` grid file, which is encoded in the UGRID convention." ] @@ -37,8 +47,8 @@ "start_time": "2023-06-09T10:04:28.705400Z" }, "collapsed": false, - "jupyter": { - "outputs_hidden": false + "pycharm": { + "name": "#%%\n" } }, "outputs": [], @@ -180,7 +190,11 @@ }, { "cell_type": "markdown", - "metadata": {}, + "metadata": { + "pycharm": { + "name": "#%% md\n" + } + }, "source": [ "For the result above, notice that the area is slightly different than the first calculation we made.\n", "\n", @@ -267,7 +281,11 @@ }, { "cell_type": "markdown", - "metadata": {}, + "metadata": { + "pycharm": { + "name": "#%% md\n" + } + }, "source": [ "Now we compare the values with actual know value and report error for each of the three cases above." ] @@ -401,7 +419,11 @@ }, { "cell_type": "markdown", - "metadata": {}, + "metadata": { + "pycharm": { + "name": "#%% md\n" + } + }, "source": [ "Additionally, if you are using a unit other than meters, you can update the units as follows" ] @@ -564,7 +586,11 @@ }, { "cell_type": "markdown", - "metadata": {}, + "metadata": { + "pycharm": { + "name": "#%% md\n" + } + }, "source": [ "## 6. Area Calculation without Grid Object\n", "\n", diff --git a/docs/examples/005-multiprecision-usage.ipynb b/docs/examples/005-multiprecision-usage.ipynb new file mode 100644 index 000000000..e412cd9a7 --- /dev/null +++ b/docs/examples/005-multiprecision-usage.ipynb @@ -0,0 +1,262 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "metadata": { + "pycharm": { + "name": "#%% md\n" + } + }, + "source": [ + "# Multiprecision Arithmetic UXArray\n", + "In this notebook we will explore the use of the UXArray which supports multiprecision arithmetic." + ] + }, + { + "cell_type": "code", + "execution_count": 1, + "metadata": { + "pycharm": { + "name": "#%%\n" + } + }, + "outputs": [], + "source": [ + "import numpy as np\n", + "import uxarray as ux\n", + "import gmpy2\n", + "from gmpy2 import mpfr, mpz" + ] + }, + { + "cell_type": "markdown", + "source": [ + "# Using GMPY2.mpfr type for Multiprecision Arithmetic\n", + "The distance between any two nodes is less than 1e-17 (the precision of the nodes is 60 bits, which is larger than the the python `float` type precision of 53 bits)" + ], + "metadata": { + "collapsed": false, + "pycharm": { + "name": "#%% md\n" + } + } + }, + { + "cell_type": "code", + "execution_count": null, + "outputs": [], + "source": [ + "# Generate the nodes on the unit sphere, and the nodes are extremely close to each other: the edge length\n", + "# between any two nodes is less than 1e-17 (the precision of the nodes is 60 bits but we will use 65 bits for\n", + "# better demonstration)\n", + "nodes = []\n", + "for i in range(4):\n", + " theta = gmpy2.div(i * 2 * gmpy2.const_pi(precision=65),\n", + " 4) # Angle in radians\n", + " x = gmpy2.cos(theta)\n", + " y = gmpy2.sin(theta)\n", + " z = mpfr(\n", + " '0') # All nodes will have z-coordinate as 0 on the unit sphere\n", + " nodes.append([x, y, z])\n", + "\n", + "# Generate the face nodes connectivity\n", + "dumb_nodes = [\n", + " ux.INT_FILL_VALUE_MPZ, ux.INT_FILL_VALUE_MPZ, ux.INT_FILL_VALUE_MPZ\n", + "]\n", + "\n", + "face_nodes_connectivity = np.array([\n", + " np.array([nodes[0], nodes[1], nodes[2], dumb_nodes], dtype=object),\n", + " np.array([nodes[0], nodes[2], nodes[3], dumb_nodes], dtype=object),\n", + " np.array([nodes[0], nodes[3], nodes[1], dumb_nodes], dtype=object),\n", + " np.array([nodes[1], nodes[2], nodes[3], nodes[0]], dtype=object)\n", + "], dtype=object)" + ], + "metadata": { + "collapsed": false, + "pycharm": { + "name": "#%%\n" + } + } + }, + { + "cell_type": "markdown", + "source": [ + "Now we can see this topology has 4 faces and 4 nodes. Three of the facess have 3 nodes and one face has 4 nodes. So we also use the\n", + "fill value `ux.INT_FILL_VALUE_MPZ` to fill the fourth node of the first three faces." + ], + "metadata": { + "collapsed": false, + "pycharm": { + "name": "#%% md\n" + } + } + }, + { + "cell_type": "code", + "execution_count": null, + "outputs": [], + "source": [ + "face_nodes_connectivity" + ], + "metadata": { + "collapsed": false, + "pycharm": { + "name": "#%%\n" + } + } + }, + { + "cell_type": "markdown", + "source": [ + "## Construct the Grid object with the multiprecision mode on" + ], + "metadata": { + "collapsed": false, + "pycharm": { + "name": "#%% md\n" + } + } + }, + { + "cell_type": "code", + "execution_count": null, + "outputs": [], + "source": [ + "# Create the grid\n", + "vgrid = ux.Grid(face_nodes_connectivity,\n", + " multi_precision=True,\n", + " precision=65,\n", + " vertices=True,\n", + " islatlon=False,\n", + " concave=False)" + ], + "metadata": { + "collapsed": false, + "pycharm": { + "name": "#%%\n" + } + } + }, + { + "cell_type": "markdown", + "source": [ + "## Check the Grid Information\n", + "We can see that even though the input nodes are extremely close to each other, the grid is still valid and not degenerated:\n", + "The face number is still 4 and the node number is still 4." + ], + "metadata": { + "collapsed": false, + "pycharm": { + "name": "#%% md\n" + } + } + }, + { + "cell_type": "code", + "execution_count": null, + "outputs": [], + "source": [ + "print(vgrid.nMesh2_face)\n", + "print(vgrid.nMesh2_node)" + ], + "metadata": { + "collapsed": false, + "pycharm": { + "name": "#%%\n" + } + } + }, + { + "cell_type": "markdown", + "source": [ + "# Using python string type for Multiprecision Arithmetic\n", + "You can also input the coordinates as string. The string will be converted to the mpfr type." + ], + "metadata": { + "collapsed": false, + "pycharm": { + "name": "#%% md\n" + } + } + }, + { + "cell_type": "code", + "execution_count": null, + "outputs": [], + "source": [ + "# Provide nodes that are extremely closed to each other\n", + "face_nodes_connectivity = [\n", + " ['1.000000000000000000', '1.000000000000000000001', '1.0000000000000000000012'],\n", + " ['1.000000000000000001', '1.000000000000000000002', '1.0000000000000000000013'],\n", + " ['1.00000000000000000001', '1.00000000000000000000', '1.0000000000000000000013']]\n", + "\n", + "uxgrid = ux.Grid(face_nodes_connectivity, multi_precision=True, precision= 100,\n", + " vertices=True,\n", + " islatlon=False,\n", + " concave=False)" + ], + "metadata": { + "collapsed": false, + "pycharm": { + "name": "#%%\n" + } + } + }, + { + "cell_type": "markdown", + "source": [ + "## Check the Grid Information\n", + "Although the input nodes are extremely close to each other, the grid is still valid and not degenerated:\n", + "The face number is still 1 and the node number is still 3." + ], + "metadata": { + "collapsed": false, + "pycharm": { + "name": "#%% md\n" + } + } + }, + { + "cell_type": "code", + "execution_count": null, + "outputs": [], + "source": [ + "print(uxgrid.nMesh2_face)\n", + "print(uxgrid.nMesh2_node)" + ], + "metadata": { + "collapsed": false, + "pycharm": { + "name": "#%%\n" + } + } + } + ], + "metadata": { + "kernelspec": { + "display_name": "Python 3.9.13 ('uxarray-docs')", + "language": "python", + "name": "python3" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 3 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython3", + "version": "3.9.13" + }, + "orig_nbformat": 4, + "vscode": { + "interpreter": { + "hash": "5df6a94ef43b22c97c277f6a9dce3b546a324b9de41c6d8e82aecf8eafd14442" + } + } + }, + "nbformat": 4, + "nbformat_minor": 2 +} diff --git a/docs/gallery.yml b/docs/gallery.yml index 2939034a5..eb04cd44d 100644 --- a/docs/gallery.yml +++ b/docs/gallery.yml @@ -13,3 +13,7 @@ - title: Working with MPAS Grids path: examples/004-working-with-mpas-grids.ipynb thumbnail: _static/examples/mpas/c-grid.png + +- title: Multiprecision Arithmetic Grid + path: examples/005-multiprecision-usage.ipynb + thumbnail: _static/thumbnails/default.svg diff --git a/docs/user_api/index.rst b/docs/user_api/index.rst index 9fe5ddbb9..19c4d4b3b 100644 --- a/docs/user_api/index.rst +++ b/docs/user_api/index.rst @@ -49,8 +49,32 @@ Methods .. autosummary:: :toctree: _autosummary - UxDataset.info - UxDataset.integrate + helpers.calculate_face_area + helpers.calculate_spherical_triangle_jacobian + helpers.calculate_spherical_triangle_jacobian_barycentric + helpers.get_all_face_area_from_coords + helpers.get_gauss_quadratureDG + helpers.get_tri_quadratureDG + helpers.grid_center_lat_lon + helpers.parse_grid_type + helpers.node_xyz_to_lonlat_rad + helpers.node_lonlat_rad_to_xyz + helpers.normalize_in_place + helpers.close_face_nodes + UxDataset.info + UxDataset.integrate + +Multi-Precision Helper Functions +------------------------------- +.. autosummary:: + :toctree: _autosummary + + multi_precision_helpers.convert_to_multiprecision + multi_precision_helpers.unique_coordinates_multiprecision + multi_precision_helpers.decimal_digits_to_precision_bits + multi_precision_helpers.precision_bits_to_decimal_digits + + UxDataArray @@ -155,3 +179,14 @@ Helpers node_lonlat_rad_to_xyz normalize_in_place parse_grid_type + + +Multi-Precision +------------------------------- +.. autosummary:: + :toctree: _autosummary + + multi_precision_helpers.convert_to_multiprecision + multi_precision_helpers.unique_coordinates_multiprecision + multi_precision_helpers.decimal_digits_to_precision_bits + multi_precision_helpers.precision_bits_to_decimal_digits diff --git a/test/test_grid.py b/test/test_grid.py index e5fc8da0b..39ff52aac 100644 --- a/test/test_grid.py +++ b/test/test_grid.py @@ -2,11 +2,16 @@ import numpy as np import numpy.testing as nt import xarray as xr +import gmpy2 +from gmpy2 import mpfr from unittest import TestCase from pathlib import Path +import timeit import uxarray as ux +import uxarray.multi_precision_helpers as mph +import numpy.testing as nt try: import constants @@ -146,6 +151,106 @@ def test_init_verts(self): assert (vgrid.nMesh2_node == 6) vgrid.encode_as("ugrid") + def test_init_verts_from_string(self): + # Provide nodes that are extremely closed to each other + face_nodes_connectivity = [[ + '1.000000000000000000', '1.000000000000000000001', + '1.0000000000000000000012' + ], + [ + '1.000000000000000001', + '1.000000000000000000002', + '1.0000000000000000000013' + ], + [ + '1.00000000000000000001', + '1.00000000000000000000', + '1.0000000000000000000013' + ]] + + uxgrid = ux.Grid(face_nodes_connectivity, + multi_precision=True, + precision=100, + vertices=True, + islatlon=False, + concave=False) + self.assertEqual(uxgrid.nMesh2_face, 1) + self.assertEqual(uxgrid.nMesh2_node, 3) + + def test_init_verts_multi_precision_triangle(self): + nodes = [] + # Generate the nodes on the unit sphere, and the nodes are extremely close to each other: the edge length + # between any two nodes is less than 1e-17 (the precision of the nodes is 60 bits but we will use 65 bits for + # testing) + for i in range(4): + theta = gmpy2.div(i * 2 * gmpy2.const_pi(precision=65), + 4) # Angle in radians + x = gmpy2.cos(theta) + y = gmpy2.sin(theta) + z = mpfr( + 0) # All nodes will have z-coordinate as 0 on the unit sphere + nodes.append([x, y, z]) + + # Generate the face nodes connectivity + face_nodes_connectivity = np.array([ + np.array([nodes[0], nodes[1], nodes[2]], dtype=object), + np.array([nodes[0], nodes[2], nodes[3]], dtype=object), + np.array([nodes[0], nodes[3], nodes[1]], dtype=object), + np.array([nodes[1], nodes[2], nodes[3]], dtype=object) + ]) + + # Create the grid + vgrid = ux.Grid(face_nodes_connectivity, + multi_precision=True, + precision=65, + vertices=True, + islatlon=False, + concave=False) + assert (vgrid.source_grid == "From vertices") + assert (vgrid.nMesh2_face == 4) + assert (vgrid.nMesh2_node == 4) + + def test_init_verts_multi_precision_fill_values(self): + # Set the desired precision + + # Generate the nodes on the unit sphere, and the nodes are extremely close to each other: the edge length + # between any two nodes is less than 1e-17 (the precision of the nodes is 60 bits but we will use 65 bits for + # testing) + + nodes = [] + for i in range(4): + theta = gmpy2.div(i * 2 * gmpy2.const_pi(precision=65), + 4) # Angle in radians + x = gmpy2.cos(theta) + y = gmpy2.sin(theta) + z = mpfr('0') + nodes.append([x, y, z]) + + # Generate the face nodes connectivity + dumb_nodes = [ + ux.INT_FILL_VALUE_MPZ, ux.INT_FILL_VALUE_MPZ, ux.INT_FILL_VALUE_MPZ + ] + face_nodes_connectivity = np.array([ + np.array([nodes[0], nodes[1], nodes[2], dumb_nodes], dtype=object), + np.array([nodes[0], nodes[2], nodes[3], dumb_nodes], dtype=object), + np.array([nodes[0], nodes[3], nodes[1], dumb_nodes], dtype=object), + np.array([nodes[1], nodes[2], nodes[3], nodes[0]], dtype=object) + ], + dtype=object) + # Create the grid + vgrid = ux.Grid(face_nodes_connectivity, + multi_precision=True, + precision=65, + vertices=True, + islatlon=False, + concave=False) + assert (vgrid.source_grid == "From vertices") + assert (vgrid.nMesh2_face == 4) + assert (vgrid.nMesh2_node == 4) + + # Test the all numbers in the vgird.ds["Mesh2_face_nodes"] are less than 4 + assert (np.all(vgrid._ds["Mesh2_face_nodes"] < 4)) + def test_init_verts_different_input_datatype(self): """Create a uxarray grid from multiple face vertices with different datatypes(ndarray, list, tuple) and saves a ugrid file. @@ -315,7 +420,7 @@ def test_calculate_total_face_area_triangle(self): islatlon=False, isconcave=False) - #calculate area + # calculate area area_gaussian = grid_verts.calculate_total_face_area( quadrature_rule="gaussian", order=5) nt.assert_almost_equal(area_gaussian, constants.TRI_AREA, decimal=3) @@ -400,7 +505,7 @@ def test_verts_calc_area(self): class TestPopulateCoordinates(TestCase): - def test_populate_cartesian_xyz_coord(self): + def test_populate_cartesian_xyz_coord_float(self): # The following testcases are generated through the matlab cart2sph/sph2cart functions # These points correspond to the eight vertices of a cube. lon_deg = [ @@ -427,8 +532,7 @@ def test_populate_cartesian_xyz_coord(self): ] verts_degree = np.stack((lon_deg, lat_deg), axis=1) - - vgrid = ux.open_grid(verts_degree, islatlon=False) + vgrid = ux.Grid([verts_degree], islatlon=True, vertices=True) vgrid._populate_cartesian_xyz_coord() for i in range(0, vgrid.nMesh2_node): @@ -442,7 +546,50 @@ def test_populate_cartesian_xyz_coord(self): cart_z[i], decimal=12) - def test_populate_lonlat_coord(self): + def test_populate_cartesian_xyz_coord_mpfr(self): + lon_deg = [ + '45.0001052295749', '45.0001052295749', '314.9998947704251', + '314.9998947704251' + ] + lat_deg = [ + '35.2655522903022', '-35.2655522903022', '35.2655522903022', + '-35.2655522903022' + ] + + cart_x = [ + '0.577340924821405', '0.577340924821405', '0.577340924821405', + '0.577340924821405' + ] + + cart_y = [ + '0.577343045516932', '0.577343045516932', '-0.577343045516932', + '-0.577343045516932' + ] + + cart_z = [ + '0.577366836872017', '-0.577366836872017', '0.577366836872017', + '-0.577366836872017' + ] + + verts_degree = np.stack((lon_deg, lat_deg), axis=1) + vgrid = ux.Grid([verts_degree], + multi_precision=True, + precision=64, + islatlon=True, + vertices=True) + vgrid._populate_cartesian_xyz_coord() + for i in range(0, vgrid.nMesh2_node): + nt.assert_almost_equal(vgrid._ds["Mesh2_node_cart_x"].values[i], + mpfr(cart_x[i]), + decimal=14) + nt.assert_almost_equal(vgrid._ds["Mesh2_node_cart_y"].values[i], + mpfr(cart_y[i]), + decimal=14) + nt.assert_almost_equal(vgrid._ds["Mesh2_node_cart_z"].values[i], + mpfr(cart_z[i]), + decimal=14) + + def test_populate_lonlat_coord_float(self): # The following testcases are generated through the matlab cart2sph/sph2cart functions # These points correspond to the 4 vertexes on a cube. @@ -481,6 +628,45 @@ def test_populate_lonlat_coord(self): lat_deg[i], decimal=12) + def test_populate_lonlat_coord_mpfr(self): + lon_deg = [ + '45.0001052295749', '45.0001052295749', '314.9998947704251', + '314.9998947704251' + ] + lat_deg = [ + '35.2655522903022', '-35.2655522903022', '35.2655522903022', + '-35.2655522903022' + ] + + cart_x = [ + '0.577340924821405', '0.577340924821405', '0.577340924821405', + '0.577340924821405' + ] + + cart_y = [ + '0.577343045516932', '0.577343045516932', '-0.577343045516932', + '-0.577343045516932' + ] + + cart_z = [ + '0.577366836872017', '-0.577366836872017', '0.577366836872017', + '-0.577366836872017' + ] + + verts_cart = np.stack((cart_x, cart_y, cart_z), axis=1) + vgrid = ux.Grid([verts_cart], + multi_precision=True, + precision=64, + islatlon=False) + vgrid._populate_lonlat_coord() + for i in range(0, vgrid.nMesh2_node): + nt.assert_almost_equal(vgrid._ds["Mesh2_node_x"].values[i], + mpfr(lon_deg[i]), + decimal=13) + nt.assert_almost_equal(vgrid._ds["Mesh2_node_y"].values[i], + mpfr(lat_deg[i]), + decimal=13) + class TestConnectivity(TestCase): mpas_filepath = current_path / "meshfiles" / "mpas" / "QU" / "mesh.QU.1920km.151026.nc" @@ -734,28 +920,6 @@ def test_build_face_edges_connectivity(self): np.array_equal(reverted_mesh2_edge_nodes[i], original_face_nodes_connectivity[i])) - def test_build_face_edges_connectivity_mpas(self): - tgrid = ux.open_grid(self.mpas_filepath) - - mesh2_face_nodes = tgrid._ds["Mesh2_face_nodes"] - - tgrid._build_face_edges_connectivity() - mesh2_face_edges = tgrid._ds.Mesh2_face_edges - mesh2_edge_nodes = tgrid._ds.Mesh2_edge_nodes - - # Assert if the mesh2_face_edges sizes are correct. - self.assertEqual(mesh2_face_edges.sizes["nMesh2_face"], - mesh2_face_nodes.sizes["nMesh2_face"]) - self.assertEqual(mesh2_face_edges.sizes["nMaxMesh2_face_edges"], - mesh2_face_nodes.sizes["nMaxMesh2_face_nodes"]) - - # Assert if the mesh2_edge_nodes sizes are correct. - # Euler formular for determining the edge numbers: n_face = n_edges - n_nodes + 2 - num_edges = mesh2_face_edges.sizes["nMesh2_face"] + tgrid._ds[ - "Mesh2_node_x"].sizes["nMesh2_node"] - 2 - size = mesh2_edge_nodes.sizes["nMesh2_edge"] - self.assertEqual(mesh2_edge_nodes.sizes["nMesh2_edge"], num_edges) - def test_build_face_edges_connectivity_fillvalues(self): verts = [ self.f0_deg, self.f1_deg, self.f2_deg, self.f3_deg, self.f4_deg, diff --git a/test/test_helpers.py b/test/test_helpers.py index fe2caee30..ea886b332 100644 --- a/test/test_helpers.py +++ b/test/test_helpers.py @@ -3,14 +3,21 @@ import numpy.testing as nt import random import xarray as xr - +import gmpy2 +from gmpy2 import mpfr, mpz +import mpmath from unittest import TestCase from pathlib import Path +import time import uxarray as ux from uxarray.utils.helpers import _replace_fill_values from uxarray.utils.constants import INT_DTYPE, INT_FILL_VALUE +from uxarray.multi_precision_helpers import convert_to_multiprecision, set_global_precision, \ + decimal_digits_to_precision_bits +from uxarray.utils.helpers import _replace_fill_values +from uxarray.utils.constants import INT_DTYPE, INT_FILL_VALUE try: import constants @@ -109,6 +116,26 @@ def test_normalize_in_place(self): self.assertLessEqual(np.absolute(np.sqrt(x * x + y * y + z * z) - 1), err_tolerance) + # Multiprecision test for places=19 + precision = decimal_digits_to_precision_bits(19) + set_global_precision(precision) + [x_mpfr, y_mpfr, + z_mpfr] = convert_to_multiprecision(np.array([ + '1.0000000000000000001', '0.0000000000000000009', + '0.0000000000000000001' + ]), + precision=precision) + normalized = ux.utils.helpers.normalize_in_place( + [x_mpfr, y_mpfr, z_mpfr]) + # Calculate the sum of squares using gmpy2.fsum() + sum_of_squares = gmpy2.fsum( + [gmpy2.square(value) for value in normalized]) + abs = gmpy2.mul(gmpy2.reldiff(mpfr('1.0'), sum_of_squares), mpfr('1.0')) + self.assertAlmostEqual(abs, 0, places=19) + + # Reset global precision to default + set_global_precision() + def test_node_xyz_to_lonlat_rad(self): [x, y, z] = ux.utils.helpers.normalize_in_place([ random.uniform(-1, 1), @@ -124,6 +151,36 @@ def test_node_xyz_to_lonlat_rad(self): self.assertLessEqual(np.absolute(new_y - y), err_tolerance) self.assertLessEqual(np.absolute(new_z - z), err_tolerance) + # Multiprecision test for places=19 + precision = decimal_digits_to_precision_bits(19) + set_global_precision(precision) + # Assign 1 at the 21th decimal place, which is beyond the precision of 19 decimal places + [x_mpfr, y_mpfr, + z_mpfr] = convert_to_multiprecision(np.array([ + '1.000000000000000000001', '0.000000000000000000001', + '0.000000000000000000001' + ]), + precision=precision) + [lon_mpfr, lat_mpfr + ] = ux.utils.helpers.node_xyz_to_lonlat_rad([x_mpfr, y_mpfr, z_mpfr]) + self.assertAlmostEqual(lon_mpfr, 0, places=19) + self.assertAlmostEqual(lat_mpfr, 0, places=19) + + # Remove 1 at the 21th decimal place, and the total digit place are 19. the results should be perfectly equal to [0,0] + [x_mpfr, y_mpfr, + z_mpfr] = convert_to_multiprecision(np.array([ + '1.0000000000000000000', '0.0000000000000000000', + '0.0000000000000000000' + ]), + precision=precision) + [lon_mpfr, lat_mpfr + ] = ux.utils.helpers.node_xyz_to_lonlat_rad([x_mpfr, y_mpfr, z_mpfr]) + self.assertTrue(gmpy2.cmp(lon_mpfr, mpfr('0')) == 0) + self.assertTrue(gmpy2.cmp(lat_mpfr, mpfr('0')) == 0) + + # Reset global precision to default + set_global_precision() + def test_node_latlon_rad_to_xyz(self): [lon, lat] = [ random.uniform(0, 2 * np.pi), @@ -137,6 +194,129 @@ def test_node_latlon_rad_to_xyz(self): self.assertLessEqual(np.absolute(new_lon - lon), err_tolerance) self.assertLessEqual(np.absolute(new_lat - lat), err_tolerance) + # Multiprecision test for places=19 + precision = decimal_digits_to_precision_bits(19) + set_global_precision(precision) + # Assign 1 at the 21th decimal place, which is beyond the precision of 19 decimal places + [lon_mpfr, lat_mpfr] = convert_to_multiprecision(np.array( + ['0.000000000000000000001', '0.000000000000000000001']), + precision=precision) + [x_mpfr, y_mpfr, + z_mpfr] = ux.utils.helpers.node_lonlat_rad_to_xyz([lon_mpfr, lat_mpfr]) + self.assertAlmostEqual(x_mpfr, 1, places=19) + self.assertAlmostEqual(y_mpfr, 0, places=19) + self.assertAlmostEqual(z_mpfr, 0, places=19) + + # Remove 1 at the 21th decimal place, and the total digit place are 19. the results should be perfectly equal to [1,0,0] + [lon_mpfr, lat_mpfr] = convert_to_multiprecision(np.array( + ['0.0000000000000000000', '0.0000000000000000000']), + precision=precision) + [x_mpfr, y_mpfr, + z_mpfr] = ux.utils.helpers.node_lonlat_rad_to_xyz([lon_mpfr, lat_mpfr]) + self.assertTrue(gmpy2.cmp(x_mpfr, mpfr('1')) == 0) + self.assertTrue(gmpy2.cmp(y_mpfr, mpfr('0')) == 0) + self.assertTrue(gmpy2.cmp(z_mpfr, mpfr('0')) == 0) + + # Reset global precision to default + set_global_precision() + + def test_precise_coordinates_conversion(self): + # Multiprecision test for places=19, And we set the global precision places to 20 + # Repeat the conversion between latitude and longitude and xyz for 1000 times + # And see if the results are the same + precision = decimal_digits_to_precision_bits(20) + set_global_precision(precision) + + # The initial coordinates + [init_x, init_y, init_z] = ux.utils.helpers.normalize_in_place([ + mpfr('0.12345678910111213149'), + mpfr('0.92345678910111213149'), + mpfr('1.72345678910111213149') + ]) + new_x = init_x + new_y = init_y + new_z = init_z + for iter in range(1000): + [new_lon, new_lat + ] = ux.utils.helpers.node_xyz_to_lonlat_rad([new_x, new_y, new_z]) + [new_x, new_y, new_z + ] = ux.utils.helpers.node_lonlat_rad_to_xyz([new_lon, new_lat]) + self.assertAlmostEqual(new_x, init_x, places=19) + self.assertAlmostEqual(new_y, init_y, places=19) + self.assertAlmostEqual(new_z, init_z, places=19) + + # Test for the longitude and latitude conversion + # The initial coordinates + [init_lon, + init_lat] = [mpfr('1.4000332309896247'), + mpfr('1.190289949682531')] + # Reset global precision to default + new_lat = init_lat + new_lon = init_lon + for iter in range(1000): + [new_x, new_y, new_z + ] = ux.utils.helpers.node_lonlat_rad_to_xyz([new_lon, new_lat]) + [new_lon, new_lat + ] = ux.utils.helpers.node_xyz_to_lonlat_rad([new_x, new_y, new_z]) + self.assertAlmostEqual(new_lon, init_lon, places=19) + self.assertAlmostEqual(new_lat, init_lat, places=19) + + # Reset global precision to default + set_global_precision() + + def test_coordinates_conversion_accumulate_error(self): + # Get the accumulated error of each function call + ux.multi_precision_helpers.set_global_precision(64) + run_time = 100 + print("\n") + + # Using the float number + new_lon = 122.987654321098765 + new_lat = 36.123456789012345 + + start_time = time.time() + for iter in range(run_time): + [new_x, new_y, new_z] = ux.utils.helpers.node_lonlat_rad_to_xyz( + np.deg2rad([new_lon, new_lat])) + [new_lon, new_lat + ] = ux.utils.helpers.node_xyz_to_lonlat_rad([new_x, new_y, new_z]) + [new_lon, new_lat] = np.rad2deg([new_lon, new_lat]) + + end_time = time.time() + diff_lat = mpfr(str(new_lat)) - mpfr('36.123456789012345') + diff_lon = mpfr(str(new_lon)) - mpfr('122.987654321098765') + print("The floating point longitude accumulated error is: " + + str(diff_lon) + "and the latitude accumulated " + "error is: " + str(diff_lat)) + print("The floating point Execution time: ", end_time - start_time, + " seconds") + + # Get the accumulated error of each function call + print("\n") + + # Using the float number + + [init_lon, init_lat + ] = [mpfr('122.987654321098765', 64), + mpfr('36.123456789012345', 64)] + new_lon = mpfr('122.987654321098765', 64) + new_lat = mpfr('36.123456789012345', 64) + start_time = time.time() + + for iter in range(run_time): + [new_x, new_y, new_z] = ux.utils.helpers.node_lonlat_rad_to_xyz( + [gmpy2.radians(val) for val in [new_lon, new_lat]]) + [new_lon, new_lat + ] = ux.utils.helpers.node_xyz_to_lonlat_rad([new_x, new_y, new_z]) + [new_lon, + new_lat] = [gmpy2.degrees(val) for val in [new_lon, new_lat]] + end_time = time.time() + diff_lat = new_lat - init_lat + diff_lon = new_lon - init_lon + print("The mpfr longitude accumulated error is: " + str(diff_lon) + + " and the latitude accumulated error is: " + str(diff_lat)) + print("The mpfr Execution time: ", end_time - start_time, " seconds") + class TestConstants(TestCase): # DTYPE as set in constants.py diff --git a/test/test_multi_precision_helpers.py b/test/test_multi_precision_helpers.py new file mode 100644 index 000000000..48ecb3a6f --- /dev/null +++ b/test/test_multi_precision_helpers.py @@ -0,0 +1,291 @@ +import os +import numpy as np +import gmpy2 +from gmpy2 import mpfr, mpz + +from unittest import TestCase +from pathlib import Path + +import uxarray as ux + +from uxarray.utils.constants import INT_DTYPE, INT_FILL_VALUE, FLOAT_PRECISION_BITS, INT_FILL_VALUE_MPZ +from uxarray.multi_precision_helpers import convert_to_multiprecision, unique_coordinates_multiprecision, precision_bits_to_decimal_digits, decimal_digits_to_precision_bits +import math + +try: + import constants +except ImportError: + from . import constants + +# Data files +current_path = Path(os.path.dirname(os.path.realpath(__file__))) + +exodus = current_path / "meshfiles" / "exodus" / "outCSne8" / "outCSne8.g" +ne8 = current_path / 'meshfiles' / "scrip" / "outCSne8" / 'outCSne8.nc' +err_tolerance = 1.0e-12 + + +class TestMultiPrecision(TestCase): + + def test_convert_to_multiprecision(self): + """Tests if the convert_to_mpfr() helper function converts a numpy + array to a numpy array of the correct dtype.""" + # test different datatypes for face_nodes + test_precision = 64 + f0_deg = np.array([ + np.array([120, -20]), + np.array([130, -10]), + np.array([120, 0]), + np.array([105, 0]), + np.array([95, -10]), + np.array([105, -20]) + ]) + f1_deg = np.array([ + np.array([120, 0]), + np.array([120, 10]), + np.array([115, 0]), + np.array([ux.INT_FILL_VALUE, ux.INT_FILL_VALUE]), + np.array([ux.INT_FILL_VALUE, ux.INT_FILL_VALUE]), + np.array([ux.INT_FILL_VALUE, ux.INT_FILL_VALUE]) + ]) + f2_deg = np.array([ + np.array([115, 0]), + np.array([120, 10]), + np.array([100, 10]), + np.array([105, 0]), + np.array([ux.INT_FILL_VALUE, ux.INT_FILL_VALUE]), + np.array([ux.INT_FILL_VALUE, ux.INT_FILL_VALUE]) + ]) + f3_deg = np.array([ + np.array([95, -10]), + np.array([105, 0]), + np.array([95, 30]), + np.array([80, 30]), + np.array([70, 0]), + np.array([75, -10]) + ]) + f4_deg = np.array([ + np.array([65, -20]), + np.array([75, -10]), + np.array([70, 0]), + np.array([55, 0]), + np.array([45, -10]), + np.array([55, -20]) + ]) + f5_deg = np.array([ + np.array([70, 0]), + np.array([80, 30]), + np.array([70, 30]), + np.array([60, 0]), + np.array([ux.INT_FILL_VALUE, ux.INT_FILL_VALUE]), + np.array([ux.INT_FILL_VALUE, ux.INT_FILL_VALUE]) + ]) + f6_deg = np.array([ + np.array([60, 0]), + np.array([70, 30]), + np.array([40, 30]), + np.array([45, 0]), + np.array([ux.INT_FILL_VALUE, ux.INT_FILL_VALUE]), + np.array([ux.INT_FILL_VALUE, ux.INT_FILL_VALUE]) + ]) + + verts = np.array( + [f0_deg, f1_deg, f2_deg, f3_deg, f4_deg, f5_deg, f6_deg]) + + verts_mpfr = convert_to_multiprecision(verts, + str_mode=False, + precision=test_precision) + + # Test if every object in the verts_mpfr array is of type mpfr + for i in range(verts_mpfr.shape[0]): + for j in range(verts_mpfr.shape[1]): + for idx, val in enumerate(verts_mpfr[i, j]): + if type(val) != mpz: + self.assertEqual(val.precision, test_precision) + # Then compare the values between verts and verts_mpfr up to the 53 bits of precision + self.assertAlmostEqual(verts[i, j][idx], + val, + places=FLOAT_PRECISION_BITS) + + else: + self.assertTrue(gmpy2.cmp(val, INT_FILL_VALUE_MPZ) == 0) + + def test_unique_coordinates_multiprecision_normal_case(self): + """The input cartesian coordinates represents 8 vertices on a cube 7. + + ---------6. + + /| /| + / | / | + 3---------2 | + | | | | + | 4------|--5 + | / | / + |/ |/ + 0---------1 + """ + + test_precision = 64 + cart_x = [ + 0.577340924821405, 0.577340924821405, 0.577340924821405, + 0.577340924821405, -0.577345166204668, -0.577345166204668, + -0.577345166204668, -0.577345166204668 + ] + cart_y = [ + 0.577343045516932, 0.577343045516932, -0.577343045516932, + -0.577343045516932, 0.577338804118089, 0.577338804118089, + -0.577338804118089, -0.577338804118089 + ] + cart_z = [ + 0.577366836872017, -0.577366836872017, 0.577366836872017, + -0.577366836872017, 0.577366836872017, -0.577366836872017, + 0.577366836872017, -0.577366836872017 + ] + + # The order of the vertexes is irrelevant, the following indexing is just for forming a face matrix + face_vertices = [ + [0, 1, 2, 3], # front face + [1, 5, 6, 2], # right face + [5, 4, 7, 6], # back face + [4, 0, 3, 7], # left face + [3, 2, 6, 7], # top face + [4, 5, 1, 0] # bottom face + ] + + # Pack the cart_x/y/z into the face matrix using the index from face_vertices + faces_coords = [] + for face in face_vertices: + face_coords = [] + for vertex_index in face: + x, y, z = cart_x[vertex_index], cart_y[vertex_index], cart_z[ + vertex_index] + face_coords.append([x, y, z]) + faces_coords.append(face_coords) + + # Now consturct the grid using the faces_coords + verts_cart = np.array(faces_coords) + verts_cart_mpfr = convert_to_multiprecision(verts_cart, + str_mode=False, + precision=test_precision) + verts_cart_mpfr_unique, unique_inverse = unique_coordinates_multiprecision( + verts_cart_mpfr.reshape(-1, verts_cart_mpfr.shape[-1]), + precision=test_precision) + recovered_verts_cart_mpfr = verts_cart_mpfr_unique[unique_inverse] + + # Compare the recovered verts_cart_mpfr with the original verts_cart_mpfr.reshape(-1, verts_cart_mpfr.shape[-1]) + expected = verts_cart_mpfr.reshape(-1, verts_cart_mpfr.shape[-1]) + for i in range(recovered_verts_cart_mpfr.shape[0]): + for j in range(recovered_verts_cart_mpfr.shape[1]): + self.assertEqual(recovered_verts_cart_mpfr[i, j].precision, + test_precision) + # Then compare the values between verts and verts_mpfr up to the 53 bits of precision + self.assertAlmostEqual(expected[i, j], + recovered_verts_cart_mpfr[i, j], + places=FLOAT_PRECISION_BITS) + + def test_unique_coordinates_multiprecision_extreme_case(self): + # Construct coordinates that are extremely closed to each other + coord1 = [ + gmpy2.mpfr('1.23456789'), + gmpy2.mpfr('2.34567890'), + gmpy2.mpfr('3.45678901') + ] + coord2 = [ + gmpy2.mpfr('1.23456789'), + gmpy2.mpfr('2.34567890'), + gmpy2.mpfr('3.45678901') + ] + # Convert the coordinates to string format + coord1_str = ["{0:.17f}".format(coord) for coord in coord1] + coord2_str = ["{0:.17f}".format(coord) for coord in coord2] + + # Create the final coordinates with the differing 64th bit + coord1_final = [coord + '0' for coord in coord1_str] + coord2_final = [coord + '1' for coord in coord2_str] + verts = np.array([coord1_final, coord2_final]) + bit_precision = decimal_digits_to_precision_bits(18) + verts_60 = convert_to_multiprecision(verts, + str_mode=True, + precision=bit_precision) + verts_57 = convert_to_multiprecision(verts, + str_mode=True, + precision=bit_precision - 1) + + verts_cart_mpfr_unique_64, unique_inverse_64 = unique_coordinates_multiprecision( + verts_60, precision=bit_precision) + verts_cart_mpfr_unique_57, unique_inverse_57 = unique_coordinates_multiprecision( + verts_57, precision=bit_precision - 1) + self.assertTrue(len(verts_cart_mpfr_unique_64) == 2) + self.assertTrue(len(verts_cart_mpfr_unique_57) == 1) + + coord1 = [ + gmpy2.mpfr('1.23456789'), + gmpy2.mpfr('2.34567890'), + gmpy2.mpfr('3.45678901') + ] + coord2 = [ + gmpy2.mpfr('1.23456789'), + gmpy2.mpfr('2.34567890'), + gmpy2.mpfr('3.45678901') + ] + # Convert the coordinates to string format + precision_bits = 200 + decimal_digits = precision_bits_to_decimal_digits(precision_bits - 1) + check = decimal_digits_to_precision_bits(decimal_digits) + + format_str = "{0:." + str(decimal_digits) + "f}" + coord1_str = [format_str.format(coord) for coord in coord1] + coord2_str = [format_str.format(coord) for coord in coord2] + + # Create the final coordinates with the differing 64th bit + coord1_final = [coord + '0' for coord in coord1_str] + coord2_final = [coord + '1' for coord in coord2_str] + verts = np.array([coord1_final, coord2_final]) + + verts_200 = convert_to_multiprecision(verts, + str_mode=True, + precision=precision_bits) + verts_199 = convert_to_multiprecision(verts, + str_mode=True, + precision=check - 1) + + verts_cart_mpfr_unique_200, unique_inverse_200 = unique_coordinates_multiprecision( + verts_200, precision=precision_bits) + verts_cart_mpfr_unique_199, unique_inverse_199 = unique_coordinates_multiprecision( + verts_199, precision=check - 1) + self.assertTrue(len(verts_cart_mpfr_unique_200) == 2) + self.assertTrue(len(verts_cart_mpfr_unique_199) == 1) + + def test_unique_coordinates_multiprecision_filled_values(self): + # The mpft_unique function should be able to handle filled values + coord1 = [ + gmpy2.mpfr('1.23456789'), + gmpy2.mpfr('2.34567890'), + gmpy2.mpfr('3.45678901') + ] + coord2 = [ + gmpy2.mpfr('1.23456789'), + gmpy2.mpfr('2.34567890'), + gmpy2.mpfr('3.45678901') + ] + # Convert the coordinates to string format + precision_bits = 200 + decimal_digits = precision_bits_to_decimal_digits(precision_bits - 1) + + format_str = "{0:." + str(decimal_digits) + "f}" + coord1_str = [format_str.format(coord) for coord in coord1] + coord2_str = [format_str.format(coord) for coord in coord2] + + # Create the final coordinates with the differing 64th bit + coord1_final = [coord + '0' for coord in coord1_str] + coord2_final = [coord + '1' for coord in coord2_str] + verts = np.array([ + coord1_final, coord2_final, + ['INT_FILL_VALUE', 'INT_FILL_VALUE', 'INT_FILL_VALUE'] + ]) + verts_200 = convert_to_multiprecision(verts, + str_mode=True, + precision=precision_bits) + unique_coords, unique_inverse = unique_coordinates_multiprecision( + verts_200, precision=precision_bits) + self.assertTrue(len(unique_coords) == 3) diff --git a/uxarray/__init__.py b/uxarray/__init__.py index b9aedd249..1ef44bb8b 100644 --- a/uxarray/__init__.py +++ b/uxarray/__init__.py @@ -1,5 +1,5 @@ from uxarray.utils.helpers import * -from uxarray.utils.constants import (INT_DTYPE, INT_FILL_VALUE) +from uxarray.utils.constants import * from uxarray.core.grid import Grid from uxarray.core.api import (open_grid, open_dataset, open_mfdataset) from uxarray.core.dataarray import UxDataArray diff --git a/uxarray/core/grid.py b/uxarray/core/grid.py index 76f76d0db..b67d556db 100644 --- a/uxarray/core/grid.py +++ b/uxarray/core/grid.py @@ -1,8 +1,12 @@ """uxarray.core.grid module.""" import xarray as xr import numpy as np +import gmpy2 +from gmpy2 import mpfr, mpz # reader and writer imports +from uxarray.multi_precision_helpers import convert_to_multiprecision, unique_coordinates_multiprecision, \ + precision_bits_to_decimal_digits, decimal_digits_to_precision_bits from uxarray.io._exodus import _read_exodus, _encode_exodus from uxarray.io._mpas import _read_mpas from uxarray.io._ugrid import _read_ugrid, _encode_ugrid @@ -11,7 +15,7 @@ from uxarray.utils.helpers import (get_all_face_area_from_coords, parse_grid_type, node_xyz_to_lonlat_rad, node_lonlat_rad_to_xyz, close_face_nodes) -from uxarray.utils.constants import INT_DTYPE, INT_FILL_VALUE +from uxarray.utils.constants import INT_DTYPE, INT_FILL_VALUE, INT_FILL_VALUE_MPZ, FLOAT_PRECISION_BITS class Grid: @@ -57,13 +61,54 @@ class Grid: >>> uxds = ux.open_dataset("filename.g") """ - def __init__(self, input_obj, **kwargs): + def __init__(self, + dataset, + multi_precision=False, + precision=FLOAT_PRECISION_BITS, + **kwargs): + """Initialize grid variables, decide if loading happens via file, verts + or gridspec. + + Parameters + ---------- + dataset : xarray.Dataset, ndarray, list, tuple, required + Input xarray.Dataset or vertex coordinates that form one face. + multi_precision : bool, optional + Specify if want to use multi-precision (mpfr) for all grid operations (default: False). + Each node coordinates will be stored as the gmpy2.mpfr/mpz(filled value) type + instead of the python ``float``/``int`` type. + Note: when using multi-precision, all grid operations are done in multi-precision using the GMPY2 library. And + the run time will significantly increase. To gain the full benefit of multi-precision, it is recommended to input + the grid coordinates in string format (eg. '0.1', '0.2', etc) and use string 'INT_FILL_VALUE' + to specify the fill value. + precision : int, optional + Specify the precision of the grid coordinates (default: FLOAT_PRECISION_BITS (python `float` bit precision)). + Other Parameters + ---------------- + islatlon : bool, optional + Specify if the grid is lat/lon based + concave: bool, optional + Specify if this grid has concave elements (internal checks for this are possible) + gridspec: bool, optional + Specifies gridspec + mesh_type: str, optional + Specify the mesh file type, eg. exo, ugrid, shp, mpas, etc + use_dual: bool, optional + Specify whether to use the primal (use_dual=False) or dual (use_dual=True) mesh if the file type is mpas + Raises + ------ + RuntimeError + If specified file not found + """ # initialize internal variable names self.__init_grid_var_names__() # initialize face_area variable self._face_areas = None + # initialize the multi-precision flag + self._multi_precision = multi_precision + self._precision = precision # TODO: fix when adding/exercising gridspec # unpack kwargs @@ -76,29 +121,53 @@ def __init__(self, input_obj, **kwargs): setattr(self, key, kwargs.get(key, None)) # check if initializing from verts: - if isinstance(input_obj, (list, tuple, np.ndarray)): - input_obj = np.asarray(input_obj) - self.mesh_type = "From vertices" + if isinstance(dataset, (list, tuple, np.ndarray)): + dataset = np.asarray(dataset) + + # Pre-process dataset if the multi-precision flag is set + if self._multi_precision: + # If the input are floats + if dataset.dtype == float: + dataset = convert_to_multiprecision( + dataset, str_mode=False, precision=self._precision) + # If the input are strings + elif np.all([ + np.issubdtype(type(element), np.str_) + for element in dataset.ravel() + ]): + dataset = convert_to_multiprecision( + dataset, str_mode=True, precision=self._precision) + else: + # If the input are not floats or strings or gmpy2.mpfr/mpz, raise an error + if ~np.any( + np.vectorize(lambda x: isinstance( + x, (gmpy2.mpfr, gmpy2.mpz)))(dataset.ravel())): + raise ValueError( + 'The input array should be either floats, strings, or gmpy2.mpfr/mpz' + ) # grid with multiple faces - if input_obj.ndim == 3: - self.__from_vert__(input_obj) + if dataset.ndim == 3: + self.__from_vert__(dataset) self.source_grid = "From vertices" # grid with a single face - elif input_obj.ndim == 2: - input_obj = np.array([input_obj]) - self.__from_vert__(input_obj) + elif dataset.ndim == 2: + dataset = np.array([dataset]) + self.__from_vert__(dataset, + multi_precision=multi_precision, + precision=precision) self.source_grid = "From vertices" else: raise RuntimeError( - f"Invalid Input Dimension: {input_obj.ndim}. Expected dimension should be " + f"Invalid Input Dimension: {dataset.ndim}. Expected dimension should be " f"3: [nMesh2_face, nMesh2_node, Two/Three] or 2 when only " f"one face is passed in.") # check if initializing from string # TODO: re-add gridspec initialization when implemented - elif isinstance(input_obj, xr.Dataset): - self.mesh_type = parse_grid_type(input_obj) - self.__from_ds__(dataset=input_obj) + elif isinstance(dataset, xr.Dataset): + self.mesh_type = parse_grid_type(dataset) + #TODO: Support multiprecision for __from_ds__ + self.__from_ds__(dataset=dataset) else: raise RuntimeError("Dataset is not a valid input type.") @@ -126,7 +195,10 @@ def __init_grid_var_names__(self): "nMaxMesh2_face_nodes": "nMaxMesh2_face_nodes" } - def __from_vert__(self, dataset): + def __from_vert__(self, + dataset, + multi_precision=False, + precision=FLOAT_PRECISION_BITS): """Create a grid with faces constructed from vertices specified by the given argument. @@ -167,19 +239,71 @@ def __from_vert__(self, dataset): z_coord = x_coord * 0.0 # Identify unique vertices and their indices - unique_verts, indices = np.unique(dataset.reshape( - -1, dataset.shape[-1]), - axis=0, - return_inverse=True) + if self._multi_precision: + # Check if the input_array is in th mpfr type + try: + # Flatten the input_array_mpfr to a 1D array so that we can check the type of each element + input_array_mpfr_copy = np.ravel(dataset) + for i in range(len(input_array_mpfr_copy)): + if type(input_array_mpfr_copy[i]) != gmpy2.mpfr and type( + input_array_mpfr_copy[i]) != gmpy2.mpz: + dataset = convert_to_multiprecision( + dataset, precision=self._precision) + except Exception as e: + raise e + + unique_verts, indices = unique_coordinates_multiprecision( + dataset.reshape(-1, dataset.shape[-1]), + precision=self._precision) + + else: + unique_verts, indices = np.unique(dataset.reshape( + -1, dataset.shape[-1]), + axis=0, + return_inverse=True) # Nodes index that contain a fill value - fill_value_mask = np.logical_or(unique_verts[:, 0] == INT_FILL_VALUE, - unique_verts[:, 1] == INT_FILL_VALUE) - if dataset[0][0].size > 2: + if self._multi_precision: + # Perform element-wise comparison using gmpy.cmp() + fill_value_mask = np.logical_or( + np.array([ + gmpy2.cmp(x, INT_FILL_VALUE_MPZ) == 0 + for x in unique_verts[:, 0] + ], + dtype=bool), + np.array([ + gmpy2.cmp(x, INT_FILL_VALUE_MPZ) == 0 + for x in unique_verts[:, 1] + ], + dtype=bool)) + + if dataset[0][0].size > 2: + fill_value_mask = np.logical_or( + np.array([ + gmpy2.cmp(x, INT_FILL_VALUE_MPZ) == 0 + for x in unique_verts[:, 0] + ], + dtype=bool), + np.array([ + gmpy2.cmp(x, INT_FILL_VALUE_MPZ) == 0 + for x in unique_verts[:, 1] + ], + dtype=bool), + np.array([ + gmpy2.cmp(x, INT_FILL_VALUE_MPZ) == 0 + for x in unique_verts[:, 2] + ], + dtype=bool)) + + else: fill_value_mask = np.logical_or( unique_verts[:, 0] == INT_FILL_VALUE, - unique_verts[:, 1] == INT_FILL_VALUE, - unique_verts[:, 2] == INT_FILL_VALUE) + unique_verts[:, 1] == INT_FILL_VALUE) + if dataset[0][0].size > 2: + fill_value_mask = np.logical_or( + unique_verts[:, 0] == INT_FILL_VALUE, + unique_verts[:, 1] == INT_FILL_VALUE, + unique_verts[:, 2] == INT_FILL_VALUE) # Get the indices of all the False values in fill_value_mask false_indices = np.where(fill_value_mask == True)[0] @@ -192,9 +316,20 @@ def __from_vert__(self, dataset): unique_verts = np.delete(unique_verts, false_indices, axis=0) # Update indices accordingly - for i, idx in enumerate(false_indices): - indices[indices == idx] = INT_FILL_VALUE - indices[(indices > idx) & (indices != INT_FILL_VALUE)] -= 1 + if self._multi_precision: + + for i, idx in enumerate(false_indices): + indices[indices == idx] = INT_FILL_VALUE_MPZ + for j in range(indices.size): + if gmpy2.cmp(mpz( + indices[j]), mpz(idx)) > 0 and gmpy2.cmp( + mpz(indices[j]), INT_FILL_VALUE_MPZ) != 0: + indices[j] -= 1 + + else: + for i, idx in enumerate(false_indices): + indices[indices == idx] = INT_FILL_VALUE + indices[(indices > idx) & (indices != INT_FILL_VALUE)] -= 1 # Create coordinate DataArrays self._ds["Mesh2_node_x"] = xr.DataArray(data=unique_verts[:, 0], @@ -215,14 +350,25 @@ def __from_vert__(self, dataset): # Create connectivity array using indices of unique vertices connectivity = indices.reshape(dataset.shape[:-1]) - self._ds["Mesh2_face_nodes"] = xr.DataArray( - data=xr.DataArray(connectivity).astype(INT_DTYPE), - dims=["nMesh2_face", "nMaxMesh2_face_nodes"], - attrs={ - "cf_role": "face_node_connectivity", - "_FillValue": INT_FILL_VALUE, - "start_index": 0 - }) + + if self._multi_precision: + self._ds["Mesh2_face_nodes"] = xr.DataArray( + data=xr.DataArray(connectivity).astype(INT_DTYPE), + dims=["nMesh2_face", "nMaxMesh2_face_nodes"], + attrs={ + "cf_role": "face_node_connectivity", + "_FillValue": INT_FILL_VALUE_MPZ, + "start_index": 0 + }) + else: + self._ds["Mesh2_face_nodes"] = xr.DataArray( + data=xr.DataArray(connectivity).astype(INT_DTYPE), + dims=["nMesh2_face", "nMaxMesh2_face_nodes"], + attrs={ + "cf_role": "face_node_connectivity", + "_FillValue": INT_FILL_VALUE, + "start_index": 0 + }) # load mesh from a file def __from_ds__(self, dataset): @@ -870,9 +1016,15 @@ def _populate_cartesian_xyz_coord(self): if "Mesh2_node_cart_x" in self._ds.keys(): return - # check for units and create Mesh2_node_cart_x/y/z set to self._ds - nodes_lon_rad = np.deg2rad(self.Mesh2_node_x.values) - nodes_lat_rad = np.deg2rad(self.Mesh2_node_y.values) + # check for units and create Mesh2_node_cart_x/y/z set to self.ds + if self._multi_precision: + nodes_lon_rad = np.array( + [gmpy2.radians(x) for x in self.Mesh2_node_x.values]) + nodes_lat_rad = np.array( + [gmpy2.radians(y) for y in self.Mesh2_node_y.values]) + else: + nodes_lon_rad = np.deg2rad(self.Mesh2_node_x.values) + nodes_lat_rad = np.deg2rad(self.Mesh2_node_y.values) nodes_rad = np.stack((nodes_lon_rad, nodes_lat_rad), axis=1) nodes_cart = np.asarray( list(map(node_lonlat_rad_to_xyz, list(nodes_rad)))) @@ -955,7 +1107,12 @@ def _populate_lonlat_coord(self): self._ds["Mesh2_node_z"].values), axis=1).tolist() nodes_rad = list(map(node_xyz_to_lonlat_rad, nodes_cart)) - nodes_degree = np.rad2deg(nodes_rad) + + if self._multi_precision: + nodes_degree = np.array( + [[gmpy2.degrees(x), gmpy2.degrees(y)] for x, y in nodes_rad]) + else: + nodes_degree = np.rad2deg(nodes_rad) self._ds["Mesh2_node_x"] = xr.DataArray( data=nodes_degree[:, 0], dims=["nMesh2_node"], diff --git a/uxarray/multi_precision_helpers.py b/uxarray/multi_precision_helpers.py new file mode 100644 index 000000000..66af712c6 --- /dev/null +++ b/uxarray/multi_precision_helpers.py @@ -0,0 +1,220 @@ +import gmpy2 +from gmpy2 import mpfr, mpz +import numpy as np +import math +from uxarray.utils.constants import FLOAT_PRECISION_BITS, INT_FILL_VALUE_MPZ + + +def set_global_precision(global_precision=FLOAT_PRECISION_BITS): + """Set the global precision of the mpfr numbers. + Important Note: + 1. To avoid arithmetic overflow, the global precision should always be higher than any other precision speicified + in the code. + 2. Modifying the precision by calling this function will modify all following codes running context until + another call to this function. + + Parameters + ---------- + global_precision : int, optional + The global precision of the expected multiprecision float. + The default precision of an mpfr is 53 bits - the same precision as Python’s `float` type. + + Returns + ------- + None + """ + + gmpy2.get_context().precision = global_precision + + +def convert_to_multiprecision(input_array, + str_mode=True, + precision=FLOAT_PRECISION_BITS): + """Convert a numpy array to a list of mpfr numbers. + + The default precision of an mpfr is 53 bits - the same precision as Python’s `float` type. + https://gmpy2.readthedocs.io/en/latest/mpfr.html + If the input array contains fill values INT_FILL_VALUE, the fill values will be converted to INT_FILL_VALUE_MPZ, + which is the multi-precision integer representation of INT_FILL_VALUE. + + Parameters + ---------- + input_array : numpy array, float/string, shape is arbitrary + The input array to be converted to mpfr. The input array should be float or string. If the input array is float, + str_mode should be False. If the input array is string, str_mode should be True. + + str_mode : bool, optional + If True, the input array should be string when passing into the function. + If False, the input array should be float when passing into the function. + str_mode is True by default and is recommended. Because to take advantage of the higher precision provided by + the mpfr type, always pass constants as strings. + precision : int, optional + The precision of the mpfr numbers. The default precision of an mpfr is 53 bits - the same precision as Python’s `float` type. + + Returns + ---------- + mpfr_array : numpy array, mpfr type, shape will be same as the input_array + The output array with mpfr type, which supports correct + rounding, selectable rounding modes, and many trigonometric, exponential, and special functions. A context + manager is used to control precision, rounding modes, and the behavior of exceptions. + + Raises + ---------- + ValueError + The input array should be string when str_mode is True, if not, raise + ValueError('The input array should be string when str_mode is True.') + """ + + # To take advantage of the higher precision provided by the mpfr type, always pass constants as strings. + # https://gmpy2.readthedocs.io/en/latest/mpfr.html + flattened_array = np.ravel(input_array) + mpfr_array = np.array(flattened_array, dtype=object) + if not str_mode: + # Cast the input 2D array to string array + for idx, val in enumerate(flattened_array): + if gmpy2.cmp(mpz(val), INT_FILL_VALUE_MPZ) == 0: + mpfr_array[idx] = INT_FILL_VALUE_MPZ + else: + decimal_digit = precision_bits_to_decimal_digits(precision) + format_str = "{0:+." + str(decimal_digit) + "f}" + val_str = format_str.format(val) + mpfr_array[idx] = mpfr(val_str, precision) + + else: + + if ~np.all([ + np.issubdtype(type(element), np.str_) + for element in flattened_array + ]): + raise ValueError( + 'The input array should be string when str_mode is True.') + # Then convert the input array to mpfr array + for idx, val in enumerate(flattened_array): + if val == "INT_FILL_VALUE": + mpfr_array[idx] = INT_FILL_VALUE_MPZ + else: + mpfr_array[idx] = mpfr(val, precision) + + mpfr_array = mpfr_array.reshape(input_array.shape) + + return mpfr_array + + +def unique_coordinates_multiprecision(input_array_mpfr, + precision=FLOAT_PRECISION_BITS): + """Find the unique coordinates in the input array with mpfr numbers. + + The default precision of an mpfr is 53 bits - the same precision as Python’s `float` type. + It can recognize the fill values INT_FILL_VALUE_MPZ, which is the multi-precision integer representation of + INT_FILL_VALUE. + + Parameters: + ---------- + input_array_mpfr : numpy.ndarray, gmpy2.mpfr type + The input array containing mpfr numbers. + + precision : int, optional + The precision in bits used for the mpfr calculations. Default is FLOAT_PRECISION_BITS. + + Returns: + ------- + unique_arr ： numpy.ndarray, gmpy2.mpfr + Array of unique coordinates in the input array. + + inverse_indices: numpy.ndarray, int + The indices to reconstruct the original array from the unique array. Only provided if return_inverse is True. + + Raises + ---------- + ValueError + The input array should be string when str_mode is True, if not, raise + ValueError('The input array should be string when str_mode is True.') + """ + + # Check if the input_array is in th mpfr type + try: + # Flatten the input_array_mpfr to a 1D array so that we can check the type of each element + input_array_mpfr_copy = np.ravel(input_array_mpfr) + for i in range(len(input_array_mpfr_copy)): + if type(input_array_mpfr_copy[i]) != gmpy2.mpfr and type( + input_array_mpfr_copy[i]) != gmpy2.mpz: + raise ValueError( + 'The input array should be in the mpfr type. You can use convert_to_mpfr() to ' + 'convert the input array to mpfr.') + except Exception as e: + raise e + + unique_arr = [] + inverse_indices = [] + m, n = input_array_mpfr.shape + unique_dict = {} + current_index = 0 + + for i in range(m): + # We only need to check the first element of each row since the elements in the same row are the same type + # (Either mpfr for valid coordinates or INT_FILL_VALUE_MPZ for fill values) + if type(input_array_mpfr[i][0]) == gmpy2.mpfr: + format_string = "{0:+." + str(precision + 1) + "Uf}" + elif type(input_array_mpfr[i][0]) == gmpy2.mpz: + format_string = "{:<+" + str(precision + 1) + "d}" + else: + raise ValueError( + 'The input array should be in the mpfr/mpz type. You can use convert_to_multiprecision() ' + 'to convert the input array to multiprecision format.') + hashable_row = tuple( + format_string.format(x) for x in input_array_mpfr[i]) + + if hashable_row not in unique_dict: + unique_dict[hashable_row] = current_index + unique_arr.append(input_array_mpfr[i]) + inverse_indices.append(current_index) + current_index += 1 + else: + inverse_indices.append(unique_dict[hashable_row]) + + unique_arr = np.array(unique_arr) + inverse_indices = np.array(inverse_indices) + + return unique_arr, inverse_indices + + +def decimal_digits_to_precision_bits(decimal_digits): + """Convert the number of decimal digits to the number of bits of precision. + + Parameters + ---------- + decimal_digits : int + The number of decimal digits of precision + + Returns + ------- + bits : int + The number of bits of precision + """ + bits = math.ceil(decimal_digits * math.log2(10)) + return bits + + +def precision_bits_to_decimal_digits(precision): + """Convert the number of bits of precision to the number of decimal digits. + + Parameters + ---------- + precision : int + The number of bits of precision + + Returns + ------- + decimal_digits : int + The number of decimal digits of precision + """ + # Compute the decimal digit count using gmpy2.log10() + log10_2 = gmpy2.log10(gmpy2.mpfr(2)) # Logarithm base 10 of 2 + log10_precision = gmpy2.div(1, + log10_2) # Logarithm base 10 of the precision + decimal_digits = gmpy2.div(precision, log10_precision) + + # Convert the result to an integer + decimal_digits = int(math.floor(decimal_digits)) + + return decimal_digits diff --git a/uxarray/utils/constants.py b/uxarray/utils/constants.py index 7a8f322f8..60ecebb24 100644 --- a/uxarray/utils/constants.py +++ b/uxarray/utils/constants.py @@ -1,5 +1,9 @@ import numpy as np +from gmpy2 import mpz # numpy indexing code is written for np.intp INT_DTYPE = np.intp INT_FILL_VALUE = np.iinfo(INT_DTYPE).min +INT_FILL_VALUE_MPZ = mpz(str(INT_FILL_VALUE)) +FLOAT_PRECISION_BITS = 53 +ERROR_TOLERANCE = 1.0e-15 diff --git a/uxarray/utils/helpers.py b/uxarray/utils/helpers.py index caf18b77c..cbbfe3235 100644 --- a/uxarray/utils/helpers.py +++ b/uxarray/utils/helpers.py @@ -1,11 +1,13 @@ import numpy as np import xarray as xr -from pathlib import PurePath +import gmpy2 +from gmpy2 import mpfr, mpz +import mpmath from .get_quadratureDG import get_gauss_quadratureDG, get_tri_quadratureDG from numba import njit, config import math -from uxarray.utils.constants import INT_DTYPE, INT_FILL_VALUE +from uxarray.utils.constants import INT_DTYPE, INT_FILL_VALUE, ERROR_TOLERANCE config.DISABLE_JIT = False @@ -96,7 +98,7 @@ def parse_grid_type(dataset): # Calculate the area of all faces. -@njit +# @njit def calculate_face_area(x, y, z, @@ -183,7 +185,7 @@ def calculate_face_area(x, return area -@njit +# @njit def get_all_face_area_from_coords(x, y, z, @@ -253,7 +255,7 @@ def get_all_face_area_from_coords(x, return area -@njit +# @njit def calculate_spherical_triangle_jacobian(node1, node2, node3, dA, dB): """Calculate Jacobian of a spherical triangle. This is a helper function for calculating face area. @@ -329,7 +331,7 @@ def calculate_spherical_triangle_jacobian(node1, node2, node3, dA, dB): return dJacobian -@njit +# @njit def calculate_spherical_triangle_jacobian_barycentric(node1, node2, node3, dA, dB): """Calculate Jacobian of a spherical triangle. This is a helper function @@ -466,19 +468,18 @@ def grid_center_lat_lon(ds): return center_lat, center_lon -@njit def node_lonlat_rad_to_xyz(node_coord): """Helper function to Convert the node coordinate from 2D longitude/latitude to normalized 3D xyz. Parameters ---------- - node: float list + node: list or np.array, python `float` or `gmpy2.mpfr` 2D coordinates[longitude, latitude] in radiance Returns ---------- - float list + node_cartesian: np.array, python `float` or `gmpy2.mpfr` the result array of the unit 3D coordinates [x, y, z] vector where :math:`x^2 + y^2 + z^2 = 1` Raises @@ -486,27 +487,42 @@ def node_lonlat_rad_to_xyz(node_coord): RuntimeError The input array doesn't have the size of 3. """ + if isinstance(node_coord, list): + node_coord = np.array(node_coord) if len(node_coord) != 2: raise RuntimeError( "Input array should have a length of 2: [longitude, latitude]") - lon = node_coord[0] - lat = node_coord[1] - return [np.cos(lon) * np.cos(lat), np.sin(lon) * np.cos(lat), np.sin(lat)] + if np.any( + np.vectorize(lambda x: isinstance(x, (gmpy2.mpfr, gmpy2.mpz)))( + node_coord)): + lon = node_coord[0] + lat = node_coord[1] + return np.array([ + gmpy2.cos(lon) * gmpy2.cos(lat), + gmpy2.sin(lon) * gmpy2.cos(lat), + gmpy2.sin(lat) + ]) + else: + lon = node_coord[0] + lat = node_coord[1] + return np.array( + [np.cos(lon) * np.cos(lat), + np.sin(lon) * np.cos(lat), + np.sin(lat)]) -@njit def node_xyz_to_lonlat_rad(node_coord): """Calculate the latitude and longitude in radiance for a node represented in the [x, y, z] 3D Cartesian coordinates. Parameters ---------- - node_coord: float list + node_coord: np.array or python list of `float` or gmpy2.mpfr 3D Cartesian Coordinates [x, y, z] of the node Returns ---------- - float list + node_coord_rad: np.array of `float` or gmpy2.mpfr the result array of longitude and latitude in radian [longitude_rad, latitude_rad] Raises @@ -514,31 +530,51 @@ def node_xyz_to_lonlat_rad(node_coord): RuntimeError The input array doesn't have the size of 3. """ + if isinstance(node_coord, list): + node_coord = np.array(node_coord) + elif not isinstance(node_coord, np.ndarray): + raise TypeError( + "Input node_coord should be either a numpy array or a list.") if len(node_coord) != 3: raise RuntimeError("Input array should have a length of 3: [x, y, z]") - reference_tolerance = 1.0e-12 - [dx, dy, dz] = normalize_in_place(node_coord) - dx /= np.absolute(dx * dx + dy * dy + dz * dz) - dy /= np.absolute(dx * dx + dy * dy + dz * dz) - dz /= np.absolute(dx * dx + dy * dy + dz * dz) - - if np.absolute(dz) < (1.0 - reference_tolerance): - d_lon_rad = math.atan2(dy, dx) - d_lat_rad = np.arcsin(dz) - - if d_lon_rad < 0.0: - d_lon_rad += 2.0 * np.pi - elif dz > 0.0: - d_lon_rad = 0.0 - d_lat_rad = 0.5 * np.pi + if np.any( + np.vectorize(lambda x: isinstance(x, (gmpy2.mpfr, gmpy2.mpz)))( + node_coord)): + [dx, dy, dz] = normalize_in_place(node_coord) + # The precision is set by the gmpy2.context through the set_global_precision() function + if gmpy2.cmp_abs(dz, mpfr('1.0')): + d_lon_rad = gmpy2.atan2(dy, dx) + d_lat_rad = gmpy2.asin(dz) + + if gmpy2.cmp(d_lon_rad, mpfr('0.0')) < 0: + d_lon_rad += mpfr('2.0') * gmpy2.const_pi() + elif gmpy2.cmp(dz, mpfr('0.0')) > 0: + d_lon_rad = mpfr('0.0') + d_lat_rad = mpfr('0.5') * gmpy2.const_pi() + else: + d_lon_rad = mpfr('0.0') + d_lat_rad = mpfr('-0.5') * gmpy2.const_pi() + + return np.array([d_lon_rad, d_lat_rad]) else: - d_lon_rad = 0.0 - d_lat_rad = -0.5 * np.pi + [dx, dy, dz] = normalize_in_place(node_coord) + + if np.absolute(dz) < (1.0 - ERROR_TOLERANCE): + d_lon_rad = math.atan2(dy, dx) + d_lat_rad = np.arcsin(dz) + + if d_lon_rad < 0.0: + d_lon_rad += 2.0 * np.pi + elif dz > 0.0: + d_lon_rad = 0.0 + d_lat_rad = 0.5 * np.pi + else: + d_lon_rad = 0.0 + d_lat_rad = -0.5 * np.pi - return [d_lon_rad, d_lat_rad] + return np.array([d_lon_rad, d_lat_rad]) -@njit def normalize_in_place(node): """Helper function to project an arbitrary node in 3D coordinates [x, y, z] on the unit sphere. It uses the `np.linalg.norm` internally to calculate @@ -546,12 +582,12 @@ def normalize_in_place(node): Parameters ---------- - node: float list + node: np.array or python list of `float` or gmpy2.mpfr 3D Cartesian Coordinates [x, y, z] Returns ---------- - float list + normalized_node: np. array of `float` or gmpy2.mpfr the result unit vector [x, y, z] where :math:`x^2 + y^2 + z^2 = 1` Raises @@ -559,10 +595,30 @@ def normalize_in_place(node): RuntimeError The input array doesn't have the size of 3. """ - if len(node) != 3: + if isinstance(node, list): + node = np.array(node) + elif not isinstance(node, np.ndarray): + raise TypeError("Input node should be either a numpy array or a list.") + + if node.size != 3: raise RuntimeError("Input array should have a length of 3: [x, y, z]") - return np.array(node) / np.linalg.norm(np.array(node), ord=2) + if np.any( + np.vectorize(lambda x: isinstance(x, (gmpy2.mpfr, gmpy2.mpz)))( + node)): + # Convert the node to mpmath array + norm = gmpy2.sqrt(gmpy2.fsum([gmpy2.square(value) for value in node])) + + # Divide each element of the node by the norm using gmpy2 + normalized_node = [gmpy2.div(element, norm) for element in node] + + # Convert the result back to numpy array or list depending on the input type + if isinstance(node, np.ndarray): + return np.array(normalized_node) + else: + return normalized_node + else: + return node / np.linalg.norm(node, ord=2) def _replace_fill_values(grid_var, original_fill, new_fill, new_dtype=None):