FluxComputeKernelBase.hpp

Go to the documentation of this file.

/*
 * ------------------------------------------------------------------------------------------------------------
 * SPDX-License-Identifier: LGPL-2.1-only
 *
 * Copyright (c) 2016-2024 Lawrence Livermore National Security LLC
 * Copyright (c) 2018-2024 TotalEnergies
 * Copyright (c) 2018-2024 The Board of Trustees of the Leland Stanford Junior University
 * Copyright (c) 2023-2024 Chevron
 * Copyright (c) 2019-     GEOS/GEOSX Contributors
 * All rights reserved
 *
 * See top level LICENSE, COPYRIGHT, CONTRIBUTORS, NOTICE, and ACKNOWLEDGEMENTS files for details.
 * ------------------------------------------------------------------------------------------------------------
 */
 
#ifndef GEOS_PHYSICSSOLVERS_FLUIDFLOW_COMPOSITIONAL_FLUXCOMPUTEKERNELBASE_HPP
#define GEOS_PHYSICSSOLVERS_FLUIDFLOW_COMPOSITIONAL_FLUXCOMPUTEKERNELBASE_HPP
 
#include "common/DataLayouts.hpp"
#include "common/DataTypes.hpp"
#include "common/GEOS_RAJA_Interface.hpp"
#include "constitutive/capillaryPressure/CapillaryPressureFields.hpp"
#include "constitutive/capillaryPressure/CapillaryPressureBase.hpp"
#include "constitutive/fluid/multifluid/MultiFluidBase.hpp"
#include "constitutive/fluid/multifluid/MultiFluidFields.hpp"
#include "constitutive/permeability/PermeabilityBase.hpp"
#include "constitutive/permeability/PermeabilityFields.hpp"
#include "mesh/ElementRegionManager.hpp"
#include "physicsSolvers/fluidFlow/FlowSolverBaseFields.hpp"
#include "physicsSolvers/fluidFlow/CompositionalMultiphaseBaseFields.hpp"
#include "physicsSolvers/fluidFlow/StencilAccessors.hpp"
 
 
namespace geos
{
 
namespace isothermalCompositionalMultiphaseFVMKernels
{
 
enum class KernelFlags
{
  CapPressure = 1 << 0, // 1
  Diffusion = 1 << 1, // 2
  Dispersion = 1 << 2, // 4
  TotalMassEquation = 1 << 3, // 8
  CheckPhasePresenceInGravity = 1 << 4, // 16
  C1PPU = 1 << 5, // 32
  IHU = 1 << 6, // 64
  HU2PH = 1 << 7 // 128
};
 
/******************************** FluxComputeKernelBase ********************************/
 
class FluxComputeKernelBase
{
public:
 
  template< typename VIEWTYPE >
  using ElementViewConst = ElementRegionManager::ElementViewConst< VIEWTYPE >;
 
  using DofNumberAccessor = ElementRegionManager::ElementViewAccessor< arrayView1d< globalIndex const > >;
 
  using CompFlowAccessors =
    StencilAccessors< fields::ghostRank,
                      fields::flow::gravityCoefficient,
                      fields::flow::pressure,
                      fields::flow::dGlobalCompFraction_dGlobalCompDensity,
                      fields::flow::phaseVolumeFraction,
                      fields::flow::dPhaseVolumeFraction,
                      fields::flow::phaseMobility,
                      fields::flow::dPhaseMobility >;
  using MultiFluidAccessors =
    StencilMaterialAccessors< constitutive::MultiFluidBase,
                              fields::multifluid::phaseDensity,
                              fields::multifluid::dPhaseDensity,
                              fields::multifluid::phaseMassDensity,
                              fields::multifluid::dPhaseMassDensity,
                              fields::multifluid::phaseCompFraction,
                              fields::multifluid::dPhaseCompFraction >;
 
  using CapPressureAccessors =
    StencilMaterialAccessors< constitutive::CapillaryPressureBase,
                              fields::cappres::phaseCapPressure,
                              fields::cappres::dPhaseCapPressure_dPhaseVolFraction >;
 
  using PermeabilityAccessors =
    StencilMaterialAccessors< constitutive::PermeabilityBase,
                              fields::permeability::permeability,
                              fields::permeability::dPerm_dPressure >;
 
  FluxComputeKernelBase( integer const numPhases,
                         globalIndex const rankOffset,
                         DofNumberAccessor const & dofNumberAccessor,
                         CompFlowAccessors const & compFlowAccessors,
                         MultiFluidAccessors const & multiFluidAccessors,
                         real64 const dt,
                         CRSMatrixView< real64, globalIndex const > const & localMatrix,
                         arrayView1d< real64 > const & localRhs,
                         BitFlags< KernelFlags > kernelFlags );
 
protected:
 
