Line data Source code
1 : // rbOOmit: An implementation of the Certified Reduced Basis method.
2 : // Copyright (C) 2009, 2010 David J. Knezevic
3 :
4 : // This file is part of rbOOmit.
5 :
6 : // rbOOmit is free software; you can redistribute it and/or
7 : // modify it under the terms of the GNU Lesser General Public
8 : // License as published by the Free Software Foundation; either
9 : // version 2.1 of the License, or (at your option) any later version.
10 :
11 : // rbOOmit is distributed in the hope that it will be useful,
12 : // but WITHOUT ANY WARRANTY; without even the implied warranty of
13 : // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
14 : // Lesser General Public License for more details.
15 :
16 : // You should have received a copy of the GNU Lesser General Public
17 : // License along with this library; if not, write to the Free Software
18 : // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
19 :
20 : // rbOOmit includes
21 : #include "libmesh/transient_rb_evaluation.h"
22 : #include "libmesh/transient_rb_theta_expansion.h"
23 :
24 : // libMesh includes
25 : #include "libmesh/getpot.h"
26 : #include "libmesh/int_range.h"
27 : #include "libmesh/libmesh_logging.h"
28 : #include "libmesh/numeric_vector.h"
29 : #include "libmesh/xdr_cxx.h"
30 :
31 : // C++ includes
32 : #include <fstream>
33 : #include <sstream>
34 : #include <iomanip>
35 :
36 : namespace libMesh
37 : {
38 :
39 12 : TransientRBEvaluation::TransientRBEvaluation(const Parallel::Communicator & comm_in) :
40 : RBEvaluation(comm_in),
41 28 : _rb_solve_data_cached(false)
42 : {
43 : // Indicate that we need to compute the RB
44 : // inner product matrix in this case
45 12 : compute_RB_inner_product = true;
46 12 : }
47 :
48 28 : TransientRBEvaluation::~TransientRBEvaluation () = default;
49 :
50 0 : void TransientRBEvaluation::clear()
51 : {
52 0 : Parent::clear();
53 :
54 0 : clear_riesz_representors();
55 0 : }
56 :
57 12 : void TransientRBEvaluation::clear_riesz_representors()
58 : {
59 12 : Parent::clear_riesz_representors();
60 :
61 : // Delete the M_q representors
62 8 : M_q_representor.clear();
63 12 : }
64 :
65 12 : void TransientRBEvaluation::resize_data_structures(const unsigned int Nmax,
66 : bool resize_error_bound_data)
67 : {
68 8 : LOG_SCOPE("resize_data_structures()", "TransientRBEvaluation");
69 :
70 12 : Parent::resize_data_structures(Nmax, resize_error_bound_data);
71 :
72 8 : RB_L2_matrix.resize(Nmax,Nmax);
73 8 : RB_LHS_matrix.resize(Nmax,Nmax);
74 8 : RB_RHS_matrix.resize(Nmax,Nmax);
75 8 : RB_RHS_save.resize(Nmax);
76 :
77 : TransientRBThetaExpansion & trans_theta_expansion =
78 12 : cast_ref<TransientRBThetaExpansion &>(get_rb_theta_expansion());
79 12 : const unsigned int Q_m = trans_theta_expansion.get_n_M_terms();
80 12 : const unsigned int Q_a = trans_theta_expansion.get_n_A_terms();
81 12 : const unsigned int Q_f = trans_theta_expansion.get_n_F_terms();
82 :
83 : // Allocate dense matrices for RB solves
84 12 : RB_M_q_vector.resize(Q_m);
85 24 : for (unsigned int q=0; q<Q_m; q++)
86 : {
87 : // Initialize the memory for the RB matrices
88 12 : RB_M_q_vector[q].resize(Nmax,Nmax);
89 : }
90 :
91 : // Initialize the initial condition storage
92 12 : RB_initial_condition_all_N.resize(Nmax);
93 252 : for (auto i : make_range(RB_initial_condition_all_N.size()))
94 : {
95 : // The i^th row holds a vector of length i+1
96 320 : RB_initial_condition_all_N[i].resize(i+1);
97 : }
98 :
99 12 : initial_L2_error_all_N.resize(Nmax, 0.);
100 :
101 :
102 12 : if (resize_error_bound_data)
103 : {
104 : // Initialize vectors for the norms of the representors
105 12 : Fq_Mq_representor_innerprods.resize(Q_f);
106 24 : for (unsigned int i=0; i<Q_f; i++)
107 : {
108 16 : Fq_Mq_representor_innerprods[i].resize(Q_m);
109 24 : for (unsigned int j=0; j<Q_m; j++)
110 : {
111 20 : Fq_Mq_representor_innerprods[i][j].resize(Nmax, 0.);
112 : }
113 : }
114 :
115 12 : unsigned int Q_m_hat = Q_m*(Q_m+1)/2;
116 12 : Mq_Mq_representor_innerprods.resize(Q_m_hat);
117 24 : for (unsigned int i=0; i<Q_m_hat; i++)
118 : {
119 16 : Mq_Mq_representor_innerprods[i].resize(Nmax);
120 252 : for (unsigned int j=0; j<Nmax; j++)
121 : {
122 400 : Mq_Mq_representor_innerprods[i][j].resize(Nmax, 0.);
123 : }
124 : }
125 :
126 12 : Aq_Mq_representor_innerprods.resize(Q_a);
127 48 : for (unsigned int i=0; i<Q_a; i++)
128 : {
129 48 : Aq_Mq_representor_innerprods[i].resize(Q_m);
130 72 : for (unsigned int j=0; j<Q_m; j++)
131 : {
132 60 : Aq_Mq_representor_innerprods[i][j].resize(Nmax);
133 756 : for (unsigned int k=0; k<Nmax; k++)
134 : {
135 1440 : Aq_Mq_representor_innerprods[i][j][k].resize(Nmax, 0.);
136 : }
137 : }
138 : }
139 :
140 : // Resize M_q_representor
141 : // This is cleared in the call to clear_riesz_representors
142 : // in Parent::resize_RB_data, so just resize here
143 12 : M_q_representor.resize(Q_m);
144 24 : for (unsigned int q_m=0; q_m<Q_m; q_m++)
145 : {
146 16 : M_q_representor[q_m].resize(Nmax);
147 : }
148 : }
149 12 : }
150 :
151 6006 : Real TransientRBEvaluation::rb_solve(unsigned int N)
152 : {
153 4004 : LOG_SCOPE("rb_solve()", "TransientRBEvaluation");
154 :
155 6006 : libmesh_error_msg_if(N > get_n_basis_functions(),
156 : "ERROR: N cannot be larger than the number of basis functions in rb_solve");
157 :
158 6006 : const RBParameters & mu = get_parameters();
159 :
160 : TransientRBThetaExpansion & trans_theta_expansion =
161 6006 : cast_ref<TransientRBThetaExpansion &>(get_rb_theta_expansion());
162 6006 : const unsigned int Q_m = trans_theta_expansion.get_n_M_terms();
163 6006 : const unsigned int Q_a = trans_theta_expansion.get_n_A_terms();
164 6006 : const unsigned int Q_f = trans_theta_expansion.get_n_F_terms();
165 :
166 6006 : const unsigned int n_time_steps = get_n_time_steps();
167 6006 : const Real dt = get_delta_t();
168 6006 : const Real euler_theta = get_euler_theta();
169 :
170 : // Resize the RB and error bound vectors
171 6006 : error_bound_all_k.resize(n_time_steps+1);
172 6006 : RB_outputs_all_k.resize(trans_theta_expansion.get_n_outputs());
173 6006 : RB_output_error_bounds_all_k.resize(trans_theta_expansion.get_n_outputs());
174 30030 : for (unsigned int n=0; n<trans_theta_expansion.get_n_outputs(); n++)
175 : {
176 32032 : RB_outputs_all_k[n].resize(n_time_steps+1, 0.);
177 32032 : RB_output_error_bounds_all_k[n].resize(n_time_steps+1, 0.);
178 : }
179 :
180 : // First assemble the mass matrix
181 10010 : DenseMatrix<Number> RB_mass_matrix_N(N,N);
182 2002 : RB_mass_matrix_N.zero();
183 10010 : DenseMatrix<Number> RB_M_q_m;
184 12012 : for (unsigned int q_m=0; q_m<Q_m; q_m++)
185 : {
186 8008 : RB_M_q_vector[q_m].get_principal_submatrix(N, RB_M_q_m);
187 6006 : RB_mass_matrix_N.add(trans_theta_expansion.eval_M_theta(q_m, mu), RB_M_q_m);
188 : }
189 :
190 6006 : RB_LHS_matrix.resize(N,N);
191 2002 : RB_LHS_matrix.zero();
192 :
193 6006 : RB_RHS_matrix.resize(N,N);
194 2002 : RB_RHS_matrix.zero();
195 :
196 6006 : RB_LHS_matrix.add(1./dt, RB_mass_matrix_N);
197 6006 : RB_RHS_matrix.add(1./dt, RB_mass_matrix_N);
198 :
199 10010 : DenseMatrix<Number> RB_Aq_a;
200 24024 : for (unsigned int q_a=0; q_a<Q_a; q_a++)
201 : {
202 24024 : RB_Aq_vector[q_a].get_principal_submatrix(N, RB_Aq_a);
203 :
204 18018 : RB_LHS_matrix.add( euler_theta*trans_theta_expansion.eval_A_theta(q_a,mu), RB_Aq_a);
205 18018 : RB_RHS_matrix.add( -(1.-euler_theta)*trans_theta_expansion.eval_A_theta(q_a,mu), RB_Aq_a);
206 : }
207 :
208 : // Add forcing terms
209 6006 : DenseVector<Number> RB_Fq_f;
210 4004 : RB_RHS_save.resize(N);
211 2002 : RB_RHS_save.zero();
212 12012 : for (unsigned int q_f=0; q_f<Q_f; q_f++)
213 : {
214 8008 : RB_Fq_vector[q_f].get_principal_subvector(N, RB_Fq_f);
215 6006 : RB_RHS_save.add(trans_theta_expansion.eval_F_theta(q_f,mu), RB_Fq_f);
216 : }
217 :
218 : // Set system time level to 0
219 6006 : set_time_step(0);
220 :
221 : // Resize/clear the solution vector
222 4004 : RB_solution.resize(N);
223 :
224 : // Load the initial condition into RB_solution
225 6006 : if (N > 0)
226 : {
227 5706 : RB_solution = RB_initial_condition_all_N[N-1];
228 : }
229 :
230 : // Resize/clear the old solution vector
231 4004 : old_RB_solution.resize(N);
232 :
233 : // Initialize the RB rhs
234 6006 : DenseVector<Number> RB_rhs(N);
235 2002 : RB_rhs.zero();
236 :
237 : // Initialize the vectors storing solution data
238 6006 : RB_temporal_solution_data.resize(n_time_steps+1);
239 612612 : for (unsigned int time_level=0; time_level<=n_time_steps; time_level++)
240 : {
241 606606 : RB_temporal_solution_data[time_level].resize(N);
242 : }
243 : // and load the initial data
244 4004 : RB_temporal_solution_data[0] = RB_solution;
245 :
246 : // Set outputs at initial time
247 : {
248 6006 : DenseVector<Number> RB_output_vector_N;
249 30030 : for (unsigned int n=0; n<trans_theta_expansion.get_n_outputs(); n++)
250 : {
251 32032 : RB_outputs_all_k[n][0] = 0.;
252 48048 : for (unsigned int q_l=0; q_l<trans_theta_expansion.get_n_output_terms(n); q_l++)
253 : {
254 40040 : RB_output_vectors[n][q_l].get_principal_subvector(N, RB_output_vector_N);
255 24024 : RB_outputs_all_k[n][0] += trans_theta_expansion.eval_output_theta(n,q_l,mu)*RB_output_vector_N.dot(RB_solution);
256 : }
257 : }
258 : }
259 :
260 : // Initialize error bounds, if necessary
261 2002 : Real error_bound_sum = 0.;
262 2002 : Real alpha_LB = 0.;
263 6006 : if (evaluate_RB_error_bound)
264 : {
265 6006 : if (N > 0)
266 : {
267 7608 : error_bound_sum += pow( initial_L2_error_all_N[N-1], 2.);
268 : }
269 :
270 : // Set error bound at the initial time
271 6006 : error_bound_all_k[get_time_step()] = std::sqrt(error_bound_sum);
272 :
273 : // Compute the outputs and associated error bounds at the initial time
274 6006 : DenseVector<Number> RB_output_vector_N;
275 30030 : for (unsigned int n=0; n<trans_theta_expansion.get_n_outputs(); n++)
276 : {
277 32032 : RB_outputs_all_k[n][0] = 0.;
278 48048 : for (unsigned int q_l=0; q_l<trans_theta_expansion.get_n_output_terms(n); q_l++)
279 : {
280 40040 : RB_output_vectors[n][q_l].get_principal_subvector(N, RB_output_vector_N);
281 24024 : RB_outputs_all_k[n][0] += trans_theta_expansion.eval_output_theta(n,q_l,mu)*RB_output_vector_N.dot(RB_solution);
282 : }
283 :
284 24024 : RB_output_error_bounds_all_k[n][0] = error_bound_all_k[0] * eval_output_dual_norm(n, nullptr);
285 : }
286 :
287 6006 : alpha_LB = get_stability_lower_bound();
288 :
289 : // Precompute time-invariant parts of the dual norm of the residual.
290 6006 : cache_online_residual_terms(N);
291 : }
292 :
293 606606 : for (unsigned int time_level=1; time_level<=n_time_steps; time_level++)
294 : {
295 600600 : set_time_step(time_level);
296 200200 : old_RB_solution = RB_solution;
297 :
298 : // Compute RB_rhs, as RB_LHS_matrix x old_RB_solution
299 600600 : RB_RHS_matrix.vector_mult(RB_rhs, old_RB_solution);
300 :
301 : // Add forcing terms
302 600600 : RB_rhs.add(get_control(time_level), RB_RHS_save);
303 :
304 600600 : if (N > 0)
305 : {
306 570600 : RB_LHS_matrix.lu_solve(RB_rhs, RB_solution);
307 : }
308 :
309 : // Save RB_solution for current time level
310 600600 : RB_temporal_solution_data[time_level] = RB_solution;
311 :
312 : // Evaluate outputs
313 600600 : DenseVector<Number> RB_output_vector_N;
314 3003000 : for (unsigned int n=0; n<trans_theta_expansion.get_n_outputs(); n++)
315 : {
316 4004000 : RB_outputs_all_k[n][time_level] = 0.;
317 4804800 : for (unsigned int q_l=0; q_l<trans_theta_expansion.get_n_output_terms(n); q_l++)
318 : {
319 4004000 : RB_output_vectors[n][q_l].get_principal_subvector(N, RB_output_vector_N);
320 6406400 : RB_outputs_all_k[n][time_level] += trans_theta_expansion.eval_output_theta(n,q_l,mu)*
321 2402400 : RB_output_vector_N.dot(RB_solution);
322 : }
323 : }
324 :
325 : // Calculate RB error bounds
326 600600 : if (evaluate_RB_error_bound)
327 : {
328 : // Evaluate the dual norm of the residual for RB_solution_vector
329 : // Real epsilon_N = uncached_compute_residual_dual_norm(N);
330 600600 : Real epsilon_N = compute_residual_dual_norm(N);
331 :
332 600600 : error_bound_sum += residual_scaling_numer(alpha_LB) * pow(epsilon_N, 2.);
333 :
334 : // store error bound at time-level _k
335 600600 : error_bound_all_k[time_level] = std::sqrt(error_bound_sum/residual_scaling_denom(alpha_LB));
336 :
337 : // Now evaluated output error bounds
338 3003000 : for (unsigned int n=0; n<trans_theta_expansion.get_n_outputs(); n++)
339 : {
340 3203200 : RB_output_error_bounds_all_k[n][time_level] = error_bound_all_k[time_level] *
341 2402400 : eval_output_dual_norm(n, nullptr);
342 : }
343 : }
344 : }
345 :
346 6006 : _rb_solve_data_cached = true ;
347 :
348 6006 : if (evaluate_RB_error_bound) // Calculate the error bounds
349 : {
350 8008 : return error_bound_all_k[n_time_steps];
351 : }
352 : else // Don't calculate the error bounds
353 : {
354 : // Just return -1. if we did not compute the error bound
355 0 : return -1.;
356 : }
357 2002 : }
358 :
359 0 : Real TransientRBEvaluation::rb_solve_again()
360 : {
361 0 : libmesh_assert(_rb_solve_data_cached);
362 :
363 0 : const unsigned int n_time_steps = get_n_time_steps();
364 : // Set system time level to 0
365 0 : set_time_step(0);
366 :
367 : // Resize/clear the solution vector
368 0 : const unsigned int N = RB_RHS_save.size();
369 0 : RB_solution.resize(N);
370 :
371 : // Load the initial condition into RB_solution
372 0 : if (N > 0)
373 0 : RB_solution = RB_initial_condition_all_N[N-1];
374 :
375 : // Resize/clear the old solution vector
376 0 : old_RB_solution.resize(N);
377 :
378 : // Initialize the RB rhs
379 0 : DenseVector<Number> RB_rhs(N);
380 0 : RB_rhs.zero();
381 :
382 0 : for (unsigned int time_level=1; time_level<=n_time_steps; time_level++)
383 : {
384 0 : set_time_step(time_level);
385 0 : old_RB_solution = RB_solution;
386 :
387 : // Compute RB_rhs, as *RB_lhs_matrix x old_RB_solution
388 0 : RB_RHS_matrix.vector_mult(RB_rhs, old_RB_solution);
389 :
390 : // Add forcing terms
391 0 : RB_rhs.add(get_control(time_level), RB_RHS_save);
392 :
393 0 : if (N > 0)
394 0 : RB_LHS_matrix.lu_solve(RB_rhs, RB_solution);
395 : }
396 :
397 : {
398 : // Just return -1. We did not compute the error bound
399 0 : return -1.;
400 : }
401 : }
402 :
403 0 : Real TransientRBEvaluation::get_error_bound_normalization()
404 : {
405 : // Just set the normalization factor to 1 in this case.
406 : // Users can override this method if specific behavior
407 : // is required.
408 :
409 0 : return 1.;
410 : }
411 :
412 600600 : Real TransientRBEvaluation::residual_scaling_numer(Real)
413 : {
414 600600 : return get_delta_t();
415 : }
416 :
417 6006 : void TransientRBEvaluation::cache_online_residual_terms(const unsigned int N)
418 : {
419 4004 : LOG_SCOPE("cache_online_residual_terms()", "TransientRBEvaluation");
420 :
421 10010 : const RBParameters mu = get_parameters();
422 :
423 : TransientRBThetaExpansion & trans_theta_expansion =
424 6006 : cast_ref<TransientRBThetaExpansion &>(get_rb_theta_expansion());
425 6006 : const unsigned int Q_m = trans_theta_expansion.get_n_M_terms();
426 6006 : const unsigned int Q_a = trans_theta_expansion.get_n_A_terms();
427 6006 : const unsigned int Q_f = trans_theta_expansion.get_n_F_terms();
428 :
429 6006 : cached_Fq_term = 0.;
430 2002 : unsigned int q=0;
431 12012 : for (unsigned int q_f1=0; q_f1<Q_f; q_f1++)
432 : {
433 6006 : Number cached_theta_q_f1 = trans_theta_expansion.eval_F_theta(q_f1,mu);
434 12012 : for (unsigned int q_f2=q_f1; q_f2<Q_f; q_f2++)
435 : {
436 6006 : Real delta = (q_f1==q_f2) ? 1. : 2.;
437 6006 : cached_Fq_term += delta*cached_theta_q_f1*trans_theta_expansion.eval_F_theta(q_f2,mu) *
438 6006 : Fq_representor_innerprods[q];
439 :
440 6006 : q++;
441 : }
442 : }
443 :
444 4004 : cached_Fq_Aq_vector.resize(N);
445 12012 : for (unsigned int q_f=0; q_f<Q_f; q_f++)
446 : {
447 6006 : Number cached_theta_q_f = trans_theta_expansion.eval_F_theta(q_f,mu);
448 24024 : for (unsigned int q_a=0; q_a<Q_a; q_a++)
449 : {
450 18018 : Number cached_theta_q_a = trans_theta_expansion.eval_A_theta(q_a,mu);
451 189378 : for (unsigned int i=0; i<N; i++)
452 : {
453 228480 : cached_Fq_Aq_vector(i) += 2.*cached_theta_q_f*cached_theta_q_a*
454 285600 : Fq_Aq_representor_innerprods[q_f][q_a][i];
455 : }
456 : }
457 : }
458 :
459 4004 : cached_Aq_Aq_matrix.resize(N,N);
460 2002 : q=0;
461 24024 : for (unsigned int q_a1=0; q_a1<Q_a; q_a1++)
462 : {
463 18018 : Number cached_theta_q_a1 = trans_theta_expansion.eval_A_theta(q_a1,mu);
464 54054 : for (unsigned int q_a2=q_a1; q_a2<Q_a; q_a2++)
465 : {
466 36036 : Number cached_theta_q_a2 = trans_theta_expansion.eval_A_theta(q_a2,mu);
467 36036 : Real delta = (q_a1==q_a2) ? 1. : 2.;
468 :
469 378756 : for (unsigned int i=0; i<N; i++)
470 : {
471 4803120 : for (unsigned int j=0; j<N; j++)
472 : {
473 4460400 : cached_Aq_Aq_matrix(i,j) += delta*
474 4460400 : cached_theta_q_a1*cached_theta_q_a2*
475 7434000 : Aq_Aq_representor_innerprods[q][i][j];
476 : }
477 : }
478 36036 : q++;
479 : }
480 : }
481 :
482 4004 : cached_Fq_Mq_vector.resize(N);
483 12012 : for (unsigned int q_f=0; q_f<Q_f; q_f++)
484 : {
485 6006 : Number cached_theta_q_f = trans_theta_expansion.eval_F_theta(q_f,mu);
486 12012 : for (unsigned int q_m=0; q_m<Q_m; q_m++)
487 : {
488 6006 : Number cached_theta_q_m = trans_theta_expansion.eval_M_theta(q_m,mu);
489 63126 : for (unsigned int i=0; i<N; i++)
490 : {
491 133280 : cached_Fq_Mq_vector(i) += 2.*cached_theta_q_f * cached_theta_q_m * Fq_Mq_representor_innerprods[q_f][q_m][i];
492 : }
493 : }
494 : }
495 :
496 4004 : cached_Aq_Mq_matrix.resize(N,N);
497 24024 : for (unsigned int q_a=0; q_a<Q_a; q_a++)
498 : {
499 18018 : Number cached_theta_q_a = trans_theta_expansion.eval_A_theta(q_a,mu);
500 :
501 36036 : for (unsigned int q_m=0; q_m<Q_m; q_m++)
502 : {
503 18018 : Number cached_theta_q_m = trans_theta_expansion.eval_M_theta(q_m,mu);
504 :
505 189378 : for (unsigned int i=0; i<N; i++)
506 : {
507 2401560 : for (unsigned int j=0; j<N; j++)
508 : {
509 5947200 : cached_Aq_Mq_matrix(i,j) += 2.*cached_theta_q_a*cached_theta_q_m*Aq_Mq_representor_innerprods[q_a][q_m][i][j];
510 : }
511 : }
512 : }
513 : }
514 :
515 4004 : cached_Mq_Mq_matrix.resize(N,N);
516 2002 : q=0;
517 12012 : for (unsigned int q_m1=0; q_m1<Q_m; q_m1++)
518 : {
519 6006 : Number cached_theta_q_m1 = trans_theta_expansion.eval_M_theta(q_m1,mu);
520 12012 : for (unsigned int q_m2=q_m1; q_m2<Q_m; q_m2++)
521 : {
522 6006 : Number cached_theta_q_m2 = trans_theta_expansion.eval_M_theta(q_m2,mu);
523 6006 : Real delta = (q_m1==q_m2) ? 1. : 2.;
524 :
525 63126 : for (unsigned int i=0; i<N; i++)
526 : {
527 800520 : for (unsigned int j=0; j<N; j++)
528 : {
529 743400 : cached_Mq_Mq_matrix(i,j) += delta*
530 743400 : cached_theta_q_m1*cached_theta_q_m2*
531 1239000 : Mq_Mq_representor_innerprods[q][i][j];
532 : }
533 : }
534 6006 : q++;
535 : }
536 : }
537 6006 : }
538 :
539 600600 : Real TransientRBEvaluation::compute_residual_dual_norm(const unsigned int N)
540 : {
541 400400 : LOG_SCOPE("compute_residual_dual_norm()", "TransientRBEvaluation");
542 :
543 : // This assembly assumes we have already called cache_online_residual_terms
544 : // and that the rb_solve parameter is constant in time
545 :
546 600600 : const Real dt = get_delta_t();
547 600600 : const Real euler_theta = get_euler_theta();
548 600600 : const Real current_control = get_control(get_time_step());
549 :
550 600600 : DenseVector<Number> RB_u_euler_theta(N);
551 600600 : DenseVector<Number> mass_coeffs(N);
552 6312600 : for (unsigned int i=0; i<N; i++)
553 : {
554 9520000 : RB_u_euler_theta(i) = euler_theta*RB_solution(i) +
555 7616000 : (1.-euler_theta)*old_RB_solution(i);
556 7616000 : mass_coeffs(i) = -(RB_solution(i) - old_RB_solution(i))/dt;
557 : }
558 :
559 600600 : Number residual_norm_sq = current_control*current_control*cached_Fq_term;
560 :
561 600600 : residual_norm_sq += current_control*RB_u_euler_theta.dot(cached_Fq_Aq_vector);
562 600600 : residual_norm_sq += current_control*mass_coeffs.dot(cached_Fq_Mq_vector);
563 :
564 6312600 : for (unsigned int i=0; i<N; i++)
565 80052000 : for (unsigned int j=0; j<N; j++)
566 : {
567 123900000 : residual_norm_sq += RB_u_euler_theta(i)*RB_u_euler_theta(j)*cached_Aq_Aq_matrix(i,j);
568 99120000 : residual_norm_sq += mass_coeffs(i)*mass_coeffs(j)*cached_Mq_Mq_matrix(i,j);
569 99120000 : residual_norm_sq += RB_u_euler_theta(i)*mass_coeffs(j)*cached_Aq_Mq_matrix(i,j);
570 : }
571 :
572 :
573 600600 : if (libmesh_real(residual_norm_sq) < 0)
574 : {
575 0 : libMesh::out << "Warning: Square of residual norm is negative "
576 0 : << "in TransientRBEvaluation::compute_residual_dual_norm()" << std::endl;
577 :
578 : // Sometimes this is negative due to rounding error,
579 : // but error is on the order of 1.e-10, so shouldn't
580 : // affect result
581 0 : residual_norm_sq = std::abs(residual_norm_sq);
582 : }
583 :
584 1201200 : return libmesh_real(std::sqrt( residual_norm_sq ));
585 : }
586 :
587 0 : Real TransientRBEvaluation::uncached_compute_residual_dual_norm(const unsigned int N)
588 : {
589 0 : LOG_SCOPE("uncached_compute_residual_dual_norm()", "TransientRBEvaluation");
590 :
591 : // Use the stored representor inner product values
592 : // to evaluate the residual norm
593 :
594 0 : const RBParameters & mu = get_parameters();
595 :
596 : TransientRBThetaExpansion & trans_theta_expansion =
597 0 : cast_ref<TransientRBThetaExpansion &>(get_rb_theta_expansion());
598 0 : const unsigned int Q_m = trans_theta_expansion.get_n_M_terms();
599 0 : const unsigned int Q_a = trans_theta_expansion.get_n_A_terms();
600 0 : const unsigned int Q_f = trans_theta_expansion.get_n_F_terms();
601 :
602 0 : const Real dt = get_delta_t();
603 0 : const Real euler_theta = get_euler_theta();
604 :
605 0 : std::vector<Number> RB_u_euler_theta(N);
606 0 : std::vector<Number> mass_coeffs(N);
607 0 : for (unsigned int i=0; i<N; i++)
608 : {
609 0 : RB_u_euler_theta[i] = euler_theta*RB_solution(i) +
610 0 : (1.-euler_theta)*old_RB_solution(i);
611 0 : mass_coeffs[i] = -(RB_solution(i) - old_RB_solution(i))/dt;
612 : }
613 :
614 0 : Number residual_norm_sq = 0.;
615 :
616 0 : unsigned int q=0;
617 0 : for (unsigned int q_f1=0; q_f1<Q_f; q_f1++)
618 : {
619 0 : Number cached_theta_q_f1 = trans_theta_expansion.eval_F_theta(q_f1,mu);
620 0 : for (unsigned int q_f2=q_f1; q_f2<Q_f; q_f2++)
621 : {
622 0 : Real delta = (q_f1==q_f2) ? 1. : 2.;
623 0 : residual_norm_sq += delta*cached_theta_q_f1*trans_theta_expansion.eval_F_theta(q_f2,mu) * Fq_representor_innerprods[q];
624 :
625 0 : q++;
626 : }
627 : }
628 :
629 0 : for (unsigned int q_f=0; q_f<Q_f; q_f++)
630 : {
631 0 : Number cached_theta_q_f = trans_theta_expansion.eval_F_theta(q_f,mu);
632 0 : for (unsigned int q_a=0; q_a<Q_a; q_a++)
633 : {
634 0 : Number cached_theta_q_a = trans_theta_expansion.eval_A_theta(q_a,mu);
635 0 : for (unsigned int i=0; i<N; i++)
636 : {
637 0 : residual_norm_sq += 2.*RB_u_euler_theta[i]*cached_theta_q_f*cached_theta_q_a*
638 0 : Fq_Aq_representor_innerprods[q_f][q_a][i];
639 : }
640 : }
641 : }
642 :
643 0 : q=0;
644 0 : for (unsigned int q_a1=0; q_a1<Q_a; q_a1++)
645 : {
646 0 : Number cached_theta_q_a1 = trans_theta_expansion.eval_A_theta(q_a1,mu);
647 0 : for (unsigned int q_a2=q_a1; q_a2<Q_a; q_a2++)
648 : {
649 0 : Number cached_theta_q_a2 = trans_theta_expansion.eval_A_theta(q_a2,mu);
650 0 : Real delta = (q_a1==q_a2) ? 1. : 2.;
651 :
652 0 : for (unsigned int i=0; i<N; i++)
653 : {
654 0 : for (unsigned int j=0; j<N; j++)
655 : {
656 0 : residual_norm_sq += delta*RB_u_euler_theta[i]*RB_u_euler_theta[j]*
657 0 : cached_theta_q_a1*cached_theta_q_a2*
658 0 : Aq_Aq_representor_innerprods[q][i][j];
659 : }
660 : }
661 0 : q++;
662 : }
663 : }
664 :
665 : // Now add the terms due to the time-derivative
666 0 : q=0;
667 0 : for (unsigned int q_m1=0; q_m1<Q_m; q_m1++)
668 : {
669 0 : Number cached_theta_q_m1 = trans_theta_expansion.eval_M_theta(q_m1,mu);
670 0 : for (unsigned int q_m2=q_m1; q_m2<Q_m; q_m2++)
671 : {
672 0 : Number cached_theta_q_m2 = trans_theta_expansion.eval_M_theta(q_m2,mu);
673 0 : Real delta = (q_m1==q_m2) ? 1. : 2.;
674 :
675 0 : for (unsigned int i=0; i<N; i++)
676 : {
677 0 : for (unsigned int j=0; j<N; j++)
678 : {
679 0 : residual_norm_sq += delta*mass_coeffs[i]*mass_coeffs[j]*
680 0 : cached_theta_q_m1*cached_theta_q_m2*
681 0 : Mq_Mq_representor_innerprods[q][i][j];
682 : }
683 : }
684 0 : q++;
685 : }
686 : }
687 :
688 0 : for (unsigned int q_f=0; q_f<Q_f; q_f++)
689 : {
690 0 : Number cached_theta_q_f = trans_theta_expansion.eval_F_theta(q_f,mu);
691 0 : for (unsigned int q_m=0; q_m<Q_m; q_m++)
692 : {
693 0 : Number cached_theta_q_m = trans_theta_expansion.eval_M_theta(q_m,mu);
694 0 : for (unsigned int i=0; i<N; i++)
695 : {
696 0 : residual_norm_sq += 2.*mass_coeffs[i]*cached_theta_q_f * cached_theta_q_m * Fq_Mq_representor_innerprods[q_f][q_m][i];
697 : }
698 : }
699 : }
700 :
701 0 : for (unsigned int q_a=0; q_a<Q_a; q_a++)
702 : {
703 0 : Number cached_theta_q_a = trans_theta_expansion.eval_A_theta(q_a,mu);
704 :
705 0 : for (unsigned int q_m=0; q_m<Q_m; q_m++)
706 : {
707 0 : Number cached_theta_q_m = trans_theta_expansion.eval_M_theta(q_m,mu);
708 :
709 0 : for (unsigned int i=0; i<N; i++)
710 : {
711 0 : for (unsigned int j=0; j<N; j++)
712 : {
713 0 : residual_norm_sq += 2.*RB_u_euler_theta[i]*mass_coeffs[j]*
714 0 : cached_theta_q_a*cached_theta_q_m*
715 0 : Aq_Mq_representor_innerprods[q_a][q_m][i][j];
716 : }
717 : }
718 : }
719 : }
720 :
721 0 : if (libmesh_real(residual_norm_sq) < 0)
722 : {
723 0 : libMesh::out << "Warning: Square of residual norm is negative "
724 0 : << "in TransientRBEvaluation::compute_residual_dual_norm()" << std::endl;
725 :
726 : // Sometimes this is negative due to rounding error,
727 : // but error is on the order of 1.e-10, so shouldn't
728 : // affect result
729 0 : residual_norm_sq = std::abs(residual_norm_sq);
730 : }
731 :
732 : // libMesh::out << "slow residual_sq = " << slow_residual_norm_sq
733 : // << ", fast residual_sq = " << residual_norm_sq << std::endl;
734 :
735 0 : return libmesh_real(std::sqrt( residual_norm_sq ));
736 : }
737 :
738 6 : void TransientRBEvaluation::legacy_write_offline_data_to_files(const std::string & directory_name,
739 : const bool write_binary_data)
740 : {
741 4 : LOG_SCOPE("legacy_write_offline_data_to_files()", "TransientRBEvaluation");
742 :
743 6 : Parent::legacy_write_offline_data_to_files(directory_name);
744 :
745 : TransientRBThetaExpansion & trans_theta_expansion =
746 6 : cast_ref<TransientRBThetaExpansion &>(get_rb_theta_expansion());
747 6 : const unsigned int Q_m = trans_theta_expansion.get_n_M_terms();
748 6 : const unsigned int Q_a = trans_theta_expansion.get_n_A_terms();
749 6 : const unsigned int Q_f = trans_theta_expansion.get_n_F_terms();
750 :
751 6 : const unsigned int n_bfs = get_n_basis_functions();
752 :
753 : // The writing mode: ENCODE for binary, WRITE for ASCII
754 6 : XdrMODE mode = write_binary_data ? ENCODE : WRITE;
755 :
756 : // The suffix to use for all the files that are written out
757 8 : const std::string suffix = write_binary_data ? ".xdr" : ".dat";
758 :
759 8 : if (this->processor_id() == 0)
760 : {
761 5 : std::ostringstream file_name;
762 :
763 : // Write out the temporal discretization data
764 5 : file_name.str("");
765 2 : file_name << directory_name << "/temporal_discretization_data" << suffix;
766 5 : Xdr temporal_discretization_data_out(file_name.str(), mode);
767 :
768 : Real real_value; unsigned int int_value;
769 3 : real_value = get_delta_t(); temporal_discretization_data_out << real_value;
770 3 : real_value = get_euler_theta(); temporal_discretization_data_out << real_value;
771 3 : int_value = get_n_time_steps(); temporal_discretization_data_out << int_value;
772 3 : int_value = get_time_step(); temporal_discretization_data_out << int_value;
773 3 : temporal_discretization_data_out.close();
774 :
775 :
776 : // Write out the L2 matrix
777 5 : file_name.str("");
778 2 : file_name << directory_name << "/RB_L2_matrix" << suffix;
779 6 : Xdr RB_L2_matrix_out(file_name.str(), mode);
780 :
781 63 : for (unsigned int i=0; i<n_bfs; i++)
782 : {
783 1260 : for (unsigned int j=0; j<n_bfs; j++)
784 : {
785 800 : RB_L2_matrix_out << RB_L2_matrix(i,j);
786 : }
787 : }
788 3 : RB_L2_matrix_out.close();
789 :
790 : // Write out the M_q matrices
791 6 : for (unsigned int q_m=0; q_m<Q_m; q_m++)
792 : {
793 5 : file_name.str("");
794 3 : file_name << directory_name << "/RB_M_";
795 : file_name << std::setw(3)
796 : << std::setprecision(0)
797 1 : << std::setfill('0')
798 1 : << std::right
799 1 : << q_m;
800 1 : file_name << suffix;
801 6 : Xdr RB_M_q_m_out(file_name.str(), mode);
802 :
803 63 : for (unsigned int i=0; i<n_bfs; i++)
804 : {
805 1260 : for (unsigned int j=0; j<n_bfs; j++)
806 : {
807 1200 : RB_M_q_m_out << RB_M_q_vector[q_m](i,j);
808 : }
809 : }
810 3 : RB_M_q_m_out.close();
811 1 : }
812 :
813 : // Write out the initial condition data
814 : // and the initial L2 error for all N
815 5 : file_name.str("");
816 2 : file_name << directory_name << "/initial_conditions" << suffix;
817 5 : Xdr initial_conditions_out(file_name.str(), mode);
818 5 : file_name.str("");
819 2 : file_name << directory_name << "/initial_L2_error" << suffix;
820 6 : Xdr initial_L2_error_out(file_name.str(), mode);
821 :
822 63 : for (unsigned int i=0; i<n_bfs; i++)
823 : {
824 80 : initial_L2_error_out << initial_L2_error_all_N[i];
825 690 : for (unsigned int j=0; j<=i; j++)
826 : {
827 630 : initial_conditions_out << RB_initial_condition_all_N[i](j);
828 : }
829 : }
830 3 : initial_conditions_out.close();
831 3 : initial_L2_error_out.close();
832 :
833 : // Next write out the Fq_Mq representor norm data
834 5 : file_name.str("");
835 2 : file_name << directory_name << "/Fq_Mq_terms" << suffix;
836 6 : Xdr RB_Fq_Mq_terms_out(file_name.str(), mode);
837 :
838 6 : for (unsigned int q_f=0; q_f<Q_f; q_f++)
839 : {
840 6 : for (unsigned int q_m=0; q_m<Q_m; q_m++)
841 : {
842 63 : for (unsigned int i=0; i<n_bfs; i++)
843 : {
844 120 : RB_Fq_Mq_terms_out << Fq_Mq_representor_innerprods[q_f][q_m][i];
845 : }
846 : }
847 : }
848 3 : RB_Fq_Mq_terms_out.close();
849 :
850 : // Next write out the Mq_Mq representor norm data
851 5 : file_name.str("");
852 2 : file_name << directory_name << "/Mq_Mq_terms" << suffix;
853 5 : Xdr RB_Mq_Mq_terms_out(file_name.str(), mode);
854 :
855 3 : unsigned int Q_m_hat = Q_m*(Q_m+1)/2;
856 6 : for (unsigned int q=0; q<Q_m_hat; q++)
857 : {
858 63 : for (unsigned int i=0; i<n_bfs; i++)
859 : {
860 1260 : for (unsigned int j=0; j<n_bfs; j++)
861 : {
862 2400 : RB_Mq_Mq_terms_out << Mq_Mq_representor_innerprods[q][i][j];
863 : }
864 : }
865 : }
866 3 : RB_Mq_Mq_terms_out.close();
867 :
868 : // Next write out the Aq_Mq representor norm data
869 5 : file_name.str("");
870 2 : file_name << directory_name << "/Aq_Mq_terms" << suffix;
871 6 : Xdr RB_Aq_Mq_terms_out(file_name.str(), mode);
872 :
873 12 : for (unsigned int q_a=0; q_a<Q_a; q_a++)
874 : {
875 18 : for (unsigned int q_m=0; q_m<Q_m; q_m++)
876 : {
877 189 : for (unsigned int i=0; i<n_bfs; i++)
878 : {
879 3780 : for (unsigned int j=0; j<n_bfs; j++)
880 : {
881 8400 : RB_Aq_Mq_terms_out << Aq_Mq_representor_innerprods[q_a][q_m][i][j];
882 : }
883 : }
884 : }
885 : }
886 3 : RB_Aq_Mq_terms_out.close();
887 1 : }
888 6 : }
889 :
890 6 : void TransientRBEvaluation::legacy_read_offline_data_from_files(const std::string & directory_name,
891 : bool read_error_bound_data,
892 : const bool read_binary_data)
893 : {
894 4 : LOG_SCOPE("legacy_read_offline_data_from_files()", "TransientRBEvaluation");
895 :
896 6 : Parent::legacy_read_offline_data_from_files(directory_name);
897 :
898 : TransientRBThetaExpansion & trans_theta_expansion =
899 6 : cast_ref<TransientRBThetaExpansion &>(get_rb_theta_expansion());
900 6 : const unsigned int Q_m = trans_theta_expansion.get_n_M_terms();
901 6 : const unsigned int Q_a = trans_theta_expansion.get_n_A_terms();
902 6 : const unsigned int Q_f = trans_theta_expansion.get_n_F_terms();
903 :
904 : // First, find out how many basis functions we had when Greedy terminated
905 : // This was set in RBSystem::read_offline_data_from_files
906 6 : unsigned int n_bfs = this->get_n_basis_functions();
907 :
908 : // The reading mode: DECODE for binary, READ for ASCII
909 6 : XdrMODE mode = read_binary_data ? DECODE : READ;
910 :
911 : // The suffix to use for all the files that are written out
912 8 : const std::string suffix = read_binary_data ? ".xdr" : ".dat";
913 :
914 : // The string stream we'll use to make the file names
915 10 : std::ostringstream file_name;
916 :
917 : // Write out the temporal discretization data
918 10 : file_name.str("");
919 4 : file_name << directory_name << "/temporal_discretization_data" << suffix;
920 10 : assert_file_exists(file_name.str());
921 :
922 12 : Xdr temporal_discretization_data_in(file_name.str(), mode);
923 :
924 : Real real_value; unsigned int int_value;
925 6 : temporal_discretization_data_in >> real_value; set_delta_t(real_value);
926 6 : temporal_discretization_data_in >> real_value; set_euler_theta(real_value);
927 6 : temporal_discretization_data_in >> int_value; set_n_time_steps(int_value);
928 6 : temporal_discretization_data_in >> int_value; set_time_step(int_value);
929 6 : temporal_discretization_data_in.close();
930 :
931 10 : file_name.str("");
932 4 : file_name << directory_name << "/RB_L2_matrix" << suffix;
933 10 : assert_file_exists(file_name.str());
934 :
935 12 : Xdr RB_L2_matrix_in(file_name.str(), mode);
936 :
937 126 : for (unsigned int i=0; i<n_bfs; i++)
938 : {
939 2520 : for (unsigned int j=0; j<n_bfs; j++)
940 : {
941 0 : Number value;
942 1600 : RB_L2_matrix_in >> value;
943 3200 : RB_L2_matrix(i,j) = value;
944 : }
945 : }
946 6 : RB_L2_matrix_in.close();
947 :
948 : // Read in the M_q matrices
949 12 : for (unsigned int q_m=0; q_m<Q_m; q_m++)
950 : {
951 10 : file_name.str("");
952 6 : file_name << directory_name << "/RB_M_";
953 : file_name << std::setw(3)
954 : << std::setprecision(0)
955 2 : << std::setfill('0')
956 2 : << std::right
957 2 : << q_m;
958 :
959 2 : file_name << suffix;
960 10 : assert_file_exists(file_name.str());
961 :
962 12 : Xdr RB_M_q_m_in(file_name.str(), mode);
963 :
964 126 : for (unsigned int i=0; i<n_bfs; i++)
965 : {
966 2520 : for (unsigned int j=0; j<n_bfs; j++)
967 : {
968 0 : Number value;
969 1600 : RB_M_q_m_in >> value;
970 3200 : RB_M_q_vector[q_m](i,j) = value;
971 : }
972 : }
973 6 : RB_M_q_m_in.close();
974 2 : }
975 :
976 :
977 : // Read in the initial condition data
978 : // and the initial L2 error for all N
979 10 : file_name.str("");
980 4 : file_name << directory_name << "/initial_conditions" << suffix;
981 10 : assert_file_exists(file_name.str());
982 :
983 10 : Xdr initial_conditions_in(file_name.str(), mode);
984 :
985 10 : file_name.str("");
986 4 : file_name << directory_name << "/initial_L2_error" << suffix;
987 10 : assert_file_exists(file_name.str());
988 :
989 12 : Xdr initial_L2_error_in(file_name.str(), mode);
990 :
991 126 : for (unsigned int i=0; i<n_bfs; i++)
992 : {
993 160 : initial_L2_error_in >> initial_L2_error_all_N[i];
994 1380 : for (unsigned int j=0; j<=i; j++)
995 : {
996 1260 : initial_conditions_in >> RB_initial_condition_all_N[i](j);
997 : }
998 : }
999 6 : initial_conditions_in.close();
1000 6 : initial_L2_error_in.close();
1001 :
1002 :
1003 6 : if (read_error_bound_data)
1004 : {
1005 : // Next read in the Fq_Mq representor norm data
1006 10 : file_name.str("");
1007 4 : file_name << directory_name << "/Fq_Mq_terms" << suffix;
1008 10 : assert_file_exists(file_name.str());
1009 :
1010 12 : Xdr RB_Fq_Mq_terms_in(file_name.str(), mode);
1011 :
1012 12 : for (unsigned int q_f=0; q_f<Q_f; q_f++)
1013 : {
1014 12 : for (unsigned int q_m=0; q_m<Q_m; q_m++)
1015 : {
1016 126 : for (unsigned int i=0; i<n_bfs; i++)
1017 : {
1018 240 : RB_Fq_Mq_terms_in >> Fq_Mq_representor_innerprods[q_f][q_m][i];
1019 : }
1020 : }
1021 : }
1022 6 : RB_Fq_Mq_terms_in.close();
1023 :
1024 : // Next read in the Mq_Mq representor norm data
1025 10 : file_name.str("");
1026 4 : file_name << directory_name << "/Mq_Mq_terms" << suffix;
1027 10 : assert_file_exists(file_name.str());
1028 :
1029 10 : Xdr RB_Mq_Mq_terms_in(file_name.str(), mode);
1030 :
1031 6 : unsigned int Q_m_hat = Q_m*(Q_m+1)/2;
1032 12 : for (unsigned int q=0; q<Q_m_hat; q++)
1033 : {
1034 126 : for (unsigned int i=0; i<n_bfs; i++)
1035 : {
1036 2520 : for (unsigned int j=0; j<n_bfs; j++)
1037 : {
1038 4800 : RB_Mq_Mq_terms_in >> Mq_Mq_representor_innerprods[q][i][j];
1039 : }
1040 : }
1041 : }
1042 6 : RB_Mq_Mq_terms_in.close();
1043 :
1044 : // Next read in the Aq_Mq representor norm data
1045 10 : file_name.str("");
1046 4 : file_name << directory_name << "/Aq_Mq_terms" << suffix;
1047 10 : assert_file_exists(file_name.str());
1048 :
1049 12 : Xdr RB_Aq_Mq_terms_in(file_name.str(), mode);
1050 :
1051 24 : for (unsigned int q_a=0; q_a<Q_a; q_a++)
1052 : {
1053 36 : for (unsigned int q_m=0; q_m<Q_m; q_m++)
1054 : {
1055 378 : for (unsigned int i=0; i<n_bfs; i++)
1056 : {
1057 7560 : for (unsigned int j=0; j<n_bfs; j++)
1058 : {
1059 16800 : RB_Aq_Mq_terms_in >> Aq_Mq_representor_innerprods[q_a][q_m][i][j];
1060 : }
1061 : }
1062 : }
1063 : }
1064 6 : RB_Aq_Mq_terms_in.close();
1065 2 : }
1066 8 : }
1067 :
1068 : }
|
