Skip to content
This repository was archived by the owner on May 18, 2019. It is now read-only.

Commit 456875e

Browse files
qichenghuaOpenModelica-Hudson
qichenghua
authored and
OpenModelica-Hudson
committed
Further Refactoring of nonlinear solver nox
1 parent 8ea3763 commit 456875e

File tree

4 files changed

+78
-117
lines changed

4 files changed

+78
-117
lines changed

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

Lines changed: 3 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -38,10 +38,10 @@ class Nox : public IAlgLoopSolver
3838
void copySolution(const Teuchos::RCP<const NOX::Solver::Generic> solver,double* const algLoopSolution);
3939
void printLogger();
4040
void divisionbyzerohandling(double const * const y0);
41-
bool checkwhethersolutionisnearby(double const * const y);
41+
bool CheckWhetherSolutionIsNearby(double const * const y);
4242
bool isdivisionbyzeroerror(const std::exception &ex);
4343
void modifySolverParameters(const Teuchos::RCP<Teuchos::ParameterList> solverParametersPtr,const int iter);
44-
44+
Teuchos::RCP<Teuchos::ParameterList> setLocaParams();
4545
bool modify_y(const int counter);
4646
void BinRep(std::vector<double> &result, const int number);
4747

@@ -96,5 +96,6 @@ class Nox : public IAlgLoopSolver
9696
bool _useDomainScaling;
9797

9898
bool _OutOfProperMethods;
99+
LogCategory _lc;
99100
};
100101
/** @} */ // end of solverNox

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

