libMesh
reduced_basis_ex7.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 7 - Acoustic Horn</h1>
25 // \author David Knezevic
26 // \date 2012
27 //
28 // In this example problem we use the Certified Reduced Basis method
29 // to solve a complex-valued Helmholtz problem for acoustics. There is only
30 // one parameter in this problem: the forcing frequency at the horn inlet.
31 // We impose a radiation boundary condition on the outer boundary.
32 //
33 // In the Online stage we write out the reflection coefficient as a function
34 // of frequency. The reflection coefficient gives the ratio of the reflected
35 // wave to the applied wave at the inlet.
36 
37 // C++ include files that we need
38 #include <iostream>
39 #include <algorithm>
40 #include <cmath>
41 #include <set>
42 
43 // Basic include file needed for the mesh functionality.
44 #include "libmesh/libmesh.h"
45 #include "libmesh/mesh.h"
46 #include "libmesh/mesh_generation.h"
47 #include "libmesh/exodusII_io.h"
48 #include "libmesh/equation_systems.h"
49 #include "libmesh/dof_map.h"
50 #include "libmesh/getpot.h"
51 #include "libmesh/elem.h"
52 #include "libmesh/rb_data_serialization.h"
53 #include "libmesh/rb_data_deserialization.h"
54 
55 // local includes
56 #include "rb_classes.h"
57 #include "assembly.h"
58 
59 // Bring in everything from the libMesh namespace
60 using namespace libMesh;
61 
62 // The main program.
63 int main (int argc, char** argv)
64 {
65  // Initialize libMesh.
66  LibMeshInit init (argc, argv);
67 
68 #ifndef LIBMESH_USE_COMPLEX_NUMBERS
69  libmesh_example_requires(false, "--enable-complex");
70 #else
71 
72 #if !defined(LIBMESH_HAVE_XDR)
73  // We need XDR support to write out reduced bases
74  libmesh_example_requires(false, "--enable-xdr");
75 #elif defined(LIBMESH_DEFAULT_SINGLE_PRECISION)
76  // XDR binary support requires double precision
77  libmesh_example_requires(false, "--disable-singleprecision");
78 #elif defined(LIBMESH_DEFAULT_TRIPLE_PRECISION)
79  // long double doesn't work on the previous examples, don't try it
80  // here
81  libmesh_example_requires(false, "double precision");
82 #endif
83  // FIXME: This example currently segfaults with Trilinos?
84  libmesh_example_requires(libMesh::default_solver_package() == PETSC_SOLVERS, "--enable-petsc");
85 
86  // Skip this 2D example if libMesh was compiled as 1D-only.
87  libmesh_example_requires(2 <= LIBMESH_DIM, "2D support");
88 
89  // Parse the input file (reduced_basis_ex7.in) using GetPot
90  std::string parameters_filename = "reduced_basis_ex7.in";
91  GetPot infile(parameters_filename);
92 
93  const unsigned int dim = 2; // The number of spatial dimensions
94 
95  bool store_basis_functions = infile("store_basis_functions", true); // Do we write the RB basis functions to disk?
96 
97  // Read the "online_mode" flag from the command line
98  GetPot command_line (argc, argv);
99  int online_mode = 0;
100  if (command_line.search(1, "-online_mode"))
101  online_mode = command_line.next(online_mode);
102 
103  // Create a mesh on the default MPI communicator.
104  Mesh mesh(init.comm(), dim);
105  mesh.read("horn.msh");
106 
107  // Create an equation systems object.
108  EquationSystems equation_systems (mesh);
109 
110  // We override RBConstruction with SimpleRBConstruction in order to
111  // specialize a few functions for this particular problem.
112  SimpleRBConstruction & rb_con =
113  equation_systems.add_system<SimpleRBConstruction> ("Acoustics");
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  equation_systems.print_info();
120  mesh.print_info();
121 
122  // Build a new RBEvaluation object which will be used to perform
123  // Reduced Basis calculations. This is required in both the
124  // "Offline" and "Online" stages.
125  SimpleRBEvaluation rb_eval(mesh.comm());
126 
127  // We need to give the RBConstruction object a pointer to
128  // our RBEvaluation object
129  rb_con.set_rb_evaluation(rb_eval);
130 
131  if (!online_mode) // Perform the Offline stage of the RB method
132  {
133  // Read in the data that defines this problem from the specified text file
134  rb_con.process_parameters_file(parameters_filename);
135 
136  // Print out info that describes the current setup of rb_con
137  rb_con.print_info();
138 
139  // Prepare rb_con for the Construction stage of the RB method.
140  // This sets up the necessary data structures and performs
141  // initial assembly of the "truth" affine expansion of the PDE.
143 
144  // Compute the reduced basis space by computing "snapshots", i.e.
145  // "truth" solves, at well-chosen parameter values and employing
146  // these snapshots as basis functions.
147  libmesh_try
148  {
149  rb_con.train_reduced_basis();
150  }
151  libmesh_catch (LogicError & e)
152  {
153  libMesh::err << "\n\n"
154  << "********************************************************************************\n"
155  << "Training reduced basis failed, this example requires a direct solver.\n"
156  << "Try running with -ksp_type preonly -pc_type lu instead.\n"
157  << "********************************************************************************"
158  << std::endl;
159  return 77;
160  }
161 
162  // Write out the data that will subsequently be required for the Evaluation stage
163 #if defined(LIBMESH_HAVE_CAPNPROTO)
165  rb_eval_writer.write_to_file("rb_eval.bin");
166 #else
168 #endif
169 
170  // If requested, write out the RB basis functions for visualization purposes
171  if (store_basis_functions)
172  {
173  // Write out the basis functions
175  }
176 
177  // Basis functions should be orthonormal, so
178  // print out the inner products to check this
180  }
181  else // Perform the Online stage of the RB method
182  {
183  // Read in the reduced basis data
184 #if defined(LIBMESH_HAVE_CAPNPROTO)
186  rb_eval_reader.read_from_file("rb_eval.bin", /*read_error_bound_data*/ true);
187 #else
188  rb_eval.legacy_read_offline_data_from_files();
189 #endif
190 
191  // Read in online_N and initialize online parameters
192  Real online_frequency = infile("online_frequency", 0.);
193  RBParameters online_mu;
194  online_mu.set_value("frequency", online_frequency);
195  rb_eval.set_parameters(online_mu);
196  rb_eval.print_parameters();
197 
198  // Now do the Online solve using the precomputed reduced basis
199  rb_eval.rb_solve(rb_eval.get_n_basis_functions());
200 
201  if (store_basis_functions)
202  {
203  // Read in the basis functions
204  rb_eval.read_in_basis_functions(rb_con);
205 
206  // Plot the solution
207  rb_con.load_rb_solution();
208 
209  //Output the variable in ExodusII format (single precision)
210  ExodusII_IO(mesh, /*single_precision=*/true).write_equation_systems ("RB_sol_float.exo", equation_systems);
211  }
212 
213  // Now do a sweep over frequencies and write out the reflection coefficient for each frequency
214  std::ofstream reflection_coeffs_out("reflection_coefficients.dat");
215 
216  Real n_frequencies = infile("n_frequencies", 0.);
217  Real delta_f = (rb_eval.get_parameter_max("frequency") - rb_eval.get_parameter_min("frequency")) / (n_frequencies-1);
218  for (unsigned int freq_i=0; freq_i<n_frequencies; freq_i++)
219  {
220  Real frequency = rb_eval.get_parameter_min("frequency") + freq_i * delta_f;
221  online_mu.set_value("frequency", frequency);
222  rb_eval.set_parameters(online_mu);
223  rb_eval.rb_solve(rb_eval.get_n_basis_functions());
224 
225  Number complex_one(1., 0.);
226  reflection_coeffs_out << frequency << " " << std::abs(rb_eval.RB_outputs[0] - complex_one) << std::endl;
227  }
228  reflection_coeffs_out.close();
229 
230  }
231 
232 #endif
233 
234  return 0;
235 }
OStreamProxy err
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.
double abs(double a)
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.
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 dim
The ExodusII_IO class implements reading meshes in the ExodusII file format from Sandia National Labs...
Definition: exodusII_io.h:52
int main(int argc, char **argv)
MeshBase & mesh
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
void print_basis_function_orthogonality()
Print out a matrix that shows the orthogonality of the RB basis functions.
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.
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 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
A class to represent the internal "this should never happen" errors, to be thrown by "libmesh_error()...
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
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...
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
virtual void read(const std::string &name, void *mesh_data=libmesh_nullptr, bool skip_renumber_nodes_and_elements=false, bool skip_find_neighbors=false)=0
Interfaces for reading/writing a mesh to/from a file.
void print_info(std::ostream &os=libMesh::out) const
Prints relevant information about the mesh.
Definition: mesh_base.C:448