Changeset b1a1225 in OpenModelica


Ignore:
Timestamp:
2016-04-07T20:24:31+02:00 (8 years ago)
Author:
Rüdiger Franke <rdgfranke@…>
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:
98c1ed01
Parents:
f90b512
Message:

Improve stop test and logging of Newton solver

File:
1 edited

Legend:

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

    r058fd5a rb1a1225  
    143143    initialize();
    144144
    145   // Get current values from system
     145  // Get current values and residuals from system
    146146  _algLoop->getReal(_y);
     147  _algLoop->evaluate();
     148  _algLoop->getRHS(_f);
    147149
    148150  // Reset status flag
     
    151153  while (_iterationStatus == CONTINUE) {
    152154    // Check stopping criterion
    153     calcFunction(_y,_f);
    154     if (totSteps) {
     155    if (!_algLoop->isLinear()) {
    155156      _iterationStatus = DONE;
    156157      for (int i=0; i<_dimSys; ++i) {
    157         if (fabs(_f[i]) > atol + fabs(_y[i]) * rtol) {
     158        if (fabs(_f[i]) > atol + rtol * fabs(_y[i])) {
    158159          _iterationStatus = CONTINUE;
    159160          break;
     
    162163    }
    163164    if (_iterationStatus == CONTINUE) {
    164       LogSysVec(_algLoop, "y" + to_string(totSteps), _y);
    165165      if (totSteps < _newtonSettings->getNewtMax()) {
    166166        // Determination of Jacobian (Fortran-format)
     
    170170          memcpy(_jac, jac, _dimSys*_dimSys*sizeof(double));
    171171          dgesv_(&_dimSys, &dimRHS, _jac, &_dimSys, _iHelp, _f, &_dimSys, &info);
    172           memcpy(_y,_f,_dimSys*sizeof(double));
     172          memcpy(_y, _f, _dimSys*sizeof(double));
    173173          _algLoop->setReal(_y);
    174174          if (info != 0)
     
    185185          const matrix_t& A_sparse = _algLoop->getSystemMatrix();
    186186          //m_t A_dense(A_sparse);
    187 
    188187          const double* jac = A_sparse.data().begin();
    189 
    190188          memcpy(_jac, jac, _dimSys*_dimSys*sizeof(double));
    191189          dgesv_(&_dimSys, &dimRHS, _jac, &_dimSys, _iHelp, _f, &_dimSys, &info);
    192           for (int i=0; i<_dimSys; i++)
    193             _y[i]=-_f[i];
     190          for (int i = 0; i < _dimSys; i++)
     191            _y[i] = -_f[i];
    194192          _algLoop->setReal(_y);
    195193          _algLoop->evaluate();
     
    201199        }
    202200        else {
     201          LogSysVec(_algLoop, "y" + to_string(totSteps), _y);
    203202          calcJacobian();
    204203
     
    214213          ++ totSteps;
    215214
    216           // New solution
    217           for (int i=0; i<_dimSys; ++i)
     215          // New iterate
     216          for (int i = 0; i < _dimSys; ++i)
    218217            _y[i] -= _newtonSettings->getDelta() * _f[i];
     218          calcFunction(_y, _f);
    219219        }
    220220      }
Note: See TracChangeset for help on using the changeset viewer.