5
$\begingroup$

I have a code (provided below) which simulates the motion of an ensemble of charged particles which are subjected to a complex static electric field which passes through some aperture in a metallic surface (N spheres, whose number is defined as numbodies in my code).

Now, I would like to run the simulation for a single particle many, many times, and plot each of their final resting places on the surface (z=-partRad) of my metalized aperture at (x(tfin), y(tfin), z= -partRad) together, in 2D. These coordinates can be displayed as simple dots in the x-y plane.

It would also be cool to superimpose the surface of my aperture onto the x-y graph of the final particle locations, to see exactly where they most frequently land. I would also like to find the average distance these particles center from some user specified point(s), after everything is said and done. I would also appreciate any cool ideas for statistical analysis regarding the outcome of the Monte Carlo simulation!!!

Each of the particles should be given a random starting position and velocity given by

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];

in my code.

The first part of the code, above the (***... MONTE CARLO CAN START...***) line in my code, doesn't need to be included in the Monte Carlo simulation (unless you find some reason I do not see).

Here is the code. Happy hunting!

Clear["Global`*"];
Needs["NDSolve`FEM`"]

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];

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

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]];

(***********************    MONTE CARLO CAN START    ******************************)

numbodies = 1;

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];

tfin = 1000;

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}];

partRad = 0.0125;

vars = Flatten[Table[{x1[j][t], y1[j][t], z1[j][t]}, {j, 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}];

ListAnimate[Rasterize[#, "Image"] & /@ frames]
Improve this question
$\endgroup$
3
  • $\begingroup$ Hi @DanielHuber. Are you running into a problem trying to execute the code? I copy/pasted the exact code provided above, and it runs without error on my end. $\endgroup$
    – Zach
    Commented Aug 1, 2022 at 14:02
  • 1
    $\begingroup$ Sorry, I must have done something wrong by copy/past. There was a Table command without iterator. $\endgroup$
    – Daniel Huber
    Commented Aug 1, 2022 at 15:16
  • 1
    $\begingroup$ @zak720 As I remember we have discussed once this mesh with two holes on mathematica.stackexchange.com/questions/269806/… $\endgroup$
    – Alex Trounev
    Commented Aug 2, 2022 at 10:03

1 Answer 1

Reset to default
4
$\begingroup$

To simulate hundreds and thousands runs we need very fast code, while at present form code for trajectories is very slow. To improve code, first, we include the electric field computation in NDSolve as follows

Clear["Global`*"];
Needs["NDSolve`FEM`"]

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]];
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}]];
pol = -1;

V0 = 1000;

{U, {Fx, Fy, Fz}} = 
  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, -Grad[V[x, y, z], {x, y, z}]}, {x, y, z} \[Element] 
    r];

Visualization

