VM2D 1.14
Vortex methods for 2D flows simulation
Loading...
Searching...
No Matches
Gmres2D.cpp
Go to the documentation of this file.
1/*--------------------------------*- VM2D -*-----------------*---------------*\
2| ## ## ## ## #### ##### | | Version 1.14 |
3| ## ## ### ### ## ## ## ## | VM2D: Vortex Method | 2026/03/06 |
4| ## ## ## # ## ## ## ## | for 2D Flow Simulation *----------------*
5| #### ## ## ## ## ## | Open Source Code |
6| ## ## ## ###### ##### | https://www.github.com/vortexmethods/VM2D |
7| |
8| Copyright (C) 2017-2026 I. Marchevsky, K. Sokol, E. Ryatina, A. Kolganova |
9*-----------------------------------------------------------------------------*
10| File name: Gmres2D.cpp |
11| Info: Source code of VM2D |
12| |
13| This file is part of VM2D. |
14| VM2D is free software: you can redistribute it and/or modify it |
15| under the terms of the GNU General Public License as published by |
16| the Free Software Foundation, either version 3 of the License, or |
17| (at your option) any later version. |
18| |
19| VM2D is distributed in the hope that it will be useful, but WITHOUT |
20| ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or |
21| FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License |
22| for more details. |
23| |
24| You should have received a copy of the GNU General Public License |
25| along with VM2D. If not, see <http://www.gnu.org/licenses/>. |
26\*---------------------------------------------------------------------------*/
27
28
42#include "Gmres2D.h"
43
44#include "Airfoil2D.h"
45#include "Boundary2D.h"
46#include "MeasureVP2D.h"
47#include "Mechanics2D.h"
48#include "Preprocessor.h"
49#include "StreamParser.h"
50#include "Velocity2D.h"
51#include "Wake2D.h"
52#include "World2D.h"
53
54
55using namespace VM2D;
56
57typedef Eigen::Map<const Eigen::VectorXd> MapVec;
58
59bool VM2D::IterRot(std::vector<std::vector<double>>& H, const double nrmRhs, double& gs, std::vector<double>& c, std::vector<double>& s, int m, int n, double epsGMRES, int iter, bool residualShow)
60{
61 bool fl;
62 double buf;
63
64 for (int i = 0; i < m - 1; ++i)
65 {
66 buf = H[i][m - 1];
67
68 H[i][m - 1] = c[i] * buf + s[i] * H[i + 1][m - 1];
69 H[i + 1][m - 1] = -s[i] * buf + c[i] * H[i + 1][m - 1];
70 }
71
72 double izn = 1.0 / sqrt(H[m - 1][m - 1] * H[m - 1][m - 1] + H[m][m - 1] * H[m][m - 1]);
73
74 c.push_back(H[m - 1][m - 1] * izn);
75 s.push_back(H[m][m - 1] * izn);
76 gs *= -s[m - 1];
77
78
79
80 H[m - 1][m - 1] = c[m - 1] * H[m - 1][m - 1] + s[m - 1] * H[m][m - 1];
81 H[m][m - 1] = 0.0;
82
83
84
85 if (residualShow)
86 std::cout << "Iteration: " << iter << ", residual = " << (fabs(gs) / nrmRhs) << std::endl;
87
88 fl = ((fabs(gs) / nrmRhs) < epsGMRES);
89 return fl;
90}
91
92
93void VM2D::SolCircleRundirect(const std::vector<double>& A, const std::vector<double>& rhs, size_t startRow, size_t startRowReg, size_t np, std::vector<double>& res)
94{
95 size_t nAllVars = rhs.size();
96 std::vector<double> alpha(np + 2), beta(np + 2), gamma(np + 2);
97 std::vector<double> u(np), v(np);
98
99 double zn = A[(startRow + np + 1) * nAllVars + (startRow + np + 1)];
100 alpha[2] = -A[(startRow + np + 1) * nAllVars + (startRow + np + 2)] / zn;
101 beta[2] = rhs[startRow + np + 1] / zn;
102 gamma[2] = -A[(startRow + np + 1) * nAllVars + (startRow + np + 0)] / zn;
103
104 for (size_t i = 2; i <= np; ++i)
105 {
106 zn = A[(startRow + np + (i % np)) * nAllVars + (startRow + np + (i % np))] + alpha[i] * A[(startRow + np + (i % np)) * nAllVars + (startRow + np + ((i - 1) % np))];
107 alpha[i + 1] = -A[(startRow + np + (i % np)) * nAllVars + (startRow + np + ((i + 1) % np))] / zn;
108 beta[i + 1] = (rhs[startRow + np + (i % np)] - beta[i] * A[(startRow + np + (i % np)) * nAllVars + (startRow + np + ((i - 1) % np))]) / zn;
109 gamma[i + 1] = -(gamma[i] * A[(startRow + np + (i % np)) * nAllVars + (startRow + np + ((i - 1) % np))]) / zn;
110 }
111
112 u[np - 1] = beta[np];
113 v[np - 1] = alpha[np] + gamma[np];
114 for (size_t i = np - 2; i >= 1; --i)
115 {
116 u[i] = alpha[i + 1] * u[i + 1] + beta[i + 1];
117 v[i] = alpha[i + 1] * v[i + 1] + gamma[i + 1];
118 }
119
120 res[startRow + np] = (beta[np + 1] + alpha[np + 1] * u[1]) / (1.0 - gamma[np + 1] - alpha[np + 1] * v[1]);
121 for (size_t i = 1; i < np; ++i)
122 res[startRow + np + i] = u[i] + res[startRow + np] * v[i];
123
124 return;
125}
126
127
128
129
130void VM2D::SolMdirect(const std::vector<double>& A, const std::vector<double>& rhs, size_t startRow, size_t startRowReg, size_t np, std::vector<double>& res, bool lin)
131{
132 std::vector<double> alpha(np), beta(np), gamma(np), delta(np), phi(np), xi(np), psi(np);
133 double yn, ynn;
134 double lam1, lam2, mu1, mu2, xi1, xi2;
135 double a;
136
137 size_t nAllVars = rhs.size();
138
139 double izn = 1.0 / A[(startRow + 0) * nAllVars + (startRow + 0)];
140 alpha[1] = -A[(startRow + 0) * nAllVars + (startRow + 1)] * izn;
141 beta[1] = -A[(startRow + 0) * nAllVars + (startRow + np - 1)] * izn;
142 gamma[1] = -A[(startRow + 0) * nAllVars + (startRowReg)] * izn;
143 delta[1] = rhs[startRow] * izn;
144 double zn;
145
146
147 for (size_t i = 1; i < np - 1; ++i)
148 {
149 a = A[(startRow + i) * nAllVars + (startRow + i - 1)];
150 zn = a * alpha[i] + A[(startRow + i) * nAllVars + (startRow + i)];
151 alpha[i + 1] = -A[(startRow + i) * nAllVars + (startRow + i + 1)] / zn;
152 beta[i + 1] = -a * beta[i] / zn;
153 gamma[i + 1] = -(A[(startRow + i) * nAllVars + (startRowReg)] + a * gamma[i]) / zn;
154 delta[i + 1] = (rhs[startRow + i] - a * delta[i]) / zn;
155 }
156
157 a = A[(startRow + (np - 2)) * nAllVars + (startRow + (np - 3))];
158 zn = alpha[np - 2] * a + A[(startRow + (np - 2)) * nAllVars + (startRow + (np - 2))];
159
160 phi[np - 1] = -(A[(startRow + np - 2) * nAllVars + (startRow + np - 1)] + a * beta[np - 2]) / zn;
161 psi[np - 1] = -(A[(startRow + np - 2) * nAllVars + (startRowReg)] + gamma[np - 2] * a) / zn;
162 xi[np - 1] = (rhs[startRow + np - 2] - delta[np - 2] * a) / zn;
163
164 for (size_t i = np - 2; i > 0; --i)
165 {
166 phi[i] = alpha[i] * phi[i + 1] + beta[i];
167 psi[i] = alpha[i] * psi[i + 1] + gamma[i];
168 xi[i] = alpha[i] * xi[i + 1] + delta[i];
169 }
170
171
172 double e = A[(startRow + np - 1) * nAllVars + (startRow + 0)];
173 a = A[(startRow + np - 1) * nAllVars + (startRow + np - 2)];
174 lam1 = e * phi[1] + a * phi[np - 1] + A[(startRow + np - 1) * nAllVars + (startRow + np - 1)];
175 mu1 = e * psi[1] + a * psi[np - 1] + A[(startRow + np - 1) * nAllVars + (startRowReg)];
176 xi1 = rhs[startRow + np - 1] - e * xi[1] - a * xi[np - 1];
177
178 lam2 = mu2 = xi2 = 0.0;
179 for (size_t j = 0; j < np - 1; ++j)
180 {
181 lam2 += phi[j + 1];
182 mu2 += psi[j + 1];
183 xi2 -= xi[j + 1];
184 }
185 lam2 += 1.0;
186
187 xi2 += rhs[startRowReg];
188
189 zn = lam1 * mu2 - lam2 * mu1;
190 yn = (xi1 * mu2 - xi2 * mu1) / zn;
191 ynn = -(xi1 * lam2 - xi2 * lam1) / zn;
192
193 res[startRowReg] = ynn;
194
195 res[startRow + np - 1] = yn;
196 for (size_t i = np - 2; i + 1 > 0; --i)
197 res[startRow + i] = phi[i + 1] * yn + psi[i + 1] * ynn + xi[i + 1];
198
199//
200 if (lin)
201 SolCircleRundirect(A, rhs, startRow, startRowReg, np, res);
202}
203
204
205
206void VM2D::GMRES_Direct(
207 const World2D& W,
208 int nAllVars,
209 int nafl,
210 const std::vector<double>& mtrDir,
211 const std::vector<double>& rhsDir,
212 const std::vector<int>& pos,
213 const std::vector<int>& vsize,
214 std::vector<std::vector<double>>& gam,
215 std::vector<double>& R)
216{
217 const size_t maxIter = nAllVars + 1;
218 std::vector<double> x(nAllVars, 0.0);
219 std::vector<double> y(nAllVars, 0.0);
220 std::vector<double> residual(nAllVars, 0.0);
221 std::vector<double> r(nAllVars, 0.0);
222 std::vector<double> w(nAllVars, 0.0);
223 std::vector<double> wOld(nAllVars, 0.0);
224 std::vector<double> g(nAllVars, 0.0);
225 std::vector<double> c(nAllVars, 0.0);
226 std::vector<double> s(nAllVars, 0.0);
227 std::vector<double> hisGs(nAllVars + 1, 0.0);
228 std::vector<double> v(nAllVars * (maxIter + 1), 0.0);
229
230 std::vector<double> vecFromV(nAllVars, 0.0);
231
232 std::vector<double> h((nAllVars + 1) * maxIter, 0.0);
233 double beta = 0.0, gs, normb = 0.0;
234 size_t m = 0;
235
236 for (size_t i = 0; i < nAllVars; ++i)
237 normb += rhsDir[i] * rhsDir[i];
238 normb = sqrt(normb);
239
240 for (size_t i = 0; i < nAllVars; ++i)
241 residual[i] = rhsDir[i];
242
243#pragma omp parallel for
244 for (int aflI = 0; aflI < nafl; ++aflI)
245 {
246 // SolMdirect(mtrDir, residual, pos[aflI], nAllVars - nafl + aflI, W.getAirfoil(aflI).getNumberOfPanels(), r, (vsize[aflI] == 2*n[aflI]));
247 r = residual;
248 }
249
250 for (size_t i = 0; i < nAllVars; ++i)
251 beta += r[i] * r[i];
252 beta = sqrt(beta);
253 gs = beta;
254 g[0] = beta;
255
256 // ,
257 if (beta == 0)
258 {
259 for (size_t aflI = 0; aflI < nafl; ++aflI)
260 {
261 for (size_t i = 0; i < vsize[aflI]; ++i)
262 gam[aflI][i] = 0.0;
263 R[aflI] = 0.0;
264 }
265 return;
266 }
267
268 for (size_t aflI = 0; aflI < nafl; ++aflI)
269 R[aflI] = x[x.size() - nafl + aflI];
270
271 for (size_t i = 0; i < nAllVars; ++i)
272 v[i * (maxIter + 1) + 0] = r[i] / beta;
273
274 for (size_t j = 0; j < nAllVars; ++j)
275 {
276 double timer1 = omp_get_wtime();
277
278 for (size_t p = 0; p < nAllVars; ++p)
279 vecFromV[p] = v[p * (maxIter + 1) + j];
280
281 //cblas_dgemv(CblasRowMajor, CblasNoTrans, nAllVars, nAllVars, 1.0, mtrDir.data(), nAllVars, &v[0 * (maxIter + 1) + j], (nAllVars + 1), 0.0, wOld.data(), 1);
282 //cblas_dgemv(CblasRowMajor, CblasNoTrans, nAllVars, nAllVars, 1.0, mtrDir.data(), nAllVars, vecFromV.data(), 1, 0.0, wOld.data(), 1);
283 for (int i = 0; i < nAllVars; ++i)
284 {
285 wOld[i] = 0.0;
286 for (int jj = 0; jj < nAllVars; ++jj)
287 {
288 wOld[i] += mtrDir[i * nAllVars + jj] * vecFromV[jj];
289 }
290 }
291
292 double timer2 = omp_get_wtime();
293
294#pragma omp parallel for
295 for (int aflI = 0; aflI < nafl; ++aflI)
296 {
297 // SolMdirect(mtrDir, wOld, pos[aflI], nAllVars - nafl + aflI, W.getAirfoil(aflI).getNumberOfPanels(), w, (vsize[aflI] == 2 * n[aflI]));
298 w = wOld;
299 }
300
301 for (size_t i = 0; i <= j; ++i)
302 {
303 h[i * maxIter + j] = 0.0;
304 for (size_t q = 0; q < nAllVars; ++q)
305 h[i * maxIter + j] += w[q] * v[q * (maxIter + 1) + i];
306
307 for (size_t q = 0; q < nAllVars; ++q)
308 w[q] -= h[i * maxIter + j] * v[q * (maxIter + 1) + i];
309 }//for i
310
311 h[(j + 1) * maxIter + j] = 0.0;
312 for (size_t p = 0; p < nAllVars; ++p)
313 h[(j + 1) * maxIter + j] += w[p] * w[p];
314 h[(j + 1) * maxIter + j] = sqrt(h[(j + 1) * maxIter + j]);
315
316 for (size_t p = 0; p < nAllVars; ++p)
317 v[p * (maxIter + 1) + (j + 1)] = w[p] / h[(j + 1) * maxIter + j];
318
319 for (size_t i = 0; i < j; ++i)
320 {
321 double buf = h[i * maxIter + j];
322 h[i * maxIter + j] = c[i] * buf + s[i] * h[(i + 1) * maxIter + j];
323 h[(i + 1) * maxIter + j] = -s[i] * buf + c[i] * h[(i + 1) * maxIter + j];
324 }
325
326 double zn = sqrt(h[j * maxIter + j] * h[j * maxIter + j] + h[(j + 1) * maxIter + j] * h[(j + 1) * maxIter + j]);
327 c[j] = fabs(h[j * maxIter + j] / zn);
328 s[j] = c[j] * h[(j + 1) * maxIter + j] / h[j * maxIter + j];
329
330 gs *= -s[j];
331
332 h[j * maxIter + j] = c[j] * h[j * maxIter + j] + s[j] * h[(j + 1) * maxIter + j];
333 h[(j + 1) * maxIter + j] = 0.0;
334 hisGs[j] = fabs(gs) / normb;
335
336 if ((fabs(gs) / normb) < W.getPassport().numericalSchemes.gmresEps)
337 {
338 m = j + 1;
339 break;
340 }
341 }//for j
342
343 //
344 for (size_t i = 0; i < m; ++i)
345 {
346 double buf = g[i];
347 g[i] *= c[i];
348 if (i < (m - 1))
349 g[i + 1] = -s[i] * buf;
350
351 }
352 y[m - 1] = g[m - 1] / h[(m - 1) * maxIter + (m - 1)];
353
354 for (size_t i = m - 2; i + 1 > 0; --i)
355 {
356 double sum = 0;
357 for (size_t j = i + 1; j < m; ++j)
358 sum += h[i * maxIter + j] * y[j];
359
360 y[i] = (g[i] - sum) / h[i * maxIter + i];
361 }
362
363 for (size_t p = 0; p < nAllVars; ++p)
364 {
365 for (size_t q = 0; q < m; ++q)
366 x[p] += v[p * (maxIter + 1) + q] * y[q];
367 }
368
369 std::cout << "#iterations: " << m << std::endl;
370
371 for (size_t aflI = 0; aflI < nafl; ++aflI)
372 {
373 for (size_t i = 0; i < vsize[aflI]; ++i)
374 gam[aflI][i] = x[pos[aflI] + i];
375 }
376
377 for (size_t aflI = 0; aflI < nafl; ++aflI)
378 R[aflI] = x[x.size() - nafl + aflI];
379}
380
381
382
383
384
385
386
387
388
390
391
392
393void VM2D::SolCircleRun(std::vector<double>& AX, const std::vector<double>& rhs, const Airfoil& afl, const std::vector<std::vector<Point2D>>& prea1, const std::vector<std::vector<Point2D>>& prec1, int p, int n, sweepVectors& sw)
394{
395 //std::vector<double> alpha(n), beta(n), delta(n), phi(n), xi(n);
396 std::vector<double>& alpha = sw.alpha;
397 std::vector<double>& beta = sw.beta;
398 std::vector<double>& delta = sw.delta;
399 std::vector<double>& phi = sw.phi;
400 std::vector<double>& xi = sw.xi;
401
402 double yn, a;
403
404 double len24 = afl.len[0] * 24.0;
405
406 alpha[1] = (afl.tau[0] & prec1[p][0]) * len24;
407 beta[1] = (afl.tau[0] & prea1[p][0]) * len24;
408 delta[1] = -len24 * rhs[n + 0];
409 double zn;
410
411 for (int i = 1; i < n - 1; ++i)
412 {
413 a = (afl.tau[i] & prea1[p][i]);
414 zn = a * alpha[i] - 1.0 / (24.0 * afl.len[i]);
415 alpha[i + 1] = -(afl.tau[i] & prec1[p][i]) / zn;
416 beta[i + 1] = -a * beta[i] / zn;
417 delta[i + 1] = (rhs[n + i] - a * delta[i]) / zn;
418 }
419
420 a = afl.tau[n - 2] & prea1[p][n - 2];
421 zn = alpha[n - 2] * a - 1.0 / (24.0 * afl.len[n - 2]);
422 phi[n - 1] = -((afl.tau[n - 2] & prec1[p][n - 2]) + beta[n - 2] * a) / zn;
423 xi[n - 1] = (rhs[n + n - 2] - delta[n - 2] * a) / zn;
424
425 for (int i = n - 2; i > 0; --i)
426 {
427 phi[i] = alpha[i] * phi[i + 1] + beta[i];
428 xi[i] = alpha[i] * xi[i + 1] + delta[i];
429 }
430
431
432 double e = (afl.tau[n - 1] & prec1[p][n - 1]);
433 a = (afl.tau[n - 1] & prea1[p][n - 1]);
434
435 yn = (rhs[n + n - 1] - e * xi[1] - a * xi[n - 1]) / (e * phi[1] + a * phi[n - 1] - 1.0 / (24.0 * afl.len[n - 1]));
436
437 AX[n + n - 1] = yn;
438 for (int i = n - 2; i >= 0; --i)
439 AX[n + i] = phi[i + 1] * yn + xi[i + 1];
440}
441
442
443
444void VM2D::SolM(std::vector<double>& AX, const std::vector<double>& rhs, const Airfoil& afl,
445 const std::vector<std::vector<Point2D>>& prea, const std::vector<std::vector<Point2D>>& prec, int p, int n, bool linScheme, sweepVectors& sw)
446{
447 //std::vector<double> alpha(n), beta(n), gamma(n), delta(n), phi(n), xi(n), psi(n);
448 std::vector<double>& alpha = sw.alpha;
449 std::vector<double>& beta = sw.beta;
450 std::vector<double>& gamma = sw.gamma;
451 std::vector<double>& delta = sw.delta;
452 std::vector<double>& phi = sw.phi;
453 std::vector<double>& xi = sw.xi;
454 std::vector<double>& psi = sw.psi;
455
456 double yn, ynn;
457 double lam1, lam2, mu1, mu2, xi1, xi2;
458 double a;
459
460 double len2 = afl.len[0] * 2.0;
461
462 alpha[1] = (afl.tau[0] & prec[p][0]) * len2;
463 beta[1] = (afl.tau[0] & prea[p][0]) * len2;
464 gamma[1] = len2;
465 delta[1] = -len2 * rhs[0];
466
467 double zn;
468
469 for (int i = 1; i < n - 1; ++i)
470 {
471 a = (afl.tau[i] & prea[p][i]);
472 zn = a * alpha[i] - 0.5 / afl.len[i];
473 alpha[i + 1] = -(afl.tau[i] & prec[p][i]) / zn;
474 beta[i + 1] = -a * beta[i] / zn;
475 gamma[i + 1] = -(1.0 + a * gamma[i]) / zn;
476 delta[i + 1] = (rhs[i] - a * delta[i]) / zn;
477 }
478
479 a = afl.tau[n - 2] & prea[p][n - 2];
480 zn = alpha[n - 2] * a - 0.5 / afl.len[n - 2];
481 phi[n - 1] = -((afl.tau[n - 2] & prec[p][n - 2]) + beta[n - 2] * a) / zn;
482 psi[n - 1] = -(1.0 + gamma[n - 2] * a) / zn;
483 xi[n - 1] = (rhs[n - 2] - delta[n - 2] * a) / zn;
484
485 for (int i = n - 2; i > 0; --i)
486 {
487 phi[i] = alpha[i] * phi[i + 1] + beta[i];
488 psi[i] = alpha[i] * psi[i + 1] + gamma[i];
489 xi[i] = alpha[i] * xi[i + 1] + delta[i];
490 }
491
492 double e = (afl.tau[n - 1] & prec[p][n - 1]);
493 a = (afl.tau[n - 1] & prea[p][n - 1]);
494 lam1 = e * phi[1] + a * phi[n - 1] - 0.5 / afl.len[n - 1];
495 mu1 = e * psi[1] + a * psi[n - 1] + 1.0;
496 xi1 = rhs[n - 1] - e * xi[1] - a * xi[n - 1];
497 lam2 = mu2 = xi2 = 0.0;
498 for (int j = 0; j < n - 1; ++j)
499 {
500 lam2 += phi[j + 1];
501 mu2 += psi[j + 1];
502 xi2 -= xi[j + 1];
503 }
504 lam2 += 1.0;
505
506 if (!linScheme)
507 xi2 += rhs[n];
508 else if (linScheme)
509 xi2 += rhs[2 * n];
510
511 zn = lam1 * mu2 - lam2 * mu1;
512 yn = (xi1 * mu2 - xi2 * mu1) / zn;
513 ynn = -(xi1 * lam2 - xi2 * lam1) / zn;
514
515 if (!linScheme)
516 AX[n] = ynn;
517 else if (linScheme)
518 AX[2 * n] = ynn;
519
520 AX[n - 1] = yn;
521 for (int i = n - 2; i >= 0; --i)
522 AX[i] = phi[i + 1] * yn + psi[i + 1] * ynn + xi[i + 1];
523}
524
525
526
527//
528void PreCalculateCoef(
529const Airfoil& aflP,
530size_t nPanelsP,
531std::vector<Point2D> & prea,
532std::vector<Point2D>& prec,
533std::vector<Point2D>& prea1, //
534std::vector<Point2D>& prec1, //
535bool linScheme)
536{
537 Point2D p1, s1, p2, s2, di, dj;
538 numvector<double, 3> alpha, lambda;
539 numvector<Point2D, 3> v00, v11;
540 Point2D i00, i11;
541
542 for (size_t i = 0; i < nPanelsP; ++i)
543 {
544 std::array<int, 2> vecV = { ((int)i - 1 >= 0) ? int(i - 1) : int(nPanelsP - 1), (i + 1 < nPanelsP) ? int(i + 1) : 0 };
545 for (int j : vecV)
546 {
547 p1 = aflP.getR(i + 1) - aflP.getR(j + 1);
548 s1 = aflP.getR(i + 1) - aflP.getR(j);
549 p2 = aflP.getR(i) - aflP.getR(j + 1);
550 s2 = aflP.getR(i) - aflP.getR(j);
551 di = aflP.getR(i + 1) - aflP.getR(i);
552 dj = aflP.getR(j + 1) - aflP.getR(j);
553
554 alpha = { \
555 (aflP.isAfter(j, i)) ? 0.0 : VMlib::Alpha(s2, s1), \
556 VMlib::Alpha(s2, p1), \
557 (aflP.isAfter(i, j)) ? 0.0 : VMlib::Alpha(p1, p2) \
558 };
559
560 lambda = { \
561 (aflP.isAfter(j, i)) ? 0.0 : VMlib::Lambda(s2, s1), \
562 VMlib::Lambda(s2, p1), \
563 (aflP.isAfter(i, j)) ? 0.0 : VMlib::Lambda(p1, p2) \
564 };
565
566 v00 = {
567 VMlib::Omega(s1, di.unit(), dj.unit()),
568 -VMlib::Omega(di, di.unit(), dj.unit()),
569 VMlib::Omega(p2, di.unit(), dj.unit())
570 };
571
572 i00 = IDPI / di.length() * (-(alpha[0] * v00[0] + alpha[1] * v00[1] + alpha[2] * v00[2]).kcross() \
573 + (lambda[0] * v00[0] + lambda[1] * v00[1] + lambda[2] * v00[2]));
574
575 if (aflP.isAfter(j, i))
576 prec[i] = (i00.kcross()) * (1.0 / dj.length());
577 if (aflP.isAfter(i, j))
578 prea[i] = (i00.kcross()) * (1.0 / dj.length());
579
580 if (linScheme) {
581
582 v11 = { 1.0 / (12.0 * di.length() * dj.length()) * \
583 (2.0 * (s1 & Omega(s1 - 3.0 * p2, di.unit(), dj.unit())) * VMlib::Omega(s1, di.unit(), dj.unit()) - \
584 s1.length2() * (s1 - 3.0 * p2)) - 0.25 * VMlib::Omega(s1, di.unit(), dj.unit()),
585 -di.length() / (12.0 * dj.length()) * Omega(di, dj.unit(), dj.unit()),
586 -dj.length() / (12.0 * di.length()) * Omega(dj, di.unit(), di.unit()) };
587
588 i11 = (IDPI / di.length()) * ((alpha[0] + alpha[2]) * v11[0] + (alpha[1] + alpha[2]) * v11[1] + alpha[2] * v11[2]\
589 + ((lambda[0] + lambda[2]) * v11[0] + (lambda[1] + lambda[2]) * v11[1] + lambda[2] * v11[2] \
590 + 1.0 / 12.0 * (dj.length() * di.unit() + di.length() * dj.unit() - 2.0 * VMlib::Omega(s1, di.unit(), dj.unit()))).kcross());
591 if (aflP.isAfter(j, i))
592 prec1[i] = (i11) * (1.0 / dj.length());
593 if (aflP.isAfter(i, j))
594 prea1[i] = (i11) * (1.0 / dj.length());
595 }
596 }
597 }
598}
599
600
601#ifdef USE_CUDA
602#include <cublas_v2.h>
603
604void VM2D::GMRES(const World2D& W,
605 std::vector<std::vector<double>>& X,
606 std::vector<double>& R,
607 const std::vector<std::vector<double>>& rhs,
608 const std::vector<double> rhsReg,
609 int& niter,
610 bool linScheme)
611{
612 VMlib::vmTimer tAll("All"), tPre("Pre"), tIter("Iter"), tPost("Post"), tIter1("Iter1"), tIterOther("IterOther"),
613 tCG("tCG"), tGC("tGC"), tWrapper("tWrapper"), tSplit("tSplit"), tPrecond("tPrecond"), tRot("tRot"),
614 tRotA("tRotA"), tRotB("tRotB"), tRotC("tRotC"), tRotD("tRotD");
615
616 cublasHandle_t cublas_handle;
617 cublasCreate(&cublas_handle);
618
619 const size_t nAfl = W.getNumberOfAirfoil();
620
621 size_t totalVsize = 0;
622 for (size_t i = 0; i < nAfl; ++i)
623 totalVsize += W.getAirfoil(i).getNumberOfPanels();
624
625 size_t nTotPan = 0;
626 for (size_t s = 0; s < nAfl; ++s)
627 nTotPan += W.getAirfoil(s).getNumberOfPanels();
628
629 if (linScheme)
630 totalVsize *= 2;
631
632
633 const size_t iterSize = 50;
634
635 double* devViBund, *devw;
636 cudaMalloc((void**)&devViBund, (totalVsize + nAfl) * sizeof(double) * iterSize);
637 cudaMalloc((void**)&devw, (totalVsize + nAfl) * sizeof(double));
638
639 double* mulptr;
640 cudaMalloc(&mulptr, iterSize * sizeof(double));
641 std::vector<double> mul(iterSize);
642
643 tAll.start();
644 tPre.start();
645
646 size_t m;
647
648 double beta;
649 std::vector<double> w(totalVsize + nAfl), c, s;
650 c.reserve(iterSize);
651 s.reserve(iterSize);
652
653 std::vector<double> Vflat;
654 Vflat.reserve((totalVsize + nAfl) * iterSize);
655
656 std::vector<std::vector<double>> H, diag(nAfl);
657 H.resize(iterSize + 1);
658 for (int i = 0; i < iterSize + 1; ++i)
659 H[i].resize(iterSize);
660
661 for (int p = 0; p < (int)nAfl; ++p)
662 {
663 size_t nPanelsP = W.getAirfoil(p).getNumberOfPanels();
664 if (!linScheme)
665 diag[p].resize(nPanelsP);
666 else
667 diag[p].resize(2 * nPanelsP);
668
669 for (int i = 0; i < nPanelsP; ++i)
670 {
671 double lenI = W.getAirfoil(p).len[i];
672 diag[p][i] = 0.5 / lenI;
673
674 if (linScheme)
675 diag[p][i + nPanelsP] = (1.0 / 24.0) / lenI;
676 }
677 }
678
679 double nrmRhs = 0.0;
680#pragma omp parallel for reduction(+:nrmRhs)
681 for (int p = 0; p < (int)nAfl; ++p)
682 nrmRhs += norm2(rhs[p]);
683 nrmRhs = sqrt(nrmRhs);
684
685 std::vector<std::vector<double>> residual(nAfl);
686 for (size_t p = 0; p < nAfl; ++p)
687 residual[p] = rhs[p];
688
689 for (size_t i = 0; i < nAfl; ++i)
690 Vflat.insert(Vflat.end(), residual[i].begin(), residual[i].end());
691
692 for (size_t i = 0; i < nAfl; ++i)
693 Vflat.push_back(rhsReg[i]);
694
695 // PRECONDITIONER
696 std::vector<std::vector<Point2D>> prea(nAfl), prec(nAfl), prea1(nAfl), prec1(nAfl);
697
698 for (size_t p = 0; p < nAfl; ++p)
699 {
700 const auto& aflP = W.getAirfoil(p);
701 size_t nPanelsP = aflP.getNumberOfPanels();
702
703 prea[p].resize(nPanelsP);
704 prec[p].resize(nPanelsP);
705 prea1[p].resize(nPanelsP);
706 prec1[p].resize(nPanelsP);
707
708 PreCalculateCoef(aflP, nPanelsP, prea[p], prec[p], prea1[p], prec1[p], linScheme);
709 }
710
711 std::vector<sweepVectors> allSW(nAfl);
712 std::vector<std::vector<double>> allBuf2(nAfl);
713 std::vector<size_t> np(nAfl, 0);
714
715 for (size_t p = 0; p < nAfl; ++p)
716 {
717 size_t nPanelsP = W.getAirfoil(p).getNumberOfPanels();
718
719 allSW[p].resize(nPanelsP);
720
721 if (!linScheme)
722 allBuf2[p].resize(nPanelsP + 1);
723 else
724 allBuf2[p].resize(2 * nPanelsP + 1);
725
726 if (p == 0)
727 continue;
728
729 np[p] = np[p - 1];
730 if (!linScheme)
731 np[p] += (/*(p == 0) ? 0 :*/ W.getAirfoil(p - 1).getNumberOfPanels());
732 else
733 np[p] += (/*(p == 0) ? 0 :*/ 2 * W.getAirfoil(p - 1).getNumberOfPanels());
734 }
735
736#pragma omp parallel for
737 for (int p = 0; p < (int)nAfl; ++p)
738 {
739 std::vector<double>& vbuf1 = allBuf2[p];
740 size_t nPanelsP = W.getAirfoil(p).getNumberOfPanels();
741
742 for (size_t i = 0; i < nPanelsP; ++i)
743 {
744 vbuf1[i] = Vflat[np[p] + i];
745 if (linScheme)
746 vbuf1[i + nPanelsP] = Vflat[np[p] + nPanelsP + i];
747 }
748
749 if (!linScheme)
750 vbuf1[nPanelsP] = Vflat[(totalVsize + nAfl) - (nAfl - p)];
751 else
752 vbuf1[2 * nPanelsP] = Vflat[(totalVsize + nAfl) - (nAfl - p)];
753
754
755 SolM(vbuf1, vbuf1, W.getAirfoil(p), prea, prec, (int)p, (int)nPanelsP, linScheme, allSW[p]);
756
757 //
758 if (linScheme)
759 SolCircleRun(vbuf1, vbuf1, W.getAirfoil(p), prea1, prec1, p, (int)nPanelsP, allSW[p]);
760
761
762 for (size_t j = np[p]; j < np[p] + nPanelsP; ++j)
763 {
764 Vflat[j] = vbuf1[j - np[p]];
765 if (linScheme)
766 Vflat[j + nPanelsP] = vbuf1[j - np[p] + nPanelsP];
767 }
768 Vflat[(totalVsize + nAfl) - (nAfl - p)] = vbuf1[vbuf1.size() - 1];
769 }
770
771 beta = 0;
772 for (int t = 0; t < (totalVsize + nAfl); ++t)
773 beta += Vflat[t] * Vflat[t];
774 beta = sqrt(beta);
775
776
777 if (beta > 0)
778 for (int t = 0; t < (totalVsize + nAfl); ++t)
779 Vflat[t] /= beta;
780
781 else
782 {
783 for (auto& xp : X)
784 for (auto& x : xp)
785 x = 0.0;
786
787 cublasDestroy(cublas_handle);
788 return;
789 }
790
791 double gs = beta;
792
793 std::vector<double> bufnewSol((linScheme ? 2 : 1) * nTotPan), bufcurrentSol(totalVsize, 0.0);
794
795 auto& treePnlInfl = *W.getNonConstCuda().inflTreePnlVortex;
796 treePnlInfl.MemoryAllocateForGMRES();
797
798 W.getNonConstCuda().AllocateSolution(W.getNonConstCuda().dev_sol, nTotPan);
799 if (linScheme)
800 W.getNonConstCuda().AllocateSolution(W.getNonConstCuda().dev_solLin, nTotPan);
801 else
802 W.getNonConstCuda().dev_solLin = nullptr;
803
804 double* dev_ptr_rhs = W.getAirfoil(0).devRhsPtr;
805 double* dev_ptr_rhsLin = (linScheme ? W.getAirfoil(0).devRhsLinPtr : nullptr);
806
807 const double zero = 0.0, unit = 1.0;
808 tPre.stop();
809
810 tIter.start();
811
812
813 //Iterations
814 for (int j = 0; j < totalVsize - 1; ++j) //+ n.size()
815 {
816 if (j == 0)
817 tIter1.start();
818 else if (j == 1)
819 tIterOther.start();
820
821 tCG.start();
822
823 size_t npred = 0;
824 if (!linScheme)
825 for (size_t i = 0; i < W.getNumberOfAirfoil(); ++i)
826 {
827#pragma omp parallel for
828 for (int p = 0; p < (int)W.getAirfoil(i).getNumberOfPanels(); ++p)
829 bufcurrentSol[npred + p] = Vflat[(totalVsize + nAfl)*(j) + npred + p] / W.getAirfoil(i).len[p];
830
831 npred += W.getAirfoil(i).getNumberOfPanels();
832 }
833 else
834 {
835 for (size_t i = 0; i < W.getNumberOfAirfoil(); ++i)
836 {
837#pragma omp parallel for
838 for (int p = 0; p < (int)W.getAirfoil(i).getNumberOfPanels(); ++p)
839 {
840 bufcurrentSol[npred + p] = Vflat[(totalVsize + nAfl) * j + (2 * npred + p)] / W.getAirfoil(i).len[p];
841 bufcurrentSol[nTotPan + npred + p] = Vflat[(totalVsize + nAfl) * j + (2 * npred + W.getAirfoil(i).getNumberOfPanels() + p)] / W.getAirfoil(i).len[p];
842 }
843 npred += W.getAirfoil(i).getNumberOfPanels();
844 }
845 W.getNonConstCuda().SetSolution(bufcurrentSol.data() + nTotPan, W.getNonConstCuda().dev_solLin, nTotPan);
846 }
847 W.getNonConstCuda().SetSolution(bufcurrentSol.data(), W.getNonConstCuda().dev_sol, nTotPan);
848
849 tCG.stop();
850
851 if (j>0)
852 tWrapper.start();
853
854 int order = W.getPassport().numericalSchemes.gmresMultipoleOrder;
855 double theta = W.getPassport().numericalSchemes.gmresTheta;
856
857 treePnlInfl.UpdatePanelFreeVortexIntensity(W.getNonConstCuda().dev_sol, W.getNonConstCuda().dev_solLin);
858 treePnlInfl.UpwardTraversal(order);
859 treePnlInfl.DownwardTraversalGMRES(dev_ptr_rhs, dev_ptr_rhsLin, theta, order, j);
860
861 if (j>0)
862 tWrapper.stop();
863
864 tGC.start();
865
866 if (!linScheme)
867 W.getCuda().CopyMemFromDev<double, 1>(nTotPan, dev_ptr_rhs, w.data(), 22);
868
869 if (linScheme)
870 {
871 W.getCuda().CopyMemFromDev<double, 1>(nTotPan, dev_ptr_rhs, bufnewSol.data(), 22);
872 W.getCuda().CopyMemFromDev<double, 1>(nTotPan, dev_ptr_rhsLin, bufnewSol.data() + nTotPan, 22);
873 }
874
875 tGC.stop();
876
877 tSplit.start();
878
879 npred = 0;
880 if (linScheme)
881 for (size_t i = 0; i < W.getNumberOfAirfoil(); ++i)
882 {
883#pragma omp parallel for
884 for (int j = 0; j < (int)W.getAirfoil(i).getNumberOfPanels(); ++j)
885 {
886 w[2 * npred + j] = bufnewSol[npred + j];
887 w[2 * npred + W.getAirfoil(i).getNumberOfPanels() + j] = bufnewSol[nTotPan + npred + j];
888 }
889 npred += W.getAirfoil(i).getNumberOfPanels();
890 }
891
892 size_t cntr = 0;
893 for (size_t pi = 0; pi < nAfl; ++pi)
894 {
895 size_t nPanelsPi = W.getAirfoil(pi).getNumberOfPanels();
896#pragma omp parallel for
897 for (int i = 0; i < (int)nPanelsPi; ++i)
898 {
899 w[cntr + i] += Vflat[(totalVsize + nAfl) * j + (totalVsize + pi)] - Vflat[(totalVsize + nAfl) * j + (cntr + i)] * diag[pi][i];
900 if (linScheme)
901 w[cntr + i + nPanelsPi] += -Vflat[(totalVsize + nAfl) * j + (cntr + i + nPanelsPi)] * diag[pi][i + nPanelsPi];
902 }
903
904 cntr += nPanelsPi;
905 if (linScheme)
906 cntr += nPanelsPi;
907 }
908
909 for (size_t i = 0; i < nAfl; ++i)
910 w[totalVsize + i] = 0.0;
911
912 cntr = 0;
913 for (size_t i = 0; i < nAfl; ++i)
914 {
915 size_t nPanelsI = W.getAirfoil(i).getNumberOfPanels();
916 for (size_t k = 0; k < nPanelsI; ++k)
917 w[totalVsize + i] += Vflat[(totalVsize + nAfl) * j + (cntr + k)];
918 cntr += nPanelsI;
919 if (linScheme)
920 cntr += nPanelsI;
921 }
922 tSplit.stop();
923
924 tPrecond.start();
925
926 // PRECONDITIONER
927
928#pragma omp parallel for
929 for (int p = 0; p < (int)nAfl; ++p)
930 {
931 std::vector<double>& vbuf2 = allBuf2[p];
932
933 size_t nPanelsP = W.getAirfoil(p).getNumberOfPanels();
934
935 for (size_t i = 0; i < nPanelsP; ++i)
936 {
937 vbuf2[i] = w[np[p] + i];
938 if (linScheme)
939 vbuf2[i + nPanelsP] = w[np[p] + nPanelsP + i];
940 }
941
942 if (!linScheme)
943 vbuf2[nPanelsP] = w[w.size() - (nAfl - p)];
944 else
945 vbuf2[2 * nPanelsP] = w[w.size() - (nAfl - p)];
946
947 SolM(vbuf2, vbuf2, W.getAirfoil(p), prea, prec, (int)p, (int)nPanelsP, linScheme, allSW[p]);
948
949 //
950 if (linScheme)
951 SolCircleRun(vbuf2, vbuf2, W.getAirfoil(p), prea1, prec1, p, (int)nPanelsP, allSW[p]);
952
953
954 for (size_t j = np[p]; j < np[p] + nPanelsP; ++j)
955 {
956 w[j] = vbuf2[j - np[p]];
957 if (linScheme)
958 w[j + nPanelsP] = vbuf2[j - np[p] + nPanelsP];
959 }
960
961 w[w.size() - (nAfl - p)] = vbuf2[vbuf2.size() - 1];
962 }
963 tPrecond.stop();
964
965
966
967 tRot.start();
968
969 tRotA.start();
970
971 if (j == H.size() - 1)
972 {
973 H.resize(j * 2 + 2);
974 for (int i = 0; i < j * 2 + 2; ++i)
975 H[i].resize(j * 2 + 1);
976
977 Vflat.reserve((totalVsize + nAfl) * j * 2);
978 c.reserve(j * 2);
979 s.reserve(j * 2);
980
981 mul.reserve(j * 2);
982 cudaFree(mulptr);
983 cudaMalloc((void**)&mulptr, j * 2 * sizeof(double));
984
985 double* devViBund_new;
986 cudaMalloc((void**)&devViBund_new, (totalVsize + nAfl) * j * 2 * sizeof(double));
987 cudaMemcpy(devViBund_new, devViBund, (totalVsize + nAfl) * j * sizeof(double), cudaMemcpyDeviceToDevice);
988 cudaFree(devViBund);
989 devViBund = devViBund_new;
990 }
991
992 tRotA.stop();
993
994 tRotB.start();
995 cudaMemcpy(devViBund + w.size() * j, Vflat.data() + w.size() * j, (1) * w.size() * sizeof(double), cudaMemcpyHostToDevice); //was
996 cudaMemcpy(devw, w.data(), w.size() * sizeof(double), cudaMemcpyHostToDevice); //was
997
998 cublasDgemv(cublas_handle, CUBLAS_OP_T, (int)w.size(), (int)Vflat.size() / (int)w.size(), &unit, devViBund, (int)w.size(), devw, 1, &zero, mulptr, 1);
999 cudaMemcpy(mul.data(), mulptr, (j+1) * sizeof(double), cudaMemcpyDeviceToHost);
1000
1001 for (int i = 0; i <= j; ++i)
1002 {
1003 double scal = 0.0;
1004 scal = mul[i];
1005 H[i][j] = scal;
1006 scal *= -1;
1007
1008 cublasDaxpy(cublas_handle, \
1009 (int)w.size(), \
1010 &scal, \
1011 devViBund + w.size() * i, 1, \
1012 devw, 1);
1013 }
1014 tRotB.stop();
1015
1016 tRotC.start();
1017
1018 cudaMemcpy(w.data(), devw, w.size() * sizeof(double), cudaMemcpyDeviceToHost);
1019
1020 m = j + 1;
1021
1022 double nrmw = 0.0;
1023
1024 //cublasDdot(cublas_handle, \
1025 // w.size(), \
1026 // devw, 1, /* w, 1 */ \
1027 // devw, 1, /* w, 1 */ \
1028 // &nrmw);
1029 //nrmw = sqrt(nrmw);
1030 //H[m][j] = nrmw;
1031
1032 nrmw = norm(w);
1033 H[m][j] = nrmw;
1034
1035 Vflat.insert(Vflat.end(), w.begin(), w.end());
1036 for (int t = 0; t < (totalVsize + nAfl); ++t)
1037 Vflat[Vflat.size() - t - 1] /= nrmw;
1038
1039 tRotC.stop();
1040
1041 tRotD.start();
1042
1043 if (IterRot(H, nrmRhs, gs, c, s, (int)m, (int)totalVsize, W.getPassport().numericalSchemes.gmresEps, (int)m, false))
1044 {
1045 std::cout << "iterations in GMRES = " << j + 1 << std::endl;
1046 if (j != 0)
1047 tIterOther.stop();
1048 tRotD.stop();
1049 tRot.stop();
1050
1051 break;
1052 }
1053
1054 tRotD.stop();
1055 tRot.stop();
1056
1057 if (j == 0)
1058 tIter1.stop();
1059 } // end of iterations
1060
1061 tIter.stop();
1062
1063 tPost.start();
1064
1065 W.getNonConstCuda().ReleaseSolution(W.getNonConstCuda().dev_sol);
1066 if (linScheme)
1067 W.getNonConstCuda().ReleaseSolution(W.getNonConstCuda().dev_solLin);
1068
1069 niter = (int)m;
1070
1071 std::vector<double> g(m + 1);
1072 g[0] = beta;
1073
1074 //GivensRotations
1075 double oldValue;
1076 for (int i = 0; i < m; i++)
1077 {
1078 oldValue = g[i];
1079 g[i] = c[i] * oldValue;
1080 g[i + 1] = -s[i] * oldValue;
1081 }
1082 //end of GivensRotations
1083
1084 // Solve HY=g
1085 std::vector<double> Y(m);
1086 Y[m - 1] = g[m - 1] / H[m - 1][m - 1];
1087
1088 double sum;
1089 for (int k = (int)m - 2; k >= 0; --k)
1090 {
1091 sum = 0.0;
1092 for (int s = k + 1; s < m; ++s)
1093 sum += H[k][s] * Y[s];
1094 Y[k] = (g[k] - sum) / H[k][k];
1095 }
1096 // end of Solve HY=g
1097
1098 size_t cntr = 0;
1099 for (size_t p = 0; p < nAfl; p++) {
1100
1101 size_t nPanelsP = W.getAirfoil(p).getNumberOfPanels();
1102
1103 if (!linScheme)
1104 for (size_t i = 0; i < nPanelsP; i++)
1105 {
1106 sum = 0.0;
1107 for (size_t j = 0; j < m; j++)
1108 sum += Vflat[(totalVsize + nAfl) * j + (i + cntr)] * Y[j];
1109 X[p][i] += sum;
1110 }
1111 else {
1112 for (size_t i = 0; i < 2 * nPanelsP; i++)
1113 {
1114 sum = 0.0;
1115 for (size_t j = 0; j < m; j++)
1116 sum += Vflat[(totalVsize + nAfl) * j + (i + cntr)] * Y[j];
1117 X[p][i] += sum;
1118 }
1119 }
1120 cntr += nPanelsP;
1121
1122 if (linScheme)
1123 cntr += nPanelsP;
1124
1125 }
1126 sum = 0.0;
1127 for (size_t p = 0; p < nAfl; p++)
1128 {
1129 for (size_t j = 0; j < m; j++)
1130 sum += Vflat[(totalVsize + nAfl) * j + (totalVsize + p)] * Y[j];
1131 //sum += V[j][totalVsize + p] * Y[j];
1132 R[p] += sum;
1133 }
1134
1135 tPost.stop();
1136
1137 tAll.stop();
1138
1139 treePnlInfl.MemoryFreeForGMRES();
1140 cudaFree(devViBund);
1141 cudaFree(devw);
1142 cudaFree(mulptr);
1143
1144 cublasDestroy(cublas_handle);
1145
1146 /*
1147 std::cout << "tAll = " << tAll.duration() << std::endl;
1148 std::cout << "tPre = " << tPre.duration() << std::endl;
1149 std::cout << "tIter = " << tIter.duration() << std::endl;
1150 std::cout << "tIter1 = " << tIter1.duration() << std::endl;
1151 if(m > 1)
1152 std::cout << "tIterN = " << tIterOther.duration() / (m-1.0) << std::endl;
1153 else
1154 std::cout << "tIterN = " << 0.0 << std::endl;
1155 std::cout << "tPost = " << tPost.duration() << std::endl;
1156
1157 std::cout << "tCG = " << tCG.duration() / (double)(m) << std::endl;
1158 std::cout << "tGC = " << tGC.duration() / (double)(m) << std::endl;
1159 std::cout << "tWrap(>0)= " << tWrapper.duration() / (double)(m - 1) << std::endl;
1160 std::cout << "tSplit = " << tSplit.duration() / (double)(m) << std::endl;
1161 std::cout << "tPrecond = " << tPrecond.duration() / (double)(m) << std::endl;
1162 std::cout << "tRot = " << (tRotA.duration() + tRotB.duration() + tRotC.duration() +tRotD.duration()) / (double)(m) << std::endl;
1163 //*/
1164
1165 //std::cout << " tRotA = " << tRotA.duration() / (double)(m) << std::endl;
1166 //std::cout << " tRotB = " << tRotB.duration() / (double)(m) << std::endl;
1167 //std::cout << " tRotC = " << tRotC.duration() / (double)(m) << std::endl;
1168 //std::cout << " tRotD = " << tRotD.duration() / (double)(m) << std::endl;
1169}
1170
1171#endif
Airfoil2D.h
Заголовочный файл с описанием класса Airfoil.
Boundary2D.h
Заголовочный файл с описанием класса Boundary.
PreCalculateCoef
void PreCalculateCoef(const Airfoil &aflP, size_t nPanelsP, std::vector< Point2D > &prea, std::vector< Point2D > &prec, std::vector< Point2D > &prea1, std::vector< Point2D > &prec1, bool linScheme)
Definition Gmres2D.cpp:528
MapVec
Eigen::Map< const Eigen::VectorXd > MapVec
Definition Gmres2D.cpp:57
Gmres2D.h
Заголовочный файл с функциями для метода GMRES.
MeasureVP2D.h
Заголовочный файл с описанием класса MeasureVP.
Mechanics2D.h
Заголовочный файл с описанием класса Mechanics.
Preprocessor.h
Заголовочный файл с описанием класса Preprocessor.
StreamParser.h
Заголовочный файл с описанием класса StreamParser.
Velocity2D.h
Заголовочный файл с описанием класса Velocity.
Wake2D.h
Заголовочный файл с описанием класса Wake.
World2D.h
Заголовочный файл с описанием класса World2D.
VM2D::AirfoilGeometry::len
std::vector< double > len
Длины панелей профиля
Definition Airfoil2D.h:94
VM2D::AirfoilGeometry::getR
const Point2D & getR(size_t q) const
Возврат константной ссылки на вершину профиля
Definition Airfoil2D.h:113
VM2D::AirfoilGeometry::tau
std::vector< Point2D > tau
Касательные к панелям профиля
Definition Airfoil2D.h:91
VM2D::AirfoilGeometry::getNumberOfPanels
size_t getNumberOfPanels() const
Возврат количества панелей на профиле
Definition Airfoil2D.h:163
VM2D::Airfoil
Абстрактный класс, определяющий обтекаемый профиль
Definition Airfoil2D.h:182
VM2D::Airfoil::isAfter
bool isAfter(size_t i, size_t j) const
Проверка, идет ли вершина i следом за вершиной j.
Definition Airfoil2D.cpp:73
VM2D::Passport::numericalSchemes
NumericalSchemes numericalSchemes
Структура с используемыми численными схемами
Definition Passport2D.h:295
VM2D::World2D
Класс, опеделяющий текущую решаемую задачу
Definition World2D.h:74
VM2D::World2D::getNumberOfAirfoil
size_t getNumberOfAirfoil() const
Возврат количества профилей в задаче
Definition World2D.h:174
VM2D::World2D::getAirfoil
const Airfoil & getAirfoil(size_t i) const
Возврат константной ссылки на объект профиля
Definition World2D.h:157
VM2D::World2D::getNonConstCuda
Gpu & getNonConstCuda() const
Возврат неконстантной ссылки на объект, связанный с видеокартой (GPU)
Definition World2D.h:266
VM2D::World2D::getCuda
const Gpu & getCuda() const
Возврат константной ссылки на объект, связанный с видеокартой (GPU)
Definition World2D.h:261
VM2D::World2D::getPassport
const Passport & getPassport() const
Возврат константной ссылки на паспорт
Definition World2D.h:251
VMlib::numvector
Шаблонный класс, определяющий вектор фиксированной длины Фактически представляет собой массив,...
Definition numvector.h:99
VMlib::numvector::kcross
numvector< T, 2 > kcross() const
Геометрический поворот двумерного вектора на 90 градусов
Definition numvector.h:510
VMlib::numvector::length2
auto length2() const -> typename std::remove_const< typename std::remove_reference< decltype(this->data[0])>::type >::type
Вычисление квадрата нормы (длины) вектора
Definition numvector.h:386
VMlib::numvector::unit
auto unit(P newlen=1) const -> numvector< typename std::remove_const< decltype(this->data[0] *newlen)>::type, n >
Вычисление орта вектора или вектора заданной длины, коллинеарного данному
Definition numvector.h:402
VMlib::numvector::length
P length() const
Вычисление 2-нормы (длины) вектора
Definition numvector.h:374
VMlib::vmTimer
Класс засекания времени
Definition TimesGen.h:59
IDPI
const double IDPI
Число .
Definition defs.h:85
VM2D::GMRES_Direct
void GMRES_Direct(const World2D &W, int nAllVars, int nafl, const std::vector< double > &mtrDir, const std::vector< double > &rhsDir, const std::vector< int > &pos, const std::vector< int > &vsize, std::vector< std::vector< double > > &gam, std::vector< double > &R)
Definition Gmres2D.cpp:206
VM2D::norm
double norm(const std::vector< T > &b)
Шаблонная функция вычисления евклидовой нормы вектора или списка
Definition Gmres2D.h:93
VM2D::SolCircleRun
void SolCircleRun(std::vector< double > &AX, const std::vector< double > &rhs, const Airfoil &afl, const std::vector< std::vector< Point2D > > &prea1, const std::vector< std::vector< Point2D > > &prec1, int p, int n, sweepVectors &sw)
Definition Gmres2D.cpp:393
VM2D::IterRot
bool IterRot(std::vector< std::vector< double > > &H, const double nrmRhs, double &gs, std::vector< double > &c, std::vector< double > &s, int m, int n, double epsGMRES, int iter, bool residualShow)
Контроль невязки после выполнения очередной итерации GMRES.
Definition Gmres2D.cpp:59
VM2D::SolMdirect
void SolMdirect(const std::vector< double > &A, const std::vector< double > &rhs, size_t startRow, size_t startRowReg, size_t np, std::vector< double > &res, bool lin)
Definition Gmres2D.cpp:130
VM2D::SolCircleRundirect
void SolCircleRundirect(const std::vector< double > &A, const std::vector< double > &rhs, size_t startRow, size_t startRowReg, size_t np, std::vector< double > &res)
Definition Gmres2D.cpp:93
VM2D::SolM
void SolM(std::vector< double > &AX, const std::vector< double > &rhs, const Airfoil &afl, const std::vector< std::vector< Point2D > > &prea, const std::vector< std::vector< Point2D > > &prec, int p, int n, bool linScheme, sweepVectors &sw)
Definition Gmres2D.cpp:444
VM2D::norm2
double norm2(const std::vector< T > &b)
Шаблонная функция вычисления евклидовой нормы вектора или списка
Definition Gmres2D.h:108
VMlib::Lambda
double Lambda(const Point2D &p, const Point2D &s)
Вспомогательная функция вычисления логарифма отношения норм векторов
Definition defs.cpp:268
VMlib::Alpha
double Alpha(const Point2D &p, const Point2D &s)
Вспомогательная функция вычисления угла между векторами
Definition defs.cpp:262
VMlib::Omega
Point2D Omega(const Point2D &a, const Point2D &b, const Point2D &c)
Вспомогательная функция вычисления величины .
Definition defs.cpp:274
VM2D::sweepVectors
Структура, определяющий необходимые массивы для рализации метода прогонки
Definition Gmres2D.h:65
VM2D::sweepVectors::phi
std::vector< double > phi
Definition Gmres2D.h:70
VM2D::sweepVectors::beta
std::vector< double > beta
Definition Gmres2D.h:67
VM2D::sweepVectors::psi
std::vector< double > psi
Definition Gmres2D.h:72
VM2D::sweepVectors::delta
std::vector< double > delta
Definition Gmres2D.h:69
VM2D::sweepVectors::alpha
std::vector< double > alpha
Definition Gmres2D.h:66
VM2D::sweepVectors::gamma
std::vector< double > gamma
Definition Gmres2D.h:68
VM2D::sweepVectors::xi
std::vector< double > xi
Definition Gmres2D.h:71