GRINS-0.8.0
inc_navier_stokes_adjoint_stab.C
Go to the documentation of this file.
1 //-----------------------------------------------------------------------bl-
2 //--------------------------------------------------------------------------
3 //
4 // GRINS - General Reacting Incompressible Navier-Stokes
5 //
6 // Copyright (C) 2014-2017 Paul T. Bauman, Roy H. Stogner
7 // Copyright (C) 2010-2013 The PECOS Development Team
8 //
9 // This library is free software; you can redistribute it and/or
10 // modify it under the terms of the Version 2.1 GNU Lesser General
11 // Public License as published by the Free Software Foundation.
12 //
13 // This library is distributed in the hope that it will be useful,
14 // but WITHOUT ANY WARRANTY; without even the implied warranty of
15 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
16 // Lesser General Public License for more details.
17 //
18 // You should have received a copy of the GNU Lesser General Public
19 // License along with this library; if not, write to the Free Software
20 // Foundation, Inc. 51 Franklin Street, Fifth Floor,
21 // Boston, MA 02110-1301 USA
22 //
23 //-----------------------------------------------------------------------el-
24 
25 
26 // This class
27 #include "grins/inc_navier_stokes_adjoint_stab.h"
28 
29 // GRINS
30 #include "grins/assembly_context.h"
31 #include "grins/inc_nav_stokes_macro.h"
32 
33 //libMesh
34 #include "libmesh/quadrature.h"
35 
36 namespace GRINS
37 {
38 
39  template<class Mu>
40  IncompressibleNavierStokesAdjointStabilization<Mu>::IncompressibleNavierStokesAdjointStabilization( const std::string& physics_name,
41  const GetPot& input )
42  : IncompressibleNavierStokesStabilizationBase<Mu>(physics_name,input)
43  {}
44 
45  template<class Mu>
46  void IncompressibleNavierStokesAdjointStabilization<Mu>::element_time_derivative
47  ( bool compute_jacobian,
48  AssemblyContext & context )
49  {
50  // The number of local degrees of freedom in each variable.
51  const unsigned int n_p_dofs = context.get_dof_indices(this->_press_var.p()).size();
52  const unsigned int n_u_dofs = context.get_dof_indices(this->_flow_vars.u()).size();
53 
54  // Element Jacobian * quadrature weights for interior integration.
55  const std::vector<libMesh::Real> &JxW =
56  context.get_element_fe(this->_flow_vars.u())->get_JxW();
57 
58  const std::vector<std::vector<libMesh::RealGradient> >& p_gradphi =
59  context.get_element_fe(this->_press_var.p())->get_dphi();
60 
61  const std::vector<std::vector<libMesh::Real> >& u_phi =
62  context.get_element_fe(this->_flow_vars.u())->get_phi();
63 
64  const std::vector<std::vector<libMesh::RealGradient> >& u_gradphi =
65  context.get_element_fe(this->_flow_vars.u())->get_dphi();
66 
67  const std::vector<std::vector<libMesh::RealTensor> >& u_hessphi =
68  context.get_element_fe(this->_flow_vars.u())->get_d2phi();
69 
70  libMesh::DenseSubVector<libMesh::Number> &Fu = context.get_elem_residual(this->_flow_vars.u()); // R_{p}
71  libMesh::DenseSubVector<libMesh::Number> &Fv = context.get_elem_residual(this->_flow_vars.v()); // R_{p}
72  libMesh::DenseSubMatrix<libMesh::Number> &Kup =
73  context.get_elem_jacobian(this->_flow_vars.u(), this->_press_var.p()); // J_{up}
74  libMesh::DenseSubMatrix<libMesh::Number> &Kuu =
75  context.get_elem_jacobian(this->_flow_vars.u(), this->_flow_vars.u()); // J_{uu}
76  libMesh::DenseSubMatrix<libMesh::Number> &Kuv =
77  context.get_elem_jacobian(this->_flow_vars.u(), this->_flow_vars.v()); // J_{uv}
78  libMesh::DenseSubMatrix<libMesh::Number> &Kvp =
79  context.get_elem_jacobian(this->_flow_vars.v(), this->_press_var.p()); // J_{vp}
80  libMesh::DenseSubMatrix<libMesh::Number> &Kvu =
81  context.get_elem_jacobian(this->_flow_vars.v(), this->_flow_vars.u()); // J_{vu}
82  libMesh::DenseSubMatrix<libMesh::Number> &Kvv =
83  context.get_elem_jacobian(this->_flow_vars.v(), this->_flow_vars.v()); // J_{vv}
84 
85  libMesh::DenseSubVector<libMesh::Number> *Fw = NULL;
86  libMesh::DenseSubMatrix<libMesh::Number> *Kuw = NULL;
87  libMesh::DenseSubMatrix<libMesh::Number> *Kvw = NULL;
88  libMesh::DenseSubMatrix<libMesh::Number> *Kwp = NULL;
89  libMesh::DenseSubMatrix<libMesh::Number> *Kwu = NULL;
90  libMesh::DenseSubMatrix<libMesh::Number> *Kwv = NULL;
91  libMesh::DenseSubMatrix<libMesh::Number> *Kww = NULL;
92 
93 
94  if(this->_flow_vars.dim() == 3)
95  {
96  Fw = &context.get_elem_residual(this->_flow_vars.w()); // R_{w}
97  Kuw = &context.get_elem_jacobian
98  (this->_flow_vars.u(), this->_flow_vars.w()); // J_{uw}
99  Kvw = &context.get_elem_jacobian
100  (this->_flow_vars.v(), this->_flow_vars.w()); // J_{vw}
101  Kwp = &context.get_elem_jacobian
102  (this->_flow_vars.w(), this->_press_var.p()); // J_{wp}
103  Kwu = &context.get_elem_jacobian
104  (this->_flow_vars.w(), this->_flow_vars.u()); // J_{wu}
105  Kwv = &context.get_elem_jacobian
106  (this->_flow_vars.w(), this->_flow_vars.v()); // J_{wv}
107  Kww = &context.get_elem_jacobian
108  (this->_flow_vars.w(), this->_flow_vars.w()); // J_{ww}
109  }
110 
111  unsigned int n_qpoints = context.get_element_qrule().n_points();
112 
113  libMesh::FEBase* fe = context.get_element_fe(this->_flow_vars.u());
114 
115  for (unsigned int qp=0; qp != n_qpoints; qp++)
116  {
117  libMesh::RealGradient g = this->_stab_helper.compute_g( fe, context, qp );
118  libMesh::RealTensor G = this->_stab_helper.compute_G( fe, context, qp );
119 
120  libMesh::RealGradient U( context.interior_value( this->_flow_vars.u(), qp ),
121  context.interior_value( this->_flow_vars.v(), qp ) );
122  if( this->_flow_vars.dim() == 3 )
123  {
124  U(2) = context.interior_value( this->_flow_vars.w(), qp );
125  }
126 
127  /*
128  libMesh::Gradient grad_u, grad_v, grad_w;
129  grad_u = context.interior_gradient(this->_flow_vars.u(), qp);
130  grad_v = context.interior_gradient(this->_flow_vars.v(), qp);
131  if (_dim == 3)
132  grad_w = context.interior_gradient(this->_flow_vars.w(), qp);
133  */
134 
135  libMesh::Real tau_M, tau_C;
136  libMesh::Real d_tau_M_d_rho, d_tau_C_d_rho;
137  libMesh::Gradient d_tau_M_dU, d_tau_C_dU;
138  libMesh::Gradient RM_s, d_RM_s_uvw_dgraduvw;
139  libMesh::Real RC;
140  libMesh::Tensor d_RC_dgradU,
141  d_RM_s_dgradp, d_RM_s_dU, d_RM_s_uvw_dhessuvw;
142 
143  // Compute the viscosity at this qp
144  libMesh::Real _mu_qp = this->_mu(context, qp);
145 
146  if (compute_jacobian)
147  {
148  this->_stab_helper.compute_tau_momentum_and_derivs
149  ( context, qp, g, G, this->_rho, U, _mu_qp,
150  tau_M, d_tau_M_d_rho, d_tau_M_dU,
151  this->_is_steady );
152  this->_stab_helper.compute_tau_continuity_and_derivs
153  ( tau_M, d_tau_M_d_rho, d_tau_M_dU,
154  g,
155  tau_C, d_tau_C_d_rho, d_tau_C_dU );
156  this->_stab_helper.compute_res_momentum_steady_and_derivs
157  ( context, qp, this->_rho, _mu_qp,
158  RM_s, d_RM_s_dgradp, d_RM_s_dU, d_RM_s_uvw_dgraduvw,
159  d_RM_s_uvw_dhessuvw);
160  this->_stab_helper.compute_res_continuity_and_derivs
161  ( context, qp, RC, d_RC_dgradU );
162  }
163  else
164  {
165  tau_M = this->_stab_helper.compute_tau_momentum( context, qp, g, G, this->_rho, U, _mu_qp, this->_is_steady );
166  tau_C = this->_stab_helper.compute_tau_continuity( tau_M, g );
167  RM_s = this->_stab_helper.compute_res_momentum_steady( context, qp, this->_rho, _mu_qp );
168  RC = this->_stab_helper.compute_res_continuity( context, qp );
169  }
170 
171  for (unsigned int i=0; i != n_u_dofs; i++)
172  {
173  libMesh::Real test_func = this->_rho*U*u_gradphi[i][qp] +
174  _mu_qp*( u_hessphi[i][qp](0,0) + u_hessphi[i][qp](1,1) + u_hessphi[i][qp](2,2) );
175  libMesh::Gradient d_test_func_dU = this->_rho*u_gradphi[i][qp];
176 
177  //libMesh::RealGradient zeroth_order_term = - this->_rho*u_phi[i][qp]*(grad_u + grad_v + grad_w);
178 
179  Fu(i) += ( -tau_M*RM_s(0)*test_func - tau_C*RC*u_gradphi[i][qp](0) )*JxW[qp];
180 
181  Fv(i) += ( -tau_M*RM_s(1)*test_func - tau_C*RC*u_gradphi[i][qp](1) )*JxW[qp];
182 
183  if(this->_flow_vars.dim() == 3)
184  {
185  (*Fw)(i) += ( -tau_M*RM_s(2)*test_func - tau_C*RC*u_gradphi[i][qp](2) )*JxW[qp];
186  }
187 
188  if (compute_jacobian)
189  {
190  const libMesh::Real fixed_deriv =
191  context.get_fixed_solution_derivative();
192  for (unsigned int j=0; j != n_p_dofs; j++)
193  {
194  Kup(i,j) += ( -tau_M*
195  ( d_RM_s_dgradp(0,0) *
196  p_gradphi[j][qp](0) +
197  d_RM_s_dgradp(0,1) *
198  p_gradphi[j][qp](1) +
199  d_RM_s_dgradp(0,2) *
200  p_gradphi[j][qp](2)
201  )*test_func)*fixed_deriv*JxW[qp];
202  Kvp(i,j) += ( -tau_M*
203  ( d_RM_s_dgradp(1,0) *
204  p_gradphi[j][qp](0) +
205  d_RM_s_dgradp(1,1) *
206  p_gradphi[j][qp](1) +
207  d_RM_s_dgradp(1,2) *
208  p_gradphi[j][qp](2)
209  )*test_func)*fixed_deriv*JxW[qp];
210  }
211  for (unsigned int j=0; j != n_u_dofs; j++)
212  {
213  Kuu(i,j) += ( -d_tau_M_dU(0)*RM_s(0)*test_func
214  -tau_M*d_RM_s_dU(0,0)*test_func
215  - d_tau_C_dU(0)*RC*u_gradphi[i][qp](0)
216  )*fixed_deriv*u_phi[j][qp]*JxW[qp];
217  Kuu(i,j) += ( -tau_M*RM_s(0)*d_test_func_dU(0)
218  )*u_phi[j][qp]*JxW[qp];
219  Kuu(i,j) += ( - tau_C*
220  (
221  d_RC_dgradU(0,0)*u_gradphi[j][qp](0) +
222  d_RC_dgradU(0,1)*u_gradphi[j][qp](1) +
223  d_RC_dgradU(0,2)*u_gradphi[j][qp](2)
224  )
225  *fixed_deriv*u_gradphi[i][qp](0)
226  )*JxW[qp];
227  Kuu(i,j) += ( -tau_M*test_func*(d_RM_s_uvw_dgraduvw*u_gradphi[j][qp])
228  )*fixed_deriv*JxW[qp];
229  Kuu(i,j) += ( -tau_M*test_func*(d_RM_s_uvw_dhessuvw.contract(u_hessphi[j][qp]))
230  )*fixed_deriv*JxW[qp];
231  Kuv(i,j) += ( -d_tau_M_dU(1)*RM_s(0)*test_func
232  -tau_M*d_RM_s_dU(0,1)*test_func
233  )*fixed_deriv*u_phi[j][qp]*JxW[qp];
234  Kuv(i,j) += ( -tau_M*RM_s(0)*d_test_func_dU(1)
235  )*u_phi[j][qp]*JxW[qp];
236  Kuv(i,j) += ( - tau_C*
237  (
238  d_RC_dgradU(1,0)*u_gradphi[j][qp](0) +
239  d_RC_dgradU(1,1)*u_gradphi[j][qp](1) +
240  d_RC_dgradU(1,2)*u_gradphi[j][qp](2)
241  )
242  *fixed_deriv*u_gradphi[i][qp](0)
243  )*JxW[qp];
244  Kvu(i,j) += ( -d_tau_M_dU(0)*RM_s(1)*test_func
245  -tau_M*d_RM_s_dU(1,0)*test_func
246  )*fixed_deriv*u_phi[j][qp]*JxW[qp];
247  Kvu(i,j) += ( -tau_M*RM_s(1)*d_test_func_dU(0)
248  )*u_phi[j][qp]*JxW[qp];
249  Kvu(i,j) += ( - tau_C*
250  (
251  d_RC_dgradU(0,0)*u_gradphi[j][qp](0) +
252  d_RC_dgradU(0,1)*u_gradphi[j][qp](1) +
253  d_RC_dgradU(0,2)*u_gradphi[j][qp](2)
254  )
255  *fixed_deriv*u_gradphi[i][qp](1)
256  )*JxW[qp];
257  Kvv(i,j) += ( -d_tau_M_dU(1)*RM_s(1)*test_func
258  -tau_M*d_RM_s_dU(1,1)*test_func
259  )*fixed_deriv*u_phi[j][qp]*JxW[qp];
260  Kvv(i,j) += ( -tau_M*RM_s(1)*d_test_func_dU(1)
261  )*u_phi[j][qp]*JxW[qp];
262  Kvv(i,j) += ( - tau_C*
263  (
264  d_RC_dgradU(1,0)*u_gradphi[j][qp](0) +
265  d_RC_dgradU(1,1)*u_gradphi[j][qp](1) +
266  d_RC_dgradU(1,2)*u_gradphi[j][qp](2)
267  )
268  *fixed_deriv*u_gradphi[i][qp](1)
269  )*JxW[qp];
270  Kvv(i,j) += ( -tau_M*test_func*(d_RM_s_uvw_dgraduvw*u_gradphi[j][qp])
271  )*fixed_deriv*JxW[qp];
272  Kvv(i,j) += ( -tau_M*test_func*(d_RM_s_uvw_dhessuvw.contract(u_hessphi[j][qp]))
273  )*fixed_deriv*JxW[qp];
274  }
275  if(this->_flow_vars.dim() == 3)
276  {
277  for (unsigned int j=0; j != n_p_dofs; j++)
278  {
279  (*Kwp)(i,j) += ( -tau_M*
280  ( d_RM_s_dgradp(2,0) *
281  p_gradphi[j][qp](0) +
282  d_RM_s_dgradp(2,1) *
283  p_gradphi[j][qp](1) +
284  d_RM_s_dgradp(2,2) *
285  p_gradphi[j][qp](2)
286  )*test_func)*fixed_deriv*JxW[qp];
287  }
288  for (unsigned int j=0; j != n_u_dofs; j++)
289  {
290  (*Kuw)(i,j) += ( -d_tau_M_dU(2)*RM_s(0)*test_func
291  -tau_M*d_RM_s_dU(0,2)*test_func
292  )*fixed_deriv*u_phi[j][qp]*JxW[qp];
293  (*Kuw)(i,j) += ( -tau_M*RM_s(0)*d_test_func_dU(2)
294  )*u_phi[j][qp]*JxW[qp];
295  (*Kuw)(i,j) += ( - tau_C*
296  (
297  d_RC_dgradU(2,0)*u_gradphi[j][qp](0) +
298  d_RC_dgradU(2,1)*u_gradphi[j][qp](1) +
299  d_RC_dgradU(2,2)*u_gradphi[j][qp](2)
300  )
301  *fixed_deriv* u_gradphi[i][qp](0)
302  )*JxW[qp];
303  (*Kvw)(i,j) += ( -d_tau_M_dU(2)*RM_s(1)*test_func
304  -tau_M*d_RM_s_dU(1,2)*test_func
305  )*fixed_deriv*u_phi[j][qp]*JxW[qp];
306  (*Kvw)(i,j) += ( -tau_M*RM_s(1)*d_test_func_dU(2)
307  )*u_phi[j][qp]*JxW[qp];
308  (*Kvw)(i,j) += ( - tau_C*
309  (
310  d_RC_dgradU(2,0)*u_gradphi[j][qp](0) +
311  d_RC_dgradU(2,1)*u_gradphi[j][qp](1) +
312  d_RC_dgradU(2,2)*u_gradphi[j][qp](2)
313  )
314  *fixed_deriv* u_gradphi[i][qp](1)
315  )*JxW[qp];
316  (*Kwu)(i,j) += ( -d_tau_M_dU(0)*RM_s(2)*test_func
317  -tau_M*d_RM_s_dU(2,0)*test_func
318  )*fixed_deriv*u_phi[j][qp]*JxW[qp];
319  (*Kwu)(i,j) += ( -tau_M*RM_s(2)*d_test_func_dU(0)
320  )*u_phi[j][qp]*JxW[qp];
321  (*Kwu)(i,j) += ( - tau_C*
322  (
323  d_RC_dgradU(0,0)*u_gradphi[j][qp](0) +
324  d_RC_dgradU(0,1)*u_gradphi[j][qp](1) +
325  d_RC_dgradU(0,2)*u_gradphi[j][qp](2)
326  )
327  *fixed_deriv* u_gradphi[i][qp](2)
328  )*JxW[qp];
329  (*Kwv)(i,j) += ( -d_tau_M_dU(1)*RM_s(2)*test_func
330  -tau_M*d_RM_s_dU(2,1)*test_func
331  )*fixed_deriv*u_phi[j][qp]*JxW[qp];
332  (*Kwv)(i,j) += ( -tau_M*RM_s(2)*d_test_func_dU(1)
333  )*u_phi[j][qp]*JxW[qp];
334  (*Kwv)(i,j) += ( - tau_C*
335  (
336  d_RC_dgradU(1,0)*u_gradphi[j][qp](0) +
337  d_RC_dgradU(1,1)*u_gradphi[j][qp](1) +
338  d_RC_dgradU(1,2)*u_gradphi[j][qp](2)
339  )
340  *fixed_deriv* u_gradphi[i][qp](2)
341  )*JxW[qp];
342  (*Kww)(i,j) += ( -d_tau_M_dU(2)*RM_s(2)*test_func
343  -tau_M*d_RM_s_dU(2,2)*test_func
344  )*fixed_deriv*u_phi[j][qp]*JxW[qp];
345  (*Kww)(i,j) += ( -tau_M*RM_s(2)*d_test_func_dU(2)
346  )*u_phi[j][qp]*JxW[qp];
347  (*Kww)(i,j) += ( - tau_C*
348  (
349  d_RC_dgradU(2,0)*u_gradphi[j][qp](0) +
350  d_RC_dgradU(2,1)*u_gradphi[j][qp](1) +
351  d_RC_dgradU(2,2)*u_gradphi[j][qp](2)
352  )
353  *fixed_deriv* u_gradphi[i][qp](2)
354  )*JxW[qp];
355  (*Kww)(i,j) += ( -tau_M*test_func*(d_RM_s_uvw_dgraduvw*u_gradphi[j][qp])
356  )*fixed_deriv*JxW[qp];
357  (*Kww)(i,j) += ( -tau_M*test_func*(d_RM_s_uvw_dhessuvw.contract(u_hessphi[j][qp]))
358  )*fixed_deriv*JxW[qp];
359 
360  }
361  }
362  }
363  }
364  }
365  }
366 
367  template<class Mu>
368  void IncompressibleNavierStokesAdjointStabilization<Mu>::element_constraint
369  ( bool compute_jacobian, AssemblyContext & context )
370  {
371  // The number of local degrees of freedom in each variable.
372  const unsigned int n_p_dofs = context.get_dof_indices(this->_press_var.p()).size();
373  const unsigned int n_u_dofs = context.get_dof_indices(this->_flow_vars.u()).size();
374 
375  // Element Jacobian * quadrature weights for interior integration.
376  const std::vector<libMesh::Real> &JxW =
377  context.get_element_fe(this->_flow_vars.u())->get_JxW();
378 
379  // The pressure shape functions at interior quadrature points.
380  const std::vector<std::vector<libMesh::RealGradient> >& p_dphi =
381  context.get_element_fe(this->_press_var.p())->get_dphi();
382 
383  // The velocity shape functions at interior quadrature points.
384  const std::vector<std::vector<libMesh::Real> >& u_phi =
385  context.get_element_fe(this->_flow_vars.u())->get_phi();
386 
387  const std::vector<std::vector<libMesh::RealGradient> >& u_dphi =
388  context.get_element_fe(this->_flow_vars.u())->get_dphi();
389 
390  const std::vector<std::vector<libMesh::RealTensor> >& u_d2phi =
391  context.get_element_fe(this->_flow_vars.u())->get_d2phi();
392 
393  libMesh::DenseSubVector<libMesh::Number> &Fp = context.get_elem_residual(this->_press_var.p()); // R_{p}
394 
395  libMesh::DenseSubMatrix<libMesh::Number> &Kpp =
396  context.get_elem_jacobian(this->_press_var.p(), this->_press_var.p()); // J_{pp}
397  libMesh::DenseSubMatrix<libMesh::Number> &Kpu =
398  context.get_elem_jacobian(this->_press_var.p(), this->_flow_vars.u()); // J_{pu}
399  libMesh::DenseSubMatrix<libMesh::Number> &Kpv =
400  context.get_elem_jacobian(this->_press_var.p(), this->_flow_vars.v()); // J_{pv}
401  libMesh::DenseSubMatrix<libMesh::Number> *Kpw = NULL;
402 
403 
404  if(this->_flow_vars.dim() == 3)
405  {
406  Kpw = &context.get_elem_jacobian
407  (this->_press_var.p(), this->_flow_vars.w()); // J_{pw}
408  }
409 
410  unsigned int n_qpoints = context.get_element_qrule().n_points();
411 
412  libMesh::FEBase* fe = context.get_element_fe(this->_flow_vars.u());
413 
414  for (unsigned int qp=0; qp != n_qpoints; qp++)
415  {
416  libMesh::RealGradient g = this->_stab_helper.compute_g( fe, context, qp );
417  libMesh::RealTensor G = this->_stab_helper.compute_G( fe, context, qp );
418 
419  libMesh::RealGradient U( context.interior_value( this->_flow_vars.u(), qp ),
420  context.interior_value( this->_flow_vars.v(), qp ) );
421  if( this->_flow_vars.dim() == 3 )
422  {
423  U(2) = context.interior_value( this->_flow_vars.w(), qp );
424  }
425 
426  libMesh::Real tau_M;
427  libMesh::Real d_tau_M_d_rho;
428  libMesh::Gradient d_tau_M_dU;
429  libMesh::Gradient RM_s, d_RM_s_uvw_dgraduvw;
430  libMesh::Tensor d_RM_s_dgradp, d_RM_s_dU, d_RM_s_uvw_dhessuvw;
431 
432  // Compute the viscosity at this qp
433  libMesh::Real _mu_qp = this->_mu(context, qp);
434 
435  if (compute_jacobian)
436  {
437  this->_stab_helper.compute_tau_momentum_and_derivs
438  ( context, qp, g, G, this->_rho, U, _mu_qp,
439  tau_M, d_tau_M_d_rho, d_tau_M_dU,
440  this->_is_steady );
441  this->_stab_helper.compute_res_momentum_steady_and_derivs
442  ( context, qp, this->_rho, _mu_qp,
443  RM_s, d_RM_s_dgradp, d_RM_s_dU, d_RM_s_uvw_dgraduvw,
444  d_RM_s_uvw_dhessuvw);
445  }
446  else
447  {
448  tau_M = this->_stab_helper.compute_tau_momentum( context, qp, g, G, this->_rho, U, _mu_qp, this->_is_steady );
449  RM_s = this->_stab_helper.compute_res_momentum_steady( context, qp, this->_rho, _mu_qp );
450  }
451 
452  // Now a loop over the pressure degrees of freedom. This
453  // computes the contributions of the continuity equation.
454  for (unsigned int i=0; i != n_p_dofs; i++)
455  {
456  Fp(i) += tau_M*(RM_s*p_dphi[i][qp])*JxW[qp];
457 
458  if (compute_jacobian)
459  {
460  const libMesh::Real fixed_deriv =
461  context.get_fixed_solution_derivative();
462 
463  const libMesh::TypeVector<libMesh::Number>
464  p_dphiiT_d_RM_s_dgradp =
465  d_RM_s_dgradp.transpose() * p_dphi[i][qp];
466 
467  for (unsigned int j=0; j != n_p_dofs; j++)
468  {
469  Kpp(i,j) += tau_M*(p_dphiiT_d_RM_s_dgradp*p_dphi[j][qp])*JxW[qp];
470  }
471 
472  for (unsigned int j=0; j != n_u_dofs; j++)
473  {
474  Kpu(i,j) += (
475  d_tau_M_dU(0)*u_phi[j][qp]*(RM_s*p_dphi[i][qp])
476  + tau_M*
477  (
478  d_RM_s_dU(0,0)*u_phi[j][qp]*p_dphi[i][qp](0) +
479  d_RM_s_dU(1,0)*u_phi[j][qp]*p_dphi[i][qp](1) +
480  d_RM_s_dU(2,0)*u_phi[j][qp]*p_dphi[i][qp](2) +
481  d_RM_s_uvw_dgraduvw*u_dphi[j][qp]*p_dphi[i][qp](0) +
482  d_RM_s_uvw_dhessuvw.contract(u_d2phi[j][qp])*p_dphi[i][qp](0)
483  )*fixed_deriv
484  )*JxW[qp];
485  Kpv(i,j) += (
486  d_tau_M_dU(1)*u_phi[j][qp]*(RM_s*p_dphi[i][qp])
487  + tau_M*
488  (
489  d_RM_s_dU(0,1)*u_phi[j][qp]*p_dphi[i][qp](0) +
490  d_RM_s_dU(1,1)*u_phi[j][qp]*p_dphi[i][qp](1) +
491  d_RM_s_dU(2,1)*u_phi[j][qp]*p_dphi[i][qp](2) +
492  d_RM_s_uvw_dgraduvw*u_dphi[j][qp]*p_dphi[i][qp](1) +
493  d_RM_s_uvw_dhessuvw.contract(u_d2phi[j][qp])*p_dphi[i][qp](0)
494  )*fixed_deriv
495  )*JxW[qp];
496  }
497 
498  if(this->_flow_vars.dim() == 3)
499  {
500  for (unsigned int j=0; j != n_u_dofs; j++)
501  {
502  (*Kpw)(i,j) += (
503  d_tau_M_dU(2)*u_phi[j][qp]*(RM_s*p_dphi[i][qp])
504  + tau_M*
505  (
506  d_RM_s_dU(0,2)*u_phi[j][qp]*p_dphi[i][qp](0) +
507  d_RM_s_dU(1,2)*u_phi[j][qp]*p_dphi[i][qp](1) +
508  d_RM_s_dU(2,2)*u_phi[j][qp]*p_dphi[i][qp](2) +
509  d_RM_s_uvw_dgraduvw*u_dphi[j][qp]*p_dphi[i][qp](2) +
510  d_RM_s_uvw_dhessuvw.contract(u_d2phi[j][qp])*p_dphi[i][qp](0)
511  )*fixed_deriv
512  )*JxW[qp];
513  }
514  }
515  }
516  }
517  }
518  }
519 
520  template<class Mu>
521  void IncompressibleNavierStokesAdjointStabilization<Mu>::mass_residual
522  ( bool compute_jacobian, AssemblyContext & context )
523  {
524  // The number of local degrees of freedom in each variable.
525  const unsigned int n_p_dofs = context.get_dof_indices(this->_press_var.p()).size();
526  const unsigned int n_u_dofs = context.get_dof_indices(this->_flow_vars.u()).size();
527 
528  // Element Jacobian * quadrature weights for interior integration.
529  const std::vector<libMesh::Real> &JxW =
530  context.get_element_fe(this->_flow_vars.u())->get_JxW();
531 
532  // The pressure shape functions at interior quadrature points.
533  const std::vector<std::vector<libMesh::RealGradient> >& p_dphi =
534  context.get_element_fe(this->_press_var.p())->get_dphi();
535 
536  const std::vector<std::vector<libMesh::Real> >& u_phi =
537  context.get_element_fe(this->_flow_vars.u())->get_phi();
538 
539  const std::vector<std::vector<libMesh::RealGradient> >& u_gradphi =
540  context.get_element_fe(this->_flow_vars.u())->get_dphi();
541 
542  const std::vector<std::vector<libMesh::RealTensor> >& u_hessphi =
543  context.get_element_fe(this->_flow_vars.u())->get_d2phi();
544 
545  libMesh::DenseSubVector<libMesh::Number> &Fu = context.get_elem_residual(this->_flow_vars.u()); // R_{p}
546  libMesh::DenseSubVector<libMesh::Number> &Fv = context.get_elem_residual(this->_flow_vars.v()); // R_{p}
547  libMesh::DenseSubVector<libMesh::Number> *Fw = NULL;
548 
549  libMesh::DenseSubVector<libMesh::Number> &Fp = context.get_elem_residual(this->_press_var.p()); // R_{p}
550 
551  libMesh::DenseSubMatrix<libMesh::Number> &Kuu =
552  context.get_elem_jacobian(this->_flow_vars.u(), this->_flow_vars.u()); // J_{uu}
553  libMesh::DenseSubMatrix<libMesh::Number> &Kuv =
554  context.get_elem_jacobian(this->_flow_vars.u(), this->_flow_vars.v()); // J_{uv}
555  libMesh::DenseSubMatrix<libMesh::Number> &Kvu =
556  context.get_elem_jacobian(this->_flow_vars.v(), this->_flow_vars.u()); // J_{vu}
557  libMesh::DenseSubMatrix<libMesh::Number> &Kvv =
558  context.get_elem_jacobian(this->_flow_vars.v(), this->_flow_vars.v()); // J_{vv}
559 
560  libMesh::DenseSubMatrix<libMesh::Number> *Kuw = NULL;
561  libMesh::DenseSubMatrix<libMesh::Number> *Kvw = NULL;
562  libMesh::DenseSubMatrix<libMesh::Number> *Kwu = NULL;
563  libMesh::DenseSubMatrix<libMesh::Number> *Kwv = NULL;
564  libMesh::DenseSubMatrix<libMesh::Number> *Kww = NULL;
565 
566  libMesh::DenseSubMatrix<libMesh::Number> &Kpu =
567  context.get_elem_jacobian(this->_press_var.p(), this->_flow_vars.u()); // J_{pu}
568  libMesh::DenseSubMatrix<libMesh::Number> &Kpv =
569  context.get_elem_jacobian(this->_press_var.p(), this->_flow_vars.v()); // J_{pv}
570  libMesh::DenseSubMatrix<libMesh::Number> *Kpw = NULL;
571 
572 
573  if(this->_flow_vars.dim() == 3)
574  {
575  Fw = &context.get_elem_residual(this->_flow_vars.w()); // R_{w}
576 
577  Kuw = &context.get_elem_jacobian
578  (this->_flow_vars.u(), this->_flow_vars.w()); // J_{uw}
579  Kvw = &context.get_elem_jacobian
580  (this->_flow_vars.v(), this->_flow_vars.w()); // J_{vw}
581  Kwu = &context.get_elem_jacobian
582  (this->_flow_vars.w(), this->_flow_vars.u()); // J_{wu}
583  Kwv = &context.get_elem_jacobian
584  (this->_flow_vars.w(), this->_flow_vars.v()); // J_{wv}
585  Kww = &context.get_elem_jacobian
586  (this->_flow_vars.w(), this->_flow_vars.w()); // J_{ww}
587 
588  Kpw = &context.get_elem_jacobian
589  (this->_press_var.p(), this->_flow_vars.w()); // J_{pw}
590  }
591 
592  unsigned int n_qpoints = context.get_element_qrule().n_points();
593 
594  for (unsigned int qp=0; qp != n_qpoints; qp++)
595  {
596  libMesh::FEBase* fe = context.get_element_fe(this->_flow_vars.u());
597 
598  libMesh::RealGradient g = this->_stab_helper.compute_g( fe, context, qp );
599  libMesh::RealTensor G = this->_stab_helper.compute_G( fe, context, qp );
600 
601  libMesh::RealGradient U( context.fixed_interior_value( this->_flow_vars.u(), qp ),
602  context.fixed_interior_value( this->_flow_vars.v(), qp ) );
603  if( this->_flow_vars.dim() == 3 )
604  {
605  U(2) = context.fixed_interior_value( this->_flow_vars.w(), qp );
606  }
607 
608  /*
609  libMesh::Gradient grad_u, grad_v, grad_w;
610  grad_u = context.interior_gradient(this->_flow_vars.u(), qp);
611  grad_v = context.interior_gradient(this->_flow_vars.v(), qp);
612  if (_dim == 3)
613  grad_w = context.interior_gradient(this->_flow_vars.w(), qp);
614  */
615 
616  libMesh::Real tau_M;
617 
618  libMesh::Real d_tau_M_d_rho;
619 
620  libMesh::Gradient d_tau_M_dU;
621 
622  libMesh::RealGradient RM_t;
623 
624  libMesh::Real d_RM_t_uvw_duvw;
625 
626  // Compute the viscosity at this qp
627  libMesh::Real _mu_qp = this->_mu(context, qp);
628 
629  if (compute_jacobian)
630  {
631  this->_stab_helper.compute_tau_momentum_and_derivs
632  ( context, qp, g, G, this->_rho, U, _mu_qp,
633  tau_M, d_tau_M_d_rho, d_tau_M_dU,
634  false );
635  this->_stab_helper.compute_res_momentum_transient_and_derivs
636  ( context, qp, this->_rho,
637  RM_t, d_RM_t_uvw_duvw );
638  }
639  else
640  {
641  tau_M = this->_stab_helper.compute_tau_momentum( context, qp, g, G, this->_rho, U, _mu_qp, false );
642  RM_t = this->_stab_helper.compute_res_momentum_transient( context, qp, this->_rho );
643  }
644 
645 
646  // Now a loop over the pressure degrees of freedom. This
647  // computes the contributions of the continuity equation.
648  for (unsigned int i=0; i != n_p_dofs; i++)
649  {
650  Fp(i) -= tau_M*RM_t*p_dphi[i][qp]*JxW[qp];
651 
652  if (compute_jacobian)
653  {
654  const libMesh::Real fixed_deriv =
655  context.get_fixed_solution_derivative();
656 
657  for (unsigned int j=0; j != n_u_dofs; j++)
658  {
659  Kpu(i,j) -= d_tau_M_dU(0)*u_phi[j][qp]*RM_t*p_dphi[i][qp]*fixed_deriv*JxW[qp];
660  Kpu(i,j) -= tau_M*d_RM_t_uvw_duvw*u_phi[j][qp]*p_dphi[i][qp](0)*fixed_deriv*JxW[qp];
661 
662  Kpv(i,j) -= d_tau_M_dU(1)*u_phi[j][qp]*RM_t*p_dphi[i][qp]*fixed_deriv*JxW[qp];
663  Kpv(i,j) -= tau_M*d_RM_t_uvw_duvw*u_phi[j][qp]*p_dphi[i][qp](1)*fixed_deriv*JxW[qp];
664  }
665 
666  if(this->_flow_vars.dim() == 3)
667  {
668  for (unsigned int j=0; j != n_u_dofs; j++)
669  {
670  (*Kpw)(i,j) -= d_tau_M_dU(2)*u_phi[j][qp]*RM_t*p_dphi[i][qp]*fixed_deriv*JxW[qp];
671  (*Kpw)(i,j) -= tau_M*d_RM_t_uvw_duvw*u_phi[j][qp]*p_dphi[i][qp](2)*fixed_deriv*JxW[qp];
672  }
673  }
674  }
675  }
676 
677  for (unsigned int i=0; i != n_u_dofs; i++)
678  {
679  libMesh::Real test_func = this->_rho*U*u_gradphi[i][qp] +
680  _mu_qp*( u_hessphi[i][qp](0,0) + u_hessphi[i][qp](1,1) + u_hessphi[i][qp](2,2) );
681  libMesh::Gradient d_test_func_dU = this->_rho*u_gradphi[i][qp];
682 
683  //libMesh::RealGradient zeroth_order_term = - this->_rho*u_phi[i][qp]*(grad_u + grad_v + grad_w);
684 
685  Fu(i) -= tau_M*RM_t(0)*test_func*JxW[qp];
686 
687  Fv(i) -= tau_M*RM_t(1)*test_func*JxW[qp];
688 
689  if(this->_flow_vars.dim() == 3)
690  {
691  (*Fw)(i) -= tau_M*RM_t(2)*test_func*JxW[qp];
692  }
693 
694  if (compute_jacobian)
695  {
696  const libMesh::Real fixed_deriv =
697  context.get_fixed_solution_derivative();
698 
699  for (unsigned int j=0; j != n_u_dofs; j++)
700  {
701  Kuu(i,j) -= d_tau_M_dU(0)*u_phi[j][qp]*RM_t(0)*test_func*fixed_deriv*JxW[qp];
702  Kuu(i,j) -= tau_M*d_RM_t_uvw_duvw*u_phi[j][qp]*test_func*fixed_deriv*JxW[qp];
703  Kuu(i,j) -= tau_M*RM_t(0)*d_test_func_dU(0)*u_phi[j][qp]*fixed_deriv*JxW[qp];
704 
705  Kuv(i,j) -= d_tau_M_dU(1)*u_phi[j][qp]*RM_t(0)*test_func*fixed_deriv*JxW[qp];
706  Kuv(i,j) -= tau_M*RM_t(0)*d_test_func_dU(1)*u_phi[j][qp]*fixed_deriv*JxW[qp];
707 
708  Kvu(i,j) -= d_tau_M_dU(0)*u_phi[j][qp]*RM_t(1)*test_func*fixed_deriv*JxW[qp];
709  Kvu(i,j) -= tau_M*RM_t(1)*d_test_func_dU(0)*u_phi[j][qp]*fixed_deriv*JxW[qp];
710 
711  Kvv(i,j) -= d_tau_M_dU(1)*u_phi[j][qp]*RM_t(1)*test_func*fixed_deriv*JxW[qp];
712  Kvv(i,j) -= tau_M*d_RM_t_uvw_duvw*u_phi[j][qp]*test_func*fixed_deriv*JxW[qp];
713  Kvv(i,j) -= tau_M*RM_t(1)*d_test_func_dU(1)*u_phi[j][qp]*fixed_deriv*JxW[qp];
714  }
715  if(this->_flow_vars.dim() == 3)
716  {
717  for (unsigned int j=0; j != n_u_dofs; j++)
718  {
719  (*Kuw)(i,j) -= d_tau_M_dU(2)*u_phi[j][qp]*RM_t(0)*test_func*fixed_deriv*JxW[qp];
720  (*Kuw)(i,j) -= tau_M*RM_t(0)*d_test_func_dU(2)*u_phi[j][qp]*fixed_deriv*JxW[qp];
721 
722  (*Kvw)(i,j) -= d_tau_M_dU(2)*u_phi[j][qp]*RM_t(1)*test_func*fixed_deriv*JxW[qp];
723  (*Kvw)(i,j) -= tau_M*RM_t(1)*d_test_func_dU(2)*u_phi[j][qp]*fixed_deriv*JxW[qp];
724 
725  (*Kwu)(i,j) -= d_tau_M_dU(0)*u_phi[j][qp]*RM_t(2)*test_func*fixed_deriv*JxW[qp];
726  (*Kwu)(i,j) -= tau_M*RM_t(2)*d_test_func_dU(0)*u_phi[j][qp]*fixed_deriv*JxW[qp];
727 
728  (*Kwv)(i,j) -= d_tau_M_dU(1)*u_phi[j][qp]*RM_t(2)*test_func*fixed_deriv*JxW[qp];
729  (*Kwv)(i,j) -= tau_M*RM_t(2)*d_test_func_dU(1)*u_phi[j][qp]*fixed_deriv*JxW[qp];
730 
731  (*Kww)(i,j) -= d_tau_M_dU(2)*u_phi[j][qp]*RM_t(2)*test_func*fixed_deriv*JxW[qp];
732  (*Kww)(i,j) -= tau_M*d_RM_t_uvw_duvw*u_phi[j][qp]*test_func*fixed_deriv*JxW[qp];
733  (*Kww)(i,j) -= tau_M*RM_t(2)*d_test_func_dU(2)*u_phi[j][qp]*fixed_deriv*JxW[qp];
734  }
735  }
736  }
737  }
738 
739  }
740  }
741 
742 } // namespace GRINS
743 
744 // Instantiate
745 INSTANTIATE_INC_NS_SUBCLASS(IncompressibleNavierStokesAdjointStabilization);
GRINS::IncompressibleNavierStokesAdjointStabilization::element_time_derivative
virtual void element_time_derivative(bool compute_jacobian, AssemblyContext &context)
Time dependent part(s) of physics for element interiors.
Definition: inc_navier_stokes_adjoint_stab.C:47
assembly_context.h
GRINS::IncompressibleNavierStokesAdjointStabilization::mass_residual
virtual void mass_residual(bool compute_jacobian, AssemblyContext &context)
Mass matrix part(s) for element interiors. All boundary terms lie within the time_derivative part...
Definition: inc_navier_stokes_adjoint_stab.C:522
GRINS
GRINS namespace.
Definition: arrhenius_catalycity.h:31
GRINS::IncompressibleNavierStokesAdjointStabilization::element_constraint
virtual void element_constraint(bool compute_jacobian, AssemblyContext &context)
Constraint part(s) of physics for element interiors.
Definition: inc_navier_stokes_adjoint_stab.C:369
inc_navier_stokes_adjoint_stab.h
INSTANTIATE_INC_NS_SUBCLASS
INSTANTIATE_INC_NS_SUBCLASS(IncompressibleNavierStokesAdjointStabilization)
inc_nav_stokes_macro.h
Generated on Tue Dec 19 2017 12:47:28 for GRINS-0.8.0 by  doxygen 1.8.9.1