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        4922 : PorousFlowFluidStateTempl<is_ad>::validParams()
      18             : {
      19        4922 :   InputParameters params = PorousFlowFluidStateBaseMaterialTempl<is_ad>::validParams();
      20        9844 :   params.addRequiredCoupledVar("gas_porepressure",
      21             :                                "Variable that is the porepressure of the gas phase");
      22        9844 :   params.addRequiredCoupledVar("z", "Total mass fraction of component i summed over all phases");
      23        9844 :   params.addCoupledVar(
      24             :       "temperature", 20, "The fluid temperature (C or K, depending on temperature_unit)");
      25        9844 :   params.addCoupledVar("xnacl", 0, "The salt mass fraction in the brine (kg/kg)");
      26        9844 :   params.addRequiredParam<UserObjectName>("fluid_state", "Name of the FluidState UserObject");
      27        4922 :   params.addClassDescription("Class for fluid state calculations using persistent primary "
      28             :                              "variables and a vapor-liquid flash");
      29        4922 :   return params;
      30           0 : }
      31             : 
      32             : template <bool is_ad>
      33        3846 : PorousFlowFluidStateTempl<is_ad>::PorousFlowFluidStateTempl(const InputParameters & parameters)
      34             :   : PorousFlowFluidStateBaseMaterialTempl<is_ad>(parameters),
      35        7692 :     _gas_porepressure(_nodal_material
      36        1461 :                           ? this->template coupledGenericDofValue<is_ad>("gas_porepressure")
      37        6231 :                           : this->template coupledGenericValue<is_ad>("gas_porepressure")),
      38        3846 :     _gas_gradp_qp(this->template coupledGenericGradient<is_ad>("gas_porepressure")),
      39        3846 :     _gas_porepressure_varnum(coupled("gas_porepressure")),
      40        7692 :     _pvar(_dictator.isPorousFlowVariable(_gas_porepressure_varnum)
      41        3846 :               ? _dictator.porousFlowVariableNum(_gas_porepressure_varnum)
      42             :               : 0),
      43        3846 :     _num_Z_vars(coupledComponents("z")),
      44        5988 :     _is_Xnacl_nodal(isCoupled("xnacl") ? getFieldVar("xnacl", 0)->isNodal() : false),
      45        1461 :     _Xnacl(_nodal_material && _is_Xnacl_nodal
      46        4389 :                ? this->template coupledGenericDofValue<is_ad>("xnacl")
      47        7149 :                : this->template coupledGenericValue<is_ad>("xnacl")),
      48        3846 :     _grad_Xnacl_qp(this->template coupledGenericGradient<is_ad>("xnacl")),
      49        3846 :     _Xnacl_varnum(coupled("xnacl")),
      50        3846 :     _Xvar(_dictator.isPorousFlowVariable(_Xnacl_varnum)
      51        3846 :               ? _dictator.porousFlowVariableNum(_Xnacl_varnum)
      52             :               : 0),
      53        3846 :     _fs(this->template getUserObject<PorousFlowFluidStateMultiComponentBase>("fluid_state")),
      54        3846 :     _aqueous_phase_number(_fs.aqueousPhaseIndex()),
      55        3846 :     _gas_phase_number(_fs.gasPhaseIndex()),
      56        3846 :     _aqueous_fluid_component(_fs.aqueousComponentIndex()),
      57        3846 :     _gas_fluid_component(_fs.gasComponentIndex()),
      58        3846 :     _salt_component(_fs.saltComponentIndex()),
      59        3846 :     _temperature(
      60        3846 :         this->template getGenericMaterialProperty<Real, is_ad>("PorousFlow_temperature" + _sfx)),
      61        3846 :     _gradT_qp(_nodal_material ? nullptr
      62        3846 :                               : &this->template getGenericMaterialProperty<RealGradient, is_ad>(
      63             :                                     "PorousFlow_grad_temperature" + _sfx)),
      64        3846 :     _dtemperature_dvar(is_ad ? nullptr
      65        3186 :                              : &this->template getMaterialProperty<std::vector<Real>>(
      66             :                                    "dPorousFlow_temperature" + _sfx + "_dvar")),
      67        3846 :     _temperature_varnum(coupled("temperature")),
      68        3846 :     _Tvar(_dictator.isPorousFlowVariable(_temperature_varnum)
      69        3846 :               ? _dictator.porousFlowVariableNum(_temperature_varnum)
      70             :               : 0),
      71        3846 :     _pidx(_fs.getPressureIndex()),
      72        3846 :     _Tidx(_fs.getTemperatureIndex()),
      73        3846 :     _Zidx(_fs.getZIndex()),
      74        7692 :     _Xidx(_fs.getXIndex())
      75             : {
      76             :   // Check that the number of phases in the fluidstate class is also provided in the Dictator
      77        3846 :   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        3846 :   _Z.resize(_num_Z_vars);
      86        3846 :   _gradZ_qp.resize(_num_Z_vars);
      87        3846 :   _Z_varnum.resize(_num_Z_vars);
      88        3846 :   _Zvar.resize(_num_Z_vars);
      89             : 
      90        7692 :   for (unsigned int i = 0; i < _num_Z_vars; ++i)
      91             :   {
      92        6231 :     _Z[i] = (_nodal_material ? &this->template coupledGenericDofValue<is_ad>("z", i)
      93        6231 :                              : &this->template coupledGenericValue<is_ad>("z", i));
      94        3846 :     _gradZ_qp[i] = &this->template coupledGenericGradient<is_ad>("z", i);
      95        3846 :     _Z_varnum[i] = coupled("z", i);
      96        7692 :     _Zvar[i] = (_dictator.isPorousFlowVariable(_Z_varnum[i])
      97        3846 :                     ? _dictator.porousFlowVariableNum(_Z_varnum[i])
      98             :                     : 0);
      99             :   }
     100        3846 : }
     101             : 
     102             : template <bool is_ad>
     103             : void
     104     1548562 : PorousFlowFluidStateTempl<is_ad>::thermophysicalProperties()
     105             : {
     106             :   // The FluidProperty objects use temperature in K
     107     1548562 :   const GenericReal<is_ad> Tk = _temperature[_qp] + _T_c2k;
     108             : 
     109     1548562 :   _fs.clearFluidStateProperties(_fsp);
     110     1548562 :   _fs.thermophysicalProperties(_gas_porepressure[_qp], Tk, _Xnacl[_qp], (*_Z[0])[_qp], _qp, _fsp);
     111     1548428 : }
     112             : 
     113             : template <bool is_ad>
     114             : void
     115       42541 : PorousFlowFluidStateTempl<is_ad>::initQpStatefulProperties()
     116             : {
     117       42541 :   PorousFlowFluidStateBaseMaterialTempl<is_ad>::initQpStatefulProperties();
     118       42541 : }
     119             : 
     120             : template <bool is_ad>
     121             : void
     122     1506021 : PorousFlowFluidStateTempl<is_ad>::computeQpProperties()
     123             : {
     124     1506021 :   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     3791811 :     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     2527874 :       if (_dictator.isPorousFlowVariable(_gas_porepressure_varnum))
     135             :       {
     136     5041672 :         (*_dporepressure_dvar)[_qp][ph][_pvar] = _fsp[ph].pressure.derivatives()[_pidx];
     137     2645478 :         (*_dsaturation_dvar)[_qp][ph][_pvar] = _fsp[ph].saturation.derivatives()[_pidx];
     138     3843575 :         (*_dfluid_density_dvar)[_qp][ph][_pvar] = _fsp[ph].density.derivatives()[_pidx];
     139     3843575 :         (*_dfluid_viscosity_dvar)[_qp][ph][_pvar] = _fsp[ph].viscosity.derivatives()[_pidx];
     140     3843575 :         (*_dfluid_enthalpy_dvar)[_qp][ph][_pvar] = _fsp[ph].enthalpy.derivatives()[_pidx];
     141     2527874 :         (*_dfluid_internal_energy_dvar)[_qp][ph][_pvar] =
     142     2527874 :             _fsp[ph].internal_energy.derivatives()[_pidx];
     143             : 
     144     8275182 :         for (unsigned int comp = 0; comp < _num_components; ++comp)
     145     5747308 :           (*_dmass_frac_dvar)[_qp][ph][comp][_pvar] =
     146     5747308 :               _fsp[ph].mass_fraction[comp].derivatives()[_pidx];
     147             :       }
     148             : 
     149             :       // If Z is a PorousFlow variable (it usually is), add derivatives wrt Z
     150     2527874 :       if (_dictator.isPorousFlowVariable(_Z_varnum[0]))
     151             :       {
     152     2586676 :         (*_dporepressure_dvar)[_qp][ph][_Zvar[0]] = _fsp[ph].pressure.derivatives()[_Zidx];
     153     2645478 :         (*_dsaturation_dvar)[_qp][ph][_Zvar[0]] = _fsp[ph].saturation.derivatives()[_Zidx];
     154     2933569 :         (*_dfluid_density_dvar)[_qp][ph][_Zvar[0]] = _fsp[ph].density.derivatives()[_Zidx];
     155     2594436 :         (*_dfluid_viscosity_dvar)[_qp][ph][_Zvar[0]] = _fsp[ph].viscosity.derivatives()[_Zidx];
     156     3777113 :         (*_dfluid_enthalpy_dvar)[_qp][ph][_Zvar[0]] = _fsp[ph].enthalpy.derivatives()[_Zidx];
     157     2527874 :         (*_dfluid_internal_energy_dvar)[_qp][ph][_Zvar[0]] =
     158     2527874 :             _fsp[ph].internal_energy.derivatives()[_Zidx];
     159             : 
     160     8275182 :         for (unsigned int comp = 0; comp < _num_components; ++comp)
     161     5747308 :           (*_dmass_frac_dvar)[_qp][ph][comp][_Zvar[0]] =
     162     5747308 :               _fsp[ph].mass_fraction[comp].derivatives()[_Zidx];
     163             :       }
     164             : 
     165             :       // If temperature is a PorousFlow variable (nonisothermal case), add derivatives wrt
     166             :       // temperature
     167     2527874 :       if (_dictator.isPorousFlowVariable(_temperature_varnum))
     168             :       {
     169     1105552 :         (*_dporepressure_dvar)[_qp][ph][_Tvar] = _fsp[ph].pressure.derivatives()[_Tidx];
     170     1132866 :         (*_dsaturation_dvar)[_qp][ph][_Tvar] = _fsp[ph].saturation.derivatives()[_Tidx];
     171     1644599 :         (*_dfluid_density_dvar)[_qp][ph][_Tvar] = _fsp[ph].density.derivatives()[_Tidx];
     172     1644599 :         (*_dfluid_viscosity_dvar)[_qp][ph][_Tvar] = _fsp[ph].viscosity.derivatives()[_Tidx];
     173     1644599 :         (*_dfluid_enthalpy_dvar)[_qp][ph][_Tvar] = _fsp[ph].enthalpy.derivatives()[_Tidx];
     174     1078238 :         (*_dfluid_internal_energy_dvar)[_qp][ph][_Tvar] =
     175     1078238 :             _fsp[ph].internal_energy.derivatives()[_Tidx];
     176             : 
     177     3747160 :         for (unsigned int comp = 0; comp < _num_components; ++comp)
     178     2668922 :           (*_dmass_frac_dvar)[_qp][ph][comp][_Tvar] =
     179     2668922 :               _fsp[ph].mass_fraction[comp].derivatives()[_Tidx];
     180             :       }
     181             : 
     182             :       // If Xnacl is a PorousFlow variable, add derivatives wrt Xnacl
     183     2527874 :       if (_dictator.isPorousFlowVariable(_Xnacl_varnum))
     184             :       {
     185      706283 :         (*_dporepressure_dvar)[_qp][ph][_Xvar] = _fsp[ph].pressure.derivatives()[_Xidx];
     186      721006 :         (*_dsaturation_dvar)[_qp][ph][_Xvar] = _fsp[ph].saturation.derivatives()[_Xidx];
     187     1037304 :         (*_dfluid_density_dvar)[_qp][ph][_Xvar] += _fsp[ph].density.derivatives()[_Xidx];
     188     1037304 :         (*_dfluid_viscosity_dvar)[_qp][ph][_Xvar] += _fsp[ph].viscosity.derivatives()[_Xidx];
     189     1037304 :         (*_dfluid_enthalpy_dvar)[_qp][ph][_Xvar] = _fsp[ph].enthalpy.derivatives()[_Xidx];
     190      691560 :         (*_dfluid_internal_energy_dvar)[_qp][ph][_Xvar] =
     191      691560 :             _fsp[ph].internal_energy.derivatives()[_Xidx];
     192             : 
     193     2766240 :         for (unsigned int comp = 0; comp < _num_components; ++comp)
     194     2074680 :           (*_dmass_frac_dvar)[_qp][ph][comp][_Xvar] =
     195     2074680 :               _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     1263937 :   if (!_nodal_material)
     204             :     if constexpr (!is_ad)
     205             :     {
     206             :       // Derivatives of capillary pressure
     207      629815 :       const Real dpc = _pc.dCapillaryPressure(_fsp[_aqueous_phase_number].saturation.value(), _qp);
     208             :       const Real d2pc =
     209      629815 :           _pc.d2CapillaryPressure(_fsp[_aqueous_phase_number].saturation.value(), _qp);
     210             : 
     211             :       // Gradients of saturation and porepressure in all phases
     212      629815 :       (*_grads_qp)[_qp][_gas_phase_number] =
     213      629815 :           (*_dsaturation_dvar)[_qp][_gas_phase_number][_pvar] * _gas_gradp_qp[_qp] +
     214      629815 :           (*_dsaturation_dvar)[_qp][_gas_phase_number][_Zvar[0]] * (*_gradZ_qp[0])[_qp] +
     215      629815 :           (*_dsaturation_dvar)[_qp][_gas_phase_number][_Tvar] * (*_gradT_qp)[_qp];
     216      629815 :       (*_grads_qp)[_qp][_aqueous_phase_number] = -(*_grads_qp)[_qp][_gas_phase_number];
     217             : 
     218      629815 :       (*_gradp_qp)[_qp][_gas_phase_number] = _gas_gradp_qp[_qp];
     219      629815 :       (*_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      629815 :       (*_grad_mass_frac_qp)[_qp][_aqueous_phase_number][_aqueous_fluid_component] =
     224      629815 :           _fsp[_aqueous_phase_number].mass_fraction[_aqueous_fluid_component].derivatives()[_pidx] *
     225             :               _gas_gradp_qp[_qp] +
     226      629815 :           _fsp[_aqueous_phase_number].mass_fraction[_aqueous_fluid_component].derivatives()[_Zidx] *
     227             :               (*_gradZ_qp[0])[_qp] +
     228      629815 :           _fsp[_aqueous_phase_number].mass_fraction[_aqueous_fluid_component].derivatives()[_Tidx] *
     229             :               (*_gradT_qp)[_qp];
     230      629815 :       (*_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      629815 :       (*_grad_mass_frac_qp)[_qp][_gas_phase_number][_aqueous_fluid_component] =
     234             :           _fsp[_gas_phase_number].mass_fraction[_aqueous_fluid_component].derivatives()[_pidx] *
     235             :               _gas_gradp_qp[_qp] +
     236             :           _fsp[_gas_phase_number].mass_fraction[_aqueous_fluid_component].derivatives()[_Zidx] *
     237             :               (*_gradZ_qp[0])[_qp] +
     238      629815 :           _fsp[_gas_phase_number].mass_fraction[_aqueous_fluid_component].derivatives()[_Tidx] *
     239             :               (*_gradT_qp)[_qp];
     240      629815 :       (*_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      629815 :       if (_dictator.isPorousFlowVariable(_gas_porepressure_varnum))
     245             :       {
     246     1889445 :         for (unsigned int ph = 0; ph < _num_phases; ++ph)
     247     1259630 :           (*_dgradp_qp_dgradv)[_qp][ph][_pvar] = 1.0;
     248             : 
     249      629815 :         (*_dgradp_qp_dgradv)[_qp][_aqueous_phase_number][_pvar] +=
     250      629815 :             -dpc * (*_dsaturation_dvar)[_qp][_aqueous_phase_number][_pvar];
     251             : 
     252      629815 :         (*_dgradp_qp_dv)[_qp][_aqueous_phase_number][_pvar] =
     253      629815 :             -d2pc * (*_grads_qp)[_qp][_aqueous_phase_number] *
     254             :             (*_dsaturation_dvar)[_qp][_aqueous_phase_number][_pvar];
     255             :       }
     256             : 
     257      629815 :       if (_dictator.isPorousFlowVariable(_Z_varnum[0]))
     258             :       {
     259      629815 :         (*_dgradp_qp_dgradv)[_qp][_aqueous_phase_number][_Zvar[0]] =
     260      629815 :             -dpc * (*_dsaturation_dvar)[_qp][_aqueous_phase_number][_Zvar[0]];
     261             : 
     262      629815 :         (*_dgradp_qp_dv)[_qp][_aqueous_phase_number][_Zvar[0]] =
     263      629815 :             -d2pc * (*_grads_qp)[_qp][_aqueous_phase_number] *
     264             :             (*_dsaturation_dvar)[_qp][_aqueous_phase_number][_Zvar[0]];
     265             :       }
     266             : 
     267      629815 :       if (_dictator.isPorousFlowVariable(_temperature_varnum))
     268             :       {
     269      268347 :         (*_dgradp_qp_dgradv)[_qp][_aqueous_phase_number][_Tvar] =
     270      268347 :             -dpc * (*_dsaturation_dvar)[_qp][_aqueous_phase_number][_Tvar];
     271             : 
     272      268347 :         (*_dgradp_qp_dv)[_qp][_aqueous_phase_number][_Tvar] =
     273      268347 :             -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      629815 :       if (_dictator.isPorousFlowVariable(_Xnacl_varnum))
     279             :       {
     280      172422 :         (*_dgradp_qp_dgradv)[_qp][_aqueous_phase_number][_Xvar] =
     281      172422 :             -dpc * (*_dsaturation_dvar)[_qp][_aqueous_phase_number][_Xvar];
     282             : 
     283      172422 :         (*_grads_qp)[_qp][_aqueous_phase_number] +=
     284      172422 :             (*_dsaturation_dvar)[_qp][_aqueous_phase_number][_Xvar] * _grad_Xnacl_qp[_qp];
     285             : 
     286      172422 :         (*_grads_qp)[_qp][_gas_phase_number] -=
     287      172422 :             (*_dsaturation_dvar)[_qp][_aqueous_phase_number][_Xvar] * _grad_Xnacl_qp[_qp];
     288             : 
     289      172422 :         (*_gradp_qp)[_qp][_aqueous_phase_number] -=
     290      172422 :             dpc * (*_dsaturation_dvar)[_qp][_aqueous_phase_number][_Xvar] * _grad_Xnacl_qp[_qp];
     291             : 
     292      172422 :         (*_dgradp_qp_dv)[_qp][_aqueous_phase_number][_Xvar] =
     293      172422 :             -d2pc * (*_grads_qp)[_qp][_aqueous_phase_number] *
     294             :             (*_dsaturation_dvar)[_qp][_aqueous_phase_number][_Xvar];
     295             : 
     296      172422 :         (*_grad_mass_frac_qp)[_qp][_aqueous_phase_number][_salt_component] = _grad_Xnacl_qp[_qp];
     297             :         (*_grad_mass_frac_qp)[_qp][_aqueous_phase_number][_aqueous_fluid_component] +=
     298             :             _fsp[_aqueous_phase_number]
     299      172422 :                 .mass_fraction[_aqueous_fluid_component]
     300      172422 :                 .derivatives()[_Xidx] *
     301             :             _grad_Xnacl_qp[_qp];
     302      172422 :         (*_grad_mass_frac_qp)[_qp][_aqueous_phase_number][_gas_fluid_component] -=
     303             :             _fsp[_aqueous_phase_number]
     304             :                 .mass_fraction[_aqueous_fluid_component]
     305      172422 :                 .derivatives()[_Xidx] *
     306             :             _grad_Xnacl_qp[_qp];
     307             :         (*_grad_mass_frac_qp)[_qp][_gas_phase_number][_aqueous_fluid_component] +=
     308      172422 :             _fsp[_gas_phase_number].mass_fraction[_aqueous_fluid_component].derivatives()[_Xidx] *
     309             :             _grad_Xnacl_qp[_qp];
     310             :         (*_grad_mass_frac_qp)[_qp][_gas_phase_number][_gas_fluid_component] -=
     311      172422 :             _fsp[_gas_phase_number].mass_fraction[_aqueous_fluid_component].derivatives()[_Xidx] *
     312             :             _grad_Xnacl_qp[_qp];
     313             :       }
     314             :     }
     315     1505887 : }
     316             : 
     317             : template class PorousFlowFluidStateTempl<false>;
     318             : template class PorousFlowFluidStateTempl<true>;