libMesh
rb_eim_evaluation.C
Go to the documentation of this file.
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 // C++ includes
21 #include <sstream>
22 #include <fstream>
23 
24 // rbOOmit includes
25 #include "libmesh/rb_eim_evaluation.h"
26 #include "libmesh/rb_eim_theta.h"
27 #include "libmesh/rb_parametrized_function.h"
28 
29 // libMesh includes
30 #include "libmesh/xdr_cxx.h"
31 #include "libmesh/libmesh_logging.h"
32 #include "libmesh/replicated_mesh.h"
33 #include "libmesh/elem.h"
34 
35 namespace libMesh
36 {
37 
39  :
40  RBEvaluation(comm_in),
41  _previous_N(0),
42  _previous_error_bound(-1),
43  _interpolation_points_mesh(comm_in)
44 {
45  // Indicate that we need to compute the RB
46  // inner product matrix in this case
48 
49  // initialize to the empty RBThetaExpansion object
51 
52  // Let's not renumber the _interpolation_points_mesh
54 }
55 
57 {
58  this->clear();
59 }
60 
62 {
63  Parent::clear();
64 
65  interpolation_points.clear();
69 
70  // Delete any RBTheta objects that were created
71  for (std::size_t i=0; i<_rb_eim_theta_objects.size(); i++)
72  delete _rb_eim_theta_objects[i];
73  _rb_eim_theta_objects.clear();
74 }
75 
76 void RBEIMEvaluation::resize_data_structures(const unsigned int Nmax,
77  bool resize_error_bound_data)
78 {
79  Parent::resize_data_structures(Nmax, resize_error_bound_data);
80 
81  // Resize the data structures relevant to the EIM system
82  interpolation_points.clear();
85  interpolation_matrix.resize(Nmax,Nmax);
86 }
87 
89 {
90  _parametrized_functions.push_back(pf);
91 }
92 
94 {
95  return cast_int<unsigned int>
96  (_parametrized_functions.size());
97 }
98 
100 {
102 }
103 
105  const Point & p,
106  const Elem & elem)
107 {
108  if (var_index >= get_n_parametrized_functions())
109  libmesh_error_msg("Error: We must have var_index < get_n_parametrized_functions() in evaluate_parametrized_function.");
110 
111  return _parametrized_functions[var_index]->evaluate(get_parameters(), p, elem);
112 }
113 
115 {
116  // Short-circuit if we are using the same parameters and value of N
117  if ((_previous_parameters == get_parameters()) && (_previous_N == N))
118  {
119  return _previous_error_bound;
120  }
121 
122  // Otherwise, update _previous parameters, _previous_N
124  _previous_N = N;
125 
126  LOG_SCOPE("rb_solve()", "RBEIMEvaluation");
127 
128  if (N > get_n_basis_functions())
129  libmesh_error_msg("ERROR: N cannot be larger than the number of basis functions in rb_solve");
130 
131  if (N==0)
132  libmesh_error_msg("ERROR: N must be greater than 0 in rb_solve");
133 
134  // Get the rhs by sampling parametrized_function
135  // at the first N interpolation_points
136  DenseVector<Number> EIM_rhs(N);
137  for (unsigned int i=0; i<N; i++)
138  {
142  }
143 
144 
145 
146  DenseMatrix<Number> interpolation_matrix_N;
147  interpolation_matrix.get_principal_submatrix(N, interpolation_matrix_N);
148 
149  interpolation_matrix_N.lu_solve(EIM_rhs, RB_solution);
150 
151  // Optionally evaluate an a posteriori error bound. The EIM error estimate
152  // recommended in the literature is based on using "next" EIM point, so
153  // we skip this if N == get_n_basis_functions()
155  {
156  // Compute the a posteriori error bound
157  // First, sample the parametrized function at x_{N+1}
161 
162  // Next, evaluate the EIM approximation at x_{N+1}
163  Number EIM_approx_at_next_x = 0.;
164  for (unsigned int j=0; j<N; j++)
165  {
166  EIM_approx_at_next_x += RB_solution(j) * interpolation_matrix(N,j);
167  }
168 
169  Real error_estimate = std::abs(g_at_next_x - EIM_approx_at_next_x);
170 
171  _previous_error_bound = error_estimate;
172  return error_estimate;
173  }
174  else // Don't evaluate an error bound
175  {
176  _previous_error_bound = -1.;
177  return -1.;
178  }
179 
180 }
181 
183 {
184  LOG_SCOPE("rb_solve()", "RBEIMEvaluation");
185 
186  if (EIM_rhs.size() > get_n_basis_functions())
187  libmesh_error_msg("ERROR: N cannot be larger than the number of basis functions in rb_solve");
188 
189  if (EIM_rhs.size()==0)
190  libmesh_error_msg("ERROR: N must be greater than 0 in rb_solve");
191 
192  const unsigned int N = EIM_rhs.size();
193  DenseMatrix<Number> interpolation_matrix_N;
194  interpolation_matrix.get_principal_submatrix(N, interpolation_matrix_N);
195 
196  interpolation_matrix_N.lu_solve(EIM_rhs, RB_solution);
197 }
198 
200 {
201  // Just set the normalization factor to 1 in this case.
202  // Users can override this method if specific behavior
203  // is required.
204 
205  return 1.;
206 }
207 
209 {
210  // Initialize the rb_theta objects that access the solution from this rb_eim_evaluation
211  _rb_eim_theta_objects.clear();
212  for (unsigned int i=0; i<get_n_basis_functions(); i++)
213  {
214  _rb_eim_theta_objects.push_back( build_eim_theta(i).release() );
215  }
216 }
217 
219 {
220  return _rb_eim_theta_objects;
221 }
222 
224 {
225  return UniquePtr<RBTheta>( new RBEIMTheta(*this, index) );
226 }
227 
228 void RBEIMEvaluation::legacy_write_offline_data_to_files(const std::string & directory_name,
229  const bool read_binary_data)
230 {
231  LOG_SCOPE("legacy_write_offline_data_to_files()", "RBEIMEvaluation");
232 
234 
235  // Get the number of basis functions
236  unsigned int n_bfs = get_n_basis_functions();
237 
238  // The writing mode: ENCODE for binary, WRITE for ASCII
239  XdrMODE mode = read_binary_data ? ENCODE : WRITE;
240 
241  // The suffix to use for all the files that are written out
242  const std::string suffix = read_binary_data ? ".xdr" : ".dat";
243 
244  if (this->processor_id() == 0)
245  {
246  std::ostringstream file_name;
247 
248  // Next write out the interpolation_matrix
249  file_name.str("");
250  file_name << directory_name << "/interpolation_matrix" << suffix;
251  Xdr interpolation_matrix_out(file_name.str(), mode);
252 
253  for (unsigned int i=0; i<n_bfs; i++)
254  {
255  for (unsigned int j=0; j<=i; j++)
256  {
257  interpolation_matrix_out << interpolation_matrix(i,j);
258  }
259  }
260 
261  // Next write out interpolation_points
262  file_name.str("");
263  file_name << directory_name << "/interpolation_points" << suffix;
264  Xdr interpolation_points_out(file_name.str(), mode);
265 
266  for (unsigned int i=0; i<n_bfs; i++)
267  {
268  interpolation_points_out << interpolation_points[i](0);
269 
270  if (LIBMESH_DIM >= 2)
271  interpolation_points_out << interpolation_points[i](1);
272 
273  if (LIBMESH_DIM >= 3)
274  interpolation_points_out << interpolation_points[i](2);
275  }
276  interpolation_points_out.close();
277 
278  // Next write out interpolation_points_var
279  file_name.str("");
280  file_name << directory_name << "/interpolation_points_var" << suffix;
281  Xdr interpolation_points_var_out(file_name.str(), mode);
282 
283  for (unsigned int i=0; i<n_bfs; i++)
284  {
285  interpolation_points_var_out << interpolation_points_var[i];
286  }
287  interpolation_points_var_out.close();
288  }
289 
290  // Write out the elements associated with the interpolation points.
291  // This uses mesh I/O, hence we have to do it on all processors.
293 }
294 
295 void RBEIMEvaluation::legacy_write_out_interpolation_points_elem(const std::string & directory_name)
296 {
298 
299  // Maintain a set of node IDs to make sure we don't insert
300  // the same node into _interpolation_points_mesh more than once
301  std::set<dof_id_type> node_ids;
302  std::map<dof_id_type, dof_id_type> node_id_map;
303 
304  unsigned int new_node_id = 0;
305  for (std::size_t i=0; i<interpolation_points_elem.size(); i++)
306  {
307  Elem * old_elem = interpolation_points_elem[i];
308 
309  for (unsigned int n=0; n<old_elem->n_nodes(); n++)
310  {
311  Node & node_ref = old_elem->node_ref(n);
312  dof_id_type old_node_id = node_ref.id();
313 
314  // Check if this node has already been added. This
315  // could happen if some of the elements are neighbors.
316  if (node_ids.find(old_node_id) == node_ids.end())
317  {
318  node_ids.insert(old_node_id);
319  _interpolation_points_mesh.add_point(node_ref, new_node_id, /* proc_id */ 0);
320 
321  node_id_map[old_node_id] = new_node_id;
322 
323  new_node_id++;
324  }
325  }
326  }
327 
328  // Maintain a map of elem IDs to make sure we don't insert
329  // the same elem into _interpolation_points_mesh more than once
330  std::map<dof_id_type,dof_id_type> elem_id_map;
331  std::vector<dof_id_type> interpolation_elem_ids(interpolation_points_elem.size());
332  dof_id_type new_elem_id = 0;
333  for (std::size_t i=0; i<interpolation_elem_ids.size(); i++)
334  {
335  Elem * old_elem = interpolation_points_elem[i];
336 
337  dof_id_type old_elem_id = old_elem->id();
338 
339  // Only insert the element into the mesh if it hasn't already been inserted
340  std::map<dof_id_type,dof_id_type>::iterator id_it = elem_id_map.find(old_elem_id);
341  if (id_it == elem_id_map.end())
342  {
343  Elem * new_elem = Elem::build(old_elem->type(), /*parent*/ libmesh_nullptr).release();
344  new_elem->subdomain_id() = old_elem->subdomain_id();
345 
346  // Assign all the nodes
347  for (unsigned int n=0; n<new_elem->n_nodes(); n++)
348  {
349  dof_id_type old_node_id = old_elem->node_id(n);
350  new_elem->set_node(n) =
351  _interpolation_points_mesh.node_ptr( node_id_map[old_node_id] );
352  }
353 
354  // Just set all proc_ids to 0
355  new_elem->processor_id() = 0;
356 
357  // Add the element to the mesh
359 
360  // Set the id of new_elem appropriately
361  new_elem->set_id(new_elem_id);
362  interpolation_elem_ids[i] = new_elem->id();
363  elem_id_map[old_elem_id] = new_elem->id();
364 
365  new_elem_id++;
366  }
367  else
368  {
369  interpolation_elem_ids[i] = id_it->second;
370  }
371 
372  }
373 
375 
376  _interpolation_points_mesh.write(directory_name + "/interpolation_points_mesh.xda");
377 
378  // Also, write out the vector that tells us which element each entry
379  // of interpolation_points_elem corresponds to. This allows us to handle
380  // the case in which elements are repeated in interpolation_points_elem.
381  if (processor_id() == 0)
382  {
383  // These are just integers, so no need for a binary format here
384  std::ofstream interpolation_elem_ids_out
385  ((directory_name + "/interpolation_elem_ids.dat").c_str(), std::ofstream::out);
386 
387  for (std::size_t i=0; i<interpolation_elem_ids.size(); i++)
388  interpolation_elem_ids_out << interpolation_elem_ids[i] << std::endl;
389 
390  interpolation_elem_ids_out.close();
391  }
392 }
393 
394 void RBEIMEvaluation::legacy_read_offline_data_from_files(const std::string & directory_name,
395  bool read_error_bound_data,
396  const bool read_binary_data)
397 {
398  LOG_SCOPE("legacy_read_offline_data_from_files()", "RBEIMEvaluation");
399 
400  Parent::legacy_read_offline_data_from_files(directory_name, read_error_bound_data);
401 
402  // First, find out how many basis functions we had when Greedy terminated
403  // This was set in RBSystem::read_offline_data_from_files
404  unsigned int n_bfs = this->get_n_basis_functions();
405 
406  // The writing mode: DECODE for binary, READ for ASCII
407  XdrMODE mode = read_binary_data ? DECODE : READ;
408 
409  // The suffix to use for all the files that are written out
410  const std::string suffix = read_binary_data ? ".xdr" : ".dat";
411 
412  // Stream for creating file names
413  std::ostringstream file_name;
414 
415  // Read in the interpolation matrix
416  file_name.str("");
417  file_name << directory_name << "/interpolation_matrix" << suffix;
418  assert_file_exists(file_name.str());
419 
420  Xdr interpolation_matrix_in(file_name.str(), mode);
421 
422  for (unsigned int i=0; i<n_bfs; i++)
423  {
424  for (unsigned int j=0; j<=i; j++)
425  {
426  Number value;
427  interpolation_matrix_in >> value;
429  }
430  }
431  interpolation_matrix_in.close();
432 
433  // Next read in interpolation_points
434  file_name.str("");
435  file_name << directory_name << "/interpolation_points" << suffix;
436  assert_file_exists(file_name.str());
437 
438  Xdr interpolation_points_in(file_name.str(), mode);
439 
440  for (unsigned int i=0; i<n_bfs; i++)
441  {
442  Real x_val, y_val, z_val = 0.;
443  interpolation_points_in >> x_val;
444 
445  if (LIBMESH_DIM >= 2)
446  interpolation_points_in >> y_val;
447 
448  if (LIBMESH_DIM >= 3)
449  interpolation_points_in >> z_val;
450 
451  Point p(x_val, y_val, z_val);
452  interpolation_points.push_back(p);
453  }
454  interpolation_points_in.close();
455 
456  // Next read in interpolation_points_var
457  file_name.str("");
458  file_name << directory_name << "/interpolation_points_var" << suffix;
459  assert_file_exists(file_name.str());
460 
461  Xdr interpolation_points_var_in(file_name.str(), mode);
462 
463  for (unsigned int i=0; i<n_bfs; i++)
464  {
465  unsigned int var;
466  interpolation_points_var_in >> var;
467  interpolation_points_var.push_back(var);
468  }
469  interpolation_points_var_in.close();
470 
471  // Read in the elements corresponding to the interpolation points
473 }
474 
475 void RBEIMEvaluation::legacy_read_in_interpolation_points_elem(const std::string & directory_name)
476 {
477  _interpolation_points_mesh.read(directory_name + "/interpolation_points_mesh.xda");
478 
479  // We have an element for each EIM basis function
480  unsigned int n_bfs = this->get_n_basis_functions();
481 
482  std::vector<dof_id_type> interpolation_elem_ids;
483  {
484  // These are just integers, so no need for a binary format here
485  std::ifstream interpolation_elem_ids_in
486  ((directory_name + "/interpolation_elem_ids.dat").c_str(), std::ifstream::in);
487 
488  if (!interpolation_elem_ids_in)
489  libmesh_error_msg("RB data missing: " + directory_name + "/interpolation_elem_ids.dat");
490 
491  for (unsigned int i=0; i<n_bfs; i++)
492  {
493  dof_id_type elem_id;
494  interpolation_elem_ids_in >> elem_id;
495  interpolation_elem_ids.push_back(elem_id);
496  }
497  interpolation_elem_ids_in.close();
498  }
499 
500  interpolation_points_elem.resize(n_bfs);
501  for (unsigned int i=0; i<n_bfs; i++)
502  {
504  _interpolation_points_mesh.elem_ptr(interpolation_elem_ids[i]);
505  }
506 }
507 
508 }
DenseMatrix< Number > interpolation_matrix
Dense matrix that stores the lower triangular interpolation matrix that can be used.
virtual void legacy_write_offline_data_to_files(const std::string &directory_name="offline_data", const bool write_binary_data=true) libmesh_override
Write out all the data to text files in order to segregate the Offline stage from the Online stage...
virtual Node * add_point(const Point &p, const dof_id_type id=DofObject::invalid_id, const processor_id_type proc_id=DofObject::invalid_processor_id) libmesh_override
functions for adding /deleting nodes elements.
bool compute_RB_inner_product
Boolean flag to indicate whether we compute the RB_inner_product_matrix.
double abs(double a)
virtual const Elem * elem_ptr(const dof_id_type i) const libmesh_override
virtual const Node * node_ptr(const dof_id_type i) const libmesh_override
The ReplicatedMesh class is derived from the MeshBase class, and is used to store identical copies of...
virtual void read(const std::string &name, void *mesh_data=libmesh_nullptr, bool skip_renumber_nodes_and_elements=false, bool skip_find_neighbors=false) libmesh_override
Reads the file specified by name.
virtual void legacy_read_offline_data_from_files(const std::string &directory_name="offline_data", bool read_error_bound_data=true, const bool read_binary_data=true) libmesh_override
Read in the saved Offline reduced basis data to initialize the system for Online solves.
virtual Node *& set_node(const unsigned int i)
Definition: elem.h:1941
A Node is like a Point, but with more information.
Definition: node.h:52
Encapsulates the MPI_Comm object.
Definition: parallel.h:657
bool evaluate_RB_error_bound
Boolean to indicate whether we evaluate a posteriori error bounds when rb_solve is called...
static UniquePtr< Elem > build(const ElemType type, Elem *p=libmesh_nullptr)
Definition: elem.C:238
virtual ElemType type() const =0
std::vector< RBTheta * > get_eim_theta_objects()
virtual void resize_data_structures(const unsigned int Nmax, bool resize_error_bound_data=true)
Resize and clear the data vectors corresponding to the value of Nmax.
virtual Elem * add_elem(Elem *e) libmesh_override
Add elem e to the end of the element array.
std::vector< RBTheta * > _rb_eim_theta_objects
The vector of RBTheta objects that are created to point to this RBEIMEvaluation.
void allow_renumbering(bool allow)
If false is passed in then this mesh will no longer be renumbered when being prepared for use...
Definition: mesh_base.h:749
DenseVector< Number > RB_solution
The RB solution vector.
virtual void resize_data_structures(const unsigned int Nmax, bool resize_error_bound_data=true) libmesh_override
Resize the data structures for storing data associated with this object.
void legacy_read_in_interpolation_points_elem(const std::string &directory_name)
Read int interpolation_points_elem from a mesh.
A simple functor class that provides a RBParameter-dependent function.
This class provides functionality required to define an RBTheta object that arises from an "Empirical...
Definition: rb_eim_theta.h:42
ReplicatedMesh & get_interpolation_points_mesh()
Get a writable reference to the interpolation points mesh.
virtual void clear() libmesh_override
Clear all internal data.
This is the base class from which all geometric element types are derived.
Definition: elem.h:89
const class libmesh_nullptr_t libmesh_nullptr
virtual void clear() libmesh_override
Clear this object.
The libMesh namespace provides an interface to certain functionality in the library.
ReplicatedMesh _interpolation_points_mesh
Mesh object that we use to store copies of the elements associated with interpolation points...
void legacy_write_out_interpolation_points_elem(const std::string &directory_name)
Write out interpolation_points_elem by putting the elements into a mesh and writing out the mesh...
virtual unsigned int size() const libmesh_override
Definition: dense_vector.h:87
void initialize_eim_theta_objects()
Build a vector of RBTheta objects that accesses the components of the RB_solution member variable of ...
dof_id_type & set_id()
Definition: dof_object.h:641
libmesh_assert(j)
std::unique_ptr< T > UniquePtr
Definition: auto_ptr.h:46
virtual dof_id_type n_elem() const libmesh_override
virtual unsigned int n_nodes() const =0
virtual UniquePtr< RBTheta > build_eim_theta(unsigned int index)
Build a theta object corresponding to EIM index index.
virtual Real get_error_bound_normalization() libmesh_override
virtual Real rb_solve(unsigned int N) libmesh_override
Calculate the EIM approximation to parametrized_function using the first N EIM basis functions...
XdrMODE
Defines an enum for read/write mode in Xdr format.
Definition: enum_xdr_mode.h:32
virtual void legacy_write_offline_data_to_files(const std::string &directory_name="offline_data", const bool write_binary_data=true)
Write out all the data to text files in order to segregate the Offline stage from the Online stage...
virtual void legacy_read_offline_data_from_files(const std::string &directory_name="offline_data", bool read_error_bound_data=true, const bool read_binary_data=true)
Read in the saved Offline reduced basis data to initialize the system for Online solves.
void get_principal_submatrix(unsigned int sub_m, unsigned int sub_n, DenseMatrix< T > &dest) const
Put the sub_m x sub_n principal submatrix into dest.
Definition: dense_matrix.C:553
virtual void write(const std::string &name) libmesh_override
Write the file specified by name.
const Node & node_ref(const unsigned int i) const
Definition: elem.h:1896
subdomain_id_type subdomain_id() const
Definition: elem.h:1951
virtual ~RBEIMEvaluation()
Destructor.
std::vector< unsigned int > interpolation_points_var
The corresponding list of variables indices at which the interpolation points were identified...
std::vector< Elem * > interpolation_points_elem
The corresponding list of elements at which the interpolation points were identified.
void set_rb_theta_expansion(RBThetaExpansion &rb_theta_expansion_in)
Set the RBThetaExpansion object.
Definition: rb_evaluation.C:84
This class implements a C++ interface to the XDR (eXternal Data Representation) format.
Definition: xdr_cxx.h:68
unsigned int _previous_N
Store the number of basis functions used for the previous solve (so we can avoid an unnecessary repea...
This class is part of the rbOOmit framework.
Definition: rb_evaluation.h:50
Number evaluate_parametrized_function(unsigned int var_index, const Point &p, const Elem &elem)
RBEIMEvaluation(const libMesh::Parallel::Communicator &comm_in LIBMESH_CAN_DEFAULT_TO_COMMWORLD)
Constructor.
virtual void clear() libmesh_override
Clear this RBEvaluation object.
Definition: rb_evaluation.C:60
void lu_solve(const DenseVector< T > &b, DenseVector< T > &x)
Solve the system Ax=b given the input vector b.
Definition: dense_matrix.C:589
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
void assert_file_exists(const std::string &file_name)
Helper function that checks if file_name exists.
OStreamProxy out
static const bool value
Definition: xdr_io.C:108
void resize(const unsigned int new_m, const unsigned int new_n)
Resize the matrix.
Definition: dense_matrix.h:776
std::vector< RBParametrizedFunction * > _parametrized_functions
This vector stores the parametrized functions that will be approximated in this EIM system...
RBThetaExpansion _empty_rb_theta_expansion
We initialize RBEIMEvaluation so that it has an "empty" RBThetaExpansion, because this isn&#39;t used at ...
dof_id_type node_id(const unsigned int i) const
Definition: elem.h:1831
const RBParameters & get_parameters() const
Get the current parameters.
virtual unsigned int get_n_basis_functions() const
Get the current number of basis functions.
dof_id_type id() const
Definition: dof_object.h:632
std::vector< Point > interpolation_points
The list of interpolation points, i.e.
Real _previous_error_bound
Store the previous error bound returned by rb_solve (so we can return it if we are avoiding an unnece...
RBParameters _previous_parameters
Store the parameters at which the previous solve was performed (so we can avoid an unnecessary repeat...
void attach_parametrized_function(RBParametrizedFunction *pf)
Attach the parametrized function that we will approximate using the Empirical Interpolation Method...
A Point defines a location in LIBMESH_DIM dimensional Real space.
Definition: point.h:38
uint8_t dof_id_type
Definition: id_types.h:64
processor_id_type processor_id() const
processor_id_type processor_id() const
Definition: dof_object.h:694
unsigned int get_n_parametrized_functions() const
Get the number of parametrized functions that have been attached to this system.