xfygogo/phasicFlow
1
0
Fork 0
mirror of https://github.com/PhasicFlow/phasicFlow.git synced 2025-07-28 03:27:05 +00:00
Code Issues Packages Projects Releases Wiki Activity

Merge pull request #92 from PhasicFlow/develop

Browse Source 
Merge (first patch of version 1.x)
This commit is contained in:
Hamidreza Norouzi
2024-03-25 01:03:20 +03:30
committed by GitHub
parent 9facc73055 c7e5b30959
commit af75d6a09c
374 changed files with 24536 additions and 13176 deletions
Download Patch File Download Diff File Expand all files Collapse all files

63
CMakeLists.txt
View File

@ -1,74 +1,73 @@
cmake_minimum_required(VERSION 3.16 FATAL_ERROR)
# set the project name and version
project(phasicFlow VERSION 0.1 )
project(phasicFlow VERSION 1.0 )
set(CMAKE_CXX_STANDARD 17 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 Release CACHE STRING "build type" FORCE)
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})
mark_as_advanced(FORCE var Kokkos_ENABLE_CUDA_LAMBDA)
mark_as_advanced(FORCE var Kokkos_ENABLE_OPENMP)
mark_as_advanced(FORCE var Kokkos_ENABLE_SERIAL)
mark_as_advanced(FORCE var Kokkos_ENABLE_CUDA_LAMBDA)
mark_as_advanced(FORCE var BUILD_SHARED_LIBS)
include(cmake/globals.cmake)
option(USE_STD_PARALLEL_ALG "Use TTB std parallel algorithms" ON)
#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} ${phasicFlow_BINARY_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)
option(pFlow_Build_Cuda "Build phasicFlow and backends for Cuda execution" OFF)
option(pFlow_Build_Double "Build phasicFlow with double precision variables" ON)
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)
set(BUILD_SHARED_LIBS ON CACHE BOOL "Build using shared libraries" FORCE)
if(pFlow_Build_Serial)
set(Kokkos_ENABLE_SERIAL ON CACHE BOOL "Serial execution" FORCE)
set(Kokkos_ENABLE_OPENMP OFF 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_ENABLE_CUDA_CONSTEXPR OFF CACHE BOOL "Enable constexpr on cuda code" FORCE)
elseif(pFlow_Build_OpenMP )
set(Kokkos_ENABLE_SERIAL ON CACHE BOOL "Serial execution" FORCE)
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_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)
set(Kokkos_ENABLE_OPENMP OFF CACHE BOOL "OpenMP execution" FORCE)
set(Kokkos_ENABLE_CUDA ON CACHE BOOL "Cuda execution" FORCE)
set(Kokkos_ENABLE_CUDA_LAMBDA ON CACHE BOOL "Cuda execution" FORCE)
set(Kokkos_ENABLE_CUDA_CONSTEXPR ON CACHE BOOL "Enable constexpr on cuda code" FORCE)
endif()
include(cmake/globals.cmake)
message(STATUS "Valid file extensions are ${validFiles}")
if(pFlow_Build_MPI)
find_package(MPI REQUIRED)
endif()
include(cmake/makeLibraryGlobals.cmake)
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}")
#main subdirectories of the code
set(Kokkos_Source_DIR)
if(DEFINED ENV{Kokkos_DIR})
set(Kokkos_Source_DIR $ENV{Kokkos_DIR})
# add_subdirectory($ENV{Kokkos_DIR} ${phasicFlow_BINARY_DIR}/kokkos)
# message(STATUS "Kokkos directory is $ENV{Kokkos_DIR}")
else()
# add_subdirectory($ENV{HOME}/Kokkos/kokkos ${phasicFlow_BINARY_DIR}/kokkos)
set(Kokkos_Source_DIR $ENV{HOME}/Kokkos/kokkos)
endif()
message(STATUS "Kokkos source directory is ${Kokkos_Source_DIR}")
add_subdirectory(${Kokkos_Source_DIR} ${phasicFlow_BINARY_DIR}/kokkos)
add_subdirectory(src)
@ -76,13 +75,11 @@ add_subdirectory(solvers)
add_subdirectory(utilities)
add_subdirectory(DEMSystems)
#add_subdirectory(testIO)
#add_subdirectory(DEMSystems)
add_subdirectory(testIO)
install(FILES "${PROJECT_BINARY_DIR}/phasicFlowConfig.H"
DESTINATION include
)
DESTINATION include)
include(InstallRequiredSystemLibraries)
set(CPACK_RESOURCE_FILE_LICENSE "${CMAKE_CURRENT_SOURCE_DIR}/LICENSE")

6
cmake/bashrc
View File

@ -1,5 +1,5 @@
export pFlow_PROJECT_VERSION=v0.1
export pFlow_PROJECT_VERSION=v-1.0
export pFlow_PROJECT=phasicFlow
@ -20,6 +20,8 @@ export pFlow_SRC_DIR="$pFlow_PROJECT_DIR/src"
export Kokkos_DIR="$kokkosDir"
export Zoltan_DIR="$projectDir/Zoltan"
# Cleanup variables (done as final statement for a clean exit code)
unset projectDir
@ -27,5 +29,7 @@ export PATH="$pFlow_BIN_DIR:$PATH"
export LD_LIBRARY_PATH="$pFlow_LIB_DIR:$LD_LIBRARY_PATH"
export LD_LIBRARY_PATH="$Zoltan_DIR/lib:$LD_LIBRARY_PATH"
#------------------------------------------------------------------------------

49
cmake/globals.cmake
View File

@ -2,3 +2,52 @@
set(validFiles)
list(APPEND validFiles *.C *.cpp *.cxx *.c *.cu *.H *.hpp *.hxx *.h *.cuh)
macro(Kokkos_cmake_settings)
#mark_as_advanced(FORCE var Kokkos_ENABLE_CUDA_LAMBDA)
mark_as_advanced(FORCE var Kokkos_CXX_STANDARD)
#mark_as_advanced(FORCE var Kokkos_ENABLE_CUDA_CONSTEXPR)
mark_as_advanced(FORCE var Kokkos_ENABLE_OPENMP)
mark_as_advanced(FORCE var Kokkos_ENABLE_SERIAL)
mark_as_advanced(FORCE var Kokkos_ENABLE_CUDA)
mark_as_advanced(FORCE var Kokkos_ENABLE_HIP)
mark_as_advanced(FORCE var Kokkos_ENABLE_AGGRESSIVE_VECTORIZATION)
mark_as_advanced(FORCE var Kokkos_ENABLE_BENCHMARKS)
mark_as_advanced(FORCE var Kokkos_ENABLE_COMPILE_AS_CMAKE_LANGUAGE)
mark_as_advanced(FORCE var Kokkos_ENABLE_CUDA_LDG_INTRINSIC)
mark_as_advanced(FORCE var Kokkos_ENABLE_CUDA_RELOCATABLE_DEVICE_CODE)
mark_as_advanced(FORCE var Kokkos_ENABLE_CUDA_UVM)
mark_as_advanced(FORCE var Kokkos_ENABLE_DEPRECATED_CODE_3)
mark_as_advanced(FORCE var Kokkos_ENABLE_DEPRECATED_CODE_4)
mark_as_advanced(FORCE var Kokkos_ENABLE_DEPRECATION_WARNINGS)
mark_as_advanced(FORCE var Kokkos_ENABLE_DESUL_ATOMICS_EXTERNAL)
mark_as_advanced(FORCE var Kokkos_ENABLE_EXAMPLES)
mark_as_advanced(FORCE var Kokkos_ENABLE_HEADER_SELF_CONTAINMENT_TESTS)
mark_as_advanced(FORCE var Kokkos_ENABLE_HIP_MULTIPLE_KERNEL_INSTANTIATIONS)
mark_as_advanced(FORCE var Kokkos_ENABLE_HIP_RELOCATABLE_DEVICE_CODE)
mark_as_advanced(FORCE var Kokkos_ENABLE_HPX)
mark_as_advanced(FORCE var Kokkos_ENABLE_HWLOC)
mark_as_advanced(FORCE var Kokkos_ENABLE_IMPL_CUDA_MALLOC_ASYNC)
mark_as_advanced(FORCE var Kokkos_ENABLE_IMPL_HPX_ASYNC_DISPATCH)
mark_as_advanced(FORCE var Kokkos_ENABLE_LARGE_MEM_TESTS)
mark_as_advanced(FORCE var Kokkos_ENABLE_DEBUG_DUALVIEW_MODIFY_CHECK)
mark_as_advanced(FORCE var Kokkos_ENABLE_LIBQUADMATH)
mark_as_advanced(FORCE var Kokkos_ENABLE_LIBRT)
mark_as_advanced(FORCE var Kokkos_ENABLE_MEMKIND)
mark_as_advanced(FORCE var Kokkos_ENABLE_ONEDPL)
mark_as_advanced(FORCE var Kokkos_ENABLE_OPENACC)
mark_as_advanced(FORCE var Kokkos_ENABLE_OPENMPTARGET)
mark_as_advanced(FORCE var Kokkos_ENABLE_ROCM)
mark_as_advanced(FORCE var Kokkos_ENABLE_SYCL)
mark_as_advanced(FORCE var Kokkos_ENABLE_TESTS)
mark_as_advanced(FORCE var Kokkos_ENABLE_THREADS)
mark_as_advanced(FORCE var Kokkos_ENABLE_TUNING)
mark_as_advanced(FORCE var Kokkos_ENABLE_UNSUPPORTED_ARCHS)
mark_as_advanced(FORCE var Kokkos_HPX_DIR)
mark_as_advanced(FORCE var Kokkos_HWLOC_DIR)
mark_as_advanced(FORCE var Kokkos_MEMKIND_DIR)
mark_as_advanced(FORCE var Kokkos_ROCM_DIR)
mark_as_advanced(FORCE var Kokkos_THREADS_DIR)
endmacro()

2
cmake/makeExecutableGlobals.cmake
View File

@ -31,7 +31,7 @@ target_include_directories(${target_name}
message(STATUS "\nCreating make file for executable ${target_name}")
message(STATUS " ${target_name} link libraries are: ${${target_link_libs}}")
message(STATUS " ${target_name} source files are: ${source_files}")
message(STATUS " ${target_name} source files are: ${${source_files}}")
message(STATUS " ${target_name} include dirs are: ${includeDirs}\n")

2
cmake/makeLibraryGlobals.cmake
View File

@ -42,7 +42,7 @@ target_include_directories(${target_name}
message(STATUS "\nCreating make file for library ${target_name}")
message(STATUS " ${target_name} link libraries are: ${${target_link_libs}}")
message(STATUS " ${target_name} source files are: ${source_files}")
message(STATUS " ${target_name} include dirs are: ${includeDirs}\n")
message(STATUS " ${target_name} include dirs are: ${includeDirs}\n \n")
install(TARGETS ${target_name} DESTINATION lib)
install(FILES ${includeFiles} DESTINATION include/${target_name})

5
phasicFlowConfig.H.in
View File

@ -5,5 +5,6 @@
#cmakedefine pFlow_Build_Serial
#cmakedefine pFlow_Build_OpenMP
#cmakedefine pFlow_Build_Cuda
#cmakedefine USE_STD_PARALLEL_ALG
#cmakedefine pFlow_Build_Double
#cmakedefine pFlow_STD_Parallel_Alg
#cmakedefine pFlow_Build_Double
#cmakedefine pFlow_Build_MPI

2
solvers/CMakeLists.txt
View File

@ -1,6 +1,6 @@
#add_subdirectory(iterateSphereParticles)
add_subdirectory(iterateSphereParticles)
add_subdirectory(iterateGeometry)

16
solvers/iterateGeometry/iterateGeometry.cpp
View File

@ -28,13 +28,14 @@ Licence:
* (https://github.com/PhasicFlow/phasicFlow/tree/main/tutorials/iterateGeometry)
* folder.
*/
#include "systemControl.hpp"
#include "geometryMotion.hpp"
#include "commandLine.hpp"
#include "readControlDict.hpp"
using pFlow::commandLine;
#include "vocabs.hpp"
#include "systemControl.hpp"
#include "geometry.hpp"
#include "commandLine.hpp"
//#include "readControlDict.hpp"
using namespace pFlow;
int main( int argc, char* argv[] )
{
@ -54,6 +55,8 @@ commandLine cmds(
// this should be palced in each main
processors::initProcessors(argc, argv);
initialize_pFlowProcessors();
#include "initialize_Control.hpp"
#include "setProperty.hpp"
@ -68,6 +71,7 @@ commandLine cmds(
// this should be palced in each main
#include "finalize.hpp"
processors::finalizeProcessors();
}

7
solvers/iterateSphereParticles/CMakeLists.txt Normal file
View File

@ -0,0 +1,7 @@
set(source_files
iterateSphereParticles.cpp
)
set(link_lib Particles)
pFlow_make_executable_install(iterateSphereParticles source_files link_lib)

25
solvers/iterateSphereParticles/createDEMComponents.hpp Executable file
View File

@ -0,0 +1,25 @@
/*------------------------------- 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.
-----------------------------------------------------------------------------*/
//
REPORT(0)<<"\nReading sphere particles . . ."<<END_REPORT;
sphereParticles sphParticles(Control, proprties);
WARNING<<"Particle insertion has not been set yet!"<<END_WARNING;

89
solvers/iterateSphereParticles/iterateSphereParticles.cpp Executable file
View File

@ -0,0 +1,89 @@
/*------------------------------- 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.
-----------------------------------------------------------------------------*/
/**
* \file sphereGranFlow.cpp
* \brief sphereGranFlow solver
*
* This solver simulate granular flow of cohesion-less, spherical particles.
* Particle insertion can be activated in the solver to gradually insert
* particles into the simulation. Geometry can be defined generally with
* built-in features of the code or through ASCII stl files or a combination
* of both. For more information refer to [tutorials/sphereGranFlow/]
* (https://github.com/PhasicFlow/phasicFlow/tree/main/tutorials/sphereGranFlow) folder.
*/
#include "vocabs.hpp"
#include "property.hpp"
#include "sphereParticles.hpp"
#include "systemControl.hpp"
#include "commandLine.hpp"
using namespace pFlow;
/**
* DEM solver for simulating granular flow of cohesion-less particles.
*
* In the root case directory just simply enter the following command to
* run the simulation. For command line options use flag -h.
*/
int main( int argc, char* argv[])
{
commandLine cmds(
"sphereGranFlow",
"DEM solver for non-cohesive spherical particles with particle insertion "
"mechanism and moving geometry");
bool isCoupling = false;
if(!cmds.parse(argc, argv)) return 0;
// this should be palced in each main
processors::initProcessors(argc, argv);
initialize_pFlowProcessors();
#include "initialize_Control.hpp"
#include "setProperty.hpp"
#include "createDEMComponents.hpp"
REPORT(0)<<"\nStart of time loop . . .\n"<<END_REPORT;
do
{
sphParticles.beforeIteration();
sphParticles.iterate();
sphParticles.afterIteration();
}while(Control++);
REPORT(0)<<"\nEnd of time loop.\n"<<END_REPORT;
// this should be palced in each main
#include "finalize.hpp"
processors::finalizeProcessors();
}

10
solvers/sphereGranFlow/createDEMComponents.hpp
View File

@ -19,11 +19,11 @@ Licence:
-----------------------------------------------------------------------------*/
//
REPORT(0)<<"\nReading sphere particles . . ."<<endREPORT;
REPORT(0)<<"\nReading sphere particles . . ."<<END_REPORT;
sphereParticles sphParticles(Control, proprties);
//
REPORT(0)<<"\nCreating particle insertion object . . ."<<endREPORT;
REPORT(0)<<"\nCreating particle insertion object . . ."<<END_REPORT;
/*auto& sphInsertion =
Control.caseSetup().emplaceObject<sphereInsertion>(
objectFile(
@ -36,12 +36,12 @@ REPORT(0)<<"\nCreating particle insertion object . . ."<<endREPORT;
sphParticles.shapes()
);*/
auto sphInsertion = sphereInsertion(
/*auto sphInsertion = sphereInsertion(
Control.caseSetup().path()+insertionFile__,
sphParticles,
sphParticles.shapes());
sphParticles.shapes());*/
REPORT(0)<<"\nCreating interaction model for sphere-sphere contact and sphere-wall contact . . ."<<endREPORT;
REPORT(0)<<"\nCreating interaction model for sphere-sphere contact and sphere-wall contact . . ."<<END_REPORT;
auto interactionPtr = interaction::create(
Control,
sphParticles,

60
solvers/sphereGranFlow/sphereGranFlow.cpp
View File

@ -30,30 +30,22 @@ Licence:
* (https://github.com/PhasicFlow/phasicFlow/tree/main/tutorials/sphereGranFlow) folder.
*/
#include "vocabs.hpp"
#include "phasicFlowKokkos.hpp"
#include "systemControl.hpp"
#include "commandLine.hpp"
#include "property.hpp"
#include "geometry.hpp"
#include "sphereParticles.hpp"
#include "Insertions.hpp"
#include "systemControl.hpp"
#include "contactSearch.hpp"
#include "sphereInteraction.hpp"
#include "commandLine.hpp"
#include "readControlDict.hpp"
using pFlow::output;
using pFlow::endl;
using pFlow::property;
using pFlow::sphereParticles;
using pFlow::objectFile;
using pFlow::sphereInsertion;
using pFlow::insertionFile__;
using pFlow::interactionFile__;
using pFlow::contactSearch;
using pFlow::interaction;
using pFlow::commandLine;
#include "interaction.hpp"
//#include "Insertions.hpp"
//#include "readControlDict.hpp"
using namespace pFlow;
/**
* DEM solver for simulating granular flow of cohesion-less particles.
*
@ -73,52 +65,58 @@ bool isCoupling = false;
if(!cmds.parse(argc, argv)) return 0;
// this should be palced in each main
processors::initProcessors(argc, argv);
initialize_pFlowProcessors();
#include "initialize_Control.hpp"
#include "setProperty.hpp"
#include "setSurfaceGeometry.hpp"
#include "createDEMComponents.hpp"
REPORT(0)<<"\nStart of time loop . . .\n"<<endREPORT;
REPORT(0)<<"\nStart of time loop . . .\n"<<END_REPORT;
do
{
if(! sphInsertion.insertParticles(
/*if(! sphInsertion.insertParticles(
Control.time().currentTime(),
Control.time().dt() ) )
{
fatalError<<
"particle insertion failed in sphereDFlow solver.\n";
return 1;
}
}*/
// set force to zero
surfGeometry.beforeIteration();
sphInteraction.beforeIteration();
// set force to zero, predict, particle deletion and etc.
sphParticles.beforeIteration();
sphInteraction.beforeIteration();
sphInteraction.iterate();
sphInteraction.iterate();
surfGeometry.iterate();
sphParticles.iterate();
surfGeometry.iterate();
sphInteraction.afterIteration();
surfGeometry.afterIteration();
sphParticles.afterIteration();
surfGeometry.afterIteration();
}while(Control++);
REPORT(0)<<"\nEnd of time loop.\n"<<endREPORT;
REPORT(0)<<"\nEnd of time loop.\n"<<END_REPORT;
// this should be palced in each main
#include "finalize.hpp"
processors::finalizeProcessors();
}

4
src/CMakeLists.txt
View File

@ -7,9 +7,9 @@ add_subdirectory(Property)
add_subdirectory(Particles)
add_subdirectory(Geometry)
add_subdirectory(Interaction)
add_subdirectory(MotionModel)
add_subdirectory(Geometry)

533
src/Geometry/geometry/geometry.cpp
View File

@ -19,33 +19,41 @@ Licence:
-----------------------------------------------------------------------------*/
#include "geometry.hpp"
#include "systemControl.hpp"
#include "vocabs.hpp"
bool pFlow::geometry::findPropertyId()
bool pFlow::geometry::createPropertyId()
{
int8Vector propId(0, surface().capacity(),RESERVE());
propId.clear();
uint32 pId;
ForAll(matI, materialName_)
if(materialName_.size() != numSurfaces() )
{
fatalErrorInFunction<<
"number of subSurface and material names do not match"<<endl;
return false;
}
uint32Vector propId(
"propId",
capacity(),
size(),
RESERVE());
ForAll(i, materialName_)
{
uint32 pIdx =0;
if( !wallProperty_.nameToIndex( materialName_[matI], pId ) )
if( !wallProperty_.nameToIndex(materialName_[i], pIdx) )
{
fatalErrorInFunction<<
"material name for the geometry is invalid: "<< materialName_[matI]<<endl;
return false;
fatalErrorInFunction<<
"Property/material name is invalid "<<materialName_[i]<<
". A list of valid names are \n"<< wallProperty_.materials()<<endl;
return false;
}
int32 surfSize = surface().surfNumTriangles(matI);
for(int32 i=0; i<surfSize; i++)
{
propId.push_back(pId);
}
}
auto triRange = subSurfaceRange(i);
propId.fill(triRange.start(), triRange.end(), pIdx);
}
propertyId_.assign(propId);
@ -53,83 +61,361 @@ bool pFlow::geometry::findPropertyId()
}
void pFlow::geometry::zeroForce()
{
contactForceWall_.fill(zero3);
}
pFlow::geometry::geometry
(
systemControl& control,
const property& prop
)
:
demGeometry(control),
multiTriSurface
(
objectFile
(
triSurfaceFile__,
"",
objectFile::READ_ALWAYS,
objectFile::WRITE_ALWAYS
),
&control.geometry()
),
demComponent
(
"geometry",
control
),
wallProperty_(prop),
geometryRepository_(control.geometry()),
triSurface_(
control.geometry().emplaceObject<multiTriSurface>(
objectFile(
"triSurface",
"",
objectFile::READ_ALWAYS,
objectFile::WRITE_ALWAYS
)
)
propertyId_
(
objectFile
(
"propertyId",
"",
objectFile::READ_NEVER,
objectFile::WRITE_NEVER
),
motionComponentName_(
control.geometry().emplaceObject<wordField>(
objectFile(
"motionComponentName",
"",
objectFile::READ_ALWAYS,
objectFile::WRITE_ALWAYS
),
"motionNamesList"
)
*this,
0u
),
contactForceWall_
(
objectFile
(
"contactForcWall",
"",
objectFile::READ_IF_PRESENT,
objectFile::WRITE_NEVER
),
materialName_(
control.geometry().emplaceObject<wordField>(
objectFile(
"materialName",
"",
objectFile::READ_ALWAYS,
objectFile::WRITE_ALWAYS
),
"materialNamesList"
)
*this,
zero3
),
normalStressWall_
(
objectFile
(
"normalStressWall",
"",
objectFile::READ_IF_PRESENT,
objectFile::WRITE_NEVER
),
propertyId_(
control.geometry().emplaceObject<int8TriSurfaceField_D>(
objectFile(
"propertyId",
"",
objectFile::READ_NEVER,
objectFile::WRITE_NEVER),
surface(),
0 ) ),
contactForceWall_(
control.geometry().emplaceObject<realx3TriSurfaceField_D>(
objectFile(
"contactForceWall",
"",
objectFile::READ_IF_PRESENT,
objectFile::WRITE_ALWAYS),
surface(),
zero3) ),
stressWall_(
control.geometry().emplaceObject<realx3TriSurfaceField_D>(
objectFile(
"stressWall",
"",
objectFile::READ_IF_PRESENT,
objectFile::WRITE_ALWAYS),
surface(),
zero3) )
*this,
zero3
),
shearStressWall_
(
objectFile
(
"shearStressWall",
"",
objectFile::READ_IF_PRESENT,
objectFile::WRITE_NEVER
),
*this,
zero3
)
{
if(!findPropertyId())
readWholeObject_ = false;
if( !IOobject::readObject() )
{
fatalErrorInFunction<<
"Error in reading from file "<<IOobject::path()<<endl;
fatalExit;
}
readWholeObject_ = true;
if( this->numSurfaces() != motionComponentName_.size() )
{
fatalErrorInFunction<<
"Number of surfaces is not equal to number of motion component names"<<endl;
fatalExit;
}
if(!createPropertyId())
{
fatalExit;
}
}
pFlow::geometry::geometry
(
systemControl &control,
const property &prop,
multiTriSurface &surf,
const wordVector &motionCompName,
const wordVector &materialName,
const dictionary& motionDict
)
:
multiTriSurface
(
objectFile
(
triSurfaceFile__,
"",
objectFile::READ_NEVER,
objectFile::WRITE_ALWAYS
),
&control.geometry(),
surf
),
demComponent
(
"geometry",
control
),
wallProperty_
(
prop
),
propertyId_
(
objectFile
(
"propertyId",
"",
objectFile::READ_NEVER,
objectFile::WRITE_NEVER
),
*this,
0u
),
contactForceWall_
(
objectFile
(
"contactForcWall",
"",
objectFile::READ_IF_PRESENT,
objectFile::WRITE_NEVER
),
*this,
zero3
),
normalStressWall_
(
objectFile
(
"normalStressWall",
"",
objectFile::READ_IF_PRESENT,
objectFile::WRITE_NEVER
),
*this,
zero3
),
shearStressWall_
(
objectFile
(
"shearStressWall",
"",
objectFile::READ_IF_PRESENT,
objectFile::WRITE_NEVER
),
*this,
zero3
)
{
motionComponentName_.assign(motionCompName);
materialName_.assign(materialName);
if( this->numSurfaces() != motionComponentName_.size() )
{
fatalErrorInFunction<<
"Number of surfaces is not equal to number of motion component names"<<endl;
fatalExit;
}
if(!createPropertyId())
{
fatalExit;
}
}
bool pFlow::geometry::beforeIteration()
{
zeroForce();
return true;
}
bool pFlow::geometry::iterate()
{
return true;
}
bool pFlow::geometry::afterIteration()
{
auto numTris = size();
auto& normalsD = normals().deviceViewAll();
auto& areaD = area().deviceViewAll();
auto& cForceD = contactForceWall_.deviceViewAll();
auto& sStressD = shearStressWall_.deviceViewAll();
auto& nStressD = normalStressWall_.deviceViewAll();
Kokkos::parallel_for(
"pFlow::geometry::afterIteration",
deviceRPolicyStatic(0, numTris),
LAMBDA_HD(uint32 i)
{
realx3 n = normalsD[i];
real A = max(areaD[i],smallValue);
realx3 nF = dot(cForceD[i],n)*n;
realx3 sF = cForceD[i] - nF;
nStressD[i] = nF/A;
sStressD[i] = sF/A;
});
Kokkos::fence();
return true;
}
bool pFlow::geometry::read(iIstream &is, const IOPattern &iop)
{
motionComponentName_.read(is, iop);
materialName_.read(is, iop);
if( readWholeObject_ )
{
return multiTriSurface::read(is, iop);
}
return true;
}
bool pFlow::geometry::write(iOstream &os, const IOPattern &iop) const
{
if( !motionComponentName_.write(os, iop) )
{
fatalErrorInFunction;
return false;
}
if( !materialName_.write(os,iop))
{
fatalErrorInFunction;
return false;
}
return multiTriSurface::write(os,iop);
}
pFlow::uniquePtr<pFlow::geometry>
pFlow::geometry::create
(
systemControl& control,
const property& prop
)
{
//
fileDictionary dict( motionModelFile__, control.time().geometry().path());
word model = dict.getVal<word>("motionModel");
auto geomModel = angleBracketsNames("geometry", model);
REPORT(1)<< "Selecting geometry model . . ."<<END_REPORT;
if( systemControlvCtorSelector_.search(geomModel) )
{
auto objPtr = systemControlvCtorSelector_[geomModel] (control, prop);
REPORT(2)<<"Model "<< Green_Text(geomModel)<<" is created.\n"<<END_REPORT;
return objPtr;
}
else
{
printKeys
(
fatalError << "Ctor Selector "<< Yellow_Text(geomModel) << " dose not exist. \n"
<<"Avaiable ones are: \n\n"
,
systemControlvCtorSelector_
);
fatalExit;
}
return nullptr;
}
pFlow::uniquePtr<pFlow::geometry>
pFlow::geometry::create
(
systemControl& control,
const property& prop,
multiTriSurface& surf,
const wordVector& motionCompName,
const wordVector& materialName,
const dictionary& motionDic
)
{
word model = motionDic.getVal<word>("motionModel");
auto geomModel = angleBracketsNames("geometry", model);
REPORT(1)<< "Selecting geometry model . . ."<<END_REPORT;
if( dictionaryvCtorSelector_.search(geomModel) )
{
auto objPtr = dictionaryvCtorSelector_[geomModel]
(
control,
prop,
surf,
motionCompName,
materialName,
motionDic
);
REPORT(2)<<"Model "<< Green_Text(geomModel)<<" is created.\n"<<END_REPORT;
return objPtr;
}
else
{
printKeys
(
fatalError << "Ctor Selector "<< Yellow_Text(geomModel) << " dose not exist. \n"
<<"Avaiable ones are: \n\n"
,
dictionaryvCtorSelector_
);
fatalExit;
}
return nullptr;
}
/*pFlow::geometry::geometry
(
systemControl& control,
const property& prop,
@ -224,52 +510,7 @@ pFlow::geometry::geometry
{}
pFlow::uniquePtr<pFlow::geometry>
pFlow::geometry::create
(
systemControl& control,
const property& prop
)
{
//motionModelFile__
auto motionDictPtr = IOobject::make<dictionary>
(
objectFile
(
motionModelFile__,
control.geometry().path(),
objectFile::READ_ALWAYS,
objectFile::WRITE_NEVER
),
motionModelFile__,
true
);
word model = motionDictPtr().getObject<dictionary>().getVal<word>("motionModel");
auto geomModel = angleBracketsNames("geometry", model);
REPORT(1)<< "Selecting geometry model . . ."<<endREPORT;
if( systemControlvCtorSelector_.search(geomModel) )
{
auto objPtr = systemControlvCtorSelector_[geomModel] (control, prop);
REPORT(2)<<"Model "<< greenText(geomModel)<<" is created.\n"<<endREPORT;
return objPtr;
}
else
{
printKeys
(
fatalError << "Ctor Selector "<< yellowText(geomModel) << " dose not exist. \n"
<<"Avaiable ones are: \n\n"
,
systemControlvCtorSelector_
);
fatalExit;
}
return nullptr;
}
bool pFlow::geometry::beforeIteration()
{
@ -296,48 +537,4 @@ bool pFlow::geometry::afterIteration()
Kokkos::fence();
return true;
}
pFlow::uniquePtr<pFlow::geometry>
pFlow::geometry::create(
systemControl& control,
const property& prop,
const dictionary& dict,
const multiTriSurface& triSurface,
const wordVector& motionCompName,
const wordVector& propName)
{
word model = dict.getVal<word>("motionModel");
auto geomModel = angleBracketsNames("geometry", model);
REPORT(1)<< "Selecting geometry model . . ."<<endREPORT;
if( dictionaryvCtorSelector_.search(geomModel) )
{
auto objPtr = dictionaryvCtorSelector_[geomModel]
(
control,
prop,
dict,
triSurface,
motionCompName,
propName
);
REPORT(2)<<"Model "<< greenText(geomModel)<<" is created.\n"<<endREPORT;
return objPtr;
}
else
{
printKeys
(
fatalError << "Ctor Selector "<< yellowText(geomModel) << " dose not exist. \n"
<<"Avaiable ones are: \n\n"
,
dictionaryvCtorSelector_
);
fatalExit;
}
return nullptr;
}
}*/

215
src/Geometry/geometry/geometry.hpp
View File

@ -23,13 +23,14 @@ Licence:
#include "virtualConstructor.hpp"
#include "demGeometry.hpp"
#include "demComponent.hpp"
#include "property.hpp"
#include "Fields.hpp"
#include "Vectors.hpp"
#include "multiTriSurface.hpp"
#include "triSurfaceFields.hpp"
#include "dictionary.hpp"
//#include "Fields.hpp"
//#include "Vectors.hpp"
namespace pFlow
{
@ -42,49 +43,50 @@ namespace pFlow
*/
class geometry
:
public demGeometry
public multiTriSurface,
public demComponent
{
protected:
private:
// - Protected members
/// Const reference to physical property of materials
const property& wallProperty_;
/// Repository to store geometry data at each simulation moment
repository& geometryRepository_;
/// All triangles in the set of wall surfaces
multiTriSurface& triSurface_;
/// The name of motion component of each wall surface
wordField& motionComponentName_;
wordField_H motionComponentName_{
"motionComponentName",
"motionComponentName"};
/// Material name of each wall surface
wordField& materialName_;
wordField_H materialName_{
"materialName",
"materialName"};
/// Property id of each triangle in the set of wall surfaces
int8TriSurfaceField_D& propertyId_;
uint32TriSurfaceField_D propertyId_;
/// Contact force on each triangle in the set of wall surfaces
realx3TriSurfaceField_D& contactForceWall_;
realx3TriSurfaceField_D contactForceWall_;
/// Stress on ech triangle in the set of wall surfaces
realx3TriSurfaceField_D& stressWall_;
/// Stress on each triangle in the set of wall surfaces
realx3TriSurfaceField_D normalStressWall_;
/// Stress on each triangle in the set of wall surfaces
realx3TriSurfaceField_D shearStressWall_;
bool readWholeObject_ = true;
// - Protected member functions
/// Find property id of each triangle based on the supplied material name
/// and the surface wall that the triangle belongs to.
bool findPropertyId();
bool createPropertyId();
/// Initialize contact force to zero
void zeroForce()
{
contactForceWall_.fill(zero3);
}
void zeroForce();
public:
/// Type info
@ -95,8 +97,15 @@ public:
/// Construct from controlSystem and property, for reading from file
geometry(systemControl& control, const property& prop);
geometry(systemControl& control,
const property& prop,
multiTriSurface& surf,
const wordVector& motionCompName,
const wordVector& materialName,
const dictionary& motionDict);
/// Construct from components
geometry(systemControl& control,
/*geometry(systemControl& control,
const property& prop,
const multiTriSurface& triSurface,
const wordVector& motionCompName,
@ -110,10 +119,10 @@ public:
const dictionary& dict,
const multiTriSurface& triSurface,
const wordVector& motionCompName,
const wordVector& propName);
const wordVector& propName);*/
/// Destructor
virtual ~geometry() = default;
~geometry()override = default;
/// Virtual constructor
create_vCtor
@ -132,156 +141,96 @@ public:
(
geometry,
dictionary,
(systemControl& control,
const property& prop,
const dictionary& dict,
const multiTriSurface& triSurface,
const wordVector& motionCompName,
const wordVector& propName),
(control, prop, dict, triSurface, motionCompName, propName)
(
systemControl& control,
const property& prop,
multiTriSurface& surf,
const wordVector& motionCompName,
const wordVector& materialName,
const dictionary& motionDic),
(control, prop, surf, motionCompName, materialName, motionDic)
);
//- Methods
/// Size of tri-surface
inline
auto size()const
//- Methods
inline
const auto& motionComponentName()const
{
return triSurface_.size();
return motionComponentName_;
}
/// Number of points in the set of surface walls
inline
auto numPoints()const
{
return triSurface_.numPoints();
}
/// Number of triangles in the set of surface walls
inline
auto numTriangles()const
{
return size();
}
/// Access to the points
inline
const auto& points()const
{
return triSurface_.points();
}
/// Access to the vertices
inline
const auto& vertices()const
{
return triSurface_.vertices();
}
/// Obtain an object for accessing triangles
inline auto getTriangleAccessor()const
{
return triSurface_.getTriangleAccessor();
}
/// Surface
inline auto& surface()
{
return triSurface_;
}
/// Surface
inline const auto& surface()const
{
return triSurface_;
}
/// Access to contact force
inline
realx3TriSurfaceField_D& contactForceWall()
auto& contactForceWall()
{
return contactForceWall_;
}
/// Access to contact force
inline
const realx3TriSurfaceField_D& contactForceWall() const
const auto& contactForceWall() const
{
return contactForceWall_;
}
/// Property ide of triangles
inline
const auto& propertyId()const
{
return propertyId_;
}
/// Access to property
inline const auto& wallProperty()const
{
return wallProperty_;
}
/// Owner repository
inline
const repository& owner()const
{
return geometryRepository_;
}
/// Owner repository
inline
repository& owner()
{
return geometryRepository_;
}
/// Path to the repository folder
inline auto path()
{
return owner().path();
}
/// The name of motion model
virtual
word motionModelTypeName()const = 0;
/// Motion model index of triangles
virtual
const int8Vector_HD& triMotionIndex() const =0;
const uint32Field_D& triMotionIndex() const =0;
/// Motion model index of points
virtual
const int8Vector_HD& pointMotionIndex()const = 0;
const uint32Field_D& pointMotionIndex()const = 0;
/// Property ide of triangles
const int8TriSurfaceField_D& propertyId() const
{
return propertyId_;
}
/// Operations before each iteration
bool beforeIteration() override;
/// Operations after each iteration
/// This is called in time loop. Perform the main calculations
/// when the component should evolve along time.
bool iterate() override;
/// This is called in time loop, after iterate.
bool afterIteration() override;
//- IO
bool read(iIstream& is, const IOPattern& iop) override;
/// write
bool write()const
{
return owner().write();
}
bool write( iOstream& os, const IOPattern& iop )const override;
//- Static members
static
uniquePtr<geometry> create(systemControl& control, const property& prop);
uniquePtr<geometry> create(
systemControl& control,
const property& prop);
static
uniquePtr<geometry> create(
systemControl& control,
const property& prop,
const dictionary& dict,
const multiTriSurface& triSurface,
const wordVector& motionCompName,
const wordVector& propName);
systemControl& control,
const property& prop,
multiTriSurface& surf,
const wordVector& motionCompName,
const wordVector& materialName,
const dictionary& motionDic);
};

294
src/Geometry/geometryMotion/geometryMotion.cpp
View File

@ -1,3 +1,4 @@
#include "geometryMotion.hpp"
/*------------------------------- phasicFlow ---------------------------------
O C enter of
O O E ngineering and
@ -21,42 +22,88 @@ Licence:
template<typename MotionModel>
bool pFlow::geometryMotion<MotionModel>::findMotionIndex()
{
motionIndex_.clear();
triMotionIndex_.reserve( this->surface().capacity() );
triMotionIndex_.clear();
ForAll( surfI, motionComponentName_)
if(motionComponentName().size() != numSurfaces() )
{
auto mName = motionComponentName_[surfI];
auto mInd = motionModel_.nameToIndex(mName);
motionIndex_.push_back(mInd);
// fill motionIndex for triangles of the surface
int32 surfSize = this->surface().surfNumTriangles(surfI);
for(int32 i=0; i<surfSize; i++)
fatalErrorInFunction<<
"size of motion component names in the triSurface is not"<<
" equal to size of number of sub-surfaces"<<endl;
return false;
}
uint32Vector surfMotionIndex("surfMotionIndex");
uint32Vector triMotionIndex("triMotionIndex");
uint32Vector pointMotionIndex("pointMotionIndex");
ForAll( surfI, motionComponentName())
{
auto mName = motionComponentName()[surfI];
uint32 mInd=0;
if( !motionModel_.nameToIndex(mName, mInd) )
{
triMotionIndex_.push_back(mInd);
fatalErrorInFunction<<
mName<< " does not exist in the list of motion names -> "<<
motionModel_.componentNames();
return false;
}
surfMotionIndex.push_back(mInd);
auto surfRange = subSurfaceRange(surfI);
for(uint32 i=0; i<surfRange.numElements(); i++)
{
triMotionIndex.push_back(mInd);
}
auto pointRange = subSurfacePointRange(surfI);
for(uint32 n=0; n<pointRange.numElements(); n++)
{
pointMotionIndex.push_back(mInd);
}
}
motionIndex_.syncViews();
triMotionIndex_.syncViews();
pointMotionIndex_.reserve(triSurface_.numPoints());
pointMotionIndex_.clear();
ForAll(surfI, motionIndex_)
{
auto nP = triSurface_.surfNumPoints(surfI);
for(int32 i=0; i<nP; i++)
{
pointMotionIndex_.push_back(motionIndex_[surfI]);
}
}
pointMotionIndex_.syncViews();
surfMotionIndex_.assign(surfMotionIndex);
triMotionIndex_.assign(triMotionIndex);
pointMotionIndex_.assign(pointMotionIndex);
return true;
}
template<typename MotionModel>
bool pFlow::geometryMotion<MotionModel>::moveGeometry()
{
uint32 iter = this->currentIter();
real t = this->currentTime();
real dt = this->dt();
auto mModel = motionModel_.getModelInterface(iter, t, dt);
auto& pointMIndexD= pointMotionIndex_.deviceViewAll();
auto& pointsD = points().deviceViewAll();
Kokkos::parallel_for(
"geometryMotion<MotionModel>::movePoints",
deviceRPolicyStatic(0, numPoints()),
LAMBDA_HD(uint32 i){
auto newPos = mModel.transferPoint(pointMIndexD[i], pointsD[i], dt);
pointsD[i] = newPos;
});
Kokkos::fence();
// move the motion components
motionModel_.move(iter, t,dt);
// end of calculations
return true;
}
template<typename MotionModel>
pFlow::geometryMotion<MotionModel>::geometryMotion
(
@ -65,144 +112,75 @@ pFlow::geometryMotion<MotionModel>::geometryMotion
)
:
geometry(control, prop),
motionModel_(
this->owner().template emplaceObject<MotionModel>(
objectFile(
motionModelFile__,
"",
objectFile::READ_ALWAYS,
objectFile::WRITE_ALWAYS
)
)
motionModel_
(
objectFile
(
motionModelFile__,
"",
objectFile::READ_ALWAYS,
objectFile::WRITE_ALWAYS
),
owner()
),
moveGeomTimer_("move geometry", &this->timers())
{
findMotionIndex();
}
template<typename MotionModel>
pFlow::geometryMotion<MotionModel>::geometryMotion
(
systemControl& control,
const property& prop,
const multiTriSurface& triSurface,
const wordVector& motionCompName,
const wordVector& propName,
const MotionModel& motionModel
)
:
geometry(
control,
prop,
triSurface,
motionCompName,
propName
),
motionModel_(
this->owner().template emplaceObject<MotionModel>(
objectFile(
motionModelFile__,
"",
objectFile::READ_NEVER,
objectFile::WRITE_ALWAYS
),
motionModel
)
),
moveGeomTimer_("move geometry", &this->timers())
{
findMotionIndex();
}
template<typename MotionModel>
pFlow::geometryMotion<MotionModel>::geometryMotion
(
systemControl& control,
const property& prop,
const dictionary& dict,
const multiTriSurface& triSurface,
const wordVector& motionCompName,
const wordVector& propName
)
:
geometry(
control,
prop,
dict,
triSurface,
motionCompName,
propName
),
motionModel_(
this->owner().template emplaceObject<MotionModel>(
objectFile(
motionModelFile__,
"",
objectFile::READ_NEVER,
objectFile::WRITE_ALWAYS
),
dict
)
),
moveGeomTimer_("move geometry", &this->timers())
{
findMotionIndex();
}
template<typename MotionModel>
bool pFlow::geometryMotion<MotionModel>::beforeIteration()
{
geometry::beforeIteration();
return true;
}
template<typename MotionModel>
bool pFlow::geometryMotion<MotionModel>::iterate()
{
if( motionModel_.isMoving() )
if(!findMotionIndex())
{
moveGeomTimer_.start();
moveGeometry();
moveGeomTimer_.end();
fatalExit;
}
return true;
}
template<typename MotionModel>
bool pFlow::geometryMotion<MotionModel>::afterIteration()
template <typename MotionModelType>
pFlow::geometryMotion<MotionModelType>::geometryMotion
(
systemControl &control,
const property &prop,
multiTriSurface &surf,
const wordVector &motionCompName,
const wordVector &materialName,
const dictionary &motionDict
)
:
geometry
(
control,
prop,
surf,
motionCompName,
materialName,
motionDict
),
motionModel_
(
objectFile
(
motionModelFile__,
"",
objectFile::READ_NEVER,
objectFile::WRITE_ALWAYS
),
motionDict,
owner()
),
moveGeomTimer_("move geometry", &this->timers())
{
geometry::afterIteration();
return true;
}
template<typename MotionModel>
bool pFlow::geometryMotion<MotionModel>::moveGeometry()
{
real dt = this->dt();
real t = this->currentTime();
auto pointMIndex= pointMotionIndex_.deviceVector();
auto mModel = motionModel_.getModel(t);
realx3* points = triSurface_.pointsData_D();
auto numPoints = triSurface_.numPoints();
if(!findMotionIndex())
{
fatalExit;
}
}
Kokkos::parallel_for(
"geometryMotion<MotionModel>::movePoints",
numPoints,
LAMBDA_HD(int32 i){
auto newPos = mModel.transferPoint(pointMIndex[i], points[i], dt);
points[i] = newPos;
});
Kokkos::fence();
// move the motion components
motionModel_.move(t,dt);
// end of calculations
moveGeomTimer_.end();
return true;
}
template<typename MotionModel>
bool pFlow::geometryMotion<MotionModel>::iterate()
{
if( motionModel_.isMoving() )
{
moveGeomTimer_.start();
moveGeometry();
this->calculateNormals();
moveGeomTimer_.end();
}
return true;
}

60
src/Geometry/geometryMotion/geometryMotion.hpp
View File

@ -20,9 +20,8 @@ Licence:
#ifndef __geometryMotion_hpp__
#define __geometryMotion_hpp__
#include "vocabs.hpp"
#include "geometry.hpp"
#include "VectorDuals.hpp"
namespace pFlow
{
@ -39,21 +38,23 @@ class geometryMotion
{
public:
using MotionModel = MotionModelType;
using MotionModel = MotionModelType;
protected:
using ModelComponent = typename MotionModelType::ModelComponent;
private:
/// Ref to motion model
MotionModel& motionModel_;
MotionModelType motionModel_;
/// motion indext mapped on each surface
int32Vector_HD motionIndex_;
uint32Field_D surfMotionIndex_{"triMotionIndex"};
/// motion index mapped on each triangle
int8Vector_HD triMotionIndex_;
uint32Field_D triMotionIndex_ {"surfMotionIndex"};
/// motion index mapped on each point
int8Vector_HD pointMotionIndex_;
uint32Field_D pointMotionIndex_{"pointMotionIndex"};
/// timer for moveGeometry
Timer moveGeomTimer_;
@ -61,32 +62,25 @@ protected:
/// determine the motion index of each triangle
bool findMotionIndex();
/// Move geometry
bool moveGeometry();
public:
/// Type info
TypeInfoTemplate("geometry", MotionModel);
TypeInfoTemplate11("geometry", ModelComponent);
// - Constructors
geometryMotion(systemControl& control, const property& prop);
geometryMotion(
systemControl& control,
systemControl& control,
const property& prop,
const multiTriSurface& triSurface,
multiTriSurface& surf,
const wordVector& motionCompName,
const wordVector& propName,
const MotionModel& motionModel);
/// Construct from components and dictionary that contains
/// motionModel
geometryMotion(systemControl& control,
const property& prop,
const dictionary& dict,
const multiTriSurface& triSurface,
const wordVector& motionCompName,
const wordVector& propName);
const wordVector& materialName,
const dictionary& motionDict);
/// Add virtual constructor
add_vCtor
@ -107,9 +101,9 @@ public:
// - Methods
/// Obtain motion model at time t
auto getModel(real t)const
auto getModel(uint32 iter, real t, real dt)const
{
return motionModel_.getModel(t);
return motionModel_.getModelInterface(iter, t, dt);
}
/// TypeName / TypeInfo of motion model
@ -119,28 +113,21 @@ public:
}
/// Access to motion model index of triangles
const int8Vector_HD& triMotionIndex()const override
const uint32Field_D& triMotionIndex()const override
{
return triMotionIndex_;
}
/// Access to motion model index of points
const int8Vector_HD& pointMotionIndex()const override
const uint32Field_D& pointMotionIndex()const override
{
return pointMotionIndex_;
}
/// Operations before each iteration
bool beforeIteration() override;
/// Iterate geometry one time step
bool iterate() override ;
/// Operations after each iteration
bool afterIteration() override;
/// Move geometry
bool moveGeometry();
};
@ -148,9 +135,6 @@ public:
#include "geometryMotion.cpp"
#ifndef BUILD_SHARED_LIBS
#include "geometryMotionsInstantiate.cpp"
#endif
#endif //__geometryMotion_hpp__

12
src/Geometry/geometryMotion/geometryMotions.cpp
View File

@ -17,9 +17,13 @@ Licence:
implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
-----------------------------------------------------------------------------*/
#include "geometryMotions.hpp"
#ifdef BUILD_SHARED_LIBS
#include "geometryMotionsInstantiate.cpp"
#endif
template class pFlow::geometryMotion<pFlow::vibratingMotion>;
template class pFlow::geometryMotion<pFlow::rotatingAxisMotion>;
template class pFlow::geometryMotion<pFlow::stationaryWall>;
//template class pFlow::geometryMotion<pFlow::multiRotatingAxisMotion>;

13
src/Geometry/geometryMotion/geometryMotions.hpp
View File

@ -22,22 +22,23 @@ Licence:
#define __geometryMotions_hpp__
#include "geometryMotion.hpp"
#include "fixedWall.hpp"
#include "stationaryWall.hpp"
#include "rotatingAxisMotion.hpp"
#include "multiRotatingAxisMotion.hpp"
//#include "multiRotatingAxisMotion.hpp"
#include "vibratingMotion.hpp"
namespace pFlow
{
typedef geometryMotion<vibratingMotion> vibratingMotionGeometry;
using vibratingMotionGeometry = geometryMotion<vibratingMotion>;
typedef geometryMotion<rotatingAxisMotion> rotationAxisMotionGeometry;
using rotationAxisMotionGeometry = geometryMotion<rotatingAxisMotion>;
typedef geometryMotion<multiRotatingAxisMotion> multiRotationAxisMotionGeometry;
using stationaryGeometry = geometryMotion<stationaryWall>;
//typedef geometryMotion<multiRotatingAxisMotion> multiRotationAxisMotionGeometry;
typedef geometryMotion<fixedWall> fixedGeometry;

163
src/Integration/AdamsBashforth2/AdamsBashforth2.cpp
View File

@ -19,59 +19,146 @@ Licence:
-----------------------------------------------------------------------------*/
#include "AdamsBashforth2.hpp"
#include "pointStructure.hpp"
#include "Time.hpp"
#include "vocabs.hpp"
namespace pFlow
{
/// Range policy for integration kernel (alias)
using rpIntegration = Kokkos::RangePolicy<
DefaultExecutionSpace,
Kokkos::Schedule<Kokkos::Static>,
Kokkos::IndexType<uint32>
>;
bool intAllActive(
real dt,
realx3Field_D& y,
realx3PointField_D& dy,
realx3PointField_D& dy1)
{
auto d_dy = dy.deviceView();
auto d_y = y.deviceView();
auto d_dy1= dy1.deviceView();
auto activeRng = dy1.activeRange();
Kokkos::parallel_for(
"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_dy1[i] = d_dy[i];
});
Kokkos::fence();
return true;
}
bool intScattered
(
real dt,
realx3Field_D& y,
realx3PointField_D& dy,
realx3PointField_D& dy1
)
{
auto d_dy = dy.deviceView();
auto d_y = y.deviceView();
auto d_dy1 = dy1.deviceView();
auto activeRng = dy1.activeRange();
const auto& activeP = dy1.activePointsMaskDevice();
Kokkos::parallel_for(
"AdamsBashforth2::correct",
rpIntegration (activeRng.start(), activeRng.end()),
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_dy1[i] = d_dy[i];
}
});
Kokkos::fence();
return true;
}
}
//const real AB2_coef[] = { 3.0 / 2.0, 1.0 / 2.0};
pFlow::AdamsBashforth2::AdamsBashforth2
(
const word& baseName,
repository& owner,
const pointStructure& pStruct,
const word& method
pointStructure& pStruct,
const word& method,
const realx3Field_D& initialValField
)
:
integration(baseName, owner, pStruct, method),
dy1_(
owner.emplaceObject<realx3PointField_D>(
objectFile(
groupNames(baseName,"dy1"),
"",
objectFile::READ_IF_PRESENT,
objectFile::WRITE_ALWAYS),
pStruct,
zero3))
{
}
integration(baseName, pStruct, method, initialValField),
realx3PointField_D
(
objectFile
(
groupNames(baseName,"dy1"),
pStruct.time().integrationFolder(),
objectFile::READ_IF_PRESENT,
objectFile::WRITE_ALWAYS
),
pStruct,
zero3,
zero3
)
{}
bool pFlow::AdamsBashforth2::predict
(
real UNUSED(dt),
realx3Vector_D& UNUSED(y),
realx3Vector_D& UNUSED(dy)
realx3PointField_D& UNUSED(y),
realx3PointField_D& UNUSED(dy)
)
{
return true;
}
bool pFlow::AdamsBashforth2::predict
(
real dt,
realx3Field_D &y,
realx3PointField_D &dy
)
{
return true;
}
bool pFlow::AdamsBashforth2::correct
(
real dt,
realx3Vector_D& y,
realx3Vector_D& dy
realx3PointField_D& y,
realx3PointField_D& dy
)
{
if(this->pStruct().allActive())
return correct(dt, y.field(), dy);
}
bool pFlow::AdamsBashforth2::correct(real dt, realx3Field_D &y, realx3PointField_D &dy)
{
auto& dy1l = dy1();
if(dy1l.isAllActive())
{
return intAll(dt, y, dy, this->pStruct().activeRange());
return intAllActive(dt, y, dy, dy1l);
}
else
{
return intRange(dt, y, dy, this->pStruct().activePointsMaskD());
return intScattered(dt, y, dy, dy1l);
}
return true;
return false;
}
bool pFlow::AdamsBashforth2::setInitialVals(
@ -81,25 +168,3 @@ bool pFlow::AdamsBashforth2::setInitialVals(
return true;
}
bool pFlow::AdamsBashforth2::intAll(
real dt,
realx3Vector_D& y,
realx3Vector_D& dy,
range activeRng)
{
auto d_dy = dy.deviceVectorAll();
auto d_y = y.deviceVectorAll();
auto d_dy1= dy1_.deviceVectorAll();
Kokkos::parallel_for(
"AdamsBashforth2::correct",
rpIntegration (activeRng.first, activeRng.second),
LAMBDA_HD(int32 i){
d_y[i] += dt*(static_cast<real>(3.0 / 2.0) * d_dy[i] - static_cast<real>(1.0 / 2.0) * d_dy1[i]);
d_dy1[i] = d_dy[i];
});
Kokkos::fence();
return true;
}

102
src/Integration/AdamsBashforth2/AdamsBashforth2.hpp
View File

@ -36,20 +36,16 @@ namespace pFlow
*/
class AdamsBashforth2
:
public integration
public integration,
public realx3PointField_D
{
protected:
private:
/// dy at t-dt
realx3PointField_D& dy1_;
/// Range policy for integration kernel (alias)
using rpIntegration = Kokkos::RangePolicy<
DefaultExecutionSpace,
Kokkos::Schedule<Kokkos::Static>,
Kokkos::IndexType<int32>
>;
auto& dy1()
{
return static_cast<realx3PointField_D&>(*this);
}
public:
/// Type info
@ -60,17 +56,12 @@ public:
/// Construct from components
AdamsBashforth2(
const word& baseName,
repository& owner,
const pointStructure& pStruct,
const word& method);
uniquePtr<integration> clone()const override
{
return makeUnique<AdamsBashforth2>(*this);
}
pointStructure& pStruct,
const word& method,
const realx3Field_D& initialValField);
/// Destructor
virtual ~AdamsBashforth2()=default;
~AdamsBashforth2()final = default;
/// Add this to the virtual constructor table
add_vCtor(
@ -80,71 +71,44 @@ public:
// - Methods
/// return integration method
word method()const override
{
return "AdamsBashforth2";
}
bool predict(
real UNUSED(dt),
realx3Vector_D& UNUSED(y),
realx3Vector_D& UNUSED(dy)) override;
realx3PointField_D& UNUSED(y),
realx3PointField_D& UNUSED(dy)) final;
bool predict(
real dt,
realx3Field_D& y,
realx3PointField_D& dy) final;
bool correct(
real dt,
realx3Vector_D& y,
realx3Vector_D& dy) override;
realx3PointField_D& y,
realx3PointField_D& dy) final;
bool correct(
real dt,
realx3Field_D& y,
realx3PointField_D& dy);
bool setInitialVals(
const int32IndexContainer& newIndices,
const realx3Vector& y) override;
const realx3Vector& y) final;
bool needSetInitialVals()const override
bool needSetInitialVals()const final
{
return false;
}
/// 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 intRange(
real dt,
realx3Vector_D& y,
realx3Vector_D& dy,
activeFunctor activeP );
};
template<typename activeFunctor>
bool pFlow::AdamsBashforth2::intRange(
real dt,
realx3Vector_D& y,
realx3Vector_D& dy,
activeFunctor activeP )
{
auto d_dy = dy.deviceVectorAll();
auto d_y = y.deviceVectorAll();
auto d_dy1= dy1_.deviceVectorAll();
auto activeRng = activeP.activeRange();
Kokkos::parallel_for(
"AdamsBashforth2::correct",
rpIntegration (activeRng.first, activeRng.second),
LAMBDA_HD(int32 i){
if( activeP(i))
{
d_y[i] += dt*(static_cast<real>(3.0 / 2.0) * d_dy[i] - static_cast<real>(1.0 / 2.0) * d_dy1[i]);
d_dy1[i] = d_dy[i];
}
});
Kokkos::fence();
return true;
}
} // pFlow

6
src/Integration/AdamsBashforth3/AdamsBashforth3.cpp
View File

@ -89,9 +89,9 @@ bool pFlow::AdamsBashforth3::intAll(
realx3Vector_D& dy,
range activeRng)
{
auto d_dy = dy.deviceVectorAll();
auto d_y = y.deviceVectorAll();
auto d_history = history_.deviceVectorAll();
auto d_dy = dy.deviceViewAll();
auto d_y = y.deviceViewAll();
auto d_history = history_.deviceViewAll();
Kokkos::parallel_for(
"AdamsBashforth3::correct",

6
src/Integration/AdamsBashforth3/AdamsBashforth3.hpp
View File

@ -160,9 +160,9 @@ bool pFlow::AdamsBashforth3::intRange(
realx3Vector_D& dy,
activeFunctor activeP )
{
auto d_dy = dy.deviceVectorAll();
auto d_y = y.deviceVectorAll();
auto d_history = history_.deviceVectorAll();
auto d_dy = dy.deviceViewAll();
auto d_y = y.deviceViewAll();
auto d_history = history_.deviceViewAll();
auto activeRng = activeP.activeRange();
Kokkos::parallel_for(

6
src/Integration/AdamsBashforth4/AdamsBashforth4.cpp
View File

@ -89,9 +89,9 @@ bool pFlow::AdamsBashforth4::intAll(
realx3Vector_D& dy,
range activeRng)
{
auto d_dy = dy.deviceVectorAll();
auto d_y = y.deviceVectorAll();
auto d_history = history_.deviceVectorAll();
auto d_dy = dy.deviceViewAll();
auto d_y = y.deviceViewAll();
auto d_history = history_.deviceViewAll();
Kokkos::parallel_for(
"AdamsBashforth4::correct",

6
src/Integration/AdamsBashforth4/AdamsBashforth4.hpp
View File

@ -165,9 +165,9 @@ bool pFlow::AdamsBashforth4::intRange(
realx3Vector_D& dy,
activeFunctor activeP )
{
auto d_dy = dy.deviceVectorAll();
auto d_y = y.deviceVectorAll();
auto d_history = history_.deviceVectorAll();
auto d_dy = dy.deviceViewAll();
auto d_y = y.deviceViewAll();
auto d_history = history_.deviceViewAll();
auto activeRng = activeP.activeRange();
Kokkos::parallel_for(

6
src/Integration/AdamsBashforth5/AdamsBashforth5.cpp
View File

@ -89,9 +89,9 @@ bool pFlow::AdamsBashforth5::intAll(
realx3Vector_D& dy,
range activeRng)
{
auto d_dy = dy.deviceVectorAll();
auto d_y = y.deviceVectorAll();
auto d_history = history_.deviceVectorAll();
auto d_dy = dy.deviceViewAll();
auto d_y = y.deviceViewAll();
auto d_history = history_.deviceViewAll();
Kokkos::parallel_for(
"AdamsBashforth5::correct",

6
src/Integration/AdamsBashforth5/AdamsBashforth5.hpp
View File

@ -165,9 +165,9 @@ bool pFlow::AdamsBashforth5::intRange(
realx3Vector_D& dy,
activeFunctor activeP )
{
auto d_dy = dy.deviceVectorAll();
auto d_y = y.deviceVectorAll();
auto d_history = history_.deviceVectorAll();
auto d_dy = dy.deviceViewAll();
auto d_y = y.deviceViewAll();
auto d_history = history_.deviceViewAll();
auto activeRng = activeP.activeRange();
Kokkos::parallel_for(

20
src/Integration/AdamsMoulton3/AdamsMoulton3.cpp
View File

@ -124,11 +124,11 @@ bool pFlow::AdamsMoulton3::predictAll(
range activeRng)
{
auto d_dy = dy.deviceVectorAll();
auto d_y = y.deviceVectorAll();
auto d_y0 = y0_.deviceVectorAll();
auto d_dy0 = dy0_.deviceVectorAll();
auto d_dy1= dy1_.deviceVectorAll();
auto d_dy = dy.deviceViewAll();
auto d_y = y.deviceViewAll();
auto d_y0 = y0_.deviceViewAll();
auto d_dy0 = dy0_.deviceViewAll();
auto d_dy1= dy1_.deviceViewAll();
Kokkos::parallel_for(
"AdamsMoulton3::predict",
@ -150,12 +150,12 @@ bool pFlow::AdamsMoulton3::intAll(
range activeRng)
{
auto d_dy = dy.deviceVectorAll();
auto d_y = y.deviceVectorAll();
auto d_dy = dy.deviceViewAll();
auto d_y = y.deviceViewAll();
auto d_dy0 = dy0_.deviceVectorAll();
auto d_y0 = y0_.deviceVectorAll();
auto d_dy1 = dy1_.deviceVectorAll();
auto d_dy0 = dy0_.deviceViewAll();
auto d_y0 = y0_.deviceViewAll();
auto d_dy1 = dy1_.deviceViewAll();
Kokkos::parallel_for(
"AdamsMoulton3::correct",

20
src/Integration/AdamsMoulton3/AdamsMoulton3.hpp
View File

@ -144,11 +144,11 @@ bool AdamsMoulton3::predictRange(
realx3Vector_D& dy,
activeFunctor activeP)
{
auto d_dy = dy.deviceVectorAll();
auto d_y = y.deviceVectorAll();
auto d_y0 = y0_.deviceVectorAll();
auto d_dy0 = dy0_.deviceVectorAll();
auto d_dy1= dy1_.deviceVectorAll();
auto d_dy = dy.deviceViewAll();
auto d_y = y.deviceViewAll();
auto d_y0 = y0_.deviceViewAll();
auto d_dy0 = dy0_.deviceViewAll();
auto d_dy1= dy1_.deviceViewAll();
auto activeRng = activeP.activeRange();
@ -182,12 +182,12 @@ bool pFlow::AdamsMoulton3::intRange(
activeFunctor activeP)
{
auto d_dy = dy.deviceVectorAll();
auto d_y = y.deviceVectorAll();
auto d_dy = dy.deviceViewAll();
auto d_y = y.deviceViewAll();
auto d_dy0 = dy0_.deviceVectorAll();
auto d_y0 = y0_.deviceVectorAll();
auto d_dy1 = dy1_.deviceVectorAll();
auto d_dy0 = dy0_.deviceViewAll();
auto d_y0 = y0_.deviceViewAll();
auto d_dy1 = dy1_.deviceViewAll();
auto activeRng = activeP.activeRange();

24
src/Integration/AdamsMoulton4/AdamsMoulton4.cpp
View File

@ -135,13 +135,13 @@ bool pFlow::AdamsMoulton4::predictAll(
range activeRng)
{
auto d_dy = dy.deviceVectorAll();
auto d_y = y.deviceVectorAll();
auto d_dy = dy.deviceViewAll();
auto d_y = y.deviceViewAll();
auto d_y0 = y0_.deviceVectorAll();
auto d_dy0 = dy0_.deviceVectorAll();
auto d_dy1 = dy1_.deviceVectorAll();
auto d_dy2 = dy2_.deviceVectorAll();
auto d_y0 = y0_.deviceViewAll();
auto d_dy0 = dy0_.deviceViewAll();
auto d_dy1 = dy1_.deviceViewAll();
auto d_dy2 = dy2_.deviceViewAll();
Kokkos::parallel_for(
"AdamsMoulton4::predict",
@ -165,13 +165,13 @@ bool pFlow::AdamsMoulton4::intAll(
range activeRng)
{
auto d_dy = dy.deviceVectorAll();
auto d_y = y.deviceVectorAll();
auto d_dy = dy.deviceViewAll();
auto d_y = y.deviceViewAll();
auto d_dy0 = dy0_.deviceVectorAll();
auto d_y0 = y0_.deviceVectorAll();
auto d_dy1 = dy1_.deviceVectorAll();
auto d_dy2 = dy2_.deviceVectorAll();
auto d_dy0 = dy0_.deviceViewAll();
auto d_y0 = y0_.deviceViewAll();
auto d_dy1 = dy1_.deviceViewAll();
auto d_dy2 = dy2_.deviceViewAll();
Kokkos::parallel_for(
"AdamsMoulton4::correct",

24
src/Integration/AdamsMoulton4/AdamsMoulton4.hpp
View File

@ -147,13 +147,13 @@ bool AdamsMoulton4::predictRange(
realx3Vector_D& dy,
activeFunctor activeP )
{
auto d_dy = dy.deviceVectorAll();
auto d_y = y.deviceVectorAll();
auto d_dy = dy.deviceViewAll();
auto d_y = y.deviceViewAll();
auto d_y0 = y0_.deviceVectorAll();
auto d_dy0 = dy0_.deviceVectorAll();
auto d_dy1 = dy1_.deviceVectorAll();
auto d_dy2 = dy2_.deviceVectorAll();
auto d_y0 = y0_.deviceViewAll();
auto d_dy0 = dy0_.deviceViewAll();
auto d_dy1 = dy1_.deviceViewAll();
auto d_dy2 = dy2_.deviceViewAll();
auto activeRng = activeP.activeRange();
@ -185,13 +185,13 @@ bool pFlow::AdamsMoulton4::intRange(
activeFunctor activeP )
{
auto d_dy = dy.deviceVectorAll();
auto d_y = y.deviceVectorAll();
auto d_dy = dy.deviceViewAll();
auto d_y = y.deviceViewAll();
auto d_dy0 = dy0_.deviceVectorAll();
auto d_y0 = y0_.deviceVectorAll();
auto d_dy1 = dy1_.deviceVectorAll();
auto d_dy2 = dy2_.deviceVectorAll();
auto d_dy0 = dy0_.deviceViewAll();
auto d_y0 = y0_.deviceViewAll();
auto d_dy1 = dy1_.deviceViewAll();
auto d_dy2 = dy2_.deviceViewAll();
auto activeRng = activeP.activeRange();

28
src/Integration/AdamsMoulton5/AdamsMoulton5.cpp
View File

@ -145,14 +145,14 @@ bool pFlow::AdamsMoulton5::predictAll(
range activeRng)
{
auto d_dy = dy.deviceVectorAll();
auto d_y = y.deviceVectorAll();
auto d_dy = dy.deviceViewAll();
auto d_y = y.deviceViewAll();
auto d_y0 = y0_.deviceVectorAll();
auto d_dy0 = dy0_.deviceVectorAll();
auto d_dy1 = dy1_.deviceVectorAll();
auto d_dy2 = dy2_.deviceVectorAll();
auto d_dy3 = dy3_.deviceVectorAll();
auto d_y0 = y0_.deviceViewAll();
auto d_dy0 = dy0_.deviceViewAll();
auto d_dy1 = dy1_.deviceViewAll();
auto d_dy2 = dy2_.deviceViewAll();
auto d_dy3 = dy3_.deviceViewAll();
Kokkos::parallel_for(
"AdamsMoulton5::predict",
@ -178,14 +178,14 @@ bool pFlow::AdamsMoulton5::intAll(
range activeRng)
{
auto d_dy = dy.deviceVectorAll();
auto d_y = y.deviceVectorAll();
auto d_dy = dy.deviceViewAll();
auto d_y = y.deviceViewAll();
auto d_dy0 = dy0_.deviceVectorAll();
auto d_y0 = y0_.deviceVectorAll();
auto d_dy1 = dy1_.deviceVectorAll();
auto d_dy2 = dy2_.deviceVectorAll();
auto d_dy3 = dy3_.deviceVectorAll();
auto d_dy0 = dy0_.deviceViewAll();
auto d_y0 = y0_.deviceViewAll();
auto d_dy1 = dy1_.deviceViewAll();
auto d_dy2 = dy2_.deviceViewAll();
auto d_dy3 = dy3_.deviceViewAll();
Kokkos::parallel_for(
"AdamsMoulton5::correct",

28
src/Integration/AdamsMoulton5/AdamsMoulton5.hpp
View File

@ -150,14 +150,14 @@ bool AdamsMoulton5::predictRange(
realx3Vector_D& dy,
activeFunctor activeP )
{
auto d_dy = dy.deviceVectorAll();
auto d_y = y.deviceVectorAll();
auto d_dy = dy.deviceViewAll();
auto d_y = y.deviceViewAll();
auto d_y0 = y0_.deviceVectorAll();
auto d_dy0 = dy0_.deviceVectorAll();
auto d_dy1 = dy1_.deviceVectorAll();
auto d_dy2 = dy2_.deviceVectorAll();
auto d_dy3 = dy3_.deviceVectorAll();
auto d_y0 = y0_.deviceViewAll();
auto d_dy0 = dy0_.deviceViewAll();
auto d_dy1 = dy1_.deviceViewAll();
auto d_dy2 = dy2_.deviceViewAll();
auto d_dy3 = dy3_.deviceViewAll();
auto activeRng = activeP.activeRange();
@ -189,14 +189,14 @@ bool pFlow::AdamsMoulton5::intRange(
activeFunctor activeP )
{
auto d_dy = dy.deviceVectorAll();
auto d_y = y.deviceVectorAll();
auto d_dy = dy.deviceViewAll();
auto d_y = y.deviceViewAll();
auto d_dy0 = dy0_.deviceVectorAll();
auto d_y0 = y0_.deviceVectorAll();
auto d_dy1 = dy1_.deviceVectorAll();
auto d_dy2 = dy2_.deviceVectorAll();
auto d_dy3 = dy3_.deviceVectorAll();
auto d_dy0 = dy0_.deviceViewAll();
auto d_y0 = y0_.deviceViewAll();
auto d_dy1 = dy1_.deviceViewAll();
auto d_dy2 = dy2_.deviceViewAll();
auto d_dy3 = dy3_.deviceViewAll();
auto activeRng = activeP.activeRange();

12
src/Integration/CMakeLists.txt
View File

@ -1,13 +1,13 @@
list(APPEND SourceFiles
integration/integration.cpp
AdamsBashforth5/AdamsBashforth5.cpp
AdamsBashforth4/AdamsBashforth4.cpp
AdamsBashforth3/AdamsBashforth3.cpp
AdamsBashforth2/AdamsBashforth2.cpp
AdamsMoulton3/AdamsMoulton3.cpp
AdamsMoulton4/AdamsMoulton4.cpp
AdamsMoulton5/AdamsMoulton5.cpp
#AdamsBashforth5/AdamsBashforth5.cpp
#AdamsBashforth4/AdamsBashforth4.cpp
#AdamsBashforth3/AdamsBashforth3.cpp
#AdamsMoulton3/AdamsMoulton3.cpp
#AdamsMoulton4/AdamsMoulton4.cpp
#AdamsMoulton5/AdamsMoulton5.cpp
)
set(link_libs Kokkos::kokkos phasicFlow)

30
src/Integration/integration/integration.cpp
View File

@ -19,33 +19,35 @@ Licence:
-----------------------------------------------------------------------------*/
#include "integration.hpp"
#include "pointStructure.hpp"
#include "repository.hpp"
pFlow::integration::integration
(
const word& baseName,
repository& owner,
const pointStructure& pStruct,
const word& method
pointStructure& pStruct,
const word&,
const realx3Field_D&
)
:
owner_(owner),
baseName_(baseName),
pStruct_(pStruct)
{
CONSUME_PARAM(method);
}
owner_(*pStruct.owner()),
pStruct_(pStruct),
baseName_(baseName)
{}
pFlow::uniquePtr<pFlow::integration>
pFlow::integration::create(
pFlow::integration::create
(
const word& baseName,
repository& owner,
const pointStructure& pStruct,
const word& method)
pointStructure& pStruct,
const word& method,
const realx3Field_D& initialValField
)
{
if( wordvCtorSelector_.search(method) )
{
return wordvCtorSelector_[method] (baseName, owner, pStruct, method);
return wordvCtorSelector_[method] (baseName, pStruct, method, initialValField);
}
else
{

62
src/Integration/integration/integration.hpp
View File

@ -23,14 +23,16 @@ Licence:
#include "virtualConstructor.hpp"
#include "Vectors.hpp"
#include "pointStructure.hpp"
#include "repository.hpp"
#include "pointFields.hpp"
namespace pFlow
{
// - Forward
class repository;
class pointStructure;
/**
* Base class for integrating the first order ODE (IVP)
*
@ -48,19 +50,19 @@ namespace pFlow
*/
class integration
{
protected:
private:
// - Protected data members
/// The owner repository that all fields are storred in
repository& owner_;
/// The base name for integration
const word baseName_;
/// A reference to pointStructure
const pointStructure& pStruct_;
/// The base name for integration
const word baseName_;
public:
/// Type info
@ -72,9 +74,9 @@ public:
/// Construct from components
integration(
const word& baseName,
repository& owner,
const pointStructure& pStruct,
const word& method);
pointStructure& pStruct,
const word& method,
const realx3Field_D& initialValField);
/// Copy constructor
integration(const integration&) = default;
@ -88,22 +90,22 @@ public:
/// Move assignment
integration& operator = (integration&&) = default;
/// Polymorphic copy/cloning
virtual
uniquePtr<integration> clone()const=0;
/// Destructor
virtual ~integration()=default;
/// Add a virtual constructor
create_vCtor(
create_vCtor
(
integration,
word,
(const word& baseName,
repository& owner,
const pointStructure& pStruct,
const word& method),
(baseName, owner, pStruct, method) );
(
const word& baseName,
pointStructure& pStruct,
const word& method,
const realx3Field_D& initialValField
),
(baseName, pStruct, method, initialValField)
);
// - Methods
@ -129,13 +131,23 @@ public:
return owner_;
}
/// return integration method
virtual
word method()const = 0 ;
/// Prediction step in integration
virtual
bool predict(real dt, realx3Vector_D& y, realx3Vector_D& dy) = 0;
bool predict(real dt, realx3PointField_D& y, realx3PointField_D& dy) = 0;
virtual
bool predict(real dt, realx3Field_D& y, realx3PointField_D& dy) = 0;
/// Correction/main integration step
virtual
bool correct(real dt, realx3Vector_D& y, realx3Vector_D& dy) = 0;
bool correct(real dt, realx3PointField_D& y, realx3PointField_D& dy) = 0;
virtual
bool correct(real dt, realx3Field_D& y, realx3PointField_D& dy) = 0;
/// Set the initial values for new indices
virtual
@ -152,9 +164,9 @@ public:
static
uniquePtr<integration> create(
const word& baseName,
repository& owner,
const pointStructure& pStruct,
const word& method);
pointStructure& pStruct,
const word& method,
const realx3Field_D& initialValField);
}; // integration

5
src/Interaction/CMakeLists.txt
View File

@ -1,5 +1,10 @@
set(SourceFiles
contactSearch/methods/cellBased/NBS/mapperNBS.cpp
contactSearch/methods/cellBased/NBS/mapperNBSKernels.cpp
contactSearch/methods/cellBased/NBS/NBSLevel0.cpp
contactSearch/methods/cellBased/NBS/NBS.cpp
contactSearch/methods/cellBased/NBS/cellsWallLevel0.cpp
contactSearch/contactSearch/contactSearch.cpp
contactSearch/ContactSearch/ContactSearchs.cpp
interaction/interaction.cpp

15
src/Interaction/Models/contactForce/linearCF.hpp
View File

@ -133,8 +133,9 @@ protected:
return false;
}
realVector etha_n(nElem);
realVector etha_t(nElem);
realVector etha_n("etha_n", nElem);
realVector etha_t("etha_t", nElem);
ForAll(i , kn)
{
etha_n[i] = -2.0*log(en[i])*sqrt(kn[i])/
@ -144,7 +145,7 @@ protected:
sqrt(pow(log(et[i]),2.0)+ pow(Pi,2.0));
}
Vector<linearProperties> prop(nElem);
Vector<linearProperties> prop("prop", nElem);
ForAll(i,kn)
{
prop[i] = {kn[i], kt[i], etha_n[i], etha_t[i], mu[i]};
@ -219,10 +220,10 @@ public:
void contactForce
(
const real dt,
const int32 i,
const int32 j,
const int32 propId_i,
const int32 propId_j,
const uint32 i,
const uint32 j,
const uint32 propId_i,
const uint32 propId_j,
const real Ri,
const real Rj,
const real ovrlp_n,

12
src/Interaction/Models/contactForce/nonLinearCF.hpp
View File

@ -121,7 +121,7 @@ protected:
}
realVector etha_n(nElem);
realVector etha_n("etha_n",nElem);
ForAll(i , en)
{
@ -137,7 +137,7 @@ protected:
}
Vector<nonLinearProperties> prop(nElem);
Vector<nonLinearProperties> prop("prop",nElem);
ForAll(i,Yeff)
{
prop[i] = {Yeff[i], Geff[i], etha_n[i], mu[i]};
@ -214,10 +214,10 @@ public:
void contactForce
(
const real dt,
const int32 i,
const int32 j,
const int32 propId_i,
const int32 propId_j,
const uint32 i,
const uint32 j,
const uint32 propId_i,
const uint32 propId_j,
const real Ri,
const real Rj,
const real ovrlp_n,

10
src/Interaction/Models/contactForce/nonLinearMod.hpp
View File

@ -121,7 +121,7 @@ protected:
}
Vector<nonLinearProperties> prop(nElem);
Vector<nonLinearProperties> prop("prop",nElem);
ForAll(i,Yeff)
{
prop[i] = {Yeff[i], Geff[i], etha_n[i], mu[i]};
@ -198,10 +198,10 @@ public:
void contactForce
(
const real dt,
const int32 i,
const int32 j,
const int32 propId_i,
const int32 propId_j,
const uint32 i,
const uint32 j,
const uint32 propId_i,
const uint32 propId_j,
const real Ri,
const real Rj,
const real ovrlp_n,

8
src/Interaction/Models/rolling/normalRolling.hpp
View File

@ -76,10 +76,10 @@ public:
void rollingFriction
(
const real dt,
const int32 i,
const int32 j,
const int32 propId_i,
const int32 propId_j,
const uint32 i,
const uint32 j,
const uint32 propId_i,
const uint32 propId_j,
const real Ri,
const real Rj,
const realx3& wi,

22
src/Interaction/contactLists/sortedContactList.hpp
View File

@ -53,7 +53,7 @@ protected:
ViewType1D<ValueType,ExecutionSpace> values_;
int32 size0_ = 0;
uint32 size0_ = 0;
ViewType1D<PairType,ExecutionSpace> sortedPairs0_;
@ -73,7 +73,7 @@ protected:
using rpReFillPairs = Kokkos::RangePolicy<
ExecutionSpace,
Kokkos::Schedule<Kokkos::Static>,
Kokkos::IndexType<int32>,
Kokkos::IndexType<uint32>,
TagReFillPairs>;
public:
@ -81,7 +81,7 @@ public:
TypeInfoNV("sortedContactList");
sortedContactList(int32 initialSize =1)
explicit sortedContactList(uint32 initialSize =1)
:
SortedPairs(initialSize),
values_("values", SortedPairs::capacity()),
@ -114,31 +114,31 @@ public:
}
INLINE_FUNCTION_HD
ValueType getValue(int32 idx)const
ValueType getValue(uint32 idx)const
{
return values_[idx];
}
INLINE_FUNCTION_HD
void setValue(int32 idx, const ValueType& val)const
void setValue(uint32 idx, const ValueType& val)const
{
values_[idx] = val;
}
INLINE_FUNCTION_HD
void operator()(TagReFillPairs, int32 idx)const
void operator()(TagReFillPairs, uint32 idx)const
{
auto searchLen = max(size0_/1000,10);
auto start = max(0,idx-searchLen);
auto end = min(size0_,idx+searchLen);
uint32 searchLen = max(size0_/1000u,10u);
uint32 start = idx-min(searchLen,idx);
uint32 end = min(size0_,idx+searchLen);
auto newPair = this->sortedPairs_[idx];
if( auto idx0 = binarySearch(
sortedPairs0_,
start,
end,
newPair);
idx0>=0)
idx0!=-1)
{
values_[idx] = values0_[idx0];
}
@ -147,7 +147,7 @@ public:
start,
end,
newPair);
idx0>=0)
idx0!=-1)
{
values_[idx] = values0_[idx0];

40
src/Interaction/contactLists/sortedPairs.hpp
View File

@ -52,24 +52,24 @@ public:
{
using PairType = typename sortedPairs::PairType;
int32 size_;
uint32 size_;
ViewType1D<PairType,ExecutionSpace> sortedParis_;
INLINE_FUNCTION_HD
int32 size()const { return size_; }
uint32 size()const { return size_; }
INLINE_FUNCTION_HD
int32 loopCount()const { return size_; }
uint32 loopCount()const { return size_; }
INLINE_FUNCTION_HD
bool isValid(int32 i)const { return i<size_; }
bool isValid(uint32 i)const { return i<size_; }
INLINE_FUNCTION_HD
PairType getPair(int i)const { return sortedParis_[i]; }
PairType getPair(uint32 i)const { return sortedParis_[i]; }
INLINE_FUNCTION_HD
bool getPair(int32 i, PairType& pair)const {
bool getPair(uint32 i, PairType& pair)const {
if(i<size_) {
pair = sortedParis_[i];
return true;
@ -85,22 +85,22 @@ public:
protected:
/// size of pair list
int32 size_ = 0;
uint32 size_ = 0;
ViewType1D<int32,ExecutionSpace> flags_;
ViewType1D<uint32,ExecutionSpace> flags_;
ViewType1D<PairType,ExecutionSpace> sortedPairs_;
using rpFillFlag = Kokkos::RangePolicy<
ExecutionSpace,
Kokkos::Schedule<Kokkos::Static>,
Kokkos::IndexType<int32>,
Kokkos::IndexType<uint32>,
TagFillFlag >;
using rpFillPairs = Kokkos::RangePolicy<
ExecutionSpace,
Kokkos::Schedule<Kokkos::Static>,
Kokkos::IndexType<int32>,
Kokkos::IndexType<uint32>,
TagFillPairs>;
public:
@ -110,7 +110,7 @@ public:
// constructors
sortedPairs(int32 initialSize =1)
explicit sortedPairs(uint32 initialSize =1)
:
UnsortedPairs(initialSize),
flags_("flags_",UnsortedPairs::capacity()+1),
@ -134,7 +134,7 @@ public:
// return the pair at index idx
// perform no check for size and existance
INLINE_FUNCTION_HD
PairType getPair(int32 idx)const
PairType getPair(uint32 idx)const
{
return sortedPairs_[idx];
}
@ -142,7 +142,7 @@ public:
// - Device/host call
// return the pair at index idx
INLINE_FUNCTION_HD
bool getPair(int32 idx, PairType& p)const
bool getPair(uint32 idx, PairType& p)const
{
if(isValid(idx))
{
@ -156,7 +156,7 @@ public:
}
INLINE_FUNCTION_HD
bool isValid(int32 idx)const
bool isValid(uint32 idx)const
{
return idx < size_;
}
@ -164,12 +164,12 @@ public:
//use this when the value of size_ is updated
INLINE_FUNCTION_H
int32 size()const
uint32 size()const
{
return size_;
}
int32 loopCount()const
uint32 loopCount()const
{
return size_;
}
@ -189,7 +189,7 @@ public:
void prepareSorted()
{
// first check the size of flags_
int32 capacity = UnsortedPairs::capacity();
uint32 capacity = UnsortedPairs::capacity();
if( capacity+1 > flags_.size() )
{
@ -218,7 +218,7 @@ public:
if( size_>sortedPairs_.size() )
{
// get more space to prevent reallocations in next iterations
int32 len = size_*1.1+1;
uint32 len = size_*1.1+1;
reallocNoInit(sortedPairs_, len);
}
@ -235,7 +235,7 @@ public:
}
INLINE_FUNCTION_HD
void operator()(TagFillFlag, int32 i)const
void operator()(TagFillFlag, uint32 i)const
{
if(this->container_.valid_at(i) )
flags_[i] = 1;
@ -244,7 +244,7 @@ public:
}
INLINE_FUNCTION_HD
void operator()(TagFillPairs, int32 i)const
void operator()(TagFillPairs, uint32 i)const
{
auto fi = flags_[i];
if(fi!=flags_[i+1])

16
src/Interaction/contactLists/unsortedContactList.hpp
View File

@ -72,7 +72,7 @@ protected:
using rpFillPairs = Kokkos::RangePolicy<
ExecutionSpace,
Kokkos::Schedule<Kokkos::Static>,
Kokkos::IndexType<int32>,
Kokkos::IndexType<uint32>,
TagReFillPairs>;
@ -80,7 +80,7 @@ public:
TypeInfoNV("unsortedContactList");
unsortedContactList(int32 capacity=1)
explicit unsortedContactList(uint32 capacity=1)
:
UnsortedPairs(capacity),
values_("values", UnsortedPairs::capacity()),
@ -122,7 +122,7 @@ public:
INLINE_FUNCTION_HD
bool getValue(const PairType& p, ValueType& val)const
{
if(auto idx = this->find(p); idx>=0)
if(auto idx = this->find(p); idx!=-1)
{
val = getValue(idx);
return true;
@ -131,7 +131,7 @@ public:
}
INLINE_FUNCTION_HD
void setValue(int32 idx, const ValueType& val)const
void setValue(uint32 idx, const ValueType& val)const
{
values_[idx] = val;
}
@ -139,7 +139,7 @@ public:
INLINE_FUNCTION_HD
bool setValue(const PairType& p, const ValueType& val)const
{
if(auto idx = this->find(p); idx>=0)
if(uint32 idx = this->find(p); idx!=-1)
{
setValue(idx, val);
return true;;
@ -148,13 +148,13 @@ public:
}
INLINE_FUNCTION_HD
void operator()(TagReFillPairs, int32 idx)const
void operator()(TagReFillPairs, uint32 idx)const
{
if( this->isValid(idx) )
{
if( int32 idx0 =
if( uint32 idx0 =
container0_.find(this->getPair(idx));
idx0>=0 )
idx0!=-1 )
{
values_[idx] = values0_[idx0];
}

34
src/Interaction/contactLists/unsortedPairs.hpp
View File

@ -41,7 +41,7 @@ public:
using memory_space = typename ExecutionSpace::memory_space;
using PairType = kPair<idType,idType>;
using PairType = Pair<idType,idType>;
using ContainerType = unorderedSet<PairType, ExecutionSpace>;
@ -52,19 +52,19 @@ public:
ContainerType Container_;
INLINE_FUNCTION_HD
int32 size()const { return Container_.size(); }
uint32 size()const { return Container_.size(); }
INLINE_FUNCTION_HD
int32 loopCount()const { return Container_.capacity(); }
uint32 loopCount()const { return Container_.capacity(); }
INLINE_FUNCTION_HD
bool isValid(int32 idx)const { return Container_.valid_at(idx); }
bool isValid(uint32 idx)const { return Container_.valid_at(idx); }
INLINE_FUNCTION_HD
PairType getPair(int idx)const { return Container_.key_at(idx); }
PairType getPair(uint32 idx)const { return Container_.key_at(idx); }
INLINE_FUNCTION_HD
bool getPair(int32 idx, PairType& pair)const {
bool getPair(uint32 idx, PairType& pair)const {
if(Container_.valid_at(idx)) {
pair = Container_.key_at(idx);
return true;
@ -84,7 +84,7 @@ public:
TypeInfoNV("unsorderedPairs");
// constructor
unsortedPairs(int32 capacity=1)
explicit unsortedPairs(uint32 capacity=1)
:
container_(capacity) // the minimum capacity would be 128
{}
@ -102,7 +102,7 @@ public:
// - Device call
INLINE_FUNCTION_HD
int32 insert(idType i, idType j)const
uint32 insert(idType i, idType j)const
{
if(auto insertResult = container_.insert(PairType(i,j)); insertResult.failed())
return -1;
@ -112,7 +112,7 @@ public:
}
INLINE_FUNCTION_HD
int32 insert(const PairType& p)const
uint32 insert(const PairType& p)const
{
if(auto insertResult = container_.insert(p); insertResult.failed())
return -1;
@ -125,7 +125,7 @@ public:
// return the pair at index idx
// perform no check for size and existance
INLINE_FUNCTION_HD
PairType getPair(int32 idx)const
PairType getPair(uint32 idx)const
{
return container_.key_at(idx);
}
@ -133,7 +133,7 @@ public:
// - Device call
// return the pair at index idx
INLINE_FUNCTION_HD
bool getPair(int32 idx, PairType& p)const
bool getPair(uint32 idx, PairType& p)const
{
if(container_.valid_at(idx))
{
@ -148,7 +148,7 @@ public:
}
INLINE_FUNCTION_HD
int32 find(const PairType & p)const
uint32 find(const PairType & p)const
{
if( auto idx = container_.find(p);
idx != Kokkos::UnorderedMapInvalidIndex )
@ -158,26 +158,26 @@ public:
}
INLINE_FUNCTION_HD
bool isValid(int32 idx)const
bool isValid(uint32 idx)const
{
return container_.valid_at(idx);
}
INLINE_FUNCTION_HD
int32 capacity() const
uint32 capacity() const
{
return container_.capacity();
}
int32 loopCount()const
uint32 loopCount()const
{
return container_.capacity();
}
//use this when the value of size_ is updated
INLINE_FUNCTION_H
int32 size()const
uint32 size()const
{
return container_.size();
}
@ -190,7 +190,7 @@ public:
/// increase the capacity of the container by at-least len
/// the content will be erased.
INLINE_FUNCTION_H
void increaseCapacityBy(int32 len)
void increaseCapacityBy(uint32 len)
{
uint newCap = container_.capacity()+len;
this->clear();

199
src/Interaction/contactSearch/ContactSearch/ContactSearch.hpp
View File

@ -22,16 +22,17 @@ Licence:
#ifndef __ContactSearch_hpp__
#define __ContactSearch_hpp__
#include "contactSearchGlobals.hpp"
#include "contactSearch.hpp"
#include "box.hpp"
#include "particles.hpp"
#include "geometry.hpp"
namespace pFlow
{
template<
template<class> class BaseMethod,
template<class> class WallMapping
class searchMethod
>
class ContactSearch
:
@ -39,32 +40,19 @@ class ContactSearch
{
public:
using IdType = typename contactSearch::IdType;
using IdType = uint32;
using IndexType = typename contactSearch::IndexType;
using ExecutionSpace = DefaultExecutionSpace;
using ExecutionSpace = typename contactSearch::ExecutionSpace;
using SearchMethodType = searchMethod;
using PairContainerType = typename contactSearch::PairContainerType;
using ParticleContactSearchType =
BaseMethod<
ExecutionSpace>;
using WallMappingType =
WallMapping<
ExecutionSpace>;
protected:
private:
uniquePtr<ParticleContactSearchType> particleContactSearch_ = nullptr;
uniquePtr<WallMappingType> wallMapping_ = nullptr;
uniquePtr<SearchMethodType> ppwContactSearch_ = nullptr;
public:
TypeInfoTemplate2("ContactSearch", ParticleContactSearchType, WallMappingType);
TypeInfoTemplate11("ContactSearch", SearchMethodType);
ContactSearch(
const dictionary& csDict,
@ -74,54 +62,41 @@ public:
Timers& timers)
:
contactSearch(csDict, domain, prtcl, geom, timers)
{
auto method = dict().getVal<word>("method");
auto wmMethod = dict().getVal<word>("wallMapping");
auto nbDict = dict().subDict(method+"Info");
real minD, maxD;
this->Particles().boundingSphereMinMax(minD, maxD);
const auto& position = this->Particles().pointPosition().deviceVectorAll();
const auto& diam = this->Particles().boundingSphere().deviceVectorAll();
const auto& position = this->Particles().pointPosition().deviceViewAll();
const auto& flags = this->Particles().dynPointStruct().activePointsMaskDevice();
const auto& diam = this->Particles().boundingSphere().deviceViewAll();
particleContactSearch_ =
makeUnique<ParticleContactSearchType>
uint32 wnPoints = this->Geometry().numPoints();
uint32 wnTri = this->Geometry().size();
const auto& wPoints = this->Geometry().points().deviceViewAll();
const auto& wVertices = this->Geometry().vertices().deviceViewAll();
ppwContactSearch_ =
makeUnique<SearchMethodType>
(
nbDict,
this->domain(),
this->domainBox(),
minD,
maxD,
position,
diam
);
REPORT(2)<<"Contact search algorithm for particle-particle is "<<
greenText(particleContactSearch_().typeName())<<endREPORT;
auto wmDict = dict().subDict(wmMethod+"Info");
int32 wnPoints = this->Geometry().numPoints();
int32 wnTri = this->Geometry().size();
const auto& wPoints = this->Geometry().points().deviceVectorAll();
const auto& wVertices = this->Geometry().vertices().deviceVectorAll();
wallMapping_ =
makeUnique<WallMappingType>(
wmDict,
particleContactSearch_().numLevels(),
particleContactSearch_().getCellsLevels(),
flags,
diam,
wnPoints,
wnTri,
wnTri,
wPoints,
wVertices
);
REPORT(2)<<"Wall mapping algorithm for particle-wall is "<<
greenText(wallMapping_().typeName())<< endREPORT;
);
REPORT(2)<<"Contact search algorithm for particle-particle is "<<
Green_Text(ppwContactSearch_().typeName())<<END_REPORT;
}
@ -130,116 +105,60 @@ public:
contactSearch,
ContactSearch,
dictionary);
bool broadSearch(
PairContainerType& ppPairs,
PairContainerType& pwPairs,
uint32 iter,
real t,
real dt,
csPairContainerType& ppPairs,
csPairContainerType& pwPairs,
bool force = false) override
{
ppTimer().start();
const auto& position = Particles().pointPosition().deviceViewAll();
const auto& flags = Particles().dynPointStruct().activePointsMaskDevice();
const auto& diam = Particles().boundingSphere().deviceViewAll();
if(particleContactSearch_)
if( !ppwContactSearch_().broadSearch(
iter,
t,
dt,
ppPairs,
pwPairs,
position,
flags,
diam,
force) )
{
auto activeRange = this->Particles().activeRange();
sphereSphereTimer_.start();
if(this->Particles().allActive())
{
particleContactSearch_().broadSearch(ppPairs, activeRange, force);
}
else
{
particleContactSearch_().broadSearch(ppPairs, activeRange, this->Particles().activePointsMaskD(), force);
}
sphereSphereTimer_.end();
}
else
fatalErrorInFunction;
return false;
if(wallMapping_)
{
sphereWallTimer_.start();
wallMapping_().broadSearch(pwPairs, particleContactSearch_(), force);
sphereWallTimer_.end();
}
else
return false;
ppTimer().end();
return true;
}
bool ppEnterBroadSearch()const override
bool enterBroadSearch(uint32 iter, real t, real dt)const override
{
if(particleContactSearch_)
if(ppwContactSearch_)
{
return particleContactSearch_().enterBoadSearch();
return ppwContactSearch_().performSearch(iter);
}
return false;
}
bool pwEnterBroadSearch()const override
bool performedBroadSearch()const override
{
if(wallMapping_)
{
return wallMapping_().enterBoadSearch();
}
return false;
return ppwContactSearch_().performedSearch();
}
bool ppPerformedBroadSearch()const override
{
if(particleContactSearch_)
{
return particleContactSearch_().performedSearch();
}
return false;
}
bool pwPerformedBroadSearch()const override
{
if(wallMapping_)
{
return wallMapping_().performedSearch();
}
return false;
}
/*bool update(const eventMessage& msg)
{
if(msg.isSizeChanged() )
{
auto newSize = this->prtcl().size();
if(!particleContactSearch_().objectSizeChanged(newSize))
{
fatalErrorInFunction<<
"erro in changing the size for particleContactSearch_ \n";
return false;
}
}
if(msg.isCapacityChanged() )
{
auto newSize = this->prtcl().capacity();
if(!particleContactSearch_().objectSizeChanged(newSize))
{
fatalErrorInFunction<<
"erro in changing the capacity for particleContactSearch_ \n";
return false;
}
}
return true;
}*/
};
}

10
src/Interaction/contactSearch/ContactSearch/ContactSearchs.cpp
View File

@ -20,11 +20,11 @@ Licence:
#include "ContactSearch.hpp"
#include "cellMapping.hpp"
//#include "cellMapping.hpp"
#include "NBS.hpp"
#include "multiGridNBS.hpp"
#include "multiGridMapping.hpp"
//#include "multiGridNBS.hpp"
//#include "multiGridMapping.hpp"
template class pFlow::ContactSearch<pFlow::NBS, pFlow::cellMapping>;
template class pFlow::ContactSearch<pFlow::multiGridNBS, pFlow::multiGridMapping>;
template class pFlow::ContactSearch<pFlow::NBS>;
//template class pFlow::ContactSearch<pFlow::multiGridNBS, pFlow::multiGridMapping>;

25
src/Interaction/contactSearch/contactSearch/contactSearch.cpp
View File

@ -19,7 +19,7 @@ Licence:
-----------------------------------------------------------------------------*/
#include "contactSearch.hpp"
#include "streams.hpp"
pFlow::contactSearch::contactSearch(
@ -29,11 +29,12 @@ pFlow::contactSearch::contactSearch(
const geometry& geom,
Timers& timers)
:
interactionBase(prtcl, geom),
domain_(domain),
dict_(dict),
sphereSphereTimer_("particle-particle contact search", &timers),
sphereWallTimer_("particle-wall contact search", &timers)
domainBox_(domain),
particles_(prtcl),
geometry_(geom),
ppTimer_("particle-particle contact search", &timers),
pwTimer_("particle-wall contact search", &timers),
dict_(dict)
{
}
@ -46,20 +47,16 @@ pFlow::uniquePtr<pFlow::contactSearch> pFlow::contactSearch::create(
const geometry& geom,
Timers& timers)
{
word baseMethName = dict.getVal<word>("method");
word baseMethName = dict.getVal<word>("method");
word wallMethod = dict.getVal<word>("wallMapping");
auto model = angleBracketsNames("ContactSearch", baseMethName);
auto model = angleBracketsNames2("ContactSearch", baseMethName, wallMethod);
REPORT(1)<<"Selecting contact search model . . ."<<endREPORT;
REPORT(1)<<"Selecting contact search model . . ."<<END_REPORT;
if( dictionaryvCtorSelector_.search(model))
{
auto objPtr = dictionaryvCtorSelector_[model] (dict, domain, prtcl, geom, timers);
REPORT(2)<<"Model "<< greenText(model)<<" is created."<<endREPORT;
REPORT(2)<<"Model "<< Green_Text(model)<<" is created."<< END_REPORT;
return objPtr;
}
else

83
src/Interaction/contactSearch/contactSearch/contactSearch.hpp
View File

@ -23,43 +23,36 @@ Licence:
#define __contactSearch_hpp__
#include "interactionBase.hpp"
#include "unsortedPairs.hpp"
#include "box.hpp"
#include "contactSearchGlobals.hpp"
#include "dictionary.hpp"
#include "virtualConstructor.hpp"
#include "Timer.hpp"
namespace pFlow
{
// - forward
class box;
class particles;
class geometry;
class contactSearch
:
public interactionBase
{
public:
using IdType = typename interactionBase::IdType;
private:
using IndexType = typename interactionBase::IndexType;
const box& domainBox_;
using ExecutionSpace = typename interactionBase::ExecutionSpace;
const particles& particles_;
using PairContainerType = unsortedPairs<ExecutionSpace, IdType>;
const geometry& geometry_;
protected:
Timer ppTimer_;
const box& domain_;
Timer pwTimer_;
dictionary dict_;
Timer sphereSphereTimer_;
Timer sphereWallTimer_;
auto& dict()
{
return dict_;
}
public:
TypeInfo("contactSearch");
@ -88,35 +81,51 @@ public:
(dict, domain, prtcl, geom, timers)
);
const auto& domain()const
{
return domain_;
}
const auto& dict()const
{
return dict_;
}
const auto& domainBox()const
{
return domainBox_;
}
const auto& Particles()const
{
return particles_;
}
const auto& Geometry()const
{
return geometry_;
}
auto& ppTimer()
{
return ppTimer_;
}
auto& pwTimer()
{
return pwTimer_;
}
virtual
bool broadSearch(
PairContainerType& ppPairs,
PairContainerType& pwPairs,
uint32 iter,
real t,
real dt,
csPairContainerType& ppPairs,
csPairContainerType& pwPairs,
bool force = false) = 0;
virtual
bool ppEnterBroadSearch()const = 0;
bool enterBroadSearch(uint32 iter, real t, real dt)const = 0;
virtual
bool pwEnterBroadSearch()const = 0;
virtual
bool ppPerformedBroadSearch()const = 0;
virtual
bool pwPerformedBroadSearch()const = 0;
bool performedBroadSearch()const = 0;
static
uniquePtr<contactSearch> create(

11
src/Interaction/contactSearch/contactSearchFunctions.hpp
View File

@ -96,17 +96,6 @@ void indexToCell(const indexType idx, const iBox<cellIndexType>& box, triple<cel
cell+= box.minPoint();
}
INLINE_FUNCTION_HD
bool sphereSphereCheck(const realx3& p1, const realx3 p2, real d1, real d2)
{
return length(p2-p1) < 0.5*(d2+d1);
}
}
#endif //__broadSearchFunctions_hpp__

23
src/phasicFlow/structuredData/trisurfaceStructure/triangleFunctions.hpp → src/Interaction/contactSearch/contactSearchGlobals.hpp
View File

@ -18,22 +18,21 @@ Licence:
-----------------------------------------------------------------------------*/
#ifndef __triangleFunctions_hpp__
#define __triangleFunctions_hpp__
#include "types.hpp"
#include "unsortedPairs.hpp"
namespace pFlow::triangleFunctions
#ifndef __contactSearchGlobals_hpp__
#define __contactSearchGlobals_hpp__
namespace pFlow
{
INLINE_FUNCTION_HD
real triangleSurface( const realx3& p1, const realx3& p2, const realx3& p3)
{
realx3 V1 = p2 - p1;
realx3 V2 = p3 - p1;
return abs((cross(V1,V2)).length() / static_cast<real>(2.0));
}
using csExecutionSpace = DefaultExecutionSpace;
using csIdType = uint32;
using csPairContainerType = unsortedPairs<DefaultExecutionSpace, uint32>;
}
#endif
#endif

209
src/Interaction/contactSearch/methods/NBS.hpp
View File

@ -1,209 +0,0 @@
/*------------------------------- 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 __NBS_hpp__
#define __NBS_hpp__
#include "NBSLevel0.hpp"
namespace pFlow
{
template<typename executionSpace>
class NBS
{
public:
using NBSLevel0Type = NBSLevel0<executionSpace>;
using cellIterator = typename NBSLevel0Type::cellIterator;
using IdType = typename NBSLevel0Type::IdType;
using IndexType = typename NBSLevel0Type::IndexType;
using Cells = typename NBSLevel0Type::Cells;
using CellType = typename Cells::CellType;
using execution_space = typename NBSLevel0Type::execution_space;
using memory_space = typename NBSLevel0Type::memory_space;
protected:
real sizeRatio_ = 1.0;
int32 updateFrequency_= 1;
int32 currentIter_ = 0;
bool performedSearch_ = false;
NBSLevel0Type NBSLevel0_;
private:
bool performSearch()
{
if(currentIter_ % updateFrequency_ == 0)
{
currentIter_++;
return true;
}else
{
currentIter_++;
return false;
}
}
public:
TypeInfoNV("NBS");
NBS(
const dictionary& dict,
const box& domain,
real minSize,
real maxSize,
const ViewType1D<realx3, memory_space>& position,
const ViewType1D<real, memory_space>& diam)
:
sizeRatio_(
max(
dict.getValOrSet<real>("sizeRatio", 1.0),
1.0
)),
updateFrequency_(
max(
dict.getValOrSet<int32>("updateFrequency", 1),
1
)),
NBSLevel0_(
domain,
maxSize*sizeRatio_,
sizeRatio_,
position,
diam)
{}
INLINE_FUNCTION_HD
NBS(const NBS&) = default;
INLINE_FUNCTION_HD
NBS& operator = (const NBS&) = default;
INLINE_FUNCTION_HD
~NBS()=default;
//// - Methods
bool enterBoadSearch()const
{
return currentIter_%updateFrequency_==0;
}
bool performedSearch()const
{
return performedSearch_;
}
Vector<cellIterator> getCellIteratorLevels()
{
return Vector<cellIterator>("cellIterator", 1, NBSLevel0_.getCellIterator());
}
auto getCellIterator(int32 lvl)const
{
return NBSLevel0_.getCellIterator();
}
int32 numLevels()const
{
return 1;
}
Vector<Cells> getCellsLevels()const
{
return Vector<Cells>("Cells", 1, NBSLevel0_.getCells());
}
auto getCells()const
{
return NBSLevel0_.getCells();
}
bool objectSizeChanged(int32 newSize)
{
NBSLevel0_.checkAllocateNext(newSize);
return true;
}
// - Perform the broad search to find pairs
// with force = true, perform broad search regardless of
// updateFrequency_ value
// on all the points in the range of [0,numPoints_)
template<typename PairsContainer>
bool broadSearch(PairsContainer& pairs, range activeRange, bool force=false)
{
if(force) currentIter_ = 0;
performedSearch_ = false;
if( !performSearch() ) return true;
NBSLevel0_.build(activeRange);
NBSLevel0_.findPairs(pairs);
performedSearch_ = true;
return true;
}
// - Perform the broad search to find pairs,
// ignore particles with incld(i) = true,
// with force = true, perform broad search regardless of
// updateFrequency_ value
template<typename PairsContainer, typename IncludeFunction>
bool broadSearch(PairsContainer& pairs, range activeRange, IncludeFunction incld, bool force = false)
{
if(force) currentIter_ = 0;
performedSearch_ = false;
if( !performSearch() ) return true;
NBSLevel0_.build(activeRange, incld);
NBSLevel0_.findPairs(pairs);
performedSearch_ = true;
return true;
}
};
}
#endif

82
src/Interaction/contactSearch/methods/NBSCrossLoop.hpp
View File

@ -1,82 +0,0 @@
/*------------------------------- 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.
-----------------------------------------------------------------------------*/
int32 m = this->head_(i,j,k);
CellType currentCell(i,j,k);
int32 n = -1;
while( m > -1 )
{
auto p_m = this->pointPosition_[m];
auto d_m = this->sizeRatio_* this->diameter_[m];
int32x3 crossIndex = mapIndexLevels(
int32x3(i,j,k),
level_,
upperLevel.level());
for(int32 ii = -1; ii<2; ii++)
{
for(int32 jj=-1; jj<2; jj++)
{
int32 kk=-1;
while( kk < 2)
{
int32x3 nghbrCI = crossIndex + int32x3(ii,jj,kk);
if( upperLevel.isInRange(nghbrCI) )
{
n = upperLevel.head_(
nghbrCI.x(),
nghbrCI.y(),
nghbrCI.z());
while( n >-1)
{
auto p_n = this->pointPosition_[n];
auto d_n = this->sizeRatio_*this->diameter_[n];
if( sphereSphereCheck(p_m, p_n, d_m, d_n) )
{
auto ln = n;
auto lm = m;
if(lm>ln) this->Swap(lm,ln);
if( auto res = pairs.insert(lm,ln); res <0)
{
getFullUpdate++;
}
}
n = this->next_[n];
}
}
kk++;
}
}
}
m = this->next_[m];
}

127
src/Interaction/contactSearch/methods/NBSLevel.hpp
View File

@ -1,127 +0,0 @@
#ifndef __NBSLevel_hpp__
#define __NBSLevel_hpp__
#include "NBSLevel0.hpp"
namespace pFlow
{
INLINE_FUNCTION_HD
int32x3 mapIndexLevels(const int32x3& ind, int32 lowerLevel, int32 upperLevel);
template<typename executionSpace>
class
NBSLevel
:
public NBSLevel0<executionSpace>
{
public:
using NBSLevel0Type = NBSLevel0<executionSpace>;
using cellIterator = typename NBSLevel0Type::cellIterator;
using IdType = typename NBSLevel0Type::IdType;
using IndexType = typename NBSLevel0Type::IndexType;
using Cells = typename NBSLevel0Type::Cells;
using CellType = typename Cells::CellType;
using execution_space = typename NBSLevel0Type::execution_space;
using memory_space = typename NBSLevel0Type::memory_space;
using mdrPolicyFindPairs = typename NBSLevel0Type::mdrPolicyFindPairs;
using HeadType = typename NBSLevel0Type::HeadType;
using NextType = typename NBSLevel0Type::NextType;
template<typename exeSpace>
friend class NBSLevels;
protected:
int32 level_ = 0;
public:
TypeInfoNV("NBSLevel0");
INLINE_FUNCTION_HD
NBSLevel(){}
NBSLevel(
int32 lvl,
const box& domain,
real cellSize,
real sizeRatio,
const ViewType1D<realx3, memory_space>& position,
const ViewType1D<real, memory_space>& diam)
:
NBSLevel0Type(
domain,
cellSize,
sizeRatio,
position,
diam,
lvl==0),
level_(lvl)
{}
INLINE_FUNCTION_HD
NBSLevel(const NBSLevel&) = default;
INLINE_FUNCTION_HD
NBSLevel& operator = (const NBSLevel&) = default;
INLINE_FUNCTION_HD
~NBSLevel() = default;
INLINE_FUNCTION_HD
auto level()const
{
return level_;
}
template<typename PairsContainer>
INLINE_FUNCTION_H
int32 findPairsCountCross(PairsContainer& pairs, NBSLevel& upperLevel)
{
mdrPolicyFindPairs
mdrPolicy(
{0,0,0},
{this->nx(),this->ny(),this->nz()} );
int32 notInsertedPairs;
Kokkos::parallel_reduce (
"NBSLevel::findPairsCountCross",
mdrPolicy,
CLASS_LAMBDA_HD(int32 i, int32 j, int32 k, int32& getFullUpdate){
#include "NBSCrossLoop.hpp"
}, notInsertedPairs);
return notInsertedPairs;
}
};
INLINE_FUNCTION_HD
int32x3 mapIndexLevels( const int32x3& ind, int32 lowerLevel, int32 upperLevel)
{
int32 a = pow(2, static_cast<int32>(upperLevel-lowerLevel));
return ind/a;
}
}
#endif

240
src/Interaction/contactSearch/methods/NBSLevel0.hpp
View File

@ -1,240 +0,0 @@
/*------------------------------- 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 __NBSLevel0_hpp__
#define __NBSLevel0_hpp__
#include "mapperNBS.hpp"
namespace pFlow
{
template<typename executionSpace>
class NBSLevel0
:
public mapperNBS<executionSpace>
{
public:
using MapperType = mapperNBS<executionSpace>;
using cellIterator = typename MapperType::cellIterator;
using IdType = typename MapperType::IdType;
using IndexType = typename MapperType::IndexType;
using Cells = typename MapperType::Cells;
using CellType = typename Cells::CellType;
using execution_space = typename MapperType::execution_space;
using memory_space = typename MapperType::memory_space;
using HeadType = typename MapperType::HeadType;
using NextType = typename MapperType::NextType;
struct TagFindPairs{};
protected:
real sizeRatio_ = 1.0;
// borrowed ownership
ViewType1D<real, memory_space> diameter_;
using mdrPolicyFindPairs =
Kokkos::MDRangePolicy<
Kokkos::Rank<3>,
Kokkos::Schedule<Kokkos::Dynamic>,
execution_space>;
static INLINE_FUNCTION_HD
void Swap(int32& x, int32& y)
{
int32 tmp = x;
x = y;
y = tmp;
}
public:
TypeInfoNV("NBSLevel0");
INLINE_FUNCTION_HD
NBSLevel0(){}
NBSLevel0(
const box& domain,
real cellSize,
const ViewType1D<realx3, memory_space>& position,
const ViewType1D<real, memory_space>& diam)
:
MapperType(domain, cellSize, position),
diameter_(diam)
{}
NBSLevel0(
const box& domain,
int32 nx,
int32 ny,
int32 nz,
const ViewType1D<realx3, memory_space>& position,
const ViewType1D<real, memory_space>& diam)
:
MapperType(domain, nx, ny, nz, position),
diameter_(diam)
{ }
NBSLevel0(
const box& domain,
real cellSize,
real sizeRatio,
const ViewType1D<realx3, memory_space>& position,
const ViewType1D<real, memory_space>& diam,
bool nextOwner = true)
:
MapperType(domain, cellSize, position, nextOwner),
sizeRatio_(sizeRatio),
diameter_(diam)
{}
INLINE_FUNCTION_HD
NBSLevel0(const NBSLevel0&) = default;
INLINE_FUNCTION_HD
NBSLevel0& operator = (const NBSLevel0&) = default;
INLINE_FUNCTION_HD
~NBSLevel0()=default;
INLINE_FUNCTION_HD
auto sizeRatio()const
{
return sizeRatio_;
}
INLINE_FUNCTION_HD
auto& diameter()
{
return diameter_;
}
// - Perform the broad search to find pairs
// with force = true, perform broad search regardless of
// updateFrequency_ value
// on all the points in the range of [0,numPoints_)
template<typename PairsContainer>
bool broadSearch(PairsContainer& pairs, range activeRange)
{
this->build(activeRange);
findPairs(pairs);
return true;
}
// - Perform the broad search to find pairs,
// ignore particles with incld(i) = true,
// with force = true, perform broad search regardless of
// updateFrequency_ value
template<typename PairsContainer, typename IncludeFunction>
bool broadSearch(PairsContainer& pairs, range activeRange, IncludeFunction incld)
{
this->build(activeRange, incld);
findPairs(pairs);
return true;
}
template<typename PairsContainer>
INLINE_FUNCTION_H
bool findPairs(PairsContainer& pairs)
{
int32 getFull = 1;
// loop until the container size fits the numebr of contact pairs
while (getFull > 0)
{
getFull = findPairsCount(pairs);
if(getFull)
{
// - resize the container
// note that getFull now shows the number of failed insertions.
uint32 len = max(getFull,500) ;
auto oldCap = pairs.capacity();
pairs.increaseCapacityBy(len);
INFORMATION<< "The contact pair container capacity increased from "<<
oldCap << " to "<<pairs.capacity()<<" in NBSLevel0."<<endINFO;
}
Kokkos::fence();
}
return true;
}
template<typename PairsContainer>
INLINE_FUNCTION_H
int32 findPairsCount(PairsContainer& pairs)
{
mdrPolicyFindPairs
mdrPolicy(
{0,0,0},
{this->nx(),this->ny(),this->nz()} );
int32 notInsertedPairs;
Kokkos::parallel_reduce (
"NBSLevel0::findPairs",
mdrPolicy,
CLASS_LAMBDA_HD(int32 i, int32 j, int32 k, int32& getFullUpdate){
#include "NBSLoop.hpp"
}, notInsertedPairs);
return notInsertedPairs;
}
};
} // pFlow
#endif // __NBSLevel0_hpp__

438
src/Interaction/contactSearch/methods/NBSLevels.hpp
View File

@ -1,438 +0,0 @@
#ifndef __NBSLevels_hpp__
#define __NBSLevels_hpp__
#include "NBSLevel.hpp"
#include "NBSLevel0.hpp"
#include "KokkosTypes.hpp"
namespace pFlow
{
template<typename executionSpace>
class NBSLevels
{
public:
using NBSLevelType = NBSLevel<executionSpace>;
using cellIterator = typename NBSLevelType::cellIterator;
using IdType = typename NBSLevelType::IdType;
using IndexType = typename NBSLevelType::IndexType;
using Cells = typename NBSLevelType::Cells;
using CellType = typename Cells::CellType;
using execution_space = typename NBSLevelType::execution_space;
using memory_space = typename NBSLevelType::memory_space;
using realRange = kPair<real,real>;
protected:
real minSize_;
real maxSize_;
int32 numLevels_=1;
Vector<NBSLevelType> nbsLevels_;
ViewType1D<realRange, memory_space> sizeRangeLevels_;
ViewType1D<realRange, HostSpace> sizeRangeLevelsHost_;
ViewType1D<real, memory_space> maxSizeLevels_;
ViewType1D<real, HostSpace> maxSizeLevelsHost_;
ViewType1D<int8, memory_space> particleLevel_;
range activeRange_{0,0};
using rangePolicyType =
Kokkos::RangePolicy<
Kokkos::IndexType<int32>,
Kokkos::Schedule<Kokkos::Static>,
execution_space>;
int32 setNumLevels()
{
int32 maxOvermin = static_cast<int32>(maxSize_/minSize_);
if (maxOvermin <=1)
return 1;
else if(maxOvermin<=3)
return 2;
else if(maxOvermin<=7)
return 3;
else if(maxOvermin<15)
return 4;
else if(maxOvermin<31)
return 5;
else if(maxOvermin<63)
return 6;
else if(maxOvermin <127)
return 7;
else
{
fatalErrorInFunction<<
"size ratio is not valid for multi-grid NBS "<< maxOvermin<<endl;
fatalExit;
}
return -1;
}
bool setDiameterRange(real sizeRatio)
{
real lvl_maxD = sizeRatio* maxSize_;
real lvl_minD = lvl_maxD/2.0;
for(int32 lvl=numLevels_-1; lvl>=0; lvl--)
{
if(lvl == 0 ) lvl_minD = 0.01*minSize_;
sizeRangeLevelsHost_[lvl] = {lvl_minD, lvl_maxD};
maxSizeLevelsHost_[lvl] = lvl_maxD;
lvl_maxD = lvl_minD;
lvl_minD /= 2.0;
}
copy(sizeRangeLevels_, sizeRangeLevelsHost_);
copy(maxSizeLevels_, maxSizeLevelsHost_);
REPORT(2)<<"Grids with "<< yellowText(numLevels_)<< " levels have been created."<<endREPORT;
for(int32 lvl=0; lvl<numLevels_; lvl++)
{
REPORT(3)<<"Cell gird No "<< yellowText(lvl)<<" with size range ("
<<sizeRangeLevelsHost_[lvl].first<<","<<sizeRangeLevelsHost_[lvl].second<<"]."<<endREPORT;
}
return true;
}
bool initLevels(
const box& domain,
real sizeRatio,
const ViewType1D<realx3, memory_space>& position,
const ViewType1D<real, memory_space>& diam)
{
for(int32 lvl = 0; lvl<numLevels_; lvl++)
{
nbsLevels_[lvl] = NBSLevelType(
lvl,
domain,
maxSizeLevelsHost_[lvl]*sizeRatio,
sizeRatio,
position,
diam );
}
auto next0 = nbsLevels_[0].next();
for(int32 lvl=1; lvl<numLevels_; lvl++)
{
nbsLevels_[lvl].setNext(next0);
}
return true;
}
void manageAllocateNext(range active)
{
activeRange_ = active;
if(activeRange_.second > nbsLevels_[0].capacity())
{
nbsLevels_[0].checkAllocateNext(activeRange_.second);
auto next0 = nbsLevels_[0].next();
for(int32 lvl=1; lvl<numLevels_; lvl++)
{
nbsLevels_[lvl].setNext(next0);
}
}
}
void nullify( range active)
{
for(int32 lvl=0; lvl<numLevels_; lvl++)
{
nbsLevels_[lvl].nullify(active);
}
}
public:
NBSLevels(
const box& domain,
real minSize,
real maxSize,
real sizeRatio,
const ViewType1D<realx3, memory_space>& position,
const ViewType1D<real, memory_space>& diam)
:
minSize_(minSize),
maxSize_(maxSize),
numLevels_(setNumLevels()),
nbsLevels_("nbsLevels", numLevels_, numLevels_, RESERVE()),
sizeRangeLevels_("sizeRangeLevels", numLevels_),
sizeRangeLevelsHost_("sizeRangeLevelsHost", numLevels_),
maxSizeLevels_("maxSizeLevels", numLevels_),
maxSizeLevelsHost_("maxSizeLevelsHost", numLevels_)
{
setDiameterRange(sizeRatio);
initLevels(domain, sizeRatio, position, diam);
}
auto getCellIterator(int32 lvl)const
{
return nbsLevels_[lvl].getCellIterator();
}
int32 numLevels()const
{
return numLevels_;
}
Cells getCells(int32 lvl)const
{
return nbsLevels_[lvl].getCells();
}
template<typename PairsContainer>
INLINE_FUNCTION_H
bool findPairs(PairsContainer& pairs)
{
int32 getFull = 1;
// loop until the container size fits the numebr of contact pairs
while (getFull > 0)
{
getFull = findPairsCount(pairs);
if(getFull)
{
// - resize the container
// note that getFull now shows the number of failed insertions.
uint32 len = max(getFull,100) ;
auto oldCap = pairs.capacity();
pairs.increaseCapacityBy(len);
INFORMATION<< "The contact pair container capacity increased from "<<
oldCap << " to "<<pairs.capacity()<<" in NBSLevels."<<endINFO;
}
Kokkos::fence();
}
return true;
}
template<typename PairsContainer>
INLINE_FUNCTION_H
int32 findPairsCount(PairsContainer& pairs)
{
int32 notInsertedCount = 0;
for(int32 lvl=0; lvl<numLevels_; lvl++)
{
// the same level
notInsertedCount+= nbsLevels_[lvl].findPairsCount(pairs);
for(int32 crsLvl = lvl+1; crsLvl<numLevels_; crsLvl++)
{
// cross levels
notInsertedCount+=
nbsLevels_[lvl].findPairsCountCross(pairs, nbsLevels_[crsLvl]);
}
}
return notInsertedCount;
}
INLINE_FUNCTION_H
void build(range activeRange)
{
// nullify next and heads
findParticleLevel(activeRange.first, activeRange.second);
manageAllocateNext(activeRange);
nullify(activeRange);
//
rangePolicyType rPolicy(activeRange.first, activeRange.second);
auto nbsLevel0 = nbsLevels_[0];
auto pointPosition = nbsLevel0.pointPosition();
auto particleLevel = particleLevel_;
auto next = nbsLevel0.next();
auto head0 = nbsLevel0.head();
typename NBSLevelType::HeadType head1, head2, head3, head4, head5, head6;
if(numLevels_>1) head1 = nbsLevels_[1].head();
if(numLevels_>2) head2 = nbsLevels_[2].head();
if(numLevels_>3) head3 = nbsLevels_[3].head();
if(numLevels_>4) head4 = nbsLevels_[4].head();
if(numLevels_>5) head5 = nbsLevels_[5].head();
if(numLevels_>6) head6 = nbsLevels_[6].head();
Kokkos::parallel_for(
"NBSLevels::build",
rPolicy,
LAMBDA_HD(int32 i){
int8 lvl = particleLevel[i];
auto ind = nbsLevel0.pointIndex(pointPosition[i]);
ind = mapIndexLevels(ind, 0, lvl);
int32 old;
if(lvl==0)
old =Kokkos::atomic_exchange(&head0(ind.x(), ind.y(), ind.z()),i);
else if(lvl==1)
old =Kokkos::atomic_exchange(&head1(ind.x(), ind.y(), ind.z()),i);
else if(lvl==2)
old =Kokkos::atomic_exchange(&head2(ind.x(), ind.y(), ind.z()),i);
else if(lvl==3)
old =Kokkos::atomic_exchange(&head3(ind.x(), ind.y(), ind.z()),i);
else if(lvl==4)
old =Kokkos::atomic_exchange(&head4(ind.x(), ind.y(), ind.z()),i);
else if(lvl==5)
old =Kokkos::atomic_exchange(&head5(ind.x(), ind.y(), ind.z()),i);
else if(lvl==6)
old =Kokkos::atomic_exchange(&head6(ind.x(), ind.y(), ind.z()),i);
next(i) = old;
});
Kokkos::fence();
}
template<typename IncludeFunction>
INLINE_FUNCTION_H
void build(range activeRange, IncludeFunction incld)
{
// nullify next and heads
findParticleLevel(activeRange.first, activeRange.second);
manageAllocateNext(activeRange);
nullify(activeRange);
rangePolicyType rPolicy(activeRange.first, activeRange.second);
auto nbsLevel0 = nbsLevels_[0];
auto pointPosition = nbsLevel0.pointPosition();
auto particleLevel = particleLevel_;
auto next = nbsLevel0.next();
auto head0 = nbsLevel0.head();
typename NBSLevelType::HeadType head1, head2, head3, head4, head5, head6;
if(numLevels_>1) head1 = nbsLevels_[1].head();
if(numLevels_>2) head2 = nbsLevels_[2].head();
if(numLevels_>3) head3 = nbsLevels_[3].head();
if(numLevels_>4) head4 = nbsLevels_[4].head();
if(numLevels_>5) head5 = nbsLevels_[5].head();
if(numLevels_>6) head6 = nbsLevels_[6].head();
Kokkos::parallel_for(
"NBSLevels::build",
rPolicy,
LAMBDA_HD(int32 i){
if(!incld(i)) return;
int8 lvl = particleLevel[i];
auto ind = nbsLevel0.pointIndex(pointPosition[i]);
ind = mapIndexLevels(ind, 0, lvl);
int32 old;
if(lvl==0)
old =Kokkos::atomic_exchange(&head0(ind.x(), ind.y(), ind.z()),i);
else if(lvl==1)
old =Kokkos::atomic_exchange(&head1(ind.x(), ind.y(), ind.z()),i);
else if(lvl==2)
old =Kokkos::atomic_exchange(&head2(ind.x(), ind.y(), ind.z()),i);
else if(lvl==3)
old =Kokkos::atomic_exchange(&head3(ind.x(), ind.y(), ind.z()),i);
else if(lvl==4)
old =Kokkos::atomic_exchange(&head4(ind.x(), ind.y(), ind.z()),i);
else if(lvl==5)
old =Kokkos::atomic_exchange(&head5(ind.x(), ind.y(), ind.z()),i);
else if(lvl==6)
old =Kokkos::atomic_exchange(&head6(ind.x(), ind.y(), ind.z()),i);
next(i) = old;
});
Kokkos::fence();
}
bool findParticleLevel(int32 first, int32 last)
{
if(last > particleLevel_.size())
{
reallocNoInit(particleLevel_,last);
}
auto diameter = nbsLevels_[0].diameter();
auto const maxSizes = maxSizeLevels_;
auto particleLevel = particleLevel_;
auto const sizeRatio = 0.999*nbsLevels_[0].sizeRatio();
int8 maxLvl = sizeRangeLevels_.size();
rangePolicyType rPolicy(first, last);
Kokkos::parallel_for(
"NBSLevels::findParticleLevel",
rPolicy,
LAMBDA_HD(int32 i)
{
for(int8 lvl = 0; lvl<maxLvl; lvl++)
{
if( sizeRatio*diameter[i]<= maxSizes[lvl] )
{
particleLevel[i] = lvl;
return;
}
}
particleLevel[i] = static_cast<int8>(-1);
});
Kokkos::fence();
return true;
}
}; //NBSLevels
}
#endif

101
src/Interaction/contactSearch/methods/NBSLoop.hpp
View File

@ -1,101 +0,0 @@
/*------------------------------- 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.
-----------------------------------------------------------------------------*/
int32 m = this->head_(i,j,k);
CellType currentCell(i,j,k);
int32 n = -1;
while( m > -1 )
{
auto p_m = this->pointPosition_[m];
auto d_m = sizeRatio_* diameter_[m];
// the same cell
n = this->next_(m);
while(n >-1)
{
auto p_n = this->pointPosition_[n];
auto d_n = sizeRatio_*diameter_[n];
if( sphereSphereCheck(p_m, p_n, d_m, d_n) )
{
auto ln = n;
auto lm = m;
if(lm>ln) Swap(lm,ln);
if( auto res = pairs.insert(lm,ln); res <0)
{
getFullUpdate++;
}
}
n = this->next_(n);
}
// neighbor cells
CellType neighborCell;
for(int32 ni=0; ni<13; ni++)
{
if(ni==0) neighborCell = currentCell + CellType( 0, 0,-1);
else if(ni==1) neighborCell = currentCell + CellType(-1, 0,-1);
else if(ni==2) neighborCell = currentCell + CellType(-1, 0, 0);
else if(ni==3) neighborCell = currentCell + CellType(-1, 0, 1);
else if(ni==4) neighborCell = currentCell + CellType( 0,-1,-1);
else if(ni==5) neighborCell = currentCell + CellType( 0,-1, 0);
else if(ni==6) neighborCell = currentCell + CellType( 0,-1, 1);
else if(ni==7) neighborCell = currentCell + CellType(-1,-1,-1);
else if(ni==8) neighborCell = currentCell + CellType(-1,-1, 0);
else if(ni==9) neighborCell = currentCell + CellType(-1,-1, 1);
else if(ni==10) neighborCell = currentCell + CellType( 1,-1,-1);
else if(ni==11) neighborCell = currentCell + CellType( 1,-1, 0);
else if(ni==12) neighborCell = currentCell + CellType( 1,-1, 1);
if( this->isInRange(neighborCell) )
{
n = this->head_(neighborCell.x(), neighborCell.y(), neighborCell.z());
while( n>-1)
{
auto p_n = this->pointPosition_[n];
auto d_n = sizeRatio_*diameter_[n];
if(sphereSphereCheck(p_m, p_n, d_m, d_n))
{
auto ln = n;
auto lm = m;
if(lm>ln) Swap(lm,ln);
if( auto res = pairs.insert(lm,ln); res <0)
{
getFullUpdate++;
}
}
n = this->next_[n];
}
}
}
m = this->next_[m];
}

67
src/Interaction/contactSearch/methods/cellBased/NBS/NBS.cpp Normal file
View File

@ -0,0 +1,67 @@
/*------------------------------- 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.
-----------------------------------------------------------------------------*/
#include "NBS.hpp"
pFlow::NBS::NBS
(
const dictionary& dict,
const box& domainBox,
real minBSSize,
real maxBSSize,
const deviceViewType1D<realx3> &position,
const pFlagTypeDevice &flags,
const deviceViewType1D<real> &diam,
uint32 nWallPoints,
uint32 nWallElements,
const ViewType1D<realx3,memory_space>& wallPoints,
const ViewType1D<uint32x3,memory_space>& wallVertices
)
:
particleWallContactSearchs<NBS>(
dict,
domainBox,
minBSSize,
maxBSSize,
position,
flags,
diam),
sizeRatio_(max(dict.getVal<real>("sizeRatio"), 1.0)),
cellExtent_(max(dict.getVal<real>("cellExtent"), 0.5)),
adjustableBox_(dict.getVal<Logical>("adjustableBox")),
NBSLevel0_
(
this->domainBox_,
maxBSSize,
sizeRatio_,
position,
flags,
adjustableBox_()
),
cellsWallLevel0_
(
cellExtent_,
nWallPoints,
nWallElements,
wallPoints,
wallVertices
)
{
}

148
src/Interaction/contactSearch/methods/cellBased/NBS/NBS.hpp Normal file
View File

@ -0,0 +1,148 @@
/*------------------------------- 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 __NBS_hpp__
#define __NBS_hpp__
#include "particleWallContactSearchs.hpp"
#include "NBSLevel0.hpp"
#include "cellsWallLevel0.hpp"
#include "Vector.hpp"
namespace pFlow
{
class NBS
:
public particleWallContactSearchs<NBS>
{
public:
using CellIterator = typename NBSLevel0::CellIterator;
private:
real sizeRatio_ = 1.0;
real cellExtent_ = 0.5;
Logical adjustableBox_;
NBSLevel0 NBSLevel0_;
cellsWallLevel0 cellsWallLevel0_;
protected:
friend particleWallContactSearchs<NBS>;
bool impl_broadSearch
(
csPairContainerType& ppPairs,
csPairContainerType& pwPairs,
const deviceViewType1D<realx3>& pointPos,
const pFlagTypeDevice& flags,
const deviceViewType1D<real>& diameter
)
{
bool searchBoxChanged;
if( !NBSLevel0_.broadSearch(
ppPairs,
pointPos,
flags,
diameter,
searchBoxChanged))
{
fatalErrorInFunction<<
"Error in broadSearch for NBS (particle-particle)"<<endl;
return false;
}
if(!cellsWallLevel0_.broadSearch(
pwPairs,
NBSLevel0_.getSearchCells(),
NBSLevel0_.getCellIterator()))
{
fatalErrorInFunction<<
"Error in broadSearch for NBS (particle-wall)"<<endl;
return false;
}
return true;
}
public:
TypeInfoNV("NBS");
NBS(
const dictionary& dict,
const box& domainBox,
real minBSSize,
real maxBSSize,
const deviceViewType1D<realx3> &position,
const pFlagTypeDevice &flags,
const deviceViewType1D<real> &diam,
uint32 nWallPoints,
uint32 nWallElements,
const ViewType1D<realx3,memory_space>& wallPoints,
const ViewType1D<uint32x3,memory_space>& wallVertices);
INLINE_FUNCTION_HD
NBS(const NBS&) = default;
INLINE_FUNCTION_HD
NBS(NBS&&) = default;
INLINE_FUNCTION_HD
NBS& operator = (const NBS&) = default;
INLINE_FUNCTION_HD
NBS& operator = (NBS&&) = default;
INLINE_FUNCTION_HD
~NBS()=default;
//// - Methods
uint32 numLevels()const
{
return 1;
}
auto getCellIterator([[maybe_unused]] uint32 lvl)const
{
return NBSLevel0_.getCellIterator();
}
Vector<cells> getDomainCellsLevels()const
{
return Vector<cells>("Cells", 1, NBSLevel0_.getDomainCells());
}
};
}
#endif

118
src/Interaction/contactSearch/methods/cellBased/NBS/NBSLevel0.cpp Normal file
View File

@ -0,0 +1,118 @@
/*------------------------------- 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.
-----------------------------------------------------------------------------*/
#include "NBSLevel0.hpp"
#include "streams.hpp"
bool pFlow::NBSLevel0::findPairs
(
csPairContainerType &pairs,
const deviceViewType1D<realx3> &pointPos,
const pFlagTypeDevice &flags,
const deviceViewType1D<real> &diameter
)
{
uint32 getFull = 1;
// loop until the container size fits the numebr of contact pairs
while (getFull > 0)
{
getFull = pFlow::NBSLevel0Kernels::findPairsCount
(
pairs,
sizeRatio_,
pointPos,
flags,
diameter,
getCellIterator()
);
if(getFull)
{
// - resize the container
// note that getFull now shows the number of failed insertions.
uint32 len = max(getFull,500u) ;
auto oldCap = pairs.capacity();
pairs.increaseCapacityBy(len);
INFORMATION<< "The contact pair container capacity increased from "<<
oldCap << " to "<<pairs.capacity()<<" in NBSLevel0."<<END_INFO;
}
Kokkos::fence();
}
return true;
}
pFlow::NBSLevel0::NBSLevel0
(
const box& domain,
real cellSize,
real sizeRatio,
const deviceViewType1D<realx3>& pointPos,
const pFlagTypeDevice& flags,
bool adjustableBox
)
:
mapperNBS
(
domain,
cellSize,
pointPos,
flags,
adjustableBox,
true
),
sizeRatio_(sizeRatio)
{
}
bool pFlow::NBSLevel0::broadSearch
(
csPairContainerType &pairs,
const deviceViewType1D<realx3> &pointPos,
const pFlagTypeDevice &flags,
const deviceViewType1D<real> &diameter,
bool& searchBoxChanged
)
{
if(!build(pointPos, flags, searchBoxChanged))
{
fatalErrorInFunction;
return false;
}
if(!findPairs(pairs, pointPos, flags, diameter))
{
fatalErrorInFunction;
return false;
}
return true;
}

98
src/Interaction/contactSearch/methods/cellBased/NBS/NBSLevel0.hpp Normal file
View File

@ -0,0 +1,98 @@
/*------------------------------- 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 __NBSLevel0_hpp__
#define __NBSLevel0_hpp__
#include "NBSLevel0Kernels.hpp"
#include "mapperNBS.hpp"
namespace pFlow
{
class NBSLevel0
:
public mapperNBS
{
public:
using MapperType = mapperNBS;
using CellIterator = typename MapperType::CellIterator;
using HeadType = typename MapperType::HeadType;
using NextType = typename MapperType::NextType;
private:
real sizeRatio_ = 1.0;
bool findPairs
(
csPairContainerType& pairs,
const deviceViewType1D<realx3>& pointPos,
const pFlagTypeDevice& flags,
const deviceViewType1D<real>& diameter
);
public:
TypeInfoNV("NBSLevel0");
INLINE_FUNCTION_HD
NBSLevel0() = default;
NBSLevel0(
const box& domain,
real cellSize,
real sizeRatio,
const deviceViewType1D<realx3>& pointPos,
const pFlagTypeDevice& flags,
bool adjustableBox);
INLINE_FUNCTION_HD
NBSLevel0(const NBSLevel0&) = default;
INLINE_FUNCTION_HD
NBSLevel0& operator = (const NBSLevel0&) = default;
INLINE_FUNCTION_HD
~NBSLevel0()=default;
bool broadSearch(
csPairContainerType& pairs,
const deviceViewType1D<realx3>& pointPos,
const pFlagTypeDevice& flags,
const deviceViewType1D<real>& diameter,
bool& searchBoxChanged);
};
} // pFlow
#endif // __NBSLevel0_hpp__

83
src/Interaction/contactSearch/methods/cellBased/NBS/NBSLevel0Kernels.hpp Normal file
View File

@ -0,0 +1,83 @@
/*------------------------------- 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.
-----------------------------------------------------------------------------*/
#include "mapperNBS.hpp"
#include "contactSearchGlobals.hpp"
#include "streams.hpp"
namespace pFlow::NBSLevel0Kernels
{
using mdrPolicyFindPairs =
Kokkos::MDRangePolicy<
Kokkos::IndexType<uint32>,
Kokkos::Rank<3>,
Kokkos::Schedule<Kokkos::Dynamic>,
DefaultExecutionSpace>;
template<typename T>
INLINE_FUNCTION_HD
void Swap(T& x, T& y)
{
T tmp = x;
x = y;
y = tmp;
}
INLINE_FUNCTION_HD
bool sphereSphereCheck(const realx3& p1, const realx3 p2, real d1, real d2)
{
return length(p2-p1) < 0.5*(d2+d1);
}
inline
uint32 findPairsCount
(
csPairContainerType& pairs,
real sizeRatio,
const deviceViewType1D<realx3>& pointPos,
const pFlagTypeDevice& flags,
const deviceViewType1D<real>& diameter,
mapperNBS::CellIterator cellIter
)
{
auto nCells = cellIter.numCells();
mdrPolicyFindPairs
mdrPolicy(
{0,0,0},
{nCells.x(), nCells.y(), nCells.z()} );
uint32 notInsertedPairs = 0u;
Kokkos::parallel_reduce (
"pFlow::NBSLevel0Kernels::findPairsCount",
mdrPolicy,
LAMBDA_HD(uint32 i, uint32 j, uint32 k, uint32& getFullUpdate){
#include "NBSLoop.hpp"
}, notInsertedPairs);
return notInsertedPairs;
}
}

100
src/Interaction/contactSearch/methods/cellBased/NBS/NBSLoop.hpp Normal file
View File

@ -0,0 +1,100 @@
/*------------------------------- 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.
-----------------------------------------------------------------------------*/
uint32 m = cellIter.start(i,j,k);
int32x3 currentCell(i,j,k);
uint32 n = mapperNBS::NoPos;
while( m != mapperNBS::NoPos)
{
auto p_m = pointPos[m];
auto d_m = sizeRatio* diameter[m];
// the same cell
n = cellIter.getNext(m);
while(n != mapperNBS::NoPos)
{
auto p_n = pointPos[n];
auto d_n = sizeRatio*diameter[n];
if( sphereSphereCheck(p_m, p_n, d_m, d_n) )
{
auto ln = n;
auto lm = m;
if(lm>ln) Swap(lm,ln);
if( pairs.insert(lm,ln) == -1)
{
getFullUpdate++;
}
}
n = cellIter.getNext(n);;
}
// neighbor cells
int32x3 neighborCell;
for(uint32 ni=0u; ni<13; ni++)
{
if(ni==0) neighborCell = currentCell + int32x3( 0, 0,-1);
else if(ni==1) neighborCell = currentCell + int32x3(-1, 0,-1);
else if(ni==2) neighborCell = currentCell + int32x3(-1, 0, 0);
else if(ni==3) neighborCell = currentCell + int32x3(-1, 0, 1);
else if(ni==4) neighborCell = currentCell + int32x3( 0,-1,-1);
else if(ni==5) neighborCell = currentCell + int32x3( 0,-1, 0);
else if(ni==6) neighborCell = currentCell + int32x3( 0,-1, 1);
else if(ni==7) neighborCell = currentCell + int32x3(-1,-1,-1);
else if(ni==8) neighborCell = currentCell + int32x3(-1,-1, 0);
else if(ni==9) neighborCell = currentCell + int32x3(-1,-1, 1);
else if(ni==10) neighborCell = currentCell + int32x3( 1,-1,-1);
else if(ni==11) neighborCell = currentCell + int32x3( 1,-1, 0);
else if(ni==12) neighborCell = currentCell + int32x3( 1,-1, 1);
if( neighborCell.x()>=0 && neighborCell.y()>=0 && neighborCell.z()>=0 &&
neighborCell.x()<nCells.x() && neighborCell.y()<nCells.y() && neighborCell.z()<nCells.z() )
{
n = cellIter.start(neighborCell.x(), neighborCell.y(), neighborCell.z());
while(n != mapperNBS::NoPos)
{
auto p_n = pointPos[n];
auto d_n = sizeRatio*diameter[n];
if(sphereSphereCheck(p_m, p_n, d_m, d_n))
{
auto ln = n;
auto lm = m;
if(lm>ln) Swap(lm,ln);
if( pairs.insert(lm,ln) == -1)
{
getFullUpdate++;
}
}
n = cellIter.next(n);
}
}
}
m = cellIter.next(m);
}

195
src/Interaction/contactSearch/methods/cellBased/NBS/cellsWallLevel0.cpp Normal file
View File

@ -0,0 +1,195 @@
/*------------------------------- 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.
-----------------------------------------------------------------------------*/
#include "cellsWallLevel0.hpp"
#include "streams.hpp"
pFlow::cellsWallLevel0::cellsWallLevel0
(
real cellExtent,
uint32 numPoints,
uint32 numElements,
const ViewType1D<realx3,
memory_space> &points,
const ViewType1D<uint32x3, memory_space> &vertices
)
:
cellExtent_( max(cellExtent, 0.5 ) ),
numElements_(numElements),
numPoints_(numPoints),
vertices_(vertices),
points_(points)
{
allocateArrays();
}
bool pFlow::cellsWallLevel0::resetElements
(
uint32 numElements,
uint32 numPoints,
ViewType1D<realx3, memory_space> &points,
ViewType1D<uint32x3, memory_space> &vertices
)
{
numElements_ = numElements;
numPoints_ = numPoints;
points_ = points;
vertices_ = vertices;
allocateArrays();
return true;
}
bool pFlow::cellsWallLevel0::broadSearch
(
csPairContainerType &pairs,
const cells& searchBox,
const mapperNBS::CellIterator &particleMap
)
{
// map walls onto the cells
this->build(searchBox);
this->particleWallFindPairs(pairs, particleMap);
return true;
}
bool pFlow::cellsWallLevel0::build(const cells & searchBox)
{
Kokkos::parallel_for(
"pFlow::cellsWallLevel0::build",
deviceRPolicyStatic(0,numElements_),
CLASS_LAMBDA_HD(uint32 i)
{
auto v = vertices_[i];
auto p1 = points_[v.x()];
auto p2 = points_[v.y()];
auto p3 = points_[v.z()];
realx3 minP;
realx3 maxP;
searchBox.extendBox(p1, p2, p3, cellExtent_, minP, maxP);
elementBox_[i] = iBoxType(searchBox.pointIndex(minP), searchBox.pointIndex(maxP));
});
Kokkos::fence();
return true;
}
bool pFlow::cellsWallLevel0::particleWallFindPairs
(
csPairContainerType &pairs,
const mapperNBS::CellIterator &particleMap
)
{
uint32 getFull = 1;
while (getFull)
{
getFull = findPairsElementRangeCount(pairs, particleMap);
if(getFull)
{
// - resize the container
// note that getFull now shows the number of failed insertions.
uint32 len = max(getFull, 50u);
auto oldCap = pairs.capacity();
pairs.increaseCapacityBy(len);
INFORMATION<<"Contact pair container capacity increased from "<<
oldCap << " to "
<< pairs.capacity() <<" in cellsWallLevel0."<<END_INFO;
Kokkos::fence();
}
}
return true;
}
pFlow::int32 pFlow::cellsWallLevel0::findPairsElementRangeCount
(
csPairContainerType &pairs,
const mapperNBS::CellIterator &particleMap
)
{
uint32 getFull =0;
//const auto pwPairs = pairs;
const auto elementBox = elementBox_;
Kokkos::parallel_reduce(
"pFlow::cellsWallLevel0::findPairsElementRangeCount",
tpPWContactSearch(numElements_, Kokkos::AUTO),
LAMBDA_HD(
const typename tpPWContactSearch::member_type & teamMember,
uint32& valueToUpdate){
const uint32 iTri = teamMember.league_rank();
const auto triBox = elementBox[iTri];
uint32 getFull2 = 0;
auto bExtent = boxExtent(triBox);
uint32 numCellBox = bExtent.x()*bExtent.y()*bExtent.z();
Kokkos::parallel_reduce(
Kokkos::TeamThreadRange( teamMember, numCellBox ),
[&] ( const uint32 linIndex, uint32 &innerUpdate )
{
int32x3 cell;
indexToCell(linIndex, triBox, cell);
uint32 n = particleMap.start(cell.x(),cell.y(),cell.z());
while( n != particleMap.NoPos)
{
// id is wall id the pair is (particle id, wall id)
if( pairs.insert(
static_cast<csIdType>(n),
static_cast<csIdType>(iTri) ) == -1 )
innerUpdate++;
n = particleMap.next(n);
}
},
getFull2
);
if ( teamMember.team_rank() == 0 ) valueToUpdate += getFull2;
},
getFull
);
return getFull;
}

139
src/Interaction/contactSearch/methods/cellBased/NBS/cellsWallLevel0.hpp Normal file
View File

@ -0,0 +1,139 @@
/*------------------------------- 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 __cellsWallLevel0_hpp__
#define __cellsWallLevel0_hpp__
#include "contactSearchGlobals.hpp"
#include "contactSearchFunctions.hpp"
#include "mapperNBS.hpp"
#include "iBox.hpp"
namespace pFlow
{
class cellsWallLevel0
{
public:
using execution_space = csExecutionSpace;
using memory_space = typename execution_space::memory_space;
using iBoxType = iBox<int32>;
class TagFindCellRange2{};
private:
// - box extent
real cellExtent_ = 0.5;
// - number of triangle elements
uint32 numElements_ = 0;
// - number of points
uint32 numPoints_ = 0;
// - ref to vectices (borrowed)
ViewType1D<uint32x3, memory_space> vertices_;
// - ref to points in the trisurface (borrowed)
ViewType1D<realx3, memory_space> points_;
// cell range of element/triangle bounding box
ViewType1D<iBoxType, memory_space> elementBox_;
using tpPWContactSearch = Kokkos::TeamPolicy<
execution_space,
Kokkos::Schedule<Kokkos::Dynamic>,
Kokkos::IndexType<uint32>>;
FUNCTION_H
void allocateArrays()
{
reallocNoInit( elementBox_, numElements_);
}
public:
TypeInfoNV("cellsWallLevel0");
INLINE_FUNCTION_HD
cellsWallLevel0(){}
FUNCTION_H
cellsWallLevel0(
real cellExtent,
uint32 numPoints,
uint32 numElements,
const ViewType1D<realx3,memory_space>& points,
const ViewType1D<uint32x3,memory_space>& vertices);
// - host call
// reset triangle elements if they have changed
bool resetElements(
uint32 numElements,
uint32 numPoints,
ViewType1D<realx3, memory_space>& points,
ViewType1D<uint32x3, memory_space>& vertices);
INLINE_FUNCTION_HD
iBoxType elementBox(uint32 i)const
{
return elementBox_[i];
}
INLINE_FUNCTION_HD
uint32 numElements()const
{
return numElements_;
}
bool broadSearch(
csPairContainerType& pairs,
const cells& searchBox,
const mapperNBS::CellIterator& particleMap);
bool build(const cells& searchBox);
bool particleWallFindPairs(
csPairContainerType& pairs,
const mapperNBS::CellIterator& particleMap);
int32 findPairsElementRangeCount(
csPairContainerType& pairs,
const mapperNBS::CellIterator& particleMap);
}; // cellsWallLevel0
} // pFlow
#endif // __cellsWallLevel0_hpp__

194
src/Interaction/contactSearch/methods/cellBased/NBS/mapperNBS.cpp Normal file
View File

@ -0,0 +1,194 @@
/*------------------------------- 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.
-----------------------------------------------------------------------------*/
#include "mapperNBS.hpp"
#include "mapperNBSKernels.hpp"
#include "streams.hpp"
pFlow::uint32 pFlow::mapperNBS::checkInterval_ = 1000;
pFlow::real pFlow::mapperNBS::enlargementFactor_ = 1.1;
bool pFlow::mapperNBS::setSearchBox
(
const deviceViewType1D<realx3> &pointPos,
const pFlagTypeDevice &flags,
real cellSize
)
{
box domainBox = domainCells_.domainBox();
if(adjustableBox_)
{
lastCheckForBox_ = buildCount_;
realx3 minP;
realx3 maxP;
pFlow::mapperNBSKernels::findPointExtends
(
pointPos,
flags,
minP, maxP
);
minP = max( minP - enlargementFactor_*cellSize, domainBox.minPoint());
maxP = min( maxP + enlargementFactor_*cellSize, domainBox.maxPoint());
box searchBox = {minP, maxP};
searchCells_ = cells(searchBox, cellSize);
INFORMATION<<"Search box for contact search has changed: "<<
"search box is ["<<searchCells_.domainBox().minPoint()<<
" "<<searchCells_.domainBox().maxPoint()<<"]"<<END_INFO;
return true;
}
else
{
searchCells_ = cells(domainBox, cellSize);
INFORMATION<<"Search box for contact search is: ["
<< domainBox.minPoint()<<" "<<domainBox.maxPoint()<<"]"<<END_INFO;
return false;
}
}
void pFlow::mapperNBS::allocateArrays(rangeU32 nextRng)
{
checkAllocateNext(nextRng);
nullifyNext(nextRng);
reallocFill(head_, searchCells_.nx(), searchCells_.ny(), searchCells_.nz(), NoPos);
}
void pFlow::mapperNBS::checkAllocateNext(rangeU32 nextRng)
{
auto newCap = nextRng.end();
if( nextCapacity_ < newCap)
{
nextCapacity_ = newCap;
if(!nextOwner_)return;
reallocNoInit(next_, nextCapacity_);
}
}
void pFlow::mapperNBS::nullifyHead()
{
fill(head_, NoPos);
}
void pFlow::mapperNBS::nullifyNext(rangeU32 nextRng)
{
if(!nextOwner_)return;
fill(next_, nextRng, NoPos);
}
pFlow::mapperNBS::mapperNBS(
const box &domain,
real cellSize,
const deviceViewType1D<realx3> &pointPos,
const pFlagTypeDevice &flags,
bool adjustableBox,
bool nextOwner
)
:
domainCells_(domain, cellSize),
searchCells_(domain, cellSize),
adjustableBox_(adjustableBox),
nextOwner_(nextOwner)
{
setSearchBox(pointPos, flags, cellSize);
allocateArrays(flags.activeRange());
}
bool pFlow::mapperNBS::build
(
const deviceViewType1D<realx3>& pointPos,
const pFlagTypeDevice & flags,
bool& searchBoxChanged
)
{
auto aRange = flags.activeRange();
buildCount_++;
if(adjustableBox_ && buildCount_%checkInterval_ == 0)
{
if(searchBoxChanged =
setSearchBox(pointPos, flags, searchCells_.cellSize());searchBoxChanged)
{
allocateArrays(aRange);
}
}
else
{
checkAllocateNext(aRange);
nullifyHead();
nullifyNext(aRange);
}
if( adjustableBox_ )
{
if(!pFlow::mapperNBSKernels::buildListsReduce(
searchCells_,
head_,
next_,
pointPos,
flags) )
{
buildCount_++;
setSearchBox(pointPos, flags, searchCells_.cellSize());
searchBoxChanged = true;
allocateArrays(flags.activeRange());
if(!pFlow::mapperNBSKernels::buildListsReduce(
searchCells_,
head_,
next_,
pointPos,
flags))
{
fatalErrorInFunction<<"failed to build list in anjustable search box mode!"<<endl;
return false;
}
}
}
else
{
pFlow::mapperNBSKernels::buildLists(
searchCells_,
head_,
next_,
pointPos,
flags
);
searchBoxChanged = false;
}
return true;
}

212
src/Interaction/contactSearch/methods/cellBased/NBS/mapperNBS.hpp Normal file
View File

@ -0,0 +1,212 @@
/*------------------------------- 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 __mapperNBS_hpp__
#define __mapperNBS_hpp__
#include "phasicFlowKokkos.hpp"
#include "pointFlag.hpp"
#include "cells.hpp"
//#include "contactSearchFunctions.hpp"
namespace pFlow
{
class mapperNBS
{
public:
using HeadType = deviceViewType3D<uint32>;
using NextType = deviceViewType1D<uint32>;
static constexpr uint32 NoPos = 0xFFFFFFFF;
class CellIterator
{
private:
HeadType head_;
NextType next_;
public:
CellIterator(const HeadType& head, const NextType& next)
:
head_(head),
next_(next)
{}
static constexpr uint32 NoPos = 0xFFFFFFFF;
INLINE_FUNCTION_HD
int32x3 numCells()const {
return int32x3(head_.extent(0), head_.extent(1), head_.extent(2));}
INLINE_FUNCTION_HD
uint32 start(int32 i, int32 j, int32 k)const {
return head_(i,j,k); }
INLINE_FUNCTION_HD
uint32 getNext(uint32 n)const {
if(n == NoPos ) return NoPos;
return next_(n); }
INLINE_FUNCTION_HD
uint32 next(uint32 n)const{
return next_(n);}
};
private:
cells domainCells_;
cells searchCells_;
HeadType head_{"NBS::head",1,1,1};
NextType next_{"NBS::next", 1};
uint32 nextCapacity_ = 0;
uint32 lastCheckForBox_ = 0;
uint32 buildCount_ = 0;
bool adjustableBox_ = false;
bool nextOwner_ = true;
static uint32 checkInterval_;
static real enlargementFactor_;
bool setSearchBox(
const deviceViewType1D<realx3>& pointPos,
const pFlagTypeDevice& flags,
real cellSize
);
void allocateArrays(rangeU32 nextRng);
void checkAllocateNext(rangeU32 nextRng);
void nullifyHead();
void nullifyNext(rangeU32 nextRng);
public:
TypeInfoNV("mapperNBS");
INLINE_FUNCTION_HD
mapperNBS() = default;
mapperNBS(
const box& domain,
real cellSize,
const deviceViewType1D<realx3>& pointPos,
const pFlagTypeDevice& flags,
bool adjustableBox,
bool nextOwner = true);
INLINE_FUNCTION_HD
mapperNBS(const mapperNBS&) = default;
INLINE_FUNCTION_HD
mapperNBS(mapperNBS&&) = default;
INLINE_FUNCTION_HD
mapperNBS& operator = (const mapperNBS&) = default;
INLINE_FUNCTION_HD
mapperNBS& operator = (mapperNBS&&) = default;
INLINE_FUNCTION_HD
~mapperNBS()=default;
//// - Methods
auto getCellIterator()const
{
return CellIterator(head_, next_);
}
const auto& getDomainCells()const
{
return domainCells_;
}
const auto& getSearchCells()const
{
return searchCells_;
}
bool build(
const deviceViewType1D<realx3>& pointPos,
const pFlagTypeDevice& flags,
bool& searchBoxChanged);
};
} // pFlow
#endif // __mapperNBS_hpp__
/*
INLINE_FUNCTION_HD
auto& head()
{
return head_;
}
INLINE_FUNCTION_HD
auto& next()
{
return next_;
}
INLINE_FUNCTION_HD
const auto& head()const
{
return head_;
}
INLINE_FUNCTION_HD
const auto& next()const
{
return next_;
}
INLINE_FUNCTION_HD
auto& pointPosition()
{
return pointPosition_;
}
*/

192
src/Interaction/contactSearch/methods/cellBased/NBS/mapperNBSKernels.cpp Normal file
View File

@ -0,0 +1,192 @@
/*------------------------------- 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.
-----------------------------------------------------------------------------*/
#include "mapperNBSKernels.hpp"
void pFlow::mapperNBSKernels::findPointExtends
(
const deviceViewType1D<realx3>& points,
const pFlagTypeDevice& flags,
realx3& minPoint,
realx3& maxPoint
)
{
if(flags.numActive() == 0)
{
minPoint = {0,0,0};
maxPoint = {0,0,0};
return;
}
real minX;
real minY;
real minZ;
real maxX;
real maxY;
real maxZ;
auto aRange = flags.activeRange();
Kokkos::parallel_reduce(
"pFlow::mapperNBSKernels::findPointExtends",
deviceRPolicyStatic(aRange.start(), aRange.end()),
LAMBDA_HD(
uint32 i,
real& minXUpdate,
real& minYUpdate,
real& minZUpdate,
real& maxXUpdate,
real& maxYUpdate,
real& maxZUpdate)
{
if(flags(i))
{
auto p = points(i);
minXUpdate = min(p.x(), minXUpdate);
minYUpdate = min(p.y(), minYUpdate);
minZUpdate = min(p.z(), minZUpdate);
maxXUpdate = max(p.x(), maxXUpdate);
maxYUpdate = max(p.y(), maxYUpdate);
maxZUpdate = max(p.z(), maxZUpdate);
}
},
Kokkos::Min<real>(minX),
Kokkos::Min<real>(minY),
Kokkos::Min<real>(minZ),
Kokkos::Max<real>(maxX),
Kokkos::Max<real>(maxY),
Kokkos::Max<real>(maxZ)
);
minPoint = {minX, minY, minZ};
maxPoint = {maxX, maxY, maxZ};
}
bool pFlow::mapperNBSKernels::buildListsReduce
(
const cells &searchCell,
const deviceViewType3D<uint32> &head,
const deviceViewType1D<uint32> &next,
const deviceViewType1D<realx3> &points,
const pFlagTypeDevice &flags
)
{
uint32 numOut = 0;
auto aRange = flags.activeRange();
if(flags.isAllActive())
{
Kokkos::parallel_reduce
(
"pFlow::mapperNBSKernels::buildListsReduce",
deviceRPolicyStatic(aRange.start(), aRange.end()),
LAMBDA_HD(uint32 i, uint32& valToUpdate)
{
int32x3 ind;
if( searchCell.pointIndexInDomain(points[i], ind) )
{
uint32 old = Kokkos::atomic_exchange(&head(ind.x(), ind.y(), ind.z()), i);
next[i] = old;
}
else
{
valToUpdate++;
}
},
numOut
);
}
else
{
Kokkos::parallel_reduce
(
"pFlow::mapperNBSKernels::buildListsReduce",
deviceRPolicyStatic(aRange.start(), aRange.end()),
LAMBDA_HD(uint32 i, uint32& valToUpdate)
{
int32x3 ind;
if( flags(i) )
{
if( searchCell.pointIndexInDomain(points[i], ind) )
{
uint32 old = Kokkos::atomic_exchange(&head(ind.x(), ind.y(), ind.z()), i);
next[i] = old;
}
else
{
valToUpdate++;
}
}
},
numOut
);
}
return numOut == 0u ;
}
bool pFlow::mapperNBSKernels::buildLists
(
const cells &searchCell,
const deviceViewType3D<uint32> &head,
const deviceViewType1D<uint32> &next,
const deviceViewType1D<realx3> &points,
const pFlagTypeDevice &flags
)
{
auto aRange = flags.activeRange();
if(flags.isAllActive() )
{
Kokkos::parallel_for
(
"pFlow::mapperNBSKernels::buildLists",
deviceRPolicyStatic(aRange.start(), aRange.end()),
LAMBDA_HD(uint32 i)
{
auto ind = searchCell.pointIndex(points[i]);
uint32 old = Kokkos::atomic_exchange(&head(ind.x(), ind.y(), ind.z()), i);
next[i] = old;
}
);
Kokkos::fence();
}
else
{
Kokkos::parallel_for
(
"pFlow::mapperNBSKernels::buildLists",
deviceRPolicyStatic(aRange.start(), aRange.end()),
LAMBDA_HD(uint32 i)
{
if( flags(i) )
{
auto ind = searchCell.pointIndex(points[i]);
uint32 old = Kokkos::atomic_exchange(&head(ind.x(), ind.y(), ind.z()), i);
next[i] = old;
}
}
);
Kokkos::fence();
}
return true;
}

49
src/Interaction/contactSearch/methods/cellBased/NBS/mapperNBSKernels.hpp Normal file
View File

@ -0,0 +1,49 @@
/*------------------------------- 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.
-----------------------------------------------------------------------------*/
#include "phasicFlowKokkos.hpp"
#include "cells.hpp"
#include "pointFlag.hpp"
namespace pFlow::mapperNBSKernels
{
void findPointExtends(
const deviceViewType1D<realx3>& points,
const pFlagTypeDevice& flags,
realx3& minPoint,
realx3& maxPoint);
bool buildListsReduce(
const cells& searchCell,
const deviceViewType3D<uint32>& head,
const deviceViewType1D<uint32>& next,
const deviceViewType1D<realx3>& points,
const pFlagTypeDevice& flags);
bool buildLists(
const cells& searchCell,
const deviceViewType3D<uint32>& head,
const deviceViewType1D<uint32>& next,
const deviceViewType1D<realx3>& points,
const pFlagTypeDevice& flags);
}

107
src/Interaction/contactSearch/cells.hpp → src/Interaction/contactSearch/methods/cellBased/cells.hpp
View File

@ -28,58 +28,40 @@ Licence:
namespace pFlow
{
template<typename indexType>
class cells
{
public:
using CellType = triple<indexType>;
protected:
private:
// - domain
box domain_{realx3(0.0), realx3(1.0)};
box domainBox_{realx3(0.0), realx3(1.0)};
// - cell size
realx3 cellSize_{1,1,1};
CellType numCells_{1,1,1};
real celldx_{1};
int32x3 numCells_{1,1,1};
// - protected methods
INLINE_FUNCTION_H
void calculate()
{
numCells_ = (domain_.maxPoint()-domain_.minPoint())/cellSize_ + realx3(1.0);
numCells_ = max( numCells_ , CellType(static_cast<indexType>(1)) );
numCells_ = (domainBox_.maxPoint()-domainBox_.minPoint())/celldx_ + realx3(1.0);
numCells_ = max( numCells_ , int32x3(1) );
}
public:
INLINE_FUNCTION_HD
cells()
{}
cells() = default;
INLINE_FUNCTION_H
cells(const box& domain, real cellSize)
:
domain_(domain),
cellSize_(cellSize)
domainBox_(domain),
celldx_(cellSize)
{
calculate();
}
INLINE_FUNCTION_H
cells(const box& domain, int32 nx, int32 ny, int32 nz)
:
domain_(domain),
cellSize_(
(domain_.maxPoint() - domain_.minPoint())/realx3(nx, ny, nz)
),
numCells_(nx, ny, nz)
{}
INLINE_FUNCTION_HD
cells(const cells&) = default;
@ -100,43 +82,36 @@ public:
INLINE_FUNCTION_H
void setCellSize(real cellSize)
{
cellSize_ = cellSize;
calculate();
}
INLINE_FUNCTION_H
void setCellSize(realx3 cellSize)
{
cellSize_ = cellSize;
celldx_ = cellSize;
calculate();
}
INLINE_FUNCTION_HD
realx3 cellSize()const
real cellSize()const
{
return cellSize_;
return celldx_;
}
INLINE_FUNCTION_HD
const CellType& numCells()const
const int32x3& numCells()const
{
return numCells_;
}
INLINE_FUNCTION_HD
indexType nx()const
int32 nx()const
{
return numCells_.x();
}
INLINE_FUNCTION_HD
indexType ny()const
int32 ny()const
{
return numCells_.y();
}
INLINE_FUNCTION_HD
indexType nz()const
int32 nz()const
{
return numCells_.z();
}
@ -149,22 +124,21 @@ public:
static_cast<int64>(numCells_.z());
}
const auto& domain()const
const auto& domainBox()const
{
return domain_;
return domainBox_;
}
INLINE_FUNCTION_HD
CellType pointIndex(const realx3& p)const
int32x3 pointIndex(const realx3& p)const
{
return CellType( (p - domain_.minPoint())/cellSize_ );
return int32x3( (p - domainBox_.minPoint())/celldx_ );
}
INLINE_FUNCTION_HD
bool pointIndexInDomain(const realx3 p, CellType& index)const
bool pointIndexInDomain(const realx3 p, int32x3& index)const
{
if( !domain_.isInside(p) ) return false;
if(!inDomain(p))return false;
index = this->pointIndex(p);
return true;
}
@ -172,11 +146,11 @@ public:
INLINE_FUNCTION_HD
bool inDomain(const realx3& p)const
{
return domain_.isInside(p);
return domainBox_.isInside(p);
}
INLINE_FUNCTION_HD
bool isInRange(const CellType& cell)const
bool inCellRange(const int32x3& cell)const
{
if(cell.x()<0)return false;
if(cell.y()<0)return false;
@ -188,7 +162,7 @@ public:
}
INLINE_FUNCTION_HD
bool isInRange(indexType i, indexType j, indexType k)const
bool inCellRange(int32 i, int32 j, int32 k)const
{
if(i<0)return false;
if(j<0)return false;
@ -199,22 +173,7 @@ public:
return true;
}
INLINE_FUNCTION_HD
void extendBox(
const CellType& p1,
const CellType& p2,
const CellType& p3,
indexType extent,
CellType& minP,
CellType& maxP)const
{
minP = min( min( p1, p2), p3)-extent;
maxP = max( max( p1, p2), p3)+extent;
minP = bound(minP);
maxP = bound(maxP);
}
INLINE_FUNCTION_HD
void extendBox(
const realx3& p1,
@ -224,17 +183,17 @@ public:
realx3& minP,
realx3& maxP)const
{
minP = min(min(p1,p2),p3) - extent*cellSize_ ;
maxP = max(max(p1,p2),p3) + extent*cellSize_ ;
minP = min(min(p1,p2),p3) - extent*celldx_ ;
maxP = max(max(p1,p2),p3) + extent*celldx_ ;
minP = bound(minP);
maxP = bound(maxP);
}
INLINE_FUNCTION_HD
CellType bound(CellType p)const
int32x3 bound(int32x3 p)const
{
return CellType(
return int32x3(
min( numCells_.x()-1, max(0,p.x())),
min( numCells_.y()-1, max(0,p.y())),
min( numCells_.z()-1, max(0,p.z()))
@ -245,9 +204,9 @@ public:
realx3 bound(realx3 p)const
{
return realx3(
min( domain_.maxPoint().x(), max(domain_.minPoint().x(),p.x())),
min( domain_.maxPoint().y(), max(domain_.minPoint().y(),p.y())),
min( domain_.maxPoint().z(), max(domain_.minPoint().z(),p.z()))
min( domainBox_.maxPoint().x(), max(domainBox_.minPoint().x(),p.x())),
min( domainBox_.maxPoint().y(), max(domainBox_.minPoint().y(),p.y())),
min( domainBox_.maxPoint().z(), max(domainBox_.minPoint().z(),p.z()))
);
}
};

389
src/Interaction/contactSearch/methods/mapperNBS.hpp
View File

@ -1,389 +0,0 @@
/*------------------------------- 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 __mapperNBS_hpp__
#define __mapperNBS_hpp__
#include "cells.hpp"
#include "contactSearchFunctions.hpp"
#include "baseAlgorithms.hpp"
#include "ViewAlgorithms.hpp"
namespace pFlow
{
template<typename executionSpace>
class mapperNBS
:
public cells<int32>
{
public:
using IdType = int32;
using IndexType = int32;
using Cells = cells<IndexType>;
using CellType = typename Cells::CellType;
using execution_space = executionSpace;
using memory_space = typename execution_space::memory_space;
using HeadType = ViewType3D<int32, memory_space>;
using NextType = ViewType1D<int32, memory_space>;
class cellIterator
{
private:
HeadType head_;
NextType next_;
public:
cellIterator(ViewType3D<int32, memory_space> head, ViewType1D<int32, memory_space> next)
:
head_(head),
next_(next)
{}
INLINE_FUNCTION_HD
Cells cellsSize()const {
return Cells(head_.extent(0), head_.extent(1), head_.extent(2));}
INLINE_FUNCTION_HD
int32 start(IndexType i, IndexType j, IndexType k)const {
return head_(i,j,k); }
INLINE_FUNCTION_HD
int32 getNext(int32 n)const {
if(n<0) return n;
return next_(n); }
};
protected:
int32 capacity_ = 1;
ViewType3D<int32, memory_space> head_;
ViewType1D<int32, memory_space> next_;
bool nextOwner_ = true;
// borrowed ownership
ViewType1D<realx3, memory_space> pointPosition_;
using rangePolicyType =
Kokkos::RangePolicy<
Kokkos::IndexType<int32>,
Kokkos::Schedule<Kokkos::Static>,
execution_space>;
INLINE_FUNCTION_H
void nullifyHead()
{
fill(
head_,
range(0,this->nx()),
range(0,this->ny()),
range(0,this->nz()),
static_cast<int32>(-1)
);
}
void nullifyNext(range nextRng)
{
if(!nextOwner_)return;
fill(
next_,
nextRng,
static_cast<int32>(-1)
);
}
void nullify()
{
nullifyHead();
nullifyNext(range(0,capacity_));
}
void nullify(range nextRng)
{
nullifyHead();
nullifyNext(nextRng);
}
void checkAllocateNext(int newCap)
{
if( capacity_ < newCap)
{
capacity_ = newCap;
if(!nextOwner_)return;
reallocNoInit(next_, capacity_);
}
}
void allocateHead()
{
reallocNoInit(head_, this->nx(), this->ny(), this->nz());
}
public:
TypeInfoNV("mapperNBS");
INLINE_FUNCTION_HD
mapperNBS(){}
mapperNBS(
const box& domain,
real cellSize,
const ViewType1D<realx3, memory_space>& position,
bool nextOwner = true)
:
Cells(domain, cellSize),
pointPosition_(position),
head_(
"mapperNBS::head_",
this->nx(),
this->ny(),
this->nz()
),
next_("mapperNBS::next_",1), //,position.size()),
nextOwner_(nextOwner)
{
checkAllocateNext(pointPosition_.size());
}
mapperNBS(
const box& domain,
int32 nx,
int32 ny,
int32 nz,
const ViewType1D<realx3, memory_space>& position,
bool nextOwner = true)
:
Cells(domain, nx, ny, nz),
pointPosition_(position),
head_("mapperNBS::head_",nx,ny,nz),
next_("mapperNBS::next_",1),
nextOwner_(nextOwner)
{
checkAllocateNext(pointPosition_.size());
}
INLINE_FUNCTION_HD
mapperNBS(const mapperNBS&) = default;
INLINE_FUNCTION_HD
mapperNBS& operator = (const mapperNBS&) = default;
INLINE_FUNCTION_HD
~mapperNBS()=default;
//// - Methods
INLINE_FUNCTION_HD
auto capacity()const
{
return capacity_;
}
cellIterator getCellIterator()const
{
return cellIterator(head_, next_);
}
bool particlesCapcityChanged(int32 newCap)
{
checkAllocateNext(newCap);
return true;
}
INLINE_FUNCTION_HD
auto& head()
{
return head_;
}
INLINE_FUNCTION_HD
auto& next()
{
return next_;
}
INLINE_FUNCTION_HD
const auto& head()const
{
return head_;
}
INLINE_FUNCTION_HD
const auto& next()const
{
return next_;
}
INLINE_FUNCTION_HD
auto& pointPosition()
{
return pointPosition_;
}
INLINE_FUNCTION_H
void setNext(ViewType1D<int32, memory_space>& next)
{
if(!nextOwner_)
{
next_ = next;
capacity_ = next.size();
}
}
// - build based on all points in active range
INLINE_FUNCTION_H
void build(range activeRange)
{
checkAllocateNext(activeRange.second);
nullify(activeRange);
Cells cellIndex = static_cast<Cells>(*this);
auto points = pointPosition_;
auto next = next_;
auto head = head_;
rangePolicyType rPolicy(activeRange.first, activeRange.second);
Kokkos::parallel_for(
"mapperNBS::build",
rPolicy,
LAMBDA_HD(int32 i){
CellType ind = cellIndex.pointIndex(points[i]);
int32 old = Kokkos::atomic_exchange(&head(ind.x(), ind.y(), ind.z()), i);
next[i] = old;
});
Kokkos::fence();
}
template<typename IncludeFunction>
INLINE_FUNCTION_H
void build(range activeRange, IncludeFunction incld)
{
checkAllocateNext(activeRange.second);
nullify(activeRange);
Cells cellIndex = static_cast<Cells>(*this);
auto points = pointPosition_;
auto next = next_;
auto head = head_;
rangePolicyType rPolicy(activeRange.first, activeRange.second);
Kokkos::parallel_for(
"mapperNBS::build_Include",
rPolicy,
LAMBDA_HD(int32 i){
if( incld(i) )
{
CellType ind = cellIndex.pointIndex(points[i]);
auto old = Kokkos::atomic_exchange(&head(ind.x(), ind.y(), ind.z()), i);
next[i] = old;
}
});
Kokkos::fence();
}
INLINE_FUNCTION_H
void buildCheckInDomain(range activeRange)
{
checkAllocateNext(activeRange.second);
nullify(activeRange);
Cells cellIndex = static_cast<Cells>(*this);
auto points = pointPosition_;
auto next = next_;
auto head = head_;
rangePolicyType rPolicy(activeRange.first, activeRange.second);
Kokkos::parallel_for(
"mapperNBS::buildCheckInDomain",
rPolicy,
LAMBDA_HD(int32 i){
CellType ind;
if( cellIndex.pointIndexInDomain(points[i], ind) )
{
int32 old = Kokkos::atomic_exchange(&head(ind.x(), ind.y(), ind.z()), i);
next[i] = old;
}
});
Kokkos::fence();
}
template<typename IncludeFunction>
INLINE_FUNCTION_H
void buildCheckInDomain(range activeRange, IncludeFunction incld)
{
checkAllocateNext(activeRange.second);
nullify(activeRange);
Cells cellIndex = static_cast<Cells>(*this);
auto points = pointPosition_;
auto next = next_;
auto head = head_;
rangePolicyType rPolicy(activeRange.first, activeRange.second);
Kokkos::parallel_for(
"mapperNBS::buildCheckInDomain_Include",
rPolicy,
LAMBDA_HD(int32 i){
CellType ind;
if( incld(i) && cellIndex.pointIndexInDomain(points[i], ind) )
{
auto old = Kokkos::atomic_exchange(&head(ind.x(), ind.y(), ind.z()), i);
next[i] = old;
}
});
Kokkos::fence();
}
};
} // pFlow
#endif // __mapperNBS_hpp__

213
src/Interaction/contactSearch/methods/multiGridNBS.hpp
View File

@ -1,213 +0,0 @@
/*------------------------------- 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 __multiGridNBS_hpp__
#define __multiGridNBS_hpp__
#include "NBSLevels.hpp"
namespace pFlow
{
template<typename executionSpace>
class multiGridNBS
{
public:
using NBSLevelsType = NBSLevels<executionSpace>;
using cellIterator = typename NBSLevelsType::cellIterator;
using IdType = typename NBSLevelsType::IdType;
using IndexType = typename NBSLevelsType::IndexType;
using Cells = typename NBSLevelsType::Cells;
using CellType = typename Cells::CellType;
using execution_space = typename NBSLevelsType::execution_space;
using memory_space = typename NBSLevelsType::memory_space;
protected:
real sizeRatio_ = 1.0;
int32 updateFrequency_= 1;
int32 currentIter_ = 0;
bool performedSearch_ = false;
NBSLevelsType NBSLevels_;
private:
bool performSearch()
{
if(currentIter_ % updateFrequency_ == 0)
{
currentIter_++;
return true;
}else
{
currentIter_++;
return false;
}
}
public:
TypeInfoNV("multiGridNBS");
multiGridNBS(
const dictionary& dict,
const box& domain,
real minSize,
real maxSize,
const ViewType1D<realx3, memory_space>& position,
const ViewType1D<real, memory_space>& diam)
:
sizeRatio_(
max(
dict.getVal<real>("sizeRatio"),
1.0
)),
updateFrequency_(
max(
dict.getVal<int32>("updateFrequency"),
1
)),
NBSLevels_(
domain,
minSize,
maxSize,
sizeRatio_,
position,
diam)
{}
INLINE_FUNCTION_HD
multiGridNBS(const multiGridNBS&) = default;
INLINE_FUNCTION_HD
multiGridNBS& operator = (const multiGridNBS&) = default;
INLINE_FUNCTION_HD
~multiGridNBS()=default;
//// - Methods
bool enterBoadSearch()const
{
return currentIter_%updateFrequency_==0;
}
bool performedSearch()const
{
return performedSearch_;
}
int32 numLevels()const
{
return NBSLevels_.numLevels();
}
auto getCellsLevels()const
{
Vector<Cells> cellsLvl("cells", numLevels(), numLevels(), RESERVE());
for(int32 lvl=0; lvl<numLevels(); lvl++)
{
cellsLvl[lvl] = NBSLevels_.getCells(lvl) ;
}
return cellsLvl;
}
auto getCells(int32 lvl)const
{
return NBSLevels_.getCells(lvl);
}
auto getCellIterator(int32 lvl)const
{
return NBSLevels_.getCellIterator(lvl);
}
bool objectSizeChanged(int32 newSize)
{
NBSLevels_.checkAllocateNext(newSize);
return true;
}
// - Perform the broad search to find pairs
// with force = true, perform broad search regardless of
// updateFrequency_ value
// on all the points in the range of [0,numPoints_)
template<typename PairsContainer>
bool broadSearch(PairsContainer& pairs, range activeRange, bool force=false)
{
if(force) currentIter_ = 0;
performedSearch_ = false;
if( !performSearch() ) return true;
NBSLevels_.build(activeRange);
NBSLevels_.findPairs(pairs);
performedSearch_ = true;
return true;
}
// - Perform the broad search to find pairs,
// ignore particles with incld(i) = true,
// with force = true, perform broad search regardless of
// updateFrequency_ value
template<typename PairsContainer, typename IncludeFunction>
bool broadSearch(PairsContainer& pairs, range activeRange, IncludeFunction incld, bool force = false)
{
if(force) currentIter_ = 0;
performedSearch_ = false;
if( !performSearch() ) return true;
NBSLevels_.build(activeRange, incld);
NBSLevels_.findPairs(pairs);
performedSearch_ = true;
return true;
}
};
}
#endif

78
src/Interaction/contactSearch/methods/particleWallContactSearchs.cpp Normal file
View File

@ -0,0 +1,78 @@
/*------------------------------- 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.
-----------------------------------------------------------------------------*/
template <typename method>
pFlow::particleWallContactSearchs<method>::particleWallContactSearchs
(
const dictionary &dict,
const box &domain,
real minSize,
real maxSize,
const ViewType1D<realx3, memory_space> &position,
const pFlagTypeDevice &flags,
const ViewType1D<real, memory_space> &diam
)
:
domainBox_(domain),
updateInterval_
(
dict.getValOrSet<uint32>("updateInterval", 1)
)
{
}
template <typename method>
bool pFlow::particleWallContactSearchs<method>::broadSearch
(
uint32 iter,
real t,
real dt,
csPairContainerType& ppPairs,
csPairContainerType& pwPairs,
const deviceViewType1D<realx3>& pointPos,
const pFlagTypeDevice& flags,
const deviceViewType1D<real>& diameter,
bool force
)
{
performedSearch_ = false;
if( !performSearch(iter, force) ) return true;
if(!getMethod().impl_broadSearch(
ppPairs,
pwPairs,
pointPos,
flags,
diameter))
{
fatalErrorInFunction<<
"Error in performing particle-particle broadSearch in method"<<
getMethod().typeName()<<endl;
return false;
}
lastUpdated_ = iter;
performedSearch_ = true;
return true;
}

123
src/Interaction/contactSearch/methods/particleWallContactSearchs.hpp Normal file
View File

@ -0,0 +1,123 @@
/*------------------------------- 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 __particleWallContactSearchs_hpp__
#define __particleWallContactSearchs_hpp__
#include "pointFlag.hpp"
#include "contactSearchGlobals.hpp"
#include "box.hpp"
namespace pFlow
{
template<typename method>
class particleWallContactSearchs
{
public:
using MethodType = method;
using IndexType = uint32;
using execution_space = DefaultExecutionSpace;
using memory_space = typename execution_space::memory_space;
private:
// friend
friend MethodType;
/// @brief box enclosing the simulation domain (local to processor)
box domainBox_;
/*/// @brief box enclosing the area for contact search (region with points)
box searchBox_;*/
/// @brief update interval in terms of iteration numebr
uint32 updateInterval_= 1;
/// @brief last iteration number which contact search has been performed
uint32 lastUpdated_ = 0;
/// @brief performed search?
bool performedSearch_ = false;
protected:
inline
auto& getMethod()
{
return static_cast<MethodType&>(*this);
}
inline
const auto& getMethod()const
{
return static_cast<const MethodType&>(*this);
}
public:
particleWallContactSearchs(
const dictionary& dict,
const box& domain,
real minSize,
real maxSize,
const ViewType1D<realx3, memory_space>& position,
const pFlagTypeDevice &flags,
const ViewType1D<real, memory_space>& diam
);
bool broadSearch
(
uint32 iter,
real t,
real dt,
csPairContainerType& ppPairs,
csPairContainerType& pwPairs,
const deviceViewType1D<realx3>& pointPos,
const pFlagTypeDevice& flags,
const deviceViewType1D<real>& diameter,
bool force = false
);
bool performedSearch()const
{
return performedSearch_;
}
bool performSearch(uint32 iter, bool force = false)const
{
if((iter-lastUpdated_) % updateInterval_ == 0 || iter == 0 || force )
{
return true;
}
return false;
}
};
} // pFlow
#include "particleWallContactSearchs.cpp"
#endif //__particleWallContactSearchs_hpp__

150
src/Interaction/contactSearch/wallMappings/cellMapping.hpp
View File

@ -1,150 +0,0 @@
/*------------------------------- 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 __cellMapping_hpp__
#define __cellMapping_hpp__
#include "cellsWallLevel0.hpp"
#include "dictionary.hpp"
namespace pFlow
{
template<
typename executionSpace
>
class cellMapping
{
public:
using cellsWallLevel0Type = cellsWallLevel0<executionSpace>;
using IdType = typename cellsWallLevel0Type::IdType;
using IndexType = typename cellsWallLevel0Type::IndexType;
using Cells = typename cellsWallLevel0Type::Cells;
using CellType = typename Cells::CellType;
using execution_space = typename cellsWallLevel0Type::execution_space;
using memory_space = typename cellsWallLevel0Type::memory_space;
using iBoxType = iBox<IndexType>;
protected:
// - update frequency
int32 updateFrequency_=1;
real cellExtent_;
int32 currentIter_ = 0;
/// a broad search has been occured during last pass?
bool performedSearch_ = false;
cellsWallLevel0Type cellsWallLevle_;
private:
bool performSearch()
{
if(currentIter_ % updateFrequency_ == 0)
{
currentIter_++;
return true;
}else
{
currentIter_++;
return false;
}
}
public:
TypeInfoNV("cellMapping");
cellMapping(
const dictionary& dict,
int32 numLevels,
const Vector<Cells>& ppCells,
int32 numPoints,
int32 numElements,
const ViewType1D<realx3,memory_space>& points,
const ViewType1D<int32x3,memory_space>& vertices
)
:
updateFrequency_(
max(
dict.getValOrSet<int32>(
"updateFrequency",
1),
1)),
cellExtent_(
max(
dict.getValOrSet<real>(
"cellExtent",
0.5),
0.5)),
cellsWallLevle_(
ppCells[0],
cellExtent_,
numPoints,
numElements,
points,
vertices
)
{}
bool enterBoadSearch()const
{
return currentIter_%updateFrequency_==0;
}
bool performedSearch()const
{
return performedSearch_;
}
template<typename PairsContainer, typename particleMapType>
bool broadSearch(PairsContainer& pairs, particleMapType& particleMap, bool force=false)
{
if(force) currentIter_ = 0;
performedSearch_= false;
if(!performSearch())return true;
cellsWallLevle_.broadSearch(pairs, particleMap);
performedSearch_ = true;
return true;
}
}; // cellMapping
} // pFlow
#endif

285
src/Interaction/contactSearch/wallMappings/cellsWallLevel0.hpp
View File

@ -1,285 +0,0 @@
/*------------------------------- 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 __cellsWallLevel0_hpp__
#define __cellsWallLevel0_hpp__
#include "types.hpp"
#include "KokkosTypes.hpp"
#include "cells.hpp"
#include "iBox.hpp"
namespace pFlow
{
template<
typename executionSpace
>
class cellsWallLevel0
:
public cells<int32>
{
public:
using IdType = int32;
using IndexType = int32;
using Cells = cells<IndexType>;
using CellType = typename Cells::CellType;
using execution_space = executionSpace;
using memory_space = typename execution_space::memory_space;
using iBoxType = iBox<IndexType>;
class TagFindCellRange2{};
protected:
// - box extent
real cellExtent_ = 0.5;
// - number of triangle elements
int32 numElements_ = 0;
// - number of points
int32 numPoints_ = 0;
// - ref to vectices (borrowed)
ViewType1D<int32x3, memory_space> vertices_;
// - ref to points in the trisurface (borrowed)
ViewType1D<realx3, memory_space> points_;
// cell range of element/triangle bounding box
ViewType1D<iBoxType, memory_space> elementBox_;
using tpPWContactSearch = Kokkos::TeamPolicy<
execution_space,
Kokkos::Schedule<Kokkos::Dynamic>,
Kokkos::IndexType<int32>
>;
using rpFindCellRange2Type =
Kokkos::RangePolicy<TagFindCellRange2, execution_space, Kokkos::IndexType<int32>>;
FUNCTION_H
void allocateArrays()
{
reallocNoInit( elementBox_, numElements_);
}
public:
TypeInfoNV("cellsWallLevel0");
INLINE_FUNCTION_HD
cellsWallLevel0(){}
FUNCTION_H
cellsWallLevel0(
const Cells& ppCells,
real cellExtent,
int32 numPoints,
int32 numElements,
const ViewType1D<realx3,memory_space>& points,
const ViewType1D<int32x3,memory_space>& vertices
)
:
Cells(ppCells),
cellExtent_( max(cellExtent, 0.5 ) ),
numElements_(numElements),
numPoints_(numPoints),
vertices_(vertices),
points_(points)
{
allocateArrays();
}
// - host call
// reset triangle elements if they have changed
bool resetElements(
int32 numElements,
int32 numPoints,
ViewType1D<realx3, memory_space>& points,
ViewType1D<int32x3, memory_space>& vertices )
{
numElements_ = numElements;
numPoints_ = numPoints;
points_ = points;
vertices_ = vertices;
allocateArrays();
return true;
}
INLINE_FUNCTION_HD
iBoxType elementBox(int32 i)const
{
return elementBox_[i];
}
INLINE_FUNCTION_HD
int32 numElements()const
{
return numElements_;
}
template<typename PairsContainer, typename particleMapType>
bool broadSearch(PairsContainer& pairs, particleMapType& particleMap)
{
// map walls onto the cells
this->build();
this->particleWallFindPairs(pairs, particleMap);
return true;
}
bool build()
{
Kokkos::parallel_for(
"cellsSimple::findcellrange2",
rpFindCellRange2Type(0,numElements_),
*this);
Kokkos::fence();
return true;
}
template<typename PairsContainer, typename particleMapType>
bool particleWallFindPairs(PairsContainer& pairs, particleMapType& particleMap)
{
int32 getFull = 1;
while (getFull)
{
getFull = findPairsElementRangeCount(pairs, particleMap.getCellIterator(0));
if(getFull)
{
// - resize the container
// note that getFull now shows the number of failed insertions.
uint32 len = max(getFull, 50);
auto oldCap = pairs.capacity();
pairs.increaseCapacityBy(len);
INFORMATION<<"Contact pair container capacity increased from "<<
oldCap << " to "
<< pairs.capacity() <<" in cellsWallLevel0."<<endINFO;
Kokkos::fence();
}
}
return true;
}
template<typename PairsContainer, typename CellIteratorType>
int32 findPairsElementRangeCount(PairsContainer& pairs, CellIteratorType cellIter)
{
int32 getFull =0;
const auto pwPairs = pairs;
const auto elementBox = elementBox_;
Kokkos::parallel_reduce(
"cellsSimple::findPairsElementRangeModified2",
tpPWContactSearch(numElements_, Kokkos::AUTO),
LAMBDA_HD(
const typename tpPWContactSearch::member_type & teamMember,
int32& valueToUpdate){
const int32 iTri = teamMember.league_rank();
const auto triBox = elementBox[iTri];
int32 getFull2 = 0;
auto bExtent = boxExtent(triBox);
int32 numCellBox = bExtent.x()*bExtent.y()*bExtent.z();
Kokkos::parallel_reduce(
Kokkos::TeamThreadRange( teamMember, numCellBox ),
[&] ( const int32 linIndex, int32 &innerUpdate )
{
CellType cell;
indexToCell(linIndex, triBox, cell);
int32 n = cellIter.start(cell.x(),cell.y(),cell.z());
while( n>-1)
{
// id is wall id the pair is (particle id, wall id)
if( pairs.insert(static_cast<IdType>(n), iTri) < 0 )
innerUpdate++;
n = cellIter.getNext(n);
}
},
getFull2
);
if ( teamMember.team_rank() == 0 ) valueToUpdate += getFull2;
},
getFull
);
return getFull;
}
INLINE_FUNCTION_HD
void operator()(TagFindCellRange2, int32 i) const
{
auto v = vertices_[i];
auto p1 = points_[v.x()];
auto p2 = points_[v.y()];
auto p3 = points_[v.z()];
realx3 minP, maxP;
this->extendBox(p1, p2, p3, cellExtent_, minP, maxP);
elementBox_[i] = iBoxType(this->pointIndex(minP), this->pointIndex(maxP));
}
}; // cellsWallLevel0
} // pFlow
#endif // __cellsWallLevel0_hpp__

152
src/Interaction/contactSearch/wallMappings/cellsWallLevels.hpp
View File

@ -1,152 +0,0 @@
/*------------------------------- 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 __cellsWallLevels_hpp__
#define __cellsWallLevels_hpp__
#include "cellsWallLevel0.hpp"
namespace pFlow
{
template<
typename executionSpace
>
class cellsWallLevels
{
public:
using cellsWallLevel0Type = cellsWallLevel0<executionSpace>;
using IdType = typename cellsWallLevel0Type::IdType;
using IndexType = typename cellsWallLevel0Type::IndexType;
using Cells = typename cellsWallLevel0Type::Cells;
using CellType = typename Cells::CellType;
using execution_space = typename cellsWallLevel0Type::execution_space;
using memory_space = typename cellsWallLevel0Type::memory_space;
using iBoxType = iBox<IndexType>;
protected:
int32 numLevles_=1;
Vector<cellsWallLevel0Type> cellsWallLevels_;
public:
TypeInfoNV("cellsWallLevels");
FUNCTION_H
cellsWallLevels(
int32 numLevels,
const Vector<Cells>& cellsLevels,
real cellExtent,
int32 numPoints,
int32 numElements,
const ViewType1D<realx3,memory_space>& points,
const ViewType1D<int32x3,memory_space>& vertices
)
:
numLevles_(numLevels),
cellsWallLevels_("cellsWallLevels",numLevels, numLevels, RESERVE())
{
for(int32 lvl=0; lvl<numLevles_; lvl++)
{
cellsWallLevels_[lvl] =
cellsWallLevel0Type(
cellsLevels[lvl],
cellExtent,
numPoints,
numElements,
points,
vertices);
}
}
template<typename PairsContainer, typename particleMapType>
bool broadSearch(PairsContainer& pairs, particleMapType& particleMap)
{
// map walls onto the cells
for(int32 lvl=0; lvl<numLevles_; lvl++)
{
cellsWallLevels_[lvl].build();
}
this->particleWallFindPairs(pairs, particleMap);
return true;
}
template<typename PairsContainer, typename particleMapType>
bool particleWallFindPairs(PairsContainer& pairs, particleMapType& particleMap)
{
int32 getFull = 1;
while (getFull)
{
getFull = 0;
for(int32 lvl=0; lvl<numLevles_; lvl++)
{
getFull +=
cellsWallLevels_[lvl].findPairsElementRangeCount(
pairs,
particleMap.getCellIterator(lvl));
}
if(getFull)
{
// - resize the container
// note that getFull now shows the number of failed insertions.
uint32 len = max(getFull, 5000);
auto oldCap = pairs.capacity();
pairs.increaseCapacityBy(len);
INFORMATION<<"Contact pair container capacity increased from "<<
oldCap << " to "
<< pairs.capacity() <<" in cellsWallLevels."<<endINFO;
Kokkos::fence();
}
}
return true;
}
}; // cellsWallLevels
} // pFlow
#endif // __cellsWallLevels_hpp__

153
src/Interaction/contactSearch/wallMappings/multiGridMapping.hpp
View File

@ -1,153 +0,0 @@
/*------------------------------- 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 __multiGridMapping_hpp__
#define __multiGridMapping_hpp__
#include "cellsWallLevels.hpp"
#include "dictionary.hpp"
namespace pFlow
{
template<
typename executionSpace
>
class multiGridMapping
{
public:
using CellsWallLevelType = cellsWallLevels<executionSpace>;
using IdType = typename CellsWallLevelType::IdType;
using IndexType = typename CellsWallLevelType::IndexType;
using Cells = typename CellsWallLevelType::Cells;
using CellType = typename Cells::CellType;
using execution_space = typename CellsWallLevelType::execution_space;
using memory_space = typename CellsWallLevelType::memory_space;
using iBoxType = iBox<IndexType>;
protected:
// - update frequency
int32 updateFrequency_=1;
real cellExtent_;
int32 currentIter_ = 0;
/// a broad search has been occured during last pass?
bool performedSearch_ = false;
CellsWallLevelType cellsWallLevle_;
private:
bool performSearch()
{
if(currentIter_ % updateFrequency_ == 0)
{
currentIter_++;
return true;
}else
{
currentIter_++;
return false;
}
}
public:
TypeInfoNV("multiGridMapping");
multiGridMapping(
const dictionary& dict,
int32 numLevels,
const Vector<Cells>& ppCells,
int32 numPoints,
int32 numElements,
const ViewType1D<realx3,memory_space>& points,
const ViewType1D<int32x3,memory_space>& vertices
)
:
updateFrequency_(
max(
dict.getVal<int32>("updateFrequency"),
1)),
cellExtent_(
max(
dict.getVal<real>("cellExtent"),
0.5)),
cellsWallLevle_(
numLevels,
ppCells,
cellExtent_,
numPoints,
numElements,
points,
vertices
)
{
REPORT(3)<<"Multi-grid wall mapping with "<<
yellowText(numLevels)<<" levels has been created."<<endREPORT;
}
bool enterBoadSearch()const
{
return currentIter_%updateFrequency_==0;
}
bool performedSearch()const
{
return performedSearch_;
}
template<typename PairsContainer, typename particleMapType>
bool broadSearch(PairsContainer& pairs, particleMapType& particleMap, bool force=false)
{
if(force) currentIter_ = 0;
performedSearch_= false;
if(!performSearch())return true;
cellsWallLevle_.broadSearch(pairs, particleMap);
performedSearch_ = true;
return true;
}
}; // multiGridMapping
} // pFlow
#endif

36
src/Interaction/interaction/interaction.cpp
View File

@ -19,7 +19,9 @@ Licence:
-----------------------------------------------------------------------------*/
#include "interaction.hpp"
#include "particles.hpp"
#include "vocabs.hpp"
#include "geometry.hpp"
pFlow::interaction::interaction
(
@ -28,26 +30,16 @@ pFlow::interaction::interaction
const geometry& geom
)
:
demInteraction(control, control.caseSetup().path()+interactionFile__),
interactionBase(prtcl, geom),
fileDict_(control.caseSetup().emplaceObject<dictionary>(
objectFile(
interactionFile__,
"",
objectFile::READ_ALWAYS,
objectFile::WRITE_NEVER),
interactionFile__,
true ))
property(interactionFile__, &control.caseSetup()),
observer
(
&prtcl.dynPointStruct(),
msg_
),
demComponent("interaction", control),
particles_(prtcl),
geometry_(geom)
{
this->subscribe(prtcl.pStruct());
contactSearch_ = contactSearch::create(
fileDict_.subDict("contactSearch"),
this->control().domain(),
prtcl,
geom,
timers()
);
}
@ -79,14 +71,14 @@ pFlow::uniquePtr<pFlow::interaction> pFlow::interaction::create
clType);
REPORT(1)<< "Selecting interaction model..."<<endREPORT;
REPORT(1)<< "Selecting interaction model..."<<END_REPORT;
if( systemControlvCtorSelector_.search(interactionModel) )
{
auto objPtr =
systemControlvCtorSelector_[interactionModel]
(control, prtcl, geom);
REPORT(2)<<"Model "<<greenText(interactionModel)<<" is created."<<endREPORT;
REPORT(2)<<"Model "<<Green_Text(interactionModel)<<" is created."<<END_REPORT;
return objPtr;
}
else

59
src/Interaction/interaction/interaction.hpp
View File

@ -21,9 +21,10 @@ Licence:
#ifndef __interaction_hpp__
#define __interaction_hpp__
#include "demInteraction.hpp"
#include "eventObserver.hpp"
#include "interactionBase.hpp"
#include "demComponent.hpp"
#include "property.hpp"
#include "observer.hpp"
#include "systemControl.hpp"
#include "contactSearch.hpp"
@ -31,27 +32,26 @@ Licence:
namespace pFlow
{
class particles;
class geometry;
class interaction
:
public demInteraction,
public eventObserver,
public interactionBase
{
public:
public property,
public observer,
public demComponent
using IdType = typename interactionBase::IdType;
{
private:
using IndexType = typename interactionBase::IndexType;
/// reference to particles object
const particles& particles_;
using ExecutionSpace = typename interactionBase::ExecutionSpace;
/// reference to geometry object
const geometry& geometry_;
protected:
/// interaction file dictionary
dictionary& fileDict_;
/// contact search object for pp and pw interactions
uniquePtr<contactSearch> contactSearch_ = nullptr;
static inline
const message msg_ = message::ITEM_REARRANGE;
public:
@ -65,11 +65,10 @@ public:
const particles& prtcl,
const geometry& geom );
~interaction() override = default;
virtual ~interaction() = default;
create_vCtor(
create_vCtor
(
interaction,
systemControl,
(
@ -78,22 +77,8 @@ public:
const geometry& geom
),
(control, prtcl, geom)
);
);
auto& contactSearchPtr()
{
return contactSearch_;
}
auto& contactSearchRef()
{
return contactSearch_();
}
const auto& fileDict()const
{
return fileDict_;
}
static
uniquePtr<interaction> create(

201
src/Interaction/sphereInteraction/sphereInteraction.cpp
View File

@ -18,22 +18,18 @@ Licence:
-----------------------------------------------------------------------------*/
template<
typename contactForceModel,
typename geometryMotionModel,
template <class, class, class> class contactListType >
bool pFlow::sphereInteraction<contactForceModel,geometryMotionModel, contactListType>::
createSphereInteraction()
template<typename cFM,typename gMM,template <class, class, class> class cLT>
bool pFlow::sphereInteraction<cFM,gMM, cLT>::createSphereInteraction()
{
realVector_D rhoD(this->densities());
realVector_D rhoD("densities", this->densities());
auto modelDict = this->fileDict().subDict("model");
auto modelDict = this->subDict("model");
REPORT(1)<<"Createing contact force model . . ."<<endREPORT;
REPORT(1)<<"Createing contact force model . . ."<<END_REPORT;
forceModel_ = makeUnique<ContactForceModel>(
this->numMaterials(),
rhoD.deviceVector(),
rhoD.deviceView(),
modelDict );
@ -47,17 +43,9 @@ bool pFlow::sphereInteraction<contactForceModel,geometryMotionModel, contactList
}
template<
typename contactForceModel,
typename geometryMotionModel,
template <class, class, class> class contactListType >
bool pFlow::sphereInteraction<contactForceModel,geometryMotionModel, contactListType>::
sphereSphereInteraction()
template<typename cFM,typename gMM,template <class, class, class> class cLT>
bool pFlow::sphereInteraction<cFM,gMM, cLT>::sphereSphereInteraction()
{
auto lastItem = ppContactList_().loopCount();
// create the kernel functor
@ -66,17 +54,17 @@ bool pFlow::sphereInteraction<contactForceModel,geometryMotionModel, contactList
this->dt(),
this->forceModel_(),
ppContactList_(), // to be read
sphParticles_.diameter().deviceVectorAll(),
sphParticles_.propertyId().deviceVectorAll(),
sphParticles_.pointPosition().deviceVectorAll(),
sphParticles_.velocity().deviceVectorAll(),
sphParticles_.rVelocity().deviceVectorAll(),
sphParticles_.contactForce().deviceVectorAll(),
sphParticles_.contactTorque().deviceVectorAll()
sphParticles_.diameter().deviceViewAll(),
sphParticles_.propertyId().deviceViewAll(),
sphParticles_.pointPosition().deviceViewAll(),
sphParticles_.velocity().deviceViewAll(),
sphParticles_.rVelocity().deviceViewAll(),
sphParticles_.contactForce().deviceViewAll(),
sphParticles_.contactTorque().deviceViewAll()
);
Kokkos::parallel_for(
"",
"ppInteraction",
rpPPInteraction(0,lastItem),
ppInteraction
);
@ -87,34 +75,32 @@ bool pFlow::sphereInteraction<contactForceModel,geometryMotionModel, contactList
}
template<
typename contactForceModel,
typename geometryMotionModel,
template <class, class, class> class contactListType >
bool pFlow::sphereInteraction<contactForceModel,geometryMotionModel, contactListType>::
sphereWallInteraction()
template<typename cFM,typename gMM,template <class, class, class> class cLT>
bool pFlow::sphereInteraction<cFM,gMM, cLT>::sphereWallInteraction()
{
int32 lastItem = pwContactList_().loopCount();
uint32 lastItem = pwContactList_().loopCount();
uint32 iter = this->currentIter();
real t = this->currentTime();
real dt = this->dt();
pFlow::sphereInteractionKernels::pwInteractionFunctor
pwInteraction(
this->dt(),
dt,
this->forceModel_(),
pwContactList_(),
geometryMotion_.getTriangleAccessor(),
geometryMotion_.getModel(t) ,
sphParticles_.diameter().deviceVectorAll() ,
sphParticles_.propertyId().deviceVectorAll(),
sphParticles_.pointPosition().deviceVectorAll(),
sphParticles_.velocity().deviceVectorAll(),
sphParticles_.rVelocity().deviceVectorAll() ,
sphParticles_.contactForce().deviceVectorAll(),
sphParticles_.contactTorque().deviceVectorAll() ,
geometryMotion_.triMotionIndex().deviceVectorAll(),
geometryMotion_.propertyId().deviceVectorAll(),
geometryMotion_.contactForceWall().deviceVectorAll()
geometryMotion_.getModel(iter, t, dt) ,
sphParticles_.diameter().deviceViewAll() ,
sphParticles_.propertyId().deviceViewAll(),
sphParticles_.pointPosition().deviceViewAll(),
sphParticles_.velocity().deviceViewAll(),
sphParticles_.rVelocity().deviceViewAll() ,
sphParticles_.contactForce().deviceViewAll(),
sphParticles_.contactTorque().deviceViewAll() ,
geometryMotion_.triMotionIndex().deviceViewAll(),
geometryMotion_.propertyId().deviceViewAll(),
geometryMotion_.contactForceWall().deviceViewAll()
);
Kokkos::parallel_for(
@ -127,4 +113,123 @@ bool pFlow::sphereInteraction<contactForceModel,geometryMotionModel, contactList
Kokkos::fence();
return true;
}
}
template<typename cFM,typename gMM,template <class, class, class> class cLT>
pFlow::sphereInteraction<cFM,gMM, cLT>::sphereInteraction
(
systemControl& control,
const particles& prtcl,
const geometry& geom
)
:
interaction(control, prtcl, geom),
geometryMotion_(dynamic_cast<const GeometryMotionModel&>(geom)),
sphParticles_(dynamic_cast<const sphereParticles&>(prtcl)),
ppInteractionTimer_("sphere-sphere interaction", &this->timers()),
pwInteractionTimer_("sphere-wall interaction", &this->timers()),
contactListTimer_("contact list management", &this->timers()),
contactListTimer0_("contact list clear", &this->timers())
{
contactSearch_ = contactSearch::create(
subDict("contactSearch"),
prtcl.thisDomain().domainBox(),
prtcl,
geom,
timers());
if(!createSphereInteraction())
{
fatalExit;
}
}
template<typename cFM,typename gMM,template <class, class, class> class cLT>
bool pFlow::sphereInteraction<cFM,gMM, cLT>::beforeIteration()
{
return true;
}
template<typename cFM,typename gMM,template <class, class, class> class cLT>
bool pFlow::sphereInteraction<cFM,gMM, cLT>::iterate()
{
auto iter = this->currentIter();
auto t = this->currentTime();
auto dt = this->dt();
//output<<"iter, t, dt "<< iter<<" "<< t << " "<<dt<<endl;
bool broadSearch = contactSearch_().enterBroadSearch(iter, t, dt);
if(broadSearch)
{
contactListTimer0_.start();
ppContactList_().beforeBroadSearch();
pwContactList_().beforeBroadSearch();
contactListTimer0_.end();
}
if( sphParticles_.numActive()<=0)return true;
if( !contactSearch_().broadSearch(
iter,
t,
dt,
ppContactList_(),
pwContactList_()) )
{
fatalErrorInFunction<<
"unable to perform broadSearch.\n";
fatalExit;
}
if(broadSearch && contactSearch_().performedBroadSearch())
{
contactListTimer_.start();
ppContactList_().afterBroadSearch();
pwContactList_().afterBroadSearch();
contactListTimer_.end();
}
ppInteractionTimer_.start();
sphereSphereInteraction();
ppInteractionTimer_.end();
pwInteractionTimer_.start();
sphereWallInteraction();
pwInteractionTimer_.end();
return true;
}
template<typename cFM,typename gMM,template <class, class, class> class cLT>
bool pFlow::sphereInteraction<cFM,gMM, cLT>::afterIteration()
{
return true;
}
template<typename cFM,typename gMM,template <class, class, class> class cLT>
bool pFlow::sphereInteraction<cFM,gMM, cLT>::hearChanges
(
real t,
real dt,
uint32 iter,
const message& msg,
const anyList& varList
)
{
if(msg.equivalentTo(message::ITEM_REARRANGE))
{
notImplementedFunction;
}
return true;
}

167
src/Interaction/sphereInteraction/sphereInteraction.hpp
View File

@ -23,7 +23,6 @@ Licence:
#include "interaction.hpp"
#include "sphereParticles.hpp"
#include "sphereInteractionKernels.hpp"
namespace pFlow
@ -43,22 +42,19 @@ public:
using ContactForceModel = contactForceModel;
using MotionModel = typename geometryMotionModel::MotionModel;
using MotionModel = typename geometryMotionModel::MotionModel;
using ModelStorage = typename ContactForceModel::contactForceStorage;
using ModelStorage = typename ContactForceModel::contactForceStorage;
using IdType = typename interaction::IdType;
using IdType = uint32;
using IndexType = typename interaction::IndexType;
using ExecutionSpace = typename interaction::ExecutionSpace;
using IndexType = uint32;
using ContactListType =
contactListType<ModelStorage, ExecutionSpace, IdType>;
contactListType<ModelStorage, DefaultExecutionSpace, IdType>;
using PairsContainerType= typename contactSearch::PairContainerType;
protected:
private:
/// const reference to geometry
const GeometryMotionModel& geometryMotion_;
@ -66,15 +62,17 @@ protected:
/// const reference to particles
const sphereParticles& sphParticles_;
/// contact search object for pp and pw interactions
uniquePtr<contactSearch> contactSearch_ = nullptr;
/// contact force model
uniquePtr<ContactForceModel> forceModel_ = nullptr;
uniquePtr<ContactForceModel> forceModel_ = nullptr;
/// contact list for particle-particle interactoins (keeps the history)
uniquePtr<ContactListType> ppContactList_ = nullptr;
uniquePtr<ContactListType> ppContactList_ = nullptr;
/// contact list for particle-wall interactions (keeps the history)
uniquePtr<ContactListType> pwContactList_ = nullptr;
uniquePtr<ContactListType> pwContactList_ = nullptr;
/// timer for particle-particle interaction computations
Timer ppInteractionTimer_;
@ -82,42 +80,39 @@ protected:
/// timer for particle-wall interaction computations
Timer pwInteractionTimer_;
Timer contactListTimer_;
Timer contactListTimer0_;
bool createSphereInteraction();
bool managePPContactLists();
bool sphereSphereInteraction();
bool managePWContactLists();
bool sphereWallInteraction();
//bool managePPContactLists();
//bool managePWContactLists();
/// range policy for p-p interaction execution
using rpPPInteraction =
Kokkos::RangePolicy<Kokkos::IndexType<int32>, Kokkos::Schedule<Kokkos::Dynamic>>;
Kokkos::RangePolicy<Kokkos::IndexType<uint32>, Kokkos::Schedule<Kokkos::Dynamic>>;
/// range policy for p-w interaction execution
using rpPWInteraction = rpPPInteraction;
public:
TypeInfoTemplate3("sphereInteraction", ContactForceModel, MotionModel, ContactListType);
// constructor
TypeInfoTemplate13("sphereInteraction", ContactForceModel, MotionModel, ContactListType);
/// Constructor from components
sphereInteraction(
systemControl& control,
const particles& prtcl,
const geometry& geom)
:
interaction(control, prtcl, geom),
geometryMotion_(dynamic_cast<const GeometryMotionModel&>(geom)),
sphParticles_(dynamic_cast<const sphereParticles&>(prtcl)),
ppInteractionTimer_("sphere-sphere interaction", &this->timers()),
pwInteractionTimer_("sphere-wall interaction", &this->timers())
{
if(!createSphereInteraction())
{
fatalExit;
}
}
const geometry& geom);
/// Add virtual constructor
add_vCtor
(
interaction,
@ -125,97 +120,25 @@ public:
systemControl
);
/// This is called in time loop, before iterate. (overriden from demComponent)
bool beforeIteration() override;
/// This is called in time loop. Perform the main calculations
/// when the component should evolve along time. (overriden from demComponent)
bool iterate() override;
/// This is called in time loop, after iterate. (overriden from demComponent)
bool afterIteration() override;
bool beforeIteration() override
{
return true;
}
bool iterate() override
{
//Info<<"before contact search"<<endInfo;
////Info<<"interaction iterrate start"<<endInfo;
if(this->contactSearch_)
{
if( this->contactSearch_().ppEnterBroadSearch())
{
//Info<<" before ppEnterBroadSearch"<<endInfo;
ppContactList_().beforeBroadSearch();
//Info<<" after ppEnterBroadSearch"<<endInfo;
}
if(this->contactSearch_().pwEnterBroadSearch())
{
//Info<<" before pwEnterBroadSearch"<<endInfo;
pwContactList_().beforeBroadSearch();
//Info<<" after pwEnterBroadSearch"<<endInfo;
}
//Info<<" before broadSearch"<<endInfo;
if( !contactSearch_().broadSearch(
ppContactList_(),
pwContactList_()) )
{
fatalErrorInFunction<<
"unable to perform broadSearch.\n";
fatalExit;
}
//Info<<" before broadSearch"<<endInfo;
if(this->contactSearch_().ppPerformedBroadSearch())
{
//Info<<" before afterBroadSearch"<<endInfo;
ppContactList_().afterBroadSearch();
//Info<<" after afterBroadSearch"<<endInfo;
}
if(this->contactSearch_().pwPerformedBroadSearch())
{
//Info<<" before pwContactList_().afterBroadSearch()"<<endInfo;
pwContactList_().afterBroadSearch();
//Info<<" after pwContactList_().afterBroadSearch()"<<endInfo;
}
}
//Info<<"after contact search"<<endInfo;
if( sphParticles_.numActive()<=0)return true;
//Info<<"before sphereSphereInteraction "<<endInfo;
ppInteractionTimer_.start();
sphereSphereInteraction();
ppInteractionTimer_.end();
//Info<<"after sphereSphereInteraction "<<endInfo;
//Info<<"before sphereWallInteraction "<<endInfo;
pwInteractionTimer_.start();
sphereWallInteraction();
pwInteractionTimer_.end();
//Info<<"after sphereWallInteraction "<<endInfo;
return true;
}
bool afterIteration() override
{
return true;
}
bool update(const eventMessage& msg)override
{
// it requires not action regarding any changes in the
// point structure
return true;
}
bool sphereSphereInteraction();
bool sphereWallInteraction();
/// Check for changes in the point structures. (overriden from observer)
bool hearChanges(
real t,
real dt,
uint32 iter,
const message& msg,
const anyList& varList)override;
};

18
src/Interaction/sphereInteraction/sphereInteractionKernels.hpp
View File

@ -42,7 +42,7 @@ struct ppInteractionFunctor
ContactListType tobeFilled_;
deviceViewType1D<real> diam_;
deviceViewType1D<int8> propId_;
deviceViewType1D<uint32> propId_;
deviceViewType1D<realx3> pos_;
deviceViewType1D<realx3> lVel_;
deviceViewType1D<realx3> rVel_;
@ -55,7 +55,7 @@ struct ppInteractionFunctor
ContactForceModel forceModel,
ContactListType tobeFilled,
deviceViewType1D<real> diam,
deviceViewType1D<int8> propId,
deviceViewType1D<uint32> propId,
deviceViewType1D<realx3> pos,
deviceViewType1D<realx3> lVel,
deviceViewType1D<realx3> rVel,
@ -75,7 +75,7 @@ struct ppInteractionFunctor
{}
INLINE_FUNCTION_HD
void operator()(const int32 n)const
void operator()(const uint32 n)const
{
if(!tobeFilled_.isValid(n))return;
@ -181,14 +181,14 @@ struct pwInteractionFunctor
MotionModel motionModel_;
deviceViewType1D<real> diam_;
deviceViewType1D<int8> propId_;
deviceViewType1D<uint32> propId_;
deviceViewType1D<realx3> pos_;
deviceViewType1D<realx3> lVel_;
deviceViewType1D<realx3> rVel_;
deviceViewType1D<realx3> cForce_;
deviceViewType1D<realx3> cTorque_;
deviceViewType1D<int8> wTriMotionIndex_;
deviceViewType1D<int8> wPropId_;
deviceViewType1D<uint32> wTriMotionIndex_;
deviceViewType1D<uint32> wPropId_;
deviceViewType1D<realx3> wCForce_;
@ -199,14 +199,14 @@ struct pwInteractionFunctor
TraingleAccessor triangles,
MotionModel motionModel ,
deviceViewType1D<real> diam ,
deviceViewType1D<int8> propId,
deviceViewType1D<uint32> propId,
deviceViewType1D<realx3> pos ,
deviceViewType1D<realx3> lVel,
deviceViewType1D<realx3> rVel,
deviceViewType1D<realx3> cForce,
deviceViewType1D<realx3> cTorque ,
deviceViewType1D<int8> wTriMotionIndex,
deviceViewType1D<int8> wPropId,
deviceViewType1D<uint32> wTriMotionIndex,
deviceViewType1D<uint32> wPropId,
deviceViewType1D<realx3> wCForce)
:
dt_(dt),

35
src/Interaction/sphereInteraction/sphereInteractions.cpp
View File

@ -28,22 +28,22 @@ Licence:
template class pFlow::sphereInteraction<
pFlow::cfModels::limitedLinearNormalRolling,
pFlow::fixedGeometry,
pFlow::stationaryGeometry,
pFlow::unsortedContactList>;
template class pFlow::sphereInteraction<
pFlow::cfModels::limitedLinearNormalRolling,
pFlow::fixedGeometry,
pFlow::stationaryGeometry,
pFlow::sortedContactList>;
template class pFlow::sphereInteraction<
pFlow::cfModels::nonLimitedLinearNormalRolling,
pFlow::fixedGeometry,
pFlow::stationaryGeometry,
pFlow::unsortedContactList>;
template class pFlow::sphereInteraction<
pFlow::cfModels::nonLimitedLinearNormalRolling,
pFlow::fixedGeometry,
pFlow::stationaryGeometry,
pFlow::sortedContactList>;
template class pFlow::sphereInteraction<
@ -88,7 +88,7 @@ template class pFlow::sphereInteraction<
pFlow::sortedContactList>;
template class pFlow::sphereInteraction<
/*template class pFlow::sphereInteraction<
pFlow::cfModels::limitedLinearNormalRolling,
pFlow::multiRotationAxisMotionGeometry,
pFlow::unsortedContactList>;
@ -106,28 +106,28 @@ template class pFlow::sphereInteraction<
template class pFlow::sphereInteraction<
pFlow::cfModels::nonLimitedLinearNormalRolling,
pFlow::multiRotationAxisMotionGeometry,
pFlow::sortedContactList>;
pFlow::sortedContactList>;*/
/// non-linear models
template class pFlow::sphereInteraction<
pFlow::cfModels::limitedNonLinearNormalRolling,
pFlow::fixedGeometry,
pFlow::stationaryGeometry,
pFlow::unsortedContactList>;
template class pFlow::sphereInteraction<
pFlow::cfModels::limitedNonLinearNormalRolling,
pFlow::fixedGeometry,
pFlow::stationaryGeometry,
pFlow::sortedContactList>;
template class pFlow::sphereInteraction<
pFlow::cfModels::nonLimitedNonLinearNormalRolling,
pFlow::fixedGeometry,
pFlow::stationaryGeometry,
pFlow::unsortedContactList>;
template class pFlow::sphereInteraction<
pFlow::cfModels::nonLimitedNonLinearNormalRolling,
pFlow::fixedGeometry,
pFlow::stationaryGeometry,
pFlow::sortedContactList>;
template class pFlow::sphereInteraction<
@ -171,7 +171,7 @@ template class pFlow::sphereInteraction<
pFlow::vibratingMotionGeometry,
pFlow::sortedContactList>;
template class pFlow::sphereInteraction<
/*template class pFlow::sphereInteraction<
pFlow::cfModels::limitedNonLinearNormalRolling,
pFlow::multiRotationAxisMotionGeometry,
pFlow::unsortedContactList>;
@ -189,28 +189,28 @@ template class pFlow::sphereInteraction<
template class pFlow::sphereInteraction<
pFlow::cfModels::nonLimitedNonLinearNormalRolling,
pFlow::multiRotationAxisMotionGeometry,
pFlow::sortedContactList>;
pFlow::sortedContactList>;*/
// - nonLinearMod models
template class pFlow::sphereInteraction<
pFlow::cfModels::limitedNonLinearModNormalRolling,
pFlow::fixedGeometry,
pFlow::stationaryGeometry,
pFlow::unsortedContactList>;
template class pFlow::sphereInteraction<
pFlow::cfModels::limitedNonLinearModNormalRolling,
pFlow::fixedGeometry,
pFlow::stationaryGeometry,
pFlow::sortedContactList>;
template class pFlow::sphereInteraction<
pFlow::cfModels::nonLimitedNonLinearModNormalRolling,
pFlow::fixedGeometry,
pFlow::stationaryGeometry,
pFlow::unsortedContactList>;
template class pFlow::sphereInteraction<
pFlow::cfModels::nonLimitedNonLinearModNormalRolling,
pFlow::fixedGeometry,
pFlow::stationaryGeometry,
pFlow::sortedContactList>;
template class pFlow::sphereInteraction<
@ -255,7 +255,7 @@ template class pFlow::sphereInteraction<
pFlow::sortedContactList>;
template class pFlow::sphereInteraction<
/*template class pFlow::sphereInteraction<
pFlow::cfModels::limitedNonLinearModNormalRolling,
pFlow::multiRotationAxisMotionGeometry,
pFlow::unsortedContactList>;
@ -274,3 +274,4 @@ template class pFlow::sphereInteraction<
pFlow::cfModels::nonLimitedNonLinearModNormalRolling,
pFlow::multiRotationAxisMotionGeometry,
pFlow::sortedContactList>;
*/

16
src/MotionModel/CMakeLists.txt
View File

@ -1,13 +1,19 @@
list(APPEND SourceFiles
entities/rotatingAxis/rotatingAxis.cpp
entities/multiRotatingAxis/multiRotatingAxis.cpp
entities/timeInterval/timeInterval.cpp
entities/vibrating/vibrating.cpp
fixedWall/fixedWall.cpp
entities/rotatingAxis/rotatingAxis.cpp
rotatingAxisMotion/rotatingAxisMotion.cpp
multiRotatingAxisMotion/multiRotatingAxisMotion.cpp
entities/vibrating/vibrating.cpp
vibratingMotion/vibratingMotion.cpp
stationaryWall/stationaryWall.cpp
entities/stationary/stationary.cpp
#entities/multiRotatingAxis/multiRotatingAxis.cpp
#multiRotatingAxisMotion/multiRotatingAxisMotion.cpp
)
set(link_libs Kokkos::kokkos phasicFlow)

147
src/MotionModel/MotionModel/MotionModel.cpp Normal file
View File

@ -0,0 +1,147 @@
/*------------------------------- 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.
-----------------------------------------------------------------------------*/
template<typename Model, typename Component>
inline
bool pFlow::MotionModel<Model, Component>::impl_nameToIndex(const word& name, uint32& indx)const
{
if( auto i = componentNames_.findi(name); i == -1)
{
return false;
}
else
{
indx = static_cast<uint32>(i);
return true;
}
}
template<typename Model, typename Component>
inline
bool pFlow::MotionModel<Model, Component>::impl_indexToName(uint32 i, word& name)const
{
if(i < numComponents_ )
{
name = componentNames_[i];
return true;
}
else
{
return false;
}
}
template<typename Model, typename Component>
inline
bool pFlow::MotionModel<Model, Component>::impl_readDictionary
(
const dictionary& dict
)
{
auto modelName = dict.getVal<word>("motionModel");
if(modelName != getTypeName<ModelComponent>())
{
fatalErrorInFunction<<
" motionModel should be "<< Yellow_Text(getTypeName<ModelComponent>())<<
", but found "<< Yellow_Text(modelName)<<endl;
return false;
}
auto& motionInfo = dict.subDict(modelName+"Info");
auto compNames = motionInfo.dictionaryKeywords();
Vector<ModelComponent> components(
"Read::modelComponent",
compNames.size()+1,
0,
RESERVE());
componentNames_.clear();
for(auto& cName: compNames)
{
auto& compDict = motionInfo.subDict(cName);
if(auto compPtr = makeUnique<ModelComponent>(compDict); compPtr)
{
components.push_back(compPtr());
componentNames_.push_back(cName);
}
}
if( !componentNames_.search("none") )
{
components.push_back
(
impl_noneComponent()
);
componentNames_.push_back("none");
}
motionComponents_.assign(components);
numComponents_ = motionComponents_.size();
return true;
}
template<typename Model, typename Component>
inline
bool pFlow::MotionModel<Model, Component>::impl_writeDictionary
(
dictionary& dict
)const
{
word modelName = getTypeName<ModelComponent>();
dict.add("motionModel", modelName );
auto modelDictName = modelName+"Info";
auto& motionInfo = dict.subDictOrCreate(modelDictName);
auto hostComponents = motionComponents_.hostView();
ForAll(i, motionComponents_)
{
auto& axDict = motionInfo.subDictOrCreate(componentNames_[i]);
if( !hostComponents[i].write(axDict))
{
fatalErrorInFunction<<
" error in writing axis "<< componentNames_[i] << " to dicrionary "
<< motionInfo.globalName()<<endl;
return false;
}
}
return true;
}
template<typename Model, typename Component>
pFlow::MotionModel<Model, Component>::MotionModel()
:
motionComponents_("motionComponents")
{}

246
src/MotionModel/MotionModel/MotionModel.hpp Normal file
View File

@ -0,0 +1,246 @@
/*------------------------------- 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 __MotionModel_hpp__
#define __MotionModel_hpp__
#include "VectorSingles.hpp"
#include "Vector.hpp"
#include "List.hpp"
#include "dictionary.hpp"
namespace pFlow
{
/**
* @brief Motion model abstract class (CRTP) for all the motion models
*/
template<typename Model, typename Component>
class MotionModel
{
public:
using ModelType = Model;
using ModelComponent = Component;
using ComponentVector_D = VectorSingle<ModelComponent>;
/** Motion model class to be passed to computational units/kernels for
* transfing points and returning velocities at various positions
*/
class ModelInterface
{
private:
deviceViewType1D<ModelComponent> components_;
uint32 numComponents_=0;
public:
INLINE_FUNCTION_HD
ModelInterface
(
deviceViewType1D<ModelComponent> Comps,
uint32 numComp
):
components_(Comps),
numComponents_(numComp)
{}
INLINE_FUNCTION_HD
ModelInterface(const ModelInterface&) = default;
INLINE_FUNCTION_HD
ModelInterface& operator=(const ModelInterface&) = default;
INLINE_FUNCTION_HD
ModelInterface(ModelInterface&&) noexcept = default;
INLINE_FUNCTION_HD
ModelInterface& operator=(ModelInterface&&) noexcept= default;
INLINE_FUNCTION_HD
~ModelInterface()=default;
INLINE_FUNCTION_HD
realx3 pointVelocity(uint32 n, const realx3& p)const
{
return components_[n].linVelocityPoint(p);
}
INLINE_FUNCTION_HD
realx3 operator()(uint32 n, const realx3& p)const
{
return pointVelocity(n,p);
}
INLINE_FUNCTION_HD
realx3 transferPoint(uint32 n, const realx3 p, real dt)const
{
return components_[n].transferPoint(p, dt);
}
INLINE_FUNCTION_HD
uint32 size()const
{
return numComponents_;
}
};
private:
// friends
friend ModelType;
/// Number of axes components
uint32 numComponents_= 0;
/// Vector to store motion components
ComponentVector_D motionComponents_;
/// Names of compoenents
wordList componentNames_;
protected:
/// @brief obtain a reference to the actual motion model
inline
auto& getModel()
{
return static_cast<ModelType&>(*this);
}
/// @brief Obtain a const reference to the actual motion model
inline
const auto& getModel()const
{
return static_cast<const ModelType&>(*this);
}
/// Return a none object for the motion model
inline
auto impl_noneComponent()const
{
return ModelType::noneComponent();
}
/// Name of the compoenent to index of the component
bool impl_nameToIndex(const word& name, uint32& idx)const;
/// Component index to motion component name
bool impl_indexToName(uint32 i, word& name)const;
/// const reference to the list of component names
inline
const wordList& impl_componentNames()const
{
return componentNames_;
}
/// Return model interface
auto impl_getModelInterface(uint32 iter, real t, real dt)const
{
getModel().impl_setTime(iter, t, dt);
return ModelInterface(
motionComponents_.deviceViewAll(),
numComponents_);
}
/// Read from dictionary
bool impl_readDictionary(const dictionary& dict);
/// Write to dictionary
bool impl_writeDictionary(dictionary& dict)const;
public:
// - Constructors
/// Empty
MotionModel();
/// Copy constructor
MotionModel(const MotionModel&) = default;
/// Move constructor
MotionModel(MotionModel&&) = default;
/// Copy assignment
MotionModel& operator=(const MotionModel&) = default;
/// No move assignment
MotionModel& operator=(MotionModel&&) = default;
/// Destructor
~MotionModel() = default;
// - Methods
/// name of the compoenent to index of the component
bool nameToIndex(const word& name, uint32& idx)const
{
return getModel().impl_nameToIndex(name, idx);
}
/// Component index to motion component name
bool indexToName(uint32 idx, word& name)const
{
return getModel().impl_indexToName(idx, name);
}
/// @brief Return a const reference to the list of compoenent names
/// @return
const wordList& componentNames()const
{
return getModel().impl_componentNames();
}
/// Is the wall assocciated to this motion component moving?
bool isMoving()const
{
return getModel().impl_isMoving();
}
/// Move the component itself
bool move(uint32 iter, real t, real dt)
{
return getModel().impl_move(iter, t, dt);
}
/// Obtain an object to model interface
auto getModelInterface(uint32 iter, real t, real dt)const
{
return getModel().impl_getModelInterface(iter, t, dt);
}
};
} // pFlow
#include "MotionModel.cpp"
#endif //__MotionModel_hpp__

4
src/MotionModel/entities/multiRotatingAxis/multiRotatingAxis.hpp
View File

@ -78,7 +78,7 @@ class multiRotatingAxis
{
protected:
/// This is either host/device pointer to all axes
/// This is device pointer to all axes
multiRotatingAxis* axisList_;
/// Index of parent axis
@ -170,7 +170,7 @@ public:
/// Set the pointer to the list of all axes.
/// This pointer is device pointer
INLINE_FUNCTION_H
void setAxisList(multiRotatingAxis* axisList)
void setAxisListPtr(multiRotatingAxis* axisList)
{
axisList_ = axisList;
}

3
src/MotionModel/entities/rotatingAxis/rotatingAxis.cpp
View File

@ -28,7 +28,6 @@ pFlow::rotatingAxis::rotatingAxis
const dictionary& dict
)
{
if(!read(dict))
{
fatalErrorInFunction<<
@ -73,7 +72,7 @@ bool pFlow::rotatingAxis::read
if(!timeInterval::read(dict))return false;
if(!line::read(dict)) return false;
real omega = dict.getValOrSet("omega", static_cast<real>(0.0));
real omega = dict.getValOrSet("omega", 0.0);
setOmega(omega);
return true;

28
src/MotionModel/entities/rotatingAxis/rotatingAxis.hpp
View File

@ -24,13 +24,12 @@ Licence:
#include "timeInterval.hpp"
#include "line.hpp"
#include "rotatingAxisFwd.hpp"
namespace pFlow
{
class dictionary;
class rotatingAxis;
#include "rotatingAxisFwd.hpp"
/**
* An axis which rotates around itself at specified speed
@ -64,7 +63,7 @@ class rotatingAxis
public timeInterval,
public line
{
protected:
private:
/// rotation speed
real omega_ = 0;
@ -76,13 +75,15 @@ public:
// - Constructor
TypeInfoNV("rotatingAxis");
/// Empty constructor
FUNCTION_HD
rotatingAxis(){}
rotatingAxis()=default;
/// Construct from dictionary
FUNCTION_H
rotatingAxis(const dictionary& dict);
explicit rotatingAxis(const dictionary& dict);
/// Construct from components
FUNCTION_HD
@ -92,9 +93,19 @@ public:
FUNCTION_HD
rotatingAxis(const rotatingAxis&) = default;
FUNCTION_HD
rotatingAxis(rotatingAxis&&) = default;
/// Copy asssignment
rotatingAxis& operator=(const rotatingAxis&) = default;
/// Copy asssignment
rotatingAxis& operator=(rotatingAxis&&) = default;
/// destructor
~rotatingAxis()=default;
/// Set omega
FUNCTION_HD
real setOmega(real omega);
@ -115,7 +126,10 @@ public:
/// Linear tangential velocity at point p
INLINE_FUNCTION_HD
realx3 linTangentialVelocityPoint(const realx3 &p)const;
realx3 linVelocityPoint(const realx3 &p)const;
INLINE_FUNCTION_HD
realx3 transferPoint(const realx3 p, real dt)const;
// - IO operation

11
src/MotionModel/entities/rotatingAxis/rotatingAxisFwd.hpp
View File

@ -18,6 +18,11 @@ Licence:
-----------------------------------------------------------------------------*/
namespace pFlow
{
class rotatingAxis;
INLINE_FUNCTION_HD
realx3 rotate(const realx3 &p, const line& ln, real theta);
@ -25,10 +30,10 @@ INLINE_FUNCTION_HD
realx3 rotate(const realx3& p, const rotatingAxis& ax, real dt);
INLINE_FUNCTION_HD
void rotate(realx3* p, size_t n, const line& ln, real theta);
void rotate(realx3* p, uint32 n, const line& ln, real theta);
INLINE_FUNCTION_HD
void rotate(realx3* p, size_t n, const rotatingAxis& ax, real dt);
void rotate(realx3* p, uint32 n, const rotatingAxis& ax, real dt);
}

17
src/MotionModel/entities/rotatingAxis/rotatingAxisI.hpp
View File

@ -1,3 +1,4 @@
#include "rotatingAxis.hpp"
/*------------------------------- phasicFlow ---------------------------------
O C enter of
O O E ngineering and
@ -19,7 +20,7 @@ Licence:
-----------------------------------------------------------------------------*/
INLINE_FUNCTION_HD
pFlow::realx3 pFlow::rotatingAxis::linTangentialVelocityPoint(const realx3 &p)const
pFlow::realx3 pFlow::rotatingAxis::linVelocityPoint(const realx3 &p)const
{
if(!inTimeRange()) return {0,0,0};
@ -28,6 +29,12 @@ pFlow::realx3 pFlow::rotatingAxis::linTangentialVelocityPoint(const realx3 &p)co
return cross(omega_*unitVector(),L);
}
INLINE_FUNCTION_HD
pFlow::realx3 pFlow::rotatingAxis::transferPoint(const realx3 p, real dt)const
{
return rotate(p, *this, dt);
}
INLINE_FUNCTION_HD
pFlow::realx3 pFlow::rotate(const realx3& p, const rotatingAxis& ax, real dt)
{
@ -97,7 +104,7 @@ pFlow::realx3 pFlow::rotate(const realx3 &p, const line& ln, real theta)
}
INLINE_FUNCTION_HD
void pFlow::rotate(realx3* p, size_t n, const line& ln, real theta)
void pFlow::rotate(realx3* p, uint32 n, const line& ln, real theta)
{
realx3 nv = ln.unitVector();
real cos_tet = cos(theta);
@ -110,7 +117,7 @@ void pFlow::rotate(realx3* p, size_t n, const line& ln, real theta)
// (a(v2+w2) - u( bv + cw - ux - vy - wz)) (1-cos_tet) + x cos_tet + (- cv + bw - wy + vz) sin_tet
realx3 res;
for(label i=0; i<n; i++ )
for(uint32 i=0; i<n; i++ )
{
res.x_ = (lp1.x_*(v2 + w2) - (nv.x_*(lp1.y_*nv.y_ + lp1.z_*nv.z_ - nv.x_*p[i].x_ - nv.y_*p[i].y_ - nv.z_*p[i].z_)))*(1 - cos_tet) +
p[i].x_ * cos_tet +
@ -133,7 +140,7 @@ void pFlow::rotate(realx3* p, size_t n, const line& ln, real theta)
}
INLINE_FUNCTION_HD
void pFlow::rotate(realx3* p, size_t n, const rotatingAxis& ax, real dt)
void pFlow::rotate(realx3* p, uint32 n, const rotatingAxis& ax, real dt)
{
if(!ax.inTimeRange()) return;
@ -148,7 +155,7 @@ void pFlow::rotate(realx3* p, size_t n, const rotatingAxis& ax, real dt)
// (a(v2+w2) - u( bv + cw - ux - vy - wz)) (1-cos_tet) + x cos_tet + (- cv + bw - wy + vz) sin_tet
realx3 res;
for(label i=0; i<n; i++ )
for(uint32 i=0; i<n; i++ )
{
res.x_ = (lp1.x_*(v2 + w2) - (nv.x_*(lp1.y_*nv.y_ + lp1.z_*nv.z_ - nv.x_*p[i].x_ - nv.y_*p[i].y_ - nv.z_*p[i].z_)))*(1 - cos_tet) +
p[i].x_ * cos_tet +

57
src/MotionModel/entities/stationary/stationary.cpp Normal file
View File

@ -0,0 +1,57 @@
/*------------------------------- 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.
-----------------------------------------------------------------------------*/
#include "stationary.hpp"
FUNCTION_H
pFlow::stationary::stationary(const dictionary& dict)
{
}
FUNCTION_H
bool pFlow::stationary::read(const dictionary& dict)
{
return true;
}
FUNCTION_H
bool pFlow::stationary::write(dictionary& dict) const
{
return true;
}
FUNCTION_H
bool pFlow::stationary::read(iIstream& is)
{
notImplementedFunction;
return true;
}
FUNCTION_H
bool pFlow::stationary::write(iOstream& os)const
{
return true;
}

104
src/MotionModel/entities/stationary/stationary.hpp Normal file
View File

@ -0,0 +1,104 @@
/*------------------------------- 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 __stationary_hpp__
#define __stationary_hpp__
#include "types.hpp"
#include "typeInfo.hpp"
#include "streams.hpp"
namespace pFlow
{
// forward
class dictionary;
/**
* stationary model for a wall
*
*/
class stationary
{
public:
TypeInfoNV("stationary");
FUNCTION_HD
stationary()=default;
FUNCTION_H
explicit stationary(const dictionary& dict);
FUNCTION_HD
stationary(const stationary&) = default;
stationary& operator=(const stationary&) = default;
INLINE_FUNCTION_HD
void setTime(real t)
{}
INLINE_FUNCTION_HD
realx3 linVelocityPoint(const realx3 &)const
{
return zero3;
}
INLINE_FUNCTION_HD
realx3 transferPoint(const realx3& p, real)const
{
return p;
}
// - IO operation
FUNCTION_H
bool read(const dictionary& dict);
FUNCTION_H
bool write(dictionary& dict) const;
FUNCTION_H
bool read(iIstream& is);
FUNCTION_H
bool write(iOstream& os)const;
};
inline iOstream& operator <<(iOstream& os, const stationary& obj)
{
return os;
}
inline iIstream& operator >>(iIstream& is, stationary& obj)
{
return is;
}
}
#endif

14
src/MotionModel/entities/vibrating/vibrating.hpp
View File

@ -67,7 +67,7 @@ class vibrating
public timeInterval
{
protected:
private:
// rotation speed
realx3 angularFreq_{0,0,0};
@ -94,11 +94,13 @@ protected:
public:
TypeInfoNV("vibrating");
FUNCTION_HD
vibrating(){}
vibrating()=default;
FUNCTION_H
vibrating(const dictionary& dict);
explicit vibrating(const dictionary& dict);
FUNCTION_HD
@ -115,16 +117,16 @@ public:
}
INLINE_FUNCTION_HD
realx3 linTangentialVelocityPoint(const realx3 &p)const
realx3 linVelocityPoint(const realx3 &)const
{
return velocity_;
}
INLINE_FUNCTION_HD
realx3 transferPoint(const realx3& p, real dt)
realx3 transferPoint(const realx3& p, real dt)const
{
if(!inTimeRange()) return p;
return p + static_cast<real>(0.5)*dt*(velocity0_+velocity_);
return p + static_cast<real>(0.5*dt)*(velocity0_+velocity_);
}
// - IO operation

110
src/MotionModel/fixedWall/fixedWall.cpp
View File

@ -1,110 +0,0 @@
/*------------------------------- 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.
-----------------------------------------------------------------------------*/
#include "fixedWall.hpp"
#include "dictionary.hpp"
#include "vocabs.hpp"
bool pFlow::fixedWall::readDictionary
(
const dictionary& dict
)
{
auto motionModel = dict.getVal<word>("motionModel");
if(motionModel != "fixedWall")
{
fatalErrorInFunction<<
" motionModel should be fixedWall, but found " << motionModel <<endl;
return false;
}
return true;
}
bool pFlow::fixedWall::writeDictionary
(
dictionary& dict
)const
{
dict.add("motionModel", "fixedWall");
auto& motionInfo = dict.subDictOrCreate("fixedWallInfo");
return true;
}
pFlow::fixedWall::fixedWall()
{}
pFlow::fixedWall::fixedWall
(
const dictionary& dict
)
{
if(! readDictionary(dict) )
{
fatalExit;
}
}
bool pFlow::fixedWall::read
(
iIstream& is
)
{
// create an empty file dictionary
dictionary motionInfo(motionModelFile__, true);
// read dictionary from stream
if( !motionInfo.read(is) )
{
ioErrorInFile(is.name(), is.lineNumber()) <<
" error in reading dictionray " << motionModelFile__ <<" from file. \n";
return false;
}
if( !readDictionary(motionInfo) ) return false;
return true;
}
bool pFlow::fixedWall::write
(
iOstream& os
)const
{
// create an empty file dictionary
dictionary motionInfo(motionModelFile__, true);
if( !writeDictionary(motionInfo))
{
return false;
}
if( !motionInfo.write(os) )
{
ioErrorInFile( os.name(), os.lineNumber() )<<
" error in writing dictionray to file. \n";
return false;
}
return true;
}

206
src/MotionModel/fixedWall/fixedWall.hpp
View File

@ -1,206 +0,0 @@
/*------------------------------- 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 __fixedWall_hpp__
#define __fixedWall_hpp__
#include "types.hpp"
#include "typeInfo.hpp"
#include "Vectors.hpp"
#include "uniquePtr.hpp"
namespace pFlow
{
class dictionary;
/**
* Fixed wall motion model for walls
*
* This class is used for geometries that are fixed during simulation.
*
*
\verbatim
// In geometryDict file, this will defines fixed walls during simulation
...
motionModel fixedWall;
...
\endverbatim
*/
class fixedWall
{
public:
/** Motion model class to be passed to computational units/kernels for
* transfing points and returning velocities at various positions
*/
class Model
{
public:
INLINE_FUNCTION_HD
Model(){}
INLINE_FUNCTION_HD
Model(const Model&) = default;
INLINE_FUNCTION_HD
Model& operator=(const Model&) = default;
INLINE_FUNCTION_HD
realx3 pointVelocity(int32 n, const realx3 p)const
{
return 0.0;
}
INLINE_FUNCTION_HD
realx3 operator()(int32 n, const realx3& p)const
{
return 0.0;
}
INLINE_FUNCTION_HD
realx3 transferPoint(int32 n, const realx3 p, real dt)const
{
return p;
}
INLINE_FUNCTION_HD int32 numComponents()const
{
return 0;
}
};
protected:
/// Name
const word name_ = "none";
/// Read from a dictionary
bool readDictionary(const dictionary& dict);
/// write to a dictionary
bool writeDictionary(dictionary& dict)const;
public:
/// Type info
TypeInfoNV("fixedWall");
// - Constructors
/// Empty
fixedWall();
/// Constructor with dictionary
fixedWall(const dictionary& dict);
/// Copy constructor
fixedWall(const fixedWall&) = default;
/// Move constructor
fixedWall(fixedWall&&) = default;
/// Copy assignment
fixedWall& operator=(const fixedWall&) = default;
/// Move assignment
fixedWall& operator=(fixedWall&&) = default;
/// Destructor
~fixedWall() = default;
// - Methods
/// Return motion model
/// t is the current simulation time
Model getModel(real t)const
{
return Model();
}
/// Name of the motion component to index
int32 nameToIndex(const word& name)const
{
return 0;
}
/// Index of motion component to name
word indexToName(label i)const
{
return name_;
}
/// Velocity at point p
INLINE_FUNCTION_HD
realx3 pointVelocity(label n, const realx3& p)const
{
return zero3;
}
/// Transfer point p for dt seconds according to motion component n
INLINE_FUNCTION_HD
realx3 transferPoint(label n, const realx3 p, real dt)const
{
return p;
}
/// Transfer a vector of point pVec for dt seconds according to motion
/// component n
INLINE_FUNCTION_HD
bool transferPoint(label n, realx3* pVec, size_t numP, real dt)
{
return true;
}
/// Are walls moving
INLINE_FUNCTION_HD
bool isMoving()const
{
return false;
}
/// Move points
bool move(real t, real dt)
{
return true;
}
// - IO operations
/// Read from input stream is
FUNCTION_H
bool read(iIstream& is);
/// Write to output stream os
FUNCTION_H
bool write(iOstream& os)const;
};
} // pFlow
#endif //__fixed_hpp__

168
src/MotionModel/rotatingAxisMotion/rotatingAxisMotion.cpp
View File

@ -19,152 +19,70 @@ Licence:
-----------------------------------------------------------------------------*/
#include "rotatingAxisMotion.hpp"
#include "dictionary.hpp"
#include "vocabs.hpp"
bool pFlow::rotatingAxisMotion::readDictionary
void pFlow::rotatingAxisMotion::impl_setTime
(
const dictionary& dict
)
{
auto motionModel = dict.getVal<word>("motionModel");
if(motionModel != "rotatingAxisMotion")
{
fatalErrorInFunction<<
" motionModel should be rotatingAxisMotion, but found "
<< motionModel <<endl;
return false;
}
auto& motionInfo = dict.subDict("rotatingAxisMotionInfo");
auto axisNames = motionInfo.dictionaryKeywords();
axis_.reserve(axisNames.size()+1);
axis_.clear();
axisName_.clear();
for(auto& aName: axisNames)
{
auto& axDict = motionInfo.subDict(aName);
if(auto axPtr = makeUnique<rotatingAxis>(axDict); axPtr)
{
axis_.push_back(axPtr());
axisName_.push_back(aName);
}
else
{
fatalErrorInFunction<<
"could not read rotating axis from "<< axDict.globalName()<<endl;
return false;
}
}
if( !axisName_.search("none") )
{
axis_.push_back
(
rotatingAxis(
realx3(0.0,0.0,0.0),
realx3(1.0,0.0,0.0),
0.0
)
);
axisName_.push_back("none");
}
axis_.syncViews();
numAxis_ = axis_.size();
return true;
}
bool pFlow::rotatingAxisMotion::writeDictionary
(
dictionary& dict
uint32 iter,
real t,
real dt
)const
{
dict.add("motionModel", "rotatingAxisMotion");
auto& motionInfo = dict.subDictOrCreate("rotatingAxisMotionInfo");
ForAll(i, axis_)
{
auto& axDict = motionInfo.subDictOrCreate(axisName_[i]);
if( !axis_.hostVectorAll()[i].write(axDict))
{
fatalErrorInFunction<<
" error in writing axis "<< axisName_[i] << " to dicrionary "
<< motionInfo.globalName()<<endl;
return false;
}
}
return true;
auto motion = motionComponents_.deviceViewAll();
Kokkos::parallel_for(
"rotatingAxisMotion::impl_setTime",
deviceRPolicyStatic(0, numComponents_),
LAMBDA_HD(uint32 i){
motion[i].setTime(t);
});
Kokkos::fence();
}
pFlow::rotatingAxisMotion::rotatingAxisMotion()
{}
pFlow::rotatingAxisMotion::rotatingAxisMotion
(
const dictionary& dict
)
pFlow::rotatingAxisMotion::rotatingAxisMotion(
const objectFile &objf,
repository *owner)
: fileDictionary(objf, owner)
{
if(! readDictionary(dict) )
if(! impl_readDictionary(*this) )
{
fatalErrorInFunction;
fatalExit;
}
}
bool pFlow::rotatingAxisMotion::read
pFlow::rotatingAxisMotion::rotatingAxisMotion
(
iIstream& is
const objectFile &objf,
const dictionary &dict,
repository *owner
)
:
fileDictionary(objf, dict, owner)
{
// create an empty file dictionary
dictionary motionInfo(motionModelFile__, true);
// read dictionary from stream
if( !motionInfo.read(is) )
if(!impl_readDictionary(*this) )
{
ioErrorInFile(is.name(), is.lineNumber()) <<
" error in reading dictionray " << motionModelFile__ <<" from file. \n";
return false;
fatalErrorInFunction;
fatalExit;
}
if( !readDictionary(motionInfo) ) return false;
return true;
}
bool pFlow::rotatingAxisMotion::write
(
iOstream& os
)const
iOstream &os,
const IOPattern &iop
) const
{
// create an empty file dictionary
dictionary motionInfo(motionModelFile__, true);
if( !writeDictionary(motionInfo))
// a global dictionary
dictionary newDict(fileDictionary::dictionary::name(), true);
if( iop.thisProcWriteData() )
{
return false;
if( !this->impl_writeDictionary(newDict) ||
!newDict.write(os))
{
fatalErrorInFunction<<
" error in writing to dictionary "<< newDict.globalName()<<endl;
return false;
}
}
if( !motionInfo.write(os) )
{
ioErrorInFile( os.name(), os.lineNumber() )<<
" error in writing dictionray to file. \n";
return false;
}
return true;
}
return true;
}

213
src/MotionModel/rotatingAxisMotion/rotatingAxisMotion.hpp
View File

@ -22,32 +22,26 @@ Licence:
#define __rotatingAxisMotion_hpp__
#include "types.hpp"
#include "typeInfo.hpp"
#include "VectorDual.hpp"
#include "Vectors.hpp"
#include "List.hpp"
#include "MotionModel.hpp"
#include "rotatingAxis.hpp"
#include "fileDictionary.hpp"
namespace pFlow
{
class dictionary;
/**
* Rotating axis motion model for walls
*
* This class is used for simulaiton that at least one wall components
* is moving according to rotatingAxis motion mode. One or more than one
* motion components can be defined in rotatingAxisMotionInfo dictionary
* motion components can be defined in rotatingAxisInfo dictionary
*
\verbatim
// In geometryDict file, this will defines rotating walls during simulation
...
motionModel rotatingAxisMotion;
motionModel rotatingAxis;
rotatingAxisMotionInfo
rotatingAxisInfo
{
rotationAxis1
{
@ -63,185 +57,50 @@ rotatingAxisMotionInfo
*
*/
class rotatingAxisMotion
:
public fileDictionary,
public MotionModel<rotatingAxisMotion, rotatingAxis>
{
public:
/** Motion model class to be passed to computational units/kernels for
* transfing points and returning velocities at various positions
*/
class Model
{
protected:
deviceViewType1D<rotatingAxis> axis_;
int32 numAxis_=0;
public:
INLINE_FUNCTION_HD
Model(deviceViewType1D<rotatingAxis> axis, int32 numAxis):
axis_(axis),
numAxis_(numAxis)
{}
INLINE_FUNCTION_HD
Model(const Model&) = default;
INLINE_FUNCTION_HD
Model& operator=(const Model&) = default;
INLINE_FUNCTION_HD
realx3 pointVelocity(int32 n, const realx3& p)const
{
return axis_[n].linTangentialVelocityPoint(p);
}
INLINE_FUNCTION_HD
realx3 operator()(int32 n, const realx3& p)const
{
return pointVelocity(n,p);
}
INLINE_FUNCTION_HD
realx3 transferPoint(int32 n, const realx3 p, real dt)const
{
return rotate(p, axis_[n], dt);
}
INLINE_FUNCTION_HD int32 numComponents()const
{
return numAxis_;
}
};
protected:
using axisVector_HD = VectorDual<rotatingAxis>;
friend MotionModel<rotatingAxisMotion, rotatingAxis>;
/// Vector to store axes
axisVector_HD axis_;
/// is the geometry attached to this component moving
bool impl_isMoving()const
{
return true;
}
/// move the component itself
bool impl_move(uint32, real, real)const
{
return true;
}
void impl_setTime(uint32 iter, real t, real dt)const;
/// Names of axes
wordList axisName_;
/// Number of axes components
label numAxis_= 0;
/// Read from dictionary
bool readDictionary(const dictionary& dict);
/// Write to dictionary
bool writeDictionary(dictionary& dict)const;
public:
/// Type info
TypeInfoNV("rotatingAxisMotion");
// - Constructors
/// Empty
FUNCTION_H
rotatingAxisMotion();
TypeInfo("rotatingAxisMotion");
/// Construct with dictionary
FUNCTION_H
rotatingAxisMotion(const dictionary& dict);
rotatingAxisMotion(const objectFile& objf, repository* owner);
/// Copy constructor
FUNCTION_H
rotatingAxisMotion(const rotatingAxisMotion&) = default;
rotatingAxisMotion(
const objectFile& objf,
const dictionary& dict,
repository* owner);
/// No move constructor
rotatingAxisMotion(rotatingAxisMotion&&) = delete;
/// Copy assignment
FUNCTION_H
rotatingAxisMotion& operator=(const rotatingAxisMotion&) = default;
/// No move assignment
rotatingAxisMotion& operator=(rotatingAxisMotion&&) = delete;
/// Destructor
FUNCTION_H
~rotatingAxisMotion() = default;
// - Methods
/// Return the motion model at time t
Model getModel(real t)
{
for(int32 i= 0; i<numAxis_; i++ )
{
axis_[i].setTime(t);
}
axis_.modifyOnHost();
axis_.syncViews();
return Model(axis_.deviceVector(), numAxis_);
}
/// Motion component name to index
INLINE_FUNCTION_H
int32 nameToIndex(const word& name)const
{
if( auto i = axisName_.findi(name); i == -1)
{
fatalErrorInFunction<<
"axis name " << name << " does not exist. \n";
fatalExit;
return i;
}
else
{
return i;
}
}
/// Motion index to motion component name
INLINE_FUNCTION_H
word indexToName(label i)const
{
if(i < numAxis_ )
return axisName_[i];
else
{
fatalErrorInFunction<<
"out of range access to the list of axes " << i <<endl<<
" size of axes_ is "<<numAxis_<<endl;
fatalExit;
return "";
}
}
/// Are walls moving
INLINE_FUNCTION_HD
bool isMoving()const
{
return true;
}
/// Move
INLINE_FUNCTION_H
bool move(real t, real dt)
{
return true;
}
// - IO operation
/// Read from input stream is
FUNCTION_H
bool read(iIstream& is);
/// Write to output stream os
FUNCTION_H
bool write(iOstream& os)const;
bool write(iOstream& os, const IOPattern& iop)const override;
static
auto noneComponent()
{
return rotatingAxis({0,0,0}, {1,0,0}, 0.0);
}
};
} // pFlow
#endif //__rotatingAxisMotion_hpp__
#endif // __rotatingAxisMotion_hpp__

Some files were not shown because too many files have changed in this diff Show More