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