PDE Solver Research Articles

A reduced order accelerator for time-dependent segregated neutronic solvers

The deterministic solution of the neutron transport problem entails the coupled solution of several partial differential equations, one for each energy group, direction and/or spherical harmonic. Several techniques have been devised for accelerating the solution of this set of equations, both for time dependent and eigenvalue calculations. This paper describes an acceleration technique based on reduced order models and applicable to the segregated solution of time dependent solutions. In this work the technique is applied to the simple case of multi-group diffusion and tested on two cases of practical interest. It shows performances that are comparable to some commonly employed acceleration techniques. Some potential advantages have been observed for transients with significant flux deformations. In addition, possibly interesting features of the proposed technique are: a relatively easy implementation in general PDE solvers and numerical libraries; its potential applicability to any kind of problem requiring the iterative solution of a system of equations; a flexible implementation with a wide margin for possible modifications.

On the approximation of functions with line singularities by ridgelets

In Grohs and Obermeier (2015), the authors discussed the existence of numerically feasible solvers for advection equations that run in optimal computational complexity. In this paper, we complete the last remaining requirement to achieve this goal — by showing that ridgelets, on which the solver is based, approximate functions with line singularities (which may appear as solutions to the advection equation) with the best possible approximation rate (up to arbitrarily small δ>0) if the functions decay sufficiently fast, and with a slight penalty if not.Structurally, the proof resembles Candès (2001), where a similar result was proved for a different ridgelet construction, which is however not well-suited for use in a PDE solver (and in particular, not suitable for the CDD-schemes (Cohen et al., 2001) we are interested in). Due to the differences between the two ridgelet constructions, we have to deal with quite a different set of issues, but are also able to relax the (support) conditions on the function being approximated. Finally, the proof employs a new convolution-type estimate that could be of independent interest due to its sharpness.

Heat Transfer in Fe3O4-H2O Nanofluid Contained in a Triangular Cavity Under a Sloping Magnetic Field

Numerical simulation is performed to explore the convective heat transfer characteristics of Fe3O4-H2O nanofluid contained in a right-angle triangular cavity considering three types of thermal boundary conditions at the bottom wall. No heat is allowed to escape through the insulated vertical wall, whereas the inclined wall is kept colder than the bottom one. A sloping magnetic field whose strength is unvarying acts upon the cavity. The physical model is converted to the mathematical form through coupled highly nonlinear partial differential equations. These equations are then transformed into the non-dimensional form with the help of a group of transformations of variables. A very robust pde solver COMSOL Multiphysics that uses the finite element method (FEM) of Galerkin type is applied to carry out the numerical calculation. Heat transfer escalation through middling Nusselt number at the lowermost cavity wall is explored for diverse model parameters and thermal circumstances. The outcomes lead us to conclude that a higher degree of heat transfer is accomplished by reducing the dimension of nanoparticles and aggregating the buoyancy force through the Rayleigh number. It is highest when there is a magnetic field leaning angle of 900 and the lowermost wall is heated homogenously.

Mitigating the influence of the boundary on PDE-based covariance operators

Gaussian random fields over infinite-dimensional Hilbert spaces require the definition of appropriate covariance operators. The use of elliptic PDE operators to construct covariance operators allows to build on fast PDE solvers for manipulations with the resulting covariance and precision operators. However, PDE operators require a choice of boundary conditions, and this choice can have a strong and usually undesired influence on the Gaussian random field. We propose two techniques that allow to ameliorate these boundary effects for large-scale problems. The first approach combines the elliptic PDE operator with a Robin boundary condition, where a varying Robin coefficient is computed from an optimization problem. The second approach normalizes the pointwise variance by rescaling the covariance operator. These approaches can be used individually or can be combined. We study properties of these approaches, and discuss their computational complexity. The performance of our approaches is studied for random fields defined over simple and complex two- and three-dimensional domains.

Evolution of permeability during the process of shale gas extraction

The sensitivity of shale permeability to effective stress has been reported widely in the literature but the relation between them is not always consistent with the principle of effective stress. This knowledge gap defines our specific research goals: (1) to resolve the inconsistency between the evolution of shale permeability and the effective stress principle, and (2) to incorporate the evolution of shale permeability into the evaluation of shale gas production rate. In this study, we define shale blocks bounded by artificially introduced hydraulic fractures. At block scale the system comprises two components represented by inorganic and organic components. Unlike previous studies, we assume that the organic system is embedded within the inorganic system and both systems are dual porosity media. We define the shale permeability as a function of strains representing the stress transfer between fractures and matrix. These permeability-strain relations (shale permeability models) are then used to couple the governing equations of shale deformation and gas flow. The composite suite of equations are then solved by a commercial PDE solver. We use the coupled model to generate typical evolution profiles of shale permeability under stress-controlled conditions for the cases of both gas injection and gas desorption. These profiles are consistent with experimental observations as reported widely in the literature. We also use the coupled model to examine how the shale permeability evolves during the extraction of shale gas. Model results show that the evolutions of both shale permeability and gas production rate at field scale are primarily controlled by the contrasts in both the transport and deformation properties of the reservoir and the interfacial dynamics of mass and effective stress transfer between matrix and fractures.

