Skip to content

[BUG] - AMR fails to migrate swarm to new mesh #135

Description

@bknight1

The code sometimes produces an index error, most likely due to particles being lost between meshes.

Here's the error:

---------------------------------------------------------------------------
IndexError                                Traceback (most recent call last)
Cell In[12], line 1
----> 1 swarm._force_migration_after_mesh_change()

File ~/Documents/Research/GitHub/underworld3/src/underworld3/swarm.py:3301, in Swarm._force_migration_after_mesh_change(self)
   3297         disabled_arrays.append(var._array_cache)
   3299 try:
   3300     # Perform standard migration
-> 3301     self.migrate(remove_sent_points=True, delete_lost_points=True)
   3302 finally:
   3303     # Re-enable variable array callbacks
   3304     for array_cache in disabled_arrays:

File ~/Documents/Research/GitHub/underworld3/src/underworld3/timing.py:310, in routine_timer_decorator.<locals>.timed(*args, **kwargs)
    308 event.begin()
    309 try:
--> 310     result = routine(*args, **kwargs)
    311     return result
    312 finally:

File ~/Documents/Research/GitHub/underworld3/src/underworld3/mpi.py:257, in collective_operation.<locals>.wrapper(*args, **kwargs)
    239         error_msg = (
    240             f"\n{'='*70}\n"
    241             f"COLLECTIVE OPERATION DEADLOCK DETECTED\n"
   (...)    253             f"{'='*70}\n"
    254         )
    255         raise CollectiveOperationError(error_msg)
--> 257 return func(*args, **kwargs)

File ~/Documents/Research/GitHub/underworld3/src/underworld3/swarm.py:3180, in Swarm.migrate(self, remove_sent_points, delete_lost_points, max_its)
   3177 swarm_coord_array = (self.dm.getField("DMSwarmPIC_coor").reshape((-1, self.dim))).copy()
   3178 self.dm.restoreField("DMSwarmPIC_coor")
-> 3180 in_or_not = self.mesh.points_in_domain(
   3181     swarm_coord_array,
   3182 )
   3184 num_points_in_domain = np.count_nonzero(in_or_not == True)
   3185 num_points_not_in_domain = np.count_nonzero(in_or_not == False)

File ~/Documents/Research/GitHub/underworld3/src/underworld3/discretisation/discretisation_mesh.py:2648, in Mesh.points_in_domain(self, points, strict_validation)
   2644 near_boundary = numpy.where(dist2 < 2 * max_radius**2)[0]
   2645 near_boundary_points = model_points[near_boundary]
   2647 in_or_not[near_boundary] = (
-> 2648     self._get_closest_local_cells_internal(near_boundary_points) != -1
   2649 )
   2651 if strict_validation:
   2652     chosen_ones = numpy.where(in_or_not == True)[0]

File ~/Documents/Research/GitHub/underworld3/src/underworld3/discretisation/discretisation_mesh.py:2748, in Mesh._get_closest_local_cells_internal(self, coords)
   2745 cells = self._indexMap[closest_points]
   2746 cStart, cEnd = self.dm.getHeightStratum(0)
-> 2748 inside = self._test_if_points_in_cells_internal(coords, cells)
   2749 cells[~inside] = -1
   2750 lost_points = np.where(inside == False)[0]

File ~/Documents/Research/GitHub/underworld3/src/underworld3/discretisation/discretisation_mesh.py:2495, in Mesh._test_if_points_in_cells_internal(self, points, cells)
   2492 insiders = numpy.ndarray(shape=(cells.shape[0], num_cell_faces), dtype=bool)
   2494 for f in range(num_cell_faces):
-> 2495     control_points_o = self.faces_outer_control_points[f, cells]
   2496     control_points_i = self.faces_inner_control_points[f, cells]
   2497     inside = (
   2498         ((control_points_o - points) ** 2).sum(axis=1)
   2499         - ((control_points_i - points) ** 2).sum(axis=1)
   2500     ) > 0

IndexError: index 147 is out of bounds for axis 1 with size 128

Code to reproduce:


import underworld3 as uw
from underworld3.systems import Stokes
import numpy as np
import sympy
import os


res = 0.2
mesh = uw.meshing.UnstructuredSimplexBox(
    minCoords=(-1.0, 0.0),
    maxCoords=(1.0, 1.0),
    cellSize= res,
    regular=False,
    qdegree=3,
)

swarm = uw.swarm.Swarm(mesh=mesh)
swarm.populate(fill_param=4)


print(mesh.X.coords[:,0].min(), mesh.X.coords[:,0].max(), mesh.X.coords[:,1].min(), mesh.X.coords[:,1].max())


x, h_min, h_max, sigma = sympy.symbols('x h_min h_max sigma', positive=True)
h_expr = h_min + (h_max - h_min)*(1 - sympy.exp(-(x/sigma)**2))


h_subbed = h_expr.subs({h_min:0.05, h_max:0.5, sigma:0.25}) ### fails
# h_subbed = h_expr.subs({h_min:0.1, h_max:0.5, sigma:0.25}) ### works
h_vals = uw.function.evaluate(h_subbed, mesh.X.coords)


metric = uw.adaptivity.create_metric(mesh, h_vals)
mesh.adapt(metric)


print(mesh.X.coords[:,0].min(), mesh.X.coords[:,0].max(), mesh.X.coords[:,1].min(), mesh.X.coords[:,1].max())


swarm._force_migration_after_mesh_change()

Metadata

Metadata

Assignees

No one assigned

    Labels

    No labels
    No labels

    Type

    No type
    No fields configured for issues without a type.

    Projects

    No projects

    Milestone

    No milestone

    Relationships

    None yet

    Development

    No branches or pull requests

    Issue actions