www.mooseframework.org
INSChorinPredictor.C
Go to the documentation of this file.
1 //* This file is part of the MOOSE framework
2 //* https://www.mooseframework.org
3 //*
4 //* All rights reserved, see COPYRIGHT for full restrictions
5 //* https://github.com/idaholab/moose/blob/master/COPYRIGHT
6 //*
7 //* Licensed under LGPL 2.1, please see LICENSE for details
8 //* https://www.gnu.org/licenses/lgpl-2.1.html
9 
10 #include "INSChorinPredictor.h"
11 #include "MooseMesh.h"
12 
13 registerMooseObject("NavierStokesApp", INSChorinPredictor);
14 
17 {
19 
20  params.addClassDescription("This class computes the 'Chorin' Predictor equation in "
21  "fully-discrete (both time and space) form.");
22  // Coupled variables
23  params.addRequiredCoupledVar("u", "x-velocity");
24  params.addCoupledVar("v", "y-velocity"); // only required in 2D and 3D
25  params.addCoupledVar("w", "z-velocity"); // only required in 3D
26 
27  // Make star also be required, even though we might not use it?
28  params.addRequiredCoupledVar("u_star", "star x-velocity");
29  params.addCoupledVar("v_star", "star y-velocity"); // only required in 2D and 3D
30  params.addCoupledVar("w_star", "star z-velocity"); // only required in 3D
31 
32  // Required parameters
33  params.addRequiredRangeCheckedParam<unsigned>(
34  "component",
35  "component>=0 & component<=2",
36  "0,1,2 depending on if we are solving the x,y,z component of the Predictor equation");
37  MooseEnum predictor_type("OLD NEW STAR");
39  "predictor_type",
40  predictor_type,
41  "One of: OLD, NEW, STAR. Indicates which velocity to use in the predictor.");
42 
43  // Optional parameters
44  params.addParam<MaterialPropertyName>("mu_name", "mu", "The name of the dynamic viscosity");
45  params.addParam<MaterialPropertyName>("rho_name", "rho", "The name of the density");
46 
47  return params;
48 }
49 
51  : Kernel(parameters),
52 
53  // Current velocities
54  _u_vel(coupledValue("u")),
55  _v_vel(_mesh.dimension() >= 2 ? coupledValue("v") : _zero),
56  _w_vel(_mesh.dimension() == 3 ? coupledValue("w") : _zero),
57 
58  // Old velocities
59  _u_vel_old(coupledValueOld("u")),
60  _v_vel_old(_mesh.dimension() >= 2 ? coupledValueOld("v") : _zero),
61  _w_vel_old(_mesh.dimension() == 3 ? coupledValueOld("w") : _zero),
62 
63  // Star velocities
64  _u_vel_star(coupledValue("u_star")),
65  _v_vel_star(_mesh.dimension() >= 2 ? coupledValue("v_star") : _zero),
66  _w_vel_star(_mesh.dimension() == 3 ? coupledValue("w_star") : _zero),
67 
68  // Velocity Gradients
69  _grad_u_vel(coupledGradient("u")),
70  _grad_v_vel(_mesh.dimension() >= 2 ? coupledGradient("v") : _grad_zero),
71  _grad_w_vel(_mesh.dimension() == 3 ? coupledGradient("w") : _grad_zero),
72 
73  // Old Velocity Gradients
74  _grad_u_vel_old(coupledGradientOld("u")),
75  _grad_v_vel_old(_mesh.dimension() >= 2 ? coupledGradientOld("v") : _grad_zero),
76  _grad_w_vel_old(_mesh.dimension() == 3 ? coupledGradientOld("w") : _grad_zero),
77 
78  // Star Velocity Gradients
79  _grad_u_vel_star(coupledGradient("u_star")),
80  _grad_v_vel_star(_mesh.dimension() >= 2 ? coupledGradient("v_star") : _grad_zero),
81  _grad_w_vel_star(_mesh.dimension() == 3 ? coupledGradient("w_star") : _grad_zero),
82 
83  // Variable numberings
84  _u_vel_var_number(coupled("u")),
85  _v_vel_var_number(_mesh.dimension() >= 2 ? coupled("v") : libMesh::invalid_uint),
86  _w_vel_var_number(_mesh.dimension() == 3 ? coupled("w") : libMesh::invalid_uint),
87 
88  // Star velocity numberings
89  _u_vel_star_var_number(coupled("u_star")),
90  _v_vel_star_var_number(_mesh.dimension() >= 2 ? coupled("v_star") : libMesh::invalid_uint),
91  _w_vel_star_var_number(_mesh.dimension() == 3 ? coupled("w_star") : libMesh::invalid_uint),
92 
93  // Required parameters
94  _component(getParam<unsigned>("component")),
95  _predictor_enum(getParam<MooseEnum>("predictor_type")),
96 
97  // Material properties
98  _mu(getMaterialProperty<Real>("mu_name")),
99  _rho(getMaterialProperty<Real>("rho_name"))
100 {
101 }
102 
103 Real
105 {
106  // Vector object for test function
107  RealVectorValue test;
108  test(_component) = _test[_i][_qp];
109 
110  // Tensor object for test function gradient
111  RealTensorValue grad_test;
112  for (unsigned k = 0; k < 3; ++k)
113  grad_test(_component, k) = _grad_test[_i][_qp](k);
114 
115  // Decide what velocity vector, gradient to use:
116  RealVectorValue U;
117  RealTensorValue grad_U;
118 
119  switch (_predictor_enum)
120  {
121  case OLD:
122  {
125  break;
126  }
127  case NEW:
128  {
131  break;
132  }
133  case STAR:
134  {
135  // Note: Donea and Huerta's book says you are supposed to use "star" velocity to make Chorin
136  // implicit, not U^{n+1}.
139  break;
140  }
141  default:
142  mooseError("Unrecognized Chorin predictor type requested.");
143  }
144 
145  //
146  // Compute the different parts
147  //
148 
149  // Note: _u is the component'th entry of "u_star" in Chorin's method.
151  Real symmetric_part = (_u[_qp] - U_old(_component)) * _test[_i][_qp];
152 
153  // Convective part. Remember to multiply by _dt!
154  Real convective_part = _dt * (grad_U * U) * test;
155 
156  // Viscous part - we are using the Laplacian form here for simplicity.
157  // Remember to multiply by _dt!
158  Real viscous_part = _dt * (_mu[_qp] / _rho[_qp]) * grad_U.contract(grad_test);
159 
160  return symmetric_part + convective_part + viscous_part;
161 }
162 
163 Real
165 {
166  // The mass matrix part is always there.
167  Real mass_part = _phi[_j][_qp] * _test[_i][_qp];
168 
169  // The on-diagonal Jacobian contribution depends on whether the predictor uses the
170  // 'new' or 'star' velocity.
171  Real other_part = 0.;
172  switch (_predictor_enum)
173  {
174  case OLD:
175  case NEW:
176  break;
177  case STAR:
178  {
180  Real convective_part =
181  _dt * ((U_star * _grad_phi[_j][_qp]) + _phi[_j][_qp] * _grad_u[_qp](_component)) *
182  _test[_i][_qp];
183  Real viscous_part =
184  _dt * ((_mu[_qp] / _rho[_qp]) * (_grad_phi[_j][_qp] * _grad_test[_i][_qp]));
185  other_part = convective_part + viscous_part;
186  break;
187  }
188  default:
189  mooseError("Unrecognized Chorin predictor type requested.");
190  }
191 
192  return mass_part + other_part;
193 }
194 
195 Real
197 {
198  switch (_predictor_enum)
199  {
200  case OLD:
201  {
202  return 0.;
203  }
204 
205  case NEW:
206  {
207  if ((jvar == _u_vel_var_number) || (jvar == _v_vel_var_number) || (jvar == _w_vel_var_number))
208  {
209  // Derivative of grad_U wrt the velocity component
210  RealTensorValue dgrad_U;
211 
212  // Initialize to invalid value, then determine correct value.
213  unsigned vel_index = 99;
214 
215  // Map jvar into the indices (0,1,2)
216  if (jvar == _u_vel_var_number)
217  vel_index = 0;
218 
219  else if (jvar == _v_vel_var_number)
220  vel_index = 1;
221 
222  else if (jvar == _w_vel_var_number)
223  vel_index = 2;
224 
225  // Fill in the vel_index'th row of dgrad_U with _grad_phi[_j][_qp]
226  for (unsigned k = 0; k < 3; ++k)
227  dgrad_U(vel_index, k) = _grad_phi[_j][_qp](k);
228 
229  // Vector object for test function
230  RealVectorValue test;
231  test(_component) = _test[_i][_qp];
232 
233  // Vector object for U
235 
236  // Tensor object for test function gradient
237  RealTensorValue grad_test;
238  for (unsigned k = 0; k < 3; ++k)
239  grad_test(_component, k) = _grad_test[_i][_qp](k);
240 
241  // Compute the convective part
242  RealVectorValue convective_jac =
243  _phi[_j][_qp] * RealVectorValue(_grad_u_vel[_qp](vel_index),
244  _grad_v_vel[_qp](vel_index),
245  _grad_w_vel[_qp](vel_index));
246 
247  // Extra contribution in vel_index component
248  convective_jac(vel_index) += U * _grad_phi[_j][_qp];
249 
250  // Be sure to scale by _dt!
251  Real convective_part = _dt * (convective_jac * test);
252 
253  // Compute the viscous part, be sure to scale by _dt. Note: the contracted
254  // value should be zero unless vel_index and _component match.
255  Real viscous_part = _dt * (_mu[_qp] / _rho[_qp]) * dgrad_U.contract(grad_test);
256 
257  // Return the result
258  return convective_part + viscous_part;
259  }
260  else
261  return 0;
262  }
263 
264  case STAR:
265  {
266  if (jvar == _u_vel_star_var_number)
267  {
268  return _dt * _phi[_j][_qp] * _grad_u[_qp](0) * _test[_i][_qp];
269  }
270 
271  else if (jvar == _v_vel_star_var_number)
272  {
273  return _dt * _phi[_j][_qp] * _grad_u[_qp](1) * _test[_i][_qp];
274  }
275 
276  else if (jvar == _w_vel_star_var_number)
277  {
278  return _dt * _phi[_j][_qp] * _grad_u[_qp](2) * _test[_i][_qp];
279  }
280 
281  else
282  return 0;
283  }
284 
285  default:
286  mooseError("Unrecognized Chorin predictor type requested.");
287  }
288 }
const MaterialProperty< Real > & _rho
const VariableValue & _w_vel
virtual Real computeQpJacobian()
const MaterialProperty< Real > & _mu
void addRequiredRangeCheckedParam(const std::string &name, const std::string &parsed_function, const std::string &doc_string)
const VariableGradient & _grad_u
static InputParameters validParams()
const unsigned int invalid_uint
This class computes the "Chorin" Predictor equation in fully-discrete (both time and space) form...
void addParam(const std::string &name, const std::initializer_list< typename T::value_type > &value, const std::string &doc_string)
const VariableValue & _u_vel
const VariablePhiGradient & _grad_phi
const VariableGradient & _grad_u_vel
virtual Real computeQpResidual()
const VariableGradient & _grad_v_vel_old
Real & _dt
The following methods are specializations for using the Parallel::packed_range_* routines for a vecto...
const VariableValue & _v_vel_old
void addRequiredParam(const std::string &name, const std::string &doc_string)
const VariableValue & _w_vel_old
const VariableValue & _u_vel_old
const VariableValue & _v_vel
TensorValue< Real > RealTensorValue
const VariableTestValue & _test
registerMooseObject("NavierStokesApp", INSChorinPredictor)
const VariableGradient & _grad_u_vel_old
const VariableGradient & _grad_w_vel
const VariableValue & _v_vel_star
const VariableGradient & _grad_v_vel_star
unsigned int _i
virtual Real computeQpOffDiagJacobian(unsigned jvar)
INSChorinPredictor(const InputParameters &parameters)
void addCoupledVar(const std::string &name, const std::string &doc_string)
void addRequiredCoupledVar(const std::string &name, const std::string &doc_string)
const VariableValue & _w_vel_star
const VariableGradient & _grad_w_vel_old
unsigned int _j
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
static InputParameters validParams()
const VariableTestGradient & _grad_test
void mooseError(Args &&... args) const
void addClassDescription(const std::string &doc_string)
const VariableGradient & _grad_w_vel_star
const VariableGradient & _grad_u_vel_star
const VariablePhiValue & _phi
const VariableGradient & _grad_v_vel
const VariableValue & _u_vel_star
static const std::string k
Definition: NS.h:124
const VariableValue & _u
unsigned int _qp