added Cvode solver source code to cpp runtime

git-svn-id: https://openmodelica.org/svn/OpenModelica/trunk@10323 f25d12d1-65f4-0310-ae8a-bbce733d8d8e

SimulationRuntime/cpp/Source/Solver/CVode/Implementation/CVode.h

@@ -0,0 +1,196 @@
+
+#pragma once
+
+#define BOOST_EXTENSION_SOLVER_DECL BOOST_EXTENSION_IMPORT_DECL
+
+#include "Solver/Implementation/SolverDefaultImplementation.h"
+#include<cvodes.h>
+#include<nvector_serial.h>
+#include<cvodes_dense.h>
+
+
+
+/*****************************************************************************/
+// Cvode aus dem SUNDIALS-Package
+// BDF-Verfahren für steife und nicht-steife ODEs
+// Dokumentation siehe offizielle Cvode Doku
+
+/*****************************************************************************
+Copyright (c) 2004, Bosch Rexroth AG, All rights reserved
+*****************************************************************************/
+class ICVodeSettings;
+
+class Cvode 
+	: public IDAESolver,  public SolverDefaultImplementation
+{
+public:
+
+	Cvode(IDAESystem* system, ISolverSettings* settings);
+
+	virtual ~Cvode();
+
+	// geerbt von Object (in SolverDefaultImplementation)
+	//---------------------------------------
+	/// Spezielle Solvereinstellungen setzten (default oder user defined)
+	virtual void init();
+
+
+	// geerbt von ISolver
+	//---------------------------------------
+	/// Setzen der Startzeit für die numerische Lösung
+	virtual void setStartTime(const double& time)
+	{
+		SolverDefaultImplementation::setStartTime(time);
+	};
+
+	/// Setzen der Endzeit für die numerische Lösung
+	virtual void setEndTime(const double& time)
+	{
+		SolverDefaultImplementation::setEndTime(time);
+	};
+
+	/// Setzen der initialen Schrittweite (z.B. auch nach Nullstelle)
+	virtual void setInitStepSize(const double& stepSize)
+	{
+		SolverDefaultImplementation::setInitStepSize(stepSize);
+	};
+
+	/// Berechung der numerischen Lösung innerhalb eines gegebenen Zeitintervalls
+	virtual void solve(const SOLVERCALL command = UNDEF_CALL);
+
+	/*/// Liefert den Zeitpunkt des letzten erfolgreichen Zeitschrittes (kann ~= tend sein)
+	virtual const double& getCurrentTime()
+	{
+	return (SolverDefaultImplementation::getCurrentTime());
+	};*/
+
+	/// Liefert den Status des Solvers nach Beendigung der Simulation
+	virtual const IDAESolver::SOLVERSTATUS getSolverStatus()
+	{
+		return (SolverDefaultImplementation::getSolverStatus());
+	};
+
+	/*/// Liefert Anzahl der Schritte im Zeitintervall
+	virtual void giveNumberOfSteps(int& totStps, int& accStps, int& rejStps)
+	{
+	SolverDefaultImplementation::giveNumberOfSteps(totStps, accStps, rejStps);
+	};
+
+	/// Liefert Anzahl der Schritte der Nullstellensuche und die Anzahl der Nullstellen im gesamten Integrationsintervall
+	virtual void giveZeroSteps(int& zeroSearchStps, int& numberOfZeros)
+	{
+	SolverDefaultImplementation::giveZeroSteps(zeroSearchStps,numberOfZeros);
+	};*/
+
+
+	/*/// Veranlasst, dass die Ausgaberoutine des Solvers (writeSolverOutput) nach jedem Schritt aufgerufen wird
+	virtual void activateSolverOutput(const bool& flag)
+	{
+	SolverDefaultImplementation::activateSolverOutput(flag);
+	};
+
+	/// Output-Stream, in den die Ausgabe des Solvers erfolgen soll, setzen
+	virtual void setOutputStream(ostream& outputStream)
+	{
+	SolverDefaultImplementation::setOutputStream(outputStream);
+	};
+	*/
+	//// Ausgabe von statistischen Informationen (wird vom SimManager nach Abschluß der Simulation aufgerufen)
+	virtual void writeSimulationInfo(ostream& outputStream);
+
+	/// Anfangszustand (und entsprechenden Zeitpunkt) des Systems (als Kopie) speichern
+	virtual void saveInitState();
+
+	/// Anfangszustand (und entsprechenden Zeitpunkt) zurück ins System kopieren
+	virtual void restoreInitState();
+
+	/// Letzten gültigen Zustand (und entsprechenden Zeitpunkt) des Systems (als Kopie) speichern
+	virtual void saveLastSuccessfullState();
+
+	/// Letzten gültigen Zustand  (und entsprechenden Zeitpunkt) zurück ins System kopieren
+	virtual void restoreLastSuccessfullState();
+
+	/// speichert den Zustand (und entsprechenden Zeitpunkt) des Systems (als Kopie) nach einem "großen Schritt" bei partitionierter Integration
+	virtual void saveLargeStepState();
+
+	/// liefert den normierten Fehler zwischen aktuellem Zustand und Zustand nach einem "großen Schritt" bei partitionierter Integration
+	virtual void giveScaledError(const double& h, double& error);
+
+	/// Approximation höherer Ordnung des Zustandes berechnen und in System kopieren
+	virtual void refineCurrentState(const double& r);
+	virtual const int reportErrorMessage(ostream& messageStream);
+private:
+
+	// Solveraufruf
+	void callCVode(); 
+
+	/// Kapselung der Berechnung der rechten Seite 
+	void calcFunction(const double& time, const double* y, double* yd);
+
+	// Callback für die rechte Seite
+	static int CV_fCallback(double t, N_Vector y, N_Vector ydot, void *user_data);
+
+	// Checks error flags of SUNDIALS
+	int check_flag(void *flagvalue, char *funcname, int opt);
+
+	// Nulltellenfunktion
+	void giveZeroVal(const double &t,const double *y,double *zeroValue);
+
+	// Callback der Nullstellenfunktion
+	static int CV_ZerofCallback(double t, N_Vector y, double *zeroval, void *user_data);
+
+
+	ICVodeSettings
+		*_cvodesettings;							///< Input			- Solver settings
+
+	void
+		*_cvodeMem,									///< Temp			- Memory for the solver
+		*_data;										///< Temp			- User data. Contains pointer to Cvode
+
+	long int
+		_dimSys,									///< Input 			- (total) Dimension of system (=number of ODE)
+		_idid,										///< Input, Output	- Status Flag
+		_locStps;									///< Output			- Number of Steps between two events
+
+	int
+		_outStps;									///< Output			- Total number of output-steps
+
+
+	double
+		*_y,										///< Output			- (Current) State vector
+		*_y0,										///< Temp			- (Old) state vector at left border of intervall (last step)
+		*_y1,										///< Temp			- (New) state vector at right border of intervall (last step)
+		*_yInit,									///< Temp			- Initial state vector
+		*_yLastSucess,								///< Temp			- State vector of last successfull step 
+		*_yLargeStep,								///< Temp			- State vector of "large step" (used by "coupling step size controller" of SimManger)
+		*_yWrite,									///< Temp			- Zustand den das System rausschreibt
+		*_f0,
+		*_f1;
+
+	double
+		_hOut,										///< Temp			- Ouput step size for dense output
+		_hZero,										///< Temp			- Downscale of step size to approach the zero crossing
+		_hUpLim,									///< Temp 			- Maximal step size
+		_hZeroCrossing,								///< Temp 			- Stores the current step size (at the time the zero crossing occurs)
+		_hUpLimZeroCrossing;						///< Temp 			- Stores the upper limit of the step size (at the time the zero crossing occurs), because it is changed for zero search
+
+
+
+	double 
+		_tOut,										///< Output			- Time for dense output
+		_tHelp,										///< Temp			- Help variable
+		_tLastZero, 								///< Temp			- Stores the time of the last zero (not last zero crossing!)
+		_tRealInitZero,								///< Temp			- Time of the very first zero in all zero functions
+		_doubleZeroDistance,						///< Temp			- In case of two zeros in one intervall (doubleZero): distance between zeros
+		_tZero,										///< Temp			- Nullstelle
+		_tLastWrite;								///< Temp			- Letzter Ausgabezeitpunkt
+
+	bool
+		_doubleZero;								///< Temp			- Flag to denote two zeros in intervall
+
+	N_Vector 
+		_CV_y0,									///< Temp			- Initial values in the Cvode Format
+		_CV_y;									///< Temp			- State in Cvode Format 
+
+};
+

SimulationRuntime/cpp/Source/Solver/CVode/Implementation/CVodeSettings.cpp

@@ -0,0 +1,29 @@
+#pragma once
+#include "stdafx.h"
+#include "CVodeSettings.h"
+
+CVodeSettings::CVodeSettings(IGlobalSettings* globalSettings)
+	: SolverSettings		(globalSettings)
+
+{
+};
+
+bool CVodeSettings::getDenseOutput()
+{
+	return _denseOutput;
+}
+void CVodeSettings::setDenseOutput(bool dense)
+{
+	_denseOutput = dense;
+}	
+
+
+bool CVodeSettings::getEventOutput()
+{
+	return _eventOutput;
+}
+
+void CVodeSettings::setEventOutput(bool eventOutput)
+{
+	_eventOutput = eventOutput;
+}

SimulationRuntime/cpp/Source/Solver/CVode/Implementation/CVodeSettings.h

@@ -0,0 +1,24 @@
+#pragma once
+#define BOOST_EXTENSION_SOLVERSETTINGS_DECL BOOST_EXTENSION_IMPORT_DECL
+#include "../../Implementation/SolverSettings.h"
+#include "../Interfaces/ICVodeSettings.h"
+
+class CVodeSettings : public ICVodeSettings, public  SolverSettings
+{
+
+public:
+	CVodeSettings(IGlobalSettings* globalSettings);
+	/**
+	Equidistant output(by interpolation polynominal) ([true,false]; default: false)
+	*/
+	 virtual bool getDenseOutput();
+	 virtual void setDenseOutput(bool);	
+	 virtual bool getEventOutput();
+	 virtual void setEventOutput(bool);	
+private:
+	 	bool
+		_denseOutput,			///< Equidistant output(by interpolation polynominal) ([true,false]; default: false)
+		_eventOutput;
+
+
+};

