Choose the right finite-difference stencil for D3Q27

pystencils/fd/derivation.py

@@ -228,9 +228,10 @@ class FiniteDifferenceStaggeredStencilDerivation:
         neighbor: the neighbor direction string or vector at whose staggered position to calculate the derivative
         dim: how many dimensions (2 or 3)
         derivative: a tuple of directions over which to perform derivatives
+        free_weights_prefix: a string to prefix to free weight symbols. If None, do not return free weights
     """
 
-    def __init__(self, neighbor, dim, derivative=tuple()):
+    def __init__(self, neighbor, dim, derivative=tuple(), free_weights_prefix=None):
         if type(neighbor) is str:
             neighbor = direction_string_to_offset(neighbor)
         if dim == 2:
@@ -281,7 +282,10 @@ class FiniteDifferenceStaggeredStencilDerivation:
 
             # if the weights are underdefined, we can choose the free symbols to find the sparsest stencil
             free_weights = set(itertools.chain(*[w.free_symbols for w in weights]))
-            if len(free_weights) > 0:
+            if free_weights_prefix is not None:
+                weights = [w.subs({fw: sp.Symbol(f"{free_weights_prefix}_{i}") for i, fw in enumerate(free_weights)})
+                           for w in weights]
+            elif len(free_weights) > 0:
                 zero_counts = defaultdict(list)
                 for values in itertools.product([-1, -sp.Rational(1, 2), 0, 1, sp.Rational(1, 2)],
                                                 repeat=len(free_weights)):

pystencils/fd/finitevolumes.py

@@ -59,7 +59,10 @@ class FVM1stOrder:
 
         assert ps.FieldType.is_staggered(flux_field)
 
+        num = 0
+
         def discretize(term, neighbor):
+            nonlocal num
             if isinstance(term, sp.Matrix):
                 nw = term.applyfunc(lambda t: discretize(t, neighbor))
                 return nw
@@ -69,7 +72,9 @@ class FVM1stOrder:
             elif isinstance(term, ps.fd.Diff):
                 access, direction = get_access_and_direction(term)
 
-                fds = FDS(neighbor, access.field.spatial_dimensions, direction)
+                fds = FDS(neighbor, access.field.spatial_dimensions, direction,
+                          free_weights_prefix=f'fvm_free_{num}' if sp.Matrix(neighbor).dot(neighbor) > 2 else None)
+                num += 1
                 return fds.apply(access)
 
             if term.args:
@@ -91,7 +96,20 @@ class FVM1stOrder:
             directional_flux = fluxes[0] * int(neighbor[0])
             for i in range(1, self.dim):
                 directional_flux += fluxes[i] * int(neighbor[i])
-            discrete_flux = discretize(directional_flux, neighbor)
+            discrete_flux = sp.simplify(discretize(directional_flux, neighbor))
+            free_weights = [s for s in discrete_flux.atoms(sp.Symbol) if s.name.startswith('fvm_free_')]
+
+            if len(free_weights) > 0:
+                discrete_flux = discrete_flux.collect(discrete_flux.atoms(ps.field.Field.Access))
+                access_counts = defaultdict(list)
+                for values in itertools.product([-1, 0, 1],
+                                                repeat=len(free_weights)):
+                    subs = {free_weight: value for free_weight, value in zip(free_weights, values)}
+                    simp = discrete_flux.subs(subs)
+                    access_count = len(simp.atoms(ps.field.Field.Access))
+                    access_counts[access_count].append(simp)
+                best_count = min(access_counts.keys())
+                discrete_flux = sum(access_counts[best_count]) / len(access_counts[best_count])
             discrete_fluxes.append(discrete_flux / sp.Matrix(neighbor).norm())
 
         if flux_field.index_dimensions > 1: