Diffusion SpecMiCP / ReactMiCP (master) Browse

reactive_transport_solver.cpp
511a07d7f37cmaster

reactive_transport_solver.cpp
View Options

			//#define EIGEN_DEFAULT_IO_FORMAT IOFormat(8)

			#include "reactive_transport_solver.hpp"
			#include "specmicp/reduced_system_solver.hpp"
			#include "boundary_conditions.hpp"
			#include "dfpm/meshes/mesh1d.hpp"
			#include "reactmicp/systems/saturated_diffusion/transport_parameters.hpp"

			#ifdef SPECMICP_USE_OPENMP
			#include <omp.h>
			#endif // SPECMICP_USE_OPENMP

			#include <iostream>

			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<mesh::Mesh1D> the_mesh,
			RawDatabasePtr the_database,
			std::shared_ptr<SaturatedDiffusionTransportParameters> transport_parameters
			):
			OptionsHandler<SIASaturatedReactiveTransportSolverOptions>(),
			PerformanceHandler<SIASaturatedReactiveTransportSolverPerformance>(),
			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<index_t> 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;
			}
			}

			void SIASaturatedReactiveTransportSolver::reset_flow()
			{
			m_transport_program.external_flow().setZero();
			}

			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
			// ------------------------
			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.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 (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 (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 << "Final residual " << residuals_storage.transport_residual << " - "
			<< residuals_storage.transport_residual/residuals_storage.transport_residual_0 << std::endl;
			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 = std::abs(m_transport_solver.velocity().maxCoeff());
			get_options().transport_options.step_tolerance = get_options().relative_update*update;

			if (get_options().disable_speciation_no_solid)
			{
			for (index_t node: m_mesh->range_nodes())
			{
			if (m_variables.moles_minerals(node) == 0.0) m_permanent_flag_compute_speciation(node) = 0;
			}
			}

			}
			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; node<m_mesh->nb_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 ((int) perf.return_code <= 0) //!= 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;
			}

			#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
			or norm_transport < get_options().absolute_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

specmicp/src/reactmicp/systems/saturated_diffusion/reactive_transport_solver.cpp511a07d7f37cmaster

reactive_transport_solver.cppView Options

specmicp/src/reactmicp/systems/saturated_diffusion/reactive_transport_solver.cpp
511a07d7f37cmaster

reactive_transport_solver.cpp
View Options