Transcript A Matrix Free Newton /Krylov Method For Coupling Complex Multi-Physics Subsystems Yunlin Xu School of Nuclear Engineering Purdue University October 23, 2006

A Matrix Free Newton /Krylov Method
For Coupling Complex Multi-Physics Subsystems
Yunlin Xu
School of Nuclear Engineering
Purdue University
October 23, 2006
1
Content
 Introduction
 MFNK and Optimal Perturbation Size Fixed Point Iteration (FPI)
for coupling subsystems
–
–
–
–
A Matrix Free Newton/Krylov method based on FPI
Local Convergence analysis of MFNK
Truncation and Round-off Error
Estimation of Optimal Finite Difference Perturbation
 Global Convergence strategies
– Line search
– Model trust region
 Numerical Examples
 Summary
2
Features of Multi-Physics Subsystems
 Multiple nonlinear subsystems are coupled together:
The solution of each subsystem depends on some
external variables which come from the other system
fi ( xi , yi )  0
xi  ni
: internal variables
yi   i ( x1, x2 ..., xi 1, xi 1..., x p )
: external variables
 Each subsystem can be solved with reliable methods
as long as they remain decoupled
3
Two General Approaches for Coupling Subsystems
 Analytic Approach: reformulate the coupled system into
a larger system of equations
– Standard Newton-type methods can be applied
 Synthetic Approach: combine the subsystem solvers for
the coupled system
– utilize the well-tested and reliable solutions of each of the
subsystems because:
• It may be too expensive to reformulate the coupled system and
forego the significant investment in developing reliable solvers for
each of the subsystems.
• One of the subsystems may be solved using commercial software
that prevents access to the source code which makes it
impossible to reformulate the coupled equations for the analytic
approach
4
Coupled Subsystem Example:
Nuclear Reactor Simulation
5
Time Advancement:
Marching vs Nest Scheme
step n+1
step n
TRAC-E
 Marching
PARCS
– The time steps must be kept small for accuracy and stability concerns
step n
 Nested
step n+1
TRAC-E
N
Converged
?
N
Y
N
Converged
?
Y
Converged
?
Y
PARCS
– Computational cost for each time step increased
– Numerical Stability and accuracy can be improved
– Time step size may be extended
6
Ringhals BWR Stability
 48 hours on 2 GHz machine for initialization!
