---

SynchronizedRungeKuttaMover.cpp

Go to the documentation of this file.

// SynchronizedRungeKuttaMover.cpp: implementation of the CSynchronizedRungeKuttaMover class.
//

#include "SynchronizedRungeKuttaMover.h"

// Construction/Destruction

CSynchronizedRungeKuttaMover::CSynchronizedRungeKuttaMover(double frequency, double stepSize, double eps)
:CQualityControlRungeKuttaMover(frequency,stepSize,eps)
{
        m_bSyncFlag = false;
}

CSynchronizedRungeKuttaMover::~CSynchronizedRungeKuttaMover()
{

}

void CSynchronizedRungeKuttaMover::Move(double xInitial, double xFinal, ReactionNetwork *pReactionNetwork)
{
        int i;
// 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 min step size to be equal to xInterval at start and max to zero
        this->SetMinStepSize(xInterval);
        this->SetMaxStepSize(0.0);
        m_iTotalStepsTaken = 0;

        // set the member pointer equal to the pointer passed in
        m_pReactionNetwork = pReactionNetwork;
        m_dTime = xInitial;
        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];

        // fill chem[] with starting chemical concentrations
        for(i = 0; i < nChem; i++)
        {
                chem[i] = m_pReactionNetwork->GetChemical(i)->GetAmount();
        }

        // the first call of Move() sets the series of stepsizes, and all subsequent calls
        // use those stepsizes
        int stepCounter = 0;

        while(m_dTime < xFinal)
        {
                // run until xFinal
                while(m_dTime < nextTime)
                {
                        if((m_bSyncFlag == false) || (m_vdSteps.size() == 0))
                        {
                                // get initial derivatives
                                ComputeDerivatives(chem,dChemdt);
                                // compute scaling vector
                                // __max is to ensure chemicals that start at ~ 0 don't get
                                // huge scaling factors, which screws up the stepper
                                for(i = 0; i < nChem; i++)
                                {
                                        m_pdScale[i] = __max(1.0,fabs(chem[i]) + fabs(m_dStepSize*dChemdt[i]) + m_scdTINY);
                                }
                        
                                // don't run over
                                if( ((m_dTime + m_dNextStep) - nextTime) > 0.0 )        
                                {
                                        m_dNextStep = nextTime-m_dTime;
                                }
                        
                                // call the stepper
                                Stepper(nChem,chem,dChemdt);
                                // put size of successful step onto list
                                m_vdSteps.push_back(m_dLastStep);

                                // 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();
                        }
                        else
                        {
                                ComputeDerivatives(chem,dChemdt);
                                m_dStepSize = m_vdSteps[stepCounter];
                                stepCounter++;
                                FourthOrderStep(nChem,chem,dChemdt);
                                m_dTime += m_dStepSize;
                        }
                }
                
                // ensure synchrony with chem and network
                for(int i = 0; i < nChem; i++)
                {
                        m_pReactionNetwork->GetChemical(i)->SetAmount(chem[i]);
                }
                
                m_iCount++;
                Notify();
                nextTime += 1.0;

        }

        delete [] chem;
        delete [] dChemdt;
        delete [] m_pdScale;
        
        // returning after first integration sets this flag to true; 
        // subsequent integrations don't affect this
        m_bSyncFlag = true;
}

---