HAZniCS Software Components for Multiphysics Problems

2025-05-06 0 0 4.94MB 31 页 10玖币

侵权投诉

HAZniCS – Soware Components for Multiphysics Problems

ANA BUDIŠA, Simula Research Laboratory, Norway

XIAOZHE HU, Department of Mathematics, Tufts University, USA

MIROSLAV KUCHTA, Simula Research Laboratory, Norway

KENT–ANDRÉ MARDAL, Department of Mathematics, University of Oslo, Norway

LUDMIL T. ZIKATANOV, Department of Mathematics, Penn State, USA

We introduce the software toolbox HAZniCS for solving interface-coupled multiphysics problems. HAZniCS

is a suite of modules that combines the well-known FEniCS framework for nite element discretization with

solver and graph library HAZmath. The focus of the paper is on the design and implementation of a pool of

robust and ecient solver algorithms which tackle issues related to the complex interfacial coupling of the

physical problems often encountered in applications in brain biomechanics. The robustness and eciency of

the numerical algorithms and methods is shown in several numerical examples, namely the Darcy-Stokes

equations that model ow of cerebrospinal uid in the human brain and the mixed-dimensional model of

electrodiusion in the brain tissue.

1 INTRODUCTION

The present paper aims to introduce a novel collection of tools for interface coupled multiphysics

problems modeled by partial dierential equations (PDEs). The interface is a main driver of the

processes in a way that strategies relying on decoupled single-physics problems typically suer

from slow convergence. Furthermore, we target multiphysics problems with geometrically complex

interfaces and slow dynamics – promoting monolithic solvers. Specically, we exploit fractional

operators and low-order interface perturbations as preconditioning techniques.

Fractional operators appear naturally on interfaces in multiphysics problems. One common ap-

proach has been using Poincaré-Steklov operators for uid-structure interaction problems [Agoshkov

1988;Deparis et al

2006;Quarteroni and Valli 1991], which exploits Dirichlet-to-Neumann map-

pings. As the Poincaré-Steklov operator takes functions in the fractional Sobolev space

𝐻1/2

functions in its dual

𝐻−1/2

, it is equivalent to a fractional Laplacian operator

(−Δ)1/2

. However, the

Poincaré-Steklov operator is not sucient for parameter-dependent problems as it is sensitive to

problem parameters, and often many sub-iterations are required. More sophisticated techniques

that include problem parameters such as Robin-to-Dirichlet, -Neumann, or -Robin maps have

been explored [Badia et al

2009], but the approach still requires tuning. We remark that the

Poincaré-Steklov operator involves the extension to a domain in a higher dimension and is, as

such, computationally expensive. However, the computational complexity is usually the same as

the involved single-physics problems.

As an alternative or generalization of Poincaré-Steklov operators, several recent papers [Boon

et al

2022a,b;Holter et al

2021;Kuchta et al

2021] have considered multiphysics problems and

derived order-optimal and parameter-robust algorithms. They are obtained by exploiting fractional

Laplacians (or sums thereof) and metric terms on the interface. Here, the fractional Laplacians

arise naturally due to trace operators appearing in the coupling conditions, e.g., conservation

of mass, that connect the unknowns of the dierent single-physics problems. We note that the

Authors’ addresses: Ana Budiša, ana@simula.no, Simula Research Laboratory, P.O. Box 134, 1325, Lysaker, Norway;

Xiaozhe Hu, Department of Mathematics, Tufts University, 503 Boston Avenue, Medford, 02155, Massachusetts, USA,

xiaozhe.hu@tufts.edu; Miroslav Kuchta, Simula Research Laboratory, P.O. Box 134, 1325, Lysaker, Norway, miroslav@

simula.no; Kent–André Mardal, Department of Mathematics, University of Oslo, P.O. Box 1053, Blindern, 0316, Oslo, Norway,

kent-and@math.uio.no; Ludmil T. Zikatanov, Department of Mathematics, Penn State, 239 McAllister Building, University

Park, 16802, Pennsylvania, USA, ludmil@psu.edu.

arXiv:2210.13274v2 [math.NA] 6 Nov 2022

