VoxelPools Class Reference

#include <VoxelPools.h>

Inheritance diagram for VoxelPools:

Collaboration diagram for VoxelPools:

Public Member Functions
void	advance (const ProcInfo *p)
	Do the numerical integration. Advance the simulation. More...
const Stoich *	getStoich ()
void	print () const
	Used for debugging. More...
void	reinit (double dt)
void	setInitDt (double dt)
	Set initial timestep to use by the solver. More...
void	setStoich (Stoich *stoich, const OdeSystem *ode)
void	setVolumeAndDependencies (double vol)
	Handles volume change and subsequent cascading updates. More...
void	updateAllRateTerms (const vector< RateTerm * > &rates, unsigned int numCoreRates)
	Updates all the rate constants from the reference rates vector. More...
void	updateRates (const double *s, double *yprime) const
void	updateRateTerms (const vector< RateTerm * > &rates, unsigned int numCoreRates, unsigned int index)
void	updateReacVelocities (const double *s, vector< double > &v) const
	VoxelPools ()
virtual	~VoxelPools ()
Public Member Functions inherited from VoxelPoolsBase
void	addProxyTransferIndex (unsigned int comptIndex, unsigned int transferIndex)
void	addProxyVoxy (unsigned int comptIndex, Id comptId, unsigned int voxel)
	Adds the index of a voxel to which proxy data should go. More...
void	backwardReacVolumeFactor (unsigned int i, double volume)
void	filterCrossRateTerms (const vector< Id > &offSolverReacs, const vector< pair< Id, Id > > &offSolverReacCompts)
void	forwardReacVolumeFactor (unsigned int i, double volume)
double	getDiffConst (unsigned int) const
double	getN (unsigned int) const
double	getNinit (unsigned int) const
double	getVolume () const
	Return the volume of the voxel. More...
double	getXreacScaleProducts (unsigned int i) const
double	getXreacScaleSubstrates (unsigned int i) const
bool	hasXfer (unsigned int comptIndex) const
	True when this voxel has data to be transferred. More...
bool	isVoxelJunctionPresent (Id i1, Id i2) const
void	print () const
	Debugging utility. More...
void	reinit ()
void	resetXreacScale (unsigned int size)
void	resizeArrays (unsigned int totNumPools)
	Allocate # of pools. More...
const double *	S () const
void	scaleVolsBufsRates (double ratio, const Stoich *stoichPtr)
void	setDiffConst (unsigned int, double v)
void	setN (unsigned int i, double v)
void	setNinit (unsigned int, double v)
void	setVolume (double vol)
	Just assigns the volume without any cascading to other values. More...
const double *	Sinit () const
unsigned int	size () const
vector< double > &	Svec ()
	Returns a handle to the mol # vector. More...
double *	varS ()
double *	varSinit ()
	VoxelPoolsBase ()
void	xferIn (const vector< unsigned int > &poolIndex, const vector< double > &values, const vector< double > &lastValues, unsigned int voxelIndex)
void	xferInOnlyProxies (const vector< unsigned int > &poolIndex, const vector< double > &values, unsigned int numProxyPools, unsigned int voxelIndex)
void	xferOut (unsigned int voxelIndex, vector< double > &values, const vector< unsigned int > &poolIndex)
	Assembles data values for sending out for x-compt reacs. More...
virtual	~VoxelPoolsBase ()

Additional Inherited Members
Protected Attributes inherited from VoxelPoolsBase
vector< RateTerm * >	rates_
const Stoich *	stoichPtr_

Detailed Description

This is the class for handling reac-diff voxels used for deterministic computations.

Definition at line 29 of file VoxelPools.h.

Constructor & Destructor Documentation

VoxelPools::VoxelPools

(

)

Definition at line 36 of file VoxelPools.cpp.

{
#ifdef USE_GSL
    driver_ = 0;
#endif
}

VoxelPools::~VoxelPools ( )

virtual

Definition at line 43 of file VoxelPools.cpp.

{
    for ( unsigned int i = 0; i < rates_.size(); ++i )
        delete( rates_[i] );
#ifdef USE_GSL
    if ( driver_ )
        gsl_odeiv2_driver_free( driver_ );
#endif
}

Member Function Documentation

void VoxelPools::advance

(

const ProcInfo *

p

)

Do the numerical integration. Advance the simulation.

