VM2D 1.14
Vortex methods for 2D flows simulation
Loading...
Searching...
No Matches
Wake2D.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: Wake2D.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#if defined(_WIN32)
41 #include <direct.h>
42#endif
43
44#include "Wake2D.h"
45
46#include "knnCPU.h"
47#include "cuKnn.cuh"
48
49#include "treeKernels.cuh"
50
51#include "Airfoil2D.h"
52#include "Boundary2D.h"
53#include "MeasureVP2D.h"
54#include "Mechanics2D.h"
55#include "Preprocessor.h"
56#include "StreamParser.h"
57#include "Velocity2D.h"
58#include "World2D.h"
59
60using namespace VM2D;
61
62bool Wake::MoveInside(const Point2D& newPos, const Point2D& oldPos, const Airfoil& afl, size_t& panThrough) const
63{
64 //const double porog_r = 1e-12;
65
66 //double minDist = 1.0E+10; //расстояние до пробиваемой панели
67 panThrough = -1;
68
69 //проверка габ. прямоугольника
70
71 if (afl.isOutsideGabarits(newPos) && afl.isOutsideGabarits(oldPos))
72 {
73 ++W.gabb;
74 return false;
75 }
76
77 //если внутри габ. прямоугольника - проводим контроль
78 bool hit = false;
79
80
81 for (size_t j = 0; j < afl.getNumberOfPanels(); ++j)
82 {
83 const Point2D& aflRj = afl.getR(j);
84 const Point2D& aflRj1 = afl.getR(j + 1);
85
86 ++W.checkPan;
87
88 if ((((aflRj - oldPos) ^ (newPos - oldPos)) * ((aflRj1 - oldPos) ^ (newPos - oldPos)) <= 0) && \
89 (((oldPos - aflRj) ^ (aflRj1 - aflRj)) * ((newPos - aflRj) ^ (aflRj1 - aflRj)) <= 0))
90 {
91 hit = true;
92 panThrough = j;
93 break;
94 }
95 }//for j
96
97
98
99 if (hit)
100 ++W.check01;
101 else
102 ++W.check02;
103
104 return hit;
105}//MoveInside(...)
106
107
108//bool Wake::MoveInside(const Point2D& newPos, const Point2D& oldPos, const Airfoil& afl, size_t& panThrough)
109//{
110// const double porog_r = 1e-12;
111//
112// double minDist = 1.0E+10; //расстояние до пробиваемой панели
113// panThrough = -1;
114//
115// //проверка габ. прямоугольника
116// if (afl.isOutsideGabarits(newPos) && afl.isOutsideGabarits(oldPos))
117// return false;
118//
119// //если внутри габ. прямоугольника - проводим контроль
120// bool hit = false;
121//
122// //Определение прямой: Ax+By+D=0 - перемещение вихря
123// double A = newPos[1] - oldPos[1];
124// double B = oldPos[0] - newPos[0];
125// double D = oldPos[1] * newPos[0] - oldPos[0] * newPos[1];
126//
127// double A1, B1, D1;
128//
129// std::pair<double, double> gabPointX = std::minmax(newPos[0], oldPos[0]);
130// std::pair<double, double> gabPointY = std::minmax(newPos[1], oldPos[1]);
131//
132// std::pair<double, double> gabPanelX, gabPanelY;
133//
134// //Проверка на пересечение
135// double r0 = 0, r1 = 0, r2 = 0, r3 = 0;
136//
137// for (size_t j = 0; j < afl.getNumberOfPanels(); ++j)
138// {
139// r0 = A * afl.getR(j)[0] + B * afl.getR(j)[1] + D;
140// r1 = A * afl.getR(j + 1)[0] + B * afl.getR(j + 1)[1] + D;
141//
142// if (fabs(r0) < porog_r) r0 = 0.0;
143// if (fabs(r1) < porog_r) r1 = 0.0;
144//
145// if (r0*r1 > 0)
146// continue;
147//
148// //Определение прямой:A1x+B1y+D1=0 - панель
149// A1 = afl.getR(j + 1)[1] - afl.getR(j)[1];
150// B1 = afl.getR(j)[0] - afl.getR(j + 1)[0];
151// D1 = afl.getR(j)[1] * afl.getR(j + 1)[0] - afl.getR(j)[0] * afl.getR(j + 1)[1];
152//
153// r2 = A1 * oldPos[0] + B1 * oldPos[1] + D1;
154// r3 = A1 * newPos[0] + B1 * newPos[1] + D1;
155//
156// if (fabs(r2) < porog_r) r2 = 0.0;
157// if (fabs(r3) < porog_r) r3 = 0.0;
158//
159// if (r2*r3 > 0)
160// continue;
161//
162// hit = true;// пробила!
163// double d2 = (oldPos[0] - (B*D1 - D * B1) / (A*B1 - B * A1))*(oldPos[0] - (B*D1 - D * B1) / (A*B1 - B * A1)) + \
164// (oldPos[1] - (A1*D - D1 * A) / (A*B1 - B * A1))*(oldPos[1] - (A1*D - D1 * A) / (A*B1 - B * A1));
165//
166// if (d2 < minDist)
167// {
168// minDist = d2;
169// panThrough = j;
170// }//if d2
171// }//for j
172//
173//
174// return hit;
175//}//MoveInside(...)
176
177
178bool Wake::MoveInsideMovingBoundary(const Point2D& newPos, const Point2D& oldPos, const AirfoilGeometry& oldAfl, const Airfoil& afl, size_t& panThrough) const
179{
180 panThrough = -1;
181
183
184 //проверка габ. прямоугольника
185 if (!afl.inverse && afl.isOutsideGabarits(newPos) && afl.isOutsideGabarits(oldPos))
186 return false;
187
188 bool hit = false;
189
190 double angle = 0;
191 double cs, sn;
192
193 double dist2 = 1000000000.0;
194
195 for (size_t i = 0; i < afl.getNumberOfPanels(); ++i)
196 {
197 Point2D v1, v2, vv;
198 v1 = afl.getR(i) - newPos;
199
200 v2 = afl.getR(i + 1) - newPos;
201
202 vv = afl.getR(i + 1) - afl.getR(i);
203
204 double dst = v1.length2() + v2.length2();
205 if (dst < dist2)
206 {
207 dist2 = dst;
208 panThrough = i;
209 }
210
211 cs = v1 & v2;
212 sn = v1 ^ v2;
213
214 angle += atan2(sn, cs);
215 }//for i
216
217 hit = ((angle > 3.14) || (angle < -3.14));
218
219 if (afl.inverse)
220 hit = !hit;
221
222 return hit;
223}//MoveInsideMovingBoundary(...)
224
225
226//bool Wake::MoveInsideMovingBoundary(const Point2D& newPos, const Point2D& oldPos, const Airfoil& oldAfl, const Airfoil& afl, size_t& panThrough)
227//{
228// panThrough = -1;
229//
230// /// \todo сравнить производительности двух inside-ов
231//
232// //проверка габ. прямоугольника
233// if (afl.isOutsideGabarits(newPos) && afl.isOutsideGabarits(oldPos))
234// return false;
235//
236// bool hit = false;
237//
238// double angle = 0;
239// double cs, sn;
240//
241// double dist2 = 1000000000.0;
242//
243// for (int i = 0; i < afl.getNumberOfPanels(); ++i)
244// {
245// Point2D v1, v2, vv;
246// v1 = afl.getR(i) - newPos;
247//
248// v2 = afl.getR(i + 1) - newPos;
249//
250// vv = afl.getR(i + 1) - afl.getR(i);
251//
252// double dst = v1.length2() + v2.length2();
253// if (dst < dist2)
254// {
255// dist2 = dst;
256// panThrough = i;
257// }
258//
259// cs = (v1 & v2);
260// sn = (v1 ^ v2);
261//
262// angle += atan2(sn, cs);
263// }//for i
264//
265// hit = ((angle > 3.00) || (angle < -3.00));
266//
267// return hit;
268//}//MoveInsideMovingBoundary(...)
269
270
271
272//Проверка пересечения вихрями следа профиля при перемещении
273void Wake::Inside(const std::vector<Point2D>& newPos, Airfoil& afl, bool isMoves, const AirfoilGeometry& oldAfl)
274{
275 std::vector<double> gamma;
276 gamma.resize(afl.getNumberOfPanels(), 0.0);
277
278 std::vector<int> through;
279 through.resize(vtx.size(), -1);
280
281//#pragma omp parallel for default(none) shared(afl, oldAfl, isMoves, through, newPos)
282 for (int i = 0; i < (int)vtx.size(); ++i)
283 {
284 size_t minN;
285
286 bool crit = isMoves ? MoveInsideMovingBoundary(newPos[i], vtx[i].r(), oldAfl, afl, minN) : MoveInside(newPos[i], vtx[i].r(), afl, minN);
287
288 if (crit)
289 through[i] = (int)minN;
290 }//for i
291
292
293 //std::stringstream sss;
294 //sss << "through_" << W.currentStep;
295 //std::ofstream of(W.getPassport().dir + "dbg/" + sss.str());
296 //for (size_t i = 0; i < gamma.size(); ++i)
297 // of << gamma[i] << std::endl;
298 //of.close();
299
300 //std::stringstream sss;
301 //sss << "through_" << W.currentStep;
302 //std::ofstream of(W.getPassport().dir + "dbg/" + sss.str());
303 //for (size_t i = 0; i < gamma.size(); ++i)
304 // of << through[i] << std::endl;
305 //of.close();
306
307 for (size_t q = 0; q < through.size(); ++q)
308 if (through[q] > -1)
309 {
310 gamma[through[q]] += vtx[q].g();
311 vtx[q].g() = 0.0;
312 }
313
314 afl.gammaThrough = gamma;
315}//Inside(...)
316
317//Поиск ближайшего соседа
318void Wake::GetPairs(int type)
319{
320 GetPairsBS(type);
321}
322
323//Поиск ближайшего соседа
324void Wake::GetPairsBS(int type)
325{
326 neighb.resize(vtx.size(), 0);
327
328 Point2D Ri, Rk;
329
331
332#pragma omp parallel for default(none) shared(type, maxG) schedule(dynamic, DYN_SCHEDULE)
333 for (int i = 0; i < vtx.size(); ++i)
334 {
335 neighb[i] = 0;
336 int s = i;
337 const Vortex2D& vtxI = vtx[i];
338
339 bool found = false;
340
341 double r2, r2test;
342
343 const double& cSP = collapseScaleParameter;
344 const double& cRBP = collapseRightBorderParameter;
345
346 while ( (!found) && ( s + 1 < (int)vtx.size() ) )
347 {
348 s++;
349 const Vortex2D& vtxK = vtx[s];
350
351 r2 = dist2(vtxI.r(), vtxK.r());
352
353 //линейное увеличение радиуса коллапса
354 double mnog = std::max(1.0, /* 2.0 * */ (vtxI.r()[0] - cRBP) / cSP);
355
356 r2test = sqr( W.getPassport().wakeDiscretizationProperties.epscol * mnog );
357
358 if (type == 1)
359 r2test *= 4.0; //Увеличение радиуса коллапса в 2 раза для коллапса вихрей разных знаков
360
361 if (r2 < r2test)
362 {
363 switch (type)
364 {
365 case 0:
366 found = ( vtxI.g()*vtxK.g() != 0.0) && (fabs(vtxI.g() + vtxK.g()) < sqr(mnog) * maxG);
367 break;
368 case 1:
369 found = (vtxI.g()*vtxK.g() < 0.0);
370 break;
371 case 2:
372 found = (vtxI.g()*vtxK.g() > 0.0) && (fabs(vtxI.g() + vtxK.g()) < sqr(mnog) * maxG);
373 break;
374 }
375 }//if r2 < r2_test
376 }//while
377
378 if (found)
379 neighb[i] = s;
380 }//for locI
381}//GetPairsBS(...)
382
383
384//Поиск ближайшего соседа
386{
387 neighb.resize(vtx.size(), 0);
388
389 Point2D Ri, Rk;
390
392
393#pragma omp parallel for default(none) shared(type, maxG) schedule(dynamic, DYN_SCHEDULE)
394 for (int i = 0; i < vtx.size(); ++i)
395 {
396 neighb[i] = 0;
397 int s = i;
398 const Vortex2D& vtxI = vtx[i];
399
400 bool found = false;
401
402 double r2, r2test = 1e+10;
403
404 const double& cSP = collapseScaleParameter;
405 const double& cRBP = collapseRightBorderParameter;
406
407 while (/*(!found) &&*/ (s + 1 < (int)vtx.size()))
408 {
409 s++;
410 const Vortex2D& vtxK = vtx[s];
411
412 r2 = dist2(vtxI.r(), vtxK.r());
413
414 if (r2 < r2test)
415 {
416 r2test = r2;
417 neighb[i] = s;
418 }//if r2 < r2_test
419 }//while
420
421 }//for locI
422}//GetPairsClosestNeib(...)
423
424
425
426#if defined(USE_CUDA)
427void Wake::GPUGetPairs(int type)
428{
429 size_t npt = vtx.size();
430
431 double tCUDASTART = 0.0, tCUDAEND = 0.0;
432
433 tCUDASTART += omp_get_wtime();
434 std::vector<int> tnei(npt, 0);
435 neighb.resize(npt, 0);
436
437 if (npt > 0)
438 {
439 cuCalculatePairs(npt, devVtxPtr, devMeshPtr, devNeiPtr, 2.0*W.getPassport().wakeDiscretizationProperties.epscol, sqr(W.getPassport().wakeDiscretizationProperties.epscol), type);
440 W.getCuda().CopyMemFromDev<int, 1>(npt, devNeiPtr, neighb.data(), 30);
441 }
442
443 tCUDAEND += omp_get_wtime();
444
445 //W.getInfo('t') << "GPU_Pairs: " << (tCUDAEND - tCUDASTART) << std::endl;
446}
447#endif
448
449#if defined(USE_CUDA)
450void Wake::GPUGetPairsClosestNeib(int type)
451{
452 size_t npt = vtx.size();
453
454 double tCUDASTART = 0.0, tCUDAEND = 0.0;
455
456 tCUDASTART += omp_get_wtime();
457 std::vector<int> tnei(npt, 0);
458 neighb.resize(npt, 0);
459
460 if (npt > 0)
461 {
462 cuCalculatePairsClosestNeib(npt, devVtxPtr, devMeshPtr, devNeiPtr, 2.0 * W.getPassport().wakeDiscretizationProperties.epscol, sqr(W.getPassport().wakeDiscretizationProperties.epscol), type);
463
464 W.getCuda().CopyMemFromDev<int, 1>(npt, devNeiPtr, neighb.data(), 30);
465 }
466
467 tCUDAEND += omp_get_wtime();
468
469 //W.getInfo('t') << "GPU_Pairs: " << (tCUDAEND - tCUDASTART) << std::endl;
470}
471#endif
472
473
474// Коллапс вихрей
475int Wake::Collaps(int type, int times)
476{
477 int nHlop = 0; //общее число убитых вихрей
478
479 //int loc_hlop = 0; //
480
481 std::vector<bool> flag; //как только вихрь сколлапсировался flag = 1
482 for (int z = 0; z < times; ++z)
483 {
484
485#if (defined(__CUDACC__) || defined(USE_CUDA)) && (defined(CU_PAIRS))
488 //if (W.getPassport().numericalSchemes.velocityComputation.second == 0)
489 {
490 const_cast<Gpu&>(W.getCuda()).RefreshWake(3);
491
492 double NeibStart = omp_get_wtime();
493 GPUGetPairs(type);
494 double NeibFinish = omp_get_wtime();
495 std::cout << "GPU_direct_closest_neib = " << NeibFinish - NeibStart << std::endl;
496 }
498
499 //else
500 //{
501 // double NeibStart = omp_get_wtime();
502 // GetPairs(type);
503 // //GetPairsClosestNeib(type);
504 // double NeibFinish = omp_get_wtime();
505 // std::cout << "CPU_direct_closest_neib = " << NeibFinish - NeibStart << std::endl;
506 //}
507#else
508 GetPairs(type);
509#endif
510
511
512 //std::ofstream neibFile(W.getPassport().dir + "neib" + std::to_string(W.currentStep));
513 //for (size_t q = 0; q < neighb.size(); ++q)
514 // neibFile << q << " " << neighb[q] << std::endl;
515 //neibFile.close();
516
517
518 //loc_hlop = 0;//число схлопнутых вихрей
519//*
520 flag.clear();
521 flag.resize(vtx.size(), false);
522
523 double sumAbsGam, iws;
524 Point2D newPos;
525
526 for (size_t vt = 0; vt + 1 < vtx.size(); ++vt)
527 {
528 Vortex2D& vtxI = vtx[vt];
529
530 if (!flag[vt])
531 {
532 int ssd = neighb[vt];
533
534 if ((ssd < 0) || (ssd >= flag.size()))
535 std::cout << "ssd = " << ssd << ", flag.size() = " << flag.size() << ", vtx.size() = " << vtx.size() << std::endl;
536
537 if ((ssd != 0) && (!flag[ssd]))
538 {
539 Vortex2D& vtxK = vtx[ssd];
540
541 flag[ssd] = true;
542
543 Vortex2D sumVtx;
544 sumVtx.g() = vtxI.g() + vtxK.g();
545
546 switch (type)
547 {
548 case 0:
549 case 2:
550 sumAbsGam = fabs(vtxI.g()) + fabs(vtxK.g());
551
552 iws = sumAbsGam > 1e-10 ? 1.0 / sumAbsGam : 1.0;
553
554 sumVtx.r() = (vtxI.r() * fabs(vtxI.g()) + vtxK.r() * fabs(vtxK.g())) * iws;
555 break;
556
557 case 1:
558 sumVtx.r() = (fabs(vtxI.g()) > fabs(vtxK.g())) ? vtxI.r() : vtxK.r();
559 break;
560 }
561
562 bool fl_hit = true;
563 //double Ch1[2];
564 size_t hitpan = -1;
565
566 for (size_t afl = 0; afl < W.getNumberOfAirfoil(); ++afl)
567 {
568 //проверим, не оказался ли новый вихрь внутри контура
569 if (MoveInside(sumVtx.r(), vtxI.r(), W.getAirfoil(afl), hitpan) || MoveInside(sumVtx.r(), vtxK.r(), W.getAirfoil(afl), hitpan))
570 fl_hit = false;
571 }//for
572
573 if (fl_hit)
574 {
575 vtxI = sumVtx;
576 vtxK.g() = 0.0;
577 nHlop++;
578 }//if (fl_hit)
579
580 }//if ((ssd!=0)&&(!flag[ssd]))
581 }//if !flag[vt]
582 }//for vt
583//*/
584 }//for z
585
586 return nHlop;
587}//Collaps(...)
588
589// Коллапс вихрей
590int Wake::CollapsNew(int type, int times)
591{
592 int nHlop = 0; //общее число убитых вихрей
593
594 const double& cSP = collapseScaleParameter;
595 const double& cRBP = collapseRightBorderParameter;
596
598
599 //int loc_hlop = 0; // neigh
600
601 std::vector<bool> flag; //как только вихрь сколлапсировался flag = 1
602
603 for (int z = 0; z < times; ++z)
604 {
605 //loc_hlop = 0;//число схлопнутых вихрей
606
607 flag.clear();
608 flag.resize(vtx.size(), false);
609
610 double sumAbsGam, iws;
611 Point2D newPos;
612
613 for (size_t vt = 0; vt + 1 < vtx.size(); ++vt)
614 {
615 Vortex2D& vtxI = vtx[vt];
616 if (!flag[vt])
617 {
618
619 for (int s = 0; s < knbForRestruct; ++s)
620 {
621 //int ssd = neighb[vt];
622 int ssd = neighbNew[vt * knbForRestruct + s];
623 if (ssd == 0)
624 continue;
625
626 Vortex2D& vtxK = vtx[ssd];
627
628 if (/*(ssd != 0) &&*/ (!flag[ssd])) //ssd==0 always true
629 {
630
631
632 Vortex2D sumVtx;
633 sumVtx.g() = vtxI.g() + vtxK.g();
634
635 switch (type)
636 {
637 case 0:
638 case 2:
639 sumAbsGam = fabs(vtxI.g()) + fabs(vtxK.g());
640 iws = sumAbsGam > 1e-10 ? 1.0 / sumAbsGam : 1.0;
641 sumVtx.r() = (vtxI.r() * fabs(vtxI.g()) + vtxK.r() * fabs(vtxK.g())) * iws;
642 break;
643
644 case 1:
645 sumVtx.r() = (fabs(vtxI.g()) > fabs(vtxK.g())) ? vtxI.r() : vtxK.r();
646 break;
647 }
648
649 bool fl_hit = true;
650 size_t hitpan = -1;
651
652 for (size_t afl = 0; afl < W.getNumberOfAirfoil(); ++afl)
653 {
654 //проверим, не оказался ли новый вихрь внутри контура
655 if (MoveInside(sumVtx.r(), vtxI.r(), W.getAirfoil(afl), hitpan) || MoveInside(sumVtx.r(), vtxK.r(), W.getAirfoil(afl), hitpan))
656 {
657 //std::cout << "HIT" << std::endl;
658 fl_hit = false;
659 }
660 }//for
661
662 if (fl_hit)
663 {
664 vtxI = sumVtx;
665 vtxK.g() = 0.0;
666 nHlop++;
667 flag[vt] = true;
668 flag[ssd] = true;
669 }//if (fl_hit)
670
671 }//if ((ssd!=0)&&(!flag[ssd]))
672
673 if (flag[vt])
674 break;
675 }// for s
676 }//if !flag[vt]
677 }//for vt
678
679 }//for z
680
681 return nHlop;
682}//CollapsNew(...)
683
684// Коллапс вихрей
685int Wake::CollapsNewFast(int type, int times, std::vector<Vortex2D>& ri, std::vector<Vortex2D>& rj, std::vector<Point2D>& rnew, std::vector<std::pair<int, int>>& rindex)
686{
687 int nHlop = 0; //общее число убитых вихрей
688
689 const double& cSP = collapseScaleParameter;
690 const double& cRBP = collapseRightBorderParameter;
691
693
694 //int loc_hlop = 0; // neigh
695
696 std::vector<bool> flag; //как только вихрь сколлапсировался flag = 1
697
698 for (int z = 0; z < times; ++z)
699 {
700 //loc_hlop = 0;//число схлопнутых вихрей
701
702 flag.clear();
703 flag.resize(vtx.size(), false);
704
705 double sumAbsGam, iws;
706 Point2D newPos;
707
708 for (size_t vt = 0; vt + 1 < vtx.size(); ++vt)
709 {
710 Vortex2D& vtxI = vtx[vt];
711 if (!flag[vt])
712 {
713
714 for (int s = 0; s < knbForRestruct; ++s)
715 {
716 //int ssd = neighb[vt];
717 int ssd = neighbNew[vt * knbForRestruct + s];
718 if (ssd == 0)
719 continue;
720
721 Vortex2D& vtxK = vtx[ssd];
722
723 if (/*(ssd != 0) &&*/ (!flag[ssd])) //ssd==0 always true
724 {
725
726
727 Vortex2D sumVtx;
728 sumVtx.g() = vtxI.g() + vtxK.g();
729
730 switch (type)
731 {
732 case 0:
733 case 2:
734 sumAbsGam = fabs(vtxI.g()) + fabs(vtxK.g());
735 iws = sumAbsGam > 1e-10 ? 1.0 / sumAbsGam : 1.0;
736 sumVtx.r() = (vtxI.r() * fabs(vtxI.g()) + vtxK.r() * fabs(vtxK.g())) * iws;
737 break;
738
739 case 1:
740 sumVtx.r() = (fabs(vtxI.g()) > fabs(vtxK.g())) ? vtxI.r() : vtxK.r();
741 break;
742 }
743
744 bool fl_hit = true;
745 size_t hitpan = -1;
746
747 //for (size_t afl = 0; afl < W.getNumberOfAirfoil(); ++afl) //кол-во профилей
748 //{
749 ri.push_back(Vortex2D(vtxI.r(), vtxI.g()));
750 rj.push_back(Vortex2D(vtxK.r(), vtxK.g()));
751 rindex.push_back({ (int)vt, ssd });
752 rnew.push_back(sumVtx.r());
753 nHlop++;
754 flag[vt] = true;
755 flag[ssd] = true;
756 //проверим, не оказался ли новый вихрь внутри контура
757 //if (MoveInside(sumVtx.r(), vtxI.r(), W.getAirfoil(afl), hitpan) || MoveInside(sumVtx.r(), vtxK.r(), W.getAirfoil(afl), hitpan))
758 //{
759 // //std::cout << "HIT" << std::endl;
760 // fl_hit = false;
761 //}
762 //}//for
763
764 //if (fl_hit)
765 //{
766 // vtxI = sumVtx;
767 // vtxK.g() = 0.0;
768 // //nHlop++;
769 // flag[vt] = true;
770 // flag[ssd] = true;
771 //}//if (fl_hit)
772
773 }//if ((ssd!=0)&&(!flag[ssd]))
774
775 if (flag[vt])
776 break;
777 }// for s
778 }//if !flag[vt]
779 }//for vt
780
781 }//for z
782
783 return nHlop;
784}//CollapsNew(...)
785
786
787
789{
790 int nFar = 0;
791 double distFar2 = sqr(W.getPassport().wakeDiscretizationProperties.distFar);
793 Point2D zerovec = { 0.0, 0.0 };
794#pragma omp parallel for default(none) shared(distFar2, zerovec) reduction(+:nFar)
795 for (int i = 0; i <static_cast<int>(vtx.size()); ++i)
796 {
797 if (dist2(vtx[i].r(), zerovec) > distFar2)
798 {
799 vtx[i].g() = 0.0;
800 nFar++;
801 }
802 }
803
804 return nFar;
805}//RemoveFar()
806
807
809{
810 const double porog_g = 1e-15;
811
812 std::vector<Vortex2D/*, VM2D::MyAlloc<VMlib::Vortex2D>*/> newWake;
813
814 newWake.reserve(vtx.size());
815
816 for (size_t q = 0; q < vtx.size(); ++q)
817 if (fabs(vtx[q].g()) > porog_g)
818 newWake.push_back(vtx[q]);
819
820 size_t delta = vtx.size() - newWake.size();
821
822 newWake.swap(vtx);
823
824 return delta;
825}//RemoveZero()
826
827
828
829
830//Реструктуризация вихревого следа
832{
833 W.getTimers().start("Restr");
834
836 {
837 double timePreKnn = -omp_get_wtime();
838
839 // Определение параметров, отвечающих за увеличение радиуса коллапса
840 std::vector<double> rightBorder, horizSpan;
841 rightBorder.reserve(W.getNumberOfAirfoil());
842 horizSpan.reserve(W.getNumberOfAirfoil());
843
844 for (size_t q = 0; q < W.getNumberOfAirfoil(); ++q)
845 {
846
847 rightBorder.emplace_back(W.getAirfoil(q).upRight[0]);
848 horizSpan.emplace_back(W.getAirfoil(q).upRight[0] - W.getAirfoil(q).lowLeft[0]);
849 }
850
851 if (W.getNumberOfAirfoil() > 0)
852 {
853 W.getNonConstWake().collapseRightBorderParameter = *std::max_element(rightBorder.begin(), rightBorder.end());
854 W.getNonConstWake().collapseScaleParameter = *std::max_element(horizSpan.begin(), horizSpan.end());
855 }
856 else
857 {
860 }
861#if defined(__CUDACC__) || defined(USE_CUDA)
863#endif
864
865 timePreKnn += omp_get_wtime();
866
867 //std::cout << "timePreKnn = " << timePreKnn * 1000 << " ms" << std::endl;
868
869
871
872
874 {
875 //CPU/GPU - прямой алгоритм
876 //double timeCollaps = -omp_get_wtime();
877 Collaps(0, 1);
878 //Collaps(1, 1);
879 //Collaps(2, 1);
880 //timeCollaps += omp_get_wtime();
881 //std::cout << "timeCollaps = " << timeCollaps * 1000 << " ms" << std::endl;
882 }
883 else
884 {
885 //быстрый алгоритм
886 for (int collapsStep = 0; collapsStep < 1; ++collapsStep)
887 {
888 //ttB = -omp_get_wtime();
889
890 double timeResize = -omp_get_wtime();
891 neighbNew.resize(vtx.size() * knbForRestruct);
892 timeResize += omp_get_wtime();
893
894 //std::cout << "timeResize = " << timeResize * 1000 << " ms" << std::endl;
895
896 const double& cSP = collapseScaleParameter;
897 const double& cRBP = collapseRightBorderParameter;
898 const double& maxG = W.getPassport().wakeDiscretizationProperties.maxGamma;
899 const double& epsCol = W.getPassport().wakeDiscretizationProperties.epscol;
900
901#ifndef USE_CUDA
902 //CPU
903 std::vector<std::vector<std::pair<double, size_t>>> initdist(vtx.size());
904 for (auto& d : initdist)
905 d.resize(2 * knbForRestruct, { -1.0, -1 });
906 double timeKnn = -omp_get_wtime();
907 WakekNNnewForCollaps(vtx, knbForRestruct, initdist, cSP, cRBP, maxG, epsCol, collapsStep);//CPU
908 timeKnn += omp_get_wtime();
909 //std::cout << "Time_Knn_CPU = " << timeKnn * 1000 << " ms" << std::endl;
910#else
911 double timeAlloc = -omp_get_wtime();
912 std::vector<std::pair<double, size_t>> initdistcuda(knbForRestruct * vtx.size()); //CUDA
913 timeAlloc += omp_get_wtime();
914 //std::cout << "timeAlloc = " << timeAlloc * 1000 << " ms" << std::endl;
915
916 double timeKnn = -omp_get_wtime();
917 W.getNonConstCuda().RefreshWake(5);
918
921 //Построение дерева для вихрей
922 BHcu::CudaTreeInfo knnTree(W.getCuda().blocks, tree_T::contr, object_T::point3, scheme_T::noScheme, false);
923 knnTree.MemoryAllocate((int)W.getCuda().n_CUDA_wake);
924 knnTree.Update((int)W.getWake().vtx.size(), W.getWake().devVtxPtr, W.getPassport().wakeDiscretizationProperties.eps);
925 knnTree.Build();
926
927 Point2D minr, maxr;
928
929 cudaMemcpy(&minr, knnTree.minrD, sizeof(Point2D), cudaMemcpyDeviceToHost);
930 cudaMemcpy(&maxr, knnTree.maxrD, sizeof(Point2D), cudaMemcpyDeviceToHost);
931
932
933 kNNcuda<knbForRestruct>(minr, maxr, W.getCuda().blocks, vtx, initdistcuda, vecForKnn, cSP, cRBP, maxG, epsCol, collapsStep); //CUDA
935 timeKnn += omp_get_wtime();
936 //std::cout << "Time_Knn_GPU = " << timeKnn * 1000 << " ms" << std::endl;
937#endif
938
939 double timeCopy = -omp_get_wtime();
940#pragma omp parallel for
941 for (int i = 0; i < vtx.size(); ++i)
942 {
943#ifndef USE_CUDA
944 for (int j = 0; j < knbForRestruct; ++j)
945 neighbNew[i * knbForRestruct + j] = (int)initdist[i][j].second;
946#else
947 for (int j = 0; j < knbForRestruct; ++j)
948 neighbNew[i * knbForRestruct + j] = (int)initdistcuda[i * knbForRestruct + j].second;
949#endif
950 }
951 timeCopy += omp_get_wtime();
952 //std::cout << "timeCopy = " << timeCopy * 1000 << " ms" << std::endl;
953
954 //
955 //std::ofstream initdistFile(W.getPassport().dir + "initdist-GPU" + std::to_string(W.currentStep));
956 //for (int i = 0; i < vtx.size(); ++i)
957 //{
958 // initdistFile << i;
959 // for (int k = 0; k < knb; ++k)
960 // initdistFile << " " << neighbNew[i * (knb)+k] << " " << (vtx[i].r() - vtx[neighbNew[i * (knb)+k]].r()).length();
961 // initdistFile << std::endl;
962 //
963 // //for (int k = 0; k < knb; ++k)
964 // //initdistFile << " " << neighbNew[i * (knb)+k] << " " << vtx[neighbNew[i * (knb)+k]].g() << " " << (vtx[neighbNew[i * (knb)+k]].r() - vtx[i].r()).length() << "; ";
965 // //initdistFile << std::endl;
966 //}
967 //initdistFile.close();
968 //
969
970
971 double timeReserve = -omp_get_wtime();
972 std::vector<Vortex2D> ri, rj;
973 std::vector<Point2D> rnew;
974 std::vector<std::pair<int, int>> rindex;
975 ri.reserve(vtx.size());
976 rj.reserve(vtx.size());
977 rnew.reserve(vtx.size());
978
979 rindex.reserve(vtx.size());
980 timeReserve += omp_get_wtime();
981 //std::cout << "timeReserve = " << timeReserve * 1000 << " ms" << std::endl;
982
983
984#ifdef USE_CUDA
985 //поиск пар для объединения
986 int nHlop = CollapsNewFast(0, 1, ri, rj, rnew, rindex); //тут заполняются ri, rj, rnew, rindex
987
988 //формирование траекторий объединямых вихрей
989 std::vector<std::pair<Point2D, Point2D>> segments(2 * ri.size());
990 for (size_t i = 0; i < ri.size(); ++i)
991 {
992 segments[2 * i + 0].first = ri[i];
993 segments[2 * i + 0].second = rnew[i];
994
995 segments[2 * i + 1].first = rj[i];
996 segments[2 * i + 1].second = rnew[i];
997 }
998
999 //контроль протыкания
1000 std::vector<int> hit(segments.size());
1001 if (ri.size() > 0)
1002 {
1003 double* devSegments_ptr;
1004
1005 cudaMalloc(&devSegments_ptr, segments.size() * sizeof(double) * 4);
1006 cudaMemcpy(devSegments_ptr, segments.data(), segments.size() * sizeof(double) * 4, cudaMemcpyHostToDevice);
1007
1009 auto& cntrTreeSeg = *W.getCuda().cntrTreeSegment;
1010 cntrTreeSeg.MemoryAllocate((int)W.getCuda().n_CUDA_wake);
1011 cntrTreeSeg.UpdatePanelGeometry((int)segments.size(), (double4*)devSegments_ptr);
1012 cntrTreeSeg.Build();
1013
1014 BHcu::treePanelsSegmentsIntersectionCalculationWrapper(*W.getNonConstCuda().auxTreePnl, cntrTreeSeg,
1015 W.getWake().devNearestPanelPtr);
1016
1017 W.timerInside.stop();
1018
1019 cudaMemcpy(hit.data(), W.getWake().devNearestPanelPtr, segments.size() * sizeof(int), cudaMemcpyDeviceToHost);
1020 cudaFree(devSegments_ptr);
1021 }//if ri.size() > 0
1022
1023 for (size_t i = 0; i < ri.size(); ++i)
1024 {
1025 if (hit[2 * i + 0] == -1 && hit[2 * i + 1] == -1)
1026 {
1027 vtx[rindex[i].first].r() = rnew[i];
1028 vtx[rindex[i].first].g() += vtx[rindex[i].second].g();
1029
1030 vtx[rindex[i].second].g() = 0.0;
1031 }
1032 }
1033
1034#else
1035 double timeCollaps;
1036 //CollapsNew(0, 1);
1037 timeCollaps = -omp_get_wtime();
1038 int nHlop = CollapsNewFast(0, 1, ri, rj, rnew, rindex);
1039
1040 timeCollaps += omp_get_wtime();
1041
1042 size_t hitA, hitB;
1043 bool ifInside;
1044#pragma omp parallel for private (hitA, hitB, ifInside)
1045 for (int i = 0; i < (int)ri.size(); ++i)
1046 {
1047 ifInside = false;
1048 for (size_t q = 0; q < W.getNumberOfAirfoil(); ++q)
1049 {
1050 MoveInside(rnew[i], ri[i], W.getAirfoil(q), hitA);
1051 MoveInside(rnew[i], rj[i], W.getAirfoil(q), hitB);
1052 if ((hitA != size_t(-1)) || (hitB != size_t(-1)))
1053 ifInside = true;
1054 }
1055 if (!ifInside)
1056 {
1057 vtx[rindex[i].first].r() = rnew[i];
1058 vtx[rindex[i].first].g() += vtx[rindex[i].second].g();
1059
1060 vtx[rindex[i].second].g() = 0.0;
1061
1062 }
1063 }
1064#endif
1065 }
1066 }
1067
1068 }
1069 double timeRemove = -omp_get_wtime();
1070 RemoveFar();
1071 RemoveZero();
1072 timeRemove += omp_get_wtime();
1073 //std::cout << "timeRemove = " << timeRemove * 1000 << std::endl;
1074
1075 W.getTimers().stop("Restr");
1076//*
1077// std::cout << "Pre-knn time = " << timePreKnn * 1000.0 << " ms" << std::endl;
1078// std::cout << "Pre-knn resize time = " << timeResize * 1000.0 << " ms" << std::endl;
1079// std::cout << "Time_Knn_GPU = " << timeKnn * 1000 << " ms" << std::endl;
1080// std::cout << "Post-knn copy time = " << timeCopy * 1000.0 << " ms" << std::endl;
1081// std::cout << "Reserve_time = " << timeReserve * 1000 << " ms" << std::endl;
1082// std::cout << "Collaps_time = " << timeCollaps * 1000 << " ms" << std::endl;
1083// std::cout << "Copy2_time = " << timeCopy2 * 1000 << " ms" << std::endl;
1084// std::cout << "Ray_time = " << timeRay * 1000 << " ms" << std::endl;
1085// std::cout << "Collaps2_time = " << timeCollaps2 * 1000 << " ms" << std::endl;
1086// std::cout << "Remove_time = " << timeRemove * 1000 << " ms" << std::endl;
1087
1088
1089//*/
1090}//Restruct()
Заголовочный файл с описанием класса Airfoil.
Заголовочный файл с описанием класса Boundary.
Заголовочный файл с описанием класса MeasureVP.
Заголовочный файл с описанием класса Mechanics.
Заголовочный файл с описанием класса Preprocessor.
Заголовочный файл с описанием класса StreamParser.
Заголовочный файл с описанием класса Velocity.
Заголовочный файл с описанием класса Wake.
Заголовочный файл с описанием класса World2D.
Класс, определяющий форму профиля
Definition Airfoil2D.h:66
const Point2D & getR(size_t q) const
Возврат константной ссылки на вершину профиля
Definition Airfoil2D.h:113
bool inverse
Признак разворота нормалей (для расчета внутренних течений)
Definition Airfoil2D.h:76
size_t getNumberOfPanels() const
Возврат количества панелей на профиле
Definition Airfoil2D.h:163
Абстрактный класс, определяющий обтекаемый профиль
Definition Airfoil2D.h:182
bool isOutsideGabarits(const Point2D &r) const
Определяет, находится ли точка с радиус-вектором вне габаритного прямоугольника профиля
Definition Airfoil2D.cpp:87
Point2D upRight
Правый верхний угол габаритного прямоугольника профиля
Definition Airfoil2D.h:271
std::vector< double > gammaThrough
Суммарные циркуляции вихрей, пересекших панели профиля на прошлом шаге
Definition Airfoil2D.h:276
Point2D lowLeft
Левый нижний угол габаритного прямоугольника профиля
Definition Airfoil2D.h:270
Класс, обеспечивающий возможность выполнения вычислений на GPU по технологии Nvidia CUDA.
Definition Gpu2D.h:69
void setCollapseCoeff(double pos_, double refLength_)
Установка правой границы самого правого профиля (для организации увеличения радиуса коллапса)
Definition Gpu2D.h:284
WakeDiscretizationProperties wakeDiscretizationProperties
Структура с параметрами дискретизации вихревого следа
Definition Passport2D.h:292
NumericalSchemes numericalSchemes
Структура с используемыми численными схемами
Definition Passport2D.h:295
std::vector< Vortex2D > vtx
Список вихревых элементов
const World2D & W
Константная ссылка на решаемую задачу
double collapseRightBorderParameter
абсцисса, правее которой происходит линейный (вправо) рост радиуса коллапса
Definition Wake2D.h:161
void GetPairsBS(int type)
Definition Wake2D.cpp:324
bool MoveInside(const Point2D &newPos, const Point2D &oldPos, const Airfoil &afl, size_t &panThrough) const
Проверка проникновения точки через границу профиля
Definition Wake2D.cpp:62
int RemoveFar()
Зануление далеко улетевших вихрей
Definition Wake2D.cpp:788
void GetPairsClosestNeib(int type)
Definition Wake2D.cpp:385
double collapseScaleParameter
характерный масштаб, на котором происходит рост радиуса коллапса
Definition Wake2D.h:164
std::vector< int > neighbNew
Definition Wake2D.h:69
std::vector< int > neighb
Вектор потенциальных соседей для будущего коллапса
Definition Wake2D.h:66
void GetPairs(int type)
Поиск ближайшего соседа
Definition Wake2D.cpp:318
size_t RemoveZero()
Исключение нулевых и мелких вихрей
Definition Wake2D.cpp:808
bool MoveInsideMovingBoundary(const Point2D &newPos, const Point2D &oldPos, const AirfoilGeometry &oldAfl, const Airfoil &afl, size_t &panThrough) const
Проверка проникновения точки через границу профиля
Definition Wake2D.cpp:178
void Restruct()
Реструктуризация вихревого следа
Definition Wake2D.cpp:831
int Collaps(int type, int times)
Коллапс вихрей
Definition Wake2D.cpp:475
void Inside(const std::vector< Point2D > &newPos, Airfoil &afl, bool isMoves, const AirfoilGeometry &oldAfl)
Проверка пересечения вихрями следа профиля при перемещении
Definition Wake2D.cpp:273
int CollapsNewFast(int type, int times, std::vector< Vortex2D > &ri, std::vector< Vortex2D > &rj, std::vector< Point2D > &rnew, std::vector< std::pair< int, int > > &rindex)
Definition Wake2D.cpp:685
int CollapsNew(int type, int times)
Definition Wake2D.cpp:590
size_t getNumberOfAirfoil() const
Возврат количества профилей в задаче
Definition World2D.h:174
VMlib::vmTimer timerMerging
Definition World2D.h:359
const Wake & getWake() const
Возврат константной ссылки на вихревой след
Definition World2D.h:226
const Airfoil & getAirfoil(size_t i) const
Возврат константной ссылки на объект профиля
Definition World2D.h:157
VMlib::vmTimer timerInside
Definition World2D.h:358
Gpu & getNonConstCuda() const
Возврат неконстантной ссылки на объект, связанный с видеокартой (GPU)
Definition World2D.h:266
VMlib::TimersGen & getTimers() const
Возврат ссылки на временную статистику выполнения шага расчета по времени
Definition World2D.h:276
const Gpu & getCuda() const
Возврат константной ссылки на объект, связанный с видеокартой (GPU)
Definition World2D.h:261
const Passport & getPassport() const
Возврат константной ссылки на паспорт
Definition World2D.h:251
Wake & getNonConstWake() const
Возврат неконстантной ссылки на вихревой след
Definition World2D.h:231
void start(const std::string &timerLabel)
Запуск счетчика
Definition TimesGen.cpp:55
Класс, опеделяющий двумерный вихревой элемент
Definition Vortex2D.h:59
HD Point2D & r()
Функция для доступа к радиус-вектору вихря
Definition Vortex2D.h:87
HD double & g()
Функция для доступа к циркуляции вихря
Definition Vortex2D.h:95
auto length2() const -> typename std::remove_const< typename std::remove_reference< decltype(this->data[0])>::type >::type
Вычисление квадрата нормы (длины) вектора
Definition numvector.h:386
const vmTimer & stop() const
Останов работающего счетчика времени
Definition TimesGen.h:125
const vmTimer & start() const
Запуск (первый или повторный) счетчика времени
Definition TimesGen.h:109
const vmTimer & reset() const
Сброс счетчика времени
Definition TimesGen.h:101
void WakekNNnewForCollaps(const std::vector< Vortex2D > &vtx, const size_t k, std::vector< std::vector< std::pair< double, size_t > > > &initdist, double cSP, double cRBP, double maxG, double epsCol, int type)
Definition knnCPU.cpp:307
knn CPU.
T sqr(T x)
Умножение a на комплексно сопряженноe к b.
Definition defsBH.h:101
std::pair< std::string, int > velocityComputation
Definition Passport2D.h:180
double eps
Радиус вихря
Definition Passport2D.h:126
double epscol
Радиус коллапса
Definition Passport2D.h:132
double maxGamma
Максимально допустимая циркуляция вихря
Definition Passport2D.h:144
double distFar
Расстояние от центра самого подветренного (правого) профиля, на котором вихри уничтожаются
Definition Passport2D.h:135