VM2D 1.14
Vortex methods for 2D flows simulation
Loading...
Searching...
No Matches
VM2D::World2D Class Reference

Класс, опеделяющий текущую решаемую задачу More...

#include <World2D.h>

Inheritance diagram for VM2D::World2D:
Collaboration diagram for VM2D::World2D:

Public Member Functions

Point2D getV0 () const
 Возврат текущей скорости набегающего потока
 
bool isAnyMovable () const
 Возврат признака того, что хотя бы один из профилей подвижный
 
bool isAnyMovableOrDeformable () const
 Возврат признака того, что хотя бы один из профилей подвижный или деформируемый
 
const Airfoil & getAirfoil (size_t i) const
 Возврат константной ссылки на объект профиля
 
const AirfoilGeometry & getOldAirfoil (size_t i) const
 Возврат константной ссылки на объект старого профиля
 
Airfoil & getNonConstAirfoil (size_t i) const
 Возврат неконстантной ссылки на объект профиля
 
size_t getNumberOfAirfoil () const
 Возврат количества профилей в задаче
 
const Boundary & getBoundary (size_t i) const
 Возврат константной ссылки на объект граничного условия
 
Boundary & getNonConstBoundary (size_t i) const
 Возврат неконстантной ссылки на объект граничного условия
 
size_t getNumberOfBoundary () const
 Возврат количества граничных условий в задаче
 
size_t getDispBoundaryInSystem (size_t i) const
 Возврат смещения в системе dispBoundaryInSystem.
 
const MeasureVP & getMeasureVP () const
 Возврат константной ссылки на measureVP.
 
MeasureVP & getNonConstMeasureVP () const
 Возврат неконстантной ссылки на measureVP.
 
const Mechanics & getMechanics (size_t i) const
 Возврат константной ссылки на объект механики
 
Mechanics & getNonConstMechanics (size_t i) const
 Возврат неконстантной ссылки на объект механики
 
const std::vector< std::unique_ptr< Mechanics > > & getMechanicsVector () const
 
const Wake & getWake () const
 Возврат константной ссылки на вихревой след
 
Wake & getNonConstWake () const
 Возврат неконстантной ссылки на вихревой след
 
CpuTreeInfo & getInflTreeWake () const
 
CpuTreeInfo & getCntrTreeWake () const
 
CpuTreeInfo & getAuxTreePnl () const
 
CpuTreeInfo & getCntrTreeVP () const
 
CpuTreeInfo & getCntrTreePnl () const
 
const WakeDataBase & getSource () const
 Возврат константной ссылки на источники в области течения
 
const Velocity & getVelocity () const
 Возврат константной ссылки на объект для вычисления скоростей
 
Velocity & getNonConstVelocity () const
 Возврат неконстантной ссылки на объект для вычисления скоростей
 
const Passport & getPassport () const
 Возврат константной ссылки на паспорт
 
Passport & getNonConstPassport () const
 Возврат неконстантной ссылки на паспорт
 
const Gpu & getCuda () const
 Возврат константной ссылки на объект, связанный с видеокартой (GPU)
 
Gpu & getNonConstCuda () const
 Возврат неконстантной ссылки на объект, связанный с видеокартой (GPU)
 
const std::pair< Eigen::MatrixXd, Eigen::MatrixXd > & getIQ (size_t i, size_t j) const
 Возврат константной ссылки на объект, связанный с матрицей интегралов от (r-xi)/|r-xi|^2.
 
VMlib::TimersGen & getTimers () const
 Возврат ссылки на временную статистику выполнения шага расчета по времени
 
bool ifDivisible (int val) const
 
void SolveLinearSystem ()
 Решение системы линейных алгебраических уравнений
 
void FillIQ ()
 Заполнение матрицы, состоящей из интегралов от (r-xi) / |r-xi|^2.
 
void FillMatrixAndRhs ()
 Заполнение матрицы системы для всех профилей
 
void ReserveMemoryForMatrixAndRhs ()
 Вычисляем размер матрицы и резервируем память под нее и под правую часть
 
void CalcVortexVelo ()
 Вычисление скоростей (и конвективных, и диффузионных) вихрей (в пелене и виртуальных), а также в точках вычисления VP.
 
void CalcPanelsVeloAndAttachedSheets ()
 Вычисление скоростей панелей и интенсивностей присоединенных слоев вихрей и источников
 
void CalcAndSolveLinearSystem ()
 Набор матрицы, правой части и решение СЛАУ
 
void MoveVortexes (std::vector< Point2D > &newPos, std::vector< double > *nutPtr=nullptr)
 Вычисляем новые положения вихрей (в пелене и виртуальных)
 
void CheckInside (std::vector< Point2D > &newPos, const std::vector< std::unique_ptr< AirfoilGeometry > > &oldAirfoil)
 Проверка проникновения вихрей внутрь профиля
 
void WakeAndAirfoilsMotion (bool dynamics, std::vector< double > *nutPtr=nullptr)
 Перемещение вихрей и профилей на шаге
 
 World2D (const VMlib::PassportGen &passport_)
 Конструктор
 
 ~World2D ()
 Деструктор
 
void GenerateMechanicsHeader (size_t mechanicsNumber)
 
virtual void Step () override
 Функция выполнения предварительного шага
 
VMlib::LogStream & getInfo () const
 Возврат ссылки на объект LogStream Используется в техничеcких целях для организации вывода
 
std::ostream & getInfo (char x) const
 Возврат ссылки на поток вывода информации Необходимо для вывода телеметрической информации, информации об ошибках и т.п.
 
size_t getCurrentStep () const
 Возврат константной ссылки на параметры распараллеливания по MPI.
 
double getCurrentTime () const
 
const PassportGen & getPassportGen () const
 
bool isFinished () const
 Функция, возвращающая признак завершения счета в решаемой задаче
 

Public Attributes

int gabb
 
int check01
 
int check02
 
int checkPan
 
std::unique_ptr< CpuTreeInfo > inflTreeWake
 Деревья для быстрого метода
 
std::unique_ptr< CpuTreeInfo > cntrTreeWake
 
std::unique_ptr< CpuTreeInfo > cntrTreeVP
 
std::unique_ptr< CpuTreeInfo > cntrTreePnl
 
std::unique_ptr< CpuTreeInfo > auxTreePnl
 
VMlib::vmTimer timerInitialBuild
 
VMlib::vmTimer timerRhs
 
VMlib::vmTimer timerFillMatrix
 
VMlib::vmTimer timerSlaeSolve
 
VMlib::vmTimer timerConvVelo
 
VMlib::vmTimer timerInside
 
VMlib::vmTimer timerMerging
 
size_t nVtxBeforeMerging
 

Protected Attributes

LogStream info
 Поток для вывода логов и сообщений об ошибках
 
const PassportGen & passportGen
 Константная ссылка на паспорт конкретного расчета
 
std::unique_ptr< TimersGen > timers
 Сведения о временах выполнения основных операций
 
size_t currentStep
 Текущий номер шага в решаемой задаче
 
double currentTime
 Текущее время в решаемой задаче
 

Private Attributes

std::vector< std::unique_ptr< Airfoil > > airfoil
 Список умных указателей на обтекаемые профили
 
std::vector< std::unique_ptr< AirfoilGeometry > > oldAirfoil
 Список умных указателей на обтекаемые профили для сохранения старого положения
 
std::vector< std::unique_ptr< Boundary > > boundary
 Список умных указателей на формирователи граничных условий на профилях
 
std::vector< size_t > dispBoundaryInSystem
 Список номеров, с которых начинаются элементы правой части (или матрицы) системы для профилей
 
std::vector< std::unique_ptr< Mechanics > > mechanics
 Список умных указателей на типы механической системы для каждого профиля
 
std::unique_ptr< Velocity > velocity
 Умный укзатель на объект, определяющий методику вычисления скоростей
 
std::unique_ptr< Wake > wake
 Умный указатель на вихревой след
 
std::unique_ptr< WakeDataBase > source
 Умный указатель на источники
 
std::unique_ptr< MeasureVP > measureVP
 Умный указатель на алгоритм вычисления полей скоростей и давления (для сохранения в файл)
 
Eigen::MatrixXd matrReord
 Матрица системы
 
std::vector< std::vector< std::pair< Eigen::MatrixXd, Eigen::MatrixXd > > > IQ
 Матрица, состоящая из пар матриц, в которых хранятся касательные и нормальные компоненты интегралов от ядра
 
Eigen::MatrixXd invMatr
 Обратная матрица
 
bool useInverseMatrix
 Признак использования обратной матрицы
 
Eigen::VectorXd rhsReord
 Правая часть системы
 
Eigen::VectorXd sol
 Решение системы
 
const Passport & passport
 Константная ссылка на паспорт конкретного расчета
 
Gpu cuda
 Объект, управляющий графическим ускорителем
 

Detailed Description

Класс, опеделяющий текущую решаемую задачу

Author
Марчевский Илья Константинович
Сокол Ксения Сергеевна
Рятина Евгения Павловна
Колганова Александра Олеговна
Version
1.14
Date
6 марта 2026 г.

Definition at line 76 of file World2D.h.

Constructor & Destructor Documentation

◆ World2D()

World2D::World2D ( const VMlib::PassportGen &  passport_)

Конструктор

Parameters
[in]passport_константная ссылка на паспорт расчета

Definition at line 73 of file World2D.cpp.

73 :
74 WorldGen(passport_),
75 passport(dynamic_cast<const Passport&>(passport_)),
76 cuda(Gpu(*this))
77{
78 std::stringstream ss;
79 ss << "#" << passport.problemNumber << " (" << passport.problemName << ")";
81
83 currentStep = 0;
84
85 std::vector<std::string> timerLabels = { "Step", "MatRhs", "Solve", "ConvVel", "Nut", /*"Nbody",*/ "DiffVel", "Force", "VelPres", "Inside", "Restr", "Save"};
86 timers = std::make_unique<VMlib::TimersGen>(*this, timerLabels);
87
88 wake.reset(new Wake(*this));
89 // загрузка пелены из файла
91 {
93 wake->ReadFromFile(passport.wakesDir, passport.wakeDiscretizationProperties.fileWake); //Считываем из каталога с пеленой
94 else
95 {
97 wake->ReadFromFile(passport.dir, passport.wakeDiscretizationProperties.fileWake); //Считываем из каталога с пеленой
98 }
99
100
101 }
102
103 source.reset(new WakeDataBase(*this));
104 // загрузка положений источников из файла
106 source->ReadFromFile(passport.dir, passport.wakeDiscretizationProperties.fileSource); //Считываем из текущего каталога
107
108
110 {
111 case 0:
112 velocity.reset(new VelocityBiotSavart(*this));
113 break;
114 case 1:
115 velocity.reset(new VelocityBarnesHut(*this));
121 else
123 break;
124 }
125
126 velocity->virtualVortexesParams.resize(passport.airfoilParams.size());
127
128 auto CreateBoundary = [this](size_t i) {
130 {
131 case 0:
132 this->boundary.emplace_back(new BoundaryVortexCollocN(*this, i));
133 //info('e') << "BoundaryMDV is not implemented now! " << std::endl;
134 //exit(1);
135 break;
136
137 case 1:
138 this->boundary.emplace_back(new BoundaryConstLayerAver(*this, i));
139 break;
140
141 case 2:
142 this->boundary.emplace_back(new BoundaryLinLayerAver(*this, i));
143 break;
144
145 default:
146 info('e') << "Unknown scheme!" << std::endl;
147 exit(1);
148 }
149 };
150
151
152 for (size_t i = 0; i < passport.airfoilParams.size(); ++i)
153 {
154
155 switch (passport.airfoilParams[i].mechanicalSystemType)
156 {
157 case 0:
158 airfoil.emplace_back(new AirfoilRigid(*this, i));
159 airfoil[i]->ReadFromFile(passport.airfoilsDir);
160 CreateBoundary(i);
161 mechanics.emplace_back(new MechanicsRigidImmovable(*this, i));
162 break;
163
164 case 1:
165 airfoil.emplace_back(new AirfoilRigid(*this, i));
166 airfoil[i]->ReadFromFile(passport.airfoilsDir);
167 CreateBoundary(i);
168 mechanics.emplace_back(new MechanicsRigidGivenLaw(*this, i));
169 break;
170
171 case 2:
172 airfoil.emplace_back(new AirfoilRigid(*this, i));
173 airfoil[i]->ReadFromFile(passport.airfoilsDir);
174 CreateBoundary(i);
175 mechanics.emplace_back(new MechanicsRigidOscillPart(*this, i));
176 break;
177
178 case 3:
179 airfoil.emplace_back(new AirfoilRigid(*this, i));
180 airfoil[i]->ReadFromFile(passport.airfoilsDir);
181 CreateBoundary(i);
182 mechanics.emplace_back(new MechanicsRigidRotatePart(*this, i));
183 break;
184
185 case 4:
186 airfoil.emplace_back(new AirfoilDeformable(*this, i));
187 airfoil[i]->ReadFromFile(passport.airfoilsDir);
188 CreateBoundary(i);
189 mechanics.emplace_back(new MechanicsDeformable(*this, i));
190 break;
191 }
192 }
193
194 if (getPassport().wakeDiscretizationProperties.sigma0 == 0)
195 {
196 double sumLength = 0.0;
197 size_t totNumPan = 0;
198 for (size_t bou = 0; bou < getNumberOfAirfoil(); ++bou)
199 {
200 for (size_t pnl = 0; pnl < getAirfoil(bou).getNumberOfPanels(); ++pnl)
201 sumLength += getAirfoil(bou).len[pnl];
202 totNumPan += getAirfoil(bou).getNumberOfPanels();
203 }
204 double eps = 0.5 * sumLength / totNumPan;
206 info('i') << "sigma0 = " << getPassport().wakeDiscretizationProperties.sigma0 << " is calculated automatically" << std::endl;
207 }
208
209 if (getPassport().wakeDiscretizationProperties.epscol == 0)
210 {
212
213 info('i') << "epscol = " << getPassport().wakeDiscretizationProperties.epscol << " is calculated automatically" << std::endl;
214 }
215
216 for (size_t afl = 0; afl < passport.airfoilParams.size(); ++afl)
217 if (passport.airfoilParams[afl].chord == 0)
218 {
219 auto& prm = getNonConstPassport().airfoilParams[afl];
220 prm.chord = (prm.initialGab.second[0] - prm.initialGab.first[0]) * prm.scale[0];
221 info('i') << "airfoil #" << afl << " chord = " << prm.chord << " is calculated automatically" << std::endl;
222 }
223
224
225#ifdef USE_CUDA
226 if (getNumberOfAirfoil() > 0)
228 cuda.Gmres.reset(new GmresSolver(*this));
229#endif
230
231 //2026-03-28
232 //считываем массив точек для подсчета и вывода поля скоростей и давлений
233 measureVP.reset(new MeasureVP(*this));
234
236 measureVP->ReadPointsFromFile(passport.dir);
237
238#ifdef OPTIMIZER
239// std::vector<Point2D> VPPoints, VPHistory;
240// for (size_t a = 0; a < getNumberOfAirfoil(); ++a)
241// {
242// const auto& afl = *airfoil[a];
243// for (size_t p = 0; p < afl.getNumberOfPanels(); ++p)
244// VPHistory.push_back(0.5 * (afl.getR(p) + afl.getR(p + 1)) + afl.nrm[p] * passport.airfoilParams[a].chord * 0.01);
245// }
246// if (passport.timeDiscretizationProperties.saveVPstep != 0)
247// measureVP->SetPoints(VPPoints, VPHistory);
248#endif
249
250 IQ.resize(passport.airfoilParams.size());
251 for (size_t i = 0; i < passport.airfoilParams.size(); ++i)
252 IQ[i].resize(passport.airfoilParams.size());
253
254
255 info.endl();
256 info.endl();
257 info('i') << "Start solving problem " << passport.problemName << std::endl;
258 info.endl();
259
260 VMlib::PrintLogoToStream(info('_') << std::endl);
261}//World2D(...)
Класс, определяющий тип обтекаемого профиля
std::vector< double > len
Длины панелей профиля
Definition Airfoil2D.h:94
size_t getNumberOfPanels() const
Возврат количества панелей на профиле
Definition Airfoil2D.h:163
Класс, определяющий тип обтекаемого профиля
Класс, определяющий способ удовлетворения граничного условия на обтекаемом профиле
Класс, определяющий способ удовлетворения граничного условия на обтекаемом профиле
Класс, определяющий способ удовлетворения граничного условия на обтекаемом профиле
Структура, хранящая данные и указатели на массивы на CPU для оптимизации итерационного решения СЛАУ н...
Definition cpuTreeInfo.h:96
Класс, обеспечивающий возможность выполнения вычислений на GPU по технологии Nvidia CUDA.
Definition Gpu2D.h:71
Класс, отвечающий за вычисление поля скорости и давления в заданых точках для вывода
Definition MeasureVP2D.h:65
Класс, определяющий вид механической системы
Класс, определяющий вид механической системы
Класс, определяющий вид механической системы
Класс, определяющий вид механической системы
Класс, определяющий вид механической системы
Класс, опеделяющий паспорт двумерной задачи
Definition Passport2D.h:256
std::string wakesDir
Каталог с файлами вихревых следов
Definition Passport2D.h:276
WakeDiscretizationProperties wakeDiscretizationProperties
Структура с параметрами дискретизации вихревого следа
Definition Passport2D.h:307
std::string airfoilsDir
Каталог с файлами профилей
Definition Passport2D.h:273
std::vector< AirfoilParams > airfoilParams
Список структур с параметрами профилей
Definition Passport2D.h:279
NumericalSchemes numericalSchemes
Структура с используемыми численными схемами
Definition Passport2D.h:310
Класс, определяющий способ вычисления скоростей
Класс, определяющий способ вычисления скоростей
Класс, опеделяющий набор вихрей
Класс, опеделяющий вихревой след (пелену)
Definition Wake2D.h:63
size_t getNumberOfAirfoil() const
Возврат количества профилей в задаче
Definition World2D.h:180
std::unique_ptr< CpuTreeInfo > inflTreeWake
Деревья для быстрого метода
Definition World2D.h:113
Gpu cuda
Объект, управляющий графическим ускорителем
Definition World2D.h:143
std::vector< std::unique_ptr< Airfoil > > airfoil
Список умных указателей на обтекаемые профили
Definition World2D.h:80
std::unique_ptr< MeasureVP > measureVP
Умный указатель на алгоритм вычисления полей скоростей и давления (для сохранения в файл)
Definition World2D.h:104
const Airfoil & getAirfoil(size_t i) const
Возврат константной ссылки на объект профиля
Definition World2D.h:163
std::unique_ptr< WakeDataBase > source
Умный указатель на источники
Definition World2D.h:101
std::vector< std::unique_ptr< Mechanics > > mechanics
Список умных указателей на типы механической системы для каждого профиля
Definition World2D.h:92
std::vector< std::unique_ptr< Boundary > > boundary
Список умных указателей на формирователи граничных условий на профилях
Definition World2D.h:86
const Passport & getPassport() const
Возврат константной ссылки на паспорт
Definition World2D.h:263
std::unique_ptr< CpuTreeInfo > cntrTreePnl
Definition World2D.h:116
Passport & getNonConstPassport() const
Возврат неконстантной ссылки на паспорт
Definition World2D.h:268
std::unique_ptr< CpuTreeInfo > cntrTreeVP
Definition World2D.h:115
std::unique_ptr< Wake > wake
Умный указатель на вихревой след
Definition World2D.h:98
std::unique_ptr< Velocity > velocity
Умный укзатель на объект, определяющий методику вычисления скоростей
Definition World2D.h:95
const Passport & passport
Константная ссылка на паспорт конкретного расчета
Definition World2D.h:140
std::vector< std::vector< std::pair< Eigen::MatrixXd, Eigen::MatrixXd > > > IQ
Матрица, состоящая из пар матриц, в которых хранятся касательные и нормальные компоненты интегралов о...
Definition World2D.h:125
std::unique_ptr< CpuTreeInfo > cntrTreeWake
Definition World2D.h:114
void endl()
Вывод в поток логов пустой строки
Definition LogStream.h:103
void assignStream(std::ostream *pStr_, const std::string &label_)
Связывание потока логов с потоком вывода
Definition LogStream.h:80
TimeDiscretizationProperties timeDiscretizationProperties
Структура с параметрами процесса интегрирования по времени
std::string problemName
Название задачи
std::string dir
Рабочий каталог задачи
size_t problemNumber
Номер задачи
std::unique_ptr< TimersGen > timers
Сведения о временах выполнения основных операций
Definition WorldGen.h:66
double currentTime
Текущее время в решаемой задаче
Definition WorldGen.h:72
VMlib::LogStream & getInfo() const
Возврат ссылки на объект LogStream Используется в техничеcких целях для организации вывода
Definition WorldGen.h:82
size_t currentStep
Текущий номер шага в решаемой задаче
Definition WorldGen.h:69
WorldGen(const VMlib::PassportGen &passport_)
Конструктор
Definition WorldGen.cpp:47
LogStream info
Поток для вывода логов и сообщений об ошибках
Definition WorldGen.h:60
bool fileExistTest(const std::string &fileName_, LogStream &info, bool exitKey=false, const std::list< std::string > &extList={})
Проверка существования файла
Definition defs.h:343
void PrintLogoToStream(std::ostream &str)
Передача в поток вывода шапки программы VM2D/VM3D.
Definition defs.cpp:108
static std::ostream * defaultWorld2DLogStream
Поток вывода логов и ошибок задачи
Definition defs.h:225
std::pair< std::string, int > boundaryCondition
Метод аппроксимации граничных условий
Definition Passport2D.h:193
std::pair< std::string, int > velocityComputation
Definition Passport2D.h:183
std::pair< std::string, int > linearSystemSolver
Definition Passport2D.h:177
std::string fileSource
Имя файла с положениями источников (без полного пути)
Definition Passport2D.h:150
double epscol
Радиус коллапса
Definition Passport2D.h:132
std::string fileWake
Имя файла с начальным состоянием вихревого следа (без полного пути)
Definition Passport2D.h:147
double sigma0
Радиус вихря
Definition Passport2D.h:126
double timeStart
Начальное время
Definition PassportGen.h:61
int saveVPstep
Шаг вычисления и сохранения скорости и давления
Definition PassportGen.h:86
Here is the call graph for this function:

