From 73734d9f53a33c34cad3fd0c583e69c28f77c055 Mon Sep 17 00:00:00 2001 From: Konstantin Kutscher <kutscher@irmb.tu-bs.de> Date: Fri, 2 Jun 2017 13:36:57 +0000 Subject: [PATCH] add MPICalculator --- .../DLR-F16/F16BombadilTest10e-6.cfg | 22 +- source/Applications/DLR-F16/f16.cpp | 86 +++-- source/Applications/LaminarTubeFlow/ltf.cfg | 8 +- source/Applications/LaminarTubeFlow/ltf.cpp | 11 +- source/Applications/mpi_benchmark/mpib.cfg | 3 + source/Applications/mpi_benchmark/mpib.cpp | 258 ++++++------- .../cmake_config_files/LOGIN07.config.cmake | 1 + .../cmake_config_files/LOGIN22.config.cmake | 1 + .../cmake_config_files/SUPERMUC.config.cmake | 23 ++ source/CMake/compilerflags/icc170.cmake | 39 ++ .../CoProcessors/MPIIORestartCoProcessor.cpp | 43 ++- .../CoProcessors/MPIIORestartCoProcessor.h | 27 +- .../WriteMacroscopicQuantitiesCoProcessor.cpp | 2 +- .../WriteMacroscopicQuantitiesCoProcessor.h | 27 +- .../Grid/CalculationManager.cpp | 27 +- .../Grid/CalculationManager.h | 4 +- source/VirtualFluidsCore/Grid/Calculator.h | 2 +- .../VirtualFluidsCore/Grid/MPICalculator.cpp | 359 ++++++++++++++++++ source/VirtualFluidsCore/Grid/MPICalculator.h | 45 +++ ...ssibleCumulantWithSpongeLayerLBMKernel.cpp | 109 +++++- ...ressibleCumulantWithSpongeLayerLBMKernel.h | 14 +- .../Parallel/MPICommunicator.cpp | 1 + 22 files changed, 872 insertions(+), 240 deletions(-) create mode 100644 source/CMake/cmake_config_files/LOGIN07.config.cmake create mode 100644 source/CMake/cmake_config_files/LOGIN22.config.cmake create mode 100644 source/CMake/cmake_config_files/SUPERMUC.config.cmake create mode 100644 source/CMake/compilerflags/icc170.cmake create mode 100644 source/VirtualFluidsCore/Grid/MPICalculator.cpp create mode 100644 source/VirtualFluidsCore/Grid/MPICalculator.h diff --git a/source/Applications/DLR-F16/F16BombadilTest10e-6.cfg b/source/Applications/DLR-F16/F16BombadilTest10e-6.cfg index 60ee8161c..257cda076 100644 --- a/source/Applications/DLR-F16/F16BombadilTest10e-6.cfg +++ b/source/Applications/DLR-F16/F16BombadilTest10e-6.cfg @@ -1,7 +1,7 @@ pathOut = d:/temp/DLR-F16 pathGeo = d:/Projects/SFB880/DLR-F16/A1_Forschungsdaten_Profilgeometrie_STL_CATIA_Rossian -fngFileWhole = f16-ascii.stl -#fngFileWhole = grundgeometrie-direkter-export.stl +#fngFileWhole = f16-ascii.stl +fngFileWhole = grundgeometrie-direkter-export.stl #fngFileWhole = grundgeometrie-mittel.stl fngFileBodyPart = f16-body-part-ascii.stl fngFileTrailingEdge = f16-trailing-edge-ascii.stl @@ -9,7 +9,7 @@ zigZagTape = 2zackenbaender0.stl numOfThreads = 4 availMem = 13e9 -refineLevel = 2 #10 +refineLevel = 4 #10 #blockNx = 7 8 8 #blockNx = 7 6 7 blockNx = 21 6 13 @@ -39,9 +39,9 @@ boundingBox = -0.90 1.20 0.035 0.065 -0.65 0.65 #deltaXfine = 0.005 #level 0 #deltaXfine = 0.0025 #level 1 -deltaXfine = 0.00125 #level 2 +#deltaXfine = 0.00125 #level 2 #deltaXfine = 0.000625 #level 3 -#deltaXfine = 0.0003125 #level 4 +deltaXfine = 0.0003125 #level 4 #deltaXfine = 0.00015625 #level 5 #deltaXfine = 0.000078125 #level 6 #deltaXfine = 0.0000390625 #level 7 @@ -50,16 +50,16 @@ deltaXfine = 0.00125 #level 2 #deltaXfine = 6.5e-6 startDistance = -1.0 -refineDistance = 0.3 +refineDistance = 5.0 #0.3 newStart = true restartStep = 10 -cpStart = 10 -cpStep = 10 +cpStart = 1000 +cpStep = 1000 -outTime = 10 -endTime = 10 +outTime = 1000 +endTime = 10000 logToFile = false @@ -68,4 +68,4 @@ porousTralingEdge = false thinWall = false -nupsStep = 10 10 10000000 \ No newline at end of file +nupsStep = 100 100 10000000 \ No newline at end of file diff --git a/source/Applications/DLR-F16/f16.cpp b/source/Applications/DLR-F16/f16.cpp index e7becc6be..82f3432dc 100644 --- a/source/Applications/DLR-F16/f16.cpp +++ b/source/Applications/DLR-F16/f16.cpp @@ -19,31 +19,57 @@ void run(string configname) ConfigurationFile config; config.load(configname); - string pathOut = config.getString("pathOut"); - string pathGeo = config.getString("pathGeo"); - string fngFileWhole = config.getString("fngFileWhole"); - string fngFileTrailingEdge = config.getString("fngFileTrailingEdge"); - string fngFileBodyPart = config.getString("fngFileBodyPart"); - string zigZagTape = config.getString("zigZagTape"); - int numOfThreads = config.getInt("numOfThreads"); + //string pathOut = config.getString("pathOut"); + //string pathGeo = config.getString("pathGeo"); + //string fngFileWhole = config.getString("fngFileWhole"); + //string fngFileTrailingEdge = config.getString("fngFileTrailingEdge"); + //string fngFileBodyPart = config.getString("fngFileBodyPart"); + //string zigZagTape = config.getString("zigZagTape"); + //int numOfThreads = config.getInt("numOfThreads"); + //vector<int> blockNx = config.getVector<int>("blockNx"); + //vector<double> boundingBox = config.getVector<double>("boundingBox"); + //double uLB = config.getDouble("uLB"); + //double restartStep = config.getDouble("restartStep"); + //double cpStart = config.getDouble("cpStart"); + //double cpStep = config.getDouble("cpStep"); + //double endTime = config.getDouble("endTime"); + //double outTime = config.getDouble("outTime"); + //double availMem = config.getDouble("availMem"); + //int refineLevel = config.getInt("refineLevel"); + //bool logToFile = config.getBool("logToFile"); + //bool porousTralingEdge = config.getBool("porousTralingEdge"); + //double deltaXfine = config.getDouble("deltaXfine")*1000.0; + //bool thinWall = config.getBool("thinWall"); + //double refineDistance = config.getDouble("refineDistance"); + //double startDistance = config.getDouble("startDistance"); + //vector<double> nupsStep = config.getVector<double>("nupsStep"); + //bool newStart = config.getBool("newStart"); + + string pathOut = config.getValue<string>("pathOut"); + string pathGeo = config.getValue<string>("pathGeo"); + string fngFileWhole = config.getValue<string>("fngFileWhole"); + string fngFileTrailingEdge = config.getValue<string>("fngFileTrailingEdge"); + string fngFileBodyPart = config.getValue<string>("fngFileBodyPart"); + string zigZagTape = config.getValue<string>("zigZagTape"); + int numOfThreads = config.getValue<int>("numOfThreads"); vector<int> blockNx = config.getVector<int>("blockNx"); vector<double> boundingBox = config.getVector<double>("boundingBox"); - double uLB = config.getDouble("uLB"); - double restartStep = config.getDouble("restartStep"); - double cpStart = config.getDouble("cpStart"); - double cpStep = config.getDouble("cpStep"); - double endTime = config.getDouble("endTime"); - double outTime = config.getDouble("outTime"); - double availMem = config.getDouble("availMem"); - int refineLevel = config.getInt("refineLevel"); - bool logToFile = config.getBool("logToFile"); - bool porousTralingEdge = config.getBool("porousTralingEdge"); - double deltaXfine = config.getDouble("deltaXfine")*1000.0; - bool thinWall = config.getBool("thinWall"); - double refineDistance = config.getDouble("refineDistance"); - double startDistance = config.getDouble("startDistance"); + double uLB = config.getValue<double>("uLB"); + double restartStep = config.getValue<double>("restartStep"); + double cpStart = config.getValue<double>("cpStart"); + double cpStep = config.getValue<double>("cpStep"); + double endTime = config.getValue<double>("endTime"); + double outTime = config.getValue<double>("outTime"); + double availMem = config.getValue<double>("availMem"); + int refineLevel = config.getValue<int>("refineLevel"); + bool logToFile = config.getValue<bool>("logToFile"); + bool porousTralingEdge = config.getValue<bool>("porousTralingEdge"); + double deltaXfine = config.getValue<double>("deltaXfine")*1000.0; + bool thinWall = config.getValue<bool>("thinWall"); + double refineDistance = config.getValue<double>("refineDistance"); + double startDistance = config.getValue<double>("startDistance"); vector<double> nupsStep = config.getVector<double>("nupsStep"); - bool newStart = config.getBool("newStart"); + bool newStart = config.getValue<bool>("newStart"); CommunicatorPtr comm = MPICommunicator::getInstance(); int myid = comm->getProcessID(); @@ -295,10 +321,10 @@ void run(string configname) fngIntrWhole = D3Q27TriFaceMeshInteractorPtr(new D3Q27TriFaceMeshInteractor(fngMeshWhole, grid, noSlipBCAdapter, Interactor3D::SOLID));//, Interactor3D::EDGES)); } - D3Q27TriFaceMeshInteractorPtr triBand1Interactor(new D3Q27TriFaceMeshInteractor(meshBand1, grid, noSlipBCAdapter, Interactor3D::SOLID));//, Interactor3D::EDGES)); - D3Q27TriFaceMeshInteractorPtr triBand2Interactor(new D3Q27TriFaceMeshInteractor(meshBand2, grid, noSlipBCAdapter, Interactor3D::SOLID));//, Interactor3D::EDGES)); - D3Q27TriFaceMeshInteractorPtr triBand3Interactor(new D3Q27TriFaceMeshInteractor(meshBand5, grid, noSlipBCAdapter, Interactor3D::SOLID));//, Interactor3D::EDGES)); - D3Q27TriFaceMeshInteractorPtr triBand4Interactor(new D3Q27TriFaceMeshInteractor(meshBand6, grid, noSlipBCAdapter, Interactor3D::SOLID));//, Interactor3D::EDGES)); + D3Q27TriFaceMeshInteractorPtr triBand1Interactor(new D3Q27TriFaceMeshInteractor(meshBand1, grid, noSlipBCAdapter, Interactor3D::SOLID, Interactor3D::EDGES)); + D3Q27TriFaceMeshInteractorPtr triBand2Interactor(new D3Q27TriFaceMeshInteractor(meshBand2, grid, noSlipBCAdapter, Interactor3D::SOLID, Interactor3D::EDGES)); + D3Q27TriFaceMeshInteractorPtr triBand3Interactor(new D3Q27TriFaceMeshInteractor(meshBand5, grid, noSlipBCAdapter, Interactor3D::SOLID, Interactor3D::EDGES)); + D3Q27TriFaceMeshInteractorPtr triBand4Interactor(new D3Q27TriFaceMeshInteractor(meshBand6, grid, noSlipBCAdapter, Interactor3D::SOLID, Interactor3D::EDGES)); if (refineLevel > 0 && myid == 0) { @@ -575,17 +601,25 @@ void run(string configname) SetKernelBlockVisitor kernelVisitor(kernel, nuLB, availMem, needMem); grid->accept(kernelVisitor); + UBLOG(logINFO, "SetKernelBlockVisitor:end"); + if (refineLevel > 0) { SetUndefinedNodesBlockVisitor undefNodesVisitor; grid->accept(undefNodesVisitor); } + UBLOG(logINFO, "SetUndefinedNodesBlockVisitor:end"); + //BC intHelper.setBC(); + UBLOG(logINFO, "intHelper.setBC():end"); + grid->accept(bcVisitor); + UBLOG(logINFO, "grid->accept(bcVisitor);:end"); + //initialization of distributions mu::Parser inflowProfileVx1, inflowProfileVx2, inflowProfileVx3; inflowProfileVx1.SetExpr("U*rangeRandom1()"); diff --git a/source/Applications/LaminarTubeFlow/ltf.cfg b/source/Applications/LaminarTubeFlow/ltf.cfg index f2f887b03..8e4777dd0 100644 --- a/source/Applications/LaminarTubeFlow/ltf.cfg +++ b/source/Applications/LaminarTubeFlow/ltf.cfg @@ -3,8 +3,8 @@ numOfThreads = 1 availMem = 10e9 #Grid -length = 71 16 32 -blocknx = 8 8 8 +length = 84 6 26 +blocknx = 21 6 13 dx = 1 refineLevel = 0 @@ -13,8 +13,8 @@ uLB = 0.001 Re = 10 -outTime = 1 -endTime = 1000 +outTime = 100 +endTime = 10000 logToFile = false diff --git a/source/Applications/LaminarTubeFlow/ltf.cpp b/source/Applications/LaminarTubeFlow/ltf.cpp index 479af1565..cc817f2df 100644 --- a/source/Applications/LaminarTubeFlow/ltf.cpp +++ b/source/Applications/LaminarTubeFlow/ltf.cpp @@ -225,8 +225,8 @@ void run(string configname) LBMKernelPtr kernel; - kernel = LBMKernelPtr(new IncompressibleCumulantLBMKernel(blocknx[0], blocknx[1], blocknx[2], IncompressibleCumulantLBMKernel::NORMAL)); - //kernel = LBMKernelPtr(new CompressibleCumulantLBMKernel(blocknx[0], blocknx[1], blocknx[2], CompressibleCumulantLBMKernel::NORMAL)); + //kernel = LBMKernelPtr(new IncompressibleCumulantLBMKernel(blocknx[0], blocknx[1], blocknx[2], IncompressibleCumulantLBMKernel::NORMAL)); + kernel = LBMKernelPtr(new CompressibleCumulantLBMKernel(blocknx[0], blocknx[1], blocknx[2], CompressibleCumulantLBMKernel::NORMAL)); // BCProcessorPtr bcProc(new BCProcessor()); @@ -255,8 +255,8 @@ void run(string configname) grid->accept(initVisitor); //set connectors - InterpolationProcessorPtr iProcessor(new IncompressibleOffsetInterpolationProcessor()); - //InterpolationProcessorPtr iProcessor(new CompressibleOffsetInterpolationProcessor()); + //InterpolationProcessorPtr iProcessor(new IncompressibleOffsetInterpolationProcessor()); + InterpolationProcessorPtr iProcessor(new CompressibleOffsetInterpolationProcessor()); SetConnectorsBlockVisitor setConnsVisitor(comm, true, D3Q27System::ENDDIR, nuLB, iProcessor); //ConnectorFactoryPtr factory(new Block3DConnectorFactory()); //ConnectorBlockVisitor setConnsVisitor(comm, nuLB, iProcessor, factory); @@ -298,7 +298,8 @@ void run(string configname) grid->accept(bcVisitor); //set connectors - InterpolationProcessorPtr iProcessor(new IncompressibleOffsetInterpolationProcessor()); + //InterpolationProcessorPtr iProcessor(new IncompressibleOffsetInterpolationProcessor()); + InterpolationProcessorPtr iProcessor(new CompressibleOffsetInterpolationProcessor()); SetConnectorsBlockVisitor setConnsVisitor(comm, true, D3Q27System::ENDDIR, nuLB, iProcessor); grid->accept(setConnsVisitor); diff --git a/source/Applications/mpi_benchmark/mpib.cfg b/source/Applications/mpi_benchmark/mpib.cfg index 39d73af8a..dabc02c8d 100644 --- a/source/Applications/mpi_benchmark/mpib.cfg +++ b/source/Applications/mpi_benchmark/mpib.cfg @@ -6,5 +6,8 @@ blockNx = 8 8 8 logToFile = false oneD = true priorityQueue = false +cpStep = 20 +restart = flase +restartStep = 200 nupsStep = 10 10 100 endTime = 10 \ No newline at end of file diff --git a/source/Applications/mpi_benchmark/mpib.cpp b/source/Applications/mpi_benchmark/mpib.cpp index a52e29f6a..d1205929b 100644 --- a/source/Applications/mpi_benchmark/mpib.cpp +++ b/source/Applications/mpi_benchmark/mpib.cpp @@ -33,6 +33,9 @@ void run(string configname) bool output = config.getBool("output"); vector<double> nupsStep = config.getVector<double>("nupsStep"); bool priorityQueue = config.getBool("priorityQueue"); + bool restart = config.getBool("restart"); + double restartStep = config.getDouble("restartStep"); + double cpStep = config.getDouble("cpStep"); if (logToFile) { @@ -64,156 +67,154 @@ void run(string configname) } - double dx = 1; - double g_minX1, g_minX2, g_minX3, g_maxX1, g_maxX2, g_maxX3; - double factor = 1.0; + LBMReal uLB = 0.05; + LBMReal Re = 20.0; + LBMReal rhoLB = 0.0; + LBMReal nueLB = 0.05842; + + Grid3DPtr grid(new Grid3D(comm)); + + ////////////////////////////////////////////////////////////////////////// + //restart + UbSchedulerPtr rSch(new UbScheduler(cpStep,cpStep)); + MPIIORestartCoProcessor rcp(grid, rSch, pathOut, comm); - if (oneD) + if (restart) { - factor = comm->getNumberOfProcesses() * numOfThreads; - g_minX1 = 0; - g_minX2 = 0; - g_minX3 = 0; - - g_maxX1 = blockNx[0]*2.0 * factor; - g_maxX2 = blockNx[1]*2.0; - g_maxX3 = blockNx[2]*2.0; + rcp.restart((int)restartStep); } else { - factor = pow(comm->getNumberOfProcesses() * numOfThreads, 1.0/3.0); - g_minX1 = 0; - g_minX2 = 0; - g_minX3 = 0; - - g_maxX1 = blockNx[0]*2.0 * factor; - g_maxX2 = blockNx[1]*2.0 * factor; - g_maxX3 = blockNx[2]*2.0 * factor; - } + double dx = 1; + double g_minX1, g_minX2, g_minX3, g_maxX1, g_maxX2, g_maxX3; + double factor = 1.0; - LBMReal uLB = 0.05; - LBMReal Re = 20.0; - LBMReal rhoLB = 0.0; - LBMReal nueLB = 0.05842; + if (oneD) + { + factor = comm->getNumberOfProcesses() * numOfThreads; + g_minX1 = 0; + g_minX2 = 0; + g_minX3 = 0; + + g_maxX1 = blockNx[0]*2.0 * factor; + g_maxX2 = blockNx[1]*2.0; + g_maxX3 = blockNx[2]*2.0; + } + else + { + factor = pow(comm->getNumberOfProcesses() * numOfThreads, 1.0/3.0); + g_minX1 = 0; + g_minX2 = 0; + g_minX3 = 0; + + g_maxX1 = blockNx[0]*2.0 * factor; + g_maxX2 = blockNx[1]*2.0 * factor; + g_maxX3 = blockNx[2]*2.0 * factor; + } - LBMUnitConverterPtr conv = LBMUnitConverterPtr(new LBMUnitConverter()); + LBMUnitConverterPtr conv = LBMUnitConverterPtr(new LBMUnitConverter()); - Grid3DPtr grid(new Grid3D(comm)); - grid->setDeltaX(dx); - grid->setBlockNX(blockNx[0], blockNx[1], blockNx[2]); + grid->setDeltaX(dx); + grid->setBlockNX(blockNx[0], blockNx[1], blockNx[2]); - GbObject3DPtr gridCube(new GbCuboid3D(g_minX1, g_minX2, g_minX3, g_maxX1, g_maxX2, g_maxX3)); - if (myid==0&&output) GbSystem3D::writeGeoObject(gridCube.get(), pathOut+"/geo/gridCube", WbWriterVtkXmlASCII::getInstance()); - GenBlocksGridVisitor genBlocks(gridCube); - grid->accept(genBlocks); + GbObject3DPtr gridCube(new GbCuboid3D(g_minX1, g_minX2, g_minX3, g_maxX1, g_maxX2, g_maxX3)); + if (myid==0&&output) GbSystem3D::writeGeoObject(gridCube.get(), pathOut+"/geo/gridCube", WbWriterVtkXmlASCII::getInstance()); + GenBlocksGridVisitor genBlocks(gridCube); + grid->accept(genBlocks); - //grid->setPeriodicX1(true); - //grid->setPeriodicX2(true); - //grid->setPeriodicX3(true); + //grid->setPeriodicX1(true); + //grid->setPeriodicX2(true); + //grid->setPeriodicX3(true); - if (myid==0) - { - UBLOG(logINFO, "//////////////////////////////////////////////////////////////////////////"); - UBLOG(logINFO, "2. PID = "<<myid<<" Total Physical Memory (RAM): "<<Utilities::getTotalPhysMem()/1073741824.0<<" GB"); - UBLOG(logINFO, "2. PID = "<<myid<<" Physical Memory currently used: "<<Utilities::getPhysMemUsed()/1073741824.0<<" GB"); - UBLOG(logINFO, "2. PID = "<<myid<<" Physical Memory currently used by current process: "<<Utilities::getPhysMemUsedByMe()/1073741824.0<<" GB"); - UBLOG(logINFO, "//////////////////////////////////////////////////////////////////////////"); - } + if (myid==0) + { + UBLOG(logINFO, "//////////////////////////////////////////////////////////////////////////"); + UBLOG(logINFO, "2. PID = "<<myid<<" Total Physical Memory (RAM): "<<Utilities::getTotalPhysMem()/1073741824.0<<" GB"); + UBLOG(logINFO, "2. PID = "<<myid<<" Physical Memory currently used: "<<Utilities::getPhysMemUsed()/1073741824.0<<" GB"); + UBLOG(logINFO, "2. PID = "<<myid<<" Physical Memory currently used by current process: "<<Utilities::getPhysMemUsedByMe()/1073741824.0<<" GB"); + UBLOG(logINFO, "//////////////////////////////////////////////////////////////////////////"); + } - if (priorityQueue) - { - if (myid==0) UBLOG(logINFO, "MetisPartitioningGridVisitor:start"); - MetisPartitioningGridVisitor metisVisitor(comm, MetisPartitioningGridVisitor::LevelBased, D3Q27System::BSW); - grid->accept(metisVisitor); - if (myid==0) UBLOG(logINFO, "MetisPartitioningGridVisitor:end"); - - //set connectors - if (myid==0) UBLOG(logINFO, "SetConnectorsBlockVisitor:start"); - InterpolationProcessorPtr iProcessor(new IncompressibleOffsetInterpolationProcessor()); - SetConnectorsBlockVisitor setConnsVisitor(comm, true, D3Q27System::ENDDIR, nueLB, iProcessor); - grid->accept(setConnsVisitor); - if (myid==0) UBLOG(logINFO, "SetConnectorsBlockVisitor:end"); - - //domain decomposition for threads - if (myid==0) UBLOG(logINFO, "PQueuePartitioningGridVisitor:start"); - PQueuePartitioningGridVisitor pqPartVisitor(numOfThreads); - grid->accept(pqPartVisitor); - if (myid==0) UBLOG(logINFO, "PQueuePartitioningGridVisitor:end"); - } - else - { - if (myid==0) UBLOG(logINFO, "MetisPartitioningGridVisitor:start"); - MetisPartitioningGridVisitor metisVisitor(comm, MetisPartitioningGridVisitor::LevelBased, D3Q27System::BSW, MetisPartitioner::KWAY, true, numOfThreads); - grid->accept(metisVisitor); - if (myid==0) UBLOG(logINFO, "MetisPartitioningGridVisitor:end"); - - //set connectors - if (myid==0) UBLOG(logINFO, "SetConnectorsBlockVisitor:start"); - InterpolationProcessorPtr iProcessor(new IncompressibleOffsetInterpolationProcessor()); - SetConnectorsBlockVisitor setConnsVisitor(comm, true, D3Q27System::ENDDIR, nueLB, iProcessor); - grid->accept(setConnsVisitor); - if (myid==0) UBLOG(logINFO, "SetConnectorsBlockVisitor:end"); - } + if (priorityQueue) + { + if (myid==0) UBLOG(logINFO, "MetisPartitioningGridVisitor:start"); + MetisPartitioningGridVisitor metisVisitor(comm, MetisPartitioningGridVisitor::LevelBased, D3Q27System::BSW); + grid->accept(metisVisitor); + if (myid==0) UBLOG(logINFO, "MetisPartitioningGridVisitor:end"); + + //domain decomposition for threads + if (myid==0) UBLOG(logINFO, "PQueuePartitioningGridVisitor:start"); + PQueuePartitioningGridVisitor pqPartVisitor(numOfThreads); + grid->accept(pqPartVisitor); + if (myid==0) UBLOG(logINFO, "PQueuePartitioningGridVisitor:end"); + } + else + { + if (myid==0) UBLOG(logINFO, "MetisPartitioningGridVisitor:start"); + MetisPartitioningGridVisitor metisVisitor(comm, MetisPartitioningGridVisitor::LevelBased, D3Q27System::BSW, MetisPartitioner::KWAY, true, numOfThreads); + grid->accept(metisVisitor); + if (myid==0) UBLOG(logINFO, "MetisPartitioningGridVisitor:end"); + } - if (output) - { - WriteBlocksCoProcessor ppblocks(grid, UbSchedulerPtr(new UbScheduler(1)), pathOut, WbWriterVtkXmlBinary::getInstance(), comm); - ppblocks.process(0); - } + if (output) + { + WriteBlocksCoProcessor ppblocks(grid, UbSchedulerPtr(new UbScheduler(1)), pathOut, WbWriterVtkXmlBinary::getInstance(), comm); + ppblocks.process(0); + } - unsigned long nob = grid->getNumberOfBlocks(); - int gl = 3; - unsigned long nodb = (blockNx[0])* (blockNx[1])* (blockNx[2]); - unsigned long nod = nob * (blockNx[0])* (blockNx[1])* (blockNx[2]); - unsigned long nodg = nob * (blockNx[0]+gl) * (blockNx[1]+gl) * (blockNx[2]+gl); - double needMemAll = double(nodg*(27*sizeof(double)+sizeof(int)+sizeof(float)*4)); - double needMem = needMemAll/double(comm->getNumberOfProcesses()); + unsigned long nob = grid->getNumberOfBlocks(); + int gl = 3; + unsigned long nodb = (blockNx[0])* (blockNx[1])* (blockNx[2]); + unsigned long nod = nob * (blockNx[0])* (blockNx[1])* (blockNx[2]); + unsigned long nodg = nob * (blockNx[0]+gl) * (blockNx[1]+gl) * (blockNx[2]+gl); + double needMemAll = double(nodg*(27*sizeof(double)+sizeof(int)+sizeof(float)*4)); + double needMem = needMemAll/double(comm->getNumberOfProcesses()); - if (myid==0) - { - UBLOG(logINFO, "//////////////////////////////////////////////////////////////////////////"); - UBLOG(logINFO, "Setup information:"); - UBLOG(logINFO, "Size of block = "<<blockNx[0]<<" x "<<blockNx[1]<<" x "<<blockNx[2] <<" nodes"); - UBLOG(logINFO, "Size of domain = "<<g_maxX1<<" x "<<g_maxX2<<" x "<<g_maxX3<<" dx "); - UBLOG(logINFO, "Number of blocks = "<<nob); - UBLOG(logINFO, "Number of nodes = "<<nod); - int minInitLevel = grid->getCoarsestInitializedLevel(); - int maxInitLevel = grid->getFinestInitializedLevel(); - for (int level = minInitLevel; level<=maxInitLevel; level++) + if (myid==0) { - int nobl = grid->getNumberOfBlocks(level); - UBLOG(logINFO, "Number of blocks for level "<<level<<" = "<<nob); - UBLOG(logINFO, "Number of nodes for level "<<level<<" = "<<nob*nodb); + UBLOG(logINFO, "//////////////////////////////////////////////////////////////////////////"); + UBLOG(logINFO, "Setup information:"); + UBLOG(logINFO, "Size of block = "<<blockNx[0]<<" x "<<blockNx[1]<<" x "<<blockNx[2]<<" nodes"); + UBLOG(logINFO, "Size of domain = "<<g_maxX1<<" x "<<g_maxX2<<" x "<<g_maxX3<<" dx "); + UBLOG(logINFO, "Number of blocks = "<<nob); + UBLOG(logINFO, "Number of nodes = "<<nod); + 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<<" = "<<nob); + UBLOG(logINFO, "Number of nodes for level "<<level<<" = "<<nob*nodb); + } + UBLOG(logINFO, "//////////////////////////////////////////////////////////////////////////"); + UBLOG(logINFO, "Necessary memory = "<<needMemAll/1073741824.0<<" GB"); + UBLOG(logINFO, "Necessary memory per process = "<<needMem/1073741824.0<<" GB"); + UBLOG(logINFO, "Available memory per process = "<<availMem/1073741824.0<<" GB"); + UBLOG(logINFO, "//////////////////////////////////////////////////////////////////////////"); } - UBLOG(logINFO, "//////////////////////////////////////////////////////////////////////////"); - UBLOG(logINFO, "Necessary memory = "<<needMemAll/1073741824.0<<" GB"); - UBLOG(logINFO, "Necessary memory per process = "<<needMem/1073741824.0<<" GB"); - UBLOG(logINFO, "Available memory per process = "<<availMem/1073741824.0<<" GB"); - UBLOG(logINFO, "//////////////////////////////////////////////////////////////////////////"); - } - LBMKernelPtr kernel; - kernel = LBMKernelPtr(new IncompressibleCumulantLBMKernel(blockNx[0], blockNx[1], blockNx[2], IncompressibleCumulantLBMKernel::NORMAL)); + LBMKernelPtr kernel; + kernel = LBMKernelPtr(new IncompressibleCumulantLBMKernel(blockNx[0], blockNx[1], blockNx[2], IncompressibleCumulantLBMKernel::NORMAL)); - BCProcessorPtr bcProc(new BCProcessor()); - kernel->setBCProcessor(bcProc); + BCProcessorPtr bcProc(new BCProcessor()); + kernel->setBCProcessor(bcProc); - SetKernelBlockVisitor kernelVisitor(kernel, nueLB, availMem, needMem); - grid->accept(kernelVisitor); + SetKernelBlockVisitor kernelVisitor(kernel, nueLB, availMem, needMem); + grid->accept(kernelVisitor); - //initialization of distributions - InitDistributionsBlockVisitor initVisitor(nueLB, rhoLB); - initVisitor.setVx1(uLB); - grid->accept(initVisitor); + //initialization of distributions + InitDistributionsBlockVisitor initVisitor(nueLB, rhoLB); + initVisitor.setVx1(uLB); + grid->accept(initVisitor); + } - ////////////////////////////////////////////////////////////////////////// - //restart - UbSchedulerPtr rSch(new UbScheduler(20, 20)); - //RestartCoProcessor rp(grid, rSch, comm, pathOut, RestartCoProcessor::BINARY); - MPIIORestartCoProcessor rcp(grid, rSch, pathOut, comm); - rcp.restart(10); - ////////////////////////////////////////////////////////////////////////// + //set connectors + if (myid==0) UBLOG(logINFO, "SetConnectorsBlockVisitor:start"); + InterpolationProcessorPtr iProcessor(new IncompressibleOffsetInterpolationProcessor()); + SetConnectorsBlockVisitor setConnsVisitor(comm, true, D3Q27System::ENDDIR, nueLB, iProcessor); + grid->accept(setConnsVisitor); + if (myid==0) UBLOG(logINFO, "SetConnectorsBlockVisitor:end"); UbSchedulerPtr nupsSch(new UbScheduler(nupsStep[0], nupsStep[1], nupsStep[2])); NUPSCounterCoProcessor npr(grid, nupsSch, numOfThreads, comm); @@ -230,7 +231,7 @@ void run(string configname) UBLOG(logINFO, "//////////////////////////////////////////////////////////////////////////"); } - CalculationManagerPtr calculation(new CalculationManager(grid, numOfThreads, endTime, visSch)); + CalculationManagerPtr calculation(new CalculationManager(grid, numOfThreads, endTime, visSch,CalculationManager::MPI)); if (myid==0) UBLOG(logINFO, "Simulation-start"); calculation->calculate(); if (myid==0) UBLOG(logINFO, "Simulation-end"); @@ -260,6 +261,7 @@ int main(int argc, char* argv[]) if (argv[1]!=NULL) { run(string(argv[1])); + UBLOG(logINFO, "run end"); } else { diff --git a/source/CMake/cmake_config_files/LOGIN07.config.cmake b/source/CMake/cmake_config_files/LOGIN07.config.cmake new file mode 100644 index 000000000..e19a785b2 --- /dev/null +++ b/source/CMake/cmake_config_files/LOGIN07.config.cmake @@ -0,0 +1 @@ +INCLUDE("CMake/cmake_config_files/SUPERMUC.config.cmake") diff --git a/source/CMake/cmake_config_files/LOGIN22.config.cmake b/source/CMake/cmake_config_files/LOGIN22.config.cmake new file mode 100644 index 000000000..e19a785b2 --- /dev/null +++ b/source/CMake/cmake_config_files/LOGIN22.config.cmake @@ -0,0 +1 @@ +INCLUDE("CMake/cmake_config_files/SUPERMUC.config.cmake") diff --git a/source/CMake/cmake_config_files/SUPERMUC.config.cmake b/source/CMake/cmake_config_files/SUPERMUC.config.cmake new file mode 100644 index 000000000..91316ea69 --- /dev/null +++ b/source/CMake/cmake_config_files/SUPERMUC.config.cmake @@ -0,0 +1,23 @@ +LIST(APPEND CAB_ADDTIONAL_COMPILER_FLAGS -D__unix__) +LIST(APPEND CAB_ADDTIONAL_COMPILER_FLAGS -D__UNIX__) + +######################################################################################## +## BOOST ALLGMEINGUELTIG ## +######################################################################################## +#standard boost +SET(BOOST_VERSION "1.61.0" CACHE STRING "std: 1.61.0") +SET(BOOST_INCLUDEDIR "/lrz/sys/libraries/boost/1.61_icc/include") +SET(BOOST_LIBRARYDIR "/lrz/sys/libraries/boost/1.61_icc/lib") + +#IF(BOOST_VERSION AND NOT BOOST_INCLUDEDIR) +# MESSAGE("${BOOST_VERSION} not found on ${CAB_MACHINE} for specified compiler") +#ENDIF() + +################################################################################# +# METIS +################################################################################# +IF(${USE_METIS}) + SET(METIS_INCLUDEDIR "/lrz/sys/libraries/metis/5.1.0/i4r4/include") + SET(METIS_DEBUG_LIBRARY "/lrz/sys/libraries/metis/5.1.0/i4r4/lib/libmetis.a") + SET(METIS_RELEASE_LIBRARY "/lrz/sys/libraries/metis/5.1.0/i4r4/lib/libmetis.a") +ENDIF() diff --git a/source/CMake/compilerflags/icc170.cmake b/source/CMake/compilerflags/icc170.cmake new file mode 100644 index 000000000..d14475465 --- /dev/null +++ b/source/CMake/compilerflags/icc170.cmake @@ -0,0 +1,39 @@ +############################################################################################################### +## +## intel170 +## +############################################################################################################### + +MACRO(SET_COMPILER_SPECIFIC_FLAGS_INTERN build_type use64BitOptions) + + #~ IF( ${use64BitOptions} ) + #~ LIST(APPEND CAB_COMPILER_ADDTIONAL_CXX_COMPILER_FLAGS "-D__amd64" ) + #~ ENDIF() + + #LIST(APPEND CAB_COMPILER_ADDTIONAL_CXX_COMPILER_FLAGS "-O") + #~ LIST(APPEND CAB_COMPILER_ADDTIONAL_CXX_COMPILER_FLAGS "-wd654") + #~ LIST(APPEND CAB_COMPILER_ADDTIONAL_CXX_COMPILER_FLAGS "-wd1125") #virtual function override intended + #~ LIST(APPEND CAB_COMPILER_ADDTIONAL_CXX_COMPILER_FLAGS "-wd1224") #warning directive: This file includes at least one deprecated or antiquated header + #~ LIST(APPEND CAB_COMPILER_ADDTIONAL_CXX_COMPILER_FLAGS "-wd377") #class "std::auto_ptr<RCF::I_ClientTransport>" has no suitable copy constructor + #~ LIST(APPEND CAB_COMPILER_ADDTIONAL_CXX_COMPILER_FLAGS "-wd327") #class "std::auto_ptr<RCF::I_ClientTransport>" has no suitable copy constructor + #~ LIST(APPEND CAB_COMPILER_ADDTIONAL_CXX_COMPILER_FLAGS "-wd327") #class "std::auto_ptr<RCF::I_ClientTransport>" has no suitable copy constructor +#~ + #~ 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") + ############################################################################################################### + ## OpenMP support + ############################################################################################################### + IF(USE_OPENMP) + LIST(APPEND CAB_COMPILER_ADDTIONAL_CXX_COMPILER_FLAGS "-qopenmp") + ENDIF() + + + ############################################################################################################### + ## mt support + ############################################################################################################### + #LIST(APPEND CAB_COMPILER_ADDTIONAL_CXX_COMPILER_FLAGS "-pthread") + #LIST(APPEND CAB_COMPILER_ADDTIONAL_C_COMPILER_FLAGS "-pthread") + +ENDMACRO(SET_COMPILER_SPECIFIC_FLAGS_INTERN build_type use64BitOptions) diff --git a/source/VirtualFluidsCore/CoProcessors/MPIIORestartCoProcessor.cpp b/source/VirtualFluidsCore/CoProcessors/MPIIORestartCoProcessor.cpp index fca82f26c..27636cc36 100644 --- a/source/VirtualFluidsCore/CoProcessors/MPIIORestartCoProcessor.cpp +++ b/source/VirtualFluidsCore/CoProcessors/MPIIORestartCoProcessor.cpp @@ -16,7 +16,8 @@ MPIIORestartCoProcessor::MPIIORestartCoProcessor(Grid3DPtr grid, UbSchedulerPtr CommunicatorPtr comm) : CoProcessor(grid, s), path(path), - comm(comm) + comm(comm), + mpiTypeFreeFlag(false) { UbSystem::makeDirectory(path + "/mpi_io_cp"); @@ -98,11 +99,15 @@ MPIIORestartCoProcessor::~MPIIORestartCoProcessor() MPI_Type_free(&blockParamType); MPI_Type_free(&block3dType); MPI_Type_free(&dataSetType); - MPI_Type_free(&dataSetDoubleType); MPI_Type_free(&boundCondType); MPI_Type_free(&boundCondType1000); MPI_Type_free(&boundCondTypeAdd); - MPI_Type_free(&bcindexmatrixType); + + if (mpiTypeFreeFlag) + { + MPI_Type_free(&dataSetDoubleType); + MPI_Type_free(&bcindexmatrixType); + } } ////////////////////////////////////////////////////////////////////////// @@ -358,13 +363,13 @@ void MPIIORestartCoProcessor::writeBlocks(int step) { if(rank == 0) { - next_view_offset = view_offset + sizeof(GridParam) + sizeof(blockParam) + blocksCount * sizeof(Block3d); + next_view_offset = view_offset + sizeof(GridParam) + sizeof(BlockParam) + blocksCount * sizeof(Block3d); MPI_Send(&next_view_offset, 1, MPI_LONG_LONG_INT, 1, 5, MPI_COMM_WORLD); } else { MPI_Recv(&view_offset, 1, MPI_LONG_LONG_INT, rank - 1, 5, MPI_COMM_WORLD, MPI_STATUS_IGNORE); - next_view_offset = view_offset + sizeof(GridParam) + sizeof(blockParam) + blocksCount * sizeof(Block3d); + next_view_offset = view_offset + sizeof(GridParam) + sizeof(BlockParam) + blocksCount * sizeof(Block3d); if(rank < size - 1) MPI_Send(&next_view_offset, 1, MPI_LONG_LONG_INT, rank + 1, 5, MPI_COMM_WORLD); } @@ -377,7 +382,7 @@ void MPIIORestartCoProcessor::writeBlocks(int step) // each process writes common parameters of a block MPI_File_write_at_all(file_handler, view_offset + sizeof(GridParam), &blockParamStr, 1, blockParamType, MPI_STATUS_IGNORE); // each process writes it's blocks - MPI_File_write_at_all(file_handler, view_offset + sizeof(GridParam) + sizeof(blockParam), &block3dArray[0], blocksCount, block3dType, MPI_STATUS_IGNORE); + MPI_File_write_at_all(file_handler, view_offset + sizeof(GridParam) + sizeof(BlockParam), &block3dArray[0], blocksCount, block3dType, MPI_STATUS_IGNORE); //MPI_File_sync(file_handler); MPI_File_close(&file_handler); @@ -388,6 +393,8 @@ void MPIIORestartCoProcessor::writeBlocks(int step) MPI_Type_contiguous(blockParamStr.bcindexmatrix_count, MPI_INT, &bcindexmatrixType); MPI_Type_commit(&bcindexmatrixType); + mpiTypeFreeFlag = true; + delete[] block3dArray; delete gridParameters; } @@ -409,7 +416,7 @@ void MPIIORestartCoProcessor::writeDataSet(int step) blocksCount += static_cast<int>(blocksVector[level].size()); } - dataSet* dataSetArray = new dataSet[blocksCount]; + DataSet* dataSetArray = new DataSet[blocksCount]; std::vector<double> doubleValuesArray; // double-values (arrays of f's) in all blocks int ic = 0; @@ -473,13 +480,13 @@ void MPIIORestartCoProcessor::writeDataSet(int step) { if(rank == 0) { - next_view_offset = view_offset + blocksCount * (sizeof(dataSet) + blockParamStr.doubleCountInBlock * sizeof(double)); + next_view_offset = view_offset + blocksCount * (sizeof(DataSet) + blockParamStr.doubleCountInBlock * sizeof(double)); MPI_Send(&next_view_offset, 1, MPI_LONG_LONG_INT, 1, 5, MPI_COMM_WORLD); } else { MPI_Recv(&view_offset, 1, MPI_LONG_LONG_INT, rank - 1, 5, MPI_COMM_WORLD, MPI_STATUS_IGNORE); - next_view_offset = view_offset + blocksCount * (sizeof(dataSet) + blockParamStr.doubleCountInBlock * sizeof(double)); + next_view_offset = view_offset + blocksCount * (sizeof(DataSet) + blockParamStr.doubleCountInBlock * sizeof(double)); if(rank < size - 1) MPI_Send(&next_view_offset, 1, MPI_LONG_LONG_INT, rank + 1, 5, MPI_COMM_WORLD); } @@ -495,7 +502,7 @@ void MPIIORestartCoProcessor::writeDataSet(int step) // each process writes data identifying blocks MPI_File_write_at_all(file_handler, view_offset, &dataSetArray[0], blocksCount, dataSetType, MPI_STATUS_IGNORE); // each process writes the dataSet arrays - MPI_File_write_at_all(file_handler, view_offset + blocksCount * sizeof(dataSet), &doubleValuesArray[0], blocksCount, dataSetDoubleType, MPI_STATUS_IGNORE); + MPI_File_write_at_all(file_handler, view_offset + blocksCount * sizeof(DataSet), &doubleValuesArray[0], blocksCount, dataSetDoubleType, MPI_STATUS_IGNORE); //MPI_File_sync(file_handler); MPI_File_close(&file_handler); @@ -687,13 +694,13 @@ void MPIIORestartCoProcessor::readBlocks(int step) { if(rank == 0) { - next_read_offset = read_offset + sizeof(GridParam) + sizeof(blockParam) + blocksCount * sizeof(Block3d); + next_read_offset = read_offset + sizeof(GridParam) + sizeof(BlockParam) + blocksCount * sizeof(Block3d); MPI_Send(&next_read_offset, 1, MPI_LONG_LONG_INT, 1, 5, MPI_COMM_WORLD); } else { MPI_Recv(&read_offset, 1, MPI_LONG_LONG_INT, rank - 1, 5, MPI_COMM_WORLD, MPI_STATUS_IGNORE); - next_read_offset = read_offset + sizeof(GridParam) + sizeof(blockParam) + blocksCount * sizeof(Block3d); + next_read_offset = read_offset + sizeof(GridParam) + sizeof(BlockParam) + blocksCount * sizeof(Block3d); if(rank < size - 1) MPI_Send(&next_read_offset, 1, MPI_LONG_LONG_INT, rank + 1, 5, MPI_COMM_WORLD); } @@ -706,7 +713,7 @@ void MPIIORestartCoProcessor::readBlocks(int step) // read parameters of a block MPI_File_read_at_all(file_handler, read_offset + sizeof(GridParam), &blockParamStr, 1, blockParamType, MPI_STATUS_IGNORE); // read all the blocks - MPI_File_read_at_all(file_handler, read_offset + sizeof(GridParam) + sizeof(blockParam), &block3dArray[0], blocksCount, block3dType, MPI_STATUS_IGNORE); + MPI_File_read_at_all(file_handler, read_offset + sizeof(GridParam) + sizeof(BlockParam), &block3dArray[0], blocksCount, block3dType, MPI_STATUS_IGNORE); //MPI_File_sync(file_handler); MPI_File_close(&file_handler); @@ -810,6 +817,8 @@ void MPIIORestartCoProcessor::readBlocks(int step) MPI_Type_contiguous(blockParamStr.bcindexmatrix_count, MPI_INT, &bcindexmatrixType); MPI_Type_commit(&bcindexmatrixType); + mpiTypeFreeFlag = true; + delete gridParameters; delete [] block3dArray; } @@ -829,7 +838,7 @@ void MPIIORestartCoProcessor::readDataSet(int step) int blocksCount = 0; MPI_File_read_at(file_handler, rank * sizeof(int), &blocksCount, 1, MPI_INT, MPI_STATUS_IGNORE); - dataSet* dataSetArray = new dataSet[blocksCount]; + DataSet* dataSetArray = new DataSet[blocksCount]; std::vector<double> doubleValuesArray(blocksCount * blockParamStr.doubleCountInBlock); // double-values in all blocks // calculate the read offset @@ -840,20 +849,20 @@ void MPIIORestartCoProcessor::readDataSet(int step) { if(rank == 0) { - next_read_offset = read_offset + blocksCount * (sizeof(dataSet) + blockParamStr.doubleCountInBlock * sizeof(double)); + next_read_offset = read_offset + blocksCount * (sizeof(DataSet) + blockParamStr.doubleCountInBlock * sizeof(double)); MPI_Send(&next_read_offset, 1, MPI_LONG_LONG_INT, 1, 5, MPI_COMM_WORLD); } else { MPI_Recv(&read_offset, 1, MPI_LONG_LONG_INT, rank - 1, 5, MPI_COMM_WORLD, MPI_STATUS_IGNORE); - next_read_offset = read_offset + blocksCount * (sizeof(dataSet) + blockParamStr.doubleCountInBlock * sizeof(double)); + next_read_offset = read_offset + blocksCount * (sizeof(DataSet) + blockParamStr.doubleCountInBlock * sizeof(double)); if(rank < size - 1) MPI_Send(&next_read_offset, 1, MPI_LONG_LONG_INT, rank + 1, 5, MPI_COMM_WORLD); } } MPI_File_read_at_all(file_handler, read_offset, &dataSetArray[0], blocksCount, dataSetType, MPI_STATUS_IGNORE); - MPI_File_read_at_all(file_handler, read_offset + blocksCount * sizeof(dataSet), &doubleValuesArray[0], blocksCount, dataSetDoubleType, MPI_STATUS_IGNORE); + MPI_File_read_at_all(file_handler, read_offset + blocksCount * sizeof(DataSet), &doubleValuesArray[0], blocksCount, dataSetDoubleType, MPI_STATUS_IGNORE); //MPI_File_sync(file_handler); MPI_File_close(&file_handler); diff --git a/source/VirtualFluidsCore/CoProcessors/MPIIORestartCoProcessor.h b/source/VirtualFluidsCore/CoProcessors/MPIIORestartCoProcessor.h index 9609d76af..a2e970b0b 100644 --- a/source/VirtualFluidsCore/CoProcessors/MPIIORestartCoProcessor.h +++ b/source/VirtualFluidsCore/CoProcessors/MPIIORestartCoProcessor.h @@ -19,7 +19,7 @@ class MPIIORestartCoProcessor: public CoProcessor //! \struct GridParam //! \brief Structure describes parameters of the grid //! \details The structure is nessasary to restore the grid correctly - typedef struct GridParam + struct GridParam { double trafoParams[33]; double deltaX; @@ -34,12 +34,12 @@ class MPIIORestartCoProcessor: public CoProcessor bool periodicX3; bool active; bool transformation; - }GridParam; + }; //! \struct blockParam //! \brief Structure describes parameters of the block that are equal in all blocks //! \details The structure used to store some parameters needed to restore dataSet arrays - typedef struct blockParam + struct BlockParam { int nx1; // to find the right block int nx2; @@ -47,13 +47,13 @@ class MPIIORestartCoProcessor: public CoProcessor int nx[9][4]; // 9 arrays x (nx1, nx2, nx3, nx4) int doubleCountInBlock; // how many double-values are in all arrays dataSet in one (any) block int bcindexmatrix_count; // how many bcindexmatrix-values are in one (any) block - }blockParam; + }; //! \struct Block3d //! \brief Structure contains information of the block //! \details The structure is used to write the data describing the block in the grid when saving the grid //! and to read it when restoring the grid - typedef struct Block3d + struct Block3d { double collFactor; double deltaT; @@ -74,24 +74,24 @@ class MPIIORestartCoProcessor: public CoProcessor bool active; bool compressible; bool withForcing; - }Block3d; + }; //! \struct dataSet //! \brief Structure containes information identifying the block //! \details The structure is used to find the needed block in the grid when restoring a dataSet - typedef struct dataSet + struct DataSet { int x1; int x2; int x3; int level; - }dataSet; + }; //! \struct BoundaryCondition //! \brief Structure containes information about boundary conditions of the block //! \details The structure is used to write data describing boundary conditions of the blocks when saving the grid //! and to read it when restoring the grid - typedef struct BoundaryCondition + struct BoundaryCondition { long long noslipBoundaryFlags; // MPI_LONG_LONG long long slipBoundaryFlags; @@ -114,14 +114,14 @@ class MPIIORestartCoProcessor: public CoProcessor float q[26]; // MPI_FLOAT char algorithmType; - }BoundaryCondition; + }; //! \struct BCAdd //! \brief Structure containes information identifying the block //! and some parameters of the arrays of boundary conditions that are equal in all blocks //! \details The structure is used to find the needed block in the grid when restoring a dataSet //! and to set common parameters - typedef struct BCAdd + struct BCAdd { int x1; // to find the right block int x2; @@ -129,7 +129,7 @@ class MPIIORestartCoProcessor: public CoProcessor int level; int boundCond_count; // how many BoundaryCondition-structures are in this block int indexContainer_count; // how many indexContainer-values are in this block - }BCAdd; + }; public: MPIIORestartCoProcessor(Grid3DPtr grid, UbSchedulerPtr s, const std::string& path, CommunicatorPtr comm); @@ -154,10 +154,11 @@ public: protected: std::string path; CommunicatorPtr comm; + bool mpiTypeFreeFlag; private: MPI_Datatype gridParamType, blockParamType, block3dType, dataSetType, dataSetDoubleType, boundCondType, boundCondType1000, boundCondTypeAdd, bcindexmatrixType; - blockParam blockParamStr; + BlockParam blockParamStr; }; diff --git a/source/VirtualFluidsCore/CoProcessors/WriteMacroscopicQuantitiesCoProcessor.cpp b/source/VirtualFluidsCore/CoProcessors/WriteMacroscopicQuantitiesCoProcessor.cpp index 4ebc36eaf..2e65951c6 100644 --- a/source/VirtualFluidsCore/CoProcessors/WriteMacroscopicQuantitiesCoProcessor.cpp +++ b/source/VirtualFluidsCore/CoProcessors/WriteMacroscopicQuantitiesCoProcessor.cpp @@ -23,7 +23,7 @@ WriteMacroscopicQuantitiesCoProcessor::WriteMacroscopicQuantitiesCoProcessor(Gri conv(conv), comm(comm) { - gridRank = Communicator::getInstance()->getProcessID(); + gridRank = comm->getProcessID(); minInitLevel = this->grid->getCoarsestInitializedLevel(); maxInitLevel = this->grid->getFinestInitializedLevel(); diff --git a/source/VirtualFluidsCore/CoProcessors/WriteMacroscopicQuantitiesCoProcessor.h b/source/VirtualFluidsCore/CoProcessors/WriteMacroscopicQuantitiesCoProcessor.h index 257023266..a5613f491 100644 --- a/source/VirtualFluidsCore/CoProcessors/WriteMacroscopicQuantitiesCoProcessor.h +++ b/source/VirtualFluidsCore/CoProcessors/WriteMacroscopicQuantitiesCoProcessor.h @@ -43,19 +43,18 @@ private: int gridRank; CommunicatorPtr comm; - friend class boost::serialization::access; - template<class Archive> - void serialize(Archive & ar, const unsigned int version) - { - ar & boost::serialization::base_object<CoProcessor>(*this); - ar & path; - ar & conv; - ar & bcInformation; - ar & blockVector; - ar & minInitLevel; - ar & maxInitLevel; - ar & gridRank; - ar & writer; - } + //friend class boost::serialization::access; + //template<class Archive> + //void serialize(Archive & ar, const unsigned int version) + //{ + // ar & boost::serialization::base_object<CoProcessor>(*this); + // ar & path; + // ar & conv; + // ar & blockVector; + // ar & minInitLevel; + // ar & maxInitLevel; + // ar & gridRank; + // ar & writer; + //} }; #endif diff --git a/source/VirtualFluidsCore/Grid/CalculationManager.cpp b/source/VirtualFluidsCore/Grid/CalculationManager.cpp index e2ca4b83d..902d081ca 100644 --- a/source/VirtualFluidsCore/Grid/CalculationManager.cpp +++ b/source/VirtualFluidsCore/Grid/CalculationManager.cpp @@ -8,6 +8,7 @@ #include <boost/foreach.hpp> #include <Calculator.h> +#include <MPICalculator.h> #include <PrePostBcCalculator.h> #include <Communicator.h> @@ -56,18 +57,26 @@ void CalculationManager::calculate() { try { - boost::thread_group threads; boost::shared_ptr<CalculationManager> this_ = shared_from_this(); - boost::exception_ptr error; + if (calcType == CalculationManager::Hybrid || calcType == CalculationManager::PrePostBc) + { + boost::thread_group threads; + boost::exception_ptr error; + + for (int i = 1; i<calcThreads.size(); i++) + { + threads.create_thread(boost::bind(&Calculator::calculate, calcThreads[i], endTime, this_, boost::ref(error))); + } + + calcThreads[0]->calculate(endTime, this_, boost::ref(error)); - for (int i=1; i<calcThreads.size(); i++) + threads.join_all(); + } + else { - threads.create_thread( boost::bind( &Calculator::calculate, calcThreads[i], endTime, this_,boost::ref(error)) ); + boost::dynamic_pointer_cast<MPICalculator>(calcThreads[0])->calculate(endTime, this_); } - calcThreads[0]->calculate(endTime, this_, boost::ref(error)); - - threads.join_all(); //if( error ) //{ @@ -182,8 +191,10 @@ CalculatorPtr CalculationManager::createCalculator(Grid3DPtr grid, SynchronizerP { switch (calcType) { + case CalculationManager::Hybrid: + return CalculatorPtr(new Calculator(grid, sync, mainThread)); case CalculationManager::MPI: - return CalculatorPtr (new Calculator(grid, sync, mainThread)); + return CalculatorPtr (new MPICalculator(grid)); #if defined VF_FETOL case CalculationManager::FETOL: return CalculatorPtr (new FETOLCalculator(grid, sync, mainThread)); diff --git a/source/VirtualFluidsCore/Grid/CalculationManager.h b/source/VirtualFluidsCore/Grid/CalculationManager.h index 2cd8155b8..f94162c11 100644 --- a/source/VirtualFluidsCore/Grid/CalculationManager.h +++ b/source/VirtualFluidsCore/Grid/CalculationManager.h @@ -12,9 +12,9 @@ typedef boost::shared_ptr<CalculationManager> CalculationManagerPtr; class CalculationManager : public boost::enable_shared_from_this<CalculationManager> { public: - enum CalculatorType{MPI, FETOL, PrePostBc}; + enum CalculatorType{Hybrid, MPI, FETOL, PrePostBc}; public: - CalculationManager(Grid3DPtr grid, int numOfThreads, double endTime, UbSchedulerPtr visScheduler, CalculatorType calcType = CalculationManager::MPI); + CalculationManager(Grid3DPtr grid, int numOfThreads, double endTime, UbSchedulerPtr visScheduler, CalculatorType calcType = CalculationManager::Hybrid); CalculationManager(Grid3DPtr grid, int numOfThreads, double endTime, UbSchedulerPtr visScheduler, CommunicatorPtr comm, int endDir, LBMReal nu, CalculatorType calcType = CalculationManager::MPI); virtual ~CalculationManager(); diff --git a/source/VirtualFluidsCore/Grid/Calculator.h b/source/VirtualFluidsCore/Grid/Calculator.h index 54c5373ec..fbb9e20e8 100644 --- a/source/VirtualFluidsCore/Grid/Calculator.h +++ b/source/VirtualFluidsCore/Grid/Calculator.h @@ -64,7 +64,7 @@ protected: int calcStep; std::vector< std::vector<Block3DPtr> > blocks; -private: + std::vector< std::vector< Block3DConnectorPtr > > localInterfaceBlockConns; std::vector< std::vector< Block3DConnectorPtr > > remoteInterfaceBlockConns; diff --git a/source/VirtualFluidsCore/Grid/MPICalculator.cpp b/source/VirtualFluidsCore/Grid/MPICalculator.cpp new file mode 100644 index 000000000..3fe8e4ca3 --- /dev/null +++ b/source/VirtualFluidsCore/Grid/MPICalculator.cpp @@ -0,0 +1,359 @@ +#include "MPICalculator.h" +#include <basics/utilities/UbException.h> +#include <boost/foreach.hpp> +#include "MathUtil.hpp" +#include "basics/writer/WbWriterVtkXmlASCII.h" + +//#define TIMING +//#define PRECOLLISIONBC + +MPICalculator::MPICalculator() +{ + +} +////////////////////////////////////////////////////////////////////////// +MPICalculator::MPICalculator(Grid3DPtr grid) : +Calculator(grid, SynchronizerPtr(), true) +{ + +} +////////////////////////////////////////////////////////////////////////// +void MPICalculator::calculate(const double& endTime, CalculationManagerPtr cm) +{ + UBLOG(logDEBUG1, "MPICalculator::calculate() - started"); + try + { + initConnectors(); + + int anzLevel = maxLevel-minLevel+1; + + int minInitLevel = minLevel; + int maxInitLevel = maxLevel-minLevel; + int straightStartLevel = minInitLevel; + int internalIterations = 1<<(maxInitLevel-minInitLevel); + int forwardStartLevel; + int threshold; + int startStep = int(grid->getTimeStep())+1; + + //UBLOG(logINFO, "startStep="<<startStep); + int anzCalcSteps = static_cast<int>(endTime); +#ifdef TIMING + UbTimer timer; + double time[6]; +#endif + + ////////////////////////////////////////////////////////////////////////// + // UBLOG(logINFO, "Number of connectors = " <<this->localConns[0].size()); + ////////////////////////////////////////////////////////////////////////// + + for (calcStep = startStep; calcStep<=anzCalcSteps+1; calcStep++) + { + + //exchange data between blocks for visualization + //sync->wait(); + ////if(visScheduler->isDue((double)(calcStep-1))) + ////{ + // //exchangeBlockData(minInitLevel, maxInitLevel, true); + ////} + + ////wait for write dump files + //write dump + grid->coProcess((double)(calcStep-1)); + + ////////////////////////////////////////////////////////////////////////// +#ifdef TIMING + UBLOG(logINFO, "calcStep = "<<calcStep); +#endif + ////////////////////////////////////////////////////////////////////////// + + for (int staggeredStep = 1; staggeredStep<=internalIterations; staggeredStep++) + { + forwardStartLevel = straightStartLevel; + if (staggeredStep==internalIterations) straightStartLevel = minInitLevel; + else + { + for (straightStartLevel = maxInitLevel, threshold = 1; + (staggeredStep&threshold)!=threshold; straightStartLevel--, threshold <<= 1); + } + ////////////////////////////////////////////////////////////////////////// +#ifdef TIMING + timer.resetAndStart(); +#endif + ////////////////////////////////////////////////////////////////////////// + + //applyPreCollisionBC(straightStartLevel, maxInitLevel); + + + + calculateBlocks(straightStartLevel, maxInitLevel); + ////calculateBlocks(minInitLevel, maxInitLevel, staggeredStep); +////////////////////////////////////////////////////////////////////////// +#ifdef TIMING + time[0] = timer.stop(); + UBLOG(logINFO, "calculateBlocks time = "<<time[0]); +#endif + ////////////////////////////////////////////////////////////////////////// + + //exchange data between blocks + //Sleep(10000); + exchangeBlockData(straightStartLevel, maxInitLevel); + ////////////////////////////////////////////////////////////////////////// +#ifdef TIMING + time[1] = timer.stop(); + UBLOG(logINFO, "exchangeBlockData time = "<<time[1]); +#endif + ////////////////////////////////////////////////////////////////////////// + //applyBCs(straightStartLevel, maxInitLevel); + applyPostCollisionBC(straightStartLevel, maxInitLevel); + + ////////////////////////////////////////////////////////////////////////// +#ifdef TIMING + time[2] = timer.stop(); + UBLOG(logINFO, "applyBCs time = "<<time[2]); +#endif + ////////////////////////////////////////////////////////////////////////// + + //swap distributions in kernel + swapDistributions(straightStartLevel, maxInitLevel); + +#ifdef PRECOLLISIONBC + exchangeBlockData(straightStartLevel, maxInitLevel); + applyPreCollisionBC(straightStartLevel, maxInitLevel); +#endif + + + ////////////////////////////////////////////////////////////////////////// +#ifdef TIMING + time[3] = timer.stop(); + UBLOG(logINFO, "swapDistributions time = "<<time[3]); +#endif + ////////////////////////////////////////////////////////////////////////// + + if (refinement) + { + // //exchange data between blocks for grid refinement + ////exchangeInterfaceBlockData(straightStartLevel, maxInitLevel, true); + //DOES NOT NEED + if (straightStartLevel<maxInitLevel) + exchangeBlockData(straightStartLevel, maxInitLevel); + // //exchangeInterfaceBlockData(straightStartLevel, maxInitLevel, true); + ////////////////////////////////////////////////////////////////////////// +#ifdef TIMING + time[4] = timer.stop(); + UBLOG(logINFO, "refinement exchangeBlockData time = "<<time[4]); +#endif + ////////////////////////////////////////////////////////////////////////// + //now ghost nodes have actual values + //interpolation of interface nodes between grid levels + interpolation(straightStartLevel, maxInitLevel); + ////////////////////////////////////////////////////////////////////////// +#ifdef TIMING + time[5] = timer.stop(); + UBLOG(logINFO, "refinement interpolation time = "<<time[5]); +#endif + ////////////////////////////////////////////////////////////////////////// + } + + if (taValuesCoProcessor) + { + taValuesCoProcessor->calculateSubtotal(calcStep-1); + } + + + } + //exchange data between blocks for visualization + if (mainThread) visScheduler->isDue((double)(calcStep-1)); + if ((int)visScheduler->getNextDueTime()==calcStep) + { + exchangeBlockData(straightStartLevel, maxInitLevel); + } + //now ghost nodes have actual values + + //dynamic load balancing + //sync->wait(); + //if (mainThread && !loadBalancingComp) + //{ + // loadBalancingComp = cm->balance(); + //} + + } + UBLOG(logDEBUG1, "MPICalculator::calculate() - stoped"); + } + catch (std::exception& e) + { + //error = boost::current_exception(); + UBLOG(logERROR, e.what()); + UBLOG(logERROR, " step = "<<calcStep); + } +} +////////////////////////////////////////////////////////////////////////// +void MPICalculator::calculateBlocks(int startLevel, int maxInitLevel) +{ + Block3DPtr blockTemp; + try + { + //startLevel bis maxInitLevel + for (int level = startLevel; level<=maxInitLevel; level++) + { + //timer.resetAndStart(); + //call LBM kernel + BOOST_FOREACH(Block3DPtr block, blocks[level]) + { + blockTemp = block; + block->getKernel()->calculate(); + } + //timer.stop(); + //UBLOG(logINFO, "level = " << level << " blocks = " << blocks[level].size() << " collision time = " << timer.getTotalTime()); + } + } + catch (std::exception& e) + { + //error = boost::current_exception(); + UBLOG(logERROR, e.what()); + UBLOG(logERROR, blockTemp->toString()<<" step = "<<calcStep); + exit(EXIT_FAILURE); + } +} +////////////////////////////////////////////////////////////////////////// +void MPICalculator::calculateBlocks(int minInitLevel, int maxInitLevel, int staggeredStep) +{ + int p, maxi, maxir, maxidp, start, end; + for (int level = minInitLevel; level<=maxInitLevel; level++) + { + p = 1<<(maxInitLevel-level); + maxi = maxir = static_cast<int>(blocks[level].size()); + maxidp = maxi/p; + if (p>maxi && maxi!=0) { + maxidp = 1; + maxi = p; + } + start = (staggeredStep-1)*maxidp; + if (start>=maxi) + start = 0; + end = start+maxidp; + if ((end+p)>=maxi) + end = maxi; + for (int i = start; i<end; i++) + { + if (i<maxir) + blocks[level][i]->getKernel()->calculate(); + } + } +} +////////////////////////////////////////////////////////////////////////// + +////////////////////////////////////////////////////////////////////////// +void MPICalculator::exchangeBlockData(int startLevel, int maxInitLevel) +{ + //startLevel bis maxInitLevel + for (int level = startLevel; level<=maxInitLevel; level++) + { + connectorsPrepare(localConns[level]); + connectorsPrepare(remoteConns[level]); + + connectorsSend(localConns[level]); + connectorsSend(remoteConns[level]); + + connectorsReceive(localConns[level]); + connectorsReceive(remoteConns[level]); + } +} +////////////////////////////////////////////////////////////////////////// +void MPICalculator::exchangeInterfaceBlockData(int startLevel, int maxInitLevel) +{ + //startLevel bis maxInitLevel + for (int level = startLevel; level<=maxInitLevel; level++) + { + connectorsPrepare(localInterfaceBlockConns[level]); + connectorsPrepare(remoteInterfaceBlockConns[level]); + + connectorsSend(localInterfaceBlockConns[level]); + connectorsSend(remoteInterfaceBlockConns[level]); + + connectorsReceive(localInterfaceBlockConns[level]); + connectorsReceive(remoteInterfaceBlockConns[level]); + } +} +////////////////////////////////////////////////////////////////////////// +void MPICalculator::swapDistributions(int startLevel, int maxInitLevel) +{ + //startLevel bis maxInitLevel + for (int level = startLevel; level<=maxInitLevel; level++) + { + BOOST_FOREACH(Block3DPtr block, blocks[level]) + { + block->getKernel()->swapDistributions(); + } + } +} +////////////////////////////////////////////////////////////////////////// +void MPICalculator::connectorsPrepare(std::vector< Block3DConnectorPtr >& connectors) +{ + BOOST_FOREACH(Block3DConnectorPtr c, connectors) + { + c->prepareForReceive(); + c->prepareForSend(); + } +} +////////////////////////////////////////////////////////////////////////// +void MPICalculator::connectorsSend(std::vector< Block3DConnectorPtr >& connectors) +{ + BOOST_FOREACH(Block3DConnectorPtr c, connectors) + { + c->fillSendVectors(); + c->sendVectors(); + } +} +////////////////////////////////////////////////////////////////////////// +void MPICalculator::connectorsReceive(std::vector< Block3DConnectorPtr >& connectors) +{ + BOOST_FOREACH(Block3DConnectorPtr c, connectors) + { + c->receiveVectors(); + c->distributeReceiveVectors(); + } +} +////////////////////////////////////////////////////////////////////////// +void MPICalculator::interpolation(int startLevel, int maxInitLevel) +{ + + + for (int level = startLevel; level<maxInitLevel; level++) + { + connectorsPrepare(localInterConns[level]); + connectorsPrepare(remoteInterConns[level]); + } + + + + for (int level = startLevel; level<maxInitLevel; level++) + { + connectorsSend(localInterConns[level]); + connectorsSend(remoteInterConns[level]); + } + + + + for (int level = startLevel; level<maxInitLevel; level++) + { + connectorsReceive(localInterConns[level]); + connectorsReceive(remoteInterConns[level]); + } + + + +} +////////////////////////////////////////////////////////////////////////// +////////////////////////////////////////////////////////////////////////// +void MPICalculator::applyPostCollisionBC(int startLevel, int maxInitLevel) +{ + //startLevel bis maxInitLevel + for (int level = startLevel; level<=maxInitLevel; level++) + { + BOOST_FOREACH(Block3DPtr block, blocks[level]) + { + block->getKernel()->getBCProcessor()->applyPostCollisionBC(); + } + } +} + diff --git a/source/VirtualFluidsCore/Grid/MPICalculator.h b/source/VirtualFluidsCore/Grid/MPICalculator.h new file mode 100644 index 000000000..48727b585 --- /dev/null +++ b/source/VirtualFluidsCore/Grid/MPICalculator.h @@ -0,0 +1,45 @@ +#ifndef MPICALCULATOR_H +#define MPICALCULATOR_H + +#include "Grid3D.h" +#include "Block3D.h" +#include "Synchronizer.h" +#include "MathUtil.hpp" +#include "basics/utilities/UbScheduler.h" +#include "basics/utilities/UbTiming.h" +#include "LoadBalancer.h" +#include "TimeAveragedValuesCoProcessor.h" + + +class MPICalculator; +typedef boost::shared_ptr<MPICalculator> MPICalculatorPtr; + +#include "CalculationManager.h" + +class MPICalculator : public Calculator +{ +public: + MPICalculator(); + MPICalculator(Grid3DPtr grid); + virtual ~MPICalculator() {} + virtual void calculate(const double& endTime, CalculationManagerPtr cm); + +protected: + void calculateBlocks(int startLevel, int maxInitLevel); + void calculateBlocks(int minInitLevel, int maxInitLevel, int staggeredStep); + void swapDistributions(int startLevel, int maxInitLevel); + virtual void exchangeBlockData(int startLevel, int maxInitLevel); + void exchangeInterfaceBlockData(int startLevel, int maxInitLevel); + virtual void connectorsPrepare(std::vector< Block3DConnectorPtr >& connectors); + virtual void connectorsSend(std::vector< Block3DConnectorPtr >& connectors); + virtual void connectorsReceive(std::vector< Block3DConnectorPtr >& connectors); + void interpolation(int startLevel, int maxInitLevel); + void deleteConnectors(std::vector< std::vector< Block3DConnectorPtr > >& conns); + void applyPostCollisionBC(int startLevel, int maxInitLevel); +private: + + +}; + +#endif + diff --git a/source/VirtualFluidsCore/LBM/IncompressibleCumulantWithSpongeLayerLBMKernel.cpp b/source/VirtualFluidsCore/LBM/IncompressibleCumulantWithSpongeLayerLBMKernel.cpp index 84ca58f08..f181acea3 100644 --- a/source/VirtualFluidsCore/LBM/IncompressibleCumulantWithSpongeLayerLBMKernel.cpp +++ b/source/VirtualFluidsCore/LBM/IncompressibleCumulantWithSpongeLayerLBMKernel.cpp @@ -28,6 +28,96 @@ void IncompressibleCumulantWithSpongeLayerLBMKernel::init() DistributionArray3DPtr d(new D3Q27EsoTwist3DSplittedVector(nx1+2, nx2+2, nx3+2, -999.0)); dataSet->setFdistributions(d); } +////////////////////////////////////////////////////////////////////////// +void IncompressibleCumulantWithSpongeLayerLBMKernel::setRelaxFactorParam(int vdir, double vL1, double vdx, double vSP) +{ + direction = vdir; + L1 = vL1; + dx = vdx; + SP = vSP; +} +////////////////////////////////////////////////////////////////////////// +void IncompressibleCumulantWithSpongeLayerLBMKernel::initRelaxFactor(int vdir, double vL1, double vdx, double vSP) +{ + direction = vdir; + L1 = vL1; + dx = vdx; + SP = vSP; + + double sizeX = L1 / dx; + double sizeSP = SP / dx; + double muX1, muX2, muX3; + + LBMReal spongeFactor; + + BCArray3D& bcArray = this->getBCProcessor()->getBCArray(); + + const int bcArrayMaxX1 = (int)bcArray.getNX1(); + const int bcArrayMaxX2 = (int)bcArray.getNX2(); + const int bcArrayMaxX3 = (int)bcArray.getNX3(); + + int minX1 = 0; + int minX2 = 0; + int minX3 = 0; + int maxX1 = bcArrayMaxX1 - ghostLayerWidth - 1; + int maxX2 = bcArrayMaxX2 - ghostLayerWidth - 1; + int maxX3 = bcArrayMaxX3 - ghostLayerWidth - 1; + + RelaxationFactorArray3DPtr relaxationFactorPtr = CbArray3D<LBMReal, IndexerX3X2X1>::CbArray3DPtr(new CbArray3D<LBMReal, IndexerX3X2X1>(maxX1, maxX2, maxX3)); + dataSet->setRelaxationFactor(relaxationFactorPtr); + + for (int x3 = minX3; x3 < maxX3; x3++) + { + for (int x2 = minX2; x2 < maxX2; x2++) + { + for (int x1 = minX1; x1 < maxX1; x1++) + { + switch (direction) + { + case D3Q27System::E: + muX1 = (double)(x1 + ix1 * maxX1); + if (muX1 >= (sizeX - sizeSP) / deltaT) + spongeFactor = (sizeX - (muX1 * deltaT + 1)) / sizeSP / 2.0 + 0.5; + else spongeFactor = 1.0; + break; + case D3Q27System::W: + muX1 = (double)(x1 + ix1 * maxX1); + if (muX1 <= sizeSP / deltaT) + spongeFactor = (sizeSP - (muX1 * deltaT + 1)) / sizeSP / 2.0 + 0.5; + else spongeFactor = 1.0; + break; + case D3Q27System::N: + muX2 = (double)(x2 + ix2 * maxX2); + if (muX2 >= (sizeX - sizeSP) / deltaT) + spongeFactor = (sizeX - (muX2 * deltaT + 1)) / sizeSP / 2.0 + 0.5; + else spongeFactor = 1.0; + break; + case D3Q27System::S: + muX2 = (double)(x2 + ix2 * maxX2); + if (muX2 <= sizeSP / deltaT) + spongeFactor = (sizeSP - (muX2 * deltaT + 1)) / sizeSP / 2.0 + 0.5; + else spongeFactor = 1.0; + break; + case D3Q27System::T: + muX3 = (double)(x3 + ix3 * maxX3); + if (muX3 >= (sizeX - sizeSP) / deltaT) + spongeFactor = (sizeX - (muX3 * deltaT + 1)) / sizeSP / 2.0 + 0.5; + else spongeFactor = 1.0; + break; + case D3Q27System::B: + muX3 = (double)(x3 + ix3 * maxX3); + if (muX3 <= sizeSP / deltaT) + spongeFactor = (sizeSP - (muX3 * deltaT + 1)) / sizeSP / 2.0 + 0.5; + else spongeFactor = 1.0; + break; + default: throw UbException(UB_EXARGS, "unknown dir"); + } + (*relaxationFactorPtr)(x1, x2, x3) = spongeFactor; + } + } + } +} + ////////////////////////////////////////////////////////////////////////// LBMKernelPtr IncompressibleCumulantWithSpongeLayerLBMKernel::clone() { @@ -53,8 +143,9 @@ LBMKernelPtr IncompressibleCumulantWithSpongeLayerLBMKernel::clone() kernel->setWithSpongeLayer(withSpongeLayer); if(withSpongeLayer) kernel->setSpongeLayer(muSpongeLayer); + boost::dynamic_pointer_cast<IncompressibleCumulantWithSpongeLayerLBMKernel>(kernel)->initRelaxFactor(direction, L1, dx, SP); return kernel; -} +} ////////////////////////////////////////////////////////////////////////// void IncompressibleCumulantWithSpongeLayerLBMKernel::calculate() { @@ -98,9 +189,9 @@ void IncompressibleCumulantWithSpongeLayerLBMKernel::collideAll() //initialization of sponge layer variables //if (withSpongeLayer) //{ - muDeltaT = deltaT; - muSpongeLayer.DefineVar("dx",&muDeltaT); - muSpongeLayer.DefineVar("x1",&muX1); muSpongeLayer.DefineVar("x2",&muX2); muSpongeLayer.DefineVar("x3",&muX3); + //muDeltaT = deltaT; + //muSpongeLayer.DefineVar("dt",&muDeltaT); + //muSpongeLayer.DefineVar("x1",&muX1); muSpongeLayer.DefineVar("x2",&muX2); muSpongeLayer.DefineVar("x3",&muX3); //} ///////////////////////////////////// @@ -109,6 +200,7 @@ void IncompressibleCumulantWithSpongeLayerLBMKernel::collideAll() zeroDistributions = boost::dynamic_pointer_cast<D3Q27EsoTwist3DSplittedVector>(dataSet->getFdistributions())->getZeroDistributions(); BCArray3D& bcArray = this->getBCProcessor()->getBCArray(); + RelaxationFactorArray3DPtr relaxationFactorPtr = dataSet->getRelaxationFactor(); const int bcArrayMaxX1 = (int)bcArray.getNX1(); const int bcArrayMaxX2 = (int)bcArray.getNX2(); @@ -122,6 +214,7 @@ void IncompressibleCumulantWithSpongeLayerLBMKernel::collideAll() int maxX3 = bcArrayMaxX3-ghostLayerWidth-1; LBMReal collFactor0 = collFactor; + LBMReal spongeFactor; for(int x3 = minX3; x3 <= maxX3; x3++) { @@ -216,12 +309,12 @@ void IncompressibleCumulantWithSpongeLayerLBMKernel::collideAll() //{ //if (!withForcing) //{ - muX1 = (double)(x1-1+ix1*maxX1); + //lk// muX1 = (double)(x1-1+ix1*maxX1); //muX2 = (double)(x2-1+ix2*maxX2); //muX3 = (double)(x3-1+ix3*maxX3); //} //spongeFactor ist von Funktion in muSpongeLayer abhängich und variiert zwischen 1 (nix tun) und 0.5 (collFactor etwa auf 1); - LBMReal spongeFactor = muSpongeLayer.Eval(); + //lk //LBMReal spongeFactor = muSpongeLayer.Eval(); //if (spongeFactor == 0.5) //{ @@ -231,8 +324,8 @@ void IncompressibleCumulantWithSpongeLayerLBMKernel::collideAll() //if(muX3 == ix3*maxX3/2 && muX2==ix2*maxX2/2) // UBLOG(logINFO," x1="<<muX1<<" spongeFactor = " << spongeFactor <<" collFactor="<<collFactor); + spongeFactor = (*relaxationFactorPtr)(x1-1, x2-1, x3-1); collFactor *= spongeFactor; - //if(muX3 == ix3*maxX3/2 && muX2==ix2*maxX2/2) // UBLOG(logINFO," x1="<<muX1<<" spongeFactor = " << spongeFactor <<" collFactor="<<collFactor); @@ -246,7 +339,7 @@ void IncompressibleCumulantWithSpongeLayerLBMKernel::collideAll() m1=mfacc+mfcaa; m2=mfaac+mfcca; oMdrho+=m0; - m1+=m2; + m1+=m2; oMdrho+=m1; m0=mfbac+mfbca; m1=mfbaa+mfbcc; diff --git a/source/VirtualFluidsCore/LBM/IncompressibleCumulantWithSpongeLayerLBMKernel.h b/source/VirtualFluidsCore/LBM/IncompressibleCumulantWithSpongeLayerLBMKernel.h index a1bb9b14f..0931848fc 100644 --- a/source/VirtualFluidsCore/LBM/IncompressibleCumulantWithSpongeLayerLBMKernel.h +++ b/source/VirtualFluidsCore/LBM/IncompressibleCumulantWithSpongeLayerLBMKernel.h @@ -25,7 +25,7 @@ typedef boost::shared_ptr<IncompressibleCumulantWithSpongeLayerLBMKernel> LBMKer //! kernel->setWithSpongeLayer(true); //! kernel->setSpongeLayer(spongeLayer); //! \endcode -//! \author K. Kucher, M. Geier +//! \author K. Kucher, M. Geier, A. Karanchuk class IncompressibleCumulantWithSpongeLayerLBMKernel : public IncompressibleCumulantLBMKernel { public: @@ -39,9 +39,19 @@ public: virtual ~IncompressibleCumulantWithSpongeLayerLBMKernel(void); LBMKernelPtr clone(); void calculate(); + void initRelaxFactor(int vdir, double vL1, double vdx, double vSP); + //! \param vdir where the sponge layer is placed + //! \param vL1 length of simulation domain + //! \param vdx subgrid space + //! \param vSP length of sponge layer + void setRelaxFactorParam(int vdir, double vL1, double vdx, double vSP); protected: void init(); - LBMReal OxyyMxzz; + LBMReal OxyyMxzz; + int direction; + double L1; + double dx; + double SP; friend class boost::serialization::access; template<class Archive> diff --git a/source/VirtualFluidsCore/Parallel/MPICommunicator.cpp b/source/VirtualFluidsCore/Parallel/MPICommunicator.cpp index dfb21b5b3..a4fcebd90 100644 --- a/source/VirtualFluidsCore/Parallel/MPICommunicator.cpp +++ b/source/VirtualFluidsCore/Parallel/MPICommunicator.cpp @@ -31,6 +31,7 @@ MPICommunicator::~MPICommunicator() if (!_mpiFinalized) { MPI_Finalize(); + //UBLOG(logINFO, "MPI_Finalize()"); } } ////////////////////////////////////////////////////////////////////////// -- GitLab