  integer const m_numPhases;
 
  globalIndex const m_rankOffset;
 
  real64 const m_dt;
 
  ElementViewConst< arrayView1d< globalIndex const > > const m_dofNumber;
 
  ElementViewConst< arrayView1d< integer const > > const m_ghostRank;
  ElementViewConst< arrayView1d< real64 const > > const m_gravCoef;
 
  // Primary and secondary variables
 
  ElementViewConst< arrayView1d< real64 const > > const m_pres;
 
  ElementViewConst< arrayView2d< real64 const, compflow::USD_PHASE > > const m_phaseVolFrac;
  ElementViewConst< arrayView3d< real64 const, compflow::USD_PHASE_DC > > const m_dPhaseVolFrac;
 
  ElementViewConst< arrayView3d< real64 const, compflow::USD_COMP_DC > > const m_dCompFrac_dCompDens;
 
  ElementViewConst< arrayView4d< real64 const, constitutive::multifluid::USD_PHASE_COMP > > const m_phaseCompFrac;
  ElementViewConst< arrayView5d< real64 const, constitutive::multifluid::USD_PHASE_COMP_DC > > const m_dPhaseCompFrac;
 
  // Residual and jacobian
 
  CRSMatrixView< real64, globalIndex const > const m_localMatrix;
  arrayView1d< real64 > const m_localRhs;
 
  BitFlags< KernelFlags > const m_kernelFlags;
};
 
namespace helpers
{
template< typename VIEWTYPE >
using ElementViewConst = ElementRegionManager::ElementViewConst< VIEWTYPE >;
 
template< localIndex numComp, localIndex numFluxSupportPoints >
GEOS_HOST_DEVICE
static void calculateMeanDensity( localIndex const ip,
                                  localIndex const (&seri)[numFluxSupportPoints],
                                  localIndex const (&sesri)[numFluxSupportPoints],
                                  localIndex const (&sei)[numFluxSupportPoints],
                                  integer const checkPhasePresenceInGravity,
                                  ElementViewConst< arrayView2d< real64 const, compflow::USD_PHASE > > const & phaseVolFrac,
                                  ElementViewConst< arrayView3d< real64 const, compflow::USD_COMP_DC > > const & dCompFrac_dCompDens,
                                  ElementViewConst< arrayView3d< real64 const, constitutive::multifluid::USD_PHASE > > const & phaseMassDens,
                                  ElementViewConst< arrayView4d< real64 const, constitutive::multifluid::USD_PHASE_DC > > const & dPhaseMassDens,
                                  real64 & densMean, real64 (& dDensMean_dPres)[numFluxSupportPoints], real64 (& dDensMean_dComp)[numFluxSupportPoints][numComp] )
{
  using Deriv = constitutive::multifluid::DerivativeOffset;
 
  densMean = 0;
  integer denom = 0;
  real64 dDens_dC[numComp]{};
  for( localIndex i = 0; i < numFluxSupportPoints; ++i )
  {
    localIndex const er = seri[i];
    localIndex const esr = sesri[i];
    localIndex const ei = sei[i];
 
    bool const phaseExists = (phaseVolFrac[er][esr][ei][ip] > 0);
    if( checkPhasePresenceInGravity && !phaseExists )
    {
      dDensMean_dPres[i] = 0.0;
      for( localIndex jc = 0; jc < numComp; ++jc )
      {
        dDensMean_dComp[i][jc] = 0.0;
      }
      continue;
    }
 
    // density
    real64 const density = phaseMassDens[er][esr][ei][0][ip];
    real64 const dDens_dPres = dPhaseMassDens[er][esr][ei][0][ip][Deriv::dP];
 
    applyChainRule( numComp,
                    dCompFrac_dCompDens[er][esr][ei],
                    dPhaseMassDens[er][esr][ei][0][ip],
                    dDens_dC,
                    Deriv::dC );
 
    // average density and derivatives
    densMean += density;
    dDensMean_dPres[i] = dDens_dPres;
    for( localIndex jc = 0; jc < numComp; ++jc )
    {
      dDensMean_dComp[i][jc] = dDens_dC[jc];
    }
    denom++;
  }
  if( denom > 1 )
  {
    densMean /= denom;
    for( localIndex i = 0; i < numFluxSupportPoints; ++i )
    {
      dDensMean_dPres[i] /= denom;
      for( integer jc = 0; jc < numComp; ++jc )
      {
        dDensMean_dComp[i][jc] /= denom;
      }
    }
  }
}
 
}
 
} // namespace isothermalCompositionalMultiphaseFVMKernels
 
} // namespace geos
 
 
#endif //GEOS_PHYSICSSOLVERS_FLUIDFLOW_COMPOSITIONAL_FLUXCOMPUTEKERNELBASE_HPP