libMesh
coupled_system.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 #include "libmesh/dirichlet_boundaries.h"
20 #include "libmesh/dof_map.h"
21 #include "libmesh/fe_base.h"
22 #include "libmesh/fe_interface.h"
23 #include "libmesh/fem_context.h"
24 #include "libmesh/getpot.h"
25 #include "libmesh/mesh.h"
26 #include "libmesh/parallel.h"
27 #include "libmesh/quadrature.h"
28 #include "libmesh/string_to_enum.h"
29 #include "libmesh/zero_function.h"
30 
31 #include "coupled_system.h"
32 
33 // Bring in everything from the libMesh namespace
34 using namespace libMesh;
35 
36 // Function to set the Dirichlet boundary function for the inlet boundary velocities
37 class BdyFunction : public FunctionBase<Number>
38 {
39 public:
40  BdyFunction (unsigned int u_var,
41  unsigned int v_var,
42  int sign) :
43  _u_var(u_var),
44  _v_var(v_var),
45  _sign(sign)
46  { this->_initialized = true; }
47 
48  virtual Number operator() (const Point &,
49  const Real = 0)
50  { libmesh_not_implemented(); }
51 
52  virtual void operator() (const Point & p,
53  const Real,
54  DenseVector<Number> & output)
55  {
56  output.resize(2);
57  output.zero();
58  const Real y=p(1);
59  // Set the parabolic inflow boundary conditions at stations 0 & 1
60  output(_u_var) = (_sign)*((y-2) * (y-3));
61  output(_v_var) = 0;
62  }
63 
65  { return UniquePtr<FunctionBase<Number>> (new BdyFunction(_u_var, _v_var, _sign)); }
66 
67 private:
68  const unsigned int _u_var, _v_var;
69  const Real _sign;
70 };
71 
72 
74 {
75  // Check the input file for Reynolds number, application type,
76  // approximation type
77  GetPot infile("coupled_system.in");
78  Peclet = infile("Peclet", 1.);
79  unsigned int pressure_p = infile("pressure_p", 1);
80  std::string fe_family = infile("fe_family", std::string("LAGRANGE"));
81 
82  // LBB needs better-than-quadratic velocities for better-than-linear
83  // pressures, and libMesh needs non-Lagrange elements for
84  // better-than-quadratic velocities.
85  libmesh_assert((pressure_p == 1) || (fe_family != "LAGRANGE"));
86 
87  FEFamily fefamily = Utility::string_to_enum<FEFamily>(fe_family);
88 
89  // Add the velocity components "u" & "v". They
90  // will be approximated using second-order approximation.
91  u_var = this->add_variable ("u", static_cast<Order>(pressure_p+1),
92  fefamily);
93  v_var = this->add_variable ("v", static_cast<Order>(pressure_p+1),
94  fefamily);
95 
96  // Add the pressure variable "p". This will
97  // be approximated with a first-order basis,
98  // providing an LBB-stable pressure-velocity pair.
99  p_var = this->add_variable ("p", static_cast<Order>(pressure_p),
100  fefamily);
101 
102  // Add the Concentration variable "C". They will
103  // be approximated using second-order approximation, the same as the velocity components
104  C_var = this->add_variable ("C", static_cast<Order>(pressure_p+1),
105  fefamily);
106 
107  // Tell the system to march velocity and concentration forward in
108  // time, with first order time derivatives, but leave pressure as a
109  // constraint only
110  this->time_evolving(u_var, 1);
111  this->time_evolving(v_var, 1);
112  this->time_evolving(C_var, 1);
113 
114  // Useful debugging options
115  this->verify_analytic_jacobians = infile("verify_analytic_jacobians", 0.);
116 
117  // Set Dirichlet boundary conditions
118  const boundary_id_type left_inlet_id = 0;
119  std::set<boundary_id_type> left_inlet_bdy;
120  left_inlet_bdy.insert(left_inlet_id);
121 
122  const boundary_id_type right_inlet_id = 1;
123  std::set<boundary_id_type> right_inlet_bdy;
124  right_inlet_bdy.insert(right_inlet_id);
125 
126  const boundary_id_type outlets_id = 2;
127  std::set<boundary_id_type> outlets_bdy;
128  outlets_bdy.insert(outlets_id);
129 
130  const boundary_id_type wall_id = 3;
131  std::set<boundary_id_type> wall_bdy;
132  wall_bdy.insert(wall_id);
133 
134  // The uv identifier for the setting the inlet and wall velocity boundary conditions
135  std::vector<unsigned int> uv(1, u_var);
136  uv.push_back(v_var);
137  // The C_only identifier for setting the concentrations at the inlets
138  std::vector<unsigned int> C_only(1, C_var);
139 
140  // The zero and constant functions
142  ConstFunction<Number> one(1);
143 
144  // We need two boundary functions for the inlets, because the signs on the velocities
145  // will be different
146  int velocity_sign = 1;
147  BdyFunction inflow_left(u_var, v_var, -velocity_sign);
148  BdyFunction inflow_right(u_var, v_var, velocity_sign);
149 
150  // On the walls we will apply the no slip and no penetration boundary condition, u=0, v=0
151  this->get_dof_map().add_dirichlet_boundary
152  (DirichletBoundary (wall_bdy, uv, &zero));
153 
154  // On the inlet (left), we apply parabolic inflow boundary conditions for the velocity, u = - (y-2)*(y-3), v=0
155  // and set C = 1
156  this->get_dof_map().add_dirichlet_boundary
157  (DirichletBoundary (left_inlet_bdy, uv, &inflow_left));
158  this->get_dof_map().add_dirichlet_boundary
159  (DirichletBoundary (left_inlet_bdy, C_only, &one));
160 
161  // On the inlet (right), we apply parabolic inflow boundary conditions for the velocity, u = (y-2)*(y-3), v=0
162  // and set C = 0
163  this->get_dof_map().add_dirichlet_boundary
164  (DirichletBoundary (right_inlet_bdy, uv, &inflow_right));
165  this->get_dof_map().add_dirichlet_boundary
166  (DirichletBoundary (right_inlet_bdy, C_only, &zero));
167 
168  // Note that the remaining boundary conditions are the natural boundary conditions for the concentration
169  // on the wall (grad(c) dot n = 0) and natural boundary conditions for the velocity and the concentration
170  // on the outlets ((grad(velocity) dot n - p n) dot t = 0, grad(C) dot n = 0)
171 
172  // Do the parent's initialization after variables and boundary constraints are defined
174 }
175 
177 {
178  FEMContext & c = cast_ref<FEMContext &>(context);
179 
180  // We should prerequest all the data
181  // we will need to build the linear system.
182  // Note that the concentration and velocity components
183  // use the same basis.
184  FEBase * u_elem_fe = libmesh_nullptr;
185  c.get_element_fe(u_var, u_elem_fe);
186  u_elem_fe->get_JxW();
187  u_elem_fe->get_phi();
188  u_elem_fe->get_dphi();
189  u_elem_fe->get_xyz();
190 
191  FEBase * p_elem_fe = libmesh_nullptr;
192  c.get_element_fe(p_var, p_elem_fe);
193  p_elem_fe->get_phi();
194 
195  FEBase * side_fe = libmesh_nullptr;
196  c.get_side_fe(u_var, side_fe);
197 
198  side_fe->get_JxW();
199  side_fe->get_phi();
200  side_fe->get_xyz();
201 }
202 
203 
204 bool CoupledSystem::element_time_derivative (bool request_jacobian,
205  DiffContext & context)
206 {
207  FEMContext & c = cast_ref<FEMContext &>(context);
208 
209  // First we get some references to cell-specific data that
210  // will be used to assemble the linear system.
211  FEBase * u_elem_fe = libmesh_nullptr;
212  c.get_element_fe(u_var, u_elem_fe);
213 
214  // Element Jacobian * quadrature weights for interior integration
215  const std::vector<Real> & JxW = u_elem_fe->get_JxW();
216 
217  // The velocity shape functions at interior quadrature points.
218  const std::vector<std::vector<Real>> & phi = u_elem_fe->get_phi();
219 
220  // The velocity shape function gradients at interior
221  // quadrature points.
222  const std::vector<std::vector<RealGradient>> & dphi = u_elem_fe->get_dphi();
223 
224  // The pressure shape functions at interior
225  // quadrature points.
226  FEBase * p_elem_fe = libmesh_nullptr;
227  c.get_element_fe(p_var, p_elem_fe);
228 
229  const std::vector<std::vector<Real>> & psi = p_elem_fe->get_phi();
230 
231  // The number of local degrees of freedom in each variable
232  const unsigned int n_p_dofs = c.get_dof_indices(p_var).size();
233  const unsigned int n_u_dofs = c.get_dof_indices(u_var).size();
234  libmesh_assert_equal_to (n_u_dofs, c.get_dof_indices(v_var).size());
235 
236  // The subvectors and submatrices we need to fill:
237  DenseSubMatrix<Number> & Kuu = c.get_elem_jacobian(u_var, u_var);
238  DenseSubMatrix<Number> & Kup = c.get_elem_jacobian(u_var, p_var);
240 
241  DenseSubMatrix<Number> & Kvv = c.get_elem_jacobian(v_var, v_var);
242  DenseSubMatrix<Number> & Kvp = c.get_elem_jacobian(v_var, p_var);
244 
245  DenseSubMatrix<Number> & KCu = c.get_elem_jacobian(C_var, u_var);
246  DenseSubMatrix<Number> & KCv = c.get_elem_jacobian(C_var, v_var);
247  DenseSubMatrix<Number> & KCC = c.get_elem_jacobian(C_var, C_var);
249 
250  // Now we will build the element Jacobian and residual.
251  // Constructing the residual requires the solution and its
252  // gradient from the previous timestep. This must be
253  // calculated at each quadrature point by summing the
254  // solution degree-of-freedom values by the appropriate
255  // weight functions.
256  unsigned int n_qpoints = c.get_element_qrule().n_points();
257 
258  for (unsigned int qp=0; qp != n_qpoints; qp++)
259  {
260  // Compute the solution & its gradient at the old Newton iterate
261  Number p = c.interior_value(p_var, qp),
262  u = c.interior_value(u_var, qp),
263  v = c.interior_value(v_var, qp);
264  Gradient grad_u = c.interior_gradient(u_var, qp),
265  grad_v = c.interior_gradient(v_var, qp),
266  grad_C = c.interior_gradient(C_var, qp);
267 
268  // Definitions for convenience. It is sometimes simpler to do a
269  // dot product if you have the full vector at your disposal.
270  NumberVectorValue U (u, v);
271  const Number C_x = grad_C(0);
272  const Number C_y = grad_C(1);
273 
274  // First, an i-loop over the velocity degrees of freedom.
275  // We know that n_u_dofs == n_v_dofs so we can compute contributions
276  // for both at the same time.
277  for (unsigned int i=0; i != n_u_dofs; i++)
278  {
279  // Stokes equations residuals
280  Fu(i) += JxW[qp] *
281  (p*dphi[i][qp](0) - // pressure term
282  (grad_u*dphi[i][qp])); // diffusion term
283 
284  Fv(i) += JxW[qp] *
285  (p*dphi[i][qp](1) - // pressure term
286  (grad_v*dphi[i][qp])); // diffusion term
287 
288  // Concentration Equation Residual
289  FC(i) += JxW[qp] *
290  ((U*grad_C)*phi[i][qp] + // convection term
291  (1./Peclet)*(grad_C*dphi[i][qp])); // diffusion term
292 
293  // Note that the Fp block is identically zero unless we are using
294  // some kind of artificial compressibility scheme...
295 
296  if (request_jacobian && c.elem_solution_derivative)
297  {
299 
300  // Matrix contributions for the uu and vv couplings.
301  for (unsigned int j=0; j != n_u_dofs; j++)
302  {
303  Kuu(i,j) += JxW[qp] * (-(dphi[i][qp]*dphi[j][qp])); // diffusion term
304 
305  Kvv(i,j) += JxW[qp] * (-(dphi[i][qp]*dphi[j][qp])); // diffusion term
306 
307  KCu(i,j) += JxW[qp] * ((phi[j][qp]*C_x)*phi[i][qp]); // convection term
308 
309  KCv(i,j) += JxW[qp] * ((phi[j][qp]*C_y)*phi[i][qp]); // convection term
310 
311  KCC(i,j) += JxW[qp]*
312  ((U*dphi[j][qp])*phi[i][qp] + // nonlinear term (convection)
313  (1./Peclet)*(dphi[j][qp]*dphi[i][qp])); // diffusion term
314  }
315 
316  // Matrix contributions for the up and vp couplings.
317  for (unsigned int j=0; j != n_p_dofs; j++)
318  {
319  Kup(i,j) += JxW[qp]*psi[j][qp]*dphi[i][qp](0);
320  Kvp(i,j) += JxW[qp]*psi[j][qp]*dphi[i][qp](1);
321  }
322  }
323  }
324  } // end of the quadrature point qp-loop
325 
326  return request_jacobian;
327 }
328 
329 
330 bool CoupledSystem::element_constraint (bool request_jacobian,
331  DiffContext & context)
332 {
333  FEMContext & c = cast_ref<FEMContext &>(context);
334 
335  // Here we define some references to cell-specific data that
336  // will be used to assemble the linear system.
337  FEBase * u_elem_fe = libmesh_nullptr;
338  c.get_element_fe(u_var, u_elem_fe);
339 
340  FEBase * p_elem_fe = libmesh_nullptr;
341  c.get_element_fe(p_var, p_elem_fe);
342 
343  // Element Jacobian * quadrature weight for interior integration
344  const std::vector<Real> & JxW = u_elem_fe->get_JxW();
345 
346  // The velocity shape function gradients at interior
347  // quadrature points.
348  const std::vector<std::vector<RealGradient>> & dphi = u_elem_fe->get_dphi();
349 
350  // The pressure shape functions at interior
351  // quadrature points.
352  const std::vector<std::vector<Real>> & psi = p_elem_fe->get_phi();
353 
354  // The number of local degrees of freedom in each variable
355  const unsigned int n_u_dofs = c.get_dof_indices(u_var).size();
356  const unsigned int n_p_dofs = c.get_dof_indices(p_var).size();
357 
358  // The subvectors and submatrices we need to fill:
359  DenseSubMatrix<Number> & Kpu = c.get_elem_jacobian(p_var, u_var);
360  DenseSubMatrix<Number> & Kpv = c.get_elem_jacobian(p_var, v_var);
362 
363  // Add the constraint given by the continuity equation
364  unsigned int n_qpoints = c.get_element_qrule().n_points();
365  for (unsigned int qp=0; qp != n_qpoints; qp++)
366  {
367  // Compute the velocity gradient at the old Newton iterate
368  Gradient
369  grad_u = c.interior_gradient(u_var, qp),
370  grad_v = c.interior_gradient(v_var, qp);
371 
372  // Now a loop over the pressure degrees of freedom. This
373  // computes the contributions of the continuity equation.
374  for (unsigned int i=0; i != n_p_dofs; i++)
375  {
376  Fp(i) += JxW[qp] * psi[i][qp] *
377  (grad_u(0) + grad_v(1));
378 
379  if (request_jacobian && c.elem_solution_derivative)
380  {
382 
383  for (unsigned int j=0; j != n_u_dofs; j++)
384  {
385  Kpu(i,j) += JxW[qp]*psi[i][qp]*dphi[j][qp](0);
386  Kpv(i,j) += JxW[qp]*psi[i][qp]*dphi[j][qp](1);
387  }
388  }
389  }
390  } // end of the quadrature point qp-loop
391 
392  return request_jacobian;
393 }
394 
396 {
397  // We need to overload the postprocess function to set the
398  // computed_QoI variable of the CoupledSystem class to the qoi value
399  // stored in System::qoi[0]
400  computed_QoI = System::qoi[0];
401 }
402 
403 // These functions supply the nonlinear weighting for the adjoint residual error estimate which
404 // arise due to the convection term in the convection-diffusion equation:
405 // ||e((u_1)_h C,1)||_{L2} ||e(C^*)||_{L2} + ||e(u1 C,1_h)||_{L2} ||e(C^*)||_{L2}
406 // ||e((u_2)_h C,2)||_{L2} ||e(C^*)||_{L2} + ||e(u2 C,2_h)||_{L2} ||e(C^*)||_{L2}
407 // These functions compute (u_1)_h or C,1_h , and (u_2)_h or C,2_h , and supply it to the weighted patch recovery error estimator
408 // In CoupledFEMFunctionsx, the object built with var = 0, returns the (u_1)_h weight, while
409 // the object built with var = 1, returns the C,1_h weight. The switch statement
410 // distinguishes the behavior of the two objects
411 // Same thing for CoupledFEMFunctionsy
413  const Point & p,
414  const Real /* time */)
415 {
416  Number weight = 0.0;
417 
418  switch(var)
419  {
420  case 0:
421  {
422  Gradient grad_C = c.point_gradient(3, p);
423 
424  weight = grad_C(0);
425  }
426  break;
427 
428  case 3:
429  {
430  Number u = c.point_value(0, p);
431 
432  weight = u;
433  }
434  break;
435 
436  default:
437  libmesh_error_msg("Wrong variable number " \
438  << var \
439  << " passed to CoupledFEMFunctionsx object! Quitting!");
440  }
441 
442  return weight;
443 }
444 
446  const Point & p,
447  const Real /* time */)
448 {
449  Number weight = 0.0;
450 
451  switch(var)
452  {
453  case 1:
454  {
455  Gradient grad_C = c.point_gradient(3, p);
456 
457  weight = grad_C(1);
458  }
459  break;
460 
461  case 3:
462  {
463  Number v = c.point_value(1, p);
464 
465  weight = v;
466  }
467  break;
468 
469  default:
470  libmesh_error_msg("Wrong variable number " \
471  << var \
472  << " passed to CoupledFEMFunctionsy object! Quitting!");
473  }
474 
475  return weight;
476 }
virtual void init_data()
Initializes the member data fields associated with the system, so that, e.g., assemble() may be used...
This class provides all data required for a physics package (e.g.
Definition: diff_context.h:54
ConstFunction that simply returns 0.
Definition: zero_function.h:35
Number point_value(unsigned int var, const Point &p) const
Definition: fem_context.C:788
void resize(const unsigned int n)
Resize the vector.
Definition: dense_vector.h:350
virtual bool element_constraint(bool request_jacobian, DiffContext &context)
Adds the constraint contribution on elem to elem_residual.
virtual void zero() libmesh_override
Set every element in the vector to 0.
Definition: dense_vector.h:374
const std::vector< dof_id_type > & get_dof_indices() const
Accessor for element dof indices.
Definition: diff_context.h:366
const class libmesh_nullptr_t libmesh_nullptr
This class allows one to associate Dirichlet boundary values with a given set of mesh boundary ids an...
virtual void postprocess()
Runs a postprocessing loop over all elements, and if postprocess_sides is true over all sides...
const unsigned int _v_var
void get_element_fe(unsigned int var, FEGenericBase< OutputShape > *&fe) const
Accessor for interior finite element object for variable var for the largest dimension in the mesh...
Definition: fem_context.h:262
The libMesh namespace provides an interface to certain functionality in the library.
const Number zero
.
Definition: libmesh.h:178
const DenseVector< Number > & get_elem_residual() const
Const accessor for element residual.
Definition: diff_context.h:248
Number interior_value(unsigned int var, unsigned int qp) const
Definition: fem_context.C:326
libmesh_assert(j)
std::unique_ptr< T > UniquePtr
Definition: auto_ptr.h:46
Defines a dense subvector for use in finite element computations.
BdyFunction(unsigned int u_var, unsigned int v_var, int sign)
const Real _sign
const DenseMatrix< Number > & get_elem_jacobian() const
Const accessor for element Jacobian.
Definition: diff_context.h:282
FEFamily
defines an enum for finite element families.
int8_t boundary_id_type
Definition: id_types.h:51
const std::vector< Real > & get_JxW() const
Definition: fe_abstract.h:237
std::vector< Number > qoi
Values of the quantities of interest.
Definition: system.h:1553
dof_id_type weight(const MeshBase &mesh, const processor_id_type pid)
Definition: mesh_tools.C:245
virtual void init_context(DiffContext &context)
Defines a dense submatrix for use in Finite Element-type computations.
const std::vector< std::vector< OutputShape > > & get_phi() const
Definition: fe_base.h:208
Real elem_solution_derivative
The derivative of elem_solution with respect to the current nonlinear solution.
Definition: diff_context.h:483
This class provides all data required for a physics package (e.g.
Definition: fem_context.h:61
virtual void operator()(const FEMContext &, const Point &, const Real, DenseVector< Number > &)
Gradient interior_gradient(unsigned int var, unsigned int qp) const
Definition: fem_context.C:381
virtual UniquePtr< FunctionBase< Number > > clone() const
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
virtual bool element_time_derivative(bool request_jacobian, DiffContext &context)
Adds the time derivative contribution on elem to elem_residual.
const QBase & get_element_qrule() const
Accessor for element interior quadrature rule for the dimension of the current _elem.
Definition: fem_context.h:765
Gradient point_gradient(unsigned int var, const Point &p) const
Definition: fem_context.C:834
Function that returns a single value that never changes.
virtual void init_data() libmesh_override
Initializes the member data fields associated with the system, so that, e.g., assemble() may be used...
Definition: fem_system.C:845
This is the base class for functor-like classes.
void get_side_fe(unsigned int var, FEGenericBase< OutputShape > *&fe) const
Accessor for edge/face (2D/3D) finite element object for variable var for the largest dimension in th...
Definition: fem_context.h:299
A Point defines a location in LIBMESH_DIM dimensional Real space.
Definition: point.h:38
unsigned int n_points() const
Definition: quadrature.h:113
sys get_dof_map()
const std::vector< Point > & get_xyz() const
Definition: fe_abstract.h:230
This class forms the foundation from which generic finite elements may be derived.
virtual void operator()(const FEMContext &, const Point &, const Real, DenseVector< Number > &)