Schedule for: 19w5189 - Women In Numerical Methods for PDEs and their Applications

A buffet dinner is served daily between 5:30pm and 7:30pm in the Vistas Dining Room, the top floor of the Sally Borden Building.

Breakfast is served daily between 7 and 9am in the Vistas Dining Room, the top floor of the Sally Borden Building.

A brief introduction to BIRS with important logistical information, technology instruction, and opportunity for participants to ask questions.

For many applications in science and engineering, the ability to efficiently and accurately approximate solutions to elliptic PDEs dictates what physical phenomena can be simulated numerically. In this talk, we present a high-order accurate discretization technique for variable coefficient PDEs with smooth coefficients. The technique comes with a nested dissection inspired direct solver that scales linearly or nearly linearly with respect to the number of unknowns. Unlike the application of nested dissection methods to classic discretization techniques, the constant prefactors do not grow with the order of the discretization. The discretization is robust even for problems with highly oscillatory solutions. For example, a problem 100 wavelengths in size can be solved to 9 digits of accuracy with 3.7 million unknowns on a desktop computer. The precomputation of the direct solver takes 6 minutes on a desktop computer. Then applying the computed solver takes 3 seconds. Since the method is naturally domain decomposing, a simple parallel implementation reduces the time for precomputation to 30 seconds. Applications of the algorithm to inverse scattering will also will be presented.

We study numerical solutions for parabolic equations with highly varying (multiscale) coefficients. Such equations typically appear when modelling heat diffusion in heterogeneous media like composite materials. For these problems classical polynomial based finite element methods fail to approximate the solution well unless the mesh width resolves the variations in the data. This leads to issues with computational cost and available memory, which calls for new approaches and methods. In this talk I will present a multiscale method based on localized orthogonal decomposition, first introduced by M\r{a}lqvist and Peterseim (2014). The focus will be on how to generalize this method to time dependent problems of parabolic type.

In seismic exploration, sources and measurements of seismic waves on the surface are used to determine model parameters representing geophysical properties of the earth. Full-waveform inversion (FWI) is a nonlinear seismic inverse technique that inverts the model parameters by minimizing the difference between the synthetic data from the forward wave propagation and the observed true data in PDE-constrained optimization. The traditional least-squares method of measuring this difference suffers from three main drawbacks including local minima trapping, sensitivity to noise, and difficulties in reconstruction below reflecting layers. Unlike the local amplitude comparison in the least-squares method, the quadratic Wasserstein distance from the optimal transport theory considers both the amplitude differences and the phase mismatches when measuring data misfit. I will briefly review our earlier development and analysis of optimal transport-based inversion and include improvements, for example, a stronger convexity proof. The main focus will be on the third "challenge" with new results on sub-reflection recovery.

Numerical analysis provides a wealth of theoretical results for time-dependent PDEs, however inverse problems pose their own difficulties that are not immediately solved using black-box solutions. This talk explores topics on the relation between backward solves vs forward solves in the time-dependent optimization loop of an unsteady problem. Special interest will be assigned to efficiency preservation in the context of analytic vs numerical treatment of the adjoint of a nonlinear unsteady PDE. Current approaches at scale in two R\&D100 codes (PETSc and Nek5000) will be presented, together with relevant use cases. Also we will outline issues of solvability of ill-posed inverse problems as well as open questions.

Lunch is served daily between 11:30am and 1:30pm in the Vistas Dining Room, the top floor of the Sally Borden Building.

Meet in the Corbett Hall Lounge for a guided tour of The Banff Centre campus.

Meet in foyer of TCPL to participate in the BIRS group photo. The photograph will be taken outdoors, so dress appropriately for the weather. Please don't be late, or you might not be in the official group photo!

Large scale numerical models of systems evolving over a wide range of temporal and spatial scales are routinely encountered in a variety of fields from fluid mechanics and plasma physics to weather prediction and chemical engineering. Many of such, so-called stiff, systems present a computational challenge and fuel continuous need to improve the fidelity, robustness and efficiency of numerical time integrators. Over the past decades, exponential integration emerged as a numerical technique that carries significant computational savings. In this talk we will explain advantages exponential methods offer and discuss theoretical and practical aspects of designing and implementing different classes of efficient exponential integrators. We will illustrate performance gains these schemes provide using test problems and examples from several applications in plasma physics and computer graphics.

