diff --git a/examples/reactmicp/leaching.cpp b/examples/reactmicp/leaching.cpp index 3267816..9c68a17 100644 --- a/examples/reactmicp/leaching.cpp +++ b/examples/reactmicp/leaching.cpp @@ -1,293 +1,309 @@ #include #include "database/database.hpp" #include "specmicp/reaction_path.hpp" #include "reactmicp/mesh.hpp" #include "reactmicp/systems/saturated_diffusion/reactive_transport_solver.hpp" +#include "reactmicp/systems/saturated_diffusion/transport_parameters.hpp" +#include "reactmicp/systems/saturated_diffusion/boundary_conditions.hpp" #include #include using namespace specmicp; using namespace specmicp::reactmicp; using namespace specmicp::reactmicp::systems::siasaturated; const double M_CaO = 56.08; const double M_SiO2 = 60.09; const double M_Al2O3 = 101.96; const double M_SO3 = 80.06; void compo_from_oxyde(Vector& compo_oxyde, Vector& compo_species) { Eigen::MatrixXd compo_in_oxyde(4, 4); Vector molar_mass(4); molar_mass << 1.0/M_CaO, 1.0/M_SiO2, 1.0/M_Al2O3, 1.0/M_SO3; compo_oxyde = compo_oxyde.cwiseProduct(molar_mass); //C3S C2S C3A Gypsum compo_in_oxyde << 3, 2, 3, 1, // CaO 1, 1, 0, 0, // SiO2 0, 0, 1, 0, // Al2O3 0, 0, 0, 1; // SO3 Eigen::FullPivLU solver(compo_in_oxyde); compo_species = solver.solve(compo_oxyde); } std::shared_ptr set_database() { specmicp::database::Database database("../data/cemdata_specmicp.js"); std::map swapping ({ {"H[+]","HO[-]"}, {"Si(OH)4", "SiO(OH)3[-]"}, {"Al[3+]","Al(OH)4[-]"} }); database.swap_components(swapping); //std::vector to_keep = {"H2O", "HO[-]", "Ca[2+]", "SiO(OH)3[-]"}; //database.keep_only_components(to_keep); return database.get_database(); } EquilibriumState get_initial_state( std::shared_ptr database, Vector& compo, double mult) { std::shared_ptr model = std::make_shared(); compo = mult*compo; const double m_c3s = compo(0); const double m_c2s = compo(1); const double m_c3a = compo(2); const double m_gypsum = compo(3); const double wc = 1.0; const double m_water = wc*( (3*M_CaO+M_SiO2)*m_c3s + (2*M_CaO+M_SiO2)*m_c2s + (3*M_CaO+M_Al2O3)*m_c3a + (M_CaO + M_SO3)*m_gypsum )*1e-3; model->amount_aqueous = { {"H2O", specmicp::reaction_amount_t(m_water/database->molar_mass_basis_si(0), 0.0)}, {"HO[-]", specmicp::reaction_amount_t(0.0, 0.0)}, }; model->amount_minerals = { {"C3S", specmicp::reaction_amount_t(m_c3s, 0.0)}, {"C2S", specmicp::reaction_amount_t(m_c2s, 0.0)}, {"C3A", specmicp::reaction_amount_t(m_c3a, 0.0)}, {"Gypsum", specmicp::reaction_amount_t(m_gypsum, 0.0)}, }; model->minerals_to_keep = { "Portlandite", "CSH,jennite", "CSH,tobermorite", "SiO2,am", "Al(OH)3,am", "Monosulfoaluminate", "Straetlingite", "Gypsum", "Ettringite" }; Eigen::VectorXd x; specmicp::ReactionPathDriver driver(model, database, x); //driver.get_options().solver_options.penalization_factor =1.0; driver.get_options().ncp_function = specmicp::micpsolver::NCPfunction::penalizedFB; //driver.get_options().ncp_function = specmicp::micpsolver::NCPfunction::min; driver.get_options().solver_options.use_scaling = false; driver.get_options().solver_options.max_factorization_step = 4; driver.get_options().solver_options.factor_gradient_search_direction = 200; driver.get_options().solver_options.maxstep = 50; driver.get_options().solver_options.maxiter_maxstep = 50; driver.get_options().solver_options.max_iter = 100; driver.get_options().allow_restart = true; specmicp::micpsolver::MiCPPerformance perf = driver.one_step(x); if (perf.return_code != specmicp::micpsolver::MiCPSolverReturnCode::ResidualMinimized) throw std::runtime_error("Error - system cannot be solved"); return driver.get_current_solution(); } EquilibriumState get_blank_state(std::shared_ptr database) { std::shared_ptr model = std::make_shared(); const double m_c3s = 1e-7; const double m_c2s = 1e-7; const double m_c3a = 1e-7; const double m_gypsum = 1e-7; const double m_water = 1; model->amount_aqueous = { {"H2O", specmicp::reaction_amount_t(m_water/database->molar_mass_basis_si(0), 0.0)}, - {"HO[-]", specmicp::reaction_amount_t(0.0, 0.0)}, + {"HO[-]", specmicp::reaction_amount_t(-0.0, 0.0)}, + {"SiO(OH)3[-]", specmicp::reaction_amount_t(0.0, 0.0)} }; model->amount_minerals = { {"C3S", specmicp::reaction_amount_t(m_c3s, 0.0)}, {"C2S", specmicp::reaction_amount_t(m_c2s, 0.0)}, {"C3A", specmicp::reaction_amount_t(m_c3a, 0.0)}, {"Gypsum", specmicp::reaction_amount_t(m_gypsum, 0.0)}, + {"SiO2,am", specmicp::reaction_amount_t(0.0001, 0.0)} }; Eigen::VectorXd x; specmicp::ReactionPathDriver driver(model, database, x); //driver.get_options().solver_options.penalization_factor =1.0; driver.get_options().ncp_function = specmicp::micpsolver::NCPfunction::penalizedFB; //driver.get_options().ncp_function = specmicp::micpsolver::NCPfunction::min; driver.get_options().solver_options.use_scaling = false; driver.get_options().solver_options.max_factorization_step = 1; driver.get_options().solver_options.factor_gradient_search_direction = 50; driver.get_options().solver_options.maxstep = 50; driver.get_options().solver_options.maxiter_maxstep = 50; driver.get_options().solver_options.max_iter = 100; driver.get_options().allow_restart = true; specmicp::micpsolver::MiCPPerformance perf = driver.one_step(x); if (perf.return_code != specmicp::micpsolver::MiCPSolverReturnCode::ResidualMinimized) throw std::runtime_error("Error - blank system cannot be solved"); return driver.get_current_solution(); } double get_dt( std::shared_ptr parameters, std::shared_ptr themesh ) { //const scalar_t max_d = parameters->diffusion_coefficient(1); //scalar_t dt = std::pow(themesh->get_dx(0), 2)/max_d; //if (dt<3600*24) dt=3600*24; //if (dt>2*3600*24) dt=2*3600*24; - return 3600*24; + return 3600*24/5; } int main() { - specmicp::stdlog::ReportLevel() = specmicp::logger::Error; + specmicp::stdlog::ReportLevel() = specmicp::logger::Warning; specmicp::logger::ErrFile::stream() = &std::cerr; Vector compo_oxyde(4), compo_species(4); compo_oxyde << 62.9, 20.6, 5.8, 3.1; compo_from_oxyde(compo_oxyde, compo_species); std::shared_ptr thedatabase = set_database(); - scalar_t radius = 7; //cm + scalar_t radius = 3.5; //cm scalar_t height = 8; //cm - index_t nb_nodes = 101; + index_t nb_nodes = 36; mesh::Mesh1DPtr themesh = mesh::axisymmetric_uniform_mesh1d(nb_nodes, radius, height); std::shared_ptr parameters = std::make_shared(nb_nodes, 1e4*std::exp(9.95*0.25-29.08), 0.25); - - EquilibriumState initial_state = get_initial_state(thedatabase, compo_species, 0.2); + parameters->diffusion_coefficient(0) = 2.2e-5; + parameters->porosity(0) = 1.0; + EquilibriumState initial_state = get_initial_state(thedatabase, compo_species, 0.5); SIABoundaryConditions bcs(nb_nodes); bcs.list_initial_states.push_back(initial_state); bcs.list_initial_states.push_back(get_blank_state(thedatabase)); bcs.bs_types[0] = BoundaryConditionType::FixedComposition; bcs.initial_states[0] = 1; database::Database dbhandler(thedatabase); index_t id_ca = dbhandler.component_label_to_id("Ca[2+]"); index_t id_si = dbhandler.component_label_to_id("SiO(OH)3[-]"); std::ofstream output_ca, output_aqca, output_si, output_poro; output_ca.open("out_ca.dat"); output_aqca.open("out_aqca.dat"); output_si.open("out_si.dat"); output_poro.open("out_poro.ca"); std::vector mineral_out(thedatabase->nb_mineral); for (int mineral: thedatabase->range_mineral()) { mineral_out[mineral].open("out_mineral_"+std::to_string(mineral)+".dat"); } SIASaturatedReactiveTransportSolver solver(themesh, thedatabase, parameters); solver.apply_boundary_conditions(bcs); solver.get_variables().scale_volume(); solver.use_sia(100); - solver.get_options().residual_tolerance = 1e-4; + solver.get_options().residual_tolerance = 1e-3; solver.get_options().threshold_update_disable_speciation = 1e-16; solver.get_options().threshold_update_transport = 1e-10; + // Transport solver.get_options().transport_options.residuals_tolerance = 1e-5; + // Speciation + solver.get_options().speciation_options.solver_options.use_crashing = false; + solver.get_options().speciation_options.solver_options.use_scaling = false; + solver.get_options().speciation_options.solver_options.max_iter = 100; + solver.get_options().speciation_options.solver_options.maxstep = 50; + + std::cout << "id SiO2,am : " << dbhandler.mineral_label_to_id("SiO2,am") << std::endl; + //thedatabase->stability_mineral( + // dbhandler.mineral_label_to_id("SiO2,am")) = database::MineralStabilityClass::cannot_dissolve; + //double dt = 5e4; double total=0 ; uindex_t k = 1; - std::cout << solver.get_variables().total_mobile_concentrations() << std::endl; - while (total < 365*24*3600) { double dt = get_dt(parameters, themesh); total += dt; std::cout << " iterations : " << k << " - dt : " << dt << std::endl; SIASaturatedReactiveTransportSolverReturnCode retcode = solver.solve_timestep(dt); if (retcode != SIASaturatedReactiveTransportSolverReturnCode::Success) throw std::runtime_error("Failed to converge : "+std::to_string((int) retcode)); SIASaturatedVariables& variables = solver.get_variables(); //std::cout << variables.total_mobile_concentrations() << std::endl; // compute porosity - for (index_t node=1; noderange_nodes()) { double volume = variables.volume_minerals(node); double cell_volume = themesh->get_volume_cell(node); specmicp_assert(volume < cell_volume); double porosity = (cell_volume - volume*1e6) / cell_volume; parameters->porosity(node) = porosity; - parameters->diffusion_coefficient(node) = porosity>0.92?2.2e-5:1e4*std::exp(9.95*porosity-29.08); + parameters->diffusion_coefficient(node) = + porosity>0.6?1e4*std::exp(9.95*0.6-29.08):1e4*std::exp(9.95*porosity-29.08); //std::cout << parameters->diffusion_coefficient(node); } + parameters->diffusion_coefficient(0) = parameters->diffusion_coefficient(1); //std::cout << variables.speciation_variables() << std::endl; if (k < 10 or k % 5 == 0) { output_ca << total << " " << total/(3600.0*24.0); output_aqca << total << " " << total/(3600.0*24.0); output_si << total << " " << total/(3600.0*24.0); output_poro << total << " " << total/(3600.0*24.0); for(index_t node: themesh->range_nodes()) { output_ca << " " << variables.nodal_component_update_total_amount(node, id_ca); output_aqca << " " << variables.total_mobile_concentration(node, id_ca); output_si << " " << variables.nodal_component_update_total_amount(node, id_si); double volume = variables.volume_minerals(node); double cell_volume = themesh->get_volume_cell(node); output_poro << " " << (cell_volume - volume*1e6) / cell_volume; } output_ca << std::endl; output_aqca << std::endl; output_si << std::endl; output_poro << std::endl; for (index_t mineral: thedatabase->range_mineral()) { mineral_out[mineral] << total << " " << total/(3600.0*24.0); for(index_t node: themesh->range_nodes()) { mineral_out[mineral] << " " << variables.mineral_amount(node, mineral); } mineral_out[mineral] << std::endl; } } ++k; } output_ca.close(); output_aqca.close(); output_si.close(); output_poro.close(); for (index_t mineral: thedatabase->range_mineral()) { mineral_out[mineral].close(); } } diff --git a/src/reactmicp/systems/diffusion/diffusion.cpp b/src/reactmicp/systems/diffusion/diffusion.cpp index 3e70f44..b0ec6de 100644 --- a/src/reactmicp/systems/diffusion/diffusion.cpp +++ b/src/reactmicp/systems/diffusion/diffusion.cpp @@ -1,623 +1,623 @@ /*------------------------------------------------------- - Module : reactmicp/systems/diffusion - File : diffusion.cpp - Author : Fabien Georget Copyright (c) 2014, Fabien 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: * Redistributions of source code must retain the above copyright notice, this list of conditions and the following disclaimer. * 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. * Neither the name of the Princeton University 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 OWNER 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 "diffusion.hpp" #include "physics/laws.hpp" #include "utils/log.hpp" #include #define EPS_J 1e-8 namespace specmicp { namespace reactmicp { namespace systems { DiffusionProgram::DiffusionProgram(std::shared_ptr themesh, std::shared_ptr thedatabase, std::shared_ptr parameters, diffusion::ListBoundaryCondition& list_bcs, const EquilibriumState& initialstate, Variables& var): DiffusionProgram(themesh, thedatabase, parameters, list_bcs, var) { initialize_no_bc_with_init(initialstate, var); } // ================================== // // // // Residuals // // // // ================================== // void DiffusionProgram::get_residuals(const Variables& variable, Vector& residual) { m_is_micp = std::vector(get_neq(), false); residual = Eigen::VectorXd::Zero(get_neq()); residual_transport(variable, residual); residual_speciation(variable, residual); } // Transport // ========= void DiffusionProgram::residual_transport(const Variables& variable, Vector& residual) { for (ind_t element: m_mesh->range_elements()) { element_residual_transport(element, variable, residual); } } void DiffusionProgram::element_residual_transport_component( ind_t element, int component, const Variables& variable, Vector& element_residual) { Eigen::Matrix2d mass, jacob; Eigen::Vector2d velocity, displacement; double mass_coeff = -(m_param->density_water() * m_param->porosity(0) * - m_mesh->crosssection()*m_mesh->get_dx(element)/2); + m_mesh->get_volume_element(element)/2); mass << 1, 0, 0, 1; mass *= mass_coeff; - double flux_coeff = -( m_mesh->crosssection() / m_mesh->get_dx(element) + double flux_coeff = -( m_mesh->get_face_area(element) / m_mesh->get_dx(element) * m_param->porosity(0) * m_param->diffusion_coefficient(element) * m_param->density_water() ); jacob << 1, -1, -1, 1; jacob *= flux_coeff; velocity << variable.velocity(m_ideq.get_dof_diffusion(m_mesh->get_node(element,0), component)), variable.velocity(m_ideq.get_dof_diffusion(m_mesh->get_node(element,1), component)); displacement << m_ideq.mobile_total_concentration(m_mesh->get_node(element,0), component, variable), m_ideq.mobile_total_concentration(m_mesh->get_node(element,1), component, variable); element_residual = mass*velocity+jacob*displacement; for (int en: m_mesh->range_nen()) { - element_residual(en) += m_mesh->crosssection()*m_mesh->get_dx(element) * + element_residual(en) += m_mesh->get_volume_element(element) * nodal_mineral_transient_term_transport( m_mesh->get_node(element, en), component, variable)/2.0; } } void DiffusionProgram::element_residual_transport(ind_t element, const Variables& variable, Vector& residual) { Eigen::VectorXd element_residual(2); for (ind_t component: m_data->range_aqueous_component()) { element_residual_transport_component(element, component, variable, element_residual); for (int en: m_mesh->range_nen()) { const ind_t node = m_mesh->get_node(element, en); const ind_t id = m_ideq.id_equation_diffusion(node, component); if (id != no_equation) {residual.coeffRef(id) += element_residual(en);} } } } double DiffusionProgram::nodal_mineral_transient_term_transport(ind_t node, ind_t component, const Variables& variable) { double mineral_term = 0; for (ind_t mineral: m_data->range_mineral()) { if (m_data->nu_mineral(mineral, component) == 0 ) continue; mineral_term -= ( m_data->nu_mineral(mineral, component) * variable.velocity(m_ideq.get_dof_mineral(node, mineral)) ); } return mineral_term; } // Speciation // ========== void DiffusionProgram::residual_speciation(const Variables& variable, Vector& residual) { for (ind_t node: m_mesh->range_nodes()) { nodal_residual_massbalance(node, variable, residual); nodal_residual_mineral(node, variable, residual); } } void DiffusionProgram::nodal_residual_massbalance(ind_t node, const Variables& variable, Vector& residual) { for (ind_t component: m_data->range_aqueous_component()) { ind_t id = m_ideq.id_equation_massbalance(node, component); if (id != no_equation) { residual(id) += nodal_component_residual_massbalance(node, component, variable); } } } double DiffusionProgram::nodal_component_residual_massbalance(ind_t node, int component, const Variables& variable) { double sum = pow10(m_ideq.component_concentration(node, component, variable)); for (ind_t aqueous: m_data->range_aqueous()) { if (m_data->nu_aqueous(aqueous, component) != 0) { sum += m_data->nu_aqueous(aqueous, component)*m_second.secondary_concentration(node, aqueous); } } return sum - m_ideq.mobile_total_concentration(node, component, variable); } void DiffusionProgram::nodal_residual_mineral(ind_t node, const Variables& variable, Vector& residual) { for (ind_t mineral: m_data->range_mineral()) { ind_t id = m_ideq.id_equation_mineral(node, mineral); if (id != no_equation) { residual(id) = nodal_mineral_residual_mineral(node, mineral, variable); m_is_micp[id] = true; } } } double DiffusionProgram::nodal_mineral_residual_mineral(ind_t node, int mineral, const Variables& variable) { double res = m_data->logk_mineral(mineral); for (ind_t component: m_data->range_aqueous_component()) { res -= m_data->nu_mineral(mineral, component) * ( m_ideq.component_concentration(node, component, variable) + m_second.loggamma_component(node, component)); } return res; } // ================================== // // // // Jacobian // // // // ================================== // void DiffusionProgram::get_jacobian(Variables &variable, Vector &residual, SparseMatrix &jacobian, double alphadt) { list_triplet_t jacob; const ind_t ncomp = m_data->nb_component-1; const ind_t nmin = m_data->nb_mineral; - const ind_t estimation = m_mesh->nnodes()*(ncomp*(3+ncomp) + ncomp*(ncomp+nmin)+nmin+ncomp*2); + const ind_t estimation = m_mesh->nb_nodes()*(ncomp*(3+ncomp) + ncomp*(ncomp+nmin)+nmin+ncomp*2); jacob.reserve(estimation); jacobian_transport(variable, residual, jacob, alphadt); jacobian_speciation(variable, jacob, alphadt); jacobian = SparseMatrix(get_neq(), get_neq()); jacobian.setFromTriplets(jacob.begin(), jacob.end()); } // Transport // ========= void DiffusionProgram::jacobian_transport(Variables& variable, Vector& residual, list_triplet_t& jacobian, double alphadt) { for (ind_t element: m_mesh->range_elements()) { element_jacobian_transport(element, variable, residual, jacobian, alphadt); } } void DiffusionProgram::element_jacobian_transport(ind_t element, Variables& variable, Vector& residual, list_triplet_t& jacobian, double alphadt) { for (int component: m_data->range_aqueous_component()) { Eigen::VectorXd element_residual_orig(Eigen::VectorXd::Zero(2)); element_residual_transport_component(element, component, variable, element_residual_orig); for (int en: m_mesh->range_nen()) { Eigen::VectorXd element_residual(Eigen::VectorXd::Zero(2)); const ind_t node = m_mesh->get_node(element, en); const ind_t idc = m_ideq.id_equation_diffusion(node, component); const ind_t dof = m_ideq.get_dof_diffusion(node, component); if (idc == no_equation) continue; const double tmp_v = variable.velocity(dof); const double tmp_d = variable.displacement(dof); double h = EPS_J*std::abs(tmp_v); if (h == 0) h = EPS_J; variable.velocity(dof) = tmp_v + h; h = variable.velocity(dof) - tmp_v; variable.displacement(dof) = tmp_d + alphadt*h; element_residual_transport_component(element, component, variable, element_residual); variable.velocity(dof) = tmp_v; variable.displacement(dof) = tmp_d; for (int enr: m_mesh->range_nen()) { const ind_t noder = m_mesh->get_node(element, enr); const ind_t idr = m_ideq.id_equation_diffusion(noder, component); if (idr == no_equation) continue; jacobian.push_back(triplet_t(idr, idc, (element_residual(enr) - element_residual_orig(enr))/h)); } // mineral -> not using finite difference for (ind_t mineral: m_data->range_mineral()) { const ind_t idm = m_ideq.id_equation_mineral(node, mineral); if (idm != no_equation) jacobian.push_back(triplet_t(idc, idm, - -alphadt*m_mesh->crosssection()*m_mesh->get_dx(element)* + -alphadt*m_mesh->get_volume_element(element)* m_data->nu_mineral(mineral, component)/(2.0))); } } } } // Speciation // ========== void DiffusionProgram::jacobian_speciation(Variables& variable, list_triplet_t& jacobian, double alphadt) { for (ind_t node: m_mesh->range_nodes()) { nodal_jacobian_speciation_fd(node, variable, jacobian, alphadt); //nodal_jacobian_massbalance(node, variable, jacobian, alphadt); //nodal_jacobian_mineral(node, variable, jacobian, alphadt); } } void DiffusionProgram::nodal_jacobian_speciation_fd( ind_t node, Variables& variable, list_triplet_t& jacobian, double alphadt) { m_second.nodal_solve_secondary_variables(node, variable); // Residuals Eigen::VectorXd residual_massbalance_orig(m_data->nb_component); for (int component: m_data->range_aqueous_component()) { //std::cout << nodal_component_residual_massbalance(node, component, variable)<< std::endl; residual_massbalance_orig(component) = nodal_component_residual_massbalance(node, component, variable); } Eigen::VectorXd residual_mineral_orig(m_data->nb_mineral); for (int mineral: m_data->range_mineral()) { residual_mineral_orig(mineral) = nodal_mineral_residual_mineral(node, mineral, variable); } for (int component: m_data->range_aqueous_component()) { const ind_t id_c = m_ideq.id_equation_massbalance(node, component); if (id_c == no_equation) continue; const ind_t dof_c = m_ideq.get_dof_massbalance(node, component); double tmp_v = variable.velocity(dof_c); double tmp_d = variable.displacement(dof_c); double h = std::abs(tmp_v)*EPS_J; if (h < EPS_J*1e-2) h = EPS_J; variable.velocity(dof_c) += h; h = variable.velocity(dof_c) - tmp_v; variable.displacement(dof_c) = update_massbalance(variable.velocity(dof_c), tmp_d, alphadt); m_second.nodal_secondary_concentrations(node, variable); //m_second.nodal_solve_secondary_variables(node, variable); // Aqueous mass balance equations for (int k: m_data->range_aqueous_component()) { const ind_t id_k = m_ideq.id_equation_massbalance(node, k); if (id_k == no_equation) continue; const double residual = nodal_component_residual_massbalance(node, k, variable); jacobian.push_back(triplet_t(id_k, id_c, (residual-residual_massbalance_orig(k))/h)); // contribution of total mobile concentration if (k == component) { const ind_t id_tc = m_ideq.id_equation_diffusion(node, component); if (id_tc != no_equation) { jacobian.push_back(triplet_t(id_c, id_tc, -alphadt)); } } } // Mineral stability for (int mineral: m_data->range_mineral()) { const ind_t id_m = m_ideq.id_equation_mineral(node, mineral); if (id_m == no_equation) continue; double residual = nodal_mineral_residual_mineral(node, mineral, variable); jacobian.push_back(triplet_t(id_m, id_c, (residual-residual_mineral_orig(mineral))/h)); } variable.velocity(dof_c) = tmp_v; variable.displacement(dof_c) = tmp_d; m_second.nodal_secondary_concentrations(node, variable); //m_second.nodal_solve_secondary_variables(node, variable); } // add dummy diagonal element for mineral => needed for reformulation for (int mineral: m_data->range_mineral()) { const ind_t id_m = m_ideq.id_equation_mineral(node, mineral); if (id_m == no_equation) continue; jacobian.push_back(triplet_t(id_m, id_m, 0)); } } void DiffusionProgram::nodal_jacobian_massbalance(ind_t node, const Variables& variable, list_triplet_t& jacobian, double alphadt) { const double logten = std::log(10.0); for (ind_t component: m_data->range_aqueous_component()) { const ind_t idp = m_ideq.id_equation_massbalance(node, component); if (idp == no_equation) continue; // total concentration const ind_t idc = m_ideq.id_equation_diffusion(node, component); if (idc != no_equation) { jacobian.push_back(triplet_t(idp, idc, -alphadt)); } for (ind_t k: m_data->range_aqueous_component()) { // concentration const ind_t ids = m_ideq.id_equation_massbalance(node, k); if (ids == no_equation) continue; double tmp_iip = 0; if (k == component) tmp_iip += logten*pow10(m_ideq.component_concentration(node, component, variable)); for (ind_t aqueous: m_data->range_aqueous()) { if (m_data->nu_aqueous(aqueous, component) == 0.0) continue; tmp_iip += (logten*m_data->nu_aqueous(aqueous, component)* m_data->nu_aqueous(aqueous, k)*m_second.secondary_concentration(node, aqueous)); // } jacobian.push_back(triplet_t(idp, ids, alphadt*tmp_iip)); } } } void DiffusionProgram::nodal_jacobian_mineral(ind_t node, const Variables& variable, list_triplet_t& jacobian, double alphadt) { for (ind_t mineral: m_data->range_mineral()) { const ind_t idm = m_ideq.id_equation_mineral(node, mineral); if (idm == no_equation) continue; for (ind_t component: m_data->range_aqueous_component()) { const ind_t idc = m_ideq.id_equation_massbalance(node, component); if (idc == no_equation) continue; jacobian.push_back(triplet_t(idm, idc, -alphadt*m_data->nu_mineral(mineral, component))); } // make place for an element jacobian.push_back(triplet_t(idm, idm, 0)); } } // Variables // ========= void DiffusionProgram::update_variables(ParabolicVariables& x, const Eigen::VectorXd& predictor, const Eigen::VectorXd& update, double factor, double alpha_dt ) { // for (ind_t node: m_mesh->range_nodes()) // { // for (ind_t edof=0; edofrange_nodes()) { for (int component: m_data->range_aqueous_component()) { const ind_t id_eq_d = m_ideq.id_equation_diffusion(node, component); if (id_eq_d != no_equation) { const ind_t dof_d = m_ideq.get_dof_diffusion(node, component); x.velocity(dof_d) += factor*update(id_eq_d); x.displacement(dof_d) = update_diffusion(x.velocity(dof_d), predictor(dof_d), alpha_dt); } const ind_t id_eq_b = m_ideq.id_equation_massbalance(node, component); if (id_eq_b != no_equation) { const ind_t dof_b = m_ideq.get_dof_massbalance(node, component); x.velocity(dof_b) += factor*update(id_eq_b); x.displacement(dof_b) = update_massbalance(x.velocity(dof_b), predictor(dof_b), alpha_dt); } } for (int mineral: m_data->range_mineral()) { const ind_t id_eq_m = m_ideq.id_equation_mineral(node, mineral); if (id_eq_m != no_equation) { const ind_t dof_m = m_ideq.get_dof_mineral(node, mineral); x.velocity(dof_m) += factor*update(id_eq_m); x.displacement(dof_m) = update_mineral(x.velocity(dof_m), predictor(dof_m), alpha_dt); } } } m_second.solve_secondary_variables(x); } void DiffusionProgram::set_predictor(ParabolicVariables& x, Eigen::VectorXd& predictor, double alpha, double dt) { - predictor.resize(m_mesh->nnodes()*get_ndf()); + predictor.resize(m_mesh->nb_nodes()*get_ndf()); predictor = x.displacement + (1-alpha)*dt*x.velocity; x.velocity.setZero(); WARNING << "Predictor size : " << predictor.rows(); } // compute secondary conc bool DiffusionProgram::hook_start_iteration(const Variables& x, double norm_residual) { return m_second.solve_secondary_variables(x); } double DiffusionProgram::maximum_lambda(const ParabolicVariables& x, const Eigen::VectorXd& update, double alphadt) { double inv_maximum_lambda = 1.0; // for (ind_t node: m_mesh->range_nodes()) // { // for (int mineral: m_data->range_mineral()) // { // if (m_ideq.id_equation_mineral(node, mineral) == no_equation // or x.displacement(m_ideq.get_dof_mineral(node, mineral)) == 0) continue; // inv_maximum_lambda = std::max(inv_maximum_lambda, // -alphadt*update(m_ideq.id_equation_mineral(node, mineral))/ // (0.9*x.displacement(m_ideq.get_dof_mineral(node, mineral)))); // } // } return 1.0/inv_maximum_lambda; } void DiffusionProgram::initialize_no_bc_with_init(const EquilibriumState& initialstate, ParabolicVariables& var) { for (ind_t node: m_mesh->range_nodes()) { if (not m_ideq.node_has_bc(node)) { for (ind_t component: m_data->range_aqueous_component()) { m_ideq.component_concentration(node, component, var) = std::log10(initialstate.molality_component(component)); m_ideq.mobile_total_concentration(node, component, var) = initialstate.total_aqueous_concentration_component(component); m_second.loggamma_component(node, component) = initialstate.loggamma_component(component); } for (ind_t aqueous: m_data->range_aqueous()) { m_second.secondary_concentration(node, aqueous) = initialstate.molality_secondary()[aqueous]; m_second.loggamma_secondary(node, aqueous) = initialstate.loggamma_secondary(aqueous); } for (ind_t mineral: m_data->range_mineral()) { - m_ideq.mole_mineral(node, mineral, var) = initialstate.moles_mineral(mineral)/m_mesh->cell_volume(1); + m_ideq.mole_mineral(node, mineral, var) = initialstate.moles_mineral(mineral)/m_mesh->get_volume_cell(1); } } } } void DiffusionProgram::copy_cp_variables(const ParabolicVariables& x, Eigen::VectorXd &moles_mineral) { moles_mineral.resize(get_neq()); moles_mineral.setZero(); for (ind_t node: m_mesh->range_nodes()) { for (int mineral: m_data->range_mineral()) { const int id_eq = m_ideq.id_equation_mineral(node, mineral); if (id_eq == no_equation) continue; moles_mineral(id_eq) = m_ideq.mole_mineral(node, mineral, x); } } } void DiffusionProgram::projection(ParabolicVariables& x) { for (int node: m_mesh->range_nodes()) { for (int mineral: m_data->range_mineral()) { if (m_ideq.mole_mineral(node, mineral, x) < 1e-10) { m_ideq.mole_mineral(node, mineral, x) = 0.0; } } } } } // end namespace systems } // end namespace reactmicp } // end namespace specmicp diff --git a/src/reactmicp/systems/diffusion/diffusion.hpp b/src/reactmicp/systems/diffusion/diffusion.hpp index 0c07edf..5a610d9 100644 --- a/src/reactmicp/systems/diffusion/diffusion.hpp +++ b/src/reactmicp/systems/diffusion/diffusion.hpp @@ -1,306 +1,306 @@ /*------------------------------------------------------- - Module : reactmicp/systems/diffusion - File : diffusion.hpp - Author : Fabien Georget Copyright (c) 2014, Fabien 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: * Redistributions of source code must retain the above copyright notice, this list of conditions and the following disclaimer. * 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. * Neither the name of the Princeton University 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 OWNER 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_DIFFUSION_HPP #define SPECMICP_REACTMICP_SYSTEMS_DIFFUSION_HPP #include #include #include "database/data_container.hpp" #include "reactmicp/common.hpp" #include "reactmicp/micpsolvers/parabolicprogram.hpp" #include "reactmicp/meshes/mesh1d.hpp" #include "diffusion_numbering.hpp" #include "diffusion_secondary.hpp" namespace specmicp { namespace reactmicp { namespace systems { //! \brief Parameters for the diffusion system class DiffusionParameter { public: DiffusionParameter(double diffusion_coefficient, double porosity): m_diff(diffusion_coefficient), m_poro(porosity) {} //! \brief Density of water (kg/m^3) double density_water() {return 1e3;} //! \brief Diffusion coefficient (element by element) (m^2/s) double diffusion_coefficient(ind_t element) {return m_diff;} //! \brief Porosity (at a node (function of composition)) double porosity(ind_t node) {return m_poro;} private: double m_diff; double m_poro; }; //! \brief Saturated diffusion in reactive porous media class DiffusionProgram: public solvers::ParabolicProgram { public: using Variables = ParabolicVariables; const int no_equation = -1; DiffusionProgram(std::shared_ptr themesh, std::shared_ptr thedatabase, std::shared_ptr parameters, diffusion::ListBoundaryCondition& list_bcs, Variables& var): m_mesh(themesh), m_data(thedatabase), m_param(parameters), - m_ideq(themesh->nnodes(), thedatabase, list_bcs, var), - m_second(themesh->nnodes(), thedatabase, m_ideq) + m_ideq(themesh->nb_nodes(), thedatabase, list_bcs, var), + m_second(themesh->nb_nodes(), thedatabase, m_ideq) { } DiffusionProgram(std::shared_ptr themesh, std::shared_ptr thedatabase, std::shared_ptr parameters, diffusion::ListBoundaryCondition& list_bcs, const EquilibriumState &initialstate, Variables& var); //! \brief Return the number of degree of freedoms ind_t get_neq() const {return m_ideq.get_neq();} ind_t get_ndf() const {return m_ideq.get_ndf();} //! \brief Return the residuals void get_residuals(const ParabolicVariables& variable, Vector& residual); //! brief Return the jacobian void get_jacobian(ParabolicVariables& variable, Vector& residual, SparseMatrix& jacobian, double alphadt); //! \brief Return the index of the MiCP equations const std::vector& micp_index() const { return m_is_micp;} //! \brief Update the variables void update_variables(ParabolicVariables& x, const Eigen::VectorXd& predictor, const Eigen::VectorXd& update, double factor, double alpha_dt ); //! \brief set the predictor void set_predictor(ParabolicVariables& x, Eigen::VectorXd& predictor, double alpha, double dt); //! \brief compute secondary concentrations bool hook_start_iteration(const ParabolicVariables& x, double norm_residual); //! \brief Return the maximum lambda allowed for the linesearch double maximum_lambda(const ParabolicVariables& x, const Eigen::VectorXd& update, double alphadt); //! \brief copy mineral amount from variables void copy_cp_variables(const ParabolicVariables& x, Eigen::VectorXd& moles_mineral); //! \brief projection step void projection(ParabolicVariables& x); // todo or not........ // double update_equation(int id_eq, // double update, // double predictor, // double alphadt // ) // { // // find _id eq .... // } private: // Residuals // ========= // transport // --------- //! \brief Compute the residuals for the diffusion equations void residual_transport(const Variables& variable, Vector& residual); //! \brief Contribution of the diffusion operator to the residuals void element_residual_transport(ind_t element, const Variables& variable, Vector& residual); //! \brief Contribution of the minerals to the residuals of the diffusion equations double nodal_mineral_transient_term_transport(ind_t node, ind_t component, const Variables& variable); void element_residual_transport_component( ind_t element, int component, const Variables& variable, Vector& element_residual); // speciation // ---------- //! \brief Compute the residual for the speciation problem void residual_speciation(const Variables& variable, Vector& residual); //! \brief Compute the residuals of the aqueous mass balances void nodal_residual_massbalance(ind_t node, const Variables& variable, Vector& residual); //! \brief compute the residual for the aqueous mass balance of 'component' at 'node' double nodal_component_residual_massbalance(ind_t node, int component, const Variables& variable); //! \brief Compute the residuals of the mineral speciation void nodal_residual_mineral(ind_t node, const Variables& variable, Vector& residual); //! \brief compute the residual for the stability equation of 'mineral' at 'node' double nodal_mineral_residual_mineral(ind_t node, int mineral, const Variables& variable); // Jacobian // ======== // transport // --------- //! \brief Contribution of the diffusion equation to the jacobian void jacobian_transport(Variables& variable, Vector& residual, list_triplet_t& jacobian, double alphadt); //! \brief Element by element assembly procedure void element_jacobian_transport(ind_t element, Variables& variable, Vector& residual, list_triplet_t& jacobian, double alphadt); // speciation // ---------- //! \brief Contribution of the speciation problem to the mass balance void jacobian_speciation(Variables &variable, list_triplet_t& jacobian, double alphadt); //! \brief Contribution of the mass balance at 'node' to the jacobian void nodal_jacobian_massbalance(ind_t node, const Variables& variable, list_triplet_t& jacobian, double alphadt); //! \brief Contribution of the mineral at equilibrium problem to the jacobian void nodal_jacobian_mineral(ind_t node, const Variables& variable, list_triplet_t& jacobian, double alphadt); //! \brief Finite difference jacobian for speciation void nodal_jacobian_speciation_fd(ind_t node, Variables &variable, list_triplet_t& jacobian, double alphadt); // Initialization // ============== void initialize_no_bc_with_init(const EquilibriumState& initialstate, ParabolicVariables& var); // Update // ====== double update_diffusion(double update, double predictor, double alphadt) { return predictor + alphadt*update; } double update_massbalance(double update, double predictor, double alphadt) { return predictor + alphadt*update; } double update_mineral(double update, double predictor, double alphadt) { return predictor + alphadt*update; } // Member // ======= std::shared_ptr m_mesh; std::shared_ptr m_data; std::shared_ptr m_param; diffusion::DiffusionNumbering m_ideq; diffusion::DiffusionSecondaryVariables m_second; std::vector m_is_micp; //SparseVector m_boundary_condition; }; } // end namespace systems } // end namespace reactmicp } // end namespace specmicp #endif // SPECMICP_REACTMICP_SYSTEMS_DIFFUSION_HPP diff --git a/src/reactmicp/systems/saturated_diffusion/reactive_transport_solver.cpp b/src/reactmicp/systems/saturated_diffusion/reactive_transport_solver.cpp index e5825d0..5828d4d 100644 --- a/src/reactmicp/systems/saturated_diffusion/reactive_transport_solver.cpp +++ b/src/reactmicp/systems/saturated_diffusion/reactive_transport_solver.cpp @@ -1,277 +1,277 @@ //#define EIGEN_DEFAULT_IO_FORMAT IOFormat(8) #include "reactive_transport_solver.hpp" #include "specmicp/reduced_system_solver.hpp" #include "boundary_conditions.hpp" #include "reactmicp/meshes/mesh1d.hpp" #ifdef SPECMICP_USE_OPENMP #include #endif // SPECMICP_USE_OPENMP #include namespace specmicp { namespace reactmicp { namespace internals { //! \brief Store residual information struct SIAResiduals { index_t nb_iterations; scalar_t transport_residual_0; scalar_t transport_residual; scalar_t transport_update; Vector speciation_residual; Vector speciation_update; SIAResiduals(index_t nb_node): nb_iterations(0), speciation_residual(Vector::Zero(nb_node)), speciation_update(Vector::Zero(nb_node)) {} }; } // end namespace internals SIASaturatedReactiveTransportSolver::SIASaturatedReactiveTransportSolver( std::shared_ptr the_mesh, RawDatabasePtr the_database, std::shared_ptr transport_parameters ): m_mesh(the_mesh), m_database(the_database), m_variables(the_database, the_mesh, transport_parameters), m_transport_program(the_mesh, the_database, transport_parameters), m_transport_solver(m_transport_program), m_bcs(the_mesh->nb_nodes()) {} void SIASaturatedReactiveTransportSolver::apply_boundary_conditions(SIABoundaryConditions bcs) { - specmicp_assert(bcs.bs_types.size() == (unsigned int) m_mesh->nnodes()); - specmicp_assert(bcs.initial_states.size() == (unsigned int) m_mesh->nnodes()); + specmicp_assert(bcs.bs_types.size() == (unsigned int) m_mesh->nb_nodes()); + specmicp_assert(bcs.initial_states.size() == (unsigned int) m_mesh->nb_nodes()); std::vector list_fixed_nodes; for (index_t node: m_mesh->range_nodes()) { switch (bcs.bs_types[node]) { case BoundaryConditionType::FixedComposition: m_bcs(node) = 1; list_fixed_nodes.push_back(node); break; default: m_bcs(node) = 0; break; } m_variables.update_composition(node, bcs.list_initial_states[bcs.initial_states[node]]); m_transport_program.number_equations(list_fixed_nodes); } } scalar_t SIASaturatedReactiveTransportSolver::residuals_transport() { Eigen::VectorXd residuals; m_transport_solver.compute_residuals(m_variables.total_mobile_concentrations(), residuals); return residuals.norm(); } SIASaturatedReactiveTransportSolverReturnCode convert_dfpm_retcode(dfpmsolver::ParabolicDriverReturnCode retcode) { switch (retcode) { case dfpmsolver::ParabolicDriverReturnCode::MaxIterations: return SIASaturatedReactiveTransportSolverReturnCode::MaxIterationsInTransport; case dfpmsolver::ParabolicDriverReturnCode::StationaryPoint: return SIASaturatedReactiveTransportSolverReturnCode::StationaryPointsInTransport; default: return SIASaturatedReactiveTransportSolverReturnCode::ErrorInTransport; } } SIASaturatedReactiveTransportSolverReturnCode SIASaturatedReactiveTransportSolver::solve_timestep(scalar_t dt) { retcode_t return_code = retcode_t::NotConvergedYet; internals::SIAResiduals residuals_storage(m_mesh->nb_nodes()); m_flag_compute_speciation = Eigen::VectorXi::Ones(m_mesh->nb_nodes()); // initialization transport program and transport solver Eigen::VectorXd transport_variable(m_variables.total_mobile_concentrations()); m_transport_solver.get_options() = get_options().transport_options; m_transport_solver.get_velocity() = Eigen::VectorXd::Zero(transport_variable.rows()); m_transport_program.external_flux().setZero(); // residuals_storage.transport_residual_0 = residuals_transport(); m_minerals_start_iteration = m_variables.speciation_variables().block( m_database->nb_component, 0, m_database->nb_mineral, m_mesh->nb_nodes()); // solve transport // --------------- dfpmsolver::ParabolicDriverReturnCode retcode = m_transport_solver.solve_timestep(dt, transport_variable); if (retcode != dfpmsolver::ParabolicDriverReturnCode::ResidualMinimized) { ERROR << "Failed to solve 1st transport problem; return code " << (int) retcode; return convert_dfpm_retcode(retcode); } residuals_storage.transport_residual = m_transport_solver.get_perfs().current_residual; residuals_storage.transport_update = m_transport_solver.get_perfs().current_update; // update m_variables.total_mobile_concentrations() = transport_variable; // Speciation // ---------- restart_speciation_iterations(dt, residuals_storage); // Convergence // ----------- return_code = check_convergence(residuals_storage); residuals_storage.nb_iterations +=1; // If SNIA, we stop if (get_options().max_iterations == 1) return retcode_t::Success; // SIA algorithm // ------------- while (return_code == retcode_t::NotConvergedYet) { //std::cout << m_flag_compute_speciation.sum() << std::endl; transport_variable = m_variables.total_mobile_concentrations(); retcode = m_transport_solver.restart_timestep(transport_variable); if (retcode != dfpmsolver::ParabolicDriverReturnCode::ResidualMinimized) { if (retcode == dfpmsolver::ParabolicDriverReturnCode::MaxIterations and m_transport_solver.get_perfs().current_residual/residuals_storage.transport_residual_0 < get_options().good_enough_transport_residual_tolerance) { restart_speciation_iterations(dt, residuals_storage); return_code = retcode_t::Success; break; } ERROR << "Failed to solve transport problem; return code " << (int) retcode; return convert_dfpm_retcode(retcode); } residuals_storage.transport_residual = m_transport_solver.get_perfs().current_residual; residuals_storage.transport_update = m_transport_solver.get_perfs().current_update; m_variables.total_mobile_concentrations() = transport_variable; residuals_storage.nb_iterations +=1; restart_speciation_iterations(dt, residuals_storage); return_code = check_convergence(residuals_storage); } std::cout << "End fixed point : return code " << (int) return_code << " - nb iter : " <nnodes(); ++node) #else for (index_t node: m_mesh->range_nodes()) #endif { if (m_bcs(node) == 0 and m_flag_compute_speciation(node)) { solve_speciation_node(dt, node, residuals_storage); } } #ifdef _OPENMP } #endif } void SIASaturatedReactiveTransportSolver::solve_speciation_node( scalar_t dt, index_t node, internals::SIAResiduals& residuals_storage ) { Eigen::VectorXd totconc; m_variables.nodal_update_total_amount(node, totconc); ReducedSystemSolver solver(m_database, totconc, get_options().speciation_options); //solver.get_options().solver_options.max_iter = 50; //solver.get_options().solver_options.factor_gradient_search_direction = 400; //solver.get_options().solver_options.use_scaling = false; Eigen::VectorXd variables = m_variables.speciation_variables(node); micpsolver::MiCPPerformance perf = solver.solve(variables); if (perf.return_code != micpsolver::MiCPSolverReturnCode::ResidualMinimized) { ERROR << "Failed to solve speciation solver for node " << node << "; return code = " << (int) perf.return_code; } residuals_storage.speciation_residual(node) = perf.current_residual; residuals_storage.speciation_update(node) = perf.current_update; if (perf.current_update < get_options().threshold_update_disable_speciation) m_flag_compute_speciation(node) = 0; EquilibriumState solution = solver.get_solution(variables); m_variables.update_composition(node, solution); // mineral const scalar_t denom = 1e-3*m_mesh->get_volume_cell(node)*dt; for (index_t component: m_database->range_aqueous_component()) { scalar_t sum = 0.0; for (index_t mineral: m_database->range_mineral()) { if (m_database->nu_mineral(mineral, component) == 0.0) continue; sum -= m_database->nu_mineral(mineral, component)*( m_variables.mineral_amount(node, mineral) -m_minerals_start_iteration(mineral, node) ); } m_transport_program.external_flux(node, component) = sum/denom; } } SIASaturatedReactiveTransportSolverReturnCode SIASaturatedReactiveTransportSolver::check_convergence( internals::SIAResiduals& residuals_storage ) { // compute residuals of transport const scalar_t norm_transport = residuals_transport(); retcode_t return_code = retcode_t::NotConvergedYet; if (norm_transport/residuals_storage.transport_residual_0 < get_options().residual_tolerance) { return_code = retcode_t::Success; } else if (residuals_storage.transport_update < get_options().threshold_update_transport) { return_code = retcode_t::Success; } else if (residuals_storage.nb_iterations >= get_options().max_iterations) { WARNING << "Maximum number of iterations (" << get_options().max_iterations << ") reached during fixed-point iterations"; return_code = retcode_t::MaxIterations; } return return_code; } } // end namespace reactmicp } // end namespace specmicp diff --git a/src/reactmicp/systems/saturated_diffusion/transport_program.cpp b/src/reactmicp/systems/saturated_diffusion/transport_program.cpp index f02212a..34b95d9 100644 --- a/src/reactmicp/systems/saturated_diffusion/transport_program.cpp +++ b/src/reactmicp/systems/saturated_diffusion/transport_program.cpp @@ -1,218 +1,220 @@ #include "transport_program.hpp" #include "transport_parameters.hpp" #include "reactmicp/meshes/mesh1d.hpp" #include #define EPS_J 1e-8 namespace specmicp { namespace reactmicp { namespace systems { namespace siasaturated { SaturatedDiffusionProgram::SaturatedDiffusionProgram( std::shared_ptr the_mesh, std::shared_ptr the_data, std::shared_ptr the_param ): m_mesh(the_mesh), m_data(the_data), m_param(the_param), m_external_flux(Vector::Zero(the_mesh->nb_nodes()*(the_data->nb_component-1))) {} void SaturatedDiffusionProgram::number_equations(std::vector list_fixed_nodes ) { m_ideq = Eigen::VectorXd::Zero(get_ndf()*m_mesh->nb_nodes()); for (auto it=list_fixed_nodes.begin(); itrange_aqueous_component()) { m_ideq(get_dof_component(*it, component)) = no_equation; } } index_t neq = 0; for (index_t dof=0; dofrange_elements()) { element_residual_transport(element, displacement, velocity, residual); } } -void SaturatedDiffusionProgram::element_residual_transport(index_t element, - const Eigen::VectorXd& displacement, - const Eigen::VectorXd& velocity, - Eigen::VectorXd& residual) +void SaturatedDiffusionProgram::element_residual_transport( + index_t element, + const Vector& displacement, + const Vector& velocity, + Vector& residual) { Eigen::VectorXd element_residual(2); for (index_t component: m_data->range_aqueous_component()) { element_residual_transport_component(element, component, displacement, velocity, element_residual); for (index_t en: m_mesh->range_nen()) { const index_t node = m_mesh->get_node(element, en); const index_t id = m_ideq(get_dof_component(node, component)); if (id != no_equation) {residual(id) += element_residual(en);} } } } -void SaturatedDiffusionProgram::element_residual_transport_component(index_t element, +void SaturatedDiffusionProgram::element_residual_transport_component( + index_t element, index_t component, const Vector &displacement, const Vector &velocity, Vector& element_residual) { Eigen::Matrix jacob; Eigen::Matrix evelocity, edisplacement; scalar_t mass_coeff = -(m_param->density_water() * m_mesh->get_volume_element(element)/2); const scalar_t porosity = (m_param->porosity(m_mesh->get_node(element, 0)) + m_param->porosity(m_mesh->get_node(element, 1)) )/2; const scalar_t inv_diff_coeff = (0.5/m_param->diffusion_coefficient(m_mesh->get_node(element, 0)) + 0.5/m_param->diffusion_coefficient(m_mesh->get_node(element, 1))); scalar_t flux_coeff = -( m_mesh->get_face_area(element) / m_mesh->get_dx(element) * porosity * m_param->density_water() * 1.0/inv_diff_coeff ); jacob << 1, -1, -1, 1; jacob *= flux_coeff; evelocity << velocity(get_dof_component(m_mesh->get_node(element, 0), component))* mass_coeff*m_param->porosity(m_mesh->get_node(element, 0)), velocity(get_dof_component(m_mesh->get_node(element, 1), component))* mass_coeff*m_param->porosity(m_mesh->get_node(element, 0)); edisplacement << displacement(get_dof_component(m_mesh->get_node(element, 0), component)), displacement(get_dof_component(m_mesh->get_node(element, 1), component)); element_residual = evelocity+jacob*edisplacement; for (index_t en: m_mesh->range_nen()) { element_residual(en) += (m_mesh->get_volume_element(element)/2 *external_flux(m_mesh->get_node(element, en), component)); } } void SaturatedDiffusionProgram::compute_jacobian( Vector& displacement, Vector& velocity, Eigen::SparseMatrix& jacobian, scalar_t alphadt ) { list_triplet_t jacob; const index_t ncomp = m_data->nb_component-1; const index_t estimation = m_mesh->nb_nodes()*(ncomp*m_mesh->nen); jacob.reserve(estimation); for (index_t element: m_mesh->range_elements()) { element_jacobian_transport(element, displacement, velocity, jacob, alphadt); } jacobian = Eigen::SparseMatrix(get_neq(), get_neq()); jacobian.setFromTriplets(jacob.begin(), jacob.end()); } // Transport // ========= void SaturatedDiffusionProgram::element_jacobian_transport( index_t element, Vector& displacement, Vector& velocity, list_triplet_t& jacobian, scalar_t alphadt) { for (index_t component: m_data->range_aqueous_component()) { Eigen::VectorXd element_residual_orig(Eigen::VectorXd::Zero(2)); element_residual_transport_component(element, component, displacement, velocity, element_residual_orig); for (index_t en: m_mesh->range_nen()) { Eigen::VectorXd element_residual(Eigen::VectorXd::Zero(2)); const index_t node = m_mesh->get_node(element, en); const index_t dof = get_dof_component(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_J*std::abs(tmp_v); if (h == 0) h = EPS_J; velocity(dof) = tmp_v + h; h = velocity(dof) - tmp_v; displacement(dof) = tmp_d + alphadt*h; element_residual_transport_component(element, component, displacement, velocity, element_residual); velocity(dof) = tmp_v; displacement(dof) = tmp_d; for (index_t enr: m_mesh->range_nen()) { const index_t noder = m_mesh->get_node(element, enr); const index_t idr = m_ideq(get_dof_component(noder, component)); if (idr == no_equation) continue; jacobian.push_back(triplet_t(idr, idc, (element_residual(enr) - element_residual_orig(enr))/h)); } } } } void SaturatedDiffusionProgram::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: m_data->range_aqueous_component()) { const index_t dof = get_dof_component(node, component); const index_t id = m_ideq(dof); if (id == no_equation) continue; velocity(dof) += lambda*update(id); } } displacement = predictor + alpha_dt*velocity; } } // end namespace siasaturated } // end namespace systems } // end namespace reactmicp } // end namespace specmicp diff --git a/src/reactmicp/systems/saturated_diffusion/transport_program.hpp b/src/reactmicp/systems/saturated_diffusion/transport_program.hpp index bcdc008..bddd32f 100644 --- a/src/reactmicp/systems/saturated_diffusion/transport_program.hpp +++ b/src/reactmicp/systems/saturated_diffusion/transport_program.hpp @@ -1,133 +1,133 @@ #ifndef SPECMICP_REACTMICP_SATURATED_DIFFUSION_TRANSPORTPROGRAM_HPP #define SPECMICP_REACTMICP_SATURATED_DIFFUSION_TRANSPORTPROGRAM_HPP #include #include "database/data_container.hpp" #include "dfpmsolver/parabolic_program.hpp" #include "reactmicp/meshes/mesh1dfwd.hpp" namespace specmicp { namespace reactmicp { namespace systems { namespace siasaturated { class SaturatedDiffusionTransportParameters; //! \brief Transport program for saturated diffusion class SaturatedDiffusionProgram: dfpmsolver::ParabolicProgram { public: using triplet_t = Eigen::Triplet; //!< Triplet type for building the transport matrix using list_triplet_t = std::vector; //!< List of triplet for building the sparse matrix SaturatedDiffusionProgram( std::shared_ptr the_mesh, std::shared_ptr the_data, std::shared_ptr the_param ); //! \brief Return the number of equations index_t get_neq() const {return m_neq;} //! \brief Return the number of degrees of freedom per node index_t get_ndf() const {return m_data->nb_component-1;} // number the equations void number_equations(std::vector list_fixed_nodes); //! \brief Compute the residuals void compute_residuals(const Vector &displacement, const Vector &velocity, Vector &residual ); //! \brief Compute the jacobian void compute_jacobian(Vector& displacement, Vector& velocity, Eigen::SparseMatrix &jacobian, scalar_t alphadt); //! \brief Update the solutions void update_solution(const Vector& update, scalar_t lambda, scalar_t alpha_dt, Vector& predictor, Vector& displacement, Vector& velocity); //! \brief Apply boundary conditions to the velocity vector //! //! by default do nothing. void apply_bc(scalar_t dt, const Vector& displacement, Vector velocity) {} //! Return the degree of freedom number of 'component' at 'node' index_t get_dof_component(index_t node, index_t component) const { specmicp_assert(component >=1 and component < m_data->nb_component); - specmicp_assert(node >= 0 and node < m_mesh->nnodes()); + //specmicp_assert(node >= 0 and node < m_mesh->nnodes()); return (component - 1)+get_ndf()*node; } //! Return the degree of freedom number of 'component' at 'node' index_t get_dof_mineral(index_t node, index_t mineral) const { specmicp_assert(mineral >= 0 and mineral < m_data->nb_mineral); - specmicp_assert(node >= 0 and node < m_mesh->nnodes()); + //specmicp_assert(node >= 0 and node < m_mesh->nnodes()); return mineral+m_data->nb_mineral*node; } //! Return the total mobile concentration scalar_t& get_total_mobile_concentration(index_t node, index_t component, Vector& concentrations) const { return concentrations.coeffRef(get_dof_component(node, component)); } //! Return the total mobile concentration scalar_t get_total_mobile_concentration(index_t node, index_t component, const Vector& concentrations) const { return concentrations.coeff(get_dof_component(node, component)); } //! Return the external flux Vector& external_flux() {return m_external_flux;} //! Return the external flux for 'component' at 'node' scalar_t& external_flux(index_t node, index_t component) {return m_external_flux(get_dof_component(node, component));} private: void element_residual_transport( index_t element, - const Eigen::VectorXd& displacement, - const Eigen::VectorXd& velocity, - Eigen::VectorXd& residual); + const Vector& displacement, + const Vector& velocity, + Vector& residual); void element_residual_transport_component( index_t element, index_t component, const Vector& displacement, const Vector& velocity, Vector& element_residual); void element_jacobian_transport( index_t element, Vector& displacement, Vector& velocity, list_triplet_t& jacobian, scalar_t alphadt); index_t m_neq; Vector m_ideq; mesh::Mesh1DPtr m_mesh; database::RawDatabasePtr m_data; std::shared_ptr m_param; Vector m_external_flux; }; } // end namespace siasaturated } // end namespace systems } // end namespace reactmicp } // end namespace specmicp #endif // SPECMICP_REACTMICP_SATURATED_DIFFUSION_TRANSPORTPROGRAM_HPP diff --git a/src/reactmicp/systems/secondary_variables/secondary.inl b/src/reactmicp/systems/secondary_variables/secondary.inl index 81bdb31..73ae5d5 100644 --- a/src/reactmicp/systems/secondary_variables/secondary.inl +++ b/src/reactmicp/systems/secondary_variables/secondary.inl @@ -1,226 +1,226 @@ #include "secondary.hpp" #include "physics/laws.hpp" namespace specmicp { namespace reactmicp { // ================================== // // // // Secondary variables wibbly-timey // // // // ================================== // //! \brief Set secondary concentrations template void SecondaryVariables::compute_secondary_concentrations() { for (index_t node=0; node void SecondaryVariables::nodal_secondary_concentrations(index_t node) { for (int species: m_data->range_aqueous()) { scalar_t sum = -m_data->logk_aqueous(species) - log_gamma_secondary(node, species); for (int component: m_data->range_aqueous_component()) { sum += m_data->nu_aqueous(species, component)*( log_gamma_component(node, component) + log_component_concentration(node, component) ); } secondary_concentration(node, species) = pow10(sum); } } template double SecondaryVariables::norm_secondary_variables() { return m_secondary_variables.block( 1, 0, m_data->nb_aqueous+(m_data->nb_component-1), m_secondary_variables.cols() ).norm(); } template double SecondaryVariables::nodal_norm_secondary_variables(index_t node) { return m_secondary_variables.block( 1, node, m_data->nb_aqueous+(m_data->nb_component-1), 1 ).norm(); } //! \brief Solve for secondary variables template int SecondaryVariables::solve_secondary_variables() { bool may_have_converged = true; for (index_t node=0; node int SecondaryVariables::nodal_solve_secondary_variables(index_t node) { bool may_have_converged = false; scalar_t previous_norm = nodal_norm_secondary_variables(node); - for (int i=0; i<6; ++i) + for (int i=0; i<10; ++i) { Vector copy = m_secondary_variables.col(node); // to compute update nodal_secondary_concentrations(node); nodal_secondary_variables(node); double new_norm = nodal_norm_secondary_variables(node); - if ( std::abs(previous_norm - new_norm)/previous_norm < 1e-6 - or (m_secondary_variables.col(node) - copy).norm() < 1e-6 ) + if ( std::abs(previous_norm - new_norm)/previous_norm < 1e-10 + or (m_secondary_variables.col(node) - copy).norm() < 1e-10 ) { may_have_converged=true; break; } previous_norm = new_norm; } //std::cout << "node : " << node << " result : " << may_have_converged << std::endl; if (may_have_converged) m_secondary_is_ok[node] = true; return may_have_converged; } //! \brief Compute secondary variables template void SecondaryVariables::compute_secondary_variables() { for (index_t node=0; node void SecondaryVariables::nodal_secondary_variables(index_t node) { nodal_ionic_strength(node); nodal_log_gamma(node); } template void SecondaryVariables::compute_ionic_strength() { for (index_t node=0; node void SecondaryVariables::nodal_ionic_strength(index_t node) { scalar_t ionics = 0; for (index_t component: m_data->range_aqueous_component()) { if (m_data->param_aq(component, 0) == 0) continue; ionics += std::pow(m_data->param_aq(component, 0), 2)*component_concentration(node, component); } for (index_t aqueous: m_data->range_aqueous()) { index_t idaq = m_data->nb_component + aqueous; if (m_data->param_aq(idaq, 0) == 0) continue; ionics += std::pow(m_data->param_aq(idaq, 0), 2)*secondary_concentration(node, aqueous); } ionic_strength(node) = ionics/2; } template void SecondaryVariables::compute_log_gamma() { for (index_t node=0; node void SecondaryVariables::nodal_log_gamma(index_t node) { nodal_log_gamma_component(node); nodal_log_gamma_aqueous(node); } //! \brief compute the logarithm of the activity coefficients for node 'node' template void SecondaryVariables::nodal_log_gamma_component(index_t node) { const scalar_t sqrti = std::sqrt(ionic_strength(node)); for (index_t component: m_data->range_aqueous_component()) { log_gamma_component(node, component) = laws::extended_debye_huckel(ionic_strength(node), sqrti, m_data->param_aq(component, 0), m_data->param_aq(component, 1), m_data->param_aq(component, 2)); } } //! \brief compute the logarithm of the activity coefficients for node 'node' template void SecondaryVariables::nodal_log_gamma_aqueous(index_t node) { const scalar_t sqrti = std::sqrt(ionic_strength(node)); for (int aqueous: m_data->range_aqueous()) { index_t idaq = m_data->nb_component+aqueous; log_gamma_secondary(node, aqueous) = laws::extended_debye_huckel(ionic_strength(node), sqrti, m_data->param_aq(idaq, 0), m_data->param_aq(idaq, 1), m_data->param_aq(idaq, 2)); } } template scalar_t SecondaryVariables::charge_balance(index_t node) { specmicp_assert(m_secondary_is_ok[node]); scalar_t balance = 0.0; for (index_t component: m_data->range_aqueous_component()) { if (m_data->charge_component(component) == 0.0) continue; balance += m_data->charge_component(component)*component_concentration(node, component); } for (index_t aqueous: m_data->range_aqueous()) { if (m_data->charge_aqueous(aqueous) == 0.0) continue; balance += m_data->charge_aqueous(aqueous)*secondary_concentration(node, aqueous); } return balance; } template scalar_t SecondaryVariables::total_mobile_concentration(index_t node, index_t component) const { specmicp_assert(m_secondary_is_ok[node]); scalar_t totconc = component_concentration(node, component); for (index_t aqueous: m_data->range_aqueous()) { if (m_data->nu_aqueous(aqueous, component) == 0.0) continue; totconc += m_data->nu_aqueous(aqueous, component)*secondary_concentration(node, aqueous); } return totconc; } } // end namespace reactmicp } // end namespace specmicp diff --git a/src/specmicp/equilibrium_data.cpp b/src/specmicp/equilibrium_data.cpp index 3280381..bbd7c73 100644 --- a/src/specmicp/equilibrium_data.cpp +++ b/src/specmicp/equilibrium_data.cpp @@ -1,157 +1,157 @@ /*------------------------------------------------------- - Module : specmicp - File : equilibrium_data.cpp - Author : Fabien Georget Copyright (c) 2014, Fabien 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: * Redistributions of source code must retain the above copyright notice, this list of conditions and the following disclaimer. * 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. * Neither the name of the Princeton University 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 OWNER 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 "equilibrium_data.hpp" #include "database/module.hpp" #include "physics/constants.hpp" #define POW10(m) std::pow(10.0, m) namespace specmicp { scalar_t EquilibriumState::logIAP(int m) const { scalar_t logiap = 0.0; for (int i=1; inb_component; ++i) logiap += m_main_variables(i)*m_data->nu_mineral(m, i); return logiap; } scalar_t EquilibriumState::logIAP(const std::string& label_mineral) const { database::DatabaseModule namer(m_data); int idm = namer.safe_label_to_id(label_mineral, m_data->labels_minerals); return logIAP(idm); } scalar_t EquilibriumState::saturation_index(const std::string& label_mineral) const { database::DatabaseModule namer(m_data); int idm = namer.safe_label_to_id(label_mineral, m_data->labels_minerals); return logIAP(idm) - m_data->logk_mineral(idm); } scalar_t EquilibriumState::logIAP_kinetic(int m) const { scalar_t logiap = 0.0; for (int i=1; inb_component; ++i) logiap += m_main_variables(i)*m_data->nu_mineral_kinetic(m, i); return logiap; } scalar_t EquilibriumState::logIAP_kinetic(const std::string& label_mineral) const { database::DatabaseModule namer(m_data); int idm = namer.safe_label_to_id(label_mineral, m_data->labels_minerals_kinetic); return logIAP_kinetic(idm); } scalar_t EquilibriumState::saturation_index_kinetic(const std::string& label_mineral) const { database::DatabaseModule namer(m_data); int idm = namer.safe_label_to_id(label_mineral, m_data->labels_minerals); return logIAP_kinetic(idm) - m_data->logk_mineral_kinetic(idm); } scalar_t EquilibriumState::total_aqueous_concentration_component(int i) const { scalar_t totconc = POW10(m_main_variables(i)); for (int j=0; jnb_aqueous; ++j) { if (m_data->nu_aqueous(j,i) != 0) totconc += m_data->nu_aqueous(j,i)*m_concentration_secondary(j); } return totconc; } void EquilibriumState::total_aqueous_concentrations(Eigen::VectorXd& totconc) const { assert(totconc.rows() == m_data->nb_component); totconc(0) = m_main_variables(0); for (int i=1; inb_component; ++i) totconc(i) = total_aqueous_concentration_component(i); } scalar_t EquilibriumState::pH() const { // search where the information about "H[+]" is stored // (i.e. was the component H[+] or HO[-] ?) database::DatabaseModule namer(m_data); scalar_t id = namer.component_label_to_id("H[+]"); if (id != no_species) return -activity_component(id); else { id = namer.aqueous_label_to_id("H[+]"); - assert(id != database::DatabaseModule::no_species); + assert(id != no_species); return -std::log10(activity_secondary(id)); } } scalar_t EquilibriumState::mass_mineral(int m) const { return m_data->molar_mass_mineral(m)*m_main_variables(m_data->nb_component+m); } scalar_t EquilibriumState::condensed_phase_volume() const { scalar_t volume = mass_water()/constants::water_density_25; for (int m=0; mnb_mineral; ++m) { volume += moles_mineral(m)*m_data->molar_volume_mineral(m); } return volume; } scalar_t EquilibriumState::gas_phase_volume() const { assert(is_fixed_volume()); return (get_volume() - condensed_phase_volume()); } scalar_t EquilibriumState::volume_minerals() const { scalar_t volume = 0; for (int m=0; mnb_mineral; ++m) { volume += moles_mineral(m)*m_data->molar_volume_mineral(m); } return volume; } void EquilibriumState::scale_condensed(scalar_t scale) { m_main_variables(0) *= scale; m_main_variables.segment(m_data->nb_component, m_data->nb_mineral) *= scale; } } // end namespace specmicp diff --git a/src/specmicp/reduced_system.cpp b/src/specmicp/reduced_system.cpp index fed4668..bfb5f53 100644 --- a/src/specmicp/reduced_system.cpp +++ b/src/specmicp/reduced_system.cpp @@ -1,525 +1,525 @@ /*------------------------------------------------------- - Module : specmicp - File : specmicp.cpp - Author : Fabien Georget Copyright (c) 2014, Fabien 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: * Redistributions of source code must retain the above copyright notice, this list of conditions and the following disclaimer. * 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. * Neither the name of the Princeton University 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 OWNER 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 #include "reduced_system.hpp" #include "utils/log.hpp" #include "equilibrium_data.hpp" #include "physics/constants.hpp" #include "physics/laws.hpp" namespace specmicp { inline double pow10(double x) { return std::pow(10.0, x); } ReducedSystem::ReducedSystem( RawDatabasePtr ptrdata, const Vector &totconc, bool conservation_water ): m_data(ptrdata), m_tot_conc(totconc), m_secondary_conc(ptrdata->nb_aqueous), m_loggamma(ptrdata->nb_component+ptrdata->nb_aqueous) { if (not m_simulbox.is_fixed_temperature() or m_simulbox.get_temperature() != units::celsius(25)) { throw std::runtime_error("Only a fixed temperature of 25°C is supported"); } number_eq(conservation_water); m_secondary_conc.setZero(); m_loggamma.setZero(); } ReducedSystem::ReducedSystem( RawDatabasePtr ptrdata, const Vector &totconc, const SimulationBox &simul_box, bool conservation_water ): m_simulbox(simul_box), m_data(ptrdata), m_tot_conc(totconc), m_secondary_conc(ptrdata->nb_aqueous), m_loggamma(ptrdata->nb_component+ptrdata->nb_aqueous) { if (not m_simulbox.is_fixed_temperature() or m_simulbox.get_temperature() != units::celsius(25)) { throw std::runtime_error("Only a fixed temperature of 25°C is supported"); } number_eq(conservation_water); m_secondary_conc.setZero(); m_loggamma.setZero(); } ReducedSystem::ReducedSystem( RawDatabasePtr ptrdata, const Vector& totconc, const SimulationBox& simul_box, const EquilibriumState& previous_solution, bool conservation_water ): m_simulbox(simul_box), m_data(ptrdata), m_tot_conc(totconc), m_secondary_conc(previous_solution.molality_secondary()), m_loggamma(previous_solution.activity_coefficient()) { if (not m_simulbox.is_fixed_temperature() or m_simulbox.get_temperature() != units::celsius(25)) { throw std::runtime_error("Only a fixed temperature of 25°C is supported"); } number_eq(conservation_water); } void ReducedSystem::number_eq(bool water_equation) { index_t neq = 0; m_ideq.reserve(m_data->nb_component+m_data->nb_mineral); if (water_equation) { m_ideq.push_back(neq); ++neq; } else { m_ideq.push_back(no_equation); } for (index_t i: m_data->range_aqueous_component()) { // Remove components that does not exist if (m_tot_conc(i) == 0 and i!= 1) { INFO << "Component " << m_data->labels_basis[i] << "has total concentration equal to zero, we desactivate it" << std::endl; m_nonactive_component.push_back(i); m_ideq.push_back(no_equation); continue; } else { m_ideq.push_back(neq); ++neq; } } m_nb_free_vars = neq; for (index_t m: m_data->range_mineral()) { bool can_precipitate = true; // Remove minerals that cannot precipitate for (auto it=m_nonactive_component.begin(); it!=m_nonactive_component.end(); ++it) { if (m_data->nu_mineral(m, *it) != 0) { can_precipitate = false; break; // this is not a mineral that can precipitate } } if (can_precipitate) { m_ideq.push_back(neq); ++neq; } else { m_ideq.push_back(no_equation); } } m_nb_tot_vars = neq; m_nb_compl_vars = m_nb_tot_vars - m_nb_free_vars; // aqueous species m_active_aqueous.reserve(m_data->nb_aqueous); for (index_t j: m_data->range_aqueous()) { bool can_exist = true; for (auto it=m_nonactive_component.begin(); it!=m_nonactive_component.end(); ++it) { if (m_data->nu_aqueous(j,*it) != 0) { can_exist = false; } } m_active_aqueous.push_back(can_exist); } } scalar_t ReducedSystem::residual_water(const Vector& x) const { scalar_t res = m_tot_conc(0) - mass_water(x)/m_data->molar_mass_basis_si(0); for (index_t j: m_data->range_aqueous()) { if (m_data->nu_aqueous(j, 0) == 0) continue; res -= mass_water(x)*m_data->nu_aqueous(j, 0)*m_secondary_conc(j); } for (index_t m: m_data->range_mineral()) { if (ideq_min(m) == no_equation or m_data->nu_mineral(m, 0) == 0) continue; res -= m_data->nu_mineral(m, 0)*x(ideq_min(m)); } return res; } double ReducedSystem::residual_component(const Vector &x, index_t i) const { const index_t idp = ideq_paq(i); scalar_t res = m_tot_conc(i) - mass_water(x)*pow10(x(idp)); for (index_t j: m_data->range_aqueous()) { if (not m_active_aqueous[j]) continue; res -= mass_water(x)*m_data->nu_aqueous(j, i)*m_secondary_conc(j); } for (index_t m: m_data->range_mineral()) { if (ideq_min(m) == no_equation) continue; res -= m_data->nu_mineral(m, i)*x(ideq_min(m)); } return res; } double ReducedSystem::residual_mineral(const Vector& x, index_t m) const { scalar_t res = m_data->logk_mineral(m); for (index_t i: m_data->range_aqueous_component()) { if (m_data->nu_mineral(m, i) != 0) res -= m_data->nu_mineral(m, i)*(x(ideq_paq(i)) + m_loggamma(i)); } return res; } void ReducedSystem::get_residuals(const Vector& x, Vector& residual) { set_secondary_concentration(x); if (ideq_w() != no_equation) residual(ideq_w()) = residual_water(x); for (index_t i: m_data->range_aqueous_component()) { if (ideq_paq(i) != no_equation) residual(ideq_paq(i)) = residual_component(x, i); } for (index_t m: m_data->range_mineral()) { if (ideq_min(m) != no_equation) residual(ideq_min(m)) = residual_mineral(x, m); } } //non-optimized Finite difference, for test only //void ReducedSystem::get_jacobian(Eigen::VectorXd& x, // Eigen::MatrixXd& jacobian) //{ // const int neq = total_variables(); // Eigen::VectorXd res(total_variables()); // Eigen::VectorXd perturbed_res(total_variables()); // get_residuals(x, res); // for (int j=0; jmolar_mass_basis_si(0); for (index_t j: m_data->range_aqueous()) { if ( m_data->nu_aqueous(j, 0) == 0 ) continue; tmp -= m_data->nu_aqueous(j, 0)*m_secondary_conc(j); } jacobian(0, 0) = tmp; for (index_t k: m_data->range_aqueous_component()) { if (ideq_paq(k) == no_equation) continue; scalar_t tmp = 0; for (index_t j: m_data->range_aqueous()) { tmp -= m_data->nu_aqueous(j,0)*m_data->nu_aqueous(j,k)*m_secondary_conc(j); } jacobian(0, ideq_paq(k)) = x(ideq_w())*log10 * tmp; } for (index_t m: m_data->range_mineral()) { if (ideq_min(m) == no_equation) continue; jacobian(0, ideq_min(m)) = -m_data->nu_mineral(m, 0); } } // aqueous component for (index_t i: m_data->range_aqueous_component()) { const index_t idp = ideq_paq(i); if (idp == no_equation) continue; // aqueous components for (index_t k: m_data->range_aqueous_component()) { if (ideq_paq(k) == no_equation) continue; scalar_t tmp_iip = 0; if (k == i) tmp_iip -= pow10(x(ideq_paq(i)))*log10; for (index_t j: m_data->range_aqueous()) { tmp_iip -= log10*m_data->nu_aqueous(j, i)*m_data->nu_aqueous(j, k)*m_secondary_conc(j); } jacobian(idp, ideq_paq(k)) = mass_water(x)*tmp_iip; } // minerals for (index_t m: m_data->range_mineral()) { if (ideq_min(m) == no_equation) continue; jacobian(idp, ideq_min(m)) = - m_data->nu_mineral(m, i); } if (ideq_w() != no_equation) // Water { scalar_t tmp_iw = -std::pow(10.0, x(idp)); for (index_t j: m_data->range_aqueous()) { if ( m_data->nu_aqueous(j, i) == 0 ) continue; tmp_iw -= m_data->nu_aqueous(j, i)*m_secondary_conc(j); } jacobian(idp, ideq_w()) = tmp_iw; } } // mineral equilibrium for (index_t m: m_data->range_mineral()) { const index_t idm = ideq_min(m); if (idm == no_equation) continue; for (index_t i: m_data->range_aqueous_component()) { if (ideq_paq(i) == no_equation) continue; jacobian(idm, ideq_paq(i)) = -m_data->nu_mineral(m, i); } } } void ReducedSystem::set_secondary_concentration(const Vector& x) { for (index_t j: m_data->range_aqueous()) { if (m_active_aqueous[j] == false) { m_secondary_conc(j) = 0; continue; } scalar_t conc = -m_data->logk_aqueous(j) - m_loggamma(j+m_data->nb_component); for (index_t k: m_data->range_aqueous_component()) { if (m_data->nu_aqueous(j, k) == 0) continue; conc += m_data->nu_aqueous(j, k)*(x(ideq_paq(k)) + m_loggamma(k)); } m_secondary_conc(j) = pow10(conc); } } void ReducedSystem::set_ionic_strength(const Vector& x) { scalar_t ionic = 0; for (index_t i: m_data->range_aqueous_component()) { if (ideq_paq(i) == no_equation or m_data->charge_component(i) == 0) continue; ionic += pow10(x(ideq_paq(i)))*std::pow(m_data->charge_component(i),2); } for (index_t j: m_data->range_aqueous()) { if (not m_active_aqueous[j] or m_data->charge_aqueous(j) == 0) continue; ionic += m_secondary_conc(j)*std::pow(m_data->charge_aqueous(j),2); } m_ionic_strength = ionic/2; } void ReducedSystem::compute_log_gamma(const Vector& x) { set_ionic_strength(x); const scalar_t sqrti = std::sqrt(m_ionic_strength); for (index_t i: m_data->range_aqueous_component()) { if (ideq_paq(i) == no_equation) { m_loggamma(i) = 0; continue; } m_loggamma(i) = laws::extended_debye_huckel(m_ionic_strength, sqrti, m_data->charge_component(i), m_data->a_debye_component(i), m_data->b_debye_component(i)); } for (index_t j: m_data->range_aqueous()) { if (not m_active_aqueous[j]) { m_loggamma(m_data->nb_component + j) = 0; continue; } m_loggamma(m_data->nb_component + j) = laws::extended_debye_huckel( m_ionic_strength, sqrti, m_data->charge_aqueous(j), m_data->a_debye_aqueous(j), m_data->b_debye_aqueous(j)); } } bool ReducedSystem::hook_start_iteration(const Vector& x, scalar_t norm_residual) { not_in_linesearch = true; scalar_t previous_norm = m_loggamma.norm(); bool may_have_converged = false; if (norm_residual < nb_free_variables()) { // Use fix point iterations for non-ideality - for (int i=0; i<6; ++i) + for (int i=0; i<10; ++i) { set_secondary_concentration(x); compute_log_gamma(x); - if (std::abs(previous_norm - m_loggamma.norm()) < 1e-6) {may_have_converged = true; break;} + if (std::abs(previous_norm - m_loggamma.norm())/previous_norm < 1e-10) {may_have_converged = true; break;} previous_norm = m_loggamma.norm(); } } return may_have_converged; } double ReducedSystem::max_lambda(const Vector& x, const Vector& update) { if (ideq_w() != no_equation) { return 1.0/std::max(1.0, -update(0)/(0.9*x(0))); } else { return 1.0; } } void ReducedSystem::reasonable_starting_guess(Vector &x, bool for_min) { x.resize(m_data->nb_component+ m_data->nb_mineral); x(0) = 2*m_tot_conc(0)*m_data->molar_mass_basis_si(0); for (index_t i: m_data->range_aqueous_component()) { x(i) = -4.0; } if (for_min) x.block(m_data->nb_component, 0, m_data->nb_mineral, 1).setZero(); else x.block(m_data->nb_component, 0, m_data->nb_mineral, 1).setZero(); m_loggamma.setZero(); m_secondary_conc.setZero(); } void ReducedSystem::reasonable_restarting_guess(Vector& x) { x(0) = m_tot_conc(0)*m_data->molar_mass_basis_si(0); for (index_t i=m_data->nb_component; inb_component+m_data->nb_mineral; ++i) { x(i) = 0.0; } m_loggamma.setZero(); m_secondary_conc.setZero(); } EquilibriumState ReducedSystem::get_solution(const Vector& xtot, const Vector& x) { double previous_norm = m_loggamma.norm(); set_secondary_concentration(x); compute_log_gamma(x); if (std::abs(previous_norm - m_loggamma.norm()) > 1e-6) { WARNING << "Activity coefficient have not converged !" << std::endl << "output can not be trusted\n Difference "+std::to_string(std::abs(previous_norm - m_loggamma.norm())); } return EquilibriumState(xtot, m_secondary_conc, m_loggamma, m_ionic_strength, m_data); } void ReducedSystem::reduce_mineral_problem(Vector& true_var) { for (index_t mineral: m_data->range_mineral()) { if (ideq_min(mineral) == no_equation) continue; if (m_data->stability_mineral(mineral) == database::MineralStabilityClass::cannot_dissolve) { if (true_var(ideq_min(mineral)) == 0.0) continue; for (index_t component: m_data->range_component()) { if (m_data->nu_mineral(mineral, component) == 0.0) continue; m_tot_conc(component) -= m_data->nu_mineral(mineral, component)*true_var(ideq_min(mineral)); } true_var(ideq_min(mineral)) = 0.0; } } } void ReducedSystem::reset_mineral_system(Vector& true_var, Vector& output_var) { for (index_t mineral: m_data->range_mineral()) { if (ideq_min(mineral) == no_equation) continue; if (m_data->stability_mineral(mineral) == database::MineralStabilityClass::cannot_dissolve) { output_var(mineral) += true_var(ideq_min(mineral)); } } } } // end namespace specmicp diff --git a/tests/reactmicp/systems/saturated_diffusion/neutrality_solver.cpp b/tests/reactmicp/systems/saturated_diffusion/neutrality_solver.cpp index 579dc5a..3b5cc3b 100644 --- a/tests/reactmicp/systems/saturated_diffusion/neutrality_solver.cpp +++ b/tests/reactmicp/systems/saturated_diffusion/neutrality_solver.cpp @@ -1,252 +1,277 @@ #include "catch.hpp" #include "utils.hpp" #include "dfpmsolver/parabolic_driver.hpp" #include "reactmicp/meshes/uniform_mesh1d.hpp" #include "reactmicp/systems/saturated_diffusion/transport_neutrality_program.hpp" #include "reactmicp/systems/saturated_diffusion/transport_neutrality_parameters.hpp" using namespace specmicp; using namespace specmicp::reactmicp; using namespace specmicp::reactmicp::systems; using namespace specmicp::reactmicp::systems::siasaturated; -#define NB_STEP 100 +#define NB_STEP 20 #define DT 1e4 TEST_CASE("dfpm solver for neutrality program", "[dfpm, solver, neutrality]") { SECTION("transport neutrality program") { std::shared_ptr thedatabase = get_test_carbo_database(); - int nb_nodes = 5; + int nb_nodes = 10; EquilibriumState initial_state = sample_carbo_composition(thedatabase); mesh::Mesh1DPtr themesh = mesh::uniform_mesh1d(nb_nodes, 0.1, 2.5); //auto parameters = std::make_shared(nb_nodes, 1e-8, 0.2); std::shared_ptr parameters = std::make_shared( nb_nodes, // nb_nodes, thedatabase->nb_component, // nb_component, thedatabase->nb_aqueous, // nb_aqueous, 2e-5, // the_diffusion_coefficient, 0.25, // porosity 1.0e-3 //reduction_factor ); //SIASaturatedVariables variables(database, themesh, parameters); SaturatedNeutralityDiffusionProgram the_program(themesh, thedatabase, parameters); Vector log_concentrations(the_program.get_ndf()*themesh->nb_nodes()); for (index_t node: themesh->range_nodes()) { for (index_t component: thedatabase->range_aqueous_component()) { log_concentrations(component-1+the_program.get_ndf()*node) = std::log10(initial_state.molality_component(component)); } } // EquilibriumState bc_state = blank_composition(thedatabase); // for (index_t component: thedatabase->range_aqueous_component()) // { // log_concentrations(component-1) = // std::log10(bc_state.molality_component(component)); // } log_concentrations.block(0, 0, the_program.get_ndf(), 1).setConstant(-9); TransportNeutralityVariables var0(nb_nodes, thedatabase); var0.component_concentrations() = log_concentrations; var0.solve_secondary_variables(); for (index_t node=1; node solver(the_program); solver.get_velocity() = Vector::Zero(log_concentrations.rows()); solver.get_options().maximum_iterations = 100; - solver.get_options().residuals_tolerance = 1e-5; - solver.get_options().quasi_newton = 4; + solver.get_options().residuals_tolerance = 1e-7; + solver.get_options().step_tolerance = 1e-6; + solver.get_options().quasi_newton = 1; solver.get_options().linesearch = dfpmsolver::ParabolicLinesearch::Strang; - + solver.set_scaling(Vector::Constant(the_program.get_neq(), 1e6)); std::cout << log_concentrations << std::endl; for (int k=0; k thedatabase = get_test_carbo_database(); database::Database dbhandler(thedatabase); - int nb_nodes = 5; + int nb_nodes = 10; EquilibriumState initial_state = sample_carbo_composition(thedatabase); mesh::Mesh1DPtr themesh = mesh::uniform_mesh1d(nb_nodes, 0.1, 2.5); //auto parameters = std::make_shared(nb_nodes, 1e-8, 0.2); std::shared_ptr parameters = std::make_shared( nb_nodes, // nb_nodes, thedatabase->nb_component, // nb_component, thedatabase->nb_aqueous, // nb_aqueous, 2e-5, // the_diffusion_coefficient, 0.25, // porosity 0.5e-3 //reduction_factor ); parameters->diff_coeff_secondary(dbhandler.aqueous_label_to_id("Si(OH)4")) = 1e-6; //SIASaturatedVariables variables(database, themesh, parameters); SaturatedNeutralityDiffusionProgram the_program(themesh, thedatabase, parameters); Vector log_concentrations(the_program.get_ndf()*themesh->nb_nodes()); for (index_t node: themesh->range_nodes()) { for (index_t component: thedatabase->range_aqueous_component()) { log_concentrations(component-1+the_program.get_ndf()*node) = std::log10(initial_state.molality_component(component)); } } log_concentrations.block(0, 0, the_program.get_ndf(), 1).setConstant(-9); the_program.number_equations({0,}); dfpmsolver::ParabolicDriver solver(the_program); solver.get_velocity() = Vector::Zero(log_concentrations.rows()); solver.get_options().maximum_iterations = 200; + solver.get_options().residuals_tolerance = 1e-6; + solver.get_options().step_tolerance = 1e-6; + solver.set_scaling(Vector::Constant(the_program.get_neq(), 1e6)); + + solver.get_options().quasi_newton = 1; + //solver.get_options().linesearch = dfpmsolver::ParabolicLinesearch::Strang; std::cout << log_concentrations << std::endl; for (int k=0; k thedatabase = get_test_carbo_database(); database::Database dbhandler(thedatabase); - int nb_nodes = 5; + int nb_nodes = 10; EquilibriumState initial_state = sample_carbo_composition(thedatabase); mesh::Mesh1DPtr themesh = mesh::uniform_mesh1d(nb_nodes, 0.1, 2.5); //auto parameters = std::make_shared(nb_nodes, 1e-8, 0.2); std::shared_ptr parameters = std::make_shared( nb_nodes, // nb_nodes, thedatabase->nb_component, // nb_component, thedatabase->nb_aqueous, // nb_aqueous, 2e-5, // the_diffusion_coefficient, 0.25, // porosity 0.5e-3 //reduction_factor ); parameters->diff_coeff_secondary(dbhandler.aqueous_label_to_id("H[+]")) = 9e-5; parameters->diff_coeff_component(dbhandler.component_label_to_id("HO[-]")) = 9e-5; parameters->diff_coeff_secondary(dbhandler.aqueous_label_to_id("Si(OH)4")) = 1e-6; parameters->diff_coeff_component(dbhandler.component_label_to_id("SiO(OH)3[-]")) = 1e-6; //SIASaturatedVariables variables(database, themesh, parameters); SaturatedNeutralityDiffusionProgram the_program(themesh, thedatabase, parameters); Vector log_concentrations(the_program.get_ndf()*themesh->nb_nodes()); for (index_t node: themesh->range_nodes()) { for (index_t component: thedatabase->range_aqueous_component()) { log_concentrations(component-1+the_program.get_ndf()*node) = std::log10(initial_state.molality_component(component)); } } log_concentrations.block(0, 0, the_program.get_ndf(), 1).setConstant(-9); the_program.number_equations({0,}); dfpmsolver::ParabolicDriver solver(the_program); solver.get_velocity() = Vector::Zero(log_concentrations.rows()); - solver.get_options().maximum_iterations = 200; + solver.get_options().maximum_iterations = 100; solver.get_options().residuals_tolerance = 1e-6; + solver.set_scaling(Vector::Constant(the_program.get_neq(), 1e6)); + + solver.get_options().quasi_newton = 1; + solver.get_options().linesearch = dfpmsolver::ParabolicLinesearch::Strang; std::cout << log_concentrations << std::endl; for (int k=0; k #include "utils.hpp" #include "reactmicp/meshes/uniform_mesh1d.hpp" #include "reactmicp/systems/saturated_diffusion/boundary_conditions.hpp" #include "reactmicp/systems/saturated_diffusion/reactive_transport_solver.hpp" +#include "reactmicp/systems/saturated_diffusion/transport_parameters.hpp" #include using namespace specmicp; using namespace specmicp::reactmicp; using namespace specmicp::reactmicp::systems; using namespace specmicp::reactmicp::systems::siasaturated; TEST_CASE("reactive transport solver", "[reactive transport, speciation, dfpm, solver]") { specmicp::stdlog::ReportLevel() = specmicp::logger::Warning; specmicp::logger::ErrFile::stream() = &std::cerr; int nb_nodes = 10; std::shared_ptr database = get_test_carbo_database(); mesh::Mesh1DPtr themesh = mesh::uniform_mesh1d(nb_nodes, 0.1, 2.5); auto parameters = std::make_shared(nb_nodes, 1e-4, 0.2); EquilibriumState initial_state = sample_carbo_composition(database, 0.02); // scaling { double volume = initial_state.volume_minerals(); volume /= (1 - parameters->porosity(1) ); double scale = themesh->get_volume_cell(1)*1e-6 / volume; initial_state.scale_condensed(scale); // } SIABoundaryConditions bcs(nb_nodes); bcs.list_initial_states.push_back(sample_carbo_composition(database)); bcs.list_initial_states.push_back(blank_composition(database)); bcs.bs_types[0] = BoundaryConditionType::FixedComposition; bcs.initial_states[0] = 1; SECTION("Initialization") { SIASaturatedReactiveTransportSolver solver(themesh, database, parameters); solver.apply_boundary_conditions(bcs); SIASaturatedVariables& variables = solver.get_variables(); REQUIRE(std::abs(variables.mineral_amount(1, 1) - bcs.list_initial_states[0].moles_mineral(1)) < 1e-10 ); REQUIRE(std::abs(variables.component_concentration(1, 1) - bcs.list_initial_states[0].molality_component(1)) < 1e-10 ); REQUIRE(std::abs(variables.mineral_amount(0, 1) - bcs.list_initial_states[1].moles_mineral(1)) < 1e-10 ); REQUIRE(std::abs(variables.component_concentration(0, 1) - bcs.list_initial_states[1].molality_component(1)) < 1e-10 ); } SECTION("Solving") { std::ofstream output; output.open("out.dat"); SIASaturatedReactiveTransportSolver solver(themesh, database, parameters); solver.apply_boundary_conditions(bcs); for (auto it=database->labels_basis.begin(); it!=database->labels_basis.end(); ++it) { std::cout << *it << " "; } std::cout << std::endl; std::cout << database->nu_mineral << std::endl; solver.use_sia(50); + //solver.get_options().transport_options.linesearch = dfpmsolver::ParabolicLinesearch::Strang; + solver.get_options().transport_options.quasi_newton = 2; + solver.get_options().transport_options.step_tolerance = 1e-10; double dt = 100; double total=0 ; - for (int k=0; k<10000; ++k) + for (int k=0; k<5; ++k) { std::cout << " iterations : " << k << std::endl; SIASaturatedReactiveTransportSolverReturnCode retcode = solver.solve_timestep(dt); - REQUIRE(retcode == SIASaturatedReactiveTransportSolverReturnCode::Success); + REQUIRE((int) retcode == (int) SIASaturatedReactiveTransportSolverReturnCode::Success); SIASaturatedVariables& variables = solver.get_variables(); output << total << " " << total/(3600.0*24.0); for(int node: themesh->range_nodes()) { output << " " << variables.nodal_component_update_total_amount(node, 3); } output << std::endl; total+=dt; } output.close(); } } diff --git a/tests/reactmicp/systems/saturated_diffusion/solver.cpp b/tests/reactmicp/systems/saturated_diffusion/solver.cpp index 2bb898b..e9c75e8 100644 --- a/tests/reactmicp/systems/saturated_diffusion/solver.cpp +++ b/tests/reactmicp/systems/saturated_diffusion/solver.cpp @@ -1,91 +1,91 @@ #include "catch.hpp" #include "utils.hpp" #include "dfpmsolver/parabolic_driver.hpp" #include "reactmicp/meshes/uniform_mesh1d.hpp" #include "reactmicp/systems/saturated_diffusion/transport_program.hpp" #include "reactmicp/systems/saturated_diffusion/variables.hpp" - +#include "reactmicp/systems/saturated_diffusion/transport_parameters.hpp" #include "specmicp/reduced_system_solver.hpp" using namespace specmicp; using namespace specmicp::reactmicp; using namespace specmicp::reactmicp::systems; using namespace specmicp::reactmicp::systems::siasaturated; TEST_CASE("dfpm solver", "[dfpm, solver]") { SECTION("transport program") { std::shared_ptr database = get_test_carbo_database(); int nb_nodes = 5; EquilibriumState initial_state = sample_carbo_composition(database); mesh::Mesh1DPtr themesh = mesh::uniform_mesh1d(nb_nodes, 1e-3, 0.001); auto parameters = std::make_shared(nb_nodes, 1e-8, 0.2); SIASaturatedVariables variables(database, themesh, parameters); SaturatedDiffusionProgram the_program(themesh, database, parameters); Eigen::VectorXd& concentrations = variables.total_mobile_concentrations(); concentrations.setOnes(); concentrations.block(0, 0, the_program.get_ndf(), 1).setZero(); the_program.number_equations({0,}); dfpmsolver::ParabolicDriver solver(the_program); solver.get_velocity() = Eigen::VectorXd::Zero(concentrations.rows()); double dt = 100; dfpmsolver::ParabolicDriverReturnCode retcode = solver.solve_timestep(dt, concentrations); REQUIRE(retcode == dfpmsolver::ParabolicDriverReturnCode::ResidualMinimized); for (int node=1; nodenb_nodes(); ++node) { for (int component: database->range_aqueous_component()) { the_program.external_flux(node, component) = parameters->porosity(0)*(1.0 - concentrations(the_program.get_dof_component(node, component))) / (parameters->density_water()*dt*themesh->get_volume_cell(node)); } } retcode = solver.restart_timestep(concentrations); REQUIRE(retcode == dfpmsolver::ParabolicDriverReturnCode::ResidualMinimized); } SECTION("micpsolver") { std::shared_ptr database = get_test_carbo_database(); database::Database database_handler(database); int nb_nodes = 5; mesh::Mesh1DPtr themesh = mesh::uniform_mesh1d(nb_nodes, 1e-3, 0.001); auto parameters = std::make_shared(nb_nodes, 1e-8, 0.2); SIASaturatedVariables variables(database, themesh, parameters); EquilibriumState initial_state = sample_carbo_composition(database); variables.update_composition(1, initial_state); Eigen::VectorXd totconc; variables.nodal_update_total_amount(1, totconc); totconc(database_handler.component_label_to_id("HCO3[-]")) += 0.05; totconc(database_handler.component_label_to_id("HO[-]")) -= 0.05; ReducedSystemSolverOptions options; options.solver_options.max_factorization_step = 1; ReducedSystemSolver solver(database, totconc, options); Eigen::VectorXd spec_var(variables.speciation_variables().col(1)); specmicp::micpsolver::MiCPPerformance perf = solver.solve(spec_var); REQUIRE((int) perf.return_code == (int) specmicp::micpsolver::MiCPSolverReturnCode::ResidualMinimized); EquilibriumState final_compo = solver.get_solution(spec_var); variables.update_composition(1, final_compo); } } diff --git a/tests/reactmicp/systems/saturated_diffusion/transport_program.cpp b/tests/reactmicp/systems/saturated_diffusion/transport_program.cpp index b7697dc..f3f9b25 100644 --- a/tests/reactmicp/systems/saturated_diffusion/transport_program.cpp +++ b/tests/reactmicp/systems/saturated_diffusion/transport_program.cpp @@ -1,27 +1,28 @@ #include "catch.hpp" #include "utils.hpp" #include "reactmicp/meshes/uniform_mesh1d.hpp" #include "reactmicp/systems/saturated_diffusion/transport_program.hpp" +#include "reactmicp/systems/saturated_diffusion/transport_parameters.hpp" using namespace specmicp; using namespace specmicp::reactmicp; using namespace specmicp::reactmicp::systems; using namespace specmicp::reactmicp::systems::siasaturated; TEST_CASE("Transport program", "[Transport, Diffusion, Finite Element, Solver]") { std::shared_ptr database = get_test_carbo_database(); int nb_nodes = 5; EquilibriumState initial_state = sample_carbo_composition(database); mesh::Mesh1DPtr themesh = mesh::uniform_mesh1d(nb_nodes, 1e-3, 0.001); auto parameters = std::make_shared(nb_nodes, 1e-6, 0.2); SECTION("Initialization") { SaturatedDiffusionProgram(themesh, database, parameters); } } diff --git a/tests/reactmicp/systems/saturated_diffusion/variables.cpp b/tests/reactmicp/systems/saturated_diffusion/variables.cpp index 1d2845d..e2ceb34 100644 --- a/tests/reactmicp/systems/saturated_diffusion/variables.cpp +++ b/tests/reactmicp/systems/saturated_diffusion/variables.cpp @@ -1,64 +1,65 @@ #include "catch.hpp" #include "utils.hpp" #include "reactmicp/systems/saturated_diffusion/variables.hpp" +#include "reactmicp/systems/saturated_diffusion/transport_parameters.hpp" #include "reactmicp/meshes/uniform_mesh1d.hpp" using namespace specmicp; using namespace specmicp::reactmicp; using namespace specmicp::reactmicp::systems::siasaturated; TEST_CASE("SIA Variables", "[Variables, Secondary variables, SIA]") { std::shared_ptr database = get_test_carbo_database(); int nb_nodes = 5; mesh::Mesh1DPtr themesh = mesh::uniform_mesh1d(nb_nodes, 1e-3, 0.001); auto parameters = std::make_shared(nb_nodes, 1e-6, 0.2); SECTION("Initialization") { SIASaturatedVariables variables(database, themesh, parameters); REQUIRE_NOTHROW( variables.component_concentration(0, 2) ); REQUIRE_NOTHROW( variables.component_concentration(2, 2) ); REQUIRE_NOTHROW( variables.component_concentration(4, 4) ); REQUIRE_NOTHROW( variables.component_concentration(3, 1) ); REQUIRE_NOTHROW( variables.log_component_concentration(3, 3) ); REQUIRE_NOTHROW( variables.log_gamma_component(4, 2) ); } SECTION("Ionic strength computation") { SIASaturatedVariables variables(database, themesh, parameters); database::Database database_handler(database); variables.log_component_concentration(1, database_handler.component_label_to_id("Ca[2+]")) = -3; variables.log_component_concentration(1, database_handler.component_label_to_id("HO[-]")) = -3; variables.log_component_concentration(1, database_handler.component_label_to_id("HCO3[-]")) = -3; variables.log_component_concentration(1, database_handler.component_label_to_id("SiO(OH)3[-]")) = -3; variables.nodal_ionic_strength(1); REQUIRE( variables.ionic_strength(1) == (4*1e-3 +3*1e-3)*0.5 ); variables.secondary_concentration(1, database_handler.aqueous_label_to_id("CO3[2-]")) = 1e-2; variables.nodal_ionic_strength(1); REQUIRE( variables.ionic_strength(1) == (4*1e-3 +3*1e-3 + 4*1e-2)*0.5 ); variables.secondary_concentration(1, database_handler.aqueous_label_to_id("CO2")) = 1e-1; variables.nodal_ionic_strength(1); REQUIRE( variables.ionic_strength(1) == (4*1e-3 +3*1e-3 + 4*1e-2)*0.5 ); } SECTION("Update") { SIASaturatedVariables variables(database, themesh, parameters); EquilibriumState composition = sample_carbo_composition(database); variables.update_composition(1, composition); - REQUIRE(variables.component_concentration(1, 1) == composition.molality_component(1)); - REQUIRE(variables.total_mobile_concentration(1, 1) == composition.total_aqueous_concentration_component(1)); + REQUIRE(variables.component_concentration(1, 1) == Approx(composition.molality_component(1)).epsilon(1e-7)); + REQUIRE(variables.total_mobile_concentration(1, 1) == Approx(composition.total_aqueous_concentration_component(1)).epsilon(1e-7)); } } diff --git a/tests/specmicp/carboalu.cpp b/tests/specmicp/carboalu.cpp index 7efb68f..008ef1a 100644 --- a/tests/specmicp/carboalu.cpp +++ b/tests/specmicp/carboalu.cpp @@ -1,199 +1,202 @@ /*------------------------------------------------------- - Module : tests/ - File : carboalu.cpp - Author : Fabien Georget Copyright (c) 2014, Fabien 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: * Redistributions of source code must retain the above copyright notice, this list of conditions and the following disclaimer. * 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. * Neither the name of the Princeton University 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 OWNER 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 #include "specmicp/reaction_path.hpp" #include "utils/log.hpp" #include double mult = 1.0; double m_c3s = mult*0.6; double m_c2s = mult*0.2; double m_c3a = mult*0.10; double m_gypsum = mult*0.10; double wc = 0.8; class Timer { public: Timer() : beg_(clock_::now()) {} void reset() { beg_ = clock_::now(); } double elapsed() const { return std::chrono::duration_cast (clock_::now() - beg_).count(); } private: typedef std::chrono::high_resolution_clock clock_; typedef std::chrono::duration > second_; std::chrono::time_point beg_; }; struct Params { specmicp::micpsolver::NCPfunction ncpfunction; double penalizationfactor; bool scaling; }; const Params CCKParam = {specmicp::micpsolver::NCPfunction::penalizedFB, 0.8, false}; const Params FBParam = {specmicp::micpsolver::NCPfunction::penalizedFB, 1.0, false}; const Params minParam = {specmicp::micpsolver::NCPfunction::min, 1.0, true}; void solve_carboalu(double delta, Params param, bool output=true, bool coldstart=false) { - specmicp::database::Database database("data/cemdata_specmicp.js"); + specmicp::database::Database database("../data/cemdata_specmicp.js"); std::shared_ptr data = database.get_database(); + std::map swapping ({ {"H[+]","HO[-]"}, {"Al[3+]","Al(OH)4[-]"}, {"Si(OH)4", "SiO(OH)3[-]"} }); database.swap_components(swapping); std::shared_ptr model = std::make_shared(); double delta_h2co3 =delta; double delta_sio2 = 0.00; double m_water = wc*1e-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)); //std::cout << m_water << std::endl; model->amount_aqueous = { {"H2O", specmicp::reaction_amount_t(m_water/data->molar_mass_basis_si(0), delta_h2co3+delta_sio2)}, {"HCO3[-]", specmicp::reaction_amount_t(0, delta_h2co3)}, {"HO[-]", specmicp::reaction_amount_t(0, -delta_h2co3-delta_sio2)}, {"SiO(OH)3[-]", specmicp::reaction_amount_t(0, delta_sio2)}, }; model->amount_minerals = { {"C3S", specmicp::reaction_amount_t(m_c3s, 0)}, {"C2S", specmicp::reaction_amount_t(m_c2s, 0)}, {"C3A", specmicp::reaction_amount_t(m_c3a, 0)}, {"Gypsum", specmicp::reaction_amount_t(m_gypsum, 0)} }; model->database_path = "data/cemdata_specmicp.js"; model->minerals_to_keep = { "Portlandite", "CSH,jennite", "CSH,tobermorite", "SiO2,am", "Calcite", "Al(OH)3,am", "Monosulfoaluminate", "Tricarboaluminate", "Monocarboaluminate", "Hemicarboaluminate", "Straetlingite", "Gypsum", "Ettringite", "Thaumasite" }; Eigen::VectorXd x; specmicp::ReactionPathDriver driver(model, data, x); Eigen::VectorXd totaq(data->nb_component); double totiter = 0; double totfact = 0; driver.get_options().solver_options.penalization_factor = param.penalizationfactor; driver.get_options().ncp_function = param.ncpfunction; driver.get_options().solver_options.use_scaling = param.scaling; - driver.get_options().solver_options.max_factorization_step = 2; + driver.get_options().solver_options.max_factorization_step = 4; driver.get_options().solver_options.factor_gradient_search_direction = 400; driver.get_options().solver_options.max_iter = 50; driver.get_options().solver_options.maxstep = 50; driver.get_options().allow_restart = true; if (output) { std::cout << "Reaction_path return_code nb_iter nb_fact pH"; for (auto it = data->labels_basis.begin(); it!= data->labels_basis.end(); ++it) { std::cout << " " << *it; } for (auto it = data->labels_minerals.begin(); it!= data->labels_minerals.end(); ++it) { std::cout << " " << *it; } std::cout << std::endl; } const int nb_step = 1 + 2.5/delta; for (int i=0; inb_component-1,1).transpose(); for (int m=0; m < data->nb_mineral; ++m) std::cout << " " << solution.moles_mineral(m); std::cout << std::endl; } totiter += perf.nb_iterations; totfact += perf.nb_factorization; + + //data->stability_mineral(database.mineral_label_to_id("SiO2,am")) = specmicp::database::MineralStabilityClass::cannot_dissolve; } if (output) { std::cout << "Average iterations : " << totiter/nb_step << std::endl; std::cout << "Average factorization : " << totfact/nb_step << std::endl; } } int main() { specmicp::stdlog::ReportLevel() = specmicp::logger::Warning; specmicp::logger::ErrFile::stream() = &std::cerr; int nbloop = 7500; double delta = 0.001; Timer tmr; - //solve_carboalu(delta, CCKParam, true); - for (int i=0; i