Loading...
Searching...
No Matches
mixed_solver.h
Go to the documentation of this file.
# ifndef _RHEO_MIXED_SOLVER_H
# define _RHEO_MIXED_SOLVER_H
// This file is part of Rheolef.
//
// Copyright (C) 2000-2009 Pierre Saramito <Pierre.Saramito@imag.fr>
//
// Rheolef is free software; you can redistribute it and/or modify
// it under the terms of the GNU General Public License as published by
// the Free Software Foundation; either version 2 of the License, or
// (at your option) any later version.
//
// Rheolef is distributed in the hope that it will be useful,
// but WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
// GNU General Public License for more details.
//
// You should have received a copy of the GNU General Public License
// along with Rheolef; if not, write to the Free Software
// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
//
// =========================================================================
#include "rheolef/uzawa.h"
#include "rheolef/cg.h"
#include "rheolef/minres.h"
namespace rheolef {
template <class Matrix, class Vector, class Solver>
struct abtbc_schur_complement { // S = B*inv(A)*B^T + C
Solver iA;
Matrix B, C;
abtbc_schur_complement (const Solver& iA1,const Matrix& B1, const Matrix& C1) : iA(iA1), B(B1), C(C1) {}
Vector operator* (const Vector& p) const {
Vector v = iA.solve(B.trans_mult(p));
return B*v + C*p;
}
};
template <class Matrix, class Vector, class Solver>
struct abtb_schur_complement { // S = B*inv(A)*B^T
Solver iA;
Matrix B;
abtb_schur_complement (const Solver& iA1,const Matrix& B1) : iA(iA1), B(B1) {}
Vector operator* (const Vector& p) const {
Vector v = iA.solve(B.trans_mult(p));
return B*v;
}
};
template <class Matrix, class Vector, class Solver, class Preconditioner>
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())
{
sopt.label = (sopt.label != "" ? sopt.label : "uzawa_abtbc");
abtbc_schur_complement<Matrix,Vector,Solver> S (inner_solver_A, B, C);
Vector v = inner_solver_A.solve (Mf);
Vector Mh = B*v - Mg;
int status = uzawa (S, p, Mh, S1, sopt);
u = inner_solver_A.solve (Mf - B.trans_mult(p));
return status;
}
template <class Matrix, class Vector, class Solver, class Preconditioner, class Real>
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())
{
sopt.label = (sopt.label != "" ? sopt.label : "uzawa_abtb");
abtb_schur_complement<Matrix,Vector,Solver> S (inner_solver_A, B);
Vector v = inner_solver_A.solve (Mf);
Vector Mh = B*v - Mg;
int status = uzawa (S, p, Mh, S1, rho, sopt);
u = inner_solver_A.solve (Mf - B.trans_mult(p));
return status;
}
/*Class:mixed_solver
NAME: @code{cg_abtb}, @code{cg_abtbc}, @code{minres_abtb}, @code{minres_abtbc} -- solvers for mixed linear problems
@findex cg
@findex minres
@findex cg\_abtb
@findex cg\_abtbc
@findex minres\_abtb
@findex minres\_abtbc
@cindex mixed linear problem
@cindex conjugate gradien algorithm
@cindex finite element method
@cindex stabilized mixed finite element method
@cindex Stokes problem
@cindex incompresible elasticity
SYNOPSIS:
@example
template <class Matrix, class Vector, class Solver, class Preconditioner>
int cg_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());
template <class Matrix, class Vector, class Solver, class Preconditioner>
int cg_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());
@end example
The synopsis is the same with the minres algorithm.
EXAMPLES:
See the user's manual for practical examples for the nearly incompressible
elasticity, the Stokes and the Navier-Stokes problems.
DESCRIPTION:
@noindent
Preconditioned conjugate gradient algorithm on the pressure p applied to
the stabilized stokes problem:
@example
[ A B^T ] [ u ] [ Mf ]
[ ] [ ] = [ ]
[ B -C ] [ p ] [ Mg ]
@end example
where A is symmetric positive definite and C is symmetric positive
and semi-definite.
Such mixed linear problems appears for instance with the discretization
of Stokes problems with stabilized P1-P1 element, or with nearly
incompressible elasticity.
Formally u = inv(A)*(Mf - B^T*p) and the reduced system writes for
all non-singular matrix S1:
@example
inv(S1)*(B*inv(A)*B^T)*p = inv(S1)*(B*inv(A)*Mf - Mg)
@end example
Uzawa or conjugate gradient algorithms are considered on the
reduced problem.
Here, S1 is some preconditioner for the Schur complement S=B*inv(A)*B^T.
Both direct or iterative solvers for S1*q = t are supported.
Application of inv(A) is performed via a call to a solver
for systems such as A*v = b.
This last system may be solved either by direct or iterative algorithms,
thus, a general matrix solver class is submitted to the algorithm.
For most applications, such as the Stokes problem,
the mass matrix for the p variable is a good S1 preconditioner
for the Schur complement.
The stopping criteria is expressed using the S1 matrix, i.e. in L2 norm
when this choice is considered.
It is scaled by the L2 norm of the right-hand side of the reduced system,
also in S1 norm.
AUTHOR:
| Pierre.Saramito@imag.fr
LJK-IMAG, 38041 Grenoble cedex 9, France
DATE: 15 april 2009
METHODS: @mixed_solver
End:
*/
template <class Matrix, class Vector, class VectorExpr1, class VectorExpr2,
class Solver, class Preconditioner>
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())
{
sopt.label = (sopt.label != "" ? sopt.label : "cg_abtbc");
abtbc_schur_complement<Matrix,Vector,Solver> S (inner_solver_A, B, C);
Vector v = inner_solver_A.solve (Mf);
Vector Mh = B*v - Mg;
int status = cg (S, p, Mh, S1, sopt);
u = inner_solver_A.solve (Mf - B.trans_mult(p));
return status;
}
template <class Matrix, class Vector, class VectorExpr1, class VectorExpr2,
class Solver, class Preconditioner>
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())
{
sopt.label = (sopt.label != "" ? sopt.label : "cg_abtb");
abtb_schur_complement<Matrix,Vector,Solver> S (inner_solver_A, B);
Vector v = inner_solver_A.solve (Mf);
Vector Mh = B*v - Mg;
int status = cg (S, p, Mh, S1, sopt);
u = inner_solver_A.solve (Mf - B.trans_mult(p));
return status;
}
template <class Matrix, class Vector, class Solver, class Preconditioner>
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())
{
sopt.label = (sopt.label != "" ? sopt.label : "minres_abtbc");
abtbc_schur_complement<Matrix,Vector,Solver> S (inner_solver_A, B, C);
Vector v = inner_solver_A.solve (Mf);
Vector Mh = B*v - Mg;
int status = minres(S, p, Mh, S1, sopt);
u = inner_solver_A.solve (Mf - B.trans_mult(p));
return status;
}
template <class Matrix, class Vector, class Solver, class Preconditioner>
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())
{
sopt.label = (sopt.label != "" ? sopt.label : "minres_abtb");
abtb_schur_complement<Matrix,Vector,Solver> S (inner_solver_A, B);
Vector v = inner_solver_A.solve (Mf);
Vector Mh = B*v - Mg;
int status = minres(S, p, Mh, S1, sopt);
u = inner_solver_A.solve (Mf - B.trans_mult(p));
return status;
}
}// namespace rheolef
# endif // _RHEO_MIXED_SOLVER_H
Float
see the Float page for the full documentation
rheolef::solver_option
see the solver_option page for the full documentation
Definition solver_option.h:155
rheolef
This file is part of Rheolef.
Definition compiler_eigen.h:39
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::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
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::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::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::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::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::uzawa
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
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
A
Definition diffusion_isotropic.h:25
p
Definition sphere.icc:25
rheolef::Vector
Definition Vector.h:176
rheolef::abtb_schur_complement
Definition mixed_solver.h:39
rheolef::abtb_schur_complement::B
Matrix B
Definition mixed_solver.h:41
rheolef::abtb_schur_complement::operator*
Vector operator*(const Vector &p) const
Definition mixed_solver.h:43
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::iA
Solver iA
Definition mixed_solver.h:40
rheolef::abtbc_schur_complement
Definition mixed_solver.h:29
rheolef::abtbc_schur_complement::B
Matrix B
Definition mixed_solver.h:31
rheolef::abtbc_schur_complement::operator*
Vector operator*(const Vector &p) const
Definition mixed_solver.h:33
rheolef::abtbc_schur_complement::C
Matrix C
Definition mixed_solver.h:31
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::abtbc_schur_complement::iA
Solver iA
Definition mixed_solver.h:30
u
Definition leveque.h:25