Homec4science
Diffusion SpecMiCP / ReactMiCP (yaml/safe) Browse

specmicp/cython/specmicp/reaction_path.pyx
f1765f058b94yaml/safe

Last Change
Raw File

reaction_path.pyx
View Options

from libcpp.pair cimport pair
from libcpp.string cimport string
from memory cimport shared_ptr, make_shared
from reaction_path cimport ReactionPathModel, ReactionPathDriver
from database cimport Database, DatabaseManager
from database import DatabaseManager
from solver_options cimport ReducedSystemSolverOptions, penFB_function, min_function
from solver_options cimport MiCPPerformance
from equilibrium_state cimport EquilibriumState, EquilibriumSystem
cdef class ReactionPathSimulation:
"""The simulation to solve"""
cdef shared_ptr[ReactionPathModel] thisptr
def __cinit__(self):
self.thisptr = make_shared[ReactionPathModel]()
cdef ReactionPathModel* _get(self):
return self.thisptr.get()
cdef shared_ptr[ReactionPathModel] get_model(self):
return self.thisptr
def add_aqueous_species(self,
bytes label,
float initial_amount,
float increase_amount):
"""Add an aqueous system in the system"""
self._get().amount_aqueous.insert(
pair[string, pair[float, float]](
label,
pair[float, float](initial_amount, increase_amount)))
def add_mineral(self,
bytes label,
float initial_amount,
float increase_amount):
"""Add a mineral in the system"""
self._get().amount_minerals.insert(
pair[string, pair[float, float]](
label,
pair[float, float](initial_amount, increase_amount)))
def set_list_minerals(self, list_mineral_to_keep):
"""Set the list of solid phases at equilibrium in the system"""
self._get().minerals_to_keep.reserve(len(list_mineral_to_keep))
for mine in list_mineral_to_keep:
self._get().minerals_to_keep.push_back(mine)
cdef class ReactionPathSolver:
"""Solve a reaction path problem"""
cdef ReactionPathDriver* thedriver
cdef VectorXd thesolution
cdef ReactionPathSimulation simulation
cdef DatabaseManager databasemanager
cdef MiCPPerformance current_perf
cdef EquilibriumState current_solution
def __cinit__(self,
ReactionPathSimulation thesimulation,
DatabaseManager databasemanager):
self.simulation = thesimulation
self.databasemanager = databasemanager
self.thedriver = new ReactionPathDriver(thesimulation.get_model(),
databasemanager.get_raw_db(),
self.thesolution
)
def one_step(self):
self.current_perf = self.thedriver.one_step(self.thesolution)
self.current_solution = self.thedriver.get_current_solution()
def get_current_solution(self):
obj = EquilibriumSystem()
obj.initialize(self.current_solution)
return obj
# --------------- #
# #
# Performance #
# #
# --------------- #
def return_code(self):
return <int>self.current_perf.return_code
def number_iterations(self):
return self.current_perf.nb_iterations
def number_factorizations(self):
return self.current_perf.nb_factorization
# --------------- #
# #
# options #
# #
# --------------- #
cdef ReducedSystemSolverOptions* _get(self):
return self.thedriver.get_ptr_options()
def use_penalized_fischer_burmeister(self, float penalization_factor=0.8):
"""Use the penalized Fischer-Burmeister function"""
self._get().ncp_function = penFB_function
self._get().solver_options.penalization_factor = penalization_factor
def use_fischer_burmeister(self):
"""Use the Fischer-Burmeister function"""
self._get().ncp_function = penFB_function
self._get().solver_options.penalization_factor = 1.0
def use_min_function(self):
"""Use the minimum function"""
self._get().ncp_function = min_function
def get_ncp_function(self):
"""Return the NCP function used by the solver.
Return a tuple (ncp_function, penalization factor)"""
if self._get().ncp_function == penFB_function:
return ("Penalized Fischer-Burmeister",
self._get().solver_options.penalization_factor )
else:
return("Minimum function", 0.0)
property allow_restart:
"""Allow the solver to restart the problem in case of failure"""
def __get__(self):
return self._get().allow_restart
def __set__(self, value):
self._get().allow_restart = value
property charge_keeper:
"""Explicit charge conservation, replace mass conservation equation of
the aqueous component 'charge keeper' by the charge neutrality equation
Expect the id of the component
#TODO : accept the label
"""
def __get__(self):
return self._get().charge_keeper
def __set__(self, value):
self._get().charge_keeper = value
property tolerance:
"""Tolerance of the solver"""
def __get__(self):
return self._get().solver_options.fvectol
def __set__(self, value):
self._get().solver_options.fvectol = value
property step_tolerance:
"""Tolerance on the minimum step length"""
def __get__(self):
return self._get().solver_options.steptol
def __set__(self, value):
self._get().solver_options.steptol = value
property max_iterations:
"""The maximum number of iterations"""
def __get__(self):
return self._get().solver_options.max_iter
def __set__(self, value):
self._get().solver_options.max_iter = value
property max_iterations_at_maxstep:
"""The maximum number of iterations at maximum step length"""
def __get__(self):
return self._get().solver_options.maxiter_maxstep
def __set__(self, value):
if (value < self.max_iterations):
raise ValueError("The maximum number of iterations at"
" maximum step length cannot be less"
" than the number of maximum iterations")
self._get().solver_options.maxiter_maxstep == value
property maximum_step_length:
"""The maximum step length"""
def __get__(self):
return self._get().solver_options.maxstep
def __set__(self, value):
self._get().solver_options.maxstep = value
property maximum_factorization_steps:
"""Maximum factorization steps in one iterations"""
def __get__(self):
return self._get().solver_options.max_factorization_step
def __set__(self, value):
self._get().solver_options.max_factorization_step = value
property condition_number_threshold:
"""The maximum condition number allowed for the jacobian."""
def __get__(self):
return self._get().solver_options.condition_limit
def __set__(self, value):
self._get().solver_options.condition_limit = value
property use_crashing:
"""Use crashing iterations to improve the starting guess."""
def __get__(self):
return self._get().solver_options.use_crashing
def __set__(self, value):
self._get().solver_options.use_crashing = value
property use_jacobian_scaling:
"""Scale the jacobian before solving the linear system."""
def __get__(self):
return self._get().solver_options.use_scaling
def __set__(self, value):
self._get().solver_options.use_scaling = value
property factor_gradient_search_direction:
"""The factor by which the steepest deschent direction is multiplied
when trying to improve the Newton direction"""
def __get__(self):
return self._get().solver_options.factor_gradient_search_direction
def __set__(self, value):
self._get().solver_options.factor_gradient_search_direction = value
property threshold_projection:
"""Threshold used to decide if a complementarity variable is zero"""
def __get__(self):
return self._get().solver_options.projection_min_variable
def __set__(self, value):
self._get().solver_options.projection_min_variable = value
property threshold_stationary_point:
"""Threshold to decided if a point is a stationary point or not"""
def __get__(self):
return self._get().solver_options.threshold_stationary_point
def __set__(self, value):
self._get().solver_options.threshold_stationary_point = value
property factor_threshold_linesearch:
"""factor to determine if linesearch is necessary or not"""
def __get__(self):
return self._get().solver_options.coeff_accept_newton_step
def __set__(self, value):
self._get().solver_options.coeff_accept_newton_step = value
c4science · Help