environment class

Every instance of the environment class is functionally a rectangular spatial domain in either two or three dimensions. The lower left corner is located at the Euclidean origin. Boundary conditions are specified with respect to the agents on each side of the domain. A fluid velocity field can be specified on a regular mesh of grid points which always includes the domain boundaries. The fluid velocity may vary in time, but the spatial mesh on which it is specified must remain constant. Analytical fluid velocity fields are also available.

Created on Tues Jan 24 2017

Author: Christopher Strickland

Email: cstric12@utk.edu

class planktos.environment(Lx=10, Ly=10, Lz=None, x_bndry='zero', y_bndry='zero', z_bndry='noflux', flow=None, flow_times=None, rho=None, mu=None, nu=None, char_L=None, U=None, init_swarms=None, units='m', ibmesh_color=None)

Rectangular environment containing fluid info, immersed meshs, and swarms.

The environment class does much of the heavy lifting of Planktos. It loads and contains information about the fluid velocity field, the dimensions of the physical environment being simulated, boundary conditions for the agents, the agent swarms that are in the environment, any immersed meshes, and all simulation time information. Additionally, it provides functions for manipulating the fluid velocity field in certain ways (e.g. extending, tiling, and interpolating), calculating vorticity and FTLE, viewing info about the fluid itself, and calling the move function on all swarms contained in the environment. It is essential to any Planktos simulation and typically the first Planktos object you create.

Parameters:

• Lx (float, default=10) – Length of domain in x and y direction, meters

• Ly (float, default=10) – Length of domain in x and y direction, meters

• Lz (float, optional) – Length of domain in z direction

• x_bndry ({'zero', 'noflux', 'periodic'} as str or [str, str], default='zero') – agent boundary condition in the x-axis (if the same on both sides), or [left bndry condition, right bndry condition]. Choices are ‘zero’, ‘noflux’, and ‘periodic’. Agents leaving a zero boundary condition will be marked as masked and cease to be updated or plotted afterward. In the noflux case, agents will undergo a sliding collision with the boundary. Movement that would have occurred through the boundary will be projected onto the boundary instead. In the periodic case, agents leaving one side of the domain will reenter on the other side.

• y_bndry (str or [str, str], default='zero') – agent boundary condition in the y-axis (if the same on both sides), or [left bndry condition, right bndry condition].

• z_bndry (str or [str, str], default='noflux') – agent boundary condition in the z-axis (if the same top and bottom), or [low bndry condition, high bndry condition].

• flow (list of ndarrays, optional) – This only needs to be specified if you already have a fluid velocity field loaded as a list of numpy arrays and wish to add it to the environment directly. Otherwise, the fluid velocity field will be assumed to be zero everywhere until something is loaded or created via a method of this class. If you specify a fluid velocity field here, it should be of the following format: [x-vel field ndarray ([t],i,j,[k]), y-vel field ndarray ([t],i,j,[k]), z-vel field ndarray if 3D]. Note! i is x index, j is y index, with the value of x and y increasing as the index increases. It is assumed that the flow mesh is equally spaced and includes values on the domain boundary. A keyword argument ‘flow_points’ must also be specified as a tuple (len==dimension) of 1D arrays specifying the mesh points along each direction. If the velocity field is time varying, the argument ‘flow_points’ must also be given, which should in this case be an interable of times at which the fluid velocity is specified (as indexed by t in the first dimension of each of the fluid ndarrays).

• flow_times ([float, float] or increasing iterable of floats, optional) – [tstart, tend] or iterable of times at which the fluid velocity is specified or scalar dt; required if flow is time-dependent.

• rho (float, optional) – fluid density of environment, kg/m**3 (optional, m here meaning length units). Auto-calculated if mu and nu are provided.

• mu (float, optional) – dynamic viscosity, kg/(m*s), Pa*s, N*s/m**2 (optional, m here meaning length units). Auto-calculated if rho and nu are provided.

• nu (float, optional) – kinematic viscosity, m**2/s (optional, m here meaning length units). Auto-calculated if rho and mu are provided.

• char_L (float, optional) – characteristic length scale. Used for calculating Reynolds number, especially in the case of immersed structures (ibmesh) and/or when simulating inertial particles

• U (float, optional) – characteristic fluid speed. Used for some calculations.

• init_swarms (swarm object or list of swarm objects, optional) – initial swarms in this environment. Can be added later.

• units (string, default='m') – length units to use, default is meters. Note that you will manually need to change self.g (accel due to gravity) if using something else. Also, none of the methods of this class use this attribute in any way, so it’s probably best to work in meters.

• ibmesh_color (matplotlib color format, optional) – color of the ibmesh. Defaults to black in 2D, ‘dimgrey’ in 3D.

Attributes

Llist of floats

length of the domain [x, y, [z]] in the stated units

unitsstr

