Reputation: 13
I'm trying to plot the dynamics of a particle moving through several cylindrical magnetic elements. I'm using solve_ivp to integrate the ODE, but I'm having trouble obtaining the expected behavior using a stop event: I'd like the integration to stop when the particle gets close to the inner radius of one of the elements (90% of the radius, in my case), provided it's within this element.
Here is how I write my stop event for the radial stop event:
def stop_event(self, t, y, collection):
"""
Event function to stop the integration when the particle's radial distance r
exceeds the diameter of any magnetic element within its x-bounds.
"""
x, y_pos, z = y[0], y[1], y[2]
r = np.sqrt(y_pos**2 + z**2) # radial distance r = sqrt(y^2 + z^2)
# Loop through the magnetic elements to check if the particle is in any element's x-range
for element in collection:
mag_element = element['magnetic_element']
x_entrance = mag_element.entrance_x # Start x of the magnetic element
x_exit = mag_element.exit_x # End x of the magnetic element
radius = mag_element.innermost_loop_diam / 2 # Critical radial boundary (90% of element radius)
# Check if particle is within the x-bounds of the magnetic element
if x_entrance <= x <= x_exit:
return (r - 0.9 * radius) # Event triggered when this crosses 0
return 1 # No event if not within the bounds of any magnetic element
And here is a shortened version of the integration method:
def solve_trajectory(self, collection, t_span, nb_points, max_step = None, h=1e-6):
y0 = np.concatenate([self.initial_position, self.initial_velocity]) # initial state [x, y, z, vx, vy, vz]
# Define the stopping event
def radial_stop_event(t, y):
return self.stop_event(t, y, collection)
radial_stop_event.terminal = True # The event should terminate the integration
radial_stop_event.direction = 1 # Detect crossing in from inside to outside
[...]
# Compute dt from the number of points
dt = (t_span[1] - t_span[0]) / (nb_points - 1)
# Call solve_ivp with Runge-Kutta method (RK45)
sol = solve_ivp(
fun=lambda t, y: self.equations_of_motion(t, y, collection, h=h),
t_span=t_span,
y0=y0,
method='RK45',
events = [radial_stop_event, longitudinal_enter_event, longitudinal_exit_event],
dense_output=True,
max_step = max_step,
rtol=1e-8,
atol=1e-10
)
[...]
return sol, times_in_elements
I unfortunately cannot provide a full working code, as this would be too long and convoluted.
As of now, the stop event works as intended if the particle does "crash" inside the magnetic element: the integration is stopped.
However, if the particle didn't crash (r
stayed smaller than 0.9 * radius
during the travel of the particle inside the magnetic element), the event triggers and stops the integration when x reaches x_exit.
I don't understand why this is happening. I've tried printing the values of the condition before and after the event triggers (by inserting print(f"t={t}, x={x}, r={r}, radius={0.9 * radius}, condition = {r - 0.9 * radius}
, and making the event not terminal), and this is what it returns just before the event is triggered :
t=9.203132978272565e-05, x=0.003033979317911664, r=1.4922512843731781e-15, radius=0.00135, condition = -0.0013499999999985078
and right after:
t=0.0002882113689884196, x=0.020175674729679623, r=7.060896498215925e-15, radius=0.00135, condition = -0.001349999999992939
so nothing that would normally cause solve_ivp to stop the integration, i.e. this is not an issue with the particle breaching the r - 0.9 * radius < 0
condition just when x > x_exit
due to the adaptative step size or whatnot.
I've also tried to set a fairly restrictive value of max_step
, used dense_output = True
, rewrote the condition with r - 0.9 * radius > 0
instead or with tolerances, nothing worked, and I'm out of ideas.
Why does reaching x_exit also trigger the event, and is there any workaround?
I could register the events and not make them terminal, only to shorten sol
a posteriori, but, as the integration is fairly costly, it'd be better if the terminal event worked as intended.
Upvotes: 1
Views: 42