1 #ifndef DUNE_PDELAB_STATIONARYLINEARPROBLEM_HH 2 #define DUNE_PDELAB_STATIONARYLINEARPROBLEM_HH 6 #include <dune/common/timer.hh> 7 #include <dune/common/parametertree.hh> 24 template<
class RFType>
37 , linear_solver_time(0.0)
38 , linear_solver_iterations(0)
43 template<
typename GO,
typename LS,
typename V>
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;
59 , _reduction(reduction)
60 , _min_defect(min_defect)
61 , _hanging_node_modifications(false)
70 , _reduction(reduction)
71 , _min_defect(min_defect)
72 , _hanging_node_modifications(false)
99 , _reduction(params.get<Real>(
"reduction"))
100 , _min_defect(params.get<Real>(
"min_defect",1
e-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))
128 , _reduction(params.get<Real>(
"reduction"))
129 , _min_defect(params.get<Real>(
"min_defect",1
e-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))
138 _hanging_node_modifications=b;
144 return _hanging_node_modifications;
164 void apply(V& x,
bool reuse_matrix =
false) {
169 void apply (
bool reuse_matrix =
false)
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;
184 assembler_time += timing;
186 else if (_go.trialGridFunctionSpace().gridView().comm().rank()==0 && _verbose>=1)
187 std::cout <<
"=== matrix setup skipped (matrix already allocated)" << std::endl;
189 if (_hanging_node_modifications)
192 _go.localAssembler().backtransform(*_x);
197 (*_jacobian) = Real(0.0);
198 _go.jacobian(*_x,*_jacobian);
201 timing = watch.elapsed();
203 if (_go.trialGridFunctionSpace().gridView().comm().rank()==0 && _verbose>=1)
206 std::cout <<
"=== matrix assembly SKIPPED" << std::endl;
208 std::cout <<
"=== matrix assembly (max) " << timing <<
" s" << std::endl;
211 assembler_time += timing;
216 W r(_go.testGridFunctionSpace(),0.0);
219 timing = watch.elapsed();
221 if (_go.trialGridFunctionSpace().gridView().comm().rank()==0 && _verbose>=1)
222 std::cout <<
"=== residual assembly (max) " << timing <<
" s" << std::endl;
223 assembler_time += timing;
226 auto defect = _ls.norm(r);
230 V z(_go.trialGridFunctionSpace(),0.0);
231 auto red = std::max(_reduction,_min_defect/
defect);
232 if (_go.trialGridFunctionSpace().gridView().comm().rank()==0)
233 std::cout <<
"=== solving (reduction: " << red <<
") " << std::endl;
234 _ls.apply(*_jacobian,z,r,red);
235 _linear_solver_result = _ls.result();
236 timing = watch.elapsed();
238 if (_go.trialGridFunctionSpace().gridView().comm().rank()==0 && _verbose>=1)
239 std::cout << timing <<
" s" << std::endl;
240 _res.linear_solver_time = timing;
242 _res.converged = _linear_solver_result.converged;
243 _res.iterations = _linear_solver_result.iterations;
244 _res.elapsed = _linear_solver_result.elapsed;
245 _res.reduction = _linear_solver_result.reduction;
246 _res.conv_rate = _linear_solver_result.conv_rate;
247 _res.first_defect =
static_cast<double>(
defect);
248 _res.defect =
static_cast<double>(
defect)*_linear_solver_result.reduction;
249 _res.linear_solver_iterations = _linear_solver_result.iterations;
252 if (_hanging_node_modifications)
255 if (_hanging_node_modifications)
256 _go.localAssembler().backtransform(*_x);
270 return _linear_solver_result;
288 shared_ptr<M> _jacobian;
293 bool _hanging_node_modifications;
bool hangingNodeModifications() const
Return whether the solver performs the necessary transformations for calculations on hanging nodes...
Definition: linearproblem.hh:142
StationaryLinearProblemSolverResult< double > Result
Definition: linearproblem.hh:53
Real reduction() const
Definition: linearproblem.hh:273
int linear_solver_iterations
Definition: linearproblem.hh:31
StationaryLinearProblemSolverResult()
Definition: linearproblem.hh:33
void setHangingNodeModifications(bool b)
Set whether the solver should apply the necessary transformations for calculations on hanging nodes...
Definition: linearproblem.hh:136
StationaryLinearProblemSolver(const GO &go, LS &ls, Real reduction, Real min_defect=1e-99, int verbose=1)
Definition: linearproblem.hh:66
Definition: linearproblem.hh:44
Definition: adaptivity.hh:27
void setReduction(Real reduction)
Definition: linearproblem.hh:278
void setKeepMatrix(bool b)
Set whether the jacobian matrix should be kept across calls to apply().
Definition: linearproblem.hh:148
void discardMatrix()
Discard the stored Jacobian matrix.
Definition: linearproblem.hh:263
void apply(bool reuse_matrix=false)
Definition: linearproblem.hh:169
const Entity & e
Definition: localfunctionspace.hh:111
RFType defect
Definition: linearproblem.hh:28
double linear_solver_time
Definition: linearproblem.hh:30
RFType reduction
Definition: solver.hh:35
StationaryLinearProblemSolver(const GO &go, LS &ls, V &x, const ParameterTree ¶ms)
Construct a StationaryLinearProblemSolver for the given objects and read parameters from a ParameterT...
Definition: linearproblem.hh:95
RFType first_defect
Definition: linearproblem.hh:27
double assembler_time
Definition: linearproblem.hh:29
StationaryLinearProblemSolver(const GO &go, LS &ls, const ParameterTree ¶ms)
Construct a StationaryLinearProblemSolver for the given objects and read parameters from a ParameterT...
Definition: linearproblem.hh:124
const Result & result() const
Definition: linearproblem.hh:159
typename impl::BackendVectorSelector< GridFunctionSpace, FieldType >::Type Vector
alias of the return type of BackendVectorSelector
Definition: backend/interface.hh:113
const Dune::PDELab::LinearSolverResult< double > & ls_result() const
Definition: linearproblem.hh:269
void apply(V &x, bool reuse_matrix=false)
Definition: linearproblem.hh:164
void set_shifted_dofs(const CG &cg, typename XG::ElementType x, XG &xg)
Definition: constraints.hh:1016
Definition: linearproblem.hh:25
bool keepMatrix() const
Return whether the jacobian matrix is kept across calls to apply().
Definition: linearproblem.hh:154
StationaryLinearProblemSolver(const GO &go, LS &ls, V &x, Real reduction, Real min_defect=1e-99, int verbose=1)
Definition: linearproblem.hh:55