ABSTRACT. Coupled 1D-3D flow models are used for a variety of applications, such as modelling fluid flow through vascularized tissue, modelling the flow of water and nutrients through soil embedded with a root system, or modelling the interaction between a well and reservoir. Veins, arteries, roots and wells all have in common that their radius is negligible compared to their length and the size of the domain as a whole. For this reason, we idealize them as being 1D geometries. The 1D structures are then endowed with a 1D flow equation, and coupled to the 3D flow equation by the use of a line source.

The main challenge associated with coupled 1D-3D flow models is that the line source makes the solution singular. This complicates both the analysis and approximation of the solution. In this talk, we show that the solution admits a splitting into two parts: (i) a term that explicitly captures the singularity and (ii) some smooth remainder. Via this splitting, we can then subtract the singularity. This yields a reformulated model in which all variables are smooth. The solution can then be approximated using any standard numerical method. We conclude by showing numerical experiments relevant to biomedical applications.

ABSTRACT. Fluid flow and solute transport in fractured porous media is the backbone of many applications including groundwater flow, underground energy harvesting, earthquake prediction, and biomedical engineering. The traditional continuous Galerkin (CG) method is not suitable for the transport equation due to lack of mass conservation. The discontinuous Galerkin (DG) method mitigates this problem; however, its computational cost is considerably more than the CG method. In this study, a robust and efficient discretization method based on the incomplete interior penalty enriched Galerkin (EG) method is proposed. This method requires fewer degrees of freedom than those of the DG method, while it achieves the same accuracy. The flow and transport models of rock matrix and fractures domains are investigated in the mixed-dimensional setting. The results of combinations of function spaces, for example, (i) CG × CG, (ii) CG × EG, and (iii) CG × DG spaces are compared. The results illustrate the superiority of the EG and DG methods in solving the flow and transport equations in fractured porous media. Furthermore, the computational burden of the EG method is two times cheaper than that of the DG method.

ABSTRACT. Non-Newtonian fluids are found in all aspects of life, from bodily fluids to engine oil. Thus advances in modeling and simulation of non-Newtonian fluids can have broad impact. Models of non-Newtonian fluids have been studied extensively for many years, but only recently have there been mathematical advances that allow models for them to be understood more completely. This understanding now allows development of numerical solution methods with a new level of reliability. Here we consider some Oldroyd models and their relation to the grade-two model. In this talk, we present a computational study of a solution method using tools from the FEniCS Project to facilitate code generation and to support correctness of the implementation.

ABSTRACT. We present an Inverse Problem PYthon library (hIPPYlib) for solving large-scale deterministic and Bayesian inverse problems governed by partial differential equations (PDEs). hIPPYlib implements state-of-the-art scalable algorithms that exploit the structure of the problem, notably the Hessian of the log posterior. The key property of the algorithms implemented in hIPPYlib is that the solution is computed at a cost, measured in forward PDE solves, that is independent of the parameter dimension. The mean of the posterior is approximated by the MAP point, which is found by minimizing the negative log posterior. This deterministic nonlinear least squares optimization problem is solved with an inexact matrix-free Newton-CG method. The posterior covariance is approximated by the inverse of the Hessian of the negative log posterior evaluated at the MAP point. This Gaussian approximation is exact when the parameter-to-observable map is linear; otherwise it can serve as a proposal for Hessian-based MCMC methods. The construction of the posterior covariance is made tractable by invoking a low-rank approximation of the Hessian of the log likelihood. hIPPYlib makes these advanced algorithms easily accessible to domain scientists and provides an environment that expedites the development of new algorithms. hIPPYlib is also a teaching tool that can be used to educate students and practitioners who are new to inverse problems and to the Bayesian inference framework.

ABSTRACT. In recent years, cloud computing has emerged as a popular new architecture for the Open Source scientific computing community, and industrial users. It has substantial benefits for end-users over traditional computing resources, for example, task-tailerd resources and negligible “queue time” with potentially low total time-to-solution.

In the MSO4SC project we have developed an Open Source multi-cloud-HPC framework allowing easy access and use of FEniCS on a diverse collection of HPC resources, both traditional supercomputers and cloud computers. Here we investigate an adaptation of MSO4SC to the Google Compute Engine (GCE).

We show performance analysis in terms of metrics comparing against a traditional supercomputer. MSO4SC-FEniCS unlocks scientific reproducibility for HPC. The Open Source code and computation are made easily reproducible, modifiable and inspectable by the Open Source MSO4SC cloud framework, here with a web interface of a familiar desktop environment.