◆ ~World2D()

VM2D::World2D::~World2D ( )
inline

Деструктор

Definition at line 354 of file World2D.h.

354{};

Member Function Documentation

◆ CalcAndSolveLinearSystem()

void World2D::CalcAndSolveLinearSystem ( )

Набор матрицы, правой части и решение СЛАУ

Вызывается в Step()

Definition at line 2108 of file World2D.cpp.

2109{
2110 if (airfoil.size() > 0)
2111 {
2113
2114#if defined(__CUDACC__) || defined(USE_CUDA)
2117
2119
2120 if ((sch == 1) || (sch == 2) || (sch == 0))
2122 else
2123 {
2124 info('e') << "schemeSwitcher is not 0, or 1, or 2! " << std::endl;
2125 exit(1);
2126 }
2127
2128 cuda.RefreshWake(1);
2129 cuda.RefreshAfls(1);
2130 cuda.RefreshVirtualWakes(1);
2131
2132#endif
2133 const int linSystemScheme = passport.numericalSchemes.linearSystemSolver.second;
2134
2135 if (linSystemScheme == 0) //Gauss
2136 {
2137 if (currentStep == 0)
2138 {
2139#ifdef BRIDGE
2140 useInverseMatrix = true;
2141#endif //BRIDGE
2142
2143#ifdef INITIAL
2145 ((mechanics.size() == 1) && (!mechanics.front()->isDeform))
2146 ||
2147 (mechanics.size() > 1 && !isAnyMovableOrDeformable())
2148 );
2149#endif //INITIAL
2150 }
2151 }//if Gauss
2152
2153 if (linSystemScheme == 0 || linSystemScheme == 1 || (linSystemScheme == 2 && isAnyMovableOrDeformable()))
2154 FillIQ();
2155
2157
2158 //{
2159 // std::stringstream ss;
2160 // ss << "IQ1-" << currentStep;
2161 // std::ofstream of(passport.dir + "dbg/" + ss.str());
2162 // for (size_t i = 0; i < (IQ[0][0]).first.rows(); ++i)
2163 // {
2164 // for (size_t j = 0; j < (IQ[0][0]).first.cols(); ++j)
2165 // of << (IQ[0][0]).first(i, j) << " ";
2166 // of << std::endl;
2167 // }
2168 // of.close();
2169 //}
2170
2171 /*
2172 {
2173 std::stringstream ss;
2174 ss << "matr-" << currentStep;
2175 std::ofstream of(passport.dir + "dbg/" + ss.str());
2176 of.precision(16);
2177 for (size_t i = 0; i < matrReord.rows(); ++i)
2178 {
2179 for (size_t j = 0; j < matrReord.cols(); ++j)
2180 of << matrReord(i, j) << " ";
2181 of << std::endl;
2182 }
2183 of.close();
2184 }
2185
2186 {
2187 std::stringstream ss;
2188 ss << "rhs-" << currentStep;
2189 std::ofstream of(passport.dir + "dbg/" + ss.str());
2190 of.precision(16);
2191 for (size_t i = 0; i < rhsReord.rows(); ++i)
2192 {
2193 of << rhsReord(i) << std::endl;
2194 }
2195 of.close();
2196 }
2197
2198 //*/
2199
2201
2202 size_t currentRow = 0;
2203 for (size_t bou = 0; bou < boundary.size(); ++bou)
2204 {
2205 size_t nVars = boundary[bou]->GetUnknownsSize();
2206 Eigen::VectorXd locSol;
2207 locSol.resize(nVars);
2208 for (size_t i = 0; i < nVars; ++i)
2209 locSol(i) = sol(currentRow + i);
2210
2211 boundary[bou]->SolutionToFreeVortexSheetAndVirtualVortex(locSol);
2212 currentRow += nVars;// +1;
2213 }
2214
2215 //14-05-2024
2216 size_t nVirtVortices = 0;
2217 for (size_t bou = 0; bou < boundary.size(); ++bou)
2218 nVirtVortices += boundary[bou]->virtualWake.vtx.size();
2219
2220 wake->vtx.reserve(wake->vtx.size() + nVirtVortices);
2221 for (size_t bou = 0; bou < boundary.size(); ++bou)
2222 for (size_t v = 0; v < boundary[bou]->virtualWake.vtx.size(); ++v)
2223 wake->vtx.push_back(Vortex2D{ boundary[bou]->virtualWake.vtx[v].r(), 0.0, getPassport().wakeDiscretizationProperties.sigma0});
2224
2225 //for (size_t v = 0; v < wake->vtx.size(); ++v)
2226 //{
2227 // if (fabs(wake->vtx[v].g()) > 1.0)
2228 // std::cout << "Error gam! " << "v = " << v << ", gam[v] = " << wake->vtx[v].g() << std::endl;
2229 //}
2230
2231 //for (size_t bou = 0; bou < boundary.size(); ++bou)
2232 // for (size_t v = 0; v < airfoil[bou]->getNumberOfPanels(); ++v)
2233 // {
2234 // if (fabs(boundary[bou]->sheets.freeVortexSheet(v, 0)) > 1000.0)
2235 // std::cout << "Error gamma sheet! " << "bou = " << bou << ", v = " << v << ", gam[v] = " << boundary[bou]->sheets.freeVortexSheet(v, 0) << std::endl;
2236 // }
2237
2238 /*
2239 if (currentStep == 0)
2240 {
2241 Point2D addMass = { 0.0, 0.0 };
2242 for (size_t q = 0; q < airfoil[0]->getNumberOfPanels(); ++q)
2243 {
2244 addMass += (sol(q) + boundary[0]->sheets.attachedVortexSheet(q,0)) * 0.5 * (airfoil[0]->getR(q) + airfoil[0]->getR(q + 1)).kcross() * airfoil[0]->len[q];
2245 }
2246 addMass *= passport.physicalProperties.rho;
2247
2248 std::cout << "AddMass = " << addMass << std::endl;
2249 //exit(-42);
2250 }
2251 */
2252 }
2253}//CalcAndSolveLinearSystem()
void setAccelCoeff(double cft_)
Установка коэффициента разгона потока
Definition Gpu2D.h:281
void setMaxGamma(double gam_)
Установка максимально допустимой циркуляции вихря
Definition Gpu2D.h:304
void setSchemeSwitcher(int schemeSwitcher_)
Установка переключателя расчетных схем
Definition Gpu2D.h:318
PhysicalProperties physicalProperties
Структура с физическими свойствами задачи
Definition Passport2D.h:304
void ReserveMemoryForMatrixAndRhs()
Вычисляем размер матрицы и резервируем память под нее и под правую часть
Definition World2D.cpp:1396
bool isAnyMovableOrDeformable() const
Возврат признака того, что хотя бы один из профилей подвижный или деформируемый
Definition World2D.cpp:2370
bool useInverseMatrix
Признак использования обратной матрицы
Definition World2D.h:131
void FillIQ()
Заполнение матрицы, состоящей из интегралов от (r-xi) / |r-xi|^2.
Definition World2D.cpp:2066
void SolveLinearSystem()
Решение системы линейных алгебраических уравнений
Definition World2D.cpp:1001
void FillMatrixAndRhs()
Заполнение матрицы системы для всех профилей
Definition World2D.cpp:1298
Eigen::VectorXd sol
Решение системы
Definition World2D.h:137
Класс, опеделяющий двумерный вихревой элемент
Definition Vortex2D.h:59
double getCurrentTime() const
Definition WorldGen.h:100
double accelCft(double currentTime) const
Функция-множитель, позволяющая моделировать разгон
double maxGamma
Максимально допустимая циркуляция вихря
Definition Passport2D.h:144
Here is the call graph for this function:
Here is the caller graph for this function:

◆ CalcPanelsVeloAndAttachedSheets()

void World2D::CalcPanelsVeloAndAttachedSheets ( )

Вычисление скоростей панелей и интенсивностей присоединенных слоев вихрей и источников

Вызывается в Step()

Definition at line 1637 of file World2D.cpp.

1638{
1639 //вычисляем скорости панелей
1640 for (size_t i = 0; i < airfoil.size(); ++i)
1642
1643 //вычисляем интенсивности присоединенных слоев
1644 for (size_t i = 0; i < airfoil.size(); ++i)
1645 boundary[i]->ComputeAttachedSheetsIntensity();
1646
1647}//CalcPanelsVeloAndAttachedSheets(...)
double dt
Шаг по времени
Definition PassportGen.h:67
Here is the caller graph for this function:

◆ CalcVortexVelo()

void World2D::CalcVortexVelo ( )

Вычисление скоростей (и конвективных, и диффузионных) вихрей (в пелене и виртуальных), а также в точках вычисления VP.

Вызывается в Step()

Definition at line 1456 of file World2D.cpp.

1457{
1458 velocity->ResizeAndZero();
1459
1460 getTimers().start("ConvVel");
1461
1462 //Конвективные скорости всех вихрей (и в пелене, и виртуальных), индуцируемые вихрями в пелене
1463 //Подготовка CUDA
1464#if defined(__CUDACC__) || defined(USE_CUDA)
1465 cuda.RefreshWake(2);
1466
1467 cuda.RefreshAfls(2);
1468 cuda.RefreshVirtualWakes(2);
1469#endif
1470 getTimers().stop("ConvVel");
1471
1472 velocity->CalcConvVelo();
1473
1474 //Расчет средних значений eps для каждой панели и их передача на видеокарту
1475
1476 getTimers().start("DiffVel");
1477
1478 for (size_t bou = 0; bou < getNumberOfBoundary(); ++bou)
1480
1481#if defined(__CUDACC__) || defined(USE_CUDA)
1482 for (size_t i = 0; i < airfoil.size(); ++i)
1483 cuda.CopyMemToDev<double, 1>(airfoil[i]->getNumberOfPanels(), airfoil[i]->meanEpsOverPanel.data(), airfoil[i]->devMeanEpsOverPanelPtr);
1484#endif
1485
1486 getTimers().stop("DiffVel");
1487
1488 //Вычисление диффузионных скоростей вихрей (и в пелене, и виртуальных)
1489 velocity->CalcDiffVelo();
1490
1491 //Обнуление вязких напряжений, если они не были вычислены
1492 for (auto& afl : airfoil)
1493 {
1494 if (afl->viscousStress.size() == 0)
1495 afl->viscousStress.resize(afl->getNumberOfPanels(), 0.0);
1496 }
1497
1498 getTimers().start("Save");
1499
1500 /*
1501 //Сохранение всех параметров для вихрей в пелене
1502 if (!(currentStep % 1))
1503 {
1504 VMlib::CreateUserDirectory(passport.dir, "dbg");
1505 std::ostringstream sss, sss2;
1506 sss << "prmWake";
1507 sss << currentStep;
1508 sss2 << "testWake";
1509 sss2 << currentStep;
1510
1511 std::ofstream prmtFile(passport.dir + "dbg/" + sss.str());
1512 prmtFile.precision(11);
1513 prmtFile << "i x y g sigma epsast convVeloX convVeloY diffVeloX diffVeloY I0 I1 I2X I2Y I3X I3Y" << std::endl;
1514 for (size_t i = 0; i < wake->vtx.size(); ++i)
1515 prmtFile << i << " " \
1516 << wake->vtx[i].r()[0] << " " << wake->vtx[i].r()[1] << " " \
1517 << wake->vtx[i].g() << " " \
1518 << wake->vtx[i].sigma() << " " \
1519 << velocity->wakeVortexesParams.epsastWake[i] << " " \
1520 << velocity->wakeVortexesParams.convVelo[i][0] << " " << velocity->wakeVortexesParams.convVelo[i][1] << " "\
1521 << velocity->wakeVortexesParams.diffVelo[i][0] << " " << velocity->wakeVortexesParams.diffVelo[i][1] << " "\
1522 //<< std::endl;
1523 << velocity->wakeVortexesParams.I0[i] << " " \
1524 << velocity->wakeVortexesParams.I1[i] << " " \
1525 << velocity->wakeVortexesParams.I2[i][0] << " " << velocity->wakeVortexesParams.I2[i][1] << " " \
1526 << velocity->wakeVortexesParams.I3[i][0] << " " << velocity->wakeVortexesParams.I3[i][1] << " " \
1527 << "\n";
1528
1529 prmtFile.close();
1530
1531 std::ofstream prmtFile2(passport.dir + "dbg/" + sss2.str());
1532 prmtFile2 << wake->vtx.size() << "\n";
1533 prmtFile2.precision(16);
1534 for (size_t i = 0; i < wake->vtx.size(); ++i)
1535 prmtFile2 << wake->vtx[i].r()[0] << " " << wake->vtx[i].r()[1] << " " << wake->vtx[i].g() << "\n";
1536 prmtFile2.close();
1537 }
1538//*/
1539
1540/*
1541 //Сохранение всех параметров для виртуальных вихрей
1542 if (!(currentStep % 1))
1543 {
1544 for (size_t b = 0; b < boundary.size(); ++b)
1545 {
1546 std::ostringstream sss;
1547 sss << "prmVirtual_";
1548 sss << b << "-";
1549 sss << currentStep;
1550 std::ofstream prmFileVirt(passport.dir + "dbg/" + sss.str());
1551 //prmFileVirt.precision(16);
1552 prmFileVirt << "i x y g sigma epsast convVeloX convVeloY diffVeloX diffVeloY I0 I1 I2X I2Y I3X I3Y" << std::endl;
1553 for (size_t i = 0; i < boundary[b]->virtualWake.vtx.size(); ++i)
1554 prmFileVirt << i << " " \
1555 << boundary[b]->virtualWake.vtx[i].r()[0] << " " << boundary[b]->virtualWake.vtx[i].r()[1] << " " \
1556 << boundary[b]->virtualWake.vtx[i].g() << " " \
1557 << boundary[b]->virtualWake.vtx[i].sigma() << " " \
1558 << velocity->virtualVortexesParams[b].epsastWake[i] << " " \
1559 << velocity->virtualVortexesParams[b].convVelo[i][0] << " " << velocity->virtualVortexesParams[b].convVelo[i][1] << " "\
1560 << velocity->virtualVortexesParams[b].diffVelo[i][0] << " " << velocity->virtualVortexesParams[b].diffVelo[i][1] << " "\
1561 //<< std::endl;
1562 << velocity->virtualVortexesParams[b].I0[i] << " " \
1563 << velocity->virtualVortexesParams[b].I1[i] << " " \
1564 << velocity->virtualVortexesParams[b].I2[i][0] << " " << velocity->virtualVortexesParams[b].I2[i][1] << " " \
1565 << velocity->virtualVortexesParams[b].I3[i][0] << " " << velocity->virtualVortexesParams[b].I3[i][1] << " " \
1566 << std::endl;
1567 prmFileVirt.close();
1568 }
1569 //if (currentStep==3) exit(-123);
1570 }
1571//*/
1572
1573 getTimers().stop("Save");
1574
1575}//CalcVortexVelo()
void calcMeanEpsOverPanel()
Вычисление средних значений eps на панелях
Definition Airfoil2D.cpp:96
VMlib::TimersGen & getTimers() const
Возврат ссылки на временную статистику выполнения шага расчета по времени
Definition World2D.h:288
size_t getNumberOfBoundary() const
Возврат количества граничных условий в задаче
Definition World2D.h:197
Airfoil & getNonConstAirfoil(size_t i) const
Возврат неконстантной ссылки на объект профиля
Definition World2D.h:175
void stop(const std::string &timerLabel)
Останов счетчика
Definition TimesGen.cpp:68
void start(const std::string &timerLabel)
Запуск счетчика
Definition TimesGen.cpp:55
Here is the call graph for this function:
Here is the caller graph for this function:

◆ CheckInside()

void World2D::CheckInside ( std::vector< Point2D > &  newPos,
const std::vector< std::unique_ptr< AirfoilGeometry > > &  oldAirfoil 
)

Проверка проникновения вихрей внутрь профиля

Вызывается в Step()

Parameters
[in]newPosновые позиции вихрей
[in]oldAirfoilконстантная ссылка на вектор из умных указателей на старые положения профилей

Definition at line 813 of file World2D.cpp.

814{
815 getTimers().start("Inside");
816
817 gabb = 0;
818 check01 = 0;
819 check02 = 0;
820 checkPan = 0;
821
822#if (!defined(USE_CUDA))
823 nVtxBeforeMerging = newPos.size();
824 for (size_t afl = 0; afl < airfoil.size(); ++afl)
825 wake->Inside(newPos, *airfoil[afl], mechanics[afl]->isMoves, *oldAirfoil[afl]);
826
827
828 //Тренируемся работать с древом панелей
829 /*
830 int nTotPan = 0;
831 for (size_t s = 0; s < getNumberOfAirfoil(); ++s)
832 nTotPan += (int)getAirfoil(s).getNumberOfPanels();
833
834 auxTreePnl.reset(new CpuTreeInfo(tree_T::aux, object_T::panel, scheme_T::noScheme));
835
836 std::vector<std::pair<Point2D, Point2D>> pnls(nTotPan);
837 int counter = 0;
838 for (size_t s = 0; s < getNumberOfAirfoil(); ++s)
839 for (size_t p = 0; p < getAirfoil(s).getNumberOfPanels(); ++p)
840 pnls[counter++] = { getAirfoil(s).getR(p), getAirfoil(s).getR(p + 1) };
841
842 auxTreePnl->UpdatePanelGeometry(pnls, 0);
843 auxTreePnl->Build();
844 auxTreePnl->UpwardTraversal(0);
845
846 std::unique_ptr<CpuTreeInfo> cntrTreePnt;
847 cntrTreePnt.reset(new CpuTreeInfo(tree_T::contr, object_T::point2, scheme_T::noScheme));
848
849 std::vector<Vortex2D> newPosVtx(newPos.size());
850 for (size_t i = 0; i < newPos.size(); ++i)
851 newPosVtx[i].r() = newPos[i];
852
853 cntrTreePnt->Update(newPosVtx);
854 float tBuild = cntrTreePnt->Build();
855 float tUpward = cntrTreePnt->UpwardTraversal(0);
856
857 std::vector<std::pair<int,double>> closeIndexDist(cntrTreePnt->object.size());
858
859
860 VMlib::vmTimer timerA;
861 timerA.start();
862 double tDownward = auxTreePnl->DownwardTraversalClosestPanelToPoints(*cntrTreePnt, closeIndexDist, false, nullptr);
863 timerA.stop();
864
865 std::ofstream treeTimeFile;
866 if (getCurrentStep() == 0)
867 {
868 treeTimeFile.open(passport.dir + "/dbg/treeTime.csv");
869 treeTimeFile << "step,time,N,tBUI,tUPW,tDNW\n";
870 }
871 else
872 treeTimeFile.open(passport.dir + "/dbg/treeTime.csv", std::ios::app);
873
874 treeTimeFile << getCurrentStep() << ',' << getCurrentTime() << ',' << newPos.size() << ',' << tBuild << ',' << tUpward << ',' << tDownward << '\n';
875
876 treeTimeFile.close();
877
878 //std::cout << "timeFastNeib = " << timerA.duration() << ", N = " << newPos.size() << "\n";
879 /*
880 std::ofstream closeFile(passport.dir + "closeFile.txt");
881 for (size_t i = 0; i < newPos.size(); ++i)
882 closeFile << newPos[i][0] << " " << newPos[i][1] << " " << closeIndexDist[i].first << "\n";
883 closeFile.close();
884 */
885
886 //*/
887
888 //std::cout << "gab: " << gabb << " check1: " << check01 << " check2: " << check02 << " checkPan: " << checkPan << std::endl;
889#else
890// //////////////////////// CUDA ///////////////////////
893
894 if ((newPos.size() > 0) && (getNumberOfAirfoil() > 0))
895 {
896 int nTotPanels = 0;
897 for (size_t afl = 0; afl < airfoil.size(); ++afl)
898 nTotPanels += (int)airfoil[afl]->getNumberOfPanels();
899
900 nVtxBeforeMerging = newPos.size();
901
902 auto& auxTree = *getNonConstCuda().auxTreePnl;
903
905 {
906 auxTree.MemoryAllocate((int)getCuda().n_CUDA_pnls);
907 auxTree.UpdatePanelGeometry(nTotPanels, (double4*)airfoil[0]->devRPtr);
908 auxTree.Build();
909 auxTree.UpwardTraversal(0);
910 }
911
912 std::vector<int> hit(newPos.size());
913
914 if (isAnyMovableOrDeformable()) //Если профили подвижны - метод псевдонормалей
915 {
916 double* devNewpos_ptr;
917 cudaMalloc(&devNewpos_ptr, newPos.size() * sizeof(double) * 2);
918 cudaMemcpy(devNewpos_ptr, newPos.data(), newPos.size() * sizeof(double) * 2, cudaMemcpyHostToDevice);
919
920 auto& cntrTreePnt = *getCuda().cntrTreePoint;
921 cntrTreePnt.MemoryAllocate((int)getCuda().n_CUDA_wake);
922 cntrTreePnt.Update((int)newPos.size(), devNewpos_ptr);
923 cntrTreePnt.Build();
924
925 BHcu::treeClosestPanelToPointsCalculationWrapper(auxTree, cntrTreePnt, getWake().devNearestPanelPtr, true, getAirfoil(0).devPsnPtr);
926
927 cudaMemcpy(hit.data(), getWake().devNearestPanelPtr, newPos.size() * sizeof(int), cudaMemcpyDeviceToHost);
928 cudaFree(devNewpos_ptr);
929
930
931 //std::ofstream vortexFile(getPassport().dir + "vortexFile.txt");
932 //for (size_t i = 0; i < getWake().vtx.size(); ++i)
933 // vortexFile << getWake().vtx[i].r()[0] << " " << getWake().vtx[i].r()[1] << " " << hit[i] << '\n';
934 //vortexFile.close();
935
936 //std::ofstream panelFile(getPassport().dir + "panelFile.txt");
937 //for (size_t i = 0; i < getAirfoil(0).getNumberOfPanels(); ++i)
938 // panelFile << getAirfoil(0).getR(i)[0] << " " << getAirfoil(0).getR(i)[1] << " " << \
939 // getAirfoil(0).getR(i + 1)[0] << " " << getAirfoil(0).getR(i + 1)[1] << " " << '\n';
940 //panelFile.close();
941
942
943 }
944 else //иначе (для неподвижного профиля - метод трассировки лучей
945 {
946 std::vector<std::pair<Point2D, Point2D>> segments(newPos.size());
947 for (size_t i = 0; i < newPos.size(); ++i)
948 {
949 segments[i].first = wake->vtx[i].r();
950 segments[i].second = newPos[i];
951 }
952 double* devSegments_ptr;
953
954 cudaMalloc(&devSegments_ptr, newPos.size() * sizeof(double) * 4);
955 cudaMemcpy(devSegments_ptr, segments.data(), newPos.size() * sizeof(double) * 4, cudaMemcpyHostToDevice);
956
957 auto& cntrTreeSeg = *getCuda().cntrTreeSegment;
958 cntrTreeSeg.MemoryAllocate((int)getCuda().n_CUDA_wake);
959 cntrTreeSeg.UpdatePanelGeometry((int)newPos.size(), (double4*)devSegments_ptr);
960 cntrTreeSeg.Build();
961
962 BHcu::treePanelsSegmentsIntersectionCalculationWrapper(auxTree, cntrTreeSeg, getWake().devNearestPanelPtr);
963 cudaMemcpy(hit.data(), getWake().devNearestPanelPtr, newPos.size() * sizeof(int), cudaMemcpyDeviceToHost);
964 cudaFree(devSegments_ptr);
965 }
966
967 size_t panCounter = 0;
968 for (size_t afl = 0; afl < airfoil.size(); ++afl)
969 {
970 std::vector<double> gamma(airfoil[afl]->getNumberOfPanels(), 0.0);
971 for (int i = 0; i < newPos.size(); ++i)
972 {
973 if (hit[i] != -1)
974 {
975 if ((hit[i] >= panCounter) && (hit[i] < panCounter + airfoil[afl]->getNumberOfPanels()))
976 {
977 if (fabs(wake->vtx[i].g()) > 1.0)
978 std::cout << "Too large gamma is through: i = " << i << ", " << hit[i] << ", " << "gamma[hit[i] - panCounter] += " << wake->vtx[i].g() << std::endl;
979
980 gamma[hit[i] - panCounter] += wake->vtx[i].g();
981 wake->vtx[i].g() = 0.0;
982 }
983 }
984 }
985 airfoil[afl]->gammaThrough = gamma;
986 panCounter += airfoil[afl]->getNumberOfPanels();
987 }//for afl
988 }
990
991#endif
992
993 getTimers().stop("Inside");
994}
const Wake & getWake() const
Возврат константной ссылки на вихревой след
Definition World2D.h:232
std::vector< std::unique_ptr< AirfoilGeometry > > oldAirfoil
Список умных указателей на обтекаемые профили для сохранения старого положения
Definition World2D.h:83
VMlib::vmTimer timerInside
Definition World2D.h:372
Gpu & getNonConstCuda() const
Возврат неконстантной ссылки на объект, связанный с видеокартой (GPU)
Definition World2D.h:278
const Gpu & getCuda() const
Возврат константной ссылки на объект, связанный с видеокартой (GPU)
Definition World2D.h:273
size_t nVtxBeforeMerging
Definition WorldGen.h:76
const vmTimer & stop() const
Останов работающего счетчика времени
Definition TimesGen.h:125
const vmTimer & start() const
Запуск (первый или повторный) счетчика времени
Definition TimesGen.h:109
const vmTimer & reset() const
Сброс счетчика времени
Definition TimesGen.h:101
Here is the call graph for this function:
Here is the caller graph for this function:

◆ FillIQ()

void World2D::FillIQ ( )

Заполнение матрицы, состоящей из интегралов от (r-xi) / |r-xi|^2.

Вызывается в Step()

Definition at line 2066 of file World2D.cpp.

2067{
2068 getTimers().start("MatRhs");
2069
2070 for (size_t bou = 0; bou < boundary.size(); ++bou)
2071 {
2072 if (currentStep == 0 || mechanics[bou]->isDeform)
2073 {
2074 boundary[bou]->FillIQSelf(IQ[bou][bou]);
2075 }
2076
2077 for (size_t oth = 0; oth < boundary.size(); ++oth)
2078 {
2079
2080#ifdef INITIAL
2081 if (currentStep == 0 || !useInverseMatrix)
2082 {
2083 //size_t nVarsOther = boundary[oth]->GetUnknownsSize();
2084 if (bou != oth)
2085 {
2086 //std::cout << "matr!" << std::endl;
2087 boundary[bou]->FillIQFromOther(*boundary[oth], IQ[bou][oth]);
2088 }
2089 }// if (currentStep == 0 || !useInverseMatrix)
2090#endif
2091
2092#ifdef BRIDGE
2093 if (currentStep == 0)
2094 {
2095 //size_t nVarsOther = boundary[oth]->GetUnknownsSize();
2096 if (bou != oth)
2097 boundary[bou]->FillIQFromOther(*boundary[oth], IQ[bou][oth]);
2098 }// if (currentStep == 0)
2099#endif
2100
2101 }// for oth
2102
2103 }// for bou
2104
2105 getTimers().stop("MatRhs");
2106}//FillIQ()
Here is the call graph for this function:
Here is the caller graph for this function:

◆ FillMatrixAndRhs()

void World2D::FillMatrixAndRhs ( )

Заполнение матрицы системы для всех профилей

Вызывается в Step()

Definition at line 1298 of file World2D.cpp.

1299{
1300 getTimers().start("MatRhs");
1301
1302 const int linSystemScheme = passport.numericalSchemes.linearSystemSolver.second;
1303
1304 if (linSystemScheme == 0 || linSystemScheme == 1)
1305 {
1308
1309 Eigen::MatrixXd locMatr;
1310 Eigen::MatrixXd otherMatr;
1311 Eigen::VectorXd locLastLine, locLastCol;
1312
1313 std::vector<std::vector<Point2D>> locIQ;
1314 std::vector<std::vector<Point2D>> locOtherIQ;
1315
1316 //обнуляем матрицу на первом шаге расчета
1317 if (currentStep == 0)
1318 {
1319 for (int i = 0; i < matrReord.rows(); ++i)
1320 for (int j = 0; j < matrReord.cols(); ++j)
1321 matrReord(i, j) = 0.0;
1322 }
1323
1324 size_t currentRow = 0;
1325 size_t currentSkosRow = 0;
1326
1327 size_t nAllVars = 0;
1328 for (size_t bou = 0; bou < boundary.size(); ++bou)
1329 nAllVars += boundary[bou]->GetUnknownsSize();
1330
1331 for (size_t bou = 0; bou < boundary.size(); ++bou)
1332 {
1333 size_t nVars = boundary[bou]->GetUnknownsSize();
1334 if (currentStep == 0 || mechanics[bou]->isDeform)
1335 {
1336 locMatr.resize(nVars, nVars);
1337 locLastLine.resize(nVars);
1338 locLastCol.resize(nVars);
1339 }
1340
1341 if (currentStep == 0 || mechanics[bou]->isDeform)
1342 boundary[bou]->FillMatrixSelf(locMatr, locLastLine, locLastCol);
1343
1344
1345 //размазываем матрицу
1346 for (size_t i = 0; i < nVars; ++i)
1347 {
1348 if (currentStep == 0 || mechanics[bou]->isDeform)
1349 {
1350 for (size_t j = 0; j < nVars; ++j)
1351 matrReord(i + currentSkosRow, j + currentSkosRow) = locMatr(i, j);
1352
1353 matrReord(nAllVars + bou, i + currentSkosRow) = locLastLine(i);
1354 matrReord(i + currentSkosRow, nAllVars + bou) = locLastCol(i);
1355 }
1356 }
1357
1358 if ((currentStep == 0) || (!useInverseMatrix))
1359 {
1360 size_t currentCol = 0;
1361 size_t currentSkosCol = 0;
1362 for (size_t oth = 0; oth < boundary.size(); ++oth)
1363 {
1364 size_t nVarsOther = boundary[oth]->GetUnknownsSize();
1365
1366 if (bou != oth)
1367 {
1368 otherMatr.resize(nVars, nVarsOther);
1369
1370 boundary[bou]->FillMatrixFromOther(*boundary[oth], otherMatr);
1371
1372 //размазываем матрицу
1373 for (size_t i = 0; i < nVars; ++i)
1374 {
1375 for (size_t j = 0; j < nVarsOther; ++j)
1376 matrReord(i + currentSkosRow, j + currentSkosCol) = otherMatr(i, j);
1377 }
1378 }// if (bou != oth)
1379 currentCol += nVarsOther + 1;
1380 currentSkosCol += nVarsOther;
1381 }// for oth
1382 }// if (currentStep == 0 || mechanics[oth]->isMoves)
1383
1384 currentRow += nVars + 1;
1385 currentSkosRow += nVars;
1386 }// for bou
1388 }
1389
1390 velocity->FillRhs(rhsReord);
1391
1392 getTimers().stop("MatRhs");
1393}//FillMatrixAndRhs()
Eigen::MatrixXd matrReord
Матрица системы
Definition World2D.h:122
Eigen::VectorXd rhsReord
Правая часть системы
Definition World2D.h:134
VMlib::vmTimer timerFillMatrix
Definition World2D.h:368
Here is the call graph for this function:
Here is the caller graph for this function:

◆ GenerateMechanicsHeader()

void World2D::GenerateMechanicsHeader ( size_t  mechanicsNumber)

Метод-обертка для вызова метода генерации заголовка файла нагрузок и заголовка файла положения (последнее — если профиль движется)

Parameters
[in]mechanicsNumberномер профиля, для которого генерируется заголовок файла

Definition at line 2059 of file World2D.cpp.

2060{
2061 mechanics[mechanicsNumber]->GenerateForcesHeader();
2062 mechanics[mechanicsNumber]->GeneratePositionHeader();
2063}//GenerateMechanicsHeader(...)
Here is the caller graph for this function:

◆ getAirfoil()

const Airfoil & VM2D::World2D::getAirfoil ( size_t  i) const
inline

Возврат константной ссылки на объект профиля

Parameters
[in]iномер профиля, константная ссылка на который возвращается
Returns
константную ссылку на i-й профиль

Definition at line 163 of file World2D.h.

163{ return *airfoil[i]; };
Here is the caller graph for this function:

◆ getAuxTreePnl()

CpuTreeInfo & VM2D::World2D::getAuxTreePnl ( ) const
inline

Definition at line 241 of file World2D.h.

241{ return *auxTreePnl; };
std::unique_ptr< CpuTreeInfo > auxTreePnl
Definition World2D.h:119

◆ getBoundary()

const Boundary & VM2D::World2D::getBoundary ( size_t  i) const
inline

Возврат константной ссылки на объект граничного условия

Parameters
[in]iномер граничного условия, константная ссылка на которое возвращается
Returns
константную ссылку на i-е граничное условие

Definition at line 186 of file World2D.h.

186{ return *boundary[i]; };
Here is the caller graph for this function:

◆ getCntrTreePnl()

CpuTreeInfo & VM2D::World2D::getCntrTreePnl ( ) const
inline

Definition at line 243 of file World2D.h.

243{ return *cntrTreePnl; }
Here is the caller graph for this function:

◆ getCntrTreeVP()

CpuTreeInfo & VM2D::World2D::getCntrTreeVP ( ) const
inline

Definition at line 242 of file World2D.h.

242{ return *cntrTreeVP; };
Here is the caller graph for this function:

◆ getCntrTreeWake()

CpuTreeInfo & VM2D::World2D::getCntrTreeWake ( ) const
inline

Definition at line 240 of file World2D.h.

240{ return *cntrTreeWake; };
Here is the caller graph for this function:

◆ getCuda()

const Gpu & VM2D::World2D::getCuda ( ) const
inline

Возврат константной ссылки на объект, связанный с видеокартой (GPU)

Returns
константную ссылку на объект, связанный с видеокартой (GPU)

Definition at line 273 of file World2D.h.

273{ return cuda; };
Here is the caller graph for this function:

◆ getCurrentStep()

size_t VMlib::WorldGen::getCurrentStep ( ) const
inlineinherited

Возврат константной ссылки на параметры распараллеливания по MPI.

Returns
константную ссылку на параметры распараллеливания по MPI const Parallel& getParallel() const { return parallel; };

Возврат номера текущего временного шага

Returns
номера текущего временного шага

Definition at line 99 of file WorldGen.h.

99{ return currentStep; };
Here is the caller graph for this function:

◆ getCurrentTime()

double VMlib::WorldGen::getCurrentTime ( ) const
inlineinherited

Definition at line 100 of file WorldGen.h.

100{ return currentTime; };
Here is the caller graph for this function:

◆ getDispBoundaryInSystem()

size_t VM2D::World2D::getDispBoundaryInSystem ( size_t  i) const
inline

Возврат смещения в системе dispBoundaryInSystem.

Parameters
[in]iномер граничного условия, константная ссылка на которое возвращается
Returns
константную ссылку на i-е граничное условие

Definition at line 203 of file World2D.h.

203{ return dispBoundaryInSystem[i]; };
std::vector< size_t > dispBoundaryInSystem
Список номеров, с которых начинаются элементы правой части (или матрицы) системы для профилей
Definition World2D.h:89

◆ getInflTreeWake()

CpuTreeInfo & VM2D::World2D::getInflTreeWake ( ) const
inline

Definition at line 239 of file World2D.h.

239{ return *inflTreeWake; };
Here is the caller graph for this function:

◆ getInfo() [1/2]

VMlib::LogStream & VMlib::WorldGen::getInfo ( ) const
inlineinherited

Возврат ссылки на объект LogStream Используется в техничеcких целях для организации вывода

Returns
ссылку на объект LogStream

Definition at line 82 of file WorldGen.h.

82{ return info; };
Here is the caller graph for this function:

◆ getInfo() [2/2]

std::ostream & VMlib::WorldGen::getInfo ( char  x) const
inlineinherited

Возврат ссылки на поток вывода информации Необходимо для вывода телеметрической информации, информации об ошибках и т.п.

Parameters
[in]xсимвол, определяющий стиль вывода сообщения
Returns
ссылку на поток вывода информации

Definition at line 89 of file WorldGen.h.

89{ return info(x); };

◆ getIQ()

const std::pair< Eigen::MatrixXd, Eigen::MatrixXd > & VM2D::World2D::getIQ ( size_t  i,
size_t  j 
) const
inline

Возврат константной ссылки на объект, связанный с матрицей интегралов от (r-xi)/|r-xi|^2.

Returns
константную ссылку на объект, связанный с матрицей интегралов от (r-xi)/|r-xi|^2

Definition at line 283 of file World2D.h.

283{ return IQ[i][j]; };
Here is the caller graph for this function:

◆ getMeasureVP()

const MeasureVP & VM2D::World2D::getMeasureVP ( ) const
inline

Возврат константной ссылки на measureVP.

Returns
константную ссылку на вихревой след

Definition at line 208 of file World2D.h.

208{ return *measureVP; };
Here is the caller graph for this function:

◆ getMechanics()

const Mechanics & VM2D::World2D::getMechanics ( size_t  i) const
inline

Возврат константной ссылки на объект механики

Parameters
[in]iномер механики, константная ссылка на который возвращается
Returns
константную ссылку на i-ю механику

Definition at line 219 of file World2D.h.

219{ return *mechanics[i]; };
Here is the caller graph for this function:

◆ getMechanicsVector()

const std::vector< std::unique_ptr< Mechanics > > & VM2D::World2D::getMechanicsVector ( ) const
inline

Definition at line 227 of file World2D.h.

227{ return mechanics; };
Here is the caller graph for this function:

◆ getNonConstAirfoil()

Airfoil & VM2D::World2D::getNonConstAirfoil ( size_t  i) const
inline

Возврат неконстантной ссылки на объект профиля

Parameters
[in]iномер профиля, неконстантная ссылка на который возвращается
Returns
неконстантную ссылку на i-й профиль

Definition at line 175 of file World2D.h.

175{ return *airfoil[i]; };
Here is the caller graph for this function:

◆ getNonConstBoundary()

Boundary & VM2D::World2D::getNonConstBoundary ( size_t  i) const
inline

Возврат неконстантной ссылки на объект граничного условия

Parameters
[in]iномер граничного условия, неконстантная ссылка на которое возвращается
Returns
неконстантную ссылку на i-е граничное условие

Definition at line 192 of file World2D.h.

192{ return *boundary[i]; };
Here is the caller graph for this function:

◆ getNonConstCuda()

Gpu & VM2D::World2D::getNonConstCuda ( ) const
inline

Возврат неконстантной ссылки на объект, связанный с видеокартой (GPU)

Returns
неконстантную ссылку на объект, связанный с видеокартой (GPU)

Definition at line 278 of file World2D.h.

278{ return cuda; };
Here is the caller graph for this function:

◆ getNonConstMeasureVP()

MeasureVP & VM2D::World2D::getNonConstMeasureVP ( ) const
inline

Возврат неконстантной ссылки на measureVP.

Returns
неконстантную ссылку на вихревой след

Definition at line 213 of file World2D.h.

213{ return *measureVP; };
Here is the caller graph for this function:

◆ getNonConstMechanics()

Mechanics & VM2D::World2D::getNonConstMechanics ( size_t  i) const
inline

Возврат неконстантной ссылки на объект механики

Parameters
[in]iномер механики, константная ссылка на который возвращается
Returns
неконстантную ссылку на i-ю механику

Definition at line 225 of file World2D.h.

225{ return *mechanics[i]; };
Here is the caller graph for this function:

◆ getNonConstPassport()

Passport & VM2D::World2D::getNonConstPassport ( ) const
inline

Возврат неконстантной ссылки на паспорт

Returns
неконстантную ссылку на паспорт

Definition at line 268 of file World2D.h.

268{ return const_cast<Passport&>(passport); };
Here is the caller graph for this function:

◆ getNonConstVelocity()

Velocity & VM2D::World2D::getNonConstVelocity ( ) const
inline

Возврат неконстантной ссылки на объект для вычисления скоростей

Returns
неконстантную ссылку на объект для вычисления скоростей

Definition at line 258 of file World2D.h.

258{ return *velocity; };

◆ getNonConstWake()

Wake & VM2D::World2D::getNonConstWake ( ) const
inline

Возврат неконстантной ссылки на вихревой след

Returns
неконстантную ссылку на вихревой след

Definition at line 237 of file World2D.h.

237{ return *wake; };
Here is the caller graph for this function:

◆ getNumberOfAirfoil()

size_t VM2D::World2D::getNumberOfAirfoil ( ) const
inline

Возврат количества профилей в задаче

Returns
количество профилей в задаче

Definition at line 180 of file World2D.h.

180{ return airfoil.size(); };
Here is the caller graph for this function:

◆ getNumberOfBoundary()

size_t VM2D::World2D::getNumberOfBoundary ( ) const
inline

Возврат количества граничных условий в задаче

Returns
количество граничных условий в задаче

Definition at line 197 of file World2D.h.

197{ return boundary.size(); };
Here is the caller graph for this function:

◆ getOldAirfoil()

const AirfoilGeometry & VM2D::World2D::getOldAirfoil ( size_t  i) const
inline

Возврат константной ссылки на объект старого профиля

Parameters
[in]iномер старого профиля, константная ссылка на который возвращается
Returns
константную ссылку на i-й старый профиль

Definition at line 169 of file World2D.h.

169{ return *oldAirfoil[i]; };
Here is the caller graph for this function:

◆ getPassport()

const Passport & VM2D::World2D::getPassport ( ) const
inline

Возврат константной ссылки на паспорт

Returns
константную ссылку на паспорт

Definition at line 263 of file World2D.h.

263{ return passport; };
Here is the caller graph for this function:

◆ getPassportGen()

const PassportGen & VMlib::WorldGen::getPassportGen ( ) const
inlineinherited

Definition at line 102 of file WorldGen.h.

102{ return passportGen; };
const PassportGen & passportGen
Константная ссылка на паспорт конкретного расчета
Definition WorldGen.h:63
Here is the caller graph for this function:

◆ getSource()

const WakeDataBase & VM2D::World2D::getSource ( ) const
inline

Возврат константной ссылки на источники в области течения

Returns
константную ссылку на источники в области течения

Definition at line 248 of file World2D.h.

248{ return *source; };
Here is the caller graph for this function:

◆ getTimers()

VMlib::TimersGen & VM2D::World2D::getTimers ( ) const
inline

Возврат ссылки на временную статистику выполнения шага расчета по времени

Returns
ссылку на временную статистику выполнения шага расчета по времени

Definition at line 288 of file World2D.h.

288{ return *timers; };
Here is the caller graph for this function:

◆ getV0()

Point2D VM2D::World2D::getV0 ( ) const
inline

Возврат текущей скорости набегающего потока

Definition at line 148 of file World2D.h.

149 {
151 }
Point2D V0(double currentTime) const
Функция скорости набегающего потока с учетом разгона
Definition Passport2D.h:93
Here is the call graph for this function:
Here is the caller graph for this function:

◆ getVelocity()

const Velocity & VM2D::World2D::getVelocity ( ) const
inline

Возврат константной ссылки на объект для вычисления скоростей

Returns
константную ссылку на объект для вычисления скоростей

Definition at line 253 of file World2D.h.

253{ return *velocity; };
Here is the caller graph for this function:

◆ getWake()

const Wake & VM2D::World2D::getWake ( ) const
inline

Возврат константной ссылки на вихревой след

Returns
константную ссылку на вихревой след

Definition at line 232 of file World2D.h.

232{ return *wake; };
Here is the caller graph for this function:

◆ ifDivisible()

bool VM2D::World2D::ifDivisible ( int  val) const
inline

Definition at line 290 of file World2D.h.

290{ return ((val > 0) && (!(currentStep % val))); };
Here is the caller graph for this function:

◆ isAnyMovable()

bool World2D::isAnyMovable ( ) const

Возврат признака того, что хотя бы один из профилей подвижный

Definition at line 2363 of file World2D.cpp.

2364{
2365 return std::any_of(getMechanicsVector().begin(), getMechanicsVector().end(),
2366 [](const std::unique_ptr<Mechanics>& m) {return (m->isMoves); });
2367}
const std::vector< std::unique_ptr< Mechanics > > & getMechanicsVector() const
Definition World2D.h:227
Here is the call graph for this function:

◆ isAnyMovableOrDeformable()

bool World2D::isAnyMovableOrDeformable ( ) const

Возврат признака того, что хотя бы один из профилей подвижный или деформируемый

Definition at line 2370 of file World2D.cpp.

2371{
2372 return std::any_of(getMechanicsVector().begin(), getMechanicsVector().end(),
2373 [](const std::unique_ptr<Mechanics>& m) {return (m->isMoves || m->isDeform); });
2374}
Here is the call graph for this function:
Here is the caller graph for this function:

◆ isFinished()

bool WorldGen::isFinished ( ) const
inherited

Функция, возвращающая признак завершения счета в решаемой задаче

true, если задача решена и выполнен признак останова; false если требуется еще выполнять шаги по времени

Definition at line 52 of file WorldGen.cpp.

53{
55}
double timeStop
Конечное время
Definition PassportGen.h:64

◆ MoveVortexes()

void World2D::MoveVortexes ( std::vector< Point2D > &  newPos,
std::vector< double > *  nutPtr = nullptr 
)

Вычисляем новые положения вихрей (в пелене и виртуальных)

Вызывается в Step()

Parameters
[out]newPosновые позиции вихрей

Definition at line 1655 of file World2D.cpp.

