---

DifferentialEquationMover.cpp

Go to the documentation of this file.

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

#include "DifferentialEquationMover.h"

// Construction/Destruction

CDifferentialEquationMover::CDifferentialEquationMover()
{

}

CDifferentialEquationMover::CDifferentialEquationMover(double frequency, double stepSize)
:CReactionMover(0,frequency)
{
        this->m_dStepSize = stepSize;
        this->m_dDefaultStepSize = stepSize;
}

CDifferentialEquationMover::~CDifferentialEquationMover()
{

}

void CDifferentialEquationMover::SetStepSize(double newSize)
{
        this->m_dStepSize = newSize;
}

void CDifferentialEquationMover::ResetStepSize()
{
        this->m_dStepSize = m_dDefaultStepSize;
}

void CDifferentialEquationMover::ComputeDerivatives(double *chem, double *dChemdt)
{
        int cNum;
        // Set dChemdt = 0.0
        // Set chemicals to values in concentration array

        int nChem = m_pReactionNetwork->GetNumberOfChemicals();
        int nReac = m_pReactionNetwork->GetNumberOfReactions();
        

        for (cNum = 0; cNum < nChem; cNum++)
        {
                dChemdt[cNum] = 0.0;
                Chemical* chemical = m_pReactionNetwork->GetChemicals()->at(cNum);
                //Chemical *chemical = m_pReactionNetwork->GetChemical(cNum);
                chemical->SetAmount(chem[cNum]);
        }

        // Give contribution to dChemdt from each reaction
        std::vector<double> *reactionRates = m_pReactionNetwork->GetReactionRates();
        std::vector<Reaction *> *reactions = m_pReactionNetwork->GetReactions();
        
        for (int rNum = 0; rNum < nReac; rNum++)
        {
                std::vector<Chemical *> *localChemicals = reactions->at(rNum)->GetChemicals();
                double reactionRate = reactionRates->at(rNum);
                int nLC = localChemicals->size();
                for (int lcNum = 0; lcNum < nLC; lcNum++)
                {
                        cNum = localChemicals->at(lcNum)->GetChemicalNumber();
                        dChemdt[cNum] += (reactions->at(rNum)->GetNumberChangedByReaction(lcNum))*reactionRate;
                }
        }
        return;
}

void CDifferentialEquationMover::ComputeJacobian(double *chem, double **dRHSdChem)
{
        int i,j;
        int nChem = m_pReactionNetwork->GetNumberOfChemicals();
        int nReac = m_pReactionNetwork->GetNumberOfReactions();

        // zero out the jacobian array
        for(i = 0; i < nChem; i++)
        {
                for(j = 0; j < nChem; j++)
                {
                        dRHSdChem[i][j] = 0.0;
                }
        }
        // make sure chemicals are current
        for (int cNum = 0; cNum < nChem; cNum++)
        {
                Chemical* chemical = m_pReactionNetwork->GetChemicals()->at(cNum);
                chemical->SetAmount(chem[cNum]);
        }

        std::vector<Reaction::JElement *> *jelements;
        std::vector<Reaction *> *reactions = m_pReactionNetwork->GetReactions();
        
        for (int rNum = 0; rNum < nReac; rNum++)
        {
                std::vector<Chemical *> *localChemicals = reactions->at(rNum)->GetChemicals();
                int nLC = localChemicals->size();
                for (int lcNum = 0; lcNum < nLC; lcNum++)
                {
                        int nChg = reactions->at(rNum)->GetNumberChangedByReaction(lcNum);
                        if(nChg != 0)
                        {
                                double sign = fabs((double) nChg)/nChg;
                                int cNum = localChemicals->at(lcNum)->GetChemicalNumber();
                                jelements = reactions->at(rNum)->GetChemicalJacobian();
                                for(int nJac = 0; nJac < jelements->size(); nJac++)
                                {
                                        int wrt = jelements->at(nJac)->GetWithRespectTo();
                                        double jvalue = jelements->at(nJac)->GetJValue();
                                        dRHSdChem[cNum][wrt] += sign*jvalue;
                                }
                        }
                        
                }
        }
}

---