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