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