Line data Source code
1 : /*
2 : * SPDX-FileCopyrightText: Copyright (c) 2014-2025: G-CSC, Goethe University Frankfurt
3 : * SPDX-License-Identifier: LicenseRef-UG4-LGPL-3.0
4 : *
5 : * Author: Arne Naegel
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 LIMEX_INTEGRATOR_HPP_
36 : #define LIMEX_INTEGRATOR_HPP_
37 : /*
38 : #define XMTHREAD_BOOST
39 : #ifdef XMTHREAD_BOOST
40 : #include <boost/thread/thread.hpp>
41 : #endif
42 : */
43 : #include <string>
44 :
45 : #include "common/stopwatch.h"
46 :
47 : #include "lib_algebra/operator/interface/operator.h"
48 : #include "lib_algebra/operator/interface/operator_inverse.h"
49 : #include "lib_algebra/operator/linear_solver/linear_solver.h"
50 : #include "lib_algebra/operator/debug_writer.h"
51 :
52 : #include "lib_disc/function_spaces/grid_function.h"
53 : #include "lib_disc/assemble_interface.h" // TODO: missing IAssemble in following file:
54 : #include "lib_disc/operator/linear_operator/assembled_linear_operator.h"
55 : #include "lib_disc/operator/non_linear_operator/assembled_non_linear_operator.h"
56 : #include "lib_disc/spatial_disc/domain_disc.h"
57 : #include "lib_disc/time_disc/time_disc_interface.h"
58 : #include "lib_disc/time_disc/theta_time_step.h"
59 : #include "lib_disc/time_disc/solution_time_series.h"
60 : #include "lib_disc/function_spaces/grid_function_util.h" // SaveVectorForConnectionViewer
61 : #include "lib_disc/function_spaces/metric_spaces.h"
62 : #include "lib_disc/io/vtkoutput.h"
63 :
64 : #include "lib_grid/refinement/refiner_interface.h"
65 :
66 :
67 : // own headers
68 : #include "time_extrapolation.h"
69 : #include "time_integrator.hpp"
70 : #include "../limex_tools.h"
71 : //#include "../multi_thread_tools.h"
72 :
73 :
74 : #ifdef UG_JSON
75 : #include <nlohmann/json.hpp>
76 : #endif
77 :
78 : #undef LIMEX_MULTI_THREAD
79 :
80 : namespace ug {
81 :
82 :
83 : /*
84 : template<class TI>
85 : class TimeIntegratorThread : public TI
86 : {
87 : TimeIntegratorThread(SmartPtr<grid_function_type> u1, number t1, ConstSmartPtr<grid_function_type> u0, number t0)
88 : {
89 :
90 : }
91 :
92 : // thre
93 : void apply()
94 : {
95 : TI::apply(u1, t1, u0, t0);
96 : }
97 :
98 :
99 : };*/
100 :
101 0 : static void MyPrintError(UGError &err)
102 : {
103 0 : for(size_t i=0;i<err.num_msg();++i)
104 : {
105 : UG_LOG("MYERROR "<<i<<":"<<err.get_msg(i)<<std::endl);
106 0 : UG_LOG(" [at "<<err.get_file(i)<<
107 : ", line "<<err.get_line(i)<<"]\n");
108 : }
109 0 : }
110 :
111 : class ILimexRefiner
112 : {
113 : virtual ~ILimexRefiner(){};
114 :
115 : protected:
116 : SmartPtr<IRefiner> m_spRefiner;
117 : };
118 :
119 : //! Abstract class for the cost of a limex stage.
120 : /*! The LimexTimeIntegrator requires a model for computing the cost per stage. */
121 : class ILimexCostStrategy
122 : {
123 : public:
124 : /// destructor
125 : virtual ~ILimexCostStrategy(){}
126 :
127 : /// provides the cost for all 0...nstages stages.
128 : virtual void update_cost(std::vector<number> &costA, const std::vector<size_t> &vSteps, const size_t nstages) = 0;
129 : };
130 :
131 :
132 : /// Cost is identical to (summation over) number of steps
133 : class LimexDefaultCost : public ILimexCostStrategy
134 : {
135 : public:
136 0 : LimexDefaultCost(){};
137 0 : void update_cost(std::vector<number> &m_costA, const std::vector<size_t> &m_vSteps, const size_t nstages)
138 : {
139 : UG_ASSERT(m_costA.size() >= nstages, "Huhh: Vectors match in size:" << m_costA.size() << "vs." << nstages);
140 :
141 0 : UG_LOG("A_0="<< m_vSteps[0] << std::endl);
142 0 : m_costA[0] = (1.0)*m_vSteps[0];
143 0 : for (size_t i=1; i<nstages; ++i)
144 : {
145 0 : m_costA[i] = m_costA[i-1] + (1.0)*m_vSteps[i];
146 0 : UG_LOG("A_i="<< m_vSteps[i] << std::endl);
147 : }
148 0 : }
149 : };
150 :
151 : /// For
152 : class LimexNonlinearCost : public ILimexCostStrategy
153 : {
154 : public:
155 0 : LimexNonlinearCost() :
156 0 : m_cAssemble (1.0), m_cMatAdd(0.0), m_cSolution(1.0), m_useGamma(1)
157 : {}
158 :
159 0 : void update_cost(std::vector<number> &m_costA, const std::vector<size_t> &nSteps, const size_t nstages)
160 : {
161 : UG_ASSERT(m_costA.size() >= nstages, "Huhh: Vectors match in size:" << m_costA.size() << "vs." << nstages);
162 :
163 : // 3 assemblies (M0, J0, Gamma0)
164 0 : m_costA[0] = (2.0+m_useGamma)*m_cAssemble + nSteps[0]*((1.0+m_useGamma)*m_cMatAdd + m_cSolution);
165 0 : for (size_t i=1; i<=nstages; ++i)
166 : {
167 : // n-1 assemblies for Mk, 2n MatAdds, n solutions
168 0 : m_costA[i] = m_costA[i-1] + (nSteps[i]-1) * m_cAssemble + nSteps[i]*((1.0+m_useGamma)*m_cMatAdd + m_cSolution);
169 : }
170 0 : }
171 :
172 : protected:
173 : double m_cAssemble; ///! Cost for matrix assembly
174 : double m_cMatAdd; ///! Cost for J=A+B
175 : double m_cSolution; ///! Cost for solving Jx=b
176 :
177 : int m_useGamma;
178 :
179 : };
180 :
181 :
182 : class LimexTimeIntegratorConfig
183 :
184 : {
185 : public:
186 : LimexTimeIntegratorConfig()
187 : : m_nstages(2),
188 : m_tol(0.01),
189 : m_rhoSafety(0.8),
190 : m_sigmaReduction(0.5),
191 : m_greedyOrderIncrease(2.0),
192 : m_max_reductions(2),
193 : m_asymptotic_order(1000),
194 : m_useCachedMatrices(false),
195 : m_conservative(0)
196 : {}
197 :
198 : LimexTimeIntegratorConfig(unsigned int nstages)
199 0 : : m_nstages(nstages),
200 0 : m_tol(0.01),
201 0 : m_rhoSafety(0.8),
202 0 : m_sigmaReduction(0.5),
203 0 : m_greedyOrderIncrease(2.0),
204 0 : m_max_reductions(2),
205 0 : m_asymptotic_order(1000),
206 0 : m_useCachedMatrices(false),
207 0 : m_conservative(0)
208 : {}
209 :
210 :
211 : protected:
212 : unsigned int m_nstages; ///< Number of Aitken-Neville stages
213 : double m_tol;
214 : double m_rhoSafety;
215 : double m_sigmaReduction;
216 :
217 : double m_greedyOrderIncrease;
218 :
219 : size_t m_max_reductions;
220 : size_t m_asymptotic_order; ///< For PDEs, we may apply an symptotic order reduction
221 :
222 : bool m_useCachedMatrices;
223 : unsigned int m_conservative; // stepping back?
224 :
225 : #ifdef UG_JSON
226 : NLOHMANN_DEFINE_TYPE_INTRUSIVE(LimexTimeIntegratorConfig,
227 : m_nstages, m_tol, m_rhoSafety,
228 : m_sigmaReduction, m_greedyOrderIncrease, m_max_reductions, m_asymptotic_order,
229 : m_useCachedMatrices, m_conservative)
230 : #endif
231 : public:
232 : std::string config_string() const {
233 : #ifdef UG_JSON
234 : try{
235 : nlohmann::json j(*this);
236 : // to_json(j, *this);
237 : return j.dump();
238 : } catch (...) {
239 : return std::string("EXCEPTION!!!");
240 : }
241 : #else
242 0 : return std::string("LimexTimeIntegratorConfig::config_string()");
243 : #endif
244 : }
245 :
246 : };
247 :
248 : //! Base class for LIMEX time integrator
249 : template<class TDomain, class TAlgebra>
250 : class LimexTimeIntegrator
251 : : public INonlinearTimeIntegrator<TDomain, TAlgebra>,
252 : public DebugWritingObject<TAlgebra>,
253 : public LimexTimeIntegratorConfig // TODO: should become a 'has a'.
254 : {
255 :
256 :
257 : public:
258 : typedef TAlgebra algebra_type;
259 : typedef typename algebra_type::matrix_type matrix_type;
260 : typedef typename algebra_type::vector_type vector_type;
261 : typedef GridFunction<TDomain, TAlgebra> grid_function_type;
262 : typedef INonlinearTimeIntegrator<TDomain, TAlgebra> base_type;
263 : typedef typename base_type::solver_type solver_type;
264 :
265 : typedef IDomainDiscretization<algebra_type> domain_discretization_type;
266 : typedef LinearImplicitEuler<algebra_type> timestep_type;
267 : typedef AitkenNevilleTimex<vector_type> timex_type;
268 : typedef INonlinearTimeIntegrator<TDomain, TAlgebra> itime_integrator_type;
269 : typedef SimpleTimeIntegrator<TDomain, TAlgebra> time_integrator_type;
270 : typedef ISubDiagErrorEst<vector_type> error_estim_type;
271 :
272 : //! Contains all data for parallel execution of time steps
273 : class ThreadData
274 : {
275 : //typedef boost::thread thread_type;
276 : public:
277 :
278 0 : ThreadData(SmartPtr<timestep_type> spTimeStep)
279 : : m_stepper(spTimeStep)
280 : {}
281 :
282 : ThreadData(){}
283 : SmartPtr<timestep_type> get_time_stepper()
284 : { return m_stepper; }
285 :
286 :
287 : void set_solver(SmartPtr<solver_type> solver)
288 0 : { m_solver = solver;}
289 :
290 : SmartPtr<solver_type> get_solver()
291 : { return m_solver;}
292 :
293 : void set_error(int e)
294 : { m_error=e; }
295 :
296 :
297 : void get_error()
298 : { /*return m_error;*/ }
299 :
300 :
301 : void set_solution(SmartPtr<grid_function_type> sol)
302 0 : { m_sol = sol;}
303 :
304 : SmartPtr<grid_function_type> get_solution()
305 : { return m_sol;}
306 :
307 : void set_derivative(SmartPtr<grid_function_type> sol)
308 0 : { m_dot = sol;}
309 :
310 : SmartPtr<grid_function_type> get_derivative()
311 : { return m_dot;}
312 :
313 : protected:
314 : // includes time step series
315 : SmartPtr<timestep_type> m_stepper;
316 : SmartPtr<solver_type> m_solver;
317 :
318 : SmartPtr<grid_function_type> m_sol;
319 : SmartPtr<grid_function_type> m_dot;
320 : int m_error;
321 : };
322 :
323 : typedef std::vector<SmartPtr<ThreadData> > thread_vector_type;
324 :
325 : public:
326 :
327 : /// forward debug info to time integrators
328 0 : void set_debug_for_timestepper(SmartPtr<IDebugWriter<algebra_type> > spDebugWriter)
329 : {
330 0 : for(size_t i=0; i<m_vThreadData.size(); ++i)
331 : {
332 0 : m_vThreadData[i].get_time_stepper()->set_debug(spDebugWriter);
333 : }
334 : UG_LOG("set_debug:" << m_vThreadData.size());
335 0 : }
336 :
337 : //using VectorDebugWritingObject<vector_type>::set_debug;
338 : protected:
339 : //using VectorDebugWritingObject<vector_type>::write_debug;
340 :
341 :
342 :
343 : public:
344 0 : LimexTimeIntegrator(int nstages):
345 : LimexTimeIntegratorConfig(nstages),
346 0 : m_gamma(m_nstages+1),
347 0 : m_costA(m_nstages+1),
348 0 : m_monitor(((m_nstages)*(m_nstages))), // TODO: wasting memory here!
349 0 : m_workload(m_nstages),
350 0 : m_lambda(m_nstages),
351 0 : m_num_reductions(m_nstages, 0),
352 0 : m_consistency_error(m_nstages),
353 0 : m_spCostStrategy(make_sp<LimexDefaultCost>(new LimexDefaultCost())),
354 0 : m_spBanachSpace(new AlgebraicSpace<grid_function_type>()), // default algebraic space
355 0 : m_bInterrupt(false),
356 0 : m_limex_step(1)
357 : {
358 0 : m_vThreadData.reserve(m_nstages);
359 0 : m_vSteps.reserve(m_nstages);
360 : init_gamma(); // init exponents (i.e. k+1, k, 2k+1, ...)
361 0 : }
362 :
363 :
364 :
365 : /// tolerance
366 0 : void set_tolerance(double tol) { m_tol = tol;}
367 0 : void set_stepsize_safety_factor(double rho) { m_rhoSafety = rho;}
368 0 : void set_stepsize_reduction_factor(double sigma) { m_sigmaReduction = sigma;}
369 0 : void set_stepsize_greedy_order_factor(double sigma) { m_greedyOrderIncrease = sigma;}
370 :
371 0 : void set_max_reductions(size_t nred) { m_max_reductions = nred;}
372 0 : void set_asymptotic_order(size_t q) { m_asymptotic_order = q;}
373 :
374 0 : void set_start_step(size_t step){m_limex_step=step;}
375 :
376 : /// add an error estimator
377 0 : void add_error_estimator(SmartPtr<error_estim_type> spErrorEstim)
378 0 : { m_spErrorEstimator = spErrorEstim; }
379 :
380 : //! add a new stage (at end of list)
381 0 : void add_stage_base(size_t nsteps, SmartPtr<solver_type> solver, SmartPtr<domain_discretization_type> spDD, SmartPtr<domain_discretization_type> spGamma=SPNULL)
382 : {
383 : UG_ASSERT(m_vThreadData.size() == m_vSteps.size(), "ERROR: m_vThreadData and m_vSteps differ in size!");
384 :
385 : UG_ASSERT(m_vThreadData.empty() || m_vSteps.back()<nsteps, "ERROR: Sequence of steps must be increasing." );
386 :
387 : // all entries have been set
388 0 : if (m_vThreadData.size() == m_nstages)
389 0 : { return; }
390 :
391 :
392 : // a) set number of steps
393 0 : m_vSteps.push_back(nsteps);
394 :
395 : // b) set time-stepper
396 : SmartPtr<timestep_type> limexStepSingleton;
397 : #ifndef LIMEX_MULTI_THREAD
398 :
399 0 : if(m_vThreadData.size()>0) {
400 : // re-use time-stepper (if applicable)
401 0 : limexStepSingleton = m_vThreadData.back().get_time_stepper();
402 : }
403 : else
404 : {
405 : // create time-stepper
406 : #endif
407 : // for mult-threading, each info object has own time-stepper
408 0 : if (spGamma.invalid()) {
409 0 : limexStepSingleton = make_sp(new timestep_type(spDD));
410 : } else {
411 0 : limexStepSingleton = make_sp(new timestep_type(spDD, spDD, spGamma));
412 : }
413 : UG_ASSERT(limexStepSingleton.valid(), "Huhh: Invalid pointer")
414 : #ifndef LIMEX_MULTI_THREAD
415 : }
416 : #endif
417 : // propagate debug info
418 0 : limexStepSingleton->set_debug(VectorDebugWritingObject<vector_type>::vector_debug_writer());
419 0 : m_vThreadData.push_back(ThreadData(limexStepSingleton));
420 :
421 : // c) set solver
422 : UG_ASSERT(solver.valid(), "Huhh: Need to supply solver!");
423 0 : m_vThreadData.back().set_solver(solver);
424 : UG_ASSERT(m_vThreadData.back().get_solver().valid(), "Huhh: Need to supply solver!");
425 : }
426 :
427 0 : void add_stage(size_t nsteps, SmartPtr<solver_type> solver, SmartPtr<domain_discretization_type> spDD)
428 0 : { add_stage_base(nsteps, solver, spDD); }
429 :
430 0 : void add_stage_ext(size_t nsteps, SmartPtr<solver_type> solver, SmartPtr<domain_discretization_type> spDD, SmartPtr<domain_discretization_type> spGamma)
431 0 : { add_stage_base(nsteps, solver, spDD, spGamma); }
432 :
433 : ///! TODO: remove this function!
434 0 : void add_stage(size_t i, size_t nsteps, SmartPtr<domain_discretization_type> spDD, SmartPtr<solver_type> solver)
435 : {
436 : UG_LOG("WARNING: add_stage(i, nsteps ,...) is deprecated. Please use 'add_stage(nsteps ,...) instead!'");
437 0 : add_stage(nsteps, solver, spDD);
438 0 : }
439 :
440 : //!
441 : SmartPtr<timestep_type> get_time_stepper(size_t i)
442 : {
443 : UG_ASSERT(i<m_vThreadData.size(), "Huhh: Invalid entry");
444 : return m_vThreadData[i].get_time_stepper();
445 : }
446 :
447 :
448 :
449 :
450 : protected:
451 : //! Initialize integrator threads (w/ solutions)
452 : void init_integrator_threads(ConstSmartPtr<grid_function_type> u);
453 :
454 : //! (Tentatively) apply integrators
455 : int apply_integrator_threads(number dtcurr, ConstSmartPtr<grid_function_type> u0, number t0, size_t nstages);
456 :
457 : //! e.g. wait for all threads to complete
458 : void join_integrator_threads();
459 :
460 : //! Override thread-wise solutions with common solution
461 : void update_integrator_threads(ConstSmartPtr<grid_function_type> ucommon, number t);
462 :
463 : //! Dispose integrator threads (w/ solutions)
464 : void dispose_integrator_threads();
465 :
466 : public:
467 : //! Integrating from t0 -> t1
468 : bool apply(SmartPtr<grid_function_type> u, number t1, ConstSmartPtr<grid_function_type> u0, number t0);
469 :
470 : number get_cost(size_t i) { return m_costA[i]; }
471 : number get_gamma(size_t i) { return m_gamma[i]; }
472 : number get_workload(size_t i) { return m_workload[i]; }
473 :
474 :
475 : protected:
476 : number& monitor(size_t k, size_t q) {
477 : UG_ASSERT(k<m_nstages, "Huhh: k mismatch");
478 : UG_ASSERT(q<m_nstages, "Huhh: q mismatch");
479 0 : return m_monitor[k+(m_nstages)*q];
480 : }
481 :
482 : /// aux: compute exponents gamma_k (for roots)
483 : void init_gamma()
484 : {
485 0 : for (size_t k=0; k<=m_nstages; ++k)
486 : {
487 0 : m_gamma[k] = 2.0+k;
488 : }
489 : }
490 :
491 : /// Updating workloads A_i for computing T_ii
492 : // (depends on m_vSteps, which must have been initialized!)
493 : void update_cost()
494 : {
495 0 : m_spCostStrategy->update_cost(m_costA, m_vSteps, m_nstages);
496 : }
497 :
498 : /// convergence monitor
499 : // (depends on cost, which must have been initialized!)
500 0 : void update_monitor()
501 : {
502 : UG_ASSERT(m_costA.size()>= m_nstages, "Cost vector too small!")
503 0 : for (size_t k=0; k<m_nstages-1; ++k)
504 : {
505 0 : UG_LOG("k= "<< k<< ", A[k]=" << m_costA[k] << ", gamma[k]=" << m_gamma[k] << "\t");
506 0 : for (size_t q=0; q<m_nstages-1; ++q)
507 : {
508 : // Deuflhard: Order and stepsize, ... eq. (3.7)
509 0 : double gamma = (m_costA[k+1] - m_costA[0] + 1.0)/(m_costA[q+1] - m_costA[0] + 1.0);
510 0 : double alpha = pow(m_tol, gamma);
511 :
512 : // for fixed order q, the monitor indicates the performance penalty compared to a strategy using k stages only
513 : // cf. eq. (4.6)
514 0 : monitor(k,q) = pow(alpha/(m_tol * m_rhoSafety), 1.0/m_gamma[k]);
515 0 : UG_LOG(monitor(k,q) << "[" << pow(alpha/(m_tol), 1.0/m_gamma[k]) <<"," << gamma<< ","<< m_costA[k+1] << "," << m_costA[q+1]<< ","<< alpha << "]" << "\t");
516 : // UG_LOG( << "\t");
517 :
518 : }
519 : UG_LOG(std::endl);
520 : }
521 :
522 0 : }
523 :
524 : // Find row k=1, ..., ntest-1 minimizing (estimated) error eps[kmin]
525 : // Also: predict column q with minimal workload W_{k+1,k} = A_{k+1} * lambda_{k+1}
526 : size_t find_optimal_solution(const std::vector<number>& eps, size_t ntest, /*size_t &kf,*/ size_t &qpred);
527 :
528 : public:
529 : /// setter for time derivative info (optional for setting $\Gamma$)
530 0 : void set_time_derivative(SmartPtr<grid_function_type> udot) {m_spDtSol = udot;}
531 :
532 : /// getter for time derivative info (optional for setting $\Gamma$)
533 0 : SmartPtr<grid_function_type> get_time_derivative() {return m_spDtSol;}
534 :
535 : /// status for time derivative info (optional for setting $\Gamma$)
536 0 : bool has_time_derivative() {return m_spDtSol!=SPNULL;}
537 :
538 :
539 0 : void enable_matrix_cache() { m_useCachedMatrices = true; } ///< Select classic LIMEX
540 0 : void disable_matrix_cache() { m_useCachedMatrices = false; } ///< Select approximate Newton (default)
541 :
542 0 : void select_cost_strategy(SmartPtr<ILimexCostStrategy> cost)
543 0 : { m_spCostStrategy = cost;}
544 :
545 :
546 : /// set banach space (e.g. for computing consistency error)
547 0 : void set_space(SmartPtr<IGridFunctionSpace<grid_function_type> > spSpace)
548 : {
549 0 : m_spBanachSpace = spSpace;
550 0 : UG_LOG("set_space:" << m_spBanachSpace->config_string());
551 0 : }
552 :
553 : /// interrupt execution of apply() by external call via observer
554 0 : void interrupt() {m_bInterrupt = true;}
555 :
556 0 : void set_conservative(bool c)
557 0 : { m_conservative = (c) ? 1 : 0; }
558 :
559 0 : std::string config_string() const
560 0 : { return LimexTimeIntegratorConfig::config_string(); }
561 :
562 : protected:
563 :
564 : SmartPtr<error_estim_type> m_spErrorEstimator; // (smart ptr for) error estimator
565 :
566 : std::vector<size_t> m_vSteps; ///< generating sequence for extrapolation
567 : std::vector<ThreadData> m_vThreadData; ///< vector with thread information
568 :
569 : std::vector<number> m_gamma; ///< gamma_i: exponent
570 :
571 : std::vector<number> m_costA; ///< Cost A_i (for completing stage i)
572 : std::vector<number> m_monitor; ///< Convergence monitor \alpha
573 :
574 : std::vector<number> m_workload;
575 : std::vector<number> m_lambda;
576 :
577 : std::vector<size_t> m_num_reductions; ///< history of reductions
578 :
579 : std::vector<number> m_consistency_error; ///<Consistency error
580 :
581 : SmartPtr<grid_function_type> m_spDtSol; ///< Time derivative
582 :
583 : SmartPtr<ILimexCostStrategy> m_spCostStrategy;
584 :
585 : /// metric space
586 : SmartPtr<IGridFunctionSpace<grid_function_type> > m_spBanachSpace;
587 :
588 : bool m_bInterrupt;
589 : int m_limex_step; ///<Current counter
590 :
591 :
592 : };
593 :
594 :
595 : /*! Create private solutions for each thread */
596 : template<class TDomain, class TAlgebra>
597 0 : void LimexTimeIntegrator<TDomain,TAlgebra>::init_integrator_threads(ConstSmartPtr<grid_function_type> u)
598 : {
599 : PROFILE_FUNC_GROUP("limex");
600 0 : const int nstages = m_vThreadData.size()-1;
601 0 : for (int i=nstages; i>=0; --i)
602 : {
603 0 : m_vThreadData[i].set_solution(u->clone());
604 0 : m_vThreadData[i].set_derivative(u->clone());
605 0 : m_vThreadData[i].get_time_stepper()->set_matrix_cache(m_useCachedMatrices);
606 : }
607 0 : }
608 :
609 : /*! Create private solutions for each thread */
610 : template<class TDomain, class TAlgebra>
611 0 : void LimexTimeIntegrator<TDomain,TAlgebra>::dispose_integrator_threads()
612 : {
613 : PROFILE_FUNC_GROUP("limex");
614 0 : const int nstages = m_vThreadData.size()-1;
615 0 : for (int i=nstages; i>=0; --i)
616 : {
617 0 : m_vThreadData[i].set_solution(SPNULL);
618 0 : m_vThreadData[i].set_derivative(SPNULL);
619 : // m_vThreadData[i].get_time_stepper()->set_matrix_cache(m_useCachedMatrices);
620 : }
621 0 : }
622 :
623 :
624 : // create (& execute) threads
625 : /*boost::thread_group g;
626 : typename thread_vector_type::reverse_iterator rit=m_vThreadData.rbegin();
627 : for (rit++; rit!= m_vThreadData.rend(); ++rit)
628 : {
629 :
630 : boost::thread *t =new boost::thread(boost::bind(&ThreadSafeTimeIntegrator::apply, *rit));
631 : //g.add_thread(t);
632 :
633 : g.create_thread(boost::bind(&ThreadSafeTimeIntegrator::apply, *rit));
634 :
635 : }*/
636 :
637 :
638 : /* TODO: PARALLEL execution?*/
639 : template<class TDomain, class TAlgebra>
640 0 : int LimexTimeIntegrator<TDomain,TAlgebra>::apply_integrator_threads(number dtcurr, ConstSmartPtr<grid_function_type> u0, number t0, size_t nstages)
641 : {
642 : PROFILE_FUNC_GROUP("limex");
643 : update_cost(); // compute cost A_i (alternative: measure times?)
644 0 : update_monitor(); // convergence monitor
645 :
646 : /*
647 : int tn = omp_get_thread_num();
648 : int nt = omp_get_num_threads();
649 : omp_set_num_threads(nstages);
650 : */
651 : int error = 0;
652 : //const int nstages = m_vThreadData.size()-1;
653 : // #pragma omp for private(i) // shared (nstages, u1) schedule(static)
654 0 : for (int i=0; i<=(int)nstages; ++i)
655 : {
656 : /*
657 : std::cerr << "I am " << tn << " of " << nt << " ("<< i<< "/" << nstages<<")!" << std::endl;
658 : UGMultiThreadEnvironment mt_env;
659 : */
660 :
661 : // copy data to private structure (one-to-many)
662 : //m_vThreadData[i].set_solution(u1->clone());
663 :
664 : // switch to "child" comm
665 : // mt_env.begin();
666 :
667 : // integrate (t0, t0+dtcurr)
668 0 : time_integrator_type integrator(m_vThreadData[i].get_time_stepper());
669 0 : integrator.set_time_step(dtcurr/m_vSteps[i]);
670 : integrator.set_dt_min(dtcurr/m_vSteps[i]);
671 : integrator.set_dt_max(dtcurr/m_vSteps[i]);
672 : integrator.set_reduction_factor(0.0); // quit immediately, if step fails
673 0 : integrator.set_solver(m_vThreadData[i].get_solver());
674 0 : integrator.set_derivative(m_vThreadData[i].get_derivative());
675 :
676 0 : if(this->debug_writer_valid())
677 : {
678 0 : integrator.set_debug(this->debug_writer());
679 : char debug_name_ext[16]; snprintf(debug_name_ext, 16, "%04d", i);
680 0 : this->enter_debug_writer_section(std::string("Stage_") + debug_name_ext);
681 : }
682 :
683 : UG_ASSERT(m_spBanachSpace.valid(), "Huhh: Need valid (default) banach space");
684 0 : integrator.set_banach_space(m_spBanachSpace);
685 0 : UG_LOG("Set space:" << m_spBanachSpace->config_string());
686 :
687 : bool exec = true;
688 : try
689 : {
690 0 : exec = integrator.apply(m_vThreadData[i].get_solution(), t0+dtcurr, u0, t0);
691 0 : m_consistency_error[i] = integrator.get_consistency_error();
692 : }
693 0 : catch(ug::UGError& err)
694 : {
695 :
696 : exec = false;
697 0 : error += (1 << i);
698 0 : UG_LOGN("Step "<< i<< " failed on stage " << i << ": " << err.get_msg());
699 0 : MyPrintError(err);
700 :
701 : }
702 :
703 0 : this->leave_debug_writer_section();
704 :
705 0 : if (!exec)
706 : {
707 : // Additional actions at failure
708 0 : UG_LOGN("Step "<< i<< " failed on stage " << i << ": reason not clear!" );
709 : }
710 :
711 : // switch to "parent" comm
712 : //mt_env.end();
713 : } /*for all stages loop*/
714 :
715 :
716 :
717 0 : return error;
718 : }
719 :
720 :
721 : template<class TDomain, class TAlgebra>
722 0 : void LimexTimeIntegrator<TDomain,TAlgebra>::join_integrator_threads()
723 : {
724 : // join all threads
725 : // g.join_all();
726 0 : const int nstages = m_vThreadData.size()-1;
727 0 : for (int i=nstages; i>=0; --i)
728 : {
729 0 : m_vThreadData[i].get_time_stepper()->invalidate();
730 : }
731 0 : }
732 :
733 :
734 : template<class TDomain, class TAlgebra>
735 0 : void LimexTimeIntegrator<TDomain,TAlgebra>::update_integrator_threads(ConstSmartPtr<grid_function_type> ucommon, number t)
736 : {
737 0 : const int nstages = m_vThreadData.size()-1;
738 0 : for (int i=nstages; i>=0; --i)
739 : {
740 : UG_ASSERT(m_vThreadData[i].get_solution()->size()==ucommon->size(), "LIMEX: Vectors must match in size!")
741 0 : *m_vThreadData[i].get_solution() = *ucommon;
742 : }
743 0 : }
744 :
745 : template<class TDomain, class TAlgebra>
746 0 : size_t LimexTimeIntegrator<TDomain,TAlgebra>::
747 : find_optimal_solution(const std::vector<number>& eps, size_t ntest, /*size_t &kf,*/ size_t &qpred)
748 : {
749 :
750 0 : const size_t qold=qpred;
751 :
752 : size_t jbest = 1;
753 0 : qpred = 1;
754 :
755 : size_t j=1;
756 : size_t k=j-1;
757 :
758 0 : m_lambda[k] = pow(m_rhoSafety*m_tol/eps[j], 1.0/m_gamma[k]); // 1/epsilon(k)
759 0 : m_workload[k] = m_costA[j]/m_lambda[k];
760 0 : UG_LOG("j=" << j << ": eps=" << eps[j] << ", lambda(j)=" <<m_lambda[k] << ", epsilon(j)=" <<1.0/m_lambda[k] << "<= alpha(k, qcurr)=" << monitor(k,qold-1) << "< alpha(k, qcurr+1)=" << monitor(k,qold) <<", A="<< m_costA[j] << ", W="<< m_workload[k] <<std::endl);
761 :
762 0 : for (j=2; j<ntest; ++j)
763 : {
764 0 : k = j-1;
765 0 : m_lambda[k] = pow(m_rhoSafety*m_tol/eps[j], 1.0/m_gamma[k]);
766 0 : m_workload[k] = m_costA[j]/m_lambda[k];
767 0 : UG_LOG("j=" << j << ": eps=" << eps[j] << ", lambda(j)=" <<m_lambda[k] << ", epsilon(j)=" <<1.0/m_lambda[k] << "<= alpha(k, qcurr)=" << monitor(k,qold-1) << "< alpha(k, qcurr+1)=" << monitor(k,qold) <<", A="<< m_costA[j] << ", W="<< m_workload[k] <<std::endl);
768 :
769 :
770 : // TODO: Convergence monitor
771 :
772 0 : qpred = (m_workload[qpred-1] > m_workload[k]) ? j : qpred;
773 0 : jbest = (eps[jbest] > eps [j]) ? j : jbest;
774 : }
775 :
776 0 : return jbest;
777 : }
778 :
779 : template<class TDomain, class TAlgebra>
780 0 : bool LimexTimeIntegrator<TDomain,TAlgebra>::
781 : apply(SmartPtr<grid_function_type> u, number t1, ConstSmartPtr<grid_function_type> u0, number t0)
782 : {
783 : PROFILE_FUNC_GROUP("limex");
784 : #ifdef UG_OPENMP
785 : // create multi-threading environment
786 : //int nt = std::min(omp_get_max_threads(), m_nstages);
787 :
788 : #endif
789 :
790 : // NOTE: we use u as common storage for future (and intermediate) solution(s)
791 0 : if (u.get() != u0.get()) // only copy if not already identical, otherwise: PST_UNDEFINED!
792 0 : *u = *u0;
793 :
794 : // initialize integrator threads
795 : // (w/ solutions)
796 0 : init_integrator_threads(u0);
797 :
798 :
799 : // write_debug
800 0 : for (unsigned int i=0; i<m_vThreadData.size(); ++i)
801 : {
802 0 : std::ostringstream ossName;
803 : ossName << std::setfill('0') << std::setw(4);
804 0 : ossName << "Limex_Init_iter" << 0 << "_stage" << i;
805 0 : this->write_debug(*m_vThreadData[i].get_solution(), ossName.str().c_str());
806 : }
807 :
808 : number t = t0;
809 0 : double dtcurr = ITimeIntegrator<TDomain, TAlgebra>::get_time_step();
810 :
811 0 : size_t kmax = m_vThreadData.size(); // maximum number of stages
812 0 : size_t qpred = kmax-1; // predicted optimal order
813 : size_t qcurr = qpred; // current order
814 :
815 : // for Gustafsson/lundh/Soederlind type PID controller
816 : /*size_t qlast = 0;
817 : double dtlast = 0.0;
818 : std::vector<double> epslast(kmax, m_rhoSafety*m_tol);
819 : */
820 :
821 : // time integration loop
822 : SmartPtr<grid_function_type> ubest = SPNULL;
823 : size_t limex_total = 1;
824 : size_t limex_success = 0;
825 : size_t ntest; ///< active number of stages <= kmax
826 : size_t jbest;
827 :
828 0 : m_bInterrupt = false;
829 : //bool bProbation = false;
830 : bool bAsymptoticReduction = false;
831 :
832 : const size_t nSwitchHistory=16;
833 : const size_t nSwitchLookBack=5;
834 0 : int vSwitchHistory[nSwitchHistory] ={ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0 };
835 :
836 0 : timex_type timex(m_vSteps);
837 0 : while ((t < t1) && ((t1-t) > base_type::m_precisionBound))
838 : {
839 : int err = 0;
840 :
841 : //UG_DLOG(LIB_LIMEX, 5, "+++ LimexTimestep +++" << limex_step << "\n");
842 0 : UG_LOG("+++ LimexTimestep +++" << m_limex_step << "\n");
843 :
844 :
845 : // save time stamp for limex step start
846 : Stopwatch stopwatch;
847 : stopwatch.start();
848 :
849 :
850 : // determine step size
851 0 : number dt = std::min(dtcurr, t1-t);
852 0 : UG_COND_THROW(dt < base_type::get_dt_min(), "Time step size below minimum. ABORTING! dt=" << dt << "; dt_min=" << base_type::get_dt_min() << "\n");
853 :
854 :
855 : // Notify init observers. (NOTE: u = u(t))¸
856 0 : itime_integrator_type::notify_init_step(u, m_limex_step, t, dt);
857 :
858 :
859 : // determine number of stages to investigate
860 0 : qcurr = qpred;
861 0 : ntest = std::min(kmax, qcurr+1);
862 : UG_LOG("ntest="<< ntest << std::endl);
863 :
864 : // checks
865 : UG_ASSERT(m_vSteps.size() >= ntest, "Huhh: sizes do not match: " << m_vSteps.size() << "<"<<ntest);
866 : UG_ASSERT(m_vThreadData.size() >= ntest, "Huhh: sizes do not match: " << m_vThreadData.size() << "< "<< ntest);
867 :
868 : ///////////////////////////////////////
869 : // PARALLEL EXECUTION: BEGIN
870 :
871 : // write_debug
872 0 : for (size_t i=0; i<ntest; ++i)
873 : {
874 0 : std::ostringstream ossName;
875 : ossName << std::setfill('0') << std::setw(4);
876 0 : ossName << "Limex_BeforeSerial_iter" << m_limex_step << "_stage" << i << "_total" << limex_total;
877 0 : this->write_debug(*m_vThreadData[i].get_solution(), ossName.str().c_str());
878 : }
879 0 : if(this->debug_writer_valid())
880 : {
881 0 : char debug_step_id_ext[16]; snprintf(debug_step_id_ext, 16, "%04d", m_limex_step);
882 0 : char debug_total_id_ext[16]; snprintf(debug_total_id_ext, 16, "%04d", (int) limex_total);
883 0 : this->enter_debug_writer_section(std::string("LimexTimeIntegrator_iter") + debug_step_id_ext + "_total" + debug_total_id_ext);
884 : }
885 :
886 : // integrate: t -> t+dt
887 0 : err = apply_integrator_threads(dt, u, t, ntest-1);
888 :
889 0 : this->leave_debug_writer_section();
890 :
891 : // write_debug
892 0 : for (size_t i=0; i<ntest; ++i)
893 : {
894 0 : std::ostringstream ossName;
895 : ossName << std::setfill('0') << std::setw(4);
896 0 : ossName << "Limex_AfterSerial_iter" << m_limex_step << "_stage" << i << "_total" << limex_total;
897 0 : this->write_debug(*m_vThreadData[i].get_solution(), ossName.str().c_str());
898 : }
899 :
900 0 : join_integrator_threads();
901 : // PARALLEL EXECUTION: END
902 : ///////////////////////////////////////
903 :
904 :
905 : ///////////////////////////////////////
906 : // SERIAL EXECUTION: BEGIN
907 :
908 : // sanity checks
909 : UG_ASSERT(m_spErrorEstimator.valid(), "Huhh: Invalid Error estimator?");
910 :
911 : double epsmin = 0.0;
912 :
913 : bool limexConverged = false;
914 0 : if (err==0)
915 : { // Compute extrapolation at t+dtcurr (SERIAL)
916 : // Consistency check.
917 0 : for (unsigned int i=1; i<ntest; ++i)
918 : {
919 0 : UG_LOG("Checking consistency:" <<m_consistency_error[i] <<"/" << m_consistency_error[i-1] << "="<< m_consistency_error[i]/m_consistency_error[i-1] << std::endl);
920 : }
921 :
922 : // Extrapolation.
923 0 : timex.set_error_estimate(m_spErrorEstimator);
924 0 : for (unsigned int i=0; i<ntest; ++i)
925 : {
926 : UG_ASSERT(m_vThreadData[i].get_solution().valid(), "Huhh: no valid solution?");
927 0 : timex.set_solution(m_vThreadData[i].get_solution(), i);
928 : }
929 0 : timex.apply(ntest);
930 :
931 : // write_debug
932 0 : for (size_t i=0; i<ntest; ++i)
933 : {
934 0 : std::ostringstream ossName;
935 : ossName << std::setfill('0') << std::setw(4);
936 0 : ossName << "Limex_Extrapolates_iter" << m_limex_step << "_stage" << i << "_total" << limex_total;
937 0 : this->write_debug(*m_vThreadData[i].get_solution(), ossName.str().c_str());
938 : }
939 0 : limex_total++;
940 :
941 : // Obtain sub-diagonal error estimates.
942 : const std::vector<number>& eps = timex.get_error_estimates();
943 : UG_ASSERT(ntest<=eps.size(), "Huhh: Not enough solutions?");
944 :
945 : // select optimal solution (w.r.t error) AND
946 : // predict optimal order (w.r.t. workload) for next step
947 0 : jbest = find_optimal_solution(eps, ntest, qpred);
948 : UG_ASSERT(jbest < ntest, "Huhh: Not enough solutions?");
949 :
950 : // best solution
951 0 : ubest = timex.get_solution(jbest-m_conservative).template cast_dynamic<grid_function_type>();
952 0 : epsmin = eps[jbest];
953 :
954 : // check for convergence
955 0 : limexConverged = (epsmin <= m_tol) ;
956 :
957 :
958 0 : if (limex_success>3)
959 : {
960 :
961 0 : vSwitchHistory[m_limex_step%nSwitchHistory] = (qpred - qcurr);
962 : UG_DLOG(LIB_LIMEX, 5, "LIMEX-ASYMPTOTIC-ORDER switch: = " << (qpred - qcurr)<< std::endl);
963 :
964 : size_t nSwitches=0;
965 0 : for (int s=nSwitchLookBack-1; s>=0; s--)
966 : {
967 0 : nSwitches += std::abs(vSwitchHistory[(m_limex_step-s)%nSwitchHistory]);
968 : UG_DLOG(LIB_LIMEX, 6, "LIMEX-ASYMPTOTIC-ORDER: s[" << s<< "] = " << vSwitchHistory[(m_limex_step-s)%nSwitchHistory] << std::endl);
969 : }
970 : UG_DLOG(LIB_LIMEX, 5,"LIMEX-ASYMPTOTIC-ORDER: nSwitches = " << nSwitches << std::endl);
971 :
972 :
973 : /*
974 : if (bProbation)
975 : {
976 :
977 : if ((qpred < qcurr) || (!limexConverged)) // consecutive (follow-up) order decrease
978 : {
979 : bProbation = true;
980 : m_num_reductions[0]++;
981 : UG_LOG("LIMEX-ASYMPTOTIC-ORDER: Decrease on parole detected: "<< qcurr << " -> " << qpred << "("<< m_num_reductions[0]<<")"<<std::endl);
982 :
983 : }
984 : else
985 : {
986 : // reset all counters
987 : bProbation = false;
988 : UG_LOG("LIMEX-ASYMPTOTIC-ORDER: Probation off!"<<std::endl);
989 :
990 : }
991 : }
992 : else
993 : {
994 : if ((qpred < qcurr) || (!limexConverged))// 1st order decrease
995 : {
996 : bProbation = true;
997 : m_num_reductions[0]++;
998 : UG_LOG("LIMEX-ASYMPTOTIC-ORDER: Decrease detected: "<< qcurr << " -> " << qpred << "("<< m_num_reductions[0]<<")"<<std::endl);
999 : }
1000 : else
1001 : {
1002 : m_num_reductions[0] = 0;
1003 : UG_LOG("LIMEX-ASYMPTOTIC-ORDER: I am totally free!"<<std::endl);
1004 : }
1005 :
1006 :
1007 : }
1008 : */
1009 :
1010 : // bAsymptoticReduction = (m_num_reductions[0] >= m_max_reductions) || bAsymptoticReduction;
1011 :
1012 : // bAsymptoticReduction = (nSwitches >= m_max_reductions) || bAsymptoticReduction;
1013 :
1014 0 : if (nSwitches >= m_max_reductions) {
1015 :
1016 :
1017 0 : m_asymptotic_order = std::min<size_t>(m_asymptotic_order-1, 2);
1018 : bAsymptoticReduction = true;
1019 :
1020 0 : for (int s=nSwitchLookBack-1; s>=0; s--)
1021 0 : { vSwitchHistory[(m_limex_step-s)%nSwitchHistory]=0; }
1022 :
1023 :
1024 : }
1025 :
1026 :
1027 :
1028 : // asymptotic order reduction
1029 : //if (m_num_reductions[0] >= m_max_reductions)
1030 0 : if (bAsymptoticReduction)
1031 : {
1032 0 : kmax = (kmax >= m_asymptotic_order ) ? m_asymptotic_order : kmax;
1033 0 : qpred = kmax - 1;
1034 : UG_DLOG(LIB_LIMEX, 5, "LIMEX-ASYMPTOTIC-ORDER: Reduction: "<< qpred );
1035 : UG_DLOG(LIB_LIMEX, 5, "(kmax=" << kmax << ", asymptotic"<< m_asymptotic_order <<") after "<<m_max_reductions<< std::endl);
1036 : }
1037 :
1038 : } // (limex_success>3)
1039 :
1040 : /*
1041 : // adjust for order bound history
1042 : if ((m_num_reductions[qpred] > m_max_reductions) && (qpred > 1))
1043 : {
1044 : UG_LOG("prohibited q:"<< qpred << "(" << m_num_reductions[qpred] << ")"<<std::endl);
1045 : //qpred--;
1046 : qpred = kmax = 2;
1047 : } else
1048 : {
1049 : UG_LOG("keeping q:"<< qpred << "(" << m_num_reductions[qpred] << ")"<<std::endl);
1050 : }
1051 : */
1052 :
1053 : //double pid = 1.0;
1054 :
1055 : // constant order?
1056 : /*if (qcurr == qlast)
1057 : {
1058 : double dtRatio =dtcurr/dtlast;
1059 : // Dd, Band 3 (DE), p. 360
1060 : int k = ntest-1;
1061 : while (k--) {
1062 : double epsRatio =eps[k+1]/epslast[k+1];
1063 :
1064 : double qopt = log(epsRatio) / log(dtRatio); // take dtcurr here, as dt may be smaller
1065 :
1066 : UG_LOG("effective p["<< k << "]="<< qopt-1.0 <<","<< eps[k+1] <<","<< epslast[k+1] <<"," << epsRatio << "("<< log(epsRatio) << ","<< pow(epsRatio, 1.0/m_gamma[k]) <<")"<< dtcurr <<","<< dtlast << "," << dtRatio << "," <<log(dtRatio) << std::endl);
1067 : }
1068 :
1069 :
1070 : double epsRatio = eps[qcurr]/epslast[qcurr];
1071 : pid = dtRatio/pow(epsRatio, 1.0/m_gamma[qcurr-1]);
1072 :
1073 : UG_LOG("pid=" << pid << std::endl);
1074 :
1075 : }*/
1076 :
1077 :
1078 : // select (predicted) order for next step//double dtpred = dtcurr*std::min(m_lambda[qpred-1], itime_integrator_type::get_increase_factor());
1079 0 : double dtpred = dtcurr*std::min(m_lambda[qpred-1], itime_integrator_type::get_increase_factor());
1080 : //double dtpred = dtcurr*m_lambda[qpred-1];
1081 : UG_LOG("+++++\nget_increase_factor() gives "<<itime_integrator_type::get_increase_factor()<<" \n+++++++")
1082 0 : UG_LOG("koptim=\t" << jbest << ",\t eps(k)=" << epsmin << ",\t q=\t" << qpred<< "("<< ntest << "), lambda(q)=" << m_lambda[qpred-1] << ", alpha(q-1,q)=" << monitor(qpred-1, qpred) << "dt(q)=" << dtpred<< std::endl);
1083 :
1084 : // EXTENSIONS: convergence model
1085 0 : if (limexConverged)
1086 : {
1087 0 : limex_success++;
1088 :
1089 : // a) aim for order increase in next step
1090 0 : if ((qpred+1==ntest) /* increase by one possible? */
1091 : // && (m_lambda[qpred-1]>1.0)
1092 : // && (m_num_reductions[qpred+1] <= m_max_reductions)
1093 0 : && (kmax>ntest)) /* still below max? */
1094 : {
1095 0 : const double alpha = monitor(qpred-1, qpred);
1096 0 : UG_LOG("CHECKING for order increase: "<< m_costA[qpred] << "*" << alpha << ">" << m_costA[qpred+1]);
1097 : // check, whether further increase could still be efficient
1098 0 : if (m_costA[qpred] * alpha > m_costA[qpred+1])
1099 : {
1100 0 : qpred++; // go for higher order
1101 0 : if (m_greedyOrderIncrease >0.0) {
1102 0 : dtpred *= (1.0-m_greedyOrderIncrease) + m_greedyOrderIncrease*alpha; // & adapt time step // TODO: check required!
1103 : }
1104 : UG_LOG("... yes.\n")
1105 :
1106 : // update history
1107 0 : vSwitchHistory[m_limex_step%nSwitchHistory] = (qpred - qcurr);
1108 0 : UG_LOG("LIMEX-ASYMPTOTIC-ORDER switch update: = " << (qpred - qcurr)<< std::endl);
1109 :
1110 : } else {
1111 : UG_LOG("... nope.\n")
1112 : }
1113 :
1114 : }
1115 :
1116 :
1117 : // b) monitor convergence (a-priori check!)
1118 :
1119 : }
1120 :
1121 : // parameters for subsequent step
1122 0 : dtcurr = std::min(dtpred, itime_integrator_type::get_dt_max());
1123 :
1124 :
1125 : }
1126 : else
1127 : {
1128 : // solver failed -> cut time step //
1129 : number base_dtmin=base_type::get_dt_min();
1130 0 : if(dtcurr <= base_dtmin)
1131 0 : base_type::set_dt_min(base_dtmin*m_sigmaReduction);
1132 :
1133 0 : dtcurr=std::max(base_type::get_dt_min(), dtcurr*m_sigmaReduction);
1134 :
1135 : //dtcurr=std:max(base_type::get_dt_min(), dtcurr);
1136 : }
1137 :
1138 : // output compute time for Limex step
1139 0 : number watchTime = stopwatch.ms()/1000.0;
1140 : UG_LOGN("Time: " << std::setprecision(3) << watchTime << "s");
1141 :
1142 0 : if ((err==0) && limexConverged)
1143 : {
1144 : // ACCEPT time step
1145 0 : UG_LOG("+++ LimexTimestep +++" << m_limex_step << " ACCEPTED"<< std::endl);
1146 : UG_LOG(" :\t time \t dt (success) \t dt (pred) \tq=\t order (curr)" << qcurr+1 << std::endl);
1147 0 : UG_LOG("LIMEX-ACCEPTING:\t" << t <<"\t"<< dt << "\t" << dtcurr << "\tq=\t" << qcurr+1 << std::endl);
1148 :
1149 : // update PID controller
1150 : /*qlast = qcurr;
1151 : epslast = timex.get_error_estimates(); // double check ???
1152 : dtlast = dt;*/
1153 : // Compute time derivatives (by extrapolation)
1154 0 : if (this->has_time_derivative())
1155 : {
1156 : UG_LOG("Computing derivative" << std::endl);
1157 0 : grid_function_type &udot = *get_time_derivative();
1158 :
1159 0 : for (size_t i=0; i<=jbest; ++i)
1160 : {
1161 0 : timex.set_solution(m_vThreadData[i].get_derivative(), i);
1162 : }
1163 0 : timex.apply(jbest+1, false); // do not compute error
1164 :
1165 0 : udot = *timex.get_solution(jbest).template cast_dynamic<grid_function_type>();
1166 :
1167 0 : std::ostringstream ossName;
1168 : ossName << std::setfill('0');
1169 0 : ossName << "Limex_Derivative_iter" << m_limex_step << "_total" << limex_total;
1170 0 : this->write_debug(udot, ossName.str().c_str());
1171 0 : }
1172 :
1173 :
1174 : // Copy best solution.
1175 : UG_ASSERT(ubest.valid(), "Huhh: Invalid error estimate?");
1176 0 : *u = *ubest;
1177 0 : t += dt;
1178 :
1179 : // Initiate post process.
1180 0 : itime_integrator_type::notify_finalize_step(u, m_limex_step++, t, dt);
1181 :
1182 : // make sure that all threads continue
1183 : // with identical initial value u(t)
1184 : // update_integrator_threads(ubest, t);
1185 :
1186 :
1187 : // working on last row => increase order
1188 : //if (ntest == q+1) ntest++;
1189 : }
1190 : else
1191 : {
1192 : // DISCARD time step
1193 0 : UG_LOG("+++ LimexTimestep +++" << m_limex_step << " FAILED" << std::endl);
1194 0 : UG_LOG(" :\t time \t dt (failed) \t dt (curr) \teps=\t" << epsmin <<"\t(tol="<<m_tol<<")\terr="<<err<< std::endl);
1195 0 : UG_LOG("LIMEX-REJECTING:\t" << t <<"\t"<< dt << "\t" << dtcurr <<std::endl);
1196 :
1197 0 : itime_integrator_type::notify_rewind_step(ubest, m_limex_step, t+dt, dt);
1198 :
1199 : }
1200 :
1201 : // interrupt if interruption flag set
1202 0 : if (m_bInterrupt)
1203 : {
1204 : UG_LOGN("Limex interrupted by external command.");
1205 0 : break;
1206 : }
1207 :
1208 : // SERIAL EXECUTION: END
1209 : ///////////////////////////////////////
1210 0 : update_integrator_threads(u, t);
1211 :
1212 :
1213 : // SOLVE
1214 :
1215 : // ESTIMATE
1216 :
1217 : // MARK
1218 :
1219 : // REFINE
1220 :
1221 :
1222 : } // time integration loop
1223 :
1224 0 : dispose_integrator_threads();
1225 :
1226 0 : return true;
1227 0 : } // apply
1228 :
1229 :
1230 : } // namespace ug
1231 :
1232 : #endif /* TIME_INTEGRATOR_HPP_ */
|
