A PDE-constrained optimization method for 3D-1D coupled problems with discontinuous solutions

A numerical method for coupled 3D-1D problems with discontinuous solutions at the interfaces is derived and discussed. This extends a previous work on the subject where only continuous solutions were considered. Thanks to properly defined function spaces a well posed 3D-1D problem is obtained from the original fully 3D problem and the solution is then found by a PDE-constrained optimization reformulation. This is a domain decomposition strategy in which unknown interface variables are introduced and a suitably defined cost functional, expressing the error in fulfilling interface conditions, is minimized constrained by the constitutive equations on the subdomains. The resulting discrete problem is robust with respect to geometrical complexity thanks to the use of independent discretizations on the various subdomains. Meshes of different sizes can be used without affecting the conditioning of the discrete linear system, and this is a peculiar aspect of the considered formulation. An efficient resolution strategy is further proposed, based on the use of a gradient based solver and yielding a method ready for parallel implementation. A numerical experiment on a problem with known analytical solution shows the accuracy of the method, and two examples on more complex configurations are proposed to address the applicability of the approach to practical problems.


Introduction
The present work deals with simulations in domains with embedded cylindrical, or nearly cylindrical inclusions, with cross-section sizes much smaller than their length and than the domain scale.This kind of problem is typical of a large variety of applications, ranging from the study of living tissues, where the inclusions are constituted, e.g., by the capillaries or by the vessels of the lymphatic system [2,17,12], to the study of the interaction between plant roots and the soil [19,10], or fibre-reinforced materials [20,16], and geological applications [7,6,3].In these cases, it is often convenient to treat the inclusions as one dimensional objects, thus actually collapsing the cross-sections on the centrelines, avoiding the complexity and the overhead of generating good quality meshes in the interior of the small inclusions.Such geometrical reduction is however non-trivial from a mathematical standpoint as it ends up in a 3D problem which is coupled to problems on the 1D domains, and non standard function spaces need to be used to allow the definition of a well posed trace operator for 3D functions on 1D manifolds.
An analysis of 3D problems with singular source terms is available in Refs.[4,5], where the solution is placed in suitable weighted Sobolev spaces.In Ref. [4] a finite element based approach is used for the resolution of the problem and optimal convergence trends for the error are observed in the used spaces.Problems with singular source terms are also considered in Ref. [21], where Dirac delta sources are replaced by regularizing terms, compactly supported.Regularizing techniques are suggested also in Ref. [9], and in Ref. [11] where line source terms are approximated by suitable kernel functions that distribute the source in a three dimensional neighbourhood of the line.A splitting technique is proposed in Ref. [8], where the solution is seen as the sum of a known low regularity term and a regular correction term that is computed solving an elliptic problem with source term and boundary data depending on a chosen extension operator for the singular source.In Ref. [13] a lifting technique is employed for the singular source term, whereas a domain decomposition approach based on the use of Lagrange multipliers is suggested in Ref. [14].The derivation of a coupled 3D-1D problem, starting from an original equidimensional formulation is proposed in Ref. [15], under suitable assumptions on the behavior of the solution in the inclusions that allow the definition of averaging operators.
Here an extension of the approach described in Ref. [1] for 3D-1D coupled elliptic problems is proposed, allowing to deal with discontinuous solution at the interfaces.This is relevant for a large variety of practical applications, as e.g. the description of biological tissues, where the boundary between the 3D and the 1D domains is a semi-permeable membrane and the Starling equation applies [2,15,11,14].A well posed 3D-1D coupled problem is derived from the original 3D-3D problem through an appropriate choice of the functional space for the solution, which allows to define extension and trace operators between spaces on 3D and 1D domains.Further, a domain decomposition approach is employed, introducing additional interface unknowns to decouple the 3D problem from the problems on the 1D inclusions, and a cost functional is designed such that the solution is obtained as the minimum of the functional constrained by the constitutive equations on the subdomains.The original work in Ref. [1] deals with continuous solution at the boundary between the 3D and the 1D domains.Here, while keeping the same structure of this original method, different interface variables are introduced for the domain decomposition process, which result in a novel setting for the PDE constrained optimization problem described.The problem is discretized using finite elements on non conforming meshes, thus providing a great flexibility in the choice of the meshes that can be independently defined on each subdomain, and a numerical scheme suitable for parallel computing is obtained thanks to the use of domain decomposition combined to a gradient scheme for the resolution of the resulting discrete problem.
The problem considered is presented in Section 2 in equi-dimensional form, and in Section 3 its 3D-1D formulation is derived in weak form in ad-hoc function spaces.The problem is re-written as a PDE constrained optimization problem in Section 4, and presented in discrete form in Section 5.The approach used to solve the obtained linear system is discussed in Section 6, whereas some numerical tests are reported in Section 7 and finally conclusions are reported in Section 8.