1656{
1657 //getTimers().start("Move");
1658
1659 size_t nvt = wake->vtx.size();
1660 size_t nVirtVortex = 0;
1661 for (size_t i = 0; i < getNumberOfBoundary(); ++i)
1662 nVirtVortex += boundary[i]->virtualWake.vtx.size();
1663
1664
1665
1666
1667 size_t counter = wake->vtx.size() - nVirtVortex;
1668 for (size_t bou = 0; bou < boundary.size(); ++bou)
1669 {
1670 for (int i = 0; i < (int)boundary[bou]->virtualWake.vtx.size(); ++i)
1671 {
1672 wake->vtx[counter].g() = boundary[bou]->virtualWake.vtx[i].g();
1673 wake->vtx[counter].sigma() = getPassport().wakeDiscretizationProperties.sigma0;
1674 ++counter;
1675 }
1676 }
1677
1678
1679
1680
1682#if defined (SPH)
1683
1684 const double c = 4.5;
1685
1686 const auto& eps = velocity->wakeVortexesParams.epsastWake;
1687
1688 double minEpsast = *std::min_element(eps.begin(), eps.end());
1689 double meanEpsast = std::accumulate(eps.begin(), eps.end(), 0.0) / eps.size();
1690
1691 double uniformH = minEpsast + 0.2 * (meanEpsast - minEpsast);
1692
1693 std::vector<std::vector<int>> neib(wake->vtx.size());
1694 for (auto& nb : neib)
1695 nb.reserve(100);
1696
1697#pragma omp parallel for schedule(dynamic, 100)
1698 for (int i = 0; i < (int)wake->vtx.size(); ++i)
1699 {
1700 const double hati = c * eps[i];
1701 const double hati2 = sqr(hati);
1702 const int dist = std::ceil(hati / uniformH);
1703
1704 Point2D vtx = wake->vtx[i].r();
1705 int cx = (int)(vtx[0] / uniformH);
1706 int cy = (int)(vtx[1] / uniformH);
1707
1708 for (int j = 0; j < (int)wake->vtx.size(); ++j)
1709 {
1710 Point2D test = wake->vtx[j].r();
1711 int tx = (int)(test[0] / uniformH);
1712 int ty = (int)(test[1] / uniformH);
1713
1714 if (abs(cx - tx) <= dist && abs(cy - ty) <= dist)
1715 if ((vtx - test).length2() < hati2)
1716 neib[i].push_back(j);
1717 }
1718 }
1719
1720 //Ядро иранца
1721 //auto W = [](double xi, double R) {double t = xi / R; return (t > 1 ? 0.0 : (8.0 / ( PI * R * R)) * (0.75 * t * t - 1.5 * t + 0.75)); };
1722 //auto gW = [](double xi, double R) {double t = xi / R; return (t > 1 ? 0.0 : (12.0 / (PI * sqr(R * R))) * (1.0 - 1.0/t)); };
1723
1724 //M4
1725 Point2D zeroVec{ 0.0, 0.0 };
1726 auto W = [](double xi, double R) {double t = xi / R; return (t > 1 ? 0.0 : 40.0 / (7.0 * PI * R * R) * (t > 0.5 ? 2.0 * cubPower(1.0 - t) : 1.0 - 6.0 * sqr(t) * (1.0 - t))); };
1727 auto WW = [zeroVec](const Point2D& v, double R) {double t = v.length() / R; return (t > 1 ? zeroVec : v * (240.0 / (7.0 * PI * sqr(sqr(R))) * (t > 0.5 ? -sqr(1.0 - t) / t : -2.0 + 3.0 * t))); };
1728
1729 //Poly6
1730 //auto W = [](double xi, double R) {double t = xi / R; return (t > 1 ? 0.0 : 4.0 / (PI * R * R) * cubPower(1.0 - t * t)); };
1731 //auto gW = [](double xi, double R) {double t = xi / R; return (t > 1 ? 0.0 : -24.0 / (PI * sqr(R * R)) * sqr(1.0 - t * t)); };
1732
1733 //Poly6 + grad Spicky
1734 //auto W = [](double xi, double R) {double t = xi / R; return (t > 1 ? 0.0 : 4.0 / (PI * R * R) * cubPower(1.0 - t * t)); };
1735 //auto gW = [](double xi, double R) {double t = xi / R; return (t > 1 ? 0.0 : -30.0 / (PI * sqr(R * R)) * sqr(1.0 - t) / t); };
1736
1737 //WendlandC6
1738 /*
1739 Point2D zeroVec{ 0.0, 0.0 };
1740 auto W = [](double xi, double R) {double t = xi / R;
1741 return (t > 1 ? 0.0 : 78.0 / (7.0 * PI * sqr(R)) * sqr(sqr(sqr((1.0 - t)))) * (32.0 * cubPower(t)+ 25.0 * sqr(t) + 8.0*t + 1.0)); };
1742 auto WW = [zeroVec](const Point2D& v, double R) -> Point2D {double t = v.length() / R;
1743 return (t > 1 ? zeroVec : v * (-1716.0 / (7.0 * PI * sqr(sqr(R))) * sqr(cubPower(1.0 - t)) * (1.0 - t) * (16.0 * sqr(t) + 7.0 * t + 1.0))); };
1744 */
1745
1746 //fp5
1747 /*
1748 auto W = [](double xi, double R) {
1749 double t = xi / R;
1750 if (t >= 1) return 0.0;
1751 if (t >= 2.0 / 3.0)
1752 return 1701.0 / (478.0 * PI) * sqr(sqr(1.0 - t)) * (1.0 - t);
1753 if (t >= 1.0 / 3.0)
1754 return 1701.0 / (478.0 * PI) * (sqr(sqr(1.0 - t)) * (1.0 - t) - 6.0 * sqr(sqr((2.0 / 3.0 - t))) * (2.0 / 3.0 - t));
1755
1756 return 1701.0 / (478.0 * PI) * (sqr(sqr(1.0 - t)) * (1.0 - t) - 6.0 * sqr(sqr((2.0 / 3.0 - t))) * (2.0 / 3.0 - t) + 15.0 * sqr(sqr((1.0 / 3.0 - t))) * (1.0 / 3.0 - t));
1757 };
1758
1759 auto gW = [](double xi, double R) {
1760 double t = xi / R;
1761 if (t >= 1) return 0.0;
1762 if (t >= 2.0 / 3.0)
1763 return -8505.0 / (478.0 * PI * R * R) * sqr(sqr(1.0 - t)) / t;
1764 if (t >= 1.0 / 3.0)
1765 return -8505.0 / (478.0 * PI * R * R) * (sqr(sqr(1.0 - t)) - 6.0 * sqr(sqr((2.0 / 3.0 - t)))) / t;
1766
1767 return -8505.0 / (478.0 * PI * R * R) * (sqr(sqr(1.0 - t)) - 6.0 * sqr(sqr((2.0 / 3.0 - t))) + 15.0 * sqr(sqr((1.0 / 3.0 - t)))) / t;
1768 };
1769 */
1770
1771
1772 std::vector<double> OIP(wake->vtx.size(), 0.0);
1773 std::vector<double> VIP(wake->vtx.size(), 0.0);
1774 std::vector<Point2D> GIP(wake->vtx.size(), { 0.0, 0.0 });
1775 std::vector<Point2D> DIFFV(wake->vtx.size(), { 0.0, 0.0 });
1776
1777
1778 const auto& vtx = wake->vtx;
1779 const auto& vel = velocity->wakeVortexesParams.convVelo;
1780
1781
1782//#ifdef TURB
1783// std::vector<Point2D> grNu(wake->vtx.size(), { 0.0, 0.0 });
1784//
1785// //std::vector<Point2D> grVx(wake->vtx.size(), { 0.0, 0.0 });
1786// //std::vector<Point2D> grVy(wake->vtx.size(), { 0.0, 0.0 });
1787// std::vector<nummatrix<double,2,2>> grV(wake->vtx.size());
1788// std::vector<Point2D> TP(wake->vtx.size(), { 0.0, 0.0 });
1789// std::vector<double> RIP(wake->vtx.size(), 0.0);
1790// std::vector<double> DELTAG(wake->vtx.size(), 0.0);
1791//#endif //TURB
1792
1793 std::vector<nummatrix<double, 2, 2>> mtrL(wake->vtx.size(), { {0.0, 0.0}, {0.0, 0.0} });
1794 std::vector<std::vector<double>> weight(wake->vtx.size());
1795
1796
1797 std::vector<Point2D> sumDIFFV(wake->vtx.size(), {0.0, 0.0});
1798 std::vector<double> sumWeight(wake->vtx.size(), 0.0);
1799
1800#pragma omp parallel for schedule(dynamic, 100)
1801 for (int i = 0; i < (int)vtx.size(); ++i)
1802 {
1803 const double hati = c * eps[i];
1804
1805 double resultOIP = 0.0;
1806 //double resultIVIP = 0.0;
1807
1808 const Point2D myPos = vtx[i].r();
1809
1810 double I1 = 0.0, I0 = 1.0;
1811
1812
1813 weight[i].reserve(neib[i].size());
1814
1815 for (auto& nbr : neib[i])
1816 {
1817 double dist2 = (myPos - vtx[nbr].r()).length2();
1818 double dist = sqrt(dist2);
1819 double w = W(dist, hati);
1820 I1 += vtx[nbr].g() * w;
1821 weight[i].push_back(w);
1822 }
1823
1824 OIP[i] = I1 / I0;
1825 }
1826
1827#pragma omp parallel for schedule(dynamic, 100)
1828 for (int i = 0; i < (int)vtx.size(); ++i)
1829 {
1830 const Point2D myPos = vtx[i].r();
1831 const double hati = c * eps[i];
1832
1833 nummatrix<double, 2, 2> M = { {0.0, 0.0}, {0.0, 0.0} };
1834 GIP[i].toZero();
1835
1836 for (auto& nbr : neib[i])
1837 {
1838 Point2D dr = myPos - vtx[nbr].r();
1839 Point2D grweight = WW(dr, hati);
1840 double vol = (vtx[nbr].g() / OIP[nbr]);
1841 if (vol != vol)
1842 vol = 0.0;
1843
1844 M += (dr | grweight) * vol;
1845 GIP[i] -= grweight * ((OIP[nbr] - OIP[i]) * vol);
1846 }
1847
1848 double detM = M[0][0] * M[1][1] - M[1][0] * M[0][1];
1849 nummatrix<double, 2, 2> L{ {0.0, 0.0}, {0.0, 0.0} };
1850 if (fabs(detM) > 1e-2)
1851 L = nummatrix<double, 2, 2>{ { M[1][1], -M[0][1] }, { -M[1][0], M[0][0] } } *(1.0 / detM);
1852
1853 GIP[i] = (L & GIP[i]);
1854 DIFFV[i] = GIP[i] * (1.0 / OIP[i]);
1855 //GIP[i] = -GIP[i];
1856 }
1857
1858 //осредняем
1859 //*
1860//#pragma omp parallel for schedule(dynamic, 100)
1861 for (int i = 0; i < (int)vtx.size(); ++i)
1862 for (int j = 0; j < neib[i].size(); ++j)
1863 if ((vtx[i].r() - vtx[neib[i][j]].r()).length() <= 1.0 * eps[i])
1864 {
1865 sumDIFFV[i] += weight[i][j] * DIFFV[neib[i][j]];
1866 sumWeight[i] += weight[i][j];
1867 }
1868
1869 //#pragma omp parallel for schedule(dynamic, 100)
1870 for (int i = 0; i < (int)vtx.size(); ++i)
1871 {
1872 DIFFV[i] = sumDIFFV[i] * (1.0 / sumWeight[i]);
1873 }
1874 //*/
1875
1876//#ifdef TURB
1877//#pragma omp parallel for schedule(dynamic, 100)
1878// for (int i = 0; i < (int)vtx.size(); ++i)
1879// {
1880// Point2D myPos = vtx[i].r();
1881// Point2D resultGradNu = { 0.0, 0.0 };
1882// Point2D resultGradVx = { 0.0, 0.0 };
1883// Point2D resultGradVy = { 0.0, 0.0 };
1884//
1885// //if (i == 0)
1886// //{
1887// // std::cout << "h = " << h << std::endl;
1888// // std::cout << "VIP = " << VIP[0] << std::endl;
1889// // std::cout << "my_coord = " << vtx[i].r() << std::endl;
1890// // std::cout << "mtrC = " << mtrC[i] << std::endl;
1891// //}
1892//
1893// for (auto& nbr : neib[i])
1894// {
1895// //if (i == 0)
1896// // std::cout << "nb_" << nbr << "_coord = " << vtx[nbr].r() << std::endl;
1897//
1898// Point2D dr = myPos - vtx[nbr].r();
1899// double dist = dr.length();
1900// double gweight = gW(dist, h);
1901//
1902// if (nbr != i)
1903// {
1904// double nutNbr = (*nutPtr)[nbr];
1905// double nutI = (*nutPtr)[i];
1906// //double nutNbr = vtx[nbr].r() & Point2D { 1.0, 1.0 };
1907// //double nutSlf = vtx[i].r() & Point2D { 1.0, 1.0 };
1908//
1909// resultGradNu += (nutNbr - nutI) * (mtrL[i] & dr) * (gweight * /*vtx[nbr].g() / OIP[nbr]*/ VIP[nbr]);
1910// resultGradVx += (vel[nbr][0] - vel[i][0]) * (mtrL[i] & dr) * (gweight * /*vtx[nbr].g() / OIP[nbr]*/ VIP[nbr]);
1911// resultGradVy += (vel[nbr][1] - vel[i][1]) * (mtrL[i] & dr) * (gweight * /*vtx[nbr].g() / OIP[nbr]*/ VIP[nbr]);
1912// }
1913// }
1914//
1915// grNu[i] = resultGradNu;
1916// grV[i][0] = resultGradVx;
1917// grV[i][1] = resultGradVy;
1918//
1919// nummatrix<double, 2, 2> regGrV = grV[i];
1920// regGrV[0][0] -= grV[i][1][1];
1921// regGrV[0][1] += grV[i][1][0];
1922// regGrV[1][0] += grV[i][0][1];
1923// regGrV[1][1] -= grV[i][0][0];
1924//
1925// regGrV *= 0.5;
1926//
1927// regGrV[0][1] += 0.5 * OIP[i];
1928// regGrV[1][0] -= 0.5 * OIP[i];
1929//
1930// TP[i] = (regGrV & grNu[i]);
1931// }//for i
1932//#endif //TURB
1933
1934
1935//#ifdef TURB
1936//
1937//#pragma omp parallel for schedule(dynamic, 100)
1938// for (int i = 0; i < (int)vtx.size(); ++i)
1939// {
1940// Point2D myPos = vtx[i].r();
1941//
1942// Point2D resultGradTPx = { 0.0, 0.0 };
1943// Point2D resultGradTPy = { 0.0, 0.0 };
1944//
1945// for (auto& nbr : neib[i])
1946// {
1947// Point2D dr = myPos - vtx[nbr].r();
1948// double dist = dr.length();
1949// double gweight = gW(dist, h);
1950//
1951// if (nbr != i)
1952// {
1953// resultGradTPx += (TP[nbr][0] - TP[i][0]) * (mtrL[i] & dr) * (gweight * /*vtx[nbr].g() / OIP[nbr]*/ VIP[nbr]);
1954// resultGradTPy += (TP[nbr][1] - TP[i][1]) * (mtrL[i] & dr) * (gweight * /*vtx[nbr].g() / OIP[nbr]*/ VIP[nbr]);
1955// }
1956// }
1957//
1958// RIP[i] = resultGradTPy[0] - resultGradTPx[1];
1960// double dG = RIP[i] * VIP[i] * getPassport().timeDiscretizationProperties.dt;
1961// //if (i == 0)
1962// // std::cout << "nVirtVortex = " << nVirtVortex << std::endl;
1963// DELTAG[i] = 0.0;//dG;// (i < (wake->vtx.size() - nVirtVortex) ? dG : 0.0);
1964//
1965// //if (fabs(dG / getWake().vtx[i].g()) <= 0.25)
1966// getNonConstWake().vtx[i].g() += DELTAG[i];
1967// //else
1968// //{
1969// // if (dG * getWake().vtx[i].g() > 0)
1970// // getNonConstWake().vtx[i].g() *= xxx;
1971// // else
1972// // getNonConstWake().vtx[i].g() /= xxx;
1973// //}
1974//
1975// getNonConstWake().vtx[i].g() = std::clamp(getWake().vtx[i].g(), -1.5 * getPassport().wakeDiscretizationProperties.maxGamma, 1.5 * getPassport().wakeDiscretizationProperties.maxGamma);
1976// }//for i
1977//
1978// wake->SaveScalarFields("SPH", VIP, OIP, *nutPtr, grNu, RIP, DELTAG);
1979//#endif //TURB
1980
1981#endif //SPH
1982
1983
1984
1985 //newPos.clear();
1986 newPos.resize(nvt);
1987
1988
1989 // 1. Инициализация генератора случайных чисел (ПСЧГ)
1990 //std::random_device rd; // Источник энтропии для инициализации
1991 //std::mt19937 gen(rd()); // Mersenne Twister, инициализированный случайным значением
1992
1993 // 2. Определение параметров нормального распределения
1994 // μ (среднее) = 0.0, σ (стандартное отклонение) = 1.0
1995 //std::normal_distribution<> d(0.0, sqrt(2.0 * getPassport().physicalProperties.nu * getPassport().timeDiscretizationProperties.dt));
1996
1997
1998 //std::ofstream checkVel (passport.dir + "/dbg/checkVel" + std::to_string(getCurrentStep()) + ".txt");
1999 //std::ofstream checkVelDyn(passport.dir + "/dbg/checkVelDyn" + std::to_string(getCurrentStep()) + ".txt");
2000
2001 //checkVel << "i x y epsast Vx Vy Wx Wy corrWx corrWy Om gOmx gOmy\n";
2002 //checkVelDyn << "i x y epsast Vx Vy Wx Wy corrWx corrWy Om gOmx gOmy\n";
2003
2004//#pragma omp parallel for
2005 for (int i = 0; i < (int)wake->vtx.size(); ++i)
2006 {
2007#ifdef SPH
2008 //Point2D W = -(passport.physicalProperties.nu / OIP[i]) * GIP[i];
2009 Point2D W = -passport.physicalProperties.nu * DIFFV[i];
2010
2011 checkVel << i << " " << wake->vtx[i].r()[0] << " " << wake->vtx[i].r()[1] << " " << velocity->wakeVortexesParams.epsastWake[i] << " " \
2012 << velocity->wakeVortexesParams.convVelo[i][0] << " " << velocity->wakeVortexesParams.convVelo[i][1] << " " \
2013 << W[0] << " " << W[1];
2014
2015 if (W.length() > 1.5 * passport.physicalProperties.vRef)
2016 {
2017 std::cout << "i = " << i << ", W was " << W << std::endl;
2018 W.toZero();
2019 //W.normalize(1.5 * passport.physicalProperties.vRef);
2020 }
2021
2022 checkVel << " " << W[0] << " " << W[1] << " " << OIP[i] << " " << GIP[i][0] << " " << GIP[i][1] << "\n";
2023
2024#else
2025 Point2D W = velocity->wakeVortexesParams.diffVelo[i] * (nutPtr ? (1.0 + (*nutPtr)[i] / passport.physicalProperties.nu) : 1.0);
2026
2027 const auto& diff = velocity->wakeVortexesParams;
2028 Point2D grad = diff.I2[i] * (1.0 / (diff.epsastWake[i] * diff.I0[i])) - diff.I3[i] * (diff.I1[i] / sqr(diff.I0[i]));
2029
2030 //checkVelDyn << i << " " << wake->vtx[i].r()[0] << " " << wake->vtx[i].r()[1] << " " << velocity->wakeVortexesParams.epsastWake[i] << " " \
2031 // << velocity->wakeVortexesParams.convVelo[i][0] << " " << velocity->wakeVortexesParams.convVelo[i][1] << " " \
2032 // << W[0] << " " << W[1] << " " << W[0] << " " << W[1] << " " << diff.I1[i]/diff.I0[i] << " " << grad[0] << " " << grad[1] << "\n";
2033#endif
2034
2035 newPos[i] = wake->vtx[i].r() + (velocity->wakeVortexesParams.convVelo[i] + W + getV0()) * passport.timeDiscretizationProperties.dt;
2036
2037 //Chorin Random Walk method
2038 //newPos[i] = wake->vtx[i].r() + (velocity->wakeVortexesParams.convVelo[i] + getV0()) * passport.timeDiscretizationProperties.dt + Point2D{ d(gen), d(gen) };
2039 }
2040 //checkVel.close();
2041 //checkVelDyn.close();
2042
2043
2044 //for (int i = 0; i < (int)wake->vtx.size(); ++i)
2045 //{
2046 // std::cout << i << " " << wake->vtx[i].r() << " " << \
2047 // velocity->wakeVortexesParams.convVelo[i] << " " << \
2048 // velocity->wakeVortexesParams.diffVelo[i] << " " << \
2049 // getV0() << "\n";
2050 //}
2051
2052
2053
2054
2055 //getTimers().stop("Move");
2056}//MoveVortexes(...)
const double PI
Число .
Definition defs.h:76
Point2D getV0() const
Возврат текущей скорости набегающего потока
Definition World2D.h:148
Шаблонный класс, определяющий матрицу фиксированного размера Фактически представляет собой массив,...
Definition nummatrix.h:66
nummatrix< T, n, m > & toZero(T val=0.0)
Definition nummatrix.h:300
numvector< T, n > & toZero(P val=0)
Установка всех компонент вектора в константу (по умолчанию — нуль)
Definition numvector.h:528
size_t size() const
Definition numvector.h:114
P length() const
Вычисление 2-нормы (длины) вектора
Definition numvector.h:374
R dist(const numvector< T, n > &x, const numvector< P, n > &y)
Вычисление расстояния между двумя точками
Definition numvector.h:852
T cubPower(T x)
Возведение числа в куб
Definition defs.h:468
T sqr(T x)
Возведение числа в квадрат
Definition defs.h:460
auto dist2(const numvector< T, n > &x, const numvector< P, n > &y) -> typename std::remove_const< decltype(x[0] - y[0])>::type
Вычисление квадрата расстояния между двумя точками
Definition numvector.h:836
double abs(const Vector3< T > &vec)
Definition utils.h:122
double vRef
Референсная скорость
Definition Passport2D.h:81
double nu
Коэффициент кинематической вязкости среды
Definition Passport2D.h:99
Here is the call graph for this function:
Here is the caller graph for this function:

◆ ReserveMemoryForMatrixAndRhs()

void World2D::ReserveMemoryForMatrixAndRhs ( )

Вычисляем размер матрицы и резервируем память под нее и под правую часть

Вызывается в Step()

Definition at line 1396 of file World2D.cpp.

1397{
1398 //getTimers().start("Mem");
1399
1400 if (currentStep == 0)
1401 {
1402 dispBoundaryInSystem.resize(boundary.size());
1403 dispBoundaryInSystem[0] = 0;
1404
1405 for (size_t i = 1; i < boundary.size(); ++i)
1406 {
1407 dispBoundaryInSystem[i] = dispBoundaryInSystem[i - 1] + boundary[i - 1]->GetUnknownsSize() + 1;
1408 }
1409
1410 size_t matrSize = boundary.size();
1411 size_t matrSkosSize = 0;
1412
1413 for (auto it = boundary.begin(); it != boundary.end(); ++it)
1414 {
1415 matrSize += (*it)->GetUnknownsSize();
1416 matrSkosSize += (*it)->GetUnknownsSize();
1417 }
1418
1419 if ((getPassport().numericalSchemes.linearSystemSolver.second == 0 || getPassport().numericalSchemes.linearSystemSolver.second == 1 || (getPassport().numericalSchemes.linearSystemSolver.second == 2 && isAnyMovableOrDeformable())))
1420 {
1421 //matr.resize(matrSize, matrSize);
1422 matrReord.resize(matrSize, matrSize);
1423 //matrSkos.resize(matrSkosSize, matrSkosSize);
1424
1425 //matr.setZero();
1426 matrReord.setZero();
1427 //matrSkos.setZero();
1428
1429
1430 for (size_t i = 0; i < getNumberOfAirfoil(); ++i)
1431 {
1432 size_t nVari = boundary[i]->GetUnknownsSize();
1433 for (size_t j = 0; j < getNumberOfAirfoil(); ++j)
1434 {
1435 size_t nVarj = boundary[j]->GetUnknownsSize();
1436 IQ[i][j].first.resize(nVari, nVarj);
1437 IQ[i][j].second.resize(nVari, nVarj);
1438 }
1439 }
1440 }
1441
1442 //rhs.resize(matrSize);
1443 rhsReord.resize(matrSize);
1444 //rhsSkos.resize(matrSkosSize);
1445 }
1446 //rhs.setZero();
1447 rhsReord.setZero();
1448 //rhsSkos.setZero();
1449
1450 //getTimers().stop("Mem");
1451}//ReserveMemoryForMatrixAndRhs()
Here is the call graph for this function:
Here is the caller graph for this function:

◆ SolveLinearSystem()

void World2D::SolveLinearSystem ( )

Решение системы линейных алгебраических уравнений

Вызывается в Step()

Definition at line 1001 of file World2D.cpp.

1002{
1003 getTimers().start("Solve");
1004/*
1005 if (currentStep == 0)
1006 {
1007 invMatr = matr;
1008 std::ifstream fileMatrix("matrix.txt");
1009 int nx, ny;
1010 fileMatrix >> nx;
1011 fileMatrix >> ny;
1012 for (int i = 0; i < matr.rows(); ++i)
1013 {
1014 for (int j = 0; j < matr.cols(); ++j)
1015 {
1016 fileMatrix >> matr(i, j);
1017 }
1018 }
1019 fileMatrix.close();
1020 }
1021*/
1022
1023//*
1024 //if (currentStep == 0)
1025 //{
1026 // std::ofstream fileMatrix(passport.dir + "/dbg/matrix.txt");
1027 // //int nx, ny;
1028 // //fileMatrix >> nx;
1029 // //fileMatrix >> ny;
1030 // fileMatrix.precision(16);
1031 // for (int i = 0; i < matrReord.rows(); ++i)
1032 // {
1033 // for (int j = 0; j < matrReord.cols(); ++j)
1034 // {
1035 // fileMatrix << matrReord(i, j) << " ";
1036 // }
1037 // fileMatrix << std::endl;
1038 // }
1039 // fileMatrix.close();
1040 //}
1041
1042
1043
1044 //{
1045 // std::ofstream fileMatrix(passport.dir + "/dbg/rhs"+std::to_string(currentStep)+".txt");
1046 // //int nx, ny;
1047 // //fileMatrix >> nx;
1048 // //fileMatrix >> ny;
1049 // fileMatrix.precision(16);
1050 // for (int i = 0; i < rhsReord.size(); ++i)
1051 // {
1052 // fileMatrix << rhsReord(i) << std::endl;
1053 // }
1054 // fileMatrix.close();
1055 //}//*/
1056
1057 //exit(-100500);
1058
1060 const int linSystemScheme = passport.numericalSchemes.linearSystemSolver.second;
1061
1062 if (linSystemScheme == 0) // Gaussian elimination
1063 {
1064 double t1 = -omp_get_wtime();
1065 if (useInverseMatrix && (currentStep == 0))
1066 {
1067 info('t') << "Inverting matrix... ";
1068
1069#if (defined(USE_CUDA))
1070 invMatr.resize(matrReord.rows(), matrReord.cols());
1071 for (int i = 0; i < (int)matrReord.rows(); ++i)
1072 for (int j = 0; j < (int)matrReord.cols(); ++j)
1073 invMatr(i, j) = (i == j) ? 1.0 : 0.0;
1074 cuInverseMatrix((int)matrReord.rows(), matrReord.data(), invMatr.data());
1075#else
1076 invMatr = matrReord.inverse();
1077#endif
1078 info('t') << "done" << std::endl;
1079 }
1080
1081 if (currentStep == 0)
1082 info('t') << "Solving system at step #0... ";
1083
1084 if (useInverseMatrix)
1085 {
1086 /*
1087 std::cout << "Solution via invMatrix" << std::endl;
1088 for (int i = 0; i < invMatr.rows(); ++i)
1089 for (int j = 0; j < invMatr.cols(); ++j)
1090 if (fabs(invMatr(i, j)) > 1e+5)
1091 std::cout << "invMatrixError: " << "invA(" << i << ", " << j << ") = " << invMatr(i, j) << std::endl;
1092
1093 for (int i = 0; i < rhsReord.size(); ++i)
1094 if (fabs(rhsReord(i)) > 1e+5)
1095 std::cout << "rhsReordError: " << "rhsReord(" << i << ") = " << rhsReord(i) << std::endl;
1096 */
1099 sol = invMatr * rhsReord;
1101 }
1102 else
1103 {
1106 sol = matrReord.partialPivLu().solve(rhsReord);
1108 }
1109
1110 if (currentStep == 0)
1111 info('t') << "done" << std::endl;
1112
1113 t1 += omp_get_wtime();
1114
1115 //std::cout << " Time in Gauss = " << t1 << std::endl;
1116/*
1117 std::ofstream solFile(getPassport().dir + "/sol" + std::to_string(currentStep) + "-Gauss.txt");
1118 solFile.precision(16);
1119 for (int i = 0; i < sol.size(); ++i)
1120 solFile << sol(i) << std::endl;
1121 solFile.close();
1122 exit(10101);
1123//*/
1124 }//Gauss
1125
1126/*
1127 if (currentStep == 0)
1128 { N
1129 invMatr = matr;
1130 std::ifstream fileMatrix("invMatrix.txt");
1131 int nx, ny;
1132 fileMatrix >> nx;
1133 fileMatrix >> ny;
1134 for (int i = 0; i < invMatr.rows(); ++i)
1135 {
1136 for (int j = 0; j < invMatr.cols(); ++j)
1137 {
1138 fileMatrix >> invMatr(i, j);
1139 }
1140 }
1141 fileMatrix.close();
1142 }
1143*/
1144
1145/*
1146 if (currentStep == 0)
1147 {
1148 Eigen::MatrixXd mmul = matr * invMatr;
1149 for (int i = 0; i < invMatr.rows(); ++i)
1150 mmul(i,i) -= 1.0;
1151
1152 double maxElem = 0.0;
1153 for (int i = 0; i < mmul.rows(); ++i)
1154 for (int j = 0; j < mmul.cols(); ++j)
1155 if (fabs(mmul(i,j)) > maxElem )
1156 maxElem = mmul(i,j);
1157
1158 std::cout << "A*A^-1: MAX ELEMENT = " << maxElem << std::endl;
1159 }
1160*/
1161
1162/*
1163 if (currentStep == 0)
1164 {
1165 std::ofstream fileMatrix("invMatrix4x3200.txt");
1166 fileMatrix << invMatr.rows() << " " << invMatr.cols() << std::endl;
1167 fileMatrix.precision(18);
1168 for (int i = 0; i < invMatr.rows(); ++i)
1169 {
1170 for (int j = 0; j < invMatr.cols(); ++j)
1171 {
1172 fileMatrix << invMatr(i, j) << " ";
1173 }
1174 fileMatrix << std::endl;
1175 }
1176 fileMatrix.close();
1177 }
1178*/
1179
1180
1181
1182 if ((linSystemScheme == 1 || linSystemScheme == 2)) // GMRES
1183 {
1184 int nFullVars = (int)getNumberOfBoundary();
1185 for (int i = 0; i < getNumberOfBoundary(); ++i)
1186 nFullVars += (int)(boundary[i]->GetUnknownsSize());
1187
1188 std::vector<std::vector<double>> Ggam(getNumberOfBoundary());
1189 for (int i = 0; i < getNumberOfBoundary(); ++i)
1190 Ggam[i].resize(boundary[i]->GetUnknownsSize());
1191
1192 std::vector<double> GR(getNumberOfBoundary());
1193
1194#ifndef USE_CUDA
1195 std::vector<double> Grhs(rhsReord.size());
1196 for (int i = 0; i < rhsReord.size(); ++i)
1197 Grhs[i] = rhsReord(i);
1198
1199 std::vector<int> Gpos(getNumberOfBoundary(), 0);
1200 for (int i = 1; i < getNumberOfBoundary(); ++i)
1201 Gpos[i] = Gpos[i - 1] + (int)(boundary[i - 1]->GetUnknownsSize());
1202
1203 std::vector<int> Gvsize(getNumberOfBoundary());
1204 for (int i = 0; i < getNumberOfBoundary(); ++i)
1205 Gvsize[i] = (int)(boundary[i]->GetUnknownsSize());
1206
1208 {
1209 //for direct GMRES
1210 std::vector<double> Gmatr(nFullVars * nFullVars);
1211 for (size_t i = 0; i < nFullVars; ++i)
1212 for (size_t j = 0; j < nFullVars; ++j)
1213 Gmatr[i * nFullVars + j] = matrReord(i, j);
1214
1215 for (int i = 0; i < getNumberOfBoundary(); ++i)
1216 for (int j = 0; j < boundary[i]->GetUnknownsSize(); ++j)
1217 auto y_test = airfoil[i]->len[j % airfoil[i]->getNumberOfPanels()];
1218
1219 //Gmres->GMRES_Direct(*this, nFullVars, (int)getNumberOfBoundary(), Gmatr, Grhs, Gpos, Gvsize, Ggam, GR);
1220 std::cout << "GMRES_Direct should be modified!" << std::endl;
1221 exit(-100);
1222
1223
1224 sol.resize(nFullVars);
1225 int cntr = 0;
1226 for (int i = 0; i < getNumberOfBoundary(); ++i)
1227 for (int j = 0; j < boundary[i]->GetUnknownsSize(); ++j)
1228 sol(cntr++) = Ggam[i][j] /*/ airfoil[i]->len[j % airfoil[i]->getNumberOfPanels()]*/;
1229 for (int i = 0; i < getNumberOfBoundary(); ++i)
1230 sol(cntr++) = GR[i];
1231
1232 //std::ofstream solFile(getPassport().dir + "/dbg/sol-direct-gmres" + std::to_string(currentStep) + ".txt");
1233 //solFile.precision(16);
1234 //for (int i = 0; i < sol.size(); ++i)
1235 // solFile << sol(i) << std::endl;
1236 //solFile.close();
1237 }
1239 {
1240 getInfo('e') << "Fast GMRES without CUDA is not implemented" << std::endl;
1241 exit(-2);
1242 }
1243#else
1245 {
1246 getInfo('e') << "Direct GMRES for CUDA is not implemented" << std::endl;
1247 exit(-2);
1248 }
1249
1250 if ( (passport.numericalSchemes.linearSystemSolver.second == 2))
1251 {
1252 // for fast GMRES
1253 std::vector<std::vector<double>> GGrhs(getNumberOfBoundary());
1254 int cntrRhs = 0;
1255 for (int i = 0; i < getNumberOfBoundary(); ++i)
1256 GGrhs[i].resize(boundary[i]->GetUnknownsSize());
1257 for (int i = 0; i < getNumberOfBoundary(); ++i)
1258 for (int j = 0; j < boundary[i]->GetUnknownsSize(); ++j)
1259 GGrhs[i][j] = rhsReord(cntrRhs++);
1260
1261 std::vector<double> GrhsReg(getNumberOfBoundary());
1262 for (int i = 0; i < getNumberOfBoundary(); ++i)
1263 GrhsReg[i] = rhsReord[nFullVars - (getNumberOfBoundary() - i)];
1264
1265
1266 int niter;
1268
1269 double time_GMRES = -omp_get_wtime();
1270 getNonConstCuda().Gmres->GMRES(Ggam, GR, GGrhs, GrhsReg, niter);//, linScheme);
1271 time_GMRES += omp_get_wtime();
1272 //std::cout << "Time_GMRES = " << time_GMRES << std::endl;
1273
1274 sol.resize(nFullVars);
1275 int cntr = 0;
1276 for (int i = 0; i < getNumberOfBoundary(); ++i)
1277 for (int j = 0; j < boundary[i]->GetUnknownsSize(); ++j)
1278 sol(cntr++) = Ggam[i][j] / airfoil[i]->len[j % airfoil[i]->getNumberOfPanels()];
1279 for (int i = 0; i < getNumberOfBoundary(); ++i)
1280 sol(cntr++) = GR[i];
1281
1282 /*
1283 std::ofstream solFile(getPassport().dir + "/sol-fast-new-gmres" + std::to_string(currentStep) + ".txt");
1284 solFile.precision(16);
1285 for (int i = 0; i < sol.size(); ++i)
1286 solFile << sol(i) << std::endl;
1287 solFile.close();
1288 exit(1010110);
1289 //*/
1290 }
1291#endif
1292 }
1293
1294 getTimers().stop("Solve");
1295}//SolveLinearSystem()
Eigen::MatrixXd invMatr
Обратная матрица
Definition World2D.h:128
VMlib::vmTimer timerSlaeSolve
Definition World2D.h:369
Here is the call graph for this function:
Here is the caller graph for this function:

◆ Step()

void World2D::Step ( )
overridevirtual

Функция выполнения предварительного шага

Основная функция выполнения одного шага по времени

Implements VMlib::WorldGen.

Definition at line 279 of file World2D.cpp.

