@unpublished{shahane24,
  author = {Shahane, S. and Chammas, S. and Bezgin, D. A. and Buhendwa, A. B. and Schmidt, S. J. and Adams, N. A. and Bryngelson, S. H. and Chen, Y.-F. and Wang, Q. and Sha, F. and Zepeda-Núñez, L.},
  title = {{Rational-WENO: A} lightweight, physically-consistent three-point weighted essentially non-oscillatory scheme},
  note = {arXiv: 2409.09217},
  doi = {10.48550/arXiv.2409.09217},
  file = {shahane-arxiv-24.pdf},
  arxiv = {arXiv.2409.09217},
  year = {2024}
}

Conventional WENO3 methods are known to be highly dissipative at lower resolutions, introducing significant errors in the pre-asymptotic regime. In this paper, we employ a rational neural network to accurately estimate the local smoothness of the solution, dynamically adapting the stencil weights based on local solution features. As rational neural networks can represent fast transitions between smooth and sharp regimes, this approach achieves a granular reconstruction with significantly reduced dissipation, improving the accuracy of the simulation. The network is trained offline on a carefully chosen dataset of analytical functions, bypassing the need for differentiable solvers. We also propose a robust model selection criterion based on estimates of the interpolation’s convergence order on a set of test functions, which correlates better with the model performance in downstream tasks. We demonstrate the effectiveness of our approach on several one-, two-, and three-dimensional fluid flow problems: our scheme generalizes across grid resolutions while handling smooth and discontinuous solutions. In most cases, our rational network-based scheme achieves higher accuracy than conventional WENO3 with the same stencil size, and in a few of them, it achieves accuracy comparable to WENO5, which uses a larger stencil.

@unpublished{chu24,
  author = {Chu, T. and Estrada, J. B. and Bryngelson, S. H.},
  title = {Bayesian optimal design accelerates discovery of material properties from bubble dynamics},
  doi = {10.48550/arXiv.2409.00011},
  note = {arXiv: 2409.00011},
  file = {chu-arxiv-24.pdf},
  arxiv = {2409.00011},
  year = {2024}
}

An optimal sequential experimental design approach is developed to characterize soft material properties at the high strain rates associated with bubble cavitation. The approach involves optimal design and model inference. The optimal design strategy maximizes the expected information gain in a Bayesian statistical setting to design experiments that provide the most informative cavitation data about unknown material properties. We infer constitutive models by characterizing the associated viscoelastic properties from measurements via a hybrid ensemble-based 4D-Var method (En4D-Var). The inertial microcavitation-based high strain-rate rheometry (IMR) method (Estrada et al. 2018) is used to simulate the bubble dynamics under laser-induced cavitation. Two sets of constitutive models demonstrate the method. We generate synthetic experimental data representing the viscoelastic behavior of stiff and soft polyacrylamide hydrogels. Accurate and efficient characterizations of the underlying models are presented.

@unpublished{kocherla24_2,
  author = {Kocherla, S. and Adams, A. and Song, Z. and Alexeev, A. and Bryngelson, S. H.},
  title = {A two-circuit approach to reducing quantum resources for the quantum lattice {B}oltzmann method},
  note = {arXiv: 2401.12248},
  file = {kocherla-arxiv-24.pdf},
  arxiv = {2401.12248},
  doi = {10.48550/arXiv.2401.12248},
  year = {2024}
}

Computational fluid dynamics (CFD) simulations often entail a large computational burden on classical computers. At present, these simulations can require up to trillions of grid points and millions of time steps. To reduce costs, novel architectures like quantum computers may be intrinsically more efficient at the appropriate computation. Current quantum algorithms for solving CFD problems use a single quantum circuit and, in some cases, lattice-based methods. We introduce the a novel multiple circuits algorithm that makes use of a quantum lattice Boltzmann method (QLBM). The two-circuit algorithm we form solves the Navier-Stokes equations with a marked reduction in CNOT gates compared to existing QLBM circuits. The problem is cast as a stream function-vorticity formulation of the 2D Navier-Stokes equations and verified and tested on a 2D lid-driven cavity flow. We show that using separate circuits for the stream function and vorticity lead to a marked CNOT reduction: 35% in total CNOT count and 16% in combined gate depth. This strategy has the additional benefit of the circuits being able to run concurrently, further halving the seen gate depth. This work is intended as a step towards practical quantum circuits for solving differential equation-based problems of scientific interest.

@article{song24,
  author = {Song, Z. and Deaton, R. and Gard, B. and Bryngelson, S. H.},
  title = {Incompressible {N}avier--{S}tokes solve on noisy quantum hardware via a hybrid quantum--classical scheme},
  journal = {Computers \& Fluids},
  pages = {106507},
  file = {song-caf-25.pdf},
  doi = {10.1016/j.compfluid.2024.106507},
  year = {2025}
}

Partial differential equation solvers are required to solve the Navier–Stokes equations for fluid flow. Recently, algorithms have been proposed to simulate fluid dynamics on quantum computers. Fault-tolerant quantum devices might enable exponential speedups over algorithms on classical computers. However, current and foreseeable quantum hardware introduce noise into computations, requiring algorithms that make judicious use of quantum resources: shallower circuit depths and fewer qubits. Under these restrictions, variational algorithms are more appropriate and robust. This work presents a hybrid quantum–classical algorithm for the incompressible Navier–Stokes equations. A classical device performs nonlinear computations, and a quantum one uses a variational solver for the pressure Poisson equation. A lid-driven cavity problem benchmarks the method. We verify the algorithm via noise-free simulation and test it on noisy IBM superconducting quantum hardware. Results show that high-fidelity results can be achieved via this approach, even on current quantum devices. Multigrid preconditioning of the Poisson problem helps avoid local minima and reduces resource requirements for the quantum device. A quantum state readout technique called HTree is used for the first time on a physical problem. HTree is appropriate for real-valued problems and achieves linear complexity in the qubit count, making the Navier–Stokes solve further tractable on current quantum devices. We compare the quantum resources required for near-term and fault-tolerant solvers to determine quantum hardware requirements for fluid simulations with complexity improvements.