Notation and formulation of the fully 3D coupled problem
Let Ω⊂R 3 be a convex domain in which a generalized cylinder Σ⊂R 3 is embedded.We denote by Λ = {λ(s), s ∈ (0, S)} the centerline of Σ, while Γ = {Γ(s), s ∈ [0, S]} is the lateral surface of Σ.In the following we assume, for simplicity, that Λ is a rectilinear segment in the three-dimensional space.We denote by Σ(s) the transversal sections of the cylinder as s ranges in the interval [0, S] and by Γ(s) their boundary.We suppose the sections to have an elliptic shape, with R(s) being the maximum axes length of the ellipse centered in λ(s).For the two extreme sections of the cylinder we adopt the compact notation Σ 0 = Σ(0) and Σ S = Σ(S).For the derivation of the model problem we assume that Σ 0 and Σ S lie on the boundary ∂Ω, but the extension to more general cases is straightforward.The portion of Ω that does not include the cylinder is denoted by D = Ω \ Σ.We define ∂D e = ∂Ω \ {Σ 0 ∪ Σ S }, referring to it as the external boundary of D, with ∂D = ∂D e ∪ Γ.In case the extreme sections of Σ were inside Ω, ∂D e would coincide with ∂Ω.
Let us now consider the following diffusion problem, with unknown pressures u in D and ũ in Σ: 3D-problem on Σ: Vectors n and ñ = −n are the outward pointing unit normal vectors to Γ, respectively for D and Σ; K, K and β are positive scalars, while f and g denote source terms.For the sake of simplicity we consider homogeneous Dirichlet boundary conditions on Σ 0 and Σ S and on ∂D e .Equations ( 3) and ( 6) allow us to couple the two problems imposing flux conservation.According to these equations, the flux across Γ is directly proportional to the jump of the pressures, with β denoting the permeability coefficient of the membrane Γ.Different coupling conditions could be considered, for example adding a pressure continuity constraint and consequently not linking the flux definition to the pressure jump, as done in Ref. [1].The choice of the interface condition depends of course on the properties of the interface, and thus on the kind of application.
Let us now suppose that R is much smaller than the size of Ω and than the longitudinal length L of the cylinder itself, in particular.This allows us to assume that the variables defined on Σ or on Γ are actually only functions of the coordinate s, considering negligible their variation on the cross-sections of the inclusion.This is the key point that allows us, in the next section, to work out a well-posed 3D-1D coupled problem from equations ( 1)-( 6).

Variational formulation of the 3D-1D problem
A 3D-1D coupled problem is obtained from problem (1)-( 6), after writing it in variational form in suitable function spaces, as here described.Let us, thus, introduce the spaces and two extension operators such that, given v ∈ H 1 0 (Λ), E Σ v(s) and E Γ v(s) are the uniform extension of v(s) respectively to Σ(s) and to Γ(s).We observe that E Γ = γ Γ • E Σ .Once denoted by V the space H 1 0 (Λ), let us further consider the spaces: such that functions in V and in H Γ are respectively the uniform extension to Σ and Γ of functions in V and functions in V D have trace on Γ belonging to H Γ .
Functions in such spaces fit the assumptions we have made on the negligible variation of the variables on the cross-sections of Σ and Γ. Denoting by (•, •) X the scalar product on a generic space X, the variational formulation of problem ( 1)-( 6) can be written as: ( K∇ũ, ∇ṽ) Let us introduce two auxiliary variables ψ D , ψ Σ ∈ H Γ , in order to formally decouple the two equations.Denoting by X the dual of the generic space X, the problem is thus rewritten as: ( K∇ũ, ∇ṽ) Let us remark that Equations ( 9)-( 10) could also be written as However, formulation (9)-( 10) is preferred, as it allows to have an empty Dirichlet boundary on either ∂D e or Σ 0 , Σ s .This is a desired property for domain decomposition purposes.Thanks to the adopted functional spaces, problem ( 9)-( 12) can be easily reduced to a 3D-1D coupled problem.Let us observe that, given η ∈ H Γ and ρ ρη dl ds.

