motion module

Library for different sorts of particle motion, including implementations of various numerical methods for relevant equations of motion. Most of these will take in a swarm object from which information about the particles and their environment can be accessed along with relevant parameters. They will then return the new particle positions or Delta x after a time dt, depending on implementation.

Created on Tues Jan 24 2017

Author: Christopher Strickland

Email: cstric12@utk.edu

planktos.motion.Euler_brownian_motion(swarm, dt, mu=None, ode=None, sigma=None)

Uses the Euler-Maruyama method to solve the Ito SDE

\[dX_t = \mu(X_t,t) dt + \sigma(t) dW_t\]

where \(\mu\) is the drift and \(\sigma\) is the diffusion. \(\mu\) can be specified directly as a constant or via an ode, or both. If both are ommited, the default is to use the mu property of the swarm object plus the local fluid drift. If an ode is given but mu is not, the the mu property of the swarm will be added to the ode velocity term before solving (however, the swarm mu is zero by default).

\(\sigma\) can be provided a number of ways (see below), and can be dependent on time or agent if passed in directly. This solver is only order 0.5 if \(\sigma\) is dependent on spatial position, so this is not directly supported.

Parameters:

• swarm (swarm object) –

• dt (float) – time step to take

• mu (1D array of length D, array of shape NxD, or array of shape 2NxD, optional) – drift velocity as a 1D array of length D (spatial dimension), an array of shape NxD (where N is the number of agents), or as an array of shape 2NxD, where N is the number of agents and D is the spatial dimension. In the last case, the first N rows give the velocity and the second N rows are the acceleration. In this case, a Taylor series method Euler step will be used. If no mu and no ode is given, a default brownian drift of swarm.get_prop(‘mu’) + the local fluid drift will be used. If mu=None but an ode is given, the default for mu will be swarm.get_prop(‘mu’) alone, as fluid interaction is assumed to be handled by the ode. Note that the default swarm value for mu is zeros, so the ode will specify the entire drift term unless mu is set to something else.

• ode (callable, optional) – an ODE function for mu with call signature func(t,x), where t is the current time and x is a 2NxD array of positions/velocities. It must return a 2NxD array of velocities/accelerations. See mu for information on default behavior if this is not specified.

• sigma (array, optional) –

  Brownian diffusion coefficient matrix. If None, use the ‘cov’ property of the swarm object, or lacking that, the ‘D’ property. For convienence, \(\sigma\) can be provided in several ways:

  • As a covariance matrix stored in swarm.get_prop(‘cov’). This is the default. The matrix given by this swarm property is assumed to be defined by \(\sigma\sigma^T\) and independent of time or spatial location. The result of using this matrix is that the integrated Wiener process over an interval dt has covariance swarm.get_prop(‘cov’)*dt, and this will be fed directly into the random number generator to produce motion with these characteristics. It should be a square matrix with the length of each side equal to the spatial dimension, and be symmetric.

  • As a diffusion tensor (matrix). The diffusion tensor is given by \(D = 0.5*\sigma\sigma^T\), so it is really just half the covariance matrix. This is the diffusion tensor as given in the Fokker-Planck equation or the heat equation. As in the case of the covariance matrix, it is assumed constant in time and space, and should be specified as a swarm property with the name ‘D’. Again, it will be fed directly into the random number generator. It should be a square matrix with the length of each side equal to the spatial dimension, and be symmetric.

  • Given directly as an argument to this function. In this case, it should be an DxD array, or an NxDxD array where N is the number of agents and D the spatial dimension. Since this is an Euler step solver, sigma is assumed constant across dt.

Returns:

NxD array (N number of agents, D spatial dimension) – New agent positions after the Euler step.

Notes

TODO: If sigma is spatially dependent, in order to maintain strong order 1 convergence we need the method due to Milstein (1974) (see Kloeden and Platen). This might be useful for turbulant models, or models where the energy of the random walk is dependent upon local fluid properties. This would be particularly useful in the case of media with different densities or in the case of a chemical concentration that excites or depresses behavior. When implementing, the details of this function will need to be changed to directly accept a spatially varying sigma. We should implement a Stratonovich function call at that time too, because with a spatially dependent sigma, Ito and Stratonovich are no longer equivalent.

planktos.motion.RK45(fun, t0, y0, tf, rtol=0.0001, atol=1e-06, h_start=None)

Runge-Kutta Dormand-Prince [1] solver (variable time-step solver).

The passed in ode function (fun) must have call signature (t,x) where x is a 2-D array with a number of columns equal to the spatial dimension. The solver will run to tf and then return, after which boundary conditions can be checked within a swarm object before restarting at the next time step.

Parameters:

• fun (callable) – Right-hand side of the ODE system. The call signature must be fun(t, y), where t is a scalar (time) and y is an NxD array where D is the dimension of the system.

• t0 (float) – initial time.

• y0 (ndarray of shape (N,D)) – initial state.

