update stride selector

readme.md files update

CMakeLists.txt

@@ -3,30 +3,17 @@ cmake_minimum_required(VERSION 3.16 FATAL_ERROR)
 # set the project name and version
 project(phasicFlow VERSION 1.0 )
 
-set(CMAKE_CXX_STANDARD 17 CACHE STRING "" FORCE)
+set(CMAKE_CXX_STANDARD 20 CACHE STRING "" FORCE)
 set(CMAKE_CXX_STANDARD_REQUIRED True)
 set(CMAKE_INSTALL_PREFIX ${phasicFlow_SOURCE_DIR} CACHE PATH "Install path of phasicFlow" FORCE)
-set(CMAKE_BUILD_TYPE Debug CACHE STRING "build type" FORCE)
+set(CMAKE_BUILD_TYPE Release CACHE STRING "build type")
 set(BUILD_SHARED_LIBS ON CACHE BOOL "Build using shared libraries" FORCE)
 mark_as_advanced(FORCE var BUILD_SHARED_LIBS)
 
-message(STATUS ${CMAKE_INSTALL_PREFIX})
+message(STATUS "Install prefix is:" ${CMAKE_INSTALL_PREFIX})
 
 include(cmake/globals.cmake)
 
-#Kokkos directory to be included 
-set(Kokkos_Source_DIR)
-
-if(DEFINED ENV{Kokkos_DIR})
-   set(Kokkos_Source_DIR $ENV{Kokkos_DIR})
-else()
-  set(Kokkos_Source_DIR $ENV{HOME}/Kokkos/kokkos)
-endif()
-message(STATUS "Kokkos source directory is ${Kokkos_Source_DIR}")
-add_subdirectory(${Kokkos_Source_DIR} ./kokkos)
-Kokkos_cmake_settings()
-
-
 option(pFlow_STD_Parallel_Alg "Use TTB std parallel algorithms" ON)
 option(pFlow_Build_Serial "Build phasicFlow and backends for serial execution" OFF)
 option(pFlow_Build_OpenMP "Build phasicFlow and backends for OpenMP execution" OFF)
