1 #ifndef DUNE_PDELAB_STATIONARYLINEARPROBLEM_HH
2 #define DUNE_PDELAB_STATIONARYLINEARPROBLEM_HH
6 #include <dune/common/timer.hh>
7 #include <dune/common/deprecated.hh>
8 #include <dune/common/parametertree.hh>
25 template<
class RFType>
44 template<
typename GO,
typename LS,
typename V>
47 typedef typename V::ElementType Real;
48 typedef typename GO::Traits::Jacobian M;
49 typedef typename GO::Traits::TrialGridFunctionSpace TrialGridFunctionSpace;
51 typedef GO GridOperator;
60 , _reduction(reduction)
61 , _min_defect(min_defect)
62 , _hanging_node_modifications(false)
71 , _reduction(reduction)
72 , _min_defect(min_defect)
73 , _hanging_node_modifications(false)
82 , _reduction(reduction)
83 , _min_defect(min_defect)
84 , _hanging_node_modifications(false)
111 , _reduction(params.get<typename V::ElementType>(
"reduction"))
112 , _min_defect(params.get<typename V::ElementType>(
"min_defect",1
e-99))
113 , _hanging_node_modifications(params.get<bool>(
"hanging_node_modifications",false))
114 , _keep_matrix(params.get<bool>(
"keep_matrix",true))
115 , _verbose(params.get<int>(
"verbosity",1))
140 , _reduction(params.get<typename V::ElementType>(
"reduction"))
141 , _min_defect(params.get<typename V::ElementType>(
"min_defect",1
e-99))
142 , _hanging_node_modifications(params.get<bool>(
"hanging_node_modifications",false))
143 , _keep_matrix(params.get<bool>(
"keep_matrix",true))
144 , _verbose(params.get<int>(
"verbosity",1))
150 _hanging_node_modifications=b;
156 return _hanging_node_modifications;
184 double timing,assembler_time=0;
191 _jacobian = make_shared<M>(_go);
192 timing = watch.elapsed();
193 if (_go.trialGridFunctionSpace().gridView().comm().rank()==0 && _verbose>=1)
194 std::cout <<
"=== matrix setup (max) " << timing <<
" s" << std::endl;
196 assembler_time += timing;
198 else if (_go.trialGridFunctionSpace().gridView().comm().rank()==0 && _verbose>=1)
199 std::cout <<
"=== matrix setup skipped (matrix already allocated)" << std::endl;
201 (*_jacobian) = Real(0.0);
202 if (_hanging_node_modifications)
205 _go.localAssembler().backtransform(*_x);
207 _go.jacobian(*_x,*_jacobian);
209 timing = watch.elapsed();
211 if (_go.trialGridFunctionSpace().gridView().comm().rank()==0 && _verbose>=1)
212 std::cout <<
"=== matrix assembly (max) " << timing <<
" s" << std::endl;
213 assembler_time += timing;
218 W r(_go.testGridFunctionSpace(),0.0);
221 timing = watch.elapsed();
223 if (_go.trialGridFunctionSpace().gridView().comm().rank()==0 && _verbose>=1)
224 std::cout <<
"=== residual assembly (max) " << timing <<
" s" << std::endl;
225 assembler_time += timing;
228 typename V::ElementType defect = _ls.norm(r);
232 V z(_go.trialGridFunctionSpace(),0.0);
233 typename V::ElementType red = std::max(_reduction,_min_defect/defect);
234 if (_go.trialGridFunctionSpace().gridView().comm().rank()==0)
235 std::cout <<
"=== solving (reduction: " << red <<
") ";
236 _ls.apply(*_jacobian,z,r,red);
237 _linear_solver_result = _ls.result();
238 timing = watch.elapsed();
240 if (_go.trialGridFunctionSpace().gridView().comm().rank()==0 && _verbose>=1)
241 std::cout << timing <<
" s" << std::endl;
250 _res.
defect =
static_cast<double>(defect)*_linear_solver_result.
reduction;
254 if (_hanging_node_modifications)
257 if (_hanging_node_modifications)
258 _go.localAssembler().backtransform(*_x);
272 return _linear_solver_result;
279 shared_ptr<M> _jacobian;
280 typename V::ElementType _reduction;
281 typename V::ElementType _min_defect;
284 bool _hanging_node_modifications;
StationaryLinearProblemSolver(const GO &go, V &x, LS &ls, typename V::ElementType reduction, typename V::ElementType min_defect=1e-99, int verbose=1)
Definition: linearproblem.hh:56
bool converged
Definition: solver.hh:32
void discardMatrix()
Discard the stored Jacobian matrix.
Definition: linearproblem.hh:265
void set_shifted_dofs(const CG &cg, typename XG::ElementType x, XG &xg)
Definition: common/constraints.hh:1009
double elapsed
Definition: solver.hh:34
void setHangingNodeModifications(bool b)
Set whether the solver should apply the necessary transformations for calculations on hanging nodes...
Definition: linearproblem.hh:148
StationaryLinearProblemSolverResult< double > Result
Definition: linearproblem.hh:54
void setKeepMatrix(bool b)
Set whether the jacobian matrix should be kept across calls to apply().
Definition: linearproblem.hh:160
bool keepMatrix() const
Return whether the jacobian matrix is kept across calls to apply().
Definition: linearproblem.hh:166
const Dune::PDELab::LinearSolverResult< double > & ls_result() const
Definition: linearproblem.hh:271
StationaryLinearProblemSolver(const GO &go, LS &ls, V &x, typename V::ElementType reduction, typename V::ElementType min_defect=1e-99, int verbose=1)
Definition: linearproblem.hh:67
bool hangingNodeModifications() const
Return whether the solver performs the necessary transformations for calculations on hanging nodes...
Definition: linearproblem.hh:154
RFType reduction
Definition: solver.hh:35
void apply(V &x)
Definition: linearproblem.hh:176
BackendVectorSelectorHelper< Backend, GridFunctionSpace, FieldType >::Type Type
Definition: backendselector.hh:14
double linear_solver_time
Definition: linearproblem.hh:31
const Result & result() const
Definition: linearproblem.hh:171
Definition: linearproblem.hh:45
StationaryLinearProblemSolverResult()
Definition: linearproblem.hh:34
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:107
StationaryLinearProblemSolver(const GO &go, LS &ls, typename V::ElementType reduction, typename V::ElementType min_defect=1e-99, int verbose=1)
Definition: linearproblem.hh:78
const E & e
Definition: interpolate.hh:172
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:136
RFType conv_rate
Definition: solver.hh:36
RFType defect
Definition: linearproblem.hh:29
void apply()
Definition: linearproblem.hh:181
RFType first_defect
Definition: linearproblem.hh:28
int linear_solver_iterations
Definition: linearproblem.hh:32
double assembler_time
Definition: linearproblem.hh:30
Definition: linearproblem.hh:26
unsigned int iterations
Definition: solver.hh:33