@article{radhakrishnan24,
  author = {Radhakrishnan, A. and {Le Berre}, H. and Wilfong, B. and Spratt, J.-S. and {Rodriguez Jr.}, M. and Colonius, T. and Bryngelson, S. H.},
  title = {Method for portable, scalable, and performant {GPU}-accelerated simulation of multiphase compressible flow},
  file = {radhakrishnan-CPC-24.pdf},
  year = {2024},
  volume = {302},
  doi = {10.1016/j.cpc.2024.109238},
  journal = {Computer Physics Communications},
  pages = {109238}
}

Multiphase compressible flows are often characterized by a broad range of space and time scales. Thus entailing large grids and small time steps, simulations of these flows on CPU-based clusters can thus take several wall-clock days. Offloading the compute kernels to GPUs appears attractive but is memory-bound for standard finite-volume and -difference methods, damping speed-ups. Even when realized, faster GPU-based kernels lead to more intrusive communication and I/O times. We present a portable strategy for GPU acceleration of multiphase compressible flow solvers that addresses these challenges and obtains large speedups at scale. We use OpenACC for portable offloading of all compute kernels while maintaining low-level control when needed. An established Fortran preprocessor and metaprogramming tool, Fypp, enables otherwise hidden compile-time optimizations. This strategy exposes compile-time optimizations and high memory reuse while retaining readable, maintainable, and compact code. Remote direct memory access, realized via CUDA-aware MPI, reduces communication times. We implement this approach in the open-source solver MFC. Metaprogramming-based preprocessing results in an 8-times speedup of the most expensive kernels, 46% of peak FLOPs on NVIDIA GPUs, and high arithmetic intensity (about 10 FLOPs/byte). In representative simulations, a single A100 GPU is 300-times faster than an Intel Xeon CPU core, corresponding to a 9-times speedup for a single A100 compared to the entire CPU die. At the same time, near-ideal (97%) weak scaling is observed for at least 13824 GPUs on Summit. A strong scaling efficiency of 84% is retained for an 8-times increase in GPU count. Collective I/O, implemented via MPI3, helps ensure negligible contribution of data transfers. Large many-GPU simulations of compressible (solid-)liquid-gas flows demonstrate the practical utility of this strategy.

