GRINS-0.7.0
velocity_drag_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_config.h"
28 #include "grins/velocity_drag_adjoint_stab.h"
29 
30 // GRINS
31 #include "grins/assembly_context.h"
32 #include "grins/inc_nav_stokes_macro.h"
33 
34 // libMesh
35 #include "libmesh/quadrature.h"
36 
37 namespace GRINS
38 {
39 
40  template<class Mu>
41  VelocityDragAdjointStabilization<Mu>::VelocityDragAdjointStabilization( const std::string& physics_name, const GetPot& input )
42  : VelocityDragBase<Mu>(physics_name, input),
43  _stab_helper( physics_name+"StabHelper", input )
44  {}
45 
46  template<class Mu>
47  VelocityDragAdjointStabilization<Mu>::~VelocityDragAdjointStabilization()
48  {
49  return;
50  }
51 
52  template<class Mu>
53  void VelocityDragAdjointStabilization<Mu>::init_context( AssemblyContext& context )
54  {
55  context.get_element_fe(this->_press_var.p())->get_dphi();
56 
57  context.get_element_fe(this->_flow_vars.u())->get_xyz();
58  context.get_element_fe(this->_flow_vars.u())->get_phi();
59  context.get_element_fe(this->_flow_vars.u())->get_dphi();
60  context.get_element_fe(this->_flow_vars.u())->get_d2phi();
61 
62  return;
63  }
64 
65  template<class Mu>
66  void VelocityDragAdjointStabilization<Mu>::element_time_derivative
67  ( bool compute_jacobian,
68  AssemblyContext& context,
69  CachedValues& /* cache */ )
70  {
71 #ifdef GRINS_USE_GRVY_TIMERS
72  this->_timer->BeginTimer("VelocityDragAdjointStabilization::element_time_derivative");
73 #endif
74 
75  libMesh::FEBase* fe = context.get_element_fe(this->_flow_vars.u());
76 
77  // Element Jacobian * quadrature weights for interior integration
78  const std::vector<libMesh::Real> &JxW = fe->get_JxW();
79 
80  // The shape functions at interior quadrature points.
81  const std::vector<std::vector<libMesh::Real> >& u_phi = fe->get_phi();
82 
83  const std::vector<std::vector<libMesh::RealGradient> >& u_gradphi =
84  context.get_element_fe(this->_flow_vars.u())->get_dphi();
85 
86  const std::vector<std::vector<libMesh::RealTensor> >& u_hessphi =
87  context.get_element_fe(this->_flow_vars.u())->get_d2phi();
88 
89  const std::vector<libMesh::Point>& u_qpoint = fe->get_xyz();
90 
91  // The number of local degrees of freedom in each variable
92  const unsigned int n_u_dofs = context.get_dof_indices(this->_flow_vars.u()).size();
93 
94  // The subvectors and submatrices we need to fill:
95  libMesh::DenseSubMatrix<libMesh::Number> &Kuu = context.get_elem_jacobian(this->_flow_vars.u(), this->_flow_vars.u()); // R_{u},{u}
96  libMesh::DenseSubMatrix<libMesh::Number> &Kuv = context.get_elem_jacobian(this->_flow_vars.u(), this->_flow_vars.v()); // R_{u},{v}
97  libMesh::DenseSubMatrix<libMesh::Number> &Kvu = context.get_elem_jacobian(this->_flow_vars.v(), this->_flow_vars.u()); // R_{v},{u}
98  libMesh::DenseSubMatrix<libMesh::Number> &Kvv = context.get_elem_jacobian(this->_flow_vars.v(), this->_flow_vars.v()); // R_{v},{v}
99 
100  libMesh::DenseSubMatrix<libMesh::Number>* Kwu = NULL;
101  libMesh::DenseSubMatrix<libMesh::Number>* Kwv = NULL;
102  libMesh::DenseSubMatrix<libMesh::Number>* Kww = NULL;
103  libMesh::DenseSubMatrix<libMesh::Number>* Kuw = NULL;
104  libMesh::DenseSubMatrix<libMesh::Number>* Kvw = NULL;
105 
106  libMesh::DenseSubVector<libMesh::Number> &Fu = context.get_elem_residual(this->_flow_vars.u()); // R_{u}
107  libMesh::DenseSubVector<libMesh::Number> &Fv = context.get_elem_residual(this->_flow_vars.v()); // R_{v}
108  libMesh::DenseSubVector<libMesh::Number>* Fw = NULL;
109 
110  if( this->_dim == 3 )
111  {
112  Kuw = &context.get_elem_jacobian(this->_flow_vars.u(), this->_flow_vars.w()); // R_{u},{w}
113  Kvw = &context.get_elem_jacobian(this->_flow_vars.v(), this->_flow_vars.w()); // R_{v},{w}
114 
115  Kwu = &context.get_elem_jacobian(this->_flow_vars.w(), this->_flow_vars.u()); // R_{w},{u}
116  Kwv = &context.get_elem_jacobian(this->_flow_vars.w(), this->_flow_vars.v()); // R_{w},{v}
117  Kww = &context.get_elem_jacobian(this->_flow_vars.w(), this->_flow_vars.w()); // R_{w},{w}
118  Fw = &context.get_elem_residual(this->_flow_vars.w()); // R_{w}
119  }
120 
121  unsigned int n_qpoints = context.get_element_qrule().n_points();
122 
123  for (unsigned int qp=0; qp != n_qpoints; qp++)
124  {
125  libMesh::RealGradient g = this->_stab_helper.compute_g( fe, context, qp );
126  libMesh::RealTensor G = this->_stab_helper.compute_G( fe, context, qp );
127 
128  libMesh::RealGradient U( context.interior_value( this->_flow_vars.u(), qp ),
129  context.interior_value( this->_flow_vars.v(), qp ) );
130  if( this->_dim == 3 )
131  {
132  U(2) = context.interior_value( this->_flow_vars.w(), qp );
133  }
134 
135  // Compute the viscosity at this qp
136  libMesh::Real mu_qp = this->_mu(context, qp);
137 
138  libMesh::Real tau_M;
139  libMesh::Real d_tau_M_d_rho;
140  libMesh::Gradient d_tau_M_dU;
141 
142  if (compute_jacobian)
143  this->_stab_helper.compute_tau_momentum_and_derivs
144  ( context, qp, g, G, this->_rho, U, mu_qp,
145  tau_M, d_tau_M_d_rho, d_tau_M_dU,
146  this->_is_steady );
147  else
148  tau_M = this->_stab_helper.compute_tau_momentum
149  ( context, qp, g, G, this->_rho, U, mu_qp,
150  this->_is_steady );
151 
152  libMesh::NumberVectorValue F;
153  libMesh::NumberTensorValue dFdU;
154  libMesh::NumberTensorValue* dFdU_ptr =
155  compute_jacobian ? &dFdU : NULL;
156  if (!this->compute_force(u_qpoint[qp], context.time, U, F, dFdU_ptr))
157  continue;
158 
159  for (unsigned int i=0; i != n_u_dofs; i++)
160  {
161  libMesh::Real test_func = this->_rho*U*u_gradphi[i][qp] +
162  mu_qp*( u_hessphi[i][qp](0,0) + u_hessphi[i][qp](1,1) + u_hessphi[i][qp](2,2) );
163  Fu(i) += tau_M*F(0)*test_func*JxW[qp];
164 
165  Fv(i) += tau_M*F(1)*test_func*JxW[qp];
166 
167  if (this->_dim == 3)
168  {
169  (*Fw)(i) += tau_M*F(2)*test_func*JxW[qp];
170  }
171 
172  if (compute_jacobian)
173  {
174  const libMesh::Real sol_deriv =
175  context.get_elem_solution_derivative();
176 
177  const libMesh::Real JxWxS = JxW[qp] * sol_deriv;
178 
179  const libMesh::Real fixed_deriv =
180  context.get_fixed_solution_derivative();
181 
182  const libMesh::Real JxWxF = JxW[qp] * fixed_deriv;
183 
184  libMesh::Gradient d_test_func_dU = this->_rho*u_gradphi[i][qp];
185 
186  for (unsigned int j=0; j != n_u_dofs; ++j)
187  {
188  Kuu(i,j) += tau_M*F(0)*d_test_func_dU(0)*u_phi[j][qp]*JxWxS;
189  Kuu(i,j) += d_tau_M_dU(0)*u_phi[j][qp]*F(0)*test_func*JxWxF;
190  Kuu(i,j) += tau_M*dFdU(0,0)*u_phi[j][qp]*test_func*JxWxS;
191  Kuv(i,j) += tau_M*F(0)*d_test_func_dU(1)*u_phi[j][qp]*JxWxS;
192  Kuv(i,j) += d_tau_M_dU(1)*u_phi[j][qp]*F(0)*test_func*JxWxF;
193  Kuv(i,j) += tau_M*dFdU(0,1)*u_phi[j][qp]*test_func*JxWxS;
194  Kvu(i,j) += tau_M*F(1)*d_test_func_dU(0)*u_phi[j][qp]*JxWxS;
195  Kvu(i,j) += d_tau_M_dU(0)*u_phi[j][qp]*F(1)*test_func*JxWxF;
196  Kvu(i,j) += tau_M*dFdU(1,0)*u_phi[j][qp]*test_func*JxWxS;
197  Kvv(i,j) += tau_M*F(1)*d_test_func_dU(1)*u_phi[j][qp]*JxWxS;
198  Kvv(i,j) += d_tau_M_dU(1)*u_phi[j][qp]*F(1)*test_func*JxWxF;
199  Kvv(i,j) += tau_M*dFdU(1,1)*u_phi[j][qp]*test_func*JxWxS;
200  }
201 
202  if (this->_dim == 3)
203  {
204  for (unsigned int j=0; j != n_u_dofs; ++j)
205  {
206  (*Kuw)(i,j) += tau_M*F(0)*d_test_func_dU(2)*u_phi[j][qp]*JxWxS;
207  (*Kuw)(i,j) += d_tau_M_dU(2)*u_phi[j][qp]*F(0)*test_func*JxWxF;
208  (*Kuw)(i,j) += tau_M*dFdU(0,2)*u_phi[j][qp]*test_func*JxWxS;
209  (*Kvw)(i,j) += tau_M*F(1)*d_test_func_dU(2)*u_phi[j][qp]*JxWxS;
210  (*Kvw)(i,j) += d_tau_M_dU(2)*u_phi[j][qp]*F(1)*test_func*JxWxF;
211  (*Kvw)(i,j) += tau_M*dFdU(1,2)*u_phi[j][qp]*test_func*JxWxS;
212  (*Kwu)(i,j) += tau_M*F(2)*d_test_func_dU(0)*u_phi[j][qp]*JxWxS;
213  (*Kwu)(i,j) += d_tau_M_dU(0)*u_phi[j][qp]*F(2)*test_func*JxWxF;
214  (*Kwu)(i,j) += tau_M*dFdU(2,0)*u_phi[j][qp]*test_func*JxWxS;
215  (*Kwv)(i,j) += tau_M*F(2)*d_test_func_dU(1)*u_phi[j][qp]*JxWxS;
216  (*Kwv)(i,j) += d_tau_M_dU(1)*u_phi[j][qp]*F(2)*test_func*JxWxF;
217  (*Kwv)(i,j) += tau_M*dFdU(2,1)*u_phi[j][qp]*test_func*JxWxS;
218  (*Kww)(i,j) += tau_M*F(2)*d_test_func_dU(2)*u_phi[j][qp]*JxWxS;
219  (*Kww)(i,j) += d_tau_M_dU(2)*u_phi[j][qp]*F(2)*test_func*JxWxF;
220  (*Kww)(i,j) += tau_M*dFdU(2,2)*u_phi[j][qp]*test_func*JxWxS;
221  }
222  }
223 
224  } // End compute_jacobian check
225 
226  } // End i dof loop
227  }
228 
229 #ifdef GRINS_USE_GRVY_TIMERS
230  this->_timer->EndTimer("VelocityDragAdjointStabilization::element_time_derivative");
231 #endif
232 
233  return;
234  }
235 
236  template<class Mu>
237  void VelocityDragAdjointStabilization<Mu>::element_constraint( bool compute_jacobian,
238  AssemblyContext& context,
239  CachedValues& /*cache*/ )
240  {
241 #ifdef GRINS_USE_GRVY_TIMERS
242  this->_timer->BeginTimer("VelocityDragAdjointStabilization::element_constraint");
243 #endif
244 
245  // The number of local degrees of freedom in each variable.
246  const unsigned int n_p_dofs = context.get_dof_indices(this->_press_var.p()).size();
247  const unsigned int n_u_dofs = context.get_dof_indices(this->_flow_vars.u()).size();
248 
249  // Element Jacobian * quadrature weights for interior integration.
250  const std::vector<libMesh::Real> &JxW =
251  context.get_element_fe(this->_flow_vars.u())->get_JxW();
252 
253  const std::vector<libMesh::Point>& u_qpoint =
254  context.get_element_fe(this->_flow_vars.u())->get_xyz();
255 
256  const std::vector<std::vector<libMesh::Real> >& u_phi =
257  context.get_element_fe(this->_flow_vars.u())->get_phi();
258 
259  const std::vector<std::vector<libMesh::RealGradient> >& p_dphi =
260  context.get_element_fe(this->_press_var.p())->get_dphi();
261 
262  libMesh::DenseSubVector<libMesh::Number> &Fp = context.get_elem_residual(this->_press_var.p()); // R_{p}
263 
264  libMesh::DenseSubMatrix<libMesh::Number> &Kpu =
265  context.get_elem_jacobian(this->_press_var.p(), this->_flow_vars.u()); // J_{pu}
266  libMesh::DenseSubMatrix<libMesh::Number> &Kpv =
267  context.get_elem_jacobian(this->_press_var.p(), this->_flow_vars.v()); // J_{pv}
268  libMesh::DenseSubMatrix<libMesh::Number> *Kpw = NULL;
269 
270  if(this->_dim == 3)
271  {
272  Kpw = &context.get_elem_jacobian
273  (this->_press_var.p(), this->_flow_vars.w()); // J_{pw}
274  }
275 
276  // Now we will build the element Jacobian and residual.
277  // Constructing the residual requires the solution and its
278  // gradient from the previous timestep. This must be
279  // calculated at each quadrature point by summing the
280  // solution degree-of-freedom values by the appropriate
281  // weight functions.
282  unsigned int n_qpoints = context.get_element_qrule().n_points();
283 
284  libMesh::FEBase* fe = context.get_element_fe(this->_flow_vars.u());
285 
286  for (unsigned int qp=0; qp != n_qpoints; qp++)
287  {
288  libMesh::RealGradient g = this->_stab_helper.compute_g( fe, context, qp );
289  libMesh::RealTensor G = this->_stab_helper.compute_G( fe, context, qp );
290 
291  libMesh::RealGradient U( context.interior_value( this->_flow_vars.u(), qp ),
292  context.interior_value( this->_flow_vars.v(), qp ) );
293  if( this->_dim == 3 )
294  {
295  U(2) = context.interior_value( this->_flow_vars.w(), qp );
296  }
297 
298  // Compute the viscosity at this qp
299  libMesh::Real mu_qp = this->_mu(context, qp);
300 
301  libMesh::Real tau_M;
302  libMesh::Real d_tau_M_d_rho;
303  libMesh::Gradient d_tau_M_dU;
304 
305  if (compute_jacobian)
306  this->_stab_helper.compute_tau_momentum_and_derivs
307  ( context, qp, g, G, this->_rho, U, mu_qp,
308  tau_M, d_tau_M_d_rho, d_tau_M_dU,
309  this->_is_steady );
310  else
311  tau_M = this->_stab_helper.compute_tau_momentum
312  ( context, qp, g, G, this->_rho, U, mu_qp,
313  this->_is_steady );
314 
315  libMesh::NumberVectorValue F;
316  libMesh::NumberTensorValue dFdU;
317  libMesh::NumberTensorValue* dFdU_ptr =
318  compute_jacobian ? &dFdU : NULL;
319  if (!this->compute_force(u_qpoint[qp], context.time, U, F, dFdU_ptr))
320  continue;
321 
322  // First, an i-loop over the velocity degrees of freedom.
323  // We know that n_u_dofs == n_v_dofs so we can compute contributions
324  // for both at the same time.
325  for (unsigned int i=0; i != n_p_dofs; i++)
326  {
327  Fp(i) += -tau_M*F*p_dphi[i][qp]*JxW[qp];
328 
329  if (compute_jacobian)
330  {
331  const libMesh::Real sol_deriv =
332  context.get_elem_solution_derivative();
333 
334  const libMesh::Real JxWxS = JxW[qp] * sol_deriv;
335 
336  const libMesh::Real fixed_deriv =
337  context.get_fixed_solution_derivative();
338 
339  const libMesh::Real JxWxF = JxW[qp] * fixed_deriv;
340 
341  for (unsigned int j=0; j != n_u_dofs; ++j)
342  {
343  Kpu(i,j) += -d_tau_M_dU(0)*u_phi[j][qp]*F*p_dphi[i][qp]*JxWxF;
344  Kpv(i,j) += -d_tau_M_dU(1)*u_phi[j][qp]*F*p_dphi[i][qp]*JxWxF;
345  for (unsigned int d=0; d != 3; ++d)
346  {
347  Kpu(i,j) += -tau_M*dFdU(d,0)*u_phi[j][qp]*p_dphi[i][qp](d)*JxWxS;
348  Kpv(i,j) += -tau_M*dFdU(d,1)*u_phi[j][qp]*p_dphi[i][qp](d)*JxWxS;
349  }
350  }
351  if( this->_dim == 3 )
352  for (unsigned int j=0; j != n_u_dofs; ++j)
353  {
354  (*Kpw)(i,j) += -d_tau_M_dU(2)*u_phi[j][qp]*F*p_dphi[i][qp]*JxWxF;
355  for (unsigned int d=0; d != 3; ++d)
356  {
357  (*Kpw)(i,j) += -tau_M*dFdU(d,2)*u_phi[j][qp]*p_dphi[i][qp](d)*JxWxS;
358  }
359  }
360  }
361  }
362  } // End quadrature loop
363 
364 #ifdef GRINS_USE_GRVY_TIMERS
365  this->_timer->EndTimer("VelocityDragAdjointStabilization::element_constraint");
366 #endif
367 
368  return;
369  }
370 
371 
372 } // namespace GRINS
373 
374 // Instantiate
375 INSTANTIATE_INC_NS_SUBCLASS(VelocityDragAdjointStabilization);
INSTANTIATE_INC_NS_SUBCLASS
INSTANTIATE_INC_NS_SUBCLASS(VelocityDragAdjointStabilization)
assembly_context.h
velocity_drag_adjoint_stab.h
GRINS::VelocityDragAdjointStabilization::element_constraint
virtual void element_constraint(bool compute_jacobian, AssemblyContext &context, CachedValues &cache)
Constraint part(s) of physics for element interiors.
Definition: velocity_drag_adjoint_stab.C:237
GRINS
GRINS namespace.
Definition: arrhenius_catalycity.h:31
GRINS::VelocityDragAdjointStabilization::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: velocity_drag_adjoint_stab.C:67
inc_nav_stokes_macro.h
GRINS::VelocityDragAdjointStabilization::init_context
virtual void init_context(AssemblyContext &context)
Initialize context for added physics variables.
Definition: velocity_drag_adjoint_stab.C:53
Generated on Thu Jun 2 2016 21:52:28 for GRINS-0.7.0 by  doxygen 1.8.10