I will present new model for micro-swimmers that takes into account the counter-rotation of the body and flagella, as seen in motile bacteria or spermatozoa. The disturbance fluid flow of one such swimmer now contains a torque-dipole singularity in addition to the leading order force-dipole singularity. This head-and-flagella counter-rotation gives rise to clock-wise circling at no-slip walls just as observed in experiments of bacteria on surfaces. I will discuss the scattering behavior of spermatozoa in a forest of cylindrical pillars, confirmed also by new experiments. Last, we show large scale and fast simulations of thousands of such swimmers that interact with each-other, surfaces, as well as immersing fluid.

Computational modeling can be used to reveal insights into the mechanisms and potential side effects of a new drug. Here we will focus on diabetes, which affects 1 in 10 people in North America. In particular, we are interested in a class of relatively novel drug treatment, the SGLT2 inhibitors (sodium-glucose co-transporter 2 inhibitors). E.g., Dapagliflozin, Canagliflozin, and Empagliflozin. We conduct simulations to better understand any side effect these drugs may have on our kidneys (which are the targets of SGLT2 inhibitors). Interestingly, these drugs may have both positive and negative side effects.

A buffet dinner is served daily between 5:30pm and 7:30pm in the Vistas Dining Room, the top floor of the Sally Borden Building.

We present new second order limiters for the discontinuous Galerkin method applied to hyperbolic conservation laws. We prove that under a suitable time step restriction, the limiters enforce the local maximum principle for linear and nonlinear scalar equations on unstructured triangular meshes. We also show that this time step size constraint is tight. We discuss under what restrictions the limiters are second order accurate and enforce the local bound. We then extend these limiters to nonconforming meshes that arise in adaptive computations. Limiting on such meshes is particularly difficult due to the lack of structure and variability in element sizes. We will discuss how limiters differ in implementation complexity, computational cost, and accuracy of computed solutions.

We develop high-order multiscale discontinuous Galerkin (DG) methods for one-dimensional stationary Schr\"{o}dinger equations with oscillating solutions. We propose two types of multiscale finite element spaces, and prove that the resulting DG methods converge optimally with respect to the mesh size $h$ in $L^2$ norm when $h$ is small enough. In the lowest order case, we prove that the second order multiscale DG method has the optimal convergence even when the mesh size is larger than the wave length. Numerically we observe that all these multiscale DG methods have at least the second-order convergence on coarse meshes and optimal high-order convergence on fine meshes.

This work was initiated because we wanted to simulate the biomechanics of the respiratory system. The main muscle that drives the respiration is the diaphragm which has an aspect ratio of approximately 1:100 of the length and thickness scales. There are several practical challenges to deal with. Creating a smooth representation of the geometry extracted from medical images; generating anisotropic scattered nodes within the thin volume; applying physically relevant boundary conditions; and solving the (non-linear) elasticity equations. In this talk, I will show results for the geometry representation as well as analysis of and results for a simplified anisotropic linear elasticity problem.

Understanding the behavior of complex fluids in biology presents mathematical, modeling, and computational challenges not encountered in classical fluid mechanics, particularly in the case of fluids with large elastic forces that interact with immersed elastic structures. I will describe recent work on micro-organism locomotion in viscoelastic fluids that highlights some of these challenges and discuss the specific modeling and numerical considerations needed for these types of problems.

Many biological fluids, like mucus and cytoplasm, have prominent viscoelastic properties, which lead to the subdiffusive behavior of immersed particles. We propose a viscoelastic generalization of the Landau-Lifschitz Navier-Stokes fluid model for particles that are passively advected by such a medium and develop a simulation techniques based on the theory of stationary Gaussian processes. In contrast to the stochastic immersed boundary method for viscous fluids, which relies on step-by-step simulation techniques exploiting the Markov property, our method is based on the numerical evaluation of the covariance associated with individual fluid modes. The numerical method is spectral, meshless and uses results from the simulations of Generalized Langevin Equations.The implementation presents many practical problems, mostly stemming from the fact that the physical regime of interest corresponds to a situation where the memory kernel has a very slow (power law) decay.