Incremental proper orthogonal decomposition for PDE simulation data

We propose an incremental algorithm to compute the proper orthogonal decomposition (POD) of simulation data for a partial differential equation. Specifically, we modify an incremental matrix SVD algorithm of Brand to accommodate data arising from Galerkin-type simulation methods for time dependent PDEs. The algorithm is applicable to data generated by many numerical methods for PDEs, including finite element and discontinuous Galerkin methods. The algorithm initializes and efficiently updates the dominant POD eigenvalues and modes during the time stepping in a PDE solver without storing the simulation data. We prove that the algorithm without truncation updates the POD exactly. We demonstrate the effectiveness of the algorithm using finite element computations for a 1D Burgers’ equation and a 2D Navier–Stokes problem.

Numerical simulation of natural convection heat transfer in a trapezoidal enclosure filled with nanoparticles

In this paper, we have investigated unsteady natural convection flow and heat transfer inside a trapezoidal enclosure filled with nine different types of nanofluids having various shapes of the nanoparticle following Tiwari and Das mathematical model. The left and right walls of the enclosure are kept at different temperatures, while the top and bottom walls of the cavity are thermally insulated. The Galerkin weighted residual based finite element method has been employed to solve the governing partial differential equations after converting them into a nondimensional form. The simulation is carried out through the pde solver COMSOL Multiphysics with Matlab interface. Comparison with the previously published result is made for a special case and an excellent agreement is found. The effects of various model parameters such as the Rayleigh number, the aspect ratio, the volume fraction and the shape factor of the nanoparticles on streamlines and isotherms have been displayed graphically and discussed. The heat transfer augmentation for various combinations of pertinent parameters has also been presented in light of the average Nusselt number on the left heated wall.

Non-Standard Finite Difference Schemes for Solving Singular Lane-Emden Equation

In this paper we construct Non-Standard finite difference schemes (NSFD) for numerical solution of nonlinear Lane-Emden type equations which are nonlinear ordinary dierential equations on semi-infinite domain. They are categorized as singular initial value problems. This equation describes a variety of phenomena in theoretical physics and astrophysics. The presented schemes are obtained by using the Non-Standard finite difference method. The use of NSFD method and its approximations play an important role for the formation of stable numerical methods. The main advantage of the schemes is that the algorithm is very simple and very easy to implement. Thus, this method may be applied as a simple and accurate solver for ODEs and PDEs and it can also be utilized as an accurate algorithm to solve linear and nonlinear equations arising in physics and other fields of applied mathematics. Illustrative examples have been discussed to demonstrate validity and applicability of the technique and the results have been compared with the exact solutions.

Hybrid PDE solver for data-driven problems and modern branching

The numerical solution of large-scale PDEs, such as those occurring in data-driven applications, unavoidably require powerful parallel computers and tailored parallel algorithms to make the best possible use of them. In fact, considerations about the parallelization and scalability of realistic problems are often critical enough to warrant acknowledgement in the modelling phase. The purpose of this paper is to spread awareness of the Probabilistic Domain Decomposition (PDD) method, a fresh approach to the parallelization of PDEs with excellent scalability properties. The idea exploits the stochastic representation of the PDE and its approximation via Monte Carlo in combination with deterministic high-performance PDE solvers. We describe the ingredients of PDD and its applicability in the scope of data science. In particular, we highlight recent advances in stochastic representations for non-linear PDEs using branching diffusions, which have significantly broadened the scope of PDD. We envision this work as a dictionary giving large-scale PDE practitioners references on the very latest algorithms and techniques of a non-standard, yet highly parallelizable, methodology at the interface of deterministic and probabilistic numerical methods. We close this work with an invitation to the fully non-linear case and open research questions.

Free Convective Boundary Layer Flow of Radiating and Reacting MHD Fluid Past a Cosinusoidally Fluctuating Heated Plate