Definition at line 89 of file VoxelPools.cpp.

References ProcInfo::currTime, ProcInfo::dt, VoxelPoolsBase::stoichPtr_, and Stoich::updateFuncs().

{
    double t = p->currTime - p->dt;
#ifdef USE_GSL
    int status = gsl_odeiv2_driver_apply( driver_, &t, p->currTime, varS());
    if ( status != GSL_SUCCESS )
    {
        cout << "Error: VoxelPools::advance: GSL integration error at time "
             << t << "\n";
        cout << "Error info: " << status << ", " <<
             gsl_strerror( status ) << endl;
        if ( status == GSL_EMAXITER )
            cout << "Max number of steps exceeded\n";
        else if ( status == GSL_ENOPROG )
            cout << "Timestep has gotten too small\n";
        else if ( status == GSL_EBADFUNC )
            cout << "Internal error\n";
        assert( 0 );
    }

#elif USE_BOOST_ODE


    // NOTE: Make sure to assing vp to BoostSys vp. In next call, it will be used by
    // updateRates func. Unlike gsl call, we can't pass extra void*  to gslFunc.
    VoxelPools* vp = reinterpret_cast< VoxelPools* >( sys_.params );
    sys_.vp = vp;
    /*-----------------------------------------------------------------------------
    NOTE: 04/21/2016 11:31:42 AM

    We need to call updateFuncs  here (unlike in GSL solver) because there
    is no way we can update const vector_type_& y in evalRatesUsingBoost
    function. In gsl implmentation one could do it, because const_cast can
    take away the constantness of double*. This probably makes the call bit
    cleaner.
     *-----------------------------------------------------------------------------*/
    vp->stoichPtr_->updateFuncs( &Svec()[0], p->currTime );

    /*-----------------------------------------------------------------------------
     * Using integrate function works with with default stepper type.
     *
     *  NOTICE to developer:
     *  If you are planning your own custom typdedef of stepper_type_ (see
     *  file BoostSystem.h), the you may run into troble. Have a look at this
     *  http://boostw.boost.org/doc/libs/1_56_0/boost/numeric/odeint/integrate/integrate.hpp
     *-----------------------------------------------------------------------------
     */

    double absTol = sys_.epsAbs;
    double relTol = sys_.epsRel;


    const double fixedDt = 0.1;

    if( sys_.method == "rk2" )
        odeint::integrate_const( rk_midpoint_stepper_type_()
                                 , sys_ , Svec()
                                 , p->currTime - p->dt, p->currTime, std::min( p->dt, fixedDt )
                               );
    else if( sys_.method == "rk4" )
        odeint::integrate_const( rk4_stepper_type_()
                                 , sys_ , Svec()
                                 , p->currTime - p->dt, p->currTime, std::min( p->dt, fixedDt )
                               );
    else if( sys_.method == "rk5")
        odeint::integrate_const( rk_karp_stepper_type_()
                                 , sys_ , Svec()
                                 , p->currTime - p->dt, p->currTime, std::min( p->dt, fixedDt )
                               );
    else if( sys_.method == "rk5a")
        odeint::integrate_adaptive(
            odeint::make_controlled<rk_karp_stepper_type_>( absTol, relTol)
            , sys_
            , Svec()
            , p->currTime - p->dt
            , p->currTime
            , p->dt
        );
    else if ("rk54" == sys_.method )
        odeint::integrate_const( rk_karp_stepper_type_()
                                 , sys_ , Svec()
                                 , p->currTime - p->dt, p->currTime, std::min( p->dt, fixedDt )
                               );
    else if ("rk54a" == sys_.method )
        odeint::integrate_adaptive(
            odeint::make_controlled<rk_karp_stepper_type_>( absTol, relTol )
            , sys_, Svec()
            , p->currTime - p->dt
            , p->currTime
            , p->dt
        );
    else if ("rk5" == sys_.method )
        odeint::integrate_const( rk_dopri_stepper_type_()
                                 , sys_ , Svec()
                                 , p->currTime - p->dt, p->currTime, std::min( p->dt, fixedDt )
                               );
    else if ("rk5a" == sys_.method )
        odeint::integrate_adaptive(
            odeint::make_controlled<rk_dopri_stepper_type_>( absTol, relTol )
            , sys_, Svec()
            , p->currTime - p->dt
            , p->currTime
            , p->dt
        );
    else if( sys_.method == "rk8" )
        odeint::integrate_const( rk_felhberg_stepper_type_()
                                 , sys_ , Svec()
                                 , p->currTime - p->dt, p->currTime, std::min( p->dt, fixedDt )
                               );
    else if( sys_.method == "rk8a" )
        odeint::integrate_adaptive(
            odeint::make_controlled<rk_felhberg_stepper_type_>( absTol, relTol )
            , sys_, Svec()
            , p->currTime - p->dt
            , p->currTime
            , p->dt
        );

    else
        odeint::integrate_adaptive(
            odeint::make_controlled<rk_karp_stepper_type_>( absTol, relTol )
            , sys_, Svec()
            , p->currTime - p->dt
            , p->currTime
            , p->dt
        );
#endif
    if ( !stoichPtr_->getAllowNegative() )   // clean out negatives
    {
        unsigned int nv = stoichPtr_->getNumVarPools();
        double* vs = varS();
        for ( unsigned int i = 0; i < nv; ++i )
        {
            if ( signbit(vs[i]) )
                vs[i] = 0.0;
        }
    }
}

Here is the call graph for this function:

const Stoich* VoxelPools::getStoich

(

)

void VoxelPools::print

(

)

const

Used for debugging.

For debugging: Print contents of voxel pool.

Definition at line 369 of file VoxelPools.cpp.

References VoxelPoolsBase::print().

{
    cout << "numAllRates = " << rates_.size() <<
         ", numLocalRates= " << stoichPtr_->getNumCoreRates() << endl;
    VoxelPoolsBase::print();
}

Here is the call graph for this function:

void VoxelPools::reinit

(

double

dt

)

Definition at line 56 of file VoxelPools.cpp.

References VoxelPoolsBase::reinit().

{
    VoxelPoolsBase::reinit();
#ifdef USE_GSL
    if ( !driver_ )
        return;
    gsl_odeiv2_driver_reset( driver_ );
    gsl_odeiv2_driver_reset_hstart( driver_, dt / 10.0 );
#endif
}

Here is the call graph for this function:

void VoxelPools::setInitDt

(

double

dt

)

Set initial timestep to use by the solver.

Definition at line 234 of file VoxelPools.cpp.

{
#ifdef USE_GSL
    gsl_odeiv2_driver_reset_hstart( driver_, dt );
#endif
}

void VoxelPools::setStoich	(	Stoich *	stoich,
		const OdeSystem *	ode
	)

setStoich: Assigns the ODE system and Stoich. Stoich may be modified to create a new RateTerm vector in case the volume of this VoxelPools is new.

Definition at line 67 of file VoxelPools.cpp.

References OdeSystem::epsAbs, OdeSystem::epsRel, OdeSystem::initStepSize, and VoxelPoolsBase::reinit().

Referenced by SteadyState::setStoich().

{
    stoichPtr_ = s;
#ifdef USE_GSL
    if ( ode )
    {
        sys_ = ode->gslSys;
        if ( driver_ )
            gsl_odeiv2_driver_free( driver_ );

        driver_ = gsl_odeiv2_driver_alloc_y_new(
                      &sys_, ode->gslStep, ode->initStepSize,
                      ode->epsAbs, ode->epsRel 
                );
    }
#elif USE_BOOST_ODE
    if( ode )
        sys_ = ode->boostSys;
#endif
    VoxelPoolsBase::reinit();
}

Here is the call graph for this function:

Here is the caller graph for this function:

void VoxelPools::setVolumeAndDependencies ( double  vol )

virtual

Handles volume change and subsequent cascading updates.

Handle volume updates. Inherited Virtual func.

Reimplemented from VoxelPoolsBase.

Definition at line 380 of file VoxelPools.cpp.

References VoxelPoolsBase::setVolumeAndDependencies().

{
    VoxelPoolsBase::setVolumeAndDependencies( vol );
    updateAllRateTerms( stoichPtr_->getRateTerms(),
                        stoichPtr_->getNumCoreRates() );
}

Here is the call graph for this function:

void VoxelPools::updateAllRateTerms ( const vector< RateTerm * > &  rates, unsigned int  numCoreRates  )

virtual

Updates all the rate constants from the reference rates vector.

Implements VoxelPoolsBase.

Definition at line 282 of file VoxelPools.cpp.

Referenced by SteadyState::setStoich().

