---

RK2ExplicitEulerHybridMover.cpp

Go to the documentation of this file.

// RK2ExplicitEulerHybridMover.cpp: implementation of the CRK2ExplicitEulerHybridMover class.
//

#include "RK2ExplicitEulerHybridMover.h"

// The details of this algorithm are described in Ashour and Hanna, 
// Computers chem. Engng., 14 (3), 267 (1990).  It is good for 
// moderately stiff problems at "sloppy tolerance."  It avoids the 
// memory and heavy operation count of most implicit solvers (it does
// not even use the Jacobian).  As coded, the method is second order.

// Static member Alpha giving combination of explicit Euler step and
// RK2 step - don't fool with this unless you know what you're doing!

const double CRK2ExplicitEulerHybridMover::m_scdAlpha=0.74;

// Construction/Destruction

CRK2ExplicitEulerHybridMover::CRK2ExplicitEulerHybridMover(double tolerance, double stepSize)
:CDifferentialEquationMover(1.0, stepSize)
{
        m_dEps = tolerance;
}

CRK2ExplicitEulerHybridMover::~CRK2ExplicitEulerHybridMover()
{

}

void CRK2ExplicitEulerHybridMover::Move(double xInitial, double xFinal, ReactionNetwork *pReactionNetwork)
{
        // do a couple of checks on the integration interval
        double xInterval = xFinal-xInitial;
        if(this->MoveTimeIsZero(xInterval)) {return;}
        if(m_dStepSize > xInterval) {m_dStepSize = xInterval;}

        // set the member pointer to the network pointer passed in
        this->m_pReactionNetwork = pReactionNetwork;
        m_dTime = xInitial;
        // may be unsafe if xInitial has a fractional part!
        m_iCount = (int) xInitial;
        double nextTime = m_dTime;
        
        int nChem = m_pReactionNetwork->GetNumberOfChemicals();
        double *chem = new double[nChem];
        double *dChemdt = new double[nChem];
        m_pdLocalError = new double[nChem];

        // set up local chemical array
        for(int cNum = 0; cNum < nChem; cNum++)
        {
                chem[cNum] = m_pReactionNetwork->GetChemical(cNum)->GetAmount();
        }

        while(m_dTime < xFinal)
        {
                while(m_dTime < nextTime)
                {
                        // make sure the stepsize isn't too small
                        if((m_dTime + m_dStepSize) == m_dTime)
                        {
                                cout << "ERROR: Stepsize too small in RK2/Euler mover." << endl;
                                exit(1);
                        }
                        // don't run over time
                        if( ((m_dTime + m_dStepSize) - nextTime) > 0.0)
                        {
                                m_dStepSize = nextTime - m_dTime;
                        }
                
                        // take a step
                        HybridStep(nChem,chem,dChemdt);
                        m_dTime += m_dStepSize;
                        // adjust the stepsize
                        m_dStepSize = m_dStepSize*SetStepScale(nChem);
                }
                // ensure synchrony with chem and network
                for (int cNum = 0; cNum < nChem; cNum++)
                {
                        m_pReactionNetwork->GetChemical(cNum)->SetAmount(chem[cNum]);
                }

                m_iCount++;
                Notify();
                nextTime += 1.0;
        }
        
        delete [] chem;
        delete [] dChemdt;
        delete [] m_pdLocalError;
}

void CRK2ExplicitEulerHybridMover::HybridStep(int nRHS, double *y, double *dydt)
{
        int i;
        double *yEuler = new double[nRHS];
        double *yRK2 = new double[nRHS];
        double *dydtx = new double[nRHS];
        double *dydtxph = new double[nRHS];

        double halfStep = 0.5*m_dStepSize;

        // take a 1st order Euler step
        ComputeDerivatives(y,dydtx);
        for(i = 0; i < nRHS; i++)
        {
                yEuler[i] = y[i] + m_dStepSize*dydtx[i];
        }
        
        // now take a 2nd order RK step using yEuler
        ComputeDerivatives(yEuler,dydtxph);
        for(i = 0; i < nRHS; i++)
        {
                yRK2[i] = y[i] + halfStep*(dydtx[i] + dydtxph[i]);
        }
        
        // weighted average to get the answer, and simulatneously compute the local
        // error
        for(i = 0; i < nRHS; i++)
        {
                y[i] = m_scdAlpha*yEuler[i] + (1 - m_scdAlpha)*yRK2[i];
                m_pdLocalError[i] = (y[i] - yRK2[i])/y[i];
        }

        // free memory
        delete [] yEuler;
        delete [] yRK2;
        delete [] dydtx;
        delete [] dydtxph;

}

//      The multiplicative factor for the next stepsize is
//                      beta = sqrt(tolerance/max_i(m_pdLocalError));


double CRK2ExplicitEulerHybridMover::SetStepScale(int nRHS)
{
        double maxErr = 0.0;
        for(int i = 0; i < nRHS; i++)
        {
                double absLocalError = fabs(m_pdLocalError[i]);
                if(absLocalError > maxErr) { maxErr = absLocalError;}
        }
        
        return sqrt(m_dEps/maxErr);
}

---