4 #ifndef DUNE_PDELAB_STATIONARY_LINEARPROBLEM_HH
5 #define DUNE_PDELAB_STATIONARY_LINEARPROBLEM_HH
9 #include <dune/common/timer.hh>
10 #include <dune/common/parametertree.hh>
27 template<
class RFType>
59 template<
typename GO,
typename LS,
typename V>
62 typedef typename Dune::template FieldTraits<typename V::ElementType >::real_type Real;
63 typedef typename GO::Traits::Jacobian M;
64 typedef typename GO::Traits::TrialGridFunctionSpace TrialGridFunctionSpace;
66 typedef GO GridOperator;
76 , _min_defect(min_defect)
77 , _hanging_node_modifications(false)
87 , _min_defect(min_defect)
88 , _hanging_node_modifications(false)
115 , _reduction(params.get<Real>(
"reduction"))
116 , _min_defect(params.get<Real>(
"min_defect",1
e-99))
117 , _hanging_node_modifications(params.get<bool>(
"hanging_node_modifications",false))
118 , _keep_matrix(params.get<bool>(
"keep_matrix",true))
119 , _verbose(params.get<int>(
"verbosity",1))
144 , _reduction(params.get<Real>(
"reduction"))
145 , _min_defect(params.get<Real>(
"min_defect",1
e-99))
146 , _hanging_node_modifications(params.get<bool>(
"hanging_node_modifications",false))
147 , _keep_matrix(params.get<bool>(
"keep_matrix",true))
148 , _verbose(params.get<int>(
"verbosity",1))
154 _hanging_node_modifications=b;
160 return _hanging_node_modifications;
182 void apply(V& x,
bool reuse_matrix =
false) {
188 void apply (
bool reuse_matrix =
false)
191 double timing,assembler_time=0;
196 if constexpr (linearSolverIsMatrixFree<LS>()){
197 if (_go.trialGridFunctionSpace().gridView().comm().rank()==0 && _verbose>=1)
198 std::cout <<
"=== matrix setup not required for matrix free solvers" << std::endl;
203 _jacobian = std::make_shared<M>(_go);
204 timing = watch.elapsed();
205 if (_go.trialGridFunctionSpace().gridView().comm().rank()==0 && _verbose>=1)
206 std::cout <<
"=== matrix setup (max) " << timing <<
" s" << std::endl;
208 assembler_time += timing;
210 else if (_go.trialGridFunctionSpace().gridView().comm().rank()==0 && _verbose>=1)
211 std::cout <<
"=== matrix setup skipped (matrix already allocated)" << std::endl;
214 if (_hanging_node_modifications)
217 _go.localAssembler().backtransform(*_x);
221 if constexpr (!linearSolverIsMatrixFree<LS>()){
224 (*_jacobian) = Real(0.0);
225 _go.jacobian(*_x,*_jacobian);
228 timing = watch.elapsed();
230 if constexpr (!linearSolverIsMatrixFree<LS>()){
231 if (_go.trialGridFunctionSpace().gridView().comm().rank()==0 && _verbose>=1)
234 std::cout <<
"=== matrix assembly SKIPPED" << std::endl;
236 std::cout <<
"=== matrix assembly (max) " << timing <<
" s" << std::endl;
240 assembler_time += timing;
245 W r(_go.testGridFunctionSpace(),0.0);
248 timing = watch.elapsed();
250 if (_go.trialGridFunctionSpace().gridView().comm().rank()==0 && _verbose>=1)
251 std::cout <<
"=== residual assembly (max) " << timing <<
" s" << std::endl;
252 assembler_time += timing;
255 auto defect = _ls.norm(r);
259 V z(_go.trialGridFunctionSpace(),0.0);
260 auto red = std::max(_reduction,_min_defect/defect);
261 if (_go.trialGridFunctionSpace().gridView().comm().rank()==0 && _verbose>=1)
263 std::cout <<
"=== solving (reduction: " << red <<
") ";
265 std::cout << std::flush;
267 std::cout << std::endl;
269 if constexpr (linearSolverIsMatrixFree<LS>()){
270 _ls.apply(z, r, red);
272 if constexpr (!linearSolverIsMatrixFree<LS>()){
273 _ls.apply(*_jacobian,z,r,red);
275 _linear_solver_result = _ls.result();
276 timing = watch.elapsed();
278 if (_go.trialGridFunctionSpace().gridView().comm().rank()==0 && _verbose>=1)
279 std::cout << timing <<
" s" << std::endl;
288 _res.
defect =
static_cast<double>(defect)*_linear_solver_result.
reduction;
292 if (_hanging_node_modifications)
295 if (_hanging_node_modifications)
296 _go.localAssembler().backtransform(*_x);
298 if constexpr (!linearSolverIsMatrixFree<LS>()){
312 return _linear_solver_result;
330 std::shared_ptr<M> _jacobian;
335 bool _hanging_node_modifications;
const Entity & e
Definition: localfunctionspace.hh:123
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
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
Definition: linearproblem.hh:29
int linear_solver_iterations
Definition: linearproblem.hh:34
double linear_solver_time
Definition: linearproblem.hh:33
StationaryLinearProblemSolverResult()
Definition: linearproblem.hh:36
RFType defect
Definition: linearproblem.hh:31
double assembler_time
Definition: linearproblem.hh:32
RFType first_defect
Definition: linearproblem.hh:30
Solve linear problems using a residual formulation.
Definition: linearproblem.hh:61
void setKeepMatrix(bool b)
Set whether the jacobian matrix should be kept across calls to apply().
Definition: linearproblem.hh:164
Real reduction() const
Definition: linearproblem.hh:315
void apply(bool reuse_matrix=false)
Solve linear problem (use initial guess that was passed at construction)
Definition: linearproblem.hh:188
StationaryLinearProblemSolver(const GO &go, LS &ls, Real reduction, Real min_defect=1e-99, int verbose=1)
Definition: linearproblem.hh:82
void setHangingNodeModifications(bool b)
Set whether the solver should apply the necessary transformations for calculations on hanging nodes.
Definition: linearproblem.hh:152
bool hangingNodeModifications() const
Return whether the solver performs the necessary transformations for calculations on hanging nodes.
Definition: linearproblem.hh:158
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:111
StationaryLinearProblemSolver(const GO &go, LS &ls, V &x, Real reduction, Real min_defect=1e-99, int verbose=1)
Definition: linearproblem.hh:71
StationaryLinearProblemSolverResult< double > Result
Definition: linearproblem.hh:69
const Dune::PDELab::LinearSolverResult< double > & ls_result() const
Definition: linearproblem.hh:311
void apply(V &x, bool reuse_matrix=false)
Solve linear problem with the provided initial guess.
Definition: linearproblem.hh:182
void discardMatrix()
Discard the stored Jacobian matrix.
Definition: linearproblem.hh:305
void setReduction(Real reduction)
Definition: linearproblem.hh:320
const Result & result() const
Return result object.
Definition: linearproblem.hh:176
bool keepMatrix() const
Return whether the jacobian matrix is kept across calls to apply().
Definition: linearproblem.hh:170
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:140