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