Search

We address the numerical challenge of solving the Hookean-type time-fractional Navier--Stokes--Fokker--Planck equation, a history-dependent system of PDEs defined on the Cartesian product of two $d$-dimensional spaces in the turbulent regime. Due to its high dimensionality, the non-locality with respect to time, and the resolution required to resolve turbulent flow, this problem is highly demanding. To overcome these challenges, we employ the Hermite spectral method for the configuration space of the Fokker--Planck equation, reducing the problem to a purely macroscopic model. Considering scenarios for available analytical solutions, we prove the existence of an optimal choice of the Hermite scaling parameter. With this choice, the macroscopic system is equivalent to solving the coupled micro-macro system. We apply second-order time integration and extrapolation of the coupling terms, achieving, for the first time, convergence rates for the fully coupled time-fractional system independent of the order of the time-fractional derivative. Our efficient implementation of the numerical scheme allows turbulent simulations of dilute polymeric fluids with memory effects in two and three dimensions. Numerical simulations show that memory effects weaken the drag-reducing effect of added polymer molecules in the turbulent flow regime.

Consensus-based optimization (CBO) is a versatile multi-particle optimization method for performing nonconvex and nonsmooth global optimizations in high dimensions. Proofs of global convergence in probability have been achieved for a broad class of objective functions in unconstrained optimizations. In this work we adapt the algorithm for solving constrained optimizations on compact and unbounded domains with boundary by leveraging emerging reflective boundary conditions. In particular, we close a relevant gap in the literature by providing a global convergence proof for the many-particle regime comprehensive of convergence rates. On the one hand, for the sake of minimizing running cost, it is desirable to keep the number of particles small. On the other hand, reducing the number of particles implies a diminished capability of exploration of the algorithm. Hence numerical heuristics are needed to ensure convergence of CBO in the few-particle regime. In this work, we also significantly improve the convergence and complexity of CBO by utilizing an adaptive region control mechanism and by choosing geometry-specific random noise. In particular, by combining a hierarchical noise structure with a multigrid finite element method, we are able to compute global minimizers for a constrained $p$-Allen-Cahn problem with obstacles, a very challenging variational problem.

Previous work shows evidence that cerebrospinal fluid (CSF) plays a crucial role in brain waste clearance processes, and that altered flow patterns are associated with various diseases of the central nervous system. In this study, we investigate the potential of deep learning to predict the distribution in human brain of a gadolinium-based CSF contrast agent (tracer) administered intrathecal. For this, T1-weighted magnetic resonance imaging (MRI) scans taken at multiple time points before and after intrathecal injection were utilized. We propose a U-net-based supervised learning model to predict pixel-wise signal increases at their peak after 24 hours. Its performance is evaluated based on different tracer distribution stages provided during training, including predictions from baseline scans taken before injection. Our findings indicate that using imaging data from just the first two hours post-injection for training yields tracer flow predictions comparable to those trained with additional later-stage scans. The model was further validated by comparing ventricular reflux gradings provided by neuroradiologists, and inter-rater grading among medical experts and the model showed excellent agreement. Our results demonstrate the potential of deep learning-based methods for CSF flow prediction, suggesting that fewer MRI scans could be sufficient for clinical analysis, which might significantly improve clinical efficiency, patient well-being, and lower healthcare costs., Comment: 13 pages, 7 figures

Biot's consolidation model is a classical model for the evolution of deformable porous media saturated by a fluid and has various interdisciplinary applications. While numerical solution methods to solve poroelasticity by typical schemes such as finite differences, finite volumes or finite elements have been intensely studied, lattice Boltzmann methods for poroelasticity have not been developed yet. In this work, we propose a novel semi-implicit coupling of lattice Boltzmann methods to solve Biot's consolidation model in two dimensions. To this end, we use a single-relaxation-time lattice Boltzmann method for reaction-diffusion equations to solve the Darcy flow and combine it with a recent pseudo-time multi-relaxation-time lattice Boltzmann scheme for quasi-static linear elasticity by Boolakee, Geier and De Lorenzis (2023, DOI: 10.1016/j.cma.2022.115756). The numerical results demonstrate that naive coupling schemes lead to instabilities when the poroelastic system is strongly coupled. However, the newly developed centered coupling scheme using fully explicit and semi-implicit contributions is stable and accurate in all considered cases, even for the Biot--Willis coefficient being one. Furthermore, the numerical results for Terzaghi's consolidation problem and a two-dimensional extension thereof highlight that the scheme is even able to capture discontinuous solutions arising from instantaneous loading., Comment: 21 pages, 11 figures

A velocity field characterized by chaotic advection induces stretching and folding processes that increase the solute-solvent interface available for diffusion. Identifying chaotic flow fields with optimized mixing enhancement is relevant for applications like groundwater remediation, microfluidics, and many more. This work uses the dilution index to quantify the temporal increase in mixing of a solute within its solvent. We introduce a new approach to select a suitable grid size for each time step in the dilution index approximation, motivated by the theory of representative elementary volumes. It preserves the central feature of the dilution index, which is monotonically increasing in time and hence leads to reliable results. Our analysis highlights the importance of a suitable choice for the grid size in the dilution index approximation. We use this approach to demonstrate the mixing enhancement for two chaotic injection-extraction systems that exhibit chaotic structures: a source-sink dipole and a rotated potential mixing. By analyzing the chaotic flow fields, we identify Kolmogorov--Arnold--Moser (KAM) islands, non-mixing regions that limit the chaotic area in the domain and, thereby, the mixing enhancement. Using our new approach, we assess the choice of design parameters of the injection-extraction systems to effectively engineer chaotic mixing. We demonstrate the important role of diffusion in filling the KAM islands and reaching complete mixing in the systems., Comment: paper: 17 pages, 15 figures. supplemental material: 4 pages, 3 figures, submitted to Physical Review E

