timestepper.hpp

Go to the documentation of this file.

#include <petsc/private/tsimpl.h>        /*I "petscts.h"  I*/
#include <petscts.h>
#include <petscksp.h>
#include "mastereq.hpp"
#include <assert.h> 
#include <iostream> 
#include "defs.hpp"
#include "output.hpp"
#include "optimtarget.hpp"
#include <deque>
#pragma once
 
class TimeStepper{
  protected:
    Vec x; 
    Vec xadj; 
    Vec xprimal; 
    std::vector<Vec> store_states; 
    std::vector<Vec> dpdm_states; 
    int mpirank_world; 
    int mpisize_petsc; 
    int mpirank_petsc; 
    PetscInt localsize_u; 
    PetscInt ilow; 
    PetscInt iupp; 
 
    bool eval_leakage; 
    bool eval_energy; 
    bool eval_dpdm; 
    bool eval_weightedcost; 
 
    double leakage_integral; 
    double energy_integral; 
    double dpdm_integral; 
    double weightedcost_integral; 
    double weightedcost_width; 
 
  public:
    MasterEq* mastereq; 
    int ntime; 
    double total_time; 
    double dt; 
    bool writeTrajectoryDataFiles;  
 
    Vec redgrad; 
 
    OptimTarget* optim_target; 
    Output* output; 
 
  public: 
    bool storeFWD; 
 
    TimeStepper(); 
 
    TimeStepper(MasterEq* mastereq_, int ntime_, double total_time_, Output* output_, bool storeFWD_); 
 
    virtual ~TimeStepper(); 
 
    Vec getState(size_t tindex);
 
    void setEvalLeakage(bool flag){ eval_leakage = flag; };
    void setEvalWeightedCost(bool flag, double width){ eval_weightedcost = flag; weightedcost_width = width; };
    void setEvalDPDM(bool flag){ eval_dpdm = flag; };
    void setEvalEnergy(bool flag){ eval_energy = flag; };
 
    double getLeakageIntegral(){ return leakage_integral; };
    double getWeightedCostIntegral(){ return weightedcost_integral; };
    double getEnergyIntegral(){ return energy_integral; };
    double getDPDMIntegral(){ return dpdm_integral; };
 
    Vec solveODE(int initid, Vec rho_t0);
 
    void solveAdjointODE(Vec rho_t0_bar, Vec finalstate, double Jbar_leakage, double Jbar_weightedcost, double Jbar_dpdm, double Jbar_energy);
 
    double evalLeakage(const Vec x);
 
    void evalLeakage_diff(const Vec x, Vec xbar, double Jbar);
 
    double evalWeightedCost(double time, const Vec x);
 
    void evalWeightedCost_diff(double time, const Vec x, Vec xbar, double Jbar);
 
    double evalDpDm(Vec x, Vec xm1, Vec xm2);
 
    void evalDpDm_diff(int n, Vec xbar, double Jbar);
    
    double evalEnergy(double time);
 
    void evalEnergy_diff(double time, double Jbar, Vec redgrad);
 
    virtual void evolveFWD(const double tstart, const double tstop, Vec x) = 0;
 
    virtual void evolveBWD(const double tstart, const double tstop, const Vec x_stop, Vec x_adj, Vec grad, bool compute_gradient);
};
 
class ExplEuler : public TimeStepper {
  protected:
  Vec stage; 
  public:
    ExplEuler(MasterEq* mastereq_, int ntime_, double total_time_, Output* output_, bool storeFWD_);
 
    ~ExplEuler();
 
    void evolveFWD(const double tstart, const double tstop, Vec x);
 
    void evolveBWD(const double tstart, const double tstop, const Vec x_stop, Vec x_adj, Vec grad, bool compute_gradient);
};
 
class ImplMidpoint : public TimeStepper {
  protected:
  Vec stage, stage_adj; 
  Vec rhs, rhs_adj; 
  KSP ksp; 
  PC  preconditioner; 
  LinearSolverType linsolve_type; 
  int linsolve_maxiter; 
  double linsolve_abstol; 
  double linsolve_reltol; 
  int linsolve_iterstaken_avg; 
  double linsolve_error_avg; 
  int linsolve_counter; 
  Vec tmp, err; 
 
  public:
    ImplMidpoint(MasterEq* mastereq_, int ntime_, double total_time_, LinearSolverType linsolve_type_, int linsolve_maxiter_, Output* output_, bool storeFWD_);
 
    ~ImplMidpoint();
 
    virtual void evolveFWD(const double tstart, const double tstop, Vec x);
 
    virtual void evolveBWD(const double tstart, const double tstop, const Vec x_stop, Vec x_adj, Vec grad, bool compute_gradient);
 
    int NeumannSolve(Mat A, Vec b, Vec x, double alpha, bool transpose);
};
 
class CompositionalImplMidpoint : public ImplMidpoint {
  protected:
 
  std::vector<double> gamma; 
  std::vector<Vec> x_stage; 
  Vec aux; 
  int order; 
 
  public:
    CompositionalImplMidpoint(int order_, MasterEq* mastereq_, int ntime_, double total_time_, LinearSolverType linsolve_type_, int linsolve_maxiter_, Output* output_, bool storeFWD_);
 
    ~CompositionalImplMidpoint();
 
    void evolveFWD(const double tstart, const double tstop, Vec x);
 
    void evolveBWD(const double tstart, const double tstop, const Vec x_stop, Vec x_adj, Vec grad, bool compute_gradient);
};