dune-pdelab 2.7-git
Loading...
Searching...
No Matches
linearproblem.hh
Go to the documentation of this file.
1#ifndef DUNE_PDELAB_STATIONARY_LINEARPROBLEM_HH
2#define DUNE_PDELAB_STATIONARY_LINEARPROBLEM_HH
3
4#include <iostream>
5
6#include <dune/common/timer.hh>
7#include <dune/common/parametertree.hh>
8
12
13namespace Dune {
14 namespace PDELab {
15
16 //===============================================================
17 // A class for solving linear stationary problems.
18 // It assembles the matrix, computes the right hand side and
19 // solves the problem.
20 // This is only a first vanilla implementation which has to be improved.
21 //===============================================================
22
23 // Status information of linear problem solver
24 template<class RFType>
26 {
27 RFType first_defect; // the first defect
28 RFType defect; // the final defect
29 double assembler_time; // Cumulative time for matrix assembly
30 double linear_solver_time; // Cumulative time for linear solver
31 int linear_solver_iterations; // Total number of linear iterations
32
40
41 };
42
43 template<typename GO, typename LS, typename V>
45 {
46 typedef typename Dune::template FieldTraits<typename V::ElementType >::real_type Real;
47 typedef typename GO::Traits::Jacobian M;
48 typedef typename GO::Traits::TrialGridFunctionSpace TrialGridFunctionSpace;
50 typedef GO GridOperator;
51
52 public:
54
55 StationaryLinearProblemSolver(const GO& go, LS& ls, V& x, Real reduction, Real min_defect = 1e-99, int verbose=1)
56 : _go(go)
57 , _ls(ls)
58 , _x(&x)
59 , _reduction(reduction)
60 , _min_defect(min_defect)
61 , _hanging_node_modifications(false)
62 , _keep_matrix(true)
63 , _verbose(verbose)
64 {}
65
66 StationaryLinearProblemSolver (const GO& go, LS& ls, Real reduction, Real min_defect = 1e-99, int verbose=1)
67 : _go(go)
68 , _ls(ls)
69 , _x()
70 , _reduction(reduction)
71 , _min_defect(min_defect)
72 , _hanging_node_modifications(false)
73 , _keep_matrix(true)
74 , _verbose(verbose)
75 {}
76
78
95 StationaryLinearProblemSolver(const GO& go, LS& ls, V& x, const ParameterTree& params)
96 : _go(go)
97 , _ls(ls)
98 , _x(&x)
99 , _reduction(params.get<Real>("reduction"))
100 , _min_defect(params.get<Real>("min_defect",1e-99))
101 , _hanging_node_modifications(params.get<bool>("hanging_node_modifications",false))
102 , _keep_matrix(params.get<bool>("keep_matrix",true))
103 , _verbose(params.get<int>("verbosity",1))
104 {}
105
107
124 StationaryLinearProblemSolver(const GO& go, LS& ls, const ParameterTree& params)
125 : _go(go)
126 , _ls(ls)
127 , _x()
128 , _reduction(params.get<Real>("reduction"))
129 , _min_defect(params.get<Real>("min_defect",1e-99))
130 , _hanging_node_modifications(params.get<bool>("hanging_node_modifications",false))
131 , _keep_matrix(params.get<bool>("keep_matrix",true))
132 , _verbose(params.get<int>("verbosity",1))
133 {}
134
137 {
138 _hanging_node_modifications=b;
139 }
140
143 {
144 return _hanging_node_modifications;
145 }
146
148 void setKeepMatrix(bool b)
149 {
150 _keep_matrix = b;
151 }
152
154 bool keepMatrix() const
155 {
156 return _keep_matrix;
157 }
158
159 const Result& result() const
160 {
161 return _res;
162 }
163
164 void apply(V& x, bool reuse_matrix = false) {
165 _x = &x;
166 apply(reuse_matrix);
167 }
168
169 void apply (bool reuse_matrix = false)
170 {
171 Dune::Timer watch;
172 double timing,assembler_time=0;
173
174 // assemble matrix; optional: assemble only on demand!
175 watch.reset();
176
177 if (!_jacobian)
178 {
179 _jacobian = std::make_shared<M>(_go);
180 timing = watch.elapsed();
181 if (_go.trialGridFunctionSpace().gridView().comm().rank()==0 && _verbose>=1)
182 std::cout << "=== matrix setup (max) " << timing << " s" << std::endl;
183 watch.reset();
184 assembler_time += timing;
185 }
186 else if (_go.trialGridFunctionSpace().gridView().comm().rank()==0 && _verbose>=1)
187 std::cout << "=== matrix setup skipped (matrix already allocated)" << std::endl;
188
189 if (_hanging_node_modifications)
190 {
191 Dune::PDELab::set_shifted_dofs(_go.localAssembler().trialConstraints(),0.0,*_x); // set hanging node DOFs to zero
192 _go.localAssembler().backtransform(*_x); // interpolate hanging nodes adjacent to Dirichlet nodes
193 }
194
195 if (!reuse_matrix)
196 {
197 (*_jacobian) = Real(0.0);
198 _go.jacobian(*_x,*_jacobian);
199 }
200
201 timing = watch.elapsed();
202 // timing = gos.trialGridFunctionSpace().gridView().comm().max(timing);
203 if (_go.trialGridFunctionSpace().gridView().comm().rank()==0 && _verbose>=1)
204 {
205 if (reuse_matrix)
206 std::cout << "=== matrix assembly SKIPPED" << std::endl;
207 else
208 std::cout << "=== matrix assembly (max) " << timing << " s" << std::endl;
209 }
210
211 assembler_time += timing;
212
213 // assemble residual
214 watch.reset();
215
216 W r(_go.testGridFunctionSpace(),0.0);
217 _go.residual(*_x,r); // residual is additive
218
219 timing = watch.elapsed();
220 // timing = gos.trialGridFunctionSpace().gridView().comm().max(timing);
221 if (_go.trialGridFunctionSpace().gridView().comm().rank()==0 && _verbose>=1)
222 std::cout << "=== residual assembly (max) " << timing << " s" << std::endl;
223 assembler_time += timing;
224 _res.assembler_time = assembler_time;
225
226 auto defect = _ls.norm(r);
227
228 // compute correction
229 watch.reset();
230 V z(_go.trialGridFunctionSpace(),0.0);
231 auto red = std::max(_reduction,_min_defect/defect);
232 if (_go.trialGridFunctionSpace().gridView().comm().rank()==0 && _verbose>=1)
233 {
234 std::cout << "=== solving (reduction: " << red << ") ";
235 if (_verbose>=1)
236 std::cout << std::flush;
237 else
238 std::cout << std::endl;
239 }
240 _ls.apply(*_jacobian,z,r,red); // solver makes right hand side consistent
241 _linear_solver_result = _ls.result();
242 timing = watch.elapsed();
243 // timing = gos.trialGridFunctionSpace().gridView().comm().max(timing);
244 if (_go.trialGridFunctionSpace().gridView().comm().rank()==0 && _verbose>=1)
245 std::cout << timing << " s" << std::endl;
246 _res.linear_solver_time = timing;
247
248 _res.converged = _linear_solver_result.converged;
249 _res.iterations = _linear_solver_result.iterations;
250 _res.elapsed = _linear_solver_result.elapsed;
251 _res.reduction = _linear_solver_result.reduction;
252 _res.conv_rate = _linear_solver_result.conv_rate;
253 _res.first_defect = static_cast<double>(defect);
254 _res.defect = static_cast<double>(defect)*_linear_solver_result.reduction;
255 _res.linear_solver_iterations = _linear_solver_result.iterations;
256
257 // and update
258 if (_hanging_node_modifications)
259 Dune::PDELab::set_shifted_dofs(_go.localAssembler().trialConstraints(),0.0,*_x); // set hanging node DOFs to zero
260 *_x -= z;
261 if (_hanging_node_modifications)
262 _go.localAssembler().backtransform(*_x); // interpolate hanging nodes adjacent to Dirichlet nodes
263
264 if (!_keep_matrix)
265 _jacobian.reset();
266 }
267
270 {
271 if(_jacobian)
272 _jacobian.reset();
273 }
274
276 return _linear_solver_result;
277 }
278
279 Real reduction() const
280 {
281 return _reduction;
282 }
283
285 {
286 _reduction = reduction;
287 }
288
289
290 private:
291 const GO& _go;
292 LS& _ls;
293 V* _x;
294 std::shared_ptr<M> _jacobian;
295 Real _reduction;
296 Real _min_defect;
297 Dune::PDELab::LinearSolverResult<double> _linear_solver_result;
298 Result _res;
299 bool _hanging_node_modifications;
300 bool _keep_matrix;
301 int _verbose;
302 };
303
304 } // namespace PDELab
305} // namespace Dune
306
307#endif // DUNE_PDELAB_STATIONARY_LINEARPROBLEM_HH
void set_shifted_dofs(const CG &cg, typename XG::ElementType x, XG &xg)
Definition constraints.hh:1014
For backward compatibility – Do not use this!
Definition adaptivity.hh:28
typename impl::BackendVectorSelector< GridFunctionSpace, FieldType >::Type Vector
alias of the return type of BackendVectorSelector
Definition backend/interface.hh:106
Definition solver.hh:31
RFType conv_rate
Definition solver.hh:36
bool converged
Definition solver.hh:32
RFType reduction
Definition solver.hh:35
unsigned int iterations
Definition solver.hh:33
double elapsed
Definition solver.hh:34
int linear_solver_iterations
Definition linearproblem.hh:31
double linear_solver_time
Definition linearproblem.hh:30
StationaryLinearProblemSolverResult()
Definition linearproblem.hh:33
RFType defect
Definition linearproblem.hh:28
double assembler_time
Definition linearproblem.hh:29
RFType first_defect
Definition linearproblem.hh:27
Definition linearproblem.hh:45
const Result & result() const
Definition linearproblem.hh:159
void setKeepMatrix(bool b)
Set whether the jacobian matrix should be kept across calls to apply().
Definition linearproblem.hh:148
Real reduction() const
Definition linearproblem.hh:279
void apply(bool reuse_matrix=false)
Definition linearproblem.hh:169
StationaryLinearProblemSolver(const GO &go, LS &ls, Real reduction, Real min_defect=1e-99, int verbose=1)
Definition linearproblem.hh:66
void setHangingNodeModifications(bool b)
Set whether the solver should apply the necessary transformations for calculations on hanging nodes.
Definition linearproblem.hh:136
bool hangingNodeModifications() const
Return whether the solver performs the necessary transformations for calculations on hanging nodes.
Definition linearproblem.hh:142
StationaryLinearProblemSolver(const GO &go, LS &ls, V &x, const ParameterTree &params)
Construct a StationaryLinearProblemSolver for the given objects and read parameters from a ParameterT...
Definition linearproblem.hh:95
StationaryLinearProblemSolver(const GO &go, LS &ls, V &x, Real reduction, Real min_defect=1e-99, int verbose=1)
Definition linearproblem.hh:55
StationaryLinearProblemSolverResult< double > Result
Definition linearproblem.hh:53
void apply(V &x, bool reuse_matrix=false)
Definition linearproblem.hh:164
const Dune::PDELab::LinearSolverResult< double > & ls_result() const
Definition linearproblem.hh:275
void discardMatrix()
Discard the stored Jacobian matrix.
Definition linearproblem.hh:269
void setReduction(Real reduction)
Definition linearproblem.hh:284
bool keepMatrix() const
Return whether the jacobian matrix is kept across calls to apply().
Definition linearproblem.hh:154
StationaryLinearProblemSolver(const GO &go, LS &ls, const ParameterTree &params)
Construct a StationaryLinearProblemSolver for the given objects and read parameters from a ParameterT...
Definition linearproblem.hh:124