Skip to content

Commit

Permalink
- removed the old Kinsol_lapack - flag of the cpp-runtime
Browse files Browse the repository at this point in the history
- KinsolLapack will now work with new and old versions of sundials
- Non linear solver can now continue simulation instead of throwing an exception, if it does not reach the given precision, by passing a flag to the executable
  • Loading branch information
Marcus Walther committed Aug 12, 2015
1 parent 3e37674 commit 12f2912
Show file tree
Hide file tree
Showing 18 changed files with 152 additions and 84 deletions.
10 changes: 0 additions & 10 deletions SimulationRuntime/cpp/CMakeLists.txt
Expand Up @@ -18,7 +18,6 @@
# if write output should be handled in parallel -DUSE_PARALLEL_OUTPUT=ON
# if ScoreP should be used for performance analysis -DUSE_SCOREP=ON
# if the boost libraries should be linked statically -DBOOST_STATIC_LINKING=ON
# if the lapack functions and data structurs of sundials should be used -DSUNDIALS_LAPACK=ON
# if boost libraries should be linked against absolute path libraries -DUSE_BOOST_REALPATHS=ON
#
# Example: "cmake -DCMAKE_BUILD_TYPE=RelWithDebInfo" to create statically linked libraries
Expand All @@ -38,7 +37,6 @@ OPTION(USE_SCOREP "USE_SCOREP" OFF)
OPTION(BOOST_STATIC_LINKING "BOOST_STATIC_LINKING" OFF)
OPTION(RUNTIME_PROFILING "RUNTIME_PROFILING" OFF)
OPTION(FMU_KINSOL "FMU_KINSOL" OFF)
OPTION(SUNDIALS_LAPACK "SUNDIALS_LAPACK" OFF)
OPTION(SUNDIALS_ROOT "SUNDIALS ROOT" "")

#Set Variables
Expand Down Expand Up @@ -104,14 +102,6 @@ ELSE(NOT BOOST_STATIC_LINKING)
ADD_DEFINITIONS(-DBOOST_STATIC_LINKING)
ENDIF(NOT BOOST_STATIC_LINKING)

# Handle sundials lapack
IF(NOT SUNDIALS_LAPACK)
MESSAGE(STATUS "Sundials lapack functions disabled")
ELSE(NOT SUNDIALS_LAPACK)
MESSAGE(STATUS "Using Sundials lapack functions")
ADD_DEFINITIONS(-DUSE_SUNDIALS_LAPACK)
ENDIF(NOT SUNDIALS_LAPACK)

# Handle parallel output
IF(USE_PARALLEL_OUTPUT)
ADD_DEFINITIONS(-DUSE_PARALLEL_OUTPUT)
Expand Down
1 change: 1 addition & 0 deletions SimulationRuntime/cpp/Core/SimController/SimController.cpp
Expand Up @@ -242,6 +242,7 @@ void SimController::Start(SimSettings simsettings, string modelKey)
global_settings->setLogSettings(simsettings.logSettings);
global_settings->setAlarmTime(simsettings.timeOut);
global_settings->setOutputPointType(simsettings.outputPointType);
global_settings->setNonLinearSolverContinueOnError(simsettings.nonLinearSolverContinueOnError);

/*boost::shared_ptr<SimManager>*/ _simMgr = boost::shared_ptr<SimManager>(new SimManager(mixedsystem, _config.get()));

