diff --git a/source/Applications/Thermoplast/config.txt b/source/Applications/Thermoplast/config.txt
new file mode 100644
index 0000000000000000000000000000000000000000..392cff19fd6b2be6a89057b918a2ab7dfa0ad51b
--- /dev/null
+++ b/source/Applications/Thermoplast/config.txt
@@ -0,0 +1,21 @@
+#simulation parameters
+#boundingBox = 0 0 0 300 1520 2320
+boundingBox = 60 1370 10 190 1530 320
+blocknx = 10 10 10
+#blocknx = 300 420 320
+endTime = 100000
+outTime = 100
+availMem = 25e9
+uLB = 0.01
+
+#PE parameters
+peMinOffset = 50 2 2
+peMaxOffset = -8 -30 -2
+sphereTime = 1000
+
+#geometry files
+pathGeo = d:/Projects/ThermoPlast/SimGeo
+michel = /michel.stl
+plexiglas = /plexiglas.stl
+
+pathOut = d:/temp/thermoplastCluster
\ No newline at end of file
diff --git a/source/Applications/Thermoplast/thermoplast.cpp b/source/Applications/Thermoplast/thermoplast.cpp
index 5c1bbfeb16b114631576477fbc43e2c91c2956e7..24f26a6184e099791a61ab0b60a8ab89c899d8f8 100644
--- a/source/Applications/Thermoplast/thermoplast.cpp
+++ b/source/Applications/Thermoplast/thermoplast.cpp
@@ -1,426 +1,465 @@
-#include <iostream>
-#include <string>
-
-#include "PointerDefinitions.h"
-
-#include <iostream>
-#include <string>
-#include <memory>
-#include <array>
-
-#include "VirtualFluids.h"
-#include <MuParser/include/muParser.h>
-#include "ForceCalculator.h"
-
-
-#include <MovableObjectInteractor.h>
-#include <DemCalculator.h>
-#include <DemCoProcessor.h>
-#include <PePartitioningGridVisitor.h>
-
-#include <PePhysicsEngineMaterialAdapter.h>
-#include <PePhysicsEngineGeometryAdapter.h>
-#include <PePhysicsEngineSolverAdapter.h>
-
-#include <VelocityBcReconstructor.h>
-#include <EquilibriumReconstructor.h>
-#include <ExtrapolationReconstructor.h>
-
-#include <DummyPhysicsEngineSolverAdapter.h>
-#include <DummyPhysicsEngineMaterialAdapter.h>
-#include <DummyPhysicsEngineGeometryAdapter.h>
-#include <WriteDemObjectsCoProcessor.h>
-
-#include "CreateDemObjectsCoProcessor.h"
-
-using namespace std;
-
-//simulation bounding box
-//double g_minX1 = -81.0;
-//double g_minX2 = -372.0;
-//double g_minX3 = -36.0;
-//
-//double g_maxX1 = 49.0;
-//double g_maxX2 = -318.0;
-//double g_maxX3 = 24;
-
-double g_minX1 = 0;
-double g_minX2 = 0;
-double g_minX3 = 0;
-
-double g_maxX1 = 100;
-double g_maxX2 = 60;
-double g_maxX3 = 60;
-
-string pathOut = "d:/temp/thermoplast3";
-string pathGeo = "d:/Projects/ThermoPlast/Geometrie";
-
-std::shared_ptr<DemCoProcessor> makePeCoProcessor(SPtr<Grid3D> grid, SPtr<Communicator> comm, const SPtr<UbScheduler> peScheduler, const std::shared_ptr<LBMUnitConverter> lbmUnitConverter, int maxpeIterations)
-{
- double peRelaxtion = 0.7;
- //int maxpeIterations = 10000;
- //Beschleunigung g
- double g = 9.81 * lbmUnitConverter->getFactorAccWToLb();
- Vector3D globalLinearAcc(0.0, 0.0, -g);
-
- std::shared_ptr<PePhysicsEngineMaterialAdapter> planeMaterial = std::make_shared<PePhysicsEngineMaterialAdapter>("granular", 1.0, 0, 0.1 / 2, 0.1 / 2, 0.5, 1, 1, 0, 0);
-
- const int gridNx = val<1>(grid->getBlockNX()) * grid->getNX1();
- const int gridNy = val<2>(grid->getBlockNX()) * grid->getNX2();
- const int gridNz = val<3>(grid->getBlockNX()) * grid->getNX3();
-
- //UbTupleInt3 simulationDomain(gridNx, gridNy, gridNz);
- //std::array<double, 6> simulationDomain = {1, 1, 1, 30, 30, 30};
- std::array<double, 6> simulationDomain = {g_minX1, g_minX2, g_minX3, g_maxX1, g_maxX2, g_maxX3};
- SPtr<GbObject3D> boxPE(new GbCuboid3D(simulationDomain[0],simulationDomain[1],simulationDomain[2],simulationDomain[3],simulationDomain[4],simulationDomain[5]));
- GbSystem3D::writeGeoObject(boxPE.get(), pathOut + "/geo/boxPE", WbWriterVtkXmlBinary::getInstance());
- UbTupleInt3 numberOfBlocks(grid->getNX1(), grid->getNX2(), grid->getNX3());
- //UbTupleInt3 numberOfBlocks((simulationDomain[3]-simulationDomain[0])/val<1>(grid->getBlockNX()), (simulationDomain[4]-simulationDomain[1])/val<2>(grid->getBlockNX()), (simulationDomain[5]-simulationDomain[2])/val<3>(grid->getBlockNX()));
- UbTupleBool3 isPeriodic(grid->isPeriodicX1(), grid->isPeriodicX2(), grid->isPeriodicX3());
- Vector3D minOffset(4,2,2);
- Vector3D maxOffset(-2,-2,-2);
- std::shared_ptr<PeParameter> peParamter = std::make_shared<PeParameter>(peRelaxtion, maxpeIterations, globalLinearAcc,
- planeMaterial, simulationDomain, numberOfBlocks, isPeriodic, minOffset, maxOffset);
- std::shared_ptr<PhysicsEngineSolverAdapter> peSolver = std::make_shared<PePhysicsEngineSolverAdapter>(peParamter);
-
- const std::shared_ptr<ForceCalculator> forceCalculator = std::make_shared<ForceCalculator>(comm);
-
- return std::make_shared<DemCoProcessor>(grid, peScheduler, comm, forceCalculator, peSolver);
-}
-
-//pipe flow with forcing
-void pf1()
-{
- SPtr<Communicator> comm = MPICommunicator::getInstance();
- int myid = comm->getProcessID();
-
- //parameters
- //string pathOut = "d:/temp/thermoplast3";
- //string pathGeo = "d:/Projects/ThermoPlast/Geometrie";
- int numOfThreads = 1;
- int blocknx[3] ={ 10,10,10 };
- double endTime = 1000000;
- double outTime = 300;
- double availMem = 8e9;
- double deltax = 1;
- double rhoLB = 0.0;
- double uLB = 0.1;
- double radius = 5;
- double Re = 100;
- double nuLB = (uLB*2.0*radius)/Re;
-
-
-
- //geometry definition
-
- ////simulation bounding box
- //double g_minX1 = -70.0;
- //double g_minX2 = -370.0;
- //double g_minX3 = -25.0;
-
- //double g_maxX1 = 25.0;
- //double g_maxX2 = -330.0;
- //double g_maxX3 = 15;
-
- double blockLength = blocknx[0]*deltax;
-
- SPtr<GbObject3D> gridCube(new GbCuboid3D(g_minX1, g_minX2, g_minX3, g_maxX1, g_maxX2, g_maxX3));
- if (myid == 0) GbSystem3D::writeGeoObject(gridCube.get(), pathOut + "/geo/gridCube", WbWriterVtkXmlBinary::getInstance());
-
- //box
- SPtr<GbObject3D> box(new GbCuboid3D(g_minX1-blockLength, g_minX2, g_minX3, g_maxX1+blockLength, g_maxX2, g_maxX3));
- GbSystem3D::writeGeoObject(box.get(), pathOut + "/geo/box", WbWriterVtkXmlBinary::getInstance());
-
-
-
- //plexiglas
- //SPtr<GbTriFaceMesh3D> plexiglasGeo;
- //if (myid==0) UBLOG(logINFO, "Read plexiglasGeo:start");
- //plexiglasGeo = SPtr<GbTriFaceMesh3D>(GbTriFaceMesh3DCreator::getInstance()->readMeshFromSTLFile2(pathGeo+"/"+"Plexiglas_perforiert1.stl", "plexiglasGeo", GbTriFaceMesh3D::KDTREE_SAHPLIT, false));
- //if (myid==0) UBLOG(logINFO, "Read plexiglasGeo:end");
- //if (myid==0) GbSystem3D::writeGeoObject(plexiglasGeo.get(), pathOut+"/geo/plexiglasGeo", WbWriterVtkXmlBinary::getInstance());
-
- if (myid == 0)
- {
- UBLOG(logINFO, "Parameters:");
- UBLOG(logINFO, "* uLB = " << uLB);
- UBLOG(logINFO, "* rhoLB = " << rhoLB);
- UBLOG(logINFO, "* nuLB = " << nuLB);
- UBLOG(logINFO, "* deltaX = " << deltax);
- UBLOG(logINFO, "* radius = " << radius);
- UBLOG(logINFO, "* Re = " << Re);
- UBLOG(logINFO, "* number of threads = "<<numOfThreads);
- UBLOG(logINFO, "* number of processes = "<<comm->getNumberOfProcesses());
- UBLOG(logINFO, "* path = "<<pathOut);
- UBLOG(logINFO, "Preprocess - start");
- }
-
- //Grid definition
- SPtr<Grid3D> grid(new Grid3D(comm));
- grid->setDeltaX(deltax);
- grid->setBlockNX(blocknx[0], blocknx[1], blocknx[2]);
- grid->setPeriodicX1(false);
- grid->setPeriodicX2(false);
- grid->setPeriodicX3(false);
-
- //blocks generating
- GenBlocksGridVisitor genBlocks(gridCube);
- grid->accept(genBlocks);
-
- //boundary conditions definition
- //boundary conditions adapters
- //////////////////////////////////////////////////////////////////////////////
- SPtr<BCAdapter> noSlipBCAdapter(new NoSlipBCAdapter());
- noSlipBCAdapter->setBcAlgorithm(SPtr<BCAlgorithm>(new NoSlipBCAlgorithm()));
-
- mu::Parser fct;
- fct.SetExpr("U");
- fct.DefineConst("U", uLB);
- SPtr<BCAdapter> inflowAdapter(new VelocityBCAdapter(true, false, false, fct, 0, BCFunction::INFCONST));
- inflowAdapter->setBcAlgorithm(SPtr<BCAlgorithm>(new VelocityBCAlgorithm()));
-
- SPtr<BCAdapter> outflowAdapter(new DensityBCAdapter(rhoLB));
- outflowAdapter->setBcAlgorithm(SPtr<BCAlgorithm>(new NonEqDensityBCAlgorithm()));
-
- //sphere
- mu::Parser fct2;
- fct2.SetExpr("U");
- fct2.DefineConst("U", 0.0);
- SPtr<BCAdapter> velocityBcParticleAdapter(new VelocityBCAdapter(true, false, false, fct2, 0, BCFunction::INFCONST));
- velocityBcParticleAdapter->setBcAlgorithm(SPtr<BCAlgorithm>(new VelocityWithDensityBCAlgorithm()));
-
- //inflow
- GbCuboid3DPtr geoInflow1(new GbCuboid3D(g_minX1-blockLength, g_minX2-radius, g_maxX3-8.0*radius, g_minX1+1, g_minX2+8.0*radius, g_maxX3+radius));
- if (myid == 0) GbSystem3D::writeGeoObject(geoInflow1.get(), pathOut + "/geo/geoInflow1", WbWriterVtkXmlASCII::getInstance());
-
- GbCuboid3DPtr geoInflow2(new GbCuboid3D(g_minX1-blockLength, g_minX2 - blockLength, g_minX3 - blockLength, g_minX1, g_maxX2 + blockLength, g_maxX3 + blockLength));
- if (myid == 0) GbSystem3D::writeGeoObject(geoInflow2.get(), pathOut + "/geo/geoInflow2", WbWriterVtkXmlASCII::getInstance());
-
- GbCuboid3DPtr geoInflow3(new GbCuboid3D(g_minX1-blockLength, g_minX2-radius, g_maxX3-4.0*radius-1, g_minX1+1, g_maxX2+radius, g_maxX3-4.0*radius));
- if (myid == 0) GbSystem3D::writeGeoObject(geoInflow3.get(), pathOut + "/geo/geoInflow3", WbWriterVtkXmlASCII::getInstance());
-
- GbCuboid3DPtr geoInflow4(new GbCuboid3D(g_minX1-blockLength, g_minX2+4.0*radius, g_maxX3-4.0*radius-1.0, g_minX1+1, g_minX2+4.0*radius+1.0, g_maxX3+radius));
- if (myid == 0) GbSystem3D::writeGeoObject(geoInflow4.get(), pathOut + "/geo/geoInflow4", WbWriterVtkXmlASCII::getInstance());
-
- //outflow
- GbCuboid3DPtr geoOutflow(new GbCuboid3D(g_maxX1, g_minX2 - blockLength, g_minX3 - blockLength, g_maxX1 + blockLength, g_maxX2 + blockLength, g_maxX3 + blockLength));
- if (myid == 0) GbSystem3D::writeGeoObject(geoOutflow.get(), pathOut + "/geo/geoOutflow", WbWriterVtkXmlASCII::getInstance());
-
- //boundary conditions visitor
- SPtr<BoundaryConditionsBlockVisitor> bcVisitor(new BoundaryConditionsBlockVisitor());
- bcVisitor->addBC(noSlipBCAdapter);
- bcVisitor->addBC(inflowAdapter);
- bcVisitor->addBC(outflowAdapter);
- bcVisitor->addBC(velocityBcParticleAdapter);
- //////////////////////////////////////////////////////////////////////////////////
-
- //set boundary conditions for blocks and create process decomposition for MPI
- SPtr<D3Q27Interactor> boxInt(new D3Q27Interactor(box, grid, noSlipBCAdapter, Interactor3D::INVERSESOLID));
-
- //sphere bc object
- //SPtr<MovableObjectInteractor> sphereInt1;
- //const std::shared_ptr<Reconstructor> velocityReconstructor = std::make_shared<VelocityBcReconstructor>();
- //const std::shared_ptr<Reconstructor> extrapolationReconstructor = std::make_shared<ExtrapolationReconstructor>(velocityReconstructor);
- //sphereInt1 = SPtr<MovableObjectInteractor>(new MovableObjectInteractor(sphere1, grid, velocityBcParticleAdapter, Interactor3D::SOLID, extrapolationReconstructor, State::UNPIN));
-
- //inflow
- SPtr<D3Q27Interactor> inflowInt1 = SPtr<D3Q27Interactor>(new D3Q27Interactor(geoInflow1, grid, inflowAdapter, Interactor3D::SOLID));
- SPtr<D3Q27Interactor> inflowInt2 = SPtr<D3Q27Interactor>(new D3Q27Interactor(geoInflow2, grid, outflowAdapter, Interactor3D::SOLID));
-
- SPtr<D3Q27Interactor> inflowInt3 = SPtr<D3Q27Interactor>(new D3Q27Interactor(geoInflow3, grid, noSlipBCAdapter, Interactor3D::SOLID));
- SPtr<D3Q27Interactor> inflowInt4 = SPtr<D3Q27Interactor>(new D3Q27Interactor(geoInflow4, grid, noSlipBCAdapter, Interactor3D::SOLID));
-
- //outflow
- SPtr<D3Q27Interactor> outflowInt = SPtr<D3Q27Interactor>(new D3Q27Interactor(geoOutflow, grid, outflowAdapter, Interactor3D::SOLID));
-
- //plexiglas
- //SPtr<Interactor3D> plexiglasInt = SPtr<D3Q27TriFaceMeshInteractor>(new D3Q27TriFaceMeshInteractor(plexiglasGeo, grid, noSlipBCAdapter, Interactor3D::SOLID));
-
- //PE
- double refLengthLb = radius*2.0;
- double refLengthWorld = refLengthLb * deltax;
- const std::shared_ptr<LBMUnitConverter> lbmUnitConverter = std::make_shared<LBMUnitConverter>(refLengthWorld, LBMUnitConverter::WORLD_MATERIAL::AIR_20C, refLengthLb);
- if (myid == 0) std::cout << lbmUnitConverter->toString() << std::endl;
- double rhoSphere = 915 * lbmUnitConverter->getFactorDensityWToLb(); // kg/m^3
- if (myid == 0) UBLOG(logINFO, "rhoSphere = "<<rhoSphere);
- SPtr<PhysicsEngineMaterialAdapter> sphereMaterial(new PePhysicsEngineMaterialAdapter("Polypropylen", rhoSphere, 0, 0.15, 0.1, 0.45, 0.5, 1, 0, 0));
- const int timestep = 2;
- const SPtr<UbScheduler> peScheduler(new UbScheduler(timestep));
- int maxpeIterations = 10;//endTime/2;
- SPtr<DemCoProcessor> demCoProcessor = makePeCoProcessor(grid, comm, peScheduler, lbmUnitConverter, maxpeIterations);
- //demCoProcessor->addInteractor(sphereInt1, sphereMaterial, Vector3D(0.0, 0.0, 0.0));
- //demCoProcessor->distributeIDs();
- demCoProcessor->setBlockVisitor(bcVisitor);
- //double g = 9.81 * lbmUnitConverter->getFactorAccWToLb();
- //double rhoFluid = lbmUnitConverter->getFactorDensityWToLb() * 830.0; // 1 // kg/m^3
- //////////////////////////////////////////////////////////////////////////
-
- SPtr<Grid3DVisitor> peVisitor(new PePartitioningGridVisitor(comm, demCoProcessor));
- InteractorsHelper intHelper(grid, peVisitor);
- intHelper.addInteractor(boxInt);
- //intHelper.addInteractor(sphereInt1);
- //intHelper.addInteractor(inflowInt2);
- intHelper.addInteractor(inflowInt1);
- intHelper.addInteractor(outflowInt);
- intHelper.addInteractor(inflowInt2);
- //intHelper.addInteractor(inflowInt3);
- //intHelper.addInteractor(inflowInt4);
-
- //intHelper.addInteractor(plexiglasInt);
-
- intHelper.selectBlocks();
-
- //write data for visualization of block grid
- SPtr<CoProcessor> ppblocks(new WriteBlocksCoProcessor(grid, SPtr<UbScheduler>(new UbScheduler(1)), pathOut, WbWriterVtkXmlBinary::getInstance(), comm));
- ppblocks->process(0);
- ppblocks.reset();
-
- unsigned long long numberOfBlocks = (unsigned long long)grid->getNumberOfBlocks();
- int ghostLayer = 3;
- unsigned long long numberOfNodesPerBlock = (unsigned long long)(blocknx[0])* (unsigned long long)(blocknx[1])* (unsigned long long)(blocknx[2]);
- unsigned long long numberOfNodes = numberOfBlocks * numberOfNodesPerBlock;
- unsigned long long numberOfNodesPerBlockWithGhostLayer = numberOfBlocks * (blocknx[0] + ghostLayer) * (blocknx[1] + ghostLayer) * (blocknx[2] + ghostLayer);
- double needMemAll = double(numberOfNodesPerBlockWithGhostLayer*(27 * sizeof(double) + sizeof(int) + sizeof(float) * 4));
- double needMem = needMemAll / double(comm->getNumberOfProcesses());
-
- if (myid == 0)
- {
- UBLOG(logINFO, "Number of blocks = " << numberOfBlocks);
- UBLOG(logINFO, "Number of nodes = " << numberOfNodes);
- int minInitLevel = grid->getCoarsestInitializedLevel();
- int maxInitLevel = grid->getFinestInitializedLevel();
- for (int level = minInitLevel; level <= maxInitLevel; level++)
- {
- int nobl = grid->getNumberOfBlocks(level);
- UBLOG(logINFO, "Number of blocks for level " << level << " = " << nobl);
- UBLOG(logINFO, "Number of nodes for level " << level << " = " << nobl*numberOfNodesPerBlock);
- }
- UBLOG(logINFO, "Necessary memory = " << needMemAll << " bytes");
- UBLOG(logINFO, "Necessary memory per process = " << needMem << " bytes");
- UBLOG(logINFO, "Available memory per process = " << availMem << " bytes");
- }
-
- //LBM kernel definition
- SPtr<LBMKernel> kernel;
- kernel = SPtr<LBMKernel>(new IncompressibleCumulantLBMKernel());
- SPtr<BCProcessor> bcProc(new BCProcessor());
- kernel->setBCProcessor(bcProc);
-
- //set forcing
- //mu::Parser fctForcingX1;
- //fctForcingX1.SetExpr("Fx1");
- //fctForcingX1.DefineConst("Fx1", 9e-7);
- //kernel->setWithForcing(true);
- //kernel->setForcingX1(fctForcingX1);
-
- //create LBM kernel
- SetKernelBlockVisitor kernelVisitor(kernel, nuLB, availMem, needMem);
- grid->accept(kernelVisitor);
-
- //set boundary conditions for nodes
- intHelper.setBC();
- grid->accept(*bcVisitor.get());
-
- //initialization of distributions
- InitDistributionsBlockVisitor initVisitor;
- initVisitor.setVx1(uLB);
- grid->accept(initVisitor);
-
- //set connectors
- InterpolationProcessorPtr iProcessor(new IncompressibleOffsetInterpolationProcessor());
- SetConnectorsBlockVisitor setConnsVisitor(comm, true, D3Q27System::ENDDIR, nuLB, iProcessor);
- grid->accept(setConnsVisitor);
-
- //write data for visualization of boundary conditions
- {
- SPtr<UbScheduler> geoSch(new UbScheduler(1));
- WriteBoundaryConditionsCoProcessor ppgeo(grid, geoSch, pathOut, WbWriterVtkXmlBinary::getInstance(), comm);
- ppgeo.process(0);
-
- WriteMacroscopicQuantitiesCoProcessor ppInit(grid, geoSch, pathOut, WbWriterVtkXmlBinary::getInstance(), SPtr<LBMUnitConverter>(new LBMUnitConverter()), comm);
- ppInit.process(0);
- }
-
- if (myid == 0) UBLOG(logINFO, "Preprocess - end");
-
- //write data for visualization of macroscopic quantities
- SPtr<UbScheduler> visSch(new UbScheduler(100));
- SPtr<WriteMacroscopicQuantitiesCoProcessor> writeMQCoProcessor(new WriteMacroscopicQuantitiesCoProcessor(grid, visSch, pathOut,
- WbWriterVtkXmlBinary::getInstance(), SPtr<LBMUnitConverter>(new LBMUnitConverter()), comm));
-
- SPtr<WriteBoundaryConditionsCoProcessor> writeBCCoProcessor(new WriteBoundaryConditionsCoProcessor(grid, visSch, pathOut,
- WbWriterVtkXmlBinary::getInstance(), comm));
-
- SPtr<WriteDemObjectsCoProcessor> writeDemObjectsCoProcessor(new WriteDemObjectsCoProcessor(grid, visSch, pathOut, SPtr<WbWriter>(WbWriterVtkXmlBinary::getInstance()), demCoProcessor, comm));
- writeDemObjectsCoProcessor->process(0);
-
- //performance control
- SPtr<UbScheduler> nupsSch(new UbScheduler(10, 30, 100));
- SPtr<NUPSCounterCoProcessor> npr(new NUPSCounterCoProcessor(grid, nupsSch, numOfThreads, comm));
-
- //////////////////////////////////////////////////////////////////////////
- //generating spheres
- SPtr<UbScheduler> sphereScheduler(new UbScheduler(100));
- Vector3D origin(g_minX1+4.0+radius, g_minX2+4.0+radius, g_maxX3-4.0-4.0*radius);
- double d = 2.0*radius;
- //std::array<double, 6> AABB={origin[0]-radius-1,origin[1]-radius,origin[2]-radius-1,origin[0]+radius+1, origin[1]+2.0*d+radius+1, origin[2]+2.0*d+radius+1};
- //SPtr<GbObject3D> boxAABB(new GbCuboid3D(AABB[0],AABB[1],AABB[2],AABB[3],AABB[4],AABB[5]));
- //GbSystem3D::writeGeoObject(boxAABB.get(), pathOut + "/geo/boxAABB", WbWriterVtkXmlBinary::getInstance());
- SPtr<CreateDemObjectsCoProcessor> createSphereCoProcessor(new CreateDemObjectsCoProcessor(grid,sphereScheduler,demCoProcessor,sphereMaterial,Vector3D(uLB, 0.0, 0.0)));
- //spheres
- for (int x3 = 0; x3 < 2; x3++)
- for (int x2 = 0; x2 < 2; x2++)
- for (int x1 = 0; x1 < 1; x1++)
- {
- //SPtr<GbObject3D> sphere(new GbSphere3D(origin[0]+x1*d, origin[1]+x2*2.0*d, origin[2]+x3*2.0*d, radius));
- SPtr<GbObject3D> sphere(new GbSphere3D(origin[0]+x1*d, origin[1]+x2*1.5*d, origin[2]+x3*1.5*d, radius));
- if (myid == 0) GbSystem3D::writeGeoObject(sphere.get(), pathOut + "/geo/sphere"+UbSystem::toString(x1)+UbSystem::toString(x2)+UbSystem::toString(x3), WbWriterVtkXmlASCII::getInstance());
- createSphereCoProcessor->addGeoObject(sphere);
- }
-
-
- //start simulation
- omp_set_num_threads(numOfThreads);
- //SPtr<UbScheduler> stepGhostLayer(new UbScheduler(outTime));
- SPtr<Calculator> calculator(new BasicCalculator(grid, peScheduler, endTime));
- //SPtr<Calculator> calculator(new DemCalculator(grid, stepGhostLayer, endTime));
- calculator->addCoProcessor(npr);
-
- calculator->addCoProcessor(createSphereCoProcessor);
- calculator->addCoProcessor(demCoProcessor);
-
- calculator->addCoProcessor(writeBCCoProcessor);
- calculator->addCoProcessor(writeDemObjectsCoProcessor);
- calculator->addCoProcessor(writeMQCoProcessor);
-
- if (myid == 0) UBLOG(logINFO, "Simulation-start");
- calculator->calculate();
- if (myid == 0) UBLOG(logINFO, "Simulation-end");
-}
-
-
-//////////////////////////////////////////////////////////////////////////
-int main(int argc, char* argv[])
-{
- try
- {
- //Sleep(30000);
- walberla::Environment env(argc, argv);
-
- pf1();
- //peFlow(std::string(argv[1]));
- return 0;
- }
- catch (std::exception& e)
- {
- cerr << e.what() << endl << flush;
- }
- catch (std::string& s)
- {
- cerr << s << endl;
- }
- catch (...)
- {
- cerr << "unknown exception" << endl;
- }
-}
+#include <iostream>
+#include <string>
+
+#include "PointerDefinitions.h"
+
+#include <iostream>
+#include <string>
+#include <memory>
+#include <array>
+
+#include "VirtualFluids.h"
+#include <MuParser/include/muParser.h>
+#include "ForceCalculator.h"
+
+
+#include <MovableObjectInteractor.h>
+#include <DemCoProcessor.h>
+#include <PePartitioningGridVisitor.h>
+
+#include <PePhysicsEngineMaterialAdapter.h>
+#include <PePhysicsEngineGeometryAdapter.h>
+#include <PePhysicsEngineSolverAdapter.h>
+
+#include <VelocityBcReconstructor.h>
+#include <EquilibriumReconstructor.h>
+#include <ExtrapolationReconstructor.h>
+
+#include <DummyPhysicsEngineSolverAdapter.h>
+#include <DummyPhysicsEngineMaterialAdapter.h>
+#include <DummyPhysicsEngineGeometryAdapter.h>
+#include <WriteDemObjectsCoProcessor.h>
+
+#include "CreateDemObjectsCoProcessor.h"
+
+using namespace std;
+
+//simulation bounding box
+double g_minX1 = 0;
+double g_minX2 = 0;
+double g_minX3 = 0;
+
+double g_maxX1 = 0;
+double g_maxX2 = 0;
+double g_maxX3 = 0;
+
+vector<double> peMinOffset;
+vector<double> peMaxOffset;
+
+string pathOut = "d:/temp/thermoplastCluster";
+string pathGeo = "d:/Projects/ThermoPlast/Geometrie";
+
+std::shared_ptr<DemCoProcessor> makePeCoProcessor(SPtr<Grid3D> grid, SPtr<Communicator> comm, const SPtr<UbScheduler> peScheduler, const std::shared_ptr<LBMUnitConverter> lbmUnitConverter, int maxpeIterations)
+{
+ double peRelaxtion = 0.7;
+ //int maxpeIterations = 10000;
+ //Beschleunigung g
+ double g = 9.81 * lbmUnitConverter->getFactorAccWToLb();
+ Vector3D globalLinearAcc(0.0, -g, 0.0);
+
+ std::shared_ptr<PePhysicsEngineMaterialAdapter> planeMaterial = std::make_shared<PePhysicsEngineMaterialAdapter>("granular", 1.0, 0, 0.1 / 2, 0.1 / 2, 0.5, 1, 1, 0, 0);
+
+ const int gridNx = val<1>(grid->getBlockNX()) * grid->getNX1();
+ const int gridNy = val<2>(grid->getBlockNX()) * grid->getNX2();
+ const int gridNz = val<3>(grid->getBlockNX()) * grid->getNX3();
+
+ //UbTupleInt3 simulationDomain(gridNx, gridNy, gridNz);
+ //std::array<double, 6> simulationDomain = {1, 1, 1, 30, 30, 30};
+ std::array<double, 6> simulationDomain = {g_minX1, g_minX2, g_minX3, g_maxX1, g_maxX2, g_maxX3};
+ UbTupleInt3 numberOfBlocks(grid->getNX1(), grid->getNX2(), grid->getNX3());
+ //UbTupleInt3 numberOfBlocks((simulationDomain[3]-simulationDomain[0])/val<1>(grid->getBlockNX()), (simulationDomain[4]-simulationDomain[1])/val<2>(grid->getBlockNX()), (simulationDomain[5]-simulationDomain[2])/val<3>(grid->getBlockNX()));
+ UbTupleBool3 isPeriodic(grid->isPeriodicX1(), grid->isPeriodicX2(), grid->isPeriodicX3());
+ Vector3D minOffset(peMinOffset[0], peMinOffset[1], peMinOffset[2]);
+ Vector3D maxOffset(peMaxOffset[0], peMaxOffset[1], peMaxOffset[2]);
+
+ SPtr<GbObject3D> boxPE(new GbCuboid3D(simulationDomain[0]+minOffset[0], simulationDomain[1]+minOffset[1], simulationDomain[2]+minOffset[2], simulationDomain[3]+maxOffset[0], simulationDomain[4]+maxOffset[1], simulationDomain[5]+maxOffset[2]));
+ GbSystem3D::writeGeoObject(boxPE.get(), pathOut + "/geo/boxPE", WbWriterVtkXmlBinary::getInstance());
+
+ std::shared_ptr<PeParameter> peParamter = std::make_shared<PeParameter>(peRelaxtion, maxpeIterations, globalLinearAcc,
+ planeMaterial, simulationDomain, numberOfBlocks, isPeriodic, minOffset, maxOffset);
+ std::shared_ptr<PhysicsEngineSolverAdapter> peSolver = std::make_shared<PePhysicsEngineSolverAdapter>(peParamter);
+
+ const std::shared_ptr<ForceCalculator> forceCalculator = std::make_shared<ForceCalculator>(comm);
+
+ return std::make_shared<DemCoProcessor>(grid, peScheduler, comm, forceCalculator, peSolver);
+}
+
+void createGridObjects()
+{
+
+}
+
+//pipe flow with forcing
+void thermoplast(string configname)
+{
+ try
+ {
+ SPtr<Communicator> comm = MPICommunicator::getInstance();
+ int myid = comm->getProcessID();
+
+ ConfigurationFile config;
+ config.load(configname);
+
+ vector<int> blocknx = config.getVector<int>("blocknx");
+ vector<double> boundingBox = config.getVector<double>("boundingBox");
+ peMinOffset = config.getVector<double>("peMinOffset");
+ peMaxOffset = config.getVector<double>("peMaxOffset");
+ int endTime = config.getValue<int>("endTime");
+ double outTime = config.getValue<double>("outTime");
+ double availMem = config.getValue<double>("availMem");
+ double uLB = config.getValue<double>("uLB");
+ string pathOut = config.getValue<string>("pathOut");
+ string pathGeo = config.getValue<string>("pathGeo");
+ string michel = config.getValue<string>("michel");
+ string plexiglas = config.getValue<string>("plexiglas");
+ double sphereTime = config.getValue<double>("sphereTime");
+
+ //parameters
+ //string pathOut = "d:/temp/thermoplast3";
+ //string pathGeo = "d:/Projects/ThermoPlast/Geometrie";
+ int numOfThreads = 1;
+ //int blocknx[3] ={ 10,10,10 };
+ //double endTime = 1000000;
+ //double outTime = 300;
+ //double availMem = 8e9;
+ double deltax = 1;
+ double rhoLB = 0.0;
+ //double uLB = 0.1;
+ double radius = 5;
+ double Re = 300;
+ double nuLB = (uLB*2.0*radius)/Re;
+
+
+
+ //geometry definition
+
+ //simulation bounding box
+ g_minX1 = boundingBox[0];
+ g_minX2 = boundingBox[1];
+ g_minX3 = boundingBox[2];
+
+ g_maxX1 = boundingBox[3];
+ g_maxX2 = boundingBox[4];
+ g_maxX3 = boundingBox[5];
+
+ double blockLength = blocknx[0]*deltax;
+
+ SPtr<GbObject3D> gridCube(new GbCuboid3D(g_minX1, g_minX2, g_minX3, g_maxX1, g_maxX2, g_maxX3));
+ if (myid == 0) GbSystem3D::writeGeoObject(gridCube.get(), pathOut + "/geo/gridCube", WbWriterVtkXmlBinary::getInstance());
+
+ //box
+ SPtr<GbObject3D> box(new GbCuboid3D(g_minX1-blockLength, g_minX2, g_minX3, g_maxX1+blockLength, g_maxX2, g_maxX3));
+ GbSystem3D::writeGeoObject(box.get(), pathOut + "/geo/box", WbWriterVtkXmlBinary::getInstance());
+
+ //michel
+ if (myid==0) UBLOG(logINFO, "Read michelGeo:start");
+ SPtr<GbTriFaceMesh3D> michelGeo = SPtr<GbTriFaceMesh3D>(GbTriFaceMesh3DCreator::getInstance()->readMeshFromSTLFile2(pathGeo+michel, "michelGeo", GbTriFaceMesh3D::KDTREE_SAHPLIT, false));
+ if (myid==0) UBLOG(logINFO, "Read michelGeo:end");
+ if (myid==0) GbSystem3D::writeGeoObject(michelGeo.get(), pathOut+"/geo/michelGeo", WbWriterVtkXmlBinary::getInstance());
+
+ //plexiglas
+ if (myid==0) UBLOG(logINFO, "Read plexiglasGeo:start");
+ SPtr<GbTriFaceMesh3D> plexiglasGeo = SPtr<GbTriFaceMesh3D>(GbTriFaceMesh3DCreator::getInstance()->readMeshFromSTLFile2(pathGeo+plexiglas, "plexiglasGeo", GbTriFaceMesh3D::KDTREE_SAHPLIT, false));
+ if (myid==0) UBLOG(logINFO, "Read plexiglasGeo:end");
+ if (myid==0) GbSystem3D::writeGeoObject(plexiglasGeo.get(), pathOut+"/geo/plexiglasGeo", WbWriterVtkXmlBinary::getInstance());
+
+ if (myid == 0)
+ {
+ UBLOG(logINFO, "Parameters:");
+ UBLOG(logINFO, "* uLB = " << uLB);
+ UBLOG(logINFO, "* rhoLB = " << rhoLB);
+ UBLOG(logINFO, "* nuLB = " << nuLB);
+ UBLOG(logINFO, "* deltaX = " << deltax);
+ UBLOG(logINFO, "* radius = " << radius);
+ UBLOG(logINFO, "* Re = " << Re);
+ UBLOG(logINFO, "* number of threads = "<<numOfThreads);
+ UBLOG(logINFO, "* number of processes = "<<comm->getNumberOfProcesses());
+ UBLOG(logINFO, "* path = "<<pathOut);
+ UBLOG(logINFO, "Preprocess - start");
+ }
+
+ //Grid definition
+ SPtr<Grid3D> grid(new Grid3D(comm));
+ grid->setDeltaX(deltax);
+ grid->setBlockNX(blocknx[0], blocknx[1], blocknx[2]);
+ grid->setPeriodicX1(false);
+ grid->setPeriodicX2(false);
+ grid->setPeriodicX3(false);
+
+ //blocks generating
+ GenBlocksGridVisitor genBlocks(gridCube);
+ grid->accept(genBlocks);
+
+ //boundary conditions definition
+ //boundary conditions adapters
+ //////////////////////////////////////////////////////////////////////////////
+ SPtr<BCAdapter> noSlipBCAdapter(new NoSlipBCAdapter());
+ noSlipBCAdapter->setBcAlgorithm(SPtr<BCAlgorithm>(new NoSlipBCAlgorithm()));
+
+ mu::Parser fct;
+ fct.SetExpr("U");
+ fct.DefineConst("U", uLB);
+ SPtr<BCAdapter> inflowAdapter(new VelocityBCAdapter(true, false, false, fct, 0, BCFunction::INFCONST));
+ inflowAdapter->setBcAlgorithm(SPtr<BCAlgorithm>(new VelocityBCAlgorithm()));
+
+ SPtr<BCAdapter> outflowAdapter(new DensityBCAdapter(rhoLB));
+ outflowAdapter->setBcAlgorithm(SPtr<BCAlgorithm>(new NonEqDensityBCAlgorithm()));
+ //outflowAdapter->setBcAlgorithm(SPtr<BCAlgorithm>(new NonReflectingOutflowBCAlgorithm()));
+
+ //sphere
+ mu::Parser fct2;
+ fct2.SetExpr("U");
+ fct2.DefineConst("U", 0.0);
+ SPtr<BCAdapter> velocityBcParticleAdapter(new VelocityBCAdapter(true, false, false, fct2, 0, BCFunction::INFCONST));
+ velocityBcParticleAdapter->setBcAlgorithm(SPtr<BCAlgorithm>(new VelocityWithDensityBCAlgorithm()));
+
+ //inflow
+ GbCuboid3DPtr geoInflow1(new GbCuboid3D(g_minX1-blockLength, g_maxX2-120.0, g_minX3+190.0, g_minX1+1, g_maxX2+20.0, g_minX3+320.0));
if (myid == 0) GbSystem3D::writeGeoObject(geoInflow1.get(), pathOut + "/geo/geoInflow1", WbWriterVtkXmlASCII::getInstance());
+
+ GbCuboid3DPtr geoInflow2(new GbCuboid3D(g_minX1-blockLength, g_minX2 - blockLength, g_minX3 - blockLength, g_minX1, g_maxX2 + blockLength, g_maxX3 + blockLength));
+ if (myid == 0) GbSystem3D::writeGeoObject(geoInflow2.get(), pathOut + "/geo/geoInflow2", WbWriterVtkXmlASCII::getInstance());
+
+ GbCuboid3DPtr geoInflow3(new GbCuboid3D(g_minX1-blockLength, g_minX2-radius, g_maxX3-4.0*radius-1, g_minX1+1, g_maxX2+radius, g_maxX3-4.0*radius));
+ if (myid == 0) GbSystem3D::writeGeoObject(geoInflow3.get(), pathOut + "/geo/geoInflow3", WbWriterVtkXmlASCII::getInstance());
+
+ GbCuboid3DPtr geoInflow4(new GbCuboid3D(g_minX1-blockLength, g_minX2+4.0*radius, g_maxX3-4.0*radius-1.0, g_minX1+1, g_minX2+4.0*radius+1.0, g_maxX3+radius));
+ if (myid == 0) GbSystem3D::writeGeoObject(geoInflow4.get(), pathOut + "/geo/geoInflow4", WbWriterVtkXmlASCII::getInstance());
+
+ //outflow
+ GbCuboid3DPtr geoOutflow(new GbCuboid3D(g_maxX1, g_minX2 - blockLength, g_minX3 - blockLength, g_maxX1 + blockLength, g_maxX2 + blockLength, g_maxX3 + blockLength));
+ if (myid == 0) GbSystem3D::writeGeoObject(geoOutflow.get(), pathOut + "/geo/geoOutflow", WbWriterVtkXmlASCII::getInstance());
+
+ //boundary conditions visitor
+ SPtr<BoundaryConditionsBlockVisitor> bcVisitor(new BoundaryConditionsBlockVisitor());
+ bcVisitor->addBC(noSlipBCAdapter);
+ bcVisitor->addBC(inflowAdapter);
+ bcVisitor->addBC(outflowAdapter);
+ bcVisitor->addBC(velocityBcParticleAdapter);
+ //////////////////////////////////////////////////////////////////////////////////
+
+ //set boundary conditions for blocks and create process decomposition for MPI
+ SPtr<D3Q27Interactor> boxInt(new D3Q27Interactor(box, grid, noSlipBCAdapter, Interactor3D::INVERSESOLID));
+
+ //sphere bc object
+ //SPtr<MovableObjectInteractor> sphereInt1;
+ //const std::shared_ptr<Reconstructor> velocityReconstructor = std::make_shared<VelocityBcReconstructor>();
+ //const std::shared_ptr<Reconstructor> extrapolationReconstructor = std::make_shared<ExtrapolationReconstructor>(velocityReconstructor);
+ //sphereInt1 = SPtr<MovableObjectInteractor>(new MovableObjectInteractor(sphere1, grid, velocityBcParticleAdapter, Interactor3D::SOLID, extrapolationReconstructor, State::UNPIN));
+
+ //inflow
+ SPtr<D3Q27Interactor> inflowInt1 = SPtr<D3Q27Interactor>(new D3Q27Interactor(geoInflow1, grid, inflowAdapter, Interactor3D::SOLID));
+ SPtr<D3Q27Interactor> inflowInt2 = SPtr<D3Q27Interactor>(new D3Q27Interactor(geoInflow2, grid, outflowAdapter, Interactor3D::SOLID));
+
+ SPtr<D3Q27Interactor> inflowInt3 = SPtr<D3Q27Interactor>(new D3Q27Interactor(geoInflow3, grid, noSlipBCAdapter, Interactor3D::SOLID));
+ SPtr<D3Q27Interactor> inflowInt4 = SPtr<D3Q27Interactor>(new D3Q27Interactor(geoInflow4, grid, noSlipBCAdapter, Interactor3D::SOLID));
+
+ //outflow
+ SPtr<D3Q27Interactor> outflowInt = SPtr<D3Q27Interactor>(new D3Q27Interactor(geoOutflow, grid, outflowAdapter, Interactor3D::SOLID));
+
+ //michel
+ SPtr<Interactor3D> michelInt = SPtr<D3Q27TriFaceMeshInteractor>(new D3Q27TriFaceMeshInteractor(michelGeo, grid, noSlipBCAdapter, Interactor3D::SOLID));
+
+ //plexiglas
+ SPtr<Interactor3D> plexiglasInt = SPtr<D3Q27TriFaceMeshInteractor>(new D3Q27TriFaceMeshInteractor(plexiglasGeo, grid, noSlipBCAdapter, Interactor3D::SOLID));
+
+ //PE
+ double refLengthLb = radius*2.0;
+ double refLengthWorld = refLengthLb * deltax;
+ const std::shared_ptr<LBMUnitConverter> lbmUnitConverter = std::make_shared<LBMUnitConverter>(refLengthWorld, LBMUnitConverter::WORLD_MATERIAL::AIR_20C, refLengthLb);
+ if (myid == 0) std::cout << lbmUnitConverter->toString() << std::endl;
+ double rhoSphere = 915 * lbmUnitConverter->getFactorDensityWToLb(); // kg/m^3
+ if (myid == 0) UBLOG(logINFO, "rhoSphere = "<<rhoSphere);
+ SPtr<PhysicsEngineMaterialAdapter> sphereMaterial(new PePhysicsEngineMaterialAdapter("Polypropylen", rhoSphere, 0, 0.15, 0.1, 0.45, 0.5, 1, 0, 0));
+ const int timestep = 2;
+ const SPtr<UbScheduler> peScheduler(new UbScheduler(timestep));
+ int maxpeIterations = 10;//endTime/2;
+ SPtr<DemCoProcessor> demCoProcessor = makePeCoProcessor(grid, comm, peScheduler, lbmUnitConverter, maxpeIterations);
+ //demCoProcessor->addInteractor(sphereInt1, sphereMaterial, Vector3D(0.0, 0.0, 0.0));
+ //demCoProcessor->distributeIDs();
+ demCoProcessor->setBlockVisitor(bcVisitor);
+ //double g = 9.81 * lbmUnitConverter->getFactorAccWToLb();
+ //double rhoFluid = lbmUnitConverter->getFactorDensityWToLb() * 830.0; // 1 // kg/m^3
+ //////////////////////////////////////////////////////////////////////////
+ SPtr<Grid3DVisitor> peVisitor(new PePartitioningGridVisitor(comm, demCoProcessor));
+ //grid->accept(peVisitor);
+
+ //SetSolidBlocksBlockVisitor solidBlocksVisitor1(michelInt);
+ //grid->accept(solidBlocksVisitor1);
+ //SetSolidBlocksBlockVisitor solidBlocksVisitor2(plexiglasInt);
+ //grid->accept(solidBlocksVisitor2);
+
+ //SetBcBlocksBlockVisitor bcBlocksVisitor1(michelInt);
+ //grid->accept(bcBlocksVisitor1);
+ //SetBcBlocksBlockVisitor bcBlocksVisitor2(plexiglasInt);
+ //grid->accept(bcBlocksVisitor2);
+
+ InteractorsHelper intHelper(grid, peVisitor);
+ intHelper.addInteractor(boxInt);
+ //intHelper.addInteractor(sphereInt1);
+ //intHelper.addInteractor(inflowInt2);
+ intHelper.addInteractor(michelInt);
+ intHelper.addInteractor(plexiglasInt);
+ intHelper.addInteractor(inflowInt1);
+ intHelper.addInteractor(outflowInt);
+ intHelper.addInteractor(inflowInt2);
+ //intHelper.addInteractor(inflowInt3);
+ //intHelper.addInteractor(inflowInt4);
+
+ //intHelper.addInteractor(plexiglasInt);
+
+ intHelper.selectBlocks();
+
+ //write data for visualization of block grid
+ SPtr<CoProcessor> ppblocks(new WriteBlocksCoProcessor(grid, SPtr<UbScheduler>(new UbScheduler(1)), pathOut, WbWriterVtkXmlBinary::getInstance(), comm));
+ ppblocks->process(0);
+ ppblocks.reset();
+
+ unsigned long long numberOfBlocks = (unsigned long long)grid->getNumberOfBlocks();
+ int ghostLayer = 3;
+ unsigned long long numberOfNodesPerBlock = (unsigned long long)(blocknx[0])* (unsigned long long)(blocknx[1])* (unsigned long long)(blocknx[2]);
+ unsigned long long numberOfNodes = numberOfBlocks * numberOfNodesPerBlock;
+ unsigned long long numberOfNodesPerBlockWithGhostLayer = numberOfBlocks * (blocknx[0] + ghostLayer) * (blocknx[1] + ghostLayer) * (blocknx[2] + ghostLayer);
+ double needMemAll = double(numberOfNodesPerBlockWithGhostLayer*(27 * sizeof(double) + sizeof(int) + sizeof(float) * 4));
+ double needMem = needMemAll / double(comm->getNumberOfProcesses());
+
+ if (myid == 0)
+ {
+ UBLOG(logINFO, "Number of blocks = " << numberOfBlocks);
+ UBLOG(logINFO, "Number of nodes = " << numberOfNodes);
+ int minInitLevel = grid->getCoarsestInitializedLevel();
+ int maxInitLevel = grid->getFinestInitializedLevel();
+ for (int level = minInitLevel; level <= maxInitLevel; level++)
+ {
+ int nobl = grid->getNumberOfBlocks(level);
+ UBLOG(logINFO, "Number of blocks for level " << level << " = " << nobl);
+ UBLOG(logINFO, "Number of nodes for level " << level << " = " << nobl*numberOfNodesPerBlock);
+ }
+ UBLOG(logINFO, "Necessary memory = " << needMemAll << " bytes");
+ UBLOG(logINFO, "Necessary memory per process = " << needMem << " bytes");
+ UBLOG(logINFO, "Available memory per process = " << availMem << " bytes");
+ }
+
+ //LBM kernel definition
+ SPtr<LBMKernel> kernel;
+ kernel = SPtr<LBMKernel>(new IncompressibleCumulantLBMKernel());
+ SPtr<BCProcessor> bcProc(new BCProcessor());
+ kernel->setBCProcessor(bcProc);
+
+ //set forcing
+ //mu::Parser fctForcingX1;
+ //fctForcingX1.SetExpr("Fx1");
+ //fctForcingX1.DefineConst("Fx1", 9e-7);
+ //kernel->setWithForcing(true);
+ //kernel->setForcingX1(fctForcingX1);
+
+ //create LBM kernel
+ SetKernelBlockVisitor kernelVisitor(kernel, nuLB, availMem, needMem);
+ grid->accept(kernelVisitor);
+
+ //set boundary conditions for nodes
+ //michelInt->initInteractor();
+ //plexiglasInt->initInteractor();
+
+ intHelper.setBC();
+ grid->accept(*bcVisitor.get());
+
+ //initialization of distributions
+ InitDistributionsBlockVisitor initVisitor;
+ initVisitor.setVx1(uLB);
+ grid->accept(initVisitor);
+
+ //set connectors
+ InterpolationProcessorPtr iProcessor(new IncompressibleOffsetInterpolationProcessor());
+ SetConnectorsBlockVisitor setConnsVisitor(comm, true, D3Q27System::ENDDIR, nuLB, iProcessor);
+ grid->accept(setConnsVisitor);
+
+ //write data for visualization of boundary conditions
+ {
+ SPtr<UbScheduler> geoSch(new UbScheduler(1));
+ WriteBoundaryConditionsCoProcessor ppgeo(grid, geoSch, pathOut, WbWriterVtkXmlBinary::getInstance(), comm);
+ ppgeo.process(0);
+
+ WriteMacroscopicQuantitiesCoProcessor ppInit(grid, geoSch, pathOut, WbWriterVtkXmlBinary::getInstance(), SPtr<LBMUnitConverter>(new LBMUnitConverter()), comm);
+ ppInit.process(0);
+ }
+
+ if (myid == 0) UBLOG(logINFO, "Preprocess - end");
+
+ //write data for visualization of macroscopic quantities
+ SPtr<UbScheduler> visSch(new UbScheduler(outTime));
+ SPtr<WriteMacroscopicQuantitiesCoProcessor> writeMQCoProcessor(new WriteMacroscopicQuantitiesCoProcessor(grid, visSch, pathOut,
+ WbWriterVtkXmlBinary::getInstance(), SPtr<LBMUnitConverter>(new LBMUnitConverter()), comm));
+
+ SPtr<WriteBoundaryConditionsCoProcessor> writeBCCoProcessor(new WriteBoundaryConditionsCoProcessor(grid, visSch, pathOut,
+ WbWriterVtkXmlBinary::getInstance(), comm));
+
+ SPtr<WriteDemObjectsCoProcessor> writeDemObjectsCoProcessor(new WriteDemObjectsCoProcessor(grid, visSch, pathOut, SPtr<WbWriter>(WbWriterVtkXmlBinary::getInstance()), demCoProcessor, comm));
+ writeDemObjectsCoProcessor->process(0);
+
+ //performance control
+ SPtr<UbScheduler> nupsSch(new UbScheduler(10, 30, 100));
+ SPtr<NUPSCounterCoProcessor> npr(new NUPSCounterCoProcessor(grid, nupsSch, numOfThreads, comm));
+
+ //////////////////////////////////////////////////////////////////////////
+ //generating spheres
+ double d = 2.0*radius;
SPtr<UbScheduler> sphereScheduler(new UbScheduler(sphereTime));
Vector3D origin(g_minX1+60.0+4.0+radius, geoInflow1->getX2Minimum()+4.0*d, geoInflow1->getX3Minimum()+2.0*d);
+ //std::array<double, 6> AABB={origin[0]-radius-1,origin[1]-radius,origin[2]-radius-1,origin[0]+radius+1, origin[1]+2.0*d+radius+1, origin[2]+2.0*d+radius+1};
+ //SPtr<GbObject3D> boxAABB(new GbCuboid3D(AABB[0],AABB[1],AABB[2],AABB[3],AABB[4],AABB[5]));
+ //GbSystem3D::writeGeoObject(boxAABB.get(), pathOut + "/geo/boxAABB", WbWriterVtkXmlBinary::getInstance());
+ SPtr<CreateDemObjectsCoProcessor> createSphereCoProcessor(new CreateDemObjectsCoProcessor(grid,sphereScheduler,demCoProcessor,sphereMaterial,Vector3D(uLB, 0.0, 0.0)));
+ //spheres
+ for (int x3 = 0; x3 < 6; x3++)
for (int x2 = 0; x2 < 4; x2++)
for (int x1 = 0; x1 < 1; x1++)
+ {
+ //SPtr<GbObject3D> sphere(new GbSphere3D(origin[0]+x1*d, origin[1]+x2*2.0*d, origin[2]+x3*2.0*d, radius));
+ SPtr<GbObject3D> sphere(new GbSphere3D(origin[0]+x1*d, origin[1]+x2*1.5*d, origin[2]+x3*1.5*d, radius));
+ if (myid == 0) GbSystem3D::writeGeoObject(sphere.get(), pathOut + "/geo/sphere"+UbSystem::toString(x1)+UbSystem::toString(x2)+UbSystem::toString(x3), WbWriterVtkXmlASCII::getInstance());
+ createSphereCoProcessor->addGeoObject(sphere);
+ }
+
+ createSphereCoProcessor->process(0);
+
+ //start simulation
+ omp_set_num_threads(numOfThreads);
+ SPtr<UbScheduler> stepGhostLayer(visSch);
+ SPtr<Calculator> calculator(new BasicCalculator(grid, peScheduler, endTime));
+ calculator->addCoProcessor(npr);
+
+ calculator->addCoProcessor(createSphereCoProcessor);
+ calculator->addCoProcessor(demCoProcessor);
+
+ calculator->addCoProcessor(writeBCCoProcessor);
+ calculator->addCoProcessor(writeDemObjectsCoProcessor);
+ calculator->addCoProcessor(writeMQCoProcessor);
+
+ if (myid == 0) UBLOG(logINFO, "Simulation-start");
+ calculator->calculate();
+ if (myid == 0) UBLOG(logINFO, "Simulation-end");
+
+ }
+ catch (std::exception& e)
{
UBLOG(logERROR, e.what());
}
catch (std::string& s)
{
UBLOG(logERROR, s);
}
catch (...)
{
UBLOG(logERROR, "unknown exception");
}
+}
+
+
+//////////////////////////////////////////////////////////////////////////
+int main(int argc, char* argv[])
+{
+ //try
+ //{
+ //Sleep(30000);
+ walberla::Environment env(argc, argv);
+
+ if (argv!=NULL)
+ {
+ if (argv[1]!=NULL)
+ {
+ thermoplast(string(argv[1]));
+ }
+ else
+ {
+ cout<<"Configuration file must be set!: "<<argv[0]<<" <config file>"<<endl<<std::flush;
+ }
+ }
+ return 0;
+ //}
+ //catch (std::exception& e)
//{
// UBLOG(logERROR, e.what());
//}
//catch (std::string& s)
//{
// UBLOG(logERROR, s);
//}
//catch (...)
//{
// UBLOG(logERROR, "unknown exception");
//}
+}
diff --git a/source/Applications/VirtualFluids.h b/source/Applications/VirtualFluids.h
index a6e7649987158160b625f69bf36757088255d713..95a3197eef0fffa99cd7d033317bb7d2bd3c80f9 100644
--- a/source/Applications/VirtualFluids.h
+++ b/source/Applications/VirtualFluids.h
@@ -272,7 +272,8 @@
#include <Visitors/SetKernelBlockVisitor.h>
#include <Visitors/SetForcingBlockVisitor.h>
#include <Visitors/SetSpongeLayerBlockVisitor.h>
-#include <Visitors/SetSolidOrBoundaryBlockVisitor.h>
+#include <Visitors/SetSolidBlocksBlockVisitor.h>
+#include <Visitors/SetBcBlocksBlockVisitor.h>
#include <Visitors/RenumberBlockVisitor.h>
#include <Visitors/ConnectorBlockVisitor.h>
#include <Visitors/ViscosityBlockVisitor.h>
diff --git a/source/CMake/cmake_config_files/PHOENIX.config.cmake b/source/CMake/cmake_config_files/PHOENIX.config.cmake
index b2f6bc0bb137cd0f0ca4754fa29614b1384d4117..ce204f104eb27aba236a2c9b5f1d563cb9978d95 100644
--- a/source/CMake/cmake_config_files/PHOENIX.config.cmake
+++ b/source/CMake/cmake_config_files/PHOENIX.config.cmake
@@ -28,6 +28,4 @@ IF(${USE_DEM_COUPLING})
SET(CORE_DEBUG_LIBRARY ${PE_BINARY_DIR}/src/core/libcore.a)
SET(CORE_RELEASE_LIBRARY ${PE_BINARY_DIR}/src/core/libcore.a)
- SET(CAB_ADDITIONAL_LINK_LIBRARIES ${CAB_ADDITIONAL_LINK_LIBRARIES} "stdc++fs")
-
ENDIF()
\ No newline at end of file
diff --git a/source/DemCoupling/CreateDemObjectsCoProcessor.cpp b/source/DemCoupling/CreateDemObjectsCoProcessor.cpp
index aa76ffc100439c78efcb45a620842aa8b9a2806f..28311218eeaa85cb8440d008c8e88ec2adf3bbb8 100644
--- a/source/DemCoupling/CreateDemObjectsCoProcessor.cpp
+++ b/source/DemCoupling/CreateDemObjectsCoProcessor.cpp
@@ -14,7 +14,7 @@
#include "PePhysicsEngineMaterialAdapter.h"
#include "muParser.h"
#include "PhysicsEngineMaterialAdapter.h"
-#include "SetSolidOrBoundaryBlockVisitor.h"
+#include "SetBcBlocksBlockVisitor.h"
#include "Grid3D.h"
CreateDemObjectsCoProcessor::CreateDemObjectsCoProcessor(SPtr<Grid3D> grid, SPtr<UbScheduler> s, SPtr<DemCoProcessor> demCoProcessor, SPtr<PhysicsEngineMaterialAdapter> demObjectMaterial, Vector3D initalVelocity) :
@@ -52,9 +52,7 @@ void CreateDemObjectsCoProcessor::process(double step)
}
SPtr<GbObject3D> geoObject((GbObject3D*)(proto->clone()));
SPtr<MovableObjectInteractor> geoObjectInt = SPtr<MovableObjectInteractor>(new MovableObjectInteractor(geoObject, grid, velocityBcParticleAdapter, Interactor3D::SOLID, extrapolationReconstructor, State::UNPIN));
- //SetSolidOrBoundaryBlockVisitor setSolidVisitor(geoObjectInt, SetSolidOrBoundaryBlockVisitor::SOLID);
- //grid->accept(setSolidVisitor);
- SetSolidOrBoundaryBlockVisitor setBcVisitor(geoObjectInt, SetSolidOrBoundaryBlockVisitor::BC);
+ SetBcBlocksBlockVisitor setBcVisitor(geoObjectInt);
grid->accept(setBcVisitor);
geoObjectInt->initInteractor();
demCoProcessor->addInteractor(geoObjectInt, demObjectMaterial, initalVelocity);
diff --git a/source/DemCoupling/DemCalculator.cpp b/source/DemCoupling/DemCalculator.cpp
deleted file mode 100644
index 5aa4690413921ef326f49eb1af8fce44fc54237d..0000000000000000000000000000000000000000
--- a/source/DemCoupling/DemCalculator.cpp
+++ /dev/null
@@ -1,134 +0,0 @@
-#include "DemCalculator.h"
-
-
-#include "Grid3D.h"
-//#include "Synchronizer.h"
-#include "UbLogger.h"
-#include "Communicator.h"
-#include "TimeAveragedValuesCoProcessor.h"
-#include "UbScheduler.h"
-
-
-DemCalculator::DemCalculator(SPtr<Grid3D> grid, SPtr<UbScheduler> additionalGhostLayerUpdateScheduler, int numberOfTimeSteps) :
- BasicCalculator(grid, additionalGhostLayerUpdateScheduler, numberOfTimeSteps)
-{
-
-}
-
-void DemCalculator::calculate()
-{
- UBLOG(logDEBUG1, "OMPCalculator::calculate() - started");
- int calcStep = 0;
- try
- {
- int minInitLevel = minLevel;
- int maxInitLevel = maxLevel - minLevel;
- int straightStartLevel = minInitLevel;
- int internalIterations = 1 << (maxInitLevel - minInitLevel);
- int forwardStartLevel;
- int threshold;
-
-#ifdef TIMING
- UbTimer timer;
- double time[6];
-#endif
-
- for (calcStep = startTimeStep; calcStep <= numberOfTimeSteps + 1; calcStep++)
- {
- coProcess((double)(calcStep - 1));
-
- exchangeBlockData(straightStartLevel, maxInitLevel);
-
-
- //////////////////////////////////////////////////////////////////////////
-#ifdef TIMING
- UBLOG(logINFO, "calcStep = " << calcStep);
-#endif
- //////////////////////////////////////////////////////////////////////////
-
- for (int staggeredStep = 1; staggeredStep <= internalIterations; staggeredStep++)
- {
- forwardStartLevel = straightStartLevel;
- if (staggeredStep == internalIterations) straightStartLevel = minInitLevel;
- else
- {
- for (straightStartLevel = maxInitLevel, threshold = 1;
- (staggeredStep&threshold) != threshold; straightStartLevel--, threshold <<= 1);
- }
- //////////////////////////////////////////////////////////////////////////
-#ifdef TIMING
- timer.resetAndStart();
-#endif
- //////////////////////////////////////////////////////////////////////////
- applyPreCollisionBC(straightStartLevel, maxInitLevel);
-
- //do collision for all blocks
- calculateBlocks(straightStartLevel, maxInitLevel, calcStep);
- //////////////////////////////////////////////////////////////////////////
-#ifdef TIMING
- time[0] = timer.stop();
- UBLOG(logINFO, "calculateBlocks time = " << time[0]);
-#endif
- //////////////////////////////////////////////////////////////////////////
- //////////////////////////////////////////////////////////////////////////
- //exchange data between blocks
- exchangeBlockData(straightStartLevel, maxInitLevel);
- //////////////////////////////////////////////////////////////////////////
-#ifdef TIMING
- time[1] = timer.stop();
- UBLOG(logINFO, "exchangeBlockData time = " << time[1]);
-#endif
- //////////////////////////////////////////////////////////////////////////
- applyPostCollisionBC(straightStartLevel, maxInitLevel);
- //////////////////////////////////////////////////////////////////////////
-#ifdef TIMING
- time[2] = timer.stop();
- UBLOG(logINFO, "applyBCs time = " << time[2]);
-#endif
- //////////////////////////////////////////////////////////////////////////
- //swap distributions in kernel
- swapDistributions(straightStartLevel, maxInitLevel);
- //////////////////////////////////////////////////////////////////////////
-#ifdef TIMING
- time[3] = timer.stop();
- UBLOG(logINFO, "swapDistributions time = " << time[3]);
-#endif
- //////////////////////////////////////////////////////////////////////////
- if (refinement)
- {
- if (straightStartLevel < maxInitLevel)
- exchangeBlockData(straightStartLevel, maxInitLevel);
- //////////////////////////////////////////////////////////////////////////
-#ifdef TIMING
- time[4] = timer.stop();
- UBLOG(logINFO, "refinement exchangeBlockData time = " << time[4]);
-#endif
- //////////////////////////////////////////////////////////////////////////
- //now ghost nodes have actual values
- //interpolation of interface nodes between grid levels
- interpolation(straightStartLevel, maxInitLevel);
- //////////////////////////////////////////////////////////////////////////
-#ifdef TIMING
- time[5] = timer.stop();
- UBLOG(logINFO, "refinement interpolation time = " << time[5]);
-#endif
- //////////////////////////////////////////////////////////////////////////
- }
- }
- //exchange data between blocks for visualization
- if ((int)additionalGhostLayerUpdateScheduler->getNextDueTime() == calcStep)
- {
- exchangeBlockData(straightStartLevel, maxInitLevel);
- }
- //now ghost nodes have actual values
- }
- UBLOG(logDEBUG1, "OMPCalculator::calculate() - stoped");
- }
- catch (std::exception& e)
- {
- UBLOG(logERROR, e.what());
- UBLOG(logERROR, " step = " << calcStep);
- //throw;
- exit(EXIT_FAILURE);
- }
-}
diff --git a/source/DemCoupling/DemCalculator.h b/source/DemCoupling/DemCalculator.h
deleted file mode 100644
index 10ce6c841e286d4afa4d5ca0a40f6d5eaa5ac75b..0000000000000000000000000000000000000000
--- a/source/DemCoupling/DemCalculator.h
+++ /dev/null
@@ -1,27 +0,0 @@
-/*
-* Author: S. Peters
-* mail: peters@irmb.tu-bs.de
-*/
-#ifndef DEM_CALCULATOR_H
-#define DEM_CALCULATOR_H
-
-#include <memory>
-
-#include "BasicCalculator.h"
-
-class Grid3D;
-class Synchronizer;
-class DemCoProcessor;
-class CalculationManager;
-
-class DemCalculator : public BasicCalculator
-{
-public:
- DemCalculator(SPtr<Grid3D> grid, SPtr<UbScheduler> additionalGhostLayerUpdateScheduler, int numberOfTimeSteps);
- virtual ~DemCalculator() {}
-
- void calculate() override;
-
-};
-
-#endif
diff --git a/source/DemCoupling/DemCoupling.cmake b/source/DemCoupling/DemCoupling.cmake
index b8480ee51bb32e8a1da2e668240482018454f42f..12762704e3cdbeaf03d5feea5d260918cc2c4869 100644
--- a/source/DemCoupling/DemCoupling.cmake
+++ b/source/DemCoupling/DemCoupling.cmake
@@ -18,4 +18,6 @@ SET(CAB_ADDITIONAL_LINK_LIBRARIES ${CAB_ADDITIONAL_LINK_LIBRARIES} ${LINK_LIBRAR
SET(LINK_LIBRARY optimized ${CORE_RELEASE_LIBRARY} debug ${CORE_DEBUG_LIBRARY})
SET(CAB_ADDITIONAL_LINK_LIBRARIES ${CAB_ADDITIONAL_LINK_LIBRARIES} ${LINK_LIBRARY})
-
+IF(${USE_GCC})
+ SET(CAB_ADDITIONAL_LINK_LIBRARIES ${CAB_ADDITIONAL_LINK_LIBRARIES} "stdc++fs")
+ENDIF()
diff --git a/source/DemCoupling/MovableObjectInteractor.cpp b/source/DemCoupling/MovableObjectInteractor.cpp
index 29e3409278e5d820afc0763e58780998d26e3fa7..91156903a42be156e3ab00b3e4ee8db063151103 100644
--- a/source/DemCoupling/MovableObjectInteractor.cpp
+++ b/source/DemCoupling/MovableObjectInteractor.cpp
@@ -11,7 +11,7 @@
#include "BCProcessor.h"
#include "ILBMKernel.h"
-#include "SetSolidOrBoundaryBlockVisitor.h"
+#include "SetBcBlocksBlockVisitor.h"
#include "BoundaryConditionsBlockVisitor.h"
#include "PhysicsEngineGeometryAdapter.h"
@@ -116,7 +116,7 @@ void MovableObjectInteractor::setBcNodesToFluid()
void MovableObjectInteractor::setBcBlocks()
{
- SetSolidOrBoundaryBlockVisitor v(shared_from_this(), SetSolidOrBoundaryBlockVisitor::BC);
+ SetBcBlocksBlockVisitor v(shared_from_this());
this->grid.lock()->accept(v);
}
diff --git a/source/DemCoupling/PePartitioningGridVisitor.cpp b/source/DemCoupling/PePartitioningGridVisitor.cpp
index b831edd04a8c9c04b83e618413ed49745436831b..df05d5ef2aa9241ddeec6ebee5d506d7304c2663 100644
--- a/source/DemCoupling/PePartitioningGridVisitor.cpp
+++ b/source/DemCoupling/PePartitioningGridVisitor.cpp
@@ -64,59 +64,59 @@ void PePartitioningGridVisitor::collectData(SPtr<Grid3D> grid)
}
}
//////////////////////////////////////////////////////////////////////////
-void PePartitioningGridVisitor::getBlocksByCuboid(double minX1, double minX2, double minX3, double maxX1, double maxX2, double maxX3, std::vector<SPtr<Block3D>>& blocks, SPtr<Grid3D> grid)
-{
- int coarsestLevel = grid->getCoarsestInitializedLevel();
- int finestLevel = grid->getFinestInitializedLevel();
-
- SPtr<CoordinateTransformation3D> trafo = grid->getCoordinateTransformator();
-
- //////////////////////////////////////////////////////////////////////////
- //MINIMALE BLOCK-INDIZES BESTIMMEN
- //
- //min:
- double dMinX1 = trafo->transformForwardToX1Coordinate(minX1, minX2, minX3)*(1<<finestLevel);
- double dMinX2 = trafo->transformForwardToX2Coordinate(minX1, minX2, minX3)*(1<<finestLevel);
- double dMinX3 = trafo->transformForwardToX3Coordinate(minX1, minX2, minX3)*(1<<finestLevel);
-
- //Achtung, wenn minX1 genau auf grenze zwischen zwei bloecken -> der "kleinere" muss genommen werden,
- //da beim Transformieren der "groessere" Index rauskommt
- int iMinX1 = (int)dMinX1; //if (UbMath::zero(dMinX1-iMinX1)) iMinX1-=1;
- int iMinX2 = (int)dMinX2; //if (UbMath::zero(dMinX2-iMinX2)) iMinX2-=1;
- int iMinX3 = (int)dMinX3; //if (UbMath::zero(dMinX3-iMinX3)) iMinX3-=1;
-
- //max (hier kann die Zusatzabfrage vernachlaessigt werden):
- int iMaxX1 = (int)(trafo->transformForwardToX1Coordinate(maxX1, maxX2, maxX3)*(1<<finestLevel));
- int iMaxX2 = (int)(trafo->transformForwardToX2Coordinate(maxX1, maxX2, maxX3)*(1<<finestLevel));
- int iMaxX3 = (int)(trafo->transformForwardToX3Coordinate(maxX1, maxX2, maxX3)*(1<<finestLevel));
-
- SPtr<Block3D> block;
-
- //set, um doppelte bloecke zu vermeiden, die u.U. bei periodic auftreten koennen
- std::set<SPtr<Block3D>> blockset;
- for (int level=coarsestLevel; level<=finestLevel; level++)
- {
- //damit bei negativen werten auch der "kleinere" genommen wird -> floor!
- int minx1 = (int)std::floor((double)iMinX1/(1<<(finestLevel-level)));
- int minx2 = (int)std::floor((double)iMinX2/(1<<(finestLevel-level)));
- int minx3 = (int)std::floor((double)iMinX3/(1<<(finestLevel-level)));
-
- int maxx1 = iMaxX1/(1<<(finestLevel-level));
- int maxx2 = iMaxX2/(1<<(finestLevel-level));
- int maxx3 = iMaxX3/(1<<(finestLevel-level));
-
- for (int ix1=minx1; ix1<maxx1; ix1++)
- for (int ix2=minx2; ix2<maxx2; ix2++)
- for (int ix3=minx3; ix3<maxx3; ix3++)
- if ((block=grid->getBlock(ix1, ix2, ix3, level)))
- {
- blockset.insert(block);
- }
- }
-
- blocks.resize(blockset.size());
- std::copy(blockset.begin(), blockset.end(), blocks.begin());
-}
+//void PePartitioningGridVisitor::getBlocksByCuboid(double minX1, double minX2, double minX3, double maxX1, double maxX2, double maxX3, std::vector<SPtr<Block3D>>& blocks, SPtr<Grid3D> grid)
+//{
+// int coarsestLevel = grid->getCoarsestInitializedLevel();
+// int finestLevel = grid->getFinestInitializedLevel();
+//
+// SPtr<CoordinateTransformation3D> trafo = grid->getCoordinateTransformator();
+//
+// //////////////////////////////////////////////////////////////////////////
+// //MINIMALE BLOCK-INDIZES BESTIMMEN
+// //
+// //min:
+// double dMinX1 = trafo->transformForwardToX1Coordinate(minX1, minX2, minX3)*(1<<finestLevel);
+// double dMinX2 = trafo->transformForwardToX2Coordinate(minX1, minX2, minX3)*(1<<finestLevel);
+// double dMinX3 = trafo->transformForwardToX3Coordinate(minX1, minX2, minX3)*(1<<finestLevel);
+//
+// //Achtung, wenn minX1 genau auf grenze zwischen zwei bloecken -> der "kleinere" muss genommen werden,
+// //da beim Transformieren der "groessere" Index rauskommt
+// int iMinX1 = (int)dMinX1; //if (UbMath::zero(dMinX1-iMinX1)) iMinX1-=1;
+// int iMinX2 = (int)dMinX2; //if (UbMath::zero(dMinX2-iMinX2)) iMinX2-=1;
+// int iMinX3 = (int)dMinX3; //if (UbMath::zero(dMinX3-iMinX3)) iMinX3-=1;
+//
+// //max (hier kann die Zusatzabfrage vernachlaessigt werden):
+// int iMaxX1 = (int)(trafo->transformForwardToX1Coordinate(maxX1, maxX2, maxX3)*(1<<finestLevel));
+// int iMaxX2 = (int)(trafo->transformForwardToX2Coordinate(maxX1, maxX2, maxX3)*(1<<finestLevel));
+// int iMaxX3 = (int)(trafo->transformForwardToX3Coordinate(maxX1, maxX2, maxX3)*(1<<finestLevel));
+//
+// SPtr<Block3D> block;
+//
+// //set, um doppelte bloecke zu vermeiden, die u.U. bei periodic auftreten koennen
+// std::set<SPtr<Block3D>> blockset;
+// for (int level=coarsestLevel; level<=finestLevel; level++)
+// {
+// //damit bei negativen werten auch der "kleinere" genommen wird -> floor!
+// int minx1 = (int)std::floor((double)iMinX1/(1<<(finestLevel-level)));
+// int minx2 = (int)std::floor((double)iMinX2/(1<<(finestLevel-level)));
+// int minx3 = (int)std::floor((double)iMinX3/(1<<(finestLevel-level)));
+//
+// int maxx1 = iMaxX1/(1<<(finestLevel-level));
+// int maxx2 = iMaxX2/(1<<(finestLevel-level));
+// int maxx3 = iMaxX3/(1<<(finestLevel-level));
+//
+// for (int ix1=minx1; ix1<maxx1; ix1++)
+// for (int ix2=minx2; ix2<maxx2; ix2++)
+// for (int ix3=minx3; ix3<maxx3; ix3++)
+// if ((block=grid->getBlock(ix1, ix2, ix3, level)))
+// {
+// blockset.insert(block);
+// }
+// }
+//
+// blocks.resize(blockset.size());
+// std::copy(blockset.begin(), blockset.end(), blocks.begin());
+//}
SPtr<Block3D> PePartitioningGridVisitor::getBlockByMinUniform(double minX1, double minX2, double minX3, SPtr<Grid3D> grid)
{
@@ -146,7 +146,7 @@ void PePartitioningGridVisitor::distributePartitionData(SPtr<Grid3D> grid)
for (int i = 0; i < totalIDs.size(); i++)
{
SPtr<Block3D> block = grid->getBlock(totalIDs[i]);
- block->setRank(totalRanks[i]);
+ if (block) block->setRank(totalRanks[i]);
}
}
//////////////////////////////////////////////////////////////////////////
diff --git a/source/DemCoupling/PePartitioningGridVisitor.h b/source/DemCoupling/PePartitioningGridVisitor.h
index df00c8a9ea9bf459b6bb91920f77cb9e9d43bc94..bce2ca0d2059b4c2915ab30fc4baa87b09cf9bf8 100644
--- a/source/DemCoupling/PePartitioningGridVisitor.h
+++ b/source/DemCoupling/PePartitioningGridVisitor.h
@@ -13,8 +13,8 @@
#include <array>
////////////////////////////////////////////////////////////////////////
-//! \brief The class implements domain decomposition with METIS library
-//! \author Kostyantyn Kucher
+//! \brief The class implements domain decomposition with PE library
+//! \author Konstantin Kutscher
//////////////////////////////////////////////////////////////////////////
class Communicator;
class Grid3D;
@@ -39,7 +39,7 @@ public:
protected:
void collectData(SPtr<Grid3D> grid);
void distributePartitionData(SPtr<Grid3D> grid);
- void getBlocksByCuboid(double minX1, double minX2, double minX3, double maxX1, double maxX2, double maxX3, std::vector<SPtr<Block3D>>& blocks, SPtr<Grid3D> grid);
+ //void getBlocksByCuboid(double minX1, double minX2, double minX3, double maxX1, double maxX2, double maxX3, std::vector<SPtr<Block3D>>& blocks, SPtr<Grid3D> grid);
SPtr<Block3D> getBlockByMinUniform(double minX1, double minX2, double minX3, SPtr<Grid3D> grid);
private:
diff --git a/source/VirtualFluidsCore/CoProcessors/WriteDemObjectsCoProcessor.cpp b/source/VirtualFluidsCore/CoProcessors/WriteDemObjectsCoProcessor.cpp
deleted file mode 100644
index 5bd5a95f42ce18463354d0ba00ac5f2e776ec068..0000000000000000000000000000000000000000
--- a/source/VirtualFluidsCore/CoProcessors/WriteDemObjectsCoProcessor.cpp
+++ /dev/null
@@ -1,82 +0,0 @@
-#include "WriteDemObjectsCoProcessor.h"
-#include "LBMKernel.h"
-#include "BCProcessor.h"
-
-#include "basics/writer/WbWriterVtkXmlBinary.h"
-#include "basics/writer/WbWriterVtkXmlASCII.h"
-
-#include "GbSphere3D.h"
-#include "Communicator.h"
-#include "UbScheduler.h"
-#include "Grid3D.h"
-#include "UbSystem.h"
-
-WriteDemObjectsCoProcessor::WriteDemObjectsCoProcessor()
-{
-
-}
-//////////////////////////////////////////////////////////////////////////
-WriteDemObjectsCoProcessor::WriteDemObjectsCoProcessor(SPtr<Grid3D> grid, SPtr<UbScheduler> s, const std::string& path, SPtr<DemCoProcessor> demCoProcessor, SPtr<Communicator> comm)
- : CoProcessor(grid, s),
- path(path),
- demCoProcessor(demCoProcessor),
- comm(comm)
-{
- this->path += "/geo/objects/sphere";
-}
-//////////////////////////////////////////////////////////////////////////
-void WriteDemObjectsCoProcessor::process(double step)
-{
- if (scheduler->isDue(step))
- {
- std::vector<UbTupleFloat3> nodes;
- std::vector<UbTupleInt3> triangles;
-
- for (size_t i = 0; i < objects.size(); i++)
- {
- objects[i]->addSurfaceTriangleSet(nodes, triangles);
-
-
- }
- int stepInt = (int)step;
- std::string outFilename = WbWriterVtkXmlBinary::getInstance()->writeTriangles(path + UbSystem::toString(stepInt), nodes, triangles);
-
-
- std::string pfilePath, partPath, subfolder, cfilePath;
-
- subfolder = "mq"+UbSystem::toString(istep);
- pfilePath = path+"/mq/"+subfolder;
- cfilePath = path+"/mq/mq_collection";
- partPath = pfilePath+"/mq"+UbSystem::toString(gridRank)+ "_" + UbSystem::toString(istep);
-
-
- std::string partName = writer->writeOctsWithNodeData(partPath, nodes, cells, datanames, data);
- size_t found=partName.find_last_of("/");
- std::string piece = partName.substr(found+1);
- piece = subfolder + "/" + piece;
-
- std::vector<std::string> cellDataNames;
- std::vector<std::string> pieces = comm->gather(piece);
- if (comm->getProcessID() == comm->getRoot())
- {
- std::string pname = WbWriterVtkXmlASCII::getInstance()->writeParallelFile(pfilePath, pieces, datanames, cellDataNames);
- found=pname.find_last_of("/");
- piece = pname.substr(found+1);
-
- std::vector<std::string> filenames;
- filenames.push_back(piece);
- if (step == CoProcessor::scheduler->getMinBegin())
- {
- WbWriterVtkXmlASCII::getInstance()->writeCollection(cfilePath, filenames, istep, false);
- }
- else
- {
- WbWriterVtkXmlASCII::getInstance()->addFilesToCollection(cfilePath, filenames, istep, false);
- }
- UBLOG(logINFO, "WriteMacroscopicQuantitiesCoProcessor step: " << istep);
- }
-
- }
-
-
-}
diff --git a/source/VirtualFluidsCore/Grid/BasicCalculator.cpp b/source/VirtualFluidsCore/Grid/BasicCalculator.cpp
index 2b08c7221354549c2cf6a6a537238c64c90a0e00..b07fde93c2089420e102aabb1c679a60acce5ff2 100644
--- a/source/VirtualFluidsCore/Grid/BasicCalculator.cpp
+++ b/source/VirtualFluidsCore/Grid/BasicCalculator.cpp
@@ -1,336 +1,366 @@
-#include "BasicCalculator.h"
-
-#include "Block3D.h"
-#include "BCProcessor.h"
-#include "LBMKernel.h"
-#include "Block3DConnector.h"
-#include "UbScheduler.h"
-#include "UbLogger.h"
-
-#ifdef _OPENMP
-#include <omp.h>
-#endif
-#define OMP_SCHEDULE guided
-
-//#define TIMING
-//#include "UbTiming.h"
-
-BasicCalculator::BasicCalculator(SPtr<Grid3D> grid, SPtr<UbScheduler> additionalGhostLayerUpdateScheduler, int numberOfTimeSteps) :
- Calculator(grid, additionalGhostLayerUpdateScheduler, numberOfTimeSteps)
-{
-
-}
-//////////////////////////////////////////////////////////////////////////
-BasicCalculator::~BasicCalculator()
-{
-
-}
-//////////////////////////////////////////////////////////////////////////
-void BasicCalculator::calculate()
-{
- UBLOG(logDEBUG1, "OMPCalculator::calculate() - started");
- int calcStep = 0;
- try
- {
- int minInitLevel = minLevel;
- int maxInitLevel = maxLevel-minLevel;
- int straightStartLevel = minInitLevel;
- int internalIterations = 1<<(maxInitLevel-minInitLevel);
- int forwardStartLevel;
- int threshold;
-
-#ifdef TIMING
- UbTimer timer;
- double time[6];
-#endif
-
- for (calcStep = startTimeStep; calcStep<=numberOfTimeSteps; calcStep++)
- {
- //////////////////////////////////////////////////////////////////////////
-#ifdef TIMING
- UBLOG(logINFO, "calcStep = "<<calcStep);
-#endif
- //////////////////////////////////////////////////////////////////////////
-
- for (int staggeredStep = 1; staggeredStep<=internalIterations; staggeredStep++)
- {
- forwardStartLevel = straightStartLevel;
- if (staggeredStep==internalIterations) straightStartLevel = minInitLevel;
- else
- {
- for (straightStartLevel = maxInitLevel, threshold = 1;
- (staggeredStep&threshold)!=threshold; straightStartLevel--, threshold <<= 1);
- }
- //////////////////////////////////////////////////////////////////////////
-#ifdef TIMING
- timer.resetAndStart();
-#endif
- //////////////////////////////////////////////////////////////////////////
- applyPreCollisionBC(straightStartLevel, maxInitLevel);
-
- //do collision for all blocks
- calculateBlocks(straightStartLevel, maxInitLevel, calcStep);
- //////////////////////////////////////////////////////////////////////////
-#ifdef TIMING
- time[0] = timer.stop();
- UBLOG(logINFO, "calculateBlocks time = "<<time[0]);
-#endif
- //////////////////////////////////////////////////////////////////////////
- //////////////////////////////////////////////////////////////////////////
- //exchange data between blocks
- exchangeBlockData(straightStartLevel, maxInitLevel);
- //////////////////////////////////////////////////////////////////////////
-#ifdef TIMING
- time[1] = timer.stop();
- UBLOG(logINFO, "exchangeBlockData time = "<<time[1]);
-#endif
- //////////////////////////////////////////////////////////////////////////
- applyPostCollisionBC(straightStartLevel, maxInitLevel);
- //////////////////////////////////////////////////////////////////////////
-#ifdef TIMING
- time[2] = timer.stop();
- UBLOG(logINFO, "applyBCs time = "<<time[2]);
-#endif
- //////////////////////////////////////////////////////////////////////////
- //swap distributions in kernel
- swapDistributions(straightStartLevel, maxInitLevel);
- //////////////////////////////////////////////////////////////////////////
-#ifdef TIMING
- time[3] = timer.stop();
- UBLOG(logINFO, "swapDistributions time = "<<time[3]);
-#endif
- //////////////////////////////////////////////////////////////////////////
- if (refinement)
- {
- if (straightStartLevel<maxInitLevel)
- exchangeBlockData(straightStartLevel, maxInitLevel);
- //////////////////////////////////////////////////////////////////////////
-#ifdef TIMING
- time[4] = timer.stop();
- UBLOG(logINFO, "refinement exchangeBlockData time = "<<time[4]);
-#endif
- //////////////////////////////////////////////////////////////////////////
- //now ghost nodes have actual values
- //interpolation of interface nodes between grid levels
- interpolation(straightStartLevel, maxInitLevel);
- //////////////////////////////////////////////////////////////////////////
-#ifdef TIMING
- time[5] = timer.stop();
- UBLOG(logINFO, "refinement interpolation time = "<<time[5]);
-#endif
- //////////////////////////////////////////////////////////////////////////
- }
- }
- //exchange data between blocks for visualization
- if (additionalGhostLayerUpdateScheduler->isDue(calcStep))
- {
- exchangeBlockData(straightStartLevel, maxInitLevel);
- }
- coProcess((double)(calcStep));
- //now ghost nodes have actual values
- }
- UBLOG(logDEBUG1, "OMPCalculator::calculate() - stoped");
- }
- catch (std::exception& e)
- {
- UBLOG(logERROR, e.what());
- UBLOG(logERROR, " step = "<<calcStep);
- //throw;
- exit(EXIT_FAILURE);
- }
-}
-//////////////////////////////////////////////////////////////////////////
-void BasicCalculator::calculateBlocks(int startLevel, int maxInitLevel, int calcStep)
-{
-#ifdef _OPENMP
-#pragma omp parallel
-#endif
- {
- SPtr<Block3D> blockTemp;
- try
- {
- //startLevel bis maxInitLevel
- for (int level = startLevel; level<=maxInitLevel; level++)
- {
- //timer.resetAndStart();
- //call LBM kernel
- int size = (int)blocks[level].size();
-#ifdef _OPENMP
-#pragma omp for schedule(OMP_SCHEDULE)
-#endif
- for (int i =0; i<size; i++)
- {
- blockTemp = blocks[level][i];
- blockTemp->getKernel()->calculate(calcStep);
- }
- //timer.stop();
- //UBLOG(logINFO, "level = " << level << " blocks = " << blocks[level].size() << " collision time = " << timer.getTotalTime());
- }
- }
- catch (std::exception& e)
- {
- UBLOG(logERROR, e.what());
- //UBLOG(logERROR, blockTemp->toString()<<" step = "<<calcStep);
- //throw;
- exit(EXIT_FAILURE);
- }
- }
-}
-//////////////////////////////////////////////////////////////////////////
-void BasicCalculator::exchangeBlockData(int startLevel, int maxInitLevel)
-{
- //startLevel bis maxInitLevel
- for (int level = startLevel; level<=maxInitLevel; level++)
- {
- //connectorsPrepareLocal(localConns[level]);
- connectorsSendLocal(localConns[level]);
- //connectorsReceiveLocal(localConns[level]);
-
- connectorsPrepareRemote(remoteConns[level]);
- connectorsSendRemote(remoteConns[level]);
- connectorsReceiveRemote(remoteConns[level]);
- }
-}
-//////////////////////////////////////////////////////////////////////////
-void BasicCalculator::swapDistributions(int startLevel, int maxInitLevel)
-{
-#ifdef _OPENMP
-#pragma omp parallel
-#endif
- {
- //startLevel bis maxInitLevel
- for (int level = startLevel; level<=maxInitLevel; level++)
- {
- int size = (int)blocks[level].size();
-#ifdef _OPENMP
-#pragma omp for schedule(OMP_SCHEDULE)
-#endif
- for (int i =0; i<size; i++)
- {
- blocks[level][i]->getKernel()->swapDistributions();
- }
- }
- }
-}
-//////////////////////////////////////////////////////////////////////////
-void BasicCalculator::connectorsPrepareLocal(std::vector< SPtr<Block3DConnector> >& connectors)
-{
- int size = (int)connectors.size();
-#ifdef _OPENMP
-#pragma omp parallel for schedule(OMP_SCHEDULE)
-#endif
- for (int i =0; i<size; i++)
- {
- connectors[i]->prepareForReceive();
- connectors[i]->prepareForSend();
- }
-}
-//////////////////////////////////////////////////////////////////////////
-void BasicCalculator::connectorsSendLocal(std::vector< SPtr<Block3DConnector> >& connectors)
-{
- int size = (int)connectors.size();
-#ifdef _OPENMP
-#pragma omp parallel for schedule(OMP_SCHEDULE)
-#endif
- for (int i =0; i<size; i++)
- {
- connectors[i]->fillSendVectors();
- connectors[i]->sendVectors();
- }
-}
-//////////////////////////////////////////////////////////////////////////
-void BasicCalculator::connectorsReceiveLocal(std::vector< SPtr<Block3DConnector> >& connectors)
-{
- int size = (int)connectors.size();
-#ifdef _OPENMP
-#pragma omp parallel for schedule(OMP_SCHEDULE)
-#endif
- for (int i =0; i<size; i++)
- {
- connectors[i]->receiveVectors();
- connectors[i]->distributeReceiveVectors();
- }
-}
-void BasicCalculator::connectorsPrepareRemote(std::vector< SPtr<Block3DConnector> >& connectors)
-{
- int size = (int)connectors.size();
- for (int i =0; i<size; i++)
- {
- connectors[i]->prepareForReceive();
- connectors[i]->prepareForSend();
- }
-}
-//////////////////////////////////////////////////////////////////////////
-void BasicCalculator::connectorsSendRemote(std::vector< SPtr<Block3DConnector> >& connectors)
-{
- int size = (int)connectors.size();
- for (int i =0; i<size; i++)
- {
- connectors[i]->fillSendVectors();
- connectors[i]->sendVectors();
- }
-}
-//////////////////////////////////////////////////////////////////////////
-void BasicCalculator::connectorsReceiveRemote(std::vector< SPtr<Block3DConnector> >& connectors)
-{
- int size = (int)connectors.size();
- for (int i =0; i<size; i++)
- {
- connectors[i]->receiveVectors();
- connectors[i]->distributeReceiveVectors();
- }
-}
-//////////////////////////////////////////////////////////////////////////
-void BasicCalculator::interpolation(int startLevel, int maxInitLevel)
-{
- for (int level = startLevel; level<maxInitLevel; level++)
- {
- connectorsPrepareLocal(localInterConns[level]);
- connectorsPrepareRemote(remoteInterConns[level]);
- }
-
- for (int level = startLevel; level<maxInitLevel; level++)
- {
- connectorsSendLocal(localInterConns[level]);
- connectorsSendRemote(remoteInterConns[level]);
- }
-
- for (int level = startLevel; level<maxInitLevel; level++)
- {
- connectorsReceiveLocal(localInterConns[level]);
- connectorsReceiveRemote(remoteInterConns[level]);
- }
-}
-//////////////////////////////////////////////////////////////////////////
-void BasicCalculator::applyPreCollisionBC(int startLevel, int maxInitLevel)
-{
- //startLevel bis maxInitLevel
- for (int level = startLevel; level<=maxInitLevel; level++)
- {
- int size = (int)blocks[level].size();
-#ifdef _OPENMP
-#pragma omp parallel for schedule(OMP_SCHEDULE)
-#endif
- for (int i =0; i<size; i++)
- {
- blocks[level][i]->getKernel()->getBCProcessor()->applyPreCollisionBC();
- }
- }
-}
-//////////////////////////////////////////////////////////////////////////
-void BasicCalculator::applyPostCollisionBC(int startLevel, int maxInitLevel)
-{
- //startLevel bis maxInitLevel
- for (int level = startLevel; level<=maxInitLevel; level++)
- {
- int size = (int)blocks[level].size();
-#ifdef _OPENMP
-#pragma omp parallel for schedule(OMP_SCHEDULE)
-#endif
- for (int i =0; i<size; i++)
- {
- blocks[level][i]->getKernel()->getBCProcessor()->applyPostCollisionBC();
- }
- }
-}
-
+#include "BasicCalculator.h"
+
+#include "Block3D.h"
+#include "BCProcessor.h"
+#include "LBMKernel.h"
+#include "Block3DConnector.h"
+#include "UbScheduler.h"
+#include "UbLogger.h"
+
+#ifdef _OPENMP
+#include <omp.h>
+#endif
+#define OMP_SCHEDULE guided
+
+//#define TIMING
+//#include "UbTiming.h"
+
+BasicCalculator::BasicCalculator(SPtr<Grid3D> grid, SPtr<UbScheduler> additionalGhostLayerUpdateScheduler, int numberOfTimeSteps) :
+ Calculator(grid, additionalGhostLayerUpdateScheduler, numberOfTimeSteps)
+{
+
+}
+//////////////////////////////////////////////////////////////////////////
+BasicCalculator::~BasicCalculator()
+{
+
+}
+//////////////////////////////////////////////////////////////////////////
+void BasicCalculator::calculate()
+{
+ UBLOG(logDEBUG1, "OMPCalculator::calculate() - started");
+ int calcStep = 0;
+ try
+ {
+ int minInitLevel = minLevel;
+ int maxInitLevel = maxLevel-minLevel;
+ int straightStartLevel = minInitLevel;
+ int internalIterations = 1<<(maxInitLevel-minInitLevel);
+ int forwardStartLevel;
+ int threshold;
+
+#ifdef TIMING
+ UbTimer timer;
+ double time[6];
+#endif
+
+ for (calcStep = startTimeStep; calcStep<=numberOfTimeSteps; calcStep++)
+ {
+ //////////////////////////////////////////////////////////////////////////
+#ifdef TIMING
+ UBLOG(logINFO, "calcStep = "<<calcStep);
+#endif
+ //////////////////////////////////////////////////////////////////////////
+
+ for (int staggeredStep = 1; staggeredStep<=internalIterations; staggeredStep++)
+ {
+ forwardStartLevel = straightStartLevel;
+ if (staggeredStep==internalIterations) straightStartLevel = minInitLevel;
+ else
+ {
+ for (straightStartLevel = maxInitLevel, threshold = 1;
+ (staggeredStep&threshold)!=threshold; straightStartLevel--, threshold <<= 1);
+ }
+ //////////////////////////////////////////////////////////////////////////
+#ifdef TIMING
+ timer.resetAndStart();
+#endif
+ //////////////////////////////////////////////////////////////////////////
+ applyPreCollisionBC(straightStartLevel, maxInitLevel);
+
+ //do collision for all blocks
+ calculateBlocks(straightStartLevel, maxInitLevel, calcStep);
+ //////////////////////////////////////////////////////////////////////////
+#ifdef TIMING
+ time[0] = timer.stop();
+ UBLOG(logINFO, "calculateBlocks time = "<<time[0]);
+#endif
+ //////////////////////////////////////////////////////////////////////////
+ //////////////////////////////////////////////////////////////////////////
+ //exchange data between blocks
+ exchangeBlockData(straightStartLevel, maxInitLevel);
+ //////////////////////////////////////////////////////////////////////////
+#ifdef TIMING
+ time[1] = timer.stop();
+ UBLOG(logINFO, "exchangeBlockData time = "<<time[1]);
+#endif
+ //////////////////////////////////////////////////////////////////////////
+ applyPostCollisionBC(straightStartLevel, maxInitLevel);
+ //////////////////////////////////////////////////////////////////////////
+#ifdef TIMING
+ time[2] = timer.stop();
+ UBLOG(logINFO, "applyBCs time = "<<time[2]);
+#endif
+ //////////////////////////////////////////////////////////////////////////
+ //swap distributions in kernel
+ swapDistributions(straightStartLevel, maxInitLevel);
+ //////////////////////////////////////////////////////////////////////////
+#ifdef TIMING
+ time[3] = timer.stop();
+ UBLOG(logINFO, "swapDistributions time = "<<time[3]);
+#endif
+ //////////////////////////////////////////////////////////////////////////
+ if (refinement)
+ {
+ if (straightStartLevel<maxInitLevel)
+ exchangeBlockData(straightStartLevel, maxInitLevel);
+ //////////////////////////////////////////////////////////////////////////
+#ifdef TIMING
+ time[4] = timer.stop();
+ UBLOG(logINFO, "refinement exchangeBlockData time = "<<time[4]);
+#endif
+ //////////////////////////////////////////////////////////////////////////
+ //now ghost nodes have actual values
+ //interpolation of interface nodes between grid levels
+ interpolation(straightStartLevel, maxInitLevel);
+ //////////////////////////////////////////////////////////////////////////
+#ifdef TIMING
+ time[5] = timer.stop();
+ UBLOG(logINFO, "refinement interpolation time = "<<time[5]);
+#endif
+ //////////////////////////////////////////////////////////////////////////
+ }
+ }
+ //exchange data between blocks for visualization
+ if (additionalGhostLayerUpdateScheduler->isDue(calcStep))
+ {
+ exchangeBlockData(straightStartLevel, maxInitLevel);
+ }
+ coProcess((double)(calcStep));
+ //now ghost nodes have actual values
+ }
+ UBLOG(logDEBUG1, "OMPCalculator::calculate() - stoped");
+ }
catch (std::exception& e)
+ {
+ UBLOG(logERROR, e.what());
+ UBLOG(logERROR, " step = "<<calcStep);
+ //throw;
+ exit(EXIT_FAILURE);
+ }
+ catch (std::string& s)
+ {
+ UBLOG(logERROR, s);
+ exit(EXIT_FAILURE);
+ }
+ catch (...)
+ {
+ UBLOG(logERROR, "unknown exception");
+ exit(EXIT_FAILURE);
+ }
+}
+//////////////////////////////////////////////////////////////////////////
+void BasicCalculator::calculateBlocks(int startLevel, int maxInitLevel, int calcStep)
+{
+#ifdef _OPENMP
+#pragma omp parallel
+#endif
+ {
+ SPtr<Block3D> blockTemp;
+ try
+ {
+ //startLevel bis maxInitLevel
+ for (int level = startLevel; level<=maxInitLevel; level++)
+ {
+ //timer.resetAndStart();
+ //call LBM kernel
+ int size = (int)blocks[level].size();
+#ifdef _OPENMP
+#pragma omp for schedule(OMP_SCHEDULE)
+#endif
+ for (int i =0; i<size; i++)
+ {
+ blockTemp = blocks[level][i];
+ blockTemp->getKernel()->calculate(calcStep);
+ }
+ //timer.stop();
+ //UBLOG(logINFO, "level = " << level << " blocks = " << blocks[level].size() << " collision time = " << timer.getTotalTime());
+ }
+ }
+ catch (std::exception& e)
+ {
+ UBLOG(logERROR, e.what());
+ //UBLOG(logERROR, blockTemp->toString()<<" step = "<<calcStep);
+ //throw;
+ exit(EXIT_FAILURE);
+ }
+ }
+}
+//////////////////////////////////////////////////////////////////////////
+void BasicCalculator::exchangeBlockData(int startLevel, int maxInitLevel)
+{
+ //startLevel bis maxInitLevel
+ for (int level = startLevel; level<=maxInitLevel; level++)
+ {
+ //connectorsPrepareLocal(localConns[level]);
+ connectorsSendLocal(localConns[level]);
+ //connectorsReceiveLocal(localConns[level]);
+
+ connectorsPrepareRemote(remoteConns[level]);
+ connectorsSendRemote(remoteConns[level]);
+ connectorsReceiveRemote(remoteConns[level]);
+ }
+}
+//////////////////////////////////////////////////////////////////////////
+void BasicCalculator::swapDistributions(int startLevel, int maxInitLevel)
+{
+#ifdef _OPENMP
+#pragma omp parallel
+#endif
+ {
+ //startLevel bis maxInitLevel
+ for (int level = startLevel; level<=maxInitLevel; level++)
+ {
+ int size = (int)blocks[level].size();
+#ifdef _OPENMP
+#pragma omp for schedule(OMP_SCHEDULE)
+#endif
+ for (int i =0; i<size; i++)
+ {
+ blocks[level][i]->getKernel()->swapDistributions();
+ }
+ }
+ }
+}
+//////////////////////////////////////////////////////////////////////////
+void BasicCalculator::connectorsPrepareLocal(std::vector< SPtr<Block3DConnector> >& connectors)
+{
+ int size = (int)connectors.size();
+#ifdef _OPENMP
+#pragma omp parallel for schedule(OMP_SCHEDULE)
+#endif
+ for (int i =0; i<size; i++)
+ {
+ connectors[i]->prepareForReceive();
+ connectors[i]->prepareForSend();
+ }
+}
+//////////////////////////////////////////////////////////////////////////
+void BasicCalculator::connectorsSendLocal(std::vector< SPtr<Block3DConnector> >& connectors)
+{
+ int size = (int)connectors.size();
+#ifdef _OPENMP
+#pragma omp parallel for schedule(OMP_SCHEDULE)
+#endif
+ for (int i =0; i<size; i++)
+ {
+ connectors[i]->fillSendVectors();
+ connectors[i]->sendVectors();
+ }
+}
+//////////////////////////////////////////////////////////////////////////
+void BasicCalculator::connectorsReceiveLocal(std::vector< SPtr<Block3DConnector> >& connectors)
+{
+ int size = (int)connectors.size();
+#ifdef _OPENMP
+#pragma omp parallel for schedule(OMP_SCHEDULE)
+#endif
+ for (int i =0; i<size; i++)
+ {
+ connectors[i]->receiveVectors();
+ connectors[i]->distributeReceiveVectors();
+ }
+}
+void BasicCalculator::connectorsPrepareRemote(std::vector< SPtr<Block3DConnector> >& connectors)
+{
+ int size = (int)connectors.size();
+ for (int i =0; i<size; i++)
+ {
+ connectors[i]->prepareForReceive();
+ connectors[i]->prepareForSend();
+ }
+}
+//////////////////////////////////////////////////////////////////////////
+void BasicCalculator::connectorsSendRemote(std::vector< SPtr<Block3DConnector> >& connectors)
+{
+ int size = (int)connectors.size();
+ for (int i =0; i<size; i++)
+ {
+ connectors[i]->fillSendVectors();
+ connectors[i]->sendVectors();
+ }
+}
+//////////////////////////////////////////////////////////////////////////
+void BasicCalculator::connectorsReceiveRemote(std::vector< SPtr<Block3DConnector> >& connectors)
+{
+ int size = (int)connectors.size();
+ for (int i =0; i<size; i++)
+ {
+ connectors[i]->receiveVectors();
+ connectors[i]->distributeReceiveVectors();
+ }
+}
+//////////////////////////////////////////////////////////////////////////
+void BasicCalculator::interpolation(int startLevel, int maxInitLevel)
+{
+ for (int level = startLevel; level<maxInitLevel; level++)
+ {
+ connectorsPrepareLocal(localInterConns[level]);
+ connectorsPrepareRemote(remoteInterConns[level]);
+ }
+
+ for (int level = startLevel; level<maxInitLevel; level++)
+ {
+ connectorsSendLocal(localInterConns[level]);
+ connectorsSendRemote(remoteInterConns[level]);
+ }
+
+ for (int level = startLevel; level<maxInitLevel; level++)
+ {
+ connectorsReceiveLocal(localInterConns[level]);
+ connectorsReceiveRemote(remoteInterConns[level]);
+ }
+}
+//////////////////////////////////////////////////////////////////////////
+void BasicCalculator::applyPreCollisionBC(int startLevel, int maxInitLevel)
+{
+ //startLevel bis maxInitLevel
+ for (int level = startLevel; level<=maxInitLevel; level++)
+ {
+ int size = (int)blocks[level].size();
+#ifdef _OPENMP
+#pragma omp parallel for schedule(OMP_SCHEDULE)
+#endif
+ for (int i =0; i<size; i++)
+ {
+ blocks[level][i]->getKernel()->getBCProcessor()->applyPreCollisionBC();
+ }
+ }
+}
+//////////////////////////////////////////////////////////////////////////
+void BasicCalculator::applyPostCollisionBC(int startLevel, int maxInitLevel)
+{
+ try{
+ //startLevel bis maxInitLevel
+ for (int level = startLevel; level<=maxInitLevel; level++)
+ {
+ int size = (int)blocks[level].size();
+#ifdef _OPENMP
+#pragma omp parallel for schedule(OMP_SCHEDULE)
+#endif
+ for (int i =0; i<size; i++)
+ {
+ blocks[level][i]->getKernel()->getBCProcessor()->applyPostCollisionBC();
+ }
+ }
+}
+ catch (std::exception& e)
+ {
+ UBLOG(logERROR, e.what());
+ //UBLOG(logERROR, " step = "<<calcStep);
+ throw;
+ //exit(EXIT_FAILURE);
+ }
+ catch (std::string& s)
+ {
+ UBLOG(logERROR, s);
+ throw;
+ //exit(EXIT_FAILURE);
+ }
+ catch (...)
+ {
+ UBLOG(logERROR, "unknown exception");
+ throw;
+ //exit(EXIT_FAILURE);
+ }
+}
+
diff --git a/source/VirtualFluidsCore/Interactors/InteractorsHelper.cpp b/source/VirtualFluidsCore/Interactors/InteractorsHelper.cpp
index 38ca4ffe62b1e9a4e1a2042ea1902505a2fd9a2a..c8aea0034857a8c2c05067e55b8924c4015b295f 100644
--- a/source/VirtualFluidsCore/Interactors/InteractorsHelper.cpp
+++ b/source/VirtualFluidsCore/Interactors/InteractorsHelper.cpp
@@ -5,6 +5,8 @@
#include <Interactor3D.h>
#include "Block3D.h"
#include "Communicator.h"
+#include "SetSolidBlocksBlockVisitor.h"
+#include "SetBcBlocksBlockVisitor.h"
InteractorsHelper::InteractorsHelper(SPtr<Grid3D> grid, SPtr<Grid3DVisitor> visitor) :
@@ -48,7 +50,9 @@ void InteractorsHelper::deleteSolidBlocks()
{
for(SPtr<Interactor3D> interactor : interactors)
{
- setBlocks(interactor, SetSolidOrBoundaryBlockVisitor::SOLID);
+ //setBlocks(interactor, SetBcBlocksBlockVisitor::SOLID);
+ SetSolidBlocksBlockVisitor v(interactor);
+ grid->accept(v);
std::vector<SPtr<Block3D>>& sb = interactor->getSolidBlockSet();
solidBlocks.insert(solidBlocks.end(), sb.begin(), sb.end());
@@ -58,16 +62,13 @@ void InteractorsHelper::deleteSolidBlocks()
updateGrid();
}
//////////////////////////////////////////////////////////////////////////
-void InteractorsHelper::setBlocks(const SPtr<Interactor3D> interactor, SetSolidOrBoundaryBlockVisitor::BlockType type) const
-{
- SetSolidOrBoundaryBlockVisitor v(interactor, type);
- grid->accept(v);
-}
-//////////////////////////////////////////////////////////////////////////
void InteractorsHelper::setBcBlocks()
{
for(const SPtr<Interactor3D> interactor : interactors)
- setBlocks(interactor, SetSolidOrBoundaryBlockVisitor::BC);
+ {
+ SetBcBlocksBlockVisitor v(interactor);
+ grid->accept(v);
+ }
}
//////////////////////////////////////////////////////////////////////////
void InteractorsHelper::updateGrid()
diff --git a/source/VirtualFluidsCore/Interactors/InteractorsHelper.h b/source/VirtualFluidsCore/Interactors/InteractorsHelper.h
index d8407768ddb648713660263f62c61d362cfb3dc8..2b10c626d9b0761274a01d35667cb9b7b46dd5b0 100644
--- a/source/VirtualFluidsCore/Interactors/InteractorsHelper.h
+++ b/source/VirtualFluidsCore/Interactors/InteractorsHelper.h
@@ -3,13 +3,12 @@
#include <vector>
#include <PointerDefinitions.h>
-#include <SetSolidOrBoundaryBlockVisitor.h>
+
class Interactor3D;
class Block3D;
class Grid3D;
class Grid3DVisitor;
-class SetSolidOrBoundaryBlockVisitor;
class InteractorsHelper
{
@@ -24,7 +23,6 @@ public:
protected:
void deleteSolidBlocks();
- void setBlocks(const SPtr<Interactor3D> interactor, SetSolidOrBoundaryBlockVisitor::BlockType type) const;
void setBcBlocks();
private:
diff --git a/source/VirtualFluidsCore/Utilities/VoxelMatrixUtil.hpp b/source/VirtualFluidsCore/Utilities/VoxelMatrixUtil.hpp
index 1ef7f15b475bae82e9ff331819544cb02c66ca1f..61a58b6396d30ef9e37eb07b4a1e4a8da6202e8c 100644
--- a/source/VirtualFluidsCore/Utilities/VoxelMatrixUtil.hpp
+++ b/source/VirtualFluidsCore/Utilities/VoxelMatrixUtil.hpp
@@ -4,7 +4,7 @@
#include "GbCuboid3D.h"
#include "NoSlipBCAdapter.h"
#include "D3Q27Interactor.h"
-#include "SetSolidOrBoundaryBlockVisitor.h"
+#include "SetBcBlocksBlockVisitor.h"
#include "Block3D.h"
#include "Grid3D.h"
@@ -50,7 +50,7 @@ namespace Utilities
// vmInt->setDifferencesToGbObject3D(block);
//}
- SetSolidOrBoundaryBlockVisitor v(vmInt, SetSolidOrBoundaryBlockVisitor::BC);
+ SetBcBlocksBlockVisitor v(vmInt);
grid->accept(v);
vmInt->initInteractor();
}
diff --git a/source/VirtualFluidsCore/Visitors/SetBcBlocksBlockVisitor.cpp b/source/VirtualFluidsCore/Visitors/SetBcBlocksBlockVisitor.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..0e8fcf3404630a6c8e2499618f2d4fac9839b702
--- /dev/null
+++ b/source/VirtualFluidsCore/Visitors/SetBcBlocksBlockVisitor.cpp
@@ -0,0 +1,24 @@
+#include "SetBcBlocksBlockVisitor.h"
+
+#include "Interactor3D.h"
+#include "Grid3DSystem.h"
+#include "Grid3D.h"
+#include "Block3D.h"
+
+SetBcBlocksBlockVisitor::SetBcBlocksBlockVisitor(SPtr<Interactor3D> interactor) :
+ Block3DVisitor(0, Grid3DSystem::MAXLEVEL), interactor(interactor)
+{
+
+}
+
+void SetBcBlocksBlockVisitor::visit(SPtr<Grid3D> grid, SPtr<Block3D> block)
+{
+ if(block->getRank() == grid->getRank())
+ {
+ if (block->isActive())
+ {
+ interactor->setBCBlock(block);
+ }
+ }
+}
+
diff --git a/source/VirtualFluidsCore/Visitors/SetBcBlocksBlockVisitor.h b/source/VirtualFluidsCore/Visitors/SetBcBlocksBlockVisitor.h
new file mode 100644
index 0000000000000000000000000000000000000000..9070d5de56d8166cddde40814e2c8ffd8f28a3c3
--- /dev/null
+++ b/source/VirtualFluidsCore/Visitors/SetBcBlocksBlockVisitor.h
@@ -0,0 +1,26 @@
+#ifndef SetBcBlocksBlockVisitor_h__
+#define SetBcBlocksBlockVisitor_h__
+
+#include <PointerDefinitions.h>
+
+#include "Block3DVisitor.h"
+
+class Grid3D;
+class Block3D;
+class Interactor3D;
+
+class SetBcBlocksBlockVisitor : public Block3DVisitor
+{
+public:
+ SetBcBlocksBlockVisitor(SPtr<Interactor3D> interactor);
+ virtual ~SetBcBlocksBlockVisitor() {}
+
+ virtual void visit(SPtr<Grid3D> grid, SPtr<Block3D> block);
+
+private:
+ SPtr<Interactor3D> interactor;
+};
+#endif // SetBcBlocksBlockVisitor_h__
+
+
+
diff --git a/source/VirtualFluidsCore/Visitors/SetSolidBlocksBlockVisitor.cpp b/source/VirtualFluidsCore/Visitors/SetSolidBlocksBlockVisitor.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..7c9342f1c6c42960bfa5d9a119fd2dc3847a80fd
--- /dev/null
+++ b/source/VirtualFluidsCore/Visitors/SetSolidBlocksBlockVisitor.cpp
@@ -0,0 +1,23 @@
+#include "SetSolidBlocksBlockVisitor.h"
+
+#include "Interactor3D.h"
+#include "Grid3DSystem.h"
+#include "Grid3D.h"
+#include "Block3D.h"
+
+SetSolidBlocksBlockVisitor::SetSolidBlocksBlockVisitor(SPtr<Interactor3D> interactor) :
+ Block3DVisitor(0, Grid3DSystem::MAXLEVEL), interactor(interactor)
+{
+
+}
+
+void SetSolidBlocksBlockVisitor::visit(SPtr<Grid3D> grid, SPtr<Block3D> block)
+{
+ if (block->getRank() == grid->getRank())
+ {
+ if (block->isActive())
+ {
+ interactor->setSolidBlock(block);
+ }
+ }
+}
\ No newline at end of file
diff --git a/source/VirtualFluidsCore/Visitors/SetSolidBlocksBlockVisitor.h b/source/VirtualFluidsCore/Visitors/SetSolidBlocksBlockVisitor.h
new file mode 100644
index 0000000000000000000000000000000000000000..1f3d549125fe3da8835e3e75c2d0c5f23fc07744
--- /dev/null
+++ b/source/VirtualFluidsCore/Visitors/SetSolidBlocksBlockVisitor.h
@@ -0,0 +1,24 @@
+#ifndef SetSolidBlocksBlockVisitor_h__
+#define SetSolidBlocksBlockVisitor_h__
+
+#include <PointerDefinitions.h>
+
+#include "Block3DVisitor.h"
+
+class Grid3D;
+class Block3D;
+class Interactor3D;
+
+class SetSolidBlocksBlockVisitor : public Block3DVisitor
+{
+public:
+ SetSolidBlocksBlockVisitor(SPtr<Interactor3D> interactor);
+ virtual ~SetSolidBlocksBlockVisitor() {}
+
+ virtual void visit(SPtr<Grid3D> grid, SPtr<Block3D> block);
+
+private:
+ SPtr<Interactor3D> interactor;
+}
+#endif // SetSolidBlocksBlockVisitor_h__
+;
diff --git a/source/VirtualFluidsCore/Visitors/SetSolidOrBoundaryBlockVisitor.cpp b/source/VirtualFluidsCore/Visitors/SetSolidOrBoundaryBlockVisitor.cpp
deleted file mode 100644
index dab7b5f0443372663444fe5ffcf5612c99ad445a..0000000000000000000000000000000000000000
--- a/source/VirtualFluidsCore/Visitors/SetSolidOrBoundaryBlockVisitor.cpp
+++ /dev/null
@@ -1,33 +0,0 @@
-#include "SetSolidOrBoundaryBlockVisitor.h"
-
-#include "Interactor3D.h"
-#include "Grid3DSystem.h"
-#include "Grid3D.h"
-#include "Block3D.h"
-
-SetSolidOrBoundaryBlockVisitor::SetSolidOrBoundaryBlockVisitor(SPtr<Interactor3D> interactor, BlockType type) :
- Block3DVisitor(0, Grid3DSystem::MAXLEVEL), interactor(interactor),
- type(type)
-{
-
-}
-
-void SetSolidOrBoundaryBlockVisitor::visit(SPtr<Grid3D> grid, SPtr<Block3D> block)
-{
- if(block->getRank() == grid->getRank())
- {
- if (block->isActive())
- {
- switch (type)
- {
- case SOLID:
- interactor->setSolidBlock(block);
- break;
- case BC:
- interactor->setBCBlock(block);
- break;
- }
- }
- }
-}
-
diff --git a/source/VirtualFluidsCore/Visitors/SetSolidOrBoundaryBlockVisitor.h b/source/VirtualFluidsCore/Visitors/SetSolidOrBoundaryBlockVisitor.h
deleted file mode 100644
index d65d12be74fb78169a23fa41d7206cddbe8f3f9c..0000000000000000000000000000000000000000
--- a/source/VirtualFluidsCore/Visitors/SetSolidOrBoundaryBlockVisitor.h
+++ /dev/null
@@ -1,28 +0,0 @@
-#ifndef SET_SOLID_OR_TRANS_BLOCK_VISITOR_H
-#define SET_SOLID_OR_TRANS_BLOCK_VISITOR_H
-
-#include <PointerDefinitions.h>
-
-#include "Block3DVisitor.h"
-
-class Grid3D;
-class Block3D;
-class Interactor3D;
-
-class SetSolidOrBoundaryBlockVisitor : public Block3DVisitor
-{
-public:
- enum BlockType { SOLID, BC };
-public:
- SetSolidOrBoundaryBlockVisitor(SPtr<Interactor3D> interactor, BlockType type);
- virtual ~SetSolidOrBoundaryBlockVisitor() {}
-
- virtual void visit(SPtr<Grid3D> grid, SPtr<Block3D> block);
-
-private:
- SPtr<Interactor3D> interactor;
- BlockType type;
-};
-
-#endif
-