heat_transfer_model.cc
No OneTemporary
Actions

Subscribers

None

File Metadata

Created: Sat, Nov 9, 21:00

heat_transfer_model.cc
View Options

	/**
	* @file heat_transfer_model.cc
	*
	* @author Guillaume Anciaux <guillaume.anciaux@epfl.ch>
	* @author Lucas Frerot <lucas.frerot@epfl.ch>
	* @author David Simon Kammer <david.kammer@epfl.ch>
	* @author Srinivasa Babu Ramisetti <srinivasa.ramisetti@epfl.ch>
	* @author Nicolas Richart <nicolas.richart@epfl.ch>
	* @author Rui Wang <rui.wang@epfl.ch>
	*
	* @date creation: Sun May 01 2011
	* @date last modification: Mon Nov 30 2015
	*
	* @brief Implementation of HeatTransferModel class
	*
	* @section LICENSE
	*
	* Copyright (©) 2010-2012, 2014, 2015 EPFL (Ecole Polytechnique Fédérale de
	* Lausanne) Laboratory (LSMS - Laboratoire de Simulation en Mécanique des
	* Solides)
	*
	* Akantu is free software: you can redistribute it and/or modify it under the
	* terms of the GNU Lesser General Public License as published by the Free
	* Software Foundation, either version 3 of the License, or (at your option) any
	* later version.
	*
	* Akantu is distributed in the hope that it will be useful, but WITHOUT ANY
	* WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
	* A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more
	* details.
	*
	* You should have received a copy of the GNU Lesser General Public License
	* along with Akantu. If not, see <http://www.gnu.org/licenses/>.
	*
	*/

	/* -------------------------------------------------------------------------- */
	#include "aka_common.hh"
	#include "heat_transfer_model.hh"
	#include "group_manager_inline_impl.cc"
	#include "dumpable_inline_impl.hh"
	#include "aka_math.hh"
	#include "aka_common.hh"
	#include "fe_engine_template.hh"
	#include "mesh.hh"
	#include "static_communicator.hh"
	#include "parser.hh"
	#include "generalized_trapezoidal.hh"

	#ifdef AKANTU_USE_MUMPS
	#include "solver_mumps.hh"
	#endif

	#ifdef AKANTU_USE_IOHELPER
	# include "dumper_paraview.hh"
	# include "dumper_elemental_field.hh"
	# include "dumper_element_partition.hh"
	# include "dumper_internal_material_field.hh"
	#endif

	/* -------------------------------------------------------------------------- */
	__BEGIN_AKANTU__

	const HeatTransferModelOptions default_heat_transfer_model_options(_explicit_lumped_capacity);

	/* -------------------------------------------------------------------------- */
	HeatTransferModel::HeatTransferModel(Mesh & mesh,
	UInt dim,
	const ID & id,
	const MemoryID & memory_id) :
	Model(mesh, dim, id, memory_id), Parsable(_st_heat, id),
	integrator(NULL),
	conductivity_matrix(NULL),
	capacity_matrix(NULL),
	jacobian_matrix(NULL),
	temperature_gradient ("temperature_gradient", id),
	temperature_on_qpoints ("temperature_on_qpoints", id),
	conductivity_on_qpoints ("conductivity_on_qpoints", id),
	k_gradt_on_qpoints ("k_gradt_on_qpoints", id),
	int_bt_k_gT ("int_bt_k_gT", id),
	bt_k_gT ("bt_k_gT", id),
	conductivity(spatial_dimension, spatial_dimension),
	thermal_energy ("thermal_energy", id),
	solver(NULL),
	pbc_synch(NULL) {
	AKANTU_DEBUG_IN();

	createSynchronizerRegistry(this);

	std::stringstream sstr; sstr << id << ":fem";
	registerFEEngineObject<MyFEEngineType>(sstr.str(), mesh,spatial_dimension);

	this->temperature= NULL;
	this->residual = NULL;
	this->blocked_dofs = NULL;

	#ifdef AKANTU_USE_IOHELPER
	this->mesh.registerDumper<DumperParaview>("paraview_all", id, true);
	this->mesh.addDumpMesh(mesh, spatial_dimension, _not_ghost, _ek_regular);
	#endif

	this->registerParam("conductivity" , conductivity , _pat_parsmod);
	this->registerParam("conductivity_variation", conductivity_variation, 0., _pat_parsmod);
	this->registerParam("temperature_reference" , T_ref , 0., _pat_parsmod);
	this->registerParam("capacity" , capacity , _pat_parsmod);
	this->registerParam("density" , density , _pat_parsmod);

	AKANTU_DEBUG_OUT();
	}

	/* -------------------------------------------------------------------------- */
	void HeatTransferModel::initModel() {
	getFEEngine().initShapeFunctions(_not_ghost);
	getFEEngine().initShapeFunctions(_ghost);
	}

	/* -------------------------------------------------------------------------- */
	void HeatTransferModel::initParallel(MeshPartition * partition,
	DataAccessor * data_accessor) {
	AKANTU_DEBUG_IN();

	if (data_accessor == NULL) data_accessor = this;
	Synchronizer & synch_parallel = createParallelSynch(partition,data_accessor);

	synch_registry->registerSynchronizer(synch_parallel, _gst_htm_capacity);
	synch_registry->registerSynchronizer(synch_parallel, _gst_htm_temperature);
	synch_registry->registerSynchronizer(synch_parallel, _gst_htm_gradient_temperature);

	AKANTU_DEBUG_OUT();
	}

	/* -------------------------------------------------------------------------- */
	void HeatTransferModel::initPBC() {
	AKANTU_DEBUG_IN();

	Model::initPBC();
	pbc_synch = new PBCSynchronizer(pbc_pair);

	synch_registry->registerSynchronizer(*pbc_synch, _gst_htm_capacity);
	synch_registry->registerSynchronizer(*pbc_synch, _gst_htm_temperature);
	changeLocalEquationNumberForPBC(pbc_pair,1);

	// as long as there are ones on the diagonal of the matrix, we can put boudandary true for slaves
	std::map<UInt, UInt>::iterator it = pbc_pair.begin();
	std::map<UInt, UInt>::iterator end = pbc_pair.end();
	while(it != end) {
	(blocked_dofs)((it).first,0) = true;
	++it;
	}



	AKANTU_DEBUG_OUT();
	}

	/* -------------------------------------------------------------------------- */
	void HeatTransferModel::initArrays() {
	AKANTU_DEBUG_IN();

	UInt nb_nodes = getFEEngine().getMesh().getNbNodes();

	std::stringstream sstr_temp; sstr_temp << Model::id << ":temperature";
	std::stringstream sstr_temp_rate; sstr_temp_rate << Model::id << ":temperature_rate";
	std::stringstream sstr_inc; sstr_inc << Model::id << ":increment";
	std::stringstream sstr_ext_flx; sstr_ext_flx << Model::id << ":external_flux";
	std::stringstream sstr_residual; sstr_residual << Model::id << ":residual";
	std::stringstream sstr_lump; sstr_lump << Model::id << ":lumped";
	std::stringstream sstr_boun; sstr_boun << Model::id << ":blocked_dofs";

	temperature = &(alloc<Real>(sstr_temp.str(), nb_nodes, 1, REAL_INIT_VALUE));
	temperature_rate = &(alloc<Real>(sstr_temp_rate.str(), nb_nodes, 1, REAL_INIT_VALUE));
	increment = &(alloc<Real>(sstr_inc.str(), nb_nodes, 1, REAL_INIT_VALUE));
	external_heat_rate = &(alloc<Real>(sstr_ext_flx.str(), nb_nodes, 1, REAL_INIT_VALUE));
	residual = &(alloc<Real>(sstr_residual.str(), nb_nodes, 1, REAL_INIT_VALUE));
	capacity_lumped = &(alloc<Real>(sstr_lump.str(), nb_nodes, 1, REAL_INIT_VALUE));
	blocked_dofs = &(alloc<bool>(sstr_boun.str(), nb_nodes, 1, false));

	Mesh::ConnectivityTypeList::const_iterator it;

	/* -------------------------------------------------------------------------- */
	// byelementtype vectors
	getFEEngine().getMesh().initElementTypeMapArray(temperature_on_qpoints,
	1,
	spatial_dimension);

	getFEEngine().getMesh().initElementTypeMapArray(temperature_gradient,
	spatial_dimension,
	spatial_dimension);

	getFEEngine().getMesh().initElementTypeMapArray(conductivity_on_qpoints,
	spatial_dimension*spatial_dimension,
	spatial_dimension);

	getFEEngine().getMesh().initElementTypeMapArray(k_gradt_on_qpoints,
	spatial_dimension,
	spatial_dimension);

	getFEEngine().getMesh().initElementTypeMapArray(bt_k_gT,
	1,
	spatial_dimension,
	true);

	getFEEngine().getMesh().initElementTypeMapArray(int_bt_k_gT,
	1,
	spatial_dimension,
	true);

	getFEEngine().getMesh().initElementTypeMapArray(thermal_energy,
	1,
	spatial_dimension);


	for(UInt g = _not_ghost; g <= _ghost; ++g) {
	GhostType gt = (GhostType) g;

	const Mesh::ConnectivityTypeList & type_list =
	getFEEngine().getMesh().getConnectivityTypeList(gt);

	for(it = type_list.begin(); it != type_list.end(); ++it) {
	if(Mesh::getSpatialDimension(*it) != spatial_dimension) continue;
	UInt nb_element = getFEEngine().getMesh().getNbElement(*it, gt);
	UInt nb_quad_points = this->getFEEngine().getNbIntegrationPoints(it, gt) nb_element;

	temperature_on_qpoints(*it, gt).resize(nb_quad_points);
	temperature_on_qpoints(*it, gt).clear();

	temperature_gradient(*it, gt).resize(nb_quad_points);
	temperature_gradient(*it, gt).clear();

	conductivity_on_qpoints(*it, gt).resize(nb_quad_points);
	conductivity_on_qpoints(*it, gt).clear();

	k_gradt_on_qpoints(*it, gt).resize(nb_quad_points);
	k_gradt_on_qpoints(*it, gt).clear();

	bt_k_gT(*it, gt).resize(nb_quad_points);
	bt_k_gT(*it, gt).clear();

	int_bt_k_gT(*it, gt).resize(nb_element);
	int_bt_k_gT(*it, gt).clear();

	thermal_energy(*it, gt).resize(nb_element);
	thermal_energy(*it, gt).clear();
	}
	}

	/* -------------------------------------------------------------------------- */
	dof_synchronizer = new DOFSynchronizer(getFEEngine().getMesh(),1);
	dof_synchronizer->initLocalDOFEquationNumbers();
	dof_synchronizer->initGlobalDOFEquationNumbers();

	AKANTU_DEBUG_OUT();
	}

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

	void HeatTransferModel::initSolver(__attribute__((unused)) SolverOptions & options) {
	#if !defined(AKANTU_USE_MUMPS) // or other solver in the future \todo add AKANTU_HAS_SOLVER in CMake
	AKANTU_DEBUG_ERROR("You should at least activate one solver.");
	#else
	UInt nb_global_nodes = mesh.getNbGlobalNodes();

	delete jacobian_matrix;
	std::stringstream sstr; sstr << Memory::id << ":jacobian_matrix";
	jacobian_matrix = new SparseMatrix(nb_global_nodes, _symmetric, sstr.str(), memory_id);

	jacobian_matrix->buildProfile(mesh, *dof_synchronizer, 1);

	delete conductivity_matrix;
	std::stringstream sstr_sti; sstr_sti << Memory::id << ":conductivity_matrix";
	conductivity_matrix = new SparseMatrix(*jacobian_matrix, sstr_sti.str(), memory_id);

	#ifdef AKANTU_USE_MUMPS
	std::stringstream sstr_solv; sstr_solv << Memory::id << ":solver";
	solver = new SolverMumps(*jacobian_matrix, sstr_solv.str());

	dof_synchronizer->initScatterGatherCommunicationScheme();
	#else
	AKANTU_DEBUG_ERROR("You should at least activate one solver.");
	#endif //AKANTU_USE_MUMPS

	if(solver)
	solver->initialize(options);
	#endif //AKANTU_HAS_SOLVER
	}

	/* -------------------------------------------------------------------------- */
	void HeatTransferModel::initImplicit(bool dynamic, SolverOptions & solver_options) {
	AKANTU_DEBUG_IN();

	method = dynamic ? _implicit_dynamic : _static;
	initSolver(solver_options);

	if(method == _implicit_dynamic) {
	if(integrator) delete integrator;
	integrator = new TrapezoidalRule1();
	}

	AKANTU_DEBUG_OUT();
	}

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

	HeatTransferModel::~HeatTransferModel()
	{
	AKANTU_DEBUG_IN();
	if (integrator) delete integrator ;
	if (conductivity_matrix) delete conductivity_matrix;
	if (capacity_matrix) delete capacity_matrix ;
	if (jacobian_matrix) delete jacobian_matrix ;
	if (solver) delete solver ;

	delete pbc_synch;

	AKANTU_DEBUG_OUT();
	}

	/* -------------------------------------------------------------------------- */
	void HeatTransferModel::assembleCapacityLumped(const GhostType & ghost_type) {
	AKANTU_DEBUG_IN();

	FEEngine & fem = getFEEngine();

	const Mesh::ConnectivityTypeList & type_list
	= fem.getMesh().getConnectivityTypeList(ghost_type);

	Mesh::ConnectivityTypeList::const_iterator it;
	for(it = type_list.begin(); it != type_list.end(); ++it)
	{
	if(Mesh::getSpatialDimension(*it) != spatial_dimension) continue;

	UInt nb_element = getFEEngine().getMesh().getNbElement(*it,ghost_type);
	UInt nb_quadrature_points = getFEEngine().getNbIntegrationPoints(*it, ghost_type);

	Array<Real> rho_1 (nb_element * nb_quadrature_points,1, capacity * density);
	fem.assembleFieldLumped(rho_1,1,*capacity_lumped,
	dof_synchronizer->getLocalDOFEquationNumbers(),
	*it, ghost_type);
	}

	AKANTU_DEBUG_OUT();
	}

	/* -------------------------------------------------------------------------- */
	void HeatTransferModel::assembleCapacityLumped() {
	AKANTU_DEBUG_IN();

	capacity_lumped->clear();

	assembleCapacityLumped(_not_ghost);
	assembleCapacityLumped(_ghost);

	getSynchronizerRegistry().synchronize(_gst_htm_capacity);

	AKANTU_DEBUG_OUT();
	}

	/* -------------------------------------------------------------------------- */
	void HeatTransferModel::updateResidual(bool compute_conductivity) {
	AKANTU_DEBUG_IN();
	/// @f$ r = q_{ext} - q_{int} - C \dot T @f$

	// start synchronization
	synch_registry->asynchronousSynchronize(_gst_htm_temperature);
	// finalize communications
	synch_registry->waitEndSynchronize(_gst_htm_temperature);

	//clear the array
	/// first @f$ r = q_{ext} @f$
	// residual->clear();
	residual->copy(*external_heat_rate);

	/// then @f$ r -= q_{int} @f$
	// update the not ghost ones

	updateResidual(_not_ghost);
	// update for the received ghosts
	updateResidual(_ghost);

	/* if (method == _explicit_lumped_capacity) {
	this->solveExplicitLumped();
	}*/


	AKANTU_DEBUG_OUT();
	}

	/* -------------------------------------------------------------------------- */
	void HeatTransferModel::assembleConductivityMatrix(bool compute_conductivity) {
	AKANTU_DEBUG_IN();

	AKANTU_DEBUG_INFO("Assemble the new stiffness matrix.");

	conductivity_matrix->clear();

	switch(mesh.getSpatialDimension()) {
	case 1: this->assembleConductivityMatrix<1>(_not_ghost,compute_conductivity); break;
	case 2: this->assembleConductivityMatrix<2>(_not_ghost,compute_conductivity); break;
	case 3: this->assembleConductivityMatrix<3>(_not_ghost,compute_conductivity); break;
	}

	AKANTU_DEBUG_OUT();
	}

	/* -------------------------------------------------------------------------- */
	template <UInt dim>
	void HeatTransferModel::assembleConductivityMatrix(const GhostType & ghost_type,bool compute_conductivity) {
	AKANTU_DEBUG_IN();

	Mesh & mesh = this->getFEEngine().getMesh();
	Mesh::type_iterator it = mesh.firstType(spatial_dimension, ghost_type);
	Mesh::type_iterator last_type = mesh.lastType(spatial_dimension, ghost_type);
	for(; it != last_type; ++it) {
	this->assembleConductivityMatrix<dim>(*it, ghost_type,compute_conductivity);
	}

	AKANTU_DEBUG_OUT();
	}

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


	/* -------------------------------------------------------------------------- */
	template <UInt dim>
	void HeatTransferModel::assembleConductivityMatrix(const ElementType & type,
	const GhostType & ghost_type,
	bool compute_conductivity) {
	AKANTU_DEBUG_IN();

	SparseMatrix & K = *conductivity_matrix;
	const Array<Real> & shapes_derivatives = this->getFEEngine().getShapesDerivatives(type, ghost_type);

	UInt nb_element = mesh.getNbElement(type, ghost_type);
	UInt nb_nodes_per_element = Mesh::getNbNodesPerElement(type);
	UInt nb_quadrature_points = getFEEngine().getNbIntegrationPoints(type, ghost_type);

	/// compute @f$\mathbf{B}^t * \mathbf{D} * \mathbf{B}@f$
	UInt bt_d_b_size = nb_nodes_per_element;

	Array<Real> * bt_d_b = new Array<Real>(nb_element * nb_quadrature_points,
	bt_d_b_size * bt_d_b_size,
	"B^tDB");

	Matrix<Real> Bt_D(nb_nodes_per_element, dim);

	Array<Real>::const_iterator< Matrix<Real> > shapes_derivatives_it = shapes_derivatives.begin(dim, nb_nodes_per_element);

	Array<Real>::iterator< Matrix<Real> > Bt_D_B_it = bt_d_b->begin(bt_d_b_size, bt_d_b_size);

	if (compute_conductivity)
	this->computeConductivityOnQuadPoints(ghost_type);
	Array<Real>::iterator< Matrix<Real> > D_it = conductivity_on_qpoints(type, ghost_type).begin(dim, dim);
	Array<Real>::iterator< Matrix<Real> > D_end = conductivity_on_qpoints(type, ghost_type).end(dim, dim);

	for (; D_it != D_end; ++D_it, ++Bt_D_B_it, ++shapes_derivatives_it) {
	Matrix<Real> & D = *D_it;
	const Matrix<Real> & B = *shapes_derivatives_it;
	Matrix<Real> & Bt_D_B = *Bt_D_B_it;

	Bt_D.mul<true, false>(B, D);
	Bt_D_B.mul<false, false>(Bt_D, B);
	}

	/// compute @f$ k_e = \int_e \mathbf{B}^t * \mathbf{D} * \mathbf{B}@f$
	Array<Real> * K_e = new Array<Real>(nb_element,
	bt_d_b_size * bt_d_b_size,
	"K_e");

	this->getFEEngine().integrate(bt_d_b, K_e,
	bt_d_b_size * bt_d_b_size,
	type, ghost_type);

	delete bt_d_b;

	this->getFEEngine().assembleMatrix(*K_e, K, 1, type, ghost_type);
	delete K_e;

	AKANTU_DEBUG_OUT();
	}

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

	void HeatTransferModel::updateResidualInternal() {

	AKANTU_DEBUG_IN();

	AKANTU_DEBUG_INFO("Update the residual");
	// f = q_ext - q_int - Mddot = q - Mddot;

	if(method != _static) {
	// f -= Mddot
	if(capacity_matrix) {
	// if full mass_matrix
	Array<Real> * Mddot = new Array<Real>(*temperature_rate, true, "Mddot");
	Mddot = *capacity_matrix;
	residual -= Mddot;
	delete Mddot;
	} else if (capacity_lumped) {

	// else lumped mass
	UInt nb_nodes = temperature_rate->getSize();
	UInt nb_degree_of_freedom = temperature_rate->getNbComponent();

	Real * capacity_val = capacity_lumped->storage();
	Real * temp_rate_val = temperature_rate->storage();
	Real * res_val = residual->storage();
	bool * blocked_dofs_val = blocked_dofs->storage();

	for (UInt n = 0; n < nb_nodes * nb_degree_of_freedom; ++n) {
	if(!(*blocked_dofs_val)) {
	res_val -= temp_rate_val * *capacity_val;
	}
	blocked_dofs_val++;
	res_val++;
	capacity_val++;
	temp_rate_val++;
	}
	} else {
	AKANTU_DEBUG_ERROR("No function called to assemble the mass matrix.");
	}
	}

	AKANTU_DEBUG_OUT();
	}

	/* -------------------------------------------------------------------------- */
	void HeatTransferModel::solveStatic() {
	AKANTU_DEBUG_IN();

	AKANTU_DEBUG_INFO("Solving Ku = f");
	AKANTU_DEBUG_ASSERT(conductivity_matrix != NULL,
	"You should first initialize the implicit solver and assemble the stiffness matrix");

	UInt nb_nodes = temperature->getSize();
	UInt nb_degree_of_freedom = temperature->getNbComponent() * nb_nodes;

	jacobian_matrix->copyContent(*conductivity_matrix);
	jacobian_matrix->applyBoundary(*blocked_dofs);

	increment->clear();

	solver->setRHS(*residual);
	solver->factorize();
	solver->solve(*increment);

	Real * increment_val = increment->storage();
	Real * temperature_val = temperature->storage();
	bool * blocked_dofs_val = blocked_dofs->storage();

	for (UInt j = 0; j < nb_degree_of_freedom;
	++j, ++temperature_val, ++increment_val, ++blocked_dofs_val) {
	if (!(*blocked_dofs_val)) {
	temperature_val += increment_val;
	}
	else {
	*increment_val = 0.0;
	}
	}

	AKANTU_DEBUG_OUT();
	}

	/* -------------------------------------------------------------------------- */
	void HeatTransferModel::computeConductivityOnQuadPoints(const GhostType & ghost_type) {
	const Mesh::ConnectivityTypeList & type_list =
	this->getFEEngine().getMesh().getConnectivityTypeList(ghost_type);
	Mesh::ConnectivityTypeList::const_iterator it;

	for(it = type_list.begin(); it != type_list.end(); ++it) {
	if(Mesh::getSpatialDimension(*it) != spatial_dimension) continue;

	Array<Real> & temperature_interpolated = temperature_on_qpoints(*it, ghost_type);

	//compute the temperature on quadrature points
	this->getFEEngine().interpolateOnIntegrationPoints(*temperature,
	temperature_interpolated,
	1 ,*it,ghost_type);

	Array<Real>::matrix_iterator C_it =
	conductivity_on_qpoints(*it, ghost_type).begin(spatial_dimension, spatial_dimension);
	Array<Real>::matrix_iterator C_end =
	conductivity_on_qpoints(*it, ghost_type).end(spatial_dimension, spatial_dimension);

	Array<Real>::iterator<Real> T_it = temperature_interpolated.begin();

	for (;C_it != C_end; ++C_it, ++T_it) {
	Matrix<Real> & C = *C_it;
	Real & T = *T_it;
	C = conductivity;

	Matrix<Real> variation(spatial_dimension, spatial_dimension, conductivity_variation * (T - T_ref));
	C += conductivity_variation;
	}
	}
	AKANTU_DEBUG_OUT();
	}
	/* -------------------------------------------------------------------------- */
	void HeatTransferModel::computeKgradT(const GhostType & ghost_type,bool compute_conductivity) {

	if (compute_conductivity)
	computeConductivityOnQuadPoints(ghost_type);

	const Mesh::ConnectivityTypeList & type_list =
	this->getFEEngine().getMesh().getConnectivityTypeList(ghost_type);
	Mesh::ConnectivityTypeList::const_iterator it;

	for(it = type_list.begin(); it != type_list.end(); ++it) {
	const ElementType & type = *it;
	if(Mesh::getSpatialDimension(*it) != spatial_dimension) continue;

	Array<Real> & gradient = temperature_gradient(*it, ghost_type);
	this->getFEEngine().gradientOnIntegrationPoints(*temperature,
	gradient,
	1 ,*it, ghost_type);

	Array<Real>::matrix_iterator C_it =
	conductivity_on_qpoints(*it, ghost_type).begin(spatial_dimension, spatial_dimension);

	Array<Real>::vector_iterator BT_it = gradient.begin(spatial_dimension);

	Array<Real>::vector_iterator k_BT_it =
	k_gradt_on_qpoints(type, ghost_type).begin(spatial_dimension);
	Array<Real>::vector_iterator k_BT_end =
	k_gradt_on_qpoints(type, ghost_type).end(spatial_dimension);

	for (;k_BT_it != k_BT_end; ++k_BT_it, ++BT_it, ++C_it) {
	Vector<Real> & k_BT = *k_BT_it;
	Vector<Real> & BT = *BT_it;
	Matrix<Real> & C = *C_it;

	k_BT.mul<false>(C, BT);
	}
	}

	AKANTU_DEBUG_OUT();
	}

	/* -------------------------------------------------------------------------- */
	void HeatTransferModel::updateResidual(const GhostType & ghost_type, bool compute_conductivity) {
	AKANTU_DEBUG_IN();

	const Mesh::ConnectivityTypeList & type_list =
	this->getFEEngine().getMesh().getConnectivityTypeList(ghost_type);
	Mesh::ConnectivityTypeList::const_iterator it;

	for(it = type_list.begin(); it != type_list.end(); ++it) {
	if(Mesh::getSpatialDimension(*it) != spatial_dimension) continue;

	Array<Real> & shapes_derivatives =
	const_cast<Array<Real> &>(getFEEngine().getShapesDerivatives(*it,ghost_type));

	UInt nb_nodes_per_element = Mesh::getNbNodesPerElement(*it);

	// compute k \grad T
	computeKgradT(ghost_type,compute_conductivity);

	Array<Real>::vector_iterator k_BT_it =
	k_gradt_on_qpoints(*it,ghost_type).begin(spatial_dimension);

	Array<Real>::matrix_iterator B_it =
	shapes_derivatives.begin(spatial_dimension, nb_nodes_per_element);

	Array<Real>::vector_iterator Bt_k_BT_it =
	bt_k_gT(*it,ghost_type).begin(nb_nodes_per_element);

	Array<Real>::vector_iterator Bt_k_BT_end =
	bt_k_gT(*it,ghost_type).end(nb_nodes_per_element);


	for (;Bt_k_BT_it != Bt_k_BT_end; ++Bt_k_BT_it, ++B_it, ++k_BT_it) {
	Vector<Real> & k_BT = *k_BT_it;
	Vector<Real> & Bt_k_BT = *Bt_k_BT_it;
	Matrix<Real> & B = *B_it;

	Bt_k_BT.mul<true>(B, k_BT);
	}

	this->getFEEngine().integrate(bt_k_gT(*it,ghost_type),
	int_bt_k_gT(*it,ghost_type),
	nb_nodes_per_element, *it,ghost_type);

	this->getFEEngine().assembleArray(int_bt_k_gT(it,ghost_type), residual,
	dof_synchronizer->getLocalDOFEquationNumbers(),
	1, *it,ghost_type, empty_filter, -1);
	}

	AKANTU_DEBUG_OUT();
	}

	/* -------------------------------------------------------------------------- */
	void HeatTransferModel::solveExplicitLumped() {
	AKANTU_DEBUG_IN();

	/// finally @f$ r -= C \dot T @f$
	// lumped C
	UInt nb_nodes = temperature_rate->getSize();
	UInt nb_degree_of_freedom = temperature_rate->getNbComponent();

	Real * capacity_val = capacity_lumped->storage();
	Real * temp_rate_val = temperature_rate->storage();
	Real * res_val = residual->storage();
	bool * blocked_dofs_val = blocked_dofs->storage();

	for (UInt n = 0; n < nb_nodes * nb_degree_of_freedom; ++n) {
	if(!(*blocked_dofs_val)) {
	res_val -= capacity_val * *temp_rate_val;
	}
	blocked_dofs_val++;
	res_val++;
	capacity_val++;
	temp_rate_val++;
	}

	#ifndef AKANTU_NDEBUG
	getSynchronizerRegistry().synchronize(akantu::_gst_htm_gradient_temperature);
	#endif

	capacity_val = capacity_lumped->storage();
	res_val = residual->storage();
	blocked_dofs_val = blocked_dofs->storage();
	Real * inc = increment->storage();

	for (UInt n = 0; n < nb_nodes * nb_degree_of_freedom; ++n) {
	if(!(*blocked_dofs_val)) {
	inc = (res_val / *capacity_val);
	}
	res_val++;
	blocked_dofs_val++;
	inc++;
	capacity_val++;
	}

	AKANTU_DEBUG_OUT();
	}

	/* -------------------------------------------------------------------------- */
	void HeatTransferModel::explicitPred() {
	AKANTU_DEBUG_IN();

	AKANTU_DEBUG_ASSERT(integrator, "integrator was not instanciated");
	integrator->integrationSchemePred(time_step,
	*temperature,
	*temperature_rate,
	*blocked_dofs);

	UInt nb_nodes = temperature->getSize();
	UInt nb_degree_of_freedom = temperature->getNbComponent();

	Real * temp = temperature->storage();
	for (UInt n = 0; n < nb_nodes * nb_degree_of_freedom; ++n, ++temp)
	if(temp < 0.) temp = 0.;

	AKANTU_DEBUG_OUT();
	}

	/* -------------------------------------------------------------------------- */
	void HeatTransferModel::explicitCorr() {
	AKANTU_DEBUG_IN();

	integrator->integrationSchemeCorrTempRate(time_step,
	*temperature,
	*temperature_rate,
	*blocked_dofs,
	*increment);

	AKANTU_DEBUG_OUT();
	}


	/* -------------------------------------------------------------------------- */
	void HeatTransferModel::implicitPred(){
	AKANTU_DEBUG_IN();

	if(method == _implicit_dynamic)
	integrator->integrationSchemePred(time_step,
	*temperature,
	*temperature_rate,
	*blocked_dofs);

	AKANTU_DEBUG_OUT();
	}
	/* -------------------------------------------------------------------------- */

	void HeatTransferModel::implicitCorr(){

	AKANTU_DEBUG_IN();

	if(method == _implicit_dynamic) {
	integrator->integrationSchemeCorrTemp(time_step,
	*temperature,
	*temperature_rate,
	*blocked_dofs,
	*increment);
	} else {
	UInt nb_nodes = temperature->getSize();
	UInt nb_degree_of_freedom = temperature->getNbComponent() * nb_nodes;

	Real * incr_val = increment->storage();
	Real * temp_val = temperature->storage();
	bool * boun_val = blocked_dofs->storage();

	for (UInt j = 0; j < nb_degree_of_freedom; ++j, ++temp_val, ++incr_val, ++boun_val){
	incr_val = (1. - *boun_val);
	temp_val += incr_val;
	}
	}

	AKANTU_DEBUG_OUT();
	}

	/* -------------------------------------------------------------------------- */
	Real HeatTransferModel::getStableTimeStep()
	{
	AKANTU_DEBUG_IN();

	Real el_size;
	Real min_el_size = std::numeric_limits<Real>::max();
	Real conductivitymax = conductivity(0, 0);

	//get the biggest parameter from k11 until k33//
	for(UInt i = 0; i < spatial_dimension; i++)
	for(UInt j = 0; j < spatial_dimension; j++)
	conductivitymax = std::max(conductivity(i, j), conductivitymax);


	const Mesh::ConnectivityTypeList & type_list =
	getFEEngine().getMesh().getConnectivityTypeList();
	Mesh::ConnectivityTypeList::const_iterator it;
	for(it = type_list.begin(); it != type_list.end(); ++it) {
	if(getFEEngine().getMesh().getSpatialDimension(*it) != spatial_dimension)
	continue;

	UInt nb_nodes_per_element = getFEEngine().getMesh().getNbNodesPerElement(*it);

	Array<Real> coord(0, nb_nodes_per_element*spatial_dimension);
	FEEngine::extractNodalToElementField(getFEEngine().getMesh(), getFEEngine().getMesh().getNodes(),
	coord, *it, _not_ghost);

	Array<Real>::matrix_iterator el_coord = coord.begin(spatial_dimension, nb_nodes_per_element);
	UInt nb_element = getFEEngine().getMesh().getNbElement(*it);

	for (UInt el = 0; el < nb_element; ++el, ++el_coord) {
	el_size = getFEEngine().getElementInradius(el_coord, it);
	min_el_size = std::min(min_el_size, el_size);
	}

	AKANTU_DEBUG_INFO("The minimum element size : " << min_el_size
	<< " and the max conductivity is : "
	<< conductivitymax);
	}

	Real min_dt = 2 * min_el_size * min_el_size * density
	* capacity / conductivitymax;

	StaticCommunicator::getStaticCommunicator().allReduce(&min_dt, 1, _so_min);

	AKANTU_DEBUG_OUT();

	return min_dt;
	}

	/* -------------------------------------------------------------------------- */
	void HeatTransferModel::readMaterials() {
	std::pair<Parser::const_section_iterator, Parser::const_section_iterator>
	sub_sect = this->parser->getSubSections(_st_heat);

	Parser::const_section_iterator it = sub_sect.first;
	const ParserSection & section = *it;

	this->parseSection(section);
	}

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

	void HeatTransferModel::initFull(const ModelOptions & options){
	Model::initFull(options);

	readMaterials();

	const HeatTransferModelOptions & my_options
	= dynamic_cast<const HeatTransferModelOptions &>(options);

	method = my_options.analysis_method;

	//initialize the vectors
	initArrays();
	temperature->clear();
	temperature_rate->clear();
	external_heat_rate->clear();


	// initialize pbc
	if(pbc_pair.size()!=0)
	initPBC();

	if (method == _explicit_lumped_capacity){
	integrator = new ForwardEuler();
	}
	if (method == _implicit_dynamic) {
	initImplicit(true);
	}
	if (method == _static) {
	initImplicit(false);
	}
	}
	/* -------------------------------------------------------------------------- */
	void HeatTransferModel::assembleCapacity() {
	AKANTU_DEBUG_IN();

	if(!capacity_matrix) {
	std::stringstream sstr; sstr << id << ":capacity_matrix";
	capacity_matrix = new SparseMatrix(*jacobian_matrix, sstr.str(), memory_id);
	}

	assembleCapacity(_not_ghost);
	// assembleMass(_ghost);

	AKANTU_DEBUG_OUT();
	}

	/* -------------------------------------------------------------------------- */
	void HeatTransferModel::assembleCapacity(GhostType ghost_type) {
	AKANTU_DEBUG_IN();

	MyFEEngineType & fem = getFEEngineClass<MyFEEngineType>();

	Array<Real> rho_1(0,1);
	//UInt nb_element;
	capacity_matrix->clear();

	Mesh::type_iterator it = mesh.firstType(spatial_dimension, ghost_type);
	Mesh::type_iterator end = mesh.lastType(spatial_dimension, ghost_type);
	for(; it != end; ++it) {
	ElementType type = *it;

	computeRho(rho_1, type, ghost_type);
	fem.assembleFieldMatrix(rho_1, 1, *capacity_matrix, type, ghost_type);
	}

	AKANTU_DEBUG_OUT();
	}

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

	void HeatTransferModel::computeRho(Array<Real> & rho,
	ElementType type,
	GhostType ghost_type) {
	AKANTU_DEBUG_IN();

	FEEngine & fem = this->getFEEngine();
	UInt nb_element = fem.getMesh().getNbElement(type,ghost_type);
	UInt nb_quadrature_points = fem.getNbIntegrationPoints(type, ghost_type);

	rho.resize(nb_element * nb_quadrature_points);
	Real * rho_1_val = rho.storage();

	/// compute @f$ rho @f$ for each nodes of each element
	for (UInt el = 0; el < nb_element; ++el) {

	for (UInt n = 0; n < nb_quadrature_points; ++n) {
	*rho_1_val++ = this->capacity;
	}
	}

	AKANTU_DEBUG_OUT();
	}

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

	template<>
	bool HeatTransferModel::testConvergence<_scc_increment>(Real tolerance, Real & error){
	AKANTU_DEBUG_IN();

	UInt nb_nodes = temperature->getSize();
	UInt nb_degree_of_freedom = temperature->getNbComponent();

	error = 0;
	Real norm[2] = {0., 0.};
	Real * increment_val = increment->storage();
	bool * blocked_dofs_val = blocked_dofs->storage();
	Real * temperature_val = temperature->storage();

	for (UInt n = 0; n < nb_nodes; ++n) {
	bool is_local_node = mesh.isLocalOrMasterNode(n);
	for (UInt d = 0; d < nb_degree_of_freedom; ++d) {
	if(!(*blocked_dofs_val) && is_local_node) {
	norm[0] += increment_val *increment_val;
	norm[1] += temperature_val *temperature_val;
	}
	blocked_dofs_val++;
	increment_val++;
	temperature_val++;
	}
	}

	StaticCommunicator::getStaticCommunicator().allReduce(norm, 2, _so_sum);

	norm[0] = sqrt(norm[0]);
	norm[1] = sqrt(norm[1]);

	AKANTU_DEBUG_ASSERT(!Math::isnan(norm[0]), "Something goes wrong in the solve phase");

	if (norm[1] < Math::getTolerance()) {
	error = norm[0];
	AKANTU_DEBUG_OUT();
	// cout<<"Error 1: "<<error<<endl;
	return error < tolerance;
	}


	AKANTU_DEBUG_OUT();
	if(norm[1] > Math::getTolerance())
	error = norm[0] / norm[1];
	else
	error = norm[0]; //In case the total displacement is zero!

	// cout<<"Error 2: "<<error<<endl;

	return (error < tolerance);
	}


	/* -------------------------------------------------------------------------- */
	void HeatTransferModel::initFEEngineBoundary(bool create_surface) {

	if(create_surface)
	MeshUtils::buildFacets(getFEEngine().getMesh());

	FEEngine & fem_boundary = getFEEngineBoundary();
	fem_boundary.initShapeFunctions();
	fem_boundary.computeNormalsOnIntegrationPoints();
	}

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

	Real HeatTransferModel::computeThermalEnergyByNode() {
	AKANTU_DEBUG_IN();

	Real ethermal = 0.;

	Array<Real>::vector_iterator heat_rate_it =
	residual->begin(residual->getNbComponent());

	Array<Real>::vector_iterator heat_rate_end =
	residual->end(residual->getNbComponent());


	UInt n = 0;
	for(;heat_rate_it != heat_rate_end; ++heat_rate_it, ++n) {
	Real heat = 0;
	bool is_local_node = mesh.isLocalOrMasterNode(n);
	bool is_not_pbc_slave_node = !isPBCSlaveNode(n);
	bool count_node = is_local_node && is_not_pbc_slave_node;

	Vector<Real> & heat_rate = *heat_rate_it;
	for (UInt i = 0; i < heat_rate.size(); ++i) {
	if (count_node)
	heat += heat_rate[i] * time_step;
	}
	ethermal += heat;
	}

	StaticCommunicator::getStaticCommunicator().allReduce(&ethermal, 1, _so_sum);

	AKANTU_DEBUG_OUT();
	return ethermal;
	}

	/* -------------------------------------------------------------------------- */
	template<class iterator>
	void HeatTransferModel::getThermalEnergy(iterator Eth,
	Array<Real>::const_iterator<Real> T_it,
	Array<Real>::const_iterator<Real> T_end) const {
	for(;T_it != T_end; ++T_it, ++Eth) {
	Eth = capacity density * *T_it;
	}
	}

	/* -------------------------------------------------------------------------- */
	Real HeatTransferModel::getThermalEnergy(const ElementType & type, UInt index) {
	AKANTU_DEBUG_IN();

	UInt nb_quadrature_points = getFEEngine().getNbIntegrationPoints(type);
	Vector<Real> Eth_on_quarature_points(nb_quadrature_points);

	Array<Real>::iterator<Real> T_it = this->temperature_on_qpoints(type).begin();
	T_it += index * nb_quadrature_points;

	Array<Real>::iterator<Real> T_end = T_it + nb_quadrature_points;

	getThermalEnergy(Eth_on_quarature_points.storage(), T_it, T_end);

	return getFEEngine().integrate(Eth_on_quarature_points, type, index);
	}

	/* -------------------------------------------------------------------------- */
	Real HeatTransferModel::getThermalEnergy() {
	Real Eth = 0;

	Mesh & mesh = getFEEngine().getMesh();
	Mesh::type_iterator it = mesh.firstType(spatial_dimension);
	Mesh::type_iterator last_type = mesh.lastType(spatial_dimension);
	for(; it != last_type; ++it) {
	UInt nb_element = getFEEngine().getMesh().getNbElement(*it, _not_ghost);
	UInt nb_quadrature_points = getFEEngine().getNbIntegrationPoints(*it, _not_ghost);
	Array<Real> Eth_per_quad(nb_element * nb_quadrature_points, 1);

	Array<Real>::iterator<Real> T_it = this->temperature_on_qpoints(*it).begin();
	Array<Real>::iterator<Real> T_end = this->temperature_on_qpoints(*it).end();
	getThermalEnergy(Eth_per_quad.begin(), T_it, T_end);

	Eth += getFEEngine().integrate(Eth_per_quad, *it);
	}

	return Eth;
	}

	/* -------------------------------------------------------------------------- */
	Real HeatTransferModel::getEnergy(const std::string & id) {
	AKANTU_DEBUG_IN();
	Real energy = 0;

	if("thermal") energy = getThermalEnergy();

	// reduction sum over all processors
	StaticCommunicator::getStaticCommunicator().allReduce(&energy, 1, _so_sum);

	AKANTU_DEBUG_OUT();
	return energy;
	}

	/* -------------------------------------------------------------------------- */
	Real HeatTransferModel::getEnergy(const std::string & energy_id, const ElementType & type, UInt index) {
	AKANTU_DEBUG_IN();

	Real energy = 0.;

	if("thermal") energy = getThermalEnergy(type, index);

	AKANTU_DEBUG_OUT();
	return energy;
	}

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

	/* -------------------------------------------------------------------------- */
	#ifdef AKANTU_USE_IOHELPER

	dumper::Field * HeatTransferModel::createNodalFieldBool(const std::string & field_name,
	const std::string & group_name,
	bool padding_flag) {

	std::map<std::string,Array<bool>* > uint_nodal_fields;
	uint_nodal_fields["blocked_dofs" ] = blocked_dofs;

	dumper::Field * field = NULL;
	field = mesh.createNodalField(uint_nodal_fields[field_name],group_name);
	return field;

	}
	/* -------------------------------------------------------------------------- */

	dumper::Field * HeatTransferModel::createNodalFieldReal(const std::string & field_name,
	const std::string & group_name,
	bool padding_flag) {

	std::map<std::string,Array<Real>* > real_nodal_fields;
	real_nodal_fields["temperature" ] = temperature;
	real_nodal_fields["temperature_rate" ] = temperature_rate;
	real_nodal_fields["external_heat_rate" ] = external_heat_rate;
	real_nodal_fields["residual" ] = residual;
	real_nodal_fields["capacity_lumped" ] = capacity_lumped;
	real_nodal_fields["increment" ] = increment;

	dumper::Field * field =
	mesh.createNodalField(real_nodal_fields[field_name],group_name);

	return field;
	}

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


	dumper::Field * HeatTransferModel
	::createElementalField(const std::string & field_name,
	const std::string & group_name,
	bool padding_flag,
	const UInt & spatial_dimension,
	const ElementKind & element_kind){


	dumper::Field * field = NULL;

	if(field_name == "partitions")
	field = mesh.createElementalField<UInt, dumper::ElementPartitionField>(mesh.getConnectivities(),group_name,this->spatial_dimension,element_kind);
	else if(field_name == "temperature_gradient"){
	ElementTypeMap<UInt> nb_data_per_elem = this->mesh.getNbDataPerElem(temperature_gradient,element_kind);

	field =
	mesh.createElementalField<Real,
	dumper::InternalMaterialField>(temperature_gradient,
	group_name,
	this->spatial_dimension,
	element_kind,
	nb_data_per_elem);
	}
	else if(field_name == "conductivity"){
	ElementTypeMap<UInt> nb_data_per_elem = this->mesh.getNbDataPerElem(conductivity_on_qpoints,element_kind);

	field =
	mesh.createElementalField<Real,
	dumper::InternalMaterialField>(conductivity_on_qpoints,
	group_name,
	this->spatial_dimension,
	element_kind,
	nb_data_per_elem);
	}

	return field;
	}
	/* -------------------------------------------------------------------------- */
	#else
	/* -------------------------------------------------------------------------- */
	dumper::Field * HeatTransferModel
	::createElementalField(const std::string & field_name,
	const std::string & group_name,
	bool padding_flag,
	const ElementKind & element_kind){

	return NULL;
	}

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

	dumper::Field * HeatTransferModel::createNodalFieldBool(const std::string & field_name,
	const std::string & group_name,
	bool padding_flag) {

	return NULL;

	}
	/* -------------------------------------------------------------------------- */

	dumper::Field * HeatTransferModel::createNodalFieldReal(const std::string & field_name,
	const std::string & group_name,
	bool padding_flag) {

	return NULL;
	}



	#endif


	__END_AKANTU__

heat_transfer_model.ccNo OneTemporaryActions

File Metadata

heat_transfer_model.ccView Options

Event Timeline

heat_transfer_model.cc
No OneTemporary
Actions

heat_transfer_model.cc
View Options