diffusion.hpp
No OneTemporary
Actions

Subscribers

None

File Metadata

Created: Sun, Nov 10, 21:43

diffusion.hpp
View Options

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

	- Module : reactmicp/systems/diffusion
	- File : diffusion.hpp
	- Author : Fabien Georget

	Copyright (c) 2014, Fabien 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:
	* Redistributions of source code must retain the above copyright
	notice, this list of conditions and the following disclaimer.
	* 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.
	* Neither the name of the Princeton University 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 OWNER 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.

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

	#ifndef SPECMICP_REACTMICP_SYSTEMS_DIFFUSION_HPP
	#define SPECMICP_REACTMICP_SYSTEMS_DIFFUSION_HPP

	#include <vector>
	#include <memory>

	#include "database/data_container.hpp"

	#include "reactmicp/common.hpp"
	#include "reactmicp/micpsolvers/parabolicprogram.hpp"
	#include "reactmicp/meshes/mesh1d.hpp"

	#include "diffusion_numbering.hpp"
	#include "diffusion_secondary.hpp"


	namespace specmicp {
	namespace reactmicp {
	namespace systems {

	//! \brief Parameters for the diffusion system
	class DiffusionParameter
	{
	public:

	DiffusionParameter(double diffusion_coefficient, double porosity):
	m_diff(diffusion_coefficient), m_poro(porosity)
	{}

	//! \brief Density of water (kg/m^3)
	double density_water() {return 1e3;}

	//! \brief Diffusion coefficient (element by element) (m^2/s)
	double diffusion_coefficient(ind_t element) {return m_diff;}

	//! \brief Porosity (at a node (function of composition))
	double porosity(ind_t node) {return m_poro;}

	private:
	double m_diff;
	double m_poro;
	};

	//! \brief Saturated diffusion in reactive porous media
	class DiffusionProgram: public solvers::ParabolicProgram<DiffusionProgram>
	{

	public:
	using Variables = ParabolicVariables;
	const int no_equation = -1;

	DiffusionProgram(std::shared_ptr<mesh::Mesh1D> themesh,
	std::shared_ptr<database::DataContainer> thedatabase,
	std::shared_ptr<DiffusionParameter> parameters,
	diffusion::ListBoundaryCondition& list_bcs,
	Variables& var):
	m_mesh(themesh), m_data(thedatabase), m_param(parameters),
	m_ideq(themesh->nnodes(), thedatabase, list_bcs, var),
	m_second(themesh->nnodes(), thedatabase, m_ideq)
	{

	}

	DiffusionProgram(std::shared_ptr<mesh::Mesh1D> themesh,
	std::shared_ptr<database::DataContainer> thedatabase,
	std::shared_ptr<DiffusionParameter> parameters,
	diffusion::ListBoundaryCondition& list_bcs,
	const EquilibriumState &initialstate,
	Variables& var);


	//! \brief Return the number of degree of freedoms
	ind_t get_neq() const {return m_ideq.get_neq();}
	ind_t get_ndf() const {return m_ideq.get_ndf();}

	//! \brief Return the residuals
	void get_residuals(const ParabolicVariables& variable,
	Vector& residual);

	//! brief Return the jacobian
	void get_jacobian(ParabolicVariables& variable,
	Vector& residual,
	SparseMatrix& jacobian,
	double alphadt);

	//! \brief Return the index of the MiCP equations
	const std::vector<bool>& micp_index() const { return m_is_micp;}


	//! \brief Update the variables
	void update_variables(ParabolicVariables& x,
	const Eigen::VectorXd& predictor,
	const Eigen::VectorXd& update,
	double factor,
	double alpha_dt
	);
	//! \brief set the predictor
	void set_predictor(ParabolicVariables& x,
	Eigen::VectorXd& predictor,
	double alpha,
	double dt);

	//! \brief compute secondary concentrations
	bool hook_start_iteration(const ParabolicVariables& x, double norm_residual);

	//! \brief Return the maximum lambda allowed for the linesearch
	double maximum_lambda(const ParabolicVariables& x, const Eigen::VectorXd& update, double alphadt);

	//! \brief copy mineral amount from variables
	void copy_cp_variables(const ParabolicVariables& x,
	Eigen::VectorXd& moles_mineral);

	//! \brief projection step
	void projection(ParabolicVariables& x);

	// todo or not........
	// double update_equation(int id_eq,
	// double update,
	// double predictor,
	// double alphadt
	// )
	// {
	// // find _id eq ....
	// }

	private:

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

	// transport
	// ---------

	//! \brief Compute the residuals for the diffusion equations
	void residual_transport(const Variables& variable,
	Vector& residual);

	//! \brief Contribution of the diffusion operator to the residuals
	void element_residual_transport(ind_t element,
	const Variables& variable,
	Vector& residual);

	//! \brief Contribution of the minerals to the residuals of the diffusion equations
	double nodal_mineral_transient_term_transport(ind_t node,
	ind_t component,
	const Variables& variable);

	void element_residual_transport_component(
	ind_t element,
	int component,
	const Variables& variable,
	Vector& element_residual);

	// speciation
	// ----------

	//! \brief Compute the residual for the speciation problem
	void residual_speciation(const Variables& variable,
	Vector& residual);

	//! \brief Compute the residuals of the aqueous mass balances
	void nodal_residual_massbalance(ind_t node,
	const Variables& variable,
	Vector& residual);

	//! \brief compute the residual for the aqueous mass balance of 'component' at 'node'
	double nodal_component_residual_massbalance(ind_t node,
	int component,
	const Variables& variable);

	//! \brief Compute the residuals of the mineral speciation
	void nodal_residual_mineral(ind_t node,
	const Variables& variable,
	Vector& residual);

	//! \brief compute the residual for the stability equation of 'mineral' at 'node'
	double nodal_mineral_residual_mineral(ind_t node,
	int mineral,
	const Variables& variable);

	// Jacobian
	// ========

	// transport
	// ---------

	//! \brief Contribution of the diffusion equation to the jacobian
	void jacobian_transport(Variables& variable,
	Vector& residual,
	list_triplet_t& jacobian,
	double alphadt);

	//! \brief Element by element assembly procedure
	void element_jacobian_transport(ind_t element,
	Variables& variable,
	Vector& residual,
	list_triplet_t& jacobian,
	double alphadt);
	// speciation
	// ----------

	//! \brief Contribution of the speciation problem to the mass balance
	void jacobian_speciation(Variables &variable,
	list_triplet_t& jacobian,
	double alphadt);

	//! \brief Contribution of the mass balance at 'node' to the jacobian
	void nodal_jacobian_massbalance(ind_t node,
	const Variables& variable,
	list_triplet_t& jacobian,
	double alphadt);

	//! \brief Contribution of the mineral at equilibrium problem to the jacobian
	void nodal_jacobian_mineral(ind_t node,
	const Variables& variable,
	list_triplet_t& jacobian,
	double alphadt);

	//! \brief Finite difference jacobian for speciation
	void nodal_jacobian_speciation_fd(ind_t node,
	Variables &variable,
	list_triplet_t& jacobian,
	double alphadt);

	// Initialization
	// ==============

	void initialize_no_bc_with_init(const EquilibriumState& initialstate,
	ParabolicVariables& var);

	// Update
	// ======

	double update_diffusion(double update,
	double predictor,
	double alphadt)
	{
	return predictor + alphadt*update;
	}
	double update_massbalance(double update,
	double predictor,
	double alphadt)
	{
	return predictor + alphadt*update;
	}
	double update_mineral(double update,
	double predictor,
	double alphadt)
	{
	return predictor + alphadt*update;
	}

	// Member
	// =======

	std::shared_ptr<mesh::Mesh1D> m_mesh;
	std::shared_ptr<database::DataContainer> m_data;
	std::shared_ptr<DiffusionParameter> m_param;

	diffusion::DiffusionNumbering m_ideq;
	diffusion::DiffusionSecondaryVariables m_second;

	std::vector<bool> m_is_micp;

	//SparseVector m_boundary_condition;

	};

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

	#endif // SPECMICP_REACTMICP_SYSTEMS_DIFFUSION_HPP

diffusion.hppNo OneTemporaryActions

File Metadata

diffusion.hppView Options

Event Timeline

diffusion.hpp
No OneTemporary
Actions

diffusion.hpp
View Options