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