ode_newton_bridge.h

Go to the documentation of this file.

 
#ifndef ODE_NEWTON_BRIDGE_H_
#define ODE_NEWTON_BRIDGE_H_
 
#include "ode_newton.h"
 
namespace ode {
 
template <class Integrator>
class OdeNewtonBridge : public OdeNewton {
 
    public:
 
        OdeNewtonBridge (unsigned long neq, unsigned long nnew, Integrator *integrator)
            : OdeNewton (nnew) {
 
            //make the default Newton iteration limit quite high
            set_iter_Newton(10000);
 
            //store pointers to the ODE solver object and its members
            integrator_ = integrator;
            sol_ = integrator_->sol_;
            Jac_ = integrator_->Jac_;
            dt_ = &(integrator_->dt_);
            //store ODE system size
            neq_ = integrator_->neq_;
            //store Newton system size;
            nnew_ = nnew;
 
            //temporary values for evaluation of ODE fun within Newton fun
            ftemp_ = new double[neq];
            //temporary solution values
            soltemp_ = new double[neq];
        }
        ~OdeNewtonBridge () {
            delete [] ftemp_;
            delete [] soltemp_;
        }
 
    protected:
 
        unsigned long neq_;
        unsigned long nnew_;
        Integrator *integrator_;
        double *sol_;
        double **Jac_;
        double *dt_;
 
        double *ftemp_;
        double *soltemp_;
 
 
        void fun(double *solin, double *fout) { integrator_->ode_fun_(solin, fout); }
 
        void jac(double *solin, double **Jout) { integrator_->ode_jac_(solin, Jout); }
};
 
 
template <class Integrator>
class OdeNewtonIRK : public OdeNewtonBridge<Integrator> {
 
    public:
 
        OdeNewtonIRK (unsigned long neq, int nk, Integrator *integrator) :
            OdeNewtonBridge<Integrator> (neq, neq*nk, integrator) {
 
            //number of stages or k vectors
            nk_ = nk;
            //pointer to stage k values
            k_ = integrator->k_;
            //pointers to tableau coefficients
            a = integrator->a;
            b = integrator->b;
        }
 
    protected:
 
        int nk_;
        double **k_;
        double **a;
        double *b;
};
 
 
template <class Integrator>
class OdeNewtonSDIRK : public OdeNewtonBridge<Integrator> {
 
    public:
 
        OdeNewtonSDIRK (unsigned long neq, Integrator *integrator) :
            OdeNewtonBridge<Integrator> (neq, neq, integrator) {
 
            //pointer to stage k values
            k_ = integrator->k_;
            //pointers to tableau coefficients
            gam = integrator->gam;
            a = integrator->a;
            b = integrator->b;
            //index of k vector being solved for
            ik_ = 0;
        }
 
        void set_ik (int ik) { ik_ = ik; }
 
    protected:
 
        double **k_;
        double gam;
        double **a;
        double *b;
        int ik_;
};
 
} // namespace ode
 
#endif