DianaContinuation.hpp

/* -----------------------------------------------------------------------------
 *  Diana process modelling, simulation and analysis software
 *  Copyright (c) 2005, Michael Krasnyk
 *  All rights reserved.

 *  This file is a part of Diana process modelling, simulation and analysis software
 *
 *  Diana is free software; you can redistribute it and/or modify it
 *  under the terms of the GNU General Public License as published
 *  by the Free Software Foundation (see accompanying file LICENSE)
 * -----------------------------------------------------------------------------
 *  $Id: DianaContinuation.hpp 8959 2009-12-10 12:46:22Z miha $
 * -----------------------------------------------------------------------------
 *  Description:
 */

#ifndef DIANA_CONTINUATION__HPP
#define DIANA_CONTINUATION__HPP

#ifdef HAVE_CONFIG_H
#include <config.h>
#endif

#include <Diana.hpp>
#include <DianaSparseSolver.hpp>
#include <Diana/DianaSignals.hpp>

#ifdef USE_NLEQ1S
#include "NLEQ1S.hpp"
#endif // USE_NLEQ1S

namespace Diana {

  enum ResidualType{
    NoneParameterisationResidual,
    LocalParameterisationResidual,
    PseudoParameterisationResidual,
  };

  enum ParameterisationType{
    ParameterisationLocal,
    ParameterisationPseudoArclength,
    ParameterisationsNumber,
  };

  enum PredictorType{
    PredictorTangent,
    PredictorChord,
    PredictorsNumber,
  };

  char* complex2str(double real, double imag, char* strbuf);
  
