equilibrium_stagger.cpp
No OneTemporary
Actions

Subscribers

None

File Metadata

Created: Tue, May 21, 10:47

equilibrium_stagger.cpp
View Options

	/* =============================================================================

	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 "equilibrium_stagger.hpp"
	#include "variables.hpp"
	#include "../../solver/staggers_base/stagger_structs.hpp"

	#include "../../../dfpm/meshes/mesh1d.hpp"
	#include "../../../database/data_container.hpp"

	#include "../../../utils/compat.hpp"

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

	#include "../../../utils/log.hpp"

	namespace specmicp {
	namespace reactmicp {
	namespace systems {
	namespace unsaturated {

	// ============================================
	//
	// Declaration of implementation details
	//
	// ============================================


	struct EquilibriumStagger::EquilibriumStaggerImpl
	{
	using OptionsVector = EquilibriumStagger::OptionsVector;

	scalar_t m_dt {-1};
	mesh::Mesh1DPtr m_mesh;
	database::RawDatabasePtr m_database;
	std::shared_ptr<BoundaryConditions> m_bcs;
	std::shared_ptr<OptionsVector> m_options;

	EquilibriumStaggerImpl(
	std::shared_ptr<UnsaturatedVariables> variables,
	std::shared_ptr<BoundaryConditions> boundary_conditions,
	std::shared_ptr<OptionsVector> options
	):
	m_mesh(variables->get_mesh()),
	m_database(variables->get_database()),
	m_bcs(boundary_conditions),
	m_options(options)
	{}

	scalar_t dt() {return m_dt;}

	int compute_one_node(index_t node, UnsaturatedVariables * const vars);
	void compute_total_concentrations(
	index_t node,
	Vector& total_concentrations,
	UnsaturatedVariables * const vars
	);
	void analyse_solution(
	index_t node,
	AdimensionalSystemSolution& solution,
	UnsaturatedVariables * const vars
	);

	void initialize_timestep_one_node(index_t node, UnsaturatedVariables* vars);

	void initialize(UnsaturatedVariables* vars);
	void initialize_timestep(scalar_t dt, UnsaturatedVariables* vars);
	StaggerReturnCode restart_timestep(UnsaturatedVariables* vars);

	units::UnitsSet& get_units() {
	return m_options->get("default").units_set;
	}


	};

	using TrueConstPtr = UnsaturatedVariables * const;

	inline TrueConstPtr cast_to_var(solver::VariablesBase * const var)
	{
	return static_cast<TrueConstPtr>(var);
	}

	// ============================================
	//
	// Implementation
	//
	// ============================================

	EquilibriumStagger::EquilibriumStagger(
	std::shared_ptr<UnsaturatedVariables> variables,
	std::shared_ptr<BoundaryConditions> boundary_conditions,
	std::shared_ptr<OptionsVector> options
	):
	m_impl(utils::make_pimpl<EquilibriumStaggerImpl>(
	variables, boundary_conditions, options))
	{

	}

	EquilibriumStagger::~EquilibriumStagger() = default;

	//! \brief Initialize the stagger at the beginning of the computation
	//!
	//! \param var a shared_ptr to the variables
	void EquilibriumStagger::initialize(VariablesBase * const var)
	{
	TrueConstPtr true_vars = cast_to_var(var);
	m_impl->initialize(true_vars);

	}

	//! \brief Initialize the stagger at the beginning of an iteration
	//!
	//! This is where the predictor can be saved, the first trivial iteration done, ...
	//!
	//! \param dt the duration of the timestep
	//! \param var a shared_ptr to the variables
	void EquilibriumStagger::initialize_timestep(
	scalar_t dt,
	VariablesBase * const var
	)
	{
	TrueConstPtr true_vars = cast_to_var(var);
	m_impl->initialize_timestep(dt, true_vars);
	}

	//! \brief Solve the equation for the timestep
	//!
	//! \param var a shared_ptr to the variables
	EquilibriumStagger::StaggerReturnCode
	EquilibriumStagger::restart_timestep(VariablesBase * const var)
	{
	TrueConstPtr true_vars = cast_to_var(var);
	return m_impl->restart_timestep(true_vars);
	}


	int
	EquilibriumStagger::EquilibriumStaggerImpl::compute_one_node(
	index_t node,
	UnsaturatedVariables * const vars
	)
	{
	AdimensionalSystemConstraints constraints = m_bcs->get_constraint(node);
	compute_total_concentrations(node, constraints.total_concentrations, vars);

	// set partial pressure model
	{
	user_model_saturation_f callable = vars->get_vapor_pressure_model();
	water_partial_pressure_f pv_model = [callable, node](scalar_t sat){
	return callable(node, sat);
	};

	constraints.set_water_partial_pressure_model(pv_model);
	}

	// We use the current value of the saturation to avoid converging
	// to the previous solution
	// This is necessary because the update to the saturation is very small
	const AdimensionalSystemSolution& previous_solution = vars->get_adim_solution(node);
	Vector x;
	AdimensionalSystemSolver the_solver;
	if (likely(previous_solution.is_valid)) // should always be true
	{
	const auto porosity = vars->get_porosity()(node);
	const auto saturation = vars->get_liquid_saturation()(node);

	the_solver = AdimensionalSystemSolver(
	vars->get_database(),
	constraints,
	previous_solution,
	m_options->operator [](node)
	);
	x = previous_solution.main_variables;
	// we use the new value of the saturation to avoid a change
	// to small to be detected by the speciaion solver
	x(0) = porosity*saturation; // x(0) = volume fraction water
	}
	else
	{
	const auto porosity = vars->get_porosity()(node);
	const auto saturation = vars->get_liquid_saturation()(node);

	the_solver = AdimensionalSystemSolver(
	vars->get_database(),
	constraints,
	m_options->operator [](node)
	);
	the_solver.initialise_variables(x, porosity*saturation, -3);
	}
	//
	micpsolver::MiCPPerformance perf = the_solver.solve(x);

	if (perf.return_code < micpsolver::MiCPSolverReturnCode::Success)
	{
	ERROR << "Failed to solve equilibrium at node " << node
	<< ". \n"
	<< "Return code : " << (int) perf.return_code;
	return -1;
	}

	AdimensionalSystemSolution solution = the_solver.get_raw_solution(x);
	analyse_solution(node, solution, vars);
	vars->set_adim_solution(node, solution);

	return 0;
	}

	void
	EquilibriumStagger::EquilibriumStaggerImpl::compute_total_concentrations(
	index_t node,
	Vector& total_concentrations,
	UnsaturatedVariables * const vars
	)
	{
	total_concentrations.resize(m_database->nb_component());
	const scalar_t saturation = vars->get_liquid_saturation()(node);
	const scalar_t porosity = vars->get_porosity()(node);
	const scalar_t rt = vars->get_rt();
	// water
	{
	const scalar_t ctilde_w = vars->get_water_aqueous_concentration()(node);
	const scalar_t cbar_w = vars->get_solid_concentration(0)(node);

	scalar_t c_w = saturationporosityctilde_w + cbar_w;

	if (vars->component_has_gas(0))
	{
	const scalar_t pv_w = vars->get_pressure_main_variables(0)(node);
	c_w += (1.0-saturation)porositypv_w/rt;
	}
	total_concentrations(0) = c_w;
	}
	total_concentrations(1) = 0.0;
	// aqueous components
	for (index_t component: m_database->range_aqueous_component())
	{
	const scalar_t ctilde_i = vars->get_aqueous_concentration(component)(node);
	const scalar_t cbar_i = vars->get_solid_concentration(component)(node);

	scalar_t c_i = saturationporosityctilde_i + cbar_i;

	if (vars->component_has_gas(component))
	{
	const scalar_t pv_i = vars->get_pressure_main_variables(component)(node);
	c_i += (1.0-saturation)porositypv_i/rt;
	}
	total_concentrations(component) = c_i;
	}
	}

	void EquilibriumStagger::EquilibriumStaggerImpl::analyse_solution(
	index_t node,
	AdimensionalSystemSolution& solution,
	UnsaturatedVariables * const vars
	)
	{
	AdimensionalSystemSolutionExtractor extr(solution,
	m_database,
	get_units());

	const scalar_t saturation = extr.saturation_water();
	const scalar_t sat_g = 1.0 - saturation;
	const scalar_t rho_l = extr.density_water();
	const scalar_t porosity = extr.porosity();
	const scalar_t rt = vars->get_rt();
	// porosity

	SecondaryTransientVariable& porosity_vars = vars->get_porosity();
	const scalar_t porosity_0 = porosity_vars.predictor(node);
	const scalar_t porosity_vel = (porosity - porosity_0)/m_dt;

	porosity_vars.variable(node) = porosity;
	porosity_vars.velocity(node) = porosity_vel;

	// water

	// saturation
	MainVariable& satvars = vars->get_liquid_saturation();
	const scalar_t saturation_0 = satvars.predictor(node);
	const scalar_t sat_vel = (saturation - saturation_0)/m_dt;

	satvars.variable(node) = saturation;
	satvars.velocity(node) = sat_vel;

	{
	// total aqueous concentration
	SecondaryTransientVariable& cwtilde_vars = vars->get_water_aqueous_concentration();
	const scalar_t cwtilde = rho_l*extr.total_aqueous_concentration(0);
	const scalar_t cwtilde_0 = cwtilde_vars.predictor(node);
	const scalar_t cwtilde_vel = (cwtilde - cwtilde_0)/m_dt;

	cwtilde_vars.variable(node) = cwtilde;
	cwtilde_vars.velocity(node) = cwtilde_vel;

	// total solid concentration
	MainVariable& solid_vars = vars->get_solid_concentration(0);
	const scalar_t cwbar = extr.total_solid_concentration(0);
	const scalar_t cwbar_0 = solid_vars.predictor(node);
	const scalar_t cwbar_vel = (cwbar - cwbar_0)/m_dt;

	solid_vars.variable(node) = cwbar;
	solid_vars.velocity(node) = cwbar_vel;
	solid_vars.chemistry_rate(node) = - cwbar_vel;

	// vapor pressure
	if (vars->component_has_gas(0))
	{
	MainVariable& pres_vars = vars->get_pressure_main_variables(0);
	const scalar_t pv = vars->get_vapor_pressure_model()(node, saturation);
	const scalar_t pv_0 = pres_vars.predictor(node);
	const scalar_t pv_vel = (pv - pv_0)/m_dt;

	pres_vars.variable(node) = pv;
	pres_vars.velocity(node) = pv_vel;

	const scalar_t transient = ( (pvporositysat_g)
	- (pv_0porosity_0(1.0-saturation_0))
	)/
	(rt*m_dt);
	const scalar_t pv_chem_rate = - transient + pres_vars.transport_fluxes(node);

	pres_vars.chemistry_rate(node) = pv_chem_rate;
	}
	}
	// aqueous components
	for (index_t component: m_database->range_aqueous_component())
	{
	// total aqueous concentration
	MainVariable& aqconc = vars->get_aqueous_concentration(component);
	const scalar_t ctildei = rho_l*extr.total_aqueous_concentration(component);
	const scalar_t ctildei_0 = aqconc.predictor(node);
	const scalar_t ctildei_vel = (ctildei - ctildei_0) / m_dt;

	aqconc.variable(node) = ctildei;
	aqconc.velocity(node) = ctildei_vel;

	// total solid concentration
	MainVariable& solconc = vars->get_solid_concentration(component);
	const scalar_t cbari = extr.total_solid_concentration(component);
	const scalar_t cbari_0 = solconc.predictor(node);
	const scalar_t cbari_vel = (cbari - cbari_0) / m_dt;

	solconc.variable(node) = cbari;
	solconc.velocity(node) = cbari_vel;
	solconc.chemistry_rate(node) = - cbari_vel;

	// partial pressure
	if (vars->component_has_gas(component))
	{
	MainVariable& pres_vars = vars->get_pressure_main_variables(component);
	index_t id_g = vars->get_id_gas(component);
	const scalar_t pi = extr.fugacity_gas(id_g)*vars->get_total_pressure();
	const scalar_t pi_0 = pres_vars.predictor(node);
	const scalar_t pi_vel = (pi - pi_0) / m_dt;

	pres_vars.variable(node) = pi;
	pres_vars.velocity(node) = pi_vel;

	const scalar_t transient = (
	(piporositysat_g)
	- (pi_0porosity_0(1.0-saturation_0))
	) / (rt * m_dt);

	const scalar_t pi_chem_rate = - transient + pres_vars.transport_fluxes(node);
	pres_vars.chemistry_rate(node) = pi_chem_rate;

	}
	}
	}

	void
	EquilibriumStagger::EquilibriumStaggerImpl::initialize_timestep_one_node(
	index_t node,
	UnsaturatedVariables * const vars
	)
	{
	}

	EquilibriumStagger::StaggerReturnCode
	EquilibriumStagger::EquilibriumStaggerImpl::restart_timestep(
	UnsaturatedVariables* vars
	)
	{
	int sum_retcode = 0;
	#if SPECMICP_USE_OPENMP
	#pragma omp parallel for \
	reduction(+: sum_retcode)
	#endif
	for (index_t node=0; node<m_mesh->nb_nodes(); ++node)
	{
	if (not m_bcs->is_normal_node(node)) continue;
	sum_retcode += compute_one_node(node, vars);
	}
	if (sum_retcode != 0)
	{
	return StaggerReturnCode::UnknownError;
	}
	return StaggerReturnCode::ResidualMinimized;
	}

	void
	EquilibriumStagger::EquilibriumStaggerImpl::initialize_timestep(
	scalar_t dt,
	UnsaturatedVariables* vars
	)
	{
	m_dt = dt;
	// initialize
	//water
	{
	MainVariable& cbar_w = vars->get_solid_concentration(0);
	cbar_w.predictor = cbar_w.variable;
	cbar_w.velocity.setZero();
	cbar_w.chemistry_rate.setZero();


	SecondaryTransientVariable& ctilde_w = vars->get_water_aqueous_concentration();
	ctilde_w.predictor = ctilde_w.variable;
	ctilde_w.velocity.setZero();

	if (vars->component_has_gas(0))
	{
	MainVariable& pres_vars = vars->get_pressure_main_variables(0);
	pres_vars.predictor = pres_vars.variable;
	pres_vars.velocity.setZero();
	pres_vars.chemistry_rate.setZero();
	}
	}
	// aqueous component
	for (index_t component: m_database->range_aqueous_component())
	{
	MainVariable& cbar_i = vars->get_solid_concentration(component);
	cbar_i.predictor = cbar_i.variable;
	cbar_i.velocity.setZero();
	cbar_i.chemistry_rate.setZero();


	if (vars->component_has_gas(component))
	{
	MainVariable& pres_vars = vars->get_pressure_main_variables(component);
	pres_vars.predictor = pres_vars.variable;
	pres_vars.velocity.setZero();
	pres_vars.chemistry_rate.setZero();
	}
	}
	// porosity
	{
	SecondaryTransientVariable& porosity = vars->get_porosity();

	porosity.predictor = porosity.variable;
	porosity.velocity.setZero();
	}

	for (auto node: m_mesh->range_nodes())
	{
	initialize_timestep_one_node(node, vars);
	}
	}

	void
	EquilibriumStagger::EquilibriumStaggerImpl::initialize(
	UnsaturatedVariables * const vars
	)
	{
	// do nothing by default
	}



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

equilibrium_stagger.cppNo OneTemporaryActions

File Metadata

equilibrium_stagger.cppView Options

Event Timeline

equilibrium_stagger.cpp
No OneTemporary
Actions

equilibrium_stagger.cpp
View Options