# Seminar: Numerical Methods for Partial Differential Equations

### Organizer:

Rodolfo Ruben Rosales, Keaton Burns

### Mailing List:

If you are interested in receiving announcements by email, please write to daisymae@math.mit.edu

### Room

2-449

## Schedule Spring 2023

### Wenesday, May 3, 2023

4:30 PM - 5:30 PM

**Room:** 2-449

Abinand Gopal (Yale)

An accelerated, high-order solver for acoustic scattering in the plane

Time-harmonic acoustic scattering in variable media can often be modeled by variable-coefficient Helmholtz equations. Such equations pose several numerical challenges. For example, they can be intrinsically ill-conditioned, necessitate the imposition of radiation conditions, and produce pollution errors when discretized with standard finite difference or finite element methods.

To avoid these issues, it is often convenient to reformulate the Helmholtz equation as an integral equation known as the Lippmann-Schwinger equation. The tradeoffs are that an integral operator with a singular kernel must be discretized and that the resulting linear system that must be inverted is dense. In this talk, I will focus on the latter issue and present a new direct solver for this purpose.

For a problem with fixed wavenumber and N degrees of freedom, the number of operations required for the solver scales as O(N^(3/2)) with very favorable constants. Moreover, the solver is highly stable in practice. The solver is particularly effective when used to compute a preconditioner for an iterative solver. In this regime, I will show that problems up to 500 by 500 wavelengths with large amounts of back scattering can be solved in mere hours on a workstation. Time permitting, I will also discuss some work in progress that extends this solver to Maxwell's equations.

This is joint work with Gunnar Martinsson (UT Austin) and the extension to Maxwell's equation is joint work with Hanwen Zhang (Yale).

### Wednesday, May 10, 2023

4:30 PM - 5:30 PM

**Room:** 2-449

Tess Smidt (MIT)

Euclidean Symmetry Equivariant Machine Learning for Atomic Systems –- Overview, Applications, and Open Questions

Atomic systems (molecules, crystals, proteins, etc.) are naturally represented by a set of coordinates in 3D space labeled by atom type. This is a challenging representation to use for machine learning because the coordinates are sensitive to 3D rotations, translations, and inversions (the symmetries of 3D Euclidean space). In this talk I’ll give an overview of Euclidean invariance and equivariance in machine learning for atomic systems. Then, I’ll share some recent applications of these methods on a variety of atomistic modeling tasks (ab initio molecular dynamics, prediction of crystal properties, and scaling of electron density predictions). Finally, I’ll explore open questions in expressivity, data-efficiency, and trainability of methods leveraging invariance and equivariance.

## Schedule Fall 2022

### Wednesday, September 7, 2022

4:30pm - 5:30pm

**Room:** 2-449

**Zoom:** https://mit.zoom.us/j/99611578020

Sidney Holden (University of Sydney)

A continuum limit for metric graphs

A metric graph is what one might expect intuitivelya spiderweb, electrical power grid, water-supply system or the British Rail. Physical continuity allows posing one-dimensional differential equations along each edge. Though rather than simple boundary conditions, solutions must satisfy conditions at the vertices which depend on the dynamics along incident edges

We study the eigenvalues and eigenfunctions of the Laplace operator on general metric graphs, which involves solving a system of nonlinear algebraic vertex conditions. In particular, we consider the continuum limit in which graphs become dense within some prescribed embedding space, e.g. a spiderweb filling in the unit disc. As the number of vertices increases, numerical solutions require a novel matrix-determinant root-finding algorithm. Alternatively, in the continuum limit, the discrete system on the vertices reduces to the eigenvalue equation of a contiguous partial differential operator. This resembles in many ways a Laplace-Beltrami operator on a Riemannian manifold but with several notable differences. Rather than a metric tensor, we derive a distinct symmetric tensor that scales differently with distance. Rather than the determinant-based volume form, we find an analogous matrix-trace-based distance form. Our findings open the possibility for a new kind of manifold geometry similar to geodesic structure, but made from underlying “graph material”. We discuss a series of examples of high-density networks, comparing PDE solutions to numerical solutions of the vertex systems. This is joint work with Geoffrey Vasil (U. Edinburgh).

### Wednesday, September 21, 2022

4:30pm - 5:30pm

**Room:** 2-449

**Zoom:** https://mit.zoom.us/j/99160816852

Miles Cranmer (Princeton University)

Interpretable Machine Learning for Physics

Would Kepler have discovered his laws if machine learning had been around in 1609? Or would he have been satisfied with the accuracy of some black box regression model, leaving Newton without the inspiration to find the law of gravitation? In this talk I will present a review of some industry-oriented machine learning algorithms, and discuss a major issue facing their use in the natural sciences: a lack of interpretability. I will then outline several approaches I have created with collaborators to help address these problems, based largely on a mix of structured deep learning and symbolic methods. This will include an introduction to the PySR/SymbolicRegression.jl software (https://astroautomata.com/PySR), a Python/Julia package for high-performance symbolic regression. I will conclude by demonstrating applications of such techniques and how we may gain new insights from such results.

### Wednesday, September 28, 2022

4:30pm - 5:30pm

**Room:** 2-449

**Zoom:** https://mit.zoom.us/j/99160816852

Anuj Kumar (UCSC)

Optimal scalar transport using steady branching pipe flows and unsteady branching blob flows

We consider the problem of "wall-to-wall optimal transport," in which we attempt to maximize the transport of a passive temperature field between hot and cold plates. Specifically, we optimize the choice of the divergence-free velocity field in the advection-diffusion equation amongst all velocities satisfying an enstrophy constraint (which can be understood as a constraint on the power required to generate the flow). Previous work established an a priori upper bound on the transport, scaling as the 1/3-power of the flow's enstrophy. Recently, Doering & Tobasco'19 constructed self-similar two-dimensional steady branching flows saturating this bound up to a logarithmic correction to scaling. This logarithmic correction appears to arise due to a topological obstruction inherent to two-dimensional steady branching flows. We present a construction of three-dimensional ``branching pipe flows" that eliminates the possibility of this logarithmic correction and therefore identifies the optimal scaling as a clean 1/3-power law. Also, using an unsteady "branching blob flow" construction, it appears that the 1/3 scaling is, in fact, optimal in two dimensions as well. We discuss the underlying physical mechanism that makes the branching flows "efficient" in transporting heat and present a design of a mechanical apparatus that can transfer heat close to the best possible scenario.

### Wednesday, October 5, 2022

4:30pm - 5:30pm

**Zoom:**https://mit.zoom.us/j/99160816852

David I. Ketcheson (King Abdullah University)

Invariant-preserving methods for dispersive wave equations

Many dispersive wave equations have a Hamiltonian structure, with important conserved quantities of both linear and nonlinear type. The accuracy of numerical discretizations over long times depends critically on the preservation of these invariants. I will present a class of invariant-preserving discretizations based on summation-by-parts in space and relaxation in time; these fully-discrete schemes conserve both linear and nonlinear invariants and can be either implicit or explicit. These schemes have been developed for a wide range of such equations, including Benjamin-Bona-Mahony (BBM), Fornberg-Whitham, Camassa-Holm, Degasperis-Procesi, Holm-Hone, and the BBM-BBM system. I will show examples demonstrating that the error for such schemes grows only linearly in time, whereas for general schemes the error grows quadratically. I will also show examples of non-dispersive hyperbolic systems where such invariant-preserving schemes lead to a similar (drastic) improvement in long-time accuracy. Finally, I will present some preliminary results of conserving multiple nonlinear invariants in multi-soliton solutions.

### Wednesday, November 16, 2022

11:00am -12:00 Noon

**Room:** 2-449

Matthew Colbrook (Cambridge University)

The mpEDMD algorithm: Structure-preserving and rigorous Koopmanism

Koopman operators globally linearize nonlinear dynamical systems, and their spectral information can be a powerful tool for the analysis and decomposition of nonlinear dynamical systems. However, Koopman operators are infinite- dimensional, and computing their spectral information is a considerable challenge. We introduce measure-preserving extended dynamic mode decomposition (mpEDMD), the first Galerkin method whose eigendecomposition converges to the spectral quantities of Koopman operators for general measure-preserving dynamical systems. mpEDMD is a data-driven and structure-preserving algorithm based on an orthogonal Procrustes problem that enforces measure-preserving truncations of Koopman operators using a general dictionary of observables. It is flexible and easy to use with any pre-existing DMD-type method and with different data types. We prove the convergence of mpEDMD for projection-valued and scalar-valued spectral measures, spectra, and Koopman mode decompositions. For the case of delay embedding (Krylov subspaces), our results include convergence rates of the approximation of spectral measures as the size of the dictionary increases. We demonstrate mpEDMD on a rang of challenging examples, its increased robustness to noise compared with other DMD-type methods, and its ability to capture the energy conservation and cascade of a turbulent boundary layer flow with Reynolds number > 60000 and state-space dimension >100000. Finally, if time permits, we discuss how this algorithm forms part of a broader program on the foundations of infinite-dimensional spectral computations.

### Wednesday, December 7, 2022

4:15pm -5:15pm

**Room:** 2-449

Raphael Pestourie (MIT)

Combining Data and Models to Accelerate Simulations and Enable Inverse Design

Inverse design is the direct optimization of a target property, it has the potential to automatize design for real-world engineering problems. However, inverse design is limited by the simulation capabilities of physical phenomena. In an iterative optimization process, a computer model of a real-world phenomenon is queried hundreds to thousands of times. On the one hand, in many applications, there exist physical models where the simulations are accurate enough to result in a meaningful design, but these simulations may be too resource-intensive to run the optimization process. This is the case for example in metasurface designoptical devices that present both subwavelength aperiodic patterns and a thousands-of-wavelengths-long diameter. On the other hand, data-driven modls are most often very fast to evaluate and they do not require a complete knowledge of the process being optimized. Unfortunately, as the number of real-world design parameters increases, these models may require unreasonable amounts of data and resources to be trained accurately. In this talk, I will present methodologies that combine data and models to accelerate simulations and enable inverse design. The center-piece of these approaches is the creation of data- and resource-efficient global surrogate models that are repeatedly called in the optimization loop. I will show theory and experimental results of surrogate-based large-scale metasurface designs. I will also show recent work in active learning and scientific machine learning (SciML) for various engineering problems, where information from both physical models and data are optimally leveraged to achieve efficient inverse design. Upon inspection, SciML may also give insights and intuition about the engineering process being simulated and optimized.

## Schedule Spring 2022

### Monday, May 16, 2022

3:30pm - 4:30pm

**Room:** 2-449

**Zoom:** https://mit.zoom.us/j/99059307672

Eric Hester (UCLA)

Modelling fluid-solid interactions

Fluid-solid interactions are a continuing challenge for numerical modelling. The main difficulty stems from imposing boundary conditions at the fluid-solid interface. This talk discusses my research on better mathematical and computational techniques for these problems.

I will start with the elegant differential geometry of the signed distance function, and show how it simplifies analysis of general boundary layers. Next, I will demonstrate how this boundary layer analysis can derive improved methods for simulating fluid-solid interactions. Finally, I will examine these techniques in two real-world problems: improving predictions of iceberg melting, and understanding the extreme boat drag of the dead water effect.

## Schedule Fall 2021

**Meeting Time:** Wednesdays, 4:30-5:30pm, Building 2, Room 449

### Wednesday, November 17th, 2021

4:30pm - 5:30pm

**Room:** 2-449

Peter Baddoo (MIT)

Integrating physical laws into data-driven modal decompositions

Incorporating partial knowledge of physical laws into data-driven architectures can improve the accuracy, generalisability, and robustness of the resulting models. In this work, we demonstrate how physical laws – such as symmetries, invariances, and conservation laws – can be integrated into the dynamic mode decomposition (DMD), which is a widely-used data analysis technique that extracts modal structures from high-dimensional measurements. Although DMD has been applied successfully to a broad range of domains, the algorithm frequently produces models that are sensitive to noise, fail to generalize outside the training regime, and violate basic physical laws. We rephrase the DMD optimization by restricting the family of admissible models to a matrix manifold that respects the physical structure of the system at hand. Termed 'physics-informed dynamic mode decomposition’ (piDMD), our formulation may be interpreted as a Procrustes problem, which allows us to leverage the substantial existing literature thereof. We focus on five of the most fundamental physical properties – conservative, self-adjoint, local, causal, and shift-invariant – and derive several closed-form solutions for the corresponding piDMD optimization problems. We demonstrate piDMD on a range of canonical problems in the physical sciences, including energy-preserving fluid flow, travelling-wave systems, the Schrödinger equation, solute advection-diffusion and systems with causal dynamics. We conclude with the challenging example of identifying the spectrum of the linearized Navier--Stokes equations from transitional channel flow data. In each case, piDMD significantly outperforms standard DMD in metrics such as spectral identification, state reconstruction and prediction, and estimation of optimal forcings and responses. If time permits, we will also discuss new solvers for elliptic PDEs based on rational approximations.

### Wednesday, November 10th, 2021

4:30pm - 5:30pm

**Room:** 2-449

**Zoom:** https://mit.zoom.us/j/98343208570

Andrew J. Horning (MIT)

Computing spectra of infinite-dimensional operators

Computing the spectrum of a differential or integral operator is often done in two steps: (1) discretize the operator to obtain a matrix eigenvalue problem and (2) compute eigenvalues of the matrix with numerical linear algebra. This “discretize-then-solve” paradigm is flexible and powerful, but tension between spectral properties of the operator and the matrix discretizations can lead to numerical artifacts that pollute computed spectra and degrade accuracy. Moreover, the “discretize-then-solve” paradigm may fail to capture infinite-dimensional phenomena like continuous spectra. In this talk, we present a robust computational framework that extracts both discrete and continuous spectral properties of infinite-dimensional operators by sampling the range of the resolvent at strategic points in the complex plane. Approximations to eigenvalues and eigenfunctions, spectral measures, and generalized eigenfunctions are constructed by (1) solving linear equations with complex shifts and (2) taking inner products in a Hilbert space. The resulting algorithms respect key structure from the operator, regardless of the underlying approximation schemes used to execute these two tasks on a computer. We illustrate the approach through a series of examples, highlighting numerical implementations for a broad class of differential and integral operators available in SpecSolve, a MATLAB package leveraging state-of-the-art spectral methods for fast and accurate spectral computations in infinite dimensions.

## Schedule Spring 2021

**Meeting Time:** Wednesdays, 4:30-5:30pm, Virtual.

### Wednesday, May 12th, 2021

4:30pm - 5:30pm

Nicholas Boffi (Harvard University)

Projection methods for quasi-static hypo-elastoplasticity

When a material deforms elastically, it returns to its original configuration after removal of an applied external force. When a material deforms plastically, it deforms irreversibly, remaining in the deformed configuration after removal of an applied load. Hypo-elastoplasticity is a flexible framework for describing hard materials exhibiting a combination of elastic and plastic deformation with a dominant plastic component. In this talk, we explore an analogy between the Navier-Stokes equations in the incompressible limit and the dynamics of a hypo-elastoplastic material in the quasi-static limit. By exploiting correspondences between 1) the hypo-elastoplastic material stress tensor and the fluid velocity field, and 2) the hypo-elastoplastic material velocity field and the fluid pressure field, we develop a projection algorithm for quasi-static hypo-elastoplasticity analogous to Chorin's projection method for incompressible fluid dynamics. After describing the basic formulation, we discuss an extension of the method, which enables the implementation of complex material deformation through an explicit coordinate transformation from a reference domain to the physical domain. We illustrate the efficacy of the methods through numerical experiments in modeling three-dimensional shear band formation in an amorphous metal.