@article{sinha24neural,
  title = {Neural networks can be {FLOP}-efficient integrators of {1D} oscillatory integrands},
  author = {Sinha, A. and Bryngelson, S. H.},
  journal = {Transactions on Machine Learning Research},
  issn = {2835-8856},
  year = {2024},
  file = {sinha-TMLR-24.pdf},
  link = {https://openreview.net/pdf?id=5psgQEHn6t}
}

We demonstrate that neural networks can be FLOP-efficient integrators of one-dimensional oscillatory integrands. We train a feed-forward neural network to compute integrals of highly oscillatory 1D functions. The training set is a parametric combination of functions with varying characters and oscillatory behavior degrees. Numerical examples show that these networks are FLOP-efficient for sufficiently oscillatory integrands with an average FLOP gain of 1000 FLOPs. The network calculates oscillatory integrals better than traditional quadrature methods under the same computational budget or number of floating point operations. We find that feed-forward networks of 5 hidden layers are satisfactory for a relative accuracy of 0.001. The computational burden of inference of the neural network is relatively small, even compared to inner-product pattern quadrature rules. We postulate that our result follows from learning latent patterns in the oscillatory integrands that are otherwise opaque to traditional numerical integrators.

@article{bryngelsonfmfm24,
  author = {{*}Bryngelson, S. H. and {*}Sch{\"a}fer, F. and Liu, J. and Mani, A.},
  title = {Fast {M}acroscopic {F}orcing {M}ethod},
  journal = {Journal of Computational Physics},
  file = {bryngelson-JCP-24.pdf},
  year = {2024},
  volume = {499},
  pages = {112721},
  doi = {10.1016/j.jcp.2023.112721},
  note = {{*}Equal contribution}
}

The macroscopic forcing method (MFM) of Mani and Park (2021) and similar methods for obtaining turbulence closure operators, such as the Green’s function-based approach of Hamba (1995), recover reduced solution operators from repeated direct numerical simulations (DNS).  MFM has already been used to successfully quantify Reynolds-averaged Navier-Stokes (RANS)-like operators for homogeneous isotropic turbulence and turbulent channel flows. Standard algorithms for MFM force each coarse-scale degree of freedom (i.e., degree of freedom in the RANS space) and conduct a corresponding fine-scale simulation (i.e., DNS), which is expensive.  We combine this method with an approach recently proposed by Schäfer and Owhadi (2023) to recover elliptic integral operators from a polylogarithmic number of matrix-vector products. The resulting Fast MFM introduced in this work applies sparse reconstruction to expose local features in the closure operator and reconstructs this coarse-grained differential operator in only a few matrix-vector products and correspondingly, a few MFM simulations.  For flows with significant nonlocality, the algorithm first “peels” long-range effects with dense matrix-vector products to expose a more local operator.  We demonstrate the algorithm’s performance for scalar transport in a laminar channel flow and momentum transport in a turbulent channel flow. For these problems, we recover eddy diffusivity- and eddy viscosity-like operators, respectively, at 1 percent of the cost of computing the exact operator via a brute-force approach for the laminar channel flow problem and 13 percent for the turbulent one. We observe that we can reconstruct these operators with an increase in accuracy by about a factor of 100 over randomized low-rank methods.  Applying these operators to compute the averaged fields of interest has visually indistinguishable behavior from the exact solution.  Our results show that a similar number of simulations are required to reconstruct the operators to the same accuracy under grid refinement.  Thus, the accuracy corresponds to the physics of the problem, not the numerics. We glean that for problems in which the RANS space is reducible to one dimension, eddy diffusivity and eddy viscosity operators can be reconstructed with reasonable accuracy using only a few simulations, regardless of simulation resolution or degrees of freedom.

@article{bati24,
  author = {Bati, A. and Bryngelson, S. H.},
  title = {{RoseNNa: A} performant, portable library for neural network inference with
  application to computational fluid dynamics},
  journal = {Computer Physics Communications},
  volume = {296},
  pages = {109052},
  file = {bati-CPC-24.pdf},
  year = {2024},
  doi = {10.1016/j.cpc.2023.109052}
}

The rise of neural network-based machine learning ushered in high-level libraries, including TensorFlow and PyTorch, to support their functionality. Computational fluid dynamics (CFD) researchers have benefited from this trend and produced powerful neural networks that promise shorter simulation times. For example, multilayer perceptrons (MLPs) and Long Short Term Memory (LSTM) recurrent-based (RNN) architectures can represent sub-grid physical effects, like turbulence. Implementing neural networks in CFD solvers is challenging because the programming languages used for machine learning and CFD are mostly non-overlapping, We present the roseNNa library, which bridges the gap between neural network inference and CFD. RoseNNa is a non-invasive, lightweight (1000 lines), and performant tool for neural network inference, with focus on the smaller networks used to augment PDE solvers, like those of CFD, which are typically written in C/C++ or Fortran. RoseNNa accomplishes this by automatically converting trained models from typical neural network training packages into a high-performance Fortran library with C and Fortran APIs. This reduces the effort needed to access trained neural networks and maintains performance in the PDE solvers that CFD researchers build and rely upon. Results show that RoseNNa reliably outperforms PyTorch (Python) and libtorch (C++) on MLPs and LSTM RNNs with less than 100 hidden layers and 100 neurons per layer, even after removing the overhead cost of API calls. Speedups range from a factor of about 10 and 2 faster than these established libraries for the smaller and larger ends of the neural network size ranges tested.

@article{kocherla24_1,
  author = {Kocherla, S. and Song, Z. and Chrit, F. E. and Gard, B. and Dumitrescu, E. F. and Alexeev, A. and Bryngelson, S. H.},
  title = {Fully quantum algorithm for mesoscale fluid simulations
  with application to partial differential equations},
  doi = {10.1116/5.0217675},
  file = {kocherla-AVS-24.pdf},
  year = {2024},
  volume = {6},
  pages = {033806},
  journal = {AVS Quantum Science}
}

Fluid flow simulations marshal our most powerful computational resources. In many cases, even this is not enough. Quantum computers provide an opportunity to speed up traditional algorithms for flow simulations. We show that lattice-based mesoscale numerical methods can be executed as efficient quantum algorithms due to their statistical features. This approach revises a quantum algorithm for lattice gas automata to reduce classical computations and state preparation at every time step. For this, the algorithm approximates the qubit relative phases and subtracts them at the end of each time step. Phases are evaluated using the iterative phase estimation algorithm and subtracted using single-qubit rotation phase gates. This method optimizes the quantum resource required and makes it more appropriate for near-term quantum hardware. We also demonstrate how the checkerboard deficiency that the D1Q2 scheme presents can be resolved using the D1Q3 scheme. The algorithm is validated by simulating two canonical partial differential equations: the diffusion and Burgers’ equations on differ

@article{liu24,
  author = {Liu, J. and Sch\"afer, F. and Bryngelson, S. H. and Zaki, T. A. and Mani, A.},
  title = {Adjoint-based computation of nonlocal eddy viscosity in turbulent channel flow},
  year = {2024},
  journal = {Physical Review Fluids},
  volume = {9},
  pages = {094606},
  doi = {10.1103/PhysRevFluids.9.094606},
  file = {liu-PRF-24.pdf}
}

Reynolds-averaged Navier–Stokes (RANS) closure operators are generally nonlocal and anisotropic for many flows, including wall-bounded turbulence. The macroscopic forcing method (MFM) and Green’s function-based approaches examine these effects using forced direct numerical simulations. For these approaches, the number of simulations required to compute the nonlocal and anisotropic eddy viscosity is equal to the number of degrees of freedom in the averaged space. We reduce the computational cost by introducing an adjoint-based formulation of MFM to obtain the nonlocal and anisotropic eddy viscosity at any Reynolds stress location using a single simulation. We then quantify the streamwise and wall-normal nonlocal eddy viscosity at near-wall locations in turbulent channel flow at shear Reynolds number 180. We demonstrate that the upstream nonlocality is not fully described by Lagrangian transport. We also quantify the significant differences in the upstream nonlocality between various components of the eddy viscosity tensor. Our results can be used to guide closure modeling, including the lengthscales and timescales used in Reynolds stress models.

@inproceedings{proceeding_25,
  title = {{OpenACC} offloading of the {MFC} compressible multiphase flow solver on {AMD} and {NVIDIA GPUs}},
  author = {Wilfong, B. and Radhakrishnan, A. and {Le Berre}, H. A. and Abbott, S. and Budiardja, R. D. and Bryngelson, S. H.},
  booktitle = {Proceedings of the SC '24 Workshops of The International Conference on High Performance Computing, Network, Storage, and Analysis},
  year = {2024},
  doi = {10.48550/arXiv.2409.10729},
  file = {wilfong-SC-24.pdf},
  keywords = {heavyref}
}

GPUs are the heart of the latest generations of supercomputers. We efficiently accelerate a compressible multiphase flow solver via OpenACC on NVIDIA and AMD Instinct GPUs. Optimization is accomplished by specifying the directive clauses ’gang vector’ and ’collapse’. Further speedups of six and ten times are achieved by packing user-defined types into coalesced multidimensional arrays and manual inlining via metaprogramming. Additional optimizations yield seven-times speedup in array packing and thirty-times speedup of select kernels on Frontier. Weak scaling efficiencies of 97% and 95% are observed when scaling to 50% of Summit and 95% of Frontier. Strong scaling efficiencies of 84% and 81% are observed when increasing the device count by a factor of 8 and 16 on V100 and MI250X hardware. The strong scaling efficiency of AMD’s MI250X increases to 92% when increasing the device count by a factor of 16 when GPU-aware MPI is used for communication.

@article{panchal23,
  author = {Panchal, A. and Bryngelson, S. H. and Menon, S.},
  title = {A seven-equation diffused interface method for resolved multiphase flows},
  journal = {Journal of Computational Physics},
  file = {panchal-JCP-23.pdf},
  volume = {475},
  pages = {111870},
  year = {2023},
  doi = {10.1016/j.jcp.2022.111870}
}

The seven-equation model is a compressible multiphase formulation that allows for phasic velocity and pressure disequilibrium. These equations are solved using a diffused interface method that models resolved multiphase flows. Novel extensions are proposed for including the effects of surface tension, viscosity, multi-species, and reactions. The allowed non-equilibrium of pressure in the seven-equation model provides numerical stability in strong shocks and allows for arbitrary and independent equations of states. A discrete equations method (DEM) models the fluxes. We show that even though stiff pressure- and velocity-relaxation solvers have been used, they are not needed for the DEM because the non-conservative fluxes are accurately modeled. An interface compression scheme controls the numerical diffusion of the interface, and its effects on the solution are discussed. Test cases are used to validate the computational method and demonstrate its applicability. They include multiphase shock tubes, shock propagation through a material interface, a surface-tension-driven oscillating droplet, an accelerating droplet in a viscous medium, and shock-detonation interacting with a deforming droplet. Simulation results are compared against exact solutions and experiments when possible.

@article{firouznia23,
  author = {Firouznia, M. and Bryngelson, S. H. and Saintillan, D.},
  title = {A spectral boundary integral method for simulating electrohydrodynamic flows in viscous drops},
  year = {2023},
  pages = {112248},
  volume = {489},
  file = {firouznia-JCP-23.pdf},
  journal = {Journal of Computational Physics},
  doi = {10.1016/j.jcp.2023.112248}
}

A weakly conducting liquid droplet immersed in another leaky dielectric liquid can exhibit rich dynamical behaviors under the effect of an applied electric field. Depending on material properties and field strength, the nonlinear coupling of interfacial charge transport and fluid flow can trigger electrohydrodynamic instabilities that lead to shape deformations and complex dynamics. We present a spectral boundary integral method to simulate droplet electrohydrodynamics in a uniform electric field. All physical variables, such as drop shape and interfacial charge density, are represented using spherical harmonic expansions. In addition to its exponential accuracy, the spectral representation affords a nondissipative dealiasing method required for numerical stability. A comprehensive charge transport model, valid under a wide range of electric field strengths, accounts for charge relaxation, Ohmic conduction, and surface charge convection by the flow. A shape reparametrization technique enables the exploration of significant droplet deformation regimes. For low-viscosity drops, the convection by the flow drives steep interfacial charge gradients near the drop equator. This introduces numerical ringing artifacts we treat via a weighted spherical harmonic expansion, resulting in solution convergence. The method and simulations are validated against experimental data and analytical predictions in the axisymmetric Taylor and Quincke electrorotation regimes.

@inproceedings{proceeding_15,
  author = {Elwasif, W. and Bastrakov, S. and Bryngelson, S. H. and Bussmann, M. and Chandrasekaran, S. and Ciorba, F. and Clark, M. A. and Debus, A. and Godoy, W. and Hagerty, N. and Hammond, J. and Hardy, D. and Harris, J. A. and Hernandez, O. and Joo, B. and Keller, S. and Kent, P. and {Le Berre}, H. and Lebrun-Grandie, D. and MacCarthy, E. and Vergara, V. G. Melesse and Messer, B. and Miller, R. and Oral, S. and Piccinali, J.-G. and Radhakrishnan, A. and Simsek, O. and Spiga, F. and Steiniger, K. and Stephan, J. and Stone, J. E. and Trott, C. and Widera, R. and Young, J.},
  title = {Early application experiences on a modern {GPU}-accelerated {A}rm-based {HPC} platform},
  booktitle = {HPC Asia '23},
  series = {International Workshop on Arm-based HPC: Practice and Experience (IWAHPCE)},
  file = {elwasif-IWAHPCE-23.pdf},
  doi = {10.1145/3581576.3581621},
  address = {Singapore},
  year = {2023},
  acceptrate = {44},
  keywords = {heavyref}
}

This paper assesses and reports the experience of eleven application teams working to build, validate, and benchmark several HPC applications on a novel GPU-accerated Arm testbed. The testbed consists of the latest, at time of writing, Arm Devkits from NVIDIA with server-class Arm CPUs and NVIDIA A100 GPUs. The applications and mini-apps are written using multiple parallel programming models, including C, CUDA, Fortran, OpenACC, and OpenMP. Each application builds extensively on the other tools available in the programming environment, including scientific libraries, compilers, and other tooling. Our goal is to evaluate application readiness for the next generation of Armand GPU-based HPC systems and determine the tooling readiness for future application developers. On both accounts, the reported case studies demonstrate that the diversity of software and tools available for GPU-accelerated Arm systems are prepared for production, even before NVIDIA deploys their next-generation such platform: Grace.

@article{bryngelson23,
  author = {Bryngelson, S. H. and Fox, R. O. and Colonius, T.},
  title = {Conditional moment methods for polydisperse cavitating flows},
  file = {bryngelson-JCP-23.pdf},
  journal = {Journal of Computational Physics},
  year = {2023},
  volume = {477},
  doi = {10.1016/j.jcp.2023.111917},
  pages = {111917}
}

The dynamics of cavitation bubbles are important in many flows, but their small sizes and high number densities often preclude direct numerical simulation. We present a computational model that averages their effect on the flow over larger spatiotemporal scales. The model is based on solving a generalized population balance equation (PBE) for nonlinear bubble dynamics and explicitly represents the evolving probability density of bubble radii and radial velocities. Conditional quadrature-based moment methods (QBMMs) are adapted to solve this PBE. A one-way-coupled bubble dynamics problem demonstrates the efficacy of different QBMMs for the evolving bubble statistics. Results show that enforcing hyperbolicity during moment inversion (CHyQMOM) provides comparable model-form accuracy to the traditional conditional method of moments and decreases computational costs by about ten times for a broad range of test cases. The CHyQMOM-based computational model is implemented in MFC, an open-source multi-phase and high-order-accurate flow solver. We assess the effect of the model and its parameters on a two-way coupled bubble screen flow problem.

@inproceedings{proceeding_14,
  author = {Zeng, Q. and Kothari, Y. and Bryngelson, S. H. and Sch{\"a}fer, F.},
  title = {Competitive physics informed networks},
  booktitle = {International Conference on Learning Representations (ICLR)},
  file = {zeng-ICLR-23.pdf},
  note = {arXiv:2204.11144},
  address = {Kigali, Rwanda},
  year = {2023},
  link = {https://openreview.net/pdf?id=z9SIj-IM7tn},
  keywords = {heavyref}
}

Physics Informed Neural Networks (PINNs) solve partial differential equations (PDEs) by representing them as neural networks. The original PINN implementation does not provide high accuracy, typically attaining about 0.1 percent relative error. We formulate and test an adversarial approach called competitive PINNs (CPINNs) to overcome this limitation. CPINNs train a discriminator that is rewarded for predicting PINN mistakes. The discriminator and PINN participate in a zero-sum game with the exact PDE solution as an optimal strategy. This approach avoids the issue of squaring the large condition numbers of PDE discretizations. Numerical experiments show that a CPINN trained with competitive gradient descent can achieve errors two orders of magnitude smaller than that of a PINN trained with Adam or stochastic gradient descent.

@article{charalampopoulos21,
  author = {Charalampopoulos, A. and Bryngelson, S. H. and Colonius, T. and Sapsis, T. P.},
  title = {Hybrid quadrature moment method for accurate and stable representation of non-{G}aussian processes and their dynamics},
  file = {charalampopoulos-RSA-21.pdf},
  journal = {Philosophical Transactions of the Royal Society A},
  year = {2022},
  volume = {380},
  number = {2229},
  doi = {10.1098/rsta.2021.0209}
}

Solving the population balance equation (PBE) for the dynamics of a dispersed phase coupled to a continuous fluid is expensive. Still, one can reduce the cost by representing the evolving particle density function in terms of its moments. In particular, quadrature-based moment methods (QBMMs) invert these moments with a quadrature rule, approximating the required statistics. QBMMs have been shown to accurately model sprays and soot with a relatively compact set of moments. However, significantly non-Gaussian processes such as bubble dynamics lead to numerical instabilities when extending their moment sets accordingly. We solve this problem by training a recurrent neural network (RNN) that adjusts the QBMM quadrature to evaluate unclosed moments with higher accuracy. The proposed method is tested on a simple model of bubbles oscillating in response to a temporally fluctuating pressure field. The approach decreases model-form error by a factor of 10 when compared to traditional QBMMs. It is both numerically stable and computationally efficient since it does not expand the baseline moment set. Additional quadrature points are also assessed, optimally placed and weighted according to an additional RNN. These points further decrease the error at low cost since the moment set is again unchanged.

@article{spratt21,
  author = {Spratt, J.-S. and Rodriguez, M. and Schmidmayer, K. and Bryngelson, S. H. and Yang, J. and Franck, C. and Colonius, T.},
  title = {Characterizing viscoelastic materials via ensemble-based data assimilation of bubble collapse observations},
  journal = {Journal of the Mechanics and Physics of Solids},
  file = {spratt-JMPS-21.pdf},
  year = {2021},
  volume = {152},
  pages = {104455},
  doi = {10.1016/j.jmps.2021.104455}
}

Bubble cavitation can induce such strain rates, and the resulting bubble dynamics are sensitive to the material properties. Thus, in principle, these properties can be inferred via measurements of the bubble dynamics. Estrada et al. (2018) demonstrated such bubble-dynamic high-strain-rate rheometry by using least-squares shooting to minimize the difference between simulated and experimental bubble radius histories. We generalize their technique to account for additional uncertainties in the model, initial conditions, and material properties needed to uniquely simulate the bubble dynamics. Ensemble-based data assimilation minimizes the computational expense associated with the bubble cavitation model. We test an ensemble Kalman filter (EnKF), an iterative ensemble Kalman smoother (IEnKS), and a hybrid ensemble-based 4D–Var method (En4D-Var) on synthetic data, assessing their estimations of the viscosity and shear modulus of a Kelvin-Voigt material. Results show that En4D–Var and IEnKS provide better moduli estimates than EnKF. Applying these methods to the experimental data of Estrada et al. (2018) yields similar material property estimates to those they obtained, but provides additional information about uncertainties. In particular, the En4D–Var yields lower viscosity estimates for some experiments, and the dynamic estimators reveal a potential mechanism that is unaccounted for in the model, whereby the viscosity is reduced in some cases due to material damage occurring at bubble collapse.

@article{bryngelson19_CPC,
  doi = {10.1016/j.cpc.2020.107396},
  volume = {266},
  year = {2021},
  pages = {107396},
  author = {Bryngelson, S. H. and Schmidmayer, K. and Coralic, V. and Maeda, K. and Meng, J. and Colonius, T.},
  title = {{MFC: A}n open-source high-order multi-component, multi-phase, and multi-scale compressible flow solver},
  journal = {Computer Physics Communications},
  file = {bryngelson-CPC-20.pdf}
}

MFC is an open-source tool for solving multi-component, multi-phase, and bubbly compressible flows. It is capable of efficiently solving a wide range of flows, including droplet atomization, shock-bubble interaction, and gas bubble cavitation. We present the 5- and 6-equation thermodynamically-consistent diffuse-interface models we use to handle such flows, which are coupled to high-order interface-capturing methods, HLL-type Riemann solvers, and TVD time-integration schemes that are capable of simulating unsteady flows with strong shocks. The numerical methods are implemented in a flexible, modular framework that is amenable to future development. The methods we employ are validated via comparisons to experimental results for shock-bubble, shock-droplet, and shock-water-cylinder interaction problems and verified to be free of spurious oscillations for material-interface advection and gas-liquid Riemann problems. For smooth solutions, such as the advection of an isentropic vortex, the methods are verified to be high-order accurate. Illustrative examples involving shock-bubble-vessel-wall and acoustic-bubble-net interactions are used to demonstrate the full capabilities of MFC.

@article{schmidmayer19,
  author = {Schmidmayer, K. and Bryngelson, S. H. and Colonius, T.},
  year = {2020},
  volume = {402},
  pages = {109080},
  journal = {Journal of Computational Physics},
  title = {An assessment of multicomponent flow models and interface capturing schemes for spherical bubble dynamics},
  file = {schmidmayer-JCP-20.pdf},
  doi = {10.1016/j.jcp.2019.109080}
}

Numerical simulation of bubble dynamics and cavitation is challenging; even the seemingly simple problem of a collapsing spherical bubble is difficult to compute accurately with a general, three-dimensional, compressible, multicomponent flow solver. Difficulties arise due to both the physical model and the numerical method chosen for its solution. We consider the 5-equation model of Allaire et al. and Massoni et al., the 5-equation model of Kapila et al., and the 6-equation model of Saurel et al. as candidate approaches for spherical bubble dynamics, and both MUSCL and WENO interface-capturing methods are implemented and compared. We demonstrate the inadequacy of the traditional 5-equation model for spherical bubble collapse problems and explain the corresponding advantages of the augmented model of Kapila et al. for representing this phenomenon. Quantitative comparisons between the augmented 5-equation and 6-equation models for three-dimensional bubble collapse problems demonstrate the versatility of the pressure-disequilibrium model. Lastly, the performance of the pressure-disequilibrium model for representing a three-dimensional spherical bubble collapse for different bubble interior/exterior pressure ratios is evaluated for different numerical methods. Pathologies associated with each factor and their origins are identified and discussed.

@article{bryngelson19_ML,
  doi = {10.1016/j.ijmultiphaseflow.2020.103262},
  year = {2020},
  volume = {127},
  pages = {103262},
  author = {Bryngelson, S. H. and Charalampopoulos, A. and Sapsis, T. P. and Colonius, T.},
  title = {A {G}aussian moment method and its augmentation via {LSTM} recurrent neural networks for the statistics of cavitating bubble populations},
  journal = {International Journal of Multiphase Flow},
  file = {bryngelson-IJMF-20.pdf}
}

Phase-averaged dilute bubbly flow models require high-order statistical moments of the bubble population. The method of classes, which directly evolve bins of bubbles in the probability space, are accurate but computationally expensive. Moment-based methods based upon a Gaussian closure present an opportunity to accelerate this approach, particularly when the bubble size distributions are broad (polydisperse). For linear bubble dynamics a Gaussian closure is exact, but for bubbles undergoing large and nonlinear oscillations, it results in a large error from misrepresented higher-order moments. Long short-term memory recurrent neural networks, trained on Monte Carlo truth data, are proposed to improve these model predictions. The networks are used to correct the low-order moment evolution equations and improve prediction of higher-order moments based upon the low-order ones. Results show that the networks can reduce model errors to less than 1 percent of their unaugmented values.

@article{bryngelson20_qbmm,
  author = {Bryngelson, S. H. and Colonius, T. and Fox, R. O.},
  title = {{QBMMlib}: A library of quadrature-based moment methods},
  year = {2020},
  journal = {SoftwareX},
  volume = {12},
  pages = {100615},
  file = {bryngelson-SoftX-20.pdf},
  doi = {10.1016/j.softx.2020.100615}
}

QBMMlib is an open source package of quadrature-based moment methods and their algorithms. Such methods are commonly used to solve fully-coupled disperse flow and combustion problems, though formulating and closing the corresponding governing equations can be complex. QBMMlib aims to make analyzing these techniques simple and more accessible. Its routines use symbolic manipulation to formulate the moment transport equations for a population balance equation and a prescribed dynamical system. However, the resulting moment transport equations are unclosed. QBMMlib trades the moments for a set of quadrature points and weights via an inversion algorithm, of which several are available. Quadratures then closes the moment transport equations. Embedded code snippets show how to use QBMMlib, with the algorithm initialization and solution spanning just 13 total lines of code. Examples are shown and analyzed for linear harmonic oscillator and bubble dynamics problems.

@article{trummler19,
  doi = {10.1017/jfm.2020.432},
  title = {Near-surface dynamics of a gas bubble collapsing above a crevice},
  year = {2020},
  volume = {899},
  pages = {A16},
  author = {Trummler, T. and Bryngelson, S. H. and Schmidmayer, K. and Schmidt, S. J. and Colonius, T. and Adams, N. A.},
  journal = {Journal of Fluid Mechanics},
  file = {trummler-JFM-20.pdf}
}

The impact of a collapsing gas bubble above rigid, notched walls is considered. Such surface crevices and imperfections often function as bubble nucleation sites, and thus have a direct relation to cavitation-induced erosion and damage structures. A generic configuration is investigated numerically using a second-order-accurate compressible multi-component flow solver in a two-dimensional axisymmetric coordinate system. Results show that the crevice geometry has a significant effect on the collapse dynamics, jet formation, subsequent wave dynamics, and interactions. The wall-pressure distribution associated with erosion potential is a direct consequence of development and intensity of these flow phenomena.

@article{bryngelson19_whales,
  doi = {10.1121/10.0000746},
  year = {2020},
  volume = {147},
  number = {2},
  pages = {1126--1135},
  author = {Bryngelson, S. H. and Colonius, T.},
  title = {Simulation of humpback whale bubble-net feeding models},
  journal = {Journal of the Acoustical Society of America},
  file = {bryngelson-JASA-20.pdf}
}

Humpback whales can generate intricate bubbly regions, called bubble nets, via blowholes. Humpback whales appear to exploit these bubble nets for feeding via loud vocalizations. A fully-coupled phase-averaging approach is used to model the flow, bubble dynamics, and corresponding acoustics. A previously hypothesized waveguiding mechanism is assessed for varying acoustic frequencies and net void fractions. Reflections within the bubbly region result in observable waveguiding for only a small range of flow parameters. A configuration of multiple whales surrounding and vocalizing towards an annular bubble net is also analyzed. For a range of flow parameters, the bubble net keeps its core region substantially quieter than the exterior. This approach appears more viable, though it relies upon the cooperation of multiple whales. A spiral bubble net configuration that circumvents this requirement is also investigated. The acoustic wave behaviors in the spiral interior vary qualitatively with the vocalization frequency and net void fraction. The competing effects of vocalization guiding and acoustic attenuation are quantified. Low void fraction cases allow low-frequency waves to partially escape the spiral region, with the remaining vocalizations still exciting the net interior. Higher void fraction nets appear preferable, guiding even low-frequency vocalizations while still maintaining a quiet net interior.

@article{bryngelson19_IJMF,
  author = {Bryngelson, S. H. and Schmidmayer, K. and Colonius, T.},
  year = {2019},
  volume = {115},
  pages = {137-143},
  journal = {International Journal of Multiphase Flow},
  title = {A quantitative comparison of phase-averaged models for bubbly, cavitating flows},
  file = {bryngelson-IJMF-19.pdf},
  doi = {10.1016/j.ijmultiphaseflow.2019.03.028}
}

We compare the computational performance of two modeling approaches for the flow of dilute cavitation bubbles in a liquid. The first approach is a deterministic model, for which bubbles are represented in a Lagrangian framework as advected features, each sampled from a distribution of equilibrium bubble sizes. The dynamic coupling to the liquid phase is modeled through local volume averaging. The second approach is stochastic; ensemble-phase averaging is used to derive mixture-averaged equations and field equations for the associated bubble properties are evolved in an Eulerian reference frame. For polydisperse mixtures, the probability density function of the equilibrium bubble radii is discretized and bubble properties are solved for each representative bin. In both cases, the equations are closed by solving Rayleigh–Plesset-like equations for the bubble dynamics as forced by the local or mixture-averaged pressure, respectively. An acoustically excited dilute bubble screen is used as a case study for comparisons. We show that observables of ensemble- and volume-averaged simulations match closely and that their convergence is first order under grid refinement. Guidelines are established for phase-averaged simulations by comparing the computational costs of methods. The primary costs are shown to be associated with stochastic closure; polydisperse ensemble-averaging requires many samples of the underlying PDF and volume-averaging requires repeated, randomized simulations to accurately represent a homogeneous bubble population. The relative sensitivities of these costs to spatial resolution and bubble void fraction are presented.

@article{bryngelson19_chaos,
  author = {Bryngelson, S. H. and Gu\'{e}niat, F. and Freund, J. B.},
  volume = {100},
  pages = {012203},
  year = {2019},
  journal = {Physical Review E},
  title = {Irregular dynamics of cellular blood flow in a model microvessel},
  file = {bryngelson-PRE-19.pdf},
  doi = {10.1103/PhysRevE.100.012203}
}

The flow of red blood cells within cylindrical vessels is complex and irregular, so long as the vessel diameter is somewhat larger than the nominal cell size. Long-time-series simulations, in which cells flow 10,000 vessel diameters, are used to characterize the chaotic kinematics, particularly to inform reduced-order models. The simulation model used includes full coupling between the elastic red blood cell membranes and surrounding viscous fluid, providing a faithful representation of the cell-scale dynamics. Results show that the flow has neither classifiable recurrent features nor a dominant frequency. jInstead, its kinematics are sensitive to the initial flow configuration in a way consistent with chaos and Lagrangian turbulence. Phase-space reconstructions show that a low-dimensional attractor does not exist, so the observed long-time dynamics are effectively stochastic. Based on this, a simple Markov chain model for the dynamics is introduced and shown to reproduce the statistics of the cell positions.

@article{bryngelson19_LAOE,
  author = {Bryngelson, S. H. and Freund, J. B.},
  volume = {77},
  pages = {171-176},
  year = {2019},
  journal = {European Journal of Mechanics B/Fluids},
  title = {Non-modal {F}loquet stability of a capsule in large amplitude oscillatory extension},
  file = {bryngelson-EJM-19.pdf},
  doi = {10.1016/j.euromechflu.2019.04.012}
}

We analyze the stability of a capsule in large-amplitude oscillatory extensional (LAOE) flow, as often used to study the rheology and dynamics of suspensions. Such a flow is typically established in a cross-slot configuration, with the particle (or particles) of interest observed in the stagnation region. However, controlling this configuration is challenging because the flow is unstable. We quantify such an instability for spherical elastic capsules suspended near the stagnation point using a non-modal global Floquet analysis, which is formulated to include full coupling of the capsule-viscous-flow dynamics. The flow is shown to be transiently, though not asymptotically, unstable. For each case considered, two predominant modes of transient amplification are identified: a predictable intra-period growth for translational capsule perturbations and period-to-period growth for certain capsule distortions. The amplitude of the intra-period growth depends linearly on the flow strength and oscillation period, which corresponds to a shift of the flow stagnation point, and the period-to-period growth saturates over several periods, commensurate with the asymptotic stability of the flow.

@article{bryngelson18_JFM,
  author = {Bryngelson, S. H. and Freund, J. B.},
  journal = {Joural of Fluid Mechanics},
  pages = {663-677},
  title = {Floquet stability analysis of capsules in viscous shear flow},
  volume = {852},
  year = {2018},
  file = {bryngelson-JFM-18.pdf},
  doi = {10.1017/jfm.2018.574}
}

Observations in experiments and simulations show that the kinematic behaviour of an elastic capsule, suspended and rotating in shear flow, depends upon the flow strength, the capsule membrane material properties and its at-rest shape. We develop a linear stability description of the periodically rotating base state of this coupled system, as represented by a boundary integral flow formulation with spherical harmonic basis functions describing the elastic capsule geometry. This yields Floquet multipliers that classify the stability of the capsule motion for elastic capillary numbers Ca ranging from 0.01 to 5. Viscous dissipation rapidly damps most perturbations. However, for all cases, a single component grows or decays slowly, depending upon Ca, over many periods of the rotation. The transitions in this stability behaviour correspond to the different classes of rotating motion observed in previous studies.

@report{report_1,
  author = {Bryngelson, S. H. and Pantano, C. and Bodony, D. and Freund, J. B.},
  title = {Adjoint-based sensitivity for flows with shocks},
  institution = {Technical Report, XPACC},
  year = {2018},
  file = {bryngelson-xpacc-18.pdf}
}

Developing a consistent solution method for discontinuous adjoint flow problems is challenging. 
The inviscid Burgers’ equation is considered as a step toward the formulation of such a method. Results are shown for linear and nonlinear numerical methods, including WENO shock capturing. Necessary and sufficient conditions for consistent and convergent solutions to the associated adjoint equation are discussed. The adjoint Euler equations are also investigated, for which no numerical schemes are provably convergent. A characteristic-based method for this problem is proposed. It transforms the adjoint equations into an uncoupled set of transport equations. These equations have the same form as the adjoint Burgers’ equation, and thus inherit their proven consistency properties.

@article{bryngelson18_PRF,
  author = {Bryngelson, S. H. and Freund, J. B.},
  journal = {Physical Review Fluids},
  number = {7},
  pages = {073101},
  title = {Global stability of flowing red blood cell trains},
  volume = {3},
  year = {2018},
  file = {bryngelson-PRF-18.pdf},
  doi = {10.1103/PhysRevFluids.3.073101}
}

A train of red blood cells flowing in a round tube will either advect steadily or break down into a complex and irregular flow, depending upon its degree of confinement. We analyze this apparent instability, including full coupling between the viscous fluid flow and the elastic cell membranes. A linear stability analysis is constructed via a complete set of orthogonal perturbations to a boundary integral formulation of the flow equations. Both transiently and asymptotically amplifying disturbances are identified. Those that amplify transiently have short-wavelength shape distortions that carry significant membrane strain energy. In contrast, asymptotic disturbances are primarily rigid-body-like tilts and translations. It is shown that an intermediate cell-cell spacing of about half a tube diameter suppresses long-time train instability, particularly when the vessel diameter is relatively small. Altering the viscosity ratio between the cytosol fluid within the cell and the suspending fluid is found to be asymptotically destabilizing for both higher and lower viscosity ratios. Altering the cytosol volume away from that of a nominally healthy discocyte alters the stability with complex dependence on train density and vessel diameter. Several of the observations are consistent with a switch from predominantly cell-cell interactions for dense trains and predominantly cell-wall interactions for less dense trains. Direct numerical simulations are used to verify the linear stability analysis and track the perturbation growth into a self-sustaining disordered regime.

@article{bryngelson16_PRF,
  author = {Bryngelson, S. H. and Freund, J. B.},
  journal = {Physical Review Fluids},
  number = {3},
  pages = {033201},
  title = {Capsule-train stability},
  volume = {1},
  year = {2016},
  file = {bryngelson-PRF-16.pdf},
  doi = {10.1103/PhysRevFluids.1.033201}
}

Elastic capsules flowing in small enough tubes, such as red blood cells in capillaries, are well known to line up into regular single-file trains. The stability of such trains in somewhat wider channels, where this organization is not observed, is studied in a two-dimensional model system that includes full coupling between the viscous flow and suspended capsules. A diverse set of linearly amplifying disturbances, both long-time asymptotic (modal) and transient (nonmodal) perturbations, is identified and analyzed. These have a range of amplification rates and their corresponding forms are wavelike, typically dominated by one of five principal perturbation classes: longitudinal and transverse translations, tilts, and symmetric and asymmetric shape distortions. Finite-amplitude transiently amplifying perturbations are shown to provide a mechanism that can bypass slower asymptotic modal linear growth and precipitate the onset of nonlinear effects. Direct numerical simulations are used to verify the linear analysis and track the subsequent transition of the regular capsule trains into an apparently chaotic flow.

@article{bryngelson16_RA,
  author = {Bryngelson, S. H. and Freund, J. B.},
  journal = {Rheologica Acta},
  number = {6},
  pages = {451-464},
  title = {Buckling and its effect on the confined flow of a model capsule suspension},
  volume = {55},
  year = {2016},
  file = {bryngelson-RA-16.pdf},
  doi = {10.1007/s00397-015-0900-9}
}

The rheology of confined flowing suspensions, such as blood, depends upon the dynamics of the components, which can be particularly rich when they are elastic capsules. Using boundary integral methods, we simulate a two-dimensional model channel through which flows a dense suspension of fluid-filled capsules. A parameter of principal interest is the equilibrium membrane perimeter, parameterized by xi_o, which ranges from round capsules with xi_o=1.0 to xi_o=3.0 capsules with a dog-bone-like equilibrium shape. It is shown that the minimum effective viscosity occurs for xi_o approx 1.6, which forms a biconcave equilibrium shape, similar to a red blood cell. The rheological behavior changes significantly over this range; transitions are linked to specific changes in the capsule dynamics. Most noteworthy is an abrupt change in behavior for xi_o = 2.0, which correlates with the onset of capsule buckling. The buckled capsules have a more varied orientation and make significant rotational (rotlet) contributions to the capsule–capsule interactions.