react_leaching.cpp
No OneTemporary
Actions

Subscribers

None

File Metadata

Created: Tue, Nov 5, 02:31

react_leaching.cpp
View Options


	#include "specmicp/adimensional/adimensional_system_solver.hpp"
	#include "specmicp/adimensional/adimensional_system_solution.hpp"
	#include "specmicp/problem_solver/dissolver.hpp"
	#include "specmicp/problem_solver/formulation.hpp"
	#include "specmicp/adimensional/adimensional_system_solution_extractor.hpp"

	#include "reactmicp/solver/staggers_base/upscaling_stagger_base.hpp"
	#include "reactmicp/solver/staggers_base/stagger_structs.hpp"
	#include "reactmicp/solver/reactive_transport_solver.hpp"
	#include "reactmicp/systems/saturated_react/variables.hpp"
	#include "reactmicp/systems/saturated_react/transport_stagger.hpp"
	#include "reactmicp/systems/saturated_react/equilibrium_stagger.hpp"
	#include "reactmicp/systems/saturated_react/init_variables.hpp"

	#include "dfpm/meshes/axisymmetric_uniform_mesh1d.hpp"
	#include "dfpm/meshes/axisymmetric_mesh1d.hpp"

	#include "utils/log.hpp"
	#include <iostream>
	#include <fstream>

	#include "utils/io/reactive_transport.hpp"

	using namespace specmicp;

	// Initial conditions
	// ==================

	specmicp::RawDatabasePtr leaching_database()
	{
	specmicp::database::Database thedatabase("../data/cemdata_specmicp.js");
	std::map<std::string, std::string> swapping ({
	{"H[+]","HO[-]"},
	{"Si(OH)4", "SiO(OH)3[-]"},
	{"Al[3+]","Al(OH)4[-]"}
	});
	thedatabase.remove_gas_phases();
	thedatabase.swap_components(swapping);
	return thedatabase.get_database();
	}


	AdimensionalSystemSolution initial_leaching_pb(
	RawDatabasePtr raw_data,
	scalar_t mult,
	units::UnitsSet the_units)
	{
	Formulation formulation;

	scalar_t m_c3s = mult*0.6;
	scalar_t m_c2s = mult*0.2;
	scalar_t m_c3a = mult*0.1;
	scalar_t m_gypsum = mult*0.1;
	scalar_t wc = 0.5;
	scalar_t m_water = wc1.0e-3(
	m_c3s(356.08+60.08)
	+ m_c2s(256.06+60.08)
	+ m_c3a(356.08+101.96)
	+ m_gypsum(56.08+80.06+218.02)
	);

	formulation.mass_solution = m_water;
	formulation.amount_minerals = {
	{"C3S", m_c3s},
	{"C2S", m_c2s},
	{"C3A", m_c3a},
	{"Gypsum", m_gypsum}
	};


	specmicp::Vector total_concentrations = specmicp::Dissolver(raw_data).dissolve(formulation);
	specmicp::AdimensionalSystemConstraints constraints(total_concentrations);
	constraints.charge_keeper = 1;
	constraints.set_saturated_system();

	specmicp::AdimensionalSystemSolverOptions options;
	options.solver_options.maxstep = 100.0;
	options.solver_options.max_iter = 200;
	options.solver_options.maxiter_maxstep = 100;
	options.solver_options.use_crashing = false;
	options.solver_options.use_scaling = true;
	options.solver_options.factor_descent_condition = -1.0;
	options.solver_options.factor_gradient_search_direction = 100;
	options.solver_options.projection_min_variable = 1e-9;
	options.solver_options.fvectol = 1e-10;
	options.solver_options.steptol = 1e-14;
	options.solver_options.non_monotone_linesearch = true;
	options.system_options.non_ideality_tolerance = 1e-10;
	options.units_set = the_units;

	specmicp::Vector x;

	specmicp::AdimensionalSystemSolver solver(raw_data, constraints, options);
	solver.initialise_variables(x, 0.8, -3.0);
	specmicp::micpsolver::MiCPPerformance perf = solver.solve(x, false);
	std::cout << "Initial return code : " << (int) perf.return_code << std::endl;
	if ((int) perf.return_code <= (int) specmicp::micpsolver::MiCPSolverReturnCode::NotConvergedYet)
	{
	throw std::runtime_error("Initial solution has not converged");
	}

	return solver.get_raw_solution(x);
	}


	AdimensionalSystemSolution initial_blank_leaching_pb(
	RawDatabasePtr raw_data,
	scalar_t mult,
	units::UnitsSet the_units)
	{
	Formulation formulation;

	scalar_t m_c3s = mult*0.6;
	scalar_t m_c2s = mult*0.2;
	scalar_t m_c3a = mult*0.1;
	scalar_t m_gypsum = mult*0.1;
	scalar_t wc = 0.8;
	scalar_t m_water = wc1.0e-3(
	m_c3s(356.08+60.08)
	+ m_c2s(256.06+60.08)
	+ m_c3a(356.08+101.96)
	+ m_gypsum(56.08+80.06+218.02)
	);

	formulation.mass_solution = m_water;
	formulation.amount_minerals = {
	{"SiO2,am", mult},
	};
	// formulation.concentration_aqueous = {
	// {"Ca[2+]", 1e-7},
	// {"Al(OH)4[-]", 1e-7},
	// {"SO4[2-]", 1e-7}
	// };
	formulation.extra_components_to_keep = {"Ca[2+]", "SO4[2-]", "Al(OH)4[-]"};

	specmicp::Vector total_concentrations = specmicp::Dissolver(raw_data).dissolve(formulation);
	specmicp::AdimensionalSystemConstraints constraints(total_concentrations);
	constraints.charge_keeper = 1;
	constraints.set_saturated_system();

	specmicp::AdimensionalSystemSolverOptions options;
	options.solver_options.maxstep = 10.0;
	options.solver_options.max_iter = 200;
	options.solver_options.maxiter_maxstep = 200;
	options.solver_options.use_crashing = false;
	options.solver_options.use_scaling = true;
	options.solver_options.factor_descent_condition = -1.0;
	options.solver_options.factor_gradient_search_direction = 100;
	options.solver_options.projection_min_variable = 1e-9;
	options.solver_options.fvectol = 1e-4;
	options.solver_options.steptol = 1e-14;
	options.solver_options.non_monotone_linesearch = true;
	options.system_options.non_ideality_tolerance = 1e-10;
	options.units_set = the_units;

	Vector x;

	AdimensionalSystemSolver solver(raw_data, constraints, options);
	solver.initialise_variables(x, 0.75, -4.0);
	micpsolver::MiCPPerformance perf = solver.solve(x, false);

	std::cout << "Initial return code : " << (int) perf.return_code << std::endl;
	if ((int) perf.return_code <= (int) specmicp::micpsolver::MiCPSolverReturnCode::NotConvergedYet)
	{
	throw std::runtime_error("Initial solution has not converged");
	}
	AdimensionalSystemSolution solution = solver.get_raw_solution(x);
	AdimensionalSystemSolutionModificator extractor(solution, raw_data, the_units);
	extractor.remove_solids();

	Vector totconc(raw_data->nb_component);
	totconc(0) = extractor.density_water()*extractor.total_aqueous_concentration(0);
	for (index_t component: raw_data->range_aqueous_component())
	{
	totconc(component) = 0.01extractor.density_water()extractor.total_aqueous_concentration(component);
	}

	constraints.total_concentrations = totconc;

	solver = AdimensionalSystemSolver(raw_data, constraints, options);
	perf = solver.solve(x, false);


	specmicp_assert((int) perf.return_code > (int) micpsolver::MiCPSolverReturnCode::NotConvergedYet);


	return solver.get_raw_solution(x);
	}

	// Upscaling
	// =========

	using namespace specmicp::reactmicp;
	using namespace specmicp::reactmicp::solver;
	using namespace specmicp::reactmicp::systems::satdiff;

	class LeachingUspcaling: public solver::UpscalingStaggerBase
	{
	public:

	LeachingUspcaling(RawDatabasePtr the_database,
	units::UnitsSet the_units):
	m_data(the_database),
	m_units_set(the_units)
	{

	}

	//! \brief Initialize the stagger at the beginning of the computation
	void initialize(VariablesBasePtr var)
	{
	SaturatedVariablesPtr true_var = cast_var_from_base(var);
	database::Database db_handler(m_data);
	m_dt = 1;
	m_id_cshj = db_handler.mineral_label_to_id("CSH,j");
	m_initial_sat_cshj = true_var->equilibrium_solution(1).main_variables(m_data->nb_component+m_id_cshj);
	//m_initial_sat_cshj = 0.36;
	std::cout << m_initial_sat_cshj << std::endl;
	true_var->upscaling_variables().setZero();
	for (index_t node=0; node<true_var->nb_nodes(); ++node)
	{
	upscaling_one_node(node, true_var);
	}


	}

	//! \brief Initialize the stagger at the beginning of an iteration
	void initialize_timestep(scalar_t dt, VariablesBasePtr var)
	{
	SaturatedVariablesPtr true_var = cast_var_from_base(var);
	m_dt = dt;

	for (index_t node=1; node<true_var->nb_nodes(); ++node)
	{
	upscaling_one_node(node, true_var);
	true_var->vel_porosity(node) = 0.0;
	}
	}

	//! \brief Solve the equation for the timestep
	StaggerReturnCode restart_timestep(VariablesBasePtr var)
	{
	SaturatedVariablesPtr true_var = cast_var_from_base(var);
	for (index_t node=1; node<true_var->nb_nodes(); ++node)
	{
	upscaling_one_node(node, true_var);
	}
	return StaggerReturnCode::ResidualMinimized;
	}

	void upscaling_one_node(index_t node, SaturatedVariablesPtr true_var)
	{
	AdimensionalSystemSolutionExtractor extractor(true_var->equilibrium_solution(node),
	m_data, m_units_set);


	scalar_t porosity = extractor.porosity();
	;//- true_var->equilibrium_solution(node).main_variables(m_data->nb_component)/m_initial_sat_cshj*0.05;
	//std::cout << "porosity : " << porosity
	// << " - saturation water" << true_var->equilibrium_solution(node).main_variables(0) << std::endl;
	true_var->vel_porosity(node) += (porosity - true_var->porosity(node))/m_dt;
	true_var->porosity(node) = porosity;

	true_var->diffusion_coefficient(node) = std::min(1e4std::exp(9.85porosity-29.08),1e-6);

	}

	private:
	RawDatabasePtr m_data;
	units::UnitsSet m_units_set;
	index_t m_id_cshj;
	scalar_t m_initial_sat_cshj;
	scalar_t m_dt;
	};

	//


	// Main
	// ====

	int main()
	{
	specmicp::logger::ErrFile::stream() = &std::cerr;
	specmicp::stdlog::ReportLevel() = specmicp::logger::Warning;

	RawDatabasePtr raw_data = leaching_database();

	units::UnitsSet the_units;
	the_units.length = units::LengthUnit::centimeter;

	AdimensionalSystemSolution initial = initial_leaching_pb(raw_data, 0.006, the_units);

	AdimensionalSystemSolutionModificator extractor(initial, raw_data, the_units);

	Database dbhandler(raw_data);
	index_t id_ca = dbhandler.component_label_to_id("Ca[2+]");
	index_t id_si = dbhandler.component_label_to_id("SiO(OH)3[-]");

	std::cout << extractor.total_solid_concentration(id_ca) << " - " << extractor.porosity() << std::endl;

	extractor.scale_total_concentration(id_ca, 0.015);
	std::cout << extractor.total_solid_concentration(id_ca) << " - " << extractor.porosity() << std::endl;

	std::cout << initial.main_variables << std::endl;

	AdimensionalSystemSolution blank = initial_blank_leaching_pb(raw_data, 0.005, the_units);

	std::cout << blank.main_variables << std::endl;


	scalar_t radius = 3.5; //cm
	scalar_t height = 8.0; //cm

	scalar_t dx = 0.0025;
	index_t additional_nodes = 1;

	radius = radius + additional_nodes*dx;

	index_t nb_nodes = 100;

	// mesh::Mesh1DPtr the_mesh = mesh::axisymmetric_uniform_mesh1d(nb_nodes, radius, dx, height);
	mesh::Mesh1DPtr the_mesh = mesh::uniform_axisymmetric_mesh1d(nb_nodes, radius, dx, height);

	std::vector<specmicp::AdimensionalSystemSolution> list_initial_composition = {initial, blank};
	std::vector<int> index_initial_state(nb_nodes, 0);
	for (int i=0; i< additional_nodes; ++i)
	index_initial_state[i] = 1;
	std::vector<specmicp::index_t> list_fixed_nodes = {0, };


	systems::satdiff::SaturatedVariablesPtr variables =
	systems::satdiff::init_variables(the_mesh, raw_data, the_units,
	list_fixed_nodes,
	list_initial_composition,
	index_initial_state);


	specmicp::AdimensionalSystemConstraints constraints;
	constraints.charge_keeper = 1;
	constraints.set_saturated_system();
	constraints.inert_volume_fraction = 0.0;

	AdimensionalSystemSolverOptions options;
	options.solver_options.maxstep = 100.0;
	options.solver_options.max_iter = 100;
	options.solver_options.maxiter_maxstep = 100;
	options.solver_options.use_crashing = false;
	options.solver_options.use_scaling = false;
	options.solver_options.factor_descent_condition = -1;
	options.solver_options.factor_gradient_search_direction = 100;
	options.solver_options.projection_min_variable = 1e-16;
	options.solver_options.fvectol = 1e-8;
	options.solver_options.steptol = 1e-12;
	options.solver_options.non_monotone_linesearch = true;
	options.system_options.non_ideality_tolerance = 1e-12;
	options.system_options.non_ideality_max_iter = 100;
	options.units_set = the_units;


	ChemistryStaggerPtr equilibrium_stagger =
	std::make_shared<reactmicp::systems::satdiff::EquilibriumStagger>(nb_nodes, constraints, options);

	std::shared_ptr<reactmicp::systems::satdiff::SaturatedTransportStagger> transport_stagger =
	std::make_shared<reactmicp::systems::satdiff::SaturatedTransportStagger>(variables, list_fixed_nodes);

	dfpmsolver::ParabolicDriverOptions& transport_options = transport_stagger->options_solver();

	transport_options.maximum_iterations = 450;
	transport_options.quasi_newton = 3;
	transport_options.residuals_tolerance = 1e-8;
	transport_options.step_tolerance = 1e-10;
	transport_options.alpha = 1;
	transport_options.linesearch = dfpmsolver::ParabolicLinesearch::Strang;
	transport_options.sparse_solver = SparseSolver::SparseLU;
	//transport_options.sparse_solver = SparseSolver::GMRES;
	transport_options.threshold_stationary_point = 1e-2;

	UpscalingStaggerPtr upscaling_stagger =
	std::make_shared<LeachingUspcaling>(raw_data, the_units);

	ReactiveTransportSolver solver(transport_stagger, equilibrium_stagger, upscaling_stagger);
	transport_stagger->initialize(variables);
	equilibrium_stagger->initialize(variables);
	upscaling_stagger->initialize(variables);

	solver.get_options().maximum_iterations = 500;
	solver.get_options().residuals_tolerance = 1e-2;

	solver.get_options().absolute_residuals_tolerance = 1e-18;
	solver.get_options().implicit_upscaling = true;

	std::ofstream output_iter("out_leaching_iter.dat");
	std::ofstream output_c_ca, output_s_ca, output_s_si;
	std::ofstream output_poro;
	std::ofstream output_loss_ca;
	output_c_ca.open("out_c_ca.dat");
	output_s_ca.open("out_s_ca.dat");
	output_s_si.open("out_s_si.dat");
	output_poro.open("out_poro.dat");
	output_loss_ca.open("out_loss_ca.dat");


	scalar_t initial_ca = 0;
	for (index_t node: the_mesh->range_nodes())
	{
	initial_ca += variables->solid_concentration(node, id_ca, variables->displacement())
	* the_mesh->get_volume_cell(node);
	}

	scalar_t total = 0.0;
	scalar_t dt=2.0;

	io::print_reactmicp_performance_long_header(&output_iter);

	int cnt =0;
	while (total < 25360024)
	{
	solver.solve_timestep(dt, variables);
	total += dt;
	io::print_reactmicp_performance_long(&output_iter, cnt, total, solver.get_perfs());

	std::cout << "time : " << total/3600/24 << std::endl
	<< "Iter : " << solver.get_perfs().nb_iterations << std::endl
	<< "Residuals : " << solver.get_perfs().residuals << std::endl;

	if (cnt % 5000 == 0)
	{
	scalar_t time = total/3600/24;
	scalar_t sqrt_time = std::sqrt(time);

	output_c_ca << time;
	output_s_ca << time;
	output_s_si << time;
	output_loss_ca << time << "\t" << sqrt_time;
	output_poro << time;
	scalar_t amount_ca = 0.0;
	for (index_t node: the_mesh->range_nodes())
	{
	amount_ca += variables->solid_concentration(node, id_ca, variables->displacement())
	* the_mesh->get_volume_cell(node);
	output_c_ca << "\t" << variables->aqueous_concentration(node, id_ca, variables->displacement());
	output_s_ca << "\t" << variables->solid_concentration( node, id_ca, variables->displacement());
	output_s_si << "\t" << variables->solid_concentration( node, id_si, variables->displacement());
	output_poro << "\t" << variables->porosity(node);
	}
	output_poro << std::endl;
	output_loss_ca << "\t" << (initial_ca - amount_ca)/1.75929e-2 << std::endl;
	output_c_ca << std::endl;
	output_s_ca << std::endl;
	output_s_si << std::endl;
	}
	++cnt;
	}

	io::print_reactmicp_timer(&std::cout, solver.get_timer());
	}

react_leaching.cppNo OneTemporaryActions

File Metadata

react_leaching.cppView Options

Event Timeline

react_leaching.cpp
No OneTemporary
Actions

react_leaching.cpp
View Options