From 9607bedb3c90f8c2af75f6775f606c1898874f6f Mon Sep 17 00:00:00 2001 From: Konstantin Kutscher <kutscher@irmb.tu-bs.de> Date: Thu, 14 Apr 2016 15:04:05 +0000 Subject: [PATCH] Fixed issue with boudary condition for thin wall --- source/Applications.cmake | 2 +- source/Applications/LaminarTubeFlow/ltf.cpp | 10 +- .../pChannel/configBombadilpChannel.cfg | 4 +- source/Applications/pChannel/pChannel.cpp | 4 +- source/Applications/sphere/CMakeLists.txt | 4 +- source/Applications/sphere/config.txt | 8 +- source/Applications/sphere/sphere.cpp | 340 ++++------ source/VirtualFluids.h | 4 +- .../BoundaryCondition/BCArray.cpp | 196 ++++++ .../BoundaryCondition/BCArray.h | 346 ++++------ .../BoundaryCondition/BCProcessor.h | 1 - .../BoundaryCondition/BoundaryCondition.cpp | 24 +- .../BoundaryCondition/BoundaryCondition.h | 18 +- .../BoundaryConditionProcessor.cpp | 58 -- .../BoundaryConditionProcessor.h | 26 - .../BoundaryCondition/D3Q27ETBCProcessor.cpp | 535 +--------------- .../BoundaryCondition/D3Q27ETBCProcessor.h | 46 +- .../D3Q27ETForThinWallBCProcessor.cpp | 57 +- .../D3Q27ETForThinWallBCProcessor.h | 1 - .../EqDensityBoundaryCondition.cpp | 5 + .../EqDensityBoundaryCondition.h | 1 + .../HighViscosityNoSlipBoundaryCondition.cpp | 5 + .../HighViscosityNoSlipBoundaryCondition.h | 1 + .../NoSlipBoundaryCondition.cpp | 5 + .../NoSlipBoundaryCondition.h | 1 + .../NonEqDensityBoundaryCondition.cpp | 5 + .../NonEqDensityBoundaryCondition.h | 1 + .../NonReflectingDensityBoundaryCondition.cpp | 5 + .../NonReflectingDensityBoundaryCondition.h | 1 + ...NonReflectingVelocityBoundaryCondition.cpp | 5 + .../NonReflectingVelocityBoundaryCondition.h | 1 + .../SlipBoundaryCondition.cpp | 5 + .../BoundaryCondition/SlipBoundaryCondition.h | 1 + .../ThinWallNoSlipBoundaryCondition.cpp | 31 +- .../ThinWallNoSlipBoundaryCondition.h | 6 +- .../VelocityBoundaryCondition.cpp | 5 + .../VelocityBoundaryCondition.h | 1 + .../MacroscopicQuantitiesCoProcessor.cpp | 12 +- .../Connectors/Block3DConnectorFactory.cpp | 4 +- .../BoostSerializationClassExportHelper.h | 2 - .../Grid/CalculationManager.cpp | 15 +- .../Grid/CalculationManager.h | 2 - source/VirtualFluidsCore/Grid/Calculator.cpp | 47 +- source/VirtualFluidsCore/Grid/Calculator.h | 5 +- source/VirtualFluidsCore/Grid/Calculator2.cpp | 546 ---------------- source/VirtualFluidsCore/Grid/Calculator2.h | 78 --- source/VirtualFluidsCore/Grid/FastSignal.hpp | 46 -- .../VirtualFluidsCore/Grid/MTCalculator.cpp | 589 ------------------ source/VirtualFluidsCore/Grid/MTCalculator.h | 82 --- .../Grid/PrePostBcCalculator.cpp | 5 - .../Grid/PrePostBcCalculator.h | 2 - .../Utilities/ConfigurationFile.hpp | 50 +- .../BoundaryConditionBlockVisitor.cpp | 4 - .../Visitors/BoundaryConditionBlockVisitor.h | 1 - 54 files changed, 663 insertions(+), 2596 deletions(-) create mode 100644 source/VirtualFluidsCore/BoundaryCondition/BCArray.cpp delete mode 100644 source/VirtualFluidsCore/BoundaryCondition/BoundaryConditionProcessor.cpp delete mode 100644 source/VirtualFluidsCore/BoundaryCondition/BoundaryConditionProcessor.h delete mode 100644 source/VirtualFluidsCore/Grid/Calculator2.cpp delete mode 100644 source/VirtualFluidsCore/Grid/Calculator2.h delete mode 100644 source/VirtualFluidsCore/Grid/FastSignal.hpp delete mode 100644 source/VirtualFluidsCore/Grid/MTCalculator.cpp delete mode 100644 source/VirtualFluidsCore/Grid/MTCalculator.h diff --git a/source/Applications.cmake b/source/Applications.cmake index 9e8705197..891ca0b1a 100644 --- a/source/Applications.cmake +++ b/source/Applications.cmake @@ -2,7 +2,7 @@ #add_subdirectory(Applications/gridRf) #add_subdirectory(Applications/greenvortex) # add_subdirectory(Applications/micropart) -# add_subdirectory(Applications/sphere) +add_subdirectory(Applications/sphere) #add_subdirectory(Applications/vfscript) #add_subdirectory(Applications/reefer) #add_subdirectory(Applications/bananas) diff --git a/source/Applications/LaminarTubeFlow/ltf.cpp b/source/Applications/LaminarTubeFlow/ltf.cpp index a3cf58ace..b3c3f842c 100644 --- a/source/Applications/LaminarTubeFlow/ltf.cpp +++ b/source/Applications/LaminarTubeFlow/ltf.cpp @@ -124,7 +124,9 @@ void run(string configname) //set connectors D3Q27InterpolationProcessorPtr iProcessor(new D3Q27IncompressibleOffsetInterpolationProcessor()); - D3Q27SetConnectorsBlockVisitor setConnsVisitor(comm, true, D3Q27System::ENDDIR, nuLB, iProcessor); + //D3Q27SetConnectorsBlockVisitor setConnsVisitor(comm, true, D3Q27System::ENDDIR, nuLB, iProcessor); + ConnectorFactoryPtr factory(new Block3DConnectorFactory()); + ConnectorBlockVisitor setConnsVisitor(comm, nuLB, iProcessor,factory); grid->accept(setConnsVisitor); //domain decomposition for threads @@ -164,11 +166,13 @@ void run(string configname) kernel = LBMKernel3DPtr(new LBMKernelETD3Q27CCLB(blocknx[0], blocknx[1], blocknx[2], LBMKernelETD3Q27CCLB::NORMAL)); // - BCProcessorPtr bcProc(new D3Q27ETBCProcessor()); + //BCProcessorPtr bcProc(new D3Q27ETBCProcessor()); + BCProcessorPtr bcProc(new D3Q27ETForThinWallBCProcessor()); BoundaryConditionPtr noSlipBC; BoundaryConditionPtr velBC; BoundaryConditionPtr denBC; - noSlipBC = BoundaryConditionPtr(new NoSlipBoundaryCondition()); + //noSlipBC = BoundaryConditionPtr(new NoSlipBoundaryCondition()); + noSlipBC = BoundaryConditionPtr(new ThinWallNoSlipBoundaryCondition()); velBC = BoundaryConditionPtr(new VelocityBoundaryCondition()); denBC = BoundaryConditionPtr(new NonEqDensityBoundaryCondition()); bcProc->addBC(noSlipBC); diff --git a/source/Applications/pChannel/configBombadilpChannel.cfg b/source/Applications/pChannel/configBombadilpChannel.cfg index 20c3befb7..31bb2040c 100644 --- a/source/Applications/pChannel/configBombadilpChannel.cfg +++ b/source/Applications/pChannel/configBombadilpChannel.cfg @@ -4,7 +4,7 @@ pathname = d:/temp/pChannel pathGeo = d:/Projects/SFB880/GeometrienPoroeseMedien/PA80-110 -numOfThreads = 1 +numOfThreads = 4 availMem = 4e9 logToFile = false @@ -91,6 +91,6 @@ timeAvStart = 1 timeAvStop = 3 endTime = 200000 -outTime = 1000 +outTime = 100 diff --git a/source/Applications/pChannel/pChannel.cpp b/source/Applications/pChannel/pChannel.cpp index 63890426b..2e202ca04 100644 --- a/source/Applications/pChannel/pChannel.cpp +++ b/source/Applications/pChannel/pChannel.cpp @@ -90,8 +90,6 @@ void run(string configname) //real velocity is 49.63 m/s double nu_LB; - BoundaryConditionProcessorPtr bcProcessor(new BoundaryConditionProcessor()); - Grid3DPtr grid(new Grid3D(comm)); ////////////////////////////////////////////////////////////////////////// @@ -388,6 +386,8 @@ void run(string configname) ////set connectors D3Q27InterpolationProcessorPtr iProcessor(new D3Q27IncompressibleOffsetInterpolationProcessor()); D3Q27SetConnectorsBlockVisitor setConnsVisitor(comm, true, D3Q27System::ENDDIR, nu_LB, iProcessor); + //ConnectorFactoryPtr factory(new Block3DConnectorFactory()); + //ConnectorBlockVisitor setConnsVisitor(comm, nu_LB, iProcessor, factory); grid->accept(setConnsVisitor); //domain decomposition for threads diff --git a/source/Applications/sphere/CMakeLists.txt b/source/Applications/sphere/CMakeLists.txt index 1d17bc866..a2ec13053 100644 --- a/source/Applications/sphere/CMakeLists.txt +++ b/source/Applications/sphere/CMakeLists.txt @@ -5,7 +5,7 @@ CMAKE_MINIMUM_REQUIRED(VERSION 2.8) ######################################################## PROJECT(sphere) -INCLUDE(${SOURCE_ROOT}/lib/IncludsList.txt) +INCLUDE(${SOURCE_ROOT}/IncludsList.cmake) ################################################################# ### LOCAL FILES ### @@ -17,7 +17,7 @@ FILE(GLOB SPECIFIC_FILES ${CMAKE_CURRENT_SOURCE_DIR}/*.h SET(ALL_SOURCES ${ALL_SOURCES} ${SPECIFIC_FILES}) SOURCE_GROUP(src FILES ${SPECIFIC_FILES}) -SET(CAB_ADDITIONAL_LINK_LIBRARIES vfluids) +SET(CAB_ADDITIONAL_LINK_LIBRARIES VirtualFluids) ################################################################# ### CREATE PROJECT ### diff --git a/source/Applications/sphere/config.txt b/source/Applications/sphere/config.txt index f6440c075..c1e511c05 100644 --- a/source/Applications/sphere/config.txt +++ b/source/Applications/sphere/config.txt @@ -1,5 +1,5 @@ #Ordner für Simulationsergebnisse -path=d:/temp/insitu_demo +path=d:/temp/sphere #Verfügbare Arbeitsspeicher in Byte memory=21e9 @@ -14,7 +14,9 @@ outstep=1 endstep=100000 #Anzahl von Threads -threads=4 +threads=1 #max refierment level (1 - 5) -level=1 \ No newline at end of file +level=1 + +test = true \ No newline at end of file diff --git a/source/Applications/sphere/sphere.cpp b/source/Applications/sphere/sphere.cpp index 0b9a3746d..4d03c01f2 100644 --- a/source/Applications/sphere/sphere.cpp +++ b/source/Applications/sphere/sphere.cpp @@ -1,11 +1,11 @@ -#include <vfluids.h> +#include <VirtualFluids.h> #include <set> #include <map> using namespace std; //////////////////////////////////////////////////////////////////////// -void chanel(const char *cstr) +void run(string configname) { try { @@ -13,8 +13,6 @@ void chanel(const char *cstr) //Sleep(30000); string machine = QUOTEME(CAB_MACHINE); - string pathname; - double availMem = 0; CommunicatorPtr comm = MPICommunicator::getInstance(); @@ -22,19 +20,33 @@ void chanel(const char *cstr) int mybundle = comm->getBundleID(); int root = comm->getRoot(); - ConfigFileReader cf(cstr); - if ( !cf.read() ) - { - std::string exceptionText = "Unable to read configuration file\n"; - throw exceptionText; - } + //ConfigFileReader cf(cstr); + //if ( !cf.read() ) + //{ + // std::string exceptionText = "Unable to read configuration file\n"; + // throw exceptionText; + //} + + //pathname = cf.getValue("path"); + //availMem = UbSystem::stringTo<double>(cf.getValue("memory")); + //string metafile = cf.getValue("metafile"); + //double outstep = UbSystem::stringTo<double>(cf.getValue("outstep")); + //double endstep = UbSystem::stringTo<double>(cf.getValue("endstep")); + //int numOfThreads = UbSystem::stringTo<int>(cf.getValue("threads")); + + ConfigurationFile config; + config.load(configname); + + string pathname = config.getValue<string>("path"); + double availMem = config.getValue<double>("memory"); + string metafile = config.getValue<string>("metafile"); + double outstep = config.getValue<double>("outstep"); + double endstep = config.getValue<double>("endstep"); + int numOfThreads = config.getValue<int>("threads"); + const int refineLevel = config.getValue<int>("level"); + + bool test = config.getValue<bool>("test"); - pathname = cf.getValue("path"); - availMem = UbSystem::stringTo<double>(cf.getValue("memory")); - string metafile = cf.getValue("metafile"); - double outstep = UbSystem::stringTo<double>(cf.getValue("outstep")); - double endstep = UbSystem::stringTo<double>(cf.getValue("endstep")); - int numOfThreads = UbSystem::stringTo<int>(cf.getValue("threads")); LBMUnitConverterPtr conv = LBMUnitConverterPtr(new LBMUnitConverter()); @@ -44,9 +56,9 @@ void chanel(const char *cstr) const int blocknx2 = 8; const int blocknx3 = 8; - const int gridNx1 = 8;// 18; - const int gridNx2 = 4;// 11; - const int gridNx3 = 4;// 11; + const int gridNx1 = 8;//18; + const int gridNx2 = 8;// 11; + const int gridNx3 = 8;// 11; //const int blocknx1 = 40; //const int blocknx2 = 40; @@ -61,15 +73,15 @@ void chanel(const char *cstr) L2 = gridNx2*blocknx1; L3 = gridNx3*blocknx1; - LBMReal radius = 2.5; - LBMReal uLB = 0.0000001; - LBMReal Re = 0.0001; + LBMReal radius = 4; + LBMReal uLB = 0.1; + LBMReal Re = 1; LBMReal rhoLB = 0.0; //LBMReal nuLB = (uLB*2.0*radius)/Re; //LBMReal nuLB = (uLB*L2)/Re; - LBMReal nuLB = 0.168666666667/100; + LBMReal nuLB = 0.168666666667 / 100; - double dp_LB =1e-6; + double dp_LB = 1e-6; double rhoLBinflow = dp_LB*3.0; Grid3DPtr grid(new Grid3D(comm)); @@ -81,14 +93,13 @@ void chanel(const char *cstr) ////////////////////////////////////////////////////////////////////////// //restart UbSchedulerPtr restartSch(new UbScheduler(100000, 100000, 100000)); - RestartPostprocessor rp(grid, restartSch, comm, pathname, RestartPostprocessor::BINARY); + RestartCoProcessor rp(grid, restartSch, comm, pathname, RestartCoProcessor::BINARY); ////////////////////////////////////////////////////////////////////////// if (grid->getTimeStep() == 0) { const int baseLevel = 0; - const int refineLevel = UbSystem::stringTo<int>(cf.getValue("level")); //bounding box double d_minX1 = 0.0; @@ -101,21 +112,21 @@ void chanel(const char *cstr) double blockLength = blocknx1*dx; - if(myid ==0) + if (myid == 0) { - UBLOG(logINFO,"Parameters:"); - UBLOG(logINFO,"uLB = " << uLB ); - UBLOG(logINFO,"rhoLB = " << rhoLB ); - UBLOG(logINFO,"nueLB = " << nuLB ); - UBLOG(logINFO,"Re = " << Re ); - UBLOG(logINFO,"dx = " << dx ); - UBLOG(logINFO,"number of levels = " << refineLevel+1 ); - UBLOG(logINFO,"numOfThreads = " << numOfThreads ); - UBLOG(logINFO,"Preprozess - start"); + UBLOG(logINFO, "Parameters:"); + UBLOG(logINFO, "uLB = " << uLB); + UBLOG(logINFO, "rhoLB = " << rhoLB); + UBLOG(logINFO, "nueLB = " << nuLB); + UBLOG(logINFO, "Re = " << Re); + UBLOG(logINFO, "dx = " << dx); + UBLOG(logINFO, "number of levels = " << refineLevel + 1); + UBLOG(logINFO, "numOfThreads = " << numOfThreads); + UBLOG(logINFO, "Preprozess - start"); } GbObject3DPtr gridCube(new GbCuboid3D(d_minX1, d_minX2, d_minX3, d_maxX1, d_maxX2, d_maxX3)); - if(myid ==0) GbSystem3D::writeGeoObject(gridCube.get(),pathname + "/geo/gridCube", WbWriterVtkXmlBinary::getInstance()); + if (myid == 0) GbSystem3D::writeGeoObject(gridCube.get(), pathname + "/geo/gridCube", WbWriterVtkXmlBinary::getInstance()); GenBlocksGridVisitor genBlocks(gridCube); grid->accept(genBlocks); @@ -124,80 +135,78 @@ void chanel(const char *cstr) //GbObject3DPtr sphereRef(new GbSphere3D(L1/4.0, L2*0.5, L3*0.5, radius+1.0)); //GbSystem3D::writeGeoObject(sphereRef.get(),pathname + "/geo/sphereRef", WbWriterVtkXmlBinary::getInstance()); + //sphere - GbObject3DPtr sphere(new GbSphere3D(L1/4.0, L2*0.5, L3*0.5, radius)); + GbObject3DPtr sphere(new GbSphere3D(L1 / 2.0, L2*0.5, L3*0.5, radius)); + //GbObject3DPtr sphere(new GbSphere3D(L1/2.0-4.0, L2*0.5+4.0, L3*0.5+4.0, radius)); //GbObject3DPtr sphere(new GbCuboid3D(L1/4.0-radius, L2/2.0-radius, L3/2.0-radius, L1/4.0+radius, L2/2.0+radius, L3/2.0+radius)); - GbSystem3D::writeGeoObject(sphere.get(),pathname + "/geo/sphere", WbWriterVtkXmlBinary::getInstance()); + GbSystem3D::writeGeoObject(sphere.get(), pathname + "/geo/sphere", WbWriterVtkXmlBinary::getInstance()); double off = 0.0; - GbObject3DPtr refCube(new GbCuboid3D(sphere->getX1Minimum()-off, sphere->getX2Minimum()-off, sphere->getX3Minimum(), - sphere->getX1Maximum()+off, sphere->getX2Maximum()+off, sphere->getX3Maximum())); - if(myid ==0) GbSystem3D::writeGeoObject(refCube.get(),pathname + "/geo/refCube", WbWriterVtkXmlBinary::getInstance()); + GbObject3DPtr refCube(new GbCuboid3D(sphere->getX1Minimum() - off, sphere->getX2Minimum() - off, sphere->getX3Minimum(), + sphere->getX1Maximum() + off, sphere->getX2Maximum() + off, sphere->getX3Maximum())); + if (myid == 0) GbSystem3D::writeGeoObject(refCube.get(), pathname + "/geo/refCube", WbWriterVtkXmlBinary::getInstance()); if (refineLevel > 0) { - if(myid == 0) UBLOG(logINFO,"Refinement - start"); + if (myid == 0) UBLOG(logINFO, "Refinement - start"); RefineCrossAndInsideGbObjectHelper refineHelper(grid, refineLevel); - //refineHelper.addGbObject(sphere, refineLevel); - refineHelper.addGbObject(refCube, refineLevel); + refineHelper.addGbObject(sphere, refineLevel); + //refineHelper.addGbObject(refCube, refineLevel); refineHelper.refine(); - if(myid == 0) UBLOG(logINFO,"Refinement - end"); + if (myid == 0) UBLOG(logINFO, "Refinement - end"); } //walls - GbCuboid3DPtr addWallYmin (new GbCuboid3D(d_minX1-4.0*blockLength, d_minX2-4.0*blockLength, d_minX3-4.0*blockLength, d_maxX1+4.0*blockLength, d_minX2, d_maxX3+4.0*blockLength)); - if(myid == 0) GbSystem3D::writeGeoObject(addWallYmin.get(), pathname+"/geo/addWallYmin", WbWriterVtkXmlASCII::getInstance()); + GbCuboid3DPtr addWallYmin(new GbCuboid3D(d_minX1 - 4.0*blockLength, d_minX2 - 4.0*blockLength, d_minX3 - 4.0*blockLength, d_maxX1 + 4.0*blockLength, d_minX2, d_maxX3 + 4.0*blockLength)); + if (myid == 0) GbSystem3D::writeGeoObject(addWallYmin.get(), pathname + "/geo/addWallYmin", WbWriterVtkXmlASCII::getInstance()); - GbCuboid3DPtr addWallZmin (new GbCuboid3D(d_minX1-4.0*blockLength, d_minX2-4.0*blockLength, d_minX3-4.0*blockLength, d_maxX1+4.0*blockLength, d_maxX2+4.0*blockLength, d_minX3)); - if(myid == 0) GbSystem3D::writeGeoObject(addWallZmin.get(), pathname+"/geo/addWallZmin", WbWriterVtkXmlASCII::getInstance()); + GbCuboid3DPtr addWallZmin(new GbCuboid3D(d_minX1 - 4.0*blockLength, d_minX2 - 4.0*blockLength, d_minX3 - 4.0*blockLength, d_maxX1 + 4.0*blockLength, d_maxX2 + 4.0*blockLength, d_minX3)); + if (myid == 0) GbSystem3D::writeGeoObject(addWallZmin.get(), pathname + "/geo/addWallZmin", WbWriterVtkXmlASCII::getInstance()); - GbCuboid3DPtr addWallYmax (new GbCuboid3D(d_minX1-4.0*blockLength, d_maxX2, d_minX3-4.0*blockLength, d_maxX1+4.0*blockLength, d_maxX2+4.0*blockLength, d_maxX3+4.0*blockLength)); - if(myid == 0) GbSystem3D::writeGeoObject(addWallYmax.get(), pathname+"/geo/addWallYmax", WbWriterVtkXmlASCII::getInstance()); + GbCuboid3DPtr addWallYmax(new GbCuboid3D(d_minX1 - 4.0*blockLength, d_maxX2, d_minX3 - 4.0*blockLength, d_maxX1 + 4.0*blockLength, d_maxX2 + 4.0*blockLength, d_maxX3 + 4.0*blockLength)); + if (myid == 0) GbSystem3D::writeGeoObject(addWallYmax.get(), pathname + "/geo/addWallYmax", WbWriterVtkXmlASCII::getInstance()); - GbCuboid3DPtr addWallZmax (new GbCuboid3D(d_minX1-4.0*blockLength, d_minX2-4.0*blockLength, d_maxX3, d_maxX1+4.0*blockLength, d_maxX2+4.0*blockLength, d_maxX3+4.0*blockLength)); - if(myid == 0) GbSystem3D::writeGeoObject(addWallZmax.get(), pathname+"/geo/addWallZmax", WbWriterVtkXmlASCII::getInstance()); + GbCuboid3DPtr addWallZmax(new GbCuboid3D(d_minX1 - 4.0*blockLength, d_minX2 - 4.0*blockLength, d_maxX3, d_maxX1 + 4.0*blockLength, d_maxX2 + 4.0*blockLength, d_maxX3 + 4.0*blockLength)); + if (myid == 0) GbSystem3D::writeGeoObject(addWallZmax.get(), pathname + "/geo/addWallZmax", WbWriterVtkXmlASCII::getInstance()); //inflow - GbCuboid3DPtr geoInflow (new GbCuboid3D(d_minX1-4.0*blockLength, d_minX2-4.0*blockLength, d_minX3-4.0*blockLength, d_minX1, d_maxX2+4.0*blockLength, d_maxX3+4.0*blockLength)); - if(myid == 0) GbSystem3D::writeGeoObject(geoInflow.get(), pathname+"/geo/geoInflow", WbWriterVtkXmlASCII::getInstance()); + GbCuboid3DPtr geoInflow(new GbCuboid3D(d_minX1 - 4.0*blockLength, d_minX2 - 4.0*blockLength, d_minX3 - 4.0*blockLength, d_minX1, d_maxX2 + 4.0*blockLength, d_maxX3 + 4.0*blockLength)); + if (myid == 0) GbSystem3D::writeGeoObject(geoInflow.get(), pathname + "/geo/geoInflow", WbWriterVtkXmlASCII::getInstance()); //outflow - GbCuboid3DPtr geoOutflow (new GbCuboid3D(d_maxX1, d_minX2-4.0*blockLength, d_minX3-4.0*blockLength, d_maxX1+4.0*blockLength, d_maxX2+4.0*blockLength, d_maxX3+4.0*blockLength)); - if(myid == 0) GbSystem3D::writeGeoObject(geoOutflow.get(), pathname+"/geo/geoOutflow", WbWriterVtkXmlASCII::getInstance()); + GbCuboid3DPtr geoOutflow(new GbCuboid3D(d_maxX1, d_minX2 - 4.0*blockLength, d_minX3 - 4.0*blockLength, d_maxX1 + 4.0*blockLength, d_maxX2 + 4.0*blockLength, d_maxX3 + 4.0*blockLength)); + if (myid == 0) GbSystem3D::writeGeoObject(geoOutflow.get(), pathname + "/geo/geoOutflow", WbWriterVtkXmlASCII::getInstance()); - BlocksPostprocessorPtr ppblocks(new BlocksPostprocessor(grid, UbSchedulerPtr(new UbScheduler(1)), pathname, WbWriterVtkXmlBinary::getInstance(), comm)); + WriteBlocksCoProcessorPtr ppblocks(new WriteBlocksCoProcessor(grid, UbSchedulerPtr(new UbScheduler(1)), pathname, WbWriterVtkXmlBinary::getInstance(), comm)); //sphere - int bbOption = 1; //0=simple Bounce Back, 1=quadr. BB - D3Q27BoundaryConditionAdapterPtr bcNoSlip(new D3Q27NoSlipBCAdapter(bbOption)); - D3Q27InteractorPtr sphereInt = D3Q27InteractorPtr ( new D3Q27Interactor(sphere, grid, bcNoSlip,Interactor3D::SOLID)); - - D3Q27BoundaryConditionAdapterPtr bcSlip(new D3Q27SlipBCAdapter(bbOption)); + D3Q27BoundaryConditionAdapterPtr bcNoSlip(new D3Q27NoSlipBCAdapter()); + D3Q27InteractorPtr sphereInt = D3Q27InteractorPtr(new D3Q27Interactor(sphere, grid, bcNoSlip, Interactor3D::SOLID)); //walls - D3Q27InteractorPtr addWallYminInt(new D3Q27Interactor(addWallYmin, grid, bcNoSlip,Interactor3D::SOLID)); - D3Q27InteractorPtr addWallZminInt(new D3Q27Interactor(addWallZmin, grid, bcNoSlip,Interactor3D::SOLID)); - D3Q27InteractorPtr addWallYmaxInt(new D3Q27Interactor(addWallYmax, grid, bcNoSlip,Interactor3D::SOLID)); - D3Q27InteractorPtr addWallZmaxInt(new D3Q27Interactor(addWallZmax, grid, bcNoSlip,Interactor3D::SOLID)); + D3Q27InteractorPtr addWallYminInt(new D3Q27Interactor(addWallYmin, grid, bcNoSlip, Interactor3D::SOLID)); + D3Q27InteractorPtr addWallZminInt(new D3Q27Interactor(addWallZmin, grid, bcNoSlip, Interactor3D::SOLID)); + D3Q27InteractorPtr addWallYmaxInt(new D3Q27Interactor(addWallYmax, grid, bcNoSlip, Interactor3D::SOLID)); + D3Q27InteractorPtr addWallZmaxInt(new D3Q27Interactor(addWallZmax, grid, bcNoSlip, Interactor3D::SOLID)); mu::Parser fct; fct.SetExpr("U"); fct.DefineConst("U", uLB); - + //inflow - //D3Q27BoundaryConditionAdapterPtr velBCAdapter(new D3Q27VelocityBCAdapter (true, false ,false ,fct, 0, D3Q27BCFunction::INFCONST)); - //velBCAdapter->setSecondaryBcOption(2); - //D3Q27InteractorPtr inflowInt = D3Q27InteractorPtr( new D3Q27Interactor(geoInflow, grid, velBCAdapter, Interactor3D::SOLID)); + D3Q27BoundaryConditionAdapterPtr velBCAdapter(new D3Q27VelocityBCAdapter(true, false, false, fct, 0, D3Q27BCFunction::INFCONST)); + D3Q27InteractorPtr inflowInt = D3Q27InteractorPtr(new D3Q27Interactor(geoInflow, grid, velBCAdapter, Interactor3D::SOLID)); - D3Q27BoundaryConditionAdapterPtr denBCAdapterInflow(new D3Q27DensityBCAdapter(rhoLBinflow)); - denBCAdapterInflow->setSecondaryBcOption(0); - D3Q27InteractorPtr inflowInt = D3Q27InteractorPtr(new D3Q27Interactor(geoInflow, grid, denBCAdapterInflow, Interactor3D::SOLID)); + //D3Q27BoundaryConditionAdapterPtr denBCAdapterInflow(new D3Q27DensityBCAdapter(rhoLBinflow)); + //denBCAdapterInflow->setSecondaryBcOption(0); + //D3Q27InteractorPtr inflowInt = D3Q27InteractorPtr(new D3Q27Interactor(geoInflow, grid, denBCAdapterInflow, Interactor3D::SOLID)); //outflow D3Q27BoundaryConditionAdapterPtr denBCAdapter(new D3Q27DensityBCAdapter(rhoLB)); - D3Q27InteractorPtr outflowInt = D3Q27InteractorPtr( new D3Q27Interactor(geoOutflow, grid, denBCAdapter,Interactor3D::SOLID)); + D3Q27InteractorPtr outflowInt = D3Q27InteractorPtr(new D3Q27Interactor(geoOutflow, grid, denBCAdapter, Interactor3D::SOLID)); Grid3DVisitorPtr metisVisitor(new MetisPartitioningGridVisitor(comm, MetisPartitioningGridVisitor::LevelBased, D3Q27System::B)); InteractorsHelper intHelper(grid, metisVisitor); @@ -217,48 +226,48 @@ void chanel(const char *cstr) //set connectors D3Q27InterpolationProcessorPtr iProcessor(new D3Q27IncompressibleOffsetInterpolationProcessor()); - D3Q27SetConnectorsBlockVisitor setConnsVisitor(comm, true, D3Q27System::ENDDIR, nuLB, iProcessor); - grid->accept( setConnsVisitor ); + //D3Q27SetConnectorsBlockVisitor setConnsVisitor(comm, true, D3Q27System::ENDDIR, nuLB, iProcessor); + //grid->accept(setConnsVisitor); Block3DConnectorFactoryPtr factory(new Block3DConnectorFactory()); - //ConnectorBlockVisitor setConnsVisitor(comm, nuLB, iProcessor, factory); - //grid->accept(setConnsVisitor); + ConnectorBlockVisitor setConnsVisitor(comm, nuLB, iProcessor, factory); + grid->accept(setConnsVisitor); - ppblocks->update(0); + ppblocks->process(0); ppblocks.reset(); unsigned long nob = grid->getNumberOfBlocks(); int gl = 3; - unsigned long nod = nob * (blocknx1+gl) * (blocknx2+gl) * (blocknx3+gl); + unsigned long nod = nob * (blocknx1 + gl) * (blocknx2 + gl) * (blocknx3 + gl); - double needMemAll = double(nod*(27*sizeof(double) + sizeof(int) + sizeof(float)*4)); - double needMem = needMemAll / double(comm->getNumberOfProcesses()); + double needMemAll = double(nod*(27 * sizeof(double) + sizeof(int) + sizeof(float) * 4)); + double needMem = needMemAll / double(comm->getNumberOfProcesses()); - if(myid == 0) + if (myid == 0) { - UBLOG(logINFO,"Number of blocks = " << nob); - UBLOG(logINFO,"Number of nodes = " << nod); - UBLOG(logINFO,"Necessary memory = " << needMemAll << " bytes"); - UBLOG(logINFO,"Necessary memory per process = " << needMem << " bytes"); - UBLOG(logINFO,"Available memory per process = " << availMem << " bytes"); - } + UBLOG(logINFO, "Number of blocks = " << nob); + UBLOG(logINFO, "Number of nodes = " << nod); + UBLOG(logINFO, "Necessary memory = " << needMemAll << " bytes"); + UBLOG(logINFO, "Necessary memory per process = " << needMem << " bytes"); + UBLOG(logINFO, "Available memory per process = " << availMem << " bytes"); + } LBMKernel3DPtr kernel(new LBMKernelETD3Q27CCLB(blocknx1, blocknx2, blocknx3, LBMKernelETD3Q27CCLB::NORMAL)); - BCProcessorPtr bcProc(new D3Q27ETBCProcessor()); - //BoundaryConditionPtr velocityBC(new VelocityBoundaryCondition(refineLevel+1)); + BCProcessorPtr bcProcessor(new D3Q27ETBCProcessor()); + BoundaryConditionPtr velocityBC(new VelocityBoundaryCondition()); BoundaryConditionPtr densityBC(new NonEqDensityBoundaryCondition()); - //BoundaryConditionPtr noSlipBC(new NoSlipBoundaryCondition()); + BoundaryConditionPtr noSlipBC(new NoSlipBoundaryCondition()); //BoundaryConditionPtr velocityBC(new NonReflectingVelocityBoundaryCondition()); //BoundaryConditionPtr densityBC(new NonReflectingDensityBoundaryCondition()); - BoundaryConditionPtr noSlipBC(new HighViscosityNoSlipBoundaryCondition()); + //BoundaryConditionPtr noSlipBC(new HighViscosityNoSlipBoundaryCondition()); - //bcProcessor->addBC(velocityBC); + bcProcessor->addBC(velocityBC); bcProcessor->addBC(densityBC); bcProcessor->addBC(noSlipBC); - kernel->setBCProcessor(bcProc); + kernel->setBCProcessor(bcProcessor); SetKernelBlockVisitor kernelVisitor(kernel, nuLB, availMem, needMem); grid->accept(kernelVisitor); @@ -271,7 +280,6 @@ void chanel(const char *cstr) intHelper.setBC(); - BoundaryConditionBlockVisitor bcVisitor; grid->accept(bcVisitor); @@ -279,69 +287,69 @@ void chanel(const char *cstr) fctRoh.SetExpr("(x1max-x1)/l*dp*3.0"); fctRoh.DefineConst("dp", dp_LB); fctRoh.DefineConst("x1max", d_maxX1); - fctRoh.DefineConst("l", d_maxX1-d_minX1); + fctRoh.DefineConst("l", d_maxX1 - d_minX1); //initialization of distributions D3Q27ETInitDistributionsBlockVisitor initVisitor(nuLB, rhoLB); - //initVisitor.setVx1(fct); - initVisitor.setRho(fctRoh); + initVisitor.setVx1(fct); + //initVisitor.setRho(fctRoh); grid->accept(initVisitor); //Postrozess UbSchedulerPtr geoSch(new UbScheduler(1)); - D3Q27MacroscopicQuantitiesPostprocessorPtr ppgeo( - new D3Q27MacroscopicQuantitiesPostprocessor(grid, geoSch, pathname, WbWriterVtkXmlBinary::getInstance(), conv, true)); - ppgeo->update(0); + MacroscopicQuantitiesCoProcessorPtr ppgeo( + new MacroscopicQuantitiesCoProcessor(grid, geoSch, pathname, WbWriterVtkXmlBinary::getInstance(), conv, true)); + ppgeo->process(0); ppgeo.reset(); - if(myid == 0) UBLOG(logINFO,"Preprozess - end"); + if (myid == 0) UBLOG(logINFO, "Preprozess - end"); } else { - UBLOG(logINFO,"SetConnectors - start, id="<<myid); + UBLOG(logINFO, "SetConnectors - start, id=" << myid); //set connectors D3Q27InterpolationProcessorPtr iProcessor(new D3Q27IncompressibleOffsetInterpolationProcessor()); - D3Q27SetConnectorsBlockVisitor setConnsVisitor(comm, true, D3Q27System::ENDDIR, nuLB, iProcessor); - //ConnectorFactoryPtr cFactory(new Block3DConnectorFactory()); - //ConnectorBlockVisitor setConnsVisitor(comm, nuLB, iProcessor, cFactory); - grid->accept( setConnsVisitor ); + //D3Q27SetConnectorsBlockVisitor setConnsVisitor(comm, true, D3Q27System::ENDDIR, nuLB, iProcessor); + ConnectorFactoryPtr cFactory(new Block3DConnectorFactory()); + ConnectorBlockVisitor setConnsVisitor(comm, nuLB, iProcessor, cFactory); + grid->accept(setConnsVisitor); - UBLOG(logINFO,"SetConnectors - end, id="<<myid); + UBLOG(logINFO, "SetConnectors - end, id=" << myid); } UbSchedulerPtr stepSch(new UbScheduler(outstep)); //stepSch->addSchedule(10000, 0, 1000000); - D3Q27MacroscopicQuantitiesPostprocessor pp(grid, stepSch, pathname, WbWriterVtkXmlBinary::getInstance(), conv); + MacroscopicQuantitiesCoProcessor pp(grid, stepSch, pathname, WbWriterVtkXmlBinary::getInstance(), conv); - //D3Q27MacroscopicQuantitiesPostprocessor ppg(grid, stepSch, pathname, WbWriterVtkXmlBinary::getInstance(), conv,true); + //D3Q27MacroscopicQuantitiesCoProcessor ppg(grid, stepSch, pathname, WbWriterVtkXmlBinary::getInstance(), conv,true); - //InSituVTKPostprocessor isp(grid, stepSch, metafile, conv); + //InSituVTKCoProcessor isp(grid, stepSch, metafile, conv); UbSchedulerPtr nupsSch(new UbScheduler(10, 30, 100)); - NUPSCounterPostprocessor npr(grid, nupsSch, numOfThreads, comm); + NUPSCounterCoProcessor npr(grid, nupsSch, numOfThreads, comm); - CalculationManagerPtr calculation(new CalculationManager(grid, numOfThreads, endstep, stepSch, bcProcessor)); + CalculationManagerPtr calculation(new CalculationManager(grid, numOfThreads, endstep, stepSch)); - if(myid == 0) - UBLOG(logINFO,"Simulation-start"); + if (myid == 0) + UBLOG(logINFO, "Simulation-start"); calculation->calculate(); - if(myid == 0) - UBLOG(logINFO,"Simulation-end"); + if (myid == 0) + UBLOG(logINFO, "Simulation-end"); } - catch(std::exception& e) + catch (std::exception& e) { cerr << e.what() << endl << flush; } - catch(std::string& s) + catch (std::string& s) { cerr << s << endl; } - catch(...) + catch (...) { cerr << "unknown exception" << endl; } @@ -351,93 +359,15 @@ void chanel(const char *cstr) ////////////////////////////////////////////////////////////////////////// int main(int argc, char* argv[]) { - //struct iNode - //{ - // int x1, x2, x3; - // LBMReal x1off, x2off, x3off; - //}; - - //typedef boost::shared_ptr<iNode> iNodePtr; - - - // - - //set <vector<int>> iNodesSetSender00; - //vector<double> offs; - //map <vector<int>> iNodesSetSender00; - - //{ - // iNodePtr inode1(new iNode()); - // inode1->x1 = 1; - // inode1->x3 = 2; - // inode1->x2 = 3; - // inode1->x1off = 0; - // inode1->x2off = 1; - // inode1->x3off = 2; - - // iNodesSetSender00.insert(inode1); - //} - - //{ - // iNodePtr inode1(new iNode()); - // inode1->x1 = 1; - // inode1->x3 = 2; - // inode1->x2 = 3; - // inode1->x1off = 0; - // inode1->x2off = 1; - // inode1->x3off = 2; - - // iNodesSetSender00.insert(inode1); - //} - - //std::pair<std::set<vector<int>>::iterator, bool> ret; - - // { - // vector<int> inv; - // inv.push_back(1); - // inv.push_back(2); - // inv.push_back(3); - // inv.push_back(1); - // inv.push_back(2); - // inv.push_back(3); - // ret=iNodesSetSender00.insert(inv); - // if (ret.second==true) - // { - // offs.push_back(0); - // offs.push_back(1); - // offs.push_back(2); - // } - // } - - // { - // vector<int> inv; - // inv.push_back(1); - // inv.push_back(1); - // inv.push_back(3); - // inv.push_back(1); - // inv.push_back(2); - // inv.push_back(3); - // ret=iNodesSetSender00.insert(inv); - // if (ret.second==true) - // { - // offs.push_back(0); - // offs.push_back(1); - // offs.push_back(2); - // } - // } - - // return 0; - - - if ( argv != NULL ) + if (argv != NULL) { - if (argc > 1) + if (argv[1] != NULL) { - chanel(argv[1]); + run(string(argv[1])); } else { - cout << "Configuration file must be set!: " << argv[0] << " <config file>" << endl << std::flush; + cout << "Configuration file is missing!" << endl; } } } diff --git a/source/VirtualFluids.h b/source/VirtualFluids.h index 317b92879..3df635fb9 100644 --- a/source/VirtualFluids.h +++ b/source/VirtualFluids.h @@ -104,7 +104,7 @@ #include <BoundaryCondition/D3Q27SlipBCAdapter.h> #include <BoundaryCondition/D3Q27VelocityBCAdapter.h> -#include <BoundaryCondition/BoundaryConditionProcessor.h> +//#include <BoundaryCondition/BoundaryConditionProcessor.h> #include <BoundaryCondition/BoundaryCondition.h> #include <BoundaryCondition/VelocityBoundaryCondition.h> #include <BoundaryCondition/NonEqDensityBoundaryCondition.h> @@ -147,7 +147,7 @@ //#include <Grid/BoostSerializationClassExportHelper.h> #include <Grid/CalculationManager.h> #include <Grid/Calculator.h> -#include <Grid/Calculator2.h> +//#include <Grid/Calculator2.h> #include <Grid/Grid3D.h> #include <Grid/Grid3DSystem.h> #include <Interactors/D3Q27Interactor.h> diff --git a/source/VirtualFluidsCore/BoundaryCondition/BCArray.cpp b/source/VirtualFluidsCore/BoundaryCondition/BCArray.cpp new file mode 100644 index 000000000..b3a66d3b2 --- /dev/null +++ b/source/VirtualFluidsCore/BoundaryCondition/BCArray.cpp @@ -0,0 +1,196 @@ +#include "BCArray.h" + + +const int BCArray::SOLID = -1; +const int BCArray::FLUID = -2; +const int BCArray::INTERFACECF = -3; +const int BCArray::INTERFACEFC = -4; +const int BCArray::UNDEFINED = -5; + +////////////////////////////////////////////////////////////////////////// +BCArray::BCArray() {} +////////////////////////////////////////////////////////////////////////// +BCArray::BCArray(std::size_t nx1, std::size_t nx2, std::size_t nx3) +{ + bcindexmatrix.resize(nx1, nx2, nx3, UNDEFINED); +} +////////////////////////////////////////////////////////////////////////// +BCArray::BCArray(std::size_t nx1, std::size_t nx2, std::size_t nx3, int val) +{ + bcindexmatrix.resize(nx1, nx2, nx3, val); +} +////////////////////////////////////////////////////////////////////////// +BCArray::~BCArray() {} +////////////////////////////////////////////////////////////////////////// +void BCArray::resize(std::size_t nx1, std::size_t nx2, std::size_t nx3) +{ + bcindexmatrix.resize(nx1, nx2, nx3); +} +////////////////////////////////////////////////////////////////////////// +void BCArray::resize(std::size_t nx1, std::size_t nx2, std::size_t nx3, int val) +{ + bcindexmatrix.resize(nx1, nx2, nx3, val); +} +////////////////////////////////////////////////////////////////////////// +bool BCArray::validIndices(std::size_t x1, std::size_t x2, std::size_t x3) const +{ + if (x1 < 0 || x1 >= this->bcindexmatrix.getNX1()) return false; + if (x2 < 0 || x2 >= this->bcindexmatrix.getNX2()) return false; + if (x3 < 0 || x3 >= this->bcindexmatrix.getNX3()) return false; + return true; +} +////////////////////////////////////////////////////////////////////////// +void BCArray::setBC(std::size_t x1, std::size_t x2, std::size_t x3, BCClassPtr const& bc) +{ + if (this->hasBC(x1, x2, x3)) + { + if (this->getBC(x1, x2, x3) == bc) return; + else this->deleteBC(x1, x2, x3); + } + + //wenn keine frei gewordene BCs vorhanden + if (indexContainer.empty()) + { + bcvector.push_back(bc); + bcindexmatrix(x1, x2, x3) = (int)bcvector.size() - 1; + } + else + { + int index = indexContainer.back(); + bcindexmatrix(x1, x2, x3) = index; + bcvector[index] = bc; + indexContainer.pop_back(); + } +} +////////////////////////////////////////////////////////////////////////// +void BCArray::setSolid(std::size_t x1, std::size_t x2, std::size_t x3) +{ + if (this->hasBC(x1, x2, x3)) this->deleteBC(x1, x2, x3); + bcindexmatrix(x1, x2, x3) = SOLID; +} +////////////////////////////////////////////////////////////////////////// +void BCArray::setFluid(std::size_t x1, std::size_t x2, std::size_t x3) +{ + if (this->hasBC(x1, x2, x3)) this->deleteBC(x1, x2, x3); + bcindexmatrix(x1, x2, x3) = FLUID; +} +////////////////////////////////////////////////////////////////////////// +void BCArray::setUndefined(std::size_t x1, std::size_t x2, std::size_t x3) +{ + if (this->hasBC(x1, x2, x3)) this->deleteBC(x1, x2, x3); + bcindexmatrix(x1, x2, x3) = UNDEFINED; +} +////////////////////////////////////////////////////////////////////////// +std::size_t BCArray::getNumberOfSolidEntries() const +{ + const std::vector<int>& data = bcindexmatrix.getDataVector(); + std::size_t counter = 0; + for (std::size_t i = 0; i < data.size(); i++) + if (data[i] == SOLID) counter++; + return counter; +} +////////////////////////////////////////////////////////////////////////// +std::size_t BCArray::getNumberOfFluidEntries() const +{ + const std::vector<int>& data = bcindexmatrix.getDataVector(); + std::size_t counter = 0; + for (std::size_t i = 0; i < data.size(); i++) + { + int tmp = data[i]; + if (tmp == FLUID || tmp >= 0) counter++; + } + return counter; +} +////////////////////////////////////////////////////////////////////////// +std::size_t BCArray::getNumberOfFluidWithoutBCEntries() const +{ + const std::vector<int>& data = bcindexmatrix.getDataVector(); + std::size_t counter = 0; + for (std::size_t i = 0; i < data.size(); i++) + if (data[i] == FLUID) counter++; + return counter; +} +////////////////////////////////////////////////////////////////////////// +std::size_t BCArray::getNumberOfBCEntries() const +{ + const std::vector<int>& data = bcindexmatrix.getDataVector(); + std::size_t counter = 0; + for (std::size_t i = 0; i < data.size(); i++) + if (data[i] >= 0) counter++; + return counter; +} +////////////////////////////////////////////////////////////////////////// +std::size_t BCArray::getNumberOfUndefinedEntries() const +{ + const std::vector<int>& data = bcindexmatrix.getDataVector(); + std::size_t counter = 0; + for (std::size_t i = 0; i < data.size(); i++) + if (data[i] == UNDEFINED) counter++; + return counter; +} +////////////////////////////////////////////////////////////////////////// +std::size_t BCArray::getBCVectorSize() const +{ + return this->bcvector.size(); +} +////////////////////////////////////////////////////////////////////////// +std::string BCArray::toString() const +{ + std::size_t solidCounter = 0; + std::size_t fluidCounter = 0; + std::size_t bcCounter = 0; + std::size_t undefCounter = 0; + + for (int x1 = 0; x1 < bcindexmatrix.getNX1(); x1++) + { + for (int x2 = 0; x2 < bcindexmatrix.getNX2(); x2++) + { + for (int x3 = 0; x3 < bcindexmatrix.getNX3(); x3++) + { + if (bcindexmatrix(x1, x2, x3) >= 0) bcCounter++; + else if (bcindexmatrix(x1, x2, x3) == FLUID) fluidCounter++; + else if (bcindexmatrix(x1, x2, x3) == SOLID) solidCounter++; + else if (bcindexmatrix(x1, x2, x3) == UNDEFINED) undefCounter++; + else throw UbException(UB_EXARGS, "invalid matrixEntry"); + } + } + } + + std::size_t unrefEntriesInBcVector = 0; + for (std::size_t i = 0; i < bcvector.size(); i++) if (!bcvector[i]) unrefEntriesInBcVector++; + + std::stringstream text; + text << "BCArray<" << typeid(BCClassPtr).name() << "," << typeid(int).name() << ">"; + text << "[ entries: " << bcindexmatrix.getNX1() << "x" << bcindexmatrix.getNX2(); + text << "x" << bcindexmatrix.getNX3() << "="; + text << bcindexmatrix.getNX1()*bcindexmatrix.getNX2()*bcindexmatrix.getNX3() << " ]:\n"; + text << " - #fluid entries : " << fluidCounter << std::endl; + text << " - #bc entries : " << bcCounter << std::endl; + text << " - #solid entries : " << solidCounter << std::endl; + text << " - #undef entries : " << undefCounter << std::endl; + text << " - bcvector-entries : " << bcvector.size() << " (empty ones: " << unrefEntriesInBcVector << ")\n"; + text << " - indexContainer-entries: " << indexContainer.size() << std::endl; + + return text.str(); +} +////////////////////////////////////////////////////////////////////////// +void BCArray::deleteBCAndSetType(std::size_t x1, std::size_t x2, std::size_t x3, int type) + { + this->deleteBC(x1, x2, x3); + + //matrix neuen Typ zuweisen + bcindexmatrix(x1, x2, x3) = type; + } +////////////////////////////////////////////////////////////////////////// +void BCArray::deleteBC(std::size_t x1, std::size_t x2, std::size_t x3) + { + //ueberpruefen, ob ueberhaupt BC vorhanden + int index = bcindexmatrix(x1, x2, x3); + if (index < 0) return; + + //frei gewordenen Index in den Indexcontainer schieben + indexContainer.push_back(index); + + //element "loeschen" + bcvector[index] = BCClassPtr(); + } \ No newline at end of file diff --git a/source/VirtualFluidsCore/BoundaryCondition/BCArray.h b/source/VirtualFluidsCore/BoundaryCondition/BCArray.h index 704d0df69..f4090c206 100644 --- a/source/VirtualFluidsCore/BoundaryCondition/BCArray.h +++ b/source/VirtualFluidsCore/BoundaryCondition/BCArray.h @@ -1,265 +1,103 @@ #ifndef BCArray_H #define BCArray_H +#include "D3Q27BoundaryCondition.h" +#include "basics/container/CbArray3D.h" + #include <typeinfo> #include <boost/serialization/serialization.hpp> #include <boost/shared_ptr.hpp> +typedef boost::shared_ptr<D3Q27BoundaryCondition> BCClassPtr; + +class BCArray; +typedef boost::shared_ptr<BCArray> BCArrayPtr; + class BCArray { public: - typedef typename boost::shared_ptr<BCArray> BCArrayPtr; - typedef std::size_t size_type; - typedef int IndexType; - -public: ////////////////////////////////////////////////////////////////////////// - BCArray() {} + BCArray(); + ////////////////////////////////////////////////////////////////////////// + BCArray(std::size_t nx1, std::size_t nx2, std::size_t nx3); + ////////////////////////////////////////////////////////////////////////// + BCArray(std::size_t nx1, std::size_t nx2, std::size_t nx3, int val); ////////////////////////////////////////////////////////////////////////// virtual ~BCArray(); ////////////////////////////////////////////////////////////////////////// - virtual inline size_type getNX1() const = 0; + inline std::size_t getNX1() const; ////////////////////////////////////////////////////////////////////////// - virtual inline size_type getNX2() const = 0; + inline std::size_t getNX2() const; ////////////////////////////////////////////////////////////////////////// - virtual inline size_type getNX3() const = 0; + inline std::size_t getNX3() const; ////////////////////////////////////////////////////////////////////////// - void resize(const size_type& nx1, const size_type& nx2, const size_type& nx3) = 0; + void resize(std::size_t nx1, std::size_t nx2, std::size_t nx3); ////////////////////////////////////////////////////////////////////////// - void resize(const size_type& nx1, const size_type& nx2, const size_type& nx3, const IndexType& val) = 0; + void resize(std::size_t nx1, std::size_t nx2, std::size_t nx3, int val); ////////////////////////////////////////////////////////////////////////// - bool validIndices(const size_type& x1, const size_type& x2, const size_type& x3) const = 0; + bool validIndices(std::size_t x1, std::size_t x2, std::size_t x3) const; ////////////////////////////////////////////////////////////////////////// - inline bool hasBC(const size_type& x1, const size_type& x2, const size_type& x3) const = 0; + inline bool hasBC(std::size_t x1, std::size_t x2, std::size_t x3) const; ////////////////////////////////////////////////////////////////////////// - void setBC(const size_type& x1, const size_type& x2, const size_type& x3, BCClassPtr const& bc) - { - if( this->hasBC(x1,x2,x3) ) - { - if( this->getBC(x1,x2,x3)==bc ) return; - else this->deleteBC(x1,x2,x3); - } - - //wenn keine frei gewordene BCs vorhanden - if( indexContainer.empty() ) - { - bcvector.push_back(bc); - bcindexmatrix(x1,x2,x3) = (IndexType)bcvector.size()-1; - } - else - { - IndexType index = indexContainer.back(); - bcindexmatrix(x1,x2,x3) = index; - bcvector[index] = bc; - indexContainer.pop_back(); - } - } + void setBC(std::size_t x1, std::size_t x2, std::size_t x3, BCClassPtr const& bc); ////////////////////////////////////////////////////////////////////////// - inline IndexType getBCVectorIndex(const size_type& x1, const size_type& x2, const size_type& x3) const - { - return bcindexmatrix(x1,x2,x3); - } + inline int getBCVectorIndex(std::size_t x1, std::size_t x2, std::size_t x3) const; ////////////////////////////////////////////////////////////////////////// - inline const BCClassPtr getBC(const size_type& x1, const size_type& x2, const size_type& x3) const - { - IndexType index = bcindexmatrix(x1,x2,x3); - if(index<0) return BCClassPtr(); //=> NULL Pointer - - return bcvector[index]; - } + inline const BCClassPtr getBC(std::size_t x1, std::size_t x2, std::size_t x3) const; ////////////////////////////////////////////////////////////////////////// - inline BCClassPtr getBC(const size_type& x1, const size_type& x2, const size_type& x3) - { - IndexType index = bcindexmatrix(x1,x2,x3); - if(index<0) return BCClassPtr(); //=> NULL Pointer - - return bcvector[index]; - } + inline BCClassPtr getBC(std::size_t x1, std::size_t x2, std::size_t x3); ////////////////////////////////////////////////////////////////////////// - void setSolid(const size_type& x1, const size_type& x2, const size_type& x3) - { - if( this->hasBC(x1,x2,x3) ) this->deleteBC(x1,x2,x3); - bcindexmatrix(x1,x2,x3)=SOLID; - } + void setSolid(std::size_t x1, std::size_t x2, std::size_t x3); ////////////////////////////////////////////////////////////////////////// - inline bool isSolid(const size_type& x1, const size_type& x2, const size_type& x3) const - { - return bcindexmatrix(x1,x2,x3)==SOLID; - } + inline bool isSolid(std::size_t x1, std::size_t x2, std::size_t x3) const; ////////////////////////////////////////////////////////////////////////// - void setFluid(const size_type& x1, const size_type& x2, const size_type& x3) - { - if( this->hasBC(x1,x2,x3) ) this->deleteBC(x1,x2,x3); - bcindexmatrix(x1,x2,x3)=FLUID; - } + void setFluid(std::size_t x1, std::size_t x2, std::size_t x3); ////////////////////////////////////////////////////////////////////////// //true : FLUID or BC //false: UNDEFINED or SOLID - inline bool isFluid(const size_type& x1, const size_type& x2, const size_type& x3) const - { - int tmp = bcindexmatrix(x1,x2,x3); - return (tmp==FLUID || tmp>=0); - } + inline bool isFluid(std::size_t x1, std::size_t x2, std::size_t x3) const; ////////////////////////////////////////////////////////////////////////// - inline bool isFluidWithoutBC(const size_type& x1, const size_type& x2, const size_type& x3) const - { - return bcindexmatrix(x1,x2,x3)==FLUID; - } + inline bool isFluidWithoutBC(std::size_t x1, std::size_t x2, std::size_t x3) const; ////////////////////////////////////////////////////////////////////////// - inline bool isUndefined(const size_type& x1, const size_type& x2, const size_type& x3) const - { - return bcindexmatrix(x1,x2,x3)==UNDEFINED; - } + inline bool isUndefined(std::size_t x1, std::size_t x2, std::size_t x3) const; ////////////////////////////////////////////////////////////////////////// - void setUndefined(const size_type& x1, const size_type& x2, const size_type& x3) - { - if( this->hasBC(x1,x2,x3) ) this->deleteBC(x1,x2,x3); - bcindexmatrix(x1,x2,x3)=UNDEFINED; - } + void setUndefined(std::size_t x1, std::size_t x2, std::size_t x3); ////////////////////////////////////////////////////////////////////////// - inline bool isInterface(const size_type& x1, const size_type& x2, const size_type& x3) const - { - return bcindexmatrix(x1,x2,x3)==INTERFACE; - } + inline bool isInterfaceCF(std::size_t x1, std::size_t x2, std::size_t x3) const; ////////////////////////////////////////////////////////////////////////// - void setInterface(const size_type& x1, const size_type& x2, const size_type& x3) - { - if( this->hasBC(x1,x2,x3) ) this->deleteBC(x1,x2,x3); - bcindexmatrix(x1,x2,x3)=INTERFACE; - } + void setInterfaceCF(std::size_t x1, std::size_t x2, std::size_t x3); ////////////////////////////////////////////////////////////////////////// - std::size_t getNumberOfSolidEntries() const - { - const std::vector<IndexType>& data = bcindexmatrix.getDataVector(); - std::size_t counter = 0; - for(std::size_t i=0; i<data.size(); i++) - if(data[i]==SOLID) counter++; - return counter; - } + inline bool isInterfaceFC(std::size_t x1, std::size_t x2, std::size_t x3) const; ////////////////////////////////////////////////////////////////////////// - std::size_t getNumberOfFluidEntries() const - { - const std::vector<IndexType>& data = bcindexmatrix.getDataVector(); - std::size_t counter = 0; - for(std::size_t i=0; i<data.size(); i++) - { - const IndexType& tmp = data[i]; - if(tmp==FLUID || tmp>=0) counter++; - } - return counter; - } + void setInterfaceFC(std::size_t x1, std::size_t x2, std::size_t x3); ////////////////////////////////////////////////////////////////////////// - std::size_t getNumberOfFluidWithoutBCEntries() const - { - const std::vector<IndexType>& data = bcindexmatrix.getDataVector(); - std::size_t counter = 0; - for(std::size_t i=0; i<data.size(); i++) - if(data[i]==FLUID) counter++; - return counter; - } + std::size_t getNumberOfSolidEntries() const; ////////////////////////////////////////////////////////////////////////// - std::size_t getNumberOfBCEntries() const - { - const std::vector<IndexType>& data = bcindexmatrix.getDataVector(); - std::size_t counter = 0; - for(std::size_t i=0; i<data.size(); i++) - if(data[i]>=0) counter++; - return counter; - } + std::size_t getNumberOfFluidEntries() const; ////////////////////////////////////////////////////////////////////////// - std::size_t getNumberOfUndefinedEntries() const - { - const std::vector<IndexType>& data = bcindexmatrix.getDataVector(); - std::size_t counter = 0; - for(std::size_t i=0; i<data.size(); i++) - if(data[i]==UNDEFINED) counter++; - return counter; - } + std::size_t getNumberOfFluidWithoutBCEntries() const; ////////////////////////////////////////////////////////////////////////// - std::size_t getBCVectorSize() const - { - return this->bcvector.size(); - } + std::size_t getNumberOfBCEntries() const; ////////////////////////////////////////////////////////////////////////// - std::string toString() const - { - std::size_t solidCounter = 0; - std::size_t fluidCounter = 0; - std::size_t bcCounter = 0; - std::size_t undefCounter = 0; - - for(int x1=0; x1<bcindexmatrix.getNX1(); x1++) - { - for(int x2=0; x2<bcindexmatrix.getNX2(); x2++) - { - for(int x3=0; x3<bcindexmatrix.getNX3(); x3++) - { - if(bcindexmatrix(x1,x2,x3)>=0 ) bcCounter++; - else if(bcindexmatrix(x1,x2,x3)==FLUID ) fluidCounter++; - else if(bcindexmatrix(x1,x2,x3)==SOLID ) solidCounter++; - else if(bcindexmatrix(x1,x2,x3)==UNDEFINED) undefCounter++; - else throw UbException(UB_EXARGS,"invalid matrixEntry"); - } - } - } - - std::size_t unrefEntriesInBcVector=0; - for(std::size_t i=0; i<bcvector.size(); i++) if(!bcvector[i]) unrefEntriesInBcVector++; - - std::stringstream text; - text<<"BCArray<"<<typeid(BCClass).name()<<","<<typeid(IndexType).name()<<">"; - text<<"[ entries: "<<bcindexmatrix.getNX1()<<"x"<<bcindexmatrix.getNX2(); - text<<"x"<<bcindexmatrix.getNX3()<<"="; - text<<bcindexmatrix.getNX1()*bcindexmatrix.getNX2()*bcindexmatrix.getNX3()<<" ]:\n"; - text<<" - #fluid entries : "<<fluidCounter<<std::endl; - text<<" - #bc entries : "<<bcCounter<<std::endl; - text<<" - #solid entries : "<<solidCounter<<std::endl; - text<<" - #undef entries : "<<undefCounter<<std::endl; - text<<" - bcvector-entries : "<<bcvector.size()<<" (empty ones: "<<unrefEntriesInBcVector<<")\n"; - text<<" - indexContainer-entries: "<<indexContainer.size()<<std::endl; - - return text.str(); - } + std::size_t getNumberOfUndefinedEntries() const; + ////////////////////////////////////////////////////////////////////////// + std::size_t getBCVectorSize() const; + ////////////////////////////////////////////////////////////////////////// + std::string toString() const; ////////////////////////////////////////////////////////////////////////// -#ifdef CAB_RCF - template<class Archive> - void SF_SERIALIZE(Archive & ar) - { - ar & bcindexmatrix; - ar & bcvector; - ar & indexContainer; - } -#endif //CAB_RCF - - static const IndexType SOLID; //definiton erfolgt ausserhalb!!! - static const IndexType FLUID; //definiton erfolgt ausserhalb!!! - static const IndexType INTERFACE; //definiton erfolgt ausserhalb!!! - static const IndexType UNDEFINED; //definiton erfolgt ausserhalb!!! + static const int SOLID; + static const int FLUID; + static const int INTERFACECF; + static const int INTERFACEFC; + static const int UNDEFINED; private: ////////////////////////////////////////////////////////////////////////// - void deleteBCAndSetType(const size_type& x1, const size_type& x2, const size_type& x3, const IndexType& type) - { - this->deleteBC(x1,x2,x3); - - //matrix neuen Typ zuweisen - bcindexmatrix(x1,x2,x3) = type; - } + void deleteBCAndSetType(std::size_t x1, std::size_t x2, std::size_t x3, int type); ////////////////////////////////////////////////////////////////////////// - void deleteBC(const size_type& x1,const size_type& x2, const size_type& x3) - { - //ueberpruefen, ob ueberhaupt BC vorhanden - IndexType index = bcindexmatrix(x1,x2,x3); - if(index<0) return; - - //frei gewordenen Index in den Indexcontainer schieben - indexContainer.push_back(index); - - //element "loeschen" - bcvector[index] = BCClassPtr(); - } + void deleteBC(std::size_t x1, std::size_t x2, std::size_t x3); friend class boost::serialization::access; template<class Archive> @@ -271,19 +109,77 @@ private: } protected: ////////////////////////////////////////////////////////////////////////// - //CbArray3D<IndexType,IndexerX1X2X3> bcindexmatrix; //-1 solid // -2 fluid -... - CbArray3D<IndexType,IndexerX3X2X1> bcindexmatrix; + //-1 solid // -2 fluid -... + CbArray3D<int, IndexerX3X2X1> bcindexmatrix; std::vector<BCClassPtr> bcvector; - std::vector<IndexType> indexContainer; + std::vector<int> indexContainer; }; -template< typename BCClass , class IndexType> -const IndexType BCArray<BCClass,IndexType>::SOLID = -1; -template< typename BCClass , class IndexType> -const IndexType BCArray<BCClass,IndexType>::FLUID = -2; -template< typename BCClass , class IndexType> -const IndexType BCArray<BCClass,IndexType>::INTERFACE = -3; -template< typename BCClass , class IndexType> -const IndexType BCArray<BCClass,IndexType>::UNDEFINED = -4; -#endif //D3Q19BCMATRIX_H +////////////////////////////////////////////////////////////////////////// +inline std::size_t BCArray::getNX1() const { return bcindexmatrix.getNX1(); } +////////////////////////////////////////////////////////////////////////// +inline std::size_t BCArray::getNX2() const { return bcindexmatrix.getNX2(); } +////////////////////////////////////////////////////////////////////////// +inline std::size_t BCArray::getNX3() const { return bcindexmatrix.getNX3(); } +////////////////////////////////////////////////////////////////////////// +inline bool BCArray::hasBC(std::size_t x1, std::size_t x2, std::size_t x3) const +{ + return bcindexmatrix(x1, x2, x3) >= 0; +} +////////////////////////////////////////////////////////////////////////// +inline int BCArray::getBCVectorIndex(std::size_t x1, std::size_t x2, std::size_t x3) const +{ + return bcindexmatrix(x1, x2, x3); +} +////////////////////////////////////////////////////////////////////////// +inline const BCClassPtr BCArray::getBC(std::size_t x1, std::size_t x2, std::size_t x3) const +{ + int index = bcindexmatrix(x1, x2, x3); + if (index < 0) return BCClassPtr(); //=> NULL Pointer + + return bcvector[index]; +} +////////////////////////////////////////////////////////////////////////// +inline BCClassPtr BCArray::getBC(std::size_t x1, std::size_t x2, std::size_t x3) +{ + int index = bcindexmatrix(x1, x2, x3); + if (index < 0) return BCClassPtr(); //=> NULL Pointer + + return bcvector[index]; +} +////////////////////////////////////////////////////////////////////////// +inline bool BCArray::isSolid(std::size_t x1, std::size_t x2, std::size_t x3) const +{ + return bcindexmatrix(x1, x2, x3) == SOLID; +} +////////////////////////////////////////////////////////////////////////// +//true : FLUID or BC +//false: UNDEFINED or SOLID +inline bool BCArray::isFluid(std::size_t x1, std::size_t x2, std::size_t x3) const +{ + int tmp = bcindexmatrix(x1, x2, x3); + return (tmp == FLUID || tmp >= 0); +} +////////////////////////////////////////////////////////////////////////// +inline bool BCArray::isFluidWithoutBC(std::size_t x1, std::size_t x2, std::size_t x3) const +{ + return bcindexmatrix(x1, x2, x3) == FLUID; +} +////////////////////////////////////////////////////////////////////////// +inline bool BCArray::isUndefined(std::size_t x1, std::size_t x2, std::size_t x3) const +{ + return bcindexmatrix(x1, x2, x3) == UNDEFINED; +} +////////////////////////////////////////////////////////////////////////// +inline bool BCArray::isInterfaceCF(std::size_t x1, std::size_t x2, std::size_t x3) const +{ + return bcindexmatrix(x1, x2, x3) == INTERFACECF; +} +////////////////////////////////////////////////////////////////////////// +inline bool BCArray::isInterfaceFC(std::size_t x1, std::size_t x2, std::size_t x3) const +{ + return bcindexmatrix(x1, x2, x3) == INTERFACEFC; +} + +#endif diff --git a/source/VirtualFluidsCore/BoundaryCondition/BCProcessor.h b/source/VirtualFluidsCore/BoundaryCondition/BCProcessor.h index b77119268..5eea70a1b 100644 --- a/source/VirtualFluidsCore/BoundaryCondition/BCProcessor.h +++ b/source/VirtualFluidsCore/BoundaryCondition/BCProcessor.h @@ -18,7 +18,6 @@ class BCProcessor public: BCProcessor(); virtual ~BCProcessor(); - virtual void applyBC() = 0; virtual void applyPreCollisionBC() = 0; virtual void applyPostCollisionBC() = 0; virtual void addBC(BoundaryConditionPtr bc) = 0; diff --git a/source/VirtualFluidsCore/BoundaryCondition/BoundaryCondition.cpp b/source/VirtualFluidsCore/BoundaryCondition/BoundaryCondition.cpp index caba2663d..4148647b7 100644 --- a/source/VirtualFluidsCore/BoundaryCondition/BoundaryCondition.cpp +++ b/source/VirtualFluidsCore/BoundaryCondition/BoundaryCondition.cpp @@ -9,14 +9,12 @@ void BoundaryCondition::apply() { if (this->compressible) { - //rho0 = REAL_CAST(0.0); calcFeqsForDirFct = &D3Q27System::getCompFeqForDirection; calcMacrosFct = &D3Q27System::calcCompMacroscopicValues; calcFeqFct = &D3Q27System::calcCompFeq; } else { - //rho0 = REAL_CAST(1.0); calcFeqsForDirFct = &D3Q27System::getIncompFeqForDirection; calcMacrosFct = &D3Q27System::calcIncompMacroscopicValues; calcFeqFct = &D3Q27System::calcIncompFeq; @@ -25,12 +23,7 @@ void BoundaryCondition::apply() int cbc = 0; int cn = 0; int j; - //int dsize = (int)distVector.size(); - //for (int i = 0; i < dsize; i++) - //{ - // distributions = distVector[i]; - //int nsize = (int)sizeVector[i] * 3 + cn; int nsize = (int)nodeVector.size(); for (j = cn; j < nsize;) @@ -45,22 +38,11 @@ void BoundaryCondition::apply() applyBC(); } cn = j; - //} } -////////////////////////////////////////////////////////////////////////// -//void BoundaryCondition::apply() -//{ -//} -////////////////////////////////////////////////////////////////////////// -void BoundaryCondition::addDistributions(EsoTwist3DPtr distributions) -{ - //distVector.push_back(distribution); - this->distributions = distributions; -} -////////////////////////////////////////////////////////////////////////// -//void BoundaryCondition::addNnode(int nnode) +//////////////////////////////////////////////////////////////////////////// +//void BoundaryCondition::addDistributions(EsoTwist3DPtr distributions) //{ -// sizeVector.push_back(nnode); +// this->distributions = distributions; //} ////////////////////////////////////////////////////////////////////////// void BoundaryCondition::addNode(int x1, int x2, int x3) diff --git a/source/VirtualFluidsCore/BoundaryCondition/BoundaryCondition.h b/source/VirtualFluidsCore/BoundaryCondition/BoundaryCondition.h index 0f88d4a05..aa46db7f6 100644 --- a/source/VirtualFluidsCore/BoundaryCondition/BoundaryCondition.h +++ b/source/VirtualFluidsCore/BoundaryCondition/BoundaryCondition.h @@ -29,18 +29,8 @@ public: BoundaryCondition(); virtual ~BoundaryCondition() {} - //void apply(int level); - - //void addDistribution(int level, EsoTwist3DPtr distribution); - //void addNnode(int level, int nnode); - //void addNode(int level, int x1, int x2, int x3); - //void addBcPointer(int level, D3Q27BoundaryConditionPtr bcPtr); - //void setCompressible(bool c); - //void setCollFactor(int level, double cf); void apply(); - - void addDistributions(EsoTwist3DPtr distributions); - //void addNnode(int nnode); + virtual void addDistributions(DistributionArray3DPtr distributions) = 0; void addNode(int x1, int x2, int x3); void addBcPointer(D3Q27BoundaryConditionPtr bcPtr); void setCompressible(bool c); @@ -51,8 +41,6 @@ public: protected: virtual void applyBC() = 0; - //std::vector <EsoTwist3DPtr> distVector; - //std::vector <int> sizeVector; std::vector <int> nodeVector; std::vector <D3Q27BoundaryConditionPtr> bcVector; @@ -61,8 +49,8 @@ protected: bool preCollision; D3Q27BoundaryConditionPtr bcPtr; - EsoTwist3DPtr distributions; - EsoTwist3DPtr distributionsTemp; + DistributionArray3DPtr distributions; + DistributionArray3DPtr distributionsTemp; LBMReal collFactor; int x1, x2, x3; diff --git a/source/VirtualFluidsCore/BoundaryCondition/BoundaryConditionProcessor.cpp b/source/VirtualFluidsCore/BoundaryCondition/BoundaryConditionProcessor.cpp deleted file mode 100644 index 5234a86aa..000000000 --- a/source/VirtualFluidsCore/BoundaryCondition/BoundaryConditionProcessor.cpp +++ /dev/null @@ -1,58 +0,0 @@ -#include "BoundaryConditionProcessor.h" -#include <boost/foreach.hpp> - -BoundaryConditionProcessor::BoundaryConditionProcessor() -{ - -} -////////////////////////////////////////////////////////////////////////// -void BoundaryConditionProcessor::addBC(BoundaryConditionPtr bc) -{ - BoundaryCondition::Type type = bc->getType(); - if (bc->isPreCollision()) - { - preBC.push_back(bc); - } - else - { - postBC.push_back(bc); - } -} -////////////////////////////////////////////////////////////////////////// -BoundaryConditionPtr BoundaryConditionProcessor::getBC(BoundaryCondition::Type type) -{ - BOOST_FOREACH(BoundaryConditionPtr bc, preBC) - { - if (bc->getType() == type) - { - return bc; - } - } - - BOOST_FOREACH(BoundaryConditionPtr bc, postBC) - { - if (bc->getType() == type) - { - return bc; - } - } - - return BoundaryConditionPtr(); -} -////////////////////////////////////////////////////////////////////////// -void BoundaryConditionProcessor::applyPreCollisionBC() -{ - BOOST_FOREACH(BoundaryConditionPtr bc, preBC) - { - bc->apply(); - } -} -////////////////////////////////////////////////////////////////////////// -void BoundaryConditionProcessor::applyPostCollisionBC() -{ - BOOST_FOREACH(BoundaryConditionPtr bc, postBC) - { - bc->apply(); - } -} - diff --git a/source/VirtualFluidsCore/BoundaryCondition/BoundaryConditionProcessor.h b/source/VirtualFluidsCore/BoundaryCondition/BoundaryConditionProcessor.h deleted file mode 100644 index 9f8bb1522..000000000 --- a/source/VirtualFluidsCore/BoundaryCondition/BoundaryConditionProcessor.h +++ /dev/null @@ -1,26 +0,0 @@ -#ifndef BoundaryConditionProcessor_h__ -#define BoundaryConditionProcessor_h__ - -#include "BoundaryCondition.h" -#include <vector> - -class BoundaryConditionProcessor; -typedef boost::shared_ptr<BoundaryConditionProcessor> BoundaryConditionProcessorPtr; - -class BoundaryConditionProcessor -{ -public: - BoundaryConditionProcessor(); - virtual ~BoundaryConditionProcessor() {} - - void addBC(BoundaryConditionPtr bc); - BoundaryConditionPtr getBC(BoundaryCondition::Type type); - void applyPreCollisionBC(); - void applyPostCollisionBC(); -protected: -private: - std::vector<BoundaryConditionPtr> preBC; - std::vector<BoundaryConditionPtr> postBC; -}; - -#endif // BoundaryConditionProcessor_h__ diff --git a/source/VirtualFluidsCore/BoundaryCondition/D3Q27ETBCProcessor.cpp b/source/VirtualFluidsCore/BoundaryCondition/D3Q27ETBCProcessor.cpp index ff567490c..8992cdecd 100644 --- a/source/VirtualFluidsCore/BoundaryCondition/D3Q27ETBCProcessor.cpp +++ b/source/VirtualFluidsCore/BoundaryCondition/D3Q27ETBCProcessor.cpp @@ -8,18 +8,18 @@ D3Q27ETBCProcessor::D3Q27ETBCProcessor() ////////////////////////////////////////////////////////////////////////// D3Q27ETBCProcessor::D3Q27ETBCProcessor(LBMKernel3DPtr kernel) { - distributions = boost::dynamic_pointer_cast<EsoTwist3D>(kernel->getDataSet()->getFdistributions()); - ghostLayerWitdh = kernel->getGhostLayerWidth(); - collFactor = kernel->getCollisionFactor(); + DistributionArray3DPtr distributions = boost::dynamic_pointer_cast<EsoTwist3D>(kernel->getDataSet()->getFdistributions()); + //ghostLayerWitdh = kernel->getGhostLayerWidth(); + //collFactor = kernel->getCollisionFactor(); bcArray.resize( distributions->getNX1(), distributions->getNX2(), distributions->getNX3(), BCArray3D<D3Q27BoundaryCondition>::FLUID); - this->compressible = kernel->getCompressible(); - - minX1 = 0; - minX2 = 0; - minX3 = 0; - maxX1 = (int)bcArray.getNX1(); - maxX2 = (int)bcArray.getNX2(); - maxX3 = (int)bcArray.getNX3(); + //this->compressible = kernel->getCompressible(); + + //minX1 = 0; + //minX2 = 0; + //minX3 = 0; + //maxX1 = (int)bcArray.getNX1(); + //maxX2 = (int)bcArray.getNX2(); + //maxX3 = (int)bcArray.getNX3(); } ////////////////////////////////////////////////////////////////////////// D3Q27ETBCProcessor::~D3Q27ETBCProcessor() @@ -27,44 +27,9 @@ D3Q27ETBCProcessor::~D3Q27ETBCProcessor() } ////////////////////////////////////////////////////////////////////////// -//void D3Q27ETBCProcessor::init(LBMKernel3DPtr kernel) -//{ -// distributions = boost::dynamic_pointer_cast<EsoTwist3D>(kernel->getDataSet()->getFdistributions()); -// ghostLayerWitdh = kernel->getGhostLayerWidth(); -// collFactor = kernel->getCollisionFactor(); -// bcArray.resize( distributions->getNX1(), distributions->getNX2(), distributions->getNX3(), BCArray3D<D3Q27BoundaryCondition>::FLUID); -// this->compressible = kernel->getCompressible(); -// CalcFeqForDirFct calcFeqsForDirFct = NULL; -// CalcMacrosFct calcMacrosFct = NULL; -// CalcFeqFct calcFeqFct = NULL; -// -// if(this->compressible) -// { -// rho0 = REAL_CAST(0.0); -// calcFeqsForDirFct = &D3Q27System::getCompFeqForDirection; -// calcMacrosFct = &D3Q27System::calcCompMacroscopicValues; -// calcFeqFct = &D3Q27System::calcCompFeq; -// } -// else -// { -// rho0 = REAL_CAST(1.0); -// calcFeqsForDirFct = &D3Q27System::getIncompFeqForDirection; -// calcMacrosFct = &D3Q27System::calcIncompMacroscopicValues; -// calcFeqFct = &D3Q27System::calcIncompFeq; -// } -// -// minX1 = 0; -// minX2 = 0; -// minX3 = 0; -// maxX1 = (int)bcArray.getNX1(); -// maxX2 = (int)bcArray.getNX2(); -// maxX3 = (int)bcArray.getNX3(); -//} -////////////////////////////////////////////////////////////////////////// BCProcessorPtr D3Q27ETBCProcessor::clone(LBMKernel3DPtr kernel) { BCProcessorPtr bcProcessor(new D3Q27ETBCProcessor(kernel)); - //boost::dynamic_pointer_cast<D3Q27ETBCProcessor>(bcProcessor)->init(kernel); BoundaryConditionPtr velocity = this->getBC(BoundaryCondition::Velocity); BoundaryConditionPtr density = this->getBC(BoundaryCondition::Density); BoundaryConditionPtr noSlip = this->getBC(BoundaryCondition::NoSlip); @@ -81,484 +46,6 @@ BCArray3D<D3Q27BoundaryCondition>& D3Q27ETBCProcessor::getBCArray() return this->bcArray; } ////////////////////////////////////////////////////////////////////////// -void D3Q27ETBCProcessor::applyBC() -{ - typedef void (*CalcMacrosFct)(const LBMReal* const& /*f[27]*/, LBMReal& /*rho*/, LBMReal& /*vx1*/, LBMReal& /*vx2*/, LBMReal& /*vx3*/); - typedef LBMReal (*CalcFeqForDirFct)(const int& /*direction*/,const LBMReal& /*(d)rho*/,const LBMReal& /*vx1*/,const LBMReal& /*vx2*/,const LBMReal& /*vx3*/); - typedef void (*CalcFeqFct)(LBMReal* const& /*feq/*[27]*/,const LBMReal& /*rho*/,const LBMReal& /*vx1*/,const LBMReal& /*vx2*/,const LBMReal& /*vx3*/); - - CalcFeqForDirFct calcFeqsForDirFct = NULL; - CalcMacrosFct calcMacrosFct = NULL; - CalcFeqFct calcFeqFct = NULL; - - if(this->compressible) - { - rho0 = REAL_CAST(0.0); - calcFeqsForDirFct = &D3Q27System::getCompFeqForDirection; - calcMacrosFct = &D3Q27System::calcCompMacroscopicValues; - calcFeqFct = &D3Q27System::calcCompFeq; - } - else - { - rho0 = REAL_CAST(1.0); - calcFeqsForDirFct = &D3Q27System::getIncompFeqForDirection; - calcMacrosFct = &D3Q27System::calcIncompMacroscopicValues; - calcFeqFct = &D3Q27System::calcIncompFeq; - } - - - - ////////////////////////////////////////////////////////////////////////// - //EsoTwist3DPtr distributionsTemp(new D3Q27EsoTwist3DSplittedVector(maxX1, maxX2, maxX2, -999.0)); - //printf("create new distributionsTemp\n"); - ////////////////////////////////////////////////////////////////////////// - - for(int x3 = minX3; x3 < maxX3; x3++) - { - for(int x2 = minX2; x2 < maxX2; x2++) - { - for(int x1 = minX1; x1 < maxX1; x1++) - { - if(!bcArray.isSolid(x1,x2,x3) && !bcArray.isUndefined(x1,x2,x3)) - { - if( (bcPtr=bcArray.getBC(x1,x2,x3)) != NULL) - { - //velocity boundary condition - if(bcPtr->hasVelocityBoundary()) - { - //distributions->getDistributionInv(f,x1,x2,x3); - //int nx1 = x1; - //int nx2 = x2; - //int nx3 = x3; - - //rho = rho0; - //if(compressible) rho += D3Q27System::getDensity(f); - - //for(int fdir=D3Q27System::FSTARTDIR; fdir<=D3Q27System::FENDDIR; fdir++) - //{ - // if(bcPtr->hasVelocityBoundaryFlag(fdir)) - // { - // int invDir = D3Q27System::INVDIR[fdir]; - // LBMReal ftemp = distributions->getDistributionInvForDirection(x1+D3Q27System::DX1[invDir], x2+D3Q27System::DX2[invDir], x3+D3Q27System::DX3[invDir], fdir); - // ftemp += rho * bcPtr->getBoundaryVelocity(fdir); - // distributions->setDistributionForDirection(ftemp, x1+D3Q27System::DX1[invDir], x2+D3Q27System::DX2[invDir], x3+D3Q27System::DX3[invDir], fdir); - // } - //} - distributions->getDistributionInv(ftemp,x1,x2,x3); - calcMacrosFct(ftemp,rho,vx1,vx2,vx3); - calcFeqFct(feq,rho,vx1,vx2,vx3); - - for(int fdir=D3Q27System::FSTARTDIR; fdir<=D3Q27System::FENDDIR; fdir++) - { - if(bcPtr->hasVelocityBoundaryFlag(fdir)) - { - const int invDir = D3Q27System::INVDIR[fdir]; - LBMReal q = 1.0;//bcPtr->getQ(invDir);// m+m q=0 stabiler - LBMReal velocity = bcPtr->getBoundaryVelocity(invDir); - //LBMReal fReturn = ((1.0-q)/(1.0+q))*((ftemp[invDir]-feq[invDir]*collFactor)/(1.0-collFactor))+((q*(ftemp[invDir]+ftemp[fdir])-velocity)/(1.0+q)); - LBMReal fReturn = ((1.0-q)/(1.0+q))*((ftemp[invDir]-feq[invDir])/(1.0-collFactor)+feq[invDir])+((q*(ftemp[invDir]+ftemp[fdir])-velocity)/(1.0+q)); - distributions->setDistributionForDirection(fReturn, x1+D3Q27System::DX1[invDir], x2+D3Q27System::DX2[invDir], x3+D3Q27System::DX3[invDir], fdir); - } - } - } -///////////////////velocity bc for non reflecting pressure bc - //if (bcPtr->hasVelocityBoundary()) - //{ - // distributions->getDistributionInv(ftemp, x1, x2, x3); - // calcMacrosFct(ftemp, rho, vx1, vx2, vx3); - // calcFeqFct(feq, rho, vx1, vx2, vx3); - - // for (int fdir = D3Q27System::FSTARTDIR; fdir <= D3Q27System::FENDDIR; fdir++) - // { - // if (bcPtr->hasVelocityBoundaryFlag(fdir)) - // { - // const int invDir = D3Q27System::INVDIR[fdir]; - // LBMReal q = 1.0;//bcPtr->getQ(invDir);// m+m q=0 stabiler - // LBMReal velocity = bcPtr->getBoundaryVelocity(invDir); - // //LBMReal fReturn = ((1.0-q)/(1.0+q))*((ftemp[invDir]-feq[invDir]*collFactor)/(1.0-collFactor))+((q*(ftemp[invDir]+ftemp[fdir])-velocity)/(1.0+q)); - // LBMReal fReturn = ((1.0 - q) / (1.0 + q))*((ftemp[invDir] - feq[invDir]) / (1.0 - collFactor) + feq[invDir]) + ((q*(ftemp[invDir] + ftemp[fdir]) - velocity) / (1.0 + q))-rho*D3Q27System::WEIGTH[invDir]; - // distributions->setDistributionForDirection(fReturn, x1 + D3Q27System::DX1[invDir], x2 + D3Q27System::DX2[invDir], x3 + D3Q27System::DX3[invDir], fdir); - // } - // } - //} - //density boundary condition - else if(bcPtr->hasDensityBoundary()) - { - distributions->getDistributionInv(f,x1,x2,x3); - int nx1 = x1; - int nx2 = x2; - int nx3 = x3; - int direction = -1; - - //flag points in direction of fluid - if ( bcPtr->hasDensityBoundaryFlag(D3Q27System::E) ) {nx1 -= 1; direction=D3Q27System::E; } - else if( bcPtr->hasDensityBoundaryFlag(D3Q27System::W) ) {nx1 += 1; direction=D3Q27System::W; } - else if( bcPtr->hasDensityBoundaryFlag(D3Q27System::N) ) {nx2 -= 1; direction=D3Q27System::N; } - else if( bcPtr->hasDensityBoundaryFlag(D3Q27System::S) ) {nx2 += 1; direction=D3Q27System::S; } - else if( bcPtr->hasDensityBoundaryFlag(D3Q27System::T) ) {nx3 -= 1; direction=D3Q27System::T; } - else if( bcPtr->hasDensityBoundaryFlag(D3Q27System::B) ) {nx3 += 1; direction=D3Q27System::B; } - else UB_THROW( UbException(UB_EXARGS,"Danger...no orthogonal BC-Flag on density boundary...") ); - - #ifdef _DEBUG - if( nx1<0 || nx1>maxX1 ) UB_THROW( UbException(UB_EXARGS,"nx1<0 || nx1>=lengthX1") ); - if( nx2<0 || nx2>maxX2 ) UB_THROW( UbException(UB_EXARGS,"nx2<0 || nx2>=lengthX2") ); - if( nx3<0 || nx3>maxX3 ) UB_THROW( UbException(UB_EXARGS,"nx3<0 || nx3>=lengthX3") ); - #endif - - //if(bcPtr->getDensitySecondaryOption(direction)!=0) - // UB_THROW( UbException(UB_EXARGS,"invalid bc->getDensitySecondaryOption()="+UbSystem::toString(bcPtr->getDensitySecondaryOption(direction))) ); - - calcMacrosFct(f,rho,vx1,vx2,vx3); - LBMReal rhoBC = bcPtr->getBoundaryDensity(); - for(int fdir=D3Q27System::STARTF; fdir<=D3Q27System::ENDF; fdir++) - { - if(bcPtr->hasDensityBoundaryFlag(fdir)) - { - short densitySO = bcPtr->getDensitySecondaryOption(fdir); - if ( densitySO == 0) - { - // Martins NEQ ADDON - ////original: 15.2.2013: - LBMReal ftemp = calcFeqsForDirFct(fdir,rho,vx1,vx2,vx3); - ftemp = calcFeqsForDirFct(fdir,rhoBC,vx1,vx2,vx3)+f[fdir]-ftemp; - distributions->setDistributionForDirection(ftemp, nx1, nx2, nx3, fdir); - } - else if(densitySO == 1) - { - //Ehsan: 15.2.2013: - LBMReal ftemp = calcFeqsForDirFct(fdir,rhoBC,vx1,vx2,vx3); - distributions->setDistributionForDirection(ftemp, nx1, nx2, nx3, fdir); - } - } - } - } -///////////////non reflecting pressure baundary condition - //else if(bcPtr->hasDensityBoundary()) - // { -// distributions->swap(); -// distributions->getDistribution(f, x1, x2, x3); -// int nx1 = x1; -// int nx2 = x2; -// int nx3 = x3; -// int direction = -1; -// -// //flag points in direction of fluid -// if (bcPtr->hasDensityBoundaryFlag(D3Q27System::E)) { nx1 += 1; direction = D3Q27System::E; } -// else if (bcPtr->hasDensityBoundaryFlag(D3Q27System::W)) { nx1 -= 1; direction = D3Q27System::W; } -// else if (bcPtr->hasDensityBoundaryFlag(D3Q27System::N)) { nx2 += 1; direction = D3Q27System::N; } -// else if (bcPtr->hasDensityBoundaryFlag(D3Q27System::S)) { nx2 -= 1; direction = D3Q27System::S; } -// else if (bcPtr->hasDensityBoundaryFlag(D3Q27System::T)) { nx3 += 1; direction = D3Q27System::T; } -// else if (bcPtr->hasDensityBoundaryFlag(D3Q27System::B)) { nx3 -= 1; direction = D3Q27System::B; } -// else UB_THROW(UbException(UB_EXARGS, "Danger...no orthogonal BC-Flag on density boundary...")); -// -//#ifdef _DEBUG -// if (nx1<0 || nx1>maxX1) UB_THROW(UbException(UB_EXARGS, "nx1<0 || nx1>=lengthX1")); -// if (nx2<0 || nx2>maxX2) UB_THROW(UbException(UB_EXARGS, "nx2<0 || nx2>=lengthX2")); -// if (nx3<0 || nx3>maxX3) UB_THROW(UbException(UB_EXARGS, "nx3<0 || nx3>=lengthX3")); -//#endif -// distributions->getDistribution(ftemp, nx1, nx2, nx3); -// calcMacrosFct(f, rho, vx1, vx2, vx3); -// -// distributions->swap(); -// -// double dim = 0.01; -// -// f[D3Q27System::W] = ftemp[D3Q27System::W] * (one_over_sqrt3 - vx1) + (1.0 - one_over_sqrt3 + vx1)*f[D3Q27System::W] - rho*dim*D3Q27System::WEIGTH[D3Q27System::W]; -// f[D3Q27System::NW] = ftemp[D3Q27System::NW] * (one_over_sqrt3 - vx1) + (1.0 - one_over_sqrt3 + vx1)*f[D3Q27System::NW] - rho*dim*D3Q27System::WEIGTH[D3Q27System::NW]; -// f[D3Q27System::SW] = ftemp[D3Q27System::SW] * (one_over_sqrt3 - vx1) + (1.0 - one_over_sqrt3 + vx1)*f[D3Q27System::SW] - rho*dim*D3Q27System::WEIGTH[D3Q27System::SW]; -// f[D3Q27System::TW] = ftemp[D3Q27System::TW] * (one_over_sqrt3 - vx1) + (1.0 - one_over_sqrt3 + vx1)*f[D3Q27System::TW] - rho*dim*D3Q27System::WEIGTH[D3Q27System::TW]; -// f[D3Q27System::BW] = ftemp[D3Q27System::BW] * (one_over_sqrt3 - vx1) + (1.0 - one_over_sqrt3 + vx1)*f[D3Q27System::BW] - rho*dim*D3Q27System::WEIGTH[D3Q27System::BW]; -// f[D3Q27System::TNW] = ftemp[D3Q27System::TNW] * (one_over_sqrt3 - vx1) + (1.0 - one_over_sqrt3 + vx1)*f[D3Q27System::TNW] - rho*dim*D3Q27System::WEIGTH[D3Q27System::TNW]; -// f[D3Q27System::TSW] = ftemp[D3Q27System::TSW] * (one_over_sqrt3 - vx1) + (1.0 - one_over_sqrt3 + vx1)*f[D3Q27System::TSW] - rho*dim*D3Q27System::WEIGTH[D3Q27System::TSW]; -// f[D3Q27System::BNW] = ftemp[D3Q27System::BNW] * (one_over_sqrt3 - vx1) + (1.0 - one_over_sqrt3 + vx1)*f[D3Q27System::BNW] - rho*dim*D3Q27System::WEIGTH[D3Q27System::BNW]; -// f[D3Q27System::BSW] = ftemp[D3Q27System::BSW] * (one_over_sqrt3 - vx1) + (1.0 - one_over_sqrt3 + vx1)*f[D3Q27System::BSW] - rho*dim*D3Q27System::WEIGTH[D3Q27System::BSW]; -// -// distributions->setDistributionInvForDirection(f[D3Q27System::W], x1, x2, x3, D3Q27System::W); -// distributions->setDistributionInvForDirection(f[D3Q27System::NW], x1, x2, x3, D3Q27System::NW); -// distributions->setDistributionInvForDirection(f[D3Q27System::SW], x1, x2, x3, D3Q27System::SW); -// distributions->setDistributionInvForDirection(f[D3Q27System::TW], x1, x2, x3, D3Q27System::TW); -// distributions->setDistributionInvForDirection(f[D3Q27System::BW], x1, x2, x3, D3Q27System::BW); -// distributions->setDistributionInvForDirection(f[D3Q27System::TNW], x1, x2, x3, D3Q27System::TNW); -// distributions->setDistributionInvForDirection(f[D3Q27System::TSW], x1, x2, x3, D3Q27System::TSW); -// distributions->setDistributionInvForDirection(f[D3Q27System::BNW], x1, x2, x3, D3Q27System::BNW); -// distributions->setDistributionInvForDirection(f[D3Q27System::BSW], x1, x2, x3, D3Q27System::BSW); - - //} - else if(bcPtr->hasNoSlipBoundary()) - { - distributions->getDistributionInv(ftemp,x1,x2,x3); - calcMacrosFct(ftemp,rho,vx1,vx2,vx3); - calcFeqFct(feq,rho,vx1,vx2,vx3); - - for(int fdir=D3Q27System::FSTARTDIR; fdir<=D3Q27System::FENDDIR; fdir++) - { - if(bcPtr->hasNoSlipBoundaryFlag(fdir)) - { - switch(bcPtr->getNoSlipSecondaryOption(fdir)) - { - case 0: - //simple bounce back - do nothing - break; - case 1: - { - //quadratic bounce back - const int invDir = D3Q27System::INVDIR[fdir]; - LBMReal q = bcPtr->getQ(invDir); - //LBMReal fReturn = ((1.0-q)/(1.0+q))*((ftemp[invDir]-feq[invDir]*collFactor)/(1.0-collFactor))+((q/(1.0+q))*(ftemp[invDir]+ftemp[fdir])); - LBMReal fReturn = ((1.0-q)/(1.0+q))*((ftemp[invDir]-feq[invDir])/(1.0-collFactor)+feq[invDir])+((q/(1.0+q))*(ftemp[invDir]+ftemp[fdir])); - distributions->setDistributionForDirection(fReturn, x1+D3Q27System::DX1[invDir], x2+D3Q27System::DX2[invDir], x3+D3Q27System::DX3[invDir], fdir); - } - break; - ////////////////////////////////////////////////////////////////////////// - ////HighViscosityNoSlipBoundaryCondition - //case 1: - //{ - // //quadratic bounce back - // const int invDir = D3Q27System::INVDIR[fdir]; - // LBMReal q = bcPtr->getQ(invDir); - // LBMReal fReturn = ((1.0-q)*ftemp[invDir]+q*((ftemp[invDir]+ftemp[fdir])*(1-collFactor)+collFactor*(feq[invDir]+feq[fdir])))/(1.0+q); - // distributions->setDistributionForDirection(fReturn, x1+D3Q27System::DX1[invDir], x2+D3Q27System::DX2[invDir], x3+D3Q27System::DX3[invDir], fdir); - //} - //break; - ////////////////////////////////////////////////////////////////////////// - //case 2: - // { - // //quadratic bounce back 2nd choice - // const int invDir = D3Q27System::INVDIR[fdir]; - // LBMReal q = bcPtr->getQ(invDir); - // LBMReal fReturn = ((1.0-q)/(1.0+q))*0.5*(ftemp[invDir]-ftemp[fdir]+(ftemp[invDir]+ftemp[fdir]-collFactor*(feq[fdir]+feq[invDir]))/(1.0-collFactor)); - // distributions->setDistributionForDirection(fReturn, x1+D3Q27System::DX1[invDir], x2+D3Q27System::DX2[invDir], x3+D3Q27System::DX3[invDir], fdir); - // } - // break; - case 2: - { - //quadratic bounce back with for thin walls - const int invDir = D3Q27System::INVDIR[fdir]; - LBMReal q = bcPtr->getQ(invDir); - //LBMReal fReturn = ((1.0-q)/(1.0+q))*((ftemp[invDir]-feq[invDir])/(1.0-collFactor)+feq[invDir])+((q/(1.0+q))*(ftemp[invDir]+ftemp[fdir])); - LBMReal fReturn = ((1.0-q)/(1.0+q))*0.5*(ftemp[invDir]-ftemp[fdir]+(ftemp[invDir]+ftemp[fdir]-collFactor*(feq[fdir]+feq[invDir]))/(1.0-collFactor)); - //printf("set distributionsTemp - start\n"); - distributionsTemp->setDistributionForDirection(fReturn, x1+D3Q27System::DX1[invDir], x2+D3Q27System::DX2[invDir], x3+D3Q27System::DX3[invDir], fdir); - //printf("set distributionsTemp - end\n"); - } - break; - } - } - } - } - else if(bcPtr->hasSlipBoundary()) - { - distributions->getDistributionInv(ftemp,x1,x2,x3); - calcMacrosFct(ftemp,rho,vx1,vx2,vx3); - calcFeqFct(feq,rho,vx1,vx2,vx3); - - UbTupleFloat3 normale = bcPtr->getNormalVector(); - LBMReal amp = vx1*val<1>(normale)+vx2*val<2>(normale)+vx3*val<3>(normale); - - vx1 = vx1 - amp * val<1>(normale); //normale zeigt von struktur weg! - vx2 = vx2 - amp * val<2>(normale); //normale zeigt von struktur weg! - vx3 = vx3 - amp * val<3>(normale); //normale zeigt von struktur weg! - - for(int fdir=D3Q27System::FSTARTDIR; fdir<=D3Q27System::FENDDIR; fdir++) - { - if(bcPtr->hasSlipBoundaryFlag(fdir)) - { - - - switch( bcPtr->getSlipSecondaryOption(fdir) ) - { - case 0: - //doesn't work - //simple bounce back - do nothing - break; - case 1: - { - //quadratic bounce back - const int invDir = D3Q27System::INVDIR[fdir]; - LBMReal q = 1.0;//bcPtr->getQ(invDir);// m+m q=0 stabiler - //vx3=0; - LBMReal velocity = 0.0; - switch(invDir) - { - case D3Q27System::E : velocity = ( UbMath::c4o9*(+vx1) ); break; //(2/cs^2)(=6)*rho_0(=1 bei imkompr)*wi*u*ei mit cs=1/sqrt(3) - case D3Q27System::W : velocity = ( UbMath::c4o9*(-vx1) ); break; //z.B. aus paper manfred MRT LB models in three dimensions (2002) - case D3Q27System::N : velocity = ( UbMath::c4o9*(+vx2) ); break; - case D3Q27System::S : velocity = ( UbMath::c4o9*(-vx2) ); break; - case D3Q27System::T : velocity = ( UbMath::c4o9*(+vx3) ); break; - case D3Q27System::B : velocity = ( UbMath::c4o9*(-vx3) ); break; - case D3Q27System::NE: velocity = ( UbMath::c1o9*(+vx1+vx2 ) ); break; - case D3Q27System::SW: velocity = ( UbMath::c1o9*(-vx1-vx2 ) ); break; - case D3Q27System::SE: velocity = ( UbMath::c1o9*(+vx1-vx2 ) ); break; - case D3Q27System::NW: velocity = ( UbMath::c1o9*(-vx1+vx2 ) ); break; - case D3Q27System::TE: velocity = ( UbMath::c1o9*(+vx1 +vx3) ); break; - case D3Q27System::BW: velocity = ( UbMath::c1o9*(-vx1 -vx3) ); break; - case D3Q27System::BE: velocity = ( UbMath::c1o9*(+vx1 -vx3) ); break; - case D3Q27System::TW: velocity = ( UbMath::c1o9*(-vx1 +vx3) ); break; - case D3Q27System::TN: velocity = ( UbMath::c1o9*( +vx2+vx3) ); break; - case D3Q27System::BS: velocity = ( UbMath::c1o9*( -vx2-vx3) ); break; - case D3Q27System::BN: velocity = ( UbMath::c1o9*( +vx2-vx3) ); break; - case D3Q27System::TS: velocity = ( UbMath::c1o9*( -vx2+vx3) ); break; - case D3Q27System::TNE: velocity = ( UbMath::c1o36*(+vx1+vx2+vx3) ); break; - case D3Q27System::BSW: velocity = ( UbMath::c1o36*(-vx1-vx2-vx3) ); break; - case D3Q27System::BNE: velocity = ( UbMath::c1o36*(+vx1+vx2-vx3) ); break; - case D3Q27System::TSW: velocity = ( UbMath::c1o36*(-vx1-vx2+vx3) ); break; - case D3Q27System::TSE: velocity = ( UbMath::c1o36*(+vx1-vx2+vx3) ); break; - case D3Q27System::BNW: velocity = ( UbMath::c1o36*(-vx1+vx2-vx3) ); break; - case D3Q27System::BSE: velocity = ( UbMath::c1o36*(+vx1-vx2-vx3) ); break; - case D3Q27System::TNW: velocity = ( UbMath::c1o36*(-vx1+vx2+vx3) ); break; - default: throw UbException(UB_EXARGS,"unknown error"); - } - //LBMReal fReturn = ((1.0-q)/(1.0+q))*((ftemp[invDir]-feq[invDir]*collFactor)/(1.0-collFactor))+((q*(ftemp[invDir]+ftemp[fdir])-velocity)/(1.0+q)); - LBMReal fReturn = ((1.0-q)/(1.0+q))*((ftemp[invDir]-feq[invDir])/(1.0-collFactor)+feq[invDir])+((q*(ftemp[invDir]+ftemp[fdir])-velocity)/(1.0+q)); - distributions->setDistributionForDirection(fReturn, x1+D3Q27System::DX1[invDir], x2+D3Q27System::DX2[invDir], x3+D3Q27System::DX3[invDir], fdir); - } - break; - default: - UB_THROW( UbException(UB_EXARGS,"unknown secondary option for SlipBC") ); - } - } - } - } - } - } - } - } - } - - -// for(int x3 = minX3; x3 < maxX3; x3++) -// { -// for(int x2 = minX2; x2 < maxX2; x2++) -// { -// for(int x1 = minX1; x1 < maxX1; x1++) -// { -// if(!bcArray.isSolid(x1,x2,x3) && !bcArray.isUndefined(x1,x2,x3)) -// { -// if( (bcPtr=bcArray.getBC(x1,x2,x3)) != NULL) -// { -// if(bcPtr->hasNoSlipBoundary()) -// { -// //distributions->getDistributionInv(ftemp,x1,x2,x3); -// //calcMacrosFct(ftemp,rho,vx1,vx2,vx3); -// //calcFeqFct(feq,rho,vx1,vx2,vx3); -// -// for(int fdir=D3Q27System::FSTARTDIR; fdir<=D3Q27System::FENDDIR; fdir++) -// { -// if(bcPtr->hasNoSlipBoundaryFlag(fdir)) -// { -// switch(bcPtr->getNoSlipSecondaryOption(fdir)) -// { -// case 2: -// { -// //quadratic bounce back with for thin walls -// const int invDir = D3Q27System::INVDIR[fdir]; -///* LBMReal q = bcPtr->getQ(invDir); -// LBMReal fReturn = ((1.0-q)/(1.0+q))*0.5*(ftemp[invDir]-ftemp[fdir]+(ftemp[invDir]+ftemp[fdir]-collFactor*(feq[fdir]+feq[invDir]))/(1.0-collFactor));*/ -// //printf("get distributionsTemp - start\n"); -// //printf("x1=%d,x2=%d,x3=%d,D3Q27System::DX1[invDir]=%d,D3Q27System::DX2[invDir]=%d,D3Q27System::DX3[invDir]%d\n",x1,x2,x3,D3Q27System::DX1[invDir],D3Q27System::DX2[invDir],D3Q27System::DX3[invDir]); -// ftemp[fdir] = distributionsTemp->getDistributionInvForDirection(x1+D3Q27System::DX1[invDir], x2+D3Q27System::DX2[invDir], x3+D3Q27System::DX3[invDir], fdir); -// //printf("get distributionsTemp - end\n"); -// distributions->setDistributionForDirection(ftemp[fdir], x1+D3Q27System::DX1[invDir], x2+D3Q27System::DX2[invDir], x3+D3Q27System::DX3[invDir], fdir); -// } -// break; -// } -// } -// } -// } -// } -// } -// } -// } -// } - -} -////////////////////////////////////////////////////////////////////////// -//void D3Q27ETBCProcessor::init() -//{ -// BoundaryConditionPtr velocity = getBC(BoundaryCondition::Velocity); -// BoundaryConditionPtr density = getBC(BoundaryCondition::Density); -// BoundaryConditionPtr noSlip = getBC(BoundaryCondition::NoSlip); -// BoundaryConditionPtr slip = getBC(BoundaryCondition::Slip); -// -// if (velocity) velocity->setCompressible(compressible); -// if (density) density->setCompressible(compressible); -// if (noSlip) noSlip->setCompressible(compressible); -// if (slip) slip->setCompressible(compressible); -// -// int vCount = 0; -// int dCount = 0; -// int nsCount = 0; -// int sCount = 0; -// -// for (int x3 = minX3; x3 < maxX3; x3++) -// { -// for (int x2 = minX2; x2 < maxX2; x2++) -// { -// for (int x1 = minX1; x1 < maxX1; x1++) -// { -// if (!bcArray.isSolid(x1, x2, x3) && !bcArray.isUndefined(x1, x2, x3)) -// { -// if ((bcPtr = bcArray.getBC(x1, x2, x3)) != NULL) -// { -// //velocity boundary condition -// if (bcPtr->hasVelocityBoundary() && velocity) -// { -// velocity->addNode(x1, x2, x3); -// velocity->addBcPointer(bcPtr); -// vCount++; -// } -// //density boundary condition -// else if (bcPtr->hasDensityBoundary() && density) -// { -// density->addNode(x1, x2, x3); -// density->addBcPointer(bcPtr); -// dCount++; -// } -// //no slip boundary condition -// else if (bcPtr->hasNoSlipBoundary() && noSlip) -// { -// noSlip->addNode(x1, x2, x3); -// noSlip->addBcPointer(bcPtr); -// nsCount++; -// } -// //slip boundary condition -// else if (bcPtr->hasSlipBoundary() && slip) -// { -// slip->addNode(x1, x2, x3); -// slip->addBcPointer(bcPtr); -// sCount++; -// } -// } -// } -// } -// } -// } -// if (vCount > 0) -// { -// velocity->addNnode(vCount); -// velocity->addDistributions(distributions); -// velocity->setCollFactor(collFactor); -// } -// if (dCount > 0) -// { -// density->addNnode(dCount); -// density->addDistributions(distributions); -// density->setCollFactor(collFactor); -// } -// if (nsCount > 0) -// { -// noSlip->addNnode(nsCount); -// noSlip->addDistributions(distributions); -// noSlip->setCollFactor(collFactor); -// } -// if (sCount > 0) -// { -// slip->addNnode(sCount); -// slip->addDistributions(distributions); -// slip->setCollFactor(collFactor); -// } -//} -////////////////////////////////////////////////////////////////////////// void D3Q27ETBCProcessor::addBC(BoundaryConditionPtr bc) { BoundaryCondition::Type type = bc->getType(); diff --git a/source/VirtualFluidsCore/BoundaryCondition/D3Q27ETBCProcessor.h b/source/VirtualFluidsCore/BoundaryCondition/D3Q27ETBCProcessor.h index 75e34f592..2599851ec 100644 --- a/source/VirtualFluidsCore/BoundaryCondition/D3Q27ETBCProcessor.h +++ b/source/VirtualFluidsCore/BoundaryCondition/D3Q27ETBCProcessor.h @@ -21,7 +21,7 @@ public: D3Q27ETBCProcessor(); D3Q27ETBCProcessor(LBMKernel3DPtr kernel); virtual ~D3Q27ETBCProcessor(); - virtual void applyBC(); + //virtual void applyBC(); virtual BCArray3D<D3Q27BoundaryCondition>& getBCArray(); virtual BCProcessorPtr clone(LBMKernel3DPtr kernel); @@ -34,20 +34,20 @@ protected: std::vector<BoundaryConditionPtr> preBC; std::vector<BoundaryConditionPtr> postBC; - int minX1, minX2, minX3, maxX1, maxX2, maxX3; - D3Q27BoundaryConditionPtr bcPtr; - LBMReal f[D3Q27System::ENDF+1]; - LBMReal ftemp[D3Q27System::ENDF+1]; - LBMReal feq[D3Q27System::ENDF+1]; - LBMReal rho, vx1, vx2, vx3, rho0; + //int minX1, minX2, minX3, maxX1, maxX2, maxX3; + //D3Q27BoundaryConditionPtr bcPtr; + //LBMReal f[D3Q27System::ENDF+1]; + //LBMReal ftemp[D3Q27System::ENDF+1]; + //LBMReal feq[D3Q27System::ENDF+1]; + //LBMReal rho, vx1, vx2, vx3, rho0; //void init(LBMKernel3DPtr kernel); BCArray3D<D3Q27BoundaryCondition> bcArray; - EsoTwist3DPtr distributions; - int ghostLayerWitdh; - LBMReal collFactor; - bool compressible; - EsoTwist3DPtr distributionsTemp; + //EsoTwist3DPtr distributions; + //int ghostLayerWitdh; + //LBMReal collFactor; + //bool compressible; + //EsoTwist3DPtr distributionsTemp; //std::vector <int> sizeVector; //std::vector <int> nodeVector; @@ -66,17 +66,17 @@ private: { ar & boost::serialization::base_object<BCProcessor>(*this); ar & bcArray; - ar & distributions; - ar & ghostLayerWitdh; - ar & collFactor; - ar & compressible; - ar & minX1; - ar & minX2; - ar & minX3; - ar & maxX1; - ar & maxX2; - ar & maxX3; - ar & distributionsTemp; + //ar & distributions; + //ar & ghostLayerWitdh; + //ar & collFactor; + //ar & compressible; + //ar & minX1; + //ar & minX2; + //ar & minX3; + //ar & maxX1; + //ar & maxX2; + //ar & maxX3; + //ar & distributionsTemp; ar & preBC; ar & postBC; } diff --git a/source/VirtualFluidsCore/BoundaryCondition/D3Q27ETForThinWallBCProcessor.cpp b/source/VirtualFluidsCore/BoundaryCondition/D3Q27ETForThinWallBCProcessor.cpp index f12da6366..e5639532f 100644 --- a/source/VirtualFluidsCore/BoundaryCondition/D3Q27ETForThinWallBCProcessor.cpp +++ b/source/VirtualFluidsCore/BoundaryCondition/D3Q27ETForThinWallBCProcessor.cpp @@ -1,6 +1,6 @@ #include "D3Q27ETForThinWallBCProcessor.h" -#include "D3Q27EsoTwist3DSplittedVector.h" +#include "ThinWallNoSlipBoundaryCondition.h" D3Q27ETForThinWallBCProcessor::D3Q27ETForThinWallBCProcessor() { @@ -9,7 +9,7 @@ D3Q27ETForThinWallBCProcessor::D3Q27ETForThinWallBCProcessor() ////////////////////////////////////////////////////////////////////////// D3Q27ETForThinWallBCProcessor::D3Q27ETForThinWallBCProcessor(LBMKernel3DPtr kernel) : D3Q27ETBCProcessor(kernel) { - //distributionsTemp = EsoTwist3DPtr(new D3Q27EsoTwist3DSplittedVector(maxX1, maxX2, maxX2, -999.0)); + } ////////////////////////////////////////////////////////////////////////// D3Q27ETForThinWallBCProcessor::~D3Q27ETForThinWallBCProcessor() @@ -20,52 +20,17 @@ D3Q27ETForThinWallBCProcessor::~D3Q27ETForThinWallBCProcessor() BCProcessorPtr D3Q27ETForThinWallBCProcessor::clone(LBMKernel3DPtr kernel) { BCProcessorPtr bcProcessor(new D3Q27ETForThinWallBCProcessor(kernel)); + BoundaryConditionPtr velocity = this->getBC(BoundaryCondition::Velocity); + BoundaryConditionPtr density = this->getBC(BoundaryCondition::Density); + BoundaryConditionPtr noSlip = this->getBC(BoundaryCondition::NoSlip); + BoundaryConditionPtr slip = this->getBC(BoundaryCondition::Slip); + if (velocity)bcProcessor->addBC(velocity->clone()); + if (density) bcProcessor->addBC(density->clone()); + if (noSlip) bcProcessor->addBC(noSlip->clone()); + if (slip) bcProcessor->addBC(slip->clone()); return bcProcessor; } ////////////////////////////////////////////////////////////////////////// -void D3Q27ETForThinWallBCProcessor::applyBC() -{ - D3Q27ETBCProcessor::applyBC(); - - LBMReal fReturn; - - for (int x3 = minX3; x3 < maxX3; x3++) - { - for (int x2 = minX2; x2 < maxX2; x2++) - { - for (int x1 = minX1; x1 < maxX1; x1++) - { - if (!bcArray.isSolid(x1, x2, x3) && !bcArray.isUndefined(x1, x2, x3)) - { - if ((bcPtr=bcArray.getBC(x1, x2, x3)) != NULL) - { - if (bcPtr->hasNoSlipBoundary()) - { - for (int fdir=D3Q27System::FSTARTDIR; fdir<=D3Q27System::FENDDIR; fdir++) - { - if (bcPtr->hasNoSlipBoundaryFlag(fdir)) - { - switch (bcPtr->getNoSlipSecondaryOption(fdir)) - { - case 2: - { - //quadratic bounce back with for thin walls - const int invDir = D3Q27System::INVDIR[fdir]; - fReturn = distributionsTemp->getDistributionInvForDirection(x1+D3Q27System::DX1[invDir], x2+D3Q27System::DX2[invDir], x3+D3Q27System::DX3[invDir], fdir); - distributions->setDistributionForDirection(fReturn, x1+D3Q27System::DX1[invDir], x2+D3Q27System::DX2[invDir], x3+D3Q27System::DX3[invDir], fdir); - } - break; - } - } - } - } - } - } - } - } - } -} -////////////////////////////////////////////////////////////////////////// void D3Q27ETForThinWallBCProcessor::applyPostCollisionBC() { D3Q27ETBCProcessor::applyPostCollisionBC(); @@ -74,7 +39,9 @@ void D3Q27ETForThinWallBCProcessor::applyPostCollisionBC() { if (bc->getType() == BoundaryCondition::NoSlip) { + boost::dynamic_pointer_cast<ThinWallNoSlipBoundaryCondition>(bc)->setPass(2); bc->apply(); + boost::dynamic_pointer_cast<ThinWallNoSlipBoundaryCondition>(bc)->setPass(1); } } } diff --git a/source/VirtualFluidsCore/BoundaryCondition/D3Q27ETForThinWallBCProcessor.h b/source/VirtualFluidsCore/BoundaryCondition/D3Q27ETForThinWallBCProcessor.h index 7acafa538..e0ff10c51 100644 --- a/source/VirtualFluidsCore/BoundaryCondition/D3Q27ETForThinWallBCProcessor.h +++ b/source/VirtualFluidsCore/BoundaryCondition/D3Q27ETForThinWallBCProcessor.h @@ -19,7 +19,6 @@ public: D3Q27ETForThinWallBCProcessor(); D3Q27ETForThinWallBCProcessor(LBMKernel3DPtr kernel); virtual ~D3Q27ETForThinWallBCProcessor(); - virtual void applyBC(); virtual BCProcessorPtr clone(LBMKernel3DPtr kernel); void applyPostCollisionBC(); protected: diff --git a/source/VirtualFluidsCore/BoundaryCondition/EqDensityBoundaryCondition.cpp b/source/VirtualFluidsCore/BoundaryCondition/EqDensityBoundaryCondition.cpp index d50f776c0..9f848d654 100644 --- a/source/VirtualFluidsCore/BoundaryCondition/EqDensityBoundaryCondition.cpp +++ b/source/VirtualFluidsCore/BoundaryCondition/EqDensityBoundaryCondition.cpp @@ -17,6 +17,11 @@ BoundaryConditionPtr EqDensityBoundaryCondition::clone() return bc; } ////////////////////////////////////////////////////////////////////////// +void EqDensityBoundaryCondition::addDistributions(DistributionArray3DPtr distributions) +{ + this->distributions = distributions; +} +////////////////////////////////////////////////////////////////////////// void EqDensityBoundaryCondition::applyBC() { LBMReal f[D3Q27System::ENDF+1]; diff --git a/source/VirtualFluidsCore/BoundaryCondition/EqDensityBoundaryCondition.h b/source/VirtualFluidsCore/BoundaryCondition/EqDensityBoundaryCondition.h index 65aaf8a08..ab1f7428f 100644 --- a/source/VirtualFluidsCore/BoundaryCondition/EqDensityBoundaryCondition.h +++ b/source/VirtualFluidsCore/BoundaryCondition/EqDensityBoundaryCondition.h @@ -12,6 +12,7 @@ public: EqDensityBoundaryCondition(); ~EqDensityBoundaryCondition(); BoundaryConditionPtr clone(); + void addDistributions(DistributionArray3DPtr distributions); protected: void applyBC(); private: diff --git a/source/VirtualFluidsCore/BoundaryCondition/HighViscosityNoSlipBoundaryCondition.cpp b/source/VirtualFluidsCore/BoundaryCondition/HighViscosityNoSlipBoundaryCondition.cpp index ad0b265d2..481ac4fc6 100644 --- a/source/VirtualFluidsCore/BoundaryCondition/HighViscosityNoSlipBoundaryCondition.cpp +++ b/source/VirtualFluidsCore/BoundaryCondition/HighViscosityNoSlipBoundaryCondition.cpp @@ -16,6 +16,11 @@ BoundaryConditionPtr HighViscosityNoSlipBoundaryCondition::clone() return bc; } ////////////////////////////////////////////////////////////////////////// +void HighViscosityNoSlipBoundaryCondition::addDistributions(DistributionArray3DPtr distributions) +{ + this->distributions = distributions; +} +////////////////////////////////////////////////////////////////////////// void HighViscosityNoSlipBoundaryCondition::applyBC() { LBMReal f[D3Q27System::ENDF+1]; diff --git a/source/VirtualFluidsCore/BoundaryCondition/HighViscosityNoSlipBoundaryCondition.h b/source/VirtualFluidsCore/BoundaryCondition/HighViscosityNoSlipBoundaryCondition.h index 4df602a09..8cc765b4c 100644 --- a/source/VirtualFluidsCore/BoundaryCondition/HighViscosityNoSlipBoundaryCondition.h +++ b/source/VirtualFluidsCore/BoundaryCondition/HighViscosityNoSlipBoundaryCondition.h @@ -12,6 +12,7 @@ public: HighViscosityNoSlipBoundaryCondition(); ~HighViscosityNoSlipBoundaryCondition(); BoundaryConditionPtr clone(); + void addDistributions(DistributionArray3DPtr distributions); protected: void applyBC(); private: diff --git a/source/VirtualFluidsCore/BoundaryCondition/NoSlipBoundaryCondition.cpp b/source/VirtualFluidsCore/BoundaryCondition/NoSlipBoundaryCondition.cpp index 5634cfad2..04235077b 100644 --- a/source/VirtualFluidsCore/BoundaryCondition/NoSlipBoundaryCondition.cpp +++ b/source/VirtualFluidsCore/BoundaryCondition/NoSlipBoundaryCondition.cpp @@ -17,6 +17,11 @@ BoundaryConditionPtr NoSlipBoundaryCondition::clone() return bc; } ////////////////////////////////////////////////////////////////////////// +void NoSlipBoundaryCondition::addDistributions(DistributionArray3DPtr distributions) +{ + this->distributions = distributions; +} +////////////////////////////////////////////////////////////////////////// void NoSlipBoundaryCondition::applyBC() { LBMReal f[D3Q27System::ENDF+1]; diff --git a/source/VirtualFluidsCore/BoundaryCondition/NoSlipBoundaryCondition.h b/source/VirtualFluidsCore/BoundaryCondition/NoSlipBoundaryCondition.h index 299e6807e..46e7fe0f8 100644 --- a/source/VirtualFluidsCore/BoundaryCondition/NoSlipBoundaryCondition.h +++ b/source/VirtualFluidsCore/BoundaryCondition/NoSlipBoundaryCondition.h @@ -12,6 +12,7 @@ public: NoSlipBoundaryCondition(); virtual ~NoSlipBoundaryCondition(); BoundaryConditionPtr clone(); + void addDistributions(DistributionArray3DPtr distributions); protected: void applyBC(); private: diff --git a/source/VirtualFluidsCore/BoundaryCondition/NonEqDensityBoundaryCondition.cpp b/source/VirtualFluidsCore/BoundaryCondition/NonEqDensityBoundaryCondition.cpp index 529c2b9bc..55dfed9d7 100644 --- a/source/VirtualFluidsCore/BoundaryCondition/NonEqDensityBoundaryCondition.cpp +++ b/source/VirtualFluidsCore/BoundaryCondition/NonEqDensityBoundaryCondition.cpp @@ -17,6 +17,11 @@ BoundaryConditionPtr NonEqDensityBoundaryCondition::clone() return bc; } ////////////////////////////////////////////////////////////////////////// +void NonEqDensityBoundaryCondition::addDistributions(DistributionArray3DPtr distributions) +{ + this->distributions = distributions; +} +////////////////////////////////////////////////////////////////////////// void NonEqDensityBoundaryCondition::applyBC() { LBMReal f[D3Q27System::ENDF+1]; diff --git a/source/VirtualFluidsCore/BoundaryCondition/NonEqDensityBoundaryCondition.h b/source/VirtualFluidsCore/BoundaryCondition/NonEqDensityBoundaryCondition.h index 59c42480d..f4a7d991b 100644 --- a/source/VirtualFluidsCore/BoundaryCondition/NonEqDensityBoundaryCondition.h +++ b/source/VirtualFluidsCore/BoundaryCondition/NonEqDensityBoundaryCondition.h @@ -12,6 +12,7 @@ public: NonEqDensityBoundaryCondition(); ~NonEqDensityBoundaryCondition(); BoundaryConditionPtr clone(); + void addDistributions(DistributionArray3DPtr distributions); protected: void applyBC(); private: diff --git a/source/VirtualFluidsCore/BoundaryCondition/NonReflectingDensityBoundaryCondition.cpp b/source/VirtualFluidsCore/BoundaryCondition/NonReflectingDensityBoundaryCondition.cpp index ff39c5320..68c6e8024 100644 --- a/source/VirtualFluidsCore/BoundaryCondition/NonReflectingDensityBoundaryCondition.cpp +++ b/source/VirtualFluidsCore/BoundaryCondition/NonReflectingDensityBoundaryCondition.cpp @@ -18,6 +18,11 @@ BoundaryConditionPtr NonReflectingDensityBoundaryCondition::clone() return bc; } ////////////////////////////////////////////////////////////////////////// +void NonReflectingDensityBoundaryCondition::addDistributions(DistributionArray3DPtr distributions) +{ + this->distributions = distributions; +} +////////////////////////////////////////////////////////////////////////// void NonReflectingDensityBoundaryCondition::applyBC() { distributions->swap(); diff --git a/source/VirtualFluidsCore/BoundaryCondition/NonReflectingDensityBoundaryCondition.h b/source/VirtualFluidsCore/BoundaryCondition/NonReflectingDensityBoundaryCondition.h index 5fcf412d5..c4f0fd14e 100644 --- a/source/VirtualFluidsCore/BoundaryCondition/NonReflectingDensityBoundaryCondition.h +++ b/source/VirtualFluidsCore/BoundaryCondition/NonReflectingDensityBoundaryCondition.h @@ -12,6 +12,7 @@ public: NonReflectingDensityBoundaryCondition(); ~NonReflectingDensityBoundaryCondition(); BoundaryConditionPtr clone(); + void addDistributions(DistributionArray3DPtr distributions); protected: void applyBC(); private: diff --git a/source/VirtualFluidsCore/BoundaryCondition/NonReflectingVelocityBoundaryCondition.cpp b/source/VirtualFluidsCore/BoundaryCondition/NonReflectingVelocityBoundaryCondition.cpp index 734545f6f..2986bdea2 100644 --- a/source/VirtualFluidsCore/BoundaryCondition/NonReflectingVelocityBoundaryCondition.cpp +++ b/source/VirtualFluidsCore/BoundaryCondition/NonReflectingVelocityBoundaryCondition.cpp @@ -17,6 +17,11 @@ BoundaryConditionPtr NonReflectingVelocityBoundaryCondition::clone() return bc; } ////////////////////////////////////////////////////////////////////////// +void NonReflectingVelocityBoundaryCondition::addDistributions(DistributionArray3DPtr distributions) +{ + this->distributions = distributions; +} +////////////////////////////////////////////////////////////////////////// void NonReflectingVelocityBoundaryCondition::applyBC() { //velocity bc for non reflecting pressure bc diff --git a/source/VirtualFluidsCore/BoundaryCondition/NonReflectingVelocityBoundaryCondition.h b/source/VirtualFluidsCore/BoundaryCondition/NonReflectingVelocityBoundaryCondition.h index cfb00fa20..cbbd66108 100644 --- a/source/VirtualFluidsCore/BoundaryCondition/NonReflectingVelocityBoundaryCondition.h +++ b/source/VirtualFluidsCore/BoundaryCondition/NonReflectingVelocityBoundaryCondition.h @@ -18,6 +18,7 @@ public: NonReflectingVelocityBoundaryCondition(); ~NonReflectingVelocityBoundaryCondition(); BoundaryConditionPtr clone(); + void addDistributions(DistributionArray3DPtr distributions); protected: void applyBC(); private: diff --git a/source/VirtualFluidsCore/BoundaryCondition/SlipBoundaryCondition.cpp b/source/VirtualFluidsCore/BoundaryCondition/SlipBoundaryCondition.cpp index fb96544cc..3aed2a494 100644 --- a/source/VirtualFluidsCore/BoundaryCondition/SlipBoundaryCondition.cpp +++ b/source/VirtualFluidsCore/BoundaryCondition/SlipBoundaryCondition.cpp @@ -17,6 +17,11 @@ BoundaryConditionPtr SlipBoundaryCondition::clone() return bc; } ////////////////////////////////////////////////////////////////////////// +void SlipBoundaryCondition::addDistributions(DistributionArray3DPtr distributions) +{ + this->distributions = distributions; +} +////////////////////////////////////////////////////////////////////////// void SlipBoundaryCondition::applyBC() { LBMReal f[D3Q27System::ENDF+1]; diff --git a/source/VirtualFluidsCore/BoundaryCondition/SlipBoundaryCondition.h b/source/VirtualFluidsCore/BoundaryCondition/SlipBoundaryCondition.h index a2a5b333b..17f6524a9 100644 --- a/source/VirtualFluidsCore/BoundaryCondition/SlipBoundaryCondition.h +++ b/source/VirtualFluidsCore/BoundaryCondition/SlipBoundaryCondition.h @@ -12,6 +12,7 @@ public: SlipBoundaryCondition(); virtual ~SlipBoundaryCondition(); BoundaryConditionPtr clone(); + void addDistributions(DistributionArray3DPtr distributions); protected: void applyBC(); private: diff --git a/source/VirtualFluidsCore/BoundaryCondition/ThinWallNoSlipBoundaryCondition.cpp b/source/VirtualFluidsCore/BoundaryCondition/ThinWallNoSlipBoundaryCondition.cpp index 21e07bf8a..d1de85397 100644 --- a/source/VirtualFluidsCore/BoundaryCondition/ThinWallNoSlipBoundaryCondition.cpp +++ b/source/VirtualFluidsCore/BoundaryCondition/ThinWallNoSlipBoundaryCondition.cpp @@ -7,6 +7,7 @@ ThinWallNoSlipBoundaryCondition::ThinWallNoSlipBoundaryCondition() { BoundaryCondition::type = BoundaryCondition::NoSlip; BoundaryCondition::preCollision = false; + pass = 1; } ////////////////////////////////////////////////////////////////////////// ThinWallNoSlipBoundaryCondition::~ThinWallNoSlipBoundaryCondition() @@ -22,23 +23,43 @@ BoundaryConditionPtr ThinWallNoSlipBoundaryCondition::clone() ////////////////////////////////////////////////////////////////////////// void ThinWallNoSlipBoundaryCondition::applyBC() { + LBMReal f[D3Q27System::ENDF + 1]; + LBMReal feq[D3Q27System::ENDF + 1]; + distributions->getDistributionInv(f, x1, x2, x3); + LBMReal rho, vx1, vx2, vx3; + calcMacrosFct(f, rho, vx1, vx2, vx3); + calcFeqFct(feq, rho, vx1, vx2, vx3); + LBMReal fReturn; for (int fdir = D3Q27System::FSTARTDIR; fdir<=D3Q27System::FENDDIR; fdir++) { if (bcPtr->hasNoSlipBoundaryFlag(fdir)) { - //quadratic bounce back with for thin walls const int invDir = D3Q27System::INVDIR[fdir]; - fReturn = distributionsTemp->getDistributionInvForDirection(x1+D3Q27System::DX1[invDir], x2+D3Q27System::DX2[invDir], x3+D3Q27System::DX3[invDir], fdir); - distributions->setDistributionForDirection(fReturn, x1+D3Q27System::DX1[invDir], x2+D3Q27System::DX2[invDir], x3+D3Q27System::DX3[invDir], fdir); + if (pass == 1) + { + LBMReal q = bcPtr->getQ(invDir); + LBMReal fReturn = ((1.0 - q) / (1.0 + q))*0.5*(f[invDir] - f[fdir] + (f[invDir] + f[fdir] - collFactor*(feq[fdir] + feq[invDir])) / (1.0 - collFactor)); + distributionsTemp->setDistributionForDirection(fReturn, x1 + D3Q27System::DX1[invDir], x2 + D3Q27System::DX2[invDir], x3 + D3Q27System::DX3[invDir], fdir); + } + else + { + //quadratic bounce back with for thin walls + fReturn = distributionsTemp->getDistributionInvForDirection(x1 + D3Q27System::DX1[invDir], x2 + D3Q27System::DX2[invDir], x3 + D3Q27System::DX3[invDir], fdir); + distributions->setDistributionForDirection(fReturn, x1 + D3Q27System::DX1[invDir], x2 + D3Q27System::DX2[invDir], x3 + D3Q27System::DX3[invDir], fdir); + } } } } ////////////////////////////////////////////////////////////////////////// -void ThinWallNoSlipBoundaryCondition::addDistributions(EsoTwist3DPtr distributions) +void ThinWallNoSlipBoundaryCondition::addDistributions(DistributionArray3DPtr distributions) { this->distributions = distributions; distributionsTemp = EsoTwist3DPtr(new D3Q27EsoTwist3DSplittedVector(distributions->getNX1(), distributions->getNX2(), distributions->getNX3(), -999.0)); } - +////////////////////////////////////////////////////////////////////////// +void ThinWallNoSlipBoundaryCondition::setPass(int pass) +{ + this->pass = pass; +} diff --git a/source/VirtualFluidsCore/BoundaryCondition/ThinWallNoSlipBoundaryCondition.h b/source/VirtualFluidsCore/BoundaryCondition/ThinWallNoSlipBoundaryCondition.h index d4d49b751..82f4f46d4 100644 --- a/source/VirtualFluidsCore/BoundaryCondition/ThinWallNoSlipBoundaryCondition.h +++ b/source/VirtualFluidsCore/BoundaryCondition/ThinWallNoSlipBoundaryCondition.h @@ -6,16 +6,18 @@ class ThinWallNoSlipBoundaryCondition; typedef boost::shared_ptr<ThinWallNoSlipBoundaryCondition> ThinWallNoSlipBoundaryConditionPtr; -class ThinWallNoSlipBoundaryCondition : public NoSlipBoundaryCondition +class ThinWallNoSlipBoundaryCondition : public BoundaryCondition { public: ThinWallNoSlipBoundaryCondition(); virtual ~ThinWallNoSlipBoundaryCondition(); BoundaryConditionPtr clone(); - void addDistributions(EsoTwist3DPtr distributions); + void addDistributions(DistributionArray3DPtr distributions); + void setPass(int pass); protected: void applyBC(); private: + int pass; friend class boost::serialization::access; template<class Archive> void serialize(Archive & ar, const unsigned int version) diff --git a/source/VirtualFluidsCore/BoundaryCondition/VelocityBoundaryCondition.cpp b/source/VirtualFluidsCore/BoundaryCondition/VelocityBoundaryCondition.cpp index fb04d2f3e..e74498925 100644 --- a/source/VirtualFluidsCore/BoundaryCondition/VelocityBoundaryCondition.cpp +++ b/source/VirtualFluidsCore/BoundaryCondition/VelocityBoundaryCondition.cpp @@ -17,6 +17,11 @@ BoundaryConditionPtr VelocityBoundaryCondition::clone() return bc; } ////////////////////////////////////////////////////////////////////////// +void VelocityBoundaryCondition::addDistributions(DistributionArray3DPtr distributions) +{ + this->distributions = distributions; +} +////////////////////////////////////////////////////////////////////////// void VelocityBoundaryCondition::applyBC() { LBMReal f[D3Q27System::ENDF+1]; diff --git a/source/VirtualFluidsCore/BoundaryCondition/VelocityBoundaryCondition.h b/source/VirtualFluidsCore/BoundaryCondition/VelocityBoundaryCondition.h index 845172548..c56c4b942 100644 --- a/source/VirtualFluidsCore/BoundaryCondition/VelocityBoundaryCondition.h +++ b/source/VirtualFluidsCore/BoundaryCondition/VelocityBoundaryCondition.h @@ -12,6 +12,7 @@ public: VelocityBoundaryCondition(); ~VelocityBoundaryCondition(); BoundaryConditionPtr clone(); + void addDistributions(DistributionArray3DPtr distributions); protected: void applyBC(); private: diff --git a/source/VirtualFluidsCore/CoProcessors/MacroscopicQuantitiesCoProcessor.cpp b/source/VirtualFluidsCore/CoProcessors/MacroscopicQuantitiesCoProcessor.cpp index 9dc9cf340..f854a0cc7 100644 --- a/source/VirtualFluidsCore/CoProcessors/MacroscopicQuantitiesCoProcessor.cpp +++ b/source/VirtualFluidsCore/CoProcessors/MacroscopicQuantitiesCoProcessor.cpp @@ -55,7 +55,7 @@ void MacroscopicQuantitiesCoProcessor::process(double step) if(scheduler->isDue(step) ) collectData(step); - UBLOG(logDEBUG3, "D3Q27MacroscopicQuantitiesCoProcessor::update:" << step); + UBLOG(logDEBUG3, "MacroscopicQuantitiesCoProcessor::update:" << step); } ////////////////////////////////////////////////////////////////////////// void MacroscopicQuantitiesCoProcessor::collectData(double step) @@ -121,7 +121,7 @@ void MacroscopicQuantitiesCoProcessor::collectData(double step) { WbWriterVtkXmlASCII::getInstance()->addFilesToCollection(cfilePath,filenames,istep,false); } - UBLOG(logINFO,"D3Q27MacroscopicQuantitiesCoProcessor step: " << istep); + UBLOG(logINFO,"MacroscopicQuantitiesCoProcessor step: " << istep); } clearData(); @@ -152,8 +152,8 @@ void MacroscopicQuantitiesCoProcessor::addDataMQ(Block3DPtr block) datanames.push_back("Vy"); datanames.push_back("Vz"); //datanames.push_back("Press"); - //datanames.push_back("Level"); - //datanames.push_back("BlockID"); + datanames.push_back("Level"); + datanames.push_back("BlockID"); @@ -255,8 +255,8 @@ void MacroscopicQuantitiesCoProcessor::addDataMQ(Block3DPtr block) //data[index++].push_back(vx2 * conv->getFactorVelocityLbToW2()); //data[index++].push_back(vx3 * conv->getFactorVelocityLbToW2()); //data[index++].push_back(press * conv->getFactorPressureLbToW2()); - //data[index++].push_back(level); - //data[index++].push_back(blockID); + data[index++].push_back(level); + data[index++].push_back(blockID); } } } diff --git a/source/VirtualFluidsCore/Connectors/Block3DConnectorFactory.cpp b/source/VirtualFluidsCore/Connectors/Block3DConnectorFactory.cpp index 44fc37843..061a87067 100644 --- a/source/VirtualFluidsCore/Connectors/Block3DConnectorFactory.cpp +++ b/source/VirtualFluidsCore/Connectors/Block3DConnectorFactory.cpp @@ -1,5 +1,5 @@ #include "Block3DConnectorFactory.h" -#include "D3Q27ETFullDirectConnector.h" +#include "D3Q27ETFullDirectConnector2.h" #include "D3Q27ETFullVectorConnector.h" #include "CoarseToFineNodeSetBlock3DConnector.h" #include "FineToCoarseNodeSetBlock3DConnector.h" @@ -14,7 +14,7 @@ Block3DConnectorFactory::~Block3DConnectorFactory() ////////////////////////////////////////////////////////////////////////// Block3DConnectorPtr Block3DConnectorFactory::createSameLevelDirectConnector(Block3DPtr from, Block3DPtr to, int sendDir) { - return Block3DConnectorPtr(new D3Q27ETFullDirectConnector(from, to, sendDir)); + return Block3DConnectorPtr(new D3Q27ETFullDirectConnector2(from, to, sendDir)); } ////////////////////////////////////////////////////////////////////////// Block3DConnectorPtr Block3DConnectorFactory::createSameLevelVectorConnector(Block3DPtr block, diff --git a/source/VirtualFluidsCore/Grid/BoostSerializationClassExportHelper.h b/source/VirtualFluidsCore/Grid/BoostSerializationClassExportHelper.h index 345447a67..854b52b76 100644 --- a/source/VirtualFluidsCore/Grid/BoostSerializationClassExportHelper.h +++ b/source/VirtualFluidsCore/Grid/BoostSerializationClassExportHelper.h @@ -5,7 +5,6 @@ #include <LBMKernelETD3Q27Cascaded.h> #include <D3Q27EsoTwist3DSplittedVector.h> #include <D3Q27ETBCProcessor.h> -#include <D3Q27ETForThinWallBCProcessor.h> #include <LBMKernelETD3Q27CascadedTI.h> #include <DataSet3D.h> #include <LBMKernelETD3Q27BGK.h> @@ -51,7 +50,6 @@ BOOST_CLASS_EXPORT(LBMKernelETD3Q27CCLB) BOOST_CLASS_EXPORT(LBMKernelETD3Q27CCLBWithSpongeLayer) BOOST_CLASS_EXPORT(D3Q27EsoTwist3DSplittedVector) BOOST_CLASS_EXPORT(D3Q27ETBCProcessor) -BOOST_CLASS_EXPORT(D3Q27ETForThinWallBCProcessor) BOOST_CLASS_EXPORT(DataSet3D) BOOST_CLASS_EXPORT(Interactor3D) BOOST_CLASS_EXPORT(D3Q27Interactor) diff --git a/source/VirtualFluidsCore/Grid/CalculationManager.cpp b/source/VirtualFluidsCore/Grid/CalculationManager.cpp index 80a86a371..a159a3726 100644 --- a/source/VirtualFluidsCore/Grid/CalculationManager.cpp +++ b/source/VirtualFluidsCore/Grid/CalculationManager.cpp @@ -43,17 +43,6 @@ CalculationManager::CalculationManager(Grid3DPtr grid, int numOfThreads, double loadBalancer = LoadBalancerPtr(new LoadBalancer(grid, comm, endDir)); } ////////////////////////////////////////////////////////////////////////// -CalculationManager::CalculationManager(Grid3DPtr grid, int numOfThreads, double endTime, UbSchedulerPtr visScheduler, BoundaryConditionProcessorPtr bcProcessor, CalculationManager::CalculatorType calcType) - : grid(grid), - numOfThreads(numOfThreads), - endTime(endTime), - visScheduler(visScheduler), - bcProcessor(bcProcessor), - calcType(calcType) -{ - init(); -} -////////////////////////////////////////////////////////////////////////// CalculationManager::~CalculationManager() { } @@ -195,13 +184,13 @@ CalculatorPtr CalculationManager::createCalculator(Grid3DPtr grid, SynchronizerP switch (calcType) { case CalculationManager::MPI: - return CalculatorPtr (new Calculator(grid, sync, bcProcessor, mainThread)); + return CalculatorPtr (new Calculator(grid, sync, mainThread)); #if defined VF_FETOL case CalculationManager::FETOL: return CalculatorPtr (new FETOLCalculator(grid, sync, mainThread)); #endif case CalculationManager::PrePostBc: - return CalculatorPtr(new PrePostBcCalculator(grid, sync, bcProcessor, mainThread)); + return CalculatorPtr(new PrePostBcCalculator(grid, sync, mainThread)); default: UB_THROW(UbException(UB_EXARGS,"This Calculator is not defined!")); } diff --git a/source/VirtualFluidsCore/Grid/CalculationManager.h b/source/VirtualFluidsCore/Grid/CalculationManager.h index 2d8c06bf2..2cd8155b8 100644 --- a/source/VirtualFluidsCore/Grid/CalculationManager.h +++ b/source/VirtualFluidsCore/Grid/CalculationManager.h @@ -17,7 +17,6 @@ public: CalculationManager(Grid3DPtr grid, int numOfThreads, double endTime, UbSchedulerPtr visScheduler, CalculatorType calcType = CalculationManager::MPI); CalculationManager(Grid3DPtr grid, int numOfThreads, double endTime, UbSchedulerPtr visScheduler, CommunicatorPtr comm, int endDir, LBMReal nu, CalculatorType calcType = CalculationManager::MPI); - CalculationManager(Grid3DPtr grid, int numOfThreads, double endTime, UbSchedulerPtr visScheduler, BoundaryConditionProcessorPtr bcProcessor, CalculatorType calcType = CalculationManager::MPI); virtual ~CalculationManager(); void calculate(); void setVisScheduler(UbSchedulerPtr s); @@ -42,7 +41,6 @@ private: LoadBalancerPtr loadBalancer; int rank; LBMReal nu; - BoundaryConditionProcessorPtr bcProcessor; }; #endif diff --git a/source/VirtualFluidsCore/Grid/Calculator.cpp b/source/VirtualFluidsCore/Grid/Calculator.cpp index e7492419f..bf0cef29d 100644 --- a/source/VirtualFluidsCore/Grid/Calculator.cpp +++ b/source/VirtualFluidsCore/Grid/Calculator.cpp @@ -34,29 +34,6 @@ Calculator::Calculator(Grid3DPtr grid, SynchronizerPtr sync, bool mainThread) : loadBalancingComp = false; } ////////////////////////////////////////////////////////////////////////// -Calculator::Calculator(Grid3DPtr grid, SynchronizerPtr sync, BoundaryConditionProcessorPtr bcProcessor, bool mainThread) : -grid(grid), -sync(sync), -bcProcessor(bcProcessor), -mainThread(mainThread), -refinement(false) -{ - minLevel = grid->getCoarsestInitializedLevel(); - maxLevel = grid->getFinestInitializedLevel(); - if (maxLevel > 0) - refinement = true; - else - refinement = false; - blocks.resize(maxLevel+1); - localConns.resize(maxLevel+1); - remoteConns.resize(maxLevel+1); - localInterfaceBlockConns.resize(maxLevel+1); - remoteInterfaceBlockConns.resize(maxLevel+1); - localInterConns.resize(maxLevel); - remoteInterConns.resize(maxLevel); - loadBalancingComp = false; -} -////////////////////////////////////////////////////////////////////////// void Calculator::calculate(const double& endTime, CalculationManagerPtr cm, boost::exception_ptr& error) { UBLOG(logDEBUG1, "Calculator::calculate() - started"); @@ -593,18 +570,18 @@ void Calculator::setTimeAveragedValuesCoProcessor(TimeAveragedValuesCoProcessorP } ////////////////////////////////////////////////////////////////////////// -void Calculator::applyBCs( int startLevel, int maxInitLevel ) -{ - //startLevel bis maxInitLevel - for(int level=startLevel; level<=maxInitLevel; level++) - { - //call LBM kernel - BOOST_FOREACH(Block3DPtr block, blocks[level]) - { - block->getKernel()->getBCProcessor()->applyBC(); - } - } -} +//void Calculator::applyBCs( int startLevel, int maxInitLevel ) +//{ +// //startLevel bis maxInitLevel +// for(int level=startLevel; level<=maxInitLevel; level++) +// { +// //call LBM kernel +// BOOST_FOREACH(Block3DPtr block, blocks[level]) +// { +// block->getKernel()->getBCProcessor()->applyBC(); +// } +// } +//} ////////////////////////////////////////////////////////////////////////// void Calculator::applyPreCollisionBC(int startLevel, int maxInitLevel) { diff --git a/source/VirtualFluidsCore/Grid/Calculator.h b/source/VirtualFluidsCore/Grid/Calculator.h index 1db546417..494cd72e3 100644 --- a/source/VirtualFluidsCore/Grid/Calculator.h +++ b/source/VirtualFluidsCore/Grid/Calculator.h @@ -8,7 +8,6 @@ #include "basics/utilities/UbScheduler.h" #include "basics/utilities/UbTiming.h" #include "LoadBalancer.h" -#include "BoundaryConditionProcessor.h" #include "TimeAveragedValuesCoProcessor.h" @@ -22,7 +21,6 @@ class Calculator public: Calculator(); Calculator(Grid3DPtr grid, SynchronizerPtr sync, bool mainThread = true); - Calculator(Grid3DPtr grid, SynchronizerPtr sync, BoundaryConditionProcessorPtr bcProcessor, bool mainThread = true); virtual ~Calculator(){} virtual void calculate(const double& endTime, CalculationManagerPtr cm, boost::exception_ptr& error); void addBlock(Block3DPtr block); @@ -49,7 +47,7 @@ protected: void connectorsSetInvStep(std::vector< Block3DConnectorPtr >& connectors, bool invStep); void interpolation(int startLevel, int maxInitLevel); void deleteConnectors(std::vector< std::vector< Block3DConnectorPtr > >& conns); - void applyBCs(int startLevel, int maxInitLevel); + //void applyBCs(int startLevel, int maxInitLevel); void applyPreCollisionBC(int startLevel, int maxInitLevel); void applyPostCollisionBC(int startLevel, int maxInitLevel); int minLevel, maxLevel; @@ -65,7 +63,6 @@ protected: Grid3DPtr grid; UbSchedulerPtr visScheduler; int calcStep; - BoundaryConditionProcessorPtr bcProcessor; std::vector< std::vector<Block3DPtr> > blocks; private: diff --git a/source/VirtualFluidsCore/Grid/Calculator2.cpp b/source/VirtualFluidsCore/Grid/Calculator2.cpp deleted file mode 100644 index 769765e67..000000000 --- a/source/VirtualFluidsCore/Grid/Calculator2.cpp +++ /dev/null @@ -1,546 +0,0 @@ -#include "Calculator2.h" -#include <basics/utilities/UbException.h> -#include <boost/foreach.hpp> -#include "SimulationParameters.h" -#include "MathUtil.hpp" -#include "basics/writer/WbWriterVtkXmlASCII.h" - -//#define TIMING - - -Calculator2::Calculator2(Grid3DPtr grid, SynchronizerPtr sync, bool mainThread) : - grid(grid), - sync(sync), - mainThread(mainThread), - refinement(false) -{ - minLevel = grid->getCoarsestInitializedLevel(); - maxLevel = grid->getFinestInitializedLevel(); - if(maxLevel > 0) - refinement = true; - else - refinement = false; - blocks.resize(maxLevel+1); - localConns.resize(maxLevel+1); - remoteConns.resize(maxLevel+1); - localInterfaceBlockConns.resize(maxLevel+1); - remoteInterfaceBlockConns.resize(maxLevel+1); - localInterConns.resize(maxLevel); - remoteInterConns.resize(maxLevel); - loadBalancingComp = false; -} -////////////////////////////////////////////////////////////////////////// -void Calculator2::calculate(const double& endTime, CalculationManagerPtr cm, boost::exception_ptr& error) -{ - UBLOG(logDEBUG1, "Calculator2::calculate() - started"); - try - { - 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; calcStep++) - { - - //exchange data between blocks for visualization - //sync->wait(); - ////if(visScheduler->isDue((double)(calcStep-1))) - ////{ - // //exchangeBlockData(minInitLevel, maxInitLevel, true); - ////} - - ////wait for write dump files - //sync->wait(); - //write dump - if (mainThread) grid->coProcess((double)(calcStep-1)); - sync->wait(); - -////////////////////////////////////////////////////////////////////////// -#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 -////////////////////////////////////////////////////////////////////////// - calculateBlocks(straightStartLevel, maxInitLevel); - //calculateBlocks(minInitLevel, maxInitLevel, staggeredStep); -////////////////////////////////////////////////////////////////////////// -#ifdef TIMING - time[0] = timer.stop(); - //UBLOG(logINFO, "calculateBlocks time = " <<time); -#endif -////////////////////////////////////////////////////////////////////////// - - //exchange data between blocks - exchangeBlockData(straightStartLevel, maxInitLevel, false); -////////////////////////////////////////////////////////////////////////// -#ifdef TIMING - time[1] = timer.stop(); - //UBLOG(logINFO, "exchangeBlockData time = " <<time); -#endif -////////////////////////////////////////////////////////////////////////// - - applyBCs(straightStartLevel, maxInitLevel); -////////////////////////////////////////////////////////////////////////// -#ifdef TIMING - time[2] = timer.stop(); - //UBLOG(logINFO, "applyBCs time = " <<time); -#endif -////////////////////////////////////////////////////////////////////////// - - //swap distributions in kernel - swapDistributions(straightStartLevel, maxInitLevel); -////////////////////////////////////////////////////////////////////////// -#ifdef TIMING - time[3] = timer.stop(); - //UBLOG(logINFO, "swapDistributions time = " <<time); -#endif -////////////////////////////////////////////////////////////////////////// - - if (refinement) - { - //exchange data between blocks for grid refinement - //exchangeInterfaceBlockData(straightStartLevel, maxInitLevel, true); - if(straightStartLevel<maxInitLevel) - //exchangeBlockData(straightStartLevel, maxInitLevel, true); - exchangeInterfaceBlockData(straightStartLevel, maxInitLevel, true); -////////////////////////////////////////////////////////////////////////// -#ifdef TIMING - time[4] = timer.stop(); - UBLOG(logINFO, "refinement exchangeBlockData time = " <<time); -#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); -#endif -////////////////////////////////////////////////////////////////////////// - } - - } - //exchange data between blocks for visualization - if((int)visScheduler->getNextDueTime() == calcStep) - { - if(mainThread) visScheduler->isDue((double)(calcStep-1)); - exchangeBlockData(straightStartLevel, maxInitLevel, true); - } - //now ghost nodes have actual values - - //dynamic load balancing - //sync->wait(); - //if (mainThread && !loadBalancingComp) - //{ - // loadBalancingComp = cm->balance(); - //} - } - error = boost::exception_ptr(); - UBLOG(logDEBUG1, "Calculator2::calculate() - stoped"); - } - catch( ... ) - { - error = boost::current_exception(); - } -} -////////////////////////////////////////////////////////////////////////// -void Calculator2::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 Calculator2::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 Calculator2::addBlock(Block3DPtr block) -{ - blocks[block->getLevel()].push_back(block); -} -////////////////////////////////////////////////////////////////////////// -void Calculator2::initConnectors() -{ - UBLOG(logDEBUG1, "Calculator2::initLocalConnectors() - start"); - - for (int l = minLevel; l <= maxLevel; l++) - { - BOOST_FOREACH(Block3DPtr block, blocks[l]) - { - block->pushBackLocalSameLevelConnectors(localConns[l]); - - if(block->hasInterpolationFlag()) - block->pushBackLocalSameLevelConnectors(localInterfaceBlockConns[l]); - if (l != maxLevel) - block->pushBackLocalInterpolationConnectorsCF(localInterConns[l]); - } - if (l != maxLevel) - { - BOOST_FOREACH(Block3DPtr block, blocks[l+1]) - { - block->pushBackLocalInterpolationConnectorsFC(localInterConns[l]); - } - } - UBLOG(logDEBUG5, "Calculator2::initConnectors()-initConnectors(localConns["<<l<<"])"); - initConnectors(localConns[l]); - - if (l != maxLevel) - { - UBLOG(logDEBUG5, "Calculator2::initConnectors()-initConnectors(localInterConns["<<l<<"])"); - initConnectors(localInterConns[l]); - } - } - - if (mainThread) - initRemoteConnectors(); - - UBLOG(logDEBUG1, "Calculator2::initLocalConnectors() - end"); -} -////////////////////////////////////////////////////////////////////////// -void Calculator2::initRemoteConnectors() -{ - std::vector< std::vector< Block3DConnectorPtr > > remoteInterConnsCF; - std::vector< std::vector< Block3DConnectorPtr > > remoteInterConnsFC; - remoteInterConnsCF.resize(maxLevel+1); - remoteInterConnsFC.resize(maxLevel+1); - - int minInitLevel = this->grid->getCoarsestInitializedLevel(); - int maxInitLevel = this->grid->getFinestInitializedLevel(); - int gridRank = grid->getRank(); - - for(int level = minInitLevel; level<=maxInitLevel;level++) - { - std::vector<Block3DPtr> blockVector; - //grid->getBlocks(level, gridRank, true, blockVector); - grid->getBlocks(level, blockVector); - BOOST_FOREACH(Block3DPtr block, blockVector) - { - int l = block->getLevel(); - block->pushBackRemoteSameLevelConnectors(remoteConns[l]); - - //if(block->isInterface()) - // block->pushBackRemoteSameLevelConnectors(remoteInterfaceBlockConns[l]); - block->pushBackRemoteInterpolationConnectorsCF(remoteInterConnsCF[l]); - block->pushBackRemoteInterpolationConnectorsFC(remoteInterConnsFC[l]); - } - } - - for (int l = minLevel; l <= maxLevel; l++) - { - UBLOG(logDEBUG5, "Calculator2::initRemoteConnectors()-initConnectors(remoteConns["<<l<<"])"); - initConnectors(remoteConns[l]); - if (l != maxLevel) - { - for(int i = 0; i < remoteInterConnsCF[l].size(); i++) - remoteInterConns[l].push_back(remoteInterConnsCF[l][i]); - for(int i = 0; i < remoteInterConnsFC[l+1].size(); i++) - remoteInterConns[l].push_back(remoteInterConnsFC[l+1][i]); - //UBLOG(logDEBUG5, "Calculator2::initRemoteConnectors()-initConnectors(remoteInterConns["<<l<<"])"); - //initConnectors(remoteInterConns[l]); - } - } - ////////////////////////////////////////////////////////////////////////// - //UBLOG(logDEBUG5, "Calculator2::initConnectors() - connectoren initialisieren - start"); - for (int l = minLevel; l <= maxLevel; l++) - { - if (l != maxLevel) - { - UBLOG(logDEBUG5, "Calculator2::initRemoteConnectors()-initConnectors(remoteInterConns["<<l<<"])"); - BOOST_FOREACH(Block3DConnectorPtr c, remoteInterConns[l] ) c->init(); - } - } - //UBLOG(logDEBUG5, "Calculator2::initConnectors() - connectoren initialisieren - end"); - ////////////////////////////////////////////////////////////////////////// - //sendTransmitterDataSize - //UBLOG(logDEBUG5, "Calculator2::initConnectors() - sendTransmitterDataSize - start"); - for (int l = minLevel; l <= maxLevel; l++) - { - if (l != maxLevel) - { - UBLOG(logDEBUG5, "Calculator2::initRemoteConnectors()-sendTransmitterDataSize(remoteInterConns["<<l<<"])"); - BOOST_FOREACH(Block3DConnectorPtr c, remoteInterConns[l] ) c->sendTransmitterDataSize(); - } - } - //UBLOG(logDEBUG5, "Calculator2::initConnectors() - sendTransmitterDataSize - end"); - ////////////////////////////////////////////////////////////////////////// - //receiveTransmitterDataSize - //wenn er hier bei verteilten berechnungen stopped, dann ist vermutlich auf einer seite ein nicht aktiver block!!! - //UBLOG(logDEBUG5, "Calculator2::initConnectors() - receiveTransmitterDataSize - start"); - for (int l = minLevel; l <= maxLevel; l++) - { - if (l != maxLevel) - { - UBLOG(logDEBUG5, "Calculator2::initRemoteConnectors()-receiveTransmitterDataSize(remoteInterConns["<<l<<"])"); - BOOST_FOREACH(Block3DConnectorPtr c, remoteInterConns[l] ) c->receiveTransmitterDataSize(); - } - } - //UBLOG(logDEBUG5, "Calculator2::initConnectors() - receiveTransmitterDataSize - end"); - ////////////////////////////////////////////////////////////////////////// -} -////////////////////////////////////////////////////////////////////////// -void Calculator2::initConnectors(std::vector<Block3DConnectorPtr>& connectors) -{ - UBLOG(logDEBUG1, "Calculator2::initConnectors() - start"); - - //initialization - ////////////////////////////////////////////////////////////////////////// - //initialize connectors - UBLOG(logDEBUG5, "Calculator2::initConnectors() - connectoren initialisieren - start"); - BOOST_FOREACH(Block3DConnectorPtr c, connectors ) c->init(); - UBLOG(logDEBUG5, "Calculator2::initConnectors() - connectoren initialisieren - end"); - ////////////////////////////////////////////////////////////////////////// - //sendTransmitterDataSize - UBLOG(logDEBUG5, "Calculator2::initConnectors() - sendTransmitterDataSize - start"); - BOOST_FOREACH(Block3DConnectorPtr c, connectors ) c->sendTransmitterDataSize(); - UBLOG(logDEBUG5, "Calculator2::initConnectors() - sendTransmitterDataSize - end"); - ////////////////////////////////////////////////////////////////////////// - //receiveTransmitterDataSize - //wenn er hier bei verteilten berechnungen stopped, dann ist vermutlich auf einer seite ein nicht aktiver block!!! - UBLOG(logDEBUG5, "Calculator2::initConnectors() - receiveTransmitterDataSize - start"); - BOOST_FOREACH(Block3DConnectorPtr c, connectors ) c->receiveTransmitterDataSize(); - UBLOG(logDEBUG5, "Calculator2::initConnectors() - receiveTransmitterDataSize - end"); - - UBLOG(logDEBUG1, "Calculator2::initConnectors() - end"); -} -////////////////////////////////////////////////////////////////////////// -void Calculator2::exchangeBlockData(int startLevel, int maxInitLevel, bool invStep) -{ - sync->wait(); - //startLevel bis maxInitLevel - for(int level=startLevel; level<=maxInitLevel; level++) - { - connectorsPrepare(localConns[level]); - connectorsPrepare(remoteConns[level]); - - connectorsSend(localConns[level], invStep); - connectorsSend(remoteConns[level], invStep); - - connectorsReceive(localConns[level], invStep); - connectorsReceive(remoteConns[level], invStep); - } - sync->wait(); -} -////////////////////////////////////////////////////////////////////////// -void Calculator2::exchangeInterfaceBlockData(int startLevel, int maxInitLevel, bool invStep) -{ - sync->wait(); - //startLevel bis maxInitLevel - for(int level=startLevel; level<=maxInitLevel; level++) - { - connectorsPrepare(localInterfaceBlockConns[level]); - connectorsPrepare(remoteInterfaceBlockConns[level]); - - connectorsSend(localInterfaceBlockConns[level], invStep); - connectorsSend(remoteInterfaceBlockConns[level], invStep); - - connectorsReceive(localInterfaceBlockConns[level], invStep); - connectorsReceive(remoteInterfaceBlockConns[level], invStep); - } - sync->wait(); -} -////////////////////////////////////////////////////////////////////////// -void Calculator2::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 Calculator2::connectorsPrepare(std::vector< Block3DConnectorPtr >& connectors) -{ - BOOST_FOREACH(Block3DConnectorPtr c, connectors) - { - c->prepareForReceive(); - c->prepareForSend(); - } -} -////////////////////////////////////////////////////////////////////////// -void Calculator2::connectorsSend(std::vector< Block3DConnectorPtr >& connectors, bool invStep) -{ - BOOST_FOREACH(Block3DConnectorPtr c, connectors) - { - c->setInvStep(invStep); - c->fillSendVectors(); - c->sendVectors(); - } -} -////////////////////////////////////////////////////////////////////////// -void Calculator2::connectorsReceive(std::vector< Block3DConnectorPtr >& connectors, bool invStep) -{ - BOOST_FOREACH(Block3DConnectorPtr c, connectors) - { - c->setInvStep(invStep); - c->receiveVectors(); - c->distributeReceiveVectors(); - } -} -////////////////////////////////////////////////////////////////////////// -void Calculator2::connectorsSetInvStep(std::vector< Block3DConnectorPtr >& connectors, bool invStep) -{ - BOOST_FOREACH(Block3DConnectorPtr c, connectors) - { - c->setInvStep(invStep); - } -} -////////////////////////////////////////////////////////////////////////// -void Calculator2::interpolation(int startLevel, int maxInitLevel) -{ - sync->wait(); - - for(int level=startLevel; level<maxInitLevel; level++) - { - connectorsPrepare(localInterConns[level]); - connectorsPrepare(remoteInterConns[level]); - } - - sync->wait(); - - for(int level=startLevel; level<maxInitLevel; level++) - { - connectorsSend(localInterConns[level], true); - connectorsSend(remoteInterConns[level], true); - } - - sync->wait(); - - for(int level=startLevel; level<maxInitLevel; level++) - { - connectorsReceive(localInterConns[level], true); - connectorsReceive(remoteInterConns[level], true); - } - - sync->wait(); - -} -////////////////////////////////////////////////////////////////////////// -void Calculator2::setVisScheduler(UbSchedulerPtr s) -{ - visScheduler = s; -} -////////////////////////////////////////////////////////////////////////// -//double Calculator2::getCallculationTime() -//{ -// return timer.getTotalTime(); -//} -////////////////////////////////////////////////////////////////////////// -std::vector< std::vector< Block3DPtr > > Calculator2::getBlocks() -{ - return blocks; -} -////////////////////////////////////////////////////////////////////////// -void Calculator2::deleteBlocks() -{ - BOOST_FOREACH(std::vector< Block3DPtr > &bs, blocks) - bs.resize(0); -} -////////////////////////////////////////////////////////////////////////// -void Calculator2::deleteConnectors() -{ - deleteConnectors(localConns); - deleteConnectors(remoteConns); - - deleteConnectors(localInterfaceBlockConns); - deleteConnectors(remoteInterfaceBlockConns); - - deleteConnectors(localInterConns); - deleteConnectors(remoteInterConns); -} -////////////////////////////////////////////////////////////////////////// -void Calculator2::deleteConnectors(std::vector< std::vector< Block3DConnectorPtr > >& conns) -{ - BOOST_FOREACH(std::vector< Block3DConnectorPtr > &c, conns) - c.resize(0); -} -////////////////////////////////////////////////////////////////////////// -void Calculator2::applyBCs( int startLevel, int maxInitLevel ) -{ - //startLevel bis maxInitLevel - for(int level=startLevel; level<=maxInitLevel; level++) - { - //call LBM kernel - BOOST_FOREACH(Block3DPtr block, blocks[level]) - { - block->getKernel()->getBCProcessor()->applyBC(); - } - } -} diff --git a/source/VirtualFluidsCore/Grid/Calculator2.h b/source/VirtualFluidsCore/Grid/Calculator2.h deleted file mode 100644 index 7c229963c..000000000 --- a/source/VirtualFluidsCore/Grid/Calculator2.h +++ /dev/null @@ -1,78 +0,0 @@ -#ifndef Calculator2_H -#define Calculator2_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" - - -class Calculator2; -typedef boost::shared_ptr<Calculator2> Calculator2Ptr; - -#include "CalculationManager.h" - -class Calculator2 -{ -public: - Calculator2(Grid3DPtr grid, SynchronizerPtr sync, bool mainThread = true); - virtual ~Calculator2(){} - void calculate(const double& endTime, CalculationManagerPtr cm, boost::exception_ptr& error); - void addBlock(Block3DPtr block); - void initConnectors(); - void setVisScheduler(UbSchedulerPtr s); - //double getCallculationTime(); - std::vector< std::vector< Block3DPtr > > getBlocks(); - void deleteBlocks(); - void deleteConnectors(); -protected: - void calculateBlocks(int startLevel, int maxInitLevel); - void calculateBlocks(int minInitLevel, int maxInitLevel, int staggeredStep); - void initConnectors(std::vector<Block3DConnectorPtr>& connectors); - void initRemoteConnectors(); - void swapDistributions(int startLevel, int maxInitLevel); - void exchangeBlockData(int startLevel, int maxInitLevel, bool invStep); - void exchangeInterfaceBlockData(int startLevel, int maxInitLevel, bool invStep); - void connectorsPrepare(std::vector< Block3DConnectorPtr >& connectors); - void connectorsSend(std::vector< Block3DConnectorPtr >& connectors, bool invStep); - void connectorsReceive(std::vector< Block3DConnectorPtr >& connectors, bool invStep); - void connectorsSetInvStep(std::vector< Block3DConnectorPtr >& connectors, bool invStep); - void interpolation(int startLevel, int maxInitLevel); - void deleteConnectors(std::vector< std::vector< Block3DConnectorPtr > >& conns); - void applyBCs(int startLevel, int maxInitLevel); -private: - boost::barrier* bar; - //double time; - SynchronizerPtr sync; - int minLevel, maxLevel; - bool mainThread; - bool refinement; - Grid3DPtr grid; - UbSchedulerPtr visScheduler; - int calcStep; - - std::vector< std::vector<Block3DPtr> > blocks; - std::vector< std::vector< Block3DConnectorPtr > > localConns; - std::vector< std::vector< Block3DConnectorPtr > > remoteConns; - - std::vector< std::vector< Block3DConnectorPtr > > localInterfaceBlockConns; - std::vector< std::vector< Block3DConnectorPtr > > remoteInterfaceBlockConns; - - //localInterConns and remoteInterConns save interpolation connectors - //every element save CF connectors for current level and FC connectors for next level - //e.g. - //localInterConns[0] = CF(0), FC(1) - //localInterConns[1] = CF(1), FC(2) - //localInterConns[2] - std::vector< std::vector< Block3DConnectorPtr > > localInterConns; - std::vector< std::vector< Block3DConnectorPtr > > remoteInterConns; - - //UbTimer timer, timer2, timer3; - bool loadBalancingComp; - -}; - -#endif diff --git a/source/VirtualFluidsCore/Grid/FastSignal.hpp b/source/VirtualFluidsCore/Grid/FastSignal.hpp deleted file mode 100644 index f8d9b0209..000000000 --- a/source/VirtualFluidsCore/Grid/FastSignal.hpp +++ /dev/null @@ -1,46 +0,0 @@ -#include <list> -#include <boost/function.hpp> - -namespace fast_signal -{ - template <class Signature> - struct signal: public std::list<boost::function<Signature> >{}; -} - -// Abkürzung für die Schleife zum Aufrufen -#define CALL_SIGNALS( SIG, REF, ARGS )\ - for( fast_signal::signal<SIG>::iterator it = REF.begin();\ - it != REF.end(); ++it )\ - (*it)ARGS; - - -//Eigentlich war das schon alles was wir brauchen um unsere sehr einfache Signal und Slot Implementation zu verwenden. -//Die Syntax ist zwar nicht ganz gleich wie bei Boost.Signals aber das sollten wir für die Zeitersparnis verkraften können. -//Jetzt also ein einfaches Verwendungsbeispiel: -//#include <iostream> -//// ... Code von vorher -// -//// Testfunktion zum Registrieren als Slots -//void test( int& i ) -//{ -// i += 1; -//} -// -//int main() -//{ -// tp_fast_signal::signal<void(int&)> test_signal; // Signal erstellen -// -// // 2 Slots registrieren: -// test_signal.push_back(&test); -// test_signal.push_back(&test); -// -// int test_int = 0; -// -// // Alle Slots aufrufen -// CALL_SIGNALS( void(int&), test_signal, (test_int) ); -// -// // Und noch testweise den aktuellen Wert von test_int ausgeben -// std::cout << "test_int=" << test_int << std::endl; -// -// return 0; -//} diff --git a/source/VirtualFluidsCore/Grid/MTCalculator.cpp b/source/VirtualFluidsCore/Grid/MTCalculator.cpp deleted file mode 100644 index 6a401effe..000000000 --- a/source/VirtualFluidsCore/Grid/MTCalculator.cpp +++ /dev/null @@ -1,589 +0,0 @@ -#include "MTCalculator.h" -#include <basics/utilities/UbException.h> -#include <boost/foreach.hpp> -#include "SimulationParameters.h" -#include "MathUtil.hpp" -#include "basics/writer/WbWriterVtkXmlASCII.h" - -#define TIMING - -MTCalculator::MTCalculator() -{ - -} -////////////////////////////////////////////////////////////////////////// -MTCalculator::MTCalculator(Grid3DPtr grid, SynchronizerPtr sync, bool mainThread) : - grid(grid), - sync(sync), - mainThread(mainThread), - refinement(false) -{ - minLevel = grid->getCoarsestInitializedLevel(); - maxLevel = grid->getFinestInitializedLevel(); - if(maxLevel > 0) - refinement = true; - else - refinement = false; - blocks.resize(maxLevel+1); - localConns.resize(maxLevel+1); - remoteConns.resize(maxLevel+1); - localInterfaceBlockConns.resize(maxLevel+1); - remoteInterfaceBlockConns.resize(maxLevel+1); - localInterConns.resize(maxLevel); - remoteInterConns.resize(maxLevel); - loadBalancingComp = false; -} -////////////////////////////////////////////////////////////////////////// -void MTCalculator::calculate(const double& endTime, CalculationManagerPtr cm, boost::exception_ptr& error) -{ - UBLOG(logDEBUG1, "Calculator::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; calcStep++) - { - - //exchange data between blocks for visualization - //sync->wait(); - ////if(visScheduler->isDue((double)(calcStep-1))) - ////{ - // //exchangeBlockData(minInitLevel, maxInitLevel, true); - ////} - - ////wait for write dump files - //sync->wait(); - //write dump - if (mainThread) grid->coProcess((double)(calcStep-1)); - sync->wait(); - - -////////////////////////////////////////////////////////////////////////// -#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 -////////////////////////////////////////////////////////////////////////// - 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, false); -////////////////////////////////////////////////////////////////////////// -#ifdef TIMING - time[1] = timer.stop(); - UBLOG(logINFO, "exchangeBlockData time = " <<time[1]); -#endif -////////////////////////////////////////////////////////////////////////// - - applyBCs(straightStartLevel, maxInitLevel); -////////////////////////////////////////////////////////////////////////// -#ifdef TIMING - time[2] = timer.stop(); - UBLOG(logINFO, "applyBCs time = " <<time[2]); -#endif -////////////////////////////////////////////////////////////////////////// - - //swap distributions in kernel - swapDistributions(straightStartLevel, maxInitLevel); -////////////////////////////////////////////////////////////////////////// -#ifdef TIMING - time[3] = timer.stop(); - UBLOG(logINFO, "swapDistributions time = " <<time[3]); -#endif -////////////////////////////////////////////////////////////////////////// - - if (refinement) - { - //exchange data between blocks for grid refinement - //exchangeInterfaceBlockData(straightStartLevel, maxInitLevel, true); - if(straightStartLevel<maxInitLevel) - exchangeBlockData(straightStartLevel, maxInitLevel, true); - //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 -////////////////////////////////////////////////////////////////////////// - } - - } - //exchange data between blocks for visualization - if(mainThread) visScheduler->isDue((double)(calcStep-1)); - if((int)visScheduler->getNextDueTime() == calcStep) - { - exchangeBlockData(straightStartLevel, maxInitLevel, true); - } - //now ghost nodes have actual values - - //dynamic load balancing - //sync->wait(); - //if (mainThread && !loadBalancingComp) - //{ - // loadBalancingComp = cm->balance(); - //} - } - error = boost::exception_ptr(); - UBLOG(logDEBUG1, "Calculator::calculate() - stoped"); - } - catch( std::exception& e ) - { - //error = boost::current_exception(); - UBLOG(logERROR, e.what()); - UBLOG(logERROR, " step = "<<calcStep); - exit(EXIT_FAILURE); - } -} -////////////////////////////////////////////////////////////////////////// -void MTCalculator::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 MTCalculator::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 MTCalculator::addBlock(Block3DPtr block) -{ - blocks[block->getLevel()].push_back(block); -} -////////////////////////////////////////////////////////////////////////// -void MTCalculator::initConnectors() -{ - UBLOG(logDEBUG1, "Calculator::initLocalConnectors() - start"); - - for (int l = minLevel; l <= maxLevel; l++) - { - BOOST_FOREACH(Block3DPtr block, blocks[l]) - { - block->pushBackLocalSameLevelConnectors(localConns[l]); - - if(block->hasInterpolationFlag()) - block->pushBackLocalSameLevelConnectors(localInterfaceBlockConns[l]); - if (l != maxLevel) - block->pushBackLocalInterpolationConnectorsCF(localInterConns[l]); - } - if (l != maxLevel) - { - BOOST_FOREACH(Block3DPtr block, blocks[l+1]) - { - block->pushBackLocalInterpolationConnectorsFC(localInterConns[l]); - } - } - UBLOG(logDEBUG5, "Calculator::initConnectors()-initConnectors(localConns["<<l<<"])"); - initConnectors(localConns[l]); - - if (l != maxLevel) - { - UBLOG(logDEBUG5, "Calculator::initConnectors()-initConnectors(localInterConns["<<l<<"])"); - initConnectors(localInterConns[l]); - } - } - - //if (mainThread) - initRemoteConnectors(); - - sync->wait(); - - UBLOG(logDEBUG1, "Calculator::initLocalConnectors() - end"); -} -////////////////////////////////////////////////////////////////////////// -void MTCalculator::initRemoteConnectors() -{ - std::vector< std::vector< Block3DConnectorPtr > > remoteInterConnsCF; - std::vector< std::vector< Block3DConnectorPtr > > remoteInterConnsFC; - remoteInterConnsCF.resize(maxLevel+1); - remoteInterConnsFC.resize(maxLevel+1); - - int minInitLevel = this->grid->getCoarsestInitializedLevel(); - int maxInitLevel = this->grid->getFinestInitializedLevel(); - int gridRank = grid->getRank(); - - - - for(int level = minInitLevel; level<=maxInitLevel;level++) - { - std::vector<Block3DPtr> blockVector; - //grid->getBlocks(level, blockVector); - //BOOST_FOREACH(Block3DPtr block, blockVector) - BOOST_FOREACH(Block3DPtr block, blocks[level]) - { - int l = block->getLevel(); - block->pushBackRemoteSameLevelConnectors(remoteConns[l]); - - //if(block->isInterface()) - // block->pushBackRemoteSameLevelConnectors(remoteInterfaceBlockConns[l]); - block->pushBackRemoteInterpolationConnectorsCF(remoteInterConnsCF[l]); - block->pushBackRemoteInterpolationConnectorsFC(remoteInterConnsFC[l]); - } - } - sync->gmtx.lock(); - for (int l = minLevel; l <= maxLevel; l++) - { - UBLOG(logDEBUG5, "Calculator::initRemoteConnectors()-initConnectors(remoteConns["<<l<<"]).size ="<<remoteConns[l].size()); - initConnectors(remoteConns[l]); - if (l != maxLevel) - { - for(int i = 0; i < remoteInterConnsCF[l].size(); i++) - remoteInterConns[l].push_back(remoteInterConnsCF[l][i]); - for(int i = 0; i < remoteInterConnsFC[l+1].size(); i++) - remoteInterConns[l].push_back(remoteInterConnsFC[l+1][i]); - } - } - sync->gmtx.unlock(); - sync->wait(); - ////////////////////////////////////////////////////////////////////////// - UBLOG(logDEBUG5, "Calculator::initConnectors() - initialize remote connectors - start"); - sync->gmtx.lock(); - for (int l = minLevel; l <= maxLevel; l++) - { - if (l != maxLevel) - { - UBLOG(logDEBUG5, "Calculator::initRemoteConnectors()-initConnectors(remoteInterConns["<<l<<"]).size ="<<remoteInterConns[l].size()); - //boost::lock_guard<boost::mutex> guard(mtx_); - BOOST_FOREACH(Block3DConnectorPtr c, remoteInterConns[l] ) c->init(); - } - } - sync->gmtx.unlock(); - - UBLOG(logDEBUG5, "Calculator::initConnectors() - initialize remote connectors - end"); - - // sync->wait(); - ////////////////////////////////////////////////////////////////////////// - //sendTransmitterDataSize - //UBLOG(logDEBUG5, "Calculator::initConnectors() - sendTransmitterDataSize - start"); - //sync->gmtx.lock(); - for (int l = minLevel; l <= maxLevel; l++) - { - if (l != maxLevel) - { - UBLOG(logDEBUG5, "Calculator::initRemoteConnectors()-sendTransmitterDataSize(remoteInterConns["<<l<<"]).size ="<<remoteInterConns[l].size()); - BOOST_FOREACH(Block3DConnectorPtr c, remoteInterConns[l] ) c->sendTransmitterDataSize(); - } - } - //sync->gmtx.unlock(); - - //sync->wait(); - //UBLOG(logDEBUG5, "Calculator::initConnectors() - sendTransmitterDataSize - end"); - ////////////////////////////////////////////////////////////////////////// - //receiveTransmitterDataSize - //wenn er hier bei verteilten berechnungen stopped, dann ist vermutlich auf einer seite ein nicht aktiver block!!! - //UBLOG(logDEBUG5, "Calculator::initConnectors() - receiveTransmitterDataSize - start"); - - //sync->gmtx.lock(); - for (int l = minLevel; l <= maxLevel; l++) - { - if (l != maxLevel) - { - UBLOG(logDEBUG5, "Calculator::initRemoteConnectors()-receiveTransmitterDataSize(remoteInterConns["<<l<<"]).size ="<<remoteInterConns[l].size()); - BOOST_FOREACH(Block3DConnectorPtr c, remoteInterConns[l] ) c->receiveTransmitterDataSize(); - } - } - //sync->gmtx.unlock(); - //sync->wait(); - //UBLOG(logDEBUG5, "Calculator::initConnectors() - receiveTransmitterDataSize - end"); - ////////////////////////////////////////////////////////////////////////// -} -////////////////////////////////////////////////////////////////////////// -void MTCalculator::initConnectors(std::vector<Block3DConnectorPtr>& connectors) -{ - UBLOG(logDEBUG1, "Calculator::initConnectors() - start"); - - //initialization - ////////////////////////////////////////////////////////////////////////// - //initialize connectors - UBLOG(logDEBUG5, "Calculator::initConnectors() - connectoren initialisieren - start"); - BOOST_FOREACH(Block3DConnectorPtr c, connectors ) c->init(); - UBLOG(logDEBUG5, "Calculator::initConnectors() - connectoren initialisieren - end"); - ////////////////////////////////////////////////////////////////////////// - //sendTransmitterDataSize - UBLOG(logDEBUG5, "Calculator::initConnectors() - sendTransmitterDataSize - start"); - BOOST_FOREACH(Block3DConnectorPtr c, connectors ) c->sendTransmitterDataSize(); - UBLOG(logDEBUG5, "Calculator::initConnectors() - sendTransmitterDataSize - end"); - ////////////////////////////////////////////////////////////////////////// - //receiveTransmitterDataSize - //wenn er hier bei verteilten berechnungen stopped, dann ist vermutlich auf einer seite ein nicht aktiver block!!! - UBLOG(logDEBUG5, "Calculator::initConnectors() - receiveTransmitterDataSize - start"); - BOOST_FOREACH(Block3DConnectorPtr c, connectors ) c->receiveTransmitterDataSize(); - UBLOG(logDEBUG5, "Calculator::initConnectors() - receiveTransmitterDataSize - end"); - - UBLOG(logDEBUG1, "Calculator::initConnectors() - end"); -} -////////////////////////////////////////////////////////////////////////// -void MTCalculator::exchangeBlockData(int startLevel, int maxInitLevel, bool invStep) -{ - UBLOG(logDEBUG5, "Calculator::exchangeBlockData() - start"); - //sync->wait(); - //startLevel bis maxInitLevel - for(int level=startLevel; level<=maxInitLevel; level++) - { - connectorsPrepare(localConns[level]); - connectorsPrepare(remoteConns[level]); - - connectorsSend(localConns[level], invStep); - connectorsSend(remoteConns[level], invStep); - - connectorsReceive(localConns[level], invStep); - connectorsReceive(remoteConns[level], invStep); - } - //sync->wait(); - UBLOG(logDEBUG5, "Calculator::exchangeBlockData() - end"); -} -////////////////////////////////////////////////////////////////////////// -void MTCalculator::exchangeInterfaceBlockData(int startLevel, int maxInitLevel, bool invStep) -{ - sync->wait(); - //startLevel bis maxInitLevel - for(int level=startLevel; level<=maxInitLevel; level++) - { - connectorsPrepare(localInterfaceBlockConns[level]); - connectorsPrepare(remoteInterfaceBlockConns[level]); - - connectorsSend(localInterfaceBlockConns[level], invStep); - connectorsSend(remoteInterfaceBlockConns[level], invStep); - - connectorsReceive(localInterfaceBlockConns[level], invStep); - connectorsReceive(remoteInterfaceBlockConns[level], invStep); - } - sync->wait(); -} -////////////////////////////////////////////////////////////////////////// -void MTCalculator::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 MTCalculator::connectorsPrepare(std::vector< Block3DConnectorPtr >& connectors) -{ - UBLOG(logDEBUG5, "Calculator::connectorsPrepare() - start"); - BOOST_FOREACH(Block3DConnectorPtr c, connectors) - { - UBLOG(logDEBUG5, "Calculator::connectorsPrepare() - c->prepareForReceive() - start"); - c->prepareForReceive(); - UBLOG(logDEBUG5, "Calculator::connectorsPrepare() - c->prepareForReceive() - end"); - UBLOG(logDEBUG5, "Calculator::connectorsPrepare() - c->prepareForSend() - start"); - c->prepareForSend(); - UBLOG(logDEBUG5, "Calculator::connectorsPrepare() - c->prepareForSend() - end"); - } - UBLOG(logDEBUG5, "Calculator::connectorsPrepare() - end"); -} -////////////////////////////////////////////////////////////////////////// -void MTCalculator::connectorsSend(std::vector< Block3DConnectorPtr >& connectors, bool invStep) -{ - UBLOG(logDEBUG5, "Calculator::connectorsSend() - start"); - BOOST_FOREACH(Block3DConnectorPtr c, connectors) - { - c->setInvStep(invStep); - c->fillSendVectors(); - c->sendVectors(); - } - UBLOG(logDEBUG5, "Calculator::connectorsSend() - end"); -} -////////////////////////////////////////////////////////////////////////// -void MTCalculator::connectorsReceive(std::vector< Block3DConnectorPtr >& connectors, bool invStep) -{ - UBLOG(logDEBUG5, "Calculator::connectorsReceive() - start"); - BOOST_FOREACH(Block3DConnectorPtr c, connectors) - { - c->setInvStep(invStep); - c->receiveVectors(); - c->distributeReceiveVectors(); - } - UBLOG(logDEBUG5, "Calculator::connectorsReceive() - end"); -} -////////////////////////////////////////////////////////////////////////// -void MTCalculator::connectorsSetInvStep(std::vector< Block3DConnectorPtr >& connectors, bool invStep) -{ - BOOST_FOREACH(Block3DConnectorPtr c, connectors) - { - c->setInvStep(invStep); - } -} -////////////////////////////////////////////////////////////////////////// -void MTCalculator::interpolation(int startLevel, int maxInitLevel) -{ - UBLOG(logDEBUG5, "Calculator::interpolation() - start"); - //sync->wait(); - - for(int level=startLevel; level<maxInitLevel; level++) - { - connectorsPrepare(localInterConns[level]); - connectorsPrepare(remoteInterConns[level]); - } - - //sync->wait(); - - for(int level=startLevel; level<maxInitLevel; level++) - { - connectorsSend(localInterConns[level], true); - connectorsSend(remoteInterConns[level], true); - } - - //sync->wait(); - - for(int level=startLevel; level<maxInitLevel; level++) - { - connectorsReceive(localInterConns[level], true); - connectorsReceive(remoteInterConns[level], true); - } - - //sync->wait(); - UBLOG(logDEBUG5, "Calculator::interpolation() - end"); -} -////////////////////////////////////////////////////////////////////////// -void MTCalculator::setVisScheduler(UbSchedulerPtr s) -{ - visScheduler = s; -} -////////////////////////////////////////////////////////////////////////// -//double Calculator::getCallculationTime() -//{ -// return timer.getTotalTime(); -//} -////////////////////////////////////////////////////////////////////////// -std::vector< std::vector< Block3DPtr > > MTCalculator::getBlocks() -{ - return blocks; -} -////////////////////////////////////////////////////////////////////////// -void MTCalculator::deleteBlocks() -{ - BOOST_FOREACH(std::vector< Block3DPtr > &bs, blocks) - bs.resize(0); -} -////////////////////////////////////////////////////////////////////////// -void MTCalculator::deleteConnectors() -{ - deleteConnectors(localConns); - deleteConnectors(remoteConns); - - deleteConnectors(localInterfaceBlockConns); - deleteConnectors(remoteInterfaceBlockConns); - - deleteConnectors(localInterConns); - deleteConnectors(remoteInterConns); -} -////////////////////////////////////////////////////////////////////////// -void MTCalculator::deleteConnectors(std::vector< std::vector< Block3DConnectorPtr > >& conns) -{ - BOOST_FOREACH(std::vector< Block3DConnectorPtr > &c, conns) - c.resize(0); -} -////////////////////////////////////////////////////////////////////////// -void MTCalculator::applyBCs( int startLevel, int maxInitLevel ) -{ - //startLevel bis maxInitLevel - for(int level=startLevel; level<=maxInitLevel; level++) - { - //call LBM kernel - BOOST_FOREACH(Block3DPtr block, blocks[level]) - { - block->getKernel()->getBCProcessor()->applyBC(); - } - } -} diff --git a/source/VirtualFluidsCore/Grid/MTCalculator.h b/source/VirtualFluidsCore/Grid/MTCalculator.h deleted file mode 100644 index a594990ef..000000000 --- a/source/VirtualFluidsCore/Grid/MTCalculator.h +++ /dev/null @@ -1,82 +0,0 @@ -#ifndef MTCALCULATOR_H -#define MTCALCULATOR_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 "Calculator.h" - -class MTCalculator; -typedef boost::shared_ptr<MTCalculator> MTCalculatorPtr; - -#include "CalculationManager.h" - -class MTCalculator : public Calculator -{ -public: - MTCalculator(); - MTCalculator(Grid3DPtr grid, SynchronizerPtr sync, bool mainThread = true); - virtual ~MTCalculator(){} - virtual void calculate(const double& endTime, CalculationManagerPtr cm, boost::exception_ptr& error); - void addBlock(Block3DPtr block); - void initConnectors(); - void setVisScheduler(UbSchedulerPtr s); - //double getCallculationTime(); - std::vector< std::vector< Block3DPtr > > getBlocks(); - void deleteBlocks(); - void deleteConnectors(); -protected: - void calculateBlocks(int startLevel, int maxInitLevel); - void calculateBlocks(int minInitLevel, int maxInitLevel, int staggeredStep); - void initConnectors(std::vector<Block3DConnectorPtr>& connectors); - virtual void initRemoteConnectors(); - void swapDistributions(int startLevel, int maxInitLevel); - virtual void exchangeBlockData(int startLevel, int maxInitLevel, bool invStep); - void exchangeInterfaceBlockData(int startLevel, int maxInitLevel, bool invStep); - virtual void connectorsPrepare(std::vector< Block3DConnectorPtr >& connectors); - virtual void connectorsSend(std::vector< Block3DConnectorPtr >& connectors, bool invStep); - virtual void connectorsReceive(std::vector< Block3DConnectorPtr >& connectors, bool invStep); - void connectorsSetInvStep(std::vector< Block3DConnectorPtr >& connectors, bool invStep); - void interpolation(int startLevel, int maxInitLevel); - void deleteConnectors(std::vector< std::vector< Block3DConnectorPtr > >& conns); - void applyBCs(int startLevel, int maxInitLevel); - int minLevel, maxLevel; - std::vector< std::vector< Block3DConnectorPtr > > localConns; - std::vector< std::vector< Block3DConnectorPtr > > remoteConns; - SynchronizerPtr sync; - - boost::barrier* bar; - //double time; - - bool mainThread; - bool refinement; - Grid3DPtr grid; - UbSchedulerPtr visScheduler; - int calcStep; - -private: - std::vector< std::vector<Block3DPtr> > blocks; - - std::vector< std::vector< Block3DConnectorPtr > > localInterfaceBlockConns; - std::vector< std::vector< Block3DConnectorPtr > > remoteInterfaceBlockConns; - - //localInterConns and remoteInterConns save interpolation connectors - //every element save CF connectors for current level and FC connectors for next level - //e.g. - //localInterConns[0] = CF(0), FC(1) - //localInterConns[1] = CF(1), FC(2) - //localInterConns[2] - std::vector< std::vector< Block3DConnectorPtr > > localInterConns; - std::vector< std::vector< Block3DConnectorPtr > > remoteInterConns; - - //UbTimer timer, timer2, timer3; - bool loadBalancingComp; - - -}; - -#endif diff --git a/source/VirtualFluidsCore/Grid/PrePostBcCalculator.cpp b/source/VirtualFluidsCore/Grid/PrePostBcCalculator.cpp index 4abdfbd44..3eaddbe33 100644 --- a/source/VirtualFluidsCore/Grid/PrePostBcCalculator.cpp +++ b/source/VirtualFluidsCore/Grid/PrePostBcCalculator.cpp @@ -16,11 +16,6 @@ PrePostBcCalculator::PrePostBcCalculator(Grid3DPtr grid, SynchronizerPtr sync, b Calculator(grid, sync, mainThread) { -} -////////////////////////////////////////////////////////////////////////// -PrePostBcCalculator::PrePostBcCalculator(Grid3DPtr grid, SynchronizerPtr sync, BoundaryConditionProcessorPtr bcProcessor, bool mainThread) : -Calculator(grid, sync, bcProcessor, mainThread) -{ } ////////////////////////////////////////////////////////////////////////// void PrePostBcCalculator::calculate(const double& endTime, CalculationManagerPtr cm, boost::exception_ptr& error) diff --git a/source/VirtualFluidsCore/Grid/PrePostBcCalculator.h b/source/VirtualFluidsCore/Grid/PrePostBcCalculator.h index 0b81c061c..3caeb94da 100644 --- a/source/VirtualFluidsCore/Grid/PrePostBcCalculator.h +++ b/source/VirtualFluidsCore/Grid/PrePostBcCalculator.h @@ -8,7 +8,6 @@ #include "basics/utilities/UbScheduler.h" #include "basics/utilities/UbTiming.h" #include "LoadBalancer.h" -#include "BoundaryConditionProcessor.h" #include "PrePostBcCalculator.h" @@ -22,7 +21,6 @@ class PrePostBcCalculator : public Calculator public: PrePostBcCalculator(); PrePostBcCalculator(Grid3DPtr grid, SynchronizerPtr sync, bool mainThread = true); - PrePostBcCalculator(Grid3DPtr grid, SynchronizerPtr sync, BoundaryConditionProcessorPtr bcProcessor, bool mainThread = true); virtual ~PrePostBcCalculator(){} virtual void calculate(const double& endTime, CalculationManagerPtr cm, boost::exception_ptr& error); protected: diff --git a/source/VirtualFluidsCore/Utilities/ConfigurationFile.hpp b/source/VirtualFluidsCore/Utilities/ConfigurationFile.hpp index 03f84c5c2..473d9a518 100644 --- a/source/VirtualFluidsCore/Utilities/ConfigurationFile.hpp +++ b/source/VirtualFluidsCore/Utilities/ConfigurationFile.hpp @@ -56,10 +56,10 @@ public: bool getBool(const std::string& key) const; std::string getString(const std::string& key) const; template<class T> - std::vector<T> getVector(const std::string& Key) const; + std::vector<T> getVector(const std::string& key) const; template<class T> - T get(const std::string& Key) const; + T getValue(const std::string& key) const; private: // the container @@ -81,7 +81,7 @@ void ConfigurationFile::clear() { data.clear(); } - +////////////////////////////////////////////////////////////////////////// bool ConfigurationFile::load(const std::string& file) { std::ifstream inFile(file.c_str()); @@ -127,12 +127,12 @@ bool ConfigurationFile::load(const std::string& file) return true; } - +////////////////////////////////////////////////////////////////////////// bool ConfigurationFile::contains(const std::string& key) const { return data.find(key) != data.end(); } - +////////////////////////////////////////////////////////////////////////// std::string ConfigurationFile::getString(const std::string& key) const { std::map<std::string, std::string>::const_iterator iter = data.find(key); @@ -147,42 +147,42 @@ std::string ConfigurationFile::getString(const std::string& key) const UB_THROW(UbException(UB_EXARGS, "The parameter \"" + key + "\" is missing!")); } } - +////////////////////////////////////////////////////////////////////////// int ConfigurationFile::getInt(const std::string& key) const { std::string str = getString(key); int value = std::atoi(str.c_str()); return value; } - +////////////////////////////////////////////////////////////////////////// long ConfigurationFile::getLong(const std::string& key) const { std::string str = getString(key); long value = std::atol(str.c_str()); return value; } - +////////////////////////////////////////////////////////////////////////// float ConfigurationFile::getFloat(const std::string& key) const { std::string str = getString(key); float value = (float)std::atof(str.c_str()); return value; } - +////////////////////////////////////////////////////////////////////////// double ConfigurationFile::getDouble(const std::string& key) const { std::string str = getString(key); double value = std::atof(str.c_str()); return value; } - +////////////////////////////////////////////////////////////////////////// bool ConfigurationFile::getBool(const std::string& key) const { std::string str = getString(key); bool value = (str == "true"); return value; } - +////////////////////////////////////////////////////////////////////////// std::string ConfigurationFile::trim(const std::string& str) { size_t first = str.find_first_not_of(" \t\n\r"); @@ -198,7 +198,7 @@ std::string ConfigurationFile::trim(const std::string& str) return ""; } } - +////////////////////////////////////////////////////////////////////////// template<class T> std::vector<T> ConfigurationFile::getVector(const std::string& key) const { @@ -215,7 +215,7 @@ std::vector<T> ConfigurationFile::getVector(const std::string& key) const } return v; } - +////////////////////////////////////////////////////////////////////////// template<class T> T ConfigurationFile::fromString(const std::string& str) const { @@ -230,5 +230,29 @@ T ConfigurationFile::fromString(const std::string& str) const stream >> t; return t; } +////////////////////////////////////////////////////////////////////////// +template<class T> +T ConfigurationFile::getValue(const std::string& key) const +{ + std::string str = getString(key); + bool bFlag = false; + if ((std::string)typeid(T).name() == (std::string)typeid(bool).name()) + { + bFlag = true; + } + + std::istringstream iss(str); + T x; + iss >> x; + if (!iss && !bFlag) + UB_THROW(UbException(UB_EXARGS, " cannot convert \"" + str + "\" to type <" + static_cast<std::string>(typeid(x).name()) + ">")); + + if (bFlag) + { + bool value = (str == "true"); + x = value; + } + return x; +} #endif // Configuration_h__ diff --git a/source/VirtualFluidsCore/Visitors/BoundaryConditionBlockVisitor.cpp b/source/VirtualFluidsCore/Visitors/BoundaryConditionBlockVisitor.cpp index 45ac0f046..b2ae5725b 100644 --- a/source/VirtualFluidsCore/Visitors/BoundaryConditionBlockVisitor.cpp +++ b/source/VirtualFluidsCore/Visitors/BoundaryConditionBlockVisitor.cpp @@ -94,25 +94,21 @@ void BoundaryConditionBlockVisitor::visit(Grid3DPtr grid, Block3DPtr block) } if (vCount > 0) { - //velocity->addNnode(vCount); velocity->addDistributions(distributions); velocity->setCollFactor(collFactor); } if (dCount > 0) { - //density->addNnode(dCount); density->addDistributions(distributions); density->setCollFactor(collFactor); } if (nsCount > 0) { - //noSlip->addNnode(nsCount); noSlip->addDistributions(distributions); noSlip->setCollFactor(collFactor); } if (sCount > 0) { - //slip->addNnode(sCount); slip->addDistributions(distributions); slip->setCollFactor(collFactor); } diff --git a/source/VirtualFluidsCore/Visitors/BoundaryConditionBlockVisitor.h b/source/VirtualFluidsCore/Visitors/BoundaryConditionBlockVisitor.h index b15c70ed7..abb3366fc 100644 --- a/source/VirtualFluidsCore/Visitors/BoundaryConditionBlockVisitor.h +++ b/source/VirtualFluidsCore/Visitors/BoundaryConditionBlockVisitor.h @@ -2,7 +2,6 @@ #define BoundaryConditionBlockVisitor_h__ #include "Block3DVisitor.h" -#include "BoundaryConditionProcessor.h" class BoundaryConditionBlockVisitor : public Block3DVisitor { -- GitLab