Homec4science
Diffusion SpecMiCP / ReactMiCP (yaml/safe) Browse

specmicp/tests/cython/thermocarbo.py
f401e32a96d2yaml/safe

Last Change
Raw File

thermocarbo.py
View Options

#!/usr/bin/env python
import unittest as ut
import sys
lib_directory="@python_module_path@"
database_path=b"@database_path@"
sys.path.insert(0, lib_directory)
import specmicp.database as database
import specmicp.reaction_path as reaction_path
if __name__=="__main__":
db = database.DatabaseManager(database_path)
swapping = {b"H[+]": b"HO[-]",
b"Si(OH)4": b"SiO(OH)3[-]"
}
db.swap_components(swapping)
model = reaction_path.ReactionPathSimulation()
m_c3s = 0.7
m_c2s = 0.3
wc = 0.5
m_water = wc*((3*56.08+60.08)*m_c3s + (2*56.08+60.08)*m_c2s)*1e-3
delta_h2co3 = 0.1
model.add_aqueous_species(b"H2O", m_water/18.015e-3, delta_h2co3)
model.add_aqueous_species(b"HCO3[-]", 0.0, delta_h2co3)
model.add_aqueous_species(b"HO[-]", 0.0, -delta_h2co3)
model.add_mineral(b"C3S", m_c3s, 0.0)
model.add_mineral(b"C2S", m_c2s, 0.0)
solver = reaction_path.ReactionPathSolver(model, db)
solver.maximum_factorization_steps = 1
solver.maximum_step_length = 50
solver.use_jacobian_scaling = False
#solver.use_min_function()
#solver.use_fischer_burmeister()
totiter = 0
totfact = 0
header = "Reaction\tReturn code\tNb iterations\tNb factorizations"
header += "\tmasswater\tpH"
for i in range(1,db.nb_component):
header += "\t"+db.component_id_to_label(i).decode()
for i in range(db.nb_mineral):
header += "\t"+db.mineral_id_to_label(i).decode()
print(header)
nb_step = 41
for i in range(nb_step):
solver.one_step()
solution = solver.get_current_solution()
output = "{0}\t{1}\t{2}\t{3}\t{4}\t{5}".format(
i*delta_h2co3,
solver.return_code(),
solver.number_iterations(),
solver.number_factorizations(),
solution.mass_water(),
solution.pH())
for i in range(1,db.nb_component):
output += "\t{0}".format(solution.total_aqueous_concentration(i))
for i in range(db.nb_mineral):
output += "\t{0}".format(solution.moles_mineral(i))
print(output)
totiter += solver.number_iterations()
totfact += solver.number_factorizations()
print("Average nb iterations : {0}".format(totiter/nb_step))
print("Average nb factorizations : {0}".format(totfact/nb_step))
c4science · Help