Multiscale Modeling & Simulation

We propose two algorithms for simulating continuous time Markov chains in the presence of metastability. We show that the algorithms correctly estimate, under the ergodicity assumption, stationary averages of the process. Both algorithms, based on the idea of the parallel replica method, use parallel computing in order to explore metastable sets more efficiently. The algorithms require no assumptions on the Markov chains beyond ergodicity and the presence of identifiable metastability. In particular, there is no assumption on reversibility. For simpler illustration of the algorithms, we assume that a synchronous architecture is used throughout the paper. We present error analyses, as well as numerical simulations on multiscale stochastic reaction network models in order to demonstrate consistency of the method and its efficiency.

We introduce a model of multiagent dynamics for self-organized motion; individuals travel at a constant speed while trying to adopt the averaged body attitude of their neighbors. The body attitudes are represented through unitary quaternions. We prove the correspondence with the model presented in [P. Degond, A. Frouvelle, and S. Merino-Aceituno, Math. Models Methods Appl. Sci., 27 (2017), pp. 1005--1049], where the body attitudes are represented by rotation matrices. Differently from this previous work, the individual-based model introduced here is based on nematic (rather than polar) alignment. From the individual-based model, the kinetic and macroscopic equations are derived. The benefit of this approach, in contrast to that of the previous one, is twofold: first, it allows for a better understanding of the macroscopic equations obtained and, second, these equations are prone to numerical studies, which is key for applications.

Multi-scale/level computer models are efficient for simulating solutions of partial differential equations on nontrivial geometries using a large number of modes. This is because fine-scale structures in the quantities of interest cannot, in practice, be resolved using a single-level finite element method (FEM) applied directly to the full model, or standard spectral basis functions are not appropriate for the nontrivial geometry and related constraints. The multi-scale/level FEM (MFEM) is based on first applying a FEM for a reduced model problem, obtained by restriction to either (i) several local subdomains or (ii) a dominant linear differential operator in the model. The resulting structure-aware precomputed solutions are then used as basis functions to simulate the full model in the MFEM framework. If these are eigenfunctions of a reduced model eigenvalue problem, such basis functions provide spectrally accurate approximations to the solution of the full model and the resulting method is called the spectral MFEM (SMFEM). In this article we develop, analyze, and implement time-space fully discrete implicit SMFEM (ISMFEM) algorithms for efficiently simulating a Ginzburg--Landau (GL) system modeling superconductivity on a class of superconducting surfaces $S$. The spatial linear differential operator in the time-dependent nonlinear GL model is the Schrödinger operator $( i \nabla + {A}_0)^2$ on $S$, with magnetic vector potential ${A}_0$. For our ISMFEM, the spectrally accurate basis functions are eigenfunctions of a reduced stationary linear model governed by $(i \nabla + {A}_0)^2$. In a recent work we developed, analyzed, and implemented an algorithm to compute these functions using an analytic-numeric hybrid method, in conjunction with a high-order FEM. Using these precomputed solutions as basis functions, the spectrally accurate algorithms developed and analyzed for our new GL computer models on nontrivial geometries are validated using both nonlinear and linearized implicit time discretizations.

Sea ice hosts a rich ecosystem of flora and fauna, from microscale to macroscale. Algae living in its porous brine microstructure, such as the diatom Melosira arctica, secrete gelatinous exopolymeric substances (EPS) which are thought to protect these communities from their cold and highly saline environment. Recent experimental work has shown significant changes in the structure and properties of young sea ice with entrained Melosira EPS, such as increased brine volume fraction, salt retention, pore tortuosity, and decreased fluid permeability. In particular, we find that the cross-sectional areas of the brine inclusions are described by a bimodal--lognormal distribution, which generalizes the classic lognormal distribution of Perovich and Gow. We propose a model for the effective fluid permeability of young, EPS-laden sea ice, consisting of a random network of pipes with cross-sectional areas chosen from this bimodal distribution. We consider an equilibrium model posed on a square lattice, incorporating only the most basic features of the geometry and connectivity of the brine microstructure, and find good agreement between our model and the observed drop in fluid permeability. Our model formulation suggests future directions for experimental work, focused on measuring the inclusion size distribution and fluid permeability of sea ice with entrained EPS as functions of brine volume fraction. The drop in fluid permeability observed in experimental work and predicted by the model is significant, and should be taken into account, for example, in physical or ecological process models involving fluid or nutrient transport.

