20 #ifndef GEOS_PHYSICSSOLVERS_FLUIDFLOW_COMPOSITIONAL_THERMALFLUXCOMPUTEKERNEL_HPP
21 #define GEOS_PHYSICSSOLVERS_FLUIDFLOW_COMPOSITIONAL_THERMALFLUXCOMPUTEKERNEL_HPP
24 #include "constitutive/fluid/multifluid/MultiFluidBase.hpp"
25 #include "constitutive/fluid/multifluid/MultiFluidFields.hpp"
26 #include "constitutive/thermalConductivity/MultiPhaseThermalConductivityBase.hpp"
27 #include "constitutive/thermalConductivity/ThermalConductivityFields.hpp"
32 namespace thermalCompositionalMultiphaseFVMKernels
44 template<
integer NUM_COMP,
integer NUM_DOF,
typename STENCILWRAPPER >
55 template<
typename VIEWTYPE >
59 using DofNumberAccessor = AbstractBase::DofNumberAccessor;
69 using AbstractBase::m_gravCoef;
71 using AbstractBase::m_dPhaseVolFrac;
73 using AbstractBase::m_dPhaseCompFrac;
85 using Base::m_dPhaseMob;
86 using Base::m_dPhaseMassDens;
87 using Base::m_dPhaseCapPressure_dPhaseVolFrac;
98 fields::multifluid::phaseEnthalpy,
99 fields::multifluid::dPhaseEnthalpy >;
103 fields::thermalconductivity::effectiveConductivity >;
126 STENCILWRAPPER
const & stencilWrapper,
127 DofNumberAccessor
const & dofNumberAccessor,
138 BitFlags< isothermalCompositionalMultiphaseFVMKernels::KernelFlags > kernelFlags )
145 capPressureAccessors,
146 permeabilityAccessors,
151 m_temp( thermalCompFlowAccessors.get( fields::flow::temperature {} ) ),
152 m_phaseEnthalpy( thermalMultiFluidAccessors.get( fields::multifluid::phaseEnthalpy {} ) ),
153 m_dPhaseEnthalpy( thermalMultiFluidAccessors.get( fields::multifluid::dPhaseEnthalpy {} ) ),
154 m_thermalConductivity( thermalConductivityAccessors.get( fields::thermalconductivity::effectiveConductivity {} ) )
189 using Deriv = constitutive::multifluid::DerivativeOffset;
211 real64 const (&dPhaseFlux_dP)[2],
223 real64 convectiveEnergyFlux = 0.0;
236 bool const phaseExists = (
m_phaseVolFrac[er_up][esr_up][ei_up][ip] > 0);
242 dDensMean_dT[i] = m_dPhaseMassDens[er][esr][ei][0][ip][Deriv::dT];
249 dDensMean_dT[i] /= denom;
267 real64 dCapPressure_dT = 0.0;
268 if( AbstractBase::m_kernelFlags.isSet( isothermalCompositionalMultiphaseFVMKernels::KernelFlags::CapPressure ) )
272 real64 const dCapPressure_dS = m_dPhaseCapPressure_dPhaseVolFrac[er][esr][ei][0][ip][jp];
273 dCapPressure_dT += dCapPressure_dS * m_dPhaseVolFrac[er][esr][ei][jp][Deriv::dT];
278 dPresGrad_dT[i] -= trans[i] * dCapPressure_dT;
279 real64 const gravD = trans[i] * m_gravCoef[er][esr][ei];
284 dGravHead_dT[j] += dDensMean_dT[j] * gravD;
299 dPhaseFlux_dT[ke] += dPresGrad_dT[ke];
303 dPhaseFlux_dT[ke] -= dGravHead_dT[ke];
307 dPhaseFlux_dT[ke] *=
m_phaseMob[er_up][esr_up][ei_up][ip];
309 dPhaseFlux_dT[k_up] += m_dPhaseMob[er_up][esr_up][ei_up][ip][Deriv::dT] * potGrad;
314 arraySlice1d<
real64 const, constitutive::multifluid::USD_PHASE_COMP - 3 > phaseCompFracSub =
316 arraySlice2d<
real64 const, constitutive::multifluid::USD_PHASE_COMP_DC - 3 > dPhaseCompFracSub =
317 m_dPhaseCompFrac[er_up][esr_up][ei_up][0][ip];
321 real64 const ycp = phaseCompFracSub[ic];
324 dCompFlux_dT[ke][ic] += dPhaseFlux_dT[ke] * ycp;
326 dCompFlux_dT[k_up][ic] += phaseFlux * dPhaseCompFracSub[ic][Deriv::dT];
344 convectiveEnergyFlux += phaseFlux * enthalpy;
348 dConvectiveEnergyFlux_dP[ke] += dPhaseFlux_dP[ke] * enthalpy;
349 dConvectiveEnergyFlux_dT[ke] += dPhaseFlux_dT[ke] * enthalpy;
353 dConvectiveEnergyFlux_dC[ke][jc] += dPhaseFlux_dC[ke][jc] * enthalpy;
357 dConvectiveEnergyFlux_dP[k_up] += phaseFlux * m_dPhaseEnthalpy[er_up][esr_up][ei_up][0][ip][Deriv::dP];
358 dConvectiveEnergyFlux_dT[k_up] += phaseFlux * m_dPhaseEnthalpy[er_up][esr_up][ei_up][0][ip][Deriv::dT];
363 m_dPhaseEnthalpy[er_up][esr_up][ei_up][0][ip],
368 dConvectiveEnergyFlux_dC[k_up][jc] += phaseFlux * dProp_dC[jc];
374 stack.
localFlux[localRowIndexEnergy0] +=
m_dt * convectiveEnergyFlux;
375 stack.
localFlux[localRowIndexEnergy1] -=
m_dt * convectiveEnergyFlux;
380 stack.
localFluxJacobian[localRowIndexEnergy0][localDofIndexPres] +=
m_dt * dConvectiveEnergyFlux_dP[ke];
381 stack.
localFluxJacobian[localRowIndexEnergy1][localDofIndexPres] -=
m_dt * dConvectiveEnergyFlux_dP[ke];
382 integer const localDofIndexTemp = localDofIndexPres +
numDof - 1;
383 stack.
localFluxJacobian[localRowIndexEnergy0][localDofIndexTemp] +=
m_dt * dConvectiveEnergyFlux_dT[ke];
384 stack.
localFluxJacobian[localRowIndexEnergy1][localDofIndexTemp] -=
m_dt * dConvectiveEnergyFlux_dT[ke];
388 integer const localDofIndexComp = localDofIndexPres + jc + 1;
389 stack.
localFluxJacobian[localRowIndexEnergy0][localDofIndexComp] +=
m_dt * dConvectiveEnergyFlux_dC[ke][jc];
390 stack.
localFluxJacobian[localRowIndexEnergy1][localDofIndexComp] -=
m_dt * dConvectiveEnergyFlux_dC[ke][jc];
407 stack.thermalTransmissibility,
419 real64 const thermalTrans[2] = { stack.thermalTransmissibility[connectionIndex][0], stack.thermalTransmissibility[connectionIndex][1] };
421 localIndex const sesri[2] = {m_sesri( iconn, k[0] ), m_sesri( iconn, k[1] )};
422 localIndex const sei[2] = {m_sei( iconn, k[0] ), m_sei( iconn, k[1] )};
424 real64 conductiveEnergyFlux = 0.0;
434 conductiveEnergyFlux += thermalTrans[ke] *
m_temp[er][esr][ei];
435 dConductiveEnergyFlux_dT[ke] += thermalTrans[ke];
441 stack.
localFlux[localRowIndexEnergy0] +=
m_dt * conductiveEnergyFlux;
442 stack.
localFlux[localRowIndexEnergy1] -=
m_dt * conductiveEnergyFlux;
447 stack.
localFluxJacobian[localRowIndexEnergy0][localDofIndexTemp] +=
m_dt * dConductiveEnergyFlux_dT[ke];
448 stack.
localFluxJacobian[localRowIndexEnergy1][localDofIndexTemp] -=
m_dt * dConductiveEnergyFlux_dT[ke];
518 template<
typename POLICY,
typename STENCILWRAPPER >
523 string const & dofKey,
525 integer const useTotalMassEquation,
526 string const & solverName,
528 STENCILWRAPPER
const & stencilWrapper,
533 isothermalCompositionalMultiphaseBaseKernels::
534 internal::kernelLaunchSelectorCompSwitch( numComps, [&](
auto NC )
536 integer constexpr NUM_COMP = NC();
537 integer constexpr NUM_DOF = NC() + 2;
541 dofNumberAccessor.setName( solverName +
"/accessors/" + dofKey );
543 BitFlags< isothermalCompositionalMultiphaseFVMKernels::KernelFlags > kernelFlags;
545 kernelFlags.set( isothermalCompositionalMultiphaseFVMKernels::KernelFlags::CapPressure );
546 if( useTotalMassEquation )
547 kernelFlags.set( isothermalCompositionalMultiphaseFVMKernels::KernelFlags::TotalMassEquation );
550 typename KernelType::CompFlowAccessors compFlowAccessors( elemManager, solverName );
551 typename KernelType::ThermalCompFlowAccessors thermalCompFlowAccessors( elemManager, solverName );
552 typename KernelType::MultiFluidAccessors multiFluidAccessors( elemManager, solverName );
553 typename KernelType::ThermalMultiFluidAccessors thermalMultiFluidAccessors( elemManager, solverName );
554 typename KernelType::CapPressureAccessors capPressureAccessors( elemManager, solverName );
555 typename KernelType::PermeabilityAccessors permeabilityAccessors( elemManager, solverName );
556 typename KernelType::ThermalConductivityAccessors thermalConductivityAccessors( elemManager, solverName );
558 KernelType kernel( numPhases, rankOffset, stencilWrapper, dofNumberAccessor,
559 compFlowAccessors, thermalCompFlowAccessors, multiFluidAccessors, thermalMultiFluidAccessors,
560 capPressureAccessors, permeabilityAccessors, thermalConductivityAccessors,
561 dt, localMatrix, localRhs, kernelFlags );
562 KernelType::template launch< POLICY >( stencilWrapper.size(), kernel );
#define GEOS_HOST_DEVICE
Marks a host-device function.
The ElementRegionManager class provides an interface to ObjectManagerBase in order to manage ElementR...
ElementViewAccessor< ArrayView< T const, NDIM, getUSD< PERM > > > constructArrayViewAccessor(string const &name, string const &neighborName=string()) const
This is a function to construct a ElementViewAccessor to access array data registered on the mesh.
array1d< array1d< VIEWTYPE > > ElementViewAccessor
The ElementViewAccessor at the ElementRegionManager level is an array of array of VIEWTYPE.
typename ElementViewAccessor< VIEWTYPE >::NestedViewTypeConst ElementViewConst
The ElementViewAccessor at the ElementRegionManager level is the type resulting from ElementViewAcces...
A struct to automatically construct and store element view accessors.
A struct to automatically construct and store element view accessors.
Base class for FluxComputeKernel that holds all data not dependent on template parameters (like stenc...
ElementViewConst< arrayView2d< real64 const, compflow::USD_PHASE > > const m_phaseVolFrac
Views on phase volume fractions.
real64 const m_dt
Time step size.
ElementViewConst< arrayView1d< globalIndex const > > const m_dofNumber
Views on dof numbers.
arrayView1d< real64 > const m_localRhs
View on the local RHS.
ElementViewConst< arrayView4d< real64 const, constitutive::multifluid::USD_PHASE_COMP > > const m_phaseCompFrac
Views on phase component fractions.
globalIndex const m_rankOffset
Offset for my MPI rank.
CRSMatrixView< real64, globalIndex const > const m_localMatrix
View on the local CRS matrix.
ElementViewConst< arrayView3d< real64 const, compflow::USD_COMP_DC > > const m_dCompFrac_dCompDens
Views on derivatives of comp fractions.
integer const m_numPhases
Number of fluid phases.
Define the interface for the assembly kernel in charge of flux terms.
static constexpr localIndex maxStencilSize
Maximum number of points in the stencil.
GEOS_HOST_DEVICE void computeFlux(localIndex const iconn, StackVariables &stack, FUNC &&compFluxKernelOp=NoOpFunc{}) const
Compute the local flux contributions to the residual and Jacobian.
GEOS_HOST_DEVICE void complete(localIndex const iconn, StackVariables &stack, FUNC &&assemblyKernelOp=NoOpFunc{}) const
Performs the complete phase for the kernel.
static constexpr integer numFluxSupportPoints
Number of flux support points (hard-coded for TFPA)
static constexpr integer numEqn
Compute time value for the number of equations (all of them, except the volume balance equation)
static constexpr integer numComp
Compile time value for the number of components.
ElementViewConst< arrayView2d< real64 const, compflow::USD_PHASE > > const m_phaseMob
Views on phase mobilities.
STENCILWRAPPER::IndexContainerViewConstType const m_seri
Connection to element maps.
static constexpr localIndex maxNumElems
Maximum number of elements at the face.
static constexpr integer numDof
Compute time value for the number of degrees of freedom.
STENCILWRAPPER const m_stencilWrapper
Reference to the stencil wrapper.
static constexpr localIndex maxNumConns
Maximum number of connections at the face.
static void createAndLaunch(integer const numComps, integer const numPhases, globalIndex const rankOffset, string const &dofKey, integer const hasCapPressure, integer const useTotalMassEquation, string const &solverName, ElementRegionManager const &elemManager, STENCILWRAPPER const &stencilWrapper, real64 const dt, CRSMatrixView< real64, globalIndex const > const &localMatrix, arrayView1d< real64 > const &localRhs)
Create a new kernel and launch.
Define the interface for the assembly kernel in charge of flux terms.
ElementViewConst< arrayView3d< real64 const, constitutive::multifluid::USD_PHASE > > const m_phaseEnthalpy
Views on phase enthalpies.
ElementViewConst< arrayView2d< real64 const, compflow::USD_PHASE > > const m_phaseVolFrac
Views on phase volume fractions.
real64 const m_dt
Time step size.
static constexpr integer numFluxSupportPoints
Number of flux support points (hard-coded for TFPA)
ElementViewConst< arrayView3d< real64 const > > const m_thermalConductivity
View on thermal conductivity.
static constexpr integer numEqn
Compute time value for the number of equations (all of them, except the volume balance equation)
static constexpr integer numComp
Compile time value for the number of components.
ElementViewConst< arrayView2d< real64 const, compflow::USD_PHASE > > const m_phaseMob
Views on phase mobilities.
GEOS_HOST_DEVICE void computeFlux(localIndex const iconn, StackVariables &stack) const
Compute the local flux contributions to the residual and Jacobian.
GEOS_HOST_DEVICE void complete(localIndex const iconn, StackVariables &stack) const
Performs the complete phase for the kernel.
ElementViewConst< arrayView4d< real64 const, constitutive::multifluid::USD_PHASE_COMP > > const m_phaseCompFrac
Views on phase component fractions.
ElementRegionManager::ElementViewConst< VIEWTYPE > ElementViewConst
The type for element-based data. Consists entirely of ArrayView's.
STENCILWRAPPER::IndexContainerViewConstType const m_seri
Connection to element maps.
static constexpr integer numDof
Compute time value for the number of degrees of freedom.
ElementViewConst< arrayView3d< real64 const, compflow::USD_COMP_DC > > const m_dCompFrac_dCompDens
Views on derivatives of comp fractions.
STENCILWRAPPER const m_stencilWrapper
Reference to the stencil wrapper.
ElementViewConst< arrayView1d< real64 const > > const m_temp
Views on temperature.
integer const m_numPhases
Number of fluid phases.
static constexpr localIndex maxNumConns
Maximum number of connections at the face.
FluxComputeKernel(integer const numPhases, globalIndex const rankOffset, STENCILWRAPPER const &stencilWrapper, DofNumberAccessor const &dofNumberAccessor, CompFlowAccessors const &compFlowAccessors, ThermalCompFlowAccessors const &thermalCompFlowAccessors, MultiFluidAccessors const &multiFluidAccessors, ThermalMultiFluidAccessors const &thermalMultiFluidAccessors, CapPressureAccessors const &capPressureAccessors, PermeabilityAccessors const &permeabilityAccessors, ThermalConductivityAccessors const &thermalConductivityAccessors, real64 const dt, CRSMatrixView< real64, globalIndex const > const &localMatrix, arrayView1d< real64 > const &localRhs, BitFlags< isothermalCompositionalMultiphaseFVMKernels::KernelFlags > kernelFlags)
Constructor for the kernel interface.
ArrayView< T, 1 > arrayView1d
Alias for 1D array view.
LvArray::CRSMatrixView< T, COL_INDEX, localIndex const, LvArray::ChaiBuffer > CRSMatrixView
Alias for CRS Matrix View.
GEOS_GLOBALINDEX_TYPE globalIndex
Global index type (for indexing objects across MPI partitions).
ArraySlice< T, 2, USD > arraySlice2d
Alias for 2D array slice.
double real64
64-bit floating point type.
GEOS_LOCALINDEX_TYPE localIndex
Local index type (for indexing objects within an MPI partition).
ArraySlice< T, 1, USD > arraySlice1d
Alias for 1D array slice.
std::int32_t integer
Signed integer type.
Kernel variables (dof numbers, jacobian and residual) located on the stack.
real64 dTrans_dPres[maxNumConns][numFluxSupportPoints]
Derivatives of transmissibility with respect to pressure.
stackArray1d< globalIndex, maxNumElems *numDof > dofColIndices
Indices of the matrix rows/columns corresponding to the dofs in this face.
stackArray2d< real64, maxNumElems *numEqn *maxStencilSize *numDof > localFluxJacobian
Storage for the face local Jacobian matrix.
stackArray1d< real64, maxNumElems *numEqn > localFlux
Storage for the face local residual vector (all equations except volume balance)
localIndex const numConnectedElems
Number of elements connected at a given connection.
real64 transmissibility[maxNumConns][numFluxSupportPoints]
Transmissibility.
localIndex const stencilSize
Stencil size for a given connection.