Skip to content

Commit d5d856e

Browse files
qichenghuaOpenModelica-Hudson
qichenghua
authored and
OpenModelica-Hudson
committed
refactored Nox
1 parent 557c017 commit d5d856e

File tree

2 files changed

+122
-113
lines changed
  • SimulationRuntime/cpp

2 files changed

+122
-113
lines changed

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

Lines changed: 4 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -33,11 +33,13 @@ class Nox : public IAlgLoopSolver
3333
void LocaHomotopySolve(const int numberofhomotopytries);
3434
NOX::StatusTest::StatusType BasicNLSsolve();
3535
void addPrintingList(const Teuchos::RCP<Teuchos::ParameterList> solverParametersPtr);
36-
void copySolution(const Teuchos::RCP<const NOX::Solver::Generic> solver,double* algLoopSolution);
36+
void copySolution(const Teuchos::RCP<const NOX::Solver::Generic> solver,double* const algLoopSolution);
3737
void printLogger();
3838
void divisionbyzerohandling(double const * const y0);
3939
bool checkwhethersolutionisnearby(double const * const y);
4040
bool isdivisionbyzeroerror(const std::exception &ex);
41+
void modifySolverParameters(const Teuchos::RCP<Teuchos::ParameterList> solverParametersPtr,const int iter);
42+
4143

4244
//void check4EventRetry(double* y)
4345

@@ -52,7 +54,7 @@ class Nox : public IAlgLoopSolver
5254
ITERATIONSTATUS
5355
_iterationStatus; ///< Output - Denotes the status of iteration
5456

55-
long int _dimSys;
57+
const long int _dimSys;
5658

5759
double
5860
*_y,

‎SimulationRuntime/cpp/Solver/Nox/Nox.cpp‎

Lines changed: 118 additions & 111 deletions
Original file line numberDiff line numberDiff line change
@@ -31,8 +31,9 @@ Nox::Nox(INonLinearAlgLoop* algLoop, INonLinSolverSettings* settings)
3131
, _generateoutput (false)
3232
, _useDomainScaling (false)
3333
, _currentIterate (NULL)
34+
, _dimSys (_algLoop->getDimReal())
3435
{
35-
_dimSys=_algLoop->getDimReal();
36+
3637
}
3738

3839
Nox::~Nox()
@@ -174,8 +175,8 @@ void Nox::solve()
174175
}
175176
catch(const std::exception &ex)
176177
{
177-
std::cout << std::endl << "sth went wrong during solver building, with error message" << ex.what() << std::endl;
178-
throw ModelicaSimulationError(ALGLOOP_SOLVER,"solving error");
178+
std::cout << std::endl << "sth went wrong during group building, with error message" << ex.what() << std::endl;
179+
throw ModelicaSimulationError(ALGLOOP_SOLVER,"solving error");//could be deleted
179180
}
180181

181182
if (_generateoutput) std::cout << "building solver" << std::endl;
@@ -276,113 +277,11 @@ void Nox::solve()
276277
// _currentIterateNorm=_statusTestNormF->getNormF();
277278
// }
278279