To characterize the Brownian motion in a bounded domain $\Omega$, it is well known that the boundary conditions of the classical diffusion equation just rely on the given information of the solution along the boundary of a domain; in contrast, for the Lévy flights or tempered Lévy flights in a bounded domain, the boundary conditions involve the information of a solution in the complementary set of $\Omega$, i.e., $\Bbb{R}^n\backslash \Omega$, with the potential reason that paths of the corresponding stochastic process are discontinuous. Guided by probability intuitions and the stochastic perspectives of anomalous diffusion, we show the reasonable ways, ensuring the clear physical meaning and well-posedness of the partial differential equations (PDEs), of specifying “boundary” conditions for space fractional PDEs modeling the anomalous diffusion. Some properties of the operators are discussed, and the well-posednesses of the PDEs with generalized boundary conditions are proved.

We study the convergence of numerical solutions of the compressible Navier--Stokes system to its incompressible limit. The numerical solution is obtained by a combined finite element--finite volume method based on the linear Crouzeix--Raviart finite element for the velocity and piecewise constant approximation for the density. The convective terms are approximated using upwinding. The distance between a numerical solution of the compressible problem and the strong solution of the incompressible Navier--Stokes equations is measured by means of a relative energy functional. For barotropic pressure exponent $\gamma \geq 3/2$ and for well-prepared initial data we obtain uniform convergence of order ${\cal O}(\sqrt{\Delta t}, h^a, \varepsilon)$, $a = \min \{ \frac{2 \gamma - 3 }{ \gamma}, 1\}$. Extensive numerical simulations confirm that the numerical solution of the compressible problem converges to the solution of the incompressible Navier--Stokes equations as the discretization parameters $\Delta t$, $h$ and the Mach number $\varepsilon$ tend to zero.

Rhythmic processes are found at all biological and ecological scales, and are fundamental to the efficient functioning of living systems in changing environments. The biochemical mechanisms underpinning these rhythms are therefore of importance, especially in the context of anthropogenic challenges such as pollution or changes in climate and land use. Here we develop and test a new method for clustering rhythmic biological data with a focus on circadian oscillations. The method combines locally stationary wavelet time series modelling with functional principal components analysis and thus extracts the time-scale patterns arising in a range of rhythmic data. We demonstrate the advantages of our methodology over alternative approaches, by means of a simulation study and real data applications, using both a published circadian dataset and a newly generated one. The new dataset records plant response to various levels of stress induced by a soil pollutant, a biological system where existing methods which assume stationarity are shown to be inappropriate. Our method successfully clusters the circadian data in an interesting way, thereby facilitating wider ranging analyses of the response of biological rhythms to environmental changes.

In this paper, we consider the harmonic extension problem, which is widely used in many applications of machine learning. We formulate the harmonic extension as solving a Laplace--Beltrami equation with Dirichlet boundary condition. We use the point integral method (PIM) proposed in [Z. Li, Z. Shi, and J. Sun, Commun. Comput. Phys., 22 (2017), pp. 228--258; Z. Shi and J. Sun, Res. Math. Sci., to appear; Z. Li and Z. Shi, Multiscale Model. Simul., 14 (2016), pp. 874--905] to solve the Laplace--Beltrami equation. The basic idea of the PIM method is to approximate the Laplace equation using an integral equation, which is easy to discretize from points. Based on the integral equation, we found that the traditional graph Laplacian method (GLM) fails to approximate the harmonic functions near the boundary. One important application of the harmonic extension in machine learning is semisupervised learning. We run a popular semisupervised learning algorithm by Zhu, Ghahramani, and Lafferty [Machine Learning, Proceedings of the Twentieth International Conference (ICML 2003), 2003, Washington, DC, 2003, ARAI Press, Menlo Park, CA pp. 912--919] over a couple of well-known datasets and compare the performance of the aforementioned approaches. Our experiments show the PIM performs the best. We also apply PIM to an image recovery problem and show it outperforms GLM. Finally, on the model problem of the Laplace--Beltrami equation with Dirichlet boundary, we prove the convergence of the point integral method.

In this paper, we present a new, optimization-based method to exhibit cyclic behavior in nonreversible stochastic processes. While our method is general, it is strongly motivated by discrete simulations of ordinary differential equations representing nonreversible biological processes, in particular, molecular simulations. Here, the discrete time steps of the simulation are often very small compared to the time scale of interest, i.e., of the whole process. In this setting, the detection of a global cyclic behavior of the process becomes difficult because transitions between individual states may appear almost reversible on the small time scale of the simulation. We address this difficulty using a mixed-integer programming model that allows us to compute a cycle of clusters with maximum net flow, i.e., large forward and small backward probability. For a synthetic genetic regulatory network consisting of a ring oscillator with three genes, we show that this approach can detect cycles that have a productivity one magnitude larger than classical spectral analysis methods. Our method applies to general nonequilibrium steady state systems such as catalytic reactions, for which the objective value computes the effectiveness of the catalyst.

