diff --git a/src/common/attached_object.hh b/src/common/attached_object.hh index 7a9c300..479b411 100644 --- a/src/common/attached_object.hh +++ b/src/common/attached_object.hh @@ -1,104 +1,121 @@ /** * @file attached_object.hh * * @author Guillaume Anciaux * * @date Mon Sep 08 23:40:22 2014 * * @brief This object is used to describe objects attached to containers of * references * * @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 . * */ #ifndef __LIBMULTISCALE_ATTCHED_OBJECT_HH__ #define __LIBMULTISCALE_ATTCHED_OBJECT_HH__ /* -------------------------------------------------------------------------- */ #include "auto_argument.hh" #include "lm_common.hh" #include "pack_buffer.hh" /* -------------------------------------------------------------------------- */ __BEGIN_LIBMULTISCALE__ /* -------------------------------------------------------------------------- */ template class ContainerArray; class AttachedObject : public AutoDispatch::ArgumentAny { public: virtual ~AttachedObject(){}; - virtual void packData(UInt index, PackBuffer &buffer, UInt toproc) = 0; - virtual void unpackData(UInt index, PackBuffer &buffer, UInt fromproc) = 0; + virtual void packData(UInt index, PackBuffer &buffer, UInt toproc, + bool verbose = false) = 0; + virtual void unpackData(UInt index, PackBuffer &buffer, UInt fromproc, + bool verbose = false) = 0; virtual void moveAttachedValues(UInt i_src, UInt i_dest) = 0; virtual void resize(UInt) { LM_TOIMPLEMENT; } + virtual LMID getID() const = 0; }; /* -------------------------------------------------------------------------- */ template class AttachedVector : public AttachedObject { public: AttachedVector(ContainerArray &vector) { *this = vector; }; virtual ~AttachedVector() = default; using AutoDispatch::ArgumentAny::operator=; - void packData(UInt index, PackBuffer &buffer, UInt) { + void packData(UInt index, PackBuffer &buffer, UInt toproc, + bool verbose = false) override { auto &v = this->getV(); for (UInt i = 0; i < v.cols(); ++i) { + if (verbose) { + DUMP("packing to proc " << toproc << ": " << v(index, i), DBG_MESSAGE); + } buffer << v(index, i); } } - void unpackData(UInt index, PackBuffer &buffer, UInt) { + void unpackData(UInt index, PackBuffer &buffer, UInt fromproc, + bool verbose = false) override { auto &v = this->getV(); if (index >= v.rows()) v.resize(index + 1); for (UInt i = 0; i < v.cols(); ++i) { buffer >> v(index, i); } + if (verbose) { + DUMP("unpacking from proc " << fromproc << ": " << v.row(index), + DBG_MESSAGE); + } } void moveAttachedValues(UInt i_src, UInt i_dest) { auto &v = this->getV(); v.row(i_dest) = v.row(i_src); } void resize(UInt sz) { auto &v = this->getV(); DUMP(&v << " resizes from " << v.rows() << " to " << sz, DBG_DETAIL); v.resize(sz); } + LMID getID() const override { return this->getV().getID(); }; + ContainerArray &getV() { return this->get>(); }; + const ContainerArray &getV() const { + return this->get>(); + }; }; /* -------------------------------------------------------------------------- */ template std::shared_ptr make_attached(ContainerArray &obj) { return std::make_shared>(obj); } /* -------------------------------------------------------------------------- */ __END_LIBMULTISCALE__ #endif /* __LIBMULTISCALE_ATTCHED_OBJECT_HH__ */ diff --git a/src/common/id_manager.hh b/src/common/id_manager.hh index 90b5a66..2d37919 100644 --- a/src/common/id_manager.hh +++ b/src/common/id_manager.hh @@ -1,299 +1,299 @@ /** * @file id_manager.hh * * @author Guillaume Anciaux * * @date Mon Sep 08 23:40:22 2014 * * @brief This is the generic manager of objects from their IDs (string) * * @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 . * */ #ifndef __LIBMULTISCALE_ID_MANAGER_HH__ #define __LIBMULTISCALE_ID_MANAGER_HH__ /* -------------------------------------------------------------------------- */ #include "lm_common.hh" #include "lm_object.hh" #include #include #include #include #include /* -------------------------------------------------------------------------- */ __BEGIN_LIBMULTISCALE__ struct UnknownID : public std::runtime_error { UnknownID(const LMID &id, const std::string &manager_type) : std::runtime_error("Unknown id: " + id + " from manager type: " + manager_type), id(id) {} LMID id; }; template class IDManager { /* ------------------------------------------------------------------------ */ /* Typedefs */ /* ------------------------------------------------------------------------ */ public: typedef std::string typeID; typedef std::map Cont; /* ------------------------------------------------------------------------ */ /* Constructors/Destructors */ /* ------------------------------------------------------------------------ */ friend struct std::default_delete>; protected: IDManager(); virtual ~IDManager(); /* ------------------------------------------------------------------------ */ /* Methods */ /* ------------------------------------------------------------------------ */ public: static void destroy(); //! add a geometry as a shared pointer typeID addObject(std::shared_ptr g); //! add a geometry as an object typeID addObject(T &g); //! generate a new id automatically typeID generateNewID(); //! remove a managed object void removeObject(const typeID &ID); //! remove a managed object void removeObject(const T &ID); //! true if singleton already allocated static bool allocated() { return static_pointer != nullptr; } //! applies a functor for each stored entry template void for_each(Func &&func) { std::for_each(this->objects.begin(), this->objects.end(), [&](auto &&pair) { func(*pair.second); }); } /* ------------------------------------------------------------------------ */ /* Accessors */ /* ------------------------------------------------------------------------ */ public: //! return a reference to a geometry from its ID T &getObject(const typeID &ID); //! return a shared_pointer from its ID std::shared_ptr getSharedObject(const typeID &ID); //! return a casted object from its ID template inline Cast &getCastedObject(const typeID &ID); //! search id of a given geometry typeID getID(T &obj); protected: //! accessor to the singleton manager template static IT &getManager(); /* ------------------------------------------------------------------------ */ /* Class Members */ /* ------------------------------------------------------------------------ */ protected: //! mapping between IDs and pointers std::map objects; //! mapping between IDs and shared pointers std::map> shared_pointer_objects; //! file to which output the object infos std::ofstream fout; static std::unique_ptr> static_pointer; std::list insert_order; }; /* -------------------------------------------------------------------------- */ template template inline IT &IDManager::getManager() { if (!static_pointer) static_pointer.reset(new IT()); IT &ref = dynamic_cast(*static_pointer); return ref; } /* -------------------------------------------------------------------------- */ template inline IDManager::IDManager() { std::string filename = "log."; filename += typeinfo(); fout.open(filename.c_str(), std::ios_base::out); } /* -------------------------------------------------------------------------- */ template inline typename IDManager::typeID IDManager::generateNewID() { UInt i = objects.size(); std::stringstream str; str << "geom" << i; return str.str(); } /* -------------------------------------------------------------------------- */ template inline typename IDManager::typeID IDManager::addObject(std::shared_ptr obj) { if (obj == nullptr) LM_FATAL("attempt to register a null pointer as " << typeid(*this).name()); typeID id = obj->getID(); if (objects.count(id) > 0 && objects[id]) LM_FATAL("ID " << id << " has already been used for " << *objects[id]); DUMP("register object: " << obj.get(), DBG_INFO); shared_pointer_objects[id] = obj; objects[id] = obj.get(); DUMPFILE(fout, id << "(" << *obj << ") added to " << typeinfo() << " manager"); insert_order.push_back(id); return id; } /* -------------------------------------------------------------------------- */ template inline typename IDManager::typeID IDManager::addObject(T &obj) { typeID id = obj.getID(); if (objects.count(id) > 0 && objects[id]) LM_FATAL("ID " << id << " has already been used for " << *objects[id]); objects[id] = &obj; DUMPFILE(fout, id << "(" << obj << ") added to " << typeinfo() << " manager"); insert_order.push_back(id); return id; } /* -------------------------------------------------------------------------- */ template inline void IDManager::removeObject(const T &obj) { typeID ID = obj.getID(); removeObject(ID); DUMPFILE(fout, "removed " << ID << " from manager " << typeinfo()); } /* -------------------------------------------------------------------------- */ template inline void IDManager::removeObject(const typeID &ID) { if (objects.count(ID) == 0) { DUMP("impossible to delete entry -> unknown id: " << ID, DBG_INFO); return; } objects.erase(ID); if (shared_pointer_objects.count(ID)) shared_pointer_objects.erase(ID); insert_order.remove(ID); } /* -------------------------------------------------------------------------- */ template inline T &IDManager::getObject(const typeID &ID) { if (objects.count(ID) == 0) { throw UnknownID(ID, typeinfo()); } return *objects[ID]; } /* -------------------------------------------------------------------------- */ template inline std::shared_ptr IDManager::getSharedObject(const typeID &ID) { if (shared_pointer_objects.count(ID) == 0) { throw UnknownID(ID, typeinfo()); } return shared_pointer_objects[ID]; } /* -------------------------------------------------------------------------- */ template template inline Cast &IDManager::getCastedObject(const typeID &ID) { std::string component_name = ID; - std::string output_name = ID; + std::string output_name = ""; auto pos = component_name.find("."); if (pos != std::string::npos) { output_name = component_name.substr(pos + 1); component_name = component_name.substr(0, pos); } if (objects.count(component_name) == 0) { throw UnknownID(ID, typeinfo()); } Cast *ptr = nullptr; if (output_name != "") ptr = &objects[component_name]->template evalOutput(output_name); else ptr = dynamic_cast(objects[component_name]); if (!ptr) { LM_FATAL("bad cast for ID: " << ID << " from " << typeinfo() << " to " << typeinfo()); } return *ptr; } /* -------------------------------------------------------------------------- */ template inline typename IDManager::typeID IDManager::getID(T &obj) { for (auto &&o : objects) { if (o.second == obj) return o.first; } DUMPFILE(fout, "Warning : auto registering of " << this->getID() << " {\n" << *obj << "\n}"); return (addObject(obj)); } /* -------------------------------------------------------------------------- */ template inline IDManager::~IDManager() { destroy(); DUMPFILE(fout, "all " << typeinfo() << " cleaned..."); fout.close(); } /* -------------------------------------------------------------------------- */ template inline void IDManager::destroy() { if (static_pointer == nullptr) return; static_pointer->shared_pointer_objects.clear(); if (static_pointer) static_pointer.release(); } /* -------------------------------------------------------------------------- */ __END_LIBMULTISCALE__ #endif /* __LIBMULTISCALE_ID_MANAGER_HH__ */ diff --git a/src/common/ref_subset.cc b/src/common/ref_subset.cc index 85a6ff8..2b72910 100644 --- a/src/common/ref_subset.cc +++ b/src/common/ref_subset.cc @@ -1,313 +1,330 @@ /** * @file ref_subset.cc * * @author Guillaume Anciaux * * @date Mon Sep 08 23:40:22 2014 * * @brief This file contains the interface of containers of references 4 the * migration system * * @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 . * */ /* -------------------------------------------------------------------------- */ #include /* -------------------------------------------------------------------------- */ #include "lib_continuum.hh" #include "lib_dd.hh" #include "lib_md.hh" #include "lm_common.hh" #include "ref_subset.hh" #include "trace_atom.hh" #include /* -------------------------------------------------------------------------- */ __BEGIN_LIBMULTISCALE__ /* -------------------------------------------------------------------------- */ template void RefSubset::declareSentRefs( MapRefToUInt &sent_refs, MapUIntToRefList &sent_reverse) { DUMP("declare refs sent for subset " << container.getID(), DBG_INFO); // clear maps sent.clear(); sent_indexes.clear(); MapRefToUInt sent_index_in_subset; // loop over subset refs if constexpr (is_subset) { for (auto &&[i, ref] : enumerate(container)) { // get the reference to the ref "i" in subset if (sent_refs.count(ref) > 0) { if (sent_refs[ref] == UINT_MAX) continue; sent_index_in_subset[ref] = i; holes.push_back(i); } } } // loop over the sent refs to build ordered buffer list UInt idx; for (auto &&[proc, listsent] : sent_reverse) { for (auto &&ref : listsent) { if (is_subset and not sent_index_in_subset.count(ref)) continue; if constexpr (not is_subset) { idx = ref.getIndex(); } else { idx = sent_index_in_subset[ref]; } sent[proc].push_back(ref); sent_indexes[proc].push_back(idx); + IF_REF_TRACED(ref, "traced ref sent to " << proc << " for container " + << container.getID()); } for (auto &&[proc, list_sent] : sent) { DUMP("I am supposed to send " << list_sent.size() << " refs to proc " << proc << " for subset " << container.getID() << "(" << subset_id << ")", DBG_INFO); } } } /* ------------------------------------------------------------------------ */ template void RefSubset::declareRecvRefs(MapRefToUInt &masks, MapUIntToRefList &recv_refs) { DUMP("declare refs received for subset " << container.getID(), DBG_INFO); UInt mask [[gnu::unused]] = 1 << subset_id; // clear maps recv.clear(); recv_indexes.clear(); auto register_newref = [&](auto proc, auto new_index, auto &ref) { - IF_REF_TRACED(at, "traced ref received from " - << proc << " for container " << &container - << " (index is now " << new_index << ")"); + IF_REF_TRACED(ref, "traced ref received from " + << proc << " for container " << container.getID() + << " (index is now " << new_index << ")"); recv_indexes[proc].push_back(new_index); recv[proc].push_back(ref); }; for (auto &&[proc, recvlist] : recv_refs) { for (auto &&ref : recvlist) { if constexpr (is_subset) { // check if concerned with ref => register if (masks[ref] & mask) { UInt new_index; // find a new index (if possible within holes) if (!holes.empty()) { new_index = holes.back(); container.get(new_index) = ref; holes.pop_back(); } else { new_index = container.size(); container.push_back(ref); } register_newref(proc, new_index, ref); } } else { register_newref(proc, ref.getIndex(), ref); } } } for (auto &&[proc, recvlist] : recv_indexes) { DUMP("I am supposed to recv " << recvlist.size() << " refs from proc " << proc << " for subset " << container.getID() << "(" << subset_id << ")", DBG_INFO); } } /* ----------------------------------------------------------------------- */ template void RefSubset::buildMasksForSentRefs(MapRefToUInt &masks) { if constexpr (not is_subset) return; UInt mask = 1 << subset_id; for (auto &&[proc, reflist] : sent) { UInt nRefs = reflist.size(); for (UInt i = 0; i < nRefs; ++i) { Ref ref = reflist[i]; if (!masks.count(ref)) masks[ref] = mask; else masks[ref] |= mask; } } } /* ----------------------------------------------------------------------- */ template void RefSubset::packAttachedData( std::map &buffer) { for (auto &&[proc, listsent] : sent_indexes) { if (proc == UINT_MAX) continue; - for (auto &&sent_index : listsent) { + for (auto &&[i, sent_index] : enumerate(listsent)) { + if (this->attached_objects.size() > 0) + DUMP(container.getID() + << " prepares to pack to proc " << proc << " subset index " + << i << " ref index " << sent_index, + DBG_MESSAGE); for (auto &&[ptr, att_obj] : this->attached_objects) { - att_obj->packData(sent_index, buffer[proc], proc); + DUMP("prepares to pack obj " << att_obj->getID(), DBG_MESSAGE); + att_obj->packData(sent_index, buffer[proc], proc, true); } } } } -/* -------------------------------------------------------------------------- - */ +/* ------------------------------------------------------------------------ */ template void RefSubset::unpackAttachedData( std::map &buffer) { for (auto &&[proc, listrecv] : recv_indexes) { if (proc == UINT_MAX) continue; - for (auto &&recv_index : listrecv) { + for (auto &&[i, recv_index] : enumerate(listrecv)) { + if (CHECK_REF_TRACED(recv[proc][i])) + DUMP(container.getID() + << " prepares to unpack from proc " << proc << " subset index " + << i << " ref index " << recv_index, + DBG_MESSAGE); + for (auto &&[ptr, att_obj] : this->attached_objects) { - att_obj->unpackData(recv_index, buffer[proc], proc); + if (CHECK_REF_TRACED(recv[proc][i])) + DUMP("prepares to pack obj " << att_obj->getID(), DBG_MESSAGE); + att_obj->unpackData(recv_index, buffer[proc], proc, + CHECK_REF_TRACED(recv[proc][i])); } } } } /* -------------------------------------------------------------------------- */ template void RefSubset::translateMovingReferences(MapRefToRef &moved) { if constexpr (not is_subset) return; // loop over the contained refs and shift references of moved refs UInt nRefs = container.size(); for (UInt i = 0; i < nRefs; ++i) { // get the reference to the ref "i" Ref &ref = container.get(i); if (moved.count(ref)) { DUMP("translating ref in subset " << subset_id << " numero " << i << " from position " << ref << " to position " << moved[ref], DBG_INFO); - IF_REF_TRACED(ref, "translate ref " << at << " to ref " << moved[ref]); + IF_REF_TRACED(ref, container.getID() << ": translate ref(" << i << ") " + << ref << " to ref " << moved[ref]); ref = moved[ref]; } } } /* -------------------------------------------------------------------------- */ template void RefSubset::fillRemainingHoles() { if constexpr (not is_subset) return; // checks if there are remaining holes for which we need to displace // the references order and attached values while (!holes.empty()) { UInt index = holes.back(); holes.pop_back(); UInt nRefs = container.size(); if (index == nRefs - 1) // container.resize(nRefs-1); resize(nRefs - 1); else moveLastReference(index); } } /* -------------------------------------------------------------------------- */ template void RefSubset::moveReference(UInt i_src, UInt i_dest) { if (i_src == i_dest) return; IF_REF_TRACED(container.get(i_src), "traced ref moved from index " << i_src << " to index " << i_dest << " for container " << &container); DUMP("ref moved from index " << i_src << " to index " << i_dest << " for container " << &container, DBG_INFO); if constexpr (is_subset) container.get(i_dest) = container.get(i_src); for (auto &&[ptr, att_obj] : this->attached_objects) { att_obj->moveAttachedValues(i_src, i_dest); } } /* ------------------------------------------------------------------------ */ template void RefSubset::moveLastReference(UInt i_dest) { if constexpr (not is_subset) return; UInt nRefs = container.size(); if (nRefs - 1 == i_dest) return; moveReference(nRefs - 1, i_dest); // container.resize(nRefs-1); resize(nRefs - 1); } /* -------------------------------------------------------------------------- */ template UInt RefSubset::resize(UInt sz) { if constexpr (is_subset) { if (sz == UInt(-1)) LM_FATAL("this should not happen"); container.resize(sz); } else { if (sz != UInt(-1)) LM_FATAL("this should not happen"); sz = container.size(); } DUMP(container.getID() << " resizes from " << container.size() << " to " << sz, DBG_INFO); for (auto &&[ptr, att_obj] : this->attached_objects) { att_obj->resize(sz); } return sz; } /* -------------------------------------------------------------------------- */ #define DECLARE_REF_SUBSET(n, data, obj) \ template class RefSubset, \ true>; \ template class RefSubset< \ typename BOOST_PP_TUPLE_ELEM(3, 0, obj)::ContainerPoints, false>; BOOST_PP_SEQ_FOR_EACH(DECLARE_REF_SUBSET, f, LIST_ATOM_MODEL) BOOST_PP_SEQ_FOR_EACH(DECLARE_REF_SUBSET, f, LIST_DD_MODEL) BOOST_PP_SEQ_FOR_EACH(DECLARE_REF_SUBSET, f, LIST_CONTINUUM_MODEL) #undef DECLARE_REF_SUBSET __END_LIBMULTISCALE__ diff --git a/src/common/ref_subset.hh b/src/common/ref_subset.hh index 3d0b483..9a729e7 100644 --- a/src/common/ref_subset.hh +++ b/src/common/ref_subset.hh @@ -1,184 +1,182 @@ /** * @file ref_subset.hh * * @author Guillaume Anciaux * * @date Mon Sep 08 23:40:22 2014 * * @brief This file contains the interface of containers of references 4 the * migration system * * @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 . * */ #ifndef __LIBMULTISCALE_REF_SUBSET_HH__ #define __LIBMULTISCALE_REF_SUBSET_HH__ /* -------------------------------------------------------------------------- */ #include #include /* -------------------------------------------------------------------------- */ #include "attached_object.hh" #include "pack_buffer.hh" /* -------------------------------------------------------------------------- */ __BEGIN_LIBMULTISCALE__ /* -------------------------------------------------------------------------- */ class ContainerInterface; /* -------------------------------------------------------------------------- */ template class RefSubsetInterface { public: using MapUIntToUIntList = std::map>; using MapRefToRef = std::map; using MapRefToUInt = std::map; using MapUIntToRefList = std::map>; //! return the container being this reference set virtual ContainerInterface &getContainer() = 0; //! declare the refs sent to other processors virtual void declareSentRefs(MapRefToUInt &sent_refs, MapUIntToRefList &sent_reverse) = 0; //! declare the refs received from other processors virtual void declareRecvRefs(MapRefToUInt &masks, MapUIntToRefList &recv_refs) = 0; //! build the mask for the sent refs virtual void buildMasksForSentRefs(MapRefToUInt &masks) = 0; //! translate the internal references virtual void translateMovingReferences(MapRefToRef &moved) = 0; //! pack attached data virtual void packAttachedData(std::map &buffer) = 0; //! unpack attached data virtual void unpackAttachedData(std::map &buffer) = 0; //! fill remaining holes so that container is one pack virtual void fillRemainingHoles() = 0; //! resize container and attached values virtual UInt resize(UInt sz = -1) = 0; //! attach a vector of data to the subset - bool attachData(std::shared_ptr &obj); + void attachData(std::shared_ptr &obj); //! detach a vector of data from the subset void detachData(std::shared_ptr &obj); virtual bool isSubset() = 0; protected: //! list of attached objects - std::map> attached_objects; + std::map> attached_objects; }; /* -------------------------------------------------------------------------- */ template class RefSubset : public RefSubsetInterface { using Ref = typename Cont::Ref; using MapRefToUInt = std::map; using MapRefToRef = std::map; using MapUIntToRefList = std::map>; using MapUIntToUIntList = std::map>; public: RefSubset(Cont &cont, UInt subset_id) : container(cont), subset_id(subset_id){}; virtual ~RefSubset(){}; bool isSubset() override { return is_subset; } //! declare the refs sent to other processors void declareSentRefs(MapRefToUInt &sent_refs, MapUIntToRefList &sent_reverse); //! declare the refs received from other processors void declareRecvRefs(MapRefToUInt &masks, MapUIntToRefList &recv_refs); //! translate the internal references void translateMovingReferences(MapRefToRef &moved); //! fill remaining holes so that container is one pack void fillRemainingHoles() override; //! build the mask for the sent refs void buildMasksForSentRefs(MapRefToUInt &masks); //! pack attached data void packAttachedData(std::map &buffer); //! unpack attached data void unpackAttachedData(std::map &buffer); //! move last reference in the array (useful to fill holes) void moveLastReference(UInt i_dest); //! move reference in the array void moveReference(UInt i_src, UInt i_dest); //! resize container and attached values UInt resize(UInt sz = -1); //! return the container being this reference set ContainerInterface &getContainer() { return container; }; private: // list of holes created by the gone refs std::vector holes; //! container of references that constitutes the subset Cont &container; //! map new proc to the list of refs that have been sent MapUIntToRefList sent; //! map old proc to the list of refs that have been received MapUIntToRefList recv; //! map new proc to the list of refs indexes that have been sent MapUIntToUIntList sent_indexes; //! map old proc to the list of refs indexes that have been received MapUIntToUIntList recv_indexes; // id of the subset UInt subset_id; }; /* -------------------------------------------------------------------------- */ // attach a vector of data to the subset template -bool RefSubsetInterface::attachData(std::shared_ptr &obj) { +void RefSubsetInterface::attachData(std::shared_ptr &obj) { DUMP("want to register obj " << &obj, DBG_DETAIL); - for (auto &&[ptr, att_obj] : attached_objects) { - if (ptr == &obj) { - DUMP("obj " << &obj << " was already attached in the past", DBG_DETAIL); - return false; + for (auto &&[name, att_obj] : attached_objects) { + if (name == obj->getID()) { + LM_FATAL("obj " << obj->getID() << " was already attached in the past"); } } - attached_objects[obj.get()] = obj; - return true; + attached_objects[obj->getID()] = obj; } /* -------------------------------------------------------------------------- */ // detach a vector of data from the subset template void RefSubsetInterface::detachData(std::shared_ptr &obj) { - attached_objects.erase(&obj); + attached_objects.erase(obj->getID()); } /* -------------------------------------------------------------------------- */ template auto make_subset(Cont &cont, UInt subset_id) { static_assert(is_subset == std::is_same_v, "This situation should not happen"); return std::make_shared>(cont, subset_id); } __END_LIBMULTISCALE__ #endif /* __LIBMULTISCALE_REF_SUBSET_HH__ */ diff --git a/src/common/reference_manager.cc b/src/common/reference_manager.cc index a05b267..56d4b60 100644 --- a/src/common/reference_manager.cc +++ b/src/common/reference_manager.cc @@ -1,467 +1,466 @@ /** * @file reference_manager.cc * * @author Guillaume Anciaux * * @date Mon Sep 08 23:40:22 2014 * * @brief This is the manager of reference and coherency with migrations * * @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 /* -------------------------------------------------------------------------- */ #include "lib_continuum.hh" #include "lib_dd.hh" #include "lib_md.hh" #include "lm_common.hh" #include "lm_communicator.hh" #include "reference_manager.hh" /* -------------------------------------------------------------------------- */ - __BEGIN_LIBMULTISCALE__ /* -------------------------------------------------------------------------- */ template ReferenceManager::ReferenceManager( typename Ref::Domain::ContainerPoints &global_container) { auto name = global_container.getID(); subsets[name] = make_subset(global_container, 0); } /* -------------------------------------------------------------------------- */ template void ReferenceManager::printPackBuffersStatus() { #ifndef LM_OPTIMIZED for (auto &&[proc, buf] : buffers_torecv) { DUMP("from " << proc << " receive buffer status: " << buf.size() << " " << buf.remainingSize(), DBG_INFO); } // same thing for sending buffers for (auto &&[proc, buf] : buffers_tosend) { DUMP("to " << proc << " send buffer status: " << buf.size() << " " << buf.remainingSize(), DBG_INFO); } #endif // LM_OPTIMIZED } /* -------------------------------------------------------------------------- */ template void ReferenceManager::clearPackBuffers() { // clear buffers buffers_tosend.clear(); buffers_torecv.clear(); // create a packbuffer for each proc I should receive from // and clear it for (auto &&[proc, ref] : newrefs) { PackBuffer &buf = buffers_torecv[proc]; buf.clear(); } // same thing for sending buffers for (auto &&[proc, ref] : sent_byproc) { PackBuffer &buf = buffers_tosend[proc]; buf.clear(); } } /* -------------------------------------------------------------------------- */ template void ReferenceManager::updateRefSubSets() { MapRefToUInt masks_send; MapRefToUInt masks_recv; //#ifndef LM_OPTIMIZED comm_group->synchronize(); // MPI_Barrier(worldCom); //#endif if (have_changed) { DUMP("have_changed = " << have_changed, DBG_INFO); } // do a bilan and checkup of coherency state // printBilan(); #ifndef LM_OPTIMIZED std::map n_tot_subset; // count number of refs in before subsets for (auto &&[name, s] : subsets) { UInt n_tot = s->getContainer().size(); comm_group->reduce(&n_tot, 1, "n_tot", OP_SUM, 0); // MPI_Reduce(&n, &n_tot, 1, MPI_INT, MPI_SUM, 0, worldCom); if (lm_my_proc_id == 0) n_tot_subset[name] = n_tot; } #endif // declare the ref gone to other processors to attached subsets for (auto &&[name, s] : subsets) s->declareSentRefs(sent, sent_byproc); // global_set.declareSentAtoms(sent, sent_byproc); // build the mask of ownership for each ref for (auto &&[name, s] : subsets) s->buildMasksForSentRefs(masks_send); // initialize the buffers clearPackBuffers(); // pack the masks to a packbuffer packMasks(masks_send); // exchange the masks exchangeBuffers(buffers_tosend, buffers_torecv); //#ifndef LM_OPTIMIZED comm_group->synchronize(); // MPI_Barrier(worldCom); //#endif // unpack the masks to a MapRefToInt unpackMasks(masks_recv); // translate references for (auto &&[name, s] : subsets) s->translateMovingReferences(moved); // declare the received refs and mask for (auto &&[name, s] : subsets) s->declareRecvRefs(masks_recv, newrefs); // global_set.declareRecvAtoms(newatoms); // initialize the buffers clearPackBuffers(); printPackBuffersStatus(); // pack attached data for (auto &&[name, s] : subsets) s->packAttachedData(buffers_tosend); // global_set.packAttachedData(buffers_tosend); // do the communication exchangeBuffers(buffers_tosend, buffers_torecv); //#ifndef LM_OPTIMIZED comm_group->synchronize(); // MPI_Barrier(worldCom); //#endif printPackBuffersStatus(); // unpack attached data for (auto &&[name, s] : subsets) s->unpackAttachedData(buffers_torecv); // global_set.unpackAttachedData(buffers_torecv); // translate references for (auto &&[name, s] : subsets) s->fillRemainingHoles(); UInt current_size [[gnu::unused]]; for (auto &&[name, s] : subsets) { if (not s->isSubset()) current_size = s->resize(); } // for (auto &&[name, s] : subsets) { // if (s->isSubset()) // s->resize(current_size); // } #ifndef LM_OPTIMIZED // count number of refs in subsets for (auto &&[name, s] : subsets) { UInt n_tot = s->getContainer().rows(); DUMP("container:" << name << " size: " << n_tot, DBG_INFO); comm_group->reduce(&n_tot, 1, "n_tot", OP_SUM, 0); if (lm_my_proc_id == 0) if (n_tot != n_tot_subset[name]) LM_FATAL("subset " << name << " mismatch of ref count " << n_tot << " " << n_tot_subset[name]); } #endif // clear the stl structures moved.clear(); sent.clear(); sent_byproc.clear(); newrefs.clear(); have_changed = false; DUMP("end update ref", DBG_INFO); } /* -------------------------------------------------------------------------- */ template void ReferenceManager::exchangeBuffers(BufferMap &toSend, BufferMap &toRecv) { std::map send_requests; /* ------------------------------------------------------------------------ */ // now send the buffers through network /* ------------------------------------------------------------------------ */ // sending of the buffers for (auto &&[proc, buf_send] : toSend) { if (proc == UINT_MAX) continue; UInt size = buf_send.size(); auto buffer = buf_send.buffer(); send_requests[proc] = comm_group->isend(buffer, size, proc, "ref_manager send"); // MPI_Isend(buffer, size, MPI_CHAR, proc, rank, worldCom, &request); } // reception async of the vectors for (auto &&[proc, buffer_to_recv] : toRecv) { if (proc == UINT_MAX) continue; DUMP("probing from proc " << proc, DBG_INFO); int size = comm_group->probe( proc, "probing ref_manager recv"); // MPI_Probe(proc, proc, // worldCom, &status); PackBuffer &buf = toRecv[proc]; buf.resize(size); auto buffer = buffer_to_recv.buffer(); DUMP("receiving from proc " << proc << " " << size << " bytes", DBG_INFO); comm_group->receive( buffer, size, proc, "recv ref_manager buffer"); //(buffer, size, MPI_CHAR, proc, proc, // worldCom, &status); DUMP("received packet from proc " << proc << " of size " << size << " bytes", DBG_DETAIL); } // now I wait for requests to be acheived for (auto &&[proc, request] : send_requests) { if (proc == int(UINT_MAX)) continue; #ifndef LM_OPTIMIZED UInt size = toRecv[proc].size(); DUMP("waiting to send a packet to proc " << proc << " of size " << size << " bytes", DBG_INFO); #endif // LM_OPTIMIZED comm_group->wait(request); DUMP("sent packet to proc " << proc << " of size " << size << " bytes", DBG_DETAIL); } } /* -------------------------------------------------------------------------- */ template void ReferenceManager::packMasks(MapRefToUInt &masks) { // pack the masks for (auto &&[proc, sendlist] : sent_byproc) { for (UInt i = 0; i < sendlist.size(); ++i) { Ref at = sendlist[i]; buffers_tosend[proc] << masks[at]; } } } /* -------------------------------------------------------------------------- */ template void ReferenceManager::unpackMasks(MapRefToUInt &masks) { for (auto &&[proc, recvlist] : newrefs) { for (UInt i = 0; i < recvlist.size(); ++i) { Ref at = recvlist[i]; UInt m; buffers_torecv[proc] >> m; masks[at] = m; } } } /* -------------------------------------------------------------------------- */ template void ReferenceManager::printBilan() { std::map nb_sent_by_proc; if (moved.size()) DUMP(moved.size() << " refs moved in memory", DBG_WARNING); for (auto &&[proc, refs] : newrefs) DUMP(refs.size() << " refs received from proc " << proc, DBG_WARNING); for (auto &&[proc, sent] : sent_byproc) DUMP(sent.size() << " refs sent to proc " << proc, DBG_WARNING); for (auto &&it : sent) ++nb_sent_by_proc[it.second]; for (auto &&[proc, sentlist] : sent_byproc) { if (sentlist.size() != nb_sent_by_proc[proc]) LM_FATAL("sent_by_proc registered " << sentlist.size() << " refs as sent but sent registered " << nb_sent_by_proc[proc] << " fatal inconstistency in send maps"); } } /* -------------------------------------------------------------------------- */ template template std::enable_if_t> ReferenceManager::addSubset(Cont &container) { auto name = container.getID(); auto it = subsets.find(name); if (it != subsets.end()) { DUMP(name << " was already registered", DBG_INFO); return; } // UInt new_subset_id = subsets.size(); auto subset = make_subset(container, subsets.size()); DUMP("adding a subset to managed subset named " << name << " " << &container, DBG_INFO_STARTUP); subsets[name] = subset; } /* -------------------------------------------------------------------------- */ template void ReferenceManager::removeSubset(const std::string &name) { subsets.erase(name); #ifndef LM_OPTIMIZED UInt size = subsets.size(); DUMP("removing a subset from managed subset at position " << size << " and named " << name, DBG_INFO); #endif } /* ------------------------------------------------------------------------ */ // template // void ReferenceManager::attachVector(ContainerArray &v, // ContainerInterface &cont) { // DUMP("try to attach vector " << &v << " to subset " << &cont, // DBG_MESSAGE); // UInt i = 0; // for (auto &&[name, subset] : subsets) { // if (&subset->getContainer() == &cont) { // auto atV = std::make_shared>(v); // if (subset->attachData(*atV)) // DUMP("attaching vector " << &v << " to subset " << cont.getID() // << " " << &cont, // DBG_MESSAGE); // break; // } // ++i; // } // if (i == subsets.size()) { // DUMP("attaching vector could not be done : did not find its subset // " // << &cont << "\nRegistered vectors are ", // DBG_MESSAGE); // for (auto &s : subsets) { // DUMP(s.first << " : " << &s.second->getContainer() << "\n", // DBG_MESSAGE); // } // LM_FATAL("abort"); // } // } /* ------------------------------------------------------------------------ */ // template // void ReferenceManager::attachVector(ContainerArray &v, // GlobalContainer &) { // AttachedVector *atV = new AttachedVector(v); // global_set.attachData(*atV); // } /* --------------------------------------------------------------------- */ // template // void ReferenceManager::detachVector(ContainerArray &, // ContainerInterface &) { // LM_TOIMPLEMENT; // } // /* -------------------------------------------------------------------------- // */ template void ReferenceManager::setCommGroup(CommGroup &group) { comm_group = &group; } // worldCom = comm.getMPIComm(); // if (worldCom != MPI_COMM_NULL) // MPI_Comm_rank(worldCom, &rank); // } #define DECLARE_REF_MANAGER(n, data, obj) \ template class ReferenceManager; /* -------------------------------------------------------------------------- */ template void ReferenceManager::acquireContext(const LMObject &obj) { this->setCommGroup(obj.getCommGroup()); } /* -------------------------------------------------------------------------- */ template void ReferenceManager::addSubset(ContainerInterface &cont) { this->addSubset(cont); } /* -------------------------------------------------------------------------- */ BOOST_PP_SEQ_FOR_EACH(DECLARE_REF_MANAGER, f, LIST_ATOM_MODEL) BOOST_PP_SEQ_FOR_EACH(DECLARE_REF_MANAGER, f, LIST_DD_MODEL) BOOST_PP_SEQ_FOR_EACH(DECLARE_REF_MANAGER, f, LIST_CONTINUUM_MODEL) #undef DECLARE_REF_MANAGER __END_LIBMULTISCALE__ diff --git a/src/coupling/duo_distributed_vector.hh b/src/coupling/duo_distributed_vector.hh index c2d621f..1109b0c 100644 --- a/src/coupling/duo_distributed_vector.hh +++ b/src/coupling/duo_distributed_vector.hh @@ -1,270 +1,275 @@ /** * @file duo_distributed_vector.hh * * @author Guillaume Anciaux * * @date Mon Oct 28 19:23:14 2013 * * @brief Migration safe representation of array of references * * @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. * */ #ifndef __DUO_DISTRIBUTED_VECTOR_H__ #define __DUO_DISTRIBUTED_VECTOR_H__ /* -------------------------------------------------------------------------- */ #include /* -------------------------------------------------------------------------- */ #include "comm_group.hh" #include "container_array.hh" #include "ref_subset.hh" /* -------------------------------------------------------------------------- */ __BEGIN_LIBMULTISCALE__ /* -------------------------------------------------------------------------- */ /** Class duodistributedvector allows redistribution of vectors, with possiblity to adapt the migration procedures (creating holes in the commuincation scheme) and eventually unfragmenting it */ class DuoDistributedVector : public AttachedObject { public: using MapUIntToUIntList = std::map>; using MapUIntToUInt = std::map; using PairUIntList = std::vector>; ~DuoDistributedVector(); DuoDistributedVector(UInt lsize, UInt tsize, const std::string &my_name); //! add an index for the communication to a given proc inline void addIndex(UInt index, UInt proc); //! remove an index for the communication : return associated duo proc inline UInt removeIndex(UInt index); //! find the real index given proc and com_index inline UInt findRealIndex(UInt com_index, UInt duo_proc); //! set the duo proc of a given index (return 1 if success 0 otherwise) inline UInt setDuoProc(UInt index, UInt proc); //! update both sides of the migrated operations UInt synchronizeMigration(CommGroup &groupAtomic, CommGroup &groupFE); //! function used to synchronize vectors (by summing them) void synchronizeVectorBySum(const std::string &name_vec, ContainerArray &vec, CommGroup &group1, CommGroup &group2); //! function used to redistribute a vector void distributeVector(const std::string &name_vec, ContainerArray &vec, CommGroup &group1, CommGroup &group2); //! print to file the current redistribution scheme (with positions info) void print(const std::string &prefix); //! print to screen summary of com scheme void printSummary(); //! set the size of the duo vector void setSize(UInt lsize, UInt tsize); //! clear the duo (for reuse) void clear(); private: //! update the communication list for the sent indexes void updateSent(); //! update the communication list for the recv indexes void updateRecv(); //! update the communication list for the moved indexes void updateMoved(); public: - inline void packData(UInt index, PackBuffer &buffer, UInt toproc); - inline void unpackData(UInt index, PackBuffer &buffer, UInt fromproc); - inline void moveAttachedValues(UInt i_src, UInt i_dest); - inline void resize(UInt sz); + inline void packData(UInt index, PackBuffer &buffer, UInt toproc, + bool verbose = false) override; + inline void unpackData(UInt index, PackBuffer &buffer, UInt fromproc, + bool verbose = false) override; + inline void moveAttachedValues(UInt i_src, UInt i_dest) override; + inline void resize(UInt sz) override; + + LMID getID() const override { LM_TOIMPLEMENT; }; private: //! map proc to communication indexes list MapUIntToUIntList com_list; //! map index to communication indexes MapUIntToUInt duo_index_map; //! map index to duo proc MapUIntToUInt duo_proc_map; //! list of indexes sent by migration MapUIntToUIntList sent; //! list of processors indexes were sent to MapUIntToUIntList sent_procs; //! list of indexes received by migration MapUIntToUIntList received; //! moved indexes in the container PairUIntList moved; //! total size of vectors that can be redistributed UInt totalsize; //! actual local size of vecotr that can be redistributed UInt localsize; //! flag used to say that migration have occured bool something_changed; //! descriptive name std::string name; }; /* -------------------------------------------------------------------------- */ inline UInt DuoDistributedVector::setDuoProc(UInt index, UInt proc) { if (duo_index_map.count(index)) { DUMP("index " << index << " was already associated with duo proc " << duo_proc_map[index], DBG_INFO); return 0; } addIndex(index, proc); return 1; } /* -------------------------------------------------------------------------- */ inline UInt DuoDistributedVector::findRealIndex(UInt com_index, UInt duo_proc) { LM_ASSERT(com_list.count(duo_proc), "processor " << duo_proc << " not part of redistribution scheme: abort"); LM_ASSERT(com_index < com_list[duo_proc].size(), "invalid com_index: abort " << com_index << " " << com_list[duo_proc].size()); return com_list[duo_proc][com_index]; } /* -------------------------------------------------------------------------- */ inline UInt DuoDistributedVector::removeIndex(UInt index) { LM_ASSERT(duo_index_map.count(index), "index not present in communication tab"); UInt com_index = duo_index_map[index]; UInt duo_proc = duo_proc_map[index]; duo_index_map.erase(index); duo_proc_map.erase(index); std::vector &coms = com_list[duo_proc]; UInt last_index = coms.back(); coms.pop_back(); // put the last com to the position that have been removed DUMP("remove index " << index << " com index " << com_index << " from proc " << duo_proc, DBG_DETAIL); if (com_index < coms.size()) { coms[com_index] = last_index; duo_index_map[last_index] = com_index; DUMP("displace index " << last_index << " from com index " << coms.size() << " to com index " << com_index << " for proc " << duo_proc << " duo is " << name, DBG_DETAIL); } return duo_proc; } /* -------------------------------------------------------------------------- */ inline void DuoDistributedVector::addIndex(UInt index, UInt duo_proc) { LM_ASSERT(!duo_index_map.count(index), "index " << index << " already present in communication tab" << " of proc " << duo_proc_map[index]); std::vector &coms = com_list[duo_proc]; UInt com_index = coms.size(); coms.push_back(index); duo_proc_map[index] = duo_proc; duo_index_map[index] = com_index; DUMP("add index " << index << " com index " << com_index << " from proc " << duo_proc << " duo is " << name, DBG_DETAIL); } /* -------------------------------------------------------------------------- */ inline void DuoDistributedVector::packData(UInt index, PackBuffer &buffer, - UInt proc) { + UInt proc, bool verbose) { LM_ASSERT(duo_index_map.count(index), "not found index : abort"); LM_ASSERT(duo_proc_map.count(index), "not found index : abort"); UInt fem_proc = duo_proc_map[index]; buffer << fem_proc; sent[fem_proc].push_back(index); sent_procs[fem_proc].push_back(proc); DUMP("register index " << index << " to sent_list" << " from proc " << fem_proc << " giving new owner as " << proc << " duo is " << name, DBG_DETAIL); } /* -------------------------------------------------------------------------- */ inline void DuoDistributedVector::unpackData(UInt index, PackBuffer &buffer, - UInt proc [[gnu::unused]]) { + UInt proc [[gnu::unused]], + bool verbose) { UInt fem_proc; buffer >> fem_proc; received[fem_proc].push_back(index); DUMP("register index " << index << " to receive_list" << " from proc " << fem_proc << " giving old owner as " << proc << " duo is " << name, DBG_DETAIL); } /* -------------------------------------------------------------------------- */ inline void DuoDistributedVector::moveAttachedValues(UInt i_src, UInt i_dest) { std::pair p(i_src, i_dest); moved.push_back(p); } /* -------------------------------------------------------------------------- */ inline void DuoDistributedVector::resize(UInt sz) { localsize = sz; } /* -------------------------------------------------------------------------- */ inline std::shared_ptr make_attached(DuoDistributedVector &obj) { LM_TOIMPLEMENT; // return std::make_shared>(obj); } /* -------------------------------------------------------------------------- */ __END_LIBMULTISCALE__ #endif /* __DUO_DISTRIBUTED_VECTOR_H__ */ diff --git a/src/filter/compute_centro_symmetry.cc b/src/filter/compute_centro_symmetry.cc index 7c168a8..2c8cc78 100644 --- a/src/filter/compute_centro_symmetry.cc +++ b/src/filter/compute_centro_symmetry.cc @@ -1,461 +1,461 @@ /** * @file compute_centro_symmetry.cc * * @author Guillaume Anciaux * * @date Wed Jul 09 21:59:47 2014 * * @brief This computes the centro symmetry criterion over a set of points * * @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 . * */ /* -------------------------------------------------------------------------- */ #include "compute_centro_symmetry.hh" #include "lib_continuum.hh" #include "lib_dd.hh" #include "lib_md.hh" #include "lm_common.hh" #include "ref_point_data.hh" #include "trace_atom.hh" /* -------------------------------------------------------------------------- */ #define SWAP(a, b) \ { \ tmp = a; \ a = b; \ b = tmp; \ } #define ISWAP(a, b) \ { \ itmp = a; \ a = b; \ b = itmp; \ } /* -------------------------------------------------------------------------- */ __BEGIN_LIBMULTISCALE__ ComputeCentroSymmetry::ComputeCentroSymmetry(const std::string &name) : LMObject(name) { this->createInput("input"); this->createOutput("output"); this->getOutput("output") = ContainerArray(this->getID() + ":" + "output"); } /* -------------------------------------------------------------------------- */ ComputeCentroSymmetry::~ComputeCentroSymmetry() {} /* -------------------------------------------------------------------------- */ inline void select(UInt k, UInt n, Real *arr) { UInt i, ir, j, l, mid; Real a, tmp; arr--; l = 1; ir = n; for (;;) { if (ir <= l + 1) { if (ir == l + 1 && arr[ir] < arr[l]) { SWAP(arr[l], arr[ir]); } return; } else { mid = (l + ir) >> 1; SWAP(arr[mid], arr[l + 1]); if (arr[l] > arr[ir]) { SWAP(arr[l], arr[ir]); } if (arr[l + 1] > arr[ir]) { SWAP(arr[l + 1], arr[ir]); } if (arr[l] > arr[l + 1]) { SWAP(arr[l], arr[l + 1]); } i = l + 1; j = ir; a = arr[l + 1]; for (;;) { do i++; while (arr[i] < a); do j--; while (arr[j] > a); if (j < i) break; SWAP(arr[i], arr[j]); } arr[l + 1] = arr[j]; arr[j] = a; if (j >= k) ir = j - 1; if (j <= k) l = i; } } } /* ---------------------------------------------------------------------- */ // function which sorts the array arr and move accordingly the iarr entries // inline void select2(UInt k, UInt n, Real *arr, UInt *iarr) { UInt i, ir, j, l, mid, ia, itmp; Real a, tmp; arr--; iarr--; l = 1; ir = n; for (;;) { if (ir <= l + 1) { if (ir == l + 1 && arr[ir] < arr[l]) { SWAP(arr[l], arr[ir]); ISWAP(iarr[l], iarr[ir]); } return; } else { mid = (l + ir) >> 1; SWAP(arr[mid], arr[l + 1]); ISWAP(iarr[mid], iarr[l + 1]); if (arr[l] > arr[ir]) { SWAP(arr[l], arr[ir]); ISWAP(iarr[l], iarr[ir]); } if (arr[l + 1] > arr[ir]) { SWAP(arr[l + 1], arr[ir]); ISWAP(iarr[l + 1], iarr[ir]); } if (arr[l] > arr[l + 1]) { SWAP(arr[l], arr[l + 1]); ISWAP(iarr[l], iarr[l + 1]); } i = l + 1; j = ir; a = arr[l + 1]; ia = iarr[l + 1]; for (;;) { do i++; while (arr[i] < a); do j--; while (arr[j] > a); if (j < i) break; SWAP(arr[i], arr[j]); ISWAP(iarr[i], iarr[j]); } arr[l + 1] = arr[j]; arr[j] = a; iarr[l + 1] = iarr[j]; iarr[j] = ia; if (j >= k) ir = j - 1; if (j <= k) l = i; } } } /* -------------------------------------------------------------------------- */ template inline Real ComputeCentroSymmetry::computeCentroSymmetry( typename Cont::Ref &at, Cont &cont, Real *distsq, UInt *nearest, UInt maxneigh [[gnu::unused]], Real *pairs, ContNeigh &neighbors) { Vector X = at.position().template cast(); UInt numneigh = 0; Cube &bbox = cont.getBoundingBox(); auto d_size = bbox.getSize(); auto d_inv_size_half = bbox.getInvSizeHalf(); UInt tested = 0; for (auto &at_neighbor : neighbors.to(X)) { auto X2 = at_neighbor.position(); if (X == X2) { continue; } Real d2 = RefPoint::template computeDistance( at.position(), X2, bbox); { - CHECK_REF_TRACED(at) { + if (CHECK_REF_TRACED(at)) { DUMP(" testing atom " << at_neighbor, DBG_MESSAGE); ++tested; } } if (d2 > rcut2) continue; // rejected DUMP(numneigh << "-th atom selected " << (X[0] - X2[0]) / 3.614999975398666 << " " << (X[1] - X2[1]) / 3.614999975398666 << " " << (X[2] - X2[2]) / 3.614999975398666 << " " << " d = " << sqrt(d2) / 3.614999975398666 << " id = " << at_neighbor.getIndex(), DBG_DETAIL); distsq[numneigh] = d2; nearest[numneigh] = at_neighbor.getIndex(); ++numneigh; } // if not 12 neighbors, centro = 0.0 if (numneigh < 12) return 0.0; LM_ASSERT(numneigh <= maxneigh, "internal error " << numneigh << " !< " << maxneigh); DUMP("found " << numneigh << " atoms as first neighbors !", DBG_DETAIL); // store 12 nearest neighs in 1st 12 locations of distsq and nearest select2(12, numneigh, distsq, nearest); { - CHECK_REF_TRACED(at) { + if (CHECK_REF_TRACED(at)) { DUMP("CONCERNING atom " << at, DBG_MESSAGE); DUMP(" I found " << numneigh << "/" << tested << " closest neighbors", DBG_MESSAGE); for (UInt j = 0; j < numneigh; j++) { UInt jj = nearest[j]; typename Cont::Ref atjj = cont.get(jj); DUMP("d = " << distsq[j] << " " << atjj, DBG_MESSAGE); } DUMP("endlist", DBG_MESSAGE); DUMP("my bbox is " << bbox, DBG_MESSAGE); } } // R = Ri + Rj for each of 66 i,j pairs among 12 neighbors // pairs = squared length of each R DUMP("looping over neighbors of atom : " << at, DBG_DETAIL); typename Cont::Ref atjj; typename Cont::Ref atkk; UInt n = 0; UInt j, jj, k, kk; for (j = 0; j < 12; j++) { jj = nearest[j]; atjj = cont.get(jj); // auto X2 = atjj.position(); for (k = j + 1; k < 12; k++) { kk = nearest[k]; atkk = cont.get(kk); pairs[n] = 0.0; for (UInt i = 0; i < Dim; ++i) { Real del1 = atjj.position()[i] - X[i]; if ((pbc_x && i == 0) || (pbc_y && i == 1) || (pbc_z && i == 2)) { int offset = static_cast(del1 * d_inv_size_half[i]); del1 -= offset * d_size[i]; } Real del2 = atkk.position()[i] - X[i]; if ((pbc_x && i == 0) || (pbc_y && i == 1) || (pbc_z && i == 2)) { int offset = static_cast(del2 * d_inv_size_half[i]); del2 -= offset * d_size[i]; } Real del = del1 + del2; pairs[n] += del * del; } ++n; } } // store 6 smallest pair distances in 1st 6 locations of pairs select(6, 66, pairs); // centrosymmetry = sum of 6 smallest squared values Real value = 0.0; for (j = 0; j < 6; j++) { DUMP("value(" << value << ")+= " << pairs[j], DBG_DETAIL); { - CHECK_REF_TRACED(at) { + if (CHECK_REF_TRACED(at)) { DUMP("add " << pairs[j] << " to " << value << " giving " << value + pairs[j], DBG_MESSAGE); } } value += pairs[j]; } DUMP("centrosymmetry computed for atom " << X << " : " << value, DBG_DETAIL); return value; } /* -------------------------------------------------------------------------- */ template inline void ComputeCentroSymmetry::computePerAtomCentroSymmetry(Cont &cont, ContNeigh &contNeigh) { UInt index = 0; UInt maxneigh = 300; Real *distsq = new Real[maxneigh]; UInt *nearest = new UInt[maxneigh]; Real pairs[66]; for (auto &at : cont) { if (index % 100000 == 0) DUMP("computing centrosymmetry for atom " << index << "/" << cont.size(), DBG_INFO); Real value = computeCentroSymmetry( at, cont, distsq, nearest, maxneigh, pairs, contNeigh); this->getOutputAsArray().push_back(value); } } /* -------------------------------------------------------------------------- */ template inline void ComputeCentroSymmetry::computePerAtomCentroSymmetry(Cont &cont, ContNeigh &contNeigh) { UInt pbc[3] = {0, 0, 0}; for (UInt i = 0; i < Dim; ++i) pbc[i] = pbc_flag[i]; if (pbc[0] == 0 && pbc[1] == 0 && pbc[2] == 0) computePerAtomCentroSymmetry(cont, contNeigh); if (pbc[0] == 1 && pbc[1] == 0 && pbc[2] == 0) computePerAtomCentroSymmetry(cont, contNeigh); if (pbc[0] == 1 && pbc[1] == 1 && pbc[2] == 0) computePerAtomCentroSymmetry(cont, contNeigh); if (pbc[0] == 1 && pbc[1] == 1 && pbc[2] == 1) computePerAtomCentroSymmetry(cont, contNeigh); if (pbc[0] == 1 && pbc[1] == 0 && pbc[2] == 1) computePerAtomCentroSymmetry(cont, contNeigh); if (pbc[0] == 0 && pbc[1] == 1 && pbc[2] == 0) computePerAtomCentroSymmetry(cont, contNeigh); if (pbc[0] == 0 && pbc[1] == 1 && pbc[2] == 1) computePerAtomCentroSymmetry(cont, contNeigh); if (pbc[0] == 0 && pbc[1] == 0 && pbc[2] == 1) computePerAtomCentroSymmetry(cont, contNeigh); } /* -------------------------------------------------------------------------- */ template void ComputeCentroSymmetry::build(Cont &cont) { constexpr UInt Dim = Cont::Dim; LM_ASSERT(Dim == 3, "cannot work if not in 3D!"); rcut2 = rcut * rcut; Cube cube = cont.getBoundingBox(); DUMP("compute centro BBOX " << cube, DBG_INFO); for (UInt i = 0; i < Dim; ++i) { cube.setXmin(i, cube.getXmin()[i] - this->skin[i]); cube.setXmax(i, cube.getXmax()[i] + this->skin[i]); } this->clear(); using ContNeigh = ContainerNeighborAtoms; ContNeigh neighbors(this->getID(), cube, rcut); neighbors.fillGrid(cont); // neighbors.setPBCFlag(pbc_flag); #ifdef TRACE_ATOM UInt indexes[3]; neighbors.coordinates2CartesianIndex(trace_ref, indexes); DUMP("box indexes containing the traced atom " << indexes[0] << " " << indexes[1] << " " << indexes[2], DBG_MESSAGE); UInt index; neighbors.coordinates2Index(trace_ref, index); DUMP("box index containing the traced atom " << indexes, DBG_MESSAGE); DUMP("block outside contains " << neighbors.getBlockOutside().size(), DBG_MESSAGE) #endif this->computePerAtomCentroSymmetry(cont, neighbors); } /* -------------------------------------------------------------------------- */ /* LMDESC CENTRO This computes the centrosymmetry for a set of points centro-symmetry criterion is detailed in Kelchner, C. L., Plimpton, S. J., Hamilton, J. C., Nov 1998. Dislocation nucleation and defect structure during surface indentation. Phys. Rev. B 58 (17), 11085-11088. */ /* LMEXAMPLE COMPUTE centro CENTRO INPUT md RCUT 3.7 */ inline void ComputeCentroSymmetry::declareParams() { ComputeInterface::declareParams(); /* LMKEYWORD RCUT Specify the cutoff radius to be used for the neighbor detection */ this->parseKeyword("RCUT", rcut); /* LMKEYWORD PBC Specify the pbc flags (one per direction) for neighbor detection */ this->parseVectorKeyword("PBC", spatial_dimension, pbc_flag, VEC_DEFAULTS(false, false, false)); /* LMKEYWORD SKIN Specify the additional tolerance to be given per spatial direction to ensure a large enough search grid */ this->parseVectorKeyword("SKIN", spatial_dimension, skin, VEC_DEFAULTS(0., 0., 0.)); } /* -------------------------------------------------------------------------- */ DECLARE_COMPUTE_MAKE_CALL(ComputeCentroSymmetry) __END_LIBMULTISCALE__ diff --git a/src/md/lammps/domain_lammps.cc b/src/md/lammps/domain_lammps.cc index 56ac11a..e8dcbd3 100644 --- a/src/md/lammps/domain_lammps.cc +++ b/src/md/lammps/domain_lammps.cc @@ -1,979 +1,979 @@ /** * @file domain_lammps.cc * * @author Guillaume Anciaux * @author Jaehyun Cho * @author Till Junge * * @date Thu Jul 31 22:41:23 2014 * * @brief This is the model wrapping LAMMPS * * @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. * */ /* -------------------------------------------------------------------------- */ //#define CHECK_STABILITY /* -------------------------------------------------------------------------- */ #include "domain_lammps.hh" #include "factory_multiscale.hh" #include "import_lammps.hh" #include "lammps_common.hh" #include "lm_common.hh" #include "reference_manager_lammps.hh" #include "trace_atom.hh" /* -------------------------------------------------------------------------- */ // LAMMPS includes #include #include #include #include #include #include #include #include #include #include #include #include #include /* -------------------------------------------------------------------------- */ #include #include #include #include /* -------------------------------------------------------------------------- */ #include "compute_lammps.hh" #include "geometry_manager.hh" #include "lm_communicator.hh" #include "lm_parser.hh" /* -------------------------------------------------------------------------- */ LAMMPS_NS::LAMMPS *lammps_main_object = NULL; /* -------------------------------------------------------------------------- */ __BEGIN_LIBMULTISCALE__ template void DomainLammps::init() { DomainAtomic, Dim>::init(); if (!this->getCommGroup().amIinGroup()) return; this->initModel(); this->update->whichflag = 1; LAMMPS_NS::LAMMPS::init(); this->update->integrate->setup(); this->initDOFs(); this->ref_manager->setState(true); this->atoms.getBoundingBox() = this->getDomainBoundingBox(); } /* -------------------------------------------------------------------------- */ template ArrayView DomainLammps::gradient() { return getField(_force); } /* -------------------------------------------------------------------------- */ template ArrayView DomainLammps::primal() { return getField(_position); } /* -------------------------------------------------------------------------- */ template auto DomainLammps::primalTimeDerivative() -> ArrayView { return getField(_velocity); } /* -------------------------------------------------------------------------- */ template auto DomainLammps::mass() -> ArrayView { return getField(_mass); } /* -------------------------------------------------------------------------- */ template auto DomainLammps::acceleration() -> ArrayView { return getField(_acceleration); } /* -------------------------------------------------------------------------- */ template auto DomainLammps::getField(FieldType field_type) -> ArrayView { switch (field_type) { case _position0: return ArrayView(&lammps_main_object->atom->x0[0][0], lammps_main_object->atom->nlocal, 3); case _velocity: return ArrayView(&lammps_main_object->atom->v[0][0], lammps_main_object->atom->nlocal, 3); case _position: return ArrayView(&lammps_main_object->atom->x[0][0], lammps_main_object->atom->nlocal, 3); case _force: return ArrayView(&lammps_main_object->atom->f[0][0], lammps_main_object->atom->nlocal, 3); case _mass: return ArrayView(&lammps_main_object->atom->mass[0], lammps_main_object->atom->ntypes, 1); case _acceleration: return ArrayView(this->acceleration_field.data(), lammps_main_object->atom->nlocal, 3); default: LM_FATAL("Not accessible field: " << field_type); } } /* -------------------------------------------------------------------------- */ template void DomainLammps::loadLammpsConfigFile() { // forward LM variables to lammps auto &vars = Parser::getAlgebraicVariables(); for (auto &&[varname, val] : vars) { std::stringstream varval; varval << std::setprecision(16) << std::scientific << val; LAMMPS_NS::Variable *lmpvar = LAMMPS_NS::LAMMPS::input->variable; char cvarname[256]; varname.copy(cvarname, 255); cvarname[varname.length()] = 0; char cvarval[30]; varval.str().copy(cvarval, 29); cvarval[varval.str().length()] = 0; if (lmpvar->find(cvarname) >= 0) LM_FATAL("cannot export LibMultiscale variable " << varname << " to LAMMPS since variable already exists"); DUMP(cvarname << " = " << cvarval, DBG_INFO); lmpvar->add_equal(cvarname, cvarval); } std::string file_to_read = lammpsfile; DUMP("lammpsfile" << " = " << lammpsfile, DBG_INFO); if (create_header_flag && this->getCommGroup().getMyRank() == 0) { std::string unit_name; if (code_unit_system == metal_unit_system) unit_name = "metal"; else if (code_unit_system == atomic_unit_system) unit_name = "real"; else if (code_unit_system == real_unit_system) unit_name = "si"; else LM_FATAL("unknown units"); file_to_read = "./generated_lammps.config"; std::ofstream file_lammps(file_to_read.c_str()); file_lammps << "units " << unit_name << std::endl; file_lammps << "dimension " << Dim << std::endl; file_lammps << "atom_style atomic" << std::endl; file_lammps << "boundary "; for (UInt i = 0; i < 3; ++i) { file_lammps << boundaries[i] << " "; } file_lammps << std::endl; file_lammps << "lattice " << lattice << " " << std::setprecision(15) << static_cast(lattice_size); file_lammps << " orient x "; for (UInt i = 0; i < 3; ++i) file_lammps << lattice_orient_x[i] << " "; if (Dim > 1) { file_lammps << "orient y "; for (UInt i = 0; i < 3; ++i) file_lammps << lattice_orient_y[i] << " "; } if (Dim == 3) { file_lammps << "orient z "; for (UInt i = 0; i < 3; ++i) file_lammps << lattice_orient_z[i] << " "; } file_lammps << "spacing "; for (UInt i = 0; i < 3; ++i) file_lammps << lattice_spacing[i] << " "; file_lammps << "origin "; for (UInt i = 0; i < Dim; ++i) file_lammps << Real(lattice_origin[i]) << " "; for (UInt i = Dim; i < 3; ++i) file_lammps << 0.0 << " "; file_lammps << std::endl; file_lammps << "region box block "; for (UInt i = 0; i < 6; ++i) file_lammps << replica[i] << " "; file_lammps << std::endl; file_lammps << "create_box 1 box" << std::endl; file_lammps << "create_atoms 1 box" << std::endl; file_lammps << std::endl << std::endl; if (triclinic == 1 && floorf(tilt[1] * 100 + 0.5) / 100 != 0.0) { // double value; std::string direction; Real tilting = floorf(tilt[1] * 100 + 0.5) / 100; if (tilt[0] == 0.0) { // value = tilt[1]/lattice_size/lattice_spacing[0]; direction = "xy"; } else if (tilt[0] == 1.0) { // value = tilt[1]/lattice_size/lattice_spacing[0]; direction = "xz"; } else if (tilt[0] == 2.0) { // value = tilt[1]/lattice_size/lattice_spacing[0]; direction = "yz"; } else LM_FATAL("tilt must be happend in xy(0), xz(1) or yz(2) axis!!!"); file_lammps << "change_box triclinic" << std::endl; file_lammps << "displace_box all " << direction << " final " << tilting << " remap none units box" << std::endl; // file_lammps << "displace_box all "< src(in.rdbuf()); std::istreambuf_iterator end; std::ostream_iterator dest(file_lammps); std::copy(src, end, dest); } DUMP("opening lammps config file " << file_to_read, DBG_INFO_STARTUP); infile = fopen(file_to_read.c_str(), "rb"); if (infile == NULL) LM_FATAL("cannot open " << file_to_read << " as a lammps config file !!!"); input->file(); DUMP("lammps file loaded", DBG_INFO); // register lammps computes to libmultiscale system UInt ncompute = this->modify->ncompute; for (UInt i = 0; i < ncompute; ++i) { auto _ptr = std::make_shared>( *(this->modify->compute[i]), *this); _ptr->acquireContext(*this); FilterManager::getManager().addObject(_ptr); } // register lammps computes to libmultiscale system UInt nfix = this->modify->nfix; for (UInt i = 0; i < nfix; ++i) { auto _ptr = std::make_shared>( *(this->modify->fix[i]), *this); _ptr->acquireContext(*this); FilterManager::getManager().addObject(_ptr); } } /* -------------------------------------------------------------------------- */ template void DomainLammps::changeBox() { auto &geom_manager = GeometryManager::getManager(); auto c = geom_manager.getGeometry(this->newGeomBox)->getBoundingBox(); this->setDomainBoundingBox(c); } /* -------------------------------------------------------------------------- */ template void DomainLammps::initModel() { lammps_main_object = this; ref_manager = std::make_shared>(this->atoms); this->atoms.setRefManager(ref_manager); ref_manager->acquireContext(*this); ImportLammpsInterface::createImportLammps(ref_manager); if (this->getGeom().getDim() == 1) LM_FATAL("LAMMPS cannot work in 1D"); ImportLammpsInterface::setGeom(&this->getGeom()); // for (auto &&_type : geom_by_type) { // Geometry &g = GeometryManager::getManager().getObject(_type.second); // UInt itype = _type.first; // ImportLammpsInterface::setGeom(&g, itype); // } loadLammpsConfigFile(); this->atoms.init(); update->whichflag = 0; update->firststep = 0; update->laststep = nb_step; update->ntimestep = 0; update->beginstep = 0; update->endstep = update->laststep; if (flag_units) { force->ftm2v = code_unit_system.ft_m2v; // force->nktv2p = 1; force->mvv2e = code_unit_system.mvv2e; } if (this->domain->dimension != Dim) LM_FATAL("mismatch dimension: check lammps config file"); if (std::string(this->atom->atom_style) != "atomic") LM_FATAL("can only work with atomic style"); this->acceleration_field.resize(atom->nlocal, 3); this->acceleration_field.setZero(); } /* -------------------------------------------------------------------------- */ template void DomainLammps::initDOFs() { { UInt nall; if (force->newton) nall = atom->nlocal + atom->nghost; else nall = atom->nlocal; Real **x = atom->x; Real **x0 = atom->x0; for (UInt i = 0; i < nall; i++) { x0[i][0] = x[i][0]; x0[i][1] = x[i][1]; x0[i][2] = x[i][2]; } } if (this->shouldRestart()) { this->restart_start = true; this->readXMLFile(this->restartfile); if (this->newGeomBox != invalidGeom && current_step >= this->when_change_box) { this->changeBox(); this->enforceCompatibility(); // perform migration after change box to // keep periodicity this->when_change_box = current_step + nb_step; // to avoid repetitions } else if (neighbor->decide()) { this->enforceCompatibility(); this->resetBox(); } this->updateGradient(); LAMMPS_NS::LAMMPS::output->setup(false); } // conversion of masses from g/mol to Kg if (flag_units) { UnitsConverter uc; uc.setOutUnits(code_unit_system); if (strcmp(update->unit_style, "real") == 0) uc.setInUnits(atomic_unit_system); else if (strcmp(update->unit_style, "metal") == 0) uc.setInUnits(metal_unit_system); else LM_FATAL("not known lammps unit system:" << " things need to be done to integrate"); uc.computeConversions(); for (int i = 1; i <= atom->ntypes; i++) atom->mass[i] *= uc.getConversion(); } -#ifdef TRACE_ATOM +#ifdef LIBMULTISCALE_TRACE_ATOM for (int i = 0; i < atom->nlocal; ++i) { - Real X[3] = {atom->x0[i][0], atom->x0[i][1], atom->x0[i][2]}; + Vector<3> X = {atom->x0[i][0], atom->x0[i][1], atom->x0[i][2]}; SET_INTERNAL_TRACE_INDEX(X, i); IF_TRACED(X, "traced atom has been detected at position " << internal_index); } VIEW_ATOM(RefLammps); #endif } /* -------------------------------------------------------------------------- */ template Real DomainLammps::potentialEnergy() { return update->minimize->pe_compute->compute_scalar(); } /* -------------------------------------------------------------------------- */ extern "C" { static UInt fargc = 1; static char prog[20] = "lammps"; static char *___tmp = &prog[0]; static char **fargv = &___tmp; } /* -------------------------------------------------------------------------- */ template DomainLammps::~DomainLammps() {} /* -------------------------------------------------------------------------- */ template DomainLammps::DomainLammps(const LMID &ID, CommGroup &group) : LMObject(ID), DomainAtomic, Dim>(ID, group), LAMMPS(fargc, fargv, group.getMPIComm()), flag_units(true), triclinic(0) { newGeomBox = invalidGeom; when_change_box = 0; lattice = ""; lattice_size = 0.; for (UInt i = 0; i < 6; ++i) replica[i] = 0; create_header_flag = false; restart_start = false; for (UInt i = 0; i < 3; ++i) { lattice_orient_x[i] = 0; lattice_orient_y[i] = 0; lattice_orient_z[i] = 0; lattice_origin[i] = 0.; } lattice_orient_x[0] = 1; lattice_orient_y[1] = 1; lattice_orient_z[2] = 1; lattice_spacing[0] = 1.0; lattice_spacing[1] = 1.0; lattice_spacing[2] = 1.0; } /* -------------------------------------------------------------------------- */ template Cube DomainLammps::getSubDomainBoundingBox() const { Cube c; c.setXmin(VectorView<3>(domain->sublo)); c.setXmax(VectorView<3>(domain->subhi)); return c; } /* -------------------------------------------------------------------------- */ template void DomainLammps::setSubDomainBoundingBox(const Cube &cube) { VectorView<3>(domain->sublo) = cube.getXmin(); VectorView<3>(domain->subhi) = cube.getXmax(); if (neighbor->decide()) this->resetBox(); } /* -------------------------------------------------------------------------- */ template Cube DomainLammps::getDomainBoundingBox() const { Cube c; c.setXmin(VectorView<3>(domain->boxlo)); c.setXmax(VectorView<3>(domain->boxhi)); return c; } /* -------------------------------------------------------------------------- */ template void DomainLammps::setDomainBoundingBox(const Cube &cube) { VectorView<3>(domain->boxlo) = cube.getXmin(); VectorView<3>(domain->boxhi) = cube.getXmax(); if (neighbor->decide()) this->resetBox(); } /* -------------------------------------------------------------------------- */ template void DomainLammps::resetBox() { atom->setup(); if (domain->triclinic) domain->x2lamda(atom->nlocal); domain->pbc(); domain->reset_box(); comm->setup(); if (neighbor->style) neighbor->setup_bins(); comm->exchange(); if (atom->sortfreq > 0) atom->sort(); comm->borders(); if (domain->triclinic) domain->lamda2x(atom->nlocal + atom->nghost); neighbor->build(); } /* -------------------------------------------------------------------------- */ template void DomainLammps::enforceCompatibility() { if (this->getCommGroup().amIinGroup()) { this->echangeAtomes(); this->acceleration_field.resize(atom->nlocal, 3); } } /* -------------------------------------------------------------------------- */ template void DomainLammps::echangeAtomes() { VIEW_ATOM(RefLammps); UInt nflag = neighbor->decide(); // if (nflag == 0 && this->restart_start==false && this->update_pad == false) // { if (nflag == 0 && this->restart_start == false) { timer->stamp(); comm->forward_comm(); timer->stamp(TIME_COMM); } else { if (modify->n_pre_exchange) modify->pre_exchange(); if (domain->triclinic) domain->x2lamda(atom->nlocal); domain->pbc(); if (domain->box_change) { domain->reset_box(); comm->setup(); if (neighbor->style) neighbor->setup_bins(); } timer->stamp(); comm->exchange(); comm->borders(); if (domain->triclinic) domain->lamda2x(atom->nlocal + atom->nghost); timer->stamp(TIME_COMM); if (modify->n_pre_neighbor) modify->pre_neighbor(); neighbor->build(); timer->stamp(TIME_NEIGHBOR); this->restart_start = false; } if (this->when_change_box == current_time) { // if one change box size, this part has to be fulfilled. if (modify->n_pre_exchange) modify->pre_exchange(); if (domain->triclinic) domain->x2lamda(atom->nlocal); domain->pbc(); if (domain->box_change) { domain->reset_box(); comm->setup(); if (neighbor->style) neighbor->setup_bins(); } timer->stamp(); comm->exchange(); comm->borders(); if (domain->triclinic) domain->lamda2x(atom->nlocal + atom->nghost); timer->stamp(TIME_COMM); if (modify->n_pre_neighbor) modify->pre_neighbor(); neighbor->build(); timer->stamp(TIME_NEIGHBOR); this->restart_start = false; } VIEW_ATOM(RefLammps); } /* -------------------------------------------------------------------------- */ template bool DomainLammps::isPeriodic(UInt dir) { switch (dir) { case 0: return this->domain->xperiodic; break; case 1: return this->domain->yperiodic; break; case 2: return this->domain->zperiodic; break; default: LM_FATAL("called on invalid-dimension"); return false; } } /* -------------------------------------------------------------------------- */ template void DomainLammps::force_clear(UInt) { UInt i; // clear global force array // nall includes ghosts only if either newton flag is set UInt nall; if (force->newton) nall = atom->nlocal + atom->nghost; else nall = atom->nlocal; Real **f = atom->f; for (i = 0; i < nall; i++) { f[i][0] = 0.0; f[i][1] = 0.0; f[i][2] = 0.0; } // if (integ.torqueflag) { // Real **torque = atom->torque; // for (i = 0; i < nall; i++) { // torque[i][0] = 0.0; // torque[i][1] = 0.0; // torque[i][2] = 0.0; // } // } // if (granflag) { // Real **phia = atom->phia; // for (i = 0; i < nall; i++) { // phia[i][0] = 0.0; // phia[i][1] = 0.0; // phia[i][2] = 0.0; // } // } // clear pair force array if using it this timestep to compute virial // if (vflag == 2 && pairflag) { // if (atom->nmax > maxpair) { // maxpair = atom->nmax; // memory->destroy_2d_Real_array(f_pair); // f_pair = memory->create_2d_Real_array(maxpair,3,"verlet:f_pair"); // } // for (i = 0; i < nall; i++) { // f_pair[i][0] = 0.0; // f_pair[i][1] = 0.0; // f_pair[i][2] = 0.0; // } // } } /* -------------------------------------------------------------------------- */ template void DomainLammps::setCurrentStep(UInt ts) { current_step = ts; // content of the Update::reset_timestep function update->eflag_global = update->vflag_global = ts; for (int i = 0; i < modify->ncompute; i++) { modify->compute[i]->invoked_scalar = -1; modify->compute[i]->invoked_vector = -1; modify->compute[i]->invoked_array = -1; modify->compute[i]->invoked_peratom = -1; modify->compute[i]->invoked_local = -1; // modify->compute[i]->eflag_global = -1; } update->ntimestep = ts; if (update->laststep < ts) update->laststep += ts; } /* -------------------------------------------------------------------------- */ template void DomainLammps::setTimeStep(Real ts) { update->dt = ts; for (int i = 0; i < modify->nfix; ++i) modify->fix[i]->init(); } /* -------------------------------------------------------------------------- */ template UInt DomainLammps::getCurrentStep() { return update->ntimestep; } /* -------------------------------------------------------------------------- */ template Real DomainLammps::getTimeStep() { return update->dt; } /* -------------------------------------------------------------------------- */ template int DomainLammps::getOrientValue(UInt axis, UInt comp) { if (axis == 0) { return this->lattice_orient_x[comp]; } else if (axis == 1) { return this->lattice_orient_y[comp]; } else if (axis == 2) { return this->lattice_orient_z[comp]; } else { LM_FATAL("wrong axis value: " << axis); } } /* -------------------------------------------------------------------------- */ template Real DomainLammps::getLatticeConstant() { return this->lattice_size; } /* -------------------------------------------------------------------------- */ template void DomainLammps::newmarkPredictor() { if (this->newGeomBox != invalidGeom && current_step >= this->when_change_box) { this->changeBox(); } VIEW_ATOM(RefLammps); STARTTIMER("AtomStepPosition"); this->update->ntimestep++; this->update->integrate->ev_set(this->update->ntimestep); this->modify->initial_integrate(0); this->echangeAtomes(); STOPTIMER("AtomStepPosition"); VIEW_ATOM(RefLammps); this->changeRelease(); } template void DomainLammps::newmarkCorrector() { VIEW_ATOM(RefLammps); DUMP("AtomStepVelocities x " << this->atom->x, DBG_INFO); STARTTIMER("AtomStepVelocities"); this->modify->final_integrate(); if (this->modify->n_end_of_step) this->modify->end_of_step(); if (LAMMPS_NS::LAMMPS::output->next == this->update->ntimestep) { this->timer->stamp(); LAMMPS_NS::LAMMPS::output->write(this->update->ntimestep); this->timer->stamp(TIME_OUTPUT); } STOPTIMER("AtomStepVelocities"); VIEW_ATOM(RefLammps); this->changeRelease(); } /* -------------------------------------------------------------------------- */ template void DomainLammps::updateAcceleration() { if (!this->getCommGroup().amIinGroup()) return; double *rmass = atom->rmass; double *mass = atom->mass; int *type = atom->type; auto f = this->gradient(); std::function compute_mass_inverse; if (rmass) compute_mass_inverse = [&](UInt i) { return code_unit_system.ft_m2v / rmass[i]; }; else compute_mass_inverse = [&](UInt i) { return code_unit_system.ft_m2v / mass[type[i]]; }; for (UInt i = 0; i < UInt(lammps_main_object->atom->nlocal); ++i) { Real mass_inverse = compute_mass_inverse(i); this->acceleration_field(i, 0) = mass_inverse * f(i, 0); this->acceleration_field(i, 1) = mass_inverse * f(i, 1); this->acceleration_field(i, 2) = mass_inverse * f(i, 2); } } /* -------------------------------------------------------------------------- */ template void DomainLammps::updateGradient() { if (!this->getCommGroup().amIinGroup()) return; VIEW_ATOM(RefLammps); STARTTIMER("AtomStepForces"); UInt eflag, vflag; eflag = this->update->integrate->eflag; vflag = this->update->integrate->vflag; this->force_clear(vflag); if (this->modify->n_pre_force) this->modify->pre_force(vflag); this->timer->stamp(); if (this->atom->molecular) { if (this->force->bond) this->force->bond->compute(eflag, vflag); if (this->force->angle) this->force->angle->compute(eflag, vflag); if (this->force->dihedral) this->force->dihedral->compute(eflag, vflag); if (this->force->improper) this->force->improper->compute(eflag, vflag); this->timer->stamp(TIME_BOND); } if (this->force->pair) { this->force->pair->compute(eflag, vflag); this->timer->stamp(TIME_PAIR); } if (this->force->kspace) { this->force->kspace->compute(eflag, vflag); this->timer->stamp(TIME_KSPACE); } if (this->force->newton) { this->comm->reverse_comm(); this->timer->stamp(TIME_COMM); } if (this->modify->n_post_force) this->modify->post_force(vflag); STOPTIMER("AtomStepForces"); VIEW_ATOM(RefLammps); this->changeRelease(); } /* -------------------------------------------------------------------------- */ /* LMDESC LAMMPS_BASE This plugin make the interface with lammps */ /* LMHERITANCE domain_md */ /* LM EXAMPLE Section MultiScale AtomsUnits\\ ...\\ GEOMETRY myGeom CUBE BBOX -1 1 -1 1 -1 1\\ MODEL LAMMPS md\\ ...\\ endSection\\ \\ Section LAMMPS:md AtomsUnits \\ LAMMPS_FILE lammps.in \\ DOMAIN_GEOMETRY 1 \\ endSection */ template void DomainLammps::declareParams() { DomainAtomic, Dim>::declareParams(); /*LMKEYWORD LAMMPS_FILE This specifies the configuration file to be used for lammps. If the create header keyword is used then only the potential definition should be into that file. Otherwise, a complete lammps configuration file can be specified. */ this->parseKeyword("LAMMPS_FILE", lammpsfile); /* LMKEYWORD CREATE_HEADER A header is automatically generated based on the information passed by keywords LATTICE, LATTICE\_SIZE, BOUNDARY and REPLICA */ this->parseTag("CREATE_HEADER", create_header_flag, false); /* LMKEYWORD LATTICE Specifies the lattice to be used to generate the lammps header. Admissible values are lammps keywords (see \url{http://lammps.sandia.gov/doc/lattice.html}) */ this->parseKeyword("LATTICE", lattice, ""); /* LMKEYWORD LATTICE_SIZE Lattice length parameter to generate the primitive lattice */ this->parseKeyword("LATTICE_SIZE", lattice_size, 0.); /* LMKEYWORD LATTICE_ORIGIN Lattice length parameter to generate the primitive lattice */ this->parseVectorKeyword("LATTICE_ORIGIN", 3, lattice_origin, VEC_DEFAULTS(0., 0., 0.)); /* LMKEYWORD LATTICE_ORIENTX Lattice length parameter to generate the primitive lattice */ this->parseVectorKeyword("LATTICE_ORIENTX", 3, lattice_orient_x, VEC_DEFAULTS(1, 0, 0)); /* LMKEYWORD LATTICE_ORIENTY Lattice length parameter to generate the primitive lattice */ this->parseVectorKeyword("LATTICE_ORIENTY", 3, lattice_orient_y, VEC_DEFAULTS(0, 1, 0)); /* LMKEYWORD LATTICE_ORIENTZ Lattice length parameter to generate the primitive lattice */ this->parseVectorKeyword("LATTICE_ORIENTZ", 3, lattice_orient_z, VEC_DEFAULTS(0, 0, 1)); /* LMKEYWORD SPACING The minimal normalized lattice sizes */ this->parseVectorKeyword("SPACING", 3, lattice_spacing, VEC_DEFAULTS(1., 1., 1.)); /* LMKEYWORD BOUNDARY Sequence of lammps code for boundary. (see \url{http://lammps.sandia.gov/doc/boundary.html}) */ this->parseVectorKeyword("BOUNDARY", 3, boundaries, VEC_DEFAULTS("", "", "")); /* LMKEYWORD REPLICA For the creation of atoms a region, specified as the number of replicas of the lattice must be provided, before calling the create atoms command. */ this->parseVectorKeyword("REPLICA", 6, replica, VEC_DEFAULTS(0, 0, 0, 0, 0, 0)); /* LMKEYWORD CHANGE_DOMAIN_BOX This allows to specify live reshaping of the bounding box. Useful to enforce stresses states or dislocation displacement field with periodic boundary conditions. */ this->parseKeyword("CHANGE_DOMAIN_BOX", newGeomBox, invalidGeom); /* LMKEYWORD CHANGE_DOMAIN_BOX_TIME Gives a timestep where the reshaping provided by CHANGE\_DOMAIN\_BOX should occur. */ this->parseKeyword("CHANGE_DOMAIN_BOX_TIME", when_change_box, 0u); /* LMKEYWORD TRICLINIC This changes the orthogonal simulation box to triclinic box */ this->parseKeyword("TRICLINIC", triclinic, 0); /* LMKEYWORD TILT This tils the simulation box by given amount value 0: xy-direction; 1: xz-direction; 2: yz-direction (TILT 0 0.25: tilt simulation box in xy-direction by 0.25) */ this->parseVectorKeyword("TILT", 2, tilt, VEC_DEFAULTS(0., 0.)); /* LMKEYWORD UNITS_CONVERSION This is for debug purpose: force libmultiscale not to touch the units conversion factors of lammps. */ this->parseTag("UNITS_CONVERSION", flag_units, true); } /* -------------------------------------------------------------------------- */ template class DomainLammps<2>; template class DomainLammps<3>; __END_LIBMULTISCALE__ diff --git a/src/md/lammps/import_lammps.cc b/src/md/lammps/import_lammps.cc index 1e0ac1a..6dbdd5a 100644 --- a/src/md/lammps/import_lammps.cc +++ b/src/md/lammps/import_lammps.cc @@ -1,741 +1,741 @@ /** * @file import_lammps.cc * * @author Guillaume Anciaux * @author Jaehyun Cho * * @date Wed Nov 20 07:29:13 2013 * * @brief The programming hooks to let LM be aware of processor migrations * * @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 "lm_common.hh" #ifdef LIBMULTISCALE_USE_LAMMPS #include "domain_lammps.hh" #include "geometry.hh" #include "reference_manager_lammps.hh" #include "trace_atom.hh" /* -------------------------------------------------------------------------- */ using namespace LAMMPS_NS; #include "atom_vec.h" #include "atom_vec_atomic.h" #include "domain.h" #include "fix.h" #include "lammps_common.hh" #include "memory.h" /* -------------------------------------------------------------------------- */ #include "import_lammps.hh" /* -------------------------------------------------------------------------- */ enum Topology { FACEXLEFT, FACEXRIGHT, FACEYLEFT, FACEYRIGHT, FACEZLEFT, FACEZRIGHT, EDGE1, EDGE2, EDGE3, EDGE4, EDGE5, EDGE6, EDGE7, EDGE8, EDGE9, EDGE10, EDGE11, EDGE12, CORNER1, CORNER2, CORNER3, CORNER4, CORNER5, CORNER6, CORNER7, CORNER8 }; /* -------------------------------------------------------------------------- */ static const int n_neighs = 26; // static const int commPairs[n_neighs][2] = { // // faces pairs // {FACEXLEFT, FACEXRIGHT}, // {FACEXRIGHT, FACEXLEFT}, // {FACEYLEFT, FACEYRIGHT}, // {FACEYRIGHT, FACEYLEFT}, // {FACEZLEFT, FACEZRIGHT}, // {FACEZRIGHT, FACEZLEFT}, // // edges pairs // {EDGE1, EDGE2}, // {EDGE2, EDGE1}, // {EDGE3, EDGE4}, // {EDGE4, EDGE3}, // {EDGE5, EDGE6}, // {EDGE6, EDGE5}, // {EDGE7, EDGE8}, // {EDGE8, EDGE7}, // {EDGE9, EDGE10}, // {EDGE10, EDGE9}, // {EDGE11, EDGE12}, // {EDGE12, EDGE11}, // // corner pairs // {CORNER1, CORNER2}, // {CORNER2, CORNER1}, // {CORNER3, CORNER4}, // {CORNER4, CORNER3}, // {CORNER5, CORNER6}, // {CORNER6, CORNER5}, // {CORNER7, CORNER8}, // {CORNER8, CORNER7}}; /* -------------------------------------------------------------------------- */ int topologyCoords[n_neighs][3] = { // faces {-1, 0, 0}, // FACEXLEFT {1, 0, 0}, // FACEXRIGHT {0, -1, 0}, // FACEYLEFT {0, 1, 0}, // FACEYRIGHT {0, 0, -1}, // FACEZLEFT {0, 0, 1}, // FACEZRIGHT // edges {1, 0, -1}, // EDGE1 {-1, 0, -1}, // EDGE2 {1, 0, 1}, // EDGE3 {-1, 0, 1}, // EDGE4 {-1, -1, 0}, // EDGE5 {1, -1, 0}, // EDGE6 {-1, 1, 0}, // EDGE7 {1, 1, 0}, // EDGE8 {0, -1, -1}, // EDGE9 {0, 1, -1}, // EDGE1 {0, -1, 1}, // EDGE1 {0, 1, 1}, // EDGE1 // corners {-1, -1, -1}, // CORNER1 {-1, 1, -1}, // CORNER2 {-1, 1, 1}, // CORNER3 {-1, -1, 1}, // CORNER4 {1, -1, -1}, // CORNER5 {1, 1, -1}, // CORNER6 {1, 1, 1}, // CORNER7 {1, -1, 1} // CORNER8 }; const char topologyStrings[n_neighs][255] = { "FACEXLEFT", "FACEXRIGHT", "FACEYLEFT", "FACEYRIGHT", "FACEZLEFT", "FACEZRIGHT", "EDGE1", "EDGE2", "EDGE3", "EDGE4", "EDGE5", "EDGE6", "EDGE7", "EDGE8", "EDGE9", "EDGE10", "EDGE11", "EDGE12", "CORNER1", "CORNER2", "CORNER3", "CORNER4", "CORNER5", "CORNER6", "CORNER7", "CORNER8"}; /* -------------------------------------------------------------------------- */ #define stamp10 1000 #define stamp20 2000 /* -------------------------------------------------------------------------- */ __BEGIN_LIBMULTISCALE__ /* -------------------------------------------------------------------------- */ void checkMPIComm(MPI_Status &status, int flag) { if (flag != MPI_SUCCESS || status.MPI_ERROR != MPI_SUCCESS) { char errStr[MPI_MAX_ERROR_STRING]; int len; MPI_Error_string(status.MPI_ERROR, errStr, &len); LM_FATAL("MPI_ERROR = " << status.MPI_ERROR << "[" << errStr); } } /* -------------------------------------------------------------------------- */ template inline void ImportLammps::init(LAMMPS_NS::Domain &domain, LAMMPS_NS::Atom &atom, LAMMPS_NS::Comm &comm, MPI_Comm &world, int me, int nprocs) { old_nlocal = atom.nlocal; VIEW_ATOM(RefLammps); this->atom = &atom; this->domain = &domain; this->comm = &comm; this->world = world; this->me = me; this->nprocs = nprocs; setPeriodic(domain.xperiodic, domain.yperiodic, domain.zperiodic); setupTopology(); MPI_Barrier(world); DUMP("before exchange nlocal = " << atom.nlocal, DBG_INFO); if (domain.triclinic == 0) { Xlo = domain.sublo; Xhi = domain.subhi; } else { Xlo = domain.sublo_lamda; Xhi = domain.subhi_lamda; } } /* -------------------------------------------------------------------------- */ template inline void ImportLammps::registerSends() { // clear the buffers for (UInt i = 0; i < buf_send.size(); ++i) buf_send[i].clear(); int nlocal = atom->nlocal; int proc[3]; for (UInt i = 0; i < 3; ++i) DUMP("box dimension " << i << " " << Xlo[i] << " " << Xhi[i], DBG_INFO); DUMP(" nlocal = " << nlocal, DBG_INFO); int at = 0; DUMP("doing loop over atoms for migration", DBG_INFO); while (at < nlocal) { Real *X = atom->x[at]; #ifndef LM_OPTIMIZED Real *X0 = atom->x0[at]; Real *V = atom->v[at]; Real *F = atom->f[at]; #endif // setting actual coordinates for (UInt k = 0; k < 3; ++k) { proc[k] = 0; // search the new processor if (X[k] < Xlo[k]) proc[k] = -1; if (X[k] >= Xhi[k]) proc[k] = 1; // special treatment for pbc if (periodic[k]) { if (domain->triclinic == 0) { // test if i am at the lower border if (Xlo[k] == domain->boxlo[k]) if (X[k] >= Xhi[k] && (X[k] - Xhi[k]) > (Xhi[k] - Xlo[k])) proc[k] = -1; // test if i am at the upper border if (Xhi[k] == domain->boxhi[k]) if (X[k] < Xlo[k] && (Xlo[k] - X[k]) > (Xhi[k] - Xlo[k])) proc[k] = 1; } else { // test if i am at the lower border if (Xlo[k] == domain->boxlo_lamda[k]) if (X[k] >= Xhi[k] && (X[k] - Xhi[k]) > (Xhi[k] - Xlo[k])) proc[k] = -1; // test if i am at the upper border if (Xhi[k] == domain->boxhi_lamda[k]) if (X[k] < Xlo[k] && (Xlo[k] - X[k]) > (Xhi[k] - Xlo[k])) proc[k] = 1; } } } // DUMP("at: "<avec->copy(nlocal - 1, at); // moveAtom(nlocal-1,at); // decrease total number of atoms nlocal--; } atom->nlocal = nlocal; DUMP("nlocal is now " << atom->nlocal, DBG_INFO); } /* -------------------------------------------------------------------------- */ template inline void ImportLammps::registerSend(UInt at, UInt toProc) { #ifdef TRACE_ATOM CHECK_TRACED(atom->x0[at]) { VIEW_ATOM(RefLammps); DUMP(ATOM_OUT(atom->x0[at]) << " atom is leaving my proc for proc " << toProc, DBG_MESSAGE); DUMP(ATOM_OUT(atom->x0[at]) << " atom is leaving my proc for proc " << toProc, DBG_WARNING); internal_index = -1; } #endif sendAtom(at, toProc); UInt nsend = buf_send[toProc].size(); buf_send[toProc].resize(nsend + packSize); DUMP("buf_send size is now " << buf_send[toProc].size(), DBG_DETAIL); UInt added = atom->avec->pack_exchange(at, &buf_send[toProc][nsend]); if (added != packSize) LM_FATAL("internal error"); } /* -------------------------------------------------------------------------- */ template inline void ImportLammps::unPackAtoms(UInt procDe) { UInt m = 0; std::vector &buf = buf_recv[procDe]; DUMP("unpacking atom " << buf_recv[procDe].size() / packSize << " atoms", DBG_INFO); while (m < buf_recv[procDe].size()) { m += atom->avec->unpack_exchange(&buf[m]); acceptAtom(atom->nlocal - 1, procDe); { UInt i = atom->nlocal - 1; Real **ptr = atom->x; Real **x0ptr = atom->x0; Real **vptr = atom->v; SET_INTERNAL_TRACE_INDEX(x0ptr[i], i); - CHECK_TRACED(x0ptr[i]) { + if (CHECK_TRACED(x0ptr[i])) { DUMP("received traced atom from " << procDe, DBG_MESSAGE); VIEW_ATOM(RefLammps); } if ((ptr[i][0] < Xlo[0] || ptr[i][0] >= Xhi[0]) || ptr[i][1] < Xlo[1] || ptr[i][1] >= Xhi[1] || ptr[i][2] < Xlo[2] || ptr[i][2] >= Xhi[2]) { LM_FATAL("pb severe in migration process (received from " << procDe << " : atom i " << i << " x: " << ptr[i][0] << " xlo: " << Xlo[0] << " xhi: " << Xhi[0] << " " << " y: " << ptr[i][1] << " ylo: " << Xlo[1] << " yhi: " << Xhi[1] << " " << " z: " << ptr[i][2] << " zlo: " << Xlo[2] << " zhi: " << Xhi[2] << " " << " big trouble an atom went far far away cannot " << "control anything now : its position is " << x0ptr[i][0] << " " << x0ptr[i][1] << " " << x0ptr[i][2] << " and its velocity is " << vptr[i][0] << " " << vptr[i][1] << " " << vptr[i][2]); } } } if (m != buf_recv[procDe].size()) LM_FATAL("not all atoms where unpacked"); buf_recv[procDe].clear(); } /* -------------------------------------------------------------------------- */ template inline void ImportLammps::checkConsistency() { MPI_Barrier(world); int test = 0; MPI_Allreduce(&(atom->nlocal), &test, 1, MPI_INT, MPI_SUM, world); DUMP("now nlocal = " << atom->nlocal << " nmax = " << atom->nmax, DBG_INFO); if (test != atom->natoms) { // ref_manager.printBilan(); LM_FATAL("atoms were lost/created " << test << " != " << atom->natoms << " " << test - atom->natoms << std::endl << "# local atoms is " << atom->nlocal << std::endl << "old # local atoms is " << old_nlocal << std::endl << "abort"); } if (atom->nmax < atom->nlocal) LM_FATAL("pb with allocation of the atoms"); } /* -------------------------------------------------------------------------- */ template inline void ImportLammps::atomCommunications() { //! pack the atoms needing a migration registerSends(); //! exchange the atoms atomCommunication(); //! unpack the data UInt ncom = proc_comm.size(); for (UInt i = 0; i < ncom; ++i) unPackAtoms(proc_comm[i]); // check constistency checkConsistency(); } /* -------------------------------------------------------------------------- */ template inline void ImportLammps::atomCommunication() { std::map requests; UInt ncom = proc_comm.size(); //! asynchronous send phase for (UInt i = 0; i < ncom; ++i) { UInt toProc = proc_comm[i]; if (int(toProc) == me) continue; int nbSend = buf_send[toProc].size(); DUMP("Send to proc " << toProc << " the " << nbSend << " Reals t:" << stamp20 + me, DBG_INFO); // if (!nbSend) continue; if (requests.count(toProc)) LM_FATAL("request to " << toProc << " already made"); MPI_Request &req = requests[toProc]; MPI_Isend(&buf_send[toProc][0], nbSend, MPI_DOUBLE, toProc, stamp20 + me, world, &req); } //! blocking receive phase for (UInt i = 0; i < ncom; ++i) { UInt fromProc = proc_comm[i]; if (int(fromProc) == me) continue; MPI_Status status; status.MPI_ERROR = MPI_SUCCESS; DUMP("probing from proc " << fromProc, DBG_INFO); MPI_Probe(fromProc, stamp20 + fromProc, world, &status); int nbRecv = 0; MPI_Get_count(&status, MPI_DOUBLE, &nbRecv); buf_recv[fromProc].resize(nbRecv); DUMP("Receive from proc " << fromProc << " the " << nbRecv << " Reals t:" << stamp20 + fromProc, DBG_INFO); // if (!nbRecv) continue; int ret = MPI_Recv(&buf_recv[fromProc][0], nbRecv, MPI_DOUBLE, fromProc, stamp20 + fromProc, world, &status); checkMPIComm(status, ret); } //! wait for sends to be complete std::map::iterator it = requests.begin(); std::map::iterator end = requests.end(); while (it != end) { MPI_Status status; status.MPI_ERROR = MPI_SUCCESS; MPI_Request req = it->second; int toProc = it->first; DUMP("wait for com to " << toProc << " to be finished", DBG_INFO); int ret = MPI_Wait(&req, &status); if (buf_send[toProc].size()) { checkMPIComm(status, ret); } DUMP("send to proc " << toProc << " finished ", DBG_INFO); ++it; } } /* -------------------------------------------------------------------------- */ template int ImportLammps::getProcNumber(int ind[3], int offset[3]) { /* Return the number of the processor located at the position i[0],i[1],i[2] */ /* It returns MPI_PROC_NULL if none */ /* It hanles periodic topology */ int index[3]; for (UInt i = 0; i < 3; ++i) index[i] = ind[i] + offset[i]; DUMP("ind " << ind[0] << " " << ind[1] << " " << ind[2], DBG_DETAIL); DUMP("offset " << offset[0] << " " << offset[1] << " " << offset[2], DBG_DETAIL); DUMP("periodic " << periodic[0] << " " << periodic[1] << " " << periodic[2], DBG_DETAIL); DUMP("nbrProc " << nbrProc[0] << " " << nbrProc[1] << " " << nbrProc[2], DBG_DETAIL); DUMP("index " << index[0] << " " << index[1] << " " << index[2], DBG_DETAIL); // exclude out of the domain faces (negative side) for (UInt i = 0; i < 3; ++i) if (index[i] < 0 && !periodic[i]) return MPI_PROC_NULL; else if (index[i] < 0) index[i] = nbrProc[i] - 1; // exclude out of the domain faces (positive side) for (UInt i = 0; i < 3; ++i) if (index[i] >= nbrProc[i] && !periodic[i]) return MPI_PROC_NULL; else if (index[i] >= nbrProc[i]) index[i] = 0; DUMP("index " << index[0] << " " << index[1] << " " << index[2], DBG_DETAIL); return proc_matrix[index]; } /* -------------------------------------------------------------------------- */ //! setting processor topology with a stamp style template void ImportLammps::setupTopology() { int reorder = 0; // get number of processors in x,y and z directions for (UInt i = 0; i < 3; ++i) nbrProc[i] = comm->procgrid[i]; // MPI communicator in cartesian coordinates MPI_Comm cartesian; // creation of the topology MPI_Cart_create(world, 3, nbrProc, periodic, reorder, &cartesian); // fill the processor repartition matrix proc_matrix.resize(nbrProc); int coords[3]; for (int n = 0; n < nprocs; n++) { MPI_Cart_coords(cartesian, n, 3, coords); proc_matrix[coords] = n; DUMP("proc " << n << " has cartesian coordinates : " << coords[0] << " " << coords[1] << " " << coords[2], DBG_INFO); } /* coordonnes du processeur */ MPI_Cart_coords(cartesian, me, 3, coords); DUMP("my coordinates are " << coords[0] << " " << coords[1] << " " << coords[2], DBG_INFO); MPI_Comm_free(&cartesian); // initialize comprocs for (int i = 0; i < n_neighs; ++i) { DUMP("get " << topologyStrings[i] << " proc number", DBG_DETAIL); comProcs[i] = getProcNumber(coords, topologyCoords[i]); DUMP(topologyStrings[i] << " " << comProcs[i], DBG_INFO); } for (UInt i = 0; i < 3; ++i) DUMP("nbrProc[" << i << "]" << nbrProc[i], DBG_INFO); // initialize com buffers array buf_send.resize(nprocs); buf_recv.resize(nprocs); // init list of processors to discuss with std::set proc_comm_with; for (int i = 0; i < n_neighs; ++i) proc_comm_with.insert(comProcs[i]); proc_comm.clear(); std::set::iterator it = proc_comm_with.begin(); std::set::iterator end = proc_comm_with.end(); while (it != end) { if (*it != MPI_PROC_NULL) { proc_comm.push_back(*it); DUMP("I should communicate with proc " << *it, DBG_DETAIL); } ++it; } DUMP("I should communicate with " << proc_comm.size() << " processors", DBG_DETAIL); // compute the packSize double tmp[1000]; packSize = atom->avec->pack_exchange(0, tmp); DUMP("packSize is " << packSize, DBG_INFO); } /* -------------------------------------------------------------------------- */ template class ImportLammps<2>; template class ImportLammps<3>; ImportLammpsInterface *ImportLammpsInterface::import = NULL; /* -------------------------------------------------------------------------- */ __END_LIBMULTISCALE__ using namespace libmultiscale; /* ---------------------------------------------------------------------- */ void AtomVecAtomic::copy(int i, int j) { tag[j] = tag[i]; type[j] = type[i]; mask[j] = mask[i]; image[j] = image[i]; x[j][0] = x[i][0]; x[j][1] = x[i][1]; x[j][2] = x[i][2]; x0[j][0] = x0[i][0]; x0[j][1] = x0[i][1]; x0[j][2] = x0[i][2]; v[j][0] = v[i][0]; v[j][1] = v[i][1]; v[j][2] = v[i][2]; if (atom->nextra_grow) for (int iextra = 0; iextra < atom->nextra_grow; iextra++) modify->fix[atom->extra_grow[iextra]]->copy_arrays(i, j); ImportLammpsInterface::getImport().moveAtom(i, j); } /* -------------------------------------------------------------------------- */ void Comm::exchange() { MPI_Barrier(world); DUMP("migration coms", DBG_INFO); ImportLammpsInterface::getImport().init(*domain, *atom, *this, world, this->me, this->nprocs); // clear global->local map since atoms move & new ghosts are created if (map_style) atom->map_clear(); DUMP("after construction of sending buffers nlocal = " << atom->nlocal << " nmax = " << atom->nmax, DBG_INFO); // send/receive communication to other procs ImportLammpsInterface::getImport().atomCommunications(); MPI_Barrier(world); ImportLammpsInterface::getImport().updateRefSubSets(); MPI_Barrier(world); DUMP("migration all done for me !", DBG_INFO); } /* -------------------------------------------------------------------------- */ // cleaned /* -------------------------------------------------------------------------- */ void AtomVecAtomic::create_atom(int itype, Real *coord) { // add the atom to my list of atoms if (!ImportLammpsInterface::getImport().getGeom(itype).contains( coord[0], coord[1], coord[2])) return; int nlocal = atom->nlocal; if (nlocal == nmax) grow(0); tag[nlocal] = 0; type[nlocal] = itype; for (UInt i = 0; i < 3; ++i) { x[nlocal][i] = coord[i]; x0[nlocal][i] = coord[i]; v[nlocal][i] = 0.0; } mask[nlocal] = 1; image[nlocal] = (512 << 20) | (512 << 10) | 512; Real *X = x[nlocal]; Real *xlo = domain->sublo; Real *xhi = domain->subhi; for (UInt i = 0; i < 3; ++i) { if (X[i] < xlo[i]) { LM_FATAL("cannot accept to create initially an atom outside" << " of the initial box" << "X[" << i << "] = " << X[i] << " < xlo = " << xlo[i]); } if (X[i] >= xhi[i]) { LM_FATAL("cannot accept to create initially an atom outside" << " of the initial box" << "X[" << i << "] = " << X[i] << " >= xhi = " << xhi[i]); } } atom->nlocal++; } /* -------------------------------------------------------------------------- */ #endif diff --git a/src/md/lammps/ref_atom_lammps.hh b/src/md/lammps/ref_atom_lammps.hh index b764e7c..134c28e 100644 --- a/src/md/lammps/ref_atom_lammps.hh +++ b/src/md/lammps/ref_atom_lammps.hh @@ -1,289 +1,291 @@ /** * @file ref_atom_lammps.hh * * @author Guillaume Anciaux * * @date Mon Oct 28 19:23:14 2013 * * @brief Reference over LAMMPS atom structure * * @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. * */ #ifndef __LIBMULTISCALE_REF_ATOM_LAMMPS_HH__ #define __LIBMULTISCALE_REF_ATOM_LAMMPS_HH__ /* -------------------------------------------------------------------------- */ #include "epot_hook.hh" #include "pair_eam.h" #include "ref_atom.hh" #include "stress_hook.hh" #include #include /* -------------------------------------------------------------------------- */ #include "atom.h" /* -------------------------------------------------------------------------- */ extern ::LAMMPS_NS::LAMMPS *lammps_main_object; /* -------------------------------------------------------------------------- */ __BEGIN_LIBMULTISCALE__ template class ReferenceManagerLammps; template class ContainerLammps; /* -------------------------------------------------------------------------- */ class ComparatorRefLammps; template class DomainLammps; /** * Class RefLammps * */ template class RefLammps : public RefAtom<_Dim, RefLammps<_Dim>> { public: static constexpr UInt Dim = _Dim; using RefComparator = ComparatorRefLammps; using Domain = DomainLammps<_Dim>; using fields = field_list<_position0, _position, _displacement, _velocity, _force, _stress, _mass>; RefLammps(UInt id) : index_atom(id), _weight(1), stress_hook(NULL), stressKin_hook(NULL), epot_hook(NULL){}; RefLammps() { RefLammps(0); } virtual ~RefLammps(){}; public: inline VectorView<_Dim> position0(); inline VectorView<_Dim> position(); inline AccessorAtomDof, _displacement> displacement(); inline VectorView<_Dim> velocity(); inline VectorView<_Dim> force(); inline MatrixView<_Dim> stress(); inline MatrixView<_Dim> kinetic_stress(); inline Real electronic_density(); inline Real getEPot(); inline UInt typeID() { return lammps_main_object->atom->type[index_atom]; }; inline UInt id(); inline Real &mass(); inline void setMass(Real mass); inline Real &charge(); inline Real &alpha(); inline UInt tag(); inline int getIndex() const { return index_atom; }; inline bool operator==(const RefLammps &ref) const { return (ref.index_atom == index_atom); }; inline void setIndex(UInt n) { index_atom = n; }; friend class ComparatorRefLammps; friend class ContainerLammps; friend class ReferenceManagerLammps; private: int index_atom; Real _weight; StressHookLammps *stress_hook; StressHookLammps *stressKin_hook; EpotHookLammps *epot_hook; }; /* -------------------------------------------------------------------------- */ template inline Real RefLammps<_Dim>::electronic_density() { #ifdef PATCHED_LAMMPS CHECK_INDEX; LM_ASSERT(dynamic_cast(lammps_main_object->force->pair), "the potential has to be EAM to get electronic density"); if (current_step == 0) return 0; LAMMPS_NS::PairEAM *ptr = dynamic_cast(lammps_main_object->force->pair); return ptr->rho[index_atom]; #else LM_FATAL("to use this option you have to patch lammps to let 'rho'" << " array within EAM strucuture to be public"); #endif return 0; } /* -------------------------------------------------------------------------- */ template inline VectorView RefLammps::position0() { return VectorView(lammps_main_object->atom->x0[index_atom]); } /* -------------------------------------------------------------------------- */ template inline AccessorAtomDof, _displacement> RefLammps::displacement() { return AccessorAtomDof, _displacement>(*this); } /* -------------------------------------------------------------------------- */ template inline VectorView RefLammps::position() { return VectorView(lammps_main_object->atom->x[index_atom]); } /* -------------------------------------------------------------------------- */ template inline UInt RefLammps::id() { return lammps_main_object->atom->type[index_atom]; } /* -------------------------------------------------------------------------- */ template inline VectorView RefLammps::velocity() { return VectorView(lammps_main_object->atom->v[index_atom]); } /* -------------------------------------------------------------------------- */ template inline VectorView RefLammps::force() { return VectorView(lammps_main_object->atom->f[index_atom]); } /* -------------------------------------------------------------------------- */ template inline UInt RefLammps::tag() { return lammps_main_object->atom->tag[index_atom]; } /* -------------------------------------------------------------------------- */ template inline MatrixView RefLammps::stress() { LM_TOIMPLEMENT; // LM_ASSERT(stress_hook, // "stress hook not initialized : stress cannot be computed"); // if (!stress_hook->Computed()) // stress_hook->ComputeStress(); // stress_hook->stress(index_atom, i); // return stress_hook->stress(index_atom, i); } /* -------------------------------------------------------------------------- */ template inline MatrixView RefLammps::kinetic_stress() { LM_TOIMPLEMENT; // LM_ASSERT(stressKin_hook, "stress hook(Kin) not initialized :" // << " stress cannot be computed"); // if (!stressKin_hook->Computed()) // stressKin_hook->ComputeStress(); // stressKin_hook->stress(index_atom, i); // return stressKin_hook->stress(index_atom, i); } /* -------------------------------------------------------------------------- */ template inline Real RefLammps::getEPot() { if (!epot_hook) { LM_ASSERT(epot_hook, "epot hook not initialized : epot cannot be computed"); } if (!epot_hook->isComputed()) epot_hook->computeEpot(); return epot_hook->epot(index_atom); } /* -------------------------------------------------------------------------- */ template inline Real &RefLammps::charge() { return lammps_main_object->atom ->dipole[lammps_main_object->atom->type[index_atom]]; } /* -------------------------------------------------------------------------- */ template inline Real &RefLammps::alpha() { return _weight; } /* -------------------------------------------------------------------------- */ template inline Real &RefLammps::mass() { return lammps_main_object->atom ->mass[lammps_main_object->atom->type[index_atom]]; } /* -------------------------------------------------------------------------- */ template inline void RefLammps::setMass(Real mass) { lammps_main_object->atom->mass[lammps_main_object->atom->type[index_atom]] = mass; } /* -------------------------------------------------------------------------- */ class ComparatorRefLammps { /* ------------------------------------------------------------------------ */ /* Methods */ /* ------------------------------------------------------------------------ */ public: template bool operator()(RefLammps &r1, RefLammps &r2) const { return r1.index_atom < r2.index_atom; } template bool operator()(const RefLammps &r1, const RefLammps &r2) const { return r1.index_atom < r2.index_atom; } }; /* -------------------------------------------------------------------------- */ template inline std::ostream &operator<<(std::ostream &os, RefLammps &ref) { UInt n = ref.getIndex(); - os << "Lammps reference : " << n << " [P] " << ref.position() << " [P0] " - << ref.position0() << " [V] " << ref.velocity() << " [F] " << ref.force(); + Eigen::IOFormat fmt(Eigen::FullPrecision, 0, ", ", "", "", "", "[", "]"); + os << "Lammps reference : " << n << " P=" << ref.position().format(fmt) + << " P0=" << ref.position0().format(fmt) + << " V=" << ref.velocity().format(fmt) << " F=" << ref.force().format(fmt); return os; } /* -------------------------------------------------------------------------- */ template inline std::ostream &operator<<(std::ostream &os, const RefLammps &ref) { LM_TOIMPLEMENT; return os; } /* -------------------------------------------------------------------------- */ __END_LIBMULTISCALE__ #endif /* __LIBMULTISCALE_REF_ATOM_LAMMPS_HH__ */ diff --git a/src/md/lammps/reference_manager_lammps.cc b/src/md/lammps/reference_manager_lammps.cc index 8d92ab3..d9780d2 100644 --- a/src/md/lammps/reference_manager_lammps.cc +++ b/src/md/lammps/reference_manager_lammps.cc @@ -1,198 +1,198 @@ /** * @file reference_manager_lammps.cc * * @author Guillaume Anciaux * * @date Wed Mar 13 16:53:59 2013 * * @brief This is the reference manager of the LAMMPS model which handles * migrations * * @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 "reference_manager_lammps.hh" #include "domain_lammps.hh" #include "lm_common.hh" #include "trace_atom.hh" /* -------------------------------------------------------------------------- */ __BEGIN_LIBMULTISCALE__ template ReferenceManagerLammps::ReferenceManagerLammps( ContainerLammps &global_container) : ReferenceManager>(global_container) {} /* -------------------------------------------------------------------------- */ template void ReferenceManagerLammps::sendAtom(UInt i, UInt proc) { if (!this->flag_need_check_movements) return; RefLammps refold(i); - IF_REF_TRACED(refold, "sending atom to proc " << proc); + IF_REF_TRACED(refold, "sending atom(" << i << ") to proc " << proc); SET_REF_INTERNAL_TRACE_INDEX(refold, -1); for (auto [temp, temp2] : this->moved) { if (refold == temp2) { DUMP("The atom to be sent was already moved away from ref" << temp << " vers " << temp2, DBG_DETAIL); refold = temp; this->moved.erase(temp); LM_ASSERT(this->moved.count(temp) == 0, "pb avec la STL"); LM_ASSERT(refold.index_atom == temp.index_atom, "pb avec la STL"); break; } } for (auto [temp, proc2] : this->sent) { if (refold == temp) { LM_FATAL( "the atom " << refold << " that is requested to be sent was already sent to " << proc2 << " and now it is to be sent to " << proc); } } this->sent[refold] = proc; this->sent_byproc[proc].push_back(refold); this->have_changed = true; } /* -------------------------------------------------------------------------- */ template void ReferenceManagerLammps::moveAtom(UInt i, UInt j) { if (!this->flag_need_check_movements) return; // pas la peine de faire un traitement si le deplacement ne fait rien if (i == j) { return; } RefLammps refold(i); RefLammps refnew(j); IF_REF_TRACED(refnew, "atom moved in memory from " << i << " to " << j << " index"); SET_REF_INTERNAL_TRACE_INDEX(refnew, j); DUMP("I move atom " << refold << " to " << refnew, DBG_INFO); LM_ASSERT(!isNaN(refnew.position0()), "l'atome de reference " << refnew << " fait de la merde : pb de deplacement des atomes"); for (auto [proc, v] : this->newrefs) { for (UInt l = 0; l < v.size(); ++l) { DUMP("proc = " << proc << " v[" << l << "]=" << v[l], DBG_ALL); LM_ASSERT(!(v[l] == refold), "houston on a un probleme " << v[l] << " was previously received from " << proc); } } for (auto &&[temp, temp2] : this->moved) { if (refold == temp2) { DUMP("previously moved from " << temp << " to " << temp2 << " : switching source from " << refold << " to " << temp, DBG_DETAIL); this->moved[temp] = refnew; this->have_changed = true; return; } } this->moved[refold] = refnew; this->have_changed = true; } /* -------------------------------------------------------------------------- */ template void ReferenceManagerLammps::acceptAtom(UInt i, UInt fromProc) { if (!this->flag_need_check_movements) return; RefLammps newref(i); IF_REF_TRACED(newref, "atom received from proc " << fromProc); SET_REF_INTERNAL_TRACE_INDEX(newref, i); this->newrefs[fromProc].push_back(newref); DUMP("I accept atom " << newref << " from " << fromProc, DBG_ALL); LM_ASSERT(!isNaN(newref.position0()), "l'atome de reference " << newref << "fait de la merde : pb dans la gestion des migrations " << "pour atome recu de " << fromProc); for (auto &&it : this->moved) { RefLammps temp = it.first; RefLammps temp2 = it.second; LM_ASSERT(!(newref == temp2), "recois un atome a un emplacement deja utilise " << "par un deplacement depuis " << temp << " vers " << temp2); } LM_ASSERT(!isNaN(newref.position0()), "l'atome de reference " << newref << "fait de la merde : pb dans la gestion des migrations " << "pour atome recu de " << fromProc); this->have_changed = true; } /* -------------------------------------------------------------------------- */ template class ReferenceManagerLammps<2>; template class ReferenceManagerLammps<3>; __END_LIBMULTISCALE__ diff --git a/src/md/trace_atom.cc b/src/md/trace_atom.cc index 9699dc2..1df2710 100644 --- a/src/md/trace_atom.cc +++ b/src/md/trace_atom.cc @@ -1,47 +1,47 @@ /** * @file trace_atom.cc * * @author Guillaume Anciaux * * @date Mon Sep 30 12:13:51 2013 * * @brief Set of macros to trace an atom through processor migrations * * @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 . * */ -#ifdef TRACE_ATOM +#ifdef LIBMULTISCALE_TRACE_ATOM #include "lm_common.hh" #include /* -------------------------------------------------------------------------- */ __BEGIN_LIBMULTISCALE__ typedef unsigned int UInt; typedef double Real; -Real trace_ref[3] = {-51.0966, -2.26547, -8.5241}; +Real trace_ref[3] = {10.4832, 16.5312, 38.7072}; Real trace_tolerance = 0.01; UInt internal_index = -1; __END_LIBMULTISCALE__ #endif diff --git a/src/md/trace_atom.hh b/src/md/trace_atom.hh index 4145bde..7e99216 100644 --- a/src/md/trace_atom.hh +++ b/src/md/trace_atom.hh @@ -1,107 +1,115 @@ /** * @file trace_atom.hh * * @author Guillaume Anciaux * * @date Mon Oct 28 19:23:14 2013 * * @brief Set of macros to trace an atom through processor migrations * * @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 . * */ #ifndef __LIBMULTISCALE_TRACE_ATOM_H__ #define __LIBMULTISCALE_TRACE_ATOM_H__ /* -------------------------------------------------------------------------- */ -#ifdef TRACE_ATOM +#ifdef LIBMULTISCALE_TRACE_ATOM /* -------------------------------------------------------------------------- */ #include "lm_common.hh" __BEGIN_LIBMULTISCALE__ /* -------------------------------------------------------------------------- */ extern Real trace_ref[3]; extern Real trace_tolerance; extern UInt internal_index; -__END_LIBMULTISCALE__ - /* -------------------------------------------------------------------------- */ -#define ATOM_OUT(X) "(" << X[0] << "," << X[1] << "," << X[2] << ")" +template bool CHECK_TRACED(V &&X) { + return (trace_ref[0] - trace_tolerance <= X[0] && + X[0] <= trace_ref[0] + trace_tolerance && + trace_ref[1] - trace_tolerance <= X[1] && + X[1] <= trace_ref[1] + trace_tolerance && + trace_ref[2] - trace_tolerance <= X[2] && + X[2] <= trace_ref[2] + trace_tolerance); +} /* -------------------------------------------------------------------------- */ -#define CHECK_TRACED(X) \ - if (trace_ref[0] - trace_tolerance <= X[0] && \ - X[0] <= trace_ref[0] + trace_tolerance && \ - trace_ref[1] - trace_tolerance <= X[1] && \ - X[1] <= trace_ref[1] + trace_tolerance && \ - trace_ref[2] - trace_tolerance <= X[2] && \ - X[2] <= trace_ref[2] + trace_tolerance) +template bool CHECK_REF_TRACED(Ref &&ref) { + return CHECK_TRACED(ref.position0()); +} /* -------------------------------------------------------------------------- */ -#define CHECK_REF_TRACED(ref) \ - Real tmp_pos[3]; \ - ref.getPositions0(tmp_pos); \ - CHECK_TRACED(tmp_pos) + +__END_LIBMULTISCALE__ + +#define ATOM_OUT(X) "(" << X[0] << "," << X[1] << "," << X[2] << ")" /* -------------------------------------------------------------------------- */ + #define IF_COORD_EQUALS(x1, x2) \ if (x1[0] - trace_tolerance <= x2[0] && x2[0] <= x1[0] + trace_tolerance && \ x1[1] - trace_tolerance <= x2[1] && x2[1] <= x1[1] + trace_tolerance && \ x1[2] - trace_tolerance <= x2[2] && x2[2] <= x1[2] + trace_tolerance) /* -------------------------------------------------------------------------- */ #define IF_TRACED(X, message) \ - CHECK_TRACED(X) { DUMP(ATOM_OUT(X) << message, DBG_MESSAGE); } + if (CHECK_TRACED(X)) { \ + DUMP(ATOM_OUT(X) << message, DBG_MESSAGE); \ + } /* -------------------------------------------------------------------------- */ #define IF_REF_TRACED(ref, message) \ - CHECK_REF_TRACED(ref) { DUMP(ref << " " << message, DBG_MESSAGE); } + if (CHECK_REF_TRACED(ref)) { \ + DUMP(ref << " " << message, DBG_MESSAGE); \ + } /* -------------------------------------------------------------------------- */ #define VIEW_ATOM(T) \ { \ if (internal_index != UINT_MAX) { \ T ref(internal_index); \ DUMP("how is my atom now ? " << ref, DBG_MESSAGE); \ } else \ DUMP("my atom is not here ", DBG_MESSAGE); \ }; /* -------------------------------------------------------------------------- */ #define SET_INTERNAL_TRACE_INDEX(X, index) \ - CHECK_TRACED(X) \ - internal_index = index; + if (CHECK_TRACED(X)) \ + internal_index = index; /* -------------------------------------------------------------------------- */ #define SET_REF_INTERNAL_TRACE_INDEX(ref, index) \ { \ - CHECK_REF_TRACED(ref) \ - internal_index = index; \ + if (CHECK_REF_TRACED(ref)) \ + internal_index = index; \ } -/* -------------------------------------------------------------------------- */ +/* -------------------------------------------------------------------------- + */ #else #define IF_TRACED(X, message) #define IF_REF_TRACED(ref, message) #define VIEW_ATOM(T) #define CHECK_TRACED(X) if (1 == 2) #define CHECK_REF_TRACED(X) if (1 == 2) #define IF_COORD_EQUALS(x1, x2) #define SET_INTERNAL_TRACE_INDEX(X, index) #define SET_REF_INTERNAL_TRACE_INDEX(ref, index) #endif -/* -------------------------------------------------------------------------- */ +/* -------------------------------------------------------------------------- + */ #endif /* __LIBMULTISCALE_TRACE_ATOM_H__ */ diff --git a/src/stimulation/stimulation_langevin.cc b/src/stimulation/stimulation_langevin.cc index 1aca1e6..500c4a9 100644 --- a/src/stimulation/stimulation_langevin.cc +++ b/src/stimulation/stimulation_langevin.cc @@ -1,400 +1,391 @@ /** * @file stimulation_langevin.cc * * @author Guillaume Anciaux * * @date Wed Jul 09 21:59:47 2014 * * @brief This stimulator applies a Langevin thermostat * * @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 . * */ #include "stimulation_langevin.hh" #include "container_array.hh" #include "factory_multiscale.hh" #include "lib_continuum.hh" #include "lib_dd.hh" #include "lib_md.hh" #include "lm_common.hh" #include "math_tools.hh" #include "reference_manager.hh" #include /* -------------------------------------------------------------------------- */ __BEGIN_LIBMULTISCALE__ /* -------------------------------------------------------------------------- */ StimulationLangevin::StimulationLangevin(const std::string &name) : LMObject(name) { damp = 100; temperature = 100; timestep = 1.; seed = 0; compute_work_flag = false; // initialize the compute auto output_names = std::list{ "temperature", "friction_work", "restitution_work", "cumulated_frictional_work", "cumulated_restitution_work"}; this->createArrayOutputs(output_names); for (auto f : output_names) { this->getOutputAsArray(f).assign(1, 0); } // initialize the per atom computes output_names = std::list{ "old_friction_per_atom", "old_restitution_per_atom", "friction_work_per_atom", "restitution_work_per_atom" }; this->createArrayOutputs(output_names); - // friction_restitution_work("friction_restitution:" + name), - // old_friction_per_atom("old_friction_per_atom:" + name), - // old_restitution_per_atom("old_restitution_per_atom:" + name), - // friction_work_per_atom("friction_work_per_atom:" + name), - // restitution_work_per_atom("restitution_work_per_atom:" + name) - cumulated_friction_work = 0; cumulated_restitution_work = 0; - - // register for computes - // FilterManager::getManager().addObject(friction_restitution_work); } /* -------------------------------------------------------------------------- */ StimulationLangevin::~StimulationLangevin() {} /* -------------------------------------------------------------------------- */ template void StimulationLangevin::stimulate(Cont &cont) { if (compute_work_flag && compute_work_peratom_flag) this->addForce(cont); else if ((!compute_work_flag) && compute_work_peratom_flag) this->addForce(cont); else if ((!compute_work_flag) && (!compute_work_peratom_flag)) this->addForce(cont); else if (compute_work_flag && (!compute_work_peratom_flag)) this->addForce(cont); else LM_FATAL("internal error: should not happend"); } /* -------------------------------------------------------------------------- */ template void StimulationLangevin::resetOldWorkPerAtom(Cont &cont) { constexpr UInt Dim = Cont::Dim; if (this->compute_work_peratom_flag) { getOutputAsArray("friction_work_per_atom").resize(cont.size(), Dim); getOutputAsArray("restitution_work_per_atom").resize(cont.size(), Dim); } auto &old_friction_per_atom = getOutputAsArray("old_friction_per_atom"); auto &old_restitution_per_atom = getOutputAsArray("old_restitution_per_atom"); DUMP("AAAAA " << cont.size() << " " << old_friction_per_atom.size() << " " << old_friction_per_atom.rows() << " " << old_friction_per_atom.cols(), DBG_DETAIL); if (old_friction_per_atom.size() == 0) { old_friction_per_atom.resize(cont.size(), Dim); old_restitution_per_atom.resize(cont.size(), Dim); cont.attachObject(old_friction_per_atom); cont.attachObject(old_restitution_per_atom); } else { if (old_friction_per_atom.size() != cont.size() * Dim) LM_FATAL("inconsistent previous storage of force vector (friction) " << old_friction_per_atom.size() << "(" << &old_friction_per_atom << ")" << " != " << cont.size() * Dim << "(" << &cont << ")"); if (old_restitution_per_atom.size() != cont.size() * Dim) LM_FATAL("inconsistent previous storage of force vector (restitution) " << old_restitution_per_atom.size() << "(" << &old_restitution_per_atom << ")" << " != " << cont.size() * Dim << "(" << &cont << ")"); } } /* -------------------------------------------------------------------------- */ template void StimulationLangevin::gatherWork(Cont &cont) { // Communicator &comm = Communicator::getCommunicator(); CommGroup &group = cont.getCommGroup(); if (compute_work_flag) { group.allReduce(&friction_work, 1, "frictional work", OP_SUM); group.allReduce(&restitution_work, 1, "restitution work", OP_SUM); friction_work *= 0.5 * timestep * code_unit_system.fd2e; restitution_work *= 0.5 * timestep * code_unit_system.fd2e; cumulated_friction_work += friction_work; cumulated_restitution_work += restitution_work; std::map toto{ {"temperature", temperature}, {"friction_work", friction_work}, {"restitution_work", restitution_work}, {"cumulated_frictional_work", cumulated_friction_work}, {"cumulated_restitution_work", cumulated_restitution_work}}; for (auto &&kv : toto) { auto &name = kv.first; auto &value = kv.second; auto &f = getOutputAsArray(name); f.clear(); f.push_back(value); } } } /* -------------------------------------------------------------------------- */ template void StimulationLangevin::addForce(Cont &cont) { if (computeWork) resetOldWorkPerAtom(cont); this->frictionCorrection(cont); this->restitutionCorrection(cont); if (computeWork) gatherWork(cont); } /* -------------------------------------------------------------------------- */ template void __attribute__((noinline)) StimulationLangevin::frictionCorrection(Cont &cont) { // constexpr UInt Dim = Cont::Dim; Real conv1 = code_unit_system.mv_t2f / damp; UInt index = 0; friction_work = 0; auto &old_friction_per_atom = getOutputAsArray("old_friction_per_atom"); auto &friction_work_per_atom = getOutputAsArray("friction_work_per_atom"); for (auto &&at : cont) { Real gamma1 = -conv1 * at.mass(); auto &&vel = at.velocity(); Vector friction_term = gamma1 * vel; if (computeWork) { Vector old_friction_term = old_friction_per_atom.row(index); friction_work += (friction_term + old_friction_term).dot(vel); old_friction_per_atom.row(index) = friction_term; if (computeWorkPerAtom) friction_work_per_atom.row(index) += friction_work; // accountFrictionWork( // friction_term, at.velocity()[i], index); } at.force() += friction_term; ++index; } } /* -------------------------------------------------------------------------- */ template void __attribute__((noinline)) StimulationLangevin::restitutionCorrection(Cont &cont) { // constexpr UInt Dim = Cont::Dim; Real conv2 = sqrt(24. * temperature * boltzmann * code_unit_system.kT2fd * code_unit_system.m_tt2f_d / damp / timestep); restitution_work = 0; auto &old_restitution_per_atom = getOutputAsArray("old_restitution_per_atom"); auto &restitution_work_per_atom = getOutputAsArray("restitution_work_per_atom"); UInt index = 0; for (auto &&at : cont) { Real gamma2 = conv2 * sqrt(at.mass()); Vector restitution_term = gamma2 * 0.5 * Vector::Random(); if (computeWork) { auto &&vel = at.velocity(); Vector old_restitution_term = old_restitution_per_atom.row(index); restitution_work += (restitution_term + old_restitution_term).dot(vel); old_restitution_term = restitution_term; if (computeWorkPerAtom) restitution_work_per_atom[index] += restitution_work; // accountRestitutionWork( // restitution_term, at.velocity()[i], index); } at.force() += restitution_term; ++index; } } /* -------------------------------------------------------------------------- */ // template // inline void StimulationLangevin::accountFrictionWork(Real friction_term, // Real velocity, // UInt index) { // Real old_friction_term = old_friction_per_atom[index]; // friction_work += (friction_term + old_friction_term) * velocity; // old_friction_per_atom[index] = friction_term; // if (computeWorkPerAtom) // friction_work_per_atom[index] += friction_work; // } // /* -------------------------------------------------------------------------- // */ // template // inline void StimulationLangevin::accountRestitutionWork(Real // restitution_term, // Real velocity, // UInt index) { // Real old_restitution_term = old_restitution_per_atom[index]; // restitution_work += (restitution_term + old_restitution_term) * velocity; // old_restitution_per_atom[index] = restitution_term; // if (computeWorkPerAtom) // restitution_work_per_atom[index] += restitution_work; // } /* -------------------------------------------------------------------------- */ /* LMDESC LANGEVIN This stimulator applies a Langevin thermostat to domain. The equation of motion becomes:\\ \begin{equation} m\ddot{x} = - \nabla \phi \underbrace{- \frac{m}{damp}{v}}_{f_{damp}} + \underbrace{\sqrt{\frac{k_{b}\,T\,m}{dt\,damp}}R}_{f_{rand}} \end{equation} with $\phi$ the classical potential for interatomic forces, $T$ the temperature of the heat bath, $damp$ the relaxation time, $v$ the velocity field, $k_b$ the boltzmann constant, $dt$ the used timestep, $m$ the particle mass and $R$ a vector of random numbers uniformly distributed between -1/2 and 1/2. \\ Explicitely, this stimulator append $f_{damp}$ and $f_{rand}$ to the force field during stage PRE\_STEP3. */ /* LMEXAMPLE STIMULATION thermostat LANGEVIN INPUT md TEMP 100 SEED 32 */ /* LMHERITANCE action_interface */ void StimulationLangevin::declareParams() { StimulationInterface::declareParams(); this->changeDefault("STAGE", PRE_STEP3); /* LMKEYWORD DAMP Set the damping parameter expressed in the units of time */ this->parseKeyword("DAMP", damp); /* LMKEYWORD TEMP Set the temperature desired to be set */ this->parseKeyword("TEMP", temperature); /* LMKEYWORD SEED Set the seed for the random generator */ this->parseKeyword("SEED", seed); /* LMKEYWORD TIMESTEP Set the time step needed for random noise contruction */ this->parseKeyword("TIMESTEP", timestep); /* LMKEYWORD WORK_PERATOM TODO */ this->parseTag("WORK_PERATOM", compute_work_peratom_flag, false); /* LMKEYWORD WORK Activates the computation of the work done by ther thermostat. A compute named friction\_restitution:ID is created and allows the visualization of the work done by the thermostat (ID should be the name given to this compute). More details available The work is computed like:\ \ $$dW^n = \int_{x^n}^{x^{n+1}} f \cdot dx$$ which for the case of the velocity verlet scheme is equivalent to $$dW^n = \frac{1}{2}(f^n + f^{n+1}) v^{n+1/2}$$ Thus when this flag is activated a compute with with the name "friction\_restitution:STIMULATOR\_NAME" and 5 entries is registered. These entries are: \begin{enumerate} \item Requested temperature \item frictional\_work: $dW_{damp}^n = \frac{1}{2}(f_{damp}^n + f_{damp}^{n+1}) v^{n+1/2}$ \item restitution\_work: $dW_{restitution}^n = \frac{1}{2}(f_{rand}^n + f_{rand}^{n+1}) v^{n+1/2}$ \item cumulated\_frictional\_work: $W_{damp}^n = \sum_0^n dW_{damp}^n$ \item cumulated\_restitution\_work: $W_{restitution}^n = \sum_0^n dW_{restitution}^n$ \end{enumerate} */ this->parseTag("WORK", compute_work_flag, false); } /* -------------------------------------------------------------------------- */ void StimulationLangevin::init() { MathTools::setSeed(seed); } /* -------------------------------------------------------------------------- */ DECLARE_STIMULATION_MAKE_CALL(StimulationLangevin) __END_LIBMULTISCALE__