libMesh
Functions
fem_system_ex1.C File Reference

Go to the source code of this file.

Functions

int main (int argc, char **argv)
 

Function Documentation

int main ( int  argc,
char **  argv 
)

Definition at line 53 of file fem_system_ex1.C.

References libMesh::DiffSolver::absolute_residual_tolerance, libMesh::EquationSystems::add_system(), libMesh::MeshTools::Generation::build_cube(), libMesh::MeshTools::Generation::build_square(), libMesh::LibMeshInit::comm(), libMesh::default_solver_package(), libMesh::DifferentiableSystem::deltat, dim, libMesh::EIGEN_SOLVERS, libMesh::ErrorEstimator::error_norm, libMesh::HEX27, libMesh::TriangleWrapper::init(), libMesh::EquationSystems::init(), libMesh::DiffSolver::initial_linear_tolerance, libMesh::INVALID_SOLVER_PACKAGE, libMesh::L2, libMesh::StatisticsVector< T >::l2_norm(), libMesh::libmesh_assert(), libMesh::DiffSolver::max_linear_iterations, libMesh::DiffSolver::max_nonlinear_iterations, libMesh::StatisticsVector< T >::maximum(), libMesh::ErrorVector::mean(), mesh, libMesh::EquationSystems::n_active_dofs(), libMesh::MeshBase::n_active_elem(), libMesh::out, NavierSystem::postprocess(), libMesh::EquationSystems::print_info(), libMesh::MeshBase::print_info(), libMesh::QUAD9, libMesh::DiffSolver::quiet, libMesh::Real, libMesh::EquationSystems::reinit(), libMesh::DiffSolver::relative_residual_tolerance, libMesh::DiffSolver::relative_step_tolerance, libMesh::FEMSystem::solve(), libMesh::solver, libMesh::System::time, libMesh::DifferentiableSystem::time_solver, libMesh::TRILINOS_SOLVERS, libMesh::DiffSolver::verbose, and libMesh::ExodusII_IO::write_timestep().

