libMesh
reduced_basis_ex6.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 6 - Heat transfer on a curved domain in 3D</h1>
25 // \author David Knezevic
26 // \date 2012
27 //
28 // In this example we consider heat transfer modeled by a Poisson
29 // equation with Robin boundary condition:
30 //
31 // -kappa \Laplacian u = 1, on \Omega
32 // -kappa du\dn = kappa Bi u, on \partial\Omega_Biot,
33 // u = 0 on \partial\Omega_Dirichlet,
34 //
35 // We consider a reference domain \Omega_hat =
36 // [-0.2,0.2]x[-0.2,0.2]x[0,3], and the physical domain is then obtain
37 // via the parametrized mapping:
38 //
39 // x = -1/mu + (1/mu+x_hat)*cos(mu*z_hat)
40 // y = y_hat
41 // z = (1/mu+x_hat)*sin(mu*z_hat)
42 //
43 // for (x_hat,y_hat,z_hat) \in \Omega_hat. (Here "hats" denotes
44 // reference domain.) Also, the "reference Dirichlet boundaries" are
45 // [-0.2,0.2]x[-0.2,0.2]x{0} and [-0.2,0.2]x[-0.2,0.2]x{3}, and the
46 // remaining boundaries are the "Biot" boundaries.
47 //
48 // Then, after putting the PDE into weak form and mapping it to the
49 // reference domain, we obtain:
50 // \kappa \int_\Omega_hat
51 // [ (1+mu*x_hat) v_x w_x + (1+mu*x_hat) v_y w_y + 1/(1+mu*x_hat) v_z w_z ]
52 // + \kappa Bi \int_\partial\Omega_hat_Biot1 (1-0.2mu) u v
53 // + \kappa Bi \int_\partial\Omega_hat_Biot2 (1+mu x_hat) u v
54 // + \kappa Bi \int_\partial\Omega_hat_Biot3 (1+0.2mu) u v
55 // = \int_\Omega_hat (1+mu x_hat) v
56 // where
57 // \partial\Omega_hat_Biot1 = [-0.2] x [-0.2,0.2] x [0,3]
58 // \partial\Omega_hat_Biot2 = [-0.2,0.2] x {-0.2} x [0,3] \UNION [-0.2,0.2] x {0.2} x [0,3]
59 // \partial\Omega_hat_Biot3 = [0.2] x [-0.2,0.2] x [0,3]
60 //
61 // The term
62 // \kappa \int_\Omega_hat 1/(1+mu*x_hat) v_z w_z
63 // is "non-affine" (in the Reduced Basis sense), since we can't
64 // express it in the form \sum theta_q(kappa,mu) a(v,w). As a result,
65 // (as in reduced_basis_ex4) we must employ the Empirical
66 // Interpolation Method (EIM) in order to apply the Reduced Basis
67 // method here.
68 //
69 // The approach we use is to construct an EIM approximation, G_EIM, to
70 // the vector-valued function G(x_hat,y_hat;mu) = (1 + mu*x_hat, 1 +
71 // mu*x_hat, 1/(1+mu*x_hat)) and then we express the "volumetric
72 // integral part" of the left-hand side operator as a(v,w;mu) =
73 // \int_\hat\Omega G_EIM(x_hat,y_hat;mu) \dot (v_x w_x, v_y w_y, v_z
74 // w_z). (We actually only need EIM for the third component of
75 // G_EIM, but it's helpful to demonstrate "vector-valued" EIM here.)
76 
77 // C++ include files that we need
78 #include <iostream>
79 #include <algorithm>
80 #include <cstdlib> // *must* precede <cmath> for proper std:abs() on PGI, Sun Studio CC
81 #include <cmath>
82 #include <set>
83 
84 // Basic include file needed for the mesh functionality.
85 #include "libmesh/libmesh.h"
86 #include "libmesh/mesh.h"
87 #include "libmesh/mesh_generation.h"
88 #include "libmesh/exodusII_io.h"
89 #include "libmesh/equation_systems.h"
90 #include "libmesh/dof_map.h"
91 #include "libmesh/getpot.h"
92 #include "libmesh/elem.h"
93 #include "libmesh/rb_data_serialization.h"
94 #include "libmesh/rb_data_deserialization.h"
95 
96 // local includes
97 #include "rb_classes.h"
98 #include "eim_classes.h"
99 #include "assembly.h"
100 
101 // Bring in everything from the libMesh namespace
102 using namespace libMesh;
103 
104 // Define a function to scale the mesh according to the parameter.
106  Real curvature,
107  const std::string & filename);
108 
109 // The main program.
110 int main (int argc, char ** argv)
111 {
112  // Initialize libMesh.
113  LibMeshInit init (argc, argv);
114 
115  // This example requires a linear solver package.
116  libmesh_example_requires(libMesh::default_solver_package() != INVALID_SOLVER_PACKAGE,
117  "--enable-petsc, --enable-trilinos, or --enable-eigen");
118 
119 #if !defined(LIBMESH_HAVE_XDR)
120  // We need XDR support to write out reduced bases
121  libmesh_example_requires(false, "--enable-xdr");
122 #elif defined(LIBMESH_DEFAULT_SINGLE_PRECISION)
123  // XDR binary support requires double precision
124  libmesh_example_requires(false, "--disable-singleprecision");
125 #elif defined(LIBMESH_DEFAULT_TRIPLE_PRECISION)
126  // I have no idea why long double isn't working here... [RHS]
127  libmesh_example_requires(false, "double precision");
128 #elif defined(LIBMESH_ENABLE_BLOCKED_STORAGE)
129  // This example dies with "New nonzero at (0,2) caused a malloc"
130  // when blocked storage is enabled.
131  libmesh_example_requires(false, "--disable-blocked-storage");
132 #endif
133 
134  // This is a 3D example
135  libmesh_example_requires(3 == LIBMESH_DIM, "3D support");
136 
137  // Parse the input file using GetPot
138  std::string eim_parameters = "eim.in";
139  std::string rb_parameters = "rb.in";
140  std::string main_parameters = "reduced_basis_ex6.in";
141  GetPot infile(main_parameters);
142 
143  unsigned int n_elem_xy = infile("n_elem_xy", 1);
144  unsigned int n_elem_z = infile("n_elem_z", 1);
145 
146  // Do we write the RB basis functions to disk?
147  bool store_basis_functions = infile("store_basis_functions", true);
148 
149  // Read the "online_mode" flag from the command line
150  GetPot command_line (argc, argv);
151  int online_mode = 0;
152  if (command_line.search(1, "-online_mode"))
153  online_mode = command_line.next(online_mode);
154 
155  // Create a mesh, with dimension to be overridden by build_cube, on
156  // the default MPI communicator. We currently have to create a
157  // ReplicatedMesh here due to a reduced_basis regression with
158  // DistributedMesh
159  ReplicatedMesh mesh(init.comm());
160 
162  n_elem_xy, n_elem_xy, n_elem_z,
163  -0.2, 0.2,
164  -0.2, 0.2,
165  0., 3.,
166  HEX8);
167 
168  // Create an equation systems object.
169  EquationSystems equation_systems (mesh);
170 
171  SimpleEIMConstruction & eim_construction =
172  equation_systems.add_system<SimpleEIMConstruction> ("EIM");
173  SimpleRBConstruction & rb_construction =
174  equation_systems.add_system<SimpleRBConstruction> ("RB");
175 
176  // Initialize the data structures for the equation system.
177  equation_systems.init ();
178 
179  // Print out some information about the "truth" discretization
180  equation_systems.print_info();
181  mesh.print_info();
182 
183  // Initialize the standard RBEvaluation object
184  SimpleRBEvaluation rb_eval(mesh.comm());
185 
186  // Initialize the EIM RBEvaluation object
187  SimpleEIMEvaluation eim_rb_eval(mesh.comm());
188 
189  // Set the rb_eval objects for the RBConstructions
190  eim_construction.set_rb_evaluation(eim_rb_eval);
191  rb_construction.set_rb_evaluation(rb_eval);
192 
193  if (!online_mode) // Perform the Offline stage of the RB method
194  {
195  // Read data from input file and print state
196  eim_construction.process_parameters_file(eim_parameters);
197  eim_construction.print_info();
198 
199  // Perform the EIM Greedy and write out the data
200  eim_construction.initialize_rb_construction();
201 
202  eim_construction.train_reduced_basis();
203 
204 #if defined(LIBMESH_HAVE_CAPNPROTO)
205  RBDataSerialization::RBEIMEvaluationSerialization rb_eim_eval_writer(eim_rb_eval);
206  rb_eim_eval_writer.write_to_file("rb_eim_eval.bin");
207 #else
208  eim_construction.get_rb_evaluation().legacy_write_offline_data_to_files("eim_data");
209 #endif
210 
211  // Read data from input file and print state
212  rb_construction.process_parameters_file(rb_parameters);
213 
214  // attach the EIM theta objects to the RBEvaluation
215  eim_rb_eval.initialize_eim_theta_objects();
216  rb_eval.get_rb_theta_expansion().attach_multiple_A_theta(eim_rb_eval.get_eim_theta_objects());
217 
218  // attach the EIM assembly objects to the RBConstruction
219  eim_construction.initialize_eim_assembly_objects();
220  rb_construction.get_rb_assembly_expansion().attach_multiple_A_assembly(eim_construction.get_eim_assembly_objects());
221 
222  // Print out the state of rb_construction now that the EIM objects have been attached
223  rb_construction.print_info();
224 
225  // Need to initialize _after_ EIM greedy so that
226  // the system knows how many affine terms there are
227  rb_construction.initialize_rb_construction();
228  rb_construction.train_reduced_basis();
229 
230 #if defined(LIBMESH_HAVE_CAPNPROTO)
231  RBDataSerialization::RBEvaluationSerialization rb_eval_writer(rb_construction.get_rb_evaluation());
232  rb_eval_writer.write_to_file("rb_eval.bin");
233 #else
234  rb_construction.get_rb_evaluation().legacy_write_offline_data_to_files("rb_data");
235 #endif
236 
237  // Write out the basis functions, if requested
238  if (store_basis_functions)
239  {
240  // Write out the basis functions
241  eim_construction.get_rb_evaluation().write_out_basis_functions(eim_construction.get_explicit_system(),
242  "eim_data");
243  rb_construction.get_rb_evaluation().write_out_basis_functions(rb_construction,
244  "rb_data");
245  }
246  }
247  else // Perform the Online stage of the RB method
248  {
249 #if defined(LIBMESH_HAVE_CAPNPROTO)
250  RBDataDeserialization::RBEIMEvaluationDeserialization rb_eim_eval_reader(eim_rb_eval);
251  rb_eim_eval_reader.read_from_file("rb_eim_eval.bin");
252 #else
253  eim_rb_eval.legacy_read_offline_data_from_files("eim_data");
254 #endif
255 
256  // attach the EIM theta objects to rb_eval objects
257  eim_rb_eval.initialize_eim_theta_objects();
258  rb_eval.get_rb_theta_expansion().attach_multiple_A_theta(eim_rb_eval.get_eim_theta_objects());
259 
260  // Read in the offline data for rb_eval
261 #if defined(LIBMESH_HAVE_CAPNPROTO)
263  rb_eval_reader.read_from_file("rb_eval.bin", /*read_error_bound_data*/ true);
264 #else
265  rb_eval.legacy_read_offline_data_from_files("rb_data");
266 #endif
267 
268  // Get the parameters at which we will do a reduced basis solve
269  Real online_curvature = infile("online_curvature", 0.);
270  Real online_Bi = infile("online_Bi", 0.);
271  Real online_kappa = infile("online_kappa", 0.);
272  RBParameters online_mu;
273  online_mu.set_value("curvature", online_curvature);
274  online_mu.set_value("Bi", online_Bi);
275  online_mu.set_value("kappa", online_kappa);
276  rb_eval.set_parameters(online_mu);
277  rb_eval.print_parameters();
278  rb_eval.rb_solve(rb_eval.get_n_basis_functions());
279 
280  // plot the solution, if requested
281  if (store_basis_functions)
282  {
283  // read in the data from files
284  eim_rb_eval.read_in_basis_functions(
285  eim_construction.get_explicit_system(), "eim_data");
286  rb_eval.read_in_basis_functions(rb_construction, "rb_data");
287 
288  eim_construction.load_rb_solution();
289  rb_construction.load_rb_solution();
290 
291  transform_mesh_and_plot(equation_systems, online_curvature, "RB_sol.e");
292  }
293  }
294 
295  return 0;
296 }
297 
299  Real curvature,
300  const std::string & filename)
301 {
302  // Loop over the mesh nodes and move them!
303  MeshBase & mesh = es.get_mesh();
304 
305  MeshBase::node_iterator node_it = mesh.nodes_begin();
306  const MeshBase::node_iterator node_end = mesh.nodes_end();
307 
308  for ( ; node_it != node_end; node_it++)
309  {
310  Node * node = *node_it;
311 
312  Real x = (*node)(0);
313  Real z = (*node)(2);
314 
315  (*node)(0) = -1./curvature + (1./curvature + x)*cos(curvature*z);
316  (*node)(2) = (1./curvature + x)*sin(curvature*z);
317  }
318 
319 #ifdef LIBMESH_HAVE_EXODUS_API
320  ExodusII_IO(mesh).write_equation_systems(filename, es);
321 #endif
322 }
This is the EquationSystems class.
The ReplicatedMesh class is derived from the MeshBase class, and is used to store identical copies of...
A Node is like a Point, but with more information.
Definition: node.h:52
This class serializes an RBEvaluation object using the Cap&#39;n Proto library.
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.
This class serializes an RBEIMEvaluation object using the Cap&#39;n Proto library.
This is the MeshBase class.
Definition: mesh_base.h:68
virtual node_iterator nodes_begin()=0
Iterate over all the nodes in the Mesh.
SolverPackage default_solver_package()
Definition: libmesh.C:995
PetscErrorCode Vec x
void transform_mesh_and_plot(EquationSystems &es, Real curvature, const std::string &filename)
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
The definition of the node_iterator struct.
Definition: mesh_base.h:1528
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
virtual node_iterator nodes_end()=0
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
const MeshBase & get_mesh() 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.
void build_cube(UnstructuredMesh &mesh, const unsigned int nx=0, const unsigned int ny=0, const unsigned int nz=0, const Real xmin=0., const Real xmax=1., const Real ymin=0., const Real ymax=1., const Real zmin=0., const Real zmax=1., const ElemType type=INVALID_ELEM, const bool gauss_lobatto_grid=false)
Builds a (elements) cube.