www.mooseframework.org
AEFVSlopeLimitingOneD.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 "AEFVSlopeLimitingOneD.h"
11 
13 
14 template <>
15 InputParameters
17 {
18  InputParameters params = validParams<SlopeLimitingBase>();
19  params.addClassDescription("One-dimensional slope limiting to get the limited slope of cell "
20  "average variable for the advection equation using a cell-centered "
21  "finite volume method.");
22  params.addRequiredCoupledVar("u", "constant monomial variable");
23  MooseEnum scheme("none minmod mc superbee", "none");
24  params.addParam<MooseEnum>("scheme", scheme, "TVD-type slope limiting scheme");
25  return params;
26 }
27 
28 AEFVSlopeLimitingOneD::AEFVSlopeLimitingOneD(const InputParameters & parameters)
29  : SlopeLimitingBase(parameters), _u(getVar("u", 0)), _scheme(getParam<MooseEnum>("scheme"))
30 {
31 }
32 
33 std::vector<RealGradient>
35 {
36  // you should know how many equations you are solving and assign this number
37  // e.g. = 1 (for the advection equation)
38  unsigned int nvars = 1;
39 
40  const Elem * elem = _current_elem;
41 
42  // index of conserved variable
43  unsigned int iv;
44 
45  // index of sides around an element
46  unsigned int is;
47 
48  // index of the neighbor element
49  unsigned int in;
50 
51  // number of sides surrounding an element = 2 in 1D
52  unsigned int nside = elem->n_sides();
53 
54  // number of reconstruction stencils = 3 in 1D
55  unsigned int nsten = nside + 1;
56 
57  // vector for the gradients of primitive variables
58  std::vector<RealGradient> ugrad(nvars, RealGradient(0., 0., 0.));
59 
60  // array to store center coordinates of this cell and its neighbor cells
61  std::vector<Real> xc(nsten, 0.);
62 
63  // the first always stores the current cell
64  xc[0] = elem->centroid()(0);
65 
66  // array for the cell-average values in the current cell and its neighbors
67  std::vector<std::vector<Real>> ucell(nsten, std::vector<Real>(nvars, 0.));
68 
69  // central slope:
70  // u_{i+1} - u {i-1}
71  // -----------------
72  // x_{i+1} - x_{i-1}
73 
74  // left-side slope:
75  // u_{i} - u {i-1}
76  // ---------------
77  // x_{i} - x_{i-1}
78 
79  // right-side slope:
80  // u_{i+1} - u {i}
81  // ---------------
82  // x_{i+1} - x_{i}
83 
84  // array to store the central and one-sided slopes, where the first should be central slope
85  std::vector<std::vector<Real>> sigma(nsten, std::vector<Real>(nvars, 0.));
86 
87  // get the cell-average variable in the central cell
88  if (_is_implicit)
89  ucell[0][0] = _u->getElementalValue(elem);
90  else
91  ucell[0][0] = _u->getElementalValueOld(elem);
92 
93  // a flag to indicate the boundary side of the current cell
94 
95  unsigned int bflag = 0;
96 
97  // loop over the sides to compute the one-sided slopes
98 
99  for (is = 0; is < nside; is++)
100  {
101  in = is + 1;
102 
103  // when the current element is an internal cell
104  if (elem->neighbor_ptr(is) != NULL)
105  {
106  const Elem * neig = elem->neighbor_ptr(is);
107  if (this->hasBlocks(neig->subdomain_id()))
108  {
109  xc[in] = neig->centroid()(0);
110 
111  // get the cell-average variable in this neighbor cell
112  if (_is_implicit)
113  ucell[in][0] = _u->getElementalValue(neig);
114  else
115  ucell[in][0] = _u->getElementalValueOld(neig);
116 
117  // calculate the one-sided slopes of primitive variables
118 
119  for (iv = 0; iv < nvars; iv++)
120  sigma[in][iv] = (ucell[0][iv] - ucell[in][iv]) / (xc[0] - xc[in]);
121  }
122  else
123  bflag = in;
124  }
125  // when the current element is at the boundary,
126  // we choose not to construct the slope in 1D just for convenience.
127  else
128  {
129  bflag = in;
130  }
131  }
132 
133  // calculate the slope of the current element
134  // according to the choice of the slope limiter algorithm
135 
136  switch (_scheme)
137  {
138  // first-order, zero slope
139  case 0:
140  break;
141 
142  // ==============
143  // minmod limiter
144  // ==============
145  case 1:
146 
147  if (bflag == 0)
148  {
149  for (iv = 0; iv < nvars; iv++)
150  {
151  if ((sigma[1][iv] * sigma[2][iv]) > 0.)
152  {
153  if (std::abs(sigma[1][iv]) < std::abs(sigma[2][iv]))
154  ugrad[iv](0) = sigma[1][iv];
155  else
156  ugrad[iv](0) = sigma[2][iv];
157  }
158  }
159  }
160  break;
161 
162  // ===========================================
163  // MC (monotonized central-difference limiter)
164  // ===========================================
165  case 2:
166 
167  if (bflag == 0)
168  {
169  for (iv = 0; iv < nvars; iv++)
170  sigma[0][iv] = (ucell[1][iv] - ucell[2][iv]) / (xc[1] - xc[2]);
171 
172  for (iv = 0; iv < nvars; iv++)
173  {
174  if (sigma[0][iv] > 0. && sigma[1][iv] > 0. && sigma[2][iv] > 0.)
175  ugrad[iv](0) = std::min(sigma[0][iv], 2. * std::min(sigma[1][iv], sigma[2][iv]));
176  else if (sigma[0][iv] < 0. && sigma[1][iv] < 0. && sigma[2][iv] < 0.)
177  ugrad[iv](0) = std::max(sigma[0][iv], 2. * std::max(sigma[1][iv], sigma[2][iv]));
178  }
179  }
180  break;
181 
182  // ================
183  // Superbee limiter
184  // ================
185  case 3:
186 
187  if (bflag == 0)
188  {
189  for (iv = 0; iv < nvars; iv++)
190  {
191  Real sigma1 = 0.;
192  Real sigma2 = 0.;
193 
194  // calculate sigma1 with minmod
195  if (sigma[2][iv] > 0. && sigma[1][iv] > 0.)
196  sigma1 = std::min(sigma[2][iv], 2. * sigma[1][iv]);
197  else if (sigma[2][iv] < 0. && sigma[1][iv] < 0.)
198  sigma1 = std::max(sigma[2][iv], 2. * sigma[1][iv]);
199 
200  // calculate sigma2 with minmod
201  if (sigma[2][iv] > 0. && sigma[1][iv] > 0.)
202  sigma2 = std::min(2. * sigma[2][iv], sigma[1][iv]);
203  else if (sigma[2][iv] < 0. && sigma[1][iv] < 0.)
204  sigma2 = std::max(2. * sigma[2][iv], sigma[1][iv]);
205 
206  // calculate sigma with maxmod
207  if (sigma1 > 0. && sigma2 > 0.)
208  ugrad[iv](0) = std::max(sigma1, sigma2);
209  else if (sigma1 < 0. && sigma2 < 0.)
210  ugrad[iv](0) = std::min(sigma1, sigma2);
211  }
212  }
213  break;
214 
215  default:
216  mooseError("Unknown 1D TVD-type slope limiter scheme");
217  break;
218  }
219 
220  return ugrad;
221 }
AEFVSlopeLimitingOneD::_u
MooseVariable * _u
the input variable
Definition: AEFVSlopeLimitingOneD.h:36
SlopeLimitingBase
Base class for slope limiting to limit the slopes of cell average variables.
Definition: SlopeLimitingBase.h:24
AEFVSlopeLimitingOneD.h
AEFVSlopeLimitingOneD::_scheme
MooseEnum _scheme
One-D slope limiting scheme.
Definition: AEFVSlopeLimitingOneD.h:39
validParams< SlopeLimitingBase >
InputParameters validParams< SlopeLimitingBase >()
Definition: SlopeLimitingBase.C:20
libMesh::RealGradient
VectorValue< Real > RealGradient
Definition: GrainForceAndTorqueInterface.h:17
validParams< AEFVSlopeLimitingOneD >
InputParameters validParams< AEFVSlopeLimitingOneD >()
Definition: AEFVSlopeLimitingOneD.C:16
AEFVSlopeLimitingOneD::limitElementSlope
virtual std::vector< RealGradient > limitElementSlope() const override
compute the limited slope of the cell
Definition: AEFVSlopeLimitingOneD.C:34
AEFVSlopeLimitingOneD
One-dimensional slope limiting to get the limited slope of cell average variable for the advection eq...
Definition: AEFVSlopeLimitingOneD.h:26
ElementLoopUserObject::_current_elem
const Elem * _current_elem
Definition: ElementLoopUserObject.h:91
registerMooseObject
registerMooseObject("RdgApp", AEFVSlopeLimitingOneD)
AEFVSlopeLimitingOneD::AEFVSlopeLimitingOneD
AEFVSlopeLimitingOneD(const InputParameters &parameters)
Definition: AEFVSlopeLimitingOneD.C:28