autostepme.h

/*************************************************************************************
 * MechSys - A C++ library to simulate (Continuum) Mechanical Systems                *
 * Copyright (C) 2005 Dorival de Moraes Pedroso <dorival.pedroso at gmail.com>       *
 * Copyright (C) 2005 Raul Dario Durand Farfan  <raul.durand at gmail.com>           *
 *                                                                                   *
 * This file is part of MechSys.                                                     *
 *                                                                                   *
 * MechSys is free software; you can redistribute it and/or modify it under the      *
 * terms of the GNU General Public License as published by the Free Software         *
 * Foundation; either version 2 of the License, or (at your option) any later        *
 * version.                                                                          *
 *                                                                                   *
 * MechSys is distributed in the hope that it will be useful, but WITHOUT ANY        *
 * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A   *
 * PARTICULAR PURPOSE. See the GNU General Public License for more details.          *
 *                                                                                   *
 * You should have received a copy of the GNU General Public License along with      *
 * MechSys; if not, write to the Free Software Foundation, Inc., 51 Franklin Street, *
 * Fifth Floor, Boston, MA 02110-1301, USA                                           *
 *************************************************************************************/

#ifndef MECHSYS_AUTOSTEPME_H
#define MECHSYS_AUTOSTEPME_H

#ifdef HAVE_CONFIG_H
  #include "config.h"
#else
  #ifndef REAL
    #define REAL double
  #endif
#endif

#include <cmath>
#ifdef DO_DEBUG
    #include <sstream>
    #include <string>
    #include <iomanip>
#endif

// MechSys
#include "numerical/integschemesctes.h"

template<typename xType, typename yType, typename zType, typename Instance>
class AutoStepME
{
public:

    typedef int (Instance::*pTgDyDz) (xType const &  x, yType const &  y, zType const &  z,
                                      xType const & dx, yType       & dy, zType       & dz) const;

    typedef REAL (Instance::*pCalcErr) (yType const & Ey, yType const & y_high,
                                        zType const & Ez, zType const & z_high) const;

    AutoStepME(IntegSchemesCtes const & ISC)
        : _maxSS (ISC.ME_maxSS()),
          _STOL  (ISC.ME_STOL ()),
          _dTini (ISC.ME_dTini()),
          _mMin  (ISC.ME_mMin ()),
          _mMax  (ISC.ME_mMax ())
    {}


    int Solve(Instance const * p2Inst, pTgDyDz tg_dy_dz, pCalcErr calc_error,
              xType       &  x,  yType &  y,  zType &  z,
              xType const & Dx)
    {
        REAL  T = 0.0;    // Pseudo time
        REAL dT = _dTini; // First increment of pseudo time
        REAL  smult;      // Multiplier for step-size
        REAL  error;      // Error estimative

        xType dx_k;  // "driver" increment
        yType dy_FE; // Forward Euler increment
        zType dz_FE; // Forward Euler increment

        xType x_FE; // FE "state"
        yType y_FE; // FE "state"
        zType z_FE; // FE "state" 

        yType dy_int; // Intermediary increment evaluate with a tangent given for the FE "state"
        zType dz_int; // Intermediary increment evaluate with a tangent given for the FE "state"

        yType y_ME; // Modified Euler "state" (high order)
        zType z_ME; // Modified Euler "state" (high order)

        yType E_y; // Estimated error on y
        zType E_z; // Estimated error on z

#ifdef DO_DEBUG
        _results << std::setw(15) << "T";
        _results << std::setw(15) << "x";
        _results << std::setw(15) << "y";
        _results << std::setw(15) << "z" << std::endl;
        _results << std::setw(15) << T;
        _results << std::setw(15) << x;
        _results << std::setw(15) << y;
        _results << std::setw(15) << z << std::endl;
#endif
        int Flag; // Flag received from tg_dy_dz: -1 => failed
                  //                               0 => continue
                  //                               1 => stop

        for (int k=1; k<=_maxSS; ++k)
        {
            if (T>=1.0) return k;

            dx_k = Dx * dT;

            Flag = (p2Inst->*tg_dy_dz)(x, y, z, dx_k, dy_FE, dz_FE); // Out: dy_FE and dz_FE
            if      (Flag==-1) return -1;
            else if (Flag==1)  return -2;

            x_FE = x + dx_k;
            y_FE = y + dy_FE;
            z_FE = z + dz_FE;

            Flag = (p2Inst->*tg_dy_dz)(x_FE, y_FE, z_FE, dx_k, dy_int, dz_int); // Out: dy_int and dz_int
            if      (Flag==-1) return -1;
            else if (Flag==1)  return -2;

            y_ME = y + 0.5*( dy_FE + dy_int );
            z_ME = z + 0.5*( dz_FE + dz_int );

            E_y = y_ME - y_FE;
            E_z = z_ME - z_FE;

            error = (p2Inst->*calc_error)(E_y, y_ME, E_z, z_ME);
            smult = 0.9*sqrt(_STOL/error);
            
            if (error<=_STOL) // Step accepted
            {
                T = T + dT; // Actualize "pseudo-time"
                x = x_FE; // Next x
                y = y_ME; // Next y
                z = z_ME; // Next z
#ifdef DO_DEBUG
            _results << std::setw(15) << T;
            _results << std::setw(15) << x;
            _results << std::setw(15) << y;
            _results << std::setw(15) << z << std::endl;
#endif
                if (smult>_mMax) smult=_mMax; // Upper bound step size
            }
            else // Step rejected
                if (smult<_mMin) smult=_mMin; // Lower bound step size

            dT = smult * dT;        // Update step size
            if (dT>1.0-T) dT=1.0-T; // Last step
        }
        return -1; // ERROR: ME did not converge!
    }

#ifdef DO_DEBUG
    std::string Results() { return _results.str(); }
#endif

private:
    int  _maxSS; 
    REAL _STOL;  
    REAL _dTini; 
    REAL _mMin;  
    REAL _mMax;  

#ifdef DO_DEBUG
    std::ostringstream _results;
#endif

}; // AutoStepME

// \cond implementation

// \endcond implementation

#endif // #define MECHSYS_AUTOSTEPME_H

---