Simulation tool developer: Solver

The code generation functionality in libCellML can be used to transform CellML models into procedural code for simulation. This example works through that process.

C++ resources: All C++ resources

• CMakeLists.txt

• example_solveGeneratedModel.cpp

• resources/

  • sineComparisonExample.c

  • sineComparisonExample.h

Python resources: All Python resources

• example_solveGeneratedModel.py

• resources/

  • sineComparisonExample.py

Create a placeholder for the solver

Interacting with procedural code means that we need a program, so the first step is to create what will soon become our solver. To start with, the program only reads and interprets command line arguments; the size and number of timesteps to simulate. This is done using a simple function, as below.

C++Python

std::map<std::string, std::string> processArguments(int argc, char **argv)
{
    if (argc == 1) {
        std::cout << "Usage:" << std::endl;
        std::cout << "  -n maxSteps -dt stepSize" << std::endl;
        std::cout << "     -n     the number of steps to take before stopping" << std::endl;
        std::cout << "     -dt    the step size to use" << std::endl;
        exit(1);
    }

    std::map<std::string, std::string> argMap;
    std::string value = argv[0];
    value.erase(0, 8); // removes the "solve_" from the start of the executable name to get back to the input
    argMap["input"] = value;

    for (size_t i = 0; i < argc - 1; ++i) {
        if (argv[i][0] == '-') {
            std::string key = argv[i];
            key.erase(0, 1);
            value = argv[i + 1];
            argMap[key] = value;
            i++;
        }
    }
    return argMap;
}

Once the function exists, we can call it from the main function to process the command line arguments.

C++Python

    // Retrieve the command line arguments.
    std::map<std::string, std::string> args = processArguments(argc, argv);
    double stepSize = std::stod(args["dt"]);
    int stepCount = std::stoi(args["n"]);
    std::string input = args["input"];

    std::cout << "-------------------------------------------------------------" << std::endl;
    std::cout << "   SIMPLE SOLVER  " << std::endl;
    std::cout << "-------------------------------------------------------------" << std::endl;
    std::cout << "       model = " << input << std::endl;
    std::cout << "       timestep = " << stepSize << std::endl;
    std::cout << "       number of steps = " << stepCount << std::endl
              << std::endl;

Connect with the generated code

The profile language will affect how you need to interact with the generated code.

For Python, the generated code must be converted into a module, and then imported. This is shown in the example code below.

For C it’s a little more complicated, as you will need to compile the solver whilst including the generated interface *.h file and linking to the implementation *.c file. A CMake file is provided which will do this for you. The build process needs to know the name of the implementation file containing the generated code, so this is given as the -DINPUT=yourBaseFileName argument (without extension). This file and the corresponding interface .h file are copied by CMake to two files named modelToSolve.cpp and modelToSolve.h respectively. This step is required so that within the solver code we can #include a known file name.

C building and linkingPython

Navigate into the folder containing the generated code sineComparisonExample.[c,h] as well as the :code:` example_solveGeneratedModel.cpp` source code from above.

cmake -DINPUT=sineComparisonExample .

You should see an output similar to this:

-- The C compiler identification is AppleClang 10.0.1.10010046
-- The CXX compiler identification is AppleClang 10.0.1.10010046
-- Check for working C compiler: /Library/Developer/CommandLineTools/usr/bin/cc
-- Check for working C compiler: /Library/Developer/CommandLineTools/usr/bin/cc -- works
-- Detecting C compiler ABI info
-- Detecting C compiler ABI info - done
-- Detecting C compile features
-- Detecting C compile features - done
-- Check for working CXX compiler: /Library/Developer/CommandLineTools/usr/bin/c++
-- Check for working CXX compiler: /Library/Developer/CommandLineTools/usr/bin/c++ -- works
-- Detecting CXX compiler ABI info
-- Detecting CXX compiler ABI info - done
-- Detecting CXX compile features
-- Detecting CXX compile features - done

1) First use 'make -j' to build the executable
2) Then solve by running: ./solve_sineComparisonExample with the arguments:
-n  step_total
-dt step_size

-- Configuring done
-- Generating done
-- Build files have been written to: your/file/location/here

Note that the combined program is now available with the prefix solve_ before the base file name you provided with the -DINPUT argument, and can be run using the instructions given in the printout above.

Retrieve the generated model information

The generated code contains model information in dictionaries, as well as functions to interface with the model’s mathematics. These are:

• VOI_INFO: a dictionary with the .name, .units, and .component attributes related to the variable of integration,

• STATE_INFO: a list of similar information for the state variables,

• VARIABLE_INFO: a list of similar information for the non-state variables.

C++Python

    // Inside the generated code are structures with information about the 
    // model and its dimensions.  These are:
    //   - VOI_INFO: a dict with the name, units, and component of the variable of integration,
    //   - STATE_INFO: a list of dicts for the state variables,
    //   - VARIABLE_INFO: a list of dicts for the non-state variables. 

    std::cout << "   VARIABLE OF INTEGRATION (units) " << std::endl;
    std::cout << "      " << VOI_INFO.name << " (" << VOI_INFO.units << ")" <<std::endl 
              << std::endl;
    std::cout << "   STATE VARIABLES (units) " << std::endl;
    std::cout << "-------------------------------------------------------------" << std::endl;
    for (size_t i = 0; i < STATE_COUNT; ++i) {
        std::cout << "      " << STATE_INFO[i].name << " (" << STATE_INFO[i].units << ")" << std::endl;
    }
    std::cout << std::endl;