In this work, we investigate a quasilinear subdiffusion model which involves a fractional derivative of order $\alpha \in (0,1)$ in time and a nonlinear diffusion coefficient. First, using smoothing properties of solution operators for linear subdiffusion and a perturbation argument, we prove several pointwise-in-time regularity estimates that are useful for numerical analysis. Then we develop a high-order time stepping scheme for solving quasilinear subdiffusion, based on convolution quadrature generated by second-order backward differentiation formula with correction at the first step. Further, we establish that the convergence order of the scheme is $O(\tau^{1+\alpha-\epsilon})$ without imposing any additional assumption on the regularity of the solution. The analysis relies on refined Sobolev regularity of the nonlinear perturbation remainder and smoothing properties of discrete solution operators. Several numerical experiments in two space dimensions show the sharpness of the error estimate., Comment: 21 pages, 1 figure

Understanding how genetically encoded rules drive and guide complex neuronal growth processes is essential to comprehending the brain's architecture, and agent-based models (ABMs) offer a powerful simulation approach to further develop this understanding. However, accurately calibrating these models remains a challenge. Here, we present a novel application of Approximate Bayesian Computation (ABC) to address this issue. ABMs are based on parametrized stochastic rules that describe the time evolution of small components -- the so-called agents -- discretizing the system, leading to stochastic simulations that require appropriate treatment. Mathematically, the calibration defines a stochastic inverse problem. We propose to address it in a Bayesian setting using ABC. We facilitate the repeated comparison between data and simulations by quantifying the morphological information of single neurons with so-called morphometrics and resort to statistical distances to measure discrepancies between populations thereof. We conduct experiments on synthetic as well as experimental data. We find that ABC utilizing Sequential Monte Carlo sampling and the Wasserstein distance finds accurate posterior parameter distributions for representative ABMs. We further demonstrate that these ABMs capture specific features of pyramidal cells of the hippocampus (CA1). Overall, this work establishes a robust framework for calibrating agent-based neuronal growth models and opens the door for future investigations using Bayesian techniques for model building, verification, and adequacy assessment., Comment: 32 pages, 12 Figures

Rupture of intracranial aneurysms results in severe subarachnoidal hemorrhage, which is associated with high morbidity and mortality. Neurointerventional occlusion of the aneurysm through coiling has evolved to a therapeutical standard. The choice of the specific coil has an important influence on secondary regrowth requiring retreatment. Aneurysm occlusion was simulated either through virtual implantation of a preshaped 3D coil or with a porous media approach. In this study, we used a recently developed numerical approach to simulate aneurysm shapes in specific challenging aneurysm anatomies and correlated these with aneurysm recurrence 6 months after treatment. The simulation showed a great variety of coil shapes depending on the variability in possible microcatheter positions. Aneurysms with a later recurrence showed a tendency for more successful coiling attempts. Results revealed further trends suggesting lower simulated packing densities in aneurysms with reoccurrence. Simulated packing densities did not correlate with those calculated by conventional software, indicating the potential for our approach to offer additional predictive value. Our study, therefore, pioneers a comprehensive numerical model for simulating aneurysm coiling, providing insights into individualized treatment strategies and outcome prediction. Future directions involve expanding the model's capabilities to simulate intraprocedural outcomes and long-term predictions, aiming to refine occlusion quality criteria and validate prediction parameters in larger patient cohorts. This simulation framework holds promise for enhancing clinical decision-making and optimizing patient outcomes in endovascular aneurysm treatment., Comment: 13 pages, 4 figures

This work highlights an approach for incorporating realistic uncertainties into scientific computing workflows based on finite elements, focusing on applications in computational mechanics and design optimization. We leverage Mat\'ern-type Gaussian random fields (GRFs) generated using the SPDE method to model aleatoric uncertainties, including environmental influences, variating material properties, and geometric ambiguities. Our focus lies on delivering practical GRF realizations that accurately capture imperfections and variations and understanding how they impact the predictions of computational models and the topology of optimized designs. We describe a numerical algorithm based on solving a generalized SPDE to sample GRFs on arbitrary meshed domains. The algorithm leverages established techniques and integrates seamlessly with the open-source finite element library MFEM and associated scientific computing workflows, like those found in industrial and national laboratory settings. Our solver scales efficiently for large-scale problems and supports various domain types, including surfaces and embedded manifolds. We showcase its versatility through biomechanics and topology optimization applications. The flexibility and efficiency of SPDE-based GRF generation empower us to run large-scale optimization problems on 2D and 3D domains, including finding optimized designs on embedded surfaces, and to generate topologies beyond the reach of conventional techniques. Moreover, these capabilities allow us to model geometric uncertainties of reconstructed submanifolds, such as the surfaces of cerebral aneurysms. In addition to offering benefits in these specific domains, the proposed techniques transcend specific applications and generalize to arbitrary forward and backward problems in uncertainty quantification involving finite elements., Comment: 38 pages, 21 figures

Intracranial aneurysms are the leading cause of hemorrhagic stroke. One of the established treatment approaches is the embolization induced by coil insertion. However, the prediction of treatment and subsequent changed flow characteristics in the aneurysm is still an open problem. In this work, we present an approach based on a patient-specific geometry and parameters including a coil representation as inhomogeneous porous medium. The model consists of the volume-averaged Navier-Stokes equations for a non-Newtonian blood rheology. We solve these equations using a problem-adapted lattice Boltzmann method and present a comparison between fully-resolved and volume-averaged simulations. The results indicate the validity of the model. Overall, this workflow allows for patient specific assessment of the flow due to potential treatment., Comment: 10 pages, 9 figures

