20 #ifndef GEOS_FINITEVOLUME_SURFACEELEMENTSTENCIL_HPP_
21 #define GEOS_FINITEVOLUME_SURFACEELEMENTSTENCIL_HPP_
47 template<
typename VIEWTYPE >
64 real64 const meanPermCoefficient );
197 {
return m_cellCenterToEdgeCenters.toViewConst(); }
223 real64 m_meanPermCoefficient;
233 virtual void move( LvArray::MemorySpace
const space )
override;
236 localIndex const *
const elementRegionIndices,
237 localIndex const *
const elementSubRegionIndices,
239 real64 const *
const weights,
249 R1Tensor const *
const cellCenterToEdgeCenter,
282 {
return m_cellCenterToEdgeCenters.toViewConst(); }
290 m_meanPermCoefficient = meanPermCoefficient;
299 real64 m_meanPermCoefficient = 1.0;
309 real64 ( & weight )[maxNumConnections][2],
310 real64 ( & dWeight_dVar )[maxNumConnections][2] )
const
320 sumOfTrans += coefficient[er][esr][ei][0][0] *
m_weights[iconn][k];
337 real64 const t0 =
m_weights[iconn][0] * coefficient[er0][esr0][ei0][0][0];
338 real64 const t1 =
m_weights[iconn][1] * coefficient[er1][esr1][ei1][0][0];
340 real64 const harmonicWeight = t0*t1 / sumOfTrans;
341 real64 const arithmeticWeight = 0.25 * (t0+t1);
343 real64 const value = m_meanPermCoefficient * harmonicWeight + (1 - m_meanPermCoefficient) * arithmeticWeight;
345 weight[connectionIndex][0] = value;
346 weight[connectionIndex][1] = -value;
348 real64 const dt0 =
m_weights[iconn][0] * dCoeff_dVar[er0][esr0][ei0][0][0];
349 real64 const dt1 =
m_weights[iconn][1] * dCoeff_dVar[er1][esr1][ei1][0][0];
352 dHarmonic[0] = ( dt0 * t1 * sumOfTrans - dt0 * t0 * t1 ) / ( sumOfTrans * sumOfTrans );
353 dHarmonic[1] = ( t0 * dt1 * sumOfTrans - dt1 * t0 * t1 ) / ( sumOfTrans * sumOfTrans );
356 dArithmetic[0] = 0.25 * dt0;
357 dArithmetic[1] = 0.25 * dt1;
359 dWeight_dVar[connectionIndex][0] = m_meanPermCoefficient * dHarmonic[0] + (1 - m_meanPermCoefficient) * dArithmetic[0];
360 dWeight_dVar[connectionIndex][1] = -( m_meanPermCoefficient * dHarmonic[1] + (1 - m_meanPermCoefficient) * dArithmetic[1] );
371 real64 ( & weight )[maxNumConnections][2],
372 real64 ( & dWeight_dVar )[maxNumConnections][2] )
const
390 real64 const harmonicWeight = t0*t1 / sumOfTrans;
391 real64 const arithmeticWeight = 0.25 * (t0+t1);
393 real64 const value = m_meanPermCoefficient * harmonicWeight + (1 - m_meanPermCoefficient) * arithmeticWeight;
395 weight[connectionIndex][0] = value;
396 weight[connectionIndex][1] = -value;
402 dHarmonic[0] = ( dt0 * t1 * sumOfTrans - dt0 * t0 * t1 ) / ( sumOfTrans * sumOfTrans );
403 dHarmonic[1] = ( t0 * dt1 * sumOfTrans - dt1 * t0 * t1 ) / ( sumOfTrans * sumOfTrans );
406 dArithmetic[0] = 0.25 * dt0;
407 dArithmetic[1] = 0.25 * dt1;
409 dWeight_dVar[connectionIndex][0] = m_meanPermCoefficient * dHarmonic[0] + (1 - m_meanPermCoefficient) * dArithmetic[0];
410 dWeight_dVar[connectionIndex][1] = -( m_meanPermCoefficient * dHarmonic[1] + (1 - m_meanPermCoefficient) * dArithmetic[1] );
426 real64 (& weight)[maxNumConnections][2],
427 real64 (& dWeight_dVar1 )[maxNumConnections][2],
428 real64 (& dWeight_dVar2 )[maxNumConnections][2][3] )
const
437 sumOfTrans += coefficient[er][esr][ei][0][0] *
m_weights[iconn][k];
454 real64 const t0 =
m_weights[iconn][0] * coefficient[er0][esr0][ei0][0][0];
455 real64 const t1 =
m_weights[iconn][1] * coefficient[er1][esr1][ei1][0][0];
457 real64 const harmonicWeight = t0*t1 / sumOfTrans;
458 real64 const arithmeticWeight = 0.25 * (t0+t1);
460 real64 const value = m_meanPermCoefficient * harmonicWeight + (1 - m_meanPermCoefficient) * arithmeticWeight;
462 weight[connectionIndex][0] = value;
463 weight[connectionIndex][1] = -value;
465 real64 const dt0_dvar1 =
m_weights[iconn][0] * dCoeff_dVar1[er0][esr0][ei0][0][0];
466 real64 const dt1_dvar1 =
m_weights[iconn][1] * dCoeff_dVar1[er1][esr1][ei1][0][0];
468 real64 dHarmonic_dvar1[2];
469 dHarmonic_dvar1[0] = ( dt0_dvar1 * t1 * sumOfTrans - dt0_dvar1 * t0 * t1 ) / ( sumOfTrans * sumOfTrans );
470 dHarmonic_dvar1[1] = ( dt0_dvar1 * t1 * sumOfTrans - dt0_dvar1 * t0 * t1 ) / ( sumOfTrans * sumOfTrans );
472 real64 dArithmetic_dvar1[2];
473 dArithmetic_dvar1[0] = 0.25 * dt0_dvar1;
474 dArithmetic_dvar1[1] = 0.25 * dt1_dvar1;
476 dWeight_dVar1[connectionIndex][0] = m_meanPermCoefficient * dHarmonic_dvar1[0] + (1 - m_meanPermCoefficient) * dArithmetic_dvar1[0];
477 dWeight_dVar1[connectionIndex][1] = -( m_meanPermCoefficient * dHarmonic_dvar1[1] + (1 - m_meanPermCoefficient) * dArithmetic_dvar1[1] );
479 real64 const dt0_dvar2 =
m_weights[iconn][0] * dCoeff_dVar2[er0][esr0][ei0][0][0][0];
480 real64 const dt1_dvar2 =
m_weights[iconn][1] * dCoeff_dVar2[er1][esr1][ei1][0][0][0];
482 real64 dHarmonic_dvar2[2];
483 dHarmonic_dvar2[0] = ( dt0_dvar2 * t1 * sumOfTrans - dt0_dvar2 * t0 * t1 ) / ( sumOfTrans * sumOfTrans );
484 dHarmonic_dvar2[1] = ( t0 * dt1_dvar2 * sumOfTrans - dt1_dvar2 * t0 * t1 ) / ( sumOfTrans * sumOfTrans );
486 real64 dArithmetic_dvar2[2];
487 dArithmetic_dvar2[0] = 0.25 * dt0_dvar2;
488 dArithmetic_dvar2[1] = 0.25 * dt1_dvar2;
490 dWeight_dVar2[connectionIndex][0][0] = ( m_meanPermCoefficient * dHarmonic_dvar2[0] + (1 - m_meanPermCoefficient) * dArithmetic_dvar2[0] );
491 dWeight_dVar2[connectionIndex][1][0] = -( m_meanPermCoefficient * dHarmonic_dvar2[1] + (1 - m_meanPermCoefficient) * dArithmetic_dvar2[1] );
505 real64 (& weight)[maxNumConnections][2] )
const
515 real64 const mult = ( LvArray::math::abs( LvArray::tensorOps::AiBi< 3 >( m_cellCenterToEdgeCenters[iconn][k], gravityVector ) ) >
MULTIPLIER_THRESHOLD )
516 ? coefficientMultiplier[er][esr][ei][0][1] : coefficientMultiplier[er][esr][ei][0][0];
518 sumOfTrans += mult * coefficient[er][esr][ei][0][0] *
m_weights[iconn][k];
536 real64 const mult0 = ( LvArray::math::abs( LvArray::tensorOps::AiBi< 3 >( m_cellCenterToEdgeCenters[iconn][k[0]], gravityVector ) ) >
MULTIPLIER_THRESHOLD )
537 ? coefficientMultiplier[er0][esr0][ei0][0][1] : coefficientMultiplier[er0][esr0][ei0][0][0];
538 real64 const mult1 = ( LvArray::math::abs( LvArray::tensorOps::AiBi< 3 >( m_cellCenterToEdgeCenters[iconn][k[1]], gravityVector ) ) >
MULTIPLIER_THRESHOLD )
539 ? coefficientMultiplier[er1][esr1][ei1][0][1] : coefficientMultiplier[er1][esr1][ei1][0][0];
541 real64 const t0 = mult0 *
m_weights[iconn][0] * coefficient[er0][esr0][ei0][0][0];
542 real64 const t1 = mult1 *
m_weights[iconn][1] * coefficient[er1][esr1][ei1][0][0];
544 real64 const harmonicWeight = t0*t1 / sumOfTrans;
545 real64 const arithmeticWeight = 0.25 * (t0+t1);
547 real64 const value = m_meanPermCoefficient * harmonicWeight + (1 - m_meanPermCoefficient) * arithmeticWeight;
549 weight[connectionIndex][0] = value;
550 weight[connectionIndex][1] = -value;
565 real64 ( & weight1 )[maxNumPointsInFlux],
566 real64 ( & weight2 )[maxNumPointsInFlux],
567 real64 ( & geometricWeight )[maxNumPointsInFlux] )
const
569 real64 sumOfGeometricWeights = 0.0;
577 real64 const cellToEdgeDistance = LvArray::tensorOps::l2Norm< 3 >( m_cellCenterToEdgeCenters[iconn][k] );
579 real64 const edgeToFaceDownDistance = -LvArray::tensorOps::AiBi< 3 >( m_cellCenterToEdgeCenters[iconn][k], unitGravityVector )
580 * edgeLength / cellToEdgeDistance;
583 ? coefficient1Multiplier[er][esr][ei][0][1] : coefficient1Multiplier[er][esr][ei][0][0];
585 weight1[k] = mult * coefficient1[er][esr][ei][0][0] *
m_weights[iconn][k];
586 weight2[k] = coefficient2[er][esr][ei] * edgeToFaceDownDistance;
588 geometricWeight[k] =
m_weights[iconn][k] / 12.0;
589 sumOfGeometricWeights += geometricWeight[k];
594 geometricWeight[k] /= sumOfGeometricWeights;
#define GEOS_HOST_DEVICE
Marks a host-device function.
#define GEOS_UNUSED_VAR(...)
Mark an unused variable and silence compiler warnings.
#define GEOS_FORCE_INLINE
Marks a function or lambda for inlining.
typename ElementViewAccessor< VIEWTYPE >::NestedViewTypeConst ElementViewConst
The ElementViewAccessor at the ElementRegionManager level is the type resulting from ElementViewAcces...
Provides management of the interior stencil points when using Two-Point flux approximation.
TRAITS::IndexContainerType m_elementRegionIndices
The container for the element region indices for each point in each stencil.
TRAITS::IndexContainerViewConstType m_elementRegionIndices
The container for the element region indices for each point in each stencil.
TRAITS::IndexContainerViewConstType m_elementSubRegionIndices
The container for the element sub region indices for each point in each stencil.
TRAITS::WeightContainerViewType m_weights
The container for the weights for each point in each stencil.
TRAITS::IndexContainerViewConstType m_elementIndices
The container for the element indices for each point in each stencil.
Provides management of the interior stencil points for a face elements when using Two-Point flux appr...
KernelWrapper createKernelWrapper() const
Create an update kernel wrapper.
virtual localIndex size() const override
Return the stencil size.
void add(localIndex const numPts, R1Tensor const *const cellCenterToEdgeCenter, localIndex const connectorIndex)
Add an entry to the stencil.
void setMeanPermCoefficient(real64 const &meanPermCoefficient)
sets the value of the mean perm conefficient
localIndex stencilSize(localIndex index) const
Give the number of stencil entries for the provided index.
virtual void move(LvArray::MemorySpace const space) override
Move the data arrays associated with the stencil to a specified memory space.
ArrayOfArraysView< R1Tensor const > getCellCenterToEdgeCenters() const
Give the array of vectors pointing from the cell center to the edge center.
virtual void add(localIndex const numPts, localIndex const *const elementRegionIndices, localIndex const *const elementSubRegionIndices, localIndex const *const elementIndices, real64 const *const weights, localIndex const connectorIndex) override
Add an entry to the stencil.
Provides access to the SurfaceElementStencil that may be called from a kernel function.
GEOS_HOST_DEVICE void computeStabilizationWeights(localIndex iconn, real64(&stabilizationWeight)[maxNumConnections][2]) const
Compute the stabilization weights.
GEOS_HOST_DEVICE GEOS_FORCE_INLINE localIndex size() const
Give the number of stencil entries.
ElementRegionManager::ElementViewConst< VIEWTYPE > CoefficientAccessor
Coefficient view accessory type.
GEOS_HOST_DEVICE GEOS_FORCE_INLINE localIndex stencilSize(localIndex index) const
Give the number of stencil entries for the provided index.
GEOS_HOST_DEVICE void addHydraulicApertureContribution(localIndex const iconn, ElementRegionManager::ElementViewConst< arrayView1d< real64 const > > hydraulicAperture) const
Add the contribution of the aperture to the weight in the stencil (done after aperture update)
static constexpr real64 MULTIPLIER_THRESHOLD
Threshold for the application of the permeability multiplier.
ArrayOfArraysView< R1Tensor const > getCellCenterToEdgeCenters() const
Accessor to the CellCenterToEdgeCenter vector.
GEOS_HOST_DEVICE void computeWeights(localIndex iconn, CoefficientAccessor< arrayView3d< real64 const > > const &coefficient, CoefficientAccessor< arrayView3d< real64 const > > const &dCoeff_dVar, real64(&weight)[maxNumConnections][2], real64(&dWeight_dVar)[maxNumConnections][2]) const
Compute weights and derivatives w.r.t to one variable.
GEOS_HOST_DEVICE GEOS_FORCE_INLINE localIndex numPointsInFlux(localIndex index) const
Give the number of points between which the flux is.
SurfaceElementStencilWrapper(IndexContainerType const &elementRegionIndices, IndexContainerType const &elementSubRegionIndices, IndexContainerType const &elementIndices, WeightContainerType const &weights, ArrayOfArrays< R1Tensor > const &cellCenterToEdgeCenters, real64 const meanPermCoefficient)
Constructor.
GEOS_HOST_DEVICE void removeHydraulicApertureContribution(localIndex const iconn, ElementRegionManager::ElementViewConst< arrayView1d< real64 const > > hydraulicAperture) const
Remove the contribution of the aperture from the weight in the stencil (done before aperture update)
ArrayView< T, 1 > arrayView1d
Alias for 1D array view.
LvArray::ArrayOfArraysView< T, INDEX_TYPE const, CONST_SIZES, LvArray::ChaiBuffer > ArrayOfArraysView
View of array of variable-sized arrays. See LvArray::ArrayOfArraysView for details.
double real64
64-bit floating point type.
GEOS_LOCALINDEX_TYPE localIndex
Local index type (for indexing objects within an MPI partition).
ArrayView< T, 4, USD > arrayView4d
Alias for 4D array view.
LvArray::ArrayOfArrays< T, INDEX_TYPE, LvArray::ChaiBuffer > ArrayOfArrays
Array of variable-sized arrays. See LvArray::ArrayOfArrays for details.
ArrayView< T, 3, USD > arrayView3d
Alias for 3D array view.
A collection of properties of a stencil type.
CONTAINER< localIndex > IndexContainerType
The array type that will be used to store the indices of the stencil contributors.
CONTAINER< real64 > WeightContainerType
The array type that is used to store the weights of the stencil contributors.
static constexpr localIndex maxNumPointsInFlux
Maximum number of points the flux.
static constexpr localIndex maxNumConnections
Maximum number of connections in a stencil.