Time marching loops

Time Marching Loops is part of a free web series, ChemPlugin Modeling with Python, by Aqueous Solutions LLC.

Controlling instances

Creating a ChemPlugin instance is a simple matter of declaring it. For example, the statement cp = ChemPlugin() creates a reference “cp” to a ChemPlugin instance. When “cp” is created, it spawns the instance to which it refers.

ChemPlugin instances run silently by default, but can be set to write console messages that trace the instance's actions and report any errors encountered. Some options include:

    cp = ChemPlugin()                            # Mute instance
    cp = ChemPlugin("stdout")                    # Messages to stdout
    cp = ChemPlugin("stderr")                    # Messages to stderr
    cp = ChemPlugin("MyMessages.txt")            # Messages to file
    cp = ChemPlugin("", "-d my_thermo.tdat")     # Command line options
                

Member functions allow a client to control ChemPlugin instances. Given a reference “cp” to a ChemPlugin instance, the statement cp.Config("pH = 5") passes a command telling instance “cp” to set initial pH to 5. Similarly, the statement cp.Initialize() solves for the initial state of an instance, given its configuration.

Note the form of a member function call: the instance in question, the name of the member function, and any arguments.

The “Console()” member function controls where an instance's console output is directed. Console output consists of routine messages an instance produces as it initializes and undertakes calculations, as well as any warning and error messages that may be generated. Depending upon your application and the number of instances you are using, you may wish to redirect the console output:

    cp.Console("stdout")                         # Send to stdout
    cp.Console("MyMessages.txt")                 # Send to file
    cp.Console()                                 # Stop messages
                

The “Config()” member function passes configuration commands to a ChemPlugin instance.

    err = cp.Config("cmd")
    err = cp.Config("cmd1; cmd2; cmd3")
                

    cmds = ("cmd1", "cmd2", "cmd3")
    for cmd in cmds:
        cp.Config(cmd)
                

The configuration commands, which are similar to those used for React, can be found in the “Configuration Commands” section of the ChemPlugin User's Guide.

The “Initialize()” member function triggers an instance to calculate its initial condition, including the distribution of mass across the various chemical species, and to prepare the instance to begin time marching.

    err = cp.Initialize()
    err = cp.Initialize(end_time, ± unit)
                

“Initialize()” is commonly called without arguments, but the client can use the call to specify the end time of the simulation (the default end time is 1 day).

Time Marching

The set of five member functions below work together to make up a time marching loop.

    deltat = cp.ReportTimeStep(± unit)
    err = cp.AdvanceTimeStep(deltat, ± unit)
    err = cp.AdvanceTransport()
    err = cp.AdvanceHeatTransport()
    err = cp.AdvanceChemicalTransport()
                

Member function “ReportTimeStep()” returns the largest time step an instance may take if it is to maintain numerical stability and honor any constraints the client program may have prescribed.

“AdvanceTimeStep()” moves the current time level carried in an instance forward by the time step specified.

“AdvanceTransport()” triggers the instance to evaluate the effect of mass transport on the instance's chemical composition over the course of the current time step.

“AdvanceHeatTransport()” triggers the instance to evaluate the effect of heat transport on the instance's temperature over the course of the current time step.

“AdvanceChemical()” causes the instance to evaluate the effect of kinetic and equilibrium reactions on an instance's chemical state over the course of the current time step.

The simple time marching loop below is composed of an infinite loop and three member functions:

    # Time Marching loop
    while True:
        deltat = cp.ReportTimeStep()
        if cp.AdvanceTimeStep(deltat): break
        if cp.AdvanceChemical(): break
        print(" Results = %4.2f" % cp.Report1(variable, unit))
                

First, a call to “ReportTimeStep()” triggers the instance to calculate an appropriate time step size. The function returns a value for Δt in seconds, and the program stores the value in variable “deltat”.

Next, a call to “AdvanceTimeStep()” passes “deltat” to the instance. Upon executing the function, the instance moves forward in reaction progress by “deltat” seconds. The function returns a value of zero, unless the client has reached the end of the simulation. In that case, a non-zero return breaks the loop and time marching ceases.

Finally, executing member function “AdvanceChemical()” causes the instance to evaluate the chemical equations at the new point in time. The function returns zero, unless an error occurs. A non-zero return marks an error, again breaking the loop.

Once the step is complete, the client writes out a result and returns to take another step.

