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 "PorousFlowFluidState.h"
      11             : 
      12             : registerMooseObject("PorousFlowApp", PorousFlowFluidState);
      13             : registerMooseObject("PorousFlowApp", ADPorousFlowFluidState);
      14             : 
      15             : template <bool is_ad>
      16             : InputParameters
      17        2434 : PorousFlowFluidStateTempl<is_ad>::validParams()
      18             : {
      19        2434 :   InputParameters params = PorousFlowFluidStateBaseMaterialTempl<is_ad>::validParams();
      20        4868 :   params.addRequiredCoupledVar("gas_porepressure",
      21             :                                "Variable that is the porepressure of the gas phase");
      22        4868 :   params.addRequiredCoupledVar("z", "Total mass fraction of component i summed over all phases");
      23        4868 :   params.addCoupledVar(
      24             :       "temperature", 20, "The fluid temperature (C or K, depending on temperature_unit)");
      25        4868 :   params.addCoupledVar("xnacl", 0, "The salt mass fraction in the brine (kg/kg)");
      26        4868 :   params.addRequiredParam<UserObjectName>("fluid_state", "Name of the FluidState UserObject");
      27        2434 :   params.addClassDescription("Class for fluid state calculations using persistent primary "
      28             :                              "variables and a vapor-liquid flash");
      29        2434 :   return params;
      30           0 : }
      31             : 
      32             : template <bool is_ad>
      33        1878 : PorousFlowFluidStateTempl<is_ad>::PorousFlowFluidStateTempl(const InputParameters & parameters)
      34             :   : PorousFlowFluidStateBaseMaterialTempl<is_ad>(parameters),
      35        3756 :     _gas_porepressure(_nodal_material
      36         729 :                           ? this->template coupledGenericDofValue<is_ad>("gas_porepressure")
      37        3027 :                           : this->template coupledGenericValue<is_ad>("gas_porepressure")),
      38        1878 :     _gas_gradp_qp(this->template coupledGenericGradient<is_ad>("gas_porepressure")),
      39        1878 :     _gas_porepressure_varnum(coupled("gas_porepressure")),
      40        3756 :     _pvar(_dictator.isPorousFlowVariable(_gas_porepressure_varnum)
      41        1878 :               ? _dictator.porousFlowVariableNum(_gas_porepressure_varnum)
      42             :               : 0),
      43        1878 :     _num_Z_vars(coupledComponents("z")),
      44        2940 :     _is_Xnacl_nodal(isCoupled("xnacl") ? getFieldVar("xnacl", 0)->isNodal() : false),
      45         729 :     _Xnacl(_nodal_material && _is_Xnacl_nodal
      46        2169 :                ? this->template coupledGenericDofValue<is_ad>("xnacl")
      47        3465 :                : this->template coupledGenericValue<is_ad>("xnacl")),
      48        1878 :     _grad_Xnacl_qp(this->template coupledGenericGradient<is_ad>("xnacl")),
      49        1878 :     _Xnacl_varnum(coupled("xnacl")),
      50        1878 :     _Xvar(_dictator.isPorousFlowVariable(_Xnacl_varnum)
      51        1878 :               ? _dictator.porousFlowVariableNum(_Xnacl_varnum)
      52             :               : 0),
      53        1878 :     _fs(this->template getUserObject<PorousFlowFluidStateMultiComponentBase>("fluid_state")),
      54        1878 :     _aqueous_phase_number(_fs.aqueousPhaseIndex()),
      55        1878 :     _gas_phase_number(_fs.gasPhaseIndex()),
      56        1878 :     _aqueous_fluid_component(_fs.aqueousComponentIndex()),
      57        1878 :     _gas_fluid_component(_fs.gasComponentIndex()),
      58        1878 :     _salt_component(_fs.saltComponentIndex()),
      59        1878 :     _temperature(
      60        1878 :         this->template getGenericMaterialProperty<Real, is_ad>("PorousFlow_temperature" + _sfx)),
      61        1878 :     _gradT_qp(_nodal_material ? nullptr
      62        1878 :                               : &this->template getGenericMaterialProperty<RealGradient, is_ad>(
      63             :                                     "PorousFlow_grad_temperature" + _sfx)),
      64        1878 :     _dtemperature_dvar(is_ad ? nullptr
      65        1578 :                              : &this->template getMaterialProperty<std::vector<Real>>(
      66             :                                    "dPorousFlow_temperature" + _sfx + "_dvar")),
      67        1878 :     _temperature_varnum(coupled("temperature")),
      68        1878 :     _Tvar(_dictator.isPorousFlowVariable(_temperature_varnum)
      69        1878 :               ? _dictator.porousFlowVariableNum(_temperature_varnum)
      70             :               : 0),
      71        1878 :     _pidx(_fs.getPressureIndex()),
      72        1878 :     _Tidx(_fs.getTemperatureIndex()),
      73        1878 :     _Zidx(_fs.getZIndex()),
      74        3756 :     _Xidx(_fs.getXIndex())
      75             : {
      76             :   // Check that the number of phases in the fluidstate class is also provided in the Dictator
      77        1878 :   if (_fs.numPhases() != _num_phases)
      78           0 :     mooseError(name(),
      79             :                ": only ",
      80           0 :                _fs.numPhases(),
      81             :                " phases are allowed. Please check the number of phases entered in the dictator is "
      82             :                "correct");
      83             : 
      84             :   // Store all total mass fractions and associated variable numbers
      85        1878 :   _Z.resize(_num_Z_vars);
      86        1878 :   _gradZ_qp.resize(_num_Z_vars);
      87        1878 :   _Z_varnum.resize(_num_Z_vars);
      88        1878 :   _Zvar.resize(_num_Z_vars);
      89             : 
      90        3756 :   for (unsigned int i = 0; i < _num_Z_vars; ++i)
      91             :   {
      92        3027 :     _Z[i] = (_nodal_material ? &this->template coupledGenericDofValue<is_ad>("z", i)
      93        3027 :                              : &this->template coupledGenericValue<is_ad>("z", i));
      94        1878 :     _gradZ_qp[i] = &this->template coupledGenericGradient<is_ad>("z", i);
      95        1878 :     _Z_varnum[i] = coupled("z", i);
      96        3756 :     _Zvar[i] = (_dictator.isPorousFlowVariable(_Z_varnum[i])
      97        1878 :                     ? _dictator.porousFlowVariableNum(_Z_varnum[i])
      98             :                     : 0);
      99             :   }
     100        1878 : }
     101             : 
     102             : template <bool is_ad>
     103             : void
     104     1075877 : PorousFlowFluidStateTempl<is_ad>::thermophysicalProperties()
     105             : {
     106             :   // The FluidProperty objects use temperature in K
     107     1075877 :   const GenericReal<is_ad> Tk = _temperature[_qp] + _T_c2k;
     108             : 
     109     1075877 :   _fs.clearFluidStateProperties(_fsp);
     110     1075877 :   _fs.thermophysicalProperties(_gas_porepressure[_qp], Tk, _Xnacl[_qp], (*_Z[0])[_qp], _qp, _fsp);
     111     1075789 : }
     112             : 
     113             : template <bool is_ad>
     114             : void
     115       26701 : PorousFlowFluidStateTempl<is_ad>::initQpStatefulProperties()
     116             : {
     117       26701 :   PorousFlowFluidStateBaseMaterialTempl<is_ad>::initQpStatefulProperties();
     118       26701 : }
     119             : 
     120             : template <bool is_ad>
     121             : void
     122     1049176 : PorousFlowFluidStateTempl<is_ad>::computeQpProperties()
     123             : {
     124     1049176 :   PorousFlowFluidStateBaseMaterialTempl<is_ad>::computeQpProperties();
     125             : 
     126             :   // If the material isn't AD, we need to compute the derivatives
     127             :   if (!is_ad)
     128             :   {
     129             :     // Derivative of properties wrt variables (calculated in fluid state class)
     130     2640681 :     for (unsigned int ph = 0; ph < _num_phases; ++ph)
     131             :     {
     132             :       // If porepressure is a PorousFlow variable (it usually is), add derivatives wrt
     133             :       // porepressure
     134     1760454 :       if (_dictator.isPorousFlowVariable(_gas_porepressure_varnum))
     135             :       {
     136     3520908 :         (*_dporepressure_dvar)[_qp][ph][_pvar] = _fsp[ph].pressure.derivatives()[_pidx];
     137     3520908 :         (*_dsaturation_dvar)[_qp][ph][_pvar] = _fsp[ph].saturation.derivatives()[_pidx];
     138     3520908 :         (*_dfluid_density_dvar)[_qp][ph][_pvar] = _fsp[ph].density.derivatives()[_pidx];
     139     3520908 :         (*_dfluid_viscosity_dvar)[_qp][ph][_pvar] = _fsp[ph].viscosity.derivatives()[_pidx];
     140     3520908 :         (*_dfluid_enthalpy_dvar)[_qp][ph][_pvar] = _fsp[ph].enthalpy.derivatives()[_pidx];
     141     1760454 :         (*_dfluid_internal_energy_dvar)[_qp][ph][_pvar] =
     142     1760454 :             _fsp[ph].internal_energy.derivatives()[_pidx];
     143             : 
     144     5799676 :         for (unsigned int comp = 0; comp < _num_components; ++comp)
     145     4039222 :           (*_dmass_frac_dvar)[_qp][ph][comp][_pvar] =
     146     4039222 :               _fsp[ph].mass_fraction[comp].derivatives()[_pidx];
     147             :       }
     148             : 
     149             :       // If Z is a PorousFlow variable (it usually is), add derivatives wrt Z
     150     1760454 :       if (_dictator.isPorousFlowVariable(_Z_varnum[0]))
     151             :       {
     152     3520908 :         (*_dporepressure_dvar)[_qp][ph][_Zvar[0]] = _fsp[ph].pressure.derivatives()[_Zidx];
     153     3520908 :         (*_dsaturation_dvar)[_qp][ph][_Zvar[0]] = _fsp[ph].saturation.derivatives()[_Zidx];
     154     3520908 :         (*_dfluid_density_dvar)[_qp][ph][_Zvar[0]] = _fsp[ph].density.derivatives()[_Zidx];
     155     3520908 :         (*_dfluid_viscosity_dvar)[_qp][ph][_Zvar[0]] = _fsp[ph].viscosity.derivatives()[_Zidx];
     156     3520908 :         (*_dfluid_enthalpy_dvar)[_qp][ph][_Zvar[0]] = _fsp[ph].enthalpy.derivatives()[_Zidx];
     157     1760454 :         (*_dfluid_internal_energy_dvar)[_qp][ph][_Zvar[0]] =
     158     1760454 :             _fsp[ph].internal_energy.derivatives()[_Zidx];
     159             : 
     160     5799676 :         for (unsigned int comp = 0; comp < _num_components; ++comp)
     161     4039222 :           (*_dmass_frac_dvar)[_qp][ph][comp][_Zvar[0]] =
     162     4039222 :               _fsp[ph].mass_fraction[comp].derivatives()[_Zidx];
     163             :       }
     164             : 
     165             :       // If temperature is a PorousFlow variable (nonisothermal case), add derivatives wrt
     166             :       // temperature
     167     1760454 :       if (_dictator.isPorousFlowVariable(_temperature_varnum))
     168             :       {
     169     1450248 :         (*_dporepressure_dvar)[_qp][ph][_Tvar] = _fsp[ph].pressure.derivatives()[_Tidx];
     170     1450248 :         (*_dsaturation_dvar)[_qp][ph][_Tvar] = _fsp[ph].saturation.derivatives()[_Tidx];
     171     1450248 :         (*_dfluid_density_dvar)[_qp][ph][_Tvar] = _fsp[ph].density.derivatives()[_Tidx];
     172     1450248 :         (*_dfluid_viscosity_dvar)[_qp][ph][_Tvar] = _fsp[ph].viscosity.derivatives()[_Tidx];
     173     1450248 :         (*_dfluid_enthalpy_dvar)[_qp][ph][_Tvar] = _fsp[ph].enthalpy.derivatives()[_Tidx];
     174      725124 :         (*_dfluid_internal_energy_dvar)[_qp][ph][_Tvar] =
     175      725124 :             _fsp[ph].internal_energy.derivatives()[_Tidx];
     176             : 
     177     2519726 :         for (unsigned int comp = 0; comp < _num_components; ++comp)
     178     1794602 :           (*_dmass_frac_dvar)[_qp][ph][comp][_Tvar] =
     179     1794602 :               _fsp[ph].mass_fraction[comp].derivatives()[_Tidx];
     180             :       }
     181             : 
     182             :       // If Xnacl is a PorousFlow variable, add derivatives wrt Xnacl
     183     1760454 :       if (_dictator.isPorousFlowVariable(_Xnacl_varnum))
     184             :       {
     185     1036628 :         (*_dporepressure_dvar)[_qp][ph][_Xvar] = _fsp[ph].pressure.derivatives()[_Xidx];
     186     1036628 :         (*_dsaturation_dvar)[_qp][ph][_Xvar] = _fsp[ph].saturation.derivatives()[_Xidx];
     187     1036628 :         (*_dfluid_density_dvar)[_qp][ph][_Xvar] += _fsp[ph].density.derivatives()[_Xidx];
     188     1036628 :         (*_dfluid_viscosity_dvar)[_qp][ph][_Xvar] += _fsp[ph].viscosity.derivatives()[_Xidx];
     189     1036628 :         (*_dfluid_enthalpy_dvar)[_qp][ph][_Xvar] = _fsp[ph].enthalpy.derivatives()[_Xidx];
     190      518314 :         (*_dfluid_internal_energy_dvar)[_qp][ph][_Xvar] =
     191      518314 :             _fsp[ph].internal_energy.derivatives()[_Xidx];
     192             : 
     193     2073256 :         for (unsigned int comp = 0; comp < _num_components; ++comp)
     194     1554942 :           (*_dmass_frac_dvar)[_qp][ph][comp][_Xvar] =
     195     1554942 :               _fsp[ph].mass_fraction[comp].derivatives()[_Xidx];
     196             :       }
     197             :     }
     198             :   }
     199             : 
     200             :   // If the material properties are being evaluated at the qps, calculate the gradients as well
     201             :   // Note: only nodal properties are evaluated in initQpStatefulProperties(), so no need to check
     202             :   // _is_initqp flag for qp properties
     203      880227 :   if (!_nodal_material)
     204             :     if constexpr (!is_ad)
     205             :     {
     206             :       // Derivatives of capillary pressure
     207      438725 :       const Real dpc = _pc.dCapillaryPressure(_fsp[_aqueous_phase_number].saturation.value(), _qp);
     208             :       const Real d2pc =
     209      438725 :           _pc.d2CapillaryPressure(_fsp[_aqueous_phase_number].saturation.value(), _qp);
     210             : 
     211             :       // Gradients of saturation and porepressure in all phases
     212      438725 :       (*_grads_qp)[_qp][_gas_phase_number] =
     213      438725 :           (*_dsaturation_dvar)[_qp][_gas_phase_number][_pvar] * _gas_gradp_qp[_qp] +
     214      438725 :           (*_dsaturation_dvar)[_qp][_gas_phase_number][_Zvar[0]] * (*_gradZ_qp[0])[_qp] +
     215      438725 :           (*_dsaturation_dvar)[_qp][_gas_phase_number][_Tvar] * (*_gradT_qp)[_qp];
     216      438725 :       (*_grads_qp)[_qp][_aqueous_phase_number] = -(*_grads_qp)[_qp][_gas_phase_number];
     217             : 
     218      438725 :       (*_gradp_qp)[_qp][_gas_phase_number] = _gas_gradp_qp[_qp];
     219      438725 :       (*_gradp_qp)[_qp][_aqueous_phase_number] =
     220             :           _gas_gradp_qp[_qp] - dpc * (*_grads_qp)[_qp][_aqueous_phase_number];
     221             : 
     222             :       // Gradients of mass fractions for each component in each phase
     223      438725 :       (*_grad_mass_frac_qp)[_qp][_aqueous_phase_number][_aqueous_fluid_component] =
     224      438725 :           _fsp[_aqueous_phase_number].mass_fraction[_aqueous_fluid_component].derivatives()[_pidx] *
     225      438725 :               _gas_gradp_qp[_qp] +
     226      438725 :           _fsp[_aqueous_phase_number].mass_fraction[_aqueous_fluid_component].derivatives()[_Zidx] *
     227      438725 :               (*_gradZ_qp[0])[_qp] +
     228      438725 :           _fsp[_aqueous_phase_number].mass_fraction[_aqueous_fluid_component].derivatives()[_Tidx] *
     229             :               (*_gradT_qp)[_qp];
     230      438725 :       (*_grad_mass_frac_qp)[_qp][_aqueous_phase_number][_gas_fluid_component] =
     231             :           -(*_grad_mass_frac_qp)[_qp][_aqueous_phase_number][_aqueous_fluid_component];
     232             : 
     233      438725 :       (*_grad_mass_frac_qp)[_qp][_gas_phase_number][_aqueous_fluid_component] =
     234      438725 :           _fsp[_gas_phase_number].mass_fraction[_aqueous_fluid_component].derivatives()[_pidx] *
     235      438725 :               _gas_gradp_qp[_qp] +
     236      438725 :           _fsp[_gas_phase_number].mass_fraction[_aqueous_fluid_component].derivatives()[_Zidx] *
     237      438725 :               (*_gradZ_qp[0])[_qp] +
     238      438725 :           _fsp[_gas_phase_number].mass_fraction[_aqueous_fluid_component].derivatives()[_Tidx] *
     239      438725 :               (*_gradT_qp)[_qp];
     240      438725 :       (*_grad_mass_frac_qp)[_qp][_gas_phase_number][_gas_fluid_component] =
     241             :           -(*_grad_mass_frac_qp)[_qp][_gas_phase_number][_aqueous_fluid_component];
     242             : 
     243             :       // Derivatives of gradients wrt variables
     244      438725 :       if (_dictator.isPorousFlowVariable(_gas_porepressure_varnum))
     245             :       {
     246     1316175 :         for (unsigned int ph = 0; ph < _num_phases; ++ph)
     247      877450 :           (*_dgradp_qp_dgradv)[_qp][ph][_pvar] = 1.0;
     248             : 
     249      438725 :         (*_dgradp_qp_dgradv)[_qp][_aqueous_phase_number][_pvar] +=
     250      438725 :             -dpc * (*_dsaturation_dvar)[_qp][_aqueous_phase_number][_pvar];
     251             : 
     252      438725 :         (*_dgradp_qp_dv)[_qp][_aqueous_phase_number][_pvar] =
     253      438725 :             -d2pc * (*_grads_qp)[_qp][_aqueous_phase_number] *
     254             :             (*_dsaturation_dvar)[_qp][_aqueous_phase_number][_pvar];
     255             :       }
     256             : 
     257      438725 :       if (_dictator.isPorousFlowVariable(_Z_varnum[0]))
     258             :       {
     259      438725 :         (*_dgradp_qp_dgradv)[_qp][_aqueous_phase_number][_Zvar[0]] =
     260      438725 :             -dpc * (*_dsaturation_dvar)[_qp][_aqueous_phase_number][_Zvar[0]];
     261             : 
     262      438725 :         (*_dgradp_qp_dv)[_qp][_aqueous_phase_number][_Zvar[0]] =
     263      438725 :             -d2pc * (*_grads_qp)[_qp][_aqueous_phase_number] *
     264             :             (*_dsaturation_dvar)[_qp][_aqueous_phase_number][_Zvar[0]];
     265             :       }
     266             : 
     267      438725 :       if (_dictator.isPorousFlowVariable(_temperature_varnum))
     268             :       {
     269      180470 :         (*_dgradp_qp_dgradv)[_qp][_aqueous_phase_number][_Tvar] =
     270      180470 :             -dpc * (*_dsaturation_dvar)[_qp][_aqueous_phase_number][_Tvar];
     271             : 
     272      180470 :         (*_dgradp_qp_dv)[_qp][_aqueous_phase_number][_Tvar] =
     273      180470 :             -d2pc * (*_grads_qp)[_qp][_aqueous_phase_number] *
     274             :             (*_dsaturation_dvar)[_qp][_aqueous_phase_number][_Tvar];
     275             :       }
     276             : 
     277             :       // If Xnacl is a PorousFlow variable, add gradients and derivatives wrt Xnacl
     278      438725 :       if (_dictator.isPorousFlowVariable(_Xnacl_varnum))
     279             :       {
     280      129279 :         (*_dgradp_qp_dgradv)[_qp][_aqueous_phase_number][_Xvar] =
     281      129279 :             -dpc * (*_dsaturation_dvar)[_qp][_aqueous_phase_number][_Xvar];
     282             : 
     283      129279 :         (*_grads_qp)[_qp][_aqueous_phase_number] +=
     284      129279 :             (*_dsaturation_dvar)[_qp][_aqueous_phase_number][_Xvar] * _grad_Xnacl_qp[_qp];
     285             : 
     286      129279 :         (*_grads_qp)[_qp][_gas_phase_number] -=
     287      129279 :             (*_dsaturation_dvar)[_qp][_aqueous_phase_number][_Xvar] * _grad_Xnacl_qp[_qp];
     288             : 
     289      129279 :         (*_gradp_qp)[_qp][_aqueous_phase_number] -=
     290      129279 :             dpc * (*_dsaturation_dvar)[_qp][_aqueous_phase_number][_Xvar] * _grad_Xnacl_qp[_qp];
     291             : 
     292      129279 :         (*_dgradp_qp_dv)[_qp][_aqueous_phase_number][_Xvar] =
     293      129279 :             -d2pc * (*_grads_qp)[_qp][_aqueous_phase_number] *
     294             :             (*_dsaturation_dvar)[_qp][_aqueous_phase_number][_Xvar];
     295             : 
     296      129279 :         (*_grad_mass_frac_qp)[_qp][_aqueous_phase_number][_salt_component] = _grad_Xnacl_qp[_qp];
     297      129279 :         (*_grad_mass_frac_qp)[_qp][_aqueous_phase_number][_aqueous_fluid_component] +=
     298             :             _fsp[_aqueous_phase_number]
     299      129279 :                 .mass_fraction[_aqueous_fluid_component]
     300      129279 :                 .derivatives()[_Xidx] *
     301             :             _grad_Xnacl_qp[_qp];
     302      129279 :         (*_grad_mass_frac_qp)[_qp][_aqueous_phase_number][_gas_fluid_component] -=
     303             :             _fsp[_aqueous_phase_number]
     304             :                 .mass_fraction[_aqueous_fluid_component]
     305      129279 :                 .derivatives()[_Xidx] *
     306      129279 :             _grad_Xnacl_qp[_qp];
     307      129279 :         (*_grad_mass_frac_qp)[_qp][_gas_phase_number][_aqueous_fluid_component] +=
     308      129279 :             _fsp[_gas_phase_number].mass_fraction[_aqueous_fluid_component].derivatives()[_Xidx] *
     309      129279 :             _grad_Xnacl_qp[_qp];
     310      129279 :         (*_grad_mass_frac_qp)[_qp][_gas_phase_number][_gas_fluid_component] -=
     311      258558 :             _fsp[_gas_phase_number].mass_fraction[_aqueous_fluid_component].derivatives()[_Xidx] *
     312      129279 :             _grad_Xnacl_qp[_qp];
     313             :       }
     314             :     }
     315     1049088 : }
     316             : 
     317             : template class PorousFlowFluidStateTempl<false>;
     318             : template class PorousFlowFluidStateTempl<true>;