Plugins
Loading...
Searching...
No Matches
time_integrator.hpp
Go to the documentation of this file.
1/*
2 * Copyright (c) 2014-2020: G-CSC, Goethe University Frankfurt
3 * Author: Arne Naegel, Andreas Kreienbuehl
4 *
5 * This file is part of UG4.
6 *
7 * UG4 is free software: you can redistribute it and/or modify it under the
8 * terms of the GNU Lesser General Public License version 3 (as published by the
9 * Free Software Foundation) with the following additional attribution
10 * requirements (according to LGPL/GPL v3 §7):
11 *
12 * (1) The following notice must be displayed in the Appropriate Legal Notices
13 * of covered and combined works: "Based on UG4 (www.ug4.org/license)".
14 *
15 * (2) The following notice must be displayed at a prominent place in the
16 * terminal output of covered works: "Based on UG4 (www.ug4.org/license)".
17 *
18 * (3) The following bibliography is recommended for citation and must be
19 * preserved in all covered files:
20 * "Reiter, S., Vogel, A., Heppner, I., Rupp, M., and Wittum, G. A massively
21 * parallel geometric multigrid solver on hierarchically distributed grids.
22 * Computing and visualization in science 16, 4 (2013), 151-164"
23 * "Vogel, A., Reiter, S., Rupp, M., Nägel, A., and Wittum, G. UG4 -- a novel
24 * flexible software system for simulating pde based models on high performance
25 * computers. Computing and visualization in science 16, 4 (2013), 165-179"
26 *
27 * This program is distributed in the hope that it will be useful,
28 * but WITHOUT ANY WARRANTY; without even the implied warranty of
29 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
30 * GNU Lesser General Public License for more details.
31 */
32
33#ifndef TIME_INTEGRATOR_HPP_
34#define TIME_INTEGRATOR_HPP_
35
36#if __cplusplus >= 201103L
37#define OVERRIDE override
38#else
39#define OVERRIDE
40#endif
41
42// std headers.
43#include <string>
44
45// UG4 headers.
50#include "lib_disc/assemble_interface.h" // TODO: missing IAssemble in following file:
59#include "lib_disc/function_spaces/grid_function_util.h" // SaveVectorForConnectionViewer
60#include "lib_disc/function_spaces/interpolate.h" //Interpolate
61#include "lib_disc/function_spaces/integrate.h" //Integral
62#include "lib_disc/function_spaces/grid_function.h" //GridFunction
64
65
66
67// My headers.
68#include "time_extrapolation.h"
69#include "../limex_tools.h"
70
71namespace ug {
72
73
75template<class TDomain, class TAlgebra>
77: public ITimeIntegratorObserver<TDomain, TAlgebra>,
78 public ITimeIntegratorStageObserver_start<TDomain, TAlgebra>,
79 public ITimeIntegratorStageObserver_finalize<TDomain, TAlgebra>,
80 public ITimeIntegratorStageObserver_end<TDomain, TAlgebra>
81{
82public:
86
89
90 VTKOutputObserver(const char *filename, SmartPtr<vtk_type> vtk)
91 : m_sp_vtk(vtk), m_filename(filename), m_startTime(0.0), m_plotStep(0.0) {}
92
93 VTKOutputObserver(const char *filename, SmartPtr<vtk_type> vtk, number plotStep)
94 : m_sp_vtk(vtk), m_filename(filename), m_startTime(0.0), m_plotStep(plotStep) {}
95
97 { m_sp_vtk = SPNULL; }
98
99
100 void set_output_scales(const std::vector<number>& vScales)
101 {
102 m_vOutputScales = vScales;
103 }
104
105
107 {
108 return true;
109 }
110
111
113 {
114 if (!m_sp_vtk.valid())
115 return false;
116
117 writeToFile(u, step, time);
118
119 m_startTime = time;
120 m_uOld = u->clone();
121
122 return true;
123 }
124
125
127 {
128 if (!m_sp_vtk.valid())
129 return false;
130
131 if (m_plotStep == 0.0)
132 {
133 writeToFile(uNew, step, time);
134 return true;
135 }
136
137 // otherwise, only plot data at multiples of given time step (interpolate if necessary)
138 number rem = fmod(time - dt - m_startTime, m_plotStep);
139 number nextPlotTimePt = time - dt - rem + m_plotStep;
140 int nextStep = (int) ((nextPlotTimePt - m_startTime + 0.5 * m_plotStep) / m_plotStep);
141
142 if (nextPlotTimePt > time)
143 {
144 m_uOld = uNew->clone();
145 return true;
146 }
147
148 SmartPtr<grid_function_type> uCur = uNew->clone_without_values();
149 while (nextPlotTimePt <= time)
150 {
151 number alpha = (time - nextPlotTimePt) / dt;
152 VecScaleAdd(static_cast<typename TAlgebra::vector_type&>(*uCur),
153 alpha, static_cast<typename TAlgebra::vector_type&>(*m_uOld),
154 1.0 - alpha, static_cast<typename TAlgebra::vector_type&>(*uNew));
155
156 writeToFile(uCur, nextStep, nextPlotTimePt);
157
158 nextPlotTimePt = (++nextStep) * m_plotStep;
159 }
160
161 m_uOld = uNew->clone();
162 return true;
163 }
164
165
167 {
168 if (!m_sp_vtk.valid())
169 return false;
170
171 m_sp_vtk->write_time_pvd(m_filename.c_str(), *u);
172 return true;
173 }
174
175
176private:
178 {
179 if (m_vOutputScales.size())
180 {
181 SmartPtr<grid_function_type> uTmp = u->clone_without_values();
182 ScaleGF<grid_function_type>(uTmp, u, m_vOutputScales);
183 m_sp_vtk->print(m_filename.c_str(), *uTmp, step, time);
184 }
185 else
186 m_sp_vtk->print(m_filename.c_str(), *u, step, time);
187 }
188
189protected:
192 std::string m_filename;
195 std::vector<number> m_vOutputScales;
196};
197
198
199
201template<class TDomain, class TAlgebra>
203: public ITimeIntegratorObserver<TDomain, TAlgebra>
204{
205public:
208
209 ConnectionViewerOutputObserver(const char *filename)
210 : m_filename(filename), m_outputTime(-1.0) {}
211
212 ConnectionViewerOutputObserver(const char *filename, number t_out)
213 : m_filename(filename), m_outputTime(t_out) {}
214
217
218 bool step_process(SmartPtr<grid_function_type> uNew, /*SmartPtr<grid_function_type> uOld, */int step, number time, number dt) OVERRIDE
219 {
220 // quit, if time does not match
221 if (m_outputTime >=0.0 && time != m_outputTime) return true;
222
223 SaveVectorForConnectionViewer<grid_function_type>(*uNew, m_filename.c_str());
224 return true;
225 }
226
227protected:
228 std::string m_filename;
230
231};
232
233
234#if 0
235template <typename TData, typename TDataIn1, typename TDataIn2>
236class LuaFunction2 // : public IFunction<TData, TDataIn1, typename TDataIn2>
237{
238 public:
240 LuaFunction2();
241 virtual ~LuaFunction2() {};
242
244
249 void set_lua_callback(const char* luaCallback, size_t numArgs);
250
252 virtual void operator() (TData& out, int numArgs1, int numArgs2,...);
253
254 protected:
256 std::string m_cbValueName;
257
259 int m_cbValueRef;
260
262 lua_State* m_L;
263
265 size_t m_numArgs;
266};
267
268
269template <typename TData, typename TDataIn1, typename TDataIn2>
270LuaFunction2<TData,TDataIn1,TDataIn2>::LuaFunction2() : m_numArgs(0)
271{
273 m_cbValueRef = LUA_NOREF;
274}
275
276template <typename TData, typename TDataIn1, typename TDataIn2>
277void LuaFunction2<TData,TDataIn1,TDataIn2>::set_lua_callback(const char* luaCallback, size_t numArgs)
278{
279// store name (string) of callback
280 m_cbValueName = luaCallback;
281
282// obtain a reference
283 lua_getglobal(m_L, m_cbValueName.c_str());
284
285// make sure that the reference is valid
286 if(lua_isnil(m_L, -1)){
287 UG_THROW("LuaFunction::set_lua_callback(...):"
288 "Specified lua callback does not exist: " << m_cbValueName);
289 }
290
291// store reference to lua function
292 m_cbValueRef = luaL_ref(m_L, LUA_REGISTRYINDEX);
293
294// remember number of arguments to be used
295 m_numArgs = numArgs;
296}
297#endif
298
299/*
300SmartUserDataWrapper* CreateNewUserData(lua_State* L, const SmartPtr<void>& ptr,
301 const char* metatableName)
302{
303// create the userdata
304 SmartUserDataWrapper* udata = (SmartUserDataWrapper*)lua_newuserdata(L,
305 sizeof(SmartUserDataWrapper));
306 new(udata) SmartUserDataWrapper;
307
308// associate the object with the userdata.
309 udata->smartPtr = ptr;
310 udata->type = SMART_POINTER;
311
312// associate the metatable (userdata is already on the stack)
313 luaL_getmetatable(L, metatableName);
314 lua_setmetatable(L, -2);
315
316 return udata;
317}
318*/
319/*
320template <typename TData, typename TDataIn1, typename TDataIn2>
321void LuaFunction2<TData,TDataIn1,TDataIn2>::operator() (TData& out, int numArgs1, SmartPtr<TDataIn1> valsArgs1[],
322 int numArgs2, ...)
323{
324 PROFILE_CALLBACK_BEGIN(operatorBracket);
325
326 UG_ASSERT((numArgs1+numArgs2) == (int)m_numArgs, "Number of arguments mismatched.");
327
328 // push the callback function on the stack
329 lua_rawgeti(m_L, LUA_REGISTRYINDEX, m_cbValueRef);
330
331 // get list of arguments
332 va_list ap;
333 va_start(ap, numArgs2);
334
335 // read all arguments and push them to the lua stack
336 for(int i = 0; i < numArgs1; ++i)
337 {
338
339 CreateNewUserData(m_L, &valArgs1[i], "");
340
341 }
342 for(int i = 0; i < numArgs2; ++i)
343 {
344 TDataIn2 val = va_arg(ap, TDataIn2);
345 lua_traits<TDataIn2>::push(m_L, val);
346 }
347
348
349 // end read in of parameters
350 va_end(ap);
351
352 // compute total args size
353 size_t argSize = lua_traits<TDataIn1>::size * numArgs1;
354 argSize += lua_traits<TDataIn2>::size * numArgs2;
355
356 // compute total return size
357 size_t retSize = lua_traits<TData>::size;
358
359 // call lua function
360 if(lua_pcall(m_L, argSize, retSize, 0) != 0)
361 UG_THROW("LuaFunction::operator(...): Error while "
362 "running callback '" << m_cbValueName << "',"
363 " lua message: "<< lua_tostring(m_L, -1));
364
365 try{
366 // read return value
367 lua_traits<TData>::read(m_L, out);
368 UG_COND_THROW(IsFiniteAndNotTooBig(out)==false, out);
369 }
370 UG_CATCH_THROW("LuaFunction::operator(...): Error while running "
371 "callback '" << m_cbValueName << "'");
372
373 // pop values
374 lua_pop(m_L, retSize);
375
376 PROFILE_CALLBACK_END();
377}
378*/
379
380template<class TDomain, class TAlgebra>
382: public ITimeIntegratorObserver<TDomain, TAlgebra>
383{
384public:
388
391
392#ifdef UG_FOR_LUA
393 PlotRefOutputObserver(const char *ExactSol)
396
397 PlotRefOutputObserver(const char *ExactSol, SmartPtr<vtk_type> vtk)
398 : m_sp_vtk(vtk)
399 { m_spReference = make_sp(new LuaUserData<number, grid_function_type::dim>(ExactSol)); }
400
401#endif
402
404 {}
405
406 // TODO: replace by call 'func (SmartPtr<G> u, int step, number dt, number t)'
407 bool step_process(SmartPtr<grid_function_type> uNew, /* SmartPtr<grid_function_type> uOld, */ int step, number time, number dt) OVERRIDE
408 {
409 UG_LOG("L2Error(\t"<< time << "\t) = \t" << L2Error(m_spReference, uNew, "c", time, 4) << std::endl);
410 if (m_sp_vtk.valid())
411 {
412 SmartPtr<grid_function_type> ref = uNew->clone();
413 Interpolate<grid_function_type> (m_spReference, ref, "c", time);
414 m_sp_vtk->print("MyReference", *ref, step, time);
415 return true;
416 }
417
418 return false;
419
420 }
421
422protected:
423 // TODO: replace by appropriate call-back
426};
427
429
432template<class TDomain, class TAlgebra>
434: public ITimeIntegratorObserver<TDomain, TAlgebra>
435{
436public:
440protected:
442 {
443 IntegralSpecs(const char* cmp, const char* subsets, int quadOrder, const char *idString) :
444 m_cmp(cmp), m_subsets(subsets), m_quadOrder(quadOrder), m_idString(idString)
445 {};
446 std::string m_cmp;
447 std::string m_subsets;
449 std::string m_idString;
450 };
451
452public:
455
458
459 // TODO: replace by call 'func (SmartPtr<G> u, int step, number dt, number t)'
460 bool step_process(SmartPtr<grid_function_type> uNew, /*SmartPtr<grid_function_type> uOld,*/ int step, number time, number dt) OVERRIDE
461 {
462
463 for (typename std::vector<IntegralSpecs>::iterator it = m_vIntegralData.begin();
464 it!=m_vIntegralData.end(); ++it)
465 {
466 number value=Integral(uNew, it->m_cmp.c_str(), it->m_subsets.c_str(), it->m_quadOrder);
467 UG_LOG("Integral(\t"<< it->m_idString << "\t"<< time << "\t)=\t" << value << std::endl);
468 }
469
470 return true;
471
472
473 }
474
475
476 void add_integral_specs(const char* cmp, const char* subsets, int quadOrder, const char* idString)
477 {
478 m_vIntegralData.push_back(IntegralSpecs(cmp, subsets, quadOrder, idString));
479 }
480
481protected:
482
483 std::vector<IntegralSpecs> m_vIntegralData;
484 //const char* cmp,
485 // const char* subsets,
486 // int quadOrder
487};
488
489
490
492template<class TDomain, class TAlgebra>
494 : public IOperator< GridFunction<TDomain, TAlgebra> >,
495 public TimeIntegratorSubject<TDomain, TAlgebra>
496{
497 public:
498
499 typedef TAlgebra algebra_type;
500 typedef typename TAlgebra::vector_type vector_type;
501 typedef typename TAlgebra::matrix_type matrix_type;
502
503
504 typedef TDomain domain_type;
506
507 protected:
508 double m_dt;
511
513
515
516
517 public:
518 // constructor
520 : m_dt(1.0), m_lower_tim(0.0), m_upper_tim(0.0), m_precisionBound(1e-10), m_bNoLogOut(false)
521 {}
522
524 virtual ~ITimeIntegrator() {};
525
527
538 virtual void init(grid_function_type const& u)
539 {
540 // UG_ASSERT(m_spDomainDisc.valid(), "TimeIntegrator<TDomain, TAlgebra>::init: m_spDomainDisc invalid.");
541 // m_spTimeDisc = make_sp(new time_disc_type(m_spDomainDisc, m_theta));
542 }
543
544
546
551 void init()
552 { UG_THROW("Please init with grid function!"); }
553
555
556 /* This method is used to prepare the in- and output vector used later in apply.
557 * It can be called, after the init method has been called at least once.
558 * The prepare method is e.g. used to set dirichlet values.
559 */
561
563
565 { UG_THROW("Fix interfaces!"); }
566
568
570 void set_time_step(double dt)
571 { m_dt = dt; }
572
574 { return m_dt; }
575
576 void set_precision_bound(double precisionBound)
577 { m_precisionBound = precisionBound; return; }
578
579 void set_no_log_out(bool bNoLogOut)
580 { m_bNoLogOut = bNoLogOut; return; }
581
582
583
584};
585
586
587
589template<class TDomain, class TAlgebra>
615
616
618template<class TAlgebra>
632
633
635template<class TDomain, class TAlgebra>
660
661
662
663template<typename TDomain, typename TAlgebra>
665{
666
668 UG_COND_THROW(!base_type::m_spLinearSolver.valid(), "Linear solver invalid");
669
670 // short-cuts
671 GridLevel const &gl = u0->grid_level();
672 time_disc_type &tdisc = *tdisc_dep_type::m_spTimeDisc;
673
674 // create solution vector & right hand side
676 SmartPtr<typename base_type::vector_type> b= u0->clone_without_values();
677
678 // solution time series
679 SmartPtr<vector_time_series_type> m_spSolTimeSeries;
680 m_spSolTimeSeries=make_sp(new vector_time_series_type());
681 m_spSolTimeSeries->clear();
682 m_spSolTimeSeries->push(uold, t0);
683
684 // create matrix operator
685 SmartPtr<typename base_type::assembled_operator_type> spAssOp=make_sp(new typename base_type::assembled_operator_type(tdisc_dep_type::m_spTimeDisc, gl));
686
687 // integrate
688 if(!base_type::m_bNoLogOut)
689 UG_LOG("+++ Integrating: ["<< t0 <<", "<< t1 <<"]\n");
690
691 double t = t0;
692 number dt_assembled = -1.0; // invalid
693 int step = 1;
694
695 number currdt = base_type::m_dt;
696
697 while((t < t1) && (t1-t > base_type::m_precisionBound))
698 {
699
700 if(!base_type::m_bNoLogOut)
701 { UG_LOG("+++ Timestep +++" << step << "\n"); }
702
703 // determine step size
704 number dt = std::min(currdt, t1-t);
705
706 // prepare step
707 tdisc.prepare_step(m_spSolTimeSeries, dt);
708 if (fabs(dt-dt_assembled) > base_type::m_precisionBound)
709 {
710 // re-assemble operator
711 if(!base_type::m_bNoLogOut)
712 UG_LOG("+++ Reassemble (t=" << t << ", dt=" << dt <<")\n");
713
714 tdisc.assemble_linear(*spAssOp, *b, gl);
715 (base_type::m_spLinearSolver)->init(spAssOp, *u1);
716 dt_assembled = dt;
717 }
718 else
719 {
720 // keep old operator
721 tdisc.assemble_rhs(*b, gl);
722 }
723
724 // execute step
725 if (base_type::m_spLinearSolver->apply(*u1, *b))
726 {
727 // ACCEPTING:
728 // push updated solution into time series
729 t += dt;
730 SmartPtr<typename base_type::vector_type> tmp = m_spSolTimeSeries->oldest();
731 VecAssign(*tmp, *u1.template cast_dynamic<typename base_type::vector_type>());
732 m_spSolTimeSeries->push_discard_oldest(tmp, t);
733
734 this->notify_finalize_step(u1, step++, t+dt, dt);
735 }
736 else
737 {
738 // DISCARDING
739 currdt *= 0.5;
740 }
741
742 }
743
744 this->notify_end(u1, step, t1, currdt);
745
746 return true;
747};
748
749
751template<class TDomain, class TAlgebra>
783
784
785
786template<typename TDomain, typename TAlgebra>
788{
790 UG_COND_THROW(!base_type::m_spLinearSolver.valid(), "Linear solver invalid");
791
792 // short-cuts
793 GridLevel const &gl = u0->grid_level();
794 time_disc_type &tdisc = *tdisc_dep_type::m_spTimeDisc;
795
796 // create solution vector & right hand side
798 SmartPtr<typename base_type::vector_type> b= u0->clone_without_values();
799
800 // solution time series
801 SmartPtr<vector_time_series_type> m_spSolTimeSeries;
802 m_spSolTimeSeries=make_sp(new vector_time_series_type());
803 m_spSolTimeSeries->push(uold, t0);
804
806
807 // select number of steps
808 double t = t0;
809 int numSteps = round((t1-t0) / base_type::m_dt);
810 number currdt = (t1-t0) / numSteps;
811
812
813
814 //std::cerr << "+++ Integrating: ["<< t0 <<", "<< t1 <<"] with dt=" << currdt << "("<< numSteps<< " iters)\n";
815 if(!base_type::m_bNoLogOut)
816 {
817 UG_LOG("+++ Integrating: [\t"<< t0 <<"\t, \t"<< t1 <<"\t] with dt=\t" << currdt << "("<< numSteps<< " iters)" << std::endl);
818 }
819
820 // integrate
821 for(int step = 1; step<=numSteps; ++step)
822 {
823 // determine step size
824 // number dt = std::min(currdt, t1-t);
825 const number dt = currdt;
826
827 if(!base_type::m_bNoLogOut)
828 {
829 UG_LOG("+++ Const timestep +++" << step<< "(t=" << t << ", dt=" << dt << ")"<< std::endl);
830 }
831 this->notify_init_step(u1, step, t, dt);
832
833 // prepare step
834 tdisc.prepare_step(m_spSolTimeSeries, dt);
835 if (spAssOp==SPNULL)
836 {
837 // Assemble operator.
838 if(!base_type::m_bNoLogOut) UG_LOG("+++ Assemble (t=" << t << ", dt=" << dt <<")" << std::endl);
839
840 spAssOp=make_sp(new typename base_type::assembled_operator_type(tdisc_dep_type::m_spTimeDisc, gl));
841 tdisc.assemble_linear(*spAssOp, *b, gl);
842 (base_type::m_spLinearSolver)->init(spAssOp, *u1);
843 }
844 else
845 {
846 // Recycle existing operator.
847 // std::cerr << "Recycling timestep " << step << "\n";
848 tdisc.assemble_rhs(*b, gl);
849 }
850
851 // execute step
852 if (base_type::m_spLinearSolver->apply(*u1, *b))
853 {
854 // ACCEPTING: push updated solution into time series
855 t += dt;
856 SmartPtr<typename base_type::vector_type> tmp = m_spSolTimeSeries->oldest();
857 VecAssign(*tmp, *u1.template cast_dynamic<typename base_type::vector_type>());
858 m_spSolTimeSeries->push_discard_oldest(tmp, t);
859 this->notify_finalize_step(u1, step, t, dt); // t updated already!
860 }
861 else
862 {
863 UG_THROW("ConstStepLinearTimeIntegrator::apply failed!!!");
864 this->notify_rewind_step(u1, step, t+dt, dt);
865 }
866
867 }
868
869 this->notify_end(u1, numSteps, t1, currdt);
870
871 return true;
872};
873
874
875
877template<class TDomain, class TAlgebra>
914
915
916template<typename TDomain, typename TAlgebra>
918{
919 // call base
920 base_type::init(u);
921 //m_spTimeDisc2 = make_sp(new typename base_type::time_disc_type(base_type::m_spDomainDisc, base_type::m_theta));
922}
923
924template<typename TDomain, typename TAlgebra>
926{
927 // short-cuts
928 GridLevel const &gl = u0->grid_level();
929 time_disc_type &tdisc = *tdisc_dep_type::m_spTimeDisc;
930 time_disc_type &tdisc2 = *m_spTimeDisc2;
931
932 // create solution vector & right hand side
934 SmartPtr<typename base_type::vector_type> b= u0->clone_without_values();
935
936 // create matrix operator
937 SmartPtr<typename base_type::assembled_operator_type> spAssOp=make_sp(new typename base_type::assembled_operator_type(tdisc_dep_type::m_spTimeDisc, gl));
938
939 // create additional solutions
942
943 // solution time series
944 SmartPtr<vector_time_series_type> m_spSolTimeSeries;
945 m_spSolTimeSeries=make_sp(new vector_time_series_type());
946 m_spSolTimeSeries->push(uold, t0);
947
948 SmartPtr<vector_time_series_type> m_spSolTimeSeries2;
949 m_spSolTimeSeries2=make_sp(new vector_time_series_type());
950 m_spSolTimeSeries2->push(u2old, t0);
951
952 // automatic selection of min/max time step
953 if (m_dtmin <= 0.0) { m_dtmin = (t1 - t0)/1.0e+5; }
954 if (m_dtmax <= 0.0) { m_dtmax = (t1 - t0)/10.0; }
955
956 // Aitken Neville extrapolation
957 const size_t tsteps[2] = {1,2};
958 std::vector<size_t> vsteps (tsteps, tsteps+2);
960
961 // integrate
962 if(!base_type::m_bNoLogOut)
963 UG_LOG("+++ Integrating: ["<< t0 <<", "<< t1 <<"]\n");
964
965 double t = t0;
966 int step = 0;
967
968 number dt = base_type::m_dt;
969 while((t < t1) && (t1-t > base_type::m_precisionBound))
970 {
971 // step: t -> t+dt
972 bool bSuccess = false;
973 while (!bSuccess){
974
975 // determine step size
976 if (dt<m_dtmin){
977 if(!base_type::m_bNoLogOut)
978 UG_LOG("Step size below minimum")
979 }
980 dt = std::min(dt, t1-t);
981
982 // basic step
983 if(!base_type::m_bNoLogOut)
984 UG_LOG("+++ Timestep: " << ++step << "\n");
985
986 tdisc.prepare_step(m_spSolTimeSeries, dt);
987 tdisc.assemble_linear(*spAssOp, *b, gl);
988 base_type::m_spLinearSolver->init(spAssOp, *u1);
989 base_type::m_spLinearSolver->apply(*u1, *b);
990
991
992 // control 1/2
993 if(!base_type::m_bNoLogOut)
994 UG_LOG("+++ Control: " << step << "\n");
995
996 tdisc2.prepare_step(m_spSolTimeSeries, 0.5*dt);
997 tdisc2.assemble_linear(*spAssOp, *b, gl);
998 base_type::m_spLinearSolver->init(spAssOp, *u2);
999 base_type::m_spLinearSolver->apply(*u2, *b);
1000
1001 // push back solution
1002 SmartPtr<typename base_type::vector_type> tmp2 = m_spSolTimeSeries2->oldest();
1003 (*tmp2)=*(u2.template cast_static<typename base_type::vector_type>());
1004 m_spSolTimeSeries2->push_discard_oldest(tmp2, t + 0.5*dt);
1005
1006 // control 2/2
1007 tdisc2.prepare_step(m_spSolTimeSeries2, 0.5*dt);
1008 tdisc2.assemble_linear(*spAssOp, *b, gl);
1009 base_type::m_spLinearSolver->init(spAssOp, *u2);
1010 base_type::m_spLinearSolver->apply(*u2, *b);
1011
1012 // obtain extrapolated solution
1013 timex.set_solution(u1, 0);
1014 timex.set_solution(u2, 1);
1015 timex.apply();
1016
1017 // predict (subsequent) time step
1018 number eps = timex.get_error_estimates()[0];
1019 number lambda = std::pow(0.8* m_tol/eps, 0.5);
1020
1021 number dtEst= dt*lambda;
1022 dtEst = std::min(dtEst, 1.5*dt);
1023 dtEst = std::min(dtEst, m_dtmax);
1024
1025 dt = dtEst;
1026 if (eps <= m_tol)
1027 {
1028 // ACCEPT STEP (and thus solution u2)
1029 if(!base_type::m_bNoLogOut)
1030 UG_LOG("ACCEPTING solution, dtnew=" << dt);
1031
1032 bSuccess = true;
1033 }
1034 else
1035 {
1036 // DISCARD step
1037 if(!base_type::m_bNoLogOut)
1038 UG_LOG("DISCARDING solution, dtnew=" << dt);
1039
1040 // => reset solutions
1041 VecAssign(*u1.template cast_dynamic<typename base_type::vector_type>(), *uold);
1042
1043 // => timeSeries2 has been updated...
1044 SmartPtr<typename base_type::vector_type> tmp = m_spSolTimeSeries2->oldest();
1045 VecAssign(*tmp, *uold);
1046 m_spSolTimeSeries2->push_discard_oldest(tmp, t);
1047
1048 }
1049
1050 }
1051
1052
1053 // prepare next loop
1054 t += dt;
1055
1056 // push updated solution into time series (and continue)
1057 SmartPtr<typename base_type::vector_type> tmp = m_spSolTimeSeries->oldest();
1058 VecAssign(*tmp, static_cast<typename base_type::vector_type> (*u2));
1059 m_spSolTimeSeries->push_discard_oldest(tmp, t);
1060
1061 }
1062
1063 return true;
1064};
1065
1067{
1068public:
1070 : m_dtMin(0.0), m_dtMax(std::numeric_limits<double>::max()), m_redFac(1.0), m_incFac(1.0)
1071 {}
1072
1073 void set_dt_min(double min) { m_dtMin = min; }
1074 double get_dt_min() { return m_dtMin; }
1075
1076 void set_dt_max(double max) { m_dtMax = max; }
1077 double get_dt_max() { return m_dtMax; }
1078
1079 void set_reduction_factor(double dec) { m_redFac = dec; }
1080 double get_reduction_factor() { return m_redFac; }
1081
1082 void set_increase_factor(double inc) { m_incFac = inc; }
1083 double get_increase_factor() { return m_incFac; }
1084
1085 void rescale(double alpha)
1086 { m_dtMin*= alpha; m_dtMax*= alpha;}
1087
1088protected:
1091};
1092
1094template<class TDomain, class TAlgebra>
1134
1135
1136} // ug
1137
1138#include "simple_integrator.hpp"
1139
1140namespace ug {
1141
1143template<class TDomain, class TAlgebra>
1145 public INonlinearTimeIntegrator<TDomain, TAlgebra>
1146{
1147public:
1148
1151
1154
1156 {
1157 int dstep = 0;
1158 auto tpoint = m_timePoints.begin();
1159 double tcurr = (tpoint == m_timePoints.end()) ? t1 : (*tpoint);
1160 double eps = 1e-8;
1161
1162 // Perform first step.
1163 //this->notify_init_step(u0, dstep, t0, tcurr-t0);
1164 bool status = m_wrappedIntegrator->apply(u1, tcurr*(1.0-eps), u0, t0);
1165 this->notify_finalize_step(u1, dstep++, tcurr, tcurr-t0);
1166
1167 // Repeat for all intermediate points.
1168 while (tpoint != m_timePoints.end())
1169 {
1170 tpoint++;
1171 double tnext = (tpoint == m_timePoints.end()) ? t1 : (*tpoint);
1172
1173 // Perform step.
1174 //this->notify_init_step(u1, dstep, tcurr, tnext-tcurr);
1175 status = status && m_wrappedIntegrator->apply(u1, tnext*(1.0-eps), u1, tcurr);
1176 this->notify_finalize_step(u1, dstep++, tnext, tnext-tcurr);
1177
1178 tcurr = tnext;
1179 }
1180 this->notify_end(u1, dstep, t1, 0.0);
1181 return status;
1182 }
1183
1184 void insert_points (std::vector<double> points)
1185 {
1186 m_timePoints = points;
1187 }
1188protected:
1190 std::vector<double> m_timePoints;
1191};
1192
1193} // namespace ug
1194
1195#endif /* TIME_INTEGRATOR_HPP_ */
Definition time_extrapolation.h:984
const std::vector< number > & get_error_estimates() const
Definition time_extrapolation.h:1028
void apply(size_t nstages, bool with_error=true)
best error estimate
Definition time_extrapolation.h:1063
void set_solution(SmartPtr< vector_type > soli, int i)
set solution (for stage i)
Definition time_extrapolation.h:1014
Sample class for integration observer: Output to VTK.
Definition time_integrator.hpp:204
ConnectionViewerOutputObserver(const char *filename)
Definition time_integrator.hpp:209
std::string m_filename
Definition time_integrator.hpp:228
ITimeIntegratorObserver< TDomain, TAlgebra > base_type
Definition time_integrator.hpp:206
virtual ~ConnectionViewerOutputObserver()
Definition time_integrator.hpp:215
number m_outputTime
Definition time_integrator.hpp:229
ConnectionViewerOutputObserver(const char *filename, number t_out)
Definition time_integrator.hpp:212
GridFunction< TDomain, TAlgebra > grid_function_type
Definition time_integrator.hpp:207
bool step_process(SmartPtr< grid_function_type > uNew, int step, number time, number dt) OVERRIDE
Definition time_integrator.hpp:218
Integrate over a given time interval (for a linear problem)
Definition time_integrator.hpp:755
ConstStepLinearTimeIntegrator(SmartPtr< time_disc_type > tDisc)
Definition time_integrator.hpp:773
base_type::linear_solver_type linear_solver_type
Definition time_integrator.hpp:761
base_type::grid_function_type grid_function_type
Definition time_integrator.hpp:762
void set_num_steps(int steps)
Definition time_integrator.hpp:779
VectorTimeSeries< typename base_type::vector_type > vector_time_series_type
Definition time_integrator.hpp:763
bool apply(SmartPtr< grid_function_type > u1, number t1, ConstSmartPtr< grid_function_type > u0, number t0)
Definition time_integrator.hpp:787
int m_numSteps
Definition time_integrator.hpp:769
ConstStepLinearTimeIntegrator(SmartPtr< time_disc_type > tDisc, SmartPtr< typename base_type::linear_solver_type > lSolver)
Definition time_integrator.hpp:776
ILinearTimeIntegrator< TDomain, TAlgebra > base_type
Definition time_integrator.hpp:759
ITimeDiscDependentObject< TAlgebra > tdisc_dep_type
Definition time_integrator.hpp:757
ITimeDiscretization< TAlgebra > time_disc_type
Definition time_integrator.hpp:760
This class integrates (t0, t1] with stops at intermediate points tk.
Definition time_integrator.hpp:1146
DiscontinuityIntegrator(SmartPtr< base_type > baseIntegrator)
Definition time_integrator.hpp:1152
void insert_points(std::vector< double > points)
Definition time_integrator.hpp:1184
SmartPtr< base_type > m_wrappedIntegrator
Definition time_integrator.hpp:1189
base_type::grid_function_type grid_function_type
Definition time_integrator.hpp:1150
INonlinearTimeIntegrator< TDomain, TAlgebra > base_type
Definition time_integrator.hpp:1149
bool apply(SmartPtr< grid_function_type > u1, number t1, ConstSmartPtr< grid_function_type > u0, number t0)
Definition time_integrator.hpp:1155
std::vector< double > m_timePoints
Definition time_integrator.hpp:1190
void assemble_linear(matrix_type &A, vector_type &b)
void assemble_rhs(vector_type &b)
integration of linear systems
Definition time_integrator.hpp:591
ILinearOperatorInverse< vector_type > linear_solver_type
Definition time_integrator.hpp:596
void set_linear_solver(SmartPtr< linear_solver_type > lSolver)
Definition time_integrator.hpp:608
AssembledLinearOperator< TAlgebra > assembled_operator_type
Definition time_integrator.hpp:597
ILinearTimeIntegrator()
Definition time_integrator.hpp:600
SmartPtr< linear_solver_type > m_spLinearSolver
Definition time_integrator.hpp:612
ITimeIntegrator< TDomain, TAlgebra > base_type
Definition time_integrator.hpp:594
ILinearTimeIntegrator(SmartPtr< linear_solver_type > lSolver)
Definition time_integrator.hpp:603
base_type::vector_type vector_type
Definition time_integrator.hpp:595
integration of non-linear systems (with bounds on dt)
Definition time_integrator.hpp:1097
double get_reduction_factor()
Definition time_integrator.hpp:1123
AssembledOperator< TAlgebra > assembled_operator_type
Definition time_integrator.hpp:1102
double get_dt_max()
Definition time_integrator.hpp:1120
void set_dt_min(double min)
Definition time_integrator.hpp:1116
base_type::vector_type vector_type
Definition time_integrator.hpp:1100
INonlinearTimeIntegrator()
Definition time_integrator.hpp:1104
ITimeIntegrator< TDomain, TAlgebra > base_type
Definition time_integrator.hpp:1099
ConstSmartPtr< solver_type > get_solver() const
Definition time_integrator.hpp:1110
TimeStepBounds m_dtBounds
Definition time_integrator.hpp:1130
IOperatorInverse< vector_type > solver_type
Definition time_integrator.hpp:1101
SmartPtr< solver_type > m_spSolver
Definition time_integrator.hpp:1129
double get_increase_factor()
Definition time_integrator.hpp:1126
void set_increase_factor(double inc)
Definition time_integrator.hpp:1125
void set_solver(SmartPtr< solver_type > solver)
Definition time_integrator.hpp:1107
void set_reduction_factor(double dec)
Definition time_integrator.hpp:1122
SmartPtr< solver_type > get_solver()
Definition time_integrator.hpp:1113
double get_dt_min()
Definition time_integrator.hpp:1117
void set_dt_max(double max)
Definition time_integrator.hpp:1119
ITimeDiscDependentObject.
Definition time_integrator.hpp:620
SmartPtr< time_disc_type > m_spTimeDisc
Definition time_integrator.hpp:630
SmartPtr< time_disc_type > get_time_disc()
Definition time_integrator.hpp:628
ITimeDiscretization< TAlgebra > time_disc_type
Definition time_integrator.hpp:622
ITimeDiscDependentObject(SmartPtr< time_disc_type > spTimeDisc)
Definition time_integrator.hpp:624
virtual void prepare_step(SmartPtr< VectorTimeSeries< vector_type > > prevSol, number dt)=0
Integrates over a given time interval [a,b] with step size dt.
Definition time_integrator.hpp:496
TDomain domain_type
Definition time_integrator.hpp:504
void set_precision_bound(double precisionBound)
Definition time_integrator.hpp:576
void prepare(grid_function_type &u)
prepares functions for application
Definition time_integrator.hpp:560
void init()
init operator
Definition time_integrator.hpp:551
double m_lower_tim
Definition time_integrator.hpp:509
virtual void init(grid_function_type const &u)
init operator depending on a function u
Definition time_integrator.hpp:538
GridFunction< TDomain, TAlgebra > grid_function_type
Definition time_integrator.hpp:505
void set_no_log_out(bool bNoLogOut)
Definition time_integrator.hpp:579
TAlgebra::matrix_type matrix_type
Definition time_integrator.hpp:501
double m_dt
Definition time_integrator.hpp:508
virtual ~ITimeIntegrator()
virtual destructor
Definition time_integrator.hpp:524
bool m_bNoLogOut
Definition time_integrator.hpp:514
void set_time_step(double dt)
Set initial time step.
Definition time_integrator.hpp:570
double m_upper_tim
Definition time_integrator.hpp:510
virtual bool apply(SmartPtr< grid_function_type > u1, number t1, ConstSmartPtr< grid_function_type > u0, number t0)=0
TAlgebra::vector_type vector_type
Definition time_integrator.hpp:500
double get_time_step()
Definition time_integrator.hpp:573
TAlgebra algebra_type
Definition time_integrator.hpp:499
void apply(grid_function_type &u1, const grid_function_type &u0)
Apply operator.
Definition time_integrator.hpp:564
double m_precisionBound
Definition time_integrator.hpp:512
ITimeIntegrator()
Definition time_integrator.hpp:519
Integration observer: Output using Lua callback.
Definition time_integrator.hpp:435
GridFunction< TDomain, TAlgebra > grid_function_type
Definition time_integrator.hpp:438
VTKOutput< TDomain::dim > vtk_type
Definition time_integrator.hpp:439
ITimeIntegratorObserver< TDomain, TAlgebra > base_type
Definition time_integrator.hpp:437
void add_integral_specs(const char *cmp, const char *subsets, int quadOrder, const char *idString)
Definition time_integrator.hpp:476
virtual ~IntegrationOutputObserver()
Definition time_integrator.hpp:456
bool step_process(SmartPtr< grid_function_type > uNew, int step, number time, number dt) OVERRIDE
Definition time_integrator.hpp:460
IntegrationOutputObserver()
Definition time_integrator.hpp:453
std::vector< IntegralSpecs > m_vIntegralData
Definition time_integrator.hpp:483
Integrate over a given time interval (for a linear problem)
Definition time_integrator.hpp:640
bool apply(SmartPtr< grid_function_type > u1, number t1, ConstSmartPtr< grid_function_type > u0, number t0)
Definition time_integrator.hpp:664
ITimeDiscretization< TAlgebra > time_disc_type
Definition time_integrator.hpp:647
ITimeDiscDependentObject< TAlgebra > tdisc_dep_type
Definition time_integrator.hpp:644
base_type::grid_function_type grid_function_type
Definition time_integrator.hpp:648
LinearTimeIntegrator(SmartPtr< time_disc_type > tDisc, SmartPtr< typename base_type::linear_solver_type > lSolver)
Definition time_integrator.hpp:655
VectorTimeSeries< typename base_type::vector_type > vector_time_series_type
Definition time_integrator.hpp:649
ILinearTimeIntegrator< TDomain, TAlgebra > base_type
Definition time_integrator.hpp:646
LinearTimeIntegrator(SmartPtr< time_disc_type > tDisc)
Definition time_integrator.hpp:652
Definition time_integrator.hpp:383
SmartPtr< UserData< number, grid_function_type::dim > > m_spReference
Definition time_integrator.hpp:424
PlotRefOutputObserver(SmartPtr< UserData< number, grid_function_type::dim > > spExactSol)
Definition time_integrator.hpp:389
virtual ~PlotRefOutputObserver()
Definition time_integrator.hpp:403
GridFunction< TDomain, TAlgebra > grid_function_type
Definition time_integrator.hpp:386
VTKOutput< TDomain::dim > vtk_type
Definition time_integrator.hpp:387
ITimeIntegratorObserver< TDomain, TAlgebra > base_type
Definition time_integrator.hpp:385
SmartPtr< vtk_type > m_sp_vtk
Definition time_integrator.hpp:425
bool step_process(SmartPtr< grid_function_type > uNew, int step, number time, number dt) OVERRIDE
Definition time_integrator.hpp:407
Integrate over a given time interval (for a linear problem)
Definition time_integrator.hpp:881
ITimeDiscDependentObject< TAlgebra > tdisc_dep_type
Definition time_integrator.hpp:883
VectorTimeSeries< typename base_type::vector_type > vector_time_series_type
Definition time_integrator.hpp:888
TimeIntegratorLinearAdaptive(SmartPtr< time_disc_type > tDisc1, SmartPtr< time_disc_type > tDisc2)
Definition time_integrator.hpp:901
double m_tol
Definition time_integrator.hpp:897
void set_time_step_max(number dt)
Definition time_integrator.hpp:912
ITimeDiscretization< TAlgebra > time_disc_type
Definition time_integrator.hpp:886
double m_dtmax
Definition time_integrator.hpp:897
bool apply(SmartPtr< grid_function_type > u1, number t1, ConstSmartPtr< grid_function_type > u0, number t0)
Definition time_integrator.hpp:925
void set_tol(double tol)
Definition time_integrator.hpp:910
void set_time_step_min(number dt)
Definition time_integrator.hpp:911
double m_dtmin
Definition time_integrator.hpp:897
SmartPtr< time_disc_type > m_spTimeDisc2
Definition time_integrator.hpp:895
ILinearTimeIntegrator< TDomain, TAlgebra > base_type
Definition time_integrator.hpp:885
base_type::grid_function_type grid_function_type
Definition time_integrator.hpp:887
Definition time_integrator.hpp:1067
void set_increase_factor(double inc)
Definition time_integrator.hpp:1082
double get_reduction_factor()
Definition time_integrator.hpp:1080
TimeStepBounds()
Definition time_integrator.hpp:1069
void set_dt_min(double min)
Definition time_integrator.hpp:1073
void set_dt_max(double max)
Definition time_integrator.hpp:1076
double m_incFac
Definition time_integrator.hpp:1090
double m_redFac
Definition time_integrator.hpp:1090
double m_dtMin
Definition time_integrator.hpp:1089
double get_dt_max()
Definition time_integrator.hpp:1077
void rescale(double alpha)
Definition time_integrator.hpp:1085
void set_reduction_factor(double dec)
Definition time_integrator.hpp:1079
double m_dtMax
Definition time_integrator.hpp:1089
double get_increase_factor()
Definition time_integrator.hpp:1083
double get_dt_min()
Definition time_integrator.hpp:1074
Sample class for integration observer: Output to VTK.
Definition time_integrator.hpp:81
bool step_process(SmartPtr< grid_function_type > uNew, int step, number time, number dt) OVERRIDE
Definition time_integrator.hpp:106
virtual ~VTKOutputObserver()
Definition time_integrator.hpp:96
SmartPtr< vtk_type > m_sp_vtk
Definition time_integrator.hpp:190
GridFunction< TDomain, TAlgebra > grid_function_type
Definition time_integrator.hpp:84
void set_output_scales(const std::vector< number > &vScales)
Definition time_integrator.hpp:100
SmartPtr< grid_function_type > m_uOld
Definition time_integrator.hpp:191
number m_startTime
Definition time_integrator.hpp:193
VTKOutputObserver(const char *filename, SmartPtr< vtk_type > vtk, number plotStep)
Definition time_integrator.hpp:93
VTKOutputObserver(const char *filename, SmartPtr< vtk_type > vtk)
Definition time_integrator.hpp:90
VTKOutputObserver()
Definition time_integrator.hpp:87
bool end_action(SmartPtr< grid_function_type > u, int step, number time, number dt) OVERRIDE
Definition time_integrator.hpp:166
std::vector< number > m_vOutputScales
Definition time_integrator.hpp:195
bool start_action(SmartPtr< grid_function_type > u, int step, number time, number dt) OVERRIDE
Definition time_integrator.hpp:112
number m_plotStep
Definition time_integrator.hpp:194
bool finalize_action(SmartPtr< grid_function_type > uNew, int step, number time, number dt) OVERRIDE
Definition time_integrator.hpp:126
std::string m_filename
Definition time_integrator.hpp:192
void writeToFile(SmartPtr< grid_function_type > u, int step, number time)
Definition time_integrator.hpp:177
ITimeIntegratorObserver< TDomain, TAlgebra > base_type
Definition time_integrator.hpp:83
VTKOutput< TDomain::dim > vtk_type
Definition time_integrator.hpp:85
void VecAssign(GPUVector< T > &dest, const GPUVector< T > &v1)
const NullSmartPtr SPNULL
#define UG_THROW(msg)
#define UG_LOG(msg)
#define UG_COND_THROW(cond, msg)
double number
#define LIMEX_PROFILE_FUNC()
Definition limex_tools.h:38
struct lua_State lua_State
lua_State * GetDefaultLuaState()
number L2Error(SmartPtr< TGridFunction > spGridFct1, const char *cmp1, SmartPtr< TGridFunction > spGridFct2, const char *cmp2, int quadOrder)
number Integral(number val, SmartPtr< TGridFunction > spGridFct)
void VecScaleAdd(double &dest, double alpha1, const double &v1, double alpha2, const double &v2)
SmartPtr< T, FreePolicy > make_sp(T *inst)
Definition time_integrator.hpp:442
int m_quadOrder
Definition time_integrator.hpp:448
std::string m_subsets
Definition time_integrator.hpp:447
std::string m_cmp
Definition time_integrator.hpp:446
std::string m_idString
Definition time_integrator.hpp:449
IntegralSpecs(const char *cmp, const char *subsets, int quadOrder, const char *idString)
Definition time_integrator.hpp:443
#define OVERRIDE
Definition time_integrator.hpp:39