Expand Down
11 changes: 11 additions & 0 deletions SimulationRuntime/cpp/Core/SimulationSettings/GlobalSettings.cpp
Expand Up @@ -19,6 +19,7 @@ GlobalSettings::GlobalSettings()
, _selected_nonlin_solver("Newton")
, _resultsfile_name("results.csv")
, _endless_sim(false)
, _nonLinSolverContinueOnError(false)
, _outputPointType(OPT_ALL)
, _alarm_time(0)
{
Expand Down Expand Up @@ -189,4 +190,14 @@ unsigned int GlobalSettings::getAlarmTime()
{
return _alarm_time;
}

void GlobalSettings::setNonLinearSolverContinueOnError(bool value)
{
_nonLinSolverContinueOnError = value;
}

bool GlobalSettings::getNonLinearSolverContinueOnError()
{
return _nonLinSolverContinueOnError;
}
/** @} */ // end of coreSimulationSettings
1 change: 1 addition & 0 deletions SimulationRuntime/cpp/Core/System/AlgLoopSolverFactory.cpp
Expand Up @@ -51,6 +51,7 @@ boost::shared_ptr<IAlgLoopSolver> AlgLoopSolverFactory::createAlgLoopSolver(IAlg

string nonlinsolver_name = _global_settings->getSelectedNonLinSolver();
boost::shared_ptr<INonLinSolverSettings> algsolversetting= createNonLinSolverSettings(nonlinsolver_name);
algsolversetting->setContinueOnError(_global_settings->getNonLinearSolverContinueOnError());
_algsolversettings.push_back(algsolversetting);

boost::shared_ptr<IAlgLoopSolver> algsolver= createNonLinSolver(algLoop,nonlinsolver_name,algsolversetting);
Expand Down
Expand Up @@ -19,6 +19,7 @@ struct SimSettings
unsigned int timeOut;
OutputPointType outputPointType;
LogSettings logSettings;
bool nonLinearSolverContinueOnError;
};

/**
Expand Down
Expand Up @@ -53,6 +53,9 @@ class GlobalSettings : public IGlobalSettings
virtual void setAlarmTime(unsigned int);
virtual unsigned int getAlarmTime();

virtual void setNonLinearSolverContinueOnError(bool);
virtual bool getNonLinearSolverContinueOnError();

private:
double
_startTime, ///< Start time of integration (default: 0.0)
Expand All @@ -61,7 +64,8 @@ class GlobalSettings : public IGlobalSettings
bool
_resultsOutput, ///< Write out results ([false,true]; default: true)
_infoOutput, ///< Write out statistical simulation infos, e.g. number of steps (at the end of simulation); [false,true]; default: true)
_endless_sim;
_endless_sim,
_nonLinSolverContinueOnError;
string
_output_path,
_selected_solver,
Expand Down
Expand Up @@ -86,5 +86,8 @@ class IGlobalSettings
virtual string getResultsFileName() = 0;
virtual void setRuntimeLibrarypath(string) = 0;
virtual string getRuntimeLibrarypath() = 0;

virtual void setNonLinearSolverContinueOnError(bool) = 0;
virtual bool getNonLinearSolverContinueOnError() = 0;
};
/** @} */ // end of coreSimulationSettings
Expand Up @@ -25,5 +25,7 @@ class INonLinSolverSettings
virtual double getDelta() = 0;
virtual void setDelta(double) = 0;
virtual void load(string) = 0;
virtual void setContinueOnError(bool) = 0;
virtual bool getContinueOnError() = 0;
};
/** @} */ // end of coreSolver
/** @} */ // end of coreSolver
14 changes: 9 additions & 5 deletions SimulationRuntime/cpp/Include/Solver/Hybrj/HybrjSettings.h
Expand Up @@ -22,11 +22,15 @@ class HybrjSettings :public INonLinSolverSettings
virtual double getDelta();
virtual void setDelta(double);
virtual void load(string);

virtual void setContinueOnError(bool);
virtual bool getContinueOnError();
private:
long int iNewt_max; ///< max. Anzahl an Newtonititerationen pro Schritt (default: 25)
long int _iNewt_max; ///< max. Anzahl an Newtonititerationen pro Schritt (default: 25)

