Skip to content

Closes #63 #64

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
nasserma merged 1 commit into from
Sep 22, 2026
Merged

Closes #63 #64

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
62 changes: 49 additions & 13 deletions opencmp/helpers/wall_func.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -91,8 +91,8 @@ def __init__(self, mesh: ngs.comp.Mesh, nu: float, C_mu: float, kappa: float,
def _mark_wall_cells(self) -> None:
"""Find the wall layer from mesh topology.

Sets ``_wall_measure``, ``_marked`` / ``wall_facet_cell_mask`` and
``mask`` to the cells that own a physical wall facet.
``_marked`` / ``wall_facet_cell_mask`` identify physical facet owners.
``mask`` additionally includes cells touching the wall at a vertex.
"""
# Assumes a simplicial mesh (one L2(0) DOF per cell). Quads/hexes parse
# but the resulting layer is untested, so refuse rather than go silently wrong.
Expand Down Expand Up @@ -121,16 +121,48 @@ def _mark_wall_cells(self) -> None:
self.wall_facet_cell_mask = ngs.GridFunction(self._fes0)
self.wall_facet_cell_mask.vec.FV().NumPy()[:] = self._marked.astype(float)

# A wall function is a first-cell treatment. Cells that do not own a
# physical wall facet use the bulk closure, even when facet-adjacent to
# a wall owner.
volume_elements = list(self._fes0.Elements(ngs.VOL))
wall_vertices = {
vertex.nr
for boundary_element in self.mesh.Elements(ngs.BND)
if boundary_element.mat == self.wall_boundary
for vertex in boundary_element.vertices
}
vertex_marked = np.zeros_like(self._marked)
for element in volume_elements:
if any(vertex.nr in wall_vertices for vertex in element.vertices):
vertex_marked[list(element.dofs)] = True

self.mask = ngs.GridFunction(self._fes0)
self.mask.vec.FV().NumPy()[:] = self._marked.astype(float)
self._wall_layer_sources = {}
self.mask.vec.FV().NumPy()[:] = vertex_marked.astype(float)
self._element_dofs = {
element.nr: tuple(element.dofs)
for element in self._fes0.Elements(ngs.VOL)
for element in volume_elements
}
owner_numbers = {
element.nr for element in volume_elements
if any(self._marked[dof] for dof in element.dofs)
}
vertex_to_owners = {}
for element in volume_elements:
if element.nr not in owner_numbers:
continue
for vertex in element.vertices:
if vertex.nr in wall_vertices:
vertex_to_owners.setdefault(vertex.nr, set()).add(element.nr)

self._wall_layer_sources = {}
for element in volume_elements:
if element.nr in owner_numbers or not any(
vertex_marked[dof] for dof in element.dofs):
continue
sources = set()
for vertex in element.vertices:
sources.update(vertex_to_owners.get(vertex.nr, ()))
if not sources:
raise ValueError(
f'Wall-vertex cell {element.nr} has no wall-facet owner source.')
self._wall_layer_sources[element.nr] = tuple(sorted(sources))
self._wall_facets = self._find_physical_wall_facets()

def _find_physical_wall_facets(self):
Expand Down Expand Up @@ -254,19 +286,23 @@ def update(self, K, U=None) -> None:
wall_nu_t = np.zeros(self._fes0.ndof)
wall_shear = np.zeros(self._fes0.ndof)
y_plus_cell = np.zeros(self._fes0.ndof)
active = self.mask.vec.FV().NumPy() > 0.5
if self.u_tau_method == 0:
projected = ngs.GridFunction(self._fes0)
projected.Set(K)
k_values = np.maximum(projected.vec.FV().NumPy(), 0.0)
u_tau[self._marked] = self.C_mu ** 0.25 * np.sqrt(
k_values[self._marked])
u_tau[active] = self.C_mu ** 0.25 * np.sqrt(k_values[active])
y_plus_cell = (self._dist_cell.vec.FV().NumPy()
* u_tau / self.nu)
else:
if U is None:
raise ValueError('Velocity-based u_tau requires the velocity iterate.')
wall_shear = self._resolved_wall_shear(U)
u_tau[self._marked] = np.sqrt(wall_shear[self._marked])
for element_number, sources in self._wall_layer_sources.items():
target_dof = self._element_dofs[element_number][0]
source_dofs = [self._element_dofs[source][0] for source in sources]
wall_shear[target_dof] = float(np.mean(wall_shear[source_dofs]))
u_tau[active] = np.sqrt(wall_shear[active])
y_plus_cell = (self._dist_cell.vec.FV().NumPy()
* u_tau / self.nu)
wall_nu_t = self._wall_viscosity_from_yplus(y_plus_cell)
Expand All @@ -281,7 +317,7 @@ def update(self, K, U=None) -> None:

def _warn_if_yplus_outside_recommended_range(self, y_plus: np.ndarray) -> None:
"""Warn once for each side of the recommended wall-function band."""
wall_values = y_plus[self._marked]
wall_values = y_plus[self.mask.vec.FV().NumPy() > 0.5]
total = wall_values.size
observed_min = float(np.min(wall_values))
observed_max = float(np.max(wall_values))
Expand Down Expand Up @@ -315,7 +351,7 @@ def _warn_if_yplus_outside_recommended_range(self, y_plus: np.ndarray) -> None:
# ------------------------------------------------------------------

def near_wall_mask(self) -> ngs.GridFunction:
"""1 only on cells that own a physical wall facet."""
"""1 on cells owning a wall facet or touching the wall at a vertex."""
return self.mask

def wall_facet_mask(self) -> ngs.GridFunction:
Expand Down
38 changes: 26 additions & 12 deletions pytests/helpers/test_k_epsilon_wall_func.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -53,6 +53,17 @@ def elements_owning_a_facet_on(mesh, wall):
if sum(1 for v in el.vertices if v.nr in wall_vertices) >= mesh.dim}


