---

QualityControlRK2TMover.cpp

Go to the documentation of this file.

// QualityControlRK2TMover.cpp: implementation of the QualityControlRK2TMover class.
//

#include "QualityControlRK2TMover.h"

// Variable stepsize RK routine - this one is for the second order 
// trapezoidal step, but changing a few names and deriving the class
// from a different base mover makes it a different order routine. 
// A factory (design pattern) is the more elegant way to have the 
// same hierarchy with less code bloat, but I'm not going to bother
// with that now. Error is computed by taking the step all at once or 
// in two half steps and computing the difference in the answer.
//                                                                                              K.S.B., 3/20/01

// NOTE - ERRCON is essentially the error tolerance; this routine
// is not designed to be terribly accurate (it's only 2nd order 
// anyway) 
const double CQualityControlRK2TMover::m_scdERRCON = 1.0e-3;

const double CQualityControlRK2TMover::m_scdFCOR = 0.06666;

const double CQualityControlRK2TMover::m_scdPGROW = -0.2;

const double CQualityControlRK2TMover::m_scdPSHRINK = -0.25;

const double CQualityControlRK2TMover::m_scdSAFETY = 0.9;

const double CQualityControlRK2TMover::m_scdTINY = 1.0e-30;


// Construction/Destruction

CQualityControlRK2TMover::CQualityControlRK2TMover(double frequency, double stepSize, double eps) 
:CRK2TMover(frequency, stepSize)
{
        this->m_dEps = eps;
        m_dLastStep = stepSize;
        m_dNextStep = stepSize;
        m_dInitialStepSize = stepSize;  // Each "move" starts with initialStepSize
}

CQualityControlRK2TMover::~CQualityControlRK2TMover()
{

}

void CQualityControlRK2TMover::Move(double xInitial, double xFinal, ReactionNetwork *pReactionNetwork)
{
        // do a couple of checks on the integration interval
        double xInterval = xFinal-xInitial;
        if(m_dInitialStepSize > xInterval) {m_dInitialStepSize = xInterval;}
        if(this->MoveTimeIsZero(xInterval)) {return;}
        this->SetMinStepSize(xInterval);
        this->SetMaxStepSize(0.0);
        m_iTotalStepsTaken = 0;

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

        for (int cNum = 0; cNum < nChem; cNum++)
        {
                chem[cNum] = m_pReactionNetwork->GetChemical(cNum)->GetAmount();
        }       

        while(m_dTime < xFinal)
        {
                while(m_dTime < nextTime)
                {
                        ComputeDerivatives(chem,dChemdt);                               
                
                        for (int i = 0; i < nChem; i++)
                        {
                                double scal1 = fabs(chem[i]) + fabs(dChemdt[i]*m_dNextStep) + m_scdTINY;
                                double scal2 = 1.0;
                                if(scal1 > scal2)
                                {
                                        m_pdScale[i] = scal1;
                                }
                                else
                                {
                                        m_pdScale[i] = scal2;
                                }
                        }
                        if( ((m_dTime + m_dNextStep) - nextTime) > 0.0 )        
                        {
                                m_dNextStep = nextTime-m_dTime;
                        }

                        // take a step
                        try
                        {
                                Stepper(nChem,chem,dChemdt);    
                        }
                        catch(std::runtime_error error)
                        {
                                cout << error.what() << endl;
                        }
                        
                        // record largest/smallest/numsteps info
                        if(m_dLastStep < this->GetMinStepSize()) {this->SetMinStepSize(m_dLastStep);}
                        if(m_dLastStep > this->GetMaxStepSize()) {this->SetMaxStepSize(m_dLastStep);}
                        this->IncrementTotalSteps();

                }
                // 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_pdScale;
}

// Takes one quality-controlled Runge-Kutta step.  
// The output (new concentrations after the move) are stored in 
// y[] on exit. 

void CQualityControlRK2TMover::Stepper(int nRHS, double *y, double *dydt)
{
        int i;
        double maxError;
        
        double *dydtTemp1 = new double[nRHS];
        double *dydtTemp2 = new double[nRHS];
        double *yTemp1 = new double[nRHS];
        double *yTemp2 = new double[nRHS];

        m_dStepSize = m_dNextStep;

        while(1)
        {
                for(i = 0; i < nRHS; i++)
                {
                        yTemp1[i] = y[i];
                        yTemp2[i] = y[i];
                        dydtTemp1[i] = dydt[i];
                        dydtTemp2[i] = dydt[i];
                }
                // cut stepSize in half and take two mini-steps
                m_dStepSize = 0.5*m_dStepSize;
                RungeKuttaStep(nRHS,yTemp1,dydtTemp1);
                RungeKuttaStep(nRHS,yTemp1,dydtTemp1);
                
                // undo halving of stepSize
                m_dStepSize = 2.0*m_dStepSize;

                if ((m_dTime + m_dStepSize) == m_dTime) 
                {
                        // need to throw exception
                        cout << "Step size too small in routine RKQC" << endl;
                        exit(1);
                }
                
                RungeKuttaStep(nRHS,yTemp2,dydtTemp2);
                
                maxError = 0.0;
                for (i = 0; i < nRHS; i++)
                {       
                        yTemp2[i] = yTemp1[i] - yTemp2[i];
                        double temp = fabs(yTemp2[i]/m_pdScale[i]);
                        if(maxError < temp) 
                        {
                                maxError = temp;
                        }
                }
                maxError = maxError/m_dEps;
                if(maxError <= 1.0)
                {
                        m_dLastStep = m_dStepSize;
                        if(maxError > m_scdERRCON)
                        {
                                m_dNextStep = m_scdSAFETY*m_dStepSize*pow(maxError,m_scdPGROW);
                        }
                        else
                        {
                                m_dNextStep = 4.0*m_dStepSize;
                        }
                        break;
                }

                m_dStepSize=m_scdSAFETY*m_dStepSize*pow(maxError,m_scdPSHRINK);
        }

        // break out of while loop; step accepted
        
        m_dTime += m_dLastStep;

        for(int i = 0; i < nRHS; i++)
        {
                y[i] = yTemp1[i] + yTemp2[i]*m_scdFCOR;
        }
        
        // free memory

        delete [] yTemp1;
        delete [] yTemp2;
        delete [] dydtTemp1;
        delete [] dydtTemp2;
}

---