libMesh
vector_fe_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 
19 
20 // <h1>FEMSystem Example 1 - Unsteady Navier-Stokes Equations with
21 // FEMSystem</h1>
22 // \author Paul Bauman
23 // \date 2012
24 //
25 // This example shows how the transient nonlinear problem from
26 // example 13 can be solved using the
27 // DifferentiableSystem class framework
28 
29 // Basic include files
30 #include "libmesh/equation_systems.h"
31 #include "libmesh/getpot.h"
32 #include "libmesh/exodusII_io.h"
33 #include "libmesh/mesh.h"
34 #include "libmesh/mesh_generation.h"
35 #include "libmesh/exact_solution.h"
36 #include "libmesh/ucd_io.h"
37 
38 // The systems and solvers we may use
39 #include "laplace_system.h"
40 #include "libmesh/diff_solver.h"
41 #include "libmesh/steady_solver.h"
42 
43 // Bring in everything from the libMesh namespace
44 using namespace libMesh;
45 
46 // The main program.
47 int main (int argc, char** argv)
48 {
49  // Initialize libMesh.
50  LibMeshInit init (argc, argv);
51 
52  // This example requires a linear solver package.
53  libmesh_example_requires(libMesh::default_solver_package() != INVALID_SOLVER_PACKAGE,
54  "--enable-petsc, --enable-trilinos, or --enable-eigen");
55 
56  // Parse the input file
57  GetPot infile("vector_fe_ex2.in");
58 
59  // Read in parameters from the input file
60  const unsigned int grid_size = infile("grid_size", 2);
61 
62  // Skip higher-dimensional examples on a lower-dimensional libMesh build
63  libmesh_example_requires(3 <= LIBMESH_DIM, "2D/3D support");
64 
65  // Create a mesh, with dimension to be overridden later, on the
66  // default MPI communicator.
67  Mesh mesh(init.comm());
68 
69  // Use the MeshTools::Generation mesh generator to create a uniform
70  // grid on the square [-1,1]^D. We instruct the mesh generator
71  // to build a mesh of 8x8 Quad9 elements in 2D, or Hex27
72  // elements in 3D. Building these higher-order elements allows
73  // us to use higher-order approximation, as in example 3.
75  grid_size,
76  grid_size,
77  grid_size,
78  -1., 1.,
79  -1., 1.,
80  -1., 1.,
81  HEX8);
82 
83  // Print information about the mesh to the screen.
84  mesh.print_info();
85 
86  // Create an equation systems object.
87  EquationSystems equation_systems (mesh);
88 
89  // Declare the system "Navier-Stokes" and its variables.
90  LaplaceSystem & system =
91  equation_systems.add_system<LaplaceSystem> ("Laplace");
92 
93  // This example only implements the steady-state problem
94  system.time_solver =
96 
97  // Initialize the system
98  equation_systems.init();
99 
100  // And the nonlinear solver options
101  DiffSolver & solver = *(system.time_solver->diff_solver().get());
102  solver.quiet = infile("solver_quiet", true);
103  solver.verbose = !solver.quiet;
104  solver.max_nonlinear_iterations = infile("max_nonlinear_iterations", 15);
105  solver.relative_step_tolerance = infile("relative_step_tolerance", 1.e-3);
106  solver.relative_residual_tolerance = infile("relative_residual_tolerance", 0.0);
107  solver.absolute_residual_tolerance = infile("absolute_residual_tolerance", 0.0);
108 
109  // And the linear solver options
110  solver.max_linear_iterations = infile("max_linear_iterations", 50000);
111  solver.initial_linear_tolerance = infile("initial_linear_tolerance", 1.e-3);
112 
113  // Print information about the system to the screen.
114  equation_systems.print_info();
115 
116  system.solve();
117 
118  ExactSolution exact_sol(equation_systems);
119 
120  SolutionFunction soln_func(system.variable_number("u"));
121  SolutionGradient soln_grad(system.variable_number("u"));
122 
123  // Build FunctionBase* containers to attach to the ExactSolution object.
124  std::vector<FunctionBase<Number> *> sols(1, &soln_func);
125  std::vector<FunctionBase<Gradient> *> grads(1, &soln_grad);
126 
127  exact_sol.attach_exact_values(sols);
128  exact_sol.attach_exact_derivs(grads);
129 
130  // Use higher quadrature order for more accurate error results
131  int extra_error_quadrature = infile("extra_error_quadrature", 2);
132  exact_sol.extra_quadrature_order(extra_error_quadrature);
133 
134  // Compute the error.
135  exact_sol.compute_error("Laplace", "u");
136 
137  // Print out the error values
138  libMesh::out << "L2-Error is: "
139  << exact_sol.l2_error("Laplace", "u")
140  << std::endl;
141  libMesh::out << "H1-Error is: "
142  << exact_sol.h1_error("Laplace", "u")
143  << std::endl;
144 
145 #ifdef LIBMESH_HAVE_EXODUS_API
146 
147  // We write the file in the ExodusII format.
148  ExodusII_IO(mesh).write_equation_systems("out.e", equation_systems);
149 
150 #endif // #ifdef LIBMESH_HAVE_EXODUS_API
151 
152  UCDIO(mesh).write_equation_systems("out.inp", equation_systems);
153 
154  // All done.
155  return 0;
156 }
This class handles the computation of the L2 and/or H1 error for the Systems in the EquationSystems o...
UniquePtr< TimeSolver > time_solver
A pointer to the solver object we&#39;re going to use.
Definition: diff_system.h:221
This is the EquationSystems class.
bool quiet
The DiffSolver should not print anything to libMesh::out unless quiet is set to false; default is tru...
Definition: diff_solver.h:162
The ExodusII_IO class implements reading meshes in the ExodusII file format from Sandia National Labs...
Definition: exodusII_io.h:52
unsigned int max_nonlinear_iterations
The DiffSolver should exit in failure if max_nonlinear_iterations is exceeded and continue_after_max_...
Definition: diff_solver.h:156
Real absolute_residual_tolerance
The DiffSolver should exit after the residual is reduced to either less than absolute_residual_tolera...
Definition: diff_solver.h:191
void attach_exact_values(const std::vector< FunctionBase< Number > * > &f)
Clone and attach arbitrary functors which compute the exact values of the EquationSystems&#39; solutions ...
MeshBase & mesh
void extra_quadrature_order(const int extraorder)
Increases or decreases the order of the quadrature rule used for numerical integration.
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.
Real initial_linear_tolerance
Any required linear solves will at first be done with this tolerance; the DiffSolver may tighten the ...
Definition: diff_solver.h:210
void compute_error(const std::string &sys_name, const std::string &unknown_name)
Computes and stores the error in the solution value e = u-u_h, the gradient grad(e) = grad(u) - grad(...
std::unique_ptr< T > UniquePtr
Definition: auto_ptr.h:46
PetscDiffSolver & solver
SolverPackage default_solver_package()
Definition: libmesh.C:995
This is a generic class that defines a solver to handle ImplicitSystem classes, including NonlinearIm...
Definition: diff_solver.h:70
unsigned short int variable_number(const std::string &var) const
Definition: system.C:1263
virtual void solve() libmesh_override
Invokes the solver associated with the system.
Definition: fem_system.C:1056
This class implements reading & writing meshes in the AVS&#39;s UCD format.
Definition: ucd_io.h:44
void print_info(std::ostream &os=libMesh::out) const
Prints information about the equation systems, by default to libMesh::out.
unsigned int max_linear_iterations
Each linear solver step should exit after max_linear_iterations is exceeded.
Definition: diff_solver.h:148
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.
This class implements a TimeSolver which does a single solve of the steady state problem.
Definition: steady_solver.h:47
const Parallel::Communicator & comm() const
Definition: libmesh.h:81
void attach_exact_derivs(const std::vector< FunctionBase< Gradient > * > &g)
Clone and attach arbitrary functors which compute the exact gradients of the EquationSystems&#39; solutio...
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
bool verbose
The DiffSolver may print a lot more to libMesh::out if verbose is set to true; default is false...
Definition: diff_solver.h:168
Real h1_error(const std::string &sys_name, const std::string &unknown_name)
OStreamProxy out
virtual void init()
Initialize all the systems.
Real l2_error(const std::string &sys_name, const std::string &unknown_name)
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 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.
Real relative_residual_tolerance
Definition: diff_solver.h:192
int main(int argc, char **argv)
Definition: vector_fe_ex2.C:47