68 "--enable-petsc, --enable-trilinos, or --enable-eigen");
71 #ifdef LIBMESH_DEFAULT_SINGLE_PRECISION 72 libmesh_example_requires(
false,
"--disable-singleprecision");
75 #ifndef LIBMESH_ENABLE_AMR 76 libmesh_example_requires(
false,
"--enable-amr");
80 #ifndef LIBMESH_ENABLE_DIRICHLET 81 libmesh_example_requires(
false,
"--enable-dirichlet");
89 GetPot infile(
"fem_system_ex1.in");
92 infile.parse_command_line(argc, argv);
95 const Real global_tolerance = infile(
"global_tolerance", 0.);
96 const unsigned int nelem_target = infile(
"n_elements", 400);
97 const bool transient = infile(
"transient",
true);
98 const Real deltat = infile(
"deltat", 0.005);
99 unsigned int n_timesteps = infile(
"n_timesteps", 20);
100 const unsigned int coarsegridsize = infile(
"coarsegridsize", 1);
101 const unsigned int coarserefinements = infile(
"coarserefinements", 0);
102 const unsigned int max_adaptivesteps = infile(
"max_adaptivesteps", 10);
103 const unsigned int dim = infile(
"dimension", 2);
104 const std::string slvr_type = infile(
"solver_type",
"newton");
105 const std::string mesh_type = infile(
"mesh_type" ,
"replicated");
106 const bool constrain_in_solver = infile(
"constrain_in_solver",
true);
109 const bool print_solutions = infile(
"print_solutions",
false);
110 const bool print_residuals = infile(
"print_residuals",
false);
111 const bool print_jacobians = infile(
"print_jacobians",
false);
113 #ifdef LIBMESH_HAVE_EXODUS_API 114 const unsigned int write_interval = infile(
"write_interval", 5);
118 libmesh_example_requires(
dim <= LIBMESH_DIM,
"2D/3D support");
125 std::shared_ptr<UnstructuredMesh>
mesh;
127 if (mesh_type ==
"distributed")
128 mesh = std::make_shared<DistributedMesh>(
init.comm());
129 else if (mesh_type ==
"replicated")
130 mesh = std::make_shared<ReplicatedMesh>(
init.comm());
132 libmesh_error_msg(
"Error: specified mesh_type not understood");
136 mesh_refinement.coarsen_by_parents() =
true;
137 mesh_refinement.absolute_global_tolerance() = global_tolerance;
138 mesh_refinement.nelem_target() = nelem_target;
139 mesh_refinement.refine_fraction() = 0.3;
140 mesh_refinement.coarsen_fraction() = 0.3;
141 mesh_refinement.coarsen_threshold() = 0.1;
165 mesh_refinement.uniformly_refine(coarserefinements);
175 equation_systems.add_system<
NavierSystem> (
"Navier-Stokes");
187 system.
time_solver = std::make_unique<EulerSolver>(system);
190 system.
time_solver = std::make_unique<SteadySolver>(system);
191 libmesh_assert_equal_to (n_timesteps, 1);
195 equation_systems.init ();
201 if (slvr_type ==
"newton")
202 system.
time_solver->diff_solver() = std::make_unique<NewtonSolver>(system);
203 else if (slvr_type ==
"petscdiff")
204 #if defined(LIBMESH_HAVE_PETSC) && defined(LIBMESH_HAVE_METAPHYSICL) 205 system.
time_solver->diff_solver() = std::make_unique<PetscDiffSolver>(system);
207 libmesh_example_requires(
false,
"--enable-petsc --enable-metaphysicl-required");
210 libmesh_error_msg(
"Error: specified solver_type not understood");
215 solver.
quiet = infile(
"solver_quiet",
true);
218 infile(
"max_nonlinear_iterations", 15);
220 infile(
"relative_step_tolerance", 1.e-3);
222 infile(
"relative_residual_tolerance", 0.0);
224 infile(
"absolute_residual_tolerance", 0.0);
230 infile(
"max_linear_iterations", 50000);
232 infile(
"initial_linear_tolerance", 1.e-3);
235 equation_systems.print_info();
239 for (
unsigned int t_step=0; t_step != n_timesteps; ++t_step)
249 unsigned int a_step = 0;
250 for (; a_step != max_adaptivesteps; ++a_step)
258 std::unique_ptr<ErrorEstimator> error_estimator;
262 if (global_tolerance != 0.)
266 libmesh_assert_equal_to (nelem_target, 0);
268 error_estimator = std::make_unique<UniformRefinementEstimator>();
272 error_estimator->error_norm =
L2;
278 libmesh_assert_greater (nelem_target, 0);
284 error_estimator = std::make_unique<KellyErrorEstimator>();
288 std::vector<Real> weights(2,1.0);
290 weights.push_back(1.0);
291 weights.push_back(0.0);
293 std::vector<FEMNormType>
294 norms(1, error_estimator->error_norm.type(0));
295 error_estimator->error_norm =
SystemNorm(norms, weights);
297 error_estimator->estimate_error(system, error);
306 if (global_tolerance != 0.)
311 if (global_tolerance != 0.)
318 if (global_tolerance != 0.)
322 if (global_error < global_tolerance)
324 mesh_refinement.flag_elements_by_error_tolerance(error);
330 if (mesh_refinement.flag_elements_by_nelem_target(error))
332 mesh_refinement.refine_and_coarsen_elements();
333 equation_systems.reinit();
334 a_step = max_adaptivesteps;
340 mesh_refinement.refine_and_coarsen_elements();
341 equation_systems.reinit();
345 <<
" active elements and " 346 << equation_systems.n_active_dofs()
351 if (a_step == max_adaptivesteps)
361 #ifdef LIBMESH_HAVE_EXODUS_API 363 if ((t_step+1)%write_interval == 0)
365 std::ostringstream file_name;
381 #endif // #ifdef LIBMESH_HAVE_EXODUS_API 383 #endif // #ifndef LIBMESH_ENABLE_AMR virtual T maximum() const
Real time
For time-dependent problems, this is the time t at the beginning of the current timestep.
This is the EquationSystems class.
virtual dof_id_type n_active_elem() const =0
bool quiet
The DiffSolver should not print anything to libMesh::out unless quiet is set to false; default is tru...
virtual Real l2_norm() const
virtual void set_constrain_in_solver(bool enable)
set_constrain_in_solver to false to apply constraints only via residual terms in the systems to be so...
bool print_jacobian_norms
Set print_jacobian_norms to true to print |J| whenever it is assembled.
The ErrorVector is a specialization of the StatisticsVector for error data computed on a finite eleme...
The ExodusII_IO class implements reading meshes in the ExodusII file format from Sandia National Labs...
unsigned int max_nonlinear_iterations
The DiffSolver should exit in failure if max_nonlinear_iterations is exceeded and continue_after_max_...
Real absolute_residual_tolerance
The DiffSolver should exit after the residual is reduced to either less than absolute_residual_tolera...
std::unique_ptr< TimeSolver > time_solver
A pointer to the solver object we're going to use.
This class defines a norm/seminorm to be applied to a NumericVector which contains coefficients in a ...
The LibMeshInit class, when constructed, initializes the dependent libraries (e.g.
bool print_jacobians
Set print_jacobians to true to print J whenever it is assembled.
virtual void init()
The initialization function.
SolverPackage default_solver_package()
This is a generic class that defines a solver to handle ImplicitSystem classes, including NonlinearIm...
bool print_solution_norms
Set print_residual_norms to true to print |U| whenever it is used in an assembly() call...
Implements (adaptive) mesh refinement algorithms for a MeshBase.
bool print_residual_norms
Set print_residual_norms to true to print |F| whenever it is assembled.
Real deltat
For time-dependent problems, this is the amount delta t to advance the solution in time...
unsigned int max_linear_iterations
Each linear solver step should exit after max_linear_iterations is exceeded.
void print_info(std::ostream &os=libMesh::out, const unsigned int verbosity=0, const bool global=true) const
Prints relevant information about the mesh.
bool print_residuals
Set print_residuals to true to print F whenever it is assembled.
void init(triangulateio &t)
Initializes the fields of t to nullptr/0 as necessary.
bool print_solutions
Set print_solutions to true to print U whenever it is used in an assembly() call. ...
bool verbose
The DiffSolver may print a lot more to libMesh::out if verbose is set to true; default is false...
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
void write_timestep(const std::string &fname, const EquationSystems &es, const int timestep, const Real time, const std::set< std::string > *system_names=nullptr)
Writes out the solution at a specific timestep.
virtual void solve() override
Invokes the solver associated with the system.
double initial_linear_tolerance
Any required linear solves will at first be done with this tolerance; the DiffSolver may tighten the ...
Real relative_residual_tolerance
Real relative_step_tolerance
virtual void postprocess()
Runs a postprocessing loop over all elements, and if postprocess_sides is true over all sides...
virtual Real mean() const override