Homec4science
Diffusion SpecMiCP / ReactMiCP (master) Browse

specmicp/examples/reactmicp/saturated_react/react_leaching.cpp
d829fc9ca5e8master

Last Change
Raw File

react_leaching.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. *
============================================================================= */
#include "reactmicp/reactmicp.hpp"
#include <iostream>
#include <fstream>
#include <Eigen/LU>
using namespace specmicp;
// Initial conditions
// ==================
specmicp::RawDatabasePtr leaching_database()
{
specmicp::database::Database thedatabase("../data/cemdata_specmicp.yaml");
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();
}
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);
}
AdimensionalSystemSolution initial_leaching_pb(
RawDatabasePtr raw_data,
scalar_t mult,
units::UnitsSet the_units)
{
Formulation formulation;
Vector oxyde_compo(4);
oxyde_compo << 62.9, 20.6, 5.8, 3.1;
Vector species_compo;
compo_from_oxyde(oxyde_compo, species_compo);
std::cout << species_compo << std::endl;
scalar_t m_c3s = mult*species_compo(0); //0.6;
scalar_t m_c2s = mult*species_compo(1); //0.2;
scalar_t m_c3a = mult*species_compo(2); //0.1;
scalar_t m_gypsum = mult*species_compo(3); //0.1;
scalar_t wc = 0.5;
scalar_t m_water = wc*1.0e-3*(
m_c3s*(3*56.08+60.08)
+ m_c2s*(2*56.06+60.08)
+ m_c3a*(3*56.08+101.96)
+ m_gypsum*(56.08+80.06+2*18.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.projection_min_variable = 1e-9;
options.solver_options.set_tolerance(1e-10, 1e-14);
options.solver_options.disable_crashing();
options.solver_options.enable_scaling();
options.solver_options.disable_descent_direction();
options.solver_options.enable_non_monotone_linesearch();
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;
Vector oxyde_compo(4);
oxyde_compo << 62.9, 20.6, 5.8, 3.1;
Vector species_compo;
compo_from_oxyde(oxyde_compo, species_compo);
std::cout << species_compo << std::endl;
scalar_t m_c3s = 0.6;
scalar_t m_c2s = 0.2;
scalar_t m_c3a = 0.1;
scalar_t m_gypsum = 0.1;
scalar_t wc = 0.8;
scalar_t m_water = wc*1.0e-3*(
m_c3s*(3*56.08+60.08)
+ m_c2s*(2*56.06+60.08)
+ m_c3a*(3*56.08+101.96)
+ m_gypsum*(56.08+80.06+2*18.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.01*extractor.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)
{
return initialize(var.get());
}
virtual void initialize(VariablesBase * const var) override
{
SaturatedVariables * const true_var = static_cast<SaturatedVariables * const>(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
virtual void initialize_timestep(scalar_t dt, VariablesBase * const var) override
{
SaturatedVariables * const true_var = static_cast<SaturatedVariables * const>(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
virtual StaggerReturnCode restart_timestep(VariablesBase * const var) override
{
SaturatedVariables* true_var = static_cast<SaturatedVariables*>(var);
for (index_t node=1; node<true_var->nb_nodes(); ++node)
{
upscaling_one_node(node, true_var);
}
return StaggerReturnCode::ResidualMinimized;
}
scalar_t diffusion_coefficient(scalar_t porosity)
{
const scalar_t factor = 1e4;
scalar_t poro = factor*std::exp(9.85*porosity-29.08);
// const scalar_t de_0 = 2.76e-12;
// const scalar_t porosity_r = 0.02;
// const scalar_t exponent = 3.32;
//if (porosity > 0.7)
// return factor*4.0*1e-10;
// else
// {
// const scalar_t ratio = ((porosity - porosity_r)/(0.25 - porosity_r));
// poro = factor*de_0*std::pow(ratio, exponent);
// }
return std::min(poro,factor*1e-10);
}
void upscaling_one_node(index_t node, SaturatedVariables * const 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) = diffusion_coefficient(porosity);
}
private:
RawDatabasePtr m_data;
units::UnitsSet m_units_set;
index_t m_id_cshj;
scalar_t m_initial_sat_cshj;
scalar_t m_dt;
};
//
void set_specmicp_options(AdimensionalSystemSolverOptions& options)
{
options.solver_options.maxstep = 10.0;
options.solver_options.max_iter = 100;
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;
options.solver_options.factor_gradient_search_direction = 100;
options.solver_options.projection_min_variable = 1e-16;
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-12;
options.system_options.non_ideality_max_iter = 100;
options.solver_options.threshold_cycling_linesearch = 1e-5;
}
void set_transport_options(dfpmsolver::ParabolicDriverOptions& transport_options)
{
transport_options.maximum_iterations = 450;
transport_options.quasi_newton = 3;
transport_options.residuals_tolerance = 1e-6;
transport_options.step_tolerance = 1e-14;
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;
transport_options.absolute_tolerance = 1e-16;
}
void set_reactive_solver_options(ReactiveTransportOptions& solver_options)
{
solver_options.maximum_iterations = 100;
solver_options.residuals_tolerance = 1e-2;
solver_options.good_enough_tolerance = 1.0;
solver_options.absolute_residuals_tolerance = 1e-16;
solver_options.implicit_upscaling = true;
}
mesh::Mesh1DPtr get_mesh(scalar_t dx, index_t nb_nodes, index_t additional_nodes)
{
mesh::UniformAxisymmetricMesh1DGeometry geometry;
geometry.dx = dx;
geometry.radius = 3.5 + additional_nodes*dx; //cm
geometry.height = 8.0;
geometry.nb_nodes = nb_nodes + additional_nodes;
return mesh::uniform_axisymmetric_mesh1d(geometry);
}
mesh::Mesh1DPtr get_mesh()
{
Vector radius(66);
const scalar_t height = 8.0;
radius <<
3.5005,
3.5000,
3.4995,
3.4990,
3.4980,
3.4970,
3.4960,
3.4950,
3.4925,
3.4900,
3.4875,
3.4850,
3.4825,
3.4800,
3.4775,
3.4750,
3.4725,
3.4700,
3.4675,
3.4650,
3.4625,
3.4600,
3.4575,
3.4550,
3.4525,
3.4500,
3.4475,
3.445,
3.4425,
3.440,
3.4375,
3.435,
3.4325,
3.430,
3.4275,
3.425,
3.4225,
3.420,
3.4175,
3.415,
3.4125,
3.410,
3.4075,
3.405,
3.4025,
3.400,
3.3975,
3.395,
3.3925,
3.390,
3.3875,
3.385,
3.3825,
3.38,
3.3775,
3.3750,
3.3725,
3.3700,
3.3675,
3.365,
3.3625,
3.36,
3.3575,
3.355,
3.325,
3.3
;
return mesh::axisymmetric_mesh1d(radius, height);
}
void output_mesh(mesh::Mesh1DPtr the_mesh)
{
std::ofstream output_mesh;
output_mesh.open("out_leaching_mesh.dat");
output_mesh << "Node\tPosition" << std::endl;
for (index_t node: the_mesh->range_nodes())
{
output_mesh << node << "\t" << the_mesh->get_position(node) << std::endl;
}
output_mesh.close();
}
void print_minerals_profile(
RawDatabasePtr the_database,
SaturatedVariablesPtr variables,
mesh::Mesh1DPtr the_mesh,
const std::string& filepath
)
{
std::ofstream ofile(filepath);
ofile << "Radius";
for (index_t mineral: the_database->range_mineral())
{
ofile << "\t" << the_database->get_label_mineral(mineral);
}
ofile << "\t" << "Porosity";
ofile << "\t" << "pH";
ofile << std::endl;
for (index_t node: the_mesh->range_nodes())
{
ofile << the_mesh->get_position(node);
AdimensionalSystemSolutionExtractor extractor(variables->equilibrium_solution(node), the_database, units::UnitsSet());
for (index_t mineral: the_database->range_mineral())
{
ofile << "\t" << extractor.volume_fraction_mineral(mineral);
}
ofile << "\t" << variables->porosity(node);
ofile << "\t" << extractor.pH();
ofile << std::endl;
}
ofile.close();
}
// Main
// ====
int main()
{
Eigen::initParallel();
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.015, the_units);
AdimensionalSystemSolutionModificator extractor(initial, raw_data, the_units);
index_t id_ca = raw_data->get_id_component("Ca[2+]");
index_t id_si = raw_data->get_id_component("SiO(OH)3[-]");
std::cout << extractor.total_solid_concentration(id_ca) << " - " << extractor.porosity() << std::endl;
extractor.scale_total_concentration(id_ca, 0.0145);
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 dx = 0.0005;
index_t nb_nodes = 150;
index_t additional_nodes = 1;
//mesh::Mesh1DPtr the_mesh = get_mesh(dx, nb_nodes, additional_nodes);
mesh::Mesh1DPtr the_mesh = get_mesh();
nb_nodes = the_mesh->nb_nodes();
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;
set_specmicp_options(options);
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();
set_transport_options(transport_options);
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);
set_reactive_solver_options(solver.get_options());
io::OutFile 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 dt=2.0;
io::print_reactmicp_performance_long_header(output_iter);
Timestepper timestepper(1.0, 500.0, 25*24*3600, 2);
int cnt =0;
while (timestepper.get_total() < timestepper.get_total_target())
{
Timer step_timer;
step_timer.start();
reactmicp::solver::ReactiveTransportReturnCode retcode = solver.solve_timestep(dt, variables);
step_timer.stop();
io::print_reactmicp_performance_long(output_iter, cnt, timestepper.get_total() , solver.get_perfs());
dt = timestepper.next_timestep(dt, retcode, solver.get_perfs().nb_iterations);
if (retcode <= reactmicp::solver::ReactiveTransportReturnCode::NotConvergedYet)
{
dt = 1.0;
variables->reset_main_variables();
retcode = solver.solve_timestep(dt, variables);
dt = timestepper.next_timestep(dt, retcode, solver.get_perfs().nb_iterations);
io::print_reactmicp_performance_long(output_iter, cnt, timestepper.get_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 % 100 == 0)
{
scalar_t time = timestepper.get_total()/3600.0/24.0;
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;
}
print_minerals_profile(raw_data, variables, the_mesh, "out_leaching_minerals_profile_25d.dat");
timestepper.set_total_target(90*24*3600);
while (timestepper.get_total() < timestepper.get_total_target())
{
Timer step_timer;
step_timer.start();
reactmicp::solver::ReactiveTransportReturnCode retcode = solver.solve_timestep(dt, variables);
step_timer.stop();
io::print_reactmicp_performance_long(output_iter, cnt, timestepper.get_total() , solver.get_perfs());
dt = timestepper.next_timestep(dt, retcode, solver.get_perfs().nb_iterations);
if (retcode <= reactmicp::solver::ReactiveTransportReturnCode::NotConvergedYet)
{
dt = 1.0;
variables->reset_main_variables();
retcode = solver.solve_timestep(dt, variables);
dt = timestepper.next_timestep(dt, retcode, solver.get_perfs().nb_iterations);
io::print_reactmicp_performance_long(output_iter, cnt, timestepper.get_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 % 100 == 0)
{
scalar_t time = timestepper.get_total()/3600.0/24.0;
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(output_iter, solver.get_timer());
}
c4science · Help