岩下 武史
北海道大学情報基盤センター
Takeshi Iwashita

Hokkaido University, Information Initiative Center
iwashita .at. iic.hokudai.ac.jp

1

Iwashita lab.

 Part I 電磁場有限要素解析の定式化，解くべき連立

一次方程式の性質
 Part II 誤差修正のフレームワーク（EEC, IEC）
 Part III 時間周期電磁場解析の高速化（TP-EEC）
 Part IV 非正則な係数行列を持つ連立一次方程式の求解高速化（Folded preconditioning）

Iwashita lab.

Time domain
Linear
Non-Linear
Frequency
domainStaticLowfrequency
50～60Hz(quasi-static)

Highfrequency
MHz, GHz(wave)

FEM
BEM

FEM
BEM

FDTD
FEM
BEM

FEM
BEM

FEM
BEM

Iwashita lab.

Maxwell equations
(1)
(2)
(3)
(4)

Constitutive equations(5)(6)

Ohm’s law

H: Magnetic field
E: Electric field
B: Magnetic flux density
D: Electric flux density
ρ: Electric charge density
J: Current density: Permeability: Permittivity: Conductivity(7)
Iwashita lab.

 Mostly linear analyses

 Time-domain analyses◦ The most popular method is FDTD (finite differential
time domain) method.

 Frequency domain analyses◦ FEM (finite element method) is often used.
◦ Basic equation: external current density◦ The linear system derived from FEM is ill-conditioned
Iwashita lab.

The displacement current is ignored.
 Non-linear problems are solved in practicalsimulations for motors, transformers, etc.
 Potential formulation is often used.


◦
◦

◦


region.
,, where A is the vector potential.
, where is the eddy current in the conductive, whereis the electric scalar potential.

Gauge condition should be payed attention to.

◦ For fixed B, A is not uniquely determined. (in non-conductive
region)◦ For fixed B and E, a set of (A, ) is not uniquely determined.
(in conductive region)
Iwashita lab.

(8)




·

0

(9)

·
0 must be satisfied.
 Without a gauge condition, multiple solutions for(A, ) exist.
 Gauge condition◦ ·0 Coulomb gauge◦ =0 (in conductive region)◦ Tree co-tree gauge (after discretization) (in non-conductive
region)Iwashita lab.

 The edge element is generally used for themagnetic vector potential.
 The unknown is set on each edge of an mesh.
◦ The unknown is given by the integral of the fieldvalue along the edge.

 The vector basis function

is used.

First order tetrahedral elementltmis the fist order scalar basis function.
Discretization of A:
Iwashita lab.

 A linear system of equations is derived fromthe weak formulation of the basic equations
and the discretization. Galerkin method isused (test function = basis function).
· grad

0
1

·

·

(10)
grad

·

grad

· grad·Iwashita lab.

 Frequency domain analysis◦→i◦ The coefficient matrix results in complex symmetric.
◦ IC-COCG method is often used to solve the linear system.

 Time domain analysis◦ When the backward time difference scheme is used,
,

→

.

→

◦ The resulting linear system:
(11)

1
Δ

The right-hand side vector b is determined by
the external current density andandin theprevious time step.

Iwashita lab.

 The coefficient matrix of the linear system (11) issymmetric and semi-positive definite.
 The standard solver is ICCG method.

◦ The application of AMG to the edge element analysis is
not easy.
◦ Other preconditioning techniques and iterative solvers
are used in some cases, for example, implementation on
GPU.

 Application of the gauge condition◦ Applying =0 to the conductive region, we get Aformulation. The resulting linear system:′.
is still singular. Applying the tree co-tree gauge◦
condition to the non-conductive region, we get a nonsingular coefficient matrix.

Iwashita lab.

Co-tree edge

 The rotation of the field in the element can be

determined by the valuable on the co-treeedge even if the valuable on the tree edges isfixed at an arbitrary value.

Iwashita lab.

When the =0 and the tree co-tree gauge conditions areused, the number of unknowns are reduced.
 However, it is known that the convergence of the ICCG
method is significantly deteriorated.
◦ Usually, the computational time is increased when these gauge
conditions are explicitly enforced.



Why?

◦ References Prof. Honma, Prof. Igarashi, and Prof. Kawaguchi’s text book in Japanese:
ISBN4-267-71641-9
 K. Fujiwara, T. Nakata, and H. Ohashi, “Improvement of convergencecharacteristic of ICCG method for the A – method using edge elements,”
