Line data    Source code

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