Page MenuHomec4science
Files F62827782

acc_carbo.cpp
No OneTemporary
Actions

File Metadata

Created
Wed, May 15, 22:42

acc_carbo.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.
------------------------------------------------------------------------------*/
//! \file acc_carbo.cpp
//! \brief Accelerated carbonation example
//!
//! This is an example to show the features of the unsaturated
//! reactive transport solver
//! It is a simple study of the accelerated carbonation of cement paste
//! (with drying)
//!
//!
#include "reactmicp_unsaturated.hpp"
#include "reactmicp/systems/unsaturated/variables.hpp"
#include <H5Cpp.h>
#include "reactmicp/io/hdf5_unsaturated.hpp"
#include "utils/io/specmicp_hdf5.hpp"
#include "dfpm/io/hdf5_mesh.hpp"
#include "database/io/hdf5_database.hpp"
#include "dfpmsolver/parabolic_structs.hpp"
#include "physics/unsaturated_laws.hpp"
#include "utils/cli/parser.hpp"
#include <iostream>
using namespace specmicp;
using namespace specmicp::reactmicp;
using namespace specmicp::reactmicp::systems::unsaturated;
using VariablesBase = solver::VariablesBase;
using StaggerReturnCode = solver::StaggerReturnCode;
class UpscalingStagger;
database::RawDatabasePtr get_database();
AdimensionalSystemSolution get_initial_condition(
database::RawDatabasePtr the_database,
const units::UnitsSet& units_set
);
// Upscaling stagger -> user models
class UpscalingStagger: public solver::UpscalingStaggerBase
{
public:
UpscalingStagger() {}
~UpscalingStagger() {}
//! \brief Initialize the stagger at the beginning of the computation
//!
//! \param var a shared_ptr to the variables
void initialize(VariablesBase* 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, VariablesBase* var);
//! \brief Solve the equation for the timestep
//!
//! \param var a shared_ptr to the variables
StaggerReturnCode restart_timestep(VariablesBase* var);
// User models
scalar_t capillary_pressure(index_t node, scalar_t saturation);
scalar_t vapor_pressure(index_t node, scalar_t saturation);
scalar_t relative_liquid_diffusivity(index_t node, scalar_t saturation);
scalar_t relative_liquid_permeability(index_t node, scalar_t saturation);
scalar_t relative_gas_diffusivity(index_t node, scalar_t saturation);
scalar_t intrinsic_permeability(scalar_t porosity);
scalar_t intrinsic_liquid_diffusivity(scalar_t porosity);
scalar_t resistance_gas_diffusivity(scalar_t porosity);
private:
scalar_t m_a {37.5479e6};
scalar_t m_b {2.1684};
scalar_t m_alpha {7.27};
scalar_t m_beta {0.440};
scalar_t m_pv_sat {3170};
scalar_t m_exp_r_l_d {3.00};
scalar_t m_exp_i_g_d {2.74};
scalar_t m_exp_r_g_d {4.20};
scalar_t m_exp_r_l_k {0.4191};
scalar_t m_k_0 {1e-21*1e4};
scalar_t m_d_0 {2.3e-13*1e4};
scalar_t m_phi_0 {0.2};
};
// implementations
int main(int argc, char* argv[])
{
cli::CommandLineParser cli_parser;
cli_parser.add_option('n', "min_dt", 0.01, "minimum timestep (s)");
cli_parser.add_option('m', "max_dt", 10.0, "maximum timestep (s)");
cli_parser.add_option('d', "dx", 0.05, "space step (cm)");
cli_parser.add_option('r', "run_until", 10.0*24.0*3600.0, "run simulation until <value> s");
cli_parser.add_option('p', "p_h20", 1800.0, "vapor pressure at the boundary");
cli_parser.set_help_message("Atmospheric carbonation");
cli_parser.parse(argc, argv);
init_logger(&std::cout, logger::Warning);
index_t nb_nodes = 50;
scalar_t dx = cli_parser.get_option<cli::ValueType::floating>("dx"); // cm
scalar_t cross_section = 5; // cm^2
units::UnitsSet units_set;
units_set.length = units::LengthUnit::centimeter;
specmicp::RawDatabasePtr the_database = get_database();
AdimensionalSystemSolution init = get_initial_condition(
the_database,
units_set);
AdimensionalSystemSolutionExtractor extr(init, the_database, units_set);
std::cout << extr.porosity() << " - " << extr.saturation_water() << std::endl;
//return EXIT_SUCCESS;
mesh::Mesh1DPtr the_mesh = mesh::uniform_mesh1d({dx, nb_nodes, cross_section});
index_t id_co2 = the_database->get_id_component("CO2");
VariablesInterface variable_interface(the_mesh, the_database, {0, id_co2}, units_set);
variable_interface.set_binary_gas_diffusivity(0, 0.240);
variable_interface.set_binary_gas_diffusivity(id_co2, 0.160);
std::shared_ptr<UpscalingStagger> upscaling_stagger =
std::make_shared<UpscalingStagger>();
auto* call_ptr = upscaling_stagger.get();
variable_interface.set_vapor_pressure_model(
[call_ptr](index_t node, scalar_t sat) -> scalar_t {
return call_ptr->vapor_pressure(node, sat);
});
for (index_t node=1; node<the_mesh->nb_nodes(); ++node)
{
variable_interface.initialize_variables(node, extr);
}
variable_interface.set_porosity(0, 1.0);
variable_interface.set_liquid_saturation(0, 0.0);
variable_interface.set_partial_pressure(id_co2, 0, 100.0);
variable_interface.set_partial_pressure(0, 0,
cli_parser.get_option<cli::ValueType::floating>("p_h20"));
std::shared_ptr<UnsaturatedVariables> vars = variable_interface.get_variables();
vars->set_aqueous_scaling(0, 1.0e3/18.3e-3*1e-6);
for (auto aqc: the_database->range_aqueous_component())
{
vars->set_aqueous_scaling(aqc, 100e-6);
}
//vars->set_gaseous_scaling(0, 1.0/vars->get_rt());
//vars->set_gaseous_scaling(id_co2, 1.0/vars->get_rt());
std::cout << "Youpla : " << vars->get_rt() << std::endl;
std::shared_ptr<EquilibriumConstraints> eq_constraints = equilibrium_constraints(the_mesh);
eq_constraints->add_fixed_node(0);
eq_constraints->set_constraint_for_all(eq_constraints->new_constraints());
eq_constraints->get_options().units_set = units_set;
micpsolver::MiCPSolverOptions& copts = eq_constraints->get_options().solver_options;
copts.set_maximum_iterations(200);
copts.set_maximum_step_length(10, 200);
copts.set_tolerance(1e-8, 1e-14);
eq_constraints->get_options().system_options.non_ideality_tolerance = 1e-12;
std::shared_ptr<EquilibriumStagger> chemistry_stagger =
equilibrium_stagger(vars, eq_constraints);
upscaling_stagger->initialize(vars.get());
vars->set_relative_variables(0);
chemistry_stagger->initialize(vars.get());
TransportConstraints transport_constraints;
transport_constraints.add_gas_node(0);
transport_constraints.add_fixed_node(0);
std::shared_ptr<UnsaturatedTransportStagger> transport_stagger =
unsaturated_transport_stagger(vars, transport_constraints, true, true);
transport_stagger->initialize(vars.get());
((*eq_constraints)[0]).set_water_partial_pressure_model(
std::bind(&UpscalingStagger::vapor_pressure,
upscaling_stagger.get(),
1,
std::placeholders::_1)
);
auto check = variable_interface.check_variables();
if (check >= VariablesValidity::error) {
throw std::runtime_error("Invalid variables");
}
solver::ReactiveTransportSolver react_solver(transport_stagger, chemistry_stagger, upscaling_stagger);
solver::ReactiveTransportOptions& opts = react_solver.get_options();
opts.residuals_tolerance = 1e-2;
opts.good_enough_tolerance = 0.99;
opts.maximum_iterations = 100;
transport_stagger->get_saturation_options()->maximum_iterations = 500;
transport_stagger->get_saturation_options()->step_tolerance = 1e-14;
transport_stagger->get_saturation_options()->residuals_tolerance = 1e-8;
//transport_stagger->get_gas_options(3)->maximum_iterations = 5000;
opts.step_tolerance = 1e-10;
io::HDF5File save_file("out_acc_carbo.hdf5", io::HDF5_OpenMode::CreateTruncate);
io::save_database_labels(save_file, "database", "/", the_database);
io::save_mesh_coordinates(save_file, "mesh", "/", the_mesh);
io::UnsaturatedHDF5Saver saver(vars);
saver.save_timepoint(save_file, "0", "/");
save_file.close();
solver::SimulationInformation simul_info("acc_carbo", 900);
auto out_pol = [&saver](
scalar_t timestep,
solver::VariablesBasePtr vars) {
io::HDF5File file("out_acc_carbo.hdf5", io::HDF5_OpenMode::OpenReadWrite);
saver.save_timepoint(file, std::to_string(timestep), "/");
};
//solver::ReactiveTransportReturnCode retcode = react_solver.solve_timestep(0.01, vars);
//retcode = react_solver.solve_timestep(1.0, vars);
//std::cout << (int) retcode << std::endl;
// retcode = react_solver.solve_timestep(5.0, vars);
// std::cout << (int) retcode << std::endl;
// retcode = react_solver.solve_timestep(10.0, vars);
// std::cout << (int) retcode << std::endl;
// retcode = react_solver.solve_timestep(10.0, vars);
// std::cout << (int) retcode << std::endl;
solver::ReactiveTransportRunner runner(
react_solver,
cli_parser.get_option<cli::ValueType::floating>("min_dt"),
cli_parser.get_option<cli::ValueType::floating>("max_dt"),
simul_info);
runner.set_output_policy(out_pol);
runner.run_until(cli_parser.get_option<cli::ValueType::floating>("run_until"),
vars);
std::cout << vars->get_liquid_saturation().variable << std::endl;
std::cout << vars->get_liquid_saturation()(1) - vars->get_liquid_saturation()(2) << std::endl;
std::cout << vars->get_pressure_main_variables(id_co2).variable << std::endl;
}
void UpscalingStagger::initialize(VariablesBase* var)
{
UnsaturatedVariables* true_var = static_cast<UnsaturatedVariables*>(var);
// set relative models
true_var->set_capillary_pressure_model(
[this](index_t x, scalar_t sat){
return this->capillary_pressure(x, sat);
});
true_var->set_relative_liquid_permeability_model(
[this](index_t x, scalar_t sat){
return this->relative_liquid_permeability(x, sat);
});
true_var->set_relative_liquid_diffusivity_model(
[this](index_t x, scalar_t sat){
return this->relative_liquid_diffusivity(x, sat);
});
true_var->set_relative_gas_diffusivity_model(
[this](index_t x, scalar_t sat){
return this->relative_gas_diffusivity(x, sat);
});
true_var->set_vapor_pressure_model(
[this](index_t x, scalar_t sat){
return this->vapor_pressure(x, sat);
});
// initialization
SecondaryTransientVariable& porosity = true_var->get_porosity();
SecondaryVariable& permeability = true_var->get_liquid_permeability();
SecondaryVariable& liq_diffusivity = true_var->get_liquid_diffusivity();
SecondaryVariable& gas_diffusivity = true_var->get_resistance_gas_diffusivity();
true_var->set_relative_variables();
gas_diffusivity(0) = resistance_gas_diffusivity(1.0);
true_var->get_relative_gas_diffusivity()(0) = relative_gas_diffusivity(0, 1.0);
for (index_t node=1; node<true_var->get_mesh()->nb_nodes(); ++node)
{
scalar_t phi = porosity(node);
permeability(node) = intrinsic_permeability(phi);
liq_diffusivity(node) = intrinsic_liquid_diffusivity(phi);
gas_diffusivity(node) = resistance_gas_diffusivity(phi);
//true_var->set_relative_variables(node);
}
}
scalar_t UpscalingStagger::capillary_pressure(index_t node, scalar_t saturation)
{
if (saturation == 0.0) return 0.0;
return laws::capillary_pressure_BaroghelBouny(saturation, m_a, m_b);
}
scalar_t UpscalingStagger::vapor_pressure(index_t node, scalar_t saturation)
{
scalar_t pv;
if (saturation == 0.0) pv = 0.0;
else if (saturation >= 1.0) pv = m_pv_sat;
else pv = m_pv_sat*laws::invert_kelvin_equation(
capillary_pressure(node, saturation));
return pv;
}
scalar_t UpscalingStagger::relative_liquid_diffusivity(
index_t node,
scalar_t saturation
)
{
if (saturation == 0.0) return 0.0;
return std::pow(saturation, m_exp_r_l_d);
}
scalar_t UpscalingStagger::relative_liquid_permeability(
index_t node,
scalar_t saturation
)
{
if (saturation == 0.0) return 0.0;
return laws::relative_liquid_permeability_Mualem(saturation, 1.0/m_b);
}
scalar_t UpscalingStagger::relative_gas_diffusivity(index_t node, scalar_t saturation)
{
if (saturation == 0.0) return 1.0;
return std::pow(1.0 - saturation, m_exp_r_g_d);
}
scalar_t UpscalingStagger::intrinsic_permeability(scalar_t porosity)
{
scalar_t tmp_1 = porosity / m_phi_0;
scalar_t tmp_2 = (1.0 - m_phi_0) / (1.0 - porosity);
tmp_1 = std::pow(tmp_1, 3);
tmp_2 = std::pow(tmp_2, 2);
return m_k_0 * tmp_1 * tmp_2;
}
scalar_t UpscalingStagger::intrinsic_liquid_diffusivity(scalar_t porosity)
{
return m_d_0 * std::exp(-9.95*porosity);
}
scalar_t UpscalingStagger::resistance_gas_diffusivity(scalar_t porosity)
{
return std::pow(porosity, m_exp_i_g_d);
}
void UpscalingStagger::initialize_timestep(scalar_t dt, VariablesBase* var)
{
}
StaggerReturnCode UpscalingStagger::restart_timestep(VariablesBase* var)
{
UnsaturatedVariables* true_var = static_cast<UnsaturatedVariables*>(var);
SecondaryTransientVariable& porosity = true_var->get_porosity();
SecondaryVariable& permeability = true_var->get_liquid_permeability();
SecondaryVariable& liq_diffusivity = true_var->get_liquid_diffusivity();
SecondaryVariable& gas_diffusivity = true_var->get_resistance_gas_diffusivity();
for (index_t node=1; node<true_var->get_mesh()->nb_nodes(); ++node)
{
scalar_t phi = porosity(node);
permeability(node) = intrinsic_permeability(phi);
liq_diffusivity(node) = intrinsic_liquid_diffusivity(phi);
gas_diffusivity(node) = resistance_gas_diffusivity(phi);
true_var->set_relative_variables(node);
}
return StaggerReturnCode::ResidualMinimized;
}
specmicp::RawDatabasePtr get_database()
{
specmicp::database::Database thedatabase(CEMDATA_PATH);
std::map<std::string, std::string> swap = {{"HCO3[-]", "CO2"},};
thedatabase.swap_components(swap);
database::RawDatabasePtr raw_data = thedatabase.get_database();
return raw_data;
}
AdimensionalSystemSolution get_initial_condition(
database::RawDatabasePtr the_database,
const units::UnitsSet& units_set
)
{
Formulation initial_formulation;
scalar_t scaling = 0.9e-3;
initial_formulation.set_mass_solution(0.6*scaling);
initial_formulation.add_mineral("C3S", 5.5*scaling);
initial_formulation.add_mineral("C2S", 1.25*scaling);
initial_formulation.add_mineral("Calcite", 0.05*scaling);
Dissolver dissolver(the_database);
UpscalingStagger upsstag_init; // just for vapor pressure model
AdimensionalSystemConstraints constraints;
constraints.set_total_concentrations(
dissolver.dissolve(initial_formulation, true)
);
//constraints.set_inert_volume_fraction(0.1);
AdimensionalSystemSolver the_solver(the_database, constraints);
the_solver.get_options().units_set = units_set;
the_solver.get_options().system_options.non_ideality_tolerance = 1e-12;
micpsolver::MiCPSolverOptions* solver_options = &the_solver.get_options().solver_options;
//solver_options.maxstep = 10.0;
solver_options->maxiter_maxstep = 200;
solver_options->set_tolerance(1e-9);
Vector x;
the_solver.initialise_variables(x, 0.8, -4);
micpsolver::MiCPPerformance perf = the_solver.solve(x);
if (perf.return_code < micpsolver::MiCPSolverReturnCode::Success)
{
throw std::runtime_error("Error : failed to solve initial conditions");
}
constraints.set_water_partial_pressure_model(
std::bind(&UpscalingStagger::vapor_pressure,
&upsstag_init,
1,
std::placeholders::_1)
);
the_solver = AdimensionalSystemSolver (the_database, constraints);
the_solver.get_options().units_set = units_set;
the_solver.get_options().system_options.non_ideality_tolerance = 1e-12;
solver_options = &the_solver.get_options().solver_options;
//solver_options.maxstep = 10.0;
solver_options->maxiter_maxstep = 200;
solver_options->set_tolerance(1e-9);
perf = the_solver.solve(x);
if (perf.return_code < micpsolver::MiCPSolverReturnCode::Success)
{
throw std::runtime_error("Error : failed to solve initial conditions");
}
return the_solver.get_raw_solution(x);
}

Event Timeline

c4science · Help