Loading...
Searching...
No Matches
ReducedUnsteadyMSR.C
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-------------------------------------------------------------------------------
12License
13 This file is part of ITHACA-FV
14 ITHACA-FV is free software: you can redistribute it and/or modify
15 it under the terms of the GNU Lesser 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 ITHACA-FV is distributed in the hope that it will be useful,
19 but WITHOUT ANY WARRANTY; without even the implied warranty of
20 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
21 GNU Lesser General Public License for more details.
22 You should have received a copy of the GNU Lesser General Public License
23 along with ITHACA-FV. If not, see <http://www.gnu.org/licenses/>.
24\*---------------------------------------------------------------------------*/
25#include "ReducedUnsteadyMSR.H"
26
27
28reducedusMSR::reducedusMSR() {}
29
30reducedusMSR::reducedusMSR(usmsrProblem& FOMproblem)
31{
32 problem = & FOMproblem;
33 N_BC = problem->inletIndex.rows();
34 N_BCt = problem->inletIndexT.rows();
35 Nphi_u = problem->B_matrix.rows();
36 Nphi_p = problem->K_matrix.cols();
37 Nphi_flux = problem->LF_matrix[0].cols();
38 Nphi_prec1 = problem->PS1_matrix.cols();
39 Nphi_prec2 = problem->PS2_matrix.cols();
40 Nphi_prec3 = problem->PS3_matrix.cols();
41 Nphi_prec4 = problem->PS4_matrix.cols();
42 Nphi_prec5 = problem->PS5_matrix.cols();
43 Nphi_prec6 = problem->PS6_matrix.cols();
44 Nphi_prec7 = problem->PS7_matrix.cols();
45 Nphi_prec8 = problem->PS8_matrix.cols();
46 Nphi_T = problem->LT_matrix.rows();
47 Nphi_dec1 = problem->LD1_matrix.rows();
48 Nphi_dec2 = problem->LD2_matrix.rows();
49 Nphi_dec3 = problem->LD3_matrix.rows();
50 Nphi_const = problem->PF_matrix[0].rows();
51 loadConstants(problem);
52
53 //load the modes function
54 for (int k = 0; k < problem->liftfield.size(); k++)
55 {
56 Umodes.append((problem->liftfield[k]).clone());
57 }
58
59 for (int k = 0; k < problem->NUmodes; k++)
60 {
61 Umodes.append((problem->Umodes[k]).clone());
62 }
63
64 for (int k = 0; k < problem->NPmodes; k++)
65 {
66 Pmodes.append((problem->Pmodes[k]).clone());
67 }
68
69 for (int k = 0; k < problem->NFluxmodes; k++)
70 {
71 Fluxmodes.append((problem->Fluxmodes[k]).clone());
72 }
73
74 for (int k = 0; k < Nphi_prec1; k++)
75 {
76 Prec1modes.append((problem->Prec1modes[k]).clone());
77 }
78
79 for (int k = 0; k < Nphi_prec2; k++)
80 {
81 Prec2modes.append((problem->Prec2modes[k]).clone());
82 }
83
84 for (int k = 0; k < Nphi_prec3; k++)
85 {
86 Prec3modes.append((problem->Prec3modes[k]).clone());
87 }
88
89 for (int k = 0; k < Nphi_prec4; k++)
90 {
91 Prec4modes.append((problem->Prec4modes[k]).clone());
92 }
93
94 for (int k = 0; k < Nphi_prec5; k++)
95 {
96 Prec5modes.append((problem->Prec5modes[k]).clone());
97 }
98
99 for (int k = 0; k < Nphi_prec6; k++)
100 {
101 Prec6modes.append((problem->Prec6modes[k]).clone());
102 }
103
104 for (int k = 0; k < Nphi_prec7; k++)
105 {
106 Prec7modes.append((problem->Prec7modes[k]).clone());
107 }
108
109 for (int k = 0; k < Nphi_prec8; k++)
110 {
111 Prec8modes.append((problem->Prec8modes[k]).clone());
112 }
113
114 for (int k = 0; k < problem->liftfieldT.size(); k++)
115 {
116 Tmodes.append((problem->liftfieldT[k]).clone());
117 }
118
119 for (int k = 0; k < problem->NTmodes; k++)
120 {
121 Tmodes.append((problem->Tmodes[k]).clone());
122 }
123
124 for (int k = 0; k < Nphi_dec1; k++)
125 {
126 Dec1modes.append((problem->Dec1modes[k]).clone());
127 }
128
129 for (int k = 0; k < Nphi_dec2; k++)
130 {
131 Dec2modes.append((problem->Dec2modes[k]).clone());
132 }
133
134 for (int k = 0; k < Nphi_dec3; k++)
135 {
136 Dec3modes.append((problem->Dec3modes[k]).clone());
137 }
138
139 for (int k = 0; k < Nphi_const; k++)
140 {
141 vmodes.append((problem->vmodes[k]).clone());
142 Dmodes.append((problem->Dmodes[k]).clone());
143 NSFmodes.append((problem->NSFmodes[k]).clone());
144 Amodes.append((problem->Amodes[k]).clone());
145 SPmodes.append((problem->SPmodes[k]).clone());
146 TXSmodes.append((problem->TXSmodes[k]).clone());
147 }
148
149 //load the snapshots
150 for (int k = 0; k < problem->Ufield.size(); k++)
151 {
152 Usnapshots.append((problem->Ufield[k]).clone());
153 Psnapshots.append((problem->Pfield[k]).clone());
154 Fluxsnapshots.append((problem->Fluxfield[k]).clone());
155 Prec1snapshots.append((problem->Prec1field[k]).clone());
156 Prec2snapshots.append((problem->Prec2field[k]).clone());
157 Prec3snapshots.append((problem->Prec3field[k]).clone());
158 Prec4snapshots.append((problem->Prec4field[k]).clone());
159 Prec5snapshots.append((problem->Prec5field[k]).clone());
160 Prec6snapshots.append((problem->Prec6field[k]).clone());
161 Prec7snapshots.append((problem->Prec7field[k]).clone());
162 Prec8snapshots.append((problem->Prec8field[k]).clone());
163 Tsnapshots.append((problem->Tfield[k]).clone());
164 Dec1snapshots.append((problem->Dec1field[k]).clone());
165 Dec2snapshots.append((problem->Dec2field[k]).clone());
166 Dec3snapshots.append((problem->Dec3field[k]).clone());
167 vsnapshots.append((problem->vFields[k]).clone());
168 Dsnapshots.append((problem->DFields[k]).clone());
169 NSFsnapshots.append((problem->NSFFields[k]).clone());
170 Asnapshots.append((problem->AFields[k]).clone());
171 SPsnapshots.append((problem->SPFields[k]).clone());
172 TXSsnapshots.append((problem->TXSFields[k]).clone());
173 }
174
175 newton_object_fd = newton_usmsr_fd(Nphi_u + Nphi_p, Nphi_u + Nphi_p,
176 FOMproblem);
177 newton_object_n = newton_usmsr_n(Nphi_flux + Nphi_prec1 + Nphi_prec2 +
178 Nphi_prec3 + Nphi_prec4 + Nphi_prec5 + Nphi_prec6 + Nphi_prec7 + Nphi_prec8,
179 Nphi_flux + Nphi_prec1 + Nphi_prec2 + Nphi_prec3 + Nphi_prec4 + Nphi_prec5 +
180 Nphi_prec6 + Nphi_prec7 + Nphi_prec8, FOMproblem);
181 newton_object_t = newton_usmsr_t(Nphi_T + Nphi_dec1 + Nphi_dec2 + Nphi_dec3,
182 Nphi_T + Nphi_dec1 + Nphi_dec2 + Nphi_dec3, FOMproblem);
183}
184
185int newton_usmsr_fd::operator()(const Eigen::VectorXd& x,
186 Eigen::VectorXd& fvec) const
187{
188 Eigen::VectorXd a_tmp(Nphi_u);
189 Eigen::VectorXd b_tmp(Nphi_p);
190 a_tmp = x.head(Nphi_u);
191 b_tmp = x.tail(Nphi_p);
192 Eigen::VectorXd a_dot(Nphi_u);
193 a_dot = (x.head(Nphi_u) - y_old.head(Nphi_u)) / dt;
195 // Convective terms
196 Eigen::MatrixXd cc(1, 1);
197 Eigen::MatrixXd gg(1, 1);
198 Eigen::MatrixXd bb(1, 1);
199 // Mom Term
200 Eigen::VectorXd M1 = problem->B_matrix * a_tmp * nu;
201 // Gradient of pressure
202 Eigen::VectorXd M2 = problem->K_matrix * b_tmp;
203 // Mass Term
204 Eigen::VectorXd M5 = problem->M_matrix * a_dot;
205 // Pressure Term
206 Eigen::VectorXd M3 = problem->D_matrix * b_tmp;
207 // BC PPE
208 Eigen::VectorXd M6 = problem->BC1_matrix * a_tmp * nu;
209 // BC PPE
210 Eigen::VectorXd M7 = problem->BC3_matrix * a_tmp * nu;
211
212 for (int i = 0; i < Nphi_u; i++)
213 {
214 cc = a_tmp.transpose() * problem->C_matrix[i] * a_tmp;
215 fvec(i) = -M5(i) + M1(i) - cc(0, 0) - M2(i);
216 }
217
218 for (int i = 0; i < Nphi_p; i++)
219 {
220 int k = i + Nphi_u;
221 gg = a_tmp.transpose() * problem->G_matrix[i] * a_tmp;
222 //bb = a_tmp.transpose() * problem->BC2_matrix[i] * a_tmp;
223 fvec(k) = M3(i, 0) + gg(0, 0) - M7(i, 0);
224 }
225
226 for (int j = 0; j < N_BC; j++)
227 {
228 fvec(j) = x(j) - BC(j);
229 }
230
231 return 0;
232}
233
234int newton_usmsr_fd::df(const Eigen::VectorXd& x,
235 Eigen::MatrixXd& fjac) const
236{
237 Eigen::NumericalDiff<newton_usmsr_fd> numDiff(* this);
238 numDiff.df(x, fjac);
239 return 0;
240}
241
242
243int newton_usmsr_n::operator()(const Eigen::VectorXd& n,
244 Eigen::VectorXd& fvecn) const
245{
246 Eigen::VectorXd c_tmp(Nphi_flux); //for flux
247 Eigen::VectorXd d1_tmp(Nphi_prec1); //for prec1
248 Eigen::VectorXd d2_tmp(Nphi_prec2); //for prec2
249 Eigen::VectorXd d3_tmp(Nphi_prec3); //for prec3
250 Eigen::VectorXd d4_tmp(Nphi_prec4); //for prec4
251 Eigen::VectorXd d5_tmp(Nphi_prec5); //for prec5
252 Eigen::VectorXd d6_tmp(Nphi_prec6); //for prec6
253 Eigen::VectorXd d7_tmp(Nphi_prec7); //for prec7
254 Eigen::VectorXd d8_tmp(Nphi_prec8); //for prec8
255 c_tmp = n.head(Nphi_flux);
256 int pos = Nphi_flux;
257 d1_tmp = n.segment(pos, Nphi_prec1);
258 pos = pos + Nphi_prec1;
259 d2_tmp = n.segment(pos, Nphi_prec2);
260 pos = pos + Nphi_prec2;
261 d3_tmp = n.segment(pos, Nphi_prec3);
262 pos = pos + Nphi_prec3;
263 d4_tmp = n.segment(pos, Nphi_prec4);
264 pos = pos + Nphi_prec4;
265 d5_tmp = n.segment(pos, Nphi_prec5);
266 pos = pos + Nphi_prec5;
267 d6_tmp = n.segment(pos, Nphi_prec6);
268 pos = pos + Nphi_prec6;
269 d7_tmp = n.segment(pos, Nphi_prec7);
270 pos = pos + Nphi_prec7;
271 d8_tmp = n.segment(pos, Nphi_prec8);
272 Eigen::VectorXd c_dot(Nphi_flux);
273 Eigen::VectorXd d1_dot(Nphi_prec1);
274 Eigen::VectorXd d2_dot(Nphi_prec2);
275 Eigen::VectorXd d3_dot(Nphi_prec3);
276 Eigen::VectorXd d4_dot(Nphi_prec4);
277 Eigen::VectorXd d5_dot(Nphi_prec5);
278 Eigen::VectorXd d6_dot(Nphi_prec6);
279 Eigen::VectorXd d7_dot(Nphi_prec7);
280 Eigen::VectorXd d8_dot(Nphi_prec8);
281 c_dot = (n.head(Nphi_flux) - w_old.head(Nphi_flux)) / dt;
282 pos = Nphi_flux;
283 d1_dot = (n.segment(pos, Nphi_prec1) - w_old.segment(pos, Nphi_prec1)) / dt;
284 pos += Nphi_prec1;
285 d2_dot = (n.segment(pos, Nphi_prec2) - w_old.segment(pos, Nphi_prec2)) / dt;
286 pos += Nphi_prec2;
287 d3_dot = (n.segment(pos, Nphi_prec3) - w_old.segment(pos, Nphi_prec3)) / dt;
288 pos += Nphi_prec3;
289 d4_dot = (n.segment(pos, Nphi_prec4) - w_old.segment(pos, Nphi_prec4)) / dt;
290 pos += Nphi_prec4;
291 d5_dot = (n.segment(pos, Nphi_prec5) - w_old.segment(pos, Nphi_prec5)) / dt;
292 pos += Nphi_prec5;
293 d6_dot = (n.segment(pos, Nphi_prec6) - w_old.segment(pos, Nphi_prec6)) / dt;
294 pos += Nphi_prec6;
295 d7_dot = (n.segment(pos, Nphi_prec7) - w_old.segment(pos, Nphi_prec7)) / dt;
296 pos += Nphi_prec7;
297 d8_dot = (n.segment(pos, Nphi_prec8) - w_old.segment(pos, Nphi_prec8)) / dt;
299 // ddt flux term
300 Eigen::VectorXd F_dot = problem->MF_matrix * c_dot * iv;
301 // Laplacian flux term
302 Eigen::MatrixXd lf(1, 1);
303 // flux production term
304 Eigen::MatrixXd pf(1, 1);
305 // flux absorption term
306 Eigen::MatrixXd af(1, 1);
307 // precursor sources
308 Eigen::VectorXd F3_1 = problem->PS1_matrix * d1_tmp * l1;
309 Eigen::VectorXd F3_2 = problem->PS2_matrix * d2_tmp * l2;
310 Eigen::VectorXd F3_3 = problem->PS3_matrix * d3_tmp * l3;
311 Eigen::VectorXd F3_4 = problem->PS4_matrix * d4_tmp * l4;
312 Eigen::VectorXd F3_5 = problem->PS5_matrix * d5_tmp * l5;
313 Eigen::VectorXd F3_6 = problem->PS6_matrix * d6_tmp * l6;
314 Eigen::VectorXd F3_7 = problem->PS7_matrix * d7_tmp * l7;
315 Eigen::VectorXd F3_8 = problem->PS8_matrix * d8_tmp * l8;
316 //ddt prec term
317 Eigen::VectorXd Pdot_1 = problem->MP1_matrix * d1_dot;
318 Eigen::VectorXd Pdot_2 = problem->MP2_matrix * d2_dot;
319 Eigen::VectorXd Pdot_3 = problem->MP3_matrix * d3_dot;
320 Eigen::VectorXd Pdot_4 = problem->MP4_matrix * d4_dot;
321 Eigen::VectorXd Pdot_5 = problem->MP5_matrix * d5_dot;
322 Eigen::VectorXd Pdot_6 = problem->MP6_matrix * d6_dot;
323 Eigen::VectorXd Pdot_7 = problem->MP7_matrix * d7_dot;
324 Eigen::VectorXd Pdot_8 = problem->MP8_matrix * d8_dot;
325 // Convective terms in prec-eq:
326 Eigen::MatrixXd pp1(1, 1);
327 Eigen::MatrixXd pp2(1, 1);
328 Eigen::MatrixXd pp3(1, 1);
329 Eigen::MatrixXd pp4(1, 1);
330 Eigen::MatrixXd pp5(1, 1);
331 Eigen::MatrixXd pp6(1, 1);
332 Eigen::MatrixXd pp7(1, 1);
333 Eigen::MatrixXd pp8(1, 1);
334 // laplacian of precursor
335 Eigen::VectorXd P1_1 = problem->LP1_matrix * d1_tmp * (nu / Sc);
336 Eigen::VectorXd P1_2 = problem->LP2_matrix * d2_tmp * (nu / Sc);
337 Eigen::VectorXd P1_3 = problem->LP3_matrix * d3_tmp * (nu / Sc);
338 Eigen::VectorXd P1_4 = problem->LP4_matrix * d4_tmp * (nu / Sc);
339 Eigen::VectorXd P1_5 = problem->LP5_matrix * d5_tmp * (nu / Sc);
340 Eigen::VectorXd P1_6 = problem->LP6_matrix * d6_tmp * (nu / Sc);
341 Eigen::VectorXd P1_7 = problem->LP7_matrix * d7_tmp * (nu / Sc);
342 Eigen::VectorXd P1_8 = problem->LP8_matrix * d8_tmp * (nu / Sc);
343 // algebric term
344 Eigen::VectorXd P2_1 = problem->MP1_matrix * d1_tmp * l1;
345 Eigen::VectorXd P2_2 = problem->MP2_matrix * d2_tmp * l2;
346 Eigen::VectorXd P2_3 = problem->MP3_matrix * d3_tmp * l3;
347 Eigen::VectorXd P2_4 = problem->MP4_matrix * d4_tmp * l4;
348 Eigen::VectorXd P2_5 = problem->MP5_matrix * d5_tmp * l5;
349 Eigen::VectorXd P2_6 = problem->MP6_matrix * d6_tmp * l6;
350 Eigen::VectorXd P2_7 = problem->MP7_matrix * d7_tmp * l7;
351 Eigen::VectorXd P2_8 = problem->MP8_matrix * d8_tmp * l8;
352 // flux source term
353 Eigen::MatrixXd fs1(1, 1);
354 Eigen::MatrixXd fs2(1, 1);
355 Eigen::MatrixXd fs3(1, 1);
356 Eigen::MatrixXd fs4(1, 1);
357 Eigen::MatrixXd fs5(1, 1);
358 Eigen::MatrixXd fs6(1, 1);
359 Eigen::MatrixXd fs7(1, 1);
360 Eigen::MatrixXd fs8(1, 1);
361
362 for (int i = 0; i < Nphi_flux; i++)
363 {
364 lf = d_c.transpose() * problem->LF_matrix[i] * c_tmp;
365 pf = nsf_c.transpose() * problem->PF_matrix[i] * c_tmp * (1 - btot);
366 af = a_c.transpose() * problem->AF_matrix[i] * c_tmp;
367 fvecn(i) = -F_dot(i) + lf(0, 0) + pf(0, 0) - af(0,
368 0) + F3_1(i) + F3_2(i) + F3_3(i) + F3_4(i) + F3_5(i) + F3_6(i) + F3_7(i) + F3_8(
369 i);
370 }
371
372 int pfvecn = Nphi_flux;
373
374 for (int i = 0; i < Nphi_prec1; i++)
375 {
376 int k = i + pfvecn;
377 pp1 = a_tmp.transpose() * problem->ST1_matrix[i] * d1_tmp;
378 fs1 = nsf_c.transpose() * problem->FS1_matrix[i] * c_tmp * b1;
379 fvecn(k) = -Pdot_1(i) - pp1(0, 0) + P1_1(i) - P2_1(i) + fs1(0, 0);
380 }
381
382 pfvecn += Nphi_prec1;
383
384 for (int i = 0; i < Nphi_prec2; i++)
385 {
386 int k = i + pfvecn;
387 pp2 = a_tmp.transpose() * problem->ST2_matrix[i] * d2_tmp;
388 fs2 = nsf_c.transpose() * problem->FS2_matrix[i] * c_tmp * b2;
389 fvecn(k) = -Pdot_2(i) - pp2(0, 0) + P1_2(i) - P2_2(i) + fs2(0, 0);
390 }
391
392 pfvecn += Nphi_prec2;
393
394 for (int i = 0; i < Nphi_prec3; i++)
395 {
396 int k = i + pfvecn;
397 pp3 = a_tmp.transpose() * problem->ST3_matrix[i] * d3_tmp;
398 fs3 = nsf_c.transpose() * problem->FS3_matrix[i] * c_tmp * b3;
399 fvecn(k) = -Pdot_3(i) - pp3(0, 0) + P1_3(i) - P2_3(i) + fs3(0, 0);
400 }
401
402 pfvecn += Nphi_prec3;
403
404 for (int i = 0; i < Nphi_prec4; i++)
405 {
406 int k = i + pfvecn;
407 pp4 = a_tmp.transpose() * problem->ST4_matrix[i] * d4_tmp;
408 fs4 = nsf_c.transpose() * problem->FS4_matrix[i] * c_tmp * b4;
409 fvecn(k) = -Pdot_4(i) - pp4(0, 0) + P1_4(i) - P2_4(i) + fs4(0, 0);
410 }
411
412 pfvecn += Nphi_prec4;
413
414 for (int i = 0; i < Nphi_prec5; i++)
415 {
416 int k = i + pfvecn;
417 pp5 = a_tmp.transpose() * problem->ST5_matrix[i] * d5_tmp;
418 fs5 = nsf_c.transpose() * problem->FS5_matrix[i] * c_tmp * b5;
419 fvecn(k) = -Pdot_5(i) - pp5(0, 0) + P1_5(i) - P2_5(i) + fs5(0, 0);
420 }
421
422 pfvecn += Nphi_prec5;
423
424 for (int i = 0; i < Nphi_prec6; i++)
425 {
426 int k = i + pfvecn;
427 pp6 = a_tmp.transpose() * problem->ST6_matrix[i] * d6_tmp;
428 fs6 = nsf_c.transpose() * problem->FS6_matrix[i] * c_tmp * b6;
429 fvecn(k) = -Pdot_6(i) - pp6(0, 0) + P1_6(i) - P2_6(i) + fs6(0, 0);
430 }
431
432 pfvecn += Nphi_prec6;
433
434 for (int i = 0; i < Nphi_prec7; i++)
435 {
436 int k = i + pfvecn;
437 pp7 = a_tmp.transpose() * problem->ST7_matrix[i] * d7_tmp;
438 fs7 = nsf_c.transpose() * problem->FS7_matrix[i] * c_tmp * b7;
439 fvecn(k) = -Pdot_7(i) - pp7(0, 0) + P1_7(i) - P2_7(i) + fs7(0, 0);
440 }
441
442 pfvecn += Nphi_prec7;
443
444 for (int i = 0; i < Nphi_prec8; i++)
445 {
446 int k = i + pfvecn;
447 pp8 = a_tmp.transpose() * problem->ST8_matrix[i] * d8_tmp;
448 fs8 = nsf_c.transpose() * problem->FS8_matrix[i] * c_tmp * b8;
449 fvecn(k) = -Pdot_8(i) - pp8(0, 0) + P1_8(i) - P2_8(i) + fs8(0, 0);
450 }
451
452 return 0;
453}
454
455int newton_usmsr_n::df(const Eigen::VectorXd& n,
456 Eigen::MatrixXd& fjacn) const
457{
458 Eigen::NumericalDiff<newton_usmsr_n> numDiff(* this);
459 numDiff.df(n, fjacn);
460 return 0;
461}
462
463int newton_usmsr_t::operator()(const Eigen::VectorXd& t,
464 Eigen::VectorXd& fvect) const
465{
466 Eigen::VectorXd e_tmp(Nphi_T); //for T
467 Eigen::VectorXd f1_tmp(Nphi_dec1); //for dec1
468 Eigen::VectorXd f2_tmp(Nphi_dec2); //for dec2
469 Eigen::VectorXd f3_tmp(Nphi_dec3); //for dec3
470 e_tmp = t.head(Nphi_T);
471 int pos = Nphi_T;
472 f1_tmp = t.segment(pos, Nphi_dec1);
473 pos += Nphi_dec1;
474 f2_tmp = t.segment(pos, Nphi_dec2);
475 pos += Nphi_dec2;
476 f3_tmp = t.segment(pos, Nphi_dec3);
477 Eigen::VectorXd e_dot(Nphi_T); //for T
478 Eigen::VectorXd f1_dot(Nphi_dec1); //for dec1
479 Eigen::VectorXd f2_dot(Nphi_dec2); //for dec2
480 Eigen::VectorXd f3_dot(Nphi_dec3); //for dec3
481 e_dot = (t.head(Nphi_T) - z_old.head(Nphi_T)) / dt;
482 pos = Nphi_T;
483 f1_dot = (t.segment(pos, Nphi_dec1) - z_old.segment(pos, Nphi_dec1)) / dt;
484 pos += Nphi_dec1;
485 f2_dot = (t.segment(pos, Nphi_dec2) - z_old.segment(pos, Nphi_dec2)) / dt;
486 pos += Nphi_dec2;
487 f3_dot = (t.segment(pos, Nphi_dec3) - z_old.segment(pos, Nphi_dec3)) / dt;
489 //ddt T term
490 Eigen::VectorXd T_dot = problem->TM_matrix * e_dot;
491 // convective term in T_eqn
492 Eigen::MatrixXd tt(1, 1);
493 // laplacian of T
494 Eigen::VectorXd T1 = problem->LT_matrix * e_tmp * nu / Pr;
495 // temp flux source
496 Eigen::MatrixXd xsf(1, 1);
497 // decay heat source term (*dli/cp)
498 Eigen::MatrixXd dhs1(1, 1);
499 Eigen::MatrixXd dhs2(1, 1);
500 Eigen::MatrixXd dhs3(1, 1);
501 //ddt dec term
502 Eigen::VectorXd DHdot_1 = problem->MD1_matrix * f1_dot;
503 Eigen::VectorXd DHdot_2 = problem->MD2_matrix * f2_dot;
504 Eigen::VectorXd DHdot_3 = problem->MD3_matrix * f3_dot;
505 // convective term in decay heat eq.
506 Eigen::MatrixXd dh1(1, 1);
507 Eigen::MatrixXd dh2(1, 1);
508 Eigen::MatrixXd dh3(1, 1);
509 //laplacian of dh
510 Eigen::VectorXd DH1_1 = problem->LD1_matrix * f1_tmp * nu / Sc;
511 Eigen::VectorXd DH1_2 = problem->LD2_matrix * f2_tmp * nu / Sc;
512 Eigen::VectorXd DH1_3 = problem->LD3_matrix * f3_tmp * nu / Sc;
513 //algebric term in dh eq
514 Eigen::VectorXd DH2_1 = problem->MD1_matrix * f1_tmp * dl1;
515 Eigen::VectorXd DH2_2 = problem->MD2_matrix * f2_tmp * dl2;
516 Eigen::VectorXd DH2_3 = problem->MD3_matrix * f3_tmp * dl3;
517 // flux source in term dh eq (*dbi)
518 Eigen::MatrixXd dfs1(1, 1);
519 Eigen::MatrixXd dfs2(1, 1);
520 Eigen::MatrixXd dfs3(1, 1);
521
522 for (int i = 0; i < Nphi_T; i++)
523 {
524 tt = a_tmp.transpose() * problem->TS_matrix[i] * e_tmp;
525 xsf = txs_c.transpose() * problem->TXS_matrix[i] * c_tmp * ((1 - dbtot) / cp);
526 dhs1 = v_c.transpose() * problem->THS1_matrix[i] * f1_tmp * (dl1 / cp);
527 dhs2 = v_c.transpose() * problem->THS2_matrix[i] * f2_tmp * (dl2 / cp);
528 dhs3 = v_c.transpose() * problem->THS3_matrix[i] * f3_tmp * (dl3 / cp);
529 fvect(i) = -T_dot(i) - tt(0, 0) + T1(i) + xsf(0, 0) + dhs1(0, 0) + dhs2(0,
530 0) + dhs3(0, 0);
531 }
532
533 int pfvect = Nphi_T;
534
535 for (int i = 0; i < Nphi_dec1; i++)
536 {
537 int k = i + pfvect;
538 dh1 = a_tmp.transpose() * problem->SD1_matrix[i] * f1_tmp;
539 dfs1 = sp_c.transpose() * problem->DFS1_matrix[i] * c_tmp * db1;
540 fvect(k) = -DHdot_1(i) - dh1(0, 0) + DH1_1(i) - DH2_1(i) + dfs1(0, 0);
541 }
542
543 pfvect += Nphi_dec1;
544
545 for (int i = 0; i < Nphi_dec2; i++)
546 {
547 int k = i + pfvect;
548 dh2 = a_tmp.transpose() * problem->SD2_matrix[i] * f2_tmp;
549 dfs2 = sp_c.transpose() * problem->DFS2_matrix[i] * c_tmp * db2;
550 fvect(k) = -DHdot_2(i) - dh2(0, 0) + DH1_2(i) - DH2_2(i) + dfs2(0, 0);
551 }
552
553 pfvect += Nphi_dec2;
554
555 for (int i = 0; i < Nphi_dec3; i++)
556 {
557 int k = i + pfvect;
558 dh3 = a_tmp.transpose() * problem->SD3_matrix[i] * f3_tmp;
559 dfs3 = sp_c.transpose() * problem->DFS3_matrix[i] * c_tmp * db3;
560 fvect(k) = -DHdot_3(i) - dh3(0, 0) + DH1_3(i) - DH2_3(i) + dfs3(0, 0);
561 }
562
563 for (int i = 0; i < N_BCt; i++)
564 {
565 fvect(i) = t(i) - BCt(i);
566 }
567
568 return 0;
569}
570
571int newton_usmsr_t::df(const Eigen::VectorXd& t,
572 Eigen::MatrixXd& fjact) const
573{
574 Eigen::NumericalDiff<newton_usmsr_t> numDiff(* this);
575 numDiff.df(t, fjact);
576 return 0;
577}
578
579
580void reducedusMSR::solveOnline(Eigen::MatrixXd vel_now,
581 Eigen::MatrixXd temp_now, Eigen::VectorXd mu_online, int startSnap)
582{
583 Info << "\n Starting online stage...\n" << endl;
584 y.resize(Nphi_u + Nphi_p, 1); //for fd
585 y.setZero();
586 w.resize(Nphi_flux + Nphi_prec1 + Nphi_prec2 + Nphi_prec3 + Nphi_prec4 +
587 Nphi_prec5 + Nphi_prec6 + Nphi_prec7 + Nphi_prec8, 1); //for n
588 w.setZero();
589 z.resize(Nphi_T + Nphi_dec1 + Nphi_dec2 + Nphi_dec3, 1); //for t
590 z.setZero();
591 int pos_w = 0;
592 int pos_z = 0;
593 y.head(Nphi_u) = ITHACAutilities::getCoeffs(Usnapshots[startSnap], Umodes);
594 y.tail(Nphi_p) = ITHACAutilities::getCoeffs(Psnapshots[startSnap], Pmodes);
595
596 for (int j = 0; j < N_BC; j++)
597 {
598 y(j) = vel_now(j, 0);
599 }
600
601 w.head(Nphi_flux) = ITHACAutilities::getCoeffs(Fluxsnapshots[startSnap],
602 Fluxmodes);
603 pos_w += Nphi_flux;
604 w.segment(pos_w, Nphi_prec1) = ITHACAutilities::getCoeffs(
605 Prec1snapshots[startSnap], Prec1modes);
606 pos_w += Nphi_prec1;
607 w.segment(pos_w, Nphi_prec2) = ITHACAutilities::getCoeffs(
608 Prec2snapshots[startSnap], Prec2modes);
609 pos_w += Nphi_prec2;
610 w.segment(pos_w, Nphi_prec3) = ITHACAutilities::getCoeffs(
611 Prec3snapshots[startSnap], Prec3modes);
612 pos_w += Nphi_prec3;
613 w.segment(pos_w, Nphi_prec4) = ITHACAutilities::getCoeffs(
614 Prec4snapshots[startSnap], Prec4modes);
615 pos_w += Nphi_prec4;
616 w.segment(pos_w, Nphi_prec5) = ITHACAutilities::getCoeffs(
617 Prec5snapshots[startSnap], Prec5modes);
618 pos_w += Nphi_prec5;
619 w.segment(pos_w, Nphi_prec6) = ITHACAutilities::getCoeffs(
620 Prec6snapshots[startSnap], Prec6modes);
621 pos_w += Nphi_prec6;
622 w.segment(pos_w, Nphi_prec7) = ITHACAutilities::getCoeffs(
623 Prec7snapshots[startSnap], Prec7modes);
624 pos_w += Nphi_prec7;
625 w.segment(pos_w, Nphi_prec8) = ITHACAutilities::getCoeffs(
626 Prec8snapshots[startSnap], Prec8modes);
627 z.head(Nphi_T) = ITHACAutilities::getCoeffs(Tsnapshots[startSnap], Tmodes);
628 pos_z += Nphi_T;
629 z.segment(pos_z, Nphi_dec1) = ITHACAutilities::getCoeffs(
630 Dec1snapshots[startSnap], Dec1modes);
631 pos_z += Nphi_dec1;
632 z.segment(pos_z, Nphi_dec2) = ITHACAutilities::getCoeffs(
633 Dec2snapshots[startSnap], Dec2modes);
634 pos_z += Nphi_dec2;
635 z.segment(pos_z, Nphi_dec3) = ITHACAutilities::getCoeffs(
636 Dec3snapshots[startSnap], Dec3modes);
637
638 for (int j = 0; j < N_BCt; j++)
639 {
640 z(j) = temp_now(j, 0);
641 }
642
643 newton_object_fd.nu = nu;
644 newton_object_fd.y_old = y;
645 newton_object_fd.dt = dt;
646 newton_object_fd.BC.resize(N_BC);
647
648 for (int j = 0; j < N_BC; j++)
649 {
650 newton_object_fd.BC(j) = vel_now(j, 0);
651 }
652
653 newton_object_n.nu = nu;
654 newton_object_n.d = d;
655 newton_object_n.m = m;
656 newton_object_n.nsf = nsf;
657 newton_object_n.Keff = Keff;
658 newton_object_n.iv = iv;
659 newton_object_n.l1 = l1;
660 newton_object_n.l2 = l2;
661 newton_object_n.l3 = l3;
662 newton_object_n.l4 = l4;
663 newton_object_n.l5 = l5;
664 newton_object_n.l6 = l6;
665 newton_object_n.l7 = l7;
666 newton_object_n.l8 = l8;
667 newton_object_n.b1 = b1;
668 newton_object_n.b2 = b2;
669 newton_object_n.b3 = b3;
670 newton_object_n.b4 = b4;
671 newton_object_n.b5 = b5;
672 newton_object_n.b6 = b6;
673 newton_object_n.b7 = b7;
674 newton_object_n.b8 = b8;
675 newton_object_n.btot = btot;
676 newton_object_n.Sc = Sc;
677 newton_object_n.w_old = w;
678 newton_object_n.dt = dt;
679 newton_object_t.nu = nu;
680 newton_object_t.Keff = Keff;
681 newton_object_t.sp = sp;
682 newton_object_t.cp = cp;
683 newton_object_t.dl1 = dl1;
684 newton_object_t.dl2 = dl2;
685 newton_object_t.dl3 = dl3;
686 newton_object_t.db1 = db1;
687 newton_object_t.db2 = db2;
688 newton_object_t.db3 = db3;
689 newton_object_t.dbtot = dbtot;
690 newton_object_t.Sc = Sc;
691 newton_object_t.Pr = Pr;
692 newton_object_t.dt = dt;
693 newton_object_t.z_old = z;
694 newton_object_t.BCt.resize(N_BCt);
695
696 for (int j = 0; j < N_BCt; j++)
697 {
698 newton_object_t.BCt(j) = temp_now(j, 0);
699 }
700
701 // Set number of online solutions
702 int Ntsteps = static_cast<int>((finalTime - tstart) / dt);
703 online_solution_fd.resize(Ntsteps);
704 online_solution_n.resize(Ntsteps);
705 online_solution_t.resize(Ntsteps);
706 online_solution_C.resize(Ntsteps);
707 // Set the initial time
708 time = tstart;
709 // Counting variable
710 int counter = 0;
711 // Create vector to store temporal solution and save initial condition as first solution
712 Eigen::MatrixXd tmp_sol_fd(Nphi_u + Nphi_p + 1, 1);
713 tmp_sol_fd(0) = time;
714 tmp_sol_fd.col(0).tail(y.rows()) = y;
715 Eigen::MatrixXd tmp_sol_n(Nphi_flux + Nphi_prec1 + Nphi_prec2 + Nphi_prec3 +
716 Nphi_prec4 + Nphi_prec5 + Nphi_prec6 + Nphi_prec7 + Nphi_prec8 + 1, 1);
717 tmp_sol_n(0) = time;
718 tmp_sol_n.col(0).tail(w.rows()) = w;
719 Eigen::MatrixXd tmp_sol_t(Nphi_T + Nphi_dec1 + Nphi_dec2 + Nphi_dec3 + 1, 1);
720 tmp_sol_t(0) = time;
721 tmp_sol_t.col(0).tail(z.rows()) = z;
722 Eigen::MatrixXd tmp_sol_C(6 * Nphi_const + 1, 1);
723 tmp_sol_C(0) = time;
724 Eigen::VectorXd v_c0;
725 Eigen::VectorXd d_c0;
726 Eigen::VectorXd nsf_c0;
727 Eigen::VectorXd a_c0;
728 Eigen::VectorXd sp_c0;
729 Eigen::VectorXd txs_c0;
730 v_c0.resize(Nphi_const);
731 d_c0.resize(Nphi_const);
732 nsf_c0.resize(Nphi_const);
733 a_c0.resize(Nphi_const);
734 sp_c0.resize(Nphi_const);
735 txs_c0.resize(Nphi_const);
736 std::vector<double> tv0;
737 tv0.resize(mu_online.size() + 1);
738 tv0[0] = time;
739
740 for (int k = 1; k < tv0.size(); k++)
741 {
742 tv0[k] = mu_online(k - 1);
743 }
744
745 for (int i = 0; i < Nphi_const; i++)
746 {
747 v_c0(i) = problem->rbfsplines_v[i]->eval(tv0);
748 d_c0(i) = problem->rbfsplines_D[i]->eval(tv0);
749 nsf_c0(i) = problem->rbfsplines_NSF[i]->eval(tv0);
750 a_c0(i) = problem->rbfsplines_A[i]->eval(tv0);
751 sp_c0(i) = problem->rbfsplines_SP[i]->eval(tv0);
752 txs_c0(i) = problem->rbfsplines_TXS[i]->eval(tv0);
753 }
754
755 int pos_c = 1;
756 tmp_sol_C.col(0).segment(pos_c, Nphi_const) = v_c0;
757 pos_c += Nphi_const;
758 tmp_sol_C.col(0).segment(pos_c, Nphi_const) = d_c0;
759 pos_c += Nphi_const;
760 tmp_sol_C.col(0).segment(pos_c, Nphi_const) = nsf_c0;
761 pos_c += Nphi_const;
762 tmp_sol_C.col(0).segment(pos_c, Nphi_const) = a_c0;
763 pos_c += Nphi_const;
764 tmp_sol_C.col(0).segment(pos_c, Nphi_const) = sp_c0;
765 pos_c += Nphi_const;
766 tmp_sol_C.col(0).segment(pos_c, Nphi_const) = txs_c0;
767 online_solution_fd[counter] = tmp_sol_fd;
768 online_solution_n[counter] = tmp_sol_n;
769 online_solution_t[counter] = tmp_sol_t;
770 online_solution_C[counter] = tmp_sol_C;
771 counter++;
772 // Create nonlinear solver object
773 Eigen::HybridNonLinearSolver<newton_usmsr_fd> hnls_fd(newton_object_fd);
774 Eigen::HybridNonLinearSolver<newton_usmsr_n> hnls_n(newton_object_n);
775 Eigen::HybridNonLinearSolver<newton_usmsr_t> hnls_t(newton_object_t);
776 Color::Modifier red(Color::FG_RED);
777 Color::Modifier green(Color::FG_GREEN);
778 Color::Modifier def(Color::FG_DEFAULT);
779
780 while (time < finalTime + dt)
781 {
782 time = time + dt;
783 std::vector<double> tv;
784 tv.resize(mu_online.size() + 1);
785 tv[0] = time;
786
787 for (int k = 1; k < tv.size(); k++)
788 {
789 tv[k] = mu_online(k - 1);
790 }
791
792 for (int i = 0; i < Nphi_const; i++)
793 {
794 newton_object_n.d_c(i) = problem->rbfsplines_D[i]->eval(tv);
795 newton_object_n.nsf_c(i) = problem->rbfsplines_NSF[i]->eval(tv);
796 newton_object_n.a_c(i) = problem->rbfsplines_A[i]->eval(tv);
797 newton_object_t.v_c(i) = problem->rbfsplines_v[i]->eval(tv);
798 newton_object_t.sp_c(i) = problem->rbfsplines_SP[i]->eval(tv);
799 newton_object_t.txs_c(i) = problem->rbfsplines_TXS[i]->eval(tv);
800 }
801
802 Eigen::VectorXd res_fd(y);
803 Eigen::VectorXd res_n(w);
804 Eigen::VectorXd res_t(z);
805 res_fd.setZero();
806 res_n.setZero();
807 res_t.setZero();
808 hnls_fd.solve(y);
809
810 for (int j = 0; j < N_BC; j++)
811 {
812 y(j) = vel_now(j, 0);
813 }
814
815 newton_object_n.a_tmp = y.head(Nphi_u);
816 hnls_n.solve(w);
817 newton_object_t.a_tmp = y.head(Nphi_u);
818 newton_object_t.c_tmp = w.head(Nphi_flux);
819 hnls_t.solve(z);
820
821 for (int j = 0; j < N_BCt; j++)
822 {
823 z(j) = temp_now(j, 0);
824 }
825
826 newton_object_fd.operator()(y, res_fd);
827 newton_object_fd.y_old = y;
828 newton_object_n.operator()(w, res_n);
829 newton_object_n.w_old = w;
830 newton_object_t.operator()(z, res_t);
831 newton_object_t.z_old = z;
832 Info << "################## Online solve N° " << count_online_solve <<
833 " ##################" << endl;
834 Info << "Time = " << time << endl;
835
836 if (res_fd.norm() / y.norm() < 1e-5)
837 {
838 Info << green << "|F_fd(x)| = " << res_fd.norm() / y.norm() <<
839 " - Minimun reached in " << hnls_fd.iter << " iterations " << def << endl
840 << endl;
841 }
842 else
843 {
844 Info << red << "|F_fd(x)| = " << res_fd.norm() / y.norm() <<
845 " - Minimun reached in " << hnls_fd.iter << " iterations " << def << endl
846 << endl;
847 }
848
849 if (res_n.norm() / w.norm() < 1e-5)
850 {
851 Info << green << "|F_n(x)| = " << res_n.norm() / w.norm() <<
852 " - Minimun reached in " << hnls_n.iter << " iterations " << def << endl <<
853 endl;
854 }
855 else
856 {
857 Info << red << "|F_n(x)| = " << res_n.norm() / w.norm() <<
858 " - Minimun reached in " << hnls_n.iter << " iterations " << def << endl <<
859 endl;
860 }
861
862 if (res_t.norm() / z.norm() < 1e-5)
863 {
864 Info << green << "|F_t(x)| = " << res_t.norm() / z.norm() <<
865 " - Minimun reached in " << hnls_t.iter << " iterations " << def << endl <<
866 endl;
867 }
868 else
869 {
870 Info << red << "|F_t(x)| = " << res_t.norm() / z.norm() <<
871 " - Minimun reached in " << hnls_t.iter << " iterations " << def << endl <<
872 endl;
873 }
874
875 count_online_solve += 1;
876 tmp_sol_fd(0) = time;
877 tmp_sol_fd.col(0).tail(y.rows()) = y;
878
879 if (counter >= online_solution_fd.size())
880 {
881 online_solution_fd.append(tmp_sol_fd);
882 }
883 else
884 {
885 online_solution_fd[counter] = tmp_sol_fd;
886 }
887
888 tmp_sol_n(0) = time;
889 tmp_sol_n.col(0).tail(w.rows()) = w;
890
891 if (counter >= online_solution_n.size())
892 {
893 online_solution_n.append(tmp_sol_n);
894 }
895 else
896 {
897 online_solution_n[counter] = tmp_sol_n;
898 }
899
900 tmp_sol_t(0) = time;
901 tmp_sol_t.col(0).tail(z.rows()) = z;
902
903 if (counter >= online_solution_t.size())
904 {
905 online_solution_t.append(tmp_sol_t);
906 }
907 else
908 {
909 online_solution_t[counter] = tmp_sol_t;
910 }
911
912 tmp_sol_C(0) = time;
913 pos_c = 1;
914 tmp_sol_C.col(0).segment(pos_c, Nphi_const) = newton_object_t.v_c;
915 pos_c += Nphi_const;
916 tmp_sol_C.col(0).segment(pos_c, Nphi_const) = newton_object_n.d_c;
917 pos_c += Nphi_const;
918 tmp_sol_C.col(0).segment(pos_c, Nphi_const) = newton_object_n.nsf_c;
919 pos_c += Nphi_const;
920 tmp_sol_C.col(0).segment(pos_c, Nphi_const) = newton_object_n.a_c;
921 pos_c += Nphi_const;
922 tmp_sol_C.col(0).segment(pos_c, Nphi_const) = newton_object_t.sp_c;
923 pos_c += Nphi_const;
924 tmp_sol_C.col(0).segment(pos_c, Nphi_const) = newton_object_t.txs_c;
925
926 if (counter >= online_solution_C.size())
927 {
928 online_solution_C.append(tmp_sol_C);
929 }
930 else
931 {
932 online_solution_C[counter] = tmp_sol_C;
933 }
934
935 counter++;
936 }
937
938 ITHACAstream::exportMatrix(online_solution_fd, "red_coeff_fd", "matlab",
939 "./ITHACAoutput/red_coeff_fd");
940 ITHACAstream::exportMatrix(online_solution_n, "red_coeff_n", "matlab",
941 "./ITHACAoutput/red_coeff_n");
942 ITHACAstream::exportMatrix(online_solution_t, "red_coeff_t", "matlab",
943 "./ITHACAoutput/red_coeff_t");
944 ITHACAstream::exportMatrix(online_solution_C, "red_coeff_C", "matlab",
945 "./ITHACAoutput/red_coeff_C");
946}
947
948void reducedusMSR::reconstructAP(fileName folder, int printevery)
949{
950 recall = true;
951 mkDir(folder);
952 ITHACAutilities::createSymLink(folder);
953 reconstruct_fd(folder, printevery);
954 reconstruct_n(folder, printevery);
955 reconstruct_C(folder, printevery);
956 reconstruct_t(folder, printevery);
957}
958
959void reducedusMSR::reconstruct_fd(fileName folder, int printevery)
960{
961 if (recall == false)
962 {
963 mkDir(folder);
964 ITHACAutilities::createSymLink(folder);
965 }
966
967 Info << "Reconstructing online solution | fluid-dynamics" << endl;
968 int counter = 0;
969 int nextwrite = 0;
970 int counter2 = 1;
971
972 for (int i = 0; i < online_solution_fd.size(); i++)
973 {
974 if (counter == nextwrite)
975 {
976 volVectorField U_rec("U", Umodes[0] * 0);
977
978 for (int j = 0; j < Nphi_u; j++)
979 {
980 U_rec += Umodes[j] * online_solution_fd[i](j + 1, 0);
981 }
982
983 ITHACAstream::exportSolution(U_rec, name(counter2), folder);
984 volScalarField P_rec("p", Pmodes[0] * 0);
985
986 for (int j = 0; j < Nphi_p; j++)
987 {
988 P_rec += Pmodes[j] * online_solution_fd[i](j + Nphi_u + 1, 0);
989 }
990
991 ITHACAstream::exportSolution(P_rec, name(counter2), folder);
992 std::ofstream of(folder + "/" + name(counter2) + "/" + name(
993 online_solution_fd[i](0)));
994 nextwrite += printevery;
995 counter2 ++;
996 UREC.append(U_rec.clone());
997 PREC.append(P_rec.clone());
998 }
999
1000 counter++;
1001 }
1002
1003 Info << "End" << endl;
1004}
1005
1006void reducedusMSR::reconstruct_n(fileName folder, int printevery)
1007{
1008 if (recall == false)
1009 {
1010 mkDir(folder);
1011 ITHACAutilities::createSymLink(folder);
1012 }
1013
1014 Info << "Reconstructing online solution | neutronics" << endl;
1015 int counter = 0;
1016 int nextwrite = 0;
1017 int counter2 = 1;
1018
1019 for (int i = 0; i < online_solution_n.size(); i++)
1020 {
1021 if (counter == nextwrite)
1022 {
1023 volScalarField Flux_rec("flux", Fluxmodes[0] * 0);
1024
1025 for (int j = 0; j < Nphi_flux; j++)
1026 {
1027 Flux_rec += Fluxmodes[j] * online_solution_n[i](j + 1, 0);
1028 }
1029
1030 ITHACAstream::exportSolution(Flux_rec, name(counter2), folder);
1031 int pos = Nphi_flux;
1032 volScalarField Prec1_rec("prec1", Prec1modes[0] * 0);
1033
1034 for (int j = 0; j < Nphi_prec1; j++)
1035 {
1036 Prec1_rec += Prec1modes[j] * online_solution_n[i](j + pos + 1, 0);
1037 }
1038
1039 ITHACAstream::exportSolution(Prec1_rec, name(counter2), folder);
1040 pos += Nphi_prec1;
1041 volScalarField Prec2_rec("prec2", Prec2modes[0] * 0);
1042
1043 for (int j = 0; j < Nphi_prec2; j++)
1044 {
1045 Prec2_rec += Prec2modes[j] * online_solution_n[i](j + pos + 1, 0);
1046 }
1047
1048 ITHACAstream::exportSolution(Prec2_rec, name(counter2), folder);
1049 pos += Nphi_prec2;
1050 volScalarField Prec3_rec("prec3", Prec3modes[0] * 0);
1051
1052 for (int j = 0; j < Nphi_prec3; j++)
1053 {
1054 Prec3_rec += Prec3modes[j] * online_solution_n[i](j + pos + 1, 0);
1055 }
1056
1057 ITHACAstream::exportSolution(Prec3_rec, name(counter2), folder);
1058 pos += Nphi_prec3;
1059 volScalarField Prec4_rec("prec4", Prec4modes[0] * 0);
1060
1061 for (int j = 0; j < Nphi_prec4; j++)
1062 {
1063 Prec4_rec += Prec4modes[j] * online_solution_n[i](j + pos + 1, 0);
1064 }
1065
1066 ITHACAstream::exportSolution(Prec4_rec, name(counter2), folder);
1067 pos += Nphi_prec4;
1068 volScalarField Prec5_rec("prec5", Prec5modes[0] * 0);
1069
1070 for (int j = 0; j < Nphi_prec5; j++)
1071 {
1072 Prec5_rec += Prec5modes[j] * online_solution_n[i](j + pos + 1, 0);
1073 }
1074
1075 ITHACAstream::exportSolution(Prec5_rec, name(counter2), folder);
1076 pos += Nphi_prec5;
1077 volScalarField Prec6_rec("prec6", Prec6modes[0] * 0);
1078
1079 for (int j = 0; j < Nphi_prec6; j++)
1080 {
1081 Prec6_rec += Prec6modes[j] * online_solution_n[i](j + pos + 1, 0);
1082 }
1083
1084 ITHACAstream::exportSolution(Prec6_rec, name(counter2), folder);
1085 pos += Nphi_prec6;
1086 volScalarField Prec7_rec("prec7", Prec7modes[0] * 0);
1087
1088 for (int j = 0; j < Nphi_prec7; j++)
1089 {
1090 Prec7_rec += Prec7modes[j] * online_solution_n[i](j + pos + 1, 0);
1091 }
1092
1093 ITHACAstream::exportSolution(Prec7_rec, name(counter2), folder);
1094 pos += Nphi_prec7;
1095 volScalarField Prec8_rec("prec8", Prec8modes[0] * 0);
1096
1097 for (int j = 0; j < Nphi_prec8; j++)
1098 {
1099 Prec8_rec += Prec8modes[j] * online_solution_n[i](j + pos + 1, 0);
1100 }
1101
1102 ITHACAstream::exportSolution(Prec8_rec, name(counter2), folder);
1103 std::ofstream of(folder + "/" + name(counter2) + "/" + name(
1104 online_solution_n[i](0)));
1105 nextwrite += printevery;
1106 counter2 ++;
1107 FLUXREC.append(Flux_rec.clone());
1108 PREC1REC.append(Prec1_rec.clone());
1109 PREC2REC.append(Prec2_rec.clone());
1110 PREC3REC.append(Prec3_rec.clone());
1111 PREC4REC.append(Prec4_rec.clone());
1112 PREC5REC.append(Prec5_rec.clone());
1113 PREC6REC.append(Prec6_rec.clone());
1114 PREC7REC.append(Prec7_rec.clone());
1115 PREC8REC.append(Prec8_rec.clone());
1116 }
1117
1118 counter++;
1119 }
1120
1121 Info << "End" << endl;
1122}
1123
1124void reducedusMSR::reconstruct_t(fileName folder, int printevery)
1125{
1126 if (recall == false)
1127 {
1128 mkDir(folder);
1129 ITHACAutilities::createSymLink(folder);
1130 }
1131
1132 Info << "Reconstructing online solution | thermal" << endl;
1133 int counter = 0;
1134 int nextwrite = 0;
1135 int counter2 = 1;
1136 dimensionedScalar decLam1("decLam1", dimensionSet(0, 0, -1, 0, 0, 0, 0), dl1);
1137 dimensionedScalar decLam2("decLam2", dimensionSet(0, 0, -1, 0, 0, 0, 0), dl2);
1138 dimensionedScalar decLam3("decLam3", dimensionSet(0, 0, -1, 0, 0, 0, 0), dl3);
1139
1140 for (int i = 0; i < online_solution_t.size(); i++)
1141 {
1142 if (counter == nextwrite)
1143 {
1144 volScalarField T_rec("T", Tmodes[0] * 0);
1145
1146 for (int j = 0; j < Nphi_T; j++)
1147 {
1148 T_rec += Tmodes[j] * online_solution_t[i](j + 1, 0);
1149 }
1150
1151 ITHACAstream::exportSolution(T_rec, name(counter2), folder);
1152 int pos = Nphi_T;
1153 volScalarField PowerDens_rec("powerDens", Dec1modes[0] * 0 * decLam1);
1154 volScalarField Dec1_rec("dec1", Dec1modes[0] * 0);
1155
1156 for (int j = 0; j < Nphi_dec1; j++)
1157 {
1158 Dec1_rec += Dec1modes[j] * online_solution_t[i](j + pos + 1, 0);
1159 PowerDens_rec += Dec1modes[j] * online_solution_t[i](j + pos + 1, 0) * decLam1;
1160 }
1161
1162 ITHACAstream::exportSolution(Dec1_rec, name(counter2), folder);
1163 pos += Nphi_dec1;
1164 volScalarField Dec2_rec("dec2", Dec2modes[0] * 0);
1165
1166 for (int j = 0; j < Nphi_dec2; j++)
1167 {
1168 Dec2_rec += Dec2modes[j] * online_solution_t[i](j + pos + 1, 0);
1169 PowerDens_rec += Dec2modes[j] * online_solution_t[i](j + pos + 1, 0) * decLam2;
1170 }
1171
1172 ITHACAstream::exportSolution(Dec2_rec, name(counter2), folder);
1173 pos += Nphi_dec2;
1174 volScalarField Dec3_rec("dec3", Dec3modes[0] * 0);
1175
1176 for (int j = 0; j < Nphi_dec3; j++)
1177 {
1178 Dec3_rec += Dec3modes[j] * online_solution_t[i](j + pos + 1, 0);
1179 PowerDens_rec += Dec3modes[j] * online_solution_t[i](j + pos + 1, 0) * decLam3;
1180 }
1181
1182 ITHACAstream::exportSolution(Dec3_rec, name(counter2), folder);
1183 PowerDens_rec += (1 - dbtot) * SPREC[counter2 - 1] * FLUXREC[counter2 - 1];
1184 ITHACAstream::exportSolution(PowerDens_rec, name(counter2), folder);
1185 std::ofstream of(folder + "/" + name(counter2) + "/" + name(
1186 online_solution_t[i](0)));
1187 nextwrite += printevery;
1188 counter2 ++;
1189 TREC.append(T_rec.clone());
1190 DEC1REC.append(Dec1_rec.clone());
1191 DEC2REC.append(Dec2_rec.clone());
1192 DEC3REC.append(Dec3_rec.clone());
1193 POWERDENSREC.append(PowerDens_rec.clone());
1194 }
1195
1196 counter++;
1197 }
1198
1199 Info << "End" << endl;
1200}
1201
1202void reducedusMSR::reconstruct_C(fileName folder, int printevery)
1203{
1204 if (recall == false)
1205 {
1206 mkDir(folder);
1207 ITHACAutilities::createSymLink(folder);
1208 }
1209
1210 Info << "Reconstructing temperature changing constants" << endl;
1211 int counter = 0;
1212 int nextwrite = 0;
1213 int counter2 = 1;
1214
1215 for (int i = 0; i < online_solution_C.size(); i++)
1216 {
1217 if (counter == nextwrite)
1218 {
1219 volScalarField v_rec("v", vmodes[0] * 0);
1220
1221 for (int j = 0; j < Nphi_const; j++)
1222 {
1223 v_rec += vmodes[j] * online_solution_C[i](j + 1, 0);
1224 }
1225
1226 ITHACAstream::exportSolution(v_rec, name(counter2), folder);
1227 int pos = Nphi_const;
1228 volScalarField D_rec("D", Dmodes[0] * 0);
1229
1230 for (int j = 0; j < Nphi_const; j++)
1231 {
1232 D_rec += Dmodes[j] * online_solution_C[i](j + pos + 1, 0);
1233 }
1234
1235 ITHACAstream::exportSolution(D_rec, name(counter2), folder);
1236 pos += Nphi_const;
1237 volScalarField NSF_rec("NSF", NSFmodes[0] * 0);
1238
1239 for (int j = 0; j < Nphi_const; j++)
1240 {
1241 NSF_rec += NSFmodes[j] * online_solution_C[i](j + pos + 1, 0);
1242 }
1243
1244 ITHACAstream::exportSolution(NSF_rec, name(counter2), folder);
1245 pos += Nphi_const;
1246 volScalarField A_rec("A", Amodes[0] * 0);
1247
1248 for (int j = 0; j < Nphi_const; j++)
1249 {
1250 A_rec += Amodes[j] * online_solution_C[i](j + pos + 1, 0);
1251 }
1252
1253 ITHACAstream::exportSolution(A_rec, name(counter2), folder);
1254 pos += Nphi_const;
1255 volScalarField SP_rec("SP", SPmodes[0] * 0);
1256
1257 for (int j = 0; j < Nphi_const; j++)
1258 {
1259 SP_rec += SPmodes[j] * online_solution_C[i](j + pos + 1, 0);
1260 }
1261
1262 ITHACAstream::exportSolution(SP_rec, name(counter2), folder);
1263 pos += Nphi_const;
1264 volScalarField TXS_rec("TXS", TXSmodes[0] * 0);
1265
1266 for (int j = 0; j < Nphi_const; j++)
1267 {
1268 TXS_rec += TXSmodes[j] * online_solution_C[i](j + pos + 1, 0);
1269 }
1270
1271 ITHACAstream::exportSolution(TXS_rec, name(counter2), folder);
1272 std::ofstream of(folder + "/" + name(counter2) + "/" + name(
1273 online_solution_C[i](0)));
1274 nextwrite += printevery;
1275 counter2 ++;
1276 vREC.append(v_rec.clone());
1277 DREC.append(D_rec.clone());
1278 NSFREC.append(NSF_rec.clone());
1279 AREC.append(A_rec.clone());
1280 SPREC.append(SP_rec.clone());
1281 TXSREC.append(TXS_rec.clone());
1282 }
1283
1284 counter++;
1285 }
1286
1287 Info << "End" << endl;
1288}
1289
reducedMSR::Usnapshots
PtrList< volVectorField > Usnapshots
List of pointers to store the snapshots for each field.
Definition ReducedMSR.H:289
reducedMSR::Umodes
PtrList< volVectorField > Umodes
List of pointers to store the modes for each field.
Definition ReducedMSR.H:265
reducedMSR::nu
scalar nu
characteristic constants of the problem
Definition ReducedMSR.H:217
reducedMSR::UREC
PtrList< volVectorField > UREC
Recontructed fields.
Definition ReducedMSR.H:312
reducedMSR::loadConstants
void loadConstants(msrProblem *problem)
Method to load all the constants needed in the ROM from ///the FOM.
Definition ReducedMSR.C:1032
reducedMSR::N_BC
int N_BC
Number of parametrized boundary conditions.
Definition ReducedMSR.H:366
reducedMSR::online_solution_fd
List< Eigen::MatrixXd > online_solution_fd
List of Eigen matrices to store the online solution.
Definition ReducedMSR.H:259
reducedMSR::count_online_solve
int count_online_solve
Counter to count the online solutions.
Definition ReducedMSR.H:370
reducedMSR::Nphi_u
int Nphi_u
Number of modes for each field.
Definition ReducedMSR.H:344
reducedMSR::y
Eigen::VectorXd y
Vector to store the solution during the Newton procedure.
Definition ReducedMSR.H:210
reducedProblem::solveOnline
virtual void solveOnline()
Virtual Method to perform and online Solve.
Definition ReducedProblem.C:56
ITHACAstream::exportSolution
void exportSolution(GeometricField< Type, PatchField, GeoMesh > &s, fileName subfolder, fileName folder, word fieldName)
Export a field to file in a certain folder and subfolder.
Definition ITHACAstream.C:992
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:349
ITHACAutilities::createSymLink
void createSymLink(word folder)
Creates symbolic links to 0, system and constant.
Definition ITHACAsystem.C:30
newton_usmsr_fd::operator()
int operator()(const Eigen::VectorXd &x, Eigen::VectorXd &fvec) const
Definition ReducedUnsteadyMSR.C:185
newton_usmsr_n::operator()
int operator()(const Eigen::VectorXd &n, Eigen::VectorXd &fvecn) const
Definition ReducedUnsteadyMSR.C:243
newton_usmsr_t::operator()
int operator()(const Eigen::VectorXd &t, Eigen::VectorXd &fvect) const
Definition ReducedUnsteadyMSR.C:463
Generated by   doxygen 1.16.1