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