Allocate space for the solution

Also within the generated code are functions to allocate space for the variables:

• create states array: to construct arrays for storage of the state variables and their rates;

• create variables array: to construct an array to store the other variables.

C++Python

    // Call module functions to construct and initialise the variable arrays.
    // Note that both the rates and the states arrays have the same dimensions,
    // so it's possible to call the createStatesArray() function for both.
    double time = 0.0;
    auto myStateVariables = createStatesArray();
    auto myRates = createStatesArray();
    auto myVariables = createVariablesArray();

Retrieve the model’s mathematical formulation

The mathematical equations that govern the model’s behaviour can be applied using functions from within the generated code:

• initialise states and constants: does what it says, sets all the initial value attributes;

• compute computed constants: computes the value of any constants which depend on others;

• compute variables: calculates those variables whose values depend on the state variables; and

• compute rates: calculates the rates of change of the state variables.

Note that all model variables which affect the rates’ values (and thereby affect the states’ values) are updated in the compute rates function. This means that you only need to call the compute variables function when you’re saving the output from a step; it does not need to be called for intermediate timesteps.

Before we begin iterating, the values of all variables are calculated.

C++Python

    // Make use of the access functions provided to initialise the variable arrays.
    initialiseStatesAndConstants(myStateVariables, myVariables);
    computeComputedConstants(myVariables);
    computeRates(time, myStateVariables, myRates, myVariables);
    computeVariables(time, myStateVariables, myRates, myVariables);

Prepare a file for the solution output

If you have an alternative way to save your solution, you can skip this step. Here we create a simple text-delimited file into which the solution can be written at each timestep.

C++Python

    // Prepare a file for writing during the solution process.
    std::cout << "   INITIAL CONDITIONS" << std::endl;
    std::cout << "-------------------------------------------------------------" << std::endl;
    for (size_t i = 0; i < STATE_COUNT; ++i) {
        std::cout << "      " << STATE_INFO[i].name << "(" << VOI_INFO.name << " = 0) = " << myStateVariables[i] << std::endl;
    }
    for (size_t i = 0; i < VARIABLE_COUNT; ++i) {
        std::cout << "      " << VARIABLE_INFO[i].name << "(" << VOI_INFO.name << " = 0) = " << myVariables[i] << std::endl;
    }

    std::cout << std::endl;
    std::string outFileName = args["input"] + "_solution.txt";
    std::ofstream outFile(outFileName);

    // Header line for output file
    outFile << "iteration";
    outFile << "\t" << VOI_INFO.name << " (" << VOI_INFO.units << ")";
    for (size_t s = 0; s < STATE_COUNT; ++s) {
        outFile << "\t" << STATE_INFO[s].name;
    }
    for (size_t s = 0; s < VARIABLE_COUNT; ++s) {
        outFile << "\t" << VARIABLE_INFO[s].name;
    }
    outFile << std::endl;

    // Initial conditions in output file
    outFile << 0;
    outFile << "\t" << 0;
    for (size_t s = 0; s < STATE_COUNT; ++s) {
        outFile << "\t" << myStateVariables[s];
    }
    for (size_t s = 0; s < VARIABLE_COUNT; ++s) {
        outFile << "\t" << myVariables[s];
    }
    outFile << std::endl;

Perform the integration steps

Finally we iterate through the timesteps, calculating the state variables, and updating the rates each step. The solution values and calculated variables are written to the output file.

C++Python

    // Numerically integrate the state variables using the Euler method to step through the solution.

    // Solution columns in output file
    for (size_t step = 1; step < stepCount; ++step) {
        time = step * stepSize;
        computeRates(time, myStateVariables, myRates, myVariables);
        outFile << step << "\t " << time;
        for (size_t s = 0; s < STATE_COUNT; ++s) {
            myStateVariables[s] = myStateVariables[s] + myRates[s] * stepSize;
            outFile << "\t" << myStateVariables[s];
        }
        // The variables in the "myVariables" array are those which do not affect the calculation
        // of rates or state variables.  They only need to be computed when outputting the 
        // results of a timestep: if you're not saving every timestep, then you can skip this
        // until you are.
        computeVariables(time, myStateVariables, myRates, myVariables); 
        for (size_t s = 0; s < VARIABLE_COUNT; ++s) {
            outFile << "\t" << myVariables[s];
        }
        outFile << "\n";
    }
    outFile.close();

    // Housekeeping.

    deleteArray(myStateVariables);
    deleteArray(myVariables);
    deleteArray(myRates);

The solution files are written in a tab-delimited format which can be read by your favourite plotting application. The plots below were generated using a step size of 0.1 for 100 iterations.

Figure 27 Plots generated from an Euler solution to the sine comparison model for a step size of 0.1.