### Wednesday, May 5th, 2021

* Special Time:* 1:00pm - 2:00pm

Felipe Vico (Universitat Politècnica de València)

High-Order Multiscale Mesh Generation And Applications To Boundary Value Problems In Electromagnetics

The scattering of electromagnetic waves in the presence of conductors and dielectrics is a fundamental problem in classical physics that has a wide range of applications in engineering. Potential theory leads to particularly powerful integral equation techniques for such problems. In this talk, we will review the main aspects of these methods: integral formulations, numerical discretization, mesh generation and fast algorithms. We will put particular emphasis on a recently developed fast mesh generation method for creating high order approximations of realistic geometries. We will show some applications and numerical examples.

### Wednesday, April 21th, 2021

4:30pm - 5:30pm

Ben Adcock (Simon Fraser University)

Approximation of high-dimensional functions via sparse polynomials and deep neural networks, with application to parametric PDEs

Driven by its various applications in scientific computing in particular, the solution of parametric and stochastic DEs arising in uncertainty quantification the approximation of smooth, high-dimensional functions via sparse polynomial expansions has received significant attention in the last decade. In the first part of this talk I will give a brief survey of recent progress in this area. In particular, I will show how the proper use of compressed sensing tools leads to algorithms for high-dimensional approximation which, unlike other approaches, provably possess near-optimal error bounds and moderate sample complexities. In particular, these techniques mitigate the curse of dimensionality to a substantial degree. The second part of the talk is devoted to most recent approaches based on deep neural networks and deep learning. Such tools are currently garnering substantial attention in the scientific computing community. Nonetheless, I will show evidence of a gap between current theory and practice. I will then discuss recent theoretical contributions showing that deep neural networks matching the performance of best-in-class schemes can be computed. This highlights the potential of deep neural networks, and sheds light on achieving robust, reliable and overall improved practical performance.

## Schedule Spring 2020

**Meeting Time:** Wednesdays, 4:30-5:30pm, Room 2-136, Refreshments will be served in Room 2-290 at 3:30pm.

### Wednesday, April 29th, 2020

9:00am - 10:00am

Albert Chern (Technical University of Berlin)

An Exact Discretization of Reflectionless Boundaries for Wave Equations

This talk concerns a classical problem in computational wave propagations: How does one truncate an infinite domain to a finite size without introducing reflection waves from the artificial boundaries? The state-of-the-art approach is attaching to those boundaries a perfectly matched layer (PML). In the continuous theory, PMLs are subject to an analytically continued wave equation that damps all incident waves without creating any interfacial reflection. However, it is believed that "numerical reflections" are unavoidable after discretization. In this talk, I will demonstrate a truly reflectionless discrete PML. The key is to uncover the geometric mechanism hidden in the differential calculus formalism; in discretizing the theory, approximations are the best one can hope for the latter, while the former often admits exact discretization.

Bio: Albert Chern; currently a postdoc in mathematics at Technical University of Berlin, and will start as an assistant professor in computer science at UC San Diego in Fall 2020. Research area: computational math, discrete differential geometry, fluid dynamics, computer graphics.

### Wednesday, January 29th, 2020

4:30pm - 5:30pm in "New Location" Room 2-132

Carlos Perez Arancibia (Institute for Mathematical and Computational Engineering, PUC Chile )

Density Interpolation Methods

In this talk, I will present present recent developments on a class of effective and simple-to-implement methods for the numerical evaluation of boundary integral operators and layer potentials associated to classical PDEs, that address longstanding efficiency, accuracy, and practical implementation issues that have hindered the applicability of boundary integral equation methods in science and engineering. The proposed techniques rely on the use of Green's third identity and local Taylor-like interpolations of density functions in terms of homogeneous solutions of the underlaying PDE. These techniques effectively regularize the singularities present in boundary integral operator and layer potentials, recasting them in terms of integrands that are bounded or even more regular depending on the order of the density interpolation. The resulting boundary integrals can then be easily, accurately, and inexpensively evaluated by means of standard quadrature rules. A variety of numerical examples demonstrate the effectiveness of these techniques in the context of Nystrom and boundary element methods for the Laplace, Helmholtz, and Maxwell equations.

This is joint work with Catalin Turc (NJIT), Luiz Faria (INRIA, ENSTA ParisTech), and Constantine Sideris (USC).

## Schedule Fall 2019

**Meeting Time:** 4:30pm - 5:30pm in Room 2-136

### Wednesday, December 11th, 2019

4:30pm - 5:30pm in Room 2-136

Daniel Fortunaton (Harvard University)

The ultraspherical spectral element method

We introduce a novel spectral element method based on the ultraspherical spectral method and the hierarchical Poincare-Steklov scheme for solving general partial differential equations on polygonal unstructured meshes. Whereas traditional finite element methods have a computational complexity that scales as O(p^6) with polynomial degree p, our method scales as O(p^4), allowing for hp-adaptivity to be based on physical considerations instead of computational ones. Properties of the ultraspherical spectral method lead to almost banded well-conditioned linear systems, allowing for the element method to be competitive in the ultra-high-polynomial regime. The hierarchical Poincare-Steklov scheme enables precomputed solution operators to be reused, allowing for fast elliptic solves in implicit and semi-implicit time-steppers. We develop an open-source software system to demonstrate the flexibility and robustness of the method. Joint work with Alex Townsend and Nick Hale.

### Tuesday, November 5th, 2019

SPECIAL TIME and PLACE

2:30pm - 3:30pm in Room 2-131

Dionisios Margetis (University of Maryland, College Park)

**NOTE: This seminar is joint with the Physical Mathematics Seminar**

Plasmonics on two-dimensional materials: Flavors of dispersion and homogenization

The advent of two-dimensional materials such as graphene and black phosphorus with a wide range of optical and electronic properties offers the promise of excellent efficiency in light-matter interaction at the nanoscale. These systems may allow for the propagation of fine-scale electromagnetic waves, called surface plasmon-polaritons (SPPs), which can defy the usual diffraction limit.

In this talk, I will discuss recent theoretical progress in understanding how the material geometry, e.g., the presence of edges and curvature as well as the formation of periodic layered structures in the presence of material anisotropy, may affect the SPP propagation. To this end, I will formulate and analyze physically inspired boundary value problems for the time-harmonic Maxwell equation.

### Wednesday, September 25th, 2019

4:30pm - 5:30pm in Room 2-136

Rasmus Ellebaek Christiansen (Technical University of Denmark)

Topology Optimization: Not Just Throwing Parameters at the Wall

The engineering of physical systems, such as the load-carrying members of an airplane or the optical components of a microscope, is a cornerstone in the improvement and development of new devices and the advancement of science. Computational design tools now allow us to rapidly explore a huge space of non-intuitive designs, but such large-scale optimization approaches require careful formulation of the problem and the design space in order to yield tractable and usable results.

This talk explains the technique of topology optimization (TO) for device design, illustrating it for a variety of engineering problems in solid mechanics and electromagnetism. In particular, it will focus on density-based TO and associated methods for imposing manufacturing constraints (e.g. minimum feature sizes) and rapid optimization (e.g. adjoint-based sensitivity analysis). A typical engineering problem cannot simply be "thrown" at TO, however - one must often reformulate the engineering objective to be well suited for TO, addressing concerns such as differentiability and density interpolation. We illustrate this process using recent results of TO for design problems such as plasmonic field localization, Raman scattering, and photonic cavity design.

## Schedule Spring 2019

**Meeting Time: Wednesdays, 4:30-5:30pm, Room 2-136, Refreshments will be served in Room 2-290 at 3:30pm.**

### Wednesday, February 13th, 2019

4:30pm - 5:30pm in Room 2-136

Thomas Anderson (California Institute of Technology)

Fast hybrid frequency/time techniques for efficient and parallelizable high-order time-domain wave scattering

We propose and discuss a frequency/time hybrid integral-equation (though other frequency-domain methods are readily usable) method for the time-dependent wave equation in two and three-dimensional spatial domains. Relying on Fourier transformation in time, the method utilizes a fixed (time-independent) number of frequency-domain integral-equation solutions to evaluate time domain solutions for arbitrarily long times. The approach relies on two main elements, namely, 1) A smooth time-windowing methodology that enables accurate band-limited representations for arbitrary long time signals, and 2) A novel Fourier transform approach which, in a time-parallel manner and without causing spurious periodicity effects, delivers numerically dispersionless spectrally accurate solutions. The algorithm can tackle complex physical structures, it enables parallelization in time in a straightforward manner, and it allows for time leapingthat is, solution sampling at any given time T at O(1)-bounded sampling cost. The proposed frequency/time hybridization strategy, which generalizes to any linear partial differential equation in the time domain for which frequency-domain solutions can be obtained (including e.g. the time-domain Maxwell equations or time domain problems posed with dispersive media) provides significant advantages over other available alternatives such as volumetric discretization and convolution-quadrature approaches.

### Tuesday, April 30, 2019

Joint with Physical Math Seminar

2:30pm - 3:30pm in Room 2-139

G. Bard Ermentrout (University of Pittsburgh)

If space turned out to be time: Resonances in the visual cortex

When subjects are exposed to full field flicker in certain frequencies, they perceive a variety of complex geometric patterns that are often called flicker hallucinations. On the other had, when looking at high contrast geometric patterns like op art, shimmering and flickering is observed. In some people, flicker or such op art can induce seizures. In this talk, I describe a simple network model of excitatory and inhibitory neurons that comprise the visual area of the brain. I show that these phenomena are reproduced and then give an explanation based on symmetry breaking bifurcations and Floquet theory. Symmetric bifurcation theory also shows why one expects a different class of patterns at high frequencies than those at low frequencies. Next, I will describe the flip side of this coin and discuss a theory of uncomfortable images. Many people exhibit visual discomfort when looking at high contrast geometric patters such as seen in op art. I'll discuss some recent results where we show that such patterns can induce global oscillations in a network similar to the one used in the flicker study.

## Schedule Spring 2018

**Meeting Time:** Wednesdays, 4:30-5:30pm, Room 2-136, Refreshments will be served in Room 2-290 at 3:30pm.

### Wednesday, April 4th, 2018

4:30pm - 5:30pm in Room 2-136

Yuval Dagan (Massachusetts Institute of Technology)

Similarity of Droplet Clouds in Vortex Flows

The dynamics of droplets and particles in rotating flows is of fundamental importance for various engineering, environmental and medical applications. The problem of calculating particles and droplet trajectories advected by the flow is challenging from the point of view of mathematical analysis. While current efforts in this field rely intensively on detailed computational analysis, such as direct numerical simulations, some of the basic underlying physics still remain hidden. It is often useful to simplify this mathematical problem, while retaining its non-trivial features.

In this talk I will present a new simplified mathematical formulation for the dynamics of particles and evaporating droplets in the vicinity of vortex flows. The complex interaction between discrete droplets and the surrounding flow is resolved using a sectional Eulerian representation of the particles. Then, new similarity solutions for droplet dynamics in reacting and non-reacting vortical flows will be presented.

Finally, I will present a new concept for the stability analysis of reacting and non-reacting shear-flows in the vicinity of polydisperse droplets and particles. Using this approach, I will show the influence of droplet evaporation on the onset of shear flow instabilities.

### Wednesday, April 11th, 2018

4:30pm - 5:30pm in Room 2-136

Catalin Turc (New Jersey Institute of Technology)

Domain decomposition for quasi-periodic scattering by layered media via robust boundary-integral equations at all frequencies

We develop a non-overlapping domain decomposition method (DDM) for the solution of quasi-periodic scalar transmission problems in layered media. Our approach relies on robust boundary-integral equation formulations of Robin-to-Robin (RtR) maps throughout the frequency spectrum, including at Rayleigh, or cutoff, frequencies. We overcome the obstacle of non-convergent quasi-periodic Green functions at these frequencies by incorporating newly introduced shifted quasi-periodic Green functions. Using the latter in the definition of our quasi-periodic boundary-integral operators leads to rigorously stable computations of RtR operators. We develop Nystrom discretizations of the RtR maps that rely on trigonometric interpolation, singularity resolution, and fast convergent windowed quasi-periodic Green functions. We solve the tridiagonal DDM system via recursive Schur complements and we establish rigorously that this procedure is always completed successfully. We present a variety of numerical results concerning Wood frequencies in two and three dimensions as well as large numbers of layers. Joint work with Stephen Shipman, Carlos Perez-Arancibia, and Stephanos Venakides.

### Wednesday, April 25th, 2018

4:30pm - 5:30pm in Room 2-136

Josue Melguizo-Gavilanes (Centre National de la Recherche Scientifique (CNRS))

Thermal Ignition

The study of thermal ignition is key to predict, prevent, mitigate and assess the risk of accidental fires and explosions. Careful experiments of canonical/simple geometries have been performed over the years to gain fundamental understanding of the phenomenon. Parametric studies analyzing the effect of size of the hot surface, mixture concentration, geometry and ignition time scales on ignition thresholds have been carried out. However, detailed numerical simulations of the experiments to study the dynamics of ignition and the most important physical and chemical processes that control ignition had not been performed. Can we reproduce the experimental ignition thresholds and reaction regimes? When/where does ignition take place? What are the key processes taking place at the ignition location? Numerical simulations of the geometries studied experimentally (i.e. moving heated spheres, rapidly and slowly heated surfaces) were performed to answer these questions.

## Schedule Fall 2017

**Meeting Time:** Wednesdays, 4:30-5:30pm, Room 2-136, Refreshments will be served in Room 2-290 at 3:30pm.