Predicting the long-term success of endovascular interventions in the clinical management of cerebral aneurysms requires detailed insight into the patient-specific physiological conditions. In this work, we not only propose numerical representations of endovascular medical devices such as coils, flow diverters or Woven EndoBridge but also outline numerical models for the prediction of blood flow patterns in the aneurysm cavity right after a surgical intervention. Detailed knowledge about the post-surgical state then lays the basis to assess the chances of a stable occlusion of the aneurysm required for a long-term treatment success. To this end, we propose mathematical and mechanical models of endovascular medical devices made out of thin metal wires. These can then be used for fully resolved flow simulations of the post-surgical blood flow, which in this work will be performed by means of a Lattice Boltzmann method applied to the incompressible Navier-Stokes equations and patient-specific geometries. To probe the suitability of homogenized models, we also investigate poro-elastic models to represent such medical devices. In particular, we examine the validity of this modeling approach for flow diverter placement across the opening of the aneurysm cavity. For both approaches, physiologically meaningful boundary conditions are provided from reduced-order models of the vascular system. The present study demonstrates our capabilities to predict the post-surgical state and lays a solid foundation to tackle the prediction of thrombus formation and, thus, the aneurysm occlusion in a next step.

The flow of cerebrospinal fluid through the perivascular spaces of the brain is believed to play a crucial role in eliminating toxic waste proteins. While the driving forces of this flow have been enigmatic, experiments have shown that arterial wall motion is central. In this work, we present a network model for simulating pulsatile fluid flow in perivascular networks. We establish the well-posedness of this model in the primal and dual mixed variational settings, and show how it can be discretized using mixed finite elements. Further, we utilize this model to investigate fundamental questions concerning the physical mechanisms governing perivascular fluid flow. Notably, our findings reveal that arterial pulsations can induce directional flow in branching perivascular networks.

We consider a Cauchy problem for the inhomogeneous differential equation given in terms of an unbounded linear operator $A$ and the Caputo fractional derivative of order $\alpha \in (0, 2)$ in time. The previously known representation of the mild solution to such a problem does not have a conventional variation-of-constants like form, with the propagator derived from the associated homogeneous problem. Instead, it relies on the existence of two propagators with different analytical properties. This fact limits theoretical and especially numerical applicability of the existing solution representation. Here, we propose an alternative representation of the mild solution to the given problem that consolidates the solution formulas for sub-parabolic, parabolic and sub-hyperbolic equations with a positive sectorial operator $A$ and non-zero initial data. The new representation is solely based on the propagator of the homogeneous problem and, therefore, can be considered as a more natural fractional extension of the solution to the classical parabolic Cauchy problem. By exploiting a trade-off between the regularity assumptions on the initial data in terms of the fractional powers of $A$ and the regularity assumptions on the right-hand side in time, we show that the proposed solution formula is strongly convergent for $t \geq 0$ under considerably weaker assumptions compared to the standard results from the literature. Crucially, the achieved relaxation of space regularity assumptions ensures that the new solution representation is practically feasible for any $\alpha \in (0, 2)$ and is amenable to the numerical evaluation using uniformly accurate quadrature-based algorithms., Comment: 20 pages

We investigate the well-posedness of a coupled Navier-Stokes-Fokker-Planck system with a time-fractional derivative. Such systems arise in the kinetic theory of dilute solutions of polymeric liquids, where the motion of noninteracting polymer chains in a Newtonian solvent is modelled by a stochastic process exhibiting power-law waiting time, in order to capture subdiffusive processes associated with non-Fickian diffusion. We outline the derivation of the model from a subordinated Langevin equation. The elastic properties of the polymer molecules immersed in the solvent are modelled by a finitely extensible nonlinear elastic (FENE) dumbbell model, and the drag term in the Fokker--Planck equation is assumed to be corotational. We prove the global-in-time existence of large-data weak solutions to this time-fractional model of order $\alpha \in (\tfrac12,1)$, and derive an energy inequality satisfied by weak solutions.

Cancer is a disease driven by random DNA mutations and the interaction of many complex phenomena. To improve the understanding and ultimately find more effective treatments, researchers leverage computer simulations mimicking the tumor growth in silico. The challenge here is to account for the many phenomena influencing the disease progression and treatment protocols. This work introduces a computational model to simulate vascular tumor growth and the response to drug treatments in 3D. It consists of two agent-based models for the tumor cells and the vasculature. Moreover, partial differential equations govern the diffusive dynamics of the nutrients, the vascular endothelial growth factor, and two cancer drugs. The model focuses explicitly on breast cancer cells over-expressing HER2 receptors and a treatment combining standard chemotherapy (Doxorubicin) and monoclonal antibodies with anti-angiogenic properties (Trastuzumab). However, large parts of the model generalize to other scenarios. We show that the model qualitatively captures the effects of the combination therapy by comparing our simulation results with previously published pre-clinical data. Furthermore, we demonstrate the scalability of the model and the associated C++ code by simulating a vascular tumor occupying a volume of 400mm3 using a total of 92.5 million agents.

