This paper describes a new software tool that has been developed for the efficient solution of systems of linear and nonlinear partial differential equations (PDEs) of parabolic type. Specifically, the software is designed to provide optimal computational performance for multiscale problems, which require highly stable, implicit, time-stepping schemes combined with a parallel implementation of adaptivity in both space and time. By combining these implicit, adaptive discretizations with an optimally efficient nonlinear multigrid solver it is possible to obtain computational solutions to a very high resolution with relatively modest computational resources. The first half of the paper describes the numerical methods that lie behind the software, along with details of their implementation, whilst the second half of the paper illustrates the flexibility and robustness of the tool by applying it to two very different example problems. These represent models of a thin film flow of a spreading viscous droplet and a multi-phase-field model of tumour growth. We conclude with a discussion of the challenges of obtaining highly scalable parallel performance for a software tool that combines both local mesh adaptivity, requiring efficient dynamic load-balancing, and a multigrid solver, requiring careful implementation of coarse grid operations and inter-grid transfer operations in parallel. (C) 2016 Elsevier Ltd. All rights reserved.

Using a phase field model, which fully couples the thermal and solute concentration field, we present simulation results in three dimensions of the rapid dendritic solidification of a class of dilute alloys at the meso scale. The key results are the prediction of steady state tip velocity and radius at varying undercooling and thermal diffusivities. Less computationally demanding 2-dimensional results are directly compared with the corresponding 3-dimensional results, where significant quantitative differences emerge. The simulations provide quantitative predictions for the range of thermal and solutal diffusivities considered and show the effectiveness and potential of the computational techniques employed. These results thus provide benchmark 3-dimensional computations, allow direct comparison with underlying analytical theory, and pave the way for further quantitative results.

The first half of the paper provides an overview of a new engineering software tool that is designed for the efficient solution of problems that may be modeled as systems of linear and nonlinear partial differential equations (PDEs) of parabolic type. Our tool is built upon the PARAMESH library, [15], which provides hierarchical mesh adaptivity in parallel in two and three dimensions. Our discretizations are based upon cell-centred finite difference schemes in space and implicit multi-step methods in time (primarily the second order backward differential formula (BDF2)). This results in the need to solve a nonlinear algebraic system at each time step, and we have implemented an optimal nonlinear multigrid method based upon full approximation scheme (FAS). The second half of the presentation illustrates the application of this new software framework to a challenging application, namely a multi-phase-field model of tumour growth [18]. We show some typical simulations for growth of the model tumours, and these results demonstrate second-order convergence in both space and time. We conclude with a discussion of the challenges of obtaining highly scalable parallel performance for a software tool that combines both local mesh adaptivity (requiring efficient dynamic load-balancing) and a multigrid solver (requiring careful implementation of coarse grid operations and inter-grid transfer operations in parallel).

We employ adaptive mesh refinement, implicit time stepping, a nonlinear multigrid solver and parallel computation to solve a multi-scale, time dependent, three dimensional, nonlinear set of coupled partial differential equations for three scalar field variables. The mathematical model represents the non-isothermal solidification of a metal alloy into a melt substantially cooled below its freezing point at the microscale. Underlying physical molecular forces are captured at this scale by a specification of the energy field. The time rate of change of the temperature, alloy concentration and an order parameter to govern the state of the material (liquid or solid) are controlled by the diffusion parameters and variational derivatives of the energy functional. The physical problem is important to material scientists for the development of solid metal alloys and, hitherto, this fully coupled thermal problem has not been simulated in three dimensions, due to its computationally demanding nature. By bringing together state of the art numerical techniques this problem is now shown here to be tractable at appropriate resolution with relatively moderate computational resources.

This paper presents an automatic locally adaptive finite element solver for the fully-coupled EHL point contact problems. The proposed algorithm uses a posteriori error estimation in the stress in order to control adaptivity in both the elasticity and lubrication domains. The implementation is based on the fact that the solution of the linear elasticity equation exhibits large variations close to the fluid domain on which the Reynolds equation is solved. Thus the local refinement in such region not only improves the accuracy of the elastic deformation solution significantly but also yield an improved accuracy in the pressure profile due to increase in the spatial resolution of fluid domain. Thus, the improved traction boundary conditions lead to even better approximation of the elastic deformation. Hence, a simple and an effective way to develop an adaptive procedure for the fully-coupled EHL problem is to apply the local refinement to the linear elasticity mesh. The proposed algorithm also seeks to improve the quality of refined meshes to ensure the best overall accuracy. It is shown that the adaptive procedure effectively refines the elements in the region(s) showing the largest local error in their solution, and reduces the overall error with optimal computational cost for a variety of EHL cases. Specifically, the computational cost of proposed adaptive algorithm is shown to be linear with respect to problem size as the number of refinement levels grows.

