diff --git a/docs/developer/subsystems/integration-point-variables.md b/docs/developer/subsystems/integration-point-variables.md index d8ee4cfc..b8edddcb 100644 --- a/docs/developer/subsystems/integration-point-variables.md +++ b/docs/developer/subsystems/integration-point-variables.md @@ -331,6 +331,117 @@ size are `nmin` and `patch_nnn` on `CellPolynomialProjector.fit`. M = uw.swarm.SwarmVariable("M", swarm, 1, proxy_location="cells", proxy_degree=2) ``` +### The swarm step at the mid time + +`swarm.advection(V_fn, dt, order=2, midtime_velocity=True)` evaluates the +RK2 mid-point velocity at the mid time, $\tfrac32 v^n - \tfrac12 v^{n-1}$, +from a `CharacteristicTrace` the swarm owns (the previous velocity is cached +by evaluation at the nodes at the end of each call), or from a solver's +shared trace passed as `characteristics=`. On a rotation whose rate ramps +linearly, ten steps of the frozen-velocity step miss 0.05 rad and the +mid-time step under 0.008 (`tests/test_0069_swarm_midtime_velocity.py`). A +steady flow is unchanged; the option is off by default and ignored when the +step is substepped. + +### Repopulation: keeping every cell fit-able + +A flow that empties cells starves the fit, and the two particle read-back +schemes both diverged on emptied corner cells before repopulation existed. +`Swarm.repopulate()` takes the per-cell census (owning cells from the strict +locator) and refills a starved cell from its own lattice, the points +`populate` uses, choosing the lattice points farthest from the particles +present. A new particle takes, for every swarm variable, the bounded Shepard +reconstruction from its nearest neighbours (`order=1` for the linear-exact +reconstruction; a starved cell is where neighbours are far, and the linear +tail extrapolated to values of 100 on a field bounded by 1), or a supplied +value (`values={var: constant or callable}`, an inflow datum for instance). +`swarm.population_control = dict(...)` makes every `advection()` end with a +repopulation, which is what the cells proxy wants: the refill runs before +the next fit. + +```python +swarm.population_control = dict() # refill to the populate() density +swarm.population_control = dict(min_per_cell=8, values={T: 0.0}) +``` + +Count is not the whole criterion. Particles the advection clamps back onto a +wall (`mesh.return_coords_to_bounds`) slide along it as a line, and the P2 +fit of a collinear set is singular whatever its count (measured: condition +number 1e300 at 92 particles in a wall cell, garbage that grew by 1e12 in +ten steps through the read-back). The fit therefore routes a cell whose Gram +matrix has condition number above `cond_max` (1e6) to the patch fit, and a +patch that is itself flat keeps only its mean. With that guard and +population control the untapered rotating box, where every wall has an +inflow and an outflow segment, runs to the same answer whether exiting +particles are clamped or deleted (`mesh.return_coords_to_bounds = None`, +the right setting for a true outflow, which also keeps the particle count +from growing). + +Measured on the rotating Gaussian (h = 0.1, C = 0.25, 10 particles per cell, +PIC, one revolution): population control takes the L2 error from 1.6e-2 to +8.8e-3, level with the integration-point history at 9.1e-3, because no cell +is ever left to the linear patch fit. A cap (`max_per_cell`) thins over-full +cells by removing the particles closest to a neighbour; measured it costs +accuracy (6.8e-2) and is off by default. + +### Viscoelastic stress history on particles + +The stress history of a viscoelastic Stokes solve is state: the stress at +the old time cannot be rebuilt from the present velocity gradient and the +rheology, and Crank-Nicolson keeps the elastic response undamped. Carried on +particles it is the Ellipsis / Underworld PIC-LIP arrangement: the particles +carry the stress along the flow, the mesh reads it at the integration points +through the cells proxy, and after each solve the new stress is evaluated at +the particles. That last step is a local ODE (the Maxwell update), so there +is no projection back to the mesh and no null space; the particle scheme's +one weakness, the re-projection of a diffused field, does not arise. + +```python +swarm = uw.swarm.Swarm(mesh) +DFDt = uw.systems.ddt.Lagrangian_Swarm( + swarm=swarm, psi_fn=sympy.Matrix.zeros(2, 2), vtype=uw.VarType.SYM_TENSOR, + degree=1, continuous=False, order=2, step_averaging=1, proxy_location="cells") +swarm.populate(fill_param=3) +swarm.population_control = dict() +stokes = uw.systems.Stokes(mesh, velocityField=v, pressureField=p, DFDt=DFDt) +stokes.constitutive_model = uw.constitutive_models.ViscoElasticPlasticFlowModel(stokes.Unknowns, order=2) +... +swarm.advection(v.sym, dt, order=2) # then +stokes.solve(timestep=dt) +``` + +The constitutive model reads the history through `psi_star[i].sym` and the +order bookkeeping only, so the swarm history slots in symbolically; the +solver assigns the stress expression to it, takes the viscoelastic order +from a supplied history, and leaves the nodal projection and shift to the +nodal history. The swarm manager evaluates every component of the new +stress at the particles before it shifts its chain, because the stress +expression reads the history it is about to overwrite. `step_averaging=1` +is required (the default 2 half-relaxes the stored stress). The ETD +integrator is not available on the swarm history. + +Maxwell shear box (`tests/test_0070_ve_stress_history_on_particles.py`): +order 1 within 5% of the analytic curve after 20 steps at dt = 0.1 t_r, +order 2 within 1%, and the particle and nodal histories agree to 0.2% of +the final stress. Uniform shear has a uniform stress, so this validates the +plumbing and the time integration. + +**Open (2026-09-09): a localised stress patch under shear.** With a +Gaussian patch in sigma_xy on the same shear box +(`~/+Simulations/integration_point_proxy/scripts/ve_stress_patch.py`), the +nodal and the particle histories each converge cleanly in h (4x per +doubling at the finest step) and in dt (first order, the same rate), but +to answers 1.2e-2 apart in L2 (peak 0.738 against 0.715), independent of +resolution, time step, particle density, proxy degree, read-back (PIC, +FLIP, an explicit P1 projection), mid-time velocity, and box width. Every +component agrees in isolation: with the particles held fixed the two +answers coincide (1.2e-3); transport alone against the exact sheared patch +puts the particles at 2e-5 and a P2 nodal history at 4e-5 (the P1 nodal +history is first order, 4e-3 at h/32), yet the coupled nodal answer does +not move when its history goes from P1 to P2; the particle proxy is +continuous across cells to 1e-8. Which limit is right is undecided and +needs a manufactured solution or an equation-level audit of both paths. + ### Why a least-squares fit and not a conservative transfer The conservative particle-to-mesh transfer solves the rule mass matrix diff --git a/src/underworld3/cython/petsc_generic_snes_solvers.pyx b/src/underworld3/cython/petsc_generic_snes_solvers.pyx index aa22e699..f7b1b195 100644 --- a/src/underworld3/cython/petsc_generic_snes_solvers.pyx +++ b/src/underworld3/cython/petsc_generic_snes_solvers.pyx @@ -2745,6 +2745,11 @@ class SolverBaseClass(uw_object): @constitutive_model.setter def constitutive_model(self, model_or_class): + # A stress history supplied by the user fixes the viscoelastic order + # (the solver's own _order is set only when it builds the history). + if self.Unknowns.DFDt is not None and self._order == 0: + self._order = getattr(self.Unknowns.DFDt, "order", 0) or 0 + ### checking if it's an instance - it will need to be reset if isinstance(model_or_class, uw.constitutive_models.Constitutive_Model): self._constitutive_model = model_or_class diff --git a/src/underworld3/swarm.py b/src/underworld3/swarm.py index ca6a93fb..7a0ff74c 100644 --- a/src/underworld3/swarm.py +++ b/src/underworld3/swarm.py @@ -2995,6 +2995,10 @@ def __init__(self, mesh, recycle_rate=0, verbose=False, clip_to_mesh=True): # a Lagrangian history registers its first sampling here so it sees # the field at the launch positions, not at the landing ones. self._pre_advection_hooks = [] + # Population control: a dict of repopulate() keyword arguments (or + # None). When set, advection() ends with repopulate(**population_control) + # so no cell is left starved before the next fit of a cells proxy. + self.population_control = None self._index = None # Particle -> proxy-node transfer operators, keyed by geometry and # stencil and shared by every proxied variable of this swarm. Entries @@ -4952,6 +4956,188 @@ def _data_layout(self, i, j=None): if self.vtype == uw.VarType.MATRIX: return i + j * self.shape[0] + @timing.routine_timer_decorator + @uw.collective_operation + def repopulate( + self, + min_per_cell=None, + max_per_cell=None, + values=None, + nnn=None, + order=0, + verbose=False, + ): + """Add particles to cells that hold too few, remove from cells that hold + too many, so every cell can support a well-posed fit of its particles. + + The trigger is the per-cell census (owning cells from the strict + locator). A starved cell is filled from its own lattice, the points + ``populate`` uses (degree ``fill_param``, cell interior), choosing the + lattice points farthest from the particles already present. A new + particle takes, for every variable, the RBF reconstruction from the + nearest existing particles at its position: bounded Shepard weights by + default (``order=0``), since a starved cell is where the neighbours + are far and a linear-exact tail extrapolates (measured: values of 100 + on a field bounded by 1 in the emptied corners of a rotating box); + ``order=1`` gives the linear-exact reconstruction. ``values`` overrides + a variable with a callable ``f(coords) -> (n, components)`` or a + constant, an inflow datum for instance. A cell above + ``max_per_cell`` loses its most redundant particles, those closest to + a neighbour in the same cell. + + Rank-local placement (a cell is filled by the rank that owns it), but + collective: every rank must call it, the domain test reduces. + + Parameters + ---------- + min_per_cell : int, optional + Particles a cell must hold; default the lattice count of + ``fill_param`` (the density ``populate`` gave). + max_per_cell : int, optional + Cap above which particles are removed; default no removal. + values : dict, optional + ``{variable or name: callable or constant}`` for new particles. + nnn : int, optional + Neighbours in the RBF reconstruction (default ``2 (dim + 1)``). + order : {0, 1}, optional + RBF reconstruction order for new particles: 0 bounded (default), + 1 linear-exact. + + Returns + ------- + (added, removed) : the counts on this rank. + """ + mesh = self.mesh + dim = self.cdim + fill = getattr(self, "fill_param", None) or 1 + lattice = np.asarray(mesh._get_coords_for_basis(fill, continuous=False)) + c0, c1 = mesh.dm.getHeightStratum(0) + ncells = c1 - c0 + n_lat = lattice.shape[0] // max(ncells, 1) + if min_per_cell is None: + min_per_cell = n_lat + if max_per_cell is not None: + min_per_cell = min(min_per_cell, max_per_cell) # a cap below the lattice count wins + + # Every rank must reach the (collective) domain test before any + # rank-local branch; the census itself is rank-local. + lat_owned = np.asarray(mesh.points_in_domain(lattice, strict_validation=True), dtype=bool) + self._flush_pending_petsc_sync() + X = np.array(self._particle_coordinates.data, copy=True) if self.local_size > 0 \ + else np.zeros((0, dim)) + cells = np.asarray(mesh._robust_owning_cells(X), dtype=np.int64) if X.shape[0] else np.zeros(0, np.int64) + npc = np.bincount(cells[cells >= 0], minlength=ncells) + lat_cells = np.asarray(mesh._robust_owning_cells(lattice), dtype=np.int64) + owned = np.zeros(ncells, dtype=bool) + owned[lat_cells[lat_owned & (lat_cells >= 0)]] = True + + added = removed = 0 + + # ---- removal: the most redundant particles of over-full cells ---------- + if max_per_cell is not None and X.shape[0] > 1: + drop = [] + over = np.nonzero(owned & (npc > max_per_cell))[0] + if over.shape[0] > 0: + # nearest-neighbour distance of every particle, one kd-tree query + d2, _ = uw.kdtree.KDTree(X).query(X, k=2) + nearest_all = np.sqrt(np.asarray(d2).reshape(X.shape[0], -1)[:, 1]) + for c in over: + idx = np.nonzero(cells == c)[0] + surplus = int(npc[c] - max_per_cell) + drop.extend(idx[np.argsort(nearest_all[idx])[:surplus]].tolist()) + if drop: + for index in sorted(drop, reverse=True): + self.dm.removePointAtIndex(int(index)) + 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() + + # ---- addition: starved cells, lattice points farthest from particles - + need = np.where(owned, np.maximum(min_per_cell - npc, 0), 0) + new_coords = [] + if need.sum() > 0: + cand_ok = lat_owned & (lat_cells >= 0) & (need[np.maximum(lat_cells, 0)] > 0) + cand = lattice[cand_ok] + cand_cells = lat_cells[cand_ok] + if X.shape[0] > 0: + dist, _ = uw.kdtree.KDTree(X).query(cand, k=1, sqr_dists=False) + dist = np.asarray(dist).reshape(-1) + else: + dist = np.zeros(cand.shape[0]) + sort_idx = np.lexsort((-dist, cand_cells)) # by cell, farthest first + cand, cand_cells, dist = cand[sort_idx], cand_cells[sort_idx], dist[sort_idx] + # rank within cell + start = np.searchsorted(cand_cells, np.arange(ncells), side="left") + rank_in_cell = np.arange(cand.shape[0]) - start[cand_cells] + take = rank_in_cell < need[cand_cells] + new_coords = cand[take] + + n_new = int(len(new_coords)) + if n_new > 0: + n_old = max(self.dm.getLocalSize(), 0) + nnn = nnn or 2 * (dim + 1) + nnn = min(nnn, max(n_old, 1)) + rbf_order = order if nnn >= dim + 2 else 0 + operator = None + if n_old > 0: + operator = uw.kdtree.KDTree(X).interpolation_matrix( + np.asarray(new_coords), nnn=nnn, p=2, order=rbf_order) + # raw values of every variable at the old particles, BEFORE the add + raw_old = {} + for name, var in self._vars.items(): + if var is self._particle_coordinates or var.clean_name in ( + "DMSwarmPIC_coor", "DMSwarm_rank", "DMSwarm_X0"): + continue + raw_old[name] = np.asarray(var.unpack_raw_data_from_petsc(squeeze=False)).reshape(n_old, -1) + + self.dm.finalizeFieldRegister() + # PETSc < 3.24 under-allocates by one on the first add to an empty + # swarm (the workaround populate() carries). + from petsc4py import PETSc + n_alloc = n_new + (1 if (n_old == 0 and PETSc.Sys.getVersion() < (3, 24, 0)) else 0) + self.dm.addNPoints(n_alloc) + coords = self.dm.getField("DMSwarmPIC_coor").reshape((-1, dim)) + coords[n_old:, :] = np.asarray(new_coords) + self.dm.restoreField("DMSwarmPIC_coor") + ranks = self.dm.getField("DMSwarm_rank") + ranks.reshape(-1)[n_old:] = uw.mpi.rank + self.dm.restoreField("DMSwarm_rank") + x0 = getattr(self, "_X0", None) + if x0 is not None: + f = self.dm.getField(x0.clean_name).reshape((-1, dim)) + f[n_old:, :] = np.asarray(new_coords) + self.dm.restoreField(x0.clean_name) + + values = values or {} + for name, var in self._vars.items(): + if name not in raw_old: + continue + spec = values.get(var, values.get(name, values.get(var.clean_name))) + ncomp = raw_old[name].shape[1] if n_old > 0 else var.num_components + if spec is not None: + vals = spec(np.asarray(new_coords)) if callable(spec) else spec + vals = np.broadcast_to(np.asarray(vals, dtype=float).reshape(n_new, -1) if np.ndim(vals) > 0 else vals, (n_new, ncomp)) + elif operator is not None: + vals = operator @ raw_old[name] + else: + vals = np.zeros((n_new, ncomp)) + f = self.dm.getField(var.clean_name).reshape((-1, ncomp)) + f[n_old:, :] = np.asarray(vals).reshape(n_new, ncomp) + self.dm.restoreField(var.clean_name) + added = n_new + self._invalidate_canonical_data() + + if added or removed: + self._population_generation += 1 + if verbose: + print(f"repopulate: rank {uw.mpi.rank} added {added}, removed {removed} " + f"(cells starved {int((need > 0).sum())})", flush=True) + return added, removed + + @timing.routine_timer_decorator def advection( self, @@ -4962,6 +5148,8 @@ def advection( restore_points_to_domain_func=None, evalf=False, step_limit=False, + midtime_velocity=False, + characteristics=None, ): r"""Advect the particle swarm through one timestep of a velocity field. @@ -5003,6 +5191,17 @@ def advection( Historical flag selecting RBF (``evalf``) velocity sampling; currently inert — velocity is always sampled with ``uw.function.global_evaluate``. Default ``False``. + midtime_velocity : bool, optional + Take the mid-point velocity of the RK2 step at the mid TIME, + :math:`\tfrac32 v^n - \tfrac12 v^{n-1}`, from a + :class:`~underworld3.systems.ddt.CharacteristicTrace` the swarm + owns (the previous velocity is cached by evaluation at the nodes + at the end of each call). Second order in an unsteady flow; a + steady flow is unchanged. Off by default. Ignored when the step + is substepped. + characteristics : CharacteristicTrace, optional + A solver's shared trace to take the mid-time velocity from + instead of a swarm-owned one (its levels are the solver's). step_limit : bool, optional If ``True``, split ``delta_t`` into substeps no larger than :meth:`estimate_dt` (roughly one element crossing per @@ -5088,6 +5287,19 @@ def advection( # Mesh.update_lvec(), and a migrate() firing there would reorder # particle rows between the coordinate array and the velocity array # captured from it. advection() performs its own migrate() at the end. + # Mid-time velocity for the RK2 mid-point stage: from the solver's + # shared trace when given, else a swarm-owned one (levels recorded at + # the end of each call). The levels, not the step cache, feed the + # expression, so a shared trace needs no step delimiting here. + trace = characteristics + owns_trace = False + if trace is None and midtime_velocity: + trace = self._characteristics_for(V_fn) + owns_trace = True + v_mid_matrix = V_fn_matrix + if trace is not None and substeps == 1: + v_mid_matrix = trace.midtime_expr() + self._deferred_migration_suspended = True # Wrap this whole thing in sub-stepping loop @@ -5098,7 +5310,8 @@ def advection( # Mid point algorithm (2nd order) if order == 2: - print(f"Advection (2nd): {self.local_size} - swarm points", flush=True) + if self.verbose: + print(f"Advection (2nd): {self.local_size} - swarm points", flush=True) # Use internal model-unit coordinates directly (no conversion needed) v_at_Vpts = np.zeros_like(self._particle_coordinates.data[...]) @@ -5128,7 +5341,7 @@ def advection( # (since the mid-points might have moved off-proc) # - v_at_Vpts[...] = uw.function.global_evaluate(V_fn_matrix, mid_pt_coords)[:, 0, :] + v_at_Vpts[...] = uw.function.global_evaluate(v_mid_matrix, mid_pt_coords)[:, 0, :] new_coords = X0.array[:, 0, :] + delta_t_model * v_at_Vpts / substeps @@ -5144,25 +5357,19 @@ def advection( # forward Euler (1st order) else: coords = self._particle_coordinates.data - print( - f"1. Advection (1st): {coords.shape} v {self.local_size} - swarm point shape", - flush=True, - ) + if self.verbose: + print(f"1. Advection (1st): {coords.shape} v {self.local_size} - swarm point shape", flush=True) v_at_Vpts = np.zeros_like(coords) v_at_Vpts[...] = uw.function.global_evaluate(V_fn_matrix, coords[...])[:, 0, :] - print( - f"2. Advection (1st): {coords.shape} v {self.local_size} - swarm point shape", - flush=True, - ) + if self.verbose: + print(f"2. Advection (1st): {coords.shape} v {self.local_size} - swarm point shape", flush=True) new_coords = coords[...] + delta_t_model * v_at_Vpts / substeps - print( - f"3. Advection (1st): {coords.shape} v {self.local_size} - swarm point shape", - flush=True, - ) + if self.verbose: + print(f"3. Advection (1st): {coords.shape} v {self.local_size} - swarm point shape", flush=True) if self.mesh.return_coords_to_bounds is not None: new_coords = self.mesh.return_coords_to_bounds(new_coords) @@ -5171,6 +5378,8 @@ def advection( ## End of substepping loop self._deferred_migration_suspended = False + if owns_trace: + trace.finish_step() # the velocity used this step becomes v^{n-1} # Re-route particles to their owning ranks and remove any that # have genuinely left the domain. Use the default max_its so that @@ -5182,8 +5391,22 @@ def advection( delete_lost_points=True, ) + if self.population_control is not None: + self.repopulate(**self.population_control) + return + def _characteristics_for(self, V_fn): + """The swarm-owned :class:`CharacteristicTrace` for ``V_fn`` (rebuilt + when the velocity expression changes).""" + from underworld3.systems.ddt import CharacteristicTrace, _matrix_of + + tr = getattr(self, "_characteristics", None) + if tr is None or not (tr.V_fn is V_fn or sympy.Matrix(_matrix_of(tr.V_fn)) == sympy.Matrix(_matrix_of(V_fn))): + tr = CharacteristicTrace(self.mesh, V_fn, midtime_velocity=True) + self._characteristics = tr + return tr + @timing.routine_timer_decorator def estimate_dt(self, V_fn): """ diff --git a/src/underworld3/systems/ddt.py b/src/underworld3/systems/ddt.py index 180e4346..23146df5 100644 --- a/src/underworld3/systems/ddt.py +++ b/src/underworld3/systems/ddt.py @@ -4117,8 +4117,12 @@ def update_pre_solve( dt: float, evalf: Optional[bool] = False, verbose: Optional[bool] = False, + **_ignored, ): - """Pre-solve: auto-initialise history on first call.""" + """Pre-solve: auto-initialise history on first call. Extra keyword + arguments (the nodal manager's ``store_result``, ``dt_physical``, + ``monotone_mode``) are accepted and ignored so a solver written for + the nodal history can drive this one.""" self._dt = dt if not self._history_initialised: @@ -4160,56 +4164,59 @@ def update_post_solve( dt: float, evalf: Optional[bool] = False, verbose: Optional[bool] = False, + **_ignored, ): - r"""Shift history chain and evaluate current :math:`\psi` on swarm.""" + r"""Evaluate the current :math:`\psi` at the particles, then shift the + history chain and store it in slot 0. + + The evaluation comes FIRST, every component of it: ``psi_fn`` may + read the history itself (a viscoelastic stress is + :math:`2\eta_{\rm eff}E_{\rm eff}(\sigma^*, \sigma^{**})`), and + writing a component of slot 0 marks its proxy stale, so an + evaluation after a partial write would read a history that is half + new (audit SWARM-06); a shift before the evaluation would hand the + stress expression the wrong levels. + """ self._dt = dt - # Record timestep history for variable-dt BDF - for i in range(self.order - 1, 0, -1): - self._dt_history[i] = self._dt_history[i - 1] - self._dt_history[0] = dt - - for h in range(self.order - 1): - i = self.order - (h + 1) - - # copy the information down the chain - if verbose: - print(f"Lagrange swarm order = {self.order}", flush=True) - print( - f"Mesh interpolant order = {self.psi_star[0]._meshVar.degree}", - flush=True, - ) - print(f"Lagrange swarm copying {i-1} to {i}", flush=True) - - self.psi_star[i].data[...] = self.psi_star[i - 1].data[...] - phi = 1 / self.step_averaging - psi_star_0 = self.psi_star[0] coords = np.asarray(self.swarm._particle_coordinates.data) if self.particle_update == "flip": # The proxy the mesh saw during this solve, at the particles: the # residual psi_p - proxy(x_p) is what the mesh never resolved. proxy_at_p = self._proxy_values_at_particles(psi_star_0, coords, evalf) - # Blend the freshly-evaluated psi into slot 0 component-by-component - # through the canonical (N, components) storage (audit SWARM-06). + updated = {} 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) - updated_psi = np.asarray( - uw.function.evaluate( - self.psi_fn[i, j], - coords, - evalf=evalf, - ) + if ij in updated: + continue # symmetric storage: one evaluation per slot + updated[ij] = np.asarray( + uw.function.evaluate(self.psi_fn[i, j], coords, evalf=evalf) ).reshape(-1) - if self.particle_update == "flip": - residual = np.asarray(psi_star_0.data[:, ij]) - proxy_at_p[:, ij] - psi_star_0.data[:, ij] = updated_psi + self.residual_retention * residual - else: - psi_star_0.data[:, ij] = ( - phi * updated_psi + (1 - phi) * psi_star_0.data[:, ij] - ) + + # Record timestep history for variable-dt BDF + for i in range(self.order - 1, 0, -1): + self._dt_history[i] = self._dt_history[i - 1] + self._dt_history[0] = dt + + for h in range(self.order - 1): + i = self.order - (h + 1) + if verbose: + print(f"Lagrange swarm order = {self.order}: copying slot {i-1} to {i}", flush=True) + self.psi_star[i].data[...] = self.psi_star[i - 1].data[...] + + # Store slot 0 component-by-component through the canonical + # (N, components) storage (audit SWARM-06). + for ij, updated_psi in updated.items(): + if self.particle_update == "flip": + residual = np.asarray(psi_star_0.data[:, ij]) - proxy_at_p[:, ij] + psi_star_0.data[:, ij] = updated_psi + self.residual_retention * residual + else: + psi_star_0.data[:, ij] = ( + phi * updated_psi + (1 - phi) * psi_star_0.data[:, ij] + ) if self._n_solves_completed < self.order: self._n_solves_completed += 1 diff --git a/src/underworld3/systems/solvers.py b/src/underworld3/systems/solvers.py index 24351f47..a9f59bbd 100644 --- a/src/underworld3/systems/solvers.py +++ b/src/underworld3/systems/solvers.py @@ -1727,37 +1727,42 @@ def solve( if uw.mpi.rank == 0 and verbose: print(f"Stokes solver - store stress and shift history", flush=True) - _advected_sigma_star = np.copy(self.DFDt.psi_star[0].array[...]) - - if getattr(self.DFDt, '_psi_star_use_multicomponent', False): - # Multi-component projection of flux → psi_star[0]. - # - # The DFDt's source-snapshot machinery (enabled once in - # _create_stress_history_ddt) intercepts psi_fn assignment - # to substitute psi_star[0] symbols with a frozen - # psi_snapshot variable, refreshed each step in - # update_pre_solve. So the projection's compiled source - # reads from psi_snapshot (not psi_star[0] itself) and is a - # true one-shot Galerkin projection — no implicit - # fixed-point iteration. - self.DFDt._psi_star_projection_solver.smoothing = 0.0 - self.DFDt._psi_star_projection_solver.solve(verbose=verbose) - # Fan flat result back to psi_star[0] tensor variable - for k, (i, j) in enumerate(self.DFDt._psi_star_indep_indices): - vals = self.DFDt._psi_star_flat_var.array[:, 0, k] - self.DFDt.psi_star[0].array[:, i, j] = vals - if i != j: - self.DFDt.psi_star[0].array[:, j, i] = vals - else: - self.DFDt._psi_star_projection_solver.uw_function = self.constitutive_model.flux - self.DFDt._psi_star_projection_solver.smoothing = 0.0 - self.DFDt._psi_star_projection_solver.solve(verbose=verbose) - - for i in range(self.DFDt.order - 1, 0, -1): - if i == 1: - self.DFDt.psi_star[i].array[...] = _advected_sigma_star + # A particle-carried history (Lagrangian_Swarm) evaluates the new + # stress at its particles and shifts its own chain in + # update_post_solve; the projection and shift below are the + # nodal semi-Lagrangian history's. + if isinstance(self.DFDt, SemiLagrangian_DDt): + _advected_sigma_star = np.copy(self.DFDt.psi_star[0].array[...]) + + if getattr(self.DFDt, '_psi_star_use_multicomponent', False): + # Multi-component projection of flux → psi_star[0]. + # + # The DFDt's source-snapshot machinery (enabled once in + # _create_stress_history_ddt) intercepts psi_fn assignment + # to substitute psi_star[0] symbols with a frozen + # psi_snapshot variable, refreshed each step in + # update_pre_solve. So the projection's compiled source + # reads from psi_snapshot (not psi_star[0] itself) and is a + # true one-shot Galerkin projection — no implicit + # fixed-point iteration. + self.DFDt._psi_star_projection_solver.smoothing = 0.0 + self.DFDt._psi_star_projection_solver.solve(verbose=verbose) + # Fan flat result back to psi_star[0] tensor variable + for k, (i, j) in enumerate(self.DFDt._psi_star_indep_indices): + vals = self.DFDt._psi_star_flat_var.array[:, 0, k] + self.DFDt.psi_star[0].array[:, i, j] = vals + if i != j: + self.DFDt.psi_star[0].array[:, j, i] = vals else: - self.DFDt.psi_star[i].array[...] = self.DFDt.psi_star[i - 1].array[...] + self.DFDt._psi_star_projection_solver.uw_function = self.constitutive_model.flux + self.DFDt._psi_star_projection_solver.smoothing = 0.0 + self.DFDt._psi_star_projection_solver.solve(verbose=verbose) + + for i in range(self.DFDt.order - 1, 0, -1): + if i == 1: + self.DFDt.psi_star[i].array[...] = _advected_sigma_star + else: + self.DFDt.psi_star[i].array[...] = self.DFDt.psi_star[i - 1].array[...] self.DFDt.update_post_solve(timestep, verbose=verbose, evalf=evalf) diff --git a/src/underworld3/utilities/cell_polynomial_projection.py b/src/underworld3/utilities/cell_polynomial_projection.py index 4fa46067..c07b04ad 100644 --- a/src/underworld3/utilities/cell_polynomial_projection.py +++ b/src/underworld3/utilities/cell_polynomial_projection.py @@ -118,7 +118,7 @@ def locate(self, coords): # -- the fit ------------------------------------------------------------ - def fit(self, coords, values, nmin=None, patch_nnn=None, old=None): + def fit(self, coords, values, nmin=None, patch_nnn=None, old=None, cond_max=1.0e6): """Fit every cell; returns nodal values shaped like ``meshVar.data``. Parameters @@ -131,6 +131,12 @@ def fit(self, coords, values, nmin=None, patch_nnn=None, old=None): old : current proxy values, shaped like ``meshVar.data``; after the first fit a cell with no particles keeps them (on the first fit, or without ``old``, it takes the linear patch fit). + cond_max : a cell whose Gram matrix has a condition number above this + is treated as thin (patch fit) however many particles it holds. + Count is not enough: particles clamped onto a wall by the + advection lie on a line, and the P2 fit of a line is singular + (measured: condition 1e300 at 92 particles, garbage that grew + by 1e12 in ten steps through the read-back). """ coords = np.asarray(coords, dtype=np.float64).reshape(-1, self.dim) values = np.asarray(values, dtype=np.float64).reshape(coords.shape[0], -1) @@ -149,6 +155,15 @@ def fit(self, coords, values, nmin=None, patch_nnn=None, old=None): U = np.zeros((self.ncells, self.Nb, nc)) nmin = nmin or self.Nb + 2 dense = npc >= nmin + self.n_ill_conditioned = 0 + if dense.any(): + ev = np.linalg.eigvalsh(G[dense]) + cond = ev[:, -1] / np.maximum(ev[:, 0], 1e-300) + ill = cond > cond_max + if ill.any(): + self.n_ill_conditioned = int(ill.sum()) + idx = np.nonzero(dense)[0][ill] + dense[idx] = False if dense.any(): ridge = 1e-10 * np.trace(G[dense], axis1=1, axis2=2)[:, None, None] / self.Nb U[dense] = np.linalg.solve(G[dense] + ridge * np.eye(self.Nb)[None], R[dense]) @@ -175,6 +190,13 @@ def fit(self, coords, values, nmin=None, patch_nnn=None, old=None): Rt = np.einsum("cpa,cpk->cak", A, psi[idx]) ridge = 1e-10 * np.trace(Gt, axis1=1, axis2=2)[:, None, None] / (self.dim + 1) + 1e-30 coef = np.linalg.solve(Gt + ridge * np.eye(self.dim + 1)[None], Rt) # (nthin, dim+1, nc) + # A patch whose particles are themselves (nearly) collinear cannot + # carry a gradient: keep only the constant term (the patch mean). + evt = np.linalg.eigvalsh(Gt) + flat = evt[:, -1] / np.maximum(evt[:, 0], 1e-300) > cond_max + if flat.any(): + coef[flat, 1:, :] = 0.0 + coef[flat, 0, :] = psi[idx[flat]].mean(axis=1) Adof = np.concatenate([np.ones((self.Nb, 1)), self.xi_dof], axis=1) # (Nb, dim+1) U[thin] = np.einsum("ba,cak->cbk", Adof, coef) diff --git a/tests/test_0067_integration_point_proxy.py b/tests/test_0067_integration_point_proxy.py index 9399f3db..89f750d7 100644 --- a/tests/test_0067_integration_point_proxy.py +++ b/tests/test_0067_integration_point_proxy.py @@ -299,13 +299,42 @@ def test_lagrangian_swarm_history_is_sampled_before_the_first_move(): swarm=swarm, psi_fn=T.sym, vtype=uw.VarType.SCALAR, degree=2, continuous=False, order=1, proxy_location="cells", ) + X0 = uw.swarm.SwarmVariable("X0", swarm, 2) # launch position, carried by the particle swarm.populate(fill_param=2) assert not lag._history_initialised - X_before = np.array(swarm._particle_coordinates.data, copy=True) + with uw.synchronised_array_update(): + X0.data[...] = np.asarray(swarm._particle_coordinates.data) swarm.advection(sympy.Matrix([[0.1, 0.0]]), 0.5, order=2) # every particle moves +0.05 in x - X_after = np.asarray(swarm._particle_coordinates.data) assert lag._history_initialised + # Particles may have changed rank: compare each against the launch + # position it carries, not against a rank-local array from before. + X_before = np.asarray(X0.data) + X_after = np.asarray(swarm._particle_coordinates.data) kept = np.abs(X_after[:, 0] - X_before[:, 0] - 0.05) < 1e-12 # particles not returned to bounds vals = np.asarray(lag.psi_star[0].data[:, 0]) assert np.allclose(vals[kept], X_before[kept, 0], atol=1e-10) # launch positions ... assert not np.allclose(vals[kept], X_after[kept, 0], atol=1e-3) # ... not landing positions + + +def test_cells_proxy_collinear_particles_do_not_blow_up(): + """Particles clamped onto a wall lie on a line; the P2 fit of a line is + singular. The fit routes such cells to the patch on the Gram condition + number, and a patch that is itself flat keeps only its mean.""" + mesh = uw.meshing.UnstructuredSimplexBox(cellSize=0.1, qdegree=2) + swarm = uw.swarm.Swarm(mesh) + M = uw.swarm.SwarmVariable("M", swarm, 1, proxy_location="cells", proxy_degree=2) + swarm.populate(fill_param=3) + X = np.array(swarm._particle_coordinates.data, copy=True) + # Pile every particle of the bottom row of cells onto the wall y = 1e-9 + bottom = X[:, 1] < 0.1 + X[bottom, 1] = 1e-9 + with uw.synchronised_array_update(): + M.data[:, 0] = 1.0 + X[:, 0] # values first: the move below migrates + with uw.synchronised_array_update(): + swarm._particle_coordinates.data[...] = X + swarm.migrate() + M._update_proxy_if_stale() + pr = M._cell_projector + assert pr.n_ill_conditioned > 0 + vals = np.asarray(M._meshVar.data[:, 0]) + assert np.isfinite(vals).all() and vals.min() > 0.5 and vals.max() < 2.5, (vals.min(), vals.max()) diff --git a/tests/test_0068_swarm_repopulation.py b/tests/test_0068_swarm_repopulation.py new file mode 100644 index 00000000..231be9f4 --- /dev/null +++ b/tests/test_0068_swarm_repopulation.py @@ -0,0 +1,110 @@ +"""Swarm.repopulate: particles added to starved cells and removed from over-full ones. + +The trigger is the per-cell census; a starved cell is filled from its own +lattice at the points farthest from the particles present; a new particle +takes the linear-exact RBF reconstruction of every variable from the nearest +existing particles (or a supplied value). ``swarm.population_control`` makes +advection() end with a repopulation, which is what keeps a cells proxy well +posed on a flow that empties cells. +""" + +import numpy as np +import pytest +import sympy + +import underworld3 as uw + +pytestmark = [pytest.mark.level_1, pytest.mark.tier_a] + + +def _census(swarm): + mesh = swarm.mesh + c0, c1 = mesh.dm.getHeightStratum(0) + X = np.asarray(swarm._particle_coordinates.data) + cells = np.asarray(mesh._robust_owning_cells(X)) + return np.bincount(cells[cells >= 0], minlength=c1 - c0) + + +def _strip_left(swarm, xmax): + """Test helper: delete every particle with x < xmax (no public removal API).""" + X = np.asarray(swarm._particle_coordinates.data) + for index in np.sort(np.nonzero(X[:, 0] < xmax)[0])[::-1]: + swarm.dm.removePointAtIndex(int(index)) + swarm._invalidate_canonical_data() + swarm._population_generation += 1 + + +def test_starved_cells_are_refilled_with_linear_exact_values(): + mesh = uw.meshing.UnstructuredSimplexBox(cellSize=0.1, qdegree=2) + x, y = mesh.X + swarm = uw.swarm.Swarm(mesh) + T = uw.swarm.SwarmVariable("T", swarm, 1, proxy_location="cells", proxy_degree=2) + V = uw.swarm.SwarmVariable("V", swarm, vtype=uw.VarType.VECTOR) + swarm.populate(fill_param=3) + n_full = swarm.local_size + before = _census(swarm) + X = np.asarray(swarm._particle_coordinates.data) + with uw.synchronised_array_update(): + T.data[:, 0] = 1.0 + 2.0 * X[:, 0] + 3.0 * X[:, 1] + V.data[:, 0] = X[:, 0] + V.data[:, 1] = -X[:, 1] + _strip_left(swarm, 0.35) + assert (_census(swarm) == 0).sum() > 0 # cells emptied + added, removed = swarm.repopulate(order=1) # linear-exact reconstruction + assert removed == 0 and added > 0 + after = _census(swarm) + # Every owned cell holds at least the lattice count again (10 for fill 3) + assert (after >= before.min()).all(), (after.min(), before.min()) + # Values: the linear fields are reconstructed exactly on the new particles + X = np.asarray(swarm._particle_coordinates.data) + assert np.allclose(np.asarray(T.data[:, 0]), 1.0 + 2.0 * X[:, 0] + 3.0 * X[:, 1], atol=1e-8) + assert np.allclose(np.asarray(V.data), np.c_[X[:, 0], -X[:, 1]], atol=1e-8) + # ... and the cells proxy reproduces the field through the weak form + err = uw.maths.Integral(mesh, (T.sym[0] - (1 + 2 * x + 3 * y)) ** 2).evaluate() + assert err < 1e-14, err + assert T._cell_projector.n_empty == 0 + # Idempotent: nothing to do on a healthy swarm + assert swarm.repopulate() == (0, 0) + assert swarm.local_size <= n_full + 0 # no more than the original lattice + + +def test_values_override_and_cap(): + mesh = uw.meshing.UnstructuredSimplexBox(cellSize=0.2, qdegree=2) + swarm = uw.swarm.Swarm(mesh) + M = uw.swarm.SwarmVariable("M", swarm, 1) + swarm.populate(fill_param=3) + with uw.synchronised_array_update(): + M.data[:, 0] = 1.0 + _strip_left(swarm, 0.5) + added, _ = swarm.repopulate(values={M: 7.0}) # default (Shepard) elsewhere + assert added > 0 + X = np.asarray(swarm._particle_coordinates.data) + vals = np.asarray(M.data[:, 0]) + assert np.allclose(vals[X[:, 0] < 0.5], 7.0) and np.allclose(vals[X[:, 0] >= 0.5], 1.0) + # A callable datum, and a cap that thins over-full cells + swarm.repopulate(min_per_cell=14, values={"M": lambda C: 2.0 * C[:, 1:2]}) + assert (_census(swarm) >= 14).all() + added, removed = swarm.repopulate(max_per_cell=8) + assert added == 0 and removed > 0 + assert (_census(swarm) <= 8).all() + + +def test_population_control_keeps_cells_filled_under_rotation(): + """Solid-body rotation of a box empties the corner cells; with population + control the cells proxy never sees an empty cell.""" + mesh = uw.meshing.UnstructuredSimplexBox(minCoords=(-1.0, -1.0), maxCoords=(1.0, 1.0), cellSize=0.2, qdegree=2) + x, y = mesh.X + V = sympy.Matrix([[-y, x]]) + swarm = uw.swarm.Swarm(mesh) + T = uw.swarm.SwarmVariable("T", swarm, 1, proxy_location="cells", proxy_degree=2) + swarm.populate(fill_param=3) + with uw.synchronised_array_update(): + T.data[:, 0] = np.asarray(swarm._particle_coordinates.data)[:, 0] + swarm.population_control = dict(min_per_cell=6) + for _ in range(4): + swarm.advection(V, 0.25, order=2) + assert (_census(swarm) >= 6).all() + T._update_proxy_if_stale() + 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 diff --git a/tests/test_0069_swarm_midtime_velocity.py b/tests/test_0069_swarm_midtime_velocity.py new file mode 100644 index 00000000..b170d31f --- /dev/null +++ b/tests/test_0069_swarm_midtime_velocity.py @@ -0,0 +1,65 @@ +"""Swarm advection with the mid-point velocity taken at the mid time. + +The RK2 step evaluates the mid-point velocity at the mid TIME, +1.5 v^n - 0.5 v^{n-1}, from a CharacteristicTrace (swarm-owned, or a +solver's shared one). On a rotation whose rate ramps linearly in time the +frozen-velocity step is first order in the rate; the mid-time step is exact +for the linear ramp after its first step. +""" + +import numpy as np +import pytest +import sympy + +import underworld3 as uw + +pytestmark = [pytest.mark.level_1, pytest.mark.tier_a] + + +def _angle_error(midtime, nsteps=10, dt=0.1): + mesh = uw.meshing.UnstructuredSimplexBox(minCoords=(-1.0, -1.0), maxCoords=(1.0, 1.0), cellSize=0.2, qdegree=2) + mesh.return_coords_to_bounds = None + x, y = mesh.X + c = uw.expression(r"c_{rate}", 1.0, "rotation rate") + V = c * sympy.Matrix([[-y, x]]) + swarm = uw.swarm.Swarm(mesh) + swarm.verbose = False + L = uw.swarm.SwarmVariable("L", swarm, 2) # launch position, carried by the particle + swarm.populate(fill_param=1) + with uw.synchronised_array_update(): + L.data[...] = np.asarray(swarm._particle_coordinates.data) + t = 0.0 + for _ in range(nsteps): + c.sym = 1.0 + t # the rate at t^n; the mid-time rate is 1 + t + dt/2 + swarm.advection(V, dt, order=2, midtime_velocity=midtime) + t += dt + X = np.asarray(swarm._particle_coordinates.data) + X0 = np.asarray(L.data) + inner = np.hypot(X0[:, 0], X0[:, 1]) < 0.6 + exact = nsteps * dt + 0.5 * (nsteps * dt) ** 2 # integral of (1 + t) + ang = np.arctan2(X[inner, 1], X[inner, 0]) - np.arctan2(X0[inner, 1], X0[inner, 0]) + ang = (ang + np.pi) % (2 * np.pi) - np.pi + return float(np.abs(ang - exact).mean()) + + +def test_midtime_velocity_makes_the_swarm_step_second_order_in_time(): + frozen = _angle_error(False) + midtime = _angle_error(True) + # Frozen rate: each step misses dt/2 of ramp over dt: 0.5 * N * dt^2 = 0.05 rad. + assert 0.03 < frozen < 0.07, frozen + # Mid-time rate: exact for the linear ramp except the first step (no v^{n-1}): 0.005 rad. + assert midtime < 0.008, midtime + assert midtime < frozen / 5 + + +def test_shared_trace_from_a_solver_is_accepted(): + mesh = uw.meshing.UnstructuredSimplexBox(cellSize=0.25, qdegree=2) + x, y = mesh.X + V = sympy.Matrix([[-y, x]]) + tr = uw.systems.ddt.CharacteristicTrace(mesh, V) + swarm = uw.swarm.Swarm(mesh) + swarm.populate(fill_param=1) + X0 = np.array(swarm._particle_coordinates.data, copy=True) + swarm.advection(V, 0.05, order=2, characteristics=tr) + assert not hasattr(swarm, "_characteristics") or swarm._characteristics is not tr + assert np.abs(np.asarray(swarm._particle_coordinates.data) - X0).max() > 1e-3 diff --git a/tests/test_0070_ve_stress_history_on_particles.py b/tests/test_0070_ve_stress_history_on_particles.py new file mode 100644 index 00000000..f5e6b88a --- /dev/null +++ b/tests/test_0070_ve_stress_history_on_particles.py @@ -0,0 +1,89 @@ +"""Viscoelastic stress history carried by particles. + +The Maxwell shear box (test_1051) with the stress history on a swarm: +``Lagrangian_Swarm`` with the cells proxy supplied as the Stokes solver's +``DFDt``. The constitutive model reads the history through the proxy's +symbol exactly as it reads the nodal one; the particles carry the stress +along the flow, the mesh integrates it at the integration points, and after +each solve the new stress is evaluated at the particles (the Maxwell update +is a local ODE, no projection back to the mesh). Uniform shear has a uniform +stress, so this checks the plumbing and the time integration against the +analytic Maxwell curve, and against the nodal history. +""" + +import numpy as np +import pytest +import sympy + +import underworld3 as uw + +pytestmark = [pytest.mark.level_2, pytest.mark.tier_a] + +ETA, MU, V0, H, W = 1.0, 1.0, 0.5, 1.0, 2.0 + + +def maxwell_stress_xy(t, gamma_dot): + return ETA * gamma_dot * (1.0 - np.exp(-t / (ETA / MU))) + + +def _run(history, order, n_steps, dt_over_tr, res=8): + dt = dt_over_tr * ETA / MU + gamma_dot = 2.0 * V0 / H + mesh = uw.meshing.StructuredQuadBox(elementRes=(2 * res, res), + minCoords=(-W / 2, -H / 2), maxCoords=(W / 2, H / 2)) + v = uw.discretisation.MeshVariable("U", mesh, mesh.dim, degree=2) + p = uw.discretisation.MeshVariable("P", mesh, 1, degree=1) + swarm = DFDt = None + if history == "particles": + swarm = uw.swarm.Swarm(mesh) + DFDt = uw.systems.ddt.Lagrangian_Swarm( + swarm=swarm, psi_fn=sympy.Matrix.zeros(2, 2), vtype=uw.VarType.SYM_TENSOR, + degree=1, continuous=False, order=order, step_averaging=1, proxy_location="cells", + ) + swarm.populate(fill_param=3) + swarm.population_control = dict() # open sides: inflow refill, outflow loss + stokes = uw.systems.Stokes(mesh, velocityField=v, pressureField=p, DFDt=DFDt) + stokes.constitutive_model = uw.constitutive_models.ViscoElasticPlasticFlowModel( + stokes.Unknowns, order=order) + stokes.constitutive_model.Parameters.shear_viscosity_0 = ETA + stokes.constitutive_model.Parameters.shear_modulus = MU + stokes.constitutive_model.Parameters.dt_elastic = dt + stokes.add_dirichlet_bc((V0, 0.0), "Top") + stokes.add_dirichlet_bc((-V0, 0.0), "Bottom") + stokes.add_dirichlet_bc((sympy.oo, 0.0), "Left") + stokes.add_dirichlet_bc((sympy.oo, 0.0), "Right") + stokes.tolerance = 1.0e-6 + stokes.petsc_options["snes_type"] = "newtonls" + stokes.petsc_options["ksp_type"] = "fgmres" + assert stokes.DFDt is DFDt if DFDt is not None else stokes.DFDt is not None + num, ana = [], [] + time = 0.0 + for _ in range(n_steps): + if swarm is not None: + swarm.advection(v.sym, dt, order=2) + stokes.solve(timestep=dt, zero_init_guess=False, evalf=False) + time += dt + num.append(float(np.asarray(uw.function.evaluate(stokes.tau.sym[0, 1], np.array([[0.0, 0.0]]))).flatten()[0])) + ana.append(maxwell_stress_xy(time, gamma_dot)) + return np.array(num), np.array(ana) + + +def test_particle_stress_history_tracks_maxwell_order1(): + num, ana = _run("particles", order=1, n_steps=20, dt_over_tr=0.1) + rel_final = abs(num[-1] - ana[-1]) / abs(ana[-1]) + assert rel_final < 0.05, rel_final + assert np.all(np.diff(num) > -1e-8) # monotone loading + + +def test_particle_and_nodal_histories_agree(): + """Uniform shear: the two histories carry the same stress, so the answers + coincide to the transport reconstruction error.""" + num_p, ana = _run("particles", order=1, n_steps=10, dt_over_tr=0.1) + num_n, _ = _run("nodal", order=1, n_steps=10, dt_over_tr=0.1) + assert np.abs(num_p - num_n).max() < 2e-3 * abs(ana[-1]), np.abs(num_p - num_n).max() + + +def test_particle_stress_history_order2(): + num, ana = _run("particles", order=2, n_steps=20, dt_over_tr=0.1) + rel_final = abs(num[-1] - ana[-1]) / abs(ana[-1]) + assert rel_final < 0.01, rel_final