ABSTRACT. We consider the numerical solution of the regularized fourth order total variation flow problem based on an implicit discretization in time and a $C^0$ Interior Penalty Discontinuous Galerkin ($C^0$ IPDG) method in space. In particular, we prove coerciveness and boundedness of the associated $C^0$ IPDG form with respect to a mesh dependent $C^0$ IPDG-norm of the underlying func- tion space. The main result is an a priori error estimate of the global dis- cretization error in the mesh dependent $C^0$ IPDG-norm. A documentation of numerical results implemented in fenics is provided illustrating the performance of the $C^0$ IPDG method and confirming the theoretical findings. The below mentioned figure suggest that the model can be used for surface regularization.

ABSTRACT. In this talk, we present compatible finite element spatial discretizations for three-dimensional equations relevant to numerical weather prediction. Our motivation revolves around developing a scalable parallel implicit solver, ideally scaling to several thousands of processors, for the full compressible system of PDEs used in operational weather forecasting. This work is a subject of ongoing development with the UK Meteorological Office, as part of the ``Gung-Ho'' dynamical core project. Traditional staggered finite difference procedures first eliminate the velocity field to produce an elliptic pressure problem. However, this is not so straightforward to accomplish within the framework of compatible finite elements. To remedy this, we use a procedure known as ''hybridization.''

The hybridization of mixed methods for second-order elliptic problems is well-known. We extend these techniques for use in the compressible equations, which permits the use of element-level static condensation. This approach mimics the elimination strategy of previous generation finite difference codes and produces a new, reduced problem defined on the mesh interfaces. We use software built on top of Firedrake to proto-type these hybridization methods, using the technology developed in this paper to facilitate both static condensation and local recovery operations.

ABSTRACT. Multigrid is an iterative method known for its efficiency for solving elliptic PDE problems. However, its performance degrades with poor quality finite element meshes, possibly leading to a failure to converge. When meshing geometrically complex domains for engineering applications, clusters of poor quality cells are not uncommon. Ideally, the mesh quality would be improved at generation time, but this may be difficult for complex geometries. Moreover, mesh generation may be performed by one team, and the analysis performed by another. The delay caused by the analyst returning a mesh to the model creator for improvement can be substantial and unacceptable in a design and analysis process. Finding techniques to overcome poor solver performance in the presence of a very small numbers of low quality cells is appealing in order to increase the robustness (and by extension acceptance) of multigrid methods in engineering practice. We propose a smoother than can overcome the degraded performance of multigrid when a small number of low quality cells are present. Regions with a small number of low quality cells are identified, and on these very small regions (usually just a handful of cells) a direct solver is used as the smoother. The approach has been tested on a range of problems, and implemented in the FEniCS/PETSc Galerkin geometric multigrid framework. Figure 1 shows results for a two-dimensional Poisson problem. A localised region containing some poorly formed cells can be observed at the centre of the domain (see Figure 1a). The relative residual against iterative count for a two-level implementation using Gauss–Seidel with the local correction and a standard Gauss–Seidel smoother are shown in Figure 1b. The convergence rate with the standard smoother on the low quality mesh is clearly poor, whereas with the local correction the residual converges must faster. With the correction, the convergence rate observed for a high quality mesh is recovered. The results indicate the smoother with localised correction has the potential to dramatically improve the speed and efficiency of multigrid for complex engineering applications. The application of robust smoothers for meshes with regions of low quality opens new possibilities for Galerkin geometric multigrid applied to complex geometries, where coarse grid representations, in particular, will inevitably have regions of low cell quality when representing geometric feature of a domain with low resolution meshes.

ABSTRACT. In this talk we will present stable discretizations and robust preconditioners for the multiple-network poroelasticity equations. In order to obtain robustness in the whole parameter space, we will present new formulations obtained through the linear combination of the fluid pressures and the associated preconditioners. Numerical simulations are performed using FEniCS.

ABSTRACT. In this talk, we will describe our implementation of the Optimized Restricted Additive Schwarz method specially tailored for finite element discretizations of the curl-curl time-harmonic Maxwell's equations with Nédélec elements. Furthermore, we will discuss the choice of transmission conditions, a key ingredient in domain decomposition methods, and provide some numerical experiments to highlight their impact on the convergence of the methods considered here.

ABSTRACT. We consider non-smooth semilinear elliptic PDEs and present a new algorithmic approach for solving optimal control problems constrained by these PDEs. In this talk, we focus exclusively on optimization problems where the non-smooth term in the PDE constraint can be represented by combinations of the Lipschitz-continuous functions abs(), max() and min().

The key idea of the optimization method under consideration is to locate stationary points by piecewise linearization and appropriate decomposition of the original problem into several smooth branch problems. Each of these branch problems can be solved by classical means. Subsequent exploitation of information given by the respective dual variables leads to the next branch and thus to a successive reduction of the function value.

Numerical results are presented for two-dimensional examples.