The work is focused on unsteady magneto-hydrodynamics free convective flow past a moving vertical porous plate in the presence of thermal radiation and first order chemical reaction with viscous dissipation. The temperature of the plate is assumed spanwise cosinusoidally fluctuating with time by taking into account of heat generation. The second order perturbation technique is employed to investigate the non-linear partial differential equations governing the fluid flow which are non-dimensionalized by introducing the similarity transformations. For the numerical computation, flex PDE Solver is used to study the effects of magnetic intensity, radiation, porous permeability, Eckert number, Schmidt number and heat generation/absorption parameters on velocity and temperature profiles. Also, the skin friction coefficients, the rate of heat transfer and the rate of mass transfer at the surface of the plate are computed numerically. The results show that within the boundary layer, the velocity and temperature are found to decrease with the increasing values of Prandtl number and radiation parameter however the trend is reverse with respect to porous permeability and heat generation/ absorption parameters.

Adaptive multilevel trust-region methods for time-dependent PDE-constrained optimization

We present a class of adaptive multilevel trust-region methods for the efficient solution of optimization problems governed by time-dependent nonlinear partial differential equations with control constraints. The algorithm is based on the ideas of the adaptive multilevel inexact SQP-method from \[26], \[27]. It is in particular well suited for problems with time-dependent PDE constraints. Instead of the quasi-normal step in a classical SQP method which results in solving the linearized PDE sufficiently well, in this algorithm a (nonlinear) solver is applied to the current discretization of the PDE. Moreover, different discretizations and solvers for the PDE and the adjoint PDE may be applied. The resulting inexactness of the reduced gradient in the current discretization is controlled within the algorithm. Thus, highly efficient PDE solvers can be coupled with the proposed optimization framework. The algorithm starts with a coarse discretization of the underlying optimization problem and provides during the optimization process implementable criteria for an adaptive refinement strategy of the current discretization based on error estimators. We prove global convergence to a stationary point of the infinite-dimensional problem. Moreover, we illustrate how the adaptive refinement strategy of the algorithm can be implemented by using a posteriori error estimators for the state and the adjoint equation. Numerical results for a semilinear parabolic PDE-constrained problem with pointwise control constraints are presented.

A resilient domain decomposition polynomial chaos solver for uncertain elliptic PDEs

A resilient method is developed for the solution of uncertain elliptic PDEs on extreme scale platforms. The method is based on a hybrid domain decomposition, polynomial chaos (PC) framework that is designed to address soft faults. Specifically, parallel and independent solves of multiple deterministic local problems are used to define PC representations of local Dirichlet boundary-to-boundary maps that are used to reconstruct the global solution. A LAD-lasso type regression is developed for this purpose. The performance of the resulting algorithm is tested on an elliptic equation with an uncertain diffusivity field. Different test cases are considered in order to analyze the impacts of correlation structure of the uncertain diffusivity field, the stochastic resolution, as well as the probability of soft faults. In particular, the computations demonstrate that, provided sufficiently many samples are generated, the method effectively overcomes the occurrence of soft faults.

A LAGRANGIAN GAUSS-NEWTON-KRYLOV SOLVER FOR MASS- AND INTENSITY-PRESERVING DIFFEOMORPHIC IMAGE REGISTRATION.

We present an efficient solver for diffeomorphic image registration problems in the framework of Large Deformations Diffeomorphic Metric Mappings (LDDMM). We use an optimal control formulation, in which the velocity field of a hyperbolic PDE needs to be found such that the distance between the final state of the system (the transformed/transported template image) and the observation (the reference image) is minimized. Our solver supports both stationary and non-stationary (i.e., transient or time-dependent) velocity fields. As transformation models, we consider both the transport equation (assuming intensities are preserved during the deformation) and the continuity equation (assuming mass-preservation). We consider the reduced form of the optimal control problem and solve the resulting unconstrained optimization problem using a discretize-then-optimize approach. A key contribution is the elimination of the PDE constraint using a Lagrangian hyperbolic PDE solver. Lagrangian methods rely on the concept of characteristic curves. We approximate these curves using a fourth-order Runge-Kutta method. We also present an efficient algorithm for computing the derivatives of the final state of the system with respect to the velocity field. This allows us to use fast Gauss-Newton based methods. We present quickly converging iterative linear solvers using spectral preconditioners that render the overall optimization efficient and scalable. Our method is embedded into the image registration framework FAIR and, thus, supports the most commonly used similarity measures and regularization functionals. We demonstrate the potential of our new approach using several synthetic and real world test problems with up to 14.7 million degrees of freedom.

A direct solver for elliptic PDEs in three dimensions based on hierarchical merging of Poincaré–Steklov operators

