1#ifndef DUNE_PDELAB_STATIONARY_LINEARPROBLEM_HH
2#define DUNE_PDELAB_STATIONARY_LINEARPROBLEM_HH
9#include <dune/pdelab/backend/interface.hh>
10#include <dune/pdelab/constraints/common/constraints.hh>
11#include <dune/pdelab/backend/solver.hh>
24 template<
class RFType>
25 struct StationaryLinearProblemSolverResult : LinearSolverResult<RFType>
29 double assembler_time;
30 double linear_solver_time;
31 int linear_solver_iterations;
33 StationaryLinearProblemSolverResult()
37 , linear_solver_time(0.0)
38 , linear_solver_iterations(0)
43 template<
typename GO,
typename LS,
typename V>
44 class StationaryLinearProblemSolver
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;
49 using W = Dune::PDELab::Backend::Vector<TrialGridFunctionSpace,typename V::ElementType>;
50 typedef GO GridOperator;
53 typedef StationaryLinearProblemSolverResult<double> Result;
55 StationaryLinearProblemSolver(
const GO& go, LS& ls, V& x, Real reduction, Real min_defect = 1e-99,
int verbose=1)
59 , _reduction(reduction)
60 , _min_defect(min_defect)
61 , _hanging_node_modifications(false)
66 StationaryLinearProblemSolver (
const GO& go, LS& ls, Real reduction, Real min_defect = 1e-99,
int verbose=1)
70 , _reduction(reduction)
71 , _min_defect(min_defect)
72 , _hanging_node_modifications(false)
95 StationaryLinearProblemSolver(
const GO& go, LS& ls, V& x,
const ParameterTree& params)
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))
124 StationaryLinearProblemSolver(
const GO& go, LS& ls,
const ParameterTree& params)
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))
136 void setHangingNodeModifications(
bool b)
138 _hanging_node_modifications=b;
142 bool hangingNodeModifications()
const
144 return _hanging_node_modifications;
148 void setKeepMatrix(
bool b)
154 bool keepMatrix()
const
159 const Result& result()
const
164 void apply(V& x,
bool reuse_matrix =
false) {
169 void apply (
bool reuse_matrix =
false)
172 double timing,assembler_time=0;
179 _jacobian = std::make_shared<M>(_go);
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);
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);
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;
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 && _verbose>=1)
234 std::cout <<
"=== solving (reduction: " << red <<
") ";
236 std::cout << std::flush;
238 std::cout << std::endl;
240 _ls.apply(*_jacobian,z,r,red);
241 _linear_solver_result = _ls.result();
244 if (_go.trialGridFunctionSpace().gridView().comm().rank()==0 && _verbose>=1)
245 std::cout << timing <<
" s" << std::endl;
246 _res.linear_solver_time = timing;
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;
258 if (_hanging_node_modifications)
261 if (_hanging_node_modifications)
262 _go.localAssembler().backtransform(*_x);
275 const Dune::PDELab::LinearSolverResult<double>& ls_result()
const{
276 return _linear_solver_result;
279 Real reduction()
const
284 void setReduction(Real reduction)
286 _reduction = reduction;
294 std::shared_ptr<M> _jacobian;
297 Dune::PDELab::LinearSolverResult<double> _linear_solver_result;
299 bool _hanging_node_modifications;
A simple stop watch.
Definition: timer.hh:41
void reset() noexcept
Reset timer while keeping the running/stopped state.
Definition: timer.hh:55
double elapsed() const noexcept
Get elapsed user-time from last reset until now/last stop in seconds.
Definition: timer.hh:75
decltype(auto) apply(F &&f, ArgTuple &&args)
Apply function with arguments given as tuple.
Definition: apply.hh:46
void set_shifted_dofs(const CG &cg, typename XG::ElementType x, XG &xg)
Definition: constraints.hh:1014
auto max(ADLTag< 0 >, const V &v1, const V &v2)
implements binary Simd::max()
Definition: defaults.hh:79
Dune namespace.
Definition: alignedallocator.hh:14
A hierarchical structure of string parameters.