ABSTRACT. Adsorbed Natural Gas (ANG) tanks are able to store gas at low pressure and high density based on the adsorption of molecules on the surfaces of porous materials. Despite advantages in comparison with other storing methods, fast adsorbing cycles are hindered by the dependence of adsorption on temperature. Adsorption is an exothermic reaction, the equilibrium adsorption density decreases as temperature increases and the thermal conductivity of typical adsorbent beds is very poor. In this work, it is sought to improve heat dissipation inside ANG tanks by distributing non-adsorbent materials presenting suitable thermal properties across adsorbent beds. Optimised placement of such materials is determined by the Topology Optimisation Method (TOM), which is implemented in FEniCS, aided by Dolfin-Adjoint. The optimisation problem consists in the maximisation of the total mass of gas admitted by a tank connected to an infinite reservoir at fixed pressure during certain time length. In order to avoid numerical instabilities and mesh dependency, a filter based on a modified Helmholtz equation is applied to the design variables field. Sections of axisymmetric and prismatic tanks are optimised considering different non-adsorbent materials. Tanks are considered to be subjected to several thermal boundary conditions, such as natural convection, forced convection and contact with cold fingers.

ABSTRACT. A trend of close integration of Computer-aided Design (CAD) and Computer-aided Engineering (CAE) has been seen for commercial and open source software packages. CAD platforms extend their functions from geometry creation to engineering simulation and product life management, while CAE software, for physical system simulation in a narrow sense, has strengthened the geometry creation and meshing capability. However, there are gaps for an automated engineering design/optimisation from geometry prototyping to Computer-aided Manufacturing (CAM). In particular, the changed geometry topology, i.e. the way surfaces forming a solid geometry, will invalidate meshing and boundary condition setup. In this paper, an automated engineering design pipeline has been developed in Python, bringing together popular open source CAE modules. Fenics and OpenFOAM solvers are utilised to solve multi-physics partial differential equations. Moreover, Fenics and OpenFOAM solvers are integrated into FreeCAD. Thereby, a one-stop design process from geometry building to analysis is provided for both GUI and scripting users. User operation in FreeCAD GUI can be recorded into Python code, which assists experienced engineers in building up automatic design pipeline in a later stage. Furthermore, two user cases of this pipeline, multiscale modelling of heat transfer in engine brush seal and optimisation of cooling performance for engine components, are illustrated. Although it is a long way towards intelligent engineering design, the pipelined process can reduce engineer’s time to numerically evaluate the optimised design based on a prototype geometry.

ABSTRACT. This work presents a methodology to assist the design of this class of innovative bistable structures that work under nonlinear buckling behavior by employing the Topology Optimization Method which satisfies the nonlinear optimization constraints. The proportional load factor function is applied as a constraint to impose structural bistability behavior as design requirements. The nonlinear buckling solution is obtained by employing the Generalized Displacement Control Method as a path-tracing algorithm.

ABSTRACT. We present recent results in both development and CFD applications of our high-performance simulation software based on FEniCS. Our computational model is based on the incompressible Navier-Stokes equations, with the option of a cheap wall model for high Reynolds number, which we solve numerically with a Galerkin least-squares stabilized finite elements method. Automation plays a fundamental role as the whole formulation is implemented in UFL, yielding more concise, more readable and less error-prone code thanks to automatic source-to-source compilation provided by FFC. Building up on automation we can foresee the development of a standard HPC flow solver where a simulation can be specified entirely through a configuration file that is input to a general binary. This design, made possible by FEniCS' high-level API, can be very appealing to scientists of other fields that are interested in exploiting the capabilities of modern HPC resources without having to deal with too much computer science in the process. Again to make our software more easily accessible, we recently simplified deployment by providing a self-contained environment in the form of a Singularity image. This container technology is available on a range of architectures, from personal computers to HPC clusters, and makes deployment and installation a matter of a one-line command. An interesting possibility is then to offer simulations as a service: this is being realized in the H2020 project MSO4SC where we contribute to the realization of a platform allowing users to shop simulation cases that use our software as a black box running on potentially a variety of different clusters. In such an agile development environment, applications can flourish. As far as applications are concerned, we will present results in two directions. On the one hand, we will discuss our recent contributions to two high-profile benchmark workshops -- NASA's HiLiftPW-3 and HiOCFD5 -- where we show that our methodology not only keeps up with, but even surpasses our main competitors in the field, especially after recent development yielded a 100 times speedup compared to our previous version. On the other hand, we will show an implementation of a multi-physics solver, built on top of Unicorn/FEniCS, that simulates a cardiac surgery procedure and shows how complex solvers can be easily built thanks to software automation and simple design techniques. All the software and the projects that we develop are based on open source software.

ABSTRACT. Compressible and nonisothermal effects are often ignored when modelling non-Newtonian fluids. How- ever, flows of viscoelastic materials often occur under extreme conditions with large temperature and pres- sure gradients where the assumption of constant density, viscosity and temperature is no longer suitable. For example, fluid compressibility can determine the stability and load-bearing capacity of lubricants used in journal bearings. Taylor-Galerkin pressure correction schemes coupled with Discrete Elastic Viscous Stress Splitting (DEVSS) stabilization enable accurate solutions to be generated for a wide range of flow parameters [2, 4]. The governing equations are derived and some numerical results for 2D flow problems are presented. In particular, the flow between eccentrically rotating cylinders (journal bearing lubrication) is considered for a range of Reynolds, Weissenberg and Mach numbers. The influence of compressibility and viscoelasticity on flow is assessed.