By linking atomistic and mesoscopic scales, we formally show how a local steric effect can hinder crystal growth and lead to a buildup of adsorbed atoms (adatoms) on a supersaturated, (1+1)-dimensional surface. Starting from a many-adatom master equation of a kinetic restricted solid-on-solid (KRSOS) model with external material deposition, we heuristically extract a coarse-grained, mesoscale description that defines the motion of a line defect (i.e., a step) in terms of statistical averages over KRSOS microstates. Near thermodynamic equilibrium, we use error estimates to show that this mesoscale picture can deviate from the standard Burton--Cabrera--Frank (BCF) step flow model in which the adatom flux at step edges is linear in the adatom supersaturation. This deviation is caused by the accumulation of adatoms near the step, which block one another from being incorporated into the crystal lattice. In the mesoscale picture, this deviation manifests as a significant contribution from many-adatom microstates to the corresponding statistical averages. We carry out kinetic Monte Carlo simulations to numerically demonstrate how certain parameters control the aforementioned deviation. From these results, we discuss empirical corrections to the BCF model that amount to a nonlinear relation for the adatom flux at the step. We also discuss how this work could be used to understand the kinetic interplay between accumulation of adatoms and step motion in recent experiments of ice surfaces.

In this paper we develop magnetic induction conforming multiscale formulations for magnetoquasistatic problems involving periodic materials. The formulations are derived using the periodic homogenization theory and applied within a heterogeneous multiscale approach. Therefore the fine-scale problem is replaced by a macroscale problem defined on a coarse mesh that covers the entire domain and many mesoscale problems defined on finely-meshed small areas around some points of interest of the macroscale mesh (e.g., numerical quadrature points). The exchange of information between these macro and meso problems is thoroughly explained in this paper. For the sake of validation, we consider a two-dimensional geometry of an idealized periodic soft magnetic composite. (An erratum is attached.)

In the paper, we present a strategy for accelerating posterior inference for unknown inputs in time fractional diffusion models. In many inference problems, the posterior may be concentrated in a small portion of the entire prior support. It will be much more efficient if we build and simulate a surrogate only over the significant region of the posterior. To this end, we construct a coarse model using the generalized multiscale finite element method (GMsFEM) and solve a least-squares problem for the coarse model with a regularizing Levenberg--Marquart algorithm. An intermediate distribution is built based on the approximate sampling distribution. For Bayesian inference, we use GMsFEM and the least-squares stochastic collocation method to obtain a reduced coarse model based on the intermediate distribution. To increase the sampling speed of Markov chain Monte Carlo, the DREAM$_{ZS}$ algorithm is used to explore the surrogate posterior density, which is based on the surrogate likelihood and the intermediate distribution. The proposed method with lower generalized polynomial chaos order gives the approximate posterior as accurately as the surrogate model directly based on the original prior. A few numerical examples for time fractional diffusion equations are carried out to demonstrate the performance of the proposed method with applications of the Bayesian inversion.

In this paper we use layer potentials and asymptotic analysis techniques to analyze the heat generation due to nanoparticles when illuminated at their plasmonic resonance. We consider arbitrary-shaped particles and the cases of both a single and multiple particles. We clarify the strong dependency of the heat generation on the geometry of the particles as it depends on the eigenvalues of the associated Neumann--Poincaré operator. For close-to-touching nanoparticles, we show that the temperature field deviates significantly from the one generated by two single nanoparticles. The results of this paper formally explain experimental results reported in the nanomedical literature. They open a door for solving the challenging problems of detecting plasmonic nanoparticles in biological media and monitoring temperature elevation in tissue generated by nanoparticle heating.

In this paper, we suggest a new heterogeneous multiscale method (HMM) for the Helmholtz equation with high contrast. The method is constructed for a setting as in Bouchitte and Felbacq [C. R. Math. Acad. Sci. Paris, 339 (2004), pp. 377-382], where the high contrast in the parameter leads to unusual effective parameters in the homogenized equation. We revisit existing homogenization approaches for this special setting and analyze the stability of the two-scale solution with respect to the wavenumber and the data. This includes a new stability result for solutions to the Helmholtz equation with discontinuous diffusion matrix. The HMM is defined as direct discretization of the two-scale limit equation. With this approach we are able to show quasi-optimality and an a priori error estimate under a resolution condition that inherits its dependence on the wavenumber from the stability constant for the analytical problem. Numerical experiments confirm our theoretical convergence results and examine the resolution condition. Moreover, the numerical simulation gives a good insight and explanation of the physical phenomenon of frequency band gaps.

