# Using WhenEvent to stop N objects at a boundary

I am attempting to implement an event into my simulation which makes any moving particle in an electric field completely stop at the instant they pass some specified threshold.

I was trying to utilize DiscreteVariable, as presented in other threads to accomplish this:

Using WhenEvent to Bound Variables

but can't seem to work it out in my particular application.

The boundary I am trying to implement is seen here under event:

eqs = Table[{x1[j]''[t] == -(1/mass)*eforceX[x, y, z] /. {x ->
x1[j][t], y -> y1[j][t], z -> z1[j][t]},
y1[j]''[t] == -(1/mass)*eforceY[x, y, z] /. {x -> x1[j][t],
y -> y1[j][t], z -> z1[j][t]},
z1[j]''[t] == -(1/mass)*eforceZ[x, y, z] /. {x -> x1[j][t],
y -> y1[j][t], z -> z1[j][t]},
x1[j][0] == pos0[[j, 1, 1]],
y1[j][0] == pos0[[j, 2, 1]],
z1[j][0] == pos0[[j, 3, 1]],
x1[j]'[0] == vel0[[j, 1, 1]],
y1[j]'[0] == vel0[[j, 2, 1]],
z1[j]'[0] == vel0[[j, 3, 1]]},
{j, numbodies}];

vars = Flatten[Table[{x1[j][t], y1[j][t], z1[j][t]}, {j, numbodies}]];

event = Table[{WhenEvent[
z1[j][t] >= -partRad , {x1[j]'[t] -> 0, y1[j]'[t] -> 0,
z1[j]'[t] -> 0}]} /. {j -> i}, {i, numbodies}];

sol1 = NDSolve[{eqs, event}, Head /@ vars, {t, 0, tfin}]


How would one go about setting up an actual boundary condition, such that if one of the N particles passes the boundary z1[j][t] >= -partRadius, they completely stop?

The full program is shown here:

Clear["Global*"];
Needs["NDSolveFEM"]

q = 1.617733*10^-18;(*Net ion Charge*)

R= Import["https://www.dropbox.com/s/dds8rm3odg2m7gu/largeAp.obj?dl=\
1"];

RegionDimension[R];
M = BoundaryMeshRegion[MeshCoordinates[R], MeshCells[R, 2]];
RegionDimension[M];
Volume[M];

r = RegionDifference[
RegionDifference[
RegionDifference[Cuboid[{0, 0, -6.5}, {2, 2, 0.5}], M],
Cuboid[{0, 0, 0.4}, {2, 2, 0.5}]],
Cuboid[{0, 0, -6.5}, {2, 2, -6.4}]];
ToElementMesh[r]["Wireframe"];
pol = -1;

V0 = 1000;
sol = NDSolveValue[{Laplacian[V[x, y, z], {x, y, z}] == 0,
DirichletCondition[
V[x, y, z] == -pol*
V0/2, (0.4 <= z <= 0.5) && (0 <= y <= 2) && (0 <= x <= 2)],
DirichletCondition[
V[x, y, z] ==
pol*V0/2, (0.0071 <= z <= 0.0072) && (0 <= y <= 2) && (0 <= x <=
2)], DirichletCondition[
V[x, y, z] ==
0, (0 <= z <= 0.0070) && (0 <= y <= 2) && (0 <= x <= 2)],
DirichletCondition[
V[x, y, z] ==
0, (-6.5 <= z <= -6.4) && (0 <= y <= 2) && (0 <= x <= 2)]},
V, {x, y, z} \[Element] r];

electricField[x_, y_, z_] = -Grad[sol[x, y, z], {x, y, z}];

v = Show[VectorPlot3D[
electricField[x, y, z], {x, 0.5, 1}, {y, 0.5, 1}, {z, -0.5, 0.1},
PlotTheme -> "Detailed", ColorFunction -> "Rainbow",
PerformanceGoal -> "Quality", VectorScale -> 0.05,
VectorPoints -> 7, PlotLegends -> Automatic], M];

vecForce =
Show[VectorPlot3D[
q*electricField[x, y, z], {x, 0.5, 1}, {y, 0.5,
1}, {z, -0.25, -0.00001}, PlotTheme -> "Detailed",
ColorFunction -> "Rainbow", PerformanceGoal -> "Quality",
VectorScale -> 0.05, VectorPoints -> 7, PlotLegends -> Automatic],
M];

numbodies = 3;
mass = 6.52*10^-11;(*particle mass in kg/m^3*)

vel0 = Table[Partition[{0, 0, 0}, 1], numbodies];
(*vel0=Table[Partition[{RandomReal[{-0.0001,0.0001}],RandomReal[{-0.\
0001,0.0001}],RandomReal[{-0.0001,0.0001}]},1],numbodies]*)
pos0 = Table[
Partition[{RandomReal[{0.5, 1}], RandomReal[{0.5, 1}],
RandomReal[{-0.4, -0.04}]}, 1], numbodies];

eforceX[x_, y_, z_] = q*Grad[sol[x, y, z], {x, y, z}][[1]];
eforceY[x_, y_, z_] = q*Grad[sol[x, y, z], {x, y, z}][[2]];
eforceZ[x_, y_, z_] = q*Grad[sol[x, y, z], {x, y, z}][[3]];

tfin = 400;

eqs = Table[{x1[j]''[t] == -(1/mass)*eforceX[x, y, z] /. {x ->
x1[j][t], y -> y1[j][t], z -> z1[j][t]},
y1[j]''[t] == -(1/mass)*eforceY[x, y, z] /. {x -> x1[j][t],
y -> y1[j][t], z -> z1[j][t]},
z1[j]''[t] == -(1/mass)*eforceZ[x, y, z] /. {x -> x1[j][t],
y -> y1[j][t], z -> z1[j][t]},
x1[j][0] == pos0[[j, 1, 1]],
y1[j][0] == pos0[[j, 2, 1]],
z1[j][0] == pos0[[j, 3, 1]],
x1[j]'[0] == vel0[[j, 1, 1]],
y1[j]'[0] == vel0[[j, 2, 1]],
z1[j]'[0] == vel0[[j, 3, 1]]},
{j, numbodies}];

vars = Flatten[Table[{x1[j][t], y1[j][t], z1[j][t]}, {j, numbodies}]];

event = Table[{WhenEvent[
z1[j][t] >= -partRad , {x1[j]'[t] -> 0, y1[j]'[t] -> 0,
z1[j]'[t] -> 0}]} /. {j -> i}, {i, numbodies}];

sol1 = NDSolve[{eqs, event}, Head /@ vars, {t, 0, tfin}]

frames = Table[
Show[v, ParametricPlot3D[
Table[{x1[j][t], y1[j][t], z1[j][t]} /. sol1, {j,
numbodies}], {t, 0, tf}, PlotRange -> All, Axes -> Off],
Graphics3D[
Table[{Hue[.35],
Sphere[{x1[j][tf], y1[j][tf], z1[j][tf]} /. sol1, 0.025]}, {j,
numbodies}]]], {tf, .025*tfin, tfin, .025*tfin}];
video = ListAnimate[frames, SaveDefinitions -> True]

• Just a quick hint, which I hope is accurate: To stop a variable, say x[t], governed by x'[t] == vx[t] for some expression vx[t], use a discrete variable x'[t] == a[t] * vx[t], with initial value a[0] == 1. Then use WhenEvent[x[t] > upperboundary, a[t] -> 0] to stop x[t]. You'll need a discrete variable for each particle, and you'll need to use it to stop each coordinate of a given particle. Commented Jul 5, 2022 at 16:30
• (You'll also need a WhenEvent for each possible boundary crossing, potentially as many as two times the number of coordinates times the number of particles.) Commented Jul 5, 2022 at 16:48
• @MichaelE2 That makes sense to me, but I suppose implementation is what is confusing. For instance, in my set of eqs, I made z1[j]''[t] == a[j][t]*(-1/mass)*eforceZ[x, y, z] and event = Table[ WhenEvent[ z1[j][t] >= -partRad , {a[j][t] -> 0} /. {j -> i}], {i, numbodies}]; sol1 = NDSolve[{eqs, event}, Head /@ vars, {t, 0, tfin}, DiscreteVariables -> a]... This throws an error "The function a[j] appears as the head of the expression a[j][t]".
– Zach
Commented Jul 5, 2022 at 17:53
• You need to set both {a[j][t] -> 0, z1[j]'[t] -> 0} since it's a 2nd order equation. Also the j -> i should go outside: Table[WhenEvent[ ...] /. {j -> i}, {i, 2}] Commented Jul 5, 2022 at 18:20
• @MichaelE2 The answer is shown below. Thank you for your help!
– Zach
Commented Jul 5, 2022 at 19:41

Here is the solution. Thanks to @Michael E2 in original question comments.

Clear["Global*"];
Needs["NDSolveFEM"]

q = 1.617733*10^-18;(*Net ion Charge*)

(*R= Import["https://www.dropbox.com/s/dds8rm3odg2m7gu/largeAp.obj?dl=\
1"];*)

R = Import["C:\\Users\\zaust\\OneDrive\\Desktop\\largeAp.obj"];
RegionDimension[R];
M = BoundaryMeshRegion[MeshCoordinates[R], MeshCells[R, 2]];
RegionDimension[M];
Volume[M];

r = RegionDifference[
RegionDifference[
RegionDifference[Cuboid[{0, 0, -6.5}, {2, 2, 0.5}], M],
Cuboid[{0, 0, 0.4}, {2, 2, 0.5}]],
Cuboid[{0, 0, -6.5}, {2, 2, -6.4}]];
ToElementMesh[r]["Wireframe"];
pol = -1;

V0 = 1000;
sol = NDSolveValue[{Laplacian[V[x, y, z], {x, y, z}] == 0,
DirichletCondition[
V[x, y, z] == -pol*
V0/2, (0.4 <= z <= 0.5) && (0 <= y <= 2) && (0 <= x <= 2)],
DirichletCondition[
V[x, y, z] ==
pol*V0/2, (0.0071 <= z <= 0.0072) && (0 <= y <= 2) && (0 <= x <=
2)], DirichletCondition[
V[x, y, z] ==
0, (0 <= z <= 0.0070) && (0 <= y <= 2) && (0 <= x <= 2)],
DirichletCondition[
V[x, y, z] ==
0, (-6.5 <= z <= -6.4) && (0 <= y <= 2) && (0 <= x <= 2)]},
V, {x, y, z} \[Element] r];

electricField[x_, y_, z_] = -Grad[sol[x, y, z], {x, y, z}];

v = Show[VectorPlot3D[
electricField[x, y, z], {x, 0.5, 1}, {y, 0.5, 1}, {z, -0.5, 0.1},
PlotTheme -> "Detailed", ColorFunction -> "Rainbow",
PerformanceGoal -> "Quality", VectorScale -> 0.05,
VectorPoints -> 7, PlotLegends -> Automatic], M];

vecForce =
Show[VectorPlot3D[
q*electricField[x, y, z], {x, 0.5, 1}, {y, 0.5,
1}, {z, -0.25, -0.00001}, PlotTheme -> "Detailed",
ColorFunction -> "Rainbow", PerformanceGoal -> "Quality",
VectorScale -> 0.05, VectorPoints -> 7, PlotLegends -> Automatic],
M];

numbodies = 3;
mass = 6.52*10^-11;(*particle mass in kg/m^3*)

vel0 = Table[Partition[{0, 0, 0}, 1], numbodies];
(*vel0=Table[Partition[{RandomReal[{-0.0001,0.0001}],RandomReal[{-0.\
0001,0.0001}],RandomReal[{-0.0001,0.0001}]},1],numbodies]*)
pos0 = Table[
Partition[{RandomReal[{0.5, 1}], RandomReal[{0.5, 1}],
RandomReal[{-0.4, -0.04}]}, 1], numbodies];
discreteStart = Table[Partition[{1, 1, 1}, 1], numbodies];

eforceX[x_, y_, z_] = q*Grad[sol[x, y, z], {x, y, z}][[1]];
eforceY[x_, y_, z_] = q*Grad[sol[x, y, z], {x, y, z}][[2]];
eforceZ[x_, y_, z_] = q*Grad[sol[x, y, z], {x, y, z}][[3]];

tfin = 600;

eqs = Table[{x1[j]''[t] == -a[t]*(1/mass)*eforceX[x, y, z] /. {x ->
x1[j][t], y -> y1[j][t], z -> z1[j][t]},
y1[j]''[t] == -a[t]*(1/mass)*eforceY[x, y, z] /. {x -> x1[j][t],
y -> y1[j][t], z -> z1[j][t]},
z1[j]''[t] == -a[t]*(1/mass)*eforceZ[x, y, z] /. {x -> x1[j][t],
y -> y1[j][t], z -> z1[j][t]},
x1[j][0] == pos0[[j, 1, 1]],
y1[j][0] == pos0[[j, 2, 1]],
z1[j][0] == pos0[[j, 3, 1]],
x1[j]'[0] == vel0[[j, 1, 1]],
y1[j]'[0] == vel0[[j, 2, 1]],
z1[j]'[0] == vel0[[j, 3, 1]],
a[0] == 1},
{j, numbodies}];

vars = Flatten[Table[{x1[j][t], y1[j][t], z1[j][t]}, {j, numbodies}]];

,{x1[j]'[t]\[Rule] 0,y1[j]'[t]\[Rule] 0,z1[j]'[t]\[Rule] 0}]}/.{j\
\[Rule]i},{i,numbodies}];*)

event = Table[
WhenEvent[
z1[j][t] >= -partRad , {a[t] -> 0, x1[j]'[t] -> 0,
y1[j]'[t] -> 0, z1[j]'[t] -> 0}] /. {j -> i}, {i, numbodies}];

sol1 = NDSolve[{eqs, event}, Head /@ vars, {t, 0, tfin},
DiscreteVariables -> a]

frames = Table[
Show[v, ParametricPlot3D[
Table[{x1[j][t], y1[j][t], z1[j][t]} /. sol1, {j,
numbodies}], {t, 0, tf}, PlotRange -> All, Axes -> Off],
Graphics3D[
Table[{Hue[.35],
Sphere[{x1[j][tf], y1[j][tf], z1[j][tf]} /. sol1, 0.025]}, {j,
numbodies}]]], {tf, .015*tfin, tfin, .015*tfin}];
video = ListAnimate[frames, SaveDefinitions -> True]
$$$$
`