Changeset e8fbdf6 in OpenModelica


Ignore:
Timestamp:
2016-04-06T13:19:25+02:00 (8 years ago)
Author:
hudson <openmodelica@…>
Branches:
Added-citation-metadata, maintenance/v1.14, maintenance/v1.15, maintenance/v1.16, maintenance/v1.17, maintenance/v1.18, maintenance/v1.19, maintenance/v1.20, maintenance/v1.21, maintenance/v1.22, master, omlib-staging
Children:
df3f1c94
Parents:
41d4cef
git-author:
Rüdiger Franke <rdgfranke@…> (04/06/16 13:19:25)
git-committer:
hudson <openmodelica@…> (04/06/16 13:19:25)
Message:

Improve stop and error tests of Newton solver

  • relate rtol to value instead of residual
  • throw nonlinear error if dgesv fails
File:
1 edited

Legend:

Unmodified
Added
Removed
  • SimulationRuntime/cpp/Solver/Newton/Newton.cpp

    r325d842 re8fbdf6  
    9898    dimRHS   = 1,        // Dimension of right hand side of linear system (=b)
    9999    info     = 0;        // Retrun-flag of Fortran code
    100 
    101100  int
    102     totStps  = 0;        // Total number of steps
     101    totSteps = 0;        // Total number of steps taken
     102  double
     103    atol = _newtonSettings->getAtol(),
     104    rtol = _newtonSettings->getRtol();
    103105
    104106  // If initialize() was not called yet
     
    106108    initialize();
    107109
    108   // Get initializeial values from system
     110  // Get current values from system
    109111  _algLoop->getReal(_y);
    110   //_algLoop->evaluate(command);
    111   _algLoop->getRHS(_f);
    112112
    113113  // Reset status flag
    114114  _iterationStatus = CONTINUE;
    115115
    116   while(_iterationStatus == CONTINUE) {
    117     _iterationStatus = DONE;
     116  while (_iterationStatus == CONTINUE) {
    118117    // Check stopping criterion
    119118    calcFunction(_y,_f);
    120     if (totStps) {
     119    if (totSteps) {
     120      _iterationStatus = DONE;
    121121      for (int i=0; i<_dimSys; ++i) {
    122         if (fabs(_f[i]) > _newtonSettings->getAtol() +_newtonSettings->getRtol() * ( fabs(_f[i]))) {
     122        if (fabs(_f[i]) > atol + fabs(_y[i]) * rtol) {
    123123          _iterationStatus = CONTINUE;
    124124          break;
     
    126126      }
    127127    }
    128     else
    129       _iterationStatus = CONTINUE;
    130 
    131     // New right hand side
    132     //calcFunction(_y,_f);
    133128
    134129    if (_iterationStatus == CONTINUE) {
    135       if (totStps < _newtonSettings->getNewtMax()) {
     130      if (totSteps < _newtonSettings->getNewtMax()) {
    136131        // Determination of Jacobian (Fortran-format)
    137         if (_algLoop->isLinear()&&!_algLoop->isLinearTearing()) {
    138           //calcFunction(_yHelp,_fHelp);
     132        if (_algLoop->isLinear() && !_algLoop->isLinearTearing()) {
    139133          const matrix_t& A = _algLoop->getSystemMatrix();
    140134          const double* jac = A.data().begin();
     
    144138          _algLoop->setReal(_y);
    145139          if (info != 0)
    146             throw ModelicaSimulationError(ALGLOOP_SOLVER,"error solving linear tearing system");
     140            throw ModelicaSimulationError(ALGLOOP_SOLVER,
     141              "error solving linear system (dgesv info: " + to_string(info) + ")");
    147142          else
    148143            _iterationStatus = DONE;
    149144        }
    150145        else if (_algLoop->isLinearTearing()) {
    151           long int dimRHS  = 1;          // Dimension of right hand side of linear system (=b)
    152 
    153           long int info  = 0;          // Retrun-flag of Fortran code
     146          long int dimRHS = 1;  // Dimension of right hand side (=b)
     147
     148          long int info = 0;    // Return-flag of Fortran code
    154149
    155150          _algLoop->setReal(_zeroVec);
     
    169164          _algLoop->evaluate();
    170165          if (info != 0)
    171             throw ModelicaSimulationError(ALGLOOP_SOLVER, "error solving linear tearing system");
     166            throw ModelicaSimulationError(ALGLOOP_SOLVER,
     167              "error solving linear tearing system (dgesv info: " + to_string(info) + ")");
    172168          else
    173169            _iterationStatus = DONE;
     
    179175          dgesv_(&_dimSys, &dimRHS, _jac, &_dimSys, _iHelp, _f, &_dimSys, &info);
    180176
    181           if (info != 0) {
    182             // TODO: Throw an error message here.
    183             _iterationStatus = SOLVERERROR;
    184             break;
    185           }
     177          if (info != 0)
     178            throw ModelicaSimulationError(ALGLOOP_SOLVER,
     179              "error solving nonlinear system (iteration: " + to_string(totSteps)
     180              + ", dgesv info: " + to_string(info) + ")");
    186181
    187182          // Increase counter
    188           ++ totStps;
     183          ++ totSteps;
    189184
    190185          // New solution
     
    193188        }
    194189      }
    195       else {
    196         _iterationStatus = SOLVERERROR;
    197         throw ModelicaSimulationError(ALGLOOP_SOLVER,"error solving non linear system");
    198 
    199       }
     190      else
     191        throw ModelicaSimulationError(ALGLOOP_SOLVER,
     192          "error solving nonlinear system (iteration limit: " + to_string(totSteps) + ")");
    200193    }
    201194  }
Note: See TracChangeset for help on using the changeset viewer.