• tf (float) – final time, e.g. time to integrate to.

• rtol (float, defaults=0.0001, 1e-06) – relative and absolute tolerance. The solver keeps the local error estimates less than atol + rtol * abs(y).

• atol (float, defaults=0.0001, 1e-06) – relative and absolute tolerance. The solver keeps the local error estimates less than atol + rtol * abs(y).

• h_start (float, optional) – time step-size to attempt first. Also the maximum step size to use. Defaults to (tf-t0)*0.5.

Returns:

y (ndarray with shape matching y0) – new state at time tf

References

[1]

Dormand, J.R., Prince, P.J. (1980). A family of embedded Runge-Kutta formulae, Journal of Computational and Applied Mathematics, 6(1), 19-26.

planktos.motion.flatten_ode(swarm)

Defines a decorator capable of converting a flattened, passed in x-variable array into a 2NxD shape for generated ODE functions, and then take the result of the ODE functions and reflatten. Need knowledge of the dimension of swarm for this, so the swarm must be passed to the decorator. 2N accounts for N equations giving the derivative of position w.r.t time and N equations giving the derivative of velocity. D is the dimension of the problem (2D or 3D).

planktos.motion.highRe_massive_drift(swarm)

Function generator for ODEs governing high Re massive drift with drag as formulated by Nathan et al. [4]. Assumes Re > 10 with neutrally buoyant particles possessing mass and a known cross-sectional area given by the property cross_sec.

\[f_D = \frac{\rho C_d A}{m}||\mathbf{u}-\mathbf{v}|| (\mathbf{u} -\mathbf{v})\]

where \(f_D\) is the drag force acting on the particle, \(C_d\) is the drag coefficient, \(A\) is the cross sectional area, \(\rho\) is the fluid density, and \(\mathbf{u}\) is the fluid velocity. The system of ODEs are specified so that the acceleration is defined to be \(f_D/m\).

Parameters:

swarm (swarm object) –

Returns:

callable, func(t,x)

Notes

Requires that the following are specified in either swarm.shared_props (if uniform across agents) or swarm.props (for individual variation):

• m: mass of each agent

• Cd: Drag coefficient acting on cross-sectional area

• cross_sec: cross-sectional area of each agent

Requires that the following are specified in the fluid environment:

• rho: fluid density

References

[4]

R. Nathan, G.G. Katul, G. Bohrer, A. Kuparinen, M.B. Soons, S.E. Thompson, A. Trakhtenbrot, H.S. Horn (2011). Mechanistic models of seed dispersal by wind. Theoretical Ecology, 4(2), 113-132

planktos.motion.inertial_particles(swarm)

Function generator for ODEs governing small, rigid, spherical particles whose dynamics can be described by the linearized Maxey-Riley equation [2] described in Haller and Sapsis (2008) [3].

\[\begin{split}\frac{d\mathbf{v}}{dt} &= \frac{3R}{2}\frac{D\mathbf{u}}{Dt} - \mu(\mathbf{v}-\mathbf{u}) + \left(1-\frac{3R}{2}\right)\mathbf{g} \\ \text{with}& \\ \frac{D\mathbf{u}}{Dt} &= \mathbf{u}_t + (\nabla\mathbf{u})\mathbf{u}\end{split}\]

Critically, it is assumed that \(\mu = R/St\) is much greater than 1, where R is the density ratio \(R=2\rho_f/(\rho_f+2\rho_p)\), and St is the Stokes number.

Parameters:

swarm (swarm object) –

Returns:

callable, func(t,x)

Notes

Requires that the following are specified in either swarm.shared_props (if uniform across agents) or swarm.props (for individual variation):

• rho or R: particle density or density ratio, respectively. if supplying R, 0<R<2/3 corresponds to aerosols, R=2/3 is neutrally buoyant, and 2/3<R<2 corresponds to bubbles.

• diam: diameter of particle

Requires that the following are specified in the fluid environment:

• char_L: characteristic length scale for Reynolds number calculation

• nu: kinematic viscosity

• U: characteristic fluid speed

• g: acceleration due to gravity (set by default to 9.80665)

• rho: fluid density (unless R is specified in swarm)

References

[2]

Maxey, M.R. and Riley, J.J. (1983). Equation of motion for a small rigid sphere in a nonuniform flow. Phys. Fluids, 26(4), 883-889.

[3]

Haller, G. and Sapsis, T. (2008). Where do inertial particles go in fluid flows? Physica D: Nonlinear Phenomena, 237(5), 573-583.

planktos.motion.tracer_particles(swarm, incl_dvdt=True)

Function generator for ODEs describing tracer particles.

Parameters:

• swarm (swarm object) –

• incl_dvdt (bool, default=True) – Whether or not to include equations for dvdt so that x has shape 2NxD matching most other ODEs (dvdt will just be given as 0). If False, will expect an NxD array for x instead of a 2NxD array.

Returns:

callable, func(t,x)