testKsolve.cpp File Reference

#include "header.h"
#include "../shell/Shell.h"
#include "RateTerm.h"
#include "muParser.h"
#include "FuncTerm.h"
#include "SparseMatrix.h"
#include "KinSparseMatrix.h"
#include "VoxelPoolsBase.h"
#include "../mesh/VoxelJunction.h"
#include "XferInfo.h"
#include "ZombiePoolInterface.h"
#include "Stoich.h"

Include dependency graph for testKsolve.cpp:

Go to the source code of this file.

Functions
Id	makeReacTest ()
void	testBuildStoich ()
void	testFuncTerm ()
void	testKsolve ()
void	testKsolveProcess ()
void	testRunGsolve ()
void	testRunKsolve ()
void	testSetupReac ()

Function Documentation

Id makeReacTest

(

)

Tab controlled by table A + Tab <===> B A + B -sumtot–> tot1 2B <===> C

C —e1Pool —> D D —e2Pool -—> E

All these are created on /kinetics, a cube compartment of vol 1e-18 m^3

Definition at line 34 of file testKsolve.cpp.

References Eref::data(), Shell::doAddMsg(), Shell::doCreate(), Shell::doSetClock(), Id::eref(), Field< A >::get(), NA, PI, Field< A >::set(), and Id::value().

Referenced by testBuildStoich(), testRunGsolve(), testRunKsolve(), testSetupReac(), and testVolScaling().

