Line data Source code
1 : /*
2 : * SPDX-FileCopyrightText: Copyright (c) 2014: G-CSC, Goethe University Frankfurt
3 : * SPDX-License-Identifier: LicenseRef-UG4-LGPL-3.0
4 : *
5 : * Author: Arne Naegel, Andreas Kreienbuehl
6 : *
7 : * This file is part of UG4.
8 : *
9 : * UG4 is free software: you can redistribute it and/or modify it under the
10 : * terms of the GNU Lesser General Public License version 3 (as published by the
11 : * Free Software Foundation) with the following additional attribution
12 : * requirements (according to LGPL/GPL v3 §7):
13 : *
14 : * (1) The following notice must be displayed in the Appropriate Legal Notices
15 : * of covered and combined works: "Based on UG4 (www.ug4.org/license)".
16 : *
17 : * (2) The following notice must be displayed at a prominent place in the
18 : * terminal output of covered works: "Based on UG4 (www.ug4.org/license)".
19 : *
20 : * (3) The following bibliography is recommended for citation and must be
21 : * preserved in all covered files:
22 : * "Reiter, S., Vogel, A., Heppner, I., Rupp, M., and Wittum, G. A massively
23 : * parallel geometric multigrid solver on hierarchically distributed grids.
24 : * Computing and visualization in science 16, 4 (2013), 151-164"
25 : * "Vogel, A., Reiter, S., Rupp, M., Nägel, A., and Wittum, G. UG4 -- a novel
26 : * flexible software system for simulating pde based models on high performance
27 : * computers. Computing and visualization in science 16, 4 (2013), 165-179"
28 : *
29 : * This program is distributed in the hope that it will be useful,
30 : * but WITHOUT ANY WARRANTY; without even the implied warranty of
31 : * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
32 : * GNU Lesser General Public License for more details.
33 : */
34 :
35 : #ifndef __H__LIMEX__TIME_INTEGRATOR_HPP__
36 : #define __H__LIMEX__TIME_INTEGRATOR_HPP__
37 :
38 : #if __cplusplus >= 201103L
39 : #define OVERRIDE override
40 : #else
41 : #define OVERRIDE
42 : #endif
43 :
44 : // std headers.
45 : #include <string>
46 :
47 : // UG4 headers
48 : #include "lib_algebra/operator/interface/operator.h"
49 : #include "lib_algebra/operator/interface/operator_inverse.h"
50 : #include "lib_algebra/operator/linear_solver/linear_solver.h"
51 : #include "lib_disc/function_spaces/grid_function.h"
52 : #include "lib_disc/assemble_interface.h" // TODO: missing IAssemble in following file:
53 : #include "lib_disc/operator/linear_operator/assembled_linear_operator.h"
54 : #include "lib_disc/operator/non_linear_operator/assembled_non_linear_operator.h"
55 : #include "lib_disc/spatial_disc/domain_disc.h"
56 : #include "lib_disc/time_disc/time_disc_interface.h"
57 : #include "lib_disc/time_disc/theta_time_step.h"
58 : #include "lib_disc/time_disc/solution_time_series.h"
59 : #include "lib_disc/time_disc/time_integrator_subject.hpp"
60 : #include "lib_disc/time_disc/time_integrator_observers/time_integrator_observer_interface.h"
61 : #include "lib_disc/function_spaces/grid_function_util.h" // SaveVectorForConnectionViewer
62 : #include "lib_disc/function_spaces/interpolate.h" //Interpolate
63 : #include "lib_disc/function_spaces/integrate.h" //Integral
64 : #include "lib_disc/function_spaces/grid_function.h" //GridFunction
65 :
66 :
67 : // Plugin headers
68 : #include "time_extrapolation.h"
69 : #include "../limex_tools.h"
70 :
71 : namespace ug {
72 :
73 : /// Integrates over a given time interval [a,b] with step size dt
74 : template<class TDomain, class TAlgebra>
75 : class ITimeIntegrator
76 : : public IOperator< GridFunction<TDomain, TAlgebra> >,
77 : public TimeIntegratorSubject<TDomain, TAlgebra>
78 : {
79 : public:
80 :
81 : typedef TAlgebra algebra_type;
82 : typedef typename TAlgebra::vector_type vector_type;
83 : typedef typename TAlgebra::matrix_type matrix_type;
84 :
85 :
86 : typedef TDomain domain_type;
87 : typedef GridFunction<TDomain, TAlgebra> grid_function_type;
88 :
89 : protected:
90 : double m_dt;
91 : double m_lower_tim;
92 : double m_upper_tim;
93 :
94 : double m_precisionBound;
95 :
96 : bool m_bNoLogOut;
97 :
98 :
99 : public:
100 : // constructor
101 0 : ITimeIntegrator()
102 0 : : m_dt(1.0), m_lower_tim(0.0), m_upper_tim(0.0), m_precisionBound(1e-10), m_bNoLogOut(false)
103 : {}
104 :
105 : /// virtual destructor
106 0 : virtual ~ITimeIntegrator() {};
107 :
108 : /// init operator depending on a function u
109 : /**
110 : * This method initializes the operator. Once initialized the 'apply'-method
111 : * can be called. The function u is passed here, since the linear operator
112 : * may be the linearization of some non-linear operator. Thus, the operator
113 : * depends on the linearization point.
114 : * If the operator is not a linearization, this method can be implemented
115 : * by simply calling init() and forgetting about the linearization point.
116 : *
117 : * \param[in] u function (linearization point)
118 : * \returns bool success flag
119 : */
120 0 : virtual void init(grid_function_type const& u)
121 : {
122 : // UG_ASSERT(m_spDomainDisc.valid(), "TimeIntegrator<TDomain, TAlgebra>::init: m_spDomainDisc invalid.");
123 : // m_spTimeDisc = make_sp(new time_disc_type(m_spDomainDisc, m_theta));
124 0 : }
125 :
126 :
127 : /// init operator
128 : /**
129 : * This method initializes the operator. Once initialized the 'apply'-method
130 : * can be called.
131 : * \returns bool success flag
132 : */
133 0 : void init()
134 0 : { UG_THROW("Please init with grid function!"); }
135 :
136 : /// prepares functions for application
137 :
138 : /* This method is used to prepare the in- and output vector used later in apply.
139 : * It can be called, after the init method has been called at least once.
140 : * The prepare method is e.g. used to set dirichlet values.
141 : */
142 0 : void prepare(grid_function_type& u) {}
143 :
144 : //! Apply operator
145 : /*! This method applies the operator, i.e, advances the time step*/
146 0 : void apply(grid_function_type& u1, const grid_function_type& u0)
147 0 : { UG_THROW("Fix interfaces!"); }
148 :
149 : virtual bool apply(SmartPtr<grid_function_type> u1, number t1, ConstSmartPtr<grid_function_type> u0, number t0) = 0;
150 :
151 : //! Set initial time step
152 0 : void set_time_step(double dt)
153 0 : { m_dt = dt; }
154 :
155 : double get_time_step()
156 0 : { return m_dt; }
157 :
158 0 : void set_precision_bound(double precisionBound)
159 0 : { m_precisionBound = precisionBound; return; }
160 :
161 0 : void set_no_log_out(bool bNoLogOut)
162 0 : { m_bNoLogOut = bNoLogOut; return; }
163 :
164 : };
165 :
166 :
167 : /// Integration of linear systems
168 : template<class TDomain, class TAlgebra>
169 : class ILinearTimeIntegrator : public ITimeIntegrator<TDomain, TAlgebra>
170 : {
171 :
172 : public:
173 : typedef ITimeIntegrator<TDomain, TAlgebra> base_type;
174 : typedef typename base_type::vector_type vector_type;
175 : typedef ILinearOperatorInverse<vector_type> linear_solver_type;
176 : typedef AssembledLinearOperator<TAlgebra> assembled_operator_type;
177 :
178 : // forward constructor
179 0 : ILinearTimeIntegrator()
180 0 : : base_type() {}
181 :
182 0 : ILinearTimeIntegrator(SmartPtr<linear_solver_type> lSolver)
183 0 : : base_type(), m_spLinearSolver(lSolver)
184 0 : {}
185 :
186 :
187 0 : void set_linear_solver(SmartPtr<linear_solver_type> lSolver)
188 0 : { m_spLinearSolver=lSolver;}
189 :
190 : protected:
191 : SmartPtr<linear_solver_type> m_spLinearSolver;
192 :
193 : };
194 :
195 :
196 : /// ITimeDiscDependentObject
197 : template<class TAlgebra>
198 0 : class ITimeDiscDependentObject
199 : {
200 : public:
201 : typedef ITimeDiscretization<TAlgebra> time_disc_type;
202 :
203 : ITimeDiscDependentObject(SmartPtr<time_disc_type> spTimeDisc) :
204 : m_spTimeDisc(spTimeDisc)
205 : {}
206 :
207 0 : SmartPtr<time_disc_type> get_time_disc() {return m_spTimeDisc;}
208 : protected:
209 : SmartPtr<time_disc_type> m_spTimeDisc;
210 : };
211 :
212 :
213 : /// Integrate over a given time interval (for a linear problem)
214 : template<class TDomain, class TAlgebra>
215 : class LinearTimeIntegrator :
216 : public ILinearTimeIntegrator<TDomain, TAlgebra>,
217 : public ITimeDiscDependentObject<TAlgebra>
218 :
219 : {
220 : private:
221 :
222 : protected:
223 : typedef ITimeDiscDependentObject<TAlgebra> tdisc_dep_type;
224 : public:
225 : typedef ILinearTimeIntegrator<TDomain, TAlgebra> base_type;
226 : typedef ITimeDiscretization<TAlgebra> time_disc_type;
227 : typedef typename base_type::grid_function_type grid_function_type;
228 : typedef VectorTimeSeries<typename base_type::vector_type> vector_time_series_type;
229 :
230 : // constructor
231 0 : LinearTimeIntegrator (SmartPtr< time_disc_type> tDisc)
232 0 : : base_type(), ITimeDiscDependentObject<TAlgebra>(tDisc) {}
233 :
234 : LinearTimeIntegrator (SmartPtr< time_disc_type> tDisc, SmartPtr<typename base_type::linear_solver_type> lSolver)
235 : : base_type(lSolver), ITimeDiscDependentObject<TAlgebra>(tDisc) {}
236 :
237 : bool apply(SmartPtr<grid_function_type> u1, number t1, ConstSmartPtr<grid_function_type> u0, number t0);
238 : };
239 :
240 :
241 : /// Integrate over a given time interval (for a linear problem)
242 : template<class TDomain, class TAlgebra>
243 : class ConstStepLinearTimeIntegrator :
244 : public ILinearTimeIntegrator<TDomain, TAlgebra>,
245 : public ITimeDiscDependentObject<TAlgebra>
246 : {
247 : protected:
248 : typedef ITimeDiscDependentObject<TAlgebra> tdisc_dep_type;
249 : public:
250 : typedef ILinearTimeIntegrator<TDomain, TAlgebra> base_type;
251 : typedef ITimeDiscretization<TAlgebra> time_disc_type;
252 : typedef typename base_type::linear_solver_type linear_solver_type;
253 : typedef typename base_type::grid_function_type grid_function_type;
254 : typedef VectorTimeSeries<typename base_type::vector_type> vector_time_series_type;
255 :
256 : private:
257 :
258 : protected:
259 :
260 : int m_numSteps;
261 : public:
262 :
263 : // constructor
264 0 : ConstStepLinearTimeIntegrator (SmartPtr<time_disc_type> tDisc)
265 0 : : base_type(), ITimeDiscDependentObject<TAlgebra>(tDisc), m_numSteps(1) {}
266 :
267 0 : ConstStepLinearTimeIntegrator (SmartPtr<time_disc_type> tDisc, SmartPtr<typename base_type::linear_solver_type> lSolver)
268 0 : : base_type(lSolver), ITimeDiscDependentObject<TAlgebra>(tDisc), m_numSteps(1) {}
269 :
270 0 : void set_num_steps(int steps) {m_numSteps = steps;}
271 :
272 : bool apply(SmartPtr<grid_function_type> u1, number t1, ConstSmartPtr<grid_function_type> u0, number t0);
273 : };
274 :
275 :
276 : /// Integrate over a given time interval (for a linear problem)
277 : template<class TDomain, class TAlgebra>
278 : class TimeIntegratorLinearAdaptive :
279 : public ILinearTimeIntegrator<TDomain, TAlgebra>,
280 : public ITimeDiscDependentObject<TAlgebra>
281 : {
282 : protected:
283 : typedef ITimeDiscDependentObject<TAlgebra> tdisc_dep_type;
284 :
285 : public:
286 : typedef ILinearTimeIntegrator<TDomain, TAlgebra> base_type;
287 : typedef ITimeDiscretization<TAlgebra> time_disc_type;
288 : typedef typename base_type::grid_function_type grid_function_type;
289 : typedef VectorTimeSeries<typename base_type::vector_type> vector_time_series_type;
290 :
291 : protected:
292 :
293 : SmartPtr<time_disc_type> m_spTimeDisc2; // set during init
294 :
295 : double m_tol, m_dtmin, m_dtmax;
296 :
297 :
298 : public:
299 0 : TimeIntegratorLinearAdaptive (SmartPtr<time_disc_type> tDisc1, SmartPtr<time_disc_type> tDisc2)
300 0 : : base_type(), ITimeDiscDependentObject<TAlgebra>(tDisc1), m_tol(1e-2), m_dtmin(-1.0), m_dtmax(-1.0)
301 : {
302 0 : m_spTimeDisc2 = tDisc2;
303 0 : }
304 :
305 0 : void init(grid_function_type const& u)
306 : {
307 : // call base
308 : base_type::init(u);
309 : //m_spTimeDisc2 = make_sp(new typename base_type::time_disc_type(base_type::m_spDomainDisc, base_type::m_theta));
310 0 : }
311 :
312 : bool apply(SmartPtr<grid_function_type> u1, number t1, ConstSmartPtr<grid_function_type> u0, number t0);
313 :
314 0 : void set_tol(double tol) {m_tol = tol;}
315 0 : void set_time_step_min(number dt) {m_dtmin = dt;}
316 0 : void set_time_step_max(number dt) {m_dtmax = dt;}
317 : };
318 :
319 :
320 : class TimeStepBounds
321 : {
322 : public:
323 : TimeStepBounds()
324 0 : : m_dtMin(0.0), m_dtMax(std::numeric_limits<double>::max()), m_redFac(1.0), m_incFac(1.0)
325 : {}
326 :
327 0 : void set_dt_min(double min) { m_dtMin = min; }
328 0 : double get_dt_min() { return m_dtMin; }
329 :
330 0 : void set_dt_max(double max) { m_dtMax = max; }
331 0 : double get_dt_max() { return m_dtMax; }
332 :
333 0 : void set_reduction_factor(double dec) { m_redFac = dec; }
334 0 : double get_reduction_factor() { return m_redFac; }
335 :
336 0 : void set_increase_factor(double inc) { m_incFac = inc; }
337 0 : double get_increase_factor() { return m_incFac; }
338 :
339 : void rescale(double alpha)
340 : { m_dtMin*= alpha; m_dtMax*= alpha;}
341 :
342 : protected:
343 : double m_dtMin, m_dtMax;
344 : double m_redFac, m_incFac;
345 : };
346 :
347 : /// Integration of non-linear systems (with bounds on dt)
348 : template<class TDomain, class TAlgebra>
349 0 : class INonlinearTimeIntegrator
350 : : public ITimeIntegrator<TDomain, TAlgebra>
351 : {
352 : public:
353 : typedef ITimeIntegrator<TDomain, TAlgebra> base_type;
354 : typedef typename base_type::vector_type vector_type;
355 : typedef IOperatorInverse<vector_type> solver_type;
356 : typedef AssembledOperator<TAlgebra> assembled_operator_type;
357 :
358 0 : INonlinearTimeIntegrator()
359 0 : : m_dtBounds() {}
360 :
361 0 : void set_solver(SmartPtr<solver_type> solver)
362 0 : { m_spSolver=solver;}
363 :
364 : ConstSmartPtr<solver_type> get_solver() const
365 : { return m_spSolver;}
366 :
367 : SmartPtr<solver_type> get_solver()
368 : { return m_spSolver;}
369 :
370 0 : void set_dt_min(double min) { m_dtBounds.set_dt_min(min); }
371 : double get_dt_min() { return m_dtBounds. get_dt_min(); }
372 :
373 0 : void set_dt_max(double max) { m_dtBounds.set_dt_max(max); }
374 : double get_dt_max() { return m_dtBounds.get_dt_max(); }
375 :
376 0 : void set_reduction_factor(double dec) { m_dtBounds.set_reduction_factor(dec); }
377 : double get_reduction_factor() { return m_dtBounds.get_reduction_factor(); }
378 :
379 0 : void set_increase_factor(double inc) { m_dtBounds.set_increase_factor(inc); }
380 : double get_increase_factor() { return m_dtBounds.get_increase_factor(); }
381 :
382 : protected:
383 : SmartPtr<solver_type> m_spSolver;
384 : TimeStepBounds m_dtBounds;
385 :
386 : };
387 :
388 :
389 : //! This class integrates (t0, t1] with stops at intermediate points tk.
390 : template<class TDomain, class TAlgebra>
391 : class DiscontinuityIntegrator :
392 : public INonlinearTimeIntegrator<TDomain, TAlgebra>
393 : {
394 : public:
395 :
396 : typedef INonlinearTimeIntegrator<TDomain, TAlgebra> base_type;
397 : typedef typename base_type::grid_function_type grid_function_type;
398 :
399 0 : DiscontinuityIntegrator(SmartPtr<base_type> baseIntegrator) :
400 0 : base_type(), m_wrappedIntegrator(baseIntegrator), m_timePoints() {};
401 :
402 0 : bool apply(SmartPtr<grid_function_type> u1, number t1, ConstSmartPtr<grid_function_type> u0, number t0)
403 : {
404 : int dstep = 0;
405 : auto tpoint = m_timePoints.begin();
406 0 : double tcurr = (tpoint == m_timePoints.end()) ? t1 : (*tpoint);
407 : double eps = 1e-8;
408 :
409 : // Perform first step.
410 : //this->notify_init_step(u0, dstep, t0, tcurr-t0);
411 0 : bool status = m_wrappedIntegrator->apply(u1, tcurr*(1.0-eps), u0, t0);
412 0 : this->notify_finalize_step(u1, dstep++, tcurr, tcurr-t0);
413 :
414 : // Repeat for all intermediate points.
415 0 : while (tpoint != m_timePoints.end())
416 : {
417 : tpoint++;
418 0 : double tnext = (tpoint == m_timePoints.end()) ? t1 : (*tpoint);
419 :
420 : // Perform step.
421 : //this->notify_init_step(u1, dstep, tcurr, tnext-tcurr);
422 0 : status = status && m_wrappedIntegrator->apply(u1, tnext*(1.0-eps), u1, tcurr);
423 0 : this->notify_finalize_step(u1, dstep++, tnext, tnext-tcurr);
424 :
425 : tcurr = tnext;
426 : }
427 0 : this->notify_end(u1, dstep, t1, 0.0);
428 0 : return status;
429 : }
430 :
431 0 : void insert_points (std::vector<double> points)
432 : {
433 0 : m_timePoints = points;
434 0 : }
435 : protected:
436 : SmartPtr<base_type> m_wrappedIntegrator;
437 : std::vector<double> m_timePoints;
438 : };
439 :
440 : } // namespace ug
441 :
442 : #include "time_integrator_impl.hpp"
443 :
444 : #endif /* __H__LIMEX__TIME_INTEGRATOR_HPP__ */
|
