From 3a0e87f7cc8bb4a39aa119d8488fcf3ce01d2f9d Mon Sep 17 00:00:00 2001
From: Helen Schottenhamml <helen.schottenhamml@fau.de>
Date: Thu, 2 May 2024 17:37:38 +0200
Subject: [PATCH] Wall function bounce
---
doc/sphinx/lbmpy.bib | 26 ++
src/lbmpy/boundaries/__init__.py | 8 +-
src/lbmpy/boundaries/boundaryconditions.py | 237 +++++++++++++++++++
src/lbmpy/boundaries/wall_function_models.py | 173 ++++++++++++++
src/lbmpy/flow_statistics.py | 82 +++++++
tests/test_wall_function_bounce.py | 116 +++++++++
6 files changed, 639 insertions(+), 3 deletions(-)
create mode 100644 src/lbmpy/boundaries/wall_function_models.py
create mode 100644 src/lbmpy/flow_statistics.py
create mode 100644 tests/test_wall_function_bounce.py
diff --git a/doc/sphinx/lbmpy.bib b/doc/sphinx/lbmpy.bib
index efce2c1a..c9d012c1 100644
--- a/doc/sphinx/lbmpy.bib
+++ b/doc/sphinx/lbmpy.bib
@@ -288,4 +288,30 @@ journal = {Communications in Computational Physics}
journal = {Physics of Fluids}
}
+@article{Han2021,
+ doi = {10.1088/1873-7005/ac1782},
+ url = {https://dx.doi.org/10.1088/1873-7005/ac1782} ,
+ year = {2021},
+ month = {aug},
+ publisher = {IOP Publishing},
+ volume = {53},
+ number = {4},
+ pages = {045506},
+ author = {Mengtao Han and Ryozo Ooka and Hideki Kikumoto},
+ title = {A wall function approach in lattice Boltzmann method: algorithm and validation using turbulent channel flow},
+ journal = {Fluid Dynamics Research}
+}
+
+@article{Maronga2020,
+ author = {Maronga, Bj{\"o}rn and Knigge, Christoph and Raasch, Siegfried},
+ year = {2020},
+ title = {{An Improved Surface Boundary Condition for Large-Eddy Simulations Based on Monin--Obukhov Similarity Theory: Evaluation and Consequences for Grid Convergence in Neutral and Stable Conditions}},
+ pages = {297--325},
+ volume = {174},
+ number = {2},
+ issn = {0006-8314},
+ journal = {{Boundary-layer meteorology}},
+ doi = {10.1007/s10546-019-00485-w}
+}
+
@Comment{jabref-meta: databaseType:bibtex;}
diff --git a/src/lbmpy/boundaries/__init__.py b/src/lbmpy/boundaries/__init__.py
index 5520f318..e522740a 100644
--- a/src/lbmpy/boundaries/__init__.py
+++ b/src/lbmpy/boundaries/__init__.py
@@ -1,10 +1,12 @@
from lbmpy.boundaries.boundaryconditions import (
- UBB, FixedDensity, DiffusionDirichlet, SimpleExtrapolationOutflow,
+ UBB, FixedDensity, DiffusionDirichlet, SimpleExtrapolationOutflow, WallFunctionBounce,
ExtrapolationOutflow, NeumannByCopy, NoSlip, NoSlipLinearBouzidi, QuadraticBounceBack, StreamInConstant, FreeSlip)
from lbmpy.boundaries.boundaryhandling import LatticeBoltzmannBoundaryHandling
+from lbmpy.boundaries.wall_function_models import MoninObukhovSimilarityTheory, LogLaw, MuskerLaw, SpaldingsLaw
-__all__ = ['NoSlip', 'NoSlipLinearBouzidi', 'QuadraticBounceBack', 'FreeSlip',
+__all__ = ['NoSlip', 'NoSlipLinearBouzidi', 'QuadraticBounceBack', 'FreeSlip', 'WallFunctionBounce',
'UBB', 'FixedDensity',
'SimpleExtrapolationOutflow', 'ExtrapolationOutflow',
'DiffusionDirichlet', 'NeumannByCopy', 'StreamInConstant',
- 'LatticeBoltzmannBoundaryHandling']
+ 'LatticeBoltzmannBoundaryHandling',
+ 'MoninObukhovSimilarityTheory', 'LogLaw', 'MuskerLaw', 'SpaldingsLaw']
diff --git a/src/lbmpy/boundaries/boundaryconditions.py b/src/lbmpy/boundaries/boundaryconditions.py
index 6dcf9d45..3d52af04 100644
--- a/src/lbmpy/boundaries/boundaryconditions.py
+++ b/src/lbmpy/boundaries/boundaryconditions.py
@@ -1,4 +1,5 @@
import abc
+from enum import Enum, auto
from warnings import warn
from pystencils import Assignment, Field
@@ -14,6 +15,7 @@ from lbmpy.maxwellian_equilibrium import discrete_equilibrium
from lbmpy.simplificationfactory import create_simplification_strategy
import sympy as sp
+import numpy as np
class LbBoundary(abc.ABC):
@@ -444,6 +446,241 @@ class FreeSlip(LbBoundary):
# end class FreeSlip
+class WallFunctionBounce(LbBoundary):
+ """
+ Wall function based on the bounce back idea, cf. :cite:`Han2021`. Its implementation is extended to the D3Q27
+ stencil, whereas different weights of the drag distribution are proposed.
+
+ Args:
+ lb_method: LB method which is used for the simulation
+ pdfs: Symbolic representation of the particle distribution functions.
+ normal_direction: Normal direction of the wall. Currently, only straight and axis-aligned walls are supported.
+ wall_function_model: Wall function that is used to retrieve the wall stress :math:`tau_w` during the simulation.
+ See :class:`lbmpy.boundaries.wall_treatment.WallFunctionModel` for more details
+ mean_velocity: Optional field or field access for the mean velocity. As wall functions are typically defined
+ in terms of the mean velocity, it is recommended to provide this variable. Per default, the
+ instantaneous velocity obtained from pdfs is used for the wall function.
+ sampling_shift: Optional sampling shift for the velocity sampling. Can be provided as symbolic variable or
+ integer. In both cases, the user must assure that the sampling shift is at least 1, as sampling
+ in boundary cells is not physical. Per default, a sampling shift of 1 is employed which
+ corresponds to a sampling in the first fluid cell normal to the wall. For lower friction
+ Reynolds numbers, choosing a sampling shift >1 has shown to improve the results for higher
+ resolutions.
+ Mutually exclusive with the Maronga sampling shift.
+ maronga_sampling_shift: Optionally, apply a correction factor to the wall shear stress proposed by Maronga et
+ al. :cite:`Maronga2020`. Has only been tested and validated for the MOST wall function.
+ No guarantee is given that it also works with other wall functions.
+ Mutually exclusive with the standard sampling shift.
+ dt: time discretisation. Usually one in LB units
+ dy: space discretisation. Usually one in LB units
+ y: distance from the wall
+ target_friction_velocity: A target friction velocity can be given if an estimate is known a priori. This target
+ friction velocity will be used as initial guess for implicit wall functions to ensure
+ convergence of the Newton algorithm.
+ weight_method: The extension of the WFB to a D3Q27 stencil is non-unique. Different weights can be chosen to
+ define the drag distribution onto the pdfs. Per default, weights corresponding to the weights
+ in the D3Q27 stencil are chosen.
+ name: Optional name of the boundary.
+ data_type: Floating-point precision. Per default, double.
+ """
+
+ class WeightMethod(Enum):
+ LATTICE_WEIGHT = auto(),
+ GEOMETRIC_WEIGHT = auto()
+
+ def __init__(self, lb_method, pdfs, normal_direction, wall_function_model,
+ mean_velocity=None, sampling_shift=1, maronga_sampling_shift=None,
+ dt=1, dy=1, y=0.5,
+ target_friction_velocity=None,
+ weight_method=WeightMethod.LATTICE_WEIGHT,
+ name=None, data_type='double'):
+ """Set an optional name here, to mark boundaries, for example for force evaluations"""
+ self.stencil = lb_method.stencil
+ if not (self.stencil.Q == 19 or self.stencil.Q == 27):
+ raise ValueError("WFB boundary is currently only defined for D3Q19 and D3Q27 stencils.")
+ self.pdfs = pdfs
+
+ self.wall_function_model = wall_function_model
+
+ if mean_velocity:
+ if isinstance(mean_velocity, Field):
+ self.mean_velocity = mean_velocity.center_vector
+ elif isinstance(mean_velocity, Field.Access):
+ self.mean_velocity = mean_velocity.field.neighbor_vector(mean_velocity.offsets)
+ else:
+ raise ValueError("Mean velocity field has to be a pystencils Field or Field.Access")
+ else:
+ self.mean_velocity = None
+
+ if not isinstance(sampling_shift, int):
+ self.sampling_shift = TypedSymbol(sampling_shift.name, np.uint32)
+ else:
+ assert sampling_shift >= 1, "The sampling shift must be greater than 1."
+ self.sampling_shift = sampling_shift
+
+ if maronga_sampling_shift:
+ assert self.mean_velocity, "Mean velocity field must be provided when using the Maronga correction"
+ if not isinstance(maronga_sampling_shift, int):
+ self.maronga_sampling_shift = TypedSymbol(maronga_sampling_shift.name, np.uint32)
+ else:
+ assert maronga_sampling_shift >= 1, "The Maronga sampling shift must be greater than 1."
+ self.maronga_sampling_shift = maronga_sampling_shift
+ else:
+ self.maronga_sampling_shift = None
+
+ if (self.sampling_shift != 1) and self.maronga_sampling_shift:
+ raise ValueError("Both sampling shift and Maronga offset are set. This is currently not supported.")
+
+ self.dt = dt
+ self.dy = dy
+ self.y = y
+ self.data_type = data_type
+
+ self.target_friction_velocity = target_friction_velocity
+
+ self.weight_method = weight_method
+
+ if len(normal_direction) - normal_direction.count(0) != 1:
+ raise ValueError("Only normal directions for straight walls are supported for example (0, 1, 0) for "
+ "a WallFunctionBounce applied to the southern boundary of the domain")
+
+ self.mirror_axis = normal_direction.index(*[direction for direction in normal_direction if direction != 0])
+
+ self.normal_direction = normal_direction
+ assert all([n in [-1, 0, 1] for n in self.normal_direction]), \
+ "Only -1, 0 and 1 allowed for defining the normal direction"
+ tangential_component = [int(not n) for n in self.normal_direction]
+ self.normal_axis = tangential_component.index(0)
+ self.tangential_axis = [0, 1, 2]
+ self.tangential_axis.remove(self.normal_axis)
+
+ self.dim = self.stencil.D
+
+ if name is None:
+ name = f"WFB : {offset_to_direction_string([-x for x in normal_direction])}"
+
+ super(WallFunctionBounce, self).__init__(name)
+
+ def get_additional_code_nodes(self, lb_method):
+ return [MirroredStencilDirections(self.stencil, self.mirror_axis),
+ NeighbourOffsetArrays(lb_method.stencil)]
+
+ def __call__(self, f_out, f_in, dir_symbol, inv_dir, lb_method, index_field):
+ # needed symbols for offsets and indices
+ # neighbour offset symbols are basically the stencil directions defined in stencils.py:L130ff.
+ neighbor_offset = NeighbourOffsetArrays.neighbour_offset(dir_symbol, lb_method.stencil)
+ tangential_offset = tuple(offset + normal for offset, normal in zip(neighbor_offset, self.normal_direction))
+ mirrored_stencil_symbol = MirroredStencilDirections._mirrored_symbol(self.mirror_axis)
+ mirrored_direction = inv_dir[sp.IndexedBase(mirrored_stencil_symbol, shape=(1,))[dir_symbol]]
+
+ name_base = "f_in_inv_offsets_"
+ offset_array_symbols = [TypedSymbol(name_base + d, mirrored_stencil_symbol.dtype) for d in ['x', 'y', 'z']]
+ mirrored_offset = sp.IndexedBase(mirrored_stencil_symbol, shape=(1,))[dir_symbol]
+ offsets = tuple(sp.IndexedBase(s, shape=(1,))[mirrored_offset] for s in offset_array_symbols)
+
+ # needed symbols in the Assignments
+ u_m = sp.Symbol("u_m")
+ tau_w = sp.Symbol("tau_w")
+ wall_stress = sp.symbols("tau_w_x tau_w_y tau_w_z")
+
+ # if the mean velocity field is not given, or the Maronga correction is applied, density and velocity values
+ # will be calculated from pdfs
+ cqc = lb_method.conserved_quantity_computation
+
+ result = []
+ if (not self.mean_velocity) or self.maronga_sampling_shift:
+ pdf_center_vector = sp.Matrix([0] * self.stencil.Q)
+
+ for i in range(self.stencil.Q):
+ pdf_center_vector[i] = self.pdfs[offsets[0] + self.normal_direction[0],
+ offsets[1] + self.normal_direction[1],
+ offsets[2] + self.normal_direction[2]](i)
+
+ eq_equations = cqc.equilibrium_input_equations_from_pdfs(pdf_center_vector)
+ result.append(eq_equations.all_assignments)
+
+ # sample velocity which will be used in the wall stress calculation
+ if self.mean_velocity:
+ if self.maronga_sampling_shift:
+ u_for_tau_wall = tuple(u_mean_i.get_shifted(
+ self.maronga_sampling_shift * self.normal_direction[0],
+ self.maronga_sampling_shift * self.normal_direction[1],
+ self.maronga_sampling_shift * self.normal_direction[2]
+ ) for u_mean_i in self.mean_velocity)
+ else:
+ u_for_tau_wall = tuple(u_mean_i.get_shifted(
+ self.sampling_shift * self.normal_direction[0],
+ self.sampling_shift * self.normal_direction[1],
+ self.sampling_shift * self.normal_direction[2]
+ ) for u_mean_i in self.mean_velocity)
+
+ rho_for_tau_wall = sp.Float(1)
+ else:
+ rho_for_tau_wall = cqc.density_symbol
+ u_for_tau_wall = cqc.velocity_symbols
+
+ # calculate Maronga factor in case of correction
+ maronga_fix = sp.Symbol("maronga_fix")
+ if self.maronga_sampling_shift:
+ inst_first_cell_vel = cqc.velocity_symbols
+ mean_first_cell_vel = tuple(u_mean_i.get_shifted(*self.normal_direction) for u_mean_i in self.mean_velocity)
+
+ mag_inst_vel_first_cell = sp.sqrt(sum([inst_first_cell_vel[i] ** 2 for i in self.tangential_axis]))
+ mag_mean_vel_first_cell = sp.sqrt(sum([mean_first_cell_vel[i] ** 2 for i in self.tangential_axis]))
+
+ result.append(Assignment(maronga_fix, mag_inst_vel_first_cell / mag_mean_vel_first_cell))
+ else:
+ maronga_fix = 1
+
+ # store which direction is tangential component (only those are used for the wall shear stress)
+ red_u_mag = sp.sqrt(sum([u_for_tau_wall[i]**2 for i in self.tangential_axis]))
+
+ u_mag = Assignment(u_m, red_u_mag)
+ result.append(u_mag)
+
+ wall_distance = self.maronga_sampling_shift if self.maronga_sampling_shift else self.sampling_shift
+
+ # using wall function model
+ wall_law_assignments = self.wall_function_model.shear_stress_assignments(
+ density_symbol=rho_for_tau_wall, velocity_symbol=u_m, shear_stress_symbol=tau_w,
+ wall_distance=(wall_distance - sp.Rational(1, 2) * self.dy),
+ u_tau_target=self.target_friction_velocity)
+ result.append(wall_law_assignments)
+
+ # calculate wall stress components and use them to calculate the drag
+ for i in self.tangential_axis:
+ result.append(Assignment(wall_stress[i], - u_for_tau_wall[i] / u_m * tau_w * maronga_fix))
+
+ weight, inv_weight_sq = sp.symbols("wfb_weight inverse_weight_squared")
+
+ if self.stencil.Q == 19:
+ result.append(Assignment(weight, sp.Rational(1, 2)))
+ elif self.stencil.Q == 27:
+ result.append(Assignment(inv_weight_sq, sum([neighbor_offset[i]**2 for i in self.tangential_axis])))
+ a, b = sp.symbols("wfb_a wfb_b")
+
+ if self.weight_method == self.WeightMethod.LATTICE_WEIGHT:
+ res_ab = sp.solve([2 * a + 4 * b - 1, a - 4 * b], [a, b]) # lattice weight scaling
+ elif self.weight_method == self.WeightMethod.GEOMETRIC_WEIGHT:
+ res_ab = sp.solve([2 * a + 4 * b - 1, a - sp.sqrt(2) * b], [a, b]) # geometric scaling
+ else:
+ raise ValueError("Unknown weighting method for the WFB D3Q27 extension. Currently, only lattice "
+ "weights and geometric weights are supported.")
+
+ result.append(Assignment(weight, sp.Piecewise((sp.Float(0), sp.Equality(inv_weight_sq, 0)),
+ (res_ab[a], sp.Equality(inv_weight_sq, 1)),
+ (res_ab[b], True))))
+
+ factor = self.dt / self.dy * weight
+ drag = sum([neighbor_offset[i] * factor * wall_stress[i] for i in self.tangential_axis])
+
+ result.append(Assignment(f_in.center(inv_dir[dir_symbol]), f_out[tangential_offset](mirrored_direction) - drag))
+
+ return result
+
+# end class WallFunctionBounce
+
+
class UBB(LbBoundary):
r"""Velocity bounce back boundary condition, enforcing specified velocity at obstacle. Furthermore, a density
at the wall can be implied. The boundary condition is implemented with the following formula:
diff --git a/src/lbmpy/boundaries/wall_function_models.py b/src/lbmpy/boundaries/wall_function_models.py
new file mode 100644
index 00000000..4f627787
--- /dev/null
+++ b/src/lbmpy/boundaries/wall_function_models.py
@@ -0,0 +1,173 @@
+import sympy as sp
+from abc import ABC, abstractmethod
+
+from pystencils import Assignment
+
+
+class WallFunctionModel(ABC):
+ def __init__(self, name):
+ self._name = name
+
+ @abstractmethod
+ def shear_stress_assignments(self, density_symbol: sp.Symbol, shear_stress_symbol: sp.Symbol,
+ velocity_symbol: sp.Symbol, wall_distance, u_tau_target):
+ """
+ Computes a symbolic representation for the log law.
+
+ Args:
+ density_symbol: symbol density, should be provided by the LB method's conserved quantity computation
+ shear_stress_symbol: symbolic wall shear stress to which the calculated shear stress will be assigned
+ velocity_symbol: symbolic velocity that is taken as a reference in the wall functions
+ wall_distance: distance to the wall, equals to 0.5 in standard cell-centered LBM
+ u_tau_target: in implicit wall functions, a target friction velocity can be provided which will be used as
+ initial guess in the Newton iteration. This target friction velocity can be obtained, e.g.,
+ from the target friction Reynolds number
+ """
+ pass
+
+
+# end class WallFunctionModel
+
+
+class ExplicitWallFunctionModel(WallFunctionModel, ABC):
+ """
+ Abstract base class for explicit wall functions that can be solved directly for the wall shear stress.
+ """
+
+ def __init__(self, name):
+ super(ExplicitWallFunctionModel, self).__init__(name=name)
+
+
+class MoninObukhovSimilarityTheory(ExplicitWallFunctionModel):
+ def __init__(self, z0, kappa=0.41, phi=0, name="MOST"):
+ self.z0 = z0
+ self.kappa = kappa
+ self.phi = phi
+
+ super(MoninObukhovSimilarityTheory, self).__init__(name=name)
+
+ def shear_stress_assignments(self, density_symbol: sp.Symbol, shear_stress_symbol: sp.Symbol,
+ velocity_symbol: sp.Symbol, wall_distance, u_tau_target=None):
+ u_tau = velocity_symbol * self.kappa / sp.ln(wall_distance / self.z0 + self.phi)
+ return [Assignment(shear_stress_symbol, u_tau ** 2 * density_symbol)]
+
+
+class ImplicitWallFunctionModel(WallFunctionModel, ABC):
+ """
+ Abstract base class for implicit wall functions that require a Newton procedure to solve for the wall shear stress.
+ """
+
+ def __init__(self, name, newton_steps, viscosity):
+ self.newton_steps = newton_steps
+ self.u_tau = sp.symbols(f"wall_function_u_tau_:{self.newton_steps + 1}")
+ self.delta = sp.symbols(f"wall_function_delta_:{self.newton_steps}")
+
+ self.viscosity = viscosity
+
+ super(ImplicitWallFunctionModel, self).__init__(name=name)
+
+ def newton_iteration(self, wall_law):
+ m = -wall_law / wall_law.diff(self.u_tau[0])
+
+ assignments = []
+ for i in range(self.newton_steps):
+ assignments.append(Assignment(self.delta[i], m.subs({self.u_tau[0]: self.u_tau[i]})))
+ assignments.append(Assignment(self.u_tau[i + 1], self.u_tau[i] + self.delta[i]))
+
+ return assignments
+
+
+class LogLaw(ImplicitWallFunctionModel):
+ """
+ Analytical model for the velocity profile inside the boundary layer, obtained from the mean velocity gradient.
+ Only valid in the log-law region.
+ """
+
+ def __init__(self, viscosity, newton_steps=5, kappa=0.41, b=5.2, name="LogLaw"):
+ self.kappa = kappa
+ self.b = b
+
+ super(LogLaw, self).__init__(name=name, newton_steps=newton_steps, viscosity=viscosity)
+
+ def shear_stress_assignments(self, density_symbol: sp.Symbol, shear_stress_symbol: sp.Symbol,
+ velocity_symbol: sp.Symbol, wall_distance, u_tau_target=None):
+ def law(u_p, y_p):
+ return 1 / self.kappa * sp.ln(y_p) + self.b - u_p
+
+ u_plus = velocity_symbol / self.u_tau[0]
+ y_plus = (wall_distance * self.u_tau[0]) / self.viscosity
+
+ u_tau_init = u_tau_target if u_tau_target else velocity_symbol / sp.Float(100)
+
+ wall_law = law(u_plus, y_plus)
+ assignments = [Assignment(self.u_tau[0], u_tau_init), # initial guess
+ *self.newton_iteration(wall_law), # newton iterations
+ Assignment(shear_stress_symbol, self.u_tau[-1] ** 2 * density_symbol)] # final result
+
+ return assignments
+
+
+class SpaldingsLaw(ImplicitWallFunctionModel):
+ """
+ Single formula for the velocity profile inside the boundary layer, proposed by Spalding :cite:`spalding1961`.
+ Valid in the inner and the outer layer.
+ """
+
+ def __init__(self, viscosity, newton_steps=5, kappa=0.41, b=5.5, name="Spalding"):
+ self.kappa = kappa
+ self.b = b
+
+ super(SpaldingsLaw, self).__init__(name=name, newton_steps=newton_steps, viscosity=viscosity)
+
+ def shear_stress_assignments(self, density_symbol: sp.Symbol, shear_stress_symbol: sp.Symbol,
+ velocity_symbol: sp.Symbol, wall_distance, u_tau_target=None):
+ def law(u_p, y_p):
+ k_times_u = self.kappa * u_p
+ frac_1 = (k_times_u ** 2) / sp.Float(2)
+ frac_2 = (k_times_u ** 3) / sp.Float(6)
+ return (u_p + sp.exp(-self.kappa * self.b) * (sp.exp(k_times_u) - sp.Float(1) - k_times_u - frac_1 - frac_2)
+ - y_p)
+
+ u_plus = velocity_symbol / self.u_tau[0]
+ y_plus = (wall_distance * self.u_tau[0]) / self.viscosity
+
+ u_tau_init = u_tau_target if u_tau_target else velocity_symbol / sp.Float(100)
+
+ wall_law = law(u_plus, y_plus)
+ assignments = [Assignment(self.u_tau[0], u_tau_init), # initial guess
+ *self.newton_iteration(wall_law), # newton iterations
+ Assignment(shear_stress_symbol, self.u_tau[-1] ** 2 * density_symbol)] # final result
+
+ return assignments
+
+
+class MuskerLaw(ImplicitWallFunctionModel):
+ """
+ Quasi-analytical model for the velocity profile inside the boundary layer, proposed by Musker. Valid in the inner
+ and the outer layer.
+ Formulation taken from :cite:`malaspinas2015`, Equation (59).
+ """
+
+ def __init__(self, viscosity, newton_steps=5, name="Musker"):
+
+ super(MuskerLaw, self).__init__(name=name, newton_steps=newton_steps, viscosity=viscosity)
+
+ def shear_stress_assignments(self, density_symbol: sp.Symbol, shear_stress_symbol: sp.Symbol,
+ velocity_symbol: sp.Symbol, wall_distance, u_tau_target=None):
+ def law(u_p, y_p):
+ arctan = sp.Float(5.424) * sp.atan(sp.Float(0.119760479041916168) * y_p - sp.Float(0.488023952095808383))
+ logarithm = (sp.Float(0.434) * sp.log((y_p + sp.Float(10.6)) ** sp.Float(9.6)
+ / (y_p ** 2 - sp.Float(8.15) * y_p + sp.Float(86)) ** 2, 10))
+ return (arctan + logarithm - sp.Float(3.50727901936264842)) - u_p
+
+ u_plus = velocity_symbol / self.u_tau[0]
+ y_plus = (wall_distance * self.u_tau[0]) / self.viscosity
+
+ u_tau_init = u_tau_target if u_tau_target else velocity_symbol / sp.Float(100)
+
+ wall_law = law(u_plus, y_plus)
+ assignments = [Assignment(self.u_tau[0], u_tau_init), # initial guess
+ *self.newton_iteration(wall_law), # newton iterations
+ Assignment(shear_stress_symbol, self.u_tau[-1] ** 2 * density_symbol)] # final result
+
+ return assignments
diff --git a/src/lbmpy/flow_statistics.py b/src/lbmpy/flow_statistics.py
new file mode 100644
index 00000000..3e5b5f5f
--- /dev/null
+++ b/src/lbmpy/flow_statistics.py
@@ -0,0 +1,82 @@
+import sympy as sp
+import pystencils as ps
+
+from pystencils.field import Field
+
+
+def welford_assignments(field, mean_field, sum_of_products_field=None):
+ r"""
+ Creates the assignments for the kernel creation of Welford's algorithm
+ (https://en.wikipedia.org/wiki/Algorithms_for_calculating_variance#Welford's_online_algorithm).
+ This algorithm is an online algorithm for the mean and variance calculation of sample data, here implemented for
+ the execution on vector fields, e.g., velocity.
+ The calculation of the variance is optional and only performed if a field for the sum of products is given.
+
+ The mean value is directly updated in the mean vector field.
+ The variance must be retrieved in a post-processing step. Let :math `M_{2,n}` denote the value of the sum of
+ products after the first :math `n` samples. According to Welford the biased sample variance
+ :math `\sigma_n^2` and the unbiased sample variance :math `s_n^2` are given by
+
+ .. math ::
+ \sigma_n^2 = \frac{M_{2,n}}{n}
+ s_n^2 = \frac{M_{2,n}}{n-1},
+
+ respectively.
+ """
+
+ if isinstance(field, Field):
+ dim = field.values_per_cell()
+ welford_field = field.center
+ elif isinstance(field, Field.Access):
+ dim = field.field.values_per_cell()
+ welford_field = field
+ else:
+ raise ValueError("Vector field has to be a pystencils Field or Field.Access")
+
+ if isinstance(mean_field, Field):
+ welford_mean_field = mean_field.center
+ elif isinstance(mean_field, Field.Access):
+ welford_mean_field = mean_field
+ else:
+ raise ValueError("Mean vector field has to be a pystencils Field or Field.Access")
+
+ if sum_of_products_field is None:
+ # sum of products will not be calculated, i.e., the variance is not retrievable
+ welford_sum_of_products_field = None
+ else:
+ if isinstance(sum_of_products_field, Field):
+ welford_sum_of_products_field = sum_of_products_field.center
+ assert sum_of_products_field.values_per_cell() == dim**2, \
+ (f"Sum of products field does not have the right layout. "
+ f"Index dimension: {sum_of_products_field.values_per_cell()}, expected: {dim**2}")
+ elif isinstance(sum_of_products_field, Field.Access):
+ welford_sum_of_products_field = sum_of_products_field
+ assert sum_of_products_field.field.values_per_cell() == dim**2, \
+ (f"Sum of products field does not have the right layout. "
+ f"Index dimension: {sum_of_products_field.field.values_per_cell()}, expected: {dim**2}")
+ else:
+ raise ValueError("Variance vector field has to be a pystencils Field or Field.Access")
+
+ counter = sp.Symbol('counter')
+ delta = sp.symbols(f"delta_:{dim}")
+
+ main_assignments = list()
+
+ for i in range(dim):
+ main_assignments.append(ps.Assignment(delta[i], welford_field.at_index(i) - welford_mean_field.at_index(i)))
+ main_assignments.append(ps.Assignment(welford_mean_field.at_index(i),
+ welford_mean_field.at_index(i) + delta[i] / counter))
+
+ if sum_of_products_field is not None:
+ delta2 = sp.symbols(f"delta2_:{dim}")
+ for i in range(dim):
+ main_assignments.append(
+ ps.Assignment(delta2[i], welford_field.at_index(i) - welford_mean_field.at_index(i)))
+ for i in range(dim):
+ for j in range(dim):
+ idx = i * dim + j
+ main_assignments.append(ps.Assignment(
+ welford_sum_of_products_field.at_index(idx),
+ welford_sum_of_products_field.at_index(idx) + delta[i] * delta2[j]))
+
+ return main_assignments
diff --git a/tests/test_wall_function_bounce.py b/tests/test_wall_function_bounce.py
new file mode 100644
index 00000000..9eeeb1f7
--- /dev/null
+++ b/tests/test_wall_function_bounce.py
@@ -0,0 +1,116 @@
+import pytest
+
+import pystencils as ps
+from lbmpy import Stencil, LBStencil, LBMConfig, Method, lattice_viscosity_from_relaxation_rate, \
+ LatticeBoltzmannStep, pdf_initialization_assignments, ForceModel
+from lbmpy.boundaries.boundaryconditions import UBB, WallFunctionBounce, NoSlip, FreeSlip
+from lbmpy.boundaries.wall_function_models import SpaldingsLaw, LogLaw, MoninObukhovSimilarityTheory, MuskerLaw
+from pystencils.slicing import slice_from_direction
+
+
+@pytest.mark.parametrize('stencil', [Stencil.D3Q19, Stencil.D3Q27])
+@pytest.mark.parametrize('wfb_type', ['wfb_i', 'wfb_ii', 'wfb_iii', 'wfb_iv'])
+def test_wfb(stencil, wfb_type):
+ stencil = LBStencil(stencil)
+ periodicity = (True, False, True)
+ domain_size = (30, 20, 15)
+ dim = len(domain_size)
+
+ omega = 1.1
+ nu = lattice_viscosity_from_relaxation_rate(omega)
+
+ # pressure gradient for laminar channel flow with u_max = 0.1
+ ext_force_density = 0.1 * 2 * nu / ((domain_size[1] - 2) / 2) ** 2
+
+ lbm_config = LBMConfig(stencil=stencil,
+ method=Method.SRT,
+ force_model=ForceModel.GUO,
+ force=(ext_force_density, 0, 0),
+ relaxation_rate=omega,
+ compressible=True)
+
+ wall_north = NoSlip()
+ normal = (0, 1, 0)
+
+ # NO-SLIP
+
+ lb_step_noslip = LatticeBoltzmannStep(domain_size=domain_size, periodicity=periodicity,
+ lbm_config=lbm_config, compute_velocity_in_every_step=True)
+
+ noslip = NoSlip()
+
+ lb_step_noslip.boundary_handling.set_boundary(noslip, slice_from_direction('S', dim))
+ lb_step_noslip.boundary_handling.set_boundary(wall_north, slice_from_direction('N', dim))
+
+ # FREE-SLIP
+
+ lb_step_freeslip = LatticeBoltzmannStep(domain_size=domain_size, periodicity=periodicity,
+ lbm_config=lbm_config, compute_velocity_in_every_step=True)
+
+ freeslip = FreeSlip(stencil=stencil, normal_direction=normal)
+
+ lb_step_freeslip.boundary_handling.set_boundary(freeslip, slice_from_direction('S', dim))
+ lb_step_freeslip.boundary_handling.set_boundary(wall_north, slice_from_direction('N', dim))
+
+ # WFB
+
+ lb_step_wfb = LatticeBoltzmannStep(domain_size=domain_size, periodicity=periodicity,
+ lbm_config=lbm_config, compute_velocity_in_every_step=True)
+
+ # pdf initialisation
+
+ init = pdf_initialization_assignments(lb_step_wfb.method, 1.0, (0.025, 0, 0),
+ lb_step_wfb.data_handling.fields[lb_step_wfb._pdf_arr_name].center_vector)
+
+ config = ps.CreateKernelConfig(target=lb_step_wfb.data_handling.default_target, cpu_openmp=False)
+ ast_init = ps.create_kernel(init, config=config)
+ kernel_init = ast_init.compile()
+
+ lb_step_wfb.data_handling.run_kernel(kernel_init)
+
+ # potential mean velocity field
+
+ mean_vel_field = lb_step_wfb.data_handling.fields[lb_step_wfb.velocity_data_name]
+ # mean_vel_field = lb_step_wfb.data_handling.add_array('mean_velocity_field', values_per_cell=stencil.D)
+ # lb_step_wfb.data_handling.fill('mean_velocity_field', 0.005, value_idx=0, ghost_layers=True)
+ lb_step_wfb.data_handling.fill(lb_step_wfb.velocity_data_name, 0.025, value_idx=0, ghost_layers=True)
+
+ # wfb arguments
+ wfb_args = {
+ 'wfb_i': {'wall_function_model': SpaldingsLaw(viscosity=nu),
+ 'weight_method': WallFunctionBounce.WeightMethod.GEOMETRIC_WEIGHT,
+ 'name': "wall"},
+ 'wfb_ii': {'wall_function_model': MuskerLaw(viscosity=nu),
+ 'weight_method': WallFunctionBounce.WeightMethod.GEOMETRIC_WEIGHT,
+ 'mean_velocity': mean_vel_field,
+ 'name': "wall"},
+ 'wfb_iii': {'wall_function_model': LogLaw(viscosity=nu),
+ 'weight_method': WallFunctionBounce.WeightMethod.LATTICE_WEIGHT,
+ 'mean_velocity': mean_vel_field.center,
+ 'sampling_shift': 2},
+ 'wfb_iv': {'wall_function_model': MoninObukhovSimilarityTheory(z0=1e-2),
+ 'weight_method': WallFunctionBounce.WeightMethod.LATTICE_WEIGHT,
+ 'mean_velocity': mean_vel_field,
+ 'maronga_sampling_shift': 2}
+ }
+
+ wall = WallFunctionBounce(lb_method=lb_step_wfb.method, normal_direction=normal,
+ pdfs=lb_step_wfb.data_handling.fields[lb_step_wfb._pdf_arr_name],
+ **wfb_args[wfb_type])
+
+ lb_step_wfb.boundary_handling.set_boundary(wall, slice_from_direction('S', dim))
+ lb_step_wfb.boundary_handling.set_boundary(wall_north, slice_from_direction('N', dim))
+
+ # rum cases
+
+ timesteps = 4000
+ lb_step_noslip.run(timesteps)
+ lb_step_freeslip.run(timesteps)
+ lb_step_wfb.run(timesteps)
+
+ noslip_velocity = lb_step_noslip.velocity[domain_size[0] // 2, :, domain_size[2] // 2, 0]
+ freeslip_velocity = lb_step_freeslip.velocity[domain_size[0] // 2, :, domain_size[2] // 2, 0]
+ wfb_velocity = lb_step_wfb.velocity[domain_size[0] // 2, :, domain_size[2] // 2, 0]
+
+ assert wfb_velocity[0] > noslip_velocity[0], f"WFB enforced velocity below no-slip velocity"
+ assert wfb_velocity[0] < freeslip_velocity[0], f"WFB enforced velocity above free-slip velocity"
--
GitLab