PARDASPKSolver.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: PARDASPKSolver.hpp 6577 2007-04-20 09:58:00Z dosta $
 * -----------------------------------------------------------------------------
 *  Description:
 */

#ifndef PAR_DASPK_SOLVER_H
#define PAR_DASPK_SOLVER_H
#include "PARDASPK_lib.h"

namespace Diana 
{
  class PARDASPKSolver : public Diana::IDianaNumericDAESolver,
                         public Diana::DianaNumericSolver    
 {
  public:
  
   PARDASPKSolver (Numeric::Solvers::Model::ICapeNumericModel* model)
      throw (Common::Error::ECapeUnknown);
   ~PARDASPKSolver();
    
   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);

     
   Numeric::Solvers::Solver::SolveReturn Solve()
      throw (Common::Error::ECapeUnknown, Common::Error::ECapeBadInvOrder,
             Common::Error::ECapeSolvingError, Common::Error::ECapeOutOfResources);          
     
   
   Common::Types::CapeArrayDouble GetSolution()
      throw (Common::Error::ECapeUnknown, Common::Error::ECapeBadInvOrder);
     
          
   void Destroy()
      throw (Common::Error::ECapeUnknown);

   void SetRelTolerance(const Common::Types::CapeArrayDouble& relTolValue)
      throw (Common::Error::ECapeUnknown, Common::Error::ECapeInvalidArgument,
             Common::Error::ECapeOutOfBounds);

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

   void SetAbsTolerance(const Common::Types::CapeArrayDouble& absTolValues)
      throw (Common::Error::ECapeUnknown, Common::Error::ECapeInvalidArgument,
             Common::Error::ECapeOutOfBounds);

   Common::Types::CapeArrayDouble GetAbsTolerance()
      throw (Common::Error::ECapeUnknown, Common::Error::ECapeBadInvOrder);
      
      
      
   Numeric::Solvers::Model::CapeArrayNumericEventInfo*
   AdvanceToNextEvent(Numeric::Solvers::Model::CapeArrayNumericEventInfo endConditions,
                       Common::Types::CapeDouble& timeBefore,
                       Common::Types::CapeDouble& timeAfter)
      throw (Common::Error::ECapeUnknown, Common::Error::ECapeInvalidArgument,
             Common::Error::ECapeOutOfBounds, Common::Error::ECapeBadInvOrder,
             Common::Error::ECapeSolvingError, Common::Error::ECapeOutOfResources,
             Common::Error::ECapeNoImpl);
         
   void setValuesFromESO();
    
   void RecalculateIC();
    
   void SolveWR()
      throw (Common::Error::ECapeUnknown, Common::Error::ECapeBadInvOrder,
             Common::Error::ECapeSolvingError, Common::Error::ECapeOutOfResources); 
                    
protected:  
   static void calculateResiduals( double* T, double Y[], double YPrime[],
      double Delta[], int IPar[] );

   static void calculateJacobian( int* NEQ, double* T, double Y[], 
      double YPRIME[], double* CJ, int IPAR[], int* RECALC );

   static void LinearSolver( int* NEQ, double B[], int IPar[] );        
   
   static void SaveInternalStates( double* X, int* NEQ,
      int* KOLD, double PHI[], double PSI[], int IPar[] );
protected:
   Common::Types::CapeString                    strName;
   Common::Types::CapeString                    strDescription;
    
   Numeric::Solvers::Model::ICapeNumericModel*  pModel; 
   Numeric::Solvers::Eso::ICapeNumericDAESO*    pESO;   
    
                                                                                                                                                
   Diana::DianaCollection                       collParameters; 
   
   Common::Types::CapeBoolean                   m_bStartSimulation; 
   Common::Types::CapeBoolean                   m_bNeedToInitialize;
   CapeBoolean                                  m_bNeedToRestoreState;
   CapeBoolean                                  m_bNeedToSaveState;
    
   Common::Types::CapeLong                  m_lVarsQuantity; 
   Common::Types::CapeLong                  m_lEquationsQuantity; 
   Common::Types::CapeDouble                    dblT; 
   Common::Types::CapeDouble                    dblT0; 
   Common::Types::CapeDouble                    dblTend; 
   Common::Types::CapeLong                      nVerboseLevel; 
    
    
   int                                          nDASKPStatus, nDASKPStatus_old; 
   CapeDouble                                   dblRTol, dblAbsTol;
   
   double                                       *vecRelTol;
   
   double                                       *vecAbsTol;
   
   int*                                         IPar; 
                                                                                                                                                
   CapeDouble                                   m_dLastInternalSolverStep;
   
   /* Structure for data exchange between fortran DASPK and
     C++ Jac, Res ans Solve functions */    
   struct sExchangeStruct
   {
      Numeric::Solvers::Eso::ICapeNumericDAESO* pESO; 
      
      CapeArrayDouble                    JacHistoryTimes;
      CapeArrayDouble                    CJHistory;
      std::vector<DianaSparseArray*>     JacHistory;
      std::vector<DianaSparseSolver*>    SolversHistory;
      CapeLong                           lCurrentIndex;      
      CapeLong                           lOffset;      
      CapeDouble                         m_dIterationNumber;
      CapeBoolean                        bNeedToSaveHistory;
            
      CapeDouble                         dTempJacCalls;
      
      CapeDouble                         dTempSolCalls;
      
      CapeDouble                         dTempResCalls;

      void*                              UpdateReceivedValues;
      void*                              pMemoryBanks;
      void*                              SaveStates;

      //this flag indicate to use backup version of jacobian
      CapeBoolean                        bUseBackupJacobian; 
      DianaSparseArray*                  BackupJacobian;
      DianaSparseSolver*                 BackupSolver;
       
      
      CapeBoolean   bTempTemp;
   } Adresses;
 
   int  nLRW, nLIW; 
   
   double                                *RWORK, *RWORK_old; 
   int                                   *IWORK, *IWORK_old; 
   double                                *Y, *Y_old;         
   double                                *YP, *YP_old;       
   double                                 T0_old;
   int                                   m_nInfo[ 30 ], m_nInfo_old[ 30 ]; 

   CapeDouble                            m_dMaxStep;
   CapeDouble                            m_bNeedToPrint;
   CapeDouble   m_dInitialStepSize;
};
};

#endif // DASPK_SOLVER_H