ABSTRACT. Mathematical models for multiphase flow in porous media are usually based on extensions of Biot’s theory; governing equations are often formulated directly at the macroscale, or obtained by upscaling of microscale conservation equations. This second option is the most suitable because it allows to properly define the connection between micro- and macro- scales and to know all assumptions and simplifications needed to obtain the final system of macroscale conservation equations. Averaging theories provide well established frameworks to develop mathematical models for multiphase systems at any scale of interest, and assure a rigorous connection between microscale and larger scale conservation equations. However, once mathematical upscaling is achieved, a number of constitutive relationships are typically required to obtain a solvable system of equations; very often, identification of relevant closure relationships (e.g. pressure-saturation relationships, phase relative permeabilities, etc.) remains therefore an open issue. With this overall aim, a microscale mathematical model for biphasic flow has been developed. The formulation consists in a set of partial differential equations coupling Cahn-Hilliard model, with Navier-Stokes equations. The model will be validated using micro-fluidic studies of literature. Once validation achieved, a set of numerical experiments will be performed to identify macroscopic closure relationships for target applications of our research team (tumor growth modeling and enhanced oil recovery (EOR)).

ABSTRACT. The application of the covariance principle to physical laws in an arbitrary reference system establishes its independence from an observer change, due to the use of time as an additional coordinate, leading to the four-dimensional framework. The covariance principle of General Relativity is formulated as follows: "the laws of nature must be written in a form which is the same in any four-dimensional system of coordinates (or, as one says, in a covariant form)" as mentioned by Landau. This new dimension is integrated in physical equations, which should transform like frame indifference equations. This is established by the use of 4D tensors and 4D operators, which will make all transformation frame indifferent. In three dimensional study, the variational principle is obtained by multiplying the partial differential equation by a test function, then integrating the resulting equation over a three-dimensional domain, which is defined clearly in FEniCS codes. In the present work to obtain the variational principle, we will use the covariant derivative instead of the partial derivatives, which will ensure the observer-independent form of our variational principle. Then we will introduce the four-dimensional momentum-energy tensor. It is widely used in relativity. It will be useful to study stress and strain in materials by writing the constitutive models in a covariant way.

ABSTRACT. In order to guarantee obtaining frame indifferent formulations, we look to express the covariance principle of physics laws. To apply this principle the theories of relativity are derived in a space-time domain. This purpose makes the 4D formalism relevant to our interest. Furthermore, to obtain behavior models, a 4D formulation of the inequation of Clausius Duhem is a way to validate the second principle of thermodynamics. This is supposed to systematically lead to obtain satisfying behavior models. The main objective is to examine the thermomechanical behaviour obtained from the 4D approach. However at this stage, the thermal model is studied as a simplied case. The need of a solver taking into account the time coordinate that we introduce is, from this point forward, essential. The use of FEniCS to fulfill this role is convenient. However, since FEniCS can only handle 3D meshes, we are led to consider the case of evolution of temperature based on Fourier and Cattaneo models with 2 dimensions dedicated to the space and 1 to the time.

ABSTRACT. The simulator Flow123d of transport processes in fractured media is shortly presented. With aim to adopt FEniCS in existing models and in further development of more complex THMC couplings we present a reduced continuum-fracture model for the poroelasticity and preliminary numerical experiments using FEniCS. We consider a poroelasticity compression problem for a material with a thin fracture modeled as an porous and elastic material with different material parameters. The Biot's model with full resolution of the fracture in 2d and a reduced model with discrete 1d fracture are compared numerically. We used the Mandel's problem to study stability of four different FEM spaces for a nearly incompressible case.

ABSTRACT. Shape-optimization with time-dependent PDE-constraints play an important role in many scientific and industrial applications, such as drag minimization over airplanes, optimizing the shape of acoustic horns and maximization of extracted energy from turbine farms. Since time dependent PDEs are computational expensive to solve, gradient based optimization methods combined with the adjoint equations offer an efficient approach for solving such problems. A well designed optimization setup can scale nearly independently of the number of design variables. However, deriving the adjoint equations for a time-dependent PDE is non-trivial.

To cope with this in FEniCS, the extension dolfin-adjoint was created. This extension uses Algorithmic Differentiation (AD), to obtain the functional sensitivity and corresponding adjoint equations. This is done by recording the solution procedure of the state-equations. With this record, it can compute the corresponding adjoint equations.

However, in dolfin-adjoint there is no support shape derivatives. These derivatives are needed to evaluate the functional sensitivity. This has been implemented in FEMORPH, which is a FEniCS-extension that computes shape-derivatives. For this tool to work, it needs to be supplied with both the state and adjoint solution. Hence, the corresponding equations needs to be solved.

This poster will present how to combine these two tools to achieve AD of shape-optimization problems with both time-dependent and time-independent PDE-constraints. We will illustrate the implementation through an example.