Swimming bacteria with helical flagella are self-propelled micro-swimmers in nature, and the swimming strategies of such bacteria vary depending on the number and the position of flagella on the cell body. In this talk, we will introduce two microorganisms, multi-flagellated E. coli and single-flagellated Vibrio A. The Kirchhoff rod theory is used to model the elastic helical flagellum and the penalty method is employed to describe the dynamics of the rigid cell body. The hydrodynamic interaction between the fluid and the cell is represented by the regularized Stokes formulation. The focus of the talk will be on how bacteria reorient swimming direction.

Cilia are hair-like appendages attached to microorganisms that allow the organisms to traverse their fluid environment. The algae Volvox are spherical swimmers with potentially thousands of individual cilia on their surface and the cilia coordination is poorly understood. In this work, we developed a mathematical model of cilia on the outer surface of a sphere submerged in a fluid. The goal was to determine if factors related to the spherical shape affected cilia synchronization. We modeled each beating cilium tip as a small rotating sphere, attached to a point just above the spherical surface by a spring. This was achieved by using a regularized image system for Stokes flow outside of a sphere. Previous models showed synchronization when cilia were attached to a sphere but this was largely because the cilia were beating in the same direction all the way around the sphere. It is known that Volvox cilia beat toward one pole, where some cilia are beating in opposite direction, which somehow helps with directed motion and rotation. By including more biologically realistic assumptions about ciliary beating in our model, we were able to simulate and understand how groups of cilia synchronize to aid in directed motion.

Conventional computational methods often create a dilemma for fluid-structure interaction problems. Typically, solids are simulated using a Lagrangian approach with a grid that moves with the material, whereas fluids are simulated using an Eulerian approach with a fixed spatial grid, requiring some type of interfacial coupling between the two different perspectives. Here, a fully Eulerian method for simulating structures immersed in a fluid will be presented. By introducing a reference map variable to model finite-deformation constitutive relations in the structures on the same grid as the fluid, the interfacial coupling problem is highly simplified. A key challenge in this technique is to extrapolate field values to new grid points as the solid moves across the fixed grid. We develop and test a least-squares regression-based algorithm that is more robust and suitable for parallel computation than the existing approaches. The reference map technique is particularly well suited for simulating soft, highly-deformable materials and many-body contact problems, which are common in nature, experiments, as well as industrial settings. We demonstrate the method by 3D simulations using a large-scale, three-dimensional, parallel implementation of it.

SUNDIALS is a suite of robust and scalable solvers for systems of ordinary differential equations, differential-algebraic equations, and nonlinear equations designed for use on computing systems ranging from desktop machines to super computers. The suite consists of six packages: CVODE(S), ARKode, IDA(S), and KINSOL, each built on common vector and solver application programming interfaces (API) allowing for application-specific and user-defined linear solvers, nonlinear solvers, data structures, encapsulated parallelism, and algorithmic flexibility. In this presentation we will overview the capabilities of the SUNDIALS suite, give a preview of new developments, and highlight applications of SUNDIALS in some DOE simulations. This work was performed under the auspices of the U.S. Department of Energy by Lawrence Livermore National Laboratory under Contract DE-AC52-07NA27344. LLNL-ABS- 771103.