Lines changed: 2 additions & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -31,10 +31,9 @@ class NoxLapackInterface : public LOCA::LAPACK::Interface {
3131
bool computeSimplifiedF(NOX::LAPACK::Vector& f, const NOX::LAPACK::Vector &x);
3232
bool computeActualF(NOX::LAPACK::Vector& f, const NOX::LAPACK::Vector &x);
3333
NOX::LAPACK::Vector applyMatrixtoVector(const NOX::LAPACK::Matrix<double> &A, const NOX::LAPACK::Vector &x);
34-
void checkdimensionof(const NOX::LAPACK::Vector &x);
3534

36-
//! Initial guess
37-
Teuchos::RCP<NOX::LAPACK::Vector> _initialGuess;
35+
//! Initial guess
36+
Teuchos::RCP<NOX::LAPACK::Vector> _initialGuess;
3837
INonLinearAlgLoop *_algLoop;///< Algebraic loop to be solved, required to obtain value of f
3938
double *_yScale, *_fScale, *_hugeabsolutevalues, *_xtemp, *_rhs;
4039
int _dimSys;

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

Lines changed: 73 additions & 99 deletions
Original file line numberDiff line numberDiff line change
@@ -28,10 +28,11 @@ Nox::Nox(INonLinearAlgLoop* algLoop, INonLinSolverSettings* settings)
2828
, _y_new (NULL)
2929
, _yScale (NULL)
3030
, _firstCall (true)
31-
, _generateoutput (true)
31+
, _generateoutput (false)
3232
, _useDomainScaling (false)
3333
, _currentIterate (NULL)
3434
, _dimSys (_algLoop->getDimReal())
35+
, _lc(LC_NLS)
3536
{
3637

3738
}
@@ -100,12 +101,13 @@ void Nox::initialize()
100101

101102
void Nox::solve()
102103
{
103-
104104
int iter=-1; //Iterationcount of proper methods
105105
NOX::StatusTest::StatusType status;
106106
bool handleddivisionbyzeroerror=false;
107107
_OutOfProperMethods=false;
108108

109+
LOGGER_WRITE_BEGIN("Nox: start solving " + to_string(_algLoop->getEquationIndex()) + " at Simulation time " + to_string(_algLoop->getSimTime()), _lc, LL_DEBUG);
110+
109111
//setup solver
110112
if (_firstCall)
111113
initialize();
@@ -140,7 +142,7 @@ void Nox::solve()
140142

141143
if (_generateoutput) std::cout << "handleddivisionbyzeroerror" << std::endl;
142144

143-
//try various methods, including homotopy
145+
//try various methods, excluding homotopy
144146

145147
while((_iterationStatus==CONTINUE) && (!_OutOfProperMethods)){
146148
iter++;
@@ -151,13 +153,7 @@ void Nox::solve()
151153
}
152154
_algLoop->setReal(_y);
153155

154-
try{
155-
_grp = Teuchos::rcp(new NOX::LAPACK::Group(*_noxLapackInterface));//this also calls the getInitialGuess-function in the NoxLapackInterface and sets the initial guess in the NOX::LAPACK::Group
156-
}
157-
catch(const std::exception &ex)
158-
{
159-
std::cout << "Here is a significant error. grp building should never throw any error. (except for maybe bad alloc or alike." << std::endl;
160-
}
156+
_grp = Teuchos::rcp(new NOX::LAPACK::Group(*_noxLapackInterface));//this also calls the getInitialGuess-function in the NoxLapackInterface and sets the initial guess in the NOX::LAPACK::Group
161157

162158
_solver = NOX::Solver::buildSolver(_grp, _statusTestsCombo, _solverParametersPtr);
163159

@@ -174,24 +170,15 @@ void Nox::solve()
174170

175171
if (_generateoutput) std::cout << "solving done" << std::endl;
176172

177-
// Get the answer
178-
copySolution(_solver,_y);
179-
180-
181-
_algLoop->setReal(_y);
182-
try{
183-
_algLoop->evaluate();
184-
}
185-
catch(const std::exception &ex)
186-
{
187-
if (status == NOX::StatusTest::Converged) std::cout << "Here is a significant error. The right hand side should be large and in particular nonzero (calling evaluate yields message " << ex.what() << "), however NOX reports success." << std::endl;
188-
}
189-
190173
if (status==NOX::StatusTest::Converged){
174+
// Get the answer
175+
copySolution(_solver,_y);
176+
_algLoop->setReal(_y);
177+
_algLoop->evaluate();
191178
_iterationStatus=DONE;
192179
}else{
193180
try{
194-
if(checkwhethersolutionisnearby(_y)){
181+
if(CheckWhetherSolutionIsNearby(_y)){
195182
_algLoop->setReal(_y);
196183
_algLoop->evaluate();
197184
_iterationStatus=DONE;
@@ -209,8 +196,8 @@ void Nox::solve()
209196
}
210197
}
211198

199+
//try homotopy
212200
int numberofhomotopytries = 0;
213-
214201
while((_iterationStatus==CONTINUE) && (numberofhomotopytries<_noxLapackInterface->getMaxNumberOfHomotopyTries())){
215202
try{
216203
LocaHomotopySolve(numberofhomotopytries);
@@ -221,7 +208,14 @@ void Nox::solve()
221208
}
222209
numberofhomotopytries++;
223210
}
224-
if (!_iterationStatus==DONE){
211+
212+
if (_iterationStatus==DONE){
213+
LOGGER_WRITE("Solution found!", _lc, LL_DEBUG);
214+
LOGGER_WRITE_VECTOR("Solution", _y, _dimSys, _lc, LL_DEBUG);
215+
LOGGER_WRITE_END(_lc, LL_DEBUG);
216+
}else{
217+
LOGGER_WRITE("Nox Failed!", _lc, LL_DEBUG);
218+
LOGGER_WRITE_END(_lc, LL_DEBUG);
225219
_iterationStatus=SOLVERERROR;
226220
throw ModelicaSimulationError(ALGLOOP_SOLVER,"Nonlinear solver failed!");
227221
}
@@ -326,7 +320,6 @@ void Nox::LocaHomotopySolve(const int numberofhomotopytries)
326320
// Create global data object
327321
Teuchos::RCP<LOCA::GlobalData> globalData = LOCA::createGlobalData(paramList, lapackFactory);
328322

329-
330323
// Set up the problem interface
331324
Teuchos::RCP<NoxLapackInterface> LocaLapackInterface = Teuchos::rcp(new NoxLapackInterface (_algLoop));//this also gets the nominal values
332325
LocaLapackInterface->setNumberOfHomotopyTries(numberofhomotopytries);
@@ -342,8 +335,6 @@ void Nox::LocaHomotopySolve(const int numberofhomotopytries)
342335
// Create the stepper
343336
LOCA::Stepper stepper(globalData, grp, _statusTestsCombo, paramList);
344337

345-
_iterationStatus=CONTINUE;
346-
347338
try{
348339
// Perform continuation run
349340
LOCA::Abstract::Iterator::IteratorStatus status = stepper.run();
@@ -366,57 +357,54 @@ void Nox::LocaHomotopySolve(const int numberofhomotopytries)
366357
if (_generateoutput) std::cout << "finishing solve!" << std::endl;
367358

368359
if(_iterationStatus==DONE){
369-
// Get the final solution from the stepper
360+
// Get the final solution from the stepper
370361
Teuchos::RCP<const LOCA::LAPACK::Group> finalGroup = Teuchos::rcp_dynamic_cast<const LOCA::LAPACK::Group>(stepper.getSolutionGroup());
371362
const NOX::LAPACK::Vector& finalSolution = dynamic_cast<const NOX::LAPACK::Vector&>(finalGroup->getX());
372363

373-
if(_generateoutput){
374-
// Print final solution
375-
std::cout << std::endl << "Final solution is " << std::endl;
376-
finalGroup->printSolution(finalSolution, finalGroup->getParam("lambda"));//does not work
377-
}
364+
if(_generateoutput){
365+
// Print final solution
366+
std::cout << std::endl << "Final solution is " << std::endl;
367+
finalGroup->printSolution(finalSolution, finalGroup->getParam("lambda"));//does not work
368+
}
378369

379-
LOCA::destroyGlobalData(globalData);
370+
LOCA::destroyGlobalData(globalData);
380371

381-
for (int i=0;i<_dimSys;i++){
382-
if (_useDomainScaling){
383-
_y[i]=finalSolution(i)/_yScale[i];
384-
}else{
385-
_y[i]=finalSolution(i);
386-
}
387-
}
372+
for (int i=0;i<_dimSys;i++){
373+
if (_useDomainScaling){
374+
_y[i]=finalSolution(i)/_yScale[i];
375+
}else{
376+
_y[i]=finalSolution(i);
377+
}
378+
}
388379

389-
_algLoop->setReal(_y);
390-
try{
391-
_algLoop->evaluate();
392-
}catch(const std::exception &ex)
393-
{
394-
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;
395-
}
380+
_algLoop->setReal(_y);
381+
try{
382+
_algLoop->evaluate();
383+
}catch(const std::exception &ex)
384+
{
385+
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;
386+
}
396387

397-
if (_generateoutput) {
398-
double * rhs = new double[_dimSys];
399-
double sum;
400-
sum=0.0;
401-
_algLoop->getRHS(rhs);
402-
for (int i=0;i<_dimSys;i++) sum+=rhs[i]*rhs[i];
403-
404-
std::cout << "solutionvector=(";
405-
for (int i=0;i<_dimSys;i++) std::cout << std::setprecision (std::numeric_limits<double>::digits10 + 8) << _y[i] << " ";
406-
std::cout << ")" << std::setprecision (6) << std::endl;
407-
std::cout << "rhs =(";
408-
for (int i=0;i<_dimSys;i++) std::cout << rhs[i] << " ";
409-
std::cout << ")" << std::endl;
410-
std::cout << "squared norm of f = " << sum << std::endl;
411-
412-
std::cout << "simtime=" << std::setprecision (std::numeric_limits<double>::digits10 + 1) << _algLoop->getSimTime() << std::endl;
413-
if (rhs) delete rhs;
414-
415-
std::cout << "Solverparamters and StatusTest at simtime " << _algLoop->getSimTime() << std::endl;
416-
paramList->print();
417-
_statusTestsCombo->print(std::cout);
418-
std::cout << "ending solve" << std::endl;
419-
}
388+
if (_generateoutput) {
389+
double * rhs = new double[_dimSys];
390+
double sum;
391+
sum=0.0;
392+
_algLoop->getRHS(rhs);
393+
for (int i=0;i<_dimSys;i++) sum+=rhs[i]*rhs[i];
394+
395+
std::cout << "solutionvector=(";
396+
for (int i=0;i<_dimSys;i++) std::cout << std::setprecision (std::numeric_limits<double>::digits10 + 8) << _y[i] << " ";
397+
std::cout << ")" << std::setprecision (6) << std::endl;
398+
std::cout << "rhs =(";
399+
for (int i=0;i<_dimSys;i++) std::cout << rhs[i] << " ";
400+
std::cout << ")" << std::endl;
401+
std::cout << "squared norm of f = " << sum << std::endl;
402+
403+
std::cout << "simtime=" << std::setprecision (std::numeric_limits<double>::digits10 + 1) << _algLoop->getSimTime() << std::endl;
404+
if (rhs) delete rhs;
405+
406+
std::cout << "ending solve" << std::endl;
407+
}
420408

421409
}
422410
}
@@ -529,9 +517,7 @@ void Nox::modifySolverParameters(const Teuchos::RCP<Teuchos::ParameterList> solv
529517
}
530518
}
531519

532-
533-
534-
bool Nox::checkwhethersolutionisnearby(double const * const y){
520+
bool Nox::CheckWhetherSolutionIsNearby(double const * const y){
535521
double* rhs=new double [_dimSys];
536522
double* rhs2=new double [_dimSys];
537523
double* ypluseps=new double [_dimSys];
@@ -540,10 +526,11 @@ bool Nox::checkwhethersolutionisnearby(double const * const y){
540526
_algLoop->evaluate();
541527
_algLoop->getRHS(rhs);
542528

529+
memcpy(ypluseps,y,_dimSys*sizeof(double));
530+
543531
for (int i=0;i<_dimSys;i++){
544532
// evaluate algLoop at ypluseps=y+eps*e_i and save in rhs2
545-
memcpy(ypluseps,y,_dimSys*sizeof(double));
546-
//ypluseps[i]*=(1.0+std::numeric_limits<double>::epsilon()*2.0);//der Faktor 2.0 ist nur zur Sicherheit und kann spaeter geloescht werden. Alternativ kann std::nextafter oder std::nexttoward verwendet werden
533+
ypluseps[i]=y[i];
547534
ypluseps[i]=std::nextafter(std::nextafter(ypluseps[i],std::numeric_limits<double>::max()),std::numeric_limits<double>::max());
548535
_algLoop->setReal(ypluseps);
549536
_algLoop->evaluate();
@@ -556,8 +543,7 @@ bool Nox::checkwhethersolutionisnearby(double const * const y){
556543
}
557544

558545
// do the same for y-eps*e_i
559-
memcpy(ypluseps,y,_dimSys*sizeof(double));
560-
// ypluseps[i]*=(1.0-std::numeric_limits<double>::epsilon()*2.0);//der Faktor 2.0 ist nur zur Sicherheit.
546+
ypluseps[i]=y[i];
561547
ypluseps[i]=std::nextafter(std::nextafter(ypluseps[i],-std::numeric_limits<double>::max()),-std::numeric_limits<double>::max());
562548
_algLoop->setReal(ypluseps);
563549
_algLoop->evaluate();
@@ -568,22 +554,10 @@ bool Nox::checkwhethersolutionisnearby(double const * const y){
568554
rhssignchange[j]= true;
569555
}
570556
}
571-
572-
//The following does not work since rhssignchange would be rewritten every time i changes.
573-
//It is also troublesome, since rhs and rhs2 are pointer to arrays whereas rhssignchange is a std::vector.
574-
//std::transform(rhs,rhs+_dimSys,rhs2,rhssignchange,[](double a, double b){return (a*b<=0.0);});
575557
}
576-
577558
delete [] ypluseps;
578559
delete [] rhs2;
579560
delete [] rhs;
580-
// just testing
581-
// bool truesol = std::all_of(rhssignchange.begin(),rhssignchange.end(),[](bool a){return a;});
582-
// bool testbool=true;
583-
// for (int i=0;i<_dimSys;i++){
584-
// if (!rhssignchange[i]) testbool=false;
585-
// }
586-
// if (testbool!=truesol) std::cout << "somehow your lambda function does not work" << std::endl;
587561
return std::all_of(rhssignchange.begin(),rhssignchange.end(),[](bool a){return a;});
588562
}
589563

@@ -696,19 +670,19 @@ void Nox::printLogger(){
696670
LOGGER_WRITE_BEGIN("NOX: ",LC_NLS,LL_DEBUG);
697671
LOGGER_WRITE((dynamic_cast<const std::stringstream &>(*_output)).str(),LC_NLS,LL_DEBUG);
698672
LOGGER_WRITE_END(LC_NLS,LL_DEBUG);
699-
//(dynamic_cast<const std::stringstream &>(*_output)).str().clear();//this does nothing. Also, using std::cout somehow miraculously disables the logging.
673+
(dynamic_cast<std::stringstream &>(*_output)).str("");
700674
}
701675
}
702676

703677

704678
//sets printing list in solverParametersPtr
705679
void Nox::addPrintingList(const Teuchos::RCP<Teuchos::ParameterList> solverParametersPtr){
706-
solverParametersPtr->sublist("Printing").set("Output Precision", 15);
707-
solverParametersPtr->sublist("Printing").set("Output Information", NOX::Utils::Error + NOX::Utils::Warning + NOX::Utils::OuterIteration + NOX::Utils::InnerIteration);
708-
if(!_generateoutput){
709-
solverParametersPtr->sublist("Printing").set("Output Stream", _output);
710-
solverParametersPtr->sublist("Printing").set("Error Stream", _output);
711-
}
680+
solverParametersPtr->sublist("Printing").set("Output Precision", 15);
681+
solverParametersPtr->sublist("Printing").set("Output Information", NOX::Utils::Error + NOX::Utils::Warning + NOX::Utils::OuterIteration + NOX::Utils::InnerIteration);
682+
if(!_generateoutput){
683+
solverParametersPtr->sublist("Printing").set("Output Stream", _output);
684+
solverParametersPtr->sublist("Printing").set("Error Stream", _output);
685+
}
712686
}
713687

714688
void Nox::copySolution(const Teuchos::RCP<const NOX::Solver::Generic> solver,double* const algLoopSolution){

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

Lines changed: 0 additions & 13 deletions
Original file line numberDiff line numberDiff line change
@@ -58,8 +58,6 @@ NoxLapackInterface::~NoxLapackInterface()
5858
//we could replace this by the x0 in Nox.cpp
5959
const NOX::LAPACK::Vector& NoxLapackInterface::getInitialGuess()
6060
{
61-
checkdimensionof(*_initialGuess);
62-
6361
if (!_computedinitialguess){
6462
double* x = new double[_dimSys];
6563

@@ -93,9 +91,6 @@ const NOX::LAPACK::Vector& NoxLapackInterface::getInitialGuess()
9391
};
9492

9593
bool NoxLapackInterface::computeActualF(NOX::LAPACK::Vector& f, const NOX::LAPACK::Vector &x){
96-
// if (_algLoop->getSimTime()>0.0) _generateoutput=true;
97-
checkdimensionof(x);
98-
9994
for (int i=0;i<_dimSys;i++){
10095
if (_useDomainScaling){
10196
_xtemp[i]=x(i)/_yScale[i];
@@ -152,9 +147,6 @@ bool NoxLapackInterface::computeActualF(NOX::LAPACK::Vector& f, const NOX::LAPAC
152147
}
153148

154149
bool NoxLapackInterface::computeJacobian(NOX::LAPACK::Matrix<double>& J, const NOX::LAPACK::Vector & x){
155-
156-
checkdimensionof(x);
157-
158150
//setting the forward difference parameters. We divide by the denominator alpha*|x_i|+beta in the computation of the difference quotient. It is similar to the Finite Difference implementation by Nox, which can be found under https://trilinos.org/docs/dev/packages/nox/doc/html/classNOX_1_1Epetra_1_1FiniteDifference.html
159151
double alpha=1.0e-11;
160152
double beta=1.0e-9;
@@ -230,7 +222,6 @@ void NoxLapackInterface::printSolution(const NOX::LAPACK::Vector &x, const doubl
230222
}
231223
}
232224

233-
//replace this function once it is implemented in Trilinos
234225
NOX::LAPACK::Vector NoxLapackInterface::applyMatrixtoVector(const NOX::LAPACK::Matrix<double> &A, const NOX::LAPACK::Vector &x){
235226
NOX::LAPACK::Vector result(A.numRows());
236227
for(int i=0;i<A.numRows();i++){
@@ -314,7 +305,3 @@ bool NoxLapackInterface::computeF(NOX::LAPACK::Vector& f, const NOX::LAPACK::Vec
314305
f.update(_lambda, g, 1.0-_lambda, h);
315306
return true;
316307
}
317-
318-
void NoxLapackInterface::checkdimensionof(const NOX::LAPACK::Vector &x){
319-
if (_dimSys!=x.length()) throw ModelicaSimulationError(ALGLOOP_SOLVER,"Dimension of solution vector is wrong in method of NoxLapackInterface!");
320-
}

0 commit comments

Comments
 (0)