VM2D 1.14
Vortex methods for 2D flows simulation
Loading...
Searching...
No Matches
VM2D::Velocity Class Referenceabstract

Абстрактный класс, определяющий способ вычисления скоростей More...

#include <Velocity2D.h>

Inheritance diagram for VM2D::Velocity:
Inheritance graph
[legend]
Collaboration diagram for VM2D::Velocity:
Collaboration graph
[legend]

Public Member Functions

 Velocity (const World2D &W_)
 Конструктор
void CalcConvVelo ()
 Вычисление конвективных скоростей вихрей и виртуальных вихрей в вихревом следе, а также в точках wakeVP.
virtual float CalcConvVeloToSetOfPointsFromWake (const WakeDataBase &pointsDb, std::vector< Point2D > &velo, std::vector< double > &domainRadius, bool calcVelo, bool calcRadius)=0
 Вычисление конвективных скоростей и радиусов вихревых доменов в заданном наборе точек от следа
virtual void CalcConvVPVeloToSetOfPointsFromWake (const WakeDataBase &pointsDb, std::vector< Point2D > &velo, std::vector< double > &domainRadius, bool calcVelo, bool calcRadius)
virtual void CalcVeloToWakeVP ()
 Вычисление скоростей в точках wakeVP.
void CalcDiffVeloI1I2ToSetOfPointsFromWake (const WakeDataBase &pointsDb, const std::vector< double > &domainRadius, const WakeDataBase &vorticesDb, std::vector< double > &I1, std::vector< Point2D > &I2)
 Вычисление числителей и знаменателей диффузионных скоростей в заданном наборе точек
void CalcDiffVeloI1I2ToWakeFromWakeFast (const WakeDataBase &wakeDb, const std::vector< double > &domainRadius, std::vector< double > &I1, std::vector< Point2D > &I2)
 Быстрый (с деревом поиска соседей) расчет I1 и I2 от вихрей следа в точках этого же следа
void CalcDiffVeloI1I2ToSetOfPointsFromSheets (const WakeDataBase &pointsDb, const std::vector< double > &domainRadius, const Boundary &bnd, std::vector< double > &I1, std::vector< Point2D > &I2)
void CalcOmegaFromWakeFast (const WakeDataBase &wakeDb, std::vector< double > &Omega)
void CalcDiffVeloI1I2ToWakeFromSheets (const WakeDataBase &pointsDb, const std::vector< double > &domainRadius, const Boundary &bnd, std::vector< double > &I1, std::vector< Point2D > &I2)
void CalcDiffVeloI1I2ToWakeFromWake (const WakeDataBase &pointsDb, const std::vector< double > &domainRadius, const WakeDataBase &vorticesDb, std::vector< double > &I1, std::vector< Point2D > &I2)
void CalcDiffVeloI1I2 ()
 Вычисление диффузионных скоростей вихрей и виртуальных вихрей в вихревом следе
void CalcDiffVeloI0I3 ()
void LimitDiffVelo (std::vector< Point2D > &diffVel)
 Контроль больших значений диффузионных скоростей
void CalcDiffVelo ()
 Вычисление диффузионных скоростей
void ResizeAndZero ()
 Очистка старых массивов под хранение скоростей, выделение новой памяти и обнуление
void SaveVisStress ()
 Сохранение вязких напряжений
void GetWakeInfluenceToRhs (const Airfoil &afl, std::vector< double > &wakeRhs) const
 Генерация вектора влияния вихревого следа на профиль
void CPUGetFASTWakeInfluenceToRhs (std::vector< double > &wakeRhs, std::vector< double > &wakeRhsLin) const
void FillRhs (Eigen::VectorXd &rhsReord) const
virtual ~Velocity ()
 Деструктор

Public Attributes

VortexesParams wakeVortexesParams
 Струтура, определяющая параметры вихрей в следе
std::vector< VortexesParams > virtualVortexesParams
 Вектор струтур, определяющий параметры виртуальных вихрей для профилей
OptimizedVelocity optimizedVelocity

Protected Attributes

const World2D & W
 Константная ссылка на решаемую задачу

Detailed Description

Абстрактный класс, определяющий способ вычисления скоростей

Author
Марчевский Илья Константинович
Сокол Ксения Сергеевна
Рятина Евгения Павловна
Колганова Александра Олеговна
Version
1.14
Date
6 марта 2026 г.

Definition at line 104 of file Velocity2D.h.

Constructor & Destructor Documentation

◆ Velocity()

VM2D::Velocity::Velocity ( const World2D & W_)
inline

Конструктор

Parameters
[in]W_константная ссылка на решаемую задачу

Definition at line 122 of file Velocity2D.h.

122 :
123 W(W_)
124 {
125 virtualVortexesParams.resize(0);
126 };
VM2D::Velocity::virtualVortexesParams
std::vector< VortexesParams > virtualVortexesParams
Вектор струтур, определяющий параметры виртуальных вихрей для профилей
Definition Velocity2D.h:115
VM2D::Velocity::W
const World2D & W
Константная ссылка на решаемую задачу
Definition Velocity2D.h:108
Here is the caller graph for this function:

◆ ~Velocity()

virtual VM2D::Velocity::~Velocity ( )
inlinevirtual

Деструктор

Definition at line 224 of file Velocity2D.h.

224{ };

Member Function Documentation

◆ CalcConvVelo()

void Velocity::CalcConvVelo ( )

Вычисление конвективных скоростей вихрей и виртуальных вихрей в вихревом следе, а также в точках wakeVP.

Warning
скорости приплюсовываются к уже имеющимся

FMM ///

Definition at line 1218 of file Velocity2D.cpp.

