DUNE-FEM (unstable)

multistep.hh
1#ifndef MULTISTEP_SOLVER_HH
2#define MULTISTEP_SOLVER_HH
3
4// inlcude all used headers before, that they don not appear in DuneODE
5
6//- system includes
7#include <iostream>
8#include <cmath>
9#include <vector>
10#include <cassert>
11
12//- Dune includes
13#include <dune/fem/solver/timeprovider.hh>
14#include <dune/fem/solver/rungekutta/explicit.hh>
15
16namespace DuneODE
17{
18
19 using namespace Dune;
20 using namespace Fem;
21 using namespace std;
22
30 template<class Operator>
31 class ExplMultiStepBase
32 {
33 public:
34 typedef typename Operator::DestinationType DestinationType;
35 typedef typename DestinationType :: DiscreteFunctionSpaceType SpaceType;
36 // typedef typename SpaceType :: GridType :: Traits :: Communication DuneCommunicatorType;
37 protected:
38 std::vector< std::vector<double> > a;
39 std::vector<double> b;
40 std::vector<double> c;
41 size_t steps_;
42 double gamma_;
43 std::vector<DestinationType*> Uj;
44 std::vector<DestinationType*> Fj;
45 std::vector<double> deltat_;
46 bool msInit,msFirst;
47 protected:
48 const int ord_;
49
50 public:
57 ExplMultiStepBase(Operator& op, TimeProviderBase& tp,
58 int pord, bool verbose = true ) :
59 a(0),b(0),c(0)
60 , steps_(pord+1)
61 , gamma_(0.5) // 1./3. -> Shu
62 , Uj(0)
63 , Fj(0)
64 , deltat_(steps_)
65 , msInit(false)
66 , msFirst(true)
67 , ord_(pord)
68 , op_(op)
69 , tp_(tp)
70 , initialized_(false)
71 {
72 a.resize(ord_);
73 for (int i=0; i<ord_; ++i)
74 {
75 a[i].resize(ord_);
76 }
77 b.resize(ord_);
78 c.resize(ord_);
79
80 switch (ord_) {
81 case 4 :
82 a[0][0]=0.; a[0][1]=0.; a[0][2]=0.; a[0][3]=0.;
83 a[1][0]=1.0; a[1][1]=0.; a[1][2]=0.; a[1][3]=0.;
84 a[2][0]=0.25; a[2][1]=0.25; a[2][2]=0.; a[2][3]=0.;
85 a[3][0]=1./6.; a[3][1]=1./6.; a[3][2]=2./3.; a[3][3]=0.;
86 b[0]=1./6.; b[1]=1./6.; b[2]=2./3.; b[3]=0.;
87 c[0]=0.; c[1]=1.0; c[2]=0.5; c[3]=1.0;
88 break;
89 case 3 :
90 a[0][0]=0.; a[0][1]=0.; a[0][2]=0.;
91 a[1][0]=1.0; a[1][1]=0.; a[1][2]=0.;
92 a[2][0]=0.25; a[2][1]=0.25; a[2][2]=0.;
93 b[0]=1./6.; b[1]=1./6.; b[2]=2./3.;
94 c[0]=0.; c[1]=1; c[2]=0.5;
95 break;
96 case 2 :
97 a[0][0]=0.; a[0][1]=0.;
98 a[1][0]=1.0; a[1][1]=0.;
99 b[0]=0.5; b[1]=0.5;
100 c[0]=0; c[1]=1;
101 break;
102 case 1:
103 a[0][0]=0.;
104 b[0]=1.;
105 c[0]=0.;
106 break;
107 default : std::cerr << "Runge-Kutta method of this order not implemented"
108 << std::endl;
109 abort();
110 }
111 // create memory for intermediate stages
112 for (size_t i=0; i<steps_; ++i)
113 {
114 Fj.push_back(new DestinationType("FMS",op_.space()) );
115 }
116 Uj.push_back(new DestinationType("UMS",op_.space()) );
117 }
118
120 ~ExplMultiStepBase()
121 {
122 for(size_t i=0; i<Fj.size(); ++i) {
123 delete Fj[i];
124 }
125 for(size_t i=0; i<Uj.size(); ++i) {
126 delete Uj[i];
127 }
128 }
129
131 void initialize(const DestinationType& U0)
132 {
133 if( ! initialized_ )
134 {
135 // Compute Steps
136 op_(U0, (*Fj[0]));
137 initialized_ = true;
138 }
139 }
140
142 void solve(DestinationType& U0)
143 {
144 if (!msInit) {
145 // Startup phase using RungeKutta
146 Uj[Uj.size()-1]->assign(U0);
147 solveRK(U0);
148 deltat_[Uj.size()-1] = tp_.deltaT();
149 Uj.push_back(new DestinationType("UMS",op_.space()) );
150 if (Uj.size()==steps_) {
151 Uj[steps_-1]->assign(U0);
152 for (size_t i=0;i<steps_-1;i++) {
153 op_(*(Uj[i]),(*Fj[i]));
154 }
155 msInit=true;
156 }
157 return;
158 }
159 // now multistep
160
161 // get cfl * timeStepEstimate
162 deltat_[steps_-1] = tp_.deltaT();
163 // get time
164
165 // Compute Steps
166 Uj[steps_-1]->assign(U0);
167
168 // set new time
169 op_.setTime( tp_.time() );
170
171 op_(*(Uj[steps_-1]), *(Fj[steps_-1]));
172
173 // provide time step estimate
174 tp_.provideTimeStepEstimate( op_.timeStepEstimate() );
175
176 // Perform Update
177 double alpha[steps_];
178 double beta[steps_];
179 switch (steps_) {
180 case 2: {
181 std::cerr << "the two step scheme of order 2 has negative beta coeffs and requires adjoint space operator!" << std::endl;
182 double w = deltat_[steps_-1]/deltat_[steps_-2];
183 double alpha2 = (1.+2.*gamma_*w)/(1.+w);
184 alpha[1] = -((1.-2.*gamma_)*w-1.)/alpha2;
185 alpha[0] = -((2.*gamma_-1.)*w*w/(1.+w))/alpha2;
186 beta[1] = (1.+gamma_*w)/alpha2*deltat_[steps_-1];
187 beta[0] = -(gamma_*w)/alpha2*deltat_[steps_-1];
188 std::cout << "# MS Coeffs : "
189 << alpha[0] << " " << alpha[1] << " "
190 << alpha2 << " "
191 << beta[0] << " " << beta[1] << " "
192 << deltat_[0] << " " << deltat_[1]
193 << std::endl;
194 } break;
195 case 3: {
196 /* Shu
197 alpha[0] = 1./4.;
198 alpha[1] = 0.;
199 alpha[2] = 3./4.;
200 beta[0] = 0.*deltat_[steps_-1];
201 beta[1] = 0.*deltat_[steps_-1];
202 beta[2] = 3./2.*deltat_[steps_-1];
203 */
204 double w1 = deltat_[1]/deltat_[0];
205 double w2 = deltat_[2]/deltat_[1];
206 double g = w1*w2/(1+w1);
207 // CFL < 1-g -> TVD
208 // therefore need g<1!
209 assert(g<1);
210 alpha[0] = g*g;
211 alpha[1] = 0.;
212 alpha[2] = 1.-g*g; // > 0 due to CFL
213 beta[0] = 0.*deltat_[steps_-1];
214 beta[1] = 0.*deltat_[steps_-1];
215 beta[2] = (1+g)*deltat_[steps_-1]; // always > 0
216 /*
217 double w1 = deltat_[1]/deltat_[0];
218 double w2 = deltat_[2]/deltat_[1];
219 double w1_1=1.+w1;
220 alpha[0] = (w1*w2-gamma_*w1*w1)/(w1_1*w1_1);
221 alpha[1] = gamma_;
222 alpha[2] = 1.-alpha[1]-alpha[0];
223 beta[0] = 0.*deltat_[steps_-1];
224 beta[1] = 0.*deltat_[steps_-1];
225 beta[2] = (1+w1+w1*w2-gamma_/w2*(1.+2.*w1))/w1_1*
226 deltat_[steps_-1]; // always > 0
227 */
228 } break;
229 case 4: {
230
231 alpha[0] = 11./27.; // 1./9.;
232 alpha[1] = 0.;
233 alpha[2] = 0.;
234 alpha[3] = 16./27.; // 8./9.;
235 beta[0] = 4./9.*deltat_[steps_-1]; // 4./9.
236 beta[1] = 0.*deltat_[steps_-1];
237 beta[2] = 0.*deltat_[steps_-1];
238 beta[3] = 16./9.*deltat_[steps_-1]; // 4./3.
239 /*
240 double k3 = deltat_[3];
241 double K2 = deltat_[0]+deltat_[1]+deltat_[2];
242 double K3 = K2 + k3;
243 double q1 = K3/K2;
244 double q2 = k3/K2;
245 double g = q2*q2*(3.+2.*q2);
246 alpha[0] = g;
247 alpha[1] = 0.;
248 alpha[2] = 0.;
249 alpha[3] = 1.-g;
250 beta[0] = q2*q1*deltat_[steps_-1]; // 0.
251 beta[1] = 0.*deltat_[steps_-1];
252 beta[2] = 0.*deltat_[steps_-1];
253 beta[3] = q1*q1*deltat_[steps_-1]; // 4./3.
254 */
255 } break;
256 }
257 #if 0 // !NDEBUG
258 { // Check order conditions
259 int p = 1;
260 double cond;
261 double M[steps_+1];
262 for (int j=0;j<=steps_;j++) {
263 M[j] = 0;
264 for (int i=0;i<j;i++)
265 M[j] += deltat_[i];
266 }
267 do {
268 cond = -1.;
269 for (int q=0;q<p;q++)
270 cond *= M[steps_];
271 for (int j=1;j<steps_;j++) {
272 double Mp=1;
273 for (int q=0;q<p-1;q++)
274 Mp *= M[j];
275 cond += (alpha[j]*M[j]+p*beta[j])*Mp;
276 }
277 if (p==1)
278 cond += beta[0];
279 std::cout << "# MS Cond for p= " << p
280 << " " << cond
281 << std::endl;
282 ++p;
283 } while (fabs(cond)<1e-7 && p<10);
284 std::cout << " # CFL and stability: ";
285 for (int i=0;i<steps_;i++) {
286 if (fabs(beta[i])>1e-10) {
287 std::cout << alpha[i]/beta[i]*
288 deltat_[steps_-1] << " ";
289 }
290 }
291 std::cout << std::endl;
292 }
293 #endif
294 // Update
295 U0 *= alpha[steps_-1];
296 U0.axpy(beta[steps_-1], *(Fj[steps_-1]) );
297 for (size_t i=0;i<steps_-1;i++) {
298 U0.axpy(alpha[i], *(Uj[i]));
299 U0.axpy(beta[i], *(Fj[i]));
300 }
301
302 DestinationType* Utmp = Uj[0];
303 DestinationType* Ftmp = Fj[0];
304 for (size_t i=1;i<steps_;i++) {
305 deltat_[i-1] = deltat_[i];
306 Uj[i-1] = Uj[i];
307 Fj[i-1] = Fj[i];
308 }
309 Uj[steps_-1] = Utmp;
310 Fj[steps_-1] = Ftmp;
311 }
313 void solveRK(DestinationType& U0)
314 {
315 DestinationType Uval("Utmp",op_.space());
316
317 // get cfl * timeStepEstimate
318 const double dt = tp_.deltaT();
319 // get time
320 const double t = tp_.time();
321
322 // set new time
323 op_.setTime( t );
324
325 // Compute Steps
326 op_(U0, *(Fj[0]));
327
328 // provide time step estimate
329 tp_.provideTimeStepEstimate( op_.timeStepEstimate() );
330
331 for (int i=1; i<ord_; ++i)
332 {
333 Uval.assign(U0);
334 for (int j=0; j<i ; ++j)
335 {
336 Uval.axpy((a[i][j]*dt), *(Fj[j]));
337 }
338
339 // set new time
340 op_.setTime( t + c[i]*dt );
341
342 // apply operator
343 op_(Uval,*(Fj[i]));
344
345 // provide time step estimate
346 tp_.provideTimeStepEstimate( op_.timeStepEstimate() );
347 }
348
349 // Perform Update
350 for (int j=0; j<ord_; ++j)
351 {
352 U0.axpy((b[j]*dt), *(Fj[j]));
353 }
354 }
355
356 protected:
357 // operator to solve for
358 const Operator& op_;
359 // time provider
360 TimeProviderBase& tp_;
361 // init flag
362 bool initialized_;
363 };
364
365#if 0
368 template<class Operator>
369 class ExplMultiStep : public TimeProvider ,
370 public ExplMultiStepBase<Operator>
371 {
372 typedef ExplMultiStepBase<Operator> BaseType;
373 public:
374 typedef typename Operator :: DestinationType DestinationType;
375 typedef typename DestinationType :: DiscreteFunctionSpaceType SpaceType;
376 typedef typename SpaceType :: GridType :: Traits :: Communication DuneCommunicatorType;
377
378 public:
385 ExplMultiStep (Operator& op,int pord,double cfl, bool verbose = true ) :
386 TimeProvider(0.0,(cfl/double(pord))*0.95),
387 BaseType(op,*this,pord,verbose),
388 tp_(op.space().gridPart().comm(),*this),
389 savetime_(0.0), savestep_(1)
390 {
391 op.timeProvider(this);
392 }
393
394 void initialize(const DestinationType& U0)
395 {
396 if(! this->initialized_)
397 {
398 // initialize
399 BaseType :: initialize(U0);
400
401 // global min of dt and reset of dtEstimate
402 this->tp_.syncTimeStep();
403 }
404 }
405
406 double solve(typename Operator::DestinationType& U0)
407 {
408 initialize( U0 );
409
410 // solve ode
411 BaseType :: solve (U0);
412
413 // calls setTime ( t + dt );
414 this->tp_.augmentTime();
415
416 // global min of dt and reset of dtEstimate
417 this->tp_.syncTimeStep();
418
419 return this->tp_.time();
420 }
421
422 void printGrid(int nr, const typename Operator::DestinationType& U)
423 {
424 if (time()>=savetime_) {
425 printSGrid(time(),savestep_*10+nr,this->op_.space(),U);
426 ++savestep_;
427 savetime_+=0.001;
428 }
429 }
430
431 void printmyInfo(string filename) const {
432 std::ostringstream filestream;
433 filestream << filename;
434 std::ofstream ofs(filestream.str().c_str(), std::ios::app);
435 ofs << "ExplMultiStep, steps: " << this->ord_ << "\n\n";
436 ofs << " cfl: " << this->tp_.cfl() << "\\\\\n\n";
437 ofs.close();
438 this->op_.printmyInfo(filename);
439 }
440
441 private:
442 // TimeProvider with communicator
443 ParallelTimeProvider<DuneCommunicatorType> tp_;
444 double savetime_;
445 int savestep_;
446 };
447#endif
448
454 template<class DestinationImp>
455 class ExplicitMultiStepSolver :
456 public OdeSolverInterface<DestinationImp> ,
457 public ExplMultiStepBase<SpaceOperatorInterface<DestinationImp> >
458 {
459 typedef DestinationImp DestinationType;
460 typedef SpaceOperatorInterface<DestinationImp> OperatorType;
461 typedef ExplMultiStepBase<OperatorType> BaseType;
462
463 private:
464 TimeProviderBase &timeProvider_;
465 double cfl_;
466
467 public:
474 ExplicitMultiStepSolver(OperatorType& op, TimeProviderBase& tp, int pord, bool verbose = false) :
475 BaseType(op,tp,pord,verbose),
476 timeProvider_(tp),
477 cfl_( std :: min( 0.45 / (2.0 * pord+1) / double(pord), 1.0 ) )
478 {
479 if(verbose)
480 std::cout << "ExplicitMultiStepSolver: cfl = " << cfl_ << "!" << std :: endl;
481 }
482
484 virtual ~ExplicitMultiStepSolver() {}
485
487 void initialize(const DestinationType& U0)
488 {
489 BaseType :: initialize(U0);
490 }
491
493 void solve(DestinationType& U0)
494 {
495 // initialize
496 if( ! this->initialized_ )
497 {
498 DUNE_THROW(InvalidStateException,"ExplicitMultiStepSolver wasn't initialized before first call!");
499 }
500
501 // solve ode
502 BaseType :: solve(U0);
503 }
504 };
505
507} // namespace DuneODE
508
509#endif // #ifndef MULTISTEP_SOLVER_HH
DuneODE::ExplMultiStepBase
Definition: multistep.hh:32
DuneODE::ExplMultiStepBase::~ExplMultiStepBase
~ExplMultiStepBase()
destructor
Definition: multistep.hh:120
DuneODE::ExplMultiStepBase::ExplMultiStepBase
ExplMultiStepBase(Operator &op, TimeProviderBase &tp, int pord, bool verbose=true)
constructor
Definition: multistep.hh:57
DuneODE::ExplMultiStepBase::solveRK
void solveRK(DestinationType &U0)
solve the system
Definition: multistep.hh:313
DuneODE::ExplMultiStepBase::initialize
void initialize(const DestinationType &U0)
apply operator once to get dt estimate
Definition: multistep.hh:131
DuneODE::ExplMultiStepBase::solve
void solve(DestinationType &U0)
solve the system
Definition: multistep.hh:142
DuneODE::ExplicitMultiStepSolver
Exlicit multi step ODE solver.
Definition: multistep.hh:458
DuneODE::ExplicitMultiStepSolver::solve
void solve(DestinationType &U0)
solve system
Definition: multistep.hh:493
DuneODE::ExplicitMultiStepSolver::ExplicitMultiStepSolver
ExplicitMultiStepSolver(OperatorType &op, TimeProviderBase &tp, int pord, bool verbose=false)
constructor
Definition: multistep.hh:474
DuneODE::ExplicitMultiStepSolver::initialize
void initialize(const DestinationType &U0)
apply operator once to get dt estimate
Definition: multistep.hh:487
DuneODE::ExplicitMultiStepSolver::~ExplicitMultiStepSolver
virtual ~ExplicitMultiStepSolver()
destructor
Definition: multistep.hh:484
DuneODE::OdeSolverInterface
Interface class for ODE Solver.
Definition: odesolverinterface.hh:21
Dune::Fem::FunctionSpace
A vector valued function space.
Definition: functionspace.hh:60
Dune::Fem::TimeProviderBase
general base for time providers
Definition: timeprovider.hh:36
Dune::Fem::TimeProviderBase::time
double time() const
obtain the current time
Definition: timeprovider.hh:94
Dune::Fem::TimeProviderBase::provideTimeStepEstimate
void provideTimeStepEstimate(const double dtEstimate)
set time step estimate to minimum of given value and internal time step estiamte
Definition: timeprovider.hh:142
Dune::Fem::TimeProviderBase::deltaT
double deltaT() const
obtain the size of the current time step
Definition: timeprovider.hh:113
Dune::InvalidStateException
Default exception if a function was called while the object is not in a valid state for that function...
Definition: exceptions.hh:375
Dune::Operator
An abstract operator Interface class for Operators. Operators are applied to Functions and the result...
Definition: operator.hh:238
DUNE_THROW
#define DUNE_THROW(E,...)
Definition: exceptions.hh:314
Dune::Hybrid::min
constexpr auto min
Function object that returns the smaller of the given values.
Definition: hybridutilities.hh:507
Dune::Simd::Overloads::cond
V cond(ADLTag< 0 >, const Mask< V > &mask, const V &ifTrue, const V &ifFalse)=delete
implements Simd::cond()
Dune
Dune namespace.
Definition: alignedallocator.hh:13
Dune::printGrid
void printGrid(const GridType &grid, std::string output_file, int size=2000, bool execute_plot=true, bool png=true, bool local_corner_indices=true, bool local_intersection_indices=true, bool outer_normals=true)
Print a grid as a gnuplot for testing and development.
Definition: printgrid.hh:72
std
STL namespace.
Dune::Fem::Operator
abstract operator
Definition: operator.hh:34
Creative Commons License   |  Legal Statements / Impressum  |  Hosted by TU Dresden & Uni Heidelberg  |  generated with Hugo v0.111.3 (Nov 3, 23:36, 2025)