FVM WIP 5

demos/fvm.py

@@ -8,9 +8,10 @@
 # v = v.at_faces(scheme='linear')
 # show(v)
 div = field.divergence(v)
-# grad = field.spatial_gradient(div)
-conv = convection(v, v)
-show(div, conv)
+grad = field.spatial_gradient(div, scheme='green-gauss')
+# conv = convection(v, v)
+# diff = diffusion(v, None, )
+show(div, grad)
 
 # field.sample(v, v.elements, 'face', v.extrapolation, scheme='linear')
 # v = v.with_values(math.random_uniform(non_channel(v.values)))

phi/field/_field_math.py

@@ -1,7 +1,7 @@
 from numbers import Number
 from typing import Callable, List, Tuple, Optional, Union, Sequence
 
-from phiml.math import Tensor, spatial, instance, tensor, channel, Shape, unstack, solve_linear, jit_compile_linear, shape, Solve, extrapolation, dual, wrap
+from phiml.math import Tensor, spatial, instance, tensor, channel, Shape, unstack, solve_linear, jit_compile_linear, shape, Solve, extrapolation, dual, wrap, rename_dims
 from phi import geom
 from phi import math
 from phi.geom import Box, Geometry, UniformGrid
@@ -103,12 +103,15 @@ def laplace(field: Field,
 
 
 def spatial_gradient(field: Field,
-                     gradient_extrapolation: Extrapolation = None,
+                     boundary: Extrapolation = None,
                      at: str = 'center',
                      dims: math.DimFilter = spatial,
                      stack_dim: Shape = channel('vector'),
                      order=2,
-                     implicit: Solve = None):
+                     implicit: Solve = None,
+                     scheme = None,
+                     interpolation = 'linear',
+                     gradient_extrapolation: Extrapolation = None):
     """
     Finite difference spatial_gradient.
 
@@ -119,7 +122,7 @@ def spatial_gradient(field: Field,
 
     Args:
         field: centered grid of any number of dimensions (scalar field, vector field, tensor field)
-        gradient_extrapolation: Extrapolation of the output
+        boundary: Boundary conditions of the gradient field.
         at: Either `'face'` or `'center'`
         dims: Along which dimensions to compute the spatial gradient. Only supported when `type==CenteredGrid`.
         stack_dim: Dimension to be added. This dimension lists the spatial_gradient w.r.t. the spatial dimensions.
@@ -129,15 +132,21 @@ def spatial_gradient(field: Field,
             Supported: 2 explicit, 4 explicit, 6 implicit.
         implicit: When a `Solve` object is passed, performs an implicit operation with the specified solver and tolerances.
             Otherwise, an explicit stencil is used.
+        gradient_extrapolation: Alias for `boundary`.
 
     Returns:
         spatial_gradient field of type `type`.
     """
     assert at in ['face', 'center']
-    if gradient_extrapolation is None:
-        gradient_extrapolation = field.extrapolation.spatial_gradient()
+    if gradient_extrapolation is not None:
+        assert boundary is None, f"Cannot specify both boundary and gradient_extrapolation"
+        boundary = gradient_extrapolation
+    if boundary is None:
+        boundary = field.extrapolation.spatial_gradient()
     if field.is_mesh and at == 'center':
-        raise NotImplementedError
+        if scheme == 'green-gauss':
+            return green_gauss_gradient(field, interpolation=interpolation, stack_dim=stack_dim, boundary=boundary)
+        raise NotImplementedError(scheme)
     extrap_map = {}
     if not implicit:
         if order == 2:
@@ -168,39 +177,49 @@ def spatial_gradient(field: Field,
     base_widths = (abs(min(needed_shifts)), max(needed_shifts))
     field.with_extrapolation(extrapolation.map(_ex_map_f(extrap_map), field.extrapolation))  # ToDo does this line do anything?
     if implicit:
-        gradient_extrapolation = extrapolation.map(_ex_map_f(extrap_map_rhs), gradient_extrapolation)
+        boundary = extrapolation.map(_ex_map_f(extrap_map_rhs), boundary)
     spatial_dims = field.shape.only(dims).names
     stack_dim = stack_dim.with_size(spatial_dims)
     if at == 'center':
         # ToDo if extrapolation == math.extrapolation.NONE, extend size by 1
         # pad = 1 if extrapolation == math.extrapolation.NONE else 0
         # bounds = Box(field.bounds.lower - field.dx, field.bounds.upper + field.dx) if extrapolation == math.extrapolation.NONE else field.bounds
         std_widths = (0, 0)
-        if gradient_extrapolation == math.extrapolation.NONE:
+        if boundary == math.extrapolation.NONE:
             base_widths = (abs(min(needed_shifts))+1, max(needed_shifts)+1)
             std_widths = (1, 1)
         padded_components = [math.pad(field.values, {dim_: base_widths if dim_ == dim else std_widths for dim_ in spatial_dims}, field.extrapolation) for dim in spatial_dims]
         shifted_components = [math.shift(padded_component, needed_shifts, stack_dim=None, padding=None, dims=dim) for padded_component, dim in zip(padded_components, spatial_dims)]
         result_components = [sum([value * shift_ for value, shift_ in zip(values, shifted_component)]) / field.dx.vector[dim] for shifted_component, dim in zip(shifted_components, field.shape.spatial.names)]
-        result = CenteredGrid(math.stack(result_components, stack_dim), gradient_extrapolation, field.bounds)
+        result = CenteredGrid(math.stack(result_components, stack_dim), boundary, field.bounds)
     elif at == 'face':
         assert spatial_dims == field.shape.spatial.names, f"spatial_gradient with type=StaggeredGrid requires dims=spatial, i.e. dims='{','.join(field.shape.spatial.names)}'"
         base_widths = (base_widths[0], base_widths[1]-1)
         padded_components = []
         for dim in spatial_dims:
-            lo, up = gradient_extrapolation.valid_outer_faces(dim)
+            lo, up = boundary.valid_outer_faces(dim)
             padding_widths = (lo + base_widths[0], up + base_widths[1])
             padded_components.append(math.pad(field.values, {dim: padding_widths}, field.extrapolation, bounds=field.bounds))
         shifted_components = [math.shift(padded_component, needed_shifts, stack_dim=None, padding=None, dims=dim) for padded_component, dim in zip(padded_components, spatial_dims)]
         result_components = [sum([value * shift_ for value, shift_ in zip(values, shifted_component)]) / field.dx.vector[dim] for shifted_component, dim in zip(shifted_components, field.shape.spatial.names)]
         assert stack_dim.name == 'vector', f"spatial_gradient with type=StaggeredGrid requires stack_dim.name == 'vector' but got '{stack_dim.name}'"
-        result = StaggeredGrid(math.stack(result_components, dual(vector=spatial_dims)), gradient_extrapolation, field.bounds)
+        result = StaggeredGrid(math.stack(result_components, dual(vector=spatial_dims)), boundary, field.bounds)
     if implicit:
         implicit.x0 = result
         result = solve_linear(_lhs_for_implicit_scheme, result, solve=implicit, values_rhs=values_rhs, needed_shifts_rhs=needed_shifts_rhs, stack_dim=stack_dim, staggered_output=type != CenteredGrid)
     return result
 
 
+def green_gauss_gradient(t: Field, interpolation='linear', stack_dim: Shape = channel('vector'), boundary: Extrapolation = None) -> Field:
+    """Computes the Green-Gauss gradient of a field at the centroids."""
+    boundary = boundary or t.boundary.spatial_gradient()
+    t = t.at_faces(interpolation=interpolation)
+    normals = rename_dims(t.face_normals, 'vector', stack_dim)
+    grad = integrate_surface(t.geometry, normals * t.values, t.boundary) / t.geometry.volume
+    grad = slice_off_constant_faces(grad, t.geometry.boundary_elements, boundary)
+    return Field(t.geometry, grad, boundary)
+
+
 def _ex_map_f(ext_dict: dict):
     def f(ext: Extrapolation):
         return ext_dict[ext.__repr__()] if ext.__repr__() in ext_dict else ext
@@ -296,7 +315,7 @@ def stagger(field: Field,
         return CenteredGrid(values, bounds=field.bounds, extrapolation=extrapolation)
 
 
-def divergence(field: Field, order=2, implicit: Solve = None, scheme: str = 'linear', upwind_vectors: Tensor = None) -> CenteredGrid:
+def divergence(field: Field, order=2, implicit: Solve = None, interpolation: str = 'linear', upwind_vectors: Tensor = None) -> CenteredGrid:
     """
     Computes the divergence of a grid using finite differences.
 
@@ -317,7 +336,7 @@ def divergence(field: Field, order=2, implicit: Solve = None, scheme: str = 'lin
         Divergence field as `CenteredGrid`
     """
     if field.is_mesh:
-        field = field.at_faces(scheme=scheme, upwind_vectors=upwind_vectors)
+        field = field.at_faces(interpolation=interpolation, upwind_vectors=upwind_vectors)
         div = integrate_flux(field.geometry, field.values, field.boundary, divide_volume=True)
         return Field(field.geometry, div, field.boundary.spatial_gradient())
     extrap_map = {}

phi/field/_resample.py

@@ -350,29 +350,30 @@ def _shift_resample(self: Field, resolution: Shape, bounds: Box, threshold=1e-5,
     return data
 
 
-def centroid_to_faces(field: Field, boundary: Extrapolation, scheme='upwind-linear', upwind_vectors: Field = None, gradient: Field = None):
+def centroid_to_faces(field: Field, boundary: Extrapolation, interpolation='upwind-linear', upwind_vectors: Field = None, gradient: Field = None):
     mesh = field.geometry
     val = field.values
     if '~neighbors' in val.shape:
         return val
     neighbor_val = mesh.pad_boundary(val, mode=field.boundary)
-    if scheme == 'upwind-linear':
+    if interpolation == 'upwind-linear':
         flows_out = (upwind_vectors or val).vector @ mesh.face_normals.vector >= 0
         if gradient is None:
-            gradient = green_gauss_gradient(field)
+            from phi.field._field_math import green_gauss_gradient
+            gradient = green_gauss_gradient(field, interpolation=interpolation)
         neighbor_grad = mesh.pad_boundary(gradient.values, mode=gradient.boundary)
         interpolated_from_self = val + gradient.values.vector @ (mesh.face_centers - mesh.center).vector
         interpolated_from_neighbor = neighbor_val + neighbor_grad.vector @ (mesh.face_centers - (mesh.center + mesh.neighbor_offsets)).vector
         # ToDo limiter
         result = math.where(flows_out, interpolated_from_self, interpolated_from_neighbor)
-        return slice_off_constant_faces(result, mesh.boundary_faces, mesh.boundary)
-    elif scheme == 'upwind':
+        return slice_off_constant_faces(result, mesh.boundary_faces, boundary)
+    elif interpolation == 'upwind':
         flows_out = (upwind_vectors or val).vector @ mesh.face_normals.vector >= 0
         return math.where(flows_out, val, neighbor_val)
-    elif scheme == 'skewfree-linear':  # does not take skewness into account
+    elif interpolation == 'skewfree-linear':  # does not take skewness into account
         relative_face_distance = slice_off_constant_faces(mesh.relative_face_distance, mesh.boundary_faces, boundary)
         return (1 - relative_face_distance) * field.values + relative_face_distance * neighbor_val
-    elif scheme == 'linear':
+    elif interpolation == 'linear':
         nb_center = math.replace_dims(mesh.center, 'cells', math.dual('~neighbors'))
         cell_deltas = math.pairwise_distances(mesh.center, format=mesh.cell_connectivity, default=None)  # x_N - x_P
         face_distance = nb_center - mesh.face_centers[mesh.interior_faces]  # x_N - x_f
@@ -385,14 +386,7 @@ def centroid_to_faces(field: Field, boundary: Extrapolation, scheme='upwind-line
         # b0 = math.tensor_like(slice_off_constant_faces(mesh.connectivity, mesh.boundary_faces, boundary), 0)
         return w * val + (1 - w) * neighbor_val
     else:
-        raise NotImplementedError(f"Scheme '{scheme}' not supported for resampling mesh values to faces")
-
-
-def green_gauss_gradient(t: Field, scheme='linear', stack_dim: Shape = channel('vector')) -> Field:
-    t = t.at_faces(scheme=scheme)
-    normals = rename_dims(t.face_normals, 'vector', stack_dim)
-    grad = integrate_surface(t.mesh, normals * t.values, t.boundary) / t.mesh.volume
-    return Field(t.geometry, grad, t.boundary.spatial_gradient())
+        raise NotImplementedError(f"Interpolation scheme '{interpolation}' not supported for resampling mesh values to faces")
 
 
 def sample_function(f: Callable, elements: Geometry, at: str, extrapolation: Extrapolation) -> Tensor:

phi/geom/_mesh.py

@@ -92,7 +92,7 @@ def face_shape(self) -> Shape:
         return self.face_areas.shape
 
     @property
-    def boundary_elements(self) -> Dict[str, Tuple[Dict[str, slice], Dict[str, slice]]]:
+    def boundary_elements(self) -> Dict[str, Dict[str, slice]]:
         return {}
 
     @property
@@ -268,9 +268,7 @@ def load_su2(file_or_mesh: str, cell_dim=instance('cells'), face_format: str = '
         face_format: Sparse storage format for cell connectivity.
 
     Returns:
-        surface: `UnstructuredMesh`
-        markers: Edges/Faces marked
-            sparse (vertices, vertices) -> int
+        `UnstructuredMesh`
     """
     if isinstance(file_or_mesh, str):
         from ezmesh import import_from_file
@@ -288,8 +286,24 @@ def load_su2(file_or_mesh: str, cell_dim=instance('cells'), face_format: str = '
 def mesh_from_numpy(points: Union[list, np.ndarray],
                     polygons: list,
                     boundaries: Dict[str, List[Sequence]],
-                    cell_dim: Shape,
-                    face_format: str) -> UnstructuredMesh:
+                    cell_dim: Shape = instance('cells'),
+                    face_format: str = 'csc') -> UnstructuredMesh:
+    """
+    Construct an unstructured mesh from vertices.
+
+    Args:
+        points: 2D numpy array of shape (num_points, point_coord).
+            The last dimension must have length 2 for 2D meshes and 3 for 3D meshes.
+        polygons: List of polygons. Each polygon is defined as a sequence of point indices mapping into `points'.
+            E.g. `[(0, 1, 2)]` denotes a single triangle connecting points 0, 1, and 2.
+        boundaries: An unstructured mesh can have multiple boundaries, each defined by a name `str` and a list of faces, defined by their vertices.
+            The `boundaries` `dict` maps boundary names to a list of edges (point pairs) in 2D and faces (3 or more points) in 3D (not yet supported).
+        cell_dim: Dimension along which to list the cells. This should be an instance dimension.
+        face_format: Sparse storage format for cell connectivity.
+
+    Returns:
+        `UnstructuredMesh`
+    """
     cell_dim = cell_dim.with_size(len(polygons))
     points = np.asarray(points)
     dim = points.shape[-1]

phi/physics/fvm.py

@@ -1,3 +1,5 @@
+from typing import Optional
+
 from phiml.math import math, Extrapolation, Tensor
 from ..field import Field
 from ..geom import UnstructuredMesh
@@ -56,7 +58,19 @@
 #     return dynamic_viscosity * mesh.integrate_surface(grad) / mesh.volume  # 1/V ∑_f ∇T ν A
 
 
-def diffusion(v: Field, prev_grad: Tensor, dynamic_viscosity=1., scheme='linear') -> Tensor:
+def diffusion(v: Field, prev_grad: Optional[Field], dynamic_viscosity=1., scheme='linear') -> Field:
+    """
+    Computes the FVM diffusion term.
+
+    Args:
+        v: Scalar or vector field sampled at the centroids of an unstructured mesh.
+        prev_grad:
+        dynamic_viscosity:
+        scheme:
+
+    Returns:
+
+    """
     neighbor_val = v.mesh.pad_boundary(v.values, mode=v.boundary)
     connecting_grad = (v.mesh.connectivity * neighbor_val - v.values) / v.mesh.neighbor_distances  # (T_N - T_P) / d_PN
     if prev_grad is not None:  # skewness correction
@@ -86,4 +100,5 @@ def convection(v: Field, prev_v: Field, density: float = 1, scheme='upwind-linea
     """
     v = v.at_faces(scheme=scheme)
     prev_v = prev_v.at_faces(scheme=scheme)
-    return density * integrate_surface(v.mesh, v * (prev_v.values.vector @ prev_v.face_normals.vector), prev_v.boundary) / v.mesh.volume
+    conv = density * integrate_surface(v.mesh, v.values * (prev_v.values.vector @ prev_v.face_normals.vector), prev_v.boundary) / v.mesh.volume
+    return Field(v.geometry, conv, v.boundary)