queso-0.57.1
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Pages
LibMeshOperatorBase.C
Go to the documentation of this file.
1 //-----------------------------------------------------------------------bl-
2 //--------------------------------------------------------------------------
3 //
4 // QUESO - a library to support the Quantification of Uncertainty
5 // for Estimation, Simulation and Optimization
6 //
7 // Copyright (C) 2008-2017 The PECOS Development Team
8 //
9 // This library is free software; you can redistribute it and/or
10 // modify it under the terms of the Version 2.1 GNU Lesser General
11 // Public License as published by the Free Software Foundation.
12 //
13 // This library is distributed in the hope that it will be useful,
14 // but WITHOUT ANY WARRANTY; without even the implied warranty of
15 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
16 // Lesser General Public License for more details.
17 //
18 // You should have received a copy of the GNU Lesser General Public
19 // License along with this library; if not, write to the Free Software
20 // Foundation, Inc. 51 Franklin Street, Fifth Floor,
21 // Boston, MA 02110-1301 USA
22 //
23 //-----------------------------------------------------------------------el-
24 
25 #include <queso/Defines.h>
26 
27 #ifdef QUESO_HAVE_LIBMESH_SLEPC
28 
29 #include <iostream>
30 #include <fstream>
31 #include <vector>
32 #include <queso/FunctionOperatorBuilder.h>
33 #include <queso/LibMeshFunction.h>
34 #include <queso/LibMeshOperatorBase.h>
35 
36 #include <libmesh/libmesh.h>
37 #include <libmesh/mesh.h>
38 #include <libmesh/mesh_generation.h>
39 #include <libmesh/equation_systems.h>
40 #include <libmesh/explicit_system.h>
41 #include <libmesh/eigen_system.h>
42 #include <libmesh/condensed_eigen_system.h>
43 #include <libmesh/numeric_vector.h>
44 #include <libmesh/exodusII_io.h>
45 #include <libmesh/parallel.h>
46 #include <libmesh/parallel_implementation.h>
47 
48 namespace QUESO {
49 
51  const FunctionOperatorBuilder & builder, libMesh::MeshBase & m)
52  : OperatorBase(),
53  equation_systems(new libMesh::EquationSystems(m)),
54  builder(builder)
55 {
56 #ifndef LIBMESH_HAVE_SLEPC
57  if (m.processor_id() == 0)
58  std::cerr << "ERROR: This example requires libMesh to be\n"
59  << "compiled with SLEPc eigen solvers support!"
60  << std::endl;
61 #else
62 
63 #ifdef LIBMESH_DEFAULT_SINGLE_PRECISION
64  // SLEPc currently gives us a nasty crash with Real==float
65  libmesh_example_assert(false, "--disable-singleprecision");
66 #endif
67 
68 #ifdef LIBMESH_USE_COMPLEX_NUMBERS
69  // SLEPc currently gives us an "inner product not well defined" with
70  // Number==complex
71  libmesh_example_assert(false, "--disable-complex");
72 #endif
73 
74  // Create a CondensedEigenSystem named "Eigensystem" and (for convenience)
75  // use a reference to the system we create.
76  this->equation_systems->add_system<libMesh::CondensedEigenSystem>("Eigensystem");
77 
78 #endif // LIBMESH_HAVE_SLEPC
79 }
80 
82 {
83 }
84 
85 void LibMeshOperatorBase::save_converged_evals(const std::string & filename) const
86 {
87  unsigned int i;
88  std::ofstream evals_file(filename.c_str());
89 
90  std::pair<libMesh::Real, libMesh::Real> eval;
91  for (i = 0; i < this->nconv; i++) {
92  eval = this->equation_systems
93  ->get_system<libMesh::EigenSystem>("Eigensystem").get_eigenpair(i);
94  if (equation_systems->processor_id() == 0) {
95  evals_file << eval.first << " " << eval.second << std::endl;
96  }
97  }
98  if (this->equation_systems->processor_id() == 0) {
99  evals_file.close();
100  }
101 }
102 
103 void LibMeshOperatorBase::save_converged_evec(const std::string & filename,
104  unsigned int i) const
105 {
106  if (i < this->nconv) {
108  es->get_system<libMesh::EigenSystem>("Eigensystem").get_eigenpair(i);
109  libMesh::ExodusII_IO(es->get_mesh()).write_equation_systems(filename, *es);
110  }
111  else {
112  std::cerr << "Warning: eigenpair " << i
113  << " did not converge. Not saving."
114  << std::endl;
115  }
116 }
117 
119  return this->nconv;
120 }
121 
122 double LibMeshOperatorBase::get_eigenvalue(unsigned int i) const
123 {
124  if (i < this->nconv) {
125  std::pair<libMesh::Real, libMesh::Real> eval;
127  eval = es->get_system<libMesh::EigenSystem>("Eigensystem").get_eigenpair(i);
128  return eval.first;
129  }
130  else {
131  return -1;
132  }
133 }
134 
136 {
137  return 1.0 / this->get_eigenvalue(i);
138 }
139 
140 libMesh::EquationSystems & LibMeshOperatorBase::get_equation_systems() const
141 {
142  return *this->equation_systems;
143 }
144 
147  double alpha) const
148 {
149  unsigned int i;
151  LibMeshFunction *kl = new LibMeshFunction(this->builder, es->get_mesh());
152 
153  // Make sure all procs in libmesh mpi communicator all have the same xi. No,
154  // I can't set the seed in QUESO. That would mess with the QUESO
155  // communicator.
156  this->equation_systems->comm().broadcast(xi);
157 
159 
160  std::pair<libMesh::Real, libMesh::Real> eval;
161  for (i = 0; i < this->get_num_converged(); i++) {
162  eval = es->get_system<libMesh::EigenSystem>("Eigensystem").get_eigenpair(i);
163  kl_eq_sys->get_system<libMesh::ExplicitSystem>("Function").solution->add(
164  xi[i] / pow(eval.first, alpha / 2.0),
165  *es->get_system<libMesh::EigenSystem>("Eigensystem").solution);
166  }
167 
169  return ap;
170 }
171 
172 } // End namespace QUESO
173 
174 #endif // QUESO_HAVE_LIBMESH_SLEPC
std::shared_ptr< T > Type
Definition: SharedPtr.h:69
virtual void save_converged_evec(const std::string &filename, unsigned int i) const
Save converged eigenfunction i to filename.
virtual void save_converged_evals(const std::string &filename) const
Save the eigenvalues to file filename.
Abstract base class for operator objects. Operators are assumed to be symmetric and positive-definite...
Definition: OperatorBase.h:45
const FunctionOperatorBuilder & builder
virtual SharedPtr< libMesh::EquationSystems >::Type get_equation_systems() const
Return the internal libmesh equation systems object.
LibMeshOperatorBase(const FunctionOperatorBuilder &builder, libMesh::MeshBase &m)
Constuct an operator on the mesh m using a builder builder.
SharedPtr< libMesh::EquationSystems >::Type equation_systems
virtual double get_inverted_eigenvalue(unsigned int i) const
Return the reciprocal of eigenvalue i.
virtual unsigned int get_num_converged() const
Return the number of converged eigenpairs.
virtual double get_eigenvalue(unsigned int i) const
Return eigenvalue i.
unsigned int nconv
The number of converged eigenpairs.
virtual SharedPtr< FunctionBase >::Type inverse_kl_transform(std::vector< double > &xi, double alpha) const
Given coefficients xi, computes the Karhunen-Loeve transform.
Function objects using libMesh for the backend.
virtual libMesh::EquationSystems & get_equation_systems() const
Return the internal libmesh equation systems object.

Generated on Tue Jun 5 2018 19:48:55 for queso-0.57.1 by  doxygen 1.8.5