1219{
1220 //VMlib::vmTimer Tmr1("a"), Tmr2("b"), Tmr3("c"), Tmr4("d"), Tmr5("e"), Tmr6("f");
1221
1222 W.getTimers().start("ConvVel");
1223
1224 //Влияние следа на след
1225#if (defined(__CUDACC__) || defined(USE_CUDA)) && (defined(CU_CONV_TOWAKE))
1226 //Вызов виртуальной функции
1227
1228 W.timerConvVelo.reset();
1229 W.timerConvVelo.start();
1230
1231 VMlib::vmTimer timerAgpu, timerBgpu;
1232
1233 //Tmr1.start();
1234 std::unique_ptr<BHcu::CudaTreeInfo>& cntrTree = W.getNonConstCuda().cntrTreeWake;
1235 auto& treeWake = *W.getNonConstCuda().inflTreeWake;
1236 if (W.getPassport().numericalSchemes.velocityComputation.second == 1)
1237 {
1238 cntrTree->MemoryAllocate((int)W.getCuda().n_CUDA_wake);
1239 cntrTree->Update((int)W.getWake().vtx.size(), W.getWake().devVtxPtr);
1240 cntrTree->Build();
1241
1242 //Перестроение дерева для вихрей //todo //временно для корректного расчета epsast
1243 treeWake.MemoryAllocate((int)W.getCuda().n_CUDA_wake);
1244 treeWake.Update((int)W.getWake().vtx.size(), W.getWake().devVtxPtr);
1245 timerAgpu.start();
1246 treeWake.Build();
1247 timerAgpu.stop();
1248
1249 timerBgpu.start();
1250 treeWake.UpwardTraversal(W.getPassport().numericalSchemes.nbodyMultipoleOrder);
1251 timerBgpu.stop();
1252 }
1253 //Tmr1.stop();
1254
1255 //Tmr2.start();
1256 double timerCgpuduration = GPUCalcConvVeloToSetOfPointsFromWake(W.getNonConstCuda().cntrTreeWake, W.getWake(), wakeVortexesParams.convVelo, wakeVortexesParams.epsastWake, true, true);
1257 //Tmr2.stop();
1258
1259 W.timerConvVelo.stop();
1260
1262 /*
1263 std::vector<fmm::particle2d> sources(W.getWake().vtx.size());
1264 #pragma omp parallel for
1265 for (int i = 0; i < (int)sources.size(); ++i)
1266 {
1267 const auto& w = W.getWake().vtx[i];
1268
1269 sources[i].center = fmm::point2d(w.r()[0], w.r()[1]);
1270 sources[i].q = w.g();
1271 sources[i].morton_code = 0;
1272 }
1273
1274 auto targets = sources;
1275 for (int i = 0; i < (int)sources.size(); ++i)
1276 targets[i].q = 0;
1277
1278 double t = omp_get_wtime();
1279 //fmm::ComputeExact<fmm::InteractionType::Potential2d>(sources, targets); // or fmm::gpu::ComputeExact... for gpu version
1280 fmm::ComputeExact<fmm::InteractionType::Force2d>(sources, targets); // or fmm::gpu::ComputeExact... for gpu version
1281 std::cout << "Compute exact time = " << omp_get_wtime() - t << std::endl;
1282
1283 t = omp_get_wtime();
1284
1285 fmm::FastMultipole fmm(sources, targets, 1.e-8, fmm::FMM_AUTO, fmm::FMM_AUTO); // or fmm::gpu::FastMultipole... for gpu version
1286 fmm::ReadError<fmm::InteractionType::Force2d>(fmm.forces);
1287 //fmm::ReadError<fmm::InteractionType::Potential2d>(fmm.potentials);
1288 std::cout << "Total time = " << omp_get_wtime() - t << std::endl;
1289 */
1290
1291
1292
1294 /*
1295 std::ofstream treeTimeFile;
1296 if (W.getCurrentStep() == 0)
1297 {
1298 treeTimeFile.open(W.getPassport().dir + "/dbg/treeTime.csv");
1299 treeTimeFile << "step,time,N,lev,averPoints,tBUI,tUPW,tDNW,tBUIgpu,tUPWgpu,tDNWgpu\n";
1300 }
1301 else
1302 treeTimeFile.open(W.getPassport().dir + "/dbg/treeTime.csv", std::ios::app);
1303
1304
1305 int estLevel = std::max(4, (int)(log2(W.getWake().vtx.size())));
1306
1307 for (int controlLevel = std::max(estLevel - 2, 1); controlLevel <= estLevel + 2; ++controlLevel)
1308 {
1309 VMlib::vmTimer timerA, timerB;
1310 if (W.getWake().vtx.size() > 0)
1311 {
1312 if (W.getPassport().numericalSchemes.velocityComputation.second == 1)
1313 {
1314 W.getInflTreeWake().Update(W.getWake().vtx, controlLevel);
1315 timerA.start();
1316 W.getInflTreeWake().Build();
1317 timerA.stop();
1318 timerB.start();
1319 W.getInflTreeWake().UpwardTraversal(W.getPassport().numericalSchemes.nbodyMultipoleOrder);
1320 timerB.stop();
1321
1322 W.getCntrTreeWake().Update(W.getWake().vtx, controlLevel);
1323 W.getCntrTreeWake().Build();
1324 W.getCntrTreeWake().UpwardTraversal(W.getPassport().numericalSchemes.nbodyMultipoleOrder);
1325 }
1326
1327 std::vector<Point2D> tempVelo(wakeVortexesParams.convVelo.size());
1328 std::vector<double> tempRadius(wakeVortexesParams.epsastWake.size());
1329
1330 //Вызов виртуальной функции
1331 double timerCduration = CalcConvVeloToSetOfPointsFromWake(W.getWake(), tempVelo, tempRadius, true, false);
1332
1333 int nObjInContr = 0;
1334 for (int c = 0; c < W.getCntrTreeWake().indexControlCells.size(); ++c)
1335 {
1336 const auto& idx = W.getCntrTreeWake().indexControlCells[c];
1337 if (idx > W.getCntrTreeWake().object.size())
1338 nObjInContr += 1;
1339 else
1340 {
1341 nObjInContr += W.getCntrTreeWake().range[idx].second - W.getCntrTreeWake().range[idx].first + 1;
1342 }
1343 }
1344
1345 treeTimeFile << W.getCurrentStep() << ',' << W.getCurrentTime() << ',' << W.getInflTreeWake().object.size() << ',' \
1346 << controlLevel << ',' << (double)nObjInContr / W.getCntrTreeWake().indexControlCells.size() << ',' \
1347 << timerA.duration() << ',' << timerB.duration() << ',' << timerCduration << ',' \
1348 << timerAgpu.duration() << ',' << timerBgpu.duration() << ',' << timerCgpuduration << '\n';
1349 }
1350 }
1351 treeTimeFile.close();
1352 //*/
1353#else
1354
1355 W.timerConvVelo.reset();
1356 W.timerConvVelo.start();
1357 if (W.getWake().vtx.size() > 0)
1358 {
1359 if (W.getPassport().numericalSchemes.velocityComputation.second == 1)
1360 {
1361
1362 //KFU
1363 //более глубокий уровень контрольных ячеек заметно ускоряет обход дерева на CPU при той же точности
1364 int estOptLevel = std::max(4, (int)(log2(W.getWake().vtx.size())));
1365 //контрольные ячейки по числу вихрей (не более 128), а не по уровню: иначе по мере удлинения следа
1366 //ячейки в плотной области у тела укрупняются и обход дерева сильно замедляется
1367 W.getInflTreeWake().controlCellSize = 128;
1368
1369 //Tmr1.start();
1370 W.getInflTreeWake().Update(W.getWake().vtx, estOptLevel);
1371 W.getInflTreeWake().Build();
1372 W.getInflTreeWake().UpwardTraversal(W.getPassport().numericalSchemes.nbodyMultipoleOrder);
1373 //Tmr1.stop();
1374 }
1375
1376 //Вызов виртуальной функции
1377 //Tmr2.start();
1378 CalcConvVeloToSetOfPointsFromWake(W.getWake(), wakeVortexesParams.convVelo, wakeVortexesParams.epsastWake, true, true);
1379 //Tmr2.stop();
1380
1381 //std::string vname = VMlib::fileNameStep("VeloBH", 0, wakeVortexesParams.convVelo.size(), "txt");
1382 //std::ofstream velofile;
1383 //velofile.open(W.getPassport().dir + vname);
1384
1385 //for (int i = 0; i < wakeVortexesParams.convVelo.size(); ++i)
1386 // velofile << wakeVortexesParams.convVelo[i][0] << " " << wakeVortexesParams.convVelo[i][1] << "\n";
1387 //velofile.close();
1388 }
1389 W.timerConvVelo.stop();
1390
1391
1393 /*
1394 VMlib::vmTimer timerExact;
1395
1396 if (W.getWake().vtx.size() > 0)
1397 {
1398 std::vector<fmm::particle2d> sources(W.getWake().vtx.size());
1399 #pragma omp parallel for
1400 for (int i = 0; i < (int)sources.size(); ++i)
1401 {
1402 const auto& w = W.getWake().vtx[i];
1403
1404 sources[i].center = fmm::point2d(w.r()[0], w.r()[1]);
1405 sources[i].q = w.g();
1406 sources[i].morton_code = 0;
1407 }
1408
1409 auto targets = sources;
1410 for (int i = 0; i < (int)sources.size(); ++i)
1411 targets[i].q = 0;
1412
1413 //double t = omp_get_wtime();
1414 //fmm::ComputeExact<fmm::InteractionType::Potential2d>(sources, targets); // or fmm::gpu::ComputeExact... for gpu version
1415 timerExact.start();
1416 fmm::ComputeExact<fmm::InteractionType::Force2d>(sources, targets); // or fmm::gpu::ComputeExact... for gpu version
1417 timerExact.stop();
1418 //std::cout << "Compute exact time = " << omp_get_wtime() - t << std::endl;
1419
1420 //t = omp_get_wtime();
1421
1422 fmm::FastMultipole fmm(sources, targets, 1.e-8, fmm::FMM_AUTO, fmm::FMM_AUTO); // or fmm::gpu::FastMultipole... for gpu version
1423 fmm::ReadError<fmm::InteractionType::Force2d>(fmm.forces);
1424 //fmm::ReadError<fmm::InteractionType::Potential2d>(fmm.potentials);
1425 //std::cout << "Total time = " << omp_get_wtime() - t << std::endl;
1426
1427 std::cout << "Hybrid time ConvVelo Wake = " << fullTimer.duration() <<
1428 //W.timerConvVelo.duration() <<
1429 " ( " << timeUpd1 + timeUpd2 + timeBld1 + timeBld2 + timeUpw1 + timeUpw2 + timeDnw << " )" << '\n';
1430
1431 fmm::ReadErrorHyb(*((std::vector<std::array<double, 2>>*) &wakeVortexesParams.convVelo));
1432 std::cout << "Exact time = " << timerExact.duration() << '\n';
1433 }
1434 //*/
1435
1436
1437
1438#endif
1439
1440 //Tmr3.start();
1441
1442 std::vector<Point2D> nullVector(0);
1443 //вычисление конвективных скоростей по закону Био-Савара от виртуальных вихрей (точнее, от слоев)
1444 for (size_t bou = 0; bou < W.getNumberOfBoundary(); ++bou)
1445 {
1446 //Влияние на след от панелей
1447#if (defined(__CUDACC__) || defined(USE_CUDA)) && (defined(CU_CONVVIRT))
1448
1449 switch (W.getPassport().numericalSchemes.velocityComputation.second)
1450 {
1451 case 0:
1452 W.getBoundary(bou).GPUCalcConvVelocityToSetOfPointsFromSheets(W.getWake(), wakeVortexesParams.convVelo);
1453 break;
1454 case 1:
1455 if (bou == 0)
1456 GPUCalcConvVelocityToSetOfPointsFromSheets(W.getNonConstCuda().cntrTreeWake, W.getWake(), wakeVortexesParams.convVelo);
1457 break;
1458 }
1459 //
1460
1461#else
1462 W.getBoundary(bou).CalcConvVelocityToSetOfPointsFromSheets(W.getWake(), wakeVortexesParams.convVelo);
1463#endif
1464 }
1465 //Tmr3.stop();
1466
1467 //Скорости только что рожденных вихрей
1468 //Tmr4.start();
1469 {
1470 size_t nVirtVortices = 0;
1471 for (size_t bou = 0; bou < W.getNumberOfBoundary(); ++bou)
1472 nVirtVortices += W.getBoundary(bou).virtualWake.vtx.size();
1473
1474 size_t counter = wakeVortexesParams.convVelo.size() - nVirtVortices;
1475
1476 for (size_t bou = 0; bou < W.getNumberOfBoundary(); ++bou)
1477 for (size_t v = 0; v < W.getBoundary(bou).virtualWake.vtx.size(); ++v)
1478 {
1479 wakeVortexesParams.convVelo[counter] = W.getBoundary(bou).afl.getV(W.getBoundary(bou).virtualWake.aflPan[v].second) \
1480 + W.getBoundary(bou).virtualWake.vecHalfGamma[v] \
1481 - W.getV0();
1482 ++counter;
1483 }
1484 }
1485 //Tmr4.stop();
1486
1487 //Заполнение структур данных виртуальных вихрей
1488 //Tmr5.start();
1489 {
1490 size_t nVirtVortices = 0;
1491 for (size_t bou = 0; bou < W.getNumberOfBoundary(); ++bou)
1492 nVirtVortices += W.getBoundary(bou).virtualWake.vtx.size();
1493
1494 size_t counter = wakeVortexesParams.convVelo.size() - nVirtVortices;
1495
1496 for (size_t bou = 0; bou < W.getNumberOfBoundary(); ++bou)
1497 for (size_t v = 0; v < W.getBoundary(bou).virtualWake.vtx.size(); ++v)
1498 {
1499 virtualVortexesParams[bou].convVelo[v] = wakeVortexesParams.convVelo[counter];
1500 virtualVortexesParams[bou].epsastWake[v] = wakeVortexesParams.epsastWake[counter];
1501 ++counter;
1502 }
1503 }
1504 //Tmr5.stop();
1505
1506 W.getTimers().stop("ConvVel");
1507 W.getTimers().start("VelPres");
1508
1509 //Tmr6.start();
1510 //вычисление скоростей в заданных точках только на соответствующем шаге по времени
1511 if ((W.getPassport().timeDiscretizationProperties.saveVPstep > 0) && (!(W.getCurrentStep() % W.getPassport().timeDiscretizationProperties.saveVPstep)))
1512#ifdef OPTIMIZER
1513 if ((W.getCurrentStep() >= OPTIMIZER_START_STEP && W.getCurrentStep() <= OPTIMIZER_STOP_STEP) || W.getCurrentStep() == 0)
1514 CalcVeloToWakeVP();
1515#else
1516 CalcVeloToWakeVP();
1517#endif
1518 //Tmr6.stop();
1519
1521 /*if ((W.getPassport().timeDiscretizationProperties.saveVPstep > 0) && (!(W.getCurrentStep() % W.getPassport().timeDiscretizationProperties.saveVPstep)))
1522 {
1523 int currentStep = W.getCurrentStep();
1524 if (currentStep >= 1100 && currentStep <= 1300)
1525 {
1526 CalcVeloToWakeVP();
1527 optimizedVelocity.AddVelocities(W.getNonConstMeasureVP().getNonConstVelocity());
1528 }
1529
1530 if (currentStep == 1200)
1531 {
1532 optimizedVelocity.PerformAveragingVelo(W.getPassport().dir);
1533 }
1534 }*/
1536
1537 //std::cout << "ConvVelo: " << Tmr1.duration() << " " << Tmr2.duration() << " " << Tmr3.duration() << " "
1538 // << Tmr4.duration() << " " << Tmr5.duration() << " " << Tmr6.duration() << "\n";
1539
1540 W.getTimers().stop("VelPres");
1541}//CalcConvVelo()
VM2D::Velocity::CalcConvVeloToSetOfPointsFromWake
virtual float CalcConvVeloToSetOfPointsFromWake(const WakeDataBase &pointsDb, std::vector< Point2D > &velo, std::vector< double > &domainRadius, bool calcVelo, bool calcRadius)=0
Вычисление конвективных скоростей и радиусов вихревых доменов в заданном наборе точек от следа
VM2D::Velocity::wakeVortexesParams
VortexesParams wakeVortexesParams
Струтура, определяющая параметры вихрей в следе
Definition Velocity2D.h:112
VM2D::Velocity::CalcVeloToWakeVP
virtual void CalcVeloToWakeVP()
Вычисление скоростей в точках wakeVP.
Definition Velocity2D.cpp:1546
VMlib::vmTimer::stop
const vmTimer & stop() const
Останов работающего счетчика времени
Definition TimesGen.h:125
VMlib::vmTimer::start
const vmTimer & start() const
Запуск (первый или повторный) счетчика времени
Definition TimesGen.h:109
Here is the call graph for this function:

◆ CalcConvVeloToSetOfPointsFromWake()

virtual float VM2D::Velocity::CalcConvVeloToSetOfPointsFromWake ( const WakeDataBase & pointsDb,
std::vector< Point2D > & velo,
std::vector< double > & domainRadius,
bool calcVelo,
bool calcRadius )
pure virtual

Вычисление конвективных скоростей и радиусов вихревых доменов в заданном наборе точек от следа

Parameters
[in]pointsDbконстантная ссылка на базу данных пелены из вихрей, в которых надо сосчитать конвективные скорости
[out]veloссылка на вектор скоростей в требуемых точках
[out]domainRadiusссылка на вектор радиусов вихревых доменов
[in]calcVeloпризнак вычисления скоростей в точках
[in]calcRadiusпризнак вычисления радиусов доменов

Implemented in VM2D::VelocityBarnesHut, and VM2D::VelocityBiotSavart.

Here is the caller graph for this function:

◆ CalcConvVPVeloToSetOfPointsFromWake()

virtual void VM2D::Velocity::CalcConvVPVeloToSetOfPointsFromWake ( const WakeDataBase & pointsDb,
std::vector< Point2D > & velo,
std::vector< double > & domainRadius,
bool calcVelo,
bool calcRadius )
inlinevirtual

Reimplemented in VM2D::VelocityBarnesHut.

Definition at line 141 of file Velocity2D.h.

141{};

◆ CalcDiffVelo()

void Velocity::CalcDiffVelo ( )

Вычисление диффузионных скоростей

Вызывается в CalcVortexVelo()

Definition at line 185 of file Velocity2D.cpp.

186{
187 W.getTimers().start("DiffVel");
188
189 if (W.getPassport().physicalProperties.nu > 0.0)
190 {
191 //VMlib::vmTimer D1, D2, D3, D4, D5;
192
193 //D1.start();
194 CalcDiffVeloI1I2();
195 //D1.stop();
196
197 //D2.start();
198 CalcDiffVeloI0I3();
199 //D2.stop();
200
201 //D3.start();
202 for (Point2D& diffV : wakeVortexesParams.diffVelo)
203 diffV *= W.getPassport().physicalProperties.nu;
204
205 //контроль застрелов диффузионной скорости
206 LimitDiffVelo(wakeVortexesParams.diffVelo);
207 //D3.stop();
208
209 //D4.start();
210 for (size_t afl = 0; afl < W.getNumberOfAirfoil(); ++afl)
211 for (size_t i = 0; i < W.getAirfoil(afl).viscousStress.size(); ++i)
212 W.getNonConstAirfoil(afl).viscousStress[i] *= W.getPassport().physicalProperties.nu;
213
214 SaveVisStress();
215 //D4.stop();
216
217 //D5.start();
218 //Заполнение структуры данных виртуальных вихрей
219 size_t nVirtVortices = 0;
220 for (size_t bou = 0; bou < W.getNumberOfBoundary(); ++bou)
221 nVirtVortices += W.getBoundary(bou).virtualWake.vtx.size();
222
223 size_t counter = wakeVortexesParams.convVelo.size() - nVirtVortices;
224
225 for (size_t bou = 0; bou < W.getNumberOfBoundary(); ++bou)
226 for (size_t v = 0; v < W.getBoundary(bou).virtualWake.vtx.size(); ++v)
227 {
228 virtualVortexesParams[bou].diffVelo[v] = wakeVortexesParams.diffVelo[counter];
229 virtualVortexesParams[bou].I0[v] = wakeVortexesParams.I0[counter];
230 virtualVortexesParams[bou].I1[v] = wakeVortexesParams.I1[counter];
231 virtualVortexesParams[bou].I2[v] = wakeVortexesParams.I2[counter];
232 virtualVortexesParams[bou].I3[v] = wakeVortexesParams.I3[counter];
233 ++counter;
234 }
235 //D5.stop();
236
237 //std::cout << "Diff: " << D1.duration() << " " << D2.duration() << " " << D3.duration() << " "
238 // << D4.duration() << " " << D5.duration() << "\n";
239
240 }
241 W.getTimers().stop("DiffVel");
242
243}// CalcDiffVelo()
VM2D::Velocity::LimitDiffVelo
void LimitDiffVelo(std::vector< Point2D > &diffVel)
Контроль больших значений диффузионных скоростей
Definition Velocity2D.cpp:173
VM2D::Velocity::CalcDiffVeloI0I3
void CalcDiffVeloI0I3()
Definition Velocity2D.cpp:115
VM2D::Velocity::SaveVisStress
void SaveVisStress()
Сохранение вязких напряжений
Definition Velocity2D.cpp:831
VM2D::Velocity::CalcDiffVeloI1I2
void CalcDiffVeloI1I2()
Вычисление диффузионных скоростей вихрей и виртуальных вихрей в вихревом следе
Definition Velocity2D.cpp:64
Here is the call graph for this function:

◆ CalcDiffVeloI0I3()

void Velocity::CalcDiffVeloI0I3 ( )

Definition at line 115 of file Velocity2D.cpp.

116{
117 //VMlib::vmTimer T03_1, T03_2, T03_3;
118
119 //T03_1.start();
120
121 for (size_t afl = 0; afl < W.getNumberOfAirfoil(); ++afl)
122 {
123#if (defined(__CUDACC__) || defined(USE_CUDA)) && (defined(CU_I0I3))
124 W.getNonConstAirfoil(afl).GPUGetDiffVelocityI0I3ToSetOfPointsAndViscousStresses(wakeVortexesParams.epsastWake, wakeVortexesParams.I0, wakeVortexesParams.I3);
125#else
126 W.getNonConstAirfoil(afl).GetDiffVelocityI0I3ToSetOfPointsAndViscousStresses(W.getWake(), wakeVortexesParams.epsastWake, wakeVortexesParams.I0, wakeVortexesParams.I3);
127#endif
128 }
129
130 //T03_1.stop();
131
132
133 //Слагаемое, связанное с движением стенки: t_tau = mu * (omega_w + 2 n * dU_w/ds),
134 //dU_w/ds на панели j --- (V_{j+1} - V_j) / len[j]; для вращения тела как целого оно дает поправку -2 mu Omega
135 //T03_2.start();
136 for (size_t afl = 0; afl < W.getNumberOfAirfoil(); ++afl)
137 {
138 auto& A = W.getNonConstAirfoil(afl);
139 for (size_t j = 0; j < A.getNumberOfPanels(); ++j)
140 A.viscousStress[j] += 2.0 * (A.nrm[j] & (A.getV(j + 1) - A.getV(j)));
141 }
142 //T03_2.stop();
143
144 //влияние поверхности
145 //T03_3.start();
146 Point2D I3;
147 double I0;
148
149 double domrad = 0.0;
150
151#pragma omp parallel for private(I0, I3, domrad)
152 for (int vt = 0; vt < (int)wakeVortexesParams.diffVelo.size(); ++vt)
153 {
154 domrad = std::max(wakeVortexesParams.epsastWake[vt], W.getPassport().wakeDiscretizationProperties.getMinEpsAst());
155
156 wakeVortexesParams.I0[vt] *= domrad;
157 wakeVortexesParams.I0[vt] += DPI * sqr(domrad);
158
159
160 I3 = wakeVortexesParams.I3[vt];
161 I0 = wakeVortexesParams.I0[vt];
162
163
164 if (fabs(I0) > 1.e-8)
165 wakeVortexesParams.diffVelo[vt] += I3 * (1.0 / I0);
166 }
167 //T03_3.stop();
168
169 //std::cout << "Diff03: " << T03_1.duration() << " " << T03_2.duration() << " " << T03_3.duration() << "\n";
170}//CalcDiffVeloI0I3()
sqr
T sqr(T x)
Возведение числа в квадрат
Definition defs.h:463
DPI
const double DPI
Число .
Definition defs.h:82
fmm::vt
void vt(const Vector3d &x, const Matrix3d *T, Matrix3d &M)
Definition multipole3d.cpp:454
Here is the call graph for this function:
Here is the caller graph for this function:

◆ CalcDiffVeloI1I2()

void Velocity::CalcDiffVeloI1I2 ( )

Вычисление диффузионных скоростей вихрей и виртуальных вихрей в вихревом следе

Вызывает 4 раза функцию CalcDiffVeloToSetOfPoints

Warning
скорости приплюсовываются к уже имеющимся

Definition at line 64 of file Velocity2D.cpp.

65{
66 // !!! пелена на пелену
67#if (defined(__CUDACC__) || defined(USE_CUDA)) && (defined(CU_I1I2))
68 if (W.getPassport().numericalSchemes.velocityComputation.second == 0)
69 GPUCalcDiffVeloI1I2ToSetOfPointsFromWake(W.getWake(), wakeVortexesParams.epsastWake, W.getWake(), wakeVortexesParams.I1, wakeVortexesParams.I2);
70 else
71 GPUDiffVeloFAST(wakeVortexesParams.epsastWake, wakeVortexesParams.I1, wakeVortexesParams.I2);
72#else
73
74 //CalcDiffVeloI1I2ToWakeFromWake(W.getWake(), wakeVortexesParams.epsastWake, W.getWake(), wakeVortexesParams.I1, wakeVortexesParams.I2);
75 //
76 // KFU
77 if (W.getPassport().numericalSchemes.velocityComputation.second == 0)
78 CalcDiffVeloI1I2ToWakeFromWake(W.getWake(), wakeVortexesParams.epsastWake, W.getWake(), wakeVortexesParams.I1, wakeVortexesParams.I2);
79 else
80 CalcDiffVeloI1I2ToWakeFromWakeFast(W.getWake(), wakeVortexesParams.epsastWake, wakeVortexesParams.I1, wakeVortexesParams.I2);
81#endif
82
83
84 for (size_t bou = 0; bou < W.getNumberOfBoundary(); ++bou)
85 {
86 //не нужно, т.к. сделано выше перед началом вычисления скоростей
87 //W.getNonConstBoundary(bou).virtualWake.WakeSynchronize();
88
89
90 //виртуальные на границе на след
91#if (defined(__CUDACC__) || defined(USE_CUDA)) && (defined(CU_I1I2))
92 GPUCalcDiffVeloI1I2ToSetOfPointsFromSheets(W.getWake(), wakeVortexesParams.epsastWake, W.getBoundary(bou), wakeVortexesParams.I1, wakeVortexesParams.I2);
93#else
94 CalcDiffVeloI1I2ToWakeFromSheets(W.getWake(), wakeVortexesParams.epsastWake, W.getBoundary(bou), wakeVortexesParams.I1, wakeVortexesParams.I2);
95#endif
96 } //for bou
97
98 Point2D I2;
99 double I1;
100
101#pragma omp parallel for private(I1, I2)
102 for (int vt = 0; vt < (int)wakeVortexesParams.diffVelo.size(); ++vt)
103 {
104 I2 = wakeVortexesParams.I2[vt];
105 I1 = wakeVortexesParams.I1[vt];
106
107 if (fabs(I1) < 1.e-8)
108 wakeVortexesParams.diffVelo[vt] = { 0.0, 0.0 };
109 else
110 wakeVortexesParams.diffVelo[vt] = I2 * (1.0 / (I1 * std::max(wakeVortexesParams.epsastWake[vt], W.getPassport().wakeDiscretizationProperties.getMinEpsAst())));
111 }
112}//CalcDiffVeloI1I2()
VM2D::Velocity::CalcDiffVeloI1I2ToWakeFromSheets
void CalcDiffVeloI1I2ToWakeFromSheets(const WakeDataBase &pointsDb, const std::vector< double > &domainRadius, const Boundary &bnd, std::vector< double > &I1, std::vector< Point2D > &I2)
Definition Velocity2D.cpp:1664
VM2D::Velocity::CalcDiffVeloI1I2ToWakeFromWakeFast
void CalcDiffVeloI1I2ToWakeFromWakeFast(const WakeDataBase &wakeDb, const std::vector< double > &domainRadius, std::vector< double > &I1, std::vector< Point2D > &I2)
Быстрый (с деревом поиска соседей) расчет I1 и I2 от вихрей следа в точках этого же следа
Definition Velocity2D.cpp:306
VM2D::Velocity::CalcDiffVeloI1I2ToWakeFromWake
void CalcDiffVeloI1I2ToWakeFromWake(const WakeDataBase &pointsDb, const std::vector< double > &domainRadius, const WakeDataBase &vorticesDb, std::vector< double > &I1, std::vector< Point2D > &I2)
Definition Velocity2D.cpp:1659
Here is the call graph for this function:
Here is the caller graph for this function:

◆ CalcDiffVeloI1I2ToSetOfPointsFromSheets()

void Velocity::CalcDiffVeloI1I2ToSetOfPointsFromSheets ( const WakeDataBase & pointsDb,
const std::vector< double > & domainRadius,
const Boundary & bnd,
std::vector< double > & I1,
std::vector< Point2D > & I2 )
Todo
Сделать переменной и синхронизировать с GPU
Todo
Учитываем пока только нулевой момент решения
Todo
Понять природу магической константы 8.0 и синхронизировать с GPU

Definition at line 449 of file Velocity2D.cpp.

450{
451 double tCPUSTART, tCPUEND;
452
453 tCPUSTART = omp_get_wtime();
454
455 std::vector<double> selfI1(pointsDb.vtx.size(), 0.0);
456 std::vector<Point2D> selfI2(pointsDb.vtx.size(), { 0.0, 0.0 });
457
458#pragma warning (push)
459#pragma warning (disable: 4101)
460 //Локальные переменные для цикла
461 Point2D Rij;
462 double rij, expr;
463 double diffRadius;
464 double left;
465 double right;
466 double posJx;
467 double domRad;
468#pragma warning (pop)
469
470
471//Оптимизированный алгоритм
472//KFU
474 const int nQuadPt = 3;
475
476 //квадратурные точки на панелях и их интенсивности вычисляются один раз
477 const size_t nPnl = bnd.afl.getNumberOfPanels();
478 std::vector<Point2D> quadPt(nPnl * nQuadPt);
479 std::vector<double> quadG(nPnl * nQuadPt);
480 NbBox quadBox;
481 quadBox.SetEmpty();
482 for (size_t j = 0; j < nPnl; ++j)
483 {
485 const double ptG = bnd.sheets.freeVortexSheet(j, 0) * bnd.afl.len[j] / nQuadPt;
486 for (int q = 0; q < nQuadPt; ++q)
487 {
488 quadPt[j * nQuadPt + q] = bnd.afl.getR(j) + bnd.afl.tau[j] * (q + 0.5) * bnd.afl.len[j] * (1.0 / nQuadPt);
489 quadG[j * nQuadPt + q] = ptG;
490 quadBox.Add(quadPt[j * nQuadPt + q]);
491 }
492 }
493
494#pragma omp parallel for default(none) shared(selfI1, selfI2, domainRadius, pointsDb, quadPt, quadG, quadBox) private(Rij, rij, expr, domRad, diffRadius, left, right, posJx)
495 for (int i = 0; i < pointsDb.vtx.size(); ++i)
496 {
497 const Vortex2D& vtxI = pointsDb.vtx[i];
498
499 domRad = std::max(domainRadius[i], W.getPassport().wakeDiscretizationProperties.getMinEpsAst());
500
502 diffRadius = 8.0 * domRad;
503
504 //точки, удаленные от всех квадратурных точек более чем на diffRadius, вклада не получают
505 if (quadBox.Dist2(vtxI.r()) >= diffRadius * diffRadius)
506 continue;
507
508 left = vtxI.r()[0] - diffRadius;
509 right = vtxI.r()[0] + diffRadius;
510
511 for (size_t jq = 0; jq < quadPt.size(); ++jq)
512 {
513 const Point2D& ptJ = quadPt[jq];
514 posJx = ptJ[0];
515
516 if ((left < posJx) && (posJx < right))
517 {
518 Rij = vtxI.r() - ptJ;
519 rij = Rij.length();
520 if (rij < diffRadius && rij > 1.e-10)
521 {
522 expr = exp(-rij / domRad);
523 selfI2[i] += (quadG[jq] * expr / rij) * Rij;
524 selfI1[i] += quadG[jq] * expr;
525 }
526 }//if (rij>1e-6)
527 }//for jq
528 } // for r
529
530
531
532// Старый "лобовой" алгоритм
533/*
534#pragma omp parallel for default(none) shared(selfI1, selfI2, domainRadius, bnd, pointsDb, std::cout) private(Rij, rij, expr, domRad, diffRadius, left, right, posJx)
535 for (int i = 0; i < pointsDb.vtx.size(); ++i)
536 {
537 const Vortex2D& vtxI = pointsDb.vtx[i];
538
539 domRad = std::max(domainRadius[i], W.getPassport().wakeDiscretizationProperties.getMinEpsAst());
540
542 diffRadius = 8.0 * domRad;
543
544 left = vtxI.r()[0] - diffRadius;
545 right = vtxI.r()[0] + diffRadius;
546
547 for (size_t j = 0; j < bnd.afl.getNumberOfPanels(); ++j)
548 {
550 const int nQuadPt = 3;
551
553 const double ptG = bnd.sheets.freeVortexSheet(j, 0) * bnd.afl.len[j] / nQuadPt;
554
555 for (int q = 0; q < nQuadPt; ++q)
556 {
557 const Point2D& ptJ = bnd.afl.getR(j) + bnd.afl.tau[j] * (q + 0.5) * bnd.afl.len[j] * (1.0 / nQuadPt); // vorticesDb.vtx[j];
558 posJx = ptJ[0];
559
560 if ((left < posJx) && (posJx < right))
561 {
562 Rij = vtxI.r() - ptJ;
563 rij = Rij.length();
564 if (rij < diffRadius && rij > 1.e-10)
565 {
566 expr = exp(-rij / domRad);
567 selfI2[i] += (ptG * expr / rij) * Rij;
568 selfI1[i] += ptG * expr;
569 }
570 }//if (rij>1e-6)
571 }
572 }//for j
573 } // for r
574*/
575
576
577
578 for (size_t i = 0; i < I1.size(); ++i)
579 {
580 I1[i] += selfI1[i];
581 I2[i] += selfI2[i];
582 }
583
584 tCPUEND = omp_get_wtime();
585 //W.getInfo('t') << "DIFF_CPU: " << tCPUEND - tCPUSTART << std::endl;
586}//CalcDiffVeloI1I2ToSetOfPointsFromSheets(...)
VM2D::AirfoilGeometry::len
std::vector< double > len
Длины панелей профиля
Definition Airfoil2D.h:94
VM2D::AirfoilGeometry::getR
const Point2D & getR(size_t q) const
Возврат константной ссылки на вершину профиля
Definition Airfoil2D.h:113
VM2D::AirfoilGeometry::tau
std::vector< Point2D > tau
Касательные к панелям профиля
Definition Airfoil2D.h:91
VM2D::AirfoilGeometry::getNumberOfPanels
size_t getNumberOfPanels() const
Возврат количества панелей на профиле
Definition Airfoil2D.h:163
VM2D::Boundary::afl
const Airfoil & afl
Definition Boundary2D.h:77
VM2D::Boundary::sheets
Sheet sheets
Слои на профиле
Definition Boundary2D.h:96
VM2D::Sheet::freeVortexSheet
const double & freeVortexSheet(size_t n, size_t moment) const
Definition Sheet2D.h:100
VM2D::WakeDataBase::vtx
std::vector< Vortex2D > vtx
Список вихревых элементов
Definition WakeDataBase2D.h:81
VMlib::Vortex2D::r
HD Point2D & r()
Функция для доступа к радиус-вектору вихря
Definition Vortex2D.h:92
numvector::length
P length() const
Вычисление 2-нормы (длины) вектора
Definition numvector.h:374
VM2D::NbBox::Add
void Add(const Point2D &p)
Definition cpuNeighborTree.h:64
VM2D::NbBox::Dist2
double Dist2(const NbBox &b) const
Квадрат расстояния до другого прямоугольника (для пустого — бесконечность).
Definition cpuNeighborTree.h:79
Here is the call graph for this function:
Here is the caller graph for this function:

◆ CalcDiffVeloI1I2ToSetOfPointsFromWake()

void Velocity::CalcDiffVeloI1I2ToSetOfPointsFromWake ( const WakeDataBase & pointsDb,
const std::vector< double > & domainRadius,
const WakeDataBase & vorticesDb,
std::vector< double > & I1,
std::vector< Point2D > & I2 )

Вычисление числителей и знаменателей диффузионных скоростей в заданном наборе точек

Parameters
[in]pointsDbконстантная ссылка на базу данных пелены из вихрей, в которых надо сосчитать диффузионные скорости
[in]domainRadiusконстантная ссылка на вектор радиусов вихревых доменов
[in]vorticesDbконстантная ссылка на на базу данных пелены из вихрей,от которых надо сосчитать влияния на points
[out]I1ссылка на вектор величин I1 (знаменателей в диффузионных скоростях) в требуемых точках
[out]I2ссылка на вектор величин I2 (числителей в диффузионных скоростях) в требуемых точках
Todo
Понять природу магической константы 8.0 и синхронизировать с GPU

Definition at line 247 of file Velocity2D.cpp.

248{
249 std::vector<double> selfI1(pointsDb.vtx.size(), 0.0);
250 std::vector<Point2D> selfI2(pointsDb.vtx.size(), { 0.0, 0.0 });
251
252#pragma warning (push)
253#pragma warning (disable: 4101)
254 //Локальные переменные для цикла
255 Point2D Rij;
256 double rij, expr;
257 double diffRadius, domRad;
258 double left;
259 double right;
260 double posJx;
261#pragma warning (pop)
262
263#pragma omp parallel for default(none) shared(selfI1, selfI2, domainRadius, vorticesDb, pointsDb) private(Rij, rij, expr, diffRadius, domRad, left, right, posJx)
264 for (int i = 0; i < pointsDb.vtx.size(); ++i)
265 {
266 const Vortex2D& vtxI = pointsDb.vtx[i];
267
268 domRad = std::max(domainRadius[i], W.getPassport().wakeDiscretizationProperties.getMinEpsAst());
269
271 diffRadius = 8.0 * domRad;
272
273 left = vtxI.r()[0] - diffRadius;
274 right = vtxI.r()[0] + diffRadius;
275
276 for (size_t j = 0; j < vorticesDb.vtx.size(); ++j)
277 {
278 const Vortex2D& vtxJ = vorticesDb.vtx[j];
279 posJx = vtxJ.r()[0];
280
281 if ((left < posJx) && (posJx < right))
282 {
283 Rij = vtxI.r() - vtxJ.r();
284 rij = Rij.length();
285 if (rij < diffRadius && rij > 1.e-10)
286 {
287 expr = exp(-rij / domRad);
288 selfI2[i] += (vtxJ.g()* expr / rij) * Rij;
289 selfI1[i] += vtxJ.g()*expr;
290 }
291 }//if (rij>1e-6)
292 }//for j
293 } // for r
294
295
296 for (size_t i = 0; i < I1.size(); ++i)
297 {
298 I1[i] += selfI1[i];
299 I2[i] += selfI2[i];
300 }
301}//CalcDiffVeloI1I2ToSetOfPointsFromWake(...)
VMlib::Vortex2D::g
HD double & g()
Функция для доступа к циркуляции вихря
Definition Vortex2D.h:100
Here is the call graph for this function:
Here is the caller graph for this function:

◆ CalcDiffVeloI1I2ToWakeFromSheets()

void Velocity::CalcDiffVeloI1I2ToWakeFromSheets ( const WakeDataBase & pointsDb,
const std::vector< double > & domainRadius,
const Boundary & bnd,
std::vector< double > & I1,
std::vector< Point2D > & I2 )

Definition at line 1664 of file Velocity2D.cpp.

1665{
1666 CalcDiffVeloI1I2ToSetOfPointsFromSheets(pointsDb, domainRadius, bnd, I1, I2);
1667}
VM2D::Velocity::CalcDiffVeloI1I2ToSetOfPointsFromSheets
void CalcDiffVeloI1I2ToSetOfPointsFromSheets(const WakeDataBase &pointsDb, const std::vector< double > &domainRadius, const Boundary &bnd, std::vector< double > &I1, std::vector< Point2D > &I2)
Definition Velocity2D.cpp:449
Here is the call graph for this function:
Here is the caller graph for this function:

◆ CalcDiffVeloI1I2ToWakeFromWake()

void Velocity::CalcDiffVeloI1I2ToWakeFromWake ( const WakeDataBase & pointsDb,
const std::vector< double > & domainRadius,
const WakeDataBase & vorticesDb,
std::vector< double > & I1,
std::vector< Point2D > & I2 )

Definition at line 1659 of file Velocity2D.cpp.

1660{
1661 CalcDiffVeloI1I2ToSetOfPointsFromWake(pointsDb, domainRadius, vorticesDb, I1, I2);
1662}//CalcDiffVeloI1I2ToWakeFromWake(...)
VM2D::Velocity::CalcDiffVeloI1I2ToSetOfPointsFromWake
void CalcDiffVeloI1I2ToSetOfPointsFromWake(const WakeDataBase &pointsDb, const std::vector< double > &domainRadius, const WakeDataBase &vorticesDb, std::vector< double > &I1, std::vector< Point2D > &I2)
Вычисление числителей и знаменателей диффузионных скоростей в заданном наборе точек
Definition Velocity2D.cpp:247
Here is the call graph for this function:
Here is the caller graph for this function:

◆ CalcDiffVeloI1I2ToWakeFromWakeFast()

void Velocity::CalcDiffVeloI1I2ToWakeFromWakeFast ( const WakeDataBase & wakeDb,
const std::vector< double > & domainRadius,
std::vector< double > & I1,
std::vector< Point2D > & I2 )

Быстрый (с деревом поиска соседей) расчет I1 и I2 от вихрей следа в точках этого же следа

Дает тот же результат, что и CalcDiffVeloI1I2ToSetOfPointsFromWake, с точностью до порядка суммирования

Todo
Понять природу магической константы 8.0 и синхронизировать с GPU (см. CalcDiffVeloI1I2ToSetOfPointsFromWake)

Definition at line 306 of file Velocity2D.cpp.