{
    double simDt = 0.1;
    double plotDt = 0.1;
    Shell* s = reinterpret_cast< Shell* >( Id().eref().data() );

    Id pools[10];
    unsigned int i = 0;
    // Make the objects.
    Id kin = s->doCreate( "CubeMesh", Id(), "kinetics", 1 );
    Id tab = s->doCreate( "StimulusTable", kin, "tab", 1 );
    Id T = pools[i++] = s->doCreate( "BufPool", kin, "T", 1 );
    Id A = pools[i++] = s->doCreate( "Pool", kin, "A", 1 );
    Id B = pools[i++] = s->doCreate( "Pool", kin, "B", 1 );
    Id C = pools[i++] = s->doCreate( "Pool", kin, "C", 1 );
    Id D = pools[i++] = s->doCreate( "Pool", kin, "D", 1 );
    Id E = pools[i++] = s->doCreate( "Pool", kin, "E", 1 );
    Id tot1 = pools[i++] = s->doCreate( "BufPool", kin, "tot1", 1 );
    Id sum = s->doCreate( "Function", tot1, "func", 1 ); // Silly that it has to have this name.
    Id sumInput( sum.value() + 1 );
    Id e1Pool = s->doCreate( "Pool", kin, "e1Pool", 1 );
    Id e2Pool = s->doCreate( "Pool", kin, "e2Pool", 1 );
    Id e1 = s->doCreate( "Enz", e1Pool, "e1", 1 );
    Id cplx = s->doCreate( "Pool", e1, "cplx", 1 );
    Id e2 = s->doCreate( "MMenz", e2Pool, "e2", 1 );
    Id r1 = s->doCreate( "Reac", kin, "r1", 1 );
    Id r2 = s->doCreate( "Reac", kin, "r2", 1 );
    Id plots = s->doCreate( "Table2", kin, "plots", 7 );

    // Connect them up
    s->doAddMsg( "Single", tab, "output", T, "setN" );

    s->doAddMsg( "Single", r1, "sub", T, "reac" );
    s->doAddMsg( "Single", r1, "sub", A, "reac" );
    s->doAddMsg( "Single", r1, "prd", B, "reac" );

    Field< unsigned int >::set( sum, "numVars", 2 );
    s->doAddMsg( "Single", A, "nOut", ObjId( sumInput, 0, 0 ), "input" );
    s->doAddMsg( "Single", B, "nOut", ObjId( sumInput, 0, 1 ), "input" );
    s->doAddMsg( "Single", sum, "valueOut", tot1, "setN" );

    s->doAddMsg( "Single", r2, "sub", B, "reac" );
    s->doAddMsg( "Single", r2, "sub", B, "reac" );
    s->doAddMsg( "Single", r2, "prd", C, "reac" );

    s->doAddMsg( "Single", e1, "sub", C, "reac" );
    s->doAddMsg( "Single", e1, "enz", e1Pool, "reac" );
    s->doAddMsg( "Single", e1, "cplx", cplx, "reac" );
    s->doAddMsg( "Single", e1, "prd", D, "reac" );

    s->doAddMsg( "Single", e2, "sub", D, "reac" );
    s->doAddMsg( "Single", e2Pool, "nOut", e2, "enzDest" );
    s->doAddMsg( "Single", e2, "prd", E, "reac" );

    // Set parameters.
    Field< double >::set( A, "concInit", 2 );
    Field< double >::set( e1Pool, "concInit", 1 );
    Field< double >::set( e2Pool, "concInit", 1 );
    Field< string >::set( sum, "expr", "x0+x1" );
    Field< double >::set( r1, "Kf", 0.2 );
    Field< double >::set( r1, "Kb", 0.1 );
    Field< double >::set( r2, "Kf", 0.1 );
    Field< double >::set( r2, "Kb", 0.0 );
    Field< double >::set( e1, "Km", 5 );
    Field< double >::set( e1, "kcat", 1 );
    Field< double >::set( e1, "ratio", 4 );
    Field< double >::set( e2, "Km", 5 );
    Field< double >::set( e2, "kcat", 1 );
    vector< double > stim( 100, 0.0 );
    double vol = Field< double >::get( kin, "volume" );
    for ( unsigned int i = 0; i< 100; ++i ) {
        stim[i] = vol * NA * (1.0 + sin( i * 2.0 * PI / 100.0 ) );
    }
    Field< vector< double > >::set( tab, "vector", stim );
    Field< double >::set( tab, "stepSize", 0.0 );
    Field< double >::set( tab, "stopTime", 10.0 );
    Field< double >::set( tab, "loopTime", 10.0 );
    Field< bool >::set( tab, "doLoop", true );


    // Connect outputs
    for ( unsigned int i = 0; i < 7; ++i )
        s->doAddMsg( "Single", ObjId( plots,i),
                        "requestOut", pools[i], "getConc" );

    // Schedule it.
    for ( unsigned int i = 11; i < 18; ++i )
        s->doSetClock( i, simDt );
    s->doSetClock( 18, plotDt );
    /*
    s->doUseClock( "/kinetics/##[ISA=Reac],/kinetics/##[ISA=EnzBase],/kinetics/##[ISA=SumFunc]",
                    "process", 4 );
    s->doUseClock( "/kinetics/##[ISA=PoolBase]", "process", 5 );
    s->doUseClock( "/kinetics/##[ISA=StimulusTable]",
                    "process", 4 );
    s->doUseClock( "/kinetics/##[ISA=Table]", "process", 8 );
    s->doSetClock( 4, simDt );
    s->doSetClock( 5, simDt );
    s->doSetClock( 8, plotDt );
    */
    return kin;
}

Here is the call graph for this function:

Here is the caller graph for this function:

void testBuildStoich

(

)

Definition at line 156 of file testKsolve.cpp.

References Eref::data(), Shell::doCreate(), Shell::doDelete(), Id::eref(), Field< A >::get(), makeReacTest(), and Field< A >::set().

Referenced by testKsolve().

