diff --git a/src/common/attached_object.hh b/src/common/attached_object.hh index 73bbd82..817dc20 100644 --- a/src/common/attached_object.hh +++ b/src/common/attached_object.hh @@ -1,120 +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.h" /* -------------------------------------------------------------------------- */ __BEGIN_LIBMULTISCALE__ /* -------------------------------------------------------------------------- */ template class ContainerArray; class AttachedObject : public AutoDispatch::ArgumentAny { public: virtual ~AttachedObject(){}; virtual void packData(UInt index, LAMMPS_NS::PackBuffer &buffer, bool verbose = false) = 0; virtual void unpackData(UInt index, LAMMPS_NS::PackBuffer &buffer, 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, LAMMPS_NS::PackBuffer &buffer, bool verbose = false) override { auto &v = this->getV(); for (UInt i = 0; i < v.cols(); ++i) { - if (verbose) { - DUMP("packing: " << v(index, i), DBG_MESSAGE); - } buffer << v(index, i); } + if (verbose) { + DUMP("packing: " << v.row(index), DBG_MESSAGE); + } } void unpackData(UInt index, LAMMPS_NS::PackBuffer &buffer, 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: " << 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/ref_subset.cc b/src/common/ref_subset.cc index f7dc19f..2c7421c 100644 --- a/src/common/ref_subset.cc +++ b/src/common/ref_subset.cc @@ -1,370 +1,458 @@ /** * @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 /* -------------------------------------------------------------------------- */ #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(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( // LAMMPS_NS::PackBuffer &buffer) { // LM_TOIMPLEMENT; // // for (auto &&[proc, listsent] : sent_indexes) { // // if (proc == UINT_MAX) // // continue; // // 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) { // // DUMP("prepares to pack obj " << att_obj->getID(), DBG_MESSAGE); // // att_obj->packData(sent_index, buffer[proc], proc, true); // // } // // } // // } // } /* ------------------------------------------------------------------------ */ // template // void RefSubset::unpackAttachedData( // LAMMPS_NS::PackBuffer &buffer) { // LM_TOIMPLEMENT; // // for (auto &&[proc, listrecv] : recv_indexes) { // // if (proc == UINT_MAX) // // continue; // // 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) { // // 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 +std::enable_if_t>, void> +printInfoAttachedValue(const std::string &mess, T1 &att_obj, T2 &container, + UInt i_src) { + DUMP(container.getID() << "[" << i_src << "] = " << container.get(i_src) + << " attached:" << att_obj.getID() << ": " + << att_obj.getV().row(i_src) << " ### " << mess, + DBG_MESSAGE); +} + +template +std::enable_if>, void> +printInfoAttachedValue(const std::string &, T1 &, T2 &, UInt, UInt) {} + +/* ----------------------------------------------------------------------- */ + template void RefSubset::packData(Ref &ref, LAMMPS_NS::PackBuffer &buffer, bool verbose) { auto index = this->getRefIndex(ref); + if (index == -1) + LM_FATAL("this should not happen"); for (auto &&[ID, att_obj] : this->attached_objects) { - if (index != -1) - att_obj->packData(index, buffer, verbose); + att_obj->packData(index, buffer, verbose); + if (true) { // CHECK_REF_TRACED(ref)) { + if (auto *ptr = dynamic_cast *>(att_obj.get())) { + printInfoAttachedValue("packing", *ptr, container, index); + } + } } } /* ----------------------------------------------------------------------- */ template void RefSubset::unpackData(Ref &ref, LAMMPS_NS::PackBuffer &buffer, bool verbose) { - for (auto &&[ID, att_obj] : this->attached_objects) { - auto index = this->getRefIndex(ref); + + auto index = this->getRefIndex(ref); + + if constexpr (is_subset) { if (index != -1) - att_obj->unpackData(index, buffer, verbose); + LM_FATAL("this should not happen: found " << container.getID() << "[" + << index << "] = " << ref); + + index = holes.back(); + this->container.get(index) = ref; + holes.pop_back(); + } + + for (auto &&[ID, att_obj] : this->attached_objects) { + att_obj->unpackData(index, buffer, verbose); + if (true) { // CHECK_REF_TRACED(container.get(index))) { + if (auto *ptr = dynamic_cast *>(att_obj.get())) { + printInfoAttachedValue("unpacking", *ptr, container, index); + } + } } } /* ----------------------------------------------------------------------- */ template void RefSubset::translateMovingReference(Ref &src, Ref &dst) { if constexpr (not is_subset) return; - // loop over the contained refs and shift references of moved refs + // get the index of the reference to be shifted in the subset auto i_src = this->getRefIndex(src); auto i_dst = this->getRefIndex(dst); // check if the moved(src) ref atom is in this subset // if not I am not concerned by this change if (i_src == -1) { return; } - container.get(i_src) = dst; + + if (CHECK_REF_TRACED(container.get(i_src))) { + for (auto &&[ptr, att_obj] : this->attached_objects) { + if (auto *ptr = dynamic_cast *>(att_obj.get())) { + printInfoAttachedValue("before translate", *ptr, container, i_src); + } + } + } // check if the ref atom of the destination is in this subset // This ref atom points certainly to an atom that left the processor - // then I have a hole to register + // if it is in this subset then I have a hole to register + container.get(i_src) = dst; + if (i_dst != -1) { holes.push_back(i_dst); } + + if (CHECK_REF_TRACED(container.get(i_src))) { + for (auto &&[ptr, att_obj] : this->attached_objects) { + if (auto *ptr = dynamic_cast *>(att_obj.get())) { + printInfoAttachedValue("after translate", *ptr, container, i_src); + } + } + } } /* -------------------------------------------------------------------------- */ 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 moveLastReferenceToHole(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); + IF_REF_TRACED(container.get(i_src), + " = " << container.getID() << "[" << i_src << "]: " + << " traced ref moved from index " << i_src + << " to index " << i_dest << " for container " + << container.getID()); - DUMP("ref moved from index " << i_src << " to index " << i_dest - << " for container " << &container, - DBG_INFO); + if (CHECK_REF_TRACED(container.get(i_src))) { + for (auto &&[ptr, att_obj] : this->attached_objects) { + if (auto *ptr = dynamic_cast *>(att_obj.get())) { + printInfoAttachedValue("before moving ", *ptr, container, i_src); + } + } + } 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); + if (CHECK_REF_TRACED(container.get(i_dest))) { + if (auto *ptr = dynamic_cast *>(att_obj.get())) { + printInfoAttachedValue("after moving attached ", *ptr, container, + i_dest); + } + } } } /* ------------------------------------------------------------------------ */ template void RefSubset::moveLastReferenceToHole(UInt i_dest) { if constexpr (not is_subset) return; UInt nRefs = container.size(); if (nRefs - 1 == i_dest) return; + + if (i_dest > nRefs - 1) { + int a = 1; + while (a) { + } + } + 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; } /* ----------------------------------------------------------------------- */ template int RefSubset::getRefIndex(Ref &ref) { for (auto &&[i, e] : enumerate(container)) { if (e == ref) return i; } return -1; } /* ----------------------------------------------------------------------- */ +template +void RefSubset::fillMasks(MapRefToUInt &masks) { + + for (auto &&e : container) { + if (masks.count(e) != 0) + masks[e] |= this->mask(); + else + masks[e] = this->mask(); + } +} + +/* ----------------------------------------------------------------------- */ + #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 3f8db5b..c4bfd5e 100644 --- a/src/common/ref_subset.hh +++ b/src/common/ref_subset.hh @@ -1,202 +1,205 @@ /** * @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 translateMovingReference(Ref &src, Ref &dst) = 0; //! pack attached data // virtual void packAttachedData(LAMMPS_NS::PackBuffer &buffer) = 0; // //! unpack attached data // virtual void unpackAttachedData(LAMMPS_NS::PackBuffer &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; //! packing all data necessary for the subsets and attached objects virtual void packData(Ref &ref, LAMMPS_NS::PackBuffer &buffer, bool verbose = false) = 0; //! packing all data necessary for the subsets and attached objects virtual void unpackData(Ref &ref, LAMMPS_NS::PackBuffer &buffer, bool verbose = false) = 0; //! attach a vector of data to the subset void attachData(std::shared_ptr &obj); //! detach a vector of data from the subset void detachData(std::shared_ptr &obj); virtual bool isSubset() = 0; virtual int mask() = 0; + virtual void fillMasks(MapRefToUInt &masks) = 0; // protected: //! list of 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 translateMovingReference(Ref &src, Ref &dst); //! fill remaining holes so that container is one pack void fillRemainingHoles() override; //! build the mask for the sent refs void buildMasksForSentRefs(MapRefToUInt &masks); //! packing all data necessary for the subsets and attached objects void packData(Ref &ref, LAMMPS_NS::PackBuffer &buffer, bool verbose = false) override; //! packing all data necessary for the subsets and attached objects void unpackData(Ref &ref, LAMMPS_NS::PackBuffer &buffer, bool verbose = false) override; int getRefIndex(Ref &ref); // //! pack attached data // void packAttachedData(LAMMPS_NS::PackBuffer &buffer) override; // //! unpack attached data // void unpackAttachedData(LAMMPS_NS::PackBuffer &buffer) override; //! move last reference in the array (useful to fill holes) void moveLastReferenceToHole(UInt i_dest); //! move reference in the subset 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; }; int mask() override { return 1 << this->subset_id; }; + void fillMasks(MapRefToUInt &masks); + 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 void RefSubsetInterface::attachData(std::shared_ptr &obj) { DUMP("want to register obj " << &obj, DBG_DETAIL); 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->getID()] = obj; } /* -------------------------------------------------------------------------- */ // detach a vector of data from the subset template void RefSubsetInterface::detachData(std::shared_ptr &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 a302514..dc4c08d 100644 --- a/src/common/reference_manager.cc +++ b/src/common/reference_manager.cc @@ -1,452 +1,464 @@ /** * @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::computeMasks() { + + masks.clear(); + for (auto &&[name, s] : this->subsets) { + s->fillMasks(masks); + } +} +/* -------------------------------------------------------------------------- */ + // 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() { LM_TOIMPLEMENT; // MapRefToUInt masks_send; // MapRefToUInt masks_recv; // comm_group->synchronize(); // 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); // 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(); // //#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(); // //#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"); - } - // 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"); - - 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"); - - 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::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"); +// } +// // 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"); + +// 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"); + +// 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; +// template void ReferenceManager::printBilan() { +// std::map nb_sent_by_proc; - // if (moved.size()) - // DUMP(moved.size() << " refs moved in memory", DBG_WARNING); +// 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, 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 &&[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 &&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"); - // } -} +// 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; } #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/common/reference_manager.hh b/src/common/reference_manager.hh index c7f7d1c..202a9d4 100644 --- a/src/common/reference_manager.hh +++ b/src/common/reference_manager.hh @@ -1,250 +1,254 @@ /** * @file reference_manager.hh * * @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. * */ #ifndef __LIBMULTISCALE_REFERENCE_MANAGER_HH__ #define __LIBMULTISCALE_REFERENCE_MANAGER_HH__ /* -------------------------------------------------------------------------- */ #include "auto_arguments.hh" #include "possible_types.hh" #include "ref_subset.hh" #include "reference_manager_interface.hh" /* -------------------------------------------------------------------------- */ #include #include /* -------------------------------------------------------------------------- */ __BEGIN_LIBMULTISCALE__ template class ReferenceManager : public ReferenceManagerInterface { public: using BufferMap = std::map; using MapRefToUInt = std::map; using MapRefToRef = std::map; using MapUIntToRefList = std::map>; ReferenceManager(typename Ref::Domain::ContainerPoints &global_container); virtual ~ReferenceManager() = default; //! function that start the updating process of all attached structures void updateRefSubSets(); //! request to manage a subset DECORATE_FUNCTION_DISPATCH(addSubset, _or) template std::enable_if_t> addSubset(Cont &cont); template std::enable_if_t> addSubset(Cont &) { LM_FATAL("wrong reference type";) } void addSubset(ContainerInterface &cont) override; //! request to remove a subset void removeSubset(const LMID &name) override; void attachObject(std::shared_ptr obj, ContainerInterface &c) override; //! request detaching a generic AttachedObject with a given container c void detachObject(std::shared_ptr obj, ContainerInterface &c) override; //! print bilan for concerned container : used only to debug - void printBilan(); + // void printBilan(); //! retreive the mpi context void acquireContext(const LMObject &obj) override; + //! compute the masks of every reference attached with subsets + void computeMasks(); protected: //! set the mpi communicator void setCommGroup(CommGroup &comm); //! generic communication routine to send/recv a bunch of buffers - void exchangeBuffers(BufferMap &toSend, BufferMap &toRecv); + // void exchangeBuffers(BufferMap &toSend, BufferMap &toRecv); + //! translate the references for moved refs // void translateMovingReferences(); //! pack masks to a BufferMap // void packMasks(MapRefToUInt &masks); //! unpack masks to a BufferMap // void unpackMasks(MapRefToUInt &masks); //! clear the buffers and/or create one entry per proc I should com with // void clearPackBuffers(); protected: //! global set // RefSet global_set; //! subset array std::map>> subsets; //! mapping between sent ref ref and new proc owner // MapRefToUInt sent; //! inverse mapping : sent refs sorted by receiving processors // MapUIntToRefList sent_byproc; //! new refs generated by migration (received) // MapUIntToRefList newrefs; //! communication buffers to be sent // BufferMap buffers_tosend; //! communication buffers to be received // BufferMap buffers_torecv; //! mapping between moved refs old ref and new ref // MapRefToRef moved; //! flag to notify if updating of attached references needed // bool have_changed = false; //! for debug purpose // void printPackBuffersStatus(); + typename ReferenceManager::MapRefToUInt masks; private: CommGroup *comm_group; }; /* -------------------------------------------------------------------------- */ template class RefPointData; template struct ReferenceManager> : public ReferenceManagerInterface { void removeSubSet(const LMID &){}; }; /* -------------------------------------------------------------------------- */ template class RefGenericElem; template struct ReferenceManager> : public ReferenceManagerInterface { void removeSubSet(const LMID &){}; }; /* -------------------------------------------------------------------------- */ template class RefGenericDDElem; template struct ReferenceManager> : public ReferenceManagerInterface { void removeSubSet(const LMID &){}; }; /* -------------------------------------------------------------------------- */ template class RefGenericDDNode; template struct ReferenceManager> : public ReferenceManagerInterface { void removeSubSet(const LMID &){}; }; /* -------------------------------------------------------------------------- */ template <> struct ReferenceManager : public ReferenceManagerInterface { void removeSubSet(const LMID &) {} }; /* -------------------------------------------------------------------------- */ template <> struct ReferenceManager : public ReferenceManagerInterface { void removeSubSet(const LMID &) {} }; /* -------------------------------------------------------------------------- */ class RefElemMeca1D; template <> struct ReferenceManager : public ReferenceManagerInterface { void removeSubSet(const LMID &) {} }; /* -------------------------------------------------------------------------- */ template class RefElemAkantu; template struct ReferenceManager> : public ReferenceManagerInterface { void removeSubSet(const LMID &) {} }; /* -------------------------------------------------------------------------- */ class RefElemParaDiS; template <> struct ReferenceManager : public ReferenceManagerInterface { void removeSubSet(const LMID &) {} }; /* -------------------------------------------------------------------------- */ template void ReferenceManager::attachObject(std::shared_ptr obj, ContainerInterface &cont) { bool found = false; for (auto &&[name, subset] : subsets) { if (&subset->getContainer() == &cont) { DUMP("attaching object " << &obj << " to subset " << cont.getID(), DBG_INFO); subset->attachData(obj); found = true; break; } } if (not found) LM_FATAL("attaching object failed : did not find its subset"); } /* ------------------------------------------------------------------------ */ template void ReferenceManager::detachObject(std::shared_ptr obj, ContainerInterface &cont) { bool found = false; for (auto &&[name, subset] : subsets) { if (&subset->getContainer() == &cont) { DUMP("detaching object " << &obj << " from subset " << cont.getID(), DBG_INFO); subset->detachData(obj); found = true; break; } } if (not found) LM_FATAL("detaching object failed : did not find its subset"); } /* -------------------------------------------------------------------------- */ __END_LIBMULTISCALE__ #endif /* __LIBMULTISCALE_REFERENCE_MANAGER_HH__ */ diff --git a/src/md/lammps/domain_lammps.cc b/src/md/lammps/domain_lammps.cc index 9c46373..83f377b 100644 --- a/src/md/lammps/domain_lammps.cc +++ b/src/md/lammps/domain_lammps.cc @@ -1,1080 +1,1081 @@ /** * @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 "atom_vec_lm.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 = nullptr; /* -------------------------------------------------------------------------- */ __BEGIN_LIBMULTISCALE__ template void DomainLammps::init() { DomainAtomic, Dim>::init(); if (!this->getCommGroup().amIinGroup()) return; this->initModel(); this->update->whichflag = 1; this->update->dt = this->timeStep; LAMMPS_NS::LAMMPS::init(); this->update->integrate->setup(0); this->initDOFs(); this->ref_manager->setState(true); this->atoms.getBoundingBox() = this->getDomainBoundingBox(); this->updateAcceleration(); } /* -------------------------------------------------------------------------- */ 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: if (not is_mass_per_particle) return ArrayView(&lammps_main_object->atom->mass[0], lammps_main_object->atom->ntypes, 1); return ArrayView(&lammps_main_object->atom->rmass[0], lammps_main_object->atom->natoms, 1); case _acceleration: return ArrayView(this->acceleration_field.data(), lammps_main_object->atom->nlocal, 3); case _angular_velocity: if (not lammps_main_object->atom->omega) { LM_FATAL("created atoms do not have angle degrees of freedom"); } return ArrayView(&lammps_main_object->atom->omega[0][0], lammps_main_object->atom->nlocal, 3); case _angular_acceleration: return ArrayView(this->angular_acceleration_field.data(), lammps_main_object->atom->nlocal, 3); case _torque: if (not lammps_main_object->atom->torque) { LM_FATAL("created atoms do not have angle degrees of freedom"); } return ArrayView(&lammps_main_object->atom->torque[0][0], lammps_main_object->atom->nlocal, 3); default: LM_FATAL("Not accessible field: " << field_type); } } /* -------------------------------------------------------------------------- */ template std::string DomainLammps::createLammpsHeader() { std::string unit_name; if (code_unit_system == metal_unit_system) unit_name = "metal"; else if (code_unit_system == real_unit_system) unit_name = "real"; else if (code_unit_system == si_unit_system) unit_name = "si"; else LM_FATAL("unknown units"); std::string file_to_generate = "./generated_lammps.config"; std::ofstream file_lammps(file_to_generate.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); return file_to_generate; } /* -------------------------------------------------------------------------- */ 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); } // ************************************************ // read lammps config file // ************************************************ std::string file_to_read = lammpsfile; DUMP("lammpsfile" << " = " << lammpsfile, DBG_INFO); // if necessary creates on the fly the lammps confif file if (create_header_flag && this->getCommGroup().getMyRank() == 0) { file_to_read = this->createLammpsHeader(); } 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 !!!"); // ************************************************ // execute 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 fixes 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(); if (auto *ptr = dynamic_cast(this->atom->avec)) { auto *avec = new AtomVecLM(*ptr); avec->setRefManager(ref_manager); this->atom->avec = avec; } 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->mvv2e = code_unit_system.mvv2e; } if (this->domain->dimension != Dim) LM_FATAL("mismatch dimension: check lammps config file"); this->acceleration_field.resize(atom->nlocal, 3); this->acceleration_field.setZero(); if (std::string(this->atom->atom_style) != "atomic") { this->is_mass_per_particle = true; this->angular_acceleration_field.resize(atom->nlocal, 3); this->angular_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(real_unit_system); else if (strcmp(update->unit_style, "metal") == 0) uc.setInUnits(metal_unit_system); else if (strcmp(update->unit_style, "si") == 0) uc.setInUnits(si_unit_system); else LM_FATAL("not known lammps unit system:" << " things need to be done to integrate"); uc.computeConversions(); if (not is_mass_per_particle) { for (int i = 1; i <= atom->ntypes; i++) atom->mass[i] *= uc.getConversion(); } else { for (int i = 0; i <= atom->natoms; i++) atom->rmass[i] *= uc.getConversion(); } } #ifdef LIBMULTISCALE_TRACE_ATOM for (int i = 0; i < atom->nlocal; ++i) { 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; if (!this->getCommGroup().amIinGroup()) return 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(false); } /* -------------------------------------------------------------------------- */ template void DomainLammps::enforceCompatibility() { if (this->getCommGroup().amIinGroup()) { this->exchangeAtoms(); this->acceleration_field.resize(atom->nlocal, 3); } } /* -------------------------------------------------------------------------- */ template void DomainLammps::exchangeAtoms() { VIEW_ATOM(RefLammps); + this->ref_manager->computeMasks(); UInt nflag = neighbor->decide(); if (nflag == 0 && this->restart_start == false) { timer->stamp(); comm->forward_comm(); timer->stamp(LAMMPS_NS::Timer::COMM); } else { if (modify->n_pre_exchange) { timer->stamp(); modify->pre_exchange(); timer->stamp(LAMMPS_NS::Timer::MODIFY); } 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(LAMMPS_NS::Timer::COMM); if (modify->n_pre_neighbor) { modify->pre_neighbor(); timer->stamp(LAMMPS_NS::Timer::MODIFY); } neighbor->build(1); timer->stamp(LAMMPS_NS::Timer::NEIGH); if (modify->n_post_neighbor) { modify->post_neighbor(); timer->stamp(LAMMPS_NS::Timer::MODIFY); } 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(LAMMPS_NS::Timer::COMM); if (modify->n_pre_neighbor) modify->pre_neighbor(); neighbor->build(false); timer->stamp(LAMMPS_NS::Timer::NEIGH); this->restart_start = false; } this->ref_manager->compactSubsets(); 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 (lammps_main_object->atom->torque) { 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->exchangeAtoms(); 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(LAMMPS_NS::Timer::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; double *radius = atom->radius; double inertia = 0.4; 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); } if (std::string(this->atom->atom_style) != "atomic") { auto torque = this->getField(_torque); for (UInt i = 0; i < UInt(lammps_main_object->atom->nlocal); ++i) { Real mass_inverse = code_unit_system.ft_m2v / inertia / (radius[i] * radius[i] * rmass[i]); this->angular_acceleration_field(i, 0) = mass_inverse * torque(i, 0); this->angular_acceleration_field(i, 1) = mass_inverse * torque(i, 1); this->angular_acceleration_field(i, 2) = mass_inverse * torque(i, 2); } } } /* -------------------------------------------------------------------------- */ template void DomainLammps::updateGradient() { if (!this->getCommGroup().amIinGroup()) return; // informing lammps of the current step this->update->ntimestep = current_step; 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); timer->stamp(LAMMPS_NS::Timer::MODIFY); } if (this->force->pair && this->force->pair->compute_flag) { this->force->pair->compute(eflag, vflag); this->timer->stamp(LAMMPS_NS::Timer::PAIR); } 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(LAMMPS_NS::Timer::BOND); } if (this->force->kspace) { this->force->kspace->compute(eflag, vflag); this->timer->stamp(LAMMPS_NS::Timer::KSPACE); } if (modify->n_pre_reverse) { modify->pre_reverse(eflag, vflag); timer->stamp(LAMMPS_NS::Timer::MODIFY); } if (this->force->newton) { this->comm->reverse_comm(); this->timer->stamp(LAMMPS_NS::Timer::COMM); } if (this->modify->n_post_force) this->modify->post_force(vflag); STOPTIMER("AtomStepForces"); VIEW_ATOM(RefLammps); this->changeRelease(); } /* -------------------------------------------------------------------------- */ /* LMDESC LAMMPS This plugin make the interface with lammps */ /* LMHERITANCE domain_md */ /* LM EXAMPLE Section MultiScale RealUnits\\ ...\\ GEOMETRY myGeom CUBE BBOX -1 1 -1 1 -1 1\\ MODEL LAMMPS md\\ ...\\ endSection\\ \\ Section LAMMPS:md RealUnits \\ 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 void RefLammps::test() { // if (this->index_atom == -1) // DUMP("toto", DBG_MESSAGE); } /* -------------------------------------------------------------------------- */ template class DomainLammps<2>; template class DomainLammps<3>; __END_LIBMULTISCALE__ diff --git a/src/md/lammps/ref_atom_lammps.hh b/src/md/lammps/ref_atom_lammps.hh index f46c346..7048ad1 100644 --- a/src/md/lammps/ref_atom_lammps.hh +++ b/src/md/lammps/ref_atom_lammps.hh @@ -1,318 +1,320 @@ /** * @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 "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, _angular_velocity, _torque, _radius>; RefLammps(UInt id) : index_atom(id), _weight(1), stress_hook(NULL), epot_hook(NULL) { this->test(); }; 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 VectorView<_Dim> angular_velocity(); inline VectorView<_Dim> angular_acceleration(); inline VectorView<_Dim> torque(); inline Real &radius(); 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; void test() override; private: int index_atom; Real _weight; StressHookLammps *stress_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() { this->test(); return VectorView(lammps_main_object->atom->x0[index_atom]); } /* -------------------------------------------------------------------------- */ template inline VectorView RefLammps::angular_velocity() { return VectorView(lammps_main_object->atom->omega[index_atom]); } /* -------------------------------------------------------------------------- */ template inline VectorView RefLammps::torque() { return VectorView(lammps_main_object->atom->torque[index_atom]); } /* -------------------------------------------------------------------------- */ template Real &RefLammps::radius() { return lammps_main_object->atom->radius[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() { LM_TOIMPLEMENT; // 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() { if (lammps_main_object->atom->mass) return lammps_main_object->atom ->mass[lammps_main_object->atom->type[index_atom]]; return lammps_main_object->atom->rmass[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(); Eigen::IOFormat fmt(Eigen::FullPrecision, 0, ", ", "", "", "", "[", "]"); - os << "Lammps reference : " << n << " P=" << ref.position().format(fmt) + os << "Ref[" << n << "] = " //<< " P=" << ref.position().format(fmt) << " P0=" << ref.position0().format(fmt) - << " V=" << ref.velocity().format(fmt) << " F=" << ref.force().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 [[gnu::unused]]) { 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 6d18239..8da3497 100644 --- a/src/md/lammps/reference_manager_lammps.cc +++ b/src/md/lammps/reference_manager_lammps.cc @@ -1,217 +1,188 @@ /** * @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(" << 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); - if (CHECK_REF_TRACED(refnew)) { - int a = 1; - while (a) { - } - } + // TRACE_BREAKPOINT(refnew) DUMP("I move atom " << refold << " to " << refnew, DBG_INFO); LM_ASSERT( !isNaN(refnew.position0()), "the atom of reference " << refnew << " is out of control : initial position is not-a-number (NaN)"); for (auto &&[name, s] : this->subsets) { s->translateMovingReference(refold, refnew); } - // 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 void ReferenceManagerLammps::compactSubsets() { for (auto &&[name, s] : this->subsets) { s->fillRemainingHoles(); } } /* -------------------------------------------------------------------------- */ template class ReferenceManagerLammps<2>; template class ReferenceManagerLammps<3>; __END_LIBMULTISCALE__ diff --git a/src/md/lammps/reference_manager_lammps.hh b/src/md/lammps/reference_manager_lammps.hh index e06086e..a74283c 100644 --- a/src/md/lammps/reference_manager_lammps.hh +++ b/src/md/lammps/reference_manager_lammps.hh @@ -1,113 +1,120 @@ /** * @file reference_manager_lammps.hh * * @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. * */ #ifndef __LIBMULTISCALE_REFERENCE_MANAGER_LAMMPS_HH__ #define __LIBMULTISCALE_REFERENCE_MANAGER_LAMMPS_HH__ /* -------------------------------------------------------------------------- */ #include "ref_atom_lammps.hh" #include "reference_manager.hh" +#include "trace_atom.hh" /* -------------------------------------------------------------------------- */ __BEGIN_LIBMULTISCALE__ /* -------------------------------------------------------------------------- */ using LMPPackBuffer = LAMMPS_NS::PackBuffer; template class ReferenceManagerLammps : public ReferenceManager> { public: ReferenceManagerLammps(ContainerLammps &global_container); // bool HaveChanged(); RefLammps &NewRef(RefLammps &ref); //! declare a migration to proc toProc // void sendAtom(UInt i, UInt toProc); //! declare a displacement in memory void moveAtom(UInt i, UInt j); //! declare a received atom from proc frmoProc // void acceptAtom(UInt i, UInt fromProc); void rebuildRefs(); void rebuildGenericRefs(); void compactSubsets(); // **************** new interface ********* void pack_subset_data(UInt i, LMPPackBuffer buf) { // build a reflammps RefLammps ref(i); + + IF_REF_TRACED(ref, "sending atom(" << ref.getIndex() << ")"); // first of all pack the masks - int mask = masks[ref]; + int mask = this->masks[ref]; buf << mask; // loop over each subset for (auto &&[name, s] : this->subsets) { if (mask & s->mask()) - s->packData(ref, buf); + s->packData(ref, buf, true); } + SET_REF_INTERNAL_TRACE_INDEX(ref, -1); } void unpack_subset_data(UInt i, LMPPackBuffer buf) { // build a reflammps RefLammps ref(i); + + DUMP("received " << ref, DBG_MESSAGE); + + IF_REF_TRACED(ref, "receiving atom(" << ref.getIndex() << ")"); // first of all pack the masks int mask; mask << buf; // loop over each subset for (auto &&[name, s] : this->subsets) { if (mask & s->mask()) - s->unpackData(ref, buf); + s->unpackData(ref, buf, true); } + SET_REF_INTERNAL_TRACE_INDEX(ref, ref.getIndex()); }; - - typename ReferenceManager>::MapRefToUInt masks; }; __END_LIBMULTISCALE__ #endif /* __LIBMULTISCALE_REFERENCE_MANAGER_LAMMPS_HH__ */ diff --git a/src/md/trace_atom.cc b/src/md/trace_atom.cc index 4b3b401..28bb2a5 100644 --- a/src/md/trace_atom.cc +++ b/src/md/trace_atom.cc @@ -1,48 +1,48 @@ /** * @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 . * */ #include "trace_atom.hh" #ifdef LIBMULTISCALE_TRACE_ATOM #include "lm_common.hh" #include /* -------------------------------------------------------------------------- */ __BEGIN_LIBMULTISCALE__ typedef unsigned int UInt; typedef double Real; -Real trace_ref[3] = {36.6912, 38.7072, -1.6128}; +Real trace_ref[3] = {-3.6288, 26.6112, -1.6128}; 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 6a59a7f..f9c4e27 100644 --- a/src/md/trace_atom.hh +++ b/src/md/trace_atom.hh @@ -1,116 +1,122 @@ /** * @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__ /* -------------------------------------------------------------------------- */ #define LIBMULTISCALE_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; /* -------------------------------------------------------------------------- */ 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); } /* -------------------------------------------------------------------------- */ template bool CHECK_REF_TRACED(Ref &&ref) { return CHECK_TRACED(ref.position0()); } /* -------------------------------------------------------------------------- */ __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) \ if (CHECK_TRACED(X)) { \ DUMP(ATOM_OUT(X) << message, DBG_MESSAGE); \ } /* -------------------------------------------------------------------------- */ #define IF_REF_TRACED(ref, 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) \ if (CHECK_TRACED(X)) \ internal_index = index; /* -------------------------------------------------------------------------- */ #define SET_REF_INTERNAL_TRACE_INDEX(ref, index) \ { \ if (CHECK_REF_TRACED(ref)) \ internal_index = index; \ } -/* -------------------------------------------------------------------------- - */ - +/* ----------------------------------------------------------------------- */ +#define TRACE_BREAKPOINT(ref) \ + if (CHECK_REF_TRACED(ref)) { \ + int a = 1; \ + while (a) { \ + } \ + } +/* ----------------------------------------------------------------------- */ #else #define IF_TRACED(X, message) #define IF_REF_TRACED(ref, message) #define VIEW_ATOM(T) #define CHECK_TRACED(X) (1 == 2) #define CHECK_REF_TRACED(X) (1 == 2) #define IF_COORD_EQUALS(x1, x2) #define SET_INTERNAL_TRACE_INDEX(X, index) #define SET_REF_INTERNAL_TRACE_INDEX(ref, index) +#define TRACE_BREAKPOINT(ref) #endif /* -------------------------------------------------------------------------- */ #endif /* __LIBMULTISCALE_TRACE_ATOM_H__ */ diff --git a/test/migration/test_migration.cc b/test/migration/test_migration.cc index e8f0d69..1952591 100644 --- a/test/migration/test_migration.cc +++ b/test/migration/test_migration.cc @@ -1,227 +1,228 @@ /* ./test/migration/test_migration.cc ********************************** author : Guillaume ANCIAUX (guillaume.anciaux@epfl.ch, g.anciaux@laposte.net) 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 is done in a collaboration between EPFL within ENAC-LSMS (http://lsms.epfl.ch/) and INRIA Bordeaux, ScAlApplix (http://www.labri.fr/projet/scalapplix/). 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 TIMER #define MAIN_SOURCE #include "algebraic_parser.hh" #include "domain_multiscale.hh" #include "factory_multiscale.hh" #include "lib_md.hh" #include "lm_common.hh" #include "reference_manager_interface.hh" #include #include #include /* -------------------------------------------------------------------------- */ using namespace libmultiscale; /* -------------------------------------------------------------------------- */ Real GaussianRand() { Real x1, x2, w; // UInt *seed = &randSeed; do { x1 = 2.0 * drand48() - 1.0; x2 = 2.0 * drand48() - 1.0; // x1 = RANDOM_abs1( seed ); // x2 = RANDOM_abs1( seed ); w = x1 * x1 + x2 * x2; } while (w >= 1.0); w = sqrt((-2.0 * log(w)) / w); return x1 * w; } /* -------------------------------------------------------------------------- */ int main(int argc, char **argv) { loadModules(argc, argv); // loading the modules macro #ifdef LIBMULTISCALE_USE_LAMMPS MPI_Barrier(MPI_COMM_WORLD); DomainMultiScale &dom = DomainMultiScale::getManager(); DOWAIT_AT_STARTUP; dom.build(argv[1]); STOPTIMER("Init"); Real time_step = 1e300; dom.for_each([&time_step](auto &d) { Real ts = d.template getParam>("TIMESTEP"); time_step = std::min(time_step, ts); }); DUMP("Multiscale timestep:" << time_step, DBG_MESSAGE); DomainLammps<3> &dL = static_cast &>(*dom.getModelByNumber(0)); ContainerLammps<3> &c = dL.getContainer(); ContainerArray> subset("everything-subset"); dL.getContainer().getRefManager()->addSubset(subset); // build the subset to be attached to reference manager for (auto at : c) { subset.push_back(at); } ContainerArray> subset2("half-subset"); dL.getContainer().getRefManager()->addSubset(subset2); // build the subset to be attached to reference manager UInt cpt = 0; for (auto at : c) { if (cpt % 2 == 0) subset2.push_back(at); ++cpt; } ContainerArray v_p0("v_p0", 0, 3); ContainerArray v2_p0("v2_p0", 0, 2); ContainerArray v3_p0("v3_p0", 0, 2); // build a vector to be migrated with the atoms cpt = 0; for (auto at : c) { Vector<3> pos0 = at.position0(); v_p0.push_back(pos0); Vector<2> pos1 = pos0.block<2, 1>(0, 0) * 2; v2_p0.push_back(pos1); if (cpt % 2 == 0) { v3_p0.push_back(pos1); } ++cpt; } dL.getContainer().getRefManager()->attachObject(v_p0, subset); dL.getContainer().getRefManager()->attachObject(v2_p0, subset); dL.getContainer().getRefManager()->attachObject(v3_p0, subset2); srand48(345); // force random migrations for (auto at : c) { at.position()[0] += GaussianRand(); at.position()[1] += GaussianRand(); at.position()[2] += GaussianRand(); } // UInt index_pb = 4; // if (lm_my_proc_id == 0){ // at = subset.Get(index_pb); // DUMP("atome " << index_pb << " coords " // << at.position0(0) << " " // << at.position0(1) << " " // << at.position0(2) << " " // << v_p0[3*index_pb],DBG_ERROR); // } dom.coupling(COUPLING_STEP3); dom.for_each([time_step](auto &d) { auto &&v = d.primalTimeDerivative(); auto &&p = d.primal(); auto &&acc = d.acceleration(); v += 0.5 * time_step * acc; p += time_step * v; d.changeRelease(); d.enforceCompatibility(); }); ++current_step; dom.coupling(COUPLING_STEP3); // if (lm_my_proc_id == 0){ // at = subset.Get(index_pb); // DUMP("atome " << index_pb << " coords " // << at.position0(0) << " " // << at.position0(1) << " " // << at.position0(2) << " " // << v_p0[3*index_pb],DBG_ERROR); // } MPI_Barrier(MPI_COMM_WORLD); // test that all attached initial positions are still correct in v_p0 + Eigen::IOFormat fmt(Eigen::FullPrecision, 0, ", ", "", "", "", "[", "]"); if (lm_my_proc_id == 0) { for (auto [i, at] : enumerate(subset)) { if (v_p0(i, 0) != at.position0()[0]) - LM_FATAL("fail on x axis " << i << " " << v_p0(i, 0) << " " - << at.position0()[0] << " " << at << " " - << v_p0.row(i)); + LM_FATAL("fail@x " << subset.getID() << "[" << i << "] " + << v_p0.row(i).format(fmt) + << " != " << at.position0().format(fmt)); if (v_p0(i, 1) != at.position0()[1]) - LM_FATAL("fail on y axis " << i << " " << v_p0(i, 1) << " " - << at.position0()[1] << " " << at << " " - << v_p0.row(i)); + LM_FATAL("fail@y " << subset.getID() << "[" << i << "] " + << v_p0.row(i).format(fmt) + << " != " << at.position0().format(fmt)); if (v_p0(i, 2) != at.position0()[2]) - LM_FATAL("fail on z axis " << i << " " << v_p0(i, 2) << " " - << at.position0()[2] << " " << at << " " - << v_p0.row(i)); + LM_FATAL("fail@z " << subset.getID() << "[" << i << "] " + << v_p0.row(i).format(fmt) + << " != " << at.position0().format(fmt)); if (v2_p0(i, 0) != at.position0()[0] * 2) - LM_FATAL("fail on x axis " << i << " " << v2_p0(i, 0) << " " - << at.position0()[0] * 2 << " " << at << " " - << v_p0.row(i)); + LM_FATAL("fail@x " << subset.getID() << "[" << i << "] " + << v2_p0.row(i).format(fmt) + << " != " << (at.position0() * 2).format(fmt)); if (v2_p0(i, 1) != at.position0()[1] * 2) - LM_FATAL("fail on y axis " << i << " " << v2_p0(i, 1) << " " - << at.position0()[1] * 2 << " " << at << " " - << v_p0.row(i)); + LM_FATAL("fail@y " << subset.getID() << "[" << i << "] " + << v2_p0.row(i).format(fmt) + << " != " << (at.position0() * 2).format(fmt)); } for (auto [i, at] : enumerate(subset2)) { if (v3_p0(i, 0) != at.position0()[0] * 2) - LM_FATAL("fail on x axis " << i << " " << v2_p0(i, 0) << " " - << at.position0()[0] * 2 << " " << at << " " - << v_p0.row(i)); + LM_FATAL("fail@ x " << subset.getID() << "[" << i << "] " + << v2_p0.row(i).format(fmt) + << " != " << (at.position0() * 2).format(fmt)); if (v3_p0(i, 1) != at.position0()[1] * 2) - LM_FATAL("fail on y axis " << i << " " << v2_p0(i, 1) << " " - << at.position0()[1] * 2 << " " << at << " " - << v_p0.row(i)); + LM_FATAL("fail@y " << subset.getID() << "[" << i << "] " + << v2_p0.row(i).format(fmt) + << " != " << (at.position0()).format(fmt)); } } #endif closeModules(); // closing the modules macro return 0; }