diff --git a/opencmp/helpers/wall_func.py b/opencmp/helpers/wall_func.py index afe2416..70ba95e 100644 --- a/opencmp/helpers/wall_func.py +++ b/opencmp/helpers/wall_func.py @@ -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. @@ -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): @@ -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) @@ -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)) @@ -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: diff --git a/pytests/helpers/test_k_epsilon_wall_func.py b/pytests/helpers/test_k_epsilon_wall_func.py index 1739402..67a92df 100644 --- a/pytests/helpers/test_k_epsilon_wall_func.py +++ b/pytests/helpers/test_k_epsilon_wall_func.py @@ -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)) @@ -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 @@ -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) @@ -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(): @@ -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) @@ -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(): @@ -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