0:2 A. Budiša, X. Hu, M. Kuchta, K.-A. Mardal and L. T. Zikatanov

fractional operators arise both when Lagrange multipliers are used to prescribe the interface

conditions, e.g., [Holter et al

2020;Layton et al

2002], and when they are avoided [Boon et al

2022a,b]. The metric terms then arise because interface conditions, such as the balance of forces,

are often expressed in terms of dierences of a quantity (e.g., displacement) across the interface

rather than the quantity itself [Boon et al

2022b;D’Angelo and Quarteroni 2008;Kuchta et al

2021]. Recently, fast solution algorithms for fractional Laplacians have been proposed based on

multilevel approaches [Bærland 2019;Bærland et al

2019;Bramble et al

2000;Führer 2022;Zhao

et al

2017] and rational approximations [Harizanov et al

2020,2018,2022]. Here, we explore the

latter for sums of fractional Laplacians. In addition to rational approximations, we will consider

multilevel algorithms that work robustly in the presence of strong metric terms at interfaces. That

is multilevel algorithms with a space decomposition aware of the metric kernel.

The software tools we developed aim to solve computational mesoscale multiphysics problems.

By computational mesoscale, in this context, we refer to problems in the range of a few hundred

thousand to tens of millions of degrees of freedom. These problems do not require parallel computing,

but they may benet signicantly from advanced algorithms. The collection of tools presented

in this paper are FEniCS [Logg et al

2012] add-ons for block assembly [Kuchta 2021] and block

preconditioning [Mardal and Haga 2012] combined with a exible algebraic multigrid (AMG)

toolbox, implemented in C, called HAZmath [Adler et al

2009]. Hence, we have named the tool

collection HAZniCS. One of the reasons for developing HAZniCS is precisely the mentioned

exibility and variety of the implementation of the AMG method in HAZmath. It allows us to easily

modify available linear solvers and preconditioners or create new model-specic solvers for the

multiphysics problems at hand. Additionally, with HAZniCS, we provide another wide range of

ecient computational methods for solving PDEs with FEniCS, but also a bridge to Python for

HAZmath to be used with other PDE simulation tools. Further in the paper, we highlight with a

series of code snippets the implementation of several solvers, namely the aggregation-based and

metric-perturbed AMG methods and the rational approximation method.

Moreover, we consider a series of examples of multiphysics problems mainly related to biome-

chanical processes. Namely, we include: (1) a simple three-dimensional example of a elliptic problem

on a regular domain, (2) Darcy-Stokes equations describing the interaction of the viscous ow

of cerebrospinal uid ow surrounding the brain and interacting with the porous media ow of

interstitial uid inside the brain, and (3) the mixed-dimensional equations representing electric

signal propagation in neurons and the surrounding matter.

The outline of the current paper is as follows: in Section 2we introduce the multiphysics models

together with the necessary mathematical concepts and numerical methods. Section 3focuses

on the implementation of those methods and the interface between the software components. In

Section 4we present the solver capabilities of our software to simulate relevant biomechanical

phenomena. Finally, we draw concluding remarks in Section 5.

2 EXAMPLES

The following three examples illustrate dierent single- and multiphysics PDE models, as well as

the relevant mathematical and computational concepts. More specically, the examples provide

an overview of iterative methods and preconditioning techniques for interface-coupled problems

that lead, e.g., to the utilization of sums of fractional operators weighted by material parameters.

Additionally, we include several code snippets in each example that highlight the most important

features of the implementation, while the full codes can be found in [Budiša et al. 2022a].

HAZniCS – Soware Components for Multiphysics Problems 0:3

2.1 Linear elliptic problem

We start with a linear elliptic problem on a three-dimensional (3

𝑑

) regular domain. This example

will serve as a baseline for the solvers in HAZniCS. Our goal is to demonstrate that our solver

performance is comparable to other established software. Additionally, the solution methods that

we use here will be incorporated and adapted to the multiphysics problems in the later examples.

Let

