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