1
2
3
4
5
6
7
8
9
1
10
11
12
S8 (12)
2
15
16
17
18
154
154
154
154
19
20
21
22
23
24
25
26
27
854
854
853
153
153
153
154
154
854
854
853
853
151
152
152
153
154
154
154
154
854
854
853
853
852
852
852
153
151
152
152
152
153
153
154
154
854
854
853
852
852
853
851
852
109
110
151
153
152
153
152
152
154
154
854
854
853
852
852
853
852
853
801
110
109
153
151
153
152
153
151
152
254
254
754
854
852
853
852
853
851
853
805
806
103
104
111
112
151
153
155
156
252
252
254
254
754
754
752
852
851
853
852
853
803
806
805
104
101
112
111
152
151
156
155
252
251
253
254
254
5
6
S7 (20)
28
30
S2 (12)
9
754
754
752
752
755
756
851
853
807
808
801
804
105
106
103
103
113
113
251
253
252
253
251
252
254
254
10
754
753
752
751
756
755
852
801
808
807
804
803
106
105
103
101
113
113
203
251
253
252
252
252
253
254
11
754
753
752
753
751
752
704
704
801
803
802
804
101
104
107
108
203
204
205
206
251
253
252
253
253
254
754
754
752
753
752
753
751
703
704
803
803
803
802
103
103
108
107
204
201
206
205
203
251
253
252
252
12
29
S1 (19)
854
4
8
14
854
854
3
7
13
854
254
254
13
754
754
753
752
751
753
712
713
711
711
708
709
801
804
102
103
201
203
202
203
203
204
207
208
251
253
252
252
254
254
14
754
753
752
752
752
751
713
712
711
710
709
708
804
803
104
102
204
202
204
202
204
201
208
207
203
251
253
252
253
254
15
754
753
751
752
714
714
701
703
706
707
702
703
705
705
901
902
305
305
301
303
306
307
303
303
314
314
351
353
353
354
16
754
753
753
751
714
714
703
703
707
706
703
701
705
705
902
901
305
305
303
302
307
306
303
301
314
314
352
351
353
354
17
654
653
652
653
651
603
607
608
601
604
602
604
602
604
502
504
403
404
308
309
310
311
312
313
351
352
352
352
353
354
18
654
654
652
652
653
651
608
607
604
603
603
602
603
601
503
502
404
401
309
308
311
311
313
312
353
351
352
353
354
354
654
354
19
654
652
652
653
651
603
605
606
601
604
507
508
503
503
402
403
403
403
304
303
351
353
352
353
352
354
20
654
653
653
652
653
651
606
605
604
603
508
507
504
501
404
402
403
401
304
304
352
351
353
352
353
354
21
654
653
652
652
652
653
651
603
513
513
501
503
505
506
403
404
407
408
401
452
355
356
351
352
353
354
22
654
654
652
651
653
652
653
651
513
513
503
503
506
505
404
401
408
407
453
451
356
355
352
352
354
354
654
654
653
651
652
555
556
551
552
511
512
501
504
405
406
403
453
452
453
451
452
352
354
354
654
654
652
652
556
555
553
551
512
511
504
503
406
405
453
451
453
452
453
452
454
354
654
654
552
551
553
552
553
551
553
509
510
401
453
452
453
452
452
453
454
454
554
554
552
552
553
552
553
551
510
509
452
451
453
452
452
453
454
454
554
554
553
553
552
552
552
551
553
452
452
452
453
453
454
454
554
554
554
554
553
552
552
551
453
453
454
454
454
454
554
554
553
553
553
453
454
454
554
554
554
554
454
454
23
24
25
26
27
28
S6 (12)
29
30
S5 (19)
128
Chans
S3 (20)
S4 (12)
7
Synthetic Approaches
 Nested Iteration
– Subsystems are chained in block Gauss-Seidel or block Jacobi iteration
– Convergence is not guaranteed.
 Matrix Free Newton/Krylov Method
– Approximate Mat-Vec by quotient:
F ' ( x k )v j 
– The system Jacobian is not constructed
F ( x k  v j )  F ( x k )

– Local Convergence guaranteed
 Problems with Direct Application of MFNK for Coupling of Subsystems
– Solvers for the subsystems are not fully utilized
– Difficult to find a good preconditioner for MFNK
– In some cases, it is not possible to obtain residuals for a subsystem if the
solver of subsystem is commercial software which can be used only as a
“black box”.
8
Objectives of Research
 Propose a general approach to implement efficient
matrix free Newton/Krylov methods for coupling
complex subsystems with their respective solvers
 Identify and address specific issues which arise in
implementing MFNK for practical applications
– Local convergence analysis of the matrix free Newton/Krylov
method
– Optimal perturbation size for the finite difference
approximation in MFNK
– Globally convergent strategies
9
Fixed Point Iteration for Coupling Subsystems
 Block Iteration
xik 1  i ( xk )  iti ( xik , yik ), 1  i  p,
k N
it denotes ti iterations, it 1 ( xik , yik )  i (it ( xik , yik ), yik )
i
i
i
 Block Iteration for coupling subsystems are fixed point
iterations
x k 1  ( x k ), x  ( x1 , x2 ,...x p ),   (1 , 2 ,... p )
x  n ,
p
 : n  n , n   n i
j 1
 The condition for convergence of FPI is ||Φ(x*)||<1.
 If ||Φ(x*)||>1, then the FPI may diverge
10
Matrix Free Newton Krylov Method
Based on a Fixed Point iteration
 Define a nonlinear system: F(x)= Φ(x) -x =0
 The solution of this system is the fixed point of function Φ(x),
which is also the solution of original coupled nonlinear system
 MFNK algorithm:
