Make the position coordinates consistent across data groups

README.md

@@ -38,7 +38,7 @@ with velocity vectors overlayed as arrows, once again using ``layers``:
 
 ```python
 ind = np.argmax(data["hydro"]["density"])
-center = data["amr"]["xyz"][ind.values]
+center = data["amr"]["position"][ind.values]
 osyris.plane({"data": data["hydro"]["density"], "norm": "log"}, # layer 1
              {"data": data["hydro"]["velocity"], "mode": "vec"}, # layer 2
              dx=2000 * osyris.units("au"),
@@ -49,7 +49,7 @@ osyris.plane({"data": data["hydro"]["density"], "norm": "log"}, # layer 1
 
 ## Have a problem or need a new feature?
 
-- Bugs reports or feature requests should be submitted by opening an [issue](https://github.com/osyris-project/osyris/issues)
+- Bug reports or feature requests should be submitted by opening an [issue](https://github.com/osyris-project/osyris/issues)
 - For general discussions or questions about how to do something with `osyris`, start a new [discussion](https://github.com/osyris-project/osyris/discussions)
 
 ## Logo credit

docs/basics.ipynb

@@ -164,7 +164,7 @@
    "outputs": [],
    "source": [
     "ind = np.argmax(data[\"hydro\"][\"density\"])\n",
-    "center = data[\"amr\"][\"xyz\"][ind.values]\n",
+    "center = data[\"amr\"][\"position\"][ind.values]\n",
     "center"
    ]
   },

docs/data_structures.ipynb

@@ -302,7 +302,7 @@
    "outputs": [],
    "source": [
     "ind = np.argmax(data[\"hydro\"][\"density\"])\n",
-    "center = data[\"amr\"][\"xyz\"][ind.values]\n",
+    "center = data[\"amr\"][\"position\"][ind.values]\n",
     "center"
    ]
   },
@@ -459,7 +459,7 @@
  ],
  "metadata": {
   "kernelspec": {
-   "display_name": "Python 3",
+   "display_name": "Python 3 (ipykernel)",
    "language": "python",
    "name": "python3"
   },
@@ -473,9 +473,9 @@
    "name": "python",
    "nbconvert_exporter": "python",
    "pygments_lexer": "ipython3",
-   "version": "3.9.1"
+   "version": "3.9.5"
   }
  },
  "nbformat": 4,
- "nbformat_minor": 2
+ "nbformat_minor": 4
 }

docs/loading_data.ipynb

@@ -108,7 +108,7 @@
    "source": [
     "data = osyris.Dataset(8, scale=\"au\", path=path).load(\n",
     "    select={\"hydro\": {\"density\": lambda d : d > 1.0e-16 * osyris.units('g/cm**3')},\n",
-    "            \"amr\": {\"xyz_x\": lambda x : x > 1500. * osyris.units('au')}})"
+    "            \"amr\": {\"x\": lambda x : x > 1500. * osyris.units('au')}})"
    ]
   },
   {
@@ -119,7 +119,7 @@
    "source": [
     "osyris.map({\"data\": data[\"hydro\"][\"density\"], \"norm\": \"log\"},\n",
     "           dx=1000 * osyris.units('au'),\n",
-    "           origin=data[\"amr\"][\"xyz\"][np.argmax(data[\"hydro\"][\"density\"]).values],\n",
+    "           origin=data[\"amr\"][\"position\"][np.argmax(data[\"hydro\"][\"density\"]).values],\n",
     "           direction='z')"
    ]
   },
@@ -252,7 +252,7 @@
    "source": [
     "osyris.map({\"data\": data[\"hydro\"][\"density\"], \"norm\": \"log\"},\n",
     "           dx=2000 * osyris.units(\"au\"),\n",
-    "           origin=data[\"amr\"][\"xyz\"][np.argmax(data[\"hydro\"][\"density\"]).values],\n",
+    "           origin=data[\"amr\"][\"position\"][np.argmax(data[\"hydro\"][\"density\"]).values],\n",
     "           direction='z')"
    ]
   },
@@ -307,7 +307,7 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "center = data['sink']['xyz'][0:1]\n",
+    "center = data['sink']['position'][0:1]\n",
     "dx = 200 * osyris.units('au')"
    ]
   },
@@ -329,9 +329,9 @@
    "outputs": [],
    "source": [
     "data.load(\n",
-    "    select={\"amr\": {\"x\": lambda x : np.logical_and(x > center.x - dx, x < center.x + dx),\n",
-    "                    \"y\": lambda y : np.logical_and(y > center.y - dx, y < center.y + dx),\n",
-    "                    \"z\": lambda z : np.logical_and(z > center.z - dx, z < center.z + dx)}\n",
+    "    select={\"amr\": {\"x\": lambda x : (x > center.x - dx) & (x < center.x + dx),\n",
+    "                    \"y\": lambda y : (y > center.y - dx) & (y < center.y + dx),\n",
+    "                    \"z\": lambda z : (z > center.z - dx) & (z < center.z + dx)}\n",
     "           })"
    ]
   },
@@ -351,7 +351,8 @@
    "outputs": [],
    "source": [
     "osyris.map(data[\"hydro\"][\"density\"],\n",
-    "           {\"data\": data[\"sink\"][\"xyz\"], \"mode\": \"scatter\", \"c\": \"white\"},\n",
+    "           {\"data\": data[\"sink\"][\"position\"], \"mode\": \"scatter\", \"c\": \"white\",\n",
+    "            \"s\": 20. * osyris.units(\"au\"), \"alpha\": 0.7},\n",
     "           norm='log', direction=\"z\", origin=center[0])"
    ]
   }

docs/plotting_1d_2d.ipynb

@@ -58,7 +58,7 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "osyris.plot(data1d['amr']['xyz_x'], data1d['hydro']['density'])"
+    "osyris.plot(data1d['amr']['position_x'], data1d['hydro']['density'])"
    ]
   },
   {
@@ -79,7 +79,7 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "osyris.plot(data1d['amr']['xyz_x'], data1d['hydro']['density'],\n",
+    "osyris.plot(data1d['amr']['position_x'], data1d['hydro']['density'],\n",
     "           marker='o', ls='None', color='red')"
    ]
   },
@@ -99,7 +99,7 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "osyris.plot(data1d['amr']['xyz_x'], data1d['hydro']['density'],\n",
+    "osyris.plot(data1d['amr']['position_x'], data1d['hydro']['density'],\n",
     "           marker='o', ls='None', logx=True, logy=True)"
    ]
   },
@@ -123,7 +123,7 @@
    "outputs": [],
    "source": [
     "try:\n",
-    "    osyris.plot(data1d['amr']['xyz_x'],\n",
+    "    osyris.plot(data1d['amr']['position_x'],\n",
     "                data1d['hydro']['density'],\n",
     "                data1d['hydro']['pressure'])\n",
     "except Exception as e:\n",
@@ -144,7 +144,7 @@
    "outputs": [],
    "source": [
     "data1d['hydro']['new_density'] = data1d['hydro']['density'] * 2\n",
-    "osyris.plot(data1d['amr']['xyz_x'],\n",
+    "osyris.plot(data1d['amr']['position_x'],\n",
     "            data1d['hydro']['density'],\n",
     "            data1d['hydro']['new_density'])"
    ]
@@ -169,8 +169,8 @@
    "source": [
     "old1d = osyris.Dataset(1, scale=\"cm\", path=path).load()\n",
     "\n",
-    "osyris.plot({'x': data1d['amr']['xyz_x'], 'y': data1d['hydro']['density']},\n",
-    "            {'x': old1d['amr']['xyz_x'], 'y': old1d['hydro']['density']},\n",
+    "osyris.plot({'x': data1d['amr']['position_x'], 'y': data1d['hydro']['density']},\n",
+    "            {'x': old1d['amr']['position_x'], 'y': old1d['hydro']['density']},\n",
     "            marker='o', ls='None')"
    ]
   },
@@ -195,13 +195,13 @@
    "source": [
     "fig, ax = plt.subplots(2, 2, figsize=(12, 9))\n",
     "\n",
-    "osyris.plot(data1d['amr']['xyz_x'], data1d['hydro']['density'],\n",
+    "osyris.plot(data1d['amr']['position_x'], data1d['hydro']['density'],\n",
     "            marker='o', ls='None', ax=ax[0, 0])\n",
-    "osyris.plot(data1d['amr']['xyz_x'], data1d['hydro']['velocity_x'],\n",
+    "osyris.plot(data1d['amr']['position_x'], data1d['hydro']['velocity_x'],\n",
     "            color='green', ax=ax[0, 1])\n",
-    "osyris.plot(data1d['amr']['xyz_x'], data1d['hydro']['pressure'],\n",
+    "osyris.plot(data1d['amr']['position_x'], data1d['hydro']['pressure'],\n",
     "            color='red', ax=ax[1, 0])\n",
-    "osyris.plot(data1d['amr']['xyz_x'], data1d['amr']['level'],\n",
+    "osyris.plot(data1d['amr']['position_x'], data1d['amr']['level'],\n",
     "            color='black', ls='dotted', ax=ax[1, 1])"
    ]
   },

docs/plotting_maps.ipynb

@@ -25,7 +25,7 @@
     "path = \"osyrisdata/starformation\"\n",
     "data = osyris.Dataset(8, scale=\"au\", path=path).load()\n",
     "ind = np.argmax(data[\"hydro\"][\"density\"])\n",
-    "center = data[\"amr\"][\"xyz\"][ind.values]"
+    "center = data[\"amr\"][\"position\"][ind.values]"
    ]
   },
   {
@@ -243,8 +243,7 @@
    "outputs": [],
    "source": [
     "osyris.map(data[\"hydro\"][\"velocity\"], mode=\"lic\",\n",
-    "           color=data[\"hydro\"][\"density\"],\n",
-    "           norm=\"log\", length=60,\n",
+    "           color=data[\"hydro\"][\"density\"], norm=\"log\",\n",
     "           dx=2000 * osyris.units(\"au\"),\n",
     "           origin=center, direction=\"z\")"
    ]
@@ -488,7 +487,7 @@
    "outputs": [],
    "source": [
     "osyris.map({\"data\": data[\"hydro\"][\"density\"], \"norm\": \"log\"}, # layer 1\n",
-    "           {\"data\": data[\"sink\"][\"xyz\"], \"mode\": \"scatter\", \"c\": \"white\"}, # layer 2\n",
+    "           {\"data\": data[\"sink\"][\"position\"], \"mode\": \"scatter\", \"c\": \"white\"}, # layer 2\n",
     "           dx=2000 * osyris.units(\"au\"),\n",
     "           origin=center,\n",
     "           direction=\"z\")"

docs/plotting_thick_maps.ipynb

@@ -49,7 +49,7 @@
     "path = \"osyrisdata/starformation\"\n",
     "data = osyris.Dataset(8, scale=\"au\", path=path).load()\n",
     "ind = np.argmax(data[\"hydro\"][\"density\"])\n",
-    "center = data[\"amr\"][\"xyz\"][ind.values]"
+    "center = data[\"amr\"][\"position\"][ind.values]"
    ]
   },
   {

docs/recipes.ipynb

@@ -94,7 +94,7 @@
    "source": [
     "# Find center\n",
     "ind = np.argmax(data8[\"hydro\"][\"density\"])\n",
-    "center = data8[\"amr\"][\"xyz\"][ind.values]\n",
+    "center = data8[\"amr\"][\"position\"][ind.values]\n",
     "\n",
     "dx = 2000 * osyris.units(\"au\")\n",
     "# Extract density slices by copying data into structures\n",
@@ -138,14 +138,14 @@
    "outputs": [],
    "source": [
     "# Re-center cell coordinates according to origin\n",
-    "data8['amr']['xyz_new'] = data8['amr']['xyz'] - center\n",
+    "data8['amr']['pos_new'] = data8['amr']['position'] - center\n",
     "\n",
     "# Create figure\n",
     "fig, ax = plt.subplots()\n",
     "\n",
     "# Make scatter plot as radial profile\n",
     "step = 100\n",
-    "osyris.scatter(data8['amr']['xyz_new'][::step],\n",
+    "osyris.scatter(data8['amr']['pos_new'][::step],\n",
     "               data8['hydro']['density'][::step],\n",
     "               c='grey', edgecolors='None', loglog=True, ax=ax)\n",
     "\n",
@@ -159,8 +159,8 @@
     "midpoints = 0.5 * (edges[1:] + edges[:-1])\n",
     "\n",
     "# Bin the data in radial bins\n",
-    "z0, _ = np.histogram(np.log10(data8['amr']['xyz_new'].norm.values), bins=edges)\n",
-    "z1, _ = np.histogram(np.log10(data8['amr']['xyz_new'].norm.values), bins=edges,\n",
+    "z0, _ = np.histogram(np.log10(data8['amr']['pos_new'].norm.values), bins=edges)\n",
+    "z1, _ = np.histogram(np.log10(data8['amr']['pos_new'].norm.values), bins=edges,\n",
     "                     weights=data8['hydro']['density'].values)\n",
     "rho_mean = z1/z0\n",
     "\n",

src/osyris/core/array.py

@@ -18,7 +18,7 @@ def _comparison_operator(lhs, rhs, op):
 
 
 class Array:
-    def __init__(self, values=None, unit=None, parent=None, name=""):
+    def __init__(self, values=0, unit=None, parent=None, name=""):
 
         if isinstance(values, np.ndarray):
             self._array = values

src/osyris/core/tools.py

@@ -44,21 +44,6 @@ def to_bin_edges(x):
     return np.append(np.insert(centers, 0, left), right)
 
 
-def perpendicular_vector(v):
-    """
-    Compute a vector perpendicular to the input vector
-    """
-
-    # x = y = z = 0 is not an acceptable solution
-    if v[0] == v[1] == v[2] == 0:
-        raise ValueError("zero-vector")
-
-    if v[2] == 0:
-        return [-v[1], v[0], 0]
-    else:
-        return [1.0, 1.0, -1.0 * (v[0] + v[1]) / v[2]]
-
-
 def value_to_string(val, precision=3):
     """
     Convert a number to a human readable string.

src/osyris/io/amr.py

@@ -16,7 +16,7 @@ def initialize(self, meta, select):
         self.initialized = False
 
         descriptor = {"level": "i", "cpu": "i", "dx": "d"}
-        descriptor.update({"xyz_{}".format(c): "d" for c in "xyz"[:meta["ndim"]]})
+        descriptor.update({f"position_{c}": "d" for c in "xyz"[:meta["ndim"]]})
 
         self.descriptor_to_variables(descriptor=descriptor, meta=meta, select=select)
 
@@ -152,7 +152,7 @@ def read_variables(self, ncache, ind, ilevel, cpuid, info):
         if self.variables["level"]["read"]:
             self.variables["level"]["buffer"]._array[:ncache, ind] = ilevel + 1
         for n in range(info["ndim"]):
-            key = "xyz_" + "xyz"[n]
+            key = "position_" + "xyz"[n]
             if self.variables[key]["read"]:
                 self.variables[key]["buffer"]._array[:ncache, ind] = (
                     self.xg[:ncache, n] + self.xcent[ind, n] - self.meta["xbound"][n]

src/osyris/io/hilbert.py

@@ -173,15 +173,15 @@ def hilbert_cpu_list(meta, scaling, select, infofile):
                         unit=1.0 * scaling.units)
     new_bbox = False
     for c in "xyz":
-        if c in select:
+        key = f"position_{c}"
+        if key in select:
             new_bbox = True
-            func_test = select[c](xyz_centers)
+            func_test = select[key](xyz_centers)
             inds = np.argwhere(func_test).ravel()
             start = xyz_centers[inds.min()] - (half_dxmin * scaling.units)
             end = xyz_centers[inds.max()] + (half_dxmin * scaling.units)
             bounding_box["{}min".format(c)] = start._array / box_size
             bounding_box["{}max".format(c)] = end._array / box_size
-            select["xyz_{}".format(c)] = select.pop(c)
 
     if new_bbox:
         return _get_cpu_list(bounding_box=bounding_box,

src/osyris/io/loader.py

@@ -70,6 +70,12 @@ def load(self, select=None, cpu_list=None, meta=None):
             for key in _select:
                 if key not in select:
                     _select[key] = False
+        # Replace aliases for x,y,z in select: x,y,x -> position_x,y,z
+        for group in _select.values():
+            if isinstance(group, dict):
+                for c in "xyz":
+                    if c in group:
+                        group[f"position_{c}"] = group.pop(c)
 
         # Take into account user specified lmax
         meta["lmax"] = meta["levelmax"]

src/osyris/io/utils.py

@@ -122,7 +122,7 @@ def make_vector_arrays(data, ndim):
                     cut = ind - 1 if key[ind - 1] == "_" else ind
                     rawkey = key[:cut] + key[ind + 1:]
                     if len(rawkey) == 0:
-                        rawkey = "xyz"
+                        rawkey = "position"
                     data[rawkey] = Array(values=np.array(
                         [data[c].values for c in comp_list]).T,
                                          unit=data[key].unit)