Banjo
Banjo

Reputation: 13

Using conditions inside events with solve_ivp triggers the event at the wrong time

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

Answers (0)

Related Questions