Diffusion SpecMiCP / ReactMiCP (reactmicp/unsaturated) Browse

unsaturated

drying.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.

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


			// ====================================================== //
			// //
			// Drying example //
			// ============== //
			// //
			// ====================================================== //
			//
			//
			// This is a simple example to show how the unsaturated
			// system of ReactMiCP is able to solve a complex drying
			// equation.
			//
			// It implements a special chemistry stagger (DryingStagger)
			// which only computes the vapor pressure.


			// Macro to define the problem
			// ===========================

			#define NB_NODES 50
			#define DX 0.05*1e-2 // cm
			#define CROSS_SECTION 5.0*1e-4 // cm^2

			#define GAMMA_WATER_GLASS 0.1 // N/m
			#define GAMMA_WATER_AIR 71.97e-3 // N/m = kg/s^2
			#define R_SMALL_BEAD 0.015*1e-2 // cm
			#define R_BIG_BEAD 0.04*1e-2 // cm

			#define M_v 18.015e-3 // kg/mol

			// Includes
			// ========

			#include "database/database.hpp"

			#include "reactmicp/systems/unsaturated/variables.hpp"
			#include "reactmicp/systems/unsaturated/variables_sub.hpp"
			#include "reactmicp/systems/unsaturated/transport_stagger.hpp"
			#include "reactmicp/systems/unsaturated/transport_constraints.hpp"

			#include "reactmicp/solver/staggers_base/chemistry_stagger_base.hpp"
			#include "reactmicp/solver/staggers_base/upscaling_stagger_base.hpp"
			#include "reactmicp/solver/staggers_base/stagger_structs.hpp"

			#include "reactmicp/solver/runner.hpp"

			#include "dfpm/meshes/generic_mesh1d.hpp"

			#include "physics/laws.hpp"
			#include "physics/units.hpp"

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

			#include <functional>
			#include <iostream>

			// Namespace
			// =========
			using namespace specmicp;
			using namespace specmicp::reactmicp;
			using namespace specmicp::reactmicp::systems;
			using namespace specmicp::reactmicp::systems::unsaturated;

			// Declarations
			// ============

			// Constants
			// ==========
			static const scalar_t rho_l_si = constants::water_density_25;
			static const scalar_t rho_l = 1e-6*constants::water_density_25;
			static const scalar_t cw_tilde_si = rho_l_si / M_v;
			static const scalar_t cw_tilde = 1e-6*cw_tilde_si;

			int main();

			database::RawDatabasePtr get_database();
			mesh::Mesh1DPtr get_mesh();

			scalar_t leverett_function(scalar_t s_eff);

			scalar_t big_bead_cap_pressure(scalar_t saturation);
			scalar_t small_bead_cap_pressure(scalar_t saturation);
			scalar_t vapor_pressure(scalar_t pc);

			// Drying Stagger
			// ==============
			//
			// This stagger computes the exchange between water vapor and liquid water
			class DryingStagger: public solver::ChemistryStaggerBase
			{
			public:
			DryingStagger() {}

			void initialize(solver::VariablesBasePtr var) {}

			//! \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 initialize_timestep(scalar_t dt, solver::VariablesBasePtr var) {
			m_dt = dt;
			}

			//! \brief Solve the equation for the timestep
			//!
			//! \param var a shared_ptr to the variables
			solver::StaggerReturnCode restart_timestep(solver::VariablesBasePtr var);

			void compute_one_node(index_t node, UnsaturatedVariables* vars);

			private:
			scalar_t m_dt {-1.0};
			};

			// Upscaling Stagger
			// =================
			//
			// Define how we compute capillary pressure, vapor pressure, ...
			class UpscalingDryingStagger: public solver::UpscalingStaggerBase
			{
			public:
			UpscalingDryingStagger(std::vector<bool>& is_big_bead):
			m_is_big_bead(is_big_bead)
			{}

			//! \brief Initialize the stagger at the beginning of the computation
			//!
			//! \param var a shared_ptr to the variables
			void initialize(solver::VariablesBasePtr var);

			//! \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 initialize_timestep(scalar_t dt, solver::VariablesBasePtr var) {}

			//! \brief Solve the equation for the timestep
			//!
			//! \param var a shared_ptr to the variables
			solver::StaggerReturnCode restart_timestep(solver::VariablesBasePtr var) {
			return solver::StaggerReturnCode::ResidualMinimized;
			}


			scalar_t capillary_pressure_model(index_t node, scalar_t saturation);

			scalar_t vapor_pressure_model(index_t node, scalar_t saturation);

			scalar_t relative_liquid_permeability_model(index_t node, scalar_t saturation);

			scalar_t relative_liquid_diffusivity_model(index_t node, scalar_t saturation);

			scalar_t relative_gas_diffusivity_model(index_t node, scalar_t saturation);

			private:
			std::vector<bool> m_is_big_bead;
			};



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

			scalar_t UpscalingDryingStagger::capillary_pressure_model(index_t node, scalar_t saturation)
			{
			scalar_t pc = 0.0;
			if (m_is_big_bead[node])
			{
			pc = big_bead_cap_pressure(saturation);
			}
			else
			pc = small_bead_cap_pressure(saturation);
			return pc;
			}

			scalar_t UpscalingDryingStagger::vapor_pressure_model(index_t node, scalar_t saturation)
			{
			return vapor_pressure(capillary_pressure_model(node, saturation));
			}


			scalar_t UpscalingDryingStagger::relative_liquid_permeability_model(index_t node, scalar_t saturation)
			{
			return std::pow(saturation, 3);
			}

			scalar_t UpscalingDryingStagger::relative_liquid_diffusivity_model(index_t node, scalar_t saturation)
			{
			return saturation;
			}

			scalar_t UpscalingDryingStagger::relative_gas_diffusivity_model(index_t node, scalar_t saturation)
			{
			return (1.0 - saturation);
			}

			void UpscalingDryingStagger::initialize(solver::VariablesBasePtr var)
			{
			UnsaturatedVariables* true_vars = static_cast<UnsaturatedVariables*>(var.get());

			true_vars->set_capillary_pressure_model(std::bind(std::mem_fn(
			&UpscalingDryingStagger::capillary_pressure_model
			), this, std::placeholders::_1, std::placeholders::_2));

			true_vars->set_vapor_pressure_model(std::bind(std::mem_fn(
			&UpscalingDryingStagger::vapor_pressure_model
			), this, std::placeholders::_1, std::placeholders::_2));

			true_vars->set_relative_liquid_permeability_model(std::bind(std::mem_fn(
			&UpscalingDryingStagger::relative_liquid_permeability_model
			), this, std::placeholders::_1, std::placeholders::_2));

			true_vars->set_relative_liquid_diffusivity_model(std::bind(std::mem_fn(
			&UpscalingDryingStagger::relative_liquid_diffusivity_model
			), this, std::placeholders::_1, std::placeholders::_2));

			true_vars->set_relative_gas_diffusivity_model(std::bind(std::mem_fn(
			&UpscalingDryingStagger::relative_gas_diffusivity_model
			), this, std::placeholders::_1, std::placeholders::_2));

			MainVariable& vapor_pressure = true_vars->get_pressure_main_variables(0);
			MainVariable& saturation = true_vars->get_liquid_saturation();

			SecondaryTransientVariable& porosity = true_vars->get_porosity();
			SecondaryVariable& liquid_diffusivity = true_vars->get_liquid_diffusivity();
			SecondaryVariable& liquid_permeability = true_vars->get_liquid_permeability();
			SecondaryVariable& resistance_gas_diffusivity = true_vars->get_resistance_gas_diffusivity();

			true_vars->get_binary_gas_diffusivity(0) = 0.282*1e-4;

			porosity(0) = 1;
			resistance_gas_diffusivity(0) = 1.0;


			for (index_t node=1; node<true_vars->get_mesh()->nb_nodes(); ++node)
			{
			vapor_pressure(node) = vapor_pressure_model(node, saturation(node));

			if (m_is_big_bead[node]) {
			porosity(node) = 0.371;
			liquid_diffusivity(node) = 0;
			liquid_permeability(node) = 3.52e-11;
			resistance_gas_diffusivity(node) = 2*0.371/(3-0.371);
			} else {
			porosity(node) = 0.387;
			liquid_diffusivity(node) = 0;
			liquid_permeability(node) = 8.41e-12;
			resistance_gas_diffusivity(node) = 2*0.387/(3-0.387);
			}
			}

			porosity.predictor = porosity.variable;

			true_vars->set_relative_variables();
			}



			solver::StaggerReturnCode DryingStagger::restart_timestep(solver::VariablesBasePtr var)
			{
			UnsaturatedVariables* true_vars = static_cast<UnsaturatedVariables*>(var.get());

			for (index_t node=1; node<true_vars->get_mesh()->nb_nodes(); ++node) {
			compute_one_node(node, true_vars);
			}

			return solver::StaggerReturnCode::ResidualMinimized;
			}

			void DryingStagger::compute_one_node(index_t node, UnsaturatedVariables *vars)
			{
			MainVariable& satvars = vars->get_liquid_saturation();
			MainVariable& presvars = vars->get_pressure_main_variables(0);


			scalar_t rt = vars->get_rt();
			scalar_t saturation = satvars(node);
			scalar_t pressure = presvars(node);
			scalar_t porosity = vars->get_porosity()(node);
			auto vapor_pressure_f = vars->get_vapor_pressure_model();

			const scalar_t tot_conc = cw_tilde_sisaturation + (1.0-saturation)pressure/rt;

			pressure = vapor_pressure_f(node, saturation);
			scalar_t cw_hat = pressure/rt;
			saturation = (tot_conc - cw_hat) / (cw_tilde_si-cw_hat);

			pressure = vapor_pressure_f(node, saturation);
			cw_hat = pressure/rt;

			scalar_t res = tot_conc - cw_tilde_sisaturation - (1.0-saturation)cw_hat;

			while (std::abs(res)/tot_conc > 1e-12) {
			saturation = (tot_conc - cw_hat) / (cw_tilde_si-cw_hat);
			pressure = vapor_pressure_f(node, saturation);
			cw_hat = pressure/rt;

			res = tot_conc - cw_tilde_sisaturation - (1.0-saturation)cw_hat;
			}

			satvars(node) = saturation;
			satvars.velocity(node) = (saturation - satvars.predictor(node))/m_dt;

			presvars(node) = pressure;
			presvars.velocity(node) = (pressure - presvars.predictor(node))/m_dt;
			presvars.chemistry_rate(node) = - porosity/(rtm_dt) (
			(pressure*(1.0-saturation)) -
			(presvars.predictor(node)*(1.0-satvars.predictor(node)))
			)
			+ presvars.transport_fluxes(node);
			}


			database::RawDatabasePtr get_database()
			{

			specmicp::database::Database thedatabase(CEMDATA_PATH);
			std::vector<std::string> to_keep = {"H2O", "H[+]"};
			thedatabase.keep_only_components(to_keep);
			thedatabase.remove_gas_phases();
			database::RawDatabasePtr raw_data = thedatabase.get_database();
			raw_data->freeze_db();
			return raw_data;
			}

			mesh::Mesh1DPtr get_mesh()
			{
			mesh::Uniform1DMeshGeometry geom;
			geom.dx = DX;
			geom.nb_nodes = NB_NODES;
			geom.section = CROSS_SECTION;
			return mesh::uniform_mesh1d(geom);
			}

			scalar_t big_bead_cap_pressure(scalar_t saturation) {
			return GAMMA_WATER_AIR/std::sqrt(3.52e-11/0.371)*leverett_function(saturation);
			}

			scalar_t small_bead_cap_pressure(scalar_t saturation) {
			return GAMMA_WATER_AIR/std::sqrt(8.41e-12/0.387)*leverett_function(saturation);
			}

			scalar_t vapor_pressure(scalar_t pc) {
			return 3e3std::exp(-M_v pc / (rho_l_si * (8.314*(25+273.15))));
			}


			scalar_t leverett_function(scalar_t s_eff) {
			return 0.325*std::pow((1.0/s_eff - 1), 0.217);
			}


			class Printer
			{
			public:
			Printer(UnsaturatedVariables* vars, std::size_t reserve_size);


			void register_timestep(scalar_t time, solver::VariablesBasePtr base_vars);

			void print(const std::string& name);

			private:
			mesh::Mesh1DPtr m_mesh;

			std::vector<std::vector<scalar_t>> m_pressure;
			std::vector<std::vector<scalar_t>> m_saturation;
			std::vector<scalar_t> m_timestep;
			};

			Printer::Printer(UnsaturatedVariables *vars, std::size_t reserve_size):
			m_mesh(vars->get_mesh()),
			m_pressure(m_mesh->nb_nodes()),
			m_saturation(m_mesh->nb_nodes())
			{
			m_timestep.reserve(reserve_size);
			for (auto node: m_mesh->range_nodes())
			{
			m_pressure[node].reserve(reserve_size);
			m_saturation[node].reserve(reserve_size);
			}
			}

			void Printer::register_timestep(scalar_t time, solver::VariablesBasePtr base_vars)
			{
			UnsaturatedVariables* vars = static_cast<UnsaturatedVariables*>(base_vars.get());

			m_timestep.push_back(time);
			MainVariable& saturation = vars->get_liquid_saturation();
			MainVariable& pressure = vars->get_pressure_main_variables(0);

			for (auto node: m_mesh->range_nodes())
			{
			m_pressure[node].push_back(pressure(node));
			m_saturation[node].push_back(saturation(node));
			}
			}

			void Printer::print(const std::string& name)
			{

			io::CSVFile pressure_file(name+"_pressure.dat");
			io::CSVFile saturation_file(name+"_saturation.dat");

			pressure_file << "Node";
			saturation_file << "Node";
			for (auto& time: m_timestep)
			{
			pressure_file.separator();
			pressure_file << time;
			saturation_file.separator();
			saturation_file << time;
			}
			pressure_file.eol();
			saturation_file.eol();

			for (auto node: m_mesh->range_nodes())
			{
			pressure_file << m_mesh->get_position(node);
			saturation_file << m_mesh->get_position(node);
			for (std::size_t ind=0; ind<m_timestep.size(); ++ind)
			{
			pressure_file.separator();
			pressure_file << m_pressure[node][ind];
			saturation_file.separator();
			saturation_file << m_saturation[node][ind];
			}
			pressure_file.eol();
			saturation_file.eol();
			}

			pressure_file.flush();
			saturation_file.flush();
			}




			int main()
			{
			init_logger(&std::cerr, logger::Warning);

			mesh::Mesh1DPtr the_mesh = get_mesh();
			database::RawDatabasePtr raw_data = get_database();
			std::vector<bool> has_gas = {true, false, false};

			std::vector<bool> is_big_bead(the_mesh->nb_nodes(), true);
			//for (int node=25; node<NB_NODES; ++node)
			// is_big_bead[node] = true;

			UnsaturatedVariablesPtr vars = std::make_shared<UnsaturatedVariables>(
			the_mesh, raw_data, has_gas, units::LengthUnit::meter);


			TransportConstraints constraints;
			constraints.add_fixed_node(0);
			constraints.add_gas_node(0);



			// variables initialisation
			MainVariable& saturation = vars->get_liquid_saturation();
			MainVariable& vapor_pressure = vars->get_pressure_main_variables(0);
			SecondaryTransientVariable& water_aqueous_concentration = vars->get_water_aqueous_concentration();
			water_aqueous_concentration.predictor = water_aqueous_concentration.variable;

			SecondaryVariable& rel_gas_d = vars->get_relative_gas_diffusivity();
			rel_gas_d(0) = 1.0;

			vapor_pressure(0) = 100;

			for (index_t node = 1; node < the_mesh->nb_nodes(); ++node) {
			saturation(node) = 0.9;
			}

			water_aqueous_concentration.set_constant(cw_tilde_si);



			std::shared_ptr<UpscalingDryingStagger> upscaling_stagger = std::make_shared<UpscalingDryingStagger>(is_big_bead);
			upscaling_stagger->initialize(vars);

			std::shared_ptr<DryingStagger> drying_stagger = std::make_shared<DryingStagger>();
			// init vapor pressure
			drying_stagger->initialize_timestep(1.0, vars);
			drying_stagger->restart_timestep(vars);


			std::shared_ptr<UnsaturatedTransportStagger> transport_stagger =
			std::make_shared<UnsaturatedTransportStagger>(vars, constraints);

			solver::ReactiveTransportSolver react_solver(transport_stagger, drying_stagger, upscaling_stagger);
			solver::ReactiveTransportOptions& opts = react_solver.get_options();
			opts.residuals_tolerance = 1e-2;
			opts.good_enough_tolerance = 0.99;

			vapor_pressure.chemistry_rate.setZero();

			//for (auto it=0; it<100; ++it)
			//{
			// auto retcode = react_solver.solve_timestep(0.1, vars);
			// std::cout << (int) retcode << std::endl;
			//}

			//std::cout << saturation.variable << std::endl;
			//std::cout << vapor_pressure.variable << std::endl;
			//std::cout << vars->get_gas_diffusivity().variable << std::endl;
			//std::cout << vars->get_relative_gas_diffusivity().variable << std::endl;

			int nb_hours = 9;
			Printer printer(vars.get(), 4*nb_hours);
			solver::SimulationInformation simul_info("drying", 900);
			solver::ReactiveTransportRunner runner(react_solver, 0.01, 10, simul_info);

			solver::output_f out_func = std::bind(&Printer::register_timestep,
			&printer,
			std::placeholders::_1, std::placeholders::_2);

			runner.set_output_policy(out_func);

			runner.run_until(nb_hours*3600, vars);

			printer.print("out_drying");

			std::cout << saturation.variable << std::endl;
			std::cout << vapor_pressure.variable << std::endl;
			}

specmicp/examples/reactmicp/unsaturated/drying.cppfb013dc5d277reactmicp/unsaturated

drying.cppView Options

specmicp/examples/reactmicp/unsaturated/drying.cpp
fb013dc5d277reactmicp/unsaturated

drying.cpp
View Options