Rheolef  7.1
an efficient C++ finite element environment
solver_abtb.cc
Go to the documentation of this file.
1 #include "rheolef/solver_abtb.h"
22 #include "rheolef/linalg.h"
23 
24 namespace rheolef {
25 
26 template <class T, class M>
28  : _opt(),
29  _a(),
30  _b(),
31  _c(),
32  _mp(),
33  _sA(),
34  _sa(),
35  _smp(),
36  _need_constraint(false)
37 {
38 }
39 template <class T, class M>
41  const csr<T,M>& a,
42  const csr<T,M>& b,
43  const csr<T,M>& mp,
44  const solver_option& opt)
45  : _opt(opt),
46  _a(a),
47  _b(b),
48  _c(),
49  _mp(mp),
50  _sA(),
51  _sa(),
52  _smp(),
53  _need_constraint(false)
54 {
55  _c.resize (_mp.row_ownership(), _mp.col_ownership());
56  init();
57 }
58 template <class T, class M>
60  const csr<T,M>& a,
61  const csr<T,M>& b,
62  const csr<T,M>& c,
63  const csr<T,M>& mp,
64  const solver_option& opt)
65  : _opt(opt),
66  _a(a),
67  _b(b),
68  _c(c),
69  _mp(mp),
70  _sA(),
71  _sa(),
72  _smp(),
73  _need_constraint(false)
74 {
75  init();
76 }
77 template <class T, class M>
78 void
80 {
81  if (_opt.iterative == solver_option::decide) {
82  _opt.iterative = (_a.pattern_dimension() > 2);
83  }
84  if (! _opt.iterative) {
85  // direct stokes solver:
86  csr<T,M> A;
87  vec<T,M> one (_b.row_ownership(), 1);
88  vec<T,M> r =_b.trans_mult (one);
89  T z = dot(r,r);
90  if (_c.dis_nnz() != 0) {
91  vec<T,M> r2 =_c*one;
92  z += dot(r2,r2);
93  }
94  _need_constraint = (fabs(z) <= std::numeric_limits<T>::epsilon());
95  if (_need_constraint) {
96  vec<T,M> zu (_a.col_ownership(), 0);
97  vec<T,M> d = _mp*one;
98  A = { { _a, trans(_b), zu },
99  { _b, - _c, d },
100  { trans(zu), trans(d), T(0)} };
101  } else {
102  A = { { _a, trans(_b)},
103  { _b, - _c } };
104  }
105  A.set_pattern_dimension (_a.pattern_dimension());
106  A.set_definite_positive (false); // A has both > 0 and < 0 eigenvalues : used by solver
107  if (_c.dis_nnz() == 0) {
108  A.set_symmetry (_a.is_symmetric());
109  } else {
110  A.set_symmetry (_a.is_symmetric() && _c.is_symmetric());
111  }
112  _sA = solver_basic<T,M>(A,_opt);
113  }
114 }
115 template <class T, class M>
116 void
118 {
119  if (_opt.iterative) {
120  // if preconditioner and inner solver are not set, use the default:
121  if (!_smp.initialized()) { _smp = solver_basic<T,M>(_mp,_opt); }
122  if (!_sa.initialized()) { _sa = solver_basic<T,M>(_a, _opt); }
123  // iterative stokes solver: mp is the preconditioner
124  check_macro (_a.is_symmetric(), "solver_abtb: iterative for unsymmetric system: not yet (HINT: use direct solver)");
125  if (_c.dis_nnz() == 0) {
126  cg_abtb (_a, _b, u, p, f, g, _smp, _sa, option());
127  } else {
128  cg_abtbc (_a, _b, _c, u, p, f, g, _smp, _sa, option());
129  }
130  check_macro (option().residue <= option().tol,
131  "solver: precision "<<option().tol<<" not reached: get "<<option().residue
132  << " after " << option().max_iter << " iterations");
133  return;
134  }
135  // direct stokes solver:
136  vec<T,M> L;
137  if (_need_constraint) {
138  L = { f, g, T(0) };
139  } else {
140  L = { f, g };
141  }
142  vec<T,M> U = _sA.solve (L);
143  u = U [range(0,u.size())];
144  p = U [range(u.size(),u.size()+p.size())];
145 
146  if (_need_constraint) {
147  // lambda no more used:
148  T lambda = (U.size() == u.size()+p.size()+1) ? U [u.size()+p.size()] : 0;
149 #ifdef _RHEOLEF_HAVE_MPI
150  mpi::broadcast (U.comm(), lambda, U.comm().size() - 1);
151 #endif // _RHEOLEF_HAVE_MPI
152  }
153 }
154 template <class T, class M>
155 bool
157 {
158  return (_a.dis_nnz() != 0);
159 }
160 // ----------------------------------------------------------------------------
161 // instanciation in library
162 // ----------------------------------------------------------------------------
164 
165 #ifdef _RHEOLEF_HAVE_MPI
167 #endif // _RHEOLEF_HAVE_MPI
168 
169 } // namespace rheolef
g
u_exact g
Definition: burgers_diffusion_exact.h:33
check_macro
check_macro(expr1.have_homogeneous_space(Xh1), "dual(expr1,expr2); expr1 should have homogeneous space. HINT: use dual(interpolate(Xh, expr1),expr2)")
residue
field residue(Float p, const field &uh)
Definition: p_laplacian_post.cc:35
rheolef::range
see the range page for the full documentation
Definition: range.h:61
rheolef::solver_basic
Definition: solver.h:188
mkgeo_ball.f
f
Definition: mkgeo_ball.sh:221
mkgeo_ball.c
c
Definition: mkgeo_ball.sh:153
rheolef::solver_abtb_basic
Definition: solver_abtb.h:115
rheolef::vec
see the vec page for the full documentation
Definition: vec.h:79
rheolef::solver_abtb_basic::solver_abtb_basic
solver_abtb_basic()
Definition: solver_abtb.cc:27
p
Definition: sphere.icc:25
a
Definition: diffusion_isotropic.h:25
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::csr< T, M >
rheolef
This file is part of Rheolef.
Definition: compiler_eigen.h:37
rheolef::solve
void solve(tiny_matrix< T > &a, tiny_vector< size_t > &piv, const tiny_vector< T > &b, tiny_vector< T > &x)
Definition: tiny_lu.h:92
u
Definition: leveque.h:25
mkgeo_obstacle.L
L
Definition: mkgeo_obstacle.sh:133
mkgeo_ball.d
d
Definition: mkgeo_ball.sh:154
mkgeo_ball.b
b
Definition: mkgeo_ball.sh:152
g
Definition: cavity_dg.h:25
epsilon
Float epsilon
Definition: transmission_error.cc:25
f
Definition: cavity_dg.h:29
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::details::dot
rheolef::details::is_vec dot
T
Expr1::float_type T
Definition: field_expr.h:218
lambda
Definition: yield_slip_circle.h:34
rheolef::solver_option
see the solver_option page for the full documentation
Definition: solver_option.h:155
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