Loading...
Searching...
No Matches
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
f
point f(const Float &u)
Definition burgers.icc:25
g
u_exact g
Definition burgers_diffusion_exact.h:33
rheolef::csr
see the csr page for the full documentation
Definition csr.h:317
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::_a
csr< T, M > _a
Definition solver_abtb.h:147
rheolef::solver_abtb_basic::_b
csr< T, M > _b
Definition solver_abtb.h:148
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::_sA
solver_basic< T, M > _sA
Definition solver_abtb.h:151
rheolef::solver_abtb_basic::_opt
solver_option _opt
Definition solver_abtb.h:146
rheolef::solver_abtb_basic::solver_abtb_basic
solver_abtb_basic()
Definition solver_abtb.cc:27
rheolef::solver_abtb_basic::_need_constraint
bool _need_constraint
Definition solver_abtb.h:154
rheolef::solver_abtb_basic::_sa
solver_basic< T, M > _sa
Definition solver_abtb.h:152
rheolef::solver_abtb_basic< T, rheo_default_memory_model >::option
const solver_option & option() const
Definition solver_abtb.h:134
rheolef::solver_abtb_basic::_c
csr< T, M > _c
Definition solver_abtb.h:149
rheolef::solver_abtb_basic::_smp
solver_basic< T, M > _smp
Definition solver_abtb.h:153
rheolef::solver_abtb_basic::initialized
bool initialized() const
Definition solver_abtb.cc:156
rheolef::solver_basic
Definition solver.h:264
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
see the vec page for the full documentation
Definition vec.h:79
u
point u(const point &x)
Definition diffusion_transport_tensor_exact.icc:25
T
Expr1::float_type T
Definition field_expr.h:230
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_obstacle.L
str L
Definition mkgeo_obstacle.sh:133
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::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::dot
rheolef::std::enable_if< details::is_vec_expr_v2_arg< Expr1 >::value &&details::is_vec_expr_v2_arg< Expr2 >::value, promote< Expr1::float_type, Expr2::float_type >::type > type dot const Expr1 expr1, const Expr2 expr2 dot(const Expr1 &expr1, const Expr2 &expr2)
Definition vec_expr_v2.h:415
d
size_t d
Definition neumann-laplace-lambda.cc:28
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