reactive_transport_solver.cpp
No OneTemporary
Actions

Subscribers

None

File Metadata

Created: Mon, May 20, 22:20

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

reactive_transport_solver.cppNo OneTemporaryActions

File Metadata

reactive_transport_solver.cppView Options

Event Timeline

reactive_transport_solver.cpp
No OneTemporary
Actions

reactive_transport_solver.cpp
View Options