A numerical method for variable coefficient elliptic PDEs on three dimensional domains is described. The method is designed for problems with smooth solutions, and is based on a multidomain spectral collocation discretization scheme. The resulting system of linear equations can very efficiently be solved using a nested dissection style direct (as opposed to iterative) solver. This makes the scheme particularly well suited to solving problems for which iterative solvers struggle; in particular for problems with oscillatory solutions. A principal feature of the scheme is that once the solution operator has been constructed, the actual solve is extremely fast. An upper bound on the asymptotic cost of the build stage of O(N4/3) is proved (for the case where the PDE is held fixed as N increases). The solve stage has close to linear complexity. The scheme requires a relatively large amount of storage per degree of freedom, but since it is a high order scheme, a small number of degrees of freedom is sufficient to achieve high accuracy. The method is presented for the case where there is no body load present, but it can with little difficulty be generalized to the non-homogeneous case. Numerical experiments demonstrate that the scheme is capable of solving Helmholtz-type equations on a domain of size 20×20×20 wavelengths to three correct digits on a modest personal workstation, with N≈2⋅106.

A fully coupled multiscale shale deformation-gas transport model for the evaluation of shale gas extraction

Horizontal drilling and hydraulic fracturing are two enabling technologies to create a shale gas reservoir. For the created reservoir scale, we define shale blocks between hydraulic fractures as matrixes. In the matrix scale, flow processes are defined in the components of inorganic minerals and kerogens, respectively. Under this framework, a set of partial differential equations are derived to define various flow and deformation processes: (1) mechanical equilibrium equation that defines the shale deformation; (2) gas flow in the kerogen system of matrix; (3) gas flow in the inorganic system of matrix; and (4) gas flow in the hydraulic fracture system. For each of gas flow systems, a permeability or diffusivity model is derived to define its evolution. All systems are fully coupled through these permeability models and mass exchanges between different systems. The fully couple PDE system was solved by using COMSOL, a popular PDE solver. The model was verified against gas production data from the Marcellus Shale and the Barnett Shale, respectively. The verified model was applied to investigate the impact of adsorption parameters, flow regimes (Knudsen number), initial permeability of the inorganic matrix, and the effective stress variations on the gas production. Model results show that the Langmuir parameters affect both the cumulative gas production and the gas extraction processes; that the impact of flow regimes is closely related to the initial permeability of the inorganic matrix; and that the impact of effective stress variations on the permeability of hydraulic fractures is more significant than that on the matrix system.

Positive Filtered P$_N$ Moment Closures for Linear Kinetic Equations

We propose a positive-preserving moment closure for linear kinetic transport equations based on a filtered spherical harmonic (FP$_N$) expansion in the angular variable. The recently proposed FP$_N$ moment equations are known to suffer from the occurrence of (unphysical) negative particle concentrations. The origin of this problem is that the FP$_N$ approximation is not always positive at the kinetic level; the new FP$^+_N$ closure is developed to address this issue. A new spherical harmonic expansion is computed via the solution of an optimization problem, with constraints that enforce positivity, but only on a finite set of preselected points. Combined with an appropriate PDE solver for the moment equations, this ensures positivity of the particle concentration at each step in the time integration. Under an additional, mild regularity assumption, we prove that FP$^+_N$ has the same consistency as FP$_N$; that is, the FP$^+_N$ approximation converges to a given target function in $L^2$ at the same rate a...

Deterministic Mean-Field Ensemble Kalman Filtering

The proof of convergence of the standard ensemble Kalman filter (EnKF) from Legland etal. (2011) is extended to non-Gaussian state space models. A density-based deterministic approximation of the mean-field limit EnKF (DMFEnKF) is proposed, consisting of a PDE solver and a quadrature rule. Given a certain minimal order of convergence $\kappa$ between the two, this extends to the deterministic filter approximation, which is therefore asymptotically superior to standard EnKF when the dimension $d<2\kappa$. The fidelity of approximation of the true distribution is also established using an extension of total variation metric to random measures. This is limited by a Gaussian bias term arising from non-linearity/non-Gaussianity of the model, which exists for both DMFEnKF and standard EnKF. Numerical results support and extend the theory.

Petascale solvers for anisotropic PDEs in atmospheric modelling on GPU clusters