279-
if (iter==0) iter = 2;
280-
281-
if (iter==1) iter=2;//skip damped Newton.
282-
//we could skip case 2 as well
283-
284-
switch (iter){
285-
case 0:
286-
break;
287-
//default Nox Line Search failed -> Try damped Newton instead
288-
case 1://should delete this. It typically does not help (verify this))
289-
_solverParametersPtr->sublist("Line Search").sublist("Full Step").set("Full Step", 0.5);
290-
break;
291-
292-
//default Nox Line Search failed -> Try Backtracking instead
293-
case 2:
294-
_solverParametersPtr->sublist("Line Search").set("Method", "Backtrack");
295-
_solverParametersPtr->sublist("Line Search").sublist("Backtrack").set("Default Step", 1024.0*1024.0);
296-
_solverParametersPtr->sublist("Line Search").sublist("Backtrack").set("Minimum Step", 1.0e-30);
297-
_solverParametersPtr->sublist("Line Search").sublist("Backtrack").set("Recovery Step", 0.0);
298-
//std::cout << "we just set the solver parameters to use backtrack. the currently used options are:" << std::endl;
299-
//_solverParametersPtr->print();
300-
//std::cout << std::endl;
301-
break;
302-
303-
//Backtracking failed -> Try Polynomial with various parameters instead
304-
case 3:
305-
_solverParametersPtr->sublist("Line Search").set("Method", "Polynomial");
306-
_solverParametersPtr->sublist("Line Search").sublist("Polynomial").set("Minimum Step", 1.0e-30);
307-
//_solverParametersPtr->sublist("Line Search").sublist("Polynomial").set("Alpha Factor", 1.0e-2);
308-
break;
309-
310-
// case 4:
311-
// _solverParametersPtr->sublist("Line Search").sublist("Polynomial").set("Interpolation Type", "Quadratic");
312-
// break;
313-
314-
// case 5:
315-
// _solverParametersPtr->sublist("Line Search").sublist("Polynomial").set("Interpolation Type", "Quadratic3");
316-
// break;
317-
318-
// case 6:
319-
// _solverParametersPtr->sublist("Line Search").sublist("Polynomial").set("Sufficient Decrease Condition", "Ared/Pred");
320-
// _solverParametersPtr->sublist("Line Search").sublist("Polynomial").set("Interpolation Type", "Cubic");
321-
// break;
322-
323-
// case 7:
324-
// _solverParametersPtr->sublist("Line Search").sublist("Polynomial").set("Interpolation Type", "Quadratic");
325-
// break;
326-
327-
// case 8:
328-
// _solverParametersPtr->sublist("Line Search").sublist("Polynomial").set("Interpolation Type", "Quadratic3");
329-
// break;
330-
331-
//Polynomial failed -> Try More'-Thuente instead
332-
case 4:
333-
_solverParametersPtr->sublist("Line Search").set("Method", "More'-Thuente");
334-
//_solverParametersPtr->sublist("Line Search").sublist("More'-Thuente").set("Sufficient Decrease", 1.0e-2);
335-
break;
336-
337-
//More'-Thuente failed -> Try Trust Region instead
338-
case 5:
339-
//_solverParametersPtr->sublist("Line Search").remove("Method", true);
340-
_solverParametersPtr->set("Nonlinear Solver", "Trust Region Based");
341-
_solverParametersPtr->sublist("Trust Region").set("Minimum Improvement Ratio", 1.0e-4);
342-
_solverParametersPtr->sublist("Trust Region").set("Minimum Trust Region Radius", 1.0e-6);
343-
break;
344-
345-
//Trust Region failed -> Try Inexact Trust Region instead
346-
case 6:
347-
//_solverParametersPtr->sublist("Trust Region").set("Use Dogleg Segment Minimization", true);
348-
_solverParametersPtr->set("Nonlinear Solver", "Inexact Trust Region Based");
349-
break;
350-
351-
//Inexact Trust Region failed -> Try Tensor instead
352-
//case 11:
353-
//_solverParametersPtr->set("Nonlinear Solver", "Tensor Based");
354-
//break;
355-
356-
//Tensor failed or other failure
357-
default:
358-
// if(_currentIterateNorm<_locTol){
359-
// memcpy(_y,_currentIterate,_dimSys*sizeof(double));
360-
// _algLoop->setReal(_y);
361-
// _algLoop->evaluate();
362-
// _algLoop->getRHS(rhs);
363-
// _iterationStatus=DONE;
364-
// break;
365-
// }
366-
// if(_currentIterateNorm<_locTol){
367-
// std::cout << "you should not see this." << std::endl;
368-
// }
369-
int numberofhomotopytries = 0;
370-
while((_iterationStatus!=DONE)){
371-
//todo: This is implemented in the worst way possible. Please fix this.
372-
//try homotopy
373-
try{
374-
if ((_generateoutput)) std::cout << "solving with numberofhomotopytries = " << numberofhomotopytries << std::endl;
375-
LocaHomotopySolve(numberofhomotopytries);
376-
}
377-
catch(const std::exception &e){
378-
std::string errorstring(e.what());
379-
_iterationStatus=SOLVERERROR;
380-
throw ModelicaSimulationError(ALGLOOP_SOLVER,"Nonlinear solver nox failed with error message " +errorstring);
381-
}
382-
numberofhomotopytries++;
383-
}
384-
break;
385-
}
280+
if ((iter==0) || (iter==1) || (iter==2)) iter = 3;
281+
282+
283+
modifySolverParameters(_solverParametersPtr,iter);
284+
386285
//comment in for debugging
387286
if ((_generateoutput)){
388287
std::cout << "solutionvector=(";
@@ -624,6 +523,114 @@ NOX::StatusTest::StatusType Nox::BasicNLSsolve(){
624523
return status;
625524
}
626525

526+
void Nox::modifySolverParameters(const Teuchos::RCP<Teuchos::ParameterList> solverParametersPtr,const int iter){
527+
528+
switch (iter){
529+
case 0:
530+
break;
531+
//default Nox Line Search failed -> Try damped Newton instead
532+
case 1://should delete this. It typically does not help (verify this))
533+
solverParametersPtr->sublist("Line Search").sublist("Full Step").set("Full Step", 0.5);
534+
break;
535+
536+
//default Nox Line Search failed -> Try Backtracking instead
537+
case 2:
538+
solverParametersPtr->sublist("Line Search").set("Method", "Backtrack");
539+
solverParametersPtr->sublist("Line Search").sublist("Backtrack").set("Default Step", 1024.0*1024.0);
540+
solverParametersPtr->sublist("Line Search").sublist("Backtrack").set("Minimum Step", 1.0e-30);
541+
solverParametersPtr->sublist("Line Search").sublist("Backtrack").set("Recovery Step", 0.0);
542+
//std::cout << "we just set the solver parameters to use backtrack. the currently used options are:" << std::endl;
543+
//solverParametersPtr->print();
544+
//std::cout << std::endl;
545+
break;
546+
547+
//Backtracking failed -> Try Polynomial with various parameters instead
548+
case 3:
549+
solverParametersPtr->sublist("Line Search").set("Method", "Polynomial");
550+
solverParametersPtr->sublist("Line Search").sublist("Polynomial").set("Minimum Step", 1.0e-30);
551+
//solverParametersPtr->sublist("Line Search").sublist("Polynomial").set("Alpha Factor", 1.0e-2);
552+
break;
553+
554+
// case 4:
555+
// solverParametersPtr->sublist("Line Search").sublist("Polynomial").set("Interpolation Type", "Quadratic");
556+
// break;
557+
558+
// case 5:
559+
// solverParametersPtr->sublist("Line Search").sublist("Polynomial").set("Interpolation Type", "Quadratic3");
560+
// break;
561+
562+
// case 6:
563+
// solverParametersPtr->sublist("Line Search").sublist("Polynomial").set("Sufficient Decrease Condition", "Ared/Pred");
564+
// solverParametersPtr->sublist("Line Search").sublist("Polynomial").set("Interpolation Type", "Cubic");
565+
// break;
566+
567+
// case 7:
568+
// solverParametersPtr->sublist("Line Search").sublist("Polynomial").set("Interpolation Type", "Quadratic");
569+
// break;
570+
571+
// case 8:
572+
// solverParametersPtr->sublist("Line Search").sublist("Polynomial").set("Interpolation Type", "Quadratic3");
573+
// break;
574+
575+
//Polynomial failed -> Try More'-Thuente instead
576+
case 4:
577+
solverParametersPtr->sublist("Line Search").set("Method", "More'-Thuente");
578+
//solverParametersPtr->sublist("Line Search").sublist("More'-Thuente").set("Sufficient Decrease", 1.0e-2);
579+
break;
580+
581+
//More'-Thuente failed -> Try Trust Region instead
582+
case 5:
583+
//solverParametersPtr->sublist("Line Search").remove("Method", true);
584+
solverParametersPtr->set("Nonlinear Solver", "Trust Region Based");
585+
solverParametersPtr->sublist("Trust Region").set("Minimum Improvement Ratio", 1.0e-4);
586+
solverParametersPtr->sublist("Trust Region").set("Minimum Trust Region Radius", 1.0e-6);
587+
break;
588+
589+
//Trust Region failed -> Try Inexact Trust Region instead
590+
case 6:
591+
//solverParametersPtr->sublist("Trust Region").set("Use Dogleg Segment Minimization", true);
592+
solverParametersPtr->set("Nonlinear Solver", "Inexact Trust Region Based");
593+
break;
594+
595+
//Inexact Trust Region failed -> Try Tensor instead
596+
//case 11:
597+
//solverParametersPtr->set("Nonlinear Solver", "Tensor Based");
598+
//break;
599+
600+
//Tensor failed or other failure
601+
default:
602+
// if(_currentIterateNorm<_locTol){
603+
// memcpy(_y,_currentIterate,_dimSys*sizeof(double));
604+
// _algLoop->setReal(_y);
605+
// _algLoop->evaluate();
606+
// _algLoop->getRHS(rhs);
607+
// _iterationStatus=DONE;
608+
// break;
609+
// }
610+
// if(_currentIterateNorm<_locTol){
611+
// std::cout << "you should not see this." << std::endl;
612+
// }
613+
int numberofhomotopytries = 0;
614+
while((_iterationStatus!=DONE)){
615+
//todo: This is implemented in the worst way possible. Please fix this.
616+
//try homotopy
617+
try{
618+
if ((_generateoutput)) std::cout << "solving with numberofhomotopytries = " << numberofhomotopytries << std::endl;
619+
LocaHomotopySolve(numberofhomotopytries);
620+
}
621+
catch(const std::exception &e){
622+
std::string errorstring(e.what());
623+
_iterationStatus=SOLVERERROR;
624+
throw ModelicaSimulationError(ALGLOOP_SOLVER,"Nonlinear solver nox failed with error message " +errorstring);
625+
}
626+
numberofhomotopytries++;
627+
}
628+
break;
629+
}
630+
}
631+
632+
633+
627634
bool Nox::checkwhethersolutionisnearby(double const * const y){
628635
double* rhs=new double [_dimSys];
629636
double* rhs2=new double [_dimSys];
@@ -768,7 +775,7 @@ void Nox::addPrintingList(const Teuchos::RCP<Teuchos::ParameterList> solverParam
768775
solverParametersPtr->sublist("Printing").set("Error Stream", _output);
769776
}
770777

771-
void Nox::copySolution(const Teuchos::RCP<const NOX::Solver::Generic> solver,double* algLoopSolution){
778+
void Nox::copySolution(const Teuchos::RCP<const NOX::Solver::Generic> solver,double* const algLoopSolution){
772779
const NOX::LAPACK::Vector& NoxSolution = dynamic_cast<const NOX::LAPACK::Vector&>((dynamic_cast<const NOX::LAPACK::Group&>(solver->getSolutionGroup())).getX());
773780
for (int i=0;i<_dimSys;i++){
774781
if (_useDomainScaling){

0 commit comments

Comments
 (0)