from ChemPlugin import *

print("Model a flow-through reactor\n")

# Configure and initialize the inlet fluid.
cp_inlet = ChemPlugin("stdout")
cmds = "Ca++ = 1 mmol/kg; HCO3- = 1 mmol/kg; pH = 1; balance on Cl-"
cp_inlet.Config(cmds)
cp_inlet.Initialize()

# Configure and initialize the stirred reactor.
cp = ChemPlugin("stdout")
cmds = ("swap Calcite for Ca++", "Calcite = 0.03 free m3", "Cl- = 1 mmol/kg",
        "swap CO2(g) for H+", "fugacity CO2(g) = 1", "balance on HCO3-",
        "volume = 1 m3", "fix f CO2(g)", "delxi = 0.01", "pluses = banner")

for c in cmds: cp.Config(c)
cp.Initialize(1.0, "day")

# Link reactor to inlet and free outlet; set rate of flow.
link1 = cp.Link(cp_inlet)
link2 = cp.Link()

link1.FlowRate(10.0, "m3/day")
link2.FlowRate(-10.0, "m3/day")

cp.PlotHeader("FlowThrough.gtp", "char")
cp.PlotBlock()

# Time marching loop.
while True:
    deltat = cp.ReportTimeStep()
    if cp.AdvanceTimeStep(deltat): break
    if cp.AdvanceTransport(): break
    if cp.AdvanceChemical(): break
    cp.PlotBlock()

input()