Skip to content

Commit

Permalink
Experimental rounding error fixing for fluid models in MSL
Browse files Browse the repository at this point in the history
  • Loading branch information
@OpenModelica-Hudson
qichenghua authored and OpenModelica-Hudson committed Jul 18, 2017
1 parent faa82e6 commit 6beb4f4
Show file tree
Hide file tree
Showing 2 changed files with 63 additions and 44 deletions.
5 changes: 4 additions & 1 deletion SimulationRuntime/cpp/Include/Solver/Nox/Nox.h
View file
Expand Up @@ -53,7 +53,10 @@ class Nox : public IAlgLoopSolver
*_y0,
*_y_old,
*_y_new,
*_yScale;
*_yScale,
_locTol,
_currentIterateNorm,
*_currentIterate;

Teuchos::RCP<NoxLapackInterface> _noxLapackInterface;

Expand Down
102 changes: 59 additions & 43 deletions SimulationRuntime/cpp/Solver/Nox/Nox.cpp
View file
Expand Up @@ -30,6 +30,7 @@ Nox::Nox(INonLinearAlgLoop* algLoop, INonLinSolverSettings* settings)
, _firstCall (true)
, _generateoutput (false)
, _useDomainScaling (false)
, _currentIterate (NULL)
{
_dimSys=_algLoop->getDimReal();
}
Expand All @@ -42,6 +43,7 @@ Nox::~Nox()
if(_y_new) delete [] _y_new;
if(_yScale) delete [] _yScale;
//if(_noxLapackInterface) delete _noxLapackInterface;
if(_currentIterate) delete [] _currentIterate;

}

Expand All @@ -64,6 +66,7 @@ void Nox::initialize()
_y_old = new double[_dimSys];
_y_new = new double[_dimSys];
_yScale = new double[_dimSys];
_currentIterate = new double[_dimSys];

_algLoop->getReal(_y);
_algLoop->getReal(_y0);
Expand Down Expand Up @@ -92,10 +95,14 @@ void Nox::solve()
double * rhs = new double[_dimSys];
double sum;

_locTol=5.0e-7;
_currentIterateNorm=1.0e2;

_output = Teuchos::rcp(new std::stringstream);
if (_firstCall) initialize();

// Set up the status tests
// can be moved to initialize.
_statusTestNormF = Teuchos::rcp(new NOX::StatusTest::NormF(1.0e-13));
_statusTestMaxIters = Teuchos::rcp(new NOX::StatusTest::MaxIters(100));
_statusTestStagnation = Teuchos::rcp(new NOX::StatusTest::Stagnation(15,0.99));
Expand Down Expand Up @@ -272,8 +279,6 @@ void Nox::solve()
}
}

int numberofdifferentnormtests=5;//this variable is used twice.

_algLoop->setReal(_y);
try{
_algLoop->evaluate();
Expand All @@ -300,53 +305,55 @@ void Nox::solve()
if (_generateoutput) std::cout << "simtime=" << _algLoop->getSimTime() << std::endl;
_iterationStatus=DONE;
}else{
//skip iter=1,2,4,6,7,9,11,12,13

if (iter==11){
iter = 14;
}
if (iter==9){
iter = 10;
}
if (iter==6){
iter = 8;
}
if (iter==4){
iter = 5;
}
if (iter==0){
iter = 3;
}
if ((iter == 1) || (iter == 2) || (iter == 4) || (iter == 6) || (iter == 7) || (iter == 9)){
std::cout << "You failed heavily at iter = " << iter << " and simtime " << _algLoop->getSimTime() << std::endl;
if(true){
try{
_algLoop->setReal(_y);
_algLoop->evaluate();
_algLoop->getRHS(rhs);

double*ypluseps=new double [_dimSys];
double*rhs2=new double [_dimSys];
memcpy(ypluseps,_y,_dimSys*sizeof(double));
for(int i=0;i<_dimSys;i++){
ypluseps[i]*=(1.0+std::numeric_limits<double>::epsilon()*8.0);//increase _ypluseps in the smallest amount possible.
}
_algLoop->setReal(ypluseps);
_algLoop->evaluate();
_algLoop->getRHS(rhs2);
double threshold=0.0;
for(int i=0;i<_dimSys;i++){
threshold+=(rhs2[i]-rhs[i])*(rhs2[i]-rhs[i]);
}
threshold=std::sqrt(threshold);

delete [] ypluseps;
delete [] rhs2;

if(_statusTestNormF->getNormF()<2*threshold){
std::cout << "threshold (this should typically be a very small value (we vary the real value by ). In case of occurence of rounding errors, the threshold should be big though.) = " << threshold << ", simtime = " << _algLoop->getSimTime() << ", algloop " << _algLoop->getEquationIndex() << std::endl;
_algLoop->setReal(_y);
_algLoop->evaluate();
_algLoop->getRHS(rhs);
_iterationStatus=DONE;
break;
}
}
catch(const std::exception &ex)
{
if (_generateoutput) std::cout << "algloop evaluation after solve failed with error message:" << std::endl << ex.what() << std::endl << "Trying to continue without. This should hopefully lead to statusTest::Failed." << std::endl;
}
}

switch(iter%(numberofdifferentnormtests)){
case 0:
*_statusTestNormF=NOX::StatusTest::NormF(1.0e-13);
break;
case 1:
*_statusTestNormF=NOX::StatusTest::NormF(1.0e-11);
break;
case 2:
*_statusTestNormF=NOX::StatusTest::NormF(1.0e-9);
break;
case 3:
*_statusTestNormF=NOX::StatusTest::NormF(1.0e-7);
break;
default:
*_statusTestNormF=NOX::StatusTest::NormF(1.0e-5);
break;
}
_statusTestsCombo = Teuchos::rcp(new NOX::StatusTest::Combo(NOX::StatusTest::Combo::OR, _statusTestNormF, _statusTestMaxIters));
_statusTestsCombo->addStatusTest(_statusTestStagnation);
_statusTestsCombo->addStatusTest(_statusTestDivergence);
// if(_statusTestNormF->getNormF()<std::min(_locTol,_currentIterateNorm)){
// memcpy(_currentIterate,_y,_dimSys*sizeof(double));
// _currentIterateNorm=_statusTestNormF->getNormF();
// }

switch (iter/numberofdifferentnormtests){
switch (iter){
case 0:
break;
//default Nox Line Search failed -> Try damped Newton instead
case 1:
case 1://should delete this. It typically does not help (verify this))
_solverParametersPtr->sublist("Line Search").sublist("Full Step").set("Full Step", 0.5);

//default Nox Line Search failed -> Try Backtracking instead
Expand Down Expand Up @@ -415,6 +422,14 @@ void Nox::solve()

//Tensor failed or other failure
default:
if(_currentIterateNorm<_locTol){
memcpy(_y,_currentIterate,_dimSys*sizeof(double));
_iterationStatus=DONE;
break;
}
if(_currentIterateNorm<_locTol){
std::cout << "you should not see this." << std::endl;
}
int numberofhomotopytries = 0;
while((_iterationStatus!=DONE)){
//todo: This is implemented in the worst way possible. Please fix this.
Expand All @@ -430,6 +445,7 @@ void Nox::solve()
}
numberofhomotopytries++;
}

break;
}
//comment in for debugging
Expand Down

0 comments on commit 6beb4f4

Please sign in to comment.