diff --git a/docs/developer/subsystems/stress-transport.md b/docs/developer/subsystems/stress-transport.md index 42242292..0459a49b 100644 --- a/docs/developer/subsystems/stress-transport.md +++ b/docs/developer/subsystems/stress-transport.md @@ -24,7 +24,7 @@ stokes.constitutive_model.Parameters.dt_elastic = dt | `semi_lagrangian` (nodal) | continuous P1 at the vertices | vertex trace-back, interpolation at the foot | any Courant number | excess stress in the first cells off a no-slip wall; on the confined cylinder that excess loses the conformation and the solve hangs | | `integration_point` | continuous P1 store, sampled at the quadrature points | trace-back of every quadrature point | Courant near one, or below one with store smoothing | a cell-scale mode of the stress that grows below Courant one when the solvent viscosity is small | | `forward` | discontinuous P1 per cell, fitted from the arrivals | fixed launch set of interior points (the integration points), one forward trajectory a step; the flux is read back at the launch points through a continuous P1 projection; an inflow cell's uncovered share is filled with the inflow value | the cylinder walls at dt 0.04; below Courant one with `flux_smoothing` at c = 0.023 (Waters-King 1/16, dt 0.0125: 0.9543 at t 1 and 0.5185 at t 6.5, against nodal 0.9622 and 0.5171) | the same cell-scale mode as the integration-point history without that smoothing (diverges at t 2.4 there); first order only; does not cross a periodic seam or follow a moving mesh | -| `lagrangian` (particles) | a swarm the solver owns and advects, one value per particle, read through a discontinuous cells proxy | the material points themselves: the constitutive flux is evaluated at the particles each step and never projected back to the mesh | any Courant number; no numerical diffusion of the history | the cost and bookkeeping of a swarm, and a proxy that needs its cells kept populated (population control refills them); the conformation check does not read a per-point tensor from it | +| `lagrangian` (particles) | a swarm the solver owns and advects, one value per particle, read through a discontinuous cells proxy | the material points themselves: the constitutive flux is evaluated at the particles each step and never projected back to the mesh; a particle that entered through an inflow takes the inflow value | any Courant number; no numerical diffusion of the history | the cost and bookkeeping of a swarm, and a proxy that needs its cells kept populated (population control refills them); the conformation check does not read a per-point tensor from it | | `eulerian` (SUPG grid) | continuous P1 | assembled transport equation with streamline upwinding | with DEVSS | without DEVSS the velocity block loses its preconditioner as the stress grows | The nodal and integration-point flavours store the stress after every solve by diff --git a/src/underworld3/swarm.py b/src/underworld3/swarm.py index 7fc16b64..b24792ce 100644 --- a/src/underworld3/swarm.py +++ b/src/underworld3/swarm.py @@ -5465,7 +5465,9 @@ def repopulate( Returns ------- - (added, removed) : the counts on this rank. + (added, removed) : the counts on this rank. New particles follow the + existing ones in storage order, so the last ``added`` rows of every + variable are the ones just created. """ mesh = self.mesh dim = self.cdim @@ -5512,12 +5514,22 @@ def repopulate( surplus = int(npc[c] - max_per_cell) drop.extend(idx[np.argsort(nearest_all[idx])[:surplus]].tolist()) if drop: + # PETSc removes a point by copying the LAST point into its + # slot (DMSwarmDataBucketRemovePointAtIndex), so the survivors + # are reordered in storage. The reconstruction below reads + # values in storage order, and the tree it is built on must + # share it: the same moves are replayed on the local copies + # (#784). Descending order, so a moved point is never one + # still to be removed. + last = X.shape[0] for index in sorted(drop, reverse=True): self.dm.removePointAtIndex(int(index)) + last -= 1 + if index != last: + X[index] = X[last] + cells[index] = cells[last] + X, cells = X[:last], cells[:last] removed = len(drop) - keep = np.ones(X.shape[0], dtype=bool) - keep[drop] = False - X, cells = X[keep], cells[keep] npc = np.bincount(cells[cells >= 0], minlength=ncells) self._invalidate_canonical_data() diff --git a/src/underworld3/systems/ddt.py b/src/underworld3/systems/ddt.py index 854a4f02..4c1aa4e7 100644 --- a/src/underworld3/systems/ddt.py +++ b/src/underworld3/systems/ddt.py @@ -1049,12 +1049,16 @@ def inflow_value(self): :class:`EulerianSUPG` compiles the value into a boundary term of its transport solve; :class:`IntegrationPointSemiLagrangian` gives it to a departure point restored to the boundary; :class:`ForwardSemiLagrangian` - fills the uncovered share of an inflow cell with it. The nodal - trace-back and the particle flavours do not use it: a departure point - or a particle that lands outside the domain is restored to the - boundary and takes the transported field's value THERE, which - constrains the inflow but is not the value set. Setting a value on - such a flavour says so once rather than dropping it in silence (#733). + fills the uncovered share of an inflow cell with it; + :class:`Lagrangian` gives it to every particle that entered through + an inflow: one whose back-trace over the step, or over one cell for a + particle the refill created, leaves the domain there (#783). The nodal + trace-back and :class:`Lagrangian_Swarm` (a swarm the caller advects) + do not use it: a departure point or a particle that lands outside the + domain is restored to the boundary and takes the transported field's + value THERE, which constrains the inflow but is not the value set. + Setting a value on such a flavour says so once rather than dropping it + in silence (#733). """ return self._inflow_value @@ -1072,14 +1076,54 @@ def inflow_value(self, value): "restores an out-of-bounds departure point to the boundary " "and reads the transported field there, which constrains " "the inflow but is not the value you set. EulerianSUPG, " - "IntegrationPointSemiLagrangian and ForwardSemiLagrangian " - "apply it (#733).", + "IntegrationPointSemiLagrangian, ForwardSemiLagrangian and " + "Lagrangian apply it (#733, #783).", stacklevel=2) self._inflow_value = value #: Whether this flavour compiles :attr:`inflow_value` into its transport. applies_inflow_value = False + def _nondim_timestep(self, dt): + r"""Reduce ``dt`` to a plain non-dimensional model-time value. + + The semi-Lagrangian trace-back is performed ENTIRELY in the mesh's + NON-DIMENSIONAL (DM) coordinate space: evaluate()/global_evaluate + treat plain arrays as DM coords and the DM point-location uses DM + values (0..L_model, NOT dimensional metres). So coords, velocity + AND dt are all reduced to non-dimensional values, whether or not + the model carries units. (Previously the has_units branch kept + dimensional coords/velocity and left dt unitless -> a 'meter' vs + 'meter/second' subtraction crash and mislocation against the ND + DM; UW3 issue #267.) + """ + if hasattr(dt, "magnitude") or hasattr(dt, "value"): + # dt carries units -> non-dimensionalise it + dt_nondim = uw.non_dimensionalise(dt, uw.get_default_model()) + if hasattr(dt_nondim, "magnitude"): + return float(dt_nondim.magnitude) + elif hasattr(dt_nondim, "value"): + return float(dt_nondim.value) + else: + return float(dt_nondim) + else: + # already non-dimensional model-time + return dt + + + def _write_inflow(self, var, coords, rows): + """Overwrite ``rows`` of ``var`` with :attr:`inflow_value` evaluated at + ``coords`` (the positions of ALL the points, so that the read, which + is collective when the value holds a field, is made on every rank; + only ``rows`` are written). Storage is non-dimensional, so the value + is reduced through the history's units. A flavour that calls this + sets ``_components`` (its stored columns) and ``_psi_units``.""" + expr = self._inflow_value + for column, (i, j) in enumerate(self._components): + vals = uw.function.evaluate(expr[i, j], coords) + var.data[rows, column] = np.asarray( + _to_nondim_ndarray(vals, units=self._psi_units)).reshape(-1)[rows] + def _unknown_shape(self): """Shape of the unknown as a matrix (``Symbolic`` stores ``_shape`` as data).""" psi = self.psi_fn @@ -3513,32 +3557,6 @@ def _record_current_field_into_history( if i != j: self.psi_star[0].array[:, j, i] = vals - def _nondim_timestep(self, dt): - r"""Reduce ``dt`` to a plain non-dimensional model-time value. - - The semi-Lagrangian trace-back is performed ENTIRELY in the mesh's - NON-DIMENSIONAL (DM) coordinate space: evaluate()/global_evaluate - treat plain arrays as DM coords and the DM point-location uses DM - values (0..L_model, NOT dimensional metres). So coords, velocity - AND dt are all reduced to non-dimensional values, whether or not - the model carries units. (Previously the has_units branch kept - dimensional coords/velocity and left dt unitless -> a 'meter' vs - 'meter/second' subtraction crash and mislocation against the ND - DM; UW3 issue #267.) - """ - if hasattr(dt, "magnitude") or hasattr(dt, "value"): - # dt carries units -> non-dimensionalise it - dt_nondim = uw.non_dimensionalise(dt, uw.get_default_model()) - if hasattr(dt_nondim, "magnitude"): - return float(dt_nondim.magnitude) - elif hasattr(dt_nondim, "value"): - return float(dt_nondim.value) - else: - return float(dt_nondim) - else: - # already non-dimensional model-time - return dt - def _trace_departure_points( self, i, node_coords_nd, dt_for_calc, evalf, subtract_v_mesh, oldframe_active ): @@ -3956,6 +3974,10 @@ class Lagrangian(_DDtBase): commits_flux_in_post_solve = True + #: A particle that entered through an inflow this step takes + #: :attr:`inflow_value` (see :meth:`_apply_inflow_value`). + applies_inflow_value = True + instances = ( 0 # count how many of these there are in order to create unique private mesh variable ids ) @@ -3989,6 +4011,13 @@ def __init__( self.V_fn = V_fn self.verbose = verbose self.order = order + # Particle storage is non-dimensional; an inflow datum with units is + # reduced through these before it is written (#783). + psi_units = uw.get_units(psi_fn) + if psi_units is not None and not uw.get_default_model().has_units(): + psi_units = None + self._psi_units = psi_units + self._components = _storage_components(vtype, tuple(sympy.Matrix(psi_fn).shape)) self._init_history_tracking(order) @@ -4022,10 +4051,13 @@ def __init__( # the swarm without bound where particles pile up against a wall): the # bounds are set from the initial occupancy, with a floor at the linear # fit minimum. + # The manager applies the control itself after each advection rather + # than leaving it on the swarm, so it knows which particles the refill + # created: those are the ones an inflow datum must reach (#783). _npart = int(uw.mpi.comm.allreduce(np.asarray(dudt_swarm._particle_coordinates.data).shape[0], op=uw.MPI.SUM)) _ncell = int(uw.mpi.comm.allreduce(mesh._centroids.shape[0], op=uw.MPI.SUM)) _mean = _npart / max(_ncell, 1) - dudt_swarm.population_control = { + self._population_control = { "min_per_cell": max(mesh.dim + 1, int(0.5 * _mean)), "max_per_cell": max(2 * (mesh.dim + 1), int(3.0 * _mean)), } @@ -4159,6 +4191,70 @@ def update_pre_solve( return + def _apply_inflow_value(self, dt, created): + r"""Give every particle that entered this step the inflow datum. + + No particle arrives from outside: the population control CREATES the + particles of an emptied inlet cell, at lattice points anywhere in the + cell, and gives them a reconstruction from the nearest old particles, + which is the wrong state for fluid that has just entered. The inlet + then carries a smear of whatever was upstream a step ago and hands it + downstream (#783). + + A particle entered if its back-trace leaves the domain and the flow + comes in where it left. The trace is :math:`x - \hat{u}\,s` with + :math:`s = |u|\,\Delta t` for a particle that moved (the test the + integration-point flavour applies to a restored departure point, + #745) and :math:`s = 2\,r_{\rm cell}` for one created this step, the + last ``created`` rows in storage, whose position inside the cell says + nothing about when it entered; :math:`r_{\rm cell}` is the RMS + vertex-to-centroid distance, so :math:`2r` is the cell's own extent + (0.9 to 1.3 of the edge on triangles). The trace of a particle beside + a wall can leave through the wall, at a corner, along a curved wall, + or where the flow separates, so the boundary velocity decides: at the + point where the trace left, the flow must cross the boundary inward + at more than 30 degrees. A no-slip wall has no velocity there and a + free-slip wall only a tangential one. Collective: the velocities and + the datum are read on every rank. + + TODO(DESIGN): a trace that wraps through a periodic seam also reads as + entered, as it does in the integration-point rule (#745). + """ + if self._inflow_value is None: + return + swarm = self.swarm + dim = self.mesh.dim + dt = self._nondim_timestep(dt) + X = np.asarray(swarm._particle_coordinates.data).reshape(-1, dim) + U = np.asarray(_to_nondim_ndarray(uw.function.evaluate(self.V_fn, X))).reshape(X.shape[0], dim) + speed = np.linalg.norm(U, axis=1) + moving = speed > 0.0 + direction = np.zeros_like(U) + direction[moving] = U[moving] / speed[moving, None] + distance = speed * dt + if created > 0: + cells = np.asarray(swarm._owning_cells())[-created:] + distance[-created:] = np.maximum( + distance[-created:], 2.0 * np.asarray(self.mesh._cell_radii)[cells]) + departure = X - direction * distance[:, None] + restored = np.asarray(self.mesh.return_coords_to_bounds(departure.copy())).reshape(departure.shape) + outward = departure - restored + left = moving & np.any(outward != 0.0, axis=1) + # The boundary velocity where each trace left; read on every rank. + U_boundary = np.asarray(_to_nondim_ndarray( + uw.function.evaluate(self.V_fn, restored[left]))).reshape(-1, dim) + crossing = np.einsum("ij,ij->i", U_boundary, outward[left]) + steep = crossing < -0.5 * np.linalg.norm(U_boundary, axis=1) * np.linalg.norm(outward[left], axis=1) + entered = np.zeros(X.shape[0], dtype=bool) + entered[np.nonzero(left)[0][steep]] = True + n_entered = int(entered.sum()) + if uw.mpi.size > 1: + n_entered = uw.mpi.comm.allreduce(n_entered, op=uw.MPI.SUM) + if n_entered == 0: + return + for slot in self.psi_star: + self._write_inflow(slot, X, entered) + def update_post_solve( self, dt: float, @@ -4192,6 +4288,10 @@ def update_post_solve( psi_star_0 = self.psi_star[0] coords = np.asarray(self.swarm._particle_coordinates.data) updated = {} + # TODO(BUG): the evaluated psi_fn is written into non-dimensional + # storage without reduction through _psi_units (see _write_inflow); a + # psi_fn carrying units lands at its physical magnitude. Same in + # initialise_history and in Lagrangian_Swarm. for i in range(psi_star_0.shape[0]): for j in range(psi_star_0.shape[1]): ij = psi_star_0._data_layout(i, j) @@ -4210,6 +4310,8 @@ def update_post_solve( delta_t=dt, restore_points_to_domain_func=self.mesh.return_coords_to_bounds, ) + created, _ = self.swarm.repopulate(**self._population_control) + self._apply_inflow_value(dt, created) if self._n_solves_completed < self.order: self._n_solves_completed += 1 @@ -5120,17 +5222,6 @@ def _write_components(self, var, expr, coords, evaluate=None, **kwargs): _to_nondim_ndarray(vals, units=self._psi_units) ).reshape(-1) - def _write_inflow(self, var, coords, rows): - """Overwrite ``rows`` of ``var`` with :attr:`inflow_value` evaluated at - ``coords`` (the restored positions of ALL the points, so that the - read, which is collective when the value holds a field, is made on - every rank; only ``rows`` are written).""" - expr = sympy.Matrix(self._inflow_value) - for column, (i, j) in enumerate(self._components): - vals = uw.function.evaluate(expr[i, j], coords) - var.data[rows, column] = np.asarray( - _to_nondim_ndarray(vals, units=self._psi_units)).reshape(-1)[rows] - def _segment_dt(self, j, dt): """Length of segment ``j`` (0 = the current step).""" if j == 0: @@ -5426,7 +5517,7 @@ def _geometry_stamp(self): """The mesh geometry the launch set was built for: the vertex count and the coordinate sum (a moved or re-meshed mesh changes one of them; adding a variable, which rebuilds the DM, changes neither).""" - coords = np.asarray(self.mesh.data) + coords = np.asarray(self.mesh.X.coords) return (coords.shape, float(coords.sum())) def _cell_measures(self): diff --git a/tests/test_0068_swarm_repopulation.py b/tests/test_0068_swarm_repopulation.py index 231be9f4..df799ce7 100644 --- a/tests/test_0068_swarm_repopulation.py +++ b/tests/test_0068_swarm_repopulation.py @@ -108,3 +108,32 @@ def test_population_control_keeps_cells_filled_under_rotation(): assert T._cell_projector.n_empty == 0 # The refilled corners carry the RBF reconstruction of x: bounded, no P2 blow-up assert np.abs(np.asarray(T.data[:, 0])).max() < 1.5 + + +def test_a_refill_after_a_removal_reads_the_right_neighbours(): + """PETSc removes a point by copying the LAST point into its slot, so a + removal reorders the surviving particles in storage. The reconstruction + for a particle created in the same call must be built on that storage + order, or its neighbours' values are read from unrelated rows (#784). + A linear field is reproduced exactly by the linear RBF, so any mismatch + shows as an error of order one.""" + mesh = uw.meshing.UnstructuredSimplexBox(minCoords=(0.0, 0.0), maxCoords=(1.0, 1.0), cellSize=0.1, qdegree=2) + swarm = uw.swarm.Swarm(mesh) + M = uw.swarm.SwarmVariable("M", swarm, 1) + swarm.populate(fill_param=3) + lattice = int(_census(swarm).max()) + # The band x < 0.2 is shifted onto [0.2, 0.4): the shifted band is empty, + # its destination holds twice the lattice, and one call must both thin + # the one and refill the other. + X = np.array(swarm._particle_coordinates.data) + X[X[:, 0] < 0.2, 0] += 0.2 + with uw.synchronised_array_update(): + swarm._particle_coordinates.data[...] = X + swarm.migrate() + X = np.asarray(swarm._particle_coordinates.data) + with uw.synchronised_array_update(): + M.data[:, 0] = X[:, 0] + 2.0 * X[:, 1] + added, removed = swarm.repopulate(max_per_cell=lattice, order=1) + assert added > 0 and removed > 0 + X = np.asarray(swarm._particle_coordinates.data) + assert np.abs(np.asarray(M.data[:, 0]) - (X[:, 0] + 2.0 * X[:, 1])).max() < 1.0e-8 diff --git a/tests/test_0075_lagrangian_history_inflow_value.py b/tests/test_0075_lagrangian_history_inflow_value.py new file mode 100644 index 00000000..c191b639 --- /dev/null +++ b/tests/test_0075_lagrangian_history_inflow_value.py @@ -0,0 +1,73 @@ +"""The particle stress history carries the inflow datum into an open channel. + +Poiseuille flow of a Maxwell fluid through a channel open at both ends, with +the stress history on the solver's own swarm (``stress_transport = +"lagrangian"``). The particles move downstream, the inlet cells empty and the +population control refills them; what the refilled particles carry is the +question. With no objective rate the fully developed stress is +:math:`\\sigma = 2\\eta\\dot{\\varepsilon}`, so the entering fluid's state is +known exactly and is what ``inflow_value`` says. A refill from the nearest +old particles instead smears the stress across the inlet cell and hands the +smear downstream (#783). +""" + +import numpy as np +import pytest +import sympy + +import underworld3 as uw + +pytestmark = [pytest.mark.level_2, pytest.mark.tier_b] + +ETA, MU = 1.0, 1.0 +LENGTH, HEIGHT, PEAK_SPEED = 4.0, 1.0, 1.0 +CELL = HEIGHT / 8 + + +def _channel(steps, dt, set_inflow): + mesh = uw.meshing.UnstructuredSimplexBox( + minCoords=(0.0, -HEIGHT / 2), maxCoords=(LENGTH, HEIGHT / 2), cellSize=CELL) + x, y = mesh.X + v = uw.discretisation.MeshVariable("U_inflow", mesh, mesh.dim, degree=2) + p = uw.discretisation.MeshVariable("P_inflow", mesh, 1, degree=1) + stokes = uw.systems.Stokes(mesh, velocityField=v, pressureField=p) + stokes.stress_transport = "lagrangian" + stokes.constitutive_model = uw.constitutive_models.ViscoElasticPlasticFlowModel( + stokes.Unknowns, order=1) + stokes.constitutive_model.Parameters.shear_viscosity_0 = ETA + stokes.constitutive_model.Parameters.shear_modulus = MU + stokes.constitutive_model.Parameters.dt_elastic = dt + + u_in = PEAK_SPEED * (1.0 - (2.0 * y / HEIGHT) ** 2) + stokes.add_dirichlet_bc((u_in, 0.0), "Left") + stokes.add_dirichlet_bc((sympy.oo, 0.0), "Right") + stokes.add_dirichlet_bc((0.0, 0.0), "Top") + stokes.add_dirichlet_bc((0.0, 0.0), "Bottom") + stokes.tolerance = 1.0e-6 + stokes.petsc_options["snes_type"] = "newtonls" + stokes.petsc_options["ksp_type"] = "fgmres" + + shear_rate = u_in.diff(y) + if set_inflow: + stokes.DFDt.inflow_value = sympy.Matrix([[0.0, ETA * shear_rate], [ETA * shear_rate, 0.0]]) + + for _ in range(steps): + stokes.solve(timestep=dt, zero_init_guess=False) + + # The inlet cell column, away from the walls: what the refilled particles carry. + probes = np.column_stack([np.full(7, 0.25 * CELL), np.linspace(-0.4, 0.4, 7)]) + carried = np.asarray(uw.function.evaluate(stokes.DFDt.psi_star[0].sym[0, 1], probes)).reshape(-1) + exact = np.asarray(uw.function.evaluate(ETA * shear_rate, probes)).reshape(-1) + return float(np.max(np.abs(carried - exact))) + + +def test_particle_history_applies_the_inflow_value(): + """After more than one transit every inlet particle has been refilled at + least once. With the inflow datum the carried shear stress matches + :math:`\\eta\\dot{\\gamma}` to the residual of the particles still loading + from the cold start; without it the inlet holds the neighbour smear.""" + assert uw.systems.ddt.Lagrangian.applies_inflow_value + error = _channel(steps=50, dt=0.1, set_inflow=True) + assert error < 0.01, error # measured 2.1e-3 (2026-09-23) + smear = _channel(steps=50, dt=0.1, set_inflow=False) + assert smear > 0.1, smear # measured 0.17