Search

In this paper we present a family of high order cut finite element methods with bound preserving properties for hyperbolic conservation laws in one space dimension. The methods are based on the discontinuous Galerkin framework and use a regular background mesh, where interior boundaries are allowed to cut through the mesh arbitrarily. Our methods include ghost penalty stabilization to handle small cut elements and a new reconstruction of the approximation on macro-elements, which are local patches consisting of cut and un-cut neighboring elements that are connected by stabilization. We show that the reconstructed solution retains conservation and order of convergence. Our lowest-order scheme results in a piecewise constant solution that satisfies a maximum principle for scalar hyperbolic conservation laws. When the lowest order scheme is applied to the Euler equations, the scheme is positivity preserving in the sense that positivity of pressure and density are retained. For the high-order schemes, suitable bound preserving limiters are applied to the reconstructed solution on macro-elements. In the scalar case, a maximum principle limiter is applied, which ensures that the limited approximation satisfies the maximum principle. Correspondingly, we use a positivity preserving limiter for the Euler equations and show that our scheme is positivity preserving. In the presence of shocks, additional limiting is needed to avoid oscillations, hence we apply a standard TVB limiter to the reconstructed solution. The time step restrictions are of the same order as for the corresponding discontinuous Galerkin methods on the background mesh. Numerical computations illustrate accuracy, bound preservation, and shock capturing capabilities of the proposed schemes., Comment: 32

We develop a hybrid spatial discretization for the wave equation in second order form, based on high-order accurate finite difference methods and discontinuous Galerkin methods. The hybridization combines computational efficiency of finite difference methods on Cartesian grids and geometrical flexibility of discontinuous Galerkin methods on unstructured meshes. The two spatial discretizations are coupled by a penalty technique at the interface such that the overall semidiscretization satisfies a discrete energy estimate to ensure stability. In addition, optimal convergence is obtained in the sense that when combining a fourth order finite difference method with a discontinuous Galerkin method using third order local polynomials, the overall convergence rate is fourth order. Furthermore, we use a novel approach to derive an error estimate for the semidiscretization by combining the energy method and the normal mode analysis for a corresponding one dimensional model problem. The stability and accuracy analysis are verified in numerical experiments.

This paper presents a stable finite element approximation for the acoustic wave equation on second-order form, with perfectly matched layers (PML) at the boundaries. Energy estimates are derived for varying PML damping for both the discrete and the continuous case. Moreover, a priori error estimates are derived for constant PML damping. Most of the analysis is performed in Laplace space. Numerical experiments in physical space validate the theoretical results., Comment: 30 pages

We develop a family of cut finite element methods of different orders based on the discontinuous Galerkin framework, for hyperbolic conservation laws with stationary interfaces in both one and two space dimensions, and for moving interfaces in one space dimension. Interface conditions are imposed weakly and so that both conservation and stability are ensured. A CutFEM with discontinuous elements in space is developed and coupled to standard explicit time-stepping schemes for linear advection problems and the acoustic wave problem with stationary interfaces. In the case of moving interfaces, we propose a space-time CutFEM based on discontinuous elements both in space and time for linear advection problems. We show that the proposed CutFEM are conservative and energy stable. For the stationary interface case an a priori error estimate is proven. Numerical computations in both one and two space dimensions support the analysis, and in addition demonstrate that the proposed methods have the expected accuracy.

It is well-known that reliable and efficient domain truncation is crucial to accurate numerical solution of most wave propagation problems. The perfectly matched layer (PML) is a method which, when stable, can provide a domain truncation scheme which is convergent with increasing layer width/damping. The difficulties in using the PML are primarily associated with stability, which can be present at the continuous level or be triggered by numerical approximations. The mathematical and numerical analysis of the PML for hyperbolic wave propagation problems has been an area of active research. It is now possible to construct stable and high order accurate numerical wave solvers by augmenting wave equations with the PML and approximating the equations using summation-by-parts finite difference methods, continuous and discontinuous Galerkin finite element methods. In this review we summarise the progress made, from mathematical, numerical and practical perspectives, point out some open problems and set the stage for future work. We also present numerical experiments of model problems corroborating the theoretical analysis, and numerical simulations of real-world wave propagation demonstrating impact. Stable and parallel implementations of the PML in the high performance computing software packages WaveQLab3D and ExaHyPE allow to sufficiently limit the computational domain of seismological problems with only a few grid points/elements around the computational boundaries where the PML is active, thus saving as much as $96\%$ of the required computational resources for a three space dimensional seismological benchmark problem.

