drying.cpp
No OneTemporary
Actions

Subscribers

None

File Metadata

Created: Sat, Nov 30, 21:35

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

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


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

	virtual void initialize(solver::VariablesBase * const var) override {}

	//! \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
	virtual void initialize_timestep(
	scalar_t dt,
	solver::VariablesBase * const var
	) override {
	m_dt = dt;
	}

	//! \brief Solve the equation for the timestep
	//!
	//! \param var a shared_ptr to the variables
	virtual solver::StaggerReturnCode restart_timestep(
	solver::VariablesBase * const var) override;

	void compute_one_node(index_t node, UnsaturatedVariables * const 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
	virtual void initialize(solver::VariablesBase * const var) override;

	//! \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
	virtual void initialize_timestep(
	scalar_t dt,
	solver::VariablesBase * const var
	) override {}

	//! \brief Solve the equation for the timestep
	//!
	//! \param var a shared_ptr to the variables
	virtual solver::StaggerReturnCode restart_timestep(
	solver::VariablesBase * const var) override
	{
	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::VariablesBase * const var)
	{
	UnsaturatedVariables * const true_vars = static_cast<UnsaturatedVariables * const>(var);

	true_vars->set_capillary_pressure_model(
	[this](index_t node, scalar_t sat) -> scalar_t {
	return this->capillary_pressure_model(node, sat);});

	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::VariablesBase * const var)
	{
	UnsaturatedVariables * const true_vars = static_cast<UnsaturatedVariables * const>(var);

	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 * const 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.get());

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


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

drying.cppNo OneTemporaryActions

File Metadata

drying.cppView Options

Event Timeline

drying.cpp
No OneTemporary
Actions

drying.cpp
View Options