ABSTRACT. In this work a numerical solution of systems of coupled non-linear parabolic equations using mixed finite elements in FeniCS is presented. The equation system is solved in space by a finite element method, whereas a backward finite difference method is applied in time, resulting in a completely implicit numerical scheme [1]. Three finite element formulations: standard (ST), mixed (MX) and dual-mixed (DM) are compared in terms of their convergence and performance [2]. The computational implementation is performed in Python based on the FEniCS software library [3].

The numerical model is applied to a case study of fine migration during a low salinity water flooding process (LSWF). This model consists of a single-phase flow and a multicomponent transport model in porous media to describe fines attachment-detachment, migration and clogging processes observed in low salinity water flooding experiments [4]. Numerical solutions obtained with FEniCS implementation are compared with those obtained with commercial software COMSOL Multiphysics [5].

References

[1] Allen, M. B., Herrera, I., and Pinder, G. F. Numerical modeling in science and engineering. John Wiley & Sons, 1988.

[2] Chen, Z., Finite Element Method and their Applications, Springer-Verlag Berlin Heidelberg, 2005.

[3] Logg, A., et. al., Automated Solution of Differential Equations by the Finite Element Method, Springer-Verlag Berlin Heidelberg, 2012.

[4] Coronado, M. and D ıı az-Viera, M. A., Modeling Fines Mobilization and Permeability Loss in Lab Cores by Low Salinity Water Flooding, Journal of Petroleum Science and Engineering vol. 150, pp. 355-365, 2017.

[5] COMSOL Multiphysics, User’s Guide Version 3.4, COMSOL AB (2007)

ABSTRACT. The multigrid method is among the most efficient methods for solving linear systems that arise from the discretisation of elliptic partial differential equations. Multigrid has a strong theoretical basis, offering optimal algorithm complexity of O(N) for linear systems with N unknowns, and is amenable to parallel implementations. However, a commonly held view in industrial communities is that multigrid is not well-suited for large-scale systems that are characterised by complex engineering geometries. It is observed that multigrid performance degrades for poor quality finite element meshes, possibly leading to a failure to converge [1]. With highly complex geometries it is inevitable that cell quality will be lower in some regions, and improving cell defects can be time-costing. Therefore, finding ways to overcome poor solver performance in presence of small number of low quality cells is of great value.

We demonstrate in the context of Galerkin geometric multigrid with simple examples how multigrid smoothers can fail, locally, to eliminate high frequency components of the residual error in regions of low cell quality. This has a severe impact on the rate convergence of the solver. Building on the understanding of the cause of degraded performance, we propose a smoother that addresses the performance and robustness issues in presence of a small number of low quality cells. The proposed multigrid smoother consists of two steps: first, a normal smoother, e.g., Gauss-Seidel, is applied on the whole domain; then, a local residual correction in regions with low quality cells is applied using a direct solver on very small systems. The approach is related to a subspace correction [2] on a small number of local subdomains. The combined smoother is a form of an additive Schwarz-type domain decomposition solver.

For geometric multigrid, a hierarchy of geometric grids are exploited. The low quality regions can be selected directly by setting a threshold on some mesh cell quality. It is found that the local correction scheme also works well in algebraic multigrid (AMG) which includes both classical AMG and smoothed aggregation AMG. In particular, the prolongation operator of smoothed aggregation AMG is reconstructed by using the shifted smoothing procedure to preserve the local property that high frequency error only appears in the regions of low quality cells.

We demonstrate the performance of the method for a range of test problems for Poisson equation and linear elasticity based on FEniCS implementation. Tests on meshes with low quality regions indicate that the local correction restores the multigrid convergence rate to that observed for high quality meshes.

The approach is promising for making multigrid methods robust with respect to mesh quality. For geometric multigrid, the coarse level meshes of complex geometries will inevitably contain some low quality cells, and the degraded convergence that this causes can be overcome by proposed approach. For algebraic multigrid, the approach is particularly appealing at extreme scale computation, e.g., contact mechanics, which is well suitable for extreme parallelisation.

References

[1] C. N. Richardson, N. Sime, and G. N. Wells. Scalable computation of thermomechanical turbomachinery problems. Finite Elements in Analysis and Design, 155:32-42, 2019.

