diff --git a/src/common/ref_subset.cc b/src/common/ref_subset.cc index a5aeda7..13f1003 100644 --- a/src/common/ref_subset.cc +++ b/src/common/ref_subset.cc @@ -1,464 +1,474 @@ /** * @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); + this->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(); + if (!this->holes.empty()) { + new_index = this->holes.back(); container.get(new_index) = ref; - holes.pop_back(); + this->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) { att_obj->packData(index, buffer, verbose); if (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) { auto index = this->getRefIndex(ref); if constexpr (is_subset) { if (index != -1) LM_FATAL("this should not happen: found " << container.getID() << "[" << index << "] = " << ref); - if (holes.size() > 0) { - index = holes.back(); + if (this->holes.size() > 0) { + index = this->holes.back(); this->container.get(index) = ref; - holes.pop_back(); + this->holes.pop_back(); } else { index = this->container.size(); this->container.resize(index + 1); } } for (auto &&[ID, att_obj] : this->attached_objects) { att_obj->unpackData(index, buffer, verbose); if (CHECK_REF_TRACED(ref)) { 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; // 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) { + DUMP("could not find " << src << " in " << this->container.getID(), + DBG_MESSAGE); return; } - 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); - } + // 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 // if it is in this subset then I have a hole to register container.get(i_src) = dst; if (i_dst != -1) { DUMP("adding hole: " << i_dst, DBG_INFO); - holes.push_back(i_dst); + this->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); - } + // 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(); + while (!this->holes.empty()) { + UInt index = this->holes.back(); + this->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), " = " << container.getID() << "[" << i_src << "]: " << " traced ref moved from index " << i_src << " to index " << i_dest << " for container " << container.getID()); 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 +bool RefSubset::is_hole(UInt index) { + auto it = std::find(this->holes.begin(), this->holes.end(), index); + bool is_hole(it != this->holes.end()); + return is_hole; +} /* ----------------------------------------------------------------------- */ template int RefSubset::getRefIndex(Ref &ref) { for (auto &&[i, e] : enumerate(container)) { - if (e == ref) + if ((not this->is_hole(i)) and 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 c4bfd5e..34cfb17 100644 --- a/src/common/ref_subset.hh +++ b/src/common/ref_subset.hh @@ -1,205 +1,216 @@ /** * @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; + template Cont &getContainer() { + return dynamic_cast(this->getContainer()); + } //! 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; + +public: + // list of holes created by the gone refs + std::vector holes; + virtual int getRefIndex(Ref &ref) = 0; + + virtual bool is_hole(UInt index) = 0; + virtual std::string getID() = 0; }; /* -------------------------------------------------------------------------- */ 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); + bool is_hole(UInt index); // //! 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); + std::string getID() { return container.getID(); }; 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/md/lammps/reference_manager_lammps.cc b/src/md/lammps/reference_manager_lammps.cc index 8da3497..d4142e7 100644 --- a/src/md/lammps/reference_manager_lammps.cc +++ b/src/md/lammps/reference_manager_lammps.cc @@ -1,188 +1,192 @@ /** * @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); // 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)"); + checkAttachedValue("before moving ", refold, this->subsets); for (auto &&[name, s] : this->subsets) { s->translateMovingReference(refold, refnew); } + checkAttachedValue("after moving ", refnew, this->subsets); } /* -------------------------------------------------------------------------- */ // 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 c99e324..e451333 100644 --- a/src/md/lammps/reference_manager_lammps.hh +++ b/src/md/lammps/reference_manager_lammps.hh @@ -1,118 +1,124 @@ /** * @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); + checkAttachedValue("before packing", ref, this->subsets); + IF_REF_TRACED(ref, "sending atom(" << ref.getIndex() << ")"); // first of all pack the masks 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); } SET_REF_INTERNAL_TRACE_INDEX(ref, -1); + checkAttachedValue("after packing", ref, this->subsets); } void unpack_subset_data(UInt i, LMPPackBuffer buf) { // build a reflammps RefLammps ref(i); + // checkAttachedValue("before unpacking", ref, this->subsets); + 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); } SET_REF_INTERNAL_TRACE_INDEX(ref, ref.getIndex()); + checkAttachedValue("after unpacking", ref, this->subsets); }; }; __END_LIBMULTISCALE__ #endif /* __LIBMULTISCALE_REFERENCE_MANAGER_LAMMPS_HH__ */ diff --git a/src/md/trace_atom.hh b/src/md/trace_atom.hh index f9c4e27..0ada948 100644 --- a/src/md/trace_atom.hh +++ b/src/md/trace_atom.hh @@ -1,122 +1,154 @@ /** * @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 "attached_object.hh" #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 false; 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()); } /* -------------------------------------------------------------------------- */ +template class ContainerArray; +template class RefLammps; + +template +void checkAttachedValue(const std::string &mess, Ref &ref, T &subsets) { + auto &s = subsets["half-subset"]; + auto &container = s->template getContainer>>(); + int index = s->getRefIndex(ref); + if (index == -1) + return; + auto &att_obj = + dynamic_cast &>(*s->attached_objects["v3_p0"]); + + bool is_hole = s->is_hole(index); + + DUMP(container.getID() << "[" << index << "] = " << ref << " attached:" + << att_obj.getID() << ": " << att_obj.getV().row(index) + << " hole=" << is_hole << " ### " << mess, + DBG_MESSAGE); + if ((not is_hole)) { + if constexpr (Ref::Dim == 3) { + Vector<3> pos = ref.position0() * 2; + Vector<2> data = att_obj.getV().row(index); + if ((data - pos.block<2, 1>(0, 0)).norm() > 1e-5) + LM_FATAL("aaaaa"); + } + } +} __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 5f2deaa..217f9b9 100644 --- a/test/migration/test_migration.cc +++ b/test/migration/test_migration.cc @@ -1,228 +1,229 @@ /* ./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@x " << subset.getID() << "[" << i << "] = data" << v_p0.row(i).format(fmt) << " != ref" << at.position0().format(fmt)); if (v_p0(i, 1) != at.position0()[1]) LM_FATAL("fail@y " << subset.getID() << "[" << i << "] = data" << v_p0.row(i).format(fmt) << " != ref" << at.position0().format(fmt)); if (v_p0(i, 2) != at.position0()[2]) LM_FATAL("fail@z " << subset.getID() << "[" << i << "] = data" << v_p0.row(i).format(fmt) << " != ref" << at.position0().format(fmt)); if (v2_p0(i, 0) != at.position0()[0] * 2) LM_FATAL("fail@x " << subset.getID() << "[" << i << "] = data" << v2_p0.row(i).format(fmt) << " != ref" << (at.position0() * 2).format(fmt)); if (v2_p0(i, 1) != at.position0()[1] * 2) LM_FATAL("fail@y " << subset.getID() << "[" << i << "] = data" << v2_p0.row(i).format(fmt) << " != ref" << (at.position0() * 2).format(fmt)); } for (auto [i, at] : enumerate(subset2)) { if (v3_p0(i, 0) != at.position0()[0] * 2) LM_FATAL("fail@ x " << subset.getID() << "[" << i << "] = data" - << v2_p0.row(i).format(fmt) << " != ref" + << v3_p0.row(i).format(fmt) << " != ref" << (at.position0() * 2).format(fmt)); if (v3_p0(i, 1) != at.position0()[1] * 2) LM_FATAL("fail@y " << subset.getID() << "[" << i << "] = data" - << v2_p0.row(i).format(fmt) << " != ref" + << v3_p0.row(i).format(fmt) << " != ref" << (at.position0()).format(fmt)); } } #endif closeModules(); // closing the modules macro return 0; }