Skip to content

Commit

Permalink
added new abort criterion for nonlinear solver based on sign changes …
Browse files Browse the repository at this point in the history
…in case of errors due to lacking precision of double precision variables implemented in commit b433b2d6ec92460fae06a9dd1b5243ee438cf0a7. Now this criterion is applied after each iteration step and not after a complete iteration cycle.
  • Loading branch information
qichenghua authored and OpenModelica-Hudson committed Jul 18, 2017
1 parent 267a1c7 commit 5914059
Show file tree
Hide file tree
Showing 6 changed files with 157 additions and 4 deletions.
@@ -0,0 +1,52 @@
#ifndef NOX_STATUSTEST_SGNCHANGE_H
#define NOX_STATUSTEST_SGNCHANGE_H

#include "NOX_Solver_Generic.H"
#include "NOX_StatusTest_Generic.H" // base class
#include "NOX_LAPACK_Vector.H"

namespace NOX{
namespace StatusTest{

class SgnChange : public Generic {

public:

//! Constructor
SgnChange() = delete;//avoid default constructor, only available in C++-11
SgnChange(double tol);

//! Destructor
~SgnChange();
// derived
virtual StatusType checkStatus(const NOX::Solver::Generic& problem, NOX::StatusTest::CheckType checkType);
// derived
virtual StatusType getStatus() const;
// derived
virtual std::ostream& print(std::ostream& stream, int indent = 0) const;

private:
//! %do initialization in checkStatus only once
void initialize(const NOX::Solver::Generic& problem);

//! %Status
StatusType _status;

//! %Tolerance
double _tol;

//! %Number of sign changes
int _numSignChanges;

//! %Only initialize once
bool _firstCall;

//! %Dimension of solution and function
unsigned int _dimSys;

//! %solution/function vector pointers
Teuchos::RCP<NOX::LAPACK::Vector> _x0,_x1,_f0,_f1;
};
}//namespace StatusTest
}//namespace NOX
#endif
3 changes: 3 additions & 0 deletions SimulationRuntime/cpp/Include/Solver/Nox/Nox.h
Expand Up @@ -5,6 +5,8 @@
*/


#include <Solver/Nox/NOX_StatusTest_SgnChange.H>


