libMesh
Functions
fem_system_ex4.C File Reference

Go to the source code of this file.

Functions

int main (int argc, char **argv)
 

Function Documentation

◆ main()

int main ( int  argc,
char **  argv 
)

Definition at line 60 of file fem_system_ex4.C.

References libMesh::DiffSolver::absolute_residual_tolerance, libMesh::BoundaryInfo::add_elements(), libMesh::EquationSystems::add_system(), libMesh::ExactSolution::attach_exact_value(), libMesh::MeshTools::Generation::build_cube(), libMesh::MeshTools::Generation::build_square(), libMesh::ExactSolution::compute_error(), libMesh::default_solver_package(), libMesh::DifferentiableSystem::deltat, dim, libMesh::err, libMesh::MeshBase::get_boundary_info(), libMesh::HEX27, libMesh::TriangleWrapper::init(), libMesh::EquationSystems::init(), libMesh::DiffSolver::initial_linear_tolerance, libMesh::INVALID_SOLVER_PACKAGE, libMesh::L2, libMesh::ExactSolution::l2_error(), 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, libMesh::FEMSystem::postprocess(), libMesh::MeshBase::prepare_for_use(), 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::DifferentiableSystem::time_solver, libMesh::DiffSolver::verbose, and libMesh::MeshOutput< MT >::write_equation_systems().