PDE-constrained optimization problem
Conditions ( 15) and ( 16) can be replaced by the minimization of a cost functional mimicking the error committed in the fulfillment of such constraints.At this aim let us define to be minimized constrained by ( 13) and (14).In order to work out the PDEconstrained optimization formulation of the problem in a compact form, let us define the linear operators A : The respective adjoints will be denoted as equations ( 13)-( 14) can be written as: Finally, the PDE-constrained optimization problem can be written as We now provide some results on the optimal control and the stepsize of the steepest descent method for Problem (26).

Proposition 1. Let us consider the trace operator γ
Then the optimal control ( ψD , ψΣ ) that provides the solution to (26) is such that where p ∈ V D and p ∈ V are the solutions respectively to Proof.Let us compute the Fréchet derivatives of J with respect to the control variables ψD and ψΣ .To this end, we introduce the increments δ ψD , δ ψΣ ∈ V and we recall that, for = D, Σ, there exists δψ ∈ H Γ : δψ = E Γ δ ψ .We have: which yield the thesis.
From the derivatives computed in Proposition 1, we now define the quantities Then the following proposition holds: Proposition 2. Given the variable X , let us increment it by a step ζδX , where δX = (δ ψD , δ ψΣ ).The steepest descent method corresponds to the stepsize and where and δp ∈ V D , δ p ∈ V are such that: Proof.It is sufficient to set to zero the derivative ∂J(X +ζδX ) ∂ζ .In the computation that follows we adopt the lighter notation: Rearranging properly the terms we get ζ = − N D with and that yields the thesis.