double dRtol; ///< Relative Toleranz für die Newtoniteration (default: 1e-6)
double dAtol; ///< Absolute Toleranz für die Newtoniteration (default: 1e-6)
double dDelta; ///< Dämpfungsfaktor (default: 0.9)
double _dRtol; ///< Relative Toleranz für die Newtoniteration (default: 1e-6)
double _dAtol; ///< Absolute Toleranz für die Newtoniteration (default: 1e-6)
double _dDelta; ///< Dämpfungsfaktor (default: 0.9)
bool _continueOnError;
};
/** @} */ // end of solverHybrj
/** @} */ // end of solverHybrj
6 changes: 5 additions & 1 deletion SimulationRuntime/cpp/Include/Solver/Kinsol/KinsolLapack.h
Expand Up @@ -7,7 +7,11 @@
int KINLapackCompletePivoting(void *kinmem, int N);
static int KINLapackCompletePivotingInit(KINMem kin_mem);
static int KINLapackCompletePivotingSetup(KINMem kin_mem);
static int KINLapackCompletePivotingSolve(KINMem kin_mem, N_Vector x, N_Vector b,realtype *res_norm/*used in new sundials: realtype *sJpnorm, realtype *sFdotJp*/);
#if SUNDIALS_MAJOR_VERSION > 2 || (SUNDIALS_MAJOR_VERSION == 2 && SUNDIALS_MINOR_VERSION > 5)
static int KINLapackCompletePivotingSolve(KINMem kin_mem, N_Vector x, N_Vector b, realtype *sJpnorm, realtype *sFdotJp*/);
#else
static int KINLapackCompletePivotingSolve(KINMem kin_mem, N_Vector x, N_Vector b, realtype *res_norm);
#endif
static void KINLapackCompletePivotingFree(KINMem kin_mem);
static int calcJacobian(KINMem kin_mem);

Expand Down
14 changes: 9 additions & 5 deletions SimulationRuntime/cpp/Include/Solver/Kinsol/KinsolSettings.h
Expand Up @@ -22,11 +22,15 @@ class KinsolSettings :public INonLinSolverSettings
virtual double getDelta();
virtual void setDelta(double);
virtual void load(string);

virtual void setContinueOnError(bool);
virtual bool getContinueOnError();
private:
long int iNewt_max; ///< max. Anzahl an Newtonititerationen pro Schritt (default: 25)
long int _iNewt_max; ///< max. Anzahl an Newtonititerationen pro Schritt (default: 25)

double dRtol; ///< Relative Toleranz für die Newtoniteration (default: 1e-6)
double dAtol; ///< Absolute Toleranz für die Newtoniteration (default: 1e-6)
double dDelta; ///< Dämpfungsfaktor (default: 0.9)
double _dRtol; ///< Relative Toleranz für die Newtoniteration (default: 1e-6)
double _dAtol; ///< Absolute Toleranz für die Newtoniteration (default: 1e-6)
double _dDelta; ///< Dämpfungsfaktor (default: 0.9)
bool _continueOnError;
};
/** @} */ // end of solverKinsol
/** @} */ // end of solverKinsol
12 changes: 8 additions & 4 deletions SimulationRuntime/cpp/Include/Solver/Newton/NewtonSettings.h
Expand Up @@ -24,11 +24,15 @@ class NewtonSettings :public INonLinSolverSettings
virtual double getDelta();
virtual void setDelta(double);
virtual void load(string);

virtual void setContinueOnError(bool);
virtual bool getContinueOnError();
private:
long int iNewt_max; ///< max. Anzahl an Newtonititerationen pro Schritt (default: 25)
long int _iNewt_max; ///< max. Anzahl an Newtonititerationen pro Schritt (default: 25)

double dRtol; ///< Relative Toleranz für die Newtoniteration (default: 1e-6)
double dAtol; ///< Absolute Toleranz für die Newtoniteration (default: 1e-6)
double dDelta; ///< Dämpfungsfaktor (default: 0.9)
double _dRtol; ///< Relative Toleranz für die Newtoniteration (default: 1e-6)
double _dAtol; ///< Absolute Toleranz für die Newtoniteration (default: 1e-6)
double _dDelta; ///< Dämpfungsfaktor (default: 0.9)
bool _continueOnError;
};
/** @} */ // end of solverNewton
Expand Up @@ -62,6 +62,7 @@ SimSettings OMCFactory::readSimulationParameter(int argc, const char* argv[])
desc.add_options()
("help", "produce help message")
("emit_protected", "emits protected variables to the result file")
("nls_continue", po::bool_switch()->default_value(false),"non linear solver will continue if it can not reach the given precision")
("runtime-library,R", po::value<string>(),"path to cpp runtime libraries")
("Modelica-system-library,M", po::value<string>(), "path to Modelica library")
("results-file,r", po::value<string>(),"name of results file")
Expand Down Expand Up @@ -106,6 +107,7 @@ SimSettings OMCFactory::readSimulationParameter(int argc, const char* argv[])
double starttime = vm["start-time"].as<double>();
double stoptime = vm["stop-time"].as<double>();
double stepsize =vm["step-size"].as<double>();
bool nlsContinueOnError = vm["nls_continue"].as<bool>();

if (!(stepsize > 0.0))
stepsize = (stoptime - starttime) / vm["number-of-intervals"].as<int>();
Expand Down Expand Up @@ -199,7 +201,7 @@ SimSettings OMCFactory::readSimulationParameter(int argc, const char* argv[])



SimSettings settings = {solver,linSolver,nonLinSolver,starttime,stoptime,stepsize,1e-24,0.01,tolerance,resultsfilename,time_out,outputPointType,logSet};
SimSettings settings = {solver,linSolver,nonLinSolver,starttime,stoptime,stepsize,1e-24,0.01,tolerance,resultsfilename,time_out,outputPointType,logSet,nlsContinueOnError};


_library_path = libraries_path;
Expand Down
37 changes: 24 additions & 13 deletions SimulationRuntime/cpp/Solver/Hybrj/HybrjSettings.cpp
Expand Up @@ -9,50 +9,61 @@
#include <Solver/Hybrj/HybrjSettings.h>

HybrjSettings::HybrjSettings()
: iNewt_max (50)
, dRtol (1e-6)
, dAtol (1.0)
, dDelta (0.9)
: _iNewt_max (50)
, _dRtol (1e-6)
, _dAtol (1.0)
, _dDelta (0.9)
, _continueOnError(false)
{
};
/*max. Anzahl an Newtonititerationen pro Schritt (default: 25)*/
long int HybrjSettings::getNewtMax()
{
return iNewt_max;
return _iNewt_max;
}
void HybrjSettings::setNewtMax(long int max)
{
iNewt_max =max;
_iNewt_max =max;
}
/* Relative Toleranz für die Newtoniteration (default: 1e-6)*/
double HybrjSettings::getRtol()
{
return dRtol;
return _dRtol;
}
void HybrjSettings::setRtol(double t)
{
dRtol=t;
_dRtol=t;
}
/*Absolute Toleranz für die Newtoniteration (default: 1e-6)*/
double HybrjSettings::getAtol()
{
return dAtol;
return _dAtol;
}
void HybrjSettings::setAtol(double t)
{
dAtol =t;
_dAtol =t;
}
/*Dämpfungsfaktor (default: 0.9)*/
double HybrjSettings::getDelta()
{
return dDelta;
return _dDelta;
}
void HybrjSettings::setDelta(double t)
{
dDelta = t;
_dDelta = t;
}

void HybrjSettings::load(string)
{
}
/** @} */ // end of solverHybrj

void HybrjSettings::setContinueOnError(bool value)
{
_continueOnError = value;
}

