diff --git a/source/Applications/DLR-F16-Porous/f16.cpp b/source/Applications/DLR-F16-Porous/f16.cpp index fe03eab1be54aae25e9c24f25c0c598f5bb5a9cc..ffaf9aaa94d7f76830409e8f3a6e2ecb217629c6 100644 --- a/source/Applications/DLR-F16-Porous/f16.cpp +++ b/source/Applications/DLR-F16-Porous/f16.cpp @@ -174,9 +174,9 @@ void initPteBC(SPtr<Grid3D> grid, vector<SPtr<Block3D>>& vectorTE, string fngFil ////////////////////////////////////////////////////////////////////////// void initPte(SPtr<Grid3D> grid, SPtr<Interactor3D> fngIntrTE, SPtr<Interactor3D> fngIntrTEmesh, string pathOut, SPtr<Communicator> comm) { - SetSolidBlockVisitor v1(fngIntrTE, BlockType::SOLID); + SetSolidOrBoundaryBlockVisitor v1(fngIntrTE, SetSolidOrBoundaryBlockVisitor::SOLID); grid->accept(v1); - SetSolidBlockVisitor v2(fngIntrTE, BlockType::BC); + SetSolidOrBoundaryBlockVisitor v2(fngIntrTE, SetSolidOrBoundaryBlockVisitor::BC); grid->accept(v2); std::vector<SPtr<Block3D>>& sb = fngIntrTE->getSolidBlockSet(); std::vector<SPtr<Block3D>>& bb = fngIntrTE->getBcBlocks(); @@ -649,7 +649,7 @@ void run(string configname) SPtr<Interactor3D> fngIntrNoTapeBody = SPtr<D3Q27TriFaceMeshInteractor>(new D3Q27TriFaceMeshInteractor(fngMeshNoTapeBody, grid, noSlipBCAdapter, Interactor3D::SOLID, (Interactor3D::Accuracy)accuracy));//, Interactor3D::POINTS)); - SetSolidBlockVisitor v(fngIntrNoTapeBody, BlockType::SOLID); + SetSolidOrBoundaryBlockVisitor v(fngIntrNoTapeBody, SetSolidOrBoundaryBlockVisitor::SOLID); grid->accept(v); std::vector<SPtr<Block3D>>& sb = fngIntrNoTapeBody->getSolidBlockSet(); for (SPtr<Block3D> block : sb) @@ -892,9 +892,9 @@ void run(string configname) 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); + SetSolidOrBoundaryBlockVisitor v1(fngIntrTE, SetSolidOrBoundaryBlockVisitor::SOLID); grid->accept(v1); - SetSolidBlockVisitor v2(fngIntrTE, BlockType::BC); + SetSolidOrBoundaryBlockVisitor v2(fngIntrTE, SetSolidOrBoundaryBlockVisitor::BC); grid->accept(v2); std::vector<SPtr<Block3D>>& vectorTE = fngIntrTE->getSolidBlockSet(); std::vector<SPtr<Block3D>>& bb = fngIntrTE->getBcBlocks(); @@ -997,7 +997,8 @@ void run(string configname) } else { - restartCoProcessor->restart((int)restartStep); + //restartCoProcessor->restart((int)restartStep); + migCoProcessor->restart((int)restartStep); grid->setTimeStep(restartStep); WriteBlocksCoProcessor ppblocks(grid, SPtr<UbScheduler>(new UbScheduler(1)), pathOut, WbWriterVtkXmlBinary::getInstance(), comm); @@ -1007,28 +1008,40 @@ void run(string configname) 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"); - + //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)); + //SetSolidOrBoundaryBlockVisitor v1(fngIntrTE, SetSolidOrBoundaryBlockVisitor::SOLID); + //grid->accept(v1); + //SetSolidOrBoundaryBlockVisitor v2(fngIntrTE, SetSolidOrBoundaryBlockVisitor::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"); +////////////////////////////////////////////////////////////////////////// + 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/fngMeshBody2", WbWriterVtkXmlBinary::getInstance()); + SPtr<Interactor3D> fngIntrBody = SPtr<D3Q27TriFaceMeshInteractor>(new D3Q27TriFaceMeshInteractor(fngMeshBody, grid, noSlipBCAdapter, Interactor3D::SOLID, (Interactor3D::Accuracy)accuracy)); + SetSolidOrBoundaryBlockVisitor v(fngIntrBody, SetSolidOrBoundaryBlockVisitor::BC); + grid->accept(v); + fngIntrBody->initInteractor(); +////////////////////////////////////////////////////////////////////////// grid->accept(bcVisitor); ////sponge layer @@ -1072,36 +1085,51 @@ void run(string configname) 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<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> 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> 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> 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> 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> 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> 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)); - //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> mic6(new IntegrateValuesHelper(grid, comm, 0.3, 0.015, -0.46+4.25*deltaXcoarse, 0.3+deltaXcoarse, 0.015+deltaXcoarse, -0.46+5.25*deltaXcoarse)); + if (myid==0) GbSystem3D::writeGeoObject(mic6->getBoundingBox().get(), pathOut+"/geo/mic6", WbWriterVtkXmlBinary::getInstance()); + SPtr<TimeseriesCoProcessor> tsp6(new TimeseriesCoProcessor(grid, stepMV, mic6, pathOut+"/mic/mic6", comm)); - //omp_set_num_threads(numOfThreads); - SPtr<Calculator> calculator(new BasicCalculator(grid, stepSch, endTime)); + omp_set_num_threads(numOfThreads); + SPtr<UbScheduler> stepGhostLayer(new UbScheduler(1)); + SPtr<Calculator> calculator(new BasicCalculator(grid, stepGhostLayer, endTime)); calculator->addCoProcessor(nupsCoProcessor); - calculator->addCoProcessor(restartCoProcessor); + //calculator->addCoProcessor(restartCoProcessor); + calculator->addCoProcessor(migCoProcessor); calculator->addCoProcessor(writeMQSelectCoProcessor); calculator->addCoProcessor(writeMQCoProcessor); - //calculator->addCoProcessor(migCoProcessor); - //calculator->addCoProcessor(tav); + calculator->addCoProcessor(tsp1); + calculator->addCoProcessor(tsp2); + calculator->addCoProcessor(tsp3); + calculator->addCoProcessor(tsp4); + calculator->addCoProcessor(tsp5); + calculator->addCoProcessor(tsp6); + calculator->addCoProcessor(tav); if (myid==0) UBLOG(logINFO, "Simulation-start"); diff --git a/source/Applications/DLR-F16-Solid/f16.cpp b/source/Applications/DLR-F16-Solid/f16.cpp index 06681a1183cfba3f160fc296d4990b4b5be91527..e5ae1722cad1d6e13207f64d91daa4ee102f20cc 100644 --- a/source/Applications/DLR-F16-Solid/f16.cpp +++ b/source/Applications/DLR-F16-Solid/f16.cpp @@ -361,7 +361,7 @@ void run(string configname) /////delete solid blocks if (myid==0) UBLOG(logINFO, "deleteSolidBlocks - start"); - SetSolidBlockVisitor v(fngIntrWhole1, BlockType::SOLID); + SetSolidOrBoundaryBlockVisitor v(fngIntrWhole1, SetSolidOrBoundaryBlockVisitor::SOLID); grid->accept(v); std::vector<SPtr<Block3D>>& sb = fngIntrWhole1->getSolidBlockSet(); for (SPtr<Block3D> block : sb) @@ -645,6 +645,19 @@ void run(string configname) SetConnectorsBlockVisitor setConnsVisitor(comm, true, D3Q27System::ENDDIR, nuLB, iProcessor); grid->accept(setConnsVisitor); +/////////////////////////////////////////////////////////////////////////////////////////////// + SPtr<GbTriFaceMesh3D> fngMeshWhole2; + if (myid==0) UBLOG(logINFO, "Read fngFileWhole2:start"); + fngMeshWhole2 = SPtr<GbTriFaceMesh3D>(GbTriFaceMesh3DCreator::getInstance()->readMeshFromSTLFile2(pathGeo+"/"+fngFileWhole2, "fngMeshWhole2", GbTriFaceMesh3D::KDTREE_SAHPLIT, false)); + if (myid==0) UBLOG(logINFO, "Read fngFileWhole2:end"); + fngMeshWhole2->rotate(0.0, 0.5, 0.0); + if (myid==0) GbSystem3D::writeGeoObject(fngMeshWhole2.get(), pathOut+"/geo/fngMeshWhole3", WbWriterVtkXmlBinary::getInstance()); + SPtr<Interactor3D> fngIntrWhole2 = SPtr<D3Q27TriFaceMeshInteractor>(new D3Q27TriFaceMeshInteractor(fngMeshWhole2, grid, noSlipBCAdapter, Interactor3D::SOLID, (Interactor3D::Accuracy)accuracy)); + SetSolidOrBoundaryBlockVisitor v(fngIntrWhole2, SetSolidOrBoundaryBlockVisitor::BC); + grid->accept(v); + fngIntrWhole2->initInteractor(); +/////////////////////////////////////////////////////////////////////////////////////////////// + grid->accept(bcVisitor); ////sponge layer @@ -688,10 +701,10 @@ void run(string configname) 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<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()); @@ -719,7 +732,7 @@ void run(string configname) SPtr<TimeseriesCoProcessor> tsp6(new TimeseriesCoProcessor(grid, stepMV, mic6, pathOut+"/mic/mic6", comm)); omp_set_num_threads(numOfThreads); - SPtr<UbScheduler> stepGhostLayer(new UbScheduler(100)); + SPtr<UbScheduler> stepGhostLayer(new UbScheduler(1)); SPtr<Calculator> calculator(new BasicCalculator(grid, stepGhostLayer, endTime)); calculator->addCoProcessor(nupsCoProcessor); calculator->addCoProcessor(restartCoProcessor); @@ -731,7 +744,7 @@ void run(string configname) calculator->addCoProcessor(tsp4); calculator->addCoProcessor(tsp5); calculator->addCoProcessor(tsp6); - //calculator->addCoProcessor(tav); + calculator->addCoProcessor(tav); if (myid==0) UBLOG(logINFO, "Simulation-start"); diff --git a/source/Applications/PoiseuilleFlow/pf1.cpp b/source/Applications/PoiseuilleFlow/pf1.cpp index a0a15394d7534002732c974c548221c11cd865a8..85582c56f2025de146a7858fe03397aaa04bc3b9 100644 --- a/source/Applications/PoiseuilleFlow/pf1.cpp +++ b/source/Applications/PoiseuilleFlow/pf1.cpp @@ -54,8 +54,6 @@ void pf1() grid->setPeriodicX2(false); grid->setPeriodicX3(false); - if (myid == 0) GbSystem3D::writeGeoObject(gridCube.get(), pathname + "/geo/gridCube", WbWriterVtkXmlBinary::getInstance()); - //blocks generating GenBlocksGridVisitor genBlocks(gridCube); grid->accept(genBlocks); @@ -79,7 +77,7 @@ void pf1() intHelper.selectBlocks(); //write data for visualization of block grid - 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)); ppblocks->process(0); ppblocks.reset(); @@ -110,7 +108,7 @@ void pf1() //LBM kernel definition 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); @@ -130,7 +128,7 @@ void pf1() grid->accept(bcVisitor); //initialization of distributions - InitDistributionsBlockVisitor initVisitor(nuLB, rhoLB); + InitDistributionsBlockVisitor initVisitor; grid->accept(initVisitor); //set connectors @@ -138,32 +136,32 @@ void pf1() SetConnectorsBlockVisitor setConnsVisitor(comm, true, D3Q27System::ENDDIR, nuLB, iProcessor); grid->accept(setConnsVisitor); - //domain decomposition for threads - PQueuePartitioningGridVisitor pqPartVisitor(numOfThreads); - grid->accept(pqPartVisitor); - //write data for visualization of boundary conditions - SPtr<UbScheduler> geoSch(new UbScheduler(1)); - WriteBoundaryConditionsSPtr<CoProcessor> ppgeo( - new WriteBoundaryConditionsCoProcessor(grid, geoSch, pathname, WbWriterVtkXmlBinary::getInstance(), SPtr<LBMUnitConverter>(new LBMUnitConverter()), comm)); - ppgeo->process(0); - ppgeo.reset(); + { + SPtr<UbScheduler> geoSch(new UbScheduler(1)); + WriteBoundaryConditionsCoProcessor ppgeo(grid, geoSch, pathname, WbWriterVtkXmlBinary::getInstance(), SPtr<LBMUnitConverter>(new LBMUnitConverter()), comm); + ppgeo.process(0); + } if (myid == 0) UBLOG(logINFO, "Preprocess - end"); //write data for visualization of macroscopic quantities SPtr<UbScheduler> visSch(new UbScheduler(outTime)); - WriteMacroscopicQuantitiesCoProcessor pp(grid, visSch, pathname, WbWriterVtkXmlASCII::getInstance(), SPtr<LBMUnitConverter>(new LBMUnitConverter()), comm); + SPtr<WriteMacroscopicQuantitiesCoProcessor> writeMQCoProcessor(new WriteMacroscopicQuantitiesCoProcessor(grid, visSch, pathname, + WbWriterVtkXmlASCII::getInstance(), SPtr<LBMUnitConverter>(new LBMUnitConverter()), comm)); //performance control SPtr<UbScheduler> nupsSch(new UbScheduler(10, 30, 100)); - NUPSCounterCoProcessor npr(grid, nupsSch, numOfThreads, comm); + SPtr<NUPSCounterCoProcessor> npr(new NUPSCounterCoProcessor(grid, nupsSch, numOfThreads, comm)); //start simulation - const SPtr<ConcreteCalculatorFactory> calculatorFactory = std::make_shared<ConcreteCalculatorFactory>(visSch); - CalculationManagerPtr calculation(new CalculationManager(grid, numOfThreads, endTime, calculatorFactory, CalculatorType::HYBRID)); + omp_set_num_threads(numOfThreads); + SPtr<UbScheduler> stepGhostLayer(new UbScheduler(outTime)); + SPtr<Calculator> calculator(new BasicCalculator(grid, stepGhostLayer, endTime)); + calculator->addCoProcessor(npr); + calculator->addCoProcessor(writeMQCoProcessor); if (myid == 0) UBLOG(logINFO, "Simulation-start"); - calculation->calculate(); + calculator->calculate(); if (myid == 0) UBLOG(logINFO, "Simulation-end"); } diff --git a/source/Applications/pChannel/configBombadilpChannel.cfg b/source/Applications/pChannel/configBombadilpChannel.cfg index 2482d5526566727c7016faed898a9d172c63a4fe..e6625606edf6530823592818fa9e1a0bad053bfb 100644 --- a/source/Applications/pChannel/configBombadilpChannel.cfg +++ b/source/Applications/pChannel/configBombadilpChannel.cfg @@ -2,15 +2,16 @@ #Simulation parameters for porous channel # -pathname = d:/temp/ChannelFlowCorr2 -pathGeo = d:/Projects/SFB880/GeometrienPoroeseMedien/PA80-110 +pathOut = d:/temp/ChannelFlow-test +pathGeo = d:/Projects/SFB880/GeometrienPoroeseMedien/isotrop/PA80-110 numOfThreads = 1 -availMem = 8e9 +availMem = 14e9 logToFile = false #porous media #rawFile = false sampleFilename = PA80-110_275x267x254_1mm.vti +#sampleFilename = PA80-110_1096x1327x303.vti #writeSample = true rawFile = false @@ -50,19 +51,17 @@ voxelDeltaX = 3.6496350365e-6 3.76789751319e-6 3.95256916996e-6 #pmL = 4e-3 0.5e-3 1e-3 -#pmL = 1e-3 1e-3 1e-3 -pmL = 4e-3 8e-3 1e-3 +pmL = 1e-3 1e-3 1e-3 +#pmL = 4e-3 8e-3 1e-3 #grid refineLevel = 0 #deltaXfine = 10e-6 -#deltaXfine = 20e-6 -deltaXfine = 40e-5 +deltaXfine = 128e-6 +#deltaXfine = 32e-6 #level 2 + +blocknx = 10 10 10 -#blocknx = 10 10 10 -blocknx = 20 20 20 -#lengthFactor = 6 -lengthFactor = 16 thinWall = false changeQs = false @@ -70,29 +69,33 @@ changeQs = false #DLR-F15 #channelHigh = 0.017 -channelHigh = 0.016 +channelHigh = 0.015 #NACA 0012 #channelHigh = 0.008 - + #channelHigh = 0.005 + +boundingBox = 0 0 0 0.024 0.016 0.016 #physic # for DLRF-15 Re = 102000/2 -Re = 51000 +Re = 15000 #real velocity is 54.95 m/s u_LB = 0.1 -nupsSteps=1000 +newStart = false +restartStep = 10 -restartStep = 20000 -restartStepStart = 20000 +cpStep = 10 +cpStart = 10 averaging = false averagingReset = false timeAvStart = 21000000 timeAvStop = 2100010000 -endTime = 5 -outTime = 1 +endTime = 100 +outTime = 10000 +nupsStep = 100 100 10000000 diff --git a/source/Applications/pChannel/configHLRNpChannel.cfg b/source/Applications/pChannel/configHLRNpChannel.cfg index 1f5ee170c92212eb947dab3c675a9c6909bf6370..03bcfb431c60b491c7ccf2c8da2dfd1536c3c4da 100644 --- a/source/Applications/pChannel/configHLRNpChannel.cfg +++ b/source/Applications/pChannel/configHLRNpChannel.cfg @@ -51,4 +51,4 @@ restartStepStart=10000 endTime = 80000 outTime = 10000 - +nupsStep = 1000 1000 10000000 diff --git a/source/Applications/pChannel/pChannel.cpp b/source/Applications/pChannel/pChannel.cpp index a557d9070c0329f7965516cbc4baa91e5ccb8a38..6d0536fd1ad2b1ecd3b7cce2c0662e580d89b1cd 100644 --- a/source/Applications/pChannel/pChannel.cpp +++ b/source/Applications/pChannel/pChannel.cpp @@ -31,37 +31,39 @@ void run(string configname) ConfigurationFile config; config.load(configname); - string pathname = config.getString("pathname"); - string pathGeo = config.getString("pathGeo"); - int numOfThreads = config.getInt("numOfThreads"); - string sampleFilename = config.getString("sampleFilename"); + string pathOut = config.getValue<string>("pathOut"); + string pathGeo = config.getValue<string>("pathGeo"); + int numOfThreads = config.getValue<int>("numOfThreads"); + string sampleFilename = config.getValue<string>("sampleFilename"); vector<int> pmNX = config.getVector<int>("pmNX"); - double lthreshold = config.getDouble("lthreshold"); - double uthreshold = config.getDouble("uthreshold"); + double lthreshold = config.getValue<double>("lthreshold"); + double uthreshold = config.getValue<double>("uthreshold"); vector<float> voxelDeltaX = config.getVector<float>("voxelDeltaX"); vector<int> blocknx = config.getVector<int>("blocknx"); - double u_LB = config.getDouble("u_LB"); - double restartStep = config.getDouble("restartStep"); - double restartStepStart = config.getDouble("restartStepStart"); - double endTime = config.getDouble("endTime"); - double outTime = config.getDouble("outTime"); - double availMem = config.getDouble("availMem"); - bool rawFile = config.getBool("rawFile"); - bool logToFile = config.getBool("logToFile"); - bool writeSample = config.getBool("writeSample"); + double u_LB = config.getValue<double>("u_LB"); + double restartStep = config.getValue<double>("restartStep"); + double cpStep = config.getValue<double>("cpStep"); + double cpStart = config.getValue<double>("cpStart"); + double endTime = config.getValue<double>("endTime"); + double outTime = config.getValue<double>("outTime"); + double availMem = config.getValue<double>("availMem"); + bool rawFile = config.getValue<bool>("rawFile"); + bool logToFile = config.getValue<bool>("logToFile"); + bool writeSample = config.getValue<bool>("writeSample"); vector<double> pmL = config.getVector<double>("pmL"); - double deltaXfine = config.getDouble("deltaXfine"); - int refineLevel = config.getInt("refineLevel"); - bool thinWall = config.getBool("thinWall"); - double Re = config.getDouble("Re"); - double channelHigh = config.getDouble("channelHigh"); - double lengthFactor = config.getDouble("lengthFactor"); - bool changeQs = config.getBool("changeQs"); - double timeAvStart = config.getDouble("timeAvStart"); - double timeAvStop = config.getDouble("timeAvStop"); - bool averaging = config.getBool("averaging"); - bool averagingReset = config.getBool("averagingReset"); - double nupsSteps = config.getDouble("nupsSteps"); + double deltaXfine = config.getValue<double>("deltaXfine"); + int refineLevel = config.getValue<int>("refineLevel"); + bool thinWall = config.getValue<bool>("thinWall"); + double Re = config.getValue<double>("Re"); + double channelHigh = config.getValue<double>("channelHigh"); + bool changeQs = config.getValue<bool>("changeQs"); + double timeAvStart = config.getValue<double>("timeAvStart"); + double timeAvStop = config.getValue<double>("timeAvStop"); + bool averaging = config.getValue<bool>("averaging"); + bool averagingReset = config.getValue<bool>("averagingReset"); + bool newStart = config.getValue<bool>("newStart"); + vector<double> nupsStep = config.getVector<double>("nupsStep"); + vector<double> boundingBox = config.getVector<double>("boundingBox"); SPtr<Communicator> comm = MPICommunicator::getInstance(); int myid = comm->getProcessID(); @@ -71,7 +73,7 @@ void run(string configname) #if defined(__unix__) if (myid == 0) { - const char* str = pathname.c_str(); + const char* str = pathOut.c_str(); mkdir(str, S_IRWXU | S_IRWXG | S_IROTH | S_IXOTH); } #endif @@ -79,12 +81,12 @@ void run(string configname) if (myid == 0) { stringstream logFilename; - logFilename << pathname + "/logfile" + UbSystem::toString(UbSystem::getTimeStamp()) + ".txt"; + logFilename << pathOut + "/logfile" + UbSystem::toString(UbSystem::getTimeStamp()) + ".txt"; UbLog::output_policy::setStream(logFilename.str()); } } - Sleep(30000); + //Sleep(30000); if (myid == 0) UBLOG(logINFO, "Testcase porous channel"); @@ -95,8 +97,7 @@ void run(string configname) const int baseLevel = 0; double deltaXcoarse = deltaXfine*(double)(1<<refineLevel); - double coord[6]; - bool restart; + //double coord[6]; vector<double> origin(3); origin[0] = 0; @@ -104,43 +105,81 @@ void run(string configname) origin[2] = 0; //real velocity is 49.63 m/s - double nu_LB; SPtr<Grid3D> grid(new Grid3D(comm)); + //BC adapters + SPtr<BCAdapter> noSlipBCAdapter(new NoSlipBCAdapter()); + noSlipBCAdapter->setBcAlgorithm(SPtr<BCAlgorithm>(new NoSlipBCAlgorithm())); + + BoundaryConditionsBlockVisitor bcVisitor; + bcVisitor.addBC(noSlipBCAdapter); + + SPtr<BCProcessor> bcProc; + bcProc = SPtr<BCProcessor>(new BCProcessor()); + + SPtr<LBMKernel> kernel = SPtr<LBMKernel>(new IncompressibleCumulantLBMKernel()); + + mu::Parser fctForcingX1; + fctForcingX1.SetExpr("Fx1"); + fctForcingX1.DefineConst("Fx1", 1.0e-6); + kernel->setWithForcing(true); + + kernel->setBCProcessor(bcProc); + ////////////////////////////////////////////////////////////////////////// //restart - SPtr<UbScheduler> rSch(new UbScheduler(restartStep, restartStepStart)); - RestartCoProcessor rp(grid, rSch, comm, pathname, RestartCoProcessor::TXT); + SPtr<UbScheduler> rSch(new UbScheduler(cpStep, cpStart)); + SPtr<MPIIORestartCoProcessor> restartCoProcessor(new MPIIORestartCoProcessor(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 (grid->getTimeStep() == 0) + //bounding box + double g_minX1 = boundingBox[0]; + double g_minX2 = boundingBox[1]; + double g_minX3 = boundingBox[2]; + + double g_maxX1 = boundingBox[3]; + double g_maxX2 = boundingBox[4]; + double g_maxX3 = boundingBox[5]; + + double blockLength = (double)blocknx[0]*deltaXcoarse; + + double channel_high = channelHigh; // g_maxX3-g_minX3; + double channel_high_LB = channel_high/deltaXcoarse; + ////////////////////////////////////////////////////////////////////////// + double nu_LB = (u_LB*channel_high_LB)/Re; + ////////////////////////////////////////////////////////////////////////// + if (myid == 0) { - if (myid == 0) UBLOG(logINFO, "new start..."); - restart = false; + UBLOG(logINFO, "Parameters:"); + UBLOG(logINFO, "Re = " << Re); + UBLOG(logINFO, "u_LB = " << u_LB); + UBLOG(logINFO, "rho_LB = " << rho_LB); + UBLOG(logINFO, "nu_LB = " << nu_LB); + UBLOG(logINFO, "dx coarse = " << deltaXcoarse << " m"); + UBLOG(logINFO, "dx fine = " << deltaXfine << " m"); + UBLOG(logINFO, "channel_high = " << channel_high << " m"); + UBLOG(logINFO, "channel_high_LB = " << channel_high_LB); + UBLOG(logINFO, "number of levels = " << refineLevel + 1); + UBLOG(logINFO, "number of processes = " << comm->getNumberOfProcesses()); + UBLOG(logINFO, "number of threads = " << numOfThreads); + UBLOG(logINFO, "path = " << pathOut); + UBLOG(logINFO, "Preprocess - start"); + } - double offsetMinX3 = pmL[2]; - - double offsetMaxX1 = pmL[0]*lengthFactor; - double offsetMaxX2 = pmL[1]*2.0; - double offsetMaxX3 = channelHigh; // pmL[2] + channelHigh; //DLR-F15 //pmL[2]*2.0; - - //bounding box - double g_minX1 = origin[0]; - double g_minX2 = origin[1]; - double g_minX3 = origin[2]; - - double g_maxX1 = origin[0] + offsetMaxX1; - double g_maxX2 = origin[1] + offsetMaxX2; - double g_maxX3 = origin[2] + offsetMaxX3; -////////////////////////////////////////////////////////////////////////// - double nx1_temp = floor((g_maxX1-g_minX1) /(deltaXcoarse*(double)blocknx[0])); - deltaXcoarse = (g_maxX1-g_minX1) /(nx1_temp*(double)blocknx[0]); + if (newStart) + { + if (myid == 0) UBLOG(logINFO, "new start..."); - g_maxX1 -= 0.5* deltaXcoarse; -////////////////////////////////////////////////////////////////////////// - double blockLength = (double)blocknx[0]*deltaXcoarse; + grid->setPeriodicX1(true); grid->setPeriodicX2(true); @@ -149,48 +188,29 @@ void run(string configname) 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(), pathname + "/geo/gridCube", WbWriterVtkXmlBinary::getInstance()); + if (myid == 0) GbSystem3D::writeGeoObject(gridCube.get(), pathOut + "/geo/gridCube", WbWriterVtkXmlBinary::getInstance()); GenBlocksGridVisitor genBlocks(gridCube); grid->accept(genBlocks); - double channel_high = channelHigh; // g_maxX3-g_minX3; - double channel_high_LB = channel_high/deltaXcoarse; -////////////////////////////////////////////////////////////////////////// - //nu_LB = 0.005; - nu_LB = (u_LB*channel_high_LB)/Re; -////////////////////////////////////////////////////////////////////////// - if (myid == 0) - { - UBLOG(logINFO, "Parameters:"); - UBLOG(logINFO, "Re = " << Re); - UBLOG(logINFO, "u_LB = " << u_LB); - UBLOG(logINFO, "rho_LB = " << rho_LB); - UBLOG(logINFO, "nu_LB = " << nu_LB); - UBLOG(logINFO, "dx coarse = " << deltaXcoarse << " m"); - UBLOG(logINFO, "dx fine = " << grid->getDeltaX(refineLevel) << " m"); - UBLOG(logINFO, "number of levels = " << refineLevel + 1); - UBLOG(logINFO, "number of processes = " << comm->getNumberOfProcesses()); - UBLOG(logINFO, "number of threads = " << numOfThreads); - UBLOG(logINFO, "path = " << pathname); - UBLOG(logINFO, "Preprocess - start"); - } + ////////////////////////////////////////////////////////////////////////// //refinement double blockLengthX3Fine = grid->getDeltaX(refineLevel) * blocknx[2]; + double refHight = 0.002; - GbCuboid3DPtr refineBoxTop(new GbCuboid3D(g_minX1-blockLength, g_minX2-blockLength, g_maxX3-blockLengthX3Fine, g_maxX1+blockLength, g_maxX2+blockLength, g_maxX3+blockLength)); - if (myid == 0) GbSystem3D::writeGeoObject(refineBoxTop.get(), pathname + "/geo/refineBoxTop", WbWriterVtkXmlASCII::getInstance()); + GbCuboid3DPtr refineBoxTop(new GbCuboid3D(g_minX1-blockLength, g_minX2-blockLength, g_maxX3-refHight, g_maxX1+blockLength, g_maxX2+blockLength, g_maxX3)); + if (myid == 0) GbSystem3D::writeGeoObject(refineBoxTop.get(), pathOut + "/geo/refineBoxTop", WbWriterVtkXmlASCII::getInstance()); //GbCuboid3DPtr refineBoxBottom(new GbCuboid3D(g_minX1-blockLength, g_minX2-blockLength, g_minX3, g_maxX1+blockLength, g_maxX2+blockLength, g_minX3+offsetMinX3+blockLengthX3Fine)); - GbCuboid3DPtr refineBoxBottom(new GbCuboid3D(g_minX1-blockLength, g_minX2-blockLength, g_minX3-blockLengthX3Fine, g_maxX1+blockLength, g_maxX2+blockLength, g_minX3+blockLengthX3Fine)); - if (myid == 0) GbSystem3D::writeGeoObject(refineBoxBottom.get(), pathname + "/geo/refineBoxBottom", WbWriterVtkXmlASCII::getInstance()); + GbCuboid3DPtr refineBoxBottom(new GbCuboid3D(g_minX1-blockLength, g_minX2-blockLength, g_minX3, g_maxX1+blockLength, g_maxX2+blockLength, g_minX3+refHight)); + if (myid == 0) GbSystem3D::writeGeoObject(refineBoxBottom.get(), pathOut + "/geo/refineBoxBottom", WbWriterVtkXmlASCII::getInstance()); if (refineLevel > 0) { if (myid == 0) UBLOG(logINFO, "Refinement - start"); - RefineCrossAndInsideGbObjectHelper refineHelper(grid, refineLevel); + RefineCrossAndInsideGbObjectHelper refineHelper(grid, refineLevel, comm); refineHelper.addGbObject(refineBoxTop, refineLevel); refineHelper.addGbObject(refineBoxBottom, refineLevel); refineHelper.refine(); @@ -200,49 +220,32 @@ void run(string configname) //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(), pathname+"/geo/addWallZmin", WbWriterVtkXmlASCII::getInstance()); + 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(), pathname+"/geo/addWallZmax", WbWriterVtkXmlASCII::getInstance()); + if (myid == 0) GbSystem3D::writeGeoObject(addWallZmax.get(), pathOut+"/geo/addWallZmax", WbWriterVtkXmlASCII::getInstance()); //wall interactors - int bbOption = 1; - D3Q27BoundaryConditionAdapterPtr bcNoSlip(new D3Q27NoSlipBCAdapter(bbOption)); - SPtr<D3Q27Interactor> addWallZminInt(new D3Q27Interactor(addWallZmin, grid, bcNoSlip, Interactor3D::SOLID)); - SPtr<D3Q27Interactor> addWallZmaxInt(new D3Q27Interactor(addWallZmax, grid, bcNoSlip, Interactor3D::SOLID)); - - //////////////////////////////////////////////// - //TEST - //SPtr<GbObject3D> testCube(new GbCuboid3D(g_minX1 + 2.0 * blockLength, g_minX2 + 2.0 * blockLength, g_minX3 + 5.0 * blockLength, - // g_minX1 + 3.0 * blockLength, g_minX2 + 3.0 * blockLength, g_minX3 + 6.0 * blockLength)); - //if (myid == 0) GbSystem3D::writeGeoObject(testCube.get(), pathname + "/geo/testCube", WbWriterVtkXmlBinary::getInstance()); - //SPtr<D3Q27Interactor> testCubeInt(new D3Q27Interactor(testCube, grid, bcNoSlip, Interactor3D::SOLID)); - /////////////////////////////////////////////// + SPtr<D3Q27Interactor> addWallZminInt(new D3Q27Interactor(addWallZmin, grid, noSlipBCAdapter, Interactor3D::SOLID)); + SPtr<D3Q27Interactor> addWallZmaxInt(new D3Q27Interactor(addWallZmax, grid, noSlipBCAdapter, Interactor3D::SOLID)); //////////////////////////////////////////// //METIS SPtr<Grid3DVisitor> metisVisitor(new MetisPartitioningGridVisitor(comm, MetisPartitioningGridVisitor::LevelBased, D3Q27System::BSW, MetisPartitioner::KWAY)); //////////////////////////////////////////// - //Zoltan - //SPtr<Grid3DVisitor> zoltanVisitor(new ZoltanPartitioningGridVisitor(comm, D3Q27System::BSW, 1)); - //grid->accept(zoltanVisitor); - /////delete solid blocks if (myid == 0) UBLOG(logINFO, "deleteSolidBlocks - start"); InteractorsHelper intHelper(grid, metisVisitor); intHelper.addInteractor(addWallZminInt); intHelper.addInteractor(addWallZmaxInt); - ////////////////////////////////////////////////////////////////////////// - //TEST - //intHelper.addInteractor(testCubeInt); - ////////////////////////////////////////////////////////////////////////// - intHelper.selectBlocks(); + intHelper.selectBlocks(); if (myid == 0) UBLOG(logINFO, "deleteSolidBlocks - end"); ////////////////////////////////////// - WriteBlocksSPtr<CoProcessor> ppblocks(new WriteBlocksCoProcessor(grid, SPtr<UbScheduler>(new UbScheduler(1)), pathname, WbWriterVtkXmlBinary::getInstance(), comm)); - ppblocks->process(0); - ppblocks.reset(); + { + WriteBlocksCoProcessor ppblocks(grid, SPtr<UbScheduler>(new UbScheduler(1)), pathOut, WbWriterVtkXmlBinary::getInstance(), comm); + ppblocks.process(0); + } unsigned long long numberOfBlocks = (unsigned long long)grid->getNumberOfBlocks(); int ghostLayer = 3; @@ -269,31 +272,6 @@ void run(string configname) UBLOG(logINFO, "Available memory per process = " << availMem << " bytes"); } - LBMKernel3DPtr kernel = LBMKernel3DPtr(new LBMKernelETD3Q27CCLB(blocknx[0], blocknx[1], blocknx[2], LBMKernelETD3Q27CCLB::NORMAL)); - - mu::Parser fctForcingX1; - fctForcingX1.SetExpr("Fx1"); - fctForcingX1.DefineConst("Fx1", 1.0e-6); - - kernel->setWithForcing(true); - - SPtr<BCProcessor> bcProc; - BoundaryConditionPtr noSlipBC; - - if (thinWall) - { - bcProc = SPtr<BCProcessor>(new D3Q27ETForThinWallBCProcessor()); - noSlipBC = BoundaryConditionPtr(new ThinWallNoSlipBoundaryCondition()); - } - else - { - bcProc = SPtr<BCProcessor>(new D3Q27ETBCProcessor()); - noSlipBC = BoundaryConditionPtr(new NoSlipBoundaryCondition()); - } - - bcProc->addBC(noSlipBC); - - kernel->setBCProcessor(bcProc); SetKernelBlockVisitor kernelVisitor(kernel, nu_LB, availMem, needMem); grid->accept(kernelVisitor); @@ -302,7 +280,7 @@ void run(string configname) //undef nodes for refinement if (refineLevel > 0) { - D3Q27SetUndefinedNodesBlockVisitor undefNodesVisitor; + SetUndefinedNodesBlockVisitor undefNodesVisitor; grid->accept(undefNodesVisitor); } @@ -311,11 +289,11 @@ void run(string configname) intHelper.setBC(); ////porous media - if(false) + if(true) { string samplePathname = pathGeo + "/" + sampleFilename; - GbVoxelMatrix3DPtr sample(new GbVoxelMatrix3D(pmNX[0], pmNX[1], pmNX[2], 0, lthreshold, uthreshold)); + SPtr<GbVoxelMatrix3D> sample(new GbVoxelMatrix3D(pmNX[0], pmNX[1], pmNX[2], 0, lthreshold, uthreshold)); if (rawFile) { sample->readBufferedMatrixFromRawFile<unsigned short>(samplePathname, GbVoxelMatrix3D::BigEndian); @@ -331,182 +309,81 @@ void run(string configname) int bounceBackOption = 1; bool vxFile = false; int i = 0; - for (int x = 0; x < lengthFactor; x+=2) - { - double offset = pmL[0] * (double)x; - //sample 0 - if (myid == 0) UBLOG(logINFO, "sample # " << i); - sample->setVoxelMatrixMininum(origin[0]+offset, origin[1], origin[2]); - Utilities::voxelMatrixDiscretisation(sample, pathname, myid, i, grid, bounceBackOption, vxFile); - i++; - - if (myid == 0) - { - 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()); - } - - //sample 1 - if (myid == 0) UBLOG(logINFO, "sample # " << i); - sample->setVoxelMatrixMininum(origin[0]+pmL[0]+offset, origin[1], origin[2]); - sample->mirrorX(); - Utilities::voxelMatrixDiscretisation(sample, pathname, myid, i, grid, bounceBackOption, vxFile); - i++; + //for (int x = 0; x < lengthFactor; x+=2) + int lenX = (int)((g_maxX1-g_minX1)/(pmL[0])); + int lenY = (int)((g_maxX2-g_minX2)/(pmL[1])); - if (myid == 0) + for (int y = 0; y < lenY; y+=2) + for (int x = 0; x < lenX; x+=2) { - 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()); + double offsetX = pmL[0] * (double)x; + double offsetY = pmL[1] * (double)y; + //sample 0 + if (myid == 0) UBLOG(logINFO, "sample # " << i); + sample->setVoxelMatrixMininum(origin[0]+offsetX, origin[1]+offsetY, origin[2]); + Utilities::voxelMatrixDiscretisation(sample, pathOut, myid, i, grid, bounceBackOption, vxFile); + i++; + + if (myid == 0) + { + 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()); + } + + //sample 1 + if (myid == 0) UBLOG(logINFO, "sample # " << i); + sample->setVoxelMatrixMininum(origin[0]+pmL[0]+offsetX, origin[1]+offsetY, origin[2]); + sample->mirrorX(); + Utilities::voxelMatrixDiscretisation(sample, pathOut, myid, i, grid, bounceBackOption, vxFile); + i++; + + if (myid == 0) + { + 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()); + } + + //sample 2 + if (myid == 0) UBLOG(logINFO, "sample # " << i); + sample->setVoxelMatrixMininum(origin[0]+pmL[0]+offsetX, origin[1]+pmL[1]+offsetY, origin[2]); + sample->mirrorY(); + Utilities::voxelMatrixDiscretisation(sample, pathOut, myid, i, grid, bounceBackOption, vxFile); + i++; + + if (myid == 0) + { + 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()); + } + + //sample 3 + if (myid == 0) UBLOG(logINFO, "sample # " << i); + sample->setVoxelMatrixMininum(origin[0]+offsetX, origin[1]+pmL[1]+offsetY, origin[2]); + sample->mirrorX(); + Utilities::voxelMatrixDiscretisation(sample, pathOut, myid, i, grid, bounceBackOption, vxFile); + sample->mirrorY(); + i++; } - //sample 2 - if (myid == 0) UBLOG(logINFO, "sample # " << i); - sample->setVoxelMatrixMininum(origin[0]+pmL[0]+offset, origin[1]+pmL[1], origin[2]); - sample->mirrorY(); - Utilities::voxelMatrixDiscretisation(sample, pathname, myid, i, grid, bounceBackOption, vxFile); - i++; - - if (myid == 0) - { - 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()); - } - - //sample 3 - if (myid == 0) UBLOG(logINFO, "sample # " << i); - sample->setVoxelMatrixMininum(origin[0]+offset, origin[1]+pmL[1], origin[2]); - sample->mirrorX(); - Utilities::voxelMatrixDiscretisation(sample, pathname, myid, i, grid, bounceBackOption, vxFile); - sample->mirrorY(); - i++; - } - } - BoundaryConditionBlockVisitor bcVisitor; + grid->accept(bcVisitor); mu::Parser inflowProfileVx1, inflowProfileVx2, inflowProfileVx3, inflowProfileRho; inflowProfileVx1.SetExpr("x3 < h ? 0.0 : uLB+1*x1"); - //inflowProfileVx1.SetExpr("x3 < h ? 0.0 : uLB+1*x1-1*x2"); - ////inflowProfile.SetExpr("uLB+1*x1-1*x2"); - // //inflowProfile.SetExpr("uLB"); inflowProfileVx1.DefineConst("uLB", u_LB); - // //inflowProfile.DefineConst("uLB", 0.0116); - // inflowProfile.DefineConst("h", pmL[2]); - - //inflowProfile = Utilities::getDuctParaboloidX(g_maxX2/2.0, g_maxX2, pmL[2]+channelHigh/2.0, channelHigh, 0.1); - - //poiseuille flow - //double Cy = g_maxX2 / 2.0; - //double Hy = g_maxX2; - //double Cz = pmL[2] + channelHigh / 2.0; - //double Hz = channelHigh; - //double V = (9.0/4.0)*u_LB; - //inflowProfile.SetExpr("x3 < h ? 0.0 : V*(((-(x2-Cy)^2.0+(Hy/2.0)^2.0)/(Hy/2.0)^2.0)*((-(x3-Cz)^2.0+(Hz/2.0)^2.0)/(Hz/2.0)^2.0))"); - //inflowProfile.DefineConst("Cy", Cy); - //inflowProfile.DefineConst("Hy", Hy); - //inflowProfile.DefineConst("Cz", Cz); - //inflowProfile.DefineConst("Hz", Hz); - //inflowProfile.DefineConst("V", V); - //inflowProfile.DefineConst("h", pmL[2]); - - //log-law - double u_tau = 0.05;//sqrt(nu_LB * (u_LB/channel_high)); //0.0013252866627413104; - double z0 = 0.0001;//nu_LB / (9.0*u_tau); - double k = 0.4; - double nois = u_LB * 0.1; - //inflowProfile.SetExpr("2.3*u_tau/k*log(x3/z0)"); - //inflowProfile.SetExpr("x3 > 0 && (zMax-x3) > 0 ? (x3 < h ? 2.3*u_tau/k*log(x3/z0) : 2.3*u_tau/k*log((zMax-x3)/z0)) : 0"); - //inflowProfile.SetExpr("x3 > 0 && (zMax-x3) > 0 ? (x3 < h ? (1.0/k)*log(9.8*x3/u_tau)*u_tau : (1.0/k)*log(9.8*(zMax-x3)/u_tau)*u_tau ) : 0"); - //inflowProfile.SetExpr("x3 < h && x3 > 0 ? 2.3*u_tau/k*log(x3/z0) : 0.0"); - - //inflowProfile.SetExpr("Uref/log((Href+z0)/z0)*log((x3-zg+z0)/z0)"); - //inflowProfileVx1.SetExpr("x3 > 0 && (zMax-x3) > 0 ? (x3 < h ? Uref/log((Href+z0)/z0)*log((x3-zg+z0)/z0) : Uref/log((Href+z0)/z0)*log((zMax-x3-zg+z0)/z0)) : 0"); - - //inflowProfileVx1.DefineFun("rangeRandom1", rangeRandom1); - //inflowProfile.SetExpr("x3 < h ? Uref/log((Href+z0)/z0)*log((x3-zg+z0)/z0)+rangeRandom(-nois, nois) : Uref/log((Href+z0)/z0)*log((zMax-x3-zg+z0)/z0)+rangeRandom(-nois, nois)"); + inflowProfileVx1.DefineConst("h", pmL[2]); - //inflowProfileVx1.SetExpr("x3 < h ? Uref/log((Href+z0)/z0)*log((x3-zg+z0)/z0)+0.1*Uref*rangeRandom1() : Uref/log((Href+z0)/z0)*log((zMax-x3-zg+z0)/z0)+0.1*Uref*rangeRandom1()"); - - //inflowProfileVx1.SetExpr("x3 < h ? Uref/log((Href+z0)/z0)*log((x3-zg+z0)/z0)+ U*cos(8.0*PI*(x1)/(L1))*sin(8.0*PI*(x3)/L3) : Uref/log((Href+z0)/z0)*log((zMax-x3-zg+z0)/z0)+U*cos(8.0*PI*(x1)/(L1))*sin(8.0*PI*(x3)/L3)"); - //inflowProfileVx1.SetExpr("U*cos(4.0*PI*(x1)/(L1))*sin(4.0*PI*(x3)/L3)"); - - //inflowProfileVx1.SetExpr("x3 < h ? Uref/log((Href+z0)/z0)*log((x3-zg+z0)/z0)+U*cos(rangeRandom(2,8)*PI*x1/L2)*sin(rangeRandom(2,8)*PI*x2/L2)*sin(rangeRandom(2,8)*PI*x3/L3) : Uref/log((Href+z0)/z0)*log((zMax-x3-zg+z0)/z0)+U*cos(rangeRandom(2,8)*PI*x1/L2)*sin(rangeRandom(2,8)*PI*x2/L2)*sin(rangeRandom(2,8)*PI*x3/L3)"); - //inflowProfileVx1.SetExpr("U*cos(2.0*PI*x1/L1)*sin(2.0*PI*x2/L2)*sin(2.0*PI*x3/L3)"); - //inflowProfileVx1.SetExpr("U*sin(2.0*PI*x1/L2)*cos(1.0*PI*x2/L2)*cos(1.0*PI*x3/L3)"); - //inflowProfileVx1.SetExpr("U*cos(2.0*PI*(x1-L1/16.0)/(L1*2.0))*sin(2.0*PI*(x3-L3/4.0)/L3)"); - //inflowProfileVx1.SetExpr("U*sin(2.0*PI*x1/L2)"); - - inflowProfileVx1.DefineConst("U", u_LB); - inflowProfileVx1.DefineConst("PI", PI); - inflowProfileVx1.DefineConst("L1", g_maxX1-g_minX1); - inflowProfileVx1.DefineConst("L2", g_maxX2-g_minX2); - inflowProfileVx1.DefineConst("L3", g_maxX3-g_minX3); - - inflowProfileVx1.DefineConst("Uref", u_LB); - inflowProfileVx1.DefineConst("Href", channelHigh); - inflowProfileVx1.DefineConst("zg", 0.0); - inflowProfileVx1.DefineConst("nois", nois); - - inflowProfileVx1.DefineConst("u_tau", u_tau); - inflowProfileVx1.DefineConst("k", k); - inflowProfileVx1.DefineConst("z0", z0); - inflowProfileVx1.DefineConst("h", channelHigh / 2.0); - inflowProfileVx1.DefineConst("zMax", channelHigh); - - inflowProfileVx2.SetExpr("0.1*U*rangeRandom1()"); - //inflowProfileVx2.SetExpr("0.0"); - //inflowProfileVx2.SetExpr("-U*cos(2.0*PI*x1/L1)*sin(2.0*PI*x2/L2)*cos(2.0*PI*x3/L3)"); - //inflowProfileVx2.SetExpr("-U/2.0*sin(2.0*PI*x1/L1)*cos(2.0*PI*x2/L2)*sin(2.0*PI*x3/L3)"); - //inflowProfileVx2.SetExpr("-U/8.0*sin(rangeRandom(2,8)*PI*x1/L2)*cos(rangeRandom(2,8)*PI*x2/L2)*sin(rangeRandom(2,8)*PI*x3/L3)"); - inflowProfileVx2.DefineConst("U", u_LB); - inflowProfileVx2.DefineConst("PI", PI); - inflowProfileVx2.DefineConst("L1", g_maxX1-g_minX1); - inflowProfileVx2.DefineConst("L2", g_maxX2-g_minX2); - inflowProfileVx2.DefineConst("L3", g_maxX3-g_minX3); - inflowProfileVx2.DefineFun("rangeRandom1", rangeRandom1); - - inflowProfileVx3.SetExpr("0.1*U*rangeRandom1()"); - //inflowProfileVx3.SetExpr("-U/2.0*sin(8.0*PI*(x1)/(L1))*cos(8.0*PI*(x3)/L3)"); - //inflowProfileVx3.SetExpr("-U/2.0*sin(2.0*PI*(x1-L1/16.0)/(L1*2.0))*cos(2.0*PI*(x3-L3/4.0)/L3)"); - //inflowProfileVx3.SetExpr("-U/2.0*sin(2.0*PI*x1/L2)*sin(2.0*PI*x3/L3)"); - //inflowProfileVx3.SetExpr("-U*cos(1.0*PI*x1/L2)*sin(1.0*PI*x2/L2)*cos(1.0*PI*x3/L3)"); - //inflowProfileVx3.SetExpr("-U/2.0*sin(2.0*PI*x1/L2)*sin(2.0*PI*x2/L2)*cos(2.0*PI*x3/L3)"); - //inflowProfileVx3.SetExpr("-U/8.0*sin(rangeRandom(2,8)*PI*x1/L2)*sin(rangeRandom(2,8)*PI*x2/L2)*cos(rangeRandom(2,8)*PI*x3/L3)"); - inflowProfileVx3.DefineConst("U", u_LB); - inflowProfileVx3.DefineConst("PI", PI); - inflowProfileVx3.DefineConst("L1", g_maxX1-g_minX1); - inflowProfileVx3.DefineConst("L2", g_maxX2-g_minX2); - inflowProfileVx3.DefineConst("L3", g_maxX3-g_minX3); - inflowProfileVx3.DefineFun("rangeRandom1", rangeRandom1); - - inflowProfileRho.SetExpr("3.0*(rho/4.0*U^2*(sin(4.0*PI*(x1-L1/4.0)/L1)+sin(4.0*PI*(x3-L3/4.0)/L3)))"); - //inflowProfileRho.SetExpr("3.0*(rho/3.0+((rho*U^2)/16.0)*(cos(2.0*PI*x1/L2)+cos(2.0*PI*x2/L2))*(cos(2.0*PI*x3/L3)+2.0))"); - //inflowProfileRho.SetExpr("3.0*(rho/3.0+((rho*U^2)/32.0)*(cos(8.0*PI*x1/L2)+cos(8.0*PI*x2/L2))*(cos(8.0*PI*x3/L3)+2.0))"); - inflowProfileRho.DefineConst("U", u_LB); - inflowProfileRho.DefineConst("PI", PI); - inflowProfileRho.DefineConst("L1", g_maxX1-g_minX1); - inflowProfileRho.DefineConst("L2", g_maxX2-g_minX2); - inflowProfileRho.DefineConst("L3", g_maxX3-g_minX3); - inflowProfileRho.DefineConst("rho", rho_LB); - - D3Q27ETInitDistributionsBlockVisitor initVisitor(nu_LB, rho_LB); + InitDistributionsBlockVisitor initVisitor; initVisitor.setVx1(inflowProfileVx1); - //initVisitor.setVx2(inflowProfileVx2); - //initVisitor.setVx3(inflowProfileVx3); - //initVisitor.setRho(inflowProfileRho); grid->accept(initVisitor); - - ////set connectors - D3Q27InterpolationProcessorPtr iProcessor(new D3Q27IncompressibleOffsetInterpolationProcessor()); - D3Q27SetConnectorsBlockVisitor setConnsVisitor(comm, true, D3Q27System::ENDDIR, nu_LB, iProcessor); - //SPtr<ConnectorFactory> factory(new Block3DConnectorFactory()); - //ConnectorBlockVisitor setConnsVisitor(comm, nu_LB, iProcessor, factory); + InterpolationProcessorPtr iProcessor(new IncompressibleOffsetInterpolationProcessor()); + SetConnectorsBlockVisitor setConnsVisitor(comm, true, D3Q27System::ENDDIR, nu_LB, iProcessor); grid->accept(setConnsVisitor); //domain decomposition for threads @@ -514,183 +391,48 @@ void run(string configname) grid->accept(pqPartVisitor); //Postrozess - SPtr<UbScheduler> geoSch(new UbScheduler(1)); - MacroscopicQuantitiesSPtr<CoProcessor> ppgeo( - new MacroscopicQuantitiesCoProcessor(grid, geoSch, pathname, WbWriterVtkXmlBinary::getInstance(), conv, true)); - ppgeo->process(0); - ppgeo.reset(); - - coord[0] = g_minX1; - coord[1] = g_minX2; - coord[2] = g_minX3 + pmL[2]; - coord[3] = g_maxX1; - coord[4] = g_maxX2; - coord[5] = g_maxX3; - - //////////////////////////////////////////////////////// - FILE * pFile; - string str = pathname + "/checkpoints/coord.txt"; - pFile = fopen(str.c_str(), "w"); - fprintf(pFile, "%g\n", deltaXcoarse); - fprintf(pFile, "%g\n", nu_LB); - fprintf(pFile, "%g\n", coord[0]); - fprintf(pFile, "%g\n", coord[1]); - fprintf(pFile, "%g\n", coord[2]); - fprintf(pFile, "%g\n", coord[3]); - fprintf(pFile, "%g\n", coord[4]); - fprintf(pFile, "%g\n", coord[5]); - fclose(pFile); - //////////////////////////////////////////////////////// + { + SPtr<UbScheduler> geoSch(new UbScheduler(1)); + WriteBoundaryConditionsCoProcessor ppgeo(grid, geoSch, pathOut, WbWriterVtkXmlBinary::getInstance(), conv, comm); + ppgeo.process(0); + } if (myid == 0) UBLOG(logINFO, "Preprozess - end"); } else { - //////////////////////////////////////////////////////// - FILE * pFile; - string str = pathname + "/checkpoints/coord.txt"; - pFile = fopen(str.c_str(), "r"); - fscanf(pFile, "%lg\n", &deltaXcoarse); - fscanf(pFile, "%lg\n", &nu_LB); - fscanf(pFile, "%lg\n", &coord[0]); - fscanf(pFile, "%lg\n", &coord[1]); - fscanf(pFile, "%lg\n", &coord[2]); - fscanf(pFile, "%lg\n", &coord[3]); - fscanf(pFile, "%lg\n", &coord[4]); - fscanf(pFile, "%lg\n", &coord[5]); - fclose(pFile); - //////////////////////////////////////////////////////// - - if (myid == 0) - { - UBLOG(logINFO, "Parameters:"); - UBLOG(logINFO, "Re = " << Re); - UBLOG(logINFO, "u_LB = " << u_LB); - UBLOG(logINFO, "rho_LB = " << rho_LB); - UBLOG(logINFO, "nu_LB = " << nu_LB); - UBLOG(logINFO, "dx coarse = " << deltaXcoarse << " m"); - UBLOG(logINFO, "dx fine = " << grid->getDeltaX(refineLevel) << " m"); - UBLOG(logINFO, "number of levels = " << refineLevel + 1); - UBLOG(logINFO, "numOfThreads = " << numOfThreads); - UBLOG(logINFO, "path = " << pathname); - } - - ///////////////////////////////////////////////////////////// - ////bounding box - //double offsetMinX3 = pmL[2]; - - //double offsetMaxX1 = pmL[0]*lengthFactor; - //double offsetMaxX2 = pmL[1]*2.0; - //double offsetMaxX3 = channelHigh; - - //double g_minX1 = origin[0]; - //double g_minX2 = origin[1]; - //double g_minX3 = origin[2]; - - //double g_maxX1 = origin[0]+offsetMaxX1; - //double g_maxX2 = origin[1]+offsetMaxX2; - //double g_maxX3 = origin[2]+offsetMaxX3; - - //double blockLength = (double)blocknx[0]*deltaXcoarse; - - //GbCuboid3DPtr addWallZmin(new GbCuboid3D(g_minX1-blockLength, g_minX2-blockLength, g_minX3-blockLength, g_maxX1+blockLength, g_maxX2+blockLength, g_minX3+offsetMinX3)); - //if (myid==0) GbSystem3D::writeGeoObject(addWallZmin.get(), pathname+"/geo/addWallZmin", WbWriterVtkXmlASCII::getInstance()); - //int bbOption = 1; - //D3Q27BoundaryConditionAdapterPtr bcNoSlip(new D3Q27NoSlipBCAdapter(bbOption)); - //SPtr<D3Q27Interactor> addWallZminInt(new D3Q27Interactor(addWallZmin, grid, bcNoSlip, Interactor3D::SOLID)); - - //SetSolidOrTransBlockVisitor v1(addWallZminInt, SetSolidOrTransBlockVisitor::SOLID); - //grid->accept(v1); - //SetSolidOrTransBlockVisitor v2(addWallZminInt, SetSolidOrTransBlockVisitor::TRANS); - //grid->accept(v2); - - //std::vector<SPtr<Block3D>> blocks; - //std::vector<SPtr<Block3D>>& sb = addWallZminInt->getSolidBlockSet(); - //if (myid==0) UBLOG(logINFO, "number of solid blocks = "<<sb.size()); - ////blocks.insert(blocks.end(), sb.begin(), sb.end()); - //std::vector<SPtr<Block3D>>& tb = addWallZminInt->getTransBlockSet(); - //if (myid==0) UBLOG(logINFO, "number of trans blocks = "<<tb.size()); - //blocks.insert(blocks.end(), tb.begin(), tb.end()); - - //if (myid==0) UBLOG(logINFO, "number of blocks = "<<blocks.size()); - - //BOOST_FOREACH(SPtr<Block3D> block, blocks) - //{ - // block->setActive(true); - // addWallZminInt->setDifferencesToGbObject3D(block); - //} - - //////////////////////////////////////////////// - //////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(addWallZminInt); - ////intHelper.selectBlocks(); - ////if (myid==0) UBLOG(logINFO, "deleteSolidBlocks - end"); - ////////////////////////////////////////// - ////intHelper.setBC(); - //////////////////////////////////////////////////////////////// - - //BoundaryConditionBlockVisitor bcVisitor; - //grid->accept(bcVisitor); - - //WriteBlocksSPtr<CoProcessor> ppblocks(new WriteBlocksCoProcessor(grid, SPtr<UbScheduler>(new UbScheduler(1)), pathname, WbWriterVtkXmlBinary::getInstance(), comm)); - //ppblocks->process(0); - //ppblocks.reset(); - - //SPtr<UbScheduler> geoSch(new UbScheduler(1)); - //MacroscopicQuantitiesSPtr<CoProcessor> ppgeo( - // new MacroscopicQuantitiesCoProcessor(grid, geoSch, pathname, WbWriterVtkXmlBinary::getInstance(), conv, true)); - //ppgeo->process(0); - //ppgeo.reset(); - ////////////////////////////////////////////////////////////////// - - //set connectors - D3Q27InterpolationProcessorPtr iProcessor(new D3Q27IncompressibleOffsetInterpolationProcessor()); - D3Q27SetConnectorsBlockVisitor setConnsVisitor(comm, true, D3Q27System::ENDDIR, nu_LB, iProcessor); + //restartCoProcessor->restart((int)restartStep); + migCoProcessor->restart((int)restartStep); + grid->setTimeStep(restartStep); + //////////////////////////////////////////////////////////////////////////// + InterpolationProcessorPtr iProcessor(new CompressibleOffsetInterpolationProcessor()); + SetConnectorsBlockVisitor setConnsVisitor(comm, true, D3Q27System::ENDDIR, nu_LB, iProcessor); grid->accept(setConnsVisitor); - //domain decomposition for threads - PQueuePartitioningGridVisitor pqPartVisitor(numOfThreads); - grid->accept(pqPartVisitor); - - restart = true; + grid->accept(bcVisitor); if (myid == 0) UBLOG(logINFO, "Restart - end"); } - SPtr<UbScheduler> nupsSch(new UbScheduler(nupsSteps)); - NUPSCounterCoProcessor npr(grid, nupsSch, numOfThreads, comm); + + 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(outTime)); - MacroscopicQuantitiesCoProcessor pp(grid, stepSch, pathname, WbWriterVtkXmlBinary::getInstance(), conv); - - double startStep = grid->getTimeStep(); + SPtr<WriteMacroscopicQuantitiesCoProcessor> writeMQCoProcessor(new WriteMacroscopicQuantitiesCoProcessor(grid, stepSch, pathOut, WbWriterVtkXmlBinary::getInstance(), conv, comm)); - //SPtr<UbScheduler> visSch(new UbScheduler()); - //visSch->addSchedule(40000,40000,40000000); - //SPtr<UbScheduler> resSchRMS(new UbScheduler()); - //resSchRMS->addSchedule(40000, startStep, 40000000); - //SPtr<UbScheduler> resSchMeans(new UbScheduler()); - //resSchMeans->addSchedule(40000, startStep, 40000000); - //SPtr<UbScheduler> stepAvSch(new UbScheduler()); - //stepAvSch->addSchedule(100, 0, 10000000); - //AverageValuesCoProcessor Avpp(grid, pathname, WbWriterVtkXmlBinary::getInstance(), - // stepSch/*wann wird rausgeschrieben*/, stepAvSch/*wann wird gemittelt*/, resSchMeans, resSchRMS/*wann wird resettet*/, restart); + 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> AdjForcSch(new UbScheduler()); AdjForcSch->addSchedule(10, 0, 10000000); - D3Q27SPtr<IntegrateValuesHelper> intValHelp(new D3Q27IntegrateValuesHelper(grid, comm, - coord[0], coord[1], coord[2], - coord[3], coord[4], coord[5])); - if (myid == 0) GbSystem3D::writeGeoObject(intValHelp->getBoundingBox().get(), pathname + "/geo/IntValHelp", WbWriterVtkXmlBinary::getInstance()); + SPtr<IntegrateValuesHelper> intValHelp(new IntegrateValuesHelper(grid, comm, g_minX1, g_minX2, g_minX3, g_maxX1, g_maxX2, g_maxX3)); + if (myid == 0) GbSystem3D::writeGeoObject(intValHelp->getBoundingBox().get(), pathOut + "/geo/IntValHelp", WbWriterVtkXmlBinary::getInstance()); double vxTarget=u_LB; - AdjustForcingCoProcessor AdjForcPPPtr(grid, AdjForcSch, pathname, intValHelp, vxTarget, comm); + AdjustForcingCoProcessor AdjForcPPPtr(grid, AdjForcSch, pathOut, intValHelp, vxTarget, comm); //mu::Parser decrViscFunc; //decrViscFunc.SetExpr("nue0+c0/(t+1)/(t+1)"); @@ -700,15 +442,15 @@ void run(string configname) //DecrViscSch->addSchedule(10, 0, 1000); //DecreaseViscosityCoProcessor decrViscPPPtr(grid, DecrViscSch, &decrViscFunc, comm); - if (changeQs) - { - double z1 = pmL[2]; - D3Q27SPtr<IntegrateValuesHelper> intValHelp2(new D3Q27IntegrateValuesHelper(grid, comm, - coord[0], coord[1], z1 - deltaXfine, - coord[3], coord[4], z1 + deltaXfine)); - if (myid == 0) GbSystem3D::writeGeoObject(intValHelp2->getBoundingBox().get(), pathname + "/geo/intValHelp2", WbWriterVtkXmlBinary::getInstance()); - Utilities::ChangeRandomQs(intValHelp2); - } + //if (changeQs) + //{ + // double z1 = pmL[2]; + // SPtr<IntegrateValuesHelper> intValHelp2(new IntegrateValuesHelper(grid, comm, + // coord[0], coord[1], z1 - deltaXfine, + // coord[3], coord[4], z1 + deltaXfine)); + // if (myid == 0) GbSystem3D::writeGeoObject(intValHelp2->getBoundingBox().get(), pathOut + "/geo/intValHelp2", WbWriterVtkXmlBinary::getInstance()); + // Utilities::ChangeRandomQs(intValHelp2); + //} std::vector<double> levelCoords; std::vector<int> levels; @@ -735,21 +477,11 @@ void run(string configname) levels.push_back(0); levelCoords.push_back(0); levelCoords.push_back(0.002); - SPtr<UbScheduler> tavSch(new UbScheduler(1, timeAvStart, timeAvStop)); - TimeAveragedValuesSPtr<CoProcessor> tav(new TimeAveragedValuesCoProcessor(grid, pathname, WbWriterVtkXmlBinary::getInstance(), tavSch, comm, - TimeAveragedValuesCoProcessor::Velocity | TimeAveragedValuesCoProcessor::Fluctuations | TimeAveragedValuesCoProcessor::Triplecorrelations, - levels, levelCoords, bounds)); + //SPtr<UbScheduler> tavSch(new UbScheduler(1, timeAvStart, timeAvStop)); + //SPtr<CoProcessor> tav(new TimeAveragedValuesCoProcessor(grid, pathOut, WbWriterVtkXmlBinary::getInstance(), tavSch, comm, + // TimeAveragedValuesCoProcessor::Velocity | TimeAveragedValuesCoProcessor::Fluctuations | TimeAveragedValuesCoProcessor::Triplecorrelations, + // levels, levelCoords, bounds)); - if (averagingReset) - { - tav->reset(); - } - - //SPtr<UbScheduler> catalystSch(new UbScheduler(1)); - //InSituCatalystCoProcessor catalyst(grid, catalystSch, "pchannel.py"); - - //SPtr<UbScheduler> exitSch(new UbScheduler(10)); - //EmergencyExitCoProcessor exitCoProc(grid, exitSch, pathname, RestartSPtr<CoProcessor>(&rp), comm); //create line time series SPtr<UbScheduler> tpcSch(new UbScheduler(1,1,3)); @@ -757,8 +489,8 @@ void run(string configname) //GbPoint3DPtr p2(new GbPoint3D(0.064,0.005,0.01)); //SPtr<GbLine3D> line(new GbLine3D(p1.get(),p2.get())); SPtr<GbLine3D> line(new GbLine3D(new GbPoint3D(0.0,0.005,0.01),new GbPoint3D(0.064,0.005,0.01))); - LineTimeSeriesCoProcessor lineTs(grid, tpcSch,pathname+"/TimeSeries/line1.csv",line, 0,comm); - if (myid==0) lineTs.writeLine(pathname+"/geo/line1"); + LineTimeSeriesCoProcessor lineTs(grid, tpcSch,pathOut+"/TimeSeries/line1.csv",line, 0,comm); + if (myid==0) lineTs.writeLine(pathOut+"/geo/line1"); if (myid == 0) { @@ -767,14 +499,16 @@ void run(string configname) UBLOG(logINFO, "PID = " << myid << " Physical Memory currently used by current process: " << Utilities::getPhysMemUsedByMe()); } - const SPtr<ConcreteCalculatorFactory> calculatorFactory = std::make_shared<ConcreteCalculatorFactory>(stepSch); - CalculationManagerPtr calculation(new CalculationManager(grid, numOfThreads, endTime, calculatorFactory, CalculatorType::HYBRID)); - if (averaging) - { - calculation->setTimeAveragedValuesCoProcessor(tav); - } + omp_set_num_threads(numOfThreads); + SPtr<UbScheduler> stepGhostLayer(new UbScheduler(1)); + SPtr<Calculator> calculator(new BasicCalculator(grid, stepGhostLayer, endTime)); + calculator->addCoProcessor(nupsCoProcessor); + calculator->addCoProcessor(migCoProcessor); + calculator->addCoProcessor(writeMQSelectCoProcessor); + calculator->addCoProcessor(writeMQCoProcessor); + if (myid == 0) UBLOG(logINFO, "Simulation-start"); - calculation->calculate(); + calculator->calculate(); if (myid == 0) UBLOG(logINFO, "Simulation-end"); } catch (exception& e) diff --git a/source/CMake/compilerflags/icc170.cmake b/source/CMake/compilerflags/icc170.cmake index 6b22171c70d08170b73e47450936d9a73bbd0f00..8e68083485c3f2ef7a20f79465daaeb3712f5b7f 100644 --- a/source/CMake/compilerflags/icc170.cmake +++ b/source/CMake/compilerflags/icc170.cmake @@ -21,7 +21,7 @@ MACRO(SET_COMPILER_SPECIFIC_FLAGS_INTERN build_type use64BitOptions) #~ LIST(APPEND CAB_COMPILER_ADDTIONAL_C_COMPILER_FLAGS "-wd266") #function "__GKfree" declared implicitly #LIST(APPEND CAB_COMPILER_ADDTIONAL_CXX_COMPILER_FLAGS "-xHOST -O3 -ip -ipo -fno-alias -mcmodel=medium -qopt-streaming-stores=always") - LIST(APPEND CAB_COMPILER_ADDTIONAL_CXX_COMPILER_FLAGS "-xHOST -O3 -ip -fno-alias -mcmodel=medium -qopt-streaming-stores=always") + LIST(APPEND CAB_COMPILER_ADDTIONAL_CXX_COMPILER_FLAGS "-xHOST -O3 -ip -fno-alias -mcmodel=medium -qopt-streaming-stores=always -wd1478") #Debug #LIST(APPEND CAB_COMPILER_ADDTIONAL_CXX_COMPILER_FLAGS "-g -traceback") diff --git a/source/ThirdParty/Library/numerics/geometry3d/GbVoxelMatrix3D.cpp b/source/ThirdParty/Library/numerics/geometry3d/GbVoxelMatrix3D.cpp index 2ae644d0afaf53e7d94c8340369f27468bbedf2b..f7c77a87b2e257ac74e58b627198746899a5fc1f 100644 --- a/source/ThirdParty/Library/numerics/geometry3d/GbVoxelMatrix3D.cpp +++ b/source/ThirdParty/Library/numerics/geometry3d/GbVoxelMatrix3D.cpp @@ -209,12 +209,12 @@ long GbVoxelMatrix3D::getNumberOfFluid() /*=======================================================*/ double GbVoxelMatrix3D::getIntersectionRaytraceFactor(const double& x1, const double& x2, const double& x3, const double& rx1, const double& rx2, const double& rx3) { - if (!((UbMath::equal(rx1, 0.0)||UbMath::equal(fabs(rx1), 1.0)||UbMath::equal(fabs(rx1), UbMath::one_over_sqrt2)) - &&(UbMath::equal(rx2, 0.0)||UbMath::equal(fabs(rx2), 1.0)||UbMath::equal(fabs(rx2), UbMath::one_over_sqrt2)) - &&(UbMath::equal(rx3, 0.0)||UbMath::equal(fabs(rx3), 1.0)||UbMath::equal(fabs(rx3), UbMath::one_over_sqrt2)))) - { - throw UbException(UB_EXARGS, "nur fuer diskrete Boltzmannrichungen implementiert!!!"); - } + //if (!((UbMath::equal(rx1, 0.0)||UbMath::equal(fabs(rx1), 1.0)||UbMath::equal(fabs(rx1), UbMath::one_over_sqrt2)) + // &&(UbMath::equal(rx2, 0.0)||UbMath::equal(fabs(rx2), 1.0)||UbMath::equal(fabs(rx2), UbMath::one_over_sqrt2)) + // &&(UbMath::equal(rx3, 0.0)||UbMath::equal(fabs(rx3), 1.0)||UbMath::equal(fabs(rx3), UbMath::one_over_sqrt2)))) + //{ + // throw UbException(UB_EXARGS, "nur fuer diskrete Boltzmannrichungen implementiert!!!"); + //} //nachbarindex ermitteln int ndx1 = 0, ndx2 = 0, ndx3 = 0; diff --git a/source/VirtualFluids.h b/source/VirtualFluids.h index 42e7461152f331063c11ddb55ac7955542cadef5..451dba3bf10a0345053f4daffaae62aa94257093 100644 --- a/source/VirtualFluids.h +++ b/source/VirtualFluids.h @@ -272,7 +272,7 @@ #include <Visitors/SetKernelBlockVisitor.h> #include <Visitors/SetForcingBlockVisitor.h> #include <Visitors/SetSpongeLayerBlockVisitor.h> -#include <Visitors/SetSolidBlockVisitor.h> +#include <Visitors/SetSolidOrBoundaryBlockVisitor.h> #include <Visitors/RenumberBlockVisitor.h> #include <Visitors/ConnectorBlockVisitor.h> #include <Visitors/ViscosityBlockVisitor.h> diff --git a/source/VirtualFluidsCore/Grid/BasicCalculator.cpp b/source/VirtualFluidsCore/Grid/BasicCalculator.cpp index 2b3505fd1cab5e2a1d50b7ec24c9f529c04348b3..285c9cdbd0ced8042276de3ab647b80a16100af8 100644 --- a/source/VirtualFluidsCore/Grid/BasicCalculator.cpp +++ b/source/VirtualFluidsCore/Grid/BasicCalculator.cpp @@ -10,7 +10,7 @@ #ifdef _OPENMP #include <omp.h> #endif -#define OMP_SCHEDULE dynamic +#define OMP_SCHEDULE guided //#define TIMING //#include "UbTiming.h" diff --git a/source/VirtualFluidsCore/Interactors/InteractorsHelper.cpp b/source/VirtualFluidsCore/Interactors/InteractorsHelper.cpp index fda7c6c24153334d0460e17a4022aa6146519694..38ca4ffe62b1e9a4e1a2042ea1902505a2fd9a2a 100644 --- a/source/VirtualFluidsCore/Interactors/InteractorsHelper.cpp +++ b/source/VirtualFluidsCore/Interactors/InteractorsHelper.cpp @@ -1,6 +1,5 @@ #include "InteractorsHelper.h" -#include <SetSolidBlockVisitor.h> #include <Grid3DVisitor.h> #include <Grid3D.h> #include <Interactor3D.h> @@ -49,7 +48,7 @@ void InteractorsHelper::deleteSolidBlocks() { for(SPtr<Interactor3D> interactor : interactors) { - setBlocks(interactor, BlockType::SOLID); + setBlocks(interactor, SetSolidOrBoundaryBlockVisitor::SOLID); std::vector<SPtr<Block3D>>& sb = interactor->getSolidBlockSet(); solidBlocks.insert(solidBlocks.end(), sb.begin(), sb.end()); @@ -59,16 +58,16 @@ void InteractorsHelper::deleteSolidBlocks() updateGrid(); } ////////////////////////////////////////////////////////////////////////// -void InteractorsHelper::setBlocks(const SPtr<Interactor3D> interactor, BlockType type) const +void InteractorsHelper::setBlocks(const SPtr<Interactor3D> interactor, SetSolidOrBoundaryBlockVisitor::BlockType type) const { - SetSolidBlockVisitor v(interactor, type); + SetSolidOrBoundaryBlockVisitor v(interactor, type); grid->accept(v); } ////////////////////////////////////////////////////////////////////////// void InteractorsHelper::setBcBlocks() { for(const SPtr<Interactor3D> interactor : interactors) - setBlocks(interactor, BlockType::BC); + setBlocks(interactor, SetSolidOrBoundaryBlockVisitor::BC); } ////////////////////////////////////////////////////////////////////////// void InteractorsHelper::updateGrid() diff --git a/source/VirtualFluidsCore/Interactors/InteractorsHelper.h b/source/VirtualFluidsCore/Interactors/InteractorsHelper.h index 5f01691e1e3e7d2da547e186d736eb2529e507e3..d8407768ddb648713660263f62c61d362cfb3dc8 100644 --- a/source/VirtualFluidsCore/Interactors/InteractorsHelper.h +++ b/source/VirtualFluidsCore/Interactors/InteractorsHelper.h @@ -3,12 +3,13 @@ #include <vector> #include <PointerDefinitions.h> +#include <SetSolidOrBoundaryBlockVisitor.h> class Interactor3D; class Block3D; class Grid3D; class Grid3DVisitor; -enum class BlockType; +class SetSolidOrBoundaryBlockVisitor; class InteractorsHelper { @@ -19,12 +20,12 @@ public: void addInteractor(SPtr<Interactor3D> interactor); void selectBlocks(); void setBC(); - void sendDomainDecompositionVisitor() const; + void sendDomainDecompositionVisitor() const; protected: void deleteSolidBlocks(); - void setBlocks(const SPtr<Interactor3D> interactor, BlockType type) const; - void setBcBlocks(); + void setBlocks(const SPtr<Interactor3D> interactor, SetSolidOrBoundaryBlockVisitor::BlockType type) const; + void setBcBlocks(); private: void updateGrid(); diff --git a/source/VirtualFluidsCore/Utilities/VoxelMatrixUtil.hpp b/source/VirtualFluidsCore/Utilities/VoxelMatrixUtil.hpp index 094e644711eb30c3a6fe9645c00bb26a5e6a68f3..1ef7f15b475bae82e9ff331819544cb02c66ca1f 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 "SetSolidBlockVisitor.h" +#include "SetSolidOrBoundaryBlockVisitor.h" #include "Block3D.h" #include "Grid3D.h" @@ -20,30 +20,39 @@ namespace Utilities { if (myid == 0) matrix->writeToVTKImageDataASCII(pathname + "/geo/vmatrix" + UbSystem::toString(fileCounter)); } - - GbCuboid3DPtr vmBox(new GbCuboid3D(matrix->getX1Minimum(), matrix->getX2Minimum(), matrix->getX3Minimum(), matrix->getX1Maximum(), matrix->getX2Maximum(), matrix->getX3Maximum())); - if (myid == 0) GbSystem3D::writeGeoObject(vmBox.get(), pathname + "/geo/vmbox" + UbSystem::toString(fileCounter), WbWriterVtkXmlASCII::getInstance()); - SPtr<D3Q27Interactor> vmBoxInt = SPtr<D3Q27Interactor>(new D3Q27Interactor(vmBox, grid, noSlipPM, Interactor3D::SOLID)); - SetSolidBlockVisitor v1(vmBoxInt, BlockType::SOLID); - grid->accept(v1); - SetSolidBlockVisitor v2(vmBoxInt, BlockType::BC); - grid->accept(v2); - - std::vector<SPtr<Block3D>> blocks; - std::vector<SPtr<Block3D>>& sb = vmBoxInt->getSolidBlockSet(); - if (myid == 0) UBLOG(logINFO, "number of solid blocks = " << sb.size()); - blocks.insert(blocks.end(), sb.begin(), sb.end()); - std::vector<SPtr<Block3D>>& tb = vmBoxInt->getBcBlocks(); - if (myid == 0) UBLOG(logINFO, "number of trans blocks = " << tb.size()); - blocks.insert(blocks.end(), tb.begin(), tb.end()); - - if (myid == 0) UBLOG(logINFO, "number of blocks = " << blocks.size()); - - for(SPtr<Block3D> block : blocks) + else { - block->setActive(true); - vmInt->setDifferencesToGbObject3D(block); + GbCuboid3DPtr vmBox(new GbCuboid3D(matrix->getX1Minimum(), matrix->getX2Minimum(), matrix->getX3Minimum(), matrix->getX1Maximum(), matrix->getX2Maximum(), matrix->getX3Maximum())); + if (myid == 0) GbSystem3D::writeGeoObject(vmBox.get(), pathname + "/geo/vmbox" + UbSystem::toString(fileCounter), WbWriterVtkXmlASCII::getInstance()); } + + //GbCuboid3DPtr vmBox(new GbCuboid3D(matrix->getX1Minimum(), matrix->getX2Minimum(), matrix->getX3Minimum(), matrix->getX1Maximum(), matrix->getX2Maximum(), matrix->getX3Maximum())); + //if (myid == 0) GbSystem3D::writeGeoObject(vmBox.get(), pathname + "/geo/vmbox" + UbSystem::toString(fileCounter), WbWriterVtkXmlASCII::getInstance()); + //SPtr<D3Q27Interactor> vmBoxInt = SPtr<D3Q27Interactor>(new D3Q27Interactor(vmBox, grid, noSlipPM, Interactor3D::SOLID)); + //SetSolidOrBoundaryBlockVisitor v1(vmBoxInt, SetSolidOrBoundaryBlockVisitor::SOLID); + //grid->accept(v1); + //SetSolidOrBoundaryBlockVisitor v2(vmBoxInt, SetSolidOrBoundaryBlockVisitor::BC); + //grid->accept(v2); + + //std::vector<SPtr<Block3D>> blocks; + //std::vector<SPtr<Block3D>>& sb = vmBoxInt->getSolidBlockSet(); + //if (myid == 0) UBLOG(logINFO, "number of solid blocks = " << sb.size()); + //blocks.insert(blocks.end(), sb.begin(), sb.end()); + //std::vector<SPtr<Block3D>>& tb = vmBoxInt->getBcBlocks(); + //if (myid == 0) UBLOG(logINFO, "number of trans blocks = " << tb.size()); + //blocks.insert(blocks.end(), tb.begin(), tb.end()); + + //if (myid == 0) UBLOG(logINFO, "number of blocks = " << blocks.size()); + + //for(SPtr<Block3D> block : blocks) + //{ + // block->setActive(true); + // vmInt->setDifferencesToGbObject3D(block); + //} + + SetSolidOrBoundaryBlockVisitor v(vmInt, SetSolidOrBoundaryBlockVisitor::BC); + grid->accept(v); + vmInt->initInteractor(); } } #endif // VoxelMatrixUtil_h__ diff --git a/source/VirtualFluidsCore/Visitors/SetSolidBlockVisitor.cpp b/source/VirtualFluidsCore/Visitors/SetSolidOrBoundaryBlockVisitor.cpp similarity index 62% rename from source/VirtualFluidsCore/Visitors/SetSolidBlockVisitor.cpp rename to source/VirtualFluidsCore/Visitors/SetSolidOrBoundaryBlockVisitor.cpp index 8312f72070e5c67a92617061d85648219175324d..dab7b5f0443372663444fe5ffcf5612c99ad445a 100644 --- a/source/VirtualFluidsCore/Visitors/SetSolidBlockVisitor.cpp +++ b/source/VirtualFluidsCore/Visitors/SetSolidOrBoundaryBlockVisitor.cpp @@ -1,18 +1,18 @@ -#include "SetSolidBlockVisitor.h" +#include "SetSolidOrBoundaryBlockVisitor.h" #include "Interactor3D.h" #include "Grid3DSystem.h" #include "Grid3D.h" #include "Block3D.h" -SetSolidBlockVisitor::SetSolidBlockVisitor(SPtr<Interactor3D> interactor, BlockType type) : +SetSolidOrBoundaryBlockVisitor::SetSolidOrBoundaryBlockVisitor(SPtr<Interactor3D> interactor, BlockType type) : Block3DVisitor(0, Grid3DSystem::MAXLEVEL), interactor(interactor), type(type) { } -void SetSolidBlockVisitor::visit(SPtr<Grid3D> grid, SPtr<Block3D> block) +void SetSolidOrBoundaryBlockVisitor::visit(SPtr<Grid3D> grid, SPtr<Block3D> block) { if(block->getRank() == grid->getRank()) { @@ -20,10 +20,10 @@ void SetSolidBlockVisitor::visit(SPtr<Grid3D> grid, SPtr<Block3D> block) { switch (type) { - case BlockType::SOLID: + case SOLID: interactor->setSolidBlock(block); break; - case BlockType::BC: + case BC: interactor->setBCBlock(block); break; } diff --git a/source/VirtualFluidsCore/Visitors/SetSolidBlockVisitor.h b/source/VirtualFluidsCore/Visitors/SetSolidOrBoundaryBlockVisitor.h similarity index 54% rename from source/VirtualFluidsCore/Visitors/SetSolidBlockVisitor.h rename to source/VirtualFluidsCore/Visitors/SetSolidOrBoundaryBlockVisitor.h index e2a21d580fb874b5b9a1176299e302b0dcc81fc7..d65d12be74fb78169a23fa41d7206cddbe8f3f9c 100644 --- a/source/VirtualFluidsCore/Visitors/SetSolidBlockVisitor.h +++ b/source/VirtualFluidsCore/Visitors/SetSolidOrBoundaryBlockVisitor.h @@ -9,18 +9,18 @@ class Grid3D; class Block3D; class Interactor3D; -enum class BlockType { SOLID, BC }; - -class SetSolidBlockVisitor : public Block3DVisitor +class SetSolidOrBoundaryBlockVisitor : public Block3DVisitor { public: - SetSolidBlockVisitor(SPtr<Interactor3D> interactor, BlockType type); - virtual ~SetSolidBlockVisitor() {} + 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; + SPtr<Interactor3D> interactor; BlockType type; };