Adds props argument to calculate_FTLE

Planktos.py

@@ -68,6 +68,7 @@ def __init__(self, Lx=10, Ly=10, Lz=None,
             char_L: characteristic length scale. Used for calculating Reynolds
                 number, especially in the case of immersed structures (ibmesh)
                 and/or when simulating inertial particles
+            U: characteristic fluid speed. Used for some calculations.
             init_swarms: initial swarms in this environment
             units: length units to use, default is meters. Note that you will
                 manually need to change self.g (accel due to gravity) if using
@@ -815,7 +816,7 @@ def read_IB2d_vtk_data(self, path, dt, print_dump,
 
         ### Convert environment dimensions and reset simulation time ###
         self.L = [self.flow_points[dim][-1] for dim in range(2)]
-        self.__reset_flow_variables(incl_rho_mu_U=True)
+        self.__reset_flow_variables()
         self.reset()
 
 
@@ -845,7 +846,7 @@ def read_IBAMR3d_vtk_data(self, filename):
 
         ### Convert environment dimensions and reset simulation time ###
         self.L = [self.flow_points[dim][-1] for dim in range(3)]
-        self.__reset_flow_variables(incl_rho_mu_U=True)
+        self.__reset_flow_variables()
         # record the original lower left corner (can be useful for later imports)
         self.fluid_domain_LLC = (mesh[0][0], mesh[1][0], mesh[2][0])
         # reset time
@@ -923,7 +924,7 @@ def read_IBAMR3d_vtk_dataset(self, path, start=None, finish=None):
 
         ### Convert environment dimensions and reset simulation time ###
         self.L = [self.flow_points[dim][-1] for dim in range(3)]
-        self.__reset_flow_variables(incl_rho_mu_U=True)
+        self.__reset_flow_variables()
         # record the original lower left corner (can be useful for later imports)
         self.fluid_domain_LLC = (mesh[0][0], mesh[1][0], mesh[2][0])
         # reset time
@@ -969,7 +970,7 @@ def read_npy_vtk_data(self, path, prefix='', flow_file='flow.npz',
 
         ### Convert environment dimensions and reset simulation time ###
         self.L = [self.flow_points[dim][-1] for dim in range(len(self.flow))]
-        self.__reset_flow_variables(incl_rho_mu_U=True)
+        self.__reset_flow_variables()
         self.reset()
 
 
@@ -1013,7 +1014,7 @@ def read_comsol_vtu_data(self, filename, vel_conv=None, grid_conv=None):
 
         ### Convert environment dimensions and reset simulation time ###
         self.L = [self.flow_points[dim][-1] for dim in range(3)]
-        self.__reset_flow_variables(incl_rho_mu_U=True)
+        self.__reset_flow_variables()
         # record the original lower left corner (can be useful for later imports)
         self.fluid_domain_LLC = (mesh[0][0], mesh[1][0], mesh[2][0])
         # reset time
@@ -1634,7 +1635,8 @@ def get_mean_fluid_speed(self):
 
 
     def calculate_FTLE(self, grid_dim=None, testdir=None, t0=0, T=0.1, dt=0.001, 
-                       ode=None, t_bound=None, swrm=None, params=None):
+                       ode_gen=None, props=None, t_bound=None, swrm=None, 
+                       params=None):
         '''Calculate the FTLE field at the given time(s) t0 with integration 
         length T on a discrete grid with given dimensions. The calculation will 
         be conducted with respect to the fluid velocity field loaded in this 
@@ -1677,14 +1679,18 @@ def calculate_FTLE(self, grid_dim=None, testdir=None, t0=0, T=0.1, dt=0.001,
                 this argument represents the length of the Euler time steps.
             t_bound: if solving ode or tracer particles, this is the bound on
                 the RK45 integration step size. Defaults to dt/100.
-            ode: [optional] function handle for an ode to be solved specifying 
-                deterministic equations of motion. Should have call signature 
-                ODEs(t,x), where t is the current time (float) and x is a 2*NxD 
-                array with the first N rows giving v=dxdt and the second N rows 
-                giving dvdt. D is the spatial dimension of the problem. See the 
-                ODE generator functions in motion.py for examples of valid functions. 
-                The passed in ODEs will be solved using a grid of initial conditions 
-                with Runge-Kutta.
+            ode_gen: [optional] function handle for an ode generator that takes
+                in a swarm object and returns an ode function handle with
+                call signature ODEs(t,x), where t is the current time (float) 
+                and x is a 2*NxD array with the first N rows giving v=dxdt and 
+                the second N rows giving dvdt. D is the spatial dimension of 
+                the problem. See the ODE generator functions in motion.py for 
+                examples of format. 
+                The ODEs will be solved using RK45 with a newly created swarm 
+                specified on a grid throughout the domain.
+            props: [optional] dictionary of properties for the swarm that will 
+                be created to solve the odes. Effectively, this passes parameter 
+                values into the ode generator.
             swarm: [optional] swarm object with user-defined movement rules as 
                 specified by the get_positions method. This allows for arbitrary 
                 FTLE calculations through subclassing and overriding this method. 
@@ -1714,7 +1720,8 @@ def calculate_FTLE(self, grid_dim=None, testdir=None, t0=0, T=0.1, dt=0.001,
             grid_dim = tuple(len(pts) for pts in self.flow_points)
 
         if swrm is None:
-            s = swarm(envir=self, init='grid', grid_dim=grid_dim, testdir=testdir)
+            s = swarm(envir=self, shared_props=props, init='grid', 
+                      grid_dim=grid_dim, testdir=testdir)
             # NOTE: swarm has been appended to this environment!
         else:
             # Get a soft copy of the swarm passed in
@@ -1747,12 +1754,12 @@ def calculate_FTLE(self, grid_dim=None, testdir=None, t0=0, T=0.1, dt=0.001,
         #   components for efficiency.
         if swrm is None:
             ### SETUP SOLVER ###
-            if ode is None:
+            if ode_gen is None:
                 ode_fun = motion.tracer_particles(s, incl_dvdt=False)
                 print("Finding {}D FTLE based on tracer particles.".format(len(grid_dim)))
             else:
-                ode_fun = ode
-                print("Finding {}D FTLE based on supplied ODE.".format(len(grid_dim)))
+                ode_fun = ode_gen(s)
+                print("Finding {}D FTLE based on supplied ODE generator.".format(len(grid_dim)))
             print(prnt_str)
 
             # keep a list of all times solved for 
@@ -1764,7 +1771,7 @@ def calculate_FTLE(self, grid_dim=None, testdir=None, t0=0, T=0.1, dt=0.001,
             while current_time < T:
                 new_time = min(current_time + dt,T)
                 ### TODO: REDO THIS SOLVER!!!!!!!!
-                if ode is None:
+                if ode_gen is None:
                     y = s.positions[~s.positions[:,0].mask,:]
                 else:
                     y = np.concatenate((s.positions[~s.positions[:,0].mask,:], 
@@ -1780,7 +1787,7 @@ def calculate_FTLE(self, grid_dim=None, testdir=None, t0=0, T=0.1, dt=0.001,
                 # Put current position in the history (maybe only do this if something exits??)
                 s.pos_history.append(s.positions.copy())
                 # pull solution into swarm object's position/velocity attributes
-                if ode is None:
+                if ode_gen is None:
                     s.positions[~s.positions[:,0].mask,:] = y_new
                 else:
                     N = round(y_new.shape[0]/2)

motion.py

@@ -347,11 +347,17 @@ def inertial_particles(swarm):
     try:
         R = swarm.get_prop('R')
     except KeyError:
-        rho_p = swarm.get_prop('rho')
-        rho_f = swarm.envir.rho
-        R = 2*rho_f/(rho_f+2*rho_p)
+        try:
+            rho_p = swarm.get_prop('rho')
+            rho_f = swarm.envir.rho
+            R = 2*rho_f/(rho_f+2*rho_p)
+        except KeyError:
+            raise KeyError("Could not find required physical property R or rho in swarm object.")
 
-    a = swarm.get_prop('diam')*0.5 # radius of particles
+    try:
+        a = swarm.get_prop('diam')*0.5 # radius of particles
+    except KeyError:
+        raise KeyError("Could not find required physical property 'diam' in swarm object.")
 
     assert swarm.envir.char_L is not None, "Characteristic length scale in envir not specified."
     L = swarm.envir.char_L

tests/visualtest_ftle_2d.py

@@ -5,18 +5,23 @@
 import numpy as np
 import sys
 sys.path.append('..')
-import Planktos
+import Planktos, motion
 
 # TODO: get a less-difficult geometry
-envir = Planktos.environment()
+envir = Planktos.environment(char_L=0.1, rho=1, mu=0.001, U=15)
 # envir.read_IB2d_vtk_data('data/channel_cyl', 5.0e-5, 1000, d_start=1)
 envir.read_IB2d_vtk_data('data/channel_cyl', 5.0e-5, 1000, d_start=20, d_finish=20)
-# envir.read_IB2d_vertex_data('data/channel_cyl/channel.vertex')
+envir.read_IB2d_vertex_data('data/channel_cyl/channel.vertex')
 
 # s = envir.add_swarm(900, init='grid', grid_dim=(30,30), testdir='x1')
 
 # Test basic tracer particles on a masked mesh.
 # envir.calculate_FTLE((102,25),T=0.1, dt=0.001, testdir='x1')
-sf, time_list, last_time = envir.calculate_FTLE((512,128),T=0.1,dt=0.001)
+# sf, time_list, last_time = envir.calculate_FTLE((512,128),T=0.1,dt=0.001)
+
+
+sf, time_list, last_time = envir.calculate_FTLE((512,128),T=0.1,dt=0.001,
+                                                ode_gen=motion.inertial_particles,
+                                                props={'R':2/3, 'diam':0.01})
 
 envir.plot_2D_FTLE()