A place to remind yourself of the spatial units you are working in. Note that none of the methods of this class use this attribute in any way, so it’s probably best to work in meters.

bndrylist of lists, each with two of {‘zero’, ‘noflux’, ‘periodic’}

Boundary conditions in each direction [x, y, [z]] for agents

swarmslist of swarm objects

The swarms that belong to this environment

timefloat

current environment time

time_historylist of floats

list of past time states

flowlist of ndarrays or fCubicSpline objects

[x-vel field ndarray ([t],i,j,[k]), y-vel field ndarray ([t],i,j,[k]),

z-vel field ndarray (if 3D)]. i is x index, j is y index, with the value of x and y increasing as the index increases. Arrays get replaced by fCubicSpline objects (if the fluid velocity is temporally varying) when they are first needed.

flow_timesndarray of floats or None

if specified, the time stamp for each index t in the flow arrays (time varying fluid velocity fields only)

flow_pointstuple (len==dimension) of 1D ndarrays

points defining the spatial grid for the fluid velocity data

fluid_domain_LLCtuple (len==dimension) of floats

original lower-left corner of domain (if from data)

tilingtuple (x,y) of floats

how much tiling was done in the x and y direction

orig_Ltuple (Lx,Ly) of floats

length of the domain in x and y direction (Lx,Ly) before tiling occured

h_pfloat

height of porous region in analytic flows, e.g. Brinkman

gfloat

accel due to gravity (length units/s**2). Only active for models of motion that include gravity (this is not the default)

rhofloat

fluid density, kg/m**3

mufloat

dynamic viscosity, kg/(m*s)

nufloat

kinematic viscosity, m**2/s

char_Lfloat

characteristic length, m

Ufloat

characteristic fluid speed, m/s

netcdfnetCDF4 Dataset object

only present if a netCDF file has been loaded for reading

ibmeshndarray

This is either an Nx2x2 array (2D, line mesh) or an Nx3x3 array (3D, triangular mesh) specifying internal mesh structures that agents will treat as solid barriers. Each value of the first index into the array specifies a mesh element, and each 1x2 or 1x3 contained in that index is a point in space (so, two for a line, three for a triangle). Avoid mesh structures that intersect with a periodic boundary - behavior related to this is not implemented.

max_meshpt_distfloat

maximum length of a mesh segment in ibmesh, typically determined automatically. This is used in search algorithms to winnow down the number of mesh elements to search for intersections of movement.

ibmesh_colormatplotlib color format

color of the ibmesh. Defaults to black in 2D, ‘dimgrey’ in 3D.

plot_structslist of function handles

List of functions that plot additional environment structures which the agents are unaware of. E.g., for visual purposes only.

plot_structs_argsList of tuples

List of argument tuples to be passed to plot_structs functions, after ax (the plot axis object) is passed

FTLE_largestndarray

FTLE field calculated using the largest eigenvalue

FTLE_smallestndarray

FTLE field calculated using the smallest eigenvalue. take the negative of this to get backward-time information

FTLE_locNx2 masked ndarray

spatial points on which the FTLE mesh was calculated

FTLE_loc_endNx2 masked ndarray

final locations for each of the FTLE mesh points

FTLE_t0float

start-time for the FTLE calculation

FTLE_Tfloat

integration time for FTLE

FTLE_grid_dimtuple of int

grid dimension for the FTLE calculation (x,y,[z])

mag_gradndarray

acts as a cache for the gradient of magnitude of the fluid velocity

mag_grad_timefloat

simulation time at which magnitude gradient above was calculated

DuDtlist of ndarrays

material derivative cache

DuDt_timefloat

simulation time at which material derivative was calculated

dt_interplist of PPoly objects

Used for temporal derivative interpolation. Set by dudt method.

Examples

Create a 3D environment that is 10x10x10 meters with fluid density and dynamic viscosity recorded. The fluid velocity is zero everywhere, but can be set to something different later.

>>> envir = planktos.environment(Lz=10, rho=1000, mu=1000)

Create a 2D 5x3 meter environment with zero fluid velocity.

>>> envir = planktos.environment(Lx=5, Ly=3)

property Re

Return the Reynolds number at the current time. Must have set self.U, self.char_L and self.nu. Reynolds number is U*char_L/nu.

add_swarm(swarm_size=100, **kwargs)

Adds a swarm into this environment.

Parameters:

• swarm_size (swarm object or int (size of swarm), default=100) – If a swarm object is given, all following parameters will be ignored (since the object is already initialized)

• init (string) – Method for initializing positions. See swarm class for options.

• seed (int) – Seed for random number generator

• kwargs – keyword arguments to be set as swarm properties (see swarm class for details)

add_vertices_to_2D_ibmesh()

