---

StiffBulirschStoerMover.cpp

Go to the documentation of this file.

#include "StiffBulirschStoerMover.h"

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

const int CStiffBulirschStoerMover::m_sciKMAXX = 7;

// this needs to be KMAXX+1
const int CStiffBulirschStoerMover::m_sciIMAXX = 8;

const double CStiffBulirschStoerMover::m_scdSAFE1 = 0.25;

const double CStiffBulirschStoerMover::m_scdSAFE2 = 0.7;

const double CStiffBulirschStoerMover::m_scdREDMAX = 1.0e-05;

const double CStiffBulirschStoerMover::m_scdREDMIN = 0.7;

const double CStiffBulirschStoerMover::m_scdSCALMX = 0.1;

CStiffBulirschStoerMover::CStiffBulirschStoerMover(double frequency, double stepSize, double eps)
:CDifferentialEquationMover(frequency,stepSize)
{
        m_dLastStep = stepSize;
        m_dNextStep = stepSize;
        m_dInitialStepSize = stepSize;
        m_dEps = eps;
}

CStiffBulirschStoerMover::~CStiffBulirschStoerMover(void)
{
}

void CStiffBulirschStoerMover::Move(double xInitial, double xFinal, ReactionNetwork *pReactionNetwork)
{
        int i;
        // do a couple of checks on the integration interval
        double xInterval = xFinal-xInitial;
        // set the member pointer equal to the pointer passed in
        m_pReactionNetwork = pReactionNetwork;
        // if there's no moving to be done and xI = xF = integer, take data there
        if(this->MoveTimeIsZero(xInterval)) 
        {
                if(xInitial > (int) xInitial)
                {
                }
                else
                {
                        m_iCount = (int) xInitial;
                        m_iCount++;
                        Notify();
                }
                // take a data point 
                return;
        }
        if(m_dInitialStepSize > xInterval) {m_dInitialStepSize = xInterval;}

        m_dTime = xInitial;
        double nextTime;
        // set up the count variable
        m_iCount = (int) xInitial;
        m_iCount++;
        nextTime = m_iCount;

        m_dStepSize = m_dInitialStepSize;
        m_dLastStep = m_dInitialStepSize;
        m_dNextStep = m_dInitialStepSize;
        
        int nChem = m_pReactionNetwork->GetNumberOfChemicals();
        double *chem = new double[nChem];
        double *dChemdt = new double[nChem];
        m_pdScale = new double[nChem];
        m_pdError = new double[nChem];

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

        while(m_dTime < xFinal)
        {
                // run until xFinal or nextTime (nextTime might be > xFinal)
                while( m_dTime < __min(xFinal,nextTime) )
                {
                        // 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] = fabs(chem[i]) + fabs(m_dStepSize*dChemdt[i]) + m_scdTINY;
                                m_pdScale[i] = __max(1.0,fabs(chem[i]));
                        }
                        
                        // don't run over - nextTime might be > xFinal
                        if( (((m_dTime + m_dNextStep) - nextTime) > 0.0) || (((m_dTime + m_dNextStep) - xFinal) > 0.0) )        
                        {
                                m_dNextStep = __min(nextTime-m_dTime,xFinal-m_dTime);
                        }

                        // call the stepper
                        StiffBSStep(nChem,chem,dChemdt);

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

        delete [] chem;
        delete [] dChemdt;
        delete [] m_pdScale;
        delete [] m_pdError;
}


void CStiffBulirschStoerMover::StiffBSStep(int nRHS, double *y, double *dydt)
{
        int i,iq,k,kk,km;
        int first = 1;
        int kmax,kopt;
        int nvold = -1;
        double epsold = -1.0;
        double xnew;
        double eps1,errmax,fact,red,scale,work,wrkmin,xest;
        int reduct = 0;
        bool exitflag = false;

        m_dStepSize = m_dNextStep;

        // memory allocation
        double *a = new double[m_sciIMAXX];
        double **alf = new double *[m_sciKMAXX];
        alf[0] = new double[(m_sciKMAXX )*(m_sciKMAXX)];
        for(i = 1; i < m_sciKMAXX; i++)
        {
                alf[i] = &(alf[0][i*(m_sciKMAXX)]);
        }
        int *nseq = new int[m_sciIMAXX];
        
        // step sequence
        nseq[0] = 2;
        nseq[1] = 6;
        nseq[2] = 10;
        nseq[3] = 14;
        nseq[4] = 22;
        nseq[5] = 34;
        nseq[6] = 50;
        nseq[7] = 70;

        double **d = new double *[nRHS];
        d[0] = new double[nRHS*(m_sciKMAXX)];
        for(i = 1; i < nRHS; i++)
        {
                d[i] = &(d[0][i*m_sciKMAXX]);
        }
        double **dfdy = new double *[nRHS];
        dfdy[0] = new double[nRHS*nRHS];
        for(i = 1; i < nRHS; i++)
        {
                dfdy[i] = &(dfdy[0][i*nRHS]);
        }
        double *x = new double[m_sciKMAXX];
        double *yerr = new double[nRHS];
        double *ysav = new double[nRHS];
        double *yseq = new double[nRHS];

        m_dStepSize = m_dNextStep;

        if(m_dEps != epsold || nRHS != nvold) // Reinitialize if either have changed
        {
                m_dNextStep = xnew = -1.0e29;
                eps1 = m_scdSAFE1*m_dEps;
                // compute work coefficients A_k
                a[0] = nseq[0] + 1;
                for(k = 0; k < m_sciKMAXX; k++)
                {
                        a[k+1] = a[k]*nseq[k+1];
                }
                
                // compute alpha(k,q)
                for(iq = 1; iq < m_sciKMAXX; iq++)
                {
                        for(k = 0; k < iq; k++)
                        {
                                double exponent = (a[k+1]-a[iq+1])/((a[iq+1] - a[0] + 1.0)*(2*k+3));
                                alf[k][iq] = pow(eps1,exponent);
                        }
                } // end iq loop
                
                epsold = m_dEps;
                nvold = nRHS;
                
                // add cost of Jacobian calculations to work coefficients
                a[0] += nRHS;
                for(k = 0; k < m_sciKMAXX; k++)
                {
                        a[k+1] = a[k] + nseq[k+1];
                }

                // determine optimal row number for convergence
                for(kopt = 1; kopt < m_sciKMAXX - 1; kopt++)
                {
                        if(a[kopt+1] > a[kopt]*alf[kopt-1][kopt])
                        {
                                break;
                        }
                }
                kmax = kopt;

        } // end eps != epsold || nv != nvold if

        // save starting values
        for(i = 0; i < nRHS; i++)
        {
                ysav[i] = y[i];
        }
        // evaluate the jacobian (at the starting point)
        this->ComputeJacobian(y,dfdy);
        // make sure derivs are current
        this->ComputeDerivatives(y,dydt);

        // a new stepsize or integration; reestablish the order window
        if(m_dTime != xnew || m_dStepSize != m_dNextStep)
        {
                first = 1;
                kopt = kmax;
        }
        reduct = 0;
        for(;;)
        {
                // evaluate the sequence of modified midpoint integrations
                for(k = 0; k <= kmax; k++)
                {
                        xnew = m_dTime + m_dStepSize;
                        // XXX debug
                        //cout << "m_dTime = " << m_dTime << ", m_dStepSize = " << m_dStepSize << ", nstep = " << nseq[k] << endl;
                        //cin.get();
                        if(xnew == m_dTime) {cout << "ERROR : Step size underflow in CStiffBulirschStoerMover::StiffBSStep()" << endl;}
                        this->SemiImplicitMidpoint(nRHS,ysav,dydt,dfdy,nseq[k],yseq);
                        xest = (m_dStepSize/nseq[k])*(m_dStepSize/nseq[k]);
                        this->PolynomialExtrapolation(nRHS,k,xest,yseq,y,yerr,x,d);
                        if(k != 0)
                        {
                                errmax = m_scdTINY;
                                for(i = 0; i < nRHS; i++)
                                {
                                        errmax = __max(errmax,fabs(yerr[i]/m_pdScale[i]));
                                }
                                errmax = errmax/m_dEps;
                                km = k-1;
                                m_pdError[km] = pow(errmax/m_scdSAFE1,1.0/(2*km + 3));
                        }
                        if(k != 0 && (k >= kopt-1 || first))
                        {
                                if(errmax < 1.0)
                                {
                                        exitflag = true;
                                        break;
                                }
                                if(k == kmax || k == kopt+1)
                                {
                                        red = m_scdSAFE2/m_pdError[km];
                                        break;
                                }
                                else if(k == kopt && alf[kopt-1][kopt] < m_pdError[km])
                                {
                                        red = 1.0/m_pdError[km];
                                        break;
                                }
                                else if(kopt == kmax && alf[km][kmax-1] < m_pdError[km])
                                {
                                        red = alf[km][kmax-1]*m_scdSAFE2/m_pdError[km];
                                        break;
                                }
                                else if(alf[km][kopt] < m_pdError[km])
                                {
                                        red = alf[km][kopt-1]/m_pdError[km];
                                        break;
                                }
                        }
                }
                if(exitflag) {break;}
                red = __min(red,m_scdREDMIN);
                red = __max(red,m_scdREDMAX);
                m_dStepSize = m_dStepSize*red;
                reduct = 1;

        }
        m_dTime = xnew;
        m_dLastStep = m_dStepSize;
        first = 0;
        wrkmin = 1.0e35;
        for(kk = 0; kk <= km; kk++)
        {
                fact = __max(m_pdError[kk],m_scdSCALMX);
                work = fact*a[kk+1];
                if(work < wrkmin)
                {
                        scale = fact;
                        wrkmin = work;
                        kopt = kk+1;
                }
        }
        m_dNextStep = m_dStepSize/scale;
        if(kopt >= k && kopt != kmax && !reduct)
        {
                fact = __max(scale/alf[kopt-1][kopt],m_scdSCALMX);
                if(a[kopt+1]*fact <= wrkmin)
                {
                        m_dNextStep = m_dStepSize/fact;
                        kopt++;
                }
        }

        // deletes
        delete [] nseq;
        delete [] a;
        delete [] alf[0];
        delete [] alf;
        delete [] d[0];
        delete [] d;
        delete [] x;
        delete [] yerr;
        delete [] ysav;
        delete [] yseq;
        delete [] dfdy[0];
        delete [] dfdy;
}


void CStiffBulirschStoerMover::SemiImplicitMidpoint(int nRHS, double *y, double *dydt, double **dfdy, int nStep, double *yout)
{
        int i,j,nn;
        double h,x;

        // memory allocation
        double **a = new double *[nRHS];
        a[0] = new double[nRHS*nRHS];
        for(i = 1; i < nRHS; i++)
        {
                a[i] = &(a[0][i*nRHS]);
        }
        double *del = new double[nRHS];
        double *ytemp = new double[nRHS];

        h = m_dStepSize/nStep;

        // set up the 1 - hf' matrix
        for(i = 0; i < nRHS; i++)
        {
                for(j = 0; j < nRHS; j++)
                {
                        a[i][j] = -h*dfdy[i][j];
                }
                ++a[i][i];
        }
        // set up right hand side for first step, use yout[] for temp. storage
        for(i = 0; i < nRHS;i++)
        {
                yout[i] = h*dydt[i];
        }
        // LU backsubstitution with yout as the rhs
        this->LUSolveLinearSystem(nRHS,a,yout);
        // first step
        for(i = 0; i < nRHS; i++)
        {
                ytemp[i] = y[i] + (del[i] = yout[i]);
        }
        x = m_dTime + h;
        // use yout[] for temporary derivative storage
        this->ComputeDerivatives(ytemp,yout);
        // general step
        for(nn = 2; nn <= nStep; nn++)
        {
                for(i = 0; i < nRHS; i++)
                {
                        yout[i] = h*yout[i] - del[i];
                }
                // Do LU Backsubstitution
                this->LUSolveLinearSystem(nRHS,a,yout);
                for(i = 0; i < nRHS; i++)
                {
                        ytemp[i] += (del[i] += 2.0*yout[i]);
                }
                x += h;
                this->ComputeDerivatives(ytemp,yout);
        }
        // set up RHS for last step
        for(i = 0; i < nRHS; i++)
        {
                yout[i] = h*yout[i] - del[i];
        }
        // Do LU backsubstitution
        this->LUSolveLinearSystem(nRHS,a,yout);
        // take last step
        for(i = 0; i < nRHS; i++)
        {
                yout[i] += ytemp[i];
        }

        delete [] a[0];
        delete [] a;
        delete [] del;
        delete [] ytemp;
}

void CStiffBulirschStoerMover::PolynomialExtrapolation(int nRHS, int iest, double xest, double *yest, double *yz, double *dy, double *x, double **d)
{
        int k1,j;
        double q,f2,f1,delta;

        double *c = new double[nRHS];
        x[iest] = xest;
        for(j = 0; j < nRHS; j++)
        {
                dy[j] = yz[j] = yest[j];
        }
        // store first estimate in first column
        if(iest == 0)
        {
                for(j = 0; j < nRHS; j++)
                {
                        d[j][0] = yest[j];
                }
        }
        else
        {
                for(j = 0; j < nRHS; j++)
                {
                        c[j] = yest[j];
                }
                for(k1 = 0; k1 < iest; k1++)
                {
                        delta = 1.0/(x[iest-k1-1] - xest);
                        f1 = xest*delta;
                        f2 = x[iest-k1-1]*delta;
                        for(j = 0; j < nRHS; j++)
                        {
                                q = d[j][k1];
                                d[j][k1] = dy[j];
                                delta = c[j] - q;
                                dy[j] = f1*delta;
                                c[j] = f2*delta;
                                yz[j] += dy[j];
                        }
                }
                for(j = 0; j < nRHS; j++)
                {
                        d[j][iest] = dy[j];
                }
        }
        delete [] c;
}

void CStiffBulirschStoerMover::LUSolveLinearSystem(int nRHS, double **a, double *b)
{
        /*
        // XXX debug
        for(int ii = 0; ii < nRHS; ii++)
        {
                for(int jj = 0; jj < nRHS; jj++)
                {
                        cout << a[ii][jj] << "    " << endl;
                }
                cout << endl;
        }
        cin.get();
        */
        // we're giving the LAPACK routine the upper triangle (though it doesn't matter)
        char UPLO = 'U';
        char TRANS = 'N';
        // can check for correct factorization
        int INFO=0; 
        // order of A
        int N = nRHS;
        int M = nRHS;
        
        // need to be >= max(1,N);
        int LDA = nRHS;
        int LDB = nRHS;
        // only one right hand side
        int NRHS  = 1;

        // allocate temporary storage space
        double *B = new double[nRHS];

        // LU DECOMPOSITION & SOLUTION (LAPACK)
        // copy alpha and beta appropriately
        double *A = new double[nRHS*nRHS];
        int *IPIV = new int[nRHS];
        int counter = 0;
        for(int j = 0; j < nRHS; j++)
        {
                B[j] = b[j];
                for(int i = 0; i < nRHS; i++)
                {
                        A[counter] = a[i][j];
                        counter++;
                }
        }
        
        // LU decompose
        DGETRF(&M,&N,A,&LDA,IPIV,&INFO);
        if(INFO > 0)
        {
                cout << "Matrix singularity suspected - StiffBulirschStoerMover::LUSolveLinearSystem failed . . . " << endl;
        }
        else
        {
                DGETRS(&TRANS,&N,&NRHS,A,&LDA,IPIV,B,&LDB,&INFO);
                // copy back the solution
                for(int i = 0; i < nRHS; i++)
                {
                        b[i] = B[i];
                }
        }

                
        delete [] A;
        delete [] B;
        delete [] IPIV;
}

---