libMesh
reduced_basis_ex3.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 3 - Transient Reduced Basis Problem</h1>
25 // \author David Knezevic
26 // \date 2011
27 //
28 // In this example problem we use the Certified Reduced Basis method
29 // to solve a transient convection-diffusion problem on the unit square.
30 // The PDE is similar to reduced_basis_ex1, except there is a time-derivative
31 // in this case.
32 
33 // Basic include file needed for the mesh functionality.
34 #include "libmesh/libmesh.h"
35 #include "libmesh/mesh.h"
36 #include "libmesh/mesh_generation.h"
37 #include "libmesh/exodusII_io.h"
38 #include "libmesh/equation_systems.h"
39 #include "libmesh/dof_map.h"
40 #include "libmesh/getpot.h"
41 #include "libmesh/elem.h"
42 #include "libmesh/rb_data_serialization.h"
43 #include "libmesh/rb_data_deserialization.h"
44 
45 // local includes
46 #include "rb_classes.h"
47 #include "assembly.h"
48 
49 // Bring in everything from the libMesh namespace
50 using namespace libMesh;
51 
52 // The main program.
53 int main (int argc, char** argv)
54 {
55  // Initialize libMesh.
56  LibMeshInit init (argc, argv);
57 
58  // This example requires SLEPc
59 #if !defined(LIBMESH_HAVE_SLEPC)
60  libmesh_example_requires(false, "--enable-slepc");
61 #else
62 
63 #if !defined(LIBMESH_HAVE_XDR)
64  // We need XDR support to write out reduced bases
65  libmesh_example_requires(false, "--enable-xdr");
66 #elif defined(LIBMESH_DEFAULT_SINGLE_PRECISION)
67  // XDR binary support requires double precision
68  libmesh_example_requires(false, "--disable-singleprecision");
69 #endif
70  // FIXME: This example currently segfaults with Trilinos?
71  libmesh_example_requires(libMesh::default_solver_package() == PETSC_SOLVERS, "--enable-petsc");
72 
73  // Skip this 2D example if libMesh was compiled as 1D-only.
74  libmesh_example_requires(2 <= LIBMESH_DIM, "2D support");
75 
76  // Parse the input file (reduced_basis_ex3.in) using GetPot
77  std::string parameters_filename = "reduced_basis_ex3.in";
78  GetPot infile(parameters_filename);
79 
80  unsigned int n_elem = infile("n_elem", 1); // Determines the number of elements in the "truth" mesh
81  const unsigned int dim = 2; // The number of spatial dimensions
82 
83  bool store_basis_functions = infile("store_basis_functions", true); // Do we write the RB basis functions to disk?
84 
85  // Read the "online_mode" flag from the command line
86  GetPot command_line (argc, argv);
87  int online_mode = 0;
88  if (command_line.search(1, "-online_mode"))
89  online_mode = command_line.next(online_mode);
90 
91  // Build a mesh on the default MPI communicator.
92  Mesh mesh (init.comm(), dim);
94  n_elem, n_elem,
95  0., 1.,
96  0., 1.,
97  QUAD4);
98 
99  // Create an equation systems object.
100  EquationSystems equation_systems (mesh);
101 
102  // We override RBConstruction with SimpleRBConstruction in order to
103  // specialize a few functions for this particular problem.
104  SimpleRBConstruction & rb_con =
105  equation_systems.add_system<SimpleRBConstruction> ("RBConvectionDiffusion");
106 
107 
108  // Initialize the data structures for the equation system.
109  equation_systems.init ();
110 
111  // Print out some information about the "truth" discretization
112  equation_systems.print_info();
113  mesh.print_info();
114 
115  // Build a new RBEvaluation object which will be used to perform
116  // Reduced Basis calculations. This is required in both the
117  // "Offline" and "Online" stages.
118  SimpleRBEvaluation rb_eval(mesh.comm());
119 
120  // Finally, we need to give the RBConstruction object a pointer to
121  // our RBEvaluation object
122  rb_con.set_rb_evaluation(rb_eval);
123 
124  if (!online_mode) // Perform the Offline stage of the RB method
125  {
126  // Read in the data that defines this problem from the specified text file
127  rb_con.process_parameters_file(parameters_filename);
128 
129  // Print out info that describes the current setup of rb_con
130  rb_con.print_info();
131 
132  // Prepare rb_con for the Construction stage of the RB method.
133  // This sets up the necessary data structures and performs
134  // initial assembly of the "truth" affine expansion of the PDE.
136 
137  // Compute the reduced basis space by computing "snapshots", i.e.
138  // "truth" solves, at well-chosen parameter values and employing
139  // these snapshots as basis functions.
140  rb_con.train_reduced_basis();
141 
142  // Write out the data that will subsequently be required for the Evaluation stage
143 #if defined(LIBMESH_HAVE_CAPNPROTO)
145  rb_eval_writer.write_to_file("trans_rb_eval.bin");
146 #else
147  rb_eval.legacy_write_offline_data_to_files();
148 #endif
149 
150  // If requested, write out the RB basis functions for visualization purposes
151  if (store_basis_functions)
152  {
153  // Write out the basis functions
155  }
156  }
157  else // Perform the Online stage of the RB method
158  {
159  // Read in the reduced basis data
160 #if defined(LIBMESH_HAVE_CAPNPROTO)
162  rb_eval_reader.read_from_file("trans_rb_eval.bin", /*read_error_bound_data*/ true);
163 #else
164  rb_eval.legacy_read_offline_data_from_files();
165 #endif
166 
167  // Read in online_N and initialize online parameters
168  unsigned int online_N = infile("online_N", 1);
169  Real online_x_vel = infile("online_x_vel", 0.);
170  Real online_y_vel = infile("online_y_vel", 0.);
171  RBParameters online_mu;
172  online_mu.set_value("x_vel", online_x_vel);
173  online_mu.set_value("y_vel", online_y_vel);
174  rb_eval.set_parameters(online_mu);
175  rb_eval.print_parameters();
176 
177  // Now do the Online solve using the precomputed reduced basis
178  Real error_bound_final_time = rb_eval.rb_solve(online_N);
179 
180  libMesh::out << "Error bound (absolute) at the final time is "
181  << error_bound_final_time << std::endl << std::endl;
182 
183  if (store_basis_functions)
184  {
185  // Read in the basis functions
186  rb_eval.read_in_basis_functions(rb_con);
187 
188  // Plot the solution at the final time level
189  rb_con.pull_temporal_discretization_data(rb_eval);
190  rb_con.set_time_step(rb_con.get_n_time_steps());
191  rb_con.load_rb_solution();
192 #ifdef LIBMESH_HAVE_EXODUS_API
193  ExodusII_IO(mesh).write_equation_systems ("RB_sol.e", equation_systems);
194 #endif
195 
196  // Plot the first basis function that was generated from the train_reduced_basis
197  // call in the Offline stage
198  rb_con.load_basis_function(0);
199 #ifdef LIBMESH_HAVE_EXODUS_API
200  ExodusII_IO(mesh).write_equation_systems ("bf0.e", equation_systems);
201 #endif
202  }
203  }
204 
205  return 0;
206 
207 #endif // LIBMESH_HAVE_SLEPC
208 }
virtual void initialize_rb_construction(bool skip_matrix_assembly=false, bool skip_vector_assembly=false)
Allocate all the data structures necessary for the construction stage of the RB method.
This is the EquationSystems class.
virtual void write_out_basis_functions(System &sys, const std::string &directory_name="offline_data", const bool write_binary_basis_functions=true)
Write out all the basis functions to file.
This class de-serializes a TransientRBEvaluation object using the Cap&#39;n Proto library.
virtual Real train_reduced_basis(const bool resize_rb_eval_data=true)
Train the reduced basis.
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
This class serializes a TransientRBEvaluation object using the Cap&#39;n Proto library.
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.
SolverPackage default_solver_package()
Definition: libmesh.C:995
virtual void load_basis_function(unsigned int i)
Load the i^th RB function into the RBConstruction solution vector.
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 void print_info()
Print out info that describes the current setup of this RBConstruction.
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.
const Parallel::Communicator & comm() const
Definition: libmesh.h:81
void write_to_file(const std::string &path)
Write the Cap&#39;n&#39;Proto buffer to disk.
This class is part of the rbOOmit framework.
Definition: rb_parameters.h:42
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
OStreamProxy out
int main(int argc, char **argv)
void read_from_file(const std::string &path, bool read_error_bound_data)
Write the Cap&#39;n&#39;Proto buffer to disk.
virtual void process_parameters_file(const std::string &parameters_filename)
Read in from the file specified by parameters_filename and set the this system&#39;s member variables acc...
RBEvaluation & get_rb_evaluation()
Get a reference to the RBEvaluation object.
virtual void init()
Initialize all the systems.
virtual void load_rb_solution()
Load the RB solution from the most recent solve with rb_eval into this system&#39;s solution vector...
void set_rb_evaluation(RBEvaluation &rb_eval_in)
Set the RBEvaluation object.
The Mesh class is a thin wrapper, around the ReplicatedMesh class by default.
Definition: mesh.h:50
void print_info(std::ostream &os=libMesh::out) const
Prints relevant information about the mesh.
Definition: mesh_base.C:448