Skip to content

Commit faa82e6

Browse files
qichenghuaOpenModelica-Hudson
qichenghua
authored and
OpenModelica-Hudson
committed
fixed error introduced in 31424fc19ae5b7e29f5a17e102eb4d6e0a29d311
1 parent 330aa9e commit faa82e6

File tree

2 files changed

+60
-143
lines changed

2 files changed

+60
-143
lines changed

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

Lines changed: 1 addition & 113 deletions
Original file line numberDiff line numberDiff line change
@@ -4,10 +4,10 @@
44
#include <Solver/Nox/FactoryExport.h>
55

66

7+
#include <Core/Utils/extension/logger.hpp>
78
#include <Solver/Nox/NoxLapackInterface.h>
89
#include "Teuchos_StandardCatchMacros.hpp"
910

10-
#include <Core/Utils/extension/logger.hpp>
1111
#include <Solver/Nox/Nox.h>
1212
#include <Solver/Nox/NoxSettings.h>
1313

@@ -730,116 +730,4 @@ NOX::StatusTest::StatusType Nox::BasicNLSsolve(){
730730
return status;
731731
}
732732

733-
NOX::StatusTest::StatusType Nox::secondBasicNLSsolve(){
734-
NOX::StatusTest::StatusType status;
735-
int iter = 0;
736-
if (_generateoutput) std::cout << "starting solving algloop " << _algLoop->getEquationIndex() << std::endl;
737-
738-
if (_firstCall) initialize();
739-
740-
_statusTestNormF = Teuchos::rcp(new NOX::StatusTest::NormF(1.0e-13));
741-
_statusTestMaxIters = Teuchos::rcp(new NOX::StatusTest::MaxIters(100));
742-
_statusTestStagnation = Teuchos::rcp(new NOX::StatusTest::Stagnation(15,0.99));
743-
_statusTestDivergence = Teuchos::rcp(new NOX::StatusTest::Divergence(1.0e13));
744-
745-
_statusTestsCombo = Teuchos::rcp(new NOX::StatusTest::Combo(NOX::StatusTest::Combo::OR, _statusTestNormF, _statusTestMaxIters));
746-
_statusTestsCombo->addStatusTest(_statusTestStagnation);
747-
_statusTestsCombo->addStatusTest(_statusTestDivergence);
748-
749-
_solverParametersPtr = Teuchos::rcp(new Teuchos::ParameterList);
750-
_solverParametersPtr->sublist("Line Search").set("Method", "Backtrack");
751-
752-
if (_generateoutput){
753-
_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.
754-
}else{
755-
_solverParametersPtr->sublist("Printing").set("Output Information", NOX::Utils::Error);
756-
}
757-
758-
_noxLapackInterface = Teuchos::rcp(new NoxLapackInterface (_algLoop,-1));//this also gets the nominal values
759-
760-
_iterationStatus=CONTINUE;
761-
762-
while(_iterationStatus==CONTINUE){
763-
iter++;
764-
765-
// Reset initial guess
766-
memcpy(_y,_y0,_dimSys*sizeof(double));
767-
_algLoop->setReal(_y);
768-
769-
_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
770-
771-
try{
772-
_solver = NOX::Solver::buildSolver(_grp, _statusTestsCombo, _solverParametersPtr);
773-
}
774-
catch(const std::exception &ex)
775-
{
776-
std::cout << std::endl << "sth went wrong during solver building, with error message" << ex.what() << std::endl;
777-
throw ModelicaSimulationError(ALGLOOP_SOLVER,"solver building error");
778-
}
779-
780-
781-
if (_generateoutput) {
782-
double* rhsssss = new double[_dimSys];//stores f(x)
783-
double* xp = new double[_dimSys];//stores x temporarily
784-
_algLoop->getReal(xp);
785-
_algLoop->getRHS(rhsssss);
786-
std::cout << "we are at position x=(";
787-
for (int i=0;i<_dimSys;i++){
788-
std::cout << xp[i] << " ";
789-
}
790-
std::cout << ")" << std::endl;
791-
std::cout << "the right hand side is given by (";
792-
for (int i=0;i<_dimSys;i++){
793-
std::cout << rhsssss[i] << " ";
794-
}
795-
std::cout << ")" << std::endl;
796-
std::cout << std::endl;
797-
delete [] rhsssss;
798-
delete [] xp;
799-
800-
std::cout << "we are solving with the following options:" <<std::endl;
801-
_solverParametersPtr->print();
802-
std::cout << std::endl;
803-
}
804-
805-
try{
806-
status = _solver->solve();
807-
}
808-
catch(const std::exception &ex)
809-
{
810-
std::cout << std::endl << "sth went wrong during solving, with error message" << ex.what() << std::endl;
811-
throw ModelicaSimulationError(ALGLOOP_SOLVER,"solving error");
812-
}
813-
814-
// Get the answer
815-
NOX::LAPACK::Group solnGrp = dynamic_cast<const NOX::LAPACK::Group&>(_solver->getSolutionGroup());
816-
NOX::LAPACK::Vector Lapacksolution=dynamic_cast<const NOX::LAPACK::Vector&>(solnGrp.getX());
817-
818-
for (int i=0;i<_dimSys;i++){
819-
if (_useDomainScaling){
820-
_y[i]=Lapacksolution(i)/_yScale[i];
821-
}else{
822-
_y[i]=Lapacksolution(i);
823-
}
824-
}
825-
826-
_algLoop->setReal(_y);
827-
_algLoop->evaluate();
828-
829-
if (_generateoutput) {
830-
std::cout << "solutionvector=(";
831-
for (int i=0;i<_dimSys;i++) std::cout << _y[i] << " ";
832-
std::cout << ")" << std::endl;
833-
}
834-
835-
if (status==NOX::StatusTest::Converged){
836-
if (_generateoutput) std::cout << "simtime=" << _algLoop->getSimTime() << std::endl;
837-
_iterationStatus=DONE;
838-
}
839-
}
840-
if (_generateoutput) std::cout << "ending solve" << std::endl;
841-
842-
return status;
843-
}
844-
845733
/** @} */ // end of solverNox

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