At kth Newton step, do
1) Find s k satisfies
k
F '( xk )sk  F ( xk )  rItk  rFD
by a Krylov iteration, in which mat-vec, F ' ( x k )v j , be approximated by:
F ' ( x )v j 
k
F ( x k  v j )  F ( x k )


( x k  v j )  ( x k )

 vj,
where rItk is the iterative residual,
k
is the residual due to the finite difference approximation.
rFD
2) Set
x k 1  x k  s k
11
Local Convergence of INM
 Inexact Newton Method (INM)
x
k 1
 x   F '( x ) 
k
k
 Local Convergence of INM
1
 F ( x )  r 
k
k
p
If r k   F ( x k ) ,   1 or r k   F ( x k ) , p>1, Then
p
 1

|| x k 1  x* ||  x k  x *  2 
 || F ' ( x*) || x k  x *
 4

2
p
 The convergence of INM depends on the inner residual,
p
assume
r k   F (xk )
– If p2, the INM has local q-quadratic convergence.
– If 1<p<2, INM converges with q-order at least p.
– If p=1 and ,   1 the INM has local q-linear convergence.
12
Local Convergence of MFNK
x k 1  x k   F '( x k ) 
 MFNK is an INM
The inner residual consists with:
– iterative residual
– finite difference residual
Em  d1  dm v1  vm 
T
1
 F ( x )  r
k
k
It
k
 rFD

rItk
k
rFD
 Em s k
d
F ( x   v)  F ( x)

 F '( x)v
 There are two conflicting sources of error in finite difference:
– Truncation error
F ( x   v)  F ( x)   F '( x)v
2
etr 
  v

– Round off error
erd 
F ( x   v)  F ( x   v)  F ( x)  F ( x)


1

 F ( x)   F ( x   v)   M x    M v
13
Local Convergence of MFNK (cont.)
 The optimal  should balance the round-off error and
truncation error
F ' ( x)
4 F ( x)
1
 opt 

k
rFD
 opt F ( x k )  rItk
v



x 

opt  4 opt  m F ' ( x k ) 1
 In theory, MFNK has local q-linear convergence, if
opt  1
 In practice, MFNK can achieve q-quadratic convergence, if
opt   || F ( x ) ||
k
14
Optimal vs Empirical Perturbation Size
 The norm of the Jacobian and  can be estimated with
information provided by the MFNK algorithm:
F ( x  v )  F ( x )
~
M  max
1 j  m
v 2
~opt
1

v
2
 F ' ( x)
2 F ( x k 1 )
F ' ( x k  0.5s k )s k  F ' ( x k )s k
~
2
 

2
k 2
s
0.5 s k
2
2
~
2 F ( x)
M
 ~ x 

~
or
~
opt
2

2
1

v
~
M
2
 ~ x  ~ x


 An empirical prescription was proposed attempt to balance the
truncation and round-off errors (Dennis)
 emp  
max  x , typ x 
v
15
Global Convergence Strategies
 Solution x* of system of nonlinear equation: F(x)=0
is also the global minimizer of optimization problem:
min f ( x) 
xR n
1
F ( x)
2
2
2
 Newton step sN is the step from current solution to
global minimizer of model problem:
c
m
(
x
min  s) 
sRn
1
F ( x c )  F '( x c )s
2
2
2
 f(xc+sN) may be larger than f(xc), due to big step sN such
that m(xc+sN) is no longer a good approximation of f(xc+sN).
In this case, we need globally convergent strategy to force
f(xc+sN)<f(xc)
16
Descent Direction
 Newton step is descent direction of both objective
function and its model:
f ( xc )  m( xc )  F '( xc )T F ( xc )
f ( xc )T s N  F ( xc )T F '( xc )s N  F ( xc )T F ( xc )  0
 For any descent direction pk, there exist λ satisfies: (1)
f ( xk   k pk )  f ( xk )   k f ( xk )T pk
f ( xk   k pk )  f ( xk )   k f ( xk )T pk
0    1
α-condition
β-condition
 A sequence {xk} generated by xk+1=xk+λkpk satisfying
