From 23f7f0e823d2a91fe5a3c889e8112b8bd4441cc4 Mon Sep 17 00:00:00 2001
From: Konstantin Kutscher <kutscher@irmb.tu-bs.de>
Date: Tue, 24 May 2016 16:26:00 +0000
Subject: [PATCH] PID controller implemented in AdjustForcingCoProcessor
---
source/Applications/FNG/fng.cpp | 74 ++++++++--------
.../pChannel/configBombadilpChannel.cfg | 14 ++--
source/Applications/pChannel/pChannel.cpp | 11 ++-
.../CoProcessors/AdjustForcingCoProcessor.cpp | 84 ++++++++++++++-----
.../CoProcessors/AdjustForcingCoProcessor.h | 22 ++++-
.../TimeAveragedValuesCoProcessor.cpp | 2 +-
6 files changed, 134 insertions(+), 73 deletions(-)
diff --git a/source/Applications/FNG/fng.cpp b/source/Applications/FNG/fng.cpp
index 8be1bb26d..620aeb85f 100644
--- a/source/Applications/FNG/fng.cpp
+++ b/source/Applications/FNG/fng.cpp
@@ -8,50 +8,56 @@
using namespace std;
-void setup(const char *cstr1, const char *cstr2)
+void run(string configname)
{
try
{
- //Sleep(30000);
-
- ConfigFileReader cf(cstr1);
- if ( !cf.read() )
- {
- std::string exceptionText = "Unable to read configuration file\n";
- throw exceptionText;
- }
-
- //parameters from config file
- string machine = cf.getValue("machine");
- string pathname = cf.getValue("path");
- string geoFile = cf.getValue("geoFile");
- int numOfThreads = UbSystem::stringTo<int>(cf.getValue("numOfThreads"));
- double availMem = UbSystem::stringTo<double>(cf.getValue("availMem"));
- int refineLevel = UbSystem::stringTo<int>(cf.getValue("refineLevel"));
- int blocknx = UbSystem::stringTo<int>(cf.getValue("blocknx"));
+ ConfigurationFile config;
+ config.load(configname);
+
+ string pathname = config.getString("pathname");
+ string pathGeo = config.getString("pathGeo");
+ string geoFile1 = config.getString("geoFile1");
+ string geoFile2 = config.getString("geoFile2");
+ int numOfThreads = config.getInt("numOfThreads");
+ 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");
+ int refineLevel = config.getInt("refineLevel");
+ bool logToFile = config.getBool("logToFile");
CommunicatorPtr comm = MPICommunicator::getInstance();
int myid = comm->getProcessID();
- if(machine == "Bombadil") int dumy=0;
- else if(machine == "Ludwig" || machine == "HLRN")
+ if (logToFile)
{
- if(myid ==0)
+#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";
+ logFilename << pathname + "/logfile" + UbSystem::toString(UbSystem::getTimeStamp()) + ".txt";
UbLog::output_policy::setStream(logFilename.str());
}
}
- else throw UbException(UB_EXARGS, "unknown machine");
- GbTriFaceMesh3DPtr geo (GbTriFaceMesh3DCreator::getInstance()->readMeshFromSTLFile(geoFile,"geo"));
+ GbTriFaceMesh3DPtr geo (GbTriFaceMesh3DCreator::getInstance()->readMeshFromSTLFile(geoFile1,"geo"));
if(myid == 0) GbSystem3D::writeGeoObject(geo.get(), pathname+"/geo/geo", WbWriterVtkXmlASCII::getInstance());
double dx = (fabs(geo->getX3Maximum()-geo->getX3Minimum())*10e-3)*(double)(1<<refineLevel);
dx /= 4.0;
- double blockLength = blocknx*dx;
+ double blockLength = blocknx[0]*dx;
double offsetX1 = fabs(geo->getX1Maximum()-geo->getX1Minimum());
double h = fabs(geo->getX3Maximum()-geo->getX3Minimum());
@@ -69,18 +75,18 @@ void setup(const char *cstr1, const char *cstr2)
//##########################################################################
//## physical parameters
//##########################################################################
- double Re = 1e6;
+ double Re = 1e6;
double rhoLB = 0.0;
- double rhoReal = 1.0;
- double nueReal = 0.000015;//0.015;
+ double rhoReal = 1.2041; //(kg/m3)
+ double nueReal = 153.5e-7; //m^2/s
- double lReal = 3.0;//<-m ;//Profile laenge in cm(! cm nicht m !)
+ double lReal = 3.0;//m
double uReal = Re*nueReal/lReal;
//##Machzahl:
//#Ma = uReal/csReal
- double Ma = 0.1;//Ma-Real!
+ double Ma = 0.15;//Ma-Real!
double csReal = uReal/Ma;
double hLB = lReal/dx;
@@ -98,7 +104,7 @@ void setup(const char *cstr1, const char *cstr2)
//////////////////////////////////////////////////////////////////////////
Grid3DPtr grid(new Grid3D(comm));
grid->setDeltaX(dx);
- grid->setBlockNX(blocknx, blocknx, blocknx);
+ grid->setBlockNX(blocknx[0], blocknx[1], blocknx[2]);
GbObject3DPtr gridCube(new GbCuboid3D(g_minX1, g_minX2, g_minX3, g_maxX1, g_maxX2, g_maxX3));
//gridCube->setCenterCoordinates(geo->getX1Centroid(), geo->getX2Centroid(), geo->getX3Centroid());
@@ -115,7 +121,7 @@ void setup(const char *cstr1, const char *cstr2)
UbSchedulerPtr rSch(new UbScheduler());
rSch->addSchedule(50,50,50);
- RestartPostprocessorPtr rp(new RestartPostprocessor(grid, rSch, comm, pathname+"/checkpoints", RestartPostprocessor::TXT));
+ RestartCoProcessorPtr rp(new RestartCoProcessor(grid, rSch, comm, pathname+"/checkpoints", RestartCoProcessor::TXT));
std::string opt;
@@ -399,9 +405,9 @@ int main(int argc, char* argv[])
if ( argv != NULL )
{
- if (argc > 1)
+ if (argv[1] != NULL)
{
- setup(argv[1], argv[2]);
+ run(string(argv[1]));
}
else
{
diff --git a/source/Applications/pChannel/configBombadilpChannel.cfg b/source/Applications/pChannel/configBombadilpChannel.cfg
index 548e8ca59..c1ae0bf88 100644
--- a/source/Applications/pChannel/configBombadilpChannel.cfg
+++ b/source/Applications/pChannel/configBombadilpChannel.cfg
@@ -2,7 +2,7 @@
#Simulation parameters for porous channel
#
-pathname = d:/temp/pChannel5
+pathname = d:/temp/pChannel6
pathGeo = d:/Projects/SFB880/GeometrienPoroeseMedien/PA80-110
numOfThreads = 4
availMem = 4e9
@@ -66,7 +66,7 @@ blocknx = 10 10 10
#blocknx = 32 40 20
lengthFactor = 4
thinWall = false
-forcing = 1e-6
+forcing = 7.83841e-09
changeQs = false
channelHigh = 0.002
@@ -82,15 +82,15 @@ channelHigh = 0.002
# for DLRF-15 Re = 102000/2
Re = 51000
#real velocity is 54.95 m/s
-u_LB = 0.1
+u_LB = 0.01
-restartStep = 3
-restartStepStart = 3
+restartStep = 500
+restartStepStart = 500
-timeAvStart = 2
+timeAvStart = 1
timeAvStop = 5
-endTime = 5
+endTime = 1000
outTime = 1000
diff --git a/source/Applications/pChannel/pChannel.cpp b/source/Applications/pChannel/pChannel.cpp
index 03b985410..6552f9866 100644
--- a/source/Applications/pChannel/pChannel.cpp
+++ b/source/Applications/pChannel/pChannel.cpp
@@ -372,11 +372,14 @@ void run(string configname)
grid->accept(bcVisitor);
mu::Parser inflowProfile;
- //inflowProfile.SetExpr("x3 < h ? 0.0 : uLB+1*x1-1*x2");
+ //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]);
+ //inflowProfile.DefineConst("uLB", u_LB);
+ inflowProfile.DefineConst("uLB", 0.0116);
+ //inflowProfile.DefineConst("h", pmL[2]);
+
+ //inflowProfile = Utilities::getDuctParaboloidX(g_maxX2/2.0, g_maxX2, pmL[2]+channelHigh/2.0, channelHigh, 0.116);
D3Q27ETInitDistributionsBlockVisitor initVisitor(nu_LB, rho_LB);
initVisitor.setVx1(inflowProfile);
@@ -490,7 +493,7 @@ void run(string configname)
UbSchedulerPtr AdjForcSch(new UbScheduler());
- AdjForcSch->addSchedule(1, 0, 10000000);
+ AdjForcSch->addSchedule(10, 0, 10000000);
D3Q27IntegrateValuesHelperPtr intValHelp(new D3Q27IntegrateValuesHelper(grid, comm,
coord[0], coord[1], coord[2],
coord[3], coord[4], coord[5]));
diff --git a/source/VirtualFluidsCore/CoProcessors/AdjustForcingCoProcessor.cpp b/source/VirtualFluidsCore/CoProcessors/AdjustForcingCoProcessor.cpp
index a233a680b..f60c7c1ef 100644
--- a/source/VirtualFluidsCore/CoProcessors/AdjustForcingCoProcessor.cpp
+++ b/source/VirtualFluidsCore/CoProcessors/AdjustForcingCoProcessor.cpp
@@ -2,7 +2,7 @@
* D3Q27AdjustForcingCoProcessor.cpp
*
*
-* Author: Sonja Uphoff
+* Author: Konstantin Kutscher
*/
#include "AdjustForcingCoProcessor.h"
@@ -19,18 +19,17 @@ using namespace std;
AdjustForcingCoProcessor::AdjustForcingCoProcessor(Grid3DPtr grid, UbSchedulerPtr s,
const std::string& path,
D3Q27IntegrateValuesHelperPtr integrateValues,
- double vTarged,
- double forcing,
+ double vTarged, double forcing,
CommunicatorPtr comm)
: CoProcessor(grid, s),
path(path),
integrateValues(integrateValues),
comm(comm),
- vTarged(vTarged),
+ vx1Targed(vTarged),
forcing(forcing)
{
- cnodes = integrateValues->getCNodes();
+ //cnodes = integrateValues->getCNodes();
if (comm->getProcessID() == comm->getRoot())
{
std::string fname = path+"/forcing/forcing.csv";
@@ -43,9 +42,37 @@ AdjustForcingCoProcessor::AdjustForcingCoProcessor(Grid3DPtr grid, UbSchedulerPt
if (path.size()>0) { UbSystem::makeDirectory(path); ostr.open(fname.c_str(), std::ios_base::out | std::ios_base::app); }
if (!ostr) throw UbException(UB_EXARGS, "couldn't open file "+fname);
}
- ostr << "step;volume;vx1average;factor;forcing\n";
+ ostr << "step;volume;vx1average;forcing\n";
ostr.close();
}
+
+ Ta = scheduler->getMaxStep();
+
+ Kpcrit = 3.0 / Ta;// 0.3;
+ Tcrit = 3.0 * Ta; // 30.0;
+ Tn = 0.5 * Tcrit;
+ Tv = 0.12 * Tcrit;
+
+ Kp = 0.6 * Kpcrit;
+ Ki = Kp / Tn;
+ Kd = Kp * Tv;
+
+ y = 0;
+ e = 0;
+ esum = 0;
+ //////////////////////////////////////////////////////////////////////////////////////////////////
+ //temporere Lösung
+ std::string fNameCfg = path + "/forcing/forcing.cfg";
+ std::ifstream istr2;
+ istr2.open(fNameCfg.c_str(), std::ios_base::in);
+ if (istr2)
+ {
+ istr2 >> forcing;
+ istr2 >> esum;
+ istr2 >> eold;
+ }
+ istr2.close();
+ ////////////////////////////////////////////////////////////////////////////////////////////////////////
}
//////////////////////////////////////////////////////////////////////////
AdjustForcingCoProcessor::~AdjustForcingCoProcessor()
@@ -60,6 +87,22 @@ void AdjustForcingCoProcessor::process(double step)
//////////////////////////////////////////////////////////////////////////
void AdjustForcingCoProcessor::collectData(double step)
{
+ //////////////////////////////////////////////////////////////////////////////////////////////////
+ //temporere Lösung
+ std::string fNameCfg = path + "/forcing/forcing.cfg";
+ std::ofstream ostr2;
+ ostr2.open(fNameCfg.c_str(), std::ios_base::out);
+ if (!ostr2)
+ {
+ ostr2.clear();
+ string path = UbSystem::getPathFromString(fNameCfg);
+ if (path.size() > 0) { UbSystem::makeDirectory(path); ostr2.open(fNameCfg.c_str(), std::ios_base::out); }
+ if (!ostr2) throw UbException(UB_EXARGS, "couldn't open file " + fNameCfg);
+ }
+ ostr2 << forcing << " " << esum << " " << eold;
+ ostr2.close();
+ ////////////////////////////////////////////////////////////////////////////////////////////////////////
+
integrateValues->calculateMQ();
UBLOG(logDEBUG3, "D3Q27AdjustForcingCoProcessor::update:" << step);
@@ -72,24 +115,17 @@ void AdjustForcingCoProcessor::collectData(double step)
double vx1 = integrateValues->getVx1();
double vx1Average = (vx1/cellsVolume);
- //double C = 5.0; //0.7 //P; //free parameter [0.1, 1]
- //double factor = 1.0 - ((vx1average - vTarged) / vTarged)*C;
-
- double factor = vTarged / vx1Average;
- //forcing *= factor;
-
- //forcing = UbMath::max(forcing, 0.0);
- //forcing = UbMath::min(forcing, 5e-3);
+//////////////////////////////////////////////////////////////////////////
+//PID-Controller (PID-Regler)
+ e = vx1Targed - vx1Average;
+ esum = esum + e;
+ y = Kp * e + Ki * Ta * esum + Kd * (e - eold) / Ta;
+ eold = e;
- if (vTarged > vx1Average)
- {
- forcing = fabs(factor*forcing);
- }
- else
- {
- forcing = -fabs(factor*forcing);
- }
+ forcing = forcing + y;
+//////////////////////////////////////////////////////////////////////////
+ //////////////////////////////////////////////////////////////////////////
mu::Parser fctForcingX1, fctForcingX2, fctForcingX3;
fctForcingX1.SetExpr("Fx1");
fctForcingX1.DefineConst("Fx1", forcing);
@@ -124,7 +160,9 @@ void AdjustForcingCoProcessor::collectData(double step)
if (!ostr) throw UbException(UB_EXARGS, "couldn't open file "+fname);
}
int istep = static_cast<int>(step);
- ostr << istep << ";" << cellsVolume << ";" << vx1Average << "; " << factor << "; " << forcing << "\n";
+
+ ostr << istep << ";" << cellsVolume << ";" << vx1Average << "; " << forcing << "\n";
ostr.close();
+
}
}
diff --git a/source/VirtualFluidsCore/CoProcessors/AdjustForcingCoProcessor.h b/source/VirtualFluidsCore/CoProcessors/AdjustForcingCoProcessor.h
index c7e99ecf4..8eb785843 100644
--- a/source/VirtualFluidsCore/CoProcessors/AdjustForcingCoProcessor.h
+++ b/source/VirtualFluidsCore/CoProcessors/AdjustForcingCoProcessor.h
@@ -10,8 +10,9 @@
class AdjustForcingCoProcessor;
typedef boost::shared_ptr<AdjustForcingCoProcessor> AdjustForcingCoProcessorPtr;
-//! \brief Computes Forcing such that a given velocity (vxZiel) is reached inside an averaging domain (h1).
-//! \details Integrate values helper, scheduler must be set in test case. Example usage: bKanal.cpp
+//! \brief Computes forcing such that a given velocity (vx1Targed) is reached inside an averaging domain (h1).
+//! \details Algorithm based on PID controller (proportional–integral–derivative controller). The parameters of PID controller estimation based on Ziegler–Nichols method.
+//! Integrate values helper, scheduler must be set in test case.
//! \author: Konstantin Kutscher
class AdjustForcingCoProcessor: public CoProcessor {
@@ -30,9 +31,22 @@ protected:
void collectData(double step);
CommunicatorPtr comm;
private:
- double vTarged; //!< target velocity.
+ double vx1Targed; //!< target velocity.
double forcing; //!< forcing at previous update step.
- std::vector<CalcNodes> cnodes;
+
+ double Kpcrit; //Kp critical
+ double Tcrit; //the oscillation period
+ double Tn;
+ double Tv;
+ double e;
+ double Ta;
+ double Kp;
+ double Ki;
+ double Kd;
+ double y;
+ double esum;
+ double eold;
+ //std::vector<CalcNodes> cnodes;
std::string path;
};
diff --git a/source/VirtualFluidsCore/CoProcessors/TimeAveragedValuesCoProcessor.cpp b/source/VirtualFluidsCore/CoProcessors/TimeAveragedValuesCoProcessor.cpp
index 910ec8721..294b6464e 100644
--- a/source/VirtualFluidsCore/CoProcessors/TimeAveragedValuesCoProcessor.cpp
+++ b/source/VirtualFluidsCore/CoProcessors/TimeAveragedValuesCoProcessor.cpp
@@ -112,7 +112,7 @@ void TimeAveragedValuesCoProcessor::process(double step)
if (scheduler->isDue(step))
{
//DEBUG/////////////////////
- UBLOG(logINFO, "step = " << step << ", lcounter = " << lcounter << ", fineStep = " << fineStep << ", maxFineStep = " << maxFineStep << ", numberOfFineSteps = " << numberOfFineSteps);
+ //UBLOG(logINFO, "step = " << step << ", lcounter = " << lcounter << ", fineStep = " << fineStep << ", maxFineStep = " << maxFineStep << ", numberOfFineSteps = " << numberOfFineSteps);
////////////////////////////
calculateSubtotal();
--
GitLab