v = Show[
  VectorPlot3D[{Fx, Fy, Fz}, {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], 
  RegionPlot3D[R, PlotStyle -> Opacity[0.5]]]

Figure 1

Second, we use Verlet integration algorithm to simulate trajectories, we have

mass = 6.52*10^-11;(*particle mass*)
SeedRandom[1234];
n = 100; X0 = Take[RandomReal[{.5, 1}, 2 n], n]; Y0 = 
 Take[RandomReal[{.5, 1}, 2 n], -n]; Z0 = 
 RandomReal[{-.4, -.04}, n]; rrv = 
 RandomReal[{-0.0001, 0.0001}, 3 n]; u0 = Take[rrv, n]; v0 = 
 Take[rrv, {n + 1, 2 n}]; w0 = Take[rrv, -n];


SymplecticLeapfrog = {"SymplecticPartitionedRungeKutta", 
  "DifferenceOrder" -> 2, "PositionVariables" :> qvars}; time = {t, 0,
   2}; tm = 10^3; partRad = 0.0125; qvars = {X[t], Y[t], Z[t]};
eqn[i_] := {X''[t] == q/mass tm^2 Fx /. {x -> X[t], y -> Y[t], 
    z -> Z[t]}, 
  Y''[t] == q/mass tm^2 Fy /. {x -> X[t], y -> Y[t], z -> Z[t]}, 
  Z''[t] == q/mass tm^2 Fz /. {x -> X[t], y -> Y[t], z -> Z[t]}, 
  X[0] == X0[[i]], Y[0] == Y0[[i]], Z[0] == Z0[[i]], X'[0] == u0[[i]],
   Y'[0] == v0[[i]], Z'[0] == w0[[i]], 
  WhenEvent[Z[t] >= -partRad, end = t; "StopIntegration"]}; h = 10^-4;
 Do[soln[i] = 
   NDSolve[eqn[i], qvars, time, StartingStepSize -> h, 
    Method -> SymplecticLeapfrog]; 
  T[i] = end;, {i, 1, n}]

It takes about 25 sec to compute 100 trajectories on my laptop. Visualization of trajectories

Show[v, Table[
  ParametricPlot3D[qvars /. soln[i], {t, 0, T[i]}, 
   ColorFunction -> Hue], {i, 100}]]

Figure 2

Finally we can compute point distribution as

fin = Table[qvars /. soln[i][[1]] /. t -> T[i], {i, n}];
Show[Graphics3D[{{Red, Point[fin]}}, 
  PlotRange -> {{.5, 1}, {.5, 1}, All}], 
 RegionPlot3D[R, PlotStyle -> Opacity[0.5]]]

Figure 3

To compute statistics we need 10^4 points as minimum. Example with 1000 points Figure 4

Example with $10^4$ points Figure 5

In the last case 33 particles of 10000 are not deposited on the plane z==-partRad, therefore we select cases and plot as follows

fin = Table[qvars /. soln[i][[1]] /. t -> T[i], {i, n}];

fin1 = Select[fin, Norm[#[[3]] + partRad] <= .01 &];

Show[Graphics3D[{{Red, Point[fin1]}}, 
  PlotRange -> {{.5, 1}, {.5, 1}, All}], 
 RegionPlot3D[R, PlotStyle -> Opacity[0.5]]]

Using fin1 we can compute mean and standard deviation on every coordinate

XYmean = Mean[fin1]

(*Out[]= {0.757097, 0.754464, -0.0125015}*)

 std = StandardDeviation[fin1]

(*Out[]= {0.0511615, 0.045477, 0.000108942}*)

Visualization of mean and std[[1]] as a sphere

Show[Graphics3D[{{Red, Point[fin1]}, {Blue, Opacity[0.25], 
    Sphere[XYmean, std[[1]]]}}, PlotRange -> {{.5, 1}, {.5, 1}, All}, 
  Axes -> {True, True, False}], 
 RegionPlot3D[R, PlotStyle -> Opacity[0.5]]]

Figure 6

Improve this answer
$\endgroup$
13
  • $\begingroup$ this is an absolute work of art. I had come up with a very crude version of this prior to seeing this answer (and it is, as you might guess, extremely slow). Fantastic work. $\endgroup$
    – Zach
    Commented Aug 5, 2022 at 1:11
  • $\begingroup$ How would I set a WhenEvent such that when a particle comes into contact with any other particle position previously simulated, its motion comes to a stop? In my crude version, I made an empty list called outCoord and filled it with the final resting coordinates of each particle simulated. I then check whether, WhenEvent[ Sqrt[(outCoord[[k, 1]][[1, 1]] - x1[j][t])^2 + (outCoord[[k, 1]][[1, 2]] - y1[j][t])^2 + (outCoord[[k, 1]][[1, 3]] - z1[j][t])^2] <= partRad for each new particle simulated, such that if this is true, it stops. How can I do it here? $\endgroup$
    – Zach
    Commented Aug 5, 2022 at 3:19
  • 1
    $\begingroup$ I don't understand last question. Do you mean collisions or something else? $\endgroup$
    – Alex Trounev
    Commented Aug 5, 2022 at 7:56
  • $\begingroup$ Yes like a collision. If a particle encounters a region of space where previously simulated particles have landed, it's motion stops. This would create a 3D structure of particles on the surface since they would now be able to "stack" on top of eachother. $\endgroup$
    – Zach
    Commented Aug 5, 2022 at 12:37
  • 1
    $\begingroup$ Ok, for this we need to define some size of particle and rule for packaging. $\endgroup$
    – Alex Trounev
    Commented Aug 5, 2022 at 12:43

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.