diff --git a/apps/cpu/Multiphase/Multiphase.cpp b/apps/cpu/Multiphase/Multiphase.cpp
index 9d91cb857b6f0716ef667b44d3742b1e31223bea..32b4295691cfb9355dda300061490d2ee9974063 100644
--- a/apps/cpu/Multiphase/Multiphase.cpp
+++ b/apps/cpu/Multiphase/Multiphase.cpp
@@ -116,12 +116,16 @@ void run(string configname)
// fctF1.SetExpr("vy1*(1-((x1-x0)^2+(x3-z0)^2)/(R^2))");
// fctF1.SetExpr("vy1*(1-(sqrt((x1-x0)^2+(x3-z0)^2)/R))^0.1");
fctF1.SetExpr("vy1");
- fctF1.DefineConst("vy1", -uLB);
+ fctF1.DefineConst("vy1", -uLB*0);
+ //fctF1.DefineConst("vy1", -uLB);
fctF1.DefineConst("R", 8.0);
fctF1.DefineConst("x0", 0.0);
fctF1.DefineConst("z0", 0.0);
- if (newStart) {
+ // mu::Parser fctFStart;
+ // fctFStart.SetExpr("0");
+
+ // if (newStart) {
// bounding box
/*double g_minX1 = 0.0;
@@ -241,8 +245,14 @@ void run(string configname)
// fct.DefineConst("U", uLB);
// BCAdapterPtr velBCAdapter(new VelocityBCAdapter(true, false, false, fct, 0, BCFunction::INFCONST));
+ //Am Anfang keine Geschwindigkeit
+ // double startGeschwindigkeit = 200;
+ // SPtr<BCAdapter> velBCAdapterFStart(
+ // new MultiphaseVelocityBCAdapter(false, true, false, fctFStart, phiH, 0.0, startGeschwindigkeit));
+
+
SPtr<BCAdapter> velBCAdapterF1(
- new MultiphaseVelocityBCAdapter(false, true, false, fctF1, phiH, 0.0, endTime));
+ new MultiphaseVelocityBCAdapter(false, true, false, fctF1, phiH, 0, endTime));
// BCAdapterPtr velBCAdapterF2_1_init(new MultiphaseVelocityBCAdapter(false, false, true, fctF2_1, phiH,
// 0.0, endTime)); BCAdapterPtr velBCAdapterF2_2_init(new MultiphaseVelocityBCAdapter(false, false, true,
@@ -253,6 +263,7 @@ void run(string configname)
// true, fctvel_F2_init, phiL, 0.0, endTime));
velBCAdapterF1->setBcAlgorithm(SPtr<BCAlgorithm>(new MultiphaseVelocityBCAlgorithm()));
+ // velBCAdapterFStart->setBcAlgorithm(SPtr<BCAlgorithm>(new MultiphaseVelocityBCAlgorithm()));
// velBCAdapterF2_1_init->setBcAlgorithm(BCAlgorithmPtr(new MultiphaseVelocityBCAlgorithm()));
// velBCAdapterF2_2_init->setBcAlgorithm(BCAlgorithmPtr(new MultiphaseVelocityBCAlgorithm()));
@@ -267,7 +278,7 @@ void run(string configname)
// BC visitor
MultiphaseBoundaryConditionsBlockVisitor bcVisitor;
bcVisitor.addBC(noSlipBCAdapter);
- bcVisitor.addBC(denBCAdapter);
+ // bcVisitor.addBC(denBCAdapter); //Ohne das BB?
bcVisitor.addBC(velBCAdapterF1);
// bcVisitor.addBC(velBCAdapterF2_1_init);
// bcVisitor.addBC(velBCAdapterF2_2_init);
@@ -283,6 +294,8 @@ void run(string configname)
SPtr<D3Q27Interactor> inflowF1Int(
new D3Q27Interactor(geoInflowF1, grid, velBCAdapterF1, Interactor3D::SOLID));
+ // inflowF1Int->addBCAdapter(velBCAdapterFStart);
+
// D3Q27InteractorPtr inflowF2_1Int_init = D3Q27InteractorPtr(new D3Q27Interactor(geoInflowF2_1, grid,
// velBCAdapterF2_1_init, Interactor3D::SOLID));
@@ -418,8 +431,8 @@ void run(string configname)
//SetConnectorsBlockVisitor setConnsVisitor(comm, true, D3Q27System::ENDDIR, nuLB, iProcessor);
// ConnectorFactoryPtr factory(new Block3DConnectorFactory());
// ConnectorBlockVisitor setConnsVisitor(comm, nuLB, iProcessor, factory);
- TwoDistributionsSetConnectorsBlockVisitor setConnsVisitor(comm);
- grid->accept(setConnsVisitor);
+ // TwoDistributionsSetConnectorsBlockVisitor setConnsVisitor(comm);
+ // grid->accept(setConnsVisitor);
// domain decomposition for threads
// PQueuePartitioningGridVisitor pqPartVisitor(numOfThreads);
@@ -436,37 +449,37 @@ void run(string configname)
if (myid == 0)
UBLOG(logINFO, "Preprocess - end");
- } else {
- if (myid == 0) {
- UBLOG(logINFO, "Parameters:");
- UBLOG(logINFO, "uLb = " << uLB);
- UBLOG(logINFO, "rho = " << rhoLB);
- UBLOG(logINFO, "nuLb = " << nuLB);
- UBLOG(logINFO, "Re = " << Re);
- UBLOG(logINFO, "dx = " << dx);
- UBLOG(logINFO, "number of levels = " << refineLevel + 1);
- UBLOG(logINFO, "numOfThreads = " << numOfThreads);
- UBLOG(logINFO, "path = " << pathname);
- }
-
- rcp.restart((int)restartStep);
- grid->setTimeStep(restartStep);
-
- // BCAdapterPtr velBCAdapter(new VelocityBCAdapter());
- // velBCAdapter->setBcAlgorithm(BCAlgorithmPtr(new VelocityBCAlgorithm()));
- // velBCAdapter->setBcAlgorithm(BCAlgorithmPtr(new VelocityWithDensityBCAlgorithm()));
- // bcVisitor.addBC(velBCAdapter);
- // grid->accept(bcVisitor);
-
- // set connectors
- // InterpolationProcessorPtr iProcessor(new IncompressibleOffsetInterpolationProcessor());
- //InterpolationProcessorPtr iProcessor(new CompressibleOffsetInterpolationProcessor());
- //SetConnectorsBlockVisitor setConnsVisitor(comm, true, D3Q27System::ENDDIR, nuLB, iProcessor);
- //grid->accept(setConnsVisitor);
-
- if (myid == 0)
- UBLOG(logINFO, "Restart - end");
- }
+ //} else {
+ // if (myid == 0) {
+ // UBLOG(logINFO, "Parameters:");
+ // UBLOG(logINFO, "uLb = " << uLB);
+ // UBLOG(logINFO, "rho = " << rhoLB);
+ // UBLOG(logINFO, "nuLb = " << nuLB);
+ // UBLOG(logINFO, "Re = " << Re);
+ // UBLOG(logINFO, "dx = " << dx);
+ // UBLOG(logINFO, "number of levels = " << refineLevel + 1);
+ // UBLOG(logINFO, "numOfThreads = " << numOfThreads);
+ // UBLOG(logINFO, "path = " << pathname);
+ // }
+
+ // rcp.restart((int)restartStep);
+ // grid->setTimeStep(restartStep);
+
+ // // BCAdapterPtr velBCAdapter(new VelocityBCAdapter());
+ // // velBCAdapter->setBcAlgorithm(BCAlgorithmPtr(new VelocityBCAlgorithm()));
+ // // velBCAdapter->setBcAlgorithm(BCAlgorithmPtr(new VelocityWithDensityBCAlgorithm()));
+ // // bcVisitor.addBC(velBCAdapter);
+ // // grid->accept(bcVisitor);
+
+ // // set connectors
+ // // InterpolationProcessorPtr iProcessor(new IncompressibleOffsetInterpolationProcessor());
+ // //InterpolationProcessorPtr iProcessor(new CompressibleOffsetInterpolationProcessor());
+ // //SetConnectorsBlockVisitor setConnsVisitor(comm, true, D3Q27System::ENDDIR, nuLB, iProcessor);
+ // //grid->accept(setConnsVisitor);
+
+ // if (myid == 0)
+ // UBLOG(logINFO, "Restart - end");
+ //}
TwoDistributionsSetConnectorsBlockVisitor setConnsVisitor(comm);
grid->accept(setConnsVisitor);
@@ -478,14 +491,37 @@ void run(string configname)
SPtr<UbScheduler> nupsSch(new UbScheduler(10, 30, 100));
SPtr<NUPSCounterCoProcessor> npr(new NUPSCounterCoProcessor(grid, nupsSch, numOfThreads, comm));
+//////
+ double fluxStart = 500.0;
SPtr<UbScheduler> stepGhostLayer(new UbScheduler(1));
- SPtr<Calculator> calculator(new BasicCalculator(grid, stepGhostLayer, endTime));
+ SPtr<Calculator> calculator(new BasicCalculator(grid, stepGhostLayer, fluxStart));
calculator->addCoProcessor(npr);
calculator->addCoProcessor(pp);
+ calculator->calculate();
+ dynamicPointerCast<MultiphaseVelocityBCAdapter>(velBCAdapterF1)->setNewVelocities(0.0,0.0,endTime, -uLB,0.0,endTime,0.0,0.0,endTime);
+ inflowF1Int->initInteractor();
+
+ // SPtr<D3Q27Interactor> inflowF1Int(
+ // new D3Q27Interactor(geoInflowF1, grid, velBCAdapterF1, Interactor3D::SOLID));
+
+//////
+
+
+
+ // SPtr<UbScheduler> stepGhostLayer(new UbScheduler(1));
+ grid->setTimeStep(fluxStart);
+ SPtr<Calculator> calculator2(new BasicCalculator(grid, stepGhostLayer, endTime));
+ calculator2->addCoProcessor(npr);
+ calculator2->addCoProcessor(pp);
+
+
+
+
+
if (myid == 0)
UBLOG(logINFO, "Simulation-start");
- calculator->calculate();
+ calculator2->calculate();
if (myid == 0)
UBLOG(logINFO, "Simulation-end");
} catch (std::exception &e) {
diff --git a/src/cpu/VirtualFluidsCore/BoundaryConditions/MultiphaseVelocityBCAlgorithm.cpp b/src/cpu/VirtualFluidsCore/BoundaryConditions/MultiphaseVelocityBCAlgorithm.cpp
index 80c1bc6518e5cfe27a4c1af9ca6bd2233d1b89d6..4d8ad0a99309c52a61be4f228f8d7e98e9528647 100644
--- a/src/cpu/VirtualFluidsCore/BoundaryConditions/MultiphaseVelocityBCAlgorithm.cpp
+++ b/src/cpu/VirtualFluidsCore/BoundaryConditions/MultiphaseVelocityBCAlgorithm.cpp
@@ -96,13 +96,15 @@ void MultiphaseVelocityBCAlgorithm::applyBC()
phiBC = bcPtr->getBoundaryPhaseField();
- D3Q27System::calcMultiphaseHeq(htemp, phiBC, vx1, vx2, vx3);
-
+ //D3Q27System::calcMultiphaseHeq(htemp, phiBC, vx1, vx2, vx3);
+ D3Q27System::calcMultiphaseHeq(htemp, phiBC, 0.0, 0.0, 0.0);//16.03.2021 dirty hack!
for (int fdir = D3Q27System::STARTF; fdir<=D3Q27System::ENDF; fdir++)
{
if (bcPtr->hasVelocityBoundaryFlag(fdir))
{
- LBMReal hReturn = htemp[fdir]+h[fdir]-heq[fdir];
+ // LBMReal hReturn = htemp[fdir]+h[fdir]-heq[fdir];
+ //17.03.2021 Let us just set the plain eq
+ LBMReal hReturn = htemp[fdir];
distributionsH->setDistributionForDirection(hReturn, nx1, nx2, nx3, fdir);
}
}
@@ -114,7 +116,9 @@ void MultiphaseVelocityBCAlgorithm::applyBC()
const int invDir = D3Q27System::INVDIR[fdir];
LBMReal q = bcPtr->getQ(invDir);// m+m q=0 stabiler
LBMReal velocity = bcPtr->getBoundaryVelocity(invDir);
- LBMReal fReturn = ((1.0-q)/(1.0+q))*((f[invDir]-feq[invDir])/(1.0-collFactor)+feq[invDir])+((q*(f[invDir]+f[fdir])-velocity)/(1.0+q));
+ //16.03.2021 quick fix for velocity BC
+ LBMReal fReturn = f[invDir] - velocity;
+ //LBMReal fReturn = ((1.0-q)/(1.0+q))*((f[invDir]-feq[invDir])/(1.0-collFactor)+feq[invDir])+((q*(f[invDir]+f[fdir])-velocity)/(1.0+q));
distributions->setDistributionForDirection(fReturn, x1+D3Q27System::DX1[invDir], x2+D3Q27System::DX2[invDir], x3+D3Q27System::DX3[invDir], fdir);
}
}
diff --git a/src/cpu/VirtualFluidsCore/LBM/MultiphaseScratchCumulantLBMKernel.cpp b/src/cpu/VirtualFluidsCore/LBM/MultiphaseScratchCumulantLBMKernel.cpp
index fec58c0a1a7d0fe15647ee2a1d5feececdab0066..a62fef298fc0c18eb053c1138949ab1fb8e58ead 100644
--- a/src/cpu/VirtualFluidsCore/LBM/MultiphaseScratchCumulantLBMKernel.cpp
+++ b/src/cpu/VirtualFluidsCore/LBM/MultiphaseScratchCumulantLBMKernel.cpp
@@ -132,6 +132,8 @@ void MultiphaseScratchCumulantLBMKernel::calculate(int step)
forcingX1 = 0.0;
forcingX2 = 0.0;
forcingX3 = 0.0;
+
+ LBMReal oneOverInterfaceScale = 1.0;// 1.0 / 3.0;
/////////////////////////////////////
localDistributionsF = dynamicPointerCast<D3Q27EsoTwist3DSplittedVector>(dataSet->getFdistributions())->getLocalDistributions();
@@ -341,6 +343,10 @@ void MultiphaseScratchCumulantLBMKernel::calculate(int step)
mfbbb = (*this->zeroDistributionsF)(x1, x2, x3) / rho;
+
+
+
+
LBMReal m0, m1, m2;
LBMReal rhoRef=c1;
@@ -380,6 +386,18 @@ void MultiphaseScratchCumulantLBMKernel::calculate(int step)
vvy += forcingX2 * deltaT * 0.5; // Y
vvz += forcingX3 * deltaT * 0.5; // Z
}
+
+ LBMReal vx2;
+ LBMReal vy2;
+ LBMReal vz2;
+ vx2 = vvx * vvx;
+ vy2 = vvy * vvy;
+ vz2 = vvz * vvz;
+
+ ///////
+
+
+
///////////////////////////////////////////////////////////////////////////////////////////
LBMReal oMdrho;
@@ -410,12 +428,7 @@ void MultiphaseScratchCumulantLBMKernel::calculate(int step)
m0 += mfbbb; //hat gefehlt
oMdrho = (rhoRef - (oMdrho + m0))/rhoRef;// 12.03.21 check derivation!!!!
- LBMReal vx2;
- LBMReal vy2;
- LBMReal vz2;
- vx2 = vvx * vvx;
- vy2 = vvy * vvy;
- vz2 = vvz * vvz;
+
////////////////////////////////////////////////////////////////////////////////////
LBMReal wadjust;
LBMReal qudricLimit = 0.01;
@@ -696,12 +709,15 @@ void MultiphaseScratchCumulantLBMKernel::calculate(int step)
//applying phase field gradients first part:
// mxxPyyPzz += c2o3 * rhoToPhi * (dX1_phi * vvx + dX2_phi * vvy + dX3_phi * vvz);
- mxxPyyPzz += c1o3 * rhoToPhi * (dX1_phi * vvx + dX2_phi * vvy + dX3_phi * vvz); // As in Hesam's code
- mxxMyy += c1o3 * rhoToPhi * (dX1_phi * vvx - dX2_phi * vvy);
- mxxMzz += c1o3 * rhoToPhi * (dX1_phi * vvx - dX3_phi * vvz);
- mfabb += c1o6 * rhoToPhi * (dX2_phi * vvz + dX3_phi * vvy);
- mfbab += c1o6 * rhoToPhi * (dX1_phi * vvz + dX3_phi * vvx);
- mfbba += c1o6 * rhoToPhi * (dX1_phi * vvy + dX2_phi * vvx);
+
+ //17.03.2021 attempt for statililization by assymptotically vanishing bias
+ LBMReal correctionScaling = rhoToPhi /rho;// +0.5;// (vx2 + vy2 + vz2) * 100;// +0.5;//(vx2 + vy2 + vz2)*1000;
+ mxxPyyPzz += (1.0/6.0) * (dX1_phi * vvx + dX2_phi * vvy + dX3_phi * vvz)* correctionScaling; // As in Hesam's code
+ mxxMyy += c1o3 * (dX1_phi * vvx - dX2_phi * vvy)* correctionScaling;
+ mxxMzz += c1o3 * (dX1_phi * vvx - dX3_phi * vvz) * correctionScaling;
+ mfabb += c1o6 * (dX2_phi * vvz + dX3_phi * vvy) * correctionScaling;
+ mfbab += c1o6 * (dX1_phi * vvz + dX3_phi * vvx) * correctionScaling;
+ mfbba += c1o6 * (dX1_phi * vvy + dX2_phi * vvx) * correctionScaling;
LBMReal dxux = 0.0;// -c1o2 * collFactorM * (mxxMyy + mxxMzz) + c1o2 * OxxPyyPzz * (/*mfaaa*/ -mxxPyyPzz);
LBMReal dyuy = 0.0;// dxux + collFactorM * c3o2 * mxxMyy;
@@ -718,18 +734,18 @@ void MultiphaseScratchCumulantLBMKernel::calculate(int step)
//applying phase field gradients second part:
//mxxPyyPzz += c2o3 * rhoToPhi * (dX1_phi * vvx + dX2_phi * vvy + dX3_phi * vvz);
- mxxPyyPzz += c1o3 * rhoToPhi * (dX1_phi * vvx + dX2_phi * vvy + dX3_phi * vvz); // As in Hesam's code
- mxxMyy += c1o3 * rhoToPhi * (dX1_phi * vvx - dX2_phi * vvy);
- mxxMzz += c1o3 * rhoToPhi * (dX1_phi * vvx - dX3_phi * vvz);
- mfabb += c1o6 * rhoToPhi * (dX2_phi * vvz + dX3_phi * vvy);
- mfbab += c1o6 * rhoToPhi * (dX1_phi * vvz + dX3_phi * vvx);
- mfbba += c1o6 * rhoToPhi * (dX1_phi * vvy + dX2_phi * vvx);
+ mxxPyyPzz += (1.0 / 6.0) * (dX1_phi * vvx + dX2_phi * vvy + dX3_phi * vvz) * correctionScaling; // As in Hesam's code
+ mxxMyy += c1o3 * (dX1_phi * vvx - dX2_phi * vvy) * correctionScaling;
+ mxxMzz += c1o3 * (dX1_phi * vvx - dX3_phi * vvz) * correctionScaling;
+ mfabb += c1o6 * (dX2_phi * vvz + dX3_phi * vvy) * correctionScaling;
+ mfbab += c1o6 * (dX1_phi * vvz + dX3_phi * vvx) * correctionScaling;
+ mfbba += c1o6 * (dX1_phi * vvy + dX2_phi * vvx) * correctionScaling;
////updated pressure
- mfaaa += rhoToPhi * (dX1_phi * vvx + dX2_phi * vvy + dX3_phi * vvz);
+ mfaaa += (dX1_phi * vvx + dX2_phi * vvy + dX3_phi * vvz) * correctionScaling;
- // mxxPyyPzz += mfaaa;//12.03.21 shifted by mfaaa
- mxxPyyPzz = mfaaa; //12.03.21 reguarized pressure !?
+ mxxPyyPzz += mfaaa;//12.03.21 shifted by mfaaa
+ // mxxPyyPzz = mfaaa; //12.03.21 reguarized pressure !?
// linear combinations back
mfcaa = c1o3 * (mxxMyy + mxxMzz + mxxPyyPzz);
mfaca = c1o3 * (-2. * mxxMyy + mxxMzz + mxxPyyPzz);
@@ -1047,7 +1063,7 @@ void MultiphaseScratchCumulantLBMKernel::calculate(int step)
+ (mfaab + mfacb + mfcab + mfccb) + (mfaba + mfabc + mfcba + mfcbc) + (mfbaa + mfbac + mfbca + mfbcc)
+ (mfabb + mfcbb) + (mfbab + mfbcb) + (mfbba + mfbbc) + mfbbb;
//LBMReal dif = fabs(drho - rho_post);
- LBMReal dif = drho+ rhoToPhi * (dX1_phi * vvx + dX2_phi * vvy + dX3_phi * vvz) - rho_post;
+ LBMReal dif = drho+ (dX1_phi * vvx + dX2_phi * vvy + dX3_phi * vvz) *correctionScaling - rho_post;
#ifdef SINGLEPRECISION
if (dif > 10.0E-7 || dif < -10.0E-7)
#else
@@ -2083,9 +2099,9 @@ void MultiphaseScratchCumulantLBMKernel::calculate(int step)
LBMReal Mccb = mfccb - mfaab * c1o9;
// collision of 1st order moments
- cx = cx * (c1 - omegaD) + omegaD * vvx * concentration + normX1 * (c1 - 0.5*omegaD) * (1.0 - phi[REST]) * (phi[REST])*c1o3;
- cy = cy * (c1 - omegaD) + omegaD * vvy * concentration + normX2 * (c1 - 0.5 * omegaD) * (1.0 - phi[REST]) * (phi[REST]) * c1o3;
- cz = cz * (c1 - omegaD) + omegaD * vvz * concentration + normX3 * (c1 - 0.5 * omegaD) * (1.0 - phi[REST]) * (phi[REST]) * c1o3;
+ cx = cx * (c1 - omegaD) + omegaD * vvx * concentration + normX1 * (c1 - 0.5 * omegaD) * (1.0 - phi[REST]) * (phi[REST]) * c1o3 * oneOverInterfaceScale;
+ cy = cy * (c1 - omegaD) + omegaD * vvy * concentration + normX2 * (c1 - 0.5 * omegaD) * (1.0 - phi[REST]) * (phi[REST]) * c1o3 * oneOverInterfaceScale;
+ cz = cz * (c1 - omegaD) + omegaD * vvz * concentration + normX3 * (c1 - 0.5 * omegaD) * (1.0 - phi[REST]) * (phi[REST]) * c1o3 * oneOverInterfaceScale;
//mhx = (ux * phi[REST] + normX1 * (tauH - 0.5) * (1.0 - phi[REST]) * (phi[REST])) / tauH + (1.0 - 1.0 / tauH) * mhx;
//mhy = (uy * phi[REST] + normX2 * (tauH - 0.5) * (1.0 - phi[REST]) * (phi[REST])) / tauH + (1.0 - 1.0 / tauH) * mhy;