Loading...
Searching...
No Matches
ReducedUnsteadyNSTurb.C
Go to the documentation of this file.
1/*---------------------------------------------------------------------------*\
2 ██╗████████╗██╗ ██╗ █████╗ ██████╗ █████╗ ███████╗██╗ ██╗
3 ██║╚══██╔══╝██║ ██║██╔══██╗██╔════╝██╔══██╗ ██╔════╝██║ ██║
4 ██║ ██║ ███████║███████║██║ ███████║█████╗█████╗ ██║ ██║
5 ██║ ██║ ██╔══██║██╔══██║██║ ██╔══██║╚════╝██╔══╝ ╚██╗ ██╔╝
6 ██║ ██║ ██║ ██║██║ ██║╚██████╗██║ ██║ ██║ ╚████╔╝
7 ╚═╝ ╚═╝ ╚═╝ ╚═╝╚═╝ ╚═╝ ╚═════╝╚═╝ ╚═╝ ╚═╝ ╚═══╝
8
9 * In real Time Highly Advanced Computational Applications for Finite Volumes
10 * Copyright (C) 2017 by the ITHACA-FV authors
11-------------------------------------------------------------------------------
12
13License
14 This file is part of ITHACA-FV
15
16 ITHACA-FV is free software: you can redistribute it and/or modify
17 it under the terms of the GNU Lesser General Public License as published by
18 the Free Software Foundation, either version 3 of the License, or
19 (at your option) any later version.
20
21 ITHACA-FV is distributed in the hope that it will be useful,
22 but WITHOUT ANY WARRANTY; without even the implied warranty of
23 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
24 GNU Lesser General Public License for more details.
25
26 You should have received a copy of the GNU Lesser General Public License
27 along with ITHACA-FV. If not, see <http://www.gnu.org/licenses/>.
28
29\*---------------------------------------------------------------------------*/
30
33
34
35#include "ReducedUnsteadyNSTurb.H"
36
37
38// * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * * * //
39
40// Constructor
44
45ReducedUnsteadyNSTurb::ReducedUnsteadyNSTurb(UnsteadyNSTurb& fomProblem)
46{
47 problem = &fomProblem;
48 N_BC = problem->inletIndex.rows();
49 Nphi_u = problem->B_matrix.rows();
50 Nphi_p = problem->K_matrix.cols();
51 nphiNut = problem->cTotalTensor.dimension(1);
52 dimA = problem->velRBF.cols();
53 interChoice = problem->interChoice;
54
55 // Create locally the velocity modes
56 for (int k = 0; k < problem->liftfield.size(); k++)
57 {
58 Umodes.append((problem->liftfield[k]).clone());
59 }
60
61 for (int k = 0; k < problem->NUmodes; k++)
62 {
63 Umodes.append((problem->Umodes[k]).clone());
64 }
65
66 for (int k = 0; k < problem->NSUPmodes; k++)
67 {
68 Umodes.append((problem->supmodes[k]).clone());
69 }
70
71 // Create locally the pressure modes
72 for (int k = 0; k < problem->NPmodes; k++)
73 {
74 Pmodes.append((problem->Pmodes[k]).clone());
75 }
76
77 // Create locally the eddy viscosity modes
78 for (int k = 0; k < problem->nNutModes; k++)
79 {
80 nutModes.append((problem->nutModes[k]).clone());
81 }
82
83 newtonObjectSUP = newtonUnsteadyNSTurbSUP(Nphi_u + Nphi_p, Nphi_u + Nphi_p,
84 fomProblem);
85 newtonObjectPPE = newtonUnsteadyNSTurbPPE(Nphi_u + Nphi_p, Nphi_u + Nphi_p,
86 fomProblem);
87 newtonObjectSUPAve = newtonUnsteadyNSTurbSUPAve(Nphi_u + Nphi_p,
88 Nphi_u + Nphi_p,
89 fomProblem);
90 newtonObjectPPEAve = newtonUnsteadyNSTurbPPEAve(Nphi_u + Nphi_p,
91 Nphi_u + Nphi_p,
92 fomProblem);
93}
94
95
96// * * * * * * * * * * * * * * * Operators supremizer * * * * * * * * * * * * * //
97
98// Operator to evaluate the residual for the supremizer approach
99int newtonUnsteadyNSTurbSUP::operator()(const Eigen::VectorXd& x,
100 Eigen::VectorXd& fvec) const
101{
102 Eigen::VectorXd a_dot(Nphi_u);
103 Eigen::VectorXd aTmp(Nphi_u);
104 Eigen::VectorXd bTmp(Nphi_p);
105 aTmp = x.head(Nphi_u);
106 bTmp = x.tail(Nphi_p);
107
108 // Choose the order of the numerical difference scheme for approximating the time derivative
109 if (problem->timeDerivativeSchemeOrder == "first")
110 {
111 a_dot = (x.head(Nphi_u) - y_old.head(Nphi_u)) / dt;
112 }
113 else
114 {
115 a_dot = (1.5 * x.head(Nphi_u) - 2 * y_old.head(Nphi_u) + 0.5 * yOldOld.head(
116 Nphi_u)) / dt;
117 }
118
119 // Convective term
120 Eigen::MatrixXd cc(1, 1);
121 // Mom Term
122 Eigen::VectorXd m1 = problem->bTotalMatrix * aTmp * nu;
123 // Gradient of pressure
124 Eigen::VectorXd m2 = problem->K_matrix * bTmp;
125 // Mass Term
126 Eigen::VectorXd m5 = problem->M_matrix * a_dot;
127 // Pressure Term
128 Eigen::VectorXd m3 = problem->P_matrix * aTmp;
129 // Penalty term
130 Eigen::MatrixXd penaltyU = Eigen::MatrixXd::Zero(Nphi_u, N_BC);
131
132 // Term for penalty method
133 if (problem->bcMethod == "penalty")
134 {
135 for (int l = 0; l < N_BC; l++)
136 {
137 penaltyU.col(l) = bc(l) * problem->bcVelVec[l] - problem->bcVelMat[l] *
138 aTmp;
139 }
140 }
141
142 for (int i = 0; i < Nphi_u; i++)
143 {
144 cc = aTmp.transpose() * Eigen::SliceFromTensor(problem->C_tensor, 0,
145 i) * aTmp - gNut.transpose() *
146 Eigen::SliceFromTensor(problem->cTotalTensor, 0, i) * aTmp;
147 fvec(i) = - m5(i) + m1(i) - cc(0, 0) - m2(i);
148
149 if (problem->bcMethod == "penalty")
150 {
151 fvec(i) += ((penaltyU * tauU)(i, 0));
152 }
153 }
154
155 for (int j = 0; j < Nphi_p; j++)
156 {
157 int k = j + Nphi_u;
158 fvec(k) = m3(j);
159 }
160
161 if (problem->bcMethod == "lift")
162 {
163 for (int j = 0; j < N_BC; j++)
164 {
165 fvec(j) = x(j) - bc(j);
166 }
167 }
168
169 return 0;
170}
171
172// Operator to evaluate the Jacobian for the supremizer approach
173int newtonUnsteadyNSTurbSUP::df(const Eigen::VectorXd& x,
174 Eigen::MatrixXd& fjac) const
175{
176 Eigen::NumericalDiff<newtonUnsteadyNSTurbSUP> numDiff(*this);
177 numDiff.df(x, fjac);
178 return 0;
179}
180
181// Operator to evaluate the residual for the supremizer approach
182int newtonUnsteadyNSTurbSUPAve::operator()(const Eigen::VectorXd& x,
183 Eigen::VectorXd& fvec) const
184{
185 Eigen::VectorXd a_dot(Nphi_u);
186 Eigen::VectorXd aTmp(Nphi_u);
187 Eigen::VectorXd bTmp(Nphi_p);
188 aTmp = x.head(Nphi_u);
189 bTmp = x.tail(Nphi_p);
190
191 // Choose the order of the numerical difference scheme for approximating the time derivative
192 if (problem->timeDerivativeSchemeOrder == "first")
193 {
194 a_dot = (x.head(Nphi_u) - y_old.head(Nphi_u)) / dt;
195 }
196 else
197 {
198 a_dot = (1.5 * x.head(Nphi_u) - 2 * y_old.head(Nphi_u) + 0.5 * yOldOld.head(
199 Nphi_u)) / dt;
200 }
201
202 // Convective term
203 Eigen::MatrixXd cc(1, 1);
204 // Mom Term
205 Eigen::VectorXd m1 = problem->bTotalMatrix * aTmp * nu;
206 // Gradient of pressure
207 Eigen::VectorXd m2 = problem->K_matrix * bTmp;
208 // Mass Term
209 Eigen::VectorXd m5 = problem->M_matrix * a_dot;
210 // Pressure Term
211 Eigen::VectorXd m3 = problem->P_matrix * aTmp;
212 // Penalty term
213 Eigen::MatrixXd penaltyU = Eigen::MatrixXd::Zero(Nphi_u, N_BC);
214
215 // Term for penalty method
216 if (problem->bcMethod == "penalty")
217 {
218 for (int l = 0; l < N_BC; l++)
219 {
220 penaltyU.col(l) = bc(l) * problem->bcVelVec[l] - problem->bcVelMat[l] *
221 aTmp;
222 }
223 }
224
225 for (int i = 0; i < Nphi_u; i++)
226 {
227 cc = aTmp.transpose() * Eigen::SliceFromTensor(problem->C_tensor, 0,
228 i) * aTmp - gNut.transpose() *
229 Eigen::SliceFromTensor(problem->cTotalTensor, 0,
230 i) * aTmp - gNutAve.transpose() *
231 Eigen::SliceFromTensor(problem->cTotalAveTensor, 0, i) * aTmp;
232 fvec(i) = - m5(i) + m1(i) - cc(0, 0) - m2(i);
233
234 if (problem->bcMethod == "penalty")
235 {
236 fvec(i) += ((penaltyU * tauU)(i, 0));
237 }
238 }
239
240 for (int j = 0; j < Nphi_p; j++)
241 {
242 int k = j + Nphi_u;
243 fvec(k) = m3(j);
244 }
245
246 if (problem->bcMethod == "lift")
247 {
248 for (int j = 0; j < N_BC; j++)
249 {
250 fvec(j) = x(j) - bc(j);
251 }
252 }
253
254 return 0;
255}
256
257// Operator to evaluate the Jacobian for the supremizer approach
258int newtonUnsteadyNSTurbSUPAve::df(const Eigen::VectorXd& x,
259 Eigen::MatrixXd& fjac) const
260{
261 Eigen::NumericalDiff<newtonUnsteadyNSTurbSUPAve> numDiff(*this);
262 numDiff.df(x, fjac);
263 return 0;
264}
265
266// * * * * * * * * * * * * * * * Operators PPE * * * * * * * * * * * * * * * //
267
268// Operator to evaluate the residual for the supremizer approach
269int newtonUnsteadyNSTurbPPE::operator()(const Eigen::VectorXd& x,
270 Eigen::VectorXd& fvec) const
271{
272 Eigen::VectorXd a_dot(Nphi_u);
273 Eigen::VectorXd aTmp(Nphi_u);
274 Eigen::VectorXd bTmp(Nphi_p);
275 aTmp = x.head(Nphi_u);
276 bTmp = x.tail(Nphi_p);
277
278 // Choose the order of the numerical difference scheme for approximating the time derivative
279 if (problem->timeDerivativeSchemeOrder == "first")
280 {
281 a_dot = (x.head(Nphi_u) - y_old.head(Nphi_u)) / dt;
282 }
283 else
284 {
285 a_dot = (1.5 * x.head(Nphi_u) - 2 * y_old.head(Nphi_u) + 0.5 * yOldOld.head(
286 Nphi_u)) / dt;
287 }
288
289 // Convective terms
290 Eigen::MatrixXd cc(1, 1);
291 Eigen::MatrixXd gg(1, 1);
292 Eigen::MatrixXd bb(1, 1);
293 // Mom Term
294 Eigen::VectorXd m1 = problem->bTotalMatrix * aTmp * nu;
295 // Gradient of pressure
296 Eigen::VectorXd m2 = problem->K_matrix * bTmp;
297 // Mass Term
298 Eigen::VectorXd m5 = problem->M_matrix * a_dot;
299 // Pressure Term
300 Eigen::VectorXd m3 = problem->D_matrix * bTmp;
301 // BC PPE
302 Eigen::VectorXd m6 = problem->BC1_matrix * aTmp * nu;
303 // BC PPE
304 Eigen::VectorXd m7 = problem->BC3_matrix * aTmp * nu;
305 // Penalty term
306 Eigen::MatrixXd penaltyU = Eigen::MatrixXd::Zero(Nphi_u, N_BC);
307
308 // Term for penalty method
309 if (problem->bcMethod == "penalty")
310 {
311 for (int l = 0; l < N_BC; l++)
312 {
313 penaltyU.col(l) = bc(l) * problem->bcVelVec[l] - problem->bcVelMat[l] *
314 aTmp;
315 }
316 }
317
318 for (int i = 0; i < Nphi_u; i++)
319 {
320 cc = aTmp.transpose() * Eigen::SliceFromTensor(problem->C_tensor, 0,
321 i) * aTmp - gNut.transpose() *
322 Eigen::SliceFromTensor(problem->cTotalTensor, 0, i) * aTmp;
323 fvec(i) = - m5(i) + m1(i) - cc(0, 0) - m2(i);
324
325 if (problem->bcMethod == "penalty")
326 {
327 fvec(i) += ((penaltyU * tauU)(i, 0));
328 }
329 }
330
331 for (int j = 0; j < Nphi_p; j++)
332 {
333 int k = j + Nphi_u;
334 gg = aTmp.transpose() * Eigen::SliceFromTensor(problem->gTensor, 0,
335 j) * aTmp;
336 bb = aTmp.transpose() * Eigen::SliceFromTensor(problem->bc2Tensor, 0,
337 j) * aTmp;
338 //fvec(k) = m3(j, 0) - gg(0, 0) - m6(j, 0) + bb(0, 0);
339 fvec(k) = m3(j, 0) + gg(0, 0) - m7(j, 0);
340 }
341
342 if (problem->bcMethod == "lift")
343 {
344 for (int j = 0; j < N_BC; j++)
345 {
346 fvec(j) = x(j) - bc(j);
347 }
348 }
349
350 return 0;
351}
352
353// Operator to evaluate the Jacobian for the supremizer approach
354int newtonUnsteadyNSTurbPPE::df(const Eigen::VectorXd& x,
355 Eigen::MatrixXd& fjac) const
356{
357 Eigen::NumericalDiff<newtonUnsteadyNSTurbPPE> numDiff(*this);
358 numDiff.df(x, fjac);
359 return 0;
360}
361
362int newtonUnsteadyNSTurbPPEAve::operator()(const Eigen::VectorXd& x,
363 Eigen::VectorXd& fvec) const
364{
365 Eigen::VectorXd a_dot(Nphi_u);
366 Eigen::VectorXd aTmp(Nphi_u);
367 Eigen::VectorXd bTmp(Nphi_p);
368 aTmp = x.head(Nphi_u);
369 bTmp = x.tail(Nphi_p);
370
371 // Choose the order of the numerical difference scheme for approximating the time derivative
372 if (problem->timeDerivativeSchemeOrder == "first")
373 {
374 a_dot = (x.head(Nphi_u) - y_old.head(Nphi_u)) / dt;
375 }
376 else
377 {
378 a_dot = (1.5 * x.head(Nphi_u) - 2 * y_old.head(Nphi_u) + 0.5 * yOldOld.head(
379 Nphi_u)) / dt;
380 }
381
382 // Convective terms
383 Eigen::MatrixXd cc(1, 1);
384 Eigen::MatrixXd gg(1, 1);
385 Eigen::MatrixXd bb(1, 1);
386 Eigen::MatrixXd nn(1, 1);
387 // Mom Term
388 Eigen::VectorXd m1 = problem->bTotalMatrix * aTmp * nu;
389 // Gradient of pressure
390 Eigen::VectorXd m2 = problem->K_matrix * bTmp;
391 // Mass Term
392 Eigen::VectorXd m5 = problem->M_matrix * a_dot;
393 // Time-derivative of the divergence Term
394 //Eigen::VectorXd m5 = problem->M_matrix * a_dot;
395 // Pressure Term
396 Eigen::VectorXd m3 = problem->D_matrix * bTmp;
397 // BC PPE
398 Eigen::VectorXd m6 = problem->BC1_matrix * aTmp * nu;
399 // BC PPE
400 Eigen::VectorXd m7 = problem->BC3_matrix * aTmp * nu;
401 // Penalty term
402 Eigen::MatrixXd penaltyU = Eigen::MatrixXd::Zero(Nphi_u, N_BC);
403
404 // Term for penalty method
405 if (problem->bcMethod == "penalty")
406 {
407 for (int l = 0; l < N_BC; l++)
408 {
409 penaltyU.col(l) = bc(l) * problem->bcVelVec[l] - problem->bcVelMat[l] *
410 aTmp;
411 }
412 }
413
414 for (int i = 0; i < Nphi_u; i++)
415 {
416 cc = aTmp.transpose() * Eigen::SliceFromTensor(problem->C_tensor, 0,
417 i) * aTmp - gNut.transpose() *
418 Eigen::SliceFromTensor(problem->cTotalTensor, 0,
419 i) * aTmp - gNutAve.transpose() *
420 Eigen::SliceFromTensor(problem->cTotalAveTensor, 0, i) * aTmp;
421 fvec(i) = - m5(i) + m1(i) - cc(0, 0) - m2(i);
422
423 if (problem->bcMethod == "penalty")
424 {
425 fvec(i) += ((penaltyU * tauU)(i, 0));
426 }
427 }
428
429 for (int j = 0; j < Nphi_p; j++)
430 {
431 int k = j + Nphi_u;
432 gg = aTmp.transpose() * Eigen::SliceFromTensor(problem->gTensor, 0,
433 j) * aTmp;
434 bb = aTmp.transpose() * Eigen::SliceFromTensor(problem->bc2Tensor, 0,
435 j) * aTmp;
436 nn = gNut.transpose() *
437 Eigen::SliceFromTensor(problem->cTotalPPETensor, 0,
438 j) * aTmp + gNutAve.transpose() *
439 Eigen::SliceFromTensor(problem->cTotalPPEAveTensor, 0, j) * aTmp;
440 //fvec(k) = m3(j, 0) + gg(0, 0) - m6(j, 0) + bb(0, 0) - nn(0, 0);
441 fvec(k) = m3(j, 0) + gg(0, 0) - m7(j, 0) - nn(0, 0);
442 //fvec(k) = m3(j, 0) + gg(0, 0) - m7(j, 0);
443 }
444
445 if (problem->bcMethod == "lift")
446 {
447 for (int j = 0; j < N_BC; j++)
448 {
449 fvec(j) = x(j) - bc(j);
450 }
451 }
452
453 return 0;
454}
455
456// Operator to evaluate the Jacobian for the PPE approach
457int newtonUnsteadyNSTurbPPEAve::df(const Eigen::VectorXd& x,
458 Eigen::MatrixXd& fjac) const
459{
460 Eigen::NumericalDiff<newtonUnsteadyNSTurbPPEAve> numDiff(*this);
461 numDiff.df(x, fjac);
462 return 0;
463}
464
465
466// * * * * * * * * * * * * * * * Solve Functions * * * * * * * * * * * * * //
467void ReducedUnsteadyNSTurb::solveOnlineSUP(Eigen::MatrixXd vel)
468{
469 M_Assert(exportEvery >= dt,
470 "The time step dt must be smaller than exportEvery.");
471 M_Assert(storeEvery >= dt,
472 "The time step dt must be smaller than storeEvery.");
473 M_Assert(ITHACAutilities::isInteger(storeEvery / dt) == true,
474 "The variable storeEvery must be an integer multiple of the time step dt.");
475 M_Assert(ITHACAutilities::isInteger(exportEvery / dt) == true,
476 "The variable exportEvery must be an integer multiple of the time step dt.");
477 M_Assert(ITHACAutilities::isInteger(exportEvery / storeEvery) == true,
478 "The variable exportEvery must be an integer multiple of the variable storeEvery.");
479 int numberOfStores = round(storeEvery / dt);
480
481 if (problem->bcMethod == "lift")
482 {
483 vel_now = setOnlineVelocity(vel);
484 }
485 else if (problem->bcMethod == "penalty")
486 {
487 vel_now = vel;
488 }
489
490 // Create and resize the solution vector
491 y.resize(Nphi_u + Nphi_p, 1);
492 y.setZero();
493 // Set Initial Conditions
494 // y.head(Nphi_u) = initCond.col(0).head(Nphi_u);
495 // y.tail(Nphi_p) = initCond.col(0).tail(Nphi_p);
496 y.head(Nphi_u) = ITHACAutilities::getCoeffs(problem->Ufield[0],
497 Umodes);
498 y.tail(Nphi_p) = ITHACAutilities::getCoeffs(problem->Pfield[0],
499 Pmodes);
500 nut0 = ITHACAutilities::getCoeffs(problem->nutFields[0],
501 nutModes);
502 int nextStore = 0;
503 int counter2 = 0;
504
505 // Change initial condition for the lifting function
506 if (problem->bcMethod == "lift")
507 {
508 for (int j = 0; j < N_BC; j++)
509 {
510 y(j) = vel_now(j, 0);
511 }
512 }
513
514 int firstRBFInd;
515
516 if (skipLift == true && problem->bcMethod == "lift")
517 {
518 firstRBFInd = N_BC;
519 }
520 else
521 {
522 firstRBFInd = 0;
523 }
524
525 // Set some properties of the newton object
526 newtonObjectSUP.nu = nu;
527 newtonObjectSUP.y_old = y;
528 newtonObjectSUP.yOldOld = newtonObjectSUP.y_old;
529 newtonObjectSUP.dt = dt;
530 newtonObjectSUP.bc.resize(N_BC);
531 newtonObjectSUP.tauU = tauU;
532 newtonObjectSUP.gNut = nut0;
533
534 for (int j = 0; j < N_BC; j++)
535 {
536 newtonObjectSUP.bc(j) = vel_now(j, 0);
537 }
538
539 // Set number of online solutions
540 int Ntsteps = static_cast<int>((finalTime - tstart) / dt);
541 int onlineSize = static_cast<int>(Ntsteps / numberOfStores);
542 online_solution.resize(onlineSize);
543 rbfCoeffMat.resize(nphiNut + 1, onlineSize + 3);
544 // Set the initial time
545 time = tstart;
546 // Counting variable
547 int counter = 0;
548 // Create vector to store temporal solution and save initial condition as first solution
549 Eigen::MatrixXd tmp_sol(Nphi_u + Nphi_p + 1, 1);
550 tmp_sol(0) = time;
551 tmp_sol.col(0).tail(y.rows()) = y;
552
553 // This part up to the while loop is just to compute the eddy viscosity field at time = startTime if time is not equal to zero
554 if ((time != 0) || (startFromZero == true))
555 {
556 online_solution[counter] = tmp_sol;
557 counter ++;
558 rbfCoeffMat(0, counter2) = time;
559 rbfCoeffMat.block(1, counter2, nphiNut, 1) = nut0;
560 counter2++;
561 nextStore += numberOfStores;
562 }
563
564 // Create nonlinear solver object
565 Eigen::HybridNonLinearSolver<newtonUnsteadyNSTurbSUP> hnls(newtonObjectSUP);
566 // Set output colors for fancy output
567 Color::Modifier red(Color::FG_RED);
568 Color::Modifier green(Color::FG_GREEN);
569 Color::Modifier def(Color::FG_DEFAULT);
570
571 // Start the time loop
572 while (time < finalTime)
573 {
574 time = time + dt;
575 Eigen::VectorXd res(y);
576 res.setZero();
577 hnls.solve(y);
578 newtonObjectSUP.operator()(y, res);
579 Eigen::VectorXd tv;
580 Eigen::VectorXd aDer;
581 aDer = (y.head(Nphi_u) - newtonObjectSUP.y_old.head(Nphi_u)) / dt;
582 tv.resize(dimA);
583
584 switch (interChoice)
585 {
586 case 1:
587 tv << y.segment(firstRBFInd, dimA);
588 break;
589
590 case 2:
591 tv << muStar, y.segment(firstRBFInd, dimA - muStar.size());
592 break;
593
594 case 3:
595 tv << y.segment(firstRBFInd, dimA / 2), aDer.segment(firstRBFInd, dimA / 2);
596 break;
597
598 case 4:
599 tv << muStar, y.segment(firstRBFInd, (dimA - muStar.size()) / 2),
600 aDer.segment(firstRBFInd, (dimA - muStar.size()) / 2);
601 break;
602
603 default:
604 tv << y.segment(firstRBFInd, dimA);
605 break;
606 }
607
608 for (int i = 0; i < nphiNut; i++)
609 {
610 newtonObjectSUP.gNut(i) = problem->rbfSplines[i]->eval(tv);
611 }
612
613 // Change initial condition for the lifting function
614 if (problem->bcMethod == "lift")
615 {
616 for (int j = 0; j < N_BC; j++)
617 {
618 y(j) = vel_now(j, 0);
619 }
620 }
621
622 newtonObjectSUP.operator()(y, res);
623 newtonObjectSUP.yOldOld = newtonObjectSUP.y_old;
624 newtonObjectSUP.y_old = y;
625 std::cout << "################## Online solve N° " << count_online_solve <<
626 " ##################" << std::endl;
627 Info << "Time = " << time << endl;
628 std::cout << "Solving for the parameter: " << vel_now << std::endl;
629
630 if (res.norm() < 1e-5)
631 {
632 std::cout << green << "|F(x)| = " << res.norm() << " - Minimun reached in " <<
633 hnls.iter << " iterations " << def << std::endl << std::endl;
634 }
635 else
636 {
637 std::cout << red << "|F(x)| = " << res.norm() << " - Minimun reached in " <<
638 hnls.iter << " iterations " << def << std::endl << std::endl;
639 }
640
641 count_online_solve += 1;
642 tmp_sol(0) = time;
643 tmp_sol.col(0).tail(y.rows()) = y;
644
645 if (counter == nextStore)
646 {
647 if (counter2 >= online_solution.size())
648 {
649 online_solution.append(tmp_sol);
650 }
651 else
652 {
653 online_solution[counter2] = tmp_sol;
654 }
655
656 rbfCoeffMat(0, counter2) = time;
657 rbfCoeffMat.block(1, counter2, nphiNut, 1) = newtonObjectSUP.gNut;
658 nextStore += numberOfStores;
659 counter2 ++;
660 }
661
662 counter ++;
663 }
664
665 // Save the solution
666 ITHACAstream::exportMatrix(online_solution, "red_coeff", "python",
667 "./ITHACAoutput/red_coeff");
668 ITHACAstream::exportMatrix(online_solution, "red_coeff", "matlab",
669 "./ITHACAoutput/red_coeff");
670 count_online_solve += 1;
671}
672
673void ReducedUnsteadyNSTurb::solveOnlineSUPAve(Eigen::MatrixXd vel)
674{
675 M_Assert(exportEvery >= dt,
676 "The time step dt must be smaller than exportEvery.");
677 M_Assert(storeEvery >= dt,
678 "The time step dt must be smaller than storeEvery.");
679 M_Assert(ITHACAutilities::isInteger(storeEvery / dt) == true,
680 "The variable storeEvery must be an integer multiple of the time step dt.");
681 M_Assert(ITHACAutilities::isInteger(exportEvery / dt) == true,
682 "The variable exportEvery must be an integer multiple of the time step dt.");
683 M_Assert(ITHACAutilities::isInteger(exportEvery / storeEvery) == true,
684 "The variable exportEvery must be an integer multiple of the variable storeEvery.");
685 int numberOfStores = round(storeEvery / dt);
686
687 if (problem->bcMethod == "lift")
688 {
689 vel_now = setOnlineVelocity(vel);
690 }
691 else if (problem->bcMethod == "penalty")
692 {
693 vel_now = vel;
694 }
695
696 // Create and resize the solution vector
697 y.resize(Nphi_u + Nphi_p, 1);
698 y.setZero();
699 // Set Initial Conditions
700 y.head(Nphi_u) = initCond.col(0).head(Nphi_u);
701 y.tail(Nphi_p) = initCond.col(0).tail(Nphi_p);
702 int nextStore = 0;
703 int counter2 = 0;
704
705 // Change initial condition for the lifting function
706 if (problem->bcMethod == "lift")
707 {
708 for (int j = 0; j < N_BC; j++)
709 {
710 y(j) = vel_now(j, 0);
711 }
712 }
713
714 int firstRBFInd;
715
716 if (skipLift == true && problem->bcMethod == "lift")
717 {
718 firstRBFInd = N_BC;
719 }
720 else
721 {
722 firstRBFInd = 0;
723 }
724
725 // Set some properties of the newton object
726 newtonObjectSUPAve.nu = nu;
727 newtonObjectSUPAve.y_old = y;
728 newtonObjectSUPAve.yOldOld = newtonObjectSUPAve.y_old;
729 newtonObjectSUPAve.dt = dt;
730 newtonObjectSUPAve.bc.resize(N_BC);
731 newtonObjectSUPAve.tauU = tauU;
732 newtonObjectSUPAve.gNut = nut0;
733 newtonObjectSUPAve.gNutAve = gNutAve;
734
735 for (int j = 0; j < N_BC; j++)
736 {
737 newtonObjectSUPAve.bc(j) = vel_now(j, 0);
738 }
739
740 // Set number of online solutions
741 int Ntsteps = static_cast<int>((finalTime - tstart) / dt);
742 int onlineSize = static_cast<int>(Ntsteps / numberOfStores);
743 online_solution.resize(onlineSize);
744 rbfCoeffMat.resize(nphiNut + 1, onlineSize + 3);
745 // Set the initial time
746 time = tstart;
747 // Counting variable
748 int counter = 0;
749 // Create vector to store temporal solution and save initial condition as first solution
750 Eigen::MatrixXd tmp_sol(Nphi_u + Nphi_p + 1, 1);
751 tmp_sol(0) = time;
752 tmp_sol.col(0).tail(y.rows()) = y;
753
754 // This part up to the while loop is just to compute the eddy viscosity field at time = startTime if time is not equal to zero
755 if ((time != 0) || (startFromZero == true))
756 {
757 online_solution[counter] = tmp_sol;
758 counter ++;
759 rbfCoeffMat(0, counter2) = time;
760 rbfCoeffMat.block(1, counter2, nphiNut, 1) = nut0;
761 counter2++;
762 nextStore += numberOfStores;
763 }
764
765 // Create nonlinear solver object
766 Eigen::HybridNonLinearSolver<newtonUnsteadyNSTurbSUPAve> hnls(
767 newtonObjectSUPAve);
768 // Set output colors for fancy output
769 Color::Modifier red(Color::FG_RED);
770 Color::Modifier green(Color::FG_GREEN);
771 Color::Modifier def(Color::FG_DEFAULT);
772
773 // Start the time loop
774 while (time < finalTime)
775 {
776 time = time + dt;
777 Eigen::VectorXd res(y);
778 res.setZero();
779 hnls.solve(y);
780 newtonObjectSUPAve.operator()(y, res);
781 Eigen::VectorXd tv;
782 Eigen::VectorXd aDer;
783 aDer = (y.head(Nphi_u) - newtonObjectSUPAve.y_old.head(Nphi_u)) / dt;
784 tv.resize(dimA);
785
786 switch (interChoice)
787 {
788 case 1:
789 tv << y.segment(firstRBFInd, dimA);
790 break;
791
792 case 2:
793 tv << muStar, y.segment(firstRBFInd, dimA - muStar.size());
794 break;
795
796 case 3:
797 tv << y.segment(firstRBFInd, dimA / 2), aDer.segment(firstRBFInd, dimA / 2);
798 break;
799
800 case 4:
801 tv << muStar, y.segment(firstRBFInd, (dimA - muStar.size()) / 2),
802 aDer.segment(firstRBFInd, (dimA - muStar.size()) / 2);
803 break;
804
805 default:
806 tv << y.segment(firstRBFInd, dimA);
807 break;
808 }
809
810 for (int i = 0; i < nphiNut; i++)
811 {
812 newtonObjectSUPAve.gNut(i) = problem->rbfSplines[i]->eval(tv);
813 }
814
815 // Change initial condition for the lifting function
816 if (problem->bcMethod == "lift")
817 {
818 for (int j = 0; j < N_BC; j++)
819 {
820 y(j) = vel_now(j, 0);
821 }
822 }
823
824 newtonObjectSUPAve.yOldOld = newtonObjectSUPAve.y_old;
825 newtonObjectSUPAve.y_old = y;
826 std::cout << "################## Online solve N° " << count_online_solve <<
827 " ##################" << std::endl;
828 Info << "Time = " << time << endl;
829 std::cout << "Solving for the parameter: " << vel_now << std::endl;
830
831 if (res.norm() < 1e-5)
832 {
833 std::cout << green << "|F(x)| = " << res.norm() << " - Minimun reached in " <<
834 hnls.iter << " iterations " << def << std::endl << std::endl;
835 }
836 else
837 {
838 std::cout << red << "|F(x)| = " << res.norm() << " - Minimun reached in " <<
839 hnls.iter << " iterations " << def << std::endl << std::endl;
840 }
841
842 count_online_solve += 1;
843 tmp_sol(0) = time;
844 tmp_sol.col(0).tail(y.rows()) = y;
845
846 if (counter == nextStore)
847 {
848 if (counter2 >= online_solution.size())
849 {
850 online_solution.append(tmp_sol);
851 }
852 else
853 {
854 online_solution[counter2] = tmp_sol;
855 }
856
857 rbfCoeffMat(0, counter2) = time;
858 rbfCoeffMat.block(1, counter2, nphiNut, 1) = newtonObjectSUPAve.gNut;
859 nextStore += numberOfStores;
860 counter2 ++;
861 }
862
863 counter ++;
864 }
865
866 // Save the solution
867 ITHACAstream::exportMatrix(online_solution, "red_coeff", "python",
868 "./ITHACAoutput/red_coeff");
869 ITHACAstream::exportMatrix(online_solution, "red_coeff", "matlab",
870 "./ITHACAoutput/red_coeff");
871 count_online_solve += 1;
872}
873
874// * * * * * * * * * * * * * * * Solve Functions * * * * * * * * * * * * * //
875void ReducedUnsteadyNSTurb::solveOnlinePPE(Eigen::MatrixXd vel)
876{
877 M_Assert(exportEvery >= dt,
878 "The time step dt must be smaller than exportEvery.");
879 M_Assert(storeEvery >= dt,
880 "The time step dt must be smaller than storeEvery.");
881 M_Assert(ITHACAutilities::isInteger(storeEvery / dt) == true,
882 "The variable storeEvery must be an integer multiple of the time step dt.");
883 M_Assert(ITHACAutilities::isInteger(exportEvery / dt) == true,
884 "The variable exportEvery must be an integer multiple of the time step dt.");
885 M_Assert(ITHACAutilities::isInteger(exportEvery / storeEvery) == true,
886 "The variable exportEvery must be an integer multiple of the variable storeEvery.");
887 int numberOfStores = round(storeEvery / dt);
888
889 if (problem->bcMethod == "lift")
890 {
891 vel_now = setOnlineVelocity(vel);
892 }
893 else if (problem->bcMethod == "penalty")
894 {
895 vel_now = vel;
896 }
897
898 // Create and resize the solution vector
899 y.resize(Nphi_u + Nphi_p, 1);
900 y.setZero();
901 // Set Initial Conditions
902 y.head(Nphi_u) = initCond.col(0).head(Nphi_u);
903 y.tail(Nphi_p) = initCond.col(0).tail(Nphi_p);
904 int nextStore = 0;
905 int counter2 = 0;
906
907 // Change initial condition for the lifting function
908 if (problem->bcMethod == "lift")
909 {
910 for (int j = 0; j < N_BC; j++)
911 {
912 y(j) = vel_now(j, 0);
913 }
914 }
915
916 int firstRBFInd;
917
918 if (skipLift == true && problem->bcMethod == "lift")
919 {
920 firstRBFInd = N_BC;
921 }
922 else
923 {
924 firstRBFInd = 0;
925 }
926
927 // Set some properties of the newton object
928 newtonObjectPPE.nu = nu;
929 newtonObjectPPE.y_old = y;
930 newtonObjectPPE.yOldOld = newtonObjectPPE.y_old;
931 newtonObjectPPE.dt = dt;
932 newtonObjectPPE.bc.resize(N_BC);
933 newtonObjectPPE.tauU = tauU;
934 newtonObjectPPE.gNut = nut0;
935
936 for (int j = 0; j < N_BC; j++)
937 {
938 newtonObjectPPE.bc(j) = vel_now(j, 0);
939 }
940
941 // Set number of online solutions
942 int Ntsteps = static_cast<int>((finalTime - tstart) / dt);
943 int onlineSize = static_cast<int>(Ntsteps / numberOfStores);
944 online_solution.resize(onlineSize);
945 rbfCoeffMat.resize(nphiNut + 1, onlineSize + 3);
946 // Set the initial time
947 time = tstart;
948 // Counting variable
949 int counter = 0;
950 // Create vectpr to store temporal solution and save initial condition as first solution
951 Eigen::MatrixXd tmp_sol(Nphi_u + Nphi_p + 1, 1);
952 tmp_sol(0) = time;
953 tmp_sol.col(0).tail(y.rows()) = y;
954
955 // This part up to the while loop is just to compute the eddy viscosity field at time = startTime if time is not equal to zero
956 if ((time != 0) || (startFromZero == true))
957 {
958 online_solution[counter] = tmp_sol;
959 counter ++;
960 rbfCoeffMat(0, counter2) = time;
961 rbfCoeffMat.block(1, counter2, nphiNut, 1) = nut0;
962 counter2++;
963 nextStore += numberOfStores;
964 }
965
966 // Create nonlinear solver object
967 Eigen::HybridNonLinearSolver<newtonUnsteadyNSTurbPPE> hnls(newtonObjectPPE);
968 // Set output colors for fancy output
969 Color::Modifier red(Color::FG_RED);
970 Color::Modifier green(Color::FG_GREEN);
971 Color::Modifier def(Color::FG_DEFAULT);
972
973 // Start the time loop
974 while (time < finalTime)
975 {
976 time = time + dt;
977 Eigen::VectorXd res(y);
978 res.setZero();
979 hnls.solve(y);
980 newtonObjectPPE.operator()(y, res);
981 Eigen::VectorXd tv;
982 Eigen::VectorXd aDer;
983 aDer = (y.head(Nphi_u) - newtonObjectPPE.y_old.head(Nphi_u)) / dt;
984 tv.resize(dimA);
985
986 switch (interChoice)
987 {
988 case 1:
989 tv << y.segment(firstRBFInd, dimA);
990 break;
991
992 case 2:
993 tv << muStar, y.segment(firstRBFInd, dimA - muStar.size());
994 break;
995
996 case 3:
997 tv << y.segment(firstRBFInd, dimA / 2), aDer.segment(firstRBFInd, dimA / 2);
998 break;
999
1000 case 4:
1001 tv << muStar, y.segment(firstRBFInd, (dimA - muStar.size()) / 2),
1002 aDer.segment(firstRBFInd, (dimA - muStar.size()) / 2);
1003 break;
1004
1005 default:
1006 tv << y.segment(firstRBFInd, dimA);
1007 break;
1008 }
1009
1010 for (int i = 0; i < nphiNut; i++)
1011 {
1012 newtonObjectPPE.gNut(i) = problem->rbfSplines[i]->eval(tv);
1013 }
1014
1015 newtonObjectPPE.operator()(y, res);
1016 newtonObjectPPE.yOldOld = newtonObjectPPE.y_old;
1017 newtonObjectPPE.y_old = y;
1018 std::cout << "################## Online solve N° " << count_online_solve <<
1019 " ##################" << std::endl;
1020 Info << "Time = " << time << endl;
1021 std::cout << "Solving for the parameter: " << vel_now << std::endl;
1022
1023 if (res.norm() < 1e-5)
1024 {
1025 std::cout << green << "|F(x)| = " << res.norm() << " - Minimun reached in " <<
1026 hnls.iter << " iterations " << def << std::endl << std::endl;
1027 }
1028 else
1029 {
1030 std::cout << red << "|F(x)| = " << res.norm() << " - Minimun reached in " <<
1031 hnls.iter << " iterations " << def << std::endl << std::endl;
1032 }
1033
1034 count_online_solve += 1;
1035 tmp_sol(0) = time;
1036 tmp_sol.col(0).tail(y.rows()) = y;
1037
1038 if (counter == nextStore)
1039 {
1040 if (counter2 >= online_solution.size())
1041 {
1042 online_solution.append(tmp_sol);
1043 }
1044 else
1045 {
1046 online_solution[counter2] = tmp_sol;
1047 }
1048
1049 rbfCoeffMat(0, counter2) = time;
1050 rbfCoeffMat.block(1, counter2, nphiNut, 1) = newtonObjectPPE.gNut;
1051 nextStore += numberOfStores;
1052 counter2 ++;
1053 }
1054
1055 counter ++;
1056 }
1057
1058 // Save the solution
1059 ITHACAstream::exportMatrix(online_solution, "red_coeff", "python",
1060 "./ITHACAoutput/red_coeff");
1061 ITHACAstream::exportMatrix(online_solution, "red_coeff", "matlab",
1062 "./ITHACAoutput/red_coeff");
1063 count_online_solve += 1;
1064}
1065
1066void ReducedUnsteadyNSTurb::solveOnlinePPEAve(Eigen::MatrixXd vel)
1067{
1068 M_Assert(exportEvery >= dt,
1069 "The time step dt must be smaller than exportEvery.");
1070 M_Assert(storeEvery >= dt,
1071 "The time step dt must be smaller than storeEvery.");
1072 M_Assert(ITHACAutilities::isInteger(storeEvery / dt) == true,
1073 "The variable storeEvery must be an integer multiple of the time step dt.");
1074 M_Assert(ITHACAutilities::isInteger(exportEvery / dt) == true,
1075 "The variable exportEvery must be an integer multiple of the time step dt.");
1076 M_Assert(ITHACAutilities::isInteger(exportEvery / storeEvery) == true,
1077 "The variable exportEvery must be an integer multiple of the variable storeEvery.");
1078 int numberOfStores = round(storeEvery / dt);
1079
1080 if (problem->bcMethod == "lift")
1081 {
1082 vel_now = setOnlineVelocity(vel);
1083 }
1084 else if (problem->bcMethod == "penalty")
1085 {
1086 vel_now = vel;
1087 }
1088
1089 // Create and resize the solution vector
1090 y.resize(Nphi_u + Nphi_p, 1);
1091 y.setZero();
1092 // Set Initial Conditions
1093 y.head(Nphi_u) = initCond.col(0).head(Nphi_u);
1094 y.tail(Nphi_p) = initCond.col(0).tail(Nphi_p);
1095 int nextStore = 0;
1096 int counter2 = 0;
1097
1098 // Change initial condition for the lifting function
1099 if (problem->bcMethod == "lift")
1100 {
1101 for (int j = 0; j < N_BC; j++)
1102 {
1103 y(j) = vel_now(j, 0);
1104 }
1105 }
1106
1107 int firstRBFInd;
1108
1109 if (skipLift == true && problem->bcMethod == "lift")
1110 {
1111 firstRBFInd = N_BC;
1112 }
1113 else
1114 {
1115 firstRBFInd = 0;
1116 }
1117
1118 // Set some properties of the newton object
1119 newtonObjectPPEAve.nu = nu;
1120 newtonObjectPPEAve.y_old = y;
1121 newtonObjectPPEAve.yOldOld = newtonObjectPPEAve.y_old;
1122 newtonObjectPPEAve.dt = dt;
1123 newtonObjectPPEAve.bc.resize(N_BC);
1124 newtonObjectPPEAve.tauU = tauU;
1125 newtonObjectPPEAve.gNut = nut0;
1126 newtonObjectPPEAve.gNutAve = gNutAve;
1127
1128 for (int j = 0; j < N_BC; j++)
1129 {
1130 newtonObjectPPEAve.bc(j) = vel_now(j, 0);
1131 }
1132
1133 // Set number of online solutions
1134 int Ntsteps = static_cast<int>((finalTime - tstart) / dt);
1135 int onlineSize = static_cast<int>(Ntsteps / numberOfStores);
1136 online_solution.resize(onlineSize);
1137 rbfCoeffMat.resize(nphiNut + 1, onlineSize + 3);
1138 // Set the initial time
1139 time = tstart;
1140 // Counting variable
1141 int counter = 0;
1142 // Create vectpr to store temporal solution and save initial condition as first solution
1143 Eigen::MatrixXd tmp_sol(Nphi_u + Nphi_p + 1, 1);
1144 tmp_sol(0) = time;
1145 tmp_sol.col(0).tail(y.rows()) = y;
1146
1147 // This part up to the while loop is just to compute the eddy viscosity field at time = startTime if time is not equal to zero
1148 if ((time != 0) || (startFromZero == true))
1149 {
1150 online_solution[counter] = tmp_sol;
1151 counter ++;
1152 rbfCoeffMat(0, counter2) = time;
1153 rbfCoeffMat.block(1, counter2, nphiNut, 1) = nut0;
1154 counter2++;
1155 nextStore += numberOfStores;
1156 }
1157
1158 // Create nonlinear solver object
1159 Eigen::HybridNonLinearSolver<newtonUnsteadyNSTurbPPEAve> hnls(
1160 newtonObjectPPEAve);
1161 // Set output colors for fancy output
1162 Color::Modifier red(Color::FG_RED);
1163 Color::Modifier green(Color::FG_GREEN);
1164 Color::Modifier def(Color::FG_DEFAULT);
1165
1166 // Start the time loop
1167 while (time < finalTime)
1168 {
1169 time = time + dt;
1170 Eigen::VectorXd res(y);
1171 res.setZero();
1172 hnls.solve(y);
1173 newtonObjectPPEAve.operator()(y, res);
1174 Eigen::VectorXd tv;
1175 Eigen::VectorXd aDer;
1176 aDer = (y.head(Nphi_u) - newtonObjectPPEAve.y_old.head(Nphi_u)) / dt;
1177 tv.resize(dimA);
1178
1179 switch (interChoice)
1180 {
1181 case 1:
1182 tv << y.segment(firstRBFInd, dimA);
1183 break;
1184
1185 case 2:
1186 tv << muStar, y.segment(firstRBFInd, dimA - muStar.size());
1187 break;
1188
1189 case 3:
1190 tv << y.segment(firstRBFInd, dimA / 2), aDer.segment(firstRBFInd, dimA / 2);
1191 break;
1192
1193 case 4:
1194 tv << muStar, y.segment(firstRBFInd, (dimA - muStar.size()) / 2),
1195 aDer.segment(firstRBFInd, (dimA - muStar.size()) / 2);
1196 break;
1197
1198 default:
1199 tv << y.segment(firstRBFInd, dimA);
1200 break;
1201 }
1202
1203 for (int i = 0; i < nphiNut; i++)
1204 {
1205 newtonObjectPPEAve.gNut(i) = problem->rbfSplines[i]->eval(tv);
1206 }
1207
1208 // Change initial condition for the lifting function
1209 if (problem->bcMethod == "lift")
1210 {
1211 for (int j = 0; j < N_BC; j++)
1212 {
1213 y(j) = vel_now(j, 0);
1214 }
1215 }
1216
1217 newtonObjectPPEAve.yOldOld = newtonObjectPPEAve.y_old;
1218 newtonObjectPPEAve.y_old = y;
1219 std::cout << "################## Online solve N° " << count_online_solve <<
1220 " ##################" << std::endl;
1221 Info << "Time = " << time << endl;
1222 std::cout << "Solving for the parameter: " << vel_now << std::endl;
1223
1224 if (res.norm() < 1e-5)
1225 {
1226 std::cout << green << "|F(x)| = " << res.norm() << " - Minimun reached in " <<
1227 hnls.iter << " iterations " << def << std::endl << std::endl;
1228 }
1229 else
1230 {
1231 std::cout << red << "|F(x)| = " << res.norm() << " - Minimun reached in " <<
1232 hnls.iter << " iterations " << def << std::endl << std::endl;
1233 }
1234
1235 count_online_solve += 1;
1236 tmp_sol(0) = time;
1237 tmp_sol.col(0).tail(y.rows()) = y;
1238
1239 if (counter == nextStore)
1240 {
1241 if (counter2 >= online_solution.size())
1242 {
1243 online_solution.append(tmp_sol);
1244 }
1245 else
1246 {
1247 online_solution[counter2] = tmp_sol;
1248 }
1249
1250 rbfCoeffMat(0, counter2) = time;
1251 rbfCoeffMat.block(1, counter2, nphiNut, 1) = newtonObjectPPEAve.gNut;
1252 nextStore += numberOfStores;
1253 counter2 ++;
1254 }
1255
1256 counter ++;
1257 }
1258
1259 // Save the solution
1260 ITHACAstream::exportMatrix(online_solution, "red_coeff", "python",
1261 "./ITHACAoutput/red_coeff");
1262 ITHACAstream::exportMatrix(online_solution, "red_coeff", "matlab",
1263 "./ITHACAoutput/red_coeff");
1264 count_online_solve += 1;
1265}
1266
1267void ReducedUnsteadyNSTurb::reconstruct(bool exportFields, fileName folder)
1268{
1269 if (exportFields)
1270 {
1271 mkDir(folder);
1272 ITHACAutilities::createSymLink(folder);
1273 }
1274
1275 int counter = 0;
1276 int nextWrite = 0;
1277 int counter2 = 1;
1278 int exportEveryIndex = round(exportEvery / storeEvery);
1279 volScalarField nutAveNow("nutAveNow", nutModes[0] * 0);
1280 List < Eigen::MatrixXd> CoeffU;
1281 List < Eigen::MatrixXd> CoeffP;
1282 List < Eigen::MatrixXd> CoeffNut;
1283 CoeffU.resize(0);
1284 CoeffP.resize(0);
1285 CoeffNut.resize(0);
1286
1287 for (int k = 0; k < problem->nutAve.size(); k++)
1288 {
1289 nutAveNow += gNutAve(k) * problem->nutAve[k];
1290 }
1291
1292 for (int i = 0; i < online_solution.size(); i++)
1293 {
1294 if (counter == nextWrite)
1295 {
1296 Eigen::MatrixXd currentUCoeff;
1297 Eigen::MatrixXd currentPCoeff;
1298 Eigen::MatrixXd currentNutCoeff;
1299 currentUCoeff = online_solution[i].block(1, 0, Nphi_u, 1);
1300 currentPCoeff = online_solution[i].bottomRows(Nphi_p);
1301 currentNutCoeff = rbfCoeffMat.block(1, i, nphiNut, 1);
1302 CoeffU.append(currentUCoeff);
1303 CoeffP.append(currentPCoeff);
1304 CoeffNut.append(currentNutCoeff);
1305 nextWrite += exportEveryIndex;
1306 }
1307
1308 counter++;
1309 }
1310
1311 volVectorField uRec("uRec", Umodes[0]);
1312 volScalarField pRec("pRec", problem->Pmodes[0]);
1313 volScalarField nutRec("nutRec", problem->nutModes[0]);
1314 uRecFields = problem->L_U_SUPmodes.reconstruct(uRec, CoeffU, "uRec");
1315 pRecFields = problem->Pmodes.reconstruct(pRec, CoeffP, "pRec");
1316 nutRecFields = problem->nutModes.reconstruct(nutRec, CoeffNut, "nutRec");
1317
1318 for (int k = 0; k < nutRecFields.size(); k++)
1319 {
1320 nutRecFields[k] += nutAveNow;
1321 }
1322
1323 if (exportFields)
1324 {
1325 ITHACAstream::exportFields(uRecFields, folder,
1326 "uRec");
1327 ITHACAstream::exportFields(pRecFields, folder,
1328 "pRec");
1329 ITHACAstream::exportFields(nutRecFields, folder,
1330 "nutRec");
1331 }
1332}
1333
1334Eigen::MatrixXd ReducedUnsteadyNSTurb::setOnlineVelocity(Eigen::MatrixXd vel)
1335{
1336 assert(problem->inletIndex.rows() == vel.rows()
1337 && "Imposed boundary conditions dimensions do not match given values matrix dimensions");
1338 Eigen::MatrixXd vel_scal;
1339 vel_scal.resize(vel.rows(), vel.cols());
1340
1341 for (int k = 0; k < problem->inletIndex.rows(); k++)
1342 {
1343 int p = problem->inletIndex(k, 0);
1344 int l = problem->inletIndex(k, 1);
1345 scalar area = gSum(problem->liftfield[0].mesh().magSf().boundaryField()[p]);
1346 scalar u_lf = gSum(problem->liftfield[k].mesh().magSf().boundaryField()[p] *
1347 problem->liftfield[k].boundaryField()[p]).component(l) / area;
1348 vel_scal(k, 0) = vel(k, 0) / u_lf;
1349 }
1350
1351 return vel_scal;
1352}
1353// ************************************************************************* //
1354
M_Assert
#define M_Assert(Expr, Msg)
Definition ITHACAassert.H:46
ReducedUnsteadyNSTurb.H
Header file of the ReducedUnsteadyNSTurb class.
Color::Modifier
Class to change color to the output stream.
Definition colormod.H:24
Modes::reconstruct
GeometricField< Type, PatchField, GeoMesh > reconstruct(GeometricField< Type, PatchField, GeoMesh > &inputField, Eigen::MatrixXd Coeff, word Name)
Function to reconstruct the solution starting from the coefficients, in this case the field is passed...
Definition Modes.C:275
ReducedUnsteadyNSTurb::solveOnlinePPEAve
void solveOnlinePPEAve(Eigen::MatrixXd velNow)
Method to perform an online solve using a PPE stabilisation method with the use of the average splitt...
Definition ReducedUnsteadyNSTurb.C:1066
ReducedUnsteadyNSTurb::ReducedUnsteadyNSTurb
ReducedUnsteadyNSTurb()
Construct Null.
Definition ReducedUnsteadyNSTurb.C:41
ReducedUnsteadyNSTurb::nphiNut
int nphiNut
Number of viscosity modes.
Definition ReducedUnsteadyNSTurb.H:214
ReducedUnsteadyNSTurb::newtonObjectSUP
newtonUnsteadyNSTurbSUP newtonObjectSUP
Function object to call the non linear solver sup approach.
Definition ReducedUnsteadyNSTurb.H:232
ReducedUnsteadyNSTurb::newtonObjectSUPAve
newtonUnsteadyNSTurbSUPAve newtonObjectSUPAve
Function object to call the non linear solver sup approach with the splitting of the eddy viscosity.
Definition ReducedUnsteadyNSTurb.H:235
ReducedUnsteadyNSTurb::nutRecFields
PtrList< volScalarField > nutRecFields
Reconstructed eddy viscosity fields list.
Definition ReducedUnsteadyNSTurb.H:256
ReducedUnsteadyNSTurb::solveOnlineSUPAve
void solveOnlineSUPAve(Eigen::MatrixXd velNow)
Method to perform an online solve using a supremizer stabilisation method with the use of the average...
Definition ReducedUnsteadyNSTurb.C:673
ReducedUnsteadyNSTurb::rbfCoeffMat
Eigen::MatrixXd rbfCoeffMat
The matrix of the eddy viscosity RBF interoplated coefficients.
Definition ReducedUnsteadyNSTurb.H:223
ReducedUnsteadyNSTurb::skipLift
bool skipLift
Interpolation boolean variable to skip lifting functions.
Definition ReducedUnsteadyNSTurb.H:253
ReducedUnsteadyNSTurb::dimA
int dimA
Dimension of the interpolation independent variable.
Definition ReducedUnsteadyNSTurb.H:217
ReducedUnsteadyNSTurb::muStar
Eigen::VectorXd muStar
Online parameter value.
Definition ReducedUnsteadyNSTurb.H:244
ReducedUnsteadyNSTurb::initCond
Eigen::MatrixXd initCond
The matrix of the initial velocity and pressure reduced coefficients.
Definition ReducedUnsteadyNSTurb.H:226
ReducedUnsteadyNSTurb::gNutAve
Eigen::VectorXd gNutAve
The reduced average vector for viscosity coefficients.
Definition ReducedUnsteadyNSTurb.H:247
ReducedUnsteadyNSTurb::solveOnlineSUP
void solveOnlineSUP(Eigen::MatrixXd velNow)
Method to perform an online solve using a supremizer stabilisation method.
Definition ReducedUnsteadyNSTurb.C:467
ReducedUnsteadyNSTurb::interChoice
int interChoice
Interpolation independent variable choice.
Definition ReducedUnsteadyNSTurb.H:250
ReducedUnsteadyNSTurb::newtonObjectPPE
newtonUnsteadyNSTurbPPE newtonObjectPPE
Function object to call the non linear solver PPE approach.
Definition ReducedUnsteadyNSTurb.H:238
ReducedUnsteadyNSTurb::reconstruct
void reconstruct(bool exportFields=false, fileName folder="./online_rec")
Method to reconstruct the solutions from an online solve with any of the two techniques SUP or the PP...
Definition ReducedUnsteadyNSTurb.C:1267
ReducedUnsteadyNSTurb::nutModes
PtrList< volScalarField > nutModes
List of pointers to store the modes for the eddy viscosity.
Definition ReducedUnsteadyNSTurb.H:220
ReducedUnsteadyNSTurb::solveOnlinePPE
void solveOnlinePPE(Eigen::MatrixXd velNow)
Method to perform an online solve using a PPE stabilisation method.
Definition ReducedUnsteadyNSTurb.C:875
ReducedUnsteadyNSTurb::problem
UnsteadyNSTurb * problem
Pointer to the FOM problem.
Definition ReducedUnsteadyNSTurb.H:211
ReducedUnsteadyNSTurb::newtonObjectPPEAve
newtonUnsteadyNSTurbPPEAve newtonObjectPPEAve
Function object to call the non linear solver PPE approach with the splitting of the eddy viscosity.
Definition ReducedUnsteadyNSTurb.H:241
ReducedUnsteadyNSTurb::setOnlineVelocity
Eigen::MatrixXd setOnlineVelocity(Eigen::MatrixXd vel)
Sets the online velocity.
Definition ReducedUnsteadyNSTurb.C:1334
ReducedUnsteadyNSTurb::nut0
Eigen::VectorXd nut0
The initial eddy viscosity reduced coefficients.
Definition ReducedUnsteadyNSTurb.H:229
UnsteadyNSTurb
Implementation of a parametrized full order unsteady NS problem and preparation of the the reduced ...
Definition UnsteadyNSTurb.H:68
UnsteadyNSTurb::cTotalTensor
Eigen::Tensor< double, 3 > cTotalTensor
Turbulent total viscosity tensor.
Definition UnsteadyNSTurb.H:125
UnsteadyNSTurb::nutAve
PtrList< volScalarField > nutAve
List of for eddy viscosity time-averaged fields.
Definition UnsteadyNSTurb.H:86
UnsteadyNSTurb::rbfSplines
std::vector< SPLINTER::RBFSpline * > rbfSplines
Create a samples for interpolation.
Definition UnsteadyNSTurb.H:92
UnsteadyNSTurb::interChoice
label interChoice
Interpolation independent variable choice.
Definition UnsteadyNSTurb.H:155
UnsteadyNSTurb::nutModes
volScalarModes nutModes
List of POD modes for eddy viscosity.
Definition UnsteadyNSTurb.H:83
UnsteadyNSTurb::cTotalPPEAveTensor
Eigen::Tensor< double, 3 > cTotalPPEAveTensor
Turbulent total average viscosity tensor for the splitting approach in the PPE equation.
Definition UnsteadyNSTurb.H:134
UnsteadyNSTurb::nutFields
PtrList< volScalarField > nutFields
List of snapshots for the solution for eddy viscosity.
Definition UnsteadyNSTurb.H:80
UnsteadyNSTurb::cTotalPPETensor
Eigen::Tensor< double, 3 > cTotalPPETensor
Turbulent total viscosity tensor in the PPE equation.
Definition UnsteadyNSTurb.H:131
UnsteadyNSTurb::bTotalMatrix
Eigen::MatrixXd bTotalMatrix
Total B Matrix.
Definition UnsteadyNSTurb.H:137
UnsteadyNSTurb::nNutModes
label nNutModes
Number of viscoisty modes used for the projection.
Definition UnsteadyNSTurb.H:152
UnsteadyNSTurb::cTotalAveTensor
Eigen::Tensor< double, 3 > cTotalAveTensor
Turbulent total average viscosity tensor for the splitting approach.
Definition UnsteadyNSTurb.H:128
UnsteadyNSTurb::velRBF
Eigen::MatrixXd velRBF
Velocity coefficients for RBF interpolation.
Definition UnsteadyNSTurb.H:143
reducedSteadyNS::y
Eigen::VectorXd y
Vector to store the solution during the Newton procedure.
Definition ReducedSteadyNS.H:118
reducedSteadyNS::Pmodes
PtrList< volScalarField > Pmodes
List of pointers to store the modes for pressure.
Definition ReducedSteadyNS.H:130
reducedSteadyNS::online_solution
List< Eigen::MatrixXd > online_solution
List of Eigen matrices to store the online solution.
Definition ReducedSteadyNS.H:124
reducedSteadyNS::Nphi_p
int Nphi_p
Number of pressure modes.
Definition ReducedSteadyNS.H:160
reducedSteadyNS::nu
scalar nu
Reduced viscosity in case of parametrized viscosity.
Definition ReducedSteadyNS.H:121
reducedSteadyNS::tauU
Eigen::MatrixXd tauU
Penalty Factor.
Definition ReducedSteadyNS.H:172
reducedSteadyNS::pRecFields
PtrList< volScalarField > pRecFields
Reconstructed pressure fields list.
Definition ReducedSteadyNS.H:145
reducedSteadyNS::count_online_solve
int count_online_solve
Counter to count the online solutions.
Definition ReducedSteadyNS.H:169
reducedSteadyNS::vel_now
Eigen::MatrixXd vel_now
Online inlet velocity vector.
Definition ReducedSteadyNS.H:106
reducedSteadyNS::Umodes
PtrList< volVectorField > Umodes
List of pointers to store the modes for velocity.
Definition ReducedSteadyNS.H:127
reducedSteadyNS::uRecFields
PtrList< volVectorField > uRecFields
Recontructed velocity fields list.
Definition ReducedSteadyNS.H:148
reducedSteadyNS::N_BC
int N_BC
Number of parametrized boundary conditions.
Definition ReducedSteadyNS.H:166
reducedSteadyNS::Nphi_u
int Nphi_u
Number of velocity modes.
Definition ReducedSteadyNS.H:157
reducedUnsteadyNS::startFromZero
bool startFromZero
A variable for starting solving the reduced system from zero or not.
Definition ReducedUnsteadyNS.H:154
reducedUnsteadyNS::exportEvery
double exportEvery
A variable for exporting the fields.
Definition ReducedUnsteadyNS.H:151
reducedUnsteadyNS::finalTime
scalar finalTime
Scalar to store the final time if the online simulation.
Definition ReducedUnsteadyNS.H:142
reducedUnsteadyNS::tstart
scalar tstart
Scalar to store the initial time if the online simulation.
Definition ReducedUnsteadyNS.H:145
reducedUnsteadyNS::time
scalar time
Scalar to store the current time.
Definition ReducedUnsteadyNS.H:136
reducedUnsteadyNS::dt
double dt
Scalar to store the time increment.
Definition ReducedUnsteadyNS.H:139
reducedUnsteadyNS::storeEvery
double storeEvery
A variable for storing the reduced coefficients.
Definition ReducedUnsteadyNS.H:148
reductionProblem::inletIndex
Eigen::MatrixXi inletIndex
Matrix that contains informations about the inlet boundaries.
Definition reductionProblem.H:137
steadyNS::bcVelVec
List< Eigen::MatrixXd > bcVelVec
Boundary term for penalty method - vector.
Definition steadyNS.H:218
steadyNS::NPmodes
label NPmodes
Number of pressure modes used for the projection.
Definition steadyNS.H:143
steadyNS::BC1_matrix
Eigen::MatrixXd BC1_matrix
PPE BC1.
Definition steadyNS.H:182
steadyNS::BC3_matrix
Eigen::MatrixXd BC3_matrix
PPE BC3.
Definition steadyNS.H:191
steadyNS::C_tensor
Eigen::Tensor< double, 3 > C_tensor
Diffusion term.
Definition steadyNS.H:167
steadyNS::Pfield
PtrList< volScalarField > Pfield
List of pointers used to form the pressure snapshots matrix.
Definition steadyNS.H:86
steadyNS::supmodes
volVectorModes supmodes
List of pointers used to form the supremizer modes.
Definition steadyNS.H:104
steadyNS::Umodes
volVectorModes Umodes
List of pointers used to form the velocity modes.
Definition steadyNS.H:101
steadyNS::Ufield
PtrList< volVectorField > Ufield
List of pointers used to form the velocity snapshots matrix.
Definition steadyNS.H:89
steadyNS::NUmodes
label NUmodes
Number of velocity modes used for the projection.
Definition steadyNS.H:140
steadyNS::liftfield
PtrList< volVectorField > liftfield
List of pointers used to form the list of lifting functions.
Definition steadyNS.H:110
steadyNS::L_U_SUPmodes
volVectorModes L_U_SUPmodes
List of pointers containing the total number of lift, supremizer and velocity modes.
Definition steadyNS.H:116
steadyNS::B_matrix
Eigen::MatrixXd B_matrix
Diffusion term.
Definition steadyNS.H:157
steadyNS::D_matrix
Eigen::MatrixXd D_matrix
Laplacian term PPE.
Definition steadyNS.H:173
steadyNS::NSUPmodes
label NSUPmodes
Number of supremizer modes used for the projection.
Definition steadyNS.H:146
steadyNS::gTensor
Eigen::Tensor< double, 3 > gTensor
Divergence of momentum PPE.
Definition steadyNS.H:179
steadyNS::K_matrix
Eigen::MatrixXd K_matrix
Gradient of pressure matrix.
Definition steadyNS.H:163
steadyNS::bcVelMat
List< Eigen::MatrixXd > bcVelMat
Boundary term for penalty method - matrix.
Definition steadyNS.H:221
steadyNS::P_matrix
Eigen::MatrixXd P_matrix
Div of velocity.
Definition steadyNS.H:170
steadyNS::M_matrix
Eigen::MatrixXd M_matrix
Mass Matrix.
Definition steadyNS.H:160
steadyNS::bc2Tensor
Eigen::Tensor< double, 3 > bc2Tensor
PPE BC2.
Definition steadyNS.H:188
steadyNS::Pmodes
volScalarModes Pmodes
List of pointers used to form the pressure modes.
Definition steadyNS.H:98
steadyNS::bcMethod
word bcMethod
Boundary Method.
Definition steadyNS.H:317
unsteadyNS::timeDerivativeSchemeOrder
word timeDerivativeSchemeOrder
Definition unsteadyNS.H:99
Color::FG_GREEN
@ FG_GREEN
Definition colormod.H:14
Color::FG_DEFAULT
@ FG_DEFAULT
Definition colormod.H:16
Color::FG_RED
@ FG_RED
Definition colormod.H:13
Eigen::SliceFromTensor
Matrix< VectorType, Dynamic, Dynamic > SliceFromTensor(Eigen::Tensor< VectorType, 3 > &tensor, label dim, label index1)
Definition EigenFunctions.H:477
ITHACAstream::exportFields
void exportFields(PtrList< GeometricField< Type, PatchField, GeoMesh > > &field, word folder, word fieldname)
Function to export a scalar of vector field.
Definition ITHACAstream.C:818
ITHACAstream::exportMatrix
void exportMatrix(Eigen::Matrix< T, -1, dim > &matrix, word Name, word type, word folder)
Export the reduced matrices in numpy (type=python), matlab (type=matlab) and txt (type=eigen) format ...
Definition ITHACAstream.C:55
ITHACAutilities::getCoeffs
Eigen::VectorXd getCoeffs(GeometricField< Type, PatchField, GeoMesh > &snapshot, PtrList< GeometricField< Type, PatchField, GeoMesh > > &modes, label Nmodes, bool consider_volumes)
Projects a snapshot on a basis function and gets the coefficients of the projection.
Definition ITHACAcoeffsMass.C:317
ITHACAutilities::isInteger
bool isInteger(double ratio)
This function checks if ratio is an integer.
Definition ITHACAutilities.C:74
ITHACAutilities::createSymLink
void createSymLink(word folder)
Creates symbolic links to 0, system and constant.
Definition ITHACAsystem.C:30
p
volScalarField & p
Definition createFields.H:41
i
label i
Definition pEqn.H:46
newtonUnsteadyNSTurbPPEAve::operator()
int operator()(const Eigen::VectorXd &x, Eigen::VectorXd &fvec) const
Definition ReducedUnsteadyNSTurb.C:362
newtonUnsteadyNSTurbPPEAve::df
int df(const Eigen::VectorXd &x, Eigen::MatrixXd &fjac) const
Definition ReducedUnsteadyNSTurb.C:457
newtonUnsteadyNSTurbPPE::operator()
int operator()(const Eigen::VectorXd &x, Eigen::VectorXd &fvec) const
Definition ReducedUnsteadyNSTurb.C:269
newtonUnsteadyNSTurbPPE::df
int df(const Eigen::VectorXd &x, Eigen::MatrixXd &fjac) const
Definition ReducedUnsteadyNSTurb.C:354
newtonUnsteadyNSTurbSUPAve::df
int df(const Eigen::VectorXd &x, Eigen::MatrixXd &fjac) const
Definition ReducedUnsteadyNSTurb.C:258
newtonUnsteadyNSTurbSUPAve::operator()
int operator()(const Eigen::VectorXd &x, Eigen::VectorXd &fvec) const
Definition ReducedUnsteadyNSTurb.C:182
newtonUnsteadyNSTurbSUP::df
int df(const Eigen::VectorXd &x, Eigen::MatrixXd &fjac) const
Definition ReducedUnsteadyNSTurb.C:173
newtonUnsteadyNSTurbSUP::operator()
int operator()(const Eigen::VectorXd &x, Eigen::VectorXd &fvec) const
Definition ReducedUnsteadyNSTurb.C:99
Generated by   doxygen 1.11.0