{
    // Clear out old rates if any
    for ( unsigned int i = 0; i < rates_.size(); ++i )
        delete( rates_[i] );

    rates_.resize( rates.size() );
    for ( unsigned int i = 0; i < numCoreRates; ++i )
        rates_[i] = rates[i]->copyWithVolScaling( getVolume(), 1, 1 );
    for ( unsigned int i = numCoreRates; i < rates.size(); ++i )
    {
        rates_[i] = rates[i]->copyWithVolScaling(  getVolume(),
                    getXreacScaleSubstrates(i - numCoreRates),
                    getXreacScaleProducts(i - numCoreRates ) );
    }
}

Here is the caller graph for this function:

void VoxelPools::updateRates	(	const double *	s,
		double *	yprime
	)		const

Core computation function. Updates the reaction velocities vector yprime given the current mol 'n' vector s.

Definition at line 318 of file VoxelPools.cpp.

References KinSparseMatrix::computeRowRate(), SparseMatrix< T >::nColumns(), and SparseMatrix< T >::nRows().

Referenced by SteadyState::classifyState().

{
    const KinSparseMatrix& N = stoichPtr_->getStoichiometryMatrix();
    vector< double > v( N.nColumns(), 0.0 );
    vector< double >::iterator j = v.begin();
    // totVar should include proxyPools only if this voxel uses them
    unsigned int totVar = stoichPtr_->getNumVarPools() +
                          stoichPtr_->getNumProxyPools();
    // totVar should include proxyPools if this voxel does not use them
    unsigned int totInvar = stoichPtr_->getNumBufPools();
    assert( N.nColumns() == 0 ||
            N.nRows() == stoichPtr_->getNumAllPools() );
    assert( N.nColumns() == rates_.size() );

    for ( vector< RateTerm* >::const_iterator
            i = rates_.begin(); i != rates_.end(); i++)
    {
        *j++ = (**i)( s );
        assert( !std::isnan( *( j-1 ) ) );
    }

    for (unsigned int i = 0; i < totVar; ++i)
        *yprime++ = N.computeRowRate( i , v );
    for (unsigned int i = 0; i < totInvar ; ++i)
        *yprime++ = 0.0;
}

Here is the call graph for this function:

Here is the caller graph for this function:

void VoxelPools::updateRateTerms ( const vector< RateTerm * > &  rates, unsigned int  numCoreRates, unsigned int  index  )

virtual

updateRateTerms updates the rate consts of a belonging to the specified index on the rates vector. It does recaling and assigning using values from the internal rates vector.

Implements VoxelPoolsBase.

Definition at line 300 of file VoxelPools.cpp.

{
    // During setup or expansion of the reac system, it is possible to
    // call this function before the rates_ term is assigned. Disable.
    if ( index >= rates_.size() )
        return;
    delete( rates_[index] );
    if ( index >= numCoreRates )
        rates_[index] = rates[index]->copyWithVolScaling(
                            getVolume(),
                            getXreacScaleSubstrates(index - numCoreRates),
                            getXreacScaleProducts(index - numCoreRates ) );
    else
        rates_[index] = rates[index]->copyWithVolScaling(
                            getVolume(), 1.0, 1.0 );
}

void VoxelPools::updateReacVelocities	(	const double *	s,
		vector< double > &	v
	)		const

updateReacVelocities computes the velocity v of each reaction from the vector s of pool #s. This is a utility function for programs like SteadyState that need to analyze velocity.

updateReacVelocities computes the velocity v of each reaction. This is a utility function for programs like SteadyState that need to analyze velocity.

Definition at line 350 of file VoxelPools.cpp.

References SparseMatrix< T >::nColumns().

{
    const KinSparseMatrix& N = stoichPtr_->getStoichiometryMatrix();
    assert( N.nColumns() == rates_.size() );

    vector< RateTerm* >::const_iterator i;
    v.clear();
    v.resize( rates_.size(), 0.0 );
    vector< double >::iterator j = v.begin();

    for ( i = rates_.begin(); i != rates_.end(); i++)
    {
        *j++ = (**i)( s );
        assert( !std::isnan( *( j-1 ) ) );
    }
}

Here is the call graph for this function:

---

The documentation for this class was generated from the following files:

• moose-core/ksolve/VoxelPools.h
• moose-core/ksolve/VoxelPools.cpp