libMesh
systems_of_equations_ex3.C
Go to the documentation of this file.
1 // The libMesh Finite Element Library.
2 // Copyright (C) 2002-2024 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>Systems Example 3 - Navier-Stokes with SCALAR Lagrange Multiplier</h1>
21 // \author David Knezevic
22 // \date 2010
23 //
24 // This example shows how the transient Navier-Stokes problem from
25 // systems_of_equations_ex2 can be solved using a scalar Lagrange multiplier
26 // formulation to constrain the integral of the pressure variable,
27 // rather than pinning the pressure at a single point.
28 
29 // C++ include files that we need
30 #include <iostream>
31 #include <algorithm>
32 #include <sstream>
33 #include <math.h>
34 
35 // Basic include file needed for the mesh functionality.
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"
58 
59 // For systems of equations the DenseSubMatrix
60 // and DenseSubVector provide convenient ways for
61 // assembling the element matrix and vector on a
62 // component-by-component basis.
63 #include "libmesh/dense_submatrix.h"
64 #include "libmesh/dense_subvector.h"
65 
66 // The definition of a geometric element
67 #include "libmesh/elem.h"
68 
69 // Bring in everything from the libMesh namespace
70 using namespace libMesh;
71 
72 // Function prototype. This function will assemble the system
73 // matrix and right-hand-side.
75  const std::string & system_name);
76 
77 // Function which sets Dirichlet BCs for the lid-driven cavity.
79 
80 // The main program.
81 int main (int argc, char ** argv)
82 {
83  // Initialize libMesh.
84  LibMeshInit init (argc, argv);
85 
86  // This example requires a linear solver package.
87  libmesh_example_requires(libMesh::default_solver_package() != INVALID_SOLVER_PACKAGE,
88  "--enable-petsc, --enable-trilinos, or --enable-eigen");
89 
90  // Skip this 2D example if libMesh was compiled as 1D-only.
91  libmesh_example_requires(2 <= LIBMESH_DIM, "2D support");
92 
93  // We use Dirichlet boundary conditions here
94 #ifndef LIBMESH_ENABLE_DIRICHLET
95  libmesh_example_requires(false, "--enable-dirichlet");
96 #endif
97 
98  // This example NaNs with the Trilinos solvers
99  libmesh_example_requires(libMesh::default_solver_package() != TRILINOS_SOLVERS, "--enable-petsc or --enable-laspack");
100 
101  // Create a mesh, with dimension to be overridden later, distributed
102  // across the default MPI communicator.
103  Mesh mesh(init.comm());
104 
105  // Get the mesh size from the command line.
106  const int n_elem = libMesh::command_line_next("-n_elem", 20);
107 
108  // Use the MeshTools::Generation mesh generator to create a uniform
109  // 2D grid on the square [-1,1]^2. We instruct the mesh generator
110  // to build a mesh of 8x8 Quad9 elements in 2D. Building these
111  // higher-order elements allows us to use higher-order
112  // approximation, as in example 3.
114  n_elem, n_elem,
115  0., 1.,
116  0., 1.,
117  QUAD9);
118 
119  // Print information about the mesh to the screen.
120  mesh.print_info();
121 
122  // Create an equation systems object.
123  EquationSystems equation_systems (mesh);
124 
125  // Declare the system and its variables.
126  // Creates a transient system named "Navier-Stokes"
128  equation_systems.add_system<TransientLinearImplicitSystem> ("Navier-Stokes");
129 
130  // Add the variables "vel_x" & "vel_y" to "Navier-Stokes". They
131  // will be approximated using second-order approximation.
132  system.add_variable ("vel_x", SECOND);
133  system.add_variable ("vel_y", SECOND);
134 
135  // Add the variable "p" to "Navier-Stokes". This will
136  // be approximated with a first-order basis,
137  // providing an LBB-stable pressure-velocity pair.
138  system.add_variable ("p", FIRST);
139 
140  // Add a scalar Lagrange multiplier to constrain the
141  // pressure to have zero mean.
142  system.add_variable ("alpha", FIRST, SCALAR);
143 
144  // Give the system a pointer to the matrix assembly
145  // function.
146  system.attach_assemble_function (assemble_stokes);
147 
148  // Set Dirichlet boundary conditions.
149  set_lid_driven_bcs(system);
150 
151  // Initialize the data structures for the equation system.
152  equation_systems.init ();
153 
154  // Prints information about the system to the screen.
155  equation_systems.print_info();
156 
157  // Create a performance-logging object for this example
158  PerfLog perf_log("Systems Example 3");
159 
160  // Get a reference to the Stokes system to use later.
161  TransientLinearImplicitSystem & navier_stokes_system =
162  equation_systems.get_system<TransientLinearImplicitSystem>("Navier-Stokes");
163 
164  // Now we begin the timestep loop to compute the time-accurate
165  // solution of the equations.
166  const Real dt = 0.1;
167  navier_stokes_system.time = 0.0;
168  const unsigned int n_timesteps = 15;
169 
170  // The number of steps and the stopping criterion are also required
171  // for the nonlinear iterations.
172  const unsigned int n_nonlinear_steps = 15;
173  const Real nonlinear_tolerance = TOLERANCE*10;
174 
175  // We also set a standard linear solver flag in the EquationSystems object
176  // which controls the maximum number of linear solver iterations allowed.
177  const int max_iter = libMesh::command_line_next("-max_iter", 250);
178 
179  equation_systems.parameters.set<unsigned int>("linear solver maximum iterations") = max_iter;
180 
181  // Tell the system of equations what the timestep is by using
182  // the set_parameter function. The matrix assembly routine can
183  // then reference this parameter.
184  equation_systems.parameters.set<Real> ("dt") = dt;
185 
186  // The kinematic viscosity, nu = mu/rho, units of length**2/time.
187  equation_systems.parameters.set<Real> ("nu") = .007;
188 
189  // The first thing to do is to get a copy of the solution at
190  // the current nonlinear iteration. This value will be used to
191  // determine if we can exit the nonlinear loop.
192  std::unique_ptr<NumericVector<Number>>
193  last_nonlinear_soln (navier_stokes_system.solution->clone());
194 
195 #ifdef LIBMESH_HAVE_EXODUS_API
196  // Since we are not doing adaptivity, write all solutions to a single Exodus file.
197  ExodusII_IO exo_io(mesh);
198 
199  // Write out the initial condition
200  exo_io.write_equation_systems ("out.e", equation_systems);
201 #endif
202 
203  for (unsigned int t_step=1; t_step<=n_timesteps; ++t_step)
204  {
205  // Increment the time counter, set the time and the
206  // time step size as parameters in the EquationSystem.
207  navier_stokes_system.time += dt;
208 
209  // A pretty update message
210  libMesh::out << "\n\n*** Solving time step "
211  << t_step
212  << ", time = "
213  << navier_stokes_system.time
214  << " ***"
215  << std::endl;
216 
217  // Now we need to update the solution vector from the
218  // previous time step. This is done directly through
219  // the reference to the Stokes system.
220  *navier_stokes_system.old_local_solution = *navier_stokes_system.current_local_solution;
221 
222  // At the beginning of each solve, reset the linear solver tolerance
223  // to a "reasonable" starting value.
224  const Real initial_linear_solver_tol = 1.e-6;
225  equation_systems.parameters.set<Real> ("linear solver tolerance") = initial_linear_solver_tol;
226 
227  // We'll set this flag when convergence is (hopefully) achieved.
228  bool converged = false;
229 
230  // Now we begin the nonlinear loop
231  for (unsigned int l=0; l<n_nonlinear_steps; ++l)
232  {
233  // Update the nonlinear solution.
234  last_nonlinear_soln->zero();
235  last_nonlinear_soln->add(*navier_stokes_system.solution);
236 
237  // Assemble & solve the linear system.
238  perf_log.push("linear solve");
239  equation_systems.get_system("Navier-Stokes").solve();
240  perf_log.pop("linear solve");
241 
242  // Compute the difference between this solution and the last
243  // nonlinear iterate.
244  last_nonlinear_soln->add (-1., *navier_stokes_system.solution);
245 
246  // Close the vector before computing its norm
247  last_nonlinear_soln->close();
248 
249  // Compute the l2 norm of the difference
250  const Real norm_delta = last_nonlinear_soln->l2_norm();
251 
252  // How many iterations were required to solve the linear system?
253  const unsigned int n_linear_iterations = navier_stokes_system.n_linear_iterations();
254 
255  // What was the final residual of the linear system?
256  const Real final_linear_residual = navier_stokes_system.final_linear_residual();
257 
258  // If the solver did no work (sometimes -ksp_converged_reason
259  // says "Linear solve converged due to CONVERGED_RTOL
260  // iterations 0") but the nonlinear residual norm is above
261  // the tolerance, we need to pick an even lower linear
262  // solver tolerance and try again. Note that the tolerance
263  // is relative to the norm of the RHS, which for this
264  // particular problem does not go to zero, since we are
265  // solving for the full solution rather than the update.
266  //
267  // Similarly, if the solver did no work and this is the 0th
268  // nonlinear step, it means that the delta between solutions
269  // is being inaccurately measured as "0" since the solution
270  // did not change. Decrease the tolerance and try again.
271  if (n_linear_iterations == 0 &&
272  (navier_stokes_system.final_linear_residual() >= nonlinear_tolerance || l==0))
273  {
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;
276  continue;
277  }
278 
279  // Print out convergence information for the linear and
280  // nonlinear iterations.
281  libMesh::out << "Linear solver converged at step: "
282  << n_linear_iterations
283  << ", final residual: "
284  << final_linear_residual
285  << " Nonlinear convergence: ||u - u_old|| = "
286  << norm_delta
287  << std::endl;
288 
289  // Terminate the solution iteration if the difference between
290  // this nonlinear iterate and the last is sufficiently small, AND
291  // if the most recent linear system was solved to a sufficient tolerance.
292  if ((norm_delta < nonlinear_tolerance) &&
293  (navier_stokes_system.final_linear_residual() < nonlinear_tolerance))
294  {
295  libMesh::out << " Nonlinear solver converged at step "
296  << l
297  << std::endl;
298  converged = true;
299  break;
300  }
301 
302  // Otherwise, decrease the linear system tolerance. For the inexact Newton
303  // method, the linear solver tolerance needs to decrease as we get closer to
304  // the solution to ensure quadratic convergence. The new linear solver tolerance
305  // is chosen (heuristically) as the square of the previous linear system residual norm.
306  //Real flr2 = final_linear_residual*final_linear_residual;
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;
309  } // end nonlinear loop
310 
311  // Don't keep going if we failed to converge.
312  libmesh_error_msg_if(!converged, "Error: Newton iterations failed to converge!");
313 
314 #ifdef LIBMESH_HAVE_EXODUS_API
315  // Write out every nth timestep to file.
316  const unsigned int write_interval = 1;
317 
318  if ((t_step+1)%write_interval == 0)
319  {
320  exo_io.write_timestep("out.e",
321  equation_systems,
322  t_step+1, // we're off by one since we wrote the IC and the Exodus numbering is 1-based.
323  navier_stokes_system.time);
324  }
325 #endif // #ifdef LIBMESH_HAVE_EXODUS_API
326  } // end timestep loop.
327 
328  // All done.
329  return 0;
330 }
331 
332 
333 
334 
335 
336 
337 // The matrix assembly function to be called at each time step to
338 // prepare for the linear solve.
340  const std::string & libmesh_dbg_var(system_name))
341 {
342  // It is a good idea to make sure we are assembling
343  // the proper system.
344  libmesh_assert_equal_to (system_name, "Navier-Stokes");
345 
346 #if LIBMESH_DIM > 1
347  // Get a constant reference to the mesh object.
348  const MeshBase & mesh = es.get_mesh();
349 
350  // The dimension that we are running
351  const unsigned int dim = mesh.mesh_dimension();
352 
353  // Get a reference to the Stokes system object.
354  TransientLinearImplicitSystem & navier_stokes_system =
355  es.get_system<TransientLinearImplicitSystem> ("Navier-Stokes");
356 
357  // Numeric ids corresponding to each variable in the system
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");
362 
363  // Get the Finite Element type for "vel_x". Note this will be
364  // the same as the type for "vel_y".
365  FEType fe_vel_type = navier_stokes_system.variable_type(u_var);
366 
367  // Get the Finite Element type for "p".
368  FEType fe_pres_type = navier_stokes_system.variable_type(p_var);
369 
370  // Build a Finite Element object of the specified type for
371  // the velocity variables.
372  std::unique_ptr<FEBase> fe_vel (FEBase::build(dim, fe_vel_type));
373 
374  // Build a Finite Element object of the specified type for
375  // the pressure variables.
376  std::unique_ptr<FEBase> fe_pres (FEBase::build(dim, fe_pres_type));
377 
378  // A Gauss quadrature rule for numerical integration.
379  // Let the FEType object decide what order rule is appropriate.
380  QGauss qrule (dim, fe_vel_type.default_quadrature_order());
381 
382  // Tell the finite element objects to use our quadrature rule.
383  fe_vel->attach_quadrature_rule (&qrule);
384  fe_pres->attach_quadrature_rule (&qrule);
385 
386  // Here we define some references to cell-specific data that
387  // will be used to assemble the linear system.
388  //
389  // The element Jacobian * quadrature weight at each integration point.
390  const std::vector<Real> & JxW = fe_vel->get_JxW();
391 
392  // The element shape functions evaluated at the quadrature points.
393  const std::vector<std::vector<Real>> & phi = fe_vel->get_phi();
394 
395  // The element shape function gradients for the velocity
396  // variables evaluated at the quadrature points.
397  const std::vector<std::vector<RealGradient>> & dphi = fe_vel->get_dphi();
398 
399  // The element shape functions for the pressure variable
400  // evaluated at the quadrature points.
401  const std::vector<std::vector<Real>> & psi = fe_pres->get_phi();
402 
403  // The value of the linear shape function gradients at the quadrature points
404  // const std::vector<std::vector<RealGradient>> & dpsi = fe_pres->get_dphi();
405 
406  // A reference to the DofMap object for this system. The DofMap
407  // object handles the index translation from node and element numbers
408  // to degree of freedom numbers. We will talk more about the DofMap
409  // in future examples.
410  const DofMap & dof_map = navier_stokes_system.get_dof_map();
411 
412  // Define data structures to contain the element matrix
413  // and right-hand-side vector contribution. Following
414  // basic finite element terminology we will denote these
415  // "Ke" and "Fe".
418 
420  Kuu(Ke), Kuv(Ke), Kup(Ke),
421  Kvu(Ke), Kvv(Ke), Kvp(Ke),
422  Kpu(Ke), Kpv(Ke), Kpp(Ke);
423  DenseSubMatrix<Number> Kalpha_p(Ke), Kp_alpha(Ke);
424 
426  Fu(Fe),
427  Fv(Fe),
428  Fp(Fe);
429 
430  // This vector will hold the degree of freedom indices for
431  // the element. These define where in the global system
432  // the element degrees of freedom get mapped.
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;
438 
439  // The global system matrix
440  SparseMatrix<Number> & matrix = navier_stokes_system.get_system_matrix();
441 
442  // Find out what the timestep size parameter is from the system, and
443  // the value of theta for the theta method. We use implicit Euler (theta=1)
444  // for this simulation even though it is only first-order accurate in time.
445  // The reason for this decision is that the second-order Crank-Nicolson
446  // method is notoriously oscillatory for problems with discontinuous
447  // initial data such as the lid-driven cavity. Therefore,
448  // we sacrifice accuracy in time for stability, but since the solution
449  // reaches steady state relatively quickly we can afford to take small
450  // timesteps. If you monitor the initial nonlinear residual for this
451  // simulation, you should see that it is monotonically decreasing in time.
452  const Real dt = es.parameters.get<Real>("dt");
453  const Real theta = 1.;
454 
455  // The kinematic viscosity, multiplies the "viscous" terms.
456  const Real nu = es.parameters.get<Real>("nu");
457 
458  // Now we will loop over all the elements in the mesh that
459  // live on the local processor. We will compute the element
460  // matrix and right-hand-side contribution. Since the mesh
461  // will be refined we want to only consider the ACTIVE elements,
462  // hence we use a variant of the active_elem_iterator.
463  for (const auto & elem : mesh.active_local_element_ptr_range())
464  {
465  // Get the degree of freedom indices for the
466  // current element. These define where in the global
467  // matrix and right-hand-side this element will
468  // contribute to.
469  dof_map.dof_indices (elem, dof_indices);
470  dof_map.dof_indices (elem, dof_indices_u, u_var);
471  dof_map.dof_indices (elem, dof_indices_v, v_var);
472  dof_map.dof_indices (elem, dof_indices_p, p_var);
473  dof_map.dof_indices (elem, dof_indices_alpha, alpha_var);
474 
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();
479 
480  // Compute the element-specific data for the current
481  // element. This involves computing the location of the
482  // quadrature points (q_point) and the shape functions
483  // (phi, dphi) for the current element.
484  fe_vel->reinit (elem);
485  fe_pres->reinit (elem);
486 
487  // Zero the element matrix and right-hand side before
488  // summing them. We use the resize member here because
489  // the number of degrees of freedom might have changed from
490  // the last element. Note that this will be the case if the
491  // element type is different (i.e. the last element was a
492  // triangle, now we are on a quadrilateral).
493  Ke.resize (n_dofs, n_dofs);
494  Fe.resize (n_dofs);
495 
496  // Reposition the submatrices... The idea is this:
497  //
498  // - - - -
499  // | Kuu Kuv Kup | | Fu |
500  // Ke = | Kvu Kvv Kvp |; Fe = | Fv |
501  // | Kpu Kpv Kpp | | Fp |
502  // - - - -
503  //
504  // The DenseSubMatrix.reposition () member takes the
505  // (row_offset, column_offset, row_size, column_size).
506  //
507  // Similarly, the DenseSubVector.reposition () member
508  // takes the (row_offset, row_size)
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);
512 
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);
516 
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);
520 
521  // Also, add a row and a column to constrain the pressure
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);
524 
525 
526  Fu.reposition (u_var*n_u_dofs, n_u_dofs);
527  Fv.reposition (v_var*n_u_dofs, n_v_dofs);
528  Fp.reposition (p_var*n_u_dofs, n_p_dofs);
529 
530  // Now we will build the element matrix and right-hand-side.
531  // Constructing the RHS requires the solution and its
532  // gradient from the previous timestep. This must be
533  // calculated at each quadrature point by summing the
534  // solution degree-of-freedom values by the appropriate
535  // weight functions.
536  for (unsigned int qp=0; qp<qrule.n_points(); qp++)
537  {
538  // Values to hold the solution & its gradient at the previous timestep.
539  Number u = 0., u_old = 0.;
540  Number v = 0., v_old = 0.;
541  Number p_old = 0.;
542  Gradient grad_u, grad_u_old;
543  Gradient grad_v, grad_v_old;
544 
545  // Compute the velocity & its gradient from the previous timestep
546  // and the old Newton iterate.
547  for (unsigned int l=0; l<n_u_dofs; l++)
548  {
549  // From the old timestep:
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]);
552  grad_u_old.add_scaled (dphi[l][qp], navier_stokes_system.old_solution (dof_indices_u[l]));
553  grad_v_old.add_scaled (dphi[l][qp], navier_stokes_system.old_solution (dof_indices_v[l]));
554 
555  // From the previous Newton iterate:
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]));
560  }
561 
562  // Compute the old pressure value at this quadrature point.
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]);
565 
566  // Definitions for convenience. It is sometimes simpler to do a
567  // dot product if you have the full vector at your disposal.
568  const NumberVectorValue U_old (u_old, v_old);
569  const NumberVectorValue U (u, v);
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);
574 
575  // First, an i-loop over the velocity degrees of freedom.
576  // We know that n_u_dofs == n_v_dofs so we can compute contributions
577  // for both at the same time.
578  for (unsigned int i=0; i<n_u_dofs; i++)
579  {
580  Fu(i) += JxW[qp]*(u_old*phi[i][qp] - // mass-matrix term
581  (1.-theta)*dt*(U_old*grad_u_old)*phi[i][qp] + // convection term
582  (1.-theta)*dt*p_old*dphi[i][qp](0) - // pressure term on rhs
583  (1.-theta)*dt*nu*(grad_u_old*dphi[i][qp]) + // diffusion term on rhs
584  theta*dt*(U*grad_u)*phi[i][qp]); // Newton term
585 
586 
587  Fv(i) += JxW[qp]*(v_old*phi[i][qp] - // mass-matrix term
588  (1.-theta)*dt*(U_old*grad_v_old)*phi[i][qp] + // convection term
589  (1.-theta)*dt*p_old*dphi[i][qp](1) - // pressure term on rhs
590  (1.-theta)*dt*nu*(grad_v_old*dphi[i][qp]) + // diffusion term on rhs
591  theta*dt*(U*grad_v)*phi[i][qp]); // Newton term
592 
593 
594  // Note that the Fp block is identically zero unless we are using
595  // some kind of artificial compressibility scheme...
596 
597  // Matrix contributions for the uu and vv couplings.
598  for (unsigned int j=0; j<n_u_dofs; j++)
599  {
600  Kuu(i,j) += JxW[qp]*(phi[i][qp]*phi[j][qp] + // mass matrix term
601  theta*dt*nu*(dphi[i][qp]*dphi[j][qp]) + // diffusion term
602  theta*dt*(U*dphi[j][qp])*phi[i][qp] + // convection term
603  theta*dt*u_x*phi[i][qp]*phi[j][qp]); // Newton term
604 
605  Kuv(i,j) += JxW[qp]*theta*dt*u_y*phi[i][qp]*phi[j][qp]; // Newton term
606 
607  Kvv(i,j) += JxW[qp]*(phi[i][qp]*phi[j][qp] + // mass matrix term
608  theta*dt*nu*(dphi[i][qp]*dphi[j][qp]) + // diffusion term
609  theta*dt*(U*dphi[j][qp])*phi[i][qp] + // convection term
610  theta*dt*v_y*phi[i][qp]*phi[j][qp]); // Newton term
611 
612  Kvu(i,j) += JxW[qp]*theta*dt*v_x*phi[i][qp]*phi[j][qp]; // Newton term
613  }
614 
615  // Matrix contributions for the up and vp couplings.
616  for (unsigned int j=0; j<n_p_dofs; j++)
617  {
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));
620  }
621  }
622 
623  // Now an i-loop over the pressure degrees of freedom. This code computes
624  // the matrix entries due to the continuity equation. Note: To maintain a
625  // symmetric matrix, we may (or may not) multiply the continuity equation by
626  // negative one. Here we do not.
627  for (unsigned int i=0; i<n_p_dofs; i++)
628  {
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++)
632  {
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);
635  }
636  }
637  } // end of the quadrature point qp-loop
638 
639  // Since we're using heterogeneous DirichletBoundary objects for
640  // the boundary conditions, we need to call a specific function
641  // to constrain the element stiffness matrix.
642  dof_map.heterogenously_constrain_element_matrix_and_vector (Ke, Fe, dof_indices);
643 
644  // The element matrix and right-hand-side are now built
645  // for this element. Add them to the global matrix and
646  // right-hand-side vector. The SparseMatrix::add_matrix()
647  // and NumericVector::add_vector() members do this for us.
648  matrix.add_matrix (Ke, dof_indices);
649  navier_stokes_system.rhs->add_vector (Fe, dof_indices);
650  } // end of element loop
651 
652  // We can set the mean of the pressure by setting Falpha. Typically
653  // a value of zero is chosen, but the value should be arbitrary.
654  navier_stokes_system.rhs->add(navier_stokes_system.rhs->size()-1, 10.);
655 #else
656  libmesh_ignore(es);
657 #endif
658 }
659 
660 
661 
663 {
664 #ifdef LIBMESH_ENABLE_DIRICHLET
665  unsigned short int
666  u_var = system.variable_number("vel_x"),
667  v_var = system.variable_number("vel_y");
668 
669  // Get a convenient reference to the System's DofMap
670  DofMap & dof_map = system.get_dof_map();
671 
672  // u=v=0 on bottom, left, right (boundaries 0,1,3)
673  dof_map.add_dirichlet_boundary(DirichletBoundary({0,1,3}, {u_var, v_var},
675  // u=1 on top (boundary 2)
676  dof_map.add_dirichlet_boundary(DirichletBoundary({2}, {u_var},
677  ConstFunction<Number>(1.)));
678  // v=0 on top
679  dof_map.add_dirichlet_boundary(DirichletBoundary({2}, {v_var},
681 #else
682  libmesh_ignore(system);
683 #endif
684 }
class FEType hides (possibly multiple) FEFamily and approximation orders, thereby enabling specialize...
Definition: fe_type.h:182
T command_line_next(std::string name, T default_value)
Use GetPot&#39;s search()/next() functions to get following arguments from the command line...
Definition: libmesh.C:1011
This is the EquationSystems class.
void pop(const char *label, const char *header="")
Pop the event label off the stack, resuming any lower event.
Definition: perf_log.C:168
void add_scaled(const TypeVector< T2 > &, const T &)
Add a scaled value to this vector without creating a temporary.
Definition: type_vector.h:651
ConstFunction that simply returns 0.
Definition: zero_function.h:38
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.
Definition: dof_map.C:1992
dof_id_type n_elem(const MeshBase::const_element_iterator &begin, const MeshBase::const_element_iterator &end)
Count up the number of elements of a specific type (as defined by an iterator range).
Definition: mesh_tools.C:850
static constexpr Real TOLERANCE
unsigned int dim
The ExodusII_IO class implements reading meshes in the ExodusII file format from Sandia National Labs...
Definition: exodusII_io.h:52
void resize(const unsigned int n)
Resize the vector.
Definition: dense_vector.h:374
void print_info(std::ostream &os=libMesh::out) const
Prints information about the equation systems, by default to libMesh::out.
MeshBase & mesh
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
Definition: fe_type.h:357
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.
This class defines a vector in LIBMESH_DIM dimensional Real or Complex space.
The LibMeshInit class, when constructed, initializes the dependent libraries (e.g.
Definition: libmesh.h:90
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.
Definition: mesh_base.h:74
Defines a dense subvector for use in finite element computations.
void set_lid_driven_bcs(TransientLinearImplicitSystem &system)
SolverPackage default_solver_package()
Definition: libmesh.C:1050
The PerfLog class allows monitoring of specific events.
Definition: perf_log.h:145
This class handles the numbering of degrees of freedom on a mesh.
Definition: dof_map.h:169
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
Definition: parameters.h:426
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
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.
Definition: perf_log.C:138
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
Definition: dof_map.h:1237
int main(int argc, char **argv)
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
T & set(const std::string &)
Definition: parameters.h:469
OStreamProxy out
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().
Definition: dense_matrix.h:895
This class implements specific orders of Gauss quadrature.
unsigned int mesh_dimension() const
Definition: mesh_base.C:324
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.
Definition: mesh.h:50