In this paper, we develop a family of high order cut discontinuous Galerkin (DG) methods for hyperbolic conservation laws in one space dimension. The ghost penalty stabilization is used to stabilize the scheme for small cut elements. The analysis shows that our proposed methods have similar stability and accuracy properties as the standard DG methods on a regular mesh. We also prove that the cut DG method with piecewise constants in space is total variation diminishing (TVD). We use the strong stability preserving Runge-Kutta method for time discretization and the time step is independent of the size of cut element. Numerical examples demonstrate that the cut DG methods are high order accurate for smooth problems and perform well for discontinuous problems., Comment: 29 pages, 10 figures

We develop a new finite difference method for the wave equation in second order form. The finite difference operators satisfy a summation-by-parts (SBP) property. With boundary conditions and material interface conditions imposed weakly by the simultaneous-approximation-term (SAT) method, we derive energy estimates for the semi-discretization. In addition, error estimates are derived by the normal mode analysis. The proposed method is termed as energy-based because of its similarity with the energy-based discontinuous Galerkin method. When imposing the Dirichlet boundary condition and material interface conditions, the traditional SBP-SAT discretization uses a penalty term with a mesh-dependent parameter, which is not needed in our method. Furthermore, numerical dissipation can be added to the discretization through the boundary and interface conditions. We present numerical experiments that verify convergence and robustness of the proposed method.

We present a stable discontinuous Galerkin (DG) method with a perfectly matched layer (PML) for three and two space dimensional linear elastodynamics, in velocity-stress formulation, subject to well-posed linear boundary conditions. First, we consider the elastodynamics equation, in a cuboidal domain, and derive an unsplit PML truncating the domain using complex coordinate stretching. Leveraging the hyperbolic structure of the underlying system, we construct continuous energy estimates, in the time domain for the elastic wave equation, and in the Laplace space for a sequence of PML model problems, with variations in one, two and three space dimensions, respectively. They correspond to PMLs normal to boundary faces, along edges and in corners. Second, we develop a DG numerical method for the linear elastodynamics equation using physically motivated numerical flux and penalty parameters, which are compatible with all well-posed, internal and external, boundary conditions. When the PML damping vanishes, by construction, our choice of penalty parameters yield an upwind scheme and a discrete energy estimate analogous to the continuous energy estimate. Third, to ensure numerical stability of the discretization when PML damping is present, it is necessary to extend the numerical DG fluxes, and the numerical inter-element and boundary procedures, to the PML auxiliary differential equations. This is crucial for deriving discrete energy estimates analogous to the continuous energy estimates. By combining the DG spatial approximation with the high order ADER time stepping scheme and the accuracy of the PML we obtain an arbitrarily accurate wave propagation solver in the time domain. Numerical experiments are presented in two and three space dimensions corroborating the theoretical results.

The conventional no-slip boundary condition leads to a non-integrable stress singularity at a contact line. This is a main challenge in numerical simulations of two-phase flows with moving contact lines. We derive a two-dimensional hydrodynamic model for the velocity field at a contact point moving with constant velocity over a curved wall. The model is a perturbation of the classical Huh and Scriven hydrodynamic solution [11], which is only valid for flow over a flat wall. The purpose of the hydrodynamic model is to investigate the macroscopic behavior of the fluids close to a contact point. We also present an idea for how the hydrodynamic solution could be used to prescribe macroscopic Dirichlet boundary conditions for the velocity in the vicinity of a moving contact point. Simulations demonstrate that the velocity field based on the non-singular boundary conditions is capable of accurately advecting the contact point.

Atomistic-continuum multiscale modelling is becoming an increasingly popular tool for simulating the behaviour of materials due to its computational efficiency and reliable accuracy. In the case of ferromagnetic materials, the atomistic approach handles the dynamics of spin magnetic moments of individual atoms, while the continuum approximations operate with volume-averaged quantities, such as magnetisation. One of the challenges for multiscale models in relation to physics of ferromagnets is the existence of the long-range dipole-dipole interactions between spins. The aim of the present paper is to demonstrate a way of including these interactions into existing atomistic-continuum coupling methods based on the partitioned-domain and the upscaling strategies. This is achieved by modelling the demagnetising field exclusively at the continuum level and coupling it to both scales. Such an approach relies on the atomistic expression for the magnetisation field converging to the continuum expression when the interatomic spacing approaches zero, which is demonstrated in this paper.

