libMesh
fem_system_ex4.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 4 - Mixed-dimension heat transfer equation
21 // with FEMSystem</h1>
22 // \author Roy Stogner
23 // \date 2015
24 //
25 // This example shows how elements of multiple dimensions can be
26 // linked and computed upon simultaneously within the
27 // DifferentiableSystem class framework
28 
29 // C++ includes
30 #include <iomanip>
31 
32 // Basic include files
33 #include "libmesh/elem.h"
34 #include "libmesh/equation_systems.h"
35 #include "libmesh/error_vector.h"
36 #include "libmesh/exact_solution.h"
37 #include "libmesh/getpot.h"
38 #include "libmesh/gmv_io.h"
39 #include "libmesh/exodusII_io.h"
40 #include "libmesh/kelly_error_estimator.h"
41 #include "libmesh/mesh.h"
42 #include "libmesh/mesh_generation.h"
43 #include "libmesh/mesh_refinement.h"
44 #include "libmesh/parsed_function.h"
45 #include "libmesh/uniform_refinement_estimator.h"
46 
47 // The systems and solvers we may use
48 #include "heatsystem.h"
49 #include "libmesh/diff_solver.h"
50 #include "libmesh/euler_solver.h"
51 #include "libmesh/steady_solver.h"
52 
53 // Bring in everything from the libMesh namespace
54 using namespace libMesh;
55 
56 // The main program.
57 int main (int argc, char ** argv)
58 {
59  // Initialize libMesh.
60  LibMeshInit init (argc, argv);
61 
62  // This example requires a linear solver package.
63  libmesh_example_requires(libMesh::default_solver_package() != INVALID_SOLVER_PACKAGE,
64  "--enable-petsc, --enable-trilinos, or --enable-eigen");
65 
66 #ifndef LIBMESH_ENABLE_AMR
67  libmesh_example_requires(false, "--enable-amr");
68 #else
69 
70  // This doesn't converge with Eigen BICGSTAB for some reason...
71  libmesh_example_requires(libMesh::default_solver_package() != EIGEN_SOLVERS, "--enable-petsc");
72 
73  // This doesn't converge without at least double precision
74  libmesh_example_requires(sizeof(Real) > 4, "--disable-singleprecision");
75 
76  // Parse the input file
77  GetPot infile("fem_system_ex4.in");
78 
79  // Read in parameters from the input file
80  const Real global_tolerance = infile("global_tolerance", 0.);
81  const unsigned int nelem_target = infile("n_elements", 400);
82  const Real deltat = infile("deltat", 0.005);
83  const unsigned int coarsegridsize = infile("coarsegridsize", 20);
84  const unsigned int coarserefinements = infile("coarserefinements", 0);
85  const unsigned int max_adaptivesteps = infile("max_adaptivesteps", 10);
86  const unsigned int dim = infile("dimension", 2);
87 
88  // Skip higher-dimensional examples on a lower-dimensional libMesh build
89  libmesh_example_requires(dim <= LIBMESH_DIM, "2D/3D support");
90 
91  // We have only defined 2 and 3 dimensional problems
92  libmesh_assert (dim == 2 || dim == 3);
93 
94  // Create a mesh, with dimension to be overridden later, distributed
95  // across the default MPI communicator.
96  Mesh mesh(init.comm());
97 
98  // And an object to refine it
99  MeshRefinement mesh_refinement(mesh);
100  mesh_refinement.coarsen_by_parents() = true;
101  mesh_refinement.absolute_global_tolerance() = global_tolerance;
102  mesh_refinement.nelem_target() = nelem_target;
103  mesh_refinement.refine_fraction() = 0.3;
104  mesh_refinement.coarsen_fraction() = 0.3;
105  mesh_refinement.coarsen_threshold() = 0.1;
106 
107  // Use the MeshTools::Generation mesh generator to create a uniform
108  // grid on the square or cube. We crop the domain at y=2/3 to allow
109  // for a homogeneous Neumann BC in our benchmark there.
110  boundary_id_type bcid = 3; // +y in 3D
111  if (dim == 2)
112  {
114  (mesh,
115  coarsegridsize,
116  coarsegridsize*2/3, // int arithmetic best we can do here
117  0., 1.,
118  0., 2./3.,
119  QUAD9);
120  bcid = 2; // +y in 2D
121  }
122  else if (dim == 3)
123  {
125  (mesh,
126  coarsegridsize,
127  coarsegridsize*2/3,
128  coarsegridsize,
129  0., 1.,
130  0., 2./3.,
131  0., 1.,
132  HEX27);
133  }
134 
135  {
136  // Add boundary elements corresponding to the +y boundary of our
137  // volume mesh
138  std::set<boundary_id_type> bcids;
139  bcids.insert(bcid);
142  }
143 
144  // To work around ExodusII file format limitations, we need elements
145  // of different dimensionality to belong to different subdomains.
146  // Our interior elements defaulted to subdomain id 0, so we'll set
147  // boundary elements to subdomain 1.
148  for (auto & elem : mesh.element_ptr_range())
149  if (elem->dim() < dim)
150  elem->subdomain_id() = 1;
151 
152  mesh_refinement.uniformly_refine(coarserefinements);
153 
154  // Print information about the mesh to the screen.
155  mesh.print_info();
156 
157  // Create an equation systems object.
158  EquationSystems equation_systems (mesh);
159 
160  // Declare the system "Heat" and its variables.
161  HeatSystem & system =
162  equation_systems.add_system<HeatSystem> ("Heat");
163 
164  // Solve this as a steady system
165  system.time_solver =
166  UniquePtr<TimeSolver>(new SteadySolver(system));
167 
168  // Initialize the system
169  equation_systems.init ();
170 
171  // Set the time stepping options
172  system.deltat = deltat;
173 
174  // And the nonlinear solver options
175  DiffSolver & solver = *(system.time_solver->diff_solver().get());
176  solver.quiet = infile("solver_quiet", true);
177  solver.verbose = !solver.quiet;
178  solver.max_nonlinear_iterations = infile("max_nonlinear_iterations", 15);
179  solver.relative_step_tolerance = infile("relative_step_tolerance", 1.e-3);
180  solver.relative_residual_tolerance = infile("relative_residual_tolerance", 0.0);
181  solver.absolute_residual_tolerance = infile("absolute_residual_tolerance", 0.0);
182 
183  // And the linear solver options
184  solver.max_linear_iterations = infile("max_linear_iterations", 50000);
185  solver.initial_linear_tolerance = infile("initial_linear_tolerance", 1.e-3);
186 
187  // Print information about the system to the screen.
188  equation_systems.print_info();
189 
190  // Adaptively solve the steady solution
191  unsigned int a_step = 0;
192  for (; a_step != max_adaptivesteps; ++a_step)
193  {
194  system.solve();
195 
196  system.postprocess();
197 
198  ErrorVector error;
199 
200  UniquePtr<ErrorEstimator> error_estimator;
201 
202  // To solve to a tolerance in this problem we
203  // need a better estimator than Kelly
204  if (global_tolerance != 0.)
205  {
206  // We can't adapt to both a tolerance and a mesh
207  // size at once
208  libmesh_assert_equal_to (nelem_target, 0);
209 
212 
213  // The lid-driven cavity problem isn't in H1, so
214  // lets estimate L2 error
215  u->error_norm = L2;
216 
217  error_estimator.reset(u);
218  }
219  else
220  {
221  // If we aren't adapting to a tolerance we need a
222  // target mesh size
223  libmesh_assert_greater (nelem_target, 0);
224 
225  // Kelly is a lousy estimator to use for a problem
226  // not in H1 - if we were doing more than a few
227  // timesteps we'd need to turn off or limit the
228  // maximum level of our adaptivity eventually
229  error_estimator.reset(new KellyErrorEstimator);
230  }
231 
232  error_estimator->estimate_error(system, error);
233 
234  // Print out status at each adaptive step.
235  Real global_error = error.l2_norm();
236  libMesh::out << "Adaptive step "
237  << a_step
238  << ": "
239  << std::endl;
240 
241  if (global_tolerance != 0.)
242  libMesh::out << "Global_error = "
243  << global_error
244  << std::endl;
245 
246  if (global_tolerance != 0.)
247  libMesh::out << "Worst element error = "
248  << error.maximum()
249  << ", mean = "
250  << error.mean()
251  << std::endl;
252 
253  if (global_tolerance != 0.)
254  {
255  // If we've reached our desired tolerance, we
256  // don't need any more adaptive steps
257  if (global_error < global_tolerance)
258  break;
259  mesh_refinement.flag_elements_by_error_tolerance(error);
260  }
261  else
262  {
263  // If flag_elements_by_nelem_target returns true, this
264  // should be our last adaptive step.
265  if (mesh_refinement.flag_elements_by_nelem_target(error))
266  {
267  mesh_refinement.refine_and_coarsen_elements();
268  equation_systems.reinit();
269  a_step = max_adaptivesteps;
270  break;
271  }
272  }
273 
274  // Carry out the adaptive mesh refinement/coarsening
275  mesh_refinement.refine_and_coarsen_elements();
276  equation_systems.reinit();
277 
278  libMesh::out << "Refined mesh to "
279  << mesh.n_active_elem()
280  << " active elements and "
281  << equation_systems.n_active_dofs()
282  << " active dofs."
283  << std::endl;
284  }
285  // Do one last solve if necessary
286  if (a_step == max_adaptivesteps)
287  {
288  system.solve();
289 
290  system.postprocess();
291  }
292 
293 
294 #ifdef LIBMESH_HAVE_EXODUS_API
296  ("out.e", equation_systems);
297 #endif // #ifdef LIBMESH_HAVE_EXODUS_API
298 
299 #ifdef LIBMESH_HAVE_GMV
301  ("out.gmv", equation_systems);
302 #endif // #ifdef LIBMESH_HAVE_GMV
303 
304 #ifdef LIBMESH_HAVE_FPARSER
305  // Check that we got close to the analytic solution
306  ExactSolution exact_sol(equation_systems);
307  const std::string exact_str = (dim == 2) ?
308  "sin(pi*x)*sin(pi*y)" : "sin(pi*x)*sin(pi*y)*sin(pi*z)";
309  ParsedFunction<Number> exact_func(exact_str);
310  exact_sol.attach_exact_value(0, &exact_func);
311  exact_sol.compute_error("Heat", "T");
312 
313  Number err = exact_sol.l2_error("Heat", "T");
314 
315  // Print out the error value
316  libMesh::out << "L2-Error is: " << err << std::endl;
317 
318  libmesh_assert_less(libmesh_real(err), 2e-3);
319 
320 #endif // #ifdef LIBMESH_HAVE_FPARSER
321 
322 #endif // #ifndef LIBMESH_ENABLE_AMR
323 
324  // All done.
325  return 0;
326 }
T libmesh_real(T a)
OStreamProxy err
This class handles the computation of the L2 and/or H1 error for the Systems in the EquationSystems o...
void add_elements(const std::set< boundary_id_type > &requested_boundary_ids, UnstructuredMesh &boundary_mesh)
Generates elements along the boundary of our _mesh, which use pre-existing nodes on the boundary_mesh...
const BoundaryInfo & get_boundary_info() const
The information about boundary ids on the mesh.
Definition: mesh_base.h:117
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.
virtual dof_id_type n_active_elem() const =0
A Function generated (via FParser) by parsing a mathematical expression.
SystemNorm error_norm
When estimating the error in a single system, the error_norm is used to control the scaling and norm ...
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
unsigned int dim
The ErrorVector is a specialization of the StatisticsVector for error data computed on a finite eleme...
Definition: error_vector.h:50
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
virtual Real l2_norm() const
Definition: statistics.C:36
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
virtual T maximum() const
Definition: statistics.C:61
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.
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
This class implements writing meshes in the GMV format.
Definition: gmv_io.h:46
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(...
libmesh_assert(j)
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
virtual void solve() libmesh_override
Invokes the solver associated with the system.
Definition: fem_system.C:1056
This is the MeshRefinement class.
int8_t boundary_id_type
Definition: id_types.h:51
int main(int argc, char **argv)
virtual SimpleRange< element_iterator > element_ptr_range()=0
void print_info(std::ostream &os=libMesh::out) const
Prints information about the equation systems, by default to libMesh::out.
Real deltat
For time-dependent problems, this is the amount delta t to advance the solution in time...
Definition: diff_system.h:248
unsigned int max_linear_iterations
Each linear solver step should exit after max_linear_iterations is exceeded.
Definition: diff_solver.h:148
virtual void reinit()
Reinitialize all the systems.
This class implements a ``brute force&#39;&#39; error estimator which integrates differences between the curr...
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
void prepare_for_use(const bool skip_renumber_nodes_and_elements=false, const bool skip_find_neighbors=false)
Prepare a newly created (or read) mesh for use.
Definition: mesh_base.C:174
const Parallel::Communicator & comm() const
Definition: libmesh.h:81
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
This class implements the Kelly error indicator which is based on the flux jumps between elements...
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
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
OStreamProxy out
void attach_exact_value(unsigned int sys_num, FunctionBase< Number > *f)
Clone and attach an arbitrary functor which computes the exact value of the system sys_num solution a...
virtual void init()
Initialize all the systems.
virtual Real mean() const libmesh_override
Definition: error_vector.C:67
Real l2_error(const std::string &sys_name, const std::string &unknown_name)
std::size_t n_active_dofs() const
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
virtual void postprocess() libmesh_override
Runs a postprocessing loop over all elements, and if postprocess_sides is true over all sides...
Definition: fem_system.C:1107