Rheolef  7.2
an efficient C++ finite element environment
 
Loading...
Searching...
No Matches
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"
26namespace rheolef {
27
28template <class Matrix, class Vector, class Solver>
29struct 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};
38template <class Matrix, class Vector, class Solver>
39struct 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};
48template <class Matrix, class Vector, class Solver, class Preconditioner>
49int 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}
62template <class Matrix, class Vector, class Solver, class Preconditioner, class Real>
63int 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");
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
77NAME: @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
90SYNOPSIS:
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
104EXAMPLES:
105 See the user's manual for practical examples for the nearly incompressible
106 elasticity, the Stokes and the Navier-Stokes problems.
107
108DESCRIPTION:
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.
142AUTHOR:
143 | Pierre.Saramito@imag.fr
144 LJK-IMAG, 38041 Grenoble cedex 9, France
145DATE: 15 april 2009
146METHODS: @mixed_solver
147End:
148*/
149
150template <class Matrix, class Vector, class VectorExpr1, class VectorExpr2,
151 class Solver, class Preconditioner>
152int 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}
164template <class Matrix, class Vector, class VectorExpr1, class VectorExpr2,
165 class Solver, class Preconditioner>
166int 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");
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}
178template <class Matrix, class Vector, class Solver, class Preconditioner>
179int 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}
191template <class Matrix, class Vector, class Solver, class Preconditioner>
192int 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");
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
see the Float page for the full documentation
see the solver_option page for the full documentation
This file is part of Rheolef.
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())
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())
int minres(const Matrix &A, Vector &x, const Vector &Mb, const Preconditioner &M, const solver_option &sopt=solver_option())
Definition minres.h:100
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())
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())
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())
int cg(const Matrix &A, Vector &x, const Vector2 &Mb, const Preconditioner &M, const solver_option &sopt=solver_option())
Definition cg.h:94
int uzawa(const Matrix &A, Vector &x, const Vector &Mb, const Preconditioner &M, const Real2 &rho, const solver_option &sopt=solver_option())
Definition uzawa.h:85
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 sphere.icc:25
Vector operator*(const Vector &p) const
abtb_schur_complement(const Solver &iA1, const Matrix &B1)
Vector operator*(const Vector &p) const
abtbc_schur_complement(const Solver &iA1, const Matrix &B1, const Matrix &C1)
Definition leveque.h:25