### Wednesday, November 15th, 2017

4:30pm - 5:30pm in Room 2-136

Ryan N. Goh (Boston University)

Vortices in Rapidly Rotating Boussinesq Convection

We study long time asymptotics in the Boussinesq approximation for rapidly rotating stably-stratified fluids in a three dimensional infinite layer with either stress-free or periodic boundary conditions. For initial data satisfying certain smallness criteria, we use self-similar variables to show that the baroclinic vorticity converges to an Oseen Vortex with algebraic rate while the baroclinic component decays to zero faster than any algebraic rate. In the case of periodic boundary conditions, we also find that the barotropic vertical velocity and thermal fluctuations converge to decaying Gaussians whose amplitudes oscillate with opposite phase of each other. We also use dispersive estimates and Lyapunov functional techniques to determine asymptotics in the rapid rotation limit for a broader class of initial data where we only require smallness in the quasi-geostrophic part of the baroclinic dynamics.

### Thursday, November 16th, 2017

4:30pm - 5:30pm in Room 2-105

Dmitry Kabanov (King Abdullah University of Science & Technology)

Numerical computation of linear Instability of detonations

We propose to study linear stability of steady-state gaseous detonations by a method combining features of normal-mode analysis and numerical simulations together. The linearized governing equations coupled with the shock-evolution equation are solved in the shock-attached frame using a fifth-order algorithm in time and space. The computed time series of the perturbation of detonation velocity are processed using the dynamic mode decomposition to extract growth rates and frequencies (stability spectra). The method is applied to the reactive Euler equations with one-step chemistry, where we compare our results with those known from normal-mode analysis and find agreement to two significant digits for stability spectra. Besides, we obtain neutral stability boundaries for several values of specific-heats ratio demonstrating that as the latter is close to the Newtonian limit, stability boundary is determined by a composition of the neutral stability curves for the fundamental mode and the first overtone, while the boundary is determined by the fundamental mode only for larger values of specific-heats ratio. Also, we apply our method to a model with two-step chemistry where the second reaction is endothermic. For this model, two types of steady-state solutions are possible with a region of identical solutions between sonic point and the shock and a region between the sonic point and the end of the reaction zone where two solutions are different. We demonstrate that the stability spectra of two types of steady-state solutions are essentially the same as stability is determined only by a region between the sonic point and the shock.

### Wednesday, November 29th, 2017

4:30pm - 5:30pm in Room 2-136

Geoffrey Vasil (The University of Sydney)

Flexible tensor calculus on domains with coordinate singularities

On a two-sphere (for example), not all smooth functions of latitude and longitude represent non-singular physical quantities; the concept of longitude breaks down at the north and south poles. Truly well-behaved functions must contain special relationships between otherwise independent variables near coordinate singularities. In applications, capturing these dependencies can become numerically very challenging. Fortunately, the famous Spherical Harmonic functions perfectly satisfy the special behaviour (for scalar functions) near the poles and can be used as a complete spectral basis. But what happens for vectors, and higher-order tensors? What happens on the unit disk? What about the coordinate singularity at the centre of the 3D ball? Or, what about the corners of triangles and tetrahedra?

This talk will discuss how to represents functions, vectors, and arbitrary tensors on domains with coordinate singularities in an accurate and numerically efficient way. The solution is based on the remarkable properties of Jacobi polynomials; which include Legendre and Chebyshev polynomials as special cases, and Hermite and Laguerre as limiting cases. Using Jacobi polynomials, it is possible to construct basis sets that account for a wide range of singular behaviour. These bases inherit simple and sparse differentiation and multiplication formula from the algebra of Jacobi polynomials, which allows for the efficient solution of PDEs.

### Wednesday, December 13th, 2017

4:30pm - 5:30pm in Room 2-136

Vsevolod Avrutsky (Moscow Institute of Physics and Technology)

Neural networks catching up with finite differences in solving partial differential equations in higher dimensions

Deep neural networks for solving Partial Differential Equations. From mathematical point of view neural network is a smooth function that depends on input vector as well as weights between its neurons, and all derivatives of the output with respect to input are available for analytical calculation. If we consider input as coordinates in 2D or 3D space we can "put" neural network into any differential equation and calculate residual error that we get for certain initial weights. By implementing an iterative process of training, we can further minimize that error and obtain a function that obeys the equation with necessary precision. This solution is defined not by values on a grid, but by parameters (weights) of a smooth function, thus it's easy to add or remove mesh points during solving. Method benefits from strong interpolating abilities of deep neural networks, and allows us to obtain solutions of linear and non linear PDEs with nearly machine precision in the whole region of space using very sparse grids. Future generalizations most likely will be able to solve equations in up to 6 dimensions.

## Schedule Spring 2017

**Meeting Time:** Wednesdays, 4:30-5:30pm, Room 2-136, Refreshments will be served in Room 2-290 at 3:30pm.

### Wednesday, February 22nd, 2017

4:30pm - 5:30pm in Room 2-136

Aminur Rahman (New Jersey Institute of Technology)

Dynamical modeling and analysis of walking droplets and chaotic logical circuits

The phenomena of wave-particle duality and logical circuits have been studied for over a hundred years. During the current century, scientists have been thinking differently about these systems with growing interest in studying hydrodynamic quantum analogs and chaotic logical circuits. Both walking droplets and chaotic logical circuits have mainly been modeled as differential and integro-differential equations. While many of these models are very good at reproducing the behavior observed in experiments, they can be quite difficult to analyze. We derive and analyze reduced models of these systems to investigate previously unobserved dynamical properties, which leads us to insights for the more complex models and even experiments. Specifically, we prove the existence of NeimarkSacker bifurcations and chaos, and show evidence of a new global bifurcation. In addition, we conduct novel experiments on circuits to show the viability of their design.

### Wednesday, March 8th, 2017

4:30pm - 5:30pm in Room 2-136

Ravi Samtaney (King Abdullah University of Science and Technology)

Discrete Exterior Calculus Discretization of the Navier-Stokes Equations

A conservative discretization of incompressible Navier-Stokes equations over surface simplicial meshes is developed using discrete exterior calculus (DEC). The DEC discretization is carried out for the exterior calculus form of Navier-Stokes equations, where the velocity field is represented by a 1-form. A distinguishing feature of our method is the use of an algebraic discretization of the interior product operator and a combinatorial discretization of the wedge product. Numerical experiments for flows over surfaces reveal a second order accuracy for the developed scheme for structured-triangular meshes, and first order accuracy for general unstructured meshes. The mimetic character of many of the DEC operators provides exact conservation of both mass and vorticity, in addition to superior kinetic energy conservation. The employment of various discrete Hodge star definitions based on both circumcentric and barycentric dual meshes is also demonstrated. The barycentric Hodge star allows the discretization to admit arbitrary simplicial meshes instead of being limited only to Delaunay meshes, as in previous DEC-based discretizations. The convergence order attained through the circumcentric Hodge operator is retained when using the barycentric Hodge. The discretization scheme is presented in detail along with numerical test cases demonstrating its numerical convergence and conservation properties. Preliminary results regarding the implementation of hybrid (circumcentric/barycentric) Hodge star operator are also presented. We conclude with some ideas for employing a similar method for magnetohydrodynamics.

### Wednesday, March 15th, 2017

4:30pm - 5:30pm in Room 2-136

Alexandria Volkening (Brown University)

Agent-based models of pattern formation on zebrafish

Zebrafish (Danio rerio) is a small fish with distinctive black and yellow stripes that form on its body and fins due to the interactions of different pigment cells. Working closely with the biological data, we present an agent-based model for these stripes that accounts for the migration, differentiation, an death of three types of pigment cells on the zebrafish body. We also explore stripe formation on the caudal fin, where bone rays and fin growth seem to help direct pigment cell placement. The development of both wild-type and mutated patterns will be discussed, as well as the non-local continuum limit associated with the model.

### Wednesday, April 5th, 2017

4:30pm - 5:30pm in Room 2-136

Ehud Yariv (Technion)

Drop electrohydrodynamics under strong electric fields

The leaky-dielectric electrohydrodynamic model was put forward by G. I. Taylor in the mid 1960's as an explanation to puzzling observations of drops deforming into an oblate spheroidal-type shape when subjected to a steady electric field. Taylor's theory neglects both fluid inertia and surface-charge convection. In addition, it assumes that the deformation from sphericity is small. These key assumptions are respectively tantamount to postulating that the Reynolds number, the electric Reynolds number and the capillary number are vanishingly small. Taken together, they eliminate the mutual decoupling between the electrostatic and flow problems, allowing for closed-form solutions where the fluid velocity scales as the square of the applied-field magnitude. Over the years, Taylor's work has been extended to situations where these numbers are asymptotically small (using regular perturbations) or even finite (using numerical simulations). The purpose of the present talk is to highlight the singular limits where the numbers become large, revealing non-conventional flow scaling and topologies.

### Wednesday, April 19th, 2017

4:30pm - 5:30pm in Room 2-136

Carlos Perez-Arancibia (MIT)

Windowed Green Function Method: An efficient high-order integral equation method for scattering in layered media

In this talk I will present a novel boundary integral equation method for the numerical solution of problems of scattering by obstacles and defects in the presence of layered media. This new approach, that we refer to as the windowed Green function (WGF) method, is based on use of smooth windowing functions and integral kernels that can be expressed directly in terms of the free-space Green function. The WGF method is fast, accurate, flexible and easy to implement. In particular straightforward modifications of existing (accelerated or unaccelerated) solvers suffice to incorporate the WGF capability. The mathematical basis of the method is simple: the method relies on a certain integral equation that is smoothly windowed by means of a low-rise windowing function, and is thus supported on the union of the obstacle and a small flat section of the interface between the two penetrable media. Various numerical experiments demonstrate that both the near- and far-field errors resulting from the proposed approach, decrease faster than any negative power of the window size. In some of those examples the proposed method is up to thousands of times faster, for a given accuracy, than an integral equation method based on use of the layer Green function and the numerical approximation of Sommerfeld integrals. Generalizations of the WGF method to problems of scattering by obstacles in layered media composed by any finite number of layers as well as wave propagation and radiation in open dielectric waveguides will also be discussed.

### Wednesday, April 26th, 2017

4:30pm - 5:30pm in Room 2-136

Sigel Gottlieb (University of Massachusetts, Dartmouth)

Error Inhibiting Block One-Step Schemes for Ordinary Differential Equations

