Diffusion SpecMiCP / ReactMiCP (yaml/safe) Browse

safe

equilibrium_stagger.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 "equilibrium_stagger.hpp"
			#include "variables.hpp"
			#include "equilibrium_constraints.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
			{
			scalar_t m_dt {-1};
			mesh::Mesh1DPtr m_mesh;
			database::RawDatabasePtr m_database;
			std::shared_ptr<EquilibriumConstraints> m_constraints;

			EquilibriumStaggerImpl(
			std::shared_ptr<UnsaturatedVariables> variables,
			std::shared_ptr<EquilibriumConstraints> constraints
			):
			m_mesh(variables->get_mesh()),
			m_database(variables->get_database()),
			m_constraints(constraints)
			{}

			scalar_t dt() {return m_dt;}

			int compute_one_node(index_t node, UnsaturatedVariables* vars);
			void compute_total_concentrations(
			index_t node,
			Vector& total_concentrations,
			UnsaturatedVariables* vars
			);
			void analyse_solution(
			index_t node,
			AdimensionalSystemSolution& solution,
			UnsaturatedVariables* 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_constraints->get_options().units_set;
			}


			};

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

			EquilibriumStagger::EquilibriumStagger(
			std::shared_ptr<UnsaturatedVariables> variables,
			std::shared_ptr<EquilibriumConstraints> constraints
			):
			m_impl(make_unique<EquilibriumStaggerImpl>(variables, constraints))
			{

			}

			EquilibriumStagger::~EquilibriumStagger() = default;

			//! \brief Initialize the stagger at the beginning of the computation
			//!
			//! \param var a shared_ptr to the variables
			void EquilibriumStagger::initialize(VariablesBasePtr var)
			{
			UnsaturatedVariables* true_vars =
			static_cast<UnsaturatedVariables*>(var.get());
			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, VariablesBasePtr var)
			{
			UnsaturatedVariables* true_vars =
			static_cast<UnsaturatedVariables*>(var.get());
			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(VariablesBasePtr var)
			{
			UnsaturatedVariables* true_vars =
			static_cast<UnsaturatedVariables*>(var.get());
			return m_impl->restart_timestep(true_vars);
			}


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

			{
			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
			AdimensionalSystemSolution& previous_solution = vars->get_adim_solution(node);
			if (previous_solution.is_valid)
			{
			const auto porosity = vars->get_porosity()(node);
			const auto saturation = vars->get_liquid_saturation()(node);
			previous_solution.main_variables(0) = porosity*saturation;
			}
			//
			AdimensionalSystemSolver the_solver(vars->get_database(),
			constraints,
			previous_solution,
			m_constraints->get_options()
			);
			Vector x = previous_solution.main_variables;
			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* 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* vars
			)
			{
			AdimensionalSystemSolutionExtractor extr(solution,
			m_database,
			get_units());

			const scalar_t saturation = extr.saturation_water();
			const scalar_t sat_g = 1-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-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* vars
			)
			{
			SecondaryTransientVariable& porosity = vars->get_porosity();

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


			}

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


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

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



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

specmicp/src/reactmicp/systems/unsaturated/equilibrium_stagger.cppe56213af5020yaml/safe

equilibrium_stagger.cppView Options

specmicp/src/reactmicp/systems/unsaturated/equilibrium_stagger.cpp
e56213af5020yaml/safe

equilibrium_stagger.cpp
View Options