Methods for auto-connecting mesh vertices into line segments may result in line segments that cross each other away from vertex points. This will cause undesirable behavior including mesh crossing. This method adds a vertex point at any such intersection to avoid such problems.

calculate_DuDt(t_indx=None, time=None)

Calculate and store material derivative of the fluid velocity with respect to time. Defaults to interpolating at the current time, given by self.time. Gradient is calculated via second order accurate central differences (using numpy) with second order accuracy at the boundaries. The material derivative is saved in case it is needed again.

The material derivative is given by .. math:: frac{Dmathbf{u}}{Dt} = mathbf{u}_t + (nablamathbf{u})mathbf{u}

Parameters:

• t_indx (int, optional) – Interpolate at a time referred to by self.envir.time_history[t_indx]

• time (float, optional) – Interpolate at a specific time. default is current time.

calculate_FTLE(grid_dim=None, testdir=None, t0=0, T=0.1, dt=0.001, 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 environment and either tracer particle movement (default), an ode specifying deterministic equations of motion, or other arbitrary particle movement as specified by a swarm object’s get_positions method and updated in discrete time intervals of length dt.

This method will set the following environment attributes: - environment.FTLE_largest - environment.FTLE_smallest - environment.FTLE_loc - environment.FTLE_t0 - environment.FTLE_T - environment.FTLE_grid_dim

All FTLE calculations will be done using a swarm object. This means that:

1. The boundary conditions specified by this environment will be respected.

2. Immersed boundaries (if any are loaded into this environment) will be treated as impassible to all particles and movement vectors crossing these boundaries will be projected onto them.

If passing in a set of ode or finding the FTLE field for tracer particles, an RK45 solver will be used. Otherwise, integration will be via the swarm object’s get_positions method.

If both ode and swarm arguments are None, the default is to calculate the FTLE based on massless tracer particles.

Parameters:

• grid_dim (tuple of int) – size of the grid in each dimension (x, y, [z]). Defaults to the fluid grid.

• testdir (str) – grid points can heuristically be removed from the interior of immersed structures. To accomplish this, a line will be drawn from each point to a domain boundary. If the number of intersections is odd, the point is considered interior and masked. See grid_init for details - this argument sets the direction that the line will be drawn in (e.g. ‘x1’ for positive x-direction, ‘y0’ for negative y-direction). If None, do not perform this check and use all gridpoints.

• t0 (float, optional) – start time for calculating FTLE. If None, default behavior is to set t0=0. TODO: Interable to calculate at many times. Default then becomes t0=0 for time invariant flows and calculate FTLE at all times the flow field was specified at (self.flow_times) for time varying flows?

• T (float, default=0.1) – integration time. Default is 1, but longer is better (up to a point), unless smallest FTLE is desired and agents are leaving the domain…

• dt (float, default=0.001) – if solving ode or tracer particles, this is the time step for checking boundary conditions. If passing in a swarm object, this argument represents the length of the Euler time steps.

• ode_gen (function handle, optional) – functional 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 (dict, optional) – dictionary of properties for the swarm that will be created to solve the odes. Effectively, this passes parameter values into the ode generator. If unspecified, will default to the values for the first agent in the props of the swarm provided.

• t_bound (float, optional) – if solving ode or tracer particles, this is the bound on the RK45 integration step size. Defaults to dt/100.

• swrm (swarm object, 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. Steps of length dt will be taken until the integration length T is reached. The swarm object itself will not be altered; a shallow copy will be created for the purpose of calculating the FTLE on a grid.

• params (dict, optional) – params to be passed to supplied swarm object’s get_positions method.

Returns:

• swarm object – used to calculuate the FTLE

• list – list of dt integration steps

• ndarray – array the same size as the point grid giving the last time in the integration before each point exited the domain

calculate_mag_gradient()

Calculate and store the gradient of the magnitude of the fluid velocity at the current time, along with the time at which it was calculated. Gradient is calculated via second order accurate central differences (using numpy) with second order accuracy at the boundaries and saved in case it is needed again.

center_cell_regrid(periodic_dim=(False, False, False))

Re-grids data that was specified at the center of cells instead of at the corners.

NOTE! This needs to be called before any immersed meshes are loaded. It will NOT look for and properly shift these meshes.

Software has a tendency to output data files where the fluid mesh is specified at the center of cells rather than at the corners. This will be readily apparent if Planktos loads your fluid velocity data and reports spacial dimensions one dx, dy, and dz smaller than you were expecting. To fix this, Planktos will interpolate/extrapolate the fluid velocity mesh using the default method to get additional grid points on the edge of the domain.

Periodicity can be enforced in specified dimensions.

Parameters:

periodic_dim (list-like of 2 or 3 bool, default=(True, True, False)) – True if that spatial dimension is periodic, otherwise False. The 3rd entry will be ignored in the 2D case.

define_pic_grid(positions, dx, dy, dz=None, return_neighbors=True)

Creates a grid across the domain with cells of size dx by dy for use in a particle-in-cell method. dx and dy should divide the domain evenly in their respective directions and represent the furthest away one needs to look from any individual agent in order to get all neighbor interactions (e.g., a characteristic distance).

Returns a dictionary of cells in which keys are (i,j) tuples indexing the cells starting at zero from the origin, and the values they point to are lists of agent indices whose positions are within that cell.

If return_neighbors is True, will also return a dictionary of cell indices in which the values are of agents located either within that cell OR in a neighboring cell. Neighbor cells are the 8 cells vertically or horizontally adjacent or diagonally adjacent (in 2D). Adjacency on the boundaries of the domain depends upon the environment boundary condition: zero or no-flux will treat the edge of the environment as a hard boundary while periodic will wrap around to find neighboring cells.

Parameters:

• positions (Nx2 or Nx3 ndarray) – ndarray of agent positions (e.g., swarm.positions)

• dx (float) – length of grid cell in the x-direction

• dy (float) – length of grid cell in the y-direction

• dz (float, optional) – length of grid cell in the z-direction

• return_neighbors (bool, default=True) – if True, return both a dictionary a dictionary of cells -> agent indices located inside cell AND a dictionary of cells -> agent indices in both cell and neighboring cells

Returns:

cells dictionary, or tuple of two dictonaries (cells, neighbors)

dudt(t_indx=None, time=None)

Return the derivative of the fluid velocity with respect to time. Defaults to interpolating at the current time, given by self.time.

Parameters:

• t_indx (int, optional) – Interpolate at a time referred to by self.envir.time_history[t_indx]

• time (float, optional) – Interpolate at a specific time. default is current time.

Returns:

List of ndarrays

extend(x_minus=0, x_plus=0, y_minus=0, y_plus=0)

Duplicate the boundary of the fluid flow a number of times in the x (+ or -) and/or y (+ or -) directions, thus extending the domain with constant fluid velocity. Good for extending domains with resolved fluid flow before/after and on the sides of a structure.

Parameters:

• x_minus (int) – number of times to duplicate bndry in the x- direction

• x_plus (int) – number of times to duplicate bndry in the x+ direction

• y_minus (int) – number of times to duplicate bndry in the y- direction

• y_plus (int) – number of times to duplicate bndry in the y+ direction

get_2D_vorticity(t_indx=None, time=None, t_n=None)

Calculuate the vorticity of the fluid velocity field at a given time.

If all time arguments are None but the flow is time-varying, the vorticity at the current time will be returned. If more than one time argument is specified, only the first will be used.

Parameters:

• t_indx (int) – integer time index into self.envir.time_history[t_indx]

• time (float) – time

• t_n (int) – integer time index into self.flow_times[t_n]

Returns:

vorticity as an ndarray

get_mean_fluid_speed()

Return the mean fluid speed at the current time, temporally interpolating the flow field if necessary.

interpolate_flow(positions, flow=None, flow_points=None, time=None, method='linear')

Spatially interpolate the fluid velocity field (or another flow field) at the supplied positions. If flow is None and self.flow is time-varying, the flow field will be interpolated in time first, using the current environmental time, or a different time if provided.

Parameters:

• positions (array) – NxD locations at which to interpolate the flow field, where D is the dimension of the system.

• flow (list-like of arrays, optional) – The fluid velocity data, with each array representing one spatial component of the velocity vector. The first dimension of each array is time if the fluid velocity field is time varying. If None, the environmental flow field is used. Interpolated in time if necessary.

• flow_points (list-like of 1-D arrays, optional) – The set of coordinates along each dimension that defines the mesh grid for the fluid velocity field. If None, the environmental flow points is used.

• time (float, optional) – if None, the present time. Otherwise, the flow field will be interpolated to the time given based on the environment flow_times. This is only supported for environmental flow fields (not ones passed in as an argument).

• method (string, default='linear') – spatial interpolation method to be passed to scipy.interpolate.interpn. Anything but splinef2d is supported.

interpolate_temporal_flow(t_indx=None, time=None)

Interpolate flow in time using a cubic spline. Defaults to interpolating at the current time, given by self.time.

Parameters:

• t_indx (int) – Interpolate at a time referred to by self.envir.time_history[t_indx]

• time (float) – Interpolate at a specific time

Returns:

interpolated flow field as a list of ndarrays

load_NetCDF(filename)

Load a NetCDF file into the netcdf attribute of the environment. Does not automatically read in any data.

Because NetCDF files can contain multiple data sets with different dimension names and associated metadata, and because it may be desirable to explore the data set first and/or load only a subset of the data, this method just loads the Dataset into the environment object. See the documentation/tutorial for netCDF4 on ways to read the metadata for the loaded Dataset. See read_NetCDF_flow for reading in data from a loaded NetCDF dataset.

Parameters:

filename (string) – path and filename of the NetCDF file, including extension

See also

read_NetCDF_flow

read in data from a loaded NetCDF dataset

move_swarms(dt=1.0, params=None, ib_collisions='sliding', silent=False)

Move all swarms in the environment.

Parameters:

• dt (float) – length of time step to move all agents

• params (any, optional) – parameters to pass along to get_positions, if necessary

• ib_collisions ({None, 'sliding' (default), 'sticky'}) – Type of interaction with immersed boundaries. If None, turn off all interaction with immersed boundaries. In sliding collisions, conduct recursive vector projection until the length of the original vector is exhausted. In sticky collisions, just return the point of intersection.

• silent (bool, default=False) – If True, suppress printing the updated time.

plot_2D_FTLE(smallest=False, clip_l=None, clip_h=None, figsize=None)

Plot the FTLE field as generated by the calculate_FTLE method. The field will be hard to visualize in 3D, so only 2D is implemented here. For 3D visualization, output the field as a vtk and visualize using VisIt, ParaView, etc.

TODO: Show a video of 2D slices as a plot of 3D FTLE

Parameters:

• smallest (bool, default=False) – If true, plot the negative, smallest, forward-time FTLE as a way of identifying attracting Lagrangian Coherent Structures (see Haller and Sapsis 2011). Otherwise, plot the largest, forward-time FTLE as a way of identifying ridges (separatrix) of LCSs.

• clip_l (float, optional) – lower clip value (below this value, mask points)

• clip_h (float, optional) – upper clip value (above this value, mask points)

• figsize (tuple, optional) – matplotlib figsize

plot_2D_vort(t=None, clip=None, interval=500, figsize=None)

Plot the vorticity of a 2D fluid at the given time t or at all times if t is None. This method will calculate the vorticity on the fly.

Clip will limit the extents of the color scale.

For time dependent velocity fields, interval is the delay between plotting of each time’s flow data, in milliseconds. Defaults to 500.

plot_envir(figsize=None)

Plot the environment without the flow, e.g. to verify ibmesh formed correctly, dimensions are correct, etc.

plot_flow(t=None, downsamp=5, interval=500, figsize=None, **kwargs)

Plot the velocity field of the fluid at a given time t or at all times if t is None. If t is not in self.flow_times, the nearest time will be shown without interpolation.

For densely sampled velocity fields, specify an int for downsamp to plot every nth vector in each direction.

For time dependent velocity fields, interval is the delay between plotting of each time’s flow data, in milliseconds. Defaults to 500.

Extra keyword arguments will be passed to pyplot’s quiver.

2D arrow lengths are scaled based on the maximum of all the data over all times.

read_IB2d_vertex_data(filename, res_factor=0.501, res=None)

Reads in 2D vertex data from a .vertex file (IB2d). Assumes that any vertices closer than res_factor (default is half + a bit for numerical stability) times the Eulerian mesh resolution are connected linearly. This method gets the Eulerian mesh resolution from the fluid velocity data, so the flow data must be imported first.

Alternatively, you can pass in the Eulerian mesh resolution directly to the res parameter to get the mesh in absence of any fluid velocity field. Calculate it from the IB2d input file by taking the domain length and dividing by the number of Eulerian grid points: Lx/Nx and Ly/Ny. The smaller of the two numbers should be used.

Avoid mesh structures that intersect with a periodic boundary; behavior related to this is not implemented.

Parameters:

• filename (string) – vertex file to load

• res_factor (float, default=0.501) – this times the Eulerian mesh resolution is the radius that will be used.

• res (float, optional) – pass the Eulerian mesh resolution in directly, instead of calculating it from a loaded fluid velocity field

read_IB2d_vtk_data(path, dt, print_dump, d_start=0, d_finish=None)

Reads in vtk flow velocity data generated by IB2d and sets environment variables accordingly.

Can read in vector data with filenames u.####.vtk or scalar data with filenames uX.####.vtk and uY.####.vtk.

IB2d is an Immersed Boundary (IB) code for solving fully coupled fluid-structure interaction models in Python and MATLAB. The code is hosted at https://github.com/nickabattista/IB2d

Parameters:

• path (str) – path to folder with vtk data

• dt (float) – dt in input2d

• print_dump (int) – print_dump in input2d

• d_start (int, default=0) – number of first vtk dump to read in

• d_finish (int, optional) – number of last vtk dump to read in, or None to read to end

• auto_regrid (bool, default=True) – IB2d always has periodic BC and returns a VTK with fluid specified at the center of the mesh cells. Regrid based on this information to make the fluid periodic within Planktos and to fill out the domain.

read_IBAMR3d_vtk_data(filename, vel_conv=None, grid_conv=None)

Reads in vtk flow data from a single source and sets environment variables accordingly. The resulting flow will be time invarient. It is assumed this data is a rectilinear grid.

All environment variables will be reset.

Parameters:

• filename (string) – filename of data to read, incl. file extension

• vel_conv (float, optional) – scalar to multiply the velocity by in order to convert units

• grid_conv (float, optional) – scalar to multiply the grid by in order to convert units

read_IBAMR3d_vtk_dataset(path, start=None, finish=None)

Reads in vtk flow data generated by VisIt from IBAMR and sets environment variables accordingly. Assumes that the vtk filenames are IBAMR_db_###.vtk where ### is the dump number, as automatically done with read_IBAMR3d_py27.py. Also assumes that the mesh is the same for each vtk.

Imported times will be translated backward so that the first time loaded corresponds to an agent environment time of 0.0.

Parameters:

• path (string) – path to vtk data

• start (int, optional) – vtk file number to start with. If None, start at first one.

• finish (int, optional) – vtk file number to end with. If None, end with last one.

read_NetCDF_flow(flow_x, flow_y=None, flow_z=None, vec_idx=None, dim_reorder=None, x_name=None, y_name=None, z_name=None, time_name=None, conv_time=False)

Read NetCDF fluid data into the environment. Must first have loaded a NetCDF dataset with load_NetCDF.

The default expectation is that the x, y, and z components of the fluid velocity data are specified in separate Variables. If that is not the case, then the index of the Variable dimension which specifies the component of the vector must be supplied in the parameter vec_idx (base 0). It will be assumed that the ordering in this Dimension is x, y, then z components.

This method assumes that mesh point values in the x, y, and z directions are given by Variables which have the same name as the Dimensions of the fluid flow Variables. If so, the method will automatically find and load them. Otherwise, strings must be specified for x_name, y_name, and z_name giving the Variable names. The same goes for the time mesh data.

Blurb about time conversion

If the coordinate variables have a Unit attribute, this will be loaded as the environment’s units.

Parameters:

• flow_x (string) – Variable name (or path, if the variable is inside group) to the fluid velocity data for the x-direction, or for all directions.

• flow_y (string, optional) – Variable name (or path, if the variable is inside group) to the fluid velocity data for the y-direction. If all directions are in the same variable, this does not need to be supplied, but vec_idx must be supplied instead.

• flow_z (string, optional) – Similar to flow_y.

• vec_idx (int, optional) – Index of the variable dimension along which the different components of the velocity are given. Leave as None if the components are given in separate variables.

• dim_reorder (tuple or list of ints, optional) – Fluid velocity data is stored within this class in the dimensional ordering ([time], x, y, [z]). If this matches the NetCDF ordering, then no adjustment is necessary. Otherwise, this must be specified as a tuple or list which contains a permutation of [0,1,..,N-1] where N is the number of spatial-temporal dimensions and the numbers correspond to where each dimension of the NetCDF variable should go in the Planktos ordering. If the NetCDF variable has a dimension specifying the component of the fluid velocity, ignore it and do not include it in the permutation list.

• x_name (string, optional) – Variable name (or path, if the variable is inside a group) for the coordinate variable corresponding with the x spatial direction. Only necessary if different from the NetCDF Dimension name.

• y_name (string, optional) – Similar to x_name.

• z_name (string, optional) – Similar to x_name.

• time_name (string, optional) – Variable name (or path, if the variable is inside a group) for the coordinate variable corresponding with the time dimension. Only necessary if different from the NetCDF Dimension name.

• conv_time (bool, default=False) – If the time variable is not numerical, set this to True in order to convert from a time format relative to a fixed date using a certain calendar to floating point numerical values. In this case, it will be expected that the time variable is a sequence of datetime objects, so conversion of this variable from a string to the python datetime object may be required first. Uses the date2num function provided by cftime. See Notes below for expected formatting in this case.

See also

load_NetCDF

Load a NetCDF file into the netcdf environment attribute.

Notes

In order for the date2num conversion to work, units and a calendar must be specified in the coordinate variable for time. Specifically, the attribute of the time variable named “units” must be a string of the form “<time units> since <reference time>”. <time units> can be days, hours, minutes, seconds, milliseconds or microseconds. <reference time> is the time origin. months since is allowed only for the 360_day calendar and common_years since is allowed only for the 365_day calendar. The attribute of the time variable named “calendar” must be a string descrbing the calendar to be used in the time calculations. All values in the CF metadata convention are supported. Valid calendars ‘standard’, ‘gregorian’, ‘proleptic_gregorian’ ‘noleap’, ‘365_day’, ‘360_day’, ‘julian’, ‘all_leap’, ‘366_day’. Default is None which means the calendar associated with the first input datetime instance will be used.

read_comsol_vtu_data(filename, vel_conv=None, grid_conv=None)

Reads in vtu flow data from a single source and sets environment variables accordingly. The resulting flow will be time invarient. It is assumed this data is on a regular grid and that a grid section is included in the data.

FOR NOW, THIS IS TIME INVARIANT ONLY.

All environment variables will be reset.

Parameters:

• filename (string) – filename of data to read, incl. file extension

• vel_conv (float, optional) – scalar to multiply the velocity by in order to convert units

• grid_conv (float, optional) – scalar to multiply the grid by in order to convert units

read_stl_mesh_data(filename, unit_conv=None)

Reads in 3D mesh data from an ascii or binary stl file. Must have the numpy-stl library installed. It is assumed that the coordinate system of the stl mesh matches the coordinate system of the flow field. Thus, the mesh will be translated using the flow LLC if necessary.

Avoid mesh structures that intersect with a periodic boundary; behavior related to this is not implemented.

Parameters:

• filename (string) – filename of data to read, incl. file extension

• unit_conv (float, optional) – scalar to multiply the mesh by in order to convert units

read_vertex_data_to_convex_hull(filename)

Reads in 2D or 3D vertex data from a vtk file or a vertex file and applies ConvexHull triangulation to get a complete boundary. This uses Qhull through Scipy under the hood http://www.qhull.org/.

regenerate_flow_data()

Regenerates the original fluid velocity data from temporal spline objects.

reset(rm_swarms=False)

Resets environment to time=0. Swarm history will be lost, and all swarms will maintain their last position and velocities. If rm_swarms=True, remove all swarms.

save_2D_vorticity(path, name, time_history=True, flow_times=False)

Save the vorticity of a 2D fluid velocity field as one or more vtk files (one for each time point).

Parameters:

• path (string) – location to save file(s)

• name (string) – prefix name for file(s)

• time_history (bool) – if True, save vorticity data for each time step in the simulation history. Only for time-varying fluid.

• flow_times (bool) – if True, save vorticity data for each time at which the fluid velocity data is explicitly specified.

save_fluid(path, name, time_history=True, flow_times=False)

Save the fluid velocity field as one or more vtk files (one for each time point).

Parameters:

• path (str) – directory in which to store the files. if it does not exist, it will be created

• name (str) – prefix to put on filenames

• time_history (bool) – if True, save for each time step in the simulation history, interpolating as needed. Only for time-varying fluid.

• flow_times (bool) – if True, save at each time for which the fluid velocity data is explicitly specified.

set_boundary_conditions(x_bndry, y_bndry, z_bndry=None)

Check and set boundary conditions. Set Z-dimension only if zdim is not None. Each boundary condition must be either a list or an iterable of length 2.

set_brinkman_flow(alpha, h_p, U, dpdx, res=101, tspan=None)

Get a fully developed Brinkman flow with a porous region.

This method sets the environment fluid velocity as a 1D Brinkman flow based on a porous layer of hight h_p in the bottom of the domain. Velocity gradient is zero in the x-direction and all flow moves parallel to the x-axis. Porous region is the lower part of the y-domain (2D) or z-domain (3D) with width h_p and an empty region above. For 3D flow, the velocity profile is the same on all slices y=c. The decision to set 2D vs. 3D flow is based on the current dimension of the environment.

After this method is successfully called, the flow property of the environment class will be set to the resulting Brinkman flow, and h_p will be set in the environment’s properties.

Parameters:

• alpha (float) – equal to 1/(hydraulic permeability). alpha=0 implies free flow (infinitely permeable)

• h_p (float) – height of porous region

• U (float or list of floats) – velocity at top of domain (v in input3d of IB2d). If a list of floats, will create a time varying fluid velocity field with Brinkman flow matching each U at a series of time points. The time points are determined by tspan.

• dpdx (float or list of floats) – dp/dx change in momentum constant. if a list, will correspond to a time varying flow field with those values of dp/dx.

• res (int) – resolution of the flow; that is, number of points at which to resolve the flow, including boundaries

• tspan ([float, float] or iterable of floats, optional) – corresponds to [tstart, tend] (start time and end time with an evenly spaced time mesh) or an iterable of times at which flow is specified in the case of a time-varying flow field. if not specified and U/dpdx are iterable, dt=1 will be used with a start time of zero.

Examples

Create a 3D environment with time varying Brinkman flow

>>> envir = planktos.environment(Lz=10, rho=1000, mu=1000)
>>> U=0.1*np.array(list(range(0,5))+list(range(5,-5,-1))+list(range(-5,8,3)))
>>> envir.set_brinkman_flow(alpha=66, h_p=1.5, U=U, dpdx=np.ones(20)*0.22306,
    res=101, tspan=[0, 20])

Notes

Brinkman’s equation [1] is written as

\[\rho(\mathbf{u}_{t}(\mathbf{x},t) +\mathbf{u}(\mathbf{x},t)\cdot\nabla \mathbf{u}(\mathbf{x},t)) = -\nabla p(\mathbf{x},t) + \mu \nabla^2 \mathbf{u}(\mathbf{x},t) - \alpha \mu \mathbf{u}(\mathbf{x},t)\]

where is \(\alpha\) is the inverse of the hydraulic permeability. We take a region of height h_p where \(\alpha>0\) with parallel shear flow above it, and we assume that the flow is steady (\(\partial u/\partial t=0\)), fully developed (\(\partial u/\partial x=0\)), and zero in all cross-stream directions. In this case, the equations can be reduced to an analytical solution, which is what we evaluate here. See [2] for more information.

References

[1]

H.C. Brinkman, (1949). “A calculation of the viscous force exerted by a flowing fluid on a dense swarm of particles,” Applied Scientific Research 1, 27(34).

[2]

C. Strickland, L.A. Miller, A. Santhanakrishnan, C. Hamlet, N.A. Battista, V. Pasour (2017). “Three-dimensional low Reynolds number flows near biological filtering and protective layers,” Fluids, 2(62).

See also

set_two_layer_channel_flow, set_canopy_flow

set_canopy_flow(h, a, u_star=None, U_h=None, beta=0.3, C=0.25, res=101, tspan=None)

Apply flow within and above a uniform homogenous canopy according to the model described in Finnigan and Belcher (2004), “Flow over a hill covered with a plant canopy” [4].

The decision to set 2D vs. 3D flow is based on the current dimension of the environment. Default values for beta and C are based on Finnigan & Belcher [4]. Must specify two of u_star, U_h, and beta, though beta has a default value of 0.3 so just giving u_star or U_h will also work. If one of u_star, U_h, or beta is given as a list-like object, the flow will vary in time.

Parameters:

• h (float) – height of canopy (m)

• a (float) – leaf area per unit volume of space $m^{-1}$. Typical values are a=1.0 for dense spruce to a=0.1 for open woodland

• u_star (float, optional (may set U_h instead)) – canopy friction velocity. u_star = U_h*beta if not set

• U_h (float, optional (may set u_star instead)) – wind speed at top of canopy. U_h = u_star/beta if not set

• beta (float, default=0.3) – mass flux through the canopy (u_star/U_h)

• C (float, default=0.25) – drag coefficient of indivudal canopy elements

• res (int) – number of points at which to resolve the flow, including boundaries

• tspan ([float, float] or iterable of floats, optional) – [tstart, tend] or iterable of times at which flow is specified if None and u_star, U_h, and/or beta are iterable, dt=1 will be used.

See also

set_brinkman_flow, set_two_layer_channel_flow

Notes

In addition to self.flow, this will set self.h_p = h

References

[4] (1,2)

J.J. Finnigan and S.E. Belcher, (2004). “Flow over a hill covered with a plant canopy,” Quarterly Journal of the Royal Meteorological Society, 130(596), 1-29.

set_two_layer_channel_flow(a, h_p, Cd, S, res=101)

Apply wide-channel flow with vegetation layer according to the two-layer model described in Defina and Bixio (2005), “Vegetated Open Channel Flow” [3].

The decision to set 2D vs. 3D flow is based on the current dimension of the environment and the fluid velocity is always time-independent.

Parameters:

• a (float) – vegetation density, given by Az*m, where Az is the frontal area of vegetation per unit depth and m the number of stems per unit area (1/m), assumed constant

• h_p (float) – plant height (m)

• Cd (float) – drag coefficient, assumed uniform (unitless)

• S (float) – bottom slope (unitless, 0-1 with 0 being no slope, resulting in no flow)

• res (int) – number of points at which to resolve the flow, including boundaries

See also

set_brinkman_flow, set_canopy_flow

Notes

In addition to self.flow, this will set self.h_p = a

References

[3]

A. Defina and A.C. Bixio, (2005). “Mean flow and turbulence in vegetated open channel flow,” Water Resources Research, 41(7).

tile_flow(x=2, y=1)

Tile fluid flow and immersed meshes a number of times in the x and/or y directions. While obviously this works best if the fluid is periodic in the direction(s) being tiled, this will not be enforced. Instead, it will just be assumed that the domain edges are equivalent, and only the right/upper domain edge will be used in tiling.

Parameters:

• x (int, default=2) – number of tiles in the x direction (counting the one already there)

• y (int, default=1) – number of tiles in the y direction (counting the one already there)

wrap_flow(periodic_dim=(True, True, False))

In some cases, software may print out fluid velocity data that omits the velocities at the right boundaries in spatial dimensions that are meant to be periodic. This helper function restores that data by copying everything over. 3rd dimension will automatically be ignored if 2D.

Parameters:

periodic_dim (list of 3 bool, default=[True, True, False]) – True if that spatial dimension is periodic, otherwise False