libMesh
reduced_basis_ex4.C
Go to the documentation of this file.
1 // The libMesh Finite Element Library.
2 // Copyright (C) 2002-2017 Benjamin S. Kirk, John W. Peterson, Roy H. Stogner
3 
4 // This library is free software; you can redistribute it and/or
5 // modify it under the terms of the GNU Lesser General Public
6 // License as published by the Free Software Foundation; either
7 // version 2.1 of the License, or (at your option) any later version.
8 
9 // This library is distributed in the hope that it will be useful,
10 // but WITHOUT ANY WARRANTY; without even the implied warranty of
11 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
12 // Lesser General Public License for more details.
13 
14 // You should have received a copy of the GNU Lesser General Public
15 // License along with this library; if not, write to the Free Software
16 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
17 
18 // rbOOmit: An implementation of the Certified Reduced Basis method.
19 // Copyright (C) 2009, 2010 David J. Knezevic
20 // This file is part of rbOOmit.
21 
22 
23 
24 // <h1>Reduced Basis Example 4 - Empirical Interpolation Method</h1>
25 // \author David Knezevic
26 // \date 2011
27 //
28 // In this example problem we develop a reduced basis approximation
29 // for a parametrized PDE that has "non-affine" parameter
30 // dependence. This requires the use of the Empirical Interpolation
31 // Method (EIM).
32 //
33 // We first use EIM to construct an affine approximation to the
34 // non-affine term, which is a parametrized function that is a
35 // Gaussian with "center" defined by the two parameters (mu_1, mu_2)
36 // \in [-1,1]^2. We then employ this EIM approximation in order to
37 // generate a reduced basis approximation for the parametrized PDE:
38 // -0.05 * Laplacian(u) = f(mu_1, mu_2), with zero Dirichlet boundary
39 // conditions.
40 
41 // Basic include file needed for the mesh functionality.
42 #include "libmesh/libmesh.h"
43 #include "libmesh/mesh.h"
44 #include "libmesh/mesh_generation.h"
45 #include "libmesh/equation_systems.h"
46 #include "libmesh/exodusII_io.h"
47 #include "libmesh/getpot.h"
48 #include "libmesh/rb_data_serialization.h"
49 #include "libmesh/rb_data_deserialization.h"
50 
51 #include "eim_classes.h"
52 #include "rb_classes.h"
53 
54 using namespace libMesh;
55 
56 
57 int main (int argc, char ** argv)
58 {
59  // Initialize libMesh.
60  LibMeshInit init (argc, argv);
61 
62  // This example requires a linear solver package.
63  libmesh_example_requires(libMesh::default_solver_package() != INVALID_SOLVER_PACKAGE,
64  "--enable-petsc, --enable-trilinos, or --enable-eigen");
65 
66 #if !defined(LIBMESH_HAVE_XDR)
67  // We need XDR support to write out reduced bases
68  libmesh_example_requires(false, "--enable-xdr");
69 #elif defined(LIBMESH_DEFAULT_SINGLE_PRECISION)
70  // XDR binary support requires double precision
71  libmesh_example_requires(false, "double precision");
72 #elif defined(LIBMESH_DEFAULT_TRIPLE_PRECISION)
73  // I have no idea why long double isn't working here... [RHS]
74  libmesh_example_requires(false, "double precision");
75 #endif
76 
77  // Skip this 2D example if libMesh was compiled as 1D-only.
78  libmesh_example_requires(2 <= LIBMESH_DIM, "2D support");
79 
80  // Define the names of the input files we will read the problem properties from
81  std::string eim_parameters = "eim.in";
82  std::string rb_parameters = "rb.in";
83  std::string main_parameters = "reduced_basis_ex4.in";
84  GetPot infile(main_parameters);
85 
86  unsigned int n_elem = infile("n_elem", 1); // Determines the number of elements in the "truth" mesh
87  const unsigned int dim = 2; // The number of spatial dimensions
88  bool store_basis_functions = infile("store_basis_functions", false); // Do we write out basis functions?
89 
90  // Read the "online_mode" flag from the command line
91  GetPot command_line (argc, argv);
92  int online_mode = 0;
93  if (command_line.search(1, "-online_mode"))
94  online_mode = command_line.next(online_mode);
95 
96  // Create a mesh (just a simple square) on the default MPI
97  // communicator. We currently have to create a ReplicatedMesh here
98  // due to a reduced_basis regression with DistributedMesh
99  ReplicatedMesh mesh (init.comm(), dim);
101  n_elem, n_elem,
102  -1., 1.,
103  -1., 1.,
104  QUAD4);
105 
106  // Initialize the EquationSystems object for this mesh and attach
107  // the EIM and RB Construction objects
108  EquationSystems equation_systems (mesh);
109 
110  SimpleEIMConstruction & eim_construction =
111  equation_systems.add_system<SimpleEIMConstruction> ("EIM");
112  SimpleRBConstruction & rb_construction =
113  equation_systems.add_system<SimpleRBConstruction> ("RB");
114 
115  // Initialize the data structures for the equation system.
116  equation_systems.init ();
117 
118  // Print out some information about the "truth" discretization
119  mesh.print_info();
120  equation_systems.print_info();
121 
122  // Initialize the standard RBEvaluation object
123  SimpleRBEvaluation rb_eval(mesh.comm());
124 
125  // Initialize the EIM RBEvaluation object
126  SimpleEIMEvaluation eim_rb_eval(mesh.comm());
127 
128  // Set the rb_eval objects for the RBConstructions
129  eim_construction.set_rb_evaluation(eim_rb_eval);
130  rb_construction.set_rb_evaluation(rb_eval);
131 
132  if (!online_mode)
133  {
134  // Read data from input file and print state
135  eim_construction.process_parameters_file(eim_parameters);
136  eim_construction.print_info();
137 
138  // Perform the EIM Greedy and write out the data
139  eim_construction.initialize_rb_construction();
140  eim_construction.train_reduced_basis();
141 
142 #if defined(LIBMESH_HAVE_CAPNPROTO)
143  RBDataSerialization::RBEIMEvaluationSerialization rb_eim_eval_writer(eim_rb_eval);
144  rb_eim_eval_writer.write_to_file("rb_eim_eval.bin");
145 #else
146  eim_construction.get_rb_evaluation().legacy_write_offline_data_to_files("eim_data");
147 #endif
148 
149  // Read data from input file and print state
150  rb_construction.process_parameters_file(rb_parameters);
151 
152  // attach the EIM theta objects to the RBConstruction and RBEvaluation objects
153  eim_rb_eval.initialize_eim_theta_objects();
154  rb_eval.get_rb_theta_expansion().attach_multiple_F_theta(eim_rb_eval.get_eim_theta_objects());
155 
156  // attach the EIM assembly objects to the RBConstruction object
157  eim_construction.initialize_eim_assembly_objects();
158  rb_construction.get_rb_assembly_expansion().attach_multiple_F_assembly(eim_construction.get_eim_assembly_objects());
159 
160  // Print out the state of rb_construction now that the EIM objects have been attached
161  rb_construction.print_info();
162 
163  // Need to initialize _after_ EIM greedy so that
164  // the system knows how many affine terms there are
165  rb_construction.initialize_rb_construction();
166  rb_construction.train_reduced_basis();
167 
168 #if defined(LIBMESH_HAVE_CAPNPROTO)
169  RBDataSerialization::RBEvaluationSerialization rb_eval_writer(rb_construction.get_rb_evaluation());
170  rb_eval_writer.write_to_file("rb_eval.bin");
171 #else
172  rb_construction.get_rb_evaluation().legacy_write_offline_data_to_files("rb_data");
173 #endif
174 
175  // Write out the basis functions, if requested
176  if (store_basis_functions)
177  {
178  // Write out the basis functions
179  eim_construction.get_rb_evaluation().write_out_basis_functions(eim_construction.get_explicit_system(),
180  "eim_data");
181 
182  rb_construction.get_rb_evaluation().write_out_basis_functions(rb_construction,
183  "rb_data");
184  }
185  }
186  else
187  {
188 #if defined(LIBMESH_HAVE_CAPNPROTO)
189  RBDataDeserialization::RBEIMEvaluationDeserialization rb_eim_eval_reader(eim_rb_eval);
190  rb_eim_eval_reader.read_from_file("rb_eim_eval.bin");
191 #else
192  eim_rb_eval.legacy_read_offline_data_from_files("eim_data");
193 #endif
194 
195  // attach the EIM theta objects to rb_eval objects
196  eim_rb_eval.initialize_eim_theta_objects();
197  rb_eval.get_rb_theta_expansion().attach_multiple_F_theta(eim_rb_eval.get_eim_theta_objects());
198 
199  // Read in the offline data for rb_eval
200 #if defined(LIBMESH_HAVE_CAPNPROTO)
202  rb_eval_reader.read_from_file("rb_eval.bin", /*read_error_bound_data*/ true);
203 #else
204  rb_eval.legacy_read_offline_data_from_files("rb_data");
205 #endif
206 
207  // Get the parameters at which we will do a reduced basis solve
208  Real online_center_x = infile("online_center_x", 0.);
209  Real online_center_y = infile("online_center_y", 0.);
210  RBParameters online_mu;
211  online_mu.set_value("center_x", online_center_x);
212  online_mu.set_value("center_y", online_center_y);
213  rb_eval.set_parameters(online_mu);
214  rb_eval.print_parameters();
215  rb_eval.rb_solve(rb_eval.get_n_basis_functions());
216 
217  // plot the solution, if requested
218  if (store_basis_functions)
219  {
220  // read in the data from files
221  eim_rb_eval.read_in_basis_functions(eim_construction.get_explicit_system(), "eim_data");
222  rb_eval.read_in_basis_functions(rb_construction, "rb_data");
223 
224  eim_construction.load_rb_solution();
225  rb_construction.load_rb_solution();
226 #ifdef LIBMESH_HAVE_EXODUS_API
227  ExodusII_IO(mesh).write_equation_systems("RB_sol.e", equation_systems);
228 #endif
229  }
230  }
231 }
This is the EquationSystems class.
The ReplicatedMesh class is derived from the MeshBase class, and is used to store identical copies of...
int main(int argc, char **argv)
This class serializes an RBEvaluation object using the Cap&#39;n Proto library.
dof_id_type n_elem(const MeshBase::const_element_iterator &begin, const MeshBase::const_element_iterator &end)
Count up the number of elements of a specific type (as defined by an iterator range).
Definition: mesh_tools.C:656
unsigned int dim
The ExodusII_IO class implements reading meshes in the ExodusII file format from Sandia National Labs...
Definition: exodusII_io.h:52
MeshBase & mesh
void build_square(UnstructuredMesh &mesh, const unsigned int nx, const unsigned int ny, const Real xmin=0., const Real xmax=1., const Real ymin=0., const Real ymax=1., const ElemType type=INVALID_ELEM, const bool gauss_lobatto_grid=false)
A specialized build_cube() for 2D meshes.
The LibMeshInit class, when constructed, initializes the dependent libraries (e.g.
Definition: libmesh.h:62
The libMesh namespace provides an interface to certain functionality in the library.
This class serializes an RBEIMEvaluation object using the Cap&#39;n Proto library.
SolverPackage default_solver_package()
Definition: libmesh.C:995
void print_info(std::ostream &os=libMesh::out) const
Prints information about the equation systems, by default to libMesh::out.
void init(triangulateio &t)
Initializes the fields of t to NULL/0 as necessary.
virtual System & add_system(const std::string &system_type, const std::string &name)
Add the system of type system_type named name to the systems array.
void read_from_file(const std::string &path, bool read_error_bound_data)
Write the Cap&#39;n&#39;Proto buffer to disk.
const Parallel::Communicator & comm() const
Definition: libmesh.h:81
This class is part of the rbOOmit framework.
Definition: rb_parameters.h:42
void read_from_file(const std::string &path)
Write the Cap&#39;n&#39;Proto buffer to disk.
void write_to_file(const std::string &path)
Write the Cap&#39;n&#39;Proto buffer to disk.
void write_to_file(const std::string &path)
Write the Cap&#39;n&#39;Proto buffer to disk.
virtual void write_equation_systems(const std::string &, const EquationSystems &, const std::set< std::string > *system_names=libmesh_nullptr)
This method implements writing a mesh with data to a specified file where the data is taken from the ...
Definition: mesh_output.C:31
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
void set_value(const std::string &param_name, Real value)
Set the value of the specified parameter.
Definition: rb_parameters.C:57
const Parallel::Communicator & comm() const
This class de-serializes an RBEvaluation object using the Cap&#39;n Proto library.
virtual void init()
Initialize all the systems.
void print_info(std::ostream &os=libMesh::out) const
Prints relevant information about the mesh.
Definition: mesh_base.C:448
This class de-serializes a RBEIMEvaluation object using the Cap&#39;n Proto library.