diff --git a/examples/reactmicp/saturated_react/carbonation.cpp b/examples/reactmicp/saturated_react/carbonation.cpp index 45ca71b..c7ce37c 100644 --- a/examples/reactmicp/saturated_react/carbonation.cpp +++ b/examples/reactmicp/saturated_react/carbonation.cpp @@ -1,692 +1,693 @@ // ====================================================== // Leaching of a cement paste in CO2-saturated water // ======================================================= // This file can be used to learn how to use ReactMiCP // Include ReactMiCP #include "reactmicp.hpp" #include "specmicp/adimensional/adimensional_system_solution_saver.hpp" #include "reactmicp/systems/saturated_react/diffusive_upscaling_stagger.hpp" // Other includes that would be needed later : #include #include "physics/laws.hpp" #include "utils/timer.hpp" #include #include #include // 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); // ============== // 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)); // 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; 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_c_hco3 << std::endl; out_c_ph << std::endl; out_s_ca << std::endl; out_calcite << std::endl; out_poro << 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; 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.js"); // 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 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 list_constraints = {cement_constraints, bc_constraints}; // index_constraints links a node to the set of constraints that will be used std::vector index_constraints(nb_nodes, 0); index_constraints[0] = 1; // Boundary conditions // This is the list of initial composition std::vector list_initial_composition = {cement_solution, bc_solution}; // This list links a node to its initial conditions std::vector 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); // Staggers // --------- // This section initialize the staggers // The equilibrium stagger // - - - - - - - - - - - - std::shared_ptr equilibrium_stagger = std::make_shared(list_constraints, index_constraints, specmicp_options); // The transport stagger // - - - - - - - - - - - std::shared_ptr transport_stagger = std::make_shared(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(upscaling_law.get_law(),the_database, units_set); // Solver // ------ // This section initialize the reactive transport solver ReactiveTransportSolver solver(transport_stagger, equilibrium_stagger, upscaling_stagger); set_reactive_solver_options(solver.get_options()); transport_stagger->initialize(variables); equilibrium_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("carbo", 3600); simul_info.output_prefix = "out_carbo_"; // 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(1.0*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"); 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, simul_info.complete_filepath("profiles_end", "dat")); + 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 = 50; //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 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.947e-2; // poro 0.47 //species_compo *= 1.08e-2; // poro 0.4 //species_compo *= 1.25e-2; // poro 0.3 species_compo *= 1.518e-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(); 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 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 solver(publi_stochio); oxyde = solver.solve(publi); Vector for_cemdata = cemdata_stochio*oxyde; return for_cemdata; } mesh::Mesh1DPtr get_mesh(scalar_t dx) { const scalar_t total_length = 0.75/2.0 + dx; //cm const scalar_t height = 20.0; index_t nb_nodes = total_length/dx + 1; return mesh::uniform_axisymmetric_mesh1d(nb_nodes, total_length, dx, height); } 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); } diff --git a/examples/reactmicp/saturated_react/carbonationfe.cpp b/examples/reactmicp/saturated_react/carbonationfe.cpp index 6229e12..73b0a7b 100644 --- a/examples/reactmicp/saturated_react/carbonationfe.cpp +++ b/examples/reactmicp/saturated_react/carbonationfe.cpp @@ -1,762 +1,763 @@ // ============================================ // Saturated carbonation of a cement past // ============================================= // This file can be used to learn how to use ReactMiCP // Include ReactMiCP #include "reactmicp.hpp" // Physical laws #include "physics/laws.hpp" // A simple timer #include "utils/timer.hpp" #include #include #include // 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); void compo_oxyde_to_mole(Vector& compo_oxyde); // 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); // Upscaling // ========= // Upscaling is the stagger in charge of finding the porosity and transport properties // This is a class that is dependant on the problem to solve and should be defined by the user // The following class is an example on how to define such a class // The Upscaling class should inherit from UpscalingStaggerBase class CarboUspcaling: public solver::UpscalingStaggerBase { public: // Initiation // This is not call automatically so can be adjusted by the user CarboUspcaling(RawDatabasePtr the_database, units::UnitsSet the_units): m_data(the_database), m_units_set(the_units), m_dt(HUGE_VAL) { } // Initialize the stagger at the beginning of the computation // The porosity and transport properties should be initialized here void initialize(VariablesBasePtr var) { SaturatedVariablesPtr true_var = cast_var_from_base(var); true_var->upscaling_variables().setZero(); for (index_t node=0; nodenb_nodes(); ++node) { upscaling_one_node(node, true_var); } } // Initialize the stagger at the beginning of a timestep // Typically do nothing but erase the "dot" variable (velocity such as the vel_porosity) void initialize_timestep(scalar_t dt, VariablesBasePtr var) { SaturatedVariablesPtr true_var = cast_var_from_base(var); m_dt = dt; for (index_t node=1; nodenb_nodes(); ++node) { upscaling_one_node(node, true_var); true_var->vel_porosity(node) = 0.0; } } //! This is the main function called during a timestep StaggerReturnCode restart_timestep(VariablesBasePtr var) { SaturatedVariablesPtr true_var = cast_var_from_base(var); for (index_t node=1; nodenb_nodes(); ++node) { upscaling_one_node(node, true_var); } return StaggerReturnCode::ResidualMinimized; // This is the value that should be returned if everything is ok } // A function that return the diffusion coefficient when porosity is equal to 'porosity' scalar_t coeff_diffusion_law(scalar_t porosity) { static constexpr scalar_t de_0 = 2.726e-12; //1e-11; static constexpr scalar_t factor = 1e4; static constexpr scalar_t porosity_r = 0.02; const scalar_t ratio = ((porosity - porosity_r)/(0.25 - porosity_r)); const scalar_t d_e = factor*de_0*std::pow(ratio, 3.32); return std::min(d_e, 1e4*1e-10); } // Compute the upscaling for 'node' void upscaling_one_node(index_t node, SaturatedVariablesPtr true_var) { // AdimensionalSystemSolutionExtractor is the class to use to // extract information from a SpecMiCP solution // To obtain correct information the correct units must be used AdimensionalSystemSolutionExtractor extractor(true_var->equilibrium_solution(node), m_data, m_units_set); // We can obtain the porosity very easily : scalar_t porosity = extractor.porosity(); true_var->vel_porosity(node) += (porosity - true_var->porosity(node))/m_dt; true_var->porosity(node) = porosity; true_var->diffusion_coefficient(node) = coeff_diffusion_law(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; }; // ============== // 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)); // 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; 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_c_hco3 << std::endl; out_c_ph << std::endl; out_s_ca << std::endl; out_calcite << std::endl; out_poro << 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; 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.js"); // 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); // Obtain the initial condition AdimensionalSystemSolution cement_solution = solve_equilibrium( the_database, cement_constraints, specmicp_options); // 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 ); // The constraints // -------------- // Boundary condition // 'list_fixed_nodes' lists the node with fixed composition std::vector 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 list_constraints = {cement_constraints, bc_constraints}; // index_constraints links a node to the set of constraints that will be used std::vector index_constraints(nb_nodes, 0); index_constraints[0] = 1; // Boundary conditions // This is the list of initial composition std::vector list_initial_composition = {cement_solution, bc_solution}; // This list links a node to its initial conditions std::vector 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); // Staggers // --------- // This section initialize the staggers std::shared_ptr equilibrium_stagger = std::make_shared(list_constraints, index_constraints, specmicp_options); std::shared_ptr transport_stagger = std::make_shared(variables, list_fixed_nodes); set_transport_options(transport_stagger->options_solver()); UpscalingStaggerPtr upscaling_stagger = std::make_shared(the_database, units_set); // Solver // ------ // This section initialize the reactive transport solver ReactiveTransportSolver solver(transport_stagger, equilibrium_stagger, upscaling_stagger); set_reactive_solver_options(solver.get_options()); transport_stagger->initialize(variables); equilibrium_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("carbo", 3600); simul_info.output_prefix = "out_carbo_"; // 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(30*24*3600, cast_var_to_base(variables)); // output information at the end of the timestep - io::print_minerals_profile(the_database, variables, the_mesh, simul_info.complete_filepath("profiles_end", "dat")); + 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 = 50.0; options.solver_options.threshold_cycling_linesearch = 1e-6; options.solver_options.set_tolerance(1e-8, 1e-16); 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.system_options.scaling_electron = 1e4; 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 = 50; //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 swapping({ {"H[+]","HO[-]"}, {"Si(OH)4", "SiO(OH)3[-]"}, {"Al[3+]","Al(OH)4[-]"}, {"Fe[2+]", "Fe(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(5); oxyde_compo << 66.98, 21.66, 2.78, 2.96, 4.41; compo_oxyde_to_mole(oxyde_compo); scalar_t mult = laws::density_water(units::celsius(25.0), units_set.length, units_set.mass); //oxyde_compo *= 0.01098; //= 1e-2; //poro 0.4 oxyde_compo *= 0.00967; //= 1e-2; //poro 0.47 Formulation formulation; formulation.mass_solution = 1.0; formulation.amount_minerals = { {"CaO(ox)", oxyde_compo(0)}, {"SiO2(ox)", oxyde_compo(1)}, {"Al2O3(ox)", oxyde_compo(2)}, {"SO3(ox)", oxyde_compo(3)}, {"Fe2O3(ox)", oxyde_compo(4)}, {"Calcite", oxyde_compo(0)*1e-4} }; 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,am", "C3AH6", "C4AH13", "C2AH8", "CAH10", "Monosulfoaluminate", "Tricarboaluminate", "Monocarboaluminate", "Hemicarboaluminate", //"Straetlingite", "Gypsum", "Ettringite", //"Thaumasite", "C3FH6", "C4FH13", "C2FH8", "Fe(OH)3(mic)", "Fe-Ettringite", "Fe-Monosulfate", "Fe-Monocarbonate", "Fe-Hemicarbonate", //"Fe-Straetlingite", }; 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(); 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_oxyde_to_mole(Vector& compo_oxyde) { constexpr double M_CaO = 56.08; constexpr double M_SiO2 = 60.09; constexpr double M_Al2O3 = 101.96; constexpr double M_SO3 = 80.06; constexpr double M_FE2O3 = 159.7; Vector molar_mass(5); molar_mass << 1.0/M_CaO, 1.0/M_SiO2, 1.0/M_Al2O3, 1.0/M_SO3, 1.0/M_FE2O3; compo_oxyde = compo_oxyde.cwiseProduct(molar_mass); } mesh::Mesh1DPtr get_mesh(scalar_t dx) { const scalar_t total_length = 0.75/2.0 + dx; //cm const scalar_t height = 20.0; index_t nb_nodes = total_length/dx + 1; return mesh::uniform_axisymmetric_mesh1d(nb_nodes, total_length, dx, height); } mesh::Mesh1DPtr get_mesh() { Vector radius(71); const scalar_t height = 20; // radius << // 3.751, // 3.750, // 3.749, // 3.748, // 3.747, // 3.745, // 3.740, // 3.735, // 3.730, // 3.720, // 3.710, // 3.700, // 3.675, // 3.650, // 3.625, // 3.600, // 3.575, // 3.550, // 3.525, // 3.500, // 3.475, // 3.450, // 3.425, // 3.400, // 3.375, // 3.350, // 3.325, // 3.300, // 3.250, // 3.200, // 3.150, // 3.100, // 3.050, // 3.000, // 2.975, // 2.950, // 2.925, // 2.900, // 2.875, // 2.850, // 2.825, // 2.800, // 2.775, // 2.750, // 2.700, // 2.650; 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); } diff --git a/examples/reactmicp/saturated_react/leaching_kinetic.cpp b/examples/reactmicp/saturated_react/leaching_kinetic.cpp index 130abe3..2c749b9 100644 --- a/examples/reactmicp/saturated_react/leaching_kinetic.cpp +++ b/examples/reactmicp/saturated_react/leaching_kinetic.cpp @@ -1,843 +1,844 @@ // ============================================ // 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.hpp" #include "specmicp/adimensional/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 #include "physics/laws.hpp" #include "utils/timer.hpp" #include #include #include // 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; // 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 void initialize_timestep(scalar_t dt, KineticSaturatedVariablesBasePtr kinetic_variables) override; //! \brief Restart a timestep void restart_timestep( index_t node, const AdimensionalSystemSolution& eq_solution, KineticSaturatedVariablesBasePtr kinetic_variables) override; //! \brief Return the volume fraction of kinetic (or inert) phases 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 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(); 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.length)); var->volume_fractions(0) = 0.0; var->amount_minerals(0) = 0.0; return std::static_pointer_cast(var); } void C3SDissolutionModel::initialize_timestep(scalar_t dt, KineticSaturatedVariablesBasePtr kinetic_variables) { m_dt = dt; C3SKineticVariablesPtr true_kin_var = std::static_pointer_cast(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(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.length); } 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(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(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(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.js"); // 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 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 list_constraints = {cement_constraints, bc_constraints}; // index_constraints links a node to the set of constraints that will be used std::vector index_constraints(nb_nodes, 0); index_constraints[0] = 1; // Boundary conditions // This is the list of initial composition std::vector list_initial_composition = {cement_solution, bc_solution}; // This list links a node to its initial conditions std::vector 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 kinetic_model = std::make_shared(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 kinetic_stagger = std::make_shared(kinetic_model, list_constraints, index_constraints, specmicp_options); // The transport stagger // - - - - - - - - - - - std::shared_ptr transport_stagger = std::make_shared(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(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"); 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, simul_info.complete_filepath("profiles_end", "dat")); + 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 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 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 solver(publi_stochio); oxyde = solver.solve(publi); Vector for_cemdata = cemdata_stochio*oxyde; return for_cemdata; } mesh::Mesh1DPtr get_mesh(scalar_t dx) { const scalar_t total_length = 0.75/2.0 + dx; //cm const scalar_t height = 20.0; index_t nb_nodes = total_length/dx + 1; return mesh::uniform_axisymmetric_mesh1d(nb_nodes, total_length, dx, height); } 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); } diff --git a/examples/reactmicp/saturated_react/momas_benchmark.cpp b/examples/reactmicp/saturated_react/momas_benchmark.cpp index 27a6fa0..366542e 100644 --- a/examples/reactmicp/saturated_react/momas_benchmark.cpp +++ b/examples/reactmicp/saturated_react/momas_benchmark.cpp @@ -1,565 +1,556 @@ #include "reactmicp.hpp" #include "database/aqueous_selector.hpp" #include #include #include #include // parameters #define DX 0.01 #define IS_ADVECTIVE true #define MIN_DT 0.01 #define MAX_DT 100.0 #define RESTART_DT MIN_DT using namespace specmicp; RawDatabasePtr get_momas_db() { database::Database thedatabase("../data/momas_benchmark.js"); return thedatabase.get_database(); } void prepare_db_for_easy_case(RawDatabasePtr raw_data) { database::Database db_handler(raw_data); db_handler.remove_components({"X5",}); if (raw_data->nb_component() != 6) { throw std::runtime_error("Not enough or Too many component"); } db_handler.remove_solid_phases(); if (not raw_data->is_valid()) { throw std::runtime_error("Database is not valid."); } //std::cout << selector.aqueous_label_to_id("C6") << " - " << selector.aqueous_label_to_id("C7") << std::endl; database::AqueousSelector selector(raw_data); selector.remove_aqueous({selector.aqueous_label_to_id("C6"), selector.aqueous_label_to_id("C7")}); std::cout << raw_data->aqueous.get_nu_matrix() << std::endl; if (raw_data->nb_aqueous() != 5) { throw std::runtime_error("Not enough or Too many aqueous."); } if (not raw_data->is_valid()) { throw std::runtime_error("Database is not valid."); } } AdimensionalSystemSolution easy_material_a( RawDatabasePtr raw_data, const AdimensionalSystemConstraints& constraints, // const AdimensionalSystemSolution& initial_solution, const AdimensionalSystemSolverOptions& options ) { specmicp::Vector variables ; specmicp::AdimensionalSystemSolver solver(raw_data, constraints, options); specmicp::micpsolver::MiCPPerformance current_perf = solver.solve(variables, true); if (current_perf.return_code != micpsolver::MiCPSolverReturnCode::ResidualMinimized) throw std::runtime_error("Failed to solve initial solution for material A."); return solver.get_raw_solution(variables); } AdimensionalSystemSolution easy_material_b( RawDatabasePtr raw_data, const AdimensionalSystemConstraints& constraints, const specmicp::AdimensionalSystemSolverOptions& options ) { specmicp::Vector variables; specmicp::AdimensionalSystemSolver solver(raw_data, constraints, options); specmicp::micpsolver::MiCPPerformance current_perf = solver.solve(variables, true); if (current_perf.return_code != micpsolver::MiCPSolverReturnCode::ResidualMinimized) throw std::runtime_error("Failed to solve initial solution for material B."); return solver.get_raw_solution(variables); } AdimensionalSystemSolution easy_bc( RawDatabasePtr raw_data, const AdimensionalSystemConstraints& constraints, const specmicp::AdimensionalSystemSolverOptions& options ) { specmicp::Vector variables; specmicp::AdimensionalSystemSolver solver(raw_data, constraints, options); solver.initialise_variables(variables, 1.0, -4.0); solver.run_pcfm(variables); specmicp::micpsolver::MiCPPerformance current_perf = solver.solve(variables, false); if (current_perf.return_code != micpsolver::MiCPSolverReturnCode::ResidualMinimized) throw std::runtime_error("Failed to solve initial solution for BC."); return solver.get_raw_solution(variables); } using namespace specmicp::reactmicp; using namespace specmicp::reactmicp::solver; using namespace specmicp::reactmicp::systems::satdiff; class MoMasUspcaling: public solver::UpscalingStaggerBase { public: MoMasUspcaling( RawDatabasePtr the_database, units::UnitsSet the_units, bool advective ): m_data(the_database), m_units_set(the_units), m_advective(advective) { } //! \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); true_var->upscaling_variables().setZero(); for (index_t node=0; nodenb_nodes(); ++node) { // if material b const scalar_t coord = true_var->get_mesh()->get_position(node); if (coord < 1.1 and coord >= 1.0) { true_var->porosity(node) = 0.5; true_var->permeability(node) = 1e-5; if (m_advective) true_var->diffusion_coefficient(node) = 3.3e-4; else true_var->diffusion_coefficient(node) = 330e-3; true_var->fluid_velocity(node) = 5.5e-3; } else // material a { true_var->porosity(node) = 0.25; true_var->permeability(node) = 1e-2; if (m_advective) true_var->diffusion_coefficient(node) = 5.5e-5; else true_var->diffusion_coefficient(node) = 55e-3; true_var->fluid_velocity(node) = 5.5e-3; } std::cout << "node : " << node << " - porosity : " << true_var->porosity(node) << std::endl; } } //! \brief Initialize the stagger at the beginning of an iteration void initialize_timestep(scalar_t dt, VariablesBasePtr var) { } //! \brief Solve the equation for the timestep StaggerReturnCode restart_timestep(VariablesBasePtr var) { return StaggerReturnCode::ResidualMinimized; } private: RawDatabasePtr m_data; units::UnitsSet m_units_set; bool m_advective; }; AdimensionalSystemConstraints get_constraints_easy_a(RawDatabasePtr& raw_data) { Vector total_concentrations(Vector::Zero(raw_data->nb_component())); total_concentrations(raw_data->get_id_component("X2")) = -2.0; total_concentrations(raw_data->get_id_component("X4")) = 2.0; AdimensionalSystemConstraints constraints_a(0.25*total_concentrations); constraints_a.disable_conservation_water(); constraints_a.surface_model.model_type = SurfaceEquationType::Equilibrium; constraints_a.surface_model.concentration = 0.25*1.0; constraints_a.inert_volume_fraction = 0.75; return constraints_a; } AdimensionalSystemConstraints get_constraints_easy_b(RawDatabasePtr& raw_data) { Vector total_concentrations(Vector::Zero(raw_data->nb_component())); total_concentrations(raw_data->get_id_component("X2")) = -2.0; total_concentrations(raw_data->get_id_component("X4")) = 2.0; AdimensionalSystemConstraints constraints_b(0.5*total_concentrations); constraints_b.disable_conservation_water(); constraints_b.surface_model.model_type = SurfaceEquationType::Equilibrium; constraints_b.surface_model.concentration = 0.5*10.0; constraints_b.inert_volume_fraction = 0.5; return constraints_b; } AdimensionalSystemConstraints get_constraints_easy_bc(RawDatabasePtr& raw_data) { Vector total_concentrations(Vector::Zero(raw_data->nb_component())); total_concentrations(raw_data->get_id_component("X1")) = 0.3; total_concentrations(raw_data->get_id_component("X2")) = 0.3; total_concentrations(raw_data->get_id_component("X3")) = 0.3; AdimensionalSystemConstraints constraints_bc(total_concentrations); constraints_bc.disable_conservation_water(); constraints_bc.surface_model.model_type = SurfaceEquationType::NoEquation; constraints_bc.surface_model.concentration = 0.0; constraints_bc.inert_volume_fraction = 0.0; return constraints_bc; } AdimensionalSystemSolverOptions get_specmicp_options() { AdimensionalSystemSolverOptions options; options.solver_options.maxstep = 1.0; options.solver_options.max_iter = 200; options.solver_options.maxiter_maxstep = 200; options.solver_options.use_crashing = false; options.solver_options.use_scaling = false; options.solver_options.non_monotone_linesearch = false; options.solver_options.factor_descent_condition = -1; 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.threshold_cycling_linesearch = 1e-8; options.system_options.non_ideality_tolerance = 1e-10; options.system_options.non_ideality = false; options.system_options.restart_concentration = -2; options.use_pcfm = true; return options; } mesh::Mesh1DPtr get_mesh(scalar_t dx) { scalar_t tot_length = 2.1+dx; index_t nb_nodes = tot_length/dx +2; return mesh::uniform_mesh1d(nb_nodes, dx, 1.0); } 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-12; transport_options.absolute_tolerance = 1e-6; 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-2; } void set_reactive_solver_options(reactmicp::solver::ReactiveTransportOptions& reactive_options) { reactive_options.maximum_iterations = 50; //500; reactive_options.residuals_tolerance = 1e-2; reactive_options.good_enough_tolerance = 1.0; reactive_options.absolute_residuals_tolerance = 1e-6; reactive_options.implicit_upscaling = false; } int main() { Eigen::initParallel(); Timer global_timer; global_timer.start(); - std::ofstream output_gen("out_momas_out.txt"); + io::OutFile output_gen("out_momas_out.txt"); - io::print_reactmicp_header(&output_gen); + io::print_reactmicp_header(output_gen); - specmicp::logger::ErrFile::stream() = &output_gen; //&std::cerr; + specmicp::logger::ErrFile::stream() = &std::cerr; specmicp::stdlog::ReportLevel() = specmicp::logger::Warning; RawDatabasePtr database = get_momas_db(); prepare_db_for_easy_case(database); units::UnitsSet the_units = units::UnitsSet(); the_units.length = units::LengthUnit::decimeter; // so that rho_w ~ 1 scalar_t dx=DX; mesh::Mesh1DPtr the_mesh = get_mesh(dx); AdimensionalSystemConstraints constraints_a = get_constraints_easy_a(database); AdimensionalSystemConstraints constraints_b = get_constraints_easy_b(database); AdimensionalSystemConstraints constraints_bc = get_constraints_easy_bc(database); AdimensionalSystemSolverOptions options = get_specmicp_options(); options.units_set = the_units; AdimensionalSystemSolution mat_bc = easy_bc(database, constraints_bc, options); AdimensionalSystemSolution mat_a = easy_material_a(database, constraints_a, options); AdimensionalSystemSolution mat_b = easy_material_b(database, constraints_b, options); options.solver_options.maxstep = 1.0; options.solver_options.use_scaling = true; options.solver_options.non_monotone_linesearch = true; index_t nb_nodes = the_mesh->nb_nodes(); std::vector list_constraints = {constraints_a, constraints_b, constraints_bc}; std::vector list_initial_composition = {mat_a, mat_b, mat_bc}; std::vector index_initial_state(nb_nodes, 0); std::vector index_constraints(nb_nodes, 0); index_t start_mat_b = 1.0/dx; index_t end_mat_b = 1.1/dx; for (index_t node=start_mat_b; node list_fixed_nodes = {0, }; std::map list_slaves_nodes = {{nb_nodes-1, nb_nodes-2}}; std::vector list_immobile_components = {0, 1}; systems::satdiff::SaturatedVariablesPtr variables = systems::satdiff::init_variables(the_mesh, database, the_units, list_fixed_nodes, list_initial_composition, index_initial_state); ChemistryStaggerPtr equilibrium_stagger = std::make_shared(list_constraints, index_constraints, options); std::shared_ptr transport_stagger = std::make_shared( variables, list_fixed_nodes, list_slaves_nodes, list_immobile_components); set_transport_options(transport_stagger->options_solver()); const bool is_advective = IS_ADVECTIVE; UpscalingStaggerPtr upscaling_stagger = std::make_shared(database, the_units, is_advective); 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()); // mesh output - std::ofstream output_mesh; - output_mesh.open("out_momas_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)-dx/2.0 << std::endl; - } - output_mesh.close(); - + io::print_mesh("out_momas_mesh.dat", the_mesh, the_units.length); //std::cout << variables->displacement() << std::endl; - std::ofstream output_iter("out_momas_iter.dat"); + io::OutFile output_iter("out_momas_iter.dat"); std::ofstream output_x1, output_x2, output_x3, output_x4; std::ofstream output_sx1, output_sx2, output_sx3, output_sx4; std::ofstream output_c1, output_c2, output_s; output_x1.open("out_momas_x1.dat"); output_x2.open("out_momas_x2.dat"); output_x3.open("out_momas_x3.dat"); output_x4.open("out_momas_x4.dat"); output_sx1.open("out_momas_sx1.dat"); output_sx2.open("out_momas_sx2.dat"); output_sx3.open("out_momas_sx3.dat"); output_sx4.open("out_momas_sx4.dat"); output_c1.open("out_momas_c1.dat"); output_c2.open("out_momas_c2.dat"); output_s.open("out_momas_s.dat"); index_t cnt = 0; scalar_t dt=RESTART_DT; reactmicp::solver::Timestepper timestepper(MIN_DT, MAX_DT, 5000, 2.0); timestepper.get_options().alpha_average = 0.3; timestepper.get_options().decrease_failure = 0.1; timestepper.get_options().increase_factor = 1.05; timestepper.get_options().decrease_factor = 0.50; timestepper.get_options().iteration_lower_target = 1.001; - io::print_reactmicp_performance_long_header(&output_iter); + io::print_reactmicp_performance_long_header(output_iter); while (timestepper.get_total() < timestepper.get_total_target()) { output_gen << "dt : " << dt << std::endl; Timer step_timer; step_timer.start(); reactmicp::solver::ReactiveTransportReturnCode retcode = solver.solve_timestep(dt, variables); step_timer.stop(); // if (retcode > reactmicp::solver::ReactiveTransportReturnCode::NotConvergedYet) //{ - io::print_reactmicp_performance_long(&output_iter, cnt, timestepper.get_total()+dt, solver.get_perfs()); + io::print_reactmicp_performance_long(output_iter, cnt, timestepper.get_total()+dt, solver.get_perfs()); //} dt = timestepper.next_timestep(dt, retcode, solver.get_perfs().nb_iterations); if (retcode <= reactmicp::solver::ReactiveTransportReturnCode::NotConvergedYet) { dt = RESTART_DT; 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()); + io::print_reactmicp_performance_long(output_iter, cnt, timestepper.get_total(), solver.get_perfs()); } output_gen << "Time : " << timestepper.get_total() << std::endl; - io::print_reactmicp_performance_short(&output_gen, solver.get_perfs()); + io::print_reactmicp_performance_short(output_gen, solver.get_perfs()); scalar_t total = timestepper.get_total(); if (cnt % 10== 0) { output_x1 << total; output_x2 << total; output_x3 << total; output_x4 << total; output_sx1 << total; output_sx2 << total; output_sx3 << total; output_sx4 << total; output_c1 << total; output_c2 << total; output_s << total; for (index_t node: the_mesh->range_nodes()) { output_x1 << "\t" << variables->aqueous_concentration(node, 1, variables->displacement()); output_x2 << "\t" << variables->aqueous_concentration(node, 2, variables->displacement()); output_x3 << "\t" << variables->aqueous_concentration(node, 3, variables->displacement()); output_x4 << "\t" << variables->aqueous_concentration(node, 4, variables->displacement()); output_sx1 << "\t" << variables->solid_concentration(node, 1, variables->displacement()); output_sx2 << "\t" << variables->solid_concentration(node, 2, variables->displacement()); output_sx3 << "\t" << variables->solid_concentration(node, 3, variables->displacement()); output_sx4 << "\t" << variables->solid_concentration(node, 4, variables->displacement()); output_c1 << "\t" << variables->equilibrium_solution(node).secondary_molalities(0); output_c2 << "\t" << variables->equilibrium_solution(node).secondary_molalities(1); output_s << "\t" << pow10(variables->equilibrium_solution(node).main_variables(5)); } output_x1 << std::endl; output_x2 << std::endl; output_x3 << std::endl; output_x4 << std::endl; output_sx1 << std::endl; output_sx2 << std::endl; output_sx3 << std::endl; output_sx4 << std::endl; output_c1 << std::endl; output_c2 << std::endl; output_s << std::endl; } ++ cnt; } variables->displacement()(database->get_id_component("X1")) = 0; variables->displacement()(database->get_id_component("X2")) = -2.0; variables->displacement()(database->get_id_component("X3")) = 0; variables->displacement()(database->get_id_component("X4")) = 2.0; timestepper.set_total_target(6000); dt = RESTART_DT; transport_stagger->options_solver().absolute_tolerance = 1e-10; solver.get_options().absolute_residuals_tolerance = 1e-6; while (timestepper.get_total() < timestepper.get_total_target()) { output_gen << "dt : " << dt << std::endl; Timer step_timer; step_timer.start(); reactmicp::solver::ReactiveTransportReturnCode retcode = solver.solve_timestep(dt, variables); step_timer.stop(); // if (retcode > reactmicp::solver::ReactiveTransportReturnCode::NotConvergedYet) //{ - io::print_reactmicp_performance_long(&output_iter, cnt, timestepper.get_total()+dt, solver.get_perfs()); + io::print_reactmicp_performance_long(output_iter, cnt, timestepper.get_total()+dt, solver.get_perfs()); //} dt = timestepper.next_timestep(dt, retcode, solver.get_perfs().nb_iterations); if (retcode <= reactmicp::solver::ReactiveTransportReturnCode::NotConvergedYet) { dt = RESTART_DT; 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()); + io::print_reactmicp_performance_long(output_iter, cnt, timestepper.get_total(), solver.get_perfs()); } output_gen << "Time : " << timestepper.get_total() << std::endl; - io::print_reactmicp_performance_short(&output_gen, solver.get_perfs()); + io::print_reactmicp_performance_short(output_gen, solver.get_perfs()); scalar_t total = timestepper.get_total(); if (cnt % 10== 0) { output_x1 << total; output_x2 << total; output_x3 << total; output_x4 << total; output_sx1 << total; output_sx2 << total; output_sx3 << total; output_sx4 << total; output_c1 << total; output_c2 << total; output_s << total; for (index_t node: the_mesh->range_nodes()) { output_x1 << "\t" << variables->aqueous_concentration(node, 1, variables->displacement()); output_x2 << "\t" << variables->aqueous_concentration(node, 2, variables->displacement()); output_x3 << "\t" << variables->aqueous_concentration(node, 3, variables->displacement()); output_x4 << "\t" << variables->aqueous_concentration(node, 4, variables->displacement()); output_sx1 << "\t" << variables->solid_concentration(node, 1, variables->displacement()); output_sx2 << "\t" << variables->solid_concentration(node, 2, variables->displacement()); output_sx3 << "\t" << variables->solid_concentration(node, 3, variables->displacement()); output_sx4 << "\t" << variables->solid_concentration(node, 4, variables->displacement()); output_c1 << "\t" << variables->equilibrium_solution(node).secondary_molalities(0); output_c2 << "\t" << variables->equilibrium_solution(node).secondary_molalities(1); output_s << "\t" << pow10(variables->equilibrium_solution(node).main_variables(5)); } output_x1 << std::endl; output_x2 << std::endl; output_x3 << std::endl; output_x4 << std::endl; output_sx1 << std::endl; output_sx2 << std::endl; output_sx3 << std::endl; output_sx4 << std::endl; output_c1 << std::endl; output_c2 << std::endl; output_s << std::endl; } ++ cnt; } reactmicp::solver::ReactiveTransportTimer timer = solver.get_timer(); - io::print_reactmicp_timer(&output_gen, timer); + io::print_reactmicp_timer(output_gen, timer); global_timer.stop(); output_gen << "Total computation time : " << global_timer.elapsed_time() << "s." << std::endl; return EXIT_SUCCESS; } diff --git a/examples/reactmicp/saturated_react/react_leaching.cpp b/examples/reactmicp/saturated_react/react_leaching.cpp index 626752e..f0c8ead 100644 --- a/examples/reactmicp/saturated_react/react_leaching.cpp +++ b/examples/reactmicp/saturated_react/react_leaching.cpp @@ -1,721 +1,721 @@ #include "reactmicp.hpp" #include #include #include using namespace specmicp; // Initial conditions // ================== specmicp::RawDatabasePtr leaching_database() { specmicp::database::Database thedatabase("../data/cemdata_specmicp.js"); std::map 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 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) { 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; nodenb_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; nodenb_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; nodenb_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, 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) = 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) { scalar_t radius = 3.5; //cm scalar_t height = 8.0; //cm radius = radius + additional_nodes*dx; return mesh::uniform_axisymmetric_mesh1d(nb_nodes+additional_nodes, radius, dx, height); } 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 list_initial_composition = {initial, blank}; std::vector index_initial_state(nb_nodes, 0); for (int i=0; i< additional_nodes; ++i) index_initial_state[i] = 1; std::vector 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(nb_nodes, constraints, options); std::shared_ptr transport_stagger = std::make_shared(variables, list_fixed_nodes); dfpmsolver::ParabolicDriverOptions& transport_options = transport_stagger->options_solver(); set_transport_options(transport_options); UpscalingStaggerPtr upscaling_stagger = std::make_shared(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()); - std::ofstream output_iter("out_leaching_iter.dat"); + 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); + 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()); + 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()); + 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()); + 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()); + 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(&std::cout, solver.get_timer()); + io::print_reactmicp_timer(output_iter, solver.get_timer()); } diff --git a/src/reactmicp.hpp b/src/reactmicp.hpp index a5cb4b5..d1e5c0c 100644 --- a/src/reactmicp.hpp +++ b/src/reactmicp.hpp @@ -1,64 +1,65 @@ /*------------------------------------------------------------------------------- Copyright (c) 2014,2015 F. Georget , 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 reactmicp.hpp //! \brief Include this file to use ReactMiCP //! //! This file includes the main headers from ReactMiCP. #include "types.hpp" #include "specmicp.hpp" #include "reactmicp/systems/saturated_react/variables.hpp" #include "reactmicp/systems/saturated_react/equilibrium_stagger.hpp" #include "reactmicp/systems/saturated_react/transport_stagger.hpp" #include "reactmicp/systems/saturated_react/init_variables.hpp" #include "reactmicp/solver/reactive_transport_solver.hpp" #include "reactmicp/solver/reactive_transport_solver_structs.hpp" #include "reactmicp/solver/staggers_base/upscaling_stagger_base.hpp" #include "reactmicp/solver/staggers_base/stagger_structs.hpp" #include "reactmicp/solver/timestepper.hpp" #include "database/database.hpp" #include "dfpm/mesh.hpp" #include "reactmicp/solver/runner.hpp" +#include "utils/io/csv_formatter.hpp" #include "dfpm/io/meshes.hpp" #include "reactmicp/io/reactive_transport.hpp" #include "reactmicp/io/saturated_react.hpp" #include "physics/io/units.hpp" #include "utils/timer.hpp" diff --git a/src/reactmicp/io/reactive_transport.cpp b/src/reactmicp/io/reactive_transport.cpp index c51bc57..d3c1c2a 100644 --- a/src/reactmicp/io/reactive_transport.cpp +++ b/src/reactmicp/io/reactive_transport.cpp @@ -1,180 +1,180 @@ /*------------------------------------------------------------------------------- Copyright (c) 2014,2015 F. Georget , 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 "reactive_transport.hpp" #include #include "../../reactmicp/solver/reactive_transport_solver_structs.hpp" #include "../../utils/timer.hpp" - +#include "../../utils/io/csv_formatter.hpp" namespace specmicp { namespace io { //! \brief Print the ReactMiCP header void print_reactmicp_header( - std::ostream* output + OutFile& output ) { - (*output) << "// ================================== //\n" + output << "// ================================== //\n" << "// //\n" << "// ReactMiCP //\n" << "// //\n" << "// ================================== //\n\n" << "A reactive transport solver based on SpecMiCP.\n" - << "------------------------------------------- \n" - << std::endl; + << "------------------------------------------- \n\n"; + output.flush(); } //! \brief Print the performance information from the ReactMiCP solver void print_reactmicp_performance_short( - std::ostream* output, + OutFile& output, const reactmicp::solver::ReactiveTransportPerformance& perf ) { - (*output) << "------------------------------------------- \n" + output << "------------------------------------------- \n" << " Performance \n " << "------------------------------------------- \n\t" << " - Return code : " << reactmicp_return_code_to_string(perf.return_code) << "\n\t" << " - Residuals : " << perf.residuals << "\n\t" << " - Number of iterations : " << perf.nb_iterations << "\n" - << "------------------------------------------- " << std::endl; + << "------------------------------------------- \n"; } //! \brief Print the time spent in the different staggers void print_reactmicp_timer( - std::ostream* output, + OutFile& output, const reactmicp::solver::ReactiveTransportTimer& timer ) { - (*output) << "------------------------------------------- \n" + output << "------------------------------------------- \n" << "Time spent in each stagger \n" << "------------------------------------------- \n\t" << "- Transport : " << timer.transport_time << " s\n\t" << "- Chemistry : " << timer.chemistry_time << " s\n\t" << "- Upscaling : " << timer.upscaling_time << " s\n" - << "-----------------------------------------" << std::endl; + << "----------------------------------------- \n"; } //! \brief Print time and timer at the end of the computation void print_reactmicp_end( - std::ostream* output, + OutFile& output, const Timer& total_timer, const reactmicp::solver::ReactiveTransportTimer& timer ) { scalar_t elapsed_s = total_timer.elapsed_time(); index_t hours = static_cast(elapsed_s) / 3600; index_t elapsed_minutes = elapsed_s - hours*3600; index_t minute = elapsed_minutes / 60; index_t seconds = elapsed_minutes - 60*minute; - (*output) << " ====================================================== \n"; - (*output) << "computation finished at " - << total_timer.get_ctime_stop(); - (*output) << " Duration of the computation : " << total_timer.elapsed_time() << " s" + output << " ====================================================== \n"; + output << "computation finished at " + << total_timer.get_ctime_stop(); + output << " Duration of the computation : " << total_timer.elapsed_time() << " s" << "( " << hours << "h " << minute << "min " << seconds << "s )\n"; print_reactmicp_timer(output, timer); + + output.flush(); } //! \brief Print the header for the performance table //! //! \sa print_reactmicp_performance_long void print_reactmicp_performance_long_header( - std::ostream* output + OutFile& output ) { - (*output) << "Id\t T\t dt\t Return_code\t Iterations \t Residuals\t Time\t Transport_time\t Chemistry_time" - << std::endl; - + output << "Id\t T\t dt\t Return_code\t Iterations \t Residuals\t Time\t Transport_time\t Chemistry_time\n"; } //! \brief Print a row in the performance table //! //! The performance table contains information about an iteration of the reactmicp solver //! //! \sa print_reactmicp_performance_long_header void print_reactmicp_performance_long( - std::ostream* output, + OutFile& output, index_t cnt, scalar_t total, const reactmicp::solver::ReactiveTransportPerformance& perfs ) { - (*output) << cnt << "\t " + output << cnt << "\t " << total << "\t " << perfs.timestep << "\t " << (int) perfs.return_code << "\t " << perfs.nb_iterations << "\t " << perfs.residuals << "\t " << perfs.total_time << "\t" << perfs.transport_time << "\t" << perfs.chemistry_time - << std::endl; + << "\n"; } //! \brief Transform a ReactiveTransportreturnCode into a string std::string reactmicp_return_code_to_string( reactmicp::solver::ReactiveTransportReturnCode retcode) { using RetCode = reactmicp::solver::ReactiveTransportReturnCode; switch (retcode) { case RetCode::ResidualMinimized: return "ResidualMinimized"; case RetCode::ErrorMinimized: return "ErrorMinimized"; default: switch (retcode) { case RetCode::StationaryPoint: return "StationaryPoint"; case RetCode::MaximumIterationsReached: return "MaximumIterationsReached"; case RetCode::GoodEnough: return "Good enough..."; case RetCode::TransportBypass: return "Transport Bypass"; case RetCode::TransportFailure: return "TransportFailure"; case RetCode::ChemistryFailure: return "ChemistryFailure"; case RetCode::UpscalingFailure: return "UpscalingFailure"; default: return "Unknow error code !"; } } } } // end namespace io } // end namespace specmicp diff --git a/src/reactmicp/io/reactive_transport.hpp b/src/reactmicp/io/reactive_transport.hpp index 79cc512..8693f39 100644 --- a/src/reactmicp/io/reactive_transport.hpp +++ b/src/reactmicp/io/reactive_transport.hpp @@ -1,110 +1,114 @@ /*------------------------------------------------------------------------------- Copyright (c) 2014,2015 F. Georget , 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. -----------------------------------------------------------------------------*/ #ifndef SPECMICP_IO_REACTMICP_HPP #define SPECMICP_IO_REACTMICP_HPP //! \file reactmicp/io/reactive_transport.hpp //! \brief Print information from the reactive transport module -#include + #include "../../types.hpp" +#include + namespace specmicp { // forward declaration class Timer; namespace reactmicp { namespace solver { struct ReactiveTransportPerformance; struct ReactiveTransportTimer; enum class ReactiveTransportReturnCode; } //end namespace solver } //end namespace reactmicp namespace io { +class OutFile; + //! \brief Print the ReactMiCP header void SPECMICP_DLL_PUBLIC print_reactmicp_header( - std::ostream* output + OutFile& output ); //! \brief Print the performance information from the ReactMiCP solver void SPECMICP_DLL_PUBLIC print_reactmicp_performance_short( - std::ostream* output, + OutFile& output, const reactmicp::solver::ReactiveTransportPerformance& perf ); //! \brief Print the time spent in the different staggers void SPECMICP_DLL_PUBLIC print_reactmicp_timer( - std::ostream* output, + OutFile& output, const reactmicp::solver::ReactiveTransportTimer& timer ); //! \brief Print time and timer at the end of the computation void SPECMICP_DLL_PUBLIC print_reactmicp_end( - std::ostream* output, + OutFile& output, const Timer& total_timer, const reactmicp::solver::ReactiveTransportTimer& timer ); //! \brief Print the header for the performance table //! //! \sa print_reactmicp_performance_long void SPECMICP_DLL_PUBLIC print_reactmicp_performance_long_header( - std::ostream* output + OutFile& output ); //! \brief Print a row in the performance table //! //! The performance table contains information about an iteration of the reactmicp solver //! //! \sa print_reactmicp_performance_long_header void SPECMICP_DLL_PUBLIC print_reactmicp_performance_long( - std::ostream* output, + OutFile& output, index_t cnt, scalar_t total, const reactmicp::solver::ReactiveTransportPerformance& perfs ); //! \brief Transform a ReactiveTransportreturnCode into a string std::string SPECMICP_DLL_PUBLIC reactmicp_return_code_to_string( reactmicp::solver::ReactiveTransportReturnCode retcode); } // end namespace io } // end namespace specmicp #endif // SPECMICP_IO_REACTMICP_HPP diff --git a/src/reactmicp/io/saturated_react.cpp b/src/reactmicp/io/saturated_react.cpp index 14b0add..d949144 100644 --- a/src/reactmicp/io/saturated_react.cpp +++ b/src/reactmicp/io/saturated_react.cpp @@ -1,206 +1,224 @@ /*------------------------------------------------------------------------------- Copyright (c) 2015 F. Georget , 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 "saturated_react.hpp" #include #include #include #include #include "../../dfpm/meshes/mesh1d.hpp" #include "../../specmicp/adimensional/adimensional_system_solution_extractor.hpp" #include "../../utils/dateandtime.hpp" #include "../../physics/io/units.hpp" #include "../systems/saturated_react/variables.hpp" +#include "../../utils/io/csv_formatter.hpp" using namespace specmicp::reactmicp::systems::satdiff; namespace specmicp { namespace io { static constexpr char field_sep = '\t'; static constexpr char comment = '#'; -//! \brief Print the time as a comment -void print_csv_header(std::ostream* ofile) +void print_csv_header(CSVFile& ofile) { - (*ofile) << comment << " - ReactMiCP - " << dateandtime::now(); + ofile.insert_comment_line("ReactMiCP"); + ofile.insert_comment_line("---------"); + ofile.insert_comment_date(); } -//! \brief Print the total aqueous concentrations in the CSV format in 'ofile' void print_components_total_aqueous_concentration( RawDatabasePtr the_database, SaturatedVariablesPtr variables, mesh::Mesh1DPtr the_mesh, const units::UnitsSet& units_set, - std::ostream* ofile + CSVFile& ofile ) { print_csv_header(ofile); - (*ofile) << comment << "units : mol/" << io::volume_unit_to_string(units_set.length) << std::endl; - (*ofile) << "Position"; + ofile.insert_comment_line("Total aqueous concentrations profiles"); + ofile.insert_comment_line( "length unit : " + io::length_unit_to_string(units_set.length)); + ofile.insert_comment_line( "concentration unit : mol/" + io::volume_unit_to_string(units_set.length)); + ofile << "Position"; for (index_t component: the_database->range_aqueous_component()) { - (*ofile) << field_sep << the_database->get_label_component(component); + ofile.separator(); + ofile << the_database->get_label_component(component); } - (*ofile) << std::endl; + ofile.eol(); for (index_t node: the_mesh->range_nodes()) { - (*ofile) << the_mesh->get_position(node); + ofile << the_mesh->get_position(node); for (index_t component: the_database->range_aqueous_component()) { - (*ofile) << field_sep << variables->aqueous_concentration(node, component); + ofile.separator(); + ofile << variables->aqueous_concentration(node, component); } - (*ofile) << std::endl; + ofile.eol(); } + ofile.flush(); } //! \brief Print the total aqueous concentrations in the CSV format in the file 'filepath' void print_components_total_aqueous_concentration( RawDatabasePtr the_database, SaturatedVariablesPtr variables, mesh::Mesh1DPtr the_mesh, const units::UnitsSet& units_set, const std::string& filepath ) { - std::ofstream ofile(filepath); - print_components_total_aqueous_concentration(the_database, variables, the_mesh, units_set, &ofile); - ofile.close(); + CSVFile ofile(filepath); + print_components_total_aqueous_concentration(the_database, variables, the_mesh, units_set, ofile); } //! \brief Print the total solid concentrations in the CSV format in 'ofile' void print_components_total_solid_concentration( RawDatabasePtr the_database, SaturatedVariablesPtr variables, mesh::Mesh1DPtr the_mesh, const units::UnitsSet& units_set, - std::ostream* ofile + CSVFile& ofile ) { print_csv_header(ofile); - (*ofile) << comment << "units : mol/" << io::volume_unit_to_string(units_set.length) << std::endl; - (*ofile) << "Position"; + ofile.insert_comment_line( "length unit : " + io::length_unit_to_string(units_set.length)); + ofile.insert_comment_line( "concentration unit : mol/" + io::volume_unit_to_string(units_set.length)); + + ofile << "Position"; for (index_t component: the_database->range_aqueous_component()) { - (*ofile) << field_sep << the_database->get_label_component(component); + ofile.separator(); + ofile << the_database->get_label_component(component); } - (*ofile) << std::endl; + ofile.eol(); for (index_t node: the_mesh->range_nodes()) { - (*ofile) << the_mesh->get_position(node); + ofile << the_mesh->get_position(node); for (index_t component: the_database->range_aqueous_component()) { - (*ofile) << field_sep << variables->solid_concentration(node, component); + ofile.separator(); + ofile << variables->solid_concentration(node, component); } - (*ofile) << std::endl; + ofile.eol(); } + ofile.flush(); } //! \brief Print the total solid concentrations in the CSV format in the file 'filepath' void print_components_total_solid_concentration( RawDatabasePtr the_database, SaturatedVariablesPtr variables, mesh::Mesh1DPtr the_mesh, const units::UnitsSet& units_set, const std::string& filepath ) { - std::ofstream ofile(filepath); - print_components_total_solid_concentration(the_database, variables, the_mesh, units_set, &ofile); - ofile.close(); + CSVFile ofile(filepath); + print_components_total_solid_concentration(the_database, variables, the_mesh, units_set, ofile); + } //! \brief Print the solid phases profiles in 'ofile' -inline void print_minerals_profile( - RawDatabasePtr the_database, +inline void print_minerals_profile(RawDatabasePtr the_database, SaturatedVariablesPtr variables, mesh::Mesh1DPtr the_mesh, - std::ostream* ofile + const units::UnitsSet& units_set, + CSVFile& ofile ) { print_csv_header(ofile); - (*ofile) << "Position"; + ofile.insert_comment_line("Solid phase profiles"); + ofile.insert_comment_line( "length unit : " + io::length_unit_to_string(units_set.length)); + ofile << "Position"; for (index_t mineral: the_database->range_mineral()) { - (*ofile) << field_sep << the_database->get_label_mineral(mineral); + ofile.separator(); + ofile << the_database->get_label_mineral(mineral); } - (*ofile) << field_sep << "Porosity"; - (*ofile) << field_sep << "pH"; - (*ofile) << std::endl; + ofile.separator(); + ofile << "Porosity"; + ofile.separator(); + ofile << "pH"; + ofile.eol(); for (index_t node: the_mesh->range_nodes()) { - (*ofile) << the_mesh->get_position(node); + 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) << field_sep << extractor.volume_fraction_mineral(mineral); + ofile.separator(); + ofile << extractor.volume_fraction_mineral(mineral); } - (*ofile) << field_sep << variables->porosity(node); - (*ofile) << field_sep << extractor.pH(); - (*ofile) << std::endl; - + ofile.separator(); + ofile << variables->porosity(node); + ofile.separator(); + ofile << extractor.pH(); + ofile.eol(); } - + ofile.flush(); } //! \brief Print the solid phases profiles in 'filepath' void print_minerals_profile( RawDatabasePtr the_database, SaturatedVariablesPtr variables, mesh::Mesh1DPtr the_mesh, + const units::UnitsSet& units_set, const std::string& filepath ) { - std::ofstream ofile(filepath); + CSVFile ofile(filepath); print_minerals_profile( the_database, variables, the_mesh, - &ofile + units_set, + ofile ); - ofile.close(); } } // end namespace io } // end namespace specmicp diff --git a/src/reactmicp/io/saturated_react.hpp b/src/reactmicp/io/saturated_react.hpp index c87a222..6e8437d 100644 --- a/src/reactmicp/io/saturated_react.hpp +++ b/src/reactmicp/io/saturated_react.hpp @@ -1,129 +1,109 @@ /*------------------------------------------------------------------------------- Copyright (c) 2014,2015 F. Georget , 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. -----------------------------------------------------------------------------*/ #ifndef SPECMICP_IO_SATURATEDREACT_HPP #define SPECMICP_IO_SATURATEDREACT_HPP //! \file saturated_react.hpp //! \brief Print information about the saturated reactive transport module #include "../../types.hpp" #include #include namespace specmicp { // forward declarations namespace database { struct DataContainer; } //end namespace database using RawDatabasePtr = std::shared_ptr; namespace units { struct UnitsSet; } //end namespace units namespace mesh { class Mesh1D; using Mesh1DPtr = std::shared_ptr; } //end namespace mesh namespace reactmicp { namespace systems { namespace satdiff { class SaturatedVariables; using SaturatedVariablesPtr = std::shared_ptr; } //end namespace satdiff } //end namespace systems } //end namespace reactmicp namespace io { using namespace specmicp::reactmicp::systems::satdiff; - -//! \brief Print the time as a comment -void SPECMICP_DLL_PUBLIC print_csv_header(std::ostream* ofile); - -//! \brief Print the total aqueous concentrations in the CSV format in 'ofile' -void SPECMICP_DLL_PUBLIC print_components_total_aqueous_concentration( - RawDatabasePtr the_database, - SaturatedVariablesPtr variables, - mesh::Mesh1DPtr the_mesh, - const units::UnitsSet& units_set, - std::ostream* ofile - ); - //! \brief Print the total aqueous concentrations in the CSV format in the file 'filepath' void SPECMICP_DLL_PUBLIC print_components_total_aqueous_concentration( RawDatabasePtr the_database, SaturatedVariablesPtr variables, mesh::Mesh1DPtr the_mesh, const units::UnitsSet& units_set, const std::string& filepath ); //! \brief Print the total solid concentrations in the CSV format in the file 'filepath' void SPECMICP_DLL_PUBLIC print_components_total_solid_concentration( RawDatabasePtr the_database, SaturatedVariablesPtr variables, mesh::Mesh1DPtr the_mesh, const units::UnitsSet& units_set, const std::string& filepath ); -//! \brief Print the solid phases profiles in 'ofile' -void SPECMICP_DLL_PUBLIC print_minerals_profile( - RawDatabasePtr the_database, - SaturatedVariablesPtr variables, - mesh::Mesh1DPtr the_mesh, - std::ostream* ofile - ); - //! \brief Print the solid phases profiles in 'filepath' void SPECMICP_DLL_PUBLIC print_minerals_profile( RawDatabasePtr the_database, SaturatedVariablesPtr variables, mesh::Mesh1DPtr the_mesh, + const units::UnitsSet& units_set, const std::string& filepath ); } // end namespace io } // end namespace specmicp #endif // SPECMICP_IO_SATURATEDREACT_HPP diff --git a/src/reactmicp/solver/runner.cpp b/src/reactmicp/solver/runner.cpp index a0cb5ad..58e52a6 100644 --- a/src/reactmicp/solver/runner.cpp +++ b/src/reactmicp/solver/runner.cpp @@ -1,112 +1,108 @@ /*------------------------------------------------------------------------------- Copyright (c) 2014,2015 F. Georget , 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 "runner.hpp" #include "../io/reactive_transport.hpp" #include "../../utils/timer.hpp" #include #include "staggers_base/variables_base.hpp" +#include "../../utils/io/csv_formatter.hpp" + namespace specmicp { namespace reactmicp { namespace solver { void ReactiveTransportRunner::run_until(scalar_t target, VariablesBasePtr variables) { Timer total_time; total_time.start(); m_timestepper.set_total_target(target); m_output_target = m_timestepper.get_total()+m_simulinfo.output_step; - std::ofstream out_iter; + io::OutFile out_iter(m_simulinfo.complete_filepath(std::string("iter"), std::string("dat"))); + io::print_reactmicp_header(out_iter); if (m_simulinfo.print_iter_info) { - out_iter.open(m_simulinfo.complete_filepath(std::string("iter"), std::string("dat"))); - io::print_reactmicp_header(&out_iter); - io::print_reactmicp_performance_long_header(&out_iter); + io::print_reactmicp_performance_long_header(out_iter); } scalar_t dt = m_timestepper.get_options().restart_timestep; ReactiveTransportReturnCode retcode = ReactiveTransportReturnCode::NotConvergedYet; while(retcode >= ReactiveTransportReturnCode::NotConvergedYet and m_timestepper.get_total() < m_timestepper.get_total_target()) { Timer step_timer; step_timer.start(); reactmicp::solver::ReactiveTransportReturnCode retcode = m_solver.solve_timestep(dt, variables); step_timer.stop(); if (m_simulinfo.print_iter_info) - io::print_reactmicp_performance_long(&out_iter, m_cnt, m_timestepper.get_total()+dt, m_solver.get_perfs()); + io::print_reactmicp_performance_long(out_iter, m_cnt, m_timestepper.get_total()+dt, m_solver.get_perfs()); dt = m_timestepper.next_timestep(dt, retcode, m_solver.get_perfs().nb_iterations); if (retcode <= reactmicp::solver::ReactiveTransportReturnCode::NotConvergedYet) { dt = m_timestepper.get_options().restart_timestep; variables->reset_main_variables(); retcode = m_solver.solve_timestep(dt, variables); if (m_simulinfo.print_iter_info) - io::print_reactmicp_performance_long(&out_iter, m_cnt, m_timestepper.get_total()+dt, m_solver.get_perfs()); + io::print_reactmicp_performance_long(out_iter, m_cnt, m_timestepper.get_total()+dt, m_solver.get_perfs()); dt = m_timestepper.next_timestep(dt, retcode, m_solver.get_perfs().nb_iterations); } ++m_cnt; // output if (m_timestepper.get_total() > m_output_target) { m_output_function(m_timestepper.get_total(), variables); m_output_target += m_simulinfo.output_step; } if (retcode == reactmicp::solver::ReactiveTransportReturnCode::UserTermination) break; } total_time.stop(); - if (m_simulinfo.print_iter_info) - { + io::print_reactmicp_end(out_iter, total_time, get_timer()); - io::print_reactmicp_end(&out_iter, total_time, get_timer()); - out_iter.close(); - - } } } // end namespace solver } // end namespace reactmicp } // end namespace specmicp diff --git a/src/utils/CMakeLists.txt b/src/utils/CMakeLists.txt index 63a9f08..d5187c4 100644 --- a/src/utils/CMakeLists.txt +++ b/src/utils/CMakeLists.txt @@ -1,83 +1,89 @@ add_custom_target(utils_inc SOURCES log.hpp # sparse solvers # -------------- sparse_solvers/sparse_solver_base.hpp sparse_solvers/sparse_solver.hpp sparse_solvers/sparse_solver_structs.hpp sparse_solvers/sparse_qr.hpp sparse_solvers/sparse_lu.hpp sparse_solvers/sparse_bicgstab.hpp sparse_solvers/sparse_gmres.hpp options_handler.hpp perfs_handler.hpp ) set(SPECMICP_COMMON_LIBRARY_FILES dateandtime.cpp timer.cpp moving_average.cpp json.cpp ../physics/laws.cpp ../physics/units.cpp ../physics/io/units.cpp + io/csv_formatter.cpp + ${JSONCPP_DIR}/jsoncpp.cpp ) add_library(specmicp_common SHARED ${SPECMICP_COMMON_LIBRARY_FILES}) install(TARGETS specmicp_common LIBRARY DESTINATION ${LIBRARY_INSTALL_DIR} ) set(UTILS_INCLUDE_LIST log.hpp options_handler.hpp perfs_handler.hpp moving_average.hpp timer.hpp dateandtime.hpp json.hpp ) set(UTILS_SPARSE_INCLUDE_LIST # sparse solvers # -------------- sparse_solvers/sparse_solver.hpp sparse_solvers/sparse_solver_base.hpp sparse_solvers/sparse_solver_structs.hpp sparse_solvers/sparse_qr.hpp sparse_solvers/sparse_lu.hpp sparse_solvers/sparse_bicgstab.hpp sparse_solvers/sparse_gmres.hpp ) install(FILES ${UTILS_INCLUDE_LIST} DESTINATION ${INCLUDE_INSTALL_DIR}/utils ) install(FILES ${UTILS_SPARSE_INCLUDE_LIST} DESTINATION ${INCLUDE_INSTALL_DIR}/utils/sparse_solvers ) +install(FILES io/csv_formatter.hpp + DESTINATION ${INCLUDE_INSTALL_DIR}/io +) + # static libraries # ---------------- if(SPECMICP_BUILD_STATIC) add_library(specmicp_common_static STATIC ${SPECMICP_COMMON_LIBRARY_FILES}) install(TARGETS specmicp_common_static ARCHIVE DESTINATION ${STATIC_LIBRARY_INSTALL_DIR} ) else() add_library(specmicp_common_static EXCLUDE_FROM_ALL STATIC ${SPECMICP_COMMON_LIBRARY_FILES}) endif() set_target_properties(specmicp_common_static PROPERTIES OUTPUT_NAME specmicp_common) diff --git a/src/reactmicp.hpp b/src/utils/io/csv_formatter.cpp similarity index 59% copy from src/reactmicp.hpp copy to src/utils/io/csv_formatter.cpp index a5cb4b5..8ec4f1e 100644 --- a/src/reactmicp.hpp +++ b/src/utils/io/csv_formatter.cpp @@ -1,64 +1,82 @@ /*------------------------------------------------------------------------------- -Copyright (c) 2014,2015 F. Georget , Princeton University +Copyright (c) 2015 F. Georget , 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 reactmicp.hpp -//! \brief Include this file to use ReactMiCP -//! -//! This file includes the main headers from ReactMiCP. +#include "csv_formatter.hpp" +#include "../dateandtime.hpp" +#include -#include "types.hpp" +namespace specmicp { +namespace io { -#include "specmicp.hpp" +std::string CSVFormatter::comment {"#"}; +std::string CSVFormatter::separator {"\t"}; +std::string CSVFormatter::space {" "}; +std::string CSVFormatter::eol {"\n"}; -#include "reactmicp/systems/saturated_react/variables.hpp" -#include "reactmicp/systems/saturated_react/equilibrium_stagger.hpp" -#include "reactmicp/systems/saturated_react/transport_stagger.hpp" -#include "reactmicp/systems/saturated_react/init_variables.hpp" -#include "reactmicp/solver/reactive_transport_solver.hpp" -#include "reactmicp/solver/reactive_transport_solver_structs.hpp" -#include "reactmicp/solver/staggers_base/upscaling_stagger_base.hpp" -#include "reactmicp/solver/staggers_base/stagger_structs.hpp" -#include "reactmicp/solver/timestepper.hpp" +void CSVFormatter::insert_comment_line(std::ostream* out, const std::string& msg) +{ + (*out) << comment << space << msg << eol; +} -#include "database/database.hpp" +void CSVFormatter::insert_comment_date(std::ostream* out) +{ + insert_comment_line(out, dateandtime::now()); +} -#include "dfpm/mesh.hpp" -#include "reactmicp/solver/runner.hpp" +OutFile::OutFile(const std::string& filepath): + m_stream(filepath) +{ + if (not m_stream) { + throw std::runtime_error("Cannot open the file : '"+filepath+"'."); + } +} -#include "dfpm/io/meshes.hpp" -#include "reactmicp/io/reactive_transport.hpp" -#include "reactmicp/io/saturated_react.hpp" -#include "physics/io/units.hpp" +void OutFile::open(const std::string& filepath) +{ + if (m_stream.is_open()) close(); + m_stream.open(filepath); +} -#include "utils/timer.hpp" +void OutFile::close() +{ + if (m_stream.is_open()) m_stream.flush(); + m_stream.close(); +} + +CSVFile::CSVFile(const std::string& filepath): + OutFile(filepath) +{} + +} //end namespace io +} //end namespace specmicp diff --git a/src/utils/io/csv_formatter.hpp b/src/utils/io/csv_formatter.hpp new file mode 100644 index 0000000..11eafff --- /dev/null +++ b/src/utils/io/csv_formatter.hpp @@ -0,0 +1,156 @@ +/*------------------------------------------------------------------------------- + +Copyright (c) 2014,2015 F. Georget , 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. + +-----------------------------------------------------------------------------*/ + +#ifndef SPECMICP_IO_CSV_FORMATTER +#define SPECMICP_IO_CSV_FORMATTER + +#include "../../types.hpp" + +#include +#include + +#include + +namespace specmicp { +namespace io { + +//! \brief A CSV Formatter +class CSVFormatter +{ +public: + static std::string comment; + static std::string separator; + static std::string space; + static std::string eol; + + //! \brief insert a comment line + static void insert_comment_line(std::ostream* out, const std::string& msg); + //! \brief insert a comment line with the date + static void insert_comment_date(std::ostream* out); + +}; + +//! \brief A output file +//! +//! This is a wrapper around a std::ofstream +//! +//! The wrapper exists to make the formatting easier +class SPECMICP_DLL_PUBLIC OutFile +{ +public: + //! \brief Open a file to write + OutFile(const std::string& filepath); + + //! \brief Close the file + ~OutFile() { + close(); + } + + //! \brief open a file + void open(const std::string& filepath); + + //! \brief Stream operator + std::ofstream& operator<< (const std::string& value) { + m_stream << value; + return m_stream; + } + + //! \brief Stream operator + std::ofstream& operator<< (scalar_t value) { + m_stream << value; + return m_stream; + } + + //! \brief Stream operator + std::ofstream& operator<< (index_t value) { + m_stream << value; + return m_stream; + } + + //! \brief Close a line and flush the stream + void endl() { + m_stream << std::endl; + } + + //! \brief Flush the stream + void flush() { + m_stream.flush(); + } + + //! \brief close the stream + void close(); + +protected: + std::ofstream m_stream; +}; + +//! \brief A CSV file +class SPECMICP_DLL_PUBLIC CSVFile: public OutFile +{ +public: + CSVFile(const std::string& filepath); + + + //! \brief insert a comment line + void insert_comment_line(const std::string& msg) { + return CSVFormatter::insert_comment_line(&m_stream, msg); + } + //! \brief insert a comment line with the date + void insert_comment_date() { + return CSVFormatter::insert_comment_date(&m_stream); + } + + //! \brief Write an end of line + void eol() { + m_stream << CSVFormatter::eol; + } + + //! \brief Write the space sequence + void space() { + m_stream << CSVFormatter::space; + } + + //! \brief Write the separator sequence + void separator() { + m_stream << CSVFormatter::separator; + } + + //! \brief Write the comment sequence + void comment() { + m_stream << CSVFormatter::comment; + } +}; + +} //end namespace io +} //end namespace specmicp + +#endif // SPECMICP_IO_CSV_FORMATTER