previous condition will converge to a minimizer of f(x). (2)
(1),(2) proofs can be found in Dennis & Schnabel’s book
17
Line Search
 Take MF Newton step as descent direction, and select λ
to minimize a model of following function
ˆf ()  f ( xk  s M )
 Quadratic model
mˆ q ( )  fˆ (0)  fˆ '(0)   fˆ (1 )  fˆ (0)  fˆ '(0)1   2 / 12
 λ predicted from quadratic model
q 
1
2  2 fˆ (1 )  fˆ (0)

   fˆ '(0) 
1
18
Information Required
in Quadratic Model
 Two function values:
fˆ (0)  f ( xk )
fˆ (1 )  f ( xk  1s M )
 One Gradient
fˆ '(0)  f ( x k )T s M  F ( x k )T F '( x k ) s M
fˆ '(0)  F ( x k )T ( F ( x k )  r )  2 fˆ (0)  F ( x k )T r
 Approximations for Gradient in MFNK
or
fˆ '(0)  2 fˆ (0)  F ( xk )T rit
fˆ '(0)  2 fˆ (0)
19
Model Trust Region
 Minimize model function in neighborhood, trust region
Min m( x c  s ) 
2
1
1
F ( x c )  F '( x c ) s  f ( x c )  m( x c ) s  sT  2 m( x c ) s
2
2
2
s 2  c
subject to
xc+s(c)
xN
sN
xc
c
20
Double Dogleg Curve
 Approximate optimal path with double dogleg curve
C.P.
xN
sN
N
xc
c
 Step along double dogleg curve.
 C .P

C .P
s
,
if



C
.
P



 C .P    C .P
s ( )   s  N
( N s N  s C .P ), if  C .P     N
C .P
 

sN ,
if    N



21
Cauchy Point
 Cauchy Point is minimizer in steepest descent direction
s C .P  
m ( x c )
2
2
m( x c )T  2 m( x c )m( x c )
JTF
m ( x c )  
T
JJ F
2
2
2
JTF
2
 Projection of Step for Cauchy Point on Krylov subspace
(Brown & Saad)
 JVm 
T
sˆC .P  Vm
2
F
2
JVm  JVm  F
JVm  Vm1Hm ,
T
JVm 
2 
T
F  Vm
dm
2
2
H mdm
2
d m   JVm  F   H mT e1 ,
T
2
2
d m  Vm
dm
2
2
Rm d m
2
dm
2
H m  QRm
22
Example Problem I: Polynomials
 Two dimensional second order polynomials
 f1 ( x1 , x2 )  a( x1  1)  b( x2  1)  c( x12  1)
F ( x1 , x2 )  
2
f
(
x
,
x
)

b
(
x

1)

a
(
x

1)

c
(
x
 2 1 2
1
2
2  1)
 Solution
( x1* , x2* )  (1,1)
 Jacobian
b 
 a  2cx1
F '( x1 , x2 )  

b
a

2
cx

2
  F '( x1* , x2* ) 2  max  a  2c  b , a  2c  b 
 1  F '( x1* , x2* )1
1
2
 min  a  2c  b , a  2c  b 
   2  F '( x1* , x2* )   
 Nonlinear level
 2 c max{ x , x } 

1
2 
  max 
2 c
2
2
x1 , x2

  x1    x2  

23
PLY 1
Truncation Error Dominated Case
Input values
a  1  2  10
10
b  0.5
c  1010
20
15
10
5
0
Key parameters
  1.5,   2,
-5
-10
  3,   2  1010 ,
-15
 opt  1.53  10 13