We present an exponentially convergent numerical method to approximate the solution of the Cauchy problem for the inhomogeneous fractional differential equation with an unbounded operator coefficient and Caputo fractional derivative in time. The numerical method is based on the newly obtained solution formula that consolidates the mild solution representations of sub-parabolic, parabolic and sub-hyperbolic equations with sectorial operator coefficient $A$ and non-zero initial data. The involved integral operators are approximated using the sinc-quadrature formulas that are tailored to the spectral parameters of $A$, fractional order $\alpha$ and the smoothness of the first initial condition, as well as to the properties of the equation's right-hand side $f(t)$. The resulting method possesses exponential convergence for positive sectorial $A$, any finite $t$, including $t = 0$ and the whole range $\alpha \in (0,2)$. It is suitable for a practically important case, when no knowledge of $f(t)$ is available outside the considered interval $t \in [0, T]$. The algorithm of the method is capable of multi-level parallelism. We provide numerical examples that confirm the theoretical error estimates., Comment: This version supersedes the official publication (https://www.mdpi.com/2227-7390/11/10/2312) at present time. Several typos were corrected and corrigendum added

Formulating tumor models that predict growth under therapy is vital for improving patient-specific treatment plans. In this context, we present our recent work on simulating non-small-scale cell lung cancer (NSCLC) in a simple, deterministic setting for two different patients receiving an immunotherapeutic treatment. At its core, our model consists of a Cahn-Hilliard-based phase-field model describing the evolution of proliferative and necrotic tumor cells. These are coupled to a simplified nutrient model that drives the growth of the proliferative cells and their decay into necrotic cells. The applied immunotherapy decreases the proliferative cell concentration. Here, we model the immunotherapeutic agent concentration in the entire lung over time by an ordinary differential equation (ODE). Finally, reaction terms provide a coupling between all these equations. By assuming spherical, symmetric tumor growth and constant nutrient inflow, we simplify this full 3D cancer simulation model to a reduced 1D model. We can then resort to patient data gathered from computed tomography (CT) scans over several years to calibrate our model. For the reduced 1D model, we show that our model can qualitatively describe observations during immunotherapy by fitting our model parameters to existing patient data. Our model covers cases in which the immunotherapy is successful and limits the tumor size, as well as cases predicting a sudden relapse, leading to exponential tumor growth. Finally, we move from the reduced model back to the full 3D cancer simulation in the lung tissue. Thereby, we show the predictive benefits a more detailed patient-specific simulation including spatial information could yield in the future.

Matrix-free techniques play an increasingly important role in large-scale simulations. Schur complement techniques and massively parallel multigrid solvers for second-order elliptic partial differential equations can significantly benefit from reduced memory traffic and consumption. The matrix-free approach often restricts solver components to purely local operations, for instance, the Jacobi- or Gauss--Seidel-Smoothers in multigrid methods. An incomplete LU (ILU) decomposition cannot be calculated from local information and is therefore not amenable to an on-the-fly computation which is typically needed for matrix-free calculations. It generally requires the storage and factorization of a sparse matrix which contradicts the low memory requirements in large scale scenarios. In this work, we propose a matrix-free ILU realization. More precisely, we introduce a memory-efficient, matrix-free ILU(0)-Smoother component for low-order conforming finite elements on tetrahedral hybrid grids. Hybrid grids consist of an unstructured macro-mesh which is subdivided into a structured micro-mesh. The ILU(0) is used for degrees-of-freedom assigned to the interior of macro-tetrahedra. This ILU(0)-Smoother can be used for the efficient matrix-free evaluation of the Steklov-Poincare operator from domain-decomposition methods. After introducing and formally defining our smoother, we investigate its performance on refined macro-tetrahedra. Secondly, the ILU(0)-Smoother on the macro-tetrahedrons is implemented via surrogate matrix polynomials in conjunction with a fast on-the-fly evaluation scheme resulting in an efficient matrix-free algorithm. The polynomial coefficients are obtained by solving a least-squares problem on a small part of the factorized ILU(0) matrices to stay memory efficient. The convergence rates of this smoother with respect to the polynomial order are thoroughly studied.

Reducing the intensity of wind excitation via aerodynamic shape modification is a major strategy to mitigate the reaction forces on supertall buildings, reduce construction and maintenance costs, and improve the comfort of future occupants. To this end, computational fluid dynamics (CFD) combined with state-of-the-art stochastic optimization algorithms is more promising than the trial and error approach adopted by the industry. The present study proposes and investigates a novel approach to risk-averse shape optimization of tall building structures that incorporates site-specific uncertainties in the wind velocity, terrain conditions, and wind flow direction. A body-fitted finite element approximation is used for the CFD with different wind directions incorporated by re-meshing the fluid domain. The bending moment at the base of the building is minimized, resulting in a building with reduced cost, material, and hence, a reduced carbon footprint. Both risk-neutral and risk-averse optimization of the twist and tapering of a representative building are presented under uncertain inflow wind conditions that have been calibrated to fit freely-available site-specific data from Basel, Switzerland. The risk-averse strategy uses the conditional value-at-risk to optimize for the low-probability high-consequence events appearing in the worst 10% of loading conditions. Adaptive sampling is used to accelerate the gradient-based stochastic optimization pipeline. The adaptive method is easy to implement and particularly helpful for compute-intensive simulations because the number of gradient samples grows only as the optimal design algorithm converges. The performance of the final risk-averse building geometry is exceptionally favorable when compared to the risk-neutral optimized geometry, thus, demonstrating the effectiveness of the risk-averse design approach in computational wind engineering.