Ω=[

be the unit cube and let

𝜕Ω

denote its boundary. Given external force

𝑓

Ω→R

and the boundary data 𝑔:𝜕Ω→R, we set to nd the solution 𝑢:Ω→Rthat satises

−Δ𝑢+𝑢=𝑓in Ω,(1a)

𝜕𝑢

𝜕𝒏=𝑔on 𝜕Ω.(1b)

To solve

(1)

computationally, we relate to

(1)

the variational formulation and the discrete problem

using the nite element method (FEM). First, let

𝐿2=𝐿2(Ω)

be the space of square-integrable

functions on

and

𝐻𝑠=𝐻𝑠(Ω)

the Sobolev spaces with

𝑠

derivatives in

𝐿2

. The corresponding

inner products and norms for any function space

𝑋

are denoted with

(·,·)𝑋

and

∥ · ∥𝑋

, respectively.

Furthermore, we let ⟨·,·⟩ denote a duality pairing between 𝑋′, the dual of 𝑋and 𝑋.

Now, let

𝑉ℎ⊂𝐻1(Ω)

be a nite element space on triangulation of

, e.g., of continuous piecewise

linear functions (P1) . A discrete variational formulation of (1) states to nd 𝑢∈𝑉ℎsuch that

𝑎(𝑢, 𝑣)=⟨𝑓 , 𝑣⟩ ∀𝑣∈𝑉ℎ,(2)

with 𝑎(𝑢, 𝑣)=(𝑢, 𝑣)𝐿2(Ω)+ (∇𝑢, ∇𝑣)𝐿2(Ω)and ⟨𝑓 , 𝑣⟩=(𝑓 , 𝑣)𝐿2(Ω)+ (𝑔, 𝑣)𝐿2(𝜕Ω).

Furthermore, it is important for the solvers to obtain matrix-vector representation of the above

system. Let the discrete operator 𝐴:𝑉ℎ→𝑉′

ℎsatisfy

⟨𝐴𝑢, 𝑣⟩=𝑎(𝑢, 𝑣), 𝑢, 𝑣 ∈𝑉ℎ(3)

with

𝑉′

ℎ

denoting the dual space of

𝑉ℎ

. Its actual implementation can be derived as follows. Let

𝜓𝑖, 𝑖 =

, . . . , 𝑚

be the nite element basis functions of

𝑉ℎ

. Dene matrix A

∈R𝑚×𝑚

and vectors

f∈R𝑚as

(A)𝑖 𝑗 =⟨𝐴𝜓𝑗, 𝜓𝑖⟩,and f𝑖=⟨𝑓 ,𝜓𝑖⟩,for 𝑖, 𝑗 =1,2, . . . , 𝑚. (4)

Consequently, we get the discrete system of equations related to

(1)

, i.e. we aim to solve for u

∈R𝑚

the algebraic system

Au =f.(5)

We remark that uand fabove are both vectors in

R𝑚

, but that uis in the so-called nodal represen-

tation, i.e.

𝑢=Í𝑖

𝑖𝜓𝑖

, while fis in the dual representation [Bramble 2019;Mardal and Winther

2011]. As such, the matrix Amaps the nodal representations of R𝑚to its dual representation.

Since Ais symmetric positive denite (SPD), we solve

(5)

with the Conjugate Gradient (CG)

method. It is well known that the number of iterations of a Krylov iterative scheme can be bounded

in terms of the condition number of the system, that is

𝜅(

)=|||A||||||A−1|||

for some matrix norm

||| · |||

. Therefore, to eciently solve the problem

(5)

we want as few iterations as possible and order

optimal scalability of the solver with regards to the number of degrees of freedom. To that aim, we

introduce a preconditioner Bsuch that

𝜅(BA)=|||BA||||||(BA)−1||| ≈ O(1),(6)

that is,

𝜅(

)

stays bounded from above independently of discretization and other problem parame-

ters. It is important to make sure that Bmaps dual representations of vectors to nodal representations

as the preconditioner Bis an approximation of the inverse of A.

The previous result is also true in the general case for symmetric operators on Hilbert spaces.

Specically, for

𝐴

(3)

we nd an operator

𝐵

𝑉′

ℎ→𝑉ℎ

such that

𝜅(𝐵𝐴)

is uniformly bounded,

0:4 A. Budiša, X. Hu, M. Kuchta, K.-A. Mardal and L. T. Zikatanov

where the matrix norm is replaced with the operator norm in the space of continuous linear

operators dened on

𝑉ℎ

. In context of operator preconditioning [Mardal and Winther 2011], a

common choice for the preconditioner operator 𝐵is the Riesz mapping, that is

(𝐵𝑓 , 𝑣)𝑉ℎ=⟨𝑓 , 𝑣⟩,∀𝑓∈𝑉′

ℎ, 𝑣 ∈𝑉ℎ.(7)

The Riesz map guarantees a uniform bound on

𝜅(𝐵𝐴)

when

𝐴

is a bounded operator that satises

the inf-sup conditions [Babuška 1971;Babuška and Aziz 1972] independent of system parameters,

such as in the case of the operator in

(3)

. Moreover, we can use any other preconditioner that gives

a uniform bound on the condition number. If we nd a spectrally equivalent operator

𝐵𝑆𝐸

such that

for parameter-independent constants 𝑐1, 𝑐2>0it satises

𝑐1∥𝑣∥2

𝐴≤ ∥𝑣∥2

𝐵−1

𝑆𝐸 ≤𝑐2∥𝑣∥2

𝐴,(8)

with

∥𝑣∥2

𝐴=⟨𝐴𝑣, 𝑣⟩

, then we retain a uniform bound on the condition number

𝜅(𝐵𝑆𝐸𝐴) ≤ 𝑐2

𝑐1𝜅(𝐵𝐴)

This is relevant when an application of

𝐵

on a function in

𝑉′

ℎ

is infeasible or inecient. We want to

replace it with a method that applies a spectrally equivalent operation. In the rest of the paper, we

will note

𝜅(𝐵𝐴)

as the condition number for both operators (

𝐴, 𝐵

) and their matrix representations

(A,B), clarifying along the way if ambiguity occurs.

In our case, it is well-known that multilevel methods, such as AMG, provide spectrally equivalent

and order optimal algorithms for the inverse of discretizations of

𝐼−Δ

. Thus, we dene the

preconditioner for (5) as B=AMG(A).

The implementation of the elliptic problem in FEniCS follows straightforwardly from the varia-

tional formulation (2) and is one of the basic examples of FEniCS software, see Listing 1.

from block.iterative i mport ConjGrad

from blo ck . alge br aic . hazm at h import AMG

from dolfin import *

mesh = UnitCubeMesh(32 , 32 , 3 2)

V = FunctionSpace( mesh , "CG", 1)

u , v = TrialFunction( V ) , TestFunction(V)

f = Expression(" s in ( p i *x [ 0] ) " , degree =4)

a = inne r (u , v ) * dx +i nn er (grad( u) , grad(v)) * dx

L = inne r (f , v ) * dx

A = assemble(a)

b = assemble(L)

B = AMG(A, parameters={"max_levels": 10 , "AMG_type": 1})

Ainv = ConjGrad( A , p re co nd = B , t ol er an c e =1 e - 10 )

x = A inv * b # So lve for the coef fi ci en t v ect or of foo in V

Listing 1. Implementation of the linear elliptic problem

(1)

. Complete code can be found in scripts

HAZniCS-examples/demo_elliptic*.py

For the preconditioner, we utilize the AMG method implemented in HAZmath, available through

our HAZniCS library. We describe the AMG method and its implementation in more detail in

Section 3.1 and showcase the performance of HAZmath AMG as compared with HYPRE [Falgout

and Yang 2002] AMG in Section 4.1.

2.2 Modeling brain clearance during sleep with Darcy-Stokes equations

We consider a multiphysics problem arising in modeling processes of waste clearance in the brain

during sleep [Eide et al

2021;Xie et al

2013] with potential links to the development of Alzheimer’s

HAZniCS – Soware Components for Multiphysics Problems 0:5

Ω

yCρ(y)

Fig. 1. (Le) Geometry and computational mesh from [Boon et al

2022a] of the Darcy-Stokes model of brain

clearance. Mesh and indicator functions for tracking subdomains making up the Darcy-(light blue) and the

Stokes domains (dark blue and orange subregions) and their interfaces are generated with SVMTK [Mardal

et al

2022]. (Right) Model reduction from 3

𝑑

-3

𝑑

to 3

𝑑

-1

𝑑

problem. Dendrites (in blue) are reduced to their

centerline while the coupling with the surrounding

is accounted for by averaging over-idealized cylindrical

surfaces with radius 𝜌.

disease. The novel model, called the glymphatic model [Ili et al

2012], states that the viscous

ow of cerebrospinal uid (CSF) is tightly coupled to the porous ow in the brain tissue and that

during sleep, in particular, it clears metabolic waste from the brain, for computational models see

e.g. [Boon et al

2022a;Holter et al

2020;Kedarasetti et al

2020]. To this end we will consider

patient-specic geometries generated from MRI images by SVTMK library [Mardal et al

2022]

used in [Boon et al

2022a], see Figure 1. Using SVMTK, the segmented brain geometry is enclosed

in a thin shell, which, together with the ventricles (the orange subregion in Figure 1), makes up the

Stokes domain. We remark that the diameter of the Stokes domain is roughly 15 cm while the shell

thickness is on average 0.8 mm.

In order to model the waste clearance, let

Ω𝐷⊂R𝑑

𝑑=

3be the domain of the porous medium

ow that represents the brain tissue

, and let

Ω𝑆⊂R𝑑

be the domain of viscous ow representing

the subarachnoid space around it saturated by CSF. Let

denote the interface between the domains,

which in this case corresponds to the surface of the brain. We then consider the Darcy-Stokes

model which seeks to nd Stokes velocity

𝒖𝑆

Ω𝑆→R𝑑

and pressure

𝑝𝑆

Ω𝑆→R

, and Darcy

velocity 𝒖𝐷:Ω𝐷→R𝑑and pressure 𝑝𝐷:Ω𝐷→Rthat satisfy

−∇ · 𝝈𝑆(𝒖𝑆, 𝑝𝑆)=𝒇𝑆in Ω𝑆,(9a)

∇ · 𝒖𝑆=0in Ω𝑆,(9b)

𝒖𝐷=−𝐾∇𝑝𝐷in Ω𝐷,(9c)

∇ · 𝒖𝐷=𝑓𝐷in Ω𝐷,(9d)

with interface conditions

𝒖𝑆·𝒏−𝒖𝐷·𝒏=0on Γ,(9e)

𝒏·𝝈𝑠(𝒖𝑆, 𝑝𝑆) · 𝒏+𝑝𝐷=0on Γ,(9f)

𝒏·𝝈𝑠(𝒖𝑆, 𝑝𝑆) ·𝝉+𝐷𝒖𝑆·𝝉=0on Γ.(9g)

1In the context of brain mechanics, the case 𝑑=2is relevant, e.g., for the slices of the brain geometry.

文档加载中……请稍候！
如果长时间未打开，您也可以点击刷新试试。

下载文档到电脑，查找使用更方便

10 玖币 0人已下载

立即下载

摘要：

HAZniCS–SoftwareComponentsforMultiphysicsProblemsANABUDIŠA,SimulaResearchLaboratory,NorwayXIAOZHEHU,DepartmentofMathematics,TuftsUniversity,USAMIROSLAVKUCHTA,SimulaResearchLaboratory,NorwayKENT–ANDRÉMARDAL,DepartmentofMathematics,UniversityofOslo,NorwayLUDMILT.ZIKATANOV,DepartmentofMathematics,PennS...

展开>> 收起<<

HAZniCS Software Components for Multiphysics Problems.pdf

共31页,预览5页

还剩页未读，继续阅读

声明：本站为文档C2C交易模式，即用户上传的文档直接被用户下载，本站只是中间服务平台，本站所有文档下载所得的收益归上传人(含作者)所有。玖贝云文库仅提供信息存储空间，仅对用户上传内容的表现方式做保护处理，对上载内容本身不做任何修改或编辑。若文档所含内容侵犯了您的版权或隐私，请立即通知玖贝云文库，我们立即给予删除！

HAZniCS Software Components for Multiphysics Problems

相关推荐

开通VIP享超值会员特权

作者详情

相关内容

热门标签

举报选择: