SteadyStateContinuation.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: SteadyStateContinuation.hpp 8160 2009-03-16 18:36:30Z miha $
 * -----------------------------------------------------------------------------
 *  Description:
 */

#ifndef STEADY_STATE_CONTINUATION__HPP
#define STEADY_STATE_CONTINUATION__HPP


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

#include <Diana.hpp>
#include <Fortran.hpp>
#include <DianaContinuation.hpp>
#include <vector>
#include <DianaSparseSolver.hpp>

namespace Diana {

  class SteadyStateContinuation : public DianaContinuation {
  public:
    SteadyStateContinuation(Numeric::Solvers::Model::ICapeNumericModel* _pModel);
    virtual ~SteadyStateContinuation();

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

  protected:
    int getModelEquations();

    int getModelVariables();

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

    void packModelVector(Common::Types::CapeArrayDouble& modelVector);

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

    void packModelResidual(Common::Types::CapeArrayDouble& modelResidual);

    void packModelJacobian(Diana::DianaSparseArray& Jaco);

    CapeBoolean calcTestCondition(CapeBoolean bFirstCall=false);

    CapeBoolean calcEigenvalues();
    CapeBoolean calcEigenvaluesLapack();
    CapeBoolean calcEigenvaluesArpack();

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

  protected:
    Common::Types::CapeBoolean      bStabilityCheck;
    Common::Types::CapeBoolean      bStability;
    Common::Types::CapeLong         nConditionTest;
    Common::Types::CapeLong         nConditionCurrent;
    const Common::Types::CapeLong   nEigMonitorMax;
    Common::Types::CapeLong         nEigMonitor;
    Common::Types::CapeLong         nEigConverged;

    Common::Parameter::ICapeParameter* prmEigSolver;

  private:
    Common::Types::CapeArrayDouble  arrReEigenValues;
    Common::Types::CapeArrayDouble  arrImEigenValues;
    Common::Types::CapeLong         nPrevReal;
    Common::Types::CapeLong         nPrevImag;

  private:
    std::vector<fortrand>           lapA;
    std::vector<fortrand>           lapB;
    std::vector<fortrand>           lapAlphaR;
    std::vector<fortrand>           lapAlphaI;
    std::vector<fortrand>           lapBeta;
    std::vector<fortrand>           lapVl;
    std::vector<fortrand>           lapVr;
    std::vector<fortrand>           lapWork;

  private:
    std::vector<fortrand>           arpResid;
    std::vector<fortrand>           arpV;
    std::vector<fortrand>           arpWorkd;
    std::vector<fortrand>           arpWorkl;
    std::vector<fortrand>           arpTmpX;
    std::vector<fortranl>           arpSelect;
    std::vector<fortrand>           arpWorkev;

    Diana::DianaSparseArray         dsaArpAsM;
    Diana::DianaSparseSolver        dssArpM;
    Diana::DianaSparseSolver        dssArpAsM;
  };

};

#endif // STEADY_STATE_CONTINUATION__H