GRINS-0.8.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-2017 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
69  ( bool compute_jacobian, AssemblyContext & context )
70  {
71  // The number of local degrees of freedom in each variable.
72  const unsigned int n_u_dofs = context.get_dof_indices(_flow_vars.u()).size();
73  const unsigned int n_T_dofs = context.get_dof_indices(_temp_vars.T()).size();
74 
75  // Element Jacobian * quadrature weights for interior integration.
76  const std::vector<libMesh::Real> &JxW =
77  context.get_element_fe(_flow_vars.u())->get_JxW();
78 
79  const std::vector<std::vector<libMesh::Real> >& T_phi =
80  context.get_element_fe(this->_temp_vars.T())->get_phi();
81 
82  const std::vector<std::vector<libMesh::Real> >& u_phi =
83  context.get_element_fe(this->_flow_vars.u())->get_phi();
84 
85  const std::vector<std::vector<libMesh::RealGradient> >& u_gradphi =
86  context.get_element_fe(this->_flow_vars.u())->get_dphi();
87 
88  const std::vector<std::vector<libMesh::RealTensor> >& u_hessphi =
89  context.get_element_fe(this->_flow_vars.u())->get_d2phi();
90 
91  // Get residuals and jacobians
92  libMesh::DenseSubVector<libMesh::Number> &Fu = context.get_elem_residual(_flow_vars.u()); // R_{u}
93  libMesh::DenseSubVector<libMesh::Number> &Fv = context.get_elem_residual(_flow_vars.v()); // R_{v}
94  libMesh::DenseSubVector<libMesh::Number> *Fw = NULL;
95 
96  libMesh::DenseSubMatrix<libMesh::Number> &KuT =
97  context.get_elem_jacobian(_flow_vars.u(), _temp_vars.T()); // J_{uT}
98  libMesh::DenseSubMatrix<libMesh::Number> &KvT =
99  context.get_elem_jacobian(_flow_vars.v(), _temp_vars.T()); // J_{vT}
100  libMesh::DenseSubMatrix<libMesh::Number> &Kuu =
101  context.get_elem_jacobian(_flow_vars.u(), _flow_vars.u()); // J_{uu}
102  libMesh::DenseSubMatrix<libMesh::Number> &Kuv =
103  context.get_elem_jacobian(_flow_vars.u(), _flow_vars.v()); // J_{uv}
104  libMesh::DenseSubMatrix<libMesh::Number> &Kvu =
105  context.get_elem_jacobian(_flow_vars.v(), _flow_vars.u()); // J_{vu}
106  libMesh::DenseSubMatrix<libMesh::Number> &Kvv =
107  context.get_elem_jacobian(_flow_vars.v(), _flow_vars.v()); // J_{vv}
108 
109  libMesh::DenseSubMatrix<libMesh::Number> *KwT = NULL;
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 
117  if(this->_flow_vars.dim() == 3)
118  {
119  Fw = &context.get_elem_residual(this->_flow_vars.w()); // R_{w}
120  KwT = &context.get_elem_jacobian
121  (_flow_vars.w(), _temp_vars.T()); // J_{wT}
122  Kuw = &context.get_elem_jacobian
123  (_flow_vars.u(), _flow_vars.w()); // J_{uw}
124  Kvw = &context.get_elem_jacobian
125  (_flow_vars.v(), _flow_vars.w()); // J_{vw}
126  Kwu = &context.get_elem_jacobian
127  (_flow_vars.w(), _flow_vars.u()); // J_{wu}
128  Kwv = &context.get_elem_jacobian
129  (_flow_vars.w(), _flow_vars.v()); // J_{wv}
130  Kww = &context.get_elem_jacobian
131  (_flow_vars.w(), _flow_vars.w()); // J_{ww}
132  }
133 
134  // Now we will build the element Jacobian and residual.
135  // Constructing the residual requires the solution and its
136  // gradient from the previous timestep. This must be
137  // calculated at each quadrature point by summing the
138  // solution degree-of-freedom values by the appropriate
139  // weight functions.
140  unsigned int n_qpoints = context.get_element_qrule().n_points();
141 
142  libMesh::FEBase* fe = context.get_element_fe(this->_flow_vars.u());
143 
144  for (unsigned int qp=0; qp != n_qpoints; qp++)
145  {
146  libMesh::RealGradient g = this->_stab_helper.compute_g( fe, context, qp );
147  libMesh::RealTensor G = this->_stab_helper.compute_G( fe, context, qp );
148 
149  libMesh::RealGradient U( context.interior_value( this->_flow_vars.u(), qp ),
150  context.interior_value( this->_flow_vars.v(), qp ) );
151  if( this->_flow_vars.dim() == 3 )
152  {
153  U(2) = context.interior_value( this->_flow_vars.w(), qp );
154  }
155 
156  // Compute the viscosity at this qp
157  libMesh::Real mu_qp = this->_mu(context, qp);
158 
159  libMesh::Real tau_M;
160  libMesh::Real d_tau_M_d_rho;
161  libMesh::Gradient d_tau_M_dU;
162 
163  if (compute_jacobian)
164  this->_stab_helper.compute_tau_momentum_and_derivs
165  ( context, qp, g, G, this->_rho, U, mu_qp,
166  tau_M, d_tau_M_d_rho, d_tau_M_dU,
167  this->_is_steady );
168  else
169  tau_M = this->_stab_helper.compute_tau_momentum
170  ( context, qp, g, G, this->_rho, U, mu_qp,
171  this->_is_steady );
172 
173  // Compute the solution & its gradient at the old Newton iterate.
174  libMesh::Number T;
175  T = context.interior_value(_temp_vars.T(), qp);
176 
177  libMesh::RealGradient d_residual_dT = _rho*_beta_T*_g;
178  // d_residual_dU = 0
179  libMesh::RealGradient residual = (T-_T_ref)*d_residual_dT;
180 
181  for (unsigned int i=0; i != n_u_dofs; i++)
182  {
183  libMesh::Real test_func = this->_rho*U*u_gradphi[i][qp] +
184  mu_qp*( u_hessphi[i][qp](0,0) + u_hessphi[i][qp](1,1) + u_hessphi[i][qp](2,2) );
185  Fu(i) += -tau_M*residual(0)*test_func*JxW[qp];
186 
187  Fv(i) += -tau_M*residual(1)*test_func*JxW[qp];
188 
189  if (this->_flow_vars.dim() == 3)
190  {
191  (*Fw)(i) += -tau_M*residual(2)*test_func*JxW[qp];
192  }
193 
194  if (compute_jacobian)
195  {
196  libMesh::Gradient d_test_func_dU = this->_rho*u_gradphi[i][qp];
197  // d_test_func_dT = 0
198 
199  for (unsigned int j=0; j != n_u_dofs; ++j)
200  {
201  Kuu(i,j) += -tau_M*residual(0)*d_test_func_dU(0)*u_phi[j][qp]*JxW[qp] * context.get_elem_solution_derivative();
202  Kuu(i,j) += -d_tau_M_dU(0)*u_phi[j][qp]*residual(0)*test_func*JxW[qp] * context.get_elem_solution_derivative();
203  Kuv(i,j) += -tau_M*residual(0)*d_test_func_dU(1)*u_phi[j][qp]*JxW[qp] * context.get_elem_solution_derivative();
204  Kuv(i,j) += -d_tau_M_dU(1)*u_phi[j][qp]*residual(0)*test_func*JxW[qp] * context.get_elem_solution_derivative();
205  Kvu(i,j) += -tau_M*residual(1)*d_test_func_dU(0)*u_phi[j][qp]*JxW[qp] * context.get_elem_solution_derivative();
206  Kvu(i,j) += -d_tau_M_dU(0)*u_phi[j][qp]*residual(1)*test_func*JxW[qp] * context.get_elem_solution_derivative();
207  Kvv(i,j) += -tau_M*residual(1)*d_test_func_dU(1)*u_phi[j][qp]*JxW[qp] * context.get_elem_solution_derivative();
208  Kvv(i,j) += -d_tau_M_dU(1)*u_phi[j][qp]*residual(1)*test_func*JxW[qp] * context.get_elem_solution_derivative();
209  }
210 
211  for (unsigned int j=0; j != n_T_dofs; ++j)
212  {
213  // KuT(i,j) += -tau_M*residual(0)*dtest_func_dT[j]*JxW[qp] * context.get_elem_solution_derivative();
214  KuT(i,j) += -tau_M*d_residual_dT(0)*T_phi[j][qp]*test_func*JxW[qp] * context.get_elem_solution_derivative();
215  // KvT(i,j) += -tau_M*residual(1)*dtest_func_dT[j]*JxW[qp] * context.get_elem_solution_derivative();
216  KvT(i,j) += -tau_M*d_residual_dT(1)*T_phi[j][qp]*test_func*JxW[qp] * context.get_elem_solution_derivative();
217  }
218  if (this->_flow_vars.dim() == 3)
219  {
220  for (unsigned int j=0; j != n_T_dofs; ++j)
221  {
222  // KwT(i,j) += -tau_M*residual(2)*dtest_func_dT[j]*JxW[qp] * context.get_elem_solution_derivative();
223  (*KwT)(i,j) += -tau_M*d_residual_dT(2)*T_phi[j][qp]*test_func*JxW[qp] * context.get_elem_solution_derivative();
224  }
225 
226  for (unsigned int j=0; j != n_u_dofs; ++j)
227  {
228  (*Kuw)(i,j) += -tau_M*residual(0)*d_test_func_dU(2)*u_phi[j][qp]*JxW[qp] * context.get_elem_solution_derivative();
229  (*Kuw)(i,j) += -d_tau_M_dU(2)*u_phi[j][qp]*residual(0)*test_func*JxW[qp] * context.get_elem_solution_derivative();
230  (*Kvw)(i,j) += -tau_M*residual(1)*d_test_func_dU(2)*u_phi[j][qp]*JxW[qp] * context.get_elem_solution_derivative();
231  (*Kvw)(i,j) += -d_tau_M_dU(2)*u_phi[j][qp]*residual(1)*test_func*JxW[qp] * context.get_elem_solution_derivative();
232  (*Kwu)(i,j) += -tau_M*residual(2)*d_test_func_dU(0)*u_phi[j][qp]*JxW[qp] * context.get_elem_solution_derivative();
233  (*Kwu)(i,j) += -d_tau_M_dU(0)*u_phi[j][qp]*residual(2)*test_func*JxW[qp] * context.get_elem_solution_derivative();
234  (*Kwv)(i,j) += -tau_M*residual(2)*d_test_func_dU(1)*u_phi[j][qp]*JxW[qp] * context.get_elem_solution_derivative();
235  (*Kwv)(i,j) += -d_tau_M_dU(1)*u_phi[j][qp]*residual(2)*test_func*JxW[qp] * context.get_elem_solution_derivative();
236  (*Kww)(i,j) += -tau_M*residual(2)*d_test_func_dU(2)*u_phi[j][qp]*JxW[qp] * context.get_elem_solution_derivative();
237  (*Kww)(i,j) += -d_tau_M_dU(2)*u_phi[j][qp]*residual(2)*test_func*JxW[qp] * context.get_elem_solution_derivative();
238  }
239  }
240 
241  } // End compute_jacobian check
242 
243  } // End i dof loop
244  } // End quadrature loop
245  }
246 
247  template<class Mu>
248  void BoussinesqBuoyancyAdjointStabilization<Mu>::element_constraint
249  ( bool compute_jacobian, AssemblyContext & context )
250  {
251  // The number of local degrees of freedom in each variable.
252  const unsigned int n_p_dofs = context.get_dof_indices(_press_var.p()).size();
253  const unsigned int n_u_dofs = context.get_dof_indices(_flow_vars.u()).size();
254  const unsigned int n_T_dofs = context.get_dof_indices(_temp_vars.T()).size();
255 
256  // Element Jacobian * quadrature weights for interior integration.
257  const std::vector<libMesh::Real> &JxW =
258  context.get_element_fe(_flow_vars.u())->get_JxW();
259 
260  const std::vector<std::vector<libMesh::Real> >& T_phi =
261  context.get_element_fe(this->_temp_vars.T())->get_phi();
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> &KpT =
272  context.get_elem_jacobian(_press_var.p(), _temp_vars.T()); // J_{pT}
273  libMesh::DenseSubMatrix<libMesh::Number> &Kpu =
274  context.get_elem_jacobian(_press_var.p(), _flow_vars.u()); // J_{pu}
275  libMesh::DenseSubMatrix<libMesh::Number> &Kpv =
276  context.get_elem_jacobian(_press_var.p(), _flow_vars.v()); // J_{pv}
277  libMesh::DenseSubMatrix<libMesh::Number> *Kpw = NULL;
278 
279  if(this->_flow_vars.dim() == 3)
280  {
281  Kpw = &context.get_elem_jacobian
282  (_press_var.p(), _flow_vars.w()); // J_{pw}
283  }
284 
285  // Now we will build the element Jacobian and residual.
286  // Constructing the residual requires the solution and its
287  // gradient from the previous timestep. This must be
288  // calculated at each quadrature point by summing the
289  // solution degree-of-freedom values by the appropriate
290  // weight functions.
291  unsigned int n_qpoints = context.get_element_qrule().n_points();
292 
293  libMesh::FEBase* fe = context.get_element_fe(this->_flow_vars.u());
294 
295  for (unsigned int qp=0; qp != n_qpoints; qp++)
296  {
297  libMesh::RealGradient g = this->_stab_helper.compute_g( fe, context, qp );
298  libMesh::RealTensor G = this->_stab_helper.compute_G( fe, context, qp );
299 
300  libMesh::RealGradient U( context.interior_value( this->_flow_vars.u(), qp ),
301  context.interior_value( this->_flow_vars.v(), qp ) );
302  if( this->_flow_vars.dim() == 3 )
303  {
304  U(2) = context.interior_value( this->_flow_vars.w(), qp );
305  }
306 
307  // Compute the viscosity at this qp
308  libMesh::Real mu_qp = this->_mu(context, qp);
309 
310  libMesh::Real tau_M;
311  libMesh::Real d_tau_M_d_rho;
312  libMesh::Gradient d_tau_M_dU;
313 
314  if (compute_jacobian)
315  this->_stab_helper.compute_tau_momentum_and_derivs
316  ( context, qp, g, G, this->_rho, U, mu_qp,
317  tau_M, d_tau_M_d_rho, d_tau_M_dU,
318  this->_is_steady );
319  else
320  tau_M = this->_stab_helper.compute_tau_momentum
321  ( context, qp, g, G, this->_rho, U, mu_qp,
322  this->_is_steady );
323 
324  // Compute the solution & its gradient at the old Newton iterate.
325  libMesh::Number T;
326  T = context.interior_value(_temp_vars.T(), qp);
327 
328  libMesh::RealGradient d_residual_dT = _rho*_beta_T*_g;
329  // d_residual_dU = 0
330  libMesh::RealGradient residual = (T-_T_ref)*d_residual_dT;
331 
332  // First, an i-loop over the velocity degrees of freedom.
333  // We know that n_u_dofs == n_v_dofs so we can compute contributions
334  // for both at the same time.
335  for (unsigned int i=0; i != n_p_dofs; i++)
336  {
337  Fp(i) += tau_M*residual*p_dphi[i][qp]*JxW[qp];
338 
339  if (compute_jacobian)
340  {
341  for (unsigned int j=0; j != n_T_dofs; ++j)
342  {
343  KpT(i,j) += tau_M*d_residual_dT*T_phi[j][qp]*p_dphi[i][qp]*JxW[qp] * context.get_elem_solution_derivative();
344  }
345 
346  for (unsigned int j=0; j != n_u_dofs; ++j)
347  {
348  Kpu(i,j) += d_tau_M_dU(0)*u_phi[j][qp]*residual*p_dphi[i][qp]*JxW[qp] * context.get_elem_solution_derivative();
349  Kpv(i,j) += d_tau_M_dU(1)*u_phi[j][qp]*residual*p_dphi[i][qp]*JxW[qp] * context.get_elem_solution_derivative();
350  }
351  if( this->_flow_vars.dim() == 3 )
352  for (unsigned int j=0; j != n_u_dofs; ++j)
353  {
354  (*Kpw)(i,j) += d_tau_M_dU(2)*u_phi[j][qp]*residual*p_dphi[i][qp]*JxW[qp] * context.get_elem_solution_derivative();
355  }
356  }
357  }
358  } // End quadrature loop
359  }
360 
361  template<class Mu>
362  void BoussinesqBuoyancyAdjointStabilization<Mu>::register_parameter
363  ( const std::string & param_name,
364  libMesh::ParameterMultiAccessor<libMesh::Number> & param_pointer )
365  const
366  {
367  ParameterUser::register_parameter(param_name, param_pointer);
368  _mu.register_parameter(param_name, param_pointer);
369  }
370 
371 
372 } // namespace GRINS
373 
374 // Instantiate
375 INSTANTIATE_INC_NS_SUBCLASS(BoussinesqBuoyancyAdjointStabilization);
INSTANTIATE_INC_NS_SUBCLASS
INSTANTIATE_INC_NS_SUBCLASS(BoussinesqBuoyancyAdjointStabilization)
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
inc_nav_stokes_macro.h
GRINS::BoussinesqBuoyancyAdjointStabilization::element_time_derivative
virtual void element_time_derivative(bool compute_jacobian, AssemblyContext &context)
Time dependent part(s) of physics for element interiors.
Definition: boussinesq_buoyancy_adjoint_stab.C:69
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::element_constraint
virtual void element_constraint(bool compute_jacobian, AssemblyContext &context)
Constraint part(s) of physics for element interiors.
Definition: boussinesq_buoyancy_adjoint_stab.C:249
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:363
Generated on Tue Dec 19 2017 12:47:28 for GRINS-0.8.0 by  doxygen 1.8.9.1