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