This paper introduces a novel framework for model adaptivity in the context of heterogeneous multiscale problems. The framework is based on the idea to interpret model adaptivity as a minimization problem of local error indicators that are derived in the general context of the dual weighted residual (DWR) method. Based on the optimization approach a postprocessing strategy is formulated that lifts the requirement of strict a priori knowledge about applicability and quality of effective models. This allows for the systematic, “goal-oriented” tuning of effective models} with respect to a quantity of interest. The framework is tested numerically on elliptic diffusion problems with different types of heterogeneous, random coefficients, as well as an advection-diffusion problem with a strong microscopic, random advection field.

To make the investigation of electronic structure of incommensurate heterostructures computationally tractable, effective alternatives to Bloch theory must be developed. In [Multiscale Model. Simul., 15(2017), pp. 476--499] we developed and analyzed a real space scheme that exploits spatial ergodicity and near-sightedness. In the present work, we present an analogous scheme formulated in momentum space, which we prove has significant computational advantages in specific incommensurate systems of physical interest, e.g., bilayers of a specified class of materials with small rotation angles. We use our theoretical analysis to obtain estimates for improved rates of convergence with respect to total CPU time for our momentum space method that are confirmed in computational experiments.

Digital Gabor filters are indispensable tools of local time-frequency analysis in signal processing. With strong orientation selectivity, discrete (tight) Gabor frames generated by two-dimensional Gabor filters also see their wide applications in image processing and volume data processing. However, owing to the lack of multiscale structures, discrete Gabor frames are less effective than multiresolution analysis (MRA) based wavelet (tight) frames when being used for modeling data composed of local structures with varying sizes. Recently, it was shown that digital Gabor filters do generate MRA-based wavelet tight frames via the unitary extension principle. However, the corresponding window function has to be a constant window, which has poor joint time-frequency resolution. In this paper, we showed that digital Gabor filters with smooth window function can generate MRA-based wavelet biframes. The MRA-based wavelet biframes generated by digital Gabor filters have both the advantages of Gabor systems on local time-frequency analysis and the advantages of wavelet systems on multiscale analysis.

In this work, we investigate the low-rank approximation of elliptic problems in heterogeneous media by means of Kolmogrov $n$-width and asymptotic expansion. This class of problems arises in many practical applications involving high-contrast media, and their efficient numerical approximation often relies crucially on certain low-rank structure of the solutions. We provide conditions on the permeability coefficient $\kappa$ that ensure a favorable low-rank approximation. These conditions are expressed in terms of the distribution of the inclusions in the coefficient $\kappa$, e.g., the values, locations, and sizes of the heterogeneous regions. Further, we provide a new asymptotic analysis for high-contrast elliptic problems based on the perfect conductivity problem and layer potential techniques, which allows deriving new estimates on the spectral gap for such high-contrast problems. These results provide theoretical underpinnings for several multiscale model reduction algorithms.

In many applications the resolution of small-scale heterogeneities remains a significant hurdle to robust and reliable predictive simulations. In particular, while material variability at the mesoscale plays a fundamental role in processes such as material failure, the resolution required to capture mechanisms at this scale is often computationally intractable. Multiscale methods aim to overcome this difficulty through judicious choice of a subscale problem and a robust manner of passing information between scales. One promising approach is the multiscale finite element method, which increases the fidelity of macroscale simulations by solving lower-scale problems that produce enriched multiscale basis functions. In this study, we present the first work toward application of the multiscale finite element method to the nonlocal peridynamic theory of solid mechanics. This is achieved within the context of a discontinuous Galerkin framework that facilitates the description of material discontinuities and does not assume the existence of spatial derivatives. Analysis of the resulting nonlocal multiscale finite element method is achieved using the ambulant Galerkin method, developed here with sufficient generality to allow for application to multiscale finite element methods for both local and nonlocal models that satisfy minimal assumptions. We conclude with preliminary results on a mixed-locality multiscale finite element method in which a nonlocal model is applied at the fine scale and a local model at the coarse scale.

In this paper we investigate the possibility of reducing the complexity of a system composed of a large number of interacting agents, whose dynamics feature a symmetry breaking. We consider first order stochastic differential equations describing the behavior of the system at the particle (i.e., Lagrangian) level and we get its continuous (i.e., Eulerian) counterpart via a kinetic description. However, the resulting continuous model alone fails to describe adequately the evolution of the system, due to the loss of granularity which prevents it from reproducing the symmetry breaking of the particle system. By suitably coupling the two models we are able to reduce considerably the necessary number of particles while still keeping the symmetry breaking and some of its large-scale statistical properties. We describe such a multiscale technique in the context of opinion dynamics, where the symmetry breaking is induced by the results of some opinion polls reported by the media.