class Nox : public IAlgLoopSolver
{
Expand Down Expand Up @@ -75,6 +77,7 @@ class Nox : public IAlgLoopSolver
Teuchos::RCP<NOX::StatusTest::MaxIters> _statusTestMaxIters;
Teuchos::RCP<NOX::StatusTest::Stagnation> _statusTestStagnation;
Teuchos::RCP<NOX::StatusTest::Divergence> _statusTestDivergence;
Teuchos::RCP<NOX::StatusTest::SgnChange> _statusTestSgnChange;

Teuchos::RCP<NOX::StatusTest::Combo> _statusTestsCombo;

Expand Down
3 changes: 2 additions & 1 deletion SimulationRuntime/cpp/Solver/Nox/CMakeLists.txt
Expand Up @@ -2,7 +2,7 @@ cmake_minimum_required(VERSION 2.8.9)

project(${NoxName})

add_library(${NoxName} Nox.cpp NoxSettings.cpp FactoryExport.cpp NoxLapackInterface.cpp)
add_library(${NoxName} Nox.cpp NoxSettings.cpp FactoryExport.cpp NoxLapackInterface.cpp NOX_StatusTest_SgnChange.C)

if(NOT BUILD_SHARED_LIBS)
set_target_properties(${NoxName} PROPERTIES COMPILE_DEFINITIONS "RUNTIME_STATIC_LINKING;ENABLE_SUNDIALS_STATIC")
Expand All @@ -18,4 +18,5 @@ install(FILES
${CMAKE_SOURCE_DIR}/Include/Solver/Nox/NoxSettings.h
${CMAKE_SOURCE_DIR}/Include/Solver/Nox/FactoryExport.h
${CMAKE_SOURCE_DIR}/Include/Solver/Nox/NoxLapackInterface.h
${CMAKE_SOURCE_DIR}/Include/Solver/Nox/NOX_StatusTest_SgnChange.H
DESTINATION include/omc/cpp/Solver/Nox)
93 changes: 93 additions & 0 deletions SimulationRuntime/cpp/Solver/Nox/NOX_StatusTest_SgnChange.C
@@ -0,0 +1,93 @@
#include <Solver/Nox/NOX_StatusTest_SgnChange.H>
#include "NOX_Common.H"
#include "NOX_LAPACK_Group.H"
#include "NOX_Abstract_Vector.H"
#include "NOX_Abstract_Group.H"

NOX::StatusTest::SgnChange::SgnChange(double tol) :
_tol(tol),
_status(NOX::StatusTest::Unevaluated),
_numSignChanges(0),
_firstCall(true),
_dimSys(0)
{

}

NOX::StatusTest::SgnChange::~SgnChange()
{

}

void NOX::StatusTest::SgnChange::initialize(const NOX::Solver::Generic& problem)
{
_dimSys = problem.getSolutionGroup().getX().length();
_x0=Teuchos::rcp(new NOX::LAPACK::Vector(_dimSys));
_x1=Teuchos::rcp(new NOX::LAPACK::Vector(_dimSys));
_f0=Teuchos::rcp(new NOX::LAPACK::Vector(_dimSys));
_f1=Teuchos::rcp(new NOX::LAPACK::Vector(_dimSys));
_firstCall=false;
}

NOX::StatusTest::StatusType NOX::StatusTest::SgnChange::checkStatus(const NOX::Solver::Generic& problem, NOX::StatusTest::CheckType checkType)
{
if(_firstCall) initialize(problem);
std::vector<bool> fSignChange(_dimSys, false);
Teuchos::RCP<NOX::LAPACK::Group> grp=Teuchos::rcp(new NOX::LAPACK::Group(dynamic_cast<const NOX::LAPACK::Group&>(problem.getSolutionGroup())));//(problem.getSolutionGroup());//solutiongroup is constant, thus we need to assign it to be able to modify it.
*_x1=*_x0=grp->getX();
grp->computeF();
*_f0=grp->getF();

for (int i=0;i<_dimSys;i++){
// compute F at x1=x0+2*eps*e_i and save in f1
(*_x1)(i)=std::nextafter(std::nextafter((*_x1)(i),std::numeric_limits<double>::max()),std::numeric_limits<double>::max());
grp->setX(*_x1);
grp->computeF();
*_f1=grp->getF();
// compare
for(int j=0;j<_dimSys;j++){
if ((*_f0)(j)*(*_f1)(j)<=0.0){
fSignChange[j]= true;
}
}
(*_x1)(i)=(*_x0)(i);

// do the same for x0-2*eps*e_i
(*_x1)(i)=std::nextafter(std::nextafter((*_x1)(i),-std::numeric_limits<double>::max()),-std::numeric_limits<double>::max());
grp->setX(*_x1);
grp->computeF();
*_f1=grp->getF();
// compare
for(int j=0;j<_dimSys;j++){
if ((*_f0)(j)*(*_f1)(j)<=0.0){
fSignChange[j]= true;
}
}
(*_x1)(i)=(*_x0)(i);
}
_numSignChanges=std::count(fSignChange.begin(), fSignChange.end(), true);

grp->setX(*_x0);
grp->computeF();

//return converged, if all entries of fSignChange are true, unconverged otherwise.//only available in C++-11//alternative: use _numSignChanges==_dimSys as criteria instead.
_status = (std::all_of(fSignChange.begin(),fSignChange.end(),[](bool a){return a;})) ? NOX::StatusTest::Converged : NOX::StatusTest::Unconverged;
return _status;
}

NOX::StatusTest::StatusType NOX::StatusTest::SgnChange::getStatus() const
{
return _status;
}

std::ostream& NOX::StatusTest::SgnChange::print(std::ostream& stream, int indent) const
{
for (int j = 0; j < indent; j ++)
stream << ' ';
stream << _status;
stream << "Found sign changes in " << _numSignChanges;
stream << " out of " << _dimSys;
stream << " components.";
stream << std::endl;
return stream;
}
6 changes: 5 additions & 1 deletion SimulationRuntime/cpp/Solver/Nox/Nox.cpp
Expand Up @@ -89,10 +89,12 @@ void Nox::initialize()
_statusTestMaxIters = Teuchos::rcp(new NOX::StatusTest::MaxIters(100));
_statusTestStagnation = Teuchos::rcp(new NOX::StatusTest::Stagnation(15,0.99));
_statusTestDivergence = Teuchos::rcp(new NOX::StatusTest::Divergence(1.0e13));
_statusTestSgnChange = Teuchos::rcp(new NOX::StatusTest::SgnChange(5.0e-7));

_statusTestsCombo = Teuchos::rcp(new NOX::StatusTest::Combo(NOX::StatusTest::Combo::OR, _statusTestNormF, _statusTestMaxIters));
_statusTestsCombo->addStatusTest(_statusTestStagnation);
_statusTestsCombo->addStatusTest(_statusTestDivergence);
_statusTestsCombo->addStatusTest(_statusTestSgnChange);

if (_generateoutput) std::cout << "ending init" << std::endl;
}
Expand Down Expand Up @@ -122,7 +124,7 @@ void Nox::solve()
// Set the level of output
addPrintingList(_solverParametersPtr);
if (_generateoutput){
_solverParametersPtr->sublist("Printing").set("Output Information", NOX::Utils::Error + NOX::Utils::Warning + NOX::Utils::OuterIteration + NOX::Utils::Details + NOX::Utils::Debug); //(there are also more options, but error and outer iteration are the ones that I commonly use.
_solverParametersPtr->sublist("Printing").set("Output Information", NOX::Utils::Error + NOX::Utils::Warning + NOX::Utils::OuterIteration + NOX::Utils::InnerIteration + NOX::Utils::Details + NOX::Utils::Debug); //(there are also more options, but error and outer iteration are the ones that I commonly use.
}

//_solverParametersPtr->sublist("Direction").set("Method", "Steepest Descent");
Expand Down Expand Up @@ -575,6 +577,7 @@ void Nox::modifySolverParameters(const Teuchos::RCP<Teuchos::ParameterList> solv
//Polynomial failed -> Try More'-Thuente instead
case 4:
solverParametersPtr->sublist("Line Search").set("Method", "More'-Thuente");
// solverParametersPtr->sublist("Line Search").sublist("More'-Thuente").set("Recovery Step", std::numeric_limits<double>::min());//I would set this to 0.0, but then More Thuente throws an error:/
//solverParametersPtr->sublist("Line Search").sublist("More'-Thuente").set("Sufficient Decrease", 1.0e-2);
break;

Expand All @@ -590,6 +593,7 @@ void Nox::modifySolverParameters(const Teuchos::RCP<Teuchos::ParameterList> solv
case 6:
//solverParametersPtr->sublist("Trust Region").set("Use Dogleg Segment Minimization", true);
solverParametersPtr->set("Nonlinear Solver", "Inexact Trust Region Based");
solverParametersPtr->sublist("Trust Region").set("Recovery Step", 0.0);
break;

//Inexact Trust Region failed -> Try Tensor instead
Expand Down
4 changes: 2 additions & 2 deletions SimulationRuntime/cpp/Solver/Nox/NoxLapackInterface.cpp
Expand Up @@ -116,7 +116,7 @@ const NOX::LAPACK::Vector& NoxLapackInterface::getInitialGuess()
};

bool NoxLapackInterface::computeActualF(NOX::LAPACK::Vector& f, const NOX::LAPACK::Vector &x){
//if (_algLoop->getSimTime()>0.088421) _generateoutput=true;
// if (_algLoop->getSimTime()>0.0) _generateoutput=true;
checkdimensionof(x);

for (int i=0;i<_dimSys;i++){
Expand Down Expand Up @@ -212,7 +212,7 @@ bool NoxLapackInterface::computeActualF(NOX::LAPACK::Vector& f, const NOX::LAPAC
// std::cout << std::endl;
// std::cout << "or also" << std::endl;
// f.print(std::cout);
}
}//_generateoutput=false;
return true;
}

Expand Down

0 comments on commit 5914059

Please sign in to comment.