www.mooseframework.org
CNSFVLeastSquaresSlopeReconstruction.C
Go to the documentation of this file.
1 /****************************************************************/
2 /* MOOSE - Multiphysics Object Oriented Simulation Environment */
3 /* */
4 /* All contents are licensed under LGPL V2.1 */
5 /* See LICENSE for full restrictions */
6 /****************************************************************/
7 
9 
10 template <>
11 InputParameters
13 {
14  InputParameters params = validParams<SlopeReconstructionMultiD>();
15 
16  params.addClassDescription("A user object that performs the least-squares slope reconstruction "
17  "to get the slopes of the P0 primitive variables.");
18 
19  params.addRequiredCoupledVar("rho", "Density at P0 (constant monomial)");
20 
21  params.addRequiredCoupledVar("rhou", "X-momentum at P0 (constant monomial)");
22 
23  params.addCoupledVar("rhov", "Y-momentum at P0 (constant monomial)");
24 
25  params.addCoupledVar("rhow", "Z-momentum at P0 (constant monomial)");
26 
27  params.addRequiredCoupledVar("rhoe", "Total energy at P0 (constant monomial)");
28 
29  params.addRequiredParam<UserObjectName>("fluid_properties",
30  "Name for fluid properties user object");
31 
32  return params;
33 }
34 
36  const InputParameters & parameters)
37  : SlopeReconstructionMultiD(parameters),
38  _rho(getVar("rho", 0)),
39  _rhou(getVar("rhou", 0)),
40  _rhov(isCoupled("rhov") ? getVar("rhov", 0) : nullptr),
41  _rhow(isCoupled("rhow") ? getVar("rhow", 0) : nullptr),
42  _rhoe(getVar("rhoe", 0)),
43  _fp(getUserObject<SinglePhaseFluidProperties>("fluid_properties"))
44 {
45 }
46 
47 void
49 {
50  const Elem * elem = _current_elem;
51 
53  dof_id_type elemID = elem->id();
54 
56  unsigned int nside = elem->n_sides();
57 
59  unsigned int nvars = 5;
60 
62  std::vector<RealGradient> ugrad(nvars, RealGradient(0., 0., 0.));
63 
65  Point snorm = Point(0., 0., 0.);
66 
68  Point scent = Point(0., 0., 0.);
69 
71  Point dcent = Point(0., 0., 0.);
72 
74  std::vector<std::vector<Real>> u(nside + 1, std::vector<Real>(nvars, 0.));
75 
78 
79  Real rhoc = _rho->getElementalValue(elem);
80  Real rhomc = 1. / rhoc;
81  Real rhouc = _rhou->getElementalValue(elem);
82  Real rhovc = 0.;
83  Real rhowc = 0.;
84  Real rhoec = _rhoe->getElementalValue(elem);
85  if (_rhov != nullptr)
86  rhovc = _rhov->getElementalValue(elem);
87  if (_rhow != nullptr)
88  rhowc = _rhow->getElementalValue(elem);
89 
90  u[0][1] = rhouc * rhomc;
91  u[0][2] = rhovc * rhomc;
92  u[0][3] = rhowc * rhomc;
93  Real eintc = rhoec * rhomc - 0.5 * (u[0][1] * u[0][1] + u[0][2] * u[0][2] + u[0][3] * u[0][3]);
94  u[0][0] = _fp.pressure(rhomc, eintc);
95  u[0][4] = _fp.temperature(rhomc, eintc);
96 
98  _avars[elemID] = u[0];
99 
101  Real A11 = 0., A12 = 0., A13 = 0., A22 = 0., A23 = 0., A33 = 0.;
102 
104  Real C11 = 0., C12 = 0., C13 = 0., C21 = 0., C22 = 0., C23 = 0., C31 = 0., C32 = 0., C33 = 0.;
105 
107  std::vector<Real> B1(nvars, 0.);
108  std::vector<Real> B2(nvars, 0.);
109  std::vector<Real> B3(nvars, 0.);
110 
112  Real lspow = 2., lswei = 0.;
113 
115  Real invdet = 0., dx = 0., dy = 0., dz = 0.;
116 
118  std::vector<Real> du(nvars, 0.);
119 
121  bool bndElem = false;
122 
124 
125  for (unsigned int is = 0; is < nside; is++)
126  {
127  unsigned int in = is + 1;
128  const Elem * neig = elem->neighbor_ptr(is);
129 
131 
132  if (neig != nullptr && this->hasBlocks(neig->subdomain_id()))
133  {
134  dof_id_type neigID = neig->id();
135 
137  {
138  _assembly.reinit(elem, is);
139  _side_centroid[std::pair<dof_id_type, dof_id_type>(elemID, neigID)] =
140  _side_elem->centroid();
141  }
142 
145 
146  Real v = 1. / _rho->getElementalValue(neig);
147 
148  u[in][1] = _rhou->getElementalValue(neig) * v;
149  if (_rhov != nullptr)
150  u[in][2] = _rhov->getElementalValue(neig) * v;
151  if (_rhow != nullptr)
152  u[in][3] = _rhow->getElementalValue(neig) * v;
153 
154  Real e = _rhoe->getElementalValue(neig) * v -
155  0.5 * (u[in][1] * u[in][1] + u[in][2] * u[in][2] + u[in][3] * u[in][3]);
156 
157  u[in][0] = _fp.pressure(v, e);
158  u[in][4] = _fp.temperature(v, e);
159 
161  _avars[neigID] = u[in];
162 
164 
165  dcent = neig->centroid() - elem->centroid();
166 
167  lswei = std::sqrt(dcent * dcent);
168 
169  lswei = std::pow(lswei, -lspow);
170 
171  for (unsigned int iv = 0; iv < nvars; iv++)
172  du[iv] = u[in][iv] - u[0][iv];
173 
174  if (_mesh.dimension() > 0)
175  {
176  dx = dcent(0);
177  A11 += lswei * dx * dx;
178  for (unsigned int iv = 0; iv < nvars; iv++)
179  B1[iv] += lswei * du[iv] * dx;
180  }
181  if (_mesh.dimension() > 1)
182  {
183  dy = dcent(1);
184  A12 += lswei * dx * dy;
185  A22 += lswei * dy * dy;
186  for (unsigned int iv = 0; iv < nvars; iv++)
187  B2[iv] += lswei * du[iv] * dy;
188  }
189  if (_mesh.dimension() > 2)
190  {
191  dz = dcent(2);
192  A13 += lswei * dx * dz;
193  A23 += lswei * dy * dz;
194  A33 += lswei * dz * dz;
195  for (unsigned int iv = 0; iv < nvars; iv++)
196  B3[iv] += lswei * du[iv] * dz;
197  }
198  }
199 
201 
202  else
203  {
204  bndElem = true;
205 
207  {
208  scent = getBoundarySideCentroid(elemID, is);
209  snorm = getBoundarySideNormal(elemID, is);
210  }
211  else
212  {
213  _assembly.reinit(elem, is);
214  scent = _side_elem->centroid();
215  snorm = _normals_face[0];
216  _bnd_side_centroid[std::pair<dof_id_type, unsigned int>(elemID, is)] = scent;
217  _bnd_side_normal[std::pair<dof_id_type, unsigned int>(elemID, is)] = snorm;
218  }
219 
221 
222  std::vector<BoundaryID> bndID = _mesh.getBoundaryIDs(elem, is);
223 
224  for (unsigned int ib = 0; ib < bndID.size(); ib++)
225  {
226  BoundaryID currentID = bndID[ib];
227 
228  std::map<BoundaryID, UserObjectName>::const_iterator pos = _bnd_uo_name_map.find(currentID);
229 
230  if (pos != _bnd_uo_name_map.end())
231  {
232  const BCUserObject & bcuo =
233  GeneralUserObject::_fe_problem.getUserObject<BCUserObject>(pos->second);
234 
235  // another question is if we can use (rho,u,v,w,p) as primitive variabls
236  // instead of the current (p,u,v,w,T)
237  // TODO test it in 1D problems first
238 
239  std::vector<Real> uvec1 = {rhoc, rhouc, rhovc, rhowc, rhoec};
240 
241  u[in] = bcuo.getGhostCellValue(is, elem->id(), uvec1, snorm);
242 
244 
245  Real rhog = u[in][0];
246  Real rhoug = u[in][1];
247  Real rhovg = u[in][2];
248  Real rhowg = u[in][3];
249  Real rhoeg = u[in][4];
250 
251  Real rhomg = 1. / rhog;
252 
253  u[in][1] = rhoug * rhomg;
254  u[in][2] = rhovg * rhomg;
255  u[in][3] = rhowg * rhomg;
256 
257  Real eintg = rhoeg * rhomg -
258  0.5 * (u[in][1] * u[in][1] + u[in][2] * u[in][2] + u[in][3] * u[in][3]);
259 
260  u[in][0] = _fp.pressure(rhomg, eintg);
261  u[in][4] = _fp.temperature(rhomg, eintg);
262  }
263  }
264 
266  _bnd_avars[std::pair<dof_id_type, unsigned int>(elemID, is)] = u[in];
267 
269 
270  dcent = 2. * (scent - elem->centroid());
271 
272  lswei = std::sqrt(dcent * dcent);
273 
274  lswei = std::pow(lswei, -lspow);
275 
276  for (unsigned int iv = 0; iv < nvars; iv++)
277  du[iv] = u[in][iv] - u[0][iv];
278 
279  if (_mesh.dimension() > 0)
280  {
281  dx = dcent(0);
282  A11 += lswei * dx * dx;
283  for (unsigned int iv = 0; iv < nvars; iv++)
284  B1[iv] += lswei * du[iv] * dx;
285  }
286  if (_mesh.dimension() > 1)
287  {
288  dy = dcent(1);
289  A12 += lswei * dx * dy;
290  A22 += lswei * dy * dy;
291  for (unsigned int iv = 0; iv < nvars; iv++)
292  B2[iv] += lswei * du[iv] * dy;
293  }
294  if (_mesh.dimension() > 2)
295  {
296  dz = dcent(2);
297  A13 += lswei * dx * dz;
298  A23 += lswei * dy * dz;
299  A33 += lswei * dz * dz;
300  for (unsigned int iv = 0; iv < nvars; iv++)
301  B3[iv] += lswei * du[iv] * dz;
302  }
303  }
304  }
305 
307 
308  if (_mesh.dimension() == 1)
309  {
310  for (unsigned int iv = 0; iv < nvars; iv++)
311  ugrad[iv](0) = B1[iv] / A11;
312  }
313  else if (_mesh.dimension() == 2)
314  {
315  invdet = 1. / (A11 * A22 - A12 * A12);
316 
317  for (unsigned int iv = 0; iv < nvars; iv++)
318  {
319  ugrad[iv](0) = invdet * (A22 * B1[iv] - A12 * B2[iv]);
320  ugrad[iv](1) = invdet * (-A12 * B1[iv] + A11 * B2[iv]);
321  }
322  }
323  else if (_mesh.dimension() == 3)
324  {
325  C11 = A22 * A33 - A23 * A23;
326  C12 = A13 * A23 - A12 * A33;
327  C13 = A12 * A23 - A13 * A22;
328 
329  C21 = A13 * A23 - A12 * A33;
330  C22 = A11 * A33 - A13 * A13;
331  C23 = A12 * A13 - A11 * A23;
332 
333  C31 = A12 * A23 - A22 * A13;
334  C32 = A12 * A13 - A11 * A23;
335  C33 = A11 * A22 - A12 * A12;
336 
337  invdet = 1. / (A11 * C11 + A12 * C12 + A13 * C12);
338 
339  for (unsigned int iv = 0; iv < nvars; iv++)
340  {
341  ugrad[iv](0) = invdet * (C11 * B1[iv] + C12 * B2[iv] + C13 * B3[iv]);
342  ugrad[iv](1) = invdet * (C21 * B1[iv] + C22 * B2[iv] + C23 * B3[iv]);
343  ugrad[iv](2) = invdet * (C31 * B1[iv] + C32 * B2[iv] + C33 * B3[iv]);
344  }
345  }
346 
347  _rslope[elemID] = ugrad;
348 }
InputParameters validParams< CNSFVLeastSquaresSlopeReconstruction >()
virtual Real temperature(Real v, Real u) const =0
Temperature as a function of specific internal energy and specific volume.
virtual void reconstructElementSlope()
compute the slope of the cell
std::map< dof_id_type, std::vector< Real > > _avars
store the average variable values into this map indexed by element ID
MooseVariable * _rhoe
the input total energy
const SinglePhaseFluidProperties & _fp
fluid properties user object
std::map< BoundaryID, UserObjectName > _bnd_uo_name_map
store the pair of boundary ID and user object name
std::map< dof_id_type, std::vector< RealGradient > > _rslope
store the reconstructed slopes into this map indexed by element ID
CNSFVLeastSquaresSlopeReconstruction(const InputParameters &parameters)
Common class for single phase fluid properties.
InputParameters validParams< SlopeReconstructionMultiD >()
std::map< std::pair< dof_id_type, unsigned int >, Point > _bnd_side_normal
store the boundary side normal into this map indexed by pair of element ID and local side ID ...
virtual const Point & getBoundarySideNormal(dof_id_type elementid, unsigned int side) const
accessor function call to get cached boundary side centroid
ExpressionBuilder::EBTerm pow(const ExpressionBuilder::EBTerm &left, T exponent)
const MooseArray< Point > & _normals_face
virtual std::vector< Real > getGhostCellValue(unsigned int iside, dof_id_type ielem, const std::vector< Real > &uvec1, const RealVectorValue &dwave) const =0
compute the ghost cell variable values
A base class of user object for calculating the variable values in ghost element according to specifi...
Definition: BCUserObject.h:42
std::map< std::pair< dof_id_type, dof_id_type >, Point > _side_centroid
store the side centroid into this map indexed by pair of element ID and neighbor ID ...
std::map< std::pair< dof_id_type, unsigned int >, Point > _bnd_side_centroid
store the boundary side centroid into this map indexed by pair of element ID and local side ID ...
virtual const Point & getBoundarySideCentroid(dof_id_type elementid, unsigned int side) const
accessor function call to get cached boundary side centroid
virtual Real pressure(Real v, Real u) const =0
Pressure as a function of specific internal energy and specific volume.
bool _side_geoinfo_cached
flag to indicated if side geometry info is cached
std::map< std::pair< dof_id_type, unsigned int >, std::vector< Real > > _bnd_avars
store the boundary average variable values into this map indexed by pair of element ID and local side...
Multi-dimensional piecewise linear slope reconstruction to get the slopes of cell average variables...