IDASolver.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: IDASolver.hpp 9727 2014-06-05 10:57:22Z andriushchenko $
 * -----------------------------------------------------------------------------
 *  Description:
 */
 
#ifndef IDA_SOLVER_H
#define IDA_SOLVER_H

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

#include <CapeOpen.hpp>
#include <Diana.hpp>

#ifdef HAVE_HSL
#include <DianaHslSolver.hpp>
#else
#include <DianaSparseSolver.hpp>
#endif

extern "C" {
#if HAVE_SUNDIALS == 1
#include "sundials/sundials_types.h"
#include "sundials/sundials_math.h"
#include "nvector_serial.h"
#include "ida.h"
#include "ida/ida_dense.h"
#include "ida/ida_impl.h"
#elif HAVE_SUNDIALS >= 2
#include "sundials/sundials_types.h"
#include "sundials/sundials_math.h"
#include "nvector/nvector_serial.h"
#include "ida/ida.h"
#include "ida/ida_dense.h"
#include "ida/ida_impl.h"
#else
#error "Unknown internal SUNDIALS version!"
#endif
}

#include <string.h>

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

  protected:

    static int errHandler(int error_code, const char *module, const char *function, char *msg, IDASolver* rdata);

    static int calculateResiduals(Common::Types::CapeDouble Tres, N_Vector Y, N_Vector Yp,
                                  N_Vector R, IDASolver* solver);
#if HAVE_SUNDIALS >= 3

    static int calculateJacobainDense(Common::Types::CapeLong Neq, Common::Types::CapeDouble T, Common::Types::CapeDouble cj, 
                                      N_Vector Y, N_Vector Yp, N_Vector Res, 
                                      DlsMat JJ, IDASolver* solver, 
                                      N_Vector tempv1, N_Vector tempv2, N_Vector tempv3);
#else

    static int calculateJacobainDense(Common::Types::CapeLong Neq, Common::Types::CapeDouble T, N_Vector Y, N_Vector Yp,
                                      N_Vector Res, Common::Types::CapeDouble cj, IDASolver* solver, DenseMat JJ,
                                      N_Vector tempv1, N_Vector tempv2, N_Vector tempv3);
#endif

    static int IDASparseSetup(IDAMem IDA_mem, N_Vector yyp, N_Vector ypp, N_Vector rrp,
                              N_Vector tmp1, N_Vector tmp2, N_Vector tmp3);
      
    static int IDASparseSolve(IDAMem IDA_mem, N_Vector b, N_Vector weight,
                              N_Vector ycur, N_Vector ypcur, N_Vector rrcur);
    
    int  calculateStructuralRank(const Diana::DianaSparseArray& dsa);

    void setValuesFromESO(bool bUpdateAbsTol=true);
  
    void startIntegrator();

    void printStep();

  protected:
    Common::Types::CapeString                    strName;
    Common::Types::CapeString                    strDescription;
    Numeric::Solvers::Model::ICapeNumericModel*  pModel; 
    Numeric::Solvers::Eso::ICapeNumericDAESO*    pESO;   
    Diana::DianaCollection                       collParameters; 
    Common::Parameter::ICapeParameter*           prmLASolver;

    void*                                        memIDA; 

    Common::Types::CapeBoolean                   bStartSimulation; 

    Common::Types::CapeBoolean                   bIntermediateResults;
    Common::Types::CapeBoolean                   bCalculateConsistentIC; 

    Common::Types::CapeLong                      nVerboseLevel;
    Common::Types::CapeLong                                      nParamVisual;
    Common::Types::CapeLong                      nStructRankLHS; 
  
    Common::Types::CapeLong                      nVars; 
    Common::Types::CapeLong                      nEqns; 

    // IDA internal variables
    N_Vector                                     vecStates; 
    N_Vector                                     vecDerivs; 
    N_Vector                                     vecId; 


    N_Vector                                     vecAbsTol; 
  
    Common::Types::CapeDouble                    dblT; 
    Common::Types::CapeDouble                    dblT0; 
    Common::Types::CapeDouble                    dblTend; 
    Common::Types::CapeDouble                    dblTout; 
    Common::Types::CapeDouble                    dblRTol; 
    Common::Types::CapeLong                      nMaxInternalSteps; 

    Common::Types::CapeBoolean                   bSingleStepping; 
    Common::Types::CapeDouble                    dblMaxStepSize; 

    Common::Types::CapeLong                      nMaxOrd; 
  private:
      Common::Types::CapeBoolean                 bIDAInitialized; 
      Common::Types::CapeBoolean                 bIDAInterruption; 
      Diana::DianaSparseArray                    dsaJacoDiff;
#ifdef HAVE_HSL
      Diana::DianaHslSolver                      dssJacoDiff;
#else
      Diana::DianaSparseSolver                   dssJacoDiff;
#endif
      IDALASolverMemRec                          stIDALASolverMem;
  };
};

#endif // IDA_SOLVER_H