IEEE Trans. Magn., vol. 32, pp. 804–807, May 1996.
 H. Igarashi and T. Honma, “On convergence of ICCG applied to finite
element equation for quasi-static fields,” IEEE Trans. Magn., vol. 38, no. 2,
pp. 565–568, Mar. 2002.

Iwashita lab.

 The large null space of discrete curl operator◦ The coefficient matrix arising in the A-formulation: The matrix C has large null space, which is sifted by
the positive-definite matrix. When the
conductivity is small, the shift generates small positive
eigenvalues, which causes worse convergence.
 The error belonging to the null space of the curl
operator may converge slowly (discussion in part II).

◦ In the linear system based on the A- formulation,
the removal of a redundant unknown causes asmaller eigenvalue in the reduced dimensionalcoefficient matrix.

Iwashita lab.

 Difficulty in AMG◦ The null space matching between fine and coarse
grids Difficulty in high-frequency analyses◦ Indefinite matrix◦ Influence upon the condition number

Iwashita lab.

Part II
Explicit / Implicit Error Correction FrameworkCollaborators

Takeshi Mifune and Masaaki Shimasaki(Kyoto Univ.)

References

Takeshi Iwashita, Takeshi Mifune and Masaaki Shimasaki, “Similarities
Between Implicit Correction Multigrid Method and A-phi Formulation inElectromagnetic Field Analysis”, IEEE Transaction on Magnetics, Vol. 44,
No. 6, (2008), pp. 946-949.
Takeshi Mifune, Soichi Moriguchi, Takeshi Iwashita and Masaaki Shimasaki,“Convergence Acceleration of Iterative Solvers for the Finite Element
Analysis Using the Implicit and Explicit Error Correction Methods”, IEEE
Transaction on Magnetics, Vol. 45, No. 3, (2009), pp. 1104-1107.
Takeshi Iwashita, Takeshi Mifune, Soichi Moriguchi and Masaaki Shimasaki,“Physical Meaning of the Advantage of A-phi Method in Convergence”, IEEE
Transaction on Magnetics, Vol. 45, No. 3, (2009), pp. 1424-1427.

Outline
We focus on some error correction methods initerative solver.
 We newly define Explicit error correction method fora class for some conventional techniques.
 We propose Implicit error correction method.
Then, we establish Explicit / Implicit error correction
framework

Introduction of two techniques described in the
proposed framework
 Numerical results confirm the effectiveness of implicit
error correction method.
Explicit error correction method(subspace correction method)

Definition Ax=b(1) : A linear system to be solved

A: n by n matrix, x, b: n-dimensional vector

If an error correction method is given by thefollowing two steps, we call it one of “explicit
error correction methods”.
 Step 1: Determination of quantity of error
correction
~
Cy
=
d
(x
)
 Solve or use(2)
C : m by m matrix, d : a function of thex.
approximate solution ~Explicit error correction method~Step 2: Error correction (update of x )


new
~
~
~
x
← x + By(3)

~y : the approximation or exact y

B : n by m matrixThe dimension m is (sufficiently) smaller than n inpractical use.

Characteristics of explicit error
correction method
Special treatment for errors thatconverge slow.
 Use of a different matrix from a
original coefficient matrixFamous examples
 Hiptmair’s hybrid smoother Coarse grid correction in a multigrid
methodImplicit error correction method

We propose a new technique “implicit error
correction method”.
The idea comes from solving the original linear
system (1) and the linear system for error correction(2) simultaneously.
Derivation of the proposed method

Paying special attention to (3), the solution
vector x is replaced by xˆ + By
ˆ in (1).
Then, we get

Axˆ + AByˆ = b (4).

Implicit error correction method

Next, we consider (2). In most explicit error
correction methods, d is given by a restricted
residual equation of (1):d = D(b − A~
x) .

~
DA
x + Cy = Db .
 Thus, we obtain

By rewriting vectors ~x and y by x̂ and ŷ ,respectively, we finally obtain A AB  xˆ   b 
 DA C  yˆ  =  Db  (6).
   Implicit error correction methodThe above linear system (6) is the formulation
based on the implicit error correction method.

The linear system (6) is solved by a preconditioned
iterative solver in the proposed method.
 It is expected that the iterative solution process
produces the effect of the error correction implicitly.
In other words, the linear system (6) has an
improved condition of the coefficient matrix overthat of the original linear system (1).
Electromagnetic field simulation
The advantage in convergence in the A-phi
method has not been explained in physicalmeaning.
 We show that the A-phi method can be
