---

QualityControlRungeKuttaMover.cpp

Go to the documentation of this file.

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

#include "QualityControlRungeKuttaMover.h"


const double CQualityControlRungeKuttaMover::m_scdERRCON = 6.0e-4;

const double CQualityControlRungeKuttaMover::m_scdFCOR = 0.06666;

const double CQualityControlRungeKuttaMover::m_scdPGROW = -0.2;

const double CQualityControlRungeKuttaMover::m_scdPSHRINK = -0.25;

const double CQualityControlRungeKuttaMover::m_scdSAFETY = 0.9;

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


// Construction/Destruction

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

CQualityControlRungeKuttaMover::~CQualityControlRungeKuttaMover()
{

}

void CQualityControlRungeKuttaMover::ResetStepSize()
{
        this->m_dInitialStepSize = m_dDefaultStepSize;
}

double CQualityControlRungeKuttaMover::GetStepSize() const
{
        return m_dInitialStepSize;
}

void CQualityControlRungeKuttaMover::SetStepSize(double newSize)
{
        this->m_dInitialStepSize = newSize;
}

void CQualityControlRungeKuttaMover::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_dInitialStepSize > xInterval) {m_dInitialStepSize = xInterval;}
        // set the increment for nextTime appropriately (all QC integrators need this line)
        double increment;
        if(xInterval < 1.0){increment = xInterval;}
        else{increment = 1.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];
        double minStep = 0.1*m_dStepSize;

        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
                        Stepper(nChem,chem,dChemdt);    

                }
                // ensure synchrony with chem and network
                for (int cNum = 0; cNum < nChem; cNum++)
                {
                        m_pReactionNetwork->GetChemical(cNum)->SetAmount(chem[cNum]);
                }

                m_iCount++;
                Notify();
                nextTime += increment;
        }
        
        delete [] chem;
        delete [] dChemdt;
        delete [] m_pdScale;
}

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

void CQualityControlRungeKuttaMover::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
                m_dStepSize = 0.5*m_dStepSize;
                FourthOrderStep(nRHS,yTemp1,dydtTemp1);
                ComputeDerivatives(yTemp1,dydtTemp1);
                FourthOrderStep(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);
                }
                
                FourthOrderStep(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;
}

---