SimulationRuntime/cpp/Source/Solver/CVode/Implementation/stdafx.h

@@ -0,0 +1,31 @@
+// stdafx.h : Includedatei für Standardsystem-Includedateien
+// oder häufig verwendete projektspezifische Includedateien,
+// die nur in unregelmäßigen Abständen geändert werden.
+//
+
+#pragma once
+
+#include <map>
+#include <boost/ref.hpp>
+#include <boost/bind.hpp>
+#include <boost/function.hpp>
+#include <boost/numeric/ublas/vector.hpp>
+#include <boost/numeric/ublas/io.hpp>
+#include <boost/numeric/ublas/matrix.hpp>
+#include <boost/numeric/ublas/matrix_proxy.hpp>
+#include <boost/archive/xml_oarchive.hpp>
+#include <boost/archive/xml_iarchive.hpp>
+#include <boost/serialization/nvp.hpp>
+#include <string>
+#include <fstream>
+#include <boost/extension/extension.hpp>
+#include <typeinfo>
+#include <boost/extension/type_map.hpp>
+#include <boost/extension/factory.hpp>
+using namespace boost::numeric;
+using namespace std;
+using std::ios;
+
+
+
+// TODO: Hier auf zusätzliche Header, die das Programm erfordert, verweisen.

SimulationRuntime/cpp/Source/Solver/CVode/Interfaces/CVodeSettings.h

@@ -0,0 +1,61 @@
+#pragma once
+
+#include "../../Implementation/SolverSettings.h"
+
+
+						/*****************************************************************************/
+						/**
+
+						Klasse zur Kapselung der Parameter (Einstellungen) für Euler.
+						Hier werden default-Einstellungen entsprechend der allgemeinen Simulations-
+						einstellugnen gemacht, diese können überprüft und ev. Fehleinstellungen korrigiert 
+						werden.
+						\date     Montag, 1. August 2005
+						\author   Daniel Kanth
+						$Id: CVodeSettings.h,v 1.1 2007/04/25 02:18:37 danikant Exp $
+						\par	Synopsis:
+						CVodeSettings(&System, &Global);
+						*/
+						/*****************************************************************************
+						Copyright (c) 2004, Bosch Rexroth AG, All rights reserved
+						*****************************************************************************/
+						class CVodeSettings : public SolverSettings
+						{
+
+						public:
+							CVodeSettings(IGlobalSettings* globalSettings)
+								: SolverSettings		(globalSettings)
+								, iMethod				(EULERFORWARD)
+								, dIterTol				(1e-8)
+								,bContinue				(false)
+								, iJacUpdate				(0)
+							{
+							};
+
+							/// Enum für Wahl der Methode
+							enum EulerMethod
+							{
+								EULERFORWARD,  			///< Euler Vorwärts
+								EULERBACKWARD,			///< Euler Rückwärts (wie Euler-Cauchy, aber impl.)
+
+							};
+
+							int 
+								iMethod;				///< Verfahrensauswahl ([0,1,2,3,4,5]; default: 0)
+
+
+							double
+								dIterTol;				///< Toleranz für die Iteration (Der Steigungen werden iteriert, bis Steigungsänderung < IterTol)
+
+							bool 
+								bContinue;             ///< Soll nach Fehler in Newtoniteration weiter gerechnet werden oder abgebrochen werden
+							int 
+								iJacUpdate;			   //Anzahl Schritte wann Jacobi Matrix aktualisiert werden soll	
+
+						};
+					}
+				}
+			}
+		}
+	}
+}

SimulationRuntime/cpp/Source/Solver/CVode/Interfaces/ICVodeSettings.h

@@ -0,0 +1,10 @@
+#pragma once
+
+class ICVodeSettings
+{
+public:
+	virtual bool getDenseOutput() =0;
+	virtual void setDenseOutput(bool) =0;	
+	  virtual bool getEventOutput() = 0;
+	 virtual void setEventOutput(bool)=0;
+};

SimulationRuntime/cpp/Source/Solver/Implementation/SolverDefaultImplementation.cpp

@@ -16,7 +16,7 @@ SolverDefaultImplementation::SolverDefaultImplementation(IDAESystem* system, ISo
 , _tEnd					(0.0)
 , _tLastSuccess			(0.0)
 , _tLastUnsucess		(0.0)
-
+, _tLargeStep			(0.0)
 , _h					(0.0)
 
 //, _firstCall			(true)
@@ -33,6 +33,7 @@ SolverDefaultImplementation::SolverDefaultImplementation(IDAESystem* system, ISo
 , _dimZeroFunc			(0)
 , _zeroVal				(NULL)
 , _zeroValLastSuccess	(NULL)
+, _zeroValLargeStep		(NULL)
 , _events				(NULL)
 , _zeroSearchActive		(false)
 
@@ -117,11 +118,12 @@ void SolverDefaultImplementation::init()
 		_zeroValLastSuccess	= new double[_dimZeroFunc];
 		_events				= new bool[_dimZeroFunc];
 		_zeroValInit			= new double[_dimZeroFunc];
-
+		_zeroValLargeStep		= new double[_dimZeroFunc];
 		event_system->giveZeroFunc(_zeroVal,_settings->getZeroTol());
 		memcpy(_zeroValLastSuccess,_zeroVal,_dimZeroFunc*sizeof(double));
 		memcpy(_zeroValInit,_zeroVal,_dimZeroFunc*sizeof(double));
 		memset(_events,false,_dimZeroFunc*sizeof(bool));
+		memcpy(_zeroValLargeStep,_zeroVal,_dimZeroFunc*sizeof(double));
 	}
      _time_events = event_system->getTimeEvents();
 	 // Set flags