-20
-20
-15
-10
-5
0
5
10
15
20
24
PLY 1 Errors
1.E+02
1.E+00
1.E-02
0
10
20
30
50
Optimal Step Size
Empirical Step Size
1.E-04
L2 norm of Error
40
1.E-06
1.E-08
1.E-10
1.E-12
1.E-14
1.E-16
1.E-18
1.E-20
Newton Step
25
PLY 1 Step Sizes
Optimal
Empirical
-6
-6
10
Computed from Jacobian
Estimated by MF
Computed from Jacobian
Empirical
finite difference Stepsize : sigma
Optimal finite difference Stepsize : sigma
10
-8
10
-10
10
-12
10
-14
10
0
-8
10
-10
10
-12
10
-14
10
20
30
Newton Step
40
50
10
0
10
20
30
Newton Step
40
50
26
PLY 1 Step Size Parametric
Relative error or residaul after
1 step Newton iteration
1.E+01
1.E+00
1.E-16
1.E-14
1.E-12
1.E-10
1.E-08
1.E-01
1.E-02
1.E-03
Error
Outer residual
Inner residual
1.E-04
Finite Difference stepsize
27
PLY 2
Round Off Error Dominated Case
Input values
20
a  108
15
b  108  0.01
10
8
5
c  10
0
-5
Key parameters
-10
  2 108 ,   100,
-15
  2 1010 ,   2 108 ,
 opt  1.772
-20
-20
-15
-10
-5
0
5
10
15
28
20
L2 norm of Error
PLY 2 Errors
1.E+14
1.E+12
1.E+10
1.E+08
1.E+06
1.E+04
1.E+02
1.E+00
1.E-02 0
1.E-04
1.E-06
1.E-08
1.E-10
1.E-12
1.E-14
1.E-16
1.E-18
1.E-20
Optimal Step Size
Empirical Step Size
10
20
30
40
50
Newton Step
29
PLY 2 Stepsizes
Optimal
Empirical
2
10
10
Finite difference Stepsize : sigma
Optimal finite difference Stepsize : sigma
10
0
10
-2
10
-4
10
Computed from Jacobian
Estimated by MF
-6
10
-8
10
1
5
10
0
10
-5
10
Computed from Jacobian
Empirical
-10
1.5
2
2.5
3
3.5
Newton Step
4
4.5
5
10
0
10
20
30
Newton Step
40
50
30
Relative error or residaul after
1 step Newton iteration
PLY 2 Stepsize Parametric
1.E+01
1.E+00
1.E-01
1.E-12
1.E-02
1.E-03
1.E-04
1.E-05
1.E-06
1.E-07
1.E-08
1.E-09
1.E-10
1.E-11
1.E-12
1.E-13
1.E-14
1.E-15
1.E-16
1.E-17
1.E-08
1.E-04
1.E+00
1.E+04
1.E+08
 opt  1.772
Error
Outer residual
Inner residual
Finite Difference stepsize
31
Numerical Examples:
Navier-Stokes-Like Problem (Goyon)
 PDE
 
 1
  
 u  (u  )u  (  u )u  q ( x, y )
2
2
Diffusion
Convection Non-physical Force function
 Boundary Condition

u  (u, v)  0 on ,
 Force function
  [0,1]  [0,1]

q ( x, y)  (q, s)
q  s   sin(2x) sin(2y)[8  3 sin 2 ( x  y)]
 Goyon, Precoditioned Newton Methods using Incremental Unknowns Methods
for Resolution of a Steady-State Navier-Stokes-Like Problem. Applied
Mathematic and Computation, 87(1997), pp. 289-311.
32
The Finite difference Equations
 2 2 vi , j 1  vi , j 1 
3u 
3u 


 
ui , j      i , j ui 1, j      i , j ui 1, j

 h 2 4h 
 h 2 4h 
 h2 h2

4hy
y
x 
x 
 x
 x
 x

 
 
vi , j 
vi , j 



ui , j 1  qi , j  0
  2
ui , j 1   2 
 h



 y 2hy 
 hy 2hy 
 2 2 ui 1, j  ui 1, j



 h2 h2
4hx
y
 x
 
3vi , j
  2 
 h
4hy
y

u 
w 
v 

u 
u


vi , j      i , j vi 1, j      i , j
 h 2 2h 
 h 2 2h

x
x
x
x






3v
vi , j 1      i , j

 h 2 4h
y
y



vi 1, j


vi , j 1  si , j  0


0  u  bu  U 
 Au (u, v)
F ( w)  A( w) w  b  
 v   b   V   0
0
A
(
u
,
v
)
v

   v   
33
Structure of the Matrices
Jacobian
Au

 Au  u u
F ' ( x)  
Av

v
 u
Au
u
v
A
Av  v
v



v

Diag block
Au

 Au  u u
B( x)  

0



A 
Av  v v
v 
0
34
NSL1: 50X50 meshes, =0.0015
 Solving (u,v) as One Nonlinear System, w0=(1-810-3) w*
Newton
I. N.(6)
MFNK
I. N.
MFNK
Preconditioner
N.A.
I
I
B
B
NNI(1)
3
18
20
3
3
NLI(2)
N.A.
3600
3991
93
94
NB(3)
0
0
0
0
0
NFE(4)
4
5
4012
4
98
NMV(5)
N.A.
3600
0
93
0
(1)NNI: Number of nonlinear iterations
(2) NLI: Number of cumulated Krylov iteration for Newton linear systems
(3) NB: Number of backtracking
(4) NFE: Number of function evaluation
(5) NMV: Number of Mat-Vecs in Krylov iteration
(6) Inexact Newton with GMRES as its linear system solver
35
NSL1: Newton Iterative error and residual
Error
Residual
1.E+00
5
10
15
20
L2 norm of error of iterative solution
1.E-02
Newton
IN(I)
MFNK(I)
IN(B)
MFNK(B)
1.E-03
1.E-04
1.E-05
1.E-06
1.E-07
1.E-08
1.E-09
1.E-10
1.E-11
1.E-01 0
L2 norm of nonlear system residual
1.E-01
1.E+00
0
5
10
15
20
1.E-02
Newton
IN(I)
MFNK(I)
IN(B)
MFNK(B)
1.E-03
1.E-04
1.E-05
1.E-06
1.E-07
1.E-08
1.E-09
1.E-10
1.E-11
1.E-12
1.E-12
1.E-13
1.E-13
Newton Step
Newton Step
36
NSL1: Coupling Subsystem
 Solving u and v as two subsystem, and coupled by
