Diffusion SpecMiCP / ReactMiCP (reactmicp/unsaturated) Browse

unsaturated

transport_stagger.cpp
View Options

			/*-------------------------------------------------------------------------------

			Copyright (c) 2015 F. Georget <fabieng@princeton.edu>, 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:

			1. Redistributions of source code must retain the above copyright notice, this
			list of conditions and the following disclaimer.

			2. 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.

			3. Neither the name of the copyright holder 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 HOLDER 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 "transport_stagger.hpp"
			#include "../../solver/staggers_base/stagger_structs.hpp"
			#include "variables.hpp"
			#include "variables_box.hpp"

			#include "saturation_equation.hpp"
			#include "pressure_equation.hpp"
			#include "aqueous_equation.hpp"

			#include "../../../database/database_holder.hpp"

			#include "../../../dfpmsolver/parabolic_driver.hpp"

			#include "../../../utils/log.hpp"

			namespace specmicp {
			namespace reactmicp {
			namespace systems {
			namespace unsaturated {

			//! \brief Type of the equation
			//! \internal
			enum class EquationType {
			Saturation,
			LiquidAqueous,
			Pressure
			};

			//! \brief Format type to a string (for message error)
			//! \internal
			static std::string to_string(EquationType eq_type);

			//! \brief Format DFPM solver retcode to a string
			static std::string to_string(dfpmsolver::ParabolicDriverReturnCode retcode);

			// forward declaration of wrapper classes
			// They are defined at the bootom of this file
			class TaskBase;
			class SaturationTask;

			using VariablesBasePtr = solver::VariablesBasePtr;


			// =================================== //
			// //
			// Declaration of the implementation //
			// //
			// =================================== //


			//! \brief Implementation class for the Unsaturated transport stagger
			//! \internal
			//!
			//! This class does all the work
			//! It was designed to be OpenMP compatible
			class UnsaturatedTransportStagger::UnsaturatedTransportStaggerImpl
			{
			public:

			// Constructor
			// -----------

			UnsaturatedTransportStaggerImpl(
			UnsaturatedVariablesPtr variables,
			const TransportConstraints& constraints
			);

			// Main functions
			// --------------

			//! \brief Return the residual
			scalar_t get_residual(UnsaturatedVariables* vars);
			//! \brief Return the first residual (start of timestep)
			scalar_t get_residual_0(UnsaturatedVariables* vars);
			//! \brief Return the update (norm of velocity vector)
			scalar_t get_update(UnsaturatedVariables* vars);

			//! \brief Initialize timestep
			void initialize_timestep(
			scalar_t dt,
			UnsaturatedVariables* vars);
			//! \brief Restart the timestep
			solver::StaggerReturnCode restart_timestep(UnsaturatedVariables* vars);

			// Timestep management
			// -------------------

			//! \brief Save the timestep
			void save_dt(scalar_t dt) {m_dt = dt;}
			//! \brief Return the timestep
			scalar_t get_dt() {return m_dt;}

			private:
			scalar_t m_norm_0 {1.0};
			// timestep management
			scalar_t m_dt {-1.0};

			//! \brief The saturation equations
			std::unique_ptr<SaturationTask> m_saturation_equation;
			// Equation and solver
			std::vector<std::unique_ptr<TaskBase>> m_equation_list;


			};

			// Main functions
			// ===============

			// call of the implementation


			UnsaturatedTransportStagger::UnsaturatedTransportStagger(
			std::shared_ptr<UnsaturatedVariables> variables,
			const TransportConstraints& constraints
			):
			m_impl(make_unique<UnsaturatedTransportStaggerImpl>(variables, constraints))
			{
			}

			UnsaturatedTransportStagger::~UnsaturatedTransportStagger() = default;


			static UnsaturatedVariables* get_var(VariablesBasePtr var) {
			return static_cast<UnsaturatedVariables*>(var.get());
			}

			void
			UnsaturatedTransportStagger::initialize_timestep(
			scalar_t dt,
			VariablesBasePtr var
			)
			{
			UnsaturatedVariables* true_var = get_var(var);
			m_impl->initialize_timestep(dt, true_var);

			}

			solver::StaggerReturnCode
			UnsaturatedTransportStagger::restart_timestep(VariablesBasePtr var)
			{
			UnsaturatedVariables* true_var = get_var(var);
			return m_impl->restart_timestep(true_var);
			}

			scalar_t
			UnsaturatedTransportStagger::get_residual(VariablesBasePtr var)
			{
			UnsaturatedVariables* true_var = get_var(var);
			return m_impl->get_residual(true_var);
			}

			scalar_t
			UnsaturatedTransportStagger::get_residual_0(VariablesBasePtr var)
			{
			UnsaturatedVariables* true_var = get_var(var);
			return m_impl->get_residual_0(true_var);
			}

			scalar_t
			UnsaturatedTransportStagger::get_update(VariablesBasePtr var)
			{
			UnsaturatedVariables* true_var = get_var(var);
			return m_impl->get_update(true_var);
			}

			// =============================== //
			// =============================== //
			// //
			// Implementation details //
			// ---------------------- //
			// //
			// =============================== //
			// =============================== //

			// 2 main sections
			// - Solver wrappers : wrapper around 1 couple equation/solver
			// - UnsaturatedTransportStaggerImpl : call the wrappers

			// =============================== //
			// //
			// Solver wrappers //
			// //
			// =============================== //

			// This wrappers are used to abstract the residual computation
			// and the methods to solve every equations


			//! \brief Base class for a task
			//!
			//! A task is how we solve governing equations,
			//! and obtain residuals and update
			//!
			//! \internal
			class SPECMICP_DLL_LOCAL TaskBase
			{
			public:
			TaskBase(index_t component, EquationType eq_type):
			m_type(eq_type),
			m_component(component)
			{}


			//! \brief Compute the squared residuals
			virtual scalar_t compute_squared_residual(UnsaturatedVariables* vars) = 0;
			//! \brief Compute the squared residuals at the beginning of the timestep
			virtual scalar_t compute_squared_residual_0(UnsaturatedVariables* vars) = 0;
			//! \brief Compute the squared update of the variables
			virtual scalar_t compute_squared_update(UnsaturatedVariables* vars) = 0;

			//! \brief Initialize the timestep
			virtual void initialize_timestep(scalar_t dt, UnsaturatedVariables* vars) = 0;
			//! \brief Solve the governing equation
			virtual dfpmsolver::ParabolicDriverReturnCode restart_timestep(UnsaturatedVariables* vars) = 0;

			//! \brief Return the component index (in the db)
			index_t component() {return m_component;}

			//! \brief The equation type
			EquationType equation_type() {return m_type;}

			private:
			EquationType m_type;
			index_t m_component;
			};


			//! \brief Base class for a equation solver
			//!
			//! \internal
			template <typename Derived, typename traits>
			class SPECMICP_DLL_LOCAL EquationTask: public TaskBase
			{
			public:
			using EqT = typename traits::EqT;
			using SolverT = typename traits::SolverT;

			EquationTask(
			index_t component,
			mesh::Mesh1DPtr the_mesh,
			typename traits::VariableBoxT& variables,
			const TransportConstraints& constraints
			):
			TaskBase(component, traits::equation_type),
			m_equation(the_mesh, variables, constraints),
			m_solver(m_equation)
			{}

			EquationTask(
			index_t component,
			mesh::Mesh1DPtr the_mesh,
			typename traits::VariableBoxT& variables,
			const TransportConstraints& constraints,
			scalar_t scaling
			):
			TaskBase(component, traits::equation_type),
			m_equation(the_mesh, variables, constraints),
			m_solver(m_equation)
			{
			m_equation.set_scaling(scaling);
			}

			Derived* derived() {return static_cast<Derived*>(this);}

			// This function must be implemented by subclass
			MainVariable& get_var(UnsaturatedVariables* vars) {
			return derived()->get_var(vars);
			}

			scalar_t compute_squared_residual(UnsaturatedVariables* vars) override {
			Vector residuals;
			const MainVariable& main = get_var(vars);
			m_equation.compute_residuals(main.variable, main.velocity, residuals, true);
			return residuals.squaredNorm()/m_norm_square_0;
			}
			scalar_t compute_squared_residual_0(UnsaturatedVariables* vars) override {
			Vector residuals;
			const MainVariable& main = get_var(vars);
			m_equation.compute_residuals(main.variable, main.velocity, residuals, false);
			m_norm_square_0 = residuals.squaredNorm();
			if (m_norm_square_0 < 1e-8) m_norm_square_0 = 1.0;
			return 1.0;
			}
			scalar_t compute_squared_update(UnsaturatedVariables *vars) override {
			const MainVariable& main = get_var(vars);
			return main.velocity.squaredNorm();
			}

			void initialize_timestep(scalar_t dt, UnsaturatedVariables* vars) override {
			m_dt = dt;
			MainVariable& main = get_var(vars);
			main.predictor = main.variable;
			main.transport_fluxes.setZero();
			main.velocity.setZero();
			m_solver.initialize_timestep(dt, get_var(vars).variable);
			}

			dfpmsolver::ParabolicDriverReturnCode restart_timestep(UnsaturatedVariables *vars) override {
			MainVariable& main = get_var(vars);
			m_solver.velocity() = main.velocity;
			dfpmsolver::ParabolicDriverReturnCode retcode = m_solver.restart_timestep(main.variable);
			if (retcode > dfpmsolver::ParabolicDriverReturnCode::NotConvergedYet)
			{
			main.velocity = m_solver.velocity();
			}
			return retcode;
			}

			dfpmsolver::ParabolicDriverOptions& get_options() {return m_solver.get_options();}

			const dfpmsolver::ParabolicDriverOptions& get_options() const {return m_solver.get_options();}

			scalar_t get_dt() {return m_dt;}

			protected:
			scalar_t m_norm_square_0 {-1};
			scalar_t m_dt {-1};
			EqT m_equation;
			SolverT m_solver;
			};

			//! \brief Traits struct for the SaturationTask class
			//!
			//! \internal
			struct SPECMICP_DLL_LOCAL SaturationTaskTraits {
			using EqT = SaturationEquation;
			using SolverT = dfpmsolver::ParabolicDriver<EqT>;
			using VariableBoxT = SaturationVariableBox;
			static constexpr EquationType equation_type {EquationType::Saturation};
			};

			//! \brief Wrapper for the saturation equation solver
			//!
			//! \internal
			class SPECMICP_DLL_LOCAL SaturationTask: public EquationTask<SaturationTask, SaturationTaskTraits>
			{
			using base = EquationTask<SaturationTask, SaturationTaskTraits>;
			public:

			SaturationTask(
			index_t component,
			mesh::Mesh1DPtr the_mesh,
			SaturationTaskTraits::VariableBoxT&& variables,
			const TransportConstraints& constraints
			):
			EquationTask<SaturationTask, SaturationTaskTraits>(
			component, the_mesh, variables, constraints)
			{}
			SaturationTask(
			index_t component,
			mesh::Mesh1DPtr the_mesh,
			SaturationTaskTraits::VariableBoxT&& variables,
			const TransportConstraints& constraints,
			scalar_t scaling
			):
			EquationTask<SaturationTask, SaturationTaskTraits>(
			component, the_mesh, variables, constraints, scaling)
			{}

			MainVariable& get_var(UnsaturatedVariables* vars) {
			return vars->get_liquid_saturation();
			}

			// Also take into account the solid total concentration
			scalar_t compute_squared_update(UnsaturatedVariables *vars) override {
			scalar_t solid_update = vars->get_solid_concentration(component()).velocity.squaredNorm();
			return solid_update + base::compute_squared_update(vars);
			}
			};

			//! \brief Traits struct for the LiquidAqueousTask class
			//!
			//! \internal
			struct SPECMICP_DLL_LOCAL LiquidAqueousTaskTraits {
			using EqT = AqueousTransportEquation;
			using SolverT = dfpmsolver::ParabolicDriver<EqT>;
			using VariableBoxT = LiquidAqueousComponentVariableBox;
			static constexpr EquationType equation_type {EquationType::LiquidAqueous};
			};

			//! \brief Wrapper for the liquid transport of aqueous component equation
			//!
			//! \internal
			class SPECMICP_DLL_LOCAL LiquidAqueousTask: public EquationTask<LiquidAqueousTask, LiquidAqueousTaskTraits>
			{
			using base = EquationTask<LiquidAqueousTask, LiquidAqueousTaskTraits>;

			public:
			LiquidAqueousTask(
			index_t component,
			mesh::Mesh1DPtr the_mesh,
			LiquidAqueousTaskTraits::VariableBoxT&& variables,
			const TransportConstraints& constraints
			):
			EquationTask<LiquidAqueousTask, LiquidAqueousTaskTraits>(
			component, the_mesh, variables, constraints)
			{}
			LiquidAqueousTask(
			index_t component,
			mesh::Mesh1DPtr the_mesh,
			LiquidAqueousTaskTraits::VariableBoxT&& variables,
			const TransportConstraints& constraints,
			scalar_t scaling
			):
			EquationTask<LiquidAqueousTask, LiquidAqueousTaskTraits>(
			component, the_mesh, variables, constraints, scaling)
			{}

			MainVariable& get_var(UnsaturatedVariables *vars) {
			return vars->get_aqueous_concentration(component());
			}

			// Also take into account the solid total concentration
			scalar_t compute_squared_update(UnsaturatedVariables *vars) override {
			scalar_t solid_update = vars->get_solid_concentration(component()).velocity.squaredNorm();
			return solid_update + base::compute_squared_update(vars);
			}
			};

			//! \brief Traits struct for the Pressure Task traits
			//!
			//! \internal
			struct SPECMICP_DLL_LOCAL PressureTaskTraits {
			using EqT = PressureEquation;
			using SolverT = dfpmsolver::ParabolicDriver<EqT>;
			using VariableBoxT = PressureVariableBox;
			static constexpr EquationType equation_type {EquationType::Pressure};
			};

			//! \brief Wrapper for the pressure equation solver
			//!
			//! \internal
			class SPECMICP_DLL_LOCAL PressureTask: public EquationTask<PressureTask, PressureTaskTraits>
			{
			using base = EquationTask<PressureTask, PressureTaskTraits>;

			public:
			PressureTask(
			index_t component,
			mesh::Mesh1DPtr the_mesh,
			PressureTaskTraits::VariableBoxT&& variables,
			const TransportConstraints& constraints
			):
			EquationTask<PressureTask, PressureTaskTraits>(
			component, the_mesh, variables, constraints)
			{}
			PressureTask(
			index_t component,
			mesh::Mesh1DPtr the_mesh,
			PressureTaskTraits::VariableBoxT&& variables,
			const TransportConstraints& constraints,
			scalar_t scaling
			):
			EquationTask<PressureTask, PressureTaskTraits>(
			component, the_mesh, variables, constraints, scaling)
			{}

			MainVariable& get_var(UnsaturatedVariables* vars) {
			return vars->get_pressure_main_variables(component());
			}
			};

			// ================================== //
			// //
			// UnsaturatedTransportStaggerImpl //
			// //
			// ================================== //

			// constructor
			// ===========

			UnsaturatedTransportStagger::UnsaturatedTransportStaggerImpl::UnsaturatedTransportStaggerImpl(
			UnsaturatedVariablesPtr variables,
			const TransportConstraints& constraints
			)
			{
			database::RawDatabasePtr raw_db = variables->get_database();
			mesh::Mesh1DPtr the_mesh = variables->get_mesh();

			// Saturation equations
			m_saturation_equation = make_unique<SaturationTask>(0, the_mesh,
			variables->get_saturation_variables(),
			constraints,
			variables->get_aqueous_scaling(0)
			);
			dfpmsolver::ParabolicDriverOptions& opts = m_saturation_equation->get_options();
			opts.residuals_tolerance = 1e-8;
			opts.absolute_tolerance = 1e-12;
			opts.step_tolerance = 1e-12;
			opts.threshold_stationary_point = 1e-12;
			opts.maximum_iterations = 500;
			opts.sparse_solver = sparse_solvers::SparseSolver::SparseLU;
			opts.quasi_newton = 3;
			opts.maximum_step_length = 0.1;
			opts.max_iterations_at_max_length = 5000;

			const index_t size = raw_db->nb_aqueous_components() + variables->nb_gas();
			m_equation_list.reserve(size);

			// Liquid aqueous diffusion-advection
			for (index_t id: raw_db->range_aqueous_component())
			{
			m_equation_list.emplace_back(
			make_unique<LiquidAqueousTask>(id, the_mesh,
			variables->get_liquid_aqueous_component_variables(id),
			constraints,
			variables->get_aqueous_scaling(id))
			);
			dfpmsolver::ParabolicDriverOptions& opts =
			static_cast<LiquidAqueousTask*>(m_equation_list[m_equation_list.size()-1].get())->get_options();
			opts.maximum_iterations = 500;
			opts.residuals_tolerance = 1e-8;
			opts.absolute_tolerance = 1e-12;
			opts.step_tolerance = 1e-12;
			opts.threshold_stationary_point = 1e-12;
			opts.sparse_solver = sparse_solvers::SparseSolver::SparseLU;
			opts.maximum_step_length= 1e6;
			}
			// Gaseous diffusion
			for (index_t id: raw_db->range_component())
			{
			if (variables->component_has_gas(id)) {
			m_equation_list.emplace_back(
			make_unique<PressureTask>(id, the_mesh,
			variables->get_pressure_variables(id),
			constraints,
			variables->get_gaseous_scaling(id)
			));
			dfpmsolver::ParabolicDriverOptions& opts =
			static_cast<PressureTask*>(m_equation_list[m_equation_list.size()-1].get())->get_options();
			opts.maximum_iterations = 500;
			opts.residuals_tolerance = 1e-8;
			opts.absolute_tolerance = 1e-12;
			opts.step_tolerance = 1e-12;
			opts.threshold_stationary_point = 1e-12;
			opts.sparse_solver = sparse_solvers::SparseSolver::SparseLU;
			opts.maximum_step_length= 1e6;
			}
			}
			}

			// Residuals
			// =========

			scalar_t
			UnsaturatedTransportStagger::UnsaturatedTransportStaggerImpl::get_residual(
			UnsaturatedVariables* vars
			)
			{
			scalar_t sum = m_saturation_equation->compute_squared_residual(vars);
			for (auto& task: m_equation_list) {
			sum += task->compute_squared_residual(vars);
			}
			auto norm = std::sqrt(sum);
			return norm;
			}

			scalar_t
			UnsaturatedTransportStagger::UnsaturatedTransportStaggerImpl::get_residual_0(
			UnsaturatedVariables* vars
			)
			{
			return m_norm_0;
			}

			scalar_t
			UnsaturatedTransportStagger::UnsaturatedTransportStaggerImpl::get_update(
			UnsaturatedVariables* vars
			)
			{
			scalar_t sum = m_saturation_equation->compute_squared_update(vars);
			for (auto& task: m_equation_list) {
			sum += task->compute_squared_update(vars);
			}
			return std::sqrt(sum);
			}

			// Solving the equations
			// =====================

			void
			UnsaturatedTransportStagger::UnsaturatedTransportStaggerImpl::initialize_timestep(
			scalar_t dt,
			UnsaturatedVariables* vars
			)
			{
			save_dt(dt);

			vars->get_advection_flux().set_zero();

			m_saturation_equation->initialize_timestep(dt, vars);
			scalar_t sum = m_saturation_equation->compute_squared_residual_0(vars);

			for (auto& task: m_equation_list) {
			task->initialize_timestep(dt, vars);
			sum += task->compute_squared_residual_0(vars);
			}
			m_norm_0 = 1.0;//std::sqrt(sum);
			}

			solver::StaggerReturnCode
			UnsaturatedTransportStagger::UnsaturatedTransportStaggerImpl::restart_timestep(
			UnsaturatedVariables* vars
			)
			{
			dfpmsolver::ParabolicDriverReturnCode retcode;

			bool flag_fail = false;
			bool flag_error_minimized = false;

			// Saturation equation
			retcode = m_saturation_equation->restart_timestep(vars);
			if (dfpmsolver::has_failed(retcode)) {
			WARNING << "Failed to solve saturation equation, return code :"
			<< to_string(retcode);
			flag_fail = true;
			}
			if (retcode == dfpmsolver::ParabolicDriverReturnCode::ErrorMinimized)
			flag_error_minimized = true;

			// Other equation
			if (not (flag_fail or flag_error_minimized))
			{
			vars->set_relative_variables();
			for (auto& task: m_equation_list) {
			retcode = task->restart_timestep(vars);
			if (dfpmsolver::has_failed(retcode)) {
			WARNING << "Equation of type '" << to_string(task->equation_type())
			<< "' for component " << task->component()
			<< " has failed with return code : "
			<< to_string(retcode);
			flag_fail = true;
			}
			flag_error_minimized = flag_error_minimized and
			(retcode == dfpmsolver::ParabolicDriverReturnCode::ErrorMinimized);
			}
			}

			// Return code
			solver::StaggerReturnCode return_code = solver::StaggerReturnCode::ResidualMinimized;
			if (flag_error_minimized)
			{
			return_code = solver::StaggerReturnCode::ErrorMinimized;
			}
			if (flag_fail)
			return_code = solver::StaggerReturnCode::UnknownError;

			return return_code;

			}

			// ================================ //
			// //
			// Helper functions //
			// //
			// ================================ //

			static std::string to_string(EquationType eq_type)
			{
			std::string str;
			switch (eq_type) {
			case EquationType::LiquidAqueous:
			str = "Liquid advection-diffusion";
			break;
			case EquationType::Pressure:
			str = "Gaseous diffusion";
			break;
			case EquationType::Saturation:
			str = "Saturation equation";
			break;
			default:
			str = "Unknown equation";
			break;
			}
			return str;
			}

			std::string to_string(dfpmsolver::ParabolicDriverReturnCode retcode)
			{
			using RetCode = dfpmsolver::ParabolicDriverReturnCode;
			std::string str;
			switch (retcode) {
			case RetCode::ResidualMinimized:
			str = "ResidualMinimized";
			break;
			case RetCode::ErrorMinimized:
			str = "ErrorMinimized";
			break;
			case RetCode::NotConvergedYet:
			str = "NotConvergedYet";
			break;
			case RetCode::StationaryPoint:
			str = "StationaryPoint";
			break;
			case RetCode::ErrorLinearSystem:
			str = "ErrorLinearSystem";
			break;
			case RetCode::MaxStepTakenTooManyTimes:
			str = "MaxStepTakenTooManyTimes";
			break;
			case RetCode::MaxIterations:
			str = "MaxIterations";
			break;
			case RetCode::LinesearchFailed:
			str = "LinesearchFailed";
			default:
			str = "Unknown return code";
			break;
			}
			return str;
			}

			} //end namespace unsaturated
			} //end namespace systems
			} //end namespace reactmicp
			} //end namespace specmicp

specmicp/src/reactmicp/systems/unsaturated/transport_stagger.cpp008df9c9036ereactmicp/unsaturated

transport_stagger.cppView Options

specmicp/src/reactmicp/systems/unsaturated/transport_stagger.cpp
008df9c9036ereactmicp/unsaturated

transport_stagger.cpp
View Options