Page MenuHomec4science
Files F71671687

leaching_kinetic.cpp
No OneTemporary
Actions

File Metadata

Created
Fri, Jul 12, 13:23

leaching_kinetic.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. *
============================================================================= */
// ============================================
// Leaching of a cement paste
// with kinetic dissolution of C3S
// =============================================
// This file can be used to learn how to use ReactMiCP
// Include ReactMiCP
#include "reactmicp/reactmicp.hpp"
#include "specmicp/io/adimensional_system_solution_saver.hpp"
#include "reactmicp/systems/saturated_react/diffusive_upscaling_stagger.hpp"
#include "reactmicp/systems/saturated_react/kinetic_stagger.hpp"
// Other includes that would be needed later :
#include <Eigen/QR>
#include "specmicp_common/physics/laws.hpp"
#include "specmicp_common/timer.hpp"
#include <functional>
#include <iostream>
#include <fstream>
// we use the following namespaces
using namespace specmicp;
using namespace reactmicp;
using namespace reactmicp::systems;
using namespace specmicp::reactmicp::systems::satdiff;
using namespace reactmicp::solver;
// We declare some functions that we will be using :
// Database
// ========
// Parse and Prepare the database
RawDatabasePtr get_database(const std::string& path_to_database);
// Boundary and initial conditions
// ===============================
// Return the chemistry constraints for the inital condition
AdimensionalSystemConstraints get_cement_initial_constraints(
RawDatabasePtr the_database,
const units::UnitsSet& units_set);
// return the chemistry constraints for the boundary condition
AdimensionalSystemConstraints get_bc_initial_constraints(
RawDatabasePtr the_database,
const units::UnitsSet& units_set);
Vector initial_compo(const Vector& publi);
void compo_from_oxyde(Vector& compo_oxyde, Vector& compo_species);
// The mesh
// ========
// There is two version
// a uniform mesh
mesh::Mesh1DPtr get_mesh(scalar_t dx);
// and a non uniform mesh
mesh::Mesh1DPtr get_mesh();
// Options
// =======
// Speciation solver options
AdimensionalSystemSolverOptions get_specmicp_options();
// Transport options
void set_transport_options(
dfpmsolver::ParabolicDriverOptions& transport_options
);
// Reactive Solver options
void set_reactive_solver_options(
reactmicp::solver::ReactiveTransportOptions& reactive_options
);
// Kinetic model and variables
// ===========================
// Kinetic variables
class C3SKineticVariables : public KineticSaturatedVariablesBase
{
public:
Vector amount_minerals; // concentration of C3S
Vector predictor; // Initial amount
Vector velocity; // dC_c3s/dt
Vector volume_fractions;
};
using C3SKineticVariablesPtr = std::shared_ptr<C3SKineticVariables>;
// The dissolution model
class C3SDissolutionModel : public KineticModelBase
{
public:
C3SDissolutionModel(RawDatabasePtr the_database, const units::UnitsSet the_units):
m_database(the_database),
m_units(the_units),
id_c3s(the_database->get_id_mineral_kinetic("C3S"))
{
if (id_c3s == no_species)
throw database::invalid_database("The database do not contain C3S !");
}
KineticSaturatedVariablesBasePtr initialize_variables(
mesh::Mesh1DPtr the_mesh,
scalar_t init_vol_fraction
);
//! \brief Initialize a timestep
virtual void initialize_timestep(
scalar_t dt,
KineticSaturatedVariablesBasePtr kinetic_variables
) override;
//! \brief Restart a timestep
virtual void restart_timestep(
index_t node,
const AdimensionalSystemSolution& eq_solution,
KineticSaturatedVariablesBasePtr kinetic_variables
) override;
//! \brief Return the volume fraction of kinetic (or inert) phases
virtual scalar_t get_volume_fraction_kinetic(
index_t node,
KineticSaturatedVariablesBasePtr kinetic_variables
) override;
//! \brief Return the velocity of the total immobile kinetic concentration for component
virtual scalar_t get_velocity_kinetic(
index_t node,
index_t component,
KineticSaturatedVariablesBasePtr kinetic_variables
) override;
private:
scalar_t volume_fraction(scalar_t concentration);
scalar_t dy_dt(scalar_t log_IAP, scalar_t Aeff);
scalar_t m_dt {-1.0};
RawDatabasePtr m_database;
units::UnitsSet m_units;
index_t id_c3s;
};
KineticSaturatedVariablesBasePtr C3SDissolutionModel::initialize_variables(
mesh::Mesh1DPtr the_mesh,
scalar_t init_vol_fraction
)
{
C3SKineticVariablesPtr var = std::make_shared<C3SKineticVariables>();
var->predictor.setZero(the_mesh->nb_nodes());
var->velocity.setZero(the_mesh->nb_nodes());
var->volume_fractions.setConstant(the_mesh->nb_nodes(), init_vol_fraction);
var->amount_minerals.setConstant(the_mesh->nb_nodes(),
init_vol_fraction/m_database->molar_volume_mineral_kinetic(id_c3s, m_units));
var->volume_fractions(0) = 0.0;
var->amount_minerals(0) = 0.0;
return std::static_pointer_cast<KineticSaturatedVariablesBase>(var);
}
void C3SDissolutionModel::initialize_timestep(
scalar_t dt,
KineticSaturatedVariablesBasePtr kinetic_variables
)
{
m_dt = dt;
C3SKineticVariablesPtr true_kin_var = std::static_pointer_cast<C3SKineticVariables>(kinetic_variables);
true_kin_var->predictor = true_kin_var->amount_minerals;
true_kin_var->velocity.setZero();
}
void C3SDissolutionModel::restart_timestep(
index_t node,
const AdimensionalSystemSolution &eq_solution,
KineticSaturatedVariablesBasePtr kinetic_variables
)
{
C3SKineticVariablesPtr true_kin_var = std::static_pointer_cast<C3SKineticVariables>(kinetic_variables);
AdimensionalSystemSolutionExtractor extr(eq_solution, m_database, m_units);
scalar_t log_IAP = 0;
for (index_t component: m_database->range_component())
{
if (m_database->nu_mineral_kinetic(id_c3s, component) == 0.0) continue;
log_IAP += m_database->nu_mineral_kinetic(id_c3s, component)*extr.log_molality_component(component);
}
scalar_t tmp = 0;
scalar_t Aeff = 0;
scalar_t porosity = extr.porosity();
if (porosity > 0.6)
{
Aeff= 1400*1e-4/1e-3*m_database->molar_mass_mineral_kinetic(id_c3s)/m_database->molar_volume_mineral_kinetic(id_c3s)*true_kin_var->volume_fractions(node);
tmp = 1e-6*dy_dt(log_IAP, Aeff);
}
true_kin_var->velocity(node) = tmp;
tmp = true_kin_var->predictor(node) + m_dt*true_kin_var->velocity(node);
if (tmp < 1e-20) {
tmp =0.0;
true_kin_var->velocity(node) = - true_kin_var->predictor(node) / m_dt;
}
true_kin_var->amount_minerals(node) = tmp;
true_kin_var->volume_fractions(node) = tmp*m_database->molar_volume_mineral_kinetic(id_c3s, m_units);
}
scalar_t C3SDissolutionModel::dy_dt(scalar_t log_IAP, scalar_t Aeff)
{
scalar_t res = (-57.0 -std::log(10)*log_IAP)/21.5;
scalar_t tmp_2 = std::pow(res, 3.73);
res = - 125.3e-6*Aeff*(1.0 - std::exp(-tmp_2));
return res;
}
scalar_t C3SDissolutionModel::get_volume_fraction_kinetic(
index_t node,
KineticSaturatedVariablesBasePtr kinetic_variables
)
{
C3SKineticVariablesPtr true_kin_var = std::static_pointer_cast<C3SKineticVariables>(kinetic_variables);
return true_kin_var->volume_fractions(node);
}
scalar_t C3SDissolutionModel::get_velocity_kinetic(
index_t node,
index_t component,
KineticSaturatedVariablesBasePtr kinetic_variables
)
{
C3SKineticVariablesPtr true_kin_var = std::static_pointer_cast<C3SKineticVariables>(kinetic_variables);
return m_database->nu_mineral_kinetic(id_c3s, component)*true_kin_var->velocity(node);
}
// ==============
// Output
// ==============
// Many information is generated throughout a computation and we can't save everything
// The reactmicp runner can call a function to output the desired information
// Here the output function is a member function for convenience
class OutputPolicy
{
public:
// Initialise the output
// SimulationInformation contain basic information from the simulation
OutputPolicy(const SimulationInformation& simul_info,
mesh::Mesh1DPtr the_mesh,
RawDatabasePtr the_database,
units::UnitsSet the_units):
m_mesh(the_mesh),
m_database(the_database),
m_units(the_units)
{
// Open the files
// Save the pH
out_c_ph.open(simul_info.complete_filepath("ph", suffix));
// Save the total aqueous concentration of HCO3[-]
out_c_hco3.open(simul_info.complete_filepath("c_hco3", suffix));
// save the total solid concentration of s_ca
out_s_ca.open(simul_info.complete_filepath("s_ca", suffix));
// Save the volume fraction of calcite
out_calcite.open(simul_info.complete_filepath("calcite", suffix));
// Save the porosity
out_poro.open(simul_info.complete_filepath("poro", suffix));
out_kin_vol_frac.open(simul_info.complete_filepath("kin_vol_frac", suffix));
// save useful indices from the database
m_id_hco3 = m_database->get_id_component("HCO3[-]");
m_id_ca = m_database->get_id_component("Ca[2+]");
m_id_calcite = m_database->get_id_mineral("Calcite");
}
void output(scalar_t time, VariablesBasePtr variables)
{
SaturatedVariablesPtr true_var = cast_var_from_base(variables);
out_c_ph << time;
out_c_hco3 << time;
out_s_ca << time;
out_calcite << time;
out_poro << time;
out_kin_vol_frac << time;
for (index_t node: m_mesh->range_nodes())
{
// Again AdimensionalSystemSolutionExtractor is used to obtain information about the speciation system
AdimensionalSystemSolutionExtractor extractor(true_var ->equilibrium_solution(node), m_database, m_units);
out_c_ph << "\t" << extractor.pH();
out_c_hco3 << "\t" << true_var->aqueous_concentration(node, m_id_hco3);
out_s_ca << "\t" << true_var->solid_concentration(node, m_id_ca);
out_calcite << "\t" << extractor.volume_fraction_mineral(m_id_calcite);
out_poro << "\t" << extractor.porosity();
out_kin_vol_frac << "\t" << std::static_pointer_cast<C3SKineticVariables>(true_var->kinetic_variables())->volume_fractions(node);
}
out_c_hco3 << std::endl;
out_c_ph << std::endl;
out_s_ca << std::endl;
out_calcite << std::endl;
out_poro << std::endl;
out_kin_vol_frac << std::endl;
}
private:
mesh::Mesh1DPtr m_mesh;
RawDatabasePtr m_database;
units::UnitsSet m_units;
std::string suffix { "dat" };
std::ofstream out_c_ph;
std::ofstream out_c_hco3;
std::ofstream out_s_ca;
std::ofstream out_calcite;
std::ofstream out_poro;
std::ofstream out_kin_vol_frac;
index_t m_id_hco3;
index_t m_id_ca;
index_t m_id_calcite;
};
// ===================
// MAIN
// ===================
int main()
{
// initialize eigen
// This is needed if OpenMP is used
Eigen::initParallel();
// Setup the log
// The output will be on cerr
specmicp::logger::ErrFile::stream() = &std::cerr;
// And only Warning and error messages will be printed
specmicp::stdlog::ReportLevel() = specmicp::logger::Warning;
// The database
// -------------
// Obtain the database
RawDatabasePtr the_database = get_database("../data/cemdata.yaml");
// SpecMiCP options
// ----------------
// Set the options
AdimensionalSystemSolverOptions specmicp_options = get_specmicp_options();
units::UnitsSet units_set = specmicp_options.units_set;
// The mesh
// --------
// Obtain the mesh
mesh::Mesh1DPtr the_mesh = get_mesh(0.0050); // uniform
//mesh::Mesh1DPtr the_mesh = get_mesh(); // non uniform
index_t nb_nodes = the_mesh->nb_nodes(); // just boilerplate to make things easier
// Print the mesh
io::print_mesh("out_carbo_mesh.dat", the_mesh, units::LengthUnit::centimeter);
// Initial and boundary conditions
// --------------------------------
// Obtain the initial constraints
// Note : the database is modified at this point
// Only the relevant components are kept at this point
AdimensionalSystemConstraints cement_constraints = get_cement_initial_constraints(the_database, units_set);
// Freeze the database
the_database->freeze_db();
// Obtain the initial condition
AdimensionalSystemSolution cement_solution = solve_equilibrium(
the_database,
cement_constraints,
specmicp_options);
if(not the_database->is_valid())
throw std::runtime_error("The database has been changed during the initial condition computation");
// Just to check that everything is ok !
AdimensionalSystemSolutionExtractor extractor(cement_solution, the_database, units_set);
std::cout << "Porosity : " << extractor.porosity() << std::endl;
std::cout << "Water volume fraction " << extractor.volume_fraction_water() << std::endl;
//return EXIT_SUCCESS;
// Obtain the boundary conditions
AdimensionalSystemConstraints bc_constraints = get_bc_initial_constraints(
the_database,
specmicp_options.units_set
);
AdimensionalSystemSolution bc_solution = solve_equilibrium(
the_database,
bc_constraints,
specmicp_options
);
if(not the_database->is_valid())
throw std::runtime_error("The database has been changed during the boundary conditions computation");
// The constraints
// --------------
// Boundary condition
// 'list_fixed_nodes' lists the node with fixed composition
std::vector<index_t> list_fixed_nodes = {0};
// list_constraints lists the different speciation constraints in the system
// Note : the total concentrations are computed from the initial solution,
// and thus different total concentrations do not constitute a reason to create different constraints
std::vector<specmicp::AdimensionalSystemConstraints> list_constraints = {cement_constraints,
bc_constraints};
// index_constraints links a node to the set of constraints that will be used
std::vector<int> index_constraints(nb_nodes, 0);
index_constraints[0] = 1; // Boundary conditions
// This is the list of initial composition
std::vector<specmicp::AdimensionalSystemSolution> list_initial_composition = {cement_solution, bc_solution};
// This list links a node to its initial conditions
std::vector<int> index_initial_state(nb_nodes, 0);
index_initial_state[0] = 1; // boundary condition
// Variables
// -----------
// Create the main set of variables
// They will be initialized using the list of initial composition
satdiff::SaturatedVariablesPtr variables =
satdiff::init_variables(the_mesh, the_database, units_set,
list_fixed_nodes,
list_initial_composition,
index_initial_state);
std::shared_ptr<C3SDissolutionModel> kinetic_model = std::make_shared<C3SDissolutionModel>(the_database, units_set);
variables->kinetic_variables() = kinetic_model->initialize_variables(the_mesh, 0.2);
// Staggers
// ---------
// This section initialize the staggers
// The equilibrium stagger
// - - - - - - - - - - - -
std::shared_ptr<KineticStagger> kinetic_stagger =
std::make_shared<satdiff::KineticStagger>(kinetic_model, list_constraints, index_constraints, specmicp_options);
// The transport stagger
// - - - - - - - - - - -
std::shared_ptr<SaturatedTransportStagger> transport_stagger =
std::make_shared<satdiff::SaturatedTransportStagger>(variables, list_fixed_nodes);
set_transport_options(transport_stagger->options_solver());
// The upscaling stagger
// - - - - - - - - - - -
CappedPowerLawParameters upscaling_param;
upscaling_param.d_eff_0 = 2.726e-12*1e4;
upscaling_param.cap = 1e10*1e6;
upscaling_param.exponent = 3.32;
upscaling_param.porosity_0 = 0.25;
upscaling_param.porosity_res = 0.02;
CappedPowerLaw upscaling_law(upscaling_param);
UpscalingStaggerPtr upscaling_stagger =
std::make_shared<DiffusiveUpscalingStagger>(upscaling_law.get_law(),the_database, units_set);
// Solver
// ------
// This section initialize the reactive transport solver
ReactiveTransportSolver solver(transport_stagger, kinetic_stagger, upscaling_stagger);
set_reactive_solver_options(solver.get_options());
transport_stagger->initialize(variables);
kinetic_stagger->initialize(variables);
upscaling_stagger->initialize(variables);
// Some basic information about the simulation
// First argument : the name of the simulation
// Second argument : the output function will be called every 3600s
SimulationInformation simul_info("kincarbo", 1800);
simul_info.output_prefix = "out_kincarbo_";
// Runner
// -------
// Create the runner : loop through the timestep
ReactiveTransportRunner runner(solver, 0.1, 100.0, simul_info);
// Create the output class
OutputPolicy output_policy(simul_info, the_mesh, the_database, specmicp_options.units_set);
// and bind it to the runner
runner.set_output_policy(std::bind(&OutputPolicy::output, &output_policy, std::placeholders::_1, std::placeholders::_2));
// Solve the problem
runner.run_until(10*24*3600, cast_var_to_base(variables));
if(not the_database->is_valid())
throw std::runtime_error("The database has been changed during the computation");
io::AdimensionalSystemSolutionSaver saver(the_database);
saver.save_solutions(simul_info.complete_filepath("solutions", "dat"), variables->equilibrium_solutions());
// output information at the end of the timestep
io::print_minerals_profile(the_database, variables, the_mesh, specmicp_options.units_set,
simul_info.complete_filepath("profiles_end", "dat"));
io::print_components_total_aqueous_concentration(the_database, variables,
the_mesh, specmicp_options.units_set,
simul_info.complete_filepath("c_profiles_end", "dat"));
io::print_components_total_solid_concentration(the_database, variables,
the_mesh, specmicp_options.units_set,
simul_info.complete_filepath("s_profiles_end", "dat"));
}
AdimensionalSystemSolverOptions get_specmicp_options()
{
AdimensionalSystemSolverOptions options;
options.solver_options.maxiter_maxstep = 100;
options.solver_options.maxstep = 20.0;
options.solver_options.threshold_cycling_linesearch = 1e-6;
options.solver_options.set_tolerance(1e-8, 1e-14);
options.solver_options.disable_condition_check();
options.solver_options.disable_descent_direction();
options.solver_options.enable_non_monotone_linesearch();
options.solver_options.enable_scaling();
options.system_options.non_ideality_tolerance = 1e-14;
options.system_options.cutoff_total_concentration = 1e-10;
options.units_set.length = units::LengthUnit::centimeter;
return options;
}
void set_transport_options(dfpmsolver::ParabolicDriverOptions& transport_options)
{
transport_options.maximum_iterations = 450;
transport_options.quasi_newton = 3;
transport_options.residuals_tolerance = 1e-8;
transport_options.step_tolerance = 1e-14;
transport_options.absolute_tolerance = 1e-18;
transport_options.alpha = 1.0;
//transport_options.linesearch = dfpmsolver::ParabolicLinesearch::Strang;
transport_options.sparse_solver = SparseSolver::SparseLU;
//transport_options.sparse_solver = SparseSolver::GMRES;
transport_options.threshold_stationary_point = 1e-4;
}
void set_reactive_solver_options(reactmicp::solver::ReactiveTransportOptions& reactive_options)
{
reactive_options.maximum_iterations = 100; //500;
reactive_options.residuals_tolerance = 1e-2;
reactive_options.good_enough_tolerance = 1.0;
reactive_options.absolute_residuals_tolerance = 1e-14;
reactive_options.implicit_upscaling = true;
}
RawDatabasePtr get_database(const std::string& path_to_database) {
database::Database the_database(path_to_database);
std::map<std::string, std::string> swapping({
{"H[+]","HO[-]"},
{"Si(OH)4", "SiO(OH)3[-]"},
{"Al[3+]","Al(OH)4[-]"}
});
the_database.swap_components(swapping);
the_database.remove_gas_phases();
return the_database.get_database();
}
AdimensionalSystemConstraints get_cement_initial_constraints(
RawDatabasePtr the_database,
const units::UnitsSet& units_set)
{
Vector oxyde_compo(4);
oxyde_compo << 66.98, 21.66, 7.19, 2.96;
Vector species_compo;
compo_from_oxyde(oxyde_compo, species_compo);
scalar_t mult = 1e-6;
species_compo *= 0.85e-2;
Formulation formulation;
formulation.mass_solution = 1.0;
formulation.amount_minerals = {
{"C3S", species_compo(0)},
{"C2S", species_compo(1)},
{"C3A", species_compo(2)},
{"Gypsum", species_compo(3)},
{"Calcite", species_compo(0)*1e-3}
};
formulation.extra_components_to_keep = {"HCO3[-]"};
formulation.concentration_aqueous =
{
{"Na[+]", mult*1.0298},
{"K[+]", mult*0.0801},
{"Cl[-]", mult*0.0001},
{"HO[-]", mult*1.8298}
},
formulation.minerals_to_keep = {
"Portlandite",
"CSH,jennite",
"CSH,tobermorite",
"SiO2(am)",
"Calcite",
"Al(OH)3(mic)",
"C3AH6",
"C3ASO_41H5_18",
"C3ASO_84H4_32",
"C4AH19",
"C4AH13",
"C2AH7_5",
"CAH10",
"Monosulfoaluminate",
"Tricarboaluminate",
"Monocarboaluminate",
"Hemicarboaluminate",
"Straetlingite",
"Gypsum",
"Ettringite",
//"Thaumasite"
};
Dissolver dissolver(the_database);
dissolver.set_units(units_set);
AdimensionalSystemConstraints constraints;
constraints.set_charge_keeper(1);
constraints.total_concentrations = dissolver.dissolve(formulation);
std::cout << constraints.total_concentrations << std::endl;
constraints.set_saturated_system();
constraints.set_inert_volume_fraction(0.2);
return constraints;
}
AdimensionalSystemConstraints get_bc_initial_constraints(
RawDatabasePtr the_database,
const units::UnitsSet& units_set)
{
const scalar_t rho_w = laws::density_water(units::celsius(25.0), units_set.length, units_set.mass);
const scalar_t nacl = 0.3*rho_w;
const scalar_t kcl = 0.28*rho_w;
Vector total_concentrations;
total_concentrations.setZero(the_database->nb_component());
total_concentrations(the_database->get_id_component("K[+]")) = kcl;
total_concentrations(the_database->get_id_component("Cl[-]")) = nacl + kcl;
total_concentrations(the_database->get_id_component("Na[+]")) = nacl;
total_concentrations(the_database->get_id_component("SiO(OH)3[-]")) = 0.0001*rho_w;
AdimensionalSystemConstraints constraints;
constraints.add_fixed_activity_component(the_database->get_id_component("HO[-]"), -10.3);
constraints.set_charge_keeper(the_database->get_id_component("HCO3[-]"));
constraints.total_concentrations = total_concentrations;
constraints.set_saturated_system();
constraints.disable_surface_model();
constraints.fixed_fugacity_cs ={};
return constraints;
}
AdimensionalSystemSolution get_bc_initial_compo(
RawDatabasePtr the_database,
const AdimensionalSystemConstraints& constraints,
const AdimensionalSystemSolverOptions& options
)
{
Vector variables;
AdimensionalSystemSolver solver(the_database,
constraints,
options);
micpsolver::MiCPPerformance perf = solver.solve(variables, true);
if (perf.return_code <= micpsolver::MiCPSolverReturnCode::NotConvergedYet)
throw std::runtime_error("Failed to solve initial composition");
std::cout << variables << std::endl;
return solver.get_raw_solution(variables);
}
void compo_from_oxyde(Vector& compo_oxyde, Vector& compo_species)
{
constexpr double M_CaO = 56.08;
constexpr double M_SiO2 = 60.09;
constexpr double M_Al2O3 = 101.96;
constexpr double M_SO3 = 80.06;
Eigen::MatrixXd compo_in_oxyde(4, 4);
Vector molar_mass(4);
molar_mass << 1.0/M_CaO, 1.0/M_SiO2, 1.0/M_Al2O3, 1.0/M_SO3;
compo_oxyde = compo_oxyde.cwiseProduct(molar_mass);
//C3S C2S C3A Gypsum
compo_in_oxyde << 3, 2, 3, 1, // CaO
1, 1, 0, 0, // SiO2
0, 0, 1, 0, // Al2O3
0, 0, 0, 1; // SO3
Eigen::ColPivHouseholderQR<Eigen::MatrixXd> solver(compo_in_oxyde);
compo_species = solver.solve(compo_oxyde);
}
Vector initial_compo(const Vector& publi)
{
Matrix publi_stochio(5,5);
publi_stochio << 1, 0, 0, 0, 0,
9, 6, 0, 0, 0,
3, 0, 1, 1, 0,
3, 0, 1, 3, 0,
0, 0, 0, 0, 1;
Matrix cemdata_stochio(5, 5);
cemdata_stochio << 1, 0, 0, 0, 0,
1.0/0.6, 1.0, 0, 0, 0,
3, 0, 1, 1, 0,
3, 0, 1, 3, 0,
0, 0, 0, 0, 1;
Vector oxyde;
Eigen::ColPivHouseholderQR<Matrix> solver(publi_stochio);
oxyde = solver.solve(publi);
Vector for_cemdata = cemdata_stochio*oxyde;
return for_cemdata;
}
mesh::Mesh1DPtr get_mesh(scalar_t dx)
{
mesh::UniformAxisymmetricMesh1DGeometry geometry;
geometry.dx = dx;
geometry.radius = 0.75/2.0 + dx; //cm
geometry.height = 20.0;
geometry.nb_nodes = geometry.radius/dx + 1;
return mesh::uniform_axisymmetric_mesh1d(geometry);
}
mesh::Mesh1DPtr get_mesh()
{
Vector radius(71);
const scalar_t height = 20;
radius <<
3.751,
3.750,
3.745,
3.740,
3.735,
3.730,
3.720,
3.710,
3.700,
3.690,
3.680,
3.670,
3.660,
3.645,
3.630,
3.615,
3.600,
3.580,
3.560,
3.540,
3.520,
3.500,
3.475,
3.450,
3.425,
3.400,
3.375,
3.350,
3.325,
3.300,
3.275,
3.250,
3.225,
3.200,
3.175,
3.150,
3.125,
3.100,
3.050,
3.025,
3.000,
2.950,
2.900,
2.850,
2.800,
2.750,
2.700,
2.650,
2.600,
2.550,
2.500,
2.450,
2.400,
2.350,
2.300,
2.250,
2.200,
2.150,
2.100,
2.050,
2.000,
1.950,
1.900,
1.850,
1.800,
1.750,
1.700,
1.650,
1.600,
1.550,
1.500
;
radius = radius/10;
return mesh::axisymmetric_mesh1d(radius, height);
}

Event Timeline

c4science · Help