Page Menu
Home
c4science
Search
Configure Global Search
Log In
Files
F91722399
leaching.cpp
No One
Temporary
Actions
Download File
Edit File
Delete File
View Transforms
Subscribe
Mute Notifications
Award Token
Subscribers
None
File Metadata
Details
File Info
Storage
Attached
Created
Wed, Nov 13, 20:42
Size
28 KB
Mime Type
text/x-c
Expires
Fri, Nov 15, 20:42 (2 d)
Engine
blob
Format
Raw Data
Handle
22313829
Attached To
rSPECMICP SpecMiCP / ReactMiCP
leaching.cpp
View Options
#include <iostream>
#include "database/database.hpp"
#include "specmicp/reaction_path.hpp"
#include "dfpm/mesh.hpp"
#include "reactmicp/systems/saturated_diffusion/reactive_transport_solver.hpp"
#include "reactmicp/systems/saturated_diffusion/transport_parameters.hpp"
#include "reactmicp/systems/saturated_diffusion/boundary_conditions.hpp"
#include "reactmicp/systems/saturated_diffusion/reactive_transport_neutrality_solver.hpp"
#include "reactmicp/systems/saturated_diffusion/transport_neutrality_parameters.hpp"
#include "utils/moving_average.hpp"
#include <fstream>
#include <Eigen/LU>
using namespace specmicp;
using namespace specmicp::reactmicp;
using namespace specmicp::reactmicp::systems::siasaturated;
const double M_CaO = 56.08;
const double M_SiO2 = 60.09;
const double M_Al2O3 = 101.96;
const double M_SO3 = 80.06;
void compo_from_oxyde(Vector& compo_oxyde, Vector& compo_species)
{
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::FullPivLU<Eigen::MatrixXd> solver(compo_in_oxyde);
compo_species = solver.solve(compo_oxyde);
}
std::shared_ptr<database::DataContainer> set_database() {
specmicp::database::Database database("../data/cemdata_specmicp.js");
std::map<std::string, std::string> swapping ({
{"H[+]","HO[-]"},
{"Si(OH)4", "SiO(OH)3[-]"},
{"Al[3+]","Al(OH)4[-]"}
});
database.swap_components(swapping);
//std::vector<std::string> to_keep = {"H2O", "HO[-]", "Ca[2+]", "SiO(OH)3[-]"};
//database.keep_only_components(to_keep);
return database.get_database();
}
EquilibriumState get_initial_state(
std::shared_ptr<specmicp::database::DataContainer> database,
Vector& compo,
double mult)
{
std::shared_ptr<specmicp::ReactionPathModel> model = std::make_shared<specmicp::ReactionPathModel>();
compo = mult*compo;
const double m_c3s = compo(0);
const double m_c2s = compo(1);
const double m_c3a = compo(2);
const double m_gypsum = compo(3);
const double wc = 1.0;
const double m_water = wc*( (3*M_CaO+M_SiO2)*m_c3s
+ (2*M_CaO+M_SiO2)*m_c2s
+ (3*M_CaO+M_Al2O3)*m_c3a
+ (M_CaO + M_SO3)*m_gypsum
)*1e-3;
model->amount_aqueous = {
{"H2O", specmicp::reaction_amount_t(m_water/database->molar_mass_basis_si(0), 0.0)},
{"HO[-]", specmicp::reaction_amount_t(0.0, 0.0)},
};
model->amount_minerals = {
{"C3S", specmicp::reaction_amount_t(m_c3s, 0.0)},
{"C2S", specmicp::reaction_amount_t(m_c2s, 0.0)},
{"C3A", specmicp::reaction_amount_t(m_c3a, 0.0)},
{"Gypsum", specmicp::reaction_amount_t(m_gypsum, 0.0)},
};
model->minerals_to_keep = {
"Portlandite",
"CSH,jennite",
"CSH,tobermorite",
"SiO2,am",
"Al(OH)3,am",
"Monosulfoaluminate",
"Straetlingite",
"Gypsum",
"Ettringite"
};
Eigen::VectorXd x;
specmicp::ReactionPathDriver driver(model, database, x);
//driver.get_options().solver_options.penalization_factor =1.0;
driver.get_options().ncp_function = specmicp::micpsolver::NCPfunction::penalizedFB;
driver.get_options().system_options.non_ideality = true;
//driver.get_options().ncp_function = specmicp::micpsolver::NCPfunction::min;
driver.get_options().solver_options.use_scaling = false;
driver.get_options().solver_options.max_factorization_step = 4;
driver.get_options().solver_options.factor_gradient_search_direction = 200;
driver.get_options().solver_options.maxstep = 50;
driver.get_options().solver_options.maxiter_maxstep = 50;
driver.get_options().solver_options.max_iter = 100;
driver.get_options().allow_restart = true;
specmicp::micpsolver::MiCPPerformance perf = driver.one_step(x);
if (perf.return_code != specmicp::micpsolver::MiCPSolverReturnCode::ResidualMinimized)
throw std::runtime_error("Error - system cannot be solved");
return driver.get_current_solution();
}
EquilibriumState get_blank_state(std::shared_ptr<specmicp::database::DataContainer> database)
{
std::shared_ptr<specmicp::ReactionPathModel> model = std::make_shared<specmicp::ReactionPathModel>();
//const double m_c3s = 1e-8;
const double m_c2s = 0.1e-8;
const double m_c3a = 0.1e-8;
const double m_gypsum = 0.1e-8;
const double m_water = 1.0;
model->amount_aqueous = {
{"H2O", specmicp::reaction_amount_t(m_water/database->molar_mass_basis_si(0), 0.0)},
{"HO[-]", specmicp::reaction_amount_t(-0.0, 0.0)},
{"SiO(OH)3[-]", specmicp::reaction_amount_t(0.0, 0.0)}
};
model->amount_minerals = {
//{"C3S", specmicp::reaction_amount_t(m_c3s, 0.0)},
{"C2S", specmicp::reaction_amount_t(m_c2s, 0.0)},
{"C3A", specmicp::reaction_amount_t(m_c3a, 0.0)},
{"Gypsum", specmicp::reaction_amount_t(m_gypsum, 0.0)},
{"SiO2,am", specmicp::reaction_amount_t(0.0, 0.0)}
};
Eigen::VectorXd x(database->nb_component+database->nb_mineral);
x.setZero();
x(0) = 1.0;
x.segment(1, database->nb_component-1).setConstant(-6);
specmicp::ReactionPathDriver driver(model, database, x);
driver.get_options().solver_options.penalization_factor =1.0;
driver.get_options().ncp_function = specmicp::micpsolver::NCPfunction::penalizedFB;
//driver.get_options().ncp_function = specmicp::micpsolver::NCPfunction::min;
driver.get_options().solver_options.use_scaling = false;
driver.get_options().solver_options.max_factorization_step = 1;
driver.get_options().solver_options.factor_gradient_search_direction = 50;
driver.get_options().solver_options.maxstep = 50;
driver.get_options().solver_options.maxiter_maxstep = 50;
driver.get_options().solver_options.max_iter = 100;
driver.get_options().allow_restart = true;
specmicp::micpsolver::MiCPPerformance perf = driver.one_step(x);
if (perf.return_code != specmicp::micpsolver::MiCPSolverReturnCode::ResidualMinimized)
throw std::runtime_error("Error - blank system cannot be solved");
EquilibriumState solution = driver.get_current_solution();
std::cout << solution.main_variables() << std::endl;
// solution.remove_minerals();
// Eigen::VectorXd variables(database->nb_component+database->nb_mineral);
// variables.setZero();
// variables(0) = 1.0;
// for (index_t component = 1; component < database->nb_component; ++component)
// variables(component) = -11;
// scalar_t ionic_strength = 0.0;
// Eigen::VectorXd secondary(database->nb_aqueous);
// Eigen::VectorXd loggamma(database->nb_component+database->nb_aqueous);
// secondary.setZero();
// loggamma.setZero();
// EquilibriumState solution(variables, secondary, loggamma, ionic_strength, database);
return solution;
}
double get_dt(
std::shared_ptr<SaturatedDiffusionTransportParameters> parameters,
std::shared_ptr<mesh::Mesh1D> themesh
)
{
const scalar_t max_d = parameters->diffusion_coefficient(1);
scalar_t dt = std::pow(themesh->get_dx(0), 2)/max_d/3;
if (dt > 5.0) return 5.0;
return 5.0;
//return 10.0;
//return dt;
//if (dt<3600*24) dt=3600*24;
//if (dt>2*3600*24) dt=2*3600*24;
//return 900;
}
int leaching()
{
Vector compo_oxyde(4), compo_species(4);
compo_oxyde << 62.9, 20.6, 5.8, 3.1;
compo_from_oxyde(compo_oxyde, compo_species);
std::shared_ptr<database::DataContainer> thedatabase = set_database();
scalar_t radius = 3.5; //cm
scalar_t height = 8.0; //cm
scalar_t dx = 0.0050;
index_t additional_nodes = 1;
radius = radius + additional_nodes*dx;
index_t nb_nodes = 50;
utils::ExponentialMovingAverage averager1(0.001, 5);
utils::ExponentialMovingAverage averager2(0.01, 5);
utils::ExponentialMovingAverage averager3(0.1, 5);
mesh::Mesh1DPtr themesh = mesh::axisymmetric_uniform_mesh1d(nb_nodes, radius, dx, height);
std::shared_ptr<SaturatedDiffusionTransportParameters> parameters =
std::make_shared<SaturatedDiffusionTransportParameters>(nb_nodes, 1e4*std::exp(9.95*0.20-29.08), 0.20);
parameters->diffusion_coefficient(0) = 2.2e-5; // cm^2/s
parameters->porosity(0) = 1.0;
EquilibriumState initial_state = get_initial_state(thedatabase, compo_species, 0.5);
SIABoundaryConditions bcs(nb_nodes);
bcs.list_initial_states.push_back(initial_state);
bcs.list_initial_states.push_back(get_blank_state(thedatabase));
bcs.bs_types[0] = BoundaryConditionType::FixedComposition;
//bcs.initial_states[0] = 1;
for (index_t node =0; node<additional_nodes; ++node)
bcs.initial_states[node] = 1;
database::Database dbhandler(thedatabase);
index_t id_ca = dbhandler.component_label_to_id("Ca[2+]");
index_t id_si = dbhandler.component_label_to_id("SiO(OH)3[-]");
index_t id_oh = dbhandler.component_label_to_id("HO[-]");
for (index_t mineral: thedatabase->range_mineral())
{
std::cout << thedatabase->labels_minerals[mineral] << std::endl;
}
std::cout << "Volume first cell " << themesh->get_volume_cell(additional_nodes) << std::endl;
//return EXIT_SUCCESS;
std::ofstream output_ca, output_aqca, output_si, output_poro, output_ph;
std::ofstream output_saca, output_sasi;
std::ofstream output_iter;
output_ca.open("out_ca.dat");
output_aqca.open("out_aqca.dat");
output_si.open("out_si.dat");
output_ph.open("out_ph.dat");
output_poro.open("out_poro.dat");
output_saca.open("out_sa_ca.dat");
output_sasi.open("out_sa_si.dat");
output_iter.open("nb_iter.dat");
std::vector<std::ofstream> mineral_out(thedatabase->nb_mineral);
for (int mineral: thedatabase->range_mineral())
{
mineral_out[mineral].open("out_mineral_"+std::to_string(mineral)+".dat");
}
SIASaturatedReactiveTransportSolver solver(themesh, thedatabase, parameters);
solver.apply_boundary_conditions(bcs);
solver.get_variables().scale_volume();
//solver.use_snia();
solver.use_sia(500);
solver.get_options().residual_tolerance = 1e-2;
solver.get_options().absolute_residual_tolerance = 1e-6;
solver.get_options().threshold_update_disable_speciation = 1e-20;
solver.get_options().threshold_update_transport = 1e-10;
solver.get_options().good_enough_transport_residual_tolerance = 1e-2;
solver.get_options().relative_update = 1e-6;
solver.get_options().use_consistent_velocity = true;
solver.get_options().use_previous_flux = false;
solver.get_options().disable_speciation_no_solid = true;
// Transport
solver.get_options().transport_options.residuals_tolerance = 1e-8;
solver.get_options().transport_options.maximum_iterations = 500;
solver.get_options().transport_options.step_tolerance = 1e-8;
solver.get_options().transport_options.threshold_stationary_point = 1e-5;
solver.get_options().transport_options.quasi_newton = 3;
solver.get_options().transport_options.sparse_solver = SparseSolver::GMRES;
//solver.get_options().transport_options.alpha = 0.0;
solver.get_options().transport_options.linesearch = dfpmsolver::ParabolicLinesearch::Strang;
// Speciation
solver.get_options().speciation_options.system_options.charge_keeper = 1;
solver.get_options().speciation_options.system_options.conservation_water = true;
solver.get_options().speciation_options.system_options.non_ideality = true;
solver.get_options().speciation_options.solver_options.fvectol = 1e-12;
solver.get_options().speciation_options.system_options.non_ideality_tolerance = 1e-14;
solver.get_options().speciation_options.solver_options.use_crashing = false;
solver.get_options().speciation_options.solver_options.use_scaling = false;
solver.get_options().speciation_options.solver_options.max_iter = 50;
solver.get_options().speciation_options.solver_options.maxstep = 20;
solver.get_options().speciation_options.solver_options.non_monotone_linesearch = true;
solver.get_options().speciation_options.solver_options.threshold_cycling_linesearch = 1e-8;
solver.get_options().speciation_options.solver_options.steptol = 1e-16;
solver.get_options().speciation_options.solver_options.threshold_stationary_point = 1e-4;
std::cout << "id SiO2,am : " << dbhandler.mineral_label_to_id("SiO2,am") << std::endl;
//thedatabase->stability_mineral(
// dbhandler.mineral_label_to_id("SiO2,am")) = database::MineralStabilityClass::cannot_dissolve;
SIASaturatedVariables& first_variables = solver.get_variables();
std::cout << first_variables.total_mobile_concentrations() << std::endl;
//double dt = 5e4;
double total=0 ;
uindex_t k = 0;
scalar_t dt = 4.0;
scalar_t dt_min = 4.0;
scalar_t dt_max = 4.0;
uindex_t max_iter;
while (total < 65*24*3600)
{
total += dt;
std::cout << " iterations : " << k << " - dt : " << dt << " - total : " << total/24/3600 << std::endl;
//start_again:
SIASaturatedReactiveTransportSolverReturnCode retcode = solver.solve_timestep(dt);
std::cout << " Return code : " << (int) retcode
<< " - Nb iter : " << solver.get_perfs().nb_iterations
<< std::endl
<< " - Nb iter transport : " << solver.get_perfs().nb_iterations_transport
<< std::endl
<< " - Nb iter chemistry : " << solver.get_perfs().nb_iterations_chemistry
<< std::endl;
if (retcode != SIASaturatedReactiveTransportSolverReturnCode::Success)
{
if (false and retcode == SIASaturatedReactiveTransportSolverReturnCode::MaxIterations)
{
max_iter += 1;
ERROR << "Maximum iterations";
}
else
{
std::cout << "Current residual : " << solver.get_perfs().current_residual << std::endl;
SIASaturatedVariables& variables = solver.get_variables();
output_ca << total << " " << total/(3600.0*24.0);
output_aqca << total << " " << total/(3600.0*24.0);
output_si << total << " " << total/(3600.0*24.0);
output_poro << total << " " << total/(3600.0*24.0);
output_ph << total << " " << total/(3600.0*24.0);
output_saca << total << " " << total/(3600.0*24.0);
output_sasi << total << " " << total/(3600.0*24.0);
for(index_t node: themesh->range_nodes())
{
output_ca << " " << variables.nodal_component_update_total_amount(node, id_ca);
output_aqca << " " << variables.total_mobile_concentration(node, id_ca);
output_si << " " << variables.nodal_component_update_total_amount(node, id_si);
output_ph << " " <<
14 + variables.log_gamma_component(node, id_oh)
+ variables.log_component_concentration(node, id_oh);
double volume = variables.volume_minerals(node);
double cell_volume = themesh->get_volume_cell(node);
output_poro << " " << (cell_volume - volume*1e6) / cell_volume;
output_saca << " " << variables.component_concentration_in_minerals(node, id_ca);
output_sasi << " " << variables.component_concentration_in_minerals(node, id_si);
}
output_ca << std::endl;
output_aqca << std::endl;
output_si << std::endl;
output_poro << std::endl;
output_ph << std::endl;
output_saca << std::endl;
output_sasi << std::endl;
output_ca.close();
output_aqca.close();
output_si.close();
output_poro.close();
output_saca.close();
output_sasi.close();
throw std::runtime_error("Failed to converge : "+std::to_string((int) retcode));
}
//dt = dt/2.0;
//solver.reset_flow();
//goto start_again;
}
//if (solver.get_perfs().nb_iterations <= 10) dt*=1.1;
SIASaturatedVariables& variables = solver.get_variables();
//std::cout << variables.total_mobile_concentrations() << std::endl;
// compute porosity
for (index_t node: themesh->range_nodes())
{
double volume = variables.volume_minerals(node);
double cell_volume = themesh->get_volume_cell(node);
specmicp_assert(volume < cell_volume);
double porosity = (cell_volume - volume*1e6) / cell_volume;
parameters->porosity(node) = porosity;
parameters->diffusion_coefficient(node) =
porosity>0.92?2.219e-5:1e4*std::exp(9.95*porosity-29.08);
}
std::cout << variables.component_concentration(1, id_ca) << std::endl;
parameters->diffusion_coefficient(0) = 2.219e-5; // parameters->diffusion_coefficient(1);
//if (solver.get_perfs().nb_iterations > 15) dt = 0.95*dt;
//else if (solver.get_perfs().nb_iterations < 2) dt = 1.05*dt;
output_iter << k
<< " " << solver.get_perfs().nb_iterations
<< " " << solver.get_perfs().current_residual
<< " " << averager1.add_point(solver.get_perfs().nb_iterations)
<< " " << averager2.add_point(solver.get_perfs().nb_iterations)
<< " " << averager3.add_point(solver.get_perfs().nb_iterations)
<< std::endl;
//std::cout << variables.speciation_variables() << std::endl;
if (k==1 or (k < 1000 and k% 100 == 0) or k % 1000 == 0)
{
output_ca << total << " " << total/(3600.0*24.0);
output_aqca << total << " " << total/(3600.0*24.0);
output_si << total << " " << total/(3600.0*24.0);
output_poro << total << " " << total/(3600.0*24.0);
output_ph << total << " " << total/(3600.0*24.0);
output_saca << total << " " << total/(3600.0*24.0);
output_sasi << total << " " << total/(3600.0*24.0);
for(index_t node: themesh->range_nodes())
{
output_ca << " " << variables.nodal_component_update_total_amount(node, id_ca);
output_aqca << " " << variables.total_mobile_concentration(node, id_ca);
output_si << " " << variables.nodal_component_update_total_amount(node, id_si);
output_ph << " " <<
14 + variables.log_gamma_component(node, id_oh)
+ variables.log_component_concentration(node, id_oh);
double volume = variables.volume_minerals(node);
double cell_volume = themesh->get_volume_cell(node);
output_poro << " " << (cell_volume - volume*1e6) / cell_volume;
output_saca << " " << variables.component_concentration_in_minerals(node, id_ca);
output_sasi << " " << variables.component_concentration_in_minerals(node, id_si);
}
output_ca << std::endl;
output_aqca << std::endl;
output_si << std::endl;
output_poro << std::endl;
output_ph << std::endl;
output_saca << std::endl;
output_sasi << std::endl;
for (index_t mineral: thedatabase->range_mineral())
{
mineral_out[mineral] << total << " " << total/(3600.0*24.0);
for(index_t node: themesh->range_nodes())
{
mineral_out[mineral] << " " << variables.mineral_amount(node, mineral);
}
mineral_out[mineral] << std::endl;
}
}
++k;
if ((averager3.current_value() - averager1.current_value()) > 10)
{
dt = std::max(dt*0.95,dt_min);
averager1.reset(averager3.current_value());
}
if ((averager2.current_value() - averager3.current_value()) > 10)
{
dt = std::min(dt*1.5, dt_max);
averager2.reset(averager3.current_value());
}
}
output_iter.close();
output_ca.close();
output_aqca.close();
output_si.close();
output_poro.close();
output_saca.close();
output_sasi.close();
for (index_t mineral: thedatabase->range_mineral())
{
mineral_out[mineral].close();
}
std::cout << "Max iterations : " << max_iter << std::endl;
return EXIT_SUCCESS;
}
int leaching_neutrality()
{
std::shared_ptr<database::DataContainer> thedatabase = set_database();
Vector compo_oxyde(4), compo_species(4);
compo_oxyde << 62.9, 20.6, 5.8, 3.1;
compo_from_oxyde(compo_oxyde, compo_species);
scalar_t radius = 3.5; //cm
scalar_t height = 8; //cm
index_t nb_nodes = 36;
database::Database dbhandler(thedatabase);
//mesh::Mesh1DPtr themesh = mesh::axisymmetric_uniform_mesh1d(nb_nodes, radius, height);
mesh::Mesh1DPtr themesh = mesh::uniform_mesh1d(nb_nodes, 0.1, 5);
//auto parameters = std::make_shared<SaturatedDiffusionTransportParameters>(nb_nodes, 1e-8, 0.2);
std::shared_ptr<SaturatedNeutralityDiffusionTransportParameters> parameters =
std::make_shared<SaturatedNeutralityDiffusionTransportParameters>(
nb_nodes, // nb_nodes,
thedatabase->nb_component, // nb_component,
thedatabase->nb_aqueous, // nb_aqueous,
2e-5, // the_diffusion_coefficient,
0.25, // porosity
0.5e-3 //reduction_factor
);
//parameters->diff_coeff_secondary(dbhandler.aqueous_label_to_id("H[+]")) = 9e-5;
//parameters->diff_coeff_component(dbhandler.component_label_to_id("HO[-]")) = 9e-5;
parameters->diff_coeff_secondary(dbhandler.aqueous_label_to_id("Si(OH)4")) = 1e-7;
parameters->diff_coeff_component(dbhandler.component_label_to_id("SiO(OH)3[-]")) = 1e-7;
parameters->porosity(0) = 1.0;
parameters->reduction_factor(0) = 1.0;
EquilibriumState initial_state = get_initial_state(thedatabase, compo_species, 0.5);
SIABoundaryConditions bcs(nb_nodes);
bcs.list_initial_states.push_back(initial_state);
bcs.list_initial_states.push_back(get_blank_state(thedatabase));
bcs.bs_types[0] = BoundaryConditionType::FixedComposition;
bcs.initial_states[0] = 1;
index_t id_ca = dbhandler.component_label_to_id("Ca[2+]");
index_t id_si = dbhandler.component_label_to_id("SiO(OH)3[-]");
std::ofstream output_ca, output_aqca, output_si, output_poro;
output_ca.open("out_ca.dat");
output_aqca.open("out_aqca.dat");
output_si.open("out_si.dat");
output_poro.open("out_poro.ca");
std::vector<std::ofstream> mineral_out(thedatabase->nb_mineral);
for (int mineral: thedatabase->range_mineral())
{
mineral_out[mineral].open("out_mineral_"+std::to_string(mineral)+".dat");
}
SIASaturatedReactiveTransportNeutralitySolver solver(themesh, thedatabase, parameters);
solver.apply_boundary_conditions(bcs);
solver.get_variables().scale_volume();
solver.use_sia(100);
solver.get_options().residual_tolerance = 1e-3;
solver.get_options().threshold_update_disable_speciation = 1e-16;
solver.get_options().threshold_update_transport = 1e-6;
// Transport
solver.get_options().transport_options.residuals_tolerance = 1e-5;
solver.get_options().transport_options.maximum_iterations = 100;
solver.get_options().transport_options.quasi_newton = 2;
solver.get_options().transport_options.maximum_step_length = 1000;
solver.get_options().transport_options.linesearch = dfpmsolver::ParabolicLinesearch::Strang;
solver.get_options().transport_options.max_iterations_at_max_length = 100;
// Speciation
solver.get_options().speciation_options.solver_options.use_crashing = false;
solver.get_options().speciation_options.solver_options.use_scaling = false;
solver.get_options().speciation_options.solver_options.max_iter = 100;
solver.get_options().speciation_options.solver_options.maxstep = 50;
std::cout << "id SiO2,am : " << dbhandler.mineral_label_to_id("SiO2,am") << std::endl;
//thedatabase->stability_mineral(
// dbhandler.mineral_label_to_id("SiO2,am")) = database::MineralStabilityClass::cannot_dissolve;
//double dt = 5e4;
double total=0 ;
uindex_t k = 1;
std::cout << "Basis :";
for (index_t component: thedatabase->range_component())
{
std:: cout << " " << thedatabase->labels_basis[component];
}
std::cout << std::endl;
std::cout << "Solid :";
for (index_t mineral: thedatabase->range_mineral())
{
std:: cout << " " << thedatabase->labels_minerals[mineral];
}
std::cout << std::endl;
while (total < 10*24*3600)
{
double dt = 3600/2;
total += dt;
std::cout << " iterations : " << k << " - dt : " << dt << std::endl;
SIASaturatedReactiveTransportSolverReturnCode retcode = solver.solve_timestep(dt);
if (retcode != SIASaturatedReactiveTransportSolverReturnCode::Success)
throw std::runtime_error("Failed to converge : "+std::to_string((int) retcode));
SIASaturatedVariables& variables = solver.get_variables();
//std::cout << variables.total_mobile_concentrations() << std::endl;
// compute porosity
for (index_t node: themesh->range_nodes())
{
double volume = variables.volume_minerals(node);
double cell_volume = themesh->get_volume_cell(node);
specmicp_assert(volume < cell_volume);
double porosity = (cell_volume - volume*1e6) / cell_volume;
parameters->porosity(node) = porosity;
parameters->reduction_factor(node) =
porosity>0.6?1e4*std::exp(9.95*0.6-29.08)/2e-5:1e4*std::exp(9.95*porosity-29.08)/2e-5;
//std::cout << parameters->diffusion_coefficient(node);
}
parameters->reduction_factor(0) = parameters->reduction_factor(1);
//std::cout << variables.speciation_variables() << std::endl;
if (k < 2 or k % 100 == 0)
{
output_ca << total << " " << total/(3600.0*24.0);
output_aqca << total << " " << total/(3600.0*24.0);
output_si << total << " " << total/(3600.0*24.0);
output_poro << total << " " << total/(3600.0*24.0);
for(index_t node: themesh->range_nodes())
{
output_ca << " " << variables.nodal_component_update_total_amount(node, id_ca);
output_aqca << " " << variables.total_mobile_concentration(node, id_ca);
output_si << " " << variables.nodal_component_update_total_amount(node, id_si);
double volume = variables.volume_minerals(node);
double cell_volume = themesh->get_volume_cell(node);
output_poro << " " << (cell_volume - volume*1e6) / cell_volume;
}
output_ca << std::endl;
output_aqca << std::endl;
output_si << std::endl;
output_poro << std::endl;
for (index_t mineral: thedatabase->range_mineral())
{
mineral_out[mineral] << total << " " << total/(3600.0*24.0);
for(index_t node: themesh->range_nodes())
{
mineral_out[mineral] << " " << variables.mineral_amount(node, mineral);
}
mineral_out[mineral] << std::endl;
}
}
++k;
}
output_ca.close();
output_aqca.close();
output_si.close();
output_poro.close();
for (index_t mineral: thedatabase->range_mineral())
{
mineral_out[mineral].close();
}
return EXIT_SUCCESS;
}
int main()
{
specmicp::stdlog::ReportLevel() = specmicp::logger::Warning;
specmicp::logger::ErrFile::stream() = &std::cerr;
return leaching();
//return leaching_neutrality();
}
Event Timeline
Log In to Comment