diff --git a/src/common/lm_globals.cc b/src/common/lm_globals.cc index 6150a29..302c205 100644 --- a/src/common/lm_globals.cc +++ b/src/common/lm_globals.cc @@ -1,124 +1,124 @@ /** * @file lm_globals.cc * * @author Guillaume Anciaux <guillaume.anciaux@epfl.ch> * @author Jaehyun Cho <jaehyun.cho@epfl.ch> * @author Till Junge <till.junge@epfl.ch> * * @date Mon Sep 08 23:40:22 2014 * * @brief This file contains all the global variables of LM * * @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 <fstream> #include <limits> /* -------------------------------------------------------------------------- */ __BEGIN_LIBMULTISCALE__ /* -------------------------------------------------------------------------- */ // global variables UInt current_step; IntegrationSchemeStage current_stage; Real current_time; UInt nb_step; UInt nb_step_next_event; UInt global_level; UInt global_mod; UInt global_proc; UInt global_proc1; UInt global_proc2; UInt start_dump; UInt end_dump; std::string my_hostname; UInt print_crap_to_file; std::ostream *global_out; std::ofstream file_out; UInt wait_on_fatal; UInt wait_at_startup; UInt create_seg_fault; UInt lm_my_proc_id; UInt lm_world_size; UInt lm_job_id; UInt spatial_dimension; UInt lm_uint_max = std::numeric_limits<UInt>::max(); Real lm_real_max = std::numeric_limits<Real>::max(); /// used to avoid infinite FATAL recursion in stage PRE_FATAL bool fatal_loop = false; #if defined(LIBMULTISCALE_USE_TIMER) std::map<std::string, struct timespec> mesT; std::map<std::string, UInt> nmes; std::map<std::string, struct timespec> tstart; std::map<std::string, struct timespec> tstop; #endif /* -------------------------------------------------------------------------- */ UInt getGlobalCurrentRelease() { UInt k = 0; switch (current_stage) { case PRE_DUMP: k = 0; break; case PRE_STEP1: k = 1; break; case PRE_STEP2: k = 2; break; case PRE_STEP3: k = 3; break; case PRE_STEP4: k = 4; break; case PRE_STEP5: k = 5; break; case PRE_STEP6: k = 6; break; case PRE_STEP7: k = 7; break; case PRE_FATAL: k = 8; break; default: LM_FATAL("error: unknown stage " << current_stage); } - return (8 * current_step + k + INITIAL_MODEL_RELEASE); + return (8 * current_step + k); } /* -------------------------------------------------------------------------- */ const int CommGroup::id_all = std::numeric_limits<int>::max(); const int CommGroup::id_none = id_all - 1; const int CommGroup::id_invalid = id_all - 2; CommGroup group_all(CommGroup::id_all); CommGroup group_none(CommGroup::id_none); CommGroup group_invalid(CommGroup::id_invalid); /* -------------------------------------------------------------------------- */ __END_LIBMULTISCALE__ diff --git a/src/common/lm_globals.hh b/src/common/lm_globals.hh index 621cc49..c0901e7 100644 --- a/src/common/lm_globals.hh +++ b/src/common/lm_globals.hh @@ -1,97 +1,95 @@ /** * @file lm_globals.hh * * @author Guillaume Anciaux <guillaume.anciaux@epfl.ch> * @author Jaehyun Cho <jaehyun.cho@epfl.ch> * @author Till Junge <till.junge@epfl.ch> * * @date Mon Sep 08 23:40:22 2014 * * @brief This file contains all the global variables of LM * * @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/>. * */ #ifndef __LIBMULTISCALE_LM_GLOBALS_HH__ #define __LIBMULTISCALE_LM_GLOBALS_HH__ /* -------------------------------------------------------------------------- */ #include "lm_macros.hh" #include "lm_types.hh" #include <iostream> /* -------------------------------------------------------------------------- */ __BEGIN_LIBMULTISCALE__ static const std::string invalidFilter = "invalid filter"; static const std::string invalidGeom = "invalid geometry"; static const std::string invalidDomain = "invalid domain"; static const std::string invalidDumper = "invalid dumper"; static const std::string invalidStimulator = "invalid stimulator"; typedef std::string GeomID; typedef std::string FilterID; typedef std::string DomainID; typedef std::string ActionID; typedef std::string DumperID; typedef std::string StimulatorID; /* -------------------------------------------------------------------------- */ // global variables EXTERN UInt current_step; EXTERN IntegrationSchemeStage current_stage; EXTERN Real current_time; EXTERN UInt nb_step; EXTERN UInt nb_step_next_event; EXTERN UInt global_level; EXTERN UInt global_mod; EXTERN UInt global_proc; EXTERN UInt global_proc1; EXTERN UInt global_proc2; EXTERN UInt start_dump; EXTERN UInt end_dump; EXTERN std::string my_hostname; EXTERN UInt print_crap_to_file; EXTERN std::ostream *global_out; EXTERN std::ofstream file_out; EXTERN UInt wait_on_fatal; EXTERN UInt wait_at_startup; EXTERN UInt create_seg_fault; EXTERN UInt lm_my_proc_id; EXTERN UInt lm_world_size; EXTERN UInt lm_job_id; EXTERN UInt spatial_dimension; EXTERN UInt lm_uint_max; EXTERN Real lm_real_max; /// used to avoid infinite FATAL recursion in stage PRE_FATAL EXTERN bool fatal_loop; /* -------------------------------------------------------------------------- */ EXTERN UInt getGlobalCurrentRelease(); /* -------------------------------------------------------------------------- */ EXTERN CommGroup group_all; EXTERN CommGroup group_none; EXTERN CommGroup group_invalid; /* -------------------------------------------------------------------------- */ -static const UInt INITIAL_MODEL_RELEASE = 100; -/* -------------------------------------------------------------------------- */ __END_LIBMULTISCALE__ #endif /* __LIBMULTISCALE_LM_GLOBALS_HH__ */ diff --git a/src/md/1d/domain_md1d.cc b/src/md/1d/domain_md1d.cc index fa0afc0..358b11a 100644 --- a/src/md/1d/domain_md1d.cc +++ b/src/md/1d/domain_md1d.cc @@ -1,534 +1,532 @@ /** * @file domain_md1d.cc * * @author Guillaume Anciaux <guillaume.anciaux@epfl.ch> * @author Srinivasa Babu Ramisetti <srinivasa.ramisetti@epfl.ch> * * @date Fri Jan 10 20:47:45 2014 * * @brief This is the domain for 1d atomic chains * * @section LICENSE * * Copyright INRIA and CEA * * The LibMultiScale is a C++ parallel framework for the multiscale * coupling methods dedicated to material simulations. This framework * provides an API which makes it possible to program coupled simulations * and integration of already existing codes. * * This Project was initiated in a collaboration between INRIA Futurs Bordeaux * within ScAlApplix team and CEA/DPTA Ile de France. * The project is now continued at the Ecole Polytechnique Fédérale de Lausanne * within the LSMS/ENAC laboratory. * * This software is governed by the CeCILL-C license under French law and * abiding by the rules of distribution of free software. You can use, * modify and/ or redistribute the software under the terms of the CeCILL-C * license as circulated by CEA, CNRS and INRIA at the following URL * "http://www.cecill.info". * * As a counterpart to the access to the source code and rights to copy, * modify and redistribute granted by the license, users are provided only * with a limited warranty and the software's author, the holder of the * economic rights, and the successive licensors have only limited * liability. * * In this respect, the user's attention is drawn to the risks associated * with loading, using, modifying and/or developing or reproducing the * software by the user in light of its specific status of free software, * that may mean that it is complicated to manipulate, and that also * therefore means that it is reserved for developers and experienced * professionals having in-depth computer knowledge. Users are therefore * encouraged to load and test the software's suitability as regards their * requirements in conditions enabling the security of their systems and/or * data to be ensured and, more generally, to use and operate it in the * same conditions as regards security. * * The fact that you are presently reading this means that you have had * knowledge of the CeCILL-C license and that you accept its terms. * */ /* -------------------------------------------------------------------------- */ #include "domain_md1d.hh" #include "lm_common.hh" #include "md_builder.hh" #include "reference_manager.hh" #include "reference_manager_md1d.hh" #include "spatial_filter_atomic.hh" /* -------------------------------------------------------------------------- */ __BEGIN_LIBMULTISCALE__ /* -------------------------------------------------------------------------- */ void DomainMD1D::correctPBCPositions() { int nb_atoms = this->getNbDofs(); Real domainSize = this->getDomainSize(); for (int i = 0; i < nb_atoms; ++i) { dofs[i].p -= domainSize * nearbyint(dofs[i].p / domainSize); } } /* -------------------------------------------------------------------------- */ inline Real DomainMD1D::computeDistance(const Real &xi, const Real &xj) { Real rij = xj - xi; Real domainSize = this->getDomainSize(); if (this->periodic_flag) rij -= domainSize * nearbyint(rij / domainSize); return rij; } /* -------------------------------------------------------------------------- */ Real DomainMD1D::getDomainSize() { if (this->periodic_flag) return this->getNbDofs() * r0; return (this->getNbDofs() - 1) * r0; } /* -------------------------------------------------------------------------- */ void DomainMD1D::init() { if (!Communicator::getCommunicator().amIinGroup(this->getGroupID())) return; ZoneMDBuilder z; Cube b = this->getGeom().getBoundingBox(); if (this->atom_file == "") z.buildMD1D(b, this->r0, this->mass, this->shift, this->periodic_flag, this->dofs); else z.readMD1D(this->atom_file, b, this->r0, this->mass, this->dofs); DUMP("Size of domain " << this->getDomainSize(), DBG_MESSAGE); DUMP("position offset " << -1.0 * this->getDomainSize() / 2 << " Center " << b.getCenter(0) << " Shift " << shift, DBG_MESSAGE); static_cast<ContainerMD1D &>(atoms).setDofs(dofs); DUMP(" MD1D R0 " << r0 << " RCUT " << rcut << " Rcut/r0 " << (Real)(rcut / r0), DBG_MESSAGE); int nb_atoms = this->getNbDofs(); stress = new Real[nb_atoms]; stiffness_parameter = 36. * epsilon * pow(2., 2. / 3.) / sigma / sigma; DUMP("Stiffness parameter " << stiffness_parameter << " epsilon " << epsilon << " sigma " << sigma, DBG_MESSAGE); this->computeNeighborLists(); - LM_TOIMPLEMENT; - //this->atoms.setRelease(INITIAL_MODEL_RELEASE); } /* -------------------------------------------------------------------------- */ Real DomainMD1D::getEcin() { LM_TOIMPLEMENT; } /* -------------------------------------------------------------------------- */ Real DomainMD1D::getEpot() { LM_TOIMPLEMENT; } /* -------------------------------------------------------------------------- */ DomainMD1D::DomainMD1D(DomainID ID, CommGroup GID) : DomainAtomic<ContainerMD1D, 1>(ID, GID), periodic_flag(false), elastic_flag(false), shift(0.0) {} /* -------------------------------------------------------------------------- */ void DomainMD1D::performStep1() { DUMP("Position step avec dt = " << timeStep, DBG_INFO); int nb_atoms = this->getNbDofs(); for (int i = 0; i < nb_atoms; ++i) { /* mise ajour des positions */ dofs[i].v += code_unit_system.ft_m2v * .5 * timeStep * dofs[i].f / dofs[i].m; dofs[i].p += timeStep * dofs[i].v; LM_ASSERT(!isinf(dofs[i].p), "one atomic position " << i << " is attributed an infinite position : " << dofs[i].v << " " << dofs[i].a << " " << timeStep); DUMP(getID() << " dof " << i << " x0 = " << dofs[i].p0 << " u = " << dofs[i].p - dofs[i].p0 << " v = " << dofs[i].v << " a = " << dofs[i].a, DBG_DETAIL); } if (periodic_flag) this->correctPBCPositions(); LM_TOIMPLEMENT; //this->atoms.incRelease(); } /* -------------------------------------------------------------------------- */ void DomainMD1D::performStep3() { DUMP("Velocity step avec dt = " << timeStep, DBG_INFO); int nb_atoms = dofs.size(); for (int i = 0; i < nb_atoms; ++i) { dofs[i].v += code_unit_system.ft_m2v * 0.5 * timeStep * dofs[i].f / dofs[i].m; DUMP(getID() << " dof " << i << " v = " << dofs[i].v << " a = " << dofs[i].a, DBG_DETAIL); } LM_TOIMPLEMENT; //this->atoms.incRelease(); } /* -------------------------------------------------------------------------- */ inline void DomainMD1D::computeNeighborLists() { UInt nb_atoms = dofs.size(); neighbor_lists.resize(nb_atoms); for (UInt i = 0; i < nb_atoms; ++i) { std::vector<UInt> neighbors; for (UInt j = i + 1; j < nb_atoms; ++j) { Real rij = std::fabs(this->computeDistance(dofs[j].p, dofs[i].p)); LM_ASSERT(rij / sigma > 1e-6, "problem in neighbor list: " << i << " " << j); if (rij <= rcut * 1.5) neighbors.push_back(j); } neighbor_lists[i] = neighbors; // for (auto j: neighbor_lists[i]) // DUMP("neighborlist " << i << ":" << j, DBG_MESSAGE); } } /* -------------------------------------------------------------------------- */ inline void DomainMD1D::computeContribution(UInt i, UInt j) { Real rij = computeDistance(dofs[i].p, dofs[j].p); Real sign = rij < 0 ? -1. : 1.; rij = std::fabs(rij); LM_ASSERT(std::fabs(rij) / sigma > 1e-6, "problem in neighbor list"); Real contribution = 24.0 * epsilon * (2.0 * MathTools::ipow(sigma / rij, 12) / rij - 1.0 * MathTools::ipow(sigma / rij, 6) / rij); if (arlequin_flag) contribution *= 0.5 * (dofs[i].alpha + dofs[j].alpha); DUMP("contribution between " << i << " and " << j << " = " << contribution << " [d=" << rij << "]", DBG_ALL); dofs[i].a -= contribution * sign; dofs[j].a += contribution * sign; } /* -------------------------------------------------------------------------- */ void DomainMD1D::performStep2() { // if (elastic_flag) { // quadraticForceComputation(); // return; // } LM_ASSERT(dofs[3].m != 0, "mass pb"); UInt nb_atoms = dofs.size(); for (UInt i = 0; i < nb_atoms; ++i) { dofs[i].f = 0.; dofs[i].a = 0.; } // loop over atoms to compute their force for (UInt i = 0; i < nb_atoms; ++i) { // loop over right neighbors for (auto &&j : neighbor_lists[i]) { Real rij = std::fabs(this->computeDistance(dofs[j].p, dofs[i].p)); if (rij <= rcut) computeContribution(i, j); } } for (UInt i = 0; i < nb_atoms; ++i) { dofs[i].f = dofs[i].a * dofs[i].m; } for (UInt i = 0; i < nb_atoms; ++i) { // if (arlequin_flag) { // dofs[i].a /= dofs[i].m * dofs[i].alpha; // dofs[i].f = dofs[i].a * dofs[i].m; // } else { // dofs[i].a /= dofs[i].m; // dofs[i].f = dofs[i].a * dofs[i].m; // } LM_ASSERT(!isinf(dofs[i].a), "one atomic " << i << " position is attributed an infinite acceleration : " << dofs[i].m); } LM_TOIMPLEMENT; //this->atoms.incRelease(); } /* -------------------------------------------------------------------------- */ // Added by Srinivasa Babu Ramisetti void DomainMD1D::quadraticForceComputation() { double contribution; double rij; double pos_right = 0, pos_left = 0; int i; int nb_atoms = dofs.size(); for (i = 0; i < nb_atoms; ++i) { dofs[i].pe = 0.0; double pos_i = dofs[i].p; int left = (i - 1 + nb_atoms) % nb_atoms; int right = (i + 1) % nb_atoms; if (periodic_flag != 1) { left = i - 1; right = i + 1; } if (right < nb_atoms) { pos_right = dofs[right].p; rij = computeDistance(pos_i, pos_right); contribution = (rij - r0) * stiffness_parameter; dofs[i].f += contribution; dofs[i].pe += (rij - r0) * contribution / 4.0; } if (left >= 0) { pos_left = dofs[left].p; rij = computeDistance(pos_i, pos_left); contribution = (rij - r0) * stiffness_parameter; dofs[i].f -= contribution; dofs[i].pe += (rij - r0) * contribution / 4.0; } dofs[i].a = dofs[i].f / dofs[i].m; } } /* -------------------------------------------------------------------------- */ // Added by Srinivasa Babu Ramisetti // void DomainMD1D::getNeighborList(RefAtom & atom, Real distance, // std::vector<RefAtom> & neigbhor_list) void DomainMD1D::getNeighborList(RefAtom &atom, UInt count, std::vector<RefAtom> &neigbhor_list) { UInt index_count = 0; neigbhor_list.clear(); RefMD1D r1d(dofs); int current_id = atom.id(); int nb_atoms = dofs.size(); if (periodic_flag == 1) { // DUMP(" Current atom id "<< atom.id() << " pos0 " << atom.position0(0) << // " nb_coefficients " << count << " nb_atoms " << nb_atoms, DBG_MESSAGE); int id = 1; int i_id1 = (current_id + id); int i_id2 = (current_id - id); // Real d = fabs(atom.position0(0) - dofs[i_id1].p0); // if(current_id == (nb_atoms - 1)) // d = fabs(atom.position0(0) - dofs[i_id2].p0); // while(d < distance){ while (index_count < count) { i_id1 = (current_id + id); if (i_id1 < 0 || i_id1 >= nb_atoms) { i_id1 = (current_id + id) % nb_atoms; // d = fabs(atom.position0(0) - dofs[current_id - id].p0); } r1d.setIndex(i_id1); neigbhor_list.push_back(r1d); i_id2 = (current_id - id); if (i_id2 < 0 || i_id2 > nb_atoms) { i_id2 = (current_id - id + nb_atoms) % nb_atoms; // d = fabs(atom.position0(0) - dofs[current_id + id].p0); } r1d.setIndex(i_id2); neigbhor_list.push_back(r1d); id++; index_count++; // DUMP(" l_i "<< i_id2 << " r_i " << i_id1 << " vel[l] " << v2 << " // vel[r] " << v1 , DBG_MESSAGE); } } else { DUMP(" Current atom id " << atom.id() << " pos0 " << atom.position0() << " nb_atoms " << nb_atoms, DBG_DETAIL); int id = 1; // Real d = fabs(atom.position0(0) - dofs[current_id + 1].p0); // while(d < distance){ while (index_count < count) { if ((current_id + id) < 0 || (current_id + id) > nb_atoms) LM_FATAL("Neighbor id's are out of bound!"); r1d.setIndex(current_id + id); neigbhor_list.push_back(r1d); if ((current_id - id) < 0 || (current_id - id) > nb_atoms) LM_FATAL("Neighbor id's are out of bound!"); r1d.setIndex(current_id - id); neigbhor_list.push_back(r1d); DUMP(" neigbhor id on left : " << current_id - id << " neigbhor id on right : " << current_id + id, DBG_DETAIL); id++; // d = fabs(atom.position0(0) - dofs[current_id + id].p0); // count = count + 2; index_count++; } } DUMP("Number of neighbor's found " << index_count * 2, DBG_DETAIL); } Real DomainMD1D::getElasticConstant() { return this->stiffness_parameter; } /* -------------------------------------------------------------------------- */ Real DomainMD1D::getR0() { return this->r0; } /* -------------------------------------------------------------------------- */ int DomainMD1D::getNbDofs() { return this->dofs.size(); } /* -------------------------------------------------------------------------- */ UInt DomainMD1D::getDim() { return 1; } /* -------------------------------------------------------------------------- */ void DomainMD1D::getDomainDimensions(Real *xmin, Real *xmax) { xmin[0] = lm_real_max; xmax[0] = -1. * lm_real_max; int nb_atoms = dofs.size(); for (int i = 0; i < nb_atoms; ++i) { if (xmin[0] > dofs[i].p0) xmin[0] = dofs[i].p0; if (xmax[0] < dofs[i].p0) xmax[0] = dofs[i].p0; } } /* -------------------------------------------------------------------------- */ void DomainMD1D::getSubDomainDimensions(Real *xmin, Real *xmax) { LM_TOIMPLEMENT; } /* -------------------------------------------------------------------------- */ void DomainMD1D::setCurrentStep(UInt ts) { LM_TOIMPLEMENT; } /* -------------------------------------------------------------------------- */ UInt DomainMD1D::getCurrentStep() { LM_TOIMPLEMENT; return 0; } /* -------------------------------------------------------------------------- */ void DomainMD1D::setTimeStep(Real ts) { LM_TOIMPLEMENT; } /* -------------------------------------------------------------------------- */ Real DomainMD1D::getTimeStep() { return timeStep; } /* -------------------------------------------------------------------------- */ void DomainMD1D::setPeriodicFlag(bool flag) { periodic_flag = flag; } /* -------------------------------------------------------------------------- */ void DomainMD1D::setReferenceManagerState(bool state) {} /* -------------------------------------------------------------------------- */ bool DomainMD1D::isPeriodic(UInt dir) { return (periodic_flag && dir == 1); } /* -------------------------------------------------------------------------- */ /* LMDESC MD1D This plugin is the implementation of a very simple one dimensional molecular dynamics class. */ /* LMHERITANCE domain_md */ /* LMEXAMPLE Section MultiScale AtomsUnits \\ ...\\ GEOMETRY myGeom CUBE BBOX -1 1 -1 1 -1 1\\ MODEL MD1D md\\ ...\\ endSection\\ \\ Section MD1D:md AtomsUnits \\ RCUT rcut \\ R0 r0 \\ MASS 39.95e-3 \\ TIMESTEP 1 \\ EPSILON 1.6567944e-21/6.94769e-21 \\ SIGMA 1.1 \\ \\ DOMAIN_GEOMETRY 1 \\ endSection \\ */ void DomainMD1D::declareParams() { DomainAtomic<ContainerMD1D, 1>::declareParams(); /* LMKEYWORD RCUT Sets the cutoff radius for force calculation */ this->parseKeyword("RCUT", rcut); /* LMKEYWORD R0 Sets the interatomic distance */ this->parseKeyword("R0", r0); /* LMKEYWORD SHIFT used in geometry creation */ this->parseKeyword("SHIFT", shift, 0.); /* LMKEYWORD MASS Sets the per atom mass */ this->parseKeyword("MASS", mass); /* LMKEYWORD TIMESTEP Sets the timestep */ this->parseKeyword("TIMESTEP", timeStep, 1.); /* LMKEYWORD EPSILON Sets the LJ energy parameter */ this->parseKeyword("EPSILON", epsilon); /* LMKEYWORD SIGMA Sets the LJ distance parameter */ this->parseKeyword("SIGMA", sigma); /* LMKEYWORD PERIODIC Sets the periodic boundary condition */ this->parseTag("PERIODIC", periodic_flag, false); /* LMKEYWORD ELASTIC Sets the flag to switch elastic interactions */ this->parseTag("ELASTIC", elastic_flag, false); /* LMKEYWORD ATOM_FILE Sets the file to read the atomic positions */ this->parseKeyword("ATOM_FILE", atom_file, ""); } /* -------------------------------------------------------------------------- */ __END_LIBMULTISCALE__