variables_interface.cpp
No OneTemporary
Actions

Subscribers

None

File Metadata

Created: Sun, Jul 28, 19:22

variables_interface.cpp
View Options

	/*-------------------------------------------------------------------------------

	Copyright (c) 2015 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 "../../solver/staggers_base/variables_base.hpp"

	#include "../../../specmicp/adimensional/adimensional_system_solution_extractor.hpp"

	#include "../../../utils/compat.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 {

	// 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 = {}
	);


	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();
	}
	};


	// =================================
	//
	// 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::VariablesInterfaceImpl::VariablesInterfaceImpl(
	mesh::Mesh1DPtr the_mesh,
	database::RawDatabasePtr the_database,
	const std::vector<index_t>& component_with_gas,
	units::UnitsSet units_set
	):
	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::~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()));
	}

	} //end namespace unsaturated
	} //end namespace systems
	} //end namespace reactmicp
	} //end namespace specmicp

variables_interface.cppNo OneTemporaryActions

File Metadata

variables_interface.cppView Options

Event Timeline

variables_interface.cpp
No OneTemporary
Actions

variables_interface.cpp
View Options