In this work, we present mixed dimensional models for simulating blood flow and transport processes in breast tissue and the vascular tree supplying it. These processes are considered, to start from the aortic inlet to the capillaries and tissue of the breast. Large variations in biophysical properties and flow conditions exist in this system necessitating the use of different flow models for different geometries and flow regimes. Large variations in biophysical properties and flow conditions exist in this system necessitating the use of different flow models for different geometries and flow regimes. In total, we consider four different model types. First, a system of 1D nonlinear hyperbolic PDEs is considered to simulate blood flow in larger arteries with highly elastic vessel walls. Second, we assign 1D linearized hyperbolic PDEs to model the smaller arteries with stiffer vessel walls. The third model type consists of ODE systems (0D models). It is used to model the arterioles and peripheral circulation. Finally, homogenized 3D porous media models are considered to simulate flow and transport in capillaries and tissue within the breast volume. Sink terms are used to account for the influence of the venous and lymphatic systems. Combining the four model types, we obtain two different 1D-0D-3D coupled models for simulating blood flow and transport processes: The first 1D-0D-3D model covers the whole path from the aorta to the breast, while the second model is a sub-model obtained by restriction to breast vasculature and tissue making possible a significant reduction in computational cost. Several numerical experiments are conducted that demonstrate realistic flow simulations compared to existing data on blood flow in human breast and vascular system.

In this work, we investigate the inverse problem of determining the kernel functions that best describe the mechanical behavior of a complex medium modeled by a general nonlocal viscoelastic wave equation. To this end, we minimize a tracking-type data misfit function under this PDE constraint. We perform the well-posedness analysis of the state and adjoint problems and, using these results, rigorously derive the first-order sensitivities. Numerical experiments in a three-dimensional setting illustrate the method.

Manufactured materials usually contain random imperfections due to the fabrication process, e.g., the 3D-printing, casting, etc. These imperfections affect significantly the effective material properties and result in uncertainties in the mechanical response. Numerical analysis of the effects of the imperfections and the uncertainty quantification (UQ) can be often done by use of digital stochastic surrogate material models. In this work, we present a new flexible class of surrogate models depending on a small number of parameters with special focus on two-phase materials. The surrogate models are constructed as the level-set of a linear combination of an intensity field representing the topological shape and a Gaussian perturbation representing the imperfections. The mathematical design parameters of the model are related to physical ones and thus easy to interpret. The calibration of the model parameters is performed using progressive batching sub-sampled quasi-Newton minimization, using a designed distance measure between the synthetic samples and the data. Then, employing a fast sampling algorithm, an arbitrary number of synthetic samples can be generated to use in Monte Carlo type methods. In particular, we illustrate the method in application to UQ of the elasto-plastic response of an imperfect octet-truss lattice which plays an important role in additive manufacturing. To this end, we study the effective material properties of the lattice unit cell under elasto-plastic deformations and investigate the sensitivity of the effective Young's modulus to the imperfections., Comment: 21 pages, 15 figures

As a mean to assess the risk dam structures are exposed to during earthquakes, we employ an abstract mathematical, three dimensional, elasto-acoustic coupled wave-propagation model taking into account (i) the dam structure itself, embedded into (ii) its surrounding topography, (iii) different material soil layers, (iv) the seismic source as well as (v) the reservoir lake filled with water treated as an acoustic medium. As a case study for extensive numerical simulations we consider the magnitude 7 seismic event of the 30$^{\rm th}$ of October 2020 taking place in the Icarian Sea (Greece) and the Tahtali dam around 30 km from there (Turkey). A challenging task is to resolve the multiple length scales that are present due to the huge differences in size between the dam building structure and the area of interest, considered for the propagation of the earthquake. Interfaces between structures and highly non-conforming meshes on different scales are resolved by means of a discontinuous Galerkin approach. The seismic source is modeled using inversion data about the real fault plane. Ultimately, we perform a real data driven, multi-scale, full source-to-site, physics based simulation based on the discontinuous Galerkin spectral element method, which allows to precisely validate the ground motion experienced along the Tahtali dam, comparing the synthetic seismograms against actually observed ones. A comparison with a more classical computational method, using a plane wave with data from a deconvolved seismogram reading as an input, is discussed., Comment: 26 pages, 16 figures

Subsurface flows, particularly hyporheic exchange fluxes, driven by streambed topography, permeability, channel gradient and dynamic flow conditions provide prominent ecological services such as nitrate removal from streams and aquifers. Stream flow dynamics cause strongly nonlinear and often episodic contributions of nutrient concentrations in river-aquifer systems. Using a fully coupled transient flow and reactive transport model, we investigated the denitrification potential of hyporheic zones during peak-flow events. The effects of streambed permeability, channel gradient and bedform amplitude on the spatio-temporal distribution of nitrate and dissolved organic carbon in streambeds and the associated denitrification potential were explored. Distinct peak-flow events with different intensity, duration and hydrograph shape were selected to represent a wide range of peak-flow scenarios. Our results indicated that the specific hydrodynamic characteristics of individual flow events largely determine the average positive or negative nitrate removal capacity of hyporheic zones, however the magnitude of this capacity is controlled by geomorphological settings (i.e. channel slope, streambed permeability and bedform amplitude). Specifically, events with longer duration and higher intensity were shown to promote higher nitrate removal efficiency with higher magnitude of removal efficiency in the scenarios with higher slope and permeability values. These results are essential for better assessment of the subsurface nitrate removal capacity under the influence of flow dynamics and particularly peak-flow events in order to provide tailored solutions for effective restoration of interconnected river-aquifer systems.

Alpine permafrost environments are highly vulnerable and sensitive to changes in regional and global climate trends. Thawing and degradation of permafrost has numerous adverse environmental, economic, and societal impacts. Mathematical modeling and numerical simulations provide powerful tools for predicting the degree of degradation and evolution of subsurface permafrost as a result of global warming. A particularly significant characteristic of alpine environments is the high variability in their topography and geomorphology which drives large lateral thermal and fluid fluxes. Additionally, harsh winds, extreme weather conditions, and various degrees of saturation have to be considered. The combination of large lateral fluxes and unsaturated ground makes alpine systems markedly different from Arctic permafrost environments and general geotechnical ground freezing applications, and therefore, alpine permafrost demands its own specialized modeling approaches. In this research work, we present a multi-physics permafrost model tailored to alpine regions. In particular, we resolve the ice-water phase transitions, unsaturated conditions, and capillary actions, and account for the impact of the evolving pore volume on fluid-matrix interactions. Moreover, the approach is multi-dimensional, and therefore, inherently resolves fluxes along topographic gradients. Through numerical cases studies based on the elevation profiles of the two prominent peaks of the Zugspitze (DE) and the Matterhorn (CH), we show the strong influence of topography driven thermal and fluid fluxes on active layer dynamics and the distribution of permafrost., Comment: 20 pages (without appendix), 9 figures

We develop a novel data-driven approach to modeling the atmospheric boundary layer. This approach leads to a nonlocal, anisotropic synthetic turbulence model which we refer to as the deep rapid distortion (DRD) model. Our approach relies on an operator regression problem which characterizes the best fitting candidate in a general family of nonlocal covariance kernels parameterized in part by a neural network. This family of covariance kernels is expressed in Fourier space and is obtained from approximate solutions to the Navier--Stokes equations at very high Reynolds numbers. Each member of the family incorporates important physical properties such as mass conservation and a realistic energy cascade. The DRD model can be calibrated with noisy data from field experiments. After calibration, the model can be used to generate synthetic turbulent velocity fields. To this end, we provide a new numerical method based on domain decomposition which delivers scalable, memory-efficient turbulence generation with the DRD model as well as others. We demonstrate the robustness of our approach with both filtered and noisy data coming from the 1968 Air Force Cambridge Research Laboratory Kansas experiments. Using this data, we witness exceptional accuracy with the DRD model, especially when compared to the International Electrotechnical Commission standard.

Time-fractional partial differential equations are nonlocal in time and show an innate memory effect. In this work, we propose an augmented energy functional which includes the history of the solution. Further, we prove the equivalence of a time-fractional gradient flow problem to an integer-order one based on our new energy. This equivalence guarantees the dissipating character of the augmented energy. The state function of the integer-order gradient flow acts on an extended domain similar to the Caffarelli-Silvestre extension for the fractional Laplacian. Additionally, we apply a numerical scheme for solving time-fractional gradient flows, which is based on kernel compressing methods. We illustrate the behavior of the original and augmented energy in the case of the Ginzburg-Landau energy functional.

In this paper, we derive the time-fractional Cahn-Hilliard equation from continuum mixture theory with a modification of Fick's law of diffusion. This model describes the process of phase separation with nonlocal memory effects. We analyze the existence, uniqueness, and regularity of weak solutions of the time-fractional Cahn-Hilliard equation. In this regard, we consider degenerating mobility functions and free energies of Landau, Flory--Huggins and double-obstacle type. We apply the Faedo-Galerkin method to the system, derive energy estimates, and use compactness theorems to pass to the limit in the discrete form. In order to compensate for the missing chain rule of fractional derivatives, we prove a fractional chain inequality for semiconvex functions. The work concludes with numerical simulations and a sensitivity analysis showing the influence of the fractional power. Here, we consider a convolution quadrature scheme for the time-fractional component, and use a mixed finite element method for the space discretization.

Due to its significance in terms of wave phenomena a considerable effort has been put into the design of preconditioners for the Helmholtz equation. One option to derive a preconditioner is to apply a multigrid method on a shifted operator. In such an approach, the wavenumber is shifted by some imaginary value. This step is motivated by the observation that the shifted problem can be more efficiently handled by iterative solvers when compared to the standard Helmholtz equation. However, up to now, it is not obvious what the best strategy for the choice of the shift parameter is. It is well known that a good shift parameter depends sensitively on the wavenumber and the discretization parameters such as the order and the mesh size. Therefore, we study the choice of a near optimal complex shift such that an FGMRES solver converges with fewer iterations. Our goal is to provide a map which returns the near optimal shift for the preconditioner depending on the wavenumber and the mesh size. In order to compute this map, a data driven approach is considered: We first generate many samples, and in a second step, we perform a nonlinear regression on this data. With this representative map, the near optimal shift can be obtained by a simple evaluation. Our preconditioner is based on a twogrid V-cycle applied to the shifted problem, allowing us to implement a semi matrix-free method. The performance of our preconditioned FGMRES solver is illustrated by several benchmark problems with heterogeneous wavenumbers in two and three space dimensions.

Fractional differential equations (FDEs) describe subdiffusion behavior of dynamical systems. Its non-local structure requires taking into account the whole evolution history during the time integration, which then possibly causes additional memory use to store the history, growing in time. An alternative to a quadrature for the history integral is to approximate the fractional kernel with the sum of exponentials, which is equivalent to considering the FDE solution as a sum of solutions to a system of ODEs. One possibility to construct this system is to approximate the Laplace spectrum of the fractional kernel with a rational function. In this paper, we use the adaptive Antoulas--Anderson (AAA) algorithm for the rational approximation of the kernel spectrum which yields only a small number of real valued poles. We propose a numerical scheme based on this idea and study its stability and convergence properties. In addition, we apply the algorithm to a time-fractional Cahn-Hilliard problem., Comment: 26 pages, 10 figures

Inspired by medical applications of high-intensity ultrasound, we study a coupled elasto-acoustic problem with general acoustic nonlinearities of quadratic type as they arise, for example, in the Westervelt and Kuznetsov equations of nonlinear acoustics. We derive convergence rates in the energy norm of a finite element approximation to the coupled problem in a setting that involves different acoustic materials and hence jumps within material parameters. A subdomain-based discontinuous Galerkin approach realizes the acoustic-acoustic coupling of different materials. At the same time, elasto-acoustic interface conditions are used for a mutual exchange of forces between the different models. Numerical simulations back up the theoretical findings in a three-dimensional setting with academic test cases as well as in an application-oriented simulation, where the modeling of human tissue as an elastic versus an acoustic medium is compared., Comment: 41 pages, 10 figures