{
        // In the updated version, we have reordered all pools and reacs by
        // Id number, modulo varPools then bufPOols.
        // Matrix looks like:
        // Reac Name    R1  R2  e1a e1b e2
        // MolName
        // D            -1  0   0   0   0
        // A            -1  0   0   0   0
        // B            +1  -2  0   0   0
        // C            0   +1  -1  0   0
        // enz1         0   0   -1  +1  0
        // e1cplx       0   0   +1  -1  0
        // E            0   0   0   +1  -1
        // F            0   0   0   0   +1
        // enz2         0   0   0   0   0
        // tot1         0   0   0   0   0
        //
        // This has been shuffled to:
        // A            -1  0   0   0   0
        // B            +1  -2  0   0   0
        // C            0   +1  -1  0   0
        // E            0   0   0   +1  -1
        // F            0   0   0   0   +1
        // enz1         0   0   -1  +1  0
        // enz2         0   0   0   0   0
        // e1cplx       0   0   +1  -1  0
        // D            -1  0   0   0   0
        // tot1         0   0   0   0   0
        //
        // (This is the output of the print command on the sparse matrix.)
        //
    Shell* s = reinterpret_cast< Shell* >( Id().eref().data() );
    Id kin = makeReacTest();
    Id ksolve = s->doCreate( "Ksolve", kin, "ksolve", 1 );
    Id stoich = s->doCreate( "Stoich", ksolve, "stoich", 1 );
    Field< Id >::set( stoich, "compartment", kin );
    Field< Id >::set( stoich, "ksolve", ksolve );

    // Used to get at the stoich matrix from gdb.
#ifndef NDEBUG
    Stoich* stoichPtr = reinterpret_cast< Stoich* >( stoich.eref().data() );
#endif

    Field< string >::set( stoich, "path", "/kinetics/##" );

    unsigned int n = Field< unsigned int >::get( stoich, "numAllPools" );
    assert( n == 10 );
    unsigned int r = Field< unsigned int >::get( stoich, "numRates" );
    assert( r == 5 ); // One each for reacs and MMenz, two for Enz.

    vector< int > entries = Field< vector< int > >::get( stoich,
                    "matrixEntry" );
    vector< unsigned int > colIndex = Field< vector< unsigned int > >::get(
                    stoich, "columnIndex" );
    vector< unsigned int > rowStart = Field< vector< unsigned int > >::get(
                    stoich, "rowStart" );

    assert( rowStart.size() == n + 1 );
    assert( entries.size() == colIndex.size() );
    assert( entries.size() == 13 );
    assert( entries[0] == -1 );
    assert( entries[1] == 1 );
    assert( entries[2] == -2 );
    assert( entries[3] == 1 );
    assert( entries[4] == -1 );
    assert( entries[5] == 1 );
    assert( entries[6] == -1 );
    assert( entries[7] == 1 );
    assert( entries[8] == -1 );
    assert( entries[9] == 1 );
    assert( entries[10] == 1 );
    assert( entries[11] == -1 );
    assert( entries[12] == -1 );

    s->doDelete( kin );
    cout << "." << flush;
}

Here is the call graph for this function:

Here is the caller graph for this function:

void testFuncTerm

(

)

Definition at line 305 of file testKsolve.cpp.

References doubleEq(), FuncTerm::setExpr(), and FuncTerm::setReactantIndex().

Referenced by testKsolve().

{
    FuncTerm ft;
    ft.setExpr( "x0 + x1 * t" );
    double args[] = { 1, 2, 3, 4, 5, 6, 7, 8, 9, 10 };

    // First check that it doesn't die even if we forget to set up anything.
    double ans = ft( args, 2.0 );

    vector< unsigned int > mol( 2, 0 );
    mol[0] = 2;
    mol[1] = 0;
    ft.setReactantIndex( mol );

    ans = ft( args, 10.0 );
    assert( doubleEq( ans, 13.0 ) );
    mol[0] = 0;
    mol[1] = 9;
    ft.setReactantIndex( mol );
    ans = ft( args, 2.0 );
    assert( doubleEq( ans, 21.0 ) );
    cout << "." << flush;
}

Here is the call graph for this function:

Here is the caller graph for this function:

void testKsolve

(

)

Definition at line 329 of file testKsolve.cpp.

References testBuildStoich(), testFuncTerm(), testRunGsolve(), testRunKsolve(), and testSetupReac().

Referenced by nonMpiTests().

{
    testSetupReac();
    testBuildStoich();
    testRunKsolve();
    testRunGsolve();
    testFuncTerm();
}

Here is the call graph for this function:

Here is the caller graph for this function:

void testKsolveProcess

(

)

Definition at line 338 of file testKsolve.cpp.

{
;
}

void testRunGsolve

(

)

Definition at line 262 of file testKsolve.cpp.

References Eref::data(), Shell::doCreate(), Shell::doDelete(), Shell::doReinit(), Shell::doSetClock(), Shell::doStart(), Shell::doUseClock(), Id::eref(), makeReacTest(), NA, PI, Field< A >::set(), and SetGet2< A1, A2 >::set().

Referenced by testKsolve().

