EXTENDING BROYDEN’S METHOD TO INTERACTION PROBLEMS

. The solution of problems involving the interaction of diﬀerent systems is a domain of ongoing research, although often a good solver already exists for each system separately. In this paper we draw our ideas from one of the best known all-round quasi-Newton methods: Broyden’s rank-one update, which we extend to algorithms using 2 approximate Jacobians. A comparison is made with the iterative substructuring method and Aitken’s acceleration method. It is shown that a Broyden method using only a single approximate Jacobian performs best.


INTRODUCTION
Simulations of coupled systems are becoming ever more important.While powerful solvers often exist for problems in a single physical domain (e.g.structural problems, flow problems, ...), development of similar tools for multiphysics problems is still ongoing, and these can be broadly put in one of the following categories: • Monolithic or simultaneous solution: the whole problem is treated as a monolithic entity and solved with a specialized solver.• Partitioned solution: the physical components are treated as isolated entities that are solved separately.The relative merits of these methods are very problem dependent: the advantage of the monolithic approach is the enhanced stability.This comes at a cost, however, as specialized software has to be written for each type of interaction problem, which can result in very large systems.The partitioned approach allows for the use of available specialized solvers for each physical component (structure, fluid,...), on condition that the coupling effects can be treated efficiently.The latter is often no problem for loosely coupled problems, which can be solved by sequentially calling the solvers for each component.Strongly coupled problems, on the other hand, still pose a real challenge; the simplest solution method is the iterative substructuring method, but this often displays very poor stability [4].
In this paper we write the coupled problem as a root-finding problem which we solve with Broyden's quasi-Newton method and an extension of this method in which two approximate Jacobians are used, one for each physical problem.This results in a partitioned approach in which each of the physical problems is solved with a specialized code that we consider to be a black box and of which the Jacobian is unknown.We assume the dominant cost lies in the use of these black boxes.
This paper is organized as follows: in section 2 we describe the type of problems we are interested in; in section 3, resp.4, we describe two wellknown methods for the solution of these problems: the iterative substructuring method, resp.Aitken's ∆ 2 method; in section 5 we briefly describe Broyden's quasi-Newton method together with two proposed variations; in section 6 we test the various algorithms on two coupled problems.

THE GENERAL MODEL
We consider coupled problems that can be stated in the following form: where Each equation describes a physical problem that is spatially decomposed.We limit p and g to values on the interface between the two physical problems.In this way the physically decoupled nature of the problem is exploited, and the method can be regarded as a special case of heterogeneous domain decomposition methods [2].
We want to solve (1) with a partitioned method that keeps the available solvers for each individual equation (1a), resp.(1b), of which we lack knowledge of the Jacobian.(1) can also be written as the root-finding problem Depending on the algorithm we will either use the formulation of (1) or (2).We assume that F , S, H and K satisfy the following properties: Remark.The choice, in equation (2), between F (S(p)) = 0 or S(F (g)) = 0, will depend on the relative dimensions of p and g; we will choose the formulation such that the (square) Jacobian of K is of the lowest dimension.

ITERATIVE SUBSTRUCTURING METHOD
The simplest strongly coupled method to solve (1) is the iterative substructuring method; it proceeds as follows: Algorithm 3.1 (Iterative substructuring method (ISM)).1. Startup.Take an initial value p o .
Depending on the strength of the coupling between the two equations in (1) this method is not always stable [4].
Another method that is often used for the solution of ( 2) is Aitken's ∆ and •, • is the classical scalar product.)

BROYDEN'S METHOD AND VARIATIONS
A well-known solution method for the nonlinear problem in ( 2) is Newton's method.However, as we suppose we don't have access to the Jacobians of either F , S, H or K, we will resort to quasi-Newton methods, which we apply either to the single equation (2) or the system (1).
5.1.Basic Broyden method.Broyden's method [1] was developed for a single equation as in (2).The resulting algorithm can be written as follows.
o is needed.As we do not have any idea what the Jacobian of F or S is, we will set them equal to zero, such that we have K o = −I as a consequence.

