Homec4science
Diffusion SpecMiCP / ReactMiCP (reactmicp/unsaturated) Browse

specmicp/cython/specmicp/constraints.pyx
b174f1ddf352reactmicp/unsaturated

Last Change
Raw File

constraints.pyx
View Options

#Copyright (c) 2014,2015 Fabien Georget <fabieng@princeton.edu>, Princeton
#University #All rights reserved.
#
#Redistribution and use in source and binary forms, with or without
#modification, are permitted provided that the following conditions are met:
#
#1. Redistributions of source code must retain the above copyright notice, this
#list of conditions and the following disclaimer.
#
#2. Redistributions in binary form must reproduce the above copyright notice,
#this list of conditions and the following disclaimer in the documentation
#and/or other materials provided with the distribution.
#
#3. Neither the name of the copyright holder nor the names of its contributors
#may be used to endorse or promote products derived from this software without
#specific prior written permission.
#
#THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND
#ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
#WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
#DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE
#FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
#DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR
#SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
#CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY,
#OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
#OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
from database cimport DatabaseManager
from constraints cimport SpecMiCPConstraints, AdimensionalSystemConstraints, SpecMiCPFormulation
from eigen cimport VectorXd, vector_setitem
from libcpp.utility cimport pair
from cython.operator cimport dereference
cdef class SpecMiCPConstraints:
"""A class managing the SpecMiCP constraints"""
def __cinit__(self, DatabaseManager database not None):
self.database = database
self.constraints = new AdimensionalSystemConstraints()
def __dealloc__(self):
del self.constraints
cdef AdimensionalSystemConstraints* _get(self):
return self.constraints
# Total concentrations
# --------------------
def set_total_concentrations(self, dictionary_total_concs):
"""Set the total concentrations
dictionary_total_concs is a dictionary mapping the labels of components
to their total concentrations.
"""
cdef VectorXd total_concs
cdef int index
total_concs.setZero(self.database.nb_component)
for (key, value) in dictionary_total_concs.items():
index = self.database.component_label_to_id(key)
vector_setitem(total_concs, index, value)
self.constraints.set_total_concentrations(total_concs)
# Water conservation
# ------------------
def enable_conservation_water(self):
"""Enable the conservation of water"""
self.constraints.enable_conservation_water()
def disable_conservation_water(self):
"""Disable the conservation of water.
No equation is solved for the water component in this case."""
self.constraints.disable_conservation_water()
def set_saturated_system(self):
"""The system is saturated with water"""
self.constraints.set_saturated_system()
# Surface model
# -------------
def enable_surface_model(self, site_concentrations):
"""Enable the surface sorption model.
The total concentration of sorption site is 'site_concentrations'"""
self.constraints.enable_surface_model(site_concentrations)
def disable_surface_model(self):
"""Disable the surface sorption model"""
self.constraints.disable_surface_model()
# components
# ----------
def set_charge_keeper(self, label):
"""Set the charge keeper to be 'label'
The charge keeper is in charge of the electroneutrality equation"""
cdef int index
index = self.database.component_label_to_id(label)
self.constraints.set_charge_keeper(index)
def add_fixed_fugacity_gas(self, gas_label, component_label, log_fugacity):
"""Set 'gas_label' to be a fixed fugacity gas"""
cdef int g_index, c_index
g_index = self.database.gas_label_to_id(gas_label)
c_index = self.database.component_label_to_id(component_label)
self.constraints.add_fixed_fugacity_gas(g_index, c_index, log_fugacity)
def add_fixed_activity_component(self, component_label, log_activity):
"""Set 'component_label' to have a fixed fugacity"""
cdef int c_index
c_index = self.database.component_label_to_id(component_label)
self.constraints.add_fixed_activity_component(c_index, log_activity)
# Inert volume
# ------------
def set_inert_volume_fraction(self, volume_fraction_inert):
"""Set the volume fraction of inert materials"""
self.constraints.set_inert_volume_fraction(volume_fraction_inert)
cdef class SpecMiCPFormulation(SpecMiCPConstraints):
"""The formulation of a problem"""
def __cinit__(self, DatabaseManager database):
self.formulation = new Formulation()
def __dealloc__(self):
del self.formulation
def set_mass_solution(self, mass_solution):
"""Set the mass of the solution"""
self.formulation.set_mass_solution(mass_solution)
def add_aqueous_species(self, label, concentration):
"""Add an aqueous species to the system.
An aqueous species can be in that case a component or a secondary
aqueous species.
'concentration' is in moles per volume
"""
self.formulation.add_aqueous_species(label, concentration)
def add_mineral(self, label, concentration):
"""Add a mineral species to the system.
The mineral can be governed by equilibrium or by kinetics.
'concentration' is in moles per volume
"""
self.formulation.add_mineral(label, concentration)
def keep_component(self, label):
"""Keep the component 'label' in the database
even if it does not exist in the initial system.
This must be used if the component will be used later in the computation
"""
self.formulation.keep_component(label)
def set_list_minerals(self, list_minerals):
"""Set the list of minerals at equilibrium
If this function is not called, or the size of the list is zero,
all minerals in the database will be used.
"""
cdef vector[string] list = list_minerals
self.formulation.set_minerals_list(list)
def initialize_system(self):
"""Initialize the system by dissolving the formulation"""
cdef Dissolver* dissolver = new Dissolver(self.database.get_raw_db())
cdef VectorXd total_conc = dissolver.dissolve(dereference(self.formulation))
self._get().set_total_concentrations(total_conc)
del dissolver
c4science · Help