3.3.0

examples/mesh/FVM_Cylinder_GMsh.ipynb

@@ -13,19 +13,26 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 1,
+   "execution_count": null,
    "metadata": {},
    "outputs": [],
    "source": [
-    "# %pip install --quiet phiflow\n",
+    "# %pip install --quiet phiflow meshio\n",
     "from phi.torch.flow import *\n",
     "# from phi.flow import *  # If JAX is not installed. You can use phi.torch or phi.tf as well.\n",
     "from tqdm.notebook import trange"
    ]
   },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "Let's download the example mesh file and visualize it!"
+   ]
+  },
   {
    "cell_type": "code",
-   "execution_count": 2,
+   "execution_count": null,
    "metadata": {},
    "outputs": [
     {
@@ -59,10 +66,19 @@
     }
    ],
    "source": [
+    "!wget https://raw.githubusercontent.com/tum-pbs/PhiFlow/master/examples/mesh/cylinder.msh -O cylinder.msh\n",
+    "\n",
     "mesh = geom.load_gmsh('cylinder.msh', ('y-', 'x+', 'y+', 'x-', 'cyl+', 'cyl-'))\n",
     "plot(Box(x=6, y=6), mesh, overlay='args', size=(4, 3), title='cylinder.msh')"
    ]
   },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "Next we can define the fluid dynamics. The momentum equation uses diffusion and advection while incompressibility is maintained by solving a separate linear system via `make_incompressible`."
+   ]
+  },
   {
    "cell_type": "code",
    "execution_count": 3,
@@ -82,6 +98,13 @@
     "    return v"
    ]
   },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "Now let's set the boundary and initial conditions and run the simulation!"