The depth-averaged shallow water wave equations are commonly used to model flows arising from natural hazards. The GeoClaw code, developed by D. George, R. J. LeVeque, M. J. Berger, K. Mandli and others is one example of a depth-averaged flow solver now widely used for modeling tsunamis, overland flooding, debris flows, storm surges and so on. Generally, depth averaged flow models show excellent large scale agreement with observations and can thus be reliably used to predict whether tsunamis will reach distant coast lines, and if, so can give vital information about arrival times. However, for other types of flows, dispersive effects missing from the SWE model can play an important role in determining localized effects such as whether waves will overtop seawalls, or whether a landslide entering a lake will trigger tsunami-like behavior on the opposite shore. Because of the importance of these dispersive effects, several depth averaged codes include dispersive corrections to the SWE. One set of equations commonly used to model these dispersive effects are the Serre-Green-Naghdi (SGN) equations. I will present my work to include dispersive correction terms into the GeoClaw extension of ForestClaw, a parallel adaptive library for Cartesian grid methods. One formulation of the SGN equations stabilizes higher order derivatives by treating them implicitly. As a result, a key component of an SGN solver is a variable coefficient Poisson solver. We will describe the SGN equations and provide an overview of their derivation, and then show preliminary results on uniform Cartesian meshes. Comparisons with the SGN solver in Basilisk (S. Popinet) and BoussClaw (J. Kim et al) will also be shown to verify our model. Preliminary results using the Hierarchical-Poincar\'e-Steklov (HPS) method developed by Gillman and Martinsson (2014) to solve the Poisson problem on adaptive meshes will also be shown.

The propagation of electromagnetic waves is modeled by time-dependent Maxwell's equations coupled with constitutive laws that describe the response of the media. In this work, we examine a nonlinear optical model that describes electromagnetic waves in linear Lorentz and nonlinear Kerr and Raman media. To design efficient, accurate, and stable computational methods, we apply high order discontinuous Galerkin discretizations and finite difference schemes in space. The challenge to achieve provable stability for fully-discrete methods lies in the temporal discretizations of the nonlinear terms. To overcome this, novel modification is proposed for the second-order leap-frog and implicit trapezoidal time integrators. The performance of the method is demonstrated via numerical examples.

Cut cells methods have been developed in recent years for computing flow around bodies with complicated geometries. They are an alternative to body fitted or unstructured grids, which may be harder to generate and more complex in the bulk of the flowfield. Cut cell methods cut the flow body out of a regular Cartesian grid resulting in so called cut cells. Cut cells can have irregular shape and may be very small. Therefore, they need special treatment. For the solution of hyperbolic conservation laws, probably the biggest issue caused by cut cells is the small cell problem: that explicit time stepping schemes are not stable on the arbitrarily small cut cells. In the context of finite volume schemes, several approaches have been suggested to successfully solve this problem. In the context of DG schemes however only very little work has been done so far. In this talk we present a new stabilization for overcoming the small cell problem in the context of piecewise linear DG schemes in one and two dimensions for the linear advection equation. Our stabilization is based on adding suitable penalty terms. In that sense it has a certain similarity to the ghost penalty stabilization [1] used for stabilizing the solution of elliptic problems on cut cells. To reflect the hyperbolic character of the problem considered however the stabilization terms that we suggest look fundamentally different and are closer in spirit to the h-box method [2]. Using the proposed space stabilization, one can use explicit time stepping even on cut cells. In one dimension, we show that the resulting method is monotone, TVD, and L 1 -stable for using piecewise constant polynomials in space and explicit Euler in time in the presence of cut cells. We also provide a TVD stability result for using piecewise linear polynomials (with limiter) and explicit second-order SSP Runge-Kutta time stepping. We conclude our talk with numerical results that confirm that both the stability considerations and the expected second-order accuracy of our scheme [3] transfer to two dimensions. References [1] E. Burman, Ghost penalty, C. R. Math. Acad. Sci. Paris, vol. 348, 2010. [2] C. Helzel, M. J. Berger, R. J. LeVeque, A high-resolution rotated grid method for conservation laws with embedded geometries, SIAM J. Sci. Comput., vol. 26, 2005. [3] C. Engwer, S. May, A. Nußing, F. Streitb\"urger. 'A stabilized DG cut-cell method for discretizing the linear transport equation. In preparation.

Stability and accuracy for a numerical method approximating an initial boundary value problem are inherently linked together. Stability means that perturbations have a bounded effect on the discrete solution, and is usually characterized by a precise estimate of norms. Such an estimate can be directly used to quantify the accuracy of the method. A very convenient and common way to investigate stability, and hence accuracy, is to use the energy method. If this approach fails one may instead attempt to get results by Laplace transforming in time and using normal mode analysis. Such analysis is usually more involved, but sharper results may follow. In this talk we will show two examples where, even though the energy method is applicable, it is rewarding to consider the problem in the Laplace domain. In the first case we get sharper accuracy results, and in the second case we get sharper temporal bounds.

