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 "RichardsSUPGstandard.h"
      11             : #include "libmesh/utility.h"
      12             : 
      13             : registerMooseObject("RichardsApp", RichardsSUPGstandard);
      14             : 
      15             : InputParameters
      16         578 : RichardsSUPGstandard::validParams()
      17             : {
      18         578 :   InputParameters params = RichardsSUPG::validParams();
      19        1156 :   params.addRequiredRangeCheckedParam<Real>(
      20             :       "p_SUPG",
      21             :       "p_SUPG > 0",
      22             :       "SUPG pressure.  This parameter controls the strength of the "
      23             :       "upwinding.  This parameter must be positive.  If you need to track "
      24             :       "advancing fronts in a problem, then set to less than your expected "
      25             :       "range of pressures in your unsaturated zone.  Otherwise, set "
      26             :       "larger, and then minimal upwinding will occur and convergence will "
      27             :       "typically be good.  If you need no SUPG, it is more efficient not "
      28             :       "to use this UserObject.");
      29         578 :   params.addClassDescription(
      30             :       "Standard SUPG relationships for Richards flow based on Appendix A of      TJR Hughes, M "
      31             :       "Mallet and A Mizukami ``A new finite element formulation for computational fluid dynamics:: "
      32             :       "II. Behond SUPG'' Computer Methods in Applied Mechanics and Engineering 54 (1986) 341--355");
      33         578 :   return params;
      34           0 : }
      35             : 
      36         285 : RichardsSUPGstandard::RichardsSUPGstandard(const InputParameters & parameters)
      37         570 :   : RichardsSUPG(parameters), _p_SUPG(getParam<Real>("p_SUPG"))
      38             : {
      39         285 : }
      40             : 
      41             : RealVectorValue
      42     6278018 : RichardsSUPGstandard::velSUPG(RealTensorValue perm,
      43             :                               RealVectorValue gradp,
      44             :                               Real density,
      45             :                               RealVectorValue gravity) const
      46             : {
      47     6278018 :   return -perm * (gradp - density * gravity); // points in direction of info propagation
      48             : }
      49             : 
      50             : RealTensorValue
      51     6278018 : RichardsSUPGstandard::dvelSUPG_dgradp(RealTensorValue perm) const
      52             : {
      53     6278018 :   return -perm;
      54             : }
      55             : 
      56             : RealVectorValue
      57     6278018 : RichardsSUPGstandard::dvelSUPG_dp(RealTensorValue perm,
      58             :                                   Real density_prime,
      59             :                                   RealVectorValue gravity) const
      60             : {
      61     6278018 :   return perm * density_prime * gravity;
      62             : }
      63             : 
      64             : Real
      65    18700854 : RichardsSUPGstandard::cosh_relation(Real alpha) const
      66             : {
      67    18700854 :   if (alpha >= 5.0 || alpha <= -5.0)
      68    12450507 :     return ((alpha > 0.0) ? 1.0 : -1.0) - 1.0 / alpha; // prevents overflows
      69     6250347 :   else if (alpha == 0)
      70             :     return 0.0;
      71     6250347 :   return 1.0 / std::tanh(alpha) - 1.0 / alpha;
      72             : }
      73             : 
      74             : Real
      75    12467236 : RichardsSUPGstandard::cosh_relation_prime(Real alpha) const
      76             : {
      77    12467236 :   if (alpha >= 5.0 || alpha <= -5.0)
      78     8300338 :     return 1.0 / (alpha * alpha); // prevents overflows
      79     4166898 :   else if (alpha == 0)
      80             :     return 1.0 / 3.0;
      81     4166898 :   return 1.0 - Utility::pow<2>(std::cosh(alpha) / std::sinh(alpha)) + 1.0 / (alpha * alpha);
      82             : }
      83             : 
      84             : // the following is physically 2*velocity/element_length
      85             : RealVectorValue
      86     6278018 : RichardsSUPGstandard::bb(RealVectorValue vel,
      87             :                          int dimen,
      88             :                          RealVectorValue xi_prime,
      89             :                          RealVectorValue eta_prime,
      90             :                          RealVectorValue zeta_prime) const
      91             : {
      92             :   RealVectorValue b;
      93     6278018 :   b(0) = xi_prime * vel;
      94     6278018 :   if (dimen >= 2)
      95     3612740 :     b(1) = eta_prime * vel;
      96     3612740 :   if (dimen == 3)
      97     1748844 :     b(2) = zeta_prime * vel;
      98     6278018 :   return b;
      99             : }
     100             : 
     101             : // following is d(bb*bb)/d(gradp)
     102             : RealVectorValue
     103     6278018 : RichardsSUPGstandard::dbb2_dgradp(RealVectorValue vel,
     104             :                                   RealTensorValue dvel_dgradp,
     105             :                                   RealVectorValue xi_prime,
     106             :                                   RealVectorValue eta_prime,
     107             :                                   RealVectorValue zeta_prime) const
     108             : {
     109             :   // if dvel_dgradp is symmetric so transpose is probably a waste of time
     110     6278018 :   return 2.0 * ((xi_prime * vel) * (dvel_dgradp.transpose() * xi_prime) +
     111     6278018 :                 (eta_prime * vel) * (dvel_dgradp.transpose() * eta_prime) +
     112     6278018 :                 (zeta_prime * vel) * (dvel_dgradp.transpose() * zeta_prime));
     113             : }
     114             : 
     115             : // following is d(bb*bb)/d(p)
     116             : Real
     117     6278018 : RichardsSUPGstandard::dbb2_dp(RealVectorValue vel,
     118             :                               RealVectorValue dvel_dp,
     119             :                               RealVectorValue xi_prime,
     120             :                               RealVectorValue eta_prime,
     121             :                               RealVectorValue zeta_prime) const
     122             : {
     123             :   return 2.0 *
     124     6278018 :          ((xi_prime * vel) * (dvel_dp * xi_prime) + (eta_prime * vel) * (dvel_dp * eta_prime) +
     125     6278018 :           (zeta_prime * vel) * (dvel_dp * zeta_prime));
     126             : }
     127             : 
     128             : Real
     129     6278018 : RichardsSUPGstandard::tauSUPG(RealVectorValue vel, Real traceperm, RealVectorValue b) const
     130             : {
     131             :   // vel = velocity, b = bb
     132     6278018 :   Real norm_v = vel.norm();
     133     6278018 :   Real norm_b = b.norm(); // Hughes et al investigate infinity-norm and 2-norm.  i just use 2-norm
     134             :                           // here.   norm_b ~ 2|a|/ele_length_in_direction_of_a
     135             : 
     136     6278018 :   if (norm_b == 0)
     137             :     return 0.0; // Only way for norm_b=0 is for zero ele size, or vel=0.  Either way we don't have
     138             :                 // to upwind.
     139             : 
     140     6233618 :   Real h = 2.0 * norm_v / norm_b; // h is a measure of the element length in the "a" direction
     141     6233618 :   Real alpha = 0.5 * norm_v * h / (traceperm * _p_SUPG); // this is the Peclet number
     142             : 
     143     6233618 :   const Real xi_tilde = RichardsSUPGstandard::cosh_relation(alpha);
     144             : 
     145     6233618 :   return xi_tilde / norm_b;
     146             : }
     147             : 
     148             : RealVectorValue
     149     6278018 : RichardsSUPGstandard::dtauSUPG_dgradp(RealVectorValue vel,
     150             :                                       RealTensorValue dvel_dgradp,
     151             :                                       Real traceperm,
     152             :                                       RealVectorValue b,
     153             :                                       RealVectorValue db2_dgradp) const
     154             : {
     155     6278018 :   Real norm_vel = vel.norm();
     156     6278018 :   if (norm_vel == 0)
     157             :     return RealVectorValue();
     158             : 
     159     6234218 :   RealVectorValue norm_vel_dgradp(dvel_dgradp * vel / norm_vel);
     160             : 
     161     6234218 :   Real norm_b = b.norm();
     162     6234218 :   if (norm_b == 0)
     163             :     return RealVectorValue();
     164             :   RealVectorValue norm_b_dgradp = db2_dgradp / 2.0 / norm_b;
     165             : 
     166     6233618 :   Real h = 2 * norm_vel / norm_b; // h is a measure of the element length in the "a" direction
     167             :   RealVectorValue h_dgradp(2 * norm_vel_dgradp / norm_b -
     168             :                            2.0 * norm_vel * norm_b_dgradp / norm_b / norm_b);
     169             : 
     170     6233618 :   Real alpha = 0.5 * norm_vel * h / traceperm / _p_SUPG; // this is the Peclet number
     171             :   RealVectorValue alpha_dgradp =
     172             :       0.5 * (norm_vel_dgradp * h + norm_vel * h_dgradp) / traceperm / _p_SUPG;
     173             : 
     174     6233618 :   Real xi_tilde = RichardsSUPGstandard::cosh_relation(alpha);
     175     6233618 :   Real xi_tilde_prime = RichardsSUPGstandard::cosh_relation_prime(alpha);
     176             :   RealVectorValue xi_tilde_dgradp = xi_tilde_prime * alpha_dgradp;
     177             : 
     178             :   RealVectorValue tau_dgradp =
     179     6233618 :       xi_tilde_dgradp / norm_b - xi_tilde * norm_b_dgradp / (norm_b * norm_b);
     180             : 
     181     6233618 :   return tau_dgradp;
     182             : }
     183             : 
     184             : Real
     185     6278018 : RichardsSUPGstandard::dtauSUPG_dp(RealVectorValue vel,
     186             :                                   RealVectorValue dvel_dp,
     187             :                                   Real traceperm,
     188             :                                   RealVectorValue b,
     189             :                                   Real db2_dp) const
     190             : {
     191     6278018 :   Real norm_vel = vel.norm();
     192     6278018 :   if (norm_vel == 0.0)
     193             :     return 0.0; // this deriv is not necessarily correct, but i can't see a better thing to do
     194             : 
     195     6234218 :   Real norm_vel_dp = dvel_dp * vel / norm_vel;
     196             : 
     197     6234218 :   Real norm_b = b.norm();
     198     6234218 :   if (norm_b == 0)
     199             :     return 0.0; // this deriv is not necessarily correct, but i can't see a better thing to do
     200     6233618 :   Real norm_b_dp = db2_dp / (2.0 * norm_b);
     201             : 
     202     6233618 :   Real h = 2.0 * norm_vel / norm_b; // h is a measure of the element length in the "a" direction
     203     6233618 :   Real h_dp = 2.0 * norm_vel_dp / norm_b - 2.0 * norm_vel * norm_b_dp / (norm_b * norm_b);
     204             : 
     205     6233618 :   Real alpha = 0.5 * norm_vel * h / (traceperm * _p_SUPG); // this is the Peclet number
     206     6233618 :   Real alpha_dp = 0.5 * (norm_vel_dp * h + norm_vel * h_dp) / (traceperm * _p_SUPG);
     207             : 
     208     6233618 :   Real xi_tilde = RichardsSUPGstandard::cosh_relation(alpha);
     209     6233618 :   Real xi_tilde_prime = RichardsSUPGstandard::cosh_relation_prime(alpha);
     210     6233618 :   Real xi_tilde_dp = xi_tilde_prime * alpha_dp;
     211             : 
     212             :   // Real tau = xi_tilde/norm_b;
     213     6233618 :   const Real tau_dp = xi_tilde_dp / norm_b - xi_tilde * norm_b_dp / (norm_b * norm_b);
     214             : 
     215     6233618 :   return tau_dp;
     216             : }
     217             : 
     218             : bool
     219     8034169 : RichardsSUPGstandard::SUPG_trivial() const
     220             : {
     221     8034169 :   return false;
     222             : }