49#include "treeKernels.cuh"
88 if ((((aflRj - oldPos) ^ (newPos - oldPos)) * ((aflRj1 - oldPos) ^ (newPos - oldPos)) <= 0) && \
89 (((oldPos - aflRj) ^ (aflRj1 - aflRj)) * ((newPos - aflRj) ^ (aflRj1 - aflRj)) <= 0))
193 double dist2 = 1000000000.0;
198 v1 = afl.
getR(i) - newPos;
200 v2 = afl.
getR(i + 1) - newPos;
214 angle += atan2(sn, cs);
217 hit = ((angle > 3.14) || (angle < -3.14));
275 std::vector<double> gamma;
278 std::vector<int> through;
279 through.resize(
vtx.size(), -1);
282 for (
int i = 0; i < (int)
vtx.size(); ++i)
289 through[i] = (int)minN;
307 for (
size_t q = 0; q < through.size(); ++q)
310 gamma[through[q]] +=
vtx[q].g();
332#pragma omp parallel for default(none) shared(type, maxG) schedule(dynamic, DYN_SCHEDULE)
333 for (
int i = 0; i <
vtx.size(); ++i)
346 while ( (!found) && ( s + 1 < (
int)
vtx.size() ) )
351 r2 = dist2(vtxI.
r(), vtxK.
r());
354 double mnog = std::max(1.0, (vtxI.
r()[0] - cRBP) / cSP);
366 found = ( vtxI.
g()*vtxK.
g() != 0.0) && (fabs(vtxI.
g() + vtxK.
g()) < sqr(mnog) * maxG);
369 found = (vtxI.
g()*vtxK.
g() < 0.0);
372 found = (vtxI.
g()*vtxK.
g() > 0.0) && (fabs(vtxI.
g() + vtxK.
g()) < sqr(mnog) * maxG);
393#pragma omp parallel for default(none) shared(type, maxG) schedule(dynamic, DYN_SCHEDULE)
394 for (
int i = 0; i <
vtx.size(); ++i)
402 double r2, r2test = 1e+10;
407 while ( (s + 1 < (
int)
vtx.size()))
412 r2 = dist2(vtxI.
r(), vtxK.
r());
427void Wake::GPUGetPairs(
int type)
429 size_t npt =
vtx.size();
431 double tCUDASTART = 0.0, tCUDAEND = 0.0;
433 tCUDASTART += omp_get_wtime();
434 std::vector<int> tnei(npt, 0);
440 W.
getCuda().CopyMemFromDev<int, 1>(npt, devNeiPtr,
neighb.data(), 30);
443 tCUDAEND += omp_get_wtime();
450void Wake::GPUGetPairsClosestNeib(
int type)
452 size_t npt =
vtx.size();
454 double tCUDASTART = 0.0, tCUDAEND = 0.0;
456 tCUDASTART += omp_get_wtime();
457 std::vector<int> tnei(npt, 0);
464 W.
getCuda().CopyMemFromDev<int, 1>(npt, devNeiPtr,
neighb.data(), 30);
467 tCUDAEND += omp_get_wtime();
481 std::vector<bool> flag;
482 for (
int z = 0; z < times; ++z)
485#if (defined(__CUDACC__) || defined(USE_CUDA)) && (defined(CU_PAIRS))
492 double NeibStart = omp_get_wtime();
494 double NeibFinish = omp_get_wtime();
495 std::cout <<
"GPU_direct_closest_neib = " << NeibFinish - NeibStart << std::endl;
521 flag.resize(
vtx.size(),
false);
523 double sumAbsGam, iws;
526 for (
size_t vt = 0; vt + 1 <
vtx.size(); ++vt)
534 if ((ssd < 0) || (ssd >= flag.size()))
535 std::cout <<
"ssd = " << ssd <<
", flag.size() = " << flag.size() <<
", vtx.size() = " <<
vtx.size() << std::endl;
537 if ((ssd != 0) && (!flag[ssd]))
544 sumVtx.
g() = vtxI.
g() + vtxK.
g();
550 sumAbsGam = fabs(vtxI.
g()) + fabs(vtxK.
g());
552 iws = sumAbsGam > 1e-10 ? 1.0 / sumAbsGam : 1.0;
554 sumVtx.
r() = (vtxI.
r() * fabs(vtxI.
g()) + vtxK.
r() * fabs(vtxK.
g())) * iws;
558 sumVtx.
r() = (fabs(vtxI.
g()) > fabs(vtxK.
g())) ? vtxI.
r() : vtxK.
r();
601 std::vector<bool> flag;
603 for (
int z = 0; z < times; ++z)
608 flag.resize(
vtx.size(),
false);
610 double sumAbsGam, iws;
613 for (
size_t vt = 0; vt + 1 <
vtx.size(); ++vt)
619 for (
int s = 0; s < knbForRestruct; ++s)
622 int ssd =
neighbNew[vt * knbForRestruct + s];
633 sumVtx.
g() = vtxI.
g() + vtxK.
g();
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;
645 sumVtx.
r() = (fabs(vtxI.
g()) > fabs(vtxK.
g())) ? vtxI.
r() : vtxK.
r();
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)
696 std::vector<bool> flag;
698 for (
int z = 0; z < times; ++z)
703 flag.resize(
vtx.size(),
false);
705 double sumAbsGam, iws;
708 for (
size_t vt = 0; vt + 1 <
vtx.size(); ++vt)
714 for (
int s = 0; s < knbForRestruct; ++s)
717 int ssd =
neighbNew[vt * knbForRestruct + s];
728 sumVtx.
g() = vtxI.
g() + vtxK.
g();
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;
740 sumVtx.
r() = (fabs(vtxI.
g()) > fabs(vtxK.
g())) ? vtxI.
r() : vtxK.
r();
751 rindex.push_back({ (int)vt, ssd });
752 rnew.push_back(sumVtx.
r());
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)
797 if (dist2(
vtx[i].r(), zerovec) > distFar2)
810 const double porog_g = 1e-15;
814 newWake.reserve(
vtx.size());
816 for (
size_t q = 0; q <
vtx.size(); ++q)
817 if (fabs(
vtx[q].g()) > porog_g)
818 newWake.push_back(
vtx[q]);
820 size_t delta =
vtx.size() - newWake.size();
837 double timePreKnn = -omp_get_wtime();
840 std::vector<double> rightBorder, horizSpan;
861#if defined(__CUDACC__) || defined(USE_CUDA)
865 timePreKnn += omp_get_wtime();
886 for (
int collapsStep = 0; collapsStep < 1; ++collapsStep)
890 double timeResize = -omp_get_wtime();
892 timeResize += omp_get_wtime();
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();
908 timeKnn += omp_get_wtime();
911 double timeAlloc = -omp_get_wtime();
912 std::vector<std::pair<double, size_t>> initdistcuda(knbForRestruct *
vtx.size());
913 timeAlloc += omp_get_wtime();
916 double timeKnn = -omp_get_wtime();
923 knnTree.MemoryAllocate((
int)
W.
getCuda().n_CUDA_wake);
929 cudaMemcpy(&minr, knnTree.minrD,
sizeof(
Point2D), cudaMemcpyDeviceToHost);
930 cudaMemcpy(&maxr, knnTree.maxrD,
sizeof(
Point2D), cudaMemcpyDeviceToHost);
933 kNNcuda<knbForRestruct>(minr, maxr,
W.
getCuda().blocks,
vtx, initdistcuda, vecForKnn, cSP, cRBP, maxG, epsCol, collapsStep);
935 timeKnn += omp_get_wtime();
939 double timeCopy = -omp_get_wtime();
940#pragma omp parallel for
941 for (
int i = 0; i <
vtx.size(); ++i)
944 for (
int j = 0; j < knbForRestruct; ++j)
945 neighbNew[i * knbForRestruct + j] = (
int)initdist[i][j].second;
947 for (
int j = 0; j < knbForRestruct; ++j)
948 neighbNew[i * knbForRestruct + j] = (
int)initdistcuda[i * knbForRestruct + j].second;
951 timeCopy += omp_get_wtime();
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());
979 rindex.reserve(
vtx.size());
980 timeReserve += omp_get_wtime();
989 std::vector<std::pair<Point2D, Point2D>> segments(2 * ri.size());
990 for (
size_t i = 0; i < ri.size(); ++i)
992 segments[2 * i + 0].first = ri[i];
993 segments[2 * i + 0].second = rnew[i];
995 segments[2 * i + 1].first = rj[i];
996 segments[2 * i + 1].second = rnew[i];
1000 std::vector<int> hit(segments.size());
1003 double* devSegments_ptr;
1005 cudaMalloc(&devSegments_ptr, segments.size() *
sizeof(
double) * 4);
1006 cudaMemcpy(devSegments_ptr, segments.data(), segments.size() *
sizeof(
double) * 4, cudaMemcpyHostToDevice);
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();
1014 BHcu::treePanelsSegmentsIntersectionCalculationWrapper(*
W.
getNonConstCuda().auxTreePnl, cntrTreeSeg,
1019 cudaMemcpy(hit.data(),
W.
getWake().devNearestPanelPtr, segments.size() *
sizeof(
int), cudaMemcpyDeviceToHost);
1020 cudaFree(devSegments_ptr);
1023 for (
size_t i = 0; i < ri.size(); ++i)
1025 if (hit[2 * i + 0] == -1 && hit[2 * i + 1] == -1)
1027 vtx[rindex[i].first].r() = rnew[i];
1028 vtx[rindex[i].first].g() +=
vtx[rindex[i].second].g();
1030 vtx[rindex[i].second].g() = 0.0;
1037 timeCollaps = -omp_get_wtime();
1040 timeCollaps += omp_get_wtime();
1044#pragma omp parallel for private (hitA, hitB, ifInside)
1045 for (
int i = 0; i < (int)ri.size(); ++i)
1052 if ((hitA !=
size_t(-1)) || (hitB !=
size_t(-1)))
1057 vtx[rindex[i].first].r() = rnew[i];
1058 vtx[rindex[i].first].g() +=
vtx[rindex[i].second].g();
1060 vtx[rindex[i].second].g() = 0.0;
1069 double timeRemove = -omp_get_wtime();
1072 timeRemove += omp_get_wtime();
1075 W.getTimers().stop(
"Restr");
Заголовочный файл с описанием класса Airfoil.
Заголовочный файл с описанием класса Boundary.
Заголовочный файл с описанием класса MeasureVP.
Заголовочный файл с описанием класса Mechanics.
Заголовочный файл с описанием класса Preprocessor.
Заголовочный файл с описанием класса StreamParser.
Заголовочный файл с описанием класса Velocity.
Заголовочный файл с описанием класса Wake.
Заголовочный файл с описанием класса World2D.
Класс, определяющий форму профиля
const Point2D & getR(size_t q) const
Возврат константной ссылки на вершину профиля
bool inverse
Признак разворота нормалей (для расчета внутренних течений)
size_t getNumberOfPanels() const
Возврат количества панелей на профиле
Абстрактный класс, определяющий обтекаемый профиль
bool isOutsideGabarits(const Point2D &r) const
Определяет, находится ли точка с радиус-вектором вне габаритного прямоугольника профиля
Point2D upRight
Правый верхний угол габаритного прямоугольника профиля
std::vector< double > gammaThrough
Суммарные циркуляции вихрей, пересекших панели профиля на прошлом шаге
Point2D lowLeft
Левый нижний угол габаритного прямоугольника профиля
Класс, обеспечивающий возможность выполнения вычислений на GPU по технологии Nvidia CUDA.
void setCollapseCoeff(double pos_, double refLength_)
Установка правой границы самого правого профиля (для организации увеличения радиуса коллапса)
WakeDiscretizationProperties wakeDiscretizationProperties
Структура с параметрами дискретизации вихревого следа
NumericalSchemes numericalSchemes
Структура с используемыми численными схемами
std::vector< Vortex2D > vtx
Список вихревых элементов
const World2D & W
Константная ссылка на решаемую задачу
double collapseRightBorderParameter
абсцисса, правее которой происходит линейный (вправо) рост радиуса коллапса
void GetPairsBS(int type)
bool MoveInside(const Point2D &newPos, const Point2D &oldPos, const Airfoil &afl, size_t &panThrough) const
Проверка проникновения точки через границу профиля
int RemoveFar()
Зануление далеко улетевших вихрей
void GetPairsClosestNeib(int type)
double collapseScaleParameter
характерный масштаб, на котором происходит рост радиуса коллапса
std::vector< int > neighbNew
std::vector< int > neighb
Вектор потенциальных соседей для будущего коллапса
void GetPairs(int type)
Поиск ближайшего соседа
size_t RemoveZero()
Исключение нулевых и мелких вихрей
bool MoveInsideMovingBoundary(const Point2D &newPos, const Point2D &oldPos, const AirfoilGeometry &oldAfl, const Airfoil &afl, size_t &panThrough) const
Проверка проникновения точки через границу профиля
void Restruct()
Реструктуризация вихревого следа
int Collaps(int type, int times)
Коллапс вихрей
void Inside(const std::vector< Point2D > &newPos, Airfoil &afl, bool isMoves, const AirfoilGeometry &oldAfl)
Проверка пересечения вихрями следа профиля при перемещении
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)
int CollapsNew(int type, int times)
size_t getNumberOfAirfoil() const
Возврат количества профилей в задаче
VMlib::vmTimer timerMerging
const Wake & getWake() const
Возврат константной ссылки на вихревой след
const Airfoil & getAirfoil(size_t i) const
Возврат константной ссылки на объект профиля
VMlib::vmTimer timerInside
Gpu & getNonConstCuda() const
Возврат неконстантной ссылки на объект, связанный с видеокартой (GPU)
VMlib::TimersGen & getTimers() const
Возврат ссылки на временную статистику выполнения шага расчета по времени
const Gpu & getCuda() const
Возврат константной ссылки на объект, связанный с видеокартой (GPU)
const Passport & getPassport() const
Возврат константной ссылки на паспорт
Wake & getNonConstWake() const
Возврат неконстантной ссылки на вихревой след
void start(const std::string &timerLabel)
Запуск счетчика
Класс, опеделяющий двумерный вихревой элемент
HD Point2D & r()
Функция для доступа к радиус-вектору вихря
HD double & g()
Функция для доступа к циркуляции вихря
auto length2() const -> typename std::remove_const< typename std::remove_reference< decltype(this->data[0])>::type >::type
Вычисление квадрата нормы (длины) вектора
const vmTimer & stop() const
Останов работающего счетчика времени
const vmTimer & start() const
Запуск (первый или повторный) счетчика времени
const vmTimer & reset() const
Сброс счетчика времени
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)
T sqr(T x)
Умножение a на комплексно сопряженноe к b.
std::pair< std::string, int > velocityComputation
double epscol
Радиус коллапса
double maxGamma
Максимально допустимая циркуляция вихря
double distFar
Расстояние от центра самого подветренного (правого) профиля, на котором вихри уничтожаются