From c9a800c3bc0d14ae5d5553dfc6105b7070bef938 Mon Sep 17 00:00:00 2001 From: Konstantin Kutscher <kutscher@irmb.tu-bs.de> Date: Thu, 26 Apr 2018 15:08:04 +0200 Subject: [PATCH] changed exption in kernel --- source/Applications.cmake | 1 + source/Applications/BoxBenchmark/bb.cfg | 13 +- source/Applications/BoxBenchmark/bb.cpp | 55 +- .../DLR-F16-Porous/CMakeLists.txt | 25 + .../DLR-F16-Porous/f16-porous.cfg | 65 + source/Applications/DLR-F16-Porous/f16.cpp | 1155 +++++++++++++++++ .../Applications/DLR-F16-Solid/f16-solid.cfg | 33 +- source/Applications/DLR-F16-Solid/f16.cpp | 77 +- .../numerics/geometry3d/GbTriFaceMesh3D.h | 12 +- .../TimeDependentBCCoProcessor.cpp | 11 +- .../CoProcessors/TimeDependentBCCoProcessor.h | 2 +- .../WriteMQFromSelectionCoProcessor.cpp | 6 +- .../Grid/BasicCalculator.cpp | 12 +- .../D3Q27TriFaceMeshInteractor.cpp | 2 + .../Interactors/Interactor3D.h | 1 + ...ibleCumulant4thOrderViscosityLBMKernel.cpp | 9 +- ...ssibleCumulant4thOrderViscosityLBMKernel.h | 2 +- .../LBM/CompressibleCumulantLBMKernel.cpp | 11 +- .../LBM/CompressibleCumulantLBMKernel.h | 2 +- source/VirtualFluidsCore/LBM/ILBMKernel.h | 2 +- .../LBM/IncompressibleCumulantLBMKernel.cpp | 10 +- .../LBM/IncompressibleCumulantLBMKernel.h | 2 +- ...ssibleCumulantWithSpongeLayerLBMKernel.cpp | 10 +- ...ressibleCumulantWithSpongeLayerLBMKernel.h | 2 +- .../LBM/InitDensityLBMKernel.cpp | 2 +- .../LBM/InitDensityLBMKernel.h | 2 +- source/VirtualFluidsCore/LBM/LBMKernel.h | 2 +- .../VirtualFluidsCore/LBM/VoidLBMKernel.cpp | 2 +- source/VirtualFluidsCore/LBM/VoidLBMKernel.h | 2 +- .../VirtualFluidsCore/Parallel/Communicator.h | 1 + .../Visitors/SpongeLayerBlockVisitor.cpp | 2 + 31 files changed, 1424 insertions(+), 109 deletions(-) create mode 100644 source/Applications/DLR-F16-Porous/CMakeLists.txt create mode 100644 source/Applications/DLR-F16-Porous/f16-porous.cfg create mode 100644 source/Applications/DLR-F16-Porous/f16.cpp diff --git a/source/Applications.cmake b/source/Applications.cmake index 3df40f374..eb161441e 100644 --- a/source/Applications.cmake +++ b/source/Applications.cmake @@ -47,6 +47,7 @@ add_subdirectory(Applications/BoxBenchmark) add_subdirectory(Applications/DHIT) add_subdirectory(Applications/DLR-F16) add_subdirectory(Applications/DLR-F16-Solid) +add_subdirectory(Applications/DLR-F16-Porous) add_subdirectory(Applications/aperm) add_subdirectory(Applications/f16test) add_subdirectory(Applications/mirror) diff --git a/source/Applications/BoxBenchmark/bb.cfg b/source/Applications/BoxBenchmark/bb.cfg index d81b5fa13..c2abee2e2 100644 --- a/source/Applications/BoxBenchmark/bb.cfg +++ b/source/Applications/BoxBenchmark/bb.cfg @@ -1,5 +1,5 @@ pathname = d:/temp/BoxBenchmark -numOfThreads = 1 +numOfThreads = 4 availMem = 11e9 #Grid @@ -14,7 +14,12 @@ uLB = 0.001 Re = 10 -outTime = 100 -endTime = 100 +outTime = 1000 +endTime = 1000 -logToFile = false \ No newline at end of file +logToFile = false + +cpStep = 10000 +cpStepStart= 100000 +restart = false +restartStep= 100000 \ No newline at end of file diff --git a/source/Applications/BoxBenchmark/bb.cpp b/source/Applications/BoxBenchmark/bb.cpp index c85de428a..aac67419c 100644 --- a/source/Applications/BoxBenchmark/bb.cpp +++ b/source/Applications/BoxBenchmark/bb.cpp @@ -14,22 +14,21 @@ void run(string configname) config.load(configname); string pathname = config.getString("pathname"); - int numOfThreads = config.getInt("numOfThreads"); + int numOfThreads = config.getValue<int>("numOfThreads"); vector<int> blocknx = config.getVector<int>("blocknx"); - double uLB = config.getDouble("uLB"); - double endTime = config.getDouble("endTime"); - double outTime = config.getDouble("outTime"); - double availMem = config.getDouble("availMem"); - int refineLevel = config.getInt("refineLevel"); - double Re = config.getDouble("Re"); - double dx = config.getDouble("dx"); + double uLB = config.getValue<double>("uLB"); + double endTime = config.getValue<double>("endTime"); + double outTime = config.getValue<double>("outTime"); + double availMem = config.getValue<double>("availMem"); + int refineLevel = config.getValue<int>("refineLevel"); + double Re = config.getValue<double>("Re"); + double dx = config.getValue<double>("dx"); vector<double> length = config.getVector<double>("length"); - bool logToFile = config.getBool("logToFile"); + bool logToFile = config.getValue<bool>("logToFile"); double cpStep = config.getValue<double>("cpStep"); double cpStepStart = config.getValue<double>("cpStepStart"); bool restart = config.getValue<bool>("restart"); - int chunk = config.getValue<int>("chunk"); double restartStep = config.getValue<double>("restartStep"); //UbLog::reportingLevel() = UbLog::logLevelFromString("DEBUG3"); @@ -86,7 +85,6 @@ void run(string configname) SPtr<UbScheduler> rSch2(new UbScheduler(cpStep, cpStepStart)); MPIIORestartCoProcessor rcp(grid, rSch2, pathname, comm); - rcp.setChunk(chunk); if (!restart) { @@ -113,7 +111,7 @@ void run(string configname) GenBlocksGridVisitor genBlocks(gridCube); grid->accept(genBlocks); - WriteBlocksSPtr<CoProcessor> ppblocks(new WriteBlocksCoProcessor(grid, SPtr<UbScheduler>(new UbScheduler(1)), pathname, WbWriterVtkXmlBinary::getInstance(), comm)); + SPtr<CoProcessor> ppblocks(new WriteBlocksCoProcessor(grid, SPtr<UbScheduler>(new UbScheduler(1)), pathname, WbWriterVtkXmlBinary::getInstance(), comm)); //int bbOption = 1; //0=simple Bounce Back, 1=quadr. BB //D3Q27BoundaryConditionAdapterPtr bcObst(new D3Q27NoSlipBCAdapter(bbOption)); @@ -128,18 +126,9 @@ void run(string configname) ppblocks.reset(); //set connectors - InterpolationProcessorPtr iProcessor(new IncompressibleOffsetInterpolationProcessor()); - SetConnectorsBlockVisitor setConnsVisitor(comm, true, D3Q27System::ENDDIR, nuLB, iProcessor); - //SPtr<ConnectorFactory> factory(new Block3DConnectorFactory()); - //ConnectorBlockVisitor setConnsVisitor(comm, nuLB, iProcessor, factory); - UBLOG(logINFO, "D3Q27SetConnectorsBlockVisitor:start"); - grid->accept(setConnsVisitor); - UBLOG(logINFO, "D3Q27SetConnectorsBlockVisitor:end"); - - //domain decomposition for threads - PQueuePartitioningGridVisitor pqPartVisitor(numOfThreads); - grid->accept(pqPartVisitor); - + //InterpolationProcessorPtr iProcessor(new IncompressibleOffsetInterpolationProcessor()); + //SetConnectorsBlockVisitor setConnsVisitor(comm, true, D3Q27System::ENDDIR, nuLB, iProcessor); + //grid->accept(setConnsVisitor); unsigned long long numberOfBlocks = (unsigned long long)grid->getNumberOfBlocks(); int ghostLayer = 3; @@ -167,9 +156,10 @@ void run(string configname) } SPtr<LBMKernel> kernel; - kernel = SPtr<LBMKernel>(new IncompressibleCumulantLBMKernel(blocknx[0], blocknx[1], blocknx[2], IncompressibleCumulantLBMKernel::NORMAL)); + kernel = SPtr<LBMKernel>(new IncompressibleCumulantLBMKernel()); SPtr<BCProcessor> bcProc(new BCProcessor()); kernel->setBCProcessor(bcProc); + kernel->setForcingX1(0.1); SetKernelBlockVisitor kernelVisitor(kernel, nuLB, availMem, needMem); grid->accept(kernelVisitor); @@ -186,7 +176,8 @@ void run(string configname) //grid->accept(bcVisitor); //initialization of distributions - InitDistributionsBlockVisitor initVisitor(nuLB, rhoLB, uLB); + InitDistributionsBlockVisitor initVisitor; + initVisitor.setVx1(0.5); grid->accept(initVisitor); if (myid == 0) UBLOG(logINFO, "Preprocess - end"); @@ -216,15 +207,17 @@ void run(string configname) } SPtr<UbScheduler> visSch(new UbScheduler(outTime)); - WriteMacroscopicQuantitiesCoProcessor pp(grid, visSch, pathname, WbWriterVtkXmlASCII::getInstance(), conv, comm); + SPtr<CoProcessor> mqCoProcessor(new WriteMacroscopicQuantitiesCoProcessor(grid, visSch, pathname, WbWriterVtkXmlASCII::getInstance(), conv, comm)); SPtr<UbScheduler> nupsSch(new UbScheduler(10, 30, 100)); - NUPSCounterCoProcessor npr(grid, nupsSch, numOfThreads, comm); + SPtr<CoProcessor> nupsCoProcessor(new NUPSCounterCoProcessor(grid, nupsSch, numOfThreads, comm)); - const SPtr<ConcreteCalculatorFactory> calculatorFactory = std::make_shared<ConcreteCalculatorFactory>(visSch); - CalculationManagerPtr calculation(new CalculationManager(grid, numOfThreads, endTime, calculatorFactory, CalculatorType::MPI)); + omp_set_num_threads(numOfThreads); + SPtr<Calculator> calculator(new BasicCalculator(grid, visSch, (int)endTime)); + calculator->addCoProcessor(nupsCoProcessor); + calculator->addCoProcessor(mqCoProcessor); if (myid == 0) UBLOG(logINFO, "Simulation-start"); - calculation->calculate(); + calculator->calculate(); if (myid == 0) UBLOG(logINFO, "Simulation-end"); } catch (std::exception& e) diff --git a/source/Applications/DLR-F16-Porous/CMakeLists.txt b/source/Applications/DLR-F16-Porous/CMakeLists.txt new file mode 100644 index 000000000..303857653 --- /dev/null +++ b/source/Applications/DLR-F16-Porous/CMakeLists.txt @@ -0,0 +1,25 @@ +CMAKE_MINIMUM_REQUIRED(VERSION 2.8) + +######################################################## +## C++ PROJECT ### +######################################################## +PROJECT(f16porous) + +INCLUDE(${SOURCE_ROOT}/IncludsList.cmake) + +################################################################# +### LOCAL FILES ### +################################################################# +FILE(GLOB SPECIFIC_FILES ${CMAKE_CURRENT_SOURCE_DIR}/*.h + ${CMAKE_CURRENT_SOURCE_DIR}/*.cpp + ${CMAKE_CURRENT_SOURCE_DIR}/*.hpp ) + +SET(ALL_SOURCES ${ALL_SOURCES} ${SPECIFIC_FILES}) +SOURCE_GROUP(src FILES ${SPECIFIC_FILES}) + +SET(CAB_ADDITIONAL_LINK_LIBRARIES VirtualFluids) + +################################################################# +### CREATE PROJECT ### +################################################################# +CREATE_CAB_PROJECT(f16porous BINARY) diff --git a/source/Applications/DLR-F16-Porous/f16-porous.cfg b/source/Applications/DLR-F16-Porous/f16-porous.cfg new file mode 100644 index 000000000..af6161cba --- /dev/null +++ b/source/Applications/DLR-F16-Porous/f16-porous.cfg @@ -0,0 +1,65 @@ +pathOut = d:/temp/DLR-F16-Porous +pathGeo = d:/Projects/SFB880/DLR-F16/Geometry + +fngFileNoTapeFull = F16_broad_Quad_noTape_full.stl +fngFileFull = F16_broad_Quad_full.stl + +fngFileNoTapeBody = F16_broad_Quad_noTape_body.stl +fngFileBody = F16_broad_Quad_body.stl + +fngFileTE = F16_broad_Quad_trailingEdge.stl +fngFileTEmesh = f16-pte-3_0.5_D_3.stl + + +accuracy = 0 + +reinit = false +pathReInit = /work/koskuche/DLR-F16_L7 +stepReInit = 10000 + +numOfThreads = 1 +availMem = 3.5e9 + +logToFile = false + +#x1min x1max x2min x2max x3min x3max [m] +boundingBox = -0.90 1.5 0.0 0.03 -0.66 0.66 + +blockNx = 10 10 10 + +refineLevel = 0 + +deltaXfine = 0.003 +#deltaXfine = 0.00075 #level 0 +#deltaXfine = 0.000375 #level 1 +#deltaXfine = 0.0001875 #level 2 +#deltaXfine = 0.00009375 #level 3 +#deltaXfine = 0.000046875 #level 4 +#deltaXfine = 0.0000234375 #level 5 + + +startDistance = -1.0e-3 +refineDistance = 0.6e-3 + +writeBlocks = false + +newStart = true +restartStep = 30 + +cpStep = 40 +cpStart = 40 + +outTimeStep = 40 +outTimeStart = 40 + +endTime = 40 + + +#Cp +pcpStart = 1000000 +pcpStop = 1000000 + +timeAvStart = 0 +timeAvStop = 40000 + +nupsStep = 1000 1000 10000000 diff --git a/source/Applications/DLR-F16-Porous/f16.cpp b/source/Applications/DLR-F16-Porous/f16.cpp new file mode 100644 index 000000000..37ec79b80 --- /dev/null +++ b/source/Applications/DLR-F16-Porous/f16.cpp @@ -0,0 +1,1155 @@ +#include <iostream> +#include <string> + +#include <PointerDefinitions.h> +#include "VirtualFluids.h" +#include <omp.h> +using namespace std; + + +////////////////////////////////////////////////////////////////////////// +void reinitGrid(SPtr<Grid3D> grid) +{ + int minInitLevel = grid->getCoarsestInitializedLevel(); + int maxInitLevel = grid->getFinestInitializedLevel(); + int gridRank = grid->getRank(); + + for (int level = minInitLevel; level<=maxInitLevel; level++) + { + std::vector<SPtr<Block3D>> blockVector; + grid->getBlocks(level, gridRank, blockVector); + for (SPtr<Block3D> block : blockVector) + { + if (block) + { + SPtr<ILBMKernel> kernel = block->getKernel(); + if (!kernel) + throw UbException(UB_EXARGS, "The LBM kernel isn't exist in block: "+block->toString()); + + SPtr<BCArray3D> bcArray = kernel->getBCProcessor()->getBCArray(); + SPtr<DistributionArray3D> distributions = kernel->getDataSet()->getFdistributions(); + + LBMReal f[D3Q27System::ENDF+1]; + LBMReal vx1,vx2,vx3,rho; + + size_t nx1 = distributions->getNX1(); + size_t nx2 = distributions->getNX2(); + size_t nx3 = distributions->getNX3(); + + for (int ix3=0; ix3<bcArray->getNX3(); ix3++) + for (int ix2=0; ix2<bcArray->getNX2(); ix2++) + for (int ix1=0; ix1<bcArray->getNX1(); ix1++) + { + distributions->getDistribution(f, ix1, ix2, ix3); + D3Q27System::calcCompMacroscopicValues(f,rho,vx1,vx2,vx3); + if (rho > 0.1 || rho < -0.1) + { + D3Q27System::calcCompFeq(f, 0, 0, 0, 0); + distributions->setDistribution(f, ix1, ix2, ix3); + distributions->setDistributionInv(f, ix1, ix2, ix3); + } + } + } + } + } +} +////////////////////////////////////////////////////////////////////////// +void voxelMatrixTransformation(vector<int> pmNX, double lthreshold, double uthreshold, vector<float> voxelDeltaX, string pathGeoTEvoxel, string pathOut, SPtr<Communicator> comm) +{ + int myid = comm->getProcessID(); + + vector<double> origin(3); + origin[0] = 0; + origin[1] = 0; + origin[2] = 0; + + SPtr<GbVoxelMatrix3D> voxelMatrix (new GbVoxelMatrix3D(pmNX[0], pmNX[1], pmNX[2], 0, lthreshold, uthreshold)); + voxelMatrix->readMatrixFromRawFile<unsigned short>(pathGeoTEvoxel, GbVoxelMatrix3D::BigEndian); + voxelMatrix->setVoxelMatrixDelta(voxelDeltaX[0], voxelDeltaX[1], voxelDeltaX[2]); + voxelMatrix->setVoxelMatrixMininum(origin[0], origin[1], origin[2]); + + voxelMatrix->rotate90aroundZ(); + voxelMatrix->rotate90aroundZ(); + voxelMatrix->rotate90aroundZ(); + voxelMatrix->rotate90aroundX(); + voxelMatrix->translate(0.2813, 0, 0.0042); + //voxelMatrix->translate(0.3, 0, 0.006); + + if (myid==0) voxelMatrix->writeToVTKImageDataAppended(pathOut+"/geo/fngTEvoxel"); + +} +////////////////////////////////////////////////////////////////////////// +void initPteBlock(SPtr<Grid3D> grid, SPtr<Block3D> block) +{ + int gridRank = grid->getRank(); + int blockRank = block->getRank(); + + if (blockRank == gridRank) + { + SPtr<ILBMKernel> kernel = block->getKernel(); + if (!kernel) + throw UbException(UB_EXARGS, "The LBM kernel isn't exist in block: "+block->toString()); + + SPtr<BCArray3D> bcArray = kernel->getBCProcessor()->getBCArray(); + SPtr<DistributionArray3D> distributions = kernel->getDataSet()->getFdistributions(); + + LBMReal f[D3Q27System::ENDF+1]; + + size_t nx1 = distributions->getNX1(); + size_t nx2 = distributions->getNX2(); + size_t nx3 = distributions->getNX3(); + + for (int ix3=0; ix3<bcArray->getNX3(); ix3++) + for (int ix2=0; ix2<bcArray->getNX2(); ix2++) + for (int ix1=0; ix1<bcArray->getNX1(); ix1++) + { + D3Q27System::calcCompFeq(f, 0, 0, 0, 0); + distributions->setDistribution(f, ix1, ix2, ix3); + distributions->setDistributionInv(f, ix1, ix2, ix3); + } + block->setActive(true); + } +} +////////////////////////////////////////////////////////////////////////// +void initPteFs(SPtr<Grid3D> grid, vector<SPtr<Block3D>>& vectorTE) +{ + for (SPtr<Block3D> block : vectorTE) + { + initPteBlock(grid, block); + } +} +////////////////////////////////////////////////////////////////////////// +void initPteBC(SPtr<Grid3D> grid, vector<SPtr<Block3D>>& vectorTE, string fngFileTEmesh, int parts, string pathGeoTE, string pathOut, SPtr<Communicator> comm) +{ + int myid = comm->getProcessID(); + + if (myid==0) UBLOG(logINFO, "vectorTE.size="<<vectorTE.size()); + + for (int i = 0; i < parts; i++) + { + SPtr<GbTriFaceMesh3D> fngMeshTEmesh; + if (myid==0) UBLOG(logINFO, "Read fngMeshTEmesh "<<fngFileTEmesh<<"_"<<UbSystem::toString(i)<<".stl"<<":start"); + fngMeshTEmesh = SPtr<GbTriFaceMesh3D>(GbTriFaceMesh3DCreator::getInstance()->readMeshFromSTLFile2(pathGeoTE+"/"+fngFileTEmesh+"_"+UbSystem::toString(i)+".stl", + "fngMeshTEmesh", GbTriFaceMesh3D::KDTREE_SAHPLIT, false)); + if (myid==0) UBLOG(logINFO, "Read fngMeshTEmesh:end"); + //if (myid==0) GbSystem3D::writeGeoObject(fngMeshTEmesh.get(), pathOut+"/geo/fngMeshTEmesh", WbWriterVtkXmlASCII::getInstance()); + + SPtr<BCAdapter> noSlipBCAdapter(new NoSlipBCAdapter()); + noSlipBCAdapter->setBcAlgorithm(SPtr<BCAlgorithm>(new NoSlipBCAlgorithm())); + SPtr<Interactor3D> fngIntrTEmesh = SPtr<D3Q27TriFaceMeshInteractor>(new D3Q27TriFaceMeshInteractor(fngMeshTEmesh, grid, noSlipBCAdapter, Interactor3D::SOLID)); + //fngMeshTEmesh.reset(); + + SPtr<GbCuboid3D> boundingBox(new GbCuboid3D(fngMeshTEmesh->getX1Minimum(), fngMeshTEmesh->getX2Minimum(), fngMeshTEmesh->getX3Minimum(), + fngMeshTEmesh->getX1Maximum(), fngMeshTEmesh->getX2Maximum(), fngMeshTEmesh->getX3Maximum())); + + for (SPtr<Block3D> block : vectorTE) + { + UbTupleDouble3 coords = grid->getBlockWorldCoordinates(block); + UbTupleDouble3 deltas = grid->getBlockLengths(block); + if (boundingBox->isCellInsideOrCuttingGbObject3D(val<1>(coords) + , val<2>(coords) + , val<3>(coords) + , val<1>(coords)+val<1>(deltas) + , val<2>(coords)+val<2>(deltas) + , val<3>(coords)+val<3>(deltas))) + { + block->setActive(true); + fngIntrTEmesh->setDifferencesToGbObject3D(block); + } + } + + //for (SPtr<Block3D> block : bb) + //{ + // block->setActive(true); + // fngIntrTEmesh->setDifferencesToGbObject3D(block); + // initPteBlock(grid, block); + //} + + dynamicPointerCast<D3Q27TriFaceMeshInteractor>(fngIntrTEmesh)->setQs(0); + } + + //fngIntrTE->removeSolidBlocks(); + //fngIntrTE->removeBcBlocks(); +} +////////////////////////////////////////////////////////////////////////// +void initPte(SPtr<Grid3D> grid, SPtr<Interactor3D> fngIntrTE, SPtr<Interactor3D> fngIntrTEmesh, string pathOut, SPtr<Communicator> comm) +{ + SetSolidBlockVisitor v1(fngIntrTE, BlockType::SOLID); + grid->accept(v1); + SetSolidBlockVisitor v2(fngIntrTE, BlockType::BC); + grid->accept(v2); + std::vector<SPtr<Block3D>>& sb = fngIntrTE->getSolidBlockSet(); + std::vector<SPtr<Block3D>>& bb = fngIntrTE->getBcBlocks(); + + for (SPtr<Block3D> block : sb) + { + block->setActive(true); + fngIntrTEmesh->setDifferencesToGbObject3D(block); + initPteBlock(grid, block); + } + + fngIntrTE->removeSolidBlocks(); + + for (SPtr<Block3D> block : bb) + { + block->setActive(true); + fngIntrTEmesh->setDifferencesToGbObject3D(block); + initPteBlock(grid, block);; + } + + dynamicPointerCast<D3Q27TriFaceMeshInteractor>(fngIntrTEmesh)->setQs(0); + + + fngIntrTE->removeBcBlocks(); +} +////////////////////////////////////////////////////////////////////////// +void setPointsTE(vector<SPtr<Block3D>>& vectorTE) +{ + for (SPtr<Block3D> block : vectorTE) + { + block->setActive(true); + + SPtr<ILBMKernel> kernel = block->getKernel(); + if (!kernel) + throw UbException(UB_EXARGS, "The LBM kernel isn't exist in block: "+block->toString()); + + SPtr<BCArray3D> bcArray = kernel->getBCProcessor()->getBCArray(); + SPtr<DistributionArray3D> distributions = kernel->getDataSet()->getFdistributions(); + + LBMReal f[D3Q27System::ENDF+1]; + + size_t nx1 = distributions->getNX1(); + size_t nx2 = distributions->getNX2(); + size_t nx3 = distributions->getNX3(); + + int minX1 = 0; + int minX2 = 0; + int minX3 = 0; + + int maxX1 = (int)bcArray->getNX1(); + int maxX2 = (int)bcArray->getNX2(); + int maxX3 = (int)bcArray->getNX3(); + + for (int ix3=minX3; ix3<maxX3; ix3++) + for (int ix2=minX2; ix2<maxX2; ix2++) + for (int ix1=minX1; ix1<maxX1; ix1++) + { + if (bcArray->isSolid(ix1, ix2, ix3) || bcArray->isUndefined(ix1,ix2,ix3)) + { + bool fFlag = true; + for (int k=-1; k<=1; k++) + for (int j=-1; j<=1; j++) + for (int i=-1; i<=1; i++) + { + if (ix1+i>minX1 && ix1+i<maxX1 && ix2+j>minX2 && ix2+j<maxX2 && ix3+k>minX3 && ix3+k<maxX3) + { + if (bcArray->isFluid(ix1+i, ix2+j, ix3+k)) + { + fFlag = fFlag && true; + } + else + { + fFlag = fFlag && false; + } + } + } + if (fFlag) + { + bcArray->setFluid(ix1, ix2, ix3); + D3Q27System::calcCompFeq(f, 0, 0, 0, 0); + distributions->setDistribution(f, ix1, ix2, ix3); + distributions->setDistributionInv(f, ix1, ix2, ix3); + } + } + /*else if (bcArray->isFluid(ix1, ix2, ix3)) + { + bool sFlag = false; + for (int k=-1; k<=1; k++) + for (int j=-1; j<=1; j++) + for (int i=-1; i<=1; i++) + { + if (ix1+i>minX1 && ix1+i<maxX1 && ix2+j>minX2 && ix2+j<maxX2 && ix3+k>minX3 && ix3+k<maxX3) + { + if (bcArray->isSolid(ix1+i, ix2+j, ix3+k)) + { + sFlag = sFlag && true; + } + else + { + sFlag = sFlag && false; + } + } + } + if (sFlag) + { + bcArray->setSolid(ix1, ix2, ix3); + } + }*/ + } + } +} + + +void run(string configname) +{ + try + { + ConfigurationFile config; + config.load(configname); + + string pathOut = config.getValue<string>("pathOut"); + string pathGeo = config.getValue<string>("pathGeo"); + //string pathGeoTE = config.getValue<string>("pathGeoTE"); + string fngFileNoTapeFull = config.getValue<string>("fngFileNoTapeFull"); + string fngFileFull = config.getValue<string>("fngFileFull"); + string fngFileNoTapeBody = config.getValue<string>("fngFileNoTapeBody"); + string fngFileBody = config.getValue<string>("fngFileBody"); + string fngFileTEmesh = config.getValue<string>("fngFileTEmesh"); + string fngFileTE = config.getValue<string>("fngFileTE"); + //string tapeFile = config.getValue<string>("tapeFile"); + int accuracy = config.getValue<int>("accuracy"); + int numOfThreads = config.getValue<int>("numOfThreads"); + vector<int> blockNx = config.getVector<int>("blockNx"); + vector<double> boundingBox = config.getVector<double>("boundingBox"); + double restartStep = config.getValue<double>("restartStep"); + double cpStart = config.getValue<double>("cpStart"); + double cpStep = config.getValue<double>("cpStep"); + int endTime = config.getValue<int>("endTime"); + double outTimeStep = config.getValue<double>("outTimeStep"); + double outTimeStart = config.getValue<double>("outTimeStart"); + double availMem = config.getValue<double>("availMem"); + int refineLevel = config.getValue<int>("refineLevel"); + bool logToFile = config.getValue<bool>("logToFile"); + double deltaXfine = config.getValue<double>("deltaXfine"); + double refineDistance = config.getValue<double>("refineDistance"); + double startDistance = config.getValue<double>("startDistance"); + vector<double> nupsStep = config.getVector<double>("nupsStep"); + bool newStart = config.getValue<bool>("newStart"); + bool writeBlocks = config.getValue<bool>("writeBlocks"); + string pathReInit = config.getValue<string>("pathReInit"); + int stepReInit = config.getValue<int>("stepReInit"); + bool reinit = config.getValue<bool>("reinit"); + + double pcpStart = config.getValue<double>("pcpStart"); + double pcpStop = config.getValue<double>("pcpStop"); + + double timeAvStart = config.getValue<double>("timeAvStart"); + double timeAvStop = config.getValue<double>("timeAvStop"); + + vector<int> pmNX = config.getVector<int>("pmNX"); + double lthreshold = config.getValue<double>("lthreshold"); + double uthreshold = config.getValue<double>("uthreshold"); + vector<float> voxelDeltaX = config.getVector<float>("voxelDeltaX"); + string pathGeoTEvoxel = config.getValue<string>("pathGeoTEvoxel"); + + + + SPtr<Communicator> comm = MPICommunicator::getInstance(); + int myid = comm->getProcessID(); + + if (logToFile) + { +#if defined(__unix__) + if (myid==0) + { + const char* str = pathOut.c_str(); + mkdir(str, S_IRWXU|S_IRWXG|S_IROTH|S_IXOTH); + } +#endif + + if (myid==0) + { + stringstream logFilename; + logFilename<<pathOut+"/logfile"+UbSystem::toString(UbSystem::getTimeStamp())+".txt"; + UbLog::output_policy::setStream(logFilename.str()); + } + } + + if (myid==0) + { + UBLOG(logINFO, "PID = "<<myid<<" Point 1"); + UBLOG(logINFO, "PID = "<<myid<<" Total Physical Memory (RAM): "<<Utilities::getTotalPhysMem()); + UBLOG(logINFO, "PID = "<<myid<<" Physical Memory currently used: "<<Utilities::getPhysMemUsed()); + UBLOG(logINFO, "PID = "<<myid<<" Physical Memory currently used by current process: "<<Utilities::getPhysMemUsedByMe()/1073741824.0<<" GB"); + } + + + //the geometry is in mm + + double g_minX1 = boundingBox[0];//*1000.0; + double g_minX2 = boundingBox[2];//*1000.0; + double g_minX3 = boundingBox[4];//*1000.0; + double g_maxX1 = boundingBox[1];//*1000.0; + double g_maxX2 = boundingBox[3];//*1000.0; + double g_maxX3 = boundingBox[5];//*1000.0; + //deltaXfine *=1000.0; + + ////////////////////////////////////////////////////////////////////////// + double deltaXcoarse = deltaXfine*(double)(1<<refineLevel); + ////////////////////////////////////////////////////////////////////////// + double blockLength = (double)blockNx[0]*deltaXcoarse; + + //########################################################################## + //## physical parameters + //########################################################################## + double Re = 1e6; + + double rhoLB = 0.0; + double rhoReal = 1.2041; //(kg/m3) + //double nueReal = 153.5e-7; //m^2/s + double uReal = 55; //m/s + double lReal = 0.3;//m + //double uReal = Re*nueReal / lReal; + double nuReal = (uReal*lReal)/Re; //m^2/s + + //##Machzahl: + //#Ma = uReal/csReal + double Ma = 0.15;//Ma-Real! + double csReal = uReal / Ma; + double hLB = lReal / deltaXcoarse; + + LBMUnitConverter unitConverter(lReal, csReal, rhoReal, hLB); + + double uLB = uReal * unitConverter.getFactorVelocityWToLb(); + double nuLB = nuReal * unitConverter.getFactorViscosityWToLb(); + double lLB = lReal/deltaXcoarse; + //double nuLB = (uLB*lLB)/Re; //0.005; + //double nuLB = 0.005; + + SPtr<LBMUnitConverter> conv = SPtr<LBMUnitConverter>(new LBMUnitConverter()); + + const int baseLevel = 0; + + //////////////////////////////////////////////////////////////////////// + //Grid + ////////////////////////////////////////////////////////////////////////// + SPtr<Grid3D> grid(new Grid3D(comm)); + + //BC adapters + SPtr<BCAdapter> noSlipBCAdapter(new NoSlipBCAdapter()); + noSlipBCAdapter->setBcAlgorithm(SPtr<BCAlgorithm>(new NoSlipBCAlgorithm())); + + //SPtr<BCAdapter> slipBCAdapter(new SlipBCAdapter()); + //slipBCAdapter->setBcAlgorithm(SPtr<BCAlgorithm>(new SlipBCAlgorithm())); + + mu::Parser fct; + fct.SetExpr("U"); + fct.DefineConst("U", uLB); + SPtr<BCAdapter> velBCAdapter(new VelocityBCAdapter(true, false, false, fct, 0, BCFunction::INFCONST)); + velBCAdapter->setBcAlgorithm(SPtr<BCAlgorithm>(new VelocityWithDensityBCAlgorithm())); + + //fct.SetExpr("U"); + //fct.DefineConst("U", 0.01); + //SPtr<BCAdapter> velBCAdapterOut(new VelocityBCAdapter(true, false, false, fct, 0, BCFunction::INFCONST)); + //velBCAdapterOut->setBcAlgorithm(SPtr<BCAlgorithm>(new VelocityBCAlgorithm())); + + SPtr<BCAdapter> outflowBCAdapter(new DensityBCAdapter(rhoLB)); + outflowBCAdapter->setBcAlgorithm(SPtr<BCAlgorithm>(new NonReflectingOutflowBCAlgorithm())); + + BoundaryConditionsBlockVisitor bcVisitor; + bcVisitor.addBC(noSlipBCAdapter); + bcVisitor.addBC(velBCAdapter); + bcVisitor.addBC(outflowBCAdapter); + + SPtr<BCProcessor> bcProc; + //bcProc = SPtr<BCProcessor>(new BCProcessor()); + bcProc = SPtr<BCProcessor>(new ThinWallBCProcessor()); + + //SPtr<LBMKernel> kernel = SPtr<LBMKernel>(new CompressibleCumulantLBMKernel()); + //dynamicPointerCast<CompressibleCumulantLBMKernel>(kernel)->setRelaxationParameter(CompressibleCumulantLBMKernel::NORMAL); + + SPtr<LBMKernel> kernel = SPtr<LBMKernel>(new CompressibleCumulant4thOrderViscosityLBMKernel()); + //dynamicPointerCast<CompressibleCumulant4thOrderViscosityLBMKernel>(kernel)->setBulkViscosity(nuLB*2.0e3); + + kernel->setBCProcessor(bcProc); + + SPtr<LBMKernel> spKernel = SPtr<LBMKernel>(new CompressibleCumulantLBMKernel()); + spKernel->setBCProcessor(bcProc); + ////////////////////////////////////////////////////////////////////////// + //restart + SPtr<UbScheduler> rSch(new UbScheduler(cpStep, cpStart)); + SPtr<MPIIORestartCoProcessor> restartCoProcessor(new MPIIORestartCoProcessor(grid, rSch, pathOut, comm)); + //SPtr<MPIIOMigrationCoProcessor> restartCoProcessor(new MPIIOMigrationCoProcessor(grid, rSch, pathOut, comm)); + restartCoProcessor->setLBMKernel(kernel); + restartCoProcessor->setBCProcessor(bcProc); + + SPtr<UbScheduler> mSch(new UbScheduler(cpStep, cpStart)); + SPtr<MPIIOMigrationCoProcessor> migCoProcessor(new MPIIOMigrationCoProcessor(grid, mSch, pathOut+"/mig", comm)); + migCoProcessor->setLBMKernel(kernel); + migCoProcessor->setBCProcessor(bcProc); + ////////////////////////////////////////////////////////////////////////// + + if (myid==0) + { + UBLOG(logINFO, "PID = "<<myid<<" Point 2"); + UBLOG(logINFO, "PID = "<<myid<<" Physical Memory currently used by current process: "<<Utilities::getPhysMemUsedByMe()/1073741824.0<<" GB"); + } + + + if (newStart) + { + //////////////////////////////////////////////////////////////////////// + //define grid + ////////////////////////////////////////////////////////////////////////// + grid->setDeltaX(deltaXcoarse); + grid->setBlockNX(blockNx[0], blockNx[1], blockNx[2]); + + 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", WbWriterVtkXmlASCII::getInstance()); + GenBlocksGridVisitor genBlocks(gridCube); + grid->accept(genBlocks); + + grid->setPeriodicX1(false); + grid->setPeriodicX2(true); + grid->setPeriodicX3(false); + + if (myid==0) + { + UBLOG(logINFO, "Parameters:"); + UBLOG(logINFO, "* Re = "<<Re); + UBLOG(logINFO, "* Ma = "<<Ma); + UBLOG(logINFO, "* velocity (uReal) = "<<uReal<<" m/s"); + UBLOG(logINFO, "* viscosity (nuReal) = "<<nuReal<<" m^2/s"); + UBLOG(logINFO, "* chord length (lReal)= "<<lReal<<" m"); + UBLOG(logINFO, "* velocity LB (uLB) = "<<uLB); + UBLOG(logINFO, "* viscosity LB (nuLB) = "<<nuLB); + UBLOG(logINFO, "* chord length (l_LB) = "<<lLB<<" dx_base"); + UBLOG(logINFO, "* dx_base = "<<deltaXcoarse<<" m"); + UBLOG(logINFO, "* dx_refine = "<<deltaXfine<<" m"); + UBLOG(logINFO, "* blocknx = "<<blockNx[0]<<"x"<<blockNx[1]<<"x"<<blockNx[2]); + UBLOG(logINFO, "* refineDistance = "<<refineDistance); + UBLOG(logINFO, "* number of levels = "<<refineLevel+1); + UBLOG(logINFO, "* number of threads = "<<numOfThreads); + UBLOG(logINFO, "* number of processes = "<<comm->getNumberOfProcesses()); + UBLOG(logINFO, "* path = "<<pathOut); + UBLOG(logINFO, "Preprozess - start"); + } + + if (myid==0) + { + UBLOG(logINFO, "PID = "<<myid<<" Point 3"); + UBLOG(logINFO, "PID = "<<myid<<" Physical Memory currently used by current process: "<<Utilities::getPhysMemUsedByMe()/1073741824.0<<" GB"); + } + ////////////////////////////////////////////////////////////////////////// + + + //voxelMatrixTransformation(pmNX, lthreshold, uthreshold, voxelDeltaX, pathGeoTEvoxel, pathOut, comm); + + //return; + + + SPtr<GbTriFaceMesh3D> fngMeshTE; + if (myid==0) UBLOG(logINFO, "Read fngMeshTE:start"); + fngMeshTE = SPtr<GbTriFaceMesh3D>(GbTriFaceMesh3DCreator::getInstance()->readMeshFromSTLFile2(pathGeo+"/"+fngFileTE, "fngMeshTE", GbTriFaceMesh3D::KDTREE_SAHPLIT, false)); + if (myid==0) UBLOG(logINFO, "Read fngMeshTE:end"); + fngMeshTE->rotate(0.0, 0.5, 0.0); + fngMeshTE->translate(0.0, 0.0, 0.0012); + if (myid==0) GbSystem3D::writeGeoObject(fngMeshTE.get(), pathOut+"/geo/fngMeshTE", WbWriterVtkXmlBinary::getInstance()); + + SPtr<Interactor3D> fngIntrTE = SPtr<D3Q27TriFaceMeshInteractor>(new D3Q27TriFaceMeshInteractor(fngMeshTE, grid, noSlipBCAdapter, Interactor3D::SOLID, (Interactor3D::Accuracy)accuracy)); + + + + if (refineLevel>0 && myid==0 && writeBlocks) + { + if (myid==0) UBLOG(logINFO, "Refinement - start"); + int rank = grid->getRank(); + grid->setRank(0); + + SPtr<GbTriFaceMesh3D> fngMeshNoTapeFull; + if (myid==0) UBLOG(logINFO, "Read fngFileNoTapeFull:start"); + fngMeshNoTapeFull = SPtr<GbTriFaceMesh3D>(GbTriFaceMesh3DCreator::getInstance()->readMeshFromSTLFile2(pathGeo+"/"+fngFileNoTapeFull, "fngMeshNoTapeBody", GbTriFaceMesh3D::KDTREE_SAHPLIT, false)); + if (myid==0) UBLOG(logINFO, "Read fngFileNoTapeFull:end"); + fngMeshNoTapeFull->rotate(0.0, 0.5, 0.0); + if (myid==0) GbSystem3D::writeGeoObject(fngMeshNoTapeFull.get(), pathOut+"/geo/fngMeshNoTapeFull", WbWriterVtkXmlBinary::getInstance()); + + SPtr<Interactor3D> fngIntrNoTapeFull = SPtr<D3Q27TriFaceMeshInteractor>(new D3Q27TriFaceMeshInteractor(fngMeshNoTapeFull, grid, noSlipBCAdapter, Interactor3D::SOLID, (Interactor3D::Accuracy)accuracy)); + + int level; + + //level 1 + level = 1; + if (refineLevel - level >= 0) + { + SPtr<GbObject3D> refBoxL1(new GbCuboid3D(0.12, -0.02625, -0.03, 1.0, 0.06, 0.03)); + if (myid==0) GbSystem3D::writeGeoObject(refBoxL1.get(), pathOut+"/geo/refBoxL1", WbWriterVtkXmlASCII::getInstance()); + RefineCrossAndInsideGbObjectBlockVisitor refVisitorBoxL1(refBoxL1, level); + grid->accept(refVisitorBoxL1); + + dynamicPointerCast<D3Q27TriFaceMeshInteractor>(fngIntrNoTapeFull)->refineBlockGridToLevel(level, startDistance, 24.0*refineDistance); + } + + //level 2 + level = 2; + if (refineLevel - level >= 0) + { + SPtr<GbObject3D> refBoxL2(new GbCuboid3D(0.15, -0.03, -0.015, 0.7, 0.06, 0.015)); + if (myid==0) GbSystem3D::writeGeoObject(refBoxL2.get(), pathOut+"/geo/refBoxL2", WbWriterVtkXmlASCII::getInstance()); + RefineCrossAndInsideGbObjectBlockVisitor refVisitorBoxL2(refBoxL2, level); + grid->accept(refVisitorBoxL2); + + dynamicPointerCast<D3Q27TriFaceMeshInteractor>(fngIntrNoTapeFull)->refineBlockGridToLevel(level, startDistance, 12.0*refineDistance); + } + + //level 3 + level = 3; + if (refineLevel - level >= 0) + { + SPtr<GbObject3D> refBoxL3(new GbCuboid3D(0.15, -0.03, -0.010, 0.4, 0.06, 0.012)); + if (myid==0) GbSystem3D::writeGeoObject(refBoxL3.get(), pathOut+"/geo/refBoxL3", WbWriterVtkXmlASCII::getInstance()); + RefineCrossAndInsideGbObjectBlockVisitor refVisitorBoxL3(refBoxL3, level); + grid->accept(refVisitorBoxL3); + + dynamicPointerCast<D3Q27TriFaceMeshInteractor>(fngIntrNoTapeFull)->refineBlockGridToLevel(level, startDistance, 6.0*refineDistance); + } + + //level 4 + level = 4; + if (refineLevel - level >= 0) + { + SPtr<GbObject3D> refBoxL4(new GbCuboid3D(0.15, -0.03, -0.005, 0.32, 0.06, 0.01)); + if (myid==0) GbSystem3D::writeGeoObject(refBoxL4.get(), pathOut+"/geo/refBoxL4", WbWriterVtkXmlASCII::getInstance()); + RefineCrossAndInsideGbObjectBlockVisitor refVisitorBoxL4(refBoxL4, level); + grid->accept(refVisitorBoxL4); + + dynamicPointerCast<D3Q27TriFaceMeshInteractor>(fngIntrNoTapeFull)->refineBlockGridToLevel(level, startDistance, 3.0*refineDistance); + } + + //level 5 + level = 5; + if (refineLevel - level >= 0) + { + dynamicPointerCast<D3Q27TriFaceMeshInteractor>(fngIntrNoTapeFull)->refineBlockGridToLevel(level, startDistance, refineDistance); + + + //SPtr<GbTriFaceMesh3D> fngMeshTE; + //if (myid==0) UBLOG(logINFO, "Read fngMeshTE:start"); + //fngMeshTE = SPtr<GbTriFaceMesh3D>(GbTriFaceMesh3DCreator::getInstance()->readMeshFromSTLFile2(pathGeo+"/"+fngFileTE, "fngMeshTE", GbTriFaceMesh3D::KDTREE_SAHPLIT, false)); + //if (myid==0) UBLOG(logINFO, "Read fngMeshTE:end"); + //fngMeshTE->rotate(0.0, 0.5, 0.0); + //fngMeshTE->translate(0.0, 0.0, 0.0012); + //if (myid==0) GbSystem3D::writeGeoObject(fngMeshTE.get(), pathOut+"/geo/fngMeshTE", WbWriterVtkXmlBinary::getInstance()); + + ////SPtr<Interactor3D> fngIntrTE = SPtr<D3Q27TriFaceMeshInteractor>(new D3Q27TriFaceMeshInteractor(fngMeshTE, grid, noSlipBCAdapter, Interactor3D::SOLID, (Interactor3D::Accuracy)accuracy)); + ////dynamicPointerCast<D3Q27TriFaceMeshInteractor>(fngIntrTE)->refineBlockGridToLevel(level, startDistance, refineDistance); + RefineCrossAndInsideGbObjectBlockVisitor refVisitorTE(fngMeshTE, level); + grid->accept(refVisitorTE); + } + + ///////delete solid blocks + if (myid==0) UBLOG(logINFO, "deleteSolidBlocks - start"); + + SPtr<GbTriFaceMesh3D> fngMeshNoTapeBody; + if (myid==0) UBLOG(logINFO, "Read fngFileNoTapeBody:start"); + fngMeshNoTapeBody = SPtr<GbTriFaceMesh3D>(GbTriFaceMesh3DCreator::getInstance()->readMeshFromSTLFile2(pathGeo+"/"+fngFileNoTapeBody, "fngMeshNoTapeBody", GbTriFaceMesh3D::KDTREE_SAHPLIT, false)); + if (myid==0) UBLOG(logINFO, "Read fngFileNoTapeBody:end"); + fngMeshNoTapeBody->rotate(0.0, 0.5, 0.0); + fngMeshNoTapeBody->translate(0.0, 0.0, -0.00011); + if (myid==0) GbSystem3D::writeGeoObject(fngMeshNoTapeBody.get(), pathOut+"/geo/fngMeshNoTapeBody", WbWriterVtkXmlBinary::getInstance()); + + SPtr<Interactor3D> fngIntrNoTapeBody = SPtr<D3Q27TriFaceMeshInteractor>(new D3Q27TriFaceMeshInteractor(fngMeshNoTapeBody, grid, noSlipBCAdapter, Interactor3D::SOLID, (Interactor3D::Accuracy)accuracy));//, Interactor3D::POINTS)); + + SetSolidBlockVisitor v(fngIntrNoTapeBody, BlockType::SOLID); + grid->accept(v); + std::vector<SPtr<Block3D>>& sb = fngIntrNoTapeBody->getSolidBlockSet(); + for (SPtr<Block3D> block : sb) + { + grid->deleteBlock(block); + } + fngIntrNoTapeBody->removeSolidBlocks(); + fngIntrNoTapeBody->removeBcBlocks(); + + + if (myid==0) UBLOG(logINFO, "deleteSolidBlocks - end"); + ////////////////////////////////////////// + + { + WriteBlocksCoProcessor ppblocks(grid, SPtr<UbScheduler>(new UbScheduler(1)), pathOut, WbWriterVtkXmlBinary::getInstance(), comm); + ppblocks.process(0); + } + + grid->setRank(rank); + + RatioBlockVisitor ratioVisitor(refineLevel); + CheckRatioBlockVisitor checkRatio(refineLevel); + int count = 0; + + do { + if (myid==0) UBLOG(logINFO, "ratioVisitor - start"); + grid->accept(ratioVisitor); + if (myid==0) UBLOG(logINFO, "ratioVisitor - end"); + if (myid==0) UBLOG(logINFO, "checkRatio - start"); + checkRatio.resetState(); + grid->accept(checkRatio); + if (myid==0) UBLOG(logINFO, "checkRatio - end"); + if (myid==0) UBLOG(logINFO, "count = "<<count++<<" state = "<<checkRatio.getState()); + } while (!checkRatio.getState()); + + + { + WriteBlocksCoProcessor ppblocks(grid, SPtr<UbScheduler>(new UbScheduler(1)), pathOut, WbWriterVtkXmlBinary::getInstance(), comm); + ppblocks.process(1); + } + + OverlapBlockVisitor overlapVisitor(refineLevel, false); + grid->accept(overlapVisitor); + + if (myid==0) UBLOG(logINFO, "Refinement - end"); + } + else if (refineLevel>0 && !writeBlocks) + { + migCoProcessor->readBlocks(0); + } + grid->updateDistributedBlocks(comm); + + std::vector<int> dirs; + for (int i = D3Q27System::E; i<=D3Q27System::TS; i++) + { + dirs.push_back(i); + } + SetInterpolationDirsBlockVisitor interDirsVisitor(dirs); + grid->accept(interDirsVisitor); + + //walls + GbCuboid3DPtr addWallZmin(new GbCuboid3D(g_minX1-blockLength, g_minX2-blockLength, g_minX3-blockLength, g_maxX1+blockLength, g_maxX2+blockLength, g_minX3)); + if (myid==0) GbSystem3D::writeGeoObject(addWallZmin.get(), pathOut+"/geo/addWallZmin", WbWriterVtkXmlASCII::getInstance()); + + GbCuboid3DPtr addWallZmax(new GbCuboid3D(g_minX1-blockLength, g_minX2-blockLength, g_maxX3, g_maxX1+blockLength, g_maxX2+blockLength, g_maxX3+blockLength)); + if (myid==0) GbSystem3D::writeGeoObject(addWallZmax.get(), pathOut+"/geo/addWallZmax", WbWriterVtkXmlASCII::getInstance()); + + //wall interactors + SPtr<D3Q27Interactor> addWallZminInt(new D3Q27Interactor(addWallZmin, grid, velBCAdapter, Interactor3D::SOLID)); + SPtr<D3Q27Interactor> addWallZmaxInt(new D3Q27Interactor(addWallZmax, grid, velBCAdapter, Interactor3D::SOLID)); + + //inflow + GbCuboid3DPtr geoInflow(new GbCuboid3D(g_minX1-blockLength, g_minX2-blockLength, g_minX3-blockLength, g_minX1, g_maxX2+blockLength, g_maxX3+blockLength)); + if (myid==0) GbSystem3D::writeGeoObject(geoInflow.get(), pathOut+"/geo/geoInflow", 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()); + + //inflow + SPtr<D3Q27Interactor> inflowIntr = SPtr<D3Q27Interactor>(new D3Q27Interactor(geoInflow, grid, velBCAdapter, Interactor3D::SOLID)); + + //outflow + SPtr<D3Q27Interactor> outflowIntr = SPtr<D3Q27Interactor>(new D3Q27Interactor(geoOutflow, grid, outflowBCAdapter, Interactor3D::SOLID)); + + //airfoil + SPtr<GbTriFaceMesh3D> fngMeshBody; + if (myid==0) UBLOG(logINFO, "Read fngFileBody:start"); + fngMeshBody = SPtr<GbTriFaceMesh3D>(GbTriFaceMesh3DCreator::getInstance()->readMeshFromSTLFile2(pathGeo+"/"+fngFileBody, "fngMeshBody", GbTriFaceMesh3D::KDTREE_SAHPLIT, false)); + if (myid==0) UBLOG(logINFO, "Read fngFileBody:end"); + fngMeshBody->rotate(0.0, 0.5, 0.0); + fngMeshBody->translate(0.0, 0.0, -0.00011); + if (myid==0) GbSystem3D::writeGeoObject(fngMeshBody.get(), pathOut+"/geo/fngMeshBody", WbWriterVtkXmlBinary::getInstance()); + + SPtr<Interactor3D> fngIntrBody = SPtr<D3Q27TriFaceMeshInteractor>(new D3Q27TriFaceMeshInteractor(fngMeshBody, grid, noSlipBCAdapter, Interactor3D::SOLID, (Interactor3D::Accuracy)accuracy)); + fngMeshBody.reset(); + + //SPtr<GbTriFaceMesh3D> fngMeshTEmesh; + //if (myid==0) UBLOG(logINFO, "Read fngMeshTEmesh:start"); + //fngMeshTEmesh = SPtr<GbTriFaceMesh3D>(GbTriFaceMesh3DCreator::getInstance()->readMeshFromSTLFile2(pathGeo+"/"+fngFileTEmesh, "fngMeshTEmesh", GbTriFaceMesh3D::KDTREE_SAHPLIT, false)); + //if (myid==0) UBLOG(logINFO, "Read fngMeshTEmesh:end"); + ////fngMeshTE->rotate(0.0, 0.5, 0.0); + ////fngMeshTE->translate(0.0, 0.0, -0.0001); + //if (myid==0) GbSystem3D::writeGeoObject(fngMeshTEmesh.get(), pathOut+"/geo/fngMeshTEmesh", WbWriterVtkXmlASCII::getInstance()); + + + //SPtr<Interactor3D> fngIntrTEmesh = SPtr<D3Q27TriFaceMeshInteractor>(new D3Q27TriFaceMeshInteractor(fngMeshTEmesh, grid, noSlipBCAdapter, Interactor3D::SOLID, (Interactor3D::Accuracy)accuracy)); + //fngMeshTEmesh.reset(); + + GbCuboid3DPtr geoAddWallP(new GbCuboid3D(0.269984375, g_minX2-blockLength, 0.0016, 0.2702421875, g_maxX2+blockLength, 0.0076)); + if (myid==0) GbSystem3D::writeGeoObject(geoAddWallP.get(), pathOut+"/geo/geoAddWallP", WbWriterVtkXmlASCII::getInstance()); + SPtr<D3Q27Interactor> addWallPIntr = SPtr<D3Q27Interactor>(new D3Q27Interactor(geoAddWallP, grid, noSlipBCAdapter, Interactor3D::SOLID)); + + ////////////////////////////////////////////////////////////////////////// + vector<double> origin(3); + origin[0] = 0; + origin[1] = 0; + origin[2] = 0; + + SPtr<GbVoxelMatrix3D> voxelMatrix1(new GbVoxelMatrix3D(pmNX[0], pmNX[1], pmNX[2], 0, lthreshold, uthreshold)); + voxelMatrix1->readMatrixFromRawFile<unsigned short>(pathGeoTEvoxel, GbVoxelMatrix3D::BigEndian); + voxelMatrix1->setVoxelMatrixDelta(voxelDeltaX[0], voxelDeltaX[1], voxelDeltaX[2]); + voxelMatrix1->setVoxelMatrixMininum(origin[0], origin[1], origin[2]); + + voxelMatrix1->rotate90aroundZ(); + voxelMatrix1->rotate90aroundZ(); + voxelMatrix1->rotate90aroundZ(); + voxelMatrix1->rotate90aroundX(); + voxelMatrix1->translate(0.2813, 0, 0.0042); + double offset = ((g_maxX2-g_minX2)/2.0 - voxelMatrix1->getLengthX2())/2.0; + voxelMatrix1->setVoxelMatrixMinX2(g_minX2+offset); + + if (myid==0) voxelMatrix1->writeToVTKImageDataAppended(pathOut+"/geo/fngTEvoxel1"); + + SPtr<D3Q27Interactor> fngIntrTEvoxel1 = SPtr<D3Q27Interactor>(new D3Q27Interactor(voxelMatrix1, grid, noSlipBCAdapter, Interactor3D::SOLID)); + + SPtr<GbVoxelMatrix3D> voxelMatrix2(new GbVoxelMatrix3D(pmNX[0], pmNX[1], pmNX[2], 0, lthreshold, uthreshold)); + voxelMatrix2->readMatrixFromRawFile<unsigned short>(pathGeoTEvoxel, GbVoxelMatrix3D::BigEndian); + voxelMatrix2->setVoxelMatrixDelta(voxelDeltaX[0], voxelDeltaX[1], voxelDeltaX[2]); + voxelMatrix2->setVoxelMatrixMininum(origin[0], origin[1], origin[2]); + + voxelMatrix2->rotate90aroundZ(); + voxelMatrix2->rotate90aroundZ(); + voxelMatrix2->rotate90aroundZ(); + voxelMatrix2->rotate90aroundX(); + voxelMatrix2->translate(0.2813, 0, 0.0042); + voxelMatrix2->mirrorY(); + voxelMatrix2->setVoxelMatrixMinX2(voxelMatrix1->getX2Maximum()); + + if (myid==0) voxelMatrix2->writeToVTKImageDataAppended(pathOut+"/geo/fngTEvoxel2"); + + SPtr<D3Q27Interactor> fngIntrTEvoxel2 = SPtr<D3Q27Interactor>(new D3Q27Interactor(voxelMatrix2, grid, noSlipBCAdapter, Interactor3D::SOLID)); + ////////////////////////////////////////////////////////////////////////// + + //////////////////////////////////////////// + //METIS + SPtr<Grid3DVisitor> metisVisitor(new MetisPartitioningGridVisitor(comm, MetisPartitioningGridVisitor::LevelBased, D3Q27System::BSW, MetisPartitioner::KWAY)); + //////////////////////////////////////////// + /////delete solid blocks + if (myid==0) UBLOG(logINFO, "deleteSolidBlocks - start"); + InteractorsHelper intHelper(grid, metisVisitor); + intHelper.addInteractor(inflowIntr); + intHelper.addInteractor(outflowIntr); + intHelper.addInteractor(addWallZminInt); + intHelper.addInteractor(addWallZmaxInt); + intHelper.addInteractor(fngIntrBody); + intHelper.addInteractor(addWallPIntr); + //intHelper.addInteractor(fngIntrTEmesh); + intHelper.addInteractor(fngIntrTEvoxel1); + intHelper.addInteractor(fngIntrTEvoxel2); + intHelper.selectBlocks(); + + if (myid==0) UBLOG(logINFO, "deleteSolidBlocks - end"); + + if (myid==0) + { + UBLOG(logINFO, "PID = "<<myid<<" Point 4"); + UBLOG(logINFO, "PID = "<<myid<<" Physical Memory currently used by current process: "<<Utilities::getPhysMemUsedByMe()/1073741824.0<<" GB"); + } + ////////////////////////////////////// + + { + WriteBlocksCoProcessor ppblocks(grid, SPtr<UbScheduler>(new UbScheduler(1)), pathOut, WbWriterVtkXmlBinary::getInstance(), comm); + ppblocks.process(2); + } + + 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"); + } + + //if (writeBlocks) + //{ + // migCoProcessor->writeBlocks(0); + //} + + SetKernelBlockVisitor kernelVisitor(kernel, nuLB, availMem, needMem); + grid->accept(kernelVisitor); + + if (myid==0) UBLOG(logINFO, "SetKernelBlockVisitor:end"); + + if (myid==0) + { + UBLOG(logINFO, "PID = "<<myid<<" Point 5"); + UBLOG(logINFO, "PID = "<<myid<<" Physical Memory currently used by current process: "<<Utilities::getPhysMemUsedByMe()/1073741824.0<<" GB"); + } + + if (refineLevel>0) + { + SetUndefinedNodesBlockVisitor undefNodesVisitor; + grid->accept(undefNodesVisitor); + } + + if (myid==0) UBLOG(logINFO, "SetUndefinedNodesBlockVisitor:end"); + + //BC + intHelper.setBC(); + + if (myid==0) UBLOG(logINFO, "intHelper.setBC():end"); + if (myid==0) UBLOG(logINFO, "PID = "<<myid<<" Physical Memory currently used by current process: "<<Utilities::getPhysMemUsedByMe()/1073741824.0<<" GB"); + + if (myid==0) UBLOG(logINFO, "initPteBC:start"); + SetSolidBlockVisitor v1(fngIntrTE, BlockType::SOLID); + grid->accept(v1); + SetSolidBlockVisitor v2(fngIntrTE, BlockType::BC); + grid->accept(v2); + std::vector<SPtr<Block3D>>& vectorTE = fngIntrTE->getSolidBlockSet(); + std::vector<SPtr<Block3D>>& bb = fngIntrTE->getBcBlocks(); + vectorTE.insert(vectorTE.end(), bb.begin(), bb.end()); + //initPteBC(grid, vectorTE, fngFileTEmesh, 4, pathGeoTE, pathOut, comm); + if (myid==0) UBLOG(logINFO, "initPteBC:end"); + + if (myid==0) + { + UBLOG(logINFO, "PID = "<<myid<<" Point 6"); + UBLOG(logINFO, "PID = "<<myid<<" Physical Memory currently used by current process: "<<Utilities::getPhysMemUsedByMe()/1073741824.0<<" GB"); + } + + //initialization of distributions + //InitDistributionsBlockVisitor initVisitor1; + //initVisitor1.setVx1(fct); + //grid->accept(initVisitor1); + + ////set connectors + InterpolationProcessorPtr iProcessor(new CompressibleOffsetInterpolationProcessor()); + SetConnectorsBlockVisitor setConnsVisitor(comm, true, D3Q27System::ENDDIR, nuLB, iProcessor); + grid->accept(setConnsVisitor); + + //initialization of distributions + if (reinit) + { + InitDistributionsBlockVisitor initVisitor1; + grid->accept(initVisitor1); + SPtr<Grid3D> oldGrid(new Grid3D(comm)); + SPtr<UbScheduler> iSch(new UbScheduler()); + SPtr<MPIIORestartCoProcessor> rcp(new MPIIORestartCoProcessor(oldGrid, iSch, pathReInit, comm)); + rcp->setLBMKernel(kernel); + rcp->setBCProcessor(bcProc); + rcp->restart(stepReInit); + InitDistributionsWithInterpolationGridVisitor initVisitor2(oldGrid, iProcessor, nuLB); + grid->accept(initVisitor2); + } + else + { + InitDistributionsBlockVisitor initVisitor; + initVisitor.setVx1(10.0/*fct*/); + grid->accept(initVisitor); + } + + initPteFs(grid, vectorTE); + //initPte(grid, fngIntrTE, fngIntrTEmesh, pathOut, comm); + + //domain decomposition for threads + //PQueuePartitioningGridVisitor pqPartVisitor(numOfThreads); + //grid->accept(pqPartVisitor); + + //bcVisitor should be accept after initialization!!!! + grid->accept(bcVisitor); + if (myid==0) UBLOG(logINFO, "grid->accept(bcVisitor):end"); + + if (myid==0) + { + UBLOG(logINFO, "PID = "<<myid<<" Point 7"); + UBLOG(logINFO, "PID = "<<myid<<" Physical Memory currently used by current process: "<<Utilities::getPhysMemUsedByMe()/1073741824.0<<" GB"); + } + + //Post process + { + SPtr<UbScheduler> geoSch(new UbScheduler(1)); + WriteBoundaryConditionsCoProcessor ppgeo(grid, geoSch, pathOut, WbWriterVtkXmlBinary::getInstance(), conv, comm); + ppgeo.process(0); + } + + if (myid==0) + { + UBLOG(logINFO, "PID = "<<myid<<" Point 8"); + UBLOG(logINFO, "PID = "<<myid<<" Physical Memory currently used by current process: "<<Utilities::getPhysMemUsedByMe()/1073741824.0<<" GB"); + } + + ////sponge layer + //////////////////////////////////////////////////////////////////////////// + + GbCuboid3DPtr spongeLayerX1max(new GbCuboid3D(g_maxX1-0.4, g_minX2-blockLength, g_minX3-blockLength, g_maxX1+blockLength, g_maxX2+blockLength, g_maxX3+blockLength)); + if (myid==0) GbSystem3D::writeGeoObject(spongeLayerX1max.get(), pathOut+"/geo/spongeLayerX1max", WbWriterVtkXmlASCII::getInstance()); + SpongeLayerBlockVisitor slVisitorX1max(spongeLayerX1max, spKernel, nuLB, D3Q27System::E); + grid->accept(slVisitorX1max); + + GbCuboid3DPtr spongeLayerX1min(new GbCuboid3D(g_minX1-blockLength, g_minX2-blockLength, g_minX3-blockLength, g_minX1+0.2, g_maxX2+blockLength, g_maxX3+blockLength)); + if (myid==0) GbSystem3D::writeGeoObject(spongeLayerX1min.get(), pathOut+"/geo/spongeLayerX1min", WbWriterVtkXmlASCII::getInstance()); + SpongeLayerBlockVisitor slVisitorX1min(spongeLayerX1min, spKernel, nuLB, D3Q27System::W); + grid->accept(slVisitorX1min); + + GbCuboid3DPtr spongeLayerX3min(new GbCuboid3D(g_minX1+0.2, g_minX2-blockLength, g_minX3-blockLength, g_maxX1-0.4, g_maxX2+blockLength, g_minX3+0.2)); + if (myid==0) GbSystem3D::writeGeoObject(spongeLayerX3min.get(), pathOut+"/geo/spongeLayerX3min", WbWriterVtkXmlASCII::getInstance()); + SpongeLayerBlockVisitor slVisitorX3min(spongeLayerX3min, spKernel, nuLB, D3Q27System::B); + grid->accept(slVisitorX3min); + + GbCuboid3DPtr spongeLayerX3max(new GbCuboid3D(g_minX1+0.2, g_minX2-blockLength, g_maxX3-0.2, g_maxX1-0.4, g_maxX2+blockLength, g_maxX3+blockLength)); + if (myid==0) GbSystem3D::writeGeoObject(spongeLayerX3max.get(), pathOut+"/geo/spongeLayerX3max", WbWriterVtkXmlASCII::getInstance()); + SpongeLayerBlockVisitor slVisitorX3max(spongeLayerX3max, spKernel, nuLB, D3Q27System::T); + grid->accept(slVisitorX3max); + + ///////////////////////////////////////////////////////////////////////////// + if (myid==0) UBLOG(logINFO, "Preprozess - end"); + } + else + { + restartCoProcessor->restart((int)restartStep); + grid->setTimeStep(restartStep); + + WriteBlocksCoProcessor ppblocks(grid, SPtr<UbScheduler>(new UbScheduler(1)), pathOut, WbWriterVtkXmlBinary::getInstance(), comm); + ppblocks.process(3); + //////////////////////////////////////////////////////////////////////////// + InterpolationProcessorPtr iProcessor(new CompressibleOffsetInterpolationProcessor()); + SetConnectorsBlockVisitor setConnsVisitor(comm, true, D3Q27System::ENDDIR, nuLB, iProcessor); + grid->accept(setConnsVisitor); + + if (myid==0) UBLOG(logINFO, "setPointsTE:start"); + SPtr<GbTriFaceMesh3D> fngMeshTE; + if (myid==0) UBLOG(logINFO, "Read fngMeshTE:start"); + fngMeshTE = SPtr<GbTriFaceMesh3D>(GbTriFaceMesh3DCreator::getInstance()->readMeshFromSTLFile2(pathGeo+"/"+fngFileTE, "fngMeshTE", GbTriFaceMesh3D::KDTREE_SAHPLIT, false)); + if (myid==0) UBLOG(logINFO, "Read fngMeshTE:end"); + fngMeshTE->rotate(0.0, 0.5, 0.0); + fngMeshTE->translate(0.0, 0.0, 0.0012); + if (myid==0) GbSystem3D::writeGeoObject(fngMeshTE.get(), pathOut+"/geo/fngMeshTE", WbWriterVtkXmlBinary::getInstance()); + SPtr<Interactor3D> fngIntrTE = SPtr<D3Q27TriFaceMeshInteractor>(new D3Q27TriFaceMeshInteractor(fngMeshTE, grid, noSlipBCAdapter, Interactor3D::SOLID, (Interactor3D::Accuracy)accuracy)); + SetSolidBlockVisitor v1(fngIntrTE, BlockType::SOLID); + grid->accept(v1); + SetSolidBlockVisitor v2(fngIntrTE, BlockType::BC); + grid->accept(v2); + std::vector<SPtr<Block3D>>& vectorTE = fngIntrTE->getSolidBlockSet(); + std::vector<SPtr<Block3D>>& bb = fngIntrTE->getBcBlocks(); + vectorTE.insert(vectorTE.end(), bb.begin(), bb.end()); + setPointsTE(vectorTE); + SPtr<UbScheduler> geoSch(new UbScheduler(1)); + WriteBoundaryConditionsCoProcessor ppgeo(grid, geoSch, pathOut, WbWriterVtkXmlBinary::getInstance(), conv, comm); + ppgeo.process(1); + if (myid==0) UBLOG(logINFO, "setPointsTE:end"); + + grid->accept(bcVisitor); + + ////sponge layer + //////////////////////////////////////////////////////////////////////////// + + GbCuboid3DPtr spongeLayerX1max(new GbCuboid3D(g_maxX1-0.4, g_minX2-blockLength, g_minX3-blockLength, g_maxX1+blockLength, g_maxX2+blockLength, g_maxX3+blockLength)); + if (myid==0) GbSystem3D::writeGeoObject(spongeLayerX1max.get(), pathOut+"/geo/spongeLayerX1max", WbWriterVtkXmlASCII::getInstance()); + SpongeLayerBlockVisitor slVisitorX1max(spongeLayerX1max, spKernel, nuLB, D3Q27System::E); + grid->accept(slVisitorX1max); + + GbCuboid3DPtr spongeLayerX1min(new GbCuboid3D(g_minX1-blockLength, g_minX2-blockLength, g_minX3-blockLength, g_minX1+0.2, g_maxX2+blockLength, g_maxX3+blockLength)); + if (myid==0) GbSystem3D::writeGeoObject(spongeLayerX1min.get(), pathOut+"/geo/spongeLayerX1min", WbWriterVtkXmlASCII::getInstance()); + SpongeLayerBlockVisitor slVisitorX1min(spongeLayerX1min, spKernel, nuLB, D3Q27System::W); + grid->accept(slVisitorX1min); + + GbCuboid3DPtr spongeLayerX3min(new GbCuboid3D(g_minX1+0.2, g_minX2-blockLength, g_minX3-blockLength, g_maxX1-0.4, g_maxX2+blockLength, g_minX3+0.2)); + if (myid==0) GbSystem3D::writeGeoObject(spongeLayerX3min.get(), pathOut+"/geo/spongeLayerX3min", WbWriterVtkXmlASCII::getInstance()); + SpongeLayerBlockVisitor slVisitorX3min(spongeLayerX3min, spKernel, nuLB, D3Q27System::B); + grid->accept(slVisitorX3min); + + GbCuboid3DPtr spongeLayerX3max(new GbCuboid3D(g_minX1+0.2, g_minX2-blockLength, g_maxX3-0.2, g_maxX1-0.4, g_maxX2+blockLength, g_maxX3+blockLength)); + if (myid==0) GbSystem3D::writeGeoObject(spongeLayerX3max.get(), pathOut+"/geo/spongeLayerX3max", WbWriterVtkXmlASCII::getInstance()); + SpongeLayerBlockVisitor slVisitorX3max(spongeLayerX3max, spKernel, nuLB, D3Q27System::T); + grid->accept(slVisitorX3max); + } + + SPtr<UbScheduler> nupsSch(new UbScheduler(nupsStep[0], nupsStep[1], nupsStep[2])); + std::shared_ptr<NUPSCounterCoProcessor> nupsCoProcessor(new NUPSCounterCoProcessor(grid, nupsSch, numOfThreads, comm)); + + SPtr<UbScheduler> stepSch(new UbScheduler(outTimeStep, outTimeStart)); + + if (myid==0) + { + UBLOG(logINFO, "PID = "<<myid<<" Point 9"); + UBLOG(logINFO, "PID = "<<myid<<" Physical Memory currently used by current process: "<<Utilities::getPhysMemUsedByMe()/1073741824.0<<" GB"); + } + + SPtr<WriteMacroscopicQuantitiesCoProcessor> writeMQCoProcessor(new WriteMacroscopicQuantitiesCoProcessor(grid, stepSch, pathOut, WbWriterVtkXmlBinary::getInstance(), conv, comm)); + + SPtr<GbObject3D> bbBox(new GbCuboid3D(g_minX1-blockLength, (g_maxX2-g_minX2)/2.0, g_minX3-blockLength, g_maxX1+blockLength, (g_maxX2-g_minX2)/2.0+deltaXcoarse, g_maxX3+blockLength)); + if (myid==0) GbSystem3D::writeGeoObject(bbBox.get(), pathOut+"/geo/bbBox", WbWriterVtkXmlASCII::getInstance()); + SPtr<WriteMQFromSelectionCoProcessor> writeMQSelectCoProcessor(new WriteMQFromSelectionCoProcessor(grid, stepSch, bbBox, pathOut, WbWriterVtkXmlBinary::getInstance(), conv, comm)); + + //SPtr<UbScheduler> tavSch(new UbScheduler(1, timeAvStart, timeAvStop)); + //SPtr<TimeAveragedValuesCoProcessor> tav(new TimeAveragedValuesCoProcessor(grid, pathOut, WbWriterVtkXmlBinary::getInstance(), tavSch, comm, + // TimeAveragedValuesCoProcessor::Density | TimeAveragedValuesCoProcessor::Velocity | TimeAveragedValuesCoProcessor::Fluctuations)); + //tav->setWithGhostLayer(true); + + //SPtr<IntegrateValuesHelper> mic1(new IntegrateValuesHelper(grid, comm, 0.3-deltaXfine, 0.015, 0.0005, 0.3, 0.015+deltaXfine, 0.0005+deltaXfine)); + //if (myid==0) GbSystem3D::writeGeoObject(mic1->getBoundingBox().get(), pathOut+"/geo/mic1", WbWriterVtkXmlBinary::getInstance()); + //SPtr<UbScheduler> stepMV(new UbScheduler(1, 0, 1000000)); + //SPtr<TimeseriesCoProcessor> tsp1(new TimeseriesCoProcessor(grid, stepMV, mic1, pathOut+"/mic/mic1", comm)); + + //SPtr<IntegrateValuesHelper> mic2(new IntegrateValuesHelper(grid, comm, 0.3+deltaXfine, 0.015, 0.001685, 0.3, 0.015+deltaXfine, 0.001685+deltaXfine)); + //if (myid==0) GbSystem3D::writeGeoObject(mic2->getBoundingBox().get(), pathOut+"/geo/mic2", WbWriterVtkXmlBinary::getInstance()); + //SPtr<TimeseriesCoProcessor> tsp2(new TimeseriesCoProcessor(grid, stepMV, mic2, pathOut+"/mic/mic2", comm)); + + //SPtr<IntegrateValuesHelper> mic3(new IntegrateValuesHelper(grid, comm, 0.3-deltaXcoarse, 0.015, -0.46+4.25*deltaXcoarse, 0.3, 0.015+deltaXcoarse, -0.46+5.25*deltaXcoarse)); + //if (myid==0) GbSystem3D::writeGeoObject(mic3->getBoundingBox().get(), pathOut+"/geo/mic3", WbWriterVtkXmlBinary::getInstance()); + //SPtr<TimeseriesCoProcessor> tsp3(new TimeseriesCoProcessor(grid, stepMV, mic3, pathOut+"/mic/mic3", comm)); + + //SPtr<IntegrateValuesHelper> mic4(new IntegrateValuesHelper(grid, comm, 0.3-deltaXcoarse, 0.015, 0.46-5.25*deltaXcoarse, 0.3, 0.015+deltaXcoarse, 0.46-4.25*deltaXcoarse)); + //if (myid==0) GbSystem3D::writeGeoObject(mic4->getBoundingBox().get(), pathOut+"/geo/mic4", WbWriterVtkXmlBinary::getInstance()); + //SPtr<TimeseriesCoProcessor> tsp4(new TimeseriesCoProcessor(grid, stepMV, mic4, pathOut+"/mic/mic4", comm)); + + //omp_set_num_threads(numOfThreads); + SPtr<Calculator> calculator(new BasicCalculator(grid, stepSch, endTime)); + calculator->addCoProcessor(nupsCoProcessor); + calculator->addCoProcessor(restartCoProcessor); + calculator->addCoProcessor(writeMQSelectCoProcessor); + calculator->addCoProcessor(writeMQCoProcessor); + //calculator->addCoProcessor(migCoProcessor); + //calculator->addCoProcessor(tav); + + + if (myid==0) UBLOG(logINFO, "Simulation-start"); + calculator->calculate(); + if (myid==0) UBLOG(logINFO, "Simulation-end"); + + if (myid==0) + { + UBLOG(logINFO, "PID = "<<myid<<" Point 10"); + UBLOG(logINFO, "PID = "<<myid<<" Physical Memory currently used by current process: "<<Utilities::getPhysMemUsedByMe()/1073741824.0<<" GB"); + } + } + catch (std::exception& e) + { + cerr<<e.what()<<endl<<flush; + } + catch (std::string& s) + { + cerr<<s<<endl; + } + catch (...) + { + cerr<<"unknown exception"<<endl; + } + +} + + +int main(int argc, char* argv[]) +{ + + if (argv!=NULL) + { + if (argv[1]!=NULL) + { + run(string(argv[1])); + } + else + { + cout<<"Configuration file must be set!: "<<argv[0]<<" <config file>"<<endl<<std::flush; + } + } + + //test_run(); + + //SuperMUC + //MPI_Finalize(); + + return 0; +} + diff --git a/source/Applications/DLR-F16-Solid/f16-solid.cfg b/source/Applications/DLR-F16-Solid/f16-solid.cfg index 0e73fa193..570058aec 100644 --- a/source/Applications/DLR-F16-Solid/f16-solid.cfg +++ b/source/Applications/DLR-F16-Solid/f16-solid.cfg @@ -1,9 +1,8 @@ -pathOut = d:/temp/DLR-F16-Porous +pathOut = d:/temp/DLR-F16-Solid-tttt pathGeo = d:/Projects/SFB880/DLR-F16/Geometry fngFileWhole1 = F16_broad_Quad_noTape_full.stl -#fngFileWhole2 = F16_broad_Quad_full.stl -fngFileWhole2 = F16_broad_Quad_body.stl +fngFileWhole2 = F16_broad_Quad_full.stl accuracy = 1 @@ -11,7 +10,7 @@ reinit = false pathReInit = /work/koskuche/DLR-F16_L7 stepReInit = 10000 -numOfThreads = 1 +numOfThreads = 4 availMem = 3.5e9 logToFile = false @@ -22,34 +21,32 @@ boundingBox = -0.90 1.5 0.0 0.03 -0.66 0.66 blockNx = 10 10 10 -refineLevel = 0 +refineLevel = 5 -deltaXfine = 0.003 +#deltaXfine = 0.003 +#deltaXfine = 0.0015 #deltaXfine = 0.00075 #level 0 #deltaXfine = 0.000375 #level 1 #deltaXfine = 0.0001875 #level 2 #deltaXfine = 0.00009375 #level 3 #deltaXfine = 0.000046875 #level 4 -#deltaXfine = 0.0000234375 #level 5 - -#startDistance = -0.0007 -#refineDistance = 0.0008 +deltaXfine = 0.0000234375 #level 5 startDistance = -1.0e-3 refineDistance = 0.6e-3 -writeBlocks = false +writeBlocks = true -newStart = true -restartStep = 100 +newStart = false +restartStep = 2000 -cpStep = 10 -cpStart = 10 +cpStep = 1000 +cpStart = 3000 -outTimeStep = 10 -outTimeStart = 10 +outTimeStep = 100 +outTimeStart = 3000 -endTime = 200 +endTime = 10000 #Cp diff --git a/source/Applications/DLR-F16-Solid/f16.cpp b/source/Applications/DLR-F16-Solid/f16.cpp index 928578940..ce6c22202 100644 --- a/source/Applications/DLR-F16-Solid/f16.cpp +++ b/source/Applications/DLR-F16-Solid/f16.cpp @@ -119,6 +119,9 @@ void run(string configname) //double nuLB = (uLB*lLB)/Re; //0.005; //double nuLB = 0.005; + UBLOG(logINFO, unitConverter.toString()); + + SPtr<LBMUnitConverter> conv = SPtr<LBMUnitConverter>(new LBMUnitConverter()); const int baseLevel = 0; @@ -170,10 +173,15 @@ void run(string configname) ////////////////////////////////////////////////////////////////////////// //restart SPtr<UbScheduler> rSch(new UbScheduler(cpStep, cpStart)); - //SPtr<MPIIORestartCoProcessor> restartCoProcessor(new MPIIORestartCoProcessor(grid, rSch, pathOut, comm)); - SPtr<MPIIOMigrationCoProcessor> restartCoProcessor(new MPIIOMigrationCoProcessor(grid, rSch, pathOut, comm)); + SPtr<MPIIORestartCoProcessor> restartCoProcessor(new MPIIORestartCoProcessor(grid, rSch, pathOut, comm)); + //SPtr<MPIIOMigrationCoProcessor> restartCoProcessor(new MPIIOMigrationCoProcessor(grid, rSch, pathOut, comm)); restartCoProcessor->setLBMKernel(kernel); restartCoProcessor->setBCProcessor(bcProc); + + SPtr<UbScheduler> mSch(new UbScheduler(cpStep, cpStart)); + SPtr<MPIIOMigrationCoProcessor> migCoProcessor(new MPIIOMigrationCoProcessor(grid, mSch, pathOut+"/mig", comm)); + migCoProcessor->setLBMKernel(kernel); + migCoProcessor->setBCProcessor(bcProc); ////////////////////////////////////////////////////////////////////////// if (myid==0) @@ -222,7 +230,7 @@ void run(string configname) UBLOG(logINFO, "Preprozess - start"); } - + //SPtr<GbObject3D> fngMeshWhole(new GbCylinder3D(15.0, 0.0, 0.0, 15.0, 100.0, 0.0, 25.0)); //GbSystem3D::writeGeoObject(fngMeshWhole.get(), pathOut + "/geo/fngMeshWholeCylinder", WbWriterVtkXmlBinary::getInstance()); @@ -503,7 +511,7 @@ void run(string configname) if (writeBlocks) { - restartCoProcessor->writeBlocks(0); + migCoProcessor->writeBlocks(0); } SetKernelBlockVisitor kernelVisitor(kernel, nuLB, availMem, needMem); @@ -545,8 +553,26 @@ void run(string configname) SetConnectorsBlockVisitor setConnsVisitor(comm, true, D3Q27System::ENDDIR, nuLB, iProcessor); grid->accept(setConnsVisitor); - //initialization of distributions //if (reinit) //{ // InitDistributionsBlockVisitor initVisitor1; // grid->accept(initVisitor1); // SPtr<Grid3D> oldGrid(new Grid3D(comm)); - // SPtr<UbScheduler> iSch(new UbScheduler()); // SPtr<MPIIORestartCoProcessor> rcp(new MPIIORestartCoProcessor(oldGrid, iSch, pathReInit, comm)); // rcp->setLBMKernel(kernel); // rcp->setBCProcessor(bcProc); // rcp->restart(stepReInit); // InitDistributionsWithInterpolationGridVisitor initVisitor2(oldGrid, iProcessor, nuLB); // grid->accept(initVisitor2); //} //else //{ // InitDistributionsBlockVisitor initVisitor; // initVisitor.setVx1(fct); // grid->accept(initVisitor); //} + //initialization of distributions + if (reinit) + { + InitDistributionsBlockVisitor initVisitor1; + grid->accept(initVisitor1); + SPtr<Grid3D> oldGrid(new Grid3D(comm)); + SPtr<UbScheduler> iSch(new UbScheduler()); + SPtr<MPIIORestartCoProcessor> rcp(new MPIIORestartCoProcessor(oldGrid, iSch, pathReInit, comm)); + rcp->setLBMKernel(kernel); + rcp->setBCProcessor(bcProc); + rcp->restart(stepReInit); + InitDistributionsWithInterpolationGridVisitor initVisitor2(oldGrid, iProcessor, nuLB); + grid->accept(initVisitor2); + } + else + { + InitDistributionsBlockVisitor initVisitor; + initVisitor.setVx1(fct); + grid->accept(initVisitor); + } //domain decomposition for threads //PQueuePartitioningGridVisitor pqPartVisitor(numOfThreads); @@ -603,6 +629,12 @@ void run(string configname) } else { + + //GbCuboid3DPtr mic5(new GbCuboid3D(0.3+deltaXfine, 0.015, 0.000517+0.00037+7.0*deltaXfine, 0.3+2.0*deltaXfine, 0.015+deltaXfine, 0.000517+0.00037+8.0*deltaXfine)); + //if (myid==0) GbSystem3D::writeGeoObject(mic5.get(), pathOut+"/geo/mic5", WbWriterVtkXmlBinary::getInstance()); + + //return; + restartCoProcessor->restart((int)restartStep); grid->setTimeStep(restartStep); //////////////////////////////////////////////////////////////////////////// @@ -658,17 +690,39 @@ void run(string configname) // TimeAveragedValuesCoProcessor::Density | TimeAveragedValuesCoProcessor::Velocity | TimeAveragedValuesCoProcessor::Fluctuations)); //tav->setWithGhostLayer(true); - SPtr<IntegrateValuesHelper> mic1(new IntegrateValuesHelper(grid, comm, 0.3-deltaXcoarse, 0.035, -0.6-deltaXcoarse, 0.3, 0.065, -0.6)); + SPtr<IntegrateValuesHelper> mic1(new IntegrateValuesHelper(grid, comm, 0.3-deltaXfine, 0.015, 0.0005, 0.3, 0.015+deltaXfine, 0.0005+deltaXfine)); if (myid==0) GbSystem3D::writeGeoObject(mic1->getBoundingBox().get(), pathOut+"/geo/mic1", WbWriterVtkXmlBinary::getInstance()); - SPtr<UbScheduler> stepMV(new UbScheduler(1)); - //TimeseriesCoProcessor tsp1(grid, stepMV, mic1, pathOut+"/mic/mic1", comm); + SPtr<UbScheduler> stepMV(new UbScheduler(1, 0, 1000000)); + SPtr<TimeseriesCoProcessor> tsp1(new TimeseriesCoProcessor(grid, stepMV, mic1, pathOut+"/mic/mic1", comm)); + + SPtr<IntegrateValuesHelper> mic2(new IntegrateValuesHelper(grid, comm, 0.3+deltaXfine, 0.015, 0.001685, 0.3, 0.015+deltaXfine, 0.001685+deltaXfine)); + if (myid==0) GbSystem3D::writeGeoObject(mic2->getBoundingBox().get(), pathOut+"/geo/mic2", WbWriterVtkXmlBinary::getInstance()); + SPtr<TimeseriesCoProcessor> tsp2(new TimeseriesCoProcessor(grid, stepMV, mic2, pathOut+"/mic/mic2", comm)); + + SPtr<IntegrateValuesHelper> mic3(new IntegrateValuesHelper(grid, comm, 0.3-deltaXcoarse, 0.015, -0.46+4.25*deltaXcoarse, 0.3, 0.015+deltaXcoarse, -0.46+5.25*deltaXcoarse)); + if (myid==0) GbSystem3D::writeGeoObject(mic3->getBoundingBox().get(), pathOut+"/geo/mic3", WbWriterVtkXmlBinary::getInstance()); + SPtr<TimeseriesCoProcessor> tsp3(new TimeseriesCoProcessor(grid, stepMV, mic3, pathOut+"/mic/mic3", comm)); + + SPtr<IntegrateValuesHelper> mic4(new IntegrateValuesHelper(grid, comm, 0.3-deltaXcoarse, 0.015, 0.46-5.25*deltaXcoarse, 0.3, 0.015+deltaXcoarse, 0.46-4.25*deltaXcoarse)); + if (myid==0) GbSystem3D::writeGeoObject(mic4->getBoundingBox().get(), pathOut+"/geo/mic4", WbWriterVtkXmlBinary::getInstance()); + SPtr<TimeseriesCoProcessor> tsp4(new TimeseriesCoProcessor(grid, stepMV, mic4, pathOut+"/mic/mic4", comm)); + + SPtr<IntegrateValuesHelper> mic5(new IntegrateValuesHelper(grid, comm,0.3+deltaXfine, 0.015, 0.000517+0.00037+7.0*deltaXfine, 0.3+2.0*deltaXfine, 0.015+deltaXfine, 0.000517+0.00037+8.0*deltaXfine)); + if (myid==0) GbSystem3D::writeGeoObject(mic5->getBoundingBox().get(), pathOut+"/geo/mic5", WbWriterVtkXmlBinary::getInstance()); + SPtr<TimeseriesCoProcessor> tsp5(new TimeseriesCoProcessor(grid, stepMV, mic5, pathOut+"/mic/mic5", comm)); omp_set_num_threads(numOfThreads); - SPtr<Calculator> calculator(new BasicCalculator(grid, stepSch, endTime)); + SPtr<UbScheduler> stepGhostLayer(new UbScheduler(100)); + SPtr<Calculator> calculator(new BasicCalculator(grid, stepGhostLayer, endTime)); calculator->addCoProcessor(nupsCoProcessor); calculator->addCoProcessor(restartCoProcessor); calculator->addCoProcessor(writeMQSelectCoProcessor); calculator->addCoProcessor(writeMQCoProcessor); + calculator->addCoProcessor(tsp1); + calculator->addCoProcessor(tsp2); + calculator->addCoProcessor(tsp3); + calculator->addCoProcessor(tsp4); + calculator->addCoProcessor(tsp5); //calculator->addCoProcessor(tav); @@ -702,7 +756,7 @@ void run(string configname) // try // { // -// SPtr<Communicator> comm = MPICommunicator::getInstance(); +// SPtr<Communicator> comm = MPICommunicator::getInstance();. // int myid = comm->getProcessID(); // // if (myid==0) @@ -790,6 +844,7 @@ void run(string configname) int main(int argc, char* argv[]) { + //Sleep(30000); if (argv!=NULL) { diff --git a/source/ThirdParty/Library/numerics/geometry3d/GbTriFaceMesh3D.h b/source/ThirdParty/Library/numerics/geometry3d/GbTriFaceMesh3D.h index 80a510b16..8abadd3a0 100644 --- a/source/ThirdParty/Library/numerics/geometry3d/GbTriFaceMesh3D.h +++ b/source/ThirdParty/Library/numerics/geometry3d/GbTriFaceMesh3D.h @@ -210,12 +210,12 @@ public: ny *= length; nz *= length; } - else - { - std::cerr<<"GbTriFaceMesh3D::TriFace - calculateNormal: nx=ny=nz=0 -> kann nich sein " - <<"(dreieck hat evtl knoten doppelt oder ist ne Linie)" - <<"->removeRedunantNodes"<<std::endl; - } + //else + //{ + // std::cerr<<"GbTriFaceMesh3D::TriFace - calculateNormal: nx=ny=nz=0 -> kann nich sein " + // <<"(dreieck hat evtl knoten doppelt oder ist ne Linie)" + // <<"->removeRedunantNodes"<<std::endl; + //} } #ifdef CAB_RCF template<class Archive> diff --git a/source/VirtualFluidsCore/CoProcessors/TimeDependentBCCoProcessor.cpp b/source/VirtualFluidsCore/CoProcessors/TimeDependentBCCoProcessor.cpp index 0aea01839..9e97fee38 100644 --- a/source/VirtualFluidsCore/CoProcessors/TimeDependentBCCoProcessor.cpp +++ b/source/VirtualFluidsCore/CoProcessors/TimeDependentBCCoProcessor.cpp @@ -4,7 +4,7 @@ #include "UbScheduler.h" #include "Grid3D.h" -TimeDependentBCCoProcessor::TimeDependentBCCoProcessor(SPtr<Grid3D> grid) : CoProcessor(grid, SPtr<UbScheduler>(new UbScheduler(1))) +TimeDependentBCCoProcessor::TimeDependentBCCoProcessor(SPtr<Grid3D> grid, SPtr<UbScheduler> s) : CoProcessor(grid, s) { } @@ -16,9 +16,12 @@ TimeDependentBCCoProcessor::~TimeDependentBCCoProcessor() ////////////////////////////////////////////////////////////////////////// void TimeDependentBCCoProcessor::process(double step) { - for(SPtr<Interactor3D> inter : interactors) - inter->updateInteractor( step ); - UBLOG(logDEBUG3, "TimeDependentBCCoProcessor::update:" << step); + if(scheduler->isDue(step) ) + { + for (SPtr<Interactor3D> inter : interactors) + inter->updateInteractor(step); + UBLOG(logDEBUG3, "TimeDependentBCCoProcessor::update:" << step); + } } ////////////////////////////////////////////////////////////////////////// void TimeDependentBCCoProcessor::addInteractor( SPtr<Interactor3D> interactor ) diff --git a/source/VirtualFluidsCore/CoProcessors/TimeDependentBCCoProcessor.h b/source/VirtualFluidsCore/CoProcessors/TimeDependentBCCoProcessor.h index 4a0372207..3e3795d6a 100644 --- a/source/VirtualFluidsCore/CoProcessors/TimeDependentBCCoProcessor.h +++ b/source/VirtualFluidsCore/CoProcessors/TimeDependentBCCoProcessor.h @@ -15,7 +15,7 @@ class Grid3D; class TimeDependentBCCoProcessor : public CoProcessor { public: - TimeDependentBCCoProcessor(SPtr<Grid3D> grid); + TimeDependentBCCoProcessor(SPtr<Grid3D> grid, SPtr<UbScheduler> s); virtual ~TimeDependentBCCoProcessor(); void process(double step) override; diff --git a/source/VirtualFluidsCore/CoProcessors/WriteMQFromSelectionCoProcessor.cpp b/source/VirtualFluidsCore/CoProcessors/WriteMQFromSelectionCoProcessor.cpp index 1d3221ddf..a27392981 100644 --- a/source/VirtualFluidsCore/CoProcessors/WriteMQFromSelectionCoProcessor.cpp +++ b/source/VirtualFluidsCore/CoProcessors/WriteMQFromSelectionCoProcessor.cpp @@ -163,9 +163,9 @@ void WriteMQFromSelectionCoProcessor::addDataMQ(SPtr<Block3D> block) calcMacros = &D3Q27System::calcIncompMacroscopicValues; } - int minX1 = 0; - int minX2 = 0; - int minX3 = 0; + int minX1 = 1; + int minX2 = 1; + int minX3 = 1; int maxX1 = (int)(distributions->getNX1()); int maxX2 = (int)(distributions->getNX2()); diff --git a/source/VirtualFluidsCore/Grid/BasicCalculator.cpp b/source/VirtualFluidsCore/Grid/BasicCalculator.cpp index 984c93579..2b3505fd1 100644 --- a/source/VirtualFluidsCore/Grid/BasicCalculator.cpp +++ b/source/VirtualFluidsCore/Grid/BasicCalculator.cpp @@ -44,9 +44,9 @@ void BasicCalculator::calculate() double time[6]; #endif - for (calcStep = startTimeStep; calcStep<=numberOfTimeSteps+1; calcStep++) + for (calcStep = startTimeStep; calcStep<=numberOfTimeSteps; calcStep++) { - coProcess((double)(calcStep-1)); + //coProcess((double)(calcStep-1)); ////////////////////////////////////////////////////////////////////////// #ifdef TIMING @@ -124,10 +124,12 @@ void BasicCalculator::calculate() } } //exchange data between blocks for visualization - if ((int)additionalGhostLayerUpdateScheduler->getNextDueTime()==calcStep) + //if ((int)additionalGhostLayerUpdateScheduler->getNextDueTime()==calcStep) + if (additionalGhostLayerUpdateScheduler->isDue(calcStep)) { exchangeBlockData(straightStartLevel, maxInitLevel); } + coProcess((double)(calcStep)); //now ghost nodes have actual values } UBLOG(logDEBUG1, "OMPCalculator::calculate() - stoped"); @@ -162,7 +164,7 @@ void BasicCalculator::calculateBlocks(int startLevel, int maxInitLevel, int calc for (int i =0; i<size; i++) { blockTemp = blocks[level][i]; - blockTemp->getKernel()->calculate(); + blockTemp->getKernel()->calculate(calcStep); } //timer.stop(); //UBLOG(logINFO, "level = " << level << " blocks = " << blocks[level].size() << " collision time = " << timer.getTotalTime()); @@ -171,7 +173,7 @@ void BasicCalculator::calculateBlocks(int startLevel, int maxInitLevel, int calc catch (std::exception& e) { UBLOG(logERROR, e.what()); - UBLOG(logERROR, blockTemp->toString()<<" step = "<<calcStep); + //UBLOG(logERROR, blockTemp->toString()<<" step = "<<calcStep); //throw; exit(EXIT_FAILURE); } diff --git a/source/VirtualFluidsCore/Interactors/D3Q27TriFaceMeshInteractor.cpp b/source/VirtualFluidsCore/Interactors/D3Q27TriFaceMeshInteractor.cpp index 53bcbf9f4..3ec38f463 100644 --- a/source/VirtualFluidsCore/Interactors/D3Q27TriFaceMeshInteractor.cpp +++ b/source/VirtualFluidsCore/Interactors/D3Q27TriFaceMeshInteractor.cpp @@ -68,6 +68,8 @@ bool D3Q27TriFaceMeshInteractor::setDifferencesToGbObject3D(const SPtr<Block3D> { if(!block) return false; + //UBLOG(logINFO, "D3Q27TriFaceMeshInteractor::setDifferencesToGbObject3D()"); + bcNodeIndicesMap[block] = set< std::vector<int> >(); set< std::vector<int> >& transNodeIndices = bcNodeIndicesMap[block]; solidNodeIndicesMap[block] = set< UbTupleInt3 >(); diff --git a/source/VirtualFluidsCore/Interactors/Interactor3D.h b/source/VirtualFluidsCore/Interactors/Interactor3D.h index 9325848ef..ebbe580f3 100644 --- a/source/VirtualFluidsCore/Interactors/Interactor3D.h +++ b/source/VirtualFluidsCore/Interactors/Interactor3D.h @@ -48,6 +48,7 @@ public: virtual SPtr<GbObject3D> getGbObject3D() const { return geoObject3D; } virtual bool setDifferencesToGbObject3D(const SPtr<Block3D> block/*, const double& x1, const double& x2, const double& x3, const double& blockLengthX1, const double& blockLengthX2, const double& blockLengthX3, const double& timestep=0*/) { + //UBLOG(logINFO, "Interactor3D::setDifferencesToGbObject3D()"); return false; } diff --git a/source/VirtualFluidsCore/LBM/CompressibleCumulant4thOrderViscosityLBMKernel.cpp b/source/VirtualFluidsCore/LBM/CompressibleCumulant4thOrderViscosityLBMKernel.cpp index 9284224cd..6ba40cdad 100644 --- a/source/VirtualFluidsCore/LBM/CompressibleCumulant4thOrderViscosityLBMKernel.cpp +++ b/source/VirtualFluidsCore/LBM/CompressibleCumulant4thOrderViscosityLBMKernel.cpp @@ -5,6 +5,7 @@ #include <math.h> #include "DataSet3D.h" #include "LBMKernel.h" +#include "Block3D.h" #define PROOF_CORRECTNESS @@ -39,6 +40,7 @@ SPtr<LBMKernel> CompressibleCumulant4thOrderViscosityLBMKernel::clone() kernel->setForcingX3(muForcingX3); kernel->setIndex(ix1, ix2, ix3); kernel->setDeltaT(deltaT); + kernel->setBlock(block.lock()); if (bulkViscosity != 0) { @@ -53,7 +55,7 @@ SPtr<LBMKernel> CompressibleCumulant4thOrderViscosityLBMKernel::clone() return kernel; } ////////////////////////////////////////////////////////////////////////// -void CompressibleCumulant4thOrderViscosityLBMKernel::calculate() +void CompressibleCumulant4thOrderViscosityLBMKernel::calculate(int step) { using namespace D3Q27System; using namespace std; @@ -999,9 +1001,8 @@ void CompressibleCumulant4thOrderViscosityLBMKernel::calculate() { UB_THROW(UbException(UB_EXARGS, "rho="+UbSystem::toString(drho)+", rho_post="+UbSystem::toString(drho_post) +" dif="+UbSystem::toString(dif) - +" rho is not correct for node "+UbSystem::toString(x1)+","+UbSystem::toString(x2)+","+UbSystem::toString(x3))); - //UBLOG(logERROR,"LBMKernelETD3Q27CCLB::collideAll(): rho is not correct for node "+UbSystem::toString(x1)+","+UbSystem::toString(x2)+","+UbSystem::toString(x3)); - //exit(EXIT_FAILURE); + +" rho is not correct for node "+UbSystem::toString(x1)+","+UbSystem::toString(x2)+","+UbSystem::toString(x3) + +" in " + block.lock()->toString()+" step = "+UbSystem::toString(step))); } #endif ////////////////////////////////////////////////////////////////////////// diff --git a/source/VirtualFluidsCore/LBM/CompressibleCumulant4thOrderViscosityLBMKernel.h b/source/VirtualFluidsCore/LBM/CompressibleCumulant4thOrderViscosityLBMKernel.h index 679d17068..61d69239d 100644 --- a/source/VirtualFluidsCore/LBM/CompressibleCumulant4thOrderViscosityLBMKernel.h +++ b/source/VirtualFluidsCore/LBM/CompressibleCumulant4thOrderViscosityLBMKernel.h @@ -19,7 +19,7 @@ public: public: CompressibleCumulant4thOrderViscosityLBMKernel(); virtual ~CompressibleCumulant4thOrderViscosityLBMKernel(void); - virtual void calculate(); + virtual void calculate(int step); virtual SPtr<LBMKernel> clone(); double getCalculationTime() override; void setBulkViscosity(LBMReal value); diff --git a/source/VirtualFluidsCore/LBM/CompressibleCumulantLBMKernel.cpp b/source/VirtualFluidsCore/LBM/CompressibleCumulantLBMKernel.cpp index 65fc41cd8..dd9725b1e 100644 --- a/source/VirtualFluidsCore/LBM/CompressibleCumulantLBMKernel.cpp +++ b/source/VirtualFluidsCore/LBM/CompressibleCumulantLBMKernel.cpp @@ -4,6 +4,8 @@ #include "D3Q27EsoTwist3DSplittedVector.h" #include <math.h> #include "DataSet3D.h" +#include "Block3D.h" + #define PROOF_CORRECTNESS ////////////////////////////////////////////////////////////////////////// @@ -41,6 +43,8 @@ SPtr<LBMKernel> CompressibleCumulantLBMKernel::clone() kernel->setForcingX3(muForcingX3); kernel->setIndex(ix1, ix2, ix3); kernel->setDeltaT(deltaT); + kernel->setBlock(block.lock()); + switch (parameter) { case NORMAL: @@ -62,7 +66,7 @@ SPtr<LBMKernel> CompressibleCumulantLBMKernel::clone() return kernel; } ////////////////////////////////////////////////////////////////////////// -void CompressibleCumulantLBMKernel::calculate() +void CompressibleCumulantLBMKernel::calculate(int step) { using namespace D3Q27System; using namespace std; @@ -980,9 +984,8 @@ void CompressibleCumulantLBMKernel::calculate() { UB_THROW(UbException(UB_EXARGS, "rho="+UbSystem::toString(drho)+", rho_post="+UbSystem::toString(drho_post) +" dif="+UbSystem::toString(dif) - +" rho is not correct for node "+UbSystem::toString(x1)+","+UbSystem::toString(x2)+","+UbSystem::toString(x3))); - //UBLOG(logERROR,"LBMKernelETD3Q27CCLB::collideAll(): rho is not correct for node "+UbSystem::toString(x1)+","+UbSystem::toString(x2)+","+UbSystem::toString(x3)); - //exit(EXIT_FAILURE); + +" rho is not correct for node "+UbSystem::toString(x1)+","+UbSystem::toString(x2)+","+UbSystem::toString(x3) + +" in " + block.lock()->toString()+" step = "+UbSystem::toString(step))); } #endif ////////////////////////////////////////////////////////////////////////// diff --git a/source/VirtualFluidsCore/LBM/CompressibleCumulantLBMKernel.h b/source/VirtualFluidsCore/LBM/CompressibleCumulantLBMKernel.h index 35a52dc71..bfededbe6 100644 --- a/source/VirtualFluidsCore/LBM/CompressibleCumulantLBMKernel.h +++ b/source/VirtualFluidsCore/LBM/CompressibleCumulantLBMKernel.h @@ -19,7 +19,7 @@ public: public: CompressibleCumulantLBMKernel(); virtual ~CompressibleCumulantLBMKernel(void); - virtual void calculate(); + virtual void calculate(int step); virtual SPtr<LBMKernel> clone(); double getCalculationTime(); void setBulkOmegaToOmega(bool value); diff --git a/source/VirtualFluidsCore/LBM/ILBMKernel.h b/source/VirtualFluidsCore/LBM/ILBMKernel.h index 2786e6a87..dcd47af52 100644 --- a/source/VirtualFluidsCore/LBM/ILBMKernel.h +++ b/source/VirtualFluidsCore/LBM/ILBMKernel.h @@ -12,7 +12,7 @@ class ILBMKernel public: virtual ~ILBMKernel() {}; - virtual void calculate() = 0; + virtual void calculate(int step) = 0; virtual double getCalculationTime() = 0; virtual void swapDistributions() = 0; diff --git a/source/VirtualFluidsCore/LBM/IncompressibleCumulantLBMKernel.cpp b/source/VirtualFluidsCore/LBM/IncompressibleCumulantLBMKernel.cpp index 86c6e5379..a5e082a2c 100644 --- a/source/VirtualFluidsCore/LBM/IncompressibleCumulantLBMKernel.cpp +++ b/source/VirtualFluidsCore/LBM/IncompressibleCumulantLBMKernel.cpp @@ -4,6 +4,7 @@ #include "D3Q27EsoTwist3DSplittedVector.h" #include "DataSet3D.h" #include <math.h> +#include "Block3D.h" #define PROOF_CORRECTNESS @@ -39,6 +40,8 @@ SPtr<LBMKernel> IncompressibleCumulantLBMKernel::clone() kernel->setForcingX3(muForcingX3); kernel->setIndex(ix1, ix2, ix3); kernel->setDeltaT(deltaT); + kernel->setBlock(block.lock()); + switch (parameter) { case NORMAL: @@ -51,7 +54,7 @@ SPtr<LBMKernel> IncompressibleCumulantLBMKernel::clone() return kernel; } ////////////////////////////////////////////////////////////////////////// -void IncompressibleCumulantLBMKernel::calculate() +void IncompressibleCumulantLBMKernel::calculate(int step) { using namespace D3Q27System; using namespace std; @@ -841,9 +844,8 @@ void IncompressibleCumulantLBMKernel::calculate() { UB_THROW(UbException(UB_EXARGS, "rho="+UbSystem::toString(rho)+", rho_post="+UbSystem::toString(rho_post) +" dif="+UbSystem::toString(dif) - +" rho is not correct for node "+UbSystem::toString(x1)+","+UbSystem::toString(x2)+","+UbSystem::toString(x3))); - //UBLOG(logERROR,"LBMKernelETD3Q27CCLB::collideAll(): rho is not correct for node "+UbSystem::toString(x1)+","+UbSystem::toString(x2)+","+UbSystem::toString(x3)); - //exit(EXIT_FAILURE); + +" rho is not correct for node "+UbSystem::toString(x1)+","+UbSystem::toString(x2)+","+UbSystem::toString(x3) + +" in " + block.lock()->toString()+" step = "+UbSystem::toString(step))); } #endif ////////////////////////////////////////////////////////////////////////// diff --git a/source/VirtualFluidsCore/LBM/IncompressibleCumulantLBMKernel.h b/source/VirtualFluidsCore/LBM/IncompressibleCumulantLBMKernel.h index bab1a42ad..c4eb162cf 100644 --- a/source/VirtualFluidsCore/LBM/IncompressibleCumulantLBMKernel.h +++ b/source/VirtualFluidsCore/LBM/IncompressibleCumulantLBMKernel.h @@ -21,7 +21,7 @@ public: public: IncompressibleCumulantLBMKernel(); virtual ~IncompressibleCumulantLBMKernel(void); - virtual void calculate(); + virtual void calculate(int step); virtual SPtr<LBMKernel> clone(); double getCalculationTime(); void setRelaxationParameter(Parameter p); diff --git a/source/VirtualFluidsCore/LBM/IncompressibleCumulantWithSpongeLayerLBMKernel.cpp b/source/VirtualFluidsCore/LBM/IncompressibleCumulantWithSpongeLayerLBMKernel.cpp index df0859c65..7c8c2d94e 100644 --- a/source/VirtualFluidsCore/LBM/IncompressibleCumulantWithSpongeLayerLBMKernel.cpp +++ b/source/VirtualFluidsCore/LBM/IncompressibleCumulantWithSpongeLayerLBMKernel.cpp @@ -4,6 +4,7 @@ #include <math.h> #include "DataSet3D.h" #include "BCArray3D.h" +#include "Block3D.h" #define PROOF_CORRECTNESS @@ -127,6 +128,8 @@ SPtr<LBMKernel> IncompressibleCumulantWithSpongeLayerLBMKernel::clone() kernel->setForcingX3(muForcingX3); kernel->setIndex(ix1, ix2, ix3); kernel->setDeltaT(deltaT); + kernel->setBlock(block.lock()); + switch (parameter) { case NORMAL: @@ -143,7 +146,7 @@ SPtr<LBMKernel> IncompressibleCumulantWithSpongeLayerLBMKernel::clone() return kernel; } ////////////////////////////////////////////////////////////////////////// -void IncompressibleCumulantWithSpongeLayerLBMKernel::calculate() +void IncompressibleCumulantWithSpongeLayerLBMKernel::calculate(int step) { using namespace D3Q27System; using namespace std; @@ -960,9 +963,8 @@ void IncompressibleCumulantWithSpongeLayerLBMKernel::calculate() { UB_THROW(UbException(UB_EXARGS,"rho="+UbSystem::toString(rho)+", rho_post="+UbSystem::toString(rho_post) +" dif="+UbSystem::toString(dif) - +" rho is not correct for node "+UbSystem::toString(x1)+","+UbSystem::toString(x2)+","+UbSystem::toString(x3))); - //UBLOG(logERROR,"LBMKernelETD3Q27CCLB::collideAll(): rho is not correct for node "+UbSystem::toString(x1)+","+UbSystem::toString(x2)+","+UbSystem::toString(x3)); - //exit(EXIT_FAILURE); + +" rho is not correct for node "+UbSystem::toString(x1)+","+UbSystem::toString(x2)+","+UbSystem::toString(x3) + +" in " + block.lock()->toString()+" step = "+UbSystem::toString(step))); } #endif ////////////////////////////////////////////////////////////////////////// diff --git a/source/VirtualFluidsCore/LBM/IncompressibleCumulantWithSpongeLayerLBMKernel.h b/source/VirtualFluidsCore/LBM/IncompressibleCumulantWithSpongeLayerLBMKernel.h index 8918e4fc3..980edc26e 100644 --- a/source/VirtualFluidsCore/LBM/IncompressibleCumulantWithSpongeLayerLBMKernel.h +++ b/source/VirtualFluidsCore/LBM/IncompressibleCumulantWithSpongeLayerLBMKernel.h @@ -23,7 +23,7 @@ public: IncompressibleCumulantWithSpongeLayerLBMKernel(); virtual ~IncompressibleCumulantWithSpongeLayerLBMKernel(void); SPtr<LBMKernel> clone(); - void calculate(); + void calculate(int step); void initRelaxFactor(int vdir, double vL1, double vdx, double vSP); //! \param vdir where the sponge layer is placed //! \param vL1 length of simulation domain diff --git a/source/VirtualFluidsCore/LBM/InitDensityLBMKernel.cpp b/source/VirtualFluidsCore/LBM/InitDensityLBMKernel.cpp index 3073d75f9..53b533508 100644 --- a/source/VirtualFluidsCore/LBM/InitDensityLBMKernel.cpp +++ b/source/VirtualFluidsCore/LBM/InitDensityLBMKernel.cpp @@ -852,7 +852,7 @@ double InitDensityLBMKernel::getCalculationTime() -void InitDensityLBMKernel::calculate() +void InitDensityLBMKernel::calculate(int step) { using namespace D3Q27System; diff --git a/source/VirtualFluidsCore/LBM/InitDensityLBMKernel.h b/source/VirtualFluidsCore/LBM/InitDensityLBMKernel.h index 6af546096..6b3866f3f 100644 --- a/source/VirtualFluidsCore/LBM/InitDensityLBMKernel.h +++ b/source/VirtualFluidsCore/LBM/InitDensityLBMKernel.h @@ -12,7 +12,7 @@ class InitDensityLBMKernel : public LBMKernel public: InitDensityLBMKernel(); ~InitDensityLBMKernel(); - void calculate(); + void calculate(int step); SPtr<LBMKernel> clone(); void setVelocity(int x1, int x2, int x3, LBMReal vvx, LBMReal vvy, LBMReal vvz); double getCalculationTime(); diff --git a/source/VirtualFluidsCore/LBM/LBMKernel.h b/source/VirtualFluidsCore/LBM/LBMKernel.h index 4854a5181..f1cc6a3d3 100644 --- a/source/VirtualFluidsCore/LBM/LBMKernel.h +++ b/source/VirtualFluidsCore/LBM/LBMKernel.h @@ -21,7 +21,7 @@ public: virtual SPtr<LBMKernel> clone() = 0; - virtual void calculate() = 0; + virtual void calculate(int step) = 0; virtual double getCalculationTime() = 0; void setBCProcessor(SPtr<BCProcessor> bcp); diff --git a/source/VirtualFluidsCore/LBM/VoidLBMKernel.cpp b/source/VirtualFluidsCore/LBM/VoidLBMKernel.cpp index c9fb04cea..1d8aef648 100644 --- a/source/VirtualFluidsCore/LBM/VoidLBMKernel.cpp +++ b/source/VirtualFluidsCore/LBM/VoidLBMKernel.cpp @@ -35,7 +35,7 @@ SPtr<LBMKernel> VoidLBMKernel::clone() return kernel; } ////////////////////////////////////////////////////////////////////////// -void VoidLBMKernel::calculate() +void VoidLBMKernel::calculate(int step) { } diff --git a/source/VirtualFluidsCore/LBM/VoidLBMKernel.h b/source/VirtualFluidsCore/LBM/VoidLBMKernel.h index b67384034..68c59a23c 100644 --- a/source/VirtualFluidsCore/LBM/VoidLBMKernel.h +++ b/source/VirtualFluidsCore/LBM/VoidLBMKernel.h @@ -9,7 +9,7 @@ public: VoidLBMKernel(); ~VoidLBMKernel(); SPtr<LBMKernel> clone(); - void calculate(); + void calculate(int step); double getCalculationTime(); void initDataSet(); protected: diff --git a/source/VirtualFluidsCore/Parallel/Communicator.h b/source/VirtualFluidsCore/Parallel/Communicator.h index 3f42d68d1..fa7b98feb 100644 --- a/source/VirtualFluidsCore/Parallel/Communicator.h +++ b/source/VirtualFluidsCore/Parallel/Communicator.h @@ -46,6 +46,7 @@ public: virtual void broadcast(std::vector<double>& values) = 0; protected: Communicator(){} + Communicator( const Communicator& ){} static SPtr<Communicator> instance; private: }; diff --git a/source/VirtualFluidsCore/Visitors/SpongeLayerBlockVisitor.cpp b/source/VirtualFluidsCore/Visitors/SpongeLayerBlockVisitor.cpp index 7d6426f20..4bbea1fd8 100644 --- a/source/VirtualFluidsCore/Visitors/SpongeLayerBlockVisitor.cpp +++ b/source/VirtualFluidsCore/Visitors/SpongeLayerBlockVisitor.cpp @@ -7,6 +7,8 @@ #include "LBMKernel.h" #include "UbException.h" #include "D3Q27System.h" +#include "BCArray3D.h" +#include "BCProcessor.h" #include <numerics/geometry3d/GbCuboid3D.h> using namespace std; -- GitLab