1#ifndef MULTISTEP_SOLVER_HH
2#define MULTISTEP_SOLVER_HH
30 template<
class Operator>
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;
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)
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);
152 for (
size_t i=0;i<
steps_-1;i++) {
181 std::cerr <<
"the two step scheme of order 2 has negative beta coeffs and requires adjoint space operator!" << std::endl;
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;
188 std::cout <<
"# MS Coeffs : "
189 << alpha[0] <<
" " << alpha[1] <<
" "
191 << beta[0] <<
" " << beta[1] <<
" "
206 double g = w1*w2/(1+w1);
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]*
291 std::cout << std::endl;
297 for (
size_t i=0;i<
steps_-1;i++) {
298 U0.axpy(alpha[i], *(
Uj[i]));
299 U0.axpy(beta[i], *(
Fj[i]));
304 for (
size_t i=1;i<
steps_;i++) {
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>
374 typedef typename Operator :: DestinationType DestinationType;
375 typedef typename DestinationType :: DiscreteFunctionSpaceType SpaceType;
376 typedef typename SpaceType :: GridType :: Traits :: CollectiveCommunication DuneCommunicatorType;
385 ExplMultiStep (
Operator& op,
int pord,
double cfl,
bool verbose =
true ) :
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;
489 BaseType :: initialize(U0);
498 DUNE_THROW(InvalidStateException,
"ExplicitMultiStepSolver wasn't initialized before first call!");
502 BaseType :: solve(U0);
Definition: bindguard.hh:11
static double min(const Double &v, const double p)
Definition: double.hh:386
Definition: multistep.hh:17
abstract operator
Definition: operator.hh:34
interface for time evolution operators
Definition: spaceoperatorif.hh:38
Definition: multistep.hh:32
std::vector< double > deltat_
Definition: multistep.hh:45
~ExplMultiStepBase()
destructor
Definition: multistep.hh:120
std::vector< double > b
Definition: multistep.hh:39
Operator::DestinationType DestinationType
Definition: multistep.hh:34
ExplMultiStepBase(Operator &op, TimeProviderBase &tp, int pord, bool verbose=true)
constructor
Definition: multistep.hh:57
size_t steps_
Definition: multistep.hh:41
void solveRK(DestinationType &U0)
solve the system
Definition: multistep.hh:313
std::vector< DestinationType * > Uj
Definition: multistep.hh:43
void initialize(const DestinationType &U0)
apply operator once to get dt estimate
Definition: multistep.hh:131
DestinationType::DiscreteFunctionSpaceType SpaceType
Definition: multistep.hh:35
bool initialized_
Definition: multistep.hh:362
std::vector< std::vector< double > > a
Definition: multistep.hh:38
const Operator & op_
Definition: multistep.hh:358
bool msFirst
Definition: multistep.hh:46
double gamma_
Definition: multistep.hh:42
bool msInit
Definition: multistep.hh:46
void solve(DestinationType &U0)
solve the system
Definition: multistep.hh:142
TimeProviderBase & tp_
Definition: multistep.hh:360
std::vector< DestinationType * > Fj
Definition: multistep.hh:44
std::vector< double > c
Definition: multistep.hh:40
const int ord_
Definition: multistep.hh:48
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
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
manager for global simulation time of time-dependent solutions
Definition: timeprovider.hh:405