diff --git a/source/Applications/DLR-F16-Solid/f16-solid.cfg b/source/Applications/DLR-F16-Solid/f16-solid.cfg index 87fdc6f55b90e34cb81e700beae0a0e9565c5558..da8cb35998b1e0a1616625c1597ffe6e0c450a61 100644 --- a/source/Applications/DLR-F16-Solid/f16-solid.cfg +++ b/source/Applications/DLR-F16-Solid/f16-solid.cfg @@ -5,12 +5,10 @@ pathGeo = d:/Projects/SFB880/DLR-F16/Geometry #fngFileWhole = grundgeometrie-direkter-export.stl fngFileWhole = F16withTapeMeshedCoarse.stl -zigZagTape = 2zackenbaender0.stl - pathReInit = /work/koskuche/DLR-F16_L7 stepReInit = 10000 -numOfThreads = 4 +numOfThreads = 4 availMem = 15e9 logToFile = false @@ -54,15 +52,15 @@ refineDistance = 0.6 writeBlocks = true newStart = true -restartStep = 10000 +restartStep = 100 cpStep = 10000 cpStart = 10000 -outTimeStep = 100 +outTimeStep = 1000 outTimeStart = 0 -endTime = 10000 +endTime = 100000 #Cp pcpStart = 1000000 diff --git a/source/Applications/DLR-F16-Solid/f16.cpp b/source/Applications/DLR-F16-Solid/f16.cpp index f7d8633fb4f420e4f3c51d4e8f6cbcb58dc99e40..7b31c4024c9d5770c2f884af5ccf7fe4c7a1803e 100644 --- a/source/Applications/DLR-F16-Solid/f16.cpp +++ b/source/Applications/DLR-F16-Solid/f16.cpp @@ -16,7 +16,6 @@ void run(string configname) string pathOut = config.getValue<string>("pathOut"); string pathGeo = config.getValue<string>("pathGeo"); string fngFileWhole = config.getValue<string>("fngFileWhole"); - string zigZagTape = config.getValue<string>("zigZagTape"); int numOfThreads = config.getValue<int>("numOfThreads"); vector<int> blockNx = config.getVector<int>("blockNx"); vector<double> boundingBox = config.getVector<double>("boundingBox"); @@ -146,13 +145,17 @@ void run(string configname) SPtr<BCAdapter> outflowBCAdapter(new DensityBCAdapter(rhoLB)); outflowBCAdapter->setBcAlgorithm(SPtr<BCAlgorithm>(new NonReflectingOutflowBCAlgorithm())); - BoundaryConditionsBlockVisitor bcVisitor; + BoundaryConditionsBlockVisitor bcVisitor; bcVisitor.addBC(noSlipBCAdapter); bcVisitor.addBC(velBCAdapter); bcVisitor.addBC(outflowBCAdapter); - SPtr<LBMKernel> kernel = SPtr<LBMKernel>(new CompressibleCumulantViscosity4thLBMKernel(blockNx[0], blockNx[1], blockNx[2], CompressibleCumulantViscosity4thLBMKernel::NORMAL)); - //SPtr<LBMKernel> kernel = SPtr<LBMKernel>(new CompressibleCumulantLBMKernel(blockNx[0], blockNx[1], blockNx[2], CompressibleCumulantLBMKernel::NORMAL)); + //SPtr<LBMKernel> kernel = SPtr<LBMKernel>(new CompressibleCumulantLBMKernel()); + //dynamicPointerCast<CompressibleCumulantLBMKernel>(kernel)->setRelaxationParameter(CompressibleCumulantLBMKernel::NORMAL); + + SPtr<LBMKernel> kernel = SPtr<LBMKernel>(new CompressibleCumulant4thOrderViscosityLBMKernel()); + + kernel->setNX(std::array<int,3>{{blockNx[0], blockNx[1], blockNx[2]}}); SPtr<BCProcessor> bcProc; bcProc = SPtr<BCProcessor>(new BCProcessor()); kernel->setBCProcessor(bcProc); @@ -453,7 +456,7 @@ void run(string configname) } //initialization of distributions - InitDistributionsBlockVisitor initVisitor1(nuLB, rhoLB); + InitDistributionsBlockVisitor initVisitor1; initVisitor1.setVx1(fct); grid->accept(initVisitor1); @@ -492,10 +495,14 @@ void run(string configname) ////sponge layer //////////////////////////////////////////////////////////////////////////// - //GbCuboid3DPtr spongeLayerX1max(new GbCuboid3D(g_maxX1-8.0*blockLength, g_minX2-blockLength, g_minX3-blockLength, g_maxX1+blockLength, g_maxX2+blockLength, g_maxX3+blockLength)); - //if (myid==0) GbSystem3D::writeGeoObject(spongeLayerX1max.get(), pathOut+"/geo/spongeLayerX1max", WbWriterVtkXmlASCII::getInstance()); - //SpongeLayerBlockVisitor slVisitorX1max(spongeLayerX1max, 1.0); - //grid->accept(slVisitorX1max); + GbCuboid3DPtr spongeLayerX1max(new GbCuboid3D(g_maxX1-8.0*blockLength, g_minX2-blockLength, g_minX3-blockLength, g_maxX1+blockLength, g_maxX2+blockLength, g_maxX3+blockLength)); + if (myid==0) GbSystem3D::writeGeoObject(spongeLayerX1max.get(), pathOut+"/geo/spongeLayerX1max", WbWriterVtkXmlASCII::getInstance()); + SpongeLayerBlockVisitor slVisitorX1max(spongeLayerX1max, 1.0); + SPtr<LBMKernel> spKernel = SPtr<LBMKernel>(new CompressibleCumulantLBMKernel()); + dynamicPointerCast<CompressibleCumulantLBMKernel>(spKernel)->setRelaxationParameter(CompressibleCumulantLBMKernel::NORMAL); + spKernel->setBCProcessor(bcProc); + slVisitorX1max.setKernel(spKernel); + grid->accept(slVisitorX1max); //GbCuboid3DPtr spongeLayerX1min(new GbCuboid3D(g_minX1-blockLength, g_minX2-blockLength, g_minX3-blockLength, g_minX1+75, g_maxX2+blockLength, g_maxX3+blockLength)); //if (myid==0) GbSystem3D::writeGeoObject(spongeLayerX1min.get(), pathOut+"/geo/spongeLayerX1min", WbWriterVtkXmlASCII::getInstance()); @@ -556,7 +563,7 @@ void run(string configname) //TimeseriesCoProcessor tsp1(grid, stepMV, mic1, pathOut+"/mic/mic1", comm); omp_set_num_threads(numOfThreads); - SPtr<Calculator> calculator(new OMPCalculator()); + SPtr<Calculator> calculator(new BasicCalculator()); calculator->setGrid(grid); calculator->setLastTimeStep(endTime); calculator->setVisScheduler(stepSch); @@ -614,7 +621,7 @@ void test_run() double g_maxX2 = 5; double g_maxX3 = 5; - double blockNx[3] ={ 5, 5, 5 }; + int blockNx[3] ={ 5, 5, 5 }; string pathOut = "d:/temp/DLR-F16-Solid-test"; @@ -641,7 +648,8 @@ void test_run() WriteBlocksCoProcessor ppblocks(grid, SPtr<UbScheduler>(new UbScheduler(1)), pathOut, WbWriterVtkXmlBinary::getInstance(), comm); ppblocks.process(0); - SPtr<LBMKernel> kernel = SPtr<LBMKernel>(new CompressibleCumulantViscosity4thLBMKernel(blockNx[0], blockNx[1], blockNx[2], CompressibleCumulantViscosity4thLBMKernel::NORMAL)); + SPtr<LBMKernel> kernel = SPtr<LBMKernel>(new CompressibleCumulant4thOrderViscosityLBMKernel()); + kernel->setNX(std::array<int,3>{{blockNx[0], blockNx[1], blockNx[2]}}); SPtr<BCProcessor> bcProc; bcProc = SPtr<BCProcessor>(new BCProcessor()); kernel->setBCProcessor(bcProc); @@ -650,7 +658,7 @@ void test_run() grid->accept(kernelVisitor); //initialization of distributions - InitDistributionsBlockVisitor initVisitor1(nuLB, rhoLB); + InitDistributionsBlockVisitor initVisitor1; initVisitor1.setVx1(0.001); grid->accept(initVisitor1); @@ -658,7 +666,7 @@ void test_run() SPtr<WriteMacroscopicQuantitiesCoProcessor> writeMQCoProcessor(new WriteMacroscopicQuantitiesCoProcessor(grid, stepSch, pathOut, WbWriterVtkXmlBinary::getInstance(), conv, comm)); //omp_set_num_threads(numOfThreads); - SPtr<Calculator> calculator(new MPICalculator()); + SPtr<Calculator> calculator(new BasicCalculator()); calculator->setGrid(grid); calculator->setLastTimeStep(2); calculator->setVisScheduler(stepSch); diff --git a/source/CMake/compilerflags/icc170.cmake b/source/CMake/compilerflags/icc170.cmake index 5b69a4249a2a1f0e29281f212d3bd3429d93c91b..6b22171c70d08170b73e47450936d9a73bbd0f00 100644 --- a/source/CMake/compilerflags/icc170.cmake +++ b/source/CMake/compilerflags/icc170.cmake @@ -22,6 +22,10 @@ MACRO(SET_COMPILER_SPECIFIC_FLAGS_INTERN build_type use64BitOptions) #LIST(APPEND CAB_COMPILER_ADDTIONAL_CXX_COMPILER_FLAGS "-xHOST -O3 -ip -ipo -fno-alias -mcmodel=medium -qopt-streaming-stores=always") LIST(APPEND CAB_COMPILER_ADDTIONAL_CXX_COMPILER_FLAGS "-xHOST -O3 -ip -fno-alias -mcmodel=medium -qopt-streaming-stores=always") + + #Debug + #LIST(APPEND CAB_COMPILER_ADDTIONAL_CXX_COMPILER_FLAGS "-g -traceback") + ############################################################################################################### ## OpenMP support ############################################################################################################### diff --git a/source/VirtualFluids.h b/source/VirtualFluids.h index c4e51224bfb70429944f6e1dad73311f3846cc34..0e0574ec8dd2f88d2f23f576be482cda1042a143 100644 --- a/source/VirtualFluids.h +++ b/source/VirtualFluids.h @@ -129,8 +129,7 @@ #include <Grid/Block3D.h> #include <Grid/Calculator.h> -#include <Grid/OMPCalculator.h> -#include <Grid/MPICalculator.h> +#include <Grid/BasicCalculator.h> #include <Grid/Grid3D.h> #include <Grid/Grid3DSystem.h> @@ -182,7 +181,7 @@ #include <LBM/LBMKernel.h> #include <LBM/IncompressibleCumulantLBMKernel.h> #include <LBM/CompressibleCumulantLBMKernel.h> -#include <LBM/CompressibleCumulantViscosity4thLBMKernel.h> +#include <LBM/CompressibleCumulant4thOrderViscosityLBMKernel.h> #include <LBM/InitDensityLBMKernel.h> #include <LBM/VoidLBMKernel.h> #include <LBM/LBMSystem.h> diff --git a/source/VirtualFluidsCore/CoProcessors/MPIIOMigrationCoProcessor.h b/source/VirtualFluidsCore/CoProcessors/MPIIOMigrationCoProcessor.h index c3248aaec649367ecfe79dc5585d1101e268f42b..e7d61504d5368dfd089d30df74514ffa6fb0cdb6 100644 --- a/source/VirtualFluidsCore/CoProcessors/MPIIOMigrationCoProcessor.h +++ b/source/VirtualFluidsCore/CoProcessors/MPIIOMigrationCoProcessor.h @@ -15,6 +15,7 @@ class LBMKernel; //! \class MPIWriteBlocksCoProcessor //! \brief Writes the grid each timestep into the files and reads the grip from the files before regenerating +//! \author Alena Karanchuk class MPIIOMigrationCoProcessor: public CoProcessor { //! \struct GridParam diff --git a/source/VirtualFluidsCore/CoProcessors/MPIIORestartCoProcessor.h b/source/VirtualFluidsCore/CoProcessors/MPIIORestartCoProcessor.h index 14b11757d1118f1ff0d7f4d3d54a585e87445b5b..f6758f9d882c2590c94495052e947df16d9789dd 100644 --- a/source/VirtualFluidsCore/CoProcessors/MPIIORestartCoProcessor.h +++ b/source/VirtualFluidsCore/CoProcessors/MPIIORestartCoProcessor.h @@ -15,6 +15,7 @@ class LBMKernel; //! \class MPIWriteBlocksCoProcessor //! \brief Writes the grid each timestep into the files and reads the grip from the files before regenerating +//! \author Alena Karanchuk class MPIIORestartCoProcessor: public CoProcessor { //! \struct GridParam @@ -64,7 +65,7 @@ class MPIIORestartCoProcessor: public CoProcessor //! \details The structure used to store some parameters needed to restore dataSet arrays struct dataSetParam { - int nx1; // to find the right block + int nx1; int nx2; int nx3; int nx[9][4]; // 9 arrays x (nx1, nx2, nx3, nx4) @@ -78,7 +79,7 @@ class MPIIORestartCoProcessor: public CoProcessor { double collFactor; double deltaT; - int x1; + int x1; // to find the right block int x2; int x3; int level; diff --git a/source/VirtualFluidsCore/Grid/OMPCalculator.cpp b/source/VirtualFluidsCore/Grid/BasicCalculator.cpp similarity index 89% rename from source/VirtualFluidsCore/Grid/OMPCalculator.cpp rename to source/VirtualFluidsCore/Grid/BasicCalculator.cpp index aee44aa60f0d47716f9b8b824e67d6c55352d549..e854ff93216528b396d66dbb378b8a4fdc616c30 100644 --- a/source/VirtualFluidsCore/Grid/OMPCalculator.cpp +++ b/source/VirtualFluidsCore/Grid/BasicCalculator.cpp @@ -1,36 +1,31 @@ -#include "OMPCalculator.h" +#include "BasicCalculator.h" -#include "Grid3D.h" #include "Block3D.h" #include "BCProcessor.h" #include "LBMKernel.h" #include "Block3DConnector.h" - -#include "MathUtil.hpp" #include "UbScheduler.h" -#include "UbTiming.h" -#include <UbException.h> +#include "UbLogger.h" #ifdef _OPENMP - #include <omp.h> +#include <omp.h> #endif #define OMP_SCHEDULE dynamic //#define TIMING +//#include "UbTiming.h" -#include "Block3DConnector.h" - -OMPCalculator::OMPCalculator() +BasicCalculator::BasicCalculator() { } ////////////////////////////////////////////////////////////////////////// -OMPCalculator::~OMPCalculator() +BasicCalculator::~BasicCalculator() { } ////////////////////////////////////////////////////////////////////////// -void OMPCalculator::calculate() +void BasicCalculator::calculate() { UBLOG(logDEBUG1, "OMPCalculator::calculate() - started"); int calcStep = 0; @@ -143,7 +138,7 @@ void OMPCalculator::calculate() } } ////////////////////////////////////////////////////////////////////////// -void OMPCalculator::calculateBlocks(int startLevel, int maxInitLevel, int calcStep) +void BasicCalculator::calculateBlocks(int startLevel, int maxInitLevel, int calcStep) { #ifdef _OPENMP #pragma omp parallel @@ -180,7 +175,7 @@ void OMPCalculator::calculateBlocks(int startLevel, int maxInitLevel, int calcSt } } ////////////////////////////////////////////////////////////////////////// -void OMPCalculator::exchangeBlockData(int startLevel, int maxInitLevel) +void BasicCalculator::exchangeBlockData(int startLevel, int maxInitLevel) { //startLevel bis maxInitLevel for (int level = startLevel; level<=maxInitLevel; level++) @@ -195,7 +190,7 @@ void OMPCalculator::exchangeBlockData(int startLevel, int maxInitLevel) } } ////////////////////////////////////////////////////////////////////////// -void OMPCalculator::swapDistributions(int startLevel, int maxInitLevel) +void BasicCalculator::swapDistributions(int startLevel, int maxInitLevel) { #ifdef _OPENMP #pragma omp parallel @@ -216,7 +211,7 @@ void OMPCalculator::swapDistributions(int startLevel, int maxInitLevel) } } ////////////////////////////////////////////////////////////////////////// -void OMPCalculator::connectorsPrepareLocal(std::vector< SPtr<Block3DConnector> >& connectors) +void BasicCalculator::connectorsPrepareLocal(std::vector< SPtr<Block3DConnector> >& connectors) { int size = (int)connectors.size(); #ifdef _OPENMP @@ -229,7 +224,7 @@ void OMPCalculator::connectorsPrepareLocal(std::vector< SPtr<Block3DConnector> > } } ////////////////////////////////////////////////////////////////////////// -void OMPCalculator::connectorsSendLocal(std::vector< SPtr<Block3DConnector> >& connectors) +void BasicCalculator::connectorsSendLocal(std::vector< SPtr<Block3DConnector> >& connectors) { int size = (int)connectors.size(); #ifdef _OPENMP @@ -242,7 +237,7 @@ void OMPCalculator::connectorsSendLocal(std::vector< SPtr<Block3DConnector> >& c } } ////////////////////////////////////////////////////////////////////////// -void OMPCalculator::connectorsReceiveLocal(std::vector< SPtr<Block3DConnector> >& connectors) +void BasicCalculator::connectorsReceiveLocal(std::vector< SPtr<Block3DConnector> >& connectors) { int size = (int)connectors.size(); #ifdef _OPENMP @@ -254,7 +249,7 @@ void OMPCalculator::connectorsReceiveLocal(std::vector< SPtr<Block3DConnector> > connectors[i]->distributeReceiveVectors(); } } -void OMPCalculator::connectorsPrepareRemote(std::vector< SPtr<Block3DConnector> >& connectors) +void BasicCalculator::connectorsPrepareRemote(std::vector< SPtr<Block3DConnector> >& connectors) { int size = (int)connectors.size(); for (int i =0; i<size; i++) @@ -264,7 +259,7 @@ void OMPCalculator::connectorsPrepareRemote(std::vector< SPtr<Block3DConnector> } } ////////////////////////////////////////////////////////////////////////// -void OMPCalculator::connectorsSendRemote(std::vector< SPtr<Block3DConnector> >& connectors) +void BasicCalculator::connectorsSendRemote(std::vector< SPtr<Block3DConnector> >& connectors) { int size = (int)connectors.size(); for (int i =0; i<size; i++) @@ -274,7 +269,7 @@ void OMPCalculator::connectorsSendRemote(std::vector< SPtr<Block3DConnector> >& } } ////////////////////////////////////////////////////////////////////////// -void OMPCalculator::connectorsReceiveRemote(std::vector< SPtr<Block3DConnector> >& connectors) +void BasicCalculator::connectorsReceiveRemote(std::vector< SPtr<Block3DConnector> >& connectors) { int size = (int)connectors.size(); for (int i =0; i<size; i++) @@ -284,7 +279,7 @@ void OMPCalculator::connectorsReceiveRemote(std::vector< SPtr<Block3DConnector> } } ////////////////////////////////////////////////////////////////////////// -void OMPCalculator::interpolation(int startLevel, int maxInitLevel) +void BasicCalculator::interpolation(int startLevel, int maxInitLevel) { for (int level = startLevel; level<maxInitLevel; level++) { @@ -305,7 +300,7 @@ void OMPCalculator::interpolation(int startLevel, int maxInitLevel) } } ////////////////////////////////////////////////////////////////////////// -void OMPCalculator::applyPreCollisionBC(int startLevel, int maxInitLevel) +void BasicCalculator::applyPreCollisionBC(int startLevel, int maxInitLevel) { //startLevel bis maxInitLevel for (int level = startLevel; level<=maxInitLevel; level++) @@ -321,7 +316,7 @@ void OMPCalculator::applyPreCollisionBC(int startLevel, int maxInitLevel) } } ////////////////////////////////////////////////////////////////////////// -void OMPCalculator::applyPostCollisionBC(int startLevel, int maxInitLevel) +void BasicCalculator::applyPostCollisionBC(int startLevel, int maxInitLevel) { //startLevel bis maxInitLevel for (int level = startLevel; level<=maxInitLevel; level++) diff --git a/source/VirtualFluidsCore/Grid/OMPCalculator.h b/source/VirtualFluidsCore/Grid/BasicCalculator.h similarity index 73% rename from source/VirtualFluidsCore/Grid/OMPCalculator.h rename to source/VirtualFluidsCore/Grid/BasicCalculator.h index 0a2c65cf9cd6975da30985d53ddda2c5788861dc..ee55d69f07bf5e4ab5778955f8b2639749707d87 100644 --- a/source/VirtualFluidsCore/Grid/OMPCalculator.h +++ b/source/VirtualFluidsCore/Grid/BasicCalculator.h @@ -1,15 +1,19 @@ -#ifndef OMPCALCULATOR_H -#define OMPCALCULATOR_H +#ifndef BasicCalculator_h__ +#define BasicCalculator_h__ #include "Calculator.h" class Block3DConnector; -class OMPCalculator : public Calculator +//! \class BasicCalculator +//! \brief Class implements basic functionality with MPI + OpenMP parallelization for main calculation loop +//! \author Konstantin Kutscher + +class BasicCalculator : public Calculator { public: - OMPCalculator(); - virtual ~OMPCalculator(); + BasicCalculator(); + virtual ~BasicCalculator(); virtual void calculate(); protected: @@ -28,5 +32,7 @@ protected: private: }; -#endif +#endif // BasicCalculator_h__ + + diff --git a/source/VirtualFluidsCore/Grid/Calculator.h b/source/VirtualFluidsCore/Grid/Calculator.h index 5ce5b3421995ee32efd45b0e33b8660dcb798022..0c33b74025a80fe244bc19bf9499887c230920cf 100644 --- a/source/VirtualFluidsCore/Grid/Calculator.h +++ b/source/VirtualFluidsCore/Grid/Calculator.h @@ -10,6 +10,10 @@ class Block3D; class Block3DConnector; class CoProcessor; +//! \class Calculator +//! \brief Abstract class for main calculation loop +//! \author Konstantin Kutscher + class Calculator { public: diff --git a/source/VirtualFluidsCore/Grid/MPICalculator.cpp b/source/VirtualFluidsCore/Grid/MPICalculator.cpp deleted file mode 100644 index 523d126165c5ec55ec2fd9b4b6ed7eeab3a8426e..0000000000000000000000000000000000000000 --- a/source/VirtualFluidsCore/Grid/MPICalculator.cpp +++ /dev/null @@ -1,268 +0,0 @@ -#include "MPICalculator.h" - -#include "Grid3D.h" -#include "Block3D.h" -#include "BCProcessor.h" -#include "LBMKernel.h" -#include "Block3DConnector.h" - -#include "MathUtil.hpp" -#include "UbScheduler.h" -#include "UbTiming.h" -#include <UbException.h> - - -//#define TIMING - -MPICalculator::MPICalculator() -{ - -} -////////////////////////////////////////////////////////////////////////// -MPICalculator::~MPICalculator() -{ - -} -////////////////////////////////////////////////////////////////////////// -void MPICalculator::calculate() -{ - UBLOG(logDEBUG1, "MPICalculator::calculate() - started"); - int calcStep = 0; - try - { - int minInitLevel = minLevel; - int maxInitLevel = maxLevel-minLevel; - int straightStartLevel = minInitLevel; - int internalIterations = 1<<(maxInitLevel-minInitLevel); - int forwardStartLevel; - int threshold; - -#ifdef TIMING - UbTimer timer; - double time[6]; -#endif - - for (calcStep = startTimeStep; calcStep<=lastTimeStep+1; calcStep++) - { - coProcess((double)(calcStep-1)); - - ////////////////////////////////////////////////////////////////////////// -#ifdef TIMING - UBLOG(logINFO, "calcStep = "<<calcStep); -#endif - ////////////////////////////////////////////////////////////////////////// - - for (int staggeredStep = 1; staggeredStep<=internalIterations; staggeredStep++) - { - forwardStartLevel = straightStartLevel; - if (staggeredStep==internalIterations) straightStartLevel = minInitLevel; - else - { - for (straightStartLevel = maxInitLevel, threshold = 1; - (staggeredStep&threshold)!=threshold; straightStartLevel--, threshold <<= 1); - } - ////////////////////////////////////////////////////////////////////////// -#ifdef TIMING - timer.resetAndStart(); -#endif - ////////////////////////////////////////////////////////////////////////// - applyPreCollisionBC(straightStartLevel, maxInitLevel); - - //do collision for all blocks - calculateBlocks(straightStartLevel, maxInitLevel, calcStep); - ////////////////////////////////////////////////////////////////////////// -#ifdef TIMING - time[0] = timer.stop(); - UBLOG(logINFO, "calculateBlocks time = "<<time[0]); -#endif - ////////////////////////////////////////////////////////////////////////// - ////////////////////////////////////////////////////////////////////////// - //exchange data between blocks - exchangeBlockData(straightStartLevel, maxInitLevel); - ////////////////////////////////////////////////////////////////////////// -#ifdef TIMING - time[1] = timer.stop(); - UBLOG(logINFO, "exchangeBlockData time = "<<time[1]); -#endif - ////////////////////////////////////////////////////////////////////////// - applyPostCollisionBC(straightStartLevel, maxInitLevel); - ////////////////////////////////////////////////////////////////////////// -#ifdef TIMING - time[2] = timer.stop(); - UBLOG(logINFO, "applyBCs time = "<<time[2]); -#endif - ////////////////////////////////////////////////////////////////////////// - //swap distributions in kernel - swapDistributions(straightStartLevel, maxInitLevel); - ////////////////////////////////////////////////////////////////////////// -#ifdef TIMING - time[3] = timer.stop(); - UBLOG(logINFO, "swapDistributions time = "<<time[3]); -#endif - ////////////////////////////////////////////////////////////////////////// - if (refinement) - { - if (straightStartLevel<maxInitLevel) - exchangeBlockData(straightStartLevel, maxInitLevel); - ////////////////////////////////////////////////////////////////////////// -#ifdef TIMING - time[4] = timer.stop(); - UBLOG(logINFO, "refinement exchangeBlockData time = "<<time[4]); -#endif - ////////////////////////////////////////////////////////////////////////// - //now ghost nodes have actual values - //interpolation of interface nodes between grid levels - interpolation(straightStartLevel, maxInitLevel); - ////////////////////////////////////////////////////////////////////////// -#ifdef TIMING - time[5] = timer.stop(); - UBLOG(logINFO, "refinement interpolation time = "<<time[5]); -#endif - ////////////////////////////////////////////////////////////////////////// - } - } - //exchange data between blocks for visualization - if ((int)visScheduler->getNextDueTime()==calcStep) - { - exchangeBlockData(straightStartLevel, maxInitLevel); - } - //now ghost nodes have actual values - } - UBLOG(logDEBUG1, "MPICalculator::calculate() - stoped"); - } - catch (std::exception& e) - { - UBLOG(logERROR, e.what()); - UBLOG(logERROR, " step = "<<calcStep); - } -} -////////////////////////////////////////////////////////////////////////// -void MPICalculator::calculateBlocks(int startLevel, int maxInitLevel, int calcStep) -{ - SPtr<Block3D> blockTemp; - try - { - //startLevel bis maxInitLevel - for (int level = startLevel; level<=maxInitLevel; level++) - { - //timer.resetAndStart(); - //call LBM kernel - for (SPtr<Block3D> 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) - { - UBLOG(logERROR, e.what()); - UBLOG(logERROR, blockTemp->toString()<<" step = "<<calcStep); - exit(EXIT_FAILURE); - } -} -////////////////////////////////////////////////////////////////////////// -void MPICalculator::exchangeBlockData(int startLevel, int maxInitLevel) -{ - //startLevel bis maxInitLevel - for (int level = startLevel; level<=maxInitLevel; level++) - { - connectorsPrepare(localConns[level]); - connectorsPrepare(remoteConns[level]); - - connectorsSend(localConns[level]); - connectorsSend(remoteConns[level]); - - connectorsReceive(localConns[level]); - connectorsReceive(remoteConns[level]); - } -} -////////////////////////////////////////////////////////////////////////// -void MPICalculator::swapDistributions(int startLevel, int maxInitLevel) -{ - //startLevel bis maxInitLevel - for (int level = startLevel; level<=maxInitLevel; level++) - { - for (SPtr<Block3D> block : blocks[level]) - { - block->getKernel()->swapDistributions(); - } - } -} -////////////////////////////////////////////////////////////////////////// -void MPICalculator::connectorsPrepare(std::vector< SPtr<Block3DConnector> >& connectors) -{ - for (SPtr<Block3DConnector> c : connectors) - { - c->prepareForReceive(); - c->prepareForSend(); - } -} -////////////////////////////////////////////////////////////////////////// -void MPICalculator::connectorsSend(std::vector< SPtr<Block3DConnector> >& connectors) -{ - for (SPtr<Block3DConnector> c : connectors) - { - c->fillSendVectors(); - c->sendVectors(); - } -} -////////////////////////////////////////////////////////////////////////// -void MPICalculator::connectorsReceive(std::vector< SPtr<Block3DConnector> >& connectors) -{ - for (SPtr<Block3DConnector> c : connectors) - { - c->receiveVectors(); - c->distributeReceiveVectors(); - } -} -////////////////////////////////////////////////////////////////////////// -void MPICalculator::interpolation(int startLevel, int maxInitLevel) -{ - - - for (int level = startLevel; level<maxInitLevel; level++) - { - connectorsPrepare(localInterConns[level]); - connectorsPrepare(remoteInterConns[level]); - } - - for (int level = startLevel; level<maxInitLevel; level++) - { - connectorsSend(localInterConns[level]); - connectorsSend(remoteInterConns[level]); - } - - for (int level = startLevel; level<maxInitLevel; level++) - { - connectorsReceive(localInterConns[level]); - connectorsReceive(remoteInterConns[level]); - } -} -////////////////////////////////////////////////////////////////////////// -void MPICalculator::applyPreCollisionBC(int startLevel, int maxInitLevel) -{ - //startLevel bis maxInitLevel - for (int level = startLevel; level<=maxInitLevel; level++) - { - for (SPtr<Block3D> block : blocks[level]) - { - block->getKernel()->getBCProcessor()->applyPreCollisionBC(); - } - } -} -////////////////////////////////////////////////////////////////////////// -void MPICalculator::applyPostCollisionBC(int startLevel, int maxInitLevel) -{ - //startLevel bis maxInitLevel - for (int level = startLevel; level<=maxInitLevel; level++) - { - for (SPtr<Block3D> block : blocks[level]) - { - block->getKernel()->getBCProcessor()->applyPostCollisionBC(); - } - } -} - diff --git a/source/VirtualFluidsCore/Grid/MPICalculator.h b/source/VirtualFluidsCore/Grid/MPICalculator.h deleted file mode 100644 index 23ed1f30feb258753ebfff2aee609c50b82e6ec7..0000000000000000000000000000000000000000 --- a/source/VirtualFluidsCore/Grid/MPICalculator.h +++ /dev/null @@ -1,31 +0,0 @@ -#ifndef MPICALCULATOR_H -#define MPICALCULATOR_H - -#include "Calculator.h" - -class Block3DConnector; - -class MPICalculator : public Calculator -{ -public: - MPICalculator(); - virtual ~MPICalculator(); - virtual void calculate(); - -protected: - void calculateBlocks(int startLevel, int maxInitLevel, int calcStep); - void swapDistributions(int startLevel, int maxInitLevel); - virtual void exchangeBlockData(int startLevel, int maxInitLevel); - virtual void connectorsPrepare(std::vector< SPtr<Block3DConnector> >& connectors); - virtual void connectorsSend(std::vector< SPtr<Block3DConnector> >& connectors); - virtual void connectorsReceive(std::vector< SPtr<Block3DConnector> >& connectors); - void interpolation(int startLevel, int maxInitLevel); - void applyPreCollisionBC(int startLevel, int maxInitLevel); - void applyPostCollisionBC(int startLevel, int maxInitLevel); -private: - - -}; - -#endif - diff --git a/source/VirtualFluidsCore/LBM/CompressibleCumulantViscosity4thLBMKernel.cpp b/source/VirtualFluidsCore/LBM/CompressibleCumulant4thOrderViscosityLBMKernel.cpp similarity index 93% rename from source/VirtualFluidsCore/LBM/CompressibleCumulantViscosity4thLBMKernel.cpp rename to source/VirtualFluidsCore/LBM/CompressibleCumulant4thOrderViscosityLBMKernel.cpp index 7e2406c942343cc4353eb1f962d543820a2d3516..63ab6e24b40cd78e3aed4bd81c991679b4e560c1 100644 --- a/source/VirtualFluidsCore/LBM/CompressibleCumulantViscosity4thLBMKernel.cpp +++ b/source/VirtualFluidsCore/LBM/CompressibleCumulant4thOrderViscosityLBMKernel.cpp @@ -1,4 +1,4 @@ -#include "CompressibleCumulantViscosity4thLBMKernel.h" +#include "CompressibleCumulant4thOrderViscosityLBMKernel.h" #include "D3Q27System.h" #include "InterpolationProcessor.h" #include "D3Q27EsoTwist3DSplittedVector.h" @@ -9,45 +9,27 @@ #define PROOF_CORRECTNESS ////////////////////////////////////////////////////////////////////////// -CompressibleCumulantViscosity4thLBMKernel::CompressibleCumulantViscosity4thLBMKernel() +CompressibleCumulant4thOrderViscosityLBMKernel::CompressibleCumulant4thOrderViscosityLBMKernel() { - this->nx1 = 0; - this->nx2 = 0; - this->nx3 = 0; - this->parameter = NORMAL; - this->OxyyMxzz = 1.0; this->compressible = true; - this->bulkOmegaToOmega = false; - this->OxxPyyPzz = 1.0; } ////////////////////////////////////////////////////////////////////////// -CompressibleCumulantViscosity4thLBMKernel::CompressibleCumulantViscosity4thLBMKernel(int nx1, int nx2, int nx3, Parameter p) -{ - this->nx1 = nx1; - this->nx2 = nx2; - this->nx3 = nx3; - this->parameter = p; - this->OxyyMxzz = 1.0; - this->compressible = true; - this->bulkOmegaToOmega = false; - this->OxxPyyPzz = 1.0; -} -////////////////////////////////////////////////////////////////////////// -CompressibleCumulantViscosity4thLBMKernel::~CompressibleCumulantViscosity4thLBMKernel(void) +CompressibleCumulant4thOrderViscosityLBMKernel::~CompressibleCumulant4thOrderViscosityLBMKernel(void) { } ////////////////////////////////////////////////////////////////////////// -void CompressibleCumulantViscosity4thLBMKernel::init() +void CompressibleCumulant4thOrderViscosityLBMKernel::initDataSet() { - SPtr<DistributionArray3D> d(new D3Q27EsoTwist3DSplittedVector(nx1+2, nx2+2, nx3+2, 0.0)); + SPtr<DistributionArray3D> d(new D3Q27EsoTwist3DSplittedVector(nx[0]+2, nx[1]+2, nx[2]+2, -999.9)); dataSet->setFdistributions(d); } ////////////////////////////////////////////////////////////////////////// -SPtr<LBMKernel> CompressibleCumulantViscosity4thLBMKernel::clone() +SPtr<LBMKernel> CompressibleCumulant4thOrderViscosityLBMKernel::clone() { - SPtr<LBMKernel> kernel(new CompressibleCumulantViscosity4thLBMKernel(nx1, nx2, nx3, parameter)); - dynamicPointerCast<CompressibleCumulantViscosity4thLBMKernel>(kernel)->init(); + SPtr<LBMKernel> kernel(new CompressibleCumulant4thOrderViscosityLBMKernel()); + kernel->setNX(nx); + dynamicPointerCast<CompressibleCumulant4thOrderViscosityLBMKernel>(kernel)->initDataSet(); kernel->setCollisionFactor(this->collFactor); kernel->setBCProcessor(bcProcessor->clone(kernel)); kernel->setWithForcing(withForcing); @@ -56,39 +38,26 @@ SPtr<LBMKernel> CompressibleCumulantViscosity4thLBMKernel::clone() kernel->setForcingX3(muForcingX3); kernel->setIndex(ix1, ix2, ix3); kernel->setDeltaT(deltaT); - switch (parameter) - { - case NORMAL: - dynamicPointerCast<CompressibleCumulantViscosity4thLBMKernel>(kernel)->OxyyMxzz = 1.0; - break; - case MAGIC: - dynamicPointerCast<CompressibleCumulantViscosity4thLBMKernel>(kernel)->OxyyMxzz = 2.0 +(-collFactor); - break; - } if (bulkOmegaToOmega) { - dynamicPointerCast<CompressibleCumulantViscosity4thLBMKernel>(kernel)->OxxPyyPzz = collFactor; + dynamicPointerCast<CompressibleCumulant4thOrderViscosityLBMKernel>(kernel)->OxxPyyPzz = collFactor; } else { - dynamicPointerCast<CompressibleCumulantViscosity4thLBMKernel>(kernel)->OxxPyyPzz = one; + dynamicPointerCast<CompressibleCumulant4thOrderViscosityLBMKernel>(kernel)->OxxPyyPzz = one; } + return kernel; } ////////////////////////////////////////////////////////////////////////// -void CompressibleCumulantViscosity4thLBMKernel::calculate() -{ - timer.resetAndStart(); - collideAll(); - timer.stop(); -} -////////////////////////////////////////////////////////////////////////// -void CompressibleCumulantViscosity4thLBMKernel::collideAll() +void CompressibleCumulant4thOrderViscosityLBMKernel::calculate() { using namespace D3Q27System; using namespace std; - + + //timer.resetAndStart(); + //initializing of forcing stuff if (withForcing) { @@ -583,16 +552,16 @@ void CompressibleCumulantViscosity4thLBMKernel::collideAll() //limiter-Scheise Teil 1 //LBMReal oxxyy,oxxzz,oxy,oxz,oyz; //LBMReal smag=0.001; - //oxxyy = omega+(one-omega)*abs(mxxMyy)/(abs(mxxMyy)+smag); - //oxxzz = omega+(one-omega)*abs(mxxMzz)/(abs(mxxMzz)+smag); - //oxy = omega+(one-omega)*abs(mfbba)/(abs(mfbba)+smag); - //oxz = omega+(one-omega)*abs(mfbab)/(abs(mfbab)+smag); - //oyz = omega+(one-omega)*abs(mfabb)/(abs(mfabb)+smag); + //oxxyy = omega+(one-omega)*fabs(mxxMyy)/(fabs(mxxMyy)+smag); + //oxxzz = omega+(one-omega)*fabs(mxxMzz)/(fabs(mxxMzz)+smag); + //oxy = omega+(one-omega)*fabs(mfbba)/(fabs(mfbba)+smag); + //oxz = omega+(one-omega)*fabs(mfbab)/(fabs(mfbab)+smag); + //oyz = omega+(one-omega)*fabs(mfabb)/(fabs(mfabb)+smag); //////////////////////////////////////////////////////////////////////////// ////Teil 1b //LBMReal constante = 1000.0; - //LBMReal nuEddi = constante * abs(mxxPyyPzz); + //LBMReal nuEddi = constante * fabs(mxxPyyPzz); //LBMReal omegaLimit = one / (one / omega + three * nuEddi); //{ @@ -675,19 +644,19 @@ void CompressibleCumulantViscosity4thLBMKernel::collideAll() //relax ////////////////////////////////////////////////////////////////////////// //das ist der limiter - wadjust = Oxyz+(one-Oxyz)*abs(mfbbb)/(abs(mfbbb)+qudricLimitD); + wadjust = Oxyz+(one-Oxyz)*fabs(mfbbb)/(fabs(mfbbb)+qudricLimitD); mfbbb += wadjust * (-mfbbb); - wadjust = OxyyPxzz+(one-OxyyPxzz)*abs(mxxyPyzz)/(abs(mxxyPyzz)+qudricLimitP); + wadjust = OxyyPxzz+(one-OxyyPxzz)*fabs(mxxyPyzz)/(fabs(mxxyPyzz)+qudricLimitP); mxxyPyzz += wadjust * (-mxxyPyzz); - wadjust = OxyyMxzz+(one-OxyyMxzz)*abs(mxxyMyzz)/(abs(mxxyMyzz)+qudricLimitM); + wadjust = OxyyMxzz+(one-OxyyMxzz)*fabs(mxxyMyzz)/(fabs(mxxyMyzz)+qudricLimitM); mxxyMyzz += wadjust * (-mxxyMyzz); - wadjust = OxyyPxzz+(one-OxyyPxzz)*abs(mxxzPyyz)/(abs(mxxzPyyz)+qudricLimitP); + wadjust = OxyyPxzz+(one-OxyyPxzz)*fabs(mxxzPyyz)/(fabs(mxxzPyyz)+qudricLimitP); mxxzPyyz += wadjust * (-mxxzPyyz); - wadjust = OxyyMxzz+(one-OxyyMxzz)*abs(mxxzMyyz)/(abs(mxxzMyyz)+qudricLimitM); + wadjust = OxyyMxzz+(one-OxyyMxzz)*fabs(mxxzMyyz)/(fabs(mxxzMyyz)+qudricLimitM); mxxzMyyz += wadjust * (-mxxzMyyz); - wadjust = OxyyPxzz+(one-OxyyPxzz)*abs(mxyyPxzz)/(abs(mxyyPxzz)+qudricLimitP); + wadjust = OxyyPxzz+(one-OxyyPxzz)*fabs(mxyyPxzz)/(fabs(mxyyPxzz)+qudricLimitP); mxyyPxzz += wadjust * (-mxyyPxzz); - wadjust = OxyyMxzz+(one-OxyyMxzz)*abs(mxyyMxzz)/(abs(mxyyMxzz)+qudricLimitM); + wadjust = OxyyMxzz+(one-OxyyMxzz)*fabs(mxyyMxzz)/(fabs(mxyyMxzz)+qudricLimitM); mxyyMxzz += wadjust * (-mxyyMxzz); ////////////////////////////////////////////////////////////////////////// //ohne limiter @@ -711,18 +680,18 @@ void CompressibleCumulantViscosity4thLBMKernel::collideAll() //4. ////////////////////////////////////////////////////////////////////////// //mit limiter - // wadjust = O4+(one-O4)*abs(CUMacc)/(abs(CUMacc)+qudricLimit); + // wadjust = O4+(one-O4)*fabs(CUMacc)/(fabs(CUMacc)+qudricLimit); //CUMacc += wadjust * (-CUMacc); - // wadjust = O4+(one-O4)*abs(CUMcac)/(abs(CUMcac)+qudricLimit); + // wadjust = O4+(one-O4)*fabs(CUMcac)/(fabs(CUMcac)+qudricLimit); //CUMcac += wadjust * (-CUMcac); - // wadjust = O4+(one-O4)*abs(CUMcca)/(abs(CUMcca)+qudricLimit); + // wadjust = O4+(one-O4)*fabs(CUMcca)/(fabs(CUMcca)+qudricLimit); //CUMcca += wadjust * (-CUMcca); - // wadjust = O4+(one-O4)*abs(CUMbbc)/(abs(CUMbbc)+qudricLimit); + // wadjust = O4+(one-O4)*fabs(CUMbbc)/(fabs(CUMbbc)+qudricLimit); //CUMbbc += wadjust * (-CUMbbc); - // wadjust = O4+(one-O4)*abs(CUMbcb)/(abs(CUMbcb)+qudricLimit); + // wadjust = O4+(one-O4)*fabs(CUMbcb)/(fabs(CUMbcb)+qudricLimit); //CUMbcb += wadjust * (-CUMbcb); - // wadjust = O4+(one-O4)*abs(CUMcbb)/(abs(CUMcbb)+qudricLimit); + // wadjust = O4+(one-O4)*fabs(CUMcbb)/(fabs(CUMcbb)+qudricLimit); //CUMcbb += wadjust * (-CUMcbb); ////////////////////////////////////////////////////////////////////////// ////////////////////////////////////////////////////////////////////////// @@ -1064,15 +1033,21 @@ void CompressibleCumulantViscosity4thLBMKernel::collideAll() } } } + //timer.stop(); } ////////////////////////////////////////////////////////////////////////// -double CompressibleCumulantViscosity4thLBMKernel::getCalculationTime() +double CompressibleCumulant4thOrderViscosityLBMKernel::getCalculationTime() { //return timer.getDuration(); return timer.getTotalTime(); } ////////////////////////////////////////////////////////////////////////// -void CompressibleCumulantViscosity4thLBMKernel::setBulkOmegaToOmega(bool value) +void CompressibleCumulant4thOrderViscosityLBMKernel::setBulkOmegaToOmega(bool value) { bulkOmegaToOmega = value; } + +//void CompressibleCumulant4thOrderViscosityLBMKernel::setViscosityFlag(bool vf) +//{ +// viscosityFlag = vf; +//} diff --git a/source/VirtualFluidsCore/LBM/CompressibleCumulantViscosity4thLBMKernel.h b/source/VirtualFluidsCore/LBM/CompressibleCumulant4thOrderViscosityLBMKernel.h similarity index 66% rename from source/VirtualFluidsCore/LBM/CompressibleCumulantViscosity4thLBMKernel.h rename to source/VirtualFluidsCore/LBM/CompressibleCumulant4thOrderViscosityLBMKernel.h index 4eefd8c3b62966bd03c8a893592be0a3ea463f12..96bd3288ed9bb39d08baca8edbbe5f99b8b81e11 100644 --- a/source/VirtualFluidsCore/LBM/CompressibleCumulantViscosity4thLBMKernel.h +++ b/source/VirtualFluidsCore/LBM/CompressibleCumulant4thOrderViscosityLBMKernel.h @@ -11,34 +11,25 @@ //! \brief compressible cumulant LBM kernel. //! \details CFD solver that use Cascaded Cumulant Lattice Boltzmann method for D3Q27 model //! \author K. Kutscher, M. Geier -class CompressibleCumulantViscosity4thLBMKernel : public LBMKernel +class CompressibleCumulant4thOrderViscosityLBMKernel : public LBMKernel { public: //! This option set relaxation parameter: NORMAL enum Parameter{NORMAL, MAGIC}; public: - CompressibleCumulantViscosity4thLBMKernel(); - //! Constructor - //! \param nx1 number of nodes in x dimension - //! \param nx2 number of nodes in y dimension - //! \param nx3 number of nodes in z dimension - //! \param p set relaxation parameter: NORMAL is OxyyMxzz = 1.0 and MAGIC is OxyyMxzz = 2.0 +(-collFactor) - CompressibleCumulantViscosity4thLBMKernel(int nx1, int nx2, int nx3, Parameter p); - virtual ~CompressibleCumulantViscosity4thLBMKernel(void); + CompressibleCumulant4thOrderViscosityLBMKernel(); + virtual ~CompressibleCumulant4thOrderViscosityLBMKernel(void); virtual void calculate(); virtual SPtr<LBMKernel> clone(); double getCalculationTime() override; void setBulkOmegaToOmega(bool value); + //void setViscosityFlag(bool vf); protected: - virtual void collideAll(); - virtual void init(); + virtual void initDataSet(); LBMReal f[D3Q27System::ENDF+1]; UbTimer timer; - LBMReal OxyyMxzz; - Parameter parameter; - CbArray4D<LBMReal,IndexerX4X3X2X1>::CbArray4DPtr localDistributions; CbArray4D<LBMReal,IndexerX4X3X2X1>::CbArray4DPtr nonLocalDistributions; CbArray3D<LBMReal,IndexerX3X2X1>::CbArray3DPtr zeroDistributions; @@ -53,6 +44,11 @@ protected: // bulk viscosity bool bulkOmegaToOmega; LBMReal OxxPyyPzz; + + //bool viscosityFlag; + //LBMReal OxyyPxzz; + //LBMReal OxyyMxzz; + //LBMReal Oxyz; }; #endif // CompressibleCumulantLBMKernel_h__ diff --git a/source/VirtualFluidsCore/LBM/CompressibleCumulantLBMKernel.cpp b/source/VirtualFluidsCore/LBM/CompressibleCumulantLBMKernel.cpp index 13ee2e1d8d12ff600d522934f8b5609b26718ce7..65fc41cd84bbea76326deb8396e6ec696e3eea9e 100644 --- a/source/VirtualFluidsCore/LBM/CompressibleCumulantLBMKernel.cpp +++ b/source/VirtualFluidsCore/LBM/CompressibleCumulantLBMKernel.cpp @@ -3,29 +3,14 @@ #include "InterpolationProcessor.h" #include "D3Q27EsoTwist3DSplittedVector.h" #include <math.h> -//#include <omp.h> #include "DataSet3D.h" #define PROOF_CORRECTNESS ////////////////////////////////////////////////////////////////////////// CompressibleCumulantLBMKernel::CompressibleCumulantLBMKernel() { - this->nx1 = 0; - this->nx2 = 0; - this->nx3 = 0; - this->parameter = NORMAL; - this->OxyyMxzz = 1.0; this->compressible = true; - this->bulkOmegaToOmega = false; - this->OxxPyyPzz = 1.0; -} -////////////////////////////////////////////////////////////////////////// -CompressibleCumulantLBMKernel::CompressibleCumulantLBMKernel(int nx1, int nx2, int nx3, Parameter p) -{ - this->nx1 = nx1; - this->nx2 = nx2; - this->nx3 = nx3; - this->parameter = p; + this->parameter = NORMAL; this->OxyyMxzz = 1.0; this->compressible = true; this->bulkOmegaToOmega = false; @@ -37,16 +22,17 @@ CompressibleCumulantLBMKernel::~CompressibleCumulantLBMKernel(void) } ////////////////////////////////////////////////////////////////////////// -void CompressibleCumulantLBMKernel::init() +void CompressibleCumulantLBMKernel::initDataSet() { - SPtr<DistributionArray3D> d(new D3Q27EsoTwist3DSplittedVector(nx1+2, nx2+2, nx3+2, -999.0)); + SPtr<DistributionArray3D> d(new D3Q27EsoTwist3DSplittedVector(nx[0]+2, nx[1]+2, nx[2]+2, -999.9)); dataSet->setFdistributions(d); } ////////////////////////////////////////////////////////////////////////// SPtr<LBMKernel> CompressibleCumulantLBMKernel::clone() { - SPtr<LBMKernel> kernel(new CompressibleCumulantLBMKernel(nx1, nx2, nx3, parameter)); - dynamicPointerCast<CompressibleCumulantLBMKernel>(kernel)->init(); + SPtr<LBMKernel> kernel(new CompressibleCumulantLBMKernel()); + kernel->setNX(nx); + dynamicPointerCast<CompressibleCumulantLBMKernel>(kernel)->initDataSet(); kernel->setCollisionFactor(this->collFactor); kernel->setBCProcessor(bcProcessor->clone(kernel)); kernel->setWithForcing(withForcing); @@ -71,23 +57,18 @@ SPtr<LBMKernel> CompressibleCumulantLBMKernel::clone() } else { - dynamicPointerCast<CompressibleCumulantLBMKernel>(kernel)->OxxPyyPzz = one; + dynamicPointerCast<CompressibleCumulantLBMKernel>(kernel)->OxxPyyPzz = one; } return kernel; } ////////////////////////////////////////////////////////////////////////// void CompressibleCumulantLBMKernel::calculate() -{ - timer.resetAndStart(); - collideAll(); - timer.stop(); -} -////////////////////////////////////////////////////////////////////////// -void CompressibleCumulantLBMKernel::collideAll() { using namespace D3Q27System; using namespace std; + //timer.resetAndStart(); + //initializing of forcing stuff if (withForcing) { @@ -247,10 +228,10 @@ void CompressibleCumulantLBMKernel::collideAll() vy2 = vvy*vvy; vz2 = vvz*vvz; //////////////////////////////////////////////////////////////////////////////////// - LBMReal wadjust; - LBMReal qudricLimitP = 0.01f;// * 0.0001f; - LBMReal qudricLimitM = 0.01f;// * 0.0001f; - LBMReal qudricLimitD = 0.01f;// * 0.001f; + //LBMReal wadjust; + //LBMReal qudricLimitP = 0.01f;// * 0.0001f; + //LBMReal qudricLimitM = 0.01f;// * 0.0001f; + //LBMReal qudricLimitD = 0.01f;// * 0.001f; //LBMReal s9 = minusomega; //test //s9 = 0.; @@ -579,16 +560,16 @@ void CompressibleCumulantLBMKernel::collideAll() //limiter-Scheise Teil 1 //LBMReal oxxyy,oxxzz,oxy,oxz,oyz; //LBMReal smag=0.001; - //oxxyy = omega+(one-omega)*abs(mxxMyy)/(abs(mxxMyy)+smag); - //oxxzz = omega+(one-omega)*abs(mxxMzz)/(abs(mxxMzz)+smag); - //oxy = omega+(one-omega)*abs(mfbba)/(abs(mfbba)+smag); - //oxz = omega+(one-omega)*abs(mfbab)/(abs(mfbab)+smag); - //oyz = omega+(one-omega)*abs(mfabb)/(abs(mfabb)+smag); + //oxxyy = omega+(one-omega)*fabs(mxxMyy)/(fabs(mxxMyy)+smag); + //oxxzz = omega+(one-omega)*fabs(mxxMzz)/(fabs(mxxMzz)+smag); + //oxy = omega+(one-omega)*fabs(mfbba)/(fabs(mfbba)+smag); + //oxz = omega+(one-omega)*fabs(mfbab)/(fabs(mfbab)+smag); + //oyz = omega+(one-omega)*fabs(mfabb)/(fabs(mfabb)+smag); //////////////////////////////////////////////////////////////////////////// ////Teil 1b //LBMReal constante = 1000.0; - //LBMReal nuEddi = constante * abs(mxxPyyPzz); + //LBMReal nuEddi = constante * fabs(mxxPyyPzz); //LBMReal omegaLimit = one / (one / omega + three * nuEddi); //{ @@ -665,19 +646,19 @@ void CompressibleCumulantLBMKernel::collideAll() //relax ////////////////////////////////////////////////////////////////////////// //das ist der limiter - //wadjust = Oxyz+(one-Oxyz)*abs(mfbbb)/(abs(mfbbb)+qudricLimitD); + //wadjust = Oxyz+(one-Oxyz)*fabs(mfbbb)/(fabs(mfbbb)+qudricLimitD); //mfbbb += wadjust * (-mfbbb); - //wadjust = OxyyPxzz+(one-OxyyPxzz)*abs(mxxyPyzz)/(abs(mxxyPyzz)+qudricLimitP); + //wadjust = OxyyPxzz+(one-OxyyPxzz)*fabs(mxxyPyzz)/(fabs(mxxyPyzz)+qudricLimitP); //mxxyPyzz += wadjust * (-mxxyPyzz); - //wadjust = OxyyMxzz+(one-OxyyMxzz)*abs(mxxyMyzz)/(abs(mxxyMyzz)+qudricLimitM); + //wadjust = OxyyMxzz+(one-OxyyMxzz)*fabs(mxxyMyzz)/(fabs(mxxyMyzz)+qudricLimitM); //mxxyMyzz += wadjust * (-mxxyMyzz); - //wadjust = OxyyPxzz+(one-OxyyPxzz)*abs(mxxzPyyz)/(abs(mxxzPyyz)+qudricLimitP); + //wadjust = OxyyPxzz+(one-OxyyPxzz)*fabs(mxxzPyyz)/(fabs(mxxzPyyz)+qudricLimitP); //mxxzPyyz += wadjust * (-mxxzPyyz); - //wadjust = OxyyMxzz+(one-OxyyMxzz)*abs(mxxzMyyz)/(abs(mxxzMyyz)+qudricLimitM); + //wadjust = OxyyMxzz+(one-OxyyMxzz)*fabs(mxxzMyyz)/(fabs(mxxzMyyz)+qudricLimitM); //mxxzMyyz += wadjust * (-mxxzMyyz); - //wadjust = OxyyPxzz+(one-OxyyPxzz)*abs(mxyyPxzz)/(abs(mxyyPxzz)+qudricLimitP); + //wadjust = OxyyPxzz+(one-OxyyPxzz)*fabs(mxyyPxzz)/(fabs(mxyyPxzz)+qudricLimitP); //mxyyPxzz += wadjust * (-mxyyPxzz); - //wadjust = OxyyMxzz+(one-OxyyMxzz)*abs(mxyyMxzz)/(abs(mxyyMxzz)+qudricLimitM); + //wadjust = OxyyMxzz+(one-OxyyMxzz)*fabs(mxyyMxzz)/(fabs(mxyyMxzz)+qudricLimitM); //mxyyMxzz += wadjust * (-mxyyMxzz); ////////////////////////////////////////////////////////////////////////// //ohne limiter @@ -701,18 +682,18 @@ void CompressibleCumulantLBMKernel::collideAll() //4. ////////////////////////////////////////////////////////////////////////// //mit limiter - // wadjust = O4+(one-O4)*abs(CUMacc)/(abs(CUMacc)+qudricLimit); + // wadjust = O4+(one-O4)*fabs(CUMacc)/(fabs(CUMacc)+qudricLimit); //CUMacc += wadjust * (-CUMacc); - // wadjust = O4+(one-O4)*abs(CUMcac)/(abs(CUMcac)+qudricLimit); + // wadjust = O4+(one-O4)*fabs(CUMcac)/(fabs(CUMcac)+qudricLimit); //CUMcac += wadjust * (-CUMcac); - // wadjust = O4+(one-O4)*abs(CUMcca)/(abs(CUMcca)+qudricLimit); + // wadjust = O4+(one-O4)*fabs(CUMcca)/(fabs(CUMcca)+qudricLimit); //CUMcca += wadjust * (-CUMcca); - // wadjust = O4+(one-O4)*abs(CUMbbc)/(abs(CUMbbc)+qudricLimit); + // wadjust = O4+(one-O4)*fabs(CUMbbc)/(fabs(CUMbbc)+qudricLimit); //CUMbbc += wadjust * (-CUMbbc); - // wadjust = O4+(one-O4)*abs(CUMbcb)/(abs(CUMbcb)+qudricLimit); + // wadjust = O4+(one-O4)*fabs(CUMbcb)/(fabs(CUMbcb)+qudricLimit); //CUMbcb += wadjust * (-CUMbcb); - // wadjust = O4+(one-O4)*abs(CUMcbb)/(abs(CUMcbb)+qudricLimit); + // wadjust = O4+(one-O4)*fabs(CUMcbb)/(fabs(CUMcbb)+qudricLimit); //CUMcbb += wadjust * (-CUMcbb); ////////////////////////////////////////////////////////////////////////// //ohne limiter @@ -1044,6 +1025,8 @@ void CompressibleCumulantLBMKernel::collideAll() } } + + //timer.stop(); } ////////////////////////////////////////////////////////////////////////// double CompressibleCumulantLBMKernel::getCalculationTime() @@ -1056,3 +1039,8 @@ void CompressibleCumulantLBMKernel::setBulkOmegaToOmega(bool value) { bulkOmegaToOmega = value; } +////////////////////////////////////////////////////////////////////////// +void CompressibleCumulantLBMKernel::setRelaxationParameter(Parameter p) +{ + parameter = p; +} diff --git a/source/VirtualFluidsCore/LBM/CompressibleCumulantLBMKernel.h b/source/VirtualFluidsCore/LBM/CompressibleCumulantLBMKernel.h index efffadcbe6acc3e5a572e5cc204199a2307f5ab1..35a52dc71471463cd4cd03cfa9304354bd5e9806 100644 --- a/source/VirtualFluidsCore/LBM/CompressibleCumulantLBMKernel.h +++ b/source/VirtualFluidsCore/LBM/CompressibleCumulantLBMKernel.h @@ -18,20 +18,14 @@ public: enum Parameter{NORMAL, MAGIC}; public: CompressibleCumulantLBMKernel(); - //! Constructor - //! \param nx1 number of nodes in x dimension - //! \param nx2 number of nodes in y dimension - //! \param nx3 number of nodes in z dimension - //! \param p set relaxation parameter: NORMAL is OxyyMxzz = 1.0 and MAGIC is OxyyMxzz = 2.0 +(-collFactor) - CompressibleCumulantLBMKernel(int nx1, int nx2, int nx3, Parameter p); virtual ~CompressibleCumulantLBMKernel(void); virtual void calculate(); virtual SPtr<LBMKernel> clone(); double getCalculationTime(); void setBulkOmegaToOmega(bool value); + void setRelaxationParameter(Parameter p); protected: - virtual void collideAll(); - virtual void init(); + virtual void initDataSet(); LBMReal f[D3Q27System::ENDF+1]; UbTimer timer; diff --git a/source/VirtualFluidsCore/LBM/IncompressibleCumulantLBMKernel.cpp b/source/VirtualFluidsCore/LBM/IncompressibleCumulantLBMKernel.cpp index 96b295477c61296694e40df8cf9d90bacb151e5a..168e37876700d382ef967255f21c55fe14f2b625 100644 --- a/source/VirtualFluidsCore/LBM/IncompressibleCumulantLBMKernel.cpp +++ b/source/VirtualFluidsCore/LBM/IncompressibleCumulantLBMKernel.cpp @@ -4,46 +4,33 @@ #include "D3Q27EsoTwist3DSplittedVector.h" #include "DataSet3D.h" #include <math.h> -//#include <omp.h> #define PROOF_CORRECTNESS ////////////////////////////////////////////////////////////////////////// IncompressibleCumulantLBMKernel::IncompressibleCumulantLBMKernel() { - this->nx1 = 0; - this->nx2 = 0; - this->nx3 = 0; this->parameter = NORMAL; this->OxyyMxzz = 1.0; this->compressible = false; } ////////////////////////////////////////////////////////////////////////// -IncompressibleCumulantLBMKernel::IncompressibleCumulantLBMKernel(int nx1, int nx2, int nx3, Parameter p) -{ - this->nx1 = nx1; - this->nx2 = nx2; - this->nx3 = nx3; - parameter = p; - this->compressible = false; -} -////////////////////////////////////////////////////////////////////////// IncompressibleCumulantLBMKernel::~IncompressibleCumulantLBMKernel(void) { } ////////////////////////////////////////////////////////////////////////// -void IncompressibleCumulantLBMKernel::init() +void IncompressibleCumulantLBMKernel::initDataSet() { - //SPtr<DistributionArray3D> d(new D3Q27EsoTwist3DSplittedVector(nx1+ghostLayerWitdh*2, nx2+ghostLayerWitdh*2, nx3+ghostLayerWitdh*2, -999.0)); - SPtr<DistributionArray3D> d(new D3Q27EsoTwist3DSplittedVector(nx1+2, nx2+2, nx3+2, -999.0)); + SPtr<DistributionArray3D> d(new D3Q27EsoTwist3DSplittedVector(nx[0]+2, nx[1]+2, nx[2]+2, -999.9)); dataSet->setFdistributions(d); } ////////////////////////////////////////////////////////////////////////// SPtr<LBMKernel> IncompressibleCumulantLBMKernel::clone() { - SPtr<LBMKernel> kernel(new IncompressibleCumulantLBMKernel(nx1, nx2, nx3, parameter)); - dynamicPointerCast<IncompressibleCumulantLBMKernel>(kernel)->init(); + SPtr<LBMKernel> kernel(new IncompressibleCumulantLBMKernel()); + kernel->setNX(nx); + dynamicPointerCast<IncompressibleCumulantLBMKernel>(kernel)->initDataSet(); kernel->setCollisionFactor(this->collFactor); kernel->setBCProcessor(bcProcessor->clone(kernel)); kernel->setWithForcing(withForcing); @@ -56,7 +43,7 @@ SPtr<LBMKernel> IncompressibleCumulantLBMKernel::clone() { case NORMAL: dynamicPointerCast<IncompressibleCumulantLBMKernel>(kernel)->OxyyMxzz = 1.0; - break; + break; case MAGIC: dynamicPointerCast<IncompressibleCumulantLBMKernel>(kernel)->OxyyMxzz = 2.0 +(-collFactor); break; @@ -65,35 +52,30 @@ SPtr<LBMKernel> IncompressibleCumulantLBMKernel::clone() } ////////////////////////////////////////////////////////////////////////// void IncompressibleCumulantLBMKernel::calculate() -{ - timer.resetAndStart(); - collideAll(); - timer.stop(); -} -////////////////////////////////////////////////////////////////////////// -void IncompressibleCumulantLBMKernel::collideAll() { using namespace D3Q27System; using namespace std; + //timer.resetAndStart(); + //initializing of forcing stuff if (withForcing) { - muForcingX1.DefineVar("x1",&muX1); muForcingX1.DefineVar("x2",&muX2); muForcingX1.DefineVar("x3",&muX3); - muForcingX2.DefineVar("x1",&muX1); muForcingX2.DefineVar("x2",&muX2); muForcingX2.DefineVar("x3",&muX3); - muForcingX3.DefineVar("x1",&muX1); muForcingX3.DefineVar("x2",&muX2); muForcingX3.DefineVar("x3",&muX3); + muForcingX1.DefineVar("x1", &muX1); muForcingX1.DefineVar("x2", &muX2); muForcingX1.DefineVar("x3", &muX3); + muForcingX2.DefineVar("x1", &muX1); muForcingX2.DefineVar("x2", &muX2); muForcingX2.DefineVar("x3", &muX3); + muForcingX3.DefineVar("x1", &muX1); muForcingX3.DefineVar("x2", &muX2); muForcingX3.DefineVar("x3", &muX3); muDeltaT = deltaT; - muForcingX1.DefineVar("dt",&muDeltaT); - muForcingX2.DefineVar("dt",&muDeltaT); - muForcingX3.DefineVar("dt",&muDeltaT); + muForcingX1.DefineVar("dt", &muDeltaT); + muForcingX2.DefineVar("dt", &muDeltaT); + muForcingX3.DefineVar("dt", &muDeltaT); muNu = (1.0/3.0)*(1.0/collFactor - 1.0/2.0); - muForcingX1.DefineVar("nu",&muNu); - muForcingX2.DefineVar("nu",&muNu); - muForcingX3.DefineVar("nu",&muNu); + muForcingX1.DefineVar("nu", &muNu); + muForcingX2.DefineVar("nu", &muNu); + muForcingX3.DefineVar("nu", &muNu); LBMReal forcingX1 = 0; LBMReal forcingX2 = 0; @@ -118,20 +100,13 @@ void IncompressibleCumulantLBMKernel::collideAll() int maxX2 = bcArrayMaxX2-ghostLayerWidth; int maxX3 = bcArrayMaxX3-ghostLayerWidth; - -//#pragma omp parallel num_threads(8) + for (int x3 = minX3; x3 < maxX3; x3++) { - // int i = omp_get_thread_num(); - // printf_s("Hello from thread %d\n", i); - //} -//#pragma omp for - for(int x3 = minX3; x3 < maxX3; x3++) - { - for(int x2 = minX2; x2 < maxX2; x2++) + for (int x2 = minX2; x2 < maxX2; x2++) { - for(int x1 = minX1; x1 < maxX1; x1++) + for (int x1 = minX1; x1 < maxX1; x1++) { - if(!bcArray->isSolid(x1,x2,x3) && !bcArray->isUndefined(x1,x2,x3)) + if (!bcArray->isSolid(x1, x2, x3) && !bcArray->isUndefined(x1, x2, x3)) { int x1p = x1 + 1; int x2p = x2 + 1; @@ -153,39 +128,39 @@ void IncompressibleCumulantLBMKernel::collideAll() //a - negative //b - null //c - positive - + // a b c //-1 0 1 - LBMReal mfcbb = (*this->localDistributions)(D3Q27System::ET_E, x1,x2,x3); - LBMReal mfbcb = (*this->localDistributions)(D3Q27System::ET_N,x1,x2,x3); - LBMReal mfbbc = (*this->localDistributions)(D3Q27System::ET_T,x1,x2,x3); - LBMReal mfccb = (*this->localDistributions)(D3Q27System::ET_NE,x1,x2,x3); - LBMReal mfacb = (*this->localDistributions)(D3Q27System::ET_NW,x1p,x2,x3); - LBMReal mfcbc = (*this->localDistributions)(D3Q27System::ET_TE,x1,x2,x3); - LBMReal mfabc = (*this->localDistributions)(D3Q27System::ET_TW, x1p,x2,x3); - LBMReal mfbcc = (*this->localDistributions)(D3Q27System::ET_TN,x1,x2,x3); - LBMReal mfbac = (*this->localDistributions)(D3Q27System::ET_TS,x1,x2p,x3); - LBMReal mfccc = (*this->localDistributions)(D3Q27System::ET_TNE,x1,x2,x3); - LBMReal mfacc = (*this->localDistributions)(D3Q27System::ET_TNW,x1p,x2,x3); - LBMReal mfcac = (*this->localDistributions)(D3Q27System::ET_TSE,x1,x2p,x3); - LBMReal mfaac = (*this->localDistributions)(D3Q27System::ET_TSW,x1p,x2p,x3); - - LBMReal mfabb = (*this->nonLocalDistributions)(D3Q27System::ET_W,x1p,x2,x3 ); - LBMReal mfbab = (*this->nonLocalDistributions)(D3Q27System::ET_S,x1,x2p,x3 ); - LBMReal mfbba = (*this->nonLocalDistributions)(D3Q27System::ET_B,x1,x2,x3p ); - LBMReal mfaab = (*this->nonLocalDistributions)(D3Q27System::ET_SW,x1p,x2p,x3 ); - LBMReal mfcab = (*this->nonLocalDistributions)(D3Q27System::ET_SE,x1,x2p,x3 ); - LBMReal mfaba = (*this->nonLocalDistributions)(D3Q27System::ET_BW,x1p,x2,x3p ); - LBMReal mfcba = (*this->nonLocalDistributions)(D3Q27System::ET_BE,x1,x2,x3p ); - LBMReal mfbaa = (*this->nonLocalDistributions)(D3Q27System::ET_BS,x1,x2p,x3p ); - LBMReal mfbca = (*this->nonLocalDistributions)(D3Q27System::ET_BN,x1,x2,x3p ); - LBMReal mfaaa = (*this->nonLocalDistributions)(D3Q27System::ET_BSW,x1p,x2p,x3p); - LBMReal mfcaa = (*this->nonLocalDistributions)(D3Q27System::ET_BSE,x1,x2p,x3p); - LBMReal mfaca = (*this->nonLocalDistributions)(D3Q27System::ET_BNW,x1p,x2,x3p); - LBMReal mfcca = (*this->nonLocalDistributions)(D3Q27System::ET_BNE,x1,x2,x3p); - - LBMReal mfbbb = (*this->zeroDistributions)(x1,x2,x3); + LBMReal mfcbb = (*this->localDistributions)(D3Q27System::ET_E, x1, x2, x3); + LBMReal mfbcb = (*this->localDistributions)(D3Q27System::ET_N, x1, x2, x3); + LBMReal mfbbc = (*this->localDistributions)(D3Q27System::ET_T, x1, x2, x3); + LBMReal mfccb = (*this->localDistributions)(D3Q27System::ET_NE, x1, x2, x3); + LBMReal mfacb = (*this->localDistributions)(D3Q27System::ET_NW, x1p, x2, x3); + LBMReal mfcbc = (*this->localDistributions)(D3Q27System::ET_TE, x1, x2, x3); + LBMReal mfabc = (*this->localDistributions)(D3Q27System::ET_TW, x1p, x2, x3); + LBMReal mfbcc = (*this->localDistributions)(D3Q27System::ET_TN, x1, x2, x3); + LBMReal mfbac = (*this->localDistributions)(D3Q27System::ET_TS, x1, x2p, x3); + LBMReal mfccc = (*this->localDistributions)(D3Q27System::ET_TNE, x1, x2, x3); + LBMReal mfacc = (*this->localDistributions)(D3Q27System::ET_TNW, x1p, x2, x3); + LBMReal mfcac = (*this->localDistributions)(D3Q27System::ET_TSE, x1, x2p, x3); + LBMReal mfaac = (*this->localDistributions)(D3Q27System::ET_TSW, x1p, x2p, x3); + + LBMReal mfabb = (*this->nonLocalDistributions)(D3Q27System::ET_W, x1p, x2, x3); + LBMReal mfbab = (*this->nonLocalDistributions)(D3Q27System::ET_S, x1, x2p, x3); + LBMReal mfbba = (*this->nonLocalDistributions)(D3Q27System::ET_B, x1, x2, x3p); + LBMReal mfaab = (*this->nonLocalDistributions)(D3Q27System::ET_SW, x1p, x2p, x3); + LBMReal mfcab = (*this->nonLocalDistributions)(D3Q27System::ET_SE, x1, x2p, x3); + LBMReal mfaba = (*this->nonLocalDistributions)(D3Q27System::ET_BW, x1p, x2, x3p); + LBMReal mfcba = (*this->nonLocalDistributions)(D3Q27System::ET_BE, x1, x2, x3p); + LBMReal mfbaa = (*this->nonLocalDistributions)(D3Q27System::ET_BS, x1, x2p, x3p); + LBMReal mfbca = (*this->nonLocalDistributions)(D3Q27System::ET_BN, x1, x2, x3p); + LBMReal mfaaa = (*this->nonLocalDistributions)(D3Q27System::ET_BSW, x1p, x2p, x3p); + LBMReal mfcaa = (*this->nonLocalDistributions)(D3Q27System::ET_BSE, x1, x2p, x3p); + LBMReal mfaca = (*this->nonLocalDistributions)(D3Q27System::ET_BNW, x1p, x2, x3p); + LBMReal mfcca = (*this->nonLocalDistributions)(D3Q27System::ET_BNE, x1, x2, x3p); + + LBMReal mfbbb = (*this->zeroDistributions)(x1, x2, x3); LBMReal m0, m1, m2; @@ -267,7 +242,7 @@ void IncompressibleCumulantLBMKernel::collideAll() m1 = mfaac - mfaaa; m0 = m2 + mfaab; mfaaa = m0; - m0 += c1o36 * oMdrho; + m0 += c1o36 * oMdrho; mfaab = m1 - m0 * vvz; mfaac = m2 - 2. * m1 * vvz + vz2 * m0; //////////////////////////////////////////////////////////////////////////////////// @@ -499,9 +474,9 @@ void IncompressibleCumulantLBMKernel::collideAll() //LBMReal CUMbcb = mfbcb - ((mfaca + c1o3 * oMdrho) * mfbab + 2. * mfbba * mfabb); // till 18.05.2015 //LBMReal CUMbbc = mfbbc - ((mfaac + c1o3 * oMdrho) * mfbba + 2. * mfbab * mfabb); // till 18.05.2015 - LBMReal CUMcbb = mfcbb - ((mfcaa + c1o3 ) * mfabb + 2. * mfbba * mfbab); - LBMReal CUMbcb = mfbcb - ((mfaca + c1o3 ) * mfbab + 2. * mfbba * mfabb); - LBMReal CUMbbc = mfbbc - ((mfaac + c1o3 ) * mfbba + 2. * mfbab * mfabb); + LBMReal CUMcbb = mfcbb - ((mfcaa + c1o3) * mfabb + 2. * mfbba * mfbab); + LBMReal CUMbcb = mfbcb - ((mfaca + c1o3) * mfbab + 2. * mfbba * mfabb); + LBMReal CUMbbc = mfbbc - ((mfaac + c1o3) * mfbba + 2. * mfbab * mfabb); LBMReal CUMcca = mfcca - ((mfcaa * mfaca + 2. * mfbba * mfbba) + c1o3 * (mfcaa + mfaca) * oMdrho + c1o9*(oMdrho-1)*oMdrho); LBMReal CUMcac = mfcac - ((mfcaa * mfaac + 2. * mfbab * mfbab) + c1o3 * (mfcaa + mfaac) * oMdrho + c1o9*(oMdrho-1)*oMdrho); @@ -513,17 +488,17 @@ void IncompressibleCumulantLBMKernel::collideAll() LBMReal CUMccb = mfccb - (mfcaa * mfacb + mfaca * mfcab + 4. * mfbba * mfbbb + 2. * (mfbab * mfbca + mfabb * mfcba)) - c1o3 * (mfacb + mfcab) * oMdrho; //Cum 6. - LBMReal CUMccc = mfccc +((-4. * mfbbb * mfbbb + LBMReal CUMccc = mfccc +((-4. * mfbbb * mfbbb - (mfcaa * mfacc + mfaca * mfcac + mfaac * mfcca) - 4. * (mfabb * mfcbb + mfbab * mfbcb + mfbba * mfbbc) - 2. * (mfbca * mfbac + mfcba * mfabc + mfcab * mfacb)) - +( 4. * (mfbab * mfbab * mfaca + mfabb * mfabb * mfcaa + mfbba * mfbba * mfaac) - + 2. * (mfcaa * mfaca * mfaac) - + 16. * mfbba * mfbab * mfabb) + +(4. * (mfbab * mfbab * mfaca + mfabb * mfabb * mfcaa + mfbba * mfbba * mfaac) + + 2. * (mfcaa * mfaca * mfaac) + + 16. * mfbba * mfbab * mfabb) - c1o3* (mfacc + mfcac + mfcca) * oMdrho -c1o9*oMdrho*oMdrho - c1o9* (mfcaa + mfaca + mfaac) * oMdrho*(1.-2.* oMdrho)- c1o27* oMdrho * oMdrho*(-2.* oMdrho) - +( 2. * (mfbab * mfbab + mfabb * mfabb + mfbba * mfbba) - + (mfaac * mfaca + mfaac * mfcaa + mfaca * mfcaa)) * c2o3*oMdrho) +c1o27*oMdrho; + +(2. * (mfbab * mfbab + mfabb * mfabb + mfbba * mfbba) + + (mfaac * mfaca + mfaac * mfcaa + mfaca * mfcaa)) * c2o3*oMdrho) +c1o27*oMdrho; //2. // linear combinations @@ -545,9 +520,9 @@ void IncompressibleCumulantLBMKernel::collideAll() mfbba += collFactor * (-mfbba); // linear combinations back - mfcaa = c1o3 * ( mxxMyy + mxxMzz + mxxPyyPzz); + mfcaa = c1o3 * (mxxMyy + mxxMzz + mxxPyyPzz); mfaca = c1o3 * (-2. * mxxMyy + mxxMzz + mxxPyyPzz); - mfaac = c1o3 * ( mxxMyy - 2. * mxxMzz + mxxPyyPzz); + mfaac = c1o3 * (mxxMyy - 2. * mxxMzz + mxxPyyPzz); //3. // linear combinations @@ -577,11 +552,11 @@ void IncompressibleCumulantLBMKernel::collideAll() mxyyMxzz += wadjust * (-mxyyMxzz); // linear combinations back - mfcba = ( mxxyMyzz + mxxyPyzz) * c1o2; + mfcba = (mxxyMyzz + mxxyPyzz) * c1o2; mfabc = (-mxxyMyzz + mxxyPyzz) * c1o2; - mfcab = ( mxxzMyyz + mxxzPyyz) * c1o2; + mfcab = (mxxzMyyz + mxxzPyyz) * c1o2; mfacb = (-mxxzMyyz + mxxzPyyz) * c1o2; - mfbca = ( mxyyMxzz + mxyyPxzz) * c1o2; + mfbca = (mxyyMxzz + mxyyPxzz) * c1o2; mfbac = (-mxyyMxzz + mxyyPxzz) * c1o2; //4. @@ -607,9 +582,9 @@ void IncompressibleCumulantLBMKernel::collideAll() //mfbcb = CUMbcb + ((mfaca + c1o3 * oMdrho) * mfbab + 2. * mfbba * mfabb); // till 18.05.2015 //mfbbc = CUMbbc + ((mfaac + c1o3 * oMdrho) * mfbba + 2. * mfbab * mfabb); // till 18.05.2015 - mfcbb = CUMcbb + ((mfcaa + c1o3 ) * mfabb + 2. * mfbba * mfbab); - mfbcb = CUMbcb + ((mfaca + c1o3 ) * mfbab + 2. * mfbba * mfabb); - mfbbc = CUMbbc + ((mfaac + c1o3 ) * mfbba + 2. * mfbab * mfabb); + mfcbb = CUMcbb + ((mfcaa + c1o3) * mfabb + 2. * mfbba * mfbab); + mfbcb = CUMbcb + ((mfaca + c1o3) * mfbab + 2. * mfbba * mfabb); + mfbbc = CUMbbc + ((mfaac + c1o3) * mfbba + 2. * mfbab * mfabb); mfcca = CUMcca + (mfcaa * mfaca + 2. * mfbba * mfbba) + c1o3 * (mfcaa + mfaca) * oMdrho + c1o9*(oMdrho-1)*oMdrho; mfcac = CUMcac + (mfcaa * mfaac + 2. * mfbab * mfbab) + c1o3 * (mfcaa + mfaac) * oMdrho + c1o9*(oMdrho-1)*oMdrho; @@ -621,17 +596,17 @@ void IncompressibleCumulantLBMKernel::collideAll() mfccb = CUMccb + (mfcaa * mfacb + mfaca * mfcab + 4. * mfbba * mfbbb + 2. * (mfbab * mfbca + mfabb * mfcba)) + c1o3 * (mfacb + mfcab) * oMdrho; //6. - mfccc = CUMccc -((-4. * mfbbb * mfbbb + mfccc = CUMccc -((-4. * mfbbb * mfbbb - (mfcaa * mfacc + mfaca * mfcac + mfaac * mfcca) - 4. * (mfabb * mfcbb + mfbac * mfbca + mfbba * mfbbc) - 2. * (mfbca * mfbac + mfcba * mfabc + mfcab * mfacb)) - +( 4. * (mfbab * mfbab * mfaca + mfabb * mfabb * mfcaa + mfbba * mfbba * mfaac) - + 2. * (mfcaa * mfaca * mfaac) - + 16. * mfbba * mfbab * mfabb) + +(4. * (mfbab * mfbab * mfaca + mfabb * mfabb * mfcaa + mfbba * mfbba * mfaac) + + 2. * (mfcaa * mfaca * mfaac) + + 16. * mfbba * mfbab * mfabb) - c1o3* (mfacc + mfcac + mfcca) * oMdrho -c1o9*oMdrho*oMdrho - c1o9* (mfcaa + mfaca + mfaac) * oMdrho*(1.-2.* oMdrho)- c1o27* oMdrho * oMdrho*(-2.* oMdrho) - +( 2. * (mfbab * mfbab + mfabb * mfabb + mfbba * mfbba) - + (mfaac * mfaca + mfaac * mfcaa + mfaca * mfcaa)) * c2o3*oMdrho) -c1o27*oMdrho; + +(2. * (mfbab * mfbab + mfabb * mfabb + mfbba * mfbba) + + (mfaac * mfaca + mfaac * mfcaa + mfaca * mfcaa)) * c2o3*oMdrho) -c1o27*oMdrho; //////////////////////////////////////////////////////////////////////////////////// //forcing @@ -646,67 +621,67 @@ void IncompressibleCumulantLBMKernel::collideAll() //mit 1, 0, 1/3, 0, 0, 0, 1/3, 0, 1/9 Konditionieren //////////////////////////////////////////////////////////////////////////////////// // Z - Dir - m0 = mfaac * c1o2 + mfaab * (vvz - c1o2) + (mfaaa + 1. * oMdrho) * ( vz2 - vvz) * c1o2; + m0 = mfaac * c1o2 + mfaab * (vvz - c1o2) + (mfaaa + 1. * oMdrho) * (vz2 - vvz) * c1o2; m1 = -mfaac - 2. * mfaab * vvz + mfaaa * (1. - vz2) - 1. * oMdrho * vz2; - m2 = mfaac * c1o2 + mfaab * (vvz + c1o2) + (mfaaa + 1. * oMdrho) * ( vz2 + vvz) * c1o2; + m2 = mfaac * c1o2 + mfaab * (vvz + c1o2) + (mfaaa + 1. * oMdrho) * (vz2 + vvz) * c1o2; mfaaa = m0; mfaab = m1; mfaac = m2; //////////////////////////////////////////////////////////////////////////////////// - m0 = mfabc * c1o2 + mfabb * (vvz - c1o2) + mfaba * ( vz2 - vvz) * c1o2; + m0 = mfabc * c1o2 + mfabb * (vvz - c1o2) + mfaba * (vz2 - vvz) * c1o2; m1 = -mfabc - 2. * mfabb * vvz + mfaba * (1. - vz2); - m2 = mfabc * c1o2 + mfabb * (vvz + c1o2) + mfaba * ( vz2 + vvz) * c1o2; + m2 = mfabc * c1o2 + mfabb * (vvz + c1o2) + mfaba * (vz2 + vvz) * c1o2; mfaba = m0; mfabb = m1; mfabc = m2; //////////////////////////////////////////////////////////////////////////////////// - m0 = mfacc * c1o2 + mfacb * (vvz - c1o2) + (mfaca + c1o3 * oMdrho) * ( vz2 - vvz) * c1o2; + m0 = mfacc * c1o2 + mfacb * (vvz - c1o2) + (mfaca + c1o3 * oMdrho) * (vz2 - vvz) * c1o2; m1 = -mfacc - 2. * mfacb * vvz + mfaca * (1. - vz2) - c1o3 * oMdrho * vz2; - m2 = mfacc * c1o2 + mfacb * (vvz + c1o2) + (mfaca + c1o3 * oMdrho) * ( vz2 + vvz) * c1o2; + m2 = mfacc * c1o2 + mfacb * (vvz + c1o2) + (mfaca + c1o3 * oMdrho) * (vz2 + vvz) * c1o2; mfaca = m0; mfacb = m1; mfacc = m2; //////////////////////////////////////////////////////////////////////////////////// //////////////////////////////////////////////////////////////////////////////////// - m0 = mfbac * c1o2 + mfbab * (vvz - c1o2) + mfbaa * ( vz2 - vvz) * c1o2; + m0 = mfbac * c1o2 + mfbab * (vvz - c1o2) + mfbaa * (vz2 - vvz) * c1o2; m1 = -mfbac - 2. * mfbab * vvz + mfbaa * (1. - vz2); - m2 = mfbac * c1o2 + mfbab * (vvz + c1o2) + mfbaa * ( vz2 + vvz) * c1o2; + m2 = mfbac * c1o2 + mfbab * (vvz + c1o2) + mfbaa * (vz2 + vvz) * c1o2; mfbaa = m0; mfbab = m1; mfbac = m2; /////////b////////////////////////////////////////////////////////////////////////// - m0 = mfbbc * c1o2 + mfbbb * (vvz - c1o2) + mfbba * ( vz2 - vvz) * c1o2; + m0 = mfbbc * c1o2 + mfbbb * (vvz - c1o2) + mfbba * (vz2 - vvz) * c1o2; m1 = -mfbbc - 2. * mfbbb * vvz + mfbba * (1. - vz2); - m2 = mfbbc * c1o2 + mfbbb * (vvz + c1o2) + mfbba * ( vz2 + vvz) * c1o2; + m2 = mfbbc * c1o2 + mfbbb * (vvz + c1o2) + mfbba * (vz2 + vvz) * c1o2; mfbba = m0; mfbbb = m1; mfbbc = m2; /////////b////////////////////////////////////////////////////////////////////////// - m0 = mfbcc * c1o2 + mfbcb * (vvz - c1o2) + mfbca * ( vz2 - vvz) * c1o2; + m0 = mfbcc * c1o2 + mfbcb * (vvz - c1o2) + mfbca * (vz2 - vvz) * c1o2; m1 = -mfbcc - 2. * mfbcb * vvz + mfbca * (1. - vz2); - m2 = mfbcc * c1o2 + mfbcb * (vvz + c1o2) + mfbca * ( vz2 + vvz) * c1o2; + m2 = mfbcc * c1o2 + mfbcb * (vvz + c1o2) + mfbca * (vz2 + vvz) * c1o2; mfbca = m0; mfbcb = m1; mfbcc = m2; //////////////////////////////////////////////////////////////////////////////////// //////////////////////////////////////////////////////////////////////////////////// - m0 = mfcac * c1o2 + mfcab * (vvz - c1o2) + (mfcaa + c1o3 * oMdrho) * ( vz2 - vvz) * c1o2; + m0 = mfcac * c1o2 + mfcab * (vvz - c1o2) + (mfcaa + c1o3 * oMdrho) * (vz2 - vvz) * c1o2; m1 = -mfcac - 2. * mfcab * vvz + mfcaa * (1. - vz2) - c1o3 * oMdrho * vz2; - m2 = mfcac * c1o2 + mfcab * (vvz + c1o2) + (mfcaa + c1o3 * oMdrho) * ( vz2 + vvz) * c1o2; + m2 = mfcac * c1o2 + mfcab * (vvz + c1o2) + (mfcaa + c1o3 * oMdrho) * (vz2 + vvz) * c1o2; mfcaa = m0; mfcab = m1; mfcac = m2; /////////c////////////////////////////////////////////////////////////////////////// - m0 = mfcbc * c1o2 + mfcbb * (vvz - c1o2) + mfcba * ( vz2 - vvz) * c1o2; + m0 = mfcbc * c1o2 + mfcbb * (vvz - c1o2) + mfcba * (vz2 - vvz) * c1o2; m1 = -mfcbc - 2. * mfcbb * vvz + mfcba * (1. - vz2); - m2 = mfcbc * c1o2 + mfcbb * (vvz + c1o2) + mfcba * ( vz2 + vvz) * c1o2; + m2 = mfcbc * c1o2 + mfcbb * (vvz + c1o2) + mfcba * (vz2 + vvz) * c1o2; mfcba = m0; mfcbb = m1; mfcbc = m2; /////////c////////////////////////////////////////////////////////////////////////// - m0 = mfccc * c1o2 + mfccb * (vvz - c1o2) + (mfcca + c1o9 * oMdrho) * ( vz2 - vvz) * c1o2; + m0 = mfccc * c1o2 + mfccb * (vvz - c1o2) + (mfcca + c1o9 * oMdrho) * (vz2 - vvz) * c1o2; m1 = -mfccc - 2. * mfccb * vvz + mfcca * (1. - vz2) - c1o9 * oMdrho * vz2; - m2 = mfccc * c1o2 + mfccb * (vvz + c1o2) + (mfcca + c1o9 * oMdrho) * ( vz2 + vvz) * c1o2; + m2 = mfccc * c1o2 + mfccb * (vvz + c1o2) + (mfcca + c1o9 * oMdrho) * (vz2 + vvz) * c1o2; mfcca = m0; mfccb = m1; mfccc = m2; @@ -715,67 +690,67 @@ void IncompressibleCumulantLBMKernel::collideAll() //mit 1/6, 2/3, 1/6, 0, 0, 0, 1/18, 2/9, 1/18 Konditionieren //////////////////////////////////////////////////////////////////////////////////// // Y - Dir - m0 = mfaca * c1o2 + mfaba * (vvy - c1o2) + (mfaaa + c1o6 * oMdrho) * ( vy2 - vvy) * c1o2; + m0 = mfaca * c1o2 + mfaba * (vvy - c1o2) + (mfaaa + c1o6 * oMdrho) * (vy2 - vvy) * c1o2; m1 = -mfaca - 2. * mfaba * vvy + mfaaa * (1. - vy2) - c1o6 * oMdrho * vy2; - m2 = mfaca * c1o2 + mfaba * (vvy + c1o2) + (mfaaa + c1o6 * oMdrho) * ( vy2 + vvy) * c1o2; + m2 = mfaca * c1o2 + mfaba * (vvy + c1o2) + (mfaaa + c1o6 * oMdrho) * (vy2 + vvy) * c1o2; mfaaa = m0; mfaba = m1; mfaca = m2; //////////////////////////////////////////////////////////////////////////////////// - m0 = mfacb * c1o2 + mfabb * (vvy - c1o2) + (mfaab + c2o3 * oMdrho) * ( vy2 - vvy) * c1o2; + m0 = mfacb * c1o2 + mfabb * (vvy - c1o2) + (mfaab + c2o3 * oMdrho) * (vy2 - vvy) * c1o2; m1 = -mfacb - 2. * mfabb * vvy + mfaab * (1. - vy2) - c2o3 * oMdrho * vy2; - m2 = mfacb * c1o2 + mfabb * (vvy + c1o2) + (mfaab + c2o3 * oMdrho) * ( vy2 + vvy) * c1o2; + m2 = mfacb * c1o2 + mfabb * (vvy + c1o2) + (mfaab + c2o3 * oMdrho) * (vy2 + vvy) * c1o2; mfaab = m0; mfabb = m1; mfacb = m2; //////////////////////////////////////////////////////////////////////////////////// - m0 = mfacc * c1o2 + mfabc * (vvy - c1o2) + (mfaac + c1o6 * oMdrho) * ( vy2 - vvy) * c1o2; + m0 = mfacc * c1o2 + mfabc * (vvy - c1o2) + (mfaac + c1o6 * oMdrho) * (vy2 - vvy) * c1o2; m1 = -mfacc - 2. * mfabc * vvy + mfaac * (1. - vy2) - c1o6 * oMdrho * vy2; - m2 = mfacc * c1o2 + mfabc * (vvy + c1o2) + (mfaac + c1o6 * oMdrho) * ( vy2 + vvy) * c1o2; + m2 = mfacc * c1o2 + mfabc * (vvy + c1o2) + (mfaac + c1o6 * oMdrho) * (vy2 + vvy) * c1o2; mfaac = m0; mfabc = m1; mfacc = m2; //////////////////////////////////////////////////////////////////////////////////// //////////////////////////////////////////////////////////////////////////////////// - m0 = mfbca * c1o2 + mfbba * (vvy - c1o2) + mfbaa * ( vy2 - vvy) * c1o2; + m0 = mfbca * c1o2 + mfbba * (vvy - c1o2) + mfbaa * (vy2 - vvy) * c1o2; m1 = -mfbca - 2. * mfbba * vvy + mfbaa * (1. - vy2); - m2 = mfbca * c1o2 + mfbba * (vvy + c1o2) + mfbaa * ( vy2 + vvy) * c1o2; + m2 = mfbca * c1o2 + mfbba * (vvy + c1o2) + mfbaa * (vy2 + vvy) * c1o2; mfbaa = m0; mfbba = m1; mfbca = m2; /////////b////////////////////////////////////////////////////////////////////////// - m0 = mfbcb * c1o2 + mfbbb * (vvy - c1o2) + mfbab * ( vy2 - vvy) * c1o2; + m0 = mfbcb * c1o2 + mfbbb * (vvy - c1o2) + mfbab * (vy2 - vvy) * c1o2; m1 = -mfbcb - 2. * mfbbb * vvy + mfbab * (1. - vy2); - m2 = mfbcb * c1o2 + mfbbb * (vvy + c1o2) + mfbab * ( vy2 + vvy) * c1o2; + m2 = mfbcb * c1o2 + mfbbb * (vvy + c1o2) + mfbab * (vy2 + vvy) * c1o2; mfbab = m0; mfbbb = m1; mfbcb = m2; /////////b////////////////////////////////////////////////////////////////////////// - m0 = mfbcc * c1o2 + mfbbc * (vvy - c1o2) + mfbac * ( vy2 - vvy) * c1o2; + m0 = mfbcc * c1o2 + mfbbc * (vvy - c1o2) + mfbac * (vy2 - vvy) * c1o2; m1 = -mfbcc - 2. * mfbbc * vvy + mfbac * (1. - vy2); - m2 = mfbcc * c1o2 + mfbbc * (vvy + c1o2) + mfbac * ( vy2 + vvy) * c1o2; + m2 = mfbcc * c1o2 + mfbbc * (vvy + c1o2) + mfbac * (vy2 + vvy) * c1o2; mfbac = m0; mfbbc = m1; mfbcc = m2; //////////////////////////////////////////////////////////////////////////////////// //////////////////////////////////////////////////////////////////////////////////// - m0 = mfcca * c1o2 + mfcba * (vvy - c1o2) + (mfcaa + c1o18 * oMdrho) * ( vy2 - vvy) * c1o2; + m0 = mfcca * c1o2 + mfcba * (vvy - c1o2) + (mfcaa + c1o18 * oMdrho) * (vy2 - vvy) * c1o2; m1 = -mfcca - 2. * mfcba * vvy + mfcaa * (1. - vy2) - c1o18 * oMdrho * vy2; - m2 = mfcca * c1o2 + mfcba * (vvy + c1o2) + (mfcaa + c1o18 * oMdrho) * ( vy2 + vvy) * c1o2; + m2 = mfcca * c1o2 + mfcba * (vvy + c1o2) + (mfcaa + c1o18 * oMdrho) * (vy2 + vvy) * c1o2; mfcaa = m0; mfcba = m1; mfcca = m2; /////////c////////////////////////////////////////////////////////////////////////// - m0 = mfccb * c1o2 + mfcbb * (vvy - c1o2) + (mfcab + c2o9 * oMdrho) * ( vy2 - vvy) * c1o2; + m0 = mfccb * c1o2 + mfcbb * (vvy - c1o2) + (mfcab + c2o9 * oMdrho) * (vy2 - vvy) * c1o2; m1 = -mfccb - 2. * mfcbb * vvy + mfcab * (1. - vy2) - c2o9 * oMdrho * vy2; - m2 = mfccb * c1o2 + mfcbb * (vvy + c1o2) + (mfcab + c2o9 * oMdrho) * ( vy2 + vvy) * c1o2; + m2 = mfccb * c1o2 + mfcbb * (vvy + c1o2) + (mfcab + c2o9 * oMdrho) * (vy2 + vvy) * c1o2; mfcab = m0; mfcbb = m1; mfccb = m2; /////////c////////////////////////////////////////////////////////////////////////// - m0 = mfccc * c1o2 + mfcbc * (vvy - c1o2) + (mfcac + c1o18 * oMdrho) * ( vy2 - vvy) * c1o2; + m0 = mfccc * c1o2 + mfcbc * (vvy - c1o2) + (mfcac + c1o18 * oMdrho) * (vy2 - vvy) * c1o2; m1 = -mfccc - 2. * mfcbc * vvy + mfcac * (1. - vy2) - c1o18 * oMdrho * vy2; - m2 = mfccc * c1o2 + mfcbc * (vvy + c1o2) + (mfcac + c1o18 * oMdrho) * ( vy2 + vvy) * c1o2; + m2 = mfccc * c1o2 + mfcbc * (vvy + c1o2) + (mfcac + c1o18 * oMdrho) * (vy2 + vvy) * c1o2; mfcac = m0; mfcbc = m1; mfccc = m2; @@ -784,67 +759,67 @@ void IncompressibleCumulantLBMKernel::collideAll() //mit 1/36, 1/9, 1/36, 1/9, 4/9, 1/9, 1/36, 1/9, 1/36 Konditionieren //////////////////////////////////////////////////////////////////////////////////// // X - Dir - m0 = mfcaa * c1o2 + mfbaa * (vvx - c1o2) + (mfaaa + c1o36 * oMdrho) * ( vx2 - vvx) * c1o2; + m0 = mfcaa * c1o2 + mfbaa * (vvx - c1o2) + (mfaaa + c1o36 * oMdrho) * (vx2 - vvx) * c1o2; m1 = -mfcaa - 2. * mfbaa * vvx + mfaaa * (1. - vx2) - c1o36 * oMdrho * vx2; - m2 = mfcaa * c1o2 + mfbaa * (vvx + c1o2) + (mfaaa + c1o36 * oMdrho) * ( vx2 + vvx) * c1o2; + m2 = mfcaa * c1o2 + mfbaa * (vvx + c1o2) + (mfaaa + c1o36 * oMdrho) * (vx2 + vvx) * c1o2; mfaaa = m0; mfbaa = m1; mfcaa = m2; //////////////////////////////////////////////////////////////////////////////////// - m0 = mfcba * c1o2 + mfbba * (vvx - c1o2) + (mfaba + c1o9 * oMdrho) * ( vx2 - vvx) * c1o2; + m0 = mfcba * c1o2 + mfbba * (vvx - c1o2) + (mfaba + c1o9 * oMdrho) * (vx2 - vvx) * c1o2; m1 = -mfcba - 2. * mfbba * vvx + mfaba * (1. - vx2) - c1o9 * oMdrho * vx2; - m2 = mfcba * c1o2 + mfbba * (vvx + c1o2) + (mfaba + c1o9 * oMdrho) * ( vx2 + vvx) * c1o2; + m2 = mfcba * c1o2 + mfbba * (vvx + c1o2) + (mfaba + c1o9 * oMdrho) * (vx2 + vvx) * c1o2; mfaba = m0; mfbba = m1; mfcba = m2; //////////////////////////////////////////////////////////////////////////////////// - m0 = mfcca * c1o2 + mfbca * (vvx - c1o2) + (mfaca + c1o36 * oMdrho) * ( vx2 - vvx) * c1o2; + m0 = mfcca * c1o2 + mfbca * (vvx - c1o2) + (mfaca + c1o36 * oMdrho) * (vx2 - vvx) * c1o2; m1 = -mfcca - 2. * mfbca * vvx + mfaca * (1. - vx2) - c1o36 * oMdrho * vx2; - m2 = mfcca * c1o2 + mfbca * (vvx + c1o2) + (mfaca + c1o36 * oMdrho) * ( vx2 + vvx) * c1o2; + m2 = mfcca * c1o2 + mfbca * (vvx + c1o2) + (mfaca + c1o36 * oMdrho) * (vx2 + vvx) * c1o2; mfaca = m0; mfbca = m1; mfcca = m2; //////////////////////////////////////////////////////////////////////////////////// //////////////////////////////////////////////////////////////////////////////////// - m0 = mfcab * c1o2 + mfbab * (vvx - c1o2) + (mfaab + c1o9 * oMdrho) * ( vx2 - vvx) * c1o2; + m0 = mfcab * c1o2 + mfbab * (vvx - c1o2) + (mfaab + c1o9 * oMdrho) * (vx2 - vvx) * c1o2; m1 = -mfcab - 2. * mfbab * vvx + mfaab * (1. - vx2) - c1o9 * oMdrho * vx2; - m2 = mfcab * c1o2 + mfbab * (vvx + c1o2) + (mfaab + c1o9 * oMdrho) * ( vx2 + vvx) * c1o2; + m2 = mfcab * c1o2 + mfbab * (vvx + c1o2) + (mfaab + c1o9 * oMdrho) * (vx2 + vvx) * c1o2; mfaab = m0; mfbab = m1; mfcab = m2; ///////////b//////////////////////////////////////////////////////////////////////// - m0 = mfcbb * c1o2 + mfbbb * (vvx - c1o2) + (mfabb + c4o9 * oMdrho) * ( vx2 - vvx) * c1o2; + m0 = mfcbb * c1o2 + mfbbb * (vvx - c1o2) + (mfabb + c4o9 * oMdrho) * (vx2 - vvx) * c1o2; m1 = -mfcbb - 2. * mfbbb * vvx + mfabb * (1. - vx2) - c4o9 * oMdrho * vx2; - m2 = mfcbb * c1o2 + mfbbb * (vvx + c1o2) + (mfabb + c4o9 * oMdrho) * ( vx2 + vvx) * c1o2; + m2 = mfcbb * c1o2 + mfbbb * (vvx + c1o2) + (mfabb + c4o9 * oMdrho) * (vx2 + vvx) * c1o2; mfabb = m0; mfbbb = m1; mfcbb = m2; ///////////b//////////////////////////////////////////////////////////////////////// - m0 = mfccb * c1o2 + mfbcb * (vvx - c1o2) + (mfacb + c1o9 * oMdrho) * ( vx2 - vvx) * c1o2; + m0 = mfccb * c1o2 + mfbcb * (vvx - c1o2) + (mfacb + c1o9 * oMdrho) * (vx2 - vvx) * c1o2; m1 = -mfccb - 2. * mfbcb * vvx + mfacb * (1. - vx2) - c1o9 * oMdrho * vx2; - m2 = mfccb * c1o2 + mfbcb * (vvx + c1o2) + (mfacb + c1o9 * oMdrho) * ( vx2 + vvx) * c1o2; + m2 = mfccb * c1o2 + mfbcb * (vvx + c1o2) + (mfacb + c1o9 * oMdrho) * (vx2 + vvx) * c1o2; mfacb = m0; mfbcb = m1; mfccb = m2; //////////////////////////////////////////////////////////////////////////////////// //////////////////////////////////////////////////////////////////////////////////// - m0 = mfcac * c1o2 + mfbac * (vvx - c1o2) + (mfaac + c1o36 * oMdrho) * ( vx2 - vvx) * c1o2; + m0 = mfcac * c1o2 + mfbac * (vvx - c1o2) + (mfaac + c1o36 * oMdrho) * (vx2 - vvx) * c1o2; m1 = -mfcac - 2. * mfbac * vvx + mfaac * (1. - vx2) - c1o36 * oMdrho * vx2; - m2 = mfcac * c1o2 + mfbac * (vvx + c1o2) + (mfaac + c1o36 * oMdrho) * ( vx2 + vvx) * c1o2; + m2 = mfcac * c1o2 + mfbac * (vvx + c1o2) + (mfaac + c1o36 * oMdrho) * (vx2 + vvx) * c1o2; mfaac = m0; mfbac = m1; mfcac = m2; ///////////c//////////////////////////////////////////////////////////////////////// - m0 = mfcbc * c1o2 + mfbbc * (vvx - c1o2) + (mfabc + c1o9 * oMdrho) * ( vx2 - vvx) * c1o2; + m0 = mfcbc * c1o2 + mfbbc * (vvx - c1o2) + (mfabc + c1o9 * oMdrho) * (vx2 - vvx) * c1o2; m1 = -mfcbc - 2. * mfbbc * vvx + mfabc * (1. - vx2) - c1o9 * oMdrho * vx2; - m2 = mfcbc * c1o2 + mfbbc * (vvx + c1o2) + (mfabc + c1o9 * oMdrho) * ( vx2 + vvx) * c1o2; + m2 = mfcbc * c1o2 + mfbbc * (vvx + c1o2) + (mfabc + c1o9 * oMdrho) * (vx2 + vvx) * c1o2; mfabc = m0; mfbbc = m1; mfcbc = m2; ///////////c//////////////////////////////////////////////////////////////////////// - m0 = mfccc * c1o2 + mfbcc * (vvx - c1o2) + (mfacc + c1o36 * oMdrho) * ( vx2 - vvx) * c1o2; + m0 = mfccc * c1o2 + mfbcc * (vvx - c1o2) + (mfacc + c1o36 * oMdrho) * (vx2 - vvx) * c1o2; m1 = -mfccc - 2. * mfbcc * vvx + mfacc * (1. - vx2) - c1o36 * oMdrho * vx2; - m2 = mfccc * c1o2 + mfbcc * (vvx + c1o2) + (mfacc + c1o36 * oMdrho) * ( vx2 + vvx) * c1o2; + m2 = mfccc * c1o2 + mfbcc * (vvx + c1o2) + (mfacc + c1o36 * oMdrho) * (vx2 + vvx) * c1o2; mfacc = m0; mfbcc = m1; mfccc = m2; @@ -855,16 +830,16 @@ void IncompressibleCumulantLBMKernel::collideAll() #ifdef PROOF_CORRECTNESS LBMReal rho_post = (mfaaa+mfaac+mfaca+mfcaa+mfacc+mfcac+mfccc+mfcca) +(mfaab+mfacb+mfcab+mfccb)+(mfaba+mfabc+mfcba+mfcbc)+(mfbaa+mfbac+mfbca+mfbcc) - +(mfabb+mfcbb)+(mfbab+mfbcb)+(mfbba+mfbbc)+mfbbb; + +(mfabb+mfcbb)+(mfbab+mfbcb)+(mfbba+mfbbc)+mfbbb; //LBMReal dif = fabs(rho - rho_post); LBMReal dif = rho - rho_post; #ifdef SINGLEPRECISION - if(dif > 10.0E-7 || dif < -10.0E-7) + if (dif > 10.0E-7 || dif < -10.0E-7) #else - if(dif > 10.0E-15 || dif < -10.0E-15) + if (dif > 10.0E-15 || dif < -10.0E-15) #endif { - UB_THROW(UbException(UB_EXARGS,"rho="+UbSystem::toString(rho)+", rho_post="+UbSystem::toString(rho_post) + UB_THROW(UbException(UB_EXARGS, "rho="+UbSystem::toString(rho)+", rho_post="+UbSystem::toString(rho_post) +" dif="+UbSystem::toString(dif) +" rho is not correct for node "+UbSystem::toString(x1)+","+UbSystem::toString(x2)+","+UbSystem::toString(x3))); //UBLOG(logERROR,"LBMKernelETD3Q27CCLB::collideAll(): rho is not correct for node "+UbSystem::toString(x1)+","+UbSystem::toString(x2)+","+UbSystem::toString(x3)); @@ -874,43 +849,42 @@ void IncompressibleCumulantLBMKernel::collideAll() ////////////////////////////////////////////////////////////////////////// //write distribution ////////////////////////////////////////////////////////////////////////// - (*this->localDistributions)(D3Q27System::ET_E,x1, x2, x3) = mfabb; - (*this->localDistributions)(D3Q27System::ET_N,x1, x2, x3) = mfbab; - (*this->localDistributions)(D3Q27System::ET_T,x1, x2, x3) = mfbba; - (*this->localDistributions)(D3Q27System::ET_NE,x1, x2, x3) = mfaab; - (*this->localDistributions)(D3Q27System::ET_NW,x1p,x2, x3) = mfcab; - (*this->localDistributions)(D3Q27System::ET_TE,x1, x2, x3) = mfaba; - (*this->localDistributions)(D3Q27System::ET_TW,x1p,x2, x3) = mfcba; - (*this->localDistributions)(D3Q27System::ET_TN,x1, x2, x3) = mfbaa; - (*this->localDistributions)(D3Q27System::ET_TS,x1, x2p,x3) = mfbca; - (*this->localDistributions)(D3Q27System::ET_TNE,x1, x2, x3) = mfaaa; - (*this->localDistributions)(D3Q27System::ET_TNW,x1p,x2, x3) = mfcaa; - (*this->localDistributions)(D3Q27System::ET_TSE,x1, x2p,x3) = mfaca; - (*this->localDistributions)(D3Q27System::ET_TSW,x1p,x2p,x3) = mfcca; - - (*this->nonLocalDistributions)(D3Q27System::ET_W,x1p,x2, x3 ) = mfcbb; - (*this->nonLocalDistributions)(D3Q27System::ET_S,x1, x2p,x3 ) = mfbcb; - (*this->nonLocalDistributions)(D3Q27System::ET_B,x1, x2, x3p ) = mfbbc; - (*this->nonLocalDistributions)(D3Q27System::ET_SW,x1p,x2p,x3 ) = mfccb; - (*this->nonLocalDistributions)(D3Q27System::ET_SE,x1, x2p,x3 ) = mfacb; - (*this->nonLocalDistributions)(D3Q27System::ET_BW,x1p,x2, x3p ) = mfcbc; - (*this->nonLocalDistributions)(D3Q27System::ET_BE,x1, x2, x3p ) = mfabc; - (*this->nonLocalDistributions)(D3Q27System::ET_BS,x1, x2p,x3p ) = mfbcc; - (*this->nonLocalDistributions)(D3Q27System::ET_BN,x1, x2, x3p ) = mfbac; - (*this->nonLocalDistributions)(D3Q27System::ET_BSW,x1p,x2p,x3p) = mfccc; - (*this->nonLocalDistributions)(D3Q27System::ET_BSE,x1, x2p,x3p) = mfacc; - (*this->nonLocalDistributions)(D3Q27System::ET_BNW,x1p,x2, x3p) = mfcac; - (*this->nonLocalDistributions)(D3Q27System::ET_BNE,x1, x2, x3p) = mfaac; - - (*this->zeroDistributions)(x1,x2,x3) = mfbbb; + (*this->localDistributions)(D3Q27System::ET_E, x1, x2, x3) = mfabb; + (*this->localDistributions)(D3Q27System::ET_N, x1, x2, x3) = mfbab; + (*this->localDistributions)(D3Q27System::ET_T, x1, x2, x3) = mfbba; + (*this->localDistributions)(D3Q27System::ET_NE, x1, x2, x3) = mfaab; + (*this->localDistributions)(D3Q27System::ET_NW, x1p, x2, x3) = mfcab; + (*this->localDistributions)(D3Q27System::ET_TE, x1, x2, x3) = mfaba; + (*this->localDistributions)(D3Q27System::ET_TW, x1p, x2, x3) = mfcba; + (*this->localDistributions)(D3Q27System::ET_TN, x1, x2, x3) = mfbaa; + (*this->localDistributions)(D3Q27System::ET_TS, x1, x2p, x3) = mfbca; + (*this->localDistributions)(D3Q27System::ET_TNE, x1, x2, x3) = mfaaa; + (*this->localDistributions)(D3Q27System::ET_TNW, x1p, x2, x3) = mfcaa; + (*this->localDistributions)(D3Q27System::ET_TSE, x1, x2p, x3) = mfaca; + (*this->localDistributions)(D3Q27System::ET_TSW, x1p, x2p, x3) = mfcca; + + (*this->nonLocalDistributions)(D3Q27System::ET_W, x1p, x2, x3) = mfcbb; + (*this->nonLocalDistributions)(D3Q27System::ET_S, x1, x2p, x3) = mfbcb; + (*this->nonLocalDistributions)(D3Q27System::ET_B, x1, x2, x3p) = mfbbc; + (*this->nonLocalDistributions)(D3Q27System::ET_SW, x1p, x2p, x3) = mfccb; + (*this->nonLocalDistributions)(D3Q27System::ET_SE, x1, x2p, x3) = mfacb; + (*this->nonLocalDistributions)(D3Q27System::ET_BW, x1p, x2, x3p) = mfcbc; + (*this->nonLocalDistributions)(D3Q27System::ET_BE, x1, x2, x3p) = mfabc; + (*this->nonLocalDistributions)(D3Q27System::ET_BS, x1, x2p, x3p) = mfbcc; + (*this->nonLocalDistributions)(D3Q27System::ET_BN, x1, x2, x3p) = mfbac; + (*this->nonLocalDistributions)(D3Q27System::ET_BSW, x1p, x2p, x3p) = mfccc; + (*this->nonLocalDistributions)(D3Q27System::ET_BSE, x1, x2p, x3p) = mfacc; + (*this->nonLocalDistributions)(D3Q27System::ET_BNW, x1p, x2, x3p) = mfcac; + (*this->nonLocalDistributions)(D3Q27System::ET_BNE, x1, x2, x3p) = mfaac; + + (*this->zeroDistributions)(x1, x2, x3) = mfbbb; ////////////////////////////////////////////////////////////////////////// } } } } - - } + //timer.stop(); } ////////////////////////////////////////////////////////////////////////// double IncompressibleCumulantLBMKernel::getCalculationTime() @@ -918,3 +892,9 @@ double IncompressibleCumulantLBMKernel::getCalculationTime() //return timer.getDuration(); return timer.getTotalTime(); } +////////////////////////////////////////////////////////////////////////// +void IncompressibleCumulantLBMKernel::setRelaxationParameter(Parameter p) +{ + parameter = p; +} + diff --git a/source/VirtualFluidsCore/LBM/IncompressibleCumulantLBMKernel.h b/source/VirtualFluidsCore/LBM/IncompressibleCumulantLBMKernel.h index 8199471268b5bd9e6ad1da6c59e507bcdf13c9be..bab1a42addf54bd5830d4fce15fb5844c4b65653 100644 --- a/source/VirtualFluidsCore/LBM/IncompressibleCumulantLBMKernel.h +++ b/source/VirtualFluidsCore/LBM/IncompressibleCumulantLBMKernel.h @@ -10,9 +10,6 @@ #include "basics/container/CbArray4D.h" #include "basics/container/CbArray3D.h" -class IncompressibleCumulantLBMKernel; -typedef SPtr<IncompressibleCumulantLBMKernel> LBMKernelETD3Q27CCLBPtr; - //! \brief Cascaded Cumulant LBM kernel. //! \details CFD solver that use Cascaded Cumulant Lattice Boltzmann method for D3Q27 model //! \author K. Kutscher, M. Geier @@ -23,20 +20,13 @@ public: enum Parameter{NORMAL, MAGIC}; public: IncompressibleCumulantLBMKernel(); - //! Constructor - //! \param nx1 number of nodes in x dimension - //! \param nx2 number of nodes in y dimension - //! \param nx3 number of nodes in z dimension - //! \param p set relaxation parameter: NORMAL is OxyyMxzz = 1.0 and MAGIC is OxyyMxzz = 2.0 +(-collFactor) - IncompressibleCumulantLBMKernel(int nx1, int nx2, int nx3, Parameter p); virtual ~IncompressibleCumulantLBMKernel(void); virtual void calculate(); virtual SPtr<LBMKernel> clone(); double getCalculationTime(); - + void setRelaxationParameter(Parameter p); protected: - virtual void collideAll(); - virtual void init(); + virtual void initDataSet(); LBMReal f[D3Q27System::ENDF+1]; UbTimer timer; diff --git a/source/VirtualFluidsCore/LBM/IncompressibleCumulantWithSpongeLayerLBMKernel.cpp b/source/VirtualFluidsCore/LBM/IncompressibleCumulantWithSpongeLayerLBMKernel.cpp index 401d65c7f07c6d16ad59950e90e31651d9dc08cc..df0859c6545d8ff9439e23e5b87cad8487415a34 100644 --- a/source/VirtualFluidsCore/LBM/IncompressibleCumulantWithSpongeLayerLBMKernel.cpp +++ b/source/VirtualFluidsCore/LBM/IncompressibleCumulantWithSpongeLayerLBMKernel.cpp @@ -11,12 +11,6 @@ IncompressibleCumulantWithSpongeLayerLBMKernel::IncompressibleCumulantWithSpongeLayerLBMKernel() { -} -////////////////////////////////////////////////////////////////////////// -IncompressibleCumulantWithSpongeLayerLBMKernel::IncompressibleCumulantWithSpongeLayerLBMKernel(int nx1, int nx2, int nx3, Parameter p) - : IncompressibleCumulantLBMKernel(nx1, nx2, nx3, p) -{ - } ////////////////////////////////////////////////////////////////////////// IncompressibleCumulantWithSpongeLayerLBMKernel::~IncompressibleCumulantWithSpongeLayerLBMKernel(void) @@ -24,10 +18,9 @@ IncompressibleCumulantWithSpongeLayerLBMKernel::~IncompressibleCumulantWithSpong } ////////////////////////////////////////////////////////////////////////// -void IncompressibleCumulantWithSpongeLayerLBMKernel::init() +void IncompressibleCumulantWithSpongeLayerLBMKernel::initDataSet() { - //SPtr<DistributionArray3D> d(new D3Q27EsoTwist3DSplittedVector(nx1+ghostLayerWitdh*2, nx2+ghostLayerWitdh*2, nx3+ghostLayerWitdh*2, -999.0)); - SPtr<DistributionArray3D> d(new D3Q27EsoTwist3DSplittedVector(nx1+2, nx2+2, nx3+2, -999.0)); + SPtr<DistributionArray3D> d(new D3Q27EsoTwist3DSplittedVector(nx[0]+2, nx[1]+2, nx[2]+2, -999.9)); dataSet->setFdistributions(d); } ////////////////////////////////////////////////////////////////////////// @@ -123,8 +116,9 @@ void IncompressibleCumulantWithSpongeLayerLBMKernel::initRelaxFactor(int vdir, d ////////////////////////////////////////////////////////////////////////// SPtr<LBMKernel> IncompressibleCumulantWithSpongeLayerLBMKernel::clone() { - SPtr<LBMKernel> kernel(new IncompressibleCumulantWithSpongeLayerLBMKernel(nx1, nx2, nx3, parameter)); - dynamicPointerCast<IncompressibleCumulantWithSpongeLayerLBMKernel>(kernel)->init(); + SPtr<LBMKernel> kernel(new IncompressibleCumulantWithSpongeLayerLBMKernel()); + kernel->setNX(nx); + dynamicPointerCast<IncompressibleCumulantWithSpongeLayerLBMKernel>(kernel)->initDataSet(); kernel->setCollisionFactor(this->collFactor); kernel->setBCProcessor(bcProcessor->clone(kernel)); kernel->setWithForcing(withForcing); @@ -150,15 +144,9 @@ SPtr<LBMKernel> IncompressibleCumulantWithSpongeLayerLBMKernel::clone() } ////////////////////////////////////////////////////////////////////////// void IncompressibleCumulantWithSpongeLayerLBMKernel::calculate() -{ - timer.resetAndStart(); - collideAll(); - timer.stop(); -} -////////////////////////////////////////////////////////////////////////// -void IncompressibleCumulantWithSpongeLayerLBMKernel::collideAll() { using namespace D3Q27System; + using namespace std; if(!withSpongeLayer) { diff --git a/source/VirtualFluidsCore/LBM/IncompressibleCumulantWithSpongeLayerLBMKernel.h b/source/VirtualFluidsCore/LBM/IncompressibleCumulantWithSpongeLayerLBMKernel.h index fe92889bac1a46a8e87e745afdea5b33a888e008..8918e4fc38741078aef47da583e34ddaee16ec2d 100644 --- a/source/VirtualFluidsCore/LBM/IncompressibleCumulantWithSpongeLayerLBMKernel.h +++ b/source/VirtualFluidsCore/LBM/IncompressibleCumulantWithSpongeLayerLBMKernel.h @@ -2,14 +2,6 @@ #define IncompressibleCumulantWithSpongeLayerLBMKernel_H #include "IncompressibleCumulantLBMKernel.h" -#include "BCProcessor.h" -#include "D3Q27System.h" -#include "basics/utilities/UbTiming.h" -#include "basics/container/CbArray4D.h" -#include "basics/container/CbArray3D.h" - -class IncompressibleCumulantWithSpongeLayerLBMKernel; -typedef SPtr<IncompressibleCumulantWithSpongeLayerLBMKernel> LBMKernelETD3Q27CCLBWithSpongeLayerPtr; //! \brief Cascaded Cumulant LBM kernel. //! \details CFD solver with sponge layer that use Cascaded Cumulant Lattice Boltzmann method for D3Q27 model <br> @@ -29,12 +21,6 @@ class IncompressibleCumulantWithSpongeLayerLBMKernel : public IncompressibleCum { public: IncompressibleCumulantWithSpongeLayerLBMKernel(); - //! Constructor - //! \param nx1 number of nodes in x dimension - //! \param nx2 number of nodes in y dimension - //! \param nx3 number of nodes in z dimension - //! \param p set relaxation parameter: NORMAL is OxyyMxzz = 1.0 and MAGIC is OxyyMxzz = 2.0 +(-collFactor) - IncompressibleCumulantWithSpongeLayerLBMKernel(int nx1, int nx2, int nx3, Parameter p); virtual ~IncompressibleCumulantWithSpongeLayerLBMKernel(void); SPtr<LBMKernel> clone(); void calculate(); @@ -45,14 +31,12 @@ public: //! \param vSP length of sponge layer void setRelaxFactorParam(int vdir, double vL1, double vdx, double vSP); protected: - void init(); + void initDataSet(); LBMReal OxyyMxzz; int direction; double L1; double dx; double SP; - - void collideAll(); }; #endif diff --git a/source/VirtualFluidsCore/LBM/InitDensityLBMKernel.cpp b/source/VirtualFluidsCore/LBM/InitDensityLBMKernel.cpp index 0edad93192153f120d97767c9ef3463608fdb918..3073d75f9d0f06a2dccd7918a611638ae7cc6a17 100644 --- a/source/VirtualFluidsCore/LBM/InitDensityLBMKernel.cpp +++ b/source/VirtualFluidsCore/LBM/InitDensityLBMKernel.cpp @@ -6,29 +6,26 @@ InitDensityLBMKernel::InitDensityLBMKernel() { + this->compressible = false; } InitDensityLBMKernel::~InitDensityLBMKernel() { } -InitDensityLBMKernel::InitDensityLBMKernel(int nx1, int nx2, int nx3) +void InitDensityLBMKernel::initDataSet() { - this->nx1 = nx1; - this->nx2 = nx2; - this->nx3 = nx3; - this->compressible = false; -} + SPtr<DistributionArray3D> d(new D3Q27EsoTwist3DSplittedVector(nx[0]+2, nx[1]+2, nx[2]+2, -999.9)); + dataSet->setFdistributions(d); + v.resize(3, nx[0]+2, nx[1]+2, nx[2]+2); -void InitDensityLBMKernel::calculate() -{ - collideAll(); } SPtr<LBMKernel> InitDensityLBMKernel::clone() { - SPtr<LBMKernel> kernel(new InitDensityLBMKernel(nx1, nx2, nx3)); - dynamicPointerCast<InitDensityLBMKernel>(kernel)->init(); + SPtr<LBMKernel> kernel(new InitDensityLBMKernel()); + kernel->setNX(nx); + dynamicPointerCast<InitDensityLBMKernel>(kernel)->initDataSet(); kernel->setCollisionFactor(this->collFactor); kernel->setBCProcessor(bcProcessor->clone(kernel)); kernel->setWithForcing(withForcing); @@ -852,16 +849,10 @@ double InitDensityLBMKernel::getCalculationTime() // //} -void InitDensityLBMKernel::init() -{ - SPtr<DistributionArray3D> d(new D3Q27EsoTwist3DSplittedVector(nx1+2, nx2+2, nx3+2, -999.0)); - dataSet->setFdistributions(d); - v.resize(3, nx1+2, nx2+2, nx3+2); -} -void InitDensityLBMKernel::collideAll() +void InitDensityLBMKernel::calculate() { using namespace D3Q27System; diff --git a/source/VirtualFluidsCore/LBM/InitDensityLBMKernel.h b/source/VirtualFluidsCore/LBM/InitDensityLBMKernel.h index 83efdeb441b3a4cfce9b36837b6722c7f4c59b60..6af546096ddbe0aef3fe3d606129ae55548d7fb0 100644 --- a/source/VirtualFluidsCore/LBM/InitDensityLBMKernel.h +++ b/source/VirtualFluidsCore/LBM/InitDensityLBMKernel.h @@ -12,14 +12,12 @@ class InitDensityLBMKernel : public LBMKernel public: InitDensityLBMKernel(); ~InitDensityLBMKernel(); - InitDensityLBMKernel(int nx1, int nx2, int nx3); void calculate(); SPtr<LBMKernel> clone(); void setVelocity(int x1, int x2, int x3, LBMReal vvx, LBMReal vvy, LBMReal vvz); double getCalculationTime(); protected: - void init(); - void collideAll(); + void initDataSet(); private: LBMReal f[D3Q27System::ENDF+1]; CbArray4D<LBMReal, IndexerX4X3X2X1>::CbArray4DPtr localDistributions; diff --git a/source/VirtualFluidsCore/LBM/LBMKernel.cpp b/source/VirtualFluidsCore/LBM/LBMKernel.cpp index 4cdf0f223e4142b7e7b7c3eb843436694106c83b..256443790170d8431699d55e79c7c455a04cf396 100644 --- a/source/VirtualFluidsCore/LBM/LBMKernel.cpp +++ b/source/VirtualFluidsCore/LBM/LBMKernel.cpp @@ -14,6 +14,9 @@ LBMKernel::LBMKernel() : ghostLayerWidth(1), this->setForcingX2(0.0); this->setForcingX3(0.0); dataSet = SPtr<DataSet3D>(new DataSet3D()); + this->nx[0] = 0; + this->nx[1] = 0; + this->nx[2] = 0; } ////////////////////////////////////////////////////////////////////////// LBMKernel::~LBMKernel() @@ -208,6 +211,16 @@ void LBMKernel::swapDistributions() dataSet->getFdistributions()->swap(); } ////////////////////////////////////////////////////////////////////////// +void LBMKernel::setNX(std::array<int, 3> nx) +{ + this->nx = nx; +} +////////////////////////////////////////////////////////////////////////// +std::array<int, 3> LBMKernel::getNX() +{ + return nx; +} +////////////////////////////////////////////////////////////////////////// bool LBMKernel::isInsideOfDomain(const int& x1, const int& x2, const int& x3) const { const SPtr<BCArray3D> bcArray = this->bcProcessor->getBCArray(); diff --git a/source/VirtualFluidsCore/LBM/LBMKernel.h b/source/VirtualFluidsCore/LBM/LBMKernel.h index 2152d7a23733013605c60f1764cfe317d264bb05..4854a5181cd9b4c592e621e323f74a7c43704f7d 100644 --- a/source/VirtualFluidsCore/LBM/LBMKernel.h +++ b/source/VirtualFluidsCore/LBM/LBMKernel.h @@ -2,12 +2,10 @@ #define LBMKERNEL_H #include <PointerDefinitions.h> - #include "LBMSystem.h" - -#include <MuParser/include/muParser.h> - #include "ILBMKernel.h" +#include <array> +#include <MuParser/include/muParser.h> class BCProcessor; class DataSet3D; @@ -72,9 +70,11 @@ public: bool isInsideOfDomain(const int &x1, const int &x2, const int &x3) const; - void swapDistributions(); + void setNX(std::array<int, 3> nx); + std::array<int, 3> getNX(); + protected: SPtr<DataSet3D> dataSet; SPtr<BCProcessor> bcProcessor; @@ -96,7 +96,7 @@ protected: WPtr<Block3D> block; - int nx1, nx2, nx3; + std::array<int, 3> nx; private: void checkFunction(mu::Parser fct); diff --git a/source/VirtualFluidsCore/LBM/VoidLBMKernel.cpp b/source/VirtualFluidsCore/LBM/VoidLBMKernel.cpp index 250664a6dfd0c218d3b9ed2e33c5b760782c7cec..c9fb04cea686c9861897c4f6f140b054704e889e 100644 --- a/source/VirtualFluidsCore/LBM/VoidLBMKernel.cpp +++ b/source/VirtualFluidsCore/LBM/VoidLBMKernel.cpp @@ -7,21 +7,23 @@ VoidLBMKernel::VoidLBMKernel() { } - -VoidLBMKernel::VoidLBMKernel(int nx1, int nx2, int nx3) : nx1(nx1), nx2(nx2), nx3(nx3) -{ - SPtr<DistributionArray3D> d(new VoidData3D(nx1+2, nx2+2, nx3+2, -999.0)); - dataSet->setFdistributions(d); -} - +////////////////////////////////////////////////////////////////////////// VoidLBMKernel::~VoidLBMKernel() { } - +////////////////////////////////////////////////////////////////////////// +void VoidLBMKernel::initDataSet() +{ + SPtr<DistributionArray3D> d(new VoidData3D(nx[0]+2, nx[1]+2, nx[2]+2, -999.9)); + dataSet->setFdistributions(d); +} +////////////////////////////////////////////////////////////////////////// SPtr<LBMKernel> VoidLBMKernel::clone() { - SPtr<LBMKernel> kernel(new VoidLBMKernel(nx1, nx2, nx3)); + SPtr<LBMKernel> kernel(new VoidLBMKernel()); + kernel->setNX(nx); + dynamicPointerCast<VoidLBMKernel>(kernel)->initDataSet(); kernel->setCollisionFactor(this->collFactor); kernel->setBCProcessor(bcProcessor->clone(kernel)); kernel->setWithForcing(withForcing); @@ -32,17 +34,12 @@ SPtr<LBMKernel> VoidLBMKernel::clone() kernel->setDeltaT(deltaT); return kernel; } - +////////////////////////////////////////////////////////////////////////// void VoidLBMKernel::calculate() { } - -void VoidLBMKernel::swapDistributions() -{ - -} - +////////////////////////////////////////////////////////////////////////// double VoidLBMKernel::getCalculationTime() { return 0.0; diff --git a/source/VirtualFluidsCore/LBM/VoidLBMKernel.h b/source/VirtualFluidsCore/LBM/VoidLBMKernel.h index 9da59a7fdaa77d1a1b8a4c8a4a8c9769637d22ef..b67384034df696b7134077105adfddbe1f183e54 100644 --- a/source/VirtualFluidsCore/LBM/VoidLBMKernel.h +++ b/source/VirtualFluidsCore/LBM/VoidLBMKernel.h @@ -7,14 +7,12 @@ class VoidLBMKernel : public LBMKernel { public: VoidLBMKernel(); - VoidLBMKernel(int nx1, int nx2, int nx3); ~VoidLBMKernel(); SPtr<LBMKernel> clone(); void calculate(); - void swapDistributions(); double getCalculationTime(); + void initDataSet(); protected: -private: - int nx1, nx2, nx3; + }; #endif // VoidLBMKernel_h__ diff --git a/source/VirtualFluidsCore/Visitors/InitDistributionsBlockVisitor.cpp b/source/VirtualFluidsCore/Visitors/InitDistributionsBlockVisitor.cpp index c1cd589dfbb88363cd8ee2fdce62189f82b211fb..2b01489faa32af6f32c8d819370969068accc7d0 100644 --- a/source/VirtualFluidsCore/Visitors/InitDistributionsBlockVisitor.cpp +++ b/source/VirtualFluidsCore/Visitors/InitDistributionsBlockVisitor.cpp @@ -17,25 +17,6 @@ InitDistributionsBlockVisitor::InitDistributionsBlockVisitor() this->setRho(0.0); } ////////////////////////////////////////////////////////////////////////// -//D3Q27ETInitDistributionsBlockVisitor::D3Q27ETInitDistributionsBlockVisitor(LBMReal rho, LBMReal vx1, LBMReal vx2, LBMReal vx3) -// : Block3DVisitor(0, Grid3DSystem::MAXLEVEL) -//{ -// this->setVx1(vx1); -// this->setVx2(vx2); -// this->setVx3(vx3); -// this->setRho(rho); -//} -////////////////////////////////////////////////////////////////////////// -InitDistributionsBlockVisitor::InitDistributionsBlockVisitor( LBMReal nu, LBMReal rho, LBMReal vx1, LBMReal vx2, LBMReal vx3) - : Block3DVisitor(0, Grid3DSystem::MAXLEVEL) -{ - this->setVx1(vx1); - this->setVx2(vx2); - this->setVx3(vx3); - this->setRho(rho); - this->setNu(nu); -} -////////////////////////////////////////////////////////////////////////// void InitDistributionsBlockVisitor::setVx1( const mu::Parser& parser) { this->checkFunction(parser); @@ -117,8 +98,6 @@ void InitDistributionsBlockVisitor::visit(const SPtr<Grid3D> grid, SPtr<Block3D> //UbTupleDouble3 blockLengths = grid->getBlockLengths(block); //UbTupleDouble3 nodeOffset = grid->getNodeOffset(block); double dx = grid->getDeltaX(block); - LBMReal o = LBMSystem::calcCollisionFactor(nu, block->getLevel()); - //define vars for functions mu::value_type x1,x2,x3; @@ -150,7 +129,9 @@ void InitDistributionsBlockVisitor::visit(const SPtr<Grid3D> grid, SPtr<Block3D> //UbTupleDouble3 org = grid->getBlockWorldCoordinates(block); SPtr<BCArray3D> bcArray = kernel->getBCProcessor()->getBCArray(); - SPtr<EsoTwist3D> distributions = dynamicPointerCast<EsoTwist3D>(kernel->getDataSet()->getFdistributions()); + SPtr<DistributionArray3D> distributions = kernel->getDataSet()->getFdistributions(); + + LBMReal o = kernel->getCollisionFactor(); LBMReal f[D3Q27System::ENDF+1]; @@ -309,8 +290,6 @@ void InitDistributionsBlockVisitor::checkFunction(mu::Parser fct) } } ////////////////////////////////////////////////////////////////////////// -void InitDistributionsBlockVisitor::setNu( LBMReal nu ) -{ - this->nu = nu; -} + + diff --git a/source/VirtualFluidsCore/Visitors/InitDistributionsBlockVisitor.h b/source/VirtualFluidsCore/Visitors/InitDistributionsBlockVisitor.h index 6a4bd3505fcb91d3b30c4af709d5df1175944e0b..a192c562e1b044ebe4ac4a0ccdb24d244e5158c5 100644 --- a/source/VirtualFluidsCore/Visitors/InitDistributionsBlockVisitor.h +++ b/source/VirtualFluidsCore/Visitors/InitDistributionsBlockVisitor.h @@ -42,14 +42,6 @@ public: public: InitDistributionsBlockVisitor(); - //D3Q27ETInitDistributionsBlockVisitor(LBMReal rho, LBMReal vx1=0.0, LBMReal vx2=0.0, LBMReal vx3=0.0); - //! Constructor - //! \param nu - viscosity - //! \param rho - density - //! \param vx1 - velocity in x - //! \param vx2 - velocity in y - //! \param vx3 - velocity in z - InitDistributionsBlockVisitor( LBMReal nu, LBMReal rho, LBMReal vx1=0.0, LBMReal vx2=0.0, LBMReal vx3=0.0); ////////////////////////////////////////////////////////////////////////// //automatic vars are: x1,x2, x3 //ussage example: setVx1("x1*0.01+x2*0.003") @@ -68,7 +60,6 @@ public: void setVx2( LBMReal vx2 ); void setVx3( LBMReal vx3 ); void setRho( LBMReal rho ); - void setNu( LBMReal nu ); void visit(SPtr<Grid3D> grid, SPtr<Block3D> block) override; @@ -80,7 +71,6 @@ private: mu::Parser muVx2; mu::Parser muVx3; mu::Parser muRho; - LBMReal nu; }; #endif //D3Q27INITDISTRIBUTIONSPATCHVISITOR_H diff --git a/source/VirtualFluidsCore/Visitors/SpongeLayerBlockVisitor.cpp b/source/VirtualFluidsCore/Visitors/SpongeLayerBlockVisitor.cpp index fe444cafc52c87bfaa9c0004fb6a44c0ed39b6e6..cd2cd0e2bd6a099d4dd853af002605964579f3ff 100644 --- a/source/VirtualFluidsCore/Visitors/SpongeLayerBlockVisitor.cpp +++ b/source/VirtualFluidsCore/Visitors/SpongeLayerBlockVisitor.cpp @@ -6,6 +6,8 @@ #include "Block3D.h" #include "ILBMKernel.h" +#include "CompressibleCumulant4thOrderViscosityLBMKernel.h" + #include <numerics/geometry3d/GbCuboid3D.h> using namespace std; @@ -39,8 +41,31 @@ void SpongeLayerBlockVisitor::visit(SPtr<Grid3D> grid, SPtr<Block3D> block) if (boundingBox->isCellInsideGbObject3D(minX1, minX2, minX3, maxX1, maxX2, maxX3)) { - SPtr<ILBMKernel> kernel = block->getKernel(); - double oldCollFactor = kernel->getCollisionFactor(); + LBMReal collFactor = LBMSystem::calcCollisionFactor(viscosity, block->getLevel()); + kernel->setCollisionFactor(collFactor); + kernel->setIndex(block->getX1(), block->getX2(), block->getX3()); + kernel->setDeltaT(LBMSystem::getDeltaT(block->getLevel())); + kernel->setBlock(block); + SPtr<LBMKernel> newKernel = kernel->clone(); + + SPtr<DataSet3D> dataSet = block->getKernel()->getDataSet(); + if (!dataSet) + { + UB_THROW(UbException(UB_EXARGS, "It is not possible to change a DataSet in kernel! Old DataSet is not exist!")); + } + + newKernel->setDataSet(dataSet); + + SPtr<BCProcessor> bcProc = block->getKernel()->getBCProcessor(); + if (!bcProc) + { + UB_THROW(UbException(UB_EXARGS, "It is not possible to change a BCProcessor in kernel! Old BCProcessor is not exist!")); + } + newKernel->setBCProcessor(bcProc); + block->setKernel(newKernel); + + //SPtr<ILBMKernel> kernel = block->getKernel(); + double oldCollFactor = newKernel->getCollisionFactor(); int ibX1 = block->getX1(); @@ -50,11 +75,21 @@ void SpongeLayerBlockVisitor::visit(SPtr<Grid3D> grid, SPtr<Block3D> block) int ibMax = val<1>(ixMax)-val<1>(ixMin)+1; double index = ibX1-val<1>(ixMin)+1; - double newCollFactor = oldCollFactor - (oldCollFactor-1.5)/ibMax*index; + double newCollFactor = oldCollFactor - (oldCollFactor-1.0)/ibMax*index; - kernel->setCollisionFactor(newCollFactor); + newKernel->setCollisionFactor(newCollFactor); + block->setKernel(newKernel); } } } - +////////////////////////////////////////////////////////////////////////// +void SpongeLayerBlockVisitor::setKernel(SPtr<LBMKernel> k) +{ + kernel = k; +} +////////////////////////////////////////////////////////////////////////// +void SpongeLayerBlockVisitor::setViscosity(LBMReal v) +{ + viscosity = v; +} diff --git a/source/VirtualFluidsCore/Visitors/SpongeLayerBlockVisitor.h b/source/VirtualFluidsCore/Visitors/SpongeLayerBlockVisitor.h index e607e3a603b330ef7d6e2382f88a5ae972ae5c6c..c53a6a53ef78adb6cfeb7149d7c6a73d821eece4 100644 --- a/source/VirtualFluidsCore/Visitors/SpongeLayerBlockVisitor.h +++ b/source/VirtualFluidsCore/Visitors/SpongeLayerBlockVisitor.h @@ -7,6 +7,7 @@ class Grid3D; class Block3D; class GbCuboid3D; +class LBMKernel; //! \brief Set sponge layer for all blocks inside boundingBox //! \details This visitor is useful if you need to set or reset sponge layer in kernels (e.g. after restart because sponge layer is not serializable). @@ -19,9 +20,14 @@ public: void visit(SPtr<Grid3D> grid, SPtr<Block3D> block) override; + void setKernel(SPtr<LBMKernel> k); + void setViscosity(LBMReal v); + private: SPtr<GbCuboid3D> boundingBox; LBMReal collFactor; + SPtr<LBMKernel> kernel; + LBMReal viscosity; }; #endif // SetSpongeLayerBlockVisitor_h__