Page MenuHomec4science
Files F75875898

compute_ekin_point.cc
No OneTemporary
Actions

File Metadata

Created
Sun, Aug 4, 20:29

compute_ekin_point.cc
View Options

/**
* @file compute_ekin_point.cc
*
* @author Guillaume Anciaux <guillaume.anciaux@epfl.ch>
*
* @date Thu Sep 18 16:13:19 2014
*
* @brief This computes the kinetic energy of a set of DOFs
*
* @section LICENSE
*
* Copyright (©) 2010-2011 EPFL (Ecole Polytechnique Fédérale de Lausanne)
* Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)
*
* LibMultiScale 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.
*
* LibMultiScale 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 LibMultiScale. If not, see <http://www.gnu.org/licenses/>.
*
*/
/* -------------------------------------------------------------------------- */
#include "lm_common.hh"
#include "lib_md.hh"
#include "lib_dd.hh"
#include "lib_continuum.hh"
#include "compute_ekin.hh"
#include "ref_point_data.hh"
#include "communicator.hh"
#include "units_converter.hh"
#include "units.hh"
/* -------------------------------------------------------------------------- */
__BEGIN_LIBMULTISCALE__
template <typename Cont>
void ComputeEKin<Cont>::build(Cont & cont){
static const UInt Dim = Cont::Dim;
Real computed_ekin[3] = {0,0,0};
UInt cpt = cont.nbElem();
typename Cont::iterator it = cont.getIterator();
typedef typename Cont::Ref RefPoint;
for (RefPoint at = it.getFirst();!it.end();at = it.getNext()) {
Real vel[Dim];
at.getVelocities(vel);
for (UInt i = 0; i < Dim; ++i) computed_ekin[i] += vel[i]*vel[i]*at.mass();
}
for (UInt i = 0; i < Dim; ++i){
computed_ekin[i] *= 0.5*code_unit_system.mvv2e;
DUMP("computed ekin[" << i << "] = " << computed_ekin[i],DBG_INFO);
}
Communicator & comm = Communicator::getCommunicator();
CommGroup group = cont.getCommGroup();
comm.allReduce(computed_ekin,Dim,group,"kinetic energy reduce", OP_SUM);
comm.allReduce(&cpt,1,group,"kinetic energy counter reduce", OP_SUM);
Real ekin = 0.0;
for (UInt i = 0; i < Dim; ++i){
DUMP("ekin[" << i << "] = " << computed_ekin[i],DBG_INFO);
ekin += computed_ekin[i];
}
//compute temperature
UnitsConverter unit;
unit.setInUnits(code_unit_system);
unit.setOutUnits(real_unit_system);
unit.computeConversions();
Real T = 2./cpt/Dim/boltzmann*(ekin*unit.getConversion<Energy>());
//convert the kinetic energies to the required unit systems
// unit.setOutUnits(this->unit_system);
// unit.computeConversions();
// ekin *= unit.getConversion<Energy>();
// for (UInt i = 0; i < Dim; ++i) computed_ekin[i] *= unit.getConversion<Energy>();
DUMP("nb_points = " << cpt,DBG_INFO);
this->empty();
Real nb_points = static_cast<Real>(cpt);
this->add(nb_points);
this->add(T);
this->add(ekin);
for (UInt i = 0; i < Dim; ++i) this->add(computed_ekin[i]);
this->name_computed.clear();
this->name_computed.push_back("nbpoints");
this->name_computed.push_back("Temp");
this->name_computed.push_back("EKin");
this->name_computed.push_back("EKin_x");
if (Dim > 1) this->name_computed.push_back("EKin_y");
if (Dim == 3) this->name_computed.push_back("EKin_z");
this->setDim(Dim + 3);
}
/* -------------------------------------------------------------------------- */
DECLARE_COMPUTE_REF(ComputeEKin,LIST_ATOM_MODEL);
DECLARE_COMPUTE_REF(ComputeEKin,LIST_DD_MODEL);
DECLARE_COMPUTE_REFPOINT(ComputeEKin);
DECLARE_COMPUTE_GENERIC_MESH(ComputeEKin);
__END_LIBMULTISCALE__

Event Timeline

c4science · Help