solver_abtb.cc
Go to the documentation of this file.
#include "rheolef/solver_abtb.h"
#include "rheolef/linalg.h"
namespace rheolef {
template <class T, class M>
solver_abtb_basic<T,M>::solver_abtb_basic ()
: _opt(),
_a(),
_b(),
_c(),
_mp(),
_sA(),
_sa(),
_smp(),
_need_constraint(false)
{
}
template <class T, class M>
solver_abtb_basic<T,M>::solver_abtb_basic(
const csr<T,M>& a,
const csr<T,M>& b,
const csr<T,M>& mp,
const solver_option& opt)
: _opt(opt),
_a(a),
_b(b),
_c(),
_mp(mp),
_sA(),
_sa(),
_smp(),
_need_constraint(false)
{
_c.resize (_mp.row_ownership(), _mp.col_ownership());
init();
}
template <class T, class M>
solver_abtb_basic<T,M>::solver_abtb_basic(
const csr<T,M>& a,
const csr<T,M>& b,
const csr<T,M>& c,
const csr<T,M>& mp,
const solver_option& opt)
: _opt(opt),
_a(a),
_b(b),
_c(c),
_mp(mp),
_sA(),
_sa(),
_smp(),
_need_constraint(false)
{
init();
}
template <class T, class M>
void
solver_abtb_basic<T,M>::init()
{
if (_opt.iterative == solver_option::decide) {
_opt.iterative = (_a.pattern_dimension() > 2);
}
if (! _opt.iterative) {
// direct stokes solver:
csr<T,M> A;
vec<T,M> one (_b.row_ownership(), 1);
vec<T,M> r =_b.trans_mult (one);
T z = dot(r,r);
if (_c.dis_nnz() != 0) {
vec<T,M> r2 =_c*one;
z += dot(r2,r2);
}
_need_constraint = (fabs(z) <= std::numeric_limits<T>::epsilon());
if (_need_constraint) {
vec<T,M> zu (_a.col_ownership(), 0);
vec<T,M> d = _mp*one;
A = { { _a, trans(_b), zu },
{ _b, - _c, d },
{ trans(zu), trans(d), T(0)} };
} else {
A = { { _a, trans(_b)},
{ _b, - _c } };
}
A.set_pattern_dimension (_a.pattern_dimension());
A.set_definite_positive (false); // A has both > 0 and < 0 eigenvalues : used by solver
if (_c.dis_nnz() == 0) {
A.set_symmetry (_a.is_symmetric());
} else {
A.set_symmetry (_a.is_symmetric() && _c.is_symmetric());
}
_sA = solver_basic<T,M>(A,_opt);
}
}
template <class T, class M>
void
solver_abtb_basic<T,M>::solve (const vec<T,M>& f, const vec<T,M>& g, vec<T,M>& u, vec<T,M>& p) const
{
if (_opt.iterative) {
// if preconditioner and inner solver are not set, use the default:
if (!_smp.initialized()) { _smp = solver_basic<T,M>(_mp,_opt); }
if (!_sa.initialized()) { _sa = solver_basic<T,M>(_a, _opt); }
// iterative stokes solver: mp is the preconditioner
check_macro (_a.is_symmetric(), "solver_abtb: iterative for unsymmetric system: not yet (HINT: use direct solver)");
if (_c.dis_nnz() == 0) {
cg_abtb (_a, _b, u, p, f, g, _smp, _sa, option());
} else {
cg_abtbc (_a, _b, _c, u, p, f, g, _smp, _sa, option());
}
check_macro (option().residue <= option().tol,
"solver: precision "<<option().tol<<" not reached: get "<<option().residue
<< " after " << option().max_iter << " iterations");
return;
}
// direct stokes solver:
vec<T,M> L;
if (_need_constraint) {
L = { f, g, T(0) };
} else {
L = { f, g };
}
vec<T,M> U = _sA.solve (L);
u = U [range(0,u.size())];
p = U [range(u.size(),u.size()+p.size())];
if (_need_constraint) {
// lambda no more used:
T lambda = (U.size() == u.size()+p.size()+1) ? U [u.size()+p.size()] : 0;
#ifdef _RHEOLEF_HAVE_MPI
mpi::broadcast (U.comm(), lambda, U.comm().size() - 1);
#endif // _RHEOLEF_HAVE_MPI
}
}
template <class T, class M>
bool
solver_abtb_basic<T,M>::initialized() const
{
return (_a.dis_nnz() != 0);
}
// ----------------------------------------------------------------------------
// instanciation in library
// ----------------------------------------------------------------------------
template class solver_abtb_basic<Float,sequential>;
#ifdef _RHEOLEF_HAVE_MPI
template class solver_abtb_basic<Float,distributed>;
#endif // _RHEOLEF_HAVE_MPI
} // namespace rheolef
g
u_exact g
Definition: burgers_diffusion_exact.h:33
rheolef::csr< T, M >
rheolef::solver_abtb_basic
Definition: solver_abtb.h:115
rheolef::solver_abtb_basic::_mp
csr< T, M > _mp
Definition: solver_abtb.h:150
rheolef::solver_abtb_basic::init
void init()
Definition: solver_abtb.cc:79
rheolef::solver_abtb_basic::solve
void solve(const vec< T, M > &f, const vec< T, M > &g, vec< T, M > &u, vec< T, M > &p) const
Definition: solver_abtb.cc:117
rheolef::solver_abtb_basic::solver_abtb_basic
solver_abtb_basic()
Definition: solver_abtb.cc:27
rheolef::solver_abtb_basic::_c
csr< T, M > _c
Definition: solver_abtb.h:149
rheolef::solver_abtb_basic::initialized
bool initialized() const
Definition: solver_abtb.cc:156
rheolef::solver_basic< T, M >
rheolef::solver_option
see the solver_option page for the full documentation
Definition: solver_option.h:155
rheolef::solver_option::decide
static const long int decide
Definition: solver_option.h:158
rheolef::vec< T, M >
T
Expr1::float_type T
Definition: field_expr.h:261
check_macro
check_macro(expr1.have_homogeneous_space(Xh1), "dual(expr1,expr2); expr1 should have homogeneous space. HINT: use dual(interpolate(Xh, expr1),expr2)")
mkgeo_ball.b
int b
Definition: mkgeo_ball.sh:152
mkgeo_ball.f
int f
Definition: mkgeo_ball.sh:221
mkgeo_ball.c
int c
Definition: mkgeo_ball.sh:153
mkgeo_ball.d
int d
Definition: mkgeo_ball.sh:154
mkgeo_obstacle.L
string L
Definition: mkgeo_obstacle.sh:133
rheolef
This file is part of Rheolef.
Definition: compiler_eigen.h:37
rheolef::trans
csr< T, sequential > trans(const csr< T, sequential > &a)
trans(a): see the form page for the full documentation
Definition: csr.h:455
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::dot
rheolef::std enable_if ::type dot const Expr1 expr1, const Expr2 expr2 dot(const Expr1 &expr1, const Expr2 &expr2)
dot(x,y): see the expression page for the full documentation
Definition: vec_expr_v2.h:415
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
residue
field residue(Float p, const field &uh)
Definition: p_laplacian_post.cc:35
a
Definition: diffusion_isotropic.h:25
f
Definition: cavity_dg.h:29
g
Definition: cavity_dg.h:25
lambda
Definition: yield_slip_circle.h:34
p
Definition: sphere.icc:25
rheolef::range
see the range page for the full documentation
Definition: range.h:61
u
Definition: leveque.h:25
epsilon
Float epsilon
Definition: transmission_error.cc:25