In this work, we present a coupled 3D-1D model of solid tumor growth within a dynamically changing vascular network to facilitate realistic simulations of angiogenesis. Additionally, the model includes erosion of the extracellular matrix, interstitial flow, and coupled flow in blood vessels and tissue. We employ continuum mixture theory with stochastic Cahn--Hilliard type phase-field models of tumor growth. The interstitial flow is governed by a mesoscale version of Darcy's law. The flow in the blood vessels is controlled by Poiseuille flow, and Starling's law is applied to model the mass transfer in and out of blood vessels. The evolution of the network of blood vessels is orchestrated by the concentration of the tumor angiogenesis factors (TAFs); blood vessels grow towards the increasing TAFs concentrations. This process is not deterministic, allowing random growth of blood vessels and, therefore, due to the coupling of nutrients in tissue and vessels, makes the growth of tumors stochastic. We demonstrate the performance of the model by applying it to a variety of scenarios. Numerical experiments illustrate the flexibility of the model and its ability to generate satellite tumors. Simulations of the effects of angiogenesis on tumor growth are presented as well as sample-independent features of cancer., Comment: 36 pages, 22 figures

A novel mathematical model for fiber-reinforced materials is proposed. It is based on a 1-dimensional beam model for the thin fiber structures, a flexible and general 3-dimensional elasticity model for the matrix and an overlapping domain decomposition approach. From a computational point of view, this is motivated by the fact that matrix and fibers can easily meshed independently. Our main interest is in fiber reinforce polymers where the Young's modulus are quite different. Thus the modeling error from the overlapping approach is of no significance. The coupling conditions acknowledge both, the forces and the moments of the beam model and transfer them to the background material. A suitable static condensation procedure is applied to remove the beam balance equations. The condensed system then forms our starting point for a numerical approximation in terms of isogeometric analysis. The choice of our discrete basis functions of higher regularity is motivated by the fact, that as a result of the static condensation, we obtain second gradient terms in fiber direction. Eventually, a series of benchmark tests demonstrate the flexibility and robustness of the proposed methodology. As a proof-of-concept, we show that our new model is able to capture bending, torsion and shear dominated situations.

We introduce adaptive sampling methods for stochastic programs with deterministic constraints. First, we propose and analyze a variant of the stochastic projected gradient method where the sample size used to approximate the reduced gradient is determined on-the-fly and updated adaptively. This method is applicable to a broad class of expectation-based risk measures and leads to a significant reduction in the individual gradient evaluations used to estimate the objective function gradient. Numerical experiments with expected risk minimization and conditional value-at-risk minimization support this conclusion and demonstrate practical performance and efficacy for both risk-neutral and risk-averse problems. Second, we propose an SQP-type method based on similar adaptive sampling principles. The benefits of this method are demonstrated in a simplified engineering design application featuring risk-averse shape optimization of a steel shell structure subject to uncertain loading conditions and model uncertainty.

This work is based on the seminar titled ``Resiliency in Numerical Algorithm Design for Extreme Scale Simulations'' held March 1-6, 2020 at Schloss Dagstuhl, that was attended by all the authors. Naive versions of conventional resilience techniques will not scale to the exascale regime: with a main memory footprint of tens of Petabytes, synchronously writing checkpoint data all the way to background storage at frequent intervals will create intolerable overheads in runtime and energy consumption. Forecasts show that the mean time between failures could be lower than the time to recover from such a checkpoint, so that large calculations at scale might not make any progress if robust alternatives are not investigated. More advanced resilience techniques must be devised. The key may lie in exploiting both advanced system features as well as specific application knowledge. Research will face two essential questions: (1) what are the reliability requirements for a particular computation and (2) how do we best design the algorithms and software to meet these requirements? One avenue would be to refine and improve on system- or application-level checkpointing and rollback strategies in the case an error is detected. Developers might use fault notification interfaces and flexible runtime systems to respond to node failures in an application-dependent fashion. Novel numerical algorithms or more stochastic computational approaches may be required to meet accuracy requirements in the face of undetectable soft errors. The goal of this Dagstuhl Seminar was to bring together a diverse group of scientists with expertise in exascale computing to discuss novel ways to make applications resilient against detected and undetected faults. In particular, participants explored the role that algorithms and applications play in the holistic approach needed to tackle this challenge., Comment: 45 pages, 3 figures, submitted to The International Journal of High Performance Computing Applications

This paper presents a class of turbulence models written in terms of fractional partial differential equations (FPDEs) with stochastic loads. Every solution of these FPDE models is an incompressible velocity field and the distribution of solutions is Gaussian. Interaction of the turbulence with solid walls is incorporated through the enforcement of various boundary conditions. The various boundary conditions deliver extensive flexibility in the near-wall statistics that can be modelled. Reproduction of both fully-developed shear-free and uniform shear boundary layer turbulence are highlighted as two simple physical applications; the first of which is also directly validated with experimental data. The rendering of inhomogeneous synthetic turbulence inlet boundary conditions is an additional application, motivated by contemporary numerical wind tunnel simulations. Calibration of model parameters and efficient numerical methods are also conferred upon.

