Source code for switchKineticSolvers

#########################################################################
## This program is part of 'MOOSE', the
## Messaging Object Oriented Simulation Environment.
##           Copyright (C) 2014 Upinder S. Bhalla. and NCBS
## It is made available under the terms of the
## GNU Lesser General Public License version 2.1
## See the file COPYING.LIB for the full notice.
#########################################################################

import moose
import pylab
import numpy
import matplotlib.pyplot as plt
import sys


def runAndSavePlots( name ):
    runtime = 20.0
    moose.reinit()
    moose.start( runtime )
    pa = moose.Neutral( '/model/graphs/' + name )
    for x in moose.wildcardFind( '/model/#graphs/conc#/#' ):
        if ( x.tick != -1 ):
            tabname = '/model/graphs/' + name + '/' + x.name + '.' + name
            y = moose.Table( tabname )
            y.vector = x.vector
            y.tick = -1

# Takes args ee, gsl, or gssa
def switchSolvers( solver ):
        if ( moose.exists( 'model/kinetics/stoich' ) ):
            moose.delete( '/model/kinetics/stoich' )
            moose.delete( '/model/kinetics/ksolve' )
        compt = moose.element( '/model/kinetics' )
        if ( solver == 'gsl' ):
            ksolve = moose.Ksolve( '/model/kinetics/ksolve' )
        if ( solver == 'gssa' ):
            ksolve = moose.Gsolve( '/model/kinetics/ksolve' )
        if ( solver != 'ee' ):
            stoich = moose.Stoich( '/model/kinetics/stoich' )
            stoich.compartment = compt
            stoich.ksolve = ksolve
            stoich.path = "/model/kinetics/##"

[docs]def main(): """ At zero order, you can select the solver you want to use within the function moose.loadModel( filename, modelpath, solver ). Having loaded in the model, you can change the solver to use on it. This example illustrates how to assign and change solvers for a kinetic model. This process is necessary in two situations: * If we want to change the numerical method employed, for example, from deterministic to stochastic. * If we are already using a solver, and we have changed the reaction network by adding or removing molecules or reactions. Note that we do not have to change the solvers if the volume or reaction rates change. In this example the model is loaded in with a gsl solver. The sequence of solver calculations is: #. gsl #. ee #. gsl #. gssa #. gsl If you're removing the solvers, you just delete the stoichiometry object and the associated ksolve/gsolve. Should there be diffusion (a dsolve)then you should delete that too. If you're building the solvers up again, then you must do the following steps in order: #. build up the ksolve/gsolve and stoich (any order) #. Assign stoich.ksolve #. Assign stoich.path. See the Reaction-diffusion section should you want to do diffusion as well. """ solver = "gsl" # Pick any of gsl, gssa, ee.. mfile = '../genesis/kkit_objects_example.g' modelId = moose.loadModel( mfile, 'model', solver ) # Increase volume so that the stochastic solver gssa # gives an interesting output compt = moose.element( '/model/kinetics' ) compt.volume = 1e-19 runAndSavePlots( 'gsl' ) ######################################################### switchSolvers( 'ee' ) runAndSavePlots( 'ee' ) ######################################################### switchSolvers( 'gsl' ) runAndSavePlots( 'gsl2' ) ######################################################### switchSolvers( 'gssa' ) runAndSavePlots( 'gssa' ) ######################################################### switchSolvers( 'gsl' ) runAndSavePlots( 'gsl3' ) ######################################################### # Display all plots. fig = plt.figure( figsize = (12, 10) ) orig = fig.add_subplot( 511 ) gsl = fig.add_subplot( 512 ) ee = fig.add_subplot( 513 ) gsl2 = fig.add_subplot( 514 ) gssa = fig.add_subplot( 515 ) plotdt = moose.element( '/clock' ).tickDt[18] for x in moose.wildcardFind( '/model/#graphs/conc#/#' ): t = numpy.arange( 0, x.vector.size, 1 ) * plotdt orig.plot( t, x.vector, label=x.name ) for x in moose.wildcardFind( '/model/graphs/gsl/#' ): t = numpy.arange( 0, x.vector.size, 1 ) * plotdt gsl.plot( t, x.vector, label=x.name ) for x in moose.wildcardFind( '/model/graphs/ee/#' ): t = numpy.arange( 0, x.vector.size, 1 ) * plotdt ee.plot( t, x.vector, label=x.name ) for x in moose.wildcardFind( '/model/graphs/gsl2/#' ): t = numpy.arange( 0, x.vector.size, 1 ) * plotdt gsl2.plot( t, x.vector, label=x.name ) for x in moose.wildcardFind( '/model/graphs/gssa/#' ): t = numpy.arange( 0, x.vector.size, 1 ) * plotdt gssa.plot( t, x.vector, label=x.name ) plt.legend() pylab.show() quit()
# Run the 'main' if this script is executed standalone. if __name__ == '__main__': main()