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