Motivated by the emerging applications of liquid-infused surfaces (LIS), we study the drag reduction and robustness of transverse flows over two-dimensional microcavities partially filled with an oily lubricant. Using separate simulations at different scales, characteristic contact line velocities at the fluid-solid intersection are first extracted from nano-scale phase field simulations and then applied to micron-scale two-phase flows, thus introducing a multiscale numerical framework to model the interface displacement and deformation within the cavities. As we explore the various effects of the lubricant-to-outer-fluid viscosity ratio $\tilde{\mu}_2/\tilde{\mu}_1$, the capillary number Ca, the static contact angle $\theta_s$, and the filling fraction of the cavity $\delta$, we find that the effective slip is most sensitive to the parameter $\delta$. The effects of $\tilde{\mu}_2/\tilde{\mu}_1$ and $\theta_s$ are generally intertwined, but weakened if $\delta < 1$. Moreover, for an initial filling fraction $\delta =0.94$, our results show that the effective slip is nearly independent of the capillary number, when it is small. Further increasing Ca to about $0.01 \tilde{\mu}_1/\tilde{\mu}_2$, we identify a possible failure mode, associated with lubricants draining from the LIS, for $\tilde{\mu}_2/\tilde{\mu}_1 \lesssim 0.1$. Very viscous lubricants (\eg $\tilde{\mu}_2/\tilde{\mu}_1 >1$), on the other hand, are immune to such failure due to their generally larger contact line velocity., Comment: Phys. Rev. Fluids (2018)

A high order cut finite element method is formulated for solving the elastic wave equation. Both a single domain problem and an interface problem are treated. The boundary or interface are allowed to cut through the background mesh. To avoid problems with small cuts, stabilizing terms are added to the bilinear forms corresponding to the mass and stiffness matrix. The stabilizing terms penalize jumps in normal derivatives over the faces of the elements cut by the boundary/interface. This ensures a stable discretization independently of how the boundary/interface cuts the mesh. Nitsche's method is used to enforce boundary and interface conditions, resulting in symmetric bilinear forms. As a result of the symmetry, an energy estimate can be made and optimal order a priori error estimates are derived for the single domain problem. Finally, numerical experiments in two dimensions are presented that verify the order of accuracy and stability with respect to small cuts.

We develop a provably energy stable discontinuous Galerkin spectral element method (DGSEM) approximation of the perfectly matched layer (PML) for the three and two space dimensional (3D and 2D) linear acoustic wave equations, in first order form, subject to well-posed linear boundary conditions. First, using the well-known complex coordinate stretching, we derive an efficient un-split modal PML for the 3D acoustic wave equation. Second, we prove asymptotic stability of the continuous PML by deriving energy estimates in the Laplace space, for the 3D PML in a heterogeneous acoustic medium, assuming piece-wise constant PML damping. Third, we develop a DGSEM for the wave equation using physically motivated numerical flux, with penalty weights, which are compatible with all well-posed, internal and external, boundary conditions. When the PML damping vanishes, by construction, our choice of penalty parameters yield an upwind scheme and a discrete energy estimate analogous to the continuous energy estimate. Fourth, to ensure numerical stability when PML damping is present, it is necessary to systematically extend the numerical numerical fluxes, and the inter-element and boundary procedures, to the PML auxiliary differential equations. This is critical for deriving discrete energy estimates analogous to the continuous energy estimates. Finally, we propose a procedure to compute PML damping coefficients such that the PML error converges to zero, at the optimal convergence rate of the underlying numerical method. Numerical experiments are presented in 2D and 3D corroborating the theoretical results.

The conventional no-slip boundary condition leads to a non-integrable stress singularity at a moving contact line. This makes numerical simulations challenging, especially when capillary effects are essential for the dynamics of the flow. This paper presents a new boundary methodology, suitable for numerical simulation of flow of two immiscible and incompressible fluids in the presence of moving contact points. The methodology is based on combining a relation between the apparent contact angle and the contact point velocity with the similarity solution for Stokes flow at a planar interface. The relation between angle and velocity can be determined by theoretical arguments, or from simulations using a more detailed model. The approach here uses the phase field model in a micro domain, with physically relevant parameters for molecular diffusion and interface thickness. The methodology is used to formulate a new boundary condition for the velocity. Numerical results illustrate the usefulness.