The commonly used one step methods and linear multi-step methods all have a global error that is of the same order as the local truncation error. In fact, this is true of the entire class of general linear methods. In practice, this means that the order of the method is typically defined solely by order conditions which are derived by studying the local truncation error. In this work we investigate the interplay between the local truncation error and the global error, and develop a methodology which defines the construction of explicit error inhibiting block one-step methods (alternatively written as explicit general linear methods in Butcher's 1993 paper). These error inhibiting schemes are constructed so that the accumulation of the local truncation error over time is controlled, which results in a global error that is one order higher than the local truncation error. In this work, we delineate how to carefully choose the coefficient matrices so that the growth of the local truncation error is inhibited. We then use this theoretical understanding to construct several methods that have higher order global error than local truncation error, and demonstrate their enhanced order of accuracy on test cases. These methods demonstrate that the error inhibiting concept is realizable. Future work will further develop new error inhibiting methods and will analyze the computational efficiency and linear stability properties of these methods.

### Wednesday, May 3rd, 2017

4:30pm - 5:30pm in Room 2-136

Roberto Camassa (University of North Carolina at Chapel Hill)

Tailoring the tails in Taylor dispersion

The interplay between fluid flow and diffusion of a solute in the fluid is a primary mechanism for transport and mixing of substances, and one of the most ubiquitous phenomena in nature. Since the seminal investigation by G.I. Taylor in 1953, it has also been the focus of much mathematical efforts to model it. Taylor's counterintuitive result -- that at long times the effective diffusivity determined by the flow scales like the inverse of the tracer's molecular diffusivity -- is a classical illustration of mathematical analysis' predictive power and arguably one of the most remarkable effects demonstrated in this context. This talk will report results that focus on the interaction of advection and diffusion with fluid boundaries, such as pipes or ducts, at early and intermediate time scales in the transport process. This can have direct applications to microfluidics. Many microfluidic systemsincluding chemical reaction, sample analysis, separation, chemotaxis, and drug development and injectionrequire precision and control of solute transport. Although concentration levels are easily specified at injection, pressure-driven transport through channels is known to spread the initial distribution, resulting in reduced concentrations downstream. By monitoring the skewness (centered, normalized third moment) of the tracer distribution in laminar, pressure driven flows an unexpected phenomenon can be revealed: The channel's cross-sectional aspect ratio alone can control the shape of the concentration profile along the channel length. Thin channels (aspect ratio « 1) deliver solutes arriving with sharp fronts and tapering tails, whereas thick channels (aspect ratio ~ 1) produce the opposite effect. Thus, it is possible to deliver solute with prescribed distributions, ranging from gradual buildup to sudden delivery, based only on the channel dimensions.

### Wednesday, May 10th, 2017

4:30pm - 5:30pm in Room 2-136

Robert Pego (Carnegie Mellon University)

An analysis of merging-splitting group dynamics by Bernstein function theory (or: How to count fish using mathematics)

We study coagulation-fragmentation equations inspired by a simple model derived in fisheries science to explain data on the size distribution of schools of pelagic fish. Although the equations lack detailed balance and admit no H-theorem, we are able to develop a rather complete description of equilibrium profiles and large-time behavior, based on complex function theory for Bernstein and Pick functions. The generating function for discrete equilibrium profiles also generates the Fuss-Catalan numbers (derived by Lambert in 1758) that count all ternary trees with n nodes. The structure of equilibrium profiles and other related sequences is explained through a new and elegant characterization of the generating functions of completely monotone sequences as those Pick functions analytic and nonnegative on (-∞,1). This is joint work with Jian-Guo Liu and Pierre Degond.

## Schedule Fall 2016

**Meeting Time:** Wednesdays, 4:30-5:30pm, Room 2-136, Refreshments will be served in Room 2-290 at 3:30pm.

### Wednesday, September 21st, 2016

NOTE - DIFFERENT ROOM

4:30pm - 5:30pm in Room 2-105

William Henshaw (Rensselaer Polytechnic Institute)

Note: This seminar is joint with the Applied Mathematics Colloquium

Over-coming fluid-structure instabilities for incompressible flows and light bodies

The added-mass instability has, for decades, plagued partitioned fluid-structure interaction (FSI) simulations of incompressible flows coupled to light solids and structures. Many current approaches require tens or hundreds of expensive sub-iterations per time-step. In this talk some new stable partitioned algorithms are described for coupling incompressible flows with (1) compressible elastic bulk solids, (2) thin structural beams and (3) rigid bodies. These added-mass partitioned (AMP) schemes require no sub-iterations, can be made fully second- or higher-order accurate, and remain stable even in the presence of strong added-mass effects. These schemes are implemented using moving and deforming overlapping grids with the Overture framework.

### Wednesday, October 12th, 2016

4:30pm - 5:30pm in Room 2-136

Marc Hodes (Tufts University)

Effect of Thermocapillary Stress on Slip Length for a Channel Textured with Parallel Ridges

Lubrication of flows in microchannels enabled by textured and superhydrophobic surfaces has received much interest in recent years. In this seminar, we will compute the apparent hydrodynamic slip length, a measure of the effectiveness of lubrication in such flows, for (fully-developed and laminar) Poiseuille flow of liquid through a heated parallel-plate channel. One side of the channel is textured with parallel (streamwise) ridges and the opposite one is smooth. On the textured side of the channel, the liquid is the Cassie (unwetted) state, i.e., a lubricating layer of gas is trapped between menisci formed between ridges and the underlying substrate of the channel. No-slip and constant heat flux boundary conditions are imposed at the solid-liquid interfaces between the ridge tips and liquid and the menisci between ridges are considered flat and adiabatic. The smooth side of the channel is subjected to no slip and adiabatic boundary conditions. We account for the streamwise and spanwise thermocapillary stresses induced along menisci. When the latter are sufficiently small, Stokes flow may be assumed. Then, our solution is based upon a conformal map, albeit a cumbersome one. When, additionally, the ratio of channel height to ridge pitch is of order 1 or larger a less cumbersome, but equally accurate, solution is derived utilizing a matched asymptotic expansion. When inertial effects are relevant the slip length is numerically computed.

### Wednesday, October 26th, 2016

4:30pm - 5:30pm in Room 2-136

David Spivak (Massachusetts Institute of Technology)

Title: The Pixel Array method for nonlinear systems, and applications to numerical PDEs

Abstract: I will introduce the Pixel Array (PA) method, a fast and elementary technique for solving nonlinear systems of equations. Unlike quasi-Newton methods, the PA method is non-iterative and does not depend on finding Jacobians or on making an initial guess. It produces all solutions in a given bounding box, for any chosen subset of variables, and I will present test cases where it is faster than quasi-Newton methods (at least Julia's NLsolve). While the PA method may indicate a solution where none exists (i.e. introduce 'false positives'), it does not introduce false negatives.

The Pixel Array method is based on the slogan "plotting is a way of linearizing": given a plot--an array of boolean pixels--for each equation in a system, basic array multiplication yields the simultaneous solution to the whole system. I will explain these ideas, as well as give applications to solving PDEs numerically.

### Wednesday, November 30th, 2016

4:30pm - 5:30pm in Room 2-136

Jim Thomas (Courant Institute of Mathematical Sciences, NYU)

Title: Wave-vortex interactions in rotating shallow water

Abstract: The rotating shallow water is one of the simplest mathematical models that can capture features of rapidly rotating and strongly stratified geophysical flows, relevant to the atmosphere and the ocean. The dynamics of this fluid model is composed of fast waves and a slow vortical component. A relatively straightforward analysis of this PDE system reveals that the slow dynamics of the system, to leading order, is captured by the vortical field alone, which evolves unaffected by the waves. In other words, to a good degree of approximation, the waves do not affect the vortical field. This has remarkable consequences in this subfield of fluid dynamics, one of them being that to capture the slow dynamics of the system we may ignore the fast waves and simply track the vortical field with the help of a reduced model. Although such reduced models have been extremely popular in the field of geophysical flows, lack of waves affecting the vortical field continues to be a surprising fact, given that the governing equations are nonlinear. This has also limited our understanding of wave induced affects on the vortical fields. In this talk, I will present some recent developments on this problem. Asymptotic analysis will be used to derive reduced models that tries to capture the effect of fast waves on vortical fields. This will be followed up with high resolution numerical simulations of the PDE system to test the validity of the asymptotic models and to investigate features of those regimes where the asymptotic analysis does not formally apply.

## Schedule Fall 2015

**Meeting Time:** Wednesdays, 4:30-5:30pm, Room E17-133, Refreshments will be served in Room E17-401A at 3:30pm.

### Wednesday, November 4th, 2015

4:30pm - 5:30pm in Room E17-133

David Shirokoff (New Jersey Institute of Technology)

Convex relaxations in material science and approximate global minimizers

A wide range of material systems exhibit energy driven pattern formation governed by an underlying non-convex energy functional. Although numerically finding and verifying local minima to these functionals is relatively straight-forward, the computation and verification of global minimizers is much more difficult. Here the verification of global minimizers is often important in understanding the material phase diagram, especially at low temperatures. In this talk I will examine two separate model functionals: those arising in non-local pairwise interaction problems, and those containing a double-well potential. In the case of the pairwise interaction problems I will present an efficient systematic approach for the computation of approximate global minimizers based on a new convex relaxation and recovery technique. The approach is sometimes exact, and also provides a numerical recovery guarantee for the approximate minimizer that is often within a few percent of the global minimum. For the double-well functionals, I will derive sufficient conditions to show that a candidate minimizer is a global minimizer, and present a numerical method for estimating the order-disorder transition (ODT) curve.

### Monday, November 9th, 2015

SPECIAL TIME and PLACE

4:30pm - 5:30pm in Room E25-117

Mitchell Luskin (University of Minnesota)

Note: This seminar is joint with the Applied Mathematics Colloquiun

A Theoretical Examination of Diffusive Molecular Dynamics

We give a mathematical formulation of diffusive molecular dynamics. The idea of diffusive molecular dynamics is to first minimize an approximate free energy of the system with respect to the mean atomic coordinates (averaging over many vibrational periods), and then to perform a diffusive step where atoms and vacancies (or a solute in a solvent) flow on a diffusive time scale. One of the tools we make use of in this analysis is relative entropy, or the Kullback-Leibler divergence (KL). Joint with Gideon Simpson and David Srolovitz.

### Wednesday, December 2nd, 2015

4:30pm - 5:30pm in Room E17-133

Snezhana Abarzhi (Carnegie Mellon University)

Scale coupling in Richtmyer-Meshkov flows

We systematically study the Richtmyer-Meshkov instability (RMI) induced by strong shocks for fluids with contrasting densities and with small and large amplitude initial perturbations imposed at the fluid interface. The Smoothed particle hydrodynamics code (SPHC) is employed to ensure accurate shock capturing, interface tracking, and accounting for the dissipation processes. Simulations results achieve good agreement with existing experiments and with the rigorous theoretical analyses including zero-order theory describing the post-shock background motion of the fluids, linear theory providing RMI growth-rate in a broad range of the Mach and Atwood numbers, weakly nonlinear theory accounting for the effect of the initial perturbation amplitude on RMI growth-rate, and highly nonlinear theory describing evolution of RM bubble front.

We report the following findings: (1) The amount of energy that can be deposited by the shock at the RM-unstable interface is restrained. Significant part of the shock energy goes into compression and background motion of the fluids. (2) The initial amplitude is key factor of RMI evolution. The initial growth-rate of RMI is a non-monotone function of the initial amplitude. (3) At late times RM flow remains laminar rather than turbulent, and RM bubbles flatten and decelerate. In the fluid bulk, the dynamics at small-scale is heterogeneous, and is characterized by reverse cumulative jets, non-uniform velocity fields, local micro-structures, high pressure regions, and "hot spots". We show that the complicated character of scale coupling in RMI suggests new possibilities for flow mitigation and control.

## Schedule Spring 2015

**Meeting Time:** Wednesdays, 4:30-5:30pm,Room E17-129,Refreshments will be served in Room E17-401A at 3:30pm.

### Friday, May 29th, 2015

2:00pm - 3:00pm in Room E17-129

Benjamin Seibold (Temple University)

Benefits of Staggered Grids and Exponential Time Integrators in Radiation Moment Methods

For the numerical approximation of linear moment models of radiative transfer, we demonstrate the structural advantages that come with using staggered grids and exponential time integrators. These include: second-order accuracy; the commutation of discretization with manipulation of the PDE and discrete systems, respectively; and most importantly: asymptotic preserving properties, meaning that the diffusive nature of radiative transfer in the regime of large scattering is reproduced automatically by the numerical scheme, even if the small scales are not resolved by the computational grid. The numerical approach is implemented in the free and open-source software StaRMAP.

### Friday, May 15th, 2015

1:00pm - 2:00pm in E17-129

Matthias Taus Institute for Computational Engineering and Sciences

The University of Texas at Austin

Isogeometric Analysis for Boundary Integral Equations

Isogeometric analysis has emerged as a framework for integrating computational geometry and numerical methods. Isogeometric analysis makes use of interpolation functions widely used in computational geometry and adopts them as basis functions for numerical methods. This approach has been shown to have several superior properties over conventional numerical methods and has been applied to numerous problems of major practical and scientific significance.

The premise of this work is that isogeometric analysis is extremely beneficial for boundary element methods. This is mainly because isogeometric analysis involves smooth basis functions and exact representations of surfaces. We develop a sophisticated numerical integration scheme for boundary integral operators and apply isogeometric boundary element methods to Laplace's equation and the equations of linear elasticity. This allows us to show the superiority of isogeometric boundary element methods over conventional ones in several aspects. In particular, we develop boundary element patch tests, prove stability of collocation methods on C2 surfaces, represent singular and hyper-singular operators in terms of weakly singular ones, construct higher-order approximations without introducing additional degrees of freedom, and formulate integral equations leading to linear algebraic systems whose conditioning is independent of the mesh size.

### Wednesday, April 1st, 2015

4:30pm - 5:30pm in E17-129

Johnny Guzman Brown University

Higher-order finite element methods for elliptic problems with interface

We present higher-order piecewise continuous finite element methods for solving a class of interface problems in two dimensions. The method is based on correction terms added to the right-hand side in the standard variational formulation of the problem. We prove optimal error estimates of the methods on general quasi-uniform and shape regular meshes in maximum norms. In addition, we apply the method to a Stokes interface problem, adding correction terms for the velocity and the pressure, obtaining optimal convergence results. This is joint work with Manuel Sanchez-Uribe and Marcus Sarkis.

### Wednesday, April 8th, 2015

4:30pm - 5:30pm in E17-129

An Energy-Minimization Finite-Element Approach for Nematic Liquid Crystals

This talk outlines an energy-minimization finite-element approach to the computational modeling of equilibrium configurations for nematic liquid crystals under free elastic effects and electric effects. The method targets minimization of the system free energy based on the Frank-Oseen free-energy model. Solutions to the intermediate discretized linearizations are shown to exist generally and are unique under certain assumptions. This requires proving continuity, coercivity, and weak coercivity for the accompanying appropriate bilinear forms within a mixed finite-element framework. Error analysis demonstrates that the method constitutes a convergent scheme. Numerical experiments are performed for problems with a range of Frank elastic constants as well as simple and patterned boundary conditions. The resulting algorithm accurately handles heterogeneous Frank constants and efficiently resolves configurations resulting from complicated boundary conditions relevant in ongoing research. Additionally, flexoelectric effects are considered that have applications in bistable devices. If time, the talk will include current work on extending these results to dynamic systems that incorporate complex fluid effects.

### Wednesday, April 15th, 2015

4:30pm - 5:30pm in E17-129

Chi-Wang Shu Brown University

Discontinuous Galerkin method for hyperbolic equations with delta-singularities

Discontinuous Galerkin (DG) methods are finite element methods with features from high resolution finite difference and finite volume methodologies and are suitable for solving hyperbolic equations with nonsmooth solutions.

In this talk we will describe our recent work on the study of DG methods for solving hyperbolic equations with $\delta$-singularities in the initial condition, in the source term, or in the solutions. For such singular solutions, many numerical techniques rely on modifications with smooth kernels and hence may severely smear such singularities, leading to large errors in the approximation. On the other hand, the DG methods are based on weak formulations and can be designed directly to solve such problems without modifications, leading to very accurate results.

We will discuss both error estimates for model linearequations, in negative norm and in strong norm afterpost-processing, and applications to nonlinear systems including the rendez-vous systems and pressureless Euler equations involving $\delta$-singularities in their solutions. For the nonlinear case a high order accuracy bound-preserving limiter is crucial to maintain nonlinear stability and to avoid blowups of the numerical solution.

This is joint work with Yang Yang, Dongming Wei and Xiangxiong Zhang.

## Schedule Fall 2014

**Meeting Time:** Wednesdays, 4:30-5:30pm, Room E17-129, Refreshments will be served in Room E17-401A at 3:30pm.

### Wednesday, November 12th, 2014

4:30pm - 5:30pm in E17-129

Yossi Cohen Massachusetts Institute of Technology

**Note: This seminar is joint with the Physical Mathematics Seminar**

Where do rivers grow? Path selection and growth in a Harmonic Field

River networks exhibit a complex ramified structure that has inspired decades of studies. Yet, an understanding of the propagation of a single stream remains elusive. Here we invoke a criterion for path selection from fracture mechanics and apply it to the growth of streams in a diffusion field. We show that a stream will follow local symmetry in order to maximize the water flux and that its trajectory is defined by the local field in its vicinity. We also study the growth of a real network. We use this principle to reconstruct the history of a network and to find a growth law associated with it. The results show that the deterministic growth of a single channel based on its local environment can be used to characterize the structure of river networks.

## Schedule Spring 2014

**Meeting Time:** Wednesdays, 4:30-5:30pm, Room E17-128, Refreshments will be served in Room E17-401A at 3:30pm.

### Wednesday, April 23rd, 2014

4:30pm - 5:30pm in E17-128

Pavel Grinfeld Drexel University

The Moving Surface Analogs of Classical Equations of Applied Mathematics

The calculus of moving surfaces (CMS) is an extension of differential geometry to deforming manifolds. The fundamental equations of applied mathematics (the Laplace equation, the heat equation and the wave equation) find intriguing CMS equivalents, in which the surface itself is the unknown quantity. I will describe the fundamental elements of the CMS and illustrate a few of its many applications in differential geometry, shape optimization and dynamics of fluid films. Along the way, I will touch on a few interesting computational questions.

### Wednesday, February 19th, 2014

4:30pm - 5:30pm in E17-128

Dong Zhou Temple University

High-order methods for a pressure Poisson equation reformulation of the incompressible Navier-Stokes equations

Popular schemes for the incompressible Navier-Stokes equations (NSE), such as projection methods, are efficient but may introduce numerical boundary layers or have limited temporal accuracy due to their fractional step nature. Pressure Poisson equation (PPE) reformulations represent a class of methods that replace the incompressibility constraint by a Poisson equation for the pressure, with a suitable choice of the boundary condition so that the incompressibility is maintained. PPE reformulations of the NSE have important advantages: the pressure is no longer implicitly coupled to the velocity, thus can be directly recovered by solving a Poisson equation, and no numerical boundary layers are generated. In this talk we focus on numerical approaches of the Shirokoff-Rosales PPE reformulation. Interestingly, the "electric" boundary conditions provided for the velocity render classical nodal finite elements non-convergent. We thus present two alternative methodologies, mixed finite element methods and meshfree finite differences, and demonstrate that these approaches allow for arbitrary order of accuracy both in space and in time.

### Wednesday, March 5th, 2014

4:30pm - 5:30pm in E17-128

**Note: This seminar is joint with the Physical Mathematics Seminar**

Jamitons and the Predictive Accuracy of Macroscopic Traffic Models

Phenomenologically, the macroscopic (i.e., fluid-dynamical) description of vehicular traffic flow is further removed from reality than its microscopic (i.e., vehicle-based) description. It is therefore an important question why one should study macroscopic models. In this talk we demonstrate how the phenomenon of “phantom traffic jams” (the occurrence of traffic waves without any discernable cause) can be understood using an anology between traffic waves (“jamitons”) and detonation waves. It turns out that the macroscopic description provides fundamental insights into the dynamics of jamitons that goes beyond what is known from the microscopic modeling of the same phenomenon. In addition to the understanding of traffic waves, we demonstrate that second-order models can model inhomogeneities in traffic flow, while remaining within a macroscopic framework. Using real traffic data, we study the predictive accuracy of various data-fitted traffic models and show that the investment into more sophisticated models does pay off.

## Schedule Fall 2013

**Meeting Time:** Wednesdays, 4-5pm, Room E17-129, Refreshments will be served in Room E17-401A at 3:30pm.

### Tuesday, October 1st, 2013

2:30pm - 3:30pm in E51-149 (MIT-Tang Center)

Jon Wilkening University of California, Berkeley

** Note: This seminar is joint with the Physical Mathematics Seminar**

Traveling-Standing Water Waves and Microseisms

We study a two-parameter family of solutions of the surface Euler equations in which solutions return to a spatial translation of their initial condition at a later time. Pure standing waves and pure traveling waves emerge as special cases at fixed values of one of the parameters. We find many examples of wave crests that nearly sharpen to a corner, with corner angles close to 120 degrees near the traveling wave of greatest height, and close to 90 degrees for large-amplitude pure standing waves. However, aside from the traveling case, we do not believe any of these solutions approach a limiting extreme wave that forms a perfect corner.

We also compute nonlinear wave packets, or breathers, which can take the form of NLS-type solitary waves or counterpropagating wave trains of nearly equal wavelength. In the latter case, an interesting phenomenon occurs in which the pressure develops a large DC component that varies in time but not space, or at least varies slowly in space compared to the wavelength of the surface waves. These large-scale pressure zones can move very rapidly since they travel at the envelope speed, and may be partially responsible for microseisms, the background noise observed in earthquake seismographs.

### Wednesday, October 16th, 2013

**Esteban Tabak** Courant Institute of Mathematical Sciences, New York University

** Note: This seminar is joint with the Physical Mathematics Seminar**

Conservation Law Modeling of Mixing in Layered Hydrostatic Flows

This talk proposes a framework to quantify the large-scale effects of fluid mixing without resolving the associated small scale motion. The equations of motion for hydrostatic flows adopt the form of a hyperbolic system of nonlinear equations that yields breaking waves. To model the shock waves that ensue, one needs to involve integral conserved quantities, such as mass and momentum. Yet in a system composed of layers that may mix, first physical principles do not provide a set of conserved quantities large enough to completely determine the flow. Our proposal is to replace the conventional conservations laws of each layer's mass and momentum, invalid after shocks form, by others, such as energy, in a way that provides a natural description of the mixing process.

On the numerical side, we develop a finite volume algorithm to simulate general systems of conservation laws with arbitrary conserved quantities. Applications range from the lock-exchange, through the mixing Rossby adjustment problem, to overturning circulations such as the atmospheric Hadley cells.

Jointly with Paul A. Milewski (Bath) and Robert E. Friel (Courant).

## Schedule Spring 2013

**Meeting Time:** Wednesdays, Wednesdays, 4-5pm, Room 2-147, Refreshments will be served in Room 2-290 at 3:30pm.

### Wednesday, February 13th, 2013

4-5pm in 2-147

**David Shirokoff** McGill University

An Active Penalty Method for the Incompressible Navier-Stokes Equations

The volume penalty method provides a simple, efficient approach for solving the incompressible Navier-Stokes equations in domains with boundaries or in the presence of moving objects. Despite the simplicity, the method suffers from poor convergence in the penalty parameter, thereby restricting accuracy of any numerical method. We demonstrate that one may achieve high order accuracy by introducing an active penalty term. We discuss how to construct the modified penalty term, and then provide 2D numerical examples demonstrating improved convergence for the heat equation and Navier-Stokes equations. In addition, we show that modifying the penalty term does not significantly alter the time step restriction from that of the conventional penalty method.

This work is joint with J-C. Nave.

### Wednesday, March 13th, 2013

4-5pm in 2-147

Benjamin Seibold Temple University

Jet Schemes and Gradient-Augmented Level Set Methods

Jet schemes are semi-Lagrangian advection approaches that evolve parts of the jet of the solution (i.e., function values and higher derivatives) along characteristic curves. Suitable Hermite interpolations give rise to methods that are high order accurate, yet optimally local, i.e., the update for the data at any grid point uses information from a single grid cell only. Jet schemes can be systematically derived from an evolve-and-project methodology in function spaces, which in particular yields stability estimates. We present a comparison of the accuracy and computational cost of jet schemes with WENO and Discontinuous Galerkin schemes.

For interface evolution problems, jet schemes give rise to gradient-augmented level set methods (GALSM). These possess sub-grid resolution and yield accurate curvature approximations. We demonstrate how the optimal locality of jet schemes gives rise to a straightforward combination with adaptive mesh refinement (AMR), and provide an outlook on jet schemes for nonlinear Hamilton-Jacobi equations.

### Wednesday, March 20th, 2013

4-5pm in 2-147

John A. Evans University of Texas at Austin

Structure-Preserving B-spline Methods for the Incompressible Navier-Stokes Equations

The incompressible Navier-Stokes equations are infused with important geometric structure, evidenced by a wide array of balance laws for momentum, energy, vorticity, enstrophy, and helicity. These balance laws are considered to be of prime importance in the evolution of laminar and turbulent flow structures, and they are even believed to play a role in the regularity of Navier-Stokes solutions. The key to unlocking much of the geometric structure of incompressible flow is precisely its volume-preserving nature, yet most numerical methods only satisfy the incompressibility constraint in an approximate sense. Consequently, such methods do not obey many fundamental laws of physics and produce unphysical results for many flow configurations of interest.

In this talk, I will discuss recently developed B-spline discretizations for viscous incompressible fluid flows which satisfy the incompressibility constraint pointwise. As incompressibility is satisfied pointwise, these discretizations replicate the geometric structure of the Navier-Stokes equations and properly balance momentum, energy, vorticity, enstrophy, and helicity. Furthermore, these discretizations provably converge to physically-relevant weak solutions of the Navier-Stokes equations which satisfy a local energy balance in space-time. The above attributes in conjunction with the geometric flexibility and resolution properties of B-splines make these discretizations an attractive candidate for reliable numerical simulation of flows on complex engineering geometries, and indeed theoretical and numerical analyses have revealed that these discretizations are more accurate and stable than classical schemes. In this talk, I will give a general overview of structure-preserving B-spline methods and present numerical results illustrating the promise of this new CFD technology in the context of low-speed flows.

### Wednesday, April 3d, 2013

4-5pm in 2-147

Jeff Aristoff Numerica Corporation

**Note: This seminar is joint with the Physical Mathematics Seminar**

Implicit Runge-Kutta Methods: Are They Practical?

A large class of methods for propagating the state of a system and its associated uncertainty (i.e., its probability density function) require the propagation of an ensemble of particles or states through nonlinear dynamics. In this talk, a new implicit Runge-Kutta-based approach for the efficient propagation of uncertainty will be described. Comparison to existing explicit-based approaches will be made in the context of space situational awareness, wherein uncertainty propagation is critical for space catalog maintenance, collision avoidance, and object characterization.

### Wednesday, April 10th, 2013

4-5pm in 2-147

Charles S. Peskin Courant, NYU

Cardiac Fluid Mechanics and Electrophysiology in a Unified Mathematical and Computational Framework

The heart is a fluid-mechanical system that is coordinated and controlled by electrical activity intrinsic to the heart itself. The immersed boundary (IB) method was introduced to study fluid-structure interaction of the heart valves, but has since been elaborated into a modeling framework for the heart as a whole, including blood flow in the cardiac chambers, elasticity of the flexible heart valve leaflets, and the active elasticity of the muscular heart walls. The IB method employs a fixed Cartesian grid for the storage of the velocity and pressure field of the entire system, and a moving collection of fibers that cut through the Cartesian grid and serve to model both the collagen fibers of the valve leaflets and the muscle fibers of the myocardium. The bidomain equations of cardiac electrophysiology separately track the spatially distributed extracellular and intracellular voltages and currents, which are coupled by capacitive and ionic transmembrane currents. The IB formulation of the bidomain equations mirrors that of the mechanical IB method. In particular, the electrical IB method employs a fixed Cartesian grid (which extends beyond the myocardium into the electrically conducting blood and surrounding tissue) for the extracellular (and extracardiac) electrical variables, and a moving system of fibers to keep track of intracellular voltages and currents, as well as membrane associated variables, in a Lagrangian manner. Both in the mechanical and also in the electrical IB method, a regularized version of the Dirac delta function is used to communicate between Eulerian and Lagrangian variables. In particular, it is used to transfer mechanical force and electrical transmembrane current from the moving cardiac fibers to the fixed Cartesian grid in which they are immersed, and to evaluate fluid velocity and extracellular voltage at the present locations of the cardiac fibers. In this way, we achieve a unified model that captures both fluid-structure interaction and electro-mechanical interaction of the beating heart.

This work is joint work with Boyce E. Griffith and David M. McQueen

### Monday, May 6th, 2013

SPECIAL TIME AND PLACE

3PM in room 2-132

Ganesh Mahadevan Colorado School of Mines

**Note: This seminar is joint with the Imaging and Computing Seminar**

Navier-Stokes equations on rotating surfaces: Regularity, algorithm, analysis, and simulation

Development of efficient algorithms with rigorous analysis for partial differential equations (PDE) on domains and surfaces requires knowledge of the regularity of solutions of the PDE. In this talk, for the Navier-Stokes equations on rotating spheres we discuss (i) the global regularity for real and complexified time; (ii) a high-order algorithm with stability and convergence analysis; (iii) long time simulation of a benchmark atmospheric flow model and justification of a turbulence theory.

## Schedule Fall 2012

**Meeting Time:** Wednesdays, 4-5pm, Room 2-135, Refreshments will be served in Room 2-290 at 3:30pm.

### Thursday, September 20, 2012

SPECIAL DAY / TIME / LOCATION

4-5pm, Room 32-124

Bill Henshaw Lawrence Livermore National Laboratory

**Jointly with MIT Distinguished Speaker Series in Computational Science and Engineering **

Solving Fluid Structure Interaction Problems on Overlapping Grids

This talk will discuss the numerical solution of fluid structure interaction (FSI) problems on overlapping grids. Overlapping grids are an efficient, flexible and accurate way to represent complex, possibly moving, geometry using a collection of overlapping structured grids. For incompressible flows with moving geometry we have been developing a scheme that is based on an approximated-factored compact scheme for the momentum equations together with a multigrid solver for the pressure equation. The overall scheme is fourth-order accurate in space and (currently) second-order accurate in time. The scheme will be described and results will be presented for some three-dimensional (parallel) computations of flows with moving rigid-bodies. In recent work, we have also been developing an FSI scheme based on the use of deforming composite grids (DCG). In the DCG approach, moving boundary-fitted grids are used near the deforming solid interfaces and these overlap non-moving grids which cover the majority of the domain. For compressible flows and elastic solids we have derived a new interface projection scheme, based on the solution of a fluid-solid Riemann problem, that overcomes the well known ``added-mass'' instability for light solids. The FSI-DCG approach is described and validated for some fluid structure problems involving high speed compressible flows coupled to rigid and elastic solids. The interesting case of a shock hitting an ellipse of zero mass is also presented.

### Wednesday, October 17th, 2012

4-5pm in 2-135

Homer F. Walker W.P.I.

Anderson acceleration for fixed-point iterations

Fixed-point iterations occur naturally and are commonly used in a broad variety of computational science and engineering applications. In practice, fixed-point iterates often converge undesirably slowly, if at all, and procedures for accelerating the convergence are desirable. This talk will focus on a particular acceleration method that originated in work of D. G. Anderson [J. Assoc. Comput. Machinery, 12 (1965), 547-560] and has been independently re-invented on several occasions. This method has enjoyed considerable success in a few applications (notably in electronic-structure computations, where it is known as Anderson mixing) but seems to have been untried or underexploited in many other important applications. Moreover, while other acceleration methods have been extensively studied by mathematicians and numerical analysts, Anderson acceleration has received relatively little attention from them until recently, despite there being many significant unanswered mathematical questions. In this talk, I will outline Anderson acceleration, discuss some of its theoretical properties, and demonstrate its performance in several PDE applications. This work is joint in part with P. Ni and in part with P. A. Lott, C. S. Woodward, and U. M. Yang.

### Wednesday, November 14th, 2012

4-5pm in 2-135

**Morris R. Flynn** U. Alberta

Buckling of a thin, viscous film in an axisymmetric geometry

By adapting the Foppl-von Karman equation, which describes the deformation of a thin elastic membrane, we present an analysis of the buckling pattern of a thin, very viscous fluid layer subject to shear in an axisymmetric geometry. A linear stability analysis yields a differential eigenvalue problem, whose solution, obtained using spectral techniques, yields the most unstable azimuthal wave-number, m*. Contrary to the discussion of Slim et al. (J. Fluid Mech., 694, pp. 5-28, 2012), it is argued that the axisymmetric problem shares the same degeneracy as its rectilinear counterpart, i.e. at the onset of instability, m* is indefinitely large. Away from this point, however, a comparison with analogue experimental results is both possible and generally favorable. In this vein, we describe the laboratory apparatus used to make new measurements of m*, the phase speed and the wave amplitude; note that no prediction concerning the latter two quantities can be made using the present theory. Experiments reveal a limited range of angular velocities wherein waves of either small or large amplitude may be excited. Transition from one to the other regime does not appear to be associated with a change in m*.

### Wednesday, November 28th, 2012

4-5pm in 2-135

Kushal Kedia MIT

An immersed-body method coupled with AMR for chemically reactive flows over heat conducting solid objects

In this talk I will discuss the development of a state-of-the-art semi-implicit numerical approach to simulate multiscale phenomena in energy conversion processes in which the coupling of chemical kinetics and molecular transport in the gas and solid phase, and the flow dynamics plays an important role. In collaboration with Sandia National Laboratories, we are developing a low-Mach number, reactive flow code for the direct simulation.of chemically reactive flows interacting with solid objects. All the relevant spatial and temporal scales are resolved incorporating detailed molecular transport and complex chemical kinetics. Adaptive Mesh Refinement (AMR) is used to achieve high resolution in regions of sharp gradients. Operator-splitting is used to separate temporal scales, such as radicals diffusion, relatively slow convection and thermal diffusion, providing a fast and efficient way to treat the stiff chemistry. Conjugate-heat exchange between embedded reactive solid objects and the reactive gas flow is incorporated. Solid walls are treated using a marker-function tracking and buffer layers to impose the correct boundary conditions. Fully resolved simulations provide a detailed insight into complex fundamental combustion processes such as flame stabilization near flame-holders, blowoff and resonant dynamics.

### Wednesday, December 12th, 2012

4-5pm in 2-135

Ugur G. Abdulla Florida Institute of Technology

On the Optimal Control of the Free Boundary Problems for the Second Order Parabolic Equations

In this talk, we describe a new variational formulation of the inverse Stefan problem, where information on the heat flux on the fixed boundary is missing and must be found along with the temperature and free boundary. This type of inverse problem arose, in particular, in a bioengineering problem on the laser ablation of biological tissues. We employ optimal control framework, where boundary heat flux and free boundary are components of the control vector, and optimality criteria consists of the minimization of the sum of $L_2$-norm declinations from the available measurement of the temperature flux on the fixed boundary and available information on the phase transition temperature on the free boundary. This approach allows one to tackle situations when the phase transition temperature is not known explicitly, and is available through measurement with possible error. It also allows for the development of iterative numerical methods of least computational cost due to the fact that for every given control vector, the parabolic PDE is solved in a fixed region instead of full free boundary problem. We prove well-posedness in Sobolev spaces framework and convergence of discrete optimal control problems to the original problem both with respect to cost functional and control.

## Schedule Spring 2012

**Announcement: This semester the NumPDE seminar is being run within the Physical Math Seminar (PMS) on:**

Tuesdays, 2:30-3:30pm, Room 2-135, Refreshments will be served in Room 2-290 at 3:30pm.

For the up-to-date schedule please check the PMS website.

### Thursday, February 2nd, 2012

SPECIAL TIME AND LOCATION

Room 2-145 from 2:30 PM - 3:30 PM

Cristina Turner University of Cordoba, FaMAF, Argentina

Shape optimization for tumor location

In non-invasive thermal diagnostics, accurate correlations between the thermal image at skin surface and interior human physiology are desired. In this work an estimation methodology to determine unknown geometrical parameters of an embedded tumor is proposed. We defined a functional that represents the mismatch between a measured experimental temperature profile, which may be obtained by infrared thermography on the skin surface, and the solution of an appropriate boundary problem. This functional is related to the geometrical parameters through the solution of the boundary problem, in such a way that finding the minimum of this functional form also means finding the unknown geometrical parameters of the embedded tumor. Sensitivity analysis techniques coupled with the adjoint method were considered to compute the shape derivative of the functional and a nonmonotone spectral projected gradient method was implemented to solve the optimization problem of finding the optimal geometric parameters.

## Schedule Fall 2011

**Announcement: This semester the NumPDE seminar is being run within the Physical Math Seminar (PMS) on:**

Tuesdays, 2:30-3:30pm, Room 2-105, Refreshments will be served in Room 2-290 at 3:30pm.

For the up-to-date schedule please check the PMS website.

## Schedule Spring 2011

### Wednesday, February 16th, 2011

4-5pm in 2-131

Enkeleida Lushi Courant, NYU

An efficient numerical method for multiple surfactant-covered interfaces in 2D Stokes Flow

I will present an efficient and highly accurate numerical method for computing the deformation of many surfactant-coated interfaces in 2D slow viscous flow. Surfactants locally alter the surface tension in interfaces and thus change the nature of their motion. We consider bubbles as well as drops covered with interface-bound or insoluble surfactant. The advection-diffusion equation for the surfactant concentration on the interface is coupled with the Stokes equations in the fluid domain through a Laplace-Young condition. The Stokes equations are first recast as an integral equation based on the complex variable theory of the biharmonic equation. The numerical method employed is spectrally accurate and uses a fast-multipole accelerated iterative procedure. The computational cost per time step is only O (N logN ) operations, with N being the number of discretization points. The interface is described by a spectral mesh and is advected according to the fluid velocity in such a manner so as to preserve equal arc length spacing of marker points. This equal arc length framework has the dual advantage of dynamically maintaining the spatial mesh and allowing efficient, implicit treatment of the stiffest terms in the dynamics. Several phenomenologically-different examples will be presented, e.g. tip streaming, bursting bubbles, interacting bubbles, etc. Last, I will talk about the numerical method extension to study an interface with soluble surfactants.

### Wednesday, March 16, 2011

4-5pm in 2-131

Burt Tilley WPI

Macroscale field effects on flows between closely spaced, thin laminates.

Fluid flows through anisotropic media are found in a wide variety of geophysical and biological systems. The macroscale behavior of these systems depend on the microstructure, which in turn may depend on local and global physical processes. Classically, geometric restrictions are needed to model these systems on the largest length scale, and we are interested in developing effective models which relax these restrictions. To explore the development of these multiscale models, we present two fundamental problems.

In the first problem, we consider an array of closely spaced, purely dielectric rigid laminates with nonuniform thickness. The laminate thickness and spacing varies on a length scale much longer than the characteristic thickness of the laminates. In the spacing between the laminates, an electrically conducting fluid is driven by an applied electric field through electroosmosis and electrophoresis along with an applied pressure gradient. Debye layers occur at the laminate-fluid interface, which are assumed to be much smaller than the laminate thickness. From a modification of the classical homogenization approach that relies on a physical microscale constraint in place of a geometric constraint, we derive an effective set of equations that describe the fluid pressure, the anion and cation concentrations in the fluid and the electric potential. Anisotropic dispersion effects in the electric field are included, and electroneutrality in the fluid is not imposed. We find that gradients in the laminate spacing can lead to charge accumulation when electroosmosis and the electrophoresis induced from the anisotropic dispersion effects balance.

Our second problem centers on the flow of an incompressible fluid that saturates an array of deformable laminates with gravity acting in the spanwise direction. The aspect ratio of the characteristic spacing between the laminates is much smaller than the characteristic scale of the laminate length, and we use this aspect ratio to find effective equations for the components of the stress tensor of the effective material. Momentum conservation on the macroscale at leading order result in a coupled set of elliptic equations for the unknown laminate spacing, orientation and the local pore pressure. At next order, effective stress-strain relations are derived based on the local geometry and material properties of the laminates and fluid. This work is done in collaboration with B. Vernescu and J.D. Plummer at WPI.

### Wednesday, March 30, 2011

4-5pm in 2-131

Hans Johnston UMASS Amherst

Numerical Study of Turbulent Thermal Convection Between Conditions of Constant Temperature and Constant Flux

We report the results of high resolution direct numerical simulations of two-dimensional Rayleigh-Benard convection for Rayleigh numbers up to Ra=10^10 in order to study the influence of temperature boundary conditions on turbulent heat transport. Specifically, we consider the extreme cases of fixed heat flux (where the top and bottom boundaries are poor thermal conductors) and fixed temperature (perfectly conducting boundaries). Both cases display identical heat transport at high Rayleigh numbers fitting a power law Nu~0.138xRa^(.285) with a scaling exponent indistinguishable from 2/7 = 0.2857... above Ra = 10^7. The overall flow dynamics for both scenarios, in particular the time averaged temperature profiles, are also indistinguishable at the highest Rayleigh numbers. The findings are compared and contrasted with results of recent three-dimensional simulations.

This is joint work with C.R. Doering (U. Michigan), Cheng Wang (UMass-Dartmouth), Jian-Guo Liu (Duke)

### Wednesday, April 6, 2011

4-5pm in 2-131

Jian-Guo Liu Duke Universtiy

Dynamics of orientational alignment and phase transition

Phase transition of directional field appears in some physical and biological systems such as ferromagnetism near Curie temperature, flocking dynamics near critical mass of self propelled particles. Dynamics of orientational alignment associated with the phase transition can be effectively described by a mean field kinetic equation. The natural free energy of the kinetic equation is non-convex with a minimum level set consisting of a sphere at super-critical case, a typical spontaneous symmetry breaking behavior in physics. In this talk, I will present some analytical results on this dynamics equation of orientational alignment and exponential convergence rate to the equilibria for both supper and sub critical cases, as well at algebraic convergence rate at the critical case. A new entropy and spontaneous symmetry breaking analysis played an important role in our analysis.

### Wednesday, April 13, 2011

4-5pm in 2-131

Michael Minion University of North Carolina at Chapel Hill

Novel Methods for Temporal Integration of PDEs

I will discuss ongoing research on the development of novel methods for the temporal integration of ODEs and PDEs. The strategy employed in the numerical methods is based on an iterative deferred corrections approach, within which one can utilize operator splitting, multirate timestepping, and parallelization in the temporal direction to achieve better computational efficiency. Much of the recent work is motivated by fluid-structure interaction problems in biological systems and I will discuss the difficulties and progress related to some target applications in this area. Recent results concerning parallelization in the temporal direction will also be presented.

### Wednesday, May 4, 2011

4-5pm in 2-131

Hyoungsu Baek MIT

Accuracy and stability enhancements of semi-implicit schemes for the Navier–Stokes equations through sub-iteration

We present an iterative semi-implicit scheme for the incompressible Navier–Stokes equations, which is stable at CFL numbers well above the nominal limit. We have implemented this scheme in conjunction with spectral discretizations, which suffer from serious time step limitations at very high resolution. Specifically, at each time level, the nonlinear convective term and the pressure boundary condition – both of which are treated explicitly in time – are updated using fixed-point iteration and Aitken relaxation. Instabilities from pressure boundary condition and convective term are suppressed by sub-iteration and relaxation of flow field, respectively. Eigenvalue analysis of Stokes problem confirms that this scheme is unconditionally stable. For Navier-Stokes flows, the proper value of the relaxation parameter is obtained as a function of the flow parameters, which leads to second- and third-order temporal accuracy at CFL number 5–14. Moreover, stability of the scheme does not depend on the order of polynomial basis employed in the spectral element method. Systematic accuracy, stability, and cost comparisons are presented against the standard semi-implicit method and a fully-implicit scheme that does not require Newton’s iterations.

## Schedule Fall 2010

### Wednesday, November 3, 2010

4-5pm in 2-132

Paul Chesler (MIT - Dept. of Physics)

Applied string theory: from gravitational collapse to heavy ion collisions

A remarkable result from heavy ion collisions at the Relativistic Heavy Ion Collider is that shortly after a collision, the medium produced behaves as a nearly ideal liquid. The system is very dynamic and evolves from a state of two colliding nuclei to a liquid in a time roughly equivalent to the time it takes light to cross a proton. Understanding the mechanisms behind the rapid approach to a liquid state is a challenging task. In recent years string theory has emerged as a powerful tool to study non-equilibrium phenomena, mapping the (challenging) dynamics of quantum systems onto the dynamics of classical gravitational systems. The creation of a liquid in a quantum theory maps onto the classical process of gravitational collapse and black hole formation. I will describe how one can use techniques borrowed from numerical relativity in astrophysics to study processes which mimic the dynamics of heavy ion collisions.

### Friday, November 5, 2010

SPECIAL DAY / TIME / LOCATION

4-5pm in 2-135

Snezhana I. Abarzhi (University of Chicago)

Turbulent mixing and beyond: problems, concepts, solutions

Turbulent mixing plays an important role in a broad variety of natural and artificial systems, spanning astrophysical to atomistic scales and low to high energy densities. Examples include inertial confinement fusion, supernovae, stellar convection, non-canonical boundary layers and optical free-space communications. Theoretical description of non-equilibrium mixing transports is a challenging problem due to singular aspects of the governing (Euler or Navier-Stokes) equations. Furthermore these processes are statistically unsteady and their fluctuating quantities are essentially time-dependent and non-Gaussian. We apply the new theoretical concept, the rate of momentum loss, to describe the transports of mass, momentum and energy in turbulent mixing flow and to capture its anisotropic and inhomogeneous character. It is shown that invariant, scaling and spectral properties of unsteady turbulent mixing differ substantially from those of isotropic and homogeneous turbulence. Time- and scale-invariance of the rate of momentum loss leads to non-dissipative momentum transfer, to 1/2 and 3/2 power-law scale-dependencies of the velocity and Reynolds number and to non-Kolmogorov spectra. Turbulent mixing exhibits more order compared to isotropic turbulence, and its viscous and dissipation scales are finite and set by flow acceleration. We suggest how to describe the random character of the unsteady turbulent flow and show that the rate of momentum loss is the statistic invariant and a robust diagnostic parameter for either sustained or time-dependent acceleration. Some criteria are outlined for the estimate of the fidelity and information capacity of the experimental and numerical data sets.

### Wednesday, November 17, 2010

4-5pm in 2-132

Paul Hand (MIT - Mathematics)

An Adaptive Multiscale Model for Simulating Cardiac Conduction

We present an adaptive multi-scale numerical method for simulating cardiac action potential propagation along a single strand of heart muscle cells. This method combines macroscale cable partial differential equations posed over the tissue with different microscale equations posed over discrete cellular geometry. The microscale equations are applied near a wavefront, where electric potential changes rapidly, and the macroscale equations are applied everywhere else. One-dimensional numerical simulations reveal that the adaptive scheme accurately reproduces the action potential waveforms and wavespeeds of the full microscale model. In particular the adaptive simulations capture the ephaptic effect that has been demonstrated by several mathematical models the literature. Our study demonstrates that capturing the behavior of action potential propagation may sometimes require simulations that are adaptive and multi-scale because (1) fully resolving all cells in a whole-heart simulation is computationally infeasable; (2) one-dimensional simulations at coarseness levels that are affordable for whole-heart simulations lead to substantial errors in conduction speed; and (3) the macroscopic cable equations are invalid at low levels of gap junctional coupling. More abstractly, this work provides an example of when and how to selectively couple the dynamics of a periodic microstructure to those of its macroscale description.

### Wednesday, December 1, 2010

4-5pm in 2-132

David Shirokoff (MIT - Mathematics)

An efficient method for the incompressible Navier-Stokes equations on irregular domains with no-slip boundary conditions, high order up to the boundary.

Common efficient schemes for the incompressible Navier-Stokes equations, such as projection or fractional step methods, have limited temporal accuracy as a result of matrix splitting errors, or introduce errors near the domain boundaries (resulting in weakly convergent solutions). In this paper we recast the incompressible Navier-Stokes equations (with the velocity prescribed at the boundary) as an equivalent system, for the primary variables velocity and pressure. We do this in the usual way away from the boundaries, by replacing the incompressibility condition on the velocity by a Poisson equation for the pressure. The key difference from the usual approaches occurs at the boundaries, where we use boundary conditions that unequivocally allow the pressure to be recovered from knowledge of the velocity at any fixed time --- avoiding the common difficulty of an, apparently, over-determined Poisson problem. Since in this alternative formulation the pressure can be accurately and efficiently recovered from the velocity, the recast equations are ideal for numerical marching methods. This new system can be discretized using a variety of methods, in principle to any desired order of accuracy. In this work we illustrate the approach with a 2-D second order finite difference scheme on a Cartesian grid, and devise an algorithm to solve the equations on domains with curved (non-conforming) boundaries, including a case with a non-trivial topology (a circular obstruction inside the domain). This algorithm achieves second order accuracy in the $L^{\infty}\/$ norm, for both the velocity and the pressure. The scheme has a natural extension to 3-D.

### Friday, December 10, 2010

SPECIAL DAY / TIME / LOCATION

4-5pm in 2-103

Benjamin Seibold (Temple University)

A new perspective on the moment closure problem in radiative transfer

Radiative transfer can be modeled by a kinetic equation that describes the evolution of the particle density function in a phase space of time, position, and the angle of flight (possibly more). While the direct simulation of this high dimensional mesoscopic equation is possible (albeit very costly), in many applications it is desirable to have a description in terms of macroscopic equations. An expansion in the angular variable yields an equivalent system of infinitely many macroscopic moment equations. The fundamental question how to best truncate this system is the moment closure problem. Many types of closure strategies exist. These are typically based on an asymptotic arguments or assume higher moments be quasi-stationary.

In this talk, we present a very different approach to derive moment closures, based on the Mori-Zwanzig formalism of irreversible statistical mechanics. Here, the influence of the truncated moments on the resolved moments is modeled by a memory term. Moment closures are then defined by a choice of a probability measure on the phase space and suitable approximations to the memory term. We demonstrate that existing closures, such as PN, SPN, and diffusion correction closures, can be re-derived with this formalism. In addition, new closures arise, such as the crescendo-diffusion closure.

### Wednesday, December 15, 2010

4-5pm in 2-132

Alex Marques (MIT - Mathematics)

High-Order Solution of the Poisson Equation with Interface Jump Conditions – Towards a Generalized Ghost Fluid Method

The Poisson equation with jump discontinuities across an interface is of central importance in multiphase diffusion phenomena. In particular, we must solve this equation to obtain the pressure distribution in multiphase flows. The difficulty in numerically solving this equation arises due to the discontinuous nature of the solution around the interface. Over the last three decades we have seen the development of several methods to address this issue. However, obtaining high order of accuracy in a systematic fashion still poses great challenges to the numerical fluid mechanics community.

In this presentation we shall introduce a new method to obtain high- order solutions to the Poisson Equation with interface jump conditions. The method is based on local smooth extrapolations of the solution field across the interface. The extrapolation procedure uses a local functional basis that satisfies the interface conditions to high order of accuracy. The overall approach is compatible with standard discretizations of the Laplace operator, and leads to minimal modifications to the final linear system to be inverted. In particular, we apply the method with cubic Hermit polynomials to develop fourth-order accurate schemes. As a result, standard Poisson solvers can be used with only minimal modifications. Details of the method and applications will be presented.

## Schedule Spring 2010

### Wednesday, March 17, 2010

4-5pm in 2-105

Maria Cameron (Courant Institute of Mathematical Science)

Computing transition paths for rare events

The overdamped Langevin equation is often used as a model in molecular dynamics. At low temperatures, a system evolving according to such an SDE spends most of the time near the potential minima and performs rare transitions between them. A number of methods have been developed to study the most likely transition paths. I will focus on one of them: the MaxFlux functional. The MaxFlux functional has been around for almost thirty years but not widely used because it is challenging to minimize. Its minimizer provides a path along which the reactive flux is maximal at a given finite temperature. I will show two ways to derive it in the framework of transition path theory: the lower bound approach and the geometrical approach. I will present an efficient way to minimize the MaxFlux functional numerically. I will demonstrate its application to the problem of finding the most likely transition paths in the Lennard-Jones-38 cluster between the face-centered-cubic and icosahedral structures.

### Thursday, April 1, 2010

4:00 PM in Building 2, Room 105

Daniela Tordella (Politecnico di Torino, Italy)

Transients of Three-dimensional Perturbations and the Role of Long Waves in the Plane Wake. Relationships with Turbulence

Results from an exploratory analysis of the transient and long-term behavior of small three-dimensional perturbations in the circular cylinder wake are presented. A few comments on the role of long waves in this stability problem and on the role of small scales in turbulence (anisotropy, possible long-term interactions) are also proposed.

### Wednesday, April 21, 2010

4:00 PM in Building 2, Room 105

John F. Gibson (Georgia Institute of Technology)

Invariant Solutions and State-space Dynamics of Low-Re Turbulence

It has recently become possible to compute precise 3D, nonlinear solutions of Navier-Stokes equations at Reynolds numbers above the onset of turbulence, for simple geometries such as pipes and channels. These solutions capture the form and dynamics of "coherent structures" and provide a starting point for understanding low-Reynolds turbulence as a dynamical system. In this talk I will present a number of equilibrium, traveling wave, and periodic orbit solutions of plane Couette flow, emphasizing visualizations of their physical structure and state-space dynamics, and comparisons to turbulent flow. Certain spatially localized solutions exhibit homoclinic snaking remarkably similar to that observed in simpler 1D PDE systems. What emerges is a picture of low-Reynolds turbulence as a walk among a set of weakly unstable invariant solutions.

### Wednesday, April 28, 2010

4:00 PM in Building 2, Room 105

Professor Anthony T Patera (Department of Mechanical Engineering, MIT)

Certified Reduced Basis Methods; Application to Continuum Mechanics and Transport

We discuss reduced basis approximation and associated a posteriori error estimation for reliable and rapid solution of parametrized partial differential equations.

The crucial ingredients are rapidly convergent Galerkin approximations over a space spanned by "snapshots" on the parametrically induced solution manifold; effective constructions for stability-constant lower bounds; rigorous and sharp a posteriori error bounds for the outputs/quantities of interest; efficient POD (in time)/Greedy (in parameter) selection of quasi-optimal samples; and Offline-Online computational procedures for very rapid response in the real-time and many-query contexts.

We consider the application of these techniques to illustrative problems from heat transfer, solid mechanics, acoustics, and fluid dynamics. We also briefly describe recent implementations on deployed platforms.

### Wednesday, May 5, 2010

4:00 PM in Building 2, Room 105

Kenneth Kamrin (Harvard University)

A Novel Finite-Difference Method for Large-Deformation Solid Mechanics

A classic problem in computational continuum mechanics is to simulate the deformation of a solid-like material to the point of very large strain. Standard approaches such as the finite-element method can lose accuracy when elements become overly distorted, and while there have been a number of mesh-repairing methods considered in the past, fully Eulerian algorithms provide a natural way to handle this issue. This talk develops and tests a fixed-grid approach to solid deformation based on a kinematic variable we call the "reference map". The method can be used in a variety of common circumstance cared about in solid mechanics: statics, quasi-statics, smooth dynamics, and shock propagation. Moreover, by the similarity of this method to finite-difference CFD algorithms, the approach also leads to a basic method for computing fluid/structure interactions on a single fixed grid.

## Schedule Fall 2009

### Wednesday, September 23, 2009

4-5pm in 2-132

**Peyré** (CNRS and Université Paris-Dauphine)

Sparse Processing of Images

In this talk, I will review recent work on the sparse representations of natural images. I will focus on the application of these emerging models for the resolution of various imaging problems, which include compression, denoising and super-resolution of images, as well as compressive sensing and compressive wave computations. Natural images exhibit a wide range of geometric regularities, such as curvilinear edges and oscillating textures. Adaptive image representations select bases from a dictionary of orthogonal or redundant frames that are parameterized by the geometry of the image. If the geometry is well estimated, the image is sparsely represented by only a few atoms in this dictionary. The resolution of ill-posed inverse problems in image processing is then regularized using sparsity constraints in these adapted representations.

### Monday, October 6th, 2009

SPECIAL DAY / TIME / LOCATION

4:30-5:30pm in 4-370

Kai Schneider (Universite de Provence, Marseille, France)

**Jointly with the Applied Math Colloquium**

Adaptive Space-Time Multiresolution Techniques for Nonlinear PDEs

We present efficient fully adaptive numerical schemes for evolutionary partial differential equations based on a finite volume (FV) discretization with explicit time discretization. A multiresolution strategy allows local grid refinement while controlling the approximation error in space. The costly fluxes are evaluated on the adaptive grid only. For time discretization we use an explicit Runge-Kutta scheme of second-order with a scale-dependent time step. On the finest scale the size of the time step is imposed by the stability condition of the explicit scheme. On larger scales, the time step can be increased without violating the stability requirement of the explicit scheme. Embedded Runge-Kutta methods of second and third order are then used to choose automatically the new time step while controlling the approximation error in time. Non-admissible choices of the time step are avoided by limiting its variation.

The implementation of the multiresolution representation uses a dynamic tree data structure, which allows memory compression and CPU time reduction. This new numerical scheme is validated using different classical test problems in one, two and three space dimensions. The gain in memory and CPU time with respect to the finite volume scheme on a regular grid is reported, which demonstrates the efficiency of the new method.

This work is joint work with M. Domingues, S. Gomes and O. Roussel.

### Wednesday, October 14, 2009

4:30-5:30pm in 2-132

Markus Schmuck (MIT - Chem. Engr.)

Modeling, Analysis, and Numerics Of The Navier-Stokes-Nernst-Planck-Poisson System

We introduce a macroscopic model which allows to describe the essential electrokinetic phenomena as electrophoresis and -osmosis. Then, we present the basic analytical results for the Navier-Stokes-Nernst-Planck-Poisson system. Next, we propose and analyze two convergent nite element discretizations which preserve all char- acteristic properties of weak solutions in the discrete setting. We begin with a scheme based on perturbation and then show how we can improve the properties and consistency of the scheme by using a suitable truncation. In the last part of the talk, we derive eective macroscopic properties of a porous solid-electrolyte composite. This leads to a huge dimensional reduction by upscaling the microstructure. The results are gained by the formal multiple-scale method in the context of homogenization.

### Wednesday, November 4, 2009

4:30-5:30pm in 2-132

Qi Qi Wang (MIT - Aero./Astro.)

Solving adjoint equations for unsteady fluid flows

Methods based on solving adjoint equations are widely used in trajectory optimization, shape design, inverse methods, optimal control and uncertainty quantification. The adjoint equation solves the sensitivity of an objective function with respect to the governing equations and its parameters. The adjoint sensitivity gradient can be used in approximating the objective function and in gradient based optimization.

The adjoint equations for unsteady incompressible Navier-Stokes equations are derived and solve. Efficient backwards time integration of the adjoint equation is performed using the dynamic checkpointing scheme. The adjoint solution reveals the sensitivity of the objective function with respect to both the geometry and velocity of solid bodies in the flow field. Applications of the adjoint solution in optimization and uncertainty quantification are discussed.

### Wednesday, November 18, 2009

4:30-5:30pm in 2-132

Leslie Greengard (NYU - CIMS)

The Nonuniform FFT, Heat Flow, and Magnetic Resonance Imaging Reconstruction

The nonuniform FFT arises is a variety of applications, from medical imaging to radio astronomy to the numerical solution of partial differential equations. In a typical problem, one is given an irregular sampling of N data points in the frequency domain with the goal of reconstructing the corresponding function at N points in the physical domain. When the sampling is uniform, the fast Fourier transform (FFT) allows this calculation to be carried out in O(N log N) operations. Unfortunately, when the data is nonuniform, the FFT does not apply. In the last few years, a number of algorithms have been developed which overcome this limitation and are often referred to as nonuniform FFTs. In this talk, we describe the basic algorithm and some of its applications.

## Schedule Spring 2009

Poster Announcement Spring 2009

### Wednesday, February 18, 2009

4-5pm in 2-146

Martin Frank (Department of Mathematics, University of Kaiserslautern)

Optimal treatment planning in radiotherapy based on Boltzmann transport equations

Treatment with high energy ionizing radiation is one of the main methods in modern cancer therapy that is in clinical use. During the last decades, two main approaches to dose calculation were used, Monte Carlo simulations and semi-empirical models based on Fermi-Eyges theory. A third way to dose calculation has only recently attracted attention in the medical physics community. This approach is based on the deterministic kinetic equations of radiative transfer.

In this work, we present a Boltzmann transport model for dose calculation in radiation therapy. We discuss simplifications of this model and additionally formulate an optimal control problem for the desired dose. Based on this formulation, we derive optimality conditions. Numerical results in one and two dimensions are presented.

### Wednesday, February 25, 2009

4-5pm in 2-146

Shing Yu Leung (Department of Mathematics, UC Irvine)

A grid based particle method for the evolution of open curves and surfaces

We present in this talk a new numerical method for modeling motion of open curves in two dimensions and open surfaces in three dimensions. Following the grid based particle method we have recently proposed, we represent the open curve or the open surface by meshless Lagrangian particles sampled according to an underlying fixed Eulerian mesh. The underlying grid is used to provide a quasi-uniform sampling and neighboring information for meshless particles. The key idea in the talk is to represent and to track end-points of the open curve and boundary points of the open surface explicitly and consistently with interior particles. We apply our algorithms to several applications including spiral crystal growth modeling and image segmentation using active contours.

### Wednesday, March 11, 2009

4-5pm in 2-146

Roland Bouffanais (Department of Mechanical Engineering, MIT)

Simulation of unsteady transitional swirling flow with a moving free surface using a spectral element method

Unsteady incompressible viscous flows of a fluid enclosed in a cylindrical container with an open top surface are discussed. These moving free-surface flows are generated by the steady rotation of the solid bottom end-wall. Such type of flows belongs to a group of recirculating lid-driven cavity flows with geometrical axisymmetry. The top surface of the cylindrical cavity is left open so that the free surface can freely deformed. The Reynolds regime corresponds to unsteady transitional flows with some incursions in the fully laminar regime. The approach taken here revealed new nonaxisymmetric flow states that are investigated based on a fully three-dimensional solution of the Navier-Stokes equations for the free-surface cylindrical swirling flow, without resorting to any symmetry property unlike all other results available in the literature. Theses solutions are obtained through direct numerical simulations based on a highly-accurate Legendre spectral element method combined with a moving-grid technique.

### Monday, April 6, 2009

4:30-5:30pm in 4-237

John Lowengrub (Department of Mathematics, UC Irvine)

**Joint seminar with Applied Mathematics Colloquium**

Multi-scale models of solid tumor growth and angiogenesis

We present and investigate models for solid tumor growth that incorporate features of the tumor microenvironment including tumor-induced angiogenesis. Using analysis and nonlinear numerical simulations, we explore the effects of the interaction between the genetic characteristics of the tumor and the tumor microenvironment on the resulting tumor progression and morphology. We account for variable cell-cell/cell-matrix adhesion in response to micro-environmental conditions (e.g. hypoxia) and to the presence of multiple tumor cell species. We focus on glioblastoma and quantify the interdependence of the tumor mass on the microenvironment and on the cellular phenotypes. The model provides resolution at various tissue physical scales, including the microvasculature, and quantifies functional links of molecular factors to phenotype that for the most part can only be tentatively established through laboratory or clinical observation. This allows observable properties of a tumor (e.g. morphology) to be used to both understand the underlying cellular physiology and to predict subsequent growth or treatment outcome, thereby providing a bridge between observable, morphologic properties of the tumor and its prognosis.

### Wednesday, April 8, 2009

4-5pm in 2-146

Jan Hesthaven (Division of Applied Mathematics, Brown University)

Discontinuous Galerkin methods for the modeling of free surface flows using high-order Boussinesq approximations

We shall discuss the modeling of free surface flows and fluid-structure interactions using high-order Boussinesq approximations. These sets of equations are characterized by being purely dispersive and strongly non-linear, with additional complications introduced by high-order spatial derivatives and cross-derivatives. The key elements of the formulation and some of the properties of the Boussinesq system will be discussed in some detail. This shall be used to argue why discontinuous Galerkin methods may be a suitable approach for the solution of these equations. We shall develop the basic elements required for solving this system, discuss a number of subtleties and address practical concerns of performance and efficient solvers. The computational approach will be extensively validated with both benchmark tests and experimental data.

This is work done in collaboration with A.P. Engsig-Karup (DTU, Denmark), P. Madsen (DTU, Denmark), H. Bingham (DTU, Denmark) and T. Warburton (Rice).

### Wednesday, April 22, 2009

4-5pm in 2-146

Pavel Grinfeld (Department of Mathematics, Drexel University)

Hamiltonian dynamic equations for fluid films

Two dimensional models for hydrodynamic systems, such as soap films, have been studied for hundreds of years. Yet there has not existed a fully nonlinear system of dynamic equations analogous to the classical Euler equations. We propose an exact nonlinear system for the dynamics of a fluid film. The system is derived in the classical Hamiltonian framework and neither the velocities nor the deviation from the equilibrium are assumed small. We discuss the properties of the proposed equations and results of numerical simulations.

### Wednesday, April 29, 2009

4-5pm in 2-146

Frederic Gibou (Department of Mechanical Engineering, UC Santa Barbara)

Sharp interface methods for moving boundary problems

In this talk, we will describe recent (and less recent) advances in the field of free boundary problems, with an emphasis on level-set and ghost-fluid methods. We will discuss their applications to two-phase flows, the Stefan problem and image guided surgery. Novel methods for adaptive mesh refinement will be discussed as well.

### Wednesday, May 6, 2009

4-5pm in 2-146

Jacob White (Department of Electrical Engineering and Computer Science, MIT)

Subtleties associated with biochemical oscillator sensitivity analysis in the presence of conservation laws

Sensitivity analysis is often dismissed as a trivial application of Taylor series, even when implicitly defined operators such as differential equations are involved. The subject maybe a bit more subtle when the Jacobian of the implicitly defined operator is singular, and the resolution of this subtlety can be quite problem dependent. In this talk we consider oscillators described by mass action kinetics, review the known material on sensitivity analysis for oscillators posed as two-point boundary value problems, and present the additional subtlety associated with the presence of conservation laws in the oscillator. In particular, we show that generalized eigenvectors play a surprising (at least to the authors) and essential role.

This is joint work with Jared Toettcher, Anya Castillo, Paul Barton and Bruce Tidor.

### Wednesday, May 13, 2009

New time: 2:30-3:30pm in 2-136

4-5pm in 2-146

Alexander Vladimirsky (Department of Mathematics, Cornell University)

Causality and efficiency: Non-iterative numerical methods

Our knowledge of the direction of information flow is fundamental for many efficient numerical methods (e.g., time-marching for evolutionary PDEs). However, for many problems (including first-order static nonlinear PDEs) the direction of information flow might be a priori unknown even if it is otherwise well-defined. This leads to a common use of iterative methods, which can be unnecessarily inefficient.

For certain systems of nonlinear equations, the "causality" present in the problem can be used to uncover the direction of information flow at runtime. Exploiting causality to effectively de-couple nonlinear systems is the fundamental idea behind Dijkstra's classical method for finding shortest paths on graphs.

We will use a continuous analogue of this principle to build efficient methods for a wide class of causal problems. We will illustrate this approach using examples from:

- continuous and hybrid optimal control (e.g., optimal traveling on foot and using the buses);
- multi-objective optimal control (e.g., quickest paths constrained by maximum pathlength & shortest paths constrained by maximum travel time);
- anisotropic front propagation (e.g., first-arrivals and multiple-arrivals in seismic imaging);
- optimal control under uncertainty (e.g., optimal traveling when the map or the terminal time is not quite known);
- Markov decision processes (e.g., stochastic shortest paths on graphs);
- dynamical systems (e.g., approximation of "geometrically stiff" invariant manifolds of vector fields).

## Schedule Fall 2008

### Monday, September 28, 2008

4:30-5:30pm in 2-231

Andrea L. Bertozzi (Department of Mathematics, University of California Los Angeles)

**Joint seminar with Applied Mathematics Colloquium**

Swarming by nature and by design

The cohesive movement of a biological population is a commonly observed natural phenomenon. With the advent of platforms of unmanned vehicles, this occurrence is attracting renewed interest from the engineering community. This talk will review recent research results on both modeling and analysis of biological swarms and also design ideas for efficient algorithms to control groups of autonomous agents. For biological models we consider two kinds of systems: driven particle systems based on force laws and continuum models based on kinematic rules. Both models involve long-range social attraction and short range dispersal and yield patterns involving clumping, mill vortices, and surface-tension-like effects. For artificial platforms we consider the problem of boundary tracking of an environmental material and consider both computer models and demonstrations on real platforms of robotic vehicles. We also consider the motion of vehicles using artificial potentials.

### Wednesday, October 1, 2008

4-5pm in 2-151

Raúl A. Radovitzky (MIT Aeronautics and Astronautics, MIT Institute for Soldier Nanotechnologies)

Discontinuous Galerkin methods applied to fracture mechanics

In this talk, I will present our progress in exploring the numerical formulation of problems in the mechanics of solid materials within the framework of discontinuous Galerkin (DG) methods. The interest in DG methods stems from its potential to address some limitations in conventional finite element formulations of problems involving complex phenomena including fracture and nonlocal material response.

Included in this talk will be the weak formulation of the boundary value problem of finite deformation elasticity, the discretized problem and its numerical properties, and the numerical implementation of the method within a conventional finite element framework. I will also discuss the extension to problems involving dynamic plastic deformations and the scalability in the case of explicit time integration.

Finally, I will discuss the application to problems involving fracture. In this respect, we discuss the connections of the method to the well-established Cohesive Zone Model (CZM) approach to crack nucleation and propagation. We show that the latter method describes the elasticity in the cohesive zone in an intrinsically inconsistent way, whereas the DG approach is inherently consistent.

Time permitting I will describe the application of DG methods to higher order problems, including shell theory and nonlocal isotropic elasticity.

### Wednesday, October 15, 2008

4-5pm in 2-151

Chi-Wang Shu (Division of Applied Mathematics, Brown University)

High order methods for convection dominated PDEs - An overview

In this talk we will give an overview of algorithm development and application, with an emphasis on recent progress, on high order methods for convection dominated partial differential equations. We will mainly discuss the finite difference weighted essentially non-oscillatory (WENO) schemes, finite volume WENO schemes, and discontinuous Galerkin (DG) finite element methods. A comparison of their relevant advantages and disadvantages will be given.

### Wednesday, October 29, 2008

4-5pm in 2-151

Benjamin Seibold (MIT, Applied Mathematics)

Particle methods - Sparsity and exact conservation

In many numerical approaches that operate on a fixed grid, a proper treatment of convection poses the largest challenge. Particle methods are a way out. Particles are moved with the flow, and thus take care of convection automatically. This advantage comes at a price: The governing equations have to be discretized on point cloud, and particle management (merging upon particle collision, and insertion of particles into holes) is required.

On the aspect of discretization, we consider meshfree finite difference approximations of the Laplace operator. A fundamental problem is how to select a small number of neighbors to a point, but guarantee the stability of the arising stencil. We present an approach, based on linear optimization, that selects optimally sparse finite difference stencils, while guaranteeing stability.

On the aspect of particle management, we present a new approach that approximates scalar 1D conservation laws only by particle motion and particle merging and insertion. The method conserved area exactly, and in addition conserves entropy exactly when no shocks are present. We outline an application that is relying heavily on such exact conservation properties.

### Wednesday, November 12, 2008

4-5pm in 2-151

Sigal Gottlieb (UMass Dartmouth, Mathematics)

Time stepping methods for numerical solution of hyperbolic PDEs with shocks

When numerically solving a hyperbolic conservation law it is important to consider the properties of the spatial discretization combined with the time discretization. If the problem is smooth, it is sufficient to linearize the problem and analyze the L2 stability properties of the resulting discretization. However, if the solution is nonsmooth, stability in the L2 norm is not sufficient. This is because for PDEs with discontinuous solutions, the presence of oscillations prevents the approximation from converging uniformly. To ensure that the method does not allow oscillations to form, we require stability of the nonlinear system in the maximum norm or in the TV semi-norm. Strong stability preserving (SSP) high order time discretizations were developed to ensure these types of nonlinear stability properties SSP methods preserve the strong stability properties -- in any norm, seminorm or convex functional -- of the spatial discretization coupled with first order Euler time stepping. I will describe the development of SSP methods and the current state-of-the-art of these methods.

### Wednesday, November 19, 2008

4-5pm in 2-151

John E. Dolbow (Duke University, Civil and Environmental Engineering)

An embedded interface finite element method: Application to modeling stimulus-responsive hydrogels

Stimulus-responsive hydrogels (SRHs) are macromolecular polymer networks immersed in a solvent, synthesized to exhibit large volumetric swelling in response to small changes in environmental stimuli. For example, SRHs have been designed to actuate in response to changes in temperature, solvent concentration, pH, and light. The unique properties of these ``soft-wet" materials make them appealing for a large range of applications. They have been used for sensors to detect trace contaminants in fuel lines and autonomous control in microfluidic systems, just to name a few. However, a lack of understanding into the relationships between gel composition, kinetics, and mechanical response has hindered designs based on SRHs and delayed the technological transfer from the laboratory to the marketplace.

This presentation will focus on our recent efforts to characterize the unique behavior of SRHs and develop robust finite element models of the same. We have developed continuum-based models for chemically and thermally-induced volume transitions in hydrogels (Dolbow et al., 2004, 2005). Consistent with experimental observations, the models allow for a sharp interface separating swelled and collapsed phases in SRHs. The models predict characteristic swelling times that are proportional to the square of the characteristic linear dimension of the specimen. Our results have also suggested several synthetic pathways that might be pursued to engineer hydrogels with optimal response times.

We discretize our models using an embedded interface finite element method, wherein the interface geometry is allowed to be independent of the underlying mesh (Dolbow and Franca, 2008; Dolbow and Harari, 2008). The method enhances the approximation basis in the vicinity of the interface to capture discontinuities in both primary and secondary fields, and is applicable to fully unstructured quad and tet meshes. Interfacial constraints are enforced weakly using a modification of a classical, variationally consistent, interior penalization method. Numerical tests indicate the accuracy of the method to be competitive with widely used finite-difference methods for elliptic interface problems.

### Wednesday, December 3, 2008

4-5pm in 2-151

Alain Karma (Northeastern University, Physics Department and Center for Interdisciplinary Research on Complex Systems)

Phase-field modeling of fracture

The phase-field method has emerged as a powerful computational tool to describe the complex evolution of interfaces in a wide range of materials science problems ranging from alloy solidification to solid-state precipitation to thin-film patterning and dislocation dynamics. The chief advantage of this method is to avoid front-tracking by the introduction of a scalar order parameter that varies smoothly in space, thereby making the interface between two phases or two crystal grains spatially diffuse. Evolution equations for these order parameters are derived variationally from a Lyapounov functional that represents the total free-energy of the system.

This talk will give an overview of the more recent application of this method to fracture where the phase-field order parameter is used as a local measure of damage. The method allows to describe both the short-scale physics of failure inside the microscopic process zone around the crack tip and macroscopic linear elasticity within a self-consistent set of partial differential equations that can be readily simulated. In addition, those equations can be analyzed in certain limits to derive, rather than to guess, laws of crack motion. Results shed light on the origin and validity of the several decade old principle of local symmetry used to determine crack paths under various loading conditions and on its extension to anisotropic materials.