VM2D 1.14
Vortex methods for 2D flows simulation
Loading...
Searching...
No Matches
World2D.cpp
Go to the documentation of this file.
1/*--------------------------------*- VM2D -*-----------------*---------------*\
2| ## ## ## ## #### ##### | | Version 1.14 |
3| ## ## ### ### ## ## ## ## | VM2D: Vortex Method | 2026/03/06 |
4| ## ## ## # ## ## ## ## | for 2D Flow Simulation *----------------*
5| #### ## ## ## ## ## | Open Source Code |
6| ## ## ## ###### ##### | https://www.github.com/vortexmethods/VM2D |
7| |
8| Copyright (C) 2017-2026 I. Marchevsky, K. Sokol, E. Ryatina, A. Kolganova |
9*-----------------------------------------------------------------------------*
10| File name: World2D.cpp |
11| Info: Source code of VM2D |
12| |
13| This file is part of VM2D. |
14| VM2D is free software: you can redistribute it and/or modify it |
15| under the terms of the GNU General Public License as published by |
16| the Free Software Foundation, either version 3 of the License, or |
17| (at your option) any later version. |
18| |
19| VM2D is distributed in the hope that it will be useful, but WITHOUT |
20| ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or |
21| FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License |
22| for more details. |
23| |
24| You should have received a copy of the GNU General Public License |
25| along with VM2D. If not, see <http://www.gnu.org/licenses/>. |
26\*---------------------------------------------------------------------------*/
27
28
40#include "World2D.h"
41
42#include "Airfoil2DRigid.h"
43#include "Airfoil2DDeformable.h"
44
45#include "Boundary2DConstLayerAver.h"
46#include "Boundary2DLinLayerAver.h"
47#include "Boundary2DVortexCollocN.h"
48
49#include "MeasureVP2D.h"
50
51#include "Mechanics2DRigidImmovable.h"
52#include "Mechanics2DRigidGivenLaw.h"
53#include "Mechanics2DRigidOscillPart.h"
54#include "Mechanics2DRigidRotatePart.h"
55#include "Mechanics2DDeformable.h"
56
57#include "StreamParser.h"
58
59#include "Velocity2DBiotSavart.h"
60#include "Velocity2DBarnesHut.h"
61
62#include "Wake2D.h"
63
64#include "Gmres2D.h"
65#include "treeKernels.cuh"
66#include <type_traits>
67#include <numeric>
68#include "nummatrix.h"
69
70using namespace VM2D;
71
72//Конструктор
73World2D::World2D(const VMlib::PassportGen& passport_) :
74 WorldGen(passport_),
75 passport(dynamic_cast<const Passport&>(passport_)),
76 cuda(Gpu(*this))
77{
78 std::stringstream ss;
79 ss << "#" << passport.problemNumber << " (" << passport.problemName << ")";
80 info.assignStream(defaults::defaultWorld2DLogStream, ss.str());
81
82 currentTime = passport.timeDiscretizationProperties.timeStart;
83 currentStep = 0;
84
85 std::vector<std::string> timerLabels = { "Step", "MatRhs", "Solve", "ConvVel", "Nut", /*"Nbody",*/ "DiffVel", "Force", "VelPres", "Inside", "Restr", "Save"};
86 timers = std::make_unique<VMlib::TimersGen>(*this, timerLabels);
87
88 wake.reset(new Wake(*this));
89 // загрузка пелены из файла
90 if (passport.wakeDiscretizationProperties.fileWake != "")
91 {
92 if (VMlib::fileExistTest(passport.wakesDir + passport.wakeDiscretizationProperties.fileWake, getInfo(), false, { "txt", "TXT" }))
93 wake->ReadFromFile(passport.wakesDir, passport.wakeDiscretizationProperties.fileWake); //Считываем из каталога с пеленой
94 else
95 {
96 if (VMlib::fileExistTest(passport.dir + passport.wakeDiscretizationProperties.fileWake, getInfo(), true, { "txt", "TXT" }))
97 wake->ReadFromFile(passport.dir, passport.wakeDiscretizationProperties.fileWake); //Считываем из каталога с пеленой
98 }
99
100
101 }
102
103 source.reset(new WakeDataBase(*this));
104 // загрузка положений источников из файла
105 if (passport.wakeDiscretizationProperties.fileSource != "")
106 source->ReadFromFile(passport.dir, passport.wakeDiscretizationProperties.fileSource); //Считываем из текущего каталога
107
108
109 switch (passport.numericalSchemes.velocityComputation.second)
110 {
111 case 0:
112 velocity.reset(new VelocityBiotSavart(*this));
113 break;
114 case 1:
115 velocity.reset(new VelocityBarnesHut(*this));
116 inflTreeWake.reset(new CpuTreeInfo(tree_T::vortex, object_T::point4, scheme_T::noScheme));
117 cntrTreeWake.reset(new CpuTreeInfo(tree_T::contr, object_T::point4, scheme_T::noScheme));
118 cntrTreeVP.reset(new CpuTreeInfo(tree_T::contr, object_T::point2, scheme_T::noScheme));
119 if(passport.numericalSchemes.boundaryCondition.second == 0)
120 cntrTreePnl.reset(new CpuTreeInfo(tree_T::contr, object_T::panel, scheme_T::constScheme));
121 else
122 cntrTreePnl.reset(new CpuTreeInfo(tree_T::contr, object_T::panel, scheme_T::linScheme));
123 break;
124 }
125
126 velocity->virtualVortexesParams.resize(passport.airfoilParams.size());
127
128 auto CreateBoundary = [this](size_t i) {
129 switch (passport.numericalSchemes.boundaryCondition.second)
130 {
131 case 0:
132 this->boundary.emplace_back(new BoundaryVortexCollocN(*this, i));
133 //info('e') << "BoundaryMDV is not implemented now! " << std::endl;
134 //exit(1);
135 break;
136
137 case 1:
138 this->boundary.emplace_back(new BoundaryConstLayerAver(*this, i));
139 break;
140
141 case 2:
142 this->boundary.emplace_back(new BoundaryLinLayerAver(*this, i));
143 break;
144
145 default:
146 info('e') << "Unknown scheme!" << std::endl;
147 exit(1);
148 }
149 };
150
151
152 for (size_t i = 0; i < passport.airfoilParams.size(); ++i)
153 {
154
155 switch (passport.airfoilParams[i].mechanicalSystemType)
156 {
157 case 0:
158 airfoil.emplace_back(new AirfoilRigid(*this, i));
159 airfoil[i]->ReadFromFile(passport.airfoilsDir);
160 CreateBoundary(i);
161 mechanics.emplace_back(new MechanicsRigidImmovable(*this, i));
162 break;
163
164 case 1:
165 airfoil.emplace_back(new AirfoilRigid(*this, i));
166 airfoil[i]->ReadFromFile(passport.airfoilsDir);
167 CreateBoundary(i);
168 mechanics.emplace_back(new MechanicsRigidGivenLaw(*this, i));
169 break;
170
171 case 2:
172 airfoil.emplace_back(new AirfoilRigid(*this, i));
173 airfoil[i]->ReadFromFile(passport.airfoilsDir);
174 CreateBoundary(i);
175 mechanics.emplace_back(new MechanicsRigidOscillPart(*this, i));
176 break;
177
178 case 3:
179 airfoil.emplace_back(new AirfoilRigid(*this, i));
180 airfoil[i]->ReadFromFile(passport.airfoilsDir);
181 CreateBoundary(i);
182 mechanics.emplace_back(new MechanicsRigidRotatePart(*this, i));
183 break;
184
185 case 4:
186 airfoil.emplace_back(new AirfoilDeformable(*this, i));
187 airfoil[i]->ReadFromFile(passport.airfoilsDir);
188 CreateBoundary(i);
189 mechanics.emplace_back(new MechanicsDeformable(*this, i));
190 break;
191 }
192 }
193
194 if (getPassport().wakeDiscretizationProperties.sigma0 == 0)
195 {
196 double sumLength = 0.0;
197 size_t totNumPan = 0;
198 for (size_t bou = 0; bou < getNumberOfAirfoil(); ++bou)
199 {
200 for (size_t pnl = 0; pnl < getAirfoil(bou).getNumberOfPanels(); ++pnl)
201 sumLength += getAirfoil(bou).len[pnl];
202 totNumPan += getAirfoil(bou).getNumberOfPanels();
203 }
204 double eps = 0.5 * sumLength / totNumPan;
205 getNonConstPassport().wakeDiscretizationProperties.sigma0 = eps;
206 info('i') << "sigma0 = " << getPassport().wakeDiscretizationProperties.sigma0 << " is calculated automatically" << std::endl;
207 }
208
209 if (getPassport().wakeDiscretizationProperties.epscol == 0)
210 {
211 getNonConstPassport().wakeDiscretizationProperties.epscol = (2.0 / 3.0) * getPassport().wakeDiscretizationProperties.sigma0;
212
213 info('i') << "epscol = " << getPassport().wakeDiscretizationProperties.epscol << " is calculated automatically" << std::endl;
214 }
215
216 for (size_t afl = 0; afl < passport.airfoilParams.size(); ++afl)
217 if (passport.airfoilParams[afl].chord == 0)
218 {
219 auto& prm = getNonConstPassport().airfoilParams[afl];
220 prm.chord = (prm.initialGab.second[0] - prm.initialGab.first[0]) * prm.scale[0];
221 info('i') << "airfoil #" << afl << " chord = " << prm.chord << " is calculated automatically" << std::endl;
222 }
223
224
225#ifdef USE_CUDA
226 if (getNumberOfAirfoil() > 0)
227 if (passport.numericalSchemes.linearSystemSolver.second == 1 || passport.numericalSchemes.linearSystemSolver.second == 2)
228 cuda.Gmres.reset(new GmresSolver(*this));
229#endif
230
231 //2026-03-28
232 //считываем массив точек для подсчета и вывода поля скоростей и давлений
233 measureVP.reset(new MeasureVP(*this));
234
235 if (passport.timeDiscretizationProperties.saveVPstep != 0)
236 measureVP->ReadPointsFromFile(passport.dir);
237
238#ifdef OPTIMIZER
239// std::vector<Point2D> VPPoints, VPHistory;
240// for (size_t a = 0; a < getNumberOfAirfoil(); ++a)
241// {
242// const auto& afl = *airfoil[a];
243// for (size_t p = 0; p < afl.getNumberOfPanels(); ++p)
244// VPHistory.push_back(0.5 * (afl.getR(p) + afl.getR(p + 1)) + afl.nrm[p] * passport.airfoilParams[a].chord * 0.01);
245// }
246// if (passport.timeDiscretizationProperties.saveVPstep != 0)
247// measureVP->SetPoints(VPPoints, VPHistory);
248#endif
249
250 IQ.resize(passport.airfoilParams.size());
251 for (size_t i = 0; i < passport.airfoilParams.size(); ++i)
252 IQ[i].resize(passport.airfoilParams.size());
253
254
255 info.endl();
256 info.endl();
257 info('i') << "Start solving problem " << passport.problemName << std::endl;
258 info.endl();
259
260 VMlib::PrintLogoToStream(info('_') << std::endl);
261}//World2D(...)
262
263
264
266//void World2D::ZeroStep()
267//{
268// getTimestat().ToZero();
269// CalcPanelsVeloAndAttachedSheets();
270//
271// getNonConstMeasureVP().Initialization();
272//
273// CalcAndSolveLinearSystem();
274//}//ZeroStep()
275
276
277
278//Основная функция выполнения одного шага по времени
279void World2D::Step() // ЮИ
280{
281 //try {
282 //Очистка статистики
283 getTimers().resetAll();
284
285#ifdef USE_CUDA
286 cuSetCurrentStep((int)currentStep, 1);
287#endif // USE_CUDA
288
289
290 //Засечка времени в начале шага
291 getTimers().start("Step");
292
293 //Обновление параметров паспорта
294 if (passport.timeDiscretizationProperties.revisePassportStep != 0 && currentStep && (currentStep % passport.timeDiscretizationProperties.revisePassportStep == 0))
295 {
296 std::unique_ptr<Passport> revisePspPtr;
297 revisePspPtr.reset(new VM2D::Passport(info, passport.dir, passport.problemNumber,
298 passport.fileFullName, passport.mechanicsFileFullName, passport.defaultsFileFullName, passport.switchersFileFullName, passport.varLine,
299 passport.timeDiscretizationProperties.reviseParameters));
300
301 this->getNonConstPassport().GetReviseParamsFromParser(*revisePspPtr, passport.timeDiscretizationProperties.reviseParameters);
302
303 }
304
305
307
308 size_t countStrongCoupling = 0;
309 for (size_t m = 0; m < mechanics.size(); ++m)
310 {
311 MechanicsRigidOscillPart* mechVar = dynamic_cast<MechanicsRigidOscillPart*>(mechanics[m].get());
312 if (mechVar && getPassport().airfoilParams[m].addedMass.length2() > 0)
313 {
314 mechVar->getStrongCoupling() = true;
315 ++countStrongCoupling;
316 }
317 }
318
319 bool semiImplicitStrategy = ((countStrongCoupling == mechanics.size()) && (mechanics.size() > 0));
320 if ((currentStep == 0) && (semiImplicitStrategy))
321 info('i') << "Strong (semi-implicit) coupling strategy" << std::endl;
322
323//НЕПОДВИЖНЫЕ ТЕЛА
324//*
325 int nTotPan = 0;
326 for (size_t s = 0; s < getNumberOfAirfoil(); ++s)
327 nTotPan += (int)getAirfoil(s).getNumberOfPanels();
328
329
330 if (!semiImplicitStrategy)
331 {
332 CalcPanelsVeloAndAttachedSheets();
333
334 timerInitialBuild.reset();
335 timerInitialBuild.start();
336 if (getPassport().numericalSchemes.velocityComputation.second == 0 && getPassport().numericalSchemes.linearSystemSolver.second == 2)
337 {
338#ifdef USE_CUDA
339 cuda.RefreshAfls(3);
340 if (nTotPan > 0)
341 {
342 auto& afl = getAirfoil(0);
343 auto& treePnlVrt = *getCuda().inflTreePnlVortex;
344 if (getCurrentStep() == 0)
345 treePnlVrt.MemoryAllocate((int)getCuda().n_CUDA_pnls);
346
347 if (getCurrentStep() == 0 || isAnyMovableOrDeformable())
348 {
349 //Построение дерева для влияющих панелей (вихревых слоев)
350 treePnlVrt.UpdatePanelGeometry(nTotPan, (double4*)afl.devRPtr);
351 treePnlVrt.Build();
352 }
353
354 treePnlVrt.UpdatePanelAttachedVortexIntensity(afl.devAttachedVortexSheetPtr, afl.devAttachedVortexSheetLinPtr);
355 treePnlVrt.UpwardTraversal(getPassport().numericalSchemes.gmresMultipoleOrder);
356 }
357#endif
358 }
359
360 if (getPassport().numericalSchemes.velocityComputation.second == 1)
361 {
362#ifdef USE_CUDA
363 cuda.RefreshWake(3);
364 cuda.RefreshAfls(3);
365
366 //Построение дерева для вихрей
367 auto& treeWake = *getCuda().inflTreeWake;
368 treeWake.MemoryAllocate((int)getCuda().n_CUDA_wake);
369 treeWake.Update((int)getWake().vtx.size(), getWake().devVtxPtr);
370 treeWake.Build();
371 treeWake.UpwardTraversal(getPassport().numericalSchemes.nbodyMultipoleOrder);
372
373 if (nTotPan > 0)
374 {
375 auto& afl = getAirfoil(0);
376 auto& treePnl = *getCuda().cntrTreePnl;
377 auto& treePnlVrt = *getCuda().inflTreePnlVortex;
378 auto& treePnlSrc = *getCuda().inflTreePnlSource;
379 auto& treePnlAux = *getCuda().auxTreePnl;
380
381 if (getCurrentStep() == 0)
382 {
383 treePnl.MemoryAllocate((int)getCuda().n_CUDA_pnls);
384 treePnlAux.MemoryAllocate((int)getCuda().n_CUDA_pnls);
385 treePnlVrt.MemoryAllocate((int)getCuda().n_CUDA_pnls);
386 if (isAnyMovableOrDeformable())
387 treePnlSrc.MemoryAllocate((int)getCuda().n_CUDA_pnls);
388 }
389
390 if (getCurrentStep() == 0 || isAnyMovableOrDeformable())
391 {
392 //Построение дерева для влияющих панелей (вихревых слоев)
393 treePnlVrt.UpdatePanelGeometry(nTotPan, (double4*)afl.devRPtr);
394 treePnlVrt.Build();
395
396 //Построение контрольного дерева панелей
397 treePnl.UpdatePanelGeometry((int)nTotPan, (double4*)afl.devRPtr);
398 treePnl.Build();
399
400 //Построение дерева для влияющих панелей (источников)
401 if (isAnyMovableOrDeformable())
402 {
403 treePnlSrc.UpdatePanelGeometry(nTotPan, (double4*)afl.devRPtr);
404 treePnlSrc.UpdatePanelAttachedSourceIntensity(afl.devAttachedSourceSheetPtr, afl.devAttachedSourceSheetLinPtr);
405 treePnlSrc.Build();
406 treePnlSrc.UpwardTraversal(getPassport().numericalSchemes.nbodyMultipoleOrder);
407 }
408 }
409
410 treePnlVrt.UpdatePanelAttachedVortexIntensity(afl.devAttachedVortexSheetPtr, afl.devAttachedVortexSheetLinPtr);
411 treePnlVrt.UpwardTraversal(getPassport().numericalSchemes.nbodyMultipoleOrder);
412
413
414 //временное дерево для отладки процедуры сохранения
415 /*
416 if (getWake().vtx.size() > 2)
417 {
418 CpuTreeInfo testTree(tree_T::vortex, object_T::point4, scheme_T::noScheme);
419 testTree.Update(getWake().vtx);
420 testTree.Build();
421 testTree.UpwardTraversal(getPassport().numericalSchemes.nbodyMultipoleOrder);
422 testTree.SaveTreeTopology(getPassport().dir + "testTree-" + std::to_string(currentStep) + ".tree");
423 }
424 */
425 }
426#else
427 if (getWake().vtx.size() > 0)
428 {
429 inflTreeWake->Update(getWake().vtx);
430 inflTreeWake->Build();
431 inflTreeWake->UpwardTraversal(getPassport().numericalSchemes.nbodyMultipoleOrder);
432 }
433 std::vector<std::pair<Point2D, Point2D>> panels;
434
435 for (size_t bou = 0; bou < getNumberOfAirfoil(); ++bou)
436 {
437 const auto& afl = getAirfoil(bou);
438
439 for (size_t i = 0; i < afl.getNumberOfPanels(); ++i)
440 panels.push_back({ afl.getR(i), afl.getR(i + 1) });
441 }
442 if (panels.size() > 0)
443 {
444 cntrTreePnl->UpdatePanelGeometry(panels, std::max(4, (int)(log2(panels.size())) - 2));
445 cntrTreePnl->Build();
446 cntrTreePnl->UpwardTraversal(getPassport().numericalSchemes.nbodyMultipoleOrder);
447 }
448#endif
449 }
450 timerInitialBuild.stop();
451
452 measureVP->Initialization();
453 CalcAndSolveLinearSystem();
454
455 //added masses
456 if (getPassport().physicalProperties.typeAccel.second == 3)
457 {
458 std::vector<Point2D> lambdaAdd(getNumberOfAirfoil(), {0.0, 0.0});
459 std::vector<double> muAdd(getNumberOfAirfoil());
460 for (size_t bou = 0; bou < getNumberOfAirfoil(); ++bou)
461 {
462 if (dynamic_cast<MechanicsRigidGivenLaw*>(&getNonConstMechanics(bou)))
463 {
464 for (size_t pnl = 0; pnl < getAirfoil(bou).getNumberOfPanels(); ++pnl)
465 {
466 double sumGam = getBoundary(bou).sheets.freeVortexSheet(pnl, 0) + getBoundary(bou).sheets.attachedVortexSheet(pnl, 0);
467 Point2D rpnl = 0.5 * (getAirfoil(bou).getR(pnl) + getAirfoil(bou).getR(pnl + 1));
468 lambdaAdd[bou] += rpnl.kcross() * sumGam * getAirfoil(bou).len[pnl];
469 muAdd[bou] += (rpnl - getAirfoil(bou).rcm).length2() * sumGam * getAirfoil(bou).len[pnl];
470 }
471 lambdaAdd[bou] *= getPassport().physicalProperties.rho;
472 muAdd[bou] *= 0.5 * getPassport().physicalProperties.rho;
473 info('i') << "Added Masses for airfoil #" << bou << " = { " << lambdaAdd[bou][0] << ", " << lambdaAdd[bou][1] << ", " << muAdd[bou] << " }" << std::endl;
474
475 char direction;
476 switch ((int)(getPassport().physicalProperties.timeAccel))
477 {
478 case 0:
479 direction = 'x';
480 break;
481 case 1:
482 direction = 'y';
483 break;
484 case 2:
485 direction = 'w';
486 break;
487 default:
488 direction = '?';
489 info('e') << "Wrong dirfection is specified!" << std::endl;
490 exit(1);
491 }
492
493 std::string addMassFileName = getPassport().dir + "addMass-" + std::to_string(getAirfoil(bou).numberInPassport);
494 std::ofstream addMassFile;
495 if (!VMlib::fileExistTest(addMassFileName, info) || getPassport().physicalProperties.timeAccel == 0)
496 {
497 addMassFile.open(addMassFileName);
498 addMassFile << "Added Masses for airfoil:" << std::endl;
499 }
500 else
501 addMassFile.open(addMassFileName, std::ios_base::app);
502
503 addMassFile << direction << "-direction: " << lambdaAdd[bou][0] << " " << lambdaAdd[bou][1] << " " << muAdd[bou] << std::endl;
504
505 addMassFile.close();
506 }
507 }
508 currentTime = getPassport().timeDiscretizationProperties.timeStop;
509 }
510
511 timerInitialBuild.start();
512 if (nTotPan > 0 && getPassport().numericalSchemes.velocityComputation.second == 1)
513 {
514 auto& afl = getAirfoil(0);
515#ifdef USE_CUDA
516 getCuda().inflTreePnlVortex->UpdatePanelFreeAndAttachedVortexIntensity(afl.devFreeVortexSheetPtr, afl.devFreeVortexSheetLinPtr, afl.devAttachedVortexSheetPtr, afl.devAttachedVortexSheetLinPtr);
517 getCuda().inflTreePnlVortex->UpwardTraversal(getPassport().numericalSchemes.nbodyMultipoleOrder);
518#endif
519 }
520 timerInitialBuild.stop();
521
522 //Вычисление скоростей вихрей: для тех, которые в следе, и виртуальных, а также в точках wakeVP
523 CalcVortexVelo();
524
525
526#ifdef TURB
527 getTimers().start("Nut");
528 std::vector<double> nut;
529 CalcVeloDifference(nut);
530 getTimers().stop("Nut");
531
532 for (size_t q = 0; q < getWake().vtx.size(); ++q)
533 getNonConstWake().vtx[q].sigma() = 4.48364 * sqrt(passport.timeDiscretizationProperties.dt * (passport.physicalProperties.nu + nut[q]));
534#endif
535
536
537//#include "gammaCirc.h"
538
539 //Расчет и сохранение поля давления
540 if (ifDivisible(passport.timeDiscretizationProperties.saveVPstep))
541
542#ifdef OPTIMIZER
543 if ((getCurrentStep() >= OPTIMIZER_START_STEP && getCurrentStep() <= OPTIMIZER_STOP_STEP) || getCurrentStep()==0) //2026-03-28
544#endif
545 {
546 const double& vRef = getPassport().physicalProperties.vRef;
547 //scaleV = 1.0 / vRef;
548 double scaleP = 1.0 / (0.5 * getPassport().physicalProperties.rho * sqr(vRef));
549
550#ifdef USE_CUDA
551 measureVP->GPUCalcPressure();
552#else
553 measureVP->CalcPressure();
554#endif
555
556 if (mechanics.size()>0)
557 {
558 //std::cout << "presForces are calculated\n";
559 //MechanicsDeformable* ptr = dynamic_cast<MechanicsDeformable*>(mechanics[0].get());
560 Mechanics* ptr = mechanics[0].get();
561
562 if (measureVP->getTotalNumberOfRealPoints() > 0)
563 measureVP->SaveVP();
564
565 //if (ptr && !ptr->beam->fsi)
566
567 //сохранение главного вектора силы, вычисленного как интеграл от давления
568 //*
569 if (measureVP->elasticPoints.size() > 0)
570 {
571 std::ofstream presForcesFile;
572 if (currentStep == 0)
573 {
574 presForcesFile.open(getPassport().dir + "presForcesFile.csv");
575 presForcesFile << "time,Fx,Fy,Py" << std::endl;
576 }
577 else
578 presForcesFile.open(getPassport().dir + "presForcesFile.csv", std::ios_base::app);
579
580 auto r = measureVP->GetVPinElasticPoints();
581 Point2D presForce = { 0.0, 0.0 };
582 double yPower = 0.0;
583
584 for (int q = 0; q < r.size(); ++q)
585 {
586 presForce -= r[q].second * getAirfoil(0).len[q] * getAirfoil(0).nrm[q];
587 yPower -= (r[q].second * getAirfoil(0).len[q] * getAirfoil(0).nrm[q])[1] * (0.5 * (getAirfoil(0).getV(q) + getAirfoil(0).getV(q + 1))[1]);
588 }
589
590 presForcesFile << currentTime << "," << scaleP * presForce[0] << "," << scaleP * presForce[1] << "," << scaleP * yPower << std::endl;
591 presForcesFile.close();
592 }
593 }
594 //*/
595
596 if (mechanics.size() > 0)
597 {
598 //Коэффициенты разложения силы по балочным функциям
599 MechanicsDeformable* ptr = dynamic_cast<MechanicsDeformable*>(mechanics[0].get());
600 if (ptr && ptr->beam->fsi)
601 {
602 auto r = measureVP->GetVPinElasticPoints();
603
604 /*
605 std::ofstream elastPresFile;
606 elastPresFile.open(getPassport().dir + "elastPresFile-" + std::to_string(currentStep) + ".txt");
607 //if (currentStep == 0)
608 // elastPresFile.open(getPassport().dir + "elastPresFile.txt");
609 //else
610 // elastPresFile.open(getPassport().dir + "elastPresFile.txt", std::ios_base::app);
611 //elastPresFile << currentStep << " ";
612
613 for (size_t q = 0; q < r.size(); ++q)
614 elastPresFile << measureVP->elasticPoints[q][0] << " " << measureVP->elasticPoints[q][1] << " " << r[q].second << std::endl;
615 //elastPresFile << std::endl;
616 elastPresFile.close();
617 */
618
619 //Сохраняем давление на последних nLastSteps шагах по дисциплине очереди
620 std::vector<double> currentPres(ptr->chord.size());
621 for (size_t j = 0; j < ptr->chord.size(); ++j)
622 {
623 //currentPres[j] = -(ptr->beam->rho * 2 * ptr->beam->F); // gravity force for g = 2.0;
624 currentPres[j] = -(r[2 * j + 0].second - r[2 * j + 1].second);
625 }
626
627 if (ptr->beam->presLastSteps.size() < ptr->beam->nLastSteps)
628 ptr->beam->presLastSteps.push_back(currentPres);
629 else
630 {
631 for (int w = 1; w < ptr->beam->nLastSteps; ++w)
632 ptr->beam->presLastSteps[w - 1] = std::move(ptr->beam->presLastSteps[w]);
633 ptr->beam->presLastSteps.back() = currentPres;
634 }
635
636 for (int q = 0; q < ptr->beam->R; ++q)
637 {
638 ptr->beam->qCoeff[q] = 0;
639
640 if (ptr->beam->presLastSteps.size() == ptr->beam->nLastSteps)
641 {
642 for (size_t j = 0; j < ptr->chord.size(); ++j)
643 {
644 double averpres = 0.0;
645
646 //осредняем давление
647 for (int i = 0; i < ptr->beam->nLastSteps; ++i)
648 averpres += ptr->beam->presLastSteps[i][j];
649 averpres /= ptr->beam->nLastSteps;
650
651 ptr->beam->qCoeff[q] += -averpres * (ptr->initialChord[j].beg - ptr->initialChord[j].end).length() * ptr->beam->shape(q, 0.5 * (ptr->initialChord[j].beg + ptr->initialChord[j].end)[0]);
652 }
653 ptr->beam->qCoeff[q] /= (ptr->beam->intSqUnitShape * ptr->beam->L);
654 }
655
656 //std::cout << "q[" << q << "] = " << ptr->beam->qCoeff[q] << std::endl;
657 }
658
659 std::ofstream phiFile;
660 if (currentStep == 0)
661 {
662 phiFile.open(getPassport().dir + "phiFile.csv");
663 phiFile << "time";
664 for (int p = 0; p < ptr->beam->R; ++p)
665 phiFile << ",phi-" << std::to_string(p + 1);
666 for (int p = 0; p < ptr->beam->R; ++p)
667 phiFile << ",q-" << std::to_string(p + 1);
668 phiFile << std::endl;
669 }
670 else
671 phiFile.open(getPassport().dir + "phiFile.csv", std::ios_base::app);
672
673 phiFile << currentTime;
674 for (int p = 0; p < ptr->beam->R; ++p)
675 phiFile << "," << ptr->beam->phi(p, 0);
676 for (int p = 0; p < ptr->beam->R; ++p)
677 phiFile << "," << ptr->beam->qCoeff[p];
678
679 phiFile << std::endl;
680
681 phiFile.close();
682 }
683 }
684//*/
685 //getInfo('e') << "Deformable airfoils are not supported now!";
686 }
687
688 //Вычисление сил, действующих на профиль и сохранение в файл
689 for (auto& mech : mechanics)
690 {
691 mech->GetHydroDynamForce();
692 mech->GenerateForcesString();
693 mech->GeneratePositionString();
694 }
695
696 //Движение вихрей (сброс вихрей + передвижение пелены)
697#ifdef TURB
698 WakeAndAirfoilsMotion(true, &nut);
699#else
700 WakeAndAirfoilsMotion(true);
701#endif
702 wake->Restruct();
703 wake->SaveKadrVtk();
704
705 //std::string fname = getPassport().dir + "/airfoil-" + std::to_string(getAirfoil(0).numberInPassport) + "-" + std::to_string(currentStep) + ".pfl";
707 //{
708 // std::ofstream of(fname);
709 // for (size_t i = 0; i < getAirfoil(0).getNumberOfPanels(); ++i)
710 // of << getAirfoil(0).getR(i)[0] << " " << getAirfoil(0).getR(i)[1] << std::endl;
711 // of.close();
712 //}
713
714 }//if (!semiImplicitStrategy)
715
716
717//СОПРЯЖЕННАЯ ПОСТАНОВКА
718 if (semiImplicitStrategy)
719 {
720 CalcPanelsVeloAndAttachedSheets();
721 measureVP->Initialization();
722
723 CalcAndSolveLinearSystem();
724
725 //Вычисление скоростей вихрей: для тех, которые в следе, и виртуальных, а также в точках wakeVP
726 CalcVortexVelo();
727
728 //for (size_t s = 0; s < mechanics.size(); ++s) {// проверка: вычисление полной силы по циркуляциям новых вихрей
729 // mechanics[s]->GetHydroDynamForce();
730 // mechanics[s]->GenerateForcesString();
731 //}
732
733 //Расчет и сохранение поля давления
734 if (ifDivisible(passport.timeDiscretizationProperties.saveVPstep))
735 {
736 measureVP->CalcPressure();
737 measureVP->SaveVP();
738 }
739
740 WakeAndAirfoilsMotion(false); //профиль - кинематически, здесь Vold := V;
741 wake->Restruct();
742
743 CalcPanelsVeloAndAttachedSheets();
744 CalcAndSolveLinearSystem();
745
746 //Вычисление сил, действующих на профиль и сохранение в файл
747
748 for (size_t m = 0; m < mechanics.size(); ++m)
749 {
750 mechanics[m]->GetHydroDynamForce();
751
752 MechanicsRigidOscillPart* mechVar = dynamic_cast<MechanicsRigidOscillPart*>(mechanics[m].get());
753 if (mechVar)
754 {
755 if (getPassport().airfoilParams[m].addedMass.length2() == 0)
756 {
757 getInfo('e') << "Added mass of the airfoil should be non-zero!" << std::endl;
758 exit(1);
759 }
760 mechVar->MoveOnlyVelo();
761 Point2D accel = (mechVar->getV() - mechVar->getVOld()) * (1.0 / getPassport().timeDiscretizationProperties.dt);
762 mechVar->hydroDynamForce -=
763 Point2D{ accel[0] * getPassport().airfoilParams[m].addedMass[0],
764 accel[1] * getPassport().airfoilParams[m].addedMass[1] };
765 mechVar->GenerateForcesString();
766 mechVar->GeneratePositionString();
767 }
768 else
769 exit(3333);
770 }
771 }
772
773
774 wake->SaveKadrVtk();
775//*/
777
778
779 //Засечка времени в конце шага
780 getTimers().stop("Step");
781 /*
782 std::cout << "nvt = " << nVtxBeforeMerging << std::endl;
783 std::cout << "tIBT = " << timerInitialBuild.duration() << std::endl;
784 std::cout << "tRHS = " << timerRhs.duration() << std::endl;
785 std::cout << "tFIL = " << timerFillMatrix.duration() << std::endl;
786 std::cout << "tLIN = " << timerSlaeSolve.duration() << std::endl;
787 std::cout << "tVEL = " << timerConvVelo.duration() << std::endl;
788 std::cout << "tINS = " << timerInside.duration() << std::endl;
789 std::cout << "tMRG = " << timerMerging.duration() << std::endl;
790 */
791
792 info('i') << "Step = " << getCurrentStep() \
793 << " PhysTime = " << getCurrentTime() \
794 << std::setprecision(3) \
795 << " StepTime = " << getTimers().durationStep() \
796 << std::setprecision(6) \
797 << std::endl;
798
799 getTimers().GenerateStatString(getCurrentStep(), getCurrentTime(), getWake().vtx.size());
800
801 currentTime += passport.timeDiscretizationProperties.dt;
802 currentStep += 1;
803 //}
804 //catch (...)
805 //{
806 // info('e') << "!!! Exception from unknown source !!!" << std::endl;
807 // exit(-1);
808 //}
809}//Step()
810
811
812//Проверка проникновения вихрей внутрь профиля
813void World2D::CheckInside(std::vector<Point2D>& newPos, const std::vector<std::unique_ptr<AirfoilGeometry>>& oldAirfoil)
814{
815 getTimers().start("Inside");
816
817 gabb = 0;
818 check01 = 0;
819 check02 = 0;
820 checkPan = 0;
821
822#if (!defined(USE_CUDA))
823 nVtxBeforeMerging = newPos.size();
824 for (size_t afl = 0; afl < airfoil.size(); ++afl)
825 wake->Inside(newPos, *airfoil[afl], mechanics[afl]->isMoves, *oldAirfoil[afl]);
826
827
828 //Тренируемся работать с древом панелей
829 /*
830 int nTotPan = 0;
831 for (size_t s = 0; s < getNumberOfAirfoil(); ++s)
832 nTotPan += (int)getAirfoil(s).getNumberOfPanels();
833
834 auxTreePnl.reset(new CpuTreeInfo(tree_T::aux, object_T::panel, scheme_T::noScheme));
835
836 std::vector<std::pair<Point2D, Point2D>> pnls(nTotPan);
837 int counter = 0;
838 for (size_t s = 0; s < getNumberOfAirfoil(); ++s)
839 for (size_t p = 0; p < getAirfoil(s).getNumberOfPanels(); ++p)
840 pnls[counter++] = { getAirfoil(s).getR(p), getAirfoil(s).getR(p + 1) };
841
842 auxTreePnl->UpdatePanelGeometry(pnls, 0);
843 auxTreePnl->Build();
844 auxTreePnl->UpwardTraversal(0);
845
846 std::unique_ptr<CpuTreeInfo> cntrTreePnt;
847 cntrTreePnt.reset(new CpuTreeInfo(tree_T::contr, object_T::point2, scheme_T::noScheme));
848
849 std::vector<Vortex2D> newPosVtx(newPos.size());
850 for (size_t i = 0; i < newPos.size(); ++i)
851 newPosVtx[i].r() = newPos[i];
852
853 cntrTreePnt->Update(newPosVtx);
854 float tBuild = cntrTreePnt->Build();
855 float tUpward = cntrTreePnt->UpwardTraversal(0);
856
857 std::vector<std::pair<int,double>> closeIndexDist(cntrTreePnt->object.size());
858
859
860 VMlib::vmTimer timerA;
861 timerA.start();
862 double tDownward = auxTreePnl->DownwardTraversalClosestPanelToPoints(*cntrTreePnt, closeIndexDist, false, nullptr);
863 timerA.stop();
864
865 std::ofstream treeTimeFile;
866 if (getCurrentStep() == 0)
867 {
868 treeTimeFile.open(passport.dir + "/dbg/treeTime.csv");
869 treeTimeFile << "step,time,N,tBUI,tUPW,tDNW\n";
870 }
871 else
872 treeTimeFile.open(passport.dir + "/dbg/treeTime.csv", std::ios::app);
873
874 treeTimeFile << getCurrentStep() << ',' << getCurrentTime() << ',' << newPos.size() << ',' << tBuild << ',' << tUpward << ',' << tDownward << '\n';
875
876 treeTimeFile.close();
877
878 //std::cout << "timeFastNeib = " << timerA.duration() << ", N = " << newPos.size() << "\n";
879 /*
880 std::ofstream closeFile(passport.dir + "closeFile.txt");
881 for (size_t i = 0; i < newPos.size(); ++i)
882 closeFile << newPos[i][0] << " " << newPos[i][1] << " " << closeIndexDist[i].first << "\n";
883 closeFile.close();
884 */
885
886 //*/
887
888 //std::cout << "gab: " << gabb << " check1: " << check01 << " check2: " << check02 << " checkPan: " << checkPan << std::endl;
889#else
890// //////////////////////// CUDA ///////////////////////
891 timerInside.reset();
892 timerInside.start();
893
894 if ((newPos.size() > 0) && (getNumberOfAirfoil() > 0))
895 {
896 int nTotPanels = 0;
897 for (size_t afl = 0; afl < airfoil.size(); ++afl)
898 nTotPanels += (int)airfoil[afl]->getNumberOfPanels();
899
900 nVtxBeforeMerging = newPos.size();
901
902 auto& auxTree = *getNonConstCuda().auxTreePnl;
903
904 if ((currentStep == 0) || isAnyMovableOrDeformable())
905 {
906 auxTree.MemoryAllocate((int)getCuda().n_CUDA_pnls);
907 auxTree.UpdatePanelGeometry(nTotPanels, (double4*)airfoil[0]->devRPtr);
908 auxTree.Build();
909 auxTree.UpwardTraversal(0);
910 }
911
912 std::vector<int> hit(newPos.size());
913
914 if (isAnyMovableOrDeformable()) //Если профили подвижны - метод псевдонормалей
915 {
916 double* devNewpos_ptr;
917 cudaMalloc(&devNewpos_ptr, newPos.size() * sizeof(double) * 2);
918 cudaMemcpy(devNewpos_ptr, newPos.data(), newPos.size() * sizeof(double) * 2, cudaMemcpyHostToDevice);
919
920 auto& cntrTreePnt = *getCuda().cntrTreePoint;
921 cntrTreePnt.MemoryAllocate((int)getCuda().n_CUDA_wake);
922 cntrTreePnt.Update((int)newPos.size(), devNewpos_ptr);
923 cntrTreePnt.Build();
924
925 BHcu::treeClosestPanelToPointsCalculationWrapper(auxTree, cntrTreePnt, getWake().devNearestPanelPtr, true, getAirfoil(0).devPsnPtr);
926
927 cudaMemcpy(hit.data(), getWake().devNearestPanelPtr, newPos.size() * sizeof(int), cudaMemcpyDeviceToHost);
928 cudaFree(devNewpos_ptr);
929
930
931 //std::ofstream vortexFile(getPassport().dir + "vortexFile.txt");
932 //for (size_t i = 0; i < getWake().vtx.size(); ++i)
933 // vortexFile << getWake().vtx[i].r()[0] << " " << getWake().vtx[i].r()[1] << " " << hit[i] << '\n';
934 //vortexFile.close();
935
936 //std::ofstream panelFile(getPassport().dir + "panelFile.txt");
937 //for (size_t i = 0; i < getAirfoil(0).getNumberOfPanels(); ++i)
938 // panelFile << getAirfoil(0).getR(i)[0] << " " << getAirfoil(0).getR(i)[1] << " " << \
939 // getAirfoil(0).getR(i + 1)[0] << " " << getAirfoil(0).getR(i + 1)[1] << " " << '\n';
940 //panelFile.close();
941
942
943 }
944 else //иначе (для неподвижного профиля - метод трассировки лучей
945 {
946 std::vector<std::pair<Point2D, Point2D>> segments(newPos.size());
947 for (size_t i = 0; i < newPos.size(); ++i)
948 {
949 segments[i].first = wake->vtx[i].r();
950 segments[i].second = newPos[i];
951 }
952 double* devSegments_ptr;
953
954 cudaMalloc(&devSegments_ptr, newPos.size() * sizeof(double) * 4);
955 cudaMemcpy(devSegments_ptr, segments.data(), newPos.size() * sizeof(double) * 4, cudaMemcpyHostToDevice);
956
957 auto& cntrTreeSeg = *getCuda().cntrTreeSegment;
958 cntrTreeSeg.MemoryAllocate((int)getCuda().n_CUDA_wake);
959 cntrTreeSeg.UpdatePanelGeometry((int)newPos.size(), (double4*)devSegments_ptr);
960 cntrTreeSeg.Build();
961
962 BHcu::treePanelsSegmentsIntersectionCalculationWrapper(auxTree, cntrTreeSeg, getWake().devNearestPanelPtr);
963 cudaMemcpy(hit.data(), getWake().devNearestPanelPtr, newPos.size() * sizeof(int), cudaMemcpyDeviceToHost);
964 cudaFree(devSegments_ptr);
965 }
966
967 size_t panCounter = 0;
968 for (size_t afl = 0; afl < airfoil.size(); ++afl)
969 {
970 std::vector<double> gamma(airfoil[afl]->getNumberOfPanels(), 0.0);
971 for (int i = 0; i < newPos.size(); ++i)
972 {
973 if (hit[i] != -1)
974 {
975 if ((hit[i] >= panCounter) && (hit[i] < panCounter + airfoil[afl]->getNumberOfPanels()))
976 {
977 if (fabs(wake->vtx[i].g()) > 1.0)
978 std::cout << "Too large gamma is through: i = " << i << ", " << hit[i] << ", " << "gamma[hit[i] - panCounter] += " << wake->vtx[i].g() << std::endl;
979
980 gamma[hit[i] - panCounter] += wake->vtx[i].g();
981 wake->vtx[i].g() = 0.0;
982 }
983 }
984 }
985 airfoil[afl]->gammaThrough = gamma;
986 panCounter += airfoil[afl]->getNumberOfPanels();
987 }//for afl
988 }
989 timerInside.stop();
990
991#endif
992
993 getTimers().stop("Inside");
994}
995
996
997
998
999
1000//Решение системы линейных алгебраических уравнений
1001void World2D::SolveLinearSystem()
1002{
1003 getTimers().start("Solve");
1004/*
1005 if (currentStep == 0)
1006 {
1007 invMatr = matr;
1008 std::ifstream fileMatrix("matrix.txt");
1009 int nx, ny;
1010 fileMatrix >> nx;
1011 fileMatrix >> ny;
1012 for (int i = 0; i < matr.rows(); ++i)
1013 {
1014 for (int j = 0; j < matr.cols(); ++j)
1015 {
1016 fileMatrix >> matr(i, j);
1017 }
1018 }
1019 fileMatrix.close();
1020 }
1021*/
1022
1023//*
1024 //if (currentStep == 0)
1025 //{
1026 // std::ofstream fileMatrix(passport.dir + "/dbg/matrix.txt");
1027 // //int nx, ny;
1028 // //fileMatrix >> nx;
1029 // //fileMatrix >> ny;
1030 // fileMatrix.precision(16);
1031 // for (int i = 0; i < matrReord.rows(); ++i)
1032 // {
1033 // for (int j = 0; j < matrReord.cols(); ++j)
1034 // {
1035 // fileMatrix << matrReord(i, j) << " ";
1036 // }
1037 // fileMatrix << std::endl;
1038 // }
1039 // fileMatrix.close();
1040 //}
1041
1042
1043
1044 //{
1045 // std::ofstream fileMatrix(passport.dir + "/dbg/rhs"+std::to_string(currentStep)+".txt");
1046 // //int nx, ny;
1047 // //fileMatrix >> nx;
1048 // //fileMatrix >> ny;
1049 // fileMatrix.precision(16);
1050 // for (int i = 0; i < rhsReord.size(); ++i)
1051 // {
1052 // fileMatrix << rhsReord(i) << std::endl;
1053 // }
1054 // fileMatrix.close();
1055 //}//*/
1056
1057 //exit(-100500);
1058
1060 const int linSystemScheme = passport.numericalSchemes.linearSystemSolver.second;
1061
1062 if (linSystemScheme == 0) // Gaussian elimination
1063 {
1064 double t1 = -omp_get_wtime();
1065 if (useInverseMatrix && (currentStep == 0))
1066 {
1067 info('t') << "Inverting matrix... ";
1068
1069#if (defined(USE_CUDA))
1070 invMatr.resize(matrReord.rows(), matrReord.cols());
1071 for (int i = 0; i < (int)matrReord.rows(); ++i)
1072 for (int j = 0; j < (int)matrReord.cols(); ++j)
1073 invMatr(i, j) = (i == j) ? 1.0 : 0.0;
1074 cuInverseMatrix((int)matrReord.rows(), matrReord.data(), invMatr.data());
1075#else
1076 invMatr = matrReord.inverse();
1077#endif
1078 info('t') << "done" << std::endl;
1079 }
1080
1081 if (currentStep == 0)
1082 info('t') << "Solving system at step #0... ";
1083
1084 if (useInverseMatrix)
1085 {
1086 /*
1087 std::cout << "Solution via invMatrix" << std::endl;
1088 for (int i = 0; i < invMatr.rows(); ++i)
1089 for (int j = 0; j < invMatr.cols(); ++j)
1090 if (fabs(invMatr(i, j)) > 1e+5)
1091 std::cout << "invMatrixError: " << "invA(" << i << ", " << j << ") = " << invMatr(i, j) << std::endl;
1092
1093 for (int i = 0; i < rhsReord.size(); ++i)
1094 if (fabs(rhsReord(i)) > 1e+5)
1095 std::cout << "rhsReordError: " << "rhsReord(" << i << ") = " << rhsReord(i) << std::endl;
1096 */
1097 timerSlaeSolve.reset();
1098 timerSlaeSolve.start();
1099 sol = invMatr * rhsReord;
1100 timerSlaeSolve.stop();
1101 }
1102 else
1103 {
1104 timerSlaeSolve.reset();
1105 timerSlaeSolve.start();
1106 sol = matrReord.partialPivLu().solve(rhsReord);
1107 timerSlaeSolve.stop();
1108 }
1109
1110 if (currentStep == 0)
1111 info('t') << "done" << std::endl;
1112
1113 t1 += omp_get_wtime();
1114
1115 //std::cout << " Time in Gauss = " << t1 << std::endl;
1116/*
1117 std::ofstream solFile(getPassport().dir + "/sol" + std::to_string(currentStep) + "-Gauss.txt");
1118 solFile.precision(16);
1119 for (int i = 0; i < sol.size(); ++i)
1120 solFile << sol(i) << std::endl;
1121 solFile.close();
1122 exit(10101);
1123//*/
1124 }//Gauss
1125
1126/*
1127 if (currentStep == 0)
1128 { N
1129 invMatr = matr;
1130 std::ifstream fileMatrix("invMatrix.txt");
1131 int nx, ny;
1132 fileMatrix >> nx;
1133 fileMatrix >> ny;
1134 for (int i = 0; i < invMatr.rows(); ++i)
1135 {
1136 for (int j = 0; j < invMatr.cols(); ++j)
1137 {
1138 fileMatrix >> invMatr(i, j);
1139 }
1140 }
1141 fileMatrix.close();
1142 }
1143*/
1144
1145/*
1146 if (currentStep == 0)
1147 {
1148 Eigen::MatrixXd mmul = matr * invMatr;
1149 for (int i = 0; i < invMatr.rows(); ++i)
1150 mmul(i,i) -= 1.0;
1151
1152 double maxElem = 0.0;
1153 for (int i = 0; i < mmul.rows(); ++i)
1154 for (int j = 0; j < mmul.cols(); ++j)
1155 if (fabs(mmul(i,j)) > maxElem )
1156 maxElem = mmul(i,j);
1157
1158 std::cout << "A*A^-1: MAX ELEMENT = " << maxElem << std::endl;
1159 }
1160*/
1161
1162/*
1163 if (currentStep == 0)
1164 {
1165 std::ofstream fileMatrix("invMatrix4x3200.txt");
1166 fileMatrix << invMatr.rows() << " " << invMatr.cols() << std::endl;
1167 fileMatrix.precision(18);
1168 for (int i = 0; i < invMatr.rows(); ++i)
1169 {
1170 for (int j = 0; j < invMatr.cols(); ++j)
1171 {
1172 fileMatrix << invMatr(i, j) << " ";
1173 }
1174 fileMatrix << std::endl;
1175 }
1176 fileMatrix.close();
1177 }
1178*/
1179
1180
1181
1182 if ((linSystemScheme == 1 || linSystemScheme == 2)) // GMRES
1183 {
1184 int nFullVars = (int)getNumberOfBoundary();
1185 for (int i = 0; i < getNumberOfBoundary(); ++i)
1186 nFullVars += (int)(boundary[i]->GetUnknownsSize());
1187
1188 std::vector<std::vector<double>> Ggam(getNumberOfBoundary());
1189 for (int i = 0; i < getNumberOfBoundary(); ++i)
1190 Ggam[i].resize(boundary[i]->GetUnknownsSize());
1191
1192 std::vector<double> GR(getNumberOfBoundary());
1193
1194#ifndef USE_CUDA
1195 std::vector<double> Grhs(rhsReord.size());
1196 for (int i = 0; i < rhsReord.size(); ++i)
1197 Grhs[i] = rhsReord(i);
1198
1199 std::vector<int> Gpos(getNumberOfBoundary(), 0);
1200 for (int i = 1; i < getNumberOfBoundary(); ++i)
1201 Gpos[i] = Gpos[i - 1] + (int)(boundary[i - 1]->GetUnknownsSize());
1202
1203 std::vector<int> Gvsize(getNumberOfBoundary());
1204 for (int i = 0; i < getNumberOfBoundary(); ++i)
1205 Gvsize[i] = (int)(boundary[i]->GetUnknownsSize());
1206
1207 if (passport.numericalSchemes.linearSystemSolver.second == 1)
1208 {
1209 //for direct GMRES
1210 std::vector<double> Gmatr(nFullVars * nFullVars);
1211 for (size_t i = 0; i < nFullVars; ++i)
1212 for (size_t j = 0; j < nFullVars; ++j)
1213 Gmatr[i * nFullVars + j] = matrReord(i, j);
1214
1215 for (int i = 0; i < getNumberOfBoundary(); ++i)
1216 for (int j = 0; j < boundary[i]->GetUnknownsSize(); ++j)
1217 auto y_test = airfoil[i]->len[j % airfoil[i]->getNumberOfPanels()];
1218
1219 //Gmres->GMRES_Direct(*this, nFullVars, (int)getNumberOfBoundary(), Gmatr, Grhs, Gpos, Gvsize, Ggam, GR);
1220 std::cout << "GMRES_Direct should be modified!" << std::endl;
1221 exit(-100);
1222
1223
1224 sol.resize(nFullVars);
1225 int cntr = 0;
1226 for (int i = 0; i < getNumberOfBoundary(); ++i)
1227 for (int j = 0; j < boundary[i]->GetUnknownsSize(); ++j)
1228 sol(cntr++) = Ggam[i][j] /*/ airfoil[i]->len[j % airfoil[i]->getNumberOfPanels()]*/;
1229 for (int i = 0; i < getNumberOfBoundary(); ++i)
1230 sol(cntr++) = GR[i];
1231
1232 //std::ofstream solFile(getPassport().dir + "/dbg/sol-direct-gmres" + std::to_string(currentStep) + ".txt");
1233 //solFile.precision(16);
1234 //for (int i = 0; i < sol.size(); ++i)
1235 // solFile << sol(i) << std::endl;
1236 //solFile.close();
1237 }
1238 if (passport.numericalSchemes.linearSystemSolver.second == 2)
1239 {
1240 getInfo('e') << "Fast GMRES without CUDA is not implemented" << std::endl;
1241 exit(-2);
1242 }
1243#else
1244 if (passport.numericalSchemes.linearSystemSolver.second == 1)
1245 {
1246 getInfo('e') << "Direct GMRES for CUDA is not implemented" << std::endl;
1247 exit(-2);
1248 }
1249
1250 if ( (passport.numericalSchemes.linearSystemSolver.second == 2))
1251 {
1252 // for fast GMRES
1253 std::vector<std::vector<double>> GGrhs(getNumberOfBoundary());
1254 int cntrRhs = 0;
1255 for (int i = 0; i < getNumberOfBoundary(); ++i)
1256 GGrhs[i].resize(boundary[i]->GetUnknownsSize());
1257 for (int i = 0; i < getNumberOfBoundary(); ++i)
1258 for (int j = 0; j < boundary[i]->GetUnknownsSize(); ++j)
1259 GGrhs[i][j] = rhsReord(cntrRhs++);
1260
1261 std::vector<double> GrhsReg(getNumberOfBoundary());
1262 for (int i = 0; i < getNumberOfBoundary(); ++i)
1263 GrhsReg[i] = rhsReord[nFullVars - (getNumberOfBoundary() - i)];
1264
1265
1266 int niter;
1267 bool linScheme = (passport.numericalSchemes.boundaryCondition.second == 2);
1268
1269 double time_GMRES = -omp_get_wtime();
1270 getNonConstCuda().Gmres->GMRES(Ggam, GR, GGrhs, GrhsReg, niter);//, linScheme);
1271 time_GMRES += omp_get_wtime();
1272 //std::cout << "Time_GMRES = " << time_GMRES << std::endl;
1273
1274 sol.resize(nFullVars);
1275 int cntr = 0;
1276 for (int i = 0; i < getNumberOfBoundary(); ++i)
1277 for (int j = 0; j < boundary[i]->GetUnknownsSize(); ++j)
1278 sol(cntr++) = Ggam[i][j] / airfoil[i]->len[j % airfoil[i]->getNumberOfPanels()];
1279 for (int i = 0; i < getNumberOfBoundary(); ++i)
1280 sol(cntr++) = GR[i];
1281
1282 /*
1283 std::ofstream solFile(getPassport().dir + "/sol-fast-new-gmres" + std::to_string(currentStep) + ".txt");
1284 solFile.precision(16);
1285 for (int i = 0; i < sol.size(); ++i)
1286 solFile << sol(i) << std::endl;
1287 solFile.close();
1288 exit(1010110);
1289 //*/
1290 }
1291#endif
1292 }
1293
1294 getTimers().stop("Solve");
1295}//SolveLinearSystem()
1296
1297//Заполнение матрицы системы для всех профилей
1298void World2D::FillMatrixAndRhs()
1299{
1300 getTimers().start("MatRhs");
1301
1302 const int linSystemScheme = passport.numericalSchemes.linearSystemSolver.second;
1303
1304 if (linSystemScheme == 0 || linSystemScheme == 1)
1305 {
1306 timerFillMatrix.reset();
1307 timerFillMatrix.start();
1308
1309 Eigen::MatrixXd locMatr;
1310 Eigen::MatrixXd otherMatr;
1311 Eigen::VectorXd locLastLine, locLastCol;
1312
1313 std::vector<std::vector<Point2D>> locIQ;
1314 std::vector<std::vector<Point2D>> locOtherIQ;
1315
1316 //обнуляем матрицу на первом шаге расчета
1317 if (currentStep == 0)
1318 {
1319 for (int i = 0; i < matrReord.rows(); ++i)
1320 for (int j = 0; j < matrReord.cols(); ++j)
1321 matrReord(i, j) = 0.0;
1322 }
1323
1324 size_t currentRow = 0;
1325 size_t currentSkosRow = 0;
1326
1327 size_t nAllVars = 0;
1328 for (size_t bou = 0; bou < boundary.size(); ++bou)
1329 nAllVars += boundary[bou]->GetUnknownsSize();
1330
1331 for (size_t bou = 0; bou < boundary.size(); ++bou)
1332 {
1333 size_t nVars = boundary[bou]->GetUnknownsSize();
1334 if (currentStep == 0 || mechanics[bou]->isDeform)
1335 {
1336 locMatr.resize(nVars, nVars);
1337 locLastLine.resize(nVars);
1338 locLastCol.resize(nVars);
1339 }
1340
1341 if (currentStep == 0 || mechanics[bou]->isDeform)
1342 boundary[bou]->FillMatrixSelf(locMatr, locLastLine, locLastCol);
1343
1344
1345 //размазываем матрицу
1346 for (size_t i = 0; i < nVars; ++i)
1347 {
1348 if (currentStep == 0 || mechanics[bou]->isDeform)
1349 {
1350 for (size_t j = 0; j < nVars; ++j)
1351 matrReord(i + currentSkosRow, j + currentSkosRow) = locMatr(i, j);
1352
1353 matrReord(nAllVars + bou, i + currentSkosRow) = locLastLine(i);
1354 matrReord(i + currentSkosRow, nAllVars + bou) = locLastCol(i);
1355 }
1356 }
1357
1358 if ((currentStep == 0) || (!useInverseMatrix))
1359 {
1360 size_t currentCol = 0;
1361 size_t currentSkosCol = 0;
1362 for (size_t oth = 0; oth < boundary.size(); ++oth)
1363 {
1364 size_t nVarsOther = boundary[oth]->GetUnknownsSize();
1365
1366 if (bou != oth)
1367 {
1368 otherMatr.resize(nVars, nVarsOther);
1369
1370 boundary[bou]->FillMatrixFromOther(*boundary[oth], otherMatr);
1371
1372 //размазываем матрицу
1373 for (size_t i = 0; i < nVars; ++i)
1374 {
1375 for (size_t j = 0; j < nVarsOther; ++j)
1376 matrReord(i + currentSkosRow, j + currentSkosCol) = otherMatr(i, j);
1377 }
1378 }// if (bou != oth)
1379 currentCol += nVarsOther + 1;
1380 currentSkosCol += nVarsOther;
1381 }// for oth
1382 }// if (currentStep == 0 || mechanics[oth]->isMoves)
1383
1384 currentRow += nVars + 1;
1385 currentSkosRow += nVars;
1386 }// for bou
1387 timerFillMatrix.stop();
1388 }
1389
1390 velocity->FillRhs(rhsReord);
1391
1392 getTimers().stop("MatRhs");
1393}//FillMatrixAndRhs()
1394
1395//вычисляем размер матрицы и резервируем память под нее и под правую часть
1396void World2D::ReserveMemoryForMatrixAndRhs()
1397{
1398 //getTimers().start("Mem");
1399
1400 if (currentStep == 0)
1401 {
1402 dispBoundaryInSystem.resize(boundary.size());
1403 dispBoundaryInSystem[0] = 0;
1404
1405 for (size_t i = 1; i < boundary.size(); ++i)
1406 {
1407 dispBoundaryInSystem[i] = dispBoundaryInSystem[i - 1] + boundary[i - 1]->GetUnknownsSize() + 1;
1408 }
1409
1410 size_t matrSize = boundary.size();
1411 size_t matrSkosSize = 0;
1412
1413 for (auto it = boundary.begin(); it != boundary.end(); ++it)
1414 {
1415 matrSize += (*it)->GetUnknownsSize();
1416 matrSkosSize += (*it)->GetUnknownsSize();
1417 }
1418
1419 if ((getPassport().numericalSchemes.linearSystemSolver.second == 0 || getPassport().numericalSchemes.linearSystemSolver.second == 1 || (getPassport().numericalSchemes.linearSystemSolver.second == 2 && isAnyMovableOrDeformable())))
1420 {
1421 //matr.resize(matrSize, matrSize);
1422 matrReord.resize(matrSize, matrSize);
1423 //matrSkos.resize(matrSkosSize, matrSkosSize);
1424
1425 //matr.setZero();
1426 matrReord.setZero();
1427 //matrSkos.setZero();
1428
1429
1430 for (size_t i = 0; i < getNumberOfAirfoil(); ++i)
1431 {
1432 size_t nVari = boundary[i]->GetUnknownsSize();
1433 for (size_t j = 0; j < getNumberOfAirfoil(); ++j)
1434 {
1435 size_t nVarj = boundary[j]->GetUnknownsSize();
1436 IQ[i][j].first.resize(nVari, nVarj);
1437 IQ[i][j].second.resize(nVari, nVarj);
1438 }
1439 }
1440 }
1441
1442 //rhs.resize(matrSize);
1443 rhsReord.resize(matrSize);
1444 //rhsSkos.resize(matrSkosSize);
1445 }
1446 //rhs.setZero();
1447 rhsReord.setZero();
1448 //rhsSkos.setZero();
1449
1450 //getTimers().stop("Mem");
1451}//ReserveMemoryForMatrixAndRhs()
1452
1453
1454
1455// Вычисление скоростей (и конвективных, и диффузионных) вихрей (в пелене и виртуальных), а также в точках вычисления VP
1456void World2D::CalcVortexVelo()
1457{
1458 velocity->ResizeAndZero();
1459
1460 getTimers().start("ConvVel");
1461
1462 //Конвективные скорости всех вихрей (и в пелене, и виртуальных), индуцируемые вихрями в пелене
1463 //Подготовка CUDA
1464#if defined(__CUDACC__) || defined(USE_CUDA)
1465 cuda.RefreshWake(2);
1466
1467 cuda.RefreshAfls(2);
1468 cuda.RefreshVirtualWakes(2);
1469#endif
1470 getTimers().stop("ConvVel");
1471
1472 velocity->CalcConvVelo();
1473
1474 //Расчет средних значений eps для каждой панели и их передача на видеокарту
1475
1476 getTimers().start("DiffVel");
1477
1478 for (size_t bou = 0; bou < getNumberOfBoundary(); ++bou)
1479 getNonConstAirfoil(bou).calcMeanEpsOverPanel();
1480
1481#if defined(__CUDACC__) || defined(USE_CUDA)
1482 for (size_t i = 0; i < airfoil.size(); ++i)
1483 cuda.CopyMemToDev<double, 1>(airfoil[i]->getNumberOfPanels(), airfoil[i]->meanEpsOverPanel.data(), airfoil[i]->devMeanEpsOverPanelPtr);
1484#endif
1485
1486 getTimers().stop("DiffVel");
1487
1488 //Вычисление диффузионных скоростей вихрей (и в пелене, и виртуальных)
1489 velocity->CalcDiffVelo();
1490
1491 //Обнуление вязких напряжений, если они не были вычислены
1492 for (auto& afl : airfoil)
1493 {
1494 if (afl->viscousStress.size() == 0)
1495 afl->viscousStress.resize(afl->getNumberOfPanels(), 0.0);
1496 }
1497
1498 getTimers().start("Save");
1499
1500 /*
1501 //Сохранение всех параметров для вихрей в пелене
1502 if (!(currentStep % 1))
1503 {
1504 VMlib::CreateUserDirectory(passport.dir, "dbg");
1505 std::ostringstream sss, sss2;
1506 sss << "prmWake";
1507 sss << currentStep;
1508 sss2 << "testWake";
1509 sss2 << currentStep;
1510
1511 std::ofstream prmtFile(passport.dir + "dbg/" + sss.str());
1512 prmtFile.precision(11);
1513 prmtFile << "i x y g sigma epsast convVeloX convVeloY diffVeloX diffVeloY I0 I1 I2X I2Y I3X I3Y" << std::endl;
1514 for (size_t i = 0; i < wake->vtx.size(); ++i)
1515 prmtFile << i << " " \
1516 << wake->vtx[i].r()[0] << " " << wake->vtx[i].r()[1] << " " \
1517 << wake->vtx[i].g() << " " \
1518 << wake->vtx[i].sigma() << " " \
1519 << velocity->wakeVortexesParams.epsastWake[i] << " " \
1520 << velocity->wakeVortexesParams.convVelo[i][0] << " " << velocity->wakeVortexesParams.convVelo[i][1] << " "\
1521 << velocity->wakeVortexesParams.diffVelo[i][0] << " " << velocity->wakeVortexesParams.diffVelo[i][1] << " "\
1522 //<< std::endl;
1523 << velocity->wakeVortexesParams.I0[i] << " " \
1524 << velocity->wakeVortexesParams.I1[i] << " " \
1525 << velocity->wakeVortexesParams.I2[i][0] << " " << velocity->wakeVortexesParams.I2[i][1] << " " \
1526 << velocity->wakeVortexesParams.I3[i][0] << " " << velocity->wakeVortexesParams.I3[i][1] << " " \
1527 << "\n";
1528
1529 prmtFile.close();
1530
1531 std::ofstream prmtFile2(passport.dir + "dbg/" + sss2.str());
1532 prmtFile2 << wake->vtx.size() << "\n";
1533 prmtFile2.precision(16);
1534 for (size_t i = 0; i < wake->vtx.size(); ++i)
1535 prmtFile2 << wake->vtx[i].r()[0] << " " << wake->vtx[i].r()[1] << " " << wake->vtx[i].g() << "\n";
1536 prmtFile2.close();
1537 }
1538//*/
1539
1540/*
1541 //Сохранение всех параметров для виртуальных вихрей
1542 if (!(currentStep % 1))
1543 {
1544 for (size_t b = 0; b < boundary.size(); ++b)
1545 {
1546 std::ostringstream sss;
1547 sss << "prmVirtual_";
1548 sss << b << "-";
1549 sss << currentStep;
1550 std::ofstream prmFileVirt(passport.dir + "dbg/" + sss.str());
1551 //prmFileVirt.precision(16);
1552 prmFileVirt << "i x y g sigma epsast convVeloX convVeloY diffVeloX diffVeloY I0 I1 I2X I2Y I3X I3Y" << std::endl;
1553 for (size_t i = 0; i < boundary[b]->virtualWake.vtx.size(); ++i)
1554 prmFileVirt << i << " " \
1555 << boundary[b]->virtualWake.vtx[i].r()[0] << " " << boundary[b]->virtualWake.vtx[i].r()[1] << " " \
1556 << boundary[b]->virtualWake.vtx[i].g() << " " \
1557 << boundary[b]->virtualWake.vtx[i].sigma() << " " \
1558 << velocity->virtualVortexesParams[b].epsastWake[i] << " " \
1559 << velocity->virtualVortexesParams[b].convVelo[i][0] << " " << velocity->virtualVortexesParams[b].convVelo[i][1] << " "\
1560 << velocity->virtualVortexesParams[b].diffVelo[i][0] << " " << velocity->virtualVortexesParams[b].diffVelo[i][1] << " "\
1561 //<< std::endl;
1562 << velocity->virtualVortexesParams[b].I0[i] << " " \
1563 << velocity->virtualVortexesParams[b].I1[i] << " " \
1564 << velocity->virtualVortexesParams[b].I2[i][0] << " " << velocity->virtualVortexesParams[b].I2[i][1] << " " \
1565 << velocity->virtualVortexesParams[b].I3[i][0] << " " << velocity->virtualVortexesParams[b].I3[i][1] << " " \
1566 << std::endl;
1567 prmFileVirt.close();
1568 }
1569 //if (currentStep==3) exit(-123);
1570 }
1571//*/
1572
1573 getTimers().stop("Save");
1574
1575}//CalcVortexVelo()
1576
1577
1578#ifdef TURB
1579void World2D::CalcVeloDifference(std::vector<double>& nut) const
1580{
1581 nut.resize(getWake().vtx.size());
1582
1583 //std::vector<Point2D> basePoints = { {0.,0.}, {0.7960681614061415,0.28974511519885826},{0.4235790213720913,0.7336603860367653},{-0.14710745031841074,0.8342878085211184},{-0.6489607110877307,0.5445426933222601},{-0.8471580427441825,0.},{-0.6489607110877307, -0.5445426933222601},{-0.14710745031841083, -0.8342878085211185},{0.42357902137209114, -0.7336603860367653},{0.7960681614061413, -0.28974511519885826},{1.5783500941642636,0.3354886691271114},{1.3054390240228053,0.9484569686703798},{0.8068056870865843,1.3974284418694807},{0.16866831725352804,1.6047718422305033},{-0.4986333369362209,1.534635612175369},{-1.0797167572280424,1.1991469430482575},{-1.474107341276333,0.6563148735601236},{-1.6136113741731684, 0.0},{-1.4741073412763332, -0.6563148735601233},{-1.0797167572280428, -1.1991469430482573},{-0.49863333693622125, -1.5346356121753688},{0.1686683172535277, -1.6047718422305033},{0.8068056870865841, -1.397428441869481},{1.3054390240228053, -0.94845696867038},{1.5783500941642634, -0.3354886691271116},{2.374139932721701,0.3578440178114266},{2.1631871809578898,1.041736042079128},{1.7600257551753844,1.6330652474247764},{1.2004783172997366,2.0792894389479355},{0.5342631127196844,2.340759639009162},{-0.1794237140005648,2.394243072889942},{-0.8771679563438638,2.234987502781785},{-1.496971976377837,1.8771434849703577},{-1.9837634669573248,1.3525070308108138},{-2.2942888678950686,0.7076943915843854},{-2.4009566345994733,0.},{-2.2942888678950686, -0.7076943915843852},{-1.983763466957325, -1.3525070308108136},{-1.4969719763778373, -1.8771434849703577},{-0.877167956343864, -2.2349875027817845},{-0.1794237140005649, -2.394243072889942},{0.5342631127196842, -2.340759639009162},{1.2004783172997364, -2.0792894389479355},{1.7600257551753842, -1.6330652474247767},{2.1631871809578898, -1.041736042079128},{2.3741399327217017, -0.3578440178114266},{3.174096170885181,0.35763474709931864},{3.0149337350056395,1.054970941725745},{2.7045899323860105,1.6994064881590272},{2.258626697172506,2.2586266971725055},{1.6994064881590272,2.7045899323860105},{1.054970941725745,3.0149337350056395},{0.35763474709931853,3.174096170885181},{-0.35763474709931864,3.174096170885181},{-1.054970941725745,3.0149337350056395},{-1.6994064881590272,2.7045899323860105},{-2.2586266971725055,2.258626697172506},{-2.7045899323860105,1.6994064881590272},{-3.0149337350056395,1.054970941725745},{-3.174096170885181,0.35763474709931853},{-3.174096170885181, -0.35763474709931864},{-3.0149337350056395, -1.054970941725745},{-2.7045899323860105, -1.6994064881590272},{-2.258626697172506, -2.2586266971725055},{-1.6994064881590272, -2.7045899323860105},{-1.054970941725745, -3.0149337350056395},{-0.35763474709931853, -3.174096170885181},{0.35763474709931864, -3.174096170885181},{1.054970941725745, -3.0149337350056395},{1.6994064881590272, -2.7045899323860105},{2.2586266971725055, -2.258626697172506},{2.7045899323860105, -1.6994064881590272},{3.0149337350056395, -1.054970941725745},{3.174096170885181, -0.35763474709931853},{3.9723533560766993,0.3680927461358311},{3.8370796268482272,1.0917432814164258},{3.5711387529972023,1.7782158086501256},{3.1835870320217006,2.404133329255034},{2.687622072962485,2.948180972772273},{2.100133367655915,3.39183184921509},{1.4411271393461673,3.719977960086802},{0.7330450547474776,3.921444683149553},{0.0, 3.9893713107821425},{-0.7330450547474779,3.921444683149553},{-1.4411271393461675,3.719977960086802},{-2.1001333676559155,3.39183184921509},{-2.687622072962485,2.948180972772273},{-3.1835870320217,2.4041333292550338},{-3.5711387529972023,1.7782158086501254},{-3.8370796268482272,1.0917432814164256},{-3.9723533560766993,0.368092746135831},{-3.9723533560766993, -0.3680927461358311},{-3.8370796268482272, -1.0917432814164258},{-3.5711387529972023, -1.7782158086501256},{-3.1835870320217006, -2.404133329255034},{-2.687622072962485, -2.948180972772273},{-2.100133367655915, -3.39183184921509},{-1.4411271393461673, -3.719977960086802},{-0.7330450547474776, -3.921444683149553},{0.0, -3.9893713107821425},{0.7330450547474779, -3.921444683149553},{1.4411271393461675, -3.719977960086802},{2.1001333676559155, -3.39183184921509},{2.687622072962485, -2.948180972772273},{3.1835870320217, -2.4041333292550338},{3.5711387529972023, -1.7782158086501254},{3.8370796268482272, -1.0917432814164256},{3.9723533560766993, -0.368092746135831} };
1584
1585 //std::vector<Point2D> basePoints = { {0.,0.}, { -2.47356, -2.00346 }, { -1.79985,1.13964 }, { -2.95238,-0.054986 }, { -1.70333,-0.910983 }, { -2.71457,1.15701 }, { -2.01605,-3.05827 }, { -1.31647,1.96372 }, { -2.99071,1.97152 }, { -0.0334834,0.767391 }, { -3.2381,0.623172 }, { -1.9359,2.42903 }, { -3.52381,-0.678158 }, { 1.05002,-2.98592 }, { 2.01605,-3.05827 }, { 1.90533,-1.62888 }, { 0.666259,0.130419 }, { 3.52381,-0.678157 }, { 1.33293,0.258719 }, { -0.678158,3.52381 }, { 1.61226,1.64718 }, { -0.699742,0.318031 }, { 0.623172,3.2381 }, { -0.376369,-0.567049 }, { 1.13573,2.72822 }, { 3.2381,0.623172 }, { 2.47356,-2.00346 }, { 2.99071,1.97152 }, { 0.351358,-1.84626 }, { -3.13329,-1.62097 }, { -1.09647,-2.98592 }, { -1.36641,0.446331 }, { -0.054986,2.95238 }, { -0.518096,-3.44188 }, { -1.90533,-1.62888 }, { 3.13329,-1.62097 }, { 0.92973,-2.31605 }, { 1.95024,3.00436 }, { 2.71457,1.15701 }, { 2.95238,-0.0549855 }, { 1.99542,-0.214456 }, { -1.62097,3.13329 }, { -0.378359,2.09655 }, { -1.99542,-0.214456 }, { -0.750952,-1.13528 }, { -0.397807,-1.84626 }, { 0.750952,-1.13528 }, { -0.976178,-2.31605 }, { 0.471648,-3.44188 }, { 0.376369,-0.567049 }, { 1.70333,-0.910983 }, { 1.76636,0.952026 }, { 0.0948167,1.43406 }, { 0.930748,1.78001 }, { 2.26517,2.3001 } };
1586
1587 std::vector<Point2D> basePoints = { {0.,0.}, { -2.10684, -1.84373 }, { 0.494872,-3.43051 }, { -1.8307,-0.772619 }, { -2.97035,-1.77037 }, { -0.869454,-2.66541 }, { -1.81226,-2.74913 }, { -0.93585,0.525783 }, { -2.84757,1.56598 }, { -0.140883,1.06415 }, { -3.2381,0.623172 }, { -1.70686,2.57839 }, { -3.36087,-0.827564 }, { 0.869454,-2.66541 }, { 1.81226,-2.74913 }, { 1.8307,-0.772618 }, { 0.794967,0.794967 }, { 2.78945,-0.204392 }, { 1.04125,-0.439927 }, { -0.366572,3.2381 }, { 1.39868,1.39868 }, { -1.04125,-0.439927 }, { 0.623172,3.2381 }, { 0.0, -1.13645 }, { 1.22688,2.50847 }, { 3.2381,0.623172 }, { 2.97035,-1.77037 }, { 2.78461,2.04139 }, { 0.0, -1.99554 }, { -2.78945,-0.204392 }, { -0.494871,-3.43051 }, { 2.10684,-1.84373 }, { 3.36087,-0.827563 }, { -0.635755,2.30225 }, { -1.87866,1.46859 }, { 2.50847,1.22688 }, { 2.04139,2.78461 } };
1588
1589
1590 //std::vector<Point2D> basePoints(101);
1591 //basePoints[0] = { 0.0, 0.0 };
1592 //for (size_t i = 0; i < 100; ++i)
1593 // basePoints[i+1] = { cos(DPI * i / 100), sin(DPI * i / 100) };
1594
1595
1596 WakeDataBase controlPoints(*this);
1597 controlPoints.vtx.resize(getWake().vtx.size() * basePoints.size());
1598
1599 for (size_t q = 0; q < getWake().vtx.size(); ++q)
1600 for (size_t i = 0; i < basePoints.size(); ++i)
1601 controlPoints.vtx[q * basePoints.size() + i].r() = getWake().vtx[q].r() + basePoints[i] * getWake().vtx[q].sigma();
1602
1603
1604 cuReserveDevMem((void*&)controlPoints.devVtxPtr, controlPoints.vtx.size() * sizeof(Vortex2D), 333);
1605 cuCopyWakeToDev(controlPoints.vtx.size(), controlPoints.vtx.data(), controlPoints.devVtxPtr, 444);
1606 cuReserveDevMem((void*&)controlPoints.devVelPtr, controlPoints.vtx.size() * 2 * sizeof(double), 555);
1607
1608 auto& treeNut = getNonConstCuda().cntrTreeNut;
1609 treeNut->MemoryAllocate((int)getCuda().n_CUDA_wake * (int)basePoints.size());
1610 treeNut->Update((int)controlPoints.vtx.size(), controlPoints.devVtxPtr);
1611 treeNut->Build();
1612
1613 std::vector<Point2D> veloVar(getWake().vtx.size() * basePoints.size(), Point2D{0.0, 0.0});
1614 std::vector<double> epsVar;
1615 velocity->GPUCalcConvVeloToSetOfPointsFromWake(treeNut, controlPoints, veloVar, epsVar, true, false);
1616 if (getNumberOfAirfoil() > 0)
1617 velocity->GPUCalcConvVelocityToSetOfPointsFromSheets(treeNut, controlPoints, veloVar);
1618
1619 cuDeleteFromDev((void*&)controlPoints.devVtxPtr, 777);
1620 cuDeleteFromDev((void*&)controlPoints.devVelPtr, 888);
1621
1622 for (size_t q = 0; q < getWake().vtx.size(); ++q)
1623 {
1624 double dv2 = 0.0;
1625 for (size_t i = 1; i < basePoints.size(); ++i)
1626 dv2 += (veloVar[q * basePoints.size()] - veloVar[q * basePoints.size() + i]).length2() * pow(getWake().vtx[q].sigma() / (controlPoints.vtx[q * basePoints.size() + i].r() - getWake().vtx[q].r()).length(), 0.6666667);
1627
1628 dv2 *= 1.0 / (basePoints.size() - 1);
1629 nut[q] = 0.105 * 0.604 * wake->vtx[q].sigma() * sqrt(dv2);
1630 }
1631}
1632#endif
1633
1634
1635
1636// Вычисление скоростей панелей и интенсивностей присоединенных слоев вихрей и источников
1637void World2D::CalcPanelsVeloAndAttachedSheets()
1638{
1639 //вычисляем скорости панелей
1640 for (size_t i = 0; i < airfoil.size(); ++i)
1641 mechanics[i]->VeloOfAirfoilPanels(currentStep * passport.timeDiscretizationProperties.dt);
1642
1643 //вычисляем интенсивности присоединенных слоев
1644 for (size_t i = 0; i < airfoil.size(); ++i)
1645 boundary[i]->ComputeAttachedSheetsIntensity();
1646
1647}//CalcPanelsVeloAndAttachedSheets(...)
1648
1649
1650
1651
1652#include <random>
1653
1654//Вычисляем новые положения вихрей (в пелене и виртуальных)
1655void World2D::MoveVortexes(std::vector<Point2D>& newPos, std::vector<double>* nutPtr)
1656{
1657 //getTimers().start("Move");
1658
1659 size_t nvt = wake->vtx.size();
1660 size_t nVirtVortex = 0;
1661 for (size_t i = 0; i < getNumberOfBoundary(); ++i)
1662 nVirtVortex += boundary[i]->virtualWake.vtx.size();
1663
1664
1665
1666
1667 size_t counter = wake->vtx.size() - nVirtVortex;
1668 for (size_t bou = 0; bou < boundary.size(); ++bou)
1669 {
1670 for (int i = 0; i < (int)boundary[bou]->virtualWake.vtx.size(); ++i)
1671 {
1672 wake->vtx[counter].g() = boundary[bou]->virtualWake.vtx[i].g();
1673 wake->vtx[counter].sigma() = getPassport().wakeDiscretizationProperties.sigma0;
1674 ++counter;
1675 }
1676 }
1677
1678
1679
1680
1682#if defined (SPH)
1683
1684 const double c = 4.5;
1685
1686 const auto& eps = velocity->wakeVortexesParams.epsastWake;
1687
1688 double minEpsast = *std::min_element(eps.begin(), eps.end());
1689 double meanEpsast = std::accumulate(eps.begin(), eps.end(), 0.0) / eps.size();
1690
1691 double uniformH = minEpsast + 0.2 * (meanEpsast - minEpsast);
1692
1693 std::vector<std::vector<int>> neib(wake->vtx.size());
1694 for (auto& nb : neib)
1695 nb.reserve(100);
1696
1697#pragma omp parallel for schedule(dynamic, 100)
1698 for (int i = 0; i < (int)wake->vtx.size(); ++i)
1699 {
1700 const double hati = c * eps[i];
1701 const double hati2 = sqr(hati);
1702 const int dist = std::ceil(hati / uniformH);
1703
1704 Point2D vtx = wake->vtx[i].r();
1705 int cx = (int)(vtx[0] / uniformH);
1706 int cy = (int)(vtx[1] / uniformH);
1707
1708 for (int j = 0; j < (int)wake->vtx.size(); ++j)
1709 {
1710 Point2D test = wake->vtx[j].r();
1711 int tx = (int)(test[0] / uniformH);
1712 int ty = (int)(test[1] / uniformH);
1713
1714 if (abs(cx - tx) <= dist && abs(cy - ty) <= dist)
1715 if ((vtx - test).length2() < hati2)
1716 neib[i].push_back(j);
1717 }
1718 }
1719
1720 //Ядро иранца
1721 //auto W = [](double xi, double R) {double t = xi / R; return (t > 1 ? 0.0 : (8.0 / ( PI * R * R)) * (0.75 * t * t - 1.5 * t + 0.75)); };
1722 //auto gW = [](double xi, double R) {double t = xi / R; return (t > 1 ? 0.0 : (12.0 / (PI * sqr(R * R))) * (1.0 - 1.0/t)); };
1723
1724 //M4
1725 Point2D zeroVec{ 0.0, 0.0 };
1726 auto W = [](double xi, double R) {double t = xi / R; return (t > 1 ? 0.0 : 40.0 / (7.0 * PI * R * R) * (t > 0.5 ? 2.0 * cubPower(1.0 - t) : 1.0 - 6.0 * sqr(t) * (1.0 - t))); };
1727 auto WW = [zeroVec](const Point2D& v, double R) {double t = v.length() / R; return (t > 1 ? zeroVec : v * (240.0 / (7.0 * PI * sqr(sqr(R))) * (t > 0.5 ? -sqr(1.0 - t) / t : -2.0 + 3.0 * t))); };
1728
1729 //Poly6
1730 //auto W = [](double xi, double R) {double t = xi / R; return (t > 1 ? 0.0 : 4.0 / (PI * R * R) * cubPower(1.0 - t * t)); };
1731 //auto gW = [](double xi, double R) {double t = xi / R; return (t > 1 ? 0.0 : -24.0 / (PI * sqr(R * R)) * sqr(1.0 - t * t)); };
1732
1733 //Poly6 + grad Spicky
1734 //auto W = [](double xi, double R) {double t = xi / R; return (t > 1 ? 0.0 : 4.0 / (PI * R * R) * cubPower(1.0 - t * t)); };
1735 //auto gW = [](double xi, double R) {double t = xi / R; return (t > 1 ? 0.0 : -30.0 / (PI * sqr(R * R)) * sqr(1.0 - t) / t); };
1736
1737 //WendlandC6
1738 /*
1739 Point2D zeroVec{ 0.0, 0.0 };
1740 auto W = [](double xi, double R) {double t = xi / R;
1741 return (t > 1 ? 0.0 : 78.0 / (7.0 * PI * sqr(R)) * sqr(sqr(sqr((1.0 - t)))) * (32.0 * cubPower(t)+ 25.0 * sqr(t) + 8.0*t + 1.0)); };
1742 auto WW = [zeroVec](const Point2D& v, double R) -> Point2D {double t = v.length() / R;
1743 return (t > 1 ? zeroVec : v * (-1716.0 / (7.0 * PI * sqr(sqr(R))) * sqr(cubPower(1.0 - t)) * (1.0 - t) * (16.0 * sqr(t) + 7.0 * t + 1.0))); };
1744 */
1745
1746 //fp5
1747 /*
1748 auto W = [](double xi, double R) {
1749 double t = xi / R;
1750 if (t >= 1) return 0.0;
1751 if (t >= 2.0 / 3.0)
1752 return 1701.0 / (478.0 * PI) * sqr(sqr(1.0 - t)) * (1.0 - t);
1753 if (t >= 1.0 / 3.0)
1754 return 1701.0 / (478.0 * PI) * (sqr(sqr(1.0 - t)) * (1.0 - t) - 6.0 * sqr(sqr((2.0 / 3.0 - t))) * (2.0 / 3.0 - t));
1755
1756 return 1701.0 / (478.0 * PI) * (sqr(sqr(1.0 - t)) * (1.0 - t) - 6.0 * sqr(sqr((2.0 / 3.0 - t))) * (2.0 / 3.0 - t) + 15.0 * sqr(sqr((1.0 / 3.0 - t))) * (1.0 / 3.0 - t));
1757 };
1758
1759 auto gW = [](double xi, double R) {
1760 double t = xi / R;
1761 if (t >= 1) return 0.0;
1762 if (t >= 2.0 / 3.0)
1763 return -8505.0 / (478.0 * PI * R * R) * sqr(sqr(1.0 - t)) / t;
1764 if (t >= 1.0 / 3.0)
1765 return -8505.0 / (478.0 * PI * R * R) * (sqr(sqr(1.0 - t)) - 6.0 * sqr(sqr((2.0 / 3.0 - t)))) / t;
1766
1767 return -8505.0 / (478.0 * PI * R * R) * (sqr(sqr(1.0 - t)) - 6.0 * sqr(sqr((2.0 / 3.0 - t))) + 15.0 * sqr(sqr((1.0 / 3.0 - t)))) / t;
1768 };
1769 */
1770
1771
1772 std::vector<double> OIP(wake->vtx.size(), 0.0);
1773 std::vector<double> VIP(wake->vtx.size(), 0.0);
1774 std::vector<Point2D> GIP(wake->vtx.size(), { 0.0, 0.0 });
1775 std::vector<Point2D> DIFFV(wake->vtx.size(), { 0.0, 0.0 });
1776
1777
1778 const auto& vtx = wake->vtx;
1779 const auto& vel = velocity->wakeVortexesParams.convVelo;
1780
1781
1782//#ifdef TURB
1783// std::vector<Point2D> grNu(wake->vtx.size(), { 0.0, 0.0 });
1784//
1785// //std::vector<Point2D> grVx(wake->vtx.size(), { 0.0, 0.0 });
1786// //std::vector<Point2D> grVy(wake->vtx.size(), { 0.0, 0.0 });
1787// std::vector<nummatrix<double,2,2>> grV(wake->vtx.size());
1788// std::vector<Point2D> TP(wake->vtx.size(), { 0.0, 0.0 });
1789// std::vector<double> RIP(wake->vtx.size(), 0.0);
1790// std::vector<double> DELTAG(wake->vtx.size(), 0.0);
1791//#endif //TURB
1792
1793 std::vector<nummatrix<double, 2, 2>> mtrL(wake->vtx.size(), { {0.0, 0.0}, {0.0, 0.0} });
1794 std::vector<std::vector<double>> weight(wake->vtx.size());
1795
1796
1797 std::vector<Point2D> sumDIFFV(wake->vtx.size(), {0.0, 0.0});
1798 std::vector<double> sumWeight(wake->vtx.size(), 0.0);
1799
1800#pragma omp parallel for schedule(dynamic, 100)
1801 for (int i = 0; i < (int)vtx.size(); ++i)
1802 {
1803 const double hati = c * eps[i];
1804
1805 double resultOIP = 0.0;
1806 //double resultIVIP = 0.0;
1807
1808 const Point2D myPos = vtx[i].r();
1809
1810 double I1 = 0.0, I0 = 1.0;
1811
1812
1813 weight[i].reserve(neib[i].size());
1814
1815 for (auto& nbr : neib[i])
1816 {
1817 double dist2 = (myPos - vtx[nbr].r()).length2();
1818 double dist = sqrt(dist2);
1819 double w = W(dist, hati);
1820 I1 += vtx[nbr].g() * w;
1821 weight[i].push_back(w);
1822 }
1823
1824 OIP[i] = I1 / I0;
1825 }
1826
1827#pragma omp parallel for schedule(dynamic, 100)
1828 for (int i = 0; i < (int)vtx.size(); ++i)
1829 {
1830 const Point2D myPos = vtx[i].r();
1831 const double hati = c * eps[i];
1832
1833 nummatrix<double, 2, 2> M = { {0.0, 0.0}, {0.0, 0.0} };
1834 GIP[i].toZero();
1835
1836 for (auto& nbr : neib[i])
1837 {
1838 Point2D dr = myPos - vtx[nbr].r();
1839 Point2D grweight = WW(dr, hati);
1840 double vol = (vtx[nbr].g() / OIP[nbr]);
1841 if (vol != vol)
1842 vol = 0.0;
1843
1844 M += (dr | grweight) * vol;
1845 GIP[i] -= grweight * ((OIP[nbr] - OIP[i]) * vol);
1846 }
1847
1848 double detM = M[0][0] * M[1][1] - M[1][0] * M[0][1];
1849 nummatrix<double, 2, 2> L{ {0.0, 0.0}, {0.0, 0.0} };
1850 if (fabs(detM) > 1e-2)
1851 L = nummatrix<double, 2, 2>{ { M[1][1], -M[0][1] }, { -M[1][0], M[0][0] } } *(1.0 / detM);
1852
1853 GIP[i] = (L & GIP[i]);
1854 DIFFV[i] = GIP[i] * (1.0 / OIP[i]);
1855 //GIP[i] = -GIP[i];
1856 }
1857
1858 //осредняем
1859 //*
1860//#pragma omp parallel for schedule(dynamic, 100)
1861 for (int i = 0; i < (int)vtx.size(); ++i)
1862 for (int j = 0; j < neib[i].size(); ++j)
1863 if ((vtx[i].r() - vtx[neib[i][j]].r()).length() <= 1.0 * eps[i])
1864 {
1865 sumDIFFV[i] += weight[i][j] * DIFFV[neib[i][j]];
1866 sumWeight[i] += weight[i][j];
1867 }
1868
1869 //#pragma omp parallel for schedule(dynamic, 100)
1870 for (int i = 0; i < (int)vtx.size(); ++i)
1871 {
1872 DIFFV[i] = sumDIFFV[i] * (1.0 / sumWeight[i]);
1873 }
1874 //*/
1875
1876//#ifdef TURB
1877//#pragma omp parallel for schedule(dynamic, 100)
1878// for (int i = 0; i < (int)vtx.size(); ++i)
1879// {
1880// Point2D myPos = vtx[i].r();
1881// Point2D resultGradNu = { 0.0, 0.0 };
1882// Point2D resultGradVx = { 0.0, 0.0 };
1883// Point2D resultGradVy = { 0.0, 0.0 };
1884//
1885// //if (i == 0)
1886// //{
1887// // std::cout << "h = " << h << std::endl;
1888// // std::cout << "VIP = " << VIP[0] << std::endl;
1889// // std::cout << "my_coord = " << vtx[i].r() << std::endl;
1890// // std::cout << "mtrC = " << mtrC[i] << std::endl;
1891// //}
1892//
1893// for (auto& nbr : neib[i])
1894// {
1895// //if (i == 0)
1896// // std::cout << "nb_" << nbr << "_coord = " << vtx[nbr].r() << std::endl;
1897//
1898// Point2D dr = myPos - vtx[nbr].r();
1899// double dist = dr.length();
1900// double gweight = gW(dist, h);
1901//
1902// if (nbr != i)
1903// {
1904// double nutNbr = (*nutPtr)[nbr];
1905// double nutI = (*nutPtr)[i];
1906// //double nutNbr = vtx[nbr].r() & Point2D { 1.0, 1.0 };
1907// //double nutSlf = vtx[i].r() & Point2D { 1.0, 1.0 };
1908//
1909// resultGradNu += (nutNbr - nutI) * (mtrL[i] & dr) * (gweight * /*vtx[nbr].g() / OIP[nbr]*/ VIP[nbr]);
1910// resultGradVx += (vel[nbr][0] - vel[i][0]) * (mtrL[i] & dr) * (gweight * /*vtx[nbr].g() / OIP[nbr]*/ VIP[nbr]);
1911// resultGradVy += (vel[nbr][1] - vel[i][1]) * (mtrL[i] & dr) * (gweight * /*vtx[nbr].g() / OIP[nbr]*/ VIP[nbr]);
1912// }
1913// }
1914//
1915// grNu[i] = resultGradNu;
1916// grV[i][0] = resultGradVx;
1917// grV[i][1] = resultGradVy;
1918//
1919// nummatrix<double, 2, 2> regGrV = grV[i];
1920// regGrV[0][0] -= grV[i][1][1];
1921// regGrV[0][1] += grV[i][1][0];
1922// regGrV[1][0] += grV[i][0][1];
1923// regGrV[1][1] -= grV[i][0][0];
1924//
1925// regGrV *= 0.5;
1926//
1927// regGrV[0][1] += 0.5 * OIP[i];
1928// regGrV[1][0] -= 0.5 * OIP[i];
1929//
1930// TP[i] = (regGrV & grNu[i]);
1931// }//for i
1932//#endif //TURB
1933
1934
1935//#ifdef TURB
1936//
1937//#pragma omp parallel for schedule(dynamic, 100)
1938// for (int i = 0; i < (int)vtx.size(); ++i)
1939// {
1940// Point2D myPos = vtx[i].r();
1941//
1942// Point2D resultGradTPx = { 0.0, 0.0 };
1943// Point2D resultGradTPy = { 0.0, 0.0 };
1944//
1945// for (auto& nbr : neib[i])
1946// {
1947// Point2D dr = myPos - vtx[nbr].r();
1948// double dist = dr.length();
1949// double gweight = gW(dist, h);
1950//
1951// if (nbr != i)
1952// {
1953// resultGradTPx += (TP[nbr][0] - TP[i][0]) * (mtrL[i] & dr) * (gweight * /*vtx[nbr].g() / OIP[nbr]*/ VIP[nbr]);
1954// resultGradTPy += (TP[nbr][1] - TP[i][1]) * (mtrL[i] & dr) * (gweight * /*vtx[nbr].g() / OIP[nbr]*/ VIP[nbr]);
1955// }
1956// }
1957//
1958// RIP[i] = resultGradTPy[0] - resultGradTPx[1];
1960// double dG = RIP[i] * VIP[i] * getPassport().timeDiscretizationProperties.dt;
1961// //if (i == 0)
1962// // std::cout << "nVirtVortex = " << nVirtVortex << std::endl;
1963// DELTAG[i] = 0.0;//dG;// (i < (wake->vtx.size() - nVirtVortex) ? dG : 0.0);
1964//
1965// //if (fabs(dG / getWake().vtx[i].g()) <= 0.25)
1966// getNonConstWake().vtx[i].g() += DELTAG[i];
1967// //else
1968// //{
1969// // if (dG * getWake().vtx[i].g() > 0)
1970// // getNonConstWake().vtx[i].g() *= xxx;
1971// // else
1972// // getNonConstWake().vtx[i].g() /= xxx;
1973// //}
1974//
1975// getNonConstWake().vtx[i].g() = std::clamp(getWake().vtx[i].g(), -1.5 * getPassport().wakeDiscretizationProperties.maxGamma, 1.5 * getPassport().wakeDiscretizationProperties.maxGamma);
1976// }//for i
1977//
1978// wake->SaveScalarFields("SPH", VIP, OIP, *nutPtr, grNu, RIP, DELTAG);
1979//#endif //TURB
1980
1981#endif //SPH
1982
1983
1984
1985 //newPos.clear();
1986 newPos.resize(nvt);
1987
1988
1989 // 1. Инициализация генератора случайных чисел (ПСЧГ)
1990 //std::random_device rd; // Источник энтропии для инициализации
1991 //std::mt19937 gen(rd()); // Mersenne Twister, инициализированный случайным значением
1992
1993 // 2. Определение параметров нормального распределения
1994 // μ (среднее) = 0.0, σ (стандартное отклонение) = 1.0
1995 //std::normal_distribution<> d(0.0, sqrt(2.0 * getPassport().physicalProperties.nu * getPassport().timeDiscretizationProperties.dt));
1996
1997
1998 //std::ofstream checkVel (passport.dir + "/dbg/checkVel" + std::to_string(getCurrentStep()) + ".txt");
1999 //std::ofstream checkVelDyn(passport.dir + "/dbg/checkVelDyn" + std::to_string(getCurrentStep()) + ".txt");
2000
2001 //checkVel << "i x y epsast Vx Vy Wx Wy corrWx corrWy Om gOmx gOmy\n";
2002 //checkVelDyn << "i x y epsast Vx Vy Wx Wy corrWx corrWy Om gOmx gOmy\n";
2003
2004//#pragma omp parallel for
2005 for (int i = 0; i < (int)wake->vtx.size(); ++i)
2006 {
2007#ifdef SPH
2008 //Point2D W = -(passport.physicalProperties.nu / OIP[i]) * GIP[i];
2009 Point2D W = -passport.physicalProperties.nu * DIFFV[i];
2010
2011 checkVel << i << " " << wake->vtx[i].r()[0] << " " << wake->vtx[i].r()[1] << " " << velocity->wakeVortexesParams.epsastWake[i] << " " \
2012 << velocity->wakeVortexesParams.convVelo[i][0] << " " << velocity->wakeVortexesParams.convVelo[i][1] << " " \
2013 << W[0] << " " << W[1];
2014
2015 if (W.length() > 1.5 * passport.physicalProperties.vRef)
2016 {
2017 std::cout << "i = " << i << ", W was " << W << std::endl;
2018 W.toZero();
2019 //W.normalize(1.5 * passport.physicalProperties.vRef);
2020 }
2021
2022 checkVel << " " << W[0] << " " << W[1] << " " << OIP[i] << " " << GIP[i][0] << " " << GIP[i][1] << "\n";
2023
2024#else
2025 Point2D W = velocity->wakeVortexesParams.diffVelo[i] * (nutPtr ? (1.0 + (*nutPtr)[i] / passport.physicalProperties.nu) : 1.0);
2026
2027 const auto& diff = velocity->wakeVortexesParams;
2028 Point2D grad = diff.I2[i] * (1.0 / (diff.epsastWake[i] * diff.I0[i])) - diff.I3[i] * (diff.I1[i] / sqr(diff.I0[i]));
2029
2030 //checkVelDyn << i << " " << wake->vtx[i].r()[0] << " " << wake->vtx[i].r()[1] << " " << velocity->wakeVortexesParams.epsastWake[i] << " " \
2031 // << velocity->wakeVortexesParams.convVelo[i][0] << " " << velocity->wakeVortexesParams.convVelo[i][1] << " " \
2032 // << W[0] << " " << W[1] << " " << W[0] << " " << W[1] << " " << diff.I1[i]/diff.I0[i] << " " << grad[0] << " " << grad[1] << "\n";
2033#endif
2034
2035 newPos[i] = wake->vtx[i].r() + (velocity->wakeVortexesParams.convVelo[i] + W + getV0()) * passport.timeDiscretizationProperties.dt;
2036
2037 //Chorin Random Walk method
2038 //newPos[i] = wake->vtx[i].r() + (velocity->wakeVortexesParams.convVelo[i] + getV0()) * passport.timeDiscretizationProperties.dt + Point2D{ d(gen), d(gen) };
2039 }
2040 //checkVel.close();
2041 //checkVelDyn.close();
2042
2043
2044 //for (int i = 0; i < (int)wake->vtx.size(); ++i)
2045 //{
2046 // std::cout << i << " " << wake->vtx[i].r() << " " << \
2047 // velocity->wakeVortexesParams.convVelo[i] << " " << \
2048 // velocity->wakeVortexesParams.diffVelo[i] << " " << \
2049 // getV0() << "\n";
2050 //}
2051
2052
2053
2054
2055 //getTimers().stop("Move");
2056}//MoveVortexes(...)
2057
2058//Метод - обертка для вызова метода генерации заголовка файла нагрузок и заголовка файла положения(последнее --- если профиль движется)
2059void World2D::GenerateMechanicsHeader(size_t mechanicsNumber)
2060{
2061 mechanics[mechanicsNumber]->GenerateForcesHeader();
2062 mechanics[mechanicsNumber]->GeneratePositionHeader();
2063}//GenerateMechanicsHeader(...)
2064
2065// Заполнение матрицы, состоящей из интегралов от (r-xi) / |r-xi|^2
2066void World2D::FillIQ()
2067{
2068 getTimers().start("MatRhs");
2069
2070 for (size_t bou = 0; bou < boundary.size(); ++bou)
2071 {
2072 if (currentStep == 0 || mechanics[bou]->isDeform)
2073 {
2074 boundary[bou]->FillIQSelf(IQ[bou][bou]);
2075 }
2076
2077 for (size_t oth = 0; oth < boundary.size(); ++oth)
2078 {
2079
2080#ifdef INITIAL
2081 if (currentStep == 0 || !useInverseMatrix)
2082 {
2083 //size_t nVarsOther = boundary[oth]->GetUnknownsSize();
2084 if (bou != oth)
2085 {
2086 //std::cout << "matr!" << std::endl;
2087 boundary[bou]->FillIQFromOther(*boundary[oth], IQ[bou][oth]);
2088 }
2089 }// if (currentStep == 0 || !useInverseMatrix)
2090#endif
2091
2092#ifdef BRIDGE
2093 if (currentStep == 0)
2094 {
2095 //size_t nVarsOther = boundary[oth]->GetUnknownsSize();
2096 if (bou != oth)
2097 boundary[bou]->FillIQFromOther(*boundary[oth], IQ[bou][oth]);
2098 }// if (currentStep == 0)
2099#endif
2100
2101 }// for oth
2102
2103 }// for bou
2104
2105 getTimers().stop("MatRhs");
2106}//FillIQ()
2107
2108void World2D::CalcAndSolveLinearSystem()
2109{
2110 if (airfoil.size() > 0)
2111 {
2112 ReserveMemoryForMatrixAndRhs();
2113
2114#if defined(__CUDACC__) || defined(USE_CUDA)
2115 cuda.setAccelCoeff(passport.physicalProperties.accelCft(getCurrentTime()));
2116 cuda.setMaxGamma(passport.wakeDiscretizationProperties.maxGamma);
2117
2118 int sch = passport.numericalSchemes.boundaryCondition.second;
2119
2120 if ((sch == 1) || (sch == 2) || (sch == 0))
2121 cuda.setSchemeSwitcher(sch);
2122 else
2123 {
2124 info('e') << "schemeSwitcher is not 0, or 1, or 2! " << std::endl;
2125 exit(1);
2126 }
2127
2128 cuda.RefreshWake(1);
2129 cuda.RefreshAfls(1);
2130 cuda.RefreshVirtualWakes(1);
2131
2132#endif
2133 const int linSystemScheme = passport.numericalSchemes.linearSystemSolver.second;
2134
2135 if (linSystemScheme == 0) //Gauss
2136 {
2137 if (currentStep == 0)
2138 {
2139#ifdef BRIDGE
2140 useInverseMatrix = true;
2141#endif //BRIDGE
2142
2143#ifdef INITIAL
2144 useInverseMatrix = (
2145 ((mechanics.size() == 1) && (!mechanics.front()->isDeform))
2146 ||
2147 (mechanics.size() > 1 && !isAnyMovableOrDeformable())
2148 );
2149#endif //INITIAL
2150 }
2151 }//if Gauss
2152
2153 if (linSystemScheme == 0 || linSystemScheme == 1 || (linSystemScheme == 2 && isAnyMovableOrDeformable()))
2154 FillIQ();
2155
2156 FillMatrixAndRhs();
2157
2158 //{
2159 // std::stringstream ss;
2160 // ss << "IQ1-" << currentStep;
2161 // std::ofstream of(passport.dir + "dbg/" + ss.str());
2162 // for (size_t i = 0; i < (IQ[0][0]).first.rows(); ++i)
2163 // {
2164 // for (size_t j = 0; j < (IQ[0][0]).first.cols(); ++j)
2165 // of << (IQ[0][0]).first(i, j) << " ";
2166 // of << std::endl;
2167 // }
2168 // of.close();
2169 //}
2170
2171 /*
2172 {
2173 std::stringstream ss;
2174 ss << "matr-" << currentStep;
2175 std::ofstream of(passport.dir + "dbg/" + ss.str());
2176 of.precision(16);
2177 for (size_t i = 0; i < matrReord.rows(); ++i)
2178 {
2179 for (size_t j = 0; j < matrReord.cols(); ++j)
2180 of << matrReord(i, j) << " ";
2181 of << std::endl;
2182 }
2183 of.close();
2184 }
2185
2186 {
2187 std::stringstream ss;
2188 ss << "rhs-" << currentStep;
2189 std::ofstream of(passport.dir + "dbg/" + ss.str());
2190 of.precision(16);
2191 for (size_t i = 0; i < rhsReord.rows(); ++i)
2192 {
2193 of << rhsReord(i) << std::endl;
2194 }
2195 of.close();
2196 }
2197
2198 //*/
2199
2200 SolveLinearSystem();
2201
2202 size_t currentRow = 0;
2203 for (size_t bou = 0; bou < boundary.size(); ++bou)
2204 {
2205 size_t nVars = boundary[bou]->GetUnknownsSize();
2206 Eigen::VectorXd locSol;
2207 locSol.resize(nVars);
2208 for (size_t i = 0; i < nVars; ++i)
2209 locSol(i) = sol(currentRow + i);
2210
2211 boundary[bou]->SolutionToFreeVortexSheetAndVirtualVortex(locSol);
2212 currentRow += nVars;// +1;
2213 }
2214
2215 //14-05-2024
2216 size_t nVirtVortices = 0;
2217 for (size_t bou = 0; bou < boundary.size(); ++bou)
2218 nVirtVortices += boundary[bou]->virtualWake.vtx.size();
2219
2220 wake->vtx.reserve(wake->vtx.size() + nVirtVortices);
2221 for (size_t bou = 0; bou < boundary.size(); ++bou)
2222 for (size_t v = 0; v < boundary[bou]->virtualWake.vtx.size(); ++v)
2223 wake->vtx.push_back(Vortex2D{ boundary[bou]->virtualWake.vtx[v].r(), 0.0, getPassport().wakeDiscretizationProperties.sigma0});
2224
2225 //for (size_t v = 0; v < wake->vtx.size(); ++v)
2226 //{
2227 // if (fabs(wake->vtx[v].g()) > 1.0)
2228 // std::cout << "Error gam! " << "v = " << v << ", gam[v] = " << wake->vtx[v].g() << std::endl;
2229 //}
2230
2231 //for (size_t bou = 0; bou < boundary.size(); ++bou)
2232 // for (size_t v = 0; v < airfoil[bou]->getNumberOfPanels(); ++v)
2233 // {
2234 // if (fabs(boundary[bou]->sheets.freeVortexSheet(v, 0)) > 1000.0)
2235 // std::cout << "Error gamma sheet! " << "bou = " << bou << ", v = " << v << ", gam[v] = " << boundary[bou]->sheets.freeVortexSheet(v, 0) << std::endl;
2236 // }
2237
2238 /*
2239 if (currentStep == 0)
2240 {
2241 Point2D addMass = { 0.0, 0.0 };
2242 for (size_t q = 0; q < airfoil[0]->getNumberOfPanels(); ++q)
2243 {
2244 addMass += (sol(q) + boundary[0]->sheets.attachedVortexSheet(q,0)) * 0.5 * (airfoil[0]->getR(q) + airfoil[0]->getR(q + 1)).kcross() * airfoil[0]->len[q];
2245 }
2246 addMass *= passport.physicalProperties.rho;
2247
2248 std::cout << "AddMass = " << addMass << std::endl;
2249 //exit(-42);
2250 }
2251 */
2252 }
2253}//CalcAndSolveLinearSystem()
2254
2255void World2D::WakeAndAirfoilsMotion(bool dynamics, std::vector<double>* nutPtr)
2256{
2257 std::vector<Point2D> newPos;
2258
2259 MoveVortexes(newPos, nutPtr);
2260
2261//#ifdef BRIDGE
2262// double totalForce = 0;
2263// for (size_t afl = 0; afl < airfoil.size(); ++afl)
2264// {
2265// totalForce += mechanics[afl]->hydroDynamForce[1];
2266// }
2267// totalForce *= 1.2;
2268//
2269// mechanics[0]->hydroDynamForce[1] = totalForce;
2270//#endif
2271
2272 oldAirfoil.resize(0);
2273 for (auto& afl : airfoil)
2274 {
2275 if (dynamic_cast<AirfoilRigid*>(afl.get()))
2276 oldAirfoil.emplace_back(new AirfoilGeometry(*afl));
2277
2278 if (dynamic_cast<AirfoilDeformable*>(afl.get()))
2279 oldAirfoil.emplace_back(new AirfoilGeometry(*afl));
2280
2281
2282#ifdef BRIDGE
2283 if (dynamics)
2284 {
2285 if (afl->numberInPassport == 0)
2286 {
2287 mechanics[afl->numberInPassport]->Move();
2288 }
2289 else
2290 {
2291 MechanicsRigidOscillPart* mechTest = dynamic_cast<MechanicsRigidOscillPart*>(mechanics[afl->numberInPassport].get());
2292 if (mechTest != nullptr)
2293 {
2294 Mechanics& mechGen = *mechanics[0];
2295 MechanicsRigidOscillPart& mech0 = dynamic_cast<MechanicsRigidOscillPart&>(mechGen);
2296
2297 Point2D dr = mech0.getR() - mech0.getROld();
2298 Point2D dv = mech0.getV() - mech0.getVOld();
2299
2300 double dphi = mech0.getPhi() - mech0.getPhiOld();
2301 double dw = mech0.getW() - mech0.getWOld();
2302
2303 Mechanics& mechGenI = *mechanics[afl->numberInPassport];
2304 MechanicsRigidOscillPart& mechI = dynamic_cast<MechanicsRigidOscillPart&>(mechGenI);
2305
2306 mechI.getROld() = mechI.getR();
2307 mechI.getVOld() = mechI.getV();
2308 mechI.getPhiOld() = mechI.getPhi();
2309 mechI.getWOld() = mechI.getW();
2310
2311 //std::cout << "afl = " << afl << ", dy = " << dy << std::endl;
2312
2313 airfoil[afl->numberInPassport]->Move(dr);
2314 airfoil[afl->numberInPassport]->Rotate(dphi);
2315
2316 mechI.getR() += dr;
2317 mechI.getV() += dv;
2318
2319 mechI.getPhi() += dphi;
2320 mechI.getW() += dw;
2321 }
2322 }
2323 }
2324#endif
2325
2326#ifdef INITIAL
2327 if (dynamics)
2328 mechanics[afl->numberInPassport]->Move();
2329 else
2330 {
2331 MechanicsRigidOscillPart* mech = dynamic_cast<MechanicsRigidOscillPart*>(mechanics[afl->numberInPassport].get());
2332 if (mech)
2333 mech->MoveKinematic();
2334 else
2335 exit(2222);
2336 }
2337#endif
2338 }//for
2339
2340#if defined(__CUDACC__) || defined(USE_CUDA)
2341 cuda.RefreshAfls(2);
2342#endif
2343
2344 for (auto& bou : boundary)
2345 bou->virtualWake.vtx.clear();
2346
2347
2348 CheckInside(newPos, oldAirfoil);
2349
2350
2351 //передача новых положений вихрей в пелену
2352 for (size_t i = 0; i < wake->vtx.size(); ++i)
2353 {
2354 wake->vtx[i].r() = newPos[i];
2355 //wake->vtx[i].sigma() = getPassport().wakeDiscretizationProperties.sigma0;
2356 }
2357
2358// getWake().SaveKadrVtk();
2359
2360}//WakeAndAirfoilsMotion()
2361
2362// Возврат признака того, что хотя бы один из профилей подвижный
2363bool World2D::isAnyMovable() const
2364{
2365 return std::any_of(getMechanicsVector().begin(), getMechanicsVector().end(),
2366 [](const std::unique_ptr<Mechanics>& m) {return (m->isMoves); });
2367}
2368
2369// Возврат признака того, что хотя бы один из профилей подвижный или деформируемый
2370bool World2D::isAnyMovableOrDeformable() const
2371{
2372 return std::any_of(getMechanicsVector().begin(), getMechanicsVector().end(),
2373 [](const std::unique_ptr<Mechanics>& m) {return (m->isMoves || m->isDeform); });
2374}
Airfoil2DDeformable.h
Заголовочный файл с описанием класса AirfoilDeformable.
Airfoil2DRigid.h
Заголовочный файл с описанием класса AirfoilRect.
Boundary2DConstLayerAver.h
Заголовочный файл с описанием класса BoundaryConstLayerAver.
Boundary2DLinLayerAver.h
Заголовочный файл с описанием класса BoundaryLinLayerAver.
Boundary2DVortexCollocN.h
Заголовочный файл с описанием класса BoundaryVortexCollocN.
Gmres2D.h
Заголовочный файл с функциями для метода GMRES.
MeasureVP2D.h
Заголовочный файл с описанием класса MeasureVP.
Mechanics2DDeformable.h
Заголовочный файл с описанием класса MechanicsDeformable.
Mechanics2DRigidGivenLaw.h
Заголовочный файл с описанием класса MechanicsRigidGivenLaw.
Mechanics2DRigidImmovable.h
Заголовочный файл с описанием класса MechanicsRigidImmovable.
Mechanics2DRigidOscillPart.h
Заголовочный файл с описанием класса MechanicsRigidOscillPart.
Mechanics2DRigidRotatePart.h
Заголовочный файл с описанием класса MechanicsRigidRotatePart.
StreamParser.h
Заголовочный файл с описанием класса StreamParser.
PI
const double PI
Число .
Definition defs.h:76
Velocity2DBarnesHut.h
Заголовочный файл с описанием класса VelocityBarnesHut.
Velocity2DBiotSavart.h
Заголовочный файл с описанием класса VelocityBiotSavart.
Wake2D.h
Заголовочный файл с описанием класса Wake.
World2D.h
Заголовочный файл с описанием класса World2D.
VM2D::AirfoilDeformable
Класс, определяющий тип обтекаемого профиля
Definition Airfoil2DDeformable.h:66
VM2D::AirfoilGeometry
Класс, определяющий форму профиля
Definition Airfoil2D.h:66
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::nrm
std::vector< Point2D > nrm
Нормали к панелям профиля
Definition Airfoil2D.h:81
VM2D::AirfoilGeometry::rcm
Point2D rcm
Положение центра масс профиля
Definition Airfoil2D.h:97
VM2D::AirfoilGeometry::getNumberOfPanels
size_t getNumberOfPanels() const
Возврат количества панелей на профиле
Definition Airfoil2D.h:163
VM2D::Airfoil::calcMeanEpsOverPanel
void calcMeanEpsOverPanel()
Вычисление средних значений eps на панелях
Definition Airfoil2D.cpp:96
VM2D::AirfoilRigid
Класс, определяющий тип обтекаемого профиля
Definition Airfoil2DRigid.h:66
VM2D::BoundaryConstLayerAver
Класс, определяющий способ удовлетворения граничного условия на обтекаемом профиле
Definition Boundary2DConstLayerAver.h:66
VM2D::Boundary::sheets
Sheet sheets
Слои на профиле
Definition Boundary2D.h:96
VM2D::BoundaryLinLayerAver
Класс, определяющий способ удовлетворения граничного условия на обтекаемом профиле
Definition Boundary2DLinLayerAver.h:66
VM2D::BoundaryVortexCollocN
Класс, определяющий способ удовлетворения граничного условия на обтекаемом профиле
Definition Boundary2DVortexCollocN.h:68
VM2D::CpuTreeInfo
Структура, хранящая данные и указатели на массивы на CPU для оптимизации итерационного решения СЛАУ н...
Definition cpuTreeInfo.h:96
VM2D::Gpu
Класс, обеспечивающий возможность выполнения вычислений на GPU по технологии Nvidia CUDA.
Definition Gpu2D.h:71
VM2D::Gpu::setAccelCoeff
void setAccelCoeff(double cft_)
Установка коэффициента разгона потока
Definition Gpu2D.h:281
VM2D::Gpu::setMaxGamma
void setMaxGamma(double gam_)
Установка максимально допустимой циркуляции вихря
Definition Gpu2D.h:304
VM2D::Gpu::setSchemeSwitcher
void setSchemeSwitcher(int schemeSwitcher_)
Установка переключателя расчетных схем
Definition Gpu2D.h:318
VM2D::MeasureVP
Класс, отвечающий за вычисление поля скорости и давления в заданых точках для вывода
Definition MeasureVP2D.h:65
VM2D::MechanicsDeformable
Класс, определяющий вид механической системы
Definition Mechanics2DDeformable.h:131
VM2D::MechanicsDeformable::initialChord
std::vector< ChordPanel > initialChord
Definition Mechanics2DDeformable.h:179
VM2D::MechanicsDeformable::chord
std::vector< ChordPanel > chord
Definition Mechanics2DDeformable.h:180
VM2D::MechanicsDeformable::beam
std::unique_ptr< Beam > beam
Definition Mechanics2DDeformable.h:184
VM2D::Mechanics
Абстрактный класс, определяющий вид механической системы
Definition Mechanics2D.h:72
VM2D::Mechanics::hydroDynamForce
Point2D hydroDynamForce
Вектор гидродинамической силы и момент, действующие на профиль
Definition Mechanics2D.h:134
VM2D::Mechanics::GeneratePositionString
void GeneratePositionString()
Сохранение строки со статистикой в файл нагрузок
Definition Mechanics2D.cpp:239
VM2D::Mechanics::GenerateForcesString
void GenerateForcesString()
Сохранение строки со статистикой в файл нагрузок
Definition Mechanics2D.cpp:184
VM2D::MechanicsRigidGivenLaw
Класс, определяющий вид механической системы
Definition Mechanics2DRigidGivenLaw.h:69
VM2D::MechanicsRigidImmovable
Класс, определяющий вид механической системы
Definition Mechanics2DRigidImmovable.h:66
VM2D::MechanicsRigidOscillPart
Класс, определяющий вид механической системы
Definition Mechanics2DRigidOscillPart.h:66
VM2D::MechanicsRigidOscillPart::getR
Point2D & getR()
текущее отклонение профиля
Definition Mechanics2DRigidOscillPart.h:100
VM2D::MechanicsRigidOscillPart::getPhi
double & getPhi()
текущий угол поворота профиля
Definition Mechanics2DRigidOscillPart.h:108
VM2D::MechanicsRigidOscillPart::getV
Point2D & getV()
текущая скорость профиля
Definition Mechanics2DRigidOscillPart.h:96
VM2D::MechanicsRigidOscillPart::getW
double & getW()
текущая угловая скорость профиля
Definition Mechanics2DRigidOscillPart.h:104
VM2D::MechanicsRigidRotatePart
Класс, определяющий вид механической системы
Definition Mechanics2DRigidRotatePart.h:67
VM2D::Passport
Класс, опеделяющий паспорт двумерной задачи
Definition Passport2D.h:256
VM2D::Passport::wakesDir
std::string wakesDir
Каталог с файлами вихревых следов
Definition Passport2D.h:276
VM2D::Passport::physicalProperties
PhysicalProperties physicalProperties
Структура с физическими свойствами задачи
Definition Passport2D.h:304
VM2D::Passport::defaultsFileFullName
std::string defaultsFileFullName
Definition Passport2D.h:284
VM2D::Passport::wakeDiscretizationProperties
WakeDiscretizationProperties wakeDiscretizationProperties
Структура с параметрами дискретизации вихревого следа
Definition Passport2D.h:307
VM2D::Passport::varLine
std::vector< std::string > varLine
Definition Passport2D.h:286
VM2D::Passport::airfoilsDir
std::string airfoilsDir
Каталог с файлами профилей
Definition Passport2D.h:273
VM2D::Passport::airfoilParams
std::vector< AirfoilParams > airfoilParams
Список структур с параметрами профилей
Definition Passport2D.h:279
VM2D::Passport::switchersFileFullName
std::string switchersFileFullName
Definition Passport2D.h:285
VM2D::Passport::mechanicsFileFullName
std::string mechanicsFileFullName
Definition Passport2D.h:283
VM2D::Passport::fileFullName
std::string fileFullName
Имена файлов
Definition Passport2D.h:282
VM2D::Passport::GetReviseParamsFromParser
void GetReviseParamsFromParser(const Passport &newPassport, const std::vector< std::string > paramList)
Считывание измененных параметров
Definition Passport2D.cpp:480
VM2D::Passport::numericalSchemes
NumericalSchemes numericalSchemes
Структура с используемыми численными схемами
Definition Passport2D.h:310
VM2D::Sheet::attachedVortexSheet
const double & attachedVortexSheet(size_t n, size_t moment) const
Definition Sheet2D.h:105
VM2D::Sheet::freeVortexSheet
const double & freeVortexSheet(size_t n, size_t moment) const
Definition Sheet2D.h:100
VM2D::VelocityBarnesHut
Класс, определяющий способ вычисления скоростей
Definition Velocity2DBarnesHut.h:65
VM2D::VelocityBiotSavart
Класс, определяющий способ вычисления скоростей
Definition Velocity2DBiotSavart.h:65
VM2D::WakeDataBase
Класс, опеделяющий набор вихрей
Definition WakeDataBase2D.h:67
VM2D::WakeDataBase::vtx
std::vector< Vortex2D > vtx
Список вихревых элементов
Definition WakeDataBase2D.h:81
VM2D::Wake
Класс, опеделяющий вихревой след (пелену)
Definition Wake2D.h:63
VM2D::World2D::getNumberOfAirfoil
size_t getNumberOfAirfoil() const
Возврат количества профилей в задаче
Definition World2D.h:180
VM2D::World2D::inflTreeWake
std::unique_ptr< CpuTreeInfo > inflTreeWake
Деревья для быстрого метода
Definition World2D.h:113
VM2D::World2D::cuda
Gpu cuda
Объект, управляющий графическим ускорителем
Definition World2D.h:143
VM2D::World2D::ReserveMemoryForMatrixAndRhs
void ReserveMemoryForMatrixAndRhs()
Вычисляем размер матрицы и резервируем память под нее и под правую часть
Definition World2D.cpp:1396
VM2D::World2D::isAnyMovable
bool isAnyMovable() const
Возврат признака того, что хотя бы один из профилей подвижный
Definition World2D.cpp:2363
VM2D::World2D::WakeAndAirfoilsMotion
void WakeAndAirfoilsMotion(bool dynamics, std::vector< double > *nutPtr=nullptr)
Перемещение вихрей и профилей на шаге
Definition World2D.cpp:2255
VM2D::World2D::airfoil
std::vector< std::unique_ptr< Airfoil > > airfoil
Список умных указателей на обтекаемые профили
Definition World2D.h:80
VM2D::World2D::timerInitialBuild
VMlib::vmTimer timerInitialBuild
Definition World2D.h:366
VM2D::World2D::isAnyMovableOrDeformable
bool isAnyMovableOrDeformable() const
Возврат признака того, что хотя бы один из профилей подвижный или деформируемый
Definition World2D.cpp:2370
VM2D::World2D::measureVP
std::unique_ptr< MeasureVP > measureVP
Умный указатель на алгоритм вычисления полей скоростей и давления (для сохранения в файл)
Definition World2D.h:104
VM2D::World2D::matrReord
Eigen::MatrixXd matrReord
Матрица системы
Definition World2D.h:122
VM2D::World2D::useInverseMatrix
bool useInverseMatrix
Признак использования обратной матрицы
Definition World2D.h:131
VM2D::World2D::getWake
const Wake & getWake() const
Возврат константной ссылки на вихревой след
Definition World2D.h:232
VM2D::World2D::invMatr
Eigen::MatrixXd invMatr
Обратная матрица
Definition World2D.h:128
VM2D::World2D::CalcPanelsVeloAndAttachedSheets
void CalcPanelsVeloAndAttachedSheets()
Вычисление скоростей панелей и интенсивностей присоединенных слоев вихрей и источников
Definition World2D.cpp:1637
VM2D::World2D::CalcVortexVelo
void CalcVortexVelo()
Вычисление скоростей (и конвективных, и диффузионных) вихрей (в пелене и виртуальных),...
Definition World2D.cpp:1456
VM2D::World2D::oldAirfoil
std::vector< std::unique_ptr< AirfoilGeometry > > oldAirfoil
Список умных указателей на обтекаемые профили для сохранения старого положения
Definition World2D.h:83
VM2D::World2D::World2D
World2D(const VMlib::PassportGen &passport_)
Конструктор
Definition World2D.cpp:73
VM2D::World2D::getAirfoil
const Airfoil & getAirfoil(size_t i) const
Возврат константной ссылки на объект профиля
Definition World2D.h:163
VM2D::World2D::getMechanicsVector
const std::vector< std::unique_ptr< Mechanics > > & getMechanicsVector() const
Definition World2D.h:227
VM2D::World2D::timerInside
VMlib::vmTimer timerInside
Definition World2D.h:372
VM2D::World2D::getNonConstCuda
Gpu & getNonConstCuda() const
Возврат неконстантной ссылки на объект, связанный с видеокартой (GPU)
Definition World2D.h:278
VM2D::World2D::source
std::unique_ptr< WakeDataBase > source
Умный указатель на источники
Definition World2D.h:101
VM2D::World2D::CalcAndSolveLinearSystem
void CalcAndSolveLinearSystem()
Набор матрицы, правой части и решение СЛАУ
Definition World2D.cpp:2108
VM2D::World2D::rhsReord
Eigen::VectorXd rhsReord
Правая часть системы
Definition World2D.h:134
VM2D::World2D::MoveVortexes
void MoveVortexes(std::vector< Point2D > &newPos, std::vector< double > *nutPtr=nullptr)
Вычисляем новые положения вихрей (в пелене и виртуальных)
Definition World2D.cpp:1655
VM2D::World2D::dispBoundaryInSystem
std::vector< size_t > dispBoundaryInSystem
Список номеров, с которых начинаются элементы правой части (или матрицы) системы для профилей
Definition World2D.h:89
VM2D::World2D::getTimers
VMlib::TimersGen & getTimers() const
Возврат ссылки на временную статистику выполнения шага расчета по времени
Definition World2D.h:288
VM2D::World2D::FillIQ
void FillIQ()
Заполнение матрицы, состоящей из интегралов от (r-xi) / |r-xi|^2.
Definition World2D.cpp:2066
VM2D::World2D::mechanics
std::vector< std::unique_ptr< Mechanics > > mechanics
Список умных указателей на типы механической системы для каждого профиля
Definition World2D.h:92
VM2D::World2D::getCuda
const Gpu & getCuda() const
Возврат константной ссылки на объект, связанный с видеокартой (GPU)
Definition World2D.h:273
VM2D::World2D::boundary
std::vector< std::unique_ptr< Boundary > > boundary
Список умных указателей на формирователи граничных условий на профилях
Definition World2D.h:86
VM2D::World2D::getPassport
const Passport & getPassport() const
Возврат константной ссылки на паспорт
Definition World2D.h:263
VM2D::World2D::cntrTreePnl
std::unique_ptr< CpuTreeInfo > cntrTreePnl
Definition World2D.h:116
VM2D::World2D::getNonConstPassport
Passport & getNonConstPassport() const
Возврат неконстантной ссылки на паспорт
Definition World2D.h:268
VM2D::World2D::GenerateMechanicsHeader
void GenerateMechanicsHeader(size_t mechanicsNumber)
Definition World2D.cpp:2059
VM2D::World2D::cntrTreeVP
std::unique_ptr< CpuTreeInfo > cntrTreeVP
Definition World2D.h:115
VM2D::World2D::CheckInside
void CheckInside(std::vector< Point2D > &newPos, const std::vector< std::unique_ptr< AirfoilGeometry > > &oldAirfoil)
Проверка проникновения вихрей внутрь профиля
Definition World2D.cpp:813
VM2D::World2D::SolveLinearSystem
void SolveLinearSystem()
Решение системы линейных алгебраических уравнений
Definition World2D.cpp:1001
VM2D::World2D::wake
std::unique_ptr< Wake > wake
Умный указатель на вихревой след
Definition World2D.h:98
VM2D::World2D::getBoundary
const Boundary & getBoundary(size_t i) const
Возврат константной ссылки на объект граничного условия
Definition World2D.h:186
VM2D::World2D::ifDivisible
bool ifDivisible(int val) const
Definition World2D.h:290
VM2D::World2D::getNumberOfBoundary
size_t getNumberOfBoundary() const
Возврат количества граничных условий в задаче
Definition World2D.h:197
VM2D::World2D::velocity
std::unique_ptr< Velocity > velocity
Умный укзатель на объект, определяющий методику вычисления скоростей
Definition World2D.h:95
VM2D::World2D::passport
const Passport & passport
Константная ссылка на паспорт конкретного расчета
Definition World2D.h:140
VM2D::World2D::FillMatrixAndRhs
void FillMatrixAndRhs()
Заполнение матрицы системы для всех профилей
Definition World2D.cpp:1298
VM2D::World2D::Step
virtual void Step() override
Функция выполнения предварительного шага
Definition World2D.cpp:279
VM2D::World2D::timerFillMatrix
VMlib::vmTimer timerFillMatrix
Definition World2D.h:368
VM2D::World2D::sol
Eigen::VectorXd sol
Решение системы
Definition World2D.h:137
VM2D::World2D::getNonConstMechanics
Mechanics & getNonConstMechanics(size_t i) const
Возврат неконстантной ссылки на объект механики
Definition World2D.h:225
VM2D::World2D::IQ
std::vector< std::vector< std::pair< Eigen::MatrixXd, Eigen::MatrixXd > > > IQ
Матрица, состоящая из пар матриц, в которых хранятся касательные и нормальные компоненты интегралов о...
Definition World2D.h:125
VM2D::World2D::getNonConstAirfoil
Airfoil & getNonConstAirfoil(size_t i) const
Возврат неконстантной ссылки на объект профиля
Definition World2D.h:175
VM2D::World2D::timerSlaeSolve
VMlib::vmTimer timerSlaeSolve
Definition World2D.h:369
VM2D::World2D::cntrTreeWake
std::unique_ptr< CpuTreeInfo > cntrTreeWake
Definition World2D.h:114
VM2D::World2D::getNonConstWake
Wake & getNonConstWake() const
Возврат неконстантной ссылки на вихревой след
Definition World2D.h:237
VMlib::LogStream::endl
void endl()
Вывод в поток логов пустой строки
Definition LogStream.h:103
VMlib::LogStream::assignStream
void assignStream(std::ostream *pStr_, const std::string &label_)
Связывание потока логов с потоком вывода
Definition LogStream.h:80
VMlib::PassportGen
Абстрактный класс, опеделяющий паспорт задачи
Definition PassportGen.h:102
VMlib::PassportGen::timeDiscretizationProperties
TimeDiscretizationProperties timeDiscretizationProperties
Структура с параметрами процесса интегрирования по времени
Definition PassportGen.h:138
VMlib::PassportGen::problemName
std::string problemName
Название задачи
Definition PassportGen.h:132
VMlib::PassportGen::dir
std::string dir
Рабочий каталог задачи
Definition PassportGen.h:129
VMlib::PassportGen::problemNumber
size_t problemNumber
Номер задачи
Definition PassportGen.h:135
VMlib::TimersGen::GenerateStatString
void GenerateStatString(size_t stepNo, double curTime, size_t N)
Формирование очередной строки файла временной статистики
Definition TimesGen.cpp:111
VMlib::TimersGen::stop
void stop(const std::string &timerLabel)
Останов счетчика
Definition TimesGen.cpp:68
VMlib::TimersGen::start
void start(const std::string &timerLabel)
Запуск счетчика
Definition TimesGen.cpp:55
VMlib::TimersGen::resetAll
void resetAll()
Сброс всех счетчиков
Definition TimesGen.cpp:81
VMlib::TimersGen::durationStep
double durationStep() const
Вывод счетчика всего шага в секундах
Definition TimesGen.h:184
VMlib::Vortex2D
Класс, опеделяющий двумерный вихревой элемент
Definition Vortex2D.h:59
VMlib::WorldGen::timers
std::unique_ptr< TimersGen > timers
Сведения о временах выполнения основных операций
Definition WorldGen.h:66
VMlib::WorldGen::currentTime
double currentTime
Текущее время в решаемой задаче
Definition WorldGen.h:72
VMlib::WorldGen::getInfo
VMlib::LogStream & getInfo() const
Возврат ссылки на объект LogStream Используется в техничеcких целях для организации вывода
Definition WorldGen.h:82
VMlib::WorldGen::getCurrentTime
double getCurrentTime() const
Definition WorldGen.h:100
VMlib::WorldGen::nVtxBeforeMerging
size_t nVtxBeforeMerging
Definition WorldGen.h:76
VMlib::WorldGen::currentStep
size_t currentStep
Текущий номер шага в решаемой задаче
Definition WorldGen.h:69
VMlib::WorldGen::getCurrentStep
size_t getCurrentStep() const
Возврат константной ссылки на параметры распараллеливания по MPI.
Definition WorldGen.h:99
VMlib::WorldGen::info
LogStream info
Поток для вывода логов и сообщений об ошибках
Definition WorldGen.h:60
VMlib::nummatrix
Шаблонный класс, определяющий матрицу фиксированного размера Фактически представляет собой массив,...
Definition nummatrix.h:66
VMlib::nummatrix::toZero
nummatrix< T, n, m > & toZero(T val=0.0)
Definition nummatrix.h:300
VMlib::numvector::kcross
numvector< T, 2 > kcross() const
Геометрический поворот двумерного вектора на 90 градусов
Definition numvector.h:511
VMlib::numvector::toZero
numvector< T, n > & toZero(P val=0)
Установка всех компонент вектора в константу (по умолчанию — нуль)
Definition numvector.h:528
VMlib::numvector::length
P length() const
Вычисление 2-нормы (длины) вектора
Definition numvector.h:374
VMlib::vmTimer::stop
const vmTimer & stop() const
Останов работающего счетчика времени
Definition TimesGen.h:125
VMlib::vmTimer::start
const vmTimer & start() const
Запуск (первый или повторный) счетчика времени
Definition TimesGen.h:109
VMlib::vmTimer::reset
const vmTimer & reset() const
Сброс счетчика времени
Definition TimesGen.h:101
VMlib::fileExistTest
bool fileExistTest(const std::string &fileName_, LogStream &info, bool exitKey=false, const std::list< std::string > &extList={})
Проверка существования файла
Definition defs.h:346
VMlib::PrintLogoToStream
void PrintLogoToStream(std::ostream &str)
Передача в поток вывода шапки программы VM2D/VM3D.
Definition defs.cpp:108
VMlib::sqr
T sqr(T x)
Возведение числа в квадрат
Definition defs.h:463
defaults::defaultWorld2DLogStream
static std::ostream * defaultWorld2DLogStream
Поток вывода логов и ошибок задачи
Definition defs.h:228
nummatrix.h
Описание класса nummatrix.
VM2D::NumericalSchemes::boundaryCondition
std::pair< std::string, int > boundaryCondition
Метод аппроксимации граничных условий
Definition Passport2D.h:193
VM2D::NumericalSchemes::velocityComputation
std::pair< std::string, int > velocityComputation
Definition Passport2D.h:183
VM2D::NumericalSchemes::linearSystemSolver
std::pair< std::string, int > linearSystemSolver
Definition Passport2D.h:177
VM2D::PhysicalProperties::accelCft
double accelCft(double currentTime) const
Функция-множитель, позволяющая моделировать разгон
Definition Passport2D.cpp:50
VM2D::PhysicalProperties::vRef
double vRef
Референсная скорость
Definition Passport2D.h:81
VM2D::PhysicalProperties::nu
double nu
Коэффициент кинематической вязкости среды
Definition Passport2D.h:99
VM2D::PhysicalProperties::rho
double rho
Плотность потока
Definition Passport2D.h:75
VM2D::WakeDiscretizationProperties::fileSource
std::string fileSource
Имя файла с положениями источников (без полного пути)
Definition Passport2D.h:150
VM2D::WakeDiscretizationProperties::epscol
double epscol
Радиус коллапса
Definition Passport2D.h:132
VM2D::WakeDiscretizationProperties::fileWake
std::string fileWake
Имя файла с начальным состоянием вихревого следа (без полного пути)
Definition Passport2D.h:147
VM2D::WakeDiscretizationProperties::sigma0
double sigma0
Радиус вихря
Definition Passport2D.h:126
VM2D::WakeDiscretizationProperties::maxGamma
double maxGamma
Максимально допустимая циркуляция вихря
Definition Passport2D.h:144
VMlib::TimeDiscretizationProperties::dt
double dt
Шаг по времени
Definition PassportGen.h:67
VMlib::TimeDiscretizationProperties::timeStart
double timeStart
Начальное время
Definition PassportGen.h:61
VMlib::TimeDiscretizationProperties::reviseParameters
std::vector< std::string > reviseParameters
Список перечитываемых параметров
Definition PassportGen.h:73
VMlib::TimeDiscretizationProperties::saveVPstep
int saveVPstep
Шаг вычисления и сохранения скорости и давления
Definition PassportGen.h:86
VMlib::TimeDiscretizationProperties::revisePassportStep
int revisePassportStep
Шаг перечитывания паспорта
Definition PassportGen.h:70
VMlib::TimeDiscretizationProperties::timeStop
double timeStop
Конечное время
Definition PassportGen.h:64