49#include "treeKernels.cuh"
89 if ((((aflRj - oldPos) ^ (newPos - oldPos)) * ((aflRj1 - oldPos) ^ (newPos - oldPos)) <= 0) && \
90 (((oldPos - aflRj) ^ (aflRj1 - aflRj)) * ((newPos - aflRj) ^ (aflRj1 - aflRj)) <= 0))
194 double dist2 = 1000000000.0;
199 v1 = afl.
getR(i) - newPos;
201 v2 = afl.
getR(i + 1) - newPos;
215 angle += atan2(sn, cs);
218 hit = ((angle > 3.14) || (angle < -3.14));
276 std::vector<double> gamma;
279 std::vector<int> through;
280 through.resize(
vtx.size(), -1);
283 for (
int i = 0; i < (int)
vtx.size(); ++i)
290 through[i] = (int)minN;
308 for (
size_t q = 0; q < through.size(); ++q)
311 gamma[through[q]] +=
vtx[q].g();
333#pragma omp parallel for default(none) shared(type, maxG) schedule(dynamic, DYN_SCHEDULE)
334 for (
int i = 0; i <
vtx.size(); ++i)
347 while ( (!found) && ( s + 1 < (
int)
vtx.size() ) )
352 r2 = dist2(vtxI.
r(), vtxK.
r());
355 double mnog = std::max(1.0, (vtxI.
r()[0] - cRBP) / cSP);
367 found = ( vtxI.
g()*vtxK.
g() != 0.0) && (fabs(vtxI.
g() + vtxK.
g()) < sqr(mnog) * maxG);
370 found = (vtxI.
g()*vtxK.
g() < 0.0);
373 found = (vtxI.
g()*vtxK.
g() > 0.0) && (fabs(vtxI.
g() + vtxK.
g()) < sqr(mnog) * maxG);
394#pragma omp parallel for default(none) shared(type, maxG) schedule(dynamic, DYN_SCHEDULE)
395 for (
int i = 0; i <
vtx.size(); ++i)
403 double r2, r2test = 1e+10;
408 while ( (s + 1 < (
int)
vtx.size()))
413 r2 = dist2(vtxI.
r(), vtxK.
r());
428void Wake::GPUGetPairs(
int type)
430 size_t npt =
vtx.size();
432 double tCUDASTART = 0.0, tCUDAEND = 0.0;
434 tCUDASTART += omp_get_wtime();
435 std::vector<int> tnei(npt, 0);
441 W.
getCuda().CopyMemFromDev<int, 1>(npt, devNeiPtr,
neighb.data(), 30);
444 tCUDAEND += omp_get_wtime();
451void Wake::GPUGetPairsClosestNeib(
int type)
453 size_t npt =
vtx.size();
455 double tCUDASTART = 0.0, tCUDAEND = 0.0;
457 tCUDASTART += omp_get_wtime();
458 std::vector<int> tnei(npt, 0);
465 W.
getCuda().CopyMemFromDev<int, 1>(npt, devNeiPtr,
neighb.data(), 30);
468 tCUDAEND += omp_get_wtime();
482 std::vector<bool> flag;
483 for (
int z = 0; z < times; ++z)
486#if (defined(__CUDACC__) || defined(USE_CUDA)) && (defined(CU_PAIRS))
493 double NeibStart = omp_get_wtime();
495 double NeibFinish = omp_get_wtime();
496 std::cout <<
"GPU_direct_closest_neib = " << NeibFinish - NeibStart << std::endl;
522 flag.resize(
vtx.size(),
false);
524 double sumAbsGam, iws;
527 for (
size_t vt = 0; vt + 1 <
vtx.size(); ++vt)
535 if ((ssd < 0) || (ssd >= flag.size()))
536 std::cout <<
"ssd = " << ssd <<
", flag.size() = " << flag.size() <<
", vtx.size() = " <<
vtx.size() << std::endl;
538 if ((ssd != 0) && (!flag[ssd]))
545 sumVtx.
g() = vtxI.
g() + vtxK.
g();
552 sumAbsGam = fabs(vtxI.
g()) + fabs(vtxK.
g());
554 iws = sumAbsGam > 1e-10 ? 1.0 / sumAbsGam : 1.0;
556 sumVtx.
r() = (vtxI.
r() * fabs(vtxI.
g()) + vtxK.
r() * fabs(vtxK.
g())) * iws;
560 sumVtx.
r() = (fabs(vtxI.
g()) > fabs(vtxK.
g())) ? vtxI.
r() : vtxK.
r();
603 std::vector<bool> flag;
605 for (
int z = 0; z < times; ++z)
610 flag.resize(
vtx.size(),
false);
612 double sumAbsGam, iws;
615 for (
size_t vt = 0; vt + 1 <
vtx.size(); ++vt)
621 for (
int s = 0; s < knbForRestruct; ++s)
624 int ssd =
neighbNew[vt * knbForRestruct + s];
635 sumVtx.
g() = vtxI.
g() + vtxK.
g();
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;
648 sumVtx.
r() = (fabs(vtxI.
g()) > fabs(vtxK.
g())) ? vtxI.
r() : vtxK.
r();
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)
699 std::vector<bool> flag;
701 for (
int z = 0; z < times; ++z)
706 flag.resize(
vtx.size(),
false);
708 double sumAbsGam, iws;
711 for (
size_t vt = 0; vt + 1 <
vtx.size(); ++vt)
717 for (
int s = 0; s < knbForRestruct; ++s)
720 int ssd =
neighbNew[vt * knbForRestruct + s];
731 sumVtx.
g() = vtxI.
g() + vtxK.
g();
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;
744 sumVtx.
r() = (fabs(vtxI.
g()) > fabs(vtxK.
g())) ? vtxI.
r() : vtxK.
r();
755 rindex.push_back({ (int)vt, ssd });
756 rnew.push_back(sumVtx.
r());
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)
802 if (dist2(
vtx[i].r(), zerovec) > distFar2)
815 const double porog_g = 1e-15;
819 newWake.reserve(
vtx.size());
821 for (
size_t q = 0; q <
vtx.size(); ++q)
822 if (fabs(
vtx[q].g()) > porog_g)
823 newWake.push_back(
vtx[q]);
825 size_t delta =
vtx.size() - newWake.size();
842 double timePreKnn = -omp_get_wtime();
845 std::vector<double> rightBorder, horizSpan;
866#if defined(__CUDACC__) || defined(USE_CUDA)
870 timePreKnn += omp_get_wtime();
878 bool useFastCollapseAlgorithm =
false;
879#if defined(__CUDACC__) || defined(USE_CUDA)
881 useFastCollapseAlgorithm =
true;
884 if (!useFastCollapseAlgorithm)
897 for (
int collapsStep = 0; collapsStep < 1; ++collapsStep)
901 double timeResize = -omp_get_wtime();
903 timeResize += omp_get_wtime();
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();
919 timeKnn += omp_get_wtime();
922 double timeAlloc = -omp_get_wtime();
923 std::vector<std::pair<double, size_t>> initdistcuda(knbForRestruct *
vtx.size());
924 timeAlloc += omp_get_wtime();
927 double timeKnn = -omp_get_wtime();
934 knnTree.MemoryAllocate((
int)
W.
getCuda().n_CUDA_wake);
940 cudaMemcpy(&minr, knnTree.minrD,
sizeof(
Point2D), cudaMemcpyDeviceToHost);
941 cudaMemcpy(&maxr, knnTree.maxrD,
sizeof(
Point2D), cudaMemcpyDeviceToHost);
944 kNNcuda<knbForRestruct>(minr, maxr,
W.
getCuda().blocks,
vtx, initdistcuda, vecForKnn, cSP, cRBP, maxG, epsCol, collapsStep);
946 timeKnn += omp_get_wtime();
950 double timeCopy = -omp_get_wtime();
951#pragma omp parallel for
952 for (
int i = 0; i <
vtx.size(); ++i)
955 for (
int j = 0; j < knbForRestruct; ++j)
956 neighbNew[i * knbForRestruct + j] = (
int)initdist[i][j].second;
958 for (
int j = 0; j < knbForRestruct; ++j)
959 neighbNew[i * knbForRestruct + j] = (
int)initdistcuda[i * knbForRestruct + j].second;
962 timeCopy += omp_get_wtime();
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());
990 rindex.reserve(
vtx.size());
991 timeReserve += omp_get_wtime();
1000 std::vector<std::pair<Point2D, Point2D>> segments(2 * ri.size());
1001 for (
size_t i = 0; i < ri.size(); ++i)
1003 segments[2 * i + 0].first = ri[i];
1004 segments[2 * i + 0].second = rnew[i];
1006 segments[2 * i + 1].first = rj[i];
1007 segments[2 * i + 1].second = rnew[i];
1011 std::vector<int> hit(segments.size());
1014 double* devSegments_ptr;
1016 cudaMalloc(&devSegments_ptr, segments.size() *
sizeof(
double) * 4);
1017 cudaMemcpy(devSegments_ptr, segments.data(), segments.size() *
sizeof(
double) * 4, cudaMemcpyHostToDevice);
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();
1025 BHcu::treePanelsSegmentsIntersectionCalculationWrapper(*
W.
getNonConstCuda().auxTreePnl, cntrTreeSeg,
1030 cudaMemcpy(hit.data(),
W.
getWake().devNearestPanelPtr, segments.size() *
sizeof(
int), cudaMemcpyDeviceToHost);
1031 cudaFree(devSegments_ptr);
1034 for (
size_t i = 0; i < ri.size(); ++i)
1036 if (hit[2 * i + 0] == -1 && hit[2 * i + 1] == -1)
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());
1042 vtx[rindex[i].second].g() = 0.0;
1049 timeCollaps = -omp_get_wtime();
1052 timeCollaps += omp_get_wtime();
1056#pragma omp parallel for private (hitA, hitB, ifInside)
1057 for (
int i = 0; i < (int)ri.size(); ++i)
1064 if ((hitA !=
size_t(-1)) || (hitB !=
size_t(-1)))
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());
1073 vtx[rindex[i].second].g() = 0.0;
1082 double timeRemove = -omp_get_wtime();
1085 timeRemove += omp_get_wtime();
1088 W.getTimers().stop(
"Restr");
Заголовочный файл с описанием класса Airfoil.
Заголовочный файл с описанием класса Boundary.
Заголовочный файл с функциями для метода GMRES.
Заголовочный файл с описанием класса 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 double & sigma()
Функция для доступа к радиусу вихря
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)
Возведение числа в квадрат
std::pair< std::string, int > velocityComputation
double epscol
Радиус коллапса
double maxGamma
Максимально допустимая циркуляция вихря
double distFar
Расстояние от центра самого подветренного (правого) профиля, на котором вихри уничтожаются