regarded as a implicit error correction method
corresponding to Hiptmair’s hybrid smoother.
 Numerical tests show the A-phi method has
similar error correction effect as the hybridsmoother.
Two formulations inelectromagnetic field simulation

A-method∂Am∇ × (ν∇ × Am ) + σ
= J0
∂t

(7)

∂ ( Am′ + gradφ )∇ × (ν∇ × Am′ ) + σ
= J0
∂t

(8)

A-phi method∂ ( Am′ + gradφ )
div(σ
)=0
∂t

(9)

Am, Am’ : the magnetic vector potential,
ν: the magnetic reluctivity, σ: the electrical conductivity,
t: time, J0: the exciting current, φ, the electric scalar potential∇ ⋅ J 0 = 0 is satisfied.

Linear system based on Aformulation
Application of the finite edge element discretization and abackward time difference method to (7) results in

Mσ 
 T
[K A ] = Cu Mν Cu + {xA } = {bA }
Δt 


(10)

where ∆t: the length of the time step, xA: the unknown vector for magnetic vector
potential, bA is determined by the previous value of the magnetic vector potential
and the exiting current.

Cu: the discrete curl operator
M v = ∫νwif ⋅ w jf dvΩ

M σ = ∫ σwie ⋅ w ej dv (11)Ω

wf: face-element basis function, we: edge-element basis function

We consider the application of Hiptmair’s
hybrid smoother to this linear system.

Relationship between A-phi method andHiptmair’s hybrid smootherWhen the hybrid smoother is applied to the linear
system (A-method) (10), its error correction
procedure can be written in a form of explicit error
correction method as follows:

C = (G M σ G ) / ∆tT~d ( x ) = G (bA − K A ~x A ) in (2) and (3)T

B=G
G : discrete gradient operator, CuG=0 is satisfied.

Corresponding implicit error correction
method…

 KA
 1 T ∆t G M σ

M σ G  xˆ A   bA   =  T T1ˆ
 G bA ∆t G M σ G  y A 
1∆tThis linear system coincides withthe linear system arising on
the A-phi formulation.
The A-phi method can be regarded as a implicit error correction
method corresponding to Hiptmair’s hybrid smoother.
Thus, it is expected that the A-phi method has a similar effect to thatof the hybrid smoother for the error of the kernel of the discrete curl
operator.
Whereas the hybrid smoother corrects these errors explicitly, the A-phi
method introduces unknown variables for the electrical scalar potential.

Numerical tests (Effectiveness ofA-phi method)
Test model: TEAM workshop problem 10
 Whitney elements, # DOFs 5968
-3 ∆t =10 (sec)Well-known property
This phenomena has been
explained by the condition
number of coefficientmatrices.

Comparison behavior (relative residual norm)

Numerical tests (Effectiveness ofA-phi method)
We focus on the error of the kernel of thediscrete curl operator K(curl) in conductive
region.
 This error is explicitly corrected in Hiptmair’s
hybrid smoother.
• The error of K(curl) in
conductive region is efficiently
removed in the A-phi method.
• Numerical tests confirm that
the A-phi method has a similareffect as that of the hybridsmoother for the error ofK(curl) in conductive region.

Numerical tests (Effectiveness ofA-phi method)
The error of K(curl) in conductive region isefficiently removed in the A-phi method.
 The slow convergence of the A-method iscaused by the error of K(curl) in conductive
region.
 Numerical tests confirm that the A-phi
method has a similar effect as that of thehybrid smoother for the error of K(curl) in
conductive region.
 We give the first physical explanationof the advantage of A-phi method.
Summary of part II
We propose an implicit error correction method that
corresponds to the explicit error correction methods such
as Hiptmair's hybrid smoother.
 It is shown that the A-phi method can be regarded as animplicit error correction version of the hybrid smoother.
 Numerical tests show that the A-phi method has a similareffect on the correction of the error of the kernel of thediscrete curl operator, which results in an advantage inconvergence.
Part III : TP-EEC Method
 Main contributors◦ Yasuhito Takahashi (Doshisha Univ.), Tadashi
Tokumasu (Toshiba Co. Ltd.) Reference◦ Y. Takahashi, T. Tokumasu, A. Kameari, H.
Kaimori, M. Fujita, T. Iwashita, and S. Wakao,“Convergence Acceleration of Time-Periodic
Electromagnetic Field Analysis by theSingularity Decomposition-Explicit Error
Correction Method,” IEEE Transactions on
Magnetics, Vol. 46, No. 8, 2010.

