yade-dev team mailing list archive
-
yade-dev team
-
Mailing list archive
-
Message #10707
[Branch ~yade-pkg/yade/git-trunk] Rev 3906: Fix conflict.
Merge authors:
Luc Sibille (luc-sibille)
------------------------------------------------------------
revno: 3906 [merge]
committer: Anton Gladky <gladky.anton@xxxxxxxxx>
timestamp: Wed 2014-04-09 21:36:24 +0200
message:
Fix conflict.
added:
pkg/lbm/
pkg/lbm/HydrodynamicsLawLBM.cpp
pkg/lbm/HydrodynamicsLawLBM.hpp
pkg/lbm/LBMbody.cpp
pkg/lbm/LBMbody.hpp
pkg/lbm/LBMlink.cpp
pkg/lbm/LBMlink.hpp
pkg/lbm/LBMnode.cpp
pkg/lbm/LBMnode.hpp
modified:
CMakeLists.txt
--
lp:yade
https://code.launchpad.net/~yade-pkg/yade/git-trunk
Your team Yade developers is subscribed to branch lp:yade.
To unsubscribe from this branch go to https://code.launchpad.net/~yade-pkg/yade/git-trunk/+edit-subscription
=== modified file 'CMakeLists.txt'
--- CMakeLists.txt 2014-04-09 09:11:49 +0000
+++ CMakeLists.txt 2014-04-09 19:36:24 +0000
@@ -15,6 +15,7 @@
# ENABLE_LINSOLV: enable LINSOLV-option (ON by default)
# ENABLE_PFVFLOW: enable PFVFLOW-option, FlowEngine (ON by default)
# ENABLE_SPH: enable SPH-option, Smoothed Particle Hydrodynamics (OFF by default, experimental)
+# ENABLE_LBMFLOW: enable LBMFLOW-option, LBM_ENGINE (OFF by default)
# runtimePREFIX: used for packaging, when install directory is not the same is runtime directory (/usr/local by default)
# CHUNKSIZE: set >1, if you want several sources to be compiled at once. Increases compilation speed and RAM-consumption during it (1 by default).
@@ -118,6 +119,7 @@
OPTION(ENABLE_LINSOLV "Enable direct solver for the flow engines (experimental)" ${DEFAULT})
OPTION(ENABLE_PFVFLOW "Enable flow engine (experimental)" ${DEFAULT})
OPTION(ENABLE_SPH "Enable SPH" OFF)
+OPTION(ENABLE_LBMFLOW "Enable LBM engine (very experimental)" OFF)
#===========================================================
# Use Eigen3 by default
@@ -317,6 +319,15 @@
INCLUDE_DIRECTORIES(${CMAKE_BINARY_DIR})
#===========================================================
+IF(ENABLE_LBMFLOW)
+ SET(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -DLBM_ENGINE")
+ SET(CONFIGURED_FEATS "${CONFIGURED_FEATS} LBMFLOW")
+ MESSAGE("LBMFLOW is still experimental, building and running LBM engine are at your own risk!")
+ELSE(ENABLE_LBMFLOW)
+ SET(DISABLED_FEATS "${DISABLED_FEATS} LBMFLOW")
+ENDIF(ENABLE_LBMFLOW)
+
+#===========================================================
IF (NOT INSTALL_PREFIX)
SET(CMAKE_INSTALL_PREFIX "/usr/local")
=== added directory 'pkg/lbm'
=== added file 'pkg/lbm/HydrodynamicsLawLBM.cpp'
--- pkg/lbm/HydrodynamicsLawLBM.cpp 1970-01-01 00:00:00 +0000
+++ pkg/lbm/HydrodynamicsLawLBM.cpp 2014-03-31 16:01:57 +0000
@@ -0,0 +1,1370 @@
+/*************************************************************************
+* Copyright (C) 2009-2012 by Franck Lominé *
+* franck.lomine@xxxxxxxxxxxxxx *
+* *
+* This program is free software; it is licensed under the terms of the *
+* GNU General Public License v2 or later. See file LICENSE for details. *
+* *
+* Luc Scholtès luc.scholtes@xxxxxxxxxxxxxxxx *
+* and Luc Sibille luc.sibille@xxxxxxxxxxxxxxx also contributed to this *
+* code. *
+* *
+* Lominé F., Scholtès L., Sibille L., Poullain P. (2013) *
+* Modelling of fluid-solid interaction in granular media with coupled *
+* LB/DE methods: application to piping erosion. International Journal *
+* for Numerical and Analytical Methods in Geomechanics, 37(6):577-596 *
+* doi: 10.1002/nag.1109 *
+* *
+* Sibille L., Lominé F., Marot D. (2012) Investigation In Modelling *
+* Piping Erosion With a Coupled «Lattice Boltzmann – Discrete Element» *
+* Numerical Method. in Proc. 6th Int. Conference on Scour and Erosion *
+* (ICSE-6), pp. 93-100. *
+* *
+*************************************************************************/
+#ifdef LBM_ENGINE
+
+#include"HydrodynamicsLawLBM.hpp"
+#include<yade/core/Omega.hpp>
+#include<yade/core/Scene.hpp>
+#include<boost/filesystem/operations.hpp>
+#include<yade/pkg/common/Box.hpp>
+#include<yade/pkg/common/Sphere.hpp>
+#include<yade/pkg/dem/CohFrictPhys.hpp>
+
+
+namespace bfs=boost::filesystem;
+using namespace std;
+
+template<class Scalar> VECTOR3_TEMPLATE(Scalar) operator*(Scalar s, const VECTOR3_TEMPLATE(Scalar)& v) {return v*s;}
+inline Vector3i vect3rToVect3i(Vector3r vect){Vector3i newvect((int)vect[0],(int)vect[1],(int)vect[2]);return(newvect);}
+
+
+HydrodynamicsLawLBM::~HydrodynamicsLawLBM() {};
+
+bool HydrodynamicsLawLBM::isActivated(){
+ DEM_ITER=scene->iter;//+1;
+ if(EngineIsActivated){
+ if((DEM_ITER % DemIterLbmIterRatio==0)&&(DEM_ITER!=DemIterLbmIterRatio)) {
+ if(DEM_ITER==0){DEMdt0 = scene->dt;scene->dt=1.e-50;}
+ return(true);}
+ else{
+ if(applyForcesAndTorques) CalculateAndApplyForcesAndTorquesOnBodies(false,true);
+ return(false);
+ }
+ }else return(false);
+}
+
+void HydrodynamicsLawLBM::action()
+{
+ timingDeltas->start();
+ NB_BODIES= scene->bodies->size();
+ int I, J, step=0;
+ NbFluidNodes=0;
+ NbSolidNodes=0;
+ Real CurMinVelOfPtc=1000000.;
+ Real CurMaxVelOfPtc=-1000000.;
+ int ErrorCriterion=2;
+ /*------------------------------------------------------------------*/
+ /* AT FIRST ITERATION */
+ /*------------------------------------------------------------------*/
+ if(firstRun){
+ // createNewFiles(); //this line is move further to create files only when the recording configuration chosen by the operator is known
+ bool initVbCutOff=false;
+ if(VbCutOff==-1) initVbCutOff=true;
+ halfWallthickness=1000000.;
+ NB_WALLS=0;
+ //NB_DYNWALLS=0;
+ NB_DYNGRAINS=0;
+ LBMbody tmpbody;
+ FOREACH(const shared_ptr<Body>& b, *scene->bodies){
+ if(!b) continue; // deleted bodies
+ if (b->shape->getClassName()=="Box"){
+ Vector3r ext(YADE_PTR_CAST<Box> ( b->shape )->extents);
+ if (ext[0]<halfWallthickness) halfWallthickness=ext[0];
+ if (ext[1]<halfWallthickness) halfWallthickness=ext[1];
+ if (ext[2]<halfWallthickness) halfWallthickness=ext[2];
+ NB_WALLS++;
+ tmpbody.setAsBox();
+ tmpbody.saveProperties=false;
+ //if(b->isDynamic()) NB_DYNWALLS++;
+ }
+ if (b->shape->getClassName()=="Sphere"){
+ const shared_ptr<Sphere>& sph = YADE_PTR_CAST<Sphere> ( b->shape );
+ if(IdFirstSphere==-1) IdFirstSphere=b->getId();
+ tmpbody.setAsPtc();
+ Real r=sph->radius;
+ if(b->isDynamic()){
+ NB_DYNGRAINS++;
+ /*--- computation of the initial volume ---*/
+ if(!strcmp(model.c_str(), "d2q9" )) Vo += Mathr::PI*(r*r);
+ else Vo += 4./3.*Mathr::PI*(r*r*r);
+ if(initVbCutOff) VbCutOff=max(VbCutOff,b->state->vel.norm()+r*b->state->angVel.norm());
+ tmpbody.saveProperties=true;
+ }else{
+ if(b->state->pos[1]>0.) tmpbody.saveProperties=true;
+ else tmpbody.saveProperties=false;
+ }
+ /*--- computation of Rmean Rmax Rmin --*/
+ MaxBodyRadius=max(r,MaxBodyRadius);
+ MinBodyRadius=min(r,MinBodyRadius);
+ MeanBodyRadius+=r;
+ }
+
+ LBbodies.push_back(tmpbody);
+ }
+ Wallthickness=2.0*halfWallthickness;
+ NB_GRAINS=NB_BODIES-NB_WALLS;
+ //Luc: is it right to count the walls as dynamic bodies?
+ //Franck why not ? Enhancement (see next line is coming ;-))
+ NB_DYNBODIES=NB_WALLS+NB_DYNGRAINS;
+ //NB_DYNBODIES=NB_DYNWALLS+NB_DYNGRAINS;
+ MeanBodyRadius=MeanBodyRadius/NB_GRAINS;
+ InitialNumberOfDynamicParticles=NB_DYNGRAINS;
+
+ /*-------------------------------------------------------------------------*/
+ /* D2Q9 model configuration */
+ /*-------------------------------------------------------------------------*/
+ if(!strcmp(model.c_str(), "d2q9" )){
+ dim=2;
+ /*--------------------------------------*/
+ /* D2Q9 model: 6 2 5 */
+ /* 3 _\|/__1 */
+ /* /|\ */
+ /* 7 4 8 */
+ /*--------------------------------------*/
+ /*----------- D2Q9 constants ---------*/
+ w.push_back(4.0/9.0);
+ for(int aa=1;aa<=4;aa++) w.push_back(1.0/9.0);
+ for(int aa=5;aa<=8;aa++) w.push_back(1.0/36.0);
+
+ /*--------- node position vectors ----*/
+ eib.push_back(Vector3r( 0., 0., 0.)); //0
+ eib.push_back(Vector3r( 1., 0., 0.)); //1
+ eib.push_back(Vector3r( 0., 1., 0.)); //2
+ eib.push_back(Vector3r(-1., 0., 0.)); //3
+ eib.push_back(Vector3r( 0.,-1., 0.)); //4
+ eib.push_back(Vector3r( 1., 1., 0.)); //5
+ eib.push_back(Vector3r(-1., 1., 0.)); //6
+ eib.push_back(Vector3r(-1.,-1., 0.)); //7
+ eib.push_back(Vector3r( 1.,-1., 0.)); //8
+ NbDir=(int) eib.size();
+
+ /*-------------- opposite nodes --------*/
+ opp.push_back(0); opp.push_back(3); opp.push_back(4);
+ opp.push_back(1); opp.push_back(2); opp.push_back(7);
+ opp.push_back(8); opp.push_back(5); opp.push_back(6);
+ }else {cerr<<"This model is not implemented yet: "<<model<<endl;exit(-1);}
+
+
+ int res=0;
+ int ll=0;
+ bool CreateLbmDir = false; //Flag to create directories if required by the recording configuration
+ bool CreateDemDir = false; //Flag to create directories if required by the recording configuration
+ bool CreateCntctDir = false; //Flag to create directories if required by the recording configuration
+ /*-------------------------------------------------------------------------*/
+ /* Recording configuration */
+ /*-------------------------------------------------------------------------*/
+ ll=LBMSavedData.size();
+ res=LBMSavedData.find("Velocity"); if(res>=0&&res<ll) {SAVE_VELOCITY =true; SAVE_VELOCITYCOMP=true; CreateLbmDir=true;}
+ res=LBMSavedData.find("velocity"); if(res>=0&&res<ll) {SAVE_VELOCITY =true; SAVE_VELOCITYCOMP=true; CreateLbmDir=true;}
+ res=LBMSavedData.find("VelXY"); if(res>=0&&res<ll) {SAVE_VELOCITYCOMP=true; CreateLbmDir=true;}
+ res=LBMSavedData.find("velXY"); if(res>=0&&res<ll) {SAVE_VELOCITYCOMP=true; CreateLbmDir=true;}
+ res=LBMSavedData.find("Rho"); if(res>=0&&res<ll) {SAVE_RHO =true; CreateLbmDir=true;}
+ res=LBMSavedData.find("rho"); if(res>=0&&res<ll) {SAVE_RHO =true; CreateLbmDir=true;}
+ res=LBMSavedData.find("Forces"); if(res>=0&&res<ll) {SAVE_FORCES =true; CreateLbmDir=true;}
+ res=LBMSavedData.find("forces"); if(res>=0&&res<ll) {SAVE_FORCES =true; CreateLbmDir=true;}
+ res=LBMSavedData.find("Bodies"); if(res>=0&&res<ll) {SAVE_BODIES =true; CreateLbmDir=true;}
+ res=LBMSavedData.find("bodies"); if(res>=0&&res<ll) {SAVE_BODIES =true; CreateLbmDir=true;}
+ res=LBMSavedData.find("NodeBD"); if(res>=0&&res<ll) {SAVE_NODEBD =true; CreateLbmDir=true;}
+ res=LBMSavedData.find("nodeBD"); if(res>=0&&res<ll) {SAVE_NODEBD =true; CreateLbmDir=true;}
+ res=LBMSavedData.find("NewNode"); if(res>=0&&res<ll) {SAVE_NODEISNEW =true; CreateLbmDir=true;}
+ res=LBMSavedData.find("newNode"); if(res>=0&&res<ll) {SAVE_NODEISNEW =true; CreateLbmDir=true;}
+ res=LBMSavedData.find("bz2"); if(res>=0&&res<ll) COMPRESS_DATA =true;
+ res=LBMSavedData.find("ObservedPtc"); if(res>=0&&res<ll) SAVE_OBSERVEDPTC =true;
+ res=LBMSavedData.find("observedptc"); if(res>=0&&res<ll) SAVE_OBSERVEDPTC =true;
+ res=LBMSavedData.find("observedPtc"); if(res>=0&&res<ll) SAVE_OBSERVEDPTC =true;
+ res=LBMSavedData.find("ObservedNode"); if(res>=0&&res<ll) SAVE_OBSERVEDNODE=true; // does not activate any recording in a a file currently
+ res=LBMSavedData.find("observednode"); if(res>=0&&res<ll) SAVE_OBSERVEDNODE=true;
+ res=LBMSavedData.find("observedNode"); if(res>=0&&res<ll) SAVE_OBSERVEDNODE=true;
+ res=LBMSavedData.find("contacts"); if(res>=0&&res<ll) {SAVE_CONTACTINFO =true; CreateCntctDir=true;} // SAVE_CONTACTINFO does not activate any (anymore) recording excepted the creation of a directory
+ res=LBMSavedData.find("Contacts"); if(res>=0&&res<ll) {SAVE_CONTACTINFO =true; CreateCntctDir=true;}
+ res=LBMSavedData.find("spheres"); if(res>=0&&res<ll) {SAVE_SPHERES =true; CreateDemDir=true;} // To save spheres_* file only if it is required by the operator
+ res=LBMSavedData.find("Spheres"); if(res>=0&&res<ll) {SAVE_SPHERES =true; CreateDemDir=true;} // To save spheres_* file only if it is required by the operator
+
+ // if(NB_DYNGRAINS==1) SAVE_OBSERVEDPTC =true; //Commented to avoid recording of observedPtc if not chosen by the operator
+ if ((ObservedPtc==-1)&&(NB_GRAINS>0)&&(SAVE_OBSERVEDPTC)) ObservedPtc=IdFirstSphere; //Condition If(SAVE_OBSERVEDPTC) is added to save observedPtc only if it is required by the operator
+ if ((SAVE_OBSERVEDPTC)&&((unsigned)ObservedPtc>=LBbodies.size())){cerr<<"Error: ObservedPtc>bodies.size()"<<endl;exit(-1);}
+ if ((SAVE_SPHERES)&&(NB_GRAINS<1)){cerr<<"Error: saving of sphere properties is switched on whereas NB_GRAINS<1"<<endl;exit(-1);}
+
+ createDirectories(CreateLbmDir,CreateDemDir,CreateCntctDir);
+ createNewFiles(); // Files are created now, since we know what the operator wants to record (ObservedPtc ... etc)
+
+
+ /*--------------------------------------------------------------*/
+ /* Periodicity configuration */
+ /*--------------------------------------------------------------*/
+ ll=periodicity.size();
+ res=periodicity.find("x"); if(res>=0&&res<ll) Xperiodicity=true;
+ res=periodicity.find("y"); if(res>=0&&res<ll) Yperiodicity=true;
+ res=periodicity.find("z"); if(res>=0&&res<ll) Zperiodicity=true;
+ cerr <<"Periodicity (XYZ): "<<Xperiodicity<<" "<<Yperiodicity<<" "<<Zperiodicity<<endl;
+
+
+
+ /*--------------------------------------------------------------*/
+ /* Lattice definition */
+ /*--------------------------------------------------------------*/
+ cerr << "---- Lattice setup -----" << endl;
+ State* sWallYm=Body::byId(WallYm_id,scene)->state.get();
+ State* sWallYp=Body::byId(WallYp_id,scene)->state.get();
+ State* sWallXm=Body::byId(WallXm_id,scene)->state.get();
+ State* sWallXp=Body::byId(WallXp_id,scene)->state.get();
+ State* sWallZp=Body::byId(WallZp_id,scene)->state.get();
+ State* sWallZm=Body::byId(WallZm_id,scene)->state.get();
+
+ height = sWallYp->se3.position.y() - sWallYm->se3.position.y();
+ width = sWallXp->se3.position.x() - sWallXm->se3.position.x();
+ depth = sWallZp->se3.position.z() - sWallZm->se3.position.z();
+
+
+ Lx0 = width-Wallthickness;
+ Ly0 = height-Wallthickness;
+ Lz0 = depth-Wallthickness;
+ Real Lx1 = width+Wallthickness;
+ Real Ly1 = height+Wallthickness;
+ Real Lz1 = depth+Wallthickness;
+ /*--------------------------------------------------------------*/
+ /* Computation of number of lattice nodes in each direction */
+ /*--------------------------------------------------------------*/
+ dx = Lx0 / (Real) (Nx-1);
+ invdx=1./dx;
+
+ //Number of nodes (after-correction)
+ Ny = ceil(invdx*Ly1)+1;
+ dx = Ly1 / (Real) (Ny-1);
+ invdx=1./dx;
+ dx2=dx*dx;
+ Nx = ceil(invdx*Lx1)+1;
+ Ny = ceil(invdx*Ly1)+1;
+ Nz=1;
+
+ cerr <<"LXYZ0= "<<Lx0<<" "<<Ly0<<" "<<Lz0<<endl;
+ cerr <<"LXYZ1= "<<Lx1<<" "<<Ly1<<" "<<Lz1<<endl;
+ cerr <<"Ny= "<<Ny<<" "<<" Nx*Ly0/Lx0="<<Nx*Ly0/Lx0<<endl;
+ cerr <<"Wallthickness= "<<Wallthickness<<" "<<" dx="<<dx<<endl;
+
+ /*-------------------------------------------------------*/
+ /* Verification of wall positionning in positive regions */
+ /*-------------------------------------------------------*/
+ if (sWallYm->se3.position.y()- halfWallthickness<0.){cerr <<"Wall position error: please avoid negatives region (Y)"<<endl;exit(-1);}
+ if (sWallXm->se3.position.x()- halfWallthickness<0.){cerr <<"Wall position error: please avoid negatives region (X)"<<endl;exit(-1);}
+ if ((dim==3)&&(sWallZm->se3.position.z()- halfWallthickness<0.)){cerr <<"Wall position error: please avoid negatives region (Z)"<<endl;exit(-1);}
+
+
+ /*--------------------------------------------------------------*/
+ /* Some tests about the validity of parameters */
+ /*--------------------------------------------------------------*/
+ if( (NB_DYNGRAINS==0)&&(use_ConvergenceCriterion)&&(ErrorCriterion==1)){
+ cerr <<"ERROR: can't use ErrorCriterion=1 when (NB_DYNGRAINS=0"<<endl;
+ exit(-1);}
+ if((ObservedNode!=-1)&&(ObservedNode>=Nx*Ny)){
+ cerr <<"Warning ObservedNode is >= Nx*Ny ... exit"<<endl;
+ exit(-1);}
+ if((SaveMode!=1)&&(SaveMode!=2)) {cerr <<"Warning unknown SaveMode."<<SaveMode<<endl; exit(-1);}
+ if((SaveMode==1)&&(IterSave<=0)) {cerr <<"Warning SaveMode==1 and IterSave<=0."<<endl; exit(-1);}
+ if((SaveMode==2)&&(TimeSave<=0)) {cerr <<"Warning SaveMode==2 and TimeSave<=0."<<endl; exit(-1);}
+ /*---------------------------------------------------------------*/
+ /*------------------ general node initialization ----------------*/
+ /*---------------------------------------------------------------*/
+ LBMnode aa;
+ for(int nidx=0; nidx<Nx*Ny; nidx++) {nodes.push_back(aa);}
+ bool j_update=false;
+ int j=0;
+ for (int nidx=0; nidx<Nx*Ny; nidx++){
+ int i=nidx-j*Nx;
+ if((nidx+1)%Nx==0) j_update=true;
+ int k=0;
+
+ nodes[nidx].SetCellIndexesAndPosition(i,j,k);
+ nodes[nidx].DispatchBoundaryConditions(Nx,Ny,Nz);
+ NbNodes++;
+ for (int dndx=0; dndx<NbDir; dndx++){
+ nodes[nidx].links_id.push_back(-1);
+ I=nodes[nidx].i+eib[dndx].x();
+ J=nodes[nidx].j+eib[dndx].y();
+ if(((I==i)&&(J==j)) || (I==-1) || (J==-1) || (I==Nx) || (J==Ny) ){
+ nodes[nidx].neighbour_id.push_back(-1);
+ }
+ else {nodes[nidx].neighbour_id.push_back(I+J*Nx);}
+ }
+ if(j_update) {j++;j_update=false;}
+ }
+
+ ////////////////////////////////////////////////////////////////////////////////////
+ ///FIXME(flomine#1): periodicity should be implemented from links to facilitate the streaming step
+ ///FIXME(flomine#1): to be optimise and bug should be corrected (cf test version)
+ LBMlink bb;
+ int link_id=-1;
+ for (int nidx=0; nidx<Nx*Ny; nidx++){
+ int I=nodes[nidx].i;
+ int J=nodes[nidx].j;
+ for (int dndx=0; dndx<NbDir; dndx++){
+ if(dndx==0) continue;
+ bb.PointingOutside=false;
+ if((!strcmp(model.c_str(), "d2q9" )) && ((dndx==1)||(dndx==2)||(dndx==5)||(dndx==6))){
+ link_id++;bb.i=dndx;bb.nid1=nidx;
+ bb.nid2=nodes[nidx].neighbour_id[dndx];
+ if(bb.nid2==-1) bb.PointingOutside=true;
+ links.push_back(bb);
+ nodes[bb.nid1].links_id[dndx]=link_id;
+ if(bb.nid2!=-1) nodes[bb.nid2].links_id[opp[dndx]]=link_id;
+ }else if(!strcmp(model.c_str(), "d2q9" )){
+ if((I==0)&&(J!=0)&&((dndx==3)||(dndx==7))){
+ link_id++;bb.i=dndx; bb.nid1=nidx;
+ bb.nid2=nodes[nidx].neighbour_id[dndx];
+ bb.PointingOutside=true;
+ if(bb.nid2!=-1) {cerr<<"ERROR: bb.id2!=-1"<<endl;exit(-1);}
+ links.push_back(bb);
+ nodes[bb.nid1].links_id[dndx]=link_id;
+ if(bb.nid2!=-1) nodes[bb.nid2].links_id[opp[dndx]]=link_id;
+ } else if((J==0)&&(I!=0)&&((dndx==4)||(dndx==7)||(dndx==8))){
+ link_id++;bb.i=dndx; bb.nid1=nidx;
+ bb.nid2=nodes[nidx].neighbour_id[dndx];
+ bb.PointingOutside=true;
+ if(bb.nid2!=-1) {cerr<<"ERROR: bb.id2!=-1"<<endl;exit(-1);}
+ links.push_back(bb);
+ nodes[bb.nid1].links_id[dndx]=link_id;
+ if(bb.nid2!=-1) nodes[bb.nid2].links_id[opp[dndx]]=link_id;
+ } else if((I==0)&&(J==0)&&((dndx==3)||(dndx==4)||(dndx==7)||(dndx==8))){
+ link_id++;bb.i=dndx; bb.nid1=nidx;
+ bb.nid2=nodes[nidx].neighbour_id[dndx];
+ bb.PointingOutside=true;
+ if(bb.nid2!=-1) {cerr<<"ERROR: bb.id2!=-1"<<endl;exit(-1);}
+ links.push_back(bb);
+ nodes[bb.nid1].links_id[dndx]=link_id;
+ if(bb.nid2!=-1) nodes[bb.nid2].links_id[opp[dndx]]=link_id;
+ }
+ }else {cerr<<"ERROR: Unknow model type: "<<model<<endl;exit(-1);}
+ }
+ }
+ ////////////////////////////////////////////////////////////////////////////////////
+
+ if((ConvergenceThreshold==-1)||(ConvergenceThreshold==0)) use_ConvergenceCriterion=false;
+ else {use_ConvergenceCriterion=true;ErrorCriterion=ConvergenceThreshold;}
+ /*------------------------------------------------------------------*/
+ /*------------------------ LBM PARAMETERS --------------------------*/
+ /*------------------------------------------------------------------*/
+ UMaxtheo = dP.norm()/(8.*Rho*Nu)*(height*height/width);
+ if(defaultLbmInitMode){
+ uMax = 0.1; c = UMaxtheo/uMax;
+ cs = c/sqrt(3.); dt = dx/c;
+ nu = dt/(dx*dx/Nu); tau = 3.*nu + 0.5;
+ }else{
+ nu = (1.0/3.0)*(tau - 0.5);
+ dt = nu*(dx*dx/Nu); c = dx/dt;
+ cs = c/sqrt(3.); uMax = UMaxtheo / c;
+ }
+ omega = 1.0/tau;
+ c2=c*c;
+ invdt=1./dt;
+ MaxBodyRadius=invdx*MaxBodyRadius;
+ MinBodyRadius=invdx*MinBodyRadius;
+ MeanBodyRadius=invdx*MeanBodyRadius;
+ outside_limit=1.5*width;
+ }//end if FirstRun
+
+ /*************************************************************************/
+ /* SOLID OBSTACLES SET-UP */
+ /*************************************************************************/
+ Real Rmin = 1000.;
+ int newFluidCells_couter=0;
+ int newObstacleCells_couter=0;
+
+ State* sWallYm=Body::byId(WallYm_id,scene)->state.get();
+ State* sWallYp=Body::byId(WallYp_id,scene)->state.get();
+ State* sWallXm=Body::byId(WallXm_id,scene)->state.get();
+ State* sWallXp=Body::byId(WallXp_id,scene)->state.get();
+ State* sWallZp=Body::byId(WallZp_id,scene)->state.get();
+ State* sWallZm=Body::byId(WallZm_id,scene)->state.get();
+
+ timingDeltas->checkpoint("Reinit:Nodes0");
+ #pragma omp parallel for
+ for (int nidx=0; nidx<Nx*Ny; nidx++){
+ /*------------------------------------------*/
+ /* Reinitialization: */
+ /*------------------------------------------*/
+ nodes[nidx].body_id=-1;
+ nodes[nidx].setAsFluid();
+ nodes[nidx].isObstacleBoundary=false;
+ nodes[nidx].isFluidBoundary=false;
+ nodes[nidx].isNewObstacle=false;
+ nodes[nidx].isNewFluid=false;
+
+ /*--- according to X+ ---*/
+ if (useWallXp&&(nodes[nidx].i>=invdx*(sWallXp->se3.position.x() - halfWallthickness))){
+ nodes[nidx].setAsObstacle();
+ nodes[nidx].isObstacleBoundary=true;
+ nodes[nidx].body_id=WallXp_id;
+ NbSolidNodes++;}
+ /*--- according to X- ---*/
+ else if (useWallXm&&(nodes[nidx].i<=invdx*(sWallXm->se3.position.x() + halfWallthickness))){
+ nodes[nidx].setAsObstacle();
+ nodes[nidx].isObstacleBoundary=true;
+ nodes[nidx].body_id=WallXm_id;
+ NbSolidNodes++;}
+ /*--- according to Y+ ---*/
+ else if (useWallYp&&(nodes[nidx].j>=invdx*(sWallYp->se3.position.y() - halfWallthickness))){
+ nodes[nidx].setAsObstacle();
+ nodes[nidx].isObstacleBoundary=true;
+ nodes[nidx].body_id=WallYp_id;
+ NbSolidNodes++;}
+ /*--- according to Y- ---*/
+ else if (useWallYm&&(nodes[nidx].j<=invdx*(sWallYm->se3.position.y() + halfWallthickness))){
+ nodes[nidx].setAsObstacle();
+ nodes[nidx].isObstacleBoundary=true;
+ nodes[nidx].body_id=WallYm_id;
+ NbSolidNodes++;}
+ /*--- according to Z+ ---*/
+ else if (useWallZp&&(nodes[nidx].k>=invdx*(sWallZp->se3.position.z() - halfWallthickness))){
+ nodes[nidx].setAsObstacle();
+ nodes[nidx].isObstacleBoundary=true;
+ nodes[nidx].body_id=WallZp_id;
+ NbSolidNodes++;}
+ /*--- according to Z- ---*/
+ else if (useWallZm&&(nodes[nidx].k<=invdx*(sWallZm->se3.position.z() + halfWallthickness))){
+ nodes[nidx].setAsObstacle();
+ nodes[nidx].isObstacleBoundary=true;
+ nodes[nidx].body_id=WallZm_id;
+ NbSolidNodes++;}
+
+ if(firstRun){nodes[nidx].wasObstacle=nodes[nidx].isObstacle;}
+ }
+
+ /*---------------------------------------------------------------*/
+ /*- Solid particle detection and recording of their properties --*/
+ /*---------------------------------------------------------------*/
+ NumberOfDynamicParticles=0;
+ if(removingCriterion!=0) IdOfNextErodedPtc.clear();
+ FOREACH(const shared_ptr<Body>& b, *scene->bodies){
+ if(!b) continue; // deleted bodies
+ State* state=b->state.get();
+ const int id=b->getId();
+ //if ((b->shape->getClassName()=="Sphere")&&(b->isDynamic())){ //ModLuc: removing of b->isDynamic() in order that non dynamic particle can be seen by the LBM
+ if (b->shape->getClassName()=="Sphere"){
+ const shared_ptr<Sphere>& sphere = YADE_PTR_CAST<Sphere> ( b->shape );
+
+ LBbodies[id].pos=invdx*state->pos;
+ LBbodies[id].vel=(state->vel)/c;
+ LBbodies[id].AVel= (state->angVel)*dt;
+ LBbodies[id].radius=invdx*RadFactor*(sphere->radius);
+
+ CurMinVelOfPtc=min(CurMinVelOfPtc,state->vel.norm());
+ CurMaxVelOfPtc=max(CurMaxVelOfPtc,state->vel.norm());
+
+ Vector3r posMax=LBbodies[id].pos+ Vector3r(LBbodies[id].radius,LBbodies[id].radius,LBbodies[id].radius);
+ Vector3r posMin=LBbodies[id].pos- Vector3r(LBbodies[id].radius,LBbodies[id].radius,LBbodies[id].radius);
+
+ Vector3r dist=Vector3r::Zero();
+ for(int ii=posMin[0]-1;ii<=posMax[0]+1;ii++)
+ for(int jj=posMin[1]-1;jj<=posMax[1]+1;jj++){
+ if((ii==-1)||(ii==Nx)||(jj==-1)||(jj==Ny)) continue;
+ if((ii<-1)||(ii>Nx)||(jj<-1)||(jj>Ny)) continue;
+ if (LBbodies[id].radius < Rmin) Rmin = LBbodies[id].radius;
+ int nidx=ii+jj*Nx;
+ dist=nodes[nidx].posb-LBbodies[id].pos;
+ if(dist.norm()<LBbodies[id].radius){
+ nodes[nidx].body_id = id;
+ nodes[nidx].setAsObstacle();
+ NbSolidNodes++;
+ NbParticleNodes++;}
+ if(firstRun){nodes[nidx].wasObstacle=nodes[nidx].isObstacle;}
+ }
+ /*-------------------------------------------------------------------*/
+ /* ///NOTE : this should be removed since it can be done with python */
+ ///Fck: pas en MODE 1
+ /*-------------------------------------------------------------------*/
+ //cerr <<"Check removing"<<endl;
+ //if(removingCriterion!=0){ //ModLuc adding of b->isDynamic() to remove only dynamic particles and not already removed particles changed into non dynamic
+ if((removingCriterion!=0)&&(b->isDynamic())){
+ switch(removingCriterion){
+ case 1:
+ /* criterion with respect to the particle postion in x direction */
+ if(LBbodies[id].pos.x()>(invdx*(sWallXp->se3.position.x())-1.05*MaxBodyRadius/RadFactor)){
+ IdOfNextErodedPtc.push_back(id);
+ }
+ break;
+ case 2:
+ /* criterion on particle velocity */
+ if((LBbodies[id].vel.norm()>VelocityThreshold)||(LBbodies[id].pos.x()>(invdx*(sWallXp->se3.position.x())-2.*LBbodies[id].radius))) {IdOfNextErodedPtc.push_back(id);}
+ break;
+ default:
+ exit(-1);
+ break;
+ }
+ }
+ //NumberOfDynamicParticles++; //ModLuc: to still count only dynamic particles and not all particles
+ if(b->isDynamic()) NumberOfDynamicParticles++;
+ }
+
+ LBbodies[id].force=Vector3r::Zero();
+ LBbodies[id].momentum=Vector3r::Zero();
+ }
+
+
+ /*------------------------------------------------------------------*/
+ /*------------------ detection of boundary nodes -------------------*/
+ /*------------------------------------------------------------------*/
+ #pragma omp parallel for
+ for (int nidx=0; nidx<Nx*Ny; nidx++)
+ if(nodes[nidx].isObstacle){
+ for(unsigned int n=0;n<nodes[nidx].neighbour_id.size();n++){
+ if(nodes[nidx].neighbour_id[n]!=-1){
+ int nidx2=nodes[nidx].neighbour_id[n];
+ if(nodes[nidx].isObstacle!=nodes[nidx2].isObstacle) {
+ nodes[nidx].isObstacleBoundary=true;
+ nodes[nidx2].isFluidBoundary=true;
+ int BodyId=nodes[nidx].body_id;
+ int lid=nodes[nidx].links_id[n];
+ links[lid].isBd=true;
+ links[lid].sid=nidx;
+ links[lid].fid=nidx2;
+ links[lid].idx_sigma_i=opp[n];
+ if(LBbodies[BodyId].isPtc()){
+ links[lid].DistMid= nodes[nidx].posb-0.5*eib[links[lid].idx_sigma_i]-LBbodies[BodyId].pos;
+ links[lid].VbMid=LBbodies[BodyId].vel+LBbodies[BodyId].AVel.cross(links[lid].DistMid);
+ if(links[lid].VbMid.norm()<VbCutOff) links[lid].VbMid=Vector3r::Zero();
+ }
+ if(LBbodies[BodyId].isBox()){
+ links[lid].DistMid= Vector3r::Zero();
+ links[lid].VbMid=Vector3r::Zero();
+ }
+ }
+ }
+ }
+ }
+ #pragma omp parallel for
+ for (int nidx=0; nidx<Nx*Ny; nidx++)
+ if((nodes[nidx].isObstacle)&&(!nodes[nidx].isObstacleBoundary)){
+ nodes[nidx].setAsFluid();
+ NbSolidNodes--;
+ if(firstRun) nodes[nidx].wasObstacle=nodes[nidx].isObstacle;
+ }
+
+ NbFluidNodes=NbNodes-NbSolidNodes;
+ /*----------------------------------------------------------------------*/
+
+ if(firstRun){
+
+ /*------------------------------------------------------------------*/
+ /*------------------------------- MODE -----------------------------*/
+ /*------------------------------------------------------------------*/
+ if((IterSubCyclingStart<=0) && (IterMax>1) ) MODE = 1;
+ if((IterSubCyclingStart<=0) && (IterMax==1) ) MODE = 2;
+ if(IterSubCyclingStart>0) MODE = 3;
+
+ /*------------------------------------------------------------------*/
+ /*---------------------------- SUBCYCLING --------------------------*/
+ /*------------------------------------------------------------------*/
+ if(MODE==3&&IterMax<IterSubCyclingStart){cerr <<"Exit because Itermax<IterSubCyclingStart"<<endl;exit(-1);}
+ if(DemIterLbmIterRatio==-1){
+ if(MODE==1) DemIterLbmIterRatio=1;
+ else DemIterLbmIterRatio=int(dt/DEMdt0);
+ newDEMdt=dt/DemIterLbmIterRatio;
+ }else{newDEMdt=dt/DemIterLbmIterRatio;}
+ scene->dt=newDEMdt;
+ if(SaveMode==2){
+ IterSave=TimeSave/(dt);
+ if(TimeSave<dt) {cerr <<"Warning SaveMode==2 and TimeSave<dt."<<endl; exit(-1);}
+ }
+ writelogfile();
+
+ /*------------------------------------------*/
+ /* Initialization of distribution functions */
+ /*------------------------------------------*/
+ for (int nidx=0; nidx<Nx*Ny; nidx++){
+ nodes[nidx].rhob=1.;
+ nodes[nidx].velb=Vector3r::Zero();
+ for (int didx=0; didx<NbDir; didx++){
+ cub = 3.0* eib[didx].dot(nodes[nidx].velb);
+ feqb = w[didx]*nodes[nidx].rhob*(1.0 + cub + 0.5*(cub*cub) - 1.5*((nodes[nidx].velb.x()*nodes[nidx].velb.x()) + (nodes[nidx].velb.y()*nodes[nidx].velb.y())));
+ nodes[nidx].f.push_back(feqb);
+ nodes[nidx].fpostcol.push_back(0.);
+ nodes[nidx].fprecol.push_back(0.);
+ }
+ }
+ firstRun = false;
+ }
+
+ /*----------------------------------------------------------------------*/
+ /* FOR ALL ITERATIONS */
+ /*----------------------------------------------------------------------*/
+ Real Error = 1.0;
+ DEMdt=scene->dt;
+ DEM_TIME=DEM_ITER*DEMdt;
+ for (iter=0; iter<IterMax; iter++){
+ step += 1;
+ LBM_ITER++;
+ LBM_TIME=LBM_ITER*dt;
+
+ /*------------------------------------------------------------------*/
+ /* REINITIALIZATION RELATIVE TO THE MODE 1 */
+ /*------------------------------------------------------------------*/
+ if(MODE==1){
+ FOREACH(const shared_ptr<Body>& b, *scene->bodies){
+ if(!b) continue;
+ const int id=b->getId();
+ LBbodies[id].force=Vector3r::Zero();
+ LBbodies[id].momentum=Vector3r::Zero();}
+ }
+
+ /*------------------------------------------------------------------*/
+ /* GENERAL REINITIALIZATION */
+ /*------------------------------------------------------------------*/
+ Vector3r WallBottomVel=Vector3r::Zero();//m s-1
+ FmoyCur=0.;
+ VmeanFluidC=0.; VmaxC=-1000000.; VminC=1000000.;
+ RhomaxC=-1000000.; RhominC=1000000.;RhoTot=0.;
+ /*------------------------------------------------------------------*/
+ /* Loop on nodes */
+ /*------------------------------------------------------------------*/
+ #pragma omp parallel for
+ for (int nidx=0; nidx<Nx*Ny; nidx++){
+ if(nodes[nidx].checkIsNewObstacle()) {newObstacleCells_couter++;}
+ else{if(nodes[nidx].checkIsNewFluid()) {newFluidCells_couter++;}}
+ if(nodes[nidx].applyBC){
+ Vector3r U=Vector3r::Zero();
+ Real density=0.;
+ /*----------- inlet ------------*/
+ if(nodes[nidx].applyXmBC){
+ if(XmBCType==1){
+ density=1.0 + dP.x()/(Rho*cs*cs);
+ U=Vector3r(1.0-((nodes[nidx].f[0]+nodes[nidx].f[2]+nodes[nidx].f[4]) + 2.0*(nodes[nidx].f[3]+nodes[nidx].f[6]+nodes[nidx].f[7]))/density,0.,0.);
+ }else if(XmBCType==2){
+ U=Vector3r::Zero();
+ density=(nodes[nidx].f[0]+nodes[nidx].f[2]+nodes[nidx].f[4]+2.*(nodes[nidx].f[3]+nodes[nidx].f[6]+nodes[nidx].f[7]))/(1.-U.x());
+ }
+ nodes[nidx].MixteBC(model,density,U,"Xm");
+ /*----------- outlet ------------*/
+ }else if( nodes[nidx].applyXpBC){
+ if(XpBCType==1){
+ density=1.0;
+ U=Vector3r(-1.0 + ((nodes[nidx].f[0]+nodes[nidx].f[2]+nodes[nidx].f[4]) + 2.0*(nodes[nidx].f[1]+nodes[nidx].f[5]+nodes[nidx].f[8]))/density,0.,0.);
+ }else if(XpBCType==2){
+ U=Vector3r::Zero();
+ density=(nodes[nidx].f[0]+nodes[nidx].f[2]+nodes[nidx].f[4]+2.*(nodes[nidx].f[1]+nodes[nidx].f[5]+nodes[nidx].f[8]))/(1.+U.x());
+ }
+ nodes[nidx].MixteBC(model,density,U,"Xp");
+ /*----------- top ------------*/
+ } else if( nodes[nidx].applyYpBC){
+ if(YpBCType==1){
+ density=1.0;
+ U=Vector3r(0.,-1.0 + ((nodes[nidx].f[0]+nodes[nidx].f[1]+nodes[nidx].f[3]) + 2.0*(nodes[nidx].f[2]+nodes[nidx].f[5]+nodes[nidx].f[6]))/density,0.);
+ }else if(YpBCType==2){
+ U=Vector3r::Zero();
+ density=(nodes[nidx].f[0]+nodes[nidx].f[1]+nodes[nidx].f[3]+2.*(nodes[nidx].f[2]+nodes[nidx].f[5]+nodes[nidx].f[6]))/(1.+U.y());
+ }
+ nodes[nidx].MixteBC(model,density,U,"Yp");
+ /*----------- bottom ------------*/
+ }else if( nodes[nidx].applyYmBC){
+ if(YmBCType==1){
+ density=1.0;
+ U=Vector3r(0.,1.0-((nodes[nidx].f[0]+nodes[nidx].f[1]+nodes[nidx].f[3]) + 2.0*(nodes[nidx].f[4]+nodes[nidx].f[7]+nodes[nidx].f[8]))/density,0.);
+ }else if(YmBCType==2){
+ U=Vector3r::Zero();
+ density=(nodes[nidx].f[0]+nodes[nidx].f[1]+nodes[nidx].f[3]+2.*(nodes[nidx].f[4]+nodes[nidx].f[7]+nodes[nidx].f[8]))/(1.-U.y());
+ }
+ nodes[nidx].MixteBC(model,density,U,"Ym");
+ /*----------- bottom-left ------------*/
+ }else if(nodes[nidx].applyYmXmBC){
+ if(XmYmZpBCType==1){
+ cerr <<"XmYmZpType=1 not implemented . Exit"<<endl;
+ exit(-1);
+ }else if(XmYmZpBCType==2){
+ U=Vector3r::Zero();
+ density=nodes[nidx+1].rhob;
+ }
+ nodes[nidx].MixteBC(model,density,U,"XmYmZp");
+ /*----------- top-left ------------*/
+ }else if(nodes[nidx].applyYpXmBC){
+ if(XmYpZpBCType==1){
+ cerr <<"XmYpZpBCType=1 not implemented . Exit"<<endl;
+ exit(-1);
+ }else if(XmYpZpBCType==2){
+ U=Vector3r::Zero();
+ density=nodes[nidx+1].rhob;
+ }
+ nodes[nidx].MixteBC(model,density,U,"XmYpZp");
+ /*----------- bottom-right ------------*/
+ }else if(nodes[nidx].applyYmXpBC){
+ if(XpYmZpBCType==1){
+ cerr <<"XpYmZpBCType=1 not implemented . Exit"<<endl;
+ exit(-1);
+ }else if(XpYmZpBCType==2){
+ U=Vector3r::Zero();
+ density=nodes[nidx-1].rhob;
+ }
+ nodes[nidx].MixteBC(model,density,U,"XpYmZp");
+ /*----------- top-right ------------*/
+ }else if(nodes[nidx].applyYpXpBC){
+ if(XpYpZpBCType==1){
+ cerr <<"XpYpZpBCType=1 not implemented . Exit"<<endl;
+ exit(-1);
+ }else if(XpYpZpBCType==2){
+ U=Vector3r::Zero();
+ density=nodes[nidx-1].rhob;
+ }
+ nodes[nidx].MixteBC(model,density,U,"XpYpZp");
+ }else{
+ cerr << "ERROR: node "<<nidx<<". Looking for a BC to apply ..."<<endl;
+ exit(-1);
+ }
+ }
+
+ nodes[nidx].rhob=0.;
+ nodes[nidx].velb=Vector3r::Zero();
+ nodes[nidx].IsolNb=0;
+ if(nodes[nidx].isFluidBoundary){nodes[nidx].IsolNb=8;}
+
+ for (int dndx=0; dndx<NbDir; dndx++){
+ nodes[nidx].fprecol[dndx] = nodes[nidx].f[dndx];
+ nodes[nidx].velb += eib[dndx]*nodes[nidx].f[dndx];
+ nodes[nidx].rhob += nodes[nidx].f[dndx];
+ if((nodes[nidx].isFluidBoundary)&&(nodes[nidx].neighbour_id[dndx]!=-1)){
+ int ns=nodes[nidx].neighbour_id[dndx];
+ if(!nodes[ns].isObstacle) nodes[nidx].IsolNb=nodes[nidx].IsolNb-1;
+ if(nodes[nidx].IsolNb<0) {cerr<<"isolNb<0"<<endl;exit(-1);}}
+
+ }
+ nodes[nidx].velb /= nodes[nidx].rhob;
+
+ Real temp0=1.5*((nodes[nidx].velb.x()*nodes[nidx].velb.x())+(nodes[nidx].velb.y()*nodes[nidx].velb.y()));
+ Real cub0 = 3.0* eib[0].dot(nodes[nidx].velb);
+ nodes[nidx].fpostcol[0]= nodes[nidx].f[0] - omega * (nodes[nidx].f[0]-(nodes[nidx].rhob* w[0]*( 1. + cub0 + 0.5*(cub0*cub0) - temp0)));
+
+
+ nodes[nidx].fpostcol[0] = nodes[nidx].fpostcol[0] + (nodes[nidx].rhob* w[0])/c2 * eib[0].dot(CstBodyForce);
+ nodes[nidx].f[0]=nodes[nidx].fpostcol[0];
+ RhoTot+=nodes[nidx].rhob;
+ if(nodes[nidx].body_id==-1)if(VmaxC<c*nodes[nidx].velb.norm()) VmaxC=c*nodes[nidx].velb.norm();
+ if(VminC>c*nodes[nidx].velb.norm()) VminC=c*nodes[nidx].velb.norm();
+ if(RhomaxC<Rho*nodes[nidx].rhob) RhomaxC=Rho*nodes[nidx].rhob;
+ if(RhominC>Rho*nodes[nidx].rhob) RhominC=Rho*nodes[nidx].rhob;
+ if(!nodes[nidx].isObstacle) VmeanFluidC+=c*nodes[nidx].velb.norm();
+ }
+
+ #pragma omp parallel for
+ for(unsigned int lid=0;lid<links.size();lid++){
+ int nidx1 = links[lid].nid1;
+ int nidx2 = links[lid].nid2;
+ int dndx1 = links[lid].i;
+ int dndx2 = opp[links[lid].i];
+
+ /*-------------------------------------------------- ---------------*/
+ /* equilibrium functions and collisions */
+ /*------------------------------------------------------------------*/
+ Real temp1=1.5*((nodes[nidx1].velb.x()*nodes[nidx1].velb.x())+(nodes[nidx1].velb.y()*nodes[nidx1].velb.y()));
+ Real cub1 = 3.0* eib[dndx1].dot(nodes[nidx1].velb);
+ nodes[nidx1].fpostcol[dndx1] = nodes[nidx1].fprecol[dndx1] - omega * (nodes[nidx1].fprecol[dndx1]-(nodes[nidx1].rhob* w[dndx1]*( 1. + cub1 + 0.5*(cub1*cub1) - temp1)));
+ nodes[nidx1].fpostcol[dndx1] = nodes[nidx1].fpostcol[dndx1] + (nodes[nidx1].rhob* w[dndx1])/c2 * eib[dndx1].dot(CstBodyForce);
+ if(!links[lid].PointingOutside){
+ Real temp2=1.5*((nodes[nidx2].velb.x()*nodes[nidx2].velb.x())+(nodes[nidx2].velb.y()*nodes[nidx2].velb.y()));
+ Real cub2 = 3.0* eib[dndx2].dot(nodes[nidx2].velb);
+ nodes[nidx2].fpostcol[dndx2] = nodes[nidx2].fprecol[dndx2] - omega * (nodes[nidx2].fprecol[dndx2]-(nodes[nidx2].rhob* w[dndx2]*( 1. + cub2 + 0.5*(cub2*cub2) - temp2)));
+ nodes[nidx2].fpostcol[dndx2] = nodes[nidx2].fpostcol[dndx2] + (nodes[nidx2].rhob* w[dndx2])/c2 * eib[dndx2].dot(CstBodyForce);
+ }
+
+ /*-------------------------------------------------- ---------------*/
+ /* Streaming */
+ /*------------------------------------------------------------------*/
+ if(links[lid].PointingOutside){
+ /// Periodicity is currently disabled until it is implemented through the link list.
+///FIXME
+ I=nodes[nidx1].i+eib[dndx1].x();
+ J=nodes[nidx1].j+eib[dndx1].y();
+ if(Xperiodicity){ if (I==Nx) {I=0;} else {if (I==-1) { I=Nx-1;}} }
+ if(Yperiodicity){ if (J==Ny) {J=0;} else {if (J==-1) { J=Ny-1;}} }
+ }else{
+ nodes[nidx1].f[dndx2]=nodes[nidx2].fpostcol[dndx2];
+ nodes[nidx2].f[dndx1]=nodes[nidx1].fpostcol[dndx1];
+ }
+
+ if(links[lid].isBd==false) continue;
+
+ int idx_sigma_i=links[lid].idx_sigma_i;
+ int sid= links[lid].sid;
+ int fid= links[lid].fid;
+ int BodyId=nodes[sid].body_id;
+
+
+ /*--- forces and momenta for this boundary link ---*/
+ links[lid].ct=3.0*w[idx_sigma_i]*nodes[sid].rhob*eib[links[lid].idx_sigma_i].dot(links[lid].VbMid);
+ Vector3r force_ij = eib[links[lid].idx_sigma_i] * (nodes[fid].fpostcol[idx_sigma_i] - links[lid].ct);
+ Vector3r lubforce_ij = Vector3r::Zero();
+ Vector3r totalforce_ij = force_ij+lubforce_ij;
+ Vector3r totalmomentum_ij = links[lid].DistMid.cross(totalforce_ij);
+
+ /* Sum over all boundary links of all boundary nodes */
+ LBbodies[BodyId].force=LBbodies[BodyId].force+totalforce_ij;
+ LBbodies[BodyId].momentum=LBbodies[BodyId].momentum+totalmomentum_ij;
+
+ /*------------------------------------------------------*/
+ /* Modified Bounce back rule */
+ if(nodes[fid].IsolNb>=5) {links[lid].VbMid=Vector3r::Zero();links[lid].ct=0.;}
+ nodes[fid].f[opp[idx_sigma_i]] = nodes[fid].fpostcol[idx_sigma_i] - 2.0*links[lid].ct;
+ nodes[sid].f[idx_sigma_i] = nodes[sid].fpostcol[opp[idx_sigma_i]]+ 2.0*links[lid].ct;
+ if( (MODE==2)||((MODE==3)&&(IterMax==1)) ) {links[lid].ReinitDynamicalProperties();}
+
+ }
+ VmeanFluidC=VmeanFluidC/NbFluidNodes;
+
+
+ /*---------------------------------------------*/
+ /* Stop criteria */
+ /*---------------------------------------------*/
+// if(use_ConvergenceCriterion){
+// switch(ErrorCriterion){
+// case 1:
+// /*--------------------------------------------------------*/
+// /* Criterion based on the mean force */
+// /*--------------------------------------------------------*/
+// if((LBM_ITER > 1000) & (LBM_ITER % 10 == 0)){
+// for (int s=NB_WALLS ; s<NB_BODIES; s++) {FmoyCur = FmoyCur + LBbodies[s].force.norm();}
+// FmoyCur = FmoyCur/(NB_DYNGRAINS);
+// if (FmoyCur!=0.){
+// Real ErrorA = abs(FmoyCur-FmoyPrev)/abs(FmoyCur);
+// Real ErrorB = abs(FmoyCur-FmoyPrevPrev)/abs(FmoyCur);
+// Error=max(ErrorA,ErrorB);
+// FmoyPrevPrev=FmoyPrev;
+// FmoyPrev=FmoyCur;
+// }
+// }
+// break;
+// case 2:
+// /*--------------------------------------------------------*/
+// /* Criterion based on the mean velocity */
+// /*--------------------------------------------------------*/
+// if((LBM_ITER > 100) & (LBM_ITER % 10 == 0)){
+// if (VmeanFluidC!=0.){
+// Real ErrorA = abs(VmeanFluidC-PrevVmeanFluidC)/abs(VmeanFluidC);
+// //Real ErrorB = abs(VmeanFluidC-PrevPrevVmeanFluidC)/abs(VmeanFluidC);
+// //Error=max(ErrorA,ErrorB);
+// Error= ErrorA;
+// PrevPrevVmeanFluidC=PrevVmeanFluidC;
+// PrevVmeanFluidC=VmeanFluidC;
+// }
+// }
+// break;
+// default:
+// cerr <<"Unknow ErrorCriterion value ! "<<endl;
+// exit(-1);
+// break;
+// }
+// }
+
+ if(LBM_ITER%IterPrint==0){
+ cerr.precision(6);
+ cerr <<"__________________________________________________________________________"<<endl;
+ cerr << "| Run in mode : "<<MODE<<endl;
+ cerr <<"| New Obstacle nodes: "<<newObstacleCells_couter<<" New Fluid nodes: "<<newFluidCells_couter<<endl;
+ cout <<"| height/width/depth (m) "<<height<<" / "<<width<<" / "<<depth<<endl;
+ cerr <<"|------------------------------------------------------------------------|"<<endl;
+ cerr <<"| \t\tDEM\t\t | \t\tLBM\t\t |"<<endl;
+ cerr <<"|------------------------------------------------------------------------|"<<endl;
+ cerr <<"| t (s)\t: "<<DEM_TIME<<"\t\t | t (s)\t\t: "<<LBM_TIME<<"\t\t "<<endl;
+ cerr <<"| Iter\t\t: "<<DEM_ITER<<"\t\t | Iter\t\t: "<<LBM_ITER<<"\t\t "<<endl;
+ cerr <<"| Nb dyn ptc\t: "<<NumberOfDynamicParticles<<"\t\t | M\t\t: "<<VmaxC/c<<"\t "<<endl;
+ cerr <<"| \t\t\t\t | "<<VminC<<"\t< <V(t)> (m/s)= "<<VmeanFluidC <<"\t< "<<VmaxC<<"\t "<<endl;
+ cerr <<"| \t\t\t\t | "<<RhominC<<"\t< rho(t) (m3/kg) <"<<RhomaxC<<"\t "<<endl;
+ cerr <<"| \t\t\t\t | RhoTot (adim) \t: "<<RhoTot<<"\t "<<endl;
+ if(ObservedPtc!=-1){
+ Vector3r tmp=2.*Rho*c2*dx*LBbodies[ObservedPtc].force; cerr <<"| \t\t\t\t\t | F (N) \t: "<<tmp<<"\t "<<endl;
+ tmp=LBbodies[ObservedPtc].pos*dx; cerr <<"| \t\t\t\t\t | pos (m) \t: "<<tmp<<"\t "<<endl;
+ tmp=LBbodies[ObservedPtc].vel*c; cerr <<"| \t\t\t\t\t | VPtc (m/s)\t: "<<tmp<<"\t \t\t\t\t "<<endl;
+ }
+ if(((MODE==2)||(MODE==3))&&(NB_DYNGRAINS>0)){
+ cerr <<"| VminPtcC (m/s)\t: "<<CurMinVelOfPtc<<"\t |\t\t\t\t "<<endl;
+ cerr <<"| VmaxPtcC (m/s)\t: "<<CurMaxVelOfPtc<<"\t |\t\t\t\t "<<endl;
+ }
+ if( ((MODE==1)&&firstRun) || (MODE==2) || (MODE==3)){
+ cerr <<"| \t\t\t\t\t | NbNodes\t: "<<NbNodes<<"\t "<<endl;
+ cerr <<"| \t\t\t\t\t | Fluid nodes\t: "<<NbFluidNodes<<"\t\t "<<endl;
+ cerr <<"| \t\t\t\t\t | Solid nodes\t: "<<NbSolidNodes<<"\t\t "<<endl;
+ }
+ if(ObservedNode!=-1){
+ cerr <<"| \t\t\t\t\t | Track node\t: "<<ObservedNode<<"\t\t "<<endl;
+ cerr <<"| \t\t\t\t\t | \t V(m/s)\t: "<<nodes[ObservedNode].velb.norm()<<"\t\t "<<endl;
+ cerr <<"| \t\t\t\t\t | \t rho\t: "<<Rho*nodes[ObservedNode].rhob<<"\t\t "<<endl;
+ }
+ if(use_ConvergenceCriterion){
+ cerr <<"| \t\t\t\t\t | Error\t: "<<Error<<"\t\t "<<endl;
+ cerr <<"| \t\t\t\t\t | with criterion "<<ErrorCriterion<<" \t\t "<<endl;
+ }
+ cerr <<"|------------------------------------------------------------------------|"<<endl;
+ }
+
+
+ /*-------------------------------------------------------------------------------*/
+ /*------------- RESULT RECORDING DURING COMPUTATION (MODE 1)--------------------*/
+ /*-------------------------------------------------------------------------------*/
+ if(((iter % (IterSave*DemIterLbmIterRatio) == 0)||(firstRun))&&MODE==1) {
+ if(((iter % (IterSave*SaveGridRatio*DemIterLbmIterRatio) == 0)||(firstRun))&&MODE==1) save(iter,dt);
+ saveStats(iter,dt);
+ if(SAVE_OBSERVEDPTC) {
+ CalculateAndApplyForcesAndTorquesOnBodies(true,false);
+ saveObservedPtc(iter,dt);
+ }
+ if(SAVE_OBSERVEDNODE) saveObservedNode(iter,dt);
+ }
+
+
+ if((IterSubCyclingStart>0)&&(iter+1>=IterSubCyclingStart)) {modeTransition();}
+ if ((Error < ConvergenceThreshold)&&(use_ConvergenceCriterion)) {
+ if(MODE==1) break;
+ if(MODE==3) modeTransition();
+ }
+ if((EndTime>0)&&(LBM_TIME>EndTime)) LbmEnd();
+ //if((InitialNumberOfDynamicParticles!=0)&&(NumberPtcEroded==InitialNumberOfDynamicParticles)) LbmEnd();
+}
+/*------------------- End of the LBM loop (iterations) --------------------*/
+
+
+if(MODE==1) {
+ cerr << "LBM ended after " << step << " iterations";
+ cerr <<" | LBM_ITER = "<< LBM_ITER<<" | Error = " << Error << endl;
+}
+
+/*--------------------------------------------------------------------------------*/
+/*-------------- APPLICATION OF HYDRODYNAMIC FORCES ON SPHERES -------------------*/
+/*--------------------------------------------------------------------------------*/
+if(applyForcesAndTorques) CalculateAndApplyForcesAndTorquesOnBodies(true,true);
+
+/*----------------------------------------------------------------------------*/
+/*----------------- SPHERES ARE MOVED OUTSIDE THE SIMULATION DOMAIN --------- */
+/*----------------------------------------------------------------------------*/
+//if(removingCriterion!=0){
+// for(unsigned int aa=0; aa<IdOfNextErodedPtc.size();aa++){
+// int IdRemoved=IdOfNextErodedPtc[aa];
+// LBbodies[IdRemoved].isEroded=true;
+// shared_ptr<Body> b=Body::byId(IdRemoved);
+// const shared_ptr<Sphere>& sphere = YADE_PTR_CAST<Sphere> ( b->shape );
+// Real r=sphere->radius;
+// outside_limit=outside_limit+1.1*r;
+// b->state->pos=Vector3r(outside_limit,0.,0.);
+// b->setDynamic(false);
+// NB_DYNGRAINS--;
+// NB_DYNBODIES--;
+// outside_limit=outside_limit+1.1*r;
+// cerr <<"Eroded Ptc: "<<IdRemoved<<endl;
+// NumberPtcEroded+=1;
+// if(!strcmp(model.c_str(), "d2q9" )) Vr += Mathr::PI*(r*r);
+// else Vr += 4./3.*Mathr::PI*(r*r*r);
+// }
+// if(IdOfNextErodedPtc.size()!=0) saveEroded(DEM_ITER,DEMdt);
+//}
+
+/*----------------------------------------------------------------------------*/
+/*------------- RESULT RECORDING AT THE END OF THE COMPUTATION ---------------*/
+/*----------------------------------------------------------------------------*/
+ if(((DEM_ITER % (IterSave*DemIterLbmIterRatio) == 0)||(firstRun))&&IterMax==1 ){
+ if(((DEM_ITER % (IterSave*SaveGridRatio*DemIterLbmIterRatio) == 0)||(firstRun))&&IterMax==1 ) save(DEM_ITER,DEMdt);
+ saveStats(DEM_ITER,DEMdt);
+ //saveEroded(DEM_ITER,DEMdt);
+ if(SAVE_OBSERVEDPTC) saveObservedPtc(DEM_ITER,DEMdt);
+ if(SAVE_OBSERVEDNODE) saveObservedNode(DEM_ITER,DEMdt);
+ }
+}
+
+void HydrodynamicsLawLBM::save(int iter_number, Real timestep)
+{
+
+ /*--------------------------------------------*/
+ /* Recording of properties of LBM nodes */
+ /*--------------------------------------------*/
+
+ bool saveLBMdata = false; // Flag to avoid to loop over all the lbm node if no lbm data are saved but only spheres data
+
+ std::stringstream Vfile_name;
+ std::stringstream Vxfile_name;
+ std::stringstream Vyfile_name;
+ std::ofstream Vfile;
+ std::ofstream Vxfile;
+ std::ofstream Vyfile;
+
+ std::stringstream FxLBMfile_name;
+ std::stringstream FyLBMfile_name;
+ std::stringstream MzLBMfile_name;
+ std::ofstream FxLBMfile;
+ std::ofstream FyLBMfile;
+ std::ofstream MzLBMfile;
+
+ std::stringstream Rhofile_name;
+ std::ofstream Rhofile;
+
+ std::stringstream Bodiesfile_name;
+ std::ofstream Bodiesfile;
+
+ std::stringstream NodeBoundaryfile_name;
+ std::ofstream NodeBoundaryfile;
+
+ std::stringstream NewNodefile_name;
+ std::ofstream NewNodefile;
+
+
+ /*--- VELOCITIES ---*/
+ if(SAVE_VELOCITY){
+ Vfile_name<<lbm_dir<<"/V"<<"_"; Vfile_name.width(10);
+ Vfile_name.fill('0'); Vfile_name<<iter_number;
+ Vfile.open(Vfile_name.str().c_str());
+
+ saveLBMdata = true;
+ }
+
+ if(SAVE_VELOCITYCOMP){
+ Vxfile_name<<lbm_dir<<"/Vx"<<"_";Vxfile_name.width(10);
+ Vyfile_name<<lbm_dir<<"/Vy"<<"_";Vyfile_name.width(10);
+
+ Vxfile_name.fill('0'); Vxfile_name<<iter_number;
+ Vyfile_name.fill('0'); Vyfile_name<<iter_number;
+
+ Vxfile.open(Vxfile_name.str().c_str());Vyfile.open(Vyfile_name.str().c_str());
+
+ saveLBMdata = true;
+ }
+
+ /*--- FORCES ---*/
+ if(SAVE_FORCES){
+ FxLBMfile_name<<lbm_dir<<"/Fx"<<"_";FxLBMfile_name.width(10);
+ FyLBMfile_name<<lbm_dir<<"/Fy"<<"_";FyLBMfile_name.width(10);
+ MzLBMfile_name<<lbm_dir<<"/Mz"<<"_";MzLBMfile_name.width(10);
+
+ FxLBMfile_name.fill('0'); FxLBMfile_name<<iter_number;
+ FyLBMfile_name.fill('0'); FyLBMfile_name<<iter_number;
+ MzLBMfile_name.fill('0'); MzLBMfile_name<<iter_number;
+
+ FxLBMfile.open(FxLBMfile_name.str().c_str());
+ FyLBMfile.open(FyLBMfile_name.str().c_str());
+ MzLBMfile.open(MzLBMfile_name.str().c_str());
+
+ saveLBMdata = true;
+ }
+
+ /*--- DENSITY ---*/
+ if(SAVE_RHO){
+ Rhofile_name<<lbm_dir<<"/Rho"<<"_";Rhofile_name.width(10);
+ Rhofile_name.fill('0'); Rhofile_name<<iter_number;
+ Rhofile.open(Rhofile_name.str().c_str());
+
+ saveLBMdata = true;
+ }
+
+ /*--- BODIES ---*/
+ if(SAVE_BODIES){
+ Bodiesfile_name<<lbm_dir<<"/Bodies"<<"_";Bodiesfile_name.width(10);
+ Bodiesfile_name.fill('0'); Bodiesfile_name<<iter_number;
+ Bodiesfile.open(Bodiesfile_name.str().c_str());
+
+ saveLBMdata = true;
+ }
+
+ /*--- BOUNDARY NODES ---*/
+ if(SAVE_NODEBD){
+ NodeBoundaryfile_name<<lbm_dir<<"/NodeBoundary"<<"_";
+ NodeBoundaryfile_name.width(10);
+ NodeBoundaryfile_name.fill('0'); NodeBoundaryfile_name<<iter_number;
+ NodeBoundaryfile.open(NodeBoundaryfile_name.str().c_str());
+
+ saveLBMdata = true;
+ }
+
+ /*--- NEW NODES ---*/
+ if(SAVE_NODEISNEW){
+ NewNodefile_name<<lbm_dir<<"/NewNode"<<"_";
+ NewNodefile_name.width(10);
+ NewNodefile_name.fill('0'); NewNodefile_name<<iter_number;
+ NewNodefile.open(NewNodefile_name.str().c_str());
+
+ saveLBMdata = true;
+ }
+
+
+
+ if(saveLBMdata){ //Condition to avoid to loop over all the lbm node if no lbm data are saved but only spheres data
+
+ cerr << "| Saving ("<<iter_number<<")"<<endl;
+
+ for (int nidx=0; nidx<Nx*Ny; nidx++){
+ if(SAVE_BODIES) Bodiesfile <<nodes[nidx].body_id<<" ";
+ if(SAVE_VELOCITY) Vfile<< c*nodes[nidx].velb.norm()<< " ";
+ if(SAVE_VELOCITYCOMP){Vxfile<< c*nodes[nidx].velb.x()<< " ";Vyfile<< c*nodes[nidx].velb.y()<< " ";}
+ if(SAVE_FORCES){
+ if(nodes[nidx].body_id>=0){
+ FxLBMfile<< 2.*Rho*c2*dx*(LBbodies[nodes[nidx].body_id].force.x()) << " ";
+ FyLBMfile<< 2.*Rho*c2*dx*(LBbodies[nodes[nidx].body_id].force.y()) << " ";
+ MzLBMfile<< 2.*Rho*c2*dx2*LBbodies[nodes[nidx].body_id].momentum.z() << " ";
+ }
+ else{FxLBMfile<< 0. << " ";FyLBMfile<< 0. << " ";MzLBMfile<< 0. << " ";}}
+ if(SAVE_RHO) Rhofile << Rho*nodes[nidx].rhob << " ";
+ if(SAVE_NODEBD) {
+ int tmp=0;
+ if(nodes[nidx].isObstacleBoundary) tmp=nodes[nidx].isObstacleBoundary;
+ if(nodes[nidx].isFluidBoundary) tmp= -1;
+ NodeBoundaryfile << tmp << " ";}
+ if(SAVE_NODEISNEW) {
+ int tmp=0;
+ if(nodes[nidx].isNewObstacle) tmp=nodes[nidx].isNewObstacle;
+ if(nodes[nidx].isNewFluid) tmp=-1;
+ NewNodefile << tmp << " ";}
+
+ if(nodes[nidx].i==Nx-1){
+ if(SAVE_BODIES) Bodiesfile<< endl;
+ if(SAVE_VELOCITY) Vfile<< endl;
+ if(SAVE_VELOCITYCOMP){Vxfile<< endl;Vyfile<< endl;}
+ if(SAVE_FORCES){FxLBMfile << endl;FyLBMfile << endl;MzLBMfile << endl;}
+ if(SAVE_RHO) Rhofile << endl;
+ if(SAVE_NODEBD) NodeBoundaryfile << endl;
+ if(SAVE_NODEISNEW) NewNodefile << endl;
+ }
+ }
+ }
+
+ std::stringstream cmd;
+ cmd<<"bzip2";
+ if(SAVE_BODIES) {Bodiesfile.close();cmd<<" "<<Bodiesfile_name.str().c_str();}
+ if(SAVE_VELOCITY) {Vfile.close();cmd<<" "<<Vfile_name.str().c_str();}
+ if(SAVE_VELOCITYCOMP) {Vxfile.close(); Vyfile.close();cmd<<" "<<Vxfile_name.str().c_str()<<" "<<Vyfile_name.str().c_str();}
+ if(SAVE_FORCES){
+ FxLBMfile.close();FyLBMfile.close();MzLBMfile.close();
+ cmd<<" "<<FxLBMfile_name.str().c_str()<<" "<<FyLBMfile_name.str().c_str()<<" "<<MzLBMfile_name.str().c_str();}
+ if(SAVE_RHO) {Rhofile.close();cmd<<" "<<Rhofile_name.str().c_str();}
+ if(SAVE_NODEBD) {NodeBoundaryfile.close();cmd<<" "<<NodeBoundaryfile_name.str().c_str();}
+ if(SAVE_NODEISNEW) {NewNodefile.close();cmd<<" "<<NewNodefile_name.str().c_str();}
+ //cmd<<"&";
+ if(COMPRESS_DATA) {if(std::system(cmd.str().c_str())) cerr<<"bzip error"<<endl;}
+
+
+ /*--------------------------------------------*/
+ /* Recording of properties of DEM objects */
+ /*--------------------------------------------*/
+ if(SAVE_SPHERES){ //condition to save spheres properties only if it is required by the operator, checking if NB_GRAINS > 0 is done earlier when recording configuration is done
+ spherefile_name.str("");
+ spherefile_name<<dem_dir<<"/spheres";
+ cerr <<" Write DEM data in "<<spherefile_name.str().c_str()<<endl;
+ spherefile_name<<"_";
+ spherefile_name.width(10);
+ spherefile_name.fill('0');
+ spherefile_name<<iter_number;
+ std::ofstream spherefile;
+ spherefile.open(spherefile_name.str().c_str());
+
+ for (unsigned int l=0;l<LBbodies.size();l++){
+ if(LBbodies[l].isPtc()){
+ if((!LBbodies[l].isEroded)&&(LBbodies[l].saveProperties)){
+ spherefile<<l<<" "; /* Id*/
+ spherefile<<dx*LBbodies[l].pos.x()<<" "<<dx*LBbodies[l].pos.y()<<" "<<dx*LBbodies[l].pos.z()<<" "; /* x, y, z*/
+ spherefile<<dx*LBbodies[l].radius<<" "; /* r*/
+ spherefile<<c*LBbodies[l].vel.x()<<" "<<c*LBbodies[l].vel.y()<<" "<<c*LBbodies[l].vel.z()<<" "; /* Vx, Vy, Vz*/
+ spherefile<<invdt*LBbodies[l].AVel.x()<<" "<<invdt*LBbodies[l].AVel.y()<<" "<<invdt*LBbodies[l].AVel.z()<<" "; /* Wx, Wy, Wz*/
+ spherefile<<LBbodies[l].Fh.x()<<" "<<LBbodies[l].Fh.y()<<" "<<LBbodies[l].Fh.z()<<" "; /* Fhx, Fhy, Fhz */
+ spherefile<<LBbodies[l].Mh.x()<<" "<<LBbodies[l].Mh.y()<<" "<<LBbodies[l].Mh.z()<<endl; /* Mhx, Mhy, Mhz */
+ }
+ }
+ }
+ spherefile.close();
+ }
+
+#ifdef LBM_VERBOSE
+cerr <<"END: HydrodynamicsLawLBM::save()"<<endl;
+#endif
+return;
+}
+
+void HydrodynamicsLawLBM::saveStats(int iter_number, Real timestep)
+{
+ #ifdef LBM_VERBOSE
+ cerr <<"START: HydrodynamicsLawLBM::saveStats()"<<endl;
+ #endif
+ cerr << "| Save stats ..."<<endl;
+ ofstream file(LBMmachFile.c_str(), ios::app);
+ file <<iter_number<<" "<<iter_number*timestep<<" "<<VmaxC<<" "<<VmaxC/c<<endl;
+ #ifdef LBM_VERBOSE
+ cerr <<"END: HydrodynamicsLawLBM::saveStats()"<<endl;
+ #endif
+ return;
+}
+//void HydrodynamicsLawLBM::saveEroded(int iter_number, Real timestep)
+//{
+//
+// cerr << "| Save Eroded Ptc ..."<<endl;
+// ofstream file(RemovedPtcFile.c_str(), ios::app);
+// file <<iter_number<<" "<<iter_number*timestep<<" "<<NumberPtcEroded<<" "<<Vr<<" "<<Vr/Vo<<" "<<FhTotale<<endl;
+//
+// return;
+//}
+void HydrodynamicsLawLBM::saveObservedPtc(int iter_number, Real timestep)
+{
+
+ cerr << "| Save Observed Ptc ..."<<endl;
+ ofstream file(ObservedPtcFile.c_str(), ios::app);
+ file <<iter_number<<" "<<iter_number*timestep<<" ";
+ file <<dx*LBbodies[ObservedPtc].pos.x()<<" "<<dx*LBbodies[ObservedPtc].pos.y()<<" "<<dx*LBbodies[ObservedPtc].pos.z()<<" ";
+ file <<dx*LBbodies[ObservedPtc].radius<<" ";
+ file <<c*LBbodies[ObservedPtc].vel.x()<<" "<<c*LBbodies[ObservedPtc].vel.y()<<" "<<c*LBbodies[ObservedPtc].vel.z()<<" ";
+ file <<invdt*LBbodies[ObservedPtc].AVel.x()<<" "<<invdt*LBbodies[ObservedPtc].AVel.y()<<" "<<invdt*LBbodies[ObservedPtc].AVel.z()<<" ";
+ file <<LBbodies[ObservedPtc].Fh.x()<<" "<<LBbodies[ObservedPtc].Fh.y()<<" "<<LBbodies[ObservedPtc].Fh.z()<<" ";
+ file <<LBbodies[ObservedPtc].Mh.x()<<" "<<LBbodies[ObservedPtc].Mh.y()<<" "<<LBbodies[ObservedPtc].Mh.z()<<endl;
+
+ return;
+}
+void HydrodynamicsLawLBM::saveObservedNode(int iter_number, Real timestep)
+{
+ return;
+}
+/*
+_________________________________________________________________________
+*/
+
+void HydrodynamicsLawLBM::createNewFiles()
+{
+
+ //spherefile_name<<dem_dir<<"/spheres"; //Useless here
+
+ ofstream file(LBMmachFile.c_str());
+ file <<"#iter_number\t time\t VmaxC\t VmaxC/c"<<endl;
+ file.close();
+
+ if(removingCriterion!=0){
+ ofstream file2(RemovedPtcFile.c_str());
+ file2 <<"#Iter time NumberPtcEroded Vr Vr/Vo FhTotale"<<endl;
+ file2.close();
+ }
+
+ if(SAVE_CONTACTINFO){
+ ofstream file3(LBMcontactsFile.c_str());
+ file3 <<"#Iter time NumberOfContact"<<endl;
+ file3.close();
+ }
+ //if(NB_GRAINS>0) {ofstream file3(spherefile_name.str().c_str());file3.close();} //For what this line is used for? It seems to work without
+
+ if(SAVE_OBSERVEDPTC){ //Condition to create observedPtc file only if the recording is required by the operator
+ ofstream file4(ObservedPtcFile.c_str());
+ file4 <<"#iter t x y z r Vx Vy Vz Wx Wy Wz Fx Fy Fz Mx My Mz"<<endl;
+ file4.close();
+ }
+
+ if(SAVE_OBSERVEDNODE){ //Condition to create observedNode file only if the recording is required by the operator
+ ofstream file5(ObservedNodeFile.c_str());file5.close();
+ }
+
+
+ return;
+}
+
+void HydrodynamicsLawLBM::createDirectories(bool dirLBM, bool dirDem, bool dirCntct)
+{
+
+ //bfs::create_directory(bfs::path(lbm_dir)); //ModLuc: to create only necessary directory
+ if(dirLBM) bfs::create_directory(bfs::path(lbm_dir));
+ //if(NB_GRAINS>0) bfs::create_directory(bfs::path(dem_dir)); //ModLuc: to create only necessary directory
+ if(dirDem) bfs::create_directory(bfs::path(dem_dir));
+ //if(SAVE_CONTACTINFO) bfs::create_directory(bfs::path(cntct_dir)); //ModLuc: to create only necessary directory
+ if(dirCntct) bfs::create_directory(bfs::path(cntct_dir));
+
+ return;
+}
+
+void HydrodynamicsLawLBM::writelogfile()
+{
+ ofstream file(LBMlogFile.c_str());
+ file <<"File format: 1"<<endl;
+ file <<"System parameters: "<<endl;
+ file <<"\t Lx0= "<<Lx0<<endl;
+ file <<"\t Ly0= "<<Ly0<<endl;
+ file <<"\t Lz0= "<<Lz0<<endl;
+ file <<"\t Wallthickness= "<<Wallthickness<<endl;
+ file <<"\t dP= "<<dP<<endl;
+ file <<"\t Nu= "<<Nu<<endl;
+ file <<"\t Rho= "<< Rho<<endl;
+ file <<"\t dx= "<< dx<<endl;
+ file <<"\t Nx= "<< Nx<<endl;
+ file <<"\t Ny= "<< Ny<<endl;
+ file <<"\t Nz= "<< Nz<<endl;
+ file <<"LBM parameters: "<<endl;
+ file <<"\t tau= "<<tau<<" omega= "<< omega <<endl;
+ file <<"\t IterMax= "<<IterMax<<endl;
+ file <<"\t SaveMode= "<<SaveMode<<endl;
+ file <<"\t IterSave= "<<IterSave<<endl;
+ file <<"\t SaveGridRatio= "<<SaveGridRatio<<endl;
+ file <<"\t DemIterLbmIterRatio= "<<DemIterLbmIterRatio<<endl;
+ file <<"\t ConvergenceThreshold= "<<ConvergenceThreshold<<endl;
+ file <<"\t Predicted Mach number (may be false)= "<<UMaxtheo/c<<endl;
+ file <<"\t LBM dt= "<<dt<<endl;
+ file <<"DEM parameters: "<<endl;
+ file <<"\t DEM dt= "<<DEMdt0<<endl;
+ file <<"\t DEM adjusted dt= "<<newDEMdt<<endl;
+ file <<"Particles: "<<endl;
+ file <<"\t InitialNumberOfDynamicParticles= "<<InitialNumberOfDynamicParticles<<endl;
+ file <<"\t NB_BODIES= "<<NB_BODIES<<" NB_GRAINS= "<<NB_GRAINS<<endl;
+ file <<"\t NB_DYNBODIES= "<<NB_DYNBODIES<<" NB_DYNGRAINS= "<<NB_DYNGRAINS<<endl;
+ file <<"\t Rmin / Rmax / Rmean = "<<dx*MinBodyRadius<<" / "<<dx*MaxBodyRadius<<" / "<<dx*MeanBodyRadius<<endl;
+ if(NB_GRAINS>0) file <<"\t NbNodePerPtc= "<<NbParticleNodes/NB_GRAINS<<endl;
+ else file <<"\t NbNodePerPtc= "<<-1<<endl;
+ file <<"\t Vo= "<<Vo<<endl;
+ file <<"Misc :"<<endl;
+ file <<"\t VbCutOff= "<<VbCutOff<<endl;
+ //file <<"\t CstBodyForceDensity= "<<CstBodyForceDensity<<endl;
+ file <<"Memory usage"<<endl;
+ file <<"\t Nodes= "<<nodes.size()<<endl;
+ file <<"\t links= "<<links.size()<<endl;
+
+ file.close();
+ return;
+}
+
+void HydrodynamicsLawLBM::modeTransition(){
+ cerr << "Mode transition "<<endl;
+ IterMax=1;
+ IterSubCyclingStart=-1;
+ LBM_ITER=1;
+ use_ConvergenceCriterion=false;
+ return;
+}
+
+void HydrodynamicsLawLBM::LbmEnd(){
+ if(MODE==1) IterMax=iter;
+ Omega::instance().stop();
+ Omega::instance().saveSimulation ("end.xml");
+}
+
+void HydrodynamicsLawLBM::CalculateAndApplyForcesAndTorquesOnBodies(bool mean,bool apply){
+ /*--------------------------------------------------------------------------------*/
+ /*---------------- APPLICATION OF HYDRODYNAMIC FORCES ON SPHERES -----------------*/
+ /*--------------------------------------------------------------------------------*/
+ if(mean) FhTotale=Vector3r::Zero();
+ FOREACH(const shared_ptr<Body>& b, *scene->bodies){
+ if(!b) continue;
+ const int id=b->getId();
+ //if ( ((b->isDynamic())&&(b->shape->getClassName()=="Sphere")) || (b->shape->getClassName()=="Box") ){ //ModLuc: remove the condition (b->isDynamic()) to be able to apply force and torque on non dynamic bodies, by this way hydrodynamic force and torque on bodies can be read through python even if bodies are non dynamic.
+ if ( (b->shape->getClassName()=="Sphere") || (b->shape->getClassName()=="Box") ){
+ if(mean){
+ LBbodies[id].fp=LBbodies[id].force;
+ LBbodies[id].force=0.5*(LBbodies[id].fp+LBbodies[id].fm);
+ LBbodies[id].fm=LBbodies[id].fp;
+ LBbodies[id].mp=LBbodies[id].momentum;
+ LBbodies[id].momentum=0.5*(LBbodies[id].mp+LBbodies[id].mm);
+ LBbodies[id].mm=LBbodies[id].mp;
+ LBbodies[id].Fh=2.*Rho*c2*dx*LBbodies[id].force;
+ LBbodies[id].Mh=2.*Rho*c2*dx2*LBbodies[id].momentum;
+ FhTotale=FhTotale+LBbodies[id].Fh;
+ }
+ if(apply){
+ scene->forces.addForce(id, LBbodies[id].Fh);
+ scene->forces.addTorque(id, LBbodies[id].Mh);
+ }
+ }
+ }
+ return;
+}
+YADE_PLUGIN((HydrodynamicsLawLBM));
+
+#endif //LBM_ENGINE
+
=== added file 'pkg/lbm/HydrodynamicsLawLBM.hpp'
--- pkg/lbm/HydrodynamicsLawLBM.hpp 1970-01-01 00:00:00 +0000
+++ pkg/lbm/HydrodynamicsLawLBM.hpp 2014-03-31 16:01:57 +0000
@@ -0,0 +1,302 @@
+/*************************************************************************
+* Copyright (C) 2009-2012 by Franck Lominé *
+* franck.lomine@xxxxxxxxxxxxxx *
+* *
+* This program is free software; it is licensed under the terms of the *
+* GNU General Public License v2 or later. See file LICENSE for details. *
+* *
+* Luc Scholtès luc.scholtes@xxxxxxxxxxxxxxxx *
+* and Luc Sibille luc.sibille@xxxxxxxxxxxxxxx also contributed to this *
+* code. *
+* *
+* Lominé F., Scholtès L., Sibille L., Poullain P. (2013) *
+* Modelling of fluid-solid interaction in granular media with coupled *
+* LB/DE methods: application to piping erosion. International Journal *
+* for Numerical and Analytical Methods in Geomechanics, 37(6):577-596 *
+* doi: 10.1002/nag.1109 *
+* *
+* Sibille L., Lominé F., Marot D. (2012) Investigation In Modelling *
+* Piping Erosion With a Coupled «Lattice Boltzmann – Discrete Element» *
+* Numerical Method. in Proc. 6th Int. Conference on Scour and Erosion *
+* (ICSE-6), pp. 93-100. *
+* *
+*************************************************************************/
+#ifdef LBM_ENGINE
+
+#pragma once
+#include<yade/pkg/lbm/LBMnode.hpp>
+#include<yade/pkg/lbm/LBMlink.hpp>
+#include<yade/pkg/lbm/LBMbody.hpp>
+#include<yade/core/GlobalEngine.hpp>
+
+
+class HydrodynamicsLawLBM : public GlobalEngine
+{
+ private :
+ std::ofstream ofile;
+
+ public :
+ bool firstRun, /*! = 1 if it is the first iteration during 1 YADE simulation*/
+ use_ConvergenceCriterion, /*! use stop condition based on the convergence criterion*/
+ SAVE_VELOCITY, /*! Switch to save node velocities*/
+ SAVE_VELOCITYCOMP, /*! Switch to save node velocities in each directions*/
+ SAVE_RHO, /*! Switch to save node densities*/
+ SAVE_FORCES, /*! Switch to save node force and momentum on grid*/
+ SAVE_BODIES, /*! Switch to save particle nodes*/
+ SAVE_NODEBD, /*! Switch to save fluid or solid boundary nodes*/
+ SAVE_NODEISNEW, /*! Switch to save new fluid/solid nodes*/
+ SAVE_DEBUGFILES, /*! Switch to save some debug data*/
+ SAVE_OBSERVEDPTC, /*! Switch to save properties of the observed particle*/
+ SAVE_OBSERVEDNODE, /*! Switch to save properties of the observed node*/
+ SAVE_CONTACTINFO, /*! Switch to save contact properties*/
+ SAVE_SPHERES, /*! Switch to save spheres properties*/
+ COMPRESS_DATA, /*! Switch to enable file compression*/
+ Xperiodicity, /*! Switch to activate lattice periodicity in x direction*/
+ Yperiodicity, /*! Switch to activate lattice periodicity in y direction*/
+ Zperiodicity; /*! Switch to activate lattice periodicity in z direction*/
+
+ int NB_BODIES, /*! Number of bodies*/
+ NB_GRAINS, /*! number of grains*/
+ NB_DYNGRAINS, /*! number of dynamic grains*/
+ NB_DYNBODIES, /*! number of dynamic bodies*/
+ NB_WALLS, /*! Number of walls*/
+ DEM_ITER, /*! Number of iteration of the DEM loop*/
+ LBM_ITER, /*! Number of iteration of the LBM loop*/
+ MODE, /*! 1->only a LBM loop, 2->lbm subcycling, 3->lbm subcycling after lbm loop*/
+ dim, /*! dimension*/
+ NbDir, /*! number of directions of the lattice model*/
+ NbNodes, /*! Total number of nodes*/
+ NbFluidNodes, /*! Number of fluid nodes*/
+ NbSolidNodes, /*! Number of solid nodes*/
+ NbParticleNodes, /*! Number of particle nodes*/
+ NbContacts, /*! Number of Contact*/
+ InitialNumberOfDynamicParticles,/*! Initial number of dynamic particles*/
+ NumberOfDynamicParticles, /*! Number of dynamic particles*/
+ Ny, /*! Number of grid divisions in y direction */
+ Nz, /*! Number of grid divisions in z direction */
+ NumberPtcEroded, /*! The bumber of eroded/removed particles*/
+ iter, /*! LBM Iteration number in current DEM loop (=1 in mode=2)*/
+ IdFirstSphere; /*! Id of the first sphere*/
+
+ Real height, /*! System height */
+ width, /*! System width */
+ depth, /*! System depth */
+ halfWallthickness, /*! Half Wall thickness */
+ Wallthickness, /*! Wall thickness */
+ cub, /*! A temporary variable to calculate equilibrium distribution function */
+ c, /*! Lattice speed */
+ c2, /*! The squared lattice speed*/
+ dx, /*! The lattice size*/
+ invdx, /*! 1 / lattice size*/
+ dx2, /*! The squared lattice size*/
+ uMax, /// TODO: PLEASE EXPLAIN uMax
+ cs, /*! c/sqrt(3) */
+ dt, /*! LBM timestep */
+ invdt, /*! one over LBM timestep */
+ nu, /*! LBM kinematic viscosity */
+ feqb, /*! Equilibrium distribution function*/
+ omega, /*! 1/tau */
+ Lx0, /*! LBM grid size in x direction*/
+ Ly0, /*! LBM grid size in y direction*/
+ Lz0, /*! LBM grid size in z direction*/
+ outside_limit, /*! the x coordinate of a point outside the system*/
+ DEMdt, /*! timestep for the DEM iteration*/
+ DEMdt0, /*! original timestep for the DEM iteration*/
+ newDEMdt, /*! the new timestep for the DEM iteration*/
+ Vr, /*! Volume of the removed particles*/
+ Vo, /*! Initial volume of dynamic particles */
+ VmeanFluidC, /*! Current mean fluid velocity */
+ PrevVmeanFluidC, /*! Previous mean fluid velocity */
+ PrevPrevVmeanFluidC, /*! Previous previous mean fluid velocity */
+ VmaxC, /*! Maximum velocity during the current time step*/
+ VminC, /*! Minimum velocity during the current time step*/
+ RhomaxC, /*! Maximum density during the current time step*/
+ RhominC, /*! Minimum density during the current time step*/
+ LBM_TIME, /*! The time ellapsed in the LB method*/
+ DEM_TIME, /*! The time ellapsed in the DE method*/
+ RhoTot, /*! Cumulative density*/
+ FmoyCur, /*! Mean force at the current LB iteration*/
+ FmoyPrev, /*! Mean force at the previous LB iteration*/
+ FmoyPrevPrev, /*! Mean force at 2 previous LB iteration*/
+ UMaxtheo, /// TODO: PLEASE EXPLAIN UMaxtheo
+ MaxBodyRadius, /*! Max radius of spheres*/
+ MinBodyRadius, /*! Min radius of spheres*/
+ MeanBodyRadius; /*! Mean radius of spheres*/
+
+ std::string LBMlogFile, /*! Name of the logfile */
+ LBMmachFile, /*! Name of the stat file */
+ LBMcontactsFile, /*! Name of the contact file */
+ RemovedPtcFile, /*! Name of the file to store removed particle informations*/
+ ObservedPtcFile, /*! Name of the file to store observed particle informations*/
+ ObservedNodeFile, /*! Name of the file to store observed particle informations*/
+ lbm_dir, /*! Directory name to save LBM files */
+ dem_dir, /*! Directory name to save DEM files */
+ cntct_dir; /*! Directory name to save contact properties */
+
+ std::stringstream spherefile_name; /*! Name of the file where sphere data are saved*/
+
+ vector<int> IdOfNextErodedPtc, /*! List of particles which will be eroded*/
+ opp; /*! opposite nodes */
+
+ vector<Real> w; /*! Weighting factor */
+
+ vector <LBMnode> nodes; /*! the LBM nodes*/
+ vector <LBMlink> links; /*! the LBM links*/
+ vector <LBMbody> LBbodies; /*! the LBM bodies*/
+
+ vector <Vector3r> eib; /*! node velocity directions*/
+
+ Vector3r FhTotale; ///Total hydrodynamic force
+
+ virtual ~HydrodynamicsLawLBM ();
+ virtual bool isActivated();
+ virtual void action();
+ void save(int iter_number, Real timestep);
+ void saveStats(int iter_number, Real timestep);
+ void saveEroded(int iter_number, Real timestep);
+ void saveContacts(int iter_number, Real timestep);
+ void saveObservedNode(int iter_number, Real timestep);
+ void saveObservedPtc(int iter_number, Real timestep);
+ void createNewFiles();
+ //void createDirectories(); // ModLuc: to create directories only if necessary
+ void createDirectories(bool dirLBM, bool dirDem, bool dirCntct);
+ void writelogfile();
+ void modeTransition();
+ void LbmEnd();
+ void CalculateAndApplyForcesAndTorquesOnBodies(bool mean,bool apply);
+
+ YADE_CLASS_BASE_DOC_ATTRS_CTOR(HydrodynamicsLawLBM,GlobalEngine,"Engine to simulate fluid flow (with the lattice Boltzmann method) with a coupling with the discrete element method.\n If you use this Engine, please cite and refer to F. Lominé et al. International Journal For Numerical and Analytical Method in Geomechanics, 2012, doi: 10.1002/nag.1109",
+
+ ((int,WallYm_id,0,,"Identifier of the Y- wall"))
+ ((bool,useWallYm,true,,"Set true if you want that the LBM see the wall in Ym"))
+ ((int,YmBCType,2,,"Boundary condition for the wall in Ym (-1: unused, 1: pressure condition, 2: velocity condition)."))
+ ((Vector3r,YmBcVel,Vector3r::Zero(),,"(!!! not fully implemented !!) The velocity imposed at the boundary"))
+ ((Real,YmBcRho,-1,,"(!!! not fully implemented !!) The density imposed at the boundary"))
+ ((int,WallYp_id,1,,"Identifier of the Y+ wall"))
+ ((bool,useWallYp,true,,"Set true if you want that the LBM see the wall in Yp"))
+ ((int,YpBCType,2,,"Boundary condition for the wall in Yp (-1: unused, 1: pressure condition, 2: velocity condition)."))
+ ((Vector3r,YpBcVel,Vector3r::Zero(),,"(!!! not fully implemented !!) The velocity imposed at the boundary"))
+ ((Real,YpBcRho,-1,,"(!!! not fully implemented !!) The density imposed at the boundary"))
+ ((int,WallXm_id,2,,"Identifier of the X- wall"))
+ ((bool,useWallXm,false,,"Set true if you want that the LBM see the wall in Xm"))
+ ((int,XmBCType,1,,"Boundary condition for the wall in Xm (-1: unused, 1: pressure condition, 2: velocity condition)."))
+ ((Vector3r,XmBcVel,Vector3r::Zero(),,"(!!! not fully implemented !!) The velocity imposed at the boundary"))
+ ((Real,XmBcRho,-1,,"(!!! not fully implemented !!) The density imposed at the boundary"))
+ ((int,WallXp_id,3,,"Identifier of the X+ wall"))
+ ((bool,useWallXp,false,,"Set true if you want that the LBM see the wall in Xp"))
+ ((int,XpBCType,1,,"Boundary condition for the wall in Xp (-1: unused, 1: pressure condition, 2: velocity condition)."))
+ ((Vector3r,XpBcVel,Vector3r::Zero(),,"(!!! not fully implemented !!) The velocity imposed at the boundary"))
+ ((Real,XpBcRho,-1,,"(!!! not fully implemented !!) The density imposed at the boundary"))
+ ((int,WallZp_id,5,,"Identifier of the Z+ wall"))
+ ((bool,useWallZp,false,,"Set true if you want that the LBM see the wall in Zp"))
+ ((int,ZpBCType,-1,,"Boundary condition for the wall in Zp (-1: unused, 1: pressure condition, 2: velocity condition)."))
+ ((Vector3r,ZpBcVel,Vector3r::Zero(),,"(!!! not fully implemented !!) The velocity imposed at the boundary"))
+ ((Real,zpBcRho,-1,,"(!!! not fully implemented !!) The density imposed at the boundary"))
+ ((int,WallZm_id,4,,"Identifier of the Z- wall"))
+ ((bool,useWallZm,false,,"Set true if you want that the LBM see the wall in Zm"))
+ ((int,ZmBCType,-1,,"Boundary condition for the wall in Zm (-1: unused, 1: pressure condition, 2: velocity condition)."))
+ ((Vector3r,ZmBcVel,Vector3r::Zero(),,"(!!! not fully implemented !!) The velocity imposed at the boundary"))
+ ((Real,ZmBcRho,-1,,"(!!! not fully implemented !!) The density imposed at the boundary"))
+ ((int,XmYmZpBCType,2,,"Boundary condition for the corner node XmYmZp (-1: unused, 1: pressure condition, 2: velocity condition)."))
+ ((int,XmYpZpBCType,2,,"Boundary condition for the corner node XmYpZp (-1: unused, 1: pressure condition, 2: velocity condition)."))
+ ((int,XpYmZpBCType,2,,"Boundary condition for the corner node XpYmZp (-1: unused, 1: pressure condition, 2: velocity condition)."))
+ ((int,XpYpZpBCType,2,,"Boundary condition for the corner node XpYpZp (-1: unused, 1: pressure condition, 2: velocity condition)."))
+ ((int,XmYmZmBCType,-1,,"Boundary condition for the corner node XmYmZm (not used with d2q9, -1: unused, 1: pressure condition, 2: velocity condition)."))
+ ((int,XmYpZmBCType,-1,,"Boundary condition for the corner node XmYpZm (not used with d2q9, -1: unused, 1: pressure condition, 2: velocity condition)."))
+ ((int,XpYmZmBCType,-1,,"Boundary condition for the corner node XpYmZm (not used with d2q9, -1: unused, 1: pressure condition, 2: velocity condition)."))
+ ((int,XpYpZmBCType,-1,,"Boundary condition for the corner node XpYpZm (not used with d2q9, -1: unused, 1: pressure condition, 2: velocity condition)."))
+
+ ((int,defaultLbmInitMode,0,,"Switch between the two initialisation methods"))
+ ((Vector3r,dP,Vector3r(0.,0.,0.),,"Pressure difference between input and output"))
+ ((Real,Rho,1000.,,"Fluid density"))
+ ((Real,Nu,0.000001,,"Fluid kinematic viscosity"))
+ ((Real,tau,0.6,,"Relaxation time"))
+ ((int,Nx,1000,,"The number of grid division in x direction"))
+ ((int,IterMax,1,,"This variable can be used to do several LBM iterations during one DEM iteration. "))
+ ((int,IterPrint,1,,"Print info on screen every IterPrint iterations"))
+ ((int,SaveMode,1,,"Save Mode (1-> default, 2-> in time (not yet implemented)"))
+ ((int,IterSave,100,,"Data are saved every IterSave LBM iteration (or see TimeSave)"))
+ ((Real,TimeSave,-1,,"Data are saved at constant time interval (or see IterSave)"))
+ ((int,SaveGridRatio,1,,"Grid data are saved every SaveGridRatio * IterSave LBM iteration (with SaveMode=1)"))
+ ((int,IterSubCyclingStart,-1,,"Iteration number when the subcycling process starts"))
+ ((int,DemIterLbmIterRatio,-1,,"Ratio between DEM and LBM iterations for subcycling"))
+ ((bool,EngineIsActivated,true,,"To activate (or not) the engine"))
+ ((bool,applyForcesAndTorques,true,,"Switch to apply forces and torques"))
+ ((int,ObservedNode,-1,,"The identifier of the node that will be observed (-1 means none)"))
+ ((int,ObservedPtc,-1,,"The identifier of the particle that will be observed (-1 means the first one)"))
+ ((Real,RadFactor,1.0,,"The radius of DEM particules seen by the LBM is the real radius of particules*RadFactor"))
+ ((Real,ConvergenceThreshold,0.000001,,""))
+ ((std::string,LBMSavedData," ",,"a list of data that will be saved. Can use velocity,velXY,forces,rho,bodies,nodeBD,newNode,observedptc,observednode,contacts,spheres,bz2"))
+ ((std::string,periodicity," ",,"periodicity"))
+ ((std::string,bc," ",,"Boundary condition"))
+ ((std::string,model,"d2q9",,"The LB model. Until now only d2q9 is implemented"))
+ ((int,removingCriterion ,0,,"Criterion to remove a sphere (1->based on particle position, 2->based on particle velocity"))
+ ((Real,VelocityThreshold,-1.,,"Velocity threshold when removingCriterion=2"))
+ ((Real,EndTime,-1,,"the time to stop the simulation"))
+ ((Vector3r,CstBodyForce,Vector3r::Zero(),,"A constant body force (=that does not vary in time or space, otherwise the implementation introduces errors)"))
+ ((Real,VbCutOff,-1,,"the minimum boundary velocity that is taken into account"))
+ ,
+ firstRun = true;
+ omega = 1.0/tau;
+ DEM_TIME = 0.;
+ LBM_TIME = 0.;
+ iter = 0;
+ use_ConvergenceCriterion = true;
+ MODE=0;
+ dim=0;
+ NbDir=0;
+ NbNodes=0;
+ NbFluidNodes=0;
+ NbSolidNodes=0;
+ NbParticleNodes=0;
+ NbContacts=0;
+ InitialNumberOfDynamicParticles=0;
+ NumberOfDynamicParticles=0;
+ NumberPtcEroded=0;
+ Vr=0.;
+ Vo=0.;
+ MaxBodyRadius=-1000000.;
+ MinBodyRadius=1000000.;
+ MeanBodyRadius=0.;
+ lbm_dir="lbm-nodes";
+ dem_dir="dem-bodies";
+ cntct_dir="contacts";
+ LBMlogFile ="LBM.log";
+ LBMmachFile ="LBM.mach";
+ LBMcontactsFile ="LBM.cntct";
+ RemovedPtcFile="eroded.dat";
+ ObservedPtcFile="observedPtc.dat";
+ ObservedNodeFile="observedNode.dat";
+ COMPRESS_DATA = false;
+ SAVE_VELOCITY = false;
+ SAVE_VELOCITYCOMP = false;
+ SAVE_RHO = false;
+ SAVE_FORCES = false;
+ SAVE_BODIES = false;
+ SAVE_NODEBD = false;
+ SAVE_NODEISNEW = false;
+ SAVE_DEBUGFILES = false;
+ SAVE_OBSERVEDPTC= false;
+ SAVE_OBSERVEDNODE=false;
+ SAVE_CONTACTINFO =false;
+ SAVE_SPHERES = false; //to save spheres_* files only if it is required by the operator
+ Xperiodicity = false;
+ Yperiodicity = false;
+ Zperiodicity = false;
+ Ny = 0;Nz = 0;
+ FmoyCur=0.;FmoyPrev=0.;
+ FmoyPrevPrev=0.;VmeanFluidC=0.;PrevVmeanFluidC=0.;PrevPrevVmeanFluidC=0.;
+ LBM_ITER=0;
+ DEM_ITER=0;
+ IdFirstSphere=-1;
+ timingDeltas=shared_ptr<TimingDeltas>(new TimingDeltas);
+
+ );
+ DECLARE_LOGGER;
+};
+
+
+REGISTER_SERIALIZABLE(HydrodynamicsLawLBM);
+
+#endif //LBM_ENGINE
+
=== added file 'pkg/lbm/LBMbody.cpp'
--- pkg/lbm/LBMbody.cpp 1970-01-01 00:00:00 +0000
+++ pkg/lbm/LBMbody.cpp 2014-03-31 16:01:57 +0000
@@ -0,0 +1,16 @@
+/*************************************************************************
+* Copyright (C) 2009-2012 by Franck Lominé *
+* franck.lomine@xxxxxxxxxxxxxx *
+* *
+* This program is free software; it is licensed under the terms of the *
+* GNU General Public License v2 or later. See file LICENSE for details. *
+* *
+*************************************************************************/
+#ifdef LBM_ENGINE
+
+#include "LBMbody.hpp"
+
+YADE_PLUGIN((LBMbody));
+LBMbody::~LBMbody(){};
+
+#endif //LBM_ENGINE
=== added file 'pkg/lbm/LBMbody.hpp'
--- pkg/lbm/LBMbody.hpp 1970-01-01 00:00:00 +0000
+++ pkg/lbm/LBMbody.hpp 2014-03-31 16:01:57 +0000
@@ -0,0 +1,46 @@
+/*************************************************************************
+* Copyright (C) 2009-2012 by Franck Lominé *
+* franck.lomine@xxxxxxxxxxxxxx *
+* *
+* This program is free software; it is licensed under the terms of the *
+* GNU General Public License v2 or later. See file LICENSE for details. *
+* *
+*************************************************************************/
+#ifdef LBM_ENGINE
+
+#pragma once
+#include<yade/lib/serialization/Serializable.hpp>
+#include<yade/lib/multimethods/Indexable.hpp>
+
+class LBMbody: public Serializable{
+ public:
+ virtual ~LBMbody();
+ //Real radius(){return ext[0];}
+ bool isBox(){if(type==1)return true; else return false;}
+ bool isPtc(){if(type==2)return true; else return false;}
+ void setAsPtc(){type=2;}
+ void setAsBox(){type=1;}
+
+ YADE_CLASS_BASE_DOC_ATTRS_CTOR(LBMbody,Serializable,
+ "Body class for Lattice Boltzmann Method ",
+ ((Vector3r,force,Vector3r::Zero(),,"Hydrodynamic force, need to be reinitialized (LB unit)"))
+ ((Vector3r,fm,Vector3r::Zero(),,"Hydrodynamic force (LB unit) at t-0.5dt"))
+ ((Vector3r,fp,Vector3r::Zero(),,"Hydrodynamic force (LB unit) at t+0.5dt"))
+ ((Vector3r,momentum,Vector3r::Zero(),,"Hydrodynamic momentum,need to be reinitialized (LB unit)"))
+ ((Vector3r,mm,Vector3r::Zero(),,"Hydrodynamic momentum (LB unit) at t-0.5dt"))
+ ((Vector3r,mp,Vector3r::Zero(),,"Hydrodynamic momentum (LB unit) at t+0.5dt"))
+ ((Vector3r,pos,Vector3r::Zero(),,"Position of body"))
+ ((Vector3r,vel,Vector3r::Zero(),,"Velocity of body"))
+ ((Vector3r,AVel,Vector3r::Zero(),,"Angular velocity of body"))
+ ((Vector3r,Fh,Vector3r::Zero(),,"Hydrodynamical force on body"))
+ ((Vector3r,Mh,Vector3r::Zero(),,"Hydrodynamical momentum on body"))
+ ((Real,radius,-1000.,,"Radius of body (for sphere)"))
+ ((bool,isEroded,false,,"Hydrodynamical force on body"))
+ ((bool,saveProperties,false,,"To save properties of the body"))
+ ((short int,type,-1,," "))
+ ,
+ );
+};
+REGISTER_SERIALIZABLE(LBMbody);
+
+#endif //LBM_ENGINE
=== added file 'pkg/lbm/LBMlink.cpp'
--- pkg/lbm/LBMlink.cpp 1970-01-01 00:00:00 +0000
+++ pkg/lbm/LBMlink.cpp 2014-03-31 16:01:57 +0000
@@ -0,0 +1,17 @@
+/*************************************************************************
+* Copyright (C) 2009-2012 by Franck Lominé *
+* franck.lomine@xxxxxxxxxxxxxx *
+* *
+* This program is free software; it is licensed under the terms of the *
+* GNU General Public License v2 or later. See file LICENSE for details. *
+* *
+*************************************************************************/
+#ifdef LBM_ENGINE
+
+#include "LBMlink.hpp"
+
+YADE_PLUGIN((LBMlink));
+LBMlink::~LBMlink(){};
+void LBMlink::ReinitDynamicalProperties(){sid=-1;fid=-1;idx_sigma_i=-1;isBd=false;VbMid=Vector3r::Zero();DistMid=Vector3r::Zero();ct=0.;return;}
+
+#endif //LBM_ENGINE
=== added file 'pkg/lbm/LBMlink.hpp'
--- pkg/lbm/LBMlink.hpp 1970-01-01 00:00:00 +0000
+++ pkg/lbm/LBMlink.hpp 2014-03-31 16:01:57 +0000
@@ -0,0 +1,37 @@
+/*************************************************************************
+* Copyright (C) 2009-2012 by Franck Lominé *
+* franck.lomine@xxxxxxxxxxxxxx *
+* *
+* This program is free software; it is licensed under the terms of the *
+* GNU General Public License v2 or later. See file LICENSE for details. *
+* *
+*************************************************************************/
+#ifdef LBM_ENGINE
+
+#pragma once
+#include<yade/lib/serialization/Serializable.hpp>
+#include<yade/lib/multimethods/Indexable.hpp>
+
+class LBMlink: public Serializable{
+ public:
+ void ReinitDynamicalProperties();
+ virtual ~LBMlink();
+
+ YADE_CLASS_BASE_DOC_ATTRS_CTOR(LBMlink,Serializable,
+ "Link class for Lattice Boltzmann Method ",
+ ((int,sid,-1,,"Solid node identifier "))
+ ((int,fid,-1,,"Fluid node identifier "))
+ ((short int,i,-1,,"direction index of the link"))
+ ((int,nid1,-1,,"fixed node identifier"))
+ ((int,nid2,-1,,"fixed node identifier or -1 if node points outside"))
+ ((short int,idx_sigma_i,-1,,"sigma_i direction index (Fluid->Solid)"))
+ ((bool,isBd,false,,"True if it is a boundary link"))
+ ((bool,PointingOutside,false,,"True if it is a link pointing outside to the system (from a fluid or solid node)"))
+ ((Vector3r,VbMid,Vector3r::Zero(),,"Velocity of boundary at midpoint"))
+ ((Vector3r,DistMid,Vector3r::Zero(),,"Distance between middle of the link and mass center of body"))
+ ((Real,ct,0.,,"Coupling term in modified bounce back rule")),
+ );
+};
+REGISTER_SERIALIZABLE(LBMlink);
+
+#endif //LBM_ENGINE
=== added file 'pkg/lbm/LBMnode.cpp'
--- pkg/lbm/LBMnode.cpp 1970-01-01 00:00:00 +0000
+++ pkg/lbm/LBMnode.cpp 2014-03-31 16:01:57 +0000
@@ -0,0 +1,120 @@
+/*************************************************************************
+* Copyright (C) 2009-2012 by Franck Lominé *
+* franck.lomine@xxxxxxxxxxxxxx *
+* *
+* This program is free software; it is licensed under the terms of the *
+* GNU General Public License v2 or later. See file LICENSE for details. *
+* *
+*************************************************************************/
+#ifdef LBM_ENGINE
+
+#include "LBMnode.hpp"
+
+YADE_PLUGIN((LBMnode));
+LBMnode::~LBMnode(){};
+
+void LBMnode::MixteBC(string lbmodel,Real density, Vector3r U, string where){
+ Real rhoVx=density*U.x();
+ Real rhoVy=density*U.y();
+ if(!strcmp(lbmodel.c_str(), "d2q9" )){
+ if(!strcmp(where.c_str(), "Xm" )){
+ f[1]=f[3]+(2./3.)*rhoVx;
+ f[5]=f[7]-0.5*(f[2]-f[4])+(1./6.)*rhoVx+ 0.5*rhoVy;
+ f[8]=f[6]+0.5*(f[2]-f[4])+(1./6.)*rhoVx- 0.5*rhoVy;
+ }
+ else if(!strcmp(where.c_str(), "Xp" )){
+ f[3]=f[1]-(2./3.)*rhoVx;
+ f[7]=f[5]+0.5*(f[2]-f[4])-(1./6.)*rhoVx - 0.5*rhoVy;
+ f[6]=f[8]-0.5*(f[2]-f[4])-(1./6.)*rhoVx + 0.5*rhoVy;
+ }
+ else if(!strcmp(where.c_str(), "Ym" )){
+ f[2]=f[4]+(2./3.)*rhoVy;
+ f[5]=f[7]-0.5*(f[1]-f[3])+0.5*rhoVx + (1./6.)*rhoVy;
+ f[6]=f[8]+0.5*(f[1]-f[3])-0.5*rhoVx + (1./6.)*rhoVy;
+ }
+ else if(!strcmp(where.c_str(), "Yp" )){
+ f[4]=f[2]-(2./3.)*rhoVy;
+ f[7]=f[5]+0.5*(f[1]-f[3])-0.5*rhoVx- (1./6.)*rhoVy;
+ f[8]=f[6]-0.5*(f[1]-f[3])+0.5*rhoVx -(1./6.)*rhoVy;
+ }
+ else if(!strcmp(where.c_str(), "XmYmZp" )){
+ f[1]=f[3]+(2./3.)*rhoVx;
+ f[2]=f[4]+(2./3.)*rhoVy;
+ f[5]=f[7] + (1./6.)*density*(U.x()+U.y());
+ f[6]=0.5*(density*(1.-U.x() -(2./3.)*U.y())-f[0]-2.*(f[3]+f[4]+f[7]));
+ f[8]=0.5*(density*(1.-(2./3.)*U.x() -U.y())-f[0]-2.*(f[3]+f[4]+f[7]));
+ }
+ else if(!strcmp(where.c_str(), "XmYpZp" )){
+ f[1]=f[3]+(2./3.)*rhoVx;
+ f[4]=f[2]-(2./3.)*rhoVy;
+ f[5]=0.5*(density*(1.-(2./3.)*U.x()+U.y())-f[0]-2.*(f[2]+f[3]+f[6]));
+ f[7]=0.5*(density*(1.-U.x()+(2./3.)*U.y())-f[0]-2.*(f[2]+f[3]+f[6]));
+ f[8]=f[6]+(1./6.)*density*(U.x()-U.y());
+ }
+ else if(!strcmp(where.c_str(), "XpYmZp" )){
+ f[2]=f[4]+(2./3.)*rhoVy;
+ f[3]=f[1]-(2./3.)*rhoVx;
+ f[5]=0.5*(density*(1.+U.x()-(2./3.)*U.y())-f[0]-2.*(f[1]+f[4]+f[8]));
+ f[6]=f[8]-(1./6.)*density*(U.x()-U.y());
+ f[7]=0.5*(density*(1.+(2./3.)*U.x()-U.y())-f[0]-2.*(f[1]+f[4]+f[8]));
+ }
+ else if(!strcmp(where.c_str(), "XpYpZp" )){
+ f[3]=f[1]-(2./3.)*rhoVx;
+ f[4]=f[2]-(2./3.)*rhoVy;
+ f[6]=0.5*(density*(1.+(2./3.)*U.x()+U.y())-f[0]-2.*(f[1]+f[2]+f[5]));
+ f[7]=f[5]-(1./6.)*density*(U.x()+U.y());
+ f[8]=0.5*(density*(1.+U.x()+(2./3.)*U.y())-f[0]-2.*(f[1]+f[2]+f[5]));
+ }
+ else {exit(-1);}
+ }else {exit(-1);}
+ return;
+}
+
+
+bool LBMnode::checkIsNewObstacle(){
+ if(isObstacle){
+ if(!wasObstacle) {isNewObstacle=true;wasObstacle=true;}
+ else {isNewObstacle=false;wasObstacle=true;}
+ return(isNewObstacle);
+ } else return(false);
+}
+
+bool LBMnode::checkIsNewFluid(){
+ if(!isObstacle){
+ if(wasObstacle) {isNewFluid=true;wasObstacle=false;}
+ else {isNewFluid=false;wasObstacle=false;}
+ return(isNewFluid);
+ } else return(false);
+}
+
+void LBMnode::DispatchBoundaryConditions(int SizeNx,int SizeNy,int SizeNz){
+ applyBC =false;
+ applyXmBC =false;
+ applyYmXmBC =false;
+ applyYpXmBC =false;
+ applyXpBC =false;
+ applyYmXpBC =false;
+ applyYpXpBC =false;
+ applyYpBC =false;
+ applyYmBC =false;
+ if((i==0)&&(j>0)&&(j<SizeNy-1)) {applyXmBC =true; applyBC=true;}
+ if((i==0)&&(j==0)) {applyYmXmBC =true; applyBC=true;}
+ if((i==0)&&(j==SizeNy-1)) {applyYpXmBC =true; applyBC=true;}
+ if((i==SizeNx-1)&&(j>0)&&(j<SizeNy-1)) {applyXpBC =true; applyBC=true;}
+ if((i==SizeNx-1)&&(j==0)) {applyYmXpBC =true; applyBC=true;}
+ if((i==SizeNx-1)&&(j==SizeNy-1)) {applyYpXpBC =true; applyBC=true;}
+ if((i>0)&&(i<SizeNx-1)&&(j==0)) {applyYmBC =true; applyBC=true;}
+ if((i>0)&&(i<SizeNx-1)&&(j==SizeNy-1)) {applyYpBC =true; applyBC=true;}
+return;
+}
+
+
+void LBMnode::SetCellIndexesAndPosition(int indI, int indJ, int indK){
+ i=indI;
+ j=indJ;
+ k=indK;
+ posb=Vector3r((Real) indI,(Real) indJ,(Real) indK);
+ return;
+}
+
+#endif //LBM_ENGINE
=== added file 'pkg/lbm/LBMnode.hpp'
--- pkg/lbm/LBMnode.hpp 1970-01-01 00:00:00 +0000
+++ pkg/lbm/LBMnode.hpp 2014-03-31 16:01:57 +0000
@@ -0,0 +1,67 @@
+/*************************************************************************
+* Copyright (C) 2009-2012 by Franck Lominé *
+* franck.lomine@xxxxxxxxxxxxxx *
+* *
+* This program is free software; it is licensed under the terms of the *
+* GNU General Public License v2 or later. See file LICENSE for details. *
+* *
+*************************************************************************/
+#ifdef LBM_ENGINE
+
+#pragma once
+#include<yade/lib/serialization/Serializable.hpp>
+#include<yade/lib/multimethods/Indexable.hpp>
+
+class LBMnode: public Serializable{
+ public:
+ int
+ i, /*! node index in x direction */
+ j, /*! node index in y direction */
+ k, /*! node index in z direction */
+ body_id; /*! the node belongs to the body indexed with body_id*/
+
+ short int IsolNb; /*! number of boundary links of a fluid boundary nodes*/
+
+ bool
+ isObstacle, /*! the node belongs to an obstacle */
+ isObstacleBoundary, /*! the node is an obstacle boundary */
+ isFluidBoundary, /*! the node is a fluid boundary */
+ wasObstacle, /*! the node was an obstacle */
+ isNewObstacle, /*! the node is an new obstacle node */
+ isNewFluid, /*! the node is an new fluid node */
+ applyBC, /*! the node is subject to the one boundary condition */
+ applyXmBC, /*! the node is subject to the left boundary condition */
+ applyXpBC, /*! the node is subject to the right boundary condition */
+ applyYpBC, /*! the node is subject to the top boundary condition */
+ applyYmBC, /*! the node is subject to the bottom boundary condition */
+ applyZpBC, /*! the node is subject to the front boundary condition NOT USED NOW*/
+ applyZmBC, /*! the node is subject to the back boundary condition NOT USED NOW*/
+ applyYmXmBC, /*! the node is subject to the bottom-left boundary condition */
+ applyYpXmBC, /*! the node is subject to the top-left boundary condition */
+ applyYmXpBC, /*! the node is subject to the bottom-right boundary condition */
+ applyYpXpBC; /*! the node is subject to the top-right boundary condition */
+
+ Vector3r posb, /*! the node position */
+ velb; /*! the node velocity */
+
+ Real rhob; /*! the node density */ /*! the node density */
+ vector<int> neighbour_id; /*! list of adjacent nodes */
+ vector<int> links_id; /*! list of links */
+ vector<Real> f;
+ vector<Real> fprecol;
+ vector<Real> fpostcol;
+
+ void DispatchBoundaryConditions(int SizeNx,int SizeNy,int SizeNz);
+ bool checkIsNewFluid();
+ bool checkIsNewObstacle();
+ void MixteBC(string lbmodel,Real density, Vector3r U, string where);
+ void SetCellIndexesAndPosition(int indI, int indJ, int indK);
+ void setAsObstacle(){isObstacle=true;}
+ void setAsFluid(){isObstacle=false;}
+
+ virtual ~LBMnode();
+ YADE_CLASS_BASE_DOC(LBMnode,Serializable,"Node class for Lattice Boltzmann Method ");
+};
+REGISTER_SERIALIZABLE(LBMnode);
+
+#endif //LBM_ENGINE
