From 9304e403b1f30e2118af82dda487170aeeeb8f08 Mon Sep 17 00:00:00 2001 From: Paul Kieckhefen Date: Fri, 16 Feb 2018 18:00:26 +0100 Subject: [PATCH 01/29] Initial commit. Works, but validation pending. --- src/lagrangian/cfdemParticle/Make/files | 2 + src/lagrangian/cfdemParticle/Make/options | 4 +- .../dataExchangeModel/dataExchangeModel.H | 8 + .../dataExchangeModel/twoWayOne2One/one2one.C | 240 +++++ .../dataExchangeModel/twoWayOne2One/one2one.H | 70 ++ .../twoWayOne2One/twoWayOne2One.C | 847 ++++++++++++++++++ .../twoWayOne2One/twoWayOne2One.H | 226 +++++ .../engineSearchMany2Many.C | 11 +- .../engineSearchMany2Many.H | 7 + 9 files changed, 1411 insertions(+), 4 deletions(-) create mode 100644 src/lagrangian/cfdemParticle/subModels/dataExchangeModel/twoWayOne2One/one2one.C create mode 100644 src/lagrangian/cfdemParticle/subModels/dataExchangeModel/twoWayOne2One/one2one.H create mode 100644 src/lagrangian/cfdemParticle/subModels/dataExchangeModel/twoWayOne2One/twoWayOne2One.C create mode 100644 src/lagrangian/cfdemParticle/subModels/dataExchangeModel/twoWayOne2One/twoWayOne2One.H diff --git a/src/lagrangian/cfdemParticle/Make/files b/src/lagrangian/cfdemParticle/Make/files index 7cc1c491..f68eab10 100644 --- a/src/lagrangian/cfdemParticle/Make/files +++ b/src/lagrangian/cfdemParticle/Make/files @@ -141,6 +141,8 @@ $(dataExchangeModels)/oneWayVTK/oneWayVTK.C $(dataExchangeModels)/twoWayFiles/twoWayFiles.C $(dataExchangeModels)/noDataExchange/noDataExchange.C $(dataExchangeModels)/twoWayMPI/twoWayMPI.C +$(dataExchangeModels)/twoWayOne2One/twoWayOne2One.C +$(dataExchangeModels)/twoWayOne2One/one2one.C $(averagingModels)/averagingModel/averagingModel.C $(averagingModels)/averagingModel/newAveragingModel.C diff --git a/src/lagrangian/cfdemParticle/Make/options b/src/lagrangian/cfdemParticle/Make/options index 9c152540..8a84bb35 100644 --- a/src/lagrangian/cfdemParticle/Make/options +++ b/src/lagrangian/cfdemParticle/Make/options @@ -34,6 +34,4 @@ LIB_LIBS = \ -lmpi_cxx \ -Wl,-rpath,$(CFDEM_LIGGGHTS_BIN_DIR) \ -L$(CFDEM_LIGGGHTS_BIN_DIR) \ - -lliggghts \ - -L$(CFDEM_Many2ManyLIB_PATH) \ - -lcoupleMany2Many + -lliggghts diff --git a/src/lagrangian/cfdemParticle/subModels/dataExchangeModel/dataExchangeModel/dataExchangeModel.H b/src/lagrangian/cfdemParticle/subModels/dataExchangeModel/dataExchangeModel/dataExchangeModel.H index e864bfee..bedfb6f1 100755 --- a/src/lagrangian/cfdemParticle/subModels/dataExchangeModel/dataExchangeModel/dataExchangeModel.H +++ b/src/lagrangian/cfdemParticle/subModels/dataExchangeModel/dataExchangeModel/dataExchangeModel.H @@ -251,11 +251,19 @@ public: for (int j=0;j<3;j++) particleCloud_.positions_[i][j]=pos[i*3+j]; } + inline void setCellIDs(label n,int* ID) const { for (int i=0;i. + +Copyright 2018- Paul Kieckhefen, TUHH + +\*---------------------------------------------------------------------------*/ + +//#define O2ODEBUG + +#ifdef O2ODEBUG +#include +#endif + +#include +#include "one2one.H" + +template +struct mpi_type_wrapper { + MPI_Datatype mpi_type; + mpi_type_wrapper(); +}; +template <> mpi_type_wrapper::mpi_type_wrapper() +: mpi_type(MPI_FLOAT) {} +template <> mpi_type_wrapper::mpi_type_wrapper() +: mpi_type(MPI_DOUBLE) {} +template <> mpi_type_wrapper::mpi_type_wrapper() +: mpi_type(MPI_INT) {} + +/* ---------------------------------------------------------------------- */ + +One2One::One2One(MPI_Comm caller) +: +ncollected_(-1), +comm_(caller), +nsrc_procs_(-1), +src_procs_(nullptr), +ndst_procs_(-1), +dst_procs_(nullptr), +nlocal_(-1), +natoms_(nullptr), +request_(nullptr), +status_(nullptr) +{ + MPI_Comm_rank(comm_,&me_); + MPI_Comm_size(comm_,&nprocs_); +} + +/* ---------------------------------------------------------------------- */ + +One2One::~One2One() +{ + deallocate(); +} + +/* ---------------------------------------------------------------------- + communicate particle ids based on processor communication pattern +------------------------------------------------------------------------- */ + +void One2One::setup +( + int nsrc_procs, + int *src_procs, + int ndst_procs, + int *dst_procs, + int nlocal +) +{ + // free any previous one2one info + deallocate(); + + src_procs_ = src_procs; + nsrc_procs_ = nsrc_procs; + dst_procs_ = dst_procs; + ndst_procs_ = ndst_procs; + nlocal_ = nlocal; + + // gather number of ids for reserving memory + natoms_ = new int[nprocs_]; + MPI_Allgather // may be replaced by send/irecv + ( + &nlocal_, + 1, + MPI_INT, + natoms_, + 1, + MPI_INT, + comm_ + ); + + ncollected_ = 0; + for (int i = 0; i < nsrc_procs_; i++) + ncollected_ += natoms_[src_procs_[i]]; + + request_ = new MPI_Request[nsrc_procs_]; + status_ = new MPI_Status[nsrc_procs_]; +} + +/* ---------------------------------------------------------------------- + src: what is present on this proc + dst: what is received from other procs + all comm according to map set up in setup(...) +------------------------------------------------------------------------- */ +template +void One2One::exchange(T *&src, T *&dst, int data_length) +{ + mpi_type_wrapper wrap; + + // post receives + int offset = 0; + for (int i = 0; i < nsrc_procs_; i++) + { + #ifdef O2ODEBUG + std::cout<< "[" << me_ << "]" + << " RCV " << i + << " of " << nsrc_procs_ + << " from: " << src_procs_[i] + << " natoms_[src_procs_[i]] " << natoms_[src_procs_[i]] + << " datalength " << data_length + << " offset " << offset + << std::endl; + #endif + MPI_Irecv + ( + dst+offset, + natoms_[src_procs_[i]]*data_length, + wrap.mpi_type, + src_procs_[i], + MPI_ANY_TAG, + comm_, + &request_[i] + ); + offset += natoms_[src_procs_[i]]*data_length; + } + + // make sure all receives are posted + MPI_Barrier(comm_); + + // blocking sends - do nonblocking instead + // since doing many-2-many here? + for (int i = 0; i < ndst_procs_; i++) + { + #ifdef O2ODEBUG + std::cout<< "[" << me_ << "]" + << " SEND to: " << dst_procs_[i] + << " nlocal_ " << nlocal_ + << " data_length " << data_length + << std::endl; + #endif + MPI_Send + ( + src, + nlocal_*data_length, + wrap.mpi_type, + dst_procs_[i], + 0, + comm_ + ); + } + + MPI_Waitall(nsrc_procs_, request_, status_); + + // copy on-proc info instead of communicating everything + // may yield a first step towards true one2one communication + //for (i = 0; i < nown; i++) + // dest[dest_own[i]] = src[src_own[i]]; + +} + +template void One2One::exchange(int*&, int*&, int); +template void One2One::exchange(double*&, double*&, int); + +// there should be a way to do this without copying data +template +void One2One::exchange(T **&src, T **&dst, int data_length) +{ + mpi_type_wrapper wrap; + + T* tmp_dst = new T[ncollected_*data_length]; + T* tmp_src = new T[nlocal_*data_length]; + + for (int i = 0; i < nlocal_; i++) + for (int j = 0; j < data_length; j++) + tmp_src[data_length*i+j] = src[i][j]; + + exchange(tmp_src, tmp_dst, data_length); + + for (int i = 0; i < ncollected_; i++) + for (int j = 0; j < data_length; j++) + dst[i][j] = tmp_dst[data_length*i+j]; + + delete [] tmp_src; + delete [] tmp_dst; +} +template void One2One::exchange(int**&, int**&, int); +template void One2One::exchange(double**&, double**&, int); + + +template +void One2One::exchange(T **&src, T *&dst, int data_length) +{ + mpi_type_wrapper wrap; + + T* tmp_src = new T[nlocal_*data_length]; + + for (int i = 0; i < nlocal_; i++) + for (int j = 0; j < data_length; j++) + tmp_src[data_length*i+j] = src[i][j]; + + exchange(tmp_src, dst, data_length); + + + delete [] tmp_src; +} +template void One2One::exchange(int**&, int*&, int); +template void One2One::exchange(double**&, double*&, int); + +/* ---------------------------------------------------------------------- */ + +void One2One::deallocate() +{ + delete [] src_procs_; + delete [] dst_procs_; + delete [] natoms_; + + delete [] request_; + delete [] status_; +} \ No newline at end of file diff --git a/src/lagrangian/cfdemParticle/subModels/dataExchangeModel/twoWayOne2One/one2one.H b/src/lagrangian/cfdemParticle/subModels/dataExchangeModel/twoWayOne2One/one2one.H new file mode 100644 index 00000000..4205e11c --- /dev/null +++ b/src/lagrangian/cfdemParticle/subModels/dataExchangeModel/twoWayOne2One/one2one.H @@ -0,0 +1,70 @@ +/*---------------------------------------------------------------------------*\ +License + + This 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. + + This code 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 this code. If not, see . + +Copyright 2018- Paul Kieckhefen, TUHH + +\*---------------------------------------------------------------------------*/ + +#ifndef ONE2ONE_H +#define ONE2ONE_H + +#include + +class One2One +{ + public: + One2One(MPI_Comm); + ~One2One(); + + void setup + ( + int nsrc_procs, + int *src_procs, + int ndst_procs, + int* dst_procs, + int nlocal + ); + + template + void exchange(T *&, T *&, int data_length=1); + template + void exchange(T **&, T **&, int data_length=1); + template + void exchange(T **&, T *&, int data_length=1); + + int ncollected_; // # of ids in from group + + protected: + + MPI_Comm comm_; + int me_, nprocs_; // rank and size + + // communication partners + int nsrc_procs_; // # of off-processor IDs + int* src_procs_; // procs I receive data from + int ndst_procs_; // # of off-processor IDs + int* dst_procs_; // procs I receive data from + + int nlocal_; // # particle ids I own + int* natoms_; + + MPI_Request* request_; + MPI_Status* status_; + + void deallocate(); +}; + +#endif diff --git a/src/lagrangian/cfdemParticle/subModels/dataExchangeModel/twoWayOne2One/twoWayOne2One.C b/src/lagrangian/cfdemParticle/subModels/dataExchangeModel/twoWayOne2One/twoWayOne2One.C new file mode 100644 index 00000000..668dc9e2 --- /dev/null +++ b/src/lagrangian/cfdemParticle/subModels/dataExchangeModel/twoWayOne2One/twoWayOne2One.C @@ -0,0 +1,847 @@ +/*---------------------------------------------------------------------------*\ + CFDEMcoupling - Open Source CFD-DEM coupling + + CFDEMcoupling is part of the CFDEMproject + www.cfdem.com + Christoph Goniva, christoph.goniva@cfdem.com + Copyright 2009-2012 JKU Linz + Copyright 2012- DCS Computing GmbH, Linz +------------------------------------------------------------------------------- +License + This file is part of CFDEMcoupling. + + CFDEMcoupling 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. + + CFDEMcoupling 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 CFDEMcoupling; if not, write to the Free Software Foundation, + Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA + +Description + This code is designed to realize coupled CFD-DEM simulations using LIGGGHTS + and OpenFOAM(R). Note: this code is not part of OpenFOAM(R) (see DISCLAIMER). + +Contributing authors + Paul Kieckhefen (TUHH) 2018- + +\*---------------------------------------------------------------------------*/ + + +#include "twoWayOne2One.H" +#include "addToRunTimeSelectionTable.H" +#include "clockModel.H" +#include "pair.h" +#include "force.h" +#include "forceModel.H" + +#include +#include + + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +namespace Foam +{ + +// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * // + +defineTypeNameAndDebug(twoWayOne2One, 0); + +addToRunTimeSelectionTable +( + dataExchangeModel, + twoWayOne2One, + dictionary +); + + +// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * // + +// Construct from components +twoWayOne2One::twoWayOne2One +( + const dictionary& dict, + cfdemCloud& sm +) +: + dataExchangeModel(dict,sm), + propsDict_(dict.subDict(typeName + "Props")), + thisLigPartner_(0), + thisFoamPartner_(0), + lig2foam_(nullptr), + foam2lig_(nullptr), + lig2foam_mask_(nullptr), + lig2foam_ids_(nullptr), + foam2lig_ids_(nullptr), + lig2foam_vec_tmp_(nullptr), + lig2foam_scl_tmp_(nullptr), + foam2lig_vec_tmp_(nullptr), + foam2lig_scl_tmp_(nullptr), + lmp(nullptr) +{ + Info<<"Starting up LIGGGHTS for first time execution"<input->file(liggghtsPath.c_str()); + + // get DEM time step size + DEMts_ = lmp->update->dt; + checkTSsize(); + + createProcMap(); +} + +void twoWayOne2One::createProcMap() +{ + const cfdemCloud& sm = this->particleCloud_; + + // calculate boundingBox of FOAM subdomain + primitivePatch tmpBoundaryFaces + ( + SubList + ( + sm.mesh().faces(), + sm.mesh().nFaces() - sm.mesh().nInternalFaces(), + sm.mesh().nInternalFaces() + ), + sm.mesh().points() + ); + typedef PrimitivePatch bPatch; + bPatch boundaryFaces + ( + tmpBoundaryFaces.localFaces(), + tmpBoundaryFaces.localPoints() + ); + List foamBoxes(Pstream::nProcs()); + treeBoundBox thisFoamBox(boundaryFaces.localPoints()); + foamBoxes[Pstream::myProcNo()] = thisFoamBox; + Pstream::gatherList(foamBoxes); + Pstream::scatterList(foamBoxes); + + // calculate bounding box of LIG subdomain + // this may have to move to couple when dynamic LB occurs + List ligBoxes(Pstream::nProcs()); + double** ligbb = o2o_liggghts_get_boundingbox(lmp); + boundBox thisLigBox + ( + point(ligbb[0][0], ligbb[0][1], ligbb[0][2]), + point(ligbb[1][0], ligbb[1][1], ligbb[1][2]) + ); + ligBoxes[Pstream::myProcNo()] = thisLigBox; + Pstream::gatherList(ligBoxes); + Pstream::scatterList(ligBoxes); + + // detect LIG subdomains which this FOAM has to interact with + forAll(ligBoxes, ligproci) + { + if (thisFoamBox.overlaps(ligBoxes[ligproci])) + { + thisLigPartner_.append(ligproci); + } + } + // detect FOAM subdomains this LIG has to interact with + // TODO: refactor to invert this list here + forAll(foamBoxes, foamproci) + { + if (thisLigBox.overlaps(foamBoxes[foamproci])) + { + thisFoamPartner_.append(foamproci); + } + } +} + +// * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * // + +twoWayOne2One::~twoWayOne2One() +{ + delete foam2lig_; + delete lig2foam_; + + destroy(lig2foam_ids_); + destroy(foam2lig_ids_); + + destroy(lig2foam_vec_tmp_); + destroy(lig2foam_scl_tmp_); + destroy(foam2lig_vec_tmp_); + destroy(foam2lig_scl_tmp_); + + delete lmp; + MPI_Comm_free(&comm_liggghts_); +} + +// * * * * * * * * * * * * * * * public Member Functions * * * * * * * * * * * * * // +void twoWayOne2One::getData +( + word name, + word type, + double ** const& field, + label /*step*/ +) const +{ + if (name == "x") // the location is transferred by couple() + { + return; + } + if (type == "vector-atom") + { + double **tmp_= static_cast(lammps_extract_atom(lmp,name.c_str())); + if (!tmp_) + { + LAMMPS_NS::Fix *fix = nullptr; + fix = lmp->modify->find_fix_property + ( + name.c_str(), + "property/atom", + "vector", + 0, + 0, + "cfd coupling", + false + ); + if (fix) + { + tmp_ = static_cast(fix)->array_atom; + } + else + { + Warning<< "coupling fix not found!" << endl; + } + + if (!tmp_) + { + FatalError<< "find_fix_property " << name + << " array_atom not found." + << abort(FatalError); + } + } + + lig2foam_->exchange + ( + tmp_, + lig2foam_vec_tmp_, + 3 + ); + extractCollected + ( + lig2foam_vec_tmp_, + const_cast(field), + 3 + ); + } + else if (type == "scalar-atom") + { + double *tmp_ = static_cast(lammps_extract_atom(lmp,name.c_str())); + if (!tmp_) + { + LAMMPS_NS::Fix *fix = nullptr; + fix = lmp->modify->find_fix_property + ( + name.c_str(), + "property/atom", + "scalar", + 0, + 0, + "cfd coupling", + true + ); + + if (fix) + { + tmp_ = static_cast(fix)->vector_atom; + } + else + { + FatalError<< "coupling fix not found!" << abort(FatalError); + } + + if (!tmp_) + { + FatalError<< "find_fix_property " << name + << " vector_atom not found." + << abort(FatalError); + } + } + lig2foam_->exchange + ( + tmp_, + lig2foam_scl_tmp_ + ); + extractCollected + ( + lig2foam_scl_tmp_, + const_cast(field) + ); + } + else + { + FatalError << "requesting type " << type << " and name " << name << abort(FatalError); + } +} + +void twoWayOne2One::getData +( + word name, + word type, + int ** const& field, + label /*step*/ +) const +{ + FatalError << "do not use this getData!!!" << abort(FatalError); +/* + o2o_data_liggghts_to_of + ( + name.c_str(), + type.c_str(), + lmp, + (void*&) field, + "int", + comm_liggghts_ + ); +*/ +} + +void twoWayOne2One::giveData +( + word name, + word type, + double ** const& field, + const char* datatype +) const +{ + if (type == "vector-atom") + { + foam2lig_->exchange + ( + const_cast(field), + foam2lig_vec_tmp_, + 3 + ); + o2o_data_of_to_liggghts + ( + name.c_str(), + type.c_str(), + lmp, + foam2lig_vec_tmp_, + datatype, + foam2lig_ids_, + foam2lig_->ncollected_ + ); + } + else if (type == "scalar-atom") + { + foam2lig_->exchange + ( + const_cast(field), + foam2lig_scl_tmp_, + 1 + ); + o2o_data_of_to_liggghts + ( + name.c_str(), + type.c_str(), + lmp, + foam2lig_scl_tmp_, + datatype, + foam2lig_ids_, + foam2lig_->ncollected_ + ); + } + else + { + FatalError<< "twoWayMany2Many::giveData requested type " << type + << " not implemented!" + << abort(FatalError); + } +} + +//============ +// double ** +void twoWayOne2One::allocateArray +( + double**& array, + double initVal, + int width, + int length +) const +{ + int len = max(length,1); + lmp->memory->grow(array, len, width, "o2o:dbl**"); + for (int i = 0; i < len; i++) + for (int j = 0; j < width; j++) + array[i][j] = initVal; +} + +void twoWayOne2One::allocateArray +( + double**& array, + double initVal, + int width, + const char* length +) const +{ + int len = max(particleCloud_.numberOfParticles(),1); + lmp->memory->grow(array, len, width, "o2o:dbl**:autolen"); + for (int i = 0; i < len; i++) + for (int j = 0; j < width; j++) + array[i][j] = initVal; +} + +void inline twoWayOne2One::destroy(double** array,int len) const +{ + lmp->memory->destroy(array); +} + +//============ +// int ** +void twoWayOne2One::allocateArray +( + int**& array, + int initVal, + int width, + int length +) const +{ + int len = max(length,1); + lmp->memory->grow(array, len, width, "o2o:int**"); + for (int i = 0; i < len; i++) + for (int j = 0; j < width; j++) + array[i][j] = initVal; +} + +void twoWayOne2One::allocateArray +( + int**& array, + int initVal, + int width, + const char* length +) const +{ + int len = max(particleCloud_.numberOfParticles(),1); + lmp->memory->grow(array, len, width, "o2o:int**:autolen"); + for (int i = 0; i < len; i++) + for (int j = 0; j < width; j++) + array[i][j] = initVal; +} + +void inline twoWayOne2One::destroy(int** array,int len) const +{ + lmp->memory->destroy(array); +} + +//============ +// double * +void twoWayOne2One::allocateArray(double*& array, double initVal, int length) const +{ + int len = max(length,1); + lmp->memory->grow(array, len, "o2o:dbl*"); + for (int i = 0; i < len; i++) + array[i] = initVal; +} + +void inline twoWayOne2One::destroy(double* array) const +{ + lmp->memory->destroy(array); +} + +//============== +// int * +void twoWayOne2One::allocateArray(int*& array, int initVal, int length) const +{ + int len = max(length,1); + lmp->memory->grow(array, len, "o2o:int*"); + for (int i = 0; i < len; i++) + array[i] = initVal; +} + +void inline twoWayOne2One::destroy(int* array) const +{ + lmp->memory->destroy(array); +} +//============== + +bool twoWayOne2One::couple(int i) const +{ + bool coupleNow = false; + if (i==0) + { + couplingStep_++; + coupleNow = true; + + + // run commands from liggghtsCommands dict + Info<< "Starting up LIGGGHTS" << endl; + particleCloud_.clockM().start(3,"LIGGGHTS"); + + // check if liggghtsCommandModels with exaxt timing are being run + bool exactTiming(false); + int runComNr = -10; + DynamicList interruptTimes(0); + DynamicList DEMstepsToInterrupt(0); + DynamicList lcModel(0); + + forAll(particleCloud_.liggghtsCommandModelList(),i) + { + // Check if exact timing is needed + // get time for execution + // store time for execution in list + if(particleCloud_.liggghtsCommand()[i]().exactTiming()) + { + exactTiming = true; + DynamicList h + = particleCloud_.liggghtsCommand()[i]().executionsWithinPeriod + ( + TSstart(), + TSend() + ); + + forAll(h,j) + { + // save interrupt times (is this necessary) + interruptTimes.append(h[j]); + + // calc stepsToInterrupt + DEMstepsToInterrupt.append(DEMstepsTillT(h[j])); + + // remember which liggghtsCommandModel to run + lcModel.append(i); + } + + // make cumulative + label len = DEMstepsToInterrupt.size(); + label ind(0); + forAll(DEMstepsToInterrupt,i) + { + ind = len - i - 1; + if(ind > 0) + { + DEMstepsToInterrupt[ind] -= DEMstepsToInterrupt[ind-1]; + } + } + + Info<< "Foam::twoWayOne2One::couple(i): interruptTimes=" << interruptTimes << nl + << "Foam::twoWayOne2One::couple(i): DEMstepsToInterrupt=" << DEMstepsToInterrupt << nl + << "Foam::twoWayOne2One::couple(i): lcModel=" << lcModel + << endl; + } + + if(particleCloud_.liggghtsCommand()[i]().type() == "runLiggghts") + { + runComNr = i; + } + } + + // models with exact timing exists + label commandLines(0); + if(exactTiming) + { + // extension for more liggghtsCommands active the same time: + // sort interrupt list within this run period + // keep track of corresponding liggghtsCommand + int DEMstepsRun(0); + + forAll(interruptTimes,j) + { + // set run command till interrupt + DEMstepsRun += DEMstepsToInterrupt[j]; + particleCloud_.liggghtsCommand()[runComNr]().set(DEMstepsToInterrupt[j]); + const char* command = particleCloud_.liggghtsCommand()[runComNr]().command(0); + Info<< "Executing run command: '"<< command <<"'"<< endl; + lmp->input->one(command); + + // run liggghts command with exact timing + command = particleCloud_.liggghtsCommand()[lcModel[j]]().command(0); + Info << "Executing command: '"<< command <<"'"<< endl; + lmp->input->one(command); + } + + // do the run + if(particleCloud_.liggghtsCommand()[runComNr]().runCommand(couplingStep())) + { + particleCloud_.liggghtsCommand()[runComNr]().set(couplingInterval() - DEMstepsRun); + const char* command = particleCloud_.liggghtsCommand()[runComNr]().command(0); + Info<< "Executing run command: '"<< command <<"'"<< endl; + lmp->input->one(command); + } + + // do the other non exact timing models + forAll(particleCloud_.liggghtsCommandModelList(),i) + { + if + ( + ! particleCloud_.liggghtsCommand()[i]().exactTiming() && + particleCloud_.liggghtsCommand()[i]().runCommand(couplingStep()) + ) + { + commandLines=particleCloud_.liggghtsCommand()[i]().commandLines(); + for(int j=0;jinput->one(command); + } + } + } + } + // no model with exact timing exists + else + { + forAll(particleCloud_.liggghtsCommandModelList(),i) + { + if(particleCloud_.liggghtsCommand()[i]().runCommand(couplingStep())) + { + commandLines=particleCloud_.liggghtsCommand()[i]().commandLines(); + for(int j=0;jinput->one(command); + } + } + } + } + + particleCloud_.clockM().stop("LIGGGHTS"); + Info<< "LIGGGHTS finished" << endl; + + setupLig2FoamCommunication(); + + locateParticles(); + + setupFoam2LigCommunication(); + } + + return coupleNow; +} + +void twoWayOne2One::setupLig2FoamCommunication() const +{ + int* src_procs = new int[thisLigPartner_.size()]; + for (int proci = 0; proci < thisLigPartner_.size(); proci++) + { + src_procs[proci] = thisLigPartner_[proci]; + } + int* dst_procs = new int[thisFoamPartner_.size()]; + for (int proci = 0; proci < thisFoamPartner_.size(); proci++) + { + dst_procs[proci] = thisFoamPartner_[proci]; + } + + delete lig2foam_; + lig2foam_ = new One2One(comm_liggghts_); + lig2foam_->setup + ( + thisLigPartner_.size(), + src_procs, + thisFoamPartner_.size(), + dst_procs, + lmp->atom->nlocal + ); + allocateArray + ( + lig2foam_vec_tmp_, + 0., + 3 * lig2foam_->ncollected_ + ); + allocateArray + ( + lig2foam_scl_tmp_, + 0., + lig2foam_->ncollected_ + ); +} + + +void twoWayOne2One::locateParticles() const +{ + // get positions for locate + double** my_positions = static_cast(lmp->atom->x); + double* my_flattened_positions = nullptr; + allocateArray(my_flattened_positions, 0., 3*lmp->atom->nlocal); + for (int atomi = 0; atomi < lmp->atom->nlocal; atomi++) + for (int coordi = 0; coordi < 3; coordi++) + my_flattened_positions[atomi*3+coordi] = my_positions[atomi][coordi]; + + double* collected_flattened_positions = nullptr; + allocateArray(collected_flattened_positions, 0., 3*lig2foam_->ncollected_); + + lig2foam_->exchange(my_flattened_positions, collected_flattened_positions, 3); + destroy(my_flattened_positions); + + if (lig2foam_mask_) + { + delete [] lig2foam_mask_; + } + lig2foam_mask_ = new bool[lig2foam_->ncollected_]; + + labelList cellIds(0); + label n_located(0); + for (int atomi = 0; atomi < lig2foam_->ncollected_; atomi++) + { + const vector position = vector + ( + collected_flattened_positions[3*atomi+0], + collected_flattened_positions[3*atomi+1], + collected_flattened_positions[3*atomi+2] + ); + + const label cellI = particleCloud_.locateM().findSingleCell(position, -1); + + lig2foam_mask_[atomi] = false; + if (cellI >= 0) // in domain + { + lig2foam_mask_[atomi] = true; + n_located++; + cellIds.append(cellI); + } + } + + setNumberOfParticles(n_located); + particleCloud_.reAllocArrays(); + + // copy positions/cellids/ids of located particles into arrays + allocateArray(lig2foam_ids_, 0, getNumberOfParticles()); + int* collected_ids = nullptr; + allocateArray(collected_ids, 0, lig2foam_->ncollected_); + lig2foam_->exchange(lmp->atom->tag, collected_ids); + extractCollected(collected_ids, lig2foam_ids_); + destroy(collected_ids); + + double* extracted_flattened_positions = new double[getNumberOfParticles()*3]; + extractCollected + ( + collected_flattened_positions, + extracted_flattened_positions, + 3 + ); + setPositions(getNumberOfParticles(), extracted_flattened_positions); + destroy(extracted_flattened_positions); + destroy(collected_flattened_positions); + + setCellIDs(cellIds); +} + +void twoWayOne2One::setupFoam2LigCommunication() const +{ + int* src_procs = new int[thisFoamPartner_.size()]; + for (int proci = 0; proci < thisFoamPartner_.size(); proci++) + { + src_procs[proci] = thisFoamPartner_[proci]; + } + + int* dst_procs = new int[thisLigPartner_.size()]; + for (int proci = 0; proci < thisLigPartner_.size(); proci++) + { + dst_procs[proci] = thisLigPartner_[proci]; + } + + delete foam2lig_; + foam2lig_ = new One2One(comm_liggghts_); + + foam2lig_->setup + ( + thisFoamPartner_.size(), + src_procs, + thisLigPartner_.size(), + dst_procs, + getNumberOfParticles() + ); + allocateArray + ( + foam2lig_ids_, + 0, + foam2lig_->ncollected_ + ); + foam2lig_->exchange(lmp->atom->tag, foam2lig_ids_); + + allocateArray + ( + foam2lig_vec_tmp_, + 1., + 3 * foam2lig_->ncollected_ + ); + allocateArray + ( + foam2lig_scl_tmp_, + 1., + foam2lig_->ncollected_ + ); +} + +template +void twoWayOne2One::extractCollected(T**& src, T**& dst, int width) const +{ + int locali = 0; + + for (int atomi = 0; atomi < lig2foam_->ncollected_; atomi++) + { + if (!lig2foam_mask_[atomi]) continue; + + for (int coordi = 0; coordi < width; coordi++) + { + dst[locali][coordi] = src[atomi][coordi] ; + } + locali++; + } +} + +template +void twoWayOne2One::extractCollected(T*& src, T*& dst, int width) const +{ + int locali = 0; + + for (int atomi = 0; atomi < lig2foam_->ncollected_; atomi++) + { + if (!lig2foam_mask_[atomi]) continue; + + for (int coordi = 0; coordi < width; coordi++) + { + dst[locali] = src[atomi*width+coordi]; + locali++; + } + } +} + +template +void twoWayOne2One::extractCollected(T*& src, T**& dst, int width) const +{ + int locali = 0; + + for (int atomi = 0; atomi < lig2foam_->ncollected_; atomi++) + { + if (!lig2foam_mask_[atomi]) continue; + + for (int coordi = 0; coordi < width; coordi++) + { + dst[locali][coordi] = src[atomi*width+coordi]; + } + locali++; + } +} + +int twoWayOne2One::getNumberOfParticles() const +{ + return particleCloud_.numberOfParticles(); +} + + + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +} // End namespace Foam + +// ************************************************************************* // diff --git a/src/lagrangian/cfdemParticle/subModels/dataExchangeModel/twoWayOne2One/twoWayOne2One.H b/src/lagrangian/cfdemParticle/subModels/dataExchangeModel/twoWayOne2One/twoWayOne2One.H new file mode 100644 index 00000000..0de31a21 --- /dev/null +++ b/src/lagrangian/cfdemParticle/subModels/dataExchangeModel/twoWayOne2One/twoWayOne2One.H @@ -0,0 +1,226 @@ +/*---------------------------------------------------------------------------*\ + CFDEMcoupling - Open Source CFD-DEM coupling + + CFDEMcoupling is part of the CFDEMproject + www.cfdem.com + Christoph Goniva, christoph.goniva@cfdem.com + Copyright 2009-2012 JKU Linz + Copyright 2012- DCS Computing GmbH, Linz +------------------------------------------------------------------------------- +License + This file is part of CFDEMcoupling. + + CFDEMcoupling 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. + + CFDEMcoupling 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 CFDEMcoupling; if not, write to the Free Software Foundation, + Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA + +Description + enhanced two way DEM-CFD coupling via MPI. + + Compared to twoWayMPI, no Allreduces are used for communication. + Instead, a geomatric map between FOAM and LIG domains is created and + subsequently used for communication. + +Class + twoWayOne2One + +SourceFiles + twoWayOne2One.C + +Contributing authors + Paul Kieckhefen (TUHH) 2018 + +\*---------------------------------------------------------------------------*/ + +#ifndef twoWayOne2One_H +#define twoWayOne2One_H + +#include "dataExchangeModel.H" +#include "liggghtsCommandModel.H" +#include "OFstream.H" +#include +#include "pair.h" +#include "force.h" +#include "forceModel.H" +#include "one2one.H" + +#include "meshSearch.H" + +//=================================// +//LAMMPS/LIGGGHTS +#include +#include +#include +#include +#include // these are LAMMPS include files +#include +#include +#include +#include +#include +#include +#include +//=================================// + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +namespace Foam +{ + +/*---------------------------------------------------------------------------*\ + Class noDrag Declaration +\*---------------------------------------------------------------------------*/ + +class twoWayOne2One +: + public dataExchangeModel +{ +private: + + // private data + dictionary propsDict_; + + mutable MPI_Comm comm_liggghts_; + + + // LIG ranks from which to retrieve particle data + labelList thisLigPartner_; + labelList thisFoamPartner_; + + mutable One2One* lig2foam_; + mutable One2One* foam2lig_; + + mutable bool* lig2foam_mask_; + + mutable int* lig2foam_ids_; + mutable int* foam2lig_ids_; + + mutable double* lig2foam_vec_tmp_; + mutable double* lig2foam_scl_tmp_; + + mutable double* foam2lig_vec_tmp_; + mutable double* foam2lig_scl_tmp_; + + // private member functions + + //- creates a geometric mapping between FOAM and LIG domains + void createProcMap(); + + //- create a One2One communicator which transfers from LIG to FOAM + void setupLig2FoamCommunication() const; + + //- locates particles received from Lig + void locateParticles() const; + + //- create a One2One communicator which transfers from FOAM to LIG + void setupFoam2LigCommunication() const; + +protected: + LAMMPS_NS::LAMMPS *lmp; + +public: + + //- Runtime type information + TypeName("twoWayOne2One"); + + + // Constructors + + //- Construct from components + twoWayOne2One + ( + const dictionary& dict, + cfdemCloud& sm + ); + + // Destructor + + ~twoWayOne2One(); + + + // Member Functions + void getData + ( + word name, + word type, + double ** const& field, + label step + ) const; + + void getData + ( + word name, + word type, + int ** const& field, + label step + ) const; + + void giveData + ( + word name, + word type, + double ** const& field, + const char* datatype + ) const; + + //============ + // double ** + void allocateArray(double**&, double, int, int) const; + void allocateArray(double**&, double, int,const char* = "nparticles") const; + void inline destroy(double**,int=0) const; + //============ + // int ** + void allocateArray(int**&, int, int, int) const; + void allocateArray(int**&, int, int,const char* = "nparticles") const; + void inline destroy(int**,int=0) const; + //============== + + //============== + // double * + void allocateArray(double*&, double, int) const; + void inline destroy(double*) const; + //============== + // int * + void allocateArray(int*&, int, int) const; + void inline destroy(int*) const; + //============== + + bool couple(int) const; + + //- extractCollected takes the collected data from Lig + // present in this Foam domain and applies the mask. + // the width parameter can be used for reshaping. + template + void extractCollected(T**&, T**&, int width=1) const; + template + void extractCollected(T*&, T*&, int width=1) const; + template + void extractCollected(T*&, T**&, int width=1) const; + + int getNumberOfParticles() const; + + word myType() const { return typeName; } + + void setCG() { particleCloud_.setCG(lmp->force->cg()); } +}; + + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +} // End namespace Foam + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +#endif + +// ************************************************************************* // diff --git a/src/lagrangian/cfdemParticle/subModels/locateModel/engineSearchMany2Many/engineSearchMany2Many.C b/src/lagrangian/cfdemParticle/subModels/locateModel/engineSearchMany2Many/engineSearchMany2Many.C index 0df1959b..7c955e90 100755 --- a/src/lagrangian/cfdemParticle/subModels/locateModel/engineSearchMany2Many/engineSearchMany2Many.C +++ b/src/lagrangian/cfdemParticle/subModels/locateModel/engineSearchMany2Many/engineSearchMany2Many.C @@ -95,7 +95,16 @@ label engineSearchMany2Many::intersection return face; } - +label engineSearchMany2Many::findCell +( + double** const& mask, + double**& positions, + int**& cellIDs, + int size +) const +{ + return 1; // locate is provided by many2may / one2one +} // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // } // End namespace Foam diff --git a/src/lagrangian/cfdemParticle/subModels/locateModel/engineSearchMany2Many/engineSearchMany2Many.H b/src/lagrangian/cfdemParticle/subModels/locateModel/engineSearchMany2Many/engineSearchMany2Many.H index 5fe4ffce..ba46b670 100755 --- a/src/lagrangian/cfdemParticle/subModels/locateModel/engineSearchMany2Many/engineSearchMany2Many.H +++ b/src/lagrangian/cfdemParticle/subModels/locateModel/engineSearchMany2Many/engineSearchMany2Many.H @@ -87,6 +87,13 @@ public: const point& pEnd ) const; + label findCell + ( + double** const& mask, + double**& positions, + int**& cellIDs, + int size + ) const; }; From e1472f5ed33f956178a95ab1135606c3609509b8 Mon Sep 17 00:00:00 2001 From: Paul Kieckhefen Date: Sat, 10 Mar 2018 17:48:25 +0100 Subject: [PATCH 02/29] moved procmap creation temporarily moved to couple(). debug output added. --- .../twoWayOne2One/twoWayOne2One.C | 31 +++++++++++++++---- .../twoWayOne2One/twoWayOne2One.H | 8 ++--- 2 files changed, 29 insertions(+), 10 deletions(-) diff --git a/src/lagrangian/cfdemParticle/subModels/dataExchangeModel/twoWayOne2One/twoWayOne2One.C b/src/lagrangian/cfdemParticle/subModels/dataExchangeModel/twoWayOne2One/twoWayOne2One.C index 668dc9e2..0754473b 100644 --- a/src/lagrangian/cfdemParticle/subModels/dataExchangeModel/twoWayOne2One/twoWayOne2One.C +++ b/src/lagrangian/cfdemParticle/subModels/dataExchangeModel/twoWayOne2One/twoWayOne2One.C @@ -102,11 +102,9 @@ twoWayOne2One::twoWayOne2One // get DEM time step size DEMts_ = lmp->update->dt; checkTSsize(); - - createProcMap(); } -void twoWayOne2One::createProcMap() +void twoWayOne2One::createProcMap() const { const cfdemCloud& sm = this->particleCloud_; @@ -146,6 +144,9 @@ void twoWayOne2One::createProcMap() Pstream::gatherList(ligBoxes); Pstream::scatterList(ligBoxes); + thisLigPartner_.clear(); + thisFoamPartner_.clear(); + // detect LIG subdomains which this FOAM has to interact with forAll(ligBoxes, ligproci) { @@ -618,6 +619,8 @@ bool twoWayOne2One::couple(int i) const particleCloud_.clockM().stop("LIGGGHTS"); Info<< "LIGGGHTS finished" << endl; + createProcMap(); + setupLig2FoamCommunication(); locateParticles(); @@ -673,8 +676,12 @@ void twoWayOne2One::locateParticles() const double* my_flattened_positions = nullptr; allocateArray(my_flattened_positions, 0., 3*lmp->atom->nlocal); for (int atomi = 0; atomi < lmp->atom->nlocal; atomi++) + { for (int coordi = 0; coordi < 3; coordi++) + { my_flattened_positions[atomi*3+coordi] = my_positions[atomi][coordi]; + } + } double* collected_flattened_positions = nullptr; allocateArray(collected_flattened_positions, 0., 3*lig2foam_->ncollected_); @@ -713,6 +720,18 @@ void twoWayOne2One::locateParticles() const setNumberOfParticles(n_located); particleCloud_.reAllocArrays(); + Pout<< "Located " << n_located << " out of " << lig2foam_->ncollected_ << " . " + << endl; + + reduce(n_located, sumOp