1#ifndef MULTISTEP_SOLVER_HH
2#define MULTISTEP_SOLVER_HH
13#include <dune/fem/solver/timeprovider.hh>
14#include <dune/fem/solver/rungekutta/explicit.hh>
30 template<
class Operator>
34 typedef typename Operator::DestinationType DestinationType;
35 typedef typename DestinationType :: DiscreteFunctionSpaceType SpaceType;
38 std::vector< std::vector<double> > a;
39 std::vector<double> b;
40 std::vector<double> c;
43 std::vector<DestinationType*> Uj;
44 std::vector<DestinationType*> Fj;
45 std::vector<double> deltat_;
58 int pord,
bool verbose =
true ) :
73 for (
int i=0; i<ord_; ++i)
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;
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;
97 a[0][0]=0.; a[0][1]=0.;
98 a[1][0]=1.0; a[1][1]=0.;
107 default : std::cerr <<
"Runge-Kutta method of this order not implemented"
112 for (
size_t i=0; i<steps_; ++i)
114 Fj.push_back(
new DestinationType(
"FMS",op_.space()) );
116 Uj.push_back(
new DestinationType(
"UMS",op_.space()) );
122 for(
size_t i=0; i<Fj.size(); ++i) {
125 for(
size_t i=0; i<Uj.size(); ++i) {
146 Uj[Uj.size()-1]->assign(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]));
162 deltat_[steps_-1] = tp_.
deltaT();
166 Uj[steps_-1]->assign(U0);
169 op_.setTime( tp_.
time() );
171 op_(*(Uj[steps_-1]), *(Fj[steps_-1]));
177 double alpha[steps_];
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] <<
" "
191 << beta[0] <<
" " << beta[1] <<
" "
192 << deltat_[0] <<
" " << deltat_[1]
204 double w1 = deltat_[1]/deltat_[0];
205 double w2 = deltat_[2]/deltat_[1];
206 double g = w1*w2/(1+w1);
213 beta[0] = 0.*deltat_[steps_-1];
214 beta[1] = 0.*deltat_[steps_-1];
215 beta[2] = (1+g)*deltat_[steps_-1];
235 beta[0] = 4./9.*deltat_[steps_-1];
236 beta[1] = 0.*deltat_[steps_-1];
237 beta[2] = 0.*deltat_[steps_-1];
238 beta[3] = 16./9.*deltat_[steps_-1];
262 for (
int j=0;j<=steps_;j++) {
264 for (
int i=0;i<j;i++)
269 for (
int q=0;q<p;q++)
271 for (
int j=1;j<steps_;j++) {
273 for (
int q=0;q<p-1;q++)
275 cond += (alpha[j]*M[j]+p*beta[j])*Mp;
279 std::cout <<
"# MS Cond for p= " << 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] <<
" ";
291 std::cout << std::endl;
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]));
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];
315 DestinationType Uval(
"Utmp",op_.space());
318 const double dt = tp_.
deltaT();
320 const double t = tp_.
time();
331 for (
int i=1; i<ord_; ++i)
334 for (
int j=0; j<i ; ++j)
336 Uval.axpy((a[i][j]*dt), *(Fj[j]));
340 op_.setTime( t + c[i]*dt );
350 for (
int j=0; j<ord_; ++j)
352 U0.axpy((b[j]*dt), *(Fj[j]));
368 template<
class Operator>
369 class ExplMultiStep :
public TimeProvider ,
370 public ExplMultiStepBase<Operator>
372 typedef ExplMultiStepBase<Operator> BaseType;
374 typedef typename Operator :: DestinationType DestinationType;
375 typedef typename DestinationType :: DiscreteFunctionSpaceType
SpaceType;
376 typedef typename SpaceType :: GridType :: Traits :: Communication DuneCommunicatorType;
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)
391 op.timeProvider(
this);
394 void initialize(
const DestinationType& U0)
396 if(! this->initialized_)
399 BaseType :: initialize(U0);
402 this->tp_.syncTimeStep();
406 double solve(
typename Operator::DestinationType& U0)
411 BaseType :: solve (U0);
414 this->tp_.augmentTime();
417 this->tp_.syncTimeStep();
419 return this->tp_.time();
422 void printGrid(
int nr,
const typename Operator::DestinationType& U)
424 if (time()>=savetime_) {
425 printSGrid(time(),savestep_*10+nr,this->op_.space(),U);
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";
438 this->op_.printmyInfo(filename);
443 ParallelTimeProvider<DuneCommunicatorType> tp_;
454 template<
class DestinationImp>
459 typedef DestinationImp DestinationType;
477 cfl_(
std ::
min( 0.45 / (2.0 * pord+1) / double(pord), 1.0 ) )
480 std::cout <<
"ExplicitMultiStepSolver: cfl = " << cfl_ <<
"!" << std :: endl;
496 if( ! this->initialized_ )
Definition: multistep.hh:32
~ExplMultiStepBase()
destructor
Definition: multistep.hh:120
ExplMultiStepBase(Operator &op, TimeProviderBase &tp, int pord, bool verbose=true)
constructor
Definition: multistep.hh:57
void solveRK(DestinationType &U0)
solve the system
Definition: multistep.hh:313
void initialize(const DestinationType &U0)
apply operator once to get dt estimate
Definition: multistep.hh:131
void solve(DestinationType &U0)
solve the system
Definition: multistep.hh:142
Exlicit multi step ODE solver.
Definition: multistep.hh:458
void solve(DestinationType &U0)
solve system
Definition: multistep.hh:493
ExplicitMultiStepSolver(OperatorType &op, TimeProviderBase &tp, int pord, bool verbose=false)
constructor
Definition: multistep.hh:474
void initialize(const DestinationType &U0)
apply operator once to get dt estimate
Definition: multistep.hh:487
virtual ~ExplicitMultiStepSolver()
destructor
Definition: multistep.hh:484
Interface class for ODE Solver.
Definition: odesolverinterface.hh:21
A vector valued function space.
Definition: functionspace.hh:60
general base for time providers
Definition: timeprovider.hh:36
double time() const
obtain the current time
Definition: timeprovider.hh:94
void provideTimeStepEstimate(const double dtEstimate)
set time step estimate to minimum of given value and internal time step estiamte
Definition: timeprovider.hh:142
double deltaT() const
obtain the size of the current time step
Definition: timeprovider.hh:113
Default exception if a function was called while the object is not in a valid state for that function...
Definition: exceptions.hh:281
An abstract operator Interface class for Operators. Operators are applied to Functions and the result...
Definition: operator.hh:238
#define DUNE_THROW(E, m)
Definition: exceptions.hh:218
constexpr auto min
Function object that returns the smaller of the given values.
Definition: hybridutilities.hh:506
V cond(ADLTag< 0 >, const Mask< V > &mask, const V &ifTrue, const V &ifFalse)=delete
implements Simd::cond()
Dune namespace.
Definition: alignedallocator.hh:13
void printGrid(const GridType &grid, const Dune::MPIHelper &helper, std::string output_file="printgrid", 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:73
abstract operator
Definition: operator.hh:34