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 "RichardsFullyUpwindFlux.h"
      11             : 
      12             : // MOOSE includes
      13             : #include "Assembly.h"
      14             : #include "MooseVariable.h"
      15             : #include "SystemBase.h"
      16             : 
      17             : #include "libmesh/quadrature.h"
      18             : 
      19             : registerMooseObject("RichardsApp", RichardsFullyUpwindFlux);
      20             : 
      21             : InputParameters
      22         194 : RichardsFullyUpwindFlux::validParams()
      23             : {
      24         194 :   InputParameters params = Kernel::validParams();
      25         388 :   params.addRequiredParam<std::vector<UserObjectName>>(
      26             :       "relperm_UO", "List of names of user objects that define relative permeability");
      27         388 :   params.addRequiredParam<std::vector<UserObjectName>>(
      28             :       "seff_UO",
      29             :       "List of name of user objects that define effective saturation as a function of "
      30             :       "pressure list");
      31         388 :   params.addRequiredParam<std::vector<UserObjectName>>(
      32             :       "density_UO", "List of names of user objects that define the fluid density");
      33         388 :   params.addRequiredParam<UserObjectName>(
      34             :       "richardsVarNames_UO", "The UserObject that holds the list of Richards variable names.");
      35         194 :   return params;
      36           0 : }
      37             : 
      38          97 : RichardsFullyUpwindFlux::RichardsFullyUpwindFlux(const InputParameters & parameters)
      39             :   : Kernel(parameters),
      40          97 :     _richards_name_UO(getUserObject<RichardsVarNames>("richardsVarNames_UO")),
      41          97 :     _num_p(_richards_name_UO.num_v()),
      42          97 :     _pvar(_richards_name_UO.richards_var_num(_var.number())),
      43          97 :     _density_UO(getUserObjectByName<RichardsDensity>(
      44         194 :         getParam<std::vector<UserObjectName>>("density_UO")[_pvar])),
      45          97 :     _seff_UO(
      46         194 :         getUserObjectByName<RichardsSeff>(getParam<std::vector<UserObjectName>>("seff_UO")[_pvar])),
      47          97 :     _relperm_UO(getUserObjectByName<RichardsRelPerm>(
      48         194 :         getParam<std::vector<UserObjectName>>("relperm_UO")[_pvar])),
      49         194 :     _viscosity(getMaterialProperty<std::vector<Real>>("viscosity")),
      50         194 :     _flux_no_mob(getMaterialProperty<std::vector<RealVectorValue>>("flux_no_mob")),
      51          97 :     _dflux_no_mob_dv(
      52          97 :         getMaterialProperty<std::vector<std::vector<RealVectorValue>>>("dflux_no_mob_dv")),
      53          97 :     _dflux_no_mob_dgradv(
      54          97 :         getMaterialProperty<std::vector<std::vector<RealTensorValue>>>("dflux_no_mob_dgradv")),
      55          97 :     _num_nodes(0),
      56          97 :     _mobility(0),
      57         194 :     _dmobility_dv(0)
      58             : {
      59          97 :   _ps_at_nodes.resize(_num_p);
      60         250 :   for (unsigned int pnum = 0; pnum < _num_p; ++pnum)
      61         153 :     _ps_at_nodes[pnum] = _richards_name_UO.nodal_var(pnum);
      62          97 : }
      63             : 
      64             : void
      65      696958 : RichardsFullyUpwindFlux::prepareNodalValues()
      66             : {
      67      696958 :   _num_nodes = (*_ps_at_nodes[_pvar]).size();
      68             : 
      69             :   Real p;
      70             :   Real density;
      71             :   Real ddensity_dp;
      72             :   Real seff;
      73             :   std::vector<Real> dseff_dp;
      74             :   Real relperm;
      75             :   Real drelperm_ds;
      76      696958 :   _mobility.resize(_num_nodes);
      77      696958 :   _dmobility_dv.resize(_num_nodes);
      78      696958 :   dseff_dp.resize(_num_p);
      79     2717996 :   for (unsigned int nodenum = 0; nodenum < _num_nodes; ++nodenum)
      80             :   {
      81             :     // retrieve and calculate basic things at the node
      82     2021038 :     p = (*_ps_at_nodes[_pvar])[nodenum];   // pressure of fluid _pvar at node nodenum
      83     2021038 :     density = _density_UO.density(p);      // density of fluid _pvar at node nodenum
      84     2021038 :     ddensity_dp = _density_UO.ddensity(p); // d(density)/dP
      85             :     seff =
      86     2021038 :         _seff_UO.seff(_ps_at_nodes, nodenum); // effective saturation of fluid _pvar at node nodenum
      87     2021038 :     _seff_UO.dseff(_ps_at_nodes, nodenum, dseff_dp); // d(seff)/d(P_ph), for ph = 0, ..., _num_p - 1
      88     2021038 :     relperm = _relperm_UO.relperm(seff); // relative permeability of fluid _pvar at node nodenum
      89     2021038 :     drelperm_ds = _relperm_UO.drelperm(seff); // d(relperm)/dseff
      90             : 
      91             :     // calculate the mobility and its derivatives wrt (variable_ph = porepressure_ph)
      92     2021038 :     _mobility[nodenum] =
      93     2021038 :         density * relperm / _viscosity[0][_pvar]; // assume viscosity is constant throughout element
      94     2021038 :     _dmobility_dv[nodenum].resize(_num_p);
      95     4991356 :     for (unsigned int ph = 0; ph < _num_p; ++ph)
      96     2970318 :       _dmobility_dv[nodenum][ph] = density * drelperm_ds * dseff_dp[ph] / _viscosity[0][_pvar];
      97     2021038 :     _dmobility_dv[nodenum][_pvar] += ddensity_dp * relperm / _viscosity[0][_pvar];
      98             :   }
      99      696958 : }
     100             : 
     101             : Real
     102     7731052 : RichardsFullyUpwindFlux::computeQpResidual()
     103             : {
     104             :   // note this is not the complete residual:
     105             :   // the upwind mobility parts get added in computeResidual
     106     7731052 :   return _grad_test[_i][_qp] * _flux_no_mob[_qp][_pvar];
     107             : }
     108             : 
     109             : void
     110      430110 : RichardsFullyUpwindFlux::computeResidual()
     111             : {
     112      430110 :   upwind(true, false, 0);
     113      430110 :   return;
     114             : }
     115             : 
     116             : void
     117      266850 : RichardsFullyUpwindFlux::computeOffDiagJacobian(const unsigned int jvar)
     118             : {
     119      266850 :   upwind(false, true, jvar);
     120      266850 :   return;
     121             : }
     122             : 
     123             : Real
     124    18334856 : RichardsFullyUpwindFlux::computeQpJac(unsigned int dvar)
     125             : {
     126             :   // this is just the derivative of the flux WITHOUT the upstream mobility terms
     127             :   // Those terms get added in during computeJacobian()
     128    18334856 :   return _grad_test[_i][_qp] * (_dflux_no_mob_dgradv[_qp][_pvar][dvar] * _grad_phi[_j][_qp] +
     129    18334856 :                                 _dflux_no_mob_dv[_qp][_pvar][dvar] * _phi[_j][_qp]);
     130             : }
     131             : 
     132             : void
     133      696960 : RichardsFullyUpwindFlux::upwind(bool compute_res, bool compute_jac, unsigned int jvar)
     134             : {
     135      696960 :   if (compute_jac && _richards_name_UO.not_richards_var(jvar))
     136           2 :     return;
     137             : 
     138             :   // calculate the mobility values and their derivatives
     139      696958 :   prepareNodalValues();
     140             : 
     141             :   // compute the residual without the mobility terms
     142             :   // Even if we are computing the jacobian we still need this
     143             :   // in order to see which nodes are upwind and which are downwind
     144      696958 :   prepareVectorTag(_assembly, _var.number());
     145     2717996 :   for (_i = 0; _i < _test.size(); _i++)
     146     9752090 :     for (_qp = 0; _qp < _qrule->n_points(); _qp++)
     147     7731052 :       _local_re(_i) += _JxW[_qp] * _coord[_qp] * computeQpResidual();
     148             : 
     149      696958 :   const unsigned int dvar = _richards_name_UO.richards_var_num(jvar);
     150      696958 :   if (compute_jac)
     151             :   {
     152      266848 :     prepareMatrixTag(_assembly, _var.number(), jvar);
     153     1127802 :     for (_i = 0; _i < _test.size(); _i++)
     154     4394526 :       for (_j = 0; _j < _phi.size(); _j++)
     155    21868428 :         for (_qp = 0; _qp < _qrule->n_points(); _qp++)
     156    18334856 :           _local_ke(_i, _j) += _JxW[_qp] * _coord[_qp] * computeQpJac(dvar);
     157             :   }
     158             : 
     159             :   // Now perform the upwinding by multiplying the residuals at the
     160             :   // upstream nodes by their mobilities
     161             :   //
     162             :   // The residual for the kernel is the darcy flux.
     163             :   // This is
     164             :   // R_i = int{mobility*flux_no_mob} = int{mobility*grad(pot)*permeability*grad(test_i)}
     165             :   // for node i.  where int is the integral over the element.
     166             :   // However, in fully-upwind, the first step is to take the mobility outside the
     167             :   // integral, which was done in the _local_re calculation above.
     168             :   //
     169             :   // NOTE: Physically _local_re(_i) is a measure of fluid flowing out of node i
     170             :   // If we had left in mobility, it would be exactly the mass flux flowing out of node i.
     171             :   //
     172             :   // This leads to the definition of upwinding:
     173             :   // ***
     174             :   // If _local_re(i) is positive then we use mobility_i.  That is
     175             :   // we use the upwind value of mobility.
     176             :   // ***
     177             :   //
     178             :   // The final subtle thing is we must also conserve fluid mass: the total mass
     179             :   // flowing out of node i must be the sum of the masses flowing
     180             :   // into the other nodes.
     181             : 
     182             :   // FIRST:
     183             :   // this is a dirty way of getting around precision loss problems
     184             :   // and problems at steadystate where upwinding oscillates from
     185             :   // node-to-node causing nonconvergence.
     186             :   // I'm not sure if i actually need to do this in moose.  Certainly
     187             :   // in cosflow it is necessary.
     188             :   // I will code a better algorithm if necessary
     189             :   bool reached_steady = true;
     190      976281 :   for (unsigned int nodenum = 0; nodenum < _num_nodes; ++nodenum)
     191             :   {
     192      962263 :     if (_local_re(nodenum) >= 1E-20)
     193             :     {
     194             :       reached_steady = false;
     195             :       break;
     196             :     }
     197             :   }
     198             : 
     199             :   // DEFINE VARIABLES USED TO ENSURE MASS CONSERVATION
     200             :   // total mass out - used for mass conservation
     201             :   Real total_mass_out = 0;
     202             :   // total flux in
     203             :   Real total_in = 0;
     204             : 
     205             :   // the following holds derivatives of these
     206             :   std::vector<Real> dtotal_mass_out;
     207             :   std::vector<Real> dtotal_in;
     208      696958 :   if (compute_jac)
     209             :   {
     210      266848 :     dtotal_mass_out.resize(_num_nodes);
     211      266848 :     dtotal_in.resize(_num_nodes);
     212     1127802 :     for (unsigned int nodenum = 0; nodenum < _num_nodes; ++nodenum)
     213             :     {
     214      860954 :       dtotal_mass_out[nodenum] = 0;
     215      860954 :       dtotal_in[nodenum] = 0;
     216             :     }
     217             :   }
     218             : 
     219             :   // PERFORM THE UPWINDING!
     220     2717996 :   for (unsigned int nodenum = 0; nodenum < _num_nodes; ++nodenum)
     221             :   {
     222     2021038 :     if (_local_re(nodenum) >= 0 || reached_steady) // upstream node
     223             :     {
     224     1030463 :       if (compute_jac)
     225             :       {
     226     2247910 :         for (_j = 0; _j < _phi.size(); _j++)
     227     1805722 :           _local_ke(nodenum, _j) *= _mobility[nodenum];
     228      442188 :         _local_ke(nodenum, nodenum) += _dmobility_dv[nodenum][dvar] * _local_re(nodenum);
     229     2247910 :         for (_j = 0; _j < _phi.size(); _j++)
     230     1805722 :           dtotal_mass_out[_j] += _local_ke(nodenum, _j);
     231             :       }
     232     1030463 :       _local_re(nodenum) *= _mobility[nodenum];
     233     1030463 :       total_mass_out += _local_re(nodenum);
     234             :     }
     235             :     else
     236             :     {
     237      990575 :       total_in -= _local_re(nodenum); // note the -= means the result is positive
     238      990575 :       if (compute_jac)
     239     2146616 :         for (_j = 0; _j < _phi.size(); _j++)
     240     1727850 :           dtotal_in[_j] -= _local_ke(nodenum, _j);
     241             :     }
     242             :   }
     243             : 
     244             :   // CONSERVE MASS
     245             :   // proportion the total_mass_out mass to the inflow nodes, weighting by their _local_re values
     246      696958 :   if (!reached_steady)
     247     2671142 :     for (unsigned int nodenum = 0; nodenum < _num_nodes; ++nodenum)
     248     1988202 :       if (_local_re(nodenum) < 0)
     249             :       {
     250      990575 :         if (compute_jac)
     251     2146616 :           for (_j = 0; _j < _phi.size(); _j++)
     252             :           {
     253     1727850 :             _local_ke(nodenum, _j) *= total_mass_out / total_in;
     254     1727850 :             _local_ke(nodenum, _j) +=
     255     1727850 :                 _local_re(nodenum) * (dtotal_mass_out[_j] / total_in -
     256     1727850 :                                       dtotal_in[_j] * total_mass_out / total_in / total_in);
     257             :           }
     258      990575 :         _local_re(nodenum) *= total_mass_out / total_in;
     259             :       }
     260             : 
     261             :   // ADD RESULTS TO RESIDUAL OR JACOBIAN
     262      696958 :   if (compute_res)
     263             :   {
     264      430110 :     accumulateTaggedLocalResidual();
     265      430110 :     if (_has_save_in)
     266           8 :       for (unsigned int i = 0; i < _save_in.size(); i++)
     267           4 :         _save_in[i]->sys().solution().add_vector(_local_re, _save_in[i]->dofIndices());
     268             :   }
     269             : 
     270      696958 :   if (compute_jac)
     271             :   {
     272      266848 :     accumulateTaggedLocalMatrix();
     273      266848 :     if (_has_diag_save_in && dvar == _pvar)
     274             :     {
     275             :       const unsigned int rows = _local_ke.m();
     276           2 :       DenseVector<Number> diag(rows);
     277           6 :       for (unsigned int i = 0; i < rows; i++)
     278           4 :         diag(i) = _local_ke(i, i);
     279             : 
     280           4 :       for (unsigned int i = 0; i < _diag_save_in.size(); i++)
     281           2 :         _diag_save_in[i]->sys().solution().add_vector(diag, _diag_save_in[i]->dofIndices());
     282             :     }
     283             :   }
     284      696958 : }