Lines changed: 59 additions & 30 deletions
Original file line numberDiff line numberDiff line change
@@ -91,8 +91,10 @@ const NOX::LAPACK::Vector& NoxLapackInterface::getInitialGuess()
9191

9292
if (_generateoutput) {
9393
std::cout << "Initial guess is given by " << std::endl;
94-
for(int i=0;i<_dimSys;i++) std::cout << (*_initialGuess)(i) << " ";
95-
std::cout << std::endl;
94+
// for(int i=0;i<_dimSys;i++) std::cout << (*_initialGuess)(i) << " ";
95+
// std::cout << std::endl;
96+
// std::cout << "or" << std::endl;
97+
_initialGuess->print(std::cout);
9698
}
9799

98100
if (_useFunctionValueScaling){
@@ -132,24 +134,40 @@ bool NoxLapackInterface::computeActualF(NOX::LAPACK::Vector& f, const NOX::LAPAC
132134
}
133135
std::cout << ")" << std::endl;
134136
std::cout << "The position seen by NOX is given by x=(";
135-
for (int i=0;i<_dimSys;i++){
136-
std::cout << std::setprecision (std::numeric_limits<double>::digits10 + 8) << x(i) << " ";
137-
}
138-
std::cout << ")" << std::endl;
137+
// for (int i=0;i<_dimSys;i++){
138+
// std::cout << std::setprecision (std::numeric_limits<double>::digits10 + 8) << x(i) << " ";
139+
// }
140+
// std::cout << ")" << std::endl;
141+
// std::cout << "or" << std::endl;
142+
x.print(std::cout);
139143
}
140144

