diff --git a/apps/benchmarks/CMakeLists.txt b/apps/benchmarks/CMakeLists.txt
index 7d5901c16064c4d3b103d5af59b40db3d9aa6403..f9c1d6e7a8feeb0d4bfbb6a1e55ab2946c2b8cee 100644
--- a/apps/benchmarks/CMakeLists.txt
+++ b/apps/benchmarks/CMakeLists.txt
@@ -1,4 +1,5 @@
add_subdirectory( AdaptiveMeshRefinementFluidParticleCoupling )
+add_subdirectory( CNT )
add_subdirectory( ComplexGeometry )
add_subdirectory( DEM )
add_subdirectory( MeshDistance )
diff --git a/apps/benchmarks/CNT/01_cnt_film.cpp b/apps/benchmarks/CNT/01_cnt_film.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..547d6f4c56214843fe29d402bace0d92075cb9fd
--- /dev/null
+++ b/apps/benchmarks/CNT/01_cnt_film.cpp
@@ -0,0 +1,376 @@
+//======================================================================================================================
+//
+// This file is part of waLBerla. waLBerla is free software: you can
+// redistribute it and/or modify it under the terms of the GNU General Public
+// License as published by the Free Software Foundation, either version 3 of
+// the License, or (at your option) any later version.
+//
+// waLBerla 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 General Public License
+// for more details.
+//
+// You should have received a copy of the GNU General Public License along
+// with waLBerla (see COPYING.txt). If not, see <http://www.gnu.org/licenses/>.
+//
+//! \file
+//! \author Igor Ostanin <i.ostanin@skoltech.ru>
+//! \author Grigorii Drozdov <drozd013@umn.edu>
+//! \author Sebastian Eibl <sebastian.eibl@fau.de>
+//
+//======================================================================================================================
+
+#include "Config.h"
+#include "InitializeCNTs.h"
+#include "SphericalSegmentAccessor.h"
+#include "SQLProperties.h"
+#include "Statistics.h"
+#include "TerminalColors.h"
+
+#include "mesa_pd/domain/BlockForestDomain.h"
+#include "mesa_pd/data/Flags.h"
+#include "mesa_pd/data/LinkedCells.h"
+#include "mesa_pd/data/ParticleStorage.h"
+#include "mesa_pd/data/SparseLinkedCells.h"
+#include <mesa_pd/kernel/AssocToBlock.h>
+#include "mesa_pd/kernel/InsertParticleIntoLinkedCells.h"
+#include "mesa_pd/kernel/InsertParticleIntoSparseLinkedCells.h"
+#include "mesa_pd/kernel/ParticleSelector.h"
+#include "mesa_pd/kernel/cnt/AnisotropicVDWContact.h"
+#include "mesa_pd/kernel/cnt/VBondContact.h"
+#include "mesa_pd/kernel/cnt/ViscousDamping.h"
+#include "mesa_pd/kernel/cnt/Parameters.h"
+#include "mesa_pd/kernel/PFCDamping.h"
+#include "mesa_pd/kernel/VelocityVerlet.h"
+#include "mesa_pd/mpi/ContactFilter.h"
+#include "mesa_pd/mpi/notifications/ForceTorqueNotification.h"
+#include "mesa_pd/mpi/ReduceProperty.h"
+#include "mesa_pd/mpi/SyncNextNeighbors.h"
+#include "mesa_pd/sorting/LinearizedCompareFunctor.h"
+#include "mesa_pd/vtk/ParticleVtkOutput.h"
+
+#include "blockforest/Initialization.h"
+#include "core/Environment.h"
+#include "core/math/Constants.h"
+#include "core/timing/TimingPool.h"
+#include "core/timing/TimingTree.h"
+#include "sqlite/SQLite.h"
+#include "vtk/VTKOutput.h"
+
+namespace walberla {
+namespace mesa_pd {
+
+int main(int argc, char **argv)
+{
+ Environment env(argc, argv);
+ walberla::mpi::MPIManager::instance()->useWorldComm();
+
+ logging::Logging::instance()->setStreamLogLevel(logging::Logging::INFO);
+ logging::Logging::instance()->setFileLogLevel(logging::Logging::INFO);
+
+ WcTimingPool globalTP;
+ WcTimingPool loopTP;
+ WcTimingTree tt;
+
+ WALBERLA_LOG_INFO_ON_ROOT("loading configuration parameters");
+ globalTP["load config"].start();
+ //---------------Generation of CNT configuration-----------------------------
+ FilmSpecimen filmSpecimen; // Data structure for CNT film
+ auto cfg = env.config();
+ if (cfg == nullptr)
+ WALBERLA_ABORT("No config specified!");
+ const Config::BlockHandle cntConf = cfg->getBlock("CNT");
+ loadFromConfig(filmSpecimen, cntConf);
+ //===========================================================================
+ globalTP["load config"].end();
+
+ WALBERLA_LOG_INFO_ON_ROOT("creating domain setup");
+ globalTP["domain setup"].start();
+ Vector3<bool> periodicity = Vector3<bool>(true, true, false);
+ if (filmSpecimen.oopp) periodicity[2] = true;
+ real_t maxZ = filmSpecimen.sizeZ;
+ real_t minZ = 0;
+ if (!filmSpecimen.oopp)
+ {
+ maxZ += real_c(filmSpecimen.numSegs) * filmSpecimen.spacing;
+ minZ -= real_c(filmSpecimen.numSegs) * filmSpecimen.spacing;
+ }
+
+ auto forest = blockforest::createBlockForest(
+ AABB(0, 0, minZ, filmSpecimen.sizeX, filmSpecimen.sizeY, maxZ), // simulation domain
+ Vector3<uint_t>(filmSpecimen.numBlocksX, filmSpecimen.numBlocksY, filmSpecimen.numBlocksZ), // blocks in each direction
+ periodicity // periodicity
+ );
+ if (!forest) WALBERLA_ABORT("No BlockForest created ... exiting!");
+
+ auto domain = std::make_shared<domain::BlockForestDomain> (forest);
+
+ auto simulationDomain = forest->getDomain();
+ WALBERLA_UNUSED(simulationDomain);
+ auto localDomain = forest->begin()->getAABB();
+ for (auto& blk : *forest)
+ {
+ localDomain.merge(blk.getAABB());
+ }
+ globalTP["domain setup"].end();
+
+ WALBERLA_LOG_INFO_ON_ROOT("creating initial particle setup");
+ globalTP["generate CNT"].start();
+ auto ps = std::make_shared<data::ParticleStorage>(10);
+ auto ac = SphericalSegmentAccessor(ps);
+ auto lc = data::SparseLinkedCells(localDomain.getExtended(4_r * kernel::cnt::outer_radius), 4_r * kernel::cnt::outer_radius );
+
+ int64_t numParticles = 0;
+ if (filmSpecimen.initialConfigurationFile.empty())
+ {
+ numParticles = generateCNTs(filmSpecimen, ps, *domain);
+ } else
+ {
+ numParticles = loadCNTs(filmSpecimen.initialConfigurationFile, ps);
+ }
+ walberla::mpi::reduceInplace(numParticles, walberla::mpi::SUM);
+ WALBERLA_LOG_INFO_ON_ROOT("#particles created: " << numParticles);
+ globalTP["generate CNT"].end();
+
+ WALBERLA_LOG_INFO_ON_ROOT("setting up VTK output");
+ auto vtkOutput = make_shared<mesa_pd::vtk::ParticleVtkOutput>(ps);
+ vtkOutput->addOutput<data::SelectParticlePosition>("position");
+ auto vtkWriter = walberla::vtk::createVTKOutput_PointData(vtkOutput,
+ "cnt",
+ 1,
+ filmSpecimen.vtkFolder,
+ "simulation_step",
+ true,
+ true,
+ true,
+ filmSpecimen.useMPIIO);
+
+ WALBERLA_LOG_INFO_ON_ROOT("setting up interaction models");
+ kernel::InsertParticleIntoSparseLinkedCells ipilc;
+ kernel::cnt::AnisotropicVDWContact anisotropicVdwContact;
+ kernel::cnt::VBondContact vBondContact;
+ kernel::cnt::ViscousDamping viscous_damping(filmSpecimen.viscousDamping * 1052.0_r,
+ filmSpecimen.viscousDamping * 1052.0_r);
+ kernel::PFCDamping pfcDamping(filmSpecimen.localDamping);
+ kernel::VelocityVerletPreForceUpdate vv_pre(kernel::cnt::dT);
+ kernel::VelocityVerletPostForceUpdate vv_post(kernel::cnt::dT);
+ mpi::ContactFilter contactFilter;
+ mpi::ReduceProperty RP;
+ mpi::SyncNextNeighbors SNN;
+
+ WALBERLA_LOG_INFO_ON_ROOT("warmup");
+ SNN(*ps, *domain);
+ sorting::LinearizedCompareFunctor linear(lc.domain_, Vector3<uint_t>(lc.numCellsPerDim_[0]));
+ ps->sort(linear);
+ vtkWriter->write(true);
+
+ WALBERLA_LOG_INFO_ON_ROOT("running simulation");
+ int64_t contactsChecked = 0;
+ int64_t contactsDetected = 0;
+ int64_t contactsTreated = 0;
+
+ std::vector<Statistics> statistics;
+
+ globalTP["sim loop"].start();
+ for (auto i = 0; i < filmSpecimen.simulationSteps; ++i)
+ {
+ Statistics stats;
+
+ loopTP["VV pre"].start();
+ ps->forEachParticle(false,
+ kernel::SelectAll(),
+ ac,
+ vv_pre,
+ ac);
+ loopTP["VV pre"].end();
+
+ loopTP["SNN"].start();
+ SNN(*ps, *domain);
+ loopTP["SNN"].end();
+
+ loopTP["GenerateSparseLinkedCells"].start();
+ lc.clear();
+ ps->forEachParticle(false, kernel::SelectAll(), ac, ipilc, ac, lc);
+ loopTP["GenerateSparseLinkedCells"].end();
+
+ loopTP["interaction"].start();
+ contactsChecked = 0;
+ contactsDetected = 0;
+ contactsTreated = 0;
+ lc.forEachParticlePairHalf(false,
+ kernel::SelectAll(),
+ ac,
+ [&](size_t p_idx1, size_t p_idx2)
+ {
+ //ensure ordering
+ if (ac.getUid(p_idx1) > ac.getUid(p_idx2))
+ std::swap(ac.getUidRef(p_idx1), ac.getUidRef(p_idx1));
+
+ ++contactsChecked;
+ auto dist = ac.getPosition(p_idx1) - ac.getPosition(p_idx2);
+ auto dist2 = math::dot(dist, dist);
+ auto cutoff = ac.getInteractionRadius(p_idx1) + ac.getInteractionRadius(p_idx2);
+ if (dist2 < cutoff * cutoff)
+ {
+ auto contactPoint = ac.getPosition(p_idx2) + 0.5_r * dist;
+ ++contactsDetected;
+ if (contactFilter(p_idx1, p_idx2, ac, contactPoint, *domain))
+ {
+ ++contactsTreated;
+ constexpr bool period_X = false;
+ constexpr bool period_Y = false;
+ constexpr int max_segsX = 0;
+ constexpr int max_segsY = 0;
+ auto sameCluster = ac.getClusterID(p_idx2) == ac.getClusterID(p_idx1);
+ auto segmentDistance = std::abs(ac.getSegmentID(p_idx2) - ac.getSegmentID(p_idx1));
+ if ((sameCluster) && ((segmentDistance == 1) ||
+ ((period_X) && (segmentDistance == max_segsX - 1)) ||
+ ((period_Y) && (segmentDistance == max_segsY - 1))
+ )
+ )
+ {
+ vBondContact(p_idx1, p_idx2, ac);
+ stats.bendingEnergy += vBondContact.getLastBendingEnergy();
+ stats.shearEnergy += vBondContact.getLastShearEnergy();
+ stats.tensileEnergy += vBondContact.getLastTensileEnergy();
+ stats.twistingEnergy += vBondContact.getLastTwistingEnergy();
+ stats.stretchEnergy += vBondContact.getLastBendingEnergy() +
+ vBondContact.getLastShearEnergy() +
+ vBondContact.getLastTensileEnergy() +
+ vBondContact.getLastTwistingEnergy();
+ }
+ else if ((sameCluster) && (segmentDistance < 20))
+ {
+ //do nothing
+ } else
+ {
+ anisotropicVdwContact(p_idx1, p_idx2, ac);
+ stats.vdwEnergy += anisotropicVdwContact.getLastEnergy();
+ stats.numAlignedPairs += anisotropicVdwContact.isParallel() ? 1 : 0;
+ }
+ viscous_damping(p_idx1, p_idx2, ac);
+ }
+ }
+ });
+ loopTP["interaction"].end();
+
+ loopTP["ReduceForce"].start();
+ RP.operator()<ForceTorqueNotification>(*ps);
+ loopTP["ReduceForce"].end();
+
+ loopTP["PFDDamping"].start();
+ ps->forEachParticle(false, kernel::SelectAll(), ac, pfcDamping, ac);
+ loopTP["PFDDamping"].end();
+
+ loopTP["VV post"].start();
+ ps->forEachParticle(false,
+ kernel::SelectAll(),
+ ac,
+ vv_post,
+ ac);
+ loopTP["VV post"].end();
+
+ if ((i + 1) % filmSpecimen.saveVTKEveryNthStep == 0)
+ {
+ loopTP["save vtk"].start();
+ WALBERLA_LOG_INFO_ON_ROOT(CYAN << " SAVING VTK " << RESET);
+ vtkWriter->write(true);
+ loopTP["save vtk"].end();
+ }
+ if ((i + 1) % filmSpecimen.saveEnergyEveryNthStep == 0)
+ {
+ loopTP["save energy"].start();
+ WALBERLA_LOG_INFO_ON_ROOT(CYAN << " SAVING ENERGY " << RESET);
+ ps->forEachParticle(false,
+ kernel::SelectLocal(),
+ ac,
+ [&](const size_t p_idx)
+ {
+ stats.kineticEnergy += 0.5_r * ac.getMass(p_idx) * ac.getLinearVelocity(p_idx) * ac.getLinearVelocity(p_idx);
+ stats.kineticEnergy += 0.5_r * ac.getAngularVelocity(p_idx) * ac.getInertiaBF(p_idx) * ac.getAngularVelocity(p_idx);
+ });
+ walberla::mpi::reduceInplace(stats, walberla::mpi::SUM);
+ WALBERLA_ROOT_SECTION()
+ {
+ statistics.push_back(stats);
+ saveStatistics(filmSpecimen.energyFolder + "/statistics.txt", statistics);
+ }
+ loopTP["save energy"].end();
+ }
+ if ((i + 1) % filmSpecimen.saveConfEveryNthStep == 0)
+ {
+ loopTP["save config"].start();
+ WALBERLA_LOG_INFO_ON_ROOT(CYAN << " SAVING CONFIG " << RESET);
+ saveConfig(ps,
+ filmSpecimen.confFolder + "/conf_" + std::to_string((i + 1) * 0.00002) + "ns",
+ IO_MODE::USE_MPI_IO);
+ loopTP["save config"].end();
+ }
+ }
+ globalTP["sim loop"].end();
+
+ auto loopTPReduced = loopTP.getReduced();
+ WALBERLA_LOG_DEVEL_ON_ROOT(*loopTPReduced);
+ auto globalTPReduced = globalTP.getReduced();
+ WALBERLA_LOG_DEVEL_ON_ROOT(*globalTPReduced);
+ auto ttReduced = tt.getReduced();
+ WALBERLA_LOG_DEVEL_ON_ROOT(ttReduced);
+
+ numParticles = 0;
+ int64_t numGhostParticles = 0;
+ ps->forEachParticle(false,
+ kernel::SelectAll(),
+ ac,
+ [&](const size_t idx)
+ {
+ if (data::particle_flags::isSet( ac.getFlagsRef(idx), data::particle_flags::GHOST))
+ {
+ ++numGhostParticles;
+ } else
+ {
+ ++numParticles;
+ }
+ });
+
+ walberla::mpi::reduceInplace(numParticles, walberla::mpi::SUM);
+ walberla::mpi::reduceInplace(numGhostParticles, walberla::mpi::SUM);
+ walberla::mpi::reduceInplace(contactsDetected, walberla::mpi::SUM);
+ walberla::mpi::reduceInplace(contactsTreated, walberla::mpi::SUM);
+
+
+ WALBERLA_ROOT_SECTION()
+ {
+ std::map<std::string, walberla::int64_t> integerProperties;
+ std::map<std::string, double> realProperties;
+ std::map<std::string, std::string> stringProperties;
+
+ addBuildInfoToSQL(integerProperties, realProperties, stringProperties);
+ addDomainPropertiesToSQL(*forest, integerProperties, realProperties, stringProperties);
+ addSlurmPropertiesToSQL(integerProperties, realProperties, stringProperties);
+
+ integerProperties["num_particles"] = numParticles;
+ integerProperties["num_ghost_particles"] = numGhostParticles;
+
+ integerProperties["contacts_detected"] = contactsDetected;
+ integerProperties["contacts_treated"] = contactsTreated;
+
+ auto runId = sqlite::storeRunInSqliteDB(filmSpecimen.sqlFile,
+ integerProperties,
+ stringProperties,
+ realProperties);
+ sqlite::storeTimingPoolInSqliteDB(filmSpecimen.sqlFile, runId, *loopTPReduced, "loop");
+ sqlite::storeTimingPoolInSqliteDB(filmSpecimen.sqlFile, runId, *globalTPReduced, "global");
+ sqlite::storeTimingTreeInSqliteDB(filmSpecimen.sqlFile, runId, tt, "tt");
+ }
+
+ return EXIT_SUCCESS;
+}
+
+} // namespace mesa_pd
+} // namespace walberla
+
+int main(int argc, char *argv[])
+{
+ return walberla::mesa_pd::main(argc, argv);
+}
+
diff --git a/apps/benchmarks/CNT/CMakeLists.txt b/apps/benchmarks/CNT/CMakeLists.txt
new file mode 100644
index 0000000000000000000000000000000000000000..41c52ccc8eb335732ecd0159128b017811da6975
--- /dev/null
+++ b/apps/benchmarks/CNT/CMakeLists.txt
@@ -0,0 +1,10 @@
+waLBerla_link_files_to_builddir( *.dat )
+
+waLBerla_add_executable( NAME 01_cnt_film
+ FILES 01_cnt_film.cpp
+ Config.cpp
+ FilmSpecimen.cpp
+ InitializeCNTs.cpp
+ SQLProperties.cpp
+ Statistics.cpp
+ DEPENDS blockforest core mesa_pd sqlite vtk )
diff --git a/apps/benchmarks/CNT/Config.cpp b/apps/benchmarks/CNT/Config.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..ee19d2a97e224628c162af49a4c8f64822305a12
--- /dev/null
+++ b/apps/benchmarks/CNT/Config.cpp
@@ -0,0 +1,108 @@
+//======================================================================================================================
+//
+// This file is part of waLBerla. waLBerla is free software: you can
+// redistribute it and/or modify it under the terms of the GNU General Public
+// License as published by the Free Software Foundation, either version 3 of
+// the License, or (at your option) any later version.
+//
+// waLBerla 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 General Public License
+// for more details.
+//
+// You should have received a copy of the GNU General Public License along
+// with waLBerla (see COPYING.txt). If not, see <http://www.gnu.org/licenses/>.
+//
+//! \file
+//! \author Igor Ostanin <i.ostanin@skoltech.ru>
+//! \author Grigorii Drozdov <drozd013@umn.edu>
+//! \author Sebastian Eibl <sebastian.eibl@fau.de>
+//
+//======================================================================================================================
+
+#include "Config.h"
+#include "TerminalColors.h"
+
+#include "core/mpi/Gatherv.h"
+#include "core/mpi/MPIIO.h"
+#include "core/mpi/RecvBuffer.h"
+#include "core/mpi/SendBuffer.h"
+
+namespace walberla {
+namespace mesa_pd {
+
+void saveConfig(const shared_ptr<data::ParticleStorage>& ps,
+ const std::string& filename,
+ const IO_MODE io_mode)
+{
+ mpi::SendBuffer sb;
+
+ for (const auto& p : *ps)
+ {
+ using namespace walberla::mesa_pd::data::particle_flags;
+ if (isSet(p.getFlags(), GHOST)) continue; //skip ghosts
+
+ sb << static_cast<int64_t> (p.getSegmentID());
+ sb << static_cast<int64_t> (p.getClusterID());
+// sb << static_cast<int64_t> (p.getGroupID());
+
+ sb << static_cast<double> (p.getPosition()[0]);
+ sb << static_cast<double> (p.getPosition()[1]);
+ sb << static_cast<double> (p.getPosition()[2]);
+
+ Vec3 orient = p.getRotation().getMatrix() * Vec3(1,0,0);
+ // Now orient has the shape (sin(th)cos(ph), sin(th)sin(ph), cos(th))
+ sb << static_cast<double> (std::acos(orient[2]));
+ sb << static_cast<double> (std::atan2(orient[1], orient[0]));
+
+ sb << static_cast<double> (p.getLinearVelocity()[0]);
+ sb << static_cast<double> (p.getLinearVelocity()[1]);
+ sb << static_cast<double> (p.getLinearVelocity()[2]);
+
+ sb << static_cast<double> (p.getAngularVelocity()[0]);
+ sb << static_cast<double> (p.getAngularVelocity()[1]);
+ sb << static_cast<double> (p.getAngularVelocity()[2]);
+ }
+
+ switch (io_mode)
+ {
+ case IO_MODE::REDUCE_TO_ROOT:
+ {
+ WALBERLA_LOG_INFO_ON_ROOT(CYAN << "Saving configuration (reduce to root): " << filename << RESET);
+ walberla::mpi::RecvBuffer rb;
+ walberla::mpi::gathervBuffer(sb, rb);
+ WALBERLA_ROOT_SECTION()
+ {
+ uint_t dataSize = sizeof(mpi::RecvBuffer::ElementType) * rb.size();
+ std::ofstream binfile;
+ binfile.open(filename + ".sav",
+ std::ios::out | std::ios::binary);
+ binfile.write(reinterpret_cast< const char * >(rb.ptr()), numeric_cast<std::streamsize>(dataSize));
+ binfile.close();
+ }
+ break;
+ }
+ case IO_MODE::USE_MPI_IO:
+ {
+ WALBERLA_LOG_INFO_ON_ROOT(CYAN << "Saving configuration (MPIIO): " << filename << RESET);
+ walberla::mpi::writeMPIIO(filename + ".mpiio", sb);
+ break;
+ }
+ case IO_MODE::ONE_FILE_PER_RANK:
+ {
+ WALBERLA_LOG_INFO_ON_ROOT(CYAN << "Saving configuration (one file per process): " << filename << RESET);
+ uint_t dataSize = sizeof(mpi::SendBuffer::ElementType) * sb.size();
+ std::ofstream binfile;
+ binfile.open(filename + "_" + std::to_string(mpi::MPIManager::instance()->rank()) + ".sav",
+ std::ios::out | std::ios::binary);
+ binfile.write(reinterpret_cast< const char * >(sb.ptr()), numeric_cast<std::streamsize>(dataSize));
+ binfile.close();
+ }
+ break;
+ default:
+ WALBERLA_ABORT("Unknown IO_MODE");
+ }
+}
+
+} //namespace mesa_pd
+} //namespace walberla
diff --git a/apps/benchmarks/CNT/Config.h b/apps/benchmarks/CNT/Config.h
new file mode 100644
index 0000000000000000000000000000000000000000..d21d2bf0c60973dfdad7c274d21251e497d1bd7d
--- /dev/null
+++ b/apps/benchmarks/CNT/Config.h
@@ -0,0 +1,37 @@
+//======================================================================================================================
+//
+// This file is part of waLBerla. waLBerla is free software: you can
+// redistribute it and/or modify it under the terms of the GNU General Public
+// License as published by the Free Software Foundation, either version 3 of
+// the License, or (at your option) any later version.
+//
+// waLBerla 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 General Public License
+// for more details.
+//
+// You should have received a copy of the GNU General Public License along
+// with waLBerla (see COPYING.txt). If not, see <http://www.gnu.org/licenses/>.
+//
+//! \file
+//! \author Igor Ostanin <i.ostanin@skoltech.ru>
+//! \author Grigorii Drozdov <drozd013@umn.edu>
+//! \author Sebastian Eibl <sebastian.eibl@fau.de>
+//
+//======================================================================================================================
+
+#pragma once
+
+#include "mesa_pd/data/ParticleStorage.h"
+
+namespace walberla {
+namespace mesa_pd {
+
+enum class IO_MODE {REDUCE_TO_ROOT, USE_MPI_IO, ONE_FILE_PER_RANK};
+
+void saveConfig(const shared_ptr<data::ParticleStorage>& ps,
+ const std::string& filename,
+ const IO_MODE io_mode);
+
+} //namespace mesa_pd
+} //namespace walberla
\ No newline at end of file
diff --git a/apps/benchmarks/CNT/FilmSpecimen.cpp b/apps/benchmarks/CNT/FilmSpecimen.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..460360e7f5002ba0fced1e62ab943da911c7393a
--- /dev/null
+++ b/apps/benchmarks/CNT/FilmSpecimen.cpp
@@ -0,0 +1,122 @@
+//======================================================================================================================
+//
+// This file is part of waLBerla. waLBerla is free software: you can
+// redistribute it and/or modify it under the terms of the GNU General Public
+// License as published by the Free Software Foundation, either version 3 of
+// the License, or (at your option) any later version.
+//
+// waLBerla 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 General Public License
+// for more details.
+//
+// You should have received a copy of the GNU General Public License along
+// with waLBerla (see COPYING.txt). If not, see <http://www.gnu.org/licenses/>.
+//
+//! \file
+//! \author Igor Ostanin <i.ostanin@skoltech.ru>
+//! \author Grigorii Drozdov <drozd013@umn.edu>
+//! \author Sebastian Eibl <sebastian.eibl@fau.de>
+//
+//======================================================================================================================
+
+//======================================================================================================================
+//
+// THIS FILE IS GENERATED - PLEASE CHANGE THE TEMPLATE !!!
+//
+//======================================================================================================================
+
+#include "FilmSpecimen.h"
+#include "TerminalColors.h"
+
+#include <core/logging/Logging.h>
+
+namespace walberla {
+namespace mesa_pd {
+
+void loadFromConfig(FilmSpecimen& params, const Config::BlockHandle& cfg)
+{
+ WALBERLA_LOG_INFO_ON_ROOT(GREEN << "Loading CNT film setup: " << RESET)
+ params.sizeX = cfg.getParameter<real_t>("sizeX", 0 );
+ WALBERLA_LOG_INFO_ON_ROOT(YELLOW << "sizeX: " << GREEN << params.sizeX << RESET);
+
+ params.sizeY = cfg.getParameter<real_t>("sizeY", 0 );
+ WALBERLA_LOG_INFO_ON_ROOT(YELLOW << "sizeY: " << GREEN << params.sizeY << RESET);
+
+ params.sizeZ = cfg.getParameter<real_t>("sizeZ", 0 );
+ WALBERLA_LOG_INFO_ON_ROOT(YELLOW << "sizeZ: " << GREEN << params.sizeZ << RESET);
+
+ params.oopp = cfg.getParameter<bool>("oopp", false );
+ WALBERLA_LOG_INFO_ON_ROOT(YELLOW << "oopp: " << GREEN << params.oopp << RESET);
+
+ params.numBlocksX = cfg.getParameter<uint_t>("numBlocksX", 1 );
+ WALBERLA_LOG_INFO_ON_ROOT(YELLOW << "numBlocksX: " << GREEN << params.numBlocksX << RESET);
+
+ params.numBlocksY = cfg.getParameter<uint_t>("numBlocksY", 1 );
+ WALBERLA_LOG_INFO_ON_ROOT(YELLOW << "numBlocksY: " << GREEN << params.numBlocksY << RESET);
+
+ params.numBlocksZ = cfg.getParameter<uint_t>("numBlocksZ", 1 );
+ WALBERLA_LOG_INFO_ON_ROOT(YELLOW << "numBlocksZ: " << GREEN << params.numBlocksZ << RESET);
+
+ params.min_OOP = cfg.getParameter<real_t>("min_OOP", 0 );
+ WALBERLA_LOG_INFO_ON_ROOT(YELLOW << "min_OOP: " << GREEN << params.min_OOP << RESET);
+
+ params.max_OOP = cfg.getParameter<real_t>("max_OOP", 0 );
+ WALBERLA_LOG_INFO_ON_ROOT(YELLOW << "max_OOP: " << GREEN << params.max_OOP << RESET);
+
+ params.numCNTs = cfg.getParameter<int>("numCNTs", 0 );
+ WALBERLA_LOG_INFO_ON_ROOT(YELLOW << "numCNTs: " << GREEN << params.numCNTs << RESET);
+
+ params.numSegs = cfg.getParameter<int>("numSegs", 0 );
+ WALBERLA_LOG_INFO_ON_ROOT(YELLOW << "numSegs: " << GREEN << params.numSegs << RESET);
+
+ params.spacing = cfg.getParameter<real_t>("spacing", 0 );
+ WALBERLA_LOG_INFO_ON_ROOT(YELLOW << "spacing: " << GREEN << params.spacing << RESET);
+
+ params.localDamping = cfg.getParameter<real_t>("localDamping", 0 );
+ WALBERLA_LOG_INFO_ON_ROOT(YELLOW << "localDamping: " << GREEN << params.localDamping << RESET);
+
+ params.viscousDamping = cfg.getParameter<real_t>("viscousDamping", 0 );
+ WALBERLA_LOG_INFO_ON_ROOT(YELLOW << "viscousDamping: " << GREEN << params.viscousDamping << RESET);
+
+ params.seed = cfg.getParameter<uint_t>("seed", 0 );
+ WALBERLA_LOG_INFO_ON_ROOT(YELLOW << "seed: " << GREEN << params.seed << RESET);
+
+ params.vdW = cfg.getParameter<int>("vdW", 0 );
+ WALBERLA_LOG_INFO_ON_ROOT(YELLOW << "vdW: " << GREEN << params.vdW << RESET);
+
+ params.simulationSteps = cfg.getParameter<int>("simulationSteps", 0 );
+ WALBERLA_LOG_INFO_ON_ROOT(YELLOW << "simulationSteps: " << GREEN << params.simulationSteps << RESET);
+
+ params.saveVTKEveryNthStep = cfg.getParameter<int>("saveVTKEveryNthStep", 0 );
+ WALBERLA_LOG_INFO_ON_ROOT(YELLOW << "saveVTKEveryNthStep: " << GREEN << params.saveVTKEveryNthStep << RESET);
+
+ params.saveEnergyEveryNthStep = cfg.getParameter<int>("saveEnergyEveryNthStep", 0 );
+ WALBERLA_LOG_INFO_ON_ROOT(YELLOW << "saveEnergyEveryNthStep: " << GREEN << params.saveEnergyEveryNthStep << RESET);
+
+ params.saveConfEveryNthStep = cfg.getParameter<int>("saveConfEveryNthStep", 0 );
+ WALBERLA_LOG_INFO_ON_ROOT(YELLOW << "saveConfEveryNthStep: " << GREEN << params.saveConfEveryNthStep << RESET);
+
+ params.vtkFolder = cfg.getParameter<std::string>("vtkFolder", "." );
+ WALBERLA_LOG_INFO_ON_ROOT(YELLOW << "vtkFolder: " << GREEN << params.vtkFolder << RESET);
+
+ params.energyFolder = cfg.getParameter<std::string>("energyFolder", "." );
+ WALBERLA_LOG_INFO_ON_ROOT(YELLOW << "energyFolder: " << GREEN << params.energyFolder << RESET);
+
+ params.confFolder = cfg.getParameter<std::string>("confFolder", "." );
+ WALBERLA_LOG_INFO_ON_ROOT(YELLOW << "confFolder: " << GREEN << params.confFolder << RESET);
+
+ params.useMPIIO = cfg.getParameter<bool>("useMPIIO", false );
+ WALBERLA_LOG_INFO_ON_ROOT(YELLOW << "useMPIIO: " << GREEN << params.useMPIIO << RESET);
+
+ params.sqlFile = cfg.getParameter<std::string>("sqlFile", "cnt.sqlite" );
+ WALBERLA_LOG_INFO_ON_ROOT(YELLOW << "sqlFile: " << GREEN << params.sqlFile << RESET);
+
+ params.initialConfigurationFile = cfg.getParameter<std::string>("initialConfigurationFile", "" );
+ WALBERLA_LOG_INFO_ON_ROOT(YELLOW << "initialConfigurationFile: " << GREEN << params.initialConfigurationFile << RESET);
+
+}
+
+} //namespace mesa_pd
+} //namespace walberla
+
diff --git a/apps/benchmarks/CNT/FilmSpecimen.h b/apps/benchmarks/CNT/FilmSpecimen.h
new file mode 100644
index 0000000000000000000000000000000000000000..5aa041f9e182d02e0c5a181998c690d20e078747
--- /dev/null
+++ b/apps/benchmarks/CNT/FilmSpecimen.h
@@ -0,0 +1,72 @@
+//======================================================================================================================
+//
+// This file is part of waLBerla. waLBerla is free software: you can
+// redistribute it and/or modify it under the terms of the GNU General Public
+// License as published by the Free Software Foundation, either version 3 of
+// the License, or (at your option) any later version.
+//
+// waLBerla 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 General Public License
+// for more details.
+//
+// You should have received a copy of the GNU General Public License along
+// with waLBerla (see COPYING.txt). If not, see <http://www.gnu.org/licenses/>.
+//
+//! \file
+//! \author Igor Ostanin <i.ostanin@skoltech.ru>
+//! \author Grigorii Drozdov <drozd013@umn.edu>
+//! \author Sebastian Eibl <sebastian.eibl@fau.de>
+//
+//======================================================================================================================
+
+//======================================================================================================================
+//
+// THIS FILE IS GENERATED - PLEASE CHANGE THE TEMPLATE !!!
+//
+//======================================================================================================================
+
+#pragma once
+
+#include "core/config/Config.h"
+
+#include <string>
+
+namespace walberla {
+namespace mesa_pd {
+
+struct FilmSpecimen
+{
+ real_t sizeX = 0; //Specimen length (x direction)
+ real_t sizeY = 0; //Specimen length (y direction)
+ real_t sizeZ = 0; //Specimen length (z direction)
+ bool oopp = false; //Specimen out-of-plane periodicity (0 - film, 1 - cnt material)
+ uint_t numBlocksX = 1; //Number of blocks in x direction
+ uint_t numBlocksY = 1; //Number of blocks in y direction
+ uint_t numBlocksZ = 1; //Number of blocks in z direction
+ real_t min_OOP = 0; //Out-of-plane angle minimum
+ real_t max_OOP = 0; //Out-of-plane angle maximum
+ int numCNTs = 0; //Number of CNTs
+ int numSegs = 0; //Number of segments in a CNT
+ real_t spacing = 0; //Segment half-spacing
+ real_t localDamping = 0; //Local damping coefficient
+ real_t viscousDamping = 0; //Viscous damping coefficient
+ uint_t seed = 0; //random generator seed
+ int vdW = 0; //type of vdW interaction model
+ int simulationSteps = 0; //Relaxation duration
+ int saveVTKEveryNthStep = 0; //timesteps between saving VTK outputs
+ int saveEnergyEveryNthStep = 0; //timesteps between saving energies
+ int saveConfEveryNthStep = 0; //timesteps between saving confs
+ std::string vtkFolder = "."; //Folder for VTK files
+ std::string energyFolder = "."; //Folder for energy files
+ std::string confFolder = "."; //Folder for conf files
+ bool useMPIIO = false; //Write a single file instead of one file per process
+ std::string sqlFile = "cnt.sqlite"; //database file
+ std::string initialConfigurationFile = ""; //restart from checkpoint
+};
+
+void loadFromConfig(FilmSpecimen& params,
+ const Config::BlockHandle& cfg);
+
+} //namespace mesa_pd
+} //namespace walberla
\ No newline at end of file
diff --git a/apps/benchmarks/CNT/FilmSpecimen.templ.cpp b/apps/benchmarks/CNT/FilmSpecimen.templ.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..1203c455c24045f502c4538145e14b56544f5c0a
--- /dev/null
+++ b/apps/benchmarks/CNT/FilmSpecimen.templ.cpp
@@ -0,0 +1,49 @@
+//======================================================================================================================
+//
+// This file is part of waLBerla. waLBerla is free software: you can
+// redistribute it and/or modify it under the terms of the GNU General Public
+// License as published by the Free Software Foundation, either version 3 of
+// the License, or (at your option) any later version.
+//
+// waLBerla 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 General Public License
+// for more details.
+//
+// You should have received a copy of the GNU General Public License along
+// with waLBerla (see COPYING.txt). If not, see <http://www.gnu.org/licenses/>.
+//
+//! \file
+//! \author Igor Ostanin <i.ostanin@skoltech.ru>
+//! \author Grigorii Drozdov <drozd013@umn.edu>
+//! \author Sebastian Eibl <sebastian.eibl@fau.de>
+//
+//======================================================================================================================
+
+//======================================================================================================================
+//
+// THIS FILE IS GENERATED - PLEASE CHANGE THE TEMPLATE !!!
+//
+//======================================================================================================================
+
+#include "FilmSpecimen.h"
+#include "TerminalColors.h"
+
+#include <core/logging/Logging.h>
+
+namespace walberla {
+namespace mesa_pd {
+
+void loadFromConfig(FilmSpecimen& params, const Config::BlockHandle& cfg)
+{
+ WALBERLA_LOG_INFO_ON_ROOT(GREEN << "Loading CNT film setup: " << RESET)
+ {%- for param in parameters %}
+ params.{{param.name}} = cfg.getParameter<{{param.type}}>("{{param.name}}", {{param.defValue}} );
+ WALBERLA_LOG_INFO_ON_ROOT(YELLOW << "{{param.name}}: " << GREEN << params.{{param.name}} << RESET);
+ {% endfor %}
+}
+
+} //namespace mesa_pd
+} //namespace walberla
+
+
diff --git a/apps/benchmarks/CNT/FilmSpecimen.templ.h b/apps/benchmarks/CNT/FilmSpecimen.templ.h
new file mode 100644
index 0000000000000000000000000000000000000000..a5eec2182f69181c67097ea8bb0d8c79a31633f4
--- /dev/null
+++ b/apps/benchmarks/CNT/FilmSpecimen.templ.h
@@ -0,0 +1,49 @@
+//======================================================================================================================
+//
+// This file is part of waLBerla. waLBerla is free software: you can
+// redistribute it and/or modify it under the terms of the GNU General Public
+// License as published by the Free Software Foundation, either version 3 of
+// the License, or (at your option) any later version.
+//
+// waLBerla 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 General Public License
+// for more details.
+//
+// You should have received a copy of the GNU General Public License along
+// with waLBerla (see COPYING.txt). If not, see <http://www.gnu.org/licenses/>.
+//
+//! \file
+//! \author Igor Ostanin <i.ostanin@skoltech.ru>
+//! \author Grigorii Drozdov <drozd013@umn.edu>
+//! \author Sebastian Eibl <sebastian.eibl@fau.de>
+//
+//======================================================================================================================
+
+//======================================================================================================================
+//
+// THIS FILE IS GENERATED - PLEASE CHANGE THE TEMPLATE !!!
+//
+//======================================================================================================================
+
+#pragma once
+
+#include "core/config/Config.h"
+
+#include <string>
+
+namespace walberla {
+namespace mesa_pd {
+
+struct FilmSpecimen
+{
+ {%- for param in parameters %}
+ {{param.type}} {{param.name}} = {{param.defValue}}; //{{param.comment}}
+ {%- endfor %}
+};
+
+void loadFromConfig(FilmSpecimen& params,
+ const Config::BlockHandle& cfg);
+
+} //namespace mesa_pd
+} //namespace walberla
diff --git a/apps/benchmarks/CNT/FilmSpecimenGenerator.py b/apps/benchmarks/CNT/FilmSpecimenGenerator.py
new file mode 100644
index 0000000000000000000000000000000000000000..4b91cd1f3a7c3b930285e7ca16c3c1d3dd174238
--- /dev/null
+++ b/apps/benchmarks/CNT/FilmSpecimenGenerator.py
@@ -0,0 +1,74 @@
+# -*- coding: utf-8 -*-
+
+from jinja2 import Environment, FileSystemLoader
+import os
+
+
+class Parameter:
+ def __init__(self, name, type, defValue="", comment=""):
+ """Propery of a data strcuture
+
+ Parameters
+ ----------
+ name : str
+ name of the property
+ type : str
+ type of the property
+ defValue : str
+ default value the property should be initialized with
+ comment : str
+ comment appended to the properties in the header file
+ """
+
+ self.name = name
+ self.type = type
+ self.defValue = defValue
+ self.comment = comment
+
+ def __str__(self):
+ return "name: {}, type: {}, defValue: {} \\\\{}".format(self.name, self.type, self.defValue, self.comment)
+
+
+class Config:
+ def __init__(self):
+ self.parameters = []
+
+ def parameter_exists(self, name):
+ for v in self.parameters:
+ if v.name == name:
+ return True
+ return False
+
+ def add_parameter(self, name, type, defValue, comment):
+ if self.parameter_exists(name):
+ print("parameters already added: " + name)
+ else:
+ self.parameters.append(Parameter(name, type, defValue, comment))
+
+ def generate_file(self, template):
+ context = dict()
+ context["parameters"] = self.parameters
+
+ path = ""
+ filename = template.replace(".templ", "")
+ dirname = os.path.dirname(__file__)
+ env = Environment(loader=FileSystemLoader(dirname))
+ print("generating: " + path + filename)
+ with open(path + filename, "wb") as fout:
+ content = env.get_template(template).render(context)
+ fout.write(content.encode('utf8'))
+
+ def generate(self):
+ print("=" * 90)
+ print("Config File:")
+ print("")
+ print("{0: <30}{1: <30}{2: <30}".format("Name", "Type", "Def. Value"))
+ print("=" * 90)
+ for param in self.parameters:
+ print("{0: <30.29}{1: <30.29}{2: <30.29}".format(param.name, param.type, param.defValue))
+ print("=" * 90)
+
+ self.generate_file("FilmSpecimen.templ.h")
+ self.generate_file("FilmSpecimen.templ.cpp")
+
+
diff --git a/apps/benchmarks/CNT/GenerateFilmSpecimen.py b/apps/benchmarks/CNT/GenerateFilmSpecimen.py
new file mode 100755
index 0000000000000000000000000000000000000000..7bc9a3868294ad08a4402bb81caf76fb436f1e00
--- /dev/null
+++ b/apps/benchmarks/CNT/GenerateFilmSpecimen.py
@@ -0,0 +1,36 @@
+#! /usr/bin/env python3
+# -*- coding: utf-8 -*-
+
+from FilmSpecimenGenerator import Config
+
+cfg = Config()
+cfg.add_parameter("sizeX", "real_t", "0", "Specimen length (x direction)")
+cfg.add_parameter("sizeY", "real_t", "0", "Specimen length (y direction)")
+cfg.add_parameter("sizeZ", "real_t", "0", "Specimen length (z direction)")
+cfg.add_parameter("oopp", "bool", "false", "Specimen out-of-plane periodicity (0 - film, 1 - cnt material)")
+cfg.add_parameter("numBlocksX", "uint_t", "1", "Number of blocks in x direction")
+cfg.add_parameter("numBlocksY", "uint_t", "1", "Number of blocks in y direction")
+cfg.add_parameter("numBlocksZ", "uint_t", "1", "Number of blocks in z direction")
+cfg.add_parameter("min_OOP", "real_t", "0", "Out-of-plane angle minimum")
+cfg.add_parameter("max_OOP", "real_t", "0", "Out-of-plane angle maximum")
+cfg.add_parameter("numCNTs", "int", "0", "Number of CNTs")
+cfg.add_parameter("numSegs", "int", "0", "Number of segments in a CNT")
+cfg.add_parameter("spacing", "real_t", "0", "Segment half-spacing")
+cfg.add_parameter("localDamping", "real_t", "0", "Local damping coefficient")
+cfg.add_parameter("viscousDamping", "real_t", "0", "Viscous damping coefficient")
+cfg.add_parameter("seed", "uint_t", "0", "random generator seed")
+cfg.add_parameter("vdW", "int", "0", "type of vdW interaction model")
+cfg.add_parameter("simulationSteps", "int", "0", "Relaxation duration")
+cfg.add_parameter("saveVTKEveryNthStep", "int", "0", "timesteps between saving VTK outputs")
+cfg.add_parameter("saveEnergyEveryNthStep", "int", "0", "timesteps between saving energies")
+cfg.add_parameter("saveConfEveryNthStep", "int", "0", "timesteps between saving confs")
+cfg.add_parameter("vtkFolder", "std::string", '"."', "Folder for VTK files")
+cfg.add_parameter("energyFolder", "std::string", '"."', "Folder for energy files")
+cfg.add_parameter("confFolder", "std::string", '"."', "Folder for conf files")
+cfg.add_parameter("useMPIIO", "bool", "false", "Write a single file instead of one file per process")
+cfg.add_parameter("sqlFile", "std::string", '"cnt.sqlite"', "database file")
+cfg.add_parameter("initialConfigurationFile", "std::string", '""', "restart from checkpoint")
+
+cfg.generate()
+
+
diff --git a/apps/benchmarks/CNT/InitializeCNTs.cpp b/apps/benchmarks/CNT/InitializeCNTs.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..da93903538c7fbb1f0da84ef60c733a2ec2fcf2c
--- /dev/null
+++ b/apps/benchmarks/CNT/InitializeCNTs.cpp
@@ -0,0 +1,174 @@
+//======================================================================================================================
+//
+// This file is part of waLBerla. waLBerla is free software: you can
+// redistribute it and/or modify it under the terms of the GNU General Public
+// License as published by the Free Software Foundation, either version 3 of
+// the License, or (at your option) any later version.
+//
+// waLBerla 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 General Public License
+// for more details.
+//
+// You should have received a copy of the GNU General Public License along
+// with waLBerla (see COPYING.txt). If not, see <http://www.gnu.org/licenses/>.
+//
+//! \file
+//! \author Igor Ostanin <i.ostanin@skoltech.ru>
+//! \author Grigorii Drozdov <drozd013@umn.edu>
+//! \author Sebastian Eibl <sebastian.eibl@fau.de>
+//
+//======================================================================================================================
+
+#include "InitializeCNTs.h"
+#include "TerminalColors.h"
+
+#include "core/math/Constants.h"
+#include "core/math/Random.h"
+#include "core/mpi/MPIManager.h"
+#include "mesa_pd/kernel/cnt/Parameters.h"
+
+#include <cmath>
+
+namespace walberla{
+namespace mesa_pd {
+
+int64_t generateCNTs(const FilmSpecimen& spec,
+ const shared_ptr<data::ParticleStorage>& ps,
+ const domain::IDomain& domain)
+{
+ auto myRank = mpi::MPIManager::instance()->rank();
+
+ auto CNT_length = 2_r * spec.spacing * real_c(spec.numSegs);
+ // Fixed random seed is necessary for coordinated generation on all MPI proc.
+ auto rand0_1 = math::RealRandom<real_t>(static_cast<std::mt19937::result_type>(spec.seed));
+ // Create an assembly of CNTs
+ int64_t numParticles = 0;
+ for (int tube = 0; tube < spec.numCNTs; tube++)
+ {
+ // This definition of theta provides uniform distribution of random orientations on a sphere when spec.OutOfPlane = 1.
+ real_t theta = 0.5_r * math::pi * spec.min_OOP + (spec.max_OOP - spec.min_OOP) * std::acos(1_r - rand0_1());
+ real_t phi = 2_r * math::pi * rand0_1();
+ real_t pos_x = spec.sizeX * rand0_1();
+ real_t pos_y = spec.sizeY * rand0_1();
+ real_t pos_z = spec.sizeZ * rand0_1();
+
+ for (int segment = 0; segment < spec.numSegs; segment++)
+ {
+ // Generate segment position
+ Vec3 pos;
+ auto r = -0.5_r * CNT_length + 2_r * spec.spacing * real_c(segment);
+ pos[0] = pos_x + r * std::sin(theta) * std::cos(phi);
+ pos[1] = pos_y + r * std::sin(theta) * std::sin(phi);
+ pos[2] = pos_z + r * std::cos(theta);
+
+ // Account for periodicity (will not work if CNTs are longer than domain size)
+ if ((pos[0] > spec.sizeX)) pos[0] -= spec.sizeX;
+ if ((pos[0] < 0)) pos[0] += spec.sizeX;
+ if ((pos[1] > spec.sizeY)) pos[1] -= spec.sizeY;
+ if ((pos[1] < 0)) pos[1] += spec.sizeY;
+ if (spec.oopp)
+ {
+ if ((pos[2] > spec.sizeZ)) pos[2] -= spec.sizeZ;
+ if ((pos[2] < 0)) pos[2] += spec.sizeZ;
+ }
+
+ //Only create segment if it is located on the process.
+ if (domain.isContainedInProcessSubdomain(uint_c(myRank), pos))
+ {
+ data::Particle &&sp = *ps->create();
+ sp.setPosition(pos);
+ sp.setOwner(myRank);
+ sp.setInteractionRadius(kernel::cnt::outer_radius);
+ sp.setSegmentID(segment);
+ sp.setClusterID(tube);
+ sp.getRotationRef().rotate( Vec3(0_r, 1_r, 0_r), -0.5_r * math::pi + theta);
+ sp.getRotationRef().rotate( Vec3(0_r, 0_r, 1_r), phi);
+ numParticles++;
+ }
+ }
+ }
+
+ return numParticles;
+}
+
+int64_t loadCNTs(const std::string& filename,
+ const shared_ptr<data::ParticleStorage>& ps)
+{
+ WALBERLA_LOG_INFO_ON_ROOT(GREEN << "Loading configuration (binary format): " << filename);
+ const auto numProcesses = mpi::MPIManager::instance()->numProcesses();
+ const auto rank = mpi::MPIManager::instance()->rank();
+ //---------Generation of WaLBerla model---------------------------------------------------------------------------------------------------------------------------------
+ // Note that walberla primitives are created below only on MPI process that is responsible
+ // for this branch of blockforest
+
+ int64_t numParticles = 0;
+
+ for (auto i = 0; i < numProcesses; ++i)
+ {
+ WALBERLA_MPI_BARRIER();
+ if (i != rank) continue; //bad practice but with the current file format we have to do this
+ std::ifstream binfile;
+ binfile.open(filename.c_str(), std::ios::in | std::ios::binary);
+ int size;
+ binfile.read((char *) &size, sizeof(int));
+ std::cout << RED << "size read form binary file is" << size << RESET << std::endl;
+ for (int id = 0; id < size; id++)
+ {
+ int ID;
+ int sID;
+ int cID;
+ int gID;
+ double x;
+ double y;
+ double z;
+ double theta;
+ double phi;
+ double vx;
+ double vy;
+ double vz;
+ double wx;
+ double wy;
+ double wz;
+ binfile.read((char *) &ID, sizeof(int));
+ binfile.read((char *) &sID, sizeof(int));
+ binfile.read((char *) &cID, sizeof(int));
+ binfile.read((char *) &gID, sizeof(int));
+ binfile.read((char *) &x, sizeof(double));
+ binfile.read((char *) &y, sizeof(double));
+ binfile.read((char *) &z, sizeof(double));
+ binfile.read((char *) &theta, sizeof(double));
+ binfile.read((char *) &phi, sizeof(double));
+ binfile.read((char *) &vx, sizeof(double));
+ binfile.read((char *) &vy, sizeof(double));
+ binfile.read((char *) &vz, sizeof(double));
+ binfile.read((char *) &wx, sizeof(double));
+ binfile.read((char *) &wy, sizeof(double));
+ binfile.read((char *) &wz, sizeof(double));
+ Vec3 pos;
+ pos[0] = real_c(x);
+ pos[1] = real_c(y);
+ pos[2] = real_c(z);
+
+ data::Particle &&sp = *ps->create();
+ sp.setPosition(pos);
+ sp.setOwner(rank);
+ sp.setInteractionRadius(kernel::cnt::outer_radius);
+ sp.setSegmentID(sID);
+ sp.setClusterID(cID);
+ sp.getRotationRef().rotate( Vec3(0_r, 1_r, 0_r), -0.5_r * math::pi + real_c(theta));
+ sp.getRotationRef().rotate( Vec3(0_r, 0_r, 1_r), real_c(phi));
+ sp.setLinearVelocity( Vec3(real_c(vx), real_c(vy), real_c(vz)) );
+ sp.setAngularVelocity( Vec3(real_c(wx), real_c(wy), real_c(wz)) );
+ numParticles++;
+ }
+ }
+
+ mpi::reduceInplace(numParticles, mpi::SUM);
+ WALBERLA_LOG_INFO_ON_ROOT("#particles created: " << numParticles);
+
+ return numParticles;
+}
+
+} //namespace mesa_pd
+} //namespace walberla
diff --git a/apps/benchmarks/CNT/InitializeCNTs.h b/apps/benchmarks/CNT/InitializeCNTs.h
new file mode 100644
index 0000000000000000000000000000000000000000..501753fb04706923dbc552b198e2a3f98d2d15d2
--- /dev/null
+++ b/apps/benchmarks/CNT/InitializeCNTs.h
@@ -0,0 +1,41 @@
+//======================================================================================================================
+//
+// This file is part of waLBerla. waLBerla is free software: you can
+// redistribute it and/or modify it under the terms of the GNU General Public
+// License as published by the Free Software Foundation, either version 3 of
+// the License, or (at your option) any later version.
+//
+// waLBerla 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 General Public License
+// for more details.
+//
+// You should have received a copy of the GNU General Public License along
+// with waLBerla (see COPYING.txt). If not, see <http://www.gnu.org/licenses/>.
+//
+//! \file
+//! \author Igor Ostanin <i.ostanin@skoltech.ru>
+//! \author Grigorii Drozdov <drozd013@umn.edu>
+//! \author Sebastian Eibl <sebastian.eibl@fau.de>
+//
+//======================================================================================================================
+
+#pragma once
+
+#include "FilmSpecimen.h"
+
+#include "mesa_pd/data/ParticleStorage.h"
+#include "mesa_pd/domain/IDomain.h"
+
+namespace walberla{
+namespace mesa_pd {
+
+int64_t generateCNTs(const FilmSpecimen& spec,
+ const shared_ptr<data::ParticleStorage>& ps,
+ const domain::IDomain& domain);
+
+int64_t loadCNTs(const std::string& filename,
+ const shared_ptr<data::ParticleStorage>& ps);
+
+} //namespace mesa_pd
+} //namespace walberla
diff --git a/apps/benchmarks/CNT/SQLProperties.cpp b/apps/benchmarks/CNT/SQLProperties.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..3daa67d5763ad5732cc6462da5f9df76e9304c23
--- /dev/null
+++ b/apps/benchmarks/CNT/SQLProperties.cpp
@@ -0,0 +1,92 @@
+//======================================================================================================================
+//
+// This file is part of waLBerla. waLBerla is free software: you can
+// redistribute it and/or modify it under the terms of the GNU General Public
+// License as published by the Free Software Foundation, either version 3 of
+// the License, or (at your option) any later version.
+//
+// waLBerla 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 General Public License
+// for more details.
+//
+// You should have received a copy of the GNU General Public License along
+// with waLBerla (see COPYING.txt). If not, see <http://www.gnu.org/licenses/>.
+//
+//! \file
+//! \author Igor Ostanin <i.ostanin@skoltech.ru>
+//! \author Grigorii Drozdov <drozd013@umn.edu>
+//! \author Sebastian Eibl <sebastian.eibl@fau.de>
+//
+//======================================================================================================================
+
+#include "SQLProperties.h"
+
+#include <blockforest/BlockForest.h>
+#include <core/DataTypes.h>
+#include <core/waLBerlaBuildInfo.h>
+#include <core/mpi/MPIManager.h>
+
+namespace walberla {
+namespace mesa_pd {
+
+void addBuildInfoToSQL(std::map<std::string, int64_t> & /*integerProperties*/,
+ std::map<std::string, double> & /*realProperties*/,
+ std::map<std::string, std::string> &stringProperties)
+{
+ stringProperties["walberla_sha"] = walberla::core::buildinfo::gitSHA1();
+ stringProperties["build_type"] = walberla::core::buildinfo::buildType();
+ stringProperties["compiler_flags"] = walberla::core::buildinfo::compilerFlags();
+ stringProperties["build_machine"] = walberla::core::buildinfo::buildMachine();
+}
+
+void addDomainPropertiesToSQL(const ::walberla::blockforest::BlockForest &forest,
+ std::map<std::string, int64_t> &integerProperties,
+ std::map<std::string, double> &realProperties,
+ std::map<std::string, std::string> & /*stringProperties*/ )
+{
+ integerProperties["numMPIProcesses"] = ::walberla::mpi::MPIManager::instance()->numProcesses();
+
+ realProperties["domainXMin"] = forest.getDomain().xMin();
+ realProperties["domainXMax"] = forest.getDomain().xMax();
+ realProperties["domainYMin"] = forest.getDomain().yMin();
+ realProperties["domainYMax"] = forest.getDomain().yMax();
+ realProperties["domainZMin"] = forest.getDomain().zMin();
+ realProperties["domainZMax"] = forest.getDomain().zMax();
+
+ integerProperties["xBlocks"] = ::walberla::int_c(forest.getXSize());
+ integerProperties["yBlocks"] = ::walberla::int_c(forest.getYSize());
+ integerProperties["zBlocks"] = ::walberla::int_c(forest.getZSize());
+
+ integerProperties["xPeriodic"] = (forest.isXPeriodic() ? 1 : 0);
+ integerProperties["yPeriodic"] = (forest.isYPeriodic() ? 1 : 0);
+ integerProperties["zPeriodic"] = (forest.isZPeriodic() ? 1 : 0);
+}
+
+std::string envToString(const char *env)
+{
+ return env != nullptr ? std::string(env) : "";
+}
+
+void addSlurmPropertiesToSQL(std::map<std::string, int64_t> & /*integerProperties*/,
+ std::map<std::string, double> & /*realProperties*/,
+ std::map<std::string, std::string> &stringProperties)
+{
+ stringProperties["SLURM_CLUSTER_NAME"] = envToString(std::getenv("SLURM_CLUSTER_NAME"));
+ stringProperties["SLURM_CPUS_ON_NODE"] = envToString(std::getenv("SLURM_CPUS_ON_NODE"));
+ stringProperties["SLURM_CPUS_PER_TASK"] = envToString(std::getenv("SLURM_CPUS_PER_TASK"));
+ stringProperties["SLURM_JOB_ACCOUNT"] = envToString(std::getenv("SLURM_JOB_ACCOUNT"));
+ stringProperties["SLURM_JOB_ID"] = envToString(std::getenv("SLURM_JOB_ID"));
+ stringProperties["SLURM_JOB_CPUS_PER_NODE"] = envToString(std::getenv("SLURM_JOB_CPUS_PER_NODE"));
+ stringProperties["SLURM_JOB_NAME"] = envToString(std::getenv("SLURM_JOB_NAME"));
+ stringProperties["SLURM_JOB_NUM_NODES"] = envToString(std::getenv("SLURM_JOB_NUM_NODES"));
+ stringProperties["SLURM_NTASKS"] = envToString(std::getenv("SLURM_NTASKS"));
+ stringProperties["SLURM_NTASKS_PER_CORE"] = envToString(std::getenv("SLURM_NTASKS_PER_CORE"));
+ stringProperties["SLURM_NTASKS_PER_NODE"] = envToString(std::getenv("SLURM_NTASKS_PER_NODE"));
+ stringProperties["SLURM_NTASKS_PER_SOCKET"] = envToString(std::getenv("SLURM_NTASKS_PER_SOCKET"));
+ stringProperties["SLURM_CPU_BIND_TYPE"] = envToString(std::getenv("SLURM_CPU_BIND_TYPE"));
+}
+
+} //namespace mesa_pd
+} //namespace walberla
+
diff --git a/apps/benchmarks/CNT/SQLProperties.h b/apps/benchmarks/CNT/SQLProperties.h
new file mode 100644
index 0000000000000000000000000000000000000000..de2523b6661aaa606d2ff8832dd67f3ca035b38a
--- /dev/null
+++ b/apps/benchmarks/CNT/SQLProperties.h
@@ -0,0 +1,48 @@
+//======================================================================================================================
+//
+// This file is part of waLBerla. waLBerla is free software: you can
+// redistribute it and/or modify it under the terms of the GNU General Public
+// License as published by the Free Software Foundation, either version 3 of
+// the License, or (at your option) any later version.
+//
+// waLBerla 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 General Public License
+// for more details.
+//
+// You should have received a copy of the GNU General Public License along
+// with waLBerla (see COPYING.txt). If not, see <http://www.gnu.org/licenses/>.
+//
+//! \file
+//! \author Igor Ostanin <i.ostanin@skoltech.ru>
+//! \author Grigorii Drozdov <drozd013@umn.edu>
+//! \author Sebastian Eibl <sebastian.eibl@fau.de>
+//
+//======================================================================================================================
+
+#pragma once
+
+#include <blockforest/BlockForest.h>
+#include <core/DataTypes.h>
+
+#include <map>
+#include <string>
+
+namespace walberla {
+namespace mesa_pd {
+
+void addBuildInfoToSQL(std::map<std::string, int64_t> &integerProperties,
+ std::map<std::string, double> &realProperties,
+ std::map<std::string, std::string> &stringProperties);
+
+void addDomainPropertiesToSQL(const walberla::blockforest::BlockForest &forest,
+ std::map<std::string, int64_t> &integerProperties,
+ std::map<std::string, double> &realProperties,
+ std::map<std::string, std::string> & /*stringProperties*/ );
+
+void addSlurmPropertiesToSQL(std::map<std::string, int64_t> &integerProperties,
+ std::map<std::string, double> &realProperties,
+ std::map<std::string, std::string> &stringProperties);
+
+} //namespace mesa_pd
+} //namespace walberla
diff --git a/apps/benchmarks/CNT/SphericalSegmentAccessor.h b/apps/benchmarks/CNT/SphericalSegmentAccessor.h
new file mode 100644
index 0000000000000000000000000000000000000000..820c4db467b63ee87d47d70cd4959e2ae64def33
--- /dev/null
+++ b/apps/benchmarks/CNT/SphericalSegmentAccessor.h
@@ -0,0 +1,60 @@
+//======================================================================================================================
+//
+// This file is part of waLBerla. waLBerla is free software: you can
+// redistribute it and/or modify it under the terms of the GNU General Public
+// License as published by the Free Software Foundation, either version 3 of
+// the License, or (at your option) any later version.
+//
+// waLBerla 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 General Public License
+// for more details.
+//
+// You should have received a copy of the GNU General Public License along
+// with waLBerla (see COPYING.txt). If not, see <http://www.gnu.org/licenses/>.
+//
+//! \file
+//! \author Igor Ostanin <i.ostanin@skoltech.ru>
+//! \author Grigorii Drozdov <drozd013@umn.edu>
+//! \author Sebastian Eibl <sebastian.eibl@fau.de>
+//
+//======================================================================================================================
+
+#pragma once
+
+#include "mesa_pd/data/ParticleAccessor.h"
+#include "mesa_pd/kernel/cnt/Parameters.h"
+
+namespace walberla {
+namespace mesa_pd {
+
+class SphericalSegmentAccessor : public data::ParticleAccessor
+{
+public:
+ SphericalSegmentAccessor(std::shared_ptr<data::ParticleStorage>& ps)
+ : ParticleAccessor(ps)
+ {}
+
+ constexpr auto getMass(const size_t /*p_idx*/) const {return kernel::cnt::mass_T;}
+ walberla::real_t& getMassRef(const size_t /*p_idx*/) {WALBERLA_ABORT("not available");}
+ void setMass(const size_t /*p_idx*/, walberla::real_t const & /*v*/) {WALBERLA_ABORT("not available");}
+
+ constexpr auto getInvMass(const size_t /*p_idx*/) const {return 1_r / kernel::cnt::mass_T;}
+ walberla::real_t& getInvMassRef(const size_t /*p_idx*/) {WALBERLA_ABORT("not available");}
+ void setInvMass(const size_t /*p_idx*/, walberla::real_t const & /*v*/) {WALBERLA_ABORT("not available");}
+
+ constexpr auto& getInertiaBF(const size_t /*p_idx*/) const {return Ia;}
+ Mat3& getInertiaBFRef(const size_t /*p_idx*/) {WALBERLA_ABORT("not available");}
+ void setInertiaBF(const size_t /*p_idx*/, Mat3 const & /*v*/) {WALBERLA_ABORT("not available");}
+
+ constexpr auto& getInvInertiaBF(const size_t /*p_idx*/) const {return invI;}
+ Mat3& getInvInertiaBFRef(const size_t /*p_idx*/) {WALBERLA_ABORT("not available");}
+ void setInvInertiaBF(const size_t /*p_idx*/, Mat3 const & /*v*/) {WALBERLA_ABORT("not available");}
+private:
+ // - sphere : I = (2/5)*mass*radius^2
+ static constexpr auto Ia = 0.4_r * kernel::cnt::mass_T * kernel::cnt::inner_radius * kernel::cnt::inner_radius;
+ static constexpr auto invI = Mat3(1_r/Ia, 0_r, 0_r, 0_r, 1_r/Ia, 0_r, 0_r, 0_r, 1_r/Ia);
+};
+
+} //namespace mesa_pd
+} //namespace walberla
diff --git a/apps/benchmarks/CNT/Statistics.cpp b/apps/benchmarks/CNT/Statistics.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..febd6f433caefe35d0811f9f0712b464d911e1ef
--- /dev/null
+++ b/apps/benchmarks/CNT/Statistics.cpp
@@ -0,0 +1,70 @@
+//======================================================================================================================
+//
+// This file is part of waLBerla. waLBerla is free software: you can
+// redistribute it and/or modify it under the terms of the GNU General Public
+// License as published by the Free Software Foundation, either version 3 of
+// the License, or (at your option) any later version.
+//
+// waLBerla 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 General Public License
+// for more details.
+//
+// You should have received a copy of the GNU General Public License along
+// with waLBerla (see COPYING.txt). If not, see <http://www.gnu.org/licenses/>.
+//
+//! \file
+//! \author Igor Ostanin <i.ostanin@skoltech.ru>
+//! \author Grigorii Drozdov <drozd013@umn.edu>
+//! \author Sebastian Eibl <sebastian.eibl@fau.de>
+//
+//======================================================================================================================
+
+#include "Statistics.h"
+#include "TerminalColors.h"
+
+#include "core/logging/Logging.h"
+
+#include <fstream>
+
+namespace walberla {
+namespace mesa_pd {
+
+void saveStatistics(const std::string& filename,
+ const std::vector<Statistics>& stats)
+{
+ WALBERLA_LOG_INFO_ON_ROOT(CYAN << "saving statistics: " << filename << RESET);
+ std::ofstream fout(filename);
+ fout.precision(10);
+ fout << "# kineticEnergy vdwEnergy stretchEnergy tensileEnergy shearEnergy bendingEnergy twistingEnergy numAlignedPairs" << std::endl;
+ for (const auto& stat : stats)
+ {
+ fout << stat.kineticEnergy << " "
+ << stat.vdwEnergy << " "
+ << stat.stretchEnergy << " "
+ << stat.tensileEnergy << " "
+ << stat.shearEnergy << " "
+ << stat.bendingEnergy << " "
+ << stat.twistingEnergy << " "
+ << stat.numAlignedPairs << std::endl;
+ }
+ fout.close();
+}
+
+} //namespace mesa_pd
+
+namespace mpi {
+void reduceInplace( mesa_pd::Statistics& stats, mpi::Operation operation, int recvRank, MPI_Comm comm )
+{
+ reduceInplace(stats.kineticEnergy, operation, recvRank, comm);
+ reduceInplace(stats.vdwEnergy, operation, recvRank, comm);
+ reduceInplace(stats.stretchEnergy, operation, recvRank, comm);
+ reduceInplace(stats.tensileEnergy, operation, recvRank, comm);
+ reduceInplace(stats.shearEnergy, operation, recvRank, comm);
+ reduceInplace(stats.bendingEnergy, operation, recvRank, comm);
+ reduceInplace(stats.twistingEnergy, operation, recvRank, comm);
+ reduceInplace(stats.numAlignedPairs, operation, recvRank, comm);
+}
+} //namespace mpi
+
+} //namespace walberla
\ No newline at end of file
diff --git a/apps/benchmarks/CNT/Statistics.h b/apps/benchmarks/CNT/Statistics.h
new file mode 100644
index 0000000000000000000000000000000000000000..51f614cc59423215dd436181e13419984dfc1f6b
--- /dev/null
+++ b/apps/benchmarks/CNT/Statistics.h
@@ -0,0 +1,54 @@
+//======================================================================================================================
+//
+// This file is part of waLBerla. waLBerla is free software: you can
+// redistribute it and/or modify it under the terms of the GNU General Public
+// License as published by the Free Software Foundation, either version 3 of
+// the License, or (at your option) any later version.
+//
+// waLBerla 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 General Public License
+// for more details.
+//
+// You should have received a copy of the GNU General Public License along
+// with waLBerla (see COPYING.txt). If not, see <http://www.gnu.org/licenses/>.
+//
+//! \file
+//! \author Igor Ostanin <i.ostanin@skoltech.ru>
+//! \author Grigorii Drozdov <drozd013@umn.edu>
+//! \author Sebastian Eibl <sebastian.eibl@fau.de>
+//
+//======================================================================================================================
+
+#pragma once
+
+#include "core/mpi/Reduce.h"
+
+#include <string>
+#include <vector>
+
+namespace walberla {
+namespace mesa_pd {
+
+struct Statistics
+{
+ real_t kineticEnergy = 0_r;
+ real_t vdwEnergy = 0_r;
+ real_t stretchEnergy = 0_r;
+ real_t tensileEnergy = 0_r;
+ real_t shearEnergy = 0_r;
+ real_t bendingEnergy = 0_r;
+ real_t twistingEnergy = 0_r;
+ int64_t numAlignedPairs = 0;
+};
+
+void saveStatistics(const std::string& filename,
+ const std::vector<Statistics>& stats);
+
+} //namespace mesa_pd
+
+namespace mpi {
+void reduceInplace( mesa_pd::Statistics& stats, mpi::Operation operation, int recvRank = 0, MPI_Comm comm = MPI_COMM_WORLD );
+} //namespace mpi
+
+} //namespace walberla
diff --git a/apps/benchmarks/CNT/TerminalColors.h b/apps/benchmarks/CNT/TerminalColors.h
new file mode 100644
index 0000000000000000000000000000000000000000..e572dffec5e357eda894ff9e93d8ba43989cbc3e
--- /dev/null
+++ b/apps/benchmarks/CNT/TerminalColors.h
@@ -0,0 +1,49 @@
+//======================================================================================================================
+//
+// This file is part of waLBerla. waLBerla is free software: you can
+// redistribute it and/or modify it under the terms of the GNU General Public
+// License as published by the Free Software Foundation, either version 3 of
+// the License, or (at your option) any later version.
+//
+// waLBerla 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 General Public License
+// for more details.
+//
+// You should have received a copy of the GNU General Public License along
+// with waLBerla (see COPYING.txt). If not, see <http://www.gnu.org/licenses/>.
+//
+//! \file
+//! \author Igor Ostanin <i.ostanin@skoltech.ru>
+//! \author Grigorii Drozdov <drozd013@umn.edu>
+//! \author Sebastian Eibl <sebastian.eibl@fau.de>
+//
+//======================================================================================================================
+
+#pragma once
+
+#include <string>
+
+namespace walberla {
+
+#ifdef PLAIN_OUTPUT
+//----Terminal color manipulators-----
+const std::string RED(""); // errors and problems
+const std::string GREEN(""); // loading operations
+const std::string YELLOW(""); // information (non-Walberla)
+const std::string PURPLE(""); // stages of program implementation
+const std::string CYAN(""); // save/write operations
+const std::string RESET("");
+//====================================
+#else
+//----Terminal color manipulators-----
+const std::string RED("\033[1;31m"); // errors and problems
+const std::string GREEN("\033[1;32m"); // loading operations
+const std::string YELLOW("\033[1;33m"); // information (non-Walberla)
+const std::string PURPLE("\033[1;35m"); // stages of program implementation
+const std::string CYAN("\033[1;36m"); // save/write operations
+const std::string RESET("\033[0m");
+//====================================
+#endif
+
+} //namespace walberla
\ No newline at end of file
diff --git a/apps/benchmarks/CNT/specimen_setup.dat b/apps/benchmarks/CNT/specimen_setup.dat
new file mode 100644
index 0000000000000000000000000000000000000000..19703ed28bf6a12cb29f40f46630bc3869c0ab8c
--- /dev/null
+++ b/apps/benchmarks/CNT/specimen_setup.dat
@@ -0,0 +1,30 @@
+//--------------CNT material generator input file (do not change anything except numerical values!)-
+// This file sets up generation of a rectangular specimen of a CNT material
+// At the generation stage, triple periodic boundary conditions are imposed.
+// Out-of-plane (xy) angle span tunes transverse anizotropy of the film (1 means isotropy)
+// Note that the numerical values should be single space separated
+
+CNT
+{
+ sizeX 20000; // A - Specimen length (x direction).
+ sizeY 20000; // A - Specimen width (y direction).
+ sizeZ 200; // A - Specimen height (z direction).
+ oopp 0; // - out of plane periodicity
+ numBlocksX 2; // - Number of blocks in x direction (blocks are distributed over available MPI processes,
+ numBlocksY 2; // - Number of blocks in y direction nx x ny x nz blocks total
+ numBlocksZ 1; // - Number of blocks in z direction (nx, ny, nz can not be less than 2 )
+ min_OOP 0.995; // (in PI/2 rads) - minimum theta angle (set 0 for isotropy)
+ max_OOP 1.0; // (in PI/2 rads) - maximum theta angle (set 1 for isotropy)
+ numCNTs 1000; // - Number of CNTs.
+ numSegs 800; // - Number of segments in a CNT.
+ spacing 6.78; // - Segment spacing
+ localDamping 0.2; // - "local" damping
+ viscousDamping 0.00; // - viscous damping
+ seed 1201; // - random seed - needs to be changed for a different realization
+ vdW 1; // - vdW model (0 - isotr.vdW (comp speed-1.0), 1 - aniz. vdW (comp speed 0.54), 2 - numer.integrated vdW (comp-speed 0.14))
+ simulationSteps 1000000; // - timesteps relaxation duration
+ saveEveryNthStep 10000; // - timesteps between checkpoints
+ saveConfEveryNthStep 100000; // - timesteps between saving confs
+ file serie_3_spec_1; // - filename to save the resulting conf (saved to ../confs/)
+}
+//==================================================================================================
diff --git a/src/mesa_pd/kernel/cnt/AnisotropicVDWContact.h b/src/mesa_pd/kernel/cnt/AnisotropicVDWContact.h
index df369216e4d5e808ab9c5567897ae13ac7c35f85..5beb8622bed219a931b8264fdf8c750fe1dc5b65 100644
--- a/src/mesa_pd/kernel/cnt/AnisotropicVDWContact.h
+++ b/src/mesa_pd/kernel/cnt/AnisotropicVDWContact.h
@@ -81,6 +81,8 @@ void AnisotropicVDWContact::operator()(const size_t p_idx1,
const size_t p_idx2,
Accessor &ac)
{
+ isParallel_ = false;
+ energy_ = 0_r;
//===Adaptation of PFC5 vdW contact model implementation====
// Getting the orientations of segments
@@ -128,7 +130,6 @@ void AnisotropicVDWContact::operator()(const size_t p_idx1,
// check that cosine belongs [-1,1]
cos_gamma = std::min(1.0_r, cos_gamma);
cos_gamma = std::max(-1.0_r, cos_gamma);
- isParallel_ = false;
if (L < 20_r && L > 16_r)
{
const auto gamma = acos(cos_gamma);
@@ -373,7 +374,6 @@ void AnisotropicVDWContact::operator()(const size_t p_idx1,
} else if (L <= 2_r * R_ * 1.2_r * TH)
{ // Small distance
//WALBERLA_LOG_DEVEL( "Small distance");
- energy_ = 0_r;
real_t F = -1_r;
addForceAtomic(p_idx1, ac, F * n);
addForceAtomic(p_idx2, ac, -F * n);
diff --git a/src/mesa_pd/kernel/cnt/Parameters.h b/src/mesa_pd/kernel/cnt/Parameters.h
index c265088bab02da9837599fbf333032b795bf2b0d..fc4b98f7f666a6fc6707ba13fbeb12bcf2c91ef5 100644
--- a/src/mesa_pd/kernel/cnt/Parameters.h
+++ b/src/mesa_pd/kernel/cnt/Parameters.h
@@ -80,10 +80,10 @@ constexpr auto M_C = 12.011_r;
constexpr auto mass_T = ro * T * M_C * 104.397_r;
/// Volume of a capsule
-double vol_capsule = (4_r/3_r) * math::pi * (R_CNT * R_CNT * R_CNT) + math::pi * R_CNT * R_CNT * T;
+const double vol_capsule = (4_r/3_r) * math::pi * (R_CNT * R_CNT * R_CNT) + math::pi * R_CNT * R_CNT * T;
/// Density of a capsule
-double dens_capsule = mass_T / vol_capsule;
+const double dens_capsule = mass_T / vol_capsule;
/// V-bond parameter
constexpr auto knorm = (En * 0.006242_r) / T;
diff --git a/src/mesa_pd/kernel/cnt/VBondContact.h b/src/mesa_pd/kernel/cnt/VBondContact.h
index 6bc5edf98e7bdb8c3adee076a712ea42bc15dccf..dc252756c77631bf3efc9b1a3cc6eaf6d5bef301 100644
--- a/src/mesa_pd/kernel/cnt/VBondContact.h
+++ b/src/mesa_pd/kernel/cnt/VBondContact.h
@@ -106,6 +106,11 @@ void VBondContact::operator()(const size_t p_idx1,
const size_t p_idx2,
Accessor &ac)
{
+ tensileEnergy = 0_r;
+ shearEnergy = 0_r;
+ bendingEnergy = 0_r;
+ twistingEnergy = 0_r;
+
Vec3 ri = ac.getPosition(p_idx1);
Vec3 rj = ac.getPosition(p_idx2);