The key to a successful recovery of seafloor characteristics from measured backscatter data generated from SONAR systems is finding a balance between the computational cost of forward modeling and the desired resolution. In the high frequency regime, many propagation models are often limited by costly simulations or unrealistic environmental assumptions, so a detailed recovery is not always possible. There is a growing need for sophisticated mathematical and computational tools that accommodate complicated scenarios, i.e., multiple seafloor layers or uncharted seafloor landscapes. To enable a rapid, remote, and accurate seafloor parameter recovery, we developed an approach a combination of localized forward modeling and deep learning. The idea is to partition environments in width into much smaller “template” domains, a few meters in width, where the sediment layer can be described using a few parameters. Machine learning is used to train a classifier using a reference library of simulations of Helmholtz equations on these domains. We investigate the potential of deep learning for classification of seafloor properties, such as sediment type, roughness, and thickness.

In this talk, I present a robust moving mesh finite difference method for the simulation of fourth order nonlinear PDEs describing elastic-electrostatic interactions in two dimensions. The parabolic Monge-Amp´ere methods from [1] are extended to solve a fourth order PDE with finite time singularity. A key feature in the implementation is the generation of a high order transformation between the computational and physical meshes that can accommodate the high order derivatives in the PDE [5]. The PDE derived from a plate contact problem develops finite time quenching singularities at discrete spatial location(s). The moving mesh method dynamically resolves these temporally forming singularities, while preserving the underlying length scales of the problem. I will show how the PMA resolves the singularities to high accuracy and gives strong evidence of self similarity near blow up. Next, I briefly discuss the prediction of the touchdown profile for given geometries using the skeleton theory from [3]. The skeleton set is numerically predicted for a variety of domains. The predictions of the skeleton method are verified using the variational moving mesh methods discretized in finite element space [2]. Accurately resolving singularities on general domains motivates recent work in extending the parabolic Monge-Amp´ere equation to problems on curved domains. Utilizing the optimal transport boundary formulation from [4, 6] creates a mapping between a fixed computational domain, on which all derivative calculations are made, to a curved physical domain. This creates an initial mesh mapping which can be paired to the PMA to adaptively resolve finescale features in the paired PDE. I give results of this method for a variety of examples including semi-linear blow-up, sharp interfaces and prescribed boundary motion on convex and select non-convex domains. [1] C.J. Budd and J.F. Williams. Moving mesh generation using the parabolic monge-ampere equation. SIAM Journal on Scientific Computing, 31(5):3438-3465, 2009. [2] K. DiPietro, R. Haynes, W. Huang, A. Lindsay, Y. Yu. Moving mesh simulation of contact sets in two dimensional models of elastic-electrostatic deflection problems. J. Compt. Phys., 375:761-782, 2018. [3] A.E. Lindsay and J. Lega, Multiple quenching solutions of a fourth order parabolic PDE with a singular nonlinearity modeling a MEMS capacitor. SIAM Journal On Applied Mathematics, 72(3):935-958, 2012. [4] J. Benamou, B. Froese, A. Oberman. Numerical solution of the optimal transportation problem using the Monge Ampere Equation. J. Compt. Phys., 260:107-126,2014. [5] K. DiPietro, A. E. Lindsay. Monge-Amp\'ere simulation of fourth order PDEs in two dimensions with application to elastic-electrostatic contact problems. J. Compt. Phys.. 349:328-350, 2017. [6] B. Froese. A numerical method for the elliptic Monge-Amp\'ere equation with transport boundary conditions. SIAM Journal on Scientific Computing, 34(3):A1432-A1459, 2012.

