{'fraction_non_alphanumeric': 0.026056492226844755, 'fraction_numerical': 0.01131304284358806, 'mean_word_length': 5.0385632437197, 'pattern_counts': {'":': 0, '<': 3, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 3, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Optimization of micropillar sequences for fluid flow sculpting

Inertial fluid flow deformation around pillars in a microchannel is a new method for controlling fluid flow. Sequences of pillars have been shown to produce a rich phase space with a wide variety of flow transformations. Previous work has successfully demonstrated manual design of pillar sequences to achieve desired transformations of the flow cross-section, with experimental validation. However, such a method is not ideal for seeking out complex sculpted shapes as the search space quickly becomes too large for efficient manual discovery. We explore fast, automated optimization methods to solve this problem. We formulate the inertial flow physics in microchannels with different micropillar configurations as a set of state transition matrix operations. These state transition matrices are constructed from experimentally validated streamtraces. This facilitates modeling the effect of a sequence of micropillars as nested matrix-matrix products, which have very efficient numerical implementations. With this new forward model, arbitrary micropillar sequences can be rapidly simulated with various inlet configurations, allowing optimization routines quick access to a large search space. We integrate this framework with the genetic algorithm and showcase its applicability by designing micropillar sequences for various useful transformations. We computationally discover micropillar sequences for complex transformations that are substantially shorter than manually designed sequences. We also determine sequences for novel transformations that were difficult to manually design. Finally, we experimentally validate these computational designs by fabricating devices and comparing predictions with the results from confocal microscopy.


I. INTRODUCTION
The physics of inertial fluid flow deformation at the microscale has seen a surge of theoretical and experimental interest in the past decade 1 . Use of the inertial flow regime (1 < Re < 100, with the Reynolds number Re = ρV D H µ , where ρ, V , and µ are the fluid density, average downstream velocity, and viscosity, and D H the hydraulic diameter) in the microfluidics community contrasts the previously held notion that most practically useful flows were in the Stokes regime 2 (Re → 0), though the effects of inertial fluid flow has been observed since the 1960s 3 . New applications have since emerged, with inertial focusing in particular, leading to new methods for high-throughput cytometry [4][5][6][7][8][9][10] .
More recently, fluid sculpting via inertial flow has been demonstrated through so-called "pillar programming" 11,12 . See Fig 1 for a general schematic of pillar programming. In pillar programming, pillars spanning the height of a microchannel create a secondary flow which can be used in sequence with additional downstream pillars to generate a net deformation to the fluid. The idea of programmability comes from the micropillars being spaced far enough apart to allow for their individual deformations to saturate. Thus, the flow sculpted by one pillar can be used as an input to another pillar, without concern for cross-talk or time dependent effects. This allows for a powerful computational shortcut for future simulations, as the deformation for particular micropillar configurations need be computed only once, with the resulting fluid displacement forming a 2D "advection map". Net deformation for a sequence of micropillars can then be rapidly simulated using the pre-computed advection maps for each micropillar in the sequence. This subverts the need to solve Navier-Stokes equations over large 3D domains for entire microchannel designs using pillar programming.
Amini, et al determined numerically that an inter-pillar spacing of 6 pillar diameters is sufficient to prevent upstream flow effects of one pillar from interacting with a previous pillar's deformation 11 . In experiments, they use a spacing of 10 pillar diameters to ensure the pillar deformations are completely isolated. Therefore, each pillar in a designed microfluidic device will extend the length of the microchannel by a fixed amount, and its pillar programming simulation is valid only in the flow conditions for the pre-computed advection map.
With a coarse but broad design space for initial exploration, a set of widely varied transformations was found and experimentally validated using this method for forward simulation 12 .  2). For a 10-pillar sequence, there are 32 10 ≈ 10 15 possible pillar sequences, in addition to multiple inlet fluid flow configurations. It becomes nontrivial (in terms of time and effort) to manually design a micropillar sequence by searching this phase space. Even if manual design is attempted, a more efficient sequence with fewer pillars might be possible. This is important in microfluidic devices, as the number of pillars in a sequence determines the length of the channel, and therefore the footprint of the device itself. Longer channels and more pillars increase the pressure required for operation, which can introduce device flexure, and exacerbate transverse mass diffusion that increases with fluid distance traveled 19 . Decreasing the channel length can not only reduce the pressure drop to improve the performance of the device, but also provide more space downstream for applications, for example, complex 3D shaped particle fabrication or particle separation and solution exchange [13][14][15]17 . Therefore, an optimization scheme is needed to enhance the capability of pillar programming while still remaining accessible to, and easily implementable by, an interested researcher with modest computing hardware. Although pillar diameter will play a role in the required pressure to drive flow through a micropillar sequence, we have determined that microchannel length has a greater overall impact on pressure drop (more information available in supplemental material 20 ). Hence, the first priority in optimization is to find a flow shape that matches the target. This is the motivation for the work presented here. Secondary to this is the minimization of micropillar sequence length  11,12,18 . Each pillar spans the height of the microchannel. Note the periodic boundary condition, for which a pillar being placed so close to the microchannel wall will have the merged volume protrude from the opposite wall. For this paper, each pillar configuration of diameter D/w and offset y/w has an integer index, which is shown next to each depiction of the micropillar. We show two steps toward a framework for automated pillar sequence design: the development of an efficient forward model for rapid evaluation of arbitrary micropillar sequences, and the formulation of the design problem as an optimization problem. To accomplish these steps, we first convert the per-pillar advection maps to sparse transition matrices, which store fluid state displacement information as transition probabilities across states. These matrices can be multiplied together to quickly propagate net deformation across many pillars. This reduction of complex three-dimensional deformation simulations to numerically efficient matrix-matrix multiplication subsequently enables effectively pairing with a genetic algorithm for optimal design 21 of micropillar sequences. To demonstrate the effectiveness of this platform, we show the design and experimental validation of optimized micropillar sequences for previously discovered deformations, as well as novel designs which have no previously known pillar sequences.

II. PROBLEM FORMULATION
The goal is to develop a computational framework that accepts a desired transformation as an input, and produces a micropillar sequence design for the nearest possible match to this transformation as an output (see Fig. 3). For this work, the simulation and fabrications parameters will remain within the constraints of pillar programming as previously defined.

A. Forward model
The design problem potentially requires evaluation of many thousands of different pillar sequences. We refer to each such evaluation as solving the forward model (in contrast to the design problem). A forward model must therefore be accurate and fast. We begin by creating advection maps from streamtrace data for various pillar configurations, which is computed from our validated in-house finite element Computational Fluid Dynamics (CFD) framework 11,12,18 . We then translate the displacement information from the maps into transition matrices, which are an effective method for fast and accurate pillar program simulation.
1. CFD data to advection maps to transition matrices We start with a dataset of 3D velocity fields calculated for a set of individual micropillar configurations at Re = 20 (See Fig 2). The Navier-Stokes equations are solved (using the finite element method) in a domain that spans six diameters upstream and downstream of the pillar. A detailed description of this procedure is provided in our earlier work 11,12 and is not repeated here for the sake of brevity. Each 3D velocity field is contracted to a 2D advection map by streamtracing uniformly distributed, neutrally buoyant particles through the velocity field (see Fig. 4(a)).
For enhanced numerical efficiency, our idea is to make pillar programming amenable to using the BLAS (Basic Linear Algebra Subprograms) libraries within the CPU and/or GPGPU (General Purpose computation on Graphics Processing Units). BLAS libraries are freely available software that are architecture-aware, utilize memory optimally, and are Final particle cells Initial particle cells Initial particle cells   We convert the displacement information contained in an advection map into a sparse "transition matrix". This is accomplished by first discretizing the advection map into a finite set of N cells, and subsequently identifying where each cell advects to (or is displaced to).
Consider an advection map constructed by using N = N Y × N Z particles that are uniformly distributed across the cross-section (with N Y particles along the y direction, and N Z particles along the z direction). We discretize the channel cross-section into N = N Y × N Z cells, with each cell center displacement given by the corresponding particle displacement from the advection map. This is illustrated in Fig. 5 results in a sparse matrix that is easily stored and manipulated. Fig. 6 illustrates the sparsity of a representative transition matrix.
Once the transition matrices for individual pillar transformations have been computed, the net deformation caused by an arbitrary sequence of pillars is easily computed as the matrix product (in sequence) of the corresponding transition matrices. This is schematically shown in Fig. 7. In this example, the inlet flow shape is represented as a vector µ inlet . The fluid is transformed by a sequence of three pillars, which have transition matrices P 1 , P 2 , and P 3 , respectively. Then, the outlet flow shape µ outlet is given by Thus, the prediction of fluid flow shapes from arbitrary micropillar sequences has been reduced from computationally expensive CFD to simple sparse matrix multiplication 22 . The key advantages of the formalism are the speed with which flow shapes can be predicted (see Fig. 8(a)), and the subsequent simplicity of integrating this with a multitude of optimization frameworks. Furthermore, sparse matrices require far less memory than a dense matrix of similar size, since only the non-zero values are stored (along with their locations in the matrix). In our case, note that each cell is displaced to a single cell (a binary transition).
Thus, the framework optimally deploys sparsity with the number of non-zero values in each row equal to 1. This reduces computational overhead. In the next section, we explore the accuracy of this transition matrix representation as a function of discretization.

B. Implementation
Consider a transition matrix constructed from a discretization of the microchannel crosssection into N Y × N Z uniform cells. The cells have dimensions w N Y × h N Z , but are effectively represented as points in the transition matrix. As such, subtle behavior of streamtraces will be truncated to displacements that align with the discretized representation. For example, fluid movements too small to leave their immediate cell space will not register in the transition matrix (see cell 1 in Fig. 5(a)). Similarly, movement just large enough to arrive in a new cell will be considered as completely displaced to this cell center (see cell 2 in Fig 5(a).
The precision of the transition matrix is therefore limited by the level of discretization of the cross-section into cells. Increasing the number of cells (larger N Y , N Z ) accounts for smaller displacements, thereby leading to a more accurate mapping onto the transition matrix (see  We next explore the representation error that these discretizations produce, especially as a function of nested matrix-matrix products. Note that experimental and mass diffusion constraints (please see Appendix A) limit the maximum number of pillars in a sequence to 10. We evaluate the representation error by comparing the predictions of nested matrix products with those of the advection maps for arbitrary pillar sequences of increasing complexity. Error is computed by first converting the transition matrices back into 2D advection maps, and comparing these to the streamtraced advection maps. Thus, for an arbitrary streamtraced advection map A and matrix-reconstructed map A P : We create many sequences of pillars and compute the error. These sequences are constructed using a quasi-random sobol number generator (which is a low discrepancy generator) for maximal coverage of the phase space of pillar combinations. The mean error over 100 real- izations is plotted in Fig. 8(b). From these plots, we chose to utilize a N Y = 801, N Z = 101 discretization for subsequent analysis, as it provides a good balance between accuracy and computational overhead. We finally validate this discretization by comparing the prediction of the transition matrix framework with experimental confocal images (see Section IV B,C) of three flow transformations. Fig. 9 illustrates this promising comparison.

III. DESIGN PROBLEM
Having established an efficient pathway for quickly evaluating arbitrary sequences of micropillars, we next formulate the problem of identifying micropillar sequences that result in user defined transformations.
We formulate the design problem as an optimization problem. Specifically, given a target outlet shape I T arget , we define a cost functional C(I T arget , I T est (s)) that returns a scalar value representing how closely a test shape, I T est (s), matches the target shape.
where the pillar sequence s that minimizes the fitness function is found.
There are several approaches to solve this optimization problem. Here, we choose to utilize a gradient free, evolutionary optimization strategy. The rationale behind this choice is motivated by the following: 1. The problem formulation as defined here is inherently discrete. This is due to the underlying experimental constraints. Fabrication of arbitrarily sized pillars at continuous locations within practical tolerance and cost is not viable. Effective manufacturing tolerances benefit from a well defined set of micropillar configurations, rather than a continuous space of micropillar diameters and offsets. Simulating a significant subset of this infinite library is also not practical. Nevertheless, the discrete set from which each pillar is chosen makes the possible search space countably large, hence, precluding an exhaustive search.
2. Though it is difficult to illustrate the phase space a 10-pillar sequence offers, a cost function (developed later in this work) evaluated for 2-pillar sequences (with each pillar chosen from 32 possible configurations as shown in Fig. 2) is shown in Fig. 10. This is a highly corrugated surface, with one global optimum but many nearly identical local minima. This precludes the utilization of gradient based methods, and instead suggests the applicability of stochastic, multistart methods that can explore the phase space efficiently.
3. The discrete problem of finding a sequence of matrices such that their product matrix has a desired structure is a variant of the minimum length generator sequence problem 23 which has been shown to be a NP-hard problem. This necessitates the utilization of (meta)heuristic optimization methods for efficient solutions.
These issues naturally suggested using gradient-free metaheuristic evolutionary search algorithms. We specifically use the genetic algorithm (GA) to locate optimal sequences. GAs are well suited to multi-modal, highly corrugated solution spaces, especially when the cost function is not easily adapted to gradient-based methods [24][25][26] .

A. Genetic Algorithm
The GA has been used with substantial success in fluid mechanics with applications ranging from efficient macro-to nano-fluid heat exchanger design [27][28][29] , bluff body flow control 30 , indoor airflow geometry 31 , and aircraft design 32 . Though there are many different flavors and modifications of the GA, the basic algorithm is generally as follows: 1. Initialize a population of randomly generated points in the search space. This is the first generation of points. The points are encoded (usually as a binary bitstring) and are called "chromosomes".
2. Evaluate each chromosome of the current generation using a fitness function (cost function) specific to the problem. 3. Select chromosomes based on their fitness to create a new generation of individuals through crossover and mutation methods. Good fitness is rewarded by increased chance of selection, which will propagate genetic material to the next generation.
There are three basic methods of creating the next generation: • Crossover combines the genetic material in individuals' chromosomes (design parameters) to create offspring that should retain some of the traits of the parents.
A crossover rate determines the proportion of each new generation that should come from crossover.
• Mutation randomly alters the bits within a chromosome. This introduces random variation into the population, which can bring candidate solutions out of local optima.
• A set of top-performing chromosomes in the population are chosen as Elites, whose chromosomes are preserved in the next generation. This retains the optimal solution throughout the GA. 4. Repeat steps 2-3 until some set of termination criteria are met, which typically include the following: • Stall generation limit: If mean or optimal fitness does not improve for a set number of generations, terminate the GA. During the evaluation stage of the GA, the previously described forward model creates a fluid flow shape for each chromosome, which is then compared to the target image in the fitness function.
The choice of the fitness function critically determines the success of the optimization procedure. Since the primary goal of this exercise is to capture the overall shape of the transformation, the fitness function should be a measure of the topology of the shape.
Extensive numerical experiments suggested that the (image) correlation coefficient (defined below) gives substantially better results than several standard pixel based norms, while remaining competitive with respect to computational speed. 33 The correlation coefficient r, is defined as This measure is extensively used to determine how similar two images are. Identical images would result in r = 1, while comparing images with precisely complementary pixels would result in r = −1. In the latter case, the result is still desirable, as the overall shape has been achieved using fluid in complementary inlets. That is, the "empty" co-flow inlets have been shaped to align with the target fluid locations in the desired fluid flow shape. The user would simply need to "flip" the inlet flow configuration in order to create their desired microfluidic device. We can treat both optimal values of r equally by squaring its value.
Thus, the fitness function C takes the form We achieved better convergence by pre-processing I T arget and I test (s) before evaluating the fitness function. Specifically, a low-pass filter applied to both I T arget , and I test (s) emphasized the large scale (topological) features of the shapes while de-emphasizing the fine scale features. This is shown in Fig. 11. When used with the correlation function, images similar in shape will generally have a better cost functional than without the low-pass filter.

A. Improved efficiency of existing designs
Three of the more complex, hierarchical designs from our earlier work showing manual designs 12 were selected to demonstrate the potential to create more efficient pillar sequences for existing transformations. These three transformations are the encapsulate, shift, and add vertex transformations (see Fig. 13(a-c,i-ii)). All three deform the same inlet flow configuration shown in Fig. 7 GA optimization found multiple designs for each target transformation, with the shortest sequences and simulated flow shapes shown in Fig. 13(e-f,i-ii). See Fig. 2 for pillar configurations corresponding to the indices in the sequence figures. The new designs resulted in 33-60% improvement to sequence length while maintaining the overall desired shape and location in the microchannel. These designs were fabricated and the pillar transformations were evaluated using confocal images. The experimental results agree well with the predictions as seen in Fig. 13(d-f,iii). The improvement by computational optimization also reduced device footprint and pressure by approximately the same percentage in experiments.
All simulations were performed using a Matlab implementation of the sparse matrix operations. We deployed the parallel GA routine available in Matlab on computing clusters.
GA parameters included a population size of 100, crossover rate of 0.8, mutation rate of 0.2, and 5 elites. Example behavior of the algorithm is shown in Fig. 12, showing optimum fitness functions and runtime for the optimization of the encapsulate flow shape. We see good support for using GAs by comparing their evolved generation's fitness functions to a computationally equivalent random searches. The GA will evolve 100 initially random pillar sequences for approximately 100 generations (with an upper limit of 200), so we evaluated 10,000 randomly generated pillar sequences. Fig 12(a) shows fits of normal distributions to the results of these searches, with the GA results having clear improvement over the random search. The final generation for the GA found an optimum fitness of C = 6.76, while the random search's best fitness was C = 21.56. The GA typically converges on the same optimal shape for smaller design spaces (see Fig 12(b), for # of pillars = 1-3). Searches using pillar sequences with a length in excess of 5 pillars find a variety of local optima, which is unsurprising given the larger, more complex search space. Runtime was as much as 30 minutes, but on average the total run of 10 trials for 10 pillar sequence lengths took about 24 hours on a 2.0 GHz 8-Core Intel E5 2650 CPU.

B. Novel designs
We next applied the framework to design novel flow transformations. These transformations change the shape of the central fluid into shapes akin to a dumbbell. The dumbbell shape was motivated by experimental particle fabrication research, where polymer precur- 29, 29, 29, 13, 9, 29, 21 13, 21, 29, 21, 29, 29, 12, 14, 6, 12 FIG. 14. Optimization targeting novel "dumbbell" shapes (a,i) and (b,i), with genetic algorithm solution sequences (ii) and post-processed flow predictions (iii). For shape (b,iii), the inlet flow shape spans a width of w 3 . Pillar indices can be found in Fig. 2. sors are shaped by a micropillar sequence and subsequently polymerized by UV light. Use of a shaped mask allows for a single section of flow to be polymerized, rather than the entire stream as in 14 . Unlike circular or ellipsoid shapes, the dumbbell offers a more complex geometry that could be applied to create particles that align along one axis in a flow 34 .

Novel sequences found via genetic algorithm
The targets used here differ with respect to the width of the dumbbell, which will show the sensitivity in the choice of a desired shape.
The GA was run with several different inlet configurations in order to allow for additional fluid. Results are in Fig. 14, with the wider dumbbell shape requiring an inlet flow shape spanning w 3 .

VI. SUMMARY
We have devised and validated an efficient forward model for simulating pillar programming via simple matrix multiplication, and implemented the model into a parallel GA. The simple nature of the new forward model is valuable for its speed and applicability. Optimization routines can now incorporate flow shape prediction directly into a computational pipeline. We have shown how the GA successfully used this new framework by optimizing for known flow shapes, and producing new micropillar sequence designs for novel flow shapes. Future work on design for pillar programming could pursue faster and more exhaustive searches by tailoring the process to more specialized hardware (e.g., the GPU), and investigating new fitness functions.
Overall, this work completes the cycle moving from a user-defined flow shape of interest to a physical implementation of a channel design that achieves this flow shape. Such an approach opens up the computer-aided design and manufacturing of shaped polymer fibers 13 and particles 14 for a range of applications. Additional uses in directing mass and heat transfer, and transferring solutions for automation of biological sample preparation should also benefit. In comparing to experimental and numerical iteration in which the full Navier-Stokes equations are solved for a set of complex channel geometries, our approach achieves orders of magnitudes improvements in time to result, making previously intractable problems that were not even attempted now possible.

VII. ACKNOWLEDGEMENTS
This research is supported in part by the National Science Foundation through NSF-1306866, NSF-1307550, and NSF-1149365. 33 Considering the importance of cost function selection with regards to optimization algorithms, we emphasize that it will be valuable to investigate additional cost functions (specifically wavelet and Fourier transform based multiscale feature functions).
34 W. E. Uspal, H. Burak Eral, and P. S. Doyle, "Engineering particle trajectories in microfluidic flows using particle shape," Nature communications 4, 2666 (2013). microchannel width w = 200µm. Pillar spacing is critical, as the deformation from one pillar must saturate before before fluid arrives at the subsequent pillar. Pillars in close proximity to each other (L < 6D P ) will encounter "cross-talk" in their respective fluid deformation, and the general premise of pillar programming will break down. Local fluid velocities in the 3D domain will vary depending geometry of the channel, but we can approximate the added time of flight by multiplying the distance travelled by a compensating coefficient f : Or, we can define the diffusive length using the Péclet number, P e = wU D : δ = wnf L P e