The study of pathological cardiac conditions such as arrhythmias, a major cause of mortality in heart failure, is becoming increasingly informed by computational simulation, numerically modelling the governing equations. This can provide insight where experimental work is constrained by technical limitations and/or ethical issues. As the models become more realistic, the construction of efficient and accurate computational models becomes increasingly challenging. In particular, recent developments have started to couple the electrophysiology models with mechanical models in order to investigate the effect of tissue deformation on arrhythmogenesis, thus introducing an element of nonlinearity into the mathematical representation. This paper outlines a biophysically-detailed computational model of coupled electromechanical cardiac activity which uses the finite element method to approximate both electrical and mechanical systems on unstructured, deforming, meshes. An ILU preconditioner is applied to improve performance of the solver. This software is used to examine the role of electrophysiology, fibrosis and mechanical deformation on the stability of spiral wave dynamics in human ventricular tissue by applying it to models of both healthy and failing tissue. The latter was simulated by modifying (i) cellular electrophysiological properties, to generate an increased action potential duration and altered intracellular calcium handling, and (ii) tissue-level properties, to simulate the gap junction remodelling, fibrosis and increased tissue stiffness seen in heart failure. The resulting numerical experiments suggest that, for the chosen mathematical models of electrophysiology and mechanical response, introducing tissue level fibrosis can have a destabilising effect on the dynamics, while the net effect of the electrophysiological remodelling stabilises the system.

Tensor Network Theory (TNT) provides efficient and highly accurate algorithms for the simulation of strongly correlated quantum systems. The corresponding numerical algorithms enable approximate descriptions of many-body states and linear operators acting on them that do not grow exponentially with system size, in contrast to exact descriptions. Whilst TNT algorithms are efficient, they are numerically demanding and require high-performance optimised and parallelised implementations. A TNT library is currently being developed at the University of Oxford to allow users to have access to the complex algorithms needed to solve these problems from their own high level codes. This project is concerned with optimising and parallelising those parts of the TNT library that involve heavy computations, to enable important quantum effects in many-body systems to be studied. The initial objectives were: Obj. 1 Improving storage of the matrices by a more efficient usage of memory and parallelisation using OpenMP. Obj. 2 Developing more efficient and scalable calculations for the core functions of the TNT algorithms which are the most computationally demanding parts. This will mainly concern the contraction and SVD operations. Obj. 3 Improvements in the performance of the code will be achieved by incorporating symmetry information that will allow a decomposition of the tensors into sub-blocks so that they may then be assigned to individual MPI processes.

A fully implicit numerical method, based upon a combination of adaptively refined hierarchical meshes and geometric multigrid, is presented for the simulation of binary alloy solidification in three space dimensions. The computational techniques are presented for a particular mathematical model, based upon the phase-field approach, however their applicability is of greater generality than for the specific phase-field model used here. In particular, an implicit second order time discretization is combined with the use of second order spatial differences to yield a large nonlinear system of algebraic equations as each time step. It is demonstrated that these equations may be solved reliably and efficiently through the use of a nonlinear multigrid scheme for locally refined grids. In effect this paper presents an extension of earlier research in two space dimensions (J. Comput. Phys., 225 (2007), pp. 1271-1287) to fully three-dimensional problems. This extension is validated against earlier two-dimensional results and against some of the limited results available in three dimensions, obtained using an explicit scheme. The efficiency of the implicit approach and the multigrid solver are then demonstrated and some sample computational results for the simulation of the growth of dendrite structures are presented.

Surface roughness on membranes has been shown to increase flux, in part because surface area was increased. However, experimental studies of the relationship between flux and surface roughness have produced contradictory results in which flux did not always increase with increased surface roughness. Increases in flux that are greater than the increase in surface area also have been reported, with hydrogen flux through palladium and palladium–copper films as an example. A mathematical model was developed to examine two-dimensional diffusion through a roughened membrane when the surface is at local equilibrium with the feed. By comparing the results to a one-dimensional diffusion model, contributions of diffusion parallel to the plane of the membrane could be separated from that of a shorter diffusion path through the thin regions. Although lateral diffusion can be significant, more often the presence of thinner regions was the dominant factor for increased flux. The model calculations predict that membrane flux can increase by more or less than the increase in surface area depending on the geometry of the surface roughness. For the geometry of the surface structures in the hydrogen permeation experiments, model calculations indicate that flux increases larger than the increase in surface area could occur.

We review the application of advanced numerical techniques such as adaptive mesh refinement, implicit time stepping, multigrid solvers and massively parallel implementations as a route to obtaining solutions to the three-dimensional phase-field problem with a domain size and interface resolution previously possible only in two dimensions. Using such techniques it is shown that such models are tractable even as the interface width approaches the solute capillary length.

This paper presents the fast preconditioned iterative solution to large sparse linear systems arising from the application of Newton and quasi-Newton methods to fully coupled elastohydrodynamic lubrication line and point contact problems. The new blockwise preconditioner that is presented combines the use of multigrid for the linear elasticity block and a separate approximation to precondition the Reynolds block. Two variants of the solver are considered, based upon the use of algebraic and geometric multigrid respectively. Numerical results are presented in order to validate the discretization and solution method and then to contrast the performance and efficiency of the proposed solution strategies compared to the use of a state-of-the-art sparse direct solver. These results demonstrate that, unlike the sparse direct solver, the preconditioned iterative approach is able to perform at computational and memory costs that both grow linearly with the number of unknowns.