@@ -34,6 +21,8 @@ option(pFlow_Build_Cuda   "Build phasicFlow and backends for Cuda execution" OFF
 option(pFlow_Build_Double "Build phasicFlow with double precision floating-oint variables" ON)
 option(pFlow_Build_MPI    "Build for MPI parallelization. This will enable multi-gpu run, CPU run on clusters (distributed memory machine). Use this combination Cuda+MPI, OpenMP + MPI or Serial+MPI " OFF)
 
+#for installing the required packages
+include(cmake/preReq.cmake)
 
 if(pFlow_Build_Serial)
 	set(Kokkos_ENABLE_SERIAL ON CACHE BOOL "Serial execution" FORCE) 
@@ -46,7 +35,8 @@ elseif(pFlow_Build_OpenMP )
 	set(Kokkos_ENABLE_OPENMP ON CACHE BOOL "OpenMP execution" FORCE)
 	set(Kokkos_ENABLE_CUDA OFF CACHE BOOL "Cuda execution" FORCE)
 	set(Kokkos_ENABLE_CUDA_LAMBDA OFF CACHE BOOL "Cuda execution" FORCE)
-	set(Kokkos_DEFAULT_HOST_PARALLEL_EXECUTION_SPACE SERIAL CACHE STRING "" FORCE)
+	set(Kokkos_DEFAULT_HOST_PARALLEL_EXECUTION_SPACE Serial CACHE STRING "" FORCE)
+	set(Kokkos_DEFAULT_DEVICE_PARALLEL_EXECUTION_SPACE OpenMP CACHE STRING "" FORCE)
 	set(Kokkos_ENABLE_CUDA_CONSTEXPR OFF CACHE BOOL "Enable constexpr on cuda code" FORCE)
 elseif(pFlow_Build_Cuda)
 	set(Kokkos_ENABLE_SERIAL ON CACHE BOOL "Serial execution" FORCE) 
@@ -65,6 +55,7 @@ include(cmake/makeExecutableGlobals.cmake)
 
 configure_file(phasicFlowConfig.H.in phasicFlowConfig.H)
 set(CMAKE_EXPORT_COMPILE_COMMANDS ON)
+
 #add a global include directory 
 include_directories(src/setHelpers src/demComponent "${PROJECT_BINARY_DIR}")
 

DEMSystems/grainDEMSystem/grainDEMSystem.cpp

@@ -78,9 +78,14 @@ pFlow::grainDEMSystem::grainDEMSystem(
 	REPORT(0)<< "\nCreating surface geometry for grainDEMSystem . . . "<<END_REPORT;
 	geometry_ = geometry::create(Control(), Property());
 
+	REPORT(0)<<"Reading shape dictionary ..."<<END_REPORT;
+	grains_ = makeUnique<grainShape>(
+		pFlow::shapeFile__,
+		&Control().caseSetup(),
+		Property() );	
 
 	REPORT(0)<<"\nReading grain particles . . ."<<END_REPORT;
-	particles_ = makeUnique<grainFluidParticles>(Control(), Property());
+	particles_ = makeUnique<grainFluidParticles>(Control(), grains_());
 
 
 	insertion_ = makeUnique<grainInsertion>( 

DEMSystems/grainDEMSystem/grainDEMSystem.hpp

@@ -46,6 +46,8 @@ protected:
 
 	uniquePtr<geometry> 			geometry_ = nullptr;
 
+	uniquePtr<grainShape> 			grains_   = nullptr;
+
 	uniquePtr<grainFluidParticles>  particles_ = nullptr;
 
 	uniquePtr<grainInsertion> 		insertion_ = nullptr;

DEMSystems/grainDEMSystem/grainFluidParticles.cpp

@@ -35,8 +35,8 @@ void pFlow::grainFluidParticles::checkHostMemory()
 
 pFlow::grainFluidParticles::grainFluidParticles(
     systemControl &control,
-    const property &prop)
+    const grainShape& grains)
-    : grainParticles(control, prop),
+    : grainParticles(control, grains),
       fluidForce_(
           objectFile(
               "fluidForce",
@@ -67,7 +67,7 @@ bool pFlow::grainFluidParticles::beforeIteration()
 
 bool pFlow::grainFluidParticles::iterate() 
 {
-
+	const auto ti = this->TimeInfo();
 	accelerationTimer().start();
 		pFlow::grainFluidParticlesKernels::acceleration(
 			control().g(),
@@ -78,16 +78,16 @@ bool pFlow::grainFluidParticles::iterate()
 			contactTorque().deviceViewAll(),
 			fluidTorque_.deviceViewAll(),
 			pStruct().activePointsMaskDevice(),
-			accelertion().deviceViewAll(),
+			acceleration().deviceViewAll(),
 			rAcceleration().deviceViewAll()
 			);
 	accelerationTimer().end();
 	
 	intCorrectTimer().start();
 	
-		dynPointStruct().correct(this->dt(), accelertion());
+		dynPointStruct().correct(ti.dt());
 	
-		rVelIntegration().correct(this->dt(), rVelocity(), rAcceleration());
+		rVelIntegration().correct(ti.dt(), rVelocity(), rAcceleration());
 	
 	intCorrectTimer().end();
 	

DEMSystems/grainDEMSystem/grainFluidParticles.hpp

@@ -59,7 +59,7 @@ protected:
 public:
 
 	/// construct from systemControl and property 
-	grainFluidParticles(systemControl &control, const property& prop);
+	grainFluidParticles(systemControl &control, const grainShape& grains);
 
 	~grainFluidParticles() override = default;
 	

DEMSystems/sphereDEMSystem/sphereDEMSystem.cpp

@@ -75,12 +75,18 @@ pFlow::sphereDEMSystem::sphereDEMSystem(
 		propertyFile__,
     	Control().caseSetup().path());
 		
+
 	REPORT(0)<< "\nCreating surface geometry for sphereDEMSystem . . . "<<END_REPORT;
 	geometry_ = geometry::create(Control(), Property());
 
+	REPORT(0)<<"Reading shapes dictionary..."<<END_REPORT;
+	spheres_ = makeUnique<sphereShape>(
+		pFlow::shapeFile__, 
+		&Control().caseSetup(), 
+		Property());
 
 	REPORT(0)<<"\nReading sphere particles . . ."<<END_REPORT;
-	particles_ = makeUnique<sphereFluidParticles>(Control(), Property());
+	particles_ = makeUnique<sphereFluidParticles>(Control(), spheres_());
 
 
 	insertion_ = makeUnique<sphereInsertion>( 

DEMSystems/sphereDEMSystem/sphereDEMSystem.hpp

@@ -46,6 +46,8 @@ protected:
 
 	uniquePtr<geometry> 			geometry_ = nullptr;
 
+	uniquePtr<sphereShape> 			spheres_  = nullptr;
+
 	uniquePtr<sphereFluidParticles> particles_ = nullptr;
 
 	uniquePtr<sphereInsertion> 		insertion_ = nullptr;

DEMSystems/sphereDEMSystem/sphereFluidParticles.cpp

@@ -32,8 +32,8 @@ void pFlow::sphereFluidParticles::checkHostMemory()
 
 pFlow::sphereFluidParticles::sphereFluidParticles(
     systemControl &control,
-    const property &prop)
+    const sphereShape& shpShape)
-    : sphereParticles(control, prop),
+    : sphereParticles(control, shpShape),
       fluidForce_(
           objectFile(
               "fluidForce",
@@ -65,7 +65,7 @@ bool pFlow::sphereFluidParticles::beforeIteration()
 
 bool pFlow::sphereFluidParticles::iterate() 
 {
-
+	const auto ti = this->TimeInfo();
 	accelerationTimer().start();
 		pFlow::sphereFluidParticlesKernels::acceleration(
 			control().g(),
@@ -76,16 +76,16 @@ bool pFlow::sphereFluidParticles::iterate()
 			contactTorque().deviceViewAll(),
 			fluidTorque_.deviceViewAll(),
 			pStruct().activePointsMaskDevice(),
-			accelertion().deviceViewAll(),
+			acceleration().deviceViewAll(),
 			rAcceleration().deviceViewAll()
 			);
 	accelerationTimer().end();
 	
 	intCorrectTimer().start();
 	
-		dynPointStruct().correct(this->dt(), accelertion());
+		dynPointStruct().correct(ti.dt());
 	
-		rVelIntegration().correct(this->dt(), rVelocity(), rAcceleration());
+		rVelIntegration().correct(ti.dt(), rVelocity(), rAcceleration());
 	
 	intCorrectTimer().end();
 	

DEMSystems/sphereDEMSystem/sphereFluidParticles.hpp

@@ -66,7 +66,7 @@ protected:
 public:
 
 	/// construct from systemControl and property 
-	sphereFluidParticles(systemControl &control, const property& prop);
+	sphereFluidParticles(systemControl &control, const sphereShape& shpShape);
 	
 	/// before iteration step 
 	bool beforeIteration() override;

README.md

@@ -1,48 +1,76 @@
-<div align ="center">
+<div align="center">
-<img src="doc/phasicFlow_logo_github.png" style="width: 400px;">
+  <img src="doc/phasicFlow_logo_github.png" style="width: 400px;" alt="PhasicFlow Logo">
 </div>
 
+## **PhasicFlow: High-Performance Discrete Element Method Simulations**
 
-**PhasicFlow** is a parallel C++ code for performing DEM simulations. It can run on shared-memory multi-core computational units such as multi-core CPUs or GPUs (for now it works on CUDA-enabled GPUs). The parallelization method mainly relies on loop-level parallelization on a shared-memory computational unit. You can build and run PhasicFlow in serial mode on regular PCs, in parallel mode for multi-core CPUs, or build it for a GPU device to off-load computations to a GPU. In its current statues you can simulate millions of particles (up to 80M particles tested) on a single desktop computer. You can see the [performance tests of PhasicFlow](https://github.com/PhasicFlow/phasicFlow/wiki/Performance-of-phasicFlow) in the wiki page.
+PhasicFlow is a robust, open-source C++ framework designed for the efficient simulation of granular materials using the Discrete Element Method (DEM). Leveraging parallel computing paradigms, PhasicFlow is capable of executing simulations on shared-memory multi-core architectures, including CPUs and NVIDIA GPUs (CUDA-enabled). The core parallelization strategy focuses on loop-level parallelism, enabling significant performance gains on modern hardware. Users can seamlessly transition between serial execution on standard PCs, parallel execution on multi-core CPUs (OpenMP), and accelerated simulations on GPUs. Currently, PhasicFlow supports simulations involving up to 80 million particles on a single desktop workstation. Detailed performance benchmarks are available on the [PhasicFlow Wiki](https://github.com/PhasicFlow/phasicFlow/wiki/Performance-of-phasicFlow).
 
-**MPI** parallelization with dynamic load balancing is under development. With this level of parallelization, PhasicFlow can leverage the computational power of **multi-gpu** workstations or clusters with distributed memory CPUs. 
+**Scalable Parallelism: MPI Integration**
-In summary PhasicFlow can have 6 execution modes:
-1. Serial on a single CPU core,
-2. Parallel on a multi-core computer/node (using OpenMP),
-3. Parallel on an nvidia-GPU (using Cuda),
-4. Parallel on distributed memory workstation (Using MPI)
-5. Parallel on distributed memory workstations with multi-core nodes (using MPI+OpenMP)
-6. Parallel on workstations with multiple GPUs (using MPI+Cuda).
-## How to build?
-You can build PhasicFlow for CPU and GPU executions. The latest release of PhasicFlow is v-0.1. [Here is a complete step-by-step procedure for building phasicFlow-v-0.1.](https://github.com/PhasicFlow/phasicFlow/wiki/How-to-Build-PhasicFlow).
 
-## Online code documentation
+Ongoing development includes the integration of MPI-based parallelization with dynamic load balancing. This enhancement will extend PhasicFlow's capabilities to distributed memory environments, such as multi-GPU workstations and high-performance computing clusters. Upon completion, PhasicFlow will offer six distinct execution modes:
-You can find a full documentation of the code, its features, and other related materials on [online documentation of the code](https://phasicflow.github.io/phasicFlow/)
 
-## How to use PhasicFlow?
+1.  **Serial Execution:** Single-core CPU.
-You can navigate into [tutorials folder](./tutorials) in the phasicFlow folder to see some simulation case setups. If you need more detailed discription, visit our [wiki page tutorials](https://github.com/PhasicFlow/phasicFlow/wiki/Tutorials).  
+2.  **Shared-Memory Parallelism:** Multi-core CPU (OpenMP).
+3.  **GPU Acceleration:** NVIDIA GPU (CUDA).
+4.  **Distributed-Memory Parallelism:** MPI.
+5.  **Hybrid Parallelism:** MPI + OpenMP.
+6.  **Multi-GPU Parallelism:** MPI + CUDA.
 
-## [PhasicFlowPlus](https://github.com/PhasicFlow/PhasicFlowPlus)
+## **Build and Installation**
-PhasicFlowPlus is and extension to PhasicFlow for simulating particle-fluid systems using resolved and unresolved CFD-DEM. [See the repository of this package.](https://github.com/PhasicFlow/PhasicFlowPlus)
 
+PhasicFlow can be compiled for both CPU and GPU execution.
 
-## Supporting packages
+* **Current Development (v-1.0):** Comprehensive build instructions are available [here](https://github.com/PhasicFlow/phasicFlow/wiki/How-to-build-PhasicFlow%E2%80%90v%E2%80%901.0).
-* [Kokkos](https://github.com/kokkos/kokkos) from National Technology & Engineering Solutions of Sandia, LLC (NTESS)
+* **Latest Release (v-0.1):** Detailed build instructions are available [here](https://github.com/PhasicFlow/phasicFlow/wiki/How-to-Build-PhasicFlow).
-* [CLI11 1.8](https://github.com/CLIUtils/CLI11) from University of Cincinnati.
+
+## **Comprehensive Documentation**
+
+In-depth documentation, including code structure, features, and usage guidelines, is accessible via the [online documentation portal](https://phasicflow.github.io/phasicFlow/).
+
+### **Tutorials and Examples**
+
+Practical examples and simulation setups are provided in the [tutorials directory](./tutorials). For detailed explanations and step-by-step guides, please refer to the [tutorial section on the PhasicFlow Wiki](https://github.com/PhasicFlow/phasicFlow/wiki/Tutorials).
+
+## Contributing to PhasicFlow
+We welcome contributions to PhasicFlow! Whether you're a developer or a new user, there are many ways to get involved. Here's how you can help:
+1. Bug Reports
+2. Suggestions for better user experience
+3. Feature request and algorithm improvements
+4. Tutorials, Simulation Case Setups and documentation
+5. Direct Code Contributions
+
+For more details on how you can contribute to PhasicFlow see [this page](https://github.com/PhasicFlow/phasicFlow/wiki/How-to-contribute-to-PhasicFlow). 
+
+## **PhasicFlowPlus: Coupled CFD-DEM Simulations**
+
+PhasicFlowPlus is an extension of PhasicFlow that facilitates the simulation of particle-fluid systems using resolved and unresolved CFD-DEM methods. The repository for PhasicFlowPlus can be found [here](https://github.com/PhasicFlow/PhasicFlowPlus).
+
+## How to cite PhasicFlow?
 
-## How to cite PhasicFlow
 If you are using PhasicFlow in your research or industrial work, cite the following [article](https://www.sciencedirect.com/science/article/pii/S0010465523001662):
+
 ```
-@article{NOROUZI2023108821,
+@article
-title = {PhasicFlow: A parallel, multi-architecture open-source code for DEM simulations},
+{
-journal = {Computer Physics Communications},
+  NOROUZI2023108821,
-volume = {291},
+  title = {PhasicFlow: A parallel, multi-architecture open-source code for DEM simulations},
-pages = {108821},
+  journal = {Computer Physics Communications},
-year = {2023},
+  volume = {291},
-issn = {0010-4655},
+  pages = {108821},
-doi = {https://doi.org/10.1016/j.cpc.2023.108821},
+  year = {2023},
-url = {https://www.sciencedirect.com/science/article/pii/S0010465523001662},
+  issn = {0010-4655},
-author = {H.R. Norouzi},
+  doi = {https://doi.org/10.1016/j.cpc.2023.108821},
-keywords = {Discrete element method, Parallel computing, CUDA, GPU, OpenMP, Granular flow}
+  url = {https://www.sciencedirect.com/science/article/pii/S0010465523001662},
+  author = {H.R. Norouzi},
+  keywords = {Discrete element method, Parallel computing, CUDA, GPU, OpenMP, Granular flow}
 }
 ```
+
+
+## **Dependencies**
+
+PhasicFlow relies on the following external libraries:
+
+* **Kokkos:** A community-led performance portability ecosystem within the Linux Foundation's High-Performance Software Foundation (HPSF). ([https://github.com/kokkos/kokkos](https://github.com/kokkos/kokkos))
+* **CLI11 1.8:** A command-line interface parser developed by the University of Cincinnati. ([https://github.com/CLIUtils/CLI11](https://github.com/CLIUtils/CLI11))

benchmarks/rotatingDrum/rotatingDrum_1mParticles/caseSetup/interaction

@@ -0,0 +1,60 @@
+/* -------------------------------*- C++ -*--------------------------------- *\ 
+|  phasicFlow File                                                            | 
+|  copyright: www.cemf.ir                                                     | 
+\* ------------------------------------------------------------------------- */ 
+
+objectName     interaction;
+objectType     dictionary;
+fileFormat     ASCII;
+
+/*---------------------------------------------------------------------------*/
+
+materials   (glassMat wallMat); // a list of materials names
+
+densities   (2500.0 2500);      // density of materials [kg/m3]
+
+contactListType   sortedContactList;
+
+contactSearch
+{
+    method         NBS;
+
+    updateInterval   20;
+
+    sizeRatio       1.1;
+
+    cellExtent     0.55;
+
+    adjustableBox   Yes;
+}
+
+model
+{
+    contactForceModel    nonLinearLimited;
+
+    rollingFrictionModel normal;
+    
+    /*
+       Property (glassMat-glassMat      glassMat-wallMat
+                                        wallMat-wallMat);
+    */
+
+    Yeff (1.0e6  1.0e6  
+                 1.0e6);                                    // Young modulus [Pa]
+   
+    Geff (0.8e6  0.8e6  
+                 0.8e6);                                    // Shear modulus [Pa]
+
+    nu   (0.25   0.25  
+                 0.25);                                     // Poisson's ratio [-]
+
+    en   (0.97   0.85   
+                 1.00);                                     // coefficient of normal restitution
+
+    mu   (0.65   0.65   
+                 0.65);                                     // dynamic friction 
+
+    mur  (0.1    0.1    
+                 0.1);                                      // rolling friction 
+}
+

benchmarks/rotatingDrum/rotatingDrum_1mParticles/caseSetup/shapes

@@ -2,12 +2,14 @@
 |  phasicFlow File                                                            | 
 |  copyright: www.cemf.ir                                                     | 
 \* ------------------------------------------------------------------------- */ 
-objectName particleInsertion;
+
-objectType dicrionary;
+objectName     shapes;
+objectType     dictionary;
 fileFormat     ASCII;
 /*---------------------------------------------------------------------------*/
-active         No;                               // is insertion active -> Yes or No
 
-collisionCheck No;                               // is checked          -> Yes or No
+names          (glassBead);    // names of shapes 
 
+diameters      (0.004);        // diameter of shapes 
 
+materials      (glassMat);     // material names for shapes 

benchmarks/rotatingDrum/rotatingDrum_1mParticles/settings/domainDict

@@ -0,0 +1,50 @@
+/* -------------------------------*- C++ -*--------------------------------- *\ 
+|  phasicFlow File                                                            | 
+|  copyright: www.cemf.ir                                                     | 
+\* ------------------------------------------------------------------------- */  
+objectName     domainDict;
+objectType     dictionary;
+fileFormat     ASCII;
+/*---------------------------------------------------------------------------*/
+globalBox                                        // Simulation domain: every particles that goes outside this domain will be deleted
+{
+    min (-0.2 -0.2  0.0);
+    max ( 0.2  0.2  0.8);
+}
+
+boundaries
+{
+    neighborListUpdateInterval 200;
+    
+    updateInterval              20;
+    
+    left
+    {
+        type    exit;                            // other options: periodic, reflective 
+    }
+
+    right
+    {
+        type     exit;                           // other options: periodic, reflective 
+    }
+
+    bottom
+    {
+        type     exit;                           // other options: periodic, reflective 
+    }
+
+    top
+    {
+        type     exit;                           // other options: periodic, reflective 
+    }
+
+    rear
+    {
+        type     exit;                           // other options: periodic, reflective 
+    }
+
+    front
+    {
+        type     exit;                           // other options: periodic, reflective 
+    }
+}

benchmarks/rotatingDrum/rotatingDrum_1mParticles/settings/geometryDict

@@ -0,0 +1,86 @@
+/* -------------------------------*- C++ -*--------------------------------- *\ 
+|  phasicFlow File                                                            | 
+|  copyright: www.cemf.ir                                                     | 
+\* ------------------------------------------------------------------------- */  
+
+objectName     geometryDict;
+objectType     dictionary;
+fileFormat     ASCII;
+
+motionModel rotatingAxis;         // motion model: rotating object around an axis
+
+rotatingAxisInfo                  // information for rotatingAxisMotion motion model 
+{
+    rotAxis 
+    {
+        p1    (0.0 0.0 0.0);      // first point for the axis of rotation    
+        
+        p2    (0.0 0.0 1.0);      // second point for the axis of rotation    
+        
+        omega 1.256;              // rotation speed (rad/s) => 12 rpm
+    }
+}
+
+surfaces
+{
+    cylinder
+    {
+        type        cylinderWall;          // type of the wall
+        
+        p1          (0.0 0.0 0.0);         // begin point of cylinder axis
+           
+        p2          (0.0 0.0 0.8);         // end point of cylinder axis
+        
+        radius1     0.2;                   // radius at p1
+        
+        radius2     0.2;                   // radius at p2
+        
+        resolution  24;                    // number of divisions
+        
+        material    wallMat;               // material name of this wall
+        
+        motion      rotAxis;               // motion component name 
+    }
+
+    /*
+        This is a plane wall at the rear end of cylinder
+    */
+    
+    wall1
+    {
+        type      planeWall;               // type of the wall  
+        
+        p1        (-0.2 -0.2 0.0);         // first point of the wall
+        
+        p2        ( 0.2 -0.2 0.0);         // second point
+        
+        p3        ( 0.2  0.2 0.0);         // third point
+        
+        p4        (-0.2  0.2 0.0);         // fourth point 
+        
+        material  wallMat;                 // material name of the wall  
+        
+        motion    rotAxis;                 // motion component name 
+    }
+
+    /*
+        This is a plane wall at the front end of cylinder
+    */
+    
+    wall2
+    {
+        type      planeWall;               // type of the wall
+        
+        p1        (-0.2 -0.2 0.8);         // first point of the wall
+        
+        p2        ( 0.2 -0.2 0.8);         // second point
+        
+        p3        ( 0.2  0.2 0.8);         // third point
+        
+        p4        (-0.2  0.2 0.8);         // fourth point 
+        
+        material  wallMat;                 // material name of the wall  
+        
+        motion    rotAxis;                 // motion component name 
+    }
+}

benchmarks/rotatingDrum/rotatingDrum_1mParticles/settings/particlesDict

@@ -0,0 +1,47 @@
+/* -------------------------------*- C++ -*--------------------------------- *\ 
+|  phasicFlow File                                                            | 
+|  copyright: www.cemf.ir                                                     | 
+\* ------------------------------------------------------------------------- */  
+
+objectName     particlesDict;
+objectType     dictionary;
+fileFormat     ASCII;
+
+setFields
+{
+    defaultValue 
+    {
+        velocity         realx3     (0 0 0);     // linear velocity (m/s)
+        acceleration     realx3     (0 0 0);     // linear acceleration (m/s2)
+        rotVelocity      realx3     (0 0 0);     // rotational velocity (rad/s)
+        shapeName        word       glassBead;   // name of the particle shape 
+    }
+
+    selectors
+    {}
+}
+
+positionParticles
+{
+    method ordered;
+
+    orderedInfo
+    {
+        diameter   0.004;           // minimum space between centers of particles
+        
+        numPoints  1000000;          // number of particles in the simulation 
+        
+        axisOrder  (z x y);         // axis order for filling the space with particles
+    }
+
+    regionType cylinder;            // other options: box and sphere    
+
+    cylinderInfo                    // cylinder for positioning particles 
+    {
+        p1     (0.0  0.0  0.01);    // lower corner point of the box 
+        
+        p2     (0.0  0.0  0.79);    // upper corner point of the box 
+        
+        radius 0.195;               // radius of cylinder 
+    }
+}

benchmarks/rotatingDrum/rotatingDrum_1mParticles/settings/settingsDict

@@ -0,0 +1,34 @@
+/* -------------------------------*- C++ -*--------------------------------- *\ 
+|  phasicFlow File                                                            | 
+|  copyright: www.cemf.ir                                                     | 
+\* ------------------------------------------------------------------------- */  
+objectName    settingsDict;
+objectType    dictionary;
+fileFormat    ASCII;
+/*---------------------------------------------------------------------------*/
+run             rotatingDrum_1mParticles;
+
+dt              0.00001;      // time step for integration (s)
+
+startTime       0;            // start time for simulation 
+
+endTime         4;           // end time for simulation 
+
+saveInterval    0.2;          // time interval for saving the simulation
+
+timePrecision   5;            // maximum number of digits for time folder 
+
+g               (0 -9.8 0);   // gravity vector (m/s2) 
+
+includeObjects (diameter);                     // save necessary (i.e., required) data on disk
+
+// exclude unnecessary data from saving on disk
+excludeObjects (rVelocity.dy1 pStructPosition.dy1 pStructVelocity.dy1); 
+
+integrationMethod         AdamsBashforth2;     // integration method 
+
+writeFormat               binary;              // data writting format (ascii or binary)
+
+timersReport              Yes;      
+
+timersReportInterval      0.01;    

benchmarks/rotatingDrum/rotatingDrum_250kParticles/caseSetup/interaction

@@ -0,0 +1,60 @@
+/* -------------------------------*- C++ -*--------------------------------- *\ 
+|  phasicFlow File                                                            | 
+|  copyright: www.cemf.ir                                                     | 
+\* ------------------------------------------------------------------------- */ 
+
+objectName     interaction;
+objectType     dictionary;
+fileFormat     ASCII;
+
+/*---------------------------------------------------------------------------*/
+
+materials   (glassMat wallMat); // a list of materials names
+
+densities   (2500.0 2500);      // density of materials [kg/m3]
+
+contactListType   sortedContactList;
+
+contactSearch
+{
+    method         NBS;
+
+    updateInterval   20;
+
+    sizeRatio       1.1;
+
+    cellExtent     0.55;
+
+    adjustableBox   Yes;
+}
+
+model
+{
+    contactForceModel    nonLinearLimited;
+
+    rollingFrictionModel normal;
+    
+    /*
+       Property (glassMat-glassMat      glassMat-wallMat
+                                        wallMat-wallMat);
+    */
+
+    Yeff (1.0e6  1.0e6  
+                 1.0e6);                                    // Young modulus [Pa]
+   
+    Geff (0.8e6  0.8e6  
+                 0.8e6);                                    // Shear modulus [Pa]
+
+    nu   (0.25   0.25  
+                 0.25);                                     // Poisson's ratio [-]
+
+    en   (0.97   0.85   
+                 1.00);                                     // coefficient of normal restitution
+
+    mu   (0.65   0.65   
+                 0.65);                                     // dynamic friction 
+
+    mur  (0.1    0.1    
+                 0.1);                                      // rolling friction 
+}
+

benchmarks/rotatingDrum/rotatingDrum_250kParticles/caseSetup/shapes

@@ -2,12 +2,14 @@
 |  phasicFlow File                                                            | 
 |  copyright: www.cemf.ir                                                     | 
 \* ------------------------------------------------------------------------- */ 
-objectName particleInsertion;
+
-objectType dicrionary;
+objectName     shapes;
+objectType     dictionary;
 fileFormat     ASCII;
 /*---------------------------------------------------------------------------*/
-active         No;                               // is insertion active -> Yes or No
 
-collisionCheck No;                               // is checked          -> Yes or No
+names          (glassBead);    // names of shapes 
 
+diameters      (0.006);        // diameter of shapes 
 
+materials      (glassMat);     // material names for shapes 

benchmarks/rotatingDrum/rotatingDrum_250kParticles/settings/domainDict

@@ -0,0 +1,50 @@
+/* -------------------------------*- C++ -*--------------------------------- *\ 
+|  phasicFlow File                                                            | 
+|  copyright: www.cemf.ir                                                     | 
+\* ------------------------------------------------------------------------- */  
+objectName     domainDict;
+objectType     dictionary;
+fileFormat     ASCII;
+/*---------------------------------------------------------------------------*/
+globalBox                                        // Simulation domain: every particles that goes outside this domain will be deleted
+{
+    min (-0.2 -0.2  0.0);
+    max ( 0.2  0.2  0.8);
+}
+
+boundaries
+{
+    neighborListUpdateInterval 200;
+    
+    updateInterval              20;
+    
+    left
+    {
+        type    exit;                            // other options: periodic, reflective 
+    }
+
+    right
+    {
+        type     exit;                           // other options: periodic, reflective 
+    }
+
+    bottom
+    {
+        type     exit;                           // other options: periodic, reflective 
+    }
+
+    top
+    {
+        type     exit;                           // other options: periodic, reflective 
+    }
+
+    rear
+    {
+        type     exit;                           // other options: periodic, reflective 
+    }
+
+    front
+    {
+        type     exit;                           // other options: periodic, reflective 
+    }
+}

benchmarks/rotatingDrum/rotatingDrum_250kParticles/settings/geometryDict

@@ -0,0 +1,86 @@
+/* -------------------------------*- C++ -*--------------------------------- *\ 
+|  phasicFlow File                                                            | 
+|  copyright: www.cemf.ir                                                     | 
+\* ------------------------------------------------------------------------- */  
+
+objectName     geometryDict;
+objectType     dictionary;
+fileFormat     ASCII;
+
+motionModel rotatingAxis;         // motion model: rotating object around an axis
+
+rotatingAxisInfo                  // information for rotatingAxisMotion motion model 
+{
+    rotAxis 
+    {
+        p1    (0.0 0.0 0.0);      // first point for the axis of rotation    
+        
+        p2    (0.0 0.0 1.0);      // second point for the axis of rotation    
+        
+        omega 1.256;              // rotation speed (rad/s) => 12 rpm
+    }
+}
+
+surfaces
+{
+    cylinder
+    {
+        type        cylinderWall;          // type of the wall
+        
+        p1          (0.0 0.0 0.0);         // begin point of cylinder axis
+           
+        p2          (0.0 0.0 0.8);         // end point of cylinder axis
+        
+        radius1     0.2;                   // radius at p1
+        
+        radius2     0.2;                   // radius at p2
+        
+        resolution  24;                    // number of divisions
+        
+        material    wallMat;               // material name of this wall
+        
+        motion      rotAxis;               // motion component name 
+    }
+
+    /*
+        This is a plane wall at the rear end of cylinder
+    */
+    
+    wall1
+    {
+        type      planeWall;               // type of the wall  
+        
+        p1        (-0.2 -0.2 0.0);         // first point of the wall
+        
+        p2        ( 0.2 -0.2 0.0);         // second point
+        
+        p3        ( 0.2  0.2 0.0);         // third point
+        
+        p4        (-0.2  0.2 0.0);         // fourth point 
+        
+        material  wallMat;                 // material name of the wall  
+        
+        motion    rotAxis;                 // motion component name 
+    }
+
+    /*
+        This is a plane wall at the front end of cylinder
+    */
+    
+    wall2
+    {
+        type      planeWall;               // type of the wall
+        
+        p1        (-0.2 -0.2 0.8);         // first point of the wall
+        
+        p2        ( 0.2 -0.2 0.8);         // second point
+        
+        p3        ( 0.2  0.2 0.8);         // third point
+        
+        p4        (-0.2  0.2 0.8);         // fourth point 
+        
+        material  wallMat;                 // material name of the wall  
+        
+        motion    rotAxis;                 // motion component name 
+    }
+}

benchmarks/rotatingDrum/rotatingDrum_250kParticles/settings/particlesDict

@@ -0,0 +1,47 @@
+/* -------------------------------*- C++ -*--------------------------------- *\ 
+|  phasicFlow File                                                            | 
+|  copyright: www.cemf.ir                                                     | 
+\* ------------------------------------------------------------------------- */  
+
+objectName     particlesDict;
+objectType     dictionary;
+fileFormat     ASCII;
+
+setFields
+{
+    defaultValue 
+    {
+        velocity         realx3     (0 0 0);     // linear velocity (m/s)
+        acceleration     realx3     (0 0 0);     // linear acceleration (m/s2)
+        rotVelocity      realx3     (0 0 0);     // rotational velocity (rad/s)
+        shapeName        word       glassBead;   // name of the particle shape 
+    }
+
+    selectors
+    {}
+}
+
+positionParticles
+{
+    method ordered;
+
+    orderedInfo
+    {
+        diameter   0.006;           // minimum space between centers of particles
+        
+        numPoints  250000;          // number of particles in the simulation 
+        
+        axisOrder  (z x y);         // axis order for filling the space with particles
+    }
+
+    regionType cylinder;            // other options: box and sphere    
+
+    cylinderInfo                    // cylinder for positioning particles 
+    {
+        p1     (0.0  0.0  0.01);    // lower corner point of the box 
+        
+        p2     (0.0  0.0  0.79);    // upper corner point of the box 
+        
+        radius 0.195;               // radius of cylinder 
+    }
+}

benchmarks/rotatingDrum/rotatingDrum_250kParticles/settings/settingsDict

@@ -0,0 +1,34 @@
+/* -------------------------------*- C++ -*--------------------------------- *\ 
+|  phasicFlow File                                                            | 
+|  copyright: www.cemf.ir                                                     | 
+\* ------------------------------------------------------------------------- */  
+objectName    settingsDict;
+objectType    dictionary;
+fileFormat    ASCII;
+/*---------------------------------------------------------------------------*/
+run             rotatingDrum_250KParticles;
+
+dt              0.00001;      // time step for integration (s)
+
+startTime       0;            // start time for simulation 
+
+endTime         4;           // end time for simulation 
+
+saveInterval    0.2;          // time interval for saving the simulation
+
+timePrecision   5;            // maximum number of digits for time folder 
+
+g               (0 -9.8 0);   // gravity vector (m/s2) 
+
+includeObjects (diameter);                     // save necessary (i.e., required) data on disk
+
+// exclude unnecessary data from saving on disk
+excludeObjects (rVelocity.dy1 pStructPosition.dy1 pStructVelocity.dy1); 
+
+integrationMethod         AdamsBashforth2;     // integration method 
+
+writeFormat               binary;              // data writting format (ascii or binary)
+
+timersReport              Yes;      
+
+timersReportInterval      0.01;    

benchmarks/rotatingDrum/rotatingDrum_2mParticles/caseSetup/interaction

@@ -0,0 +1,60 @@
+/* -------------------------------*- C++ -*--------------------------------- *\ 
+|  phasicFlow File                                                            | 
+|  copyright: www.cemf.ir                                                     | 
+\* ------------------------------------------------------------------------- */ 
+
+objectName     interaction;
+objectType     dictionary;
+fileFormat     ASCII;
+
+/*---------------------------------------------------------------------------*/
+
+materials   (glassMat wallMat); // a list of materials names
+
+densities   (2500.0 2500);      // density of materials [kg/m3]
+
+contactListType   sortedContactList;
+
+contactSearch
+{
+    method         NBS;
+
+    updateInterval   20;
+
+    sizeRatio       1.1;
+
+    cellExtent     0.55;
+
+    adjustableBox   Yes;
+}
+
+model
+{
+    contactForceModel    nonLinearLimited;
+
+    rollingFrictionModel normal;
+    
+    /*
+       Property (glassMat-glassMat      glassMat-wallMat
+                                        wallMat-wallMat);
+    */
+
+    Yeff (1.0e6  1.0e6  
+                 1.0e6);                                    // Young modulus [Pa]
+   
+    Geff (0.8e6  0.8e6  
+                 0.8e6);                                    // Shear modulus [Pa]
+
+    nu   (0.25   0.25  
+                 0.25);                                     // Poisson's ratio [-]
+
+    en   (0.97   0.85   
+                 1.00);                                     // coefficient of normal restitution
+
+    mu   (0.65   0.65   
+                 0.65);                                     // dynamic friction 
+
+    mur  (0.1    0.1    
+                 0.1);                                      // rolling friction 
+}
+

benchmarks/rotatingDrum/rotatingDrum_2mParticles/caseSetup/shapes

@@ -2,11 +2,14 @@
 |  phasicFlow File                                                            | 
 |  copyright: www.cemf.ir                                                     | 
 \* ------------------------------------------------------------------------- */ 
-objectName 	sphereDict;
+
-objectType 	sphereShape;
+objectName     shapes;
+objectType     dictionary;
 fileFormat     ASCII;
 /*---------------------------------------------------------------------------*/
 
-names 		(sphere1); 	// names of shapes 
+names          (glassBead);    // names of shapes 
-diameters 	(0.004);	// diameter of shapes 
+
-materials	(prop1);	// material names for shapes 
+diameters      (0.003);        // diameter of shapes 
+
+materials      (glassMat);     // material names for shapes 

benchmarks/rotatingDrum/rotatingDrum_2mParticles/cleanThisCase

@@ -0,0 +1,7 @@
+#!/bin/sh
+cd ${0%/*} || exit 1    # Run from this directory
+
+ls | grep -P "^(([0-9]+\.?[0-9]*)|(\.[0-9]+))$" | xargs -d"\n" rm -rf
+rm -rf VTK 
+
+#------------------------------------------------------------------------------

benchmarks/rotatingDrum/rotatingDrum_2mParticles/runThisCase

@@ -0,0 +1,21 @@
+#!/bin/sh
+cd ${0%/*} || exit 1    # Run from this directory
+echo "\n<--------------------------------------------------------------------->"
+echo "1) Creating particles"
+echo "<--------------------------------------------------------------------->\n"
+particlesPhasicFlow
+
+echo "\n<--------------------------------------------------------------------->"
+echo "2) Creating geometry"
+echo "<--------------------------------------------------------------------->\n"
+geometryPhasicFlow
+
+echo "\n<--------------------------------------------------------------------->"
+echo "3) Running the case"
+echo "<--------------------------------------------------------------------->\n"
+sphereGranFlow
+
+
+
+
+#------------------------------------------------------------------------------

benchmarks/rotatingDrum/rotatingDrum_2mParticles/settings/domainDict

@@ -0,0 +1,50 @@
+/* -------------------------------*- C++ -*--------------------------------- *\ 
+|  phasicFlow File                                                            | 
+|  copyright: www.cemf.ir                                                     | 
+\* ------------------------------------------------------------------------- */  
+objectName     domainDict;
+objectType     dictionary;
+fileFormat     ASCII;
+/*---------------------------------------------------------------------------*/
+globalBox                                        // Simulation domain: every particles that goes outside this domain will be deleted
+{
+    min (-0.2 -0.2  0.0);
+    max ( 0.2  0.2  1.2);
+}
+
+boundaries
+{
+    neighborListUpdateInterval 200;
+    
+    updateInterval              20;
+    
+    left
+    {
+        type    exit;                            // other options: periodic, reflective 
+    }
+
+    right
+    {
+        type     exit;                           // other options: periodic, reflective 
+    }
+
+    bottom
+    {
+        type     exit;                           // other options: periodic, reflective 
+    }
+
+    top
+    {
+        type     exit;                           // other options: periodic, reflective 
+    }
+
+    rear
+    {
+        type     exit;                           // other options: periodic, reflective 
+    }
+
+    front
+    {
+        type     exit;                           // other options: periodic, reflective 
+    }
+}

benchmarks/rotatingDrum/rotatingDrum_2mParticles/settings/geometryDict

@@ -0,0 +1,86 @@
+/* -------------------------------*- C++ -*--------------------------------- *\ 
+|  phasicFlow File                                                            | 
+|  copyright: www.cemf.ir                                                     | 
+\* ------------------------------------------------------------------------- */  
+
+objectName     geometryDict;
+objectType     dictionary;
+fileFormat     ASCII;
+
+motionModel rotatingAxis;         // motion model: rotating object around an axis
+
+rotatingAxisInfo                  // information for rotatingAxisMotion motion model 
+{
+    rotAxis 
+    {
+        p1    (0.0 0.0 0.0);      // first point for the axis of rotation    
+        
+        p2    (0.0 0.0 1.0);      // second point for the axis of rotation    
+        
+        omega 1.256;              // rotation speed (rad/s) => 12 rpm
+    }
+}
+
+surfaces
+{
+    cylinder
+    {
+        type        cylinderWall;          // type of the wall
+        
+        p1          (0.0 0.0 0.0);         // begin point of cylinder axis
+           
+        p2          (0.0 0.0 1.2);         // end point of cylinder axis
+        
+        radius1     0.2;                   // radius at p1
+        
+        radius2     0.2;                   // radius at p2
+        
+        resolution  24;                    // number of divisions
+        
+        material    wallMat;               // material name of this wall
+        
+        motion      rotAxis;               // motion component name 
+    }
+
+    /*
+        This is a plane wall at the rear end of cylinder
+    */
+    
+    wall1
+    {
+        type      planeWall;               // type of the wall  
+        
+        p1        (-0.2 -0.2 0.0);         // first point of the wall
+        
+        p2        ( 0.2 -0.2 0.0);         // second point
+        
+        p3        ( 0.2  0.2 0.0);         // third point
+        
+        p4        (-0.2  0.2 0.0);         // fourth point 
+        
+        material  wallMat;                 // material name of the wall  
+        
+        motion    rotAxis;                 // motion component name 
+    }
+
+    /*
+        This is a plane wall at the front end of cylinder
+    */
+    
+    wall2
+    {
+        type      planeWall;               // type of the wall
+        
+        p1        (-0.2 -0.2 1.2);         // first point of the wall
+        
+        p2        ( 0.2 -0.2 1.2);         // second point
+        
+        p3        ( 0.2  0.2 1.2);         // third point
+        
+        p4        (-0.2  0.2 1.2);         // fourth point 
+        
+        material  wallMat;                 // material name of the wall  
+        
+        motion    rotAxis;                 // motion component name 
+    }
+}

benchmarks/rotatingDrum/rotatingDrum_2mParticles/settings/particlesDict

@@ -0,0 +1,47 @@
+/* -------------------------------*- C++ -*--------------------------------- *\ 
+|  phasicFlow File                                                            | 
+|  copyright: www.cemf.ir                                                     | 
+\* ------------------------------------------------------------------------- */  
+
+objectName     particlesDict;
+objectType     dictionary;
+fileFormat     ASCII;
+
+setFields
+{
+    defaultValue 
+    {
+        velocity         realx3     (0 0 0);     // linear velocity (m/s)
+        acceleration     realx3     (0 0 0);     // linear acceleration (m/s2)
+        rotVelocity      realx3     (0 0 0);     // rotational velocity (rad/s)
+        shapeName        word       glassBead;   // name of the particle shape 
+    }
+
+    selectors
+    {}
+}
+
+positionParticles
+{
+    method ordered;
+
+    orderedInfo
+    {
+        diameter   0.003;           // minimum space between centers of particles
+        
+        numPoints  2000000;          // number of particles in the simulation 
+        
+        axisOrder  (z x y);         // axis order for filling the space with particles
+    }
+
+    regionType cylinder;            // other options: box and sphere    
+
+    cylinderInfo                    // cylinder for positioning particles 
+    {
+        p1     (0.0  0.0  0.01);    // lower corner point of the box 
+        
+        p2     (0.0  0.0  1.19);    // upper corner point of the box 
+        
+        radius 0.195;               // radius of cylinder 
+    }
+}

benchmarks/rotatingDrum/rotatingDrum_2mParticles/settings/settingsDict

@@ -0,0 +1,34 @@
+/* -------------------------------*- C++ -*--------------------------------- *\ 
+|  phasicFlow File                                                            | 
+|  copyright: www.cemf.ir                                                     | 
+\* ------------------------------------------------------------------------- */  
+objectName    settingsDict;
+objectType    dictionary;
+fileFormat    ASCII;
+/*---------------------------------------------------------------------------*/
+run             rotatingDrum_2mParticles;
+
+dt              0.00001;      // time step for integration (s)
+
+startTime       0;            // start time for simulation 
+
+endTime         4;           // end time for simulation 
+
+saveInterval    0.2;          // time interval for saving the simulation
+
+timePrecision   5;            // maximum number of digits for time folder 
+
+g               (0 -9.8 0);   // gravity vector (m/s2) 
+
+includeObjects (diameter);                     // save necessary (i.e., required) data on disk
+
+// exclude unnecessary data from saving on disk
+excludeObjects (rVelocity.dy1 pStructPosition.dy1 pStructVelocity.dy1); 
+
+integrationMethod         AdamsBashforth2;     // integration method 
+
+writeFormat               binary;              // data writting format (ascii or binary)
+
+timersReport              Yes;      
+
+timersReportInterval      0.01;    

benchmarks/rotatingDrum/rotatingDrum_4mParticles/caseSetup/interaction

@@ -0,0 +1,60 @@
+/* -------------------------------*- C++ -*--------------------------------- *\ 
+|  phasicFlow File                                                            | 
+|  copyright: www.cemf.ir                                                     | 
+\* ------------------------------------------------------------------------- */ 
+
+objectName     interaction;
+objectType     dictionary;
+fileFormat     ASCII;
+
+/*---------------------------------------------------------------------------*/
+
+materials   (glassMat wallMat); // a list of materials names
+
+densities   (2500.0 2500);      // density of materials [kg/m3]
+
+contactListType   sortedContactList;
+
+contactSearch
+{
+    method         NBS;
+
+    updateInterval   20;
+
+    sizeRatio       1.1;
+
+    cellExtent     0.55;
+
+    adjustableBox   Yes;
+}
+
+model
+{
+    contactForceModel    nonLinearLimited;
+
+    rollingFrictionModel normal;
+    
+    /*
+       Property (glassMat-glassMat      glassMat-wallMat
+                                        wallMat-wallMat);
+    */
+
+    Yeff (1.0e6  1.0e6  
+                 1.0e6);                                    // Young modulus [Pa]
+   
+    Geff (0.8e6  0.8e6  
+                 0.8e6);                                    // Shear modulus [Pa]
+
+    nu   (0.25   0.25  
+                 0.25);                                     // Poisson's ratio [-]
+
+    en   (0.97   0.85   
+                 1.00);                                     // coefficient of normal restitution
+
+    mu   (0.65   0.65   
+                 0.65);                                     // dynamic friction 
+
+    mur  (0.1    0.1    
+                 0.1);                                      // rolling friction 
+}
+

benchmarks/rotatingDrum/rotatingDrum_4mParticles/caseSetup/shapes

@@ -0,0 +1,15 @@
+/* -------------------------------*- C++ -*--------------------------------- *\ 
+|  phasicFlow File                                                            | 
+|  copyright: www.cemf.ir                                                     | 
+\* ------------------------------------------------------------------------- */ 
+
+objectName     shapes;
+objectType     dictionary;
+fileFormat     ASCII;
+/*---------------------------------------------------------------------------*/
+
+names          (glassBead);    // names of shapes 
+
+diameters      (0.003);        // diameter of shapes 
+
+materials      (glassMat);     // material names for shapes 

benchmarks/rotatingDrum/rotatingDrum_4mParticles/cleanThisCase

@@ -0,0 +1,7 @@
+#!/bin/sh
+cd ${0%/*} || exit 1    # Run from this directory
+
+ls | grep -P "^(([0-9]+\.?[0-9]*)|(\.[0-9]+))$" | xargs -d"\n" rm -rf
+rm -rf VTK 
+
+#------------------------------------------------------------------------------

benchmarks/rotatingDrum/rotatingDrum_4mParticles/runThisCase

@@ -0,0 +1,21 @@
+#!/bin/sh
+cd ${0%/*} || exit 1    # Run from this directory
+echo "\n<--------------------------------------------------------------------->"
+echo "1) Creating particles"
+echo "<--------------------------------------------------------------------->\n"
+particlesPhasicFlow
+
+echo "\n<--------------------------------------------------------------------->"
+echo "2) Creating geometry"
+echo "<--------------------------------------------------------------------->\n"
+geometryPhasicFlow
+
+echo "\n<--------------------------------------------------------------------->"
+echo "3) Running the case"
+echo "<--------------------------------------------------------------------->\n"
+sphereGranFlow
+
+
+
+
+#------------------------------------------------------------------------------

benchmarks/rotatingDrum/rotatingDrum_4mParticles/settings/domainDict

@@ -0,0 +1,50 @@
+/* -------------------------------*- C++ -*--------------------------------- *\ 
+|  phasicFlow File                                                            | 
+|  copyright: www.cemf.ir                                                     | 
+\* ------------------------------------------------------------------------- */  
+objectName     domainDict;
+objectType     dictionary;
+fileFormat     ASCII;
+/*---------------------------------------------------------------------------*/
+globalBox                                        // Simulation domain: every particles that goes outside this domain will be deleted
+{
+    min (-0.2 -0.2  0.0);
+    max ( 0.2  0.2  1.6);
+}
+
+boundaries
+{
+    neighborListUpdateInterval 200;
+    
+    updateInterval              20;
+    
+    left
+    {
+        type    exit;                            // other options: periodic, reflective 
+    }
+
+    right
+    {
+        type     exit;                           // other options: periodic, reflective 
+    }
+
+    bottom
+    {
+        type     exit;                           // other options: periodic, reflective 
+    }
+
+    top
+    {
+        type     exit;                           // other options: periodic, reflective 
+    }
+
+    rear
+    {
+        type     exit;                           // other options: periodic, reflective 
+    }
+
+    front
+    {
+        type     exit;                           // other options: periodic, reflective 
+    }
+}

benchmarks/rotatingDrum/rotatingDrum_4mParticles/settings/geometryDict

@@ -0,0 +1,86 @@
+/* -------------------------------*- C++ -*--------------------------------- *\ 
+|  phasicFlow File                                                            | 
+|  copyright: www.cemf.ir                                                     | 
+\* ------------------------------------------------------------------------- */  
+
+objectName     geometryDict;
+objectType     dictionary;
+fileFormat     ASCII;
+
+motionModel rotatingAxis;         // motion model: rotating object around an axis
+
+rotatingAxisInfo                  // information for rotatingAxisMotion motion model 
+{
+    rotAxis 
+    {
+        p1    (0.0 0.0 0.0);      // first point for the axis of rotation    
+        
+        p2    (0.0 0.0 1.0);      // second point for the axis of rotation    
+        
+        omega 1.256;              // rotation speed (rad/s) => 12 rpm
+    }
+}
+
+surfaces
+{
+    cylinder
+    {
+        type        cylinderWall;          // type of the wall
+        
+        p1          (0.0 0.0 0.0);         // begin point of cylinder axis
+           
+        p2          (0.0 0.0 1.6);         // end point of cylinder axis
+        
+        radius1     0.2;                   // radius at p1
+        
+        radius2     0.2;                   // radius at p2
+        
+        resolution  24;                    // number of divisions
+        
+        material    wallMat;               // material name of this wall
+        
+        motion      rotAxis;               // motion component name 
+    }
+
+    /*
+        This is a plane wall at the rear end of cylinder
+    */
+    
+    wall1
+    {
+        type      planeWall;               // type of the wall  
+        
+        p1        (-0.2 -0.2 0.0);         // first point of the wall
+        
+        p2        ( 0.2 -0.2 0.0);         // second point
+        
+        p3        ( 0.2  0.2 0.0);         // third point
+        
+        p4        (-0.2  0.2 0.0);         // fourth point 
+        
+        material  wallMat;                 // material name of the wall  
+        
+        motion    rotAxis;                 // motion component name 
+    }
+
+    /*
+        This is a plane wall at the front end of cylinder
+    */
+    
+    wall2
+    {
+        type      planeWall;               // type of the wall
+        
+        p1        (-0.2 -0.2 1.6);         // first point of the wall
+        
+        p2        ( 0.2 -0.2 1.6);         // second point
+        
+        p3        ( 0.2  0.2 1.6);         // third point
+        
+        p4        (-0.2  0.2 1.6);         // fourth point 
+        
+        material  wallMat;                 // material name of the wall  
+        
+        motion    rotAxis;                 // motion component name 
+    }
+}

benchmarks/rotatingDrum/rotatingDrum_4mParticles/settings/particlesDict

@@ -0,0 +1,47 @@
+/* -------------------------------*- C++ -*--------------------------------- *\ 
+|  phasicFlow File                                                            | 
+|  copyright: www.cemf.ir                                                     | 
+\* ------------------------------------------------------------------------- */  
+
+objectName     particlesDict;
+objectType     dictionary;
+fileFormat     ASCII;
+
+setFields
+{
+    defaultValue 
+    {
+        velocity         realx3     (0 0 0);     // linear velocity (m/s)
+        acceleration     realx3     (0 0 0);     // linear acceleration (m/s2)
+        rotVelocity      realx3     (0 0 0);     // rotational velocity (rad/s)
+        shapeName        word       glassBead;   // name of the particle shape 
+    }
+
+    selectors
+    {}
+}
+
+positionParticles
+{
+    method ordered;
+
+    orderedInfo
+    {
+        diameter   0.003;           // minimum space between centers of particles
+        
+        numPoints  4000000;          // number of particles in the simulation 
+        
+        axisOrder  (z x y);         // axis order for filling the space with particles
+    }
+
+    regionType cylinder;            // other options: box and sphere    
+
+    cylinderInfo                    // cylinder for positioning particles 
+    {
+        p1     (0.0  0.0  0.01);    // lower corner point of the box 
+        
+        p2     (0.0  0.0  1.59);    // upper corner point of the box 
+        
+        radius 0.195;               // radius of cylinder 
+    }
+}

benchmarks/rotatingDrum/rotatingDrum_4mParticles/settings/settingsDict

@@ -0,0 +1,34 @@
+/* -------------------------------*- C++ -*--------------------------------- *\ 
+|  phasicFlow File                                                            | 
+|  copyright: www.cemf.ir                                                     | 
+\* ------------------------------------------------------------------------- */  
+objectName    settingsDict;
+objectType    dictionary;
+fileFormat    ASCII;
+/*---------------------------------------------------------------------------*/
+run             rotatingDrum_4mParticles;
+
+dt              0.00001;      // time step for integration (s)
+
+startTime       0;            // start time for simulation 
+
+endTime         4;           // end time for simulation 
+
+saveInterval    0.2;          // time interval for saving the simulation
+
+timePrecision   5;            // maximum number of digits for time folder 
+
+g               (0 -9.8 0);   // gravity vector (m/s2) 
+
+includeObjects (diameter);                     // save necessary (i.e., required) data on disk
+
+// exclude unnecessary data from saving on disk
+excludeObjects (rVelocity.dy1 pStructPosition.dy1 pStructVelocity.dy1); 
+
+integrationMethod         AdamsBashforth2;     // integration method 
+
+writeFormat               binary;              // data writting format (ascii or binary)
+
+timersReport              Yes;      
+
+timersReportInterval      0.01;    

benchmarks/rotatingDrum/rotatingDrum_500kParticles/caseSetup/interaction

@@ -0,0 +1,60 @@
+/* -------------------------------*- C++ -*--------------------------------- *\ 
+|  phasicFlow File                                                            | 
+|  copyright: www.cemf.ir                                                     | 
+\* ------------------------------------------------------------------------- */ 
+
+objectName     interaction;
+objectType     dictionary;
+fileFormat     ASCII;
+
+/*---------------------------------------------------------------------------*/
+
+materials   (glassMat wallMat); // a list of materials names
+
+densities   (2500.0 2500);      // density of materials [kg/m3]
+
+contactListType   sortedContactList;
+
+contactSearch
+{
+    method         NBS;
+
+    updateInterval   20;
+
+    sizeRatio       1.1;
+
+    cellExtent     0.55;
+
+    adjustableBox   Yes;
+}
+
+model
+{
+    contactForceModel    nonLinearLimited;
+
+    rollingFrictionModel normal;
+    
+    /*
+       Property (glassMat-glassMat      glassMat-wallMat
+                                        wallMat-wallMat);
+    */
+
+    Yeff (1.0e6  1.0e6  
+                 1.0e6);                                    // Young modulus [Pa]
+   
+    Geff (0.8e6  0.8e6  
+                 0.8e6);                                    // Shear modulus [Pa]
+
+    nu   (0.25   0.25  
+                 0.25);                                     // Poisson's ratio [-]
+
+    en   (0.97   0.85   
+                 1.00);                                     // coefficient of normal restitution
+
+    mu   (0.65   0.65   
+                 0.65);                                     // dynamic friction 
+
+    mur  (0.1    0.1    
+                 0.1);                                      // rolling friction 
+}
+

benchmarks/rotatingDrum/rotatingDrum_500kParticles/caseSetup/shapes

@@ -0,0 +1,15 @@
+/* -------------------------------*- C++ -*--------------------------------- *\ 
+|  phasicFlow File                                                            | 
+|  copyright: www.cemf.ir                                                     | 
+\* ------------------------------------------------------------------------- */ 
+
+objectName     shapes;
+objectType     dictionary;
+fileFormat     ASCII;
+/*---------------------------------------------------------------------------*/
+
+names          (glassBead);    // names of shapes 
+
+diameters      (0.005);        // diameter of shapes 
+
+materials      (glassMat);     // material names for shapes 

benchmarks/rotatingDrum/rotatingDrum_500kParticles/cleanThisCase

@@ -0,0 +1,7 @@
+#!/bin/sh
+cd ${0%/*} || exit 1    # Run from this directory
+
+ls | grep -P "^(([0-9]+\.?[0-9]*)|(\.[0-9]+))$" | xargs -d"\n" rm -rf
+rm -rf VTK 
+
+#------------------------------------------------------------------------------

benchmarks/rotatingDrum/rotatingDrum_500kParticles/runThisCase

@@ -0,0 +1,21 @@
+#!/bin/sh
+cd ${0%/*} || exit 1    # Run from this directory
+echo "\n<--------------------------------------------------------------------->"
+echo "1) Creating particles"
+echo "<--------------------------------------------------------------------->\n"
+particlesPhasicFlow
+
+echo "\n<--------------------------------------------------------------------->"
+echo "2) Creating geometry"
+echo "<--------------------------------------------------------------------->\n"
+geometryPhasicFlow
+
+echo "\n<--------------------------------------------------------------------->"
+echo "3) Running the case"
+echo "<--------------------------------------------------------------------->\n"
+sphereGranFlow
+
+
+
+
+#------------------------------------------------------------------------------

benchmarks/rotatingDrum/rotatingDrum_500kParticles/settings/domainDict

@@ -0,0 +1,50 @@
+/* -------------------------------*- C++ -*--------------------------------- *\ 
+|  phasicFlow File                                                            | 
+|  copyright: www.cemf.ir                                                     | 
+\* ------------------------------------------------------------------------- */  
+objectName     domainDict;
+objectType     dictionary;
+fileFormat     ASCII;
+/*---------------------------------------------------------------------------*/
+globalBox                                        // Simulation domain: every particles that goes outside this domain will be deleted
+{
+    min (-0.2 -0.2  0.0);
+    max ( 0.2  0.2  0.8);
+}
+
+boundaries
+{
+    neighborListUpdateInterval 200;
+    
+    updateInterval              20;
+    
+    left
+    {
+        type    exit;                            // other options: periodic, reflective 
+    }
+
+    right
+    {
+        type     exit;                           // other options: periodic, reflective 
+    }
+
+    bottom
+    {
+        type     exit;                           // other options: periodic, reflective 
+    }
+
+    top
+    {
+        type     exit;                           // other options: periodic, reflective 
+    }
+
+    rear
+    {
+        type     exit;                           // other options: periodic, reflective 
+    }
+
+    front
+    {
+        type     exit;                           // other options: periodic, reflective 
+    }
+}

benchmarks/rotatingDrum/rotatingDrum_500kParticles/settings/geometryDict

@@ -0,0 +1,86 @@
+/* -------------------------------*- C++ -*--------------------------------- *\ 
+|  phasicFlow File                                                            | 
+|  copyright: www.cemf.ir                                                     | 
+\* ------------------------------------------------------------------------- */  
+
+objectName     geometryDict;
+objectType     dictionary;
+fileFormat     ASCII;
+
+motionModel rotatingAxis;         // motion model: rotating object around an axis
+
+rotatingAxisInfo                  // information for rotatingAxisMotion motion model 
+{
+    rotAxis 
+    {
+        p1    (0.0 0.0 0.0);      // first point for the axis of rotation    
+        
+        p2    (0.0 0.0 1.0);      // second point for the axis of rotation    
+        
+        omega 1.256;              // rotation speed (rad/s) => 12 rpm
+    }
+}
+
+surfaces
+{
+    cylinder
+    {
+        type        cylinderWall;          // type of the wall
+        
+        p1          (0.0 0.0 0.0);         // begin point of cylinder axis
+           
+        p2          (0.0 0.0 0.8);         // end point of cylinder axis
+        
+        radius1     0.2;                   // radius at p1
+        
+        radius2     0.2;                   // radius at p2
+        
+        resolution  24;                    // number of divisions
+        
+        material    wallMat;               // material name of this wall
+        
+        motion      rotAxis;               // motion component name 
+    }
+
+    /*
+        This is a plane wall at the rear end of cylinder
+    */
+    
+    wall1
+    {
+        type      planeWall;               // type of the wall  
+        
+        p1        (-0.2 -0.2 0.0);         // first point of the wall
+        
+        p2        ( 0.2 -0.2 0.0);         // second point
+        
+        p3        ( 0.2  0.2 0.0);         // third point
+        
+        p4        (-0.2  0.2 0.0);         // fourth point 
+        
+        material  wallMat;                 // material name of the wall  
+        
+        motion    rotAxis;                 // motion component name 
+    }
+
+    /*
+        This is a plane wall at the front end of cylinder
+    */
+    
+    wall2
+    {
+        type      planeWall;               // type of the wall
+        
+        p1        (-0.2 -0.2 0.8);         // first point of the wall
+        
+        p2        ( 0.2 -0.2 0.8);         // second point
+        
+        p3        ( 0.2  0.2 0.8);         // third point
+        
+        p4        (-0.2  0.2 0.8);         // fourth point 
+        
+        material  wallMat;                 // material name of the wall  
+        
+        motion    rotAxis;                 // motion component name 
+    }
+}

benchmarks/rotatingDrum/rotatingDrum_500kParticles/settings/particlesDict

@@ -0,0 +1,47 @@
+/* -------------------------------*- C++ -*--------------------------------- *\ 
+|  phasicFlow File                                                            | 
+|  copyright: www.cemf.ir                                                     | 
+\* ------------------------------------------------------------------------- */  
+
+objectName     particlesDict;
+objectType     dictionary;
+fileFormat     ASCII;
+
+setFields
+{
+    defaultValue 
+    {
+        velocity         realx3     (0 0 0);     // linear velocity (m/s)
+        acceleration     realx3     (0 0 0);     // linear acceleration (m/s2)
+        rotVelocity      realx3     (0 0 0);     // rotational velocity (rad/s)
+        shapeName        word       glassBead;   // name of the particle shape 
+    }
+
+    selectors
+    {}
+}
+
+positionParticles
+{
+    method ordered;
+
+    orderedInfo
+    {
+        diameter   0.005;           // minimum space between centers of particles
+        
+        numPoints  500000;          // number of particles in the simulation 
+        
+        axisOrder  (z x y);         // axis order for filling the space with particles
+    }
+
+    regionType cylinder;            // other options: box and sphere    
+
+    cylinderInfo                    // cylinder for positioning particles 
+    {
+        p1     (0.0  0.0  0.01);    // lower corner point of the box 
+        
+        p2     (0.0  0.0  0.79);    // upper corner point of the box 
+        
+        radius 0.195;               // radius of cylinder 
+    }
+}

benchmarks/rotatingDrum/rotatingDrum_500kParticles/settings/settingsDict

@@ -0,0 +1,34 @@
+/* -------------------------------*- C++ -*--------------------------------- *\ 
+|  phasicFlow File                                                            | 
+|  copyright: www.cemf.ir                                                     | 
+\* ------------------------------------------------------------------------- */  
+objectName    settingsDict;
+objectType    dictionary;
+fileFormat    ASCII;
+/*---------------------------------------------------------------------------*/
+run             rotatingDrum_500KParticles;
+
+dt              0.00001;      // time step for integration (s)
+
+startTime       0;            // start time for simulation 
+
+endTime         4;           // end time for simulation 
+
+saveInterval    0.2;          // time interval for saving the simulation
+
+timePrecision   5;            // maximum number of digits for time folder 
+
+g               (0 -9.8 0);   // gravity vector (m/s2) 
+
+includeObjects (diameter);                     // save necessary (i.e., required) data on disk
+
+// exclude unnecessary data from saving on disk
+excludeObjects (rVelocity.dy1 pStructPosition.dy1 pStructVelocity.dy1); 
+
+integrationMethod         AdamsBashforth2;     // integration method 
+
+writeFormat               binary;              // data writting format (ascii or binary)
+
+timersReport              Yes;      
+
+timersReportInterval      0.01;    

benchmarks/rotatingDrum/rotatingDrum_8mParticles/caseSetup/interaction

@@ -0,0 +1,60 @@
+/* -------------------------------*- C++ -*--------------------------------- *\ 
+|  phasicFlow File                                                            | 
+|  copyright: www.cemf.ir                                                     | 
+\* ------------------------------------------------------------------------- */ 
+
+objectName     interaction;
+objectType     dictionary;
+fileFormat     ASCII;
+
+/*---------------------------------------------------------------------------*/
+
+materials   (glassMat wallMat); // a list of materials names
+
+densities   (2500.0 2500);      // density of materials [kg/m3]
+
+contactListType   sortedContactList;
+
+contactSearch
+{
+    method         NBS;
+
+    updateInterval   20;
+
+    sizeRatio       1.1;
+
+    cellExtent     0.55;
+
+    adjustableBox   Yes;
+}
+
+model
+{
+    contactForceModel    nonLinearLimited;
+
+    rollingFrictionModel normal;
+    
+    /*
+       Property (glassMat-glassMat      glassMat-wallMat
+                                        wallMat-wallMat);
+    */
+
+    Yeff (1.0e6  1.0e6  
+                 1.0e6);                                    // Young modulus [Pa]
+   
+    Geff (0.8e6  0.8e6  
+                 0.8e6);                                    // Shear modulus [Pa]
+
+    nu   (0.25   0.25  
+                 0.25);                                     // Poisson's ratio [-]
+
+    en   (0.97   0.85   
+                 1.00);                                     // coefficient of normal restitution
+
+    mu   (0.65   0.65   
+                 0.65);                                     // dynamic friction 
+
+    mur  (0.1    0.1    
+                 0.1);                                      // rolling friction 
+}
+

benchmarks/rotatingDrum/rotatingDrum_8mParticles/caseSetup/shapes

@@ -0,0 +1,15 @@
+/* -------------------------------*- C++ -*--------------------------------- *\ 
+|  phasicFlow File                                                            | 
+|  copyright: www.cemf.ir                                                     | 
+\* ------------------------------------------------------------------------- */ 
+
+objectName     shapes;
+objectType     dictionary;
+fileFormat     ASCII;
+/*---------------------------------------------------------------------------*/
+
+names          (glassBead);    // names of shapes 
+
+diameters      (0.002);        // diameter of shapes 
+
+materials      (glassMat);     // material names for shapes 

benchmarks/rotatingDrum/rotatingDrum_8mParticles/cleanThisCase

@@ -0,0 +1,7 @@
+#!/bin/sh
+cd ${0%/*} || exit 1    # Run from this directory
+
+ls | grep -P "^(([0-9]+\.?[0-9]*)|(\.[0-9]+))$" | xargs -d"\n" rm -rf
+rm -rf VTK 
+
+#------------------------------------------------------------------------------

benchmarks/rotatingDrum/rotatingDrum_8mParticles/runThisCase

@@ -0,0 +1,21 @@
+#!/bin/sh
+cd ${0%/*} || exit 1    # Run from this directory
+echo "\n<--------------------------------------------------------------------->"
+echo "1) Creating particles"
+echo "<--------------------------------------------------------------------->\n"
+particlesPhasicFlow
+
+echo "\n<--------------------------------------------------------------------->"
+echo "2) Creating geometry"
+echo "<--------------------------------------------------------------------->\n"
+geometryPhasicFlow
+
+echo "\n<--------------------------------------------------------------------->"
+echo "3) Running the case"
+echo "<--------------------------------------------------------------------->\n"
+sphereGranFlow
+
+
+
+
+#------------------------------------------------------------------------------

benchmarks/rotatingDrum/rotatingDrum_8mParticles/settings/domainDict

@@ -0,0 +1,50 @@
+/* -------------------------------*- C++ -*--------------------------------- *\ 
+|  phasicFlow File                                                            | 
+|  copyright: www.cemf.ir                                                     | 
+\* ------------------------------------------------------------------------- */  
+objectName     domainDict;
+objectType     dictionary;
+fileFormat     ASCII;
+/*---------------------------------------------------------------------------*/
+globalBox                                        // Simulation domain: every particles that goes outside this domain will be deleted
+{
+    min (-0.2 -0.2  0.0);
+    max ( 0.2  0.2  1.6);
+}
+
+boundaries
+{
+    neighborListUpdateInterval 200;
+    
+    updateInterval              20;
+    
+    left
+    {
+        type    exit;                            // other options: periodic, reflective 
+    }
+
+    right
+    {
+        type     exit;                           // other options: periodic, reflective 
+    }
+
+    bottom
+    {
+        type     exit;                           // other options: periodic, reflective 
+    }
+
+    top
+    {
+        type     exit;                           // other options: periodic, reflective 
+    }
+
+    rear
+    {
+        type     exit;                           // other options: periodic, reflective 
+    }
+
+    front
+    {
+        type     exit;                           // other options: periodic, reflective 
+    }
+}

benchmarks/rotatingDrum/rotatingDrum_8mParticles/settings/geometryDict

@@ -0,0 +1,86 @@
+/* -------------------------------*- C++ -*--------------------------------- *\ 
+|  phasicFlow File                                                            | 
+|  copyright: www.cemf.ir                                                     | 
+\* ------------------------------------------------------------------------- */  
+
+objectName     geometryDict;
+objectType     dictionary;
+fileFormat     ASCII;
+
+motionModel rotatingAxis;         // motion model: rotating object around an axis
+
+rotatingAxisInfo                  // information for rotatingAxisMotion motion model 
+{
+    rotAxis 
+    {
+        p1    (0.0 0.0 0.0);      // first point for the axis of rotation    
+        
+        p2    (0.0 0.0 1.0);      // second point for the axis of rotation    
+        
+        omega 1.256;              // rotation speed (rad/s) => 12 rpm
+    }
+}
+
+surfaces
+{
+    cylinder
+    {
+        type        cylinderWall;          // type of the wall
+        
+        p1          (0.0 0.0 0.0);         // begin point of cylinder axis
+           
+        p2          (0.0 0.0 1.6);         // end point of cylinder axis
+        
+        radius1     0.2;                   // radius at p1
+        
+        radius2     0.2;                   // radius at p2
+        
+        resolution  24;                    // number of divisions
+        
+        material    wallMat;               // material name of this wall
+        
+        motion      rotAxis;               // motion component name 
+    }
+
+    /*
+        This is a plane wall at the rear end of cylinder
+    */
+    
+    wall1
+    {
+        type      planeWall;               // type of the wall  
+        
+        p1        (-0.2 -0.2 0.0);         // first point of the wall
+        
+        p2        ( 0.2 -0.2 0.0);         // second point
+        
+        p3        ( 0.2  0.2 0.0);         // third point
+        
+        p4        (-0.2  0.2 0.0);         // fourth point 
+        
+        material  wallMat;                 // material name of the wall  
+        
+        motion    rotAxis;                 // motion component name 
+    }
+
+    /*
+        This is a plane wall at the front end of cylinder
+    */
+    
+    wall2
+    {
+        type      planeWall;               // type of the wall
+        
+        p1        (-0.2 -0.2 1.6);         // first point of the wall
+        
+        p2        ( 0.2 -0.2 1.6);         // second point
+        
+        p3        ( 0.2  0.2 1.6);         // third point
+        
+        p4        (-0.2  0.2 1.6);         // fourth point 
+        
+        material  wallMat;                 // material name of the wall  
+        
+        motion    rotAxis;                 // motion component name 
+    }
+}

benchmarks/rotatingDrum/rotatingDrum_8mParticles/settings/particlesDict

@@ -0,0 +1,47 @@
+/* -------------------------------*- C++ -*--------------------------------- *\ 
+|  phasicFlow File                                                            | 
+|  copyright: www.cemf.ir                                                     | 
+\* ------------------------------------------------------------------------- */  
+
+objectName     particlesDict;
+objectType     dictionary;
+fileFormat     ASCII;
+
+setFields
+{
+    defaultValue 
+    {
+        velocity         realx3     (0 0 0);     // linear velocity (m/s)
+        acceleration     realx3     (0 0 0);     // linear acceleration (m/s2)
+        rotVelocity      realx3     (0 0 0);     // rotational velocity (rad/s)
+        shapeName        word       glassBead;   // name of the particle shape 
+    }
+
+    selectors
+    {}
+}
+
+positionParticles
+{
+    method ordered;
+
+    orderedInfo
+    {
+        diameter   0.003;           // minimum space between centers of particles
+        
+        numPoints  6000000;          // number of particles in the simulation 
+        
+        axisOrder  (z x y);         // axis order for filling the space with particles
+    }
+
+    regionType cylinder;            // other options: box and sphere    
+
+    cylinderInfo                    // cylinder for positioning particles 
+    {
+        p1     (0.0  0.0  0.01);    // lower corner point of the box 
+        
+        p2     (0.0  0.0  1.59);    // upper corner point of the box 
+        
+        radius 0.195;               // radius of cylinder 
+    }
+}

benchmarks/rotatingDrum/rotatingDrum_8mParticles/settings/settingsDict

@@ -0,0 +1,34 @@
+/* -------------------------------*- C++ -*--------------------------------- *\ 
+|  phasicFlow File                                                            | 
+|  copyright: www.cemf.ir                                                     | 
+\* ------------------------------------------------------------------------- */  
+objectName    settingsDict;
+objectType    dictionary;
+fileFormat    ASCII;
+/*---------------------------------------------------------------------------*/
+run             rotatingDrum_4mParticles;
+
+dt              0.00001;      // time step for integration (s)
+
+startTime       0;            // start time for simulation 
+
+endTime         4;           // end time for simulation 
+
+saveInterval    0.2;          // time interval for saving the simulation
+
+timePrecision   5;            // maximum number of digits for time folder 
+
+g               (0 -9.8 0);   // gravity vector (m/s2) 
+
+includeObjects (diameter);                     // save necessary (i.e., required) data on disk
+
+// exclude unnecessary data from saving on disk
+excludeObjects (rVelocity.dy1 pStructPosition.dy1 pStructVelocity.dy1); 
+
+integrationMethod         AdamsBashforth2;     // integration method 
+
+writeFormat               binary;              // data writting format (ascii or binary)
+
+timersReport              Yes;      
+
+timersReportInterval      0.01;    

benchmarks/rotatingDrum_4MParticles/caseSetup/interaction

@@ -1,59 +0,0 @@
-/* -------------------------------*- C++ -*--------------------------------- *\ 
-|  phasicFlow File                                                            | 
-|  copyright: www.cemf.ir                                                     | 
-\* ------------------------------------------------------------------------- */ 
-
-objectName interaction;
-objectType dicrionary;
-
-materials   (glassMat wallMat); // a list of materials names
-densities   (2500.0 2500);      // density of materials [kg/m3]
-
-contactListType   sortedContactList;
-
-model
-{
-   contactForceModel nonLinearLimited;
-   rollingFrictionModel normal;
-
-   Yeff (1.0e6 1.0e6  // Young modulus [Pa]
-               1.0e6);
-   
-   Geff (0.8e6 0.8e6  // Shear modulus [Pa]
-               0.8e6);
-
-   nu    (0.25 0.25  // Poisson's ratio [-]
-               0.25);
-
-   en (0.97  0.85   // coefficient of normal restitution
-             1.00);
-
-   et (1.0   1.0    // coefficient of tangential restitution 
-             1.0);
-
-   mu (0.65  0.65   // dynamic friction 
-             0.65);
-
-   mur (0.1  0.1    // rolling friction 
-             0.1);
-        
-}
-
-contactSearch
-{
-   method         NBS;           
-   wallMapping    cellMapping;   
-
-   NBSInfo
-   {
-      updateFrequency 10;        // each 20 timesteps, update neighbor list 
-      sizeRatio      1.05;       // bounding box size to particle diameter (max)
-   }
-
-   cellMappingInfo
-   {
-      updateFrequency 10;        // each 20 timesteps, update neighbor list  
-      cellExtent     0.6;        // bounding box for particle-wall search (> 0.5)
-   }
-
-}

benchmarks/rotatingDrum_4MParticles/caseSetup/particleInsertion

@@ -1,14 +0,0 @@
-/* -------------------------------*- C++ -*--------------------------------- *\ 
-|  phasicFlow File                                                            | 
-|  copyright: www.cemf.ir                                                     | 
-\* ------------------------------------------------------------------------- */ 
-
-objectName particleInsertion;
-objectType dicrionary;
-
-
-active  no;          // is insertion active?
-
-collisionCheck No;   // not implemented for yes
-
-

benchmarks/rotatingDrum_4MParticles/caseSetup/sphereShape

@@ -1,11 +0,0 @@
-/* -------------------------------*- C++ -*--------------------------------- *\ 
-|  phasicFlow File                                                            | 
-|  copyright: www.cemf.ir                                                     | 
-\* ------------------------------------------------------------------------- */ 
-
-objectName 	sphereDict;
-objectType 	sphereShape;
-
-names 		(glassBead); 	// names of shapes 
-diameters 	(0.003);		// diameter of shapes 
-materials	(glassMat);		// material names for shapes 

benchmarks/rotatingDrum_4MParticles/settings/geometryDict

@@ -1,63 +0,0 @@
-/* -------------------------------*- C++ -*--------------------------------- *\ 
-|  phasicFlow File                                                            | 
-|  copyright: www.cemf.ir                                                     | 
-\* ------------------------------------------------------------------------- */  
-
-objectName 	geometryDict;
-objectType 	dictionary;
-
-motionModel rotatingAxisMotion; 
-
-surfaces
-{
-	
-	cylinder
-	{
-		type cylinderWall;  
-		p1 (0.0 0.0 0.0);   
-		p2 (0.0 0.0 1.6);  
-		radius1 0.2;		
-		radius2 0.2;		
-		resolution 24;      
-		material wallMat;   
-		motion rotAxis;
-	}
-
-	
-	wall1
-	{
-		type planeWall;			
-		p1 (-0.2 -0.2 0.0);	
-		p2 ( 0.2 -0.2 0.0);   
-		p3 ( 0.2  0.2 0.0);   
-		p4 (-0.2  0.2 0.0);   
-		material wallMat;     
-		motion rotAxis;		
-	}
-
-	/*
-	This is a plane wall at the front end of cylinder
-	*/
-	wall2
-	{
-		type planeWall;
-		p1 (-0.2 -0.2 1.6);
-		p2 ( 0.2 -0.2 1.6);
-		p3 ( 0.2  0.2 1.6);
-		p4 (-0.2  0.2 1.6);
-		material wallMat;
-		motion rotAxis;
-	}
-		
-}
-
-// information for rotatingAxisMotion motion model 
-rotatingAxisMotionInfo
-{
-	rotAxis 
-	{
-		p1 (0.0 0.0 0.0);	
-		p2 (0.0 0.0 1.0);	
-		omega 1.256; 		// rotation speed (rad/s) => 12 rpm
-	}
-}

benchmarks/rotatingDrum_4MParticles/settings/particlesDict

@@ -1,44 +0,0 @@
-/* -------------------------------*- C++ -*--------------------------------- *\ 
-|  phasicFlow File                                                            | 
-|  copyright: www.cemf.ir                                                     | 
-\* ------------------------------------------------------------------------- */  
-
-objectName 	particlesDict;
-objectType 	dictionary;
-
-setFields
-{
-	
-	defaultValue 
-	{
-		velocity 		realx3 	(0 0 0); // linear velocity (m/s)
-		acceleration 	realx3 	(0 0 0); // linear acceleration (m/s2)
-		rotVelocity 	realx3 	(0 0 0); // rotational velocity (rad/s)
-		shapeName 		word	glassBead; // name of the particle shape 
-	}
-
-	selectors
-	{}
-}
-
-positionParticles
-{
-	method positionOrdered;     		
-
-	maxNumberOfParticles 	4000001; 
-	mortonSorting 			Yes;    
-
-	cylinder  // box for positioning particles 
-	{
-		p1 ( 0.0  0.0 0.01);    // lower corner point of the box 
-		p2 ( 0.0  0.0 1.59);	// upper corner point of the box 
-		radius 0.195;
-	}
-	
-	positionOrderedInfo
-	{
-		diameter 0.003; 	// minimum space between centers of particles
-		numPoints 4000000; 	// number of particles in the simulation 
-		axisOrder (z x y);  // axis order for filling the space with particles
-	}
-}

benchmarks/rotatingDrum_4MParticles/settings/settingsDict

@@ -1,32 +0,0 @@
-/* -------------------------------*- C++ -*--------------------------------- *\ 
-|  phasicFlow File                                                            | 
-|  copyright: www.cemf.ir                                                     | 
-\* ------------------------------------------------------------------------- */  
-objectName 	settingsDict;
-objectType 	dictionary;;
-
-run 			rotatingDrum_1;
-
-dt 				0.00001; 	// time step for integration (s)
-
-startTime 		0; 			// start time for simulation 
-
-endTime 		10; 		// end time for simulation 
-
-saveInterval 	0.2; 		// time interval for saving the simulation
-
-timePrecision   5;			// maximum number of digits for time folder 
-
-g 				(0 -9.8 0); // gravity vector (m/s2) 
-
-domain 
-{
-	min (-0.2 -0.2 -0.0);
-	max ( 0.2  0.2  1.6);
-}
-
-integrationMethod 		AdamsBashforth3; 	// integration method 
-
-timersReport 			Yes;  	
-
-timersReportInterval   	0.01;	

cmake/autoComplete

@@ -1,17 +1,96 @@
 PF_cFlags="--description --help --version"
+AllTimeFolders=
+__getAllTime(){
+    # Initialize empty array for time folders
+    local time_folders=()
+    
+    # Loop through all directories in current folder
+    for dir in */; do
+        # Remove trailing slash
+        dir=${dir%/}
+        
+        # Check if directory name is a valid floating point number
+        # This pattern matches integers and floating point numbers
+        if [[ $dir =~ ^[0-9]+(\.[0-9]+)?$ ]]; then
+            time_folders+=("$dir")
+        fi
+    done
+    
+    # Set completion reply to the time folders
+    COMPREPLY=("${time_folders[@]}")
+    AllTimeFolders="${time_folders[@]}"
+}
+
+__getFields(){
+    __getAllTime
+    local -A unique_files=()
+    # Files to exclude from suggestions
+    local exclude_files=("shapeHash" "pStructure" "particleInsertion" "p" "alpha" "U" "Sp" "Su" "phi")
+    declare -A exclude_dict
+    
+    # Build exclude dictionary for faster lookups
+    for file in "${exclude_files[@]}"; do
+        exclude_dict["$file"]=1
+    done
+
+    for dir in $AllTimeFolders; do
+        # Skip if not a directory
+        [ ! -d "$dir" ] && continue
+        
+        # Find all files in this time directory
+        while IFS= read -r filename; do
+            # Skip empty lines and excluded files
+            [ -z "$filename" ] || [ "${exclude_dict[$filename]+exists}" ] && continue
+            
+            # Add to unique files
+            unique_files["$filename"]=1
+        done < <(find "$dir" -maxdepth 1 -type f -printf '%f\n')
+    done
+    
+    # Set completion reply to the unique filenames
+    COMPREPLY=(${!unique_files[@]})
+    
+    # Clear global variable
+    AllTimeFolders=
+}
 
 _pFlowToVTK(){
-    if [ "$3" == "--time" ]; then
+    local cur="${COMP_WORDS[COMP_CWORD]}"
-      COMPREPLY=( $(ls) )
+    local prev="${COMP_WORDS[COMP_CWORD-1]}"
+    
+    # Check if we're completing a field
+    local is_field=0
+    for ((i=1; i<COMP_CWORD; i++)); do
+        if [[ "${COMP_WORDS[i]}" == "--fields" ]]; then
+            is_field=1
+            break
+        fi
+    done
+    
+    if [ "$prev" == "--time" ]; then
+        __getAllTime
+    elif [ "$prev" == "--fields" ] || [ $is_field -eq 1 ]; then 
+        # We're completing field names
+        __getFields
+        # Filter the results based on the current word prefix
+        if [ -n "$cur" ]; then
+            local filtered=()
+            for item in "${COMPREPLY[@]}"; do
+                if [[ "$item" == "$cur"* ]]; then
+                    filtered+=("$item")
+                fi
+            done
+            COMPREPLY=("${filtered[@]}")
+        fi
     else
-      COMPREPLY=( $(compgen -W "$PF_cFlags --binary --no-geometry --no-particles --out-folder --time --separate-surfaces --fields" -- "$2") )
+        COMPREPLY=( $(compgen -W "$PF_cFlags --binary --no-geometry --no-particles --out-folder --time --separate-surfaces --fields" -- "$cur") )
     fi
 }
 complete -F _pFlowToVTK pFlowToVTK
 
 _postprocessPhasicFlow(){
     if [ "$3" == "--time" ]; then
-      COMPREPLY=( $(ls) )
+      __getAllTime
     else
       COMPREPLY=( $(compgen -W "$PF_cFlags --out-folder --time --zeroFolder" -- "$2") )
     fi

cmake/preReq.cmake

@@ -0,0 +1,54 @@
+if(pFlow_STD_Parallel_Alg)
+    # Check if libtbb-dev is installed
+    execute_process(
+        COMMAND dpkg -s libtbb-dev
+        RESULT_VARIABLE TBB_IS_INSTALLED
+        OUTPUT_QUIET
+        ERROR_QUIET)
+
+    if(NOT TBB_IS_INSTALLED EQUAL 0)
+        message(STATUS "libtbb-dev not found. Installing libtbb-dev...")
+        execute_process(
+            COMMAND sudo apt-get update
+            COMMAND sudo apt-get install -y libtbb-dev
+            RESULT_VARIABLE TBB_INSTALL_RESULT)
+            
+        if(NOT TBB_INSTALL_RESULT EQUAL 0)
+            message(FATAL_ERROR "Failed to install libtbb-dev")
+        endif()
+    else()
+        message(STATUS "libtbb-dev is already installed.")
+    endif()
+endif()
+
+
+# Kokkos folder creation 
+set(Kokkos_Source_DIR $ENV{HOME}/Kokkos/kokkos)
+
+if(NOT EXISTS "${Kokkos_Source_DIR}/CMakeLists.txt")
+
+    # Check CMake version and set policy CMP0169 if CMake version is 3.30 or higher
+    if(${CMAKE_VERSION} VERSION_GREATER_EQUAL "3.30")
+        cmake_policy(SET CMP0169 OLD)
+    endif()
+
+    include(FetchContent)
+    FetchContent_Declare(
+        kokkos
+        GIT_REPOSITORY https://github.com/kokkos/kokkos.git
+        GIT_TAG 4.4.01
+    )
+    
+    FetchContent_GetProperties(kokkos)
+    if(NOT kokkos_POPULATED)
+        message(STATUS "Kokkos source directory not found. Downloading Kokkos version 4.4.1 ...")
+        FetchContent_Populate(kokkos)
+        set(Kokkos_Source_DIR ${kokkos_SOURCE_DIR})
+    endif()
+
+endif()
+
+message(STATUS "Kokkos source directory is ${Kokkos_Source_DIR}")
+add_subdirectory(${Kokkos_Source_DIR} ./kokkos)
+#Kokkos_cmake_settings()
+

doc/mdDocs/codingStyle.md

@@ -0,0 +1,231 @@
+# PhasicFlow Coding Style Guidelines
+This document outlines the coding style guidelines for the PhasicFlow codebase. 
+Adhering to these guidelines ensures consistency, readability, and maintainability of the project.
+## 1. FormattingIndentation: 
+* Use 4 spaces for every logical level and block.
+* Line Spacing: Leave two empty lines between sections (e.g., between functions in a .cpp file, between class members).
+
+## 2. Naming ConventionsGeneral Naming: 
+* All names should start with lowercase letters, except for special names (e.g., Ergun, Hertz).
+* Macro names start with Upper case or all the letters are in UPPER case.
+* Compound Names: For compound names, the first word starts with a lowercase letter, and subsequent words start with an uppercase letter (e.g., boundaryBase, motionModel).
+
+## 3. File Structure 
+* Header Files: Use the .hpp extension for header files.
+* Source Files: Use the .cpp extension for source files.
+* Header and Source File Headers: All header and source files must include a standardized header that describes the project's intention and licensing information.
+* File Naming: Header and source file names should correspond to the class they contain. Aim for one class per file.
+* Inline Functions: Place inline functions in a separate classNameI.hpp file to avoid cluttering the main header file.
+
+## 4. Class DesignClass Member Order:
+* Private members and methods
+* Private static members and methods
+* Public methods
+* Public static methods
+* Enumerations and Nested Classes: Declare enumerations and nested classes before all class members and methods.
+* Special Functions: Each class must explicitly define all special functions:Constructor, Copy constructor and assignment operator, Move constructor and assignment operator
+* Destructor: Each class must have an explicit destructor declaration:`~className() = default`; or `~className();`
+* Interface classes or classes with virtual methods must have a virtual destructor.
+* Virtual Method Overrides: When implementing a `virtual` method from a base class in a derived class, use the `override` keyword. The same applies to derived class destructors.
+## 5. NamespacesOfficial Namespace: 
+The official namespace for the codebase is pFlow. All entities should be defined within this namespace.
+
+### Example File Structure
+```
+src/
+├── componentName1/
+│   ├── componentName1.hpp
+│   ├── componentName1.cpp
+│   ├── componentName1I.hpp
+│   └── ...
+└── componentName2/
+    ├── componentName2.hpp
+    ├── componentName2.cpp
+    └── ...
+```
+### Example Class Structure
+```C++
+namespace pFlow 
+{
+
+class MyClass 
+{
+public:
+    enum class MyEnum 
+    {
+        Value1,
+        Value2
+    };
+
+    class NestedClass 
+    {
+        // ...
+    };
+
+private:
+    int privateMember_;
+    
+    void privateMethod();
+
+    
+    static int privateStaticMember_;
+    
+    static void privateStaticMethod();
+
+public:
+    
+    MyClass();
+
+    MyClass(const MyClass& other);
+    
+    MyClass(MyClass&& other);
+    
+    MyClass& operator=(const MyClass& other);
+    
+    MyClass& operator=(MyClass&& other);
+    
+    ~MyClass();
+
+    
+    void publicMethod();
+    
+    static void publicStaticMethod();
+
+};
+
+// assuming base class has virtual methods 
+class DerivedClass 
+: 
+    public BaseClass 
+{
+public:
+     
+     ...
+
+    ~DerivedClass() override;
+    
+     void virtualMethod() override;
+};
+
+} // namespace pFlow
+``` 
+
+## 6. Doxygen Documentation
+
+### 6.1. Ruls
+provide the documentations in the header files only. In rare cases where you are generating documentations for executables, the doxygen documentation can be supplied in the .cpp file.  
+
+* **General Doxygen Style:**
+
+    * Use `///` for short documentations for methods and members.
+
+    * Use `/** */` for classes and main methods which play a significant role in the class or code.
+
+    * Place Doxygen comments *before* the declaration of the entity being documented (e.g., class, function, variable).
+
+    * Use `@param` to document function parameters, `@return` for return values, `@brief` for a short description, and `@details` for a more in-depth explanation.
+
+    * Use Markdown syntax within Doxygen comments for formatting.
+
+* **File Headers:** Each file should contain a Doxygen comment at the top, including:
+
+    * `@file` : The name of the file.
+
+    * `@brief`: A brief description of the file's purpose.
+
+    * `@author`: The author(s) of the file.
+
+    * `@date` : The date of creation or last modification.
+
+* **Class Documentation:**
+
+    * Use `/** */` for class documentation.
+
+    * Provide a `@brief` description of the class.
+
+    * Use `@tparam` to document template parameters.
+
+    * Document the purpose of the class, its invariants, and how it should be used.
+
+* **Function/Method Documentation:**
+
+    * Use `///` for short documentations.
+
+    * Use `/** */` for main methods which play a significant role.
+
+    * Provide a `@brief` description of the function.
+
+    * Use `@param` to describe each parameter, including its purpose and whether it is an input, output, or input/output parameter.
+
+    * Use `@return` to describe the return value, including its meaning and possible values.
+
+    * Use `@pre` to document any preconditions that must be met before calling the function.
+
+    * Use `@post` to document any postconditions that will be true after the function returns.
+
+    * Use `@throws` to document any exceptions that the function may throw.
+
+    * Use `@details` for a more detailed explanation of the function's behavior, algorithms, or any other relevant information.
+
+* **Variable Documentation:**
+
+    * Use `///<` for single-line documentation of variables.
+
+### 6.2. Doxygen Documentation Examples
+
+* **Class example**
+
+```cpp
+/**
+ * @brief Represents a particle in the simulation.
+ * @details This class stores the position, velocity, and other physical
+ * properties of a particle.
+ */
+class Particle
+{
+private:
+
+    Point     position_;   ///< The current position of the particle.
+
+    Vector    velocity_;   ///< The current velocity of the particle.
+
+    double     mass_;      ///< The mass of the particle.
+
+public:
+
+    /// Constructs a particle with default values.
+    Particle();
+
+    /**
+     * @brief Updates the position of the particle based on its velocity
+     * and the given time step.
+     * @param deltaTime The time elapsed since the last update, in seconds.
+     */
+    void updatePosition(const timeInfo& ti );
+
+    /// Gets the current position of the particle.
+    Point getPosition() const;
+};
+```
+
+* **Function Example**
+
+```cpp
+/**
+ * @brief Calculates the distance between two points.
+ * @param p1 The first point.
+ * @param p2 The second point.
+ * @return The distance between the two points.
+ */
+double calculateDistance(const Point& p1, const Point& p2)
+{
+    // Implementation
+    return 0.0;
+}
+
+/// Returns the velocity of the particle.
+Vector getVelocity() const
+{
+    return velocity_;
+}
+``` 

solvers/grainGranFlow/createDEMComponents.hpp

@@ -17,10 +17,17 @@ Licence:
   implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
 
 -----------------------------------------------------------------------------*/
+REPORT(0)<<"Reading shape dictionary ..."<<END_REPORT;
+pFlow::grainShape grains
+(
+	pFlow::shapeFile__,
+	&Control.caseSetup(),
+	proprties
+);
 
 //
 REPORT(0)<<"\nReading sphere particles . . ."<<END_REPORT;
-pFlow::grainParticles grnParticles(Control, proprties);
+pFlow::grainParticles grnParticles(Control, grains);
 
 //
 REPORT(0)<<"\nCreating particle insertion object . . ."<<END_REPORT;

solvers/iterateSphereParticles/createDEMComponents.hpp

@@ -18,8 +18,16 @@ Licence:
 
 -----------------------------------------------------------------------------*/
 
+REPORT(0)<<"Reading shape dictionary ..."<<END_REPORT;
+pFlow::sphereShape spheres
+(
+	pFlow::shapeFile__,
+	&Control.caseSetup(),
+	proprties
+);
+
 //
 REPORT(0)<<"\nReading sphere particles . . ."<<END_REPORT;
-pFlow::sphereParticles sphParticles(Control, proprties);
+pFlow::sphereParticles sphParticles(Control, spheres);
 
 WARNING<<"Particle insertion has not been set yet!"<<END_WARNING;

solvers/sphereGranFlow/createDEMComponents.hpp

@@ -17,10 +17,18 @@ Licence:
   implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
 
 -----------------------------------------------------------------------------*/
-
+REPORT(0)<<"Reading shapes dictionary..."<<END_REPORT;
+pFlow::sphereShape spheres
+(
+	pFlow::shapeFile__,
+	&Control.caseSetup(),
+	proprties
+);
 //
 REPORT(0)<<"\nReading sphere particles . . ."<<END_REPORT;
-pFlow::sphereParticles sphParticles(Control, proprties);
+
+
+pFlow::sphereParticles sphParticles(Control, spheres);
 
 //
 REPORT(0)<<"\nCreating particle insertion object . . ."<<END_REPORT;

solvers/sphereGranFlow/sphereGranFlow.cpp

@@ -74,7 +74,7 @@ pFlow::initialize_pFlowProcessors();
 
 	do 
 	{
-		
+		//Ping;
 		if(! sphInsertion.insertParticles( 
 			Control.time().currentIter(),
 			Control.time().currentTime(),
@@ -90,21 +90,25 @@ pFlow::initialize_pFlowProcessors();
 
 		// set force to zero, predict, particle deletion and etc. 
 		sphParticles.beforeIteration();
-
+		//Ping;
 		sphInteraction.beforeIteration();
 		
 		sphInteraction.iterate();
 			
 		surfGeometry.iterate();
 
+		//Ping;
 		sphParticles.iterate();
 		
+		//Ping;
 		sphInteraction.afterIteration();
 
+		//Ping;
 		surfGeometry.afterIteration();
 
+		//Ping;
 		sphParticles.afterIteration();
-				
+		//Ping;	
 
 	}while(Control++);
 

src/CMakeLists.txt

@@ -13,3 +13,5 @@ add_subdirectory(Interaction)
 
 add_subdirectory(MotionModel)
 
+add_subdirectory(PostprocessData)
+

src/Geometry/geometryMotion/geometryMotion.cpp

@@ -21,7 +21,6 @@ Licence:
 template<typename MotionModel>
 bool pFlow::geometryMotion<MotionModel>::findMotionIndex()
 {
-	
 	if(motionComponentName().size() != numSurfaces() )
 	{
 		fatalErrorInFunction<<

src/Geometry/geometryMotion/geometryMotions.cpp

@@ -28,4 +28,4 @@ template class pFlow::geometryMotion<pFlow::stationaryWall>;
 
 template class pFlow::geometryMotion<pFlow::conveyorBeltMotion>;
 
-//template class pFlow::geometryMotion<pFlow::multiRotatingAxisMotion>;
+template class pFlow::geometryMotion<pFlow::multiRotatingAxisMotion>;

src/Geometry/geometryMotion/geometryMotions.hpp

@@ -25,7 +25,7 @@ Licence:
 #include "stationaryWall.hpp"
 #include "rotatingAxisMotion.hpp"
 #include "conveyorBeltMotion.hpp"
-//#include "multiRotatingAxisMotion.hpp"
+#include "multiRotatingAxisMotion.hpp"
 #include "vibratingMotion.hpp"
 
 
@@ -40,10 +40,7 @@ using stationaryGeometry = geometryMotion<stationaryWall>;
 
 using conveyorBeltMotionGeometry = geometryMotion<conveyorBeltMotion>;
 
-//typedef geometryMotion<multiRotatingAxisMotion> multiRotationAxisMotionGeometry;
+using multiRotationAxisMotionGeometry = geometryMotion<multiRotatingAxisMotion>;
-
-
-
 
 }
 

src/Integration/AdamsBashforth2/AdamsBashforth2.cpp

@@ -37,7 +37,8 @@ bool intAllActive(
 	real dt, 
 	realx3Field_D& y, 
 	realx3PointField_D& dy,
-	realx3PointField_D& dy1)
+	realx3PointField_D& dy1,
+	real damping = 1.0)
 {
 
 	auto d_dy = dy.deviceView();
@@ -49,7 +50,7 @@ bool intAllActive(
 		"AdamsBashforth2::correct",
 		rpIntegration (activeRng.start(), activeRng.end()),
 		LAMBDA_HD(uint32 i){
-			d_y[i] +=  dt*(static_cast<real>(1.5) * d_dy[i] - static_cast<real>(0.5) * d_dy1[i]);
+			d_y[i] +=  damping * dt*(static_cast<real>(1.5) * d_dy[i] - static_cast<real>(0.5) * d_dy1[i]);
 			d_dy1[i] = d_dy[i];
 		});
 	Kokkos::fence();
@@ -62,7 +63,8 @@ bool intScattered
 	real dt, 
 	realx3Field_D& y,
 	realx3PointField_D& dy,
-	realx3PointField_D& dy1
+	realx3PointField_D& dy1,
+	real damping = 1.0
 )
 {
 
@@ -78,7 +80,7 @@ bool intScattered
 		LAMBDA_HD(uint32 i){
 			if( activeP(i))
 			{
-				d_y[i] +=  dt*(static_cast<real>(1.5) * d_dy[i] - static_cast<real>(0.5) * d_dy1[i]);
+				d_y[i] +=  damping * dt*(static_cast<real>(1.5) * d_dy[i] - static_cast<real>(0.5) * d_dy1[i]);
 				d_dy1[i] = d_dy[i];
 			}
 		});
@@ -95,10 +97,11 @@ pFlow::AdamsBashforth2::AdamsBashforth2
 	const word& baseName,
 	pointStructure& pStruct,
 	const word& method,
-	const realx3Field_D& initialValField
+	const realx3Field_D& initialValField,
+	bool  keepHistory
 )
 :
-	integration(baseName, pStruct, method, initialValField),
+	integration(baseName, pStruct, method, initialValField, keepHistory),
 	realx3PointField_D
 	(
 		objectFile
@@ -106,14 +109,17 @@ pFlow::AdamsBashforth2::AdamsBashforth2
 			groupNames(baseName,"dy1"),
 			pStruct.time().integrationFolder(),
 			objectFile::READ_IF_PRESENT,
-			objectFile::WRITE_ALWAYS
+			keepHistory?objectFile::WRITE_ALWAYS:objectFile::WRITE_NEVER
 		),
 		pStruct,
 		zero3,
 		zero3
 	),
+	initialValField_(initialValField),
 	boundaryList_(pStruct, method, *this)
-{}
+{
+	realx3PointField_D::addEvent(message::ITEMS_INSERT);
+}
 
 void pFlow::AdamsBashforth2::updateBoundariesSlaveToMasterIfRequested()
 {
@@ -142,18 +148,19 @@ bool pFlow::AdamsBashforth2::correct
 (
 	real dt,
 	realx3PointField_D& y,
-	realx3PointField_D& dy
+	realx3PointField_D& dy,
+	real damping
 )
 {
 	auto& dy1l = dy1();
 	bool success = false;
 	if(dy1l.isAllActive())
 	{
-		success = intAllActive(dt, y.field(), dy, dy1l);
+		success = intAllActive(dt, y.field(), dy, dy1(), damping);
 	}
 	else
 	{
-		success = intScattered(dt, y.field(), dy, dy1l);
+		success = intScattered(dt, y.field(), dy, dy1(), damping);
 	}
 
 	success = success && boundaryList_.correct(dt, y, dy);
@@ -171,11 +178,11 @@ bool pFlow::AdamsBashforth2::correctPStruct(
 	bool success = false;
 	if(dy1l.isAllActive())
 	{
-		success = intAllActive(dt, pStruct.pointPosition(), vel, dy1l);
+		success = intAllActive(dt, pStruct.pointPosition(), vel,  dy1());
 	}
 	else
 	{
-		success = intScattered(dt, pStruct.pointPosition(), vel, dy1l);
+		success = intScattered(dt, pStruct.pointPosition(), vel,  dy1());
 	}
 
 	success = success && boundaryList_.correctPStruct(dt, pStruct, vel);
@@ -183,11 +190,21 @@ bool pFlow::AdamsBashforth2::correctPStruct(
 	return success;
 }
 
-
+/*bool pFlow::AdamsBashforth2::hearChanges
-bool pFlow::AdamsBashforth2::setInitialVals(
+(
-	const int32IndexContainer& newIndices,
+	const timeInfo &ti, 
-	const realx3Vector& y)
+	const message &msg, 
+	const anyList &varList
+)
 {
-	return true;
+	if(msg.equivalentTo(message::ITEMS_INSERT))
-}
+	{
 		
+		return insertValues(varList, initialValField_.deviceViewAll(), dy1());
+	}
+	else
+	{
+		return realx3PointField_D::hearChanges(ti, msg, varList);
+	}
+    
+}*/

src/Integration/AdamsBashforth2/AdamsBashforth2.hpp

@@ -41,10 +41,14 @@ class AdamsBashforth2
 {
 private:
 
+	const realx3Field_D& initialValField_;
+
 	boundaryIntegrationList boundaryList_;
 
 	friend class processorAB2BoundaryIntegration;
 
+protected:
+	
 	const auto& dy1()const
 	{
 		return static_cast<const realx3PointField_D&>(*this);
@@ -55,6 +59,16 @@ private:
 		return static_cast<realx3PointField_D&>(*this);
 	}
 
+	auto& initialValField()
+	{
+		return initialValField_;
+	}
+
+	boundaryIntegrationList& boundaryList()
+	{
+		return boundaryList_;
+	}
+	
 public:
 
 	/// Class info
@@ -67,10 +81,11 @@ public:
 			const word& baseName,
 			pointStructure& pStruct,
 			const word& method,
-			const realx3Field_D& initialValField);
+			const realx3Field_D& initialValField,
+			bool  keepHistory);
 
 		/// Destructor 
-		~AdamsBashforth2()final = default;
+		~AdamsBashforth2()override = default;
 
 		/// Add this to the virtual constructor table 
 		add_vCtor(
@@ -102,32 +117,21 @@ public:
 		bool correct(
 			real dt, 
 			realx3PointField_D& y, 
-			realx3PointField_D& dy) final;
+			realx3PointField_D& dy,
+			real damping = 1.0) override;
 
 		bool correctPStruct(
 			real dt, 
 			pointStructure& pStruct, 
-			realx3PointField_D& vel) final;
+			realx3PointField_D& vel) override;
-			
 		
 		/*bool hearChanges
 		(
-			real t,
+			const timeInfo& ti,
-			real dt,
-			uint32 iter,
 			const message& msg, 
 			const anyList& varList
 		) override;*/
 
-		bool setInitialVals(
-			const int32IndexContainer& newIndices,
-			const realx3Vector& y) final;
-
-		bool needSetInitialVals()const final
-		{
-			return false;
-		}
-
 };
 
 

src/Integration/AdamsBashforth3/AdamsBashforth3.cpp

@@ -20,90 +20,185 @@ Licence:
 
 #include "AdamsBashforth3.hpp"
 
-//const real AB3_coef[] = { 23.0 / 12.0, 16.0 / 12.0, 5.0 / 12.0 };
+namespace pFlow
-
-pFlow::AdamsBashforth3::AdamsBashforth3
-(
-	const word& baseName,
-	repository& owner,
-	const pointStructure& pStruct,
-	const word& method
-)
-:
-	integration(baseName, owner, pStruct, method),
-	history_(
-		owner.emplaceObject<pointField<VectorSingle,AB3History>>(
-			objectFile(
-				groupNames(baseName,"AB3History"),
-				"",
-				objectFile::READ_IF_PRESENT,
-				objectFile::WRITE_ALWAYS),
-			pStruct,
-			AB3History({zero3,zero3})))
-
 {
 
-}
+/// Range policy for integration kernel (alias)
+using rpIntegration = Kokkos::RangePolicy<
+		DefaultExecutionSpace,
+		Kokkos::Schedule<Kokkos::Static>,
+		Kokkos::IndexType<uint32>
+		>;
 
-bool pFlow::AdamsBashforth3::predict
+bool intAllActive(
-(
-	real UNUSED(dt),
-	realx3Vector_D& UNUSED(y),
-	realx3Vector_D& UNUSED(dy)
-)
-{
-
-	return true;
-}
-
-bool pFlow::AdamsBashforth3::correct
-(
 	real dt, 
-	realx3Vector_D& y,
+	realx3Field_D& y, 
-	realx3Vector_D& dy
+	realx3PointField_D& dy,
-)
+	realx3PointField_D& dy1,
+	realx3PointField_D& dy2,
+	real damping = 1.0)
 {
 
-	if(this->pStruct().allActive())
+	auto d_dy = dy.deviceView();
-	{
+	auto d_y  = y.deviceView();
-		return intAll(dt, y, dy, this->pStruct().activeRange());
+	auto d_dy1= dy1.deviceView();
-	}
+	auto d_dy2= dy2.deviceView();
-	else
+	auto activeRng = dy1.activeRange();
-	{
-		return intRange(dt, y, dy, this->pStruct().activePointsMaskD());
-	}
-
-	return true;
-}
-
-bool pFlow::AdamsBashforth3::setInitialVals(
-	const int32IndexContainer& newIndices,
-	const realx3Vector& y)
-{
-	return true;
-}
-
-bool pFlow::AdamsBashforth3::intAll(
-	real dt, 
-	realx3Vector_D& y, 
-	realx3Vector_D& dy, 
-	range activeRng)
-{
-	auto d_dy = dy.deviceViewAll();
-	auto d_y  = y.deviceViewAll();
-	auto d_history = history_.deviceViewAll();
 
 	Kokkos::parallel_for(
 		"AdamsBashforth3::correct",
-		rpIntegration (activeRng.first, activeRng.second),
+		rpIntegration (activeRng.start(), activeRng.end()),
-		LAMBDA_HD(int32 i){
+		LAMBDA_HD(uint32 i){
-			auto ldy = d_dy[i];
+			d_y[i] += damping* dt*( static_cast<real>(23.0 / 12.0) * d_dy[i] 
-			d_y[i] += dt*( static_cast<real>(23.0 / 12.0) * ldy 
+						- static_cast<real>(16.0 / 12.0) * d_dy1[i] 
-				- static_cast<real>(16.0 / 12.0) * d_history[i].dy1_ 
+						+ static_cast<real>(5.0 / 12.0) * d_dy2[i]) ;
-				+ static_cast<real>(5.0 / 12.0) * d_history[i].dy2_);
+			
-			d_history[i] = {ldy ,d_history[i].dy1_};
+			d_dy2[i] = d_dy1[i];
+			d_dy1[i] = d_dy[i];
+
 		});
 	Kokkos::fence();
 
 	return true;	
 }
+
+bool intScattered
+(
+	real dt, 
+	realx3Field_D& y,
+	realx3PointField_D& dy,
+	realx3PointField_D& dy1,
+	realx3PointField_D& dy2,
+	real damping = 1.0
+)
+{
+
+	auto d_dy 		= dy.deviceView();
+	auto d_y  		= y.deviceView();
+	auto d_dy1		= dy1.deviceView();
+	auto d_dy2		= dy2.deviceView();
+	auto activeRng 	= dy1.activeRange();
+	const auto& activeP = dy1.activePointsMaskDevice();
+
+	Kokkos::parallel_for(
+		"AdamsBashforth3::correct",
+		rpIntegration (activeRng.start(), activeRng.end()),
+		LAMBDA_HD(uint32 i){
+			if( activeP(i))
+			{
+				d_y[i] += damping * dt*( static_cast<real>(23.0 / 12.0) * d_dy[i] 
+						- static_cast<real>(16.0 / 12.0) * d_dy1[i] 
+						+ static_cast<real>(5.0 / 12.0) * d_dy2[i]);
+			
+				d_dy2[i] = d_dy1[i];
+				d_dy1[i] = d_dy[i];
+			}
+		});
+	Kokkos::fence();
+
+
+	return true;
+}
+
+}
+
+//const real AB3_coef[] = { 23.0 / 12.0, 16.0 / 12.0, 5.0 / 12.0 };
+
+pFlow::AdamsBashforth3::AdamsBashforth3
+(
+	const word& baseName,
+	pointStructure& pStruct,
+	const word& method,
+	const realx3Field_D& initialValField,
+	bool  keepHistory
+)
+:
+	AdamsBashforth2(baseName, pStruct, method, initialValField, keepHistory),
+	dy2_
+	(
+		objectFile
+		(
+			groupNames(baseName,"dy2"),
+			pStruct.time().integrationFolder(),
+			objectFile::READ_IF_PRESENT,
+			keepHistory ? objectFile::WRITE_ALWAYS : objectFile::WRITE_NEVER
+		),
+		pStruct,
+		zero3,
+		zero3
+	)
+
+{
+	
+}
+
+void pFlow::AdamsBashforth3::updateBoundariesSlaveToMasterIfRequested()
+{
+	AdamsBashforth2::updateBoundariesSlaveToMasterIfRequested();
+	dy2_.updateBoundariesSlaveToMasterIfRequested();
+}
+
+
+
+bool pFlow::AdamsBashforth3::correct
+(
+	real dt, 
+	realx3PointField_D& y, 
+	realx3PointField_D& dy,
+	real damping
+)
+{
+	
+	bool success = false;
+	if(y.isAllActive())
+	{
+		success = intAllActive(dt, y.field(), dy, dy1(), dy2(), damping);
+	}
+	else
+	{
+		success = intScattered(dt, y.field(), dy, dy1(), dy2(), damping);
+	}
+
+	success = success && boundaryList().correct(dt, y, dy);
+
+	return success;
+}
+
+bool pFlow::AdamsBashforth3::correctPStruct(
+	real dt, 
+	pointStructure &pStruct, 
+	realx3PointField_D &vel)
+{
+	
+	bool success = false;
+	if(dy2().isAllActive())
+	{
+		success = intAllActive(dt, pStruct.pointPosition(), vel, dy1(), dy2());
+	}
+	else
+	{
+		success = intScattered(dt, pStruct.pointPosition(), vel, dy1(), dy2());
+	}
+
+	success = success && boundaryList().correctPStruct(dt, pStruct, vel);
+
+	return success;
+}
+
+/*bool pFlow::AdamsBashforth3::hearChanges
+(
+	const timeInfo &ti, 
+	const message &msg, 
+	const anyList &varList
+)
+{
+	if(msg.equivalentTo(message::ITEMS_INSERT))
+	{
+		return insertValues(varList, initialValField().deviceViewAll(), dy1()) &&
+		insertValues(varList, initialValField().deviceViewAll(), dy2());
+	}
+	else
+	{
+		return realx3PointField_D::hearChanges(ti, msg, varList);
+	}
+}*/

src/Integration/AdamsBashforth3/AdamsBashforth3.hpp

@@ -22,48 +22,14 @@ Licence:
 #define __AdamsBashforth3_hpp__
 
 
-#include "integration.hpp"
+#include "AdamsBashforth2.hpp"
 #include "pointFields.hpp"
+#include "boundaryIntegrationList.hpp"
+
 
 namespace pFlow
 {
 
-struct AB3History
-{
-	TypeInfoNV("AB3History");
-
-	realx3 dy1_={0,0,0};
-	realx3 dy2_={0,0,0};
-};
-
-
-INLINE_FUNCTION
-iIstream& operator>>(iIstream& str, AB3History& ab3)
-{
-	str.readBegin("AB3History");
-
-    str >> ab3.dy1_;
-    str >> ab3.dy2_;
-    
-    str.readEnd("AB3History");
-
-    str.check(FUNCTION_NAME);
-
-	return str;
-	
-}
-
-INLINE_FUNCTION
-iOstream& operator<<(iOstream& str, const AB3History& ab3)
-{
-	str << token::BEGIN_LIST << ab3.dy1_
-	    << token::SPACE << ab3.dy2_
-		<< token::END_LIST;	    
-
-	str.check(FUNCTION_NAME);
-
-	return str;
-}
 
 /**
  * Third order Adams-Bashforth integration method for solving ODE
@@ -72,19 +38,26 @@ iOstream& operator<<(iOstream& str, const AB3History& ab3)
  */
 class AdamsBashforth3
 :
-	public integration
+	public AdamsBashforth2
 {
+private:
+	
+	friend class processorAB3BoundaryIntegration;
+
+	realx3PointField_D 		dy2_;
+
 protected:
 
-	/// Integration history
+	const auto& dy2()const
-	pointField<VectorSingle,AB3History>& history_;
+	{
+		return dy2_;
+	}
+
+	auto& dy2()
+	{
+		return dy2_;
+	}
 	
-	/// Range policy for integration kernel
-	using rpIntegration = Kokkos::RangePolicy<
-			DefaultExecutionSpace,
-			Kokkos::Schedule<Kokkos::Static>,
-			Kokkos::IndexType<int32>
-			>;
 
 public:
 
@@ -96,17 +69,14 @@ public:
 		/// Construct from components
 		AdamsBashforth3(
 			const word& baseName,
-			repository& owner,
+			pointStructure& pStruct,
-			const pointStructure& pStruct,
+			const word& method,
-			const word& method);
+			const realx3Field_D& initialValField,
+			bool  keepHistory);
 		
-		uniquePtr<integration> clone()const override
-		{
-			return makeUnique<AdamsBashforth3>(*this);
-		}
 
 		/// Destructor 
-		virtual ~AdamsBashforth3()=default;
+		~AdamsBashforth3() override =default;
 
 		/// Add this to the virtual constructor table 
 		add_vCtor(
@@ -117,43 +87,39 @@ public:
 
 	// - Methods
 
-		bool predict(
+		void updateBoundariesSlaveToMasterIfRequested()override;
-			real UNUSED(dt),
-			realx3Vector_D & UNUSED(y),
-			realx3Vector_D& UNUSED(dy)) override;
 
-		bool correct(real dt, 
+		/// return integration method 
-			realx3Vector_D & y,
+		word method()const override
-			realx3Vector_D& dy) override;
-
-		bool setInitialVals(
-			const int32IndexContainer& newIndices,
-			const realx3Vector& y) override;
-
-		bool needSetInitialVals()const override
 		{
-			return false;
+			return "AdamsBashforth3";
 		}
 
-		/// Integrate on all points in the active range 
-		bool intAll(
-			real dt, 
-			realx3Vector_D& y, 
-			realx3Vector_D& dy, 
-			range activeRng);
 		
-		/// Integrate on active points in the active range 
+
-		template<typename activeFunctor>
+		bool correct(
-		bool intRange(
 			real dt, 
-			realx3Vector_D& y,
+			realx3PointField_D& y, 
-			realx3Vector_D& dy,
+			realx3PointField_D& dy,
-			activeFunctor activeP );
+			real damping = 1.0) override;
+
+		bool correctPStruct(
+			real dt, 
+			pointStructure& pStruct, 
+			realx3PointField_D& vel) override;
+			
+		
+		/*bool hearChanges
+		(
+			const timeInfo& ti,
+			const message& msg, 
+			const anyList& varList
+		) override;*/
 
 };
 
 
-template<typename activeFunctor>
+/*template<typename activeFunctor>
 bool pFlow::AdamsBashforth3::intRange(
 	real dt, 
 	realx3Vector_D& y,
@@ -181,7 +147,7 @@ bool pFlow::AdamsBashforth3::intRange(
 	Kokkos::fence();
 
 	return true;
-}
+}*/
 
 } // pFlow
 

src/Integration/AdamsBashforth4/AdamsBashforth4.cpp

@@ -20,96 +20,187 @@ Licence:
 
 #include "AdamsBashforth4.hpp"
 
-
+namespace pFlow
-
-pFlow::AdamsBashforth4::AdamsBashforth4
-(
-	const word& baseName,
-	repository& owner,
-	const pointStructure& pStruct,
-	const word& method
-)
-:
-	integration(baseName, owner, pStruct, method),
-	history_(
-		owner.emplaceObject<pointField<VectorSingle,AB4History>>(
-			objectFile(
-				groupNames(baseName,"AB4History"),
-				"",
-				objectFile::READ_IF_PRESENT,
-				objectFile::WRITE_ALWAYS),
-			pStruct,
-			AB4History({zero3,zero3, zero3})))
-
 {
 
-}
+/// Range policy for integration kernel (alias)
+using rpIntegration = Kokkos::RangePolicy<
+		DefaultExecutionSpace,
+		Kokkos::Schedule<Kokkos::Static>,
+		Kokkos::IndexType<uint32>
+		>;
 
-bool pFlow::AdamsBashforth4::predict
+bool intAllActive(
-(
-	real UNUSED(dt),
-	realx3Vector_D& UNUSED(y),
-	realx3Vector_D& UNUSED(dy)
-)
-{
-
-	return true;
-}
-
-bool pFlow::AdamsBashforth4::correct
-(
 	real dt, 
-	realx3Vector_D& y,
+	realx3Field_D& y, 
-	realx3Vector_D& dy
+	realx3PointField_D& dy,
-)
+	realx3PointField_D& dy1,
+	realx3PointField_D& dy2,
+	realx3PointField_D& dy3,
+	real damping = 1.0)
 {
 
-	if(this->pStruct().allActive())
+	auto d_dy = dy.deviceView();
-	{
+	auto d_y  = y.deviceView();
-		return intAll(dt, y, dy, this->pStruct().activeRange());
+	auto d_dy1= dy1.deviceView();
-	}
+	auto d_dy2= dy2.deviceView();
-	else
+	auto d_dy3= dy3.deviceView();
-	{
+	auto activeRng = dy1.activeRange();
-		return intRange(dt, y, dy, this->pStruct().activePointsMaskD());
-	}
-
-	return true;
-}
-
-bool pFlow::AdamsBashforth4::setInitialVals(
-	const int32IndexContainer& newIndices,
-	const realx3Vector& y)
-{
-	return true;
-}
-
-bool pFlow::AdamsBashforth4::intAll(
-	real dt, 
-	realx3Vector_D& y, 
-	realx3Vector_D& dy, 
-	range activeRng)
-{
-	auto d_dy = dy.deviceViewAll();
-	auto d_y  = y.deviceViewAll();
-	auto d_history = history_.deviceViewAll();
 
 	Kokkos::parallel_for(
 		"AdamsBashforth4::correct",
-		rpIntegration (activeRng.first, activeRng.second),
+		rpIntegration (activeRng.start(), activeRng.end()),
-		LAMBDA_HD(int32 i){				
+		LAMBDA_HD(uint32 i){
-			d_y[i] += dt*( 
+			d_y[i] += damping * dt*( 
 				  static_cast<real>(55.0 / 24.0) * d_dy[i]
-					- static_cast<real>(59.0 / 24.0) * d_history[i].dy1_ 
+				- static_cast<real>(59.0 / 24.0) * d_dy1[i] 
-					+ static_cast<real>(37.0 / 24.0) * d_history[i].dy2_
+				+ static_cast<real>(37.0 / 24.0) * d_dy2[i]
-					- static_cast<real>( 9.0 / 24.0) * d_history[i].dy3_
+				- static_cast<real>( 9.0 / 24.0) * d_dy3[i]);
-					);
-			d_history[i].dy3_ = d_history[i].dy2_;
-			d_history[i].dy2_ = d_history[i].dy1_;
-			d_history[i].dy1_ = d_dy[i];
-			
 
+			d_dy3[i] = d_dy2[i];
+			d_dy2[i] = d_dy1[i];
+			d_dy1[i] = d_dy[i];
 		});
 	Kokkos::fence();
 
 	return true;	
 }
+
+bool intScattered
+(
+	real dt, 
+	realx3Field_D& y,
+	realx3PointField_D& dy,
+	realx3PointField_D& dy1,
+	realx3PointField_D& dy2,
+	realx3PointField_D& dy3,
+	real damping = 1.0
+)
+{
+
+	auto d_dy 		= dy.deviceView();
+	auto d_y  		= y.deviceView();
+	auto d_dy1		= dy1.deviceView();
+	auto d_dy2		= dy2.deviceView();
+	auto d_dy3		= dy3.deviceView();
+	auto activeRng 	= dy1.activeRange();
+	const auto& activeP = dy1.activePointsMaskDevice();
+
+	Kokkos::parallel_for(
+		"AdamsBashforth4::correct",
+		rpIntegration (activeRng.start(), activeRng.end()),
+		LAMBDA_HD(uint32 i){
+			if( activeP(i))
+			{
+				d_y[i] += damping*  dt*( 
+				  static_cast<real>(55.0 / 24.0) * d_dy[i]
+				- static_cast<real>(59.0 / 24.0) * d_dy1[i] 
+				+ static_cast<real>(37.0 / 24.0) * d_dy2[i]
+				- static_cast<real>( 9.0 / 24.0) * d_dy3[i] );
+				d_dy3[i] = d_dy2[i];
+				d_dy2[i] = d_dy1[i];
+				d_dy1[i] = d_dy[i];
+			}
+		});
+	Kokkos::fence();
+
+
+	return true;
+}
+
+}
+
+pFlow::AdamsBashforth4::AdamsBashforth4
+(
+	const word& baseName,
+	pointStructure& pStruct,
+	const word& method,
+	const realx3Field_D& initialValField,
+	bool keepHistory
+)
+:
+	AdamsBashforth3(baseName, pStruct, method, initialValField, keepHistory),
+	dy3_
+	(
+		objectFile
+		(
+			groupNames(baseName,"dy3"),
+			pStruct.time().integrationFolder(),
+			objectFile::READ_IF_PRESENT,
+			keepHistory?objectFile::WRITE_ALWAYS:objectFile::WRITE_NEVER
+		),
+		pStruct,
+		zero3,
+		zero3
+	)
+{}
+
+void pFlow::AdamsBashforth4::updateBoundariesSlaveToMasterIfRequested()
+{
+	AdamsBashforth3::updateBoundariesSlaveToMasterIfRequested();
+	dy3_.updateBoundariesSlaveToMasterIfRequested();
+}
+
+bool pFlow::AdamsBashforth4::correct
+(
+	real dt, 
+	realx3PointField_D& y, 
+	realx3PointField_D& dy,
+	real damping
+)
+{
+	
+	bool success = false;
+	if(y.isAllActive())
+	{
+		success = intAllActive(dt, y.field(), dy, dy1(), dy2(), dy3(), damping);
+	}
+	else
+	{
+		success = intScattered(dt, y.field(), dy, dy1(), dy2(), dy3(), damping);
+	}
+
+	success = success && boundaryList().correct(dt, y, dy);
+
+	return success;
+}
+
+bool pFlow::AdamsBashforth4::correctPStruct(
+	real dt, 
+	pointStructure &pStruct, 
+	realx3PointField_D &vel)
+{
+	
+	bool success = false;
+	if(dy2().isAllActive())
+	{
+		success = intAllActive(dt, pStruct.pointPosition(), vel, dy1(), dy2(), dy3());
+	}
+	else
+	{
+		success = intScattered(dt, pStruct.pointPosition(), vel, dy1(), dy2(), dy3());
+	}
+
+	success = success && boundaryList().correctPStruct(dt, pStruct, vel);
+
+	return success;
+}
+
+/*bool pFlow::AdamsBashforth4::hearChanges
+(
+	const timeInfo &ti, 
+	const message &msg, 
+	const anyList &varList
+)
+{
+	if(msg.equivalentTo(message::ITEMS_INSERT))
+	{
+		return insertValues(varList, initialValField().deviceViewAll(), dy1()) &&
+		insertValues(varList, initialValField().deviceViewAll(), dy2()) && 
+		insertValues(varList, initialValField().deviceViewAll(), dy3());
+	}
+	else
+	{
+		return realx3PointField_D::hearChanges(ti, msg, varList);
+	}
+}*/

src/Integration/AdamsBashforth4/AdamsBashforth4.hpp

@@ -22,53 +22,14 @@ Licence:
 #define __AdamsBashforth4_hpp__
 
 
-#include "integration.hpp"
+#include "AdamsBashforth3.hpp"
 #include "pointFields.hpp"
+#include "boundaryIntegrationList.hpp"
+
 
 namespace pFlow
 {
 
-struct AB4History
-{
-	TypeInfoNV("AB4History");
-
-	realx3 dy1_={0,0,0};
-	realx3 dy2_={0,0,0};
-	realx3 dy3_={0,0,0};
-	
-};
-
-
-INLINE_FUNCTION
-iIstream& operator>>(iIstream& str, AB4History& ab4)
-{
-	str.readBegin("AB4History");
-
-    str >> ab4.dy1_;
-    str >> ab4.dy2_;
-    str >> ab4.dy3_;
-    
-    str.readEnd("AB4History");
-
-    str.check(FUNCTION_NAME);
-
-	return str;
-	
-}
-
-INLINE_FUNCTION
-iOstream& operator<<(iOstream& str, const AB4History& ab4)
-{
-	str << token::BEGIN_LIST << ab4.dy1_
-	    << token::SPACE << ab4.dy2_
-	    << token::SPACE << ab4.dy3_
-		<< token::END_LIST;	    
-
-	str.check(FUNCTION_NAME);
-
-	return str;
-}
-
 /**
  * Fourth order Adams-Bashforth integration method for solving ODE
  * 
@@ -76,19 +37,25 @@ iOstream& operator<<(iOstream& str, const AB4History& ab4)
  */
 class AdamsBashforth4
 :
-	public integration
+	public AdamsBashforth3
 {
+private:
+
+	friend class processorAB4BoundaryIntegration;
+
+	realx3PointField_D 		dy3_;
+
 protected:
 
-	/// Integration history
+	const auto& dy3()const
-	pointField<VectorSingle,AB4History>& history_;
+	{
+		return dy3_;
+	}
 
-	/// Range policy for integration kernel
+	auto& dy3()
-	using rpIntegration = Kokkos::RangePolicy<
+	{
-			DefaultExecutionSpace,
+		return dy3_;
-			Kokkos::Schedule<Kokkos::Static>,
+	}
-			Kokkos::IndexType<int32>
-			>;
 
 public:
 
@@ -100,17 +67,15 @@ public:
 		/// Construct from components
 		AdamsBashforth4(
 			const word& baseName,
-			repository& owner,
+			pointStructure& pStruct,
-			const pointStructure& pStruct,
+			const word& method,
-			const word& method);
+			const realx3Field_D& initialValField,
+			bool  keepHistory);
+		
 		
-		uniquePtr<integration> clone()const override
-		{
-			return makeUnique<AdamsBashforth4>(*this);
-		}
 
 		/// Destructor
-		virtual ~AdamsBashforth4()=default;
+		~AdamsBashforth4() override =default;
 
 		/// Add a this to the virtual constructor table 
 		add_vCtor(
@@ -121,77 +86,36 @@ public:
 
 	// - Methods
 		
-		bool predict(
+		void updateBoundariesSlaveToMasterIfRequested()override;
-			real UNUSED(dt),
+
-			realx3Vector_D & UNUSED(y),
+		/// return integration method 
-			realx3Vector_D& UNUSED(dy)) override;
+		word method()const override
+		{
+			return "AdamsBashforth4";
+		}
 
 		bool correct(
 			real dt, 
-			realx3Vector_D & y,
+			realx3PointField_D& y, 
-			realx3Vector_D& dy) override;
+			realx3PointField_D& dy,
+			real damping = 1.0) override;
 
-		bool setInitialVals(
+		bool correctPStruct(
-			const int32IndexContainer& newIndices,
-			const realx3Vector& y) override;
-
-		bool needSetInitialVals()const override
-		{
-			return false;
-		}
-
-		/// Integrate on all points in the active range 
-		bool intAll(
 			real dt, 
-			realx3Vector_D& y, 
+			pointStructure& pStruct, 
-			realx3Vector_D& dy, 
+			realx3PointField_D& vel) override;
-			range activeRng);
 		
-		/// Integrate on active points in the active range 
+		/*bool hearChanges
-		template<typename activeFunctor>
+		(
-		bool intRange(
+			const timeInfo& ti,
-			real dt,
+			const message& msg, 
-			realx3Vector_D& y,
+			const anyList& varList
-			realx3Vector_D& dy,
+		) override;*/
-			activeFunctor activeP );
 
 };
 
 
-template<typename activeFunctor>
-bool pFlow::AdamsBashforth4::intRange(
-	real dt, 
-	realx3Vector_D& y,
-	realx3Vector_D& dy,
-	activeFunctor activeP )
-{
-	auto d_dy = dy.deviceViewAll();
-	auto d_y  = y.deviceViewAll();
-	auto d_history = history_.deviceViewAll();
-	auto activeRng = activeP.activeRange();
 
-	Kokkos::parallel_for(
-		"AdamsBashforth4::correct",
-		rpIntegration (activeRng.first, activeRng.second),
-		LAMBDA_HD(int32 i){
-			if( activeP(i))
-			{
-				
-				d_y[i] += dt*( 
-							static_cast<real>(55.0 / 24.0) * d_dy[i]
-						- static_cast<real>(59.0 / 24.0) * d_history[i].dy1_ 
-						+ static_cast<real>(37.0 / 24.0) * d_history[i].dy2_
-						- static_cast<real>( 9.0 / 24.0) * d_history[i].dy3_
-						);
-				d_history[i].dy3_ = d_history[i].dy2_;
-				d_history[i].dy2_ = d_history[i].dy1_;
-				d_history[i].dy1_ = d_dy[i];
-			}
-		});
-	Kokkos::fence();
-
-	return true;
-}
 
 } // pFlow
 

src/Integration/AdamsBashforth5/AdamsBashforth5.cpp

@@ -20,93 +20,179 @@ Licence:
 
 #include "AdamsBashforth5.hpp"
 
-
+namespace pFlow
-
-pFlow::AdamsBashforth5::AdamsBashforth5
-(
-	const word& baseName,
-	repository& owner,
-	const pointStructure& pStruct,
-	const word& method
-)
-:
-	integration(baseName, owner, pStruct, method),
-	history_(
-		owner.emplaceObject<pointField<VectorSingle,AB5History>>(
-			objectFile(
-				groupNames(baseName,"AB5History"),
-				"",
-				objectFile::READ_IF_PRESENT,
-				objectFile::WRITE_ALWAYS),
-			pStruct,
-			AB5History({zero3,zero3, zero3})))
-
 {
 
-}
+/// Range policy for integration kernel (alias)
+using rpIntegration = Kokkos::RangePolicy<
+		DefaultExecutionSpace,
+		Kokkos::Schedule<Kokkos::Static>,
+		Kokkos::IndexType<uint32>
+		>;
 
-bool pFlow::AdamsBashforth5::predict
+bool intAllActive(
-(
-	real UNUSED(dt),
-	realx3Vector_D& UNUSED(y),
-	realx3Vector_D& UNUSED(dy)
-)
-{
-
-	return true;
-}
-
-bool pFlow::AdamsBashforth5::correct
-(
 	real dt, 
-	realx3Vector_D& y,
+	realx3Field_D& y, 
-	realx3Vector_D& dy
+	realx3PointField_D& dy,
-)
+	realx3PointField_D& dy1,
+	realx3PointField_D& dy2,
+	realx3PointField_D& dy3,
+	realx3PointField_D& dy4,
+	real damping = 1.0)
 {
 
-	if(this->pStruct().allActive())
+	auto d_dy = dy.deviceView();
-	{
+	auto d_y  = y.deviceView();
-		return intAll(dt, y, dy, this->pStruct().activeRange());
+	auto d_dy1= dy1.deviceView();
-	}
+	auto d_dy2= dy2.deviceView();
-	else
+	auto d_dy3= dy3.deviceView();
-	{
+	auto d_dy4= dy4.deviceView();
-		return intRange(dt, y, dy, this->pStruct().activePointsMaskD());
+	auto activeRng = dy1.activeRange();
-	}
-
-	return true;
-}
-
-bool pFlow::AdamsBashforth5::setInitialVals(
-	const int32IndexContainer& newIndices,
-	const realx3Vector& y)
-{
-	return true;
-}
-
-bool pFlow::AdamsBashforth5::intAll(
-	real dt, 
-	realx3Vector_D& y, 
-	realx3Vector_D& dy, 
-	range activeRng)
-{
-	auto d_dy = dy.deviceViewAll();
-	auto d_y  = y.deviceViewAll();
-	auto d_history = history_.deviceViewAll();
 
 	Kokkos::parallel_for(
 		"AdamsBashforth5::correct",
-		rpIntegration (activeRng.first, activeRng.second),
+		rpIntegration (activeRng.start(), activeRng.end()),
-		LAMBDA_HD(int32 i){				
+		LAMBDA_HD(uint32 i){
-			d_y[i] += dt*( 
+			d_y[i] += damping*  dt*( 
 				  static_cast<real>(1901.0 / 720.0) * d_dy[i]
-				- static_cast<real>(2774.0 / 720.0) * d_history[i].dy1_ 
+				- static_cast<real>(2774.0 / 720.0) * d_dy1[i]
-				+ static_cast<real>(2616.0 / 720.0) * d_history[i].dy2_
+				+ static_cast<real>(2616.0 / 720.0) * d_dy2[i]
-				- static_cast<real>(1274.0 / 720.0) * d_history[i].dy3_
+				- static_cast<real>(1274.0 / 720.0) * d_dy3[i]
-				+ static_cast<real>( 251.0 / 720.0) * d_history[i].dy4_
+				+ static_cast<real>( 251.0 / 720.0) * d_dy4[i]);
-				);
+				
-			d_history[i] = {d_dy[i] ,d_history[i].dy1_, d_history[i].dy2_, d_history[i].dy3_};
+			d_dy4[i] = d_dy3[i];
+			d_dy3[i] = d_dy2[i];
+			d_dy2[i] = d_dy1[i];
+			d_dy1[i] = d_dy[i];
 		});
 	Kokkos::fence();
 
 	return true;	
 }
+
+bool intScattered
+(
+	real dt, 
+	realx3Field_D& y,
+	realx3PointField_D& dy,
+	realx3PointField_D& dy1,
+	realx3PointField_D& dy2,
+	realx3PointField_D& dy3,
+	realx3PointField_D& dy4,
+	real damping = 1.0
+)
+{
+
+	auto d_dy 		= dy.deviceView();
+	auto d_y  		= y.deviceView();
+	auto d_dy1		= dy1.deviceView();
+	auto d_dy2		= dy2.deviceView();
+	auto d_dy3		= dy3.deviceView();
+	auto d_dy4		= dy4.deviceView();
+	auto activeRng 	= dy1.activeRange();
+	const auto& activeP = dy1.activePointsMaskDevice();
+
+	Kokkos::parallel_for(
+		"AdamsBashforth5::correct",
+		rpIntegration (activeRng.start(), activeRng.end()),
+		LAMBDA_HD(uint32 i){
+			if( activeP(i))
+			{
+				d_y[i] += damping*  dt*( 
+				  static_cast<real>(1901.0 / 720.0) * d_dy[i]
+				- static_cast<real>(2774.0 / 720.0) * d_dy1[i]
+				+ static_cast<real>(2616.0 / 720.0) * d_dy2[i]
+				- static_cast<real>(1274.0 / 720.0) * d_dy3[i]
+				+ static_cast<real>( 251.0 / 720.0) * d_dy4[i]);
+
+				d_dy4[i] = d_dy3[i];
+				d_dy3[i] = d_dy2[i];
+				d_dy2[i] = d_dy1[i];
+				d_dy1[i] = d_dy[i];
+			}
+		});
+	Kokkos::fence();
+
+	return true;
+}
+
+}
+
+pFlow::AdamsBashforth5::AdamsBashforth5
+(
+	const word &baseName, 
+	pointStructure &pStruct, 
+	const word &method, 
+	const realx3Field_D &initialValField,
+	bool  keepHistory
+)
+:
+	AdamsBashforth4(baseName, pStruct, method, initialValField, keepHistory),
+	dy4_
+	(
+		objectFile
+		(
+			groupNames(baseName,"dy4"),
+			pStruct.time().integrationFolder(),
+			objectFile::READ_IF_PRESENT,
+			keepHistory?objectFile::WRITE_ALWAYS:objectFile::WRITE_NEVER
+		),
+		pStruct,
+		zero3,
+		zero3
+	)
+{}
+
+void pFlow::AdamsBashforth5::updateBoundariesSlaveToMasterIfRequested()
+{
+	AdamsBashforth4::updateBoundariesSlaveToMasterIfRequested();
+	dy4_.updateBoundariesSlaveToMasterIfRequested();
+}
+
+bool pFlow::AdamsBashforth5::correct
+(
+	real dt, 
+	realx3PointField_D &y, 
+	realx3PointField_D &dy, 
+	real damping
+)
+{
+	
+	bool success = false;
+	if(y.isAllActive())
+	{
+		success = intAllActive(dt, y.field(), dy, dy1(), dy2(), dy3(), dy4(), damping);
+	}
+	else
+	{
+		success = intScattered(dt, y.field(), dy, dy1(), dy2(), dy3(), dy4(), damping);
+	}
+
+	success = success && boundaryList().correct(dt, y, dy);
+
+	return success;
+}
+
+bool pFlow::AdamsBashforth5::correctPStruct
+(
+	real dt, 
+	pointStructure &pStruct, 
+	realx3PointField_D &vel
+)
+{
+bool success = false;
+	if(dy2().isAllActive())
+	{
+		success = intAllActive(dt, pStruct.pointPosition(), vel, dy1(), dy2(), dy3(), dy4());
+	}
+	else
+	{
+		success = intScattered(dt, pStruct.pointPosition(), vel, dy1(), dy2(), dy3(), dy4());
+	}
+
+	success = success && boundaryList().correctPStruct(dt, pStruct, vel);
+
+	return success;
+}
+
+

src/Integration/AdamsBashforth5/AdamsBashforth5.hpp

@@ -22,54 +22,12 @@ Licence:
 #define __AdamsBashforth5_hpp__
 
 
-#include "integration.hpp"
+#include "AdamsBashforth4.hpp"
 #include "pointFields.hpp"
 
 namespace pFlow
 {
 
-struct AB5History
-{
-	TypeInfoNV("AB5History");
-
-	realx3 dy1_={0,0,0};
-	realx3 dy2_={0,0,0};
-	realx3 dy3_={0,0,0};
-	realx3 dy4_={0,0,0};
-};
-
-
-INLINE_FUNCTION
-iIstream& operator>>(iIstream& str, AB5History& ab5)
-{
-	str.readBegin("AB5History");
-
-    str >> ab5.dy1_;
-    str >> ab5.dy2_;
-    str >> ab5.dy3_;
-    str >> ab5.dy4_;
-    
-    str.readEnd("AB5History");
-
-    str.check(FUNCTION_NAME);
-
-	return str;
-	
-}
-
-INLINE_FUNCTION
-iOstream& operator<<(iOstream& str, const AB5History& ab5)
-{
-	str << token::BEGIN_LIST << ab5.dy1_
-	    << token::SPACE << ab5.dy2_
-	    << token::SPACE << ab5.dy3_
-	    << token::SPACE << ab5.dy4_
-		<< token::END_LIST;	    
-
-	str.check(FUNCTION_NAME);
-
-	return str;
-}
 
 /**
  * Fifth order Adams-Bashforth integration method for solving ODE
@@ -78,19 +36,25 @@ iOstream& operator<<(iOstream& str, const AB5History& ab5)
  */
 class AdamsBashforth5
 :
-	public integration
+	public AdamsBashforth4
 {
+private:
+
+	friend class processorAB4BoundaryIntegration;
+
+	realx3PointField_D 		dy4_;
+
 protected:
 
-	/// Integration history
+	const auto& dy4()const
-	pointField<VectorSingle,AB5History>& history_;
+	{
+		return dy4_;
+	}
 
-	/// Range policy for integration kernel
+	auto& dy4()
-	using rpIntegration = Kokkos::RangePolicy<
+	{
-			DefaultExecutionSpace,
+		return dy4_;
-			Kokkos::Schedule<Kokkos::Static>,
+	}
-			Kokkos::IndexType<int32>
-			>;
 
 public:
 
@@ -99,20 +63,20 @@ public:
 
 	// - Constructors
 		
+		/// Construct from components
 		AdamsBashforth5(
 			const word& baseName,
-			repository& owner,
+			pointStructure& pStruct,
-			const pointStructure& pStruct,
+			const word& method,
-			const word& method);
+			const realx3Field_D& initialValField,
+			bool  keepHistory);
 		
-		uniquePtr<integration> clone()const override
-		{
-			return makeUnique<AdamsBashforth5>(*this);
-		}
 		
-		virtual ~AdamsBashforth5()=default;
 
-		/// Add this to the virtual constructor table 
+		/// Destructor
+		~AdamsBashforth5() override =default;
+
+		/// Add a this to the virtual constructor table 
 		add_vCtor(
 			integration,
 			AdamsBashforth5,
@@ -121,44 +85,29 @@ public:
 
 	// - Methods
 		
-		bool predict(
+		void updateBoundariesSlaveToMasterIfRequested()override;
-			real UNUSED(dt),
+
-			realx3Vector_D & UNUSED(y),
+		/// return integration method 
-			realx3Vector_D& UNUSED(dy)) override;
+		word method()const override
+		{
+			return "AdamsBashforth5";
+		}
 
 		bool correct(
 			real dt, 
-			realx3Vector_D & y,
+			realx3PointField_D& y, 
-			realx3Vector_D& dy) override;
+			realx3PointField_D& dy,
+			real damping = 1.0) override;
 
-		bool setInitialVals(
+		bool correctPStruct(
-			const int32IndexContainer& newIndices,
-			const realx3Vector& y) override;
-
-		bool needSetInitialVals()const override
-		{
-			return false;
-		}
-
-		/// Integrate on all points in the active range 
-		bool intAll(
 			real dt, 
-			realx3Vector_D& y, 
+			pointStructure& pStruct, 
-			realx3Vector_D& dy, 
+			realx3PointField_D& vel) override;
-			range activeRng);
-
-		/// Integrate on active points in the active range 
-		template<typename activeFunctor>
-		bool intRange(
-			real dt,
-			realx3Vector_D& y,
-			realx3Vector_D& dy,
-			activeFunctor activeP );
 
 };
 
 
-template<typename activeFunctor>
+/*template<typename activeFunctor>
 bool pFlow::AdamsBashforth5::intRange(
 	real dt, 
 	realx3Vector_D& y,
@@ -190,7 +139,7 @@ bool pFlow::AdamsBashforth5::intRange(
 	Kokkos::fence();
 
 	return true;
-}
+}*/
 
 } // pFlow
 

src/Integration/CMakeLists.txt

@@ -4,9 +4,10 @@ integration/integration.cpp
 boundaries/boundaryIntegration.cpp
 boundaries/boundaryIntegrationList.cpp 
 AdamsBashforth2/AdamsBashforth2.cpp
-#AdamsBashforth5/AdamsBashforth5.cpp
+AdamsBashforth3/AdamsBashforth3.cpp
-#AdamsBashforth4/AdamsBashforth4.cpp
+AdamsBashforth4/AdamsBashforth4.cpp
-#AdamsBashforth3/AdamsBashforth3.cpp
+AdamsBashforth5/AdamsBashforth5.cpp
+
 #AdamsMoulton3/AdamsMoulton3.cpp
 #AdamsMoulton4/AdamsMoulton4.cpp
 #AdamsMoulton5/AdamsMoulton5.cpp

src/Integration/boundaries/boundaryIntegration.hpp

@@ -51,9 +51,7 @@ public:
     );
 
     bool hearChanges(
-        real t,
+        const timeInfo& ti,
-        real dt,
-        uint32 iter,
         const message &msg,
         const anyList &varList) override
     {
@@ -77,6 +75,11 @@ public:
         return true;
     }
 
+    bool isActive()const override
+    {
+        return false;
+    }
+    
     static
     uniquePtr<boundaryIntegration> create(
         const boundaryBase& boundary,

src/Integration/boundaries/boundaryIntegrationList.cpp

@@ -7,11 +7,11 @@ pFlow::boundaryIntegrationList::boundaryIntegrationList(
     integration &intgrtn
 )
 :
-    ListPtr<boundaryIntegration>(6),
+    boundaryListPtr<boundaryIntegration>(),
 	boundaries_(pStruct.boundaries())
 {
 	
-    for(uint32 i=0; i<6; i++)
+	ForAllBoundariesPtr(i, this)
 	{
 		this->set(
 			i,
@@ -23,20 +23,18 @@ pFlow::boundaryIntegrationList::boundaryIntegrationList(
 			)
 		);
 	}	
-
 }
 
 bool pFlow::boundaryIntegrationList::correct(real dt, realx3PointField_D &y, realx3PointField_D &dy)
 {
-	for(auto& bndry:*this)
+	ForAllActiveBoundariesPtr(i,this)
 	{
-		if(!bndry->correct(dt, y, dy))
+		if(!boundaryPtr(i)->correct(dt, y, dy))
 		{
 			fatalErrorInFunction<<"Error in correcting boundary "<<
-			bndry->boundaryName()<<endl;
+			boundaryPtr(i)->boundaryName()<<endl;
 			return false;
 		}
-		
 	}
     return true;
 }
@@ -46,14 +44,15 @@ bool pFlow::boundaryIntegrationList::correctPStruct(
 	pointStructure &pStruct, 
 	const realx3PointField_D &vel)
 {
-	for(auto& bndry:*this)
+	ForAllActiveBoundariesPtr(i,this)
 	{
-		if(!bndry->correctPStruct(dt, vel))
+		if(!boundaryPtr(i)->correctPStruct(dt, vel))
 		{
 			fatalErrorInFunction<<"Error in correcting boundary "<<
-			bndry->boundaryName()<<" in pointStructure."<<endl;
+			boundaryPtr(i)->boundaryName()<<" in pointStructure."<<endl;
 			return false;
 		}
 	}
+	
     return true;
 }

src/Integration/boundaries/boundaryIntegrationList.hpp

@@ -4,7 +4,7 @@
 
 
 #include "boundaryList.hpp"
-#include "ListPtr.hpp"
+#include "boundaryListPtr.hpp"
 #include "boundaryIntegration.hpp"
 
 
@@ -15,7 +15,7 @@ class integration;
 
 class boundaryIntegrationList
 :
-    public ListPtr<boundaryIntegration>
+    public boundaryListPtr<boundaryIntegration>
 {
 private:
 

src/Integration/integration/integration.cpp

@@ -22,17 +22,41 @@ Licence:
 #include "pointStructure.hpp"
 #include "repository.hpp"
 
-pFlow::integration::integration
+bool pFlow::integration::insertValues
 (
-	const word& baseName,
+	const anyList& varList, 
-	pointStructure& pStruct,
+	const deviceViewType1D<realx3>& srcVals,
-	const word&,
+	realx3PointField_D& dstFeild
-	const realx3Field_D&
 )
-:
+{
-	owner_(*pStruct.owner()),
+	const word eventName = message::eventName(message::ITEMS_INSERT);
+
+	if(!varList.checkObjectType<uint32IndexContainer>(eventName))
+	{
+		fatalErrorInFunction<<"Insertion failed in integration for "<< baseName_<<
+		", variable "<< eventName <<
+		" does not exist or the type is incorrect"<<endl;
+		return false;
+	}
+
+	const auto& indices = varList.getObject<uint32IndexContainer>(
+		eventName);
+	
+	dstFeild.field().insertSetElement(indices, srcVals);
+
+    return true;
+}
+
+pFlow::integration::integration(
+    const word &baseName,
+    pointStructure &pStruct,
+    const word &,
+    const realx3Field_D &,
+	bool keepHistory)
+    : owner_(*pStruct.owner()),
       pStruct_(pStruct),
-	baseName_(baseName)	
+      baseName_(baseName),
+      keepHistory_(keepHistory)
 {}
 
 
@@ -42,12 +66,13 @@ pFlow::uniquePtr<pFlow::integration>
 		const word& baseName,
 		pointStructure& pStruct,
 		const word& method,
-		const realx3Field_D& initialValField
+		const realx3Field_D& initialValField,
+		bool  keepHistory
 )
 {
 	if( wordvCtorSelector_.search(method) )
 	{
-		return wordvCtorSelector_[method] (baseName, pStruct, method, initialValField);
+		return wordvCtorSelector_[method] (baseName, pStruct, method, initialValField, keepHistory);
 	}
 	else
 	{

src/Integration/integration/integration.hpp

@@ -24,6 +24,7 @@ Licence:
 
 #include "virtualConstructor.hpp"
 #include "pointFields.hpp"
+#include "Logical.hpp"
 
 
 namespace pFlow
@@ -63,6 +64,15 @@ private:
 		/// The base name for integration 
 		const word baseName_;
 
+		bool 	keepHistory_;
+
+protected:
+
+	bool insertValues(
+		const anyList& varList, 
+		const deviceViewType1D<realx3>& srcVals,
+		realx3PointField_D& dstFeild);
+
 public:
 
 	/// Type info
@@ -76,7 +86,8 @@ public:
 			const word& baseName,
 			pointStructure& pStruct,
 			const word& method,
-			const realx3Field_D& initialValField);
+			const realx3Field_D& initialValField,
+			bool  keepHistory);
 
 		/// Copy constructor 
 		integration(const integration&) = default;
@@ -102,9 +113,10 @@ public:
 				const word& baseName,
 				pointStructure& pStruct,
 				const word& method,
-				const realx3Field_D& initialValField
+				const realx3Field_D& initialValField,
+				bool keepHistory
 			),
-			(baseName, pStruct, method, initialValField)
+			(baseName, pStruct, method, initialValField, keepHistory)
 		);
 
 
@@ -131,6 +143,11 @@ public:
 			return owner_;
 		}
 
+		bool keepHistory()const
+		{
+			return keepHistory_;
+		}
+
 		virtual 
 		void updateBoundariesSlaveToMasterIfRequested() = 0;
 		/// return integration method 
@@ -146,29 +163,19 @@ public:
 
 		/// Correction/main integration step
 		virtual 
-		bool correct(real dt, realx3PointField_D& y, realx3PointField_D& dy) = 0;
+		bool correct(real dt, realx3PointField_D& y, realx3PointField_D& dy, real damping = 1.0) = 0;
 
 		virtual 
 		bool correctPStruct(real dt, pointStructure& pStruct, realx3PointField_D& vel) = 0;
 
-		/// Set the initial values for new indices 
-		virtual 
-		bool setInitialVals(
-				const int32IndexContainer& newIndices,
-				const realx3Vector& y) = 0;
-
-		/// Check if the method requires any set initial vals
-		virtual 
-		bool needSetInitialVals()const = 0;
-
-
 	/// Create the polymorphic object based on inputs
 	static
 		uniquePtr<integration> create(
 			const word& baseName,
 			pointStructure& pStruct,
 			const word& method,
-			const realx3Field_D& initialValField);
+			const realx3Field_D& initialValField,
+			bool  keepHistory);
 
 }; // integration
 

src/Interaction/Models/contactForce/cGAbsoluteLinearCF.hpp

@@ -25,12 +25,6 @@ Licence:
 #include "symArrays.hpp"
 
 
-
-
-
-
-
-
 namespace pFlow::cfModels
 {
 
@@ -159,7 +153,7 @@ protected:
 		{
 
 			etha_t[i] = -2.0*log( et[i])*sqrt(kt[i]*2/7) /
-					sqrt(log(pow(et[i],2.0))+ pow(Pi,2.0));
+					sqrt(log(pow(et[i],2))+ pow(Pi,2));
 		}
 
 		Vector<linearProperties> prop("prop", nElem);
@@ -285,8 +279,8 @@ public:
 
 		history.overlap_t_ += Vt*dt;
 
-		real mi = 3*Pi/4*pow(Ri,3.0)*rho_[propId_i];
+		real mi = 3*Pi/4*pow(Ri,3)*rho_[propId_i];
-		real mj = 3*Pi/4*pow(Rj,3.0)*rho_[propId_j];
+		real mj = 3*Pi/4*pow(Rj,3)*rho_[propId_j];
 
 		real sqrt_meff = sqrt((mi*mj)/(mi+mj));
 
@@ -299,29 +293,24 @@ public:
 		else if (addDissipationModel_==3)
 		{
 			auto pie =3.14;
-			prop.en_ = exp((pow(f_,1.5)*log(prop.en_)*sqrt( (1-((pow(log(prop.en_),2))/(pow(log(prop.en_),2)+pow(pie,2))))/(1-(pow(f_,3)*(pow(log(prop.en_),2))/(pow(log(prop.en_),2)+pow(pie,2)))) ) ));
+			prop.en_ = exp((pow(f_,static_cast<real>(1.5))*log(prop.en_)*sqrt( (1-((pow(log(prop.en_),2))/(pow(log(prop.en_),2)+pow(pie,2))))/(1-(pow(f_,3)*(pow(log(prop.en_),2))/(pow(log(prop.en_),2)+pow(pie,2)))) ) ));
 		}
 
-
-
-
 		real ethan_  = -2.0*log(prop.en_)*sqrt(prop.kn_)/
-					sqrt(pow(log(prop.en_),2.0)+ pow(Pi,2.0));
+					sqrt(pow(log(prop.en_),2)+ pow(Pi,2));
 
 			//REPORT(0)<<"\n en n is  :  "<<END_REPORT;
 			//REPORT(0)<< prop.en_ <<END_REPORT;
 
-		FCn = (  -pow(f_,3.0)*prop.kn_ * ovrlp_n - sqrt_meff * pow(f_,1.5) * ethan_ * vrn)*Nij;
+		FCn = (  -pow(f_,3)*prop.kn_ * ovrlp_n - sqrt_meff * pow(f_,static_cast<real>(1.5)) * ethan_ * vrn)*Nij;
-		FCt = ( -pow(f_,3.0)*prop.kt_ * history.overlap_t_ - sqrt_meff * pow(f_,1.5) *  prop.ethat_*Vt);
+		FCt = ( -pow(f_,3)*prop.kt_ * history.overlap_t_ - sqrt_meff * pow(f_,static_cast<real>(1.5)) *  prop.ethat_*Vt);
-		
-
 				
 		real ft = length(FCt);
 		real ft_fric = prop.mu_ * length(FCn);
 		
 		if(ft > ft_fric)
 		{
-			if( length(history.overlap_t_) >static_cast<real>(0.0))
+			if( length(history.overlap_t_) >zero)
 			{
 				if constexpr (limited)
 				{

src/Interaction/Models/contactForce/cGNonLinearCF.hpp

@@ -0,0 +1,307 @@
+/*------------------------------- phasicFlow ---------------------------------
+      O        C enter of
+     O O       E ngineering and
+    O   O      M ultiscale modeling of
+   OOOOOOO     F luid flow       
+------------------------------------------------------------------------------
+  Copyright (C): www.cemf.ir
+  email: hamid.r.norouzi AT gmail.com
+------------------------------------------------------------------------------  
+Licence:
+  This file is part of phasicFlow code. It is a free software for simulating 
+  granular and multiphase flows. You can redistribute it and/or modify it under
+  the terms of GNU General Public License v3 or any other later versions. 
+ 
+  phasicFlow is distributed to help others in their research in the field of 
+  granular and multiphase flows, but WITHOUT ANY WARRANTY; without even the
+  implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
+
+-----------------------------------------------------------------------------*/
+
+#ifndef __cGNonLinearCF_hpp__
+#define __cGNonLinearCF_hpp__
+
+#include "types.hpp"
+
+namespace pFlow::cfModels
+{
+
+template<bool limited=true>
+class cGNonLinear
+{
+public:
+
+	struct contactForceStorage
+	{
+		realx3 overlap_t_ = 0.0;
+	};
+
+	struct cGNonLinearProperties
+	{		
+		real 	Yeff_	= 1000000.0;
+		real 	Geff_	= 8000000.0;
+		real 	en_ 	= 1.0;
+		real  	mu_ 	= 0.00001;
+
+		INLINE_FUNCTION_HD
+		cGNonLinearProperties(){}
+
+		INLINE_FUNCTION_HD
+		cGNonLinearProperties(real Yeff, real Geff, real en, real mu ):
+			Yeff_(Yeff), Geff_(Geff), en_(en), mu_(mu)
+		{}		
+
+		INLINE_FUNCTION_HD
+		cGNonLinearProperties(const cGNonLinearProperties&)=default;
+
+		INLINE_FUNCTION_HD
+		cGNonLinearProperties& operator=(const cGNonLinearProperties&)=default;
+
+		INLINE_FUNCTION_HD
+		~cGNonLinearProperties() = default;
+	};
+
+protected:
+
+	using cGNonLinearArrayType = 	symArray<cGNonLinearProperties>;
+
+	int32 						numMaterial_ = 0;
+
+	ViewType1D<real>  		 	rho_;
+
+	int32 						addDissipationModel_;
+
+	cGNonLinearArrayType  		nonlinearProperties_;
+
+	bool readNonLinearDictionary(const dictionary& dict)
+	{
+		auto Yeff = dict.getVal<realVector>("Yeff");
+		auto Geff = dict.getVal<realVector>("Geff");
+		auto nu = dict.getVal<realVector>("nu");
+		auto en = dict.getVal<realVector>("en");
+		auto mu = dict.getVal<realVector>("mu");
+
+		auto nElem = Yeff.size();
+
+		if(nElem != nu.size())
+		{
+			fatalErrorInFunction<<
+			"sizes of Yeff("<<nElem<<") and nu("<<nu.size()<<") do not match.\n";
+			return false;
+		}
+
+		if(nElem != en.size())
+		{
+			fatalErrorInFunction<<
+			"sizes of Yeff("<<nElem<<") and en("<<en.size()<<") do not match.\n";
+			return false;
+		}
+
+
+		if(nElem != mu.size())
+		{
+			fatalErrorInFunction<<
+			"sizes of Yeff("<<nElem<<") and mu("<<mu.size()<<") do not match.\n";
+			return false;
+		}
+
+		
+		// check if size of vector matchs a symetric array
+		uint32 nMat;
+		if( !cGNonLinearArrayType::getN(nElem, nMat) )
+		{
+			fatalErrorInFunction<<
+			"sizes of properties do not match a symetric array with size ("<<
+			numMaterial_<<"x"<<numMaterial_<<").\n";
+			return false;
+		}
+		else if( numMaterial_ != nMat)
+		{
+			fatalErrorInFunction<<
+			"size mismatch for porperties. \n";
+			return false;
+		}
+
+		Vector<cGNonLinearProperties> prop("prop",nElem);
+		ForAll(i,Yeff)
+		{
+			prop[i] = {Yeff[i], Geff[i], en[i], mu[i]};
+		}
+
+		nonlinearProperties_.assign(prop);
+
+		auto adm = dict.getVal<word>("additionalDissipationModel");
+		if(adm == "none")
+		{
+			addDissipationModel_ = 1;
+		}
+		else if(adm == "LU")
+		{
+			addDissipationModel_ = 2;
+		}
+		else if (adm == "GB")
+		{
+			addDissipationModel_ = 3;
+		}
+		else
+		{
+			addDissipationModel_ = 1;
+		}
+
+		return true;
+
+	}
+
+	static const char* modelName()
+	{
+		if constexpr (limited)
+		{
+			return "cGNonLinearLimited";
+		}
+		else
+		{
+			return "cGNonLinearNonLimited";
+		}
+		return "";
+	}
+
+public:
+
+	TypeInfoNV(modelName());
+
+	INLINE_FUNCTION_HD
+	cGNonLinear(){}
+
+	cGNonLinear(
+		int32 nMaterial,
+		const ViewType1D<real>& rho,
+		const dictionary& dict)
+	:
+		numMaterial_(nMaterial),
+		rho_("rho",nMaterial),
+		nonlinearProperties_("cGNonLinearProperties",nMaterial)
+	{
+
+		Kokkos::deep_copy(rho_,rho);
+		if(!readNonLinearDictionary(dict))
+		{
+			fatalExit;
+		}
+	}
+	
+	INLINE_FUNCTION_HD
+	cGNonLinear(const cGNonLinear&) = default;
+
+	INLINE_FUNCTION_HD
+	cGNonLinear(cGNonLinear&&) = default;
+
+	INLINE_FUNCTION_HD
+	cGNonLinear& operator=(const cGNonLinear&) = default;
+
+	INLINE_FUNCTION_HD
+	cGNonLinear& operator=(cGNonLinear&&) = default;
+
+
+	INLINE_FUNCTION_HD
+	~cGNonLinear()=default;
+
+	INLINE_FUNCTION_HD
+	int32 numMaterial()const
+	{
+		return numMaterial_;
+	}
+
+	//// - Methods
+
+	INLINE_FUNCTION_HD
+	void contactForce
+	(
+		const real dt,
+		const uint32 i,
+		const uint32 j,
+		const uint32 propId_i,
+		const uint32 propId_j,
+		const real Ri,
+		const real Rj,
+		const real cGFi,
+		const real cGFj,
+		const real ovrlp_n,
+		const realx3& Vr,
+		const realx3& Nij,
+		contactForceStorage& history,
+		realx3& FCn,
+		realx3& FCt
+	)const
+	{
+
+		const auto prop = nonlinearProperties_(propId_i,propId_j);
+
+		const real f = 2/( 1/cGFi + 1/cGFj );
+
+		real vrn = dot(Vr, Nij);	
+		realx3 Vt = Vr - vrn*Nij;
+
+		history.overlap_t_ += Vt*dt;
+
+		real mi = 3*Pi/4*pow(Ri,static_cast<real>(3))*rho_[propId_i];
+		real mj = 3*Pi/4*pow(Rj,static_cast<real>(3))*rho_[propId_j];
+		real Reff = 1/(1/Ri + 1/Rj);
+
+		real K_hertz = 4.0/3.0*prop.Yeff_*sqrt(Reff);
+		real sqrt_meff_K_hertz = sqrt((mi*mj)/(mi+mj) * K_hertz);
+
+		real en = prop.en_;
+		if (addDissipationModel_==2)
+		{
+			en = sqrt(1+((pow(prop.en_,2)-1)*f));
+		}
+		else if (addDissipationModel_==3)
+		{
+			
+			en = exp((pow(f,static_cast<real>(1.5))*log(prop.en_)*sqrt( (1-((pow(log(prop.en_),2))/(pow(log(prop.en_),2)+pow(Pi,2))))/(1-(pow(f,3)*(pow(log(prop.en_),2))/(pow(log(prop.en_),2)+pow(Pi,2)))) ) ));
+		}
+		
+		real Kn = static_cast<real>(4.0/3.0) * prop.Yeff_ * sqrt(Reff*ovrlp_n);
+
+		real ethan_  = -2.0*log(en)*sqrt(Kn)/
+			sqrt(pow(log(en),2)+ pow(Pi,2));
+
+		FCn = ( - Kn*ovrlp_n - 
+			   sqrt_meff_K_hertz*ethan_*pow(ovrlp_n,static_cast<real>(0.25))*vrn)*Nij;
+
+		FCt = (f*f)*(- static_cast<real>(8.0) * prop.Geff_ * sqrt(Reff*ovrlp_n) ) * history.overlap_t_;
+
+		real ft = length(FCt);
+		real ft_fric = prop.mu_ * length(FCn);
+
+		// apply friction 
+		if(ft > ft_fric)
+		{
+			if( length(history.overlap_t_) >0.0)
+			{
+				if constexpr (limited)
+				{
+					real kt = static_cast<real>(8.0) * prop.Geff_ * sqrt(Reff*ovrlp_n);
+					FCt *= (ft_fric/ft);
+                	history.overlap_t_ = - FCt/(f*f*kt);
+				}
+				else
+				{
+					FCt = (FCt/ft)*ft_fric;	
+				}
+				
+			}
+			else
+			{
+				FCt = 0.0;
+			}
+		}
+
+	}
+   
+
+};
+
+} //pFlow::CFModels
+
+#endif

src/Interaction/Models/contactForce/cGRelativeLinearCF.hpp

@@ -26,11 +26,6 @@ Licence:
 
 
 
-
-
-
-
-
 namespace pFlow::cfModels
 {
 
@@ -49,16 +44,15 @@ public:
 	{		
 		real 	kn_ 	= 1000.0;
 		real 	kt_	= 800.0;
-		real 	en_ 	= 0.0;
+		real 	en_ 	= 1.0;
-		real 	ethat_ 	= 0.0;
 		real  mu_ 	= 0.00001;
 
 		INLINE_FUNCTION_HD
 		linearProperties(){}
 
 		INLINE_FUNCTION_HD
-		linearProperties(real kn, real kt, real en, real etha_t, real mu  ):
+		linearProperties(real kn, real kt, real en, real mu  ):
-			kn_(kn), kt_(kt), en_(en),ethat_(etha_t), mu_(mu)
+			kn_(kn), kt_(kt), en_(en), mu_(mu)
 		{}		
 
 		INLINE_FUNCTION_HD
@@ -93,7 +87,6 @@ protected:
 		auto kn = dict.getVal<realVector>("kn");
 		auto kt = dict.getVal<realVector>("kt");
 		auto en = dict.getVal<realVector>("en");
-		auto et = dict.getVal<realVector>("et");
 		auto mu = dict.getVal<realVector>("mu");
 		
 
@@ -114,12 +107,6 @@ protected:
 			return false;
 		}
 
-		if(nElem != et.size())
-		{
-			fatalErrorInFunction<<
-			"sizes of kn("<<nElem<<") and et("<<et.size()<<") do not match.\n";
-			return false;
-		}
 
 		if(nElem != mu.size())
 		{
@@ -148,30 +135,16 @@ protected:
 		}
 
 
-		realVector etha_t("etha_t", nElem);
-
-
-		ForAll(i , kn)
-		{
-
-
-			etha_t[i] = -2.0*log( et[i])*sqrt(kt[i]*2/7) /
-					sqrt(log(pow(et[i],2.0))+ pow(Pi,2.0));
-					
-					
-		}
-
 		Vector<linearProperties> prop("prop", nElem);
 		ForAll(i,kn)
 		{
-			prop[i] = {kn[i], kt[i], en[i], etha_t[i], mu[i]  };
+			prop[i] = {kn[i], kt[i], en[i], mu[i]  };
 		}
 
 		linearProperties_.assign(prop);
 		
+
 		auto adm = dict.getVal<word>("additionalDissipationModel");
-
-		
 		if(adm == "none")
 		{
 			addDissipationModel_ = 1;
@@ -282,27 +255,40 @@ public:
 
 		history.overlap_t_ += Vt*dt;
 
-		real mi = 3*Pi/4*pow(Ri,3.0)*rho_[propId_i];
+		real mi = 3*Pi/4*pow(Ri,3)*rho_[propId_i];
-		real mj = 3*Pi/4*pow(Rj,3.0)*rho_[propId_j];
+		real mj = 3*Pi/4*pow(Rj,3)*rho_[propId_j];
 
 		real sqrt_meff = sqrt((mi*mj)/(mi+mj));
 		
-		
+		real en = prop.en_;
 		if (addDissipationModel_==2)
 		{
-			prop.en_ = sqrt(1+((pow(prop.en_,2)-1)*f_));
+			en = sqrt(1+((pow(prop.en_,2)-1)*f_));
 		}
 		else if (addDissipationModel_==3)
 		{
-			auto pie =3.14;
+			
-			prop.en_ = exp((pow(f_,1.5)*log(prop.en_)*sqrt( (1-((pow(log(prop.en_),2))/(pow(log(prop.en_),2)+pow(pie,2))))/(1-(pow(f_,3)*(pow(log(prop.en_),2))/(pow(log(prop.en_),2)+pow(pie,2)))) ) ));
+			en = exp(
+				(
+					pow(f_,static_cast<real>(1.5))*
+					log(prop.en_)*
+					sqrt( 
+						(1-((pow(log(prop.en_),static_cast<real>(2.0))
+						)/
+						(
+							pow(log(prop.en_),static_cast<real>(2.0))+
+							pow(Pi,static_cast<real>(2.0)))))/
+							(1-(pow(f_,3)*(pow(log(prop.en_),2))/
+							(pow(log(prop.en_),static_cast<real>(2.0))+
+							pow(Pi,static_cast<real>(2.0))))) ) ));
 		}
-		real ethan_  = -2.0*log(prop.en_)*sqrt(prop.kn_)/
+
-			sqrt(pow(log(prop.en_),2.0)+ pow(Pi,2.0));
+		real ethan_  = -2.0*log(en)*sqrt(prop.kn_)/
+			sqrt(pow(log(en),2)+ pow(Pi,2));
 		
 
-		FCn = (  -pow(f_,1.0)*prop.kn_ * ovrlp_n - sqrt_meff * pow(f_,0.5) * ethan_ * vrn)*Nij;
+		FCn = ( -f_*prop.kn_ * ovrlp_n - sqrt_meff * pow(f_,half) * ethan_ * vrn)*Nij;
-		FCt = ( -pow(f_,1.0)*prop.kt_ * history.overlap_t_ - sqrt_meff * pow(f_,0.5) *  prop.ethat_*Vt);
+		FCt = ( -f_*prop.kt_ * history.overlap_t_);
 		
 
 		

src/Interaction/Models/contactForce/linearCF.hpp

@@ -139,10 +139,10 @@ protected:
 		ForAll(i , kn)
 		{
 			etha_n[i] = -2.0*log(en[i])*sqrt(kn[i])/
-					sqrt(pow(log(en[i]),2.0)+ pow(Pi,2.0));
+					sqrt(pow(log(en[i]),static_cast<real>(2.0))+ pow(Pi,static_cast<real>(2.0)));
 
 			etha_t[i] = -2.0*log( et[i]*sqrt(kt[i]) )/
-					sqrt(pow(log(et[i]),2.0)+ pow(Pi,2.0));
+					sqrt(pow(log(et[i]),static_cast<real>(2.0))+ pow(Pi,static_cast<real>(2.0)));
 		}
 
 		Vector<linearProperties> prop("prop", nElem);
@@ -243,8 +243,8 @@ public:
 
 		history.overlap_t_ += Vt*dt;
 
-		real mi = 3*Pi/4*pow(Ri,3.0)*rho_[propId_i];
+		real mi = 3*Pi/4*pow(Ri,static_cast<real>(3.0))*rho_[propId_i];
-		real mj = 3*Pi/4*pow(Rj,3.0)*rho_[propId_j];
+		real mj = 3*Pi/4*pow(Rj,static_cast<real>(3.0))*rho_[propId_j];
 
 		real sqrt_meff = sqrt((mi*mj)/(mi+mj));
 

src/Interaction/Models/contactForce/nonLinearCF.hpp

@@ -131,7 +131,7 @@ protected:
 			// we take out sqrt(meff*K_hertz) here and then consider this term 
 			// when calculating damping part. 
 			etha_n[i] = -2.2664*log(en[i])/
-					sqrt(pow(log(en[i]),2.0)+ pow(Pi,2.0));
+					sqrt(pow(log(en[i]),static_cast<real>(2.0))+ pow(Pi,static_cast<real>(2.0)));
 
 			// no damping for tangential part 
 
@@ -255,7 +255,7 @@ public:
 		// apply friction 
 		if(ft > ft_fric)
 		{
-			if( length(history.overlap_t_) >0.0)
+			if( length(history.overlap_t_) >zero)
 			{
 				if constexpr (limited)
 				{

src/Interaction/Models/grainContactForceModels.hpp

@@ -23,6 +23,7 @@ Licence:
 
 #include "cGAbsoluteLinearCF.hpp"
 #include "cGRelativeLinearCF.hpp"
+#include "cGNonLinearCF.hpp"
 
 #include "grainRolling.hpp"
 
@@ -38,6 +39,10 @@ using nonLimitedCGAbsoluteLinearGrainRolling = grainRolling<cGAbsoluteLinear<fal
 using limitedCGRelativeLinearGrainRolling   = grainRolling<cGRelativeLinear<true>>;
 using nonLimitedCGRelativeLinearGrainRolling = grainRolling<cGRelativeLinear<false>>;
 
+using limitedCGNonLinearGrainRolling   = grainRolling<cGNonLinear<true>>;
+using nonLimitedCGNonLinearGrainRolling = grainRolling<cGNonLinear<false>>;
+
+
 }
 
 

src/Interaction/Models/rolling/grainRolling.hpp

@@ -96,10 +96,10 @@ public:
 		realx3 w_hat = wi-wj;
 		real w_hat_mag = length(w_hat);
 
-		if( !equal(w_hat_mag,0.0) )
+		if( !equal(w_hat_mag,static_cast<real>(0.0)) )
 			w_hat /= w_hat_mag;
 		else
-			w_hat = 0.0;
+			w_hat = static_cast<real>(0.0);
 
 		auto Reff = (Ri*Rj)/(Ri+Rj);
 

src/Interaction/Models/rolling/normalRolling.hpp

@@ -94,10 +94,10 @@ public:
 		realx3 w_hat = wi-wj;
 		real w_hat_mag = length(w_hat);
 
-		if( !equal(w_hat_mag,0.0) )
+		if( !equal(w_hat_mag,static_cast<real>(0.0)) )
 			w_hat /= w_hat_mag;
 		else
-			w_hat = 0.0;
+			w_hat = static_cast<real>(0.0);
 
 		auto Reff = (Ri*Rj)/(Ri+Rj);
 

src/Interaction/contactLists/sortedContactList.hpp

@@ -80,13 +80,17 @@ public:
 
 	TypeInfoNV("sortedContactList");
 
-
+	sortedContactList(uint32 initialSize =1)
-	explicit sortedContactList(uint32 initialSize =1)
 	:
-		SortedPairs(initialSize),
+		sortedContactList("sortedContactList", initialSize)
-		values_("values", SortedPairs::capacity()),
+	{}
-		sortedPairs0_("sortedPairs0", SortedPairs::capacity()),
+	
-		values0_("values0", SortedPairs::capacity())
+	sortedContactList(const word& name, uint32 initialSize =1)
+	:
+		SortedPairs(name, initialSize),
+		values_(groupNames(name, "values"), SortedPairs::capacity()),
+		sortedPairs0_(groupNames(name, "sortedPairs0"), SortedPairs::capacity()),
+		values0_(groupNames(name, "values0"), SortedPairs::capacity())
 	{}
 
 	bool beforeBroadSearch()

src/Interaction/contactLists/sortedPairs.hpp

@@ -110,11 +110,11 @@ public:
 
 
 	// constructors
-	explicit sortedPairs(uint32 initialSize =1)
+	explicit sortedPairs(const word& name, uint32 initialSize =1)
 	:
 		UnsortedPairs(initialSize),
-		flags_("flags_",UnsortedPairs::capacity()+1),
+		flags_( groupNames(name, "flags_"), UnsortedPairs::capacity()+1),
-		sortedPairs_("sortedPairs_",UnsortedPairs::capacity())
+		sortedPairs_(groupNames(name, "sortedPairs_"), UnsortedPairs::capacity())
 	{}
 
 	
@@ -193,7 +193,7 @@ public:
 
 		if( capacity+1 > flags_.size() )
 		{
-			reallocNoInit(flags_, capacity+1);
+			reallocInit(flags_, capacity+1);
 		}
 
 		// fill the flags
@@ -219,7 +219,7 @@ public:
 		{
 			// get more space to prevent reallocations in next iterations
 			uint32 len = size_*1.1+1; 
-			reallocNoInit(sortedPairs_, len);
+			reallocInit(sortedPairs_, len);
 		}
 		
 		Kokkos::parallel_for(
@@ -232,6 +232,7 @@ public:
 		sort(sortedPairs_, 0, size_ );
 		
 		
+
 	}
 
 	INLINE_FUNCTION_HD

src/Interaction/contactLists/unsortedContactList.hpp

@@ -82,11 +82,16 @@ public:
 	TypeInfoNV("unsortedContactList");
 
 	explicit unsortedContactList(uint32 capacity=1)
+	:
+	unsortedContactList("unsortedContactList", capacity)
+	{}
+	
+	unsortedContactList(const word& name, uint32 capacity=1)
 	:
 		UnsortedPairs(capacity),
-		values_("values", UnsortedPairs::capacity()),
+		values_(groupNames(name, "values"), UnsortedPairs::capacity()),
 		container0_(capacity),
-		values0_("values0",container0_.capacity())
+		values0_(groupNames(name, "values0"),container0_.capacity())
 	{}