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