Elastohydrodynamic lubrication modelling plays an important role in engineering design and analysis, since a number of important mechanical components operate under elastohydrodynamic lubrication conditions. In this article, methods are presented for solving both line and point contact cases using multiphysics software. The advantages, and the overheads, of using such an approach over developing highly specialised, bespoke software are highlighted. In order to calculate the deformation of the contacts three different methods are developed and their relative performance is assessed. The advantage of using a nested solution strategy has also been examined. The flexibility of the multiphysics software approach is highlighted in results involving a complex transient case modelling an involute gear.

The fully coupled approach for the solution of elastohydrodynamic lubrication point contact problems requires the numerical solution of the elasticity problem on a large 3D domain. A sufficiently fine computational mesh is required to obtain the elastic deformation solution to the necessary accuracy, and hence, a sufficiently accurate elastohydrodynamic lubrication point contact solution. This article discusses the accuracy of an elastohydrodynamic lubrication point contact solution over different finite element meshes. We present a family of efficient 3D meshes which can be used to calculate accurate elastohydrodynamic lubrication point contact solutions at a minimal computational cost. This article also illustrates the fact that the unstructured hierarchical meshes can lead to a poor quality elastohydrodynamic lubrication solution unless an appropriate post-processing, smoothing technique is applied.

This paper describes an adaptive implementation of a high order Discontinuous Galerkin (DG) method for the solution of Elastohydrodynamic Lubrication (EHL) point contact problems. These problems arise when modelling the thin lubricating film between contacts which are under sufficiently high pressure that the elastic deformation of the contacting elements cannot be neglected. The governing equations are highly non-linear and include a second order partial differential equation that is derived via the thin-film approximation. Furthermore, the problem features a free boundary, which models where cavitation occurs, and this is automatically captured as part of the solution process. The need for spatial adaptivity stems from the highly variable length scales that are present in typical solutions. Results are presented which demonstrate both the effectiveness and the limitations of the proposed adaptive algorithm.

Ionic conductivity in nanocomposite electrolytes is examined through the use of a numerical model. A rigorous description of the space charge layer and its impact on conductivity are developed for a composite system consisting of insulating spheres dispersed within an ion conducting material. Model simulations are performed to understand how the effective conductivity, which can exceed the conductivity of the bulk material, depends on the volume fraction, size, configuration, and particle size distribution of the nanoparticles in the bulk material. Several deliberately chosen regular particle configurations are used to establish the lower and upper bounds for conductivity enhancement. A simple cubic array of particles is demonstrated to provide a reasonable estimate for the behavior expected from a random distribution of particles. Finally, conductivity modulation is shown to be significant when the particle radius is comparable to or smaller than the thickness of the space charge layer. (C) 2012 Elsevier Ltd. All rights reserved.

Potential conductivity enhancement due to formation of space charge layers in nanoionic composites was computed numerically for several structures of non-conducting nanoparticles in a bulk ionic conductor. Optimum loading fractions were extracted from simulation results and found to depend strongly on the thickness of the space charge layer relative to the size of the particles. This behavior agreed well with an approximate analytical expression derived herein. In certain cases, significant conductivity enhancement was predicted, even for nanoparticle loadings as small as 0.5 volume %. The model was also applied to a space charge layer depletion scenario, and found to be in good agreement with results from a recent experimental study. (C) 2011 Elsevier Ltd. All rights reserved.

The study of cardiac arrhythmias is a major focus of computational biology, and undertaking biophysically detailed simulations is computationally demanding. An efficient coupled electromechanical solver to model cardiac tissue has been developed. This provides features to model fibre direction, and utilises computationally efficient techniques to reduce the simulation times. In this paper the break up of human re-entrant arrhythmias has been simulated. The results suggest that tissue deformation is a contributory factor in the break up of stable re-entrant spiral waves.

The use of an adjoint technique for goal‐based error estimation described by Hartit et al. ( Int. J. Numer. Meth. Fluids 2005; 47 :1069–1074) is extended to the numerical solution of free boundary problems that arise in elastohydrodynamic lubrication (EHL). EHL systems are highly nonlinear and consist of a thin‐film approximation of the flow of a non‐Newtonian lubricant which separates two bodies that are forced together by an applied load, coupled with a linear elastic model for the deformation of the bodies. A finite difference discretization of the line contact flow problem is presented, along with the numerical evaluation of an exact solution for the elastic deformation, and a moving grid representation of the free boundary that models cavitation at the outflow in this one‐dimensional case. The application of a goal‐based error estimate for this problem is then described. This estimate relies on the solution of an adjoint problem; its effectiveness is demonstrated for the physically important goal of the total friction through the contact. Finally, the application of this error estimate to drive local mesh refinement is demonstrated. Copyright © 2010 John Wiley & Sons, Ltd.