Discrete matrix formulation
In this section we work out the discrete matrix formulation of problem (26).In general, the 3D-1D coupling does not present particular issues in the discrete framework.Nonetheless our approach has the additional advantage of allowing for the use of non conforming meshes: thanks to the optimization framework, the partitions of the 1D inclusions can be defined in a completely independent manner from the surrounding 3D mesh, without any theoretical or practical constraint on mesh sizes.Further, the proposed formulation provides the direct computation of interface variables, and it allows to decouple the 3D problem from the 1D problems, thus paving the way to the use of possibly different constitutive equations and to the application of efficient solvers based on parallel computing techniques.
For the sake of generality, we consider I segments of different length and orientation crossing the domain Ω.The segments are defined as Λ i = {λ i (s), s ∈ (0, S i )}, i = 1, ..., I and they represent the centerlines of I cylindrical inclusions Σ i .
The proposed approach can easily handle intersections among inclusions centrelines.Intersecting segments are split into sub-segments in correspondence of their intersection point q.In this way, q always corresponds to a segment endpoint, in which pressure continuity and flux conservation are constrained.It is to remark that a variety of intersection modes is possible for the original three dimensional inclusions.As an example, 3D inclusions might partially overlap whereas the corresponding centrelines might not intersect.By considering here only intersections between centrelines, we implicitly assume that the intersection volume of the corresponding three dimensional inclusions is small and can be reduced to a point in the scale of the domain.A deeper investigation on the treatment of different intersection models is out of the scope of the present work.
After having extended the domain D to the whole Ω, let us consider a tetrahedral mesh T of domain Ω, on which we define Lagrangian finite element basis functions {ϕ k } N k=1 , such that U = N k=1 U k ϕ k is the discrete approximation of pressure u.On each segment Λ i we build three different partitions, independent from each other and from T .We denote them by Ti , τ D i and τ Σ i and we define the basis functions on τ Σ i , with Ni , N D i and N Σ i denoting the number of DOFs of the discrete approximations of the variables ûi , ψD ,i and ψΣ,i respectively.Such approximations are defined as We then define the following matrices: and the vectors we can group the matrices as follows: where matrix Q simply equates the DOFs at the extrema of connected subsegments and allow us to enforce continuity through Lagrange multipliers.Let us observe how Â = Â in case no intersections occur among segments.Finally we can write with The discrete functional is derived from ( 17) replacing the norms in H Γ with norms in L 2 (Λ) and summing over the I inclusions.First we define matrices and then The discrete cost functional then reads: Finally, the discrete matrix formulation of the 3D-1D problem can be written as: First order optimality conditions for problem (36) are collected in the saddlepoint system Proposition 3. Matrix K in (37) is non-singular and the unique solution of (37) is equivalent to the solution of the optimization problem (36).
The following lemma is used to prove Proposition 3.
and let G ∈ R (N + N +N D +NΣ)×(N + N +N D +NΣ) be defined as Then matrix A is full rank and matrix G is symmetric positive definite on ker(A).
Proof.Matrix A is full rank for the ellipticity of operators A in (18) and Â in (19), whereas matrix G is symmetric positive semi-definite as be the k-th element of the canonical basis of R N D +NΣ and let z k ∈ ker (A) be defined as: Thus it is either e D k = 0 and z Σ k = 0, either e Σ k = 0 and z D k = 0, and consequently for a certain segment Λ, depending on k.As a consequence, z k ∈ ker (G) for any k = 1, . . ., N D + N Σ .The vector space ker (A) = span{z 1 , . . ., z N D +N Σ } is a subspace of Im(G), and ker (G) ∩ ker (A) = {0}.
The proof of Proposition 3 derives from classical arguments of quadratic programming, observing that 6 Solving strategies Solving system (37) is equivalent to solve the optimum problem (36).A different resolution approach is however proposed, based on an iterative solver and allowing to take full advantage of the decoupling introduced by the proposed method.
Let us formally replace in the cost functional (35) the expressions Matrix M is symmetric positive definite, given the equivalence of this formulation with the saddle-point system (37).This allows us to perform the minimization of the unconstrained functional (38) via a gradient based scheme, looking for the minimum as the solution of A preconditioner P can be defined for the resolution of system (39).In particular we set where Â = diag Â1 , ..., ÂI .This means that the top-left block of matrix P corresponds exactly to the top-left block of matrix M in case no intersections among the segment occur.Otherwise, ( Â ) −1 is an approximation of the inverse of matrix Â which can be built inverting independently the matrices related to the single segments and which maintains a block-diagonal structure, i.e. ( Â For what concerns the bottom-right block, only matrix M Σ is kept with respect to the same block of matrix M, so that even this portion of the preconditioner can be built block-diagonalwise, assembling matrices which are independently related to each single 1D inclusion.The conjugate gradient scheme which is employed to solve system (39) is reported in Algorithm 1.The quantity MδX , whose computation is required at each iteration in the algorithm, can actually be obtained without explicitly building matrix M. In fact, after some computations we obtain where δX = [δΨ T D , δΨ T Σ ], and δU , δ Û , δP , δ P are the solutions of the linear systems which require the resolution of local sub-problems on the 1D segments and on the 3D domain.

Numerical results
In this section we present three numerical examples to better highlight the characteristics of the proposed approach.The simulations are performed using linear finite elements on the 3D and 1D non-conforming meshes, independently Algorithm 1: Conjugate gradient method for MX + d = 0 Solve Pz k+1 = r k+1 ; 10 13 end generated on the sub-domains.Parameter h denotes the maximum diameter of the tetrahedra for the 3D mesh, while other three parameters, namely δu,i , δ D,i and δ Σ,i express the refinement level of the 1D meshes Ti , τ D i , τ Σ i , i = 1, ...I, respectively.Each of these three parameters represents the ratio between the number of nodes in the 1D mesh and the number of intersections of the segment Λ i with the faces of the tetrahedra in T .In the simulations, for simplicity, we adopt unique, but possibly different, values of δu , δ D and δ Σ for the various segments: for this reason we drop, in the following, the segment index i = 1, ...I for these parameters.In all cases, linear Lagrangian finite elements on tethrahedra are used in the 3D domains and piecewise continuous linear basis functions on equally spaced meshes are chosen for the 1D functions.

Test Problem 1 (TP1)
Let us consider a cube Ω of edge l = 2 centered in the axes origin and whose faces are parallel to the coordinate axes.Let us further consider a cylinder Σ of radius R = 10 −2 and height h = 2 whose centreline Λ lies on the z axis (see Figure 1, on the left).Let us denote by ∂Ω l , ∂Ω + and ∂Ω − respectively the lateral, the top and the bottom faces of the cube.We aim at solving a problem in the form of ( 13)-( 16), obtained by reducing Σ to its centerline, with the following data: The problem is completed with appropriate boundary conditions such that the exact solution is: In particular we consider Neumann boundary conditions on ∂Ω + and ∂Ω − , whereas Dirichlet boundary conditions are imposed on ∂Ω l .Dirichlet boundary conditions equal to 1 are imposed at segment endpoints.The obtained solution is shown in Figure 1: on the left, the 3D solution is shown on a portion of the domain, whereas, on the right, the solution U on the y − z plane containing the z-axis is plotted along with the solution Û , with solution values reported along the x-axis.The computational mesh used for this solution has parameters h = 0.083, δu = 1, δ D = 0.5 and δ Σ = 0.5, corresponding to N = 4155 DOFs in the cube and N = 57 DOFs on the segment.
Errors indicators E L 2 , E H 1 are chosen for the 3D solution and E L 2 and E H 1 for the 1D problem, defined as: ., . where The values of these errors indicators on the same meshes considered before, are reported in Figure 3, on the left for E D ψ and, on the right, for E Σ ψ .Again optimal convergence curves are obtained, if it is considered that Ψ D is the approximation of a quantity in H 1 2 (Λ) and Ψ Σ the approximation of a function in H 1 (Λ).
Figure 4 shows the trend of the condition number of matrix K, defined in (37), under the variation of the 1D-mesh parameters.On the left the conditioning is plotted under the variation of δ D and for different values of δu , while a constant δ Σ = 0.5 is used; on the right δ Σ varies instead, while δ D = 0.5 is fixed, again for δu ranging between 0.05 and 2. In both cases we can observe how, in general, at the increase δu , slightly higher conditioning values are registered, with some exceptions for very small values of the parameter.In all cases, however, the impact of this parameter is quite marginal on the conditioning of the system.Looking at the left plot we can see that the value of δ D has no impact on the conditioning.Looking instead at Figure 4 on the right, it can be noticed that δ Σ has a larger impact on the conditioning, but only if very small values are used, and, at the same time, a value δu > δ Σ is chosen; in these cases an increase of up to two orders of magnitude in the conditioning is observed.However, for δ Σ > 0.2 the effect of δ Σ on the conditioning becomes almost irrelevant, independently from the choice of the other parameters.The behaviour here observed is quite different from the one observed in Ref. [1], where a larger effect of the mesh parameters on the conditioning was instead observed.Further, system (39) is known to be even better conditioned than the corresponding system (37), [18].For this simple example it is possible to explicitly compute matrix M. Its conditioning is plotted in Figure 5, against variations of the 1D mesh parameters.Trends similar to the ones of Figures 4 are observed and the behaviour is almost the same for all the values of δu considered, but the conditioning of M is, in general, between 4 and 5 orders of magnitude smaller than the one of K.This is expected to have a positive impact on the number of iterations of Algorithm 1, with few iterations required to reach the prescribed tolerance even without the use of a preconditioner.The analysis of the performances of the proposed iterative solver is deferred to the last example here proposed, in which a more complex setting is considered.

Test Problem 2 (TP2)
For this numerical example we consider a set of 19 inclusions of radius Ř = 10 −2 , whose centerlines intersect in 9 points.The resulting network is embedded in the same cube of edge l = 2 considered in Test Problem 1.We impose homogeneous Dirichlet boundary conditions on all the faces of the cube and at the dead ends of the network which intersect the cube at its top and bottom faces, as shown in Figure 6.Homogeneous Neumann boundary conditions are imposed at segment endpoints lying inside the cube.For what concerns problem coefficients we consider K = 1, f = 0 and Ki = 100, g i = 100, β i = 5e − 2, ∀i = 1, ..., 19.
The problem is solved both with the method proposed in this article and with a different approach in which no auxiliary variables are introduced.The  following 3D-1D coupled problem is derived from ( 7)-( 8 which, after discretization yields the following global system where the nomenclature is the same of Section 5 and the new matrix B is defined as We will refer to this method as "coupled" and we will use it as a comparison term for our approach, which, instead, will be labelled as OPT (optimization based).The results shown in the following are obtained by considering a 3D mesh, non conforming to the inclusions, with h = 0.083 and N = 3320 and a 1D mesh with δu = 1, corresponding to N = 234 DOFs.Parameters δ D = δ Σ = 0.5 are used in the OPT approach.Figure 6 reports the solutions obtained with both approaches on the network of segments, while Figure 7 on the right proposes a comparison of the solutions on two selected segments (marked in Figure 6), showing an almost perfect agreement.Figure 7 on the left, instead, reports

Conjugate gradient test (CGtest)
We now consider a more complex numerical example, characterized by the presence of multiple intersecting inclusions.The setting of this example might be considered as realistic of a living tissue with a network of vessels.The purpose of the present example is to test the performances of the proposed resolution strategy and preconditioner in a realistic setting.In particular we consider the domain of Figure 8, where 873 segments organized into two connected clusters are immersed in a cubic domain Ω of edge l = 2, as considered in the previous examples.On the faces of the cube we consider Neumann boundary conditions, namely K∇u • n = 2 • 10 −5 , with n denoting, in this case, the outward pointing unit normal vector to ∂Ω.At the inlets of the two networks, i.e. the segment endpoints lying on the face z = −1, we impose Dirichlet boundary conditions equal to 5 • 10 −3 , and at all other segment endpoints we consider homogeneous Neumann conditions.Problem data are as follows: The cluster of segments of this example is composed by many segments of small length, compared to domain size.Thus, a set of simulations is performed by varying both parameter h of the 3D mesh and parameters δ D and δ Σ of the 1D meshes, in order to control the meshsize independently, being, instead, δu = 1 fixed.The values of the parameters used in the simulations are reported in the first three columns of Table 7.3, whereas the fourth column reports the number of DOFs of the variables Ψ D and Ψ Σ , corresponding to the size of the system (39).The coarsest and the finest mesh combination considered are shown in  ).It can be seen that, in all cases, the number of iterations is small compared to the number of unknowns, and it only marginally grows as the stopping tolerance is reduced.The use of the preconditioner allows to further reduce the number of iteration, the obtained reduction ranging between 15% and 30%.It is to be highlighted that the proposed preconditioner can be obtained and applied at a very low computational cost, as it only involves the resolution of local 1D problems and can be computed and used in parallel.The effectiveness of the proposed resolution approach reflects the good conditioning of the obtained system, as it was pointed out in Test Problem 1.

Conclusions
A PDE constrained formulation for 3D-1D coupled problems with discontinuous solution at the interfaces has been derived and proposed.The approach is based on the introduction of unknown interface variables to decouple the 3D and 1D problems and on the minimization of a cost functional to enforce interface conditions.The problem is discretized resorting to standard finite elements on non conforming meshes independently set on each subdomain.Well posedness results for the discrete problem are obtained independently of the choice of the mesh parameters of the various domains.The proposed test on a problem with known analytical solution shows that optimal convergence trends of the error are obtained for both the 3D and 1D solution.Also the linear system corresponding to the application of the method appears to be well conditioned for a wide range of choices of the mesh parameters.The examples on more complex domains reveal the applicability of the method to realistic configurations and also the good performances of the proposed gradient-based solver.

Figure 1 :
Figure 1: TP1: left: view of the numerical solution inside the cube; right: solution obtained on the segment and on a section of the cube parallel to the z-axis and containing Λ. Parameters h = 0.13, δu = 1, δ D = δ Σ = 0.5.

Figure 2 :
Figure 2: TP1: trend of the L 2 and H 1 -norms of the relative errors under mesh refinement.On the left: error on the cube with respect to (41); on the right: error on the segment with respect to (42).Other parameters: δu = 1, δ D = δ Σ = 0.5.

Figure 2
Figure 2 displays the convergence trends for the above quantities against mesh refinement.Four meshes are considered, obtained by choosing h = 0.208, 0.131, 0.083, 0.052, which correspond to N = 257, 1026, 4155, 16545 DOFs and N = 15, 29, 57, 88 DOFs, respectively.The parameters δu = 1 and δ D = δ Σ = 0.5 are fixed for all cases, hence a 3D mesh refinement induces a refinement of all the 1D meshes.Optimal convergence rates are obtained for the considered indicators in relation to the regularity of the exact solution, as reported in the picture.Two additional error indicators are instead considered for the interface variables Ψ D and Ψ Σ as:

Figure 4 :Figure 5 :Figure 6 :
Figure 4: TP1: trend of the conditioning of the KKT system under the variation of the 1D mesh parameters.On the left variable δ D and different values of δu while δ Σ = 0.5.On the right variable δ Σ and δ D = 0.5.In both cases h = 0.083

Figure 7 :
Figure 7: TP2: On the left: solutions with OPT and coupled methods obtained inside the cube on three different planes parallel to the x − y plane and located at z = −0.5, z = 0 and z = 0.5.Solution amplified by a factor 50 with respect to domain size; on the right: comparison of the solution on two selected segments with the OPT and the coupled methods.

Figure 8 :
Figure 8: CGtest: Representation of the problem geometry.The blue spheres highlight the Dirichlet boundary conditions at the inlets of the two networks.

Figure 9 ,
Figure9, whereas the solution on the finest mesh is reported in Figure10, on the left for the 1D network and on the right for the 3D solution on three planes orthogonal to the z-axis.The remaining columns of Table7.3 report the number of iterations required by Algorithm 1 to solve system (39) up to a relative residual of 10 −6 or 10 −9 , without using a preconditioner (columns CG(10 −6 ) it and CG (10 −9 ) it ) or using preconditioner (40) (columns P CG (10 −6 ) it