LCOV - code coverage report
Current view: top level - src/kernels - PorousFlowDarcyBase.C (source / functions) Hit Total Coverage
Test: porous_flow Test Coverage Lines: 195 215 90.7 %
Date: 2017-11-17 17:48:31 Functions: 13 17 76.5 %
Legend: Lines: hit not hit

          Line data    Source code
       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             : 
       8             : #include "PorousFlowDarcyBase.h"
       9             : 
      10             : // MOOSE includes
      11             : #include "Assembly.h"
      12             : #include "MooseMesh.h"
      13             : #include "MooseVariable.h"
      14             : #include "SystemBase.h"
      15             : 
      16             : #include "libmesh/quadrature.h"
      17             : 
      18             : template <>
      19             : InputParameters
      20         127 : validParams<PorousFlowDarcyBase>()
      21             : {
      22         127 :   InputParameters params = validParams<Kernel>();
      23         381 :   params.addRequiredParam<RealVectorValue>("gravity",
      24         127 :                                            "Gravitational acceleration vector downwards (m/s^2)");
      25         381 :   params.addRequiredParam<UserObjectName>(
      26         127 :       "PorousFlowDictator", "The UserObject that holds the list of PorousFlow variable names");
      27         381 :   params.addParam<unsigned>("full_upwind_threshold",
      28             :                             5,
      29             :                             "If, for each timestep, the number of "
      30             :                             "upwind-downwind swaps in an element is less than "
      31             :                             "this quantity, then full upwinding is used for that element.  "
      32         127 :                             "Otherwise the fallback scheme is employed.");
      33         381 :   MooseEnum fallback_enum("quick harmonic", "quick");
      34         381 :   params.addParam<MooseEnum>("fallback_scheme",
      35             :                              fallback_enum,
      36             :                              "quick: use nodal mobility without "
      37             :                              "preserving mass.  harmonic: use a "
      38             :                              "harmonic mean of nodal mobilities "
      39         127 :                              "and preserve fluid mass");
      40         254 :   params.addClassDescription("Fully-upwinded advective Darcy flux");
      41         127 :   return params;
      42             : }
      43             : 
      44         128 : PorousFlowDarcyBase::PorousFlowDarcyBase(const InputParameters & parameters)
      45             :   : Kernel(parameters),
      46         256 :     _permeability(getMaterialProperty<RealTensorValue>("PorousFlow_permeability_qp")),
      47             :     _dpermeability_dvar(
      48         256 :         getMaterialProperty<std::vector<RealTensorValue>>("dPorousFlow_permeability_qp_dvar")),
      49             :     _dpermeability_dgradvar(getMaterialProperty<std::vector<std::vector<RealTensorValue>>>(
      50         256 :         "dPorousFlow_permeability_qp_dgradvar")),
      51             :     _fluid_density_node(
      52         256 :         getMaterialProperty<std::vector<Real>>("PorousFlow_fluid_phase_density_nodal")),
      53             :     _dfluid_density_node_dvar(getMaterialProperty<std::vector<std::vector<Real>>>(
      54         256 :         "dPorousFlow_fluid_phase_density_nodal_dvar")),
      55         256 :     _fluid_density_qp(getMaterialProperty<std::vector<Real>>("PorousFlow_fluid_phase_density_qp")),
      56             :     _dfluid_density_qp_dvar(getMaterialProperty<std::vector<std::vector<Real>>>(
      57         256 :         "dPorousFlow_fluid_phase_density_qp_dvar")),
      58         256 :     _fluid_viscosity(getMaterialProperty<std::vector<Real>>("PorousFlow_viscosity_nodal")),
      59             :     _dfluid_viscosity_dvar(
      60         256 :         getMaterialProperty<std::vector<std::vector<Real>>>("dPorousFlow_viscosity_nodal_dvar")),
      61         256 :     _pp(getMaterialProperty<std::vector<Real>>("PorousFlow_porepressure_nodal")),
      62         256 :     _grad_p(getMaterialProperty<std::vector<RealGradient>>("PorousFlow_grad_porepressure_qp")),
      63             :     _dgrad_p_dgrad_var(getMaterialProperty<std::vector<std::vector<Real>>>(
      64         256 :         "dPorousFlow_grad_porepressure_qp_dgradvar")),
      65             :     _dgrad_p_dvar(getMaterialProperty<std::vector<std::vector<RealGradient>>>(
      66         256 :         "dPorousFlow_grad_porepressure_qp_dvar")),
      67         256 :     _porousflow_dictator(getUserObject<PorousFlowDictator>("PorousFlowDictator")),
      68         128 :     _num_phases(_porousflow_dictator.numPhases()),
      69         256 :     _gravity(getParam<RealVectorValue>("gravity")),
      70         384 :     _full_upwind_threshold(getParam<unsigned>("full_upwind_threshold")),
      71         256 :     _fallback_scheme(getParam<MooseEnum>("fallback_scheme").getEnum<FallbackEnum>()),
      72         128 :     _proto_flux(_num_phases),
      73         128 :     _jacobian(_num_phases),
      74             :     _num_upwinds(),
      75        2816 :     _num_downwinds()
      76             : {
      77             : #ifdef LIBMESH_HAVE_TBB_API
      78         128 :   if (libMesh::n_threads() > 1)
      79             :     mooseWarning("PorousFlowDarcyBase: num_upwinds and num_downwinds may not be computed "
      80             :                  "accurately when using TBB and greater than 1 thread");
      81             : #endif
      82         128 : }
      83             : 
      84             : void
      85        1538 : PorousFlowDarcyBase::timestepSetup()
      86             : {
      87        1538 :   Kernel::timestepSetup();
      88        3076 :   _num_upwinds = std::unordered_map<unsigned, std::vector<std::vector<unsigned>>>();
      89        3076 :   _num_downwinds = std::unordered_map<unsigned, std::vector<std::vector<unsigned>>>();
      90        1538 : }
      91             : 
      92             : Real
      93    18494136 : PorousFlowDarcyBase::darcyQp(unsigned int ph) const
      94             : {
      95    18494136 :   return _grad_test[_i][_qp] *
      96    92470680 :          (_permeability[_qp] * (_grad_p[_qp][ph] - _fluid_density_qp[_qp][ph] * _gravity));
      97             : }
      98             : 
      99             : Real
     100    45592216 : PorousFlowDarcyBase::darcyQpJacobian(unsigned int jvar, unsigned int ph) const
     101             : {
     102    45592216 :   if (_porousflow_dictator.notPorousFlowVariable(jvar))
     103             :     return 0.0;
     104             : 
     105    45592216 :   const unsigned int pvar = _porousflow_dictator.porousFlowVariableNum(jvar);
     106   136776648 :   RealVectorValue deriv = _dpermeability_dvar[_qp][pvar] * _phi[_j][_qp] *
     107   182368864 :                           (_grad_p[_qp][ph] - _fluid_density_qp[_qp][ph] * _gravity);
     108   319145512 :   for (unsigned i = 0; i < LIBMESH_DIM; ++i)
     109   410329944 :     deriv += _dpermeability_dgradvar[_qp][i][pvar] * _grad_phi[_j][_qp](i) *
     110   136776648 :              (_grad_p[_qp][ph] - _fluid_density_qp[_qp][ph] * _gravity);
     111   182368864 :   deriv += _permeability[_qp] * (_grad_phi[_j][_qp] * _dgrad_p_dgrad_var[_qp][ph][pvar] -
     112    91184432 :                                  _phi[_j][_qp] * _dfluid_density_qp_dvar[_qp][ph][pvar] * _gravity);
     113    91184432 :   deriv += _permeability[_qp] * (_dgrad_p_dvar[_qp][ph][pvar] * _phi[_j][_qp]);
     114    45592216 :   return _grad_test[_i][_qp] * deriv;
     115             : }
     116             : 
     117             : Real
     118           0 : PorousFlowDarcyBase::computeQpResidual()
     119             : {
     120           0 :   mooseError("PorousFlowDarcyBase: computeQpResidual called");
     121             :   return 0.0;
     122             : }
     123             : 
     124             : void
     125      713915 : PorousFlowDarcyBase::computeResidual()
     126             : {
     127      713915 :   computeResidualAndJacobian(JacRes::CALCULATE_RESIDUAL, 0.0);
     128      713915 : }
     129             : 
     130             : void
     131           0 : PorousFlowDarcyBase::computeJacobian()
     132             : {
     133           0 :   computeResidualAndJacobian(JacRes::CALCULATE_JACOBIAN, _var.number());
     134           0 : }
     135             : 
     136             : void
     137      476218 : PorousFlowDarcyBase::computeOffDiagJacobian(unsigned int jvar)
     138             : {
     139      476218 :   computeResidualAndJacobian(JacRes::CALCULATE_JACOBIAN, jvar);
     140      476218 : }
     141             : 
     142             : void
     143     1190133 : PorousFlowDarcyBase::computeResidualAndJacobian(JacRes res_or_jac, unsigned int jvar)
     144             : {
     145     1666351 :   if ((res_or_jac == JacRes::CALCULATE_JACOBIAN) &&
     146      476218 :       _porousflow_dictator.notPorousFlowVariable(jvar))
     147           0 :     return;
     148             : 
     149             :   /// The PorousFlow variable index corresponding to the variable number jvar
     150             :   const unsigned int pvar =
     151     1190133 :       ((res_or_jac == JacRes::CALCULATE_JACOBIAN) ? _porousflow_dictator.porousFlowVariableNum(jvar)
     152             :                                                   : 0);
     153             : 
     154     1190133 :   DenseMatrix<Number> & ke = _assembly.jacobianBlock(_var.number(), jvar);
     155     1190133 :   if ((ke.n() == 0) &&
     156             :       (res_or_jac == JacRes::CALCULATE_JACOBIAN)) // this removes a problem encountered in
     157             :                                                   // the initial timestep when
     158             :                                                   // use_displaced_mesh=true
     159             :     return;
     160             : 
     161             :   /// The number of nodes in the element
     162     1190133 :   const unsigned int num_nodes = _test.size();
     163             : 
     164             :   /// Compute the residual and jacobian without the mobility terms. Even if we are computing the Jacobian
     165             :   /// we still need this in order to see which nodes are upwind and which are downwind.
     166     4321563 :   for (unsigned ph = 0; ph < _num_phases; ++ph)
     167             :   {
     168     3131430 :     _proto_flux[ph].assign(num_nodes, 0);
     169    12757566 :     for (_qp = 0; _qp < _qrule->n_points(); _qp++)
     170             :     {
     171    23307204 :       for (_i = 0; _i < num_nodes; ++_i)
     172    73976544 :         _proto_flux[ph][_i] += _JxW[_qp] * _coord[_qp] * darcyQp(ph);
     173             :     }
     174             :   }
     175             : 
     176             :   // for this element, record whether each node is "upwind" or "downwind" (or neither)
     177     1190133 :   const unsigned elem = _current_elem->id();
     178     1190133 :   if (_num_upwinds.find(elem) == _num_upwinds.end())
     179             :   {
     180      248046 :     _num_upwinds[elem] = std::vector<std::vector<unsigned>>(_num_phases);
     181      248046 :     _num_downwinds[elem] = std::vector<std::vector<unsigned>>(_num_phases);
     182      416227 :     for (unsigned ph = 0; ph < _num_phases; ++ph)
     183             :     {
     184      292204 :       _num_upwinds[elem][ph].assign(num_nodes, 0);
     185      292204 :       _num_downwinds[elem][ph].assign(num_nodes, 0);
     186             :     }
     187             :   }
     188             :   // record the information once per nonlinear iteration
     189     1190133 :   if (res_or_jac == JacRes::CALCULATE_JACOBIAN && jvar == _var.number())
     190             :   {
     191     1108823 :     for (unsigned ph = 0; ph < _num_phases; ++ph)
     192             :     {
     193     2929521 :       for (unsigned nod = 0; nod < num_nodes; ++nod)
     194             :       {
     195     2525156 :         if (_proto_flux[ph][nod] > 0)
     196      635696 :           _num_upwinds[elem][ph][nod]++;
     197      626882 :         else if (_proto_flux[ph][nod] < 0)
     198      564476 :           _num_downwinds[elem][ph][nod]++;
     199             :       }
     200             :     }
     201             :   }
     202             : 
     203             :   // based on _num_upwinds and _num_downwinds, calculate the maximum number
     204             :   // of upwind-downwind swaps that have been encountered in this timestep
     205             :   // for this element
     206     1190133 :   std::vector<unsigned> max_swaps(_num_phases, 0);
     207     4321563 :   for (unsigned ph = 0; ph < _num_phases; ++ph)
     208             :   {
     209    11191851 :     for (unsigned nod = 0; nod < num_nodes; ++nod)
     210     9626136 :       max_swaps[ph] = std::max(
     211     4813068 :           max_swaps[ph], std::min(_num_upwinds[elem][ph][nod], _num_downwinds[elem][ph][nod]));
     212             :   }
     213             : 
     214             :   // size the _jacobian correctly and calculate it for the case residual = _proto_flux
     215     1190133 :   if (res_or_jac == JacRes::CALCULATE_JACOBIAN)
     216             :   {
     217     1845862 :     for (unsigned ph = 0; ph < _num_phases; ++ph)
     218             :     {
     219     1369644 :       _jacobian[ph].resize(ke.m());
     220     2831240 :       for (_i = 0; _i < _test.size(); _i++)
     221             :       {
     222     4292836 :         _jacobian[ph][_i].assign(ke.n(), 0.0);
     223    10857966 :         for (_j = 0; _j < _phi.size(); _j++)
     224   108607528 :           for (_qp = 0; _qp < _qrule->n_points(); _qp++)
     225   182368864 :             _jacobian[ph][_i][_j] += _JxW[_qp] * _coord[_qp] * darcyQpJacobian(jvar, ph);
     226             :       }
     227             :     }
     228             :   }
     229             : 
     230             :   // Loop over all the phases, computing the mass flux, which
     231             :   // gets placed into _proto_flux, and the derivative of this
     232             :   // which gets placed into _jacobian
     233     4321563 :   for (unsigned int ph = 0; ph < _num_phases; ++ph)
     234             :   {
     235     3131430 :     if (max_swaps[ph] < _full_upwind_threshold)
     236     1564808 :       fullyUpwind(res_or_jac, ph, pvar);
     237             :     else
     238             :     {
     239         907 :       switch (_fallback_scheme)
     240             :       {
     241             :         case FallbackEnum::QUICK:
     242         347 :           quickUpwind(res_or_jac, ph, pvar);
     243             :           break;
     244             :         case FallbackEnum::HARMONIC:
     245         560 :           harmonicMean(res_or_jac, ph, pvar);
     246             :           break;
     247             :       }
     248             :     }
     249             :   }
     250             : 
     251             :   /// Add results to the Residual or Jacobian
     252     1190133 :   if (res_or_jac == JacRes::CALCULATE_RESIDUAL)
     253             :   {
     254      713915 :     DenseVector<Number> & re = _assembly.residualBlock(_var.number());
     255             : 
     256      713915 :     _local_re.resize(re.size());
     257             :     _local_re.zero();
     258     2954453 :     for (_i = 0; _i < _test.size(); _i++)
     259     7573838 :       for (unsigned int ph = 0; ph < _num_phases; ++ph)
     260     7999950 :         _local_re(_i) += _proto_flux[ph][_i];
     261             : 
     262             :     re += _local_re;
     263             : 
     264      713915 :     if (_has_save_in)
     265             :     {
     266             :       Threads::spin_mutex::scoped_lock lock(Threads::spin_mtx);
     267           0 :       for (unsigned int i = 0; i < _save_in.size(); i++)
     268           0 :         _save_in[i]->sys().solution().add_vector(_local_re, _save_in[i]->dofIndices());
     269             :     }
     270             :   }
     271             : 
     272     1190133 :   if (res_or_jac == JacRes::CALCULATE_JACOBIAN)
     273             :   {
     274      476218 :     _local_ke.resize(ke.m(), ke.n());
     275             :     _local_ke.zero();
     276             : 
     277     2094380 :     for (_i = 0; _i < _test.size(); _i++)
     278     8829006 :       for (_j = 0; _j < _phi.size(); _j++)
     279    24633940 :         for (unsigned int ph = 0; ph < _num_phases; ++ph)
     280    26134644 :           _local_ke(_i, _j) += _jacobian[ph][_i][_j];
     281             : 
     282             :     ke += _local_ke;
     283             : 
     284      476218 :     if (_has_diag_save_in && jvar == _var.number())
     285             :     {
     286             :       unsigned int rows = ke.m();
     287             :       DenseVector<Number> diag(rows);
     288           0 :       for (unsigned int i = 0; i < rows; i++)
     289           0 :         diag(i) = _local_ke(i, i);
     290             : 
     291             :       Threads::spin_mutex::scoped_lock lock(Threads::spin_mtx);
     292           0 :       for (unsigned int i = 0; i < _diag_save_in.size(); i++)
     293           0 :         _diag_save_in[i]->sys().solution().add_vector(diag, _diag_save_in[i]->dofIndices());
     294             :     }
     295             :   }
     296             : }
     297             : 
     298             : void
     299     1564808 : PorousFlowDarcyBase::fullyUpwind(JacRes res_or_jac, unsigned int ph, unsigned int pvar)
     300             : {
     301             :   /**
     302             :    * Perform the full upwinding by multiplying the residuals at the upstream nodes by their
     303             :    * mobilities.
     304             :    * Mobility is different for each phase, and in each situation:
     305             :    *   mobility = density / viscosity    for single-component Darcy flow
     306             :    *   mobility = mass_fraction * density * relative_perm / viscosity    for multi-component,
     307             :    *multiphase flow
     308             :    *   mobility = enthalpy * density * relative_perm / viscosity    for heat convection
     309             :    *
     310             :    * The residual for the kernel is the sum over Darcy fluxes for each phase.
     311             :    * The Darcy flux for a particular phase is
     312             :    * R_i = int{mobility*flux_no_mob} = int{mobility*grad(pot)*permeability*grad(test_i)}
     313             :    * for node i.  where int is the integral over the element.
     314             :    * However, in fully-upwind, the first step is to take the mobility outside the integral,
     315             :    * which was done in the _proto_flux calculation above.
     316             :    *
     317             :    * NOTE: Physically _proto_flux[i][ph] is a measure of fluid of phase ph flowing out of node i.
     318             :    * If we had left in mobility, it would be exactly the component mass flux flowing out of node
     319             :    *i.
     320             :    *
     321             :    * This leads to the definition of upwinding:
     322             :    *
     323             :    * If _proto_flux(i)[ph] is positive then we use mobility_i.  That is we use the upwind value of
     324             :    * mobility.
     325             :    *
     326             :    * The final subtle thing is we must also conserve fluid mass: the total component mass flowing
     327             :    * out of node i must be the sum of the masses flowing into the other nodes.
     328             :   **/
     329             : 
     330             :   /// The number of nodes in the element
     331     1564808 :   const unsigned int num_nodes = _test.size();
     332             : 
     333             :   Real mob;
     334             :   Real dmob;
     335             :   /// Define variables used to ensure mass conservation
     336             :   Real total_mass_out = 0.0;
     337             :   Real total_in = 0.0;
     338             : 
     339             :   /// The following holds derivatives of these
     340             :   std::vector<Real> dtotal_mass_out;
     341             :   std::vector<Real> dtotal_in;
     342     1564808 :   if (res_or_jac == JacRes::CALCULATE_JACOBIAN)
     343             :   {
     344     1369382 :     dtotal_mass_out.assign(num_nodes, 0.0);
     345     1369382 :     dtotal_in.assign(num_nodes, 0.0);
     346             :   }
     347             : 
     348             :   /// Perform the upwinding using the mobility
     349     3129616 :   std::vector<bool> upwind_node(num_nodes);
     350    11183688 :   for (unsigned int n = 0; n < num_nodes; ++n)
     351             :   {
     352     9618880 :     if (_proto_flux[ph][n] >= 0.0) // upstream node
     353             :     {
     354             :       upwind_node[n] = true;
     355             :       /// The mobility at the upstream node
     356     2665307 :       mob = mobility(n, ph);
     357     2665307 :       if (res_or_jac == JacRes::CALCULATE_JACOBIAN)
     358             :       {
     359             :         /// The derivative of the mobility wrt the PorousFlow variable
     360     1159920 :         dmob = dmobility(n, ph, pvar);
     361             : 
     362     5777668 :         for (_j = 0; _j < _phi.size(); _j++)
     363     9235496 :           _jacobian[ph][n][_j] *= mob;
     364             : 
     365     1159920 :         if (_test.size() == _phi.size())
     366             :           /* mobility at node=n depends only on the variables at node=n, by construction.  For
     367             :            * linear-lagrange variables, this means that Jacobian entries involving the derivative
     368             :            * of mobility will only be nonzero for derivatives wrt variables at node=n.  Hence the
     369             :            * [n][n] in the line below.  However, for other variable types (eg constant monomials)
     370             :            * I cannot tell what variable number contributes to the derivative.  However, in all
     371             :            * cases I can possibly imagine, the derivative is zero anyway, since in the full
     372             :            * upwinding scheme, mobility shouldn't depend on these other sorts of variables.
     373             :            */
     374     2319776 :           _jacobian[ph][n][n] += dmob * _proto_flux[ph][n];
     375             : 
     376     5777668 :         for (_j = 0; _j < _phi.size(); _j++)
     377    13853244 :           dtotal_mass_out[_j] += _jacobian[ph][n][_j];
     378             :       }
     379     2665307 :       _proto_flux[ph][n] *= mob;
     380     2665307 :       total_mass_out += _proto_flux[ph][n];
     381             :     }
     382             :     else
     383             :     {
     384             :       upwind_node[n] = false;
     385     2144133 :       total_in -= _proto_flux[ph][n]; /// note the -= means the result is positive
     386     2144133 :       if (res_or_jac == JacRes::CALCULATE_JACOBIAN)
     387     5077678 :         for (_j = 0; _j < _phi.size(); _j++)
     388    12275112 :           dtotal_in[_j] -= _jacobian[ph][n][_j];
     389             :     }
     390             :   }
     391             : 
     392             :   /// Conserve mass over all phases by proportioning the total_mass_out mass to the inflow nodes, weighted by their proto_flux values
     393    11183688 :   for (unsigned int n = 0; n < num_nodes; ++n)
     394             :   {
     395     9618880 :     if (!upwind_node[n]) // downstream node
     396             :     {
     397     2144133 :       if (res_or_jac == JacRes::CALCULATE_JACOBIAN)
     398     5077678 :         for (_j = 0; _j < _phi.size(); _j++)
     399             :         {
     400     8183408 :           _jacobian[ph][n][_j] *= total_mass_out / total_in;
     401     4091704 :           _jacobian[ph][n][_j] +=
     402    12275112 :               _proto_flux[ph][n] * (dtotal_mass_out[_j] / total_in -
     403     8183408 :                                     dtotal_in[_j] * total_mass_out / total_in / total_in);
     404             :         }
     405     4288266 :       _proto_flux[ph][n] *= total_mass_out / total_in;
     406             :     }
     407             :   }
     408     1564808 : }
     409             : 
     410             : void
     411         347 : PorousFlowDarcyBase::quickUpwind(JacRes res_or_jac, unsigned int ph, unsigned int pvar)
     412             : {
     413             :   /// The number of nodes in the element
     414         347 :   const unsigned int num_nodes = _test.size();
     415             : 
     416             :   Real mob;
     417             :   Real dmob;
     418             : 
     419             :   /// Use the raw nodal mobility
     420        3123 :   for (unsigned int n = 0; n < num_nodes; ++n)
     421             :   {
     422             :     /// The mobility at the node
     423        1388 :     mob = mobility(n, ph);
     424        1388 :     if (res_or_jac == JacRes::CALCULATE_JACOBIAN)
     425             :     {
     426             :       /// The derivative of the mobility wrt the PorousFlow variable
     427         236 :       dmob = dmobility(n, ph, pvar);
     428             : 
     429        1180 :       for (_j = 0; _j < _phi.size(); _j++)
     430        1888 :         _jacobian[ph][n][_j] *= mob;
     431             : 
     432         236 :       if (_test.size() == _phi.size())
     433             :         /* mobility at node=n depends only on the variables at node=n, by construction.  For
     434             :          * linear-lagrange variables, this means that Jacobian entries involving the derivative
     435             :          * of mobility will only be nonzero for derivatives wrt variables at node=n.  Hence the
     436             :          * [n][n] in the line below.  However, for other variable types (eg constant monomials)
     437             :          * I cannot tell what variable number contributes to the derivative.  However, in all
     438             :          * cases I can possibly imagine, the derivative is zero anyway, since in the full
     439             :          * upwinding scheme, mobility shouldn't depend on these other sorts of variables.
     440             :          */
     441         708 :         _jacobian[ph][n][n] += dmob * _proto_flux[ph][n];
     442             :     }
     443        2776 :     _proto_flux[ph][n] *= mob;
     444             :   }
     445         347 : }
     446             : 
     447             : void
     448         560 : PorousFlowDarcyBase::harmonicMean(JacRes res_or_jac, unsigned int ph, unsigned int pvar)
     449             : {
     450             :   /// The number of nodes in the element
     451         560 :   const unsigned int num_nodes = _test.size();
     452             : 
     453         560 :   std::vector<Real> mob(num_nodes);
     454             :   unsigned num_zero = 0;
     455             :   unsigned zero_mobility_node;
     456             :   Real harmonic_mob = 0;
     457        5040 :   for (unsigned n = 0; n < num_nodes; ++n)
     458             :   {
     459        4480 :     mob[n] = mobility(n, ph);
     460        2240 :     if (mob[n] == 0.0)
     461             :     {
     462             :       zero_mobility_node = n;
     463         704 :       num_zero++;
     464             :     }
     465             :     else
     466        1536 :       harmonic_mob += 1.0 / mob[n];
     467             :   }
     468         560 :   if (num_zero > 0)
     469             :     harmonic_mob = 0.0;
     470             :   else
     471         376 :     harmonic_mob = (1.0 * num_nodes) / harmonic_mob;
     472             : 
     473             :   /*
     474             :   if (res_or_jac == JacRes::CALCULATE_RESIDUAL)
     475             :   {
     476             :     Moose::out << "RES ";
     477             :     for (unsigned n = 0; n < num_nodes; ++n)
     478             :       Moose::out << std::setprecision(16) << _pp[n][0] << " ";
     479             :     Moose::out << "\n";
     480             :     for (unsigned n = 0; n < num_nodes; ++n)
     481             :       Moose::out << mob[n] << " ";
     482             :     Moose::out << "harmonic_mob = " << harmonic_mob << "\n";
     483             :   }
     484             :   if (res_or_jac == JacRes::CALCULATE_JACOBIAN)
     485             :   {
     486             :     Moose::out << "JAC ";
     487             :     for (unsigned n = 0; n < num_nodes; ++n)
     488             :       Moose::out << mob[n] << " ";
     489             :     Moose::out << "harmonic_mob = " << harmonic_mob << "\n";
     490             :   }
     491             :   */
     492             : 
     493             :   // d(harmonic_mob)/d(PorousFlow variable at node n)
     494         560 :   std::vector<Real> dharmonic_mob(num_nodes, 0.0);
     495         560 :   if (res_or_jac == JacRes::CALCULATE_JACOBIAN)
     496             :   {
     497          72 :     const Real harm2 = std::pow(harmonic_mob, 2) / (1.0 * num_nodes);
     498          72 :     if (num_zero == 0)
     499         432 :       for (unsigned n = 0; n < num_nodes; ++n)
     500         768 :         dharmonic_mob[n] = dmobility(n, ph, pvar) * harm2 / std::pow(mob[n], 2);
     501          24 :     else if (num_zero == 1)
     502           0 :       dharmonic_mob[zero_mobility_node] =
     503           0 :           num_nodes * dmobility(zero_mobility_node, ph, pvar); // other derivs are zero
     504             :     // if num_zero > 1 then all dharmonic_mob = 0.0
     505             :   }
     506             : 
     507         560 :   if (res_or_jac == JacRes::CALCULATE_JACOBIAN)
     508         648 :     for (unsigned n = 0; n < num_nodes; ++n)
     509        1440 :       for (_j = 0; _j < _phi.size(); _j++)
     510             :       {
     511        2304 :         _jacobian[ph][n][_j] *= harmonic_mob;
     512        1152 :         if (_test.size() == _phi.size())
     513        2304 :           _jacobian[ph][n][_j] += dharmonic_mob[_j] * _proto_flux[ph][n];
     514             :       }
     515             : 
     516         560 :   if (res_or_jac == JacRes::CALCULATE_RESIDUAL)
     517        4392 :     for (unsigned n = 0; n < num_nodes; ++n)
     518        3904 :       _proto_flux[ph][n] *= harmonic_mob;
     519         560 : }
     520             : 
     521             : Real
     522           0 : PorousFlowDarcyBase::mobility(unsigned nodenum, unsigned phase) const
     523             : {
     524           0 :   return _fluid_density_node[nodenum][phase] / _fluid_viscosity[nodenum][phase];
     525             : }
     526             : 
     527             : Real
     528           0 : PorousFlowDarcyBase::dmobility(unsigned nodenum, unsigned phase, unsigned pvar) const
     529             : {
     530           0 :   Real dm = _dfluid_density_node_dvar[nodenum][phase][pvar] / _fluid_viscosity[nodenum][phase];
     531           0 :   dm -= _fluid_density_node[nodenum][phase] * _dfluid_viscosity_dvar[nodenum][phase][pvar] /
     532             :         std::pow(_fluid_viscosity[nodenum][phase], 2);
     533           0 :   return dm;
     534        2499 : }

Generated by: LCOV version 1.11