DuaSPHysics Cell Division Code Structure

Try to straight up the cell division part of DuaSPHysics's tricky code. A structure map is shown.

Class inheriting

JSphGpuSingle::Run

ConfigCellDivision();

void JSph::ConfigCellDivision()

Map_Size: Maximum number of cells within case limits, =KernelSize/ScellDiv (KernelSize or KernelSize/2).

Map_Cells: Maximum number of cells within case limits. =TUint3(unsigned(ceil(Map_Size.xyz/Scell)

Output CfgInit_MapCells.vtk looks like blow:

SelecDomain(TUint3(0,0,0),Map_Cells);

void JSph::SelecDomain(tuint3 celini,tuint3 celfin)

sets local domain of simulation within Map_Cells

Cal these parms below (defined in class JSph):

DomCells: actual Number of cells in each direction

DomPosMin: Lower limit of simulation + edge (KernelSize) if periodic conditions. DomPosMin=Map_PosMin+(DomCelIni*Scell);

DomPosMax: Upper limit of simulation + edge (KernelSize) if periodic conditions.

DomPosMax=min(Map_PosMax,Map_PosMin+(DomCelFin*Scell));

DomRealPosMin/DomRealPosMax: actual limits of local simulation according to DomCelIni/Fin (without periodic condition borders)

DomRealPosMin=max(DomPosMin,MapRealPosMin)

DomRealPosMax=min(DomPosMax,MapRealPosMax)

DomCellCode: Indicator of group of particles & other special markers.

LoadDcellParticles(Np,Code,AuxPos,Dcell);

void JSph::LoadDcellParticles(unsigned n,const typecode *code,const tdouble3 *pos,unsigned *dcell)const{

traverse all the particles, computes the initial cell, store in dcell

checks if there are unexpected excluded particles

cx,cy,cz: the cell index of the particle belongs to,

dx=pos[p].x-DomPosMin.x, dy=pos[p].y-DomPosMin.y, double dz=pos[p].z-DomPosMin.

cx=dx/Scell, cy=dy/Scell, cz=dz/Scell;

dcell[p]=PC__Cell(DomCellCode,cx,cy,cz);

PC__Cell packs all three indexes into a single integer in different method according to different DomCellCode

#define PC__Cell(cc,cx,cy,cz) ((cx<<((cc>>10)&31))|(cy<<((cc>>15)&31))|cz)  //-Cell value for cx,cy and cz.

CellDivSingle = new JCellDivGpuSingle

Note that CellDivSingle is a member variable of class JSphGpuSingle. Initilized in the constructor. Deleted in the destructor.

It news at the above line.

And all the configs should be passed in because of the class inheriting

class JSphGpuSingle : public JSphGpu
{
protected:
  JCellDivGpuSingle* CellDivSingle;
  // ...
    
public:
  JSphGpuSingle();
  ~JSphGpuSingle();
  void Run(std::string appname,const JSphCfgRun *cfg,JLog2 *log);

  // ...
}

JSphGpuSingle::JSphGpuSingle():JSphGpu(false){
  ClassName="JSphGpuSingle";
  CellDivSingle=NULL;
}

JSphGpuSingle::~JSphGpuSingle(){
  DestructorActive=true;
  delete CellDivSingle; CellDivSingle=NULL;
}

CellDivSingle->Divide(Npb,Np-Npb-NpbPer-NpfPer,NpbPer,NpfPer,BoundChanged,Dcellg,Codeg,Timers,Posxyg,Poszg,Idpg);

void JCellDivGpuSingle::Divide(
    unsigned npb1  		// Npb;   *///< Number of boundary particles*
    unsigned npf1 		// Np-Npb-NpbPer-NpfPer
                        // unsigned Np;    *///< Total number of particles*
                        // unsigned NpbPer;   *///< Number of periodic boundary particles.*
                        // unsigned NpfPer;   *///< Number of periodic particles (fluid-floating).* 
    unsigned npb2 		// unsigned NpbPer;   *///< Number of periodic boundary particles.*
    unsigned npf2 		// unsigned NpfPer;   *///< Number of periodic particles (fluid-floating).* 
    bool boundchanged; 	// bool BoundChanged   *///< Indicates if a selected boundary particle has changed since the last time step.*
    const unsigned*dcellg 	// unsigned *Dcellg; *///< Cells inside DomCells coded with DomCellCode.*
    const typecode*codeg 	// typecode *Codeg;  *///< Indicator of group of particles & other special markers.*
    TimersGpu timers
    const double2 *posxy  	// double2 *Posxyg;
    const double *posz   	// double *Poszg;
    const unsigned *idp 	// unsigned *Idpg;   *///< Identifier of particle*
)

Initial processing of Divide: Calculate the limits of the domain and compute the new position of each particle (SortPart).

The value for np includes periodic boundary and fluid particles (npbper and npfper).

The floating bodies are treated as fluids (both to be ignored as excluded).

cudiv::Sort

void Sort(unsigned* keys,unsigned* values,unsigned size,bool stable){
  if(size){
    thrust::device_ptr<unsigned> dev_keysg(keys);
    thrust::device_ptr<unsigned> dev_valuesg(values);
    if(stable)thrust::stable_sort_by_key(dev_keysg,dev_keysg+size,dev_valuesg);
    else thrust::sort_by_key(dev_keysg,dev_keysg+size,dev_valuesg);
  }
}

Last arg always stable

Different key, values, list length for if(DivideFull)