+   ]
+  },
   {
    "cell_type": "code",
    "execution_count": 11,

phi/VERSION

@@ -1 +1 @@
-3.2.0
+3.3.0

phi/field/_field.py

@@ -9,6 +9,7 @@
 from phi.math import Shape, Tensor, channel, non_batch, expand, instance, spatial, wrap, dual, non_dual
 from phi.math.extrapolation import Extrapolation
 from phi.math.magic import BoundDim, slicing_dict
+from phiml.dataclasses import sliceable
 from phiml.math import batch, Solve, DimFilter, unstack, concat_shapes, pack_dims, shape
 from phiml.math.extrapolation import domain_slice
 
@@ -45,6 +46,7 @@ def __call__(cls,
         return result
 
 
+@sliceable
 @dataclass(frozen=True)
 class Field(metaclass=_FieldType):
     """
@@ -174,7 +176,7 @@ def numpy(self, order: DimFilter = None):
             assert order is not None, f"order must be specified for non-uniform Field values"
             order = self.values.shape.only(order, reorder=True)
             stack_dims = order.non_uniform_shape
-            inner_order = order.without(stack_dims)
+            inner_order = order.without(stack_dims).names
             return [v.numpy(inner_order) for v in unstack(self.values, stack_dims)]
 
     def uniform_values(self):
@@ -206,7 +208,7 @@ def shape(self) -> Shape:
         if self.is_grid and '~vector' in self.values.shape:
             return batch(self.geometry) & self.resolution & non_dual(self.values).without(self.resolution) & self.geometry.shape['vector']
         set_shape = self.geometry.sets[self.sampled_at]
-        return batch(self.geometry) & (channel(self.geometry) - 'vector') & set_shape & self.values
+        return batch(self.geometry) & (channel(self.geometry) - 'vector') & set_shape & self.values.shape
 
     @property
     def resolution(self):
@@ -232,7 +234,7 @@ def bounds(self) -> BaseBox:
 
         Fields whose spatial rank is determined only during sampling return an empty `Box`.
         """
-        if isinstance(self.geometry.bounds, BaseBox):
+        if hasattr(self.geometry, 'bounds') and isinstance(self.geometry.bounds, BaseBox):
             return self.geometry.bounds
         extent = self.geometry.bounding_half_extent().vector.as_dual('_extent')
         points = self.geometry.center + extent
@@ -242,6 +244,10 @@ def bounds(self) -> BaseBox:
 
     box = bounds
 
+    @property
+    def boundary_names(self):
+        return tuple(self.geometry.boundary_faces)
+
     @property
     def is_grid(self):
         """A Field represents grid data if its `geometry` is a `phi.geom.UniformGrid` instance."""
@@ -662,7 +668,7 @@ def __getitem__(self, item) -> 'Field':
         item = slicing_dict(self, item)
         if not item:
             return self
-        boundary = domain_slice(self.boundary, item, self.resolution)
+        boundary = domain_slice(self.boundary, item, domain_dims=self.boundary_names)
         item_without_vec = {dim: selection for dim, selection in item.items() if dim != 'vector'}
         geometry = self.geometry[item_without_vec]
         if self.is_staggered and 'vector' in item and '~vector' in self.geometry.face_shape:
@@ -682,9 +688,6 @@ def __getitem__(self, item) -> 'Field':
         values = self.values[item]
         return Field(geometry, values, boundary)
 
-    def __getattr__(self, name: str) -> BoundDim:
-        return BoundDim(self, name)
-
     def dimension(self, name: str):
         """
         Returns a reference to one of the dimensions of this field.

phi/field/_field_io.py

@@ -3,6 +3,7 @@
 import numpy as np
 
 from phi import geom, math
+from phiml.math._shape import from_dict
 from ._field import Field
 from ._grid import unstack_staggered_tensor, CenteredGrid, StaggeredGrid
 from ._field_math import stack
@@ -55,7 +56,7 @@ def write_single_field(field: Field, file: str):
         field_type = 'StaggeredGrid' if field.is_staggered else 'CenteredGrid'
         np.savez_compressed(file,
                             dim_names=dim_names,
-                            dim_types=field.shape.types,
+                            dim_types=field.shape.dim_types,
                             dim_item_names=np.asarray(field.shape.item_names, dtype=object),
                             field_type=field_type,
                             lower=lower,
@@ -103,7 +104,8 @@ def read_single_field(file: str, convert_to_backend=True) -> Field:
         raise NotImplementedError(f"{ftype} not implemented")
     data_arr = stored['data']
     dim_item_names = stored.get('dim_item_names', (None,) * len(data_arr.shape))
-    data = tensor(data_arr, Shape(data_arr.shape, tuple(stored['dim_names']), tuple(stored['dim_types']), tuple(dim_item_names)), convert=convert_to_backend)
+    shape_spec = {'names': tuple(stored['dim_names']), 'sizes': data_arr.shape, 'types': tuple(stored['dim_types']), 'item_names': tuple(dim_item_names)}
+    data = tensor(data_arr, from_dict(shape_spec), convert=convert_to_backend)
     bounds_item_names = stored.get('bounds_item_names', None)
     if bounds_item_names is None or bounds_item_names.shape == ():  # None or empty array
         bounds_item_names = spatial(data).names

phi/field/_field_math.py

@@ -2,8 +2,9 @@
 from typing import Callable, List, Tuple, Optional, Union, Sequence
 
 import numpy as np
-from phiml.math import Tensor, spatial, instance, tensor, channel, batch, Shape, unstack, solve_linear, jit_compile_linear, \
-    shape, Solve, extrapolation, dual, wrap, rename_dims, factorial, concat, zeros, ones
+from phiml.math import Tensor, spatial, instance, tensor, channel, batch, Shape, unstack, solve_linear, \
+    jit_compile_linear, \
+    shape, Solve, extrapolation, dual, wrap, rename_dims, factorial, concat, zeros, ones, neighbor_mean
 from phi import geom
 from phi import math
 from phi.geom import Box, Geometry, UniformGrid
@@ -50,7 +51,8 @@ def laplace(u: Field,
             implicitness: int = None,
             weights: Union[Tensor, Field] = None,
             upwind: Field = None,
-            correct_skew=True) -> Field:
+            correct_skew=True,
+            wide_stencil=False) -> Field:
     """
     Spatial Laplace operator for scalar grid.
 
@@ -77,7 +79,6 @@ def laplace(u: Field,
     Returns:
         laplacian field as `CenteredGrid`
     """
-
     if implicitness is None:
         implicitness = 0 if implicit is None else 2
     elif implicitness != 0:
@@ -94,6 +95,11 @@ def laplace(u: Field,
             raise NotImplementedError(f"laplace on meshes is not yet supported with vector-valued weights")
         neighbor_val = u.mesh.pad_boundary(u.values, mode=u.boundary)
         nb_distances = u.mesh.neighbor_distances
+        if wide_stencil:
+            assert weights is None
+            grad_p = spatial_gradient(u, order=order, scheme='green-gauss', upwind=upwind)
+            div_grad_p = grad_p.divergence(order=order, upwind=upwind)
+            return div_grad_p
         connecting_grad = (u.mesh.connectivity * neighbor_val - u.values) / nb_distances  # (T_N - T_P) / d_PN
         if correct_skew and gradient is not None:  # skewness correction
             assert dual(gradient).names == ('~vector',), f"gradient must contain one dual dim '~vector' listing the gradient components but got {gradient.shape}"
@@ -109,41 +115,33 @@ def laplace(u: Field,
         laplace_values = u.mesh.integrate_surface(grad) / u.mesh.volume  # 1/V ∑_f ∇T ν A
         result = weights * laplace_values if weights is not None else laplace_values
         return Field(u.mesh, result, u.boundary - u.boundary)
-
     # --- Grid ---
-
     laplace_ext = u.extrapolation.spatial_gradient().spatial_gradient()
     laplace_dims = u.shape.only(axes).names
-
-    if u.vector.exists and (u.is_centered or order > 2):
+    if 'vector' in u.shape and (u.is_centered or order > 2):
         fields = [f for f in u.vector]
     else:
         fields = [u]
-
     result = []
     for f in fields:
         if order == 2:
             result.append(math.map_d2c(math.laplace)(f.values, dx=f.dx, padding=f.extrapolation, dims=axes, weights=weights, padding_kwargs={'bounds': f.bounds}))  # uses ghost cells
         else:
-            result_components = [perform_finite_difference_operation(f.values, dim, 2, f.dx.vector[dim], f.extrapolation,
-                                                                         laplace_ext, 'center', order, implicit,
-                                                                         implicitness) for dim in laplace_dims]
+            result_components = [perform_finite_difference_operation(f.values, dim, 2, f.dx.vector[dim], f.extrapolation, laplace_ext, 'center', order, implicit, implicitness) for dim in laplace_dims]
             if weights is not None:
                 if channel(weights):
                     result_components = [c * weights[ax] for c, ax in zip(result_components, axes_names)]
                 else:
                     result_components = [c * weights for c in result_components]
 
             result.append(sum(result_components))
-
-    if u.vector.exists and (u.is_centered or order > 2):
+    if 'vector' in u.shape and (u.is_centered or order > 2):
         if u.is_staggered:
             result = math.stack(result, dual(vector=u.vector.item_names))
         else:
             result = math.stack(result, channel(vector=u.vector.item_names))
     else:
         result = result[0]
-
     return u.with_values(result).with_extrapolation(laplace_ext)
 
 
@@ -205,17 +203,14 @@ def spatial_gradient(field: Field,
             return least_squares_gradient(field, stack_dim=stack_dim, boundary=boundary)
         raise NotImplementedError(scheme)
 
-
-    if field.vector.exists:
+    if 'vector' in field.shape:
         assert stack_dim.name != 'vector', "`stack_dim=vector` is inadmissible if the input is a vector grid"
         if field == StaggeredGrid:
             assert at == 'faces', "for a `StaggeredGrid` input only `type == StaggeredGrid` is possible"
 
     if at == 'faces':
         assert stack_dim.name == 'vector', f"spatial_gradient with type=StaggeredGrid requires stack_dim.name == 'vector' but got '{stack_dim.name}'"
 
-
-
     if gradient_extrapolation is None:
         gradient_extrapolation = field.extrapolation.spatial_gradient()
 
@@ -657,8 +652,8 @@ def curl(field: Field, at='corner'):
     if field.is_grid and field.is_staggered and field.spatial_rank == 2 and at == 'corner':
         x, y = field.vector.item_names
         values = field.with_boundary(None).values
-        vx = math.pad(values.vector.dual[x], {y: (1, 1)}, field.extrapolation[{'vector': y}])
-        vy = math.pad(values.vector.dual[y], {x: (1, 1)}, field.extrapolation[{'vector': x}])
+        vx = math.pad(values.vector.dual[x], {y: (1, 1)}, field.boundary[{'vector': y}])
+        vy = math.pad(values.vector.dual[y], {x: (1, 1)}, field.boundary[{'vector': x}])
         vy_dx = math.spatial_gradient(vy, dims=x, dx=field.dx[x], padding=None, stack_dim=None, difference='forward')
         vx_dy = math.spatial_gradient(vx, dims=y, dx=field.dx[y], padding=None, stack_dim=None, difference='forward')
         curl_val = vy_dx - vx_dy
@@ -699,6 +694,36 @@ def curl(field: Field, at='corner'):
     #         vy_dx = math.spatial_gradient(y_padded, field.dx, 'forward', None, dims='x', stack_dim=None)
     #         result = vy_dx - vx_dy
     #         return CenteredGrid(result, field.extrapolation.spatial_gradient(), field.bounds)
+    elif field.is_grid and field.is_staggered and field.spatial_rank == 3 and at == 'corner':
+        x, y, z = field.vector.item_names
+        values = field.with_boundary(None).values
+        vx = math.pad(values.vector.dual[x], {y: (1, 1), z: (1, 1)}, field.boundary[{'vector': y}])
+        vy = math.pad(values.vector.dual[y], {x: (1, 1), z: (1, 1)}, field.boundary[{'vector': x}])
+        vz = math.pad(values.vector.dual[z], {x: (1, 1), y: (1, 1)}, field.boundary[{'vector': z}])
+        vx_dy = neighbor_mean(math.spatial_gradient(vx, dims=y, dx=field.dx[y], padding=None, stack_dim=None, difference='forward'), z)
+        vx_dz = neighbor_mean(math.spatial_gradient(vx, dims=z, dx=field.dx[z], padding=None, stack_dim=None, difference='forward'), y)
+        vy_dx = neighbor_mean(math.spatial_gradient(vy, dims=x, dx=field.dx[x], padding=None, stack_dim=None, difference='forward'), z)
+        vy_dz = neighbor_mean(math.spatial_gradient(vy, dims=z, dx=field.dx[z], padding=None, stack_dim=None, difference='forward'), x)
+        vz_dx = neighbor_mean(math.spatial_gradient(vz, dims=x, dx=field.dx[x], padding=None, stack_dim=None, difference='forward'), y)
+        vz_dy = neighbor_mean(math.spatial_gradient(vz, dims=y, dx=field.dx[y], padding=None, stack_dim=None, difference='forward'), x)
+        curl_val = math.stack([vz_dy-vy_dz, vx_dz-vz_dx, vy_dx-vx_dy], field.shape['vector'])
+        corners = UniformGrid(field.resolution + 1, Box(field.bounds.lower - field.dx / 2, field.bounds.upper + field.dx / 2))
+        return Field(corners, curl_val, field.boundary.spatial_gradient())
+    elif field.is_grid and field.is_centered and field.spatial_rank == 3 and at == 'corner':
+        raise NotImplementedError
+        x, y, z = field.vector.item_names
+        values = pad(field, 1).values
+        # ToDo 8 diag offset vectors, account for cell stretching
+        # Then sum (offset x v) / |offset|^2 ??
+        diag_basis = wrap([(1, 1), (1, -1)], channel(diag='pos,neg'), dual(vector=[x, y]))
+        diag_comp = diag_basis @ values
+        ll = diag_comp[{x: slice(-1), y: slice(-1), 'diag': 'neg'}]
+        ul = diag_comp[{x: slice(-1), y: slice(1, None), 'diag': 'pos'}]
+        lr = diag_comp[{x: slice(1, None), y: slice(-1), 'diag': 'pos'}]
+        ur = diag_comp[{x: slice(1, None), y: slice(1, None), 'diag': 'neg'}]
+        curl_val = ll - ul + lr - ur
+        corners = UniformGrid(field.resolution + 1, Box(field.bounds.lower - field.dx / 2, field.bounds.upper + field.dx / 2))
+        return Field(corners, curl_val, field.boundary.spatial_gradient())
     raise NotImplementedError("Only 2D curl at corner currently supported")
 
 
@@ -762,7 +787,10 @@ def mean(field: Field, dim=lambda s: s.non_channel.non_batch) -> Tensor:
     Returns:
         `phi.Tensor`
     """
-    result = math.mean(field.values, dim=dim)
+    if field.is_grid:
+        result = math.mean(field.values, dim=dim)
+    else:
+        result = math.mean(field.values, dim=dim, weight=field.geometry.volume)
     if (instance(field.geometry) & spatial(field.geometry)) in result.shape:
         return field.with_values(result)
     return result
@@ -948,7 +976,10 @@ def stack(fields: Sequence[Field], dim: Shape, dim_bounds: Box = None):
             return fields[0].with_values(values).with_boundary(boundary)
     else:
         values = math.stack([f.values for f in fields], dim)
-        geometry = fields[0].geometry if all(f.geometry == fields[0].geometry for f in fields) else math.stack([f.geometry for f in fields], dim)
+        geometry = fields[0].geometry if all(f.geometry == fields[0].geometry for f in fields) else math.stack([f.geometry for f in fields], dim, layout_non_matching=True)
+        if isinstance(geometry, Tensor):
+            from phi.geom._geom_ops import GeometryStack
+            geometry = GeometryStack(geometry)
         return Field(geometry, values, boundary)
 
 

phi/field/_grid.py

@@ -158,7 +158,7 @@ def StaggeredGrid(values: Any = 0.,
             values = sample(values, elements, at='face', boundary=extrapolation, dot_face_normal=elements)
         elif callable(values):
             values = sample_function(values, elements, 'face', extrapolation)
-            if elements.shape.shape.rank > 1:  # Different number of X and Y faces
+            if elements.shape.is_non_uniform:  # Different number of X and Y faces
                 assert isinstance(values, TensorStack), f"values function must return a staggered Tensor but returned {type(values)}"
             assert '~vector' in values.shape
             if 'vector' in values.shape:
@@ -203,7 +203,7 @@ def resolution_from_staggered_tensor(values: Tensor, extrapolation: Extrapolatio
     x_shape = values.shape.after_gather({'vector': any_dim, '~vector': any_dim})
     ext_lower, ext_upper = extrapolation.valid_outer_faces(any_dim)
     delta = int(ext_lower) + int(ext_upper) - 1
-    resolution = x_shape.spatial._replace_single_size(any_dim, x_shape.get_size(any_dim) - delta)
+    resolution = x_shape.spatial.with_dim_size(any_dim, x_shape.get_size(any_dim) - delta)
     return resolution
 
 

phi/field/_resample.py

@@ -7,7 +7,7 @@
 from phi.math.extrapolation import Extrapolation, ConstantExtrapolation, PERIODIC
 from phiml.math import unstack, channel, rename_dims, batch, extrapolation
 from ._field import Field, FieldInitializer, as_boundary, slice_off_constant_faces
-from phiml.math._tensors import may_vary_along
+from phiml.math._tensors import may_vary_along, wrap
 
 
 def resample(value: Union[Field, Geometry, Tensor, float, FieldInitializer], to: Union[Field, Geometry], keep_boundary=False, **kwargs):
@@ -48,7 +48,9 @@ def resample(value: Union[Field, Geometry, Tensor, float, FieldInitializer], to:
         >>> field.resample(grid, to=grid) == grid
         True
     """
-    assert isinstance(to, (Field, Geometry)), f"'to' must be a Field or Geomoetry but got {to}"
+    assert isinstance(to, (Field, Geometry)), f"'to' must be a Field or Geometry but got {to}"
+    if isinstance(to, Geometry):
+        to = Field(to, wrap(0.), value.extrapolation if isinstance(value, Field) else 0.)
     if not isinstance(value, (Field, Geometry, FieldInitializer)):
         return to.with_values(value)
     if isinstance(value, Field) and keep_boundary: