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