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 "Q2PPorepressureFlux.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             : // C++ includes
      20             : #include <iostream>
      21             : 
      22             : registerMooseObject("RichardsApp", Q2PPorepressureFlux);
      23             : 
      24             : InputParameters
      25          28 : Q2PPorepressureFlux::validParams()
      26             : {
      27          28 :   InputParameters params = Kernel::validParams();
      28          56 :   params.addRequiredParam<UserObjectName>(
      29             :       "fluid_density",
      30             :       "A RichardsDensity UserObject that defines the fluid density as a function of pressure.");
      31          56 :   params.addRequiredCoupledVar("saturation_variable", "The variable representing fluid saturation");
      32          56 :   params.addRequiredParam<UserObjectName>(
      33             :       "fluid_relperm",
      34             :       "A RichardsRelPerm UserObject (eg RichardsRelPermPower) that defines the "
      35             :       "fluid relative permeability as a function of the saturation variable");
      36          56 :   params.addRequiredParam<Real>("fluid_viscosity", "The fluid dynamic viscosity");
      37          28 :   params.addClassDescription("Flux according to Darcy-Richards flow.  The Variable for this Kernel "
      38             :                              "should be the porepressure.");
      39          28 :   return params;
      40           0 : }
      41             : 
      42          14 : Q2PPorepressureFlux::Q2PPorepressureFlux(const InputParameters & parameters)
      43             :   : Kernel(parameters),
      44          14 :     _density(getUserObject<RichardsDensity>("fluid_density")),
      45          14 :     _sat(coupledDofValues("saturation_variable")),
      46          14 :     _sat_var(coupled("saturation_variable")),
      47          14 :     _relperm(getUserObject<RichardsRelPerm>("fluid_relperm")),
      48          28 :     _viscosity(getParam<Real>("fluid_viscosity")),
      49          28 :     _gravity(getMaterialProperty<RealVectorValue>("gravity")),
      50          28 :     _permeability(getMaterialProperty<RealTensorValue>("permeability")),
      51          14 :     _num_nodes(0),
      52          14 :     _mobility(0),
      53          14 :     _dmobility_dp(0),
      54          28 :     _dmobility_ds(0)
      55             : {
      56          14 : }
      57             : 
      58             : void
      59       12006 : Q2PPorepressureFlux::prepareNodalValues()
      60             : {
      61       12006 :   _num_nodes = _sat.size();
      62             : 
      63             :   Real density;
      64             :   Real ddensity_dp;
      65             :   Real relperm;
      66             :   Real drelperm_ds;
      67             : 
      68       12006 :   _mobility.resize(_num_nodes);
      69       12006 :   _dmobility_dp.resize(_num_nodes);
      70       12006 :   _dmobility_ds.resize(_num_nodes);
      71             : 
      72       40398 :   for (unsigned int nodenum = 0; nodenum < _num_nodes; ++nodenum)
      73             :   {
      74       28392 :     density = _density.density(_var.dofValues()[nodenum]);      // fluid density at the node
      75       28392 :     ddensity_dp = _density.ddensity(_var.dofValues()[nodenum]); // d(fluid density)/dP at the node
      76       28392 :     relperm = _relperm.relperm(_sat[nodenum]); // relative permeability of the fluid at node nodenum
      77       28392 :     drelperm_ds = _relperm.drelperm(_sat[nodenum]); // d(relperm)/dsat
      78             : 
      79             :     // calculate the mobility and its derivatives wrt P and S
      80       28392 :     _mobility[nodenum] = density * relperm / _viscosity;
      81       28392 :     _dmobility_dp[nodenum] = ddensity_dp * relperm / _viscosity;
      82       28392 :     _dmobility_ds[nodenum] = density * drelperm_ds / _viscosity;
      83             :   }
      84       12006 : }
      85             : 
      86             : Real
      87       88032 : Q2PPorepressureFlux::computeQpResidual()
      88             : {
      89             :   // note this is not the complete residual:
      90             :   // the upwind mobility parts get added in computeResidual
      91       88032 :   return _grad_test[_i][_qp] *
      92       88032 :          (_permeability[_qp] * (_grad_u[_qp] - _density.density(_u[_qp]) * _gravity[_qp]));
      93             : }
      94             : 
      95             : void
      96        5535 : Q2PPorepressureFlux::computeResidual()
      97             : {
      98        5535 :   upwind(true, false, 0);
      99        5535 : }
     100             : 
     101             : void
     102          33 : Q2PPorepressureFlux::computeJacobian()
     103             : {
     104          33 :   upwind(false, true, _var.number());
     105          33 : }
     106             : 
     107             : void
     108        6439 : Q2PPorepressureFlux::computeOffDiagJacobian(const unsigned int jvar)
     109             : {
     110        6439 :   upwind(false, true, jvar);
     111        6439 : }
     112             : 
     113             : Real
     114      211648 : Q2PPorepressureFlux::computeQpJac(unsigned int dvar)
     115             : {
     116             :   // this is just the derivative of the flux WITHOUT the upstream mobility terms
     117             :   // Those terms get added in during computeJacobian()
     118      211648 :   if (dvar == _var.number())
     119      106880 :     return _grad_test[_i][_qp] *
     120      106880 :            (_permeability[_qp] *
     121      106880 :             (_grad_phi[_j][_qp] - _density.ddensity(_u[_qp]) * _gravity[_qp] * _phi[_j][_qp]));
     122             :   else
     123             :     return 0;
     124             : }
     125             : 
     126             : void
     127       12007 : Q2PPorepressureFlux::upwind(bool compute_res, bool compute_jac, unsigned int jvar)
     128             : {
     129       12007 :   if (compute_jac && !(jvar == _var.number() || jvar == _sat_var))
     130           1 :     return;
     131             : 
     132             :   // calculate the mobility values and their derivatives
     133       12006 :   prepareNodalValues();
     134             : 
     135             :   // compute the residual without the mobility terms
     136             :   // Even if we are computing the jacobian we still need this
     137             :   // in order to see which nodes are upwind and which are downwind
     138       12006 :   prepareVectorTag(_assembly, _var.number());
     139             : 
     140       40398 :   for (_i = 0; _i < _test.size(); _i++)
     141      116424 :     for (_qp = 0; _qp < _qrule->n_points(); _qp++)
     142       88032 :       _local_re(_i) += _JxW[_qp] * _coord[_qp] * computeQpResidual();
     143             : 
     144       12006 :   if (compute_jac)
     145             :   {
     146        6471 :     prepareMatrixTag(_assembly, _var.number(), jvar);
     147       21547 :     for (_i = 0; _i < _test.size(); _i++)
     148       60916 :       for (_j = 0; _j < _phi.size(); _j++)
     149      257488 :         for (_qp = 0; _qp < _qrule->n_points(); _qp++)
     150      211648 :           _local_ke(_i, _j) += _JxW[_qp] * _coord[_qp] * computeQpJac(jvar);
     151             :   }
     152             : 
     153             :   // Now perform the upwinding by multiplying the residuals at the
     154             :   // upstream nodes by their mobilities
     155             :   //
     156             :   // The residual for the kernel is the darcy flux.
     157             :   // This is
     158             :   // R_i = int{mobility*flux_no_mob} = int{mobility*grad(pot)*permeability*grad(test_i)}
     159             :   // for node i.  where int is the integral over the element.
     160             :   // However, in fully-upwind, the first step is to take the mobility outside the
     161             :   // integral, which was done in the _local_re calculation above.
     162             :   //
     163             :   // NOTE: Physically _local_re(_i) is a measure of fluid flowing out of node i
     164             :   // If we had left in mobility, it would be exactly the mass flux flowing out of node i.
     165             :   //
     166             :   // This leads to the definition of upwinding:
     167             :   // ***
     168             :   // If _local_re(i) is positive then we use mobility_i.  That is
     169             :   // we use the upwind value of mobility.
     170             :   // ***
     171             :   //
     172             :   // The final subtle thing is we must also conserve fluid mass: the total mass
     173             :   // flowing out of node i must be the sum of the masses flowing
     174             :   // into the other nodes.
     175             : 
     176             :   // FIRST:
     177             :   // this is a dirty way of getting around precision loss problems
     178             :   // and problems at steadystate where upwinding oscillates from
     179             :   // node-to-node causing nonconvergence.
     180             :   // I'm not sure if i actually need to do this in moose.  Certainly
     181             :   // in cosflow it is necessary.
     182             :   // I will code a better algorithm if necessary
     183             :   bool reached_steady = true;
     184       19712 :   for (unsigned int nodenum = 0; nodenum < _num_nodes; ++nodenum)
     185             :   {
     186       19220 :     if (_local_re(nodenum) >= 1E-20)
     187             :     {
     188             :       reached_steady = false;
     189             :       break;
     190             :     }
     191             :   }
     192             :   reached_steady = false;
     193             : 
     194             :   // DEFINE VARIABLES USED TO ENSURE MASS CONSERVATION
     195             :   // total mass out - used for mass conservation
     196             :   Real total_mass_out = 0;
     197             :   // total flux in
     198             :   Real total_in = 0;
     199             : 
     200             :   // the following holds derivatives of these
     201             :   std::vector<Real> dtotal_mass_out;
     202             :   std::vector<Real> dtotal_in;
     203       12006 :   if (compute_jac)
     204             :   {
     205        6471 :     dtotal_mass_out.assign(_num_nodes, 0);
     206        6471 :     dtotal_in.assign(_num_nodes, 0);
     207             :   }
     208             : 
     209             :   // PERFORM THE UPWINDING!
     210       40398 :   for (unsigned int nodenum = 0; nodenum < _num_nodes; ++nodenum)
     211             :   {
     212       28392 :     if (_local_re(nodenum) >= 0 || reached_steady) // upstream node
     213             :     {
     214       14326 :       if (compute_jac)
     215             :       {
     216       31062 :         for (_j = 0; _j < _phi.size(); _j++)
     217       23448 :           _local_ke(nodenum, _j) *= _mobility[nodenum];
     218        7614 :         if (jvar == _var.number())
     219             :           // deriv wrt P
     220        3842 :           _local_ke(nodenum, nodenum) += _dmobility_dp[nodenum] * _local_re(nodenum);
     221             :         else
     222             :           // deriv wrt S
     223        3772 :           _local_ke(nodenum, nodenum) += _dmobility_ds[nodenum] * _local_re(nodenum);
     224       31062 :         for (_j = 0; _j < _phi.size(); _j++)
     225       23448 :           dtotal_mass_out[_j] += _local_ke(nodenum, _j);
     226             :       }
     227       14326 :       _local_re(nodenum) *= _mobility[nodenum];
     228       14326 :       total_mass_out += _local_re(nodenum);
     229             :     }
     230             :     else
     231             :     {
     232       14066 :       total_in -= _local_re(nodenum); // note the -= means the result is positive
     233       14066 :       if (compute_jac)
     234       29854 :         for (_j = 0; _j < _phi.size(); _j++)
     235       22392 :           dtotal_in[_j] -= _local_ke(nodenum, _j);
     236             :     }
     237             :   }
     238             : 
     239             :   // CONSERVE MASS
     240             :   // proportion the total_mass_out mass to the inflow nodes, weighting by their _local_re values
     241             :   if (!reached_steady)
     242       40398 :     for (unsigned int nodenum = 0; nodenum < _num_nodes; ++nodenum)
     243       28392 :       if (_local_re(nodenum) < 0)
     244             :       {
     245       14066 :         if (compute_jac)
     246       29854 :           for (_j = 0; _j < _phi.size(); _j++)
     247             :           {
     248       22392 :             _local_ke(nodenum, _j) *= total_mass_out / total_in;
     249       22392 :             _local_ke(nodenum, _j) +=
     250       22392 :                 _local_re(nodenum) * (dtotal_mass_out[_j] / total_in -
     251       22392 :                                       dtotal_in[_j] * total_mass_out / total_in / total_in);
     252             :           }
     253       14066 :         _local_re(nodenum) *= total_mass_out / total_in;
     254             :       }
     255             : 
     256             :   // ADD RESULTS TO RESIDUAL OR JACOBIAN
     257       12006 :   if (compute_res)
     258             :   {
     259        5535 :     accumulateTaggedLocalResidual();
     260        5535 :     if (_has_save_in)
     261         602 :       for (unsigned int i = 0; i < _save_in.size(); i++)
     262         301 :         _save_in[i]->sys().solution().add_vector(_local_re, _save_in[i]->dofIndices());
     263             :   }
     264             : 
     265       12006 :   if (compute_jac)
     266             :   {
     267        6471 :     accumulateTaggedLocalMatrix();
     268        6471 :     if (_has_diag_save_in && jvar == _var.number())
     269             :     {
     270             :       const unsigned int rows = _local_ke.m();
     271         103 :       DenseVector<Number> diag(rows);
     272         515 :       for (unsigned int i = 0; i < rows; i++)
     273         412 :         diag(i) = _local_ke(i, i);
     274             : 
     275         206 :       for (unsigned int i = 0; i < _diag_save_in.size(); i++)
     276         103 :         _diag_save_in[i]->sys().solution().add_vector(diag, _diag_save_in[i]->dofIndices());
     277             :     }
     278             :   }
     279       12006 : }