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