In this work, we present and analyse a system of coupled partial differential equations, which models tumour growth under the influence of subdiffusion, mechanical effects, nutrient supply, and chemotherapy. The subdiffusion of the system is modelled by a time fractional derivative in the equation governing the volume fraction of the tumour cells. The mass densities of the nutrients and the chemotherapeutic agents are modelled by reaction diffusion equations. We prove the existence and uniqueness of a weak solution to the model via the Faedo--Galerkin method and the application of appropriate compactness theorems. Lastly, we propose a fully discretised system and illustrate the effects of the fractional derivative and the influence of the fractional parameter in numerical examples.

In this work, we present and analyze a mathematical model for tumor growth incorporating ECM erosion, interstitial flow, and the effect of vascular flow and nutrient transport. The model is of phase-field or diffused-interface type in which multiple phases of cell species and other constituents are separated by smooth evolving interfaces. The model involves a mesoscale version of Darcy's law to capture the flow mechanism in the tissue matrix. Modeling flow and transport processes in the vasculature supplying the healthy and cancerous tissue, one-dimensional (1D) equations are considered. Since the models governing the transport and flow processes are defined together with cell species models on a three-dimensional (3D) domain, we obtain a 3D-1D coupled model. We show some mathematical results on the existence of weak solutions. Furthermore, simulation results are presented illustrating the evolution of tumors and the effects of ECM erosion.

Complex geometries as common in industrial applications consist of multiple patches, if spline based parametrizations are used. The requirements for the generation of analysis-suitable models are increasing dramatically since isogeometric analysis is directly based on the spline parametrization and nowadays used for the calculation of higher-order partial differential equations. The computational, or more general, the engineering analysis necessitates suitable coupling techniques between the different patches. Mortar methods have been successfully applied for coupling of patches and for contact mechanics in recent years to resolve the arising issues within the interface. We present here current achievements in the design of mortar technologies in isogeometric analysis within the Priority Program SPP 1748, Reliable Simulation Techniques in Solid Mechanics. Development of Non-standard Discretisation Methods, Mechanical and Mathematical Analysis.

The surrogate matrix methodology delivers low-cost approximations of matrices (i.e., surrogate matrices) which are normally computed in Galerkin methods via element-scale quadrature formulas. In this paper, the methodology is applied to a number of model problems in wave mechanics treated in the Galerkin isogeometic setting. Herein, the resulting surrogate methods are shown to significantly reduce the assembly time in high frequency wave propagation problems. In particular, the assembly time is reduced with negligible loss in solution accuracy. This paper also extends the scope of previous articles in its series by considering multi-patch discretizations of time-harmonic, transient, and nonlinear PDEs as particular use cases of the methodology. Our a priori error analysis for the Helmholtz equation demonstrates that the additional consistency error introduced by the presence of surrogate matrices is independent of the wave number. In addition, our floating point analysis establishes that the computational complexity of the methodology compares favorably to other contemporary fast assembly techniques for isogeometric methods. Our numerical experiments demonstrate clear performance gains for time-harmonic problems, both with and without the presence of perfectly matched layers. Notable speed-ups are also presented for a transient problem with a compressible neo-Hookean material.

In this work, we introduce an algorithmic approach to generate microvascular networks starting from larger vessels that can be reconstructed without noticeable segmentation errors. Contrary to larger vessels, the reconstruction of fine-scale components of microvascular networks shows significant segmentation errors, and an accurate mapping is time and cost intense. Thus there is a need for fast and reliable reconstruction algorithms yielding surrogate networks having similar stochastic properties as the original ones. The microvascular networks are constructed in a marching way by adding vessels to the outlets of the vascular tree from the previous step. To optimise the structure of the vascular trees, we use Murray's law to determine the radii of the vessels and bifurcation angles. In each step, we compute the local gradient of the partial pressure of oxygen and adapt the orientation of the new vessels to this gradient. At the same time, we use the partial pressure of oxygen to check whether the considered tissue block is supplied sufficiently with oxygen. Computing the partial pressure of oxygen, we use a 3D-1D coupled model for blood flow and oxygen transport. To decrease the complexity of a fully coupled 3D model, we reduce the blood vessel network to a 1D graph structure and use a bi-directional coupling with the tissue which is described by a 3D homogeneous porous medium. The resulting surrogate networks are analysed with respect to morphological and physiological aspects.

We propose a high-order discontinuous Galerkin scheme for nonlinear acoustic waves on polytopic meshes. To model sound propagation with and without losses, we use Westervelt's nonlinear wave equation with and without strong damping. Challenges in the numerical analysis lie in handling the nonlinearity in the model, which involves the derivatives in time of the acoustic velocity potential, and in preventing the equation from degenerating. We rely in our approach on the Banach fixed-point theorem combined with a stability and convergence analysis of a linear wave equation with a variable coefficient in front of the second time derivative. By doing so, we derive an a priori error estimate for Westervelt's equation in a suitable energy norm for the polynomial degree $p \geq 2$. Numerical experiments carried out in two-dimensional settings illustrate the theoretical convergence results. In addition, we demonstrate efficiency of the method in a three-dimensional domain with varying medium parameters, where we use the discontinuous Galerkin approach in a hybrid way.

In this article, we consider scenarios in which traditional estimates for the active subspace method based on probabilistic Poincar\'e inequalities are not valid due to unbounded Poincar\'e constants. Consequently, we propose a framework that allows to derive generalized estimates in the sense that it enables to control the trade-off between the size of the Poincar\'e constant and a weaker order of the final error bound. In particular, we investigate independently exponentially distributed random variables in dimension two or larger and give explicit expressions for corresponding Poincar\'e constants showing their dependence on the dimension of the problem. Finally, we suggest possibilities for future work that aim for extending the class of distributions applicable to the active subspace method as we regard this as an opportunity to enlarge its usability., Comment: 27 pages, 6 figures