61 {
62  // Initialize libMesh.
63  LibMeshInit init (argc, argv);
64 
65  // This example requires a linear solver package.
66  libmesh_example_requires(libMesh::default_solver_package() != INVALID_SOLVER_PACKAGE,
67  "--enable-petsc, --enable-trilinos, or --enable-eigen");
68 
69 #ifndef LIBMESH_ENABLE_AMR
70  libmesh_example_requires(false, "--enable-amr");
71 #else
72 
73  // We use Dirichlet boundary conditions here
74 #ifndef LIBMESH_ENABLE_DIRICHLET
75  libmesh_example_requires(false, "--enable-dirichlet");
76 #endif
77 
78  // This doesn't converge without at least double precision
79  libmesh_example_requires(sizeof(Real) > 4, "--disable-singleprecision");
80 
81  // Parse the input file
82  GetPot infile("fem_system_ex4.in");
83 
84  // But allow the command line to override it.
85  infile.parse_command_line(argc, argv);
86 
87  // Read in parameters from the input file
88  const Real global_tolerance = infile("global_tolerance", 0.);
89  const unsigned int nelem_target = infile("n_elements", 400);
90  const Real deltat = infile("deltat", 0.005);
91  const unsigned int coarsegridsize = infile("coarsegridsize", 20);
92  const unsigned int coarserefinements = infile("coarserefinements", 0);
93  const unsigned int max_adaptivesteps = infile("max_adaptivesteps", 10);
94  const unsigned int dim = infile("dimension", 2);
95 
96  // Skip higher-dimensional examples on a lower-dimensional libMesh build
97  libmesh_example_requires(dim <= LIBMESH_DIM, "2D/3D support");
98 
99  // We have only defined 2 and 3 dimensional problems
100  libmesh_assert (dim == 2 || dim == 3);
101 
102  // Create a mesh, with dimension to be overridden later, distributed
103  // across the default MPI communicator.
104  Mesh mesh(init.comm());
105 
106  // And an object to refine it
107  MeshRefinement mesh_refinement(mesh);
108  mesh_refinement.coarsen_by_parents() = true;
109  mesh_refinement.absolute_global_tolerance() = global_tolerance;
110  mesh_refinement.nelem_target() = nelem_target;
111  mesh_refinement.refine_fraction() = 0.3;
112  mesh_refinement.coarsen_fraction() = 0.3;
113  mesh_refinement.coarsen_threshold() = 0.1;
114 
115  // Use the MeshTools::Generation mesh generator to create a uniform
116  // grid on the square or cube. We crop the domain at y=2/3 to allow
117  // for a homogeneous Neumann BC in our benchmark there.
118  boundary_id_type bcid = 3; // +y in 3D
119  if (dim == 2)
120  {
122  (mesh,
123  coarsegridsize,
124  coarsegridsize*2/3, // int arithmetic best we can do here
125  0., 1.,
126  0., 2./3.,
127  QUAD9);
128  bcid = 2; // +y in 2D
129  }
130  else if (dim == 3)
131  {
133  (mesh,
134  coarsegridsize,
135  coarsegridsize*2/3,
136  coarsegridsize,
137  0., 1.,
138  0., 2./3.,
139  0., 1.,
140  HEX27);
141  }
142 
143  {
144  // Add boundary elements corresponding to the +y boundary of our
145  // volume mesh
146  std::set<boundary_id_type> bcids;
147  bcids.insert(bcid);
150  }
151 
152  // To work around ExodusII file format limitations, we need elements
153  // of different dimensionality to belong to different subdomains.
154  // Our interior elements defaulted to subdomain id 0, so we'll set
155  // boundary elements to subdomain 1.
156  for (auto & elem : mesh.element_ptr_range())
157  if (elem->dim() < dim)
158  elem->subdomain_id() = 1;
159 
160  mesh_refinement.uniformly_refine(coarserefinements);
161 
162  // Print information about the mesh to the screen.
163  mesh.print_info();
164 
165  // Create an equation systems object.
166  EquationSystems equation_systems (mesh);
167 
168  // Declare the system "Heat" and its variables.
169  HeatSystem & system =
170  equation_systems.add_system<HeatSystem> ("Heat");
171 
172  // Solve this as a steady system
173  system.time_solver = std::make_unique<SteadySolver>(system);
174 
175  // Initialize the system
176  equation_systems.init ();
177 
178  // Set the time stepping options
179  system.deltat = deltat;
180 
181  // And the nonlinear solver options
182  DiffSolver & solver = *(system.time_solver->diff_solver().get());
183  solver.quiet = infile("solver_quiet", true);
184  solver.verbose = !solver.quiet;
185  solver.max_nonlinear_iterations = infile("max_nonlinear_iterations", 15);
186  solver.relative_step_tolerance = infile("relative_step_tolerance", 1.e-3);
187  solver.relative_residual_tolerance = infile("relative_residual_tolerance", 0.0);
188  solver.absolute_residual_tolerance = infile("absolute_residual_tolerance", 0.0);
189 
190  // And the linear solver options
191  solver.max_linear_iterations = infile("max_linear_iterations", 50000);
192  solver.initial_linear_tolerance = infile("initial_linear_tolerance", 1.e-3);
193 
194  // Print information about the system to the screen.
195  equation_systems.print_info();
196 
197  // Adaptively solve the steady solution
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  std::unique_ptr<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 
217  auto u = std::make_unique<UniformRefinementEstimator>();
218 
219  // The lid-driven cavity problem isn't in H1, so
220  // lets estimate L2 error
221  u->error_norm = L2;
222  error_estimator = std::move(u);
223  }
224  else
225  {
226  // If we aren't adapting to a tolerance we need a
227  // target mesh size
228  libmesh_assert_greater (nelem_target, 0);
229 
230  // Kelly is a lousy estimator to use for a problem
231  // not in H1 - if we were doing more than a few
232  // timesteps we'd need to turn off or limit the
233  // maximum level of our adaptivity eventually
234  error_estimator = std::make_unique<KellyErrorEstimator>();
235  }
236 
237  error_estimator->estimate_error(system, error);
238 
239  // Print out status at each adaptive step.
240  Real global_error = error.l2_norm();
241  libMesh::out << "Adaptive step "
242  << a_step
243  << ": "
244  << std::endl;
245 
246  if (global_tolerance != 0.)
247  libMesh::out << "Global_error = "
248  << global_error
249  << std::endl;
250 
251  if (global_tolerance != 0.)
252  libMesh::out << "Worst element error = "
253  << error.maximum()
254  << ", mean = "
255  << error.mean()
256  << std::endl;
257 
258  if (global_tolerance != 0.)
259  {
260  // If we've reached our desired tolerance, we
261  // don't need any more adaptive steps
262  if (global_error < global_tolerance)
263  break;
264  mesh_refinement.flag_elements_by_error_tolerance(error);
265  }
266  else
267  {
268  // If flag_elements_by_nelem_target returns true, this
269  // should be our last adaptive step.
270  if (mesh_refinement.flag_elements_by_nelem_target(error))
271  {
272  mesh_refinement.refine_and_coarsen_elements();
273  equation_systems.reinit();
274  a_step = max_adaptivesteps;
275  break;
276  }
277  }
278 
279  // Carry out the adaptive mesh refinement/coarsening
280  mesh_refinement.refine_and_coarsen_elements();
281  equation_systems.reinit();
282 
283  libMesh::out << "Refined mesh to "
284  << mesh.n_active_elem()
285  << " active elements and "
286  << equation_systems.n_active_dofs()
287  << " active dofs."
288  << std::endl;
289  }
290  // Do one last solve if necessary
291  if (a_step == max_adaptivesteps)
292  {
293  system.solve();
294 
295  system.postprocess();
296  }
297 
298 
299 #ifdef LIBMESH_HAVE_EXODUS_API
301  ("out.e", equation_systems);
302 #endif // #ifdef LIBMESH_HAVE_EXODUS_API
303 
304 #ifdef LIBMESH_HAVE_GMV
306  ("out.gmv", equation_systems);
307 #endif // #ifdef LIBMESH_HAVE_GMV
308 
309 #ifdef LIBMESH_HAVE_FPARSER
310  // Check that we got close to the analytic solution
311  ExactSolution exact_sol(equation_systems);
312  const std::string exact_str = (dim == 2) ?
313  "sin(pi*x)*sin(pi*y)" : "sin(pi*x)*sin(pi*y)*sin(pi*z)";
314  ParsedFunction<Number> exact_func(exact_str);
315  exact_sol.attach_exact_value(0, &exact_func);
316  exact_sol.compute_error("Heat", "T");
317 
318  Real err = exact_sol.l2_error("Heat", "T");
319 
320  // Print out the error value
321  libMesh::out << "L2-Error is: " << err << std::endl;
322 
323  libmesh_assert_less(err, 2e-3);
324 
325 #endif // #ifdef LIBMESH_HAVE_FPARSER
326 
327 #endif // #ifndef LIBMESH_ENABLE_AMR
328 
329  // All done.
330  return 0;
331 }
virtual T maximum() const
Definition: statistics.C:62
OStreamProxy err
This class handles the computation of the L2 and/or H1 error for the Systems in the EquationSystems o...
This is the EquationSystems class.
virtual dof_id_type n_active_elem() const =0
void add_elements(const std::set< boundary_id_type > &requested_boundary_ids, UnstructuredMesh &boundary_mesh, bool store_parent_side_ids=false)
Generates elements along the boundary of our _mesh, which use pre-existing nodes on the boundary_mesh...
A Function generated (via FParser) by parsing a mathematical expression.
virtual void write_equation_systems(const std::string &, const EquationSystems &, const std::set< std::string > *system_names=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 quiet
The DiffSolver should not print anything to libMesh::out unless quiet is set to false; default is tru...
Definition: diff_solver.h:160
virtual Real l2_norm() const
Definition: statistics.C:37
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:154
void prepare_for_use(const bool skip_renumber_nodes_and_elements, const bool skip_find_neighbors)
Prepare a newly ecreated (or read) mesh for use.
Definition: mesh_base.C:710
Real absolute_residual_tolerance
The DiffSolver should exit after the residual is reduced to either less than absolute_residual_tolera...
Definition: diff_solver.h:189
std::unique_ptr< TimeSolver > time_solver
A pointer to the solver object we&#39;re going to use.
Definition: diff_system.h:251
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:90
This class implements writing meshes in the GMV format.
Definition: gmv_io.h:46
const BoundaryInfo & get_boundary_info() const
The information about boundary ids on the mesh.
Definition: mesh_base.h:159
SolverPackage default_solver_package()
Definition: libmesh.C:1050
This is a generic class that defines a solver to handle ImplicitSystem classes, including NonlinearIm...
Definition: diff_solver.h:68
Implements (adaptive) mesh refinement algorithms for a MeshBase.
int8_t boundary_id_type
Definition: id_types.h:51
Real deltat
For time-dependent problems, this is the amount delta t to advance the solution in time...
Definition: diff_system.h:278
unsigned int max_linear_iterations
Each linear solver step should exit after max_linear_iterations is exceeded.
Definition: diff_solver.h:146
void print_info(std::ostream &os=libMesh::out, const unsigned int verbosity=0, const bool global=true) const
Prints relevant information about the mesh.
Definition: mesh_base.C:1489
void init(triangulateio &t)
Initializes the fields of t to nullptr/0 as necessary.
libmesh_assert(ctx)
virtual void postprocess() override
Runs a postprocessing loop over all elements, and if postprocess_sides is true over all sides...
Definition: fem_system.C:1127
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:166
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
OStreamProxy out
virtual void solve() override
Invokes the solver associated with the system.
Definition: fem_system.C:1076
The Mesh class is a thin wrapper, around the ReplicatedMesh class by default.
Definition: mesh.h:50
double initial_linear_tolerance
Any required linear solves will at first be done with this tolerance; the DiffSolver may tighten the ...
Definition: diff_solver.h:208
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:190
virtual Real mean() const override
Definition: error_vector.C:69