280{
281 //try {
282 //Очистка статистики
284
285#ifdef USE_CUDA
286 cuSetCurrentStep((int)currentStep, 1);
287#endif // USE_CUDA
288
289
290 //Засечка времени в начале шага
291 getTimers().start("Step");
292
293 //Обновление параметров паспорта
295 {
296 std::unique_ptr<Passport> revisePspPtr;
297 revisePspPtr.reset(new VM2D::Passport(info, passport.dir, passport.problemNumber,
300
302
303 }
304
305
307
308 size_t countStrongCoupling = 0;
309 for (size_t m = 0; m < mechanics.size(); ++m)
310 {
311 MechanicsRigidOscillPart* mechVar = dynamic_cast<MechanicsRigidOscillPart*>(mechanics[m].get());
312 if (mechVar && getPassport().airfoilParams[m].addedMass.length2() > 0)
313 {
314 mechVar->getStrongCoupling() = true;
315 ++countStrongCoupling;
316 }
317 }
318
319 bool semiImplicitStrategy = ((countStrongCoupling == mechanics.size()) && (mechanics.size() > 0));
320 if ((currentStep == 0) && (semiImplicitStrategy))
321 info('i') << "Strong (semi-implicit) coupling strategy" << std::endl;
322
323//НЕПОДВИЖНЫЕ ТЕЛА
324//*
325 int nTotPan = 0;
326 for (size_t s = 0; s < getNumberOfAirfoil(); ++s)
327 nTotPan += (int)getAirfoil(s).getNumberOfPanels();
328
329
330 if (!semiImplicitStrategy)
331 {
333
336 if (getPassport().numericalSchemes.velocityComputation.second == 0 && getPassport().numericalSchemes.linearSystemSolver.second == 2)
337 {
338#ifdef USE_CUDA
339 cuda.RefreshAfls(3);
340 if (nTotPan > 0)
341 {
342 auto& afl = getAirfoil(0);
343 auto& treePnlVrt = *getCuda().inflTreePnlVortex;
344 if (getCurrentStep() == 0)
345 treePnlVrt.MemoryAllocate((int)getCuda().n_CUDA_pnls);
346
348 {
349 //Построение дерева для влияющих панелей (вихревых слоев)
350 treePnlVrt.UpdatePanelGeometry(nTotPan, (double4*)afl.devRPtr);
351 treePnlVrt.Build();
352 }
353
354 treePnlVrt.UpdatePanelAttachedVortexIntensity(afl.devAttachedVortexSheetPtr, afl.devAttachedVortexSheetLinPtr);
355 treePnlVrt.UpwardTraversal(getPassport().numericalSchemes.gmresMultipoleOrder);
356 }
357#endif
358 }
359
360 if (getPassport().numericalSchemes.velocityComputation.second == 1)
361 {
362#ifdef USE_CUDA
363 cuda.RefreshWake(3);
364 cuda.RefreshAfls(3);
365
366 //Построение дерева для вихрей
367 auto& treeWake = *getCuda().inflTreeWake;
368 treeWake.MemoryAllocate((int)getCuda().n_CUDA_wake);
369 treeWake.Update((int)getWake().vtx.size(), getWake().devVtxPtr);
370 treeWake.Build();
371 treeWake.UpwardTraversal(getPassport().numericalSchemes.nbodyMultipoleOrder);
372
373 if (nTotPan > 0)
374 {
375 auto& afl = getAirfoil(0);
376 auto& treePnl = *getCuda().cntrTreePnl;
377 auto& treePnlVrt = *getCuda().inflTreePnlVortex;
378 auto& treePnlSrc = *getCuda().inflTreePnlSource;
379 auto& treePnlAux = *getCuda().auxTreePnl;
380
381 if (getCurrentStep() == 0)
382 {
383 treePnl.MemoryAllocate((int)getCuda().n_CUDA_pnls);
384 treePnlAux.MemoryAllocate((int)getCuda().n_CUDA_pnls);
385 treePnlVrt.MemoryAllocate((int)getCuda().n_CUDA_pnls);
387 treePnlSrc.MemoryAllocate((int)getCuda().n_CUDA_pnls);
388 }
389
391 {
392 //Построение дерева для влияющих панелей (вихревых слоев)
393 treePnlVrt.UpdatePanelGeometry(nTotPan, (double4*)afl.devRPtr);
394 treePnlVrt.Build();
395
396 //Построение контрольного дерева панелей
397 treePnl.UpdatePanelGeometry((int)nTotPan, (double4*)afl.devRPtr);
398 treePnl.Build();
399
400 //Построение дерева для влияющих панелей (источников)
402 {
403 treePnlSrc.UpdatePanelGeometry(nTotPan, (double4*)afl.devRPtr);
404 treePnlSrc.UpdatePanelAttachedSourceIntensity(afl.devAttachedSourceSheetPtr, afl.devAttachedSourceSheetLinPtr);
405 treePnlSrc.Build();
406 treePnlSrc.UpwardTraversal(getPassport().numericalSchemes.nbodyMultipoleOrder);
407 }
408 }
409
410 treePnlVrt.UpdatePanelAttachedVortexIntensity(afl.devAttachedVortexSheetPtr, afl.devAttachedVortexSheetLinPtr);
411 treePnlVrt.UpwardTraversal(getPassport().numericalSchemes.nbodyMultipoleOrder);
412
413
414 //временное дерево для отладки процедуры сохранения
415 /*
416 if (getWake().vtx.size() > 2)
417 {
418 CpuTreeInfo testTree(tree_T::vortex, object_T::point4, scheme_T::noScheme);
419 testTree.Update(getWake().vtx);
420 testTree.Build();
421 testTree.UpwardTraversal(getPassport().numericalSchemes.nbodyMultipoleOrder);
422 testTree.SaveTreeTopology(getPassport().dir + "testTree-" + std::to_string(currentStep) + ".tree");
423 }
424 */
425 }
426#else
427 if (getWake().vtx.size() > 0)
428 {
429 inflTreeWake->Update(getWake().vtx);
430 inflTreeWake->Build();
431 inflTreeWake->UpwardTraversal(getPassport().numericalSchemes.nbodyMultipoleOrder);
432 }
433 std::vector<std::pair<Point2D, Point2D>> panels;
434
435 for (size_t bou = 0; bou < getNumberOfAirfoil(); ++bou)
436 {
437 const auto& afl = getAirfoil(bou);
438
439 for (size_t i = 0; i < afl.getNumberOfPanels(); ++i)
440 panels.push_back({ afl.getR(i), afl.getR(i + 1) });
441 }
442 if (panels.size() > 0)
443 {
444 cntrTreePnl->UpdatePanelGeometry(panels, std::max(4, (int)(log2(panels.size())) - 2));
445 cntrTreePnl->Build();
446 cntrTreePnl->UpwardTraversal(getPassport().numericalSchemes.nbodyMultipoleOrder);
447 }
448#endif
449 }
451
452 measureVP->Initialization();
454
455 //added masses
456 if (getPassport().physicalProperties.typeAccel.second == 3)
457 {
458 std::vector<Point2D> lambdaAdd(getNumberOfAirfoil(), {0.0, 0.0});
459 std::vector<double> muAdd(getNumberOfAirfoil());
460 for (size_t bou = 0; bou < getNumberOfAirfoil(); ++bou)
461 {
462 if (dynamic_cast<MechanicsRigidGivenLaw*>(&getNonConstMechanics(bou)))
463 {
464 for (size_t pnl = 0; pnl < getAirfoil(bou).getNumberOfPanels(); ++pnl)
465 {
466 double sumGam = getBoundary(bou).sheets.freeVortexSheet(pnl, 0) + getBoundary(bou).sheets.attachedVortexSheet(pnl, 0);
467 Point2D rpnl = 0.5 * (getAirfoil(bou).getR(pnl) + getAirfoil(bou).getR(pnl + 1));
468 lambdaAdd[bou] += rpnl.kcross() * sumGam * getAirfoil(bou).len[pnl];
469 muAdd[bou] += (rpnl - getAirfoil(bou).rcm).length2() * sumGam * getAirfoil(bou).len[pnl];
470 }
471 lambdaAdd[bou] *= getPassport().physicalProperties.rho;
472 muAdd[bou] *= 0.5 * getPassport().physicalProperties.rho;
473 info('i') << "Added Masses for airfoil #" << bou << " = { " << lambdaAdd[bou][0] << ", " << lambdaAdd[bou][1] << ", " << muAdd[bou] << " }" << std::endl;
474
475 char direction;
476 switch ((int)(getPassport().physicalProperties.timeAccel))
477 {
478 case 0:
479 direction = 'x';
480 break;
481 case 1:
482 direction = 'y';
483 break;
484 case 2:
485 direction = 'w';
486 break;
487 default:
488 direction = '?';
489 info('e') << "Wrong dirfection is specified!" << std::endl;
490 exit(1);
491 }
492
493 std::string addMassFileName = getPassport().dir + "addMass-" + std::to_string(getAirfoil(bou).numberInPassport);
494 std::ofstream addMassFile;
495 if (!VMlib::fileExistTest(addMassFileName, info) || getPassport().physicalProperties.timeAccel == 0)
496 {
497 addMassFile.open(addMassFileName);
498 addMassFile << "Added Masses for airfoil:" << std::endl;
499 }
500 else
501 addMassFile.open(addMassFileName, std::ios_base::app);
502
503 addMassFile << direction << "-direction: " << lambdaAdd[bou][0] << " " << lambdaAdd[bou][1] << " " << muAdd[bou] << std::endl;
504
505 addMassFile.close();
506 }
507 }
509 }
510
512 if (nTotPan > 0 && getPassport().numericalSchemes.velocityComputation.second == 1)
513 {
514 auto& afl = getAirfoil(0);
515#ifdef USE_CUDA
516 getCuda().inflTreePnlVortex->UpdatePanelFreeAndAttachedVortexIntensity(afl.devFreeVortexSheetPtr, afl.devFreeVortexSheetLinPtr, afl.devAttachedVortexSheetPtr, afl.devAttachedVortexSheetLinPtr);
517 getCuda().inflTreePnlVortex->UpwardTraversal(getPassport().numericalSchemes.nbodyMultipoleOrder);
518#endif
519 }
521
522 //Вычисление скоростей вихрей: для тех, которые в следе, и виртуальных, а также в точках wakeVP
524
525
526#ifdef TURB
527 getTimers().start("Nut");
528 std::vector<double> nut;
529 CalcVeloDifference(nut);
530 getTimers().stop("Nut");
531
532 for (size_t q = 0; q < getWake().vtx.size(); ++q)
533 getNonConstWake().vtx[q].sigma() = 4.48364 * sqrt(passport.timeDiscretizationProperties.dt * (passport.physicalProperties.nu + nut[q]));
534#endif
535
536
537//#include "gammaCirc.h"
538
539 //Расчет и сохранение поля давления
541
542#ifdef OPTIMIZER
543 if ((getCurrentStep() >= OPTIMIZER_START_STEP && getCurrentStep() <= OPTIMIZER_STOP_STEP) || getCurrentStep()==0) //2026-03-28
544#endif
545 {
546 const double& vRef = getPassport().physicalProperties.vRef;
547 //scaleV = 1.0 / vRef;
548 double scaleP = 1.0 / (0.5 * getPassport().physicalProperties.rho * sqr(vRef));
549
550#ifdef USE_CUDA
551 measureVP->GPUCalcPressure();
552#else
553 measureVP->CalcPressure();
554#endif
555
556 if (mechanics.size()>0)
557 {
558 //std::cout << "presForces are calculated\n";
559 //MechanicsDeformable* ptr = dynamic_cast<MechanicsDeformable*>(mechanics[0].get());
560 Mechanics* ptr = mechanics[0].get();
561
562 if (measureVP->getTotalNumberOfRealPoints() > 0)
563 measureVP->SaveVP();
564
565 //if (ptr && !ptr->beam->fsi)
566
567 //сохранение главного вектора силы, вычисленного как интеграл от давления
568 //*
569 if (measureVP->elasticPoints.size() > 0)
570 {
571 std::ofstream presForcesFile;
572 if (currentStep == 0)
573 {
574 presForcesFile.open(getPassport().dir + "presForcesFile.csv");
575 presForcesFile << "time,Fx,Fy,Py" << std::endl;
576 }
577 else
578 presForcesFile.open(getPassport().dir + "presForcesFile.csv", std::ios_base::app);
579
580 auto r = measureVP->GetVPinElasticPoints();
581 Point2D presForce = { 0.0, 0.0 };
582 double yPower = 0.0;
583
584 for (int q = 0; q < r.size(); ++q)
585 {
586 presForce -= r[q].second * getAirfoil(0).len[q] * getAirfoil(0).nrm[q];
587 yPower -= (r[q].second * getAirfoil(0).len[q] * getAirfoil(0).nrm[q])[1] * (0.5 * (getAirfoil(0).getV(q) + getAirfoil(0).getV(q + 1))[1]);
588 }
589
590 presForcesFile << currentTime << "," << scaleP * presForce[0] << "," << scaleP * presForce[1] << "," << scaleP * yPower << std::endl;
591 presForcesFile.close();
592 }
593 }
594 //*/
595
596 if (mechanics.size() > 0)
597 {
598 //Коэффициенты разложения силы по балочным функциям
599 MechanicsDeformable* ptr = dynamic_cast<MechanicsDeformable*>(mechanics[0].get());
600 if (ptr && ptr->beam->fsi)
601 {
602 auto r = measureVP->GetVPinElasticPoints();
603
604 /*
605 std::ofstream elastPresFile;
606 elastPresFile.open(getPassport().dir + "elastPresFile-" + std::to_string(currentStep) + ".txt");
607 //if (currentStep == 0)
608 // elastPresFile.open(getPassport().dir + "elastPresFile.txt");
609 //else
610 // elastPresFile.open(getPassport().dir + "elastPresFile.txt", std::ios_base::app);
611 //elastPresFile << currentStep << " ";
612
613 for (size_t q = 0; q < r.size(); ++q)
614 elastPresFile << measureVP->elasticPoints[q][0] << " " << measureVP->elasticPoints[q][1] << " " << r[q].second << std::endl;
615 //elastPresFile << std::endl;
616 elastPresFile.close();
617 */
618
619 //Сохраняем давление на последних nLastSteps шагах по дисциплине очереди
620 std::vector<double> currentPres(ptr->chord.size());
621 for (size_t j = 0; j < ptr->chord.size(); ++j)
622 {
623 //currentPres[j] = -(ptr->beam->rho * 2 * ptr->beam->F); // gravity force for g = 2.0;
624 currentPres[j] = -(r[2 * j + 0].second - r[2 * j + 1].second);
625 }
626
627 if (ptr->beam->presLastSteps.size() < ptr->beam->nLastSteps)
628 ptr->beam->presLastSteps.push_back(currentPres);
629 else
630 {
631 for (int w = 1; w < ptr->beam->nLastSteps; ++w)
632 ptr->beam->presLastSteps[w - 1] = std::move(ptr->beam->presLastSteps[w]);
633 ptr->beam->presLastSteps.back() = currentPres;
634 }
635
636 for (int q = 0; q < ptr->beam->R; ++q)
637 {
638 ptr->beam->qCoeff[q] = 0;
639
640 if (ptr->beam->presLastSteps.size() == ptr->beam->nLastSteps)
641 {
642 for (size_t j = 0; j < ptr->chord.size(); ++j)
643 {
644 double averpres = 0.0;
645
646 //осредняем давление
647 for (int i = 0; i < ptr->beam->nLastSteps; ++i)
648 averpres += ptr->beam->presLastSteps[i][j];
649 averpres /= ptr->beam->nLastSteps;
650
651 ptr->beam->qCoeff[q] += -averpres * (ptr->initialChord[j].beg - ptr->initialChord[j].end).length() * ptr->beam->shape(q, 0.5 * (ptr->initialChord[j].beg + ptr->initialChord[j].end)[0]);
652 }
653 ptr->beam->qCoeff[q] /= (ptr->beam->intSqUnitShape * ptr->beam->L);
654 }
655
656 //std::cout << "q[" << q << "] = " << ptr->beam->qCoeff[q] << std::endl;
657 }
658
659 std::ofstream phiFile;
660 if (currentStep == 0)
661 {
662 phiFile.open(getPassport().dir + "phiFile.csv");
663 phiFile << "time";
664 for (int p = 0; p < ptr->beam->R; ++p)
665 phiFile << ",phi-" << std::to_string(p + 1);
666 for (int p = 0; p < ptr->beam->R; ++p)
667 phiFile << ",q-" << std::to_string(p + 1);
668 phiFile << std::endl;
669 }
670 else
671 phiFile.open(getPassport().dir + "phiFile.csv", std::ios_base::app);
672
673 phiFile << currentTime;
674 for (int p = 0; p < ptr->beam->R; ++p)
675 phiFile << "," << ptr->beam->phi(p, 0);
676 for (int p = 0; p < ptr->beam->R; ++p)
677 phiFile << "," << ptr->beam->qCoeff[p];
678
679 phiFile << std::endl;
680
681 phiFile.close();
682 }
683 }
684//*/
685 //getInfo('e') << "Deformable airfoils are not supported now!";
686 }
687
688 //Вычисление сил, действующих на профиль и сохранение в файл
689 for (auto& mech : mechanics)
690 {
691 mech->GetHydroDynamForce();
692 mech->GenerateForcesString();
693 mech->GeneratePositionString();
694 }
695
696 //Движение вихрей (сброс вихрей + передвижение пелены)
697#ifdef TURB
698 WakeAndAirfoilsMotion(true, &nut);
699#else
701#endif
702 wake->Restruct();
703 wake->SaveKadrVtk();
704
705 //std::string fname = getPassport().dir + "/airfoil-" + std::to_string(getAirfoil(0).numberInPassport) + "-" + std::to_string(currentStep) + ".pfl";
707 //{
708 // std::ofstream of(fname);
709 // for (size_t i = 0; i < getAirfoil(0).getNumberOfPanels(); ++i)
710 // of << getAirfoil(0).getR(i)[0] << " " << getAirfoil(0).getR(i)[1] << std::endl;
711 // of.close();
712 //}
713
714 }//if (!semiImplicitStrategy)
715
716
717//СОПРЯЖЕННАЯ ПОСТАНОВКА
718 if (semiImplicitStrategy)
719 {
721 measureVP->Initialization();
722
724
725 //Вычисление скоростей вихрей: для тех, которые в следе, и виртуальных, а также в точках wakeVP
727
728 //for (size_t s = 0; s < mechanics.size(); ++s) {// проверка: вычисление полной силы по циркуляциям новых вихрей
729 // mechanics[s]->GetHydroDynamForce();
730 // mechanics[s]->GenerateForcesString();
731 //}
732
733 //Расчет и сохранение поля давления
735 {
736 measureVP->CalcPressure();
737 measureVP->SaveVP();
738 }
739
740 WakeAndAirfoilsMotion(false); //профиль - кинематически, здесь Vold := V;
741 wake->Restruct();
742
745
746 //Вычисление сил, действующих на профиль и сохранение в файл
747
748 for (size_t m = 0; m < mechanics.size(); ++m)
749 {
750 mechanics[m]->GetHydroDynamForce();
751
752 MechanicsRigidOscillPart* mechVar = dynamic_cast<MechanicsRigidOscillPart*>(mechanics[m].get());
753 if (mechVar)
754 {
755 if (getPassport().airfoilParams[m].addedMass.length2() == 0)
756 {
757 getInfo('e') << "Added mass of the airfoil should be non-zero!" << std::endl;
758 exit(1);
759 }
760 mechVar->MoveOnlyVelo();
761 Point2D accel = (mechVar->getV() - mechVar->getVOld()) * (1.0 / getPassport().timeDiscretizationProperties.dt);
762 mechVar->hydroDynamForce -=
763 Point2D{ accel[0] * getPassport().airfoilParams[m].addedMass[0],
764 accel[1] * getPassport().airfoilParams[m].addedMass[1] };
765 mechVar->GenerateForcesString();
766 mechVar->GeneratePositionString();
767 }
768 else
769 exit(3333);
770 }
771 }
772
773
774 wake->SaveKadrVtk();
775//*/
777
778
779 //Засечка времени в конце шага
780 getTimers().stop("Step");
781 /*
782 std::cout << "nvt = " << nVtxBeforeMerging << std::endl;
783 std::cout << "tIBT = " << timerInitialBuild.duration() << std::endl;
784 std::cout << "tRHS = " << timerRhs.duration() << std::endl;
785 std::cout << "tFIL = " << timerFillMatrix.duration() << std::endl;
786 std::cout << "tLIN = " << timerSlaeSolve.duration() << std::endl;
787 std::cout << "tVEL = " << timerConvVelo.duration() << std::endl;
788 std::cout << "tINS = " << timerInside.duration() << std::endl;
789 std::cout << "tMRG = " << timerMerging.duration() << std::endl;
790 */
791
792 info('i') << "Step = " << getCurrentStep() \
793 << " PhysTime = " << getCurrentTime() \
794 << std::setprecision(3) \
795 << " StepTime = " << getTimers().durationStep() \
796 << std::setprecision(6) \
797 << std::endl;
798
800
802 currentStep += 1;
803 //}
804 //catch (...)
805 //{
806 // info('e') << "!!! Exception from unknown source !!!" << std::endl;
807 // exit(-1);
808 //}
809}//Step()
const Point2D & getR(size_t q) const
Возврат константной ссылки на вершину профиля
Definition Airfoil2D.h:113
std::vector< Point2D > nrm
Нормали к панелям профиля
Definition Airfoil2D.h:81
Point2D rcm
Положение центра масс профиля
Definition Airfoil2D.h:97
Sheet sheets
Слои на профиле
Definition Boundary2D.h:96
std::vector< ChordPanel > initialChord
std::vector< ChordPanel > chord
std::unique_ptr< Beam > beam
Абстрактный класс, определяющий вид механической системы
Definition Mechanics2D.h:72
Point2D hydroDynamForce
Вектор гидродинамической силы и момент, действующие на профиль
void GeneratePositionString()
Сохранение строки со статистикой в файл нагрузок
void GenerateForcesString()
Сохранение строки со статистикой в файл нагрузок
Point2D & getV()
текущая скорость профиля
std::string defaultsFileFullName
Definition Passport2D.h:284
std::vector< std::string > varLine
Definition Passport2D.h:286
std::string switchersFileFullName
Definition Passport2D.h:285
std::string mechanicsFileFullName
Definition Passport2D.h:283
std::string fileFullName
Имена файлов
Definition Passport2D.h:282
void GetReviseParamsFromParser(const Passport &newPassport, const std::vector< std::string > paramList)
Считывание измененных параметров
const double & attachedVortexSheet(size_t n, size_t moment) const
Definition Sheet2D.h:105
const double & freeVortexSheet(size_t n, size_t moment) const
Definition Sheet2D.h:100
std::vector< Vortex2D > vtx
Список вихревых элементов
void WakeAndAirfoilsMotion(bool dynamics, std::vector< double > *nutPtr=nullptr)
Перемещение вихрей и профилей на шаге
Definition World2D.cpp:2255
VMlib::vmTimer timerInitialBuild
Definition World2D.h:366
void CalcPanelsVeloAndAttachedSheets()
Вычисление скоростей панелей и интенсивностей присоединенных слоев вихрей и источников
Definition World2D.cpp:1637
void CalcVortexVelo()
Вычисление скоростей (и конвективных, и диффузионных) вихрей (в пелене и виртуальных),...
Definition World2D.cpp:1456
void CalcAndSolveLinearSystem()
Набор матрицы, правой части и решение СЛАУ
Definition World2D.cpp:2108
const Boundary & getBoundary(size_t i) const
Возврат константной ссылки на объект граничного условия
Definition World2D.h:186
bool ifDivisible(int val) const
Definition World2D.h:290
Mechanics & getNonConstMechanics(size_t i) const
Возврат неконстантной ссылки на объект механики
Definition World2D.h:225
Wake & getNonConstWake() const
Возврат неконстантной ссылки на вихревой след
Definition World2D.h:237
void GenerateStatString(size_t stepNo, double curTime, size_t N)
Формирование очередной строки файла временной статистики
Definition TimesGen.cpp:111
void resetAll()
Сброс всех счетчиков
Definition TimesGen.cpp:81
double durationStep() const
Вывод счетчика всего шага в секундах
Definition TimesGen.h:184
size_t getCurrentStep() const
Возврат константной ссылки на параметры распараллеливания по MPI.
Definition WorldGen.h:99
numvector< T, 2 > kcross() const
Геометрический поворот двумерного вектора на 90 градусов
Definition numvector.h:511
double rho
Плотность потока
Definition Passport2D.h:75
std::vector< std::string > reviseParameters
Список перечитываемых параметров
Definition PassportGen.h:73
int revisePassportStep
Шаг перечитывания паспорта
Definition PassportGen.h:70
Here is the call graph for this function:

◆ WakeAndAirfoilsMotion()

void World2D::WakeAndAirfoilsMotion ( bool  dynamics,
std::vector< double > *  nutPtr = nullptr 
)

Перемещение вихрей и профилей на шаге

Вызывается в Step()

Definition at line 2255 of file World2D.cpp.

2256{
2257 std::vector<Point2D> newPos;
2258
2259 MoveVortexes(newPos, nutPtr);
2260
2261//#ifdef BRIDGE
2262// double totalForce = 0;
2263// for (size_t afl = 0; afl < airfoil.size(); ++afl)
2264// {
2265// totalForce += mechanics[afl]->hydroDynamForce[1];
2266// }
2267// totalForce *= 1.2;
2268//
2269// mechanics[0]->hydroDynamForce[1] = totalForce;
2270//#endif
2271
2272 oldAirfoil.resize(0);
2273 for (auto& afl : airfoil)
2274 {
2275 if (dynamic_cast<AirfoilRigid*>(afl.get()))
2276 oldAirfoil.emplace_back(new AirfoilGeometry(*afl));
2277
2278 if (dynamic_cast<AirfoilDeformable*>(afl.get()))
2279 oldAirfoil.emplace_back(new AirfoilGeometry(*afl));
2280
2281
2282#ifdef BRIDGE
2283 if (dynamics)
2284 {
2285 if (afl->numberInPassport == 0)
2286 {
2287 mechanics[afl->numberInPassport]->Move();
2288 }
2289 else
2290 {
2291 MechanicsRigidOscillPart* mechTest = dynamic_cast<MechanicsRigidOscillPart*>(mechanics[afl->numberInPassport].get());
2292 if (mechTest != nullptr)
2293 {
2294 Mechanics& mechGen = *mechanics[0];
2295 MechanicsRigidOscillPart& mech0 = dynamic_cast<MechanicsRigidOscillPart&>(mechGen);
2296
2297 Point2D dr = mech0.getR() - mech0.getROld();
2298 Point2D dv = mech0.getV() - mech0.getVOld();
2299
2300 double dphi = mech0.getPhi() - mech0.getPhiOld();
2301 double dw = mech0.getW() - mech0.getWOld();
2302
2303 Mechanics& mechGenI = *mechanics[afl->numberInPassport];
2304 MechanicsRigidOscillPart& mechI = dynamic_cast<MechanicsRigidOscillPart&>(mechGenI);
2305
2306 mechI.getROld() = mechI.getR();
2307 mechI.getVOld() = mechI.getV();
2308 mechI.getPhiOld() = mechI.getPhi();
2309 mechI.getWOld() = mechI.getW();
2310
2311 //std::cout << "afl = " << afl << ", dy = " << dy << std::endl;
2312
2313 airfoil[afl->numberInPassport]->Move(dr);
2314 airfoil[afl->numberInPassport]->Rotate(dphi);
2315
2316 mechI.getR() += dr;
2317 mechI.getV() += dv;
2318
2319 mechI.getPhi() += dphi;
2320 mechI.getW() += dw;
2321 }
2322 }
2323 }
2324#endif
2325
2326#ifdef INITIAL
2327 if (dynamics)
2328 mechanics[afl->numberInPassport]->Move();
2329 else
2330 {
2331 MechanicsRigidOscillPart* mech = dynamic_cast<MechanicsRigidOscillPart*>(mechanics[afl->numberInPassport].get());
2332 if (mech)
2333 mech->MoveKinematic();
2334 else
2335 exit(2222);
2336 }
2337#endif
2338 }//for
2339
2340#if defined(__CUDACC__) || defined(USE_CUDA)
2341 cuda.RefreshAfls(2);
2342#endif
2343
2344 for (auto& bou : boundary)
2345 bou->virtualWake.vtx.clear();
2346
2347
2348 CheckInside(newPos, oldAirfoil);
2349
2350
2351 //передача новых положений вихрей в пелену
2352 for (size_t i = 0; i < wake->vtx.size(); ++i)
2353 {
2354 wake->vtx[i].r() = newPos[i];
2355 //wake->vtx[i].sigma() = getPassport().wakeDiscretizationProperties.sigma0;
2356 }
2357
2358// getWake().SaveKadrVtk();
2359
2360}//WakeAndAirfoilsMotion()
Класс, определяющий форму профиля
Definition Airfoil2D.h:66
Point2D & getR()
текущее отклонение профиля
double & getPhi()
текущий угол поворота профиля
double & getW()
текущая угловая скорость профиля
void MoveVortexes(std::vector< Point2D > &newPos, std::vector< double > *nutPtr=nullptr)
Вычисляем новые положения вихрей (в пелене и виртуальных)
Definition World2D.cpp:1655
void CheckInside(std::vector< Point2D > &newPos, const std::vector< std::unique_ptr< AirfoilGeometry > > &oldAirfoil)
Проверка проникновения вихрей внутрь профиля
Definition World2D.cpp:813
Here is the call graph for this function:
Here is the caller graph for this function:

Member Data Documentation

◆ airfoil

std::vector<std::unique_ptr<Airfoil> > VM2D::World2D::airfoil
private

Список умных указателей на обтекаемые профили

Definition at line 80 of file World2D.h.

◆ auxTreePnl

std::unique_ptr<CpuTreeInfo> VM2D::World2D::auxTreePnl

Definition at line 119 of file World2D.h.

◆ boundary

std::vector<std::unique_ptr<Boundary> > VM2D::World2D::boundary
private

Список умных указателей на формирователи граничных условий на профилях

Definition at line 86 of file World2D.h.

◆ check01

int VM2D::World2D::check01
mutable

Definition at line 108 of file World2D.h.

◆ check02

int VM2D::World2D::check02
mutable

Definition at line 109 of file World2D.h.

◆ checkPan

int VM2D::World2D::checkPan
mutable

Definition at line 110 of file World2D.h.

◆ cntrTreePnl

std::unique_ptr<CpuTreeInfo> VM2D::World2D::cntrTreePnl

Definition at line 116 of file World2D.h.

◆ cntrTreeVP

std::unique_ptr<CpuTreeInfo> VM2D::World2D::cntrTreeVP

Definition at line 115 of file World2D.h.

◆ cntrTreeWake

std::unique_ptr<CpuTreeInfo> VM2D::World2D::cntrTreeWake

Definition at line 114 of file World2D.h.

◆ cuda

Gpu VM2D::World2D::cuda
mutableprivate

Объект, управляющий графическим ускорителем

Definition at line 143 of file World2D.h.

◆ currentStep

size_t VMlib::WorldGen::currentStep
protectedinherited

Текущий номер шага в решаемой задаче

Definition at line 69 of file WorldGen.h.

◆ currentTime

double VMlib::WorldGen::currentTime
protectedinherited

Текущее время в решаемой задаче

Definition at line 72 of file WorldGen.h.

◆ dispBoundaryInSystem

std::vector<size_t> VM2D::World2D::dispBoundaryInSystem
private

Список номеров, с которых начинаются элементы правой части (или матрицы) системы для профилей

Definition at line 89 of file World2D.h.

◆ gabb

int VM2D::World2D::gabb
mutable

Definition at line 107 of file World2D.h.

◆ inflTreeWake

std::unique_ptr<CpuTreeInfo> VM2D::World2D::inflTreeWake

Деревья для быстрого метода

Definition at line 113 of file World2D.h.

◆ info

LogStream VMlib::WorldGen::info
mutableprotectedinherited

Поток для вывода логов и сообщений об ошибках

Definition at line 60 of file WorldGen.h.

◆ invMatr

Eigen::MatrixXd VM2D::World2D::invMatr
private

Обратная матрица

Definition at line 128 of file World2D.h.

◆ IQ

std::vector<std::vector<std::pair<Eigen::MatrixXd, Eigen::MatrixXd> > > VM2D::World2D::IQ
private

Матрица, состоящая из пар матриц, в которых хранятся касательные и нормальные компоненты интегралов от ядра

Definition at line 125 of file World2D.h.

◆ matrReord

Eigen::MatrixXd VM2D::World2D::matrReord
private

Матрица системы

Definition at line 122 of file World2D.h.

◆ measureVP

std::unique_ptr<MeasureVP> VM2D::World2D::measureVP
private

Умный указатель на алгоритм вычисления полей скоростей и давления (для сохранения в файл)

Definition at line 104 of file World2D.h.

◆ mechanics

std::vector<std::unique_ptr<Mechanics> > VM2D::World2D::mechanics
private

Список умных указателей на типы механической системы для каждого профиля

Definition at line 92 of file World2D.h.

◆ nVtxBeforeMerging

size_t VMlib::WorldGen::nVtxBeforeMerging
inherited

Definition at line 76 of file WorldGen.h.

◆ oldAirfoil

std::vector<std::unique_ptr<AirfoilGeometry> > VM2D::World2D::oldAirfoil
private

Список умных указателей на обтекаемые профили для сохранения старого положения

Definition at line 83 of file World2D.h.

◆ passport

const Passport& VM2D::World2D::passport
private

Константная ссылка на паспорт конкретного расчета

Definition at line 140 of file World2D.h.

◆ passportGen

const PassportGen& VMlib::WorldGen::passportGen
protectedinherited

Константная ссылка на паспорт конкретного расчета

Definition at line 63 of file WorldGen.h.

◆ rhsReord

Eigen::VectorXd VM2D::World2D::rhsReord
private

Правая часть системы

Definition at line 134 of file World2D.h.

◆ sol

Eigen::VectorXd VM2D::World2D::sol
private

Решение системы

Definition at line 137 of file World2D.h.

◆ source

std::unique_ptr<WakeDataBase> VM2D::World2D::source
private

Умный указатель на источники

Definition at line 101 of file World2D.h.

◆ timerConvVelo

VMlib::vmTimer VM2D::World2D::timerConvVelo

Definition at line 371 of file World2D.h.

◆ timerFillMatrix

VMlib::vmTimer VM2D::World2D::timerFillMatrix

Definition at line 368 of file World2D.h.

◆ timerInitialBuild

VMlib::vmTimer VM2D::World2D::timerInitialBuild

Definition at line 366 of file World2D.h.

◆ timerInside

VMlib::vmTimer VM2D::World2D::timerInside

Definition at line 372 of file World2D.h.

◆ timerMerging

VMlib::vmTimer VM2D::World2D::timerMerging

Definition at line 373 of file World2D.h.

◆ timerRhs

VMlib::vmTimer VM2D::World2D::timerRhs

Definition at line 367 of file World2D.h.

◆ timers

std::unique_ptr<TimersGen> VMlib::WorldGen::timers
protectedinherited

Сведения о временах выполнения основных операций

Definition at line 66 of file WorldGen.h.

◆ timerSlaeSolve

VMlib::vmTimer VM2D::World2D::timerSlaeSolve

Definition at line 369 of file World2D.h.

◆ useInverseMatrix

bool VM2D::World2D::useInverseMatrix
private

Признак использования обратной матрицы

Definition at line 131 of file World2D.h.

◆ velocity

std::unique_ptr<Velocity> VM2D::World2D::velocity
private

Умный укзатель на объект, определяющий методику вычисления скоростей

Definition at line 95 of file World2D.h.

◆ wake

std::unique_ptr<Wake> VM2D::World2D::wake
private

Умный указатель на вихревой след

Definition at line 98 of file World2D.h.


The documentation for this class was generated from the following files: