diff --git a/examples/CMakeLists.txt b/examples/CMakeLists.txt index 542d65c..ebd2eb2 100644 --- a/examples/CMakeLists.txt +++ b/examples/CMakeLists.txt @@ -1,128 +1,136 @@ ###### Installation ####################### set(EXAMPLES_INSTALL_DIR ${SHARE_INSTALL_DIR}/examples) # SpecMiCP # ======== set (SPECMICP_EXAMPLES specmicp/adimensional/thermocarbo.cpp specmicp/adimensional/equilibrium_curve.cpp specmicp/adimensional/momas_thermo.cpp specmicp/adimensional/carbofe.cpp ) install(FILES ${SPECMICP_EXAMPLES} DESTINATION ${EXAMPLES_INSTALL_DIR}/specmicp ) # ReactMiCP # ========= set (REACTMICP_SATURATED_EXAMPLES reactmicp/equilibrium_curve.cpp reactmicp/saturated_react/react_leaching.cpp reactmicp/saturated_react/carbonation.cpp reactmicp/saturated_react/carbonationfe.cpp + reactmicp/saturated_react/leaching_kinetic.cpp ) install(FILES ${REACTMICP_SATURATED_EXAMPLES} DESTINATION ${EXAMPLES_INSTALL_DIR}/reactmicp/saturated ) ###### Build (optional) #################### if (NOT SPECMICP_BUILD_EXAMPLE) set(EXAMPLE_IS_OPTIONAL EXCLUDE_FROM_ALL) else() set(EXAMPLE_IS_OPTIONAL ) endif() # SpecMiCP # ======== # thermocarbo add_executable(ex_adim_thermocarbo ${EXAMPLE_IS_OPTIONAL} specmicp/adimensional/thermocarbo.cpp ) target_link_libraries(ex_adim_thermocarbo ${SPECMICP_LIBS}) # carbofe add_executable(ex_adim_carbofe ${EXAMPLE_IS_OPTIONAL} specmicp/adimensional/carbofe.cpp ) target_link_libraries(ex_adim_carbofe ${SPECMICP_LIBS}) # Momas thermodynamic example add_executable(ex_adim_momas ${EXAMPLE_IS_OPTIONAL} specmicp/adimensional/momas_thermo.cpp ) target_link_libraries(ex_adim_momas ${SPECMICP_LIBS}) # equilibrium curve add_executable(ex_adim_equilibriumcurve ${EXAMPLE_IS_OPTIONAL} specmicp/adimensional/equilibrium_curve.cpp ) target_link_libraries(ex_adim_equilibriumcurve ${SPECMICP_LIBS}) # ReactMiCP # ========= # EquilibriumCurve add_executable(ex_reactmicp_equilibriumcurve ${EXAMPLE_IS_OPTIONAL} reactmicp/equilibrium_curve.cpp ) target_link_libraries(ex_reactmicp_equilibriumcurve ${REACTMICP_LIBS}) -# Leaching +# Leaching add_executable(ex_reactmicp_leaching ${EXAMPLE_IS_OPTIONAL} reactmicp/saturated_react/react_leaching.cpp ) target_link_libraries(ex_reactmicp_leaching ${REACTMICP_LIBS}) -# Momas benchmark +# Momas benchmark add_executable(ex_reactmicp_momas ${EXAMPLE_IS_OPTIONAL} reactmicp/saturated_react/momas_benchmark.cpp ) target_link_libraries(ex_reactmicp_momas ${REACTMICP_LIBS}) -# Carbonation +# Carbonation add_executable(ex_reactmicp_carbo ${EXAMPLE_IS_OPTIONAL} reactmicp/saturated_react/carbonation.cpp ) target_link_libraries(ex_reactmicp_carbo ${REACTMICP_LIBS}) add_executable(ex_reactmicp_carbo_static EXCLUDE_FROM_ALL reactmicp/saturated_react/carbonation.cpp ) target_link_libraries(ex_reactmicp_carbo_static ${REACTMICP_STATIC_LIBS}) -# Carbonation with Fe +# Carbonation with Fe add_executable(ex_reactmicp_carbofe ${EXAMPLE_IS_OPTIONAL} reactmicp/saturated_react/carbonationfe.cpp ) target_link_libraries(ex_reactmicp_carbofe ${REACTMICP_LIBS}) +# Leaching with hydration + +add_executable(ex_reactmicp_leaching_kinetic + ${EXAMPLE_IS_OPTIONAL} + reactmicp/saturated_react/leaching_kinetic.cpp +) +target_link_libraries(ex_reactmicp_leaching_kinetic ${REACTMICP_STATIC_LIBS}) diff --git a/examples/reactmicp/saturated_react/leaching_kinetic.cpp b/examples/reactmicp/saturated_react/leaching_kinetic.cpp new file mode 100644 index 0000000..2e2d66b --- /dev/null +++ b/examples/reactmicp/saturated_react/leaching_kinetic.cpp @@ -0,0 +1,843 @@ +// ============================================ + +// 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); + //! \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_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/src/reactmicp/CMakeLists.txt b/src/reactmicp/CMakeLists.txt index a14893f..64e345a 100644 --- a/src/reactmicp/CMakeLists.txt +++ b/src/reactmicp/CMakeLists.txt @@ -1,126 +1,131 @@ # ReactMiCP # ============= set(REACTMICP_SOLVER_DIR solver) set(REACTMICP_SYSTEMS_DIR systems) set(REACTMICP_EQUILIBRIUMCURVE equilibrium_curve) add_custom_target(reactmicp_incl SOURCES ${REACTMICP_SOLVER_DIR}/reactive_transport_solver_structs.hpp ${REACTMICP_SOLVER_DIR}/staggers_base/variables_base.hpp ${REACTMICP_SOLVER_DIR}/staggers_base/transport_stagger_base.hpp ${REACTMICP_SOLVER_DIR}/staggers_base/chemistry_stagger_base.hpp ${REACTMICP_SOLVER_DIR}/staggers_base/upscaling_stagger_base.hpp ${REACTMICP_SOLVER_DIR}/staggers_base/stagger_structs.hpp ${REACTMICP_SOLVER_DIR}/staggers_base/decl.inl ${REACTMICP_SOLVER_DIR}/staggers_base/staggers_base.hpp ${REACTMICP_SYSTEMS_DIR}/saturated_react/variablesfwd.hpp ) set(REACTMICP_LIBFILES ${REACTMICP_EQUILIBRIUMCURVE}/chemistry.cpp ${REACTMICP_EQUILIBRIUMCURVE}/eqcurve_extractor.cpp ${REACTMICP_EQUILIBRIUMCURVE}/eqcurve_coupler.cpp ${REACTMICP_EQUILIBRIUMCURVE}/eqcurve_solid_transport.cpp # Flexible reactive transport solver # ---------------------------------- ${REACTMICP_SOLVER_DIR}/reactive_transport_solver.cpp ${REACTMICP_SOLVER_DIR}/timestepper.cpp ${REACTMICP_SOLVER_DIR}/runner.cpp # ${REACTMICP_SYSTEMS_DIR}/saturated_react/variables.cpp ${REACTMICP_SYSTEMS_DIR}/saturated_react/transport_program.cpp ${REACTMICP_SYSTEMS_DIR}/saturated_react/transport_stagger.cpp ${REACTMICP_SYSTEMS_DIR}/saturated_react/equilibrium_stagger.cpp + ${REACTMICP_SYSTEMS_DIR}/saturated_react/kinetic_stagger.cpp ${REACTMICP_SYSTEMS_DIR}/saturated_react/react_solver.cpp ${REACTMICP_SYSTEMS_DIR}/saturated_react/init_variables.cpp ${REACTMICP_SYSTEMS_DIR}/saturated_react/diffusive_upscaling_stagger.cpp ) if(OPENMP_FOUND) if(SPECMICP_USE_OPENMP) set_source_files_properties(${REACTMICP_SYSTEMS_DIR}/saturated_react/equilibrium_stagger.cpp PROPERTIES COMPILE_DEFINITIONS "SPECMICP_USE_OPENMP" ) + set_source_files_properties(${REACTMICP_SYSTEMS_DIR}/saturated_react/kinetic_stagger.cpp + PROPERTIES COMPILE_DEFINITIONS "SPECMICP_USE_OPENMP" + ) endif() endif() add_library(reactmicp SHARED ${REACTMICP_LIBFILES}) target_link_libraries(reactmicp dfpm) install(TARGETS reactmicp LIBRARY DESTINATION ${LIBRARY_INSTALL_DIR} ) # include files # ------------- set(REACTMICP_STAGGERBASE_INCLUDE_LIST ${REACTMICP_SOLVER_DIR}/staggers_base/variables_base.hpp ${REACTMICP_SOLVER_DIR}/staggers_base/transport_stagger_base.hpp ${REACTMICP_SOLVER_DIR}/staggers_base/chemistry_stagger_base.hpp ${REACTMICP_SOLVER_DIR}/staggers_base/upscaling_stagger_base.hpp ${REACTMICP_SOLVER_DIR}/staggers_base/stagger_structs.hpp ${REACTMICP_SOLVER_DIR}/staggers_base/decl.inl ${REACTMICP_SOLVER_DIR}/staggers_base/staggers_base.hpp ) set(REACTMICP_SOLVER_INCLUDE_LIST ${REACTMICP_SOLVER_DIR}/reactive_transport_solver_structs.hpp ${REACTMICP_SOLVER_DIR}/reactive_transport_solver.hpp ${REACTMICP_SOLVER_DIR}/timestepper.hpp ${REACTMICP_SOLVER_DIR}/runner.hpp ) set(REACTMICP_SATURATEDREACT_INCLUDE_LIST ${REACTMICP_SYSTEMS_DIR}/saturated_react/variablesfwd.hpp ${REACTMICP_SYSTEMS_DIR}/saturated_react/variables.hpp ${REACTMICP_SYSTEMS_DIR}/saturated_react/transport_program.hpp ${REACTMICP_SYSTEMS_DIR}/saturated_react/transport_stagger.hpp ${REACTMICP_SYSTEMS_DIR}/saturated_react/equilibrium_stagger.hpp + ${REACTMICP_SYSTEMS_DIR}/saturated_react/kinetic_stagger.hpp ${REACTMICP_SYSTEMS_DIR}/saturated_react/react_solver.hpp ${REACTMICP_SYSTEMS_DIR}/saturated_react/init_variables.hpp ${REACTMICP_SYSTEMS_DIR}/saturated_react/diffusive_upscaling_stagger.hpp ) set(REACTMICP_EQUILIBRIUMCURVE_INCLUDE_LIST ${REACTMICP_EQUILIBRIUMCURVE}/chemistry.cpp ${REACTMICP_EQUILIBRIUMCURVE}/eqcurve_extractor.cpp ${REACTMICP_EQUILIBRIUMCURVE}/eqcurve_coupler.cpp ${REACTMICP_EQUILIBRIUMCURVE}/eqcurve_solid_transport.cpp ) install(FILES ${REACTMICP_STAGGERBASE_INCLUDE_LIST} DESTINATION ${INCLUDE_INSTALL_DIR}/reactmicp/solver/staggers_base ) install(FILES ${REACTMICP_SOLVER_INCLUDE_LIST} DESTINATION ${INCLUDE_INSTALL_DIR}/reactmicp/solver ) install(FILES ${REACTMICP_SATURATEDREACT_INCLUDE_LIST} DESTINATION ${INCLUDE_INSTALL_DIR}/reactmicp/systems/saturated_react ) install(FILES ${REACTMICP_EQUILIBRIUMCURVE_INCLUDE_LIST} DESTINATION ${INCLUDE_INSTALL_DIR}/reactmicp/equilibrium_curve ) # static libraries # ---------------- if(SPECMICP_BUILD_STATIC) add_library(reactmicp_static STATIC ${REACTMICP_LIBFILES}) install(TARGETS reactmicp_static ARCHIVE DESTINATION ${STATIC_LIBRARY_INSTALL_DIR} ) else() add_library(reactmicp_static EXCLUDE_FROM_ALL STATIC ${REACTMICP_LIBFILES}) endif() set_target_properties(reactmicp_static PROPERTIES OUTPUT_NAME reactmicp) target_link_libraries(reactmicp_static dfpm) diff --git a/src/reactmicp/systems/saturated_react/kinetic_stagger.cpp b/src/reactmicp/systems/saturated_react/kinetic_stagger.cpp new file mode 100644 index 0000000..b6f2ebf --- /dev/null +++ b/src/reactmicp/systems/saturated_react/kinetic_stagger.cpp @@ -0,0 +1,200 @@ +/*------------------------------------------------------------------------------- + +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 "kinetic_stagger.hpp" +#include "variables.hpp" + +#include "../../../specmicp/adimensional/adimensional_system_solver.hpp" +#include "../../../specmicp/adimensional/adimensional_system_solution_extractor.hpp" +#include "../../../utils/log.hpp" +#include "../../../reactmicp/solver/staggers_base/stagger_structs.hpp" + +#include + +namespace specmicp { +namespace reactmicp { +namespace systems { +namespace satdiff { + + +//! \brief Initialize the stagger at the beginning of an iteration +void KineticStagger::initialize_timestep(scalar_t dt, VariablesBasePtr var) +{ + m_dt = dt; + + SaturatedVariablesPtr true_var = cast_var_from_base(var); + + // Initialize velocity using values from previous timestep + for (index_t node=0; nodenb_nodes(); ++node) + { + if (true_var->is_fixed_composition(node)) continue; + scalar_t alpha = 1.0; + for (index_t component: true_var->get_database()->range_aqueous_component()) + { + alpha = std::max(alpha, + 0.9*dt*true_var->solid_concentration(node, component, true_var->chemistry_rate()) + /(true_var->solid_concentration(node, component, true_var->displacement())) + ); + } + auto solid_velocity = true_var->velocity().segment( + true_var->offset_node(node)+true_var->offset_solid_concentration(), + true_var->nb_component()); + auto solid_chemistry_rate = true_var->chemistry_rate().segment( + true_var->offset_node(node)+true_var->offset_solid_concentration(), + true_var->nb_component()); + solid_velocity = 1/alpha * solid_chemistry_rate; + } + + m_model->initialize_timestep(dt, true_var->kinetic_variables()); + for (index_t node=0; nodenb_nodes(); ++node) + { + if (true_var->is_fixed_composition(node)) continue; + m_model->restart_timestep(node, true_var->equilibrium_solution(node), true_var->kinetic_variables()); + + for (index_t component=0; componentnb_component(); ++component) + { + + true_var->aqueous_concentration(node, component, true_var->chemistry_rate()) = - + m_model->get_velocity_kinetic(node, component, true_var->kinetic_variables()); + true_var->aqueous_concentration(node, component, true_var->velocity()) += ( + //true_var->aqueous_concentration(node, component, true_var->transport_rate()) + + true_var->aqueous_concentration(node, component, true_var->chemistry_rate()) + //- true_var->solid_concentration(node, component, true_var->chemistry_rate()) + ); + + + true_var->aqueous_concentration(node, component, true_var->displacement()) = + true_var->aqueous_concentration(node, component, true_var->predictor()) + + m_dt*true_var->aqueous_concentration(node, component, true_var->velocity()); + } + } +} + +//! \brief Solve the equation for the timestep +solver::StaggerReturnCode KineticStagger::restart_timestep(VariablesBasePtr var) +{ + SaturatedVariablesPtr true_var = cast_var_from_base(var); + int failed_chemistry = 0; +#ifdef SPECMICP_USE_OPENMP +#pragma omp parallel default(none) shared(true_var, failed_chemistry) + { +#pragma omp for schedule(dynamic, 5) + for (index_t node=0; nodenb_nodes(); ++node) + { + // only solve if necessary + if (true_var->is_fixed_composition(node) or failed_chemistry > 0) continue; + const auto retcode = solve_one_node(node, true_var); + if (retcode > 0) + { + ++failed_chemistry; + } + } + } +#else + { + for (index_t node=0; nodenb_nodes(); ++node) + { + if (true_var->is_fixed_composition(node)) continue; + const auto retcode = solve_one_node(node, true_var); + if (retcode > 0) + { + ++failed_chemistry; + break; + } + } + } +#endif // SPECMICP_USE_OPENMP + if (failed_chemistry > 0) + return solver::StaggerReturnCode::UnknownError; + return solver::StaggerReturnCode::ResidualMinimized; +} + +//! +int KineticStagger::solve_one_node(index_t node, SaturatedVariablesPtr true_var) +{ + + + + + AdimensionalSystemConstraints constraints(get_constraints(node)); + + constraints.total_concentrations = true_var->total_concentrations(node); + constraints.set_inert_volume_fraction(m_model->get_volume_fraction_kinetic(node, true_var->kinetic_variables())); + + AdimensionalSystemSolver adim_solver(true_var->get_database(), + constraints, + true_var->equilibrium_solution(node), + m_options); + Vector variables(true_var->equilibrium_solution(node).main_variables); + micpsolver::MiCPPerformance perf = adim_solver.solve(variables); + + micpsolver::MiCPSolverReturnCode retcode = perf.return_code; + if (retcode <= micpsolver::MiCPSolverReturnCode::NotConvergedYet) + { + + ERROR << "Failed to solve chemistry problem at node " << node + << ", return code = " << static_cast(retcode) + << ", residual = " << perf.current_residual; + ERROR << "Total concentration : \n" << constraints.total_concentrations; + return 1; + } + + true_var->equilibrium_solution(node) = adim_solver.get_raw_solution(variables); + AdimensionalSystemSolutionExtractor extractor(true_var->equilibrium_solution(node), + true_var->get_database(), + m_options.units_set); + + for (index_t component=0; componentnb_component(); ++component) + { + const scalar_t c_aq = extractor.density_water()*extractor.total_aqueous_concentration(component); + true_var->aqueous_concentration(node, component, true_var->displacement()) = c_aq; + + const scalar_t c_aq_0 = true_var->aqueous_concentration(node, component, true_var->predictor()); + const scalar_t vel_aq = (c_aq - c_aq_0)/m_dt; + true_var->aqueous_concentration(node, component, true_var->velocity()) = vel_aq; + + const scalar_t c_sol = extractor.total_immobile_concentration(component); + true_var->solid_concentration(node, component, true_var->displacement()) = c_sol; + + const scalar_t c_sol_0 = true_var->solid_concentration(node, component, true_var->predictor()); + const scalar_t vel_sol = (c_sol - c_sol_0)/m_dt; + true_var->solid_concentration(node, component, true_var->velocity()) = vel_sol; + true_var->solid_concentration(node, component, true_var->chemistry_rate()) = vel_sol; + } + + return 0; +} + +} // end namespace satdiff +} // end namespace systems +} // end namespace reactmicp +} // end namespace specmicp diff --git a/src/reactmicp/systems/saturated_react/kinetic_stagger.hpp b/src/reactmicp/systems/saturated_react/kinetic_stagger.hpp new file mode 100644 index 0000000..e92e1d0 --- /dev/null +++ b/src/reactmicp/systems/saturated_react/kinetic_stagger.hpp @@ -0,0 +1,137 @@ +/*------------------------------------------------------------------------------- + +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_REACTMICP_SATDIFF_KINETICSTAGGER_HPP +#define SPECMICP_REACTMICP_SATDIFF_KINETICSTAGGER_HPP + +//! \file saturated_react/kinetic_stagger.hpp +//! \brief The kinetic stagger for the saturated reactive transport solver + +#include "../../solver/staggers_base/chemistry_stagger_base.hpp" +#include "variablesfwd.hpp" + +#include "../../../specmicp/adimensional/adimensional_system_solver_structs.hpp" +#include "../../../specmicp/adimensional/adimensional_system_solution.hpp" + +#include + +namespace specmicp { +namespace reactmicp { +namespace systems { +namespace satdiff { + +using VariablesBasePtr = solver::VariablesBasePtr; + +//! \brief The Kinetic model +//! +//! This is an abstract class. +class KineticModelBase +{ +public: + virtual ~KineticModelBase() {} + //! \brief Initialize a timestep + virtual void initialize_timestep(scalar_t dt, KineticSaturatedVariablesBasePtr kinetic_variables) = 0; + //! \brief Restart a timestep + virtual void restart_timestep( + index_t node, + const AdimensionalSystemSolution& eq_solution, + KineticSaturatedVariablesBasePtr kinetic_variables) = 0; + //! \brief Return the volume fraction of kinetic (or inert) phases + virtual scalar_t get_volume_fraction_kinetic(index_t node, KineticSaturatedVariablesBasePtr kinetic_variables) = 0; + //! \brief Return the velocity of the total immobile kinetic concentration for component + virtual scalar_t get_velocity_kinetic(index_t node, index_t component, KineticSaturatedVariablesBasePtr kinetic_variables) = 0; +}; + +using KineticModelBasePtr = std::shared_ptr; + +//! \brief Solve the kinetic and equilibrium problems +class KineticStagger: public solver::ChemistryStaggerBase +{ +public: + KineticStagger( + index_t nb_nodes, + KineticModelBasePtr the_model, + AdimensionalSystemConstraints constraints, + AdimensionalSystemSolverOptions options + ): + m_model(the_model), + m_id_constraints(nb_nodes, 0), + m_list_constraints({constraints,}), + m_options(options) + {} + KineticStagger( + KineticModelBasePtr the_model, + std::vector list_constraints, + std::vector index_constraints, + AdimensionalSystemSolverOptions options + ): + m_model(the_model), + m_id_constraints(index_constraints), + m_list_constraints(list_constraints), + m_options(options) + {} + //! \brief Initialize the stagger at the beginning of the computation + void initialize(VariablesBasePtr var) override {} + + //! \brief Initialize the stagger at the beginning of an iteration + void initialize_timestep(scalar_t dt, VariablesBasePtr var) override; + + //! \brief Solve the equation for the timestep + solver::StaggerReturnCode restart_timestep(VariablesBasePtr var) override; + + //! \brief Solve the speciation problem at one node + int solve_one_node(index_t node, SaturatedVariablesPtr var); + + //! \brief Return the constraints for 'node' + //! + //! \param node Index of the node + AdimensionalSystemConstraints& get_constraints(index_t node) { + return m_list_constraints[m_id_constraints[node]]; + } + +private: + scalar_t m_dt {-1.0}; + + KineticModelBasePtr m_model; + + std::vector m_id_constraints; + std::vector m_list_constraints; + AdimensionalSystemSolverOptions m_options; +}; + +} // end namespace satdiff +} // end namespace systems +} // end namespace reactmicp +} // end namespace specmicp + + +#endif // SPECMICP_REACTMICP_SATDIFF_KINETICSTAGGER_HPP diff --git a/src/reactmicp/systems/saturated_react/transport_program.cpp b/src/reactmicp/systems/saturated_react/transport_program.cpp index 01182ed..dcf1390 100644 --- a/src/reactmicp/systems/saturated_react/transport_program.cpp +++ b/src/reactmicp/systems/saturated_react/transport_program.cpp @@ -1,317 +1,321 @@ /*------------------------------------------------------------------------------- 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 "transport_program.hpp" #include "../../../dfpm/meshes/mesh1d.hpp" #include "variables.hpp" namespace specmicp { namespace reactmicp { namespace systems { namespace satdiff { //class SaturatedDiffusion:: SaturatedDiffusion::SaturatedDiffusion(SaturatedVariablesPtr variables, std::vector list_fixed_nodes): m_ndf(2*variables->nb_component()), m_tot_ndf(2*variables->nb_component()*variables->get_mesh()->nb_nodes()), m_mesh(variables->get_mesh()), m_variables(variables), m_is_in_residual_computation(false) { number_equations(list_fixed_nodes, {}, {0, 1}); } SaturatedDiffusion::SaturatedDiffusion(SaturatedVariablesPtr variables, std::vector list_fixed_nodes, std::map list_slave_nodes, std::vector list_immobile_components): m_ndf(2*variables->nb_component()), m_tot_ndf(2*variables->nb_component()*variables->get_mesh()->nb_nodes()), m_mesh(variables->get_mesh()), m_variables(variables), m_is_in_residual_computation(false) { number_equations(list_fixed_nodes, list_slave_nodes, list_immobile_components); } void SaturatedDiffusion::number_equations(std::vector list_fixed_nodes, std::map list_slave_nodes, std::vector list_immobile_components) { m_ideq.resizeLike(m_variables->displacement()); m_ideq.setZero(); // flag fixed nodes for (index_t node: list_fixed_nodes) { for (index_t component=2; componentnb_component(); ++component) { m_ideq(m_variables->dof_aqueous_concentration(node, component)) = no_equation; } } // flag slaves nodes // we flag them by making their ideq more negative than no_equation for (auto slave_pair: list_slave_nodes) { for (index_t component=2; componentnb_component(); ++component) { m_ideq(m_variables->dof_aqueous_concentration(slave_pair.first, component)) = no_equation-1; } } // set equation numbers index_t neq = 0; for (index_t node=0; nodenb_nodes(); ++node) { for (index_t component: list_immobile_components) { m_ideq(m_variables->dof_aqueous_concentration(node, component)) = no_equation; m_ideq(m_variables->dof_aqueous_concentration(node, component)) = no_equation; } for (index_t component=0; componentnb_component(); ++component) { index_t dof = m_variables->dof_aqueous_concentration(node, component); if (m_ideq(dof) > no_equation) // attribute an equation number if it is NOT a slave nor a fixed node { m_ideq(dof) = neq; ++neq; } dof = m_variables->dof_solid_concentration(node, component); m_ideq(dof) = no_equation; } } // slave nodes // attribute the correct equation number for (auto slave_pair: list_slave_nodes) { for (index_t component=1; componentnb_component(); ++component) { m_ideq(m_variables->dof_aqueous_concentration(slave_pair.first, component)) = m_ideq(m_variables->dof_aqueous_concentration(slave_pair.second, component)); } } m_neq = neq; } void SaturatedDiffusion::compute_residuals( const Vector& displacement, const Vector& velocity, Vector& residual) { residual = Vector::Zero(get_neq()); m_is_in_residual_computation = true; for (index_t element: m_mesh->range_elements()) { for (index_t component=2; componentnb_component(); ++component) { Eigen::Vector2d element_residual; element_residual.setZero(); residuals_element_component(element, component, displacement, velocity, element_residual); for (index_t en=0; en<2; ++en) { const index_t node = m_mesh->get_node(element, en); const index_t id = m_ideq(m_variables->dof_aqueous_concentration(node, component)); if (id != no_equation) {residual(id) += element_residual(en);} } } } m_is_in_residual_computation = false; } void SaturatedDiffusion::residuals_element_component( index_t element, index_t component, const Vector& displacement, const Vector& velocity, Eigen::Vector2d& element_residual ) { const scalar_t mass_coeff_0 = -m_mesh->get_volume_cell_element(element, 0); const scalar_t mass_coeff_1 = -m_mesh->get_volume_cell_element(element, 1); const index_t node_0 = m_mesh->get_node(element, 0); const index_t node_1 = m_mesh->get_node(element, 1); const scalar_t diff_coeff = 1.0/(0.5/m_variables->diffusion_coefficient(node_0) + 0.5/m_variables->diffusion_coefficient(node_1)); scalar_t flux_coeff = -( m_mesh->get_face_area(element) / m_mesh->get_dx(element) * diff_coeff ) ; const index_t dof_0 = m_variables->dof_aqueous_concentration(node_0, component); const index_t dof_1 = m_variables->dof_aqueous_concentration(node_1, component); // diffusion scalar_t flux_diffusion = flux_coeff*(displacement(dof_0) - displacement(dof_1)); element_residual(0) = flux_diffusion; element_residual(1) = - flux_diffusion; // advection if (m_variables->fluid_velocity(element) != 0.0) { scalar_t flux_advection = (m_mesh->get_face_area(element)) *m_variables->fluid_velocity(element); if (m_variables->fluid_velocity(element) > 0) { flux_advection *= (displacement(dof_0) - displacement(dof_1)); element_residual(1) += flux_advection; } else { flux_advection *= (displacement(dof_1) - displacement(dof_0)); element_residual(0) -= flux_advection; } } if (m_is_in_residual_computation) { m_variables->aqueous_concentration(node_0, component, m_variables->transport_rate()) += element_residual(0); m_variables->aqueous_concentration(node_1, component, m_variables->transport_rate()) += element_residual(1); } // velocity element_residual(0) += mass_coeff_0*(velocity(dof_0)*m_variables->porosity(node_0) +m_variables->vel_porosity(node_0)*displacement(dof_0)); element_residual(1) += mass_coeff_1*(velocity(dof_1)*m_variables->porosity(node_1) + m_variables->vel_porosity(node_1)*displacement(dof_1)); // external rate - element_residual(0) += mass_coeff_0*m_variables->solid_concentration(node_0, component, - m_variables->chemistry_rate()); - element_residual(1) += mass_coeff_1*m_variables->solid_concentration(node_1, component, - m_variables->chemistry_rate()); + element_residual(0) -= mass_coeff_0*( + - m_variables->solid_concentration(node_0, component,m_variables->chemistry_rate()) + + m_variables->aqueous_concentration(node_0, component,m_variables->chemistry_rate()) + ); + element_residual(1) -= mass_coeff_1*( + - m_variables->solid_concentration(node_1, component,m_variables->chemistry_rate()) + + m_variables->aqueous_concentration(node_1, component,m_variables->chemistry_rate()) + ); } void SaturatedDiffusion::compute_jacobian( Vector& displacement, Vector& velocity, Eigen::SparseMatrix& jacobian, scalar_t alphadt ) { dfpm::list_triplet_t jacob; const index_t ncomp = m_variables->nb_component(); const index_t estimation = m_mesh->nb_nodes()*(ncomp*m_mesh->nen); jacob.reserve(estimation); for (index_t element: m_mesh->range_elements()) { jacobian_element(element, displacement, velocity, jacob, alphadt); } jacobian = Eigen::SparseMatrix(get_neq(), get_neq()); jacobian.setFromTriplets(jacob.begin(), jacob.end()); } void SaturatedDiffusion::jacobian_element( index_t element, Vector& displacement, Vector& velocity, dfpm::list_triplet_t& jacobian, scalar_t alphadt) { for (index_t component=1; componentnb_component(); ++component) { Eigen::Vector2d element_residual_orig(Eigen::Vector2d::Zero()); residuals_element_component(element, component, displacement, velocity, element_residual_orig); for (index_t en=0; en<2; ++en) { Eigen::Vector2d element_residual(Eigen::Vector2d::Zero()); const index_t node = m_mesh->get_node(element, en); const index_t dof = m_variables->dof_aqueous_concentration(node, component); const index_t idc = m_ideq(dof); if (idc == no_equation) continue; const scalar_t tmp_v = velocity(dof); const scalar_t tmp_d = displacement(dof); scalar_t h = eps_jacobian*std::abs(tmp_v); if (h < 1e-4*eps_jacobian) h = eps_jacobian; velocity(dof) = tmp_v + h; h = velocity(dof) - tmp_v; displacement(dof) = tmp_d + alphadt*h; residuals_element_component(element, component, displacement, velocity, element_residual); velocity(dof) = tmp_v; displacement(dof) = tmp_d; for (index_t enr=0; enr<2; ++enr) { const index_t noder = m_mesh->get_node(element, enr); const index_t idr = m_ideq(m_variables->dof_aqueous_concentration(noder, component)); if (idr == no_equation) continue; jacobian.push_back(dfpm::triplet_t( idr, idc, (element_residual(enr) - element_residual_orig(enr))/h )); } } } } //! \brief Update the solutions void SaturatedDiffusion::update_solution(const Vector& update, scalar_t lambda, scalar_t alpha_dt, Vector& predictor, Vector& displacement, Vector& velocity) { for (index_t node: m_mesh->range_nodes()) { for (index_t component=1; component< m_variables->nb_component(); ++component) { const index_t dof = m_variables->dof_aqueous_concentration(node, component); const index_t id = m_ideq(dof); if (id == no_equation) continue; velocity(dof) += lambda*update(id); } } //displacement = m_variables->predictor() + alpha_dt*velocity; displacement = predictor + alpha_dt*velocity; } } // end namespace satdiff } // end namespace systems } // end namespace reactmicp } // end namespace specmicp diff --git a/src/reactmicp/systems/saturated_react/variables.hpp b/src/reactmicp/systems/saturated_react/variables.hpp index 1003aa4..b7a5c11 100644 --- a/src/reactmicp/systems/saturated_react/variables.hpp +++ b/src/reactmicp/systems/saturated_react/variables.hpp @@ -1,262 +1,280 @@ /*------------------------------------------------------------------------------- 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_REACTMICP_SYSTEMS_SATURATED_VARIABLES_HPP #define SPECMICP_REACTMICP_SYSTEMS_SATURATED_VARIABLES_HPP #include "../../../types.hpp" #include "../../../database.hpp" #include "../../solver/staggers_base/variables_base.hpp" #include "../../../specmicp/adimensional/adimensional_system_solution.hpp" #include // forward declaration // =================== #include "../../../dfpm/meshes/mesh1dfwd.hpp" namespace specmicp { namespace reactmicp { namespace solver { using VariablesBasePtr = std::shared_ptr; } namespace systems { namespace satdiff { class SaturatedVariablesFactory; } // end namespace satdiff } // end namespace systems } // end namespace reactmicp } // end namespace specmicp // Class declaration // ================= namespace specmicp { namespace reactmicp { namespace systems { namespace satdiff { +//! \brief Base class for the kinetic variables +class KineticSaturatedVariablesBase +{ +}; + +//! \brief The kinetic variables are +using KineticSaturatedVariablesBasePtr = std::shared_ptr; + +class DummyKineticSaturatedVariables: public KineticSaturatedVariablesBase +{}; + //! \brief Variables for the saturated reactive transport system //! //! Contain all the variables that need to be shared between the staggers class SaturatedVariables: public solver::VariablesBase { // SaturatedVariablesFactory should be the class to use to inialize // the variables correctly friend class SaturatedVariablesFactory; public: SaturatedVariables(mesh::Mesh1DPtr the_mesh, RawDatabasePtr the_database); //! \brief Return the mesh mesh::Mesh1DPtr get_mesh() {return m_mesh;} //! \brief Return the database RawDatabasePtr get_database() {return m_database;} //! \brief Return the number of components index_t nb_component() {return m_database->nb_component();} //! \brief Return the number of nodes index_t nb_nodes() {return m_is_fixed_composition.size();} //! \brief Return true if 'node' has a fixed composition index_t is_fixed_composition(index_t node) {return m_is_fixed_composition[node];} // Main variables // ============== //! \brief Return the main variable vector Vector& displacement() {return m_displacement;} //! \brief Return the main variable vector at the beginning of the timestep Vector& predictor() {return m_predictor;} //! \brief Return the velocity of the main variables Vector& velocity() {return m_velocity;} //! \brief Return the rate of change of the main variables due to the transport operator Vector& transport_rate() {return m_transport_rate;} //! \brief Return the rate of change of the main variables due to the chemistry operator Vector& chemistry_rate() {return m_chemistry_rate;} // Access to main variables // ======================== //! \brief Return the number of degree of freedom (per node) in the main variables vector index_t ndf() {return 2*m_database->nb_component();} //! \brief Return the offset of 'node' in the main variables vector index_t offset_node(index_t node) {return node*ndf();} //! \brief Return the offset of the aqueous concentration variables in the main variables vector index_t offset_aqueous_concentration() {return 0;} //! \brief Return the offset of the aqueous concentrations variables in the main variables vector index_t offset_aqueous_concentration(index_t node) { return offset_aqueous_concentration()+offset_node(node);} //! \brief Return the offset of the solid concentration variables in the main variables vector index_t offset_solid_concentration() {return m_database->nb_component();} //! \brief Return the offset of the solid concentrations variables in the main variables vector index_t offset_solid_concentration(index_t node) { return offset_solid_concentration()+offset_node(node);} //! \brief Return the degree of freedom number for the aqueous concentration of 'component' at 'node' index_t dof_aqueous_concentration(index_t node, index_t component) { return (component + offset_aqueous_concentration(node)); } //! \brief Return the degree of freedom number for the solid concentration of 'component' at 'node' index_t dof_solid_concentration(index_t node, index_t component) { return (component + offset_solid_concentration(node)); } //! \brief Return the aqueous concentration of 'component' at 'node' in 'var' //! //! 'var' is any of the main variables vector, it may be a velocity vector scalar_t& aqueous_concentration(index_t node, index_t component, Vector& var) { return var(dof_aqueous_concentration(node, component)); } //! \brief Return the aqueous concentration of 'component' at 'node' in main variables //! //! 'var' is any of the main variables vector, it may be a velocity vector scalar_t& aqueous_concentration(index_t node, index_t component) { return m_displacement(dof_aqueous_concentration(node, component)); } //! \brief Return the solid concentration of 'component' at 'node' in 'var' //! //! 'var' is any of the main variables vector, it may be a velocity vector scalar_t& solid_concentration(index_t node, index_t component, Vector& var){ return var(dof_solid_concentration(node, component)); } //! \brief Return the solid concentration of 'component' at 'node' in main variables //! //! 'var' is any of the main variables vector, it may be a velocity vector scalar_t& solid_concentration(index_t node, index_t component){ return m_displacement(dof_solid_concentration(node, component)); } //! \brief Return a vector containing the total concentrations computed from the main variables //! //! This is to be used to restart the chemistry computation Vector total_concentrations(index_t node); // Equilibrium // =========== //! \brief Return the solution of the speciation solver at 'node' AdimensionalSystemSolution& equilibrium_solution(index_t node) { return m_equilibrium_solutions[node]; } //! \brief Return a reference to the equilibrium solutions solution const std::vector& equilibrium_solutions() { return m_equilibrium_solutions; } // Upscaling // ========= //! \brief Return the offset for 'node' in the upscaling variables vector index_t offset_node_upscaling(index_t node) {return ndf_upscaling()*node;} //! \brief Return the number fo degree of freedom (per node) for the upscaling vector index_t ndf_upscaling() {return 5;} //! \brief Return the degree of freedom for the porosity at 'node' index_t dof_porosity(index_t node) {return 0 + offset_node_upscaling(node);} //! \brief Return the degree of freedom of the porosity velocity at 'node' index_t dof_vel_porosity(index_t node) {return 1 + offset_node_upscaling(node);} //! \brief Return the degree of freedom of the diffusion coefficient at 'node' index_t dof_diffusion_coefficient(index_t node) {return 2 + offset_node_upscaling(node);} //! \brief Return the degree of freedom of the permeability at 'node' index_t dof_permeability(index_t node) {return 3 + offset_node_upscaling(node);} //! \brief Return the fluid velocity index_t dof_fluid_velocity(index_t node) {return 4 + offset_node_upscaling(node);} //! \brief Return the porosity at 'node' scalar_t& porosity(index_t node) {return m_upscaling(dof_porosity(node));} //! \brief Return the rate of change of the porosity at 'node' scalar_t& vel_porosity(index_t node) {return m_upscaling(dof_vel_porosity(node));} //! \brief Return the diffusion coefficient at 'node' scalar_t& diffusion_coefficient(index_t node) {return m_upscaling(dof_diffusion_coefficient(node));} //! \brief Return the permeability at 'node' scalar_t& permeability(index_t node) {return m_upscaling(dof_permeability(node));} //! \brief Return the fluid velocity at 'node' scalar_t& fluid_velocity(index_t node) {return m_upscaling(dof_fluid_velocity(node));} //! \brief Return the vector of upscaling variables Vector& upscaling_variables() {return m_upscaling;} + KineticSaturatedVariablesBasePtr& kinetic_variables() {return m_kinetic_variables;} + + //! \brief Reset the main variables void reset_main_variables() override; private: // ############ // // Attributes // // ############ // mesh::Mesh1DPtr m_mesh; RawDatabasePtr m_database; std::vector m_is_fixed_composition; // Main variables // ============== Vector m_displacement; Vector m_predictor; Vector m_velocity; Vector m_transport_rate; Vector m_chemistry_rate; // Equilibrium // =========== std::vector m_equilibrium_solutions; + // Kinetic + // ======= + KineticSaturatedVariablesBasePtr m_kinetic_variables; + // Upscaling // ========= Vector m_upscaling; }; //! \brief typedef of a shared pointer of a SaturatedVariables using SaturatedVariablesPtr = std::shared_ptr; // Casting function // ================= //! \brief Static cast to a SaturatedVariablesPtr inline SaturatedVariablesPtr cast_var_from_base(solver::VariablesBasePtr var) { return std::static_pointer_cast(var); } //! \brief Static cast from a SaturatedVariablesPtr inline solver::VariablesBasePtr cast_var_to_base(SaturatedVariablesPtr var) { return std::static_pointer_cast(var); } } // end namespace satdiff } // end namespace systems } // end namespace reactmicp } // end namespace specmicp #endif // SPECMICP_REACTMICP_SYSTEMS_SATURATED_VARIABLES_HPP diff --git a/src/reactmicp/systems/saturated_react/variablesfwd.hpp b/src/reactmicp/systems/saturated_react/variablesfwd.hpp index 7cf95b9..be33d62 100644 --- a/src/reactmicp/systems/saturated_react/variablesfwd.hpp +++ b/src/reactmicp/systems/saturated_react/variablesfwd.hpp @@ -1,51 +1,57 @@ /*------------------------------------------------------------------------------- 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_REACTMICP_SYSTEMS_SATURATED_VARIABLESFWD_HPP #define SPECMICP_REACTMICP_SYSTEMS_SATURATED_VARIABLESFWD_HPP #include namespace specmicp { namespace reactmicp { namespace systems { namespace satdiff { class SaturatedVariables; using SaturatedVariablesPtr = std::shared_ptr; +//! \brief Base class for the kinetic variables +class KineticSaturatedVariablesBase; + +//! \brief The kinetic variables are +using KineticSaturatedVariablesBasePtr = std::shared_ptr; + } // end namespace satdiff } // end namespace systems } // end namespace reactmicp } // end namespace specmicp #endif // SPECMICP_REACTMICP_SYSTEMS_SATURATED_VARIABLESFWD_HPP