Homec4science
Diffusion akantu (bugfixes/all_bugs_found_by_working_on_tutorials) Browse

rAKA/src/model/heat_transfer/heat_transfer_model.cc
66aebf94f133bugfixes/all_bugs_found…

Last Change
Raw File

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 Srinivasa Babu Ramisetti <srinivasa.ramisetti@epfl.ch>
* @author David Simon Kammer <david.kammer@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 Sep 15 2014
*
* @brief Implementation of HeatTransferModel class
*
* @section LICENSE
*
* Copyright (©) 2014 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_material_internal_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) {
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();
PBCSynchronizer * synch = new PBCSynchronizer(pbc_pair);
synch_registry->registerSynchronizer(*synch, _gst_htm_capacity);
synch_registry->registerSynchronizer(*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().getNbQuadraturePoints(*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 (solver) delete solver;
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().getNbQuadraturePoints(*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().getNbQuadraturePoints(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^t*D*B");
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().interpolateOnQuadraturePoints(*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().gradientOnQuadraturePoints(*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.getNbQuadraturePoints(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.computeNormalsOnControlPoints();
}
/* -------------------------------------------------------------------------- */
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().getNbQuadraturePoints(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().getNbQuadraturePoints(*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__
c4science · Help