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