def elements_touching_wall_vertices(mesh, wall):
wall_vertices = {
vertex.nr
for boundary_element in mesh.Elements(ngs.BND)
if boundary_element.mat == wall
for vertex in boundary_element.vertices
}
return {element.nr for element in mesh.Elements(ngs.VOL)
if any(vertex.nr in wall_vertices for vertex in element.vertices)}


def add_one_face_connected_layer(mesh, element_numbers):
"""Ground truth for one topological dilation through volume-cell facets."""
elements = list(mesh.Elements(ngs.VOL))
Expand All @@ -74,20 +85,20 @@ def add_one_face_connected_layer(mesh, element_numbers):
# ----------------------------------------------------------------------

@pytest.mark.parametrize('meshfile, wall', [(SQUARE, 'bottom'), (CHANNEL, 'wall')])
def test_mask_contains_only_wall_facet_owners(meshfile, wall):
def test_mask_contains_every_cell_touching_a_wall_vertex(meshfile, wall):
mesh = ngs.Mesh(meshfile)
wf = build(mesh, wall)
wall_cells = elements_owning_a_facet_on(mesh, wall)
wall_cells = elements_touching_wall_vertices(mesh, wall)
assert marked_element_numbers(wf) == wall_cells


def test_near_wall_mask_matches_wall_facet_owners():
def test_near_wall_mask_matches_wall_vertex_cells():
mesh = ngs.Mesh(CHANNEL)
wf = build(mesh, 'wall')
values = wf.near_wall_mask().vec.FV().NumPy()
actual = {element.nr for element in wf._fes0.Elements(ngs.VOL)
if any(values[dof] > 0.5 for dof in element.dofs)}
wall_cells = elements_owning_a_facet_on(mesh, 'wall')
wall_cells = elements_touching_wall_vertices(mesh, 'wall')
assert actual == wall_cells


Expand Down Expand Up @@ -117,8 +128,7 @@ def test_wall_measure_sums_to_the_exact_boundary_measure():
assert wf._wall_measure.sum() == pytest.approx(exact, rel=1e-12)


def test_cell_touching_the_wall_only_at_a_vertex_is_not_marked():
"""A vertex-only touch is not a wall-function cell."""
def test_cell_touching_the_wall_only_at_a_vertex_is_marked():
mesh = ngs.Mesh(SQUARE)
wf = build(mesh, 'bottom')
marked = marked_element_numbers(wf)
Expand All @@ -131,7 +141,7 @@ def test_cell_touching_the_wall_only_at_a_vertex_is_not_marked():
vertex_only = {el.nr for el in mesh.Elements(ngs.VOL)
if sum(1 for v in el.vertices if v.nr in bottom_vertices) == 1}
assert vertex_only, 'mesh exercises no vertex-only touch; test is vacuous'
assert not vertex_only & marked
assert vertex_only <= marked


def test_missing_wall_marker_raises_a_clear_error():
Expand Down Expand Up @@ -360,11 +370,12 @@ def test_bulk_viscosity_matches_the_plain_k_epsilon_formula(channel_wf):
assert nu_t[~marked] == pytest.approx(0.09 * k ** 2 / eps, rel=1e-6)


def test_wall_owner_neighbours_are_not_in_the_wall_function_mask(channel_wf):
def test_only_vertex_connected_owner_neighbours_join_the_wall_mask(channel_wf):
mesh, wf = channel_wf
owners = elements_owning_a_facet_on(mesh, 'wall')
expanded = add_one_face_connected_layer(mesh, owners)
neighbours = expanded - owners
vertex_cells = elements_touching_wall_vertices(mesh, 'wall') - owners
assert neighbours

k = ngs.GridFunction(wf._fes0)
Expand All @@ -374,8 +385,11 @@ def test_wall_owner_neighbours_are_not_in_the_wall_function_mask(channel_wf):
u_tau = wf.u_tau_cell.vec.FV().NumPy()
for element in wf._fes0.Elements(ngs.VOL):
if element.nr in neighbours:
assert mask[list(element.dofs)] == pytest.approx(0.0, abs=1e-14)
assert u_tau[list(element.dofs)] == pytest.approx(0.0, abs=1e-14)
expected = 1.0 if element.nr in vertex_cells else 0.0
assert mask[list(element.dofs)] == pytest.approx(expected, abs=1e-14)
expected_utau = 0.09 ** 0.25 if expected else 0.0
assert u_tau[list(element.dofs)] == pytest.approx(
expected_utau, abs=1e-14)


def test_velocity_method_uses_molecular_tangential_not_normal_stress():
Expand Down Expand Up @@ -455,9 +469,9 @@ def test_tet_mesh_attributes_every_wall_face_to_exactly_one_cell(cube_wf):
assert wf._wall_measure.sum() == pytest.approx(exact, rel=1e-12)


def test_tet_mask_contains_only_wall_face_owners(cube_wf):
def test_tet_mask_contains_every_wall_vertex_cell(cube_wf):
mesh, wf = cube_wf
wall_cells = elements_owning_a_facet_on(mesh, 'bottom')
wall_cells = elements_touching_wall_vertices(mesh, 'bottom')
assert wall_cells, 'no cell owns a wall face; test is vacuous'
assert marked_element_numbers(wf) == wall_cells

Expand Down
Loading