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_construction.h"
22 : #include "libmesh/transient_rb_evaluation.h"
23 : #include "libmesh/transient_rb_theta_expansion.h"
24 : #include "libmesh/transient_rb_assembly_expansion.h"
25 :
26 : // LibMesh includes
27 : #include "libmesh/numeric_vector.h"
28 : #include "libmesh/sparse_matrix.h"
29 : #include "libmesh/dof_map.h"
30 : #include "libmesh/libmesh_logging.h"
31 : #include "libmesh/linear_solver.h"
32 : #include "libmesh/equation_systems.h"
33 : #include "libmesh/exodusII_io.h"
34 : #include "libmesh/getpot.h"
35 : #include "libmesh/dense_matrix.h"
36 : #include "libmesh/dense_vector.h"
37 : #include "libmesh/xdr_cxx.h"
38 : #include "libmesh/timestamp.h"
39 :
40 : // TIMPI includes
41 : #include "timpi/communicator.h"
42 :
43 : // For checking for the existence of files
44 : #include <sys/stat.h>
45 :
46 : #include <fstream>
47 : #include <sstream>
48 :
49 : // Need SLEPc to get the POD eigenvalues
50 : #if defined(LIBMESH_HAVE_SLEPC)
51 : // LAPACK include (via SLEPc)
52 : #include "libmesh/ignore_warnings.h"
53 : #include <petscsys.h>
54 : #include <slepcblaslapack.h>
55 : #include "libmesh/restore_warnings.h"
56 : #endif // LIBMESH_HAVE_SLEPC
57 :
58 : namespace libMesh
59 : {
60 :
61 140 : TransientRBConstruction::TransientRBConstruction (EquationSystems & es,
62 : const std::string & name_in,
63 140 : const unsigned int number_in)
64 : : Parent(es, name_in, number_in),
65 132 : L2_matrix(SparseMatrix<Number>::build(es.comm())),
66 132 : non_dirichlet_L2_matrix(SparseMatrix<Number>::build(es.comm())),
67 132 : nonzero_initialization(false),
68 132 : compute_truth_projection_error(false),
69 132 : init_filename(""),
70 132 : POD_tol(-1.),
71 132 : max_truth_solves(-1),
72 272 : L2_assembly(nullptr)
73 : {
74 : // Indicate that we need to compute the RB
75 : // inner product matrix in this case
76 140 : compute_RB_inner_product = true;
77 :
78 : // We should not necessarily exit the greedy due to repeated parameters in
79 : // the transient case
80 140 : exit_on_repeated_greedy_parameters = false;
81 140 : }
82 :
83 396 : TransientRBConstruction::~TransientRBConstruction () = default;
84 :
85 0 : void TransientRBConstruction::clear()
86 : {
87 0 : Parent::clear();
88 :
89 : // clear the mass matrices
90 0 : M_q_vector.clear();
91 :
92 0 : if (store_non_dirichlet_operators)
93 0 : non_dirichlet_M_q_vector.clear();
94 :
95 : // clear the temporal_data
96 0 : temporal_data.clear();
97 0 : }
98 :
99 70 : void TransientRBConstruction::initialize_rb_construction(bool skip_matrix_assembly,
100 : bool skip_vector_assembly)
101 : {
102 : // Check that the theta and assembly objects are consistently sized
103 : #ifndef NDEBUG
104 : TransientRBThetaExpansion & trans_theta_expansion =
105 2 : cast_ref<TransientRBThetaExpansion &>(get_rb_theta_expansion());
106 :
107 : TransientRBAssemblyExpansion & trans_assembly_expansion =
108 2 : cast_ref<TransientRBAssemblyExpansion &>(get_rb_assembly_expansion());
109 : #endif
110 : // This assert only gets called if DEBUG is on
111 2 : libmesh_assert_equal_to (trans_theta_expansion.get_n_M_terms(), trans_assembly_expansion.get_n_M_terms());
112 :
113 70 : Parent::initialize_rb_construction(skip_matrix_assembly, skip_vector_assembly);
114 70 : }
115 :
116 70 : void TransientRBConstruction::process_parameters_file (const std::string & parameters_filename)
117 : {
118 70 : Parent::process_parameters_file(parameters_filename);
119 :
120 : // Read in data from parameters_filename
121 144 : GetPot infile(parameters_filename);
122 :
123 : // Read in the generic temporal discretization data
124 70 : process_temporal_parameters_file(parameters_filename);
125 :
126 : // Read in the data specific to Construction
127 70 : nonzero_initialization = infile("nonzero_initialization",nonzero_initialization);
128 70 : init_filename = infile("init_filename",init_filename);
129 :
130 70 : const Real POD_tol_in = infile("POD_tol", POD_tol);
131 70 : const int max_truth_solves_in = infile("max_truth_solves", max_truth_solves);
132 70 : const unsigned int delta_N_in = infile("delta_N", delta_N);
133 :
134 2 : set_POD_tol(POD_tol_in);
135 2 : set_max_truth_solves(max_truth_solves_in);
136 2 : set_delta_N(delta_N_in);
137 :
138 : // Pass the temporal discretization data to the RBEvaluation
139 70 : TransientRBEvaluation & trans_rb_eval = cast_ref<TransientRBEvaluation &>(get_rb_evaluation());
140 70 : trans_rb_eval.pull_temporal_discretization_data( *this );
141 70 : }
142 :
143 70 : void TransientRBConstruction::print_info() const
144 : {
145 70 : Parent::print_info();
146 :
147 2 : libMesh::out << std::endl << "TransientRBConstruction parameters:" << std::endl;
148 :
149 70 : if (is_rb_eval_initialized())
150 : {
151 : // Print out info that describes the current setup
152 : auto & trans_theta_expansion =
153 70 : cast_ref<const TransientRBThetaExpansion &>(get_rb_theta_expansion());
154 70 : libMesh::out << "Q_m: " << trans_theta_expansion.get_n_M_terms() << std::endl;
155 : }
156 : else
157 : {
158 0 : libMesh::out << "RBThetaExpansion member is not set yet" << std::endl;
159 : }
160 70 : libMesh::out << "Number of time-steps: " << get_n_time_steps() << std::endl;
161 70 : libMesh::out << "dt: " << get_delta_t() << std::endl;
162 70 : libMesh::out << "euler_theta (time discretization parameter): " << get_euler_theta() << std::endl;
163 70 : if (get_POD_tol() > 0.)
164 0 : libMesh::out << "POD_tol: " << get_POD_tol() << std::endl;
165 70 : if (max_truth_solves > 0)
166 0 : libMesh::out << "Maximum number of truth solves: " << max_truth_solves << std::endl;
167 4 : libMesh::out << "delta_N (number of basis functions to add each POD-Greedy step): " << get_delta_N() << std::endl;
168 70 : if (nonzero_initialization)
169 : {
170 0 : libMesh::out << "Reading initial condition from " << init_filename << std::endl;
171 : }
172 : else
173 : {
174 2 : libMesh::out << "Using zero initial condition" << std::endl;
175 : }
176 2 : libMesh::out << std::endl;
177 70 : }
178 :
179 70 : void TransientRBConstruction::allocate_data_structures()
180 : {
181 70 : Parent::allocate_data_structures();
182 :
183 : TransientRBThetaExpansion & trans_theta_expansion =
184 70 : cast_ref<TransientRBThetaExpansion &>(get_rb_theta_expansion());
185 70 : const unsigned int Q_m = trans_theta_expansion.get_n_M_terms();
186 70 : const unsigned int n_outputs = trans_theta_expansion.get_n_outputs();
187 :
188 : // Resize and allocate vectors for storing mesh-dependent data
189 70 : const unsigned int n_time_levels = get_n_time_steps()+1;
190 70 : temporal_data.resize(n_time_levels);
191 :
192 : // Resize vectors for storing mesh-dependent data but only
193 : // initialize if initialize_mesh_dependent_data == true
194 70 : M_q_vector.resize(Q_m);
195 :
196 : // Only initialize the mass matrices if we
197 : // are not in single-matrix mode
198 : {
199 2 : DofMap & dof_map = this->get_dof_map();
200 :
201 70 : dof_map.attach_matrix(*L2_matrix);
202 70 : L2_matrix->init();
203 70 : L2_matrix->zero();
204 :
205 140 : for (unsigned int q=0; q<Q_m; q++)
206 : {
207 : // Initialize the memory for the matrices
208 70 : M_q_vector[q] = SparseMatrix<Number>::build(this->comm());
209 72 : dof_map.attach_matrix(*M_q_vector[q]);
210 72 : M_q_vector[q]->init();
211 72 : M_q_vector[q]->zero();
212 : }
213 :
214 : // We also need to initialize a second set of non-Dirichlet operators
215 70 : if (store_non_dirichlet_operators)
216 : {
217 0 : dof_map.attach_matrix(*non_dirichlet_L2_matrix);
218 0 : non_dirichlet_L2_matrix->init();
219 0 : non_dirichlet_L2_matrix->zero();
220 :
221 0 : non_dirichlet_M_q_vector.resize(Q_m);
222 0 : for (unsigned int q=0; q<Q_m; q++)
223 : {
224 : // Initialize the memory for the matrices
225 0 : non_dirichlet_M_q_vector[q] = SparseMatrix<Number>::build(this->comm());
226 0 : dof_map.attach_matrix(*non_dirichlet_M_q_vector[q]);
227 0 : non_dirichlet_M_q_vector[q]->init();
228 0 : non_dirichlet_M_q_vector[q]->zero();
229 : }
230 : }
231 : }
232 :
233 7140 : for (unsigned int i=0; i<n_time_levels; i++)
234 : {
235 7070 : temporal_data[i] = NumericVector<Number>::build(this->comm());
236 7272 : temporal_data[i]->init (this->n_dofs(), this->n_local_dofs(), false, PARALLEL);
237 : }
238 :
239 : // and the truth output vectors
240 70 : truth_outputs_all_k.resize(n_outputs);
241 350 : for (unsigned int n=0; n<n_outputs; n++)
242 : {
243 288 : truth_outputs_all_k[n].resize(n_time_levels);
244 : }
245 :
246 : // This vector is for storing rhs entries for
247 : // computing the projection of the initial condition
248 : // into the RB space
249 70 : RB_ic_proj_rhs_all_N.resize(Nmax);
250 70 : }
251 :
252 70 : void TransientRBConstruction::assemble_affine_expansion(bool skip_matrix_assembly,
253 : bool skip_vector_assembly)
254 : {
255 : // Call parent's assembly functions
256 70 : Parent::assemble_affine_expansion(skip_matrix_assembly, skip_vector_assembly);
257 :
258 : // Now update RB_ic_proj_rhs_all_N if necessary.
259 : // This allows us to compute the L2 projection
260 : // of the initial condition into the RB space
261 : // so that we can continue to enrich a given RB
262 : // space.
263 70 : if (get_rb_evaluation().get_n_basis_functions() > 0)
264 : {
265 : // Load the initial condition into the solution vector
266 0 : initialize_truth();
267 :
268 0 : std::unique_ptr<NumericVector<Number>> temp1 = NumericVector<Number>::build(this->comm());
269 0 : temp1->init (this->n_dofs(), this->n_local_dofs(), false, PARALLEL);
270 :
271 : // First compute the right-hand side vector for the L2 projection
272 0 : L2_matrix->vector_mult(*temp1, *solution);
273 :
274 0 : for (unsigned int i=0; i<get_rb_evaluation().get_n_basis_functions(); i++)
275 : {
276 0 : RB_ic_proj_rhs_all_N(i) = temp1->dot(get_rb_evaluation().get_basis_function(i));
277 : }
278 0 : }
279 70 : }
280 :
281 70 : Real TransientRBConstruction::train_reduced_basis(const bool resize_rb_eval_data)
282 : {
283 138 : libmesh_error_msg_if(get_RB_training_type() == "POD",
284 : "POD RB training is not supported with TransientRBConstruction");
285 :
286 70 : compute_truth_projection_error = true;
287 70 : Real value = Parent::train_reduced_basis(resize_rb_eval_data);
288 70 : compute_truth_projection_error = false;
289 :
290 70 : return value;
291 : }
292 :
293 294770 : SparseMatrix<Number> * TransientRBConstruction::get_M_q(unsigned int q)
294 : {
295 : TransientRBThetaExpansion & trans_theta_expansion =
296 294770 : cast_ref<TransientRBThetaExpansion &>(get_rb_theta_expansion());
297 :
298 294770 : libmesh_error_msg_if(q >= trans_theta_expansion.get_n_M_terms(),
299 : "Error: We must have q < Q_m in get_M_q.");
300 :
301 303192 : return M_q_vector[q].get();
302 : }
303 :
304 0 : SparseMatrix<Number> * TransientRBConstruction::get_non_dirichlet_M_q(unsigned int q)
305 : {
306 0 : libmesh_error_msg_if(!store_non_dirichlet_operators,
307 : "Error: Must have store_non_dirichlet_operators==true to access non_dirichlet_M_q.");
308 :
309 : TransientRBThetaExpansion & trans_theta_expansion =
310 0 : cast_ref<TransientRBThetaExpansion &>(get_rb_theta_expansion());
311 :
312 0 : libmesh_error_msg_if(q >= trans_theta_expansion.get_n_M_terms(),
313 : "Error: We must have q < Q_m in get_M_q.");
314 :
315 0 : return non_dirichlet_M_q_vector[q].get();
316 : }
317 :
318 0 : void TransientRBConstruction::get_all_matrices(std::map<std::string, SparseMatrix<Number> *> & all_matrices)
319 : {
320 0 : Parent::get_all_matrices(all_matrices);
321 :
322 0 : all_matrices["L2_matrix"] = L2_matrix.get();
323 :
324 0 : if (store_non_dirichlet_operators)
325 0 : all_matrices["L2_matrix_non_dirichlet"] = non_dirichlet_L2_matrix.get();
326 :
327 : TransientRBThetaExpansion & trans_theta_expansion =
328 0 : cast_ref<TransientRBThetaExpansion &>(get_rb_theta_expansion());
329 0 : const unsigned int Q_m = trans_theta_expansion.get_n_M_terms();
330 :
331 0 : for (unsigned int q_m=0; q_m<Q_m; q_m++)
332 : {
333 0 : std::stringstream matrix_name;
334 0 : matrix_name << "M" << q_m;
335 0 : all_matrices[matrix_name.str()] = get_M_q(q_m);
336 :
337 0 : if (store_non_dirichlet_operators)
338 : {
339 0 : matrix_name << "_non_dirichlet";
340 0 : all_matrices[matrix_name.str()] = get_non_dirichlet_M_q(q_m);
341 : }
342 0 : }
343 0 : }
344 :
345 70 : void TransientRBConstruction::assemble_L2_matrix(SparseMatrix<Number> * input_matrix, bool apply_dirichlet_bc)
346 : {
347 70 : input_matrix->zero();
348 70 : add_scaled_matrix_and_vector(1.,
349 : L2_assembly,
350 : input_matrix,
351 : nullptr,
352 : false, /* symmetrize */
353 : apply_dirichlet_bc);
354 70 : }
355 :
356 0 : void TransientRBConstruction::assemble_mass_matrix(SparseMatrix<Number> * input_matrix)
357 : {
358 0 : input_matrix->zero();
359 0 : add_scaled_mass_matrix(1., input_matrix);
360 0 : }
361 :
362 140000 : void TransientRBConstruction::add_scaled_mass_matrix(Number scalar, SparseMatrix<Number> * input_matrix)
363 : {
364 140000 : const RBParameters & mu = get_parameters();
365 :
366 : TransientRBThetaExpansion & trans_theta_expansion =
367 140000 : cast_ref<TransientRBThetaExpansion &>(get_rb_theta_expansion());
368 :
369 140000 : const unsigned int Q_m = trans_theta_expansion.get_n_M_terms();
370 :
371 280000 : for (unsigned int q=0; q<Q_m; q++)
372 140000 : input_matrix->add(scalar * trans_theta_expansion.eval_M_theta(q,mu), *get_M_q(q));
373 140000 : }
374 :
375 140000 : void TransientRBConstruction::mass_matrix_scaled_matvec(Number scalar,
376 : NumericVector<Number> & dest,
377 : NumericVector<Number> & arg)
378 : {
379 8000 : LOG_SCOPE("mass_matrix_scaled_matvec()", "TransientRBConstruction");
380 :
381 140000 : dest.zero();
382 :
383 140000 : const RBParameters & mu = get_parameters();
384 :
385 : TransientRBThetaExpansion & trans_theta_expansion =
386 140000 : cast_ref<TransientRBThetaExpansion &>(get_rb_theta_expansion());
387 :
388 140000 : const unsigned int Q_m = trans_theta_expansion.get_n_M_terms();
389 :
390 144000 : std::unique_ptr<NumericVector<Number>> temp_vec = NumericVector<Number>::build(this->comm());
391 140000 : temp_vec->init (this->n_dofs(), this->n_local_dofs(), false, PARALLEL);
392 :
393 280000 : for (unsigned int q=0; q<Q_m; q++)
394 : {
395 140000 : get_M_q(q)->vector_mult(*temp_vec, arg);
396 144000 : dest.add(scalar * trans_theta_expansion.eval_M_theta(q,mu), *temp_vec);
397 : }
398 140000 : }
399 :
400 140000 : void TransientRBConstruction::truth_assembly()
401 : {
402 8000 : LOG_SCOPE("truth_assembly()", "TransientRBConstruction");
403 :
404 140000 : this->matrix->close();
405 :
406 140000 : this->matrix->zero();
407 140000 : this->rhs->zero();
408 :
409 140000 : const RBParameters & mu = get_parameters();
410 :
411 : TransientRBThetaExpansion & trans_theta_expansion =
412 140000 : cast_ref<TransientRBThetaExpansion &>(get_rb_theta_expansion());
413 :
414 140000 : const unsigned int Q_a = trans_theta_expansion.get_n_A_terms();
415 140000 : const unsigned int Q_f = trans_theta_expansion.get_n_F_terms();
416 :
417 140000 : const Real dt = get_delta_t();
418 140000 : const Real euler_theta = get_euler_theta();
419 :
420 : {
421 : // We should have already assembled the matrices
422 : // and vectors in the affine expansion, so
423 : // just use them
424 :
425 140000 : add_scaled_mass_matrix(1./dt, matrix);
426 140000 : mass_matrix_scaled_matvec(1./dt, *rhs, *current_local_solution);
427 :
428 144000 : std::unique_ptr<NumericVector<Number>> temp_vec = NumericVector<Number>::build(this->comm());
429 140000 : temp_vec->init (this->n_dofs(), this->n_local_dofs(), false, PARALLEL);
430 :
431 560000 : for (unsigned int q_a=0; q_a<Q_a; q_a++)
432 : {
433 420000 : matrix->add(euler_theta*trans_theta_expansion.eval_A_theta(q_a,mu), *get_Aq(q_a));
434 :
435 420000 : get_Aq(q_a)->vector_mult(*temp_vec, *current_local_solution);
436 420000 : temp_vec->scale( -(1.-euler_theta)*trans_theta_expansion.eval_A_theta(q_a,mu) );
437 432000 : rhs->add(*temp_vec);
438 : }
439 :
440 280000 : for (unsigned int q_f=0; q_f<Q_f; q_f++)
441 : {
442 140000 : *temp_vec = *get_Fq(q_f);
443 140000 : temp_vec->scale( get_control(get_time_step())*trans_theta_expansion.eval_F_theta(q_f,mu) );
444 144000 : rhs->add(*temp_vec);
445 : }
446 :
447 132000 : }
448 :
449 140000 : this->matrix->close();
450 140000 : this->rhs->close();
451 140000 : }
452 :
453 140 : void TransientRBConstruction::set_L2_assembly(ElemAssembly & L2_assembly_in)
454 : {
455 140 : L2_assembly = &L2_assembly_in;
456 140 : }
457 :
458 0 : ElemAssembly & TransientRBConstruction::get_L2_assembly()
459 : {
460 0 : libmesh_error_msg_if(!L2_assembly, "Error: L2_assembly hasn't been initialized yet");
461 :
462 0 : return *L2_assembly;
463 : }
464 :
465 70 : void TransientRBConstruction::assemble_Mq_matrix(unsigned int q, SparseMatrix<Number> * input_matrix, bool apply_dirichlet_bc)
466 : {
467 : TransientRBThetaExpansion & trans_theta_expansion =
468 70 : cast_ref<TransientRBThetaExpansion &>(get_rb_theta_expansion());
469 :
470 : TransientRBAssemblyExpansion & trans_assembly_expansion =
471 70 : cast_ref<TransientRBAssemblyExpansion &>(get_rb_assembly_expansion());
472 :
473 70 : libmesh_error_msg_if(q >= trans_theta_expansion.get_n_M_terms(),
474 : "Error: We must have q < Q_m in assemble_Mq_matrix.");
475 :
476 70 : input_matrix->zero();
477 72 : add_scaled_matrix_and_vector(1.,
478 70 : &trans_assembly_expansion.get_M_assembly(q),
479 : input_matrix,
480 : nullptr,
481 : false, /* symmetrize */
482 : apply_dirichlet_bc);
483 70 : }
484 :
485 70 : void TransientRBConstruction::assemble_all_affine_operators()
486 : {
487 70 : Parent::assemble_all_affine_operators();
488 :
489 : TransientRBThetaExpansion & trans_theta_expansion =
490 70 : cast_ref<TransientRBThetaExpansion &>(get_rb_theta_expansion());
491 :
492 140 : for (unsigned int q=0; q<trans_theta_expansion.get_n_M_terms(); q++)
493 70 : assemble_Mq_matrix(q, get_M_q(q));
494 :
495 70 : if (store_non_dirichlet_operators)
496 : {
497 0 : for (unsigned int q=0; q<trans_theta_expansion.get_n_M_terms(); q++)
498 0 : assemble_Mq_matrix(q, get_non_dirichlet_M_q(q), false);
499 : }
500 70 : }
501 :
502 70 : void TransientRBConstruction::assemble_misc_matrices()
503 : {
504 2 : libMesh::out << "Assembling L2 matrix" << std::endl;
505 70 : assemble_L2_matrix(L2_matrix.get());
506 :
507 70 : if (store_non_dirichlet_operators)
508 : {
509 0 : libMesh::out << "Assembling non-Dirichlet L2 matrix" << std::endl;
510 0 : assemble_L2_matrix(non_dirichlet_L2_matrix.get(), /* apply_dirichlet_bc = */ false);
511 : }
512 :
513 70 : Parent::assemble_misc_matrices();
514 70 : }
515 :
516 1400 : Real TransientRBConstruction::truth_solve(int write_interval)
517 : {
518 40 : LOG_SCOPE("truth_solve()", "TransientRBConstruction");
519 :
520 1400 : const RBParameters & mu = get_parameters();
521 1400 : const unsigned int n_time_steps = get_n_time_steps();
522 :
523 : // // NumericVector for computing true L2 error
524 : // std::unique_ptr<NumericVector<Number>> temp = NumericVector<Number>::build();
525 : // temp->init (this->n_dofs(), this->n_local_dofs(), false, PARALLEL);
526 :
527 : // Apply initial condition again.
528 1400 : initialize_truth();
529 1400 : set_time_step(0);
530 :
531 : // Now compute the truth outputs
532 7000 : for (unsigned int n=0; n<get_rb_theta_expansion().get_n_outputs(); n++)
533 : {
534 5760 : truth_outputs_all_k[n][0] = 0.;
535 11200 : for (unsigned int q_l=0; q_l<get_rb_theta_expansion().get_n_output_terms(n); q_l++)
536 : {
537 11360 : truth_outputs_all_k[n][0] += get_rb_theta_expansion().eval_output_theta(n,q_l,mu)*
538 5600 : get_output_vector(n,q_l)->dot(*solution);
539 : }
540 : }
541 :
542 : // Load initial projection error into temporal_data dense matrix
543 1400 : if (compute_truth_projection_error)
544 1400 : set_error_temporal_data();
545 :
546 141400 : for (unsigned int time_level=1; time_level<=n_time_steps; time_level++)
547 : {
548 140000 : set_time_step(time_level);
549 :
550 144000 : *old_local_solution = *current_local_solution;
551 :
552 : // We assume that the truth assembly has been attached to the system
553 140000 : truth_assembly();
554 :
555 : // truth_assembly assembles into matrix and rhs, so use those for the solve
556 140000 : solve_for_matrix_and_rhs(*get_linear_solver(), *matrix, *rhs);
557 :
558 : // The matrix doesn't change at each timestep, so we
559 : // can set reuse_preconditioner == true
560 140000 : linear_solver->reuse_preconditioner(true);
561 :
562 140000 : if (assert_convergence)
563 : {
564 140000 : check_convergence(*get_linear_solver());
565 : }
566 :
567 : // Now compute the truth outputs
568 700000 : for (unsigned int n=0; n<get_rb_theta_expansion().get_n_outputs(); n++)
569 : {
570 592000 : truth_outputs_all_k[n][time_level] = 0.;
571 1120000 : for (unsigned int q_l=0; q_l<get_rb_theta_expansion().get_n_output_terms(n); q_l++)
572 : {
573 576000 : truth_outputs_all_k[n][time_level] +=
574 560000 : get_rb_theta_expansion().eval_output_theta(n,q_l,mu)*get_output_vector(n,q_l)->dot(*solution);
575 : }
576 : }
577 :
578 : // load projection error into column _k of temporal_data matrix
579 140000 : if (compute_truth_projection_error)
580 140000 : set_error_temporal_data();
581 :
582 140000 : if ((write_interval > 0) && (time_level%write_interval == 0))
583 : {
584 0 : libMesh::out << std::endl << "Truth solve, plotting time step " << time_level << std::endl;
585 :
586 0 : std::ostringstream file_name;
587 :
588 0 : file_name << "truth.e.";
589 : file_name << std::setw(3)
590 : << std::setprecision(0)
591 0 : << std::setfill('0')
592 0 : << std::right
593 0 : << time_level;
594 :
595 : #ifdef LIBMESH_HAVE_EXODUS_API
596 0 : ExodusII_IO(get_mesh()).write_equation_systems (file_name.str(),
597 0 : this->get_equation_systems());
598 : #endif
599 0 : }
600 : }
601 :
602 : // Set reuse_preconditioner back to false for subsequent solves.
603 1400 : linear_solver->reuse_preconditioner(false);
604 :
605 : // Get the L2 norm of the truth solution at time-level _K
606 : // Useful for normalizing our true error data
607 1400 : L2_matrix->vector_mult(*inner_product_storage_vector, *solution);
608 1440 : Real final_truth_L2_norm = libmesh_real(std::sqrt(inner_product_storage_vector->dot(*solution)));
609 :
610 :
611 1440 : return final_truth_L2_norm;
612 : }
613 :
614 : bool
615 1400 : TransientRBConstruction::greedy_termination_test(Real abs_greedy_error,
616 : Real initial_greedy_error,
617 : int count)
618 : {
619 1400 : if ((get_max_truth_solves()>0) && (count >= get_max_truth_solves()))
620 : {
621 0 : libMesh::out << "Maximum number of truth solves reached: max = "
622 0 : << count << std::endl;
623 0 : return true;
624 : }
625 :
626 1400 : return Parent::greedy_termination_test(abs_greedy_error, initial_greedy_error, count);
627 : }
628 :
629 141400 : Number TransientRBConstruction::set_error_temporal_data()
630 : {
631 8080 : LOG_SCOPE("set_error_temporal_data()", "TransientRBConstruction");
632 :
633 : // first compute the projection of solution onto the current
634 : // RB space
635 :
636 141400 : const unsigned int time_step = get_time_step();
637 :
638 141400 : if (get_rb_evaluation().get_n_basis_functions() == 0)
639 : {
640 : // If the basis is empty, then the error is the solution itself
641 7272 : temporal_data[time_step]->zero();
642 7474 : temporal_data[time_step]->add(1., *solution);
643 : }
644 : else
645 : {
646 134330 : unsigned int RB_size = get_rb_evaluation().get_n_basis_functions();
647 :
648 138168 : std::unique_ptr<NumericVector<Number>> temp = NumericVector<Number>::build(this->comm());
649 134330 : temp->init (this->n_dofs(), this->n_local_dofs(), false, PARALLEL);
650 :
651 : // First compute the right-hand side vector for the projection
652 134330 : get_non_dirichlet_inner_product_matrix_if_avail()->vector_mult(*temp, *solution);
653 :
654 : // zero_dirichlet_dofs_on_vector(*temp);
655 :
656 : // Do not assume that RB_stiffness matrix is diagonal,
657 : // diagonality degrades as N increases
658 :
659 : // Get an appropriately sized copy of RB_inner_product_matrix
660 142006 : DenseMatrix<Number> RB_inner_product_matrix_N(RB_size,RB_size);
661 1477630 : for (unsigned int i=0; i<RB_size; i++)
662 18806200 : for (unsigned int j=0; j<RB_size; j++)
663 : {
664 17462900 : RB_inner_product_matrix_N(i,j) = get_rb_evaluation().RB_inner_product_matrix(i,j);
665 : }
666 :
667 : // Compute the projection RHS
668 134330 : DenseVector<Number> RB_proj_rhs(RB_size);
669 1477630 : for (unsigned int i=0; i<RB_size; i++)
670 : {
671 1343300 : RB_proj_rhs(i) = temp->dot(get_rb_evaluation().get_basis_function(i));
672 : }
673 :
674 134330 : DenseVector<Number> RB_proj(RB_size);
675 :
676 : // Now solve the linear system
677 134330 : RB_inner_product_matrix_N.lu_solve(RB_proj_rhs, RB_proj);
678 :
679 : // Load the RB projection into temp
680 134330 : temp->zero();
681 1477630 : for (unsigned int i=0; i<RB_size; i++)
682 : {
683 1343300 : temp->add(RB_proj(i), get_rb_evaluation().get_basis_function(i));
684 : }
685 :
686 138168 : temp->add(-1., *solution);
687 :
688 : // Now temp holds the projection error, store in temporal_data
689 142006 : *(temporal_data[time_step]) = *temp;
690 126654 : }
691 :
692 : // return the square of the X norm of the truth solution
693 141400 : get_non_dirichlet_inner_product_matrix_if_avail()->vector_mult(*inner_product_storage_vector,*solution);
694 :
695 153520 : return solution->dot(*inner_product_storage_vector);
696 : }
697 :
698 0 : const NumericVector<Number> & TransientRBConstruction::get_error_temporal_data()
699 : {
700 0 : LOG_SCOPE("get_error_temporal_data()", "TransientRBConstruction");
701 :
702 0 : const unsigned int time_step = get_time_step();
703 :
704 0 : return *temporal_data[time_step];
705 : }
706 :
707 2800 : void TransientRBConstruction::initialize_truth ()
708 : {
709 2800 : if (nonzero_initialization)
710 : {
711 : // Use System::read_serialized_data to read the initial condition
712 : // into this->solution
713 0 : Xdr IC_data(init_filename, READ);
714 0 : read_serialized_data(IC_data, false);
715 0 : }
716 : else
717 : {
718 : // Otherwise zero out the solution as a default
719 2800 : this->solution->zero();
720 : }
721 2800 : this->solution->close();
722 2800 : this->update();
723 2800 : }
724 :
725 0 : void TransientRBConstruction::add_IC_to_RB_space()
726 : {
727 0 : LOG_SCOPE("add_IC_to_RB_space()", "TransientRBConstruction");
728 :
729 0 : libmesh_error_msg_if(get_rb_evaluation().get_n_basis_functions() > 0,
730 : "Error: Should not call TransientRBConstruction::add_IC_to_RB_space() "
731 : "on a system that already contains basis functions.");
732 :
733 0 : libmesh_error_msg_if(!nonzero_initialization,
734 : "Error: Should not call TransientRBConstruction::add_IC_to_RB_space() "
735 : "when nonzero_initialization==false.");
736 :
737 0 : initialize_truth();
738 :
739 : // load the new basis function into the basis_functions vector.
740 0 : get_rb_evaluation().basis_functions.emplace_back(NumericVector<Number>::build(this->comm()));
741 0 : NumericVector<Number> & current_bf = get_rb_evaluation().get_basis_function(get_rb_evaluation().get_n_basis_functions()-1);
742 0 : current_bf.init (this->n_dofs(), this->n_local_dofs(), false, PARALLEL);
743 0 : current_bf = *solution;
744 :
745 : // We can just set the norm to 1.
746 0 : get_non_dirichlet_inner_product_matrix_if_avail()->vector_mult(*inner_product_storage_vector,*solution);
747 :
748 0 : Real current_bf_norm = libmesh_real(std::sqrt( current_bf.dot(*inner_product_storage_vector) ));
749 0 : current_bf.scale(1./current_bf_norm);
750 :
751 0 : unsigned int saved_delta_N = get_delta_N();
752 0 : set_delta_N(1);
753 0 : update_system();
754 0 : set_delta_N(saved_delta_N);
755 0 : }
756 :
757 1400 : void TransientRBConstruction::enrich_RB_space()
758 : {
759 : // Need SLEPc to get the POD eigenvalues
760 80 : LOG_SCOPE("enrich_RB_space()", "TransientRBConstruction");
761 :
762 : // With the "method of snapshots", the size of
763 : // the eigenproblem is determined by the number
764 : // of time-steps (rather than the number of spatial dofs).
765 1400 : unsigned int n_snapshots = temporal_data.size();
766 1480 : DenseMatrix<Number> correlation_matrix(n_snapshots,n_snapshots);
767 142800 : for (unsigned int i=0; i<n_snapshots; i++)
768 : {
769 145440 : get_non_dirichlet_inner_product_matrix_if_avail()->vector_mult(
770 141400 : *inner_product_storage_vector, *temporal_data[i]);
771 :
772 7352800 : for (unsigned int j=0; j<=i; j++)
773 : {
774 : // Scale the inner products by the number of time-steps to normalize the
775 : // POD energy norm appropriately
776 7623480 : Number inner_prod = (temporal_data[j]->dot(*inner_product_storage_vector)) /
777 7211400 : (Real)(get_n_time_steps()+1);
778 :
779 7211400 : correlation_matrix(i,j) = inner_prod;
780 7211400 : if(i != j)
781 : {
782 7070000 : correlation_matrix(j,i) = libmesh_conj(inner_prod);
783 : }
784 : }
785 : }
786 :
787 : // The POD can be formulated in terms of either the SVD or an eigenvalue problem.
788 : // Here we use the SVD of the correlation matrix to obtain the POD eigenvalues and
789 : // eigenvectors.
790 1400 : DenseVector<Real> sigma( n_snapshots );
791 1480 : DenseMatrix<Number> U( n_snapshots, n_snapshots );
792 1480 : DenseMatrix<Number> VT( n_snapshots, n_snapshots );
793 1400 : correlation_matrix.svd(sigma, U, VT );
794 :
795 40 : libMesh::out << std::endl << "POD singular values:" << std::endl;
796 4200 : for (unsigned int i=0; i<=1; i++)
797 : {
798 160 : libMesh::out << "singular value " << i << " = " << sigma(i) << std::endl;
799 : }
800 40 : libMesh::out << "..." << std::endl;
801 1440 : libMesh::out << "last singular value = " << sigma(n_snapshots-1) << std::endl;
802 40 : libMesh::out << std::endl;
803 :
804 : // Now load the new basis functions
805 40 : unsigned int j = 0;
806 : while (true)
807 : {
808 : // load the new basis function into the basis_functions vector.
809 1400 : get_rb_evaluation().basis_functions.emplace_back(NumericVector<Number>::build(this->comm()));
810 : NumericVector<Number> & current_bf =
811 1400 : get_rb_evaluation().get_basis_function(get_rb_evaluation().get_n_basis_functions()-1);
812 1400 : current_bf.init (this->n_dofs(), this->n_local_dofs(), false, PARALLEL);
813 1400 : current_bf.zero();
814 :
815 : // Perform the matrix multiplication of temporal data with
816 : // the next POD eigenvector
817 142800 : for (unsigned int i=0; i<n_snapshots; i++)
818 : {
819 145440 : current_bf.add( U.el(i,j), *temporal_data[i] );
820 : }
821 :
822 : // We just set the norm to 1.
823 1400 : get_non_dirichlet_inner_product_matrix_if_avail()->vector_mult(*inner_product_storage_vector,current_bf);
824 :
825 1440 : Real current_bf_norm = std::abs( std::sqrt( current_bf.dot(*inner_product_storage_vector) ) );
826 1400 : current_bf.scale(1./current_bf_norm);
827 :
828 : // Increment j here since we use the incremented counter
829 : // in the if clauses below
830 1400 : j++;
831 :
832 : // If positive POD_tol, we use it to determine the number of basis functions
833 : // to add, and then break the loop when POD_tol is satisfied, or after Nmax
834 : // basis functions have been added. Else we break the loop after delta_N
835 : // (or Nmax) new basis functions.
836 1400 : if (POD_tol > 0.)
837 : {
838 0 : set_delta_N(1);
839 :
840 : // We need to define the updated RB system matrices before the RB solve
841 0 : update_system();
842 0 : Real error_bound = get_rb_evaluation().rb_solve(get_rb_evaluation().get_n_basis_functions());
843 :
844 0 : if ((error_bound <= POD_tol) || (get_rb_evaluation().get_n_basis_functions()==get_Nmax()))
845 : {
846 0 : set_delta_N(0);
847 0 : break;
848 : }
849 : }
850 : else
851 : {
852 1400 : if (j == get_delta_N())
853 : {
854 40 : break;
855 : }
856 : else
857 0 : if (get_rb_evaluation().get_n_basis_functions()==get_Nmax())
858 : {
859 0 : set_delta_N(j);
860 0 : break;
861 : }
862 : }
863 0 : }
864 2720 : }
865 :
866 :
867 1400 : void TransientRBConstruction::update_system()
868 : {
869 : // If delta_N is set to zero, there is nothing to update
870 1400 : if (get_delta_N() == 0)
871 0 : return;
872 :
873 1400 : Parent::update_system();
874 :
875 40 : libMesh::out << "Updating RB initial conditions" << std::endl;
876 1400 : update_RB_initial_condition_all_N();
877 : }
878 :
879 560 : SparseMatrix<Number> & TransientRBConstruction::get_matrix_for_output_dual_solves()
880 : {
881 560 : return *L2_matrix;
882 : }
883 :
884 70 : void TransientRBConstruction::load_rb_solution()
885 : {
886 4 : LOG_SCOPE("load_rb_solution()", "TransientRBConstruction");
887 :
888 70 : solution->zero();
889 :
890 70 : const unsigned int time_step = get_time_step();
891 :
892 70 : TransientRBEvaluation & trans_rb_eval = cast_ref<TransientRBEvaluation &>(get_rb_evaluation());
893 72 : DenseVector<Number> RB_solution_vector_k = trans_rb_eval.RB_temporal_solution_data[time_step];
894 :
895 70 : libmesh_error_msg_if(RB_solution_vector_k.size() > get_rb_evaluation().get_n_basis_functions(),
896 : "ERROR: rb_eval object contains "
897 : << get_rb_evaluation().get_n_basis_functions()
898 : << " basis functions. RB_solution vector contains "
899 : << RB_solution_vector_k.size()
900 : << " entries. RB_solution in TransientRBConstruction::load_rb_solution is too long!");
901 :
902 1470 : for (unsigned int i=0; i<RB_solution_vector_k.size(); i++)
903 1400 : solution->add(RB_solution_vector_k(i), get_rb_evaluation().get_basis_function(i));
904 :
905 70 : update();
906 70 : }
907 :
908 1400 : void TransientRBConstruction::update_RB_system_matrices()
909 : {
910 80 : LOG_SCOPE("update_RB_system_matrices()", "TransientRBConstruction");
911 :
912 1400 : Parent::update_RB_system_matrices();
913 :
914 1400 : TransientRBEvaluation & trans_rb_eval = cast_ref<TransientRBEvaluation &>(get_rb_evaluation());
915 :
916 : TransientRBThetaExpansion & trans_theta_expansion =
917 1400 : cast_ref<TransientRBThetaExpansion &>(get_rb_theta_expansion());
918 1400 : const unsigned int Q_m = trans_theta_expansion.get_n_M_terms();
919 :
920 1400 : unsigned int RB_size = get_rb_evaluation().get_n_basis_functions();
921 :
922 1440 : std::unique_ptr<NumericVector<Number>> temp = NumericVector<Number>::build(this->comm());
923 1400 : temp->init (this->n_dofs(), this->n_local_dofs(), false, PARALLEL);
924 :
925 2800 : for (unsigned int i=(RB_size-delta_N); i<RB_size; i++)
926 : {
927 16100 : for (unsigned int j=0; j<RB_size; j++)
928 : {
929 : // Compute reduced L2 matrix
930 14700 : temp->zero();
931 14700 : L2_matrix->vector_mult(*temp, get_rb_evaluation().get_basis_function(j));
932 :
933 14700 : Number value = get_rb_evaluation().get_basis_function(i).dot(*temp);
934 14700 : trans_rb_eval.RB_L2_matrix(i,j) = value;
935 14700 : if (i!=j)
936 : {
937 : // The L2 matrix is assumed
938 : // to be symmetric
939 13300 : trans_rb_eval.RB_L2_matrix(j,i) = value;
940 : }
941 :
942 29400 : for (unsigned int q_m=0; q_m<Q_m; q_m++)
943 : {
944 : // Compute reduced M_q matrix
945 14700 : temp->zero();
946 14700 : get_M_q(q_m)->vector_mult(*temp, get_rb_evaluation().get_basis_function(j));
947 :
948 14700 : value = (get_rb_evaluation().get_basis_function(i)).dot(*temp);
949 14700 : trans_rb_eval.RB_M_q_vector[q_m](i,j) = value;
950 :
951 14700 : if (i!=j)
952 : {
953 : // Each mass matrix term is assumed
954 : // to be symmetric
955 13300 : trans_rb_eval.RB_M_q_vector[q_m](j,i) = value;
956 : }
957 : }
958 :
959 : }
960 : }
961 1400 : }
962 :
963 :
964 :
965 :
966 1330 : void TransientRBConstruction::update_residual_terms(bool compute_inner_products)
967 : {
968 76 : LOG_SCOPE("update_residual_terms()", "TransientRBConstruction");
969 :
970 1330 : Parent::update_residual_terms(compute_inner_products);
971 :
972 1330 : TransientRBEvaluation & trans_rb_eval = cast_ref<TransientRBEvaluation &>(get_rb_evaluation());
973 :
974 : TransientRBThetaExpansion & trans_theta_expansion =
975 1330 : cast_ref<TransientRBThetaExpansion &>(get_rb_theta_expansion());
976 :
977 1330 : const unsigned int Q_m = trans_theta_expansion.get_n_M_terms();
978 1330 : const unsigned int Q_a = trans_theta_expansion.get_n_A_terms();
979 1330 : const unsigned int Q_f = trans_theta_expansion.get_n_F_terms();
980 :
981 1330 : unsigned int RB_size = get_rb_evaluation().get_n_basis_functions();
982 :
983 2660 : for (unsigned int q_m=0; q_m<Q_m; q_m++)
984 : {
985 2660 : for (unsigned int i=(RB_size-delta_N); i<RB_size; i++)
986 : {
987 : // Initialize the vectors when we need them
988 1406 : if (!trans_rb_eval.M_q_representor[q_m][i])
989 : {
990 2584 : trans_rb_eval.M_q_representor[q_m][i] = NumericVector<Number>::build(this->comm());
991 1406 : trans_rb_eval.M_q_representor[q_m][i]->init (this->n_dofs(), this->n_local_dofs(), false, PARALLEL);
992 : }
993 :
994 38 : libmesh_assert(trans_rb_eval.M_q_representor[q_m][i]->size() == this->n_dofs() &&
995 : trans_rb_eval.M_q_representor[q_m][i]->local_size() == this->n_local_dofs() );
996 :
997 1330 : rhs->zero();
998 1368 : M_q_vector[q_m]->vector_mult(*rhs, get_rb_evaluation().get_basis_function(i));
999 :
1000 1330 : if (!is_quiet())
1001 0 : libMesh::out << "Starting solve i="
1002 0 : << i << " in TransientRBConstruction::update_residual_terms() at "
1003 0 : << Utility::get_timestamp() << std::endl;
1004 :
1005 1368 : solve_for_matrix_and_rhs(*inner_product_solver, *inner_product_matrix, *rhs);
1006 :
1007 1330 : if (assert_convergence)
1008 1330 : check_convergence(*inner_product_solver);
1009 :
1010 1330 : if (!is_quiet())
1011 : {
1012 0 : libMesh::out << "Finished solve i="
1013 0 : << i << " in TransientRBConstruction::update_residual_terms() at "
1014 0 : << Utility::get_timestamp() << std::endl;
1015 :
1016 0 : libMesh::out << this->n_linear_iterations()
1017 0 : << " iterations, final residual "
1018 0 : << this->final_linear_residual() << std::endl;
1019 : }
1020 :
1021 1444 : *trans_rb_eval.M_q_representor[q_m][i] = *solution;
1022 : }
1023 : }
1024 :
1025 : // Now compute and store the inner products if requested
1026 1330 : if (compute_inner_products)
1027 : {
1028 2660 : for (unsigned int q_f=0; q_f<Q_f; q_f++)
1029 : {
1030 1330 : get_non_dirichlet_inner_product_matrix_if_avail()->vector_mult(*inner_product_storage_vector, *Fq_representor[q_f]);
1031 :
1032 2660 : for (unsigned int i=(RB_size-delta_N); i<RB_size; i++)
1033 : {
1034 2660 : for (unsigned int q_m=0; q_m<Q_m; q_m++)
1035 : {
1036 1444 : trans_rb_eval.Fq_Mq_representor_innerprods[q_f][q_m][i] =
1037 1444 : trans_rb_eval.M_q_representor[q_m][i]->dot(*inner_product_storage_vector);
1038 : } // end for q_m
1039 : } // end for i
1040 : } // end for q_f
1041 :
1042 38 : unsigned int q=0;
1043 2660 : for (unsigned int q_m1=0; q_m1<Q_m; q_m1++)
1044 : {
1045 2660 : for (unsigned int q_m2=q_m1; q_m2<Q_m; q_m2++)
1046 : {
1047 2660 : for (unsigned int i=(RB_size-delta_N); i<RB_size; i++)
1048 : {
1049 14630 : for (unsigned int j=0; j<RB_size; j++)
1050 : {
1051 13300 : get_non_dirichlet_inner_product_matrix_if_avail()->vector_mult(*inner_product_storage_vector, *trans_rb_eval.M_q_representor[q_m2][j]);
1052 :
1053 14060 : trans_rb_eval.Mq_Mq_representor_innerprods[q][i][j] =
1054 14440 : trans_rb_eval.M_q_representor[q_m1][i]->dot(*inner_product_storage_vector);
1055 :
1056 13300 : if (i != j)
1057 : {
1058 12312 : get_non_dirichlet_inner_product_matrix_if_avail()->vector_mult(*inner_product_storage_vector,
1059 1026 : *trans_rb_eval.M_q_representor[q_m2][i]);
1060 :
1061 12996 : trans_rb_eval.Mq_Mq_representor_innerprods[q][j][i] =
1062 12996 : trans_rb_eval.M_q_representor[q_m1][j]->dot(*inner_product_storage_vector);
1063 : }
1064 : } // end for j
1065 : } // end for i
1066 1330 : q++;
1067 : } // end for q_m2
1068 : } // end for q_m1
1069 :
1070 :
1071 2660 : for (unsigned int i=(RB_size-delta_N); i<RB_size; i++)
1072 : {
1073 14630 : for (unsigned int j=0; j<RB_size; j++)
1074 : {
1075 53200 : for (unsigned int q_a=0; q_a<Q_a; q_a++)
1076 : {
1077 79800 : for (unsigned int q_m=0; q_m<Q_m; q_m++)
1078 : {
1079 41040 : get_non_dirichlet_inner_product_matrix_if_avail()->vector_mult(*inner_product_storage_vector,
1080 41040 : *trans_rb_eval.M_q_representor[q_m][j]);
1081 :
1082 43320 : trans_rb_eval.Aq_Mq_representor_innerprods[q_a][q_m][i][j] =
1083 43320 : trans_rb_eval.Aq_representor[q_a][i]->dot(*inner_product_storage_vector);
1084 :
1085 39900 : if (i != j)
1086 : {
1087 36936 : get_non_dirichlet_inner_product_matrix_if_avail()->vector_mult(*inner_product_storage_vector,
1088 3078 : *trans_rb_eval.M_q_representor[q_m][i]);
1089 :
1090 40014 : trans_rb_eval.Aq_Mq_representor_innerprods[q_a][q_m][j][i] =
1091 38988 : trans_rb_eval.Aq_representor[q_a][j]->dot(*inner_product_storage_vector);
1092 : }
1093 : } // end for q_m
1094 : } // end for q_a
1095 : } // end for j
1096 : } // end for i
1097 : } // end if (compute_inner_products)
1098 1330 : }
1099 :
1100 :
1101 1400 : void TransientRBConstruction::update_RB_initial_condition_all_N()
1102 : {
1103 80 : LOG_SCOPE("update_RB_initial_condition_all_N()", "TransientRBConstruction");
1104 :
1105 1400 : TransientRBEvaluation & trans_rb_eval = cast_ref<TransientRBEvaluation &>(get_rb_evaluation());
1106 :
1107 : // Load the initial condition into the solution vector
1108 1400 : initialize_truth();
1109 :
1110 1440 : std::unique_ptr<NumericVector<Number>> temp1 = NumericVector<Number>::build(this->comm());
1111 1400 : temp1->init (this->n_dofs(), this->n_local_dofs(), false, PARALLEL);
1112 :
1113 1440 : std::unique_ptr<NumericVector<Number>> temp2 = NumericVector<Number>::build(this->comm());
1114 1400 : temp2->init (this->n_dofs(), this->n_local_dofs(), false, PARALLEL);
1115 :
1116 :
1117 1400 : unsigned int RB_size = get_rb_evaluation().get_n_basis_functions();
1118 :
1119 : // First compute the right-hand side vector for the L2 projection
1120 1400 : L2_matrix->vector_mult(*temp1, *solution);
1121 :
1122 2800 : for (unsigned int i=(RB_size-delta_N); i<RB_size; i++)
1123 : {
1124 1400 : RB_ic_proj_rhs_all_N(i) = temp1->dot(get_rb_evaluation().get_basis_function(i));
1125 : }
1126 :
1127 :
1128 : // Now compute the projection for each N
1129 1480 : DenseMatrix<Number> RB_L2_matrix_N;
1130 1400 : DenseVector<Number> RB_rhs_N;
1131 2800 : for (unsigned int N=(RB_size-delta_N); N<RB_size; N++)
1132 : {
1133 : // We have to index here by N+1 since the loop index is zero-based.
1134 1400 : trans_rb_eval.RB_L2_matrix.get_principal_submatrix(N+1, RB_L2_matrix_N);
1135 :
1136 1400 : RB_ic_proj_rhs_all_N.get_principal_subvector(N+1, RB_rhs_N);
1137 :
1138 1360 : DenseVector<Number> RB_ic_N(N+1);
1139 :
1140 : // Now solve the linear system
1141 1400 : RB_L2_matrix_N.lu_solve(RB_rhs_N, RB_ic_N);
1142 :
1143 : // Load RB_ic_N into RB_initial_condition_all_N
1144 1400 : trans_rb_eval.RB_initial_condition_all_N[N] = RB_ic_N;
1145 :
1146 : // Compute the L2 error for the RB initial condition
1147 : // This part is dependent on the truth space.
1148 :
1149 : // load the RB solution into temp1
1150 1400 : temp1->zero();
1151 16100 : for (unsigned int i=0; i<N+1; i++)
1152 : {
1153 14700 : temp1->add(RB_ic_N(i), get_rb_evaluation().get_basis_function(i));
1154 : }
1155 :
1156 : // subtract truth initial condition from RB_ic_N
1157 1440 : temp1->add(-1., *solution);
1158 :
1159 : // Compute L2 norm error, i.e. sqrt(M(solution,solution))
1160 1400 : temp2->zero();
1161 1400 : L2_matrix->vector_mult(*temp2, *temp1);
1162 :
1163 1440 : trans_rb_eval.initial_L2_error_all_N[N] = libmesh_real(std::sqrt(temp2->dot(*temp1)));
1164 : }
1165 1400 : }
1166 :
1167 : //Real TransientRBConstruction::uncached_compute_residual_dual_norm(const unsigned int N)
1168 : //{
1169 : // LOG_SCOPE("uncached_compute_residual_dual_norm()", "TransientRBConstruction");
1170 : //
1171 : // // This is the "slow" way of computing the residual, but it is useful
1172 : // // for validating the "fast" way.
1173 : // // Note that this only works in serial since otherwise each processor will
1174 : // // have a different parameter value during the Greedy training.
1175 : //
1176 : // // Assemble the right-hand side to find the Reisz representor
1177 : // // of the residual in the X norm
1178 : // std::unique_ptr<NumericVector<Number>> RB_sol = NumericVector<Number>::build();
1179 : // RB_sol->init (this->n_dofs(), this->n_local_dofs(), false, PARALLEL);
1180 : // RB_sol->zero();
1181 : //
1182 : // std::unique_ptr<NumericVector<Number>> ghosted_temp = NumericVector<Number>::build();
1183 : // ghosted_temp->init (this->n_dofs(), this->n_local_dofs(),
1184 : // this->get_dof_map().get_send_list(), false,
1185 : // GHOSTED);
1186 : //
1187 : // std::unique_ptr<NumericVector<Number>> parallel_temp = NumericVector<Number>::build();
1188 : // parallel_temp->init (this->n_dofs(), this->n_local_dofs(), false, PARALLEL);
1189 : //
1190 : // // Store current_local_solution, since we don't want to corrupt it
1191 : // *ghosted_temp = *current_local_solution;
1192 : //
1193 : // for (unsigned int i=0; i<N; i++)
1194 : // {
1195 : // RB_sol->add(RB_solution(i), rb_eval->get_basis_function(i));
1196 : // parallel_temp->add(old_RB_solution(i), rb_eval->get_basis_function(i));
1197 : // }
1198 : //
1199 : // // Load parallel_temp into current_local_solution in order to do assembly
1200 : // const std::vector<unsigned int> & send_list = this->get_dof_map().get_send_list();
1201 : // parallel_temp->localize (*current_local_solution, send_list);
1202 : //
1203 : // // Load the system_matrix
1204 : // this->truth_assembly();
1205 : //
1206 : // // Restore current_local_solution
1207 : // *current_local_solution = *ghosted_temp;
1208 : //
1209 : // matrix->vector_mult(*parallel_temp, *RB_sol);
1210 : // rhs->add(-1., *parallel_temp);
1211 : // rhs->close();
1212 : //
1213 : // zero_dirichlet_dofs_on_rhs();
1214 : //
1215 : // // Then solve the system to get the Reisz representor
1216 : // matrix->zero();
1217 : // {
1218 : // matrix->add(1., *inner_product_matrix);
1219 : // if (constrained_problem)
1220 : // matrix->add(1., *constraint_matrix);
1221 : // }
1222 : //
1223 : //
1224 : // solution->zero();
1225 : // solve();
1226 : // // Make sure we didn't max out the number of iterations
1227 : // if ((this->n_linear_iterations() >=
1228 : // this->get_equation_systems().parameters.get<unsigned int>("linear solver maximum iterations")) &&
1229 : // (this->final_linear_residual() >
1230 : // this->get_equation_systems().parameters.get<Real>("linear solver tolerance")))
1231 : // {
1232 : // libmesh_error_msg("Warning: Linear solver may not have converged! Final linear residual = "
1233 : // << this->final_linear_residual() << ", number of iterations = "
1234 : // << this->n_linear_iterations());
1235 : // }
1236 : //
1237 : // get_non_dirichlet_inner_product_matrix_if_avail()->vector_mult(*inner_product_storage_vector, *solution);
1238 : //
1239 : // Real slow_residual_norm_sq = solution->dot(*inner_product_storage_vector);
1240 : //
1241 : // return libmesh_real(std::sqrt( slow_residual_norm_sq ));
1242 : //}
1243 :
1244 0 : void TransientRBConstruction::write_riesz_representors_to_files(const std::string & riesz_representors_dir,
1245 : const bool write_binary_residual_representors)
1246 : {
1247 0 : LOG_SCOPE("write_riesz_representors_to_files()", "TransientRBConstruction");
1248 :
1249 : // Write out the M_q_representors. These are useful to have when restarting,
1250 : // so you don't have to recompute them all over again. There should be
1251 : // this->rb_eval->get_n_basis_functions() of these.
1252 0 : libMesh::out << "Writing out the M_q_representors..." << std::endl;
1253 :
1254 0 : std::ostringstream file_name;
1255 0 : const std::string riesz_representor_suffix = (write_binary_residual_representors ? ".xdr" : ".dat");
1256 :
1257 0 : TransientRBEvaluation & trans_rb_eval = cast_ref<TransientRBEvaluation &>(get_rb_evaluation());
1258 :
1259 0 : const unsigned int istop = trans_rb_eval.get_n_basis_functions();
1260 0 : const unsigned int istart = istop-get_delta_N();
1261 :
1262 0 : for (std::size_t q=0; q<trans_rb_eval.M_q_representor.size(); ++q)
1263 0 : for (unsigned int i=istart; i<istop; ++i)
1264 : {
1265 0 : libMesh::out << "Writing out M_q_representor[" << q << "][" << i << "]..." << std::endl;
1266 0 : libmesh_assert(trans_rb_eval.M_q_representor[q][i]);
1267 :
1268 0 : file_name.str(""); // reset filename
1269 0 : file_name << riesz_representors_dir << "/M_q_representor" << i << riesz_representor_suffix;
1270 :
1271 : {
1272 : // No need to copy!
1273 : //*solution = *(M_q_representor[q][i]);
1274 0 : trans_rb_eval.M_q_representor[q][i]->swap(*solution);
1275 :
1276 0 : Xdr mr_data(file_name.str(),
1277 0 : write_binary_residual_representors ? ENCODE : WRITE);
1278 :
1279 0 : write_serialized_data(mr_data, false);
1280 :
1281 : // Synchronize before moving on
1282 0 : this->comm().barrier();
1283 :
1284 : // Swap back.
1285 0 : trans_rb_eval.M_q_representor[q][i]->swap(*solution);
1286 :
1287 : // TODO: bzip the resulting file? See $LIBMESH_DIR/src/mesh/unstructured_mesh.C
1288 : // for the system call, be sure to do it only on one processor, etc.
1289 0 : }
1290 : }
1291 0 : }
1292 :
1293 0 : void TransientRBConstruction::read_riesz_representors_from_files(const std::string & riesz_representors_dir,
1294 : const bool read_binary_residual_representors)
1295 : {
1296 0 : LOG_SCOPE("read_riesz_representors_from_files()", "TransientRBConstruction");
1297 :
1298 : const std::string riesz_representor_suffix =
1299 0 : (read_binary_residual_representors ? ".xdr" : ".dat");
1300 :
1301 0 : std::ostringstream file_name;
1302 : struct stat stat_info;
1303 :
1304 0 : TransientRBEvaluation & trans_rb_eval = cast_ref<TransientRBEvaluation &>(get_rb_evaluation());
1305 :
1306 0 : libMesh::out << "Reading in the M_q_representors..." << std::endl;
1307 :
1308 : // Read in the Aq representors. The class makes room for [Q_m][Nmax] of these. We are going to
1309 : // read in [Q_m][this->rb_eval->get_n_basis_functions()]. FIXME:
1310 : // should we be worried about leaks in the locations where we're about to fill entries?
1311 0 : for (std::size_t i=0; i<trans_rb_eval.M_q_representor.size(); ++i)
1312 0 : for (std::size_t j=0; j<trans_rb_eval.M_q_representor[i].size(); ++j)
1313 0 : libmesh_error_msg_if(trans_rb_eval.M_q_representor[i][j] != nullptr,
1314 : "Error, must delete existing M_q_representor before reading in from file.");
1315 :
1316 : // Now ready to read them in from file!
1317 0 : for (std::size_t i=0; i<trans_rb_eval.M_q_representor.size(); ++i)
1318 0 : for (std::size_t j=0; j<trans_rb_eval.get_n_basis_functions(); ++j)
1319 : {
1320 0 : file_name.str(""); // reset filename
1321 : file_name << riesz_representors_dir
1322 0 : << "/M_q_representor" << i << "_" << j << riesz_representor_suffix;
1323 :
1324 : // On processor zero check to be sure the file exists
1325 0 : if (this->processor_id() == 0)
1326 : {
1327 0 : int stat_result = stat(file_name.str().c_str(), &stat_info);
1328 :
1329 0 : libmesh_error_msg_if(stat_result != 0, "File does not exist: " << file_name.str());
1330 : }
1331 :
1332 0 : Xdr aqr_data(file_name.str(),
1333 0 : read_binary_residual_representors ? DECODE : READ);
1334 :
1335 0 : read_serialized_data(aqr_data, false);
1336 :
1337 0 : trans_rb_eval.M_q_representor[i][j] = NumericVector<Number>::build(this->comm());
1338 0 : trans_rb_eval.M_q_representor[i][j]->init (n_dofs(), n_local_dofs(),
1339 0 : false, PARALLEL);
1340 :
1341 : // No need to copy, just swap
1342 : //*M_q_representor[i][j] = *solution;
1343 0 : trans_rb_eval.M_q_representor[i][j]->swap(*solution);
1344 0 : }
1345 0 : }
1346 :
1347 : } // namespace libMesh
|