[2] J. Xu. The method of subspace corrections. Journal of Computational and Applied Mathematics, 128(1):335- 362, 2001.

ABSTRACT. Mixed-dimensional partial differential equations (PDEs) are equations coupling unknown fields defined over domains of differing topological dimension. Such mixed-dimensional PDEs naturally arise in a wide range of fields including geology, bio-medicine, and fracture mechanics. Mixed-dimensional models can also be used to impose non-standard conditions over a subspace of lower dimension such as a part of the boundary or an interface between two domains, through a Lagrange multiplier. Finite element discretizations of mixed-dimensional PDEs involve nested meshes of heterogeneous topological dimension. The assembly of such systems is non-standard and non-trivial especially with regard to the terms involved in the interactions between the different domains. In other words, automated solution of mixed-dimensional PDEs requires the design of both generic high level software abstractions and lower level algorithms.

While external FEniCS-based packages have been developed for simulating mixed-dimensional problems, there is also an important need for embedding these features as an intrinsic part of the FEniCS library. We introduce an automated framework dedicated to mixed-dimensional problems addressing this gap, which has already been used with diverse applications ranging from reservoir simulations to ion concentration dynamics models.

In this talk, we give an overview of the abstractions and algorithms involved, with a focus on newly implemented features including parallel submeshes extraction and interface formulation handling. The introduced tools will be illustrated by concrete examples of applications in biomedicine.

ABSTRACT. Electrical conduction in brain tissue is commonly modeled using classical bidomain models. These models fundamentally assume that the discrete nature of brain tissue can be represented by homogenized equations where the extracellular space, the cell membrane, and the intracellular spare are continuous and exist everywhere. Consequently, they do not allow simulations highlighting the effect of a nonuniform distribution of ion channels along the cell membrane or the complex morphology of the cells.

In this talk, we present a more accurate model for cerebral electrodiffusion with an explicit representation of the geometry of the cell, the cell membrane and the extracellular space. To take full advantage of this framework, a numerical solution scheme capable of efficiently handling three-dimensional, complicated geometries is required. We propose a novel numerical solution scheme using a mortar finite element method, allowing for the coupling of variational problems posed over the non-overlapping intra and extracellular domains by weakly enforcing interface conditions on the cell membrane. This solution algorithm flexibly allows for arbitrary geometries and efficient solution of the separate subproblems. Finally, we implement the numerical scheme in the new mixed dimensional framework allowing for 3D- 3D coupling in FEniCS.

ABSTRACT. I will present the Unified Continuum model that we have been developing and implementing in recent years as part of our effort toward applying our stabilized finite elements method to solids and fluids altogether. Recently, we have focused more on applications in biomedical engineering and biomedicine, an activity that, internally, we refer to as \emph{Digital Uman - Unified human modeling}, so my presentation will highlight progress in this direction. The idea that we are pushing is that of developing a framework that seamlessly applies to several different applications with minor adaptations, tackling problems related to the flow of blood and the deformations typical of biological tissue. Alongside fluid flows and fluid-structure interaction, we have also been working on multi-physics applications and coupling of FEniCS with other software. According to the plan, we are involved in a number of projects targeting applications in several areas of biomedicine, often very close to clinical applications. I will present at least two of them: simulation of cardiac radiofrequency ablation and coupling of FEniCS with a 0D model of the circulatory system. In biomedical applications, modelling is mission-critical and an error can yield dramatic consequences; for this reason, we are promoting a culture of reproducibility and falsifiability in our work that allows other scientists to easily inspect and test our results. On this topic, I will present our framework developed as part of the MSO4SC EU project that is currently based on Google Cloud technology.

ABSTRACT. When complex mechanical behaviors are involved, the high-fidelity mesoscale/multiscale methods for soft tissue constitutive models are predictive but are also computationally very demanding. Traditional phenomenological models have simple forms based on physical insight, but they often lack the ability to be predictive. In this work, we considered a recently developed predictive structural constitutive models for native and exogenously crosslinked collagenous soft tissues. We have constructed a small and shallow neural network (NN) as an efficient surrogate model for the detailed but computationally costly structural model. The NN model can faithfully represent the response of the planar soft tissues with a range of meso-scale fiber structures. To integrate tissue geometry more readily, the finite element simulation with the NN model was implemented in tIGAr, a python library for isogeometric analysis (IGA) using FEniCS. Thanks to the expressiveness of the NN model, it paves the way for applying machine learning methods for simulation of soft tissues and identification its mechanical properties to much broader applications.

ABSTRACT. Heart-on-a-chip devices can be used to control and monitor the cellular microenvironment of human cardiac microtissue, showing great potential as powerful cardiotoxicity screening platforms. Human induced pluripotent stem cell derived cardiomyocytes are combined with microfluidic technology to derive a tissue mimicking adult human ventricular dimensions, physiology and uniaxial beating properties, which can be used for pharmalogical studies and disease modeling. Within such microphysiological systems, indicators of cardiac function, such as electrophysiology (action potential, calcium transient) and beating physiology (beat rate, contraction velocity) can be readily measured using image analysis techniques. However, important mechanical properties such as myocyte contractility, tissue stress and tissue stiffness, cannot easily be revealed using current techniques. We thus propose a numerical model simulating the dynamic mechanics of the microtissue.

For our model we combined previously published models for cardiac electro-chemical and mechanical properties, embedded with a geometrical representation of the tissue chamber, to quantify the passive and active properties of the system. We used a zero-dimensional point model of cardiac myofilaments developed by Rice et. al. to calculate the force generated by single myocytes. This was coupled with a three-dimensional continuum model for myocardial mechanics given by Guccione et. al. Together this gave a set of non-linear partial differential equations describing the stress and deformation, which we set up to find an equilibrium solution and calculate the resulting Cauchy stress. We imposed appropriate boundary conditions and solve the system using FEniCS – both the framework for spatially dependent nonlinear boundary value problems, and for time-dependent initial value problems – and evaluated the system using data provided from experimental pharmacology studies.

Our results provide a quantification of the Cauchy stress developed through active myocardial contraction through the domain, and relative changes under differing simulated conditions. We compared these results with experimental measurements of contraction forces, using the same points for measurements. In this dual analysis we observed a comparable relative increase in Cauchy stress and microscopic measurements of contraction force, respectively, correlated with increasing drug doses. This combined in silico / in vitro system give us a platform for the quantification of internal mechanical properties, providing a framework for understanding the physiology complexity of micro-organ tissues. Such a system can give new insight in cardiac mechanics in itself, and might also be important for further development of appropriate mathematical models.

ABSTRACT. Reaction-driven cracking is a coupled process whereby fluid-induced reactions drive large volume changes in the host rock which produce stresses leading to crack propagation and failure. This in turn generates new surface area and fluid-flow pathways for subsequent reaction in a potentially self-sustaining system. This mechanism has has been proposed for the pervasive serpentinization and carbonation of peridotite, as well as applications to mineral carbon sequestration and hydrocarbon extraction.

The key computational issue in this problem is implementing algorithms that adequately model the formation of discrete fractures. Here we present models using a discontinuous finite element method for modeling fracture formation (Radovitsky et al., 2011). Cracks are introduced along facets of the mesh by the relaxation of penalty parameters once a failure criterion is met. It is fully described in the weak form of the equations, requiring no modification of the underlying mesh structure and allowing fluid properties to be easily adjusted along cracked facets. To develop and test the method, we start by implementing the algorithm for the simplified Biot equations for poro-elasticity using our model assembler TerraFERMA, built on the finite element library FEniCS.

We consider hydro-fracking around a borehole (Grassl et al., 2015), where elevated fluid pressure in the poro-elastic solid causes it to fail radially in tension. We discuss issues arising from this method, including the formation of null spaces and appropriate preconditioning and solution strategies. Initial results suggest that this method provides a promising way to incorporate cracking into our reactive fluid flow models and future work aims to integrate the mechanical and chemical aspects of this process.

ABSTRACT. Pneumatic storm surge barriers promise several economic and environmental advantages over traditional coastal structures, such as sea walls. Such barriers, made of synthetic rubber, can be stowed and deployed, and designed for a shorter lifespan considering uncertainties surrounding future storm behavior. For this type of barrier to become a feasible option in coastal protection, reliable analysis methods must be developed that can capture the nonlinear finite deformation and the fluid-structure interaction (FSI) between the storm surge and the membrane.

This paper addresses the finite-element modeling of the pneumatic barrier subject to hydrodynamic loading. After a review of appropriate numerical modeling techniques, a preliminary FSI model of a pneumatic storm-surge barrier in a single-phase flow is presented. The pressurized barrier, which represents a deformable portion of the boundary of the fluid domain, is modeled following the work of [1] as a deformable membrane enclosing an ideal gas. The fluid is modeled via the incompressible Navier-Stokes equations using a Continuous Galerkin discretization with SUPG stabilization. Arbitrary Lagrangian-Eulerian strategies are employed to ensure conformity of the fluid mesh with the structural deformation. The end goal of the modeling effort is to study the response of the membrane to loads imparted by free surface flows (e.g. breaking waves). To that end, the paper also investigates coupling the model with interface capturing schemes such as the conservative level set method [2] and more sophisticated FEniCS packages such as the two-phase Navier Stokes solver Ocellaris [3]. This paper is of relevance to the scholar interested in pneumatic membrane structures and fluid-structure interaction.

ABSTRACT. We present an efficient preconditioner for solving the Helmholtz equation with high-order Lagrangian Finite Element Methods (FEM). The proposed method is based on a two-level Schwarz Domain Decomposition Method (DDM) with plane-wave enrichment at the interfaces, combined with a signal processing tool to estimate dominant wave directions.

The solution of the Helmholtz equation is susceptible to the 'pollution effect, which induces a rapid increase in numerical error as the frequency grows. This creates a stringent requirement for the mesh size so that the accuracy of the solution is acceptable. The use of high-order methods can alleviate this requirement, leading to smaller matrices for the same precision. However, as the polynomial degree increases, the spectrum of the resulting matrix rapidly deteriorates. Therefore, the need for an effective preconditioner becomes even more evident.

Because of the suitability for working on parallel computers and demonstrated convergence capabilities, domain decomposition based preconditioners appear as an appealing alternative. However one-level methods are not scalable with respect to the number of subdomains; information is exchanged only locally between adjacent subdomains. Also, standard coarse-grid methods are not suitable for time-harmonic wave problems; the coarsest mesh is still fine to resolve the oscillatory nature of the waves. In this work, hence, an alternative means of creating the second level is investigated.

The construction of our coarse space is motivated by the geometric optics ansatz, and it is carried out during a pre-processing step. We assume that a wave field is a weighted superposition of plane-waves, and their directions change only slightly with frequency. Then we solve a cheaper low-frequency problem to extract the directions and the weights, using the 2nd order numerical micro-local analysis (NMLA) algorithm. We perform a careful selection of the dominant plane waves to enrich the interfaces of each subdomain, improving the scalability of the original problem and ensuring load balancing. While the actual problem is discretized with high-order Lagrangian finite elements, the low-frequency problem, on the same mesh, is discretized by linear Lagrangian elements.

We perform numerical experiments to validate our proposal. For two-dimensional scattering problems, we observe a nearly constant iteration count for the preconditioned GMRES, with an increasing number of subdomains ($N_s = 16 \cdots 128$). All numerical experiments have been performed using the \textit{FEniCS} library.

ABSTRACT. The behaviour and segregation of magma from the mantle into volcanoes is an open question. Modelling the liquid magma in the permeable network of pores between grains gives rise to four coupled nonlinear partial differential equations, conserving solid deformation, compaction pressure, magma porous flow and porosity of the mantle.

We demonstrate how the new features of DOLFIN-X facilitate the scalable solution of the finite element formulation of this system. In particular: high order geometry representation, block matrix assembly in the MatNest data structure of PETSc and automatic detection of block reassembly requirements.

Finally we show results from preliminary simulations inspired by experiments demonstrated in laboratory. We highlight the scalability of the system assembly and computation of the numerical solution. Additionally we monitor the impact of a weakly enforced free slip boundary condition on the conditioning of the underlying linear system.

ABSTRACT. In the field of tissue engineering and regenerative medicine, electromagnetic stimulation of biological samples has received much attention in the last decades. Some applications such as electrical muscle stimulation are already used for medical therapy. Still, the interaction of different kinds of cells with external electromagnetic fields is a current topic of research. Examples are the stimulation of stem cells to guide their differentiation or the stimulation of bone or cartilage cells to enhance the production of extracellular matrix proteins. Besides a few commercial solutions most of the experimental devices are custom-made and operated by medical scientists. In order to link the medical/biological findings to the applied stimulus, numerical computations of the electric field distribution and other quantities of interest are required. The simulation results also serve to optimize experimental parameters such as electrode spacing and applied voltage. To date, most of the researchers use commercial packages such as COMSOL Multiphysics R and do not reveal details of their modeling approach that are crucial for reproducibility. We present our open-source package EMStimTools [1] that automates the solution of relevant partial differential equations (PDEs) and summarizes the modeling assumptions in a comprehensible YAML input file. A geometry including subdomains as well as boundaries is built and meshed in SALOME and converted to a FEniCS data format to solve relevant PDEs on it. The result is exported in the xdmf file format and can be visualized by ParaView.