  class DianaContinuation : virtual public IDianaContinuation,
                            virtual public DianaNumericSolver {

  public:
    DianaContinuation(Numeric::Solvers::Model::ICapeNumericModel* _pModel);

    virtual ~DianaContinuation();

    const Common::Types::CapeString& GetComponentName() const
      throw (Common::Error::ECapeUnknown);

    void SetComponentName(const Common::Types::CapeString& _name)
      throw (Common::Error::ECapeUnknown, Common::Error::ECapeInvalidArgument);

    const Common::Types::CapeString& GetComponentDescription() const
      throw (Common::Error::ECapeUnknown);

    void SetComponentDescription(const Common::Types::CapeString& _description)
      throw (Common::Error::ECapeUnknown, Common::Error::ECapeInvalidArgument);

    Common::Collection::ICapeCollection* GetParameters()
      throw (Common::Error::ECapeUnknown, Common::Error::ECapeFailedInitialisation,
             Common::Error::ECapeNoImpl);

    void SetSimulationContext(Common::Identification::ICapeIdentification* simContextManager)
      throw (Common::Error::ECapeUnknown, Common::Error::ECapeInvalidArgument,
             Common::Error::ECapeFailedInitialisation, Common::Error::ECapeNoImpl);

    void Initialize()
      throw (Common::Error::ECapeUnknown, Common::Error::ECapeFailedInitialisation,
             Common::Error::ECapeOutOfResources, Common::Error::ECapeLicenceError,
             Common::Error::ECapeBadInvOrder);

    void Terminate()
      throw (Common::Error::ECapeUnknown, Common::Error::ECapeOutOfResources,
             Common::Error::ECapeBadInvOrder);

    Common::Types::CapeArrayDouble GetSolution()
      throw (Common::Error::ECapeUnknown, Common::Error::ECapeBadInvOrder);

    void Destroy()
      throw (Common::Error::ECapeUnknown);

    Numeric::Solvers::Solver::SolveReturn Solve()
      throw (Common::Error::ECapeUnknown, Common::Error::ECapeBadInvOrder,
             Common::Error::ECapeSolvingError, Common::Error::ECapeOutOfResources);

    ContiReturn Continuate()
      throw (Common::Error::ECapeUnknown, Common::Error::ECapeSolvingError);

    void AddFreeParameter(const Common::Types::CapeString& strParName,
                          Common::Types::CapeDouble minValue=0.0,
                          Common::Types::CapeDouble maxValue=0.0)
      throw (Common::Error::ECapeUnknown, Common::Error::ECapeInvalidArgument);

    void RemoveFreeParameter(const Common::Types::CapeString& strParName)
      throw (Common::Error::ECapeUnknown, Common::Error::ECapeInvalidArgument);

    Common::Types::CapeString GetVariableName(Common::Types::CapeLong index)
      throw (Common::Error::ECapeUnknown, Common::Error::ECapeInvalidArgument);

#ifdef DEBUG
    double TestJacobian();
#endif

  protected:
    void             initSystemVariables();
    void             destroySystemVariables();

  protected:
    ContiReturn makeContinuation(const Common::Types::CapeArrayDouble& vecYI,
                                                    Common::Types::CapeDouble initialDirection)
      throw (Common::Error::ECapeUnknown, Common::Error::ECapeSolvingError);

    Numeric::Solvers::Solver::SolveReturn doNewtonSolve(Common::Types::CapeArrayDouble& Y, int nMaxIterations,
                                                        const DianaINTSignalHandler& inthandler)
                        throw (Common::Error::ECapeUnknown, Common::Error::ECapeSolvingError);

    Numeric::Solvers::Solver::SolveReturn doNewtonSolve_(Common::Types::CapeArrayDouble& Y, int nMaxIterations,
                                                         const DianaINTSignalHandler& inthandler)
      throw (Common::Error::ECapeUnknown, Common::Error::ECapeSolvingError);
    
    void restoreInitialValues()
      throw (Common::Error::ECapeUnknown);

    virtual int getModelEquations() = 0;

    virtual int getModelVariables() = 0;

    virtual CapeBoolean calcTestCondition(CapeBoolean bFirstCall=false);

    CapeBoolean checkIsolaCondition(const CapeArrayDouble& vecY0, const CapeArrayDouble& vecY1,
                             const CapeArrayDouble& vecY2);

    void formFunction(const Common::Types::CapeArrayDouble& Y, Common::Types::CapeArrayDouble& Res);

  private:
    void packVector(Common::Types::CapeArrayDouble& work_vector);

    void unpackVector(const Common::Types::CapeArrayDouble& work_vector);

    void packFreeParameters(int nFreeParamsPos, Common::Types::CapeArrayDouble& modelVector);

    void unpackFreeParameters(int nFreeParamsPos, const Common::Types::CapeArrayDouble& modelVector);

    void getParameterBoundaries(Common::Parameter::ICapeParameter* param,
                                Common::Types::CapeDouble& dblMin,
                                Common::Types::CapeDouble& dblMax);

  protected:

    virtual void packModelVector(Common::Types::CapeArrayDouble& modelVector) = 0;

    virtual void packModelBoundaries(Common::Types::CapeArrayDouble& arrMinY,
                                     Common::Types::CapeArrayDouble& arrMaxY) = 0;

    virtual void unpackModelVector(const Common::Types::CapeArrayDouble& modelVector) = 0;

    virtual void packModelResidual(Common::Types::CapeArrayDouble& modelResidual) = 0;

    virtual void packModelJacobian(Diana::DianaSparseArray& Jaco);

    virtual void addDerivativeByParameter(Diana::DianaSparseArray& Jaco, int nJacCol, int nParamIndex);

#ifdef USE_NLEQ1S
  private:
    Numeric::Solvers::Solver::SolveReturn doNewtonSolve_NLEQ1S(Common::Types::CapeArrayDouble& Y, int nMaxIterations,
                                                               const DianaINTSignalHandler& inthandler)
      throw (Common::Error::ECapeUnknown, Common::Error::ECapeSolvingError);
    static void NLEQ1SRes(fortrani& n, fortrand* x, fortrand* f, fortrani& ifail, DianaContinuation* param);
    static void NLEQ1SJac(fortrani& n, fortrand* x, fortrand* dfdx, fortrani* irow, fortrani* icol, fortrani& nfill, 
                          fortrani& ifail, DianaContinuation* param);
#endif // USE_NLEQ1S

  protected:
    int              checkVariableFiniteness(const Common::Types::CapeArrayDouble& Y);

    void             calcTangentVector(const Common::Types::CapeArrayDouble& Y, Common::Types::CapeArrayDouble& T,
                                       PredictorType type);

    CapeBoolean      checkVariableBoundaries(const Common::Types::CapeArrayDouble& Y);

    int              makeBoundaryCorrection(Common::Types::CapeArrayDouble& Y, const DianaINTSignalHandler& inthandler);

  private:
    /* parametrisation methods */
    void             setNoneParameterisation();

    /* local parameterisation */
    void             setLocalParameterisation(CapeLong idx, CapeDouble val, CapeDouble sign);
    Common::Types::CapeDouble calcLocalParameterisation(const Common::Types::CapeArrayDouble& dblY);
    void             addLocalParameterisationJacRow(Diana::DianaSparseArray& Jaco, int nRow);

    /* pseudo-arclength parameterisation */
    void             setPseudoParameterisation();
    Common::Types::CapeDouble calcPseudoParameterisation(const Common::Types::CapeArrayDouble& dblY);
    void             addPseudoParameterisationJacRow(Diana::DianaSparseArray& Jaco, int nRow);

  private:
    // Linear algebra
    void             formJacobianMatrix(const Common::Types::CapeArrayDouble& Y);

  protected:
    Common::Types::CapeDouble  vecNorm2(const Common::Types::CapeArrayDouble& dblY);
    Common::Types::CapeBoolean isZeroValue(Common::Types::CapeDouble value);
    Common::Types::CapeDouble  calcCosAngle(const Common::Types::CapeArrayDouble& Y1, const Common::Types::CapeArrayDouble& Y0,
                                            const Common::Types::CapeArrayDouble& Ym1);

    void printVector(const Common::Types::CapeString& name, const Common::Types::CapeArrayDouble& vec);

    void calcArpackSMEigenpair(Diana::DianaSparseArray&        dsaA,
                               Common::Types::CapeArrayDouble& arrEigenvector);

    std::vector<std::complex<Common::Types::CapeDouble> > GetEigenvalues()
      throw (Common::Error::ECapeUnknown);

    ContiReturn init(const Common::Types::CapeArrayDouble& vecYI,
                     Common::Types::CapeDouble initialDirection,
                     const int nInitialCorrectionDir,
                     const DianaINTSignalHandler& inthandler);
    ContiReturn one_step(const DianaINTSignalHandler& inthandler);

  protected:
    Common::Types::CapeString                   strName;
    Common::Types::CapeString                   strDescription;
    Numeric::Solvers::Model::ICapeNumericModel* pModel;               
    Diana::IDianaDAESO*                         pESO;                 
    Diana::DianaCollection*                     pESORealParameters;
    Diana::DianaCollection                      collSolverParameters; 
    CapeBoolean                                 bInitialized;         

    CapeBoolean                                 bstart;
    CapeBoolean                                 bintermediate;


  protected:
    // parameters lists
    Common::Types::CapeArrayLong          arrFreeParams;     
    Common::Types::CapeArrayDouble        arrFreeParamsMin;  
    Common::Types::CapeArrayDouble        arrFreeParamsMax;  

  protected:
    Common::Types::CapeArrayDouble        dblMinY;      
    Common::Types::CapeArrayDouble        dblMaxY;      
  private:
    Common::Types::CapeArrayDouble        dblScaleY;    

  private:
    // model variables
    int                            nVars;               
    int                            nModelVars;          
    int                            nModelEqns;          
    Common::Types::CapeArrayDouble vecY;                
    Common::Types::CapeArrayDouble vecYBC;              
    Common::Types::CapeArrayDouble vecRes;              
    Common::Types::CapeArrayDouble vecEye;              
    Common::Types::CapeArrayDouble vecY0;               
    Common::Types::CapeArrayDouble vecT0;               
    Common::Types::CapeArrayDouble vecYm1;              
    Common::Types::CapeArrayDouble vecYStart;           
    Diana::DianaSparseArray        matJacobian;         

    Diana::DianaSparseSolver       spmNewton;

  private:
    // variables for predictor and parameterisation functions
    ResidualType                nResidualType;          

    Common::Types::CapeLong     nParameterisationType;  
    Common::Types::CapeLong     nPredictorType;         
    Common::Parameter::ICapeParameter* prmParameterisationType;
    Common::Parameter::ICapeParameter* prmPredictorType;
    Common::Parameter::ICapeParameter* prmInitialCorrection;

    // local parameterisation internal parameters
    CapeDouble                  dblLocalXk;             
    CapeLong                    indLocalXk;             

  protected:
    // continuation paramameters (access with help of GetParameters method)
    Common::Types::CapeLong     nVerboseLevel;          
    Common::Types::CapeBoolean  bBidirectional;         
    Common::Types::CapeDouble   dblInitialDirection;    

    Common::Types::CapeDouble   dblStepSize;            
    Common::Types::CapeDouble   dblInitialStepSize;     
    Common::Types::CapeDouble   dblFailedStepScale;     
    Common::Types::CapeDouble   dblSucceededStepScale;  
    Common::Types::CapeDouble   dblMinStepSize;         
    Common::Types::CapeDouble   dblMaxStepSize;         
    Common::Types::CapeLong     nMaxStepsAmount;        
    Common::Types::CapeLong     nStepNumber;            
    Common::Types::CapeLong     nMaxNewtonIts;          
    Common::Types::CapeDouble   dblTol;                 
    Common::Types::CapeDouble   dblTolF;                
    Common::Types::CapeLong     nSuccessIts;            

    Common::Types::CapeLong     nBoundaryVariableIdx;

    Common::Types::CapeDouble   dblFDEpsilon;           
    Common::Types::CapeDouble   dblUround;              
    Common::Types::CapeDouble   dblIsolaTrustInterval;  

    Common::Types::CapeBoolean  bOrthogonalBranch;
    Common::Types::CapeBoolean  bCheckVariablesBoundaries; 
  };
};

#endif // DIANA_CONTINUATION__HPP