ABSTRACT. A stabilized finite element method for compressible flow has been implemented within the FEniCS framework. This work describes the analytical techniques and numerical methods that have been used to develop the solution algorithm. Results are presented for a validation test case.

The governing equations addressed by this solver are the compressible Euler equations. These equations describe the adiabatic, inviscid flow of a gas. The solver implementation is based on the Streamline-Upwing Petrov-Galerkin (SUPG) method applied to compressible flows. This allows equal-order representation of the solution variables. An additional term is added to the variational formulation to provide the extra stability needed for shock-capturing. This technique is known as the YZ-beta method. The extra term introduces an artificial viscosity in the vicinity of shock waves.

A feature of modelling flows governed by the Euler equations is that a slip condition must be applied at wall boundaries. This is different from flows governed by the Navier-Stokes equations, in which a no-slip condition must be applied. In the new solver the slip condition is imposed as part of the weak formulation, in the style of the Nitsche method. The advantage of the Nitsche method in this context lies in the avoidance of local co-ordinate transformations, which would be needed along the wall surfaces if the slip condition were to be imposed in strong (Dirichlet) form.

A test case has been studied to validate the new implementation. The geometry is that of a diverging nozzle. There exists an exact one-dimensional solution for this problem.

The flow of air within the geometry is simulated. For a particular combination of inlet and outlet boundary conditions a plane shock wave forms at the nozzle midpoint. The flow under these conditions is modelled by the new solver and the location of the shock wave is compared with the analytical solution.

ABSTRACT. The charged particles that constitutes a plasma are subject to electromagnetic forces while themselves inducing electromagnetic fields [1]. In order to study plasma phenomena, one often employs Particle–In–Cell (PIC) simulations, where simulation particles are numerically integrated in time through electro- magnetic fields represented on a mesh [2, 3]. The source terms of the field equations in turn depends on the particles’ positions.

One phenomena which is studied using PIC simulations is the charging of electrically conducting objects such as spacecrafts and their instruments in ionospheric plasma. Figure 1 shows an example of a Langmuir needle probe. The more accurate geometrical representation achieved by using an unstructured mesh is clearly advantageous to the structured meshes traditionally used in PIC simulations. While a couple of FEM–PIC codes already exists [4, 5], we have created a novel code using FEniCS [6] for solving the field equations. The flexibility inherited from FEniCS offers rapid testing and prototyping of new ideas. Moreover, while other FEM–PIC codes uses ad-hoc techniques to obtain the unknown boundary potential on the objects, we treat it as an integral part of the variational problem by imposing constraints using Lagrange multipliers.

For verification the code is run for a spherical object and compared to the well known results by Laframboise [9].