Extension of Broyden's method.
We recall that K(p) = F (S(p)) − p and that as such We therefore propose to build approximate Jacobians F s and Ŝ s along the lines of the Broyden update: It is easy to see that this construction respects the following secant equations:1 This now allows us to write the following algorithm: Algorithm 5.2 (Broyden's method-2 Jacobians-Variant 1 (BrVar1)).
Remark.While Broyden's original update also respects (4c), it can be shown that the approximate Jacobians differ between algorithms 5.1 and 5.2.
In [5] another approach was proposed to couple (1) with two approximate Jacobians, which we now apply to Broyden's method.In this approach we look at each equation in (1) separately and write a Taylor expansion in which we insert the approximate Jacobians: We require the solution of the coupled problem to satisfy (1) at iteration s + 1; it follows that Solving (5) for g and p gives We assign the value of p in equation (6a) to p s+1 , and use it as an input for S, We can then update the Jacobian ( Ŝ ) in equation (6b) and assign the value of g to g s+1 : ) .Combining this with Broyden's method results in the following algorithm: Algorithm 5.3 (Broyden's method-2 Jacobians-Variant 2 (BrVar2)).
Remark.It is obvious from the way we construct F s and Ŝ s and equations (5a) and (5b) that they respect the secant equations (4a) and (4b).For reasons of comparison with the other algorithms, we start from F o = Ŝ o = 0.

TEST-CASES
We propose two test-cases: the heat equation and a fluid-structure interaction problem.The convergence requirement will be defined based on a relative reduction of the residual, which is defined by K(ps) K(po) ≤ 10 −5 for algorithms 4.1, 5.1 and 5.2, and on F (gs)−ps F (go)−po ≤ 10 −5 for algorithms 3.1 and 5.3, as K(p s ) is not explicitly computed in these.The performance parameter is the number of times the function F (or S) is called, as we have assumed this is the dominant cost of the algorithm; we call this value FC.We break off the algorithm if FC > 100 if it has not converged at that point.
where T is temperature, ρ is density, C is heat capacity, k is thermal conductivity.We assume we only have a linear solver to solve (7) with an imposed value of ρ, C and k that may vary in space.
If we want to take dependency of ρ, C and k on T into account, we need to do this outside the solver for (7).They can be obtained by the following interpolation polynomials (based on values between 0 and 300 • C in [6]): We use a finite differencing scheme on N equidistant nodes for the solver of (7), with spacing ∆x, and an implicit time discretization with time step ∆t.This gives the N linear equations (i = 1, . . ., N ): We will use g to denote the ensemble of discretized constants ρ i , C i and k i , which will be defined by a node-to-node relation with T i according to (8).
We only consider one time-step and solve (9) for the temperature at the next time-step.We write this as F (g s ) = T s+1 .The parameters are computed for a given temperature, which is written as S(T s ) = g s .
As an initial condition we take the values at time-level n, which is a uniform temperature of T o = 150K, the right boundary condition is a temperature of 150K, while the left boundary condition is a function of time, defined as 2 sin πt 10 .We use an initial relaxation factor ω = 0.1.
6.1.2.Results and discussion.Various combinations of the ratio ∆t ∆x 2 were tried for N = 100 and N = 1000 (tables 1 and 2).For small values of ∆t ∆x 2 the three Broyden algorithms and the ISM performed equally well; Aitken's method trailed behind.As the ratio ∆t ∆x 2 grows, the problem becomes harder to solve and the number of iterations becomes larger.The performance of the ISM initially degrades faster than the Broyden methods, which perform equally well.For N = 1000 and ∆t ∆x 2 = 10 9 the ISM method catches up with the basic Broyden method, while the other three methods lag behind.For even higher values of ∆t ∆x 2 convergence could not be obtained.It is thus seen that for this mildly non-linear problem there is not much difference between the methods, but the basic Broyden method is always the best.This is perhaps surprising as it only uses one Jacobian and thus less information than the Broyden variants.6.2.Fluid-structure interaction problem.In this problem both physical problems are vector-based and non-linear.6.2.1.The model.We consider one-dimensional unsteady flow in an elastic tube (figure 1).(For a two-dimensional simulation, see [5]).The fluid is incompressible and inviscid and gravity is neglected.The governing equations are the conservation of mass and momentum which can be written in conservative to which suitable boundary conditions need to be added (velocity at the inlet and pressure at the outlet).g represents the cross-sectional area of the tube, u the velocity along the axis of the tube, ρ is the density of the fluid and p the static pressure.
We introduce the kinematic pressure p = p /ρ.If the elastic wall of the tube has a constitutive law of the form g = g(p(x, t)), with its inertia neglected, then (10) can be rewritten in the following form: The wave speed c is defined by For the constitutive law of the structure we take the non-linear relation defined by: where c mk is the Moens-Korteweg wave speed c mk = Eh 2ρro , E Young's modulus, h the thickness of the tube wall and r o a transversal reference length of the tube.The length of the tube is L.
The velocity is imposed at the inlet as where u o is the initial velocity and T the period (= 1).At the outlet, a non-reflecting boundary condition is imposed: We use a one-dimensional mesh with N nodes of length ∆x.Fluid velocity, pressure and diameter of the tube are stored at the cell centers.Central discretization is used for all terms in (11), except for the convective term in (11b) where a first-order upwind scheme is used; pressure stabilization is added [5].Pressure and velocity at the inlet, resp.outlet are linearly interpolated from nearby mesh-points.For the time-discretization we use backward Euler with time-step ∆t.
As in §6.1, we only consider one time-step and use those values as initial values for the iteration.We solve the discretized fluid problem for the pressure at the next time-step and a fixed geometry, which we write as F (g s ) = p s+1 .For the structural problem we compute g for a given pressure, which we write as S(p s ) = g s .
The initial relaxation parameter is ω = 10 −2 , except for (κ, τ ) = (10, 10 −2 ) and (10, 10 −3 ), where it is 10 −4 , resp. 10 −5 .6.2.2.Results and discussion.Various combinations of parameters were tried for N = 100 and N = 1000.When the governing equations are non-dimensionalized it is seen that they are characterized by two non-dimensional parameters κ = co uo , being a dimensionless structural stiffness, and τ = uo∆t L , being a dimensionless time-step; c o , u o , are the linearizing, i.e. reference, values of the wave speed, resp.fluid velocity.Various combinations of these were tested (tables 3 and 4).It is seen that a more flexible tube (lower value of κ) and a smaller time-step (lower value of τ ) represent a more difficult problem to solve.The three Broyden methods steadily outperform the ISM and Aitken's method but share similar performance for the "easier" test-cases.When the problem becomes harder, the basic Broyden method appears the best, followed by the second variant of Broyden's method.Again this is somewhat surprising as all the information of the two systems S and F is lumped into a single approximate Jacobian K s .The basic Broyden method was tested on interaction problems and compared with two proposed variants as well as the iterative substructuring method and the Aitken ∆ 2 method.For weakly non-linear problems and/or problems requiring few iterations to converge there was little to chose from between the three Broyden methods, although they were better than the remaining two.For more difficult problems not only did the three Broyden-based methods outperform the other two, but it appeared that the basic Broyden method came out best.This is surprising as, compared with the two variants, it only uses a single Jacobian to capture all the information that is available of the two constituent systems, while the other two use two approximate Jacobians, one for each system.

( 1 )
F , resp S, H, is twice continuously differentiable in an open set D F , resp D S , D H . (2) There is a p o such that the set C(p o ) = {p; K(p) ≤ K(p o ) } is compact.(3) K(p) = 0 has a single solution p * in C(p o ) (4) (K (p)) −1 exists and is continuous in an open set containing C(p o ).

6. 1 .
Heat equation.6.1.1.The model.The heat equation in one dimension without heat source is given by:

Table 1 .
FC required for convergence for the heat equation.N = 100.