diff --git a/src/common/lm_log.cc b/src/common/lm_log.cc
index 49b55b2..f3cc804 100644
--- a/src/common/lm_log.cc
+++ b/src/common/lm_log.cc
@@ -1,175 +1,368 @@
/**
* @file lm_log.cc
*
* @author Guillaume Anciaux <guillaume.anciaux@epfl.ch>
*
* @date Mon Sep 08 23:40:22 2014
*
* @brief Log/Debug system of LibMultiScale
*
* @section LICENSE
*
* Copyright (©) 2010-2011 EPFL (Ecole Polytechnique Fédérale de Lausanne)
* Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)
*
* LibMultiScale is free software: you can redistribute it and/or modify it
* under the
* terms of the GNU Lesser General Public License as published by the Free
* Software Foundation, either version 3 of the License, or (at your option) any
* later version.
*
* LibMultiScale is distributed in the hope that it will be useful, but
* WITHOUT ANY
* WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
* A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more
* details.
*
* You should have received a copy of the GNU Lesser General Public License
* along with LibMultiScale. If not, see <http://www.gnu.org/licenses/>.
*
*/
#include "action_interface.hh"
#include "factory_multiscale.hh"
#include "lm_common.hh"
#include "lm_globals.hh"
+#include <cxxabi.h>
+#include <execinfo.h>
#include <fstream>
+#include <iomanip>
#include <iostream>
#include <mpi.h>
/* -------------------------------------------------------------------------- */
__BEGIN_LIBMULTISCALE__
/* -------------------------------------------------------------------------- */
#define CREATE_SEGFAULT \
if (create_seg_fault) { \
UInt *a = NULL; \
*a = 1; \
}
/* -------------------------------------------------------------------------- */
void doWait(bool condition) {
if (!condition)
return;
UInt pid = lm_getpid();
std::cerr << "Proc " << lm_my_proc_id << " Waiting for gdb to attach: pid is "
<< pid << std::endl;
UInt a = 1;
while (a) {
};
}
/* -------------------------------------------------------------------------- */
void doWaitOnFatal() { doWait(wait_on_fatal); }
/* -------------------------------------------------------------------------- */
void doWaitAtStartup() {
doWait((wait_at_startup != UINT_MAX) && (lm_my_proc_id == wait_at_startup));
}
/* -------------------------------------------------------------------------- */
+std::string exec(const std::string &cmd) {
+ FILE *pipe = popen(cmd.c_str(), "r");
+ if (!pipe)
+ return "";
+ char buffer[1024];
+ std::string result = "";
+ while (!feof(pipe)) {
+ if (fgets(buffer, 128, pipe) != nullptr)
+ result += buffer;
+ }
+
+ result = result.substr(0, result.size() - 1);
+ pclose(pipe);
+ return result;
+}
+
+/* -------------------------------------------------------------------------- */
+
+/** Print a demangled stack backtrace of the caller function to a string */
+static inline std::string print_stacktrace(UInt max_frames = 10) {
+
+ std::stringstream mesg_stream;
+ mesg_stream << "\n==== C++ stack trace ==== \n";
+
+ std::string me = "";
+ char buf[1024];
+ /* The manpage says it won't null terminate. Let's zero the buffer. */
+ memset(buf, 0, sizeof(buf));
+ /* Note we use sizeof(buf)-1 since we may need an extra char for NUL. */
+ if (readlink("/proc/self/exe", buf, sizeof(buf) - 1))
+ me = std::string(buf);
+
+ std::ifstream inmaps;
+ inmaps.open("/proc/self/maps");
+ std::map<std::string, size_t> addr_map;
+ std::string line;
+ while (inmaps.good()) {
+ std::getline(inmaps, line);
+ std::stringstream sstr(line);
+
+ size_t first = line.find('-');
+ std::stringstream sstra(line.substr(0, first));
+ size_t addr;
+ sstra >> std::hex >> addr;
+
+ std::string lib;
+ sstr >> lib;
+ sstr >> lib;
+ sstr >> lib;
+ sstr >> lib;
+ sstr >> lib;
+ sstr >> lib;
+ if (lib != "" && addr_map.find(lib) == addr_map.end()) {
+ addr_map[lib] = addr;
+ }
+ }
+
+ if (me != "")
+ addr_map[me] = 0;
+
+ size_t stack_depth;
+ void *stack_addrs[max_frames];
+ char **stack_strings;
+
+ stack_depth = backtrace(stack_addrs, max_frames);
+ stack_strings = backtrace_symbols(stack_addrs, stack_depth);
+
+ //auto w = size_t(std::floor(log(double(stack_depth)) / std::log(10.)) + 1);
+
+ for (UInt i = stack_depth - 1; i > 0; i--) {
+ //mesg_stream << std::dec << " [" << std::setw(w) << i << "] ";
+ std::string bt_line(stack_strings[i]);
+ size_t first, second;
+
+ if ((first = bt_line.find('(')) != std::string::npos &&
+ (second = bt_line.find('+')) != std::string::npos) {
+ std::string location = bt_line.substr(0, first);
+#if defined(READLINK_COMMAND)
+ std::string location_cmd =
+ std::string("/usr/bin/readlink") + std::string(" -f ") + location;
+ location = exec(location_cmd);
+#endif
+ std::string call =
+ demangle(bt_line.substr(first + 1, second - first - 1).c_str());
+ size_t f = bt_line.find('[');
+ size_t s = bt_line.find(']');
+ std::string address = bt_line.substr(f + 1, s - f - 1);
+ std::stringstream sstra(address);
+ size_t addr;
+ sstra >> std::hex >> addr;
+
+ //mesg_stream << location << " [" << call << "]";
+
+ auto it = addr_map.find(location);
+ if (it != addr_map.end()) {
+ std::stringstream syscom;
+ syscom << "/usr/bin/addr2line"
+ << " 0x" << std::hex << (addr - it->second) << " -i -e "
+ << location;
+ std::string line = exec(syscom.str());
+ std::replace(line.begin(), line.end(), ' ', ':');
+ mesg_stream << line << ": " << std::endl;
+ } else {
+ mesg_stream << location << ":0x" << std::hex << addr << std::endl;
+ }
+ } else {
+ mesg_stream << bt_line << std::endl;
+ }
+ }
+
+ free(stack_strings);
+
+ return mesg_stream.str();
+}
+
+/* -------------------------------------------------------------------------- */
+/** Print a demangled stack backtrace of the caller function to a string */
+static inline std::string print_stacktrace2(UInt max_frames = 10) {
+
+ std::string mesg;
+ mesg += "\n==== C++ stack trace ==== \n";
+
+ // storage array for stack trace address data
+ std::vector<void *> addrlist(max_frames);
+ // retrieve current stack addresses
+ int addrlen = backtrace(addrlist.data(), addrlist.size());
+
+ if (addrlen == 0) {
+ mesg += " <empty, possibly corrupt>\n";
+ return mesg;
+ }
+
+ addrlist.resize(addrlen);
+
+ // resolve addresses into strings containing "filename(function+address)",
+ // this array must be free()-ed
+ auto symbollist = backtrace_symbols(addrlist.data(), addrlist.size());
+
+ // iterate over the returned symbol lines. skip the first, it is the
+ // address of this function.
+ for (UInt i = addrlist.size() - 1; i > 0; i--) {
+ char *begin_name = 0, *begin_offset = 0, *end_offset = 0;
+
+ // find parentheses and +address offset surrounding the mangled name:
+ // ./module(function+0x15c) [0x8048a6d]
+ for (char *p = symbollist[i]; *p; ++p) {
+ if (*p == '(')
+ begin_name = p;
+ else if (*p == '+')
+ begin_offset = p;
+ else if (*p == ')' && begin_offset) {
+ end_offset = p;
+ break;
+ }
+ }
+
+ if (begin_name && begin_offset && end_offset && begin_name < begin_offset) {
+ *begin_name++ = '\0';
+ *begin_offset++ = '\0';
+ *end_offset = '\0';
+
+ // mangled name is now in [begin_name, begin_offset) and caller
+ // offset in [begin_offset, end_offset). now apply
+ // __cxa_demangle():
+
+ int status;
+ decltype(auto) ret = abi::__cxa_demangle(begin_name, NULL, NULL, &status);
+
+ std::string funcname;
+
+ if (status == 0) {
+ funcname = ret;
+ free(ret);
+ } else
+ funcname = begin_name;
+
+ mesg += symbollist[i];
+ mesg += " : " + funcname + "+" + begin_offset + "\n";
+
+ } else {
+ // couldn't parse the line? print the whole line.
+ mesg += symbollist[i];
+ mesg += "\n";
+ }
+ }
+ free(symbollist);
+ return mesg;
+}
+
+/* -------------------------------------------------------------------------- */
+
void lmFatal(std::stringstream &x, UInt line, const std::string &file,
const std::string &function, const std::string &pretty_function) {
if (!fatal_loop) {
fatal_loop = true;
try {
auto &actions = ActionManager::getManager();
current_stage = PRE_FATAL;
actions.action();
} catch (...) {
}
}
std::stringstream sstr;
- sstr << "\n"
+ sstr << print_stacktrace() << "\n\n"
<< file << ":" << line << ":FATAL: (" << current_step << ")"
<< "[" << lm_my_proc_id << "]:" << x.str() << std::endl;
(*global_out) << sstr.str();
(*global_out).flush();
CREATE_SEGFAULT;
DOWAIT_ON_FATAL;
throw LibMultiScaleException(sstr.str());
CREATE_SEGFAULT;
// ::libmultiscale::lm_exit(LM_EXIT_FAILURE);
}
/* -------------------------------------------------------------------------- */
void openGlobal() {
if (print_crap_to_file) {
std::stringstream sstr;
sstr << my_hostname << "-" << lm_my_proc_id << ".log";
file_out.open(sstr.str().c_str(), std::ios_base::out);
global_out = &file_out;
} else {
global_out = &std::cerr;
}
}
/* -------------------------------------------------------------------------- */
inline std::ostream &operator<<(std::ostream &stream, const dbgLevel &_this) {
switch (_this) {
case dbg_message:;
break;
case dbg_warning:
stream << "WARNING:";
break;
case dbg_info_startup:
stream << "INFO_STARTUP:";
break;
case dbg_info:
stream << "INFO:";
break;
case dbg_detail:
stream << "DETAIL:";
break;
case dbg_all:
stream << "ALL:";
break;
}
return stream;
}
/* -------------------------------------------------------------------------- */
#define INFO_PROCESS \
"s=" << current_step << ":" << current_stage << ":proc=" << lm_my_proc_id \
<< ": "
#define INFO_PROCESS_FULL \
"s=" << current_step << ":" << current_stage << ":proc=" << lm_my_proc_id \
<< ":host=" << my_hostname << "," << lm_getpid() << " : "
#define INFO_SOURCE function << ":" << line << ":"
-#define INFO_SOURCE_FULL \
- file << ":" << line << ":" << pretty_function << ":"
+#define INFO_SOURCE_FULL file << ":" << line << ":" << pretty_function << ":"
/* ----------------- */
// no info
//#define APPEND_FILE ""
/* ------------------*/
// small info
#define APPEND_FILE INFO_SOURCE << INFO_PROCESS
/* ------------------*/
// full info
//#define APPEND_FILE INFO_SOURCE_FULL << INFO_PROCESS_FULL
/* -------------------------------------------------------------------------- */
void lmPrintOut(std::stringstream &x, dbgLevel level, UInt line,
const std::string &file, const std::string &function,
const std::string &pretty_function) {
if (global_level >= level)
if ((global_proc == lm_uint_max || global_proc == lm_my_proc_id ||
global_proc1 == lm_my_proc_id || global_proc2 == lm_my_proc_id) &&
(current_step >= start_dump &&
(end_dump == lm_uint_max || current_step <= end_dump))) {
if (global_out)
(*global_out) << APPEND_FILE << x.str() << std::endl << std::flush;
}
}
/* -------------------------------------------------------------------------- */
__END_LIBMULTISCALE__
diff --git a/src/container/container_array.hh b/src/container/container_array.hh
index 2508580..a22d34a 100644
--- a/src/container/container_array.hh
+++ b/src/container/container_array.hh
@@ -1,368 +1,379 @@
/**
* @file container_array.hh
*
* @author Guillaume Anciaux <guillaume.anciaux@epfl.ch>
*
* @date Thu Jul 24 14:21:58 2014
*
* @brief This is the contiguous array-like container for LM
*
* @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_CONTAINER_ARRAY_HH__
#define __LIBMULTISCALE_CONTAINER_ARRAY_HH__
/* -------------------------------------------------------------------------- */
#include "container.hh"
#include <algorithm>
#include <vector>
template <typename Ref> class ReferenceManager;
/* -------------------------------------------------------------------------- */
#ifdef TRACE_ATOM
#include "lammps/ref_atom_lammps.hh"
#include "trace_atom.hh"
#endif
/* -------------------------------------------------------------------------- */
__BEGIN_LIBMULTISCALE__
/**
* Class ContainerArray<T>
* array implementation of a generic container
*/
/* -------------------------------------------------------------------------- */
template <class T> struct GetDim {
static constexpr UInt Dim = T::Dim;
inline UInt getDim() { return T::Dim; }
};
template <> struct GetDim<Real> {};
template <> struct GetDim<UInt> {};
/* -------------------------------------------------------------------------- */
template <typename T> class ContainerArrayIterator;
template <typename T> class ContainerArray;
template <typename T>
std::enable_if_t<!std::is_arithmetic<T>::value, UInt>
getContainerArrayDim(ContainerArray<T> &cont) {
return T::Dim;
};
template <typename T>
std::enable_if_t<std::is_arithmetic<T>::value, UInt>
getContainerArrayDim(ContainerArray<T> &cont) {
return cont.cols();
};
/* -------------------------------------------------------------------------- */
template <class T>
class ContainerArray : public Container_base<T>,
public GetDim<T>,
public virtual LMObject {
/* ------------------------------------------------------------------------ */
/* Typedefs */
/* ------------------------------------------------------------------------ */
public:
using ContainerSubSet = ContainerArray<T>;
using iterator = ContainerArrayIterator<T>;
using EigenArray = Eigen::Array<T, Eigen::Dynamic, Eigen::Dynamic>;
/* ------------------------------------------------------------------------ */
/* Constructors/Destructors */
/* ------------------------------------------------------------------------ */
ContainerArray(const LMID &id);
ContainerArray();
ContainerArray(UInt size);
ContainerArray(ContainerArray &&);
ContainerArray(ContainerArray &);
ContainerArray &operator=(const ContainerArray &);
ContainerArray &operator=(const EigenArray &);
~ContainerArray();
/* ------------------------------------------------------------------------ */
/* Methods */
/* ------------------------------------------------------------------------ */
UInt getDim() { return getContainerArrayDim(*this); }
//! return an associated iterator
virtual iterator begin(DOFType dt = dt_local);
virtual iterator end(DOFType dt = dt_local);
//! sort items
virtual void sort() override;
//! search for a given item
virtual UInt search(T &el) override;
//! add a particular item
virtual void push_back(const T &el) override;
+ //! add a particular item
+ template <UInt Dim> void push_back(const Array1D<Dim> &el);
// //! search for a given item in a previously sorted array
// UInt searchInSorted(T &el);
// //! add an entry in a sorted array
// void addInSorted(T &el);
//! resize the container (by clearing the additional values)
virtual void resize(UInt size);
//! assign all values up to size to a single value
virtual void assign(UInt size, T val);
//! get an item from its index
inline T &get(UInt index) override;
//! get an item from its index
inline T &operator[](UInt index) override;
//! return the number of contained items
virtual UInt size(DOFType dt = dt_local) override;
//! return the number of columns in the array
UInt cols() { return array->cols(); };
//! function that clear the container
virtual void clear();
// ! copy the container info
virtual void copyContainerInfo(ContainerInterface &cont) override;
// //! get an item from its computed name
// virtual Real &operator()(const std::string &name, UInt index = 0) {
// std::vector<std::string>::iterator it;
// it = find(name_computed.begin(), name_computed.end(), name);
// if (it == name_computed.end()) {
// LM_FATAL("Name computed '" << name << "' not found");
// } else {
// UInt position = it - name_computed.begin();
// LM_FATAL("temp " << position);
// }
// }
// virtual UInt getDim() { return name_computed.size(); };
// std::vector<std::string> name_computed;
private:
EigenArray array_storage;
std::unique_ptr<Eigen::Map<EigenArray>> array;
};
/* -------------------------------------------------------------------------- */
__END_LIBMULTISCALE__
#include "iterator_array.hh"
#include "reference_manager.hh"
__BEGIN_LIBMULTISCALE__
/* -------------------------------------------------------------------------- */
template <typename T>
ContainerArray<T>::ContainerArray(const LMID &id)
: LMObject(id), Container_base<T>() {
array = std::make_unique<Eigen::Map<EigenArray>>(array_storage.data(), 0, 1);
}
/* -------------------------------------------------------------------------- */
template <typename T>
ContainerArray<T>::ContainerArray()
: LMObject("anonymous"), Container_base<T>() {
array = std::make_unique<Eigen::Map<EigenArray>>(array_storage.data(), 0, 1);
}
/* -------------------------------------------------------------------------- */
template <typename T>
ContainerArray<T>::ContainerArray(UInt size)
: LMObject("anonymous"), Container_base<T>() {
array =
std::make_unique<Eigen::Map<EigenArray>>(array_storage.data(), size, 1);
}
/* -------------------------------------------------------------------------- */
template <typename T>
ContainerArray<T>::ContainerArray(ContainerArray &&other)
: LMObject(other.getID()) {
this->array_storage = other.array_storage;
array = std::make_unique<Eigen::Map<EigenArray>>(
array_storage.data(), other.array->rows(), other.array->cols());
}
/* -------------------------------------------------------------------------- */
template <typename T>
ContainerArray<T>::ContainerArray(ContainerArray &other)
: LMObject(other.getID()) {
this->array_storage = other.array_storage;
array = std::make_unique<Eigen::Map<EigenArray>>(
array_storage.data(), other.array->rows(), other.array->cols());
}
/* -------------------------------------------------------------------------- */
template <typename T>
ContainerArray<T> &ContainerArray<T>::operator=(const ContainerArray &other) {
this->array_storage = other.array_storage;
array = std::make_unique<Eigen::Map<EigenArray>>(
array_storage.data(), other.array->rows(), other.array->cols());
return *this;
}
/* -------------------------------------------------------------------------- */
template <typename T>
ContainerArray<T> &ContainerArray<T>::operator=(const EigenArray &other) {
this->array_storage = other;
array = std::make_unique<Eigen::Map<EigenArray>>(array_storage.data(),
other.rows(), other.cols());
return *this;
}
/* -------------------------------------------------------------------------- */
template <typename T>
inline void
copy_ref_manager(std::weak_ptr<ReferenceManagerInterface> refmanager,
std::weak_ptr<ReferenceManagerInterface> other_refmanager) {
try {
__attribute__((unused)) auto &test =
dynamic_cast<ReferenceManager<T> &>(*other_refmanager.lock());
refmanager = other_refmanager;
} catch (...) {
}
}
template <>
inline void copy_ref_manager<Real>(
std::weak_ptr<ReferenceManagerInterface> refmanager,
std::weak_ptr<ReferenceManagerInterface> other_refmanager) {}
template <typename T>
void ContainerArray<T>::copyContainerInfo(ContainerInterface &cont) {
ContainerInterface::copyContainerInfo(cont);
if (cont.hasRefManager()) {
copy_ref_manager<T>(this->refmanager, cont.getRefManager());
}
}
/* -------------------------------------------------------------------------- */
template <typename T> ContainerArray<T>::~ContainerArray() {
clear();
if (this->hasRefManager()) {
this->getRefManager()->removeSubSet(this->getID(), *this);
}
}
/* ------------------------------------------------------------------------ */
template <class T> void ContainerArray<T>::sort() {
// std::sort(array.begin(),array.end());
LM_TOIMPLEMENT;
}
/* ------------------------------------------------------------------------ */
template <class T> inline UInt ContainerArray<T>::size(DOFType dt) {
return array->size();
}
/* ------------------------------------------------------------------------ */
template <class T> inline void ContainerArray<T>::clear() {
auto cols = array_storage.cols();
array_storage.resize(0, cols);
array =
std::make_unique<Eigen::Map<EigenArray>>(array_storage.data(), 0, cols);
}
/* ------------------------------------------------------------------------ */
template <class T> inline void ContainerArray<T>::resize(UInt size) {
UInt ncols = array_storage.cols();
if (ncols == 0)
ncols = 1;
if (size % ncols != 0)
LM_FATAL("invalid size");
UInt nrows = size / ncols;
while (nrows > array_storage.rows()) {
array_storage.conservativeResize(1 + array_storage.rows() * 2, ncols);
}
array = std::make_unique<Eigen::Map<EigenArray>>(array_storage.data(), nrows,
ncols);
}
/* ------------------------------------------------------------------------ */
template <class T> inline void ContainerArray<T>::assign(UInt size, T val) {
this->resize(size);
*array = val;
}
/* ------------------------------------------------------------------------ */
template <class T> inline void ContainerArray<T>::push_back(const T &el) {
auto row = array->rows();
this->resize(row + 1);
array->row(row) = el;
}
/* ------------------------------------------------------------------------ */
+template <class T>
+template <UInt Dim>
+inline void ContainerArray<T>::push_back(const Array1D<Dim> &arr) {
+ for (UInt i = 0; i < Dim; ++i)
+ this->push_back(arr[i]);
+}
+/* ------------------------------------------------------------------------ */
+
template <class T> inline T &ContainerArray<T>::get(UInt index) {
- LM_ASSERT(index < size(),
- "out of range access : " << index << " over " << size());
+ LM_ASSERT(index < size(), "out of range access : " << index << " over "
+ << size() << " for "
+ << this->getID());
using EigenArray1 = Eigen::Array<T, Eigen::Dynamic, 1u>;
Eigen::Map<EigenArray1> view(array->data(), this->size());
return view(index);
}
/* ------------------------------------------------------------------------ */
template <class T> inline T &ContainerArray<T>::operator[](UInt index) {
LM_ASSERT(index < size(),
"out of range access : " << index << " over " << size());
return this->get(index);
}
/* ------------------------------------------------------------------------ */
template <class T> UInt ContainerArray<T>::search(T &el) {
for (UInt i = 0; i < size(); ++i) {
T &e = this->get(i);
if (el == e)
return i;
}
return UINT_MAX;
}
/* ------------------------------------------------------------------------ */
template <class T>
typename ContainerArray<T>::iterator ContainerArray<T>::begin(DOFType dt) {
return iterator(array->data());
}
/* ------------------------------------------------------------------------ */
template <class T>
typename ContainerArray<T>::iterator ContainerArray<T>::end(DOFType dt) {
return iterator(array->data() + this->size());
}
/* ------------------------------------------------------------------------ */
__END_LIBMULTISCALE__
#endif /* __LIBMULTISCALE_CONTAINER_ARRAY_HH__ */
diff --git a/src/filter/compute_ekin.cc b/src/filter/compute_ekin.cc
index cf72f27..2419bf5 100644
--- a/src/filter/compute_ekin.cc
+++ b/src/filter/compute_ekin.cc
@@ -1,233 +1,232 @@
/**
* @file compute_ekin.cc
*
* @author Guillaume Anciaux <guillaume.anciaux@epfl.ch>
*
* @date Wed Jul 09 21:59:47 2014
*
* @brief This computes the kinetic energy of a set of DOFs
*
* @section LICENSE
*
* Copyright (©) 2010-2011 EPFL (Ecole Polytechnique Fédérale de Lausanne)
* Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)
*
* LibMultiScale is free software: you can redistribute it and/or modify it
* under the
* terms of the GNU Lesser General Public License as published by the Free
* Software Foundation, either version 3 of the License, or (at your option) any
* later version.
*
* LibMultiScale is distributed in the hope that it will be useful, but
* WITHOUT ANY
* WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
* A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more
* details.
*
* You should have received a copy of the GNU Lesser General Public License
* along with LibMultiScale. If not, see <http://www.gnu.org/licenses/>.
*
*/
/* -------------------------------------------------------------------------- */
#include "compute_ekin.hh"
#include "communicator.hh"
#include "lib_continuum.hh"
#include "lib_dd.hh"
#include "lib_md.hh"
#include "lm_common.hh"
#include "ref_point_data.hh"
#include "units.hh"
#include "units_converter.hh"
/* -------------------------------------------------------------------------- */
__BEGIN_LIBMULTISCALE__
ComputeEKin::ComputeEKin(const std::string &name)
: LMObject(name), based_on_nodes_flag(false) {
this->createInput("velocities");
auto output_names =
std::list<std::string>{"nb_points", "Temperature", "EKin", "EKin_vector"};
for (auto f : output_names) {
this->createOutput(f);
- this->getOutput("nbpoints") =
+ this->getOutput(f) =
make_argument(ContainerArray<Real>(this->getID() + ":" + f));
}
};
/* -------------------------------------------------------------------------- */
ComputeEKin::~ComputeEKin(){};
/* -------------------------------------------------------------------------- */
template <typename Cont> enable_if_md<Cont> ComputeEKin::build(Cont &cont) {
this->copyContainerInfo(cont);
constexpr UInt Dim = Cont::Dim;
Array1D<Dim> computed_ekin;
computed_ekin = 0.;
UInt cpt = cont.size();
for (auto &&at : cont) {
Array1D<Dim> vel = at.velocity();
computed_ekin += vel * vel * at.mass();
}
computed_ekin *= 0.5 * code_unit_system.mvv2e;
for (UInt i = 0; i < Dim; ++i) {
DUMP("computed ekin[" << i << "] = " << computed_ekin[i], DBG_INFO);
}
CommGroup &group = cont.getCommGroup();
group.allReduce(computed_ekin.data(), Dim, "kinetic energy reduce", OP_SUM);
group.allReduce(&cpt, 1, "kinetic energy counter reduce", OP_SUM);
Real ekin = 0.0;
for (UInt i = 0; i < Dim; ++i) {
DUMP("ekin[" << i << "] = " << computed_ekin[i], DBG_INFO);
ekin += computed_ekin[i];
}
// compute temperature
UnitsConverter unit;
unit.setInUnits(code_unit_system);
unit.setOutUnits(real_unit_system);
unit.computeConversions();
Real T = 2. / cpt / Dim / boltzmann * (ekin * unit.getConversion<Energy>());
// convert the kinetic energies to the required unit systems
// unit.setOutUnits(this->unit_system);
// unit.computeConversions();
// ekin *= unit.getConversion<Energy>();
// for (UInt i = 0; i < Dim; ++i) computed_ekin[i] *=
// unit.getConversion<Energy>();
DUMP("nb_points = " << cpt, DBG_INFO);
this->clear();
Real nb_points = static_cast<Real>(cpt);
this->getArray("nb_points").push_back(nb_points);
this->getArray("Temperature").push_back(T);
- this->getArray().push_back(ekin);
- for (UInt i = 0; i < Dim; ++i)
- this->getArray().push_back(computed_ekin[i]);
+ this->getArray("EKin").push_back(ekin);
+ this->getArray("EKin_vector").push_back(computed_ekin);
// this->name_computed.clear();
// this->name_computed.push_back("nbpoints");
// this->name_computed.push_back("Temp");
// this->name_computed.push_back("EKin");
// this->name_computed.push_back("EKin_x");
// if (Dim > 1)
// this->name_computed.push_back("EKin_y");
// if (Dim == 3)
// this->name_computed.push_back("EKin_z");
}
/* -------------------------------------------------------------------------- */
template <typename Cont> enable_if_mesh<Cont> ComputeEKin::build(Cont &cont) {
constexpr UInt Dim = Cont::Dim;
if (based_on_nodes_flag) {
buildBasedOnNodes<Dim>(cont);
return;
}
Real computed_ekin = 0.0;
for (auto &&el : cont.getContainerElems()) {
computed_ekin += el.getEKin();
}
CommGroup &group = cont.getCommGroup();
group.allReduce(&computed_ekin, 1, " kinetic energy reduce", OP_SUM);
this->clear();
this->getArray().push_back(computed_ekin);
// this->name_computed.clear();
// this->name_computed.push_back("EKin");
}
/* -------------------------------------------------------------------------- */
template <UInt Dim, typename Cont>
void ComputeEKin::buildBasedOnNodes(Cont &cont) {
Real computed_ekin[3] = {0, 0, 0};
UInt cpt = cont.size();
for (auto &&at : cont) {
auto vel = at.velocity();
for (UInt i = 0; i < Dim; ++i)
computed_ekin[i] +=
0.5 * at.mass() * vel[i] * vel[i] * code_unit_system.mvv2e;
}
// Communicator &comm = Communicator::getCommunicator();
CommGroup &group = cont.getCommGroup();
group.allReduce(computed_ekin, Dim, "kinetic energy reduce", OP_SUM);
group.allReduce(&cpt, 1, "kinetic energy counter reduce", OP_SUM);
Real ekin = 0.0;
for (UInt i = 0; i < Dim; ++i) {
DUMP("ekin[" << i << "] = " << computed_ekin[i], DBG_INFO);
ekin += computed_ekin[i];
}
DUMP("nb_points = " << cpt, DBG_INFO);
this->clear();
Real nb_points = static_cast<Real>(cpt);
this->getArray("nb_points").push_back(nb_points);
this->getArray("EKin").push_back(ekin);
for (UInt i = 0; i < Dim; ++i)
this->getArray().push_back(computed_ekin[i]);
// this->name_computed.clear();
// this->name_computed.push_back("nbpoints");
// this->name_computed.push_back("EKin");
// this->name_computed.push_back("EKin_x");
// if (Dim > 1)
// this->name_computed.push_back("EKin_y");
// if (Dim == 3)
// this->name_computed.push_back("EKin_z");
DUMP("Kinetic energy based on nodes could be higher than the one computed "
"from elements.",
DBG_WARNING);
}
/* -------------------------------------------------------------------------- */
/* LMDESC EKIN
This computes the kinetic energy of the region provided as input
*/
/* LMEXAMPLE FILTER ekin EKIN INPUT md */
inline void ComputeEKin::declareParams() {
ComputeInterface::declareParams();
/* LMKEYWORD BASED_ON_NODES
Specify that the selection should be based on the
nodes of the element. Valid for only FE.
*/
this->parseTag("BASED_ON_NODES", based_on_nodes_flag, false);
};
/* -------------------------------------------------------------------------- */
DECLARE_COMPUTE_MAKE_CALL(ComputeEKin);
__END_LIBMULTISCALE__
diff --git a/src/filter/compute_interface.hh b/src/filter/compute_interface.hh
index fb562ed..6f1244e 100644
--- a/src/filter/compute_interface.hh
+++ b/src/filter/compute_interface.hh
@@ -1,122 +1,127 @@
/**
* @file compute_interface.hh
*
* @author Guillaume Anciaux <guillaume.anciaux@epfl.ch>
*
* @date Mon Oct 28 19:23:14 2013
*
* @brief This is the interface of all computes
*
* @section LICENSE
*
* Copyright (©) 2010-2011 EPFL (Ecole Polytechnique Fédérale de Lausanne)
* Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)
*
* LibMultiScale is free software: you can redistribute it and/or modify it
* under the
* terms of the GNU Lesser General Public License as published by the Free
* Software Foundation, either version 3 of the License, or (at your option) any
* later version.
*
* LibMultiScale is distributed in the hope that it will be useful, but
* WITHOUT ANY
* WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
* A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more
* details.
*
* You should have received a copy of the GNU Lesser General Public License
* along with LibMultiScale. If not, see <http://www.gnu.org/licenses/>.
*
*/
#ifndef __LIBMULTISCALE_COMPUTE_INTERFACE_HH__
#define __LIBMULTISCALE_COMPUTE_INTERFACE_HH__
/* -------------------------------------------------------------------------- */
#include "container_array.hh"
#include "filter_interface.hh"
/* -------------------------------------------------------------------------- */
__BEGIN_LIBMULTISCALE__
/* -------------------------------------------------------------------------- */
class ComputeInterface : public FilterInterface {
public:
ComputeInterface();
virtual ~ComputeInterface();
//! gather result routine
ContainerArray<Real> &gatherData(UInt source_rank, UInt root_rank = 0);
//! gather all other processor result routine
ContainerArray<Real> &gatherAllData(UInt root_rank = 0);
//! all gather result routine
ContainerArray<Real> &allGatherAllData();
//! return the total number of elements (parallel safe)
UInt getTotalNbData(UInt root_rank = 0);
inline virtual void compute_make_call() override { LM_TOIMPLEMENT; };
ContainerArray<Real> &getArray(const std::string &output_name = "") {
auto name = output_name;
if (name == "") {
if (this->getOutputs().size() != 1) {
- LM_FATAL("Attached array is not alone: need a name");
+ std::string mesg = "Attached array is not alone: need a name";
+ mesg += "\n\nPossible names: \n\n";
+ for (auto &&pair : this->getOutputs()) {
+ mesg += "\t" + pair.first + "\n";
+ }
+ LM_FATAL(mesg);
}
name = this->getOutputs().begin()->first;
}
return this->getCastedOutput<ContainerArray<Real>>(name, false);
}
CommGroup &getCommGroup() {
if (comm_group == nullptr)
LM_FATAL("invalid group");
return *comm_group;
}
//! set the communication group
virtual void setCommGroup(CommGroup &group) {
comm_group = &group;
for (auto &&pair : this->getOutputs()) {
this->getArray(pair.first).setCommGroup(group);
}
};
//! set the communication group
virtual void copyContainerInfo(ContainerInterface &cont) {
for (auto &&pair : this->getOutputs()) {
this->getArray(pair.first).copyContainerInfo(cont);
}
};
//! set the communication group
virtual void clear() {
for (auto &&pair : this->getOutputs()) {
this->getArray(pair.first).clear();
}
};
protected:
//! vector for gather routine usage
ContainerArray<Real> data_gather;
//! vector for receive counts on allGatherAllData.
std::vector<int> rcnts;
//! flag specifying that input data was gathered(for python)
bool gather_flag;
//! processor number compute_python use to gather data
UInt gather_root_proc;
private:
CommGroup *comm_group;
};
__END_LIBMULTISCALE__
#define DECLARE_COMPUTE(class_name) \
class_name(const std::string &name); \
virtual ~class_name(); \
\
void declareParams() override; \
void compute_make_call() override
#endif /* __LIBMULTISCALE_COMPUTE_INTERFACE_HH__ */
diff --git a/src/filter/compute_time_delta.cc b/src/filter/compute_time_delta.cc
index 9857780..47835bb 100644
--- a/src/filter/compute_time_delta.cc
+++ b/src/filter/compute_time_delta.cc
@@ -1,159 +1,159 @@
/**
* @file compute_time_delta.cc
*
* @author Guillaume Anciaux <guillaume.anciaux@epfl.ch>
*
* @date Tue Dec 10 16:00:49 2013
*
* @brief This compute maintains a time history in order to compute the time
* variation
*
* @section LICENSE
*
* Copyright (©) 2010-2011 EPFL (Ecole Polytechnique Fédérale de Lausanne)
* Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)
*
* LibMultiScale is free software: you can redistribute it and/or modify it
* under the
* terms of the GNU Lesser General Public License as published by the Free
* Software Foundation, either version 3 of the License, or (at your option) any
* later version.
*
* LibMultiScale is distributed in the hope that it will be useful, but
* WITHOUT ANY
* WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
* A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more
* details.
*
* You should have received a copy of the GNU Lesser General Public License
* along with LibMultiScale. If not, see <http://www.gnu.org/licenses/>.
*
*/
/* -------------------------------------------------------------------------- */
#include "compute_time_delta.hh"
#include "factory_multiscale.hh"
#include "lm_common.hh"
#include "math_tools.hh"
/* -------------------------------------------------------------------------- */
__BEGIN_LIBMULTISCALE__
/* -------------------------------------------------------------------------- */
ComputeTimeDelta::ComputeTimeDelta(const std::string &name) : LMObject(name) {
delta_size = 1;
// ActionManager::getManager().addObject(*this);
}
/* -------------------------------------------------------------------------- */
ComputeTimeDelta::~ComputeTimeDelta() {}
/* -------------------------------------------------------------------------- */
void ComputeTimeDelta::action() { LM_TOIMPLEMENT; }
/* -------------------------------------------------------------------------- */
template <typename Cont> void ComputeTimeDelta::account(Cont &cont) {
int k = (current_step + delta_size) % this->frequency;
LM_ASSERT(k >= 0 || k == (int)(delta_size % this->frequency),
"this should not happen");
//const UInt Dim = cont.getDim();
previous_data.resize(cont.size());
this->getArray().resize(cont.size());
if (k == (int)(delta_size % this->frequency)) {
for (UInt i = 0; i < cont.size(); ++i)
this->getArray()[i] = cont[i] - previous_data[i];
}
if (k == 0) {
for (UInt i = 0; i < cont.size(); ++i)
previous_data[i] = cont[i];
}
this->last_accounted_step = current_step;
this->last_accounted_k = k;
}
/* -------------------------------------------------------------------------- */
template <typename Cont> void ComputeTimeDelta::build(Cont &cont) {
#ifndef LM_OPTIMIZED
int k = (current_step + delta_size) % this->frequency;
#endif // LM_OPTIMIZED
LM_ASSERT(k == (int)(delta_size % this->frequency),
"time delta is not tallied to correct timesteps");
LM_ASSERT(last_accounted_k == (delta_size % this->frequency),
"time delta is not tallied to correct timesteps");
LM_ASSERT(last_accounted_step == current_step,
"time delta is not tallied to correct timesteps");
- const UInt Dim = cont.getDim();
+ //const UInt Dim = cont.getDim();
// this->name_computed.clear();
// for (UInt i = 0; i < Dim; ++i) {
// std::stringstream str;
// if (cont.name_computed.size() > i)
// str << "D-" << cont.name_computed[i];
// else
// str << "D-field" << i;
// this->name_computed.push_back(str.str());
// }
}
/* -------------------------------------------------------------------------- */
/* LMDESC TIMEDELTA
This computes a time delta from another input compute
*/
/* LMEXAMPLE
COMPUTE time_delta TIMEDELTA INPUT md FIELD displacement\\
COMPUTE sum TIMEAVG INPUT disp
*/
/* LMHERITANCE action_interface*/
void ComputeTimeDelta::declareParams() {
ActionInterface::declareParams();
/* LMKEYWORD SIZE
The time(multiple of timestep size)
over which to compute the delta. The default is 1
*/
this->parseKeyword("SIZE", delta_size);
}
/* -------------------------------------------------------------------------- */
bool ComputeTimeDelta::shouldMakeAction() {
bool sM = this->doesStageMaskMatch();
if (!sM)
return false;
this->setOneStep();
if (this->end_step >= current_step && this->start_step <= current_step)
if (this->frequency != UINT_MAX) {
int k = (current_step + delta_size) % this->frequency;
if (k == 0 || k == (int)delta_size) {
DUMP("should proceed action "
<< this->getID() << " current_step = " << current_step
<< " start " << this->start_step << " end " << this->end_step,
DBG_INFO);
return true;
}
}
return false;
}
/* -------------------------------------------------------------------------- */
DECLARE_COMPUTE_MAKE_CALL(ComputeTimeDelta);
__END_LIBMULTISCALE__
diff --git a/src/stimulation/stimulation_temperature.cc b/src/stimulation/stimulation_temperature.cc
index 0e0553d..820228e 100644
--- a/src/stimulation/stimulation_temperature.cc
+++ b/src/stimulation/stimulation_temperature.cc
@@ -1,136 +1,136 @@
/**
* @file stimulation_temperature.cc
*
* @author Guillaume Anciaux <guillaume.anciaux@epfl.ch>
* @author Srinivasa Babu Ramisetti <srinivasa.ramisetti@epfl.ch>
*
* @date Wed Jul 09 21:59:47 2014
*
* @brief Set a given temperature to a set of DOFs
*
* @section LICENSE
*
* Copyright (©) 2010-2011 EPFL (Ecole Polytechnique Fédérale de Lausanne)
* Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)
*
* LibMultiScale is free software: you can redistribute it and/or modify it
* under the
* terms of the GNU Lesser General Public License as published by the Free
* Software Foundation, either version 3 of the License, or (at your option) any
* later version.
*
* LibMultiScale is distributed in the hope that it will be useful, but
* WITHOUT ANY
* WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
* A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more
* details.
*
* You should have received a copy of the GNU Lesser General Public License
* along with LibMultiScale. If not, see <http://www.gnu.org/licenses/>.
*
*/
/* -------------------------------------------------------------------------- */
/* -------------------------------------------------------------------------- */
#include "stimulation_temperature.hh"
#include "compute_ekin.hh"
#include "lib_continuum.hh"
#include "lib_dd.hh"
#include "lib_md.hh"
#include "lm_common.hh"
#include <random>
/* -------------------------------------------------------------------------- */
__BEGIN_LIBMULTISCALE__
/* -------------------------------------------------------------------------- */
/* LMDESC TEMPERATURE
This stimulator sets initial velocity to input a temperature.
*/
/* LMEXAMPLE STIMULATION TEMPERATURE TEMP 100 SEED 32 */
/* LMEXAMPLE STIMULATION TEMPERATURE FLUX_FLAG FLUX 1 SEED 32 */
/* LMHERITANCE action_interface */
inline void StimulationTemperature::declareParams() {
StimulationInterface::declareParams();
/* LMKEYWORD TEMP
Set the temperature desired to be set
*/
this->parseKeyword("TEMP", temperature);
/* LMKEYWORD SEED
Set the seed for the random generator
*/
this->parseKeyword("SEED", seed);
/* LMKEYWORD HEAT_RATE
Set the heat rate per node (only for finite elements)
*/
this->parseKeyword("HEAT_RATE", heat_rate_per_node, 0.);
/* LMKEYWORD FLUX_FLAG
Set the flux per node (only for finite elements)
*/
this->parseTag("FLUX_FLAG", flux_flag, false);
}
/* -------------------------------------------------------------------------- */
StimulationTemperature::StimulationTemperature(const std::string &name)
: LMObject(name), flux_flag(false){};
/* -------------------------------------------------------------------------- */
StimulationTemperature::~StimulationTemperature(){};
/* -------------------------------------------------------------------------- */
template <typename Cont>
enable_if_md<Cont> StimulationTemperature::stimulate(Cont &cont) {
constexpr UInt Dim = Cont::Dim;
std::default_random_engine generator(seed);
std::normal_distribution<double> distribution;
for (auto &&at : cont) {
auto vel = at.velocity();
for (UInt i = 0; i < Dim; ++i)
vel[i] = distribution(generator);
}
ComputeEKin compute_ekin("ComputeEkin:" + this->getID());
compute_ekin.build(cont);
- Real T_ = compute_ekin.getArray().get(1);
+ Real T_ = compute_ekin.getArray("Temperature").get(0);
Real alpha = sqrt(temperature / T_);
for (auto &&at : cont) {
DUMP("procede a la stimulation de l'atome a la position " << at, DBG_INFO);
at.velocity() *= alpha;
}
compute_ekin.build(cont);
- T_ = compute_ekin.getArray().get(1);
+ T_ = compute_ekin.getArray("EKin").get(0);
}
/* -------------------------------------------------------------------------- */
template <typename Cont>
enable_if_mesh<Cont> StimulationTemperature::stimulate(Cont &cont) {
if (!flux_flag) {
for (auto &&n : cont) {
n.temperature() = this->temperature;
}
} else {
for (auto &&n : cont) {
n.externalHeatRate() = this->heat_rate_per_node;
}
}
}
/* -------------------------------------------------------------------------- */
DECLARE_STIMULATION_MAKE_CALL(StimulationTemperature);
__END_LIBMULTISCALE__