diff --git a/source/Applications.cmake b/source/Applications.cmake
index 41187a45eaa2c7f765b2d26f7e436117538d7182..7d3766f694ae5baf6c3a13dd48c98bb7f128a500 100644
--- a/source/Applications.cmake
+++ b/source/Applications.cmake
@@ -44,4 +44,5 @@ add_subdirectory(Applications/pChannel)
#add_subdirectory(Applications/pDisk)
add_subdirectory(Applications/BoxBenchmark)
add_subdirectory(Applications/DHIT)
-add_subdirectory(Applications/FNG)
\ No newline at end of file
+add_subdirectory(Applications/FNG)
+add_subdirectory(Applications/aperm)
\ No newline at end of file
diff --git a/source/Applications/DHIT/dhit.cfg b/source/Applications/DHIT/dhit.cfg
index 768c4a3df70988afc1b9589ba5358af3474dbd08..1585f232dd94f9312dff9127c1ae1d452c95dc6a 100644
--- a/source/Applications/DHIT/dhit.cfg
+++ b/source/Applications/DHIT/dhit.cfg
@@ -1,18 +1,20 @@
-pathname = d:/temp/DHIT
+pathname = d:/temp/DHIT_128_2
numOfThreads = 4
availMem = 11e9
#Grid
-length = 128 128 128
+length = 128 128 128
blocknx = 32 32 32
+initTime = 10000
+
outTime = 100
-endTime = 100
+endTime = 1000000
logToFile = false
#Simulation
-initFile = d:/Projects/DHIT/F_128.txt
+initFile = d:/Projects/DHIT/Velocities_128_2.txt
nuLB = 1.2395e-2
uRMS = 0.0234
lambda = 0.1
diff --git a/source/Applications/DHIT/dhit.cpp b/source/Applications/DHIT/dhit.cpp
index 622d0a8c583b7b42752792b141fd28a630e0f20d..850c34019ef19d9d1aacaaf9948078c655484a7b 100644
--- a/source/Applications/DHIT/dhit.cpp
+++ b/source/Applications/DHIT/dhit.cpp
@@ -25,6 +25,7 @@ void run(string configname)
double nuLB = config.getDouble("nuLB");
double uRMS = config.getDouble("uRMS");
double lambda = config.getDouble("lambda");
+ double initTime = config.getDouble("initTime");
CommunicatorPtr comm = MPICommunicator::getInstance();
int myid = comm->getProcessID();
@@ -60,9 +61,9 @@ void run(string configname)
double g_minX2 = 0.0;
double g_minX3 = 0.0;
- double g_maxX1 = length[0];
- double g_maxX2 = length[1];
- double g_maxX3 = length[2];
+ double g_maxX1 = length[0];//-1.0;
+ double g_maxX2 = length[1];//-1.0;
+ double g_maxX3 = length[2];//-1.0;
//geometry
GbObject3DPtr box(new GbCuboid3D(g_minX1, g_minX2, g_minX3, g_maxX1, g_maxX2, g_maxX3));
@@ -149,7 +150,8 @@ void run(string configname)
LBMKernel3DPtr kernel;
- kernel = LBMKernel3DPtr(new LBMKernelETD3Q27CCLB(blocknx[0], blocknx[1], blocknx[2], LBMKernelETD3Q27CCLB::NORMAL));
+ //kernel = LBMKernel3DPtr(new LBMKernelETD3Q27CCLB(blocknx[0], blocknx[1], blocknx[2], LBMKernelETD3Q27CCLB::NORMAL));
+ kernel = LBMKernel3DPtr(new InitDensityLBMKernel(blocknx[0], blocknx[1], blocknx[2]));
BCProcessorPtr bcProc(new D3Q27ETBCProcessor());
kernel->setBCProcessor(bcProc);
@@ -181,16 +183,31 @@ void run(string configname)
UBLOG(logINFO, "PID = " << myid << " Physical Memory currently used by current process: " << Utilities::getPhysMemUsedByMe());
}
- UbSchedulerPtr visSch(new UbScheduler(outTime));
- MacroscopicQuantitiesCoProcessor pp(grid, visSch, pathname, WbWriterVtkXmlBinary::getInstance(), conv);
+ UbSchedulerPtr outputSch(new UbScheduler(outTime));
+ MacroscopicQuantitiesCoProcessor pp(grid, outputSch, pathname, WbWriterVtkXmlBinary::getInstance(), conv);
UbSchedulerPtr nupsSch(new UbScheduler(10, 30, 100));
NUPSCounterCoProcessor npr(grid, nupsSch, numOfThreads, comm);
+ CalculationManagerPtr initialisation(new CalculationManager(grid, numOfThreads, initTime, outputSch));
+ if (myid == 0) UBLOG(logINFO, "Initialisation-start");
+ initialisation->calculate();
+ if (myid == 0) UBLOG(logINFO, "Initialisation-end");
+
+
+ kernel = LBMKernel3DPtr(new LBMKernelETD3Q27CCLB(blocknx[0], blocknx[1], blocknx[2], LBMKernelETD3Q27CCLB::NORMAL));
+ kernel->setBCProcessor(bcProc);
+ SetKernelBlockVisitor kernelVisitor2(kernel, nuLB, availMem, needMem, SetKernelBlockVisitor::Change);
+ grid->accept(kernelVisitor2);
+
+ UbSchedulerPtr visSch(new UbScheduler(outTime));
+
+ if (myid==0) UBLOG(logINFO, "Simulation-start");
+ grid->setTimeStep(initTime+1.0);
CalculationManagerPtr calculation(new CalculationManager(grid, numOfThreads, endTime, visSch));
- if (myid == 0) UBLOG(logINFO, "Simulation-start");
calculation->calculate();
- if (myid == 0) UBLOG(logINFO, "Simulation-end");
+ if (myid==0) UBLOG(logINFO, "Simulation-end");
+
}
catch (std::exception& e)
{
diff --git a/source/Applications/FNG/fng.cpp b/source/Applications/FNG/fng.cpp
index 0f683374477e8e85f2d8e99521a42b8a87b57c2c..9e3448ed12aa03553507ae10cbf498f2e1fa76ad 100644
--- a/source/Applications/FNG/fng.cpp
+++ b/source/Applications/FNG/fng.cpp
@@ -136,17 +136,17 @@ void run(string configname)
if (myid == 0)
{
UBLOG(logINFO, "Parameters:");
- UBLOG(logINFO, "* Re =" << Re);
- UBLOG(logINFO, "* Ma =" << Ma);
- UBLOG(logINFO, "* uReal =" << uReal);
- UBLOG(logINFO, "* nuReal =" << nueReal);
- UBLOG(logINFO, "* nuLB =" << nuLB);
- UBLOG(logINFO, "* velocity =" << uLB);
- UBLOG(logINFO, "* dx_base =" << deltaXcoarse/1000.0 << "m");
- UBLOG(logINFO, "* dx_refine =" << deltaXfine/1000.0 << "m");
-
- UBLOG(logINFO, "number of levels = " << refineLevel + 1);
- UBLOG(logINFO, "numOfThreads = " << numOfThreads);
+ UBLOG(logINFO, "* Re = "<<Re);
+ UBLOG(logINFO, "* Ma = "<<Ma);
+ UBLOG(logINFO, "* velocity (uReal) = "<<uReal<<" m/s");
+ UBLOG(logINFO, "* viscosity (nuReal) = "<<nueReal<<" m^2/s");
+ UBLOG(logINFO, "* velocity LB (uLB) = "<<uLB);
+ UBLOG(logINFO, "* viscosity LB (nuLB) = "<<nuLB);
+ UBLOG(logINFO, "* dx_base = "<<deltaXcoarse/1000.0<<" m");
+ UBLOG(logINFO, "* dx_refine = "<<deltaXfine/1000.0<<" m");
+ UBLOG(logINFO, "* number of levels = " << refineLevel + 1);
+ UBLOG(logINFO, "* number of threads = " << numOfThreads);
+ UBLOG(logINFO, "* number of processes = " << comm->getNumberOfProcesses());
UBLOG(logINFO, "Preprozess - start");
}
@@ -297,16 +297,16 @@ void run(string configname)
if (porousTralingEdge)
{
- //grid->setRank(0);
- //boost::dynamic_pointer_cast<D3Q27TriFaceMeshInteractor>(fngIntrTrailingEdge)->refineBlockGridToLevel(refineLevel, -1.0, 10.0);
- //grid->setRank(rank);
+ grid->setRank(0);
+ boost::dynamic_pointer_cast<D3Q27TriFaceMeshInteractor>(fngIntrTrailingEdge)->refineBlockGridToLevel(refineLevel, -2.0, refineDistance);
+ grid->setRank(rank);
- GbObject3DPtr trailingEdgeCube(new GbCuboid3D(fngMeshTrailingEdge->getX1Minimum()-blockLength, fngMeshTrailingEdge->getX2Minimum(), fngMeshTrailingEdge->getX3Minimum()-blockLength/2.0,
- fngMeshTrailingEdge->getX1Maximum()+blockLength, fngMeshTrailingEdge->getX2Maximum(), fngMeshTrailingEdge->getX3Maximum()+blockLength/2.0));
- if (myid == 0) GbSystem3D::writeGeoObject(trailingEdgeCube.get(), pathOut + "/geo/trailingEdgeCube", WbWriterVtkXmlASCII::getInstance());
+ //GbObject3DPtr trailingEdgeCube(new GbCuboid3D(fngMeshTrailingEdge->getX1Minimum()-blockLength, fngMeshTrailingEdge->getX2Minimum(), fngMeshTrailingEdge->getX3Minimum()-blockLength/2.0,
+ // fngMeshTrailingEdge->getX1Maximum()+blockLength, fngMeshTrailingEdge->getX2Maximum(), fngMeshTrailingEdge->getX3Maximum()+blockLength/2.0));
+ //if (myid == 0) GbSystem3D::writeGeoObject(trailingEdgeCube.get(), pathOut + "/geo/trailingEdgeCube", WbWriterVtkXmlASCII::getInstance());
- RefineCrossAndInsideGbObjectBlockVisitor refVisitor(trailingEdgeCube, refineLevel);
- grid->accept(refVisitor);
+ //RefineCrossAndInsideGbObjectBlockVisitor refVisitor(trailingEdgeCube, refineLevel);
+ //grid->accept(refVisitor);
}
RatioBlockVisitor ratioVisitor(refineLevel);
@@ -397,7 +397,7 @@ void run(string configname)
if (porousTralingEdge)
{
intHelper.addInteractor(fngIntrBodyPart);
- intHelper.addInteractor(fngIntrTrailingEdge);
+ //intHelper.addInteractor(fngIntrTrailingEdge);
}
else
{
diff --git a/source/Applications/FNG/fng15Bombadil.cfg b/source/Applications/FNG/fng15Bombadil.cfg
index 393ca8a897cab1fc5c8dbd3f5372c746ba811805..a27fbb5cf7eaddab5ea9becd3f8f18e4fe6428a2 100644
--- a/source/Applications/FNG/fng15Bombadil.cfg
+++ b/source/Applications/FNG/fng15Bombadil.cfg
@@ -1,4 +1,4 @@
-pathOut = d:/temp/fng2
+pathOut = d:/temp/fngPorous
pathGeo = d:/Projects/SFB880/FNG/A1_Forschungsdaten_Profilgeometrie_STL_CATIA_Rossian
#fngFileWhole = f16-ascii.stl
fngFileWhole = grundgeometrie-direkter-export.stl
@@ -11,7 +11,7 @@ availMem = 20e9
refineLevel = 8
#blockNx = 8 4 8
blockNx = 21 6 13
-uLB = 0.01
+uLB = 0.1
#x1min x1max x2min x2max x3min x3max [m]
boundingBox = -0.90 1.20 0.035 0.065 -0.65 0.65
@@ -34,6 +34,6 @@ endTime = 2000
logToFile = false
-porousTralingEdge = false
+porousTralingEdge = true
thinWall = false
diff --git a/source/Applications/aperm/CMakeLists.txt b/source/Applications/aperm/CMakeLists.txt
new file mode 100644
index 0000000000000000000000000000000000000000..673364ed8283715b7644df46b6fce1bb967ec44b
--- /dev/null
+++ b/source/Applications/aperm/CMakeLists.txt
@@ -0,0 +1,25 @@
+CMAKE_MINIMUM_REQUIRED(VERSION 2.8)
+
+########################################################
+## C++ PROJECT ###
+########################################################
+PROJECT(aperm)
+
+INCLUDE(${SOURCE_ROOT}/IncludsList.cmake)
+
+#################################################################
+### LOCAL FILES ###
+#################################################################
+FILE(GLOB SPECIFIC_FILES ${CMAKE_CURRENT_SOURCE_DIR}/*.h
+ ${CMAKE_CURRENT_SOURCE_DIR}/*.cpp
+ ${CMAKE_CURRENT_SOURCE_DIR}/*.hpp )
+
+SET(ALL_SOURCES ${ALL_SOURCES} ${SPECIFIC_FILES})
+SOURCE_GROUP(src FILES ${SPECIFIC_FILES})
+
+SET(CAB_ADDITIONAL_LINK_LIBRARIES VirtualFluids)
+
+#################################################################
+### CREATE PROJECT ###
+#################################################################
+CREATE_CAB_PROJECT(aperm BINARY)
diff --git a/source/Applications/aperm/aperm.cpp b/source/Applications/aperm/aperm.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..beae1f1e04387b05c1d513ebfbf559a1cff9d525
--- /dev/null
+++ b/source/Applications/aperm/aperm.cpp
@@ -0,0 +1,470 @@
+#include <iostream>
+#include <string>
+#include <VirtualFluids.h>
+
+using namespace std;
+
+void changeDP()
+{
+}
+//////////////////////////////////////////////////////////////////////////
+void run(string configname)
+{
+ try
+ {
+ ConfigurationFile config;
+ config.load(configname);
+
+ string pathname = config.getString("pathname");
+ string pathGeo = config.getString("pathGeo");
+ int numOfThreads = config.getInt("numOfThreads");
+ string sampleFilename = config.getString("sampleFilename");
+ int pmNX1 = config.getInt("pmNX1");
+ int pmNX2 = config.getInt("pmNX2");
+ int pmNX3 = config.getInt("pmNX3");
+ double lthreshold = config.getDouble("lthreshold");
+ double uthreshold = config.getDouble("uthreshold");
+ double pmL1 = config.getDouble("pmL1");
+ double pmL2 = config.getDouble("pmL2");
+ double pmL3 = config.getDouble("pmL3");
+ int blocknx = config.getInt("blocknx");
+ double nx3 = config.getDouble("nx3");
+ double dp_LB = config.getDouble("dp_LB");
+ double nu_LB = config.getDouble("nu_LB");
+ string timeSeriesFile = config.getString("timeSeriesFile");
+ double restartStep = config.getDouble("restartStep");
+ double restartStepStart = config.getDouble("restartStepStart");
+ double endTime = config.getDouble("endTime");
+ double outTime = config.getDouble("outTime");
+ double availMem = config.getDouble("availMem");
+ bool rawFile = config.getBool("rawFile");
+ double timeSeriesOutTime = config.getDouble("timeSeriesOutTime");
+ bool logToFile = config.getBool("logToFile");
+ bool spongeLayer = config.getBool("spongeLayer");
+
+
+ CommunicatorPtr comm = MPICommunicator::getInstance();
+ int myid = comm->getProcessID();
+
+ if (logToFile)
+ {
+#if defined(__unix__)
+ if (myid == 0)
+ {
+ const char* str = pathname.c_str();
+ int status = mkdir(str, S_IRWXU | S_IRWXG | S_IROTH | S_IXOTH);
+ }
+#endif
+
+ if (myid == 0)
+ {
+ stringstream logFilename;
+ logFilename << pathname + "/logfile" + UbSystem::toString(UbSystem::getTimeStamp()) + ".txt";
+ UbLog::output_policy::setStream(logFilename.str());
+ }
+ }
+
+ //Sleep(30000);
+
+ if (myid == 0) UBLOG(logINFO, "Testcase permeability");
+
+ string machinename = UbSystem::getMachineName();
+ UBLOG(logINFO, "PID = " << myid << " Hostname: " << machinename);
+ UBLOG(logINFO, "PID = " << myid << " Total Physical Memory (RAM): " << Utilities::getTotalPhysMem());
+ UBLOG(logINFO, "PID = " << myid << " Physical Memory currently used: " << Utilities::getPhysMemUsed());
+ UBLOG(logINFO, "PID = " << myid << " Physical Memory currently used by current process: " << Utilities::getPhysMemUsedByMe());
+
+ int blocknx1 = blocknx;
+ int blocknx2 = blocknx;
+ int blocknx3 = blocknx;
+
+ LBMReal rho_LB = 0.0;
+ double rhoLBinflow = dp_LB*3.0;
+
+ LBMUnitConverterPtr conv = LBMUnitConverterPtr(new LBMUnitConverter());
+
+ const int baseLevel = 0;
+
+ double coord[6];
+ double deltax;
+
+ Grid3DPtr grid(new Grid3D(comm));
+
+ //////////////////////////////////////////////////////////////////////////
+ //restart
+ UbSchedulerPtr rSch(new UbScheduler(restartStep, restartStepStart));
+ RestartCoProcessor rp(grid, rSch, comm, pathname, RestartCoProcessor::BINARY);
+ //////////////////////////////////////////////////////////////////////////
+
+ if (grid->getTimeStep() == 0)
+ {
+ if (myid == 0) UBLOG(logINFO, "new start..");
+
+ UBLOG(logINFO, "new start PID = " << myid << " Hostname: " << machinename);
+ UBLOG(logINFO, "new start PID = " << myid << " Total Physical Memory (RAM): " << Utilities::getTotalPhysMem());
+ UBLOG(logINFO, "new start PID = " << myid << " Physical Memory currently used: " << Utilities::getPhysMemUsed());
+ UBLOG(logINFO, "new start PID = " << myid << " Physical Memory currently used by current process: " << Utilities::getPhysMemUsedByMe());
+
+ string samplePathname = pathGeo + sampleFilename;
+
+ double deltaVoxelX1 = pmL1/(double)pmNX1;
+ double deltaVoxelX2 = pmL2/(double)pmNX2;
+ double deltaVoxelX3 = pmL3/(double)pmNX3;
+
+ GbVoxelMatrix3DPtr sample(new GbVoxelMatrix3D(pmNX1, pmNX2, pmNX3, 0, lthreshold, uthreshold));
+ if (rawFile)
+ {
+ sample->readMatrixFromRawFile<unsigned short>(samplePathname, GbVoxelMatrix3D::BigEndian);
+ }
+ else
+ {
+ sample->readMatrixFromVtiASCIIFile(samplePathname);
+ }
+
+ sample->setVoxelMatrixDelta((float)deltaVoxelX1, (float)deltaVoxelX2, (float)deltaVoxelX3);
+ sample->setVoxelMatrixMininum(0.0, 0.0, 0.0);
+
+ if (myid == 0) sample->writeToVTKImageDataASCII(pathname + "/geo/sample");
+
+ ///////////////////////////////////////////////////////
+
+ ////////////////////////////////////////////////////////////////////////
+
+ double offset1 = sample->getLengthX1()/10.0;
+ double offset2 = 2.0*offset1;
+ //bounding box
+ double g_minX1 = sample->getX1Minimum() - offset1;
+ double g_minX2 = sample->getX2Minimum();
+ double g_minX3 = sample->getX3Minimum();
+
+ double g_maxX1 = sample->getX1Maximum() + offset2;
+ double g_maxX2 = sample->getX2Maximum();
+ double g_maxX3 = sample->getX3Maximum();
+
+ deltax = (g_maxX3-g_minX3) /(nx3*blocknx3);
+
+ double blockLength = (double)blocknx1*deltax;
+
+ grid->setPeriodicX1(false);
+ grid->setPeriodicX2(false);
+ grid->setPeriodicX3(false);
+ grid->setDeltaX(deltax);
+ grid->setBlockNX(blocknx1, blocknx2, blocknx3);
+
+ GbObject3DPtr gridCube(new GbCuboid3D(g_minX1, g_minX2, g_minX3, g_maxX1, g_maxX2, g_maxX3));
+ if (myid == 0) GbSystem3D::writeGeoObject(gridCube.get(), pathname + "/geo/gridCube", WbWriterVtkXmlBinary::getInstance());
+
+ GenBlocksGridVisitor genBlocks(gridCube);
+ grid->accept(genBlocks);
+
+ if (myid == 0)
+ {
+ UBLOG(logINFO, "Parameters:");
+ UBLOG(logINFO, "rho_LB = " << rho_LB);
+ UBLOG(logINFO, "nu_LB = " << nu_LB);
+ UBLOG(logINFO, "dp_LB = " << dp_LB);
+ UBLOG(logINFO, "dx = " << deltax << " m");
+ UBLOG(logINFO, "numOfThreads = " << numOfThreads);
+ UBLOG(logINFO, "path = " << pathname);
+ UBLOG(logINFO, "Preprozess - start");
+ }
+
+ //walls
+ GbCuboid3DPtr addWallYmin(new GbCuboid3D(g_minX1-blockLength, g_minX2-blockLength, g_minX3-blockLength, g_maxX1+blockLength, g_minX2, g_maxX3+blockLength));
+ if (myid == 0) GbSystem3D::writeGeoObject(addWallYmin.get(), pathname+"/geo/addWallYmin", WbWriterVtkXmlASCII::getInstance());
+
+ GbCuboid3DPtr addWallZmin(new GbCuboid3D(g_minX1-blockLength, g_minX2-blockLength, g_minX3-blockLength, g_maxX1+blockLength, g_maxX2+blockLength, g_minX3));
+ if (myid == 0) GbSystem3D::writeGeoObject(addWallZmin.get(), pathname+"/geo/addWallZmin", WbWriterVtkXmlASCII::getInstance());
+
+ GbCuboid3DPtr addWallYmax(new GbCuboid3D(g_minX1-blockLength, g_maxX2, g_minX3-blockLength, g_maxX1+blockLength, g_maxX2+blockLength, g_maxX3+blockLength));
+ if (myid == 0) GbSystem3D::writeGeoObject(addWallYmax.get(), pathname+"/geo/addWallYmax", WbWriterVtkXmlASCII::getInstance());
+
+ GbCuboid3DPtr addWallZmax(new GbCuboid3D(g_minX1-blockLength, g_minX2-blockLength, g_maxX3, g_maxX1+blockLength, g_maxX2+blockLength, g_maxX3+blockLength));
+ if (myid == 0) GbSystem3D::writeGeoObject(addWallZmax.get(), pathname+"/geo/addWallZmax", WbWriterVtkXmlASCII::getInstance());
+
+
+ //inflow
+ GbCuboid3DPtr geoInflow(new GbCuboid3D(g_minX1-blockLength, g_minX2-blockLength, g_minX3-blockLength, g_minX1, g_maxX2+blockLength, g_maxX3+blockLength));
+ if (myid == 0) GbSystem3D::writeGeoObject(geoInflow.get(), pathname + "/geo/geoInflow", WbWriterVtkXmlASCII::getInstance());
+
+ //outflow
+ GbCuboid3DPtr geoOutflow(new GbCuboid3D(g_maxX1, g_minX2-blockLength, g_minX3-blockLength, g_maxX1+blockLength, g_maxX2+blockLength, g_maxX3+blockLength));
+ if (myid == 0) GbSystem3D::writeGeoObject(geoOutflow.get(), pathname + "/geo/geoOutflow", WbWriterVtkXmlASCII::getInstance());
+
+ WriteBlocksCoProcessorPtr ppblocks(new WriteBlocksCoProcessor(grid, UbSchedulerPtr(new UbScheduler(1)), pathname, WbWriterVtkXmlBinary::getInstance(), comm));
+
+ //bone interactor
+ int bcOptionNoSlip = 1; //0=simple Bounce Back, 1=quadr. BB, 2=thin wall
+ D3Q27BoundaryConditionAdapterPtr bcNoSlip(new D3Q27NoSlipBCAdapter(bcOptionNoSlip));
+ D3Q27InteractorPtr sampleInt(new D3Q27Interactor(sample, grid, bcNoSlip, Interactor3D::SOLID));
+
+ //wall interactors
+ D3Q27InteractorPtr addWallYminInt(new D3Q27Interactor(addWallYmin, grid, bcNoSlip, Interactor3D::SOLID));
+ D3Q27InteractorPtr addWallZminInt(new D3Q27Interactor(addWallZmin, grid, bcNoSlip, Interactor3D::SOLID));
+ D3Q27InteractorPtr addWallYmaxInt(new D3Q27Interactor(addWallYmax, grid, bcNoSlip, Interactor3D::SOLID));
+ D3Q27InteractorPtr addWallZmaxInt(new D3Q27Interactor(addWallZmax, grid, bcNoSlip, Interactor3D::SOLID));
+
+ D3Q27BoundaryConditionAdapterPtr denBCAdapterInflow(new D3Q27DensityBCAdapter(rhoLBinflow));
+ denBCAdapterInflow->setSecondaryBcOption(0);
+ D3Q27InteractorPtr inflowInt = D3Q27InteractorPtr(new D3Q27Interactor(geoInflow, grid, denBCAdapterInflow, Interactor3D::SOLID));
+
+ //outflow
+ D3Q27BoundaryConditionAdapterPtr denBCAdapterOutflow(new D3Q27DensityBCAdapter(rho_LB));
+ denBCAdapterOutflow->setSecondaryBcOption(0);
+ D3Q27InteractorPtr outflowInt = D3Q27InteractorPtr(new D3Q27Interactor(geoOutflow, grid, denBCAdapterOutflow, Interactor3D::SOLID));
+
+
+ UBLOG(logINFO, "PID = " << myid << " Hostname: " << machinename);
+ UBLOG(logINFO, "PID = " << myid << " Total Physical Memory (RAM): " << Utilities::getTotalPhysMem());
+ UBLOG(logINFO, "PID = " << myid << " Physical Memory currently used: " << Utilities::getPhysMemUsed());
+ UBLOG(logINFO, "PID = " << myid << " Physical Memory currently used by current process: " << Utilities::getPhysMemUsedByMe());
+
+ ////////////////////////////////////////////
+ //METIS
+ Grid3DVisitorPtr metisVisitor(new MetisPartitioningGridVisitor(comm, MetisPartitioningGridVisitor::LevelBased, D3Q27System::BSW, MetisPartitioner::RECURSIVE));
+ ////////////////////////////////////////////
+ //Zoltan
+ //Grid3DVisitorPtr zoltanVisitor(new ZoltanPartitioningGridVisitor(comm, D3Q27System::BSW, 1));
+ //grid->accept(zoltanVisitor);
+ /////delete solid blocks
+ if (myid == 0) UBLOG(logINFO, "deleteSolidBlocks - start");
+ InteractorsHelper intHelper(grid, metisVisitor);
+ intHelper.addInteractor(addWallYminInt);
+ intHelper.addInteractor(addWallZminInt);
+ intHelper.addInteractor(addWallYmaxInt);
+ intHelper.addInteractor(addWallZmaxInt);
+ intHelper.addInteractor(inflowInt);
+ intHelper.addInteractor(outflowInt);
+ intHelper.addInteractor(sampleInt);
+ intHelper.selectBlocks();
+ if (myid == 0) UBLOG(logINFO, "deleteSolidBlocks - end");
+ //////////////////////////////////////
+
+ //set connectors
+ D3Q27InterpolationProcessorPtr iProcessor(new D3Q27IncompressibleOffsetInterpolationProcessor());
+ D3Q27SetConnectorsBlockVisitor setConnsVisitor(comm, true, D3Q27System::ENDDIR, nu_LB, iProcessor);
+ grid->accept(setConnsVisitor);
+
+ //domain decomposition for threads
+ PQueuePartitioningGridVisitor pqPartVisitor(numOfThreads);
+ grid->accept(pqPartVisitor);
+
+ ppblocks->process(0);
+ ppblocks.reset();
+
+ unsigned long nob = grid->getNumberOfBlocks();
+ int gl = 3;
+ unsigned long nodb = (blocknx1)* (blocknx2)* (blocknx3);
+ unsigned long nod = nob * (blocknx1)* (blocknx2)* (blocknx3);
+ unsigned long nodg = nob * (blocknx1 + gl) * (blocknx2 + gl) * (blocknx3 + gl);
+ double needMemAll = double(nodg*(27 * sizeof(double) + sizeof(int) + sizeof(float) * 4));
+ double needMem = needMemAll / double(comm->getNumberOfProcesses());
+
+ if (myid == 0)
+ {
+ UBLOG(logINFO, "Number of blocks = " << nob);
+ UBLOG(logINFO, "Number of nodes = " << nod);
+ int minInitLevel = grid->getCoarsestInitializedLevel();
+ int maxInitLevel = grid->getFinestInitializedLevel();
+ for (int level = minInitLevel; level <= maxInitLevel; level++)
+ {
+ int nobl = grid->getNumberOfBlocks(level);
+ UBLOG(logINFO, "Number of blocks for level " << level << " = " << nobl);
+ UBLOG(logINFO, "Number of nodes for level " << level << " = " << nobl*nodb);
+ }
+ UBLOG(logINFO, "Necessary memory = " << needMemAll << " bytes");
+ UBLOG(logINFO, "Necessary memory per process = " << needMem << " bytes");
+ UBLOG(logINFO, "Available memory per process = " << availMem << " bytes");
+ }
+
+ LBMKernel3DPtr kernel;
+ kernel = LBMKernel3DPtr(new LBMKernelETD3Q27CCLB(blocknx1, blocknx2, blocknx3, LBMKernelETD3Q27CCLB::NORMAL));
+
+ BCProcessorPtr bcProc(new D3Q27ETBCProcessor());
+
+ BoundaryConditionPtr densityBC(new NonEqDensityBoundaryCondition());
+ BoundaryConditionPtr noSlipBC(new NoSlipBoundaryCondition());
+
+ bcProc->addBC(densityBC);
+ bcProc->addBC(noSlipBC);
+
+ kernel->setBCProcessor(bcProc);
+
+ SetKernelBlockVisitor kernelVisitor(kernel, nu_LB, availMem, needMem);
+ grid->accept(kernelVisitor);
+
+ //BC
+ intHelper.setBC();
+
+ BoundaryConditionBlockVisitor bcVisitor;
+ grid->accept(bcVisitor);
+
+ //Press*1.6e8+(14.76-coordsX)/3.5*5000
+ //initialization of distributions
+ mu::Parser fct;
+ fct.SetExpr("(x1max-x1)/l*dp*3.0");
+ fct.DefineConst("dp", dp_LB);
+ fct.DefineConst("x1max", g_maxX1);
+ fct.DefineConst("l", g_maxX1-g_minX1);
+
+ D3Q27ETInitDistributionsBlockVisitor initVisitor(nu_LB, rho_LB);
+ initVisitor.setRho(fct);
+ grid->accept(initVisitor);
+
+ //Postrozess
+ UbSchedulerPtr geoSch(new UbScheduler(1));
+ MacroscopicQuantitiesCoProcessorPtr ppgeo(
+ new MacroscopicQuantitiesCoProcessor(grid, geoSch, pathname, WbWriterVtkXmlBinary::getInstance(), conv, true));
+ ppgeo->process(0);
+ ppgeo.reset();
+
+ coord[0] = sample->getX1Minimum();
+ coord[1] = sample->getX2Minimum();
+ coord[2] = sample->getX3Minimum();
+ coord[3] = sample->getX1Maximum();
+ coord[4] = sample->getX2Maximum();
+ coord[5] = sample->getX3Maximum();
+
+ ////////////////////////////////////////////////////////
+ FILE * pFile;
+ string str = pathname + "/checkpoints/coord.txt";
+ pFile = fopen(str.c_str(), "w");
+ fprintf(pFile, "%g\n", deltax);
+ fprintf(pFile, "%g\n", coord[0]);
+ fprintf(pFile, "%g\n", coord[1]);
+ fprintf(pFile, "%g\n", coord[2]);
+ fprintf(pFile, "%g\n", coord[3]);
+ fprintf(pFile, "%g\n", coord[4]);
+ fprintf(pFile, "%g\n", coord[5]);
+ fclose(pFile);
+ ////////////////////////////////////////////////////////
+
+ grid->addInteractor(inflowInt);
+
+ if (myid == 0) UBLOG(logINFO, "Preprozess - end");
+ }
+ else
+ {
+ ////////////////////////////////////////////////////////
+ FILE * pFile;
+ string str = pathname + "/checkpoints/coord.txt";
+ pFile = fopen(str.c_str(), "r");
+ fscanf(pFile, "%lg\n", &deltax);
+ fscanf(pFile, "%lg\n", &coord[0]);
+ fscanf(pFile, "%lg\n", &coord[1]);
+ fscanf(pFile, "%lg\n", &coord[2]);
+ fscanf(pFile, "%lg\n", &coord[3]);
+ fscanf(pFile, "%lg\n", &coord[4]);
+ fscanf(pFile, "%lg\n", &coord[5]);
+ fclose(pFile);
+ ////////////////////////////////////////////////////////
+
+ //new nu
+ //ViscosityBlockVisitor nuVisitor(nu_LB);
+ //grid->accept(nuVisitor);
+
+ //new dp
+ Grid3D::Interactor3DSet interactors = grid->getInteractors();
+ interactors[0]->setGrid3D(grid);
+ boost::dynamic_pointer_cast<D3Q27Interactor>(interactors[0])->deleteBCAdapter();
+ D3Q27BoundaryConditionAdapterPtr denBCAdapterFront(new D3Q27DensityBCAdapter(rhoLBinflow));
+ boost::dynamic_pointer_cast<D3Q27Interactor>(interactors[0])->addBCAdapter(denBCAdapterFront);
+ interactors[0]->updateInteractor();
+
+ if (myid == 0)
+ {
+ UBLOG(logINFO, "Parameters:");
+ UBLOG(logINFO, "rho_LB = " << rho_LB);
+ UBLOG(logINFO, "nu_LB = " << nu_LB);
+ UBLOG(logINFO, "dp_LB = " << dp_LB);
+ UBLOG(logINFO, "dx = " << deltax << " m");
+ }
+
+ BoundaryConditionBlockVisitor bcVisitor;
+ grid->accept(bcVisitor);
+
+ //set connectors
+ D3Q27InterpolationProcessorPtr iProcessor(new D3Q27IncompressibleOffsetInterpolationProcessor());
+ D3Q27SetConnectorsBlockVisitor setConnsVisitor(comm, true, D3Q27System::ENDDIR, nu_LB, iProcessor);
+ grid->accept(setConnsVisitor);
+
+ if (myid == 0) UBLOG(logINFO, "Restart - end");
+ }
+ UbSchedulerPtr nupsSch(new UbScheduler(10, 30, 100));
+ //nupsSch->addSchedule(nupsStep[0], nupsStep[1], nupsStep[2]);
+ NUPSCounterCoProcessor npr(grid, nupsSch, numOfThreads, comm);
+
+ UbSchedulerPtr stepSch(new UbScheduler(outTime));
+
+ MacroscopicQuantitiesCoProcessor pp(grid, stepSch, pathname, WbWriterVtkXmlBinary::getInstance(), conv);
+
+ deltax = grid->getDeltaX(baseLevel);
+ double dxd2 = deltax / 2.0;
+
+ D3Q27IntegrateValuesHelperPtr ih1(new D3Q27IntegrateValuesHelper(grid, comm, coord[0] - dxd2*10.0, coord[1] - dxd2, coord[2] - dxd2,
+ coord[0] - dxd2*10.0 - 2.0*dxd2, coord[4] + dxd2, coord[5] + dxd2));
+
+ //D3Q27IntegrateValuesHelperPtr ih2(new D3Q27IntegrateValuesHelper(grid, comm, coord[3]/2.0, coord[1] - dxd2, coord[2] - dxd2,
+ // coord[3]/2.0 + 2.0*dxd2, coord[4] + dxd2, coord[5] + dxd2));
+ D3Q27IntegrateValuesHelperPtr ih2(new D3Q27IntegrateValuesHelper(grid, comm, coord[0], coord[1], coord[2], coord[3], coord[4], coord[5]));
+
+ D3Q27IntegrateValuesHelperPtr ih3(new D3Q27IntegrateValuesHelper(grid, comm, coord[3] + dxd2*10.0, coord[1] - dxd2, coord[2] - dxd2,
+ coord[3] + dxd2*10.0 + 2.0*dxd2, coord[4] + dxd2, coord[5] + dxd2));
+
+ //D3Q27IntegrateValuesHelperPtr ih1(new D3Q27IntegrateValuesHelper(grid, comm, coord[0], coord[1], coord[2], coord[3], coord[4], coord[5]));
+ if (myid == 0) GbSystem3D::writeGeoObject(ih1->getBoundingBox().get(), pathname + "/geo/ih1", WbWriterVtkXmlBinary::getInstance());
+ if (myid == 0) GbSystem3D::writeGeoObject(ih2->getBoundingBox().get(), pathname + "/geo/ih2", WbWriterVtkXmlBinary::getInstance());
+ if (myid == 0) GbSystem3D::writeGeoObject(ih3->getBoundingBox().get(), pathname + "/geo/ih3", WbWriterVtkXmlBinary::getInstance());
+
+ double factorp = 1; // dp_real / dp_LB;
+ double factorv = 1;// dx / dt;
+ UbSchedulerPtr stepMV(new UbScheduler(timeSeriesOutTime));
+
+ TimeseriesCoProcessor tsp1(grid, stepMV, ih1, pathname+timeSeriesFile+"_1", comm);
+ TimeseriesCoProcessor tsp2(grid, stepMV, ih2, pathname+timeSeriesFile+"_2", comm);
+ TimeseriesCoProcessor tsp3(grid, stepMV, ih3, pathname+timeSeriesFile+"_3", comm);
+
+ if (myid == 0)
+ {
+ UBLOG(logINFO, "PID = " << myid << " Total Physical Memory (RAM): " << Utilities::getTotalPhysMem());
+ UBLOG(logINFO, "PID = " << myid << " Physical Memory currently used: " << Utilities::getPhysMemUsed());
+ UBLOG(logINFO, "PID = " << myid << " Physical Memory currently used by current process: " << Utilities::getPhysMemUsedByMe());
+ }
+
+ CalculationManagerPtr calculation(new CalculationManager(grid, numOfThreads, endTime, stepMV));
+ if (myid == 0) UBLOG(logINFO, "Simulation-start");
+ calculation->calculate();
+ if (myid == 0) UBLOG(logINFO, "Simulation-end");
+ }
+ catch (exception& e)
+ {
+ cerr << e.what() << endl << flush;
+ }
+ catch (string& s)
+ {
+ cerr << s << endl;
+ }
+ catch (...)
+ {
+ cerr << "unknown exception" << endl;
+ }
+
+}
+//////////////////////////////////////////////////////////////////////////
+int main(int argc, char* argv[])
+{
+
+ if (argv!=NULL)
+ {
+ if (argv[1]!=NULL)
+ {
+ run(string(argv[1]));
+ }
+ else
+ {
+ cout<<"Configuration file is missing!"<<endl;
+ }
+ }
+
+ return 0;
+}
diff --git a/source/Applications/aperm/configBombadil2.txt b/source/Applications/aperm/configBombadil2.txt
new file mode 100644
index 0000000000000000000000000000000000000000..9a29f30c232af059660accf738a4259403d9d860
--- /dev/null
+++ b/source/Applications/aperm/configBombadil2.txt
@@ -0,0 +1,44 @@
+pathname = d:/temp/aperm
+pathGeo = d:/Projects/SFB880/GeometrienPoroeseMedien/isotrop/PA80-110
+numOfThreads = 4
+availMem = 3e9
+logToFile = false
+
+#poroeses Medium
+rawFile = false
+sampleFilename = /alu_80-110.vti
+
+#Diminsion in Voxel
+pmNX1 = 200
+pmNX2 = 200
+pmNX3 = 200
+
+#Threshold
+lthreshold = 29041
+uthreshold = 65535
+
+#Diminsion in m
+pmL1 = 0.726e-3
+pmL2 = 0.75e-3
+pmL3 = 0.786e-3
+
+#grid
+blocknx = 30
+nx3 = 5
+spongeLayer=true
+
+dp_LB = 0.001
+#nu_LB = 0.168666666667
+#nu_LB = 0.168666666667e-1
+#nu_LB = 0.168666666667e-2
+#nu_LB = 0.168666666667e-3
+nu_LB = 0.168666666667e-4
+
+timeSeriesFile = /timeseries/simAlu80_5
+timeSeriesOutTime = 10
+
+restartStep = 20000
+restartStepStart=20000
+
+endTime = 60000
+outTime = 10
\ No newline at end of file
diff --git a/source/Applications/aperm/configBombadilSBP120s500.txt b/source/Applications/aperm/configBombadilSBP120s500.txt
new file mode 100644
index 0000000000000000000000000000000000000000..4fb351a97ce8bb1aa040c3d1e9a904f6ca9a7e11
--- /dev/null
+++ b/source/Applications/aperm/configBombadilSBP120s500.txt
@@ -0,0 +1,50 @@
+#
+#Simulation parameters for determitatoin of permeability
+#SBP120
+
+pathname = d:/temp/perm
+pathGeo = d:/Projects/SFB880/GeometrienPoroeseMedien/isotrop/SBP120
+numOfThreads = 4
+availMem = 3e9
+logToFile = false
+
+#porous media
+rawFile = false
+sampleFilename = /SBP120s500_center_closed.vti
+
+#diminsions [voxel]
+pmNX1 = 500
+pmNX2 = 500
+pmNX3 = 500
+
+#threshold
+#lthreshold = 38370
+#uthreshold = 65535
+lthreshold = 1
+uthreshold = 1
+
+
+#diminsions [m]
+pmL1 = 1.87e-3
+pmL2 = 1.87e-3
+pmL3 = 1.87e-3
+
+#grid
+#blocknx = 30
+#nx3 = 5
+blocknx = 50
+nx3 = 10
+spongeLayer=true
+
+#physic
+dp_LB = 1e-7
+nu_LB = 0.01
+
+timeSeriesFile = /timeseries/simSBP120_1
+timeSeriesOutTime = 10
+
+restartStep = 20000
+restartStepStart=20000
+
+endTime = 60000
+outTime = 100
diff --git a/source/Applications/aperm/config_HLRS_SBP120.cfg b/source/Applications/aperm/config_HLRS_SBP120.cfg
new file mode 100644
index 0000000000000000000000000000000000000000..54ceee42016bd06904a472fedccb35c68bbb06be
--- /dev/null
+++ b/source/Applications/aperm/config_HLRS_SBP120.cfg
@@ -0,0 +1,43 @@
+#HLRS
+#Simulation parameters for determitatoin of permeability
+#SBP120
+
+pathname = /univ_1/ws1/ws/xrmkuchr-perm-0/SBP120
+pathGeo = /univ_1/ws1/ws/xrmkuchr-perm-0/SBP120/Data
+numOfThreads = 24
+availMem = 128e9
+logToFile = true
+#porous media
+rawFile = false
+sampleFilename = /Sinterbronze_SBP120_1358x1376x1572.raw
+
+#diminsions [voxel]
+pmNX1 = 1358
+pmNX2 = 1376
+pmNX3 = 1572
+
+#threshold
+lthreshold = 38370
+uthreshold = 65535
+
+#diminsions [m]
+pmL1 = 5092499.73e-9
+pmL2 = 5159999.85e-9
+pmL3 = 5894999.98e-9
+
+#grid
+blocknx = 64
+nx3 = 22
+
+#physic
+dp_LB = 1e-7
+nu_LB = 0.01
+
+timeSeriesFile = /timeseries/simSBP120_1
+timeSeriesOutTime = 100
+
+restartStep = 1000
+restartStepStart=1000
+
+endTime = 2000
+outTime = 1000
diff --git a/source/Applications/pChannel/configBombadilpChannel.cfg b/source/Applications/pChannel/configBombadilpChannel.cfg
index 8bd0735882a22985395dc17a9034a26830ebebb9..2482d5526566727c7016faed898a9d172c63a4fe 100644
--- a/source/Applications/pChannel/configBombadilpChannel.cfg
+++ b/source/Applications/pChannel/configBombadilpChannel.cfg
@@ -2,9 +2,9 @@
#Simulation parameters for porous channel
#
-pathname = d:/temp/ChannelFlowC2
+pathname = d:/temp/ChannelFlowCorr2
pathGeo = d:/Projects/SFB880/GeometrienPoroeseMedien/PA80-110
-numOfThreads = 4
+numOfThreads = 1
availMem = 8e9
logToFile = false
@@ -51,7 +51,7 @@ voxelDeltaX = 3.6496350365e-6 3.76789751319e-6 3.95256916996e-6
#pmL = 4e-3 0.5e-3 1e-3
#pmL = 1e-3 1e-3 1e-3
-pmL = 4e-3 4e-3 1e-3
+pmL = 4e-3 8e-3 1e-3
#grid
refineLevel = 0
@@ -60,7 +60,7 @@ refineLevel = 0
deltaXfine = 40e-5
#blocknx = 10 10 10
-blocknx = 20 40 20
+blocknx = 20 20 20
#lengthFactor = 6
lengthFactor = 16
thinWall = false
@@ -84,15 +84,15 @@ u_LB = 0.1
nupsSteps=1000
-restartStep = 21000000
-restartStepStart = 21000000
+restartStep = 20000
+restartStepStart = 20000
averaging = false
averagingReset = false
timeAvStart = 21000000
timeAvStop = 2100010000
-endTime = 1000000
-outTime = 1000
+endTime = 5
+outTime = 1
diff --git a/source/Applications/pChannel/pChannel.cpp b/source/Applications/pChannel/pChannel.cpp
index c7636c3c5709e01592978ba53f82f6046d788c0f..c6aab304ee6a33e8c3bf140c4fd577edf4cd5992 100644
--- a/source/Applications/pChannel/pChannel.cpp
+++ b/source/Applications/pChannel/pChannel.cpp
@@ -84,7 +84,7 @@ void run(string configname)
}
}
- //Sleep(30000);
+ Sleep(30000);
if (myid == 0) UBLOG(logINFO, "Testcase porous channel");
@@ -389,10 +389,11 @@ void run(string configname)
grid->accept(bcVisitor);
mu::Parser inflowProfileVx1, inflowProfileVx2, inflowProfileVx3, inflowProfileRho;
- // inflowProfile.SetExpr("x3 < h ? 0.0 : uLB+1*x1-1*x2");
+ inflowProfileVx1.SetExpr("x3 < h ? 0.0 : uLB+1*x1");
+ //inflowProfileVx1.SetExpr("x3 < h ? 0.0 : uLB+1*x1-1*x2");
////inflowProfile.SetExpr("uLB+1*x1-1*x2");
// //inflowProfile.SetExpr("uLB");
- // inflowProfile.DefineConst("uLB", u_LB);
+ inflowProfileVx1.DefineConst("uLB", u_LB);
// //inflowProfile.DefineConst("uLB", 0.0116);
// inflowProfile.DefineConst("h", pmL[2]);
@@ -425,10 +426,10 @@ void run(string configname)
//inflowProfile.SetExpr("Uref/log((Href+z0)/z0)*log((x3-zg+z0)/z0)");
//inflowProfileVx1.SetExpr("x3 > 0 && (zMax-x3) > 0 ? (x3 < h ? Uref/log((Href+z0)/z0)*log((x3-zg+z0)/z0) : Uref/log((Href+z0)/z0)*log((zMax-x3-zg+z0)/z0)) : 0");
- inflowProfileVx1.DefineFun("rangeRandom1", rangeRandom1);
+ //inflowProfileVx1.DefineFun("rangeRandom1", rangeRandom1);
//inflowProfile.SetExpr("x3 < h ? Uref/log((Href+z0)/z0)*log((x3-zg+z0)/z0)+rangeRandom(-nois, nois) : Uref/log((Href+z0)/z0)*log((zMax-x3-zg+z0)/z0)+rangeRandom(-nois, nois)");
- inflowProfileVx1.SetExpr("x3 < h ? Uref/log((Href+z0)/z0)*log((x3-zg+z0)/z0)+0.1*Uref*rangeRandom1() : Uref/log((Href+z0)/z0)*log((zMax-x3-zg+z0)/z0)+0.1*Uref*rangeRandom1()");
+ //inflowProfileVx1.SetExpr("x3 < h ? Uref/log((Href+z0)/z0)*log((x3-zg+z0)/z0)+0.1*Uref*rangeRandom1() : Uref/log((Href+z0)/z0)*log((zMax-x3-zg+z0)/z0)+0.1*Uref*rangeRandom1()");
//inflowProfileVx1.SetExpr("x3 < h ? Uref/log((Href+z0)/z0)*log((x3-zg+z0)/z0)+ U*cos(8.0*PI*(x1)/(L1))*sin(8.0*PI*(x3)/L3) : Uref/log((Href+z0)/z0)*log((zMax-x3-zg+z0)/z0)+U*cos(8.0*PI*(x1)/(L1))*sin(8.0*PI*(x3)/L3)");
//inflowProfileVx1.SetExpr("U*cos(4.0*PI*(x1)/(L1))*sin(4.0*PI*(x3)/L3)");
@@ -494,8 +495,8 @@ void run(string configname)
D3Q27ETInitDistributionsBlockVisitor initVisitor(nu_LB, rho_LB);
initVisitor.setVx1(inflowProfileVx1);
- initVisitor.setVx2(inflowProfileVx2);
- initVisitor.setVx3(inflowProfileVx3);
+ //initVisitor.setVx2(inflowProfileVx2);
+ //initVisitor.setVx3(inflowProfileVx3);
//initVisitor.setRho(inflowProfileRho);
grid->accept(initVisitor);
@@ -747,8 +748,17 @@ void run(string configname)
//UbSchedulerPtr catalystSch(new UbScheduler(1));
//InSituCatalystCoProcessor catalyst(grid, catalystSch, "pchannel.py");
- UbSchedulerPtr exitSch(new UbScheduler(10));
- EmergencyExitCoProcessor exitCoProc(grid, exitSch, pathname, RestartCoProcessorPtr(&rp), comm);
+ //UbSchedulerPtr exitSch(new UbScheduler(10));
+ //EmergencyExitCoProcessor exitCoProc(grid, exitSch, pathname, RestartCoProcessorPtr(&rp), comm);
+
+ //create line time series
+ UbSchedulerPtr tpcSch(new UbScheduler(1,1,3));
+ //GbPoint3DPtr p1(new GbPoint3D(0.0,0.005,0.01));
+ //GbPoint3DPtr p2(new GbPoint3D(0.064,0.005,0.01));
+ //GbLine3DPtr line(new GbLine3D(p1.get(),p2.get()));
+ GbLine3DPtr line(new GbLine3D(new GbPoint3D(0.0,0.005,0.01),new GbPoint3D(0.064,0.005,0.01)));
+ LineTimeSeriesCoProcessor lineTs(grid, tpcSch,pathname+"/TimeSeries/line1.csv",line, 0,comm);
+ if (myid==0) lineTs.writeLine(pathname+"/geo/line1");
if (myid == 0)
{
diff --git a/source/VirtualFluids.h b/source/VirtualFluids.h
index d56163928d7ab8e04ae37328c581d7448205be63..7549b401dc6fff84b2e79c41606481df04da9474 100644
--- a/source/VirtualFluids.h
+++ b/source/VirtualFluids.h
@@ -157,6 +157,7 @@
//#include <CoProcessors/MeanValuesCoProcessor.h>
#include <CoProcessors/TimeAveragedValuesCoProcessor.h>
#include <CoProcessors/InSituCatalystCoProcessor.h>
+#include <LineTimeSeriesCoProcessor.h>
#include <LBM/D3Q27CompactInterpolationProcessor.h>
#include <LBM/D3Q27IncompressibleOffsetInterpolationProcessor.h>
#include <LBM/D3Q27IntegrateValuesHelper.h>
@@ -171,6 +172,7 @@
#include <LBM/LBMKernelETD3Q27.h>
#include <LBM/LBMKernelETD3Q27CCLB.h>
#include <LBM/CompressibleCumulantLBMKernel.h>
+#include <LBM/InitDensityLBMKernel.h>
#include <LBM/LBMSystem.h>
#include <LBM/LBMSystems.h>
#include <LBM/LBMUnitConverter.h>
diff --git a/source/VirtualFluidsCore/BoundaryCondition/D3Q27ETBCProcessor.cpp b/source/VirtualFluidsCore/BoundaryCondition/D3Q27ETBCProcessor.cpp
index 8992cdecd4cf249b407f831771fac3069980b77b..1626c5b907207823d31f1cbcd60e6c4fced33911 100644
--- a/source/VirtualFluidsCore/BoundaryCondition/D3Q27ETBCProcessor.cpp
+++ b/source/VirtualFluidsCore/BoundaryCondition/D3Q27ETBCProcessor.cpp
@@ -9,17 +9,7 @@ D3Q27ETBCProcessor::D3Q27ETBCProcessor()
D3Q27ETBCProcessor::D3Q27ETBCProcessor(LBMKernel3DPtr kernel)
{
DistributionArray3DPtr distributions = boost::dynamic_pointer_cast<EsoTwist3D>(kernel->getDataSet()->getFdistributions());
- //ghostLayerWitdh = kernel->getGhostLayerWidth();
- //collFactor = kernel->getCollisionFactor();
bcArray.resize( distributions->getNX1(), distributions->getNX2(), distributions->getNX3(), BCArray3D<D3Q27BoundaryCondition>::FLUID);
- //this->compressible = kernel->getCompressible();
-
- //minX1 = 0;
- //minX2 = 0;
- //minX3 = 0;
- //maxX1 = (int)bcArray.getNX1();
- //maxX2 = (int)bcArray.getNX2();
- //maxX3 = (int)bcArray.getNX3();
}
//////////////////////////////////////////////////////////////////////////
D3Q27ETBCProcessor::~D3Q27ETBCProcessor()
diff --git a/source/VirtualFluidsCore/BoundaryCondition/D3Q27ETBCProcessor.h b/source/VirtualFluidsCore/BoundaryCondition/D3Q27ETBCProcessor.h
index 2599851ec7093ce3a51c83f0202af6f3528a67d0..f1f459be7f20b186dd8afed4e5fed1f1ab8febb9 100644
--- a/source/VirtualFluidsCore/BoundaryCondition/D3Q27ETBCProcessor.h
+++ b/source/VirtualFluidsCore/BoundaryCondition/D3Q27ETBCProcessor.h
@@ -33,32 +33,8 @@ public:
protected:
std::vector<BoundaryConditionPtr> preBC;
std::vector<BoundaryConditionPtr> postBC;
-
- //int minX1, minX2, minX3, maxX1, maxX2, maxX3;
- //D3Q27BoundaryConditionPtr bcPtr;
- //LBMReal f[D3Q27System::ENDF+1];
- //LBMReal ftemp[D3Q27System::ENDF+1];
- //LBMReal feq[D3Q27System::ENDF+1];
- //LBMReal rho, vx1, vx2, vx3, rho0;
-
- //void init(LBMKernel3DPtr kernel);
BCArray3D<D3Q27BoundaryCondition> bcArray;
- //EsoTwist3DPtr distributions;
- //int ghostLayerWitdh;
- //LBMReal collFactor;
- //bool compressible;
- //EsoTwist3DPtr distributionsTemp;
-
- //std::vector <int> sizeVector;
- //std::vector <int> nodeVector;
- //std::vector <D3Q27BoundaryConditionPtr> bcVector;
- //typedef void(*CalcMacrosFct)(const LBMReal* const& /*f[27]*/, LBMReal& /*rho*/, LBMReal& /*vx1*/, LBMReal& /*vx2*/, LBMReal& /*vx3*/);
- //typedef LBMReal(*CalcFeqForDirFct)(const int& /*direction*/, const LBMReal& /*(d)rho*/, const LBMReal& /*vx1*/, const LBMReal& /*vx2*/, const LBMReal& /*vx3*/);
- //typedef void(*CalcFeqFct)(LBMReal* const& /*feq/*[27]*/, const LBMReal& /*rho*/, const LBMReal& /*vx1*/, const LBMReal& /*vx2*/, const LBMReal& /*vx3*/);
- //CalcFeqForDirFct calcFeqsForDirFct;
- //CalcMacrosFct calcMacrosFct;
- //CalcFeqFct calcFeqFct;
private:
friend class boost::serialization::access;
template<class Archive>
@@ -66,17 +42,6 @@ private:
{
ar & boost::serialization::base_object<BCProcessor>(*this);
ar & bcArray;
- //ar & distributions;
- //ar & ghostLayerWitdh;
- //ar & collFactor;
- //ar & compressible;
- //ar & minX1;
- //ar & minX2;
- //ar & minX3;
- //ar & maxX1;
- //ar & maxX2;
- //ar & maxX3;
- //ar & distributionsTemp;
ar & preBC;
ar & postBC;
}
diff --git a/source/VirtualFluidsCore/CoProcessors/AdjustForcingCoProcessor.cpp b/source/VirtualFluidsCore/CoProcessors/AdjustForcingCoProcessor.cpp
index 98c20a4ac2af225a49a80432a3d1e4cab4782eca..3ae67b046856c2aa5e700fd8d12c9f214c11030c 100644
--- a/source/VirtualFluidsCore/CoProcessors/AdjustForcingCoProcessor.cpp
+++ b/source/VirtualFluidsCore/CoProcessors/AdjustForcingCoProcessor.cpp
@@ -72,8 +72,8 @@ AdjustForcingCoProcessor::AdjustForcingCoProcessor(Grid3DPtr grid, UbSchedulerPt
if (istr2)
{
istr2 >> forcing;
- istr2 >> esum;
- istr2 >> eold;
+ //istr2 >> esum;
+ //istr2 >> eold;
}
istr2.close();
}
@@ -156,6 +156,7 @@ void AdjustForcingCoProcessor::collectData(double step)
//UBLOG(logINFO, "D3Q27AdjustForcingCoProcessor step: " << static_cast<int>(step));
//UBLOG(logINFO, "new forcing is: " << forcing);
std::string fname = path + "/forcing/forcing.csv";
+ //std::string fname = path + "/forcing/forcing_"+UbSystem::toString(comm->getProcessID())+".csv";
std::ofstream ostr;
ostr.open(fname.c_str(), std::ios_base::out | std::ios_base::app);
if (!ostr)
diff --git a/source/VirtualFluidsCore/CoProcessors/CalculateForcesCoProcessor.cpp b/source/VirtualFluidsCore/CoProcessors/CalculateForcesCoProcessor.cpp
index cfa597bd52942a636a3063d77c3e0857319933be..52ad99b2b8cb4e21162d0d06eaaa8f37d6edccfe 100644
--- a/source/VirtualFluidsCore/CoProcessors/CalculateForcesCoProcessor.cpp
+++ b/source/VirtualFluidsCore/CoProcessors/CalculateForcesCoProcessor.cpp
@@ -102,12 +102,13 @@ void CalculateForcesCoProcessor::calculateForces()
DistributionArray3DPtr distributions = kernel->getDataSet()->getFdistributions();
distributions->swap();
- int minX1 = kernel->getGhostLayerWidth();
- int maxX1 = (int)bcArray.getNX1() - 1 - kernel->getGhostLayerWidth();
- int minX2 = kernel->getGhostLayerWidth();
- int maxX2 = (int)bcArray.getNX2() - 1 - kernel->getGhostLayerWidth();
- int minX3 = kernel->getGhostLayerWidth();
- int maxX3 = (int)bcArray.getNX3() - 1 - kernel->getGhostLayerWidth();
+ int ghostLayerWidth = kernel->getGhostLayerWidth();
+ int minX1 = ghostLayerWidth;
+ int maxX1 = (int)bcArray.getNX1() - 1 - ghostLayerWidth;
+ int minX2 = ghostLayerWidth;
+ int maxX2 = (int)bcArray.getNX2() - 1 - ghostLayerWidth;
+ int minX3 = ghostLayerWidth;
+ int maxX3 = (int)bcArray.getNX3() - 1 - ghostLayerWidth;
BOOST_FOREACH(std::vector<int> node, transNodeIndicesSet)
{
diff --git a/source/VirtualFluidsCore/CoProcessors/LineTimeSeriesCoProcessor.cpp b/source/VirtualFluidsCore/CoProcessors/LineTimeSeriesCoProcessor.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..53a05c5e094ef338c021113321ccad80cb677e64
--- /dev/null
+++ b/source/VirtualFluidsCore/CoProcessors/LineTimeSeriesCoProcessor.cpp
@@ -0,0 +1,237 @@
+#include "LineTimeSeriesCoProcessor.h"
+#include "D3Q27ETBCProcessor.h"
+#include "WbWriterVtkXmlASCII.h"
+
+LineTimeSeriesCoProcessor::LineTimeSeriesCoProcessor(Grid3DPtr grid, UbSchedulerPtr s, const std::string& path, GbLine3DPtr line, int level, CommunicatorPtr comm) :
+ CoProcessor(grid, s),
+ path(path),
+ length(0),
+ ix1(0),
+ ix2(0),
+ ix3(0),
+ level(level),
+ line(line)
+{
+ root = comm->isRoot();
+ fname = path;
+
+ mpi_comm = *((MPI_Comm*)comm->getNativeCommunicator());
+ numOfProc = comm->getNumberOfProcesses();
+ gridRank = comm->getProcessID();
+
+ double dx = CoProcessor::grid->getDeltaX(level);
+
+ CoordinateTransformation3DPtr trafo = grid->getCoordinateTransformator();
+ double orgX1 = trafo->getX1CoordinateOffset();
+ double orgX2 = trafo->getX2CoordinateOffset();
+ double orgX3 = trafo->getX3CoordinateOffset();
+
+
+ int x1min = (line->getX1Minimum()-orgX1)/dx;
+ int x1max = (line->getX1Maximum()-orgX1)/dx;
+ int x2min = (line->getX2Minimum()-orgX2)/dx;
+ int x2max = (line->getX2Maximum()-orgX2)/dx;
+ int x3min = (line->getX3Minimum()-orgX3)/dx;
+ int x3max = (line->getX3Maximum()-orgX3)/dx;
+
+ UbTupleInt3 blockNx = grid->getBlockNX();
+
+ if (x1min!=x1max)
+ {
+ dir = X1;
+ blocknx = val<1>(blockNx);
+ length = x1max;
+ }
+ else if (x2min!=x2max)
+ {
+ dir = X2;
+ blocknx = val<2>(blockNx);
+ length = x2max;
+ }
+ else if (x3min!=x3max)
+ {
+ dir = X3;
+ blocknx = val<3>(blockNx);
+ length = x3max;
+ }
+ blockix1 = x1min/val<1>(blockNx);
+ blockix2 = x2min/val<2>(blockNx);
+ blockix3 = x3min/val<3>(blockNx);
+
+ ix1 = x1min%val<1>(blockNx)+1;
+ ix2 = x2min%val<2>(blockNx)+1;
+ ix3 = x3min%val<3>(blockNx)+1;
+}
+//////////////////////////////////////////////////////////////////////////
+void LineTimeSeriesCoProcessor::process(double step)
+{
+ if (scheduler->isDue(step))
+ {
+ collectData();
+ }
+
+ UBLOG(logDEBUG3, "MacroscopicQuantitiesCoProcessor::update:"<<step);
+}
+//////////////////////////////////////////////////////////////////////////
+void LineTimeSeriesCoProcessor::writeLine(const std::string& path)
+{
+ std::vector<UbTupleFloat3 > nodes(2);
+ std::vector<UbTupleInt2 > lines(1);
+ val<1>(nodes[0]) = line->getX1Minimum();
+ val<2>(nodes[0]) = line->getX2Minimum();
+ val<3>(nodes[0]) = line->getX3Minimum();
+ val<1>(nodes[1]) = line->getX1Maximum();
+ val<2>(nodes[1]) = line->getX2Maximum();
+ val<3>(nodes[1]) = line->getX3Maximum();
+ val<1>(lines[0]) = 0;
+ val<1>(lines[0]) = 1;
+ WbWriterVtkXmlASCII *writer = WbWriterVtkXmlASCII::getInstance();
+ writer->writeLines(path, nodes, lines);
+}
+//////////////////////////////////////////////////////////////////////////
+void LineTimeSeriesCoProcessor::collectData()
+{
+ LBMReal f[27];
+ LBMReal vx1, vx2, vx3, rho;
+ MPI_Status status;
+ std::vector<double> v1(length, 0);
+ std::vector<double> v2(length, 0);
+ std::vector<double> v3(length, 0);
+ std::vector<double> p(length, 0);
+ for (int x = 0; x<length; x += blocknx)
+ {
+ if (dir == X1)
+ {
+ blockix1 = x/blocknx;
+ }
+ else if (dir == X2)
+ {
+ blockix2 = x/blocknx;
+ }
+ else if (dir == X3)
+ {
+ blockix3 = x/blocknx;
+ }
+
+ Block3DPtr block = CoProcessor::grid->getBlock(blockix1, blockix2, blockix3, level);
+ if (block)
+ {
+ if (block->getRank()==gridRank)
+ {
+ LBMKernel3DPtr kernel = block->getKernel();
+ calcMacros = NULL;
+ if (kernel->getCompressible())
+ {
+ calcMacros = &D3Q27System::calcCompMacroscopicValues;
+ }
+ else
+ {
+ calcMacros = &D3Q27System::calcIncompMacroscopicValues;
+ }
+ DistributionArray3DPtr distributions = kernel->getDataSet()->getFdistributions();
+
+ for (int ix = 1; ix<=blocknx; ix++)
+ {
+ if (dir==X1)
+ {
+ ix1 = ix;
+ }
+ else if (dir==X2)
+ {
+ ix2 = ix;
+ }
+ else if (dir==X3)
+ {
+ ix3 = ix;
+ }
+ distributions->getDistribution(f, ix1, ix2, ix3);
+ calcMacros(f, rho, vx1, vx2, vx3);
+ v1[x+(ix-1)] = vx1;
+ v2[x+(ix-1)] = vx2;
+ v3[x+(ix-1)] = vx3;
+ p[x+(ix-1)] = rho;
+ }
+ }
+ }
+
+ }
+
+ if (root)
+ {
+ for (int i = 1; i<numOfProc; i++)
+ {
+ std::vector<double> v1temp(length, 0);
+ std::vector<double> v2temp(length, 0);
+ std::vector<double> v3temp(length, 0);
+ std::vector<double> ptemp(length, 0);
+ MPI_Recv(&v1temp[0], length, MPI_DOUBLE, i, 1, mpi_comm, &status);
+ MPI_Recv(&v2temp[0], length, MPI_DOUBLE, i, 2, mpi_comm, &status);
+ MPI_Recv(&v3temp[0], length, MPI_DOUBLE, i, 3, mpi_comm, &status);
+ MPI_Recv( &ptemp[0], length, MPI_DOUBLE, i, 4, mpi_comm, &status);
+ for (int j = 0; j<length; j++)
+ {
+ v1[j] += v1temp[j];
+ v2[j] += v2temp[j];
+ v3[j] += v3temp[j];
+ p[j] += ptemp[j];
+ }
+ }
+
+ std::ofstream ostr;
+ ostr.open(fname.c_str(), std::ios_base::out|std::ios_base::app);
+ if (!ostr)
+ {
+ ostr.clear();
+ std::string path = UbSystem::getPathFromString(fname);
+ if (path.size()>0) { UbSystem::makeDirectory(path); ostr.open(fname.c_str(), std::ios_base::out|std::ios_base::app); }
+ if (!ostr) throw UbException(UB_EXARGS, "couldn't open file "+fname);
+ }
+
+ for (int x = 0; x<length; x++)
+ {
+ ostr<<v1[x];
+ if (x<length-1)
+ {
+ ostr<<";";
+ }
+ }
+ ostr<<"\n";
+ for (int x = 0; x<length; x++)
+ {
+ ostr<<v2[x];
+ if (x<length-1)
+ {
+ ostr<<";";
+ }
+ }
+ ostr<<"\n";
+ for (int x = 0; x<length; x++)
+ {
+ ostr<<v3[x];
+ if (x<length-1)
+ {
+ ostr<<";";
+ }
+ }
+ ostr<<"\n";
+ for (int x = 0; x<length; x++)
+ {
+ ostr<<p[x];
+ if (x<length-1)
+ {
+ ostr<<";";
+ }
+ }
+ ostr<<"\n";
+ ostr.close();
+ }
+ else
+ {
+ MPI_Send(&v1[0], length, MPI_DOUBLE, 0, 1, mpi_comm);
+ MPI_Send(&v2[0], length, MPI_DOUBLE, 0, 2, mpi_comm);
+ MPI_Send(&v3[0], length, MPI_DOUBLE, 0, 3, mpi_comm);
+ MPI_Send( &p[0], length, MPI_DOUBLE, 0, 4, mpi_comm);
+ }
+}
+
+
diff --git a/source/VirtualFluidsCore/CoProcessors/LineTimeSeriesCoProcessor.h b/source/VirtualFluidsCore/CoProcessors/LineTimeSeriesCoProcessor.h
new file mode 100644
index 0000000000000000000000000000000000000000..9cedcfddea7ebbff079208db82b28769d51cbfcb
--- /dev/null
+++ b/source/VirtualFluidsCore/CoProcessors/LineTimeSeriesCoProcessor.h
@@ -0,0 +1,47 @@
+#ifndef LineTimeSeriesCoProcessor_h__
+#define LineTimeSeriesCoProcessor_h__
+
+#include "CoProcessor.h"
+#include "GbLine3D.h"
+#include "MPICommunicator.h"
+
+//! \brief Writes to .csv file time series for a line in x1 direction.
+//! \details It can be used to compute for given time range the time averaged two-point correlations for a line. <br>
+//! \f$ R_{ij}(x_{a},x{b},t) = <u_{i}(x_{a},t)u_{j}(x_{a}+r,t)> \f$ <br>
+//
+//! \author Konstantin Kutscher
+
+class LineTimeSeriesCoProcessor : public CoProcessor
+{
+public:
+enum Direction {X1, X2, X3};
+public:
+ LineTimeSeriesCoProcessor(Grid3DPtr grid, UbSchedulerPtr s, const std::string& path, GbLine3DPtr line, int level, CommunicatorPtr comm);
+ ~LineTimeSeriesCoProcessor(){}
+ void process(double step);
+ void writeLine(const std::string& path);
+protected:
+ void collectData();
+private:
+ std::string path;
+ std::string fname;
+ bool root;
+ GbLine3DPtr line;
+ //function pointer
+ typedef void(*CalcMacrosFct)(const LBMReal* const& /*feq[27]*/, LBMReal& /*(d)rho*/, LBMReal& /*vx1*/, LBMReal& /*vx2*/, LBMReal& /*vx3*/);
+ CalcMacrosFct calcMacros;
+ int blocknx;
+ int blockix1;
+ int blockix2;
+ int blockix3;
+ int level;
+ int ix1;
+ int ix2;
+ int ix3;
+ int length;
+ MPI_Comm mpi_comm;
+ int numOfProc;
+ int gridRank;
+ Direction dir;
+};
+#endif // LineTimeSeriesCoProcessor_h__
diff --git a/source/VirtualFluidsCore/CoProcessors/MacroscopicQuantitiesCoProcessor.h b/source/VirtualFluidsCore/CoProcessors/MacroscopicQuantitiesCoProcessor.h
index 2fa8c3a5fa57940b454560427f4bbdb4d020a09d..2f25a109a5d3f81e23897cedf804dd393cdee1ca 100644
--- a/source/VirtualFluidsCore/CoProcessors/MacroscopicQuantitiesCoProcessor.h
+++ b/source/VirtualFluidsCore/CoProcessors/MacroscopicQuantitiesCoProcessor.h
@@ -20,6 +20,7 @@ public:
const std::string& path, WbWriter* const writer,
LBMUnitConverterPtr conv,
bool nodesInformation = false);
+ ~MacroscopicQuantitiesCoProcessor(){}
void process(double step);
protected:
void collectData(double step);
diff --git a/source/VirtualFluidsCore/CoProcessors/TimeAveragedValuesCoProcessor.cpp b/source/VirtualFluidsCore/CoProcessors/TimeAveragedValuesCoProcessor.cpp
index 3304a2626e72503d4d34e2843774649e1f88be11..157fd55f0bd38c1b2b6ef1197d59148b0784eb32 100644
--- a/source/VirtualFluidsCore/CoProcessors/TimeAveragedValuesCoProcessor.cpp
+++ b/source/VirtualFluidsCore/CoProcessors/TimeAveragedValuesCoProcessor.cpp
@@ -105,7 +105,8 @@ void TimeAveragedValuesCoProcessor::init(UbSchedulerPtr s)
levelFactor = 1 << maxInitLevel;
maxStep = scheduler->getMaxEnd();
- numberOfFineSteps = (maxStep - scheduler->getMinBegin()) * levelFactor;
+ numberOfFineSteps = int(maxStep - scheduler->getMinBegin()) * levelFactor;
+ numberOfSteps = int(maxStep - scheduler->getMinBegin());
//function pointer
using namespace D3Q27System;
@@ -122,13 +123,16 @@ void TimeAveragedValuesCoProcessor::init(UbSchedulerPtr s)
//////////////////////////////////////////////////////////////////////////
void TimeAveragedValuesCoProcessor::process(double step)
{
+ calculateSubtotal(step);
+
if (step == maxStep)
{
//DEBUG/////////////////////
//UBLOG(logINFO, "process::step = " << step << ", maxStep = " << maxStep << ", levelFactor = " << levelFactor << ", numberOfFineSteps = " << numberOfFineSteps);
////////////////////////////
- calculateAverageValues((double)numberOfFineSteps);
+ //calculateAverageValues((double)numberOfFineSteps);
+ calculateAverageValues((double)numberOfSteps);
if (timeAveraging)
{
diff --git a/source/VirtualFluidsCore/CoProcessors/TimeAveragedValuesCoProcessor.h b/source/VirtualFluidsCore/CoProcessors/TimeAveragedValuesCoProcessor.h
index 861f96dc50dddba61eb4dc54c0b1a4dce7462b7a..939aaab07ec196bf5cd875a54cd6bd0a0f5a3d3a 100644
--- a/source/VirtualFluidsCore/CoProcessors/TimeAveragedValuesCoProcessor.h
+++ b/source/VirtualFluidsCore/CoProcessors/TimeAveragedValuesCoProcessor.h
@@ -20,10 +20,6 @@ typedef boost::shared_ptr<TimeAveragedValuesCoProcessor> TimeAveragedValuesCoPro
//! \author Konstantin Kutscher
// \f$ u_{mean}=\frac{1}{N}\sum\limits_{i=1}^n u_{i} \f$
-//struct plotZ
-//{
-//
-//};
class TimeAveragedValuesCoProcessor : public CoProcessor
{
@@ -89,6 +85,7 @@ private:
int options;
int lcounter;
int numberOfFineSteps;
+ int numberOfSteps;
int fineStep;
int minFineStep;
int maxFineStep;
diff --git a/source/VirtualFluidsCore/Data/D3Q27EsoTwist3DSplittedVector.cpp b/source/VirtualFluidsCore/Data/D3Q27EsoTwist3DSplittedVector.cpp
index 68173279730e22414810bfa7354980077ad8ca41..e8a1185fa59f1d98f9267b39d6cedb10d2562f96 100644
--- a/source/VirtualFluidsCore/Data/D3Q27EsoTwist3DSplittedVector.cpp
+++ b/source/VirtualFluidsCore/Data/D3Q27EsoTwist3DSplittedVector.cpp
@@ -62,15 +62,15 @@ void D3Q27EsoTwist3DSplittedVector::getDistribution(LBMReal* const f, size_t x1,
//////////////////////////////////////////////////////////////////////////
void D3Q27EsoTwist3DSplittedVector::setDistribution(const LBMReal* const f, size_t x1, size_t x2, size_t x3)
{
- (*this->localDistributions)(D3Q27System::ET_E,x1, x2, x3) = f[D3Q27System::INV_E];
- (*this->localDistributions)(D3Q27System::ET_N,x1, x2, x3) = f[D3Q27System::INV_N];
- (*this->localDistributions)(D3Q27System::ET_T,x1, x2, x3) = f[D3Q27System::INV_T];
- (*this->localDistributions)(D3Q27System::ET_NE,x1, x2, x3) = f[D3Q27System::INV_NE];
- (*this->localDistributions)(D3Q27System::ET_NW,x1+1,x2, x3) = f[D3Q27System::INV_NW];
- (*this->localDistributions)(D3Q27System::ET_TE,x1, x2, x3) = f[D3Q27System::INV_TE];
- (*this->localDistributions)(D3Q27System::ET_TW,x1+1,x2, x3) = f[D3Q27System::INV_TW];
- (*this->localDistributions)(D3Q27System::ET_TN,x1, x2, x3) = f[D3Q27System::INV_TN];
- (*this->localDistributions)(D3Q27System::ET_TS,x1, x2+1,x3) = f[D3Q27System::INV_TS];
+ (*this->localDistributions)(D3Q27System::ET_E,x1, x2, x3) = f[D3Q27System::INV_E];
+ (*this->localDistributions)(D3Q27System::ET_N,x1, x2, x3) = f[D3Q27System::INV_N];
+ (*this->localDistributions)(D3Q27System::ET_T,x1, x2, x3) = f[D3Q27System::INV_T];
+ (*this->localDistributions)(D3Q27System::ET_NE,x1, x2, x3) = f[D3Q27System::INV_NE];
+ (*this->localDistributions)(D3Q27System::ET_NW,x1+1,x2, x3) = f[D3Q27System::INV_NW];
+ (*this->localDistributions)(D3Q27System::ET_TE,x1, x2, x3) = f[D3Q27System::INV_TE];
+ (*this->localDistributions)(D3Q27System::ET_TW,x1+1,x2, x3) = f[D3Q27System::INV_TW];
+ (*this->localDistributions)(D3Q27System::ET_TN,x1, x2, x3) = f[D3Q27System::INV_TN];
+ (*this->localDistributions)(D3Q27System::ET_TS,x1, x2+1,x3) = f[D3Q27System::INV_TS];
(*this->localDistributions)(D3Q27System::ET_TNE,x1, x2, x3) = f[D3Q27System::INV_TNE];
(*this->localDistributions)(D3Q27System::ET_TNW,x1+1,x2, x3) = f[D3Q27System::INV_TNW];
(*this->localDistributions)(D3Q27System::ET_TSE,x1, x2+1,x3) = f[D3Q27System::INV_TSE];
diff --git a/source/VirtualFluidsCore/Data/DataSet3D.h b/source/VirtualFluidsCore/Data/DataSet3D.h
index 1cc45b619a9964126f3b84fa127f114a8696a24e..888ac68dba4ae89936105f608444762ad4ad7dea 100644
--- a/source/VirtualFluidsCore/Data/DataSet3D.h
+++ b/source/VirtualFluidsCore/Data/DataSet3D.h
@@ -13,20 +13,11 @@ typedef boost::shared_ptr<DataSet3D> DataSet3DPtr;
typedef CbArray4D<LBMReal,IndexerX4X3X2X1> AverageValuesArray3D;
typedef boost::shared_ptr< AverageValuesArray3D > AverageValuesArray3DPtr;
-//typedef CbArray4D<LBMReal, IndexerX4X3X2X1> AverageVelocityArray3D;
-//typedef boost::shared_ptr< AverageValuesArray3D > AverageVelocityArray3DPtr;
-//
-//typedef CbArray4D<LBMReal, IndexerX4X3X2X1> AverageFluctuationsArray3D;
-//typedef boost::shared_ptr< AverageFluctuationsArray3D > AverageFluctuationsArray3DPtr;
-//
-//typedef CbArray4D<LBMReal, IndexerX4X3X2X1> AverageTriplecorrelationsArray3D;
-//typedef boost::shared_ptr< AverageValuesArray3D > AverageTriplecorrelationsArray3DPtr;
-
typedef CbArray4D<LBMReal,IndexerX4X3X2X1> ShearStressValuesArray3D;
typedef boost::shared_ptr< ShearStressValuesArray3D > ShearStressValuesArray3DPtr;
-typedef CbArray3D<LBMReal, IndexerX3X2X1> ViscosityArray3D;
-typedef boost::shared_ptr< ViscosityArray3D > ViscosityArray3DPtr;
+typedef CbArray3D<LBMReal, IndexerX3X2X1> RelaxationFactorArray3D;
+typedef boost::shared_ptr< RelaxationFactorArray3D > RelaxationFactorArray3DPtr;
class DataSet3D
{
@@ -42,96 +33,109 @@ public:
AverageValuesArray3DPtr getAverageTriplecorrelations() const;
void setAverageTriplecorrelations(AverageValuesArray3DPtr values);
-
-
-
AverageValuesArray3DPtr getAverageValues() const;
void setAverageValues(AverageValuesArray3DPtr values);
ShearStressValuesArray3DPtr getShearStressValues() const;
void setShearStressValues(ShearStressValuesArray3DPtr values);
+
+ RelaxationFactorArray3DPtr getRelaxationFactor() const;
+ void setRelaxationFactor(RelaxationFactorArray3DPtr values);
protected:
private:
- DistributionArray3DPtr mFdistributions;
- AverageValuesArray3DPtr mAverageValues;
+ DistributionArray3DPtr fdistributions;
+ AverageValuesArray3DPtr averageValues;
+
+ AverageValuesArray3DPtr averageVelocity;
+ AverageValuesArray3DPtr averageFluktuations;
+ AverageValuesArray3DPtr averageTriplecorrelations;
- AverageValuesArray3DPtr mAverageVelocity;
- AverageValuesArray3DPtr mAverageFluktuations;
- AverageValuesArray3DPtr mAverageTriplecorrelations;
+ ShearStressValuesArray3DPtr shearStressValues;
- ShearStressValuesArray3DPtr mShearStressValues;
+ RelaxationFactorArray3DPtr relaxationFactor;
friend class boost::serialization::access;
template<class Archive>
void serialize(Archive & ar, const unsigned int version)
{
- ar & mFdistributions;
- ar & mAverageValues;
- ar & mShearStressValues;
- ar & mAverageVelocity;
- ar & mAverageFluktuations;
- ar & mAverageTriplecorrelations;
+ ar & fdistributions;
+ ar & averageValues;
+ ar & shearStressValues;
+ ar & averageVelocity;
+ ar & averageFluktuations;
+ ar & averageTriplecorrelations;
+ ar & relaxationFactor;
}
};
inline DistributionArray3DPtr DataSet3D::getFdistributions() const
{
- return mFdistributions;
+ return fdistributions;
}
inline void DataSet3D::setFdistributions(DistributionArray3DPtr distributions)
{
- mFdistributions = distributions;
+ fdistributions = distributions;
}
inline AverageValuesArray3DPtr DataSet3D::getAverageValues() const
{
- return mAverageValues;
+ return averageValues;
}
inline void DataSet3D::setAverageValues(AverageValuesArray3DPtr values)
{
- mAverageValues = values;
+ averageValues = values;
}
inline AverageValuesArray3DPtr DataSet3D::getAverageVelocity() const
{
- return mAverageVelocity;
+ return averageVelocity;
}
inline void DataSet3D::setAverageVelocity(AverageValuesArray3DPtr values)
{
- mAverageVelocity = values;
+ averageVelocity = values;
}
inline AverageValuesArray3DPtr DataSet3D::getAverageFluctuations() const
{
- return mAverageFluktuations;
+ return averageFluktuations;
}
inline void DataSet3D::setAverageFluctuations(AverageValuesArray3DPtr values)
{
- mAverageFluktuations = values;
+ averageFluktuations = values;
}
inline AverageValuesArray3DPtr DataSet3D::getAverageTriplecorrelations() const
{
- return mAverageTriplecorrelations;
+ return averageTriplecorrelations;
}
inline void DataSet3D::setAverageTriplecorrelations(AverageValuesArray3DPtr values)
{
- mAverageTriplecorrelations = values;
+ averageTriplecorrelations = values;
}
inline ShearStressValuesArray3DPtr DataSet3D::getShearStressValues() const
{
- return mShearStressValues;
+ return shearStressValues;
}
inline void DataSet3D::setShearStressValues(ShearStressValuesArray3DPtr values)
{
- mShearStressValues = values;
+ shearStressValues = values;
+}
+
+inline RelaxationFactorArray3DPtr DataSet3D::getRelaxationFactor() const
+{
+ return relaxationFactor;
+}
+
+inline void DataSet3D::setRelaxationFactor(RelaxationFactorArray3DPtr values)
+{
+ relaxationFactor = values;
}
#endif
diff --git a/source/VirtualFluidsCore/Interactors/D3Q27Interactor.cpp b/source/VirtualFluidsCore/Interactors/D3Q27Interactor.cpp
index b56bddbb3d3894b9f3f2bdf8eb30ac6b8724b806..ae529f9fa7062e37e6b4f9e71801a5e589d8604c 100644
--- a/source/VirtualFluidsCore/Interactors/D3Q27Interactor.cpp
+++ b/source/VirtualFluidsCore/Interactors/D3Q27Interactor.cpp
@@ -652,6 +652,7 @@ vector< pair<GbPoint3D,GbPoint3D> > D3Q27Interactor::getQsLineSet()
{
LBMKernel3DPtr kernel = block->getKernel();
BCArray3D<D3Q27BoundaryCondition>& bcMatrix = boost::dynamic_pointer_cast<D3Q27ETBCProcessor>(kernel->getBCProcessor())->getBCArray();
+ UbTupleDouble3 nodeOffset = grid.lock()->getNodeOffset(block);
//double collFactor = ((LbD3Q27Calculator*)grid->getCalculator())->getCollisionsFactors()[block->getLevel()];
//checken ob obere reihe doppelt im system vorhanden oder nicht
@@ -701,9 +702,9 @@ vector< pair<GbPoint3D,GbPoint3D> > D3Q27Interactor::getQsLineSet()
{
if( !bcMatrix.hasBC(ix1,ix2,ix3) ) continue;
D3Q27BoundaryConditionPtr bc = bcMatrix.getBC(ix1,ix2,ix3);
- double x1a = x1+dx * ix1;
- double x2a = x2+dx * ix2;
- double x3a = x3+dx * ix3;
+ double x1a = x1-val<1>(nodeOffset)+dx * ix1;
+ double x2a = x2-val<2>(nodeOffset)+dx * ix2;
+ double x3a = x3-val<3>(nodeOffset)+dx * ix3;
pointpair.first.setX1(x1a);
pointpair.first.setX2(x2a);
pointpair.first.setX3(x3a);
diff --git a/source/VirtualFluidsCore/LBM/D3Q27IntegrateValuesHelper.cpp b/source/VirtualFluidsCore/LBM/D3Q27IntegrateValuesHelper.cpp
index 62929fd669f29bf55465f777c6460eda6f817aee..fd1a3ef209a09551c53a32aaaa11d0f93b66adcd 100644
--- a/source/VirtualFluidsCore/LBM/D3Q27IntegrateValuesHelper.cpp
+++ b/source/VirtualFluidsCore/LBM/D3Q27IntegrateValuesHelper.cpp
@@ -2,6 +2,7 @@
#include <boost/foreach.hpp>
#include <numerics/geometry3d/GbCuboid3D.h>
+#include <numerics/geometry3d/CoordinateTransformation3D.h>
#include <vector>
#include "LBMKernelETD3Q27.h"
@@ -135,6 +136,104 @@ void D3Q27IntegrateValuesHelper::init(int level)
}
}
+}
+//////////////////////////////////////////////////////////////////////////
+void D3Q27IntegrateValuesHelper::prepare2DMatrix(int level)
+{
+ root = comm->isRoot();
+
+ double orgX1, orgX2, orgX3;
+ int gridRank = grid->getRank();
+ int minInitLevel, maxInitLevel;
+ if (level<0)
+ {
+ minInitLevel = this->grid->getCoarsestInitializedLevel();
+ maxInitLevel = this->grid->getFinestInitializedLevel();
+ }
+ else
+ {
+ minInitLevel = level;
+ maxInitLevel = level;
+ }
+ double dx = grid->getDeltaX(level);
+ CoordinateTransformation3D trafo(boundingBox->getX1Minimum(),boundingBox->getX2Minimum(),boundingBox->getX3Minimum(),dx,dx,dx);
+ cnodes2DMatrix.resize(UbMath::integerRounding<double>(boundingBox->getX1Maximum()/dx),UbMath::integerRounding<double>(boundingBox->getX2Maximum()/dx));
+
+ double numSolids = 0.0;
+ double numFluids = 0.0;
+ for (int level = minInitLevel; level<=maxInitLevel; level++)
+ {
+ vector<Block3DPtr> blockVector;
+ grid->getBlocks(level, gridRank, blockVector);
+ BOOST_FOREACH(Block3DPtr block, blockVector)
+ {
+ Node cn;
+ cn.block = block;
+ //Koords bestimmen
+ UbTupleDouble3 org = grid->getBlockWorldCoordinates(block);
+
+ orgX1 = val<1>(org);
+ orgX2 = val<2>(org);
+ orgX3 = val<3>(org);
+
+ LBMKernelETD3Q27Ptr kernel = boost::dynamic_pointer_cast<LBMKernelETD3Q27>(block->getKernel());
+ BCArray3D<D3Q27BoundaryCondition>& bcArray = boost::dynamic_pointer_cast<D3Q27ETBCProcessor>(kernel->getBCProcessor())->getBCArray();
+ int ghostLayerWitdh = kernel->getGhostLayerWidth();
+ DistributionArray3DPtr distributions = kernel->getDataSet()->getFdistributions();
+ double internX1, internX2, internX3;
+
+ double dx = grid->getDeltaX(block);
+ UbTupleDouble3 orgDelta = grid->getNodeOffset(block);
+
+ for (int ix3 = ghostLayerWitdh; ix3<(int)distributions->getNX3()-ghostLayerWitdh; ix3++)
+ {
+ for (int ix2 = ghostLayerWitdh; ix2<(int)distributions->getNX2()-ghostLayerWitdh; ix2++)
+ {
+ for (int ix1 = ghostLayerWitdh; ix1<(int)distributions->getNX1()-ghostLayerWitdh; ix1++)
+ {
+ internX1 = orgX1-val<1>(orgDelta)+ix1 * dx;
+ internX2 = orgX2-val<2>(orgDelta)+ix2 * dx;
+ internX3 = orgX3-val<3>(orgDelta)+ix3 * dx;
+ if (boundingBox->isPointInGbObject3D(internX1, internX2, internX3))
+ {
+ if (!bcArray.isSolid(ix1, ix2, ix3)&&!bcArray.isUndefined(ix1, ix2, ix3))
+ {
+ cn.node = UbTupleInt3(ix1, ix2, ix3);
+ int x1 = (int)trafo.transformForwardToX1Coordinate(internX1, internX2, internX3);
+ int x2 = (int)trafo.transformForwardToX2Coordinate(internX1, internX2, internX3);
+ cnodes2DMatrix(x1,x2)=cn;
+ numFluids++;
+ }
+ else if (bcArray.isSolid(ix1, ix2, ix3))
+ {
+ numSolids++;
+ }
+ }
+ }
+ }
+ }
+
+ }
+ }
+ vector<double> rvalues;
+ vector<double> values;
+ values.push_back(numSolids);
+ values.push_back(numFluids);
+ rvalues = comm->gather(values);
+
+ if (root)
+ {
+ numberOfSolidNodes = 0.0;
+ numberOfFluidsNodes = 0.0;
+ int rsize = (int)rvalues.size();
+ int vsize = (int)values.size();
+ for (int i = 0; i<rsize; i += vsize)
+ {
+ numberOfSolidNodes += rvalues[i];
+ numberOfFluidsNodes += rvalues[i+1];
+ }
+ }
+
}
// calculation conventional rho, velocity and averaged data
void D3Q27IntegrateValuesHelper::calculateAV()
@@ -354,6 +453,45 @@ void D3Q27IntegrateValuesHelper::calculateAV2()
}
}
}
+//////////////////////////////////////////////////////////////////////////
+void D3Q27IntegrateValuesHelper::calculateTwoPointCorrelations()
+{
+ int x1max = cnodes2DMatrix.getNX1();
+ int x2max = cnodes2DMatrix.getNX2();
+ LBMReal f[27];
+
+ for (int x2; x2 < x2max; x2++)
+ {
+ for (int x1; x1 < x1max; x1++)
+ {
+ Block3DPtr block = cnodes2DMatrix(x1,x2).block;
+ UbTupleInt3 node = cnodes2DMatrix(x1,x2).node;
+ LBMKernel3DPtr kernel = block->getKernel();
+ DistributionArray3DPtr distributions = kernel->getDataSet()->getFdistributions();
+ distributions->getDistribution(f, val<1>(node), val<2>(node), val<3>(node));
+
+ //Funktionszeiger
+ typedef void(*CalcMacrosFct)(const LBMReal* const& /*feq[27]*/, LBMReal& /*(d)rho*/, LBMReal& /*vx1*/, LBMReal& /*vx2*/, LBMReal& /*vx3*/);
+
+ CalcMacrosFct calcMacros = NULL;
+ if (kernel->getCompressible())
+ {
+ calcMacros = &D3Q27System::calcCompMacroscopicValues;
+ }
+ else
+ {
+ calcMacros = &D3Q27System::calcIncompMacroscopicValues;
+ }
+
+ //////////////////////////////////////////////////////////////////////////
+ //compute velocity
+ //////////////////////////////////////////////////////////////////////////
+ LBMReal vx, vy, vz, rho;
+ calcMacros(f, rho, vx, vy, vz);
+
+ }
+ }
+}
//////////////////////////////////////////////////////////////////////////
void D3Q27IntegrateValuesHelper::calculateMQ()
@@ -455,7 +593,7 @@ GbCuboid3DPtr D3Q27IntegrateValuesHelper::getBoundingBox()
return this->boundingBox;
}
//////////////////////////////////////////////////////////////////////////
-std::vector<CalcNodes> D3Q27IntegrateValuesHelper::getCNodes()
+std::vector<D3Q27IntegrateValuesHelper::CalcNodes> D3Q27IntegrateValuesHelper::getCNodes()
{
return cnodes;
}
diff --git a/source/VirtualFluidsCore/LBM/D3Q27IntegrateValuesHelper.h b/source/VirtualFluidsCore/LBM/D3Q27IntegrateValuesHelper.h
index 034bf5da978e2b054d9df5771cee5cee322fcd70..f223b07a04558dd5b8b0063d16bcfd7868ed9971 100644
--- a/source/VirtualFluidsCore/LBM/D3Q27IntegrateValuesHelper.h
+++ b/source/VirtualFluidsCore/LBM/D3Q27IntegrateValuesHelper.h
@@ -5,18 +5,37 @@
#include "D3Q27System.h"
#include "Communicator.h"
#include "GbCuboid3D.h"
+#include "CbArray2D.h"
-struct CalcNodes
-{
- Block3DPtr block;
- std::vector<UbTupleInt3> nodes;
-};
+//struct CalcNodes
+//{
+// Block3DPtr block;
+// std::vector<UbTupleInt3> nodes;
+//};
+//
+//struct Nodes
+//{
+// Block3DPtr block;
+// UbTupleInt3 nodes;
+//};
class D3Q27IntegrateValuesHelper;
typedef boost::shared_ptr<D3Q27IntegrateValuesHelper> D3Q27IntegrateValuesHelperPtr;
class D3Q27IntegrateValuesHelper
{
+public:
+ struct CalcNodes
+ {
+ Block3DPtr block;
+ std::vector<UbTupleInt3> nodes;
+ };
+
+ struct Node
+ {
+ Block3DPtr block;
+ UbTupleInt3 node;
+ };
public:
D3Q27IntegrateValuesHelper(Grid3DPtr grid, CommunicatorPtr comm,
double minX1, double minX2, double minX3,
@@ -29,6 +48,8 @@ public:
void calculateMQ();
void calculateAV();
void calculateAV2();
+ void prepare2DMatrix(int level);
+ void calculateTwoPointCorrelations();
void clearData();
double getRho() {return sRho;}
@@ -86,6 +107,7 @@ private:
std::vector<CalcNodes> cnodes;
GbCuboid3DPtr boundingBox;
CommunicatorPtr comm;
+ CbArray2D<Node> cnodes2DMatrix;
enum Values{AvVx = 0, AvVy = 1, AvVz = 2, AvVxx = 3, AvVyy = 4, AvVzz = 5, AvVxy = 6, AvVyz = 7, AvVxz = 8};
double saVx, saVy, saVz;
diff --git a/source/VirtualFluidsCore/LBM/InitDensityLBMKernel.cpp b/source/VirtualFluidsCore/LBM/InitDensityLBMKernel.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..f06861bff40e200794a756fa4750d4e369912165
--- /dev/null
+++ b/source/VirtualFluidsCore/LBM/InitDensityLBMKernel.cpp
@@ -0,0 +1,1074 @@
+#include "InitDensityLBMKernel.h"
+#include "D3Q27EsoTwist3DSplittedVector.h"
+
+InitDensityLBMKernel::InitDensityLBMKernel()
+{
+}
+
+InitDensityLBMKernel::~InitDensityLBMKernel()
+{
+}
+
+InitDensityLBMKernel::InitDensityLBMKernel(int nx1, int nx2, int nx3) : LBMKernelETD3Q27(nx1, nx2, nx3)
+{
+ this->compressible = false;
+}
+
+void InitDensityLBMKernel::calculate()
+{
+ collideAll();
+}
+
+LBMKernel3DPtr InitDensityLBMKernel::clone()
+{
+ LBMKernel3DPtr kernel(new InitDensityLBMKernel(nx1, nx2, nx3));
+ boost::dynamic_pointer_cast<InitDensityLBMKernel>(kernel)->init();
+ kernel->setCollisionFactor(this->collFactor);
+ kernel->setBCProcessor(bcProcessor->clone(kernel));
+ kernel->setWithForcing(withForcing);
+ kernel->setForcingX1(muForcingX1);
+ kernel->setForcingX2(muForcingX2);
+ kernel->setForcingX3(muForcingX3);
+ kernel->setIndex(ix1, ix2, ix3);
+ kernel->setDeltaT(deltaT);
+ boost::dynamic_pointer_cast<InitDensityLBMKernel>(kernel)->OxyyMxzz = 1.0;
+ return kernel;
+}
+
+void InitDensityLBMKernel::setVelocity(int x1, int x2, int x3, LBMReal vvx, LBMReal vvy, LBMReal vvz)
+{
+ v(0, x1, x2, x3) = vvx;
+ v(1, x1, x2, x3) = vvy;
+ v(2, x1, x2, x3) = vvz;
+}
+
+double InitDensityLBMKernel::getCallculationTime()
+{
+ return 0;
+}
+
+//void InitDensityLBMKernel::collideAll()
+//{
+// using namespace D3Q27System;
+//
+// localDistributions = boost::dynamic_pointer_cast<D3Q27EsoTwist3DSplittedVector>(dataSet->getFdistributions())->getLocalDistributions();
+// nonLocalDistributions = boost::dynamic_pointer_cast<D3Q27EsoTwist3DSplittedVector>(dataSet->getFdistributions())->getNonLocalDistributions();
+// zeroDistributions = boost::dynamic_pointer_cast<D3Q27EsoTwist3DSplittedVector>(dataSet->getFdistributions())->getZeroDistributions();
+//
+// BCArray3D<D3Q27BoundaryCondition>& bcArray = boost::dynamic_pointer_cast<D3Q27ETBCProcessor>(this->getBCProcessor())->getBCArray();
+//
+// const int bcArrayMaxX1 = (int)bcArray.getNX1();
+// const int bcArrayMaxX2 = (int)bcArray.getNX2();
+// const int bcArrayMaxX3 = (int)bcArray.getNX3();
+//
+// int minX1 = ghostLayerWidth;
+// int minX2 = ghostLayerWidth;
+// int minX3 = ghostLayerWidth;
+// int maxX1 = bcArrayMaxX1-ghostLayerWidth;
+// int maxX2 = bcArrayMaxX2-ghostLayerWidth;
+// int maxX3 = bcArrayMaxX3-ghostLayerWidth;
+//
+//
+// for (int x3 = minX3; x3<maxX3; x3++)
+// {
+// for (int x2 = minX2; x2<maxX2; x2++)
+// {
+// for (int x1 = minX1; x1<maxX1; x1++)
+// {
+// if (!bcArray.isSolid(x1, x2, x3)&&!bcArray.isUndefined(x1, x2, x3))
+// {
+// int x1p = x1+1;
+// int x2p = x2+1;
+// int x3p = x3+1;
+// //////////////////////////////////////////////////////////////////////////
+// //read distribution
+// ////////////////////////////////////////////////////////////////////////////
+// //////////////////////////////////////////////////////////////////////////
+//
+// //E N T
+// //c c c
+// //////////
+// //W S B
+// //a a a
+//
+// //Rest ist b
+//
+// //mfxyz
+// //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 m0, m1, m2;
+//
+// LBMReal rho = (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;
+//
+// //LBMReal vvx = ((((mfccc-mfaaa)+(mfcac-mfaca))+((mfcaa-mfacc)+(mfcca-mfaac)))+
+// // (((mfcba-mfabc)+(mfcbc-mfaba))+((mfcab-mfacb)+(mfccb-mfaab)))+
+// // (mfcbb-mfabb));
+// //LBMReal vvy = ((((mfccc-mfaaa)+(mfaca-mfcac))+((mfacc-mfcaa)+(mfcca-mfaac)))+
+// // (((mfbca-mfbac)+(mfbcc-mfbaa))+((mfacb-mfcab)+(mfccb-mfaab)))+
+// // (mfbcb-mfbab));
+// //LBMReal vvz = ((((mfccc-mfaaa)+(mfcac-mfaca))+((mfacc-mfcaa)+(mfaac-mfcca)))+
+// // (((mfbac-mfbca)+(mfbcc-mfbaa))+((mfabc-mfcba)+(mfcbc-mfaba)))+
+// // (mfbbc-mfbba));
+//
+// LBMReal vvx = v(0,x1,x2,x3);
+// LBMReal vvy = v(1,x1,x2,x3);
+// LBMReal vvz = v(2,x1,x2,x3);
+// //LBMReal rho = v(3,x1,x2,x3);
+//
+// LBMReal oMdrho;
+//
+// oMdrho = mfccc+mfaaa;
+// m0 = mfaca+mfcac;
+// m1 = mfacc+mfcaa;
+// m2 = mfaac+mfcca;
+// oMdrho += m0;
+// m1 += m2;
+// oMdrho += m1;
+// m0 = mfbac+mfbca;
+// m1 = mfbaa+mfbcc;
+// m0 += m1;
+// m1 = mfabc+mfcba;
+// m2 = mfaba+mfcbc;
+// m1 += m2;
+// m0 += m1;
+// m1 = mfacb+mfcab;
+// m2 = mfaab+mfccb;
+// m1 += m2;
+// m0 += m1;
+// oMdrho += m0;
+// m0 = mfabb+mfcbb;
+// m1 = mfbab+mfbcb;
+// m2 = mfbba+mfbbc;
+// m0 += m1+m2;
+// m0 += mfbbb; //hat gefehlt
+// oMdrho = 1.-(oMdrho+m0);
+//
+// LBMReal vx2;
+// LBMReal vy2;
+// LBMReal vz2;
+// vx2 = vvx*vvx;
+// vy2 = vvy*vvy;
+// vz2 = vvz*vvz;
+// ////////////////////////////////////////////////////////////////////////////////////
+// LBMReal wadjust;
+// LBMReal qudricLimit = 0.01;
+// ////////////////////////////////////////////////////////////////////////////////////
+// //Hin
+// ////////////////////////////////////////////////////////////////////////////////////
+// // mit 1/36, 1/9, 1/36, 1/9, 4/9, 1/9, 1/36, 1/9, 1/36 Konditionieren
+// ////////////////////////////////////////////////////////////////////////////////////
+// // Z - Dir
+// m2 = mfaaa+mfaac;
+// m1 = mfaac-mfaaa;
+// m0 = m2+mfaab;
+// mfaaa = m0;
+// m0 += c1o36 * oMdrho;
+// mfaab = m1-m0 * vvz;
+// mfaac = m2-2. * m1 * vvz+vz2 * m0;
+// ////////////////////////////////////////////////////////////////////////////////////
+// m2 = mfaba+mfabc;
+// m1 = mfabc-mfaba;
+// m0 = m2+mfabb;
+// mfaba = m0;
+// m0 += c1o9 * oMdrho;
+// mfabb = m1-m0 * vvz;
+// mfabc = m2-2. * m1 * vvz+vz2 * m0;
+// ////////////////////////////////////////////////////////////////////////////////////
+// m2 = mfaca+mfacc;
+// m1 = mfacc-mfaca;
+// m0 = m2+mfacb;
+// mfaca = m0;
+// m0 += c1o36 * oMdrho;
+// mfacb = m1-m0 * vvz;
+// mfacc = m2-2. * m1 * vvz+vz2 * m0;
+// ////////////////////////////////////////////////////////////////////////////////////
+// ////////////////////////////////////////////////////////////////////////////////////
+// m2 = mfbaa+mfbac;
+// m1 = mfbac-mfbaa;
+// m0 = m2+mfbab;
+// mfbaa = m0;
+// m0 += c1o9 * oMdrho;
+// mfbab = m1-m0 * vvz;
+// mfbac = m2-2. * m1 * vvz+vz2 * m0;
+// ////////////////////////////////////////////////////////////////////////////////////
+// m2 = mfbba+mfbbc;
+// m1 = mfbbc-mfbba;
+// m0 = m2+mfbbb;
+// mfbba = m0;
+// m0 += c4o9 * oMdrho;
+// mfbbb = m1-m0 * vvz;
+// mfbbc = m2-2. * m1 * vvz+vz2 * m0;
+// ////////////////////////////////////////////////////////////////////////////////////
+// m2 = mfbca+mfbcc;
+// m1 = mfbcc-mfbca;
+// m0 = m2+mfbcb;
+// mfbca = m0;
+// m0 += c1o9 * oMdrho;
+// mfbcb = m1-m0 * vvz;
+// mfbcc = m2-2. * m1 * vvz+vz2 * m0;
+// ////////////////////////////////////////////////////////////////////////////////////
+// ////////////////////////////////////////////////////////////////////////////////////
+// m2 = mfcaa+mfcac;
+// m1 = mfcac-mfcaa;
+// m0 = m2+mfcab;
+// mfcaa = m0;
+// m0 += c1o36 * oMdrho;
+// mfcab = m1-m0 * vvz;
+// mfcac = m2-2. * m1 * vvz+vz2 * m0;
+// ////////////////////////////////////////////////////////////////////////////////////
+// m2 = mfcba+mfcbc;
+// m1 = mfcbc-mfcba;
+// m0 = m2+mfcbb;
+// mfcba = m0;
+// m0 += c1o9 * oMdrho;
+// mfcbb = m1-m0 * vvz;
+// mfcbc = m2-2. * m1 * vvz+vz2 * m0;
+// ////////////////////////////////////////////////////////////////////////////////////
+// m2 = mfcca+mfccc;
+// m1 = mfccc-mfcca;
+// m0 = m2+mfccb;
+// mfcca = m0;
+// m0 += c1o36 * oMdrho;
+// mfccb = m1-m0 * vvz;
+// mfccc = m2-2. * m1 * vvz+vz2 * m0;
+// ////////////////////////////////////////////////////////////////////////////////////
+// ////////////////////////////////////////////////////////////////////////////////////
+// // mit 1/6, 0, 1/18, 2/3, 0, 2/9, 1/6, 0, 1/18 Konditionieren
+// ////////////////////////////////////////////////////////////////////////////////////
+// // Y - Dir
+// m2 = mfaaa+mfaca;
+// m1 = mfaca-mfaaa;
+// m0 = m2+mfaba;
+// mfaaa = m0;
+// m0 += c1o6 * oMdrho;
+// mfaba = m1-m0 * vvy;
+// mfaca = m2-2. * m1 * vvy+vy2 * m0;
+// ////////////////////////////////////////////////////////////////////////////////////
+// m2 = mfaab+mfacb;
+// m1 = mfacb-mfaab;
+// m0 = m2+mfabb;
+// mfaab = m0;
+// mfabb = m1-m0 * vvy;
+// mfacb = m2-2. * m1 * vvy+vy2 * m0;
+// ////////////////////////////////////////////////////////////////////////////////////
+// m2 = mfaac+mfacc;
+// m1 = mfacc-mfaac;
+// m0 = m2+mfabc;
+// mfaac = m0;
+// m0 += c1o18 * oMdrho;
+// mfabc = m1-m0 * vvy;
+// mfacc = m2-2. * m1 * vvy+vy2 * m0;
+// ////////////////////////////////////////////////////////////////////////////////////
+// ////////////////////////////////////////////////////////////////////////////////////
+// m2 = mfbaa+mfbca;
+// m1 = mfbca-mfbaa;
+// m0 = m2+mfbba;
+// mfbaa = m0;
+// m0 += c2o3 * oMdrho;
+// mfbba = m1-m0 * vvy;
+// mfbca = m2-2. * m1 * vvy+vy2 * m0;
+// ////////////////////////////////////////////////////////////////////////////////////
+// m2 = mfbab+mfbcb;
+// m1 = mfbcb-mfbab;
+// m0 = m2+mfbbb;
+// mfbab = m0;
+// mfbbb = m1-m0 * vvy;
+// mfbcb = m2-2. * m1 * vvy+vy2 * m0;
+// ////////////////////////////////////////////////////////////////////////////////////
+// m2 = mfbac+mfbcc;
+// m1 = mfbcc-mfbac;
+// m0 = m2+mfbbc;
+// mfbac = m0;
+// m0 += c2o9 * oMdrho;
+// mfbbc = m1-m0 * vvy;
+// mfbcc = m2-2. * m1 * vvy+vy2 * m0;
+// ////////////////////////////////////////////////////////////////////////////////////
+// ////////////////////////////////////////////////////////////////////////////////////
+// m2 = mfcaa+mfcca;
+// m1 = mfcca-mfcaa;
+// m0 = m2+mfcba;
+// mfcaa = m0;
+// m0 += c1o6 * oMdrho;
+// mfcba = m1-m0 * vvy;
+// mfcca = m2-2. * m1 * vvy+vy2 * m0;
+// ////////////////////////////////////////////////////////////////////////////////////
+// m2 = mfcab+mfccb;
+// m1 = mfccb-mfcab;
+// m0 = m2+mfcbb;
+// mfcab = m0;
+// mfcbb = m1-m0 * vvy;
+// mfccb = m2-2. * m1 * vvy+vy2 * m0;
+// ////////////////////////////////////////////////////////////////////////////////////
+// m2 = mfcac+mfccc;
+// m1 = mfccc-mfcac;
+// m0 = m2+mfcbc;
+// mfcac = m0;
+// m0 += c1o18 * oMdrho;
+// mfcbc = m1-m0 * vvy;
+// mfccc = m2-2. * m1 * vvy+vy2 * m0;
+// ////////////////////////////////////////////////////////////////////////////////////
+// ////////////////////////////////////////////////////////////////////////////////////
+// // mit 1, 0, 1/3, 0, 0, 0, 1/3, 0, 1/9 Konditionieren
+// ////////////////////////////////////////////////////////////////////////////////////
+// // X - Dir
+// m2 = mfaaa+mfcaa;
+// m1 = mfcaa-mfaaa;
+// m0 = m2+mfbaa;
+// mfaaa = m0;
+// m0 += 1. * oMdrho;
+// mfbaa = m1-m0 * vvx;
+// mfcaa = m2-2. * m1 * vvx+vx2 * m0;
+// ////////////////////////////////////////////////////////////////////////////////////
+// m2 = mfaba+mfcba;
+// m1 = mfcba-mfaba;
+// m0 = m2+mfbba;
+// mfaba = m0;
+// mfbba = m1-m0 * vvx;
+// mfcba = m2-2. * m1 * vvx+vx2 * m0;
+// ////////////////////////////////////////////////////////////////////////////////////
+// m2 = mfaca+mfcca;
+// m1 = mfcca-mfaca;
+// m0 = m2+mfbca;
+// mfaca = m0;
+// m0 += c1o3 * oMdrho;
+// mfbca = m1-m0 * vvx;
+// mfcca = m2-2. * m1 * vvx+vx2 * m0;
+// ////////////////////////////////////////////////////////////////////////////////////
+// ////////////////////////////////////////////////////////////////////////////////////
+// m2 = mfaab+mfcab;
+// m1 = mfcab-mfaab;
+// m0 = m2+mfbab;
+// mfaab = m0;
+// mfbab = m1-m0 * vvx;
+// mfcab = m2-2. * m1 * vvx+vx2 * m0;
+// ////////////////////////////////////////////////////////////////////////////////////
+// m2 = mfabb+mfcbb;
+// m1 = mfcbb-mfabb;
+// m0 = m2+mfbbb;
+// mfabb = m0;
+// mfbbb = m1-m0 * vvx;
+// mfcbb = m2-2. * m1 * vvx+vx2 * m0;
+// ////////////////////////////////////////////////////////////////////////////////////
+// m2 = mfacb+mfccb;
+// m1 = mfccb-mfacb;
+// m0 = m2+mfbcb;
+// mfacb = m0;
+// mfbcb = m1-m0 * vvx;
+// mfccb = m2-2. * m1 * vvx+vx2 * m0;
+// ////////////////////////////////////////////////////////////////////////////////////
+// ////////////////////////////////////////////////////////////////////////////////////
+// m2 = mfaac+mfcac;
+// m1 = mfcac-mfaac;
+// m0 = m2+mfbac;
+// mfaac = m0;
+// m0 += c1o3 * oMdrho;
+// mfbac = m1-m0 * vvx;
+// mfcac = m2-2. * m1 * vvx+vx2 * m0;
+// ////////////////////////////////////////////////////////////////////////////////////
+// m2 = mfabc+mfcbc;
+// m1 = mfcbc-mfabc;
+// m0 = m2+mfbbc;
+// mfabc = m0;
+// mfbbc = m1-m0 * vvx;
+// mfcbc = m2-2. * m1 * vvx+vx2 * m0;
+// ////////////////////////////////////////////////////////////////////////////////////
+// m2 = mfacc+mfccc;
+// m1 = mfccc-mfacc;
+// m0 = m2+mfbcc;
+// mfacc = m0;
+// m0 += c1o9 * oMdrho;
+// mfbcc = m1-m0 * vvx;
+// mfccc = m2-2. * m1 * vvx+vx2 * m0;
+// ////////////////////////////////////////////////////////////////////////////////////
+// // Cumulants
+// ////////////////////////////////////////////////////////////////////////////////////
+// LBMReal OxxPyyPzz = 1.; //omega2 or bulk viscosity
+// LBMReal OxyyPxzz = 1.;//-s9;//2+s9;//
+// //LBMReal OxyyMxzz = 1.;//2+s9;//
+// LBMReal O4 = 1.;
+// LBMReal O5 = 1.;
+// LBMReal O6 = 1.;
+//
+// //Cum 4.
+// //LBMReal CUMcbb = mfcbb - ((mfcaa + c1o3 * oMdrho) * mfabb + 2. * mfbba * mfbab); // till 18.05.2015
+// //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 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);
+// LBMReal CUMacc = mfacc-((mfaac * mfaca+2. * mfabb * mfabb)+c1o3 * (mfaac+mfaca) * oMdrho+c1o9*(oMdrho-1)*oMdrho);
+//
+// //Cum 5.
+// LBMReal CUMbcc = mfbcc-(mfaac * mfbca+mfaca * mfbac+4. * mfabb * mfbbb+2. * (mfbab * mfacb+mfbba * mfabc))-c1o3 * (mfbca+mfbac) * oMdrho;
+// LBMReal CUMcbc = mfcbc-(mfaac * mfcba+mfcaa * mfabc+4. * mfbab * mfbbb+2. * (mfabb * mfcab+mfbba * mfbac))-c1o3 * (mfcba+mfabc) * oMdrho;
+// 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
+// -(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)
+// -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.
+// // linear combinations
+// LBMReal mxxPyyPzz = mfcaa+mfaca+mfaac;
+// LBMReal mxxMyy = mfcaa-mfaca;
+// LBMReal mxxMzz = mfcaa-mfaac;
+//
+// LBMReal dxux = -c1o2 * collFactor *(mxxMyy+mxxMzz)+c1o2 * OxxPyyPzz*(mfaaa-mxxPyyPzz);
+// LBMReal dyuy = dxux+collFactor * c3o2 * mxxMyy;
+// LBMReal dzuz = dxux+collFactor * c3o2 * mxxMzz;
+//
+// //relax
+// mxxPyyPzz += OxxPyyPzz*(mfaaa-mxxPyyPzz)-3. * (1.-c1o2 * OxxPyyPzz) * (vx2 * dxux+vy2 * dyuy+vz2 * dzuz);
+// mxxMyy += collFactor * (-mxxMyy)-3. * (1.-c1o2 * collFactor) * (vx2 * dxux-vy2 * dyuy);
+// mxxMzz += collFactor * (-mxxMzz)-3. * (1.-c1o2 * collFactor) * (vx2 * dxux-vz2 * dzuz);
+//
+// mfabb += collFactor * (-mfabb);
+// mfbab += collFactor * (-mfbab);
+// mfbba += collFactor * (-mfbba);
+//
+// // linear combinations back
+// mfcaa = c1o3 * (mxxMyy+mxxMzz+mxxPyyPzz);
+// mfaca = c1o3 * (-2. * mxxMyy+mxxMzz+mxxPyyPzz);
+// mfaac = c1o3 * (mxxMyy-2. * mxxMzz+mxxPyyPzz);
+//
+// //3.
+// // linear combinations
+// LBMReal mxxyPyzz = mfcba+mfabc;
+// LBMReal mxxyMyzz = mfcba-mfabc;
+//
+// LBMReal mxxzPyyz = mfcab+mfacb;
+// LBMReal mxxzMyyz = mfcab-mfacb;
+//
+// LBMReal mxyyPxzz = mfbca+mfbac;
+// LBMReal mxyyMxzz = mfbca-mfbac;
+//
+// //relax
+// wadjust = OxyyMxzz+(1.-OxyyMxzz)*fabs(mfbbb)/(fabs(mfbbb)+qudricLimit);
+// mfbbb += wadjust * (-mfbbb);
+// wadjust = OxyyPxzz+(1.-OxyyPxzz)*fabs(mxxyPyzz)/(fabs(mxxyPyzz)+qudricLimit);
+// mxxyPyzz += wadjust * (-mxxyPyzz);
+// wadjust = OxyyMxzz+(1.-OxyyMxzz)*fabs(mxxyMyzz)/(fabs(mxxyMyzz)+qudricLimit);
+// mxxyMyzz += wadjust * (-mxxyMyzz);
+// wadjust = OxyyPxzz+(1.-OxyyPxzz)*fabs(mxxzPyyz)/(fabs(mxxzPyyz)+qudricLimit);
+// mxxzPyyz += wadjust * (-mxxzPyyz);
+// wadjust = OxyyMxzz+(1.-OxyyMxzz)*fabs(mxxzMyyz)/(fabs(mxxzMyyz)+qudricLimit);
+// mxxzMyyz += wadjust * (-mxxzMyyz);
+// wadjust = OxyyPxzz+(1.-OxyyPxzz)*fabs(mxyyPxzz)/(fabs(mxyyPxzz)+qudricLimit);
+// mxyyPxzz += wadjust * (-mxyyPxzz);
+// wadjust = OxyyMxzz+(1.-OxyyMxzz)*fabs(mxyyMxzz)/(fabs(mxyyMxzz)+qudricLimit);
+// mxyyMxzz += wadjust * (-mxyyMxzz);
+//
+// // linear combinations back
+// mfcba = (mxxyMyzz+mxxyPyzz) * c1o2;
+// mfabc = (-mxxyMyzz+mxxyPyzz) * c1o2;
+// mfcab = (mxxzMyyz+mxxzPyyz) * c1o2;
+// mfacb = (-mxxzMyyz+mxxzPyyz) * c1o2;
+// mfbca = (mxyyMxzz+mxyyPxzz) * c1o2;
+// mfbac = (-mxyyMxzz+mxyyPxzz) * c1o2;
+//
+// //4.
+// CUMacc += O4 * (-CUMacc);
+// CUMcac += O4 * (-CUMcac);
+// CUMcca += O4 * (-CUMcca);
+//
+// CUMbbc += O4 * (-CUMbbc);
+// CUMbcb += O4 * (-CUMbcb);
+// CUMcbb += O4 * (-CUMcbb);
+//
+// //5.
+// CUMbcc += O5 * (-CUMbcc);
+// CUMcbc += O5 * (-CUMcbc);
+// CUMccb += O5 * (-CUMccb);
+//
+// //6.
+// CUMccc += O6 * (-CUMccc);
+//
+// //back cumulants to central moments
+// //4.
+// //mfcbb = CUMcbb + ((mfcaa + c1o3 * oMdrho) * mfabb + 2. * mfbba * mfbab); // till 18.05.2015
+// //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);
+//
+// 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;
+// mfacc = CUMacc+(mfaac * mfaca+2. * mfabb * mfabb)+c1o3 * (mfaac+mfaca) * oMdrho+c1o9*(oMdrho-1)*oMdrho;
+//
+// //5.
+// mfbcc = CUMbcc+(mfaac * mfbca+mfaca * mfbac+4. * mfabb * mfbbb+2. * (mfbab * mfacb+mfbba * mfabc))+c1o3 * (mfbca+mfbac) * oMdrho;
+// mfcbc = CUMcbc+(mfaac * mfcba+mfcaa * mfabc+4. * mfbab * mfbbb+2. * (mfabb * mfcab+mfbba * mfbac))+c1o3 * (mfcba+mfabc) * oMdrho;
+// mfccb = CUMccb+(mfcaa * mfacb+mfaca * mfcab+4. * mfbba * mfbbb+2. * (mfbab * mfbca+mfabb * mfcba))+c1o3 * (mfacb+mfcab) * oMdrho;
+//
+// //6.
+// 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)
+// -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;
+//
+// mfaab = 0.0;
+// mfaba = 0.0;
+// mfbaa = 0.0;
+//
+// //mfaab *= -0.5;
+// //mfaba *= -0.5;
+// //mfbaa *= -0.5;
+//
+//
+// ////////////////////////////////////////////////////////////////////////////////////
+// //back
+// ////////////////////////////////////////////////////////////////////////////////////
+// //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;
+// m1 = -mfaac-2. * mfaab * vvz+mfaaa * (1.-vz2)-1. * oMdrho * vz2;
+// 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;
+// m1 = -mfabc-2. * mfabb * vvz+mfaba * (1.-vz2);
+// 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;
+// m1 = -mfacc-2. * mfacb * vvz+mfaca * (1.-vz2)-c1o3 * oMdrho * vz2;
+// 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;
+// m1 = -mfbac-2. * mfbab * vvz+mfbaa * (1.-vz2);
+// 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;
+// m1 = -mfbbc-2. * mfbbb * vvz+mfbba * (1.-vz2);
+// 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;
+// m1 = -mfbcc-2. * mfbcb * vvz+mfbca * (1.-vz2);
+// 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;
+// m1 = -mfcac-2. * mfcab * vvz+mfcaa * (1.-vz2)-c1o3 * oMdrho * vz2;
+// 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;
+// m1 = -mfcbc-2. * mfcbb * vvz+mfcba * (1.-vz2);
+// 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;
+// m1 = -mfccc-2. * mfccb * vvz+mfcca * (1.-vz2)-c1o9 * oMdrho * vz2;
+// m2 = mfccc * c1o2+mfccb * (vvz+c1o2)+(mfcca+c1o9 * oMdrho) * (vz2+vvz) * c1o2;
+// mfcca = m0;
+// mfccb = m1;
+// mfccc = m2;
+// ////////////////////////////////////////////////////////////////////////////////////
+// ////////////////////////////////////////////////////////////////////////////////////
+// //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;
+// m1 = -mfaca-2. * mfaba * vvy+mfaaa * (1.-vy2)-c1o6 * oMdrho * vy2;
+// 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;
+// m1 = -mfacb-2. * mfabb * vvy+mfaab * (1.-vy2)-c2o3 * oMdrho * vy2;
+// 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;
+// m1 = -mfacc-2. * mfabc * vvy+mfaac * (1.-vy2)-c1o6 * oMdrho * vy2;
+// 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;
+// m1 = -mfbca-2. * mfbba * vvy+mfbaa * (1.-vy2);
+// 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;
+// m1 = -mfbcb-2. * mfbbb * vvy+mfbab * (1.-vy2);
+// 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;
+// m1 = -mfbcc-2. * mfbbc * vvy+mfbac * (1.-vy2);
+// 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;
+// m1 = -mfcca-2. * mfcba * vvy+mfcaa * (1.-vy2)-c1o18 * oMdrho * vy2;
+// 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;
+// m1 = -mfccb-2. * mfcbb * vvy+mfcab * (1.-vy2)-c2o9 * oMdrho * vy2;
+// 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;
+// m1 = -mfccc-2. * mfcbc * vvy+mfcac * (1.-vy2)-c1o18 * oMdrho * vy2;
+// m2 = mfccc * c1o2+mfcbc * (vvy+c1o2)+(mfcac+c1o18 * oMdrho) * (vy2+vvy) * c1o2;
+// mfcac = m0;
+// mfcbc = m1;
+// mfccc = m2;
+// ////////////////////////////////////////////////////////////////////////////////////
+// ////////////////////////////////////////////////////////////////////////////////////
+// //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;
+// m1 = -mfcaa-2. * mfbaa * vvx+mfaaa * (1.-vx2)-c1o36 * oMdrho * vx2;
+// 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;
+// m1 = -mfcba-2. * mfbba * vvx+mfaba * (1.-vx2)-c1o9 * oMdrho * vx2;
+// 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;
+// m1 = -mfcca-2. * mfbca * vvx+mfaca * (1.-vx2)-c1o36 * oMdrho * vx2;
+// 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;
+// m1 = -mfcab-2. * mfbab * vvx+mfaab * (1.-vx2)-c1o9 * oMdrho * vx2;
+// 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;
+// m1 = -mfcbb-2. * mfbbb * vvx+mfabb * (1.-vx2)-c4o9 * oMdrho * vx2;
+// 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;
+// m1 = -mfccb-2. * mfbcb * vvx+mfacb * (1.-vx2)-c1o9 * oMdrho * vx2;
+// 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;
+// m1 = -mfcac-2. * mfbac * vvx+mfaac * (1.-vx2)-c1o36 * oMdrho * vx2;
+// 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;
+// m1 = -mfcbc-2. * mfbbc * vvx+mfabc * (1.-vx2)-c1o9 * oMdrho * vx2;
+// 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;
+// m1 = -mfccc-2. * mfbcc * vvx+mfacc * (1.-vx2)-c1o36 * oMdrho * vx2;
+// m2 = mfccc * c1o2+mfbcc * (vvx+c1o2)+(mfacc+c1o36 * oMdrho) * (vx2+vvx) * c1o2;
+// mfacc = m0;
+// mfbcc = m1;
+// mfccc = m2;
+//
+// //////////////////////////////////////////////////////////////////////////
+// //proof correctness
+// //////////////////////////////////////////////////////////////////////////
+//#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;
+// //LBMReal dif = fabs(rho - rho_post);
+// LBMReal dif = rho-rho_post;
+//#ifdef SINGLEPRECISION
+// if (dif>10.0E-7||dif<-10.0E-7)
+//#else
+// 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)
+// +" dif="+UbSystem::toString(dif)
+// +" rho is not correct for node "+UbSystem::toString(x1)+","+UbSystem::toString(x2)+","+UbSystem::toString(x3)));
+// //UBLOG(logERROR,"LBMKernelETD3Q27CCLB::collideAll(): rho is not correct for node "+UbSystem::toString(x1)+","+UbSystem::toString(x2)+","+UbSystem::toString(x3));
+// //exit(EXIT_FAILURE);
+// }
+//#endif
+// ////////////////////////////////////////////////////////////////////////////
+// ////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;
+// ////////////////////////////////////////////////////////////////////////////
+//
+// }
+// }
+// }
+// }
+//
+//}
+
+void InitDensityLBMKernel::init()
+{
+ DistributionArray3DPtr 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()
+{
+ using namespace D3Q27System;
+
+ localDistributions = boost::dynamic_pointer_cast<D3Q27EsoTwist3DSplittedVector>(dataSet->getFdistributions())->getLocalDistributions();
+ nonLocalDistributions = boost::dynamic_pointer_cast<D3Q27EsoTwist3DSplittedVector>(dataSet->getFdistributions())->getNonLocalDistributions();
+ zeroDistributions = boost::dynamic_pointer_cast<D3Q27EsoTwist3DSplittedVector>(dataSet->getFdistributions())->getZeroDistributions();
+
+ BCArray3D<D3Q27BoundaryCondition>& bcArray = boost::dynamic_pointer_cast<D3Q27ETBCProcessor>(this->getBCProcessor())->getBCArray();
+ D3Q27BoundaryConditionPtr bcPtr;
+ LBMReal f[D3Q27System::ENDF+1];
+ LBMReal feq[D3Q27System::ENDF+1];
+ LBMReal drho, vx1, vx2, vx3;
+ const int bcArrayMaxX1 = (int)bcArray.getNX1();
+ const int bcArrayMaxX2 = (int)bcArray.getNX2();
+ const int bcArrayMaxX3 = (int)bcArray.getNX3();
+
+ int minX1 = ghostLayerWidth;
+ int minX2 = ghostLayerWidth;
+ int minX3 = ghostLayerWidth;
+ int maxX1 = bcArrayMaxX1-ghostLayerWidth;
+ int maxX2 = bcArrayMaxX2-ghostLayerWidth;
+ int maxX3 = bcArrayMaxX3-ghostLayerWidth;
+
+ //collFactor = 1.0/(1.0/2.0+1.0/sqrt(6.0));
+ collFactor = 1.0;
+
+ for (int x3 = minX3; x3<maxX3; x3++)
+ {
+ for (int x2 = minX2; x2<maxX2; x2++)
+ {
+ for (int x1 = minX1; x1<maxX1; x1++)
+ {
+ if (!bcArray.isSolid(x1, x2, x3)&&!bcArray.isUndefined(x1, x2, x3))
+ {
+ int x1p = x1+1;
+ int x2p = x2+1;
+ int x3p = x3+1;
+ //////////////////////////////////////////////////////////////////////////
+ //read distribution
+ ////////////////////////////////////////////////////////////////////////////
+ f[ZERO] = (*this->zeroDistributions)(x1, x2, x3);
+
+ f[E] = (*this->localDistributions)(D3Q27System::ET_E, x1, x2, x3);
+ f[N] = (*this->localDistributions)(D3Q27System::ET_N, x1, x2, x3);
+ f[T] = (*this->localDistributions)(D3Q27System::ET_T, x1, x2, x3);
+ f[NE] = (*this->localDistributions)(D3Q27System::ET_NE, x1, x2, x3);
+ f[NW] = (*this->localDistributions)(D3Q27System::ET_NW, x1p, x2, x3);
+ f[TE] = (*this->localDistributions)(D3Q27System::ET_TE, x1, x2, x3);
+ f[TW] = (*this->localDistributions)(D3Q27System::ET_TW, x1p, x2, x3);
+ f[TN] = (*this->localDistributions)(D3Q27System::ET_TN, x1, x2, x3);
+ f[TS] = (*this->localDistributions)(D3Q27System::ET_TS, x1, x2p, x3);
+ f[TNE] = (*this->localDistributions)(D3Q27System::ET_TNE, x1, x2, x3);
+ f[TNW] = (*this->localDistributions)(D3Q27System::ET_TNW, x1p, x2, x3);
+ f[TSE] = (*this->localDistributions)(D3Q27System::ET_TSE, x1, x2p, x3);
+ f[TSW] = (*this->localDistributions)(D3Q27System::ET_TSW, x1p, x2p, x3);
+
+ f[W] = (*this->nonLocalDistributions)(D3Q27System::ET_W, x1p, x2, x3);
+ f[S] = (*this->nonLocalDistributions)(D3Q27System::ET_S, x1, x2p, x3);
+ f[B] = (*this->nonLocalDistributions)(D3Q27System::ET_B, x1, x2, x3p);
+ f[SW] = (*this->nonLocalDistributions)(D3Q27System::ET_SW, x1p, x2p, x3);
+ f[SE] = (*this->nonLocalDistributions)(D3Q27System::ET_SE, x1, x2p, x3);
+ f[BW] = (*this->nonLocalDistributions)(D3Q27System::ET_BW, x1p, x2, x3p);
+ f[BE] = (*this->nonLocalDistributions)(D3Q27System::ET_BE, x1, x2, x3p);
+ f[BS] = (*this->nonLocalDistributions)(D3Q27System::ET_BS, x1, x2p, x3p);
+ f[BN] = (*this->nonLocalDistributions)(D3Q27System::ET_BN, x1, x2, x3p);
+ f[BSW] = (*this->nonLocalDistributions)(D3Q27System::ET_BSW, x1p, x2p, x3p);
+ f[BSE] = (*this->nonLocalDistributions)(D3Q27System::ET_BSE, x1, x2p, x3p);
+ f[BNW] = (*this->nonLocalDistributions)(D3Q27System::ET_BNW, x1p, x2, x3p);
+ f[BNE] = (*this->nonLocalDistributions)(D3Q27System::ET_BNE, x1, x2, x3p);
+ //////////////////////////////////////////////////////////////////////////
+
+ drho = ((f[TNE]+f[BSW])+(f[TSE]+f[BNW]))+((f[BSE]+f[TNW])+(f[TSW]+f[BNE]))
+ +(((f[NE]+f[SW])+(f[SE]+f[NW]))+((f[TE]+f[BW])+(f[BE]+f[TW]))
+ +((f[BN]+f[TS])+(f[TN]+f[BS])))+((f[E]+f[W])+(f[N]+f[S])
+ +(f[T]+f[B]))+f[ZERO];
+
+ //vx1 = ((((f[TNE]-f[BSW])+(f[TSE]-f[BNW]))+((f[BSE]-f[TNW])+(f[BNE]-f[TSW])))+
+ // (((f[BE]-f[TW])+(f[TE]-f[BW]))+((f[SE]-f[NW])+(f[NE]-f[SW])))+
+ // (f[E]-f[W]));
+
+ //vx2 = ((((f[TNE]-f[BSW])+(f[BNW]-f[TSE]))+((f[TNW]-f[BSE])+(f[BNE]-f[TSW])))+
+ // (((f[BN]-f[TS])+(f[TN]-f[BS]))+((f[NW]-f[SE])+(f[NE]-f[SW])))+
+ // (f[N]-f[S]));
+
+ //vx3 = ((((f[TNE]-f[BSW])+(f[TSE]-f[BNW]))+((f[TNW]-f[BSE])+(f[TSW]-f[BNE])))+
+ // (((f[TS]-f[BN])+(f[TN]-f[BS]))+((f[TW]-f[BE])+(f[TE]-f[BW])))+
+ // (f[T]-f[B]));
+
+ vx1 = v(0,x1,x2,x3);
+ vx2 = v(1,x1,x2,x3);
+ vx3 = v(2,x1,x2,x3);
+
+ //LBMReal vvx = v(0,x1,x2,x3);
+ //LBMReal vvy = v(1,x1,x2,x3);
+ //LBMReal vvz = v(2,x1,x2,x3);
+
+ //vx1 = vx1+(vvx-vx1);
+ //vx2 = vx2+(vvy-vx2);
+ //vx3 = vx3+(vvz-vx3);
+
+ LBMReal cu_sq = 1.5*(vx1*vx1+vx2*vx2+vx3*vx3);
+
+ feq[ZERO] = c8o27*(drho-cu_sq);
+ feq[E] = c2o27*(drho+3.0*(vx1)+c9o2*(vx1)*(vx1)-cu_sq);
+ feq[W] = c2o27*(drho+3.0*(-vx1)+c9o2*(-vx1)*(-vx1)-cu_sq);
+ feq[N] = c2o27*(drho+3.0*(vx2)+c9o2*(vx2)*(vx2)-cu_sq);
+ feq[S] = c2o27*(drho+3.0*(-vx2)+c9o2*(-vx2)*(-vx2)-cu_sq);
+ feq[T] = c2o27*(drho+3.0*(vx3)+c9o2*(vx3)*(vx3)-cu_sq);
+ feq[B] = c2o27*(drho+3.0*(-vx3)+c9o2*(-vx3)*(-vx3)-cu_sq);
+ feq[NE] = c1o54*(drho+3.0*(vx1+vx2)+c9o2*(vx1+vx2)*(vx1+vx2)-cu_sq);
+ feq[SW] = c1o54*(drho+3.0*(-vx1-vx2)+c9o2*(-vx1-vx2)*(-vx1-vx2)-cu_sq);
+ feq[SE] = c1o54*(drho+3.0*(vx1-vx2)+c9o2*(vx1-vx2)*(vx1-vx2)-cu_sq);
+ feq[NW] = c1o54*(drho+3.0*(-vx1+vx2)+c9o2*(-vx1+vx2)*(-vx1+vx2)-cu_sq);
+ feq[TE] = c1o54*(drho+3.0*(vx1+vx3)+c9o2*(vx1+vx3)*(vx1+vx3)-cu_sq);
+ feq[BW] = c1o54*(drho+3.0*(-vx1-vx3)+c9o2*(-vx1-vx3)*(-vx1-vx3)-cu_sq);
+ feq[BE] = c1o54*(drho+3.0*(vx1-vx3)+c9o2*(vx1-vx3)*(vx1-vx3)-cu_sq);
+ feq[TW] = c1o54*(drho+3.0*(-vx1+vx3)+c9o2*(-vx1+vx3)*(-vx1+vx3)-cu_sq);
+ feq[TN] = c1o54*(drho+3.0*(vx2+vx3)+c9o2*(vx2+vx3)*(vx2+vx3)-cu_sq);
+ feq[BS] = c1o54*(drho+3.0*(-vx2-vx3)+c9o2*(-vx2-vx3)*(-vx2-vx3)-cu_sq);
+ feq[BN] = c1o54*(drho+3.0*(vx2-vx3)+c9o2*(vx2-vx3)*(vx2-vx3)-cu_sq);
+ feq[TS] = c1o54*(drho+3.0*(-vx2+vx3)+c9o2*(-vx2+vx3)*(-vx2+vx3)-cu_sq);
+ feq[TNE] = c1o216*(drho+3.0*(vx1+vx2+vx3)+c9o2*(vx1+vx2+vx3)*(vx1+vx2+vx3)-cu_sq);
+ feq[BSW] = c1o216*(drho+3.0*(-vx1-vx2-vx3)+c9o2*(-vx1-vx2-vx3)*(-vx1-vx2-vx3)-cu_sq);
+ feq[BNE] = c1o216*(drho+3.0*(vx1+vx2-vx3)+c9o2*(vx1+vx2-vx3)*(vx1+vx2-vx3)-cu_sq);
+ feq[TSW] = c1o216*(drho+3.0*(-vx1-vx2+vx3)+c9o2*(-vx1-vx2+vx3)*(-vx1-vx2+vx3)-cu_sq);
+ feq[TSE] = c1o216*(drho+3.0*(vx1-vx2+vx3)+c9o2*(vx1-vx2+vx3)*(vx1-vx2+vx3)-cu_sq);
+ feq[BNW] = c1o216*(drho+3.0*(-vx1+vx2-vx3)+c9o2*(-vx1+vx2-vx3)*(-vx1+vx2-vx3)-cu_sq);
+ feq[BSE] = c1o216*(drho+3.0*(vx1-vx2-vx3)+c9o2*(vx1-vx2-vx3)*(vx1-vx2-vx3)-cu_sq);
+ feq[TNW] = c1o216*(drho+3.0*(-vx1+vx2+vx3)+c9o2*(-vx1+vx2+vx3)*(-vx1+vx2+vx3)-cu_sq);
+
+ //Relaxation
+ f[ZERO] += (feq[ZERO]-f[ZERO])*collFactor;
+ f[E] += (feq[E]-f[E])*collFactor;
+ f[W] += (feq[W]-f[W])*collFactor;
+ f[N] += (feq[N]-f[N])*collFactor;
+ f[S] += (feq[S]-f[S])*collFactor;
+ f[T] += (feq[T]-f[T])*collFactor;
+ f[B] += (feq[B]-f[B])*collFactor;
+ f[NE] += (feq[NE]-f[NE])*collFactor;
+ f[SW] += (feq[SW]-f[SW])*collFactor;
+ f[SE] += (feq[SE]-f[SE])*collFactor;
+ f[NW] += (feq[NW]-f[NW])*collFactor;
+ f[TE] += (feq[TE]-f[TE])*collFactor;
+ f[BW] += (feq[BW]-f[BW])*collFactor;
+ f[BE] += (feq[BE]-f[BE])*collFactor;
+ f[TW] += (feq[TW]-f[TW])*collFactor;
+ f[TN] += (feq[TN]-f[TN])*collFactor;
+ f[BS] += (feq[BS]-f[BS])*collFactor;
+ f[BN] += (feq[BN]-f[BN])*collFactor;
+ f[TS] += (feq[TS]-f[TS])*collFactor;
+
+ f[TNE] += (feq[TNE]-f[TNE])*collFactor;
+ f[BSW] += (feq[BSW]-f[BSW])*collFactor;
+ f[BNE] += (feq[BNE]-f[BNE])*collFactor;
+ f[TSW] += (feq[TSW]-f[TSW])*collFactor;
+ f[TSE] += (feq[TSE]-f[TSE])*collFactor;
+ f[BNW] += (feq[BNW]-f[BNW])*collFactor;
+ f[BSE] += (feq[BSE]-f[BSE])*collFactor;
+ f[TNW] += (feq[TNW]-f[TNW])*collFactor;
+
+ //////////////////////////////////////////////////////////////////////////
+#ifdef PROOF_CORRECTNESS
+ LBMReal rho_post = f[ZERO]+f[E]+f[W]+f[N]+f[S]+f[T]+f[B]
+ +f[NE]+f[SW]+f[SE]+f[NW]+f[TE]+f[BW]+f[BE]
+ +f[TW]+f[TN]+f[BS]+f[BN]+f[TS]+f[TNE]+f[TSW]
+ +f[TSE]+f[TNW]+f[BNE]+f[BSW]+f[BSE]+f[BNW];
+ LBMReal dif = drho-rho_post;
+#ifdef SINGLEPRECISION
+ if (dif>10.0E-7||dif<-10.0E-7)
+#else
+ if (dif>10.0E-15||dif<-10.0E-15)
+#endif
+ {
+ UB_THROW(UbException(UB_EXARGS, "rho is not correct"));
+ }
+#endif
+ //////////////////////////////////////////////////////////////////////////
+ //write distribution
+ //////////////////////////////////////////////////////////////////////////
+ (*this->localDistributions)(D3Q27System::ET_E, x1, x2, x3) = f[D3Q27System::INV_E];
+ (*this->localDistributions)(D3Q27System::ET_N, x1, x2, x3) = f[D3Q27System::INV_N];
+ (*this->localDistributions)(D3Q27System::ET_T, x1, x2, x3) = f[D3Q27System::INV_T];
+ (*this->localDistributions)(D3Q27System::ET_NE, x1, x2, x3) = f[D3Q27System::INV_NE];
+ (*this->localDistributions)(D3Q27System::ET_NW, x1p, x2, x3) = f[D3Q27System::INV_NW];
+ (*this->localDistributions)(D3Q27System::ET_TE, x1, x2, x3) = f[D3Q27System::INV_TE];
+ (*this->localDistributions)(D3Q27System::ET_TW, x1p, x2, x3) = f[D3Q27System::INV_TW];
+ (*this->localDistributions)(D3Q27System::ET_TN, x1, x2, x3) = f[D3Q27System::INV_TN];
+ (*this->localDistributions)(D3Q27System::ET_TS, x1, x2p, x3) = f[D3Q27System::INV_TS];
+ (*this->localDistributions)(D3Q27System::ET_TNE, x1, x2, x3) = f[D3Q27System::INV_TNE];
+ (*this->localDistributions)(D3Q27System::ET_TNW, x1p, x2, x3) = f[D3Q27System::INV_TNW];
+ (*this->localDistributions)(D3Q27System::ET_TSE, x1, x2p, x3) = f[D3Q27System::INV_TSE];
+ (*this->localDistributions)(D3Q27System::ET_TSW, x1p, x2p, x3) = f[D3Q27System::INV_TSW];
+
+ (*this->nonLocalDistributions)(D3Q27System::ET_W, x1p, x2, x3) = f[D3Q27System::INV_W];
+ (*this->nonLocalDistributions)(D3Q27System::ET_S, x1, x2p, x3) = f[D3Q27System::INV_S];
+ (*this->nonLocalDistributions)(D3Q27System::ET_B, x1, x2, x3p) = f[D3Q27System::INV_B];
+ (*this->nonLocalDistributions)(D3Q27System::ET_SW, x1p, x2p, x3) = f[D3Q27System::INV_SW];
+ (*this->nonLocalDistributions)(D3Q27System::ET_SE, x1, x2p, x3) = f[D3Q27System::INV_SE];
+ (*this->nonLocalDistributions)(D3Q27System::ET_BW, x1p, x2, x3p) = f[D3Q27System::INV_BW];
+ (*this->nonLocalDistributions)(D3Q27System::ET_BE, x1, x2, x3p) = f[D3Q27System::INV_BE];
+ (*this->nonLocalDistributions)(D3Q27System::ET_BS, x1, x2p, x3p) = f[D3Q27System::INV_BS];
+ (*this->nonLocalDistributions)(D3Q27System::ET_BN, x1, x2, x3p) = f[D3Q27System::INV_BN];
+ (*this->nonLocalDistributions)(D3Q27System::ET_BSW, x1p, x2p, x3p) = f[D3Q27System::INV_BSW];
+ (*this->nonLocalDistributions)(D3Q27System::ET_BSE, x1, x2p, x3p) = f[D3Q27System::INV_BSE];
+ (*this->nonLocalDistributions)(D3Q27System::ET_BNW, x1p, x2, x3p) = f[D3Q27System::INV_BNW];
+ (*this->nonLocalDistributions)(D3Q27System::ET_BNE, x1, x2, x3p) = f[D3Q27System::INV_BNE];
+
+ (*this->zeroDistributions)(x1, x2, x3) = f[D3Q27System::ZERO];
+ //////////////////////////////////////////////////////////////////////////
+
+
+ }
+ }
+ }
+ }
+}
\ No newline at end of file
diff --git a/source/VirtualFluidsCore/LBM/InitDensityLBMKernel.h b/source/VirtualFluidsCore/LBM/InitDensityLBMKernel.h
new file mode 100644
index 0000000000000000000000000000000000000000..e73a85c4da77e51ff51d190ce2855c481f9fabf2
--- /dev/null
+++ b/source/VirtualFluidsCore/LBM/InitDensityLBMKernel.h
@@ -0,0 +1,30 @@
+#ifndef InitDensityLBMKernel_h__
+#define InitDensityLBMKernel_h__
+
+#include "LBMKernelETD3Q27.h"
+#include "basics/utilities/UbTiming.h"
+
+class InitDensityLBMKernel : public LBMKernelETD3Q27
+{
+public:
+ InitDensityLBMKernel();
+ ~InitDensityLBMKernel();
+ InitDensityLBMKernel(int nx1, int nx2, int nx3);
+ void calculate();
+ LBMKernel3DPtr clone();
+ void setVelocity(int x1, int x2, int x3, LBMReal vvx, LBMReal vvy, LBMReal vvz);
+ double getCallculationTime();
+protected:
+ void init();
+ void collideAll();
+private:
+ LBMReal f[D3Q27System::ENDF+1];
+ CbArray4D<LBMReal, IndexerX4X3X2X1>::CbArray4DPtr localDistributions;
+ CbArray4D<LBMReal, IndexerX4X3X2X1>::CbArray4DPtr nonLocalDistributions;
+ CbArray3D<LBMReal, IndexerX3X2X1>::CbArray3DPtr zeroDistributions;
+ LBMReal OxyyMxzz;
+ CbArray4D<LBMReal, IndexerX4X3X2X1> v;
+};
+
+#endif // InitDensityLBMKernel_h__
+
diff --git a/source/VirtualFluidsCore/LBM/LBMKernel3D.h b/source/VirtualFluidsCore/LBM/LBMKernel3D.h
index 011a6bd037b3d175cbe3c510b56f5702e1728938..e90da3ee6f5feba76b9264f9bbb0619a88023533 100644
--- a/source/VirtualFluidsCore/LBM/LBMKernel3D.h
+++ b/source/VirtualFluidsCore/LBM/LBMKernel3D.h
@@ -33,14 +33,14 @@ public:
virtual void swapDistributions() = 0;
virtual double getCallculationTime() = 0;
- virtual void setBCProcessor(BCProcessorPtr bcp);
- virtual BCProcessorPtr getBCProcessor();
+ void setBCProcessor(BCProcessorPtr bcp);
+ BCProcessorPtr getBCProcessor();
- virtual void setCollisionFactor(double collFactor);
- virtual double getCollisionFactor() const;
+ void setCollisionFactor(double collFactor);
+ double getCollisionFactor() const;
- virtual void setGhostLayerWidth(int witdh);
- virtual int getGhostLayerWidth() const;
+ void setGhostLayerWidth(int witdh);
+ int getGhostLayerWidth() const;
void setDataSet(DataSet3DPtr dataSet);
DataSet3DPtr getDataSet() const;
diff --git a/source/VirtualFluidsCore/Utilities/ChangeRandomQs.hpp b/source/VirtualFluidsCore/Utilities/ChangeRandomQs.hpp
index a4adf5a266d0ffd9f38eba5200aa51aa0eb29620..9a91f578fc913f470c0d877e14a4ddc5df3fbf40 100644
--- a/source/VirtualFluidsCore/Utilities/ChangeRandomQs.hpp
+++ b/source/VirtualFluidsCore/Utilities/ChangeRandomQs.hpp
@@ -11,9 +11,9 @@ namespace Utilities
{
void ChangeRandomQs(D3Q27IntegrateValuesHelperPtr integrateValues)
{
- std::vector<CalcNodes> cnodes = integrateValues->getCNodes();
+ std::vector<D3Q27IntegrateValuesHelper::CalcNodes> cnodes = integrateValues->getCNodes();
- BOOST_FOREACH(CalcNodes cn, cnodes)
+ BOOST_FOREACH(D3Q27IntegrateValuesHelper::CalcNodes cn, cnodes)
{
LBMKernel3DPtr kernel = cn.block->getKernel();
BCArray3D<D3Q27BoundaryCondition>& bcArray = boost::dynamic_pointer_cast<D3Q27ETBCProcessor>(kernel->getBCProcessor())->getBCArray();
diff --git a/source/VirtualFluidsCore/Visitors/InitDistributionsFromFileBlockVisitor.cpp b/source/VirtualFluidsCore/Visitors/InitDistributionsFromFileBlockVisitor.cpp
index 9e9236c0b02d40b5765f68de26638aaa789f140c..98b3dd7d1e8a987be6e4fcc14a019c8d3db5da42 100644
--- a/source/VirtualFluidsCore/Visitors/InitDistributionsFromFileBlockVisitor.cpp
+++ b/source/VirtualFluidsCore/Visitors/InitDistributionsFromFileBlockVisitor.cpp
@@ -3,14 +3,15 @@
#include "LBMKernel3D.h"
#include "D3Q27ETBCProcessor.h"
#include "Grid3DSystem.h"
+#include "InitDensityLBMKernel.h"
-InitDistributionsFromFileBlockVisitor::InitDistributionsFromFileBlockVisitor(LBMReal nu, LBMReal rho, std::string filename)
+InitDistributionsFromFileBlockVisitor::InitDistributionsFromFileBlockVisitor(LBMReal nu, LBMReal rho, std::string filename)
: Block3DVisitor(0, Grid3DSystem::MAXLEVEL), nu(nu), rho(rho)
{
UbFileInputASCII in(filename);
if (!in)
{
- throw UbException(UB_EXARGS, "could not open file " + filename);
+ throw UbException(UB_EXARGS, "could not open file "+filename);
}
int nodesX1 = in.readInteger();
@@ -19,9 +20,9 @@ InitDistributionsFromFileBlockVisitor::InitDistributionsFromFileBlockVisitor(LBM
matrix = CbArray4D<LBMReal, IndexerX4X3X2X1>(3, nodesX1, nodesX2, nodesX3, 0);
- for (int x3 = 0; x3 < nodesX3; x3++)
- for (int x2 = 0; x2 < nodesX2; x2++)
- for (int x1 = 0; x1 < nodesX1; x1++)
+ for (int x3 = 0; x3<nodesX3; x3++)
+ for (int x2 = 0; x2<nodesX2; x2++)
+ for (int x1 = 0; x1<nodesX1; x1++)
{
matrix(Vx1, x1, x2, x3) = in.readDouble();
matrix(Vx2, x1, x2, x3) = in.readDouble();
@@ -49,16 +50,16 @@ void InitDistributionsFromFileBlockVisitor::visit(const Grid3DPtr grid, Block3DP
typedef void(*CalcFeqsFct)(LBMReal* const& /*feq[27]*/, const LBMReal& /*(d)rho*/, const LBMReal& /*vx1*/, const LBMReal& /*vx2*/, const LBMReal& /*vx3*/);
CalcFeqsFct calcFeqsFct = NULL;
- LBMReal vx1, vx2, vx3, rho;
+ LBMReal vx1, vx2, vx3;
int gridRank = grid->getRank();
int blockRank = block->getRank();
- if (blockRank == gridRank && block->isActive())
+ if (blockRank==gridRank && block->isActive())
{
LBMKernel3DPtr kernel = block->getKernel();
if (!kernel)
- throw UbException(UB_EXARGS, "The LBM kernel isn't exist in block: " + block->toString());
+ throw UbException(UB_EXARGS, "The LBM kernel isn't exist in block: "+block->toString());
if (kernel->getCompressible())
calcFeqsFct = &D3Q27System::calcCompFeq;
@@ -70,7 +71,7 @@ void InitDistributionsFromFileBlockVisitor::visit(const Grid3DPtr grid, Block3DP
BCArray3D<D3Q27BoundaryCondition> bcArray = boost::dynamic_pointer_cast<D3Q27ETBCProcessor>(kernel->getBCProcessor())->getBCArray();
EsoTwist3DPtr distributions = boost::dynamic_pointer_cast<EsoTwist3D>(kernel->getDataSet()->getFdistributions());
- LBMReal f[D3Q27System::ENDF + 1];
+ LBMReal f[D3Q27System::ENDF+1];
size_t nx1 = distributions->getNX1();
size_t nx2 = distributions->getNX2();
@@ -84,128 +85,164 @@ void InitDistributionsFromFileBlockVisitor::visit(const Grid3DPtr grid, Block3DP
int maxX2 = (int)bcArray.getNX2();
int maxX3 = (int)bcArray.getNX3();
- for (int ix3 = minX1; ix3 < maxX3; ix3++)
- for (int ix2 = minX2; ix2 < maxX2; ix2++)
- for (int ix1 = minX1; ix1 < maxX1; ix1++)
+ int maxMX1 = matrix.getNX2();
+ int maxMX2 = matrix.getNX3();
+ int maxMX3 = matrix.getNX4();
+
+ int blockix1 = block->getX1();
+ int blockix2 = block->getX2();
+ int blockix3 = block->getX3();
+
+ UbTupleInt3 blockNx = grid->getBlockNX();
+
+ for (int ix3 = minX1; ix3<maxX3; ix3++)
+ for (int ix2 = minX2; ix2<maxX2; ix2++)
+ for (int ix1 = minX1; ix1<maxX1; ix1++)
{
- UbTupleInt3 coords = grid->getNodeIndexes(block, ix1, ix2, ix3);
- int x1 = val<1>(coords);
- int x2 = val<2>(coords);
- int x3 = val<3>(coords);
+ int x1 = blockix1*val<1>(blockNx)+ix1-1;
+ int x2 = blockix2*val<2>(blockNx)+ix2-1;
+ int x3 = blockix3*val<3>(blockNx)+ix3-1;
+
+ if (x1==-1)
+ {
+ x1 = maxMX1-1;
+ }
+ if (x2==-1)
+ {
+ x2 = maxMX2-1;
+ }
+ if (x3==-1)
+ {
+ x3 = maxMX3-1;
+ }
+
+ if (x1==maxMX1)
+ {
+ x1 = 1;
+ }
+ if (x2==maxMX2)
+ {
+ x2 = 1;
+ }
+ if (x3==maxMX3)
+ {
+ x3 = 1;
+ }
vx1 = matrix(Vx1, x1, x2, x3);
vx2 = matrix(Vx2, x1, x2, x3);
vx3 = matrix(Vx3, x1, x2, x3);
-
- //x-derivative
- //double deltaX = dx*0.5;
- int deltaX = 1;
- x1 = val<1>(coords) + deltaX;
- if (x1 > maxX1) x1 = val<1>(coords);
- double vx1Plusx1 = matrix(Vx1, x1, x2, x3);
- double vx2Plusx1 = matrix(Vx2, x1, x2, x3);
- double vx3Plusx1 = matrix(Vx3, x1, x2, x3);
-
- x1 = val<1>(coords) - deltaX;
- if (x1 < minX1) x1 = val<1>(coords);
- double vx1Minusx1 = matrix(Vx1, x1, x2, x3);
- double vx2Minusx1 = matrix(Vx2, x1, x2, x3);
- double vx3Minusx1 = matrix(Vx3, x1, x2, x3);
-
- //y-derivative
- x1 = val<1>(coords);
- x2 = val<2>(coords) + deltaX;
- if (x2 > maxX2) x2 = val<2>(coords);
- double vx1Plusx2 = matrix(Vx1, x1, x2, x3);
- double vx2Plusx2 = matrix(Vx2, x1, x2, x3);
- double vx3Plusx2 = matrix(Vx3, x1, x2, x3);
-
- x2 = val<2>(coords) - deltaX;
- if (x2 < minX2) x2 = val<2>(coords);
- double vx1Minusx2 = matrix(Vx1, x1, x2, x3);
- double vx2Minusx2 = matrix(Vx2, x1, x2, x3);
- double vx3Minusx2 = matrix(Vx3, x1, x2, x3);
-
- //z-derivative
- x2 = val<2>(coords);
- x3 = val<3>(coords) + deltaX;
- if (x3 > maxX3) x3 = val<3>(coords);
- double vx1Plusx3 = matrix(Vx1, x1, x2, x3);
- double vx2Plusx3 = matrix(Vx2, x1, x2, x3);
- double vx3Plusx3 = matrix(Vx3, x1, x2, x3);
-
- x3 = val<3>(coords) - deltaX;
- if (x3 < minX3) x3 = val<3>(coords);
- double vx1Minusx3 = matrix(Vx1, x1, x2, x3);
- double vx2Minusx3 = matrix(Vx2, x1, x2, x3);
- double vx3Minusx3 = matrix(Vx3, x1, x2, x3);
-
- double ax = (vx1Plusx1 - vx1Minusx1) / (2.0*deltaX);
- double bx = (vx2Plusx1 - vx2Minusx1) / (2.0*deltaX);
- double cx = (vx3Plusx1 - vx3Minusx1) / (2.0*deltaX);
-
- double ay = (vx1Plusx2 - vx1Minusx2) / (2.0*deltaX);
- double by = (vx2Plusx2 - vx2Minusx2) / (2.0*deltaX);
- double cy = (vx3Plusx2 - vx3Minusx2) / (2.0*deltaX);
-
- double az = (vx1Plusx3 - vx1Minusx3) / (2.0*deltaX);
- double bz = (vx2Plusx3 - vx2Minusx3) / (2.0*deltaX);
- double cz = (vx3Plusx3 - vx3Minusx3) / (2.0*deltaX);
- double eps_new = 1.0;
- LBMReal op = 1.;
-
- LBMReal feq[27];
-
- calcFeqsFct(feq, rho, vx1, vx2, vx3);
-
- double f_E = eps_new *((5.*ax*o + 5.*by*o + 5.*cz*o - 8.*ax*op + 4.*by*op + 4.*cz*op) / (54.*o*op));
- double f_N = f_E + eps_new *((2.*(ax - by)) / (9.*o));
- double f_T = f_E + eps_new *((2.*(ax - cz)) / (9.*o));
- double f_NE = eps_new *(-(5.*cz*o + 3.*(ay + bx)*op - 2.*cz*op + ax*(5.*o + op) + by*(5.*o + op)) / (54.*o*op));
- double f_SE = f_NE + eps_new *((ay + bx) / (9.*o));
- double f_TE = eps_new *(-(5.*cz*o + by*(5.*o - 2.*op) + 3.*(az + cx)*op + cz*op + ax*(5.*o + op)) / (54.*o*op));
- double f_BE = f_TE + eps_new *((az + cx) / (9.*o));
- double f_TN = eps_new *(-(5.*ax*o + 5.*by*o + 5.*cz*o - 2.*ax*op + by*op + 3.*bz*op + 3.*cy*op + cz*op) / (54.*o*op));
- double f_BN = f_TN + eps_new *((bz + cy) / (9.*o));
- double f_ZERO = eps_new *((5.*(ax + by + cz)) / (9.*op));
- double f_TNE = eps_new *(-(ay + az + bx + bz + cx + cy) / (72.*o));
- double f_TSW = -eps_new *((ay + bx) / (36.*o)) - f_TNE;
- double f_TSE = -eps_new *((az + cx) / (36.*o)) - f_TNE;
- double f_TNW = -eps_new *((bz + cy) / (36.*o)) - f_TNE;
-
-
- f[E] = f_E + feq[E];
- f[W] = f_E + feq[W];
- f[N] = f_N + feq[N];
- f[S] = f_N + feq[S];
- f[T] = f_T + feq[T];
- f[B] = f_T + feq[B];
- f[NE] = f_NE + feq[NE];
- f[SW] = f_NE + feq[SW];
- f[SE] = f_SE + feq[SE];
- f[NW] = f_SE + feq[NW];
- f[TE] = f_TE + feq[TE];
- f[BW] = f_TE + feq[BW];
- f[BE] = f_BE + feq[BE];
- f[TW] = f_BE + feq[TW];
- f[TN] = f_TN + feq[TN];
- f[BS] = f_TN + feq[BS];
- f[BN] = f_BN + feq[BN];
- f[TS] = f_BN + feq[TS];
- f[TNE] = f_TNE + feq[TNE];
- f[TNW] = f_TNW + feq[TNW];
- f[TSE] = f_TSE + feq[TSE];
- f[TSW] = f_TSW + feq[TSW];
- f[BNE] = f_TSW + feq[BNE];
- f[BNW] = f_TSE + feq[BNW];
- f[BSE] = f_TNW + feq[BSE];
- f[BSW] = f_TNE + feq[BSW];
- f[ZERO] = f_ZERO + feq[ZERO];
+ //int x1p, x2p, x3p;
+
+ ////x-derivative
+ //if (x1+1 >= maxMX1) x1p = x1;
+ //else x1p = x1+1;
+ //double vx1Plusx1 = matrix(Vx1, x1p, x2, x3);
+ //double vx2Plusx1 = matrix(Vx2, x1p, x2, x3);
+ //double vx3Plusx1 = matrix(Vx3, x1p, x2, x3);
+
+ //if (x1-1 < minX1) x1p = x1;
+ //else x1p = x1-1;
+ //double vx1Minusx1 = matrix(Vx1, x1p, x2, x3);
+ //double vx2Minusx1 = matrix(Vx2, x1p, x2, x3);
+ //double vx3Minusx1 = matrix(Vx3, x1p, x2, x3);
+
+ ////y-derivative
+ //if (x2+1 >= maxMX2) x2p = x2;
+ //else x2p = x2+1;
+ //double vx1Plusx2 = matrix(Vx1, x1, x2p, x3);
+ //double vx2Plusx2 = matrix(Vx2, x1, x2p, x3);
+ //double vx3Plusx2 = matrix(Vx3, x1, x2p, x3);
+
+ //if (x2-1 < minX2) x2p = x2;
+ //else x2p = x2-1;
+ //double vx1Minusx2 = matrix(Vx1, x1, x2p, x3);
+ //double vx2Minusx2 = matrix(Vx2, x1, x2p, x3);
+ //double vx3Minusx2 = matrix(Vx3, x1, x2p, x3);
+
+ ////z-derivative
+ //if (x3+1 >= maxMX3) x3p = x3;
+ //else x3p = x3+1;
+ //double vx1Plusx3 = matrix(Vx1, x1, x2, x3p);
+ //double vx2Plusx3 = matrix(Vx2, x1, x2, x3p);
+ //double vx3Plusx3 = matrix(Vx3, x1, x2, x3p);
+
+ //if (x3-1 < minX3) x3p = x3;
+ //else x3p = x3-1;
+ //double vx1Minusx3 = matrix(Vx1, x1, x2, x3);
+ //double vx2Minusx3 = matrix(Vx2, x1, x2, x3);
+ //double vx3Minusx3 = matrix(Vx3, x1, x2, x3);
+
+ //double ax = (vx1Plusx1 - vx1Minusx1) / (2.0);
+ //double bx = (vx2Plusx1 - vx2Minusx1) / (2.0);
+ //double cx = (vx3Plusx1 - vx3Minusx1) / (2.0);
+
+ //double ay = (vx1Plusx2 - vx1Minusx2) / (2.0);
+ //double by = (vx2Plusx2 - vx2Minusx2) / (2.0);
+ //double cy = (vx3Plusx2 - vx3Minusx2) / (2.0);
+
+ //double az = (vx1Plusx3 - vx1Minusx3) / (2.0);
+ //double bz = (vx2Plusx3 - vx2Minusx3) / (2.0);
+ //double cz = (vx3Plusx3 - vx3Minusx3) / (2.0);
+ //double eps_new = 1.0;
+ //LBMReal op = 1.;
+
+ //LBMReal feq[27];
+
+ //calcFeqsFct(feq, rho, vx1, vx2, vx3);
+
+ //double f_E = eps_new *((5.*ax*o + 5.*by*o + 5.*cz*o - 8.*ax*op + 4.*by*op + 4.*cz*op) / (54.*o*op));
+ //double f_N = f_E + eps_new *((2.*(ax - by)) / (9.*o));
+ //double f_T = f_E + eps_new *((2.*(ax - cz)) / (9.*o));
+ //double f_NE = eps_new *(-(5.*cz*o + 3.*(ay + bx)*op - 2.*cz*op + ax*(5.*o + op) + by*(5.*o + op)) / (54.*o*op));
+ //double f_SE = f_NE + eps_new *((ay + bx) / (9.*o));
+ //double f_TE = eps_new *(-(5.*cz*o + by*(5.*o - 2.*op) + 3.*(az + cx)*op + cz*op + ax*(5.*o + op)) / (54.*o*op));
+ //double f_BE = f_TE + eps_new *((az + cx) / (9.*o));
+ //double f_TN = eps_new *(-(5.*ax*o + 5.*by*o + 5.*cz*o - 2.*ax*op + by*op + 3.*bz*op + 3.*cy*op + cz*op) / (54.*o*op));
+ //double f_BN = f_TN + eps_new *((bz + cy) / (9.*o));
+ //double f_ZERO = eps_new *((5.*(ax + by + cz)) / (9.*op));
+ //double f_TNE = eps_new *(-(ay + az + bx + bz + cx + cy) / (72.*o));
+ //double f_TSW = -eps_new *((ay + bx) / (36.*o)) - f_TNE;
+ //double f_TSE = -eps_new *((az + cx) / (36.*o)) - f_TNE;
+ //double f_TNW = -eps_new *((bz + cy) / (36.*o)) - f_TNE;
+
+
+ //f[E] = f_E + feq[E];
+ //f[W] = f_E + feq[W];
+ //f[N] = f_N + feq[N];
+ //f[S] = f_N + feq[S];
+ //f[T] = f_T + feq[T];
+ //f[B] = f_T + feq[B];
+ //f[NE] = f_NE + feq[NE];
+ //f[SW] = f_NE + feq[SW];
+ //f[SE] = f_SE + feq[SE];
+ //f[NW] = f_SE + feq[NW];
+ //f[TE] = f_TE + feq[TE];
+ //f[BW] = f_TE + feq[BW];
+ //f[BE] = f_BE + feq[BE];
+ //f[TW] = f_BE + feq[TW];
+ //f[TN] = f_TN + feq[TN];
+ //f[BS] = f_TN + feq[BS];
+ //f[BN] = f_BN + feq[BN];
+ //f[TS] = f_BN + feq[TS];
+ //f[TNE] = f_TNE + feq[TNE];
+ //f[TNW] = f_TNW + feq[TNW];
+ //f[TSE] = f_TSE + feq[TSE];
+ //f[TSW] = f_TSW + feq[TSW];
+ //f[BNE] = f_TSW + feq[BNE];
+ //f[BNW] = f_TSE + feq[BNW];
+ //f[BSE] = f_TNW + feq[BSE];
+ //f[BSW] = f_TNE + feq[BSW];
+ //f[ZERO] = f_ZERO + feq[ZERO];
+
+ calcFeqsFct(f, rho, vx1, vx2, vx3);
distributions->setDistribution(f, ix1, ix2, ix3);
distributions->setDistributionInv(f, ix1, ix2, ix3);
+ boost::dynamic_pointer_cast<InitDensityLBMKernel>(kernel)->setVelocity(ix1, ix2, ix3, vx1, vx2, vx3);
}
+
}
}
//////////////////////////////////////////////////////////////////////////