In another example of a time marching loop

    # Time Marching loop
    while True:
        deltat = cp.ReportTimeStep()
        if cp.AdvanceTimeStep(deltat): break
        if cp.AdvanceChemical(): break
        pH = cp.Report1("pH")
        print(" pH = %4.2f" % pH)
        if pH > 8.0: break
                

a client can be made to break out of the loop upon exceeding a threshold for some variable, like pH.

Task 1: Titration model

In this exercise, we'll create, configure, and initialize a ChemPlugin instance, then go through a time marching loop to simulate a simple titration. The client program to do this is mostly complete, but you'll have to insert a few lines of working code.

To get started, right-click on Titration1_skel.py and open it in Notepad++.

The client program is a simple titration model.

First, find where the instance is configured and replace the “# >>>” line with code to apply the configuration commands we've added to the “cmds” list

    # Configure the instance
    cmds = ("Ca++ = 1 mmol/kg",
            "Na+ = 1 mmol/kg",
            "Cl- = 3 mmol/kg",
            "HCO3- = 2 mmol/kg",
            "pH = 4")
            
    cmds += ("react 3 mmol/kg NaOH",
             "delxi = 0.01")
             
    # >>> use "cmds" to configure "cp" <<< 
                

    for cmd in cmds:
        cp.Config(cmd)
                    

Scroll down and replace the next “# >>>” line with code to initialize the instance:

    # Initialize the instance
    # >>> initialize "cp" <<<
    print(" Titrated %4.2f mmol, pH = %4.2f" % 
           (cp.Report1("mass_reacted", "mmol"), cp.Report1("pH")))
                

    cp.Initialize()
                    

Finally, replace the “# >>>” lines within the time marching loop with code to figure the largest allowable time step, move the current time level forward by that amount, and then evaluate the effects of chemical reactions:

    # Time Marching loop
    while True:
        # >>> Report time step <<< 
        # >>> Advance time step <<< 
        # >>> Advance chemical equations <<< 
        print(" Titrated %4.2f mmol, pH = %4.2f" % 
           (cp.Report1("mass_reacted", "mmol"), cp.Report1("pH")))
                

        deltat = cp.ReportTimeStep()
        if cp.AdvanceTimeStep(deltat): break
        if cp.AdvanceChemical(): break
                    

Save the script.

To see how it works, double-click the modified Python script. The Command Prompt will launch:

Upon initialization, the ChemPlugin instance loads 17 aqueous species and 16 minerals to consider in the calculation, based on our choice of basis entries Ca++, Na+, Cl-, HCO3-, and H+. At that point, no reactant has been added; the instance reports its initial pH of 4.

Subsequent lines of output show how the pH changes as NaOH is added to the system. The extent of a ChemPlugin instance is 1 kg solvent by default, so the configuration commands

    "react 3 mmol/kg NaOH"
    "delxi = 0.01"
                

set up a reaction path in which 3 mmoles of NaOH are added to the system in increments of 0.03 mmoles. If “delxi” was instead set to 0.1, the same mass of NaOH would be added in increments of 0.3 mmoles.

Scroll down to see the point in the simulation where calcite reaches saturation. At that point, it gets swapped into the Basis so that it can precipitate:

Finally, the simulation reaches its endpoint

and the instance reports a successful completion of the reaction path.

Task 2: Another time marching loop

Building upon the previous example, set up the time marching loop so that the titration stops at a specified point. In this case, let's stop titrating NaOH once the pH exceeds 8.

Returning to our script in Notepad++, replace the last two lines in the time marching loop with code to retrieve the pH and break out of the loop once the target pH is exceeded.

    pH = cp.Report1("pH")
    print(" Titrated %4.2f mmol, pH = %4.2f" %
           (cp.Report1("mass_reacted", "mmol"), pH))
    if pH > 8.0: break
                    

Save your script and rerun. Your results should look like this:

How much NaOH was required to reach a pH of 8?

Authors

Craig M. Bethke and Brian Farrell. © Copyright 2016–2024 Aqueous Solutions LLC. This lesson may be reproduced and modified freely to support any licensed use of The Geochemist's Workbench® software, provided that any derived materials acknowledge original authorship.

References

Bethke, C.M., 2022, Geochemical and Biogeochemical Reaction Modeling, 3^rd ed. Cambridge University Press, New York, 520 pp.

Bethke, C.M., 2024, The Geochemist's Workbench®, Release 17: ChemPlugin™ User's Guide. Aqueous Solutions LLC, Champaign, IL, 303 pp.