fixed issue that occured if evaluating the algloop at the initial guess throws a division by zero error and added homotopy solver

Author: qichenghua

SimulationRuntime/cpp/Include/Solver/Nox/Nox.h

@@ -30,6 +30,9 @@ class Nox : public IAlgLoopSolver
   void solverinit();
   void createStatusTests();
   void createSolverParameters();
+  void LocaHomotopySolve(int numberofhomotopytries);
+  NOX::StatusTest::StatusType BasicNLSsolve();
+  NOX::StatusTest::StatusType secondBasicNLSsolve();
 
   //void check4EventRetry(double* y)
 
@@ -60,6 +63,9 @@ class Nox : public IAlgLoopSolver
   //used for status tests
   Teuchos::RCP<NOX::StatusTest::NormF> _statusTestNormF;
   Teuchos::RCP<NOX::StatusTest::MaxIters> _statusTestMaxIters;
+  Teuchos::RCP<NOX::StatusTest::Stagnation> _statusTestStagnation;
+  Teuchos::RCP<NOX::StatusTest::Divergence> _statusTestDivergence;
+
   Teuchos::RCP<NOX::StatusTest::Combo> _statusTestsCombo;
 
   //list of solver parameters
@@ -71,5 +77,6 @@ class Nox : public IAlgLoopSolver
   bool _firstCall;
   bool _generateoutput;
   bool _useDomainScaling;
+
 };
 /** @} */ // end of solverNox

SimulationRuntime/cpp/Include/Solver/Nox/NoxLapackInterface.h

@@ -1,11 +1,13 @@
-#include "NOX_LAPACK_Group.H"
+//somehow this is not good...
+#include <LOCA.H>
+#include <LOCA_LAPACK.H>
 
-class NoxLapackInterface : public NOX::LAPACK::Interface {
+class NoxLapackInterface : public LOCA::LAPACK::Interface {
 
   public:
 
     //! Constructor
-    NoxLapackInterface(INonLinearAlgLoop *algLoop);
+    NoxLapackInterface(INonLinearAlgLoop *algLoop, int numberofhomotopytries);
 
     //! Destructor
     ~NoxLapackInterface();
@@ -16,7 +18,19 @@ class NoxLapackInterface : public NOX::LAPACK::Interface {
 
     bool computeJacobian(NOX::LAPACK::Matrix<double>& J, const NOX::LAPACK::Vector & x);
 
+	//! Sets parameters
+	void setParams(const LOCA::ParameterVector& p);
+
+	//! Prints solution after successful step
+	void printSolution(const NOX::LAPACK::Vector &x, const double conParam);
+
   private:
+
+	bool computeSimplifiedF(NOX::LAPACK::Vector& f, const NOX::LAPACK::Vector &x);
+	bool computeActualF(NOX::LAPACK::Vector& f, const NOX::LAPACK::Vector &x);
+	NOX::LAPACK::Vector applyMatrixtoVector(const NOX::LAPACK::Matrix<double> &A, const NOX::LAPACK::Vector &x);
+	void checkdimensionof(const NOX::LAPACK::Vector &x);
+
     //! Initial guess
     Teuchos::RCP<NOX::LAPACK::Vector> _initialGuess;
 	INonLinearAlgLoop *_algLoop;///< Algebraic loop to be solved, required to obtain value of f
@@ -25,4 +39,9 @@ class NoxLapackInterface : public NOX::LAPACK::Interface {
 	bool _generateoutput;
 	bool _useDomainScaling;
 	bool _useFunctionValueScaling;
+	double _lambda;//homotopy parameter
+	bool _computedinitialguess;
+	int _numberofhomotopytries;
+	bool _evaluatedJacobianAtInitialGuess;
+	Teuchos::RCP<NOX::LAPACK::Matrix<double>> _J;//F'(x_0)
 };

SimulationRuntime/cpp/Solver/Kinsol/Kinsol.cpp

@@ -718,6 +718,7 @@ void Kinsol::solveNLS()
 		idid = KINSol(_kinMem, _Kin_y, method, _Kin_yScale, _Kin_fScale);
 
 		KINGetFuncNorm(_kinMem, &_fnorm);
+		if(!_fValid && (idid==KIN_FIRST_SYSFUNC_ERR)) throw ModelicaSimulationError(ALGLOOP_SOLVER,"Algloop could not be evaluated! Evaluation failed at the first call.");
 		//if(_fnorm/euclidNorm(_dimSys,_yScale) < 1e-4)
 		if(idid != KIN_SUCCESS &&  _fnorm < locTol && _fnorm < _currentIterateNorm)
 		{