Rheolef  7.1
an efficient C++ finite element environment
mixed_solver.h
Go to the documentation of this file.
1 # ifndef _RHEO_MIXED_SOLVER_H
2 # define _RHEO_MIXED_SOLVER_H
3 // This file is part of Rheolef.
4 //
5 // Copyright (C) 2000-2009 Pierre Saramito <Pierre.Saramito@imag.fr>
6 //
7 // Rheolef is free software; you can redistribute it and/or modify
8 // it under the terms of the GNU General Public License as published by
9 // the Free Software Foundation; either version 2 of the License, or
10 // (at your option) any later version.
11 //
12 // Rheolef is distributed in the hope that it will be useful,
13 // but WITHOUT ANY WARRANTY; without even the implied warranty of
14 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
15 // GNU General Public License for more details.
16 //
17 // You should have received a copy of the GNU General Public License
18 // along with Rheolef; if not, write to the Free Software
19 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
20 //
21 // =========================================================================
22 
23 #include "rheolef/uzawa.h"
24 #include "rheolef/cg.h"
25 #include "rheolef/minres.h"
26 namespace rheolef {
27 
28 template <class Matrix, class Vector, class Solver>
29 struct abtbc_schur_complement { // S = B*inv(A)*B^T + C
30  Solver iA;
31  Matrix B, C;
32  abtbc_schur_complement (const Solver& iA1,const Matrix& B1, const Matrix& C1) : iA(iA1), B(B1), C(C1) {}
33  Vector operator* (const Vector& p) const {
34  Vector v = iA.solve(B.trans_mult(p));
35  return B*v + C*p;
36  }
37 };
38 template <class Matrix, class Vector, class Solver>
39 struct abtb_schur_complement { // S = B*inv(A)*B^T
40  Solver iA;
41  Matrix B;
42  abtb_schur_complement (const Solver& iA1,const Matrix& B1) : iA(iA1), B(B1) {}
43  Vector operator* (const Vector& p) const {
44  Vector v = iA.solve(B.trans_mult(p));
45  return B*v;
46  }
47 };
48 template <class Matrix, class Vector, class Solver, class Preconditioner>
49 int uzawa_abtbc (const Matrix& A, const Matrix& B, const Matrix& C, Vector& u, Vector& p,
50  const Vector& Mf, const Vector& Mg, const Preconditioner& S1,
51  const Solver& inner_solver_A, const Float& rho,
52  const solver_option& sopt = solver_option())
53 {
54  sopt.label = (sopt.label != "" ? sopt.label : "uzawa_abtbc");
55  abtbc_schur_complement<Matrix,Vector,Solver> S (inner_solver_A, B, C);
56  Vector v = inner_solver_A.solve (Mf);
57  Vector Mh = B*v - Mg;
58  int status = uzawa (S, p, Mh, S1, sopt);
59  u = inner_solver_A.solve (Mf - B.trans_mult(p));
60  return status;
61 }
62 template <class Matrix, class Vector, class Solver, class Preconditioner, class Real>
63 int uzawa_abtb (const Matrix& A, const Matrix& B, Vector& u, Vector& p,
64  const Vector& Mf, const Vector& Mg, const Preconditioner& S1,
65  const Solver& inner_solver_A, const Real& rho,
66  const solver_option& sopt = solver_option())
67 {
68  sopt.label = (sopt.label != "" ? sopt.label : "uzawa_abtb");
69  abtb_schur_complement<Matrix,Vector,Solver> S (inner_solver_A, B);
70  Vector v = inner_solver_A.solve (Mf);
71  Vector Mh = B*v - Mg;
72  int status = uzawa (S, p, Mh, S1, rho, sopt);
73  u = inner_solver_A.solve (Mf - B.trans_mult(p));
74  return status;
75 }
76 /*Class:mixed_solver
77 NAME: @code{cg_abtb}, @code{cg_abtbc}, @code{minres_abtb}, @code{minres_abtbc} -- solvers for mixed linear problems
78 @findex cg
79 @findex minres
80 @findex cg\_abtb
81 @findex cg\_abtbc
82 @findex minres\_abtb
83 @findex minres\_abtbc
84 @cindex mixed linear problem
85 @cindex conjugate gradien algorithm
86 @cindex finite element method
87 @cindex stabilized mixed finite element method
88 @cindex Stokes problem
89 @cindex incompresible elasticity
90 SYNOPSIS:
91  @example
92  template <class Matrix, class Vector, class Solver, class Preconditioner>
93  int cg_abtb (const Matrix& A, const Matrix& B, Vector& u, Vector& p,
94  const Vector& Mf, const Vector& Mg, const Preconditioner& S1,
95  const Solver& inner_solver_A, const solver_option& sopt = solver_option());
96 
97  template <class Matrix, class Vector, class Solver, class Preconditioner>
98  int cg_abtbc (const Matrix& A, const Matrix& B, const Matrix& C, Vector& u, Vector& p,
99  const Vector& Mf, const Vector& Mg, const Preconditioner& S1,
100  const Solver& inner_solver_A, const solver_option& sopt = solver_option());
101  @end example
102  The synopsis is the same with the minres algorithm.
103 
104 EXAMPLES:
105  See the user's manual for practical examples for the nearly incompressible
106  elasticity, the Stokes and the Navier-Stokes problems.
107 
108 DESCRIPTION:
109  @noindent
110  Preconditioned conjugate gradient algorithm on the pressure p applied to
111  the stabilized stokes problem:
112  @example
113  [ A B^T ] [ u ] [ Mf ]
114  [ ] [ ] = [ ]
115  [ B -C ] [ p ] [ Mg ]
116  @end example
117  where A is symmetric positive definite and C is symmetric positive
118  and semi-definite.
119  Such mixed linear problems appears for instance with the discretization
120  of Stokes problems with stabilized P1-P1 element, or with nearly
121  incompressible elasticity.
122  Formally u = inv(A)*(Mf - B^T*p) and the reduced system writes for
123  all non-singular matrix S1:
124  @example
125  inv(S1)*(B*inv(A)*B^T)*p = inv(S1)*(B*inv(A)*Mf - Mg)
126  @end example
127  Uzawa or conjugate gradient algorithms are considered on the
128  reduced problem.
129  Here, S1 is some preconditioner for the Schur complement S=B*inv(A)*B^T.
130  Both direct or iterative solvers for S1*q = t are supported.
131  Application of inv(A) is performed via a call to a solver
132  for systems such as A*v = b.
133  This last system may be solved either by direct or iterative algorithms,
134  thus, a general matrix solver class is submitted to the algorithm.
135  For most applications, such as the Stokes problem,
136  the mass matrix for the p variable is a good S1 preconditioner
137  for the Schur complement.
138  The stopping criteria is expressed using the S1 matrix, i.e. in L2 norm
139  when this choice is considered.
140  It is scaled by the L2 norm of the right-hand side of the reduced system,
141  also in S1 norm.
142 AUTHOR:
143  | Pierre.Saramito@imag.fr
144  LJK-IMAG, 38041 Grenoble cedex 9, France
145 DATE: 15 april 2009
146 METHODS: @mixed_solver
147 End:
148 */
149 
150 template <class Matrix, class Vector, class VectorExpr1, class VectorExpr2,
151  class Solver, class Preconditioner>
152 int cg_abtbc (const Matrix& A, const Matrix& B, const Matrix& C, Vector& u, Vector& p,
153  const VectorExpr1& Mf, const VectorExpr2& Mg, const Preconditioner& S1,
154  const Solver& inner_solver_A, const solver_option& sopt = solver_option())
155 {
156  sopt.label = (sopt.label != "" ? sopt.label : "cg_abtbc");
157  abtbc_schur_complement<Matrix,Vector,Solver> S (inner_solver_A, B, C);
158  Vector v = inner_solver_A.solve (Mf);
159  Vector Mh = B*v - Mg;
160  int status = cg (S, p, Mh, S1, sopt);
161  u = inner_solver_A.solve (Mf - B.trans_mult(p));
162  return status;
163 }
164 template <class Matrix, class Vector, class VectorExpr1, class VectorExpr2,
165  class Solver, class Preconditioner>
166 int cg_abtb (const Matrix& A, const Matrix& B, Vector& u, Vector& p,
167  const VectorExpr1& Mf, const VectorExpr2& Mg, const Preconditioner& S1,
168  const Solver& inner_solver_A, const solver_option& sopt = solver_option())
169 {
170  sopt.label = (sopt.label != "" ? sopt.label : "cg_abtb");
171  abtb_schur_complement<Matrix,Vector,Solver> S (inner_solver_A, B);
172  Vector v = inner_solver_A.solve (Mf);
173  Vector Mh = B*v - Mg;
174  int status = cg (S, p, Mh, S1, sopt);
175  u = inner_solver_A.solve (Mf - B.trans_mult(p));
176  return status;
177 }
178 template <class Matrix, class Vector, class Solver, class Preconditioner>
179 int minres_abtbc (const Matrix& A, const Matrix& B, const Matrix& C, Vector& u, Vector& p,
180  const Vector& Mf, const Vector& Mg, const Preconditioner& S1,
181  const Solver& inner_solver_A, const solver_option& sopt = solver_option())
182 {
183  sopt.label = (sopt.label != "" ? sopt.label : "minres_abtbc");
184  abtbc_schur_complement<Matrix,Vector,Solver> S (inner_solver_A, B, C);
185  Vector v = inner_solver_A.solve (Mf);
186  Vector Mh = B*v - Mg;
187  int status = minres(S, p, Mh, S1, sopt);
188  u = inner_solver_A.solve (Mf - B.trans_mult(p));
189  return status;
190 }
191 template <class Matrix, class Vector, class Solver, class Preconditioner>
192 int minres_abtb (const Matrix& A, const Matrix& B, Vector& u, Vector& p,
193  const Vector& Mf, const Vector& Mg, const Preconditioner& S1,
194  const Solver& inner_solver_A, const solver_option& sopt = solver_option())
195 {
196  sopt.label = (sopt.label != "" ? sopt.label : "minres_abtb");
197  abtb_schur_complement<Matrix,Vector,Solver> S (inner_solver_A, B);
198  Vector v = inner_solver_A.solve (Mf);
199  Vector Mh = B*v - Mg;
200  int status = minres(S, p, Mh, S1, sopt);
201  u = inner_solver_A.solve (Mf - B.trans_mult(p));
202  return status;
203 }
204 }// namespace rheolef
205 # endif // _RHEO_MIXED_SOLVER_H
rheolef::cg
int cg(const Matrix &A, Vector &x, const Vector2 &Mb, const Preconditioner &M, const solver_option &sopt=solver_option())
Definition: cg.h:94
rheolef::Vector
Definition: Vector.h:176
rheolef::abtb_schur_complement::operator*
Vector operator*(const Vector &p) const
Definition: mixed_solver.h:43
rheolef::minres_abtb
int minres_abtb(const Matrix &A, const Matrix &B, Vector &u, Vector &p, const Vector &Mf, const Vector &Mg, const Preconditioner &S1, const Solver &inner_solver_A, const solver_option &sopt=solver_option())
Definition: mixed_solver.h:192
rheolef::abtbc_schur_complement
Definition: mixed_solver.h:29
rheolef::minres
int minres(const Matrix &A, Vector &x, const Vector &Mb, const Preconditioner &M, const solver_option &sopt=solver_option())
Definition: minres.h:100
rheolef::abtbc_schur_complement::operator*
Vector operator*(const Vector &p) const
Definition: mixed_solver.h:33
p
Definition: sphere.icc:25
rheolef::abtbc_schur_complement::iA
Solver iA
Definition: mixed_solver.h:30
rheolef::cg_abtbc
int cg_abtbc(const Matrix &A, const Matrix &B, const Matrix &C, Vector &u, Vector &p, const VectorExpr1 &Mf, const VectorExpr2 &Mg, const Preconditioner &S1, const Solver &inner_solver_A, const solver_option &sopt=solver_option())
Definition: mixed_solver.h:152
rheolef::uzawa_abtb
int uzawa_abtb(const Matrix &A, const Matrix &B, Vector &u, Vector &p, const Vector &Mf, const Vector &Mg, const Preconditioner &S1, const Solver &inner_solver_A, const Real &rho, const solver_option &sopt=solver_option())
Definition: mixed_solver.h:63
rheolef::abtb_schur_complement::iA
Solver iA
Definition: mixed_solver.h:40
rheolef
This file is part of Rheolef.
Definition: compiler_eigen.h:37
rheolef::uzawa_abtbc
int uzawa_abtbc(const Matrix &A, const Matrix &B, const Matrix &C, Vector &u, Vector &p, const Vector &Mf, const Vector &Mg, const Preconditioner &S1, const Solver &inner_solver_A, const Float &rho, const solver_option &sopt=solver_option())
Definition: mixed_solver.h:49
u
Definition: leveque.h:25
Float
see the Float page for the full documentation
rheolef::abtbc_schur_complement::abtbc_schur_complement
abtbc_schur_complement(const Solver &iA1, const Matrix &B1, const Matrix &C1)
Definition: mixed_solver.h:32
rheolef::minres_abtbc
int minres_abtbc(const Matrix &A, const Matrix &B, const Matrix &C, Vector &u, Vector &p, const Vector &Mf, const Vector &Mg, const Preconditioner &S1, const Solver &inner_solver_A, const solver_option &sopt=solver_option())
Definition: mixed_solver.h:179
rheolef::abtbc_schur_complement::B
Matrix B
Definition: mixed_solver.h:31
rheolef::abtbc_schur_complement::C
Matrix C
Definition: mixed_solver.h:31
rheolef::abtb_schur_complement::abtb_schur_complement
abtb_schur_complement(const Solver &iA1, const Matrix &B1)
Definition: mixed_solver.h:42
rheolef::abtb_schur_complement
Definition: mixed_solver.h:39
mkgeo_contraction.status
status
Definition: mkgeo_contraction.sh:290
rheolef::abtb_schur_complement::B
Matrix B
Definition: mixed_solver.h:41
rheolef::solver_option
see the solver_option page for the full documentation
Definition: solver_option.h:155
rheolef::cg_abtb
int cg_abtb(const Matrix &A, const Matrix &B, Vector &u, Vector &p, const VectorExpr1 &Mf, const VectorExpr2 &Mg, const Preconditioner &S1, const Solver &inner_solver_A, const solver_option &sopt=solver_option())
Definition: mixed_solver.h:166