{
    double simDt = 0.1;
    // double plotDt = 0.1;
    Shell* s = reinterpret_cast< Shell* >( Id().eref().data() );
    Id kin = makeReacTest();
    double volume = 1e-21;
    Field< double >::set( kin, "volume", volume );
    Field< double >::set( ObjId( "/kinetics/A" ), "concInit", 2 );
    Field< double >::set( ObjId( "/kinetics/e1Pool" ), "concInit", 1 );
    Field< double >::set( ObjId( "/kinetics/e2Pool" ), "concInit", 1 );
    Id e1( "/kinetics/e1Pool/e1" );
    Field< double >::set( e1, "Km", 5 );
    Field< double >::set( e1, "kcat", 1 );
    vector< double > stim( 100, 0.0 );
    for ( unsigned int i = 0; i< 100; ++i ) {
        stim[i] = volume * NA * (1.0 + sin( i * 2.0 * PI / 100.0 ) );
    }
    Field< vector< double > >::set( ObjId( "/kinetics/tab" ), "vector", stim );


    Id gsolve = s->doCreate( "Gsolve", kin, "gsolve", 1 );
    Id stoich = s->doCreate( "Stoich", gsolve, "stoich", 1 );
    Field< Id >::set( stoich, "compartment", kin );
    Field< Id >::set( stoich, "ksolve", gsolve );

    Field< string >::set( stoich, "path", "/kinetics/##" );
    s->doUseClock( "/kinetics/gsolve", "process", 4 );
    s->doSetClock( 4, simDt );

    s->doReinit();
    s->doStart( 20.0 );
    Id plots( "/kinetics/plots" );
    for ( unsigned int i = 0; i < 7; ++i ) {
        stringstream ss;
        ss << "plot." << i;
        SetGet2< string, string >::set( ObjId( plots, i ), "xplot",
                        "tsr3.plot", ss.str() );
    }
    s->doDelete( kin );
    cout << "." << flush;
}

Here is the call graph for this function:

Here is the caller graph for this function:

void testRunKsolve

(

)

Definition at line 235 of file testKsolve.cpp.

References Eref::data(), Shell::doCreate(), Shell::doDelete(), Shell::doReinit(), Shell::doSetClock(), Shell::doStart(), Shell::doUseClock(), Id::eref(), makeReacTest(), Field< A >::set(), and SetGet2< A1, A2 >::set().

Referenced by testKsolve().

{
    double simDt = 0.1;
    // double plotDt = 0.1;
    Shell* s = reinterpret_cast< Shell* >( Id().eref().data() );
    Id kin = makeReacTest();
    Id ksolve = s->doCreate( "Ksolve", kin, "ksolve", 1 );
    Id stoich = s->doCreate( "Stoich", ksolve, "stoich", 1 );
    Field< Id >::set( stoich, "compartment", kin );
    Field< Id >::set( stoich, "ksolve", ksolve );
    Field< string >::set( stoich, "path", "/kinetics/##" );
    s->doUseClock( "/kinetics/ksolve", "process", 4 );
    s->doSetClock( 4, simDt );

    s->doReinit();
    s->doStart( 20.0 );
    Id plots( "/kinetics/plots" );
    for ( unsigned int i = 0; i < 7; ++i ) {
        stringstream ss;
        ss << "plot." << i;
        SetGet2< string, string >::set( ObjId( plots, i ), "xplot",
                        "tsr2.plot", ss.str() );
    }
    s->doDelete( kin );
    cout << "." << flush;
}

Here is the call graph for this function:

Here is the caller graph for this function:

void testSetupReac

(

)

Definition at line 137 of file testKsolve.cpp.

References Eref::data(), Shell::doDelete(), Shell::doReinit(), Shell::doStart(), Id::eref(), and makeReacTest().

Referenced by testKsolve().

{
    Shell* s = reinterpret_cast< Shell* >( Id().eref().data() );
    Id kin = makeReacTest();
    s->doReinit();
    s->doStart( 20.0 );
    Id plots( "/kinetics/plots" );
    /*
    for ( unsigned int i = 0; i < 7; ++i ) {
        stringstream ss;
        ss << "plot." << i;
        SetGet2< string, string >::set( ObjId( plots, i ), "xplot",
                        "tsr.plot", ss.str() );
    }
    */
    s->doDelete( kin );
    cout << "." << flush;
}

Here is the call graph for this function:

Here is the caller graph for this function: