36 #include "libmesh/libmesh.h" 37 #include "libmesh/mesh.h" 38 #include "libmesh/mesh_generation.h" 39 #include "libmesh/exodusII_io.h" 40 #include "libmesh/equation_systems.h" 41 #include "libmesh/fe.h" 42 #include "libmesh/quadrature_gauss.h" 43 #include "libmesh/dof_map.h" 44 #include "libmesh/sparse_matrix.h" 45 #include "libmesh/numeric_vector.h" 46 #include "libmesh/dense_matrix.h" 47 #include "libmesh/dense_vector.h" 48 #include "libmesh/linear_implicit_system.h" 49 #include "libmesh/transient_system.h" 50 #include "libmesh/perf_log.h" 51 #include "libmesh/boundary_info.h" 52 #include "libmesh/utility.h" 53 #include "libmesh/dirichlet_boundaries.h" 54 #include "libmesh/zero_function.h" 55 #include "libmesh/const_function.h" 56 #include "libmesh/enum_solver_package.h" 57 #include "libmesh/getpot.h" 63 #include "libmesh/dense_submatrix.h" 64 #include "libmesh/dense_subvector.h" 67 #include "libmesh/elem.h" 75 const std::string & system_name);
81 int main (
int argc,
char ** argv)
88 "--enable-petsc, --enable-trilinos, or --enable-eigen");
91 libmesh_example_requires(2 <= LIBMESH_DIM,
"2D support");
94 #ifndef LIBMESH_ENABLE_DIRICHLET 95 libmesh_example_requires(
false,
"--enable-dirichlet");
132 system.add_variable (
"vel_x",
SECOND);
133 system.add_variable (
"vel_y",
SECOND);
138 system.add_variable (
"p",
FIRST);
152 equation_systems.
init ();
158 PerfLog perf_log(
"Systems Example 3");
167 navier_stokes_system.time = 0.0;
168 const unsigned int n_timesteps = 15;
172 const unsigned int n_nonlinear_steps = 15;
179 equation_systems.
parameters.
set<
unsigned int>(
"linear solver maximum iterations") = max_iter;
192 std::unique_ptr<NumericVector<Number>>
193 last_nonlinear_soln (navier_stokes_system.solution->clone());
195 #ifdef LIBMESH_HAVE_EXODUS_API 200 exo_io.write_equation_systems (
"out.e", equation_systems);
203 for (
unsigned int t_step=1; t_step<=n_timesteps; ++t_step)
207 navier_stokes_system.time += dt;
213 << navier_stokes_system.time
220 *navier_stokes_system.
old_local_solution = *navier_stokes_system.current_local_solution;
224 const Real initial_linear_solver_tol = 1.e-6;
225 equation_systems.
parameters.
set<
Real> (
"linear solver tolerance") = initial_linear_solver_tol;
228 bool converged =
false;
231 for (
unsigned int l=0; l<n_nonlinear_steps; ++l)
234 last_nonlinear_soln->zero();
235 last_nonlinear_soln->add(*navier_stokes_system.solution);
238 perf_log.
push(
"linear solve");
239 equation_systems.
get_system(
"Navier-Stokes").solve();
240 perf_log.
pop(
"linear solve");
244 last_nonlinear_soln->add (-1., *navier_stokes_system.solution);
247 last_nonlinear_soln->close();
250 const Real norm_delta = last_nonlinear_soln->l2_norm();
253 const unsigned int n_linear_iterations = navier_stokes_system.n_linear_iterations();
256 const Real final_linear_residual = navier_stokes_system.final_linear_residual();
271 if (n_linear_iterations == 0 &&
272 (navier_stokes_system.final_linear_residual() >= nonlinear_tolerance || l==0))
274 Real old_linear_solver_tolerance = equation_systems.
parameters.
get<
Real> (
"linear solver tolerance");
275 equation_systems.
parameters.
set<
Real> (
"linear solver tolerance") = 1.e-3 * old_linear_solver_tolerance;
282 << n_linear_iterations
283 <<
", final residual: " 284 << final_linear_residual
285 <<
" Nonlinear convergence: ||u - u_old|| = " 292 if ((norm_delta < nonlinear_tolerance) &&
293 (navier_stokes_system.final_linear_residual() < nonlinear_tolerance))
307 Real new_linear_solver_tolerance = std::min(Utility::pow<2>(final_linear_residual), initial_linear_solver_tol);
308 equation_systems.
parameters.
set<
Real> (
"linear solver tolerance") = new_linear_solver_tolerance;
312 libmesh_error_msg_if(!converged,
"Error: Newton iterations failed to converge!");
314 #ifdef LIBMESH_HAVE_EXODUS_API 316 const unsigned int write_interval = 1;
318 if ((t_step+1)%write_interval == 0)
320 exo_io.write_timestep(
"out.e",
323 navier_stokes_system.time);
325 #endif // #ifdef LIBMESH_HAVE_EXODUS_API 340 const std::string & libmesh_dbg_var(system_name))
344 libmesh_assert_equal_to (system_name,
"Navier-Stokes");
358 const unsigned int u_var = navier_stokes_system.variable_number (
"vel_x");
359 const unsigned int v_var = navier_stokes_system.variable_number (
"vel_y");
360 const unsigned int p_var = navier_stokes_system.variable_number (
"p");
361 const unsigned int alpha_var = navier_stokes_system.variable_number (
"alpha");
365 FEType fe_vel_type = navier_stokes_system.variable_type(u_var);
368 FEType fe_pres_type = navier_stokes_system.variable_type(p_var);
383 fe_vel->attach_quadrature_rule (&qrule);
384 fe_pres->attach_quadrature_rule (&qrule);
390 const std::vector<Real> & JxW = fe_vel->get_JxW();
393 const std::vector<std::vector<Real>> & phi = fe_vel->get_phi();
397 const std::vector<std::vector<RealGradient>> & dphi = fe_vel->get_dphi();
401 const std::vector<std::vector<Real>> & psi = fe_pres->get_phi();
410 const DofMap & dof_map = navier_stokes_system.get_dof_map();
420 Kuu(Ke), Kuv(Ke), Kup(Ke),
421 Kvu(Ke), Kvv(Ke), Kvp(Ke),
422 Kpu(Ke), Kpv(Ke), Kpp(Ke);
433 std::vector<dof_id_type> dof_indices;
434 std::vector<dof_id_type> dof_indices_u;
435 std::vector<dof_id_type> dof_indices_v;
436 std::vector<dof_id_type> dof_indices_p;
437 std::vector<dof_id_type> dof_indices_alpha;
453 const Real theta = 1.;
463 for (
const auto & elem :
mesh.active_local_element_ptr_range())
473 dof_map.
dof_indices (elem, dof_indices_alpha, alpha_var);
475 const unsigned int n_dofs = dof_indices.size();
476 const unsigned int n_u_dofs = dof_indices_u.size();
477 const unsigned int n_v_dofs = dof_indices_v.size();
478 const unsigned int n_p_dofs = dof_indices_p.size();
484 fe_vel->reinit (elem);
485 fe_pres->reinit (elem);
493 Ke.
resize (n_dofs, n_dofs);
509 Kuu.reposition (u_var*n_u_dofs, u_var*n_u_dofs, n_u_dofs, n_u_dofs);
510 Kuv.reposition (u_var*n_u_dofs, v_var*n_u_dofs, n_u_dofs, n_v_dofs);
511 Kup.reposition (u_var*n_u_dofs, p_var*n_u_dofs, n_u_dofs, n_p_dofs);
513 Kvu.reposition (v_var*n_v_dofs, u_var*n_v_dofs, n_v_dofs, n_u_dofs);
514 Kvv.reposition (v_var*n_v_dofs, v_var*n_v_dofs, n_v_dofs, n_v_dofs);
515 Kvp.reposition (v_var*n_v_dofs, p_var*n_v_dofs, n_v_dofs, n_p_dofs);
517 Kpu.reposition (p_var*n_u_dofs, u_var*n_u_dofs, n_p_dofs, n_u_dofs);
518 Kpv.reposition (p_var*n_u_dofs, v_var*n_u_dofs, n_p_dofs, n_v_dofs);
519 Kpp.
reposition (p_var*n_u_dofs, p_var*n_u_dofs, n_p_dofs, n_p_dofs);
522 Kp_alpha.
reposition (p_var*n_u_dofs, p_var*n_u_dofs+n_p_dofs, n_p_dofs, 1);
523 Kalpha_p.reposition (p_var*n_u_dofs+n_p_dofs, p_var*n_u_dofs, 1, n_p_dofs);
526 Fu.reposition (u_var*n_u_dofs, n_u_dofs);
527 Fv.reposition (v_var*n_u_dofs, n_v_dofs);
536 for (
unsigned int qp=0; qp<qrule.n_points(); qp++)
539 Number u = 0., u_old = 0.;
540 Number v = 0., v_old = 0.;
547 for (
unsigned int l=0; l<n_u_dofs; l++)
550 u_old += phi[l][qp]*navier_stokes_system.
old_solution (dof_indices_u[l]);
551 v_old += phi[l][qp]*navier_stokes_system.
old_solution (dof_indices_v[l]);
556 u += phi[l][qp]*navier_stokes_system.current_solution (dof_indices_u[l]);
557 v += phi[l][qp]*navier_stokes_system.current_solution (dof_indices_v[l]);
558 grad_u.
add_scaled (dphi[l][qp], navier_stokes_system.current_solution (dof_indices_u[l]));
559 grad_v.
add_scaled (dphi[l][qp], navier_stokes_system.current_solution (dof_indices_v[l]));
563 for (
unsigned int l=0; l<n_p_dofs; l++)
564 p_old += psi[l][qp]*navier_stokes_system.
old_solution (dof_indices_p[l]);
570 const Number u_x = grad_u(0);
571 const Number u_y = grad_u(1);
572 const Number v_x = grad_v(0);
573 const Number v_y = grad_v(1);
578 for (
unsigned int i=0; i<n_u_dofs; i++)
580 Fu(i) += JxW[qp]*(u_old*phi[i][qp] -
581 (1.-theta)*dt*(U_old*grad_u_old)*phi[i][qp] +
582 (1.-theta)*dt*p_old*dphi[i][qp](0) -
583 (1.-theta)*dt*nu*(grad_u_old*dphi[i][qp]) +
584 theta*dt*(U*grad_u)*phi[i][qp]);
587 Fv(i) += JxW[qp]*(v_old*phi[i][qp] -
588 (1.-theta)*dt*(U_old*grad_v_old)*phi[i][qp] +
589 (1.-theta)*dt*p_old*dphi[i][qp](1) -
590 (1.-theta)*dt*nu*(grad_v_old*dphi[i][qp]) +
591 theta*dt*(U*grad_v)*phi[i][qp]);
598 for (
unsigned int j=0; j<n_u_dofs; j++)
600 Kuu(i,j) += JxW[qp]*(phi[i][qp]*phi[j][qp] +
601 theta*dt*nu*(dphi[i][qp]*dphi[j][qp]) +
602 theta*dt*(U*dphi[j][qp])*phi[i][qp] +
603 theta*dt*u_x*phi[i][qp]*phi[j][qp]);
605 Kuv(i,j) += JxW[qp]*theta*dt*u_y*phi[i][qp]*phi[j][qp];
607 Kvv(i,j) += JxW[qp]*(phi[i][qp]*phi[j][qp] +
608 theta*dt*nu*(dphi[i][qp]*dphi[j][qp]) +
609 theta*dt*(U*dphi[j][qp])*phi[i][qp] +
610 theta*dt*v_y*phi[i][qp]*phi[j][qp]);
612 Kvu(i,j) += JxW[qp]*theta*dt*v_x*phi[i][qp]*phi[j][qp];
616 for (
unsigned int j=0; j<n_p_dofs; j++)
618 Kup(i,j) += JxW[qp]*(-theta*dt*psi[j][qp]*dphi[i][qp](0));
619 Kvp(i,j) += JxW[qp]*(-theta*dt*psi[j][qp]*dphi[i][qp](1));
627 for (
unsigned int i=0; i<n_p_dofs; i++)
629 Kp_alpha(i,0) += JxW[qp]*psi[i][qp];
630 Kalpha_p(0,i) += JxW[qp]*psi[i][qp];
631 for (
unsigned int j=0; j<n_u_dofs; j++)
633 Kpu(i,j) += JxW[qp]*psi[i][qp]*dphi[j][qp](0);
634 Kpv(i,j) += JxW[qp]*psi[i][qp]*dphi[j][qp](1);
649 navier_stokes_system.rhs->add_vector (Fe, dof_indices);
654 navier_stokes_system.rhs->add(navier_stokes_system.rhs->size()-1, 10.);
664 #ifdef LIBMESH_ENABLE_DIRICHLET 666 u_var = system.variable_number(
"vel_x"),
667 v_var = system.variable_number(
"vel_y");
670 DofMap & dof_map = system.get_dof_map();
class FEType hides (possibly multiple) FEFamily and approximation orders, thereby enabling specialize...
T command_line_next(std::string name, T default_value)
Use GetPot's search()/next() functions to get following arguments from the command line...
This is the EquationSystems class.
void pop(const char *label, const char *header="")
Pop the event label off the stack, resuming any lower event.
void add_scaled(const TypeVector< T2 > &, const T &)
Add a scaled value to this vector without creating a temporary.
ConstFunction that simply returns 0.
void dof_indices(const Elem *const elem, std::vector< dof_id_type > &di) const
Fills the vector di with the global degree of freedom indices for the element.
static constexpr Real TOLERANCE
The ExodusII_IO class implements reading meshes in the ExodusII file format from Sandia National Labs...
void resize(const unsigned int n)
Resize the vector.
void print_info(std::ostream &os=libMesh::out) const
Prints information about the equation systems, by default to libMesh::out.
This class allows one to associate Dirichlet boundary values with a given set of mesh boundary ids an...
Manages storage and variables for transient systems.
Order default_quadrature_order() const
This class defines a vector in LIBMESH_DIM dimensional Real or Complex space.
The LibMeshInit class, when constructed, initializes the dependent libraries (e.g.
The libMesh namespace provides an interface to certain functionality in the library.
const T_sys & get_system(std::string_view name) const
NumericVector< Number > * old_local_solution
All the values I need to compute my contribution to the simulation at hand.
This is the MeshBase class.
Defines a dense subvector for use in finite element computations.
void set_lid_driven_bcs(TransientLinearImplicitSystem &system)
SolverPackage default_solver_package()
The PerfLog class allows monitoring of specific events.
This class handles the numbering of degrees of freedom on a mesh.
void reposition(const unsigned int ioff, const unsigned int n)
Changes the location of the subvector in the parent vector.
void libmesh_ignore(const Args &...)
virtual void add_matrix(const DenseMatrix< T > &dm, const std::vector< numeric_index_type > &rows, const std::vector< numeric_index_type > &cols)=0
Add the full matrix dm to the SparseMatrix.
const T & get(std::string_view) const
void print_info(std::ostream &os=libMesh::out, const unsigned int verbosity=0, const bool global=true) const
Prints relevant information about the mesh.
Defines a dense submatrix for use in Finite Element-type computations.
void init(triangulateio &t)
Initializes the fields of t to nullptr/0 as necessary.
static std::unique_ptr< FEGenericBase > build(const unsigned int dim, const FEType &type)
Builds a specific finite element type.
void push(const char *label, const char *header="")
Push the event label onto the stack, pausing any active event.
void assemble_stokes(EquationSystems &es, const std::string &system_name)
Number old_solution(const dof_id_type global_dof_number) const
void heterogenously_constrain_element_matrix_and_vector(DenseMatrix< Number > &matrix, DenseVector< Number > &rhs, std::vector< dof_id_type > &elem_dofs, bool asymmetric_constraint_rows=true, int qoi_index=-1) const
int main(int argc, char **argv)
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
T & set(const std::string &)
const MeshBase & get_mesh() const
void resize(const unsigned int new_m, const unsigned int new_n)
Resizes the matrix to the specified size and calls zero().
This class implements specific orders of Gauss quadrature.
unsigned int mesh_dimension() const
Parameters parameters
Data structure holding arbitrary parameters.
Function that returns a single value that never changes.
void add_dirichlet_boundary(const DirichletBoundary &dirichlet_boundary)
Adds a copy of the specified Dirichlet boundary to the system.
void reposition(const unsigned int ioff, const unsigned int joff, const unsigned int new_m, const unsigned int new_n)
Changes the location of the submatrix in the parent matrix.
virtual void init()
Initialize all the systems.
virtual System & add_system(std::string_view system_type, std::string_view name)
Add the system of type system_type named name to the systems array.
The Mesh class is a thin wrapper, around the ReplicatedMesh class by default.