Rheolef  7.2
an efficient C++ finite element environment
 
Loading...
Searching...
No Matches
solver_abtb.cc
Go to the documentation of this file.
1
21#include "rheolef/solver_abtb.h"
22#include "rheolef/linalg.h"
23
24namespace rheolef {
25
26template <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}
39template <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}
58template <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}
77template <class T, class M>
78void
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}
115template <class T, class M>
116void
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}
154template <class T, class M>
155bool
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
see the csr page for the full documentation
Definition csr.h:317
void solve(const vec< T, M > &f, const vec< T, M > &g, vec< T, M > &u, vec< T, M > &p) const
see the solver_option page for the full documentation
static const long int decide
see the vec page for the full documentation
Definition vec.h:79
Expr1::float_type T
Definition field_expr.h:230
check_macro(expr1.have_homogeneous_space(Xh1), "dual(expr1,expr2); expr1 should have homogeneous space. HINT: use dual(interpolate(Xh, expr1),expr2)")
This file is part of Rheolef.
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())
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())
rheolef::std enable_if ::type dot const Expr1 expr1, const Expr2 expr2 dot(const Expr1 &expr1, const Expr2 &expr2)
csr< T, sequential > trans(const csr< T, sequential > &a)
trans(a): see the form page for the full documentation
Definition csr.h:455
field residue(Float p, const field &uh)
Definition cavity_dg.h:29
Definition cavity_dg.h:25
Definition sphere.icc:25
see the range page for the full documentation
Definition range.h:61
Definition leveque.h:25