The extrapolation method known as Anderson acceleration has been used for decades to speed the convergence of nonlinear solvers in many applications. A mathematical justification of the improved convergence rate however has remained elusive. Here, we provide theory to establish the improved convergence rate. The key ideas of the analysis are relating the difference of consecutive iterates to residuals based on performing the inner-optimization in a Hilbert space setting, and explicitly defining the gain in the optimization stage to be the ratio of improvement over a step of the unaccelerated fixed point iteration. The main result we prove is this method of acceleration improves the convergence rate of a fixed point iteration to first order by a factor of the gain at each step as the method converges.

We've all looked at stability domains for ODE time-steppers. At the most basic level, these are found by studying how the time-stepper handles the ODE x' = sigma x where sigma is a complex number with negative real part. This leads to a stability domain that has a continuous boundary. The underlying analysis generalizes to systems of ODEs if the linearized system is diagonalizable. In this talk, I'll discuss an implicit-explicit time-stepping scheme for which the linearized system is not diagonalizable; standard stability theory doesn't apply. I'll demonstrate that an adaptive time-stepper can be used to explore the stability domain and I'll give an example of a system for which the stability domain can have a discontinuous boundary; a small change in a parameter can lead to a jump in the stability threshold of the time-step size. This is joint work with my former PhD student, Dave Yan.

A large toolbox of numerical schemes for the nonlinear Schrödinger equation has been established, based on different discretization techniques such as discretizing the variation-of-constants formula (e.g., exponential integrators) or splitting the full equation into a series of simpler subproblems (e.g., splitting methods). In many situations these classical schemes allow a precise and efficient approximation. This, however, drastically changes whenever “non-smooth’’ phenomena enter the scene such as for problems at low-regularity and high oscillations. Classical schemes fail to capture the oscillatory parts within the solution which leads to severe instabilities and loss of convergence. In this talk I present a new class of Fourier integrators for the nonlinear Schrödinger equation at low-regularity. The key idea in the construction of the new schemes is to tackle and hardwire the underlying structure of resonances into the numerical discretization.​ These terms are the cornerstones of theoretical analysis of the long time behaviour of differential equations and their numerical discretizations and offer the new schemes strong geometric structure at low regularity.

Low Mach number equation sets approximate the equations of motion of a compressible fluid by filtering out the sound waves, which allows the system to evolve on the advective rather than the acoustic time scale. Depending on the degree of approximation, low Mach number models retain some subset of possible compressible effects. In this talk I will discuss low Mach number models for reacting and stratified flows arising in combustion, astrophysics and atmospheric modeling.

Accurate evaluation of layer potentials near boundaries and interfaces are needed in many applications, including fluid-structure interaction problems. A classical method to approximate the solution everywhere in the domain consists of using the same quadrature rule (Nyström method) used to solve the underlying boundary integral equation. This method is problematic for evaluations close to boundaries and interfaces. For a fixed number, N, of quadrature points, this method incurs a non-uniform error with O(1) errors in a boundary layer of thickness O(1/N). We have developed new asymptotic methods to remove this error. To demonstrate this method, we consider the Laplace problem and show the methods extended to the Stokes equations.

We present a numerical method for computing the Stokeslet and stresslet integrals in Stokes flow, motivated by applications such as swimming of microorganisms, particle and drop motion, or biomembrane and red blood cell mechanics. Evaluating the integrals accurately for points near the boundary, e.g., when two interfaces are close together, is the most difficult case. The accurate solution is obtained by regularizing the kernels and adding analytically derived correction terms to eliminate the largest error. On the surface, high order regularizations are designed so that corrections are not required. To evaluate the resulting sums efficiently, we developed a kernel-independent treecode based on barycentric Lagrange interpolation.

We study 3rd order accurate finite volume methods for hyperbolic problems, which are inspired by the active flux method of Eymann and Roe. For one-dimensional linear problems we show that the unlimited active flux method can be interpreted as an ADER method. This interpretation motivates the construction of new third order accurate methods. The active flux method is based on a continuous reconstruction of the conserved quantities. This simplifies the construction of Cartesian grid cut cell methods.

5-day workshop participants are welcome to use BIRS facilities (BIRS Coffee Lounge, TCPL and Reading Room) until 3 pm on Friday, although participants are still required to checkout of the guest rooms by 12 noon.