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