ABSTRACT. I intend to give a brief description of my MSc final project research status, supervised by Johan Jansson at KTH, which includes an illustration of the chosen methodology and structure of the research stages. One of the current principal challenges in Computational Fluid Dynamics (CFD) is the predictability of turbulent separated flow, in particular for a complete aircraft. In [1] the authors describe the results of their adaptive Direct FEM methodology in Unicorn/FEniCS for time-resolved simulation of aerodynamics. They simulate the Navier-Stokes equations (NSE) without including an explicit turbulence or boundary layer model, but they are able to directly incorporate the effects of turbulence and boundary layer in the final outputs. This becomes possible by adding a parametrization of the wall shear stress in terms of a skin friction. For very high Reynolds numbers ( Re > 1e7 ), this leads to an approximation of the skin friction to zero (free slip boundary condition), which in other words means that the dissipative effect of the boundary layer disappears. This is the starting point of my MSc thesis project. However, for slightly lower Reynolds numbers (1e5

ABSTRACT. We consider the numerical solution of singularly perturbed differential equations by finite element methods, and implemented in FEniCS. The simplest model problem is the one dimensional linear reaction-diffusion problem:

−ε^2u’’(x) + b(x)u(x) = f(x) on Ω = (0, 1), (1a)

with the boundary conditions

u(0) = u(1) = 0. (1b)

It is singularly perturbed in the sense that the positive real parameter ε may be arbitrarily small, but, if we formally set ε = 0, then it is ill-posed. When ε is small, the solution may exhibit boundary layers (and, perhaps, also interior layers, depending on the smoothness of b and f ). The numerical solution of (1), and its generalisations to higher-dimensional problems, and to the closely related convection-diffusion problem, presents numerous mathematical and computational challenges, particularly as ε → 0. The development of algorithms for robust solution of (1) is the subject of intense mathematical investigation. Here “robust” means two things:

1. The algorithm should yield a “reasonable” solution for all ranges of ε, including resolving any layers present;

2. The mathematical analysis of the method should be valid for all ranges of ε.

(This is a rather informal explanation of the concept of parameter-robustness. It is sufficient for this presentation. For more precise definitions, see, e.g., [2, 4]). Historically, most of the research on numerical solution of reaction-diffusion equations has focused on finite difference methods, and provable ε-uniform convergence, see, e.g., the canonical references [3, 5]. Of course, there have been many studies of finite element methods, but they are far fewer in number. One reason for this is that finite element methods are usually analysed with respect to the energy norm, which, arguably, underweights the layer terms in the solution [1]. In this presentation, we will describe the numerical solution of (1), and its generalisations, by FEMs in FEniCS. Full consideration will be given to the issue of balanced and non-balanced norms, and to the generation of layer resolving meshes.

References

[1] R. Lin and M. Stynes. A balanced finite element method for singularly perturbed reaction diffusion problems. SIAM J. Numer. Anal., 50(5):2729–2743, 2012.

[2] T. Linß. Layer-adapted meshes for reaction-convection-diffusion problems, volume 1985 of Lecture Notes in Mathematics. Springer-Verlag, Berlin, 2010.

[3] J. J. H. Miller, E. O’Riordan, and G. I. Shishkin. Fitted numerical methods for singular perturbation problems. World Scientific Publishing Co. Pte. Ltd., Hackensack, NJ, revised edition, 2012. Error estimates in the maximum norm for linear problems in one and two dimensions.

[4] H.-G. Roos, M. Stynes, and L. Tobiska. Robust Numerical Methods for Singularly Perturbed Differential Equations, volume 24 of Springer Series in Computational Mathematics. Springer-Verlag, Berlin, 2nd edition, 2008.

[5] G. I. Shishkin and L. P. Shishkina. Difference methods for singular perturbation problems, volume 140 of Monographs and Surveys in Pure and Applied Mathematics. CRC Press, Boca Raton, FL, 2009.

ABSTRACT. Stimulating and recording electrodes are a popular means to record the electrical activity of neural networks as well as to stimulate brain tissue to alleviate the effects of neurological disorders such as Parkinson’s disease. Electrical probes could be improved in terms of the shape and size of the contact sites and their geometrical configuration, for this the electrical fields generated by neural activity and the corresponding signal on contact sites need to be modeled accurately. Volume conductor models provide a valuable tool to estimate the electric fields resulting from neurostimulation or bioelectric sources.

In order to develop an accurate model, properties of the different materials and their interactions need to be considered. Typically, the electrode-electrolyte interface and the capacitive and dispersive tissue properties are disregarded. Specifically, models of recording electrodes assume purely resistive volume conduction.

This study aims to give a comprehensive overview of and comparison between the results obtained by using different models for the electrode-electrolyte interfaces and dielectric properties of the extacellular medium [1,2]. Furthermore, we expand upon earlier implementations of recording electrodes and present an easy-to-use pipeline centered around FEniCS [3] that can be applied to any experimental set-up.

References [1] D.A. Robinson, The electrical properties of metal microelectrodes, Proceedings of the IEEE, 56 (1968), 1065-1071. [2] D.R. Cantrell, S. Inayat, A. Taove, R.S. Ruoff, J.B. Troy, Incorporation of the electrode{electrolyte interface into finite-element models of metal microelectrodes, Journal of neural engineering, 5 (2007), 54-67. [3] M. Alnaes, J. Blechta, J. Hake, A. Johansson, B. Kehlet, A. Logg, C. Richardson, J. Ring, M.E. Rognes, G.N. Wells, The FEniCS Project version 1.5, Archive of Numerical Software, 3 (2015).

ABSTRACT. Computational solution of partial differential equations with interface conditions is a challenge that appears in many applications. Our aim is to develop a general high performance computing framework for mesh conforming interface conditions. In this paper we consider two concrete examples: diffusion in biological tissues and turbulent flow fluid-structure interaction models.

ABSTRACT. A posteriori error estimation is a crucial tool for measuring the discretisation error in finite element methods. A posteriori error estimators allows us to bound the error by computable quantities. Moreover when the estimator is local we can refine only the regions where the error is highest. Consequently we improve the convergence while maintaining a reasonable computational cost.

Here we will study the convergence properties of the Bank-Weiser a posteriori error estimator. Our goal is to show that the Bank-Weiser estimator has convergence properties comparable to that of the residual a posteriori error estimator, regardless of the degree of polynomials in the finite element discretisation. Moreover, due to some pathological cases of problems where the solution can be rough, we want to avoid extra-regularity assumptions (e.g. H^2) on the solution of the PDE.

ABSTRACT. Recent advances in the field of laser-driven ion acceleration have made it possible to resolve the material modification upon irradiation with an unprecedented time resolution. Using picosecond proton pulses, the early stages of track formation and generation of reactive species can now be directly probed. As described in [1], when amorphous SiO2 samples are irradiated with a burst of protons, transient opacity depending on the sample depth is observed. In collaboration with the experimentalists at Queens University Belfast, we are developing a multiscale framework to describe the early stages of ion irradiation of dielectrics. In this poster, we present the finite element implementation of a drift-diffusion model of the non-equilibrium carrier dynamics in SiO2 following irradiation. The transient opacity which follows the proton irradiation is then obtained by considering the temperature dependent free carrier absorption. The qualitative agreement between the model and the experimental results is finally discussed.

References [1] Dromey, B., et al. Picosecond metrology of laser-driven proton bursts Nature communications 7 (2016): 10642. [2] Klaumunzer, S. Thermal-spike models for ion track physics: A critical examination. Med. Phys.. (2008): 293-323. The Really Good Journal, 6 (1969), 501-510.

ABSTRACT. In this work the development of a single-phase slightly compressible flow and multicomponent transport model in porous media at laboratory scale is presented. The development methodology consisted on first establishing a conceptual model where the hypothesis, goals, objectives and limitations of the model are considered. Applying the systematic formulation approach of continuum medium modeling the mathematical model is derived, resulting a partial differential equation system with initial and boundary conditions. The mathematical model is discretized using a mixed finite element method. Finally, its computational implementation is performed on the open source software FeniCS project. The numerical solutions were compared with the commercial software COMSOL Multiphysics. Two case studies are presented: one on the effect of the salinity reduction over the fines migration and permeability loss and the other on dynamic porosity and permeability modification due to the injection of water with microorganisms and nutrients.

ABSTRACT. A finite element discrete fracture model is presented for modeling a fractured porous media and to simulate two-phase flow through multiple fractures and rock matrix, with drastically different petrophysical properties, at a scale such that the fractures could be modeled explicitly. To represent the porous matrix a continuum approach is employed whereas for the fractures a discrete one. By applying an averaging procedure to the equation that governs the fluid flow in the fractures, and based on the concept of crossflow equilibrium between the faces of the walls of the fracture and the matrix, the necessity of representing the thickness of each of the fractures in the computational domain is eliminated, given that it is explicitly considered in the mathematical model, solving the problem of a geometry with a high aspect ratio compared to equi-dimensional models. Thus, the fractured porous media is modeled with mixed dimensions elements, representing the fractures as elements of n − 1 dimensions immersed in a porous matrix of n dimensions by isolated internal boundaries lying on the element edges, where the equations that govern the flow of fluids in the matrix and the fractures are coupled by means of jump and average equations, taking into account interactions between the fractures and the surrounding porous media by an additional source/sink term. The derived flow model in fractured porous media is two-phase, based on the saturation and total velocity. The model is based on a previous model [1] and a conforming discrete fracture model which provide a comprehensive derivation of the lower-dimensional formulation [2]. For the numerical solution is applied a Finite Element Method and its computational implementation is carried out in Python using FeniCS project. Finally, the quality of the results obtained by the discrete fracture model is studied by comparison to an equi-dimensional discretization on a case study in two dimensions.

ABSTRACT. In this poster, we present a domain-specific language (DSL) to facilitate the automatic generation of element-local dense linear algebra kernels for finite element methods. The language, known as Slate, is a compact, Python-embedded DSL for code-generation systems using the unified form language (UFL) as the primary user-facing abstraction for describing finite element problems. Slate takes the discretized equations, expressed in UFL, and provides high-level representations of the local element tensors associated with their respective forms. Equipped with a standard set of operations, expressions describing symbolic operators for localized algebra can be constructed from these representations. A specialized compiler for linear algebra is then used to generate low-level code from these mathematical descriptions.

Key examples are presented in this poster demonstrating both syntax and use. Model applications include the static condensation of both continuous Galerkin problems and hybridized finite element systems, local recovery of discontinuous fields, and localized post-processing techniques. Slate is a developing project and continues to expand in its features as various applications arise. All of its components are fully-integrated in Firedrake, providing access to hybridization and static condensation procedures for several ongoing applications. Future integration within FEniCS is possible and we warmly welcome discussion and ideas.

ABSTRACT. Within the endfoot barrier of astrocytes near the perivascular spaces surrounding the blood vessels, there is a high concentration of the channel membrane protein aquaporin-4. The AQP4 proteins form structures in the astrocytic membranes allowing for highly efficient water transport, and play an important role in mechanisms underlying volume and water homeostasis in the brain. The water transport through AQP4 is driven by osmotic pressure gradients, primarily regulated by movement of ions in the brain tissue.

Typically, astrocyte cell dynamics are modeled using electrodiffusive frameworks, addressing both electrical variables such as potentials, and chemical variables such as ion concentrations. Such models comprise nonlinear numerically stiff systems of partial- and ordinary differential equations. Models that additionally account for mechanical variables such as microscopic fluid flow and hydrostatic pressure are not well explored numerically, leaving modelling of glial dynamics and their role in interstitial fluid flow largely open.

I will present ongoing work on a new finite element scheme for such electrodiffusive models connecting mechanical and electrochemical variables at cellular level, as well as details from the FEniCS implementation. In addition, I will sketch out how the scheme will be applied to computational studies on the effect of glial water dynamics on fluid and interstitial transport.

ABSTRACT. % This template is losely based on the abstract template from Conference % on Geometry: Theory and Applications (CGTA 2015), and uses code % listings style ``anslistings'' from the Archive of Numerical Software. \documentclass[10pt,a4paper]{article} %\usepackage[square]{natbib} \usepackage{amsmath} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{graphicx} \usepackage{algorithm} \usepackage{algorithmicx} \usepackage{algpseudocode} \usepackage{anslistings} % <-- should be provided with this template

\pagestyle{empty} \usepackage[left=25mm, right=25mm, top=15mm, bottom=20mm, noheadfoot]{geometry} % please don't change geometry settings!

%\usepackage[T1]{fontenc} %% get hyphenation and accented letters right %\usepackage{mathptmx} %% use fitting times fonts also in formulas

%%%

\begin{document} \thispagestyle{empty}

% Insert title here: \title{\texttt{miXFEM} - an XFEM toolbox to tackle multiphysics problem with \texttt{FEniCS}}

% Insert authors and afiliations here, speaker in bold: \author{{\bf Mischa Jahn}, University of Bremen, mischa@math.uni-bremen.de}

\date{} % please leave date empty \maketitle\thispagestyle{empty}

% Insert keywords here Keywords: \emph{FEniCS, XFEM, Levelset method}\\

% Insert text here: On the \texttt{FEniCS}'18 workshop, we present a poster summarizing the current status of our work on the library \texttt{miXFEM} - an XFEM toolbox for \texttt{FEniCS}. Mathematically, our approach is based on introducing a physically motivated order of discontinuous features and, thereto, transform multiphysics problems involving discontinuities into a setting suitable to be tackled using automated code generation. Therefore, all arising discontinuities are represented by zero levels of indicator functions, and the usually used finite element spaces \begin{align} V_{{\rm cg},h}^m = \{v_h \in {C}^0 (\Omega) : {v_h}|_S \in \mathcal{P}^m(S),\, S \in \mathcal{S}_h\} \end{align} or \begin{align} V_{{\rm dg},h}^m = \{v_h \in {L}^2 (\Omega) : {v_h}|_S \in \mathcal{P}^m(S),\, S \in \mathcal{S}_h\}, \end{align} with $m \ge 1$, are enriched accordingly to the discontinuity's position. Heaviside functions are used as enrichment, introducing $\mathcal{N}_{\rm enr} = \{1,\ldots,N_{\rm enr}\}$ strong discontinuous features to the approximation space $V_h^{\rm XFEM}$ and Nitsche's method \cite{nitsche} is used to impose interface and boundary conditions. Enriched functions $u_h \in V_h^{\rm XFEM}$ then have the representation \begin{equation} \begin{aligned} u_h = \sum_{j \in \mathcal{N}} u_j v_j + \sum_{i \in \mathcal{N}_{\rm enr}} \left(\sum_{k \in \mathcal{N}_i} u_{i,k} v_{i,k}\right) = {\boldsymbol u}_h \cdot {\boldsymbol v}_h \end{aligned} \end{equation} with \begin{equation} \begin{aligned} {\boldsymbol v}_h= [\underbrace{v_1,\ldots,v_{N_\mathcal{B}}}_{\rm std \ basis \ functions}, \underbrace{v_{1,1},\ldots,v_{1,|\mathcal{N}_1|}}_{{\rm add. \ basis \ functions \ for \ } \Gamma_{1,h}}, \ldots, \underbrace{v_{N_{\rm enr},1},\ldots,v_{N_{\rm enr},|\mathcal{N}_{N_{\rm enr}}|}}_{{\rm add. \ basis \ functions \ for \ } \Gamma_{N_{\rm enr},h} }]^T\\ \end{aligned} \end{equation} and \begin{equation} \begin{aligned} {\boldsymbol u}_h= [\underbrace{u_1,\ldots,u_{N_\mathcal{B}}}_{\rm std.\ coefficients},\underbrace{u_{1,1},\ldots,u_{1,|\mathcal{N}_1|}}_{{\rm coefficients\ for\ }\Gamma_{1,h}}, \ldots, \underbrace{u_{N_{\rm enr},1},\ldots,u_{N_{\rm enr},|\mathcal{N}_{N_{\rm enr}}|}}_{{\rm coefficients\ for\ }\Gamma_{N_{\rm enr},h}}]^T \end{aligned} \end{equation} whereas $\mathcal{N}_i$ denotes the respective set of enriched basis functions and $|\mathcal{N}_i|$ its cardinality. In regards to the linear system, the corresponding matrix has the same order which results in a block scheme \begin{equation} \begin{aligned} A = \begin{pmatrix} A_{{\rm STD} \times {\rm STD}} & A_{{\Gamma_1} \times {\rm STD}} & A_{{\Gamma_2} \times {\rm STD}} & \ldots& \ldots & A_{{\Gamma_{N_{\rm enr}}} \times {\rm STD}} \\ A_{{\rm STD} \times {\Gamma_1}} & A_{{\Gamma_1} \times {\Gamma_1}} & A_{{\Gamma_2} \times {\Gamma_1}} & \ldots &\ldots & A_{{\Gamma_{N_{\rm enr}}} \times {\Gamma_1}} \\ \vdots & \vdots & \vdots & \ldots &\ldots & \vdots\\ A_{{\rm STD} \times {\Gamma_{N_{\rm enr}}}} & A_{{\Gamma_1} \times {\Gamma_{N_{\rm enr}}}} & A_{{\Gamma_2} \times {\Gamma_{N_{\rm enr}}}} & \ldots&\ldots & A_{{\Gamma_{N_{\rm enr}}} \times {\Gamma_{N_{\rm enr}}}} \\ \end{pmatrix} \end{aligned} \end{equation}

with

\begin{equation} \begin{aligned} A_{{\rm STD} \times {\rm STD}} &= \begin{pmatrix} a_h(v_j,v_l) \end{pmatrix}_{j,l\in\mathcal{N}},\\ A_{{\Gamma_i} \times {\rm STD}} &= \begin{pmatrix} a_h(v_{i,k},v_j) \end{pmatrix}_{j\in\mathcal{N},\, k\in \mathcal{N}_i,\, i\in\mathcal{N}_{\rm enr}},\\ A_{{\rm STD}\times {\Gamma_i} } &= \begin{pmatrix} a_h(v_j,v_{i,k}) \end{pmatrix}_{j\in\mathcal{N},\, k\in \mathcal{N}_i,\, i\in\mathcal{N}_{\rm enr}},\\ A_{{\Gamma_i} \times{\Gamma_l}} &= \begin{pmatrix} a_h(v_{i,k},v_{l,\tilde{k}}) \end{pmatrix}_{k\in \mathcal{N}_i,\, \tilde{k}\in \mathcal{N}_l,\, i,l\in\mathcal{N}_{\rm enr}}.\\ \end{aligned} \end{equation}

The implementation of our toolbox \texttt{miXFEM} \cite{mixfem1,mixfem2}, partly based on the \texttt{PUM library} by Nikbaht and Wells \cite{pum}, consists of a level set toolbox \cite{levelset}, an extended \texttt{FEniCS Form Compiler}, and a \texttt{C++} implementation. It can handle an arbitrary (but a priori known) number of evolving discontinuities and has been validated by several steady-state and also time-dependent examples with known solutions from different fields and is currently used to tackle problems such as two-phase flows (Figure \ref{fig:two-phase}), Stefan problems, and laser-based keyhole welding (Figure \ref{fig:keyhole}).

\begin{figure} \centering \begin{minipage}[c]{0.245\textwidth} \centering \includegraphics[width=0.97\textwidth]{two-phase-flow.png} \caption{Two-phase flow} \label{fig:two-phase} \end{minipage} \begin{minipage}[c]{0.005\textwidth} \ \end{minipage} \begin{minipage}[c]{0.735\textwidth} \centering \includegraphics[width=0.99\textwidth]{keyhole.png} \caption{Laser-based keyhole welding} \label{fig:keyhole} \end{minipage} \end{figure}

\begin{thebibliography}{00} \addcontentsline{toc}{chapter}{References} % FIXNE: Correct order necessary! \bibitem{levelset} M. Jahn, T. Klock. A level set toolbox including reinitialization and mass correction algorithms for FEniCS. Technical Report 16-01, University of Bremen, 2016. \bibitem{mixfem1} M. Jahn, T. Klock. Numerical solution of the Stefan problem in level set formulation with the eXtended finite element method in FEniCS. Technical Report 17-01, University of Bremen, 2017. \bibitem{mixfem2} M. Jahn, A. Luttmann. Solving the Stefan problem with prescribed interface using an XFEM toolbox for FEniCS. Technical Report 16-03, University of Bremen, 2016. \bibitem{nitsche} J. Nitsche. {\"U}ber ein Variationsprinzip zur L{\"o}sung von Dirichlet-Problemen bei Verwendung von Teilr{\"a}umen, die keinen Randbedingungen unterworfen sind. Abhandlungen aus dem Mathematischen Seminar der Universit{\"a}t Hamburg, 36(1):9--15, 1971 \bibitem{pum} M. Nikbakht and G. N. Wells. Automated modeling of evolving discontinuities. Algorithms, 2(3):1008--1030, 2009. \end{thebibliography}

\end{document}

ABSTRACT. In this contribution we present an efficient strategy for numerical solution of viscoelastic flow problems the mathematical formulation of which involves the Oldroyd-B type constitutive models. The viscoelastic problem is discretized using the DEVSS method, originally developed by Guénette and Fortin (1995), based on the mixed finite element formulation involving three unknown fields (velocity, pressure and viscous polymer stress). The arising algebraic systems of equations are tackled using iterative solvers of the Krylov subspace type with a certain type of block preconditioning strategy following the ideas from Elman et al. (2014) and Hwang (2011).

Motivated by the large scale problems encountered in rubber industry, we aim to provide the implementation suitable for deployment on HPC systems. We build on top of the FEniCS Project (fenicsproject.org), revisiting the implementation of the DEVSS approach in Rheagen (launchpad.net/rheagen), and we implement the preconditioning strategy using the state-of-the-art algorithms available in PETSc (www.mcs.anl.gov/petsc) by extending the ideas previously applied in FENaPack (github.com/blechta/fenapack).

ABSTRACT. We consider a multiscale method to solve parabolic equations with a highly varying spatial coefficients. These equations can be used to model physical behavior in a heterogeneous powder medium of the selective laser melting process. In this study a localized orthogonal decomposition (LOD) technique is used to construct prolongation and restriction operators for generating coarse scale system with the aid of FEniCS. LOD is used to solve the heat equation with a highly varying thermal conductivity. Based on numerical examples, we show that optimal convergence rate of the method is independent of the variations in the coefficient.