The present study concerns the numerical homogenization of second order hyperbolic equations in non-divergence form, where the model problem includes a rapidly oscillating coefficient function. These small scales influence the large scale behavior, hence their effects should be accurately modelled in a numerical simulation. A direct numerical simulation is prohibitively expensive since a minimum of two points per wavelength are needed to resolve the small scales. A multiscale method, under the equation free methodology, is proposed to approximate the coarse scale behaviour of the exact solution at a cost independent of the small scales in the problem. We prove convergence rates for the upscaled quantities in one as well as in multi-dimensional periodic settings. Moreover, numerical results in one and two dimensions are provided to support the theory.

In this work, we discuss and compare three methods for the numerical approximation of constant- and variable-coefficient diffusion equations in both single and composite domains with possible discontinuity in the solution/flux at interfaces, considering (i) the Cut Finite Element Method; (ii) the Difference Potentials Method; and (iii) the summation-by-parts Finite Difference Method. First we give a brief introduction for each of the three methods. Next, we propose benchmark problems, and consider numerical tests-with respect to accuracy and convergence-for linear parabolic problems on a single domain, and continue with similar tests for linear parabolic problems on a composite domain (with the interface defined either explicitly or implicitly). Lastly, a comparative discussion of the methods and numerical results will be given., Comment: 45 pages, 12 figures, in revision for Journal of Scientific Computing

In this article, a few problems related to multiscale modelling of magnetic materials at finite temperatures and possible ways of solving these problems are discussed. The discussion is mainly centred around two established multiscale concepts: the partitioned domain and the upscaling-based methodologies. The major challenge for both multiscale methods is to capture the correct value of magnetisation length accurately, which is affected by a random temperature-dependent force. Moreover, general limitations of these multiscale techniques in application to spin systems are discussed., Comment: 30 pages

When using a finite difference method to solve an initial--boundary--value problem, the truncation error is often of lower order at a few grid points near boundaries than in the interior. Normal mode analysis is a powerful tool to analyze the effect of the large truncation error near boundaries on the overall convergence rate, and has been used in many previous literatures for different equations. However, existing work only concerns problems in one space dimension. In this paper, we extend the analysis to problems in two space dimensions. The two dimensional analysis is based on a diagonalization procedure that decomposes a two dimensional problem to many one dimensional problems of the same type. We present a general framework of analyzing convergence for such one dimensional problems, and explain how to obtain the result for the corresponding two dimensional problem. In particular, we consider two kinds of truncation errors in two space dimensions: the truncation error along an entire boundary, and the truncation error localized at a few grid points close to a corner of the computational domain. The accuracy analysis is in a general framework, here applied to the second order wave equation. Numerical experiments corroborate our accuracy analysis.

The scalar wave equation is solved using higher order immersed finite elements. We demonstrate that higher order convergence can be obtained. Small cuts with the background mesh are stabilized by adding penalty terms to the weak formulation. This ensures that the condition numbers of the mass and stiffness matrix are independent of how the boundary cuts the mesh. The penalties consist of jumps in higher order derivatives integrated over the interior faces of the elements cut by the boundary. The dependence on the polynomial degree of the condition number of the stabilized mass matrix is estimated. We conclude that the condition number grows extremely fast when increasing the polynomial degree of the finite element space. The time step restriction of the resulting system is investigated numerically and is concluded not to be worse than for a standard (non-immersed) finite element method.

We consider a multiscale strategy addressing the disparate scales in the Landau-Lifschitz equations in micro-magnetism. At the microscopic scale, the dynamics of magnetic moments are driven by a high frequency field. On the macroscopic scale we are interested in simulating the dynamics of the magnetisation without fully resolving the microscopic scales. The method follows the framework of heterogeneous multiscale methods and it has two main ingredients: a micro- and a macroscale model. The microscopic model is assumed to be known exactly whereas the macro model is incomplete as it lacks effective quantities. The two models use different temporal and spatial scales and effective parameter values for the macro model are computed on the fly, allowing for improved efficiency over traditional one-scale schemes. For the analysis, we consider a single spin under a high frequency field and show that effective quantities can be obtained accurately with step-sizes much larger than the size of the microscopic scales required to resolve the microscopic features. Numerical results both for a single magnetic particle as well as a chain of interacting magnetic particles are given to validate the theory., Comment: 23 pages

