diff --git a/src/dfpmsolver/parabolic_driver.hpp b/src/dfpmsolver/parabolic_driver.hpp index c7e7a1d..2b5fd80 100644 --- a/src/dfpmsolver/parabolic_driver.hpp +++ b/src/dfpmsolver/parabolic_driver.hpp @@ -1,143 +1,147 @@ #ifndef SPECMICP_DFPMSOLVER_PARABOLICDRIVER_HPP #define SPECMICP_DFPMSOLVER_PARABOLICDRIVER_HPP #include "common.hpp" #include #include "driver.hpp" #include "parabolic_structs.hpp" namespace specmicp { namespace dfpmsolver { //! \brief The parabolic driver //! //! Parabolic driver for finite element programs template class ParabolicDriver: public Driver { public: using base=Driver; using base::program; using base::get_neq; using base::get_options; using base::get_perfs; using base::scaling; using base::initialize_scaling; ParabolicDriver(Program& the_program): Driver(the_program), m_velocity(Vector::Zero(the_program.get_tot_ndf())) {} //! \brief Solve a timestep of length dt ParabolicDriverReturnCode solve_timestep(scalar_t dt, Vector& displacement); //! \brief Restart the current timestep ParabolicDriverReturnCode restart_timestep(Vector& displacement); //! \brief Check if the solution has been found ParabolicDriverReturnCode check_convergence(); //! \brief Compute the residuals void compute_residuals(const Vector& displacement, Vector& residual) { compute_residuals(displacement, m_velocity, residual); } //! \brief Compute the residuals void compute_residuals(const Vector& displacement, const Vector& velocity, Vector& residual) { program().compute_residuals(displacement, velocity, residual); get_perfs().nb_call_residuals += 1; } //! \brief Compute the jacobian void compute_jacobian(Vector& displacement, Vector& velocity, Eigen::SparseMatrix& jacobian ); //! \brief Return the norm of the current residuals double norm_residuals() {return m_residuals.norm();} //! \brief Read/Write reference to the velocity vector const Vector& get_velocity() const {return m_velocity;} Vector& velocity() {return m_velocity;} void set_velocity(Vector& velocity_vector); -private: - void reset_velocity(); + + //! \brief Initialize the computation + void initialize_timestep(scalar_t dt, Eigen::VectorXd& displacement); //! \brief Strang Linesearch //! //! ref : //! - Matthies et al. (1979) //! - JHP course notes ParabolicLinesearchReturnCode strang_linesearch( Vector& update, Vector& displacements, scalar_t& lambda ); + scalar_t residuals_0() {return m_norm_0;} +private: + + void reset_velocity(); + //! \brief Backtracking Linesearch //! //! ref : //! - Algo A6.3.1 : Dennis and Schnabel (1983) //! - Nocedal & Wrigth (2006) ParabolicLinesearchReturnCode backtracking_linesearch( Vector& update, Vector& displacements, scalar_t& lambda_out ); //! Update the variables in displacement void update_variable( const Vector& update, scalar_t lambda, Vector& displacement ); //! \brief Set the predictor void set_predictor(Vector& displacement); //! \brief solve the linear system SparseSolverReturnCode solve_linear_system(Eigen::VectorXd& update); //! \brief perform the linesearch ParabolicDriverReturnCode linesearch( Vector &update, Vector &displacements ); //! Compute the residuals for the linesearch double compute_residuals_linesearch( Vector& update, scalar_t lambda, Vector& displacement ); double compute_residuals_strang_linesearch( Vector &update, scalar_t lambda, Vector &displacement ); scalar_t update_norm(const Vector& update); - //! \brief Initialize the computation - void initialize_timestep(scalar_t dt, Eigen::VectorXd& displacement); double m_norm_0; double m_current_dt; Eigen::VectorXd m_gradient; Eigen::VectorXd m_velocity; Eigen::VectorXd m_residuals; Eigen::VectorXd m_predictor; Eigen::SparseMatrix m_jacobian; bool m_velocity_is_initialized; }; } // end namespace dfpmsolver } // end namespace specmicp // implementation // ============== #include "parabolic_driver.inl" #endif // SPECMICP_DFPMSOLVER_PARABOLICDRIVER_HPP diff --git a/src/dfpmsolver/parabolic_driver.inl b/src/dfpmsolver/parabolic_driver.inl index 23d09f9..6f629c8 100644 --- a/src/dfpmsolver/parabolic_driver.inl +++ b/src/dfpmsolver/parabolic_driver.inl @@ -1,471 +1,471 @@ #include "parabolic_driver.hpp" // for syntaxic coloration only #include "utils/log.hpp" #include "utils/sparse_solver.hpp" namespace specmicp { namespace dfpmsolver { template void ParabolicDriver::compute_jacobian(Vector& displacement, Vector& velocity, Eigen::SparseMatrix& jacobian ) { program().compute_jacobian(displacement, velocity, jacobian, get_options().alpha*m_current_dt); jacobian = jacobian*scaling().asDiagonal(); get_perfs().nb_call_jacobian += 1; } template ParabolicDriverReturnCode ParabolicDriver::solve_timestep(scalar_t dt, Eigen::VectorXd& displacement) { initialize_timestep(dt, displacement); - compute_residuals(displacement, m_velocity, m_residuals); - m_norm_0 = norm_residuals(); - get_perfs().nb_iterations = 0; return restart_timestep(displacement); } template void ParabolicDriver::initialize_timestep(scalar_t dt, Eigen::VectorXd& displacement) { initialize_scaling(); m_residuals = Eigen::VectorXd::Zero(get_neq()); m_current_dt = dt; set_predictor(displacement); reset_velocity(); program().apply_bc(dt, displacement, m_velocity); + + compute_residuals(displacement, m_velocity, m_residuals); + m_norm_0 = norm_residuals(); + get_perfs().nb_iterations = 0; } template void ParabolicDriver::set_predictor(Vector& displacement) { if (get_options().alpha < 1) m_predictor = displacement + (1-get_options().alpha)*m_current_dt*m_velocity; else m_predictor = displacement; } template scalar_t ParabolicDriver::update_norm(const Vector& update) { // l-∞ scaled norm scalar_t norm = 0.0; for (index_t dof=0; dof ParabolicDriverReturnCode ParabolicDriver::restart_timestep(Vector& displacement) { ParabolicDriverReturnCode return_code = ParabolicDriverReturnCode::NotConvergedYet; Eigen::VectorXd update(get_neq()); update.setZero(); get_perfs().current_update = 0; while (return_code == ParabolicDriverReturnCode::NotConvergedYet) { compute_residuals(displacement, m_velocity, m_residuals); get_perfs().current_residual = m_residuals.norm()/m_norm_0; get_perfs().current_update = update_norm(update); - //std::cout << " NB iterations : " << get_perfs().nb_iterations - // << " - res : " << get_perfs().current_residual - // << " - update : " << get_perfs().current_update - // << std::endl; + DEBUG << " NB iterations : " << get_perfs().nb_iterations + << " - res : " << get_perfs().current_residual + << " - update : " << get_perfs().current_update; return_code = check_convergence(); if (return_code != ParabolicDriverReturnCode::NotConvergedYet) break; if (get_perfs().nb_iterations % get_options().quasi_newton == 0) { // Todo : avoid computation of decomposition if it is the same matrix compute_jacobian(displacement, m_velocity, m_jacobian); } get_perfs().nb_iterations += 1; m_gradient = m_jacobian.transpose()*m_residuals; SparseSolverReturnCode retcode = solve_linear_system(update); if (retcode != SparseSolverReturnCode::Success) { ERROR << "Error when solving linear system : " << (int) retcode << std::endl; - return_code = ParabolicDriverReturnCode::ErrorLinearSystem; + return ParabolicDriverReturnCode::ErrorLinearSystem; } //if (update.norm() < get_options().step_tolerance) return_code = ParabolicDriverReturnCode::ErrorMinimized; //else return_code = linesearch(update, displacement); } return return_code; } template ParabolicDriverReturnCode ParabolicDriver::check_convergence() { ParabolicDriverReturnCode termcode = ParabolicDriverReturnCode::NotConvergedYet; const int nb_iterations = get_perfs().nb_iterations; const scalar_t norm_residuals = get_perfs().current_residual; const scalar_t norm_update = get_perfs().current_update; //std::cout << "Residuals : " << nb_iterations << " - " << norm_residuals/m_norm_0 << std::endl; DEBUG << "Residuals : " << nb_iterations << " - " << norm_residuals/m_norm_0; if (norm_residuals < get_options().residuals_tolerance) { termcode = ParabolicDriverReturnCode::ResidualMinimized; } else if (nb_iterations > 0 and norm_update > 0 and norm_update < 1.01*get_options().step_tolerance) { if (norm_residuals > get_options().threshold_stationary_point) { ERROR << "Stationary point detected !"; termcode = ParabolicDriverReturnCode::StationaryPoint; } WARNING << "Error is minimized - may indicate a stationnary point"; termcode = ParabolicDriverReturnCode::ErrorMinimized; } else if (nb_iterations > get_options().maximum_iterations) { ERROR << "Maximum number of iteration reached (" << get_options().maximum_iterations << ")"; termcode = ParabolicDriverReturnCode::MaxIterations; } else if (get_perfs().maximum_step_taken) { get_perfs().nb_consecutive_max_step_taken += 1; get_perfs().nb_max_step_taken += 1; if (get_perfs().nb_consecutive_max_step_taken == get_options().max_iterations_at_max_length) { ERROR << "Divergence detected - Maximum step length taken two many times"; termcode = ParabolicDriverReturnCode::MaxStepTakenTooManyTimes; } } else { get_perfs().nb_consecutive_max_step_taken = 0; } get_perfs().return_code = termcode; return termcode; } template SparseSolverReturnCode ParabolicDriver::solve_linear_system(Eigen::VectorXd& update) { return solve_sparse_linear_system(scaling().asDiagonal()*m_residuals, m_jacobian, update, get_options().sparse_solver); } template double ParabolicDriver::compute_residuals_linesearch( Vector &update, scalar_t lambda, Vector &displacement ) { Eigen::VectorXd velocity(m_velocity); Eigen::VectorXd residual = Eigen::VectorXd::Zero(get_neq()); program().update_solution(update, lambda, get_options().alpha*m_current_dt, m_predictor, displacement, velocity); compute_residuals(displacement, velocity, residual); return 0.5*residual.squaredNorm(); } template double ParabolicDriver::compute_residuals_strang_linesearch( Vector &update, scalar_t lambda, Vector &displacement ) { Eigen::VectorXd velocity(m_velocity); Eigen::VectorXd residual = Eigen::VectorXd::Zero(get_neq()); program().update_solution(update, lambda, get_options().alpha*m_current_dt, m_predictor, displacement, velocity); compute_residuals(displacement, velocity, residual); return update.dot(residual); } template void ParabolicDriver::update_variable( const Vector& update, scalar_t lambda, Vector& displacement ) { program().update_solution(update, lambda, get_options().alpha*m_current_dt, m_predictor, displacement, m_velocity); } template ParabolicDriverReturnCode ParabolicDriver::linesearch( Vector &update, Vector &displacements ) { base::is_step_too_long(update); get_perfs().maximum_step_taken = false; scalar_t lambda; ParabolicLinesearchReturnCode retcode; switch (get_options().linesearch) { case ParabolicLinesearch::Bactracking: retcode = backtracking_linesearch(update, displacements, lambda); break; case ParabolicLinesearch::Strang: retcode = strang_linesearch(update, displacements, lambda); break; default: throw std::runtime_error("Linesearch type for Parabolic driver is not recognized"); break; } if (retcode != ParabolicLinesearchReturnCode::Success) { return ParabolicDriverReturnCode::LinesearchFailed; } update_variable(update, lambda, displacements); update *= lambda; return ParabolicDriverReturnCode::NotConvergedYet; } namespace internal { //! \brief Return true if a and b have the same sign inline bool have_same_sign(double a, double b) { return (std::copysign(1., a) * std::copysign(1., b) > 0); } } // end namespace internal template ParabolicLinesearchReturnCode ParabolicDriver::strang_linesearch( Vector& update, Vector& displacements, scalar_t& lambda ) { DEBUG << "Strang linesearch"; Eigen::VectorXd xp(displacements); const scalar_t s_tol = 0.5; const scalar_t s_max = 2.0; const int lin_max = 10; scalar_t s_b = 0.0; scalar_t g_b = 0.5*m_residuals.squaredNorm(); scalar_t s_a = 1.0; scalar_t g_a = compute_residuals_strang_linesearch(update, 1.0, xp); scalar_t newtlen = (scaling().asDiagonal()*update).norm(); // const scalar_t s_r = s_a; const scalar_t g_r = g_a; if (std::abs(g_a) <= s_tol*std::abs(g_b)) { DEBUG << "Skip linesearch "; lambda = 1.0; if (lambda == 1.0 and (newtlen > 0.99 * get_options().maximum_step_length)) { get_perfs().maximum_step_taken = true; } return ParabolicLinesearchReturnCode::Success; } while (internal::have_same_sign(g_a, g_b) and s_a < s_max) { s_b = s_a; s_a = 2*s_a; g_b = g_a; g_a = compute_residuals_strang_linesearch(update, s_a, xp); } scalar_t g_ = g_a; scalar_t g_0 = g_a; scalar_t s = s_a; int l; for (l=0; l= g_r) { WARNING << "Failed to find better update in Strang linesearch"; //lambda = 0.1*s_r; return backtracking_linesearch(update, displacements, lambda); } lambda = s_a; if (lambda == 1.0 and (newtlen > 0.99 * get_options().maximum_step_length)) { get_perfs().maximum_step_taken = true; } return ParabolicLinesearchReturnCode::Success; } template ParabolicLinesearchReturnCode ParabolicDriver::backtracking_linesearch( Vector& update, Vector& displacements, scalar_t& lambda_out ) { // References // ---------- // - Algo A6.3.1 : Dennis and Schnabel (1983) // - Nocedal & Wrigth (2006) DEBUG << "Linesearch"; Eigen::VectorXd xp(displacements); double fcp; int retcode = 2; // 2 not converged, 1 problem, 0 success const scalar_t alpha = 1e-4; scalar_t newtlen = (scaling().asDiagonal()*update).norm(); scalar_t init_slope = m_gradient.dot(update); const scalar_t rellength = update_norm(update); const scalar_t minlambda = get_options().step_tolerance / rellength; scalar_t lambda = 1.0; scalar_t lambda_prev = lambda; scalar_t merit_value = 0.5*m_residuals.squaredNorm(); // new residual fcp = compute_residuals_linesearch(update, lambda, xp); // Skip linesearch if enough progress is done // ------------------------------------------ if (fcp < get_options().coeff_accept_newton_step *merit_value) { DEBUG << "Skip linesearch "; lambda_out = 1.0; return ParabolicLinesearchReturnCode::Success; } // The linesearch // -------------- scalar_t fc = merit_value; scalar_t fcp_prev; int cnt = 0; do { SPAM << "cnt : " < 0.99 * get_options().maximum_step_length)) { get_perfs().maximum_step_taken = true; } break; } else if (lambda < minlambda) { retcode = 0; lambda = minlambda; //retcode = 1; break; } else { //WARNING << "denom : " << fcp - fc -init_slope; // Select a new step length // - - - - - - - - - - - - double lambdatmp; if (cnt == 0) { // only a quadratic at the first lambdatmp = - init_slope / (2*(fcp - fc -init_slope)); } else { const scalar_t factor = 1.0 /(lambda - lambda_prev); const scalar_t x1 = fcp - fc - lambda*init_slope; const scalar_t x2 = fcp_prev - fc - lambda_prev*init_slope; const scalar_t a = factor * ( x1/(lambda*lambda) - x2/(lambda_prev*lambda_prev)); const scalar_t b = factor * ( -x1*lambda_prev/(lambda*lambda) + x2*lambda/(lambda_prev*lambda_prev)); if (a == 0) { // cubic interpolation is in fact a quadratic interpolation lambdatmp = - init_slope/(2*b); } else { const scalar_t disc = b*b-3*a*init_slope; lambdatmp = (-b+std::sqrt(disc))/(3*a); } if (lambdatmp > 0.5*lambda ) lambdatmp = 0.5*lambda; } //WARNING << "lambdatmp : " << lambdatmp; lambda_prev = lambda; fcp_prev = fcp; if (lambdatmp < 0.1*lambda) { lambda = 0.1 * lambda; } else { lambda = lambdatmp; } } if (not std::isfinite(lambda)) { ERROR << "Lambda is non finite - we stop"; return ParabolicLinesearchReturnCode::Divergence; } fcp = compute_residuals_linesearch(update, lambda, xp); //xp.velocity = x.velocity; ++cnt; } while(retcode == 2 and cnt < 50); //WARNING << "Lambda : " << lambda << " - iterations : " << cnt; if (cnt == 50) { ERROR << "Too much linesearch iterations ! We stop"; return ParabolicLinesearchReturnCode::MaximumIterations; } lambda_out = lambda; switch (retcode) { case 0: return ParabolicLinesearchReturnCode::Success; break; case 1: return ParabolicLinesearchReturnCode::LambdaTooSmall; default: return ParabolicLinesearchReturnCode::NotSupposedToHappen; break; } } template void ParabolicDriver::set_velocity(Vector& velocity_vector) { m_velocity.resize(program().get_tot_ndf()); m_velocity.setZero(); for (index_t dof=0; dof void ParabolicDriver::reset_velocity() { for (index_t dof=0; dof #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 ): OptionsHandler(), PerformanceHandler(), m_mesh(the_mesh), m_database(the_database), m_variables(the_database, the_mesh, transport_parameters), m_transport_parameters(transport_parameters), m_transport_program(the_mesh, the_database, transport_parameters), m_transport_solver(m_transport_program), m_bcs(the_mesh->nb_nodes()), m_permanent_flag_compute_speciation(Eigen::VectorXi::Ones(the_mesh->nb_nodes())) {} void SIASaturatedReactiveTransportSolver::apply_boundary_conditions(SIABoundaryConditions bcs) { 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); } set_scaling(); } void SIASaturatedReactiveTransportSolver::set_scaling() { Vector scaling(m_transport_program.get_neq()); for (index_t node: m_mesh->range_nodes()) { scalar_t typical_v = m_transport_parameters->diffusion_coefficient(node); index_t element; if (node != m_mesh->nb_nodes()-1) element = node; else element = node-1; typical_v *= m_mesh->get_face_area(element)/(m_mesh->get_dx(element)*m_mesh->get_volume_cell(node)); for (index_t component: m_database->range_aqueous_component()) { const index_t id_eq = m_transport_program.id_equation(m_transport_program.get_dof_component(node, component)); if (id_eq == no_equation) continue; scaling(id_eq) = 1.0/typical_v; } } m_transport_solver.set_scaling(scaling); } scalar_t SIASaturatedReactiveTransportSolver::residuals_transport(SIASaturatedVariables& variables) { Eigen::VectorXd residuals; m_transport_solver.compute_residuals(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) { // Initialization // ############## // Secondary variables necessary to the algorithm // ---------------------------------------------- retcode_t return_code = retcode_t::NotConvergedYet; internals::SIAResiduals residuals_storage(m_mesh->nb_nodes()); reset_perfs(); m_flag_compute_speciation = m_permanent_flag_compute_speciation; // Copy main variables // -------------------- // We copy main variables to recover from a failure SIASaturatedVariables tmp_variables(m_variables); // initialization transport program and transport solver Eigen::VectorXd transport_variable(tmp_variables.total_mobile_concentrations()); m_transport_solver.get_options() = get_options().transport_options; if (not get_options().use_previous_flux) m_transport_program.external_flow().setZero(); // Solid phases, used to compute ΔB_i/Δt m_minerals_start_iteration = tmp_variables.speciation_variables().block( m_database->nb_component, 0, m_database->nb_mineral, m_mesh->nb_nodes()); // First Iteration // ############### // Compute initial residual // ------------------------ - residuals_storage.transport_residual_0 = residuals_transport(tmp_variables); + m_transport_solver.initialize_timestep(dt, transport_variable); + residuals_storage.transport_residual_0 = m_transport_solver.residuals_0(); // Solve transport operator // ------------------------ - dfpmsolver::ParabolicDriverReturnCode retcode = m_transport_solver.solve_timestep(dt, transport_variable); + dfpmsolver::ParabolicDriverReturnCode retcode = m_transport_solver.restart_timestep(transport_variable); // If we cannot solve the linear system, we try SparseQR solver if (retcode == dfpmsolver::ParabolicDriverReturnCode::ErrorLinearSystem) { get_perfs().error_transport_linearsystem += 1; - if (m_transport_solver.get_options().sparse_solver != SparseSolver::SparseQR) + if (get_options().sparse_qr_fallback + and m_transport_solver.get_options().sparse_solver != SparseSolver::SparseQR) { m_transport_solver.get_options().sparse_solver = SparseSolver::SparseQR; retcode = m_transport_solver.restart_timestep(transport_variable); } } if (retcode != dfpmsolver::ParabolicDriverReturnCode::ResidualMinimized) { ERROR << "Failed to solve 1st transport problem; return code " << (int) retcode; return convert_dfpm_retcode(retcode); } // save informations about the transport solver residuals_storage.transport_residual = m_transport_solver.get_perfs().current_residual; residuals_storage.transport_update = m_transport_solver.get_perfs().current_update; get_perfs().nb_iterations_transport = m_transport_solver.get_perfs().nb_iterations; // Update // ------ tmp_variables.total_mobile_concentrations() = transport_variable; // Speciation // ---------- restart_speciation_iterations(dt, residuals_storage, tmp_variables); // speciation solver update directly tmp_variables (node by node) // If we use the consitent velocity if (get_options().use_consistent_velocity) m_transport_solver.velocity() = ( tmp_variables.total_mobile_concentrations() - m_variables.total_mobile_concentrations() )/dt; - + //std::cout << "norm velocity 0 : " << m_transport_solver.get_velocity().norm() << std::endl; + //std::cout << "norm flow 0 : " << m_transport_program.external_flow().norm() << std::endl; // Check convergence // ----------------- return_code = check_convergence(tmp_variables, residuals_storage); residuals_storage.nb_iterations +=1; get_perfs().nb_iterations += 1; // If SNIA is used, we stop // ------------------------ if (get_options().snia()) { // Copy the tmp variables into the true variables m_variables = tmp_variables; // Return success return retcode_t::Success; } // SIA algorithm // ============= // 2nd iterations and more while (return_code == retcode_t::NotConvergedYet) { // update variables for the transport transport_variable = tmp_variables.total_mobile_concentrations(); retcode = m_transport_solver.restart_timestep(transport_variable); if ((int) retcode <= 0) // if (error) { // If we cannot solve the transport linear system, we try the SparseQR solver if (retcode == dfpmsolver::ParabolicDriverReturnCode::ErrorLinearSystem) { get_perfs().error_transport_linearsystem += 1; - if (m_transport_solver.get_options().sparse_solver != SparseSolver::SparseQR) + if (get_options().sparse_qr_fallback + and m_transport_solver.get_options().sparse_solver != SparseSolver::SparseQR) { m_transport_solver.get_options().sparse_solver = SparseSolver::SparseQR; retcode = m_transport_solver.restart_timestep(transport_variable); } } // If maximum number of iterations is reached, check if the solution is good enough 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) { // compute equilibrium one last time restart_speciation_iterations(dt, residuals_storage, tmp_variables); m_variables = tmp_variables; return_code = retcode_t::Success; break; } else { ERROR << "Failed to solve transport problem; return code " << (int) retcode; return convert_dfpm_retcode(retcode); } } // save informations about the transport solver residuals_storage.transport_residual = m_transport_solver.get_perfs().current_residual; residuals_storage.transport_update = m_transport_solver.get_perfs().current_update; get_perfs().nb_iterations_transport = m_transport_solver.get_perfs().nb_iterations; // Update // ------ tmp_variables.total_mobile_concentrations() = transport_variable; // Speciation // ---------- restart_speciation_iterations(dt, residuals_storage, tmp_variables); // This function also do the update, node by node // If we use the consitent velocity, update velocity if (get_options().use_consistent_velocity) m_transport_solver.velocity() = ( tmp_variables.total_mobile_concentrations() - m_variables.total_mobile_concentrations() )/dt; + //std::cout << "norm velocity : " << m_transport_solver.get_velocity().norm() << std::endl; + //std::cout << "norm flow : " << m_transport_program.external_flow().norm() << std::endl; + // Iteration is finished residuals_storage.nb_iterations +=1; get_perfs().nb_iterations += 1; // Check convergence // ----------------- return_code = check_convergence(tmp_variables, residuals_storage); } // copy the solution if it is good // ------------------------------- std::cout << (int) return_code << " ?== " << (int) SIASaturatedReactiveTransportSolverReturnCode::Success << std::endl; if (return_code == SIASaturatedReactiveTransportSolverReturnCode::Success) { std::cout << "replace variables" << std::endl; m_variables = tmp_variables; // compute new update threshold // ---------------------------- scalar_t update = m_transport_solver.velocity().maxCoeff(); get_options().transport_options.step_tolerance = get_options().relative_update*update; } return return_code; } void SIASaturatedReactiveTransportSolver::restart_transport_iterations( internals::SIAResiduals& residuals_storage ) { Eigen::VectorXd transport_variable(m_variables.total_mobile_concentrations()); dfpmsolver::ParabolicDriverReturnCode retcode = m_transport_solver.restart_timestep(transport_variable); if (retcode != dfpmsolver::ParabolicDriverReturnCode::ResidualMinimized) { ERROR << "Failed to solve transport problem; return code " << (int) retcode; throw std::runtime_error("Failed to solve transport program"); } residuals_storage.transport_residual = m_transport_solver.get_perfs().current_residual; residuals_storage.transport_update = m_transport_solver.get_perfs().current_update; } void SIASaturatedReactiveTransportSolver::restart_speciation_iterations( scalar_t dt, internals::SIAResiduals& residuals_storage, SIASaturatedVariables& variables ) { #ifdef _OPENMP #pragma omp parallel { #pragma omp for schedule(runtime) for (index_t node=0; nodenb_nodes(); ++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, variables); } } #ifdef _OPENMP } #endif } void SIASaturatedReactiveTransportSolver::solve_speciation_node( scalar_t dt, index_t node, internals::SIAResiduals& residuals_storage, SIASaturatedVariables& variables ) { Eigen::VectorXd totconc; variables.nodal_update_total_amount(node, totconc); ReducedSystemSolver solver(m_database, totconc, get_options().speciation_options); Eigen::VectorXd speciation_variables = variables.speciation_variables(node); micpsolver::MiCPPerformance perf = solver.solve(speciation_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(speciation_variables); variables.update_composition(node, solution); // mineral const scalar_t denom = 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)*( variables.mineral_amount(node, mineral) -m_minerals_start_iteration(mineral, node) ); } m_transport_program.external_flow(node, component) = sum/denom; } if (solution.moles_minerals() == 0) m_permanent_flag_compute_speciation(node) = 0; #ifdef _OPENMP #pragma omp critical get_perfs().nb_iterations_chemistry += perf.nb_iterations; #else get_perfs().nb_iterations_chemistry += perf.nb_iterations; #endif } SIASaturatedReactiveTransportSolverReturnCode SIASaturatedReactiveTransportSolver::check_convergence( SIASaturatedVariables& variables, internals::SIAResiduals& residuals_storage ) { // compute residuals of transport const scalar_t norm_transport = residuals_transport(variables); get_perfs().current_residual = norm_transport/residuals_storage.transport_residual_0; + DEBUG << "RESIDUAL : " << get_perfs().current_residual << " - " << norm_transport; retcode_t return_code = retcode_t::NotConvergedYet; if (norm_transport/residuals_storage.transport_residual_0 < get_options().residual_tolerance) { + std::cout << "Residual minimized" << std::endl; return_code = retcode_t::Success; } else if (residuals_storage.transport_update < get_options().threshold_update_transport) { + std::cout << "Error minimized" << std::endl; 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"; if (norm_transport/residuals_storage.transport_residual_0 < get_options().good_enough_transport_residual_tolerance ) return_code = retcode_t::Success; return_code = retcode_t::MaxIterations; } return return_code; } } // end namespace reactmicp } // end namespace specmicp diff --git a/src/reactmicp/systems/saturated_diffusion/reactive_transport_solver_structs.hpp b/src/reactmicp/systems/saturated_diffusion/reactive_transport_solver_structs.hpp index 698da94..ade9f8d 100644 --- a/src/reactmicp/systems/saturated_diffusion/reactive_transport_solver_structs.hpp +++ b/src/reactmicp/systems/saturated_diffusion/reactive_transport_solver_structs.hpp @@ -1,97 +1,99 @@ #ifndef SPECMICP_REACTMICP_SATURATED_DIFFUSION_OPTIONS_HPP #define SPECMICP_REACTMICP_SATURATED_DIFFUSION_OPTIONS_HPP #include "common.hpp" #include "dfpmsolver/parabolic_structs.hpp" #include "micpsolver/micpsolver_structs.hpp" #include "specmicp/reduced_system_solver.hpp" namespace specmicp { namespace reactmicp { namespace systems { namespace siasaturated { //! \brief Return codes for the SIASaturatedReactiveTransportSolver enum class SIASaturatedReactiveTransportSolverReturnCode { LolThatsNotSuposedToHappen = -30, //!< Well, have fun debugging that... // Transport error codes ErrorInTransport = -20, //!< Error in the transport solver MaxIterationsInTransport = -21, //!< Maximum iterations reached in the transport solver StationaryPointsInTransport = -22, //!< Stationnary points in the transport solver // Speciation error codes ErrorInSpeciation = -10, //!< Error in the speciation solver MaxIterationsInSpeciation = -11, //!< Maximum iterations reached in the speciation solver StationaryPointsInSpeciation = -12, //!< Stationnary points in the speciation solver // Coupling algorithm MaxIterations = -1, //!< Maximum iterations of the coupling solver StationaryPoints = -2, //!< Stationary points reached in the coupling solver // The good ones NotConvergedYet = 0, //!< Problem is not converged yet Success = 1 //!< The problem is succesfully solved }; //! \brief Options for the SIASaturatedReactiveTransportSolver struct SIASaturatedReactiveTransportSolverOptions { int max_iterations; //!< Maximum number of fixed-point iterations scalar_t residual_tolerance; //!< The tolerance for fixed-point convergence scalar_t good_enough_transport_residual_tolerance; //!< If transport failed, this value is used to check if it is 'good' enough scalar_t threshold_update_transport; //!< Threshold to accept a step if the update is too small scalar_t threshold_update_disable_speciation; //!< Threshold to disable next computation of speciation at a node scalar_t relative_update; //!< To compute a new update threshold bool use_consistent_velocity; //!< True to compute residuals using speciation velocity bool use_previous_flux; //!< Use flux of previous timestep to accelerate the computation when possible bool speciation_for_solid_only; //!< Compute speciation only when there are some solids left + bool sparse_qr_fallback; //!< use sparse qr if linear system could not be solved dfpmsolver::ParabolicDriverOptions transport_options; //!< Options to use for the transport solver specmicp::ReducedSystemSolverOptions speciation_options; //!< Options for the speciation solver SIASaturatedReactiveTransportSolverOptions(): max_iterations(50), residual_tolerance(1e-3), good_enough_transport_residual_tolerance(1e-2), threshold_update_transport(1e-10), threshold_update_disable_speciation(1e-10), use_consistent_velocity(false), use_previous_flux(false), - speciation_for_solid_only(false) + speciation_for_solid_only(false), + sparse_qr_fallback(true) {} //! \brief Use a sequential non-iterative algorithm void use_snia() {max_iterations=1;} //! \brief Return true if snia is used bool snia() {return max_iterations == 1;} //! \brief Use a sequential iterative algorithm void use_sia(int nb_iter=50) {max_iterations=nb_iter;} //! \brief Set the tolerance void set_tolerance(scalar_t tol) {residual_tolerance = tol;} }; struct SIASaturatedReactiveTransportSolverPerformance { index_t nb_iterations; index_t nb_iterations_transport; index_t nb_iterations_chemistry; scalar_t current_residual; scalar_t previous_update; index_t error_transport_linearsystem; SIASaturatedReactiveTransportSolverPerformance(): nb_iterations(0), nb_iterations_transport(0), nb_iterations_chemistry(0), current_residual(-1.0), previous_update(0), error_transport_linearsystem(0) {} }; } // end namespace siasaturated } // end namespace systems } // end namespace reactmicp } // end namespace specmicp #endif // SPECMICP_REACTMICP_SATURATED_DIFFUSION_OPTIONS_HPP diff --git a/src/reactmicp/systems/saturated_diffusion/transport_program.cpp b/src/reactmicp/systems/saturated_diffusion/transport_program.cpp index 87ab1cd..cd62430 100644 --- a/src/reactmicp/systems/saturated_diffusion/transport_program.cpp +++ b/src/reactmicp/systems/saturated_diffusion/transport_program.cpp @@ -1,236 +1,236 @@ #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 { index_t SaturatedDiffusionProgram::get_tot_ndf() const { return m_mesh->nb_nodes()*get_ndf(); } 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_flow(Vector::Zero(the_mesh->nb_nodes()*(the_data->nb_component-1))), m_internal_flow(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 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, 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.0); const index_t node_0 = m_mesh->get_node(element, 0); const index_t node_1 = m_mesh->get_node(element, 1); const scalar_t porosity = (m_param->porosity(node_0) + m_param->porosity(node_1) )/2.0; const scalar_t diff_coeff = 1.0/(0.5/m_param->diffusion_coefficient(node_0) + 0.5/m_param->diffusion_coefficient(node_1)); const index_t dof_0 = get_dof_component(node_0, component); const index_t dof_1 = get_dof_component(node_1, component); scalar_t flux_coeff = -( m_mesh->get_face_area(element) / m_mesh->get_dx(element) * porosity * m_param->density_water() * diff_coeff ); jacob << 1.0, -1.0, -1.0, 1.0; jacob *= flux_coeff; evelocity << velocity(dof_0)* mass_coeff*m_param->porosity(node_0), velocity(dof_1)* mass_coeff*m_param->porosity(node_1); edisplacement << displacement(dof_0), displacement(dof_1); element_residual = jacob*edisplacement; m_internal_flow(dof_0) += element_residual(0); m_internal_flow(dof_1) += element_residual(1); for (index_t en: m_mesh->range_nen()) { element_residual(en) += evelocity(en); element_residual(en) += 1e6*(m_mesh->get_volume_element(element)/2.0 *external_flow(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; + if (h < 1e-4*EPS_J) 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