Target problem
 Motors and transformers (practicalsimulations) Low frequency field problem (50-60Hz)◦ Eddy current problem
 Non-linear problem◦ Non-linear magnetic property
 Calculation for time-periodic steady statesolution

Basic equations
 A-phi formulation∂ ( A + gradφ )∇ × (ν∇ × A) + σ= J0∂t∂ ( A + gradφ ))=0div (σ∂tin conductive region

A : the magnetic vector potential, ν: the magnetic reluctivity, σ: the electricalconductivity, t: time, J0: the exiting current, φ, the electric scalar potential
PDE problem -> FEM
ν has non-linear property: ν = ν (rot A) -> Newton Raphson method
Transient analysis should be performed.
Treatment of the time derivative -> (1) Step by step method(2) TPFEM

Conventional
methodTransient analysis∂
S ( x) + C
x= f∂t

A system of nonlinear equations derived
from the finite element discretization

Adopting the θ method for the time integration scheme~~~
T ( x i ) − C ( x i −1 ) = f i(The subscripts indicate time step.)

Time periodic conditionC~T ( xi ) = θ S ( xi ) +
xi
∆t
C
~
C ( x i ) = −(1 − θ ) S ( x i ) +xi∆t~f i = θ f i + (1 − θ ) f i −1+

 Ordinary time-periodic condition: x(t+T) = x(t) Half time-periodic condition: x(t+T/2) = −x(t)T is the period.

Step by step method
In the step-by-step method, a series of the following non-linear
equations are repeatedly solved until the steady state solution is
obtained.

~ ~~Gi = fi − T (x i ) + C (x i −1 ) = 0(i = 1, 2,, nt )

 0 Initial guess (only used in the first time step)x0 =  xnt Ordinary time-periodic condition case− x Half time-periodic condition case nt
Each non-linear equation has n unknowns, where n is the number of unknowns of
the analyzed model (discrete vector and scalar potentials).

Overview of the simulation code
based on the step-by-step method The simulation code consists of three nested loop (fortime integration, Newton method, linear solver).
Loop for time-dependent calculations
Loop for Newton Raphson method
Loop for linear iterative solver
ICCG, ILU-GMRES, ILU-BiCGSTAB etc.

Computational cost: C x n x Nite x Nr x NT
C : a constant value determined by the linear iterative solver
Nite : the average number of iterations of iterative solver
Nr : the average number of non linear iterations
NT : the number of total time steps

A problem of the step by step
methodThe transient analysis is a kind of initial value problems of ODE.
If the time constant of the problem is very large, many time steps are needed
to obtain the steady state solution.

Coil current [A]In practical simulations, the analyzed model involves an external circuit,which often leads a large time constant.

25
20
15
10
5
0

One period

I

0

0.05
Time [s]

0.1

0.15

(

)

0.2

FEM V = 1000 1 − cos ωtregionTime-periodic finite-element
method (TPFEM)
 All linear systems are combined with thetime periodic condition and are solved
simultaneously.
~~ T ( x1 )  C ( xn )   f1 
  
 ~
~
 T ( x2 ) − C ( x1 )   f2  T~( x3 ) − C~( x2 )  =  f3 

  

  

T~( xn ) − C~( xn −1 )  fn tt
  t


Solve this large system of
non-linear equations bymeans of Newton Raphson
method.

However, the convergence of the Newton method is not good, and thus
TPFEM is usually inferior to the step by step method in the computational time.
But, TPFEM has an advantage in parallel computation.

Main idea Acceleration of convergence in the transient

analysis by using EEC (Explicit Error
Correction) (Subspace correction) method
The problem: linear system of equations Kx=b
The approximation vector 𝑥𝑥� is updated as follows:r = b − Kx~

B KBy = B r is (approximately) solved.
~~~
xnew = x + By
T

T

EEC method for non-linear
equationK ( x ) = b : the nonlinear equation to solve
K ( x~ + e ) = b

Linearization
Equation for the error

~)
∂
K
(
x
K ( x~ ) +e=b∂x∂K ( x~ )e = b − K ( x~ )
∂x

~)
(
x
∂
K
T
T
for error
B
BeH = B ( b − K ( x~ )) Equation
correction (lower-rank)∂x
newEquation for error~~x← x + BeH

correction (lower-lank)

EEC ProcedureWhat kind of errors are slow inconvergence?
 Tokumasu and Takahashi’s idea

◦ The DC component of the error can have a
low convergence rate.
◦ How is it extracted ?
 Use TPFEM formulation (for K(x)=b)
 The DC component of the error has the identicalvalue of potentials in all time steps.
The DC error component is written as
𝑒𝑒1 , 𝑒𝑒2 , ⋯ , 𝑒𝑒𝑛𝑛 , 𝑒𝑒1 , 𝑒𝑒2 , ⋯ , 𝑒𝑒𝑛𝑛 , ⋯ ⋯ , 𝑒𝑒1 , 𝑒𝑒2 , ⋯ , 𝑒𝑒𝑛𝑛 T

in TPFEM formulation.

It is rewritten as

𝑰𝑰 𝑰𝑰 ⋯ 𝑰𝑰 𝑇𝑇 𝑒𝑒1 , 𝑒𝑒2 , ⋯ , 𝑒𝑒𝑛𝑛 T

BeH

I : n x n unit matrix

TP-EEC method(1) Perform the step by step method for nt

time steps(2) Perform the following error correction forthe current approximation vector xnt
~)
∂
K
(
x
T
Solve B T
BeH = B ( b − K ( x~ ))
∂x
T
T
T
T
~
x = ( x0 , x1 ,, xnt )B = ( I , I ,, I )T (n x nt) by n matrix
K : the coefficient matrix in TPFEM
Update xnt to xnt + BeH(3) If the convergence condition is notsatisfied, return to (1) with x0 = ± xnt .

Simplified TP-EECThe technique is available only for a half time periodic problem.

()

 nt ~ ~~ ~~ ∑ T i − C i + 2 C nt eH = −C ( x0 ) − C ( xnt ) The equation for error
correction
 i =1• Assume that e H = αq = α ( − x0 − xnt )• Adopt the backward Euler method (θ=1)• Multiply both sides by q T

α=1(∆t q (∑ S )q q C q + 2)
T

n

i =1

T

i

When the influence of the matrix C is larger than Si, which corresponds to thecase that the time constant is large, α = 1/2

1e H = (− x0 − xnt )2

The simplified TP-EEC method, which is
unnecessary to solve the linear system forthe error.

Numerical tests (Iron core model)
SS400(σ=7.505×106)

200 µFI

Coil:3000 Turn

Half periodicity: V=1000 sinωt
FEM region

Methods

f =50 Hzor

Ordinary periodicity: V=1000(1-cos ωt)

Elapsed time (s)Half periodic problem

Periodic problem

Normal step-by-step

3360.4 (207)6848.7 (432)

TP-EEC method648.7 (39)1999.3 (121)

Simplified method

1101.8 (68)

Numbers in parenthesis indicate the number of time steps required for
obtaining steady state solutions in one period

3.0

Half timeperiodic problem

Coil current [A]2.0

Coil current [A]

TP-EEC
NormalSimplified
Solution
TPFEM

1.0
0.0
-1.0
-2.0
7 0
6
5
4
3
2
1
0
-1
-2
-3
0

Converge to steady state

0.02

0.04
0.06
Ordinary timeperiodic problem

0.08
TP-EEC
NormalSolution
TPFEM

0.1

0.02

0.04

0.08

0.1

Time [sec]0.06

Converge to steady stateRelative error of eddy-current loss [%]Half time-periodic problem

1.0E+2
1.0E+1
1.0E+01.0E-1
1.0E-2
1.0E-3
1.0E-4
1.0E-5
1.0E-6

TP-EEC
Simplified
Normal

1.0E+3 0

20

60 problem 80
Ordinary40time-periodic

100

1.0E+2
1.0E+1
1.0E+01.0E-1TP-EEC

Normal

1.0E-2

00

0.02
20

0.04
40
Time [s]

0.06
60

0.08
80

0.1
100

Numerical test (IPM motor)
permanent
magnetcoil Mesh of IPM motor (1/6 model)（#Elements：17,556，#Unknowns：49,798）
The stator is satisfied with half timeperiodicity

+
shaft

The rotor has ordinary time-periodicity.
simultaneously considered in the TP-EEC method
The simplified method is only applied to the stator

Methods

Elapsed time (s)Normal step-by-step for 30 periods

Not converged

TP-EEC method

4167 (244)

Simplified TP-EEC method4079 (242)Eddy-current loss [W]Torque [Nm]

TPFEM
Solution

Normal

12000
10000
8000
6000
4000
2000
0
-2000

Simplified

TP-EEC

Converge to steady state
0

180

360

540

720

900

1080

0

180

360

540

720

900

1080

6000
5000
4000
3000
2000
1000
0

Electric angle [deg]Summary for part III
 The TP-EEC method was introduced toaccelerate the convergence to the steadystate solution in a transient analysis.
◦ The TP-EEC method is based on the error
correction method for the non-linear equation.
◦ It focuses on the error correction of the DC
component of the error.
◦ Numerical tests confirms that the TP-EEC
method improve the convergence.

Related work
 Parallel TP-EEC method was proposed.

Each process/thread independently performthe step-by-step method. (Parallel in time)(2) The continuity of the solution between
processes/threads and the periodic boundary
condition should be satisfied. The EEC
framework for non-linear equation is applied
for reducing the gap.
(3) The technique can be applied to any time
evolution analysis (non time periodic problems).
(4) Parallel TP-EEC can be regarded as one of
PARAREAL method.
(1)

Y. Takahashi, T. Tokumasu, K. Fujiwara, T. Iwashita, and H. Nakashima,“Parallel TP-EEC Method Based on Phase Conversion for
Time-Periodic Nonlinear Magnetic Field Problems,” IEEE Transactions
on Magnetics, Vol. 51, No. 3, 2015.

主たる研究者： 美舩健（京大）
 共同研究者： 髙橋康人（同志社大）
 References◦ Takeshi Mifune, Yasuhito Takahashi and Takeshi Iwashita, “Folded
Preconditioner: a New Class of Preconditioners for Krylov
Subspace Methods to Solve Redundancy-Reduced Linear Systemsof Equations”, IEEE Transactions on Magnetics, Vol. 45, No. 5,(2009), pp. 2068-2075.
◦ Takeshi Mifune, Yasuhito Takahashi, and Takeshi Iwashita, “NewPreconditioning Technique to Avoid Convergence Deterioration
due to the Zero-Tree Gauge Condition in Magnetostatic Analysis”,
IEEE Trans. Magn., Vol. 46, No. 7, (2010), pp. 2579-2584.
◦ Yasuhito Takahashi, Takeshi Mifune, Takeshi Iwashita, Koji
Fujiwara and Yoshiyuki Ishihara, “Folded IC Preconditioning in
Quasi-Static Field Analysis Taking Account of Both Tree-Cotree
and φ=0 Gauge Conditions”, IEEE Trans. Magn., Vol. 47, No. 5,(2011), pp. 1342-1345.

Iwashita lab.

Ax = b
N: 方程式の未知数の数
A : 係数行列 (N×N)
x : N次の解ベクトル（未知）b : N次の右辺ベクトル（既知）※本資料では A が正方(N×N)の場合のみに限定

Iwashita lab.

Ax = b(A ∈ C

N ×N,

x ∈C N , b∈C N)

 第 i ステップの近似解： x(i)

※簡単のため、本資料では x(0) = 0 を仮定
 決められた手順で x(i) を更新する( i = 1, 2, …)
 残差 r(i) = b – A x(i) が十分小さくなったら収束とみなす

Iwashita lab.

 クリロフ部分空間：

K

(i )

( A; b ) = span{b, Ab, , A b}i −1

 第 i ステップの近似解を K(i)(A, b) から探索

x ∈K
(i )

（※本資料では(i )( A; b )x(0) = 0 を仮定している）{}
x ( 2 ) ∈ span{b, Ab}

つまり、 x (1) ∈ span b{}

x (3) ∈ span b, Ab, A2 b…
Iwashita lab.

 前処理：収束に必要な反復数を減らし、計算を高速化するための工夫
元の方程式（右）前処理された方程式

Ax = b

M : 前処理行列(N×N)
AM −1Mx = b
A′ = A M －1
x′ = M x

A′x ′ = b

Iwashita lab.

 各ステップの近似解を K(i)(A′, b) から探索{}

x ′(i ) ∈ K (i ) ( A′; b ) = span b, A′b,  , A′i −1b

A′ = AM-1, 元の方程式の解は x(i) = M-1 x′(i) で求まるので,

{

−1

−1

−1

x ∈ span M b, M AM b,  , M(i )(∈ K (i ) M −1 A; M −1b)−1(AM ) b}−1 i −1

M が適切に選ばれていれば、反復収束性が改善される（通常、M ≃ A を目指す）

Iwashita lab.

 1 1  x1   1  解は x1 + x2 = 1 を
det (A) = 0
例：   =  満たすもの全て
x
2
2
2




2


(rank(A) < N)
 解は無数に存在し一意に決まらない（解が存在しない場合はここでは扱わない）
 電磁界解析では，しばしば現れる（無数の解のうち、ある一つの解を得ればよい）◦ curl curl の形（静磁界）の式の離散化◦ 辺・節点要素を併用する A-φ法

Iwashita lab.

Ax = b, rank(A) = L < N
未知数, 式の順序を適切に選べば, 一般に
L
L
N-L

 Ar

＊r
 CA

N-L

A＊
 br 
r B  x r   = CA＊r B  xex   C
＊br 

L
N-L

Ar ∈C L×L は正則, B ∈CL×(N – L), C ∈C(N – L)×L

CE Ar BE x = CE br
BE = (EE
B ), CE =   E ∈CL×L は 単位行列
C 
Iwashita lab.

Ar x rCAr x r+ Ar Bxex = br
+ CAr Bxex = Cbr
Ar BE x = br

下の式は上の式に C を乗じただけなので不要( BE x = xr + Bxex )y = BE x ∈CL とおく

Ar y = br L 次に縮小された方程式（ y は一意）
y が求まれば, 元の式の解の一つが ,xr ← y, xex ← 0 で求まる
Iwashita lab.

ベクトルポテンシャルとスカラポテンシャルを未知数として定式化（辺要素・節点要素を併用した有限要素法）

波動方程式
離散化

 Ar

 CAr

ベクトルポテンシャルに対応する未知数

Ar B  x r   br   = 
CAr B  xex   Cbr 

スカラポテンシャルに対応する未知数

Ar y = br
未知数の数を削減（A法）

Iwashita lab.

 電磁波解析の例題

A-φ 法

A法

未知数の数

483,887

363,107行列の非零要素数（上三角部）

13,930,984

6,006,199

前処理－クリロフ部分空間法

IC－COCG

IC－COCG

メモリ消費 (MB)

663

356

収束までの反復回数302

求解時間 (s)

98.8

1772
265.5

Iwashita lab.

不定な方程式 Ax = b, rank(A) = L < N

A = CE Ar BE , b = CE br ………(※１)
ただし, BE = (EB )∈ C

L× N

E, CE =   ∈ C N × LC 

Ar ∈C L×L は正則, E ∈C L×L は単位行列縮小した方程式 Ar y = br ( y = BE x )反復の収束が悪く, 計算時間面でかえって不利に（メモリは節約できる ）Iwashita lab.

 FP (Folded Preconditioning)

不定な方程式に対する
前処理
IC
SOR
MG
…

縮小した方程式に対する
前処理
IC
SOR
MG
…
FIC
FSOR
FMG
…

FP

Iwashita lab.

A = CE Ar BE についての前処理 Mに対して,
以下で定義される L × L 行列 Mf を使う(Arについての)前処理を FP と呼ぶ

M

−1
f

= BE M CE ……………(※２)−1

Iwashita lab.

y = BE x

不定な方程式( N 次)
Ax = b

ベースとする前処理付 (i)(i)
きクリロフ部分空間法 y = BE x(x (i ) ∈ K (i ) M −1 A; M −1b)

縮小した方程式( L 次)
Ar y = bry(i) は, 縮小した方程式の近似解とみなせる

x(i) を求めるアルゴリズムと y(i) を求めるアルゴリズムは,元の方程式の残差の観点から, 同等 （次スライド）
2. y(i) を求める反復は, 縮小した方程式に対する前処理付
きクリロフ部分空間法の形で得られる。ただしベースとした前処理に対応する FP を用いる （次々スライド）
1.

Iwashita lab.

r (i ) = b − Ax (i )= b − C E Ar BE x (i ) BE x (i ) = b − C E Ar BE 0 y (i ) = b − A0( BE = (EB ))………(※３)

Iwashita lab.

(x (i ) ∈ K (i ) M −1 A; M −1b()(※１)よりx (i ) ∈ K (i ) M −1C E Ar BE ; M −1C E br()

y ( i ) = BE x (i ) ∈ K (i ) BE M −1C E Ar ; BE M −1C E br(※２)より(y (i ) ∈ K (i ) M f−1 Ar ; M f−1br))

縮小した方程式に Mf を前処理として用いたときの探索空間に等しい
Iwashita lab.

Ax = b に対する（任意の）前処理付きク
リロフ部分空間法をベースとして, （※３）の意味で同等のアルゴリズムを, Ar y = br に対する前処理付きクリロフ部分空間法の形で設計できる。ただし ベースとした前処理に対応する FP を用いる。

Iwashita lab.

Ax = b
前処理付き
クリロフ部分空間法
IC-CG
SOR-BiCG
MG-BiCGSTAB

Ar y = br
前処理付き
クリロフ部分空間法
同等（同じ収束性）

FIC-XXXFSOR-YYY
FMG-ZZZ

MG-WWW

Iwashita lab.

(1)ベースとした前処理付きクリロフ部分空間法と同等になる探索アルゴリズム（前スライドのXXX, YYY,
ZZZ）を設計する(2)既知の効率的に設計されたアルゴリズム（CG,
BiCGSTABなど）を用いる （推奨）
 特別な場合◦ ベースとする探索アルゴリズムがCG, BiCGなどの時、(１)と(２)が同じ結果になることもある

Iwashita lab.

x(i) は, 双直交条件を満たすよう決定されるはシャドウ(r , (M A ) b ) = 0, 0 ≤ j < i （b*ベクトル）(C r , (M A ) b ) = 0, 0 ≤ j < i
(i )

−H

H

−H

(i )
E r

j

*

H

j

*

(※１) (※２)より(r , (M A ) C b ) = 0, 0 ≤ j < i
(i )
r

−H
f

H
r

j

H
E

*

Iwashita lab.

Ax = b
前処理付き
クリロフ部分空間法
IC-CG
SOR-BiCG
MG-BiCGSTAB

Ar y = br
前処理付き
クリロフ部分空間法
同等（同じ収束性）

FIC-CGFSORBiCG
FMG-ZZZFMGBiCGSTAB（“同じ収束性”は保証されな
いが、性能は期待できる）
Iwashita lab.

探索アルゴリズム(CGなど)の部分については, 未知数を減らしたことで, 計算コストの削減が期待できる
前処理の実装に必要なのは, 行列ベクトル積： z = Mf-1
w のみ

M f−1 = BE M −1CE・ M -1 についての積－ ベースとした前処理のコストと同じ
・ BE, CE についての積－ 効率よく計算できれば成功
Iwashita lab.

A-φ 法の方程式では,
 E BE = (E G ), CE =  T 
G 

なので, G, GT に関する積が速く計算できれば良い
ここで G は, 離散勾配行列（次スライド）

Iwashita lab.

有限要素メッシュ
G : （辺数）× （節点数）行列
節点 n1
辺e
節点 n2

G の 第 e 行の非零成分は(e, n1), (e, n2) 成分の２つのみ
G, GT とベクトルの積は十分高速に行える

Iwashita lab.

 電磁波解析
cavity surface

ε r = 6.0-j1.0
dielectric
magnetic wall

a

要素数
節点数
ポート数
周波数

121,500
129,0762
約6.2GHz(k0a = 2.6)

PC(Xeon X5472, 8GB RAM)

Iwashita lab.

Formulations

A-φ

A

未知数の数

483,887

363,107係数行列の非零要素数(上三角部分)13,930,984

6,006,199

A-φ
解法

A

IC－COCG IC－COCG FIC－COCG

FIC(改良版)－
COCG

メモリ消費
(MB)

663

356

542

524

反復回数

302

1772

302

302

求解時間 (s)

98.8

265.5

89.4

81.8

Iwashita lab.

静磁界解析

直方体辺要素70×70×70 分割（辺数は 1,058,610）

z [m]

0.2

Coil

J

yz, zx 平面は固定境界
μr = 1000

0.12

Tree-Zero ゲージについて
は, x 軸に平行な辺全てをTree として選択

0.11
0.1

0.075 0.085

0.101

0.2

x [m]

PC(Core2 X9770, 8GB RAM)

0.1

Iwashita lab.

Formulations

ゲージ無し

Tree-Zero ゲージ

未知数の数

1,058,610

690,900係数行列の非零要素数(上三角部分)17,297,701

7,443,450

ゲージ無し

Tree-Zero ゲージ

解法

IC－CG

IC－CG

FIC－CG

メモリ消費
(MB)

496

240

378

反復回数

673

10000 以上

679

求解時間 (s)

134.1

－

111.6
Iwashita lab.

 非正則な係数行列を持つ連立一次方程式に対する

Folded preconditioning（FP）を提案
 同前処理をクリロフ部分空間法に適用した場合，収束性の劣化を伴わずに冗長な未知変数を削除できる◦ 一部の反復法では等価な収束性が得られる

 FPを使用する際のコストは反復前に評価できる

Iwashita lab.