From bf84b0e58aebe3b3a39e5f2326ff61f42d363c0e Mon Sep 17 00:00:00 2001 From: Christopher Strickland Date: Wed, 2 Jun 2021 14:35:52 -0400 Subject: [PATCH] Adds props argument to calculate_FTLE --- Planktos.py | 47 +++++++++++++++++++++---------------- motion.py | 14 +++++++---- tests/visualtest_ftle_2d.py | 13 ++++++---- 3 files changed, 46 insertions(+), 28 deletions(-) diff --git a/Planktos.py b/Planktos.py index 30b0656..3c31f51 100644 --- a/Planktos.py +++ b/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) diff --git a/motion.py b/motion.py index 741b461..983932f 100644 --- a/motion.py +++ b/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 diff --git a/tests/visualtest_ftle_2d.py b/tests/visualtest_ftle_2d.py index dd6f11e..638e767 100644 --- a/tests/visualtest_ftle_2d.py +++ b/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() \ No newline at end of file