LSODA with diffusion

ksolve/Ksolve.cpp

@@ -659,8 +659,10 @@ void Ksolve::reinit( const Eref& e, ProcPtr p )
 
     if ( isBuilt_ )
     {
-        for ( unsigned int i = 0 ; i < pools_.size(); ++i )
+        for ( unsigned int i = 0 ; i < pools_.size(); ++i ) {
+            pools_[i].setNumVoxels( pools_.size() );
             pools_[i].reinit( p->dt );
+		}
     }
     else
     {

ksolve/VoxelPools.cpp

@@ -36,6 +36,7 @@ using namespace boost::numeric;
 
 VoxelPools::VoxelPools() : pLSODA(nullptr)
 {
+	lsodaState_ = 1;
 #ifdef USE_GSL
     driver_ = 0;
 #endif
@@ -57,6 +58,7 @@ VoxelPools::~VoxelPools()
 void VoxelPools::reinit( double dt )
 {
     VoxelPoolsBase::reinit();
+	lsodaState_ = 1;
 #ifdef USE_GSL
     if ( !driver_ )
         return;
@@ -106,27 +108,19 @@ const string VoxelPools::getMethod( )
 void VoxelPools::advance( const ProcInfo* p )
 {
     double t = p->currTime - p->dt;
-	if ( p->isStart() ) // True if first step or restart.
-		lsodaState = 1;
-	// lsodaState = 1;
-	/**
-	if (t < p->dt )
-		lsodaState = 1;
-	else
-		lsodaState = 1; // Earlier this was set to 3. This made it 
-		// impossible to alter non-buffered molecules from the script after
-		// the simulation had started, even if the control returned to the
-		// parser.
-	**/
     Ksolve* k = reinterpret_cast<Ksolve*>( stoichPtr_->getKsolve().eref().data() );
 
     if( getMethod() == "lsoda" )
     {
+		// True if first step or restart, or if diffusion. Tells LSODA to 
+		// recalculate using new pool n values, which slows it down a bit. 
+		if ( p->isStart() || (numVoxels_ > 1) ) 
+			lsodaState_ = 1;
     	size_t totVar = stoichPtr_->getNumVarPools() + stoichPtr_->getNumProxyPools();
         vector<double> yout(size()+1);
         pLSODA->lsoda_update( &VoxelPools::lsodaSys, size()
                 , Svec(), yout , &t
-                , p->currTime, &lsodaState, this
+                , p->currTime, &lsodaState_, this
                 );
 
         // Now update the y from yout. This is different thant normal GSL or
@@ -135,7 +129,7 @@ void VoxelPools::advance( const ProcInfo* p )
         for (size_t i = 0; i < totVar; i++)
             varS()[i] = yout[i+1];
 
-        if( lsodaState == 0 )
+        if( lsodaState_ == 0 )
         {
             cerr << "Error: VoxelPools::advance: LSODA integration error at time "
                  << t << "\n";

ksolve/VoxelPools.h

@@ -101,7 +101,7 @@ class VoxelPools: public VoxelPoolsBase
 
     std::shared_ptr<LSODA> pLSODA;
     LSODA_ODE_SYSTEM_TYPE lsodaSystem;
-    int lsodaState = 1;
+    int lsodaState_;
 
 #ifdef USE_GSL
     gsl_odeiv2_driver* driver_;

ksolve/VoxelPoolsBase.cpp

@@ -139,6 +139,11 @@ void VoxelPoolsBase::scaleVolsBufsRates(double ratio, const Stoich* stoichPtr)
     }
 }
 
+void VoxelPoolsBase::setNumVoxels( unsigned int n )
+{
+	numVoxels_ = n;
+}
+
 //////////////////////////////////////////////////////////////
 // Zombie Pool Access functions
 //////////////////////////////////////////////////////////////

ksolve/VoxelPoolsBase.h

@@ -205,12 +205,19 @@ class VoxelPoolsBase
 
     void scaleVolsBufsRates( double ratio, const Stoich* stoichPtr );
 
+	void setNumVoxels( unsigned int );
+
     /// Debugging utility
     void print() const;
 
 protected:
     const Stoich* stoichPtr_;
     vector< RateTerm* > rates_;
+	/**
+	 * Number of voxels. If > 1, set flag for LSODA to handle molecule
+	 * flux on each timestep during diffusion, which slows it down.
+	 */
+	unsigned int numVoxels_;
 
 private:
     /**

python/rdesigneur/rdesigneur.py

@@ -133,7 +133,7 @@ def __init__(self,
             moogList = [], 
             outputFileList = [], # List of all file save specifications.
             modelFileNameList = [], # List of any files used to build.
-            ode_method = "gsl",  # gsl, lsoda, gssa, gillespie
+            ode_method = "lsoda",  # gsl, lsoda, gssa, gillespie
             isLegacyMethod = False,
             params = None
         ):
@@ -1588,10 +1588,12 @@ def _buildNeuroMesh( self ):
         self.chemid.name = 'temp_chem'
         newChemid = moose.Neutral( self.model.path + '/chem' )
         comptlist = self._assignComptNamesFromKkit_SBML()
-        if len( comptlist ) == 1 and comptlist[0].name == 'kinetics':
+        #if len( comptlist ) == 1 and comptlist[0].name == 'kinetics':
+        if len( comptlist ) == 1:
             comptlist[0].name = 'dend'
         comptdict = { i.name:i for i in comptlist }
         if len(comptdict) == 1 or 'dend' in comptdict:
+            #print( "COMPTDICT = ", comptdict )
             self.dendCompt = moose.NeuroMesh( newChemid.path + '/dend' )
             self.dendCompt.geometryPolicy = 'cylinder'
             self.dendCompt.separateSpines = 0