bool HybrjSettings::getContinueOnError()
{
return _continueOnError;
}
/** @} */ // end of solverHybrj
27 changes: 13 additions & 14 deletions SimulationRuntime/cpp/Solver/Kinsol/Kinsol.cpp
Expand Up @@ -9,13 +9,8 @@
#include <nvector/nvector_serial.h>
#include <kinsol/kinsol.h>

#ifdef USE_SUNDIALS_LAPACK
#include <kinsol/kinsol_lapack.h>
#else
#include <kinsol/kinsol_spgmr.h>
#include <kinsol/kinsol_dense.h>
#endif //USE_SUNDIALS_LAPACK

#include <kinsol/kinsol_spgmr.h>
#include <kinsol/kinsol_dense.h>

#include <kinsol/kinsol_spbcgs.h>
#include <kinsol/kinsol_sptfqmr.h>
Expand Down Expand Up @@ -255,11 +250,7 @@ void Kinsol::initialize()
if (check_flag(&idid, (char *)"KINSetUserData", 1))
throw ModelicaSimulationError(ALGLOOP_SOLVER,"Kinsol::initialize()");

#ifdef USE_SUNDIALS_LAPACK
KINLapackDense(_kinMem, _dimSys);
#else
KINDense(_kinMem, _dimSys);
#endif //USE_SUNDIALS_LAPACK

/*will be used with new sundials version
if(_algLoop->isLinearTearing())
Expand Down Expand Up @@ -558,11 +549,19 @@ void Kinsol::solve()
return;
}

if(_iterationStatus == SOLVERERROR && !_eventRetry)
if(_iterationStatus == SOLVERERROR && !_eventRetry)
{
if(_kinsolSettings->getContinueOnError())
{
if(!_solverErrorNotificationGiven)
{
throw ModelicaSimulationError(ALGLOOP_SOLVER,"Nonlinear solver failed!");
Logger::write("Kinsol: Solver error detected. The simulation will continue, but the results may be incorrect.",LC_NLS,LL_WARNING);
_solverErrorNotificationGiven = true;
}

}
else
throw ModelicaSimulationError(ALGLOOP_SOLVER,"Nonlinear solver failed!");
}
}
}

Expand Down
13 changes: 9 additions & 4 deletions SimulationRuntime/cpp/Solver/Kinsol/KinsolLapack.cpp
Expand Up @@ -104,18 +104,23 @@ static void KINLapackCompletePivotingFree(KINMem kin_mem)
* \return Return_Description
* \details Details
*/
static int KINLapackCompletePivotingSolve(KINMem kin_mem, N_Vector x, N_Vector b,realtype *res_norm/*used in new sundials: realtype *sJpnorm, realtype *sFdotJp*/)
#if SUNDIALS_MAJOR_VERSION > 2 || (SUNDIALS_MAJOR_VERSION == 2 && SUNDIALS_MINOR_VERSION > 5)
static int KINLapackCompletePivotingSolve(KINMem kin_mem, N_Vector x, N_Vector b, realtype *sJpnorm, realtype *sFdotJp)
#else
static int KINLapackCompletePivotingSolve(KINMem kin_mem, N_Vector x, N_Vector b, realtype *res_norm)
#endif
{
int flag=0;
linSysData* data = (linSysData*)kin_mem->kin_lmem;

dgesc2_(&data->n,data->jac,&data->n, NV_DATA_S(b), data->ihelpArray, data->jhelpArray, data->scale);


memcpy(NV_DATA_S(x),NV_DATA_S(b),data->n*sizeof(double));
#if SUNDIALS_MAJOR_VERSION > 2 || (SUNDIALS_MAJOR_VERSION == 2 && SUNDIALS_MINOR_VERSION > 5)
*sFdotJp = N_VDotProd(kin_mem->kin_fval, b);
#else
N_VProd(b, kin_mem->kin_fscale, b);

//*used in new sundials: *sFdotJp = N_VDotProd(kin_mem->kin_fval, b);
#endif

return(flag);
}
Expand Down

0 comments on commit 12f2912

Please sign in to comment.