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