FPI or MFNK
FPI
Steps/tolerance
FPI
MFNK
schemes for Subsystems NNI NB(1) time(s) (2)
NNI NB NFE time(s)
1
60
/ 0
386
460
46.64
7
0
194
24.05
2
70
/ 0
388
2232
61.58
4
0
112
17.90
3
80
/ 0
450
2768
89.65
4
0
111
21.45
4
70
/ 10-4
358
2030
44.05
5
0
135
19.12
5
70
/ 10-5
547
3509
78.83
4
0
113
17.89
6
Direct solver
186
484
-(3)
4
0
135
- (3)
(1) Maximum 8 times halving the increment of solution for each Fixed point iteration
(2) Run on Purdue SP2 cluster, the CPUs are 375 MHz PWR3.
(3) Run on PC with Matlab linear system direct solver.
37
NSL1: Global convergence
w0=(1-810-3) w*
Global strategy
NO GS
LS
MTR
NNI
500
28
23
NLI
3534667
1875
1399
0
61
33
50501
1965
1456
NMV
3534667
135029
99993
Time(s)
7919.09
308.21
221.08
NB
NFE
38
L2 norm of Residual
NSL1 Residual
1.E+02
1.E+01
1.E+00
1.E-01 0
1.E-02
1.E-03
1.E-04
1.E-05
1.E-06
1.E-07
1.E-08
1.E-09
1.E-10
1.E-11
1.E-12
1.E-13
5
10
15
20
25
30
NO GS
LS
MTR
Newton Step
39
NSL1: First Backtracking
The first backtracking occurred after the first
Newton iteration.
The L2 norm of residual Lambda
before
1.65150
No GS
2.42510
LS
1.07990
0.316829
MTR
1.03694
0.316829
40
NSL2: 1000X1000 meshes
Test
NSL2a
NSL2a
-3
cases
(=10 ) (=10-3)
Scheme
FPI(1)
MFR(2)
(5)
(5)
NNI
500
500
NLI
0
25000
(6)
(7)
NB
4964
1273
NFE
5465
26774
NMV
25000
0
Time(s)
100422
69292
(1) FPI with 50 steps for solving
NSL2a
NSL2b
-3
(=10 )
(=510-4)
Opt(3) Emp(4) Opt
Emp
9
165
0
175
3500
9211
100
1945
261
2307
46140
65690
11
212
6
230
4600
12423
100
1877
271
2249
44980
100534
NSL2c
(=210-4)
Opt
Emp
11
420
8
440
9000
24731
100
2000
253
2354
47080
96791
linear system of subsystems
(2)
MFNK based on original nonlinear system, 50 Krylov steps for each Newton step
(3)
MFNK based on FPI, with optimal perturbation.
(4)
MFNK based on FPI, with empirical perturbation.
(5)
500 is limit of nonlinear iterations, NNI=500 indicates not converge
(6)
At most 10 times halving for each FPI
(7)
At most 3 backtrackings for each Newton step
41
NSL2: Residuals for FPI and MFR
FPI
MFR
9
10
9
8
7
L2 norm of Residual
L2 norm of Residual
8
NSL2a
6
NSL2b
5
NSL2c
4
3
NSL2a
NSL2b
NSL2c
8
7
7
6
2
6
1
0
5
0
10
20
30
FPI Step
40
50
0
10
20
30
40
50
Newton Step
42
NSL 2:
Optimal vs Empirical Finite Difference
Optimal
Empirical
1.E+01
1.E+01
1.E+00
5
10
1.E-01
15
20
NSL2a
NSL2b
NSL2c
1.E-02
1.E-03
1.E+00
L2 norm of error
L2 norm of error
0
0
20
40
60
80
100
NSL2a
NSL2b
NSL2c
1.E-01
1.E-02
1.E-04
1.E-05
1.E-03
Newton Step
Newton Step
43
NSL2 Step Sizes
1.E-04
Estimated Optimal Finite Difference Size
0
5
10
15
20
1.E-05
Empirical finite
difference size for NSL2
1.E-06
1.E-07
1.E-08
NSL2a
NSL2b
NSL2c
1.E-09
Newton Step
44
NSL2 Step Size parametric
L2-norm of error after 1 Newton step
5.E-05
4.E-05
Empirical
3.E-05
2.E-05
1.E-05
Estimated Optimal
0.E+00
1.E-12 1.E-11 1.E-10 1.E-09 1.E-08 1.E-07 1.E-06 1.E-05 1.E-04
Finite Difference Size
45
Summary
 A general approach, MFNK, was presented here for coupling
subsystems with respective solvers.
 Based on any FPI, a corresponding MFNK method can be constructed.
 MFNK provides a more efficient method than FPI for coupling
subsystems.
 MFNK can converge for several cases in which the corresponding FPI
diverges.
 Locally, MFNK converges at least q-linearly and in many cases qquadratically.
 A more sophisticated FPI scheme provides a more efficient nonlinear
system for the corresponding MFNK.
46
Summary (Cont.)
 A method was proposed to estimate the optimal perturbation size for
matrix free Newton/Krylov methods.
 The method was based on an analysis of the truncation error and the
round-off error introduced by the finite difference approximation.
 The optimal perturbation size can be accurately estimated in the MFNK
algorithm with almost no additional computational cost.
 Numerical examples shows that the optimum perturbation size leads to
improved convergence of the MFNK method compared to the
perturbation determined by empirical formulas.
47
Summary (Cont.)
 Line Search and Model Trust Region, were implemented within
the framework of MFNK
– For the line search method, a quadratic or higher order model was used with
an approximation for the gradient
– For the model trust region strategy, a double dogleg approach was
implemented using the projection of a Newton step and Cauchy point within
the Krylov subspace
– the model trust region strategy showed better local performance than the line
search strategy
 Peer-to-Peer parallel MFNK algorithm was implemented
48