54 {
55  // Initialize libMesh.
56  LibMeshInit init (argc, argv);
57 
58  // This example requires a linear solver package.
59  libmesh_example_requires(libMesh::default_solver_package() != INVALID_SOLVER_PACKAGE,
60  "--enable-petsc, --enable-trilinos, or --enable-eigen");
61 
62  // This example fails without at least double precision FP
63 #ifdef LIBMESH_DEFAULT_SINGLE_PRECISION
64  libmesh_example_requires(false, "--disable-singleprecision");
65 #endif
66 
67 #ifndef LIBMESH_ENABLE_AMR
68  libmesh_example_requires(false, "--enable-amr");
69 #else
70 
71  // Trilinos and Eigen solvers NaN by default here.
72  // We'll skip this example for now.
73  libmesh_example_requires(libMesh::default_solver_package() != EIGEN_SOLVERS, "--enable-petsc");
74  libmesh_example_requires(libMesh::default_solver_package() != TRILINOS_SOLVERS, "--enable-petsc");
75 
76  // Parse the input file
77  GetPot infile("fem_system_ex1.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 bool transient = infile("transient", true);
83  const Real deltat = infile("deltat", 0.005);
84  unsigned int n_timesteps = infile("n_timesteps", 20);
85  const unsigned int coarsegridsize = infile("coarsegridsize", 1);
86  const unsigned int coarserefinements = infile("coarserefinements", 0);
87  const unsigned int max_adaptivesteps = infile("max_adaptivesteps", 10);
88  const unsigned int dim = infile("dimension", 2);
89 
90 #ifdef LIBMESH_HAVE_EXODUS_API
91  const unsigned int write_interval = infile("write_interval", 5);
92 #endif
93 
94  // Skip higher-dimensional examples on a lower-dimensional libMesh build
95  libmesh_example_requires(dim <= LIBMESH_DIM, "2D/3D support");
96 
97  // We have only defined 2 and 3 dimensional problems
98  libmesh_assert (dim == 2 || dim == 3);
99 
100  // Create a mesh, with dimension to be overridden later, distributed
101  // across the default MPI communicator.
102  Mesh mesh(init.comm());
103 
104  // And an object to refine it
105  MeshRefinement mesh_refinement(mesh);
106  mesh_refinement.coarsen_by_parents() = true;
107  mesh_refinement.absolute_global_tolerance() = global_tolerance;
108  mesh_refinement.nelem_target() = nelem_target;
109  mesh_refinement.refine_fraction() = 0.3;
110  mesh_refinement.coarsen_fraction() = 0.3;
111  mesh_refinement.coarsen_threshold() = 0.1;
112 
113  // Use the MeshTools::Generation mesh generator to create a uniform
114  // grid on the square [-1,1]^D. We instruct the mesh generator
115  // to build a mesh of 8x8 Quad9 elements in 2D, or Hex27
116  // elements in 3D. Building these higher-order elements allows
117  // us to use higher-order approximation, as in example 3.
118  if (dim == 2)
120  coarsegridsize,
121  coarsegridsize,
122  0., 1.,
123  0., 1.,
124  QUAD9);
125  else if (dim == 3)
127  coarsegridsize,
128  coarsegridsize,
129  coarsegridsize,
130  0., 1.,
131  0., 1.,
132  0., 1.,
133  HEX27);
134 
135  mesh_refinement.uniformly_refine(coarserefinements);
136 
137  // Print information about the mesh to the screen.
138  mesh.print_info();
139 
140  // Create an equation systems object.
141  EquationSystems equation_systems (mesh);
142 
143  // Declare the system "Navier-Stokes" and its variables.
144  NavierSystem & system =
145  equation_systems.add_system<NavierSystem> ("Navier-Stokes");
146 
147  // Solve this as a time-dependent or steady system
148  if (transient)
149  system.time_solver =
150  UniquePtr<TimeSolver>(new EulerSolver(system));
151  else
152  {
153  system.time_solver =
154  UniquePtr<TimeSolver>(new SteadySolver(system));
155  libmesh_assert_equal_to (n_timesteps, 1);
156  }
157 
158  // Initialize the system
159  equation_systems.init ();
160 
161  // Set the time stepping options
162  system.deltat = deltat;
163 
164  // And the nonlinear solver options
165  DiffSolver & solver = *(system.time_solver->diff_solver().get());
166  solver.quiet = infile("solver_quiet", true);
167  solver.verbose = !solver.quiet;
168  solver.max_nonlinear_iterations =
169  infile("max_nonlinear_iterations", 15);
170  solver.relative_step_tolerance =
171  infile("relative_step_tolerance", 1.e-3);
173  infile("relative_residual_tolerance", 0.0);
175  infile("absolute_residual_tolerance", 0.0);
176 
177  // And the linear solver options
178  solver.max_linear_iterations =
179  infile("max_linear_iterations", 50000);
180  solver.initial_linear_tolerance =
181  infile("initial_linear_tolerance", 1.e-3);
182 
183  // Print information about the system to the screen.
184  equation_systems.print_info();
185 
186  // Now we begin the timestep loop to compute the time-accurate
187  // solution of the equations.
188  for (unsigned int t_step=0; t_step != n_timesteps; ++t_step)
189  {
190  // A pretty update message
191  libMesh::out << "\n\nSolving time step "
192  << t_step
193  << ", time = "
194  << system.time
195  << std::endl;
196 
197  // Adaptively solve the timestep
198  unsigned int a_step = 0;
199  for (; a_step != max_adaptivesteps; ++a_step)
200  {
201  system.solve();
202 
203  system.postprocess();
204 
205  ErrorVector error;
206 
207  UniquePtr<ErrorEstimator> error_estimator;
208 
209  // To solve to a tolerance in this problem we
210  // need a better estimator than Kelly
211  if (global_tolerance != 0.)
212  {
213  // We can't adapt to both a tolerance and a mesh
214  // size at once
215  libmesh_assert_equal_to (nelem_target, 0);
216 
218 
219  // The lid-driven cavity problem isn't in H1, so
220  // lets estimate L2 error
221  u->error_norm = L2;
222 
223  error_estimator.reset(u);
224  }
225  else
226  {
227  // If we aren't adapting to a tolerance we need a
228  // target mesh size
229  libmesh_assert_greater (nelem_target, 0);
230 
231  // Kelly is a lousy estimator to use for a problem
232  // not in H1 - if we were doing more than a few
233  // timesteps we'd need to turn off or limit the
234  // maximum level of our adaptivity eventually
235  error_estimator.reset(new KellyErrorEstimator);
236  }
237 
238  // Calculate error based on u and v (and w?) but not p
239  std::vector<Real> weights(2,1.0); // u, v
240  if (dim == 3)
241  weights.push_back(1.0); // w
242  weights.push_back(0.0); // p
243  // Keep the same default norm type.
244  std::vector<FEMNormType>
245  norms(1, error_estimator->error_norm.type(0));
246  error_estimator->error_norm = SystemNorm(norms, weights);
247 
248  error_estimator->estimate_error(system, error);
249 
250  // Print out status at each adaptive step.
251  Real global_error = error.l2_norm();
252  libMesh::out << "Adaptive step "
253  << a_step
254  << ": "
255  << std::endl;
256 
257  if (global_tolerance != 0.)
258  libMesh::out << "Global_error = "
259  << global_error
260  << std::endl;
261 
262  if (global_tolerance != 0.)
263  libMesh::out << "Worst element error = "
264  << error.maximum()
265  << ", mean = "
266  << error.mean()
267  << std::endl;
268 
269  if (global_tolerance != 0.)
270  {
271  // If we've reached our desired tolerance, we
272  // don't need any more adaptive steps
273  if (global_error < global_tolerance)
274  break;
275  mesh_refinement.flag_elements_by_error_tolerance(error);
276  }
277  else
278  {
279  // If flag_elements_by_nelem_target returns true, this
280  // should be our last adaptive step.
281  if (mesh_refinement.flag_elements_by_nelem_target(error))
282  {
283  mesh_refinement.refine_and_coarsen_elements();
284  equation_systems.reinit();
285  a_step = max_adaptivesteps;
286  break;
287  }
288  }
289 
290  // Carry out the adaptive mesh refinement/coarsening
291  mesh_refinement.refine_and_coarsen_elements();
292  equation_systems.reinit();
293 
294  libMesh::out << "Refined mesh to "
295  << mesh.n_active_elem()
296  << " active elements and "
297  << equation_systems.n_active_dofs()
298  << " active dofs."
299  << std::endl;
300  }
301  // Do one last solve if necessary
302  if (a_step == max_adaptivesteps)
303  {
304  system.solve();
305 
306  system.postprocess();
307  }
308 
309  // Advance to the next timestep in a transient problem
310  system.time_solver->advance_timestep();
311 
312 #ifdef LIBMESH_HAVE_EXODUS_API
313  // Write out this timestep if we're requested to
314  if ((t_step+1)%write_interval == 0)
315  {
316  std::ostringstream file_name;
317 
318  // We write the file in the ExodusII format.
319  file_name << "out_"
320  << std::setw(3)
321  << std::setfill('0')
322  << std::right
323  << t_step+1
324  << ".e";
325 
326  ExodusII_IO(mesh).write_timestep(file_name.str(),
327  equation_systems,
328  1, // This number indicates how many time steps
329  // are being written to the file
330  system.time);
331  }
332 #endif // #ifdef LIBMESH_HAVE_EXODUS_API
333  }
334 #endif // #ifndef LIBMESH_ENABLE_AMR
335 
336  // All done.
337  return 0;
338 }
Real time
For time-dependent problems, this is the time t at the beginning of the current timestep.
Definition: system.h:1545
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
void write_timestep(const std::string &fname, const EquationSystems &es, const int timestep, const Real time)
Writes out the solution at a specific timestep.
Definition: exodusII_io.C:773
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
This class defines a norm/seminorm to be applied to a NumericVector which contains coefficients in a ...
Definition: system_norm.h:43
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
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
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.
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
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.
This class implements a TimeSolver which does a single solve of the steady state problem.
Definition: steady_solver.h:47
This class defines a theta-method Euler (defaulting to Backward Euler with theta = 1...
Definition: euler_solver.h:43
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
virtual Real mean() const libmesh_override
Definition: error_vector.C:67
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()
Runs a postprocessing loop over all elements, and if postprocess_sides is true over all sides...
Definition: naviersystem.C:618