libMesh
reduced_basis_ex2.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 2 - Successive Constraint Method</h1>
25 // \author David Knezevic
26 // \date 2011
27 //
28 // In this example we extend reduced_basis_ex1 to solve a steady
29 // convection-diffusion problem on the unit square via the Reduced
30 // Basis Method. In this case, we modify the PDE so that it no longer
31 // has a parameter-independent coercivity constant. Therefore, in
32 // order to obtain an error bound, we need to employ the Successive
33 // Constraint Method (SCM) implemented in
34 // RBSCMConstruction/RBSCMEvaluation to obtain a parameter-dependent
35 // lower bound for the coercivity constant.
36 //
37 // The PDE being solved is div(k*grad(u)) + Beta*grad(u) = f
38 // k is the diffusion coefficient :
39 // - constant in the domain 0<=x<0.5 , its value is given by the first parameter mu[0]
40 // - constant in the domain 0.5<=x<=1 , its value is given by the second parameter mu[1]
41 // Beta is the convection velocity :
42 // - constant in the whole domain
43 // - equal to zero in the y-direction
44 // - its value in the x-direction is given by the third (and last) parameter mu[2]
45 // Boundary conditions :
46 // - dyu=0 on top and bottom
47 // - u=0 on the left side
48 // - dxu + Beta*u = 0 on the right side
49 
50 // C++ include files that we need
51 #include <iostream>
52 #include <algorithm>
53 #include <cstdlib> // *must* precede <cmath> for proper std:abs() on PGI, Sun Studio CC
54 #include <cmath>
55 #include <set>
56 
57 // Basic include file needed for the mesh functionality.
58 #include "libmesh/libmesh.h"
59 #include "libmesh/mesh.h"
60 #include "libmesh/mesh_generation.h"
61 #include "libmesh/exodusII_io.h"
62 #include "libmesh/equation_systems.h"
63 #include "libmesh/dof_map.h"
64 #include "libmesh/getpot.h"
65 #include "libmesh/elem.h"
66 #include "libmesh/rb_data_serialization.h"
67 #include "libmesh/rb_data_deserialization.h"
68 
69 // local includes
70 #include "rb_classes.h"
71 #include "assembly.h"
72 
73 // Bring in everything from the libMesh namespace
74 using namespace libMesh;
75 
76 // The main program.
77 int main (int argc, char ** argv)
78 {
79  // Initialize libMesh.
80  LibMeshInit init (argc, argv);
81 
82  // This example requires SLEPc and GLPK
83 #if !defined(LIBMESH_HAVE_SLEPC) || !defined(LIBMESH_HAVE_GLPK)
84  libmesh_example_requires(false, "--enable-slepc --enable-glpk");
85 #else
86 
87 #if !defined(LIBMESH_HAVE_XDR)
88  // We need XDR support to write out reduced bases
89  libmesh_example_requires(false, "--enable-xdr");
90 #elif defined(LIBMESH_DEFAULT_SINGLE_PRECISION)
91  // XDR binary support requires double precision
92  libmesh_example_requires(false, "--disable-singleprecision");
93 #endif
94  // FIXME: This example currently segfaults with Trilinos?
95  libmesh_example_requires(libMesh::default_solver_package() == PETSC_SOLVERS, "--enable-petsc");
96 
97  // FIXME: This example needs lapack, which is serial only
98  libmesh_example_requires(init.comm().size() == 1, "mpirun -np 1");
99 
100  // Skip this 2D example if libMesh was compiled as 1D-only.
101  libmesh_example_requires(2 <= LIBMESH_DIM, "2D support");
102 
103  // Parse the input file (reduced_basis_ex2.in) using GetPot
104  std::string parameters_filename = "reduced_basis_ex2.in";
105  GetPot infile(parameters_filename);
106 
107  unsigned int n_elem = infile("n_elem", 1); // Determines the number of elements in the "truth" mesh
108  const unsigned int dim = 2; // The number of spatial dimensions
109 
110  bool store_basis_functions = infile("store_basis_functions", true); // Do we write the RB basis functions to disk?
111 
112  // Read the "online_mode" flag from the command line
113  GetPot command_line (argc, argv);
114  int online_mode = 0;
115  if (command_line.search(1, "-online_mode"))
116  online_mode = command_line.next(online_mode);
117 
118  // Build a mesh on the default MPI communicator.
119  Mesh mesh (init.comm(), dim);
121  n_elem, n_elem,
122  0., 1.,
123  0., 1.,
124  QUAD4);
125 
126  // Create an equation systems object.
127  EquationSystems equation_systems (mesh);
128 
129  // We override RBConstruction with SimpleRBConstruction in order to
130  // specialize a few functions for this particular problem.
131  SimpleRBConstruction & rb_con =
132  equation_systems.add_system<SimpleRBConstruction> ("RBConvectionDiffusion");
133 
134  // Initialize the SCM Construction object
135  RBSCMConstruction & rb_scm_con =
136  equation_systems.add_system<RBSCMConstruction> ("RBSCMConvectionDiffusion");
137  rb_scm_con.set_RB_system_name("RBConvectionDiffusion");
138  rb_scm_con.add_variable("p", FIRST);
139 
140  // Initialize the data structures for the equation system.
141  equation_systems.init ();
142 
143  // Print out some information about the "truth" discretization
144  equation_systems.print_info();
145  mesh.print_info();
146 
147  // Set parameters for the eigenvalue problems that will be solved by rb_scm_con
148  equation_systems.parameters.set<unsigned int>("eigenpairs") = 1;
149  equation_systems.parameters.set<unsigned int>("basis vectors") = 3;
150  equation_systems.parameters.set<unsigned int>
151  ("linear solver maximum iterations") = 1000;
152 
153  // Build a new RBEvaluation object which will be used to perform
154  // Reduced Basis calculations. This is required in both the
155  // "Offline" and "Online" stages.
156  SimpleRBEvaluation rb_eval(mesh.comm());
157 
158  // We need to give the RBConstruction object a pointer to
159  // our RBEvaluation object
160  rb_con.set_rb_evaluation(rb_eval);
161 
162  // We also need a SCM evaluation object to perform SCM calculations
163  RBSCMEvaluation rb_scm_eval(mesh.comm());
164  rb_scm_eval.set_rb_theta_expansion(rb_eval.get_rb_theta_expansion());
165 
166  // Tell rb_eval about rb_scm_eval
167  rb_eval.rb_scm_eval = &rb_scm_eval;
168 
169  // Finally, need to give rb_scm_con and rb_eval a pointer to the
170  // SCM evaluation object, rb_scm_eval
171  rb_scm_con.set_rb_scm_evaluation(rb_scm_eval);
172 
173  if (!online_mode) // Perform the Offline stage of the RB method
174  {
175  // Read in the data that defines this problem from the specified text file
176  rb_con.process_parameters_file(parameters_filename);
177  rb_scm_con.process_parameters_file(parameters_filename);
178 
179  // Print out info that describes the current setup of rb_con
180  rb_con.print_info();
181  rb_scm_con.print_info();
182 
183  // Prepare rb_con for the Construction stage of the RB method.
184  // This sets up the necessary data structures and performs
185  // initial assembly of the "truth" affine expansion of the PDE.
187 
188  // Perform the SCM Greedy algorithm to derive the data required
189  // for rb_scm_eval to provide a coercivity lower bound.
190  rb_scm_con.perform_SCM_greedy();
191 
192  // Compute the reduced basis space by computing "snapshots", i.e.
193  // "truth" solves, at well-chosen parameter values and employing
194  // these snapshots as basis functions.
195  rb_con.train_reduced_basis();
196 
197  // Write out the data that will subsequently be required for the Evaluation stage
198 #if defined(LIBMESH_HAVE_CAPNPROTO)
200  rb_eval_writer.write_to_file("rb_eval.bin");
201  RBDataSerialization::RBSCMEvaluationSerialization rb_scm_eval_writer(rb_scm_con.get_rb_scm_evaluation());
202  rb_scm_eval_writer.write_to_file("rb_scm_eval.bin");
203 #else
205  rb_scm_con.get_rb_scm_evaluation().legacy_write_offline_data_to_files("scm_data");
206 #endif
207 
208  // If requested, write out the RB basis functions for visualization purposes
209  if (store_basis_functions)
210  {
211  // Write out the basis functions
212  rb_con.get_rb_evaluation().write_out_basis_functions(rb_con, "rb_data");
213  }
214  }
215  else // Perform the Online stage of the RB method
216  {
217 
218  // Read in the reduced basis data
219 #if defined(LIBMESH_HAVE_CAPNPROTO)
221  rb_eval_reader.read_from_file("rb_eval.bin", /*read_error_bound_data*/ true);
222  RBDataDeserialization::RBSCMEvaluationDeserialization rb_scm_eval_reader(rb_scm_eval);
223  rb_scm_eval_reader.read_from_file("rb_scm_eval.bin");
224 #else
225  rb_eval.legacy_read_offline_data_from_files("rb_data");
226  rb_scm_eval.legacy_read_offline_data_from_files("scm_data");
227 #endif
228 
229  // Read in online_N and initialize online parameters
230  unsigned int online_N = infile("online_N", 1);
231  Real online_mu_0 = infile("online_mu_0", 0.);
232  Real online_mu_1 = infile("online_mu_1", 0.);
233  Real online_mu_2 = infile("online_mu_2", 0.);
234  RBParameters online_mu;
235  online_mu.set_value("mu_0", online_mu_0);
236  online_mu.set_value("mu_1", online_mu_1);
237  online_mu.set_value("mu_2", online_mu_2);
238  rb_eval.set_parameters(online_mu);
239  rb_eval.print_parameters();
240 
241  // Now do the Online solve using the precomputed reduced basis
242  rb_eval.rb_solve(online_N);
243 
244  // Print out outputs as well as the corresponding output error bounds.
245  libMesh::out << "output 1, value = " << rb_eval.RB_outputs[0]
246  << ", bound = " << rb_eval.RB_output_error_bounds[0]
247  << std::endl;
248  libMesh::out << "output 2, value = " << rb_eval.RB_outputs[1]
249  << ", bound = " << rb_eval.RB_output_error_bounds[1]
250  << std::endl;
251  libMesh::out << "output 3, value = " << rb_eval.RB_outputs[2]
252  << ", bound = " << rb_eval.RB_output_error_bounds[2]
253  << std::endl;
254  libMesh::out << "output 4, value = " << rb_eval.RB_outputs[3]
255  << ", bound = " << rb_eval.RB_output_error_bounds[3]
256  << std::endl << std::endl;
257 
258  if (store_basis_functions)
259  {
260  // Read in the basis functions
261  rb_eval.read_in_basis_functions(rb_con, "rb_data");
262 
263  // Plot the solution
264  rb_con.load_rb_solution();
265 #ifdef LIBMESH_HAVE_EXODUS_API
266  ExodusII_IO(mesh).write_equation_systems ("RB_sol.e", equation_systems);
267 #endif
268 
269  // Plot the first basis function that was generated from the train_reduced_basis
270  // call in the Offline stage
271  rb_con.load_basis_function(0);
272 #ifdef LIBMESH_HAVE_EXODUS_API
273  ExodusII_IO(mesh).write_equation_systems ("bf0.e", equation_systems);
274 #endif
275  }
276  }
277 
278  return 0;
279 
280 #endif // LIBMESH_HAVE_SLEPC && LIBMESH_HAVE_GLPK
281 }
int main(int argc, char **argv)
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.
void set_rb_theta_expansion(RBThetaExpansion &rb_theta_expansion_in)
Set the RBThetaExpansion object.
virtual Real train_reduced_basis(const bool resize_rb_eval_data=true)
Train the reduced basis.
This class serializes an RBEvaluation object using the Cap&#39;n Proto library.
unsigned int size() const
Definition: parallel.h:726
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.
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.
This class is part of the rbOOmit framework.
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...
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.
void write_to_file(const std::string &path)
Write the Cap&#39;n&#39;Proto buffer to disk.
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
This class de-serializes a RBSCMEvaluation object using the Cap&#39;n Proto library.
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
T & set(const std::string &)
Definition: parameters.h:469
const Parallel::Communicator & comm() const
OStreamProxy out
Parameters parameters
Data structure holding arbitrary parameters.
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...
This class de-serializes an RBEvaluation object using the Cap&#39;n Proto library.
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...
This class serializes an RBSCMEvaluation object using the Cap&#39;n Proto library.
void set_rb_evaluation(RBEvaluation &rb_eval_in)
Set the RBEvaluation object.
void set_RB_system_name(const std::string &new_name)
Set the name of the associated RB system — we need this to load the (symmetrized) affine operators...
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
void read_from_file(const std::string &path)
Write the Cap&#39;n&#39;Proto buffer to disk.
This class is part of the rbOOmit framework.