From 437cf3df178ef900ca1d1541560ae0b5870c9237 Mon Sep 17 00:00:00 2001
From: kutscher <kutscher@irmb.tu-bs.de>
Date: Mon, 4 Jun 2018 13:22:47 +0200
Subject: [PATCH] add ForceCalculator, shift setTimeStep to CoProcessor
---
.../pChannel/configBombadilpChannel.cfg | 13 +-
source/Applications/pChannel/pChannel.cpp | 2 +-
.../CoProcessors/ForceCalculator.cpp | 142 ++++++++++++++++++
.../CoProcessors/ForceCalculator.h | 43 ++++++
.../MPIIOMigrationCoProcessor.cpp | 2 +
.../CoProcessors/MPIIORestartCoProcessor.cpp | 2 +
6 files changed, 197 insertions(+), 7 deletions(-)
create mode 100644 source/VirtualFluidsCore/CoProcessors/ForceCalculator.cpp
create mode 100644 source/VirtualFluidsCore/CoProcessors/ForceCalculator.h
diff --git a/source/Applications/pChannel/configBombadilpChannel.cfg b/source/Applications/pChannel/configBombadilpChannel.cfg
index 69f20a844..96f3f4ef9 100644
--- a/source/Applications/pChannel/configBombadilpChannel.cfg
+++ b/source/Applications/pChannel/configBombadilpChannel.cfg
@@ -11,6 +11,7 @@ logToFile = false
#porous media
#rawFile = false
sampleFilename = PA80-110_275x267x254_1mm.vti
+#sampleFilename = PA80-110_1096x1062x254_4x4x1mm.vti
#sampleFilename = PA80-110_1096x1327x303.vti
#writeSample = true
@@ -52,12 +53,12 @@ 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 8e-3 1e-3
+#pmL = 4e-3 4e-3 1e-3
#grid
refineLevel = 0
#deltaXfine = 10e-6
-deltaXfine = 128e-6
+deltaXfine = 20e-6
#deltaXfine = 32e-6 #level 2
blocknx = 10 10 10
@@ -69,21 +70,21 @@ changeQs = false
#DLR-F15
#channelHigh = 0.017
-channelHigh = 0.015
+channelHigh = 0.002
#NACA 0012
#channelHigh = 0.008
#channelHigh = 0.005
-boundingBox = 0 0 0 0.024 0.016 0.016
+boundingBox = 0 0 0 0.006 0.003 0.003
#physic
# for DLRF-15 Re = 102000/2
-Re = 15000
+Re = 25000
#real velocity is 54.95 m/s
u_LB = 0.1
-newStart = false
+newStart = true
restartStep = 230000
cpStep = 100
diff --git a/source/Applications/pChannel/pChannel.cpp b/source/Applications/pChannel/pChannel.cpp
index 2dad90343..56a519953 100644
--- a/source/Applications/pChannel/pChannel.cpp
+++ b/source/Applications/pChannel/pChannel.cpp
@@ -428,7 +428,7 @@ void run(string configname)
SPtr<UbScheduler> AdjForcSch(new UbScheduler());
AdjForcSch->addSchedule(10, 0, 10000000);
- SPtr<IntegrateValuesHelper> intValHelp(new IntegrateValuesHelper(grid, comm, g_minX1, g_minX2, g_minX3, g_maxX1, g_maxX2, g_maxX3));
+ SPtr<IntegrateValuesHelper> intValHelp(new IntegrateValuesHelper(grid, comm, g_minX1, g_minX2, g_minX3+pmL[2], g_maxX1, g_maxX2, g_maxX3));
if (myid == 0) GbSystem3D::writeGeoObject(intValHelp->getBoundingBox().get(), pathOut + "/geo/IntValHelp", WbWriterVtkXmlBinary::getInstance());
double vxTarget=u_LB;
diff --git a/source/VirtualFluidsCore/CoProcessors/ForceCalculator.cpp b/source/VirtualFluidsCore/CoProcessors/ForceCalculator.cpp
new file mode 100644
index 000000000..a96372748
--- /dev/null
+++ b/source/VirtualFluidsCore/CoProcessors/ForceCalculator.cpp
@@ -0,0 +1,142 @@
+#include "ForceCalculator.h"
+#include "BCProcessor.h"
+
+#include "Communicator.h"
+#include "D3Q27Interactor.h"
+#include "DataSet3D.h"
+#include "LBMKernel.h"
+#include "Block3D.h"
+#include "BoundaryConditions.h"
+#include "BCArray3D.h"
+#include "Communicator.h"
+#include "DistributionArray3D.h"
+#include "BoundaryConditions.h"
+#include "D3Q27Interactor.h"
+
+ForceCalculator::ForceCalculator(SPtr<Communicator> comm) : comm(comm), forceX1global(0), forceX2global(0), forceX3global(0)
+{
+
+}
+
+ForceCalculator::~ForceCalculator()
+{
+
+}
+
+
+
+
+Vector3D ForceCalculator::getForces(int x1, int x2, int x3, SPtr<DistributionArray3D> distributions, SPtr<BoundaryConditions> bc, const Vector3D& boundaryVelocity) const
+{
+ double forceX1 = 0;
+ double forceX2 = 0;
+ double forceX3 = 0;
+ if (bc)
+ {
+ for (int fdir = D3Q27System::FSTARTDIR; fdir <= D3Q27System::FENDDIR; fdir++)
+ {
+ if (bc->hasNoSlipBoundaryFlag(fdir) || bc->hasVelocityBoundaryFlag(fdir))
+ {
+ const int invDir = D3Q27System::INVDIR[fdir];
+ const double f = distributions->getDistributionInvForDirection(x1, x2, x3, invDir);
+ const double fnbr = distributions->getDistributionInvForDirection(x1 + D3Q27System::DX1[invDir], x2 + D3Q27System::DX2[invDir], x3 + D3Q27System::DX3[invDir], fdir);
+
+ double correction[3] = { 0.0, 0.0, 0.0 };
+ if(bc->hasVelocityBoundaryFlag(fdir))
+ {
+ const double forceTerm = f - fnbr;
+ correction[0] = forceTerm * boundaryVelocity[0];
+ correction[1] = forceTerm * boundaryVelocity[1];
+ correction[2] = forceTerm * boundaryVelocity[2];
+ }
+
+ //UBLOG(logINFO, "c, c * bv(x,y,z): " << correction << ", " << correction * val<1>(boundaryVelocity) << ", " << correction * val<2>(boundaryVelocity) << ", " << correction * val<3>(boundaryVelocity));
+
+ // force consists of the MEM part and the galilean invariance correction including the boundary velocity
+ forceX1 += (f + fnbr) * D3Q27System::DX1[invDir] - correction[0];
+ forceX2 += (f + fnbr) * D3Q27System::DX2[invDir] - correction[1];
+ forceX3 += (f + fnbr) * D3Q27System::DX3[invDir] - correction[2];
+ }
+ }
+ }
+ return Vector3D(forceX1, forceX2, forceX3);
+}
+
+void ForceCalculator::calculateForces(std::vector<SPtr<D3Q27Interactor> > interactors)
+{
+ forceX1global = 0.0;
+ forceX2global = 0.0;
+ forceX3global = 0.0;
+
+ for (const SPtr<D3Q27Interactor> interactor : interactors)
+ {
+ for (BcNodeIndicesMap::value_type t : interactor->getBcNodeIndicesMap())
+ {
+ double forceX1 = 0.0;
+ double forceX2 = 0.0;
+ double forceX3 = 0.0;
+
+ SPtr<Block3D>block = t.first;
+ SPtr<ILBMKernel> kernel = block->getKernel();
+ SPtr<BCArray3D> bcArray = kernel->getBCProcessor()->getBCArray();
+ SPtr<DistributionArray3D> distributions = kernel->getDataSet()->getFdistributions();
+ distributions->swap();
+
+ std::set< std::vector<int> >& transNodeIndices = t.second;
+ for (std::vector<int> node : transNodeIndices)
+ {
+ int x1 = node[0];
+ int x2 = node[1];
+ int x3 = node[2];
+
+ if (kernel->isInsideOfDomain(x1, x2, x3) && bcArray->isFluid(x1, x2, x3))
+ {
+ SPtr<BoundaryConditions> bc = bcArray->getBC(x1, x2, x3);
+ Vector3D forceVec = getForces(x1, x2, x3, distributions, bc);
+ forceX1 += forceVec[0];
+ forceX2 += forceVec[1];
+ forceX3 += forceVec[2];
+ }
+ }
+ //if we have got discretization with more level
+ // deltaX is LBM deltaX and equal LBM deltaT
+ double deltaX = LBMSystem::getDeltaT(block->getLevel()); //grid->getDeltaT(block);
+ double deltaXquadrat = deltaX*deltaX;
+ forceX1 *= deltaXquadrat;
+ forceX2 *= deltaXquadrat;
+ forceX3 *= deltaXquadrat;
+
+ distributions->swap();
+
+ forceX1global += forceX1;
+ forceX2global += forceX2;
+ forceX3global += forceX3;
+ }
+ }
+ gatherGlobalForces();
+}
+
+void ForceCalculator::gatherGlobalForces()
+{
+ std::vector<double> values{ forceX1global , forceX2global, forceX3global };
+ std::vector<double> rvalues = comm->gather(values);
+
+ if (comm->isRoot())
+ {
+ forceX1global = 0.0;
+ forceX2global = 0.0;
+ forceX3global = 0.0;
+
+ for (int i = 0; i < (int)rvalues.size(); i += 3)
+ {
+ forceX1global += rvalues[i];
+ forceX2global += rvalues[i + 1];
+ forceX3global += rvalues[i + 2];
+ }
+ }
+}
+
+Vector3D ForceCalculator::getGlobalForces() const
+{
+ return Vector3D(forceX1global, forceX2global, forceX3global);
+}
diff --git a/source/VirtualFluidsCore/CoProcessors/ForceCalculator.h b/source/VirtualFluidsCore/CoProcessors/ForceCalculator.h
new file mode 100644
index 000000000..adcb66669
--- /dev/null
+++ b/source/VirtualFluidsCore/CoProcessors/ForceCalculator.h
@@ -0,0 +1,43 @@
+/*
+ * ForceCalculator.h
+ *
+ * Created on: 25.10.2017
+ * Author: S. Peters
+ */
+
+#ifndef ForceCalculator_H
+#define ForceCalculator_H
+
+#include <vector>
+#include <memory>
+
+#include "VirtualFluidsBasics/basics/utilities/Vector3D.h"
+
+class D3Q27Interactor;
+class Communicator;
+class DistributionArray3D;
+class BoundaryConditions;
+
+class ForceCalculator
+{
+public:
+ ForceCalculator(std::shared_ptr<Communicator> comm);
+ virtual ~ForceCalculator();
+
+ void calculateForces(std::vector<std::shared_ptr<D3Q27Interactor> > interactors);
+ Vector3D getForces(int x1, int x2, int x3, std::shared_ptr<DistributionArray3D> distributions, std::shared_ptr<BoundaryConditions> bc, const Vector3D& boundaryVelocity = Vector3D(0.0, 0.0, 0.0)) const;
+
+ Vector3D getGlobalForces() const;
+
+private:
+ void gatherGlobalForces();
+
+ std::shared_ptr<Communicator> comm;
+
+ double forceX1global;
+ double forceX2global;
+ double forceX3global;
+};
+
+
+#endif
diff --git a/source/VirtualFluidsCore/CoProcessors/MPIIOMigrationCoProcessor.cpp b/source/VirtualFluidsCore/CoProcessors/MPIIOMigrationCoProcessor.cpp
index cdea328e6..68e1cde14 100644
--- a/source/VirtualFluidsCore/CoProcessors/MPIIOMigrationCoProcessor.cpp
+++ b/source/VirtualFluidsCore/CoProcessors/MPIIOMigrationCoProcessor.cpp
@@ -1475,6 +1475,8 @@ void MPIIOMigrationCoProcessor::restart(int step)
readDataSet(step);
readBoundaryConds(step);
+ grid->setTimeStep(step);
+
if (comm->isRoot()) UBLOG(logINFO, "Load check point - end");
//this->reconnect(grid);
}
diff --git a/source/VirtualFluidsCore/CoProcessors/MPIIORestartCoProcessor.cpp b/source/VirtualFluidsCore/CoProcessors/MPIIORestartCoProcessor.cpp
index 0188db222..f5a50b6a9 100644
--- a/source/VirtualFluidsCore/CoProcessors/MPIIORestartCoProcessor.cpp
+++ b/source/VirtualFluidsCore/CoProcessors/MPIIORestartCoProcessor.cpp
@@ -1651,6 +1651,8 @@ void MPIIORestartCoProcessor::restart(int step)
readDataSet(step);
readBoundaryConds(step);
+ grid->setTimeStep(step);
+
if (comm->isRoot()) UBLOG(logINFO, "Load check point - end");
}
--
GitLab