[1] https://github.com/j-zimmermann/EMStimTools/releases/tag/v0.1.3.dev0

ABSTRACT. We study the stability of a class of three-field mixed finite element methods for nonlinear elasticity called CSFEMs near regular solutions. To this end, we write suitable material-independent and material-dependent inf-sup conditions. The material-independent condition is only a necessary condition for the stability while the material-depend one is sufficient as well. By considering several combinations of first-order and second-order elements, we study the validity of these inf-sup conditions by using singular values associated to the matrix form of the inf-sup conditions.

ABSTRACT. In this work we ran a suite of finite element simulations to estimate the anisotropic elastic moduli of a nickel alloy from high energy x-ray diffraction measurements. In the experiment, a polycrystal specimen was scanned in an unloaded state to determine a spatial map of crystal orientations on the interior. The specimen was then loaded and scanned in situ to determine grain centroids and average grain strain tensors at several loads. To estimate the moduli, we ran a series of finite element simulations of the polycrystal specimen. Each grain was modeled as a subdomain of the mesh, and its orientation was initialized from the x-ray data. The material was modeled using anisotropic linear elasticity, and we varied the moduli systematically. The simulation results were processed to compute grain-averaged strains and then compared with experimental data. We present the optimal moduli and show how the error varies near the peak.

ABSTRACT. A new computational approach for the general numerical simulation of fluid flows is presented and demonstrated for compressible natural convection flows with real fluid properties. The Characteristic-Based Split finite element method is generalized, and abbreviated as GCBS, in the present work to introduce both pressure-expansion effects in addition to thermal-expansion effects on density within the projection step, in contrast to recent works that preserve only one of those two expansion effects. Moreover, we keep the fluid in thermodynamic equilibrium to support a real equation of state and fluid model. The method retains the full Navier-Stokes equations, transitioning smoothly between compressible and incompressible flows, without special artificial-compressibility methods, special low-Mach-number methods, or other modeling treatments. Neither isothermal or isobaric assumptions are made. The non-linear solution obtained is validated by the linearized solution during the early time period of the simulation.

Natural convection in a cavity is explored for its convenience in testing across known benchmark flow regimes, from classical Boussinesq flow to the full Navier-Stokes equations with real fluid models. Results for natural convection in a cavity across Rayleigh numbers, non-dimensional adiabatic temperature gradients, and temperature differences, are obtained across flow regimes with the GCBS method. Classical Boussinesq, Thermodynamic Boussinesq, low-Mach number, and real-gas flows are presented, demonstrating the generalization of the method across flow regimes. Thus, our contribution is to generalize the Characteristic-Based Split method to capture the full Navier-Stokes equations while demonstrating the capability with compressible natural convection.

The FEniCS framework was actively employed in for the theoretical development of the GCBS scheme, as well as for the computational development of a scalable CFD code. The flexibility provided by FEniCS and its implementation of the Unified Form Language (UFL) were instrumental for developing the GCBS scheme through incremental trial-and-error exploration of various possible versions of operator-splitting into a series of linear predictor and corrector steps to solve the nonlinear compressible Navier-Stokes equations. The theoretical development of the GCBS scheme, with the aid of FEniCS, began with the Chorin-style operator split as a foundation upon which, through gradual generalization, we built the original CBS scheme, and then the GCBS scheme. The FEniCS framework enabled rapid implementation of parallelization through MPI, abstraction of two- and three-dimensional flows, and easy exploration of the effects of various element types and interpolating polynomial orders. Moreover, the FEniCS framework also facilitated the validation of the GCBS scheme through the computation of various challenging fluid benchmark problems.

ABSTRACT. Here we provide an overview of global-scale mantle convection models constructed with FEniCS. We explain some of the challenges in modeling geodynamic phenomena, including discontinuous boundaries, highly non-linear rheology, and extreme variability in length and time scales of the Earth system. We demonstrate how FEniCS can be used to address such challenges and show that our simulations produce geologically reasonable solutions. We conclude by suggesting areas of development where FEniCS may prove useful in creating the next generation of geodynamic models.