Page MenuHomec4science

variables_interface.cpp
No OneTemporary

File Metadata

Created
Fri, Aug 23, 05:59

variables_interface.cpp

/* =============================================================================
Copyright (c) 2014 - 2016
F. 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. *
============================================================================= */
#include "variables_interface.hpp"
#include "variables.hpp"
#include "variables_box.hpp"
#include "reactmicp/solver/staggers_base/variables_base.hpp"
#include "specmicp/adimensional/adimensional_system_solution_extractor.hpp"
#include "specmicp_common/compat.hpp"
#include "specmicp_common/log.hpp"
#include "specmicp_common/eigen/value_checker.hpp"
#include "specmicp_common/eigen/vector_checker.hpp"
#include <stdexcept>
// This is a boring class, with only small differences between
// all the functions and variables. But, they make all the
// difference at the end. The main purpose of the class is to
// hide all the variables structure and only present a unified
// interface to the user.
// None of the structure defined in "variables_sub.hpp" should
// be leaked outside this class
namespace specmicp {
namespace reactmicp {
namespace systems {
namespace unsaturated {
using namespace specmicp::vector_checker;
// typedef
// -------
using VariablesBase = solver::VariablesBase;
using VariablesBasePtr = std::shared_ptr<VariablesBase>;
// =================================
//
// Declaration of the implementation
//
// =================================
// This class just store and retrieve variables group
struct VariablesInterface::VariablesInterfaceImpl
{
UnsaturatedVariablesPtr m_vars;
mesh::Mesh1DPtr m_mesh;
database::RawDatabasePtr m_database;
VariablesInterfaceImpl(
mesh::Mesh1DPtr the_mesh,
database::RawDatabasePtr the_database,
const std::vector<index_t>& component_with_gas,
units::UnitsSet units_set = {}
);
VariablesInterfaceImpl(
std::shared_ptr<UnsaturatedVariables> vars
);
MainVariable& liquid_saturation() {
return m_vars->get_liquid_saturation();
}
const MainVariable& liquid_saturation() const {
return m_vars->get_liquid_saturation();
}
MainVariable& aqueous_concentration(index_t component) {
return m_vars->get_aqueous_concentration(component);
}
const MainVariable& aqueous_concentration(index_t component) const {
return m_vars->get_aqueous_concentration(component);
}
MainVariable& partial_pressure(index_t component) {
if (not m_vars->component_has_gas(component)) {
throw std::logic_error("Component has no associated gas");
}
return m_vars->get_pressure_main_variables(component);
}
const MainVariable& partial_pressure(index_t component) const {
return m_vars->get_pressure_main_variables(component);
}
MainVariable& solid_concentration(index_t component) {
return m_vars->get_solid_concentration(component);
}
const MainVariable& solid_concentration(index_t component) const {
return m_vars->get_solid_concentration(component);
}
SecondaryTransientVariable& porosity() {
return m_vars->get_porosity();
}
const SecondaryTransientVariable& porosity() const {
return m_vars->get_porosity();
}
SecondaryTransientVariable& water_aqueous_concentration() {
return m_vars->get_water_aqueous_concentration();
}
const SecondaryTransientVariable& water_aqueous_concentration() const {
return m_vars->get_water_aqueous_concentration();
}
SecondaryVariable& liquid_diffusivity() {
return m_vars->get_liquid_diffusivity();
}
const SecondaryVariable& liquid_diffusivity() const {
return m_vars->get_liquid_diffusivity();
}
SecondaryVariable& liquid_permeability() {
return m_vars->get_liquid_permeability();
}
const SecondaryVariable& liquid_permeability() const {
return m_vars->get_liquid_permeability();
}
scalar_t& binary_gas_diffusivity(index_t component) {
return m_vars->get_binary_gas_diffusivity(component);
}
scalar_t binary_gas_diffusivity(index_t component) const {
return m_vars->get_binary_gas_diffusivity(component);
}
SecondaryVariable& resistance_gas_diffusivity() {
return m_vars->get_resistance_gas_diffusivity();
}
const SecondaryVariable& resistance_gas_diffusivity() const {
return m_vars->get_resistance_gas_diffusivity();
}
SecondaryVariable& advection_flux() {
return m_vars->get_advection_flux();
}
const SecondaryVariable& advection_flux() const {
return m_vars->get_advection_flux();
}
VariablesValidity check_variables();
bool check_bounds(const Vector& to_check, scalar_t lower, scalar_t upper);
VariablesValidity check_saturation();
VariablesValidity check_relative();
VariablesValidity check_porosity();
VariablesValidity check_transport_param();
};
// =================================
//
// Implementation
//
// =================================
VariablesInterface::VariablesInterface(
mesh::Mesh1DPtr the_mesh,
database::RawDatabasePtr the_database,
const std::vector<index_t>& component_with_gas):
m_impl(make_unique<VariablesInterfaceImpl>(
the_mesh, the_database, component_with_gas))
{}
VariablesInterface::VariablesInterface(
mesh::Mesh1DPtr the_mesh,
database::RawDatabasePtr the_database,
const std::vector<index_t>& component_with_gas,
const units::UnitsSet& units_set
):
m_impl(make_unique<VariablesInterfaceImpl>(
the_mesh, the_database, component_with_gas, units_set))
{}
VariablesInterface::VariablesInterface(
std::shared_ptr<UnsaturatedVariables> vars
):
m_impl(make_unique<VariablesInterfaceImpl>(vars))
{}
VariablesInterface::VariablesInterfaceImpl::VariablesInterfaceImpl(
mesh::Mesh1DPtr the_mesh,
database::RawDatabasePtr the_database,
const std::vector<index_t>& component_with_gas,
units::UnitsSet units_set
):
m_vars(nullptr),
m_mesh(the_mesh),
m_database(the_database)
{
std::vector<bool> comp_has_gas(the_database->nb_component(), false);
for (auto component: component_with_gas) {
comp_has_gas[component] = true;
}
m_vars = std::make_shared<UnsaturatedVariables>(
the_mesh, the_database, comp_has_gas, units_set.length);
}
VariablesInterface::VariablesInterfaceImpl::VariablesInterfaceImpl(
std::shared_ptr<UnsaturatedVariables> vars
):
m_vars(vars),
m_mesh(m_vars->get_mesh()),
m_database(m_vars->get_database())
{}
VariablesInterface::~VariablesInterface() = default;
std::shared_ptr<VariablesBase> VariablesInterface::get_base_variables()
{
return std::static_pointer_cast<VariablesBase>(m_impl->m_vars);
}
std::shared_ptr<UnsaturatedVariables> VariablesInterface::get_variables()
{
return m_impl->m_vars;
}
UnsaturatedVariables* VariablesInterface::get_raw_variables()
{
return m_impl->m_vars.get();
}
// Saturation
const Vector& VariablesInterface::get_liquid_saturation() const
{
return m_impl->liquid_saturation().variable;
}
void VariablesInterface::set_liquid_saturation(index_t node, scalar_t value)
{
m_impl->liquid_saturation().variable(node) = value;
m_impl->liquid_saturation().predictor(node) = value;
}
void VariablesInterface::set_liquid_saturation(const Vector& value)
{
m_impl->liquid_saturation().variable = value;
m_impl->liquid_saturation().predictor = value;
}
void VariablesInterface::set_liquid_saturation(scalar_t value)
{
m_impl->liquid_saturation().variable.setConstant(value);
m_impl->liquid_saturation().predictor.setConstant(value);
}
// Aqueous component concentration
const Vector& VariablesInterface::get_aqueous_concentration(index_t component) const
{
return m_impl->aqueous_concentration(component).variable;
}
void VariablesInterface::set_aqueous_concentration(index_t component, index_t node, scalar_t value)
{
MainVariable& aq_conc = m_impl->aqueous_concentration(component);
aq_conc.variable(node) = value;
aq_conc.predictor(node) = value;
}
void VariablesInterface::set_aqueous_concentration(index_t component, const Vector& value)
{
MainVariable& aq_conc = m_impl->aqueous_concentration(component);
aq_conc.variable = value;
aq_conc.predictor = value;
}
void VariablesInterface::set_aqueous_concentration(index_t component, scalar_t value)
{
MainVariable& aq_conc = m_impl->aqueous_concentration(component);
aq_conc.variable.setConstant(value);
aq_conc.predictor.setConstant(value);
}
// Partial pressure
const Vector& VariablesInterface::get_partial_pressure(index_t component) const
{
return m_impl->partial_pressure(component).variable;
}
void VariablesInterface::set_partial_pressure(index_t component, index_t node, scalar_t value)
{
MainVariable& comp_pressure = m_impl->partial_pressure(component);
comp_pressure.variable(node) = value;
comp_pressure.predictor(node) = value;
}
void VariablesInterface::set_partial_pressure(index_t component, const Vector& value)
{
MainVariable& comp_pressure = m_impl->partial_pressure(component);
comp_pressure.variable = value;
comp_pressure.predictor = value;
}
void VariablesInterface::set_partial_pressure(index_t component, scalar_t value)
{
MainVariable& comp_pressure = m_impl->partial_pressure(component);
comp_pressure.variable.setConstant(value);
comp_pressure.predictor.setConstant(value);
}
// Solid concentration
const Vector& VariablesInterface::get_solid_concentration(index_t component) const
{
return m_impl->solid_concentration(component).variable;
}
void VariablesInterface::set_solid_concentration(index_t component, index_t node, scalar_t value)
{
MainVariable& solid_conc = m_impl->solid_concentration(component);
solid_conc.variable(node) = value;
solid_conc.predictor(node) = value;
}
void VariablesInterface::set_solid_concentration(index_t component, const Vector& value)
{
MainVariable& solid_conc = m_impl->solid_concentration(component);
solid_conc.variable = value;
solid_conc.predictor = value;
}
void VariablesInterface::set_solid_concentration(index_t component, scalar_t value)
{
MainVariable& solid_conc = m_impl->solid_concentration(component);
solid_conc.variable.setConstant(value);
solid_conc.predictor.setConstant(value);
}
// Porosity
const Vector& VariablesInterface::get_porosity() const
{
return m_impl->porosity().variable;
}
void VariablesInterface::set_porosity(index_t node, scalar_t value)
{
m_impl->porosity().variable(node) = value;
m_impl->porosity().predictor(node) = value;
}
void VariablesInterface::set_porosity(const Vector& value)
{
m_impl->porosity().variable = value;
m_impl->porosity().predictor = value;
}
void VariablesInterface::set_porosity(scalar_t value)
{
m_impl->porosity().variable.setConstant(value);
m_impl->porosity().predictor.setConstant(value);
}
// Water aqueous concentration
const Vector& VariablesInterface::get_water_aqueous_concentration() const
{
return m_impl->water_aqueous_concentration().variable;
}
void VariablesInterface::set_water_aqueous_concentration(index_t node, scalar_t value)
{
m_impl->water_aqueous_concentration().variable(node) = value;
m_impl->water_aqueous_concentration().predictor(node) = value;
}
void VariablesInterface::set_water_aqueous_concentration(const Vector& value)
{
m_impl->water_aqueous_concentration().variable = value;
m_impl->water_aqueous_concentration().predictor = value;
}
void VariablesInterface::set_water_aqueous_concentration(scalar_t value)
{
m_impl->water_aqueous_concentration().variable.setConstant(value);
m_impl->water_aqueous_concentration().predictor.setConstant(value);
}
// Liquid diffusion coefficient
const Vector& VariablesInterface::get_liquid_diffusivity() const
{
return m_impl->liquid_diffusivity().variable;
}
void VariablesInterface::set_liquid_diffusivity(index_t node, scalar_t value)
{
m_impl->liquid_diffusivity().variable(node) = value;
}
void VariablesInterface::set_liquid_diffusivity(const Vector& value)
{
m_impl->liquid_diffusivity().variable = value;
}
void VariablesInterface::set_liquid_diffusivity(scalar_t value)
{
m_impl->liquid_diffusivity().variable.setConstant(value);
}
// Liquid permeability
const Vector& VariablesInterface::get_liquid_permeability() const
{
return m_impl->liquid_permeability().variable;
}
void VariablesInterface::set_liquid_permeability(index_t node, scalar_t value)
{
m_impl->liquid_permeability().variable(node) = value;
}
void VariablesInterface::set_liquid_permeability(const Vector& value)
{
m_impl->liquid_permeability().variable = value;
}
void VariablesInterface::set_liquid_permeability(scalar_t value)
{
m_impl->liquid_permeability().variable.setConstant(value);
}
// Gas diffusivity
scalar_t
VariablesInterface::get_binary_gas_diffusivity(
index_t component
) const
{
return m_impl->binary_gas_diffusivity(component);
}
void
VariablesInterface::set_binary_gas_diffusivity(
index_t component,
scalar_t value
)
{
m_impl->binary_gas_diffusivity(component) = value;
}
const Vector&
VariablesInterface::get_resistance_gas_diffusivity() const
{
return m_impl->resistance_gas_diffusivity().variable;
}
void
VariablesInterface::set_resistance_gas_diffusivity(
index_t node,
scalar_t value
)
{
m_impl->resistance_gas_diffusivity().variable(node) = value;
}
void
VariablesInterface::set_resistance_gas_diffusivity(
const Vector& value)
{
m_impl->resistance_gas_diffusivity().variable = value;
}
void
VariablesInterface::set_resistance_gas_diffusivity(
scalar_t value)
{
m_impl->resistance_gas_diffusivity().variable.setConstant(value);
}
const Vector& VariablesInterface::get_advection_flux() const
{
return m_impl->advection_flux().variable;
}
void VariablesInterface::set_advection_flux(index_t node, scalar_t value)
{
m_impl->advection_flux().variable(node) = value;
}
void VariablesInterface::set_advection_flux(const Vector& value)
{
m_impl->advection_flux().variable = value;
}
void VariablesInterface::set_advection_flux(scalar_t value)
{
m_impl->advection_flux().variable.setConstant(value);
}
// User Models
const Vector& VariablesInterface::get_capillary_pressure() const
{
return m_impl->m_vars->get_capillary_pressure().variable;
}
void VariablesInterface::set_capillary_pressure_model(
user_model_saturation_f capillary_pressure_model
)
{
m_impl->m_vars->set_capillary_pressure_model(capillary_pressure_model);
const Vector& sat = get_liquid_saturation();
Vector& cap_pressure = m_impl->m_vars->get_capillary_pressure().variable;
for (auto node: m_impl->m_mesh->range_nodes())
{
cap_pressure(node) = capillary_pressure_model(node, sat(node));
}
}
const Vector& VariablesInterface::get_vapor_pressure() const
{
return m_impl->m_vars->get_pressure_main_variables(0).variable;
}
void VariablesInterface::set_vapor_pressure_model(
user_model_saturation_f vapor_pressure_model
)
{
m_impl->m_vars->set_vapor_pressure_model(vapor_pressure_model);
}
const Vector& VariablesInterface::get_relative_liquid_diffusivity() const
{
return m_impl->m_vars->get_relative_liquid_diffusivity().variable;
}
void VariablesInterface::set_relative_liquid_diffusivity_model(
user_model_saturation_f relative_liquid_diffusivity_model
)
{
m_impl->m_vars->set_relative_liquid_diffusivity_model(relative_liquid_diffusivity_model);
const Vector& sat = get_liquid_saturation();
Vector& rel_diff = m_impl->m_vars->get_relative_liquid_diffusivity().variable;
for (auto node: m_impl->m_mesh->range_nodes())
{
rel_diff(node) = relative_liquid_diffusivity_model(node, sat(node));
}
}
const Vector& VariablesInterface::get_relative_liquid_permeability() const
{
return m_impl->m_vars->get_relative_liquid_permeability().variable;
}
void VariablesInterface::set_relative_liquid_permeability_model(
user_model_saturation_f relative_liquid_permeability_model
)
{
m_impl->m_vars->set_relative_liquid_permeability_model(relative_liquid_permeability_model);
const Vector& sat = get_liquid_saturation();
Vector& rel_diff = m_impl->m_vars->get_relative_liquid_permeability().variable;
for (auto node: m_impl->m_mesh->range_nodes())
{
rel_diff(node) = relative_liquid_permeability_model(node, sat(node));
}
}
const Vector& VariablesInterface::get_relative_gas_diffusivity() const
{
return m_impl->m_vars->get_relative_gas_diffusivity().variable;
}
void VariablesInterface::set_relative_gas_diffusivity_model(
user_model_saturation_f relative_gas_diffusivity_model
)
{
m_impl->m_vars->set_relative_gas_diffusivity_model(relative_gas_diffusivity_model);
const Vector& sat = get_liquid_saturation();
Vector& rel_diff = m_impl->m_vars->get_relative_gas_diffusivity().variable;
for (auto node: m_impl->m_mesh->range_nodes())
{
rel_diff(node) = relative_gas_diffusivity_model(node, sat(node));
}
}
// chemistry
void VariablesInterface::initialize_variables(index_t node,
const AdimensionalSystemSolutionExtractor& extractor
)
{
set_porosity(node, extractor.porosity());
scalar_t saturation = extractor.saturation_water();
set_liquid_saturation(node, extractor.saturation_water());
scalar_t rho_l = extractor.density_water();
set_water_aqueous_concentration(node,
rho_l*extractor.total_aqueous_concentration(0));
set_solid_concentration(0, node,
extractor.total_solid_concentration(0));
if (m_impl->m_vars->component_has_gas(0))
{
set_partial_pressure(0, node,
m_impl->m_vars->get_vapor_pressure_model()(node, saturation)
);
}
for (index_t component: m_impl->m_database->range_aqueous_component())
{
set_aqueous_concentration(component, node,
rho_l*extractor.total_aqueous_concentration(component));
set_solid_concentration(component, node,
extractor.total_solid_concentration(component));
if (m_impl->m_vars->component_has_gas(component))
{
set_partial_pressure(component, node,
extractor.fugacity_gas(m_impl->m_vars->get_id_gas(component))
*m_impl->m_vars->get_total_pressure());
}
}
m_impl->m_vars->set_adim_solution(node,
std::move(extractor.get_solution()));
}
VariablesValidity VariablesInterface::check_variables()
{
return m_impl->check_variables();
}
bool VariablesInterface::VariablesInterfaceImpl::check_bounds(
const Vector& to_check,
scalar_t lower,
scalar_t upper
)
{
bool flag = true;
for (index_t ind=0; ind< to_check.rows(); ++ind)
{
const scalar_t& value = to_check(ind);
if (value < lower or value > upper) {
flag = false;
break;
}
}
return flag;
}
VariablesValidity max_error_level(VariablesValidity err1, VariablesValidity err2)
{
return (err1>err2)?err1:err2;
}
VariablesValidity VariablesInterface::VariablesInterfaceImpl::check_variables()
{
VariablesValidity return_code = VariablesValidity::good;
auto tmpret = check_saturation();
return_code = max_error_level(return_code, tmpret);
tmpret = check_porosity();
return_code = max_error_level(return_code, tmpret);
tmpret = check_relative();
return_code = max_error_level(return_code, tmpret);
return return_code;
}
static bool SPECMICP_DLL_LOCAL saturation_check(const Vector& x)
{
return (is_bounded(x, 0.0, 1.0) && (any(x) > 0.0));
}
VariablesValidity
VariablesInterface::VariablesInterfaceImpl::check_saturation()
{
using namespace vector_checker;
const auto& x = m_vars->get_liquid_saturation().variable;
if (not saturation_check(x))
{
ERROR << "Saturation vector is incorrect !";
return VariablesValidity::error;
}
return VariablesValidity::good;
}
static bool SPECMICP_DLL_LOCAL porosity_check(const Vector& x)
{
return (is_bounded(x, 0.0, 1.0) && (any(x) > 0.0));
}
VariablesValidity
VariablesInterface::VariablesInterfaceImpl::check_porosity()
{
using namespace vector_checker;
const auto& x = m_vars->get_porosity().variable;
if (not porosity_check(x)) {
ERROR << "Porosity is incorrect";
return VariablesValidity::error;
}
return VariablesValidity::good;
}
static bool SPECMICP_DLL_LOCAL relative_variable_check(const Vector& x)
{
return (is_bounded(x, 0.0, 1.0) && (any(x) > 0.0));
}
VariablesValidity
VariablesInterface::VariablesInterfaceImpl::check_relative()
{
auto is_valid = VariablesValidity::good;
if (not relative_variable_check(
m_vars->get_relative_liquid_diffusivity().variable)) {
ERROR << "Relative liquid diffusivity is incorrect";
is_valid = VariablesValidity::error;
}
if (not relative_variable_check(
m_vars->get_relative_liquid_permeability().variable)) {
ERROR << "Relative liquid permeability is incorrect";
is_valid = VariablesValidity::error;
}
if (not relative_variable_check(
m_vars->get_relative_gas_diffusivity().variable)) {
ERROR << "Relative gas diffusivity is incorrect";
is_valid = VariablesValidity::error;
}
if (not relative_variable_check(
m_vars->get_resistance_gas_diffusivity().variable)) {
ERROR << "Resistance gas diffusivity is incorrect";
is_valid = VariablesValidity::error;
}
return is_valid;
}
static bool SPECMICP_DLL_LOCAL transport_parameters_check(const Vector& x)
{
return (is_lower_bounded(x, 0.0) && (any(x) > 0.0));
}
VariablesValidity
VariablesInterface::VariablesInterfaceImpl::check_transport_param()
{
auto is_valid = VariablesValidity::good;
if (not transport_parameters_check(
m_vars->get_liquid_diffusivity().variable)) {
ERROR << "Liquid diffusivity is incorrect";
is_valid = VariablesValidity::error;
}
if (not transport_parameters_check(
m_vars->get_liquid_permeability().variable)) {
ERROR << "Liquid permeability is incorrect";
is_valid = VariablesValidity::error;
}
return is_valid;
}
} //end namespace unsaturated
} //end namespace systems
} //end namespace reactmicp
} //end namespace specmicp

Event Timeline