[field] Minor FVM fixes

phi/field/_field_math.py

@@ -252,15 +252,15 @@ def spatial_gradient(field: Field,
     return result
 
 
-def green_gauss_gradient(t: Field, boundary: Extrapolation = None, stack_dim: Shape = channel('vector'), order=2, upwind: Field = None) -> Field:
+def green_gauss_gradient(u: Field, boundary: Extrapolation = None, stack_dim: Shape = channel('vector'), order=2, upwind: Field = None) -> Field:
     """Computes the Green-Gauss gradient of a field at the centroids."""
-    assert stack_dim not in t.shape, f"Gradient dimension is already part of field {t.shape}. Please use a different dimension"
-    boundary = boundary or t.boundary.spatial_gradient()
-    t = t.at_faces(boundary=NONE, order=order, upwind=upwind)
-    normals = rename_dims(t.geometry.face_normals, 'vector', stack_dim)
-    grad = t.geometry.integrate_surface(normals * t.values) / t.geometry.volume
-    grad = slice_off_constant_faces(grad, t.geometry.boundary_elements, boundary)
-    return Field(t.geometry, grad, boundary)
+    assert stack_dim not in u.shape, f"Gradient dimension is already part of field {u.shape}. Please use a different dimension"
+    boundary = boundary or u.boundary.spatial_gradient()
+    u = u.at_faces(boundary=NONE, order=order, upwind=upwind)
+    normals = rename_dims(u.geometry.face_normals, 'vector', stack_dim)
+    grad = u.geometry.integrate_surface(normals * u.values) / u.geometry.volume
+    grad = slice_off_constant_faces(grad, u.geometry.boundary_elements, boundary)
+    return Field(u.geometry, grad, boundary)
 
 
 def _ex_map_f(ext_dict: dict):

phi/physics/advect.py

@@ -87,6 +87,7 @@ def differential(u: Field,
         u: Scalar or vector-valued `Field` sampled on a `CenteredGrid`, `StaggeredGrid` or `UnstructuredMesh`.
         velocity: `Field` that can be sampled at the elements of `u`.
             For FVM, the advection term is typically linearized by setting `velocity = previous_velocity`.
+            Passing `velocity=u` yields non-linear terms which cannot be traced inside linear functions.
         order: Spatial order of accuracy.
             Higher orders entail larger stencils and more computation time but result in more accurate results assuming a large enough resolution.
             Supported for grids: 2 explicit, 4 explicit, 6 implicit (inherited from `phi.field.spatial_gradient()` and resampling).
@@ -115,7 +116,7 @@ def differential(u: Field,
         amounts = velocity * grad
         return - sum(amounts.gradient)
     elif u.is_mesh:
-        u = u.at_faces(boundary=NONE, order=order, upwind=u if upwind is True else upwind)
+        u = u.at_faces(boundary=NONE, order=order, upwind=velocity if upwind is True else upwind)
         velocity = velocity.at_faces(boundary=NONE, order=order, upwind=velocity if upwind is True else upwind)
         conv = density * u.mesh.integrate_surface(u.values * (velocity.values.vector @ velocity.face_normals.vector)) / u.mesh.volume
         return Field(u.geometry, conv, u.boundary)