Multi-GPU parallelisation of scalable CG and multigrid solvers for anisotropic PDEs.Efficient matrix-free CUDA implementation minimises global memory access on GPUs.Excellent weak scaling on up to 16,384 nVidia K20X cards (44 mio. cores, Titan, OLCF).PDEs with 5x1011 unknowns can be solved in 1?s to an accuracy of 10 - 5 .Achieved performance of 0.78PFLOPs and memory bandwidth utilisation of more than 40%. Memory bound applications such as solvers for large sparse systems of equations remain a challenge for GPUs. Fast solvers should be based on numerically efficient algorithms and implemented such that global memory access is minimised. To solve systems with trillions ( O ( 10 12 ) ) unknowns the code has to make efficient use of several million individual processor cores on large GPU clusters.We describe the multi-GPU implementation of two algorithmically optimal iterative solvers for anisotropic PDEs which are encountered in (semi-) implicit time stepping procedures in atmospheric modelling. In this application the condition number is large but independent of the grid resolution and both methods are asymptotically optimal, albeit with different absolute performance. In particular, an important constant in the discretisation is the CFL number; only the multigrid solver is robust to changes in this constant. We parallelise the solvers and adapt them to the specific features of GPU architectures, paying particular attention to efficient global memory access. We achieve a performance of up to 0.78 PFLOPs when solving an equation with 0.55 ? 1012unknowns on 16384 GPUs; this corresponds to about 3% of the theoretical peak performance of the machine and we use more than 40% of the peak memory bandwidth with a Conjugate Gradient (CG) solver. Although the other solver, a geometric multigrid algorithm, has a slightly worse performance in terms of FLOPs per second, overall it is faster as it needs less iterations to converge; the multigrid algorithm can solve a linear PDE with half a trillion unknowns in about one second.

Toward Predictive Multiscale Modeling of Vascular Tumor Growth

New directions in medical and biomedical sciences have gradually emerged over recent years that will change the way diseases are diagnosed and treated and are leading to the redirection of medicine toward patient-specific treatments. We refer to these new approaches for studying biomedical systems as predictive medicine, a new version of medical science that involves the use of advanced computer models of biomedical phenomena, high-performance computing, new experimental methods for model data calibration, modern imaging technologies, cutting-edge numerical algorithms for treating large stochastic systems, modern methods for model selection, calibration, validation, verification, and uncertainty quantification, and new approaches for drug design and delivery, all based on predictive models. The methodologies are designed to study events at multiple scales, from genetic data, to sub-cellular signaling mechanisms, to cell interactions, to tissue physics and chemistry, to organs in living human subjects. The present document surveys work on the development and implementation of predictive models of vascular tumor growth, covering aspects of what might be called modeling-and-experimentally based computational oncology. The work described is that of a multi-institutional team, centered at ICES with strong participation by members at M. D. Anderson Cancer Center and University of Texas at San Antonio. This exposition covers topics on signaling models, cell and cell-interaction models, tissue models based on multi-species mixture theories, models of angiogenesis, and beginning work of drug effects. A number of new parallel computer codes for implementing finite-element methods, multi-level Markov Chain Monte Carlo sampling methods, data classification methods, stochastic PDE solvers, statistical inverse algorithms for model calibration and validation, models of events at different spatial and temporal scales is presented. Importantly, new methods for model selection in the presence of uncertainties fundamental to predictive medical science, are described which are based on the notion of Bayesian model plausibilities. Also, as part of this general approach, new codes for determining the sensitivity of model outputs to variations in model parameters are described that provide a basis for assessing the importance of model parameters and controlling and reducing the number of relevant model parameters. Model specific data is to be accessible through careful and model-specific platforms in the Tumor Engineering Laboratory. We describe parallel computer platforms on which large-scale calculations are run as well as specific time-marching algorithms needed to treat stiff systems encountered in some phase-field mixture models. We also cover new non-invasive imaging and data classification methods that provide in vivo data for model validation. The study concludes with a brief discussion of future work and open challenges.

Nonlinear multigrid methods for second order differential operators with nonlinear diffusion coefficient

Nonlinear multigrid methods such as the Full Approximation Scheme (FAS) and Newton-multigrid (Newton-MG) are well established as fast solvers for nonlinear PDEs of elliptic and parabolic type. In this paper we consider Newton-MG and FAS iterations applied to second order differential operators with nonlinear diffusion coefficient. Under mild assumptions arising in practical applications, an approximation (shown to be sharp) of the execution time of the algorithms is derived, which demonstrates that Newton-MG can be expected to be a faster iteration than a standard FAS iteration for a finite element discretisation. Results are provided for elliptic and parabolic problems, demonstrating a faster execution time as well as greater stability of the Newton-MG iteration. Results are explained using current theory for the convergence of multigrid methods, giving a qualitative insight into how the nonlinear multigrid methods can be expected to perform in practice.