diff --git a/source/Applications/aperm/aperm.cpp b/source/Applications/aperm/aperm.cpp
index beae1f1e04387b05c1d513ebfbf559a1cff9d525..327d63564f95f2c7120634b42c8d5e21f6cb8915 100644
--- a/source/Applications/aperm/aperm.cpp
+++ b/source/Applications/aperm/aperm.cpp
@@ -41,7 +41,10 @@ void run(string configname)
double timeSeriesOutTime = config.getDouble("timeSeriesOutTime");
bool logToFile = config.getBool("logToFile");
bool spongeLayer = config.getBool("spongeLayer");
-
+ vector<double> nupsStep = config.getVector<double>("nupsStep");
+ int numOfParts = config.getInt("numOfParts");
+ bool gridPrepare = config.getBool("gridPrepare");
+ double deltax = config.getDouble("deltax");
CommunicatorPtr comm = MPICommunicator::getInstance();
int myid = comm->getProcessID();
@@ -64,7 +67,7 @@ void run(string configname)
}
}
- //Sleep(30000);
+ Sleep(30000);
if (myid == 0) UBLOG(logINFO, "Testcase permeability");
@@ -86,17 +89,17 @@ void run(string configname)
const int baseLevel = 0;
double coord[6];
- double deltax;
+ //double deltax;
Grid3DPtr grid(new Grid3D(comm));
//////////////////////////////////////////////////////////////////////////
//restart
UbSchedulerPtr rSch(new UbScheduler(restartStep, restartStepStart));
- RestartCoProcessor rp(grid, rSch, comm, pathname, RestartCoProcessor::BINARY);
+ RestartCoProcessor rp(grid, rSch, comm, pathname, RestartCoProcessor::TXT);
//////////////////////////////////////////////////////////////////////////
- if (grid->getTimeStep() == 0)
+ if (gridPrepare)
{
if (myid == 0) UBLOG(logINFO, "new start..");
@@ -141,7 +144,15 @@ void run(string configname)
double g_maxX2 = sample->getX2Maximum();
double g_maxX3 = sample->getX3Maximum();
- deltax = (g_maxX3-g_minX3) /(nx3*blocknx3);
+ //////////////////////////////////////////////////////////////////////////
+ double nx3_temp = floor((g_maxX3-g_minX3)/(deltax*(double)blocknx3));
+
+ deltax = (g_maxX3-g_minX3)/(nx3_temp*(double)blocknx3);
+
+ //g_maxX3 -= 0.5* deltax;
+ //////////////////////////////////////////////////////////////////////////
+
+ //deltax = (g_maxX3-g_minX3) /(nx3*blocknx3);
double blockLength = (double)blocknx1*deltax;
@@ -219,6 +230,7 @@ void run(string configname)
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));
@@ -226,6 +238,7 @@ void run(string configname)
//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);
@@ -251,6 +264,7 @@ void run(string configname)
ppblocks->process(0);
ppblocks.reset();
+
unsigned long nob = grid->getNumberOfBlocks();
int gl = 3;
@@ -278,8 +292,8 @@ void run(string configname)
}
LBMKernel3DPtr kernel;
- kernel = LBMKernel3DPtr(new LBMKernelETD3Q27CCLB(blocknx1, blocknx2, blocknx3, LBMKernelETD3Q27CCLB::NORMAL));
-
+ //kernel = LBMKernel3DPtr(new LBMKernelETD3Q27CCLB(blocknx1, blocknx2, blocknx3, LBMKernelETD3Q27CCLB::NORMAL));
+ kernel = LBMKernel3DPtr(new VoidLBMKernel(blocknx1, blocknx2, blocknx3));
BCProcessorPtr bcProc(new D3Q27ETBCProcessor());
BoundaryConditionPtr densityBC(new NonEqDensityBoundaryCondition());
@@ -296,27 +310,27 @@ void run(string configname)
//BC
intHelper.setBC();
- BoundaryConditionBlockVisitor bcVisitor;
- grid->accept(bcVisitor);
+ //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);
+ //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);
+ //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();
+ //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();
@@ -326,39 +340,71 @@ void run(string configname)
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);
+ UbFileOutputASCII outf(pathname + "/checkpoints/coord.txt");
+ outf.writeDouble(deltax);
+ outf.writeDouble(coord[0]);
+ outf.writeDouble(coord[1]);
+ outf.writeDouble(coord[2]);
+ outf.writeDouble(coord[3]);
+ outf.writeDouble(coord[4]);
+ outf.writeDouble(coord[5]);
+ outf.writeDouble(g_minX1);
+ outf.writeDouble(g_maxX1);
+ outf.writeDouble(availMem);
+ outf.writeDouble(needMem);
////////////////////////////////////////////////////////
grid->addInteractor(inflowInt);
+
if (myid == 0) UBLOG(logINFO, "Preprozess - end");
+
+ if (comm->getNumberOfProcesses() == 1)
+ {
+ UBLOG(logINFO, "Prepare grid - start");
+ rp.writeDistributedGrid(grid, numOfParts);
+ UBLOG(logINFO, "Prepare grid - end");
+ return;
+ }
}
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);
+ UbFileInputASCII inf(pathname+"/checkpoints/coord.txt");
+ deltax = inf.readDouble();
+ coord[0] = inf.readDouble();
+ coord[1] = inf.readDouble();
+ coord[2] = inf.readDouble();
+ coord[3] = inf.readDouble();
+ coord[4] = inf.readDouble();
+ coord[5] = inf.readDouble();
+ double g_minX1 = inf.readDouble();
+ double g_maxX1 = inf.readDouble();
+ double availMem = inf.readDouble();
+ double needMem = inf.readDouble();
////////////////////////////////////////////////////////
+ LBMKernel3DPtr kernel = LBMKernel3DPtr(new LBMKernelETD3Q27CCLB(blocknx1, blocknx2, blocknx3, LBMKernelETD3Q27CCLB::NORMAL));
+ BCProcessorPtr bcProc(new D3Q27ETBCProcessor());
+ kernel->setBCProcessor(bcProc);
+ SetKernelBlockVisitor kernelVisitor(kernel, nu_LB, availMem, needMem, SetKernelBlockVisitor::ChangeBC);
+ grid->accept(kernelVisitor);
+
+ 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);
+
//new nu
//ViscosityBlockVisitor nuVisitor(nu_LB);
//grid->accept(nuVisitor);
@@ -380,17 +426,21 @@ void run(string configname)
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);
+ //Postrozess
+ UbSchedulerPtr geoSch(new UbScheduler(1));
+ MacroscopicQuantitiesCoProcessorPtr ppgeo(
+ new MacroscopicQuantitiesCoProcessor(grid, geoSch, pathname, WbWriterVtkXmlBinary::getInstance(), conv, true));
+ ppgeo->process(0);
+ ppgeo.reset();
+
if (myid == 0) UBLOG(logINFO, "Restart - end");
}
- UbSchedulerPtr nupsSch(new UbScheduler(10, 30, 100));
+ UbSchedulerPtr nupsSch(new UbScheduler(nupsStep[0], nupsStep[1], nupsStep[2]));
//nupsSch->addSchedule(nupsStep[0], nupsStep[1], nupsStep[2]);
NUPSCounterCoProcessor npr(grid, nupsSch, numOfThreads, comm);
diff --git a/source/Applications/aperm/configBombadil2.txt b/source/Applications/aperm/configBombadil2.txt
index 9a29f30c232af059660accf738a4259403d9d860..9528e9aa9c0ab4ae425931e4c903327fdf3c2c97 100644
--- a/source/Applications/aperm/configBombadil2.txt
+++ b/source/Applications/aperm/configBombadil2.txt
@@ -1,6 +1,6 @@
pathname = d:/temp/aperm
pathGeo = d:/Projects/SFB880/GeometrienPoroeseMedien/isotrop/PA80-110
-numOfThreads = 4
+numOfThreads = 1
availMem = 3e9
logToFile = false
@@ -25,7 +25,8 @@ pmL3 = 0.786e-3
#grid
blocknx = 30
nx3 = 5
-spongeLayer=true
+deltax = 10e-6
+spongeLayer=false
dp_LB = 0.001
#nu_LB = 0.168666666667
@@ -37,8 +38,13 @@ nu_LB = 0.168666666667e-4
timeSeriesFile = /timeseries/simAlu80_5
timeSeriesOutTime = 10
-restartStep = 20000
-restartStepStart=20000
+restartStep = 100
+restartStepStart=100
-endTime = 60000
-outTime = 10
\ No newline at end of file
+endTime = 100
+outTime = 10
+
+nupsStep = 10 30 100
+
+gridPrepare = true
+numOfParts = 4
\ No newline at end of file
diff --git a/source/Applications/pChannel/pChannel.cpp.hlrn b/source/Applications/pChannel/pChannel.cpp.hlrn
new file mode 100644
index 0000000000000000000000000000000000000000..04274a4a24fdb42156d54e56baff7560a323df18
--- /dev/null
+++ b/source/Applications/pChannel/pChannel.cpp.hlrn
@@ -0,0 +1,729 @@
+#include <iostream>
+#include <string>
+#include "VirtualFluids.h"
+
+
+
+//#include <thread>
+
+using namespace std;
+
+//////////////////////////////////////////////////////////////////////////
+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");
+ vector<int> pmNX = config.getVector<int>("pmNX");
+ double lthreshold = config.getDouble("lthreshold");
+ double uthreshold = config.getDouble("uthreshold");
+ vector<float> voxelDeltaX = config.getVector<float>("voxelDeltaX");
+ vector<int> blocknx = config.getVector<int>("blocknx");
+ double u_LB = config.getDouble("u_LB");
+ double restartStep = config.getDouble("restartStep");
+ double restartStepStart = config.getDouble("restartStepStart");
+ double endTime = config.getDouble("endTime");
+ double outTime = config.getDouble("outTime");
+ double availMem = config.getDouble("availMem");
+ bool rawFile = config.getBool("rawFile");
+ bool logToFile = config.getBool("logToFile");
+ bool writeSample = config.getBool("writeSample");
+ vector<double> pmL = config.getVector<double>("pmL");
+ double deltaXfine = config.getDouble("deltaXfine");
+ int refineLevel = config.getInt("refineLevel");
+ bool thinWall = config.getBool("thinWall");
+ double Re = config.getDouble("Re");
+ double channelHigh = config.getDouble("channelHigh");
+ double lengthFactor = config.getDouble("lengthFactor");
+ double forcing = config.getDouble("forcing");
+ bool changeQs = config.getBool("changeQs");
+ double timeAvStart = config.getDouble("timeAvStart");
+ double timeAvStop = config.getDouble("timeAvStop");
+ bool averaging = config.getBool("averaging");
+ bool averagingReset = config.getBool("averagingReset");
+ double nupsSteps = config.getDouble("nupsSteps");
+ double timeLineTsStart = config.getDouble("timeLineTsStart");
+ double timeLineTsStop = config.getDouble("timeLineTsStop");
+
+
+ CommunicatorPtr comm = MPICommunicator::getInstance();
+ int myid = comm->getProcessID();
+
+ if (logToFile)
+ {
+#if defined(__unix__)
+ if (myid == 0)
+ {
+ const char* str = pathname.c_str();
+ 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 porous channel");
+
+ LBMReal rho_LB = 0.0;
+
+ LBMUnitConverterPtr conv = LBMUnitConverterPtr(new LBMUnitConverter());
+
+ const int baseLevel = 0;
+ double deltaXcoarse = deltaXfine*(double)(1<<refineLevel);
+
+ double coord[6];
+ bool restart;
+
+ vector<double> origin(3);
+ origin[0] = 0;
+ origin[1] = 0;
+ origin[2] = 0;
+
+ //real velocity is 49.63 m/s
+ double nu_LB;
+
+ 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...");
+ restart = false;
+
+ double offsetMinX3 = pmL[2];
+
+ double offsetMaxX1 = pmL[0]*lengthFactor;
+ double offsetMaxX2 = pmL[1]*2.0;
+ double offsetMaxX3 = pmL[2] + channelHigh; //DLR-F15 //pmL[2]*2.0;
+
+ //bounding box
+ double g_minX1 = origin[0];
+ double g_minX2 = origin[1];
+ double g_minX3 = origin[2];
+
+ double g_maxX1 = origin[0] + offsetMaxX1;
+ double g_maxX2 = origin[1] + offsetMaxX2;
+ double g_maxX3 = origin[2] + offsetMaxX3;
+//////////////////////////////////////////////////////////////////////////
+ double nx1_temp = floor((g_maxX1-g_minX1) /(deltaXcoarse*(double)blocknx[0]));
+
+ deltaXcoarse = (g_maxX1-g_minX1) /(nx1_temp*(double)blocknx[0]);
+
+ g_maxX1 -= 0.5* deltaXcoarse;
+//////////////////////////////////////////////////////////////////////////
+ double blockLength = (double)blocknx[0]*deltaXcoarse;
+
+ grid->setPeriodicX1(true);
+ grid->setPeriodicX2(true);
+ grid->setPeriodicX3(false);
+ grid->setDeltaX(deltaXcoarse);
+ grid->setBlockNX(blocknx[0], blocknx[1], blocknx[2]);
+
+ GbObject3DPtr gridCube(new GbCuboid3D(g_minX1, g_minX2, g_minX3, g_maxX1, g_maxX2, g_maxX3));
+ if (myid == 0) GbSystem3D::writeGeoObject(gridCube.get(), pathname + "/geo/gridCube", WbWriterVtkXmlBinary::getInstance());
+
+
+ GenBlocksGridVisitor genBlocks(gridCube);
+ grid->accept(genBlocks);
+ double channel_high = channelHigh; // g_maxX3-g_minX3;
+ double channel_high_LB = channel_high/deltaXcoarse;
+//////////////////////////////////////////////////////////////////////////
+ //nu_LB = 0.005;
+ nu_LB = (u_LB*channel_high_LB)/Re;
+//////////////////////////////////////////////////////////////////////////
+ if (myid == 0)
+ {
+ UBLOG(logINFO, "Parameters:");
+ UBLOG(logINFO, "Re = " << Re);
+ UBLOG(logINFO, "u_LB = " << u_LB);
+ UBLOG(logINFO, "rho_LB = " << rho_LB);
+ UBLOG(logINFO, "nu_LB = " << nu_LB);
+ UBLOG(logINFO, "dx coarse = " << deltaXcoarse << " m");
+ UBLOG(logINFO, "dx fine = " << grid->getDeltaX(refineLevel) << " m");
+ UBLOG(logINFO, "number of levels = " << refineLevel + 1);
+ UBLOG(logINFO, "number of processes = " << comm->getNumberOfProcesses());
+ UBLOG(logINFO, "number of threads = " << numOfThreads);
+ UBLOG(logINFO, "path = " << pathname);
+ UBLOG(logINFO, "Preprocess - start");
+ }
+
+ //////////////////////////////////////////////////////////////////////////
+ //refinement
+ double blockLengthX3Fine = grid->getDeltaX(refineLevel) * blocknx[2];
+
+ GbCuboid3DPtr refineBoxTop(new GbCuboid3D(g_minX1-blockLength, g_minX2-blockLength, g_maxX3-blockLengthX3Fine, g_maxX1+blockLength, g_maxX2+blockLength, g_maxX3+blockLength));
+ if (myid == 0) GbSystem3D::writeGeoObject(refineBoxTop.get(), pathname + "/geo/refineBoxTop", WbWriterVtkXmlASCII::getInstance());
+
+ GbCuboid3DPtr refineBoxBottom(new GbCuboid3D(g_minX1-blockLength, g_minX2-blockLength, g_minX3, g_maxX1+blockLength, g_maxX2+blockLength, g_minX3+offsetMinX3+blockLengthX3Fine));
+ //GbCuboid3DPtr refineBoxBottom(new GbCuboid3D(g_minX1-blockLength, g_minX2-blockLength, g_minX3-blockLengthX3Fine, g_maxX1+blockLength, g_maxX2+blockLength, g_minX3+blockLengthX3Fine));
+ if (myid == 0) GbSystem3D::writeGeoObject(refineBoxBottom.get(), pathname + "/geo/refineBoxBottom", WbWriterVtkXmlASCII::getInstance());
+
+ if (refineLevel > 0)
+ {
+ if (myid == 0) UBLOG(logINFO, "Refinement - start");
+ RefineCrossAndInsideGbObjectHelper refineHelper(grid, refineLevel);
+ refineHelper.addGbObject(refineBoxTop, refineLevel);
+ refineHelper.addGbObject(refineBoxBottom, refineLevel);
+ refineHelper.refine();
+ if (myid == 0) UBLOG(logINFO, "Refinement - end");
+ }
+ //////////////////////////////////////////////////////////////////////////
+
+ //walls
+ GbCuboid3DPtr addWallZmin(new GbCuboid3D(g_minX1-blockLength, g_minX2-blockLength, g_minX3-blockLength, g_maxX1+blockLength, g_maxX2+blockLength, g_minX3));
+ if (myid == 0) GbSystem3D::writeGeoObject(addWallZmin.get(), pathname+"/geo/addWallZmin", WbWriterVtkXmlASCII::getInstance());
+
+ 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());
+
+
+ //wall interactors
+ int bbOption = 1;
+ D3Q27BoundaryConditionAdapterPtr bcNoSlip(new D3Q27NoSlipBCAdapter(bbOption));
+ D3Q27InteractorPtr addWallZminInt(new D3Q27Interactor(addWallZmin, grid, bcNoSlip, Interactor3D::SOLID));
+ D3Q27InteractorPtr addWallZmaxInt(new D3Q27Interactor(addWallZmax, grid, bcNoSlip, Interactor3D::SOLID));
+
+ ////////////////////////////////////////////////
+ //TEST
+ //GbObject3DPtr testCube(new GbCuboid3D(g_minX1 + 2.0 * blockLength, g_minX2 + 2.0 * blockLength, g_minX3 + 5.0 * blockLength,
+ // g_minX1 + 3.0 * blockLength, g_minX2 + 3.0 * blockLength, g_minX3 + 6.0 * blockLength));
+ //if (myid == 0) GbSystem3D::writeGeoObject(testCube.get(), pathname + "/geo/testCube", WbWriterVtkXmlBinary::getInstance());
+ //D3Q27InteractorPtr testCubeInt(new D3Q27Interactor(testCube, grid, bcNoSlip, Interactor3D::SOLID));
+ ///////////////////////////////////////////////
+
+ ////////////////////////////////////////////
+ //METIS
+ Grid3DVisitorPtr metisVisitor(new MetisPartitioningGridVisitor(comm, MetisPartitioningGridVisitor::LevelBased, D3Q27System::BSW, MetisPartitioner::KWAY));
+ ////////////////////////////////////////////
+ //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(addWallZminInt);
+ intHelper.addInteractor(addWallZmaxInt);
+ //////////////////////////////////////////////////////////////////////////
+ //TEST
+ //intHelper.addInteractor(testCubeInt);
+ //////////////////////////////////////////////////////////////////////////
+ intHelper.selectBlocks();
+ if (myid == 0) UBLOG(logINFO, "deleteSolidBlocks - end");
+ //////////////////////////////////////
+
+ WriteBlocksCoProcessorPtr ppblocks(new WriteBlocksCoProcessor(grid, UbSchedulerPtr(new UbScheduler(1)), pathname, WbWriterVtkXmlBinary::getInstance(), comm));
+ ppblocks->process(0);
+ ppblocks.reset();
+
+ unsigned long long numberOfBlocks = (unsigned long long)grid->getNumberOfBlocks();
+ int ghostLayer = 3;
+ unsigned long long numberOfNodesPerBlock = (unsigned long long)(blocknx[0])* (unsigned long long)(blocknx[1])* (unsigned long long)(blocknx[2]);
+ unsigned long long numberOfNodes = numberOfBlocks * numberOfNodesPerBlock;
+ unsigned long long numberOfNodesPerBlockWithGhostLayer = numberOfBlocks * (blocknx[0] + ghostLayer) * (blocknx[1] + ghostLayer) * (blocknx[2] + ghostLayer);
+ double needMemAll = double(numberOfNodesPerBlockWithGhostLayer*(27 * sizeof(double) + sizeof(int) + sizeof(float) * 4));
+ double needMem = needMemAll / double(comm->getNumberOfProcesses());
+
+
+ if (myid == 0)
+ {
+ UBLOG(logINFO, "Number of blocks = " << numberOfBlocks);
+ UBLOG(logINFO, "Number of nodes = " << numberOfNodes);
+ 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*numberOfNodesPerBlock);
+ }
+ UBLOG(logINFO, "Necessary memory = " << needMemAll << " bytes");
+ UBLOG(logINFO, "Necessary memory per process = " << needMem << " bytes");
+ UBLOG(logINFO, "Available memory per process = " << availMem << " bytes");
+ }
+
+ LBMKernel3DPtr kernel = LBMKernel3DPtr(new LBMKernelETD3Q27CCLB(blocknx[0], blocknx[1], blocknx[2], LBMKernelETD3Q27CCLB::NORMAL));
+
+ mu::Parser fctForcingX1;
+ fctForcingX1.SetExpr("Fx1");
+ fctForcingX1.DefineConst("Fx1", 1.0e-6);
+
+ kernel->setWithForcing(true);
+
+ BCProcessorPtr bcProc;
+ BoundaryConditionPtr noSlipBC;
+
+ if (thinWall)
+ {
+ bcProc = BCProcessorPtr(new D3Q27ETForThinWallBCProcessor());
+ noSlipBC = BoundaryConditionPtr(new ThinWallNoSlipBoundaryCondition());
+ }
+ else
+ {
+ bcProc = BCProcessorPtr(new D3Q27ETBCProcessor());
+ noSlipBC = BoundaryConditionPtr(new NoSlipBoundaryCondition());
+ }
+
+ bcProc->addBC(noSlipBC);
+
+ kernel->setBCProcessor(bcProc);
+
+ SetKernelBlockVisitor kernelVisitor(kernel, nu_LB, availMem, needMem);
+ grid->accept(kernelVisitor);
+
+ //////////////////////////////////
+ //undef nodes for refinement
+ if (refineLevel > 0)
+ {
+ D3Q27SetUndefinedNodesBlockVisitor undefNodesVisitor;
+ grid->accept(undefNodesVisitor);
+ }
+
+
+ //BC
+ intHelper.setBC();
+
+ ////porous media
+ {
+ string samplePathname = pathGeo + "/" + sampleFilename;
+
+ GbVoxelMatrix3DPtr sample(new GbVoxelMatrix3D(pmNX[0], pmNX[1], pmNX[2], 0, lthreshold, uthreshold));
+ if (rawFile)
+ {
+ sample->readBufferedMatrixFromRawFile<unsigned short>(samplePathname, GbVoxelMatrix3D::BigEndian);
+ }
+ else
+ {
+ sample->readMatrixFromVtiASCIIFile(samplePathname);
+ }
+
+ sample->setVoxelMatrixDelta(voxelDeltaX[0], voxelDeltaX[1], voxelDeltaX[2]);
+ sample->setVoxelMatrixMininum(origin[0], origin[1], origin[2]);
+
+ int bounceBackOption = 1;
+ bool vxFile = false;
+ int i = 0;
+ for (int x = 0; x < lengthFactor; x+=2)
+ {
+ double offset = pmL[0] * (double)x;
+ //sample 0
+ if (myid == 0) UBLOG(logINFO, "sample # " << i);
+ sample->setVoxelMatrixMininum(origin[0]+offset, origin[1], origin[2]);
+ Utilities::voxelMatrixDiscretisation(sample, pathname, myid, i, grid, bounceBackOption, vxFile);
+ i++;
+
+ if (myid == 0)
+ {
+ UBLOG(logINFO, "PID = " << myid << " Total Physical Memory (RAM): " << Utilities::getTotalPhysMem());
+ UBLOG(logINFO, "PID = " << myid << " Physical Memory currently used: " << Utilities::getPhysMemUsed());
+ UBLOG(logINFO, "PID = " << myid << " Physical Memory currently used by current process: " << Utilities::getPhysMemUsedByMe());
+ }
+
+ //sample 1
+ if (myid == 0) UBLOG(logINFO, "sample # " << i);
+ sample->setVoxelMatrixMininum(origin[0]+pmL[0]+offset, origin[1], origin[2]);
+ sample->mirrorX();
+ Utilities::voxelMatrixDiscretisation(sample, pathname, myid, i, grid, bounceBackOption, vxFile);
+ i++;
+
+ if (myid == 0)
+ {
+ UBLOG(logINFO, "PID = " << myid << " Total Physical Memory (RAM): " << Utilities::getTotalPhysMem());
+ UBLOG(logINFO, "PID = " << myid << " Physical Memory currently used: " << Utilities::getPhysMemUsed());
+ UBLOG(logINFO, "PID = " << myid << " Physical Memory currently used by current process: " << Utilities::getPhysMemUsedByMe());
+ }
+
+ //sample 2
+ if (myid == 0) UBLOG(logINFO, "sample # " << i);
+ sample->setVoxelMatrixMininum(origin[0]+pmL[0]+offset, origin[1]+pmL[1], origin[2]);
+ sample->mirrorY();
+ Utilities::voxelMatrixDiscretisation(sample, pathname, myid, i, grid, bounceBackOption, vxFile);
+ i++;
+
+ if (myid == 0)
+ {
+ UBLOG(logINFO, "PID = " << myid << " Total Physical Memory (RAM): " << Utilities::getTotalPhysMem());
+ UBLOG(logINFO, "PID = " << myid << " Physical Memory currently used: " << Utilities::getPhysMemUsed());
+ UBLOG(logINFO, "PID = " << myid << " Physical Memory currently used by current process: " << Utilities::getPhysMemUsedByMe());
+ }
+
+ //sample 3
+ if (myid == 0) UBLOG(logINFO, "sample # " << i);
+ sample->setVoxelMatrixMininum(origin[0]+offset, origin[1]+pmL[1], origin[2]);
+ sample->mirrorX();
+ Utilities::voxelMatrixDiscretisation(sample, pathname, myid, i, grid, bounceBackOption, vxFile);
+ sample->mirrorY();
+ i++;
+ }
+
+ }
+ BoundaryConditionBlockVisitor bcVisitor;
+ grid->accept(bcVisitor);
+
+ mu::Parser inflowProfile;
+ inflowProfile.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);
+ inflowProfile.DefineConst("h", pmL[2]);
+
+ D3Q27ETInitDistributionsBlockVisitor initVisitor(nu_LB, rho_LB);
+ initVisitor.setVx1(inflowProfile);
+ //initVisitor.setVx1(u_LB);
+ grid->accept(initVisitor);
+
+ ////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);
+
+ //Postrozess
+ UbSchedulerPtr geoSch(new UbScheduler(1));
+ MacroscopicQuantitiesCoProcessorPtr ppgeo(
+ new MacroscopicQuantitiesCoProcessor(grid, geoSch, pathname, WbWriterVtkXmlBinary::getInstance(), conv, true));
+ ppgeo->process(0);
+ ppgeo.reset();
+
+ coord[0] = g_minX1;
+ coord[1] = g_minX2;
+ coord[2] = g_minX3;
+ coord[3] = g_maxX1;
+ coord[4] = g_maxX2;
+ coord[5] = g_maxX3;
+
+ ////////////////////////////////////////////////////////
+ FILE * pFile;
+ string str = pathname + "/checkpoints/coord.txt";
+ pFile = fopen(str.c_str(), "w");
+ fprintf(pFile, "%g\n", deltaXcoarse);
+ fprintf(pFile, "%g\n", nu_LB);
+ fprintf(pFile, "%g\n", coord[0]);
+ fprintf(pFile, "%g\n", coord[1]);
+ fprintf(pFile, "%g\n", coord[2]);
+ fprintf(pFile, "%g\n", coord[3]);
+ fprintf(pFile, "%g\n", coord[4]);
+ fprintf(pFile, "%g\n", coord[5]);
+ fclose(pFile);
+ ////////////////////////////////////////////////////////
+
+ if (myid == 0) UBLOG(logINFO, "Preprozess - end");
+ }
+ else
+ {
+ ////////////////////////////////////////////////////////
+ FILE * pFile;
+ string str = pathname + "/checkpoints/coord.txt";
+ pFile = fopen(str.c_str(), "r");
+ fscanf(pFile, "%lg\n", &deltaXcoarse);
+ fscanf(pFile, "%lg\n", &nu_LB);
+ fscanf(pFile, "%lg\n", &coord[0]);
+ fscanf(pFile, "%lg\n", &coord[1]);
+ fscanf(pFile, "%lg\n", &coord[2]);
+ fscanf(pFile, "%lg\n", &coord[3]);
+ fscanf(pFile, "%lg\n", &coord[4]);
+ fscanf(pFile, "%lg\n", &coord[5]);
+ fclose(pFile);
+ ////////////////////////////////////////////////////////
+
+ if (myid == 0)
+ {
+ UBLOG(logINFO, "Parameters:");
+ UBLOG(logINFO, "Re = " << Re);
+ UBLOG(logINFO, "u_LB = " << u_LB);
+ UBLOG(logINFO, "rho_LB = " << rho_LB);
+ UBLOG(logINFO, "nu_LB = " << nu_LB);
+ UBLOG(logINFO, "dx coarse = " << deltaXcoarse << " m");
+ UBLOG(logINFO, "dx fine = " << grid->getDeltaX(refineLevel) << " m");
+ UBLOG(logINFO, "number of levels = " << refineLevel + 1);
+ UBLOG(logINFO, "numOfThreads = " << numOfThreads);
+ UBLOG(logINFO, "path = " << pathname);
+ }
+
+ /////////////////////////////////////////////////////////////
+ /////////////////////////////////////////////////////////////
+ //bounding box
+// double offsetMinX3 = pmL[2];
+//
+// double offsetMaxX1 = pmL[0]*lengthFactor;
+// double offsetMaxX2 = pmL[1]*2.0;
+// double offsetMaxX3 = channelHigh;
+//
+// double g_minX1 = origin[0];
+// double g_minX2 = origin[1];
+// double g_minX3 = origin[2];
+//
+// double g_maxX1 = origin[0]+offsetMaxX1;
+// double g_maxX2 = origin[1]+offsetMaxX2;
+// double g_maxX3 = origin[2]+offsetMaxX3;
+//
+// double blockLength = (double)blocknx[0]*deltaXcoarse;
+//
+// GbCuboid3DPtr addWallZmin(new GbCuboid3D(g_minX1-blockLength, g_minX2-blockLength, g_minX3-blockLength, g_maxX1+blockLength, g_maxX2+blockLength, g_minX3+offsetMinX3));
+// if (myid==0) GbSystem3D::writeGeoObject(addWallZmin.get(), pathname+"/geo/addWallZmin", WbWriterVtkXmlASCII::getInstance());
+// int bbOption = 1;
+// D3Q27BoundaryConditionAdapterPtr bcNoSlip(new D3Q27NoSlipBCAdapter(bbOption));
+// D3Q27InteractorPtr addWallZminInt(new D3Q27Interactor(addWallZmin, grid, bcNoSlip, Interactor3D::SOLID));
+//
+// SetSolidOrTransBlockVisitor v1(addWallZminInt, SetSolidOrTransBlockVisitor::SOLID);
+// grid->accept(v1);
+// SetSolidOrTransBlockVisitor v2(addWallZminInt, SetSolidOrTransBlockVisitor::TRANS);
+// grid->accept(v2);
+//
+// std::vector<Block3DPtr> blocks;
+// std::vector<Block3DPtr>& sb = addWallZminInt->getSolidBlockSet();
+// if (myid==0) UBLOG(logINFO, "number of solid blocks = "<<sb.size());
+// blocks.insert(blocks.end(), sb.begin(), sb.end());
+// std::vector<Block3DPtr>& tb = addWallZminInt->getTransBlockSet();
+// if (myid==0) UBLOG(logINFO, "number of trans blocks = "<<tb.size());
+// blocks.insert(blocks.end(), tb.begin(), tb.end());
+//
+// if (myid==0) UBLOG(logINFO, "number of blocks = "<<blocks.size());
+//
+// BOOST_FOREACH(Block3DPtr block, blocks)
+// {
+// block->setActive(true);
+// addWallZminInt->setDifferencesToGbObject3D(block);
+// }
+//
+// //////////////////////////////////////////////
+// ////METIS
+// //Grid3DVisitorPtr metisVisitor(new MetisPartitioningGridVisitor(comm, MetisPartitioningGridVisitor::LevelBased, D3Q27System::BSW, MetisPartitioner::KWAY));
+// //////////////////////////////////////////////
+// ///////delete solid blocks
+// //if (myid==0) UBLOG(logINFO, "deleteSolidBlocks - start");
+// //InteractorsHelper intHelper(grid, metisVisitor);
+// //intHelper.addInteractor(addWallZminInt);
+// //intHelper.selectBlocks();
+// //if (myid==0) UBLOG(logINFO, "deleteSolidBlocks - end");
+// ////////////////////////////////////////
+// //intHelper.setBC();
+// //////////////////////////////////////////////////////////////
+//
+// BoundaryConditionBlockVisitor bcVisitor;
+// grid->accept(bcVisitor);
+//
+// WriteBlocksCoProcessorPtr ppblocks(new WriteBlocksCoProcessor(grid, UbSchedulerPtr(new UbScheduler(1)), pathname, WbWriterVtkXmlBinary::getInstance(), comm));
+// ppblocks->process(0);
+// ppblocks.reset();
+//
+// UbSchedulerPtr geoSch(new UbScheduler(1));
+// MacroscopicQuantitiesCoProcessorPtr ppgeo(
+// new MacroscopicQuantitiesCoProcessor(grid, geoSch, pathname, WbWriterVtkXmlBinary::getInstance(), conv, true));
+// ppgeo->process(0);
+// ppgeo.reset();
+ ////////////////////////////////////////////////////////////////
+
+ //set connectors
+ D3Q27InterpolationProcessorPtr iProcessor(new D3Q27IncompressibleOffsetInterpolationProcessor());
+ D3Q27SetConnectorsBlockVisitor setConnsVisitor(comm, true, D3Q27System::ENDDIR, nu_LB, iProcessor);
+ grid->accept(setConnsVisitor);
+
+ restart = true;
+
+ if (myid == 0) UBLOG(logINFO, "Restart - end");
+ }
+ UbSchedulerPtr nupsSch(new UbScheduler(nupsSteps));
+ NUPSCounterCoProcessor npr(grid, nupsSch, numOfThreads, comm);
+
+ UbSchedulerPtr stepSch(new UbScheduler(outTime));
+
+ MacroscopicQuantitiesCoProcessor pp(grid, stepSch, pathname, WbWriterVtkXmlBinary::getInstance(), conv);
+
+ double startStep = grid->getTimeStep();
+
+ //UbSchedulerPtr visSch(new UbScheduler());
+ //visSch->addSchedule(40000,40000,40000000);
+ //UbSchedulerPtr resSchRMS(new UbScheduler());
+ //resSchRMS->addSchedule(40000, startStep, 40000000);
+ //UbSchedulerPtr resSchMeans(new UbScheduler());
+ //resSchMeans->addSchedule(40000, startStep, 40000000);
+ //UbSchedulerPtr stepAvSch(new UbScheduler());
+ //stepAvSch->addSchedule(100, 0, 10000000);
+ //AverageValuesPostprocessor Avpp(grid, pathname, WbWriterVtkXmlBinary::getInstance(),
+ // stepSch/*wann wird rausgeschrieben*/, stepAvSch/*wann wird gemittelt*/, resSchMeans, resSchRMS/*wann wird resettet*/, restart);
+
+
+ UbSchedulerPtr AdjForcSch(new UbScheduler());
+ AdjForcSch->addSchedule(10, 0, 10000000);
+ D3Q27IntegrateValuesHelperPtr intValHelp(new D3Q27IntegrateValuesHelper(grid, comm,
+ coord[0], coord[1], coord[2],
+ coord[3], coord[4], coord[5]));
+ if (myid == 0) GbSystem3D::writeGeoObject(intValHelp->getBoundingBox().get(), pathname + "/geo/IntValHelp", WbWriterVtkXmlBinary::getInstance());
+
+ double vxTarget=u_LB;
+ AdjustForcingCoProcessor AdjForcPPPtr(grid, AdjForcSch, pathname, intValHelp, vxTarget, comm);
+
+ //mu::Parser decrViscFunc;
+ //decrViscFunc.SetExpr("nue0+c0/(t+1)/(t+1)");
+ //decrViscFunc.DefineConst("nue0", nu_LB*4.0);
+ //decrViscFunc.DefineConst("c0", 0.1);
+ //UbSchedulerPtr DecrViscSch(new UbScheduler());
+ //DecrViscSch->addSchedule(10, 0, 1000);
+ //DecreaseViscosityPostprocessor decrViscPPPtr(grid, DecrViscSch, &decrViscFunc, comm);
+
+ //if (changeQs)
+ //{
+ // double z1 = pmL[2];
+ // D3Q27IntegrateValuesHelperPtr intValHelp2(new D3Q27IntegrateValuesHelper(grid, comm,
+ // coord[0], coord[1], z1 - deltaXfine,
+ // coord[3], coord[4], z1 + deltaXfine));
+ // if (myid == 0) GbSystem3D::writeGeoObject(intValHelp2->getBoundingBox().get(), pathname + "/geo/intValHelp2", WbWriterVtkXmlBinary::getInstance());
+ // Utilities::ChangeRandomQs(intValHelp2);
+ //}
+
+ std::vector<double> levelCoords;
+ std::vector<int> levels;
+ std::vector<double> bounds;
+ bounds.push_back(0);
+ bounds.push_back(0);
+ bounds.push_back(1e-3);
+ bounds.push_back(0.064);
+ bounds.push_back(0.008);
+ bounds.push_back(0.018);
+ levels.push_back(3);
+ levels.push_back(2);
+ levels.push_back(1);
+ levels.push_back(0);
+ levels.push_back(1);
+ levels.push_back(2);
+ levels.push_back(3);
+ levelCoords.push_back(0);
+ levelCoords.push_back(0.0016-6.0*deltaXfine);
+ levelCoords.push_back(0.0016);
+ levelCoords.push_back(0.0024-6.0*deltaXfine*2.0);
+ levelCoords.push_back(0.0024);
+ levelCoords.push_back(0.0048-6.0*deltaXfine*4.0);
+ levelCoords.push_back(0.0048);
+ levelCoords.push_back(0.0144);
+ levelCoords.push_back(0.0144+6.0*deltaXfine*4.0);
+ levelCoords.push_back(0.0168);
+ levelCoords.push_back(0.0168+6.0*deltaXfine*2.0);
+ levelCoords.push_back(0.0176);
+ levelCoords.push_back(0.0176+6.0*deltaXfine);
+ levelCoords.push_back(0.018);
+ UbSchedulerPtr tavSch(new UbScheduler(1, timeAvStart, timeAvStop));
+ TimeAveragedValuesCoProcessorPtr tav(new TimeAveragedValuesCoProcessor(grid, pathname, WbWriterVtkXmlBinary::getInstance(), tavSch, comm,
+ TimeAveragedValuesCoProcessor::Velocity | TimeAveragedValuesCoProcessor::Fluctuations | TimeAveragedValuesCoProcessor::Triplecorrelations//));
+ ,levels, levelCoords, bounds, false));
+ if (averagingReset)
+ {
+ tav->reset();
+ }
+ //UbSchedulerPtr catalystSch(new UbScheduler(1));
+ //InSituCatalystCoProcessor catalyst(grid, catalystSch, "pchannel.py");
+
+ UbSchedulerPtr exitSch(new UbScheduler(10));
+ EmergencyExitCoProcessor exitCoProc(grid, exitSch, pathname, RestartCoProcessorPtr(&rp), comm);
+
+ //create line time series
+ UbSchedulerPtr tpcSch(new UbScheduler(1,timeLineTsStart,timeLineTsStop));
+
+ GbLine3DPtr line1(new GbLine3D(new GbPoint3D(0.0,0.004,0.00078),new GbPoint3D(0.064,0.004,0.00078)));
+ LineTimeSeriesCoProcessor lineTs1(grid, tpcSch,pathname+"/TimeSeries/line1.csv",line1, 3,comm);
+ if (myid==0) lineTs1.writeLine(pathname+"/geo/line1");
+
+ GbLine3DPtr line2(new GbLine3D(new GbPoint3D(0.0,0.004,0.001+deltaXfine*8.0),new GbPoint3D(0.064,0.004,0.001+deltaXfine*8.0)));
+ LineTimeSeriesCoProcessor lineTs2(grid, tpcSch,pathname+"/TimeSeries/line2.csv",line2, 3,comm);
+ if (myid==0) lineTs2.writeLine(pathname+"/geo/line2");
+
+ GbLine3DPtr line3(new GbLine3D(new GbPoint3D(0.03,0.0,0.00078),new GbPoint3D(0.03,0.008,0.00078)));
+ LineTimeSeriesCoProcessor lineTs3(grid, tpcSch,pathname+"/TimeSeries/line3.csv",line3, 3,comm);
+ if (myid==0) lineTs3.writeLine(pathname+"/geo/line3");
+
+ GbLine3DPtr line4(new GbLine3D(new GbPoint3D(0.03,0.0,0.001+deltaXfine*8.0),new GbPoint3D(0.03,0.008,0.001+deltaXfine*8.0)));
+ LineTimeSeriesCoProcessor lineTs4(grid, tpcSch,pathname+"/TimeSeries/line4.csv",line4, 3,comm);
+ if (myid==0) lineTs4.writeLine(pathname+"/geo/line4");
+
+ GbLine3DPtr line5(new GbLine3D(new GbPoint3D(0.0,0.004,0.002),new GbPoint3D(0.064,0.004,0.002)));
+ LineTimeSeriesCoProcessor lineTs5(grid, tpcSch,pathname+"/TimeSeries/line5.csv",line5,2,comm);
+ if (myid==0) lineTs5.writeLine(pathname+"/geo/line5");
+
+ GbLine3DPtr line6(new GbLine3D(new GbPoint3D(0.0,0.004,0.0035),new GbPoint3D(0.064,0.004,0.0035)));
+ LineTimeSeriesCoProcessor lineTs6(grid, tpcSch,pathname+"/TimeSeries/line6.csv",line6,1,comm);
+ if (myid==0) lineTs6.writeLine(pathname+"/geo/line6");
+
+ GbLine3DPtr line7(new GbLine3D(new GbPoint3D(0.0,0.004,0.009),new GbPoint3D(0.064,0.004,0.009)));
+ LineTimeSeriesCoProcessor lineTs7(grid, tpcSch,pathname+"/TimeSeries/line7.csv",line7,0,comm);
+ if (myid==0) lineTs7.writeLine(pathname+"/geo/line7");
+
+ GbLine3DPtr line8(new GbLine3D(new GbPoint3D(0.0,0.004,0.015),new GbPoint3D(0.064,0.004,0.015)));
+ LineTimeSeriesCoProcessor lineTs8(grid, tpcSch,pathname+"/TimeSeries/line8.csv",line8,1,comm);
+ if (myid==0) lineTs8.writeLine(pathname+"/geo/line8");
+
+ GbLine3DPtr line9(new GbLine3D(new GbPoint3D(0.0,0.004,0.017),new GbPoint3D(0.064,0.004,0.017)));
+ LineTimeSeriesCoProcessor lineTs9(grid, tpcSch,pathname+"/TimeSeries/line9.csv",line9,2,comm);
+ if (myid==0) lineTs9.writeLine(pathname+"/geo/line9");
+
+ GbLine3DPtr line10(new GbLine3D(new GbPoint3D(0.0,0.004,0.018-deltaXfine*14.0),new GbPoint3D(0.064,0.004,0.018-deltaXfine*14.0)));
+ LineTimeSeriesCoProcessor lineTs10(grid, tpcSch,pathname+"/TimeSeries/line10.csv",line10,3,comm);
+ if (myid==0) lineTs10.writeLine(pathname+"/geo/line10");
+
+ 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, stepSch));
+ if (averaging)
+ {
+ calculation->setTimeAveragedValuesCoProcessor(tav);
+ }
+ 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/CMake/cmake_config_files/BOMBADIL.config.cmake b/source/CMake/cmake_config_files/BOMBADIL.config.cmake
index dabf1253e8b67732010af9bf4818101d5bbf7470..c410eae506d52fb1291967096811dd1a399d10ae 100644
--- a/source/CMake/cmake_config_files/BOMBADIL.config.cmake
+++ b/source/CMake/cmake_config_files/BOMBADIL.config.cmake
@@ -26,7 +26,7 @@ SET(BOOST_LIBRARYDIR ${BOOST_ROOT}"/stageMSVC64/lib")
#################################################################################
# VTK
#################################################################################
-#set(VTK_DIR "c:/Tools/VTK/build/VTK-6.1.0")
+set(VTK_DIR "c:/Tools/VTK/build/VTK-6.1.0")
#################################################################################
# ZOLTAN
#################################################################################
diff --git a/source/ThirdParty/Library/numerics/geometry3d/GbVoxelMatrix3D.cpp b/source/ThirdParty/Library/numerics/geometry3d/GbVoxelMatrix3D.cpp
index 87271118d828ddb95e33cbc7f315ebc603e93445..c7823d0daff5633f8775a799bf38baa83689ae15 100644
--- a/source/ThirdParty/Library/numerics/geometry3d/GbVoxelMatrix3D.cpp
+++ b/source/ThirdParty/Library/numerics/geometry3d/GbVoxelMatrix3D.cpp
@@ -138,7 +138,7 @@ void GbVoxelMatrix3D::setClosedVoidSpaceToSolid()
x1NbrTemp.clear();
x2NbrTemp.clear();
x3NbrTemp.clear();
- size = x1Nbr.size();
+ size = (int)x1Nbr.size();
}
UBLOG(logINFO, "setClosedVoidSpaceToSolid:end");
voxelMatrix = voxelMatrixTemp;
diff --git a/source/ThirdParty/Library/numerics/geometry3d/GbVoxelMatrix3D.h b/source/ThirdParty/Library/numerics/geometry3d/GbVoxelMatrix3D.h
index 2c26ca05a2a6b008fe7391e8b27aed68803ee93b..7eaace198edf1547cee5e6333b6e0c91bd7cc94d 100644
--- a/source/ThirdParty/Library/numerics/geometry3d/GbVoxelMatrix3D.h
+++ b/source/ThirdParty/Library/numerics/geometry3d/GbVoxelMatrix3D.h
@@ -236,7 +236,7 @@ void GbVoxelMatrix3D::readMatrixFromRawFile(std::string filename, GbVoxelMatrix3
//UBLOG(logINFO,"number of nodes = "<<nodesX1*nodesX2*nodesX3*sizeof(T)<<" file size = "<<(long)length);
//if( (nodesX1*nodesX2*nodesX3)*sizeof(float) != (long)length )
- unsigned long long nofn = nodesX1*nodesX2*nodesX3*sizeof(T);
+ unsigned long long nofn = (unsigned long long)nodesX1*(unsigned long long)nodesX2*(unsigned long long)nodesX3*(unsigned long long)sizeof(T);
if (nofn != (unsigned long long)length)
{
throw UbException(UB_EXARGS,"number of nodes("+UbSystem::toString(nofn)+") doesn't match file size("+UbSystem::toString((long)length)+")");
diff --git a/source/VirtualFluids.h b/source/VirtualFluids.h
index 7549b401dc6fff84b2e79c41606481df04da9474..6fb6880b1bb68fb6327df99bdaaba5ea5d7a36f0 100644
--- a/source/VirtualFluids.h
+++ b/source/VirtualFluids.h
@@ -124,6 +124,7 @@
#include <Data/EsoTwist3D.h>
#include <Data/EsoTwistD3Q27System.h>
#include <Data/EsoTwistD3Q27SparseData.h>
+#include <Data/VoidData3D.h>
#include <Grid/Block3D.h>
//#include <Grid/BoostSerializationClassExportHelper.h>
#include <Grid/CalculationManager.h>
@@ -173,6 +174,7 @@
#include <LBM/LBMKernelETD3Q27CCLB.h>
#include <LBM/CompressibleCumulantLBMKernel.h>
#include <LBM/InitDensityLBMKernel.h>
+#include <LBM/VoidLBMKernel.h>
#include <LBM/LBMSystem.h>
#include <LBM/LBMSystems.h>
#include <LBM/LBMUnitConverter.h>
diff --git a/source/VirtualFluidsCore/BoundaryCondition/D3Q27BoundaryCondition.h b/source/VirtualFluidsCore/BoundaryCondition/D3Q27BoundaryCondition.h
index 99d8e63c533a3fa1b57ff4a76c8bc9b72ba13b5f..219d1d835c601bdf4a115a338bc46d8823d5a479 100644
--- a/source/VirtualFluidsCore/BoundaryCondition/D3Q27BoundaryCondition.h
+++ b/source/VirtualFluidsCore/BoundaryCondition/D3Q27BoundaryCondition.h
@@ -125,7 +125,7 @@ public:
{
switch(direction)
{
- case D3Q27System::E : return (float)( UbMath::c4o9*(+bcVelocityX1) ); //(2/cs^2)(=6)*rho_0(=1 bei imkompr)*wi*u*ei mit cs=1/sqrt(3)
+ case D3Q27System::E : return (float)( UbMath::c4o9*(+bcVelocityX1) ); //(2/cs^2)(=6)*rho_0(=1 bei inkompr)*wi*u*ei mit cs=1/sqrt(3)
case D3Q27System::W : return (float)( UbMath::c4o9*(-bcVelocityX1) ); //z.B. aus paper manfred MRT LB models in three dimensions (2002)
case D3Q27System::N : return (float)( UbMath::c4o9*(+bcVelocityX2) );
case D3Q27System::S : return (float)( UbMath::c4o9*(-bcVelocityX2) );
diff --git a/source/VirtualFluidsCore/BoundaryCondition/D3Q27ETBCProcessor.cpp b/source/VirtualFluidsCore/BoundaryCondition/D3Q27ETBCProcessor.cpp
index 1626c5b907207823d31f1cbcd60e6c4fced33911..4b33a7114c153bc21ac762f28dd0b08a81b24292 100644
--- a/source/VirtualFluidsCore/BoundaryCondition/D3Q27ETBCProcessor.cpp
+++ b/source/VirtualFluidsCore/BoundaryCondition/D3Q27ETBCProcessor.cpp
@@ -85,3 +85,8 @@ void D3Q27ETBCProcessor::applyPostCollisionBC()
bc->apply();
}
}
+//////////////////////////////////////////////////////////////////////////
+//void D3Q27ETBCProcessor::resizeBcArray(int nx1, int nx2, int nx3)
+//{
+// bcArray.resize( nx1, nx2, nx3, BCArray3D<D3Q27BoundaryCondition>::FLUID);
+//}
diff --git a/source/VirtualFluidsCore/BoundaryCondition/D3Q27ETBCProcessor.h b/source/VirtualFluidsCore/BoundaryCondition/D3Q27ETBCProcessor.h
index f1f459be7f20b186dd8afed4e5fed1f1ab8febb9..927637d34775e0358562d168f8b017d765d273e8 100644
--- a/source/VirtualFluidsCore/BoundaryCondition/D3Q27ETBCProcessor.h
+++ b/source/VirtualFluidsCore/BoundaryCondition/D3Q27ETBCProcessor.h
@@ -29,6 +29,7 @@ public:
BoundaryConditionPtr getBC(BoundaryCondition::Type type);
void applyPreCollisionBC();
void applyPostCollisionBC();
+ //void resizeBcArray(int nx1, int nx2, int nx3);
//void init();
protected:
std::vector<BoundaryConditionPtr> preBC;
diff --git a/source/VirtualFluidsCore/CoProcessors/LineTimeSeriesCoProcessor.cpp b/source/VirtualFluidsCore/CoProcessors/LineTimeSeriesCoProcessor.cpp
index 53a05c5e094ef338c021113321ccad80cb677e64..c91876a37df54b8533d38ee6b28603cee80e6366 100644
--- a/source/VirtualFluidsCore/CoProcessors/LineTimeSeriesCoProcessor.cpp
+++ b/source/VirtualFluidsCore/CoProcessors/LineTimeSeriesCoProcessor.cpp
@@ -27,12 +27,12 @@ LineTimeSeriesCoProcessor::LineTimeSeriesCoProcessor(Grid3DPtr grid, UbScheduler
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;
+ int x1min = (int)(line->getX1Minimum()-orgX1)/dx;
+ int x1max = (int)(line->getX1Maximum()-orgX1)/dx;
+ int x2min = (int)(line->getX2Minimum()-orgX2)/dx;
+ int x2max = (int)(line->getX2Maximum()-orgX2)/dx;
+ int x3min = (int)(line->getX3Minimum()-orgX3)/dx;
+ int x3max = (int)(line->getX3Maximum()-orgX3)/dx;
UbTupleInt3 blockNx = grid->getBlockNX();
@@ -77,12 +77,12 @@ 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>(nodes[0]) = (float)line->getX1Minimum();
+ val<2>(nodes[0]) = (float)line->getX2Minimum();
+ val<3>(nodes[0]) = (float)line->getX3Minimum();
+ val<1>(nodes[1]) = (float)line->getX1Maximum();
+ val<2>(nodes[1]) = (float)line->getX2Maximum();
+ val<3>(nodes[1]) = (float)line->getX3Maximum();
val<1>(lines[0]) = 0;
val<1>(lines[0]) = 1;
WbWriterVtkXmlASCII *writer = WbWriterVtkXmlASCII::getInstance();
diff --git a/source/VirtualFluidsCore/CoProcessors/RestartCoProcessor.cpp b/source/VirtualFluidsCore/CoProcessors/RestartCoProcessor.cpp
index 5165f778ff4fa65e6df99e912e8cfaf5df2d92b8..c2866958cece78596c1965e6f4f2e78ad5adc242 100644
--- a/source/VirtualFluidsCore/CoProcessors/RestartCoProcessor.cpp
+++ b/source/VirtualFluidsCore/CoProcessors/RestartCoProcessor.cpp
@@ -8,9 +8,12 @@
#include <basics/utilities/UbFileOutputASCII.h>
#include <basics/utilities/UbFileInputASCII.h>
+#include "MetisPartitioningGridVisitor.h"
+
#include "BoostSerializationClassExportHelper.h"
#include <MemoryUtil.h>
+#include <vector>
RestartCoProcessor::RestartCoProcessor(Grid3DPtr& grid, UbSchedulerPtr s, CommunicatorPtr comm, const std::string& path, ArchiveType type) :
CoProcessor(grid, s),
@@ -89,11 +92,11 @@ void RestartCoProcessor::doCheckPoint(int step)
if (archiveType == TXT)
{
- saveTxtArchive(filename + ".txt");
+ saveTxtArchive(filename + ".txt", grid);
}
else if(archiveType == BINARY)
{
- saveBinArchive(filename + ".bin");
+ saveBinArchive(filename + ".bin", grid);
}
comm->barrier();
@@ -167,7 +170,7 @@ void RestartCoProcessor::addBlockVisitor( Block3DVisitorPtr v )
blockVisitors.push_back(v);
}
//////////////////////////////////////////////////////////////////////////
-void RestartCoProcessor::saveTxtArchive(std::string filename)
+void RestartCoProcessor::saveTxtArchive( std::string filename, Grid3DPtr grid )
{
std::ofstream file(filename.c_str());
if(!file)
@@ -214,7 +217,7 @@ void RestartCoProcessor::loadTxtArchive( std::string filename )
file.close();
}
//////////////////////////////////////////////////////////////////////////
-void RestartCoProcessor::saveBinArchive( std::string filename )
+void RestartCoProcessor::saveBinArchive( std::string filename, Grid3DPtr grid )
{
//impotent for binary archive add std::ios::binary
std::ofstream file(filename.c_str(), std::ios::binary);
@@ -285,6 +288,74 @@ void RestartCoProcessor::checkMetafile()
}
file.close();
}
+//////////////////////////////////////////////////////////////////////////
+void RestartCoProcessor::writeDistributedGrid(Grid3DPtr sgrid, int numberOfProcesses)
+{
+ using namespace std;
+
+ Grid3DVisitorPtr metisVisitor(new MetisPartitioningGridVisitor(comm, MetisPartitioningGridVisitor::LevelBased, D3Q27System::BSW, MetisPartitioner::RECURSIVE));
+ boost::dynamic_pointer_cast<MetisPartitioningGridVisitor>(metisVisitor)->setNumberOfProcesses(numberOfProcesses);
+ sgrid->accept(metisVisitor);
+
+ int minInitLevel = sgrid->getCoarsestInitializedLevel();
+ int maxInitLevel = sgrid->getFinestInitializedLevel();
+
+ for (int i = 0; i<numberOfProcesses; i++)
+ {
+ UBLOG(logINFO, "Create dump file for rank " << i <<" - start");
+ Grid3DPtr newGrid(new Grid3D());
+ newGrid->setRank(i);
+ newGrid->setDeltaX(sgrid->getDeltaX(0));
+ newGrid->setNX1(sgrid->getNX1());
+ newGrid->setNX2(sgrid->getNX2());
+ newGrid->setNX3(sgrid->getNX3());
+ newGrid->setCoordinateTransformator(sgrid->getCoordinateTransformator());
+ UbTupleInt3 blockNX = sgrid->getBlockNX();
+ newGrid->setBlockNX(val<1>(blockNX), val<2>(blockNX), val<3>(blockNX));
+ Grid3D::Interactor3DSet interactors = sgrid->getInteractors();
+ for(int inter=0; inter < interactors.size(); inter++)
+ {
+ newGrid->addInteractor(interactors[inter]);
+ }
+
+ for (int level = minInitLevel; level<=maxInitLevel; level++)
+ {
+ vector<Block3DPtr> blockVector;
+ grid->getBlocks(level, blockVector);
+ BOOST_FOREACH(Block3DPtr block, blockVector)
+ {
+ if (block)
+ {
+ if (block->getRank() == i)
+ {
+ newGrid->addBlock(block);
+ }
+ else
+ {
+ Block3DPtr newBlock(new Block3D(block->getX1(), block->getX2(), block->getX3(), block->getLevel()));
+ newBlock->setRank(block->getRank());
+ newGrid->addBlock(newBlock);
+ }
+ }
+ }
+ }
+
+ std::string filename = path+"/checkpoint"+UbSystem::toString(1)+"/checkpoint"+UbSystem::toString(i)+"_"+UbSystem::toString(1);
+
+ if (archiveType==TXT)
+ {
+ saveTxtArchive(filename+".txt", newGrid);
+ }
+ else if (archiveType==BINARY)
+ {
+ saveBinArchive(filename+".bin", newGrid);
+ }
+
+ UBLOG(logINFO, "Create dump file for rank " << i <<" - end");
+ }
+ writeMetafile(1);
+}
+
diff --git a/source/VirtualFluidsCore/CoProcessors/RestartCoProcessor.h b/source/VirtualFluidsCore/CoProcessors/RestartCoProcessor.h
index 541166af86b9aa3972f705525cd9cf3c971c21c2..d8c2f5a4c447ed294af81fd7f035ad5b56576245 100644
--- a/source/VirtualFluidsCore/CoProcessors/RestartCoProcessor.h
+++ b/source/VirtualFluidsCore/CoProcessors/RestartCoProcessor.h
@@ -26,12 +26,13 @@ public:
void addBlockVisitor(Block3DVisitorPtr v);
void doCheckPoint(int step);
Grid3DPtr restart();
+ void writeDistributedGrid(Grid3DPtr grid, int numberOfProcesses);
protected:
void acceptGridVisitors();
void acceptBlockVisitors();
- void saveTxtArchive(std::string filename);
+ void saveTxtArchive(std::string filename, Grid3DPtr grid);
void loadTxtArchive(std::string filename);
- void saveBinArchive(std::string filename);
+ void saveBinArchive(std::string filename, Grid3DPtr grid);
void loadBinArchive(std::string filename);
void writeMetafile(int step);
int readMetafile();
diff --git a/source/VirtualFluidsCore/CoProcessors/TimeAveragedValuesCoProcessor.cpp b/source/VirtualFluidsCore/CoProcessors/TimeAveragedValuesCoProcessor.cpp
index 157fd55f0bd38c1b2b6ef1197d59148b0784eb32..107e2330f16050eecf84045a0788224b6ec58ed8 100644
--- a/source/VirtualFluidsCore/CoProcessors/TimeAveragedValuesCoProcessor.cpp
+++ b/source/VirtualFluidsCore/CoProcessors/TimeAveragedValuesCoProcessor.cpp
@@ -430,16 +430,16 @@ void TimeAveragedValuesCoProcessor::calculateAverageValues(double timeSteps)
if ((options&Triplecorrelations) == Triplecorrelations)
{
//triple-correlations
- (*at)(Vxxx, ix1, ix2, ix3) = (*at)(Vxxx, ix1, ix2, ix3) / timeSteps - 3 * uxx*ux + 2 * ux*ux*ux;
- (*at)(Vxxy, ix1, ix2, ix3) = (*at)(Vxxy, ix1, ix2, ix3) / timeSteps - 2 * uxy*ux - uxx*uy + 2 * ux*ux*uy;
- (*at)(Vxxz, ix1, ix2, ix3) = (*at)(Vxxz, ix1, ix2, ix3) / timeSteps - 2 * uxz*ux - uxx*uz + 2 * ux*ux*uz;
- (*at)(Vyyy, ix1, ix2, ix3) = (*at)(Vyyy, ix1, ix2, ix3) / timeSteps - 3 * uyy*uy + 2 * uy*uy*uy;
- (*at)(Vyyx, ix1, ix2, ix3) = (*at)(Vyyx, ix1, ix2, ix3) / timeSteps - 2 * uxy*uy - uyy*ux + 2 * uy*uy*ux;
- (*at)(Vyyz, ix1, ix2, ix3) = (*at)(Vyyz, ix1, ix2, ix3) / timeSteps - 2 * uyz*uy - uyy*uz + 2 * uy*uy*uz;
- (*at)(Vzzz, ix1, ix2, ix3) = (*at)(Vzzz, ix1, ix2, ix3) / timeSteps - 3 * uzz*uz + 2 * uz*uz*uz;
- (*at)(Vzzx, ix1, ix2, ix3) = (*at)(Vzzx, ix1, ix2, ix3) / timeSteps - 2 * uxz*uz - uzz*ux + 2 * uz*uz*ux;
- (*at)(Vzzy, ix1, ix2, ix3) = (*at)(Vzzy, ix1, ix2, ix3) / timeSteps - 2 * uyz*uz - uzz*uy + 2 * uz*uz*uy;
- (*at)(Vxyz, ix1, ix2, ix3) = (*at)(Vxyz, ix1, ix2, ix3) / timeSteps - uxy*uz - uxz*uy - uyz*ux + 2 * ux*uy*uz;
+ (*at)(Vxxx, ix1, ix2, ix3) = (*at)(Vxxx, ix1, ix2, ix3) / timeSteps - 3.0 * uxx*ux + 2.0 * ux*ux*ux;
+ (*at)(Vxxy, ix1, ix2, ix3) = (*at)(Vxxy, ix1, ix2, ix3) / timeSteps - 2.0 * uxy*ux - uxx*uy + 2.0 * ux*ux*uy;
+ (*at)(Vxxz, ix1, ix2, ix3) = (*at)(Vxxz, ix1, ix2, ix3) / timeSteps - 2.0 * uxz*ux - uxx*uz + 2.0 * ux*ux*uz;
+ (*at)(Vyyy, ix1, ix2, ix3) = (*at)(Vyyy, ix1, ix2, ix3) / timeSteps - 3.0 * uyy*uy + 2.0 * uy*uy*uy;
+ (*at)(Vyyx, ix1, ix2, ix3) = (*at)(Vyyx, ix1, ix2, ix3) / timeSteps - 2.0 * uxy*uy - uyy*ux + 2.0 * uy*uy*ux;
+ (*at)(Vyyz, ix1, ix2, ix3) = (*at)(Vyyz, ix1, ix2, ix3) / timeSteps - 2.0 * uyz*uy - uyy*uz + 2.0 * uy*uy*uz;
+ (*at)(Vzzz, ix1, ix2, ix3) = (*at)(Vzzz, ix1, ix2, ix3) / timeSteps - 3.0 * uzz*uz + 2.0 * uz*uz*uz;
+ (*at)(Vzzx, ix1, ix2, ix3) = (*at)(Vzzx, ix1, ix2, ix3) / timeSteps - 2.0 * uxz*uz - uzz*ux + 2.0 * uz*uz*ux;
+ (*at)(Vzzy, ix1, ix2, ix3) = (*at)(Vzzy, ix1, ix2, ix3) / timeSteps - 2.0 * uyz*uz - uzz*uy + 2.0 * uz*uz*uy;
+ (*at)(Vxyz, ix1, ix2, ix3) = (*at)(Vxyz, ix1, ix2, ix3) / timeSteps - uxy*uz - uxz*uy - uyz*ux + 2.0 * ux*uy*uz;
}
//////////////////////////////////////////////////////////////////////////
}
diff --git a/source/VirtualFluidsCore/Data/VoidData3D.h b/source/VirtualFluidsCore/Data/VoidData3D.h
new file mode 100644
index 0000000000000000000000000000000000000000..482277eb8ce3badc28437a2d989cf81e0c23d9bc
--- /dev/null
+++ b/source/VirtualFluidsCore/Data/VoidData3D.h
@@ -0,0 +1,47 @@
+#ifndef VoidData3D_H
+#define VoidData3D_H
+
+#include "EsoTwist3D.h"
+#include <boost/serialization/serialization.hpp>
+#include <boost/smart_ptr/shared_ptr.hpp>
+
+class VoidData3D;
+typedef boost::shared_ptr<VoidData3D> VoidData3DPtr;
+
+class VoidData3D : public EsoTwist3D
+{
+public:
+ VoidData3D() {};
+ VoidData3D (size_t nx1, size_t nx2, size_t nx3, LBMReal value)
+ {
+ this->NX1 = nx1;
+ this->NX2 = nx2;
+ this->NX3 = nx3;
+ }
+ ~VoidData3D() {};
+ size_t getNX1() const { return NX1;}
+ size_t getNX2() const { return NX2;}
+ size_t getNX3() const { return NX3;}
+ void getDistribution(LBMReal* const f, size_t x1, size_t x2, size_t x3) {}
+ void setDistribution(const LBMReal* const f, size_t x1, size_t x2, size_t x3){}
+ void getDistributionInv(LBMReal* const f, size_t x1, size_t x2, size_t x3) {}
+ void setDistributionInv(const LBMReal* const f, size_t x1, size_t x2, size_t x3) {}
+ void setDistributionForDirection(const LBMReal* const f, size_t x1, size_t x2, size_t x3, unsigned long int direction) {}
+ void setDistributionForDirection(LBMReal f, size_t x1, size_t x2, size_t x3, int direction) {}
+ LBMReal getDistributionInvForDirection(size_t x1, size_t x2, size_t x3, int direction) {return 0.0;}
+ void setDistributionInvForDirection(const LBMReal* const f, size_t x1, size_t x2, size_t x3, unsigned long int direction) {}
+ void setDistributionInvForDirection(LBMReal f, size_t x1, size_t x2, size_t x3, unsigned long int direction) {}
+ LBMReal getDistributionForDirection(size_t x1, size_t x2, size_t x3, int direction) {return 0.0;}
+ void swap() {}
+protected:
+private:
+ size_t NX1, NX2, NX3;
+ friend class boost::serialization::access;
+ template<class Archive>
+ void serialize(Archive & ar, const unsigned int version)
+ {
+ ar & boost::serialization::base_object< EsoTwist3D >(*this);
+ }
+};
+
+#endif
diff --git a/source/VirtualFluidsCore/Grid/BoostSerializationClassExportHelper.h b/source/VirtualFluidsCore/Grid/BoostSerializationClassExportHelper.h
index 97fcd522f5fe7e73b0198dc8e71118948767d0ae..ef0073a3f53e8618726c7e06a12a2494efb65bb0 100644
--- a/source/VirtualFluidsCore/Grid/BoostSerializationClassExportHelper.h
+++ b/source/VirtualFluidsCore/Grid/BoostSerializationClassExportHelper.h
@@ -38,12 +38,15 @@
#include <D3Q27VelocityBCAdapter.h>
#include "ThinWallNoSlipBoundaryCondition.h"
#include "D3Q27ETForThinWallBCProcessor.h"
-
+#include "VoidLBMKernel.h"
+#include "VoidData3D.h"
#include <boost/serialization/export.hpp>
BOOST_CLASS_EXPORT(LBMKernelETD3Q27)
BOOST_CLASS_EXPORT(LBMKernelETD3Q27CCLB)
BOOST_CLASS_EXPORT(LBMKernelETD3Q27CCLBWithSpongeLayer)
+BOOST_CLASS_EXPORT(VoidLBMKernel)
+BOOST_CLASS_EXPORT(VoidData3D)
BOOST_CLASS_EXPORT(D3Q27EsoTwist3DSplittedVector)
BOOST_CLASS_EXPORT(BCProcessor)
BOOST_CLASS_EXPORT(D3Q27ETBCProcessor)
diff --git a/source/VirtualFluidsCore/LBM/D3Q27IntegrateValuesHelper.cpp b/source/VirtualFluidsCore/LBM/D3Q27IntegrateValuesHelper.cpp
index fd1a3ef209a09551c53a32aaaa11d0f93b66adcd..b8c798c34055aacf8913c83a5e72850031a57c46 100644
--- a/source/VirtualFluidsCore/LBM/D3Q27IntegrateValuesHelper.cpp
+++ b/source/VirtualFluidsCore/LBM/D3Q27IntegrateValuesHelper.cpp
@@ -453,46 +453,6 @@ 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()
{
diff --git a/source/VirtualFluidsCore/LBM/D3Q27IntegrateValuesHelper.h b/source/VirtualFluidsCore/LBM/D3Q27IntegrateValuesHelper.h
index f223b07a04558dd5b8b0063d16bcfd7868ed9971..0ba14b1aca727ef5de6a88f73aff43a1db1ecf2f 100644
--- a/source/VirtualFluidsCore/LBM/D3Q27IntegrateValuesHelper.h
+++ b/source/VirtualFluidsCore/LBM/D3Q27IntegrateValuesHelper.h
@@ -49,7 +49,6 @@ public:
void calculateAV();
void calculateAV2();
void prepare2DMatrix(int level);
- void calculateTwoPointCorrelations();
void clearData();
double getRho() {return sRho;}
diff --git a/source/VirtualFluidsCore/LBM/LBMKernel3D.cpp b/source/VirtualFluidsCore/LBM/LBMKernel3D.cpp
index 4a82ef7cd6c4e4bc4bd61586ca6743bef75ef8c1..c503464df8bab9a038ef17e5ed70a1453034bceb 100644
--- a/source/VirtualFluidsCore/LBM/LBMKernel3D.cpp
+++ b/source/VirtualFluidsCore/LBM/LBMKernel3D.cpp
@@ -4,7 +4,8 @@
LBMKernel3D::LBMKernel3D() : ghostLayerWidth(1),
deltaT(1.0),
withForcing(false),
- withSpongeLayer(false)
+ withSpongeLayer(false),
+ compressible(false)
{
this->setForcingX1(0.0);
this->setForcingX2(0.0);
diff --git a/source/VirtualFluidsCore/LBM/VoidLBMKernel.cpp b/source/VirtualFluidsCore/LBM/VoidLBMKernel.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..35b188a7a016b2bed10c3a15135e0d028c889ada
--- /dev/null
+++ b/source/VirtualFluidsCore/LBM/VoidLBMKernel.cpp
@@ -0,0 +1,47 @@
+#include "VoidLBMKernel.h"
+#include "VoidData3D.h"
+
+VoidLBMKernel::VoidLBMKernel()
+{
+
+}
+
+VoidLBMKernel::VoidLBMKernel(int nx1, int nx2, int nx3) : nx1(nx1), nx2(nx2), nx3(nx3)
+{
+ DistributionArray3DPtr d(new VoidData3D(nx1+2, nx2+2, nx3+2, -999.0));
+ dataSet->setFdistributions(d);
+}
+
+VoidLBMKernel::~VoidLBMKernel()
+{
+
+}
+
+LBMKernel3DPtr VoidLBMKernel::clone()
+{
+ LBMKernel3DPtr kernel(new VoidLBMKernel(nx1, nx2, nx3));
+ 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);
+ return kernel;
+}
+
+void VoidLBMKernel::calculate()
+{
+
+}
+
+void VoidLBMKernel::swapDistributions()
+{
+
+}
+
+double VoidLBMKernel::getCallculationTime()
+{
+ return 0.0;
+}
diff --git a/source/VirtualFluidsCore/LBM/VoidLBMKernel.h b/source/VirtualFluidsCore/LBM/VoidLBMKernel.h
new file mode 100644
index 0000000000000000000000000000000000000000..0165e711a8cb216e2969ecc7d5d3e1de79761337
--- /dev/null
+++ b/source/VirtualFluidsCore/LBM/VoidLBMKernel.h
@@ -0,0 +1,28 @@
+#ifndef VoidLBMKernel_h__
+#define VoidLBMKernel_h__
+
+#include "LBMKernel3D.h"
+
+class VoidLBMKernel : public LBMKernel3D
+{
+public:
+ VoidLBMKernel();
+ VoidLBMKernel(int nx1, int nx2, int nx3);
+ ~VoidLBMKernel();
+ LBMKernel3DPtr clone();
+ void calculate();
+ void swapDistributions();
+ double getCallculationTime();
+protected:
+private:
+ int nx1, nx2, nx3;
+ friend class boost::serialization::access;
+ template<class Archive>
+ void serialize(Archive & ar, const unsigned int version)
+ {
+ ar & boost::serialization::base_object<LBMKernel3D>(*this);
+ ar & nx1 & nx2 & nx3;
+ }
+
+};
+#endif // VoidLBMKernel_h__
diff --git a/source/VirtualFluidsCore/Visitors/InitDistributionsFromFileBlockVisitor.cpp b/source/VirtualFluidsCore/Visitors/InitDistributionsFromFileBlockVisitor.cpp
index 98b3dd7d1e8a987be6e4fcc14a019c8d3db5da42..9df3d1d371945fc8fa3b9a7a7fb7a38efc988abb 100644
--- a/source/VirtualFluidsCore/Visitors/InitDistributionsFromFileBlockVisitor.cpp
+++ b/source/VirtualFluidsCore/Visitors/InitDistributionsFromFileBlockVisitor.cpp
@@ -85,9 +85,9 @@ void InitDistributionsFromFileBlockVisitor::visit(const Grid3DPtr grid, Block3DP
int maxX2 = (int)bcArray.getNX2();
int maxX3 = (int)bcArray.getNX3();
- int maxMX1 = matrix.getNX2();
- int maxMX2 = matrix.getNX3();
- int maxMX3 = matrix.getNX4();
+ int maxMX1 = (int)matrix.getNX2();
+ int maxMX2 = (int)matrix.getNX3();
+ int maxMX3 = (int)matrix.getNX4();
int blockix1 = block->getX1();
int blockix2 = block->getX2();
diff --git a/source/VirtualFluidsCore/Visitors/MetisPartitioningGridVisitor.cpp b/source/VirtualFluidsCore/Visitors/MetisPartitioningGridVisitor.cpp
index 2485bf1de29d897dedf772315a4b3aa1441528b6..e5c5aab95a8a32f0b43ea8e337189f378a4a6bb1 100644
--- a/source/VirtualFluidsCore/Visitors/MetisPartitioningGridVisitor.cpp
+++ b/source/VirtualFluidsCore/Visitors/MetisPartitioningGridVisitor.cpp
@@ -15,7 +15,7 @@ MetisPartitioningGridVisitor::MetisPartitioningGridVisitor(CommunicatorPtr comm,
graphType(graphType),
partType(partType)
{
-
+ numberOfProcesses = comm->getNumberOfProcesses();
}
//////////////////////////////////////////////////////////////////////////
MetisPartitioningGridVisitor::~MetisPartitioningGridVisitor()
@@ -44,7 +44,7 @@ void MetisPartitioningGridVisitor::visit(Grid3DPtr grid)
processRoot = comm->getProcessRoot();
processID = comm->getProcessID();
- int numberOfProcesses = comm->getNumberOfProcesses();
+ /*int numberOfProcesses = comm->getNumberOfProcesses();*/
if (numberOfProcesses > 1)
{
int temp = bundleID;
@@ -53,7 +53,7 @@ void MetisPartitioningGridVisitor::visit(Grid3DPtr grid)
if (bundleRoot == bundleID && processRoot == processID)
{
bundleID = i;
- numberOfProcesses = comm->getNumberOfProcessesInBundle(i);
+ //numberOfProcesses = comm->getNumberOfProcessesInBundle(i);
collectData(grid, numberOfProcesses, PROCESS);
bundleID = temp;
}
@@ -294,4 +294,10 @@ void MetisPartitioningGridVisitor::clear()
parts.clear();
}
//////////////////////////////////////////////////////////////////////////
+void MetisPartitioningGridVisitor::setNumberOfProcesses(int np)
+{
+ numberOfProcesses = np;
+}
+
+
#endif //VF_METIS
diff --git a/source/VirtualFluidsCore/Visitors/MetisPartitioningGridVisitor.h b/source/VirtualFluidsCore/Visitors/MetisPartitioningGridVisitor.h
index 34bfbe82f46cde9e663a83488d1f70b3df240a9f..4f00c17caee41c2392f45b8bbab99b61b7ef11e8 100644
--- a/source/VirtualFluidsCore/Visitors/MetisPartitioningGridVisitor.h
+++ b/source/VirtualFluidsCore/Visitors/MetisPartitioningGridVisitor.h
@@ -34,6 +34,7 @@ public:
MetisPartitioningGridVisitor(CommunicatorPtr comm, GraphType graphType, int numOfDirs, MetisPartitioner::PartType partType = MetisPartitioner::KWAY, bool threads = false, int numberOfThreads = 0);
virtual ~MetisPartitioningGridVisitor();
void visit(Grid3DPtr grid);
+ void setNumberOfProcesses(int np);
protected:
enum PartLevel {BUNDLE, PROCESS, THREAD};
void collectData(Grid3DPtr grid, int nofSegments, PartLevel level);
@@ -56,6 +57,7 @@ protected:
bool threads;
GraphType graphType;
MetisPartitioner::PartType partType;
+ int numberOfProcesses;
};
#endif //VF_MPI
diff --git a/source/VirtualFluidsCore/Visitors/SetKernelBlockVisitor.cpp b/source/VirtualFluidsCore/Visitors/SetKernelBlockVisitor.cpp
index 4da48ba57a90c38e34696484b4d849d3b696eb00..608b14e36b2c2f88cc23d8c8c5376a9b42160718 100644
--- a/source/VirtualFluidsCore/Visitors/SetKernelBlockVisitor.cpp
+++ b/source/VirtualFluidsCore/Visitors/SetKernelBlockVisitor.cpp
@@ -18,7 +18,7 @@
//}
//////////////////////////////////////////////////////////////////////////
SetKernelBlockVisitor::SetKernelBlockVisitor(LBMKernel3DPtr kernel, LBMReal nue, double availMem, double needMem, SetKernelBlockVisitor::Action action /*= SetKernelBlockVisitor::New*/) :
- Block3DVisitor(0, Grid3DSystem::MAXLEVEL), kernel(kernel), nue(nue), action(action)
+ Block3DVisitor(0, Grid3DSystem::MAXLEVEL), kernel(kernel), nue(nue), action(action), dataSetFlag(true)
{
if (needMem > availMem)
{
@@ -43,12 +43,15 @@ void SetKernelBlockVisitor::visit(Grid3DPtr grid, Block3DPtr block)
block->setKernel(newKernel);
break;
case SetKernelBlockVisitor::Change:
+ {
DataSet3DPtr dataSet = block->getKernel()->getDataSet();
if (!dataSet)
{
UB_THROW(UbException(UB_EXARGS, "It is not possible to change a DataSet in kernel! Old DataSet is not exist!"));
}
+
newKernel->setDataSet(dataSet);
+
BCProcessorPtr bcProc = block->getKernel()->getBCProcessor();
if (!bcProc)
{
@@ -56,9 +59,27 @@ void SetKernelBlockVisitor::visit(Grid3DPtr grid, Block3DPtr block)
}
newKernel->setBCProcessor(bcProc);
block->setKernel(newKernel);
+ }
+ break;
+
+ case SetKernelBlockVisitor::ChangeBC:
+ {
+ BCProcessorPtr bcProc = block->getKernel()->getBCProcessor();
+ if (!bcProc)
+ {
+ UB_THROW(UbException(UB_EXARGS, "It is not possible to change a BCProcessor in kernel! Old BCProcessor is not exist!"));
+ }
+ newKernel->setBCProcessor(bcProc);
+ block->setKernel(newKernel);
+ }
break;
}
}
}
+void SetKernelBlockVisitor::setNoDataSetFlag(bool flag)
+{
+ dataSetFlag = flag;
+}
+
diff --git a/source/VirtualFluidsCore/Visitors/SetKernelBlockVisitor.h b/source/VirtualFluidsCore/Visitors/SetKernelBlockVisitor.h
index 2b8bcca338a23ba5fa4f44b3c11f8c157bc923dd..d2bdff5fda07d0d85c98154c4a0ef5ff237906e2 100644
--- a/source/VirtualFluidsCore/Visitors/SetKernelBlockVisitor.h
+++ b/source/VirtualFluidsCore/Visitors/SetKernelBlockVisitor.h
@@ -8,7 +8,7 @@
class SetKernelBlockVisitor : public Block3DVisitor
{
public:
- enum Action { New, Change };
+ enum Action { New, Change, ChangeBC};
public:
//SetKernelBlockVisitor(LBMKernel3DPtr kernel, LBMReal nue);
@@ -20,10 +20,13 @@ public:
virtual void visit(Grid3DPtr grid, Block3DPtr block);
+ void setNoDataSetFlag(bool flag);
+
private:
LBMKernel3DPtr kernel;
LBMReal nue;
Action action;
+ bool dataSetFlag;
};
#endif