Diffusion SpecMiCP / ReactMiCP (cython) Browse

reactive_transport_neutrality_solver.cpp
6d08d8744732cython

reactive_transport_neutrality_solver.cpp
View Options

			//#define EIGEN_DEFAULT_IO_FORMAT IOFormat(8)

			#include "reactive_transport_neutrality_solver.hpp"
			#include "specmicp/reduced_system_solver.hpp"
			#include "boundary_conditions.hpp"
			#include "reactmicp/meshes/mesh1d.hpp"
			#include "reactmicp/systems/saturated_diffusion/transport_neutrality_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


			SIASaturatedReactiveTransportNeutralitySolver::SIASaturatedReactiveTransportNeutralitySolver(
			std::shared_ptr<mesh::Mesh1D> the_mesh,
			RawDatabasePtr the_database,
			std::shared_ptr<SaturatedNeutralityDiffusionTransportParameters> transport_parameters
			):
			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())
			{}

			void SIASaturatedReactiveTransportNeutralitySolver::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.extern_update_nodal_composition(node,
			bcs.list_initial_states[bcs.initial_states[node]]);
			m_transport_program.number_equations(list_fixed_nodes);
			}
			set_scaling();
			}

			void SIASaturatedReactiveTransportNeutralitySolver::set_scaling()
			{
			Vector scaling(m_transport_program.get_neq());
			for (index_t node: m_mesh->range_nodes())
			{

			scalar_t typical_v = 1.0;
			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) = 1000.0/(typical_v*m_transport_parameters->diffusion_component(node, component));
			}
			}
			m_transport_solver.set_scaling(scaling);
			}

			scalar_t SIASaturatedReactiveTransportNeutralitySolver::residuals_transport(const Vector& transport_variables)
			{
			Eigen::VectorXd residuals;
			m_transport_solver.compute_residuals(transport_variables, residuals);
			return residuals.norm();
			}

			namespace internal {

			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;
			}
			}
			} // end namespace internal

			SIASaturatedReactiveTransportSolverReturnCode SIASaturatedReactiveTransportNeutralitySolver::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_variables(m_transport_program.get_tot_ndf());
			for (index_t node: m_mesh->range_nodes())
			{
			for (index_t component: m_database->range_aqueous_component())
			{
			transport_variables(m_transport_program.get_dof_component(node, component)) =
			m_variables.log_component_concentration(node, component);
			}
			}
			m_transport_solver.get_options() = get_options().transport_options;
			m_transport_program.external_flux().setZero();
			//
			residuals_storage.transport_residual_0 = residuals_transport(transport_variables);
			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_variables);
			//if (retcode != dfpmsolver::ParabolicDriverReturnCode::ResidualMinimized)
			if ((int) retcode <= 0)
			{
			ERROR << "Failed to solve 1st transport problem; return code " << (int) retcode;
			return internal::convert_dfpm_retcode(retcode);
			}
			residuals_storage.transport_residual = m_transport_solver.get_perfs().current_residual;
			residuals_storage.transport_update = 1.; // don't check update the first time
			//residuals_storage.transport_update = m_transport_solver.get_perfs().current_update;
			// update
			Vector total_concentrations(m_transport_program.get_tot_ndf());
			m_transport_program.get_total_mobile_concentrations(transport_variables, total_concentrations);
			m_variables.total_mobile_concentrations() = total_concentrations;

			for (index_t node: m_mesh->range_nodes())
			{
			scalar_t charge = m_variables.charge_balance(node);
			std::cout << " " << charge;
			if (std::abs(charge) > 1e-6)
			m_variables.log_component_concentration(node, 2) = std::log10(m_variables.component_concentration(node, 2) + charge);
			}
			std::cout << std::endl;
			// Speciation
			// ----------
			restart_speciation_iterations(dt, residuals_storage, transport_variables);
			// Convergence
			// -----------
			return_code = check_convergence(residuals_storage, transport_variables);
			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;
			retcode = m_transport_solver.restart_timestep(transport_variables);
			//if (retcode != dfpmsolver::ParabolicDriverReturnCode::ResidualMinimized)
			if ((int) retcode <= 0)
			{
			if (retcode == dfpmsolver::ParabolicDriverReturnCode::MaxIterations and
			m_transport_solver.get_perfs().current_residual
			< get_options().good_enough_transport_residual_tolerance)
			{
			restart_speciation_iterations(dt, residuals_storage, transport_variables);
			return_code = retcode_t::Success;
			break;
			}

			ERROR << "Failed to solve transport problem; return code " << (int) retcode;
			return internal::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;

			Vector total_concentrations(m_transport_program.get_tot_ndf());
			m_transport_program.get_total_mobile_concentrations(transport_variables, total_concentrations);
			m_variables.total_mobile_concentrations() = total_concentrations;

			residuals_storage.nb_iterations +=1;
			restart_speciation_iterations(dt, residuals_storage, transport_variables);

			return_code = check_convergence(residuals_storage, transport_variables);
			}
			std::cout << "End fixed point : return code " << (int) return_code
			<< " - nb iter : " <<residuals_storage.nb_iterations
			<< std::endl;

			set_scaling();

			return return_code;
			}

			void SIASaturatedReactiveTransportNeutralitySolver::restart_transport_iterations(
			internals::SIAResiduals& residuals_storage
			)
			{
			Eigen::VectorXd transport_variable(m_transport_program.get_tot_ndf());
			for (index_t node: m_mesh->range_nodes())
			{
			for (index_t component: m_database->range_aqueous_component())
			{
			transport_variable(m_transport_program.get_dof_component(node, component)) =
			m_variables.log_component_concentration(node, component);
			}
			}
			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;

			Vector total_concentrations(m_transport_program.get_tot_ndf());
			m_transport_program.get_total_mobile_concentrations(transport_variable, total_concentrations);
			m_variables.total_mobile_concentrations() = total_concentrations;
			}

			void SIASaturatedReactiveTransportNeutralitySolver::restart_speciation_iterations(
			scalar_t dt,
			internals::SIAResiduals& residuals_storage,
			Vector& transport_variables
			)
			{
			#ifdef _OPENMP
			#pragma omp parallel
			{
			#pragma omp for schedule(static)
			for (index_t node=0; node<m_mesh->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, transport_variables);
			}
			}
			#ifdef _OPENMP
			}
			#endif
			}


			void SIASaturatedReactiveTransportNeutralitySolver::solve_speciation_node(
			scalar_t dt,
			index_t node,
			internals::SIAResiduals& residuals_storage,
			Vector& transport_variables
			)
			{
			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

			m_transport_program.extern_update_nodal_composition(node, solution, transport_variables);


			const scalar_t denom = 1e-3m_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 SIASaturatedReactiveTransportNeutralitySolver::check_convergence(
			internals::SIAResiduals& residuals_storage,
			const Vector& transport_variables
			)
			{
			// compute residuals of transport
			const scalar_t norm_transport = residuals_transport(transport_variables);
			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

specmicp/src/reactmicp/systems/saturated_diffusion/reactive_transport_neutrality_solver.cpp6d08d8744732cython

reactive_transport_neutrality_solver.cppView Options

specmicp/src/reactmicp/systems/saturated_diffusion/reactive_transport_neutrality_solver.cpp
6d08d8744732cython

reactive_transport_neutrality_solver.cpp
View Options