libMesh
systems_of_equations_ex2.C
Go to the documentation of this file.
1 // The libMesh Finite Element Library.
2 // Copyright (C) 2002-2017 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 2 - Unsteady Nonlinear Navier-Stokes</h1>
21 // \author John W. Peterson
22 // \date 2004
23 //
24 // This example shows how a simple, unsteady, nonlinear system of equations
25 // can be solved in parallel. The system of equations are the familiar
26 // Navier-Stokes equations for low-speed incompressible fluid flow. This
27 // example introduces the concept of the inner nonlinear loop for each
28 // timestep, and requires a good deal of linear algebra number-crunching
29 // at each step. If you have a ExodusII viewer such as ParaView installed,
30 // the script movie.sh in this directory will also take appropriate screen
31 // shots of each of the solution files in the time sequence. These rgb files
32 // can then be animated with the "animate" utility of ImageMagick if it is
33 // installed on your system. On a PIII 1GHz machine in debug mode, this
34 // example takes a little over a minute to run. If you would like to see
35 // a more detailed time history, or compute more timesteps, that is certainly
36 // possible by changing the n_timesteps and dt variables below.
37 
38 // C++ include files that we need
39 #include <iostream>
40 #include <algorithm>
41 #include <sstream>
42 #include <math.h>
43 
44 // Basic include file needed for the mesh functionality.
45 #include "libmesh/libmesh.h"
46 #include "libmesh/mesh.h"
47 #include "libmesh/mesh_generation.h"
48 #include "libmesh/exodusII_io.h"
49 #include "libmesh/equation_systems.h"
50 #include "libmesh/fe.h"
51 #include "libmesh/quadrature_gauss.h"
52 #include "libmesh/dof_map.h"
53 #include "libmesh/sparse_matrix.h"
54 #include "libmesh/numeric_vector.h"
55 #include "libmesh/dense_matrix.h"
56 #include "libmesh/dense_vector.h"
57 #include "libmesh/linear_implicit_system.h"
58 #include "libmesh/transient_system.h"
59 #include "libmesh/perf_log.h"
60 #include "libmesh/boundary_info.h"
61 #include "libmesh/utility.h"
62 #include "libmesh/dirichlet_boundaries.h"
63 #include "libmesh/zero_function.h"
64 #include "libmesh/const_function.h"
65 #include "libmesh/parsed_function.h"
66 
67 // For systems of equations the DenseSubMatrix
68 // and DenseSubVector provide convenient ways for
69 // assembling the element matrix and vector on a
70 // component-by-component basis.
71 #include "libmesh/dense_submatrix.h"
72 #include "libmesh/dense_subvector.h"
73 
74 // The definition of a geometric element
75 #include "libmesh/elem.h"
76 
77 // Bring in everything from the libMesh namespace
78 using namespace libMesh;
79 
80 // Function prototype. This function will assemble the system
81 // matrix and right-hand-side.
83  const std::string & system_name);
84 
85 // Functions which set Dirichlet BCs corresponding to different problems.
89 
90 // The main program.
91 int main (int argc, char** argv)
92 {
93  // Initialize libMesh.
94  LibMeshInit init (argc, argv);
95 
96  // This example requires a linear solver package.
97  libmesh_example_requires(libMesh::default_solver_package() != INVALID_SOLVER_PACKAGE,
98  "--enable-petsc, --enable-trilinos, or --enable-eigen");
99 
100  // Skip this 2D example if libMesh was compiled as 1D-only.
101  libmesh_example_requires(2 <= LIBMESH_DIM, "2D support");
102 
103  // This example NaNs with the Eigen sparse linear solvers and
104  // Trilinos solvers, but should work OK with either PETSc or
105  // Laspack.
106  libmesh_example_requires(libMesh::default_solver_package() != EIGEN_SOLVERS, "--enable-petsc or --enable-laspack");
107  libmesh_example_requires(libMesh::default_solver_package() != TRILINOS_SOLVERS, "--enable-petsc or --enable-laspack");
108 
109  // Create a mesh, with dimension to be overridden later, distributed
110  // across the default MPI communicator.
111  Mesh mesh(init.comm());
112 
113  // Use the MeshTools::Generation mesh generator to create a uniform
114  // 2D grid on the square [-1,1]^2. We instruct the mesh generator
115  // to build a mesh of 8x8 Quad9 elements in 2D. Building these
116  // higher-order elements allows us to use higher-order polynomial
117  // approximations for the velocity.
119  20, 20,
120  0., 1.,
121  0., 1.,
122  QUAD9);
123 
124  // Print information about the mesh to the screen.
125  mesh.print_info();
126 
127  // Create an equation systems object.
128  EquationSystems equation_systems (mesh);
129 
130  // Declare the system and its variables.
131  // Creates a transient system named "Navier-Stokes"
133  equation_systems.add_system<TransientLinearImplicitSystem> ("Navier-Stokes");
134 
135  // Add the variables "vel_x" & "vel_y" to "Navier-Stokes". They
136  // will be approximated using second-order approximation.
137  system.add_variable ("vel_x", SECOND);
138  system.add_variable ("vel_y", SECOND);
139 
140  // Add the variable "p" to "Navier-Stokes". This will
141  // be approximated with a first-order basis,
142  // providing an LBB-stable pressure-velocity pair.
143  system.add_variable ("p", FIRST);
144 
145  // Give the system a pointer to the matrix assembly
146  // function.
147  system.attach_assemble_function (assemble_stokes);
148 
149  // Note: only pick one set of BCs!
150  set_lid_driven_bcs(system);
151  // set_stagnation_bcs(system);
152  // set_poiseuille_bcs(system);
153 
154  // Initialize the data structures for the equation system.
155  equation_systems.init ();
156 
157  // Prints information about the system to the screen.
158  equation_systems.print_info();
159 
160  // Create a performance-logging object for this example
161  PerfLog perf_log("Systems Example 2");
162 
163  // Get a reference to the Stokes system to use later.
164  TransientLinearImplicitSystem & navier_stokes_system =
165  equation_systems.get_system<TransientLinearImplicitSystem>("Navier-Stokes");
166 
167  // Now we begin the timestep loop to compute the time-accurate
168  // solution of the equations.
169  const Real dt = 0.1;
170  navier_stokes_system.time = 0.0;
171  const unsigned int n_timesteps = 15;
172 
173  // The number of steps and the stopping criterion are also required
174  // for the nonlinear iterations.
175  const unsigned int n_nonlinear_steps = 15;
176  const Real nonlinear_tolerance = 1.e-5;
177 
178  // We also set a standard linear solver flag in the EquationSystems object
179  // which controls the maximum number of linear solver iterations allowed.
180  equation_systems.parameters.set<unsigned int>("linear solver maximum iterations") = 250;
181 
182  // Tell the system of equations what the timestep is by using
183  // the set_parameter function. The matrix assembly routine can
184  // then reference this parameter.
185  equation_systems.parameters.set<Real> ("dt") = dt;
186 
187  // The kinematic viscosity, nu = mu/rho, units of length**2/time.
188  equation_systems.parameters.set<Real> ("nu") = .007;
189 
190  // The first thing to do is to get a copy of the solution at
191  // the current nonlinear iteration. This value will be used to
192  // determine if we can exit the nonlinear loop.
194  last_nonlinear_soln (navier_stokes_system.solution->clone());
195 
196  // Since we are not doing adaptivity, write all solutions to a single Exodus file.
197  ExodusII_IO exo_io(mesh);
198 
199  for (unsigned int t_step=1; t_step<=n_timesteps; ++t_step)
200  {
201  // Increment the time counter, set the time step size as
202  // a parameter in the EquationSystem.
203  navier_stokes_system.time += dt;
204 
205  // A pretty update message
206  libMesh::out << "\n\n*** Solving time step "
207  << t_step
208  << ", time = "
209  << navier_stokes_system.time
210  << " ***"
211  << std::endl;
212 
213  // Now we need to update the solution vector from the
214  // previous time step. This is done directly through
215  // the reference to the Stokes system.
216  *navier_stokes_system.old_local_solution = *navier_stokes_system.current_local_solution;
217 
218  // At the beginning of each solve, reset the linear solver tolerance
219  // to a "reasonable" starting value.
220  const Real initial_linear_solver_tol = 1.e-6;
221  equation_systems.parameters.set<Real> ("linear solver tolerance") = initial_linear_solver_tol;
222 
223  // We'll set this flag when convergence is (hopefully) achieved.
224  bool converged = false;
225 
226  // Now we begin the nonlinear loop
227  for (unsigned int l=0; l<n_nonlinear_steps; ++l)
228  {
229  // Update the nonlinear solution.
230  last_nonlinear_soln->zero();
231  last_nonlinear_soln->add(*navier_stokes_system.solution);
232 
233  // Assemble & solve the linear system.
234  perf_log.push("linear solve");
235  equation_systems.get_system("Navier-Stokes").solve();
236  perf_log.pop("linear solve");
237 
238  // Compute the difference between this solution and the last
239  // nonlinear iterate.
240  last_nonlinear_soln->add (-1., *navier_stokes_system.solution);
241 
242  // Close the vector before computing its norm
243  last_nonlinear_soln->close();
244 
245  // Compute the l2 norm of the difference
246  const Real norm_delta = last_nonlinear_soln->l2_norm();
247 
248  // How many iterations were required to solve the linear system?
249  const unsigned int n_linear_iterations = navier_stokes_system.n_linear_iterations();
250 
251  // What was the final residual of the linear system?
252  const Real final_linear_residual = navier_stokes_system.final_linear_residual();
253 
254  // If the solver did no work (sometimes -ksp_converged_reason
255  // says "Linear solve converged due to CONVERGED_RTOL
256  // iterations 0") but the nonlinear residual norm is above
257  // the tolerance, we need to pick an even lower linear
258  // solver tolerance and try again. Note that the tolerance
259  // is relative to the norm of the RHS, which for this
260  // particular problem does not go to zero, since we are
261  // solving for the full solution rather than the update.
262  //
263  // Similarly, if the solver did no work and this is the 0th
264  // nonlinear step, it means that the delta between solutions
265  // is being inaccurately measured as "0" since the solution
266  // did not change. Decrease the tolerance and try again.
267  if (n_linear_iterations == 0 &&
268  (navier_stokes_system.final_linear_residual() >= nonlinear_tolerance || l==0))
269  {
270  Real old_linear_solver_tolerance = equation_systems.parameters.get<Real> ("linear solver tolerance");
271  equation_systems.parameters.set<Real> ("linear solver tolerance") = 1.e-3 * old_linear_solver_tolerance;
272  continue;
273  }
274 
275  // Print out convergence information for the linear and
276  // nonlinear iterations.
277  libMesh::out << "Linear solver converged at step: "
278  << n_linear_iterations
279  << ", final residual: "
280  << final_linear_residual
281  << " Nonlinear convergence: ||u - u_old|| = "
282  << norm_delta
283  << std::endl;
284 
285  // Terminate the solution iteration if the difference between
286  // this nonlinear iterate and the last is sufficiently small, AND
287  // if the most recent linear system was solved to a sufficient tolerance.
288  if ((norm_delta < nonlinear_tolerance) &&
289  (navier_stokes_system.final_linear_residual() < nonlinear_tolerance))
290  {
291  libMesh::out << " Nonlinear solver converged at step "
292  << l
293  << std::endl;
294  converged = true;
295  break;
296  }
297 
298  // Otherwise, decrease the linear system tolerance. For the inexact Newton
299  // method, the linear solver tolerance needs to decrease as we get closer to
300  // the solution to ensure quadratic convergence. The new linear solver tolerance
301  // is chosen (heuristically) as the square of the previous linear system residual norm.
302  //Real flr2 = final_linear_residual*final_linear_residual;
303  Real new_linear_solver_tolerance = std::min(Utility::pow<2>(final_linear_residual), initial_linear_solver_tol);
304  equation_systems.parameters.set<Real> ("linear solver tolerance") = new_linear_solver_tolerance;
305  } // end nonlinear loop
306 
307  // Don't keep going if we failed to converge.
308  if (!converged)
309  libmesh_error_msg("Error: Newton iterations failed to converge!");
310 
311 #ifdef LIBMESH_HAVE_EXODUS_API
312  // Write out every nth timestep to file.
313  const unsigned int write_interval = 1;
314 
315  if ((t_step+1)%write_interval == 0)
316  {
317  exo_io.write_timestep("out.e",
318  equation_systems,
319  t_step+1, // we're off by one since we wrote the IC and the Exodus numbering is 1-based.
320  navier_stokes_system.time);
321  }
322 #endif // #ifdef LIBMESH_HAVE_EXODUS_API
323  } // end timestep loop.
324 
325  // All done.
326  return 0;
327 }
328 
329 
330 
331 
332 
333 
334 // The matrix assembly function to be called at each time step to
335 // prepare for the linear solve.
337  const std::string & libmesh_dbg_var(system_name))
338 {
339  // It is a good idea to make sure we are assembling
340  // the proper system.
341  libmesh_assert_equal_to (system_name, "Navier-Stokes");
342 
343  // Get a constant reference to the mesh object.
344  const MeshBase & mesh = es.get_mesh();
345 
346  // The dimension that we are running
347  const unsigned int dim = mesh.mesh_dimension();
348 
349  // Get a reference to the Stokes system object.
350  TransientLinearImplicitSystem & navier_stokes_system =
351  es.get_system<TransientLinearImplicitSystem> ("Navier-Stokes");
352 
353  // Numeric ids corresponding to each variable in the system
354  const unsigned int u_var = navier_stokes_system.variable_number ("vel_x");
355  const unsigned int v_var = navier_stokes_system.variable_number ("vel_y");
356  const unsigned int p_var = navier_stokes_system.variable_number ("p");
357 
358  // Get the Finite Element type for "u". Note this will be
359  // the same as the type for "v".
360  FEType fe_vel_type = navier_stokes_system.variable_type(u_var);
361 
362  // Get the Finite Element type for "p".
363  FEType fe_pres_type = navier_stokes_system.variable_type(p_var);
364 
365  // Build a Finite Element object of the specified type for
366  // the velocity variables.
367  UniquePtr<FEBase> fe_vel (FEBase::build(dim, fe_vel_type));
368 
369  // Build a Finite Element object of the specified type for
370  // the pressure variables.
371  UniquePtr<FEBase> fe_pres (FEBase::build(dim, fe_pres_type));
372 
373  // A Gauss quadrature rule for numerical integration.
374  // Let the FEType object decide what order rule is appropriate.
375  QGauss qrule (dim, fe_vel_type.default_quadrature_order());
376 
377  // Tell the finite element objects to use our quadrature rule.
378  fe_vel->attach_quadrature_rule (&qrule);
379  fe_pres->attach_quadrature_rule (&qrule);
380 
381  // Here we define some references to cell-specific data that
382  // will be used to assemble the linear system.
383  //
384  // The element Jacobian * quadrature weight at each integration point.
385  const std::vector<Real> & JxW = fe_vel->get_JxW();
386 
387  // The element shape functions evaluated at the quadrature points.
388  const std::vector<std::vector<Real>> & phi = fe_vel->get_phi();
389 
390  // The element shape function gradients for the velocity
391  // variables evaluated at the quadrature points.
392  const std::vector<std::vector<RealGradient>> & dphi = fe_vel->get_dphi();
393 
394  // The element shape functions for the pressure variable
395  // evaluated at the quadrature points.
396  const std::vector<std::vector<Real>> & psi = fe_pres->get_phi();
397 
398  // The value of the linear shape function gradients at the quadrature points
399  // const std::vector<std::vector<RealGradient>> & dpsi = fe_pres->get_dphi();
400 
401  // A reference to the DofMap object for this system. The DofMap
402  // object handles the index translation from node and element numbers
403  // to degree of freedom numbers. We will talk more about the DofMap
404  // in future examples.
405  const DofMap & dof_map = navier_stokes_system.get_dof_map();
406 
407  // Define data structures to contain the element matrix
408  // and right-hand-side vector contribution. Following
409  // basic finite element terminology we will denote these
410  // "Ke" and "Fe".
413 
415  Kuu(Ke), Kuv(Ke), Kup(Ke),
416  Kvu(Ke), Kvv(Ke), Kvp(Ke),
417  Kpu(Ke), Kpv(Ke), Kpp(Ke);
418 
420  Fu(Fe),
421  Fv(Fe),
422  Fp(Fe);
423 
424  // This vector will hold the degree of freedom indices for
425  // the element. These define where in the global system
426  // the element degrees of freedom get mapped.
427  std::vector<dof_id_type> dof_indices;
428  std::vector<dof_id_type> dof_indices_u;
429  std::vector<dof_id_type> dof_indices_v;
430  std::vector<dof_id_type> dof_indices_p;
431 
432  // Find out what the timestep size parameter is from the system, and
433  // the value of theta for the theta method. We use implicit Euler (theta=1)
434  // for this simulation even though it is only first-order accurate in time.
435  // The reason for this decision is that the second-order Crank-Nicolson
436  // method is notoriously oscillatory for problems with discontinuous
437  // initial data such as the lid-driven cavity. Therefore,
438  // we sacrifice accuracy in time for stability, but since the solution
439  // reaches steady state relatively quickly we can afford to take small
440  // timesteps. If you monitor the initial nonlinear residual for this
441  // simulation, you should see that it is monotonically decreasing in time.
442  const Real dt = es.parameters.get<Real>("dt");
443  const Real theta = 1.;
444 
445  // The kinematic viscosity, multiplies the "viscous" terms.
446  const Real nu = es.parameters.get<Real>("nu");
447 
448  // The system knows whether or not we need to do a pressure pin.
449  // This is only required for problems with all-Dirichlet boundary
450  // conditions on the velocity.
451  const bool pin_pressure = es.parameters.get<bool>("pin_pressure");
452 
453  // Now we will loop over all the elements in the mesh that
454  // live on the local processor. We will compute the element
455  // matrix and right-hand-side contribution. Since the mesh
456  // will be refined we want to only consider the ACTIVE elements,
457  // hence we use a variant of the active_elem_iterator.
458  for (const auto & elem : mesh.active_local_element_ptr_range())
459  {
460  // Get the degree of freedom indices for the
461  // current element. These define where in the global
462  // matrix and right-hand-side this element will
463  // contribute to.
464  dof_map.dof_indices (elem, dof_indices);
465  dof_map.dof_indices (elem, dof_indices_u, u_var);
466  dof_map.dof_indices (elem, dof_indices_v, v_var);
467  dof_map.dof_indices (elem, dof_indices_p, p_var);
468 
469  const unsigned int n_dofs = dof_indices.size();
470  const unsigned int n_u_dofs = dof_indices_u.size();
471  const unsigned int n_v_dofs = dof_indices_v.size();
472  const unsigned int n_p_dofs = dof_indices_p.size();
473 
474  // Compute the element-specific data for the current
475  // element. This involves computing the location of the
476  // quadrature points (q_point) and the shape functions
477  // (phi, dphi) for the current element.
478  fe_vel->reinit (elem);
479  fe_pres->reinit (elem);
480 
481  // Zero the element matrix and right-hand side before
482  // summing them. We use the resize member here because
483  // the number of degrees of freedom might have changed from
484  // the last element. Note that this will be the case if the
485  // element type is different (i.e. the last element was a
486  // triangle, now we are on a quadrilateral).
487  Ke.resize (n_dofs, n_dofs);
488  Fe.resize (n_dofs);
489 
490  // Reposition the submatrices... The idea is this:
491  //
492  // - - - -
493  // | Kuu Kuv Kup | | Fu |
494  // Ke = | Kvu Kvv Kvp |; Fe = | Fv |
495  // | Kpu Kpv Kpp | | Fp |
496  // - - - -
497  //
498  // The DenseSubMatrix.reposition () member takes the
499  // (row_offset, column_offset, row_size, column_size).
500  //
501  // Similarly, the DenseSubVector.reposition () member
502  // takes the (row_offset, row_size)
503  Kuu.reposition (u_var*n_u_dofs, u_var*n_u_dofs, n_u_dofs, n_u_dofs);
504  Kuv.reposition (u_var*n_u_dofs, v_var*n_u_dofs, n_u_dofs, n_v_dofs);
505  Kup.reposition (u_var*n_u_dofs, p_var*n_u_dofs, n_u_dofs, n_p_dofs);
506 
507  Kvu.reposition (v_var*n_v_dofs, u_var*n_v_dofs, n_v_dofs, n_u_dofs);
508  Kvv.reposition (v_var*n_v_dofs, v_var*n_v_dofs, n_v_dofs, n_v_dofs);
509  Kvp.reposition (v_var*n_v_dofs, p_var*n_v_dofs, n_v_dofs, n_p_dofs);
510 
511  Kpu.reposition (p_var*n_u_dofs, u_var*n_u_dofs, n_p_dofs, n_u_dofs);
512  Kpv.reposition (p_var*n_u_dofs, v_var*n_u_dofs, n_p_dofs, n_v_dofs);
513  Kpp.reposition (p_var*n_u_dofs, p_var*n_u_dofs, n_p_dofs, n_p_dofs);
514 
515  Fu.reposition (u_var*n_u_dofs, n_u_dofs);
516  Fv.reposition (v_var*n_u_dofs, n_v_dofs);
517  Fp.reposition (p_var*n_u_dofs, n_p_dofs);
518 
519  // Now we will build the element matrix and right-hand-side.
520  // Constructing the RHS requires the solution and its
521  // gradient from the previous timestep. This must be
522  // calculated at each quadrature point by summing the
523  // solution degree-of-freedom values by the appropriate
524  // weight functions.
525  for (unsigned int qp=0; qp<qrule.n_points(); qp++)
526  {
527  // Values to hold the solution & its gradient at the previous timestep.
528  Number u = 0., u_old = 0.;
529  Number v = 0., v_old = 0.;
530  Number p_old = 0.;
531  Gradient grad_u, grad_u_old;
532  Gradient grad_v, grad_v_old;
533 
534  // Compute the velocity & its gradient from the previous timestep
535  // and the old Newton iterate.
536  for (unsigned int l=0; l<n_u_dofs; l++)
537  {
538  // From the old timestep:
539  u_old += phi[l][qp]*navier_stokes_system.old_solution (dof_indices_u[l]);
540  v_old += phi[l][qp]*navier_stokes_system.old_solution (dof_indices_v[l]);
541  grad_u_old.add_scaled (dphi[l][qp],navier_stokes_system.old_solution (dof_indices_u[l]));
542  grad_v_old.add_scaled (dphi[l][qp],navier_stokes_system.old_solution (dof_indices_v[l]));
543 
544  // From the previous Newton iterate:
545  u += phi[l][qp]*navier_stokes_system.current_solution (dof_indices_u[l]);
546  v += phi[l][qp]*navier_stokes_system.current_solution (dof_indices_v[l]);
547  grad_u.add_scaled (dphi[l][qp],navier_stokes_system.current_solution (dof_indices_u[l]));
548  grad_v.add_scaled (dphi[l][qp],navier_stokes_system.current_solution (dof_indices_v[l]));
549  }
550 
551  // Compute the old pressure value at this quadrature point.
552  for (unsigned int l=0; l<n_p_dofs; l++)
553  p_old += psi[l][qp]*navier_stokes_system.old_solution (dof_indices_p[l]);
554 
555  // Definitions for convenience. It is sometimes simpler to do a
556  // dot product if you have the full vector at your disposal.
557  const NumberVectorValue U_old (u_old, v_old);
558  const NumberVectorValue U (u, v);
559  const Number u_x = grad_u(0);
560  const Number u_y = grad_u(1);
561  const Number v_x = grad_v(0);
562  const Number v_y = grad_v(1);
563 
564  // First, an i-loop over the velocity degrees of freedom.
565  // We know that n_u_dofs == n_v_dofs so we can compute contributions
566  // for both at the same time.
567  for (unsigned int i=0; i<n_u_dofs; i++)
568  {
569  Fu(i) += JxW[qp]*(u_old*phi[i][qp] - // mass-matrix term
570  (1.-theta)*dt*(U_old*grad_u_old)*phi[i][qp] + // convection term
571  (1.-theta)*dt*p_old*dphi[i][qp](0) - // pressure term on rhs
572  (1.-theta)*dt*nu*(grad_u_old*dphi[i][qp]) + // diffusion term on rhs
573  theta*dt*(U*grad_u)*phi[i][qp]); // Newton term
574 
575 
576  Fv(i) += JxW[qp]*(v_old*phi[i][qp] - // mass-matrix term
577  (1.-theta)*dt*(U_old*grad_v_old)*phi[i][qp] + // convection term
578  (1.-theta)*dt*p_old*dphi[i][qp](1) - // pressure term on rhs
579  (1.-theta)*dt*nu*(grad_v_old*dphi[i][qp]) + // diffusion term on rhs
580  theta*dt*(U*grad_v)*phi[i][qp]); // Newton term
581 
582 
583  // Note that the Fp block is identically zero unless we are using
584  // some kind of artificial compressibility scheme...
585 
586  // Matrix contributions for the uu and vv couplings.
587  for (unsigned int j=0; j<n_u_dofs; j++)
588  {
589  Kuu(i,j) += JxW[qp]*(phi[i][qp]*phi[j][qp] + // mass matrix term
590  theta*dt*nu*(dphi[i][qp]*dphi[j][qp]) + // diffusion term
591  theta*dt*(U*dphi[j][qp])*phi[i][qp] + // convection term
592  theta*dt*u_x*phi[i][qp]*phi[j][qp]); // Newton term
593 
594  Kuv(i,j) += JxW[qp]*theta*dt*u_y*phi[i][qp]*phi[j][qp]; // Newton term
595 
596  Kvv(i,j) += JxW[qp]*(phi[i][qp]*phi[j][qp] + // mass matrix term
597  theta*dt*nu*(dphi[i][qp]*dphi[j][qp]) + // diffusion term
598  theta*dt*(U*dphi[j][qp])*phi[i][qp] + // convection term
599  theta*dt*v_y*phi[i][qp]*phi[j][qp]); // Newton term
600 
601  Kvu(i,j) += JxW[qp]*theta*dt*v_x*phi[i][qp]*phi[j][qp]; // Newton term
602  }
603 
604  // Matrix contributions for the up and vp couplings.
605  for (unsigned int j=0; j<n_p_dofs; j++)
606  {
607  Kup(i,j) += JxW[qp]*(-theta*dt*psi[j][qp]*dphi[i][qp](0));
608  Kvp(i,j) += JxW[qp]*(-theta*dt*psi[j][qp]*dphi[i][qp](1));
609  }
610  }
611 
612  // Now an i-loop over the pressure degrees of freedom. This code computes
613  // the matrix entries due to the continuity equation. Note: To maintain a
614  // symmetric matrix, we may (or may not) multiply the continuity equation by
615  // negative one. Here we do not.
616  for (unsigned int i=0; i<n_p_dofs; i++)
617  for (unsigned int j=0; j<n_u_dofs; j++)
618  {
619  Kpu(i,j) += JxW[qp]*psi[i][qp]*dphi[j][qp](0);
620  Kpv(i,j) += JxW[qp]*psi[i][qp]*dphi[j][qp](1);
621  }
622  } // end of the quadrature point qp-loop
623 
624 
625  // At this point the interior element integration has been
626  // completed. We now need to pin the pressure to zero at global
627  // node number "pressure_node". This effectively removes the
628  // non-trivial null space of constant pressure solutions. The
629  // pressure pin is not necessary in problems that have "outflow"
630  // BCs, like Poiseuille flow with natural BCs. In fact it is
631  // actually wrong to do so, since the pressure is not
632  // under-specified in that situation.
633  if (pin_pressure)
634  {
635  const Real penalty = 1.e10;
636  const unsigned int pressure_node = 0;
637  const Real p_value = 0.0;
638  for (auto c : elem->node_index_range())
639  if (elem->node_id(c) == pressure_node)
640  {
641  Kpp(c,c) += penalty;
642  Fp(c) += penalty*p_value;
643  }
644  }
645 
646  // Since we're using heterogeneous DirichletBoundary objects for
647  // the boundary conditions, we need to call a specific function
648  // to constrain the element stiffness matrix.
649  dof_map.heterogenously_constrain_element_matrix_and_vector (Ke, Fe, dof_indices);
650 
651  // The element matrix and right-hand-side are now built
652  // for this element. Add them to the global matrix and
653  // right-hand-side vector. The SparseMatrix::add_matrix()
654  // and NumericVector::add_vector() members do this for us.
655  navier_stokes_system.matrix->add_matrix (Ke, dof_indices);
656  navier_stokes_system.rhs->add_vector (Fe, dof_indices);
657  } // end of element loop
658 }
659 
660 
661 
663 {
664  unsigned short int
665  u_var = system.variable_number("vel_x"),
666  v_var = system.variable_number("vel_y");
667 
668  // This problem *does* require a pressure pin, there are Dirichlet
669  // boundary conditions for u and v on the entire boundary.
670  system.get_equation_systems().parameters.set<bool>("pin_pressure") = true;
671 
672  // Get a convenient reference to the System's DofMap
673  DofMap & dof_map = system.get_dof_map();
674 
675  {
676  // u=v=0 on bottom, left, right
677  std::set<boundary_id_type> boundary_ids;
678  boundary_ids.insert(0);
679  boundary_ids.insert(1);
680  boundary_ids.insert(3);
681 
682  std::vector<unsigned int> variables;
683  variables.push_back(u_var);
684  variables.push_back(v_var);
685 
686  dof_map.add_dirichlet_boundary(DirichletBoundary(boundary_ids,
687  variables,
689  }
690  {
691  // u=1 on top
692  std::set<boundary_id_type> boundary_ids;
693  boundary_ids.insert(2);
694 
695  std::vector<unsigned int> variables;
696  variables.push_back(u_var);
697 
698  dof_map.add_dirichlet_boundary(DirichletBoundary(boundary_ids,
699  variables,
700  ConstFunction<Number>(1.)));
701  }
702  {
703  // v=0 on top
704  std::set<boundary_id_type> boundary_ids;
705  boundary_ids.insert(2);
706 
707  std::vector<unsigned int> variables;
708  variables.push_back(v_var);
709 
710  dof_map.add_dirichlet_boundary(DirichletBoundary(boundary_ids,
711  variables,
713  }
714 }
715 
716 
717 
719 {
720  unsigned short int
721  u_var = system.variable_number("vel_x"),
722  v_var = system.variable_number("vel_y");
723 
724  // This problem does not require a pressure pin, the Neumann outlet
725  // BCs are sufficient to set the value of the pressure.
726  system.get_equation_systems().parameters.set<bool>("pin_pressure") = false;
727 
728  // Get a convenient reference to the System's DofMap
729  DofMap & dof_map = system.get_dof_map();
730 
731  {
732  // u=v=0 on bottom
733  std::set<boundary_id_type> boundary_ids;
734  boundary_ids.insert(0);
735 
736  std::vector<unsigned int> variables;
737  variables.push_back(u_var);
738  variables.push_back(v_var);
739 
740  dof_map.add_dirichlet_boundary(DirichletBoundary(boundary_ids,
741  variables,
743  }
744  {
745  // u=0 on left (symmetry)
746  std::set<boundary_id_type> boundary_ids;
747  boundary_ids.insert(3);
748 
749  std::vector<unsigned int> variables;
750  variables.push_back(u_var);
751 
752  dof_map.add_dirichlet_boundary(DirichletBoundary(boundary_ids,
753  variables,
755  }
756  {
757  // u = k*x on top
758  std::set<boundary_id_type> boundary_ids;
759  boundary_ids.insert(2);
760 
761  std::vector<unsigned int> variables;
762  variables.push_back(u_var);
763 
764  // Set up ParsedFunction parameters
765  std::vector<std::string> additional_vars;
766  additional_vars.push_back("k");
767  std::vector<Number> initial_vals;
768  initial_vals.push_back(1.);
769 
770  dof_map.add_dirichlet_boundary(DirichletBoundary(boundary_ids,
771  variables,
773  &additional_vars,
774  &initial_vals)));
775  }
776  {
777  // v = -k*y on top
778  std::set<boundary_id_type> boundary_ids;
779  boundary_ids.insert(2);
780 
781  std::vector<unsigned int> variables;
782  variables.push_back(v_var);
783 
784  // Set up ParsedFunction parameters
785  std::vector<std::string> additional_vars;
786  additional_vars.push_back("k");
787  std::vector<Number> initial_vals;
788  initial_vals.push_back(1.);
789 
790  // Note: we have to specify LOCAL_VARIABLE_ORDER here, since we're
791  // using a ParsedFunction to set the value of v_var, which is
792  // actually the second variable in the system.
793  dof_map.add_dirichlet_boundary(DirichletBoundary(boundary_ids,
794  variables,
795  ParsedFunction<Number>("-k*y",
796  &additional_vars,
797  &initial_vals),
799  }
800 }
801 
802 
803 
805 {
806  unsigned short int
807  u_var = system.variable_number("vel_x"),
808  v_var = system.variable_number("vel_y");
809 
810  // This problem does not require a pressure pin, the Neumann outlet
811  // BCs are sufficient to set the value of the pressure.
812  system.get_equation_systems().parameters.set<bool>("pin_pressure") = false;
813 
814  // Get a convenient reference to the System's DofMap
815  DofMap & dof_map = system.get_dof_map();
816 
817  {
818  // u=v=0 on top, bottom
819  std::set<boundary_id_type> boundary_ids;
820  boundary_ids.insert(0);
821  boundary_ids.insert(2);
822 
823  std::vector<unsigned int> variables;
824  variables.push_back(u_var);
825  variables.push_back(v_var);
826 
827  dof_map.add_dirichlet_boundary(DirichletBoundary(boundary_ids,
828  variables,
830  }
831  {
832  // u=quadratic on left
833  std::set<boundary_id_type> boundary_ids;
834  boundary_ids.insert(3);
835 
836  std::vector<unsigned int> variables;
837  variables.push_back(u_var);
838 
839  dof_map.add_dirichlet_boundary(DirichletBoundary(boundary_ids,
840  variables,
841  ParsedFunction<Number>("4*y*(1-y)")));
842  }
843  {
844  // v=0 on left
845  std::set<boundary_id_type> boundary_ids;
846  boundary_ids.insert(3);
847 
848  std::vector<unsigned int> variables;
849  variables.push_back(v_var);
850 
851  dof_map.add_dirichlet_boundary(DirichletBoundary(boundary_ids,
852  variables,
854  }
855 }
class FEType hides (possibly multiple) FEFamily and approximation orders, thereby enabling specialize...
Definition: fe_type.h:178
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:162
UniquePtr< NumericVector< Number > > old_local_solution
All the values I need to compute my contribution to the simulation at hand.
ConstFunction that simply returns 0.
Definition: zero_function.h:35
A Function generated (via FParser) by parsing a mathematical expression.
void set_lid_driven_bcs(TransientLinearImplicitSystem &system)
void add_scaled(const TypeVector< T2 > &, const T)
Add a scaled value to this vector without creating a temporary.
Definition: type_vector.h:624
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:350
MeshBase & mesh
This class allows one to associate Dirichlet boundary values with a given set of mesh boundary ids an...
This class provides a specific system class.
const T & get(const std::string &) const
Definition: parameters.h:430
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
The libMesh namespace provides an interface to certain functionality in the library.
This is the MeshBase class.
Definition: mesh_base.h:68
std::unique_ptr< T > UniquePtr
Definition: auto_ptr.h:46
Defines a dense subvector for use in finite element computations.
SolverPackage default_solver_package()
Definition: libmesh.C:995
virtual SimpleRange< element_iterator > active_local_element_ptr_range()=0
The PerfLog class allows monitoring of specific events.
Definition: perf_log.h:125
This class handles the numbering of degrees of freedom on a mesh.
Definition: dof_map.h:167
void set_poiseuille_bcs(TransientLinearImplicitSystem &system)
void reposition(const unsigned int ioff, const unsigned int n)
Changes the location of the subvector in the parent vector.
void assemble_stokes(EquationSystems &es, const std::string &system_name)
void print_info(std::ostream &os=libMesh::out) const
Prints information about the equation systems, by default to libMesh::out.
Order default_quadrature_order() const
Definition: fe_type.h:332
Defines a dense submatrix for use in Finite Element-type computations.
void init(triangulateio &t)
Initializes the fields of t to NULL/0 as necessary.
int main(int argc, char **argv)
virtual System & add_system(const std::string &system_type, const std::string &name)
Add the system of type system_type named name to the systems array.
void push(const char *label, const char *header="")
Push the event label onto the stack, pausing any active event.
Definition: perf_log.C:132
const Parallel::Communicator & comm() const
Definition: libmesh.h:81
PetscErrorCode Vec Mat libmesh_dbg_var(j)
void set_stagnation_bcs(TransientLinearImplicitSystem &system)
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
Constrains the element matrix and vector.
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
T & set(const std::string &)
Definition: parameters.h:469
OStreamProxy out
unsigned int mesh_dimension() const
Definition: mesh_base.C:148
void resize(const unsigned int new_m, const unsigned int new_n)
Resize the matrix.
Definition: dense_matrix.h:776
This class implements specific orders of Gauss quadrature.
Parameters parameters
Data structure holding arbitrary parameters.
Function that returns a single value that never changes.
const MeshBase & get_mesh() const
Number old_solution(const dof_id_type global_dof_number) const
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.
const T_sys & get_system(const std::string &name) const
virtual void init()
Initialize all the systems.
The Mesh class is a thin wrapper, around the ReplicatedMesh class by default.
Definition: mesh.h:50
long double min(long double a, double b)
void print_info(std::ostream &os=libMesh::out) const
Prints relevant information about the mesh.
Definition: mesh_base.C:448
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:1917
static UniquePtr< FEGenericBase > build(const unsigned int dim, const FEType &type)
Builds a specific finite element type.