We compare two high order finite-difference methods that solve the elastic wave equation in two dimensional domains with curved boundaries and material discontinuities. Two numerical experiments are designed with focus on wave boundary interaction, the response of a pressure wave impinging on a circular cavity and the wave field generated by an explosive impact on the wall an underground tunnel. Qualitative comparisons of the results are made where similarities and differences are pointed out., Comment: 12 pages, 5 figures

A main challenge in numerical simulations of moving contact line problems is that the adherence, or no-slip boundary condition leads to a non-integrable stress singularity at the contact line. In this report we perform the first steps in developing the macroscopic part of an accurate multiscale model for a moving contact line problem in two space dimensions. We assume that a micro model has been used to determine a relation between the contact angle and the contact line velocity. An intermediate region is introduced where an analytical expression for the velocity exists. This expression is used to implement boundary conditions for the moving contact line at a macroscopic scale, along a fictitious boundary located a small distance away from the physical boundary. Model problems where the shape of the interface is constant thought the simulation are introduced. For these problems, experiments show that the errors in the resulting contact line velocities converge with the grid size $h$ at a rate of convergence $p\approx 2$.

We use high order finite difference methods to solve the wave equation in the second order form. The spatial discretization is performed by finite difference operators satisfying a summation-by-parts property. The focus of this work is on the numerical treatment of non-conforming grid interfaces. The interface conditions are imposed weakly by the simultaneous approximation term technique in combination with interface operators, which move the discrete solutions between the grids on the interface. In particular, we consider interpolation operators and projection operators. A norm-compatibility condition, which leads to stability for first order hyperbolic systems, does not suffice for second order wave equations. An extra constraint on the interface operators must be satisfied to derive an energy estimate for stability. We carry out eigenvalue analyses to investigate the additional constraint and how it is related to stability, and find that the projection operators have better stability properties than the interpolation operators. In addition, a truncation error analysis is performed to study the convergence property of the numerical schemes. In the numerical experiments, the stability and accuracy properties of the numerical schemes are further explored, and the practical usefulness of non-conforming grid interfaces is presented and discussed in two efficiency studies., Comment: 27 pages, 15 figures

In this paper, we consider finite difference approximations of the second order wave equation. We use finite difference operators satisfying the summation-by-parts property to discretize the equation in space. Boundary conditions and grid interface conditions are imposed by the simultaneous-approximation-term technique. Typically, the truncation error is larger at the grid points near a boundary or grid interface than that in the interior. Normal mode analysis can be used to analyze how the large truncation error affects the convergence rate of the underlying stable numerical scheme. If the semi-discretized equation satisfies a determinant condition, two orders are gained from the large truncation error. However, many interesting second order equations do not satisfy the determinant condition. We then carefully analyze the solution of the boundary system to derive a sharp estimate for the error in the solution and acquire the gain in convergence rate. The result shows that stability does not automatically yield a gain of two orders in convergence rate. The accuracy analysis is verified by numerical experiments., Comment: In version 2, we have added a new section on the convergence analysis of the Neumann problem, and have improved formulations in many places

A recent study shows that the classical theory concerning accuracy and points per wavelength is not valid for surface waves in almost incompressible elastic materials. The grid size must instead be proportional to $(\frac{\mu}{\lambda})^{(1/p)}$ to achieve a certain accuracy. Here $p$ is the order of accuracy the scheme and $\mu$ and $\lambda$ are the Lame parameters. This accuracy requirement becomes very restrictive close to the incompressible limit where $\frac{\mu}{\lambda} \ll 1$, especially for low order methods. We present results concerning how to choose the number of grid points for 4th, 6th and 8th order summation-by-parts finite difference schemes. The result is applied to Lambs problem in an almost incompressible material., Comment: 14 pages, 5 figures

We prove nonlinear stability for finite amplitude perturbations of plane Couette flow. A bound of the solution of the resolvent equation in the unstable complex half-plane is used to estimate the solution of the full nonlinear problem.The result is a lower bound, including Reynolds number dependence,of the threshold amplitude below which all perturbations are stable. Our result is an improvement of the corresponding bound derived in \cite{KLH:1}., Comment: arxiv version is already official