141145
_algLoop->setReal(_xtemp);
142-
//_algLoop->getRHS(_rhs);
146+
_algLoop->getRHS(_rhs);
143147
try{
144148
_algLoop->evaluate();
145149
_algLoop->getRHS(_rhs);
146150
}catch(const std::exception &ex)
147151
{
148152
if (_generateoutput) std::cout << "calculating right hand side failed with error message:" << std::endl << ex.what() << std::endl;
149153
//the following should be done when some to be implemented flag like "continue if function evaluation fails" is activated.
154+
155+
156+
//even newer, experimental version. Delete comments, when this is tested.
157+
// for(int i=0;i<_dimSys;i++){
158+
// _rhs[i]=_fScale[i]*(std::abs(x(i)-(getInitialGuess())(i))+1)*((_rhs[i]<0.0) ? -1.0e6 : 1.0e6);
159+
// this has conical form.
160+
// it is based on the assumption that getInitialGuess is contained in the domain.
161+
// We have two goals: One is to ensure, that the rhs is at least 1.0e6.
162+
// The other is to make sure, that if there is an open ball around the area where we cannot evaluate, that the next newton iterate is the initial guess.
163+
// That is actually not good, the initial guess is where we started...
164+
// _rhs[i]=_fScale[i]*(std::abs(x(i))+1)*((_rhs[i]<0.0) ? -1.0e6 : 1.0e6);
165+
// this has conical form with center 0.
166+
// }
167+
150168
//new, experimental version. Delete comments, when this is tested.
151169
for(int i=0;i<_dimSys;i++){
152-
_rhs[i]= (_rhs[i]<0.0) ? -1.0e6 : 1.0e6;
170+
_rhs[i]= ((_rhs[i]<0.0) ? -1.0e6 : 1.0e6);
153171
}
154172
//Maybe this should have conical form, ie. in case of high values, use 1.0e6*x+1.0e6 instead.
155173

@@ -173,7 +191,7 @@ bool NoxLapackInterface::computeActualF(NOX::LAPACK::Vector& f, const NOX::LAPAC
173191
std::cout << _rhs[i] << " ";
174192
}
175193
std::cout << ")" << std::endl;
176-
std::cout << "the right hand side seen by NOX is given by (";
194+
// std::cout << "the right hand side seen by NOX is given by (";
177195
}
178196
for (int i=0;i<_dimSys;i++){
179197

@@ -183,15 +201,17 @@ bool NoxLapackInterface::computeActualF(NOX::LAPACK::Vector& f, const NOX::LAPAC
183201
f(i)=_rhs[i];
184202
}
185203
//checking for infinity.
186-
if (f(i)>=std::numeric_limits<double>::max()) f(i)=1.0e6;
187-
if (f(i)<=-std::numeric_limits<double>::max()) f(i)=-1.0e6;
204+
if (f(i)>=std::numeric_limits<double>::max()) f(i)=1.0e12;
205+
if (f(i)<=-std::numeric_limits<double>::max()) f(i)=-1.0e12;
188206
//checking for NaN. Do NOT delete the next line, it makes sense.
189-
if (!(f(i)==f(i))) f(i)=1.0e6;
190-
if (_generateoutput) std::cout << f(i) << " ";
207+
if (!(f(i)==f(i))) f(i)=1.0e12;
208+
// if (_generateoutput) std::cout << f(i) << " ";
191209
}
192210
if (_generateoutput){
193-
std::cout << ")" << std::endl;
194-
std::cout << std::endl;
211+
// std::cout << ")" << std::endl;
212+
// std::cout << std::endl;
213+
// std::cout << "or also" << std::endl;
214+
// f.print(std::cout);
195215
}
196216
return true;
197217
}
@@ -211,13 +231,16 @@ bool NoxLapackInterface::computeJacobian(NOX::LAPACK::Matrix<double>& J, const N
211231

212232
if (_generateoutput){
213233
std::cout << "we are at simtime " << _algLoop->getSimTime() << " and at position (seen by NOX) x=(";
214-
for (int i=0;i<_dimSys;i++){
215-
std::cout << x(i) << " ";
216-
}
217-
std::cout << ")" << std::endl;
218-
std::cout << std::endl;
234+
// for (int i=0;i<_dimSys;i++){
235+
// std::cout << x(i) << " ";
236+
// }
237+
// std::cout << ")" << std::endl;
238+
// std::cout << std::endl;
239+
// std::cout << "or" << std::endl;
240+
x.print(std::cout);
219241
}
220242

243+
221244
for (int i=0;i<_dimSys;i++){
222245
//adding the denominator of the difference quotient
223246
xplushei(i)+=alpha*std::abs(xplushei(i))+beta;
@@ -228,15 +251,19 @@ bool NoxLapackInterface::computeJacobian(NOX::LAPACK::Matrix<double>& J, const N
228251
xplushei(i)=x(i);
229252
}
230253

254+
//if (_generateoutput){
255+
// std::cout << "the Jacobian is given by (the transpose of) " << std::endl;
256+
// for (int i=0;i<_dimSys;i++){
257+
// for (int j=0;j<_dimSys;j++){
258+
// std::cout << J(j,i) << " ";
259+
// }
260+
// std::cout << std::endl;
261+
// }
262+
// std::cout << std::endl << "done computing Jacobian" << std::endl;
263+
//}
231264
if (_generateoutput){
232-
std::cout << "the Jacobian is given by (the transpose of) " << std::endl;
233-
for (int i=0;i<_dimSys;i++){
234-
for (int j=0;j<_dimSys;j++){
235-
std::cout << J(j,i) << " ";
236-
}
237-
std::cout << std::endl;
238-
}
239-
std::cout << std::endl << "done computing Jacobian" << std::endl;
265+
std::cout << "the Jacobian is given by" << std::endl;
266+
J.print(std::cout);//is correct, no transposing necessary.
240267
}
241268
return true;
242269
}
@@ -250,8 +277,10 @@ void NoxLapackInterface::printSolution(const NOX::LAPACK::Vector &x, const doubl
250277
{
251278
if(_generateoutput){
252279
std::cout << "At parameter value: " << std::setprecision(8) << conParam << " the solution vector (norm=" << x.norm() << ") is" << std::endl;
253-
for (int i=0; i<_dimSys; i++) std::cout << " " << x(i);
254-
std::cout << ")" << std::endl;
280+
// for (int i=0; i<_dimSys; i++) std::cout << " " << x(i);
281+
// std::cout << ")" << std::endl;
282+
// std::cout << "or" << std::endl;
283+
x.print(std::cout);
255284
std::cout << "Simtime: " << _algLoop->getSimTime() << std::endl;
256285
}
257286
}

0 commit comments

Comments
 (0)