307{
308 const std::vector<Vortex2D>& vtx = wakeDb.vtx;
309 const int n = (int)vtx.size();
310 const double minEpsAst = W.getPassport().wakeDiscretizationProperties.getMinEpsAst();
311
312 NeighborTree tree;
313 tree.Build(n, [&vtx](int i) -> const Point2D& { return vtx[i].r(); });
314
315 std::vector<double> gSorted(n);
316#pragma omp parallel for
317 for (int k = 0; k < n; ++k)
318 gSorted[k] = vtx[tree.idx[k]].g();
319
320#pragma omp parallel
321 {
322 std::vector<int> cand;
323 cand.reserve(256);
324
325#pragma omp for schedule(dynamic, 4)
326 for (int leaf = 0; leaf < tree.nLeaves; ++leaf)
327 {
329 double maxDiffRadius = 0.0;
330 for (int k = tree.LeafBegin(leaf); k < tree.LeafEnd(leaf); ++k)
331 maxDiffRadius = std::max(maxDiffRadius, 8.0 * std::max(domainRadius[tree.idx[k]], minEpsAst));
332
333 tree.CollectLeaves(tree.LeafBox(leaf), maxDiffRadius, cand);
334
335 for (int k = tree.LeafBegin(leaf); k < tree.LeafEnd(leaf); ++k)
336 {
337 const int i = tree.idx[k];
338 const Point2D& posI = tree.pos[k];
339 const double domRad = std::max(domainRadius[i], minEpsAst);
340 const double iDomRad = 1.0 / domRad;
341 const double diffRadius2 = sqr(8.0 * domRad);
342
343 double selfI1 = 0.0;
344 Point2D selfI2 = { 0.0, 0.0 };
345
346 for (int c : cand)
347 {
348 if (tree.LeafBox(c).Dist2(posI) >= diffRadius2)
349 continue;
350
351 for (int m = tree.LeafBegin(c); m < tree.LeafEnd(c); ++m)
352 {
353 const Point2D Rij = posI - tree.pos[m];
354 const double rij2 = Rij.length2();
355 if (rij2 < diffRadius2 && rij2 > 1.e-20)
356 {
357 const double rij = sqrt(rij2);
358 const double expr = gSorted[m] * exp(-rij * iDomRad);
359 selfI2 += (expr / rij) * Rij;
360 selfI1 += expr;
361 }
362 }
363 }
364
365 I1[i] += selfI1;
366 I2[i] += selfI2;
367 }
368 }
369 }
370}//CalcDiffVeloI1I2ToWakeFromWakeFast(...)
VM2D::NeighborTree::idx
std::vector< int > idx
исходные номера точек в порядке сортировки
Definition cpuNeighborTree.h:109
VM2D::NeighborTree::LeafBegin
int LeafBegin(int leaf) const
Definition cpuNeighborTree.h:113
VM2D::NeighborTree::LeafEnd
int LeafEnd(int leaf) const
Definition cpuNeighborTree.h:114
VM2D::NeighborTree::Build
void Build(int n_, F getPos)
Построение дерева; getPos(i) возвращает координаты i-й точки
Definition cpuNeighborTree.h:119
VM2D::NeighborTree::CollectLeaves
void CollectLeaves(const NbBox &q, double R, std::vector< int > &leaves) const
Список листьев, прямоугольники которых находятся на расстоянии меньше R от прямоугольника q.
Definition cpuNeighborTree.h:185
VM2D::NeighborTree::LeafBox
const NbBox & LeafBox(int leaf) const
Definition cpuNeighborTree.h:115
VM2D::NeighborTree::pos
std::vector< Point2D > pos
координаты точек в порядке сортировки
Definition cpuNeighborTree.h:110
numvector::length2
auto length2() const -> typename std::remove_const< typename std::remove_reference< decltype(this->data[0])>::type >::type
Вычисление квадрата нормы (длины) вектора
Definition numvector.h:386
Here is the call graph for this function:
Here is the caller graph for this function:

◆ CalcOmegaFromWakeFast()

void Velocity::CalcOmegaFromWakeFast ( const WakeDataBase & wakeDb,
std::vector< double > & Omega )
Todo
Понять природу магической константы 8.0 и синхронизировать с GPU

Definition at line 376 of file Velocity2D.cpp.

377{
378 const std::vector<Vortex2D>& vtx = wakeDb.vtx;
379 const int n = (int)vtx.size();
380 const double minEpsAst = W.getPassport().wakeDiscretizationProperties.getMinEpsAst();
381
382 NeighborTree tree;
383 tree.Build(n, [&vtx](int i) -> const Point2D& { return vtx[i].r(); });
384
385 std::vector<double> gSorted(n), sigma(n);
386#pragma omp parallel for
387 for (int k = 0; k < n; ++k)
388 {
389 gSorted[k] = vtx[tree.idx[k]].g();
390 sigma[k] = vtx[tree.idx[k]].sigma();
391 Omega[k] = 0.0;
392 }
393
394#pragma omp parallel
395 {
396 std::vector<int> cand;
397 cand.reserve(256);
398
399#pragma omp for schedule(dynamic, 4)
400 for (int leaf = 0; leaf < tree.nLeaves; ++leaf)
401 {
403 double maxDiffRadius = 0.0;
404 for (int k = tree.LeafBegin(leaf); k < tree.LeafEnd(leaf); ++k)
405 maxDiffRadius = std::max(maxDiffRadius, 8.0 * std::max(sigma[k], minEpsAst));
406
407 tree.CollectLeaves(tree.LeafBox(leaf), maxDiffRadius, cand);
408
409 for (int k = tree.LeafBegin(leaf); k < tree.LeafEnd(leaf); ++k)
410 {
411 const int i = tree.idx[k];
412 const Point2D& posI = tree.pos[k];
413 const double domRad = std::max(sigma[k], minEpsAst);
414 const double iDomRad = 1.0 / domRad;
415 const double diffRadius2 = sqr(8.0 * domRad);
416
417 double selfOmega = 0.0;
418
419 for (int c : cand)
420 {
421 if (tree.LeafBox(c).Dist2(posI) >= diffRadius2)
422 continue;
423
424 for (int m = tree.LeafBegin(c); m < tree.LeafEnd(c); ++m)
425 {
426 const Point2D Rij = posI - tree.pos[m];
427 const double rij2 = Rij.length2();
428 if (rij2 < diffRadius2)
429 {
430 const double rij = sqrt(rij2);
431 const double expr = gSorted[m] * exp(-rij * iDomRad);
432 selfOmega += expr;
433 }
434 }
435 }
436
437 Omega[i] = selfOmega * IDPI * sqr(iDomRad);
438 }
439 }
440 }
441}//CalcOmegaFromWakeFast(...)
IDPI
const double IDPI
Число .
Definition defs.h:79
VMlib::Omega
Point2D Omega(const Point2D &a, const Point2D &b, const Point2D &c)
Вспомогательная функция вычисления величины .
Definition defs.cpp:274
Here is the call graph for this function:

◆ CalcVeloToWakeVP()

void Velocity::CalcVeloToWakeVP ( )
virtual

Вычисление скоростей в точках wakeVP.

Definition at line 1546 of file Velocity2D.cpp.

1547{
1548 //VMlib::vmTimer T1, T2, T3, T4, T5, T6;
1549
1550 //T1.start();
1551 std::vector<Point2D> velConvWake;
1552 std::vector<std::vector<Point2D>> velConvBou;
1553
1554 int addWSize = (int)W.getMeasureVP().getWakeVP().vtx.size();
1555
1556 velConvWake.resize(addWSize, { 0.0, 0.0 });
1557
1558 velConvBou.resize(W.getNumberOfBoundary());
1559 for (size_t i = 0; i < W.getNumberOfBoundary(); ++i)
1560 velConvBou[i].resize(addWSize, { 0.0, 0.0 });
1561
1562 std::vector<Point2D>& velocityRef = W.getNonConstMeasureVP().getNonConstVelocity();
1563 velocityRef.assign(addWSize, W.getV0());
1564
1565 //T1.stop();
1566
1567 //T2.start();
1568
1569#if (defined(USE_CUDA))
1570 W.getNonConstCuda().RefreshVP();
1571#endif
1572 //T2.stop();
1573
1574#if (defined(__CUDACC__) || defined(USE_CUDA)) && (defined(CU_VP))
1575 std::unique_ptr<BHcu::CudaTreeInfo>& cntrTree = W.getNonConstCuda().cntrTreeVP;
1576
1577 switch(W.getPassport().numericalSchemes.velocityComputation.second)
1578 {
1579 case 1:
1580
1581 if (W.getCurrentStep() == 0 || W.isAnyMovableOrDeformable())
1582 {
1583 cntrTree->MemoryAllocate((int)W.getCuda().n_CUDA_velVP);
1584 cntrTree->Update((int)W.getMeasureVP().getWakeVP().vtx.size(), W.getMeasureVP().getWakeVP().devVtxPtr);
1585 cntrTree->Build();
1586 }
1587 //T3.start();
1588 GPUCalcConvVeloToSetOfPointsFromWake(cntrTree, W.getMeasureVP().getWakeVP(), velocityRef, W.getNonConstMeasureVP().getNonConstDomainRadius(), true, false);
1589 //T3.stop();
1590 break;
1591 case 0:
1592 //T4.start();
1593 GPUCalcConvVeloToSetOfPointsFromWake(cntrTree, W.getMeasureVP().getWakeVP(), velocityRef, W.getNonConstMeasureVP().getNonConstDomainRadius(), true, false);
1594 //T4.stop();
1595 break;
1596 }
1597
1598
1599#else
1600 switch (W.getPassport().numericalSchemes.velocityComputation.second)
1601 {
1602 case 0:
1603 CalcConvVeloToSetOfPointsFromWake(W.getMeasureVP().getWakeVP(), velConvWake, W.getNonConstMeasureVP().getNonConstDomainRadius(), true, false);
1604 break;
1605 case 1:
1606 int order = W.getPassport().numericalSchemes.nbodyMultipoleOrder;
1607 double theta = W.getPassport().numericalSchemes.nbodyTheta;
1608
1609 W.getCntrTreeVP().Update(W.getMeasureVP().getWakeVP().vtx);
1610 W.getCntrTreeVP().Build();
1611 W.getCntrTreeVP().UpwardTraversal(W.getPassport().numericalSchemes.nbodyMultipoleOrder);
1612 W.getInflTreeWake().DownwardTraversalVorticesToPoints(W.getCntrTreeVP(), velConvWake, W.getNonConstMeasureVP().getNonConstDomainRadius(), theta, order, false, W.getPassport().wakeDiscretizationProperties.coefEpsAst);
1613 }
1614#endif
1615
1616
1617 for (size_t i = 0; i < W.getNumberOfBoundary(); ++i)
1618#if (defined(__CUDACC__) || defined(USE_CUDA)) && (defined(CU_VP))
1619 switch (W.getPassport().numericalSchemes.velocityComputation.second)
1620 {
1621 case 0:
1622 //T5.start();
1623 W.getNonConstBoundary(i).GPUCalcConvVelocityToSetOfPointsFromSheets(W.getMeasureVP().getWakeVP(), velConvBou[i]);
1624
1625 for (int s = 0; s < addWSize; ++s)
1626 velocityRef[s] += velConvWake[s];
1627
1628 for (size_t bou = 0; bou < velConvBou.size(); ++bou)
1629 for (int j = 0; j < addWSize; ++j)
1630 velocityRef[j] += velConvBou[bou][j];
1631 //T5.stop();
1632 break;
1633 case 1:
1634 if (i == 0)
1635 {
1636 //T6.start();
1637 GPUCalcConvVelocityToSetOfPointsFromSheets(cntrTree, W.getMeasureVP().getWakeVP(), velocityRef);
1638 //T6.stop();
1639 }
1640 break;
1641 }
1642#else
1643 W.getNonConstBoundary(i).CalcConvVelocityToSetOfPointsFromSheets(W.getMeasureVP().getWakeVP(), velConvBou[i]);
1644 for (int i = 0; i < addWSize; ++i)
1645 velocityRef[i] += velConvWake[i];
1646
1647 for (size_t bou = 0; bou < velConvBou.size(); ++bou)
1648 for (int j = 0; j < addWSize; ++j)
1649 velocityRef[j] += velConvBou[bou][j];
1650#endif
1651
1652 //std::cout << "calcV: " << T1.duration() << " " << T2.duration() << " " << T3.duration() << " "
1653 // << T4.duration() << " " << T5.duration() << " " << T6.duration() << "\n";
1654
1655
1656}// CalcVeloToWakeVP()
CU_VP
#define CU_VP
Definition Gpudefs.h:92
Here is the call graph for this function:
Here is the caller graph for this function:

◆ CPUGetFASTWakeInfluenceToRhs()

void Velocity::CPUGetFASTWakeInfluenceToRhs ( std::vector< double > & wakeRhs,
std::vector< double > & wakeRhsLin ) const

Definition at line 902 of file Velocity2D.cpp.

903{
904 const size_t nvt = W.getWake().vtx.size();
905
906
907 size_t nTotPan = 0;
908
909 for (size_t s = 0; s < W.getNumberOfAirfoil(); ++s)
910 nTotPan += W.getAirfoil(s).getNumberOfPanels();
911
912 wakeRhs.assign(nTotPan, 0.0);
913 wakeRhsLin.assign(nTotPan, 0.0);
914
915 if (nvt > 0)
916 {
917 auto& inflTree = W.getInflTreeWake();
918
919 auto& cntrTree = W.getCntrTreePnl();
920
921 inflTree.DownwardTraversalVorticesToPanels(
922 cntrTree,
923 wakeRhs,
924 wakeRhsLin,
925 W.getPassport().numericalSchemes.nbodyTheta,
926 W.getPassport().numericalSchemes.nbodyMultipoleOrder);
927 }
928}//CPUGetFASTWakeInfluenceToRhs(...)
Here is the caller graph for this function:

◆ FillRhs()

void Velocity::FillRhs ( Eigen::VectorXd & rhsReord) const

!!!!!!!!!!влияние присоединенных слоев от самого себя и от других профилей

Definition at line 1055 of file Velocity2D.cpp.

1056{
1057 Eigen::VectorXd locRhs;
1058 std::vector<double> lastRhs(W.getNumberOfBoundary());
1059
1060 size_t currentRow = 0;
1061 size_t currentSkosRow = 0;
1062
1063 size_t nAllVars = 0;
1064 for (size_t bou = 0; bou < W.getNumberOfBoundary(); ++bou)
1065 nAllVars += W.getBoundary(bou).GetUnknownsSize();
1066
1067 //Временные массивы для хранения правой части по быстрому методу (размеры задаются внутри CPUGetFASTWakeInfluenceToRhs)
1068 std::vector<double> fastWakeRhs;
1069 std::vector<double> fastWakeRhsLin;
1070
1071#ifndef USE_CUDA
1072 if (W.getPassport().numericalSchemes.velocityComputation.second == 1) //Быстрый метод на CPU
1073 CPUGetFASTWakeInfluenceToRhs(fastWakeRhs, fastWakeRhsLin);
1074#endif
1075
1076 size_t curGlobPnl = 0;
1077
1078 W.timerRhs.reset();
1079
1080 for (size_t bou = 0; bou < W.getNumberOfBoundary(); ++bou)
1081 {
1082 const Airfoil& afl = W.getAirfoil(bou);
1083 size_t np = afl.getNumberOfPanels();
1084
1085 size_t nVars;
1086
1087 nVars = W.getBoundary(bou).GetUnknownsSize();
1088 locRhs.resize(nVars);
1089
1090 std::vector<double> wakeRhs;
1091
1092 double tt1 = omp_get_wtime();
1093
1094#if (defined(__CUDACC__) || defined(USE_CUDA)) && (defined(CU_RHS))
1095 W.timerRhs.start();
1096 if (W.getPassport().numericalSchemes.velocityComputation.second == 1)
1097 GPUFASTGetWakeInfluenceToRhs(afl, wakeRhs);
1098 else
1099 GPUGetWakeInfluenceToRhs(afl, wakeRhs);
1100 W.timerRhs.stop();
1101#else
1102 if (W.getPassport().numericalSchemes.velocityComputation.second == 1) //Собираем результат из быстрого метода
1103 {
1104 size_t shDim = W.getBoundary(afl.numberInPassport).sheetDim;
1105
1106 wakeRhs.assign(np * shDim, 0.0);
1107
1108 for (size_t i = 0; i < np; ++i)
1109 {
1110 wakeRhs[i] = fastWakeRhs[curGlobPnl + i];
1111
1112 if (shDim != 1)
1113 wakeRhs[np + i] = fastWakeRhsLin[curGlobPnl + i];
1114 }
1115 curGlobPnl += np;
1116 }
1117 else
1118 GetWakeInfluenceToRhs(afl, wakeRhs); //Прямой расчет
1119
1120 //std::string vname = VMlib::fileNameStep("rhsLinBH", 0, wakeRhs.size(), "txt");
1121 //std::ofstream velofile;
1122 //velofile.open(W.getPassport().dir + vname);
1123
1124 //for (int i = 0; i < wakeRhs.size(); ++i)
1125 // velofile << wakeRhs[i] << "\n";
1126 //velofile.close();
1127 //exit(-124);
1128#endif
1129
1130 std::vector<double> vInfRhs;
1131 afl.GetInfluenceFromVInfToPanel(vInfRhs);
1132
1133#pragma omp parallel for \
1134 default(none) \
1135 shared(locRhs, afl, bou, wakeRhs, vInfRhs, np) \
1136 schedule(dynamic, DYN_SCHEDULE)
1137 for (int i = 0; i < (int)afl.getNumberOfPanels(); ++i)
1138 {
1139 locRhs(i) = -vInfRhs[i] - wakeRhs[i] + 0.25 * ((afl.getV(i) + afl.getV(i + 1)) & afl.tau[i]); //0.25 * (afl.getV(i) + afl.getV(i + 1))*afl.tau[i] - прямолинейные
1140 if (W.getBoundary(bou).sheetDim > 1)
1141 locRhs(np + i) = -vInfRhs[np + i] - wakeRhs[np + i];
1142
1144 //Если включен быстрый метод и есть подвижные тела, то пока что все равно вычисляем IQ
1145 if (W.isAnyMovableOrDeformable())
1146 {
1147 for (size_t q = 0; q < W.getNumberOfBoundary(); ++q)
1148 {
1149 const auto& sht = W.getBoundary(q).sheets;
1150 const auto& iq = W.getIQ(bou, q);
1151
1152 const Airfoil& aflOther = W.getAirfoil(q);
1153 if (W.getMechanics(q).isMoves)
1154 {
1155 for (size_t j = 0; j < aflOther.getNumberOfPanels(); ++j)
1156 {
1157 if ((i != j) || (bou != q))
1158 {
1159 locRhs(i) += -iq.first(i, j) * sht.attachedVortexSheet(j, 0);
1160 locRhs(i) += -iq.second(i, j) * sht.attachedSourceSheet(j, 0);
1161 }//if (i != j)
1162
1163 }//for j
1164 }
1165 }//for q
1166 }
1167 }//for i
1168
1169 lastRhs[bou] = 0.0;
1170
1171#pragma omp parallel for
1172 for (int q = 0; q < (int)afl.gammaThrough.size(); ++q)
1173 lastRhs[bou] += afl.gammaThrough[q];
1174
1176
1177 const double currT = W.getCurrentTime();
1178 const double dt = W.getPassport().timeDiscretizationProperties.dt;
1179
1180 lastRhs[bou] += (W.getMechanics(bou).circulation - W.getMechanics(bou).circulationOld) * (W.getAirfoil(bou).inverse ? 1.0 : -1.0);
1181
1183
1184 //размазываем правую часть
1185 for (size_t i = 0; i < nVars; ++i)
1186 rhsReord(i + currentSkosRow) = locRhs(i);
1187
1188 rhsReord(nAllVars + bou) = lastRhs[bou];
1189
1190 currentRow += nVars + 1;
1191 currentSkosRow += nVars;
1192 }// for bou
1193
1194
1196 /*
1197 std::ofstream treeTimeFile;
1198 if (W.getCurrentStep() == 0)
1199 {
1200 treeTimeFile.open(W.getPassport().dir + "/dbg/rhsTime.csv");
1201 treeTimeFile << "step,time,Nvtx,Npnl,tDNW,tDNWgpu\n";
1202 }
1203 else
1204 treeTimeFile.open(W.getPassport().dir + "/dbg/rhsTime.csv", std::ios::app);
1205
1206
1207 treeTimeFile << W.getCurrentStep() << ',' << W.getCurrentTime() << ',' << W.getInflTreeWake().object.size() << ',' << W.getCntrTreePnl().object.size() << ',' \
1208 << timerRhsTest.duration() << ',' << W.timerRhs.duration() << '\n';
1209
1210 treeTimeFile.close();
1211 //*/
1212}
VM2D::AirfoilGeometry::getV
const Point2D & getV(size_t q) const
Возврат константной ссылки на скорость вершины профиля
Definition Airfoil2D.h:137
VM2D::Airfoil::gammaThrough
std::vector< double > gammaThrough
Суммарные циркуляции вихрей, пересекших панели профиля на прошлом шаге
Definition Airfoil2D.h:277
VM2D::Airfoil::GetInfluenceFromVInfToPanel
virtual void GetInfluenceFromVInfToPanel(std::vector< double > &vInfRhs) const
Вычисление влияния набегающего потока на панель для правой части
Definition Airfoil2D.cpp:1557
VM2D::Airfoil::numberInPassport
const size_t numberInPassport
Номер профиля в паспорте
Definition Airfoil2D.h:188
VM2D::Velocity::GetWakeInfluenceToRhs
void GetWakeInfluenceToRhs(const Airfoil &afl, std::vector< double > &wakeRhs) const
Генерация вектора влияния вихревого следа на профиль
Definition Velocity2D.cpp:865
VM2D::Velocity::CPUGetFASTWakeInfluenceToRhs
void CPUGetFASTWakeInfluenceToRhs(std::vector< double > &wakeRhs, std::vector< double > &wakeRhsLin) const
Definition Velocity2D.cpp:902
Here is the call graph for this function:

◆ GetWakeInfluenceToRhs()

void Velocity::GetWakeInfluenceToRhs ( const Airfoil & afl,
std::vector< double > & wakeRhs ) const

Генерация вектора влияния вихревого следа на профиль

Генерирует вектор влияния вихревого следа на профиль, используемый затем для расчета вектора правой части.

Parameters
[in]aflконстантная ссылка на профиль, правая часть для которого вычисляется
[out]wakeRhsссылка на вектор влияния вихревого следа на ОДИН профиль

Definition at line 865 of file Velocity2D.cpp.

866{
867 size_t np = afl.getNumberOfPanels();
868 size_t shDim = W.getBoundary(afl.numberInPassport).sheetDim;
869
870 wakeRhs.resize(W.getBoundary(afl.numberInPassport).GetUnknownsSize());
871
872 //локальные переменные для цикла
873 std::vector<double> velI(shDim, 0.0);
874
875#pragma omp parallel for default(none) shared(shDim, afl, np, wakeRhs, IDPI) private(velI)
876 for (int i = 0; i < np; ++i)
877 {
878 velI.assign(shDim, 0.0);
879
880 if (W.getWake().vtx.size() > 0)
881 {
882 //Учет влияния следа
883 afl.GetInfluenceFromVorticesToPanel(i, W.getWake().vtx.data(), W.getWake().vtx.size(), velI);
884 }
885
886 if (W.getSource().vtx.size() > 0)
887 {
888 //Учет влияния источников
889 afl.GetInfluenceFromSourcesToPanel(i, W.getSource().vtx.data(), W.getSource().vtx.size(), velI);
890 }
891
892 for (size_t j = 0; j < shDim; ++j)
893 velI[j] *= IDPI / afl.len[i];
894
895 wakeRhs[i] = velI[0];
896
897 if (shDim != 1)
898 wakeRhs[np + i] = velI[1];
899 }//for i
900}//GetWakeInfluenceToRhs(...)
VM2D::Airfoil::GetInfluenceFromSourcesToPanel
virtual void GetInfluenceFromSourcesToPanel(size_t panel, const Vortex2D *ptr, ptrdiff_t count, std::vector< double > &panelRhs) const
Вычисление влияния части подряд идущих источников из области течения на панель для правой части
Definition Airfoil2D.cpp:1540
VM2D::Airfoil::GetInfluenceFromVorticesToPanel
virtual void GetInfluenceFromVorticesToPanel(size_t panel, const Vortex2D *ptr, ptrdiff_t count, std::vector< double > &panelRhs) const
Вычисление влияния части подряд идущих вихрей из вихревого следа на панель для правой части
Definition Airfoil2D.cpp:1533
Here is the call graph for this function:
Here is the caller graph for this function:

◆ LimitDiffVelo()

void Velocity::LimitDiffVelo ( std::vector< Point2D > & diffVel)

Контроль больших значений диффузионных скоростей

Parameters
[in,out]diffVelссылка на вектор диффузионных скоростей

Definition at line 173 of file Velocity2D.cpp.

174{
175 for (size_t i = 0; i < diffVel.size(); ++i)
176 {
177 Point2D& diffV = diffVel[i];
178
179 if (diffV.length() > 1.5 * W.getPassport().physicalProperties.vRef)
180 diffV.normalize(1.5 * W.getPassport().physicalProperties.vRef);
181 }
182}
numvector::normalize
void normalize(P newlen=1.0)
Нормирование вектора на заданную длину
Definition numvector.h:419
Here is the call graph for this function:
Here is the caller graph for this function:

◆ ResizeAndZero()

void Velocity::ResizeAndZero ( )

Очистка старых массивов под хранение скоростей, выделение новой памяти и обнуление

Вызывается в CalcVortexVelo() на каждом шаге расчета перед непосредственно расчетом скоростей

Definition at line 767 of file Velocity2D.cpp.

768{
769 W.getTimers().start("ConvVel");
770 wakeVortexesParams.convVelo.clear();
771 wakeVortexesParams.convVelo.resize(W.getWake().vtx.size(), { 0.0, 0.0 });
772 W.getTimers().stop("ConvVel");
773
774 W.getTimers().start("DiffVel");
775 wakeVortexesParams.I0.clear();
776 wakeVortexesParams.I0.resize(W.getWake().vtx.size(), 0.0);
777
778 wakeVortexesParams.I1.clear();
779 wakeVortexesParams.I1.resize(W.getWake().vtx.size(), 0.0);
780
781 wakeVortexesParams.I2.clear();
782 wakeVortexesParams.I2.resize(W.getWake().vtx.size(), { 0.0, 0.0 });
783
784 wakeVortexesParams.I3.clear();
785 wakeVortexesParams.I3.resize(W.getWake().vtx.size(), { 0.0, 0.0 });
786
787 wakeVortexesParams.diffVelo.clear();
788 wakeVortexesParams.diffVelo.resize(W.getWake().vtx.size(), { 0.0, 0.0 });
789
790 wakeVortexesParams.epsastWake.clear();
791 wakeVortexesParams.epsastWake.resize(W.getWake().vtx.size(), 0.0);
792 W.getTimers().stop("DiffVel");
793
794 //Создаем массивы под виртуальные вихри
795 W.getTimers().start("ConvVel");
796 virtualVortexesParams.clear();
797 virtualVortexesParams.resize(W.getNumberOfBoundary());
798 W.getTimers().stop("ConvVel");
799
800 for (size_t bou = 0; bou < W.getNumberOfBoundary(); ++bou)
801 {
802 W.getTimers().start("ConvVel");
803 virtualVortexesParams[bou].convVelo.clear();
804 virtualVortexesParams[bou].convVelo.resize(W.getBoundary(bou).virtualWake.vtx.size(), { 0.0, 0.0 });
805 W.getTimers().stop("ConvVel");
806
807 W.getTimers().start("DiffVel");
808 virtualVortexesParams[bou].I0.clear();
809 virtualVortexesParams[bou].I0.resize(W.getBoundary(bou).virtualWake.vtx.size(), 0.0);
810
811 virtualVortexesParams[bou].I1.clear();
812 virtualVortexesParams[bou].I1.resize(W.getBoundary(bou).virtualWake.vtx.size(), 0.0);
813
814 virtualVortexesParams[bou].I2.clear();
815 virtualVortexesParams[bou].I2.resize(W.getBoundary(bou).virtualWake.vtx.size(), { 0.0, 0.0 });
816
817 virtualVortexesParams[bou].I3.clear();
818 virtualVortexesParams[bou].I3.resize(W.getBoundary(bou).virtualWake.vtx.size(), { 0.0, 0.0 });
819
820 virtualVortexesParams[bou].diffVelo.clear();
821 virtualVortexesParams[bou].diffVelo.resize(W.getBoundary(bou).virtualWake.vtx.size(), { 0.0, 0.0 });
822
823 virtualVortexesParams[bou].epsastWake.clear();
824 virtualVortexesParams[bou].epsastWake.resize(W.getBoundary(bou).virtualWake.vtx.size(), 0.0);
825
826 W.getTimers().stop("DiffVel");
827 }
828}//ResizeAndZero()

◆ SaveVisStress()

void Velocity::SaveVisStress ( )

Сохранение вязких напряжений

Вызывается в CalcDiffVelo()

Definition at line 831 of file Velocity2D.cpp.

832{
833 if ((W.getPassport().timeDiscretizationProperties.saveVtxStep > 0) && (W.ifDivisible(W.getPassport().timeDiscretizationProperties.saveVisStress)))
834 // if ((W.getPassport().timeDiscretizationProperties.saveVTK > 0) && (W.ifDivisible(10)) && (W.getNumberOfAirfoil() > 0))
835 {
836 for (size_t q = 0; q < W.getNumberOfAirfoil(); ++q)
837 {
838 if (q == 0)
839 VMlib::CreateUserDirectory(W.getPassport().dir, "visStress");
840
841 std::stringstream ss;
842 ss << "VisStress_" << q << "-";
843 std::string fname = VMlib::fileNameStep(ss.str(), W.getPassport().timeDiscretizationProperties.nameLength, W.getCurrentStep(), "txt");
844 std::ofstream outfile;
845 outfile.open(W.getPassport().dir + "visStress/" + fname);
846
847 outfile << W.getAirfoil(q).viscousStress.size() << std::endl; //Сохранение числа вихрей в пелене
848
849 for (size_t i = 0; i < W.getAirfoil(q).viscousStress.size(); ++i)
850 {
851 const Point2D& r = 0.5 * (W.getAirfoil(q).getR(i + 1) + W.getAirfoil(q).getR(i));
852
853 //Делим обратно на длину панели, чтобы получить размерность напряжения
854 double gi = W.getAirfoil(q).viscousStress[i] / W.getAirfoil(q).len[i];
855 outfile << static_cast<int>(i) << " " << r[0] << " " << r[1] << " " << gi << std::endl;
856 }//for i
857 outfile.close();
858 }
859 }
860}//SaveVisStress()
VMlib::CreateUserDirectory
void CreateUserDirectory(const std::string &dir, const std::string &name)
Создание каталога
Definition defs.h:447
VMlib::fileNameStep
std::string fileNameStep(const std::string &name, int length, size_t number, const std::string &ext)
Формирование имени файла
Definition defs.h:387
Here is the call graph for this function:
Here is the caller graph for this function:

Member Data Documentation

◆ optimizedVelocity

OptimizedVelocity VM2D::Velocity::optimizedVelocity

Definition at line 117 of file Velocity2D.h.

◆ virtualVortexesParams

std::vector<VortexesParams> VM2D::Velocity::virtualVortexesParams

Вектор струтур, определяющий параметры виртуальных вихрей для профилей

Definition at line 115 of file Velocity2D.h.

◆ W

const World2D& VM2D::Velocity::W
protected

Константная ссылка на решаемую задачу

Definition at line 108 of file Velocity2D.h.

◆ wakeVortexesParams

VortexesParams VM2D::Velocity::wakeVortexesParams

Струтура, определяющая параметры вихрей в следе

Definition at line 112 of file Velocity2D.h.

The documentation for this class was generated from the following files: