[refactor] Remove unnecessary copy in axisCentroid and simply code.

We stored 3 copies of each cell centroid in axisCentroid. This seemed
to be a waste of memory and also made the function distanceVector_ hard
to understand.

With this change we omit this copy of information and simplify distanceVector_
This commit is contained in:
Markus Blatt 2024-06-27 15:15:33 +02:00
parent 54d303ae5d
commit d7c869d01a
2 changed files with 17 additions and 62 deletions

View File

@ -257,10 +257,8 @@ protected:
const DimVector& distance, const DimVector& distance,
const Scalar& poro) const; const Scalar& poro) const;
DimVector distanceVector_(const DimVector& center, DimVector distanceVector_(const DimVector& faceCenter,
int faceIdx, // in the reference element that contains the intersection const std::array<double,dimWorld>& cellCenter) const;
unsigned elemIdx,
const std::array<std::vector<DimVector>, dimWorld>& axisCentroids) const;
void applyMultipliers_(Scalar& trans, void applyMultipliers_(Scalar& trans,
unsigned faceIdx, unsigned faceIdx,

View File

@ -182,25 +182,6 @@ update(bool global, const TransUpdateQuantities update_quantities,
else else
extractPermeability_(); extractPermeability_();
// calculate the axis specific centroids of all elements
std::array<std::vector<DimVector>, dimWorld> axisCentroids;
for (unsigned dimIdx = 0; dimIdx < dimWorld; ++dimIdx)
axisCentroids[dimIdx].resize(numElements);
for (const auto& elem : elements(gridView_)) {
unsigned elemIdx = elemMapper.index(elem);
// compute the axis specific "centroids" used for the transmissibilities. for
// consistency with the flow simulator, we use the element centers as
// computed by opm-parser's Opm::EclipseGrid class for all axes.
std::array<double, dimWorld> centroid = centroids_(elemIdx);
for (unsigned axisIdx = 0; axisIdx < dimWorld; ++axisIdx)
for (unsigned dimIdx = 0; dimIdx < dimWorld; ++dimIdx)
axisCentroids[axisIdx][elemIdx][dimIdx] = centroid[dimIdx];
}
// reserving some space in the hashmap upfront saves quite a bit of time because // reserving some space in the hashmap upfront saves quite a bit of time because
// resizes are costly for hashmaps and there would be quite a few of them if we // resizes are costly for hashmaps and there would be quite a few of them if we
// would not have a rough idea of how large the final map will be (the rough idea // would not have a rough idea of how large the final map will be (the rough idea
@ -272,9 +253,7 @@ update(bool global, const TransUpdateQuantities update_quantities,
faceAreaNormal, faceAreaNormal,
intersection.indexInInside(), intersection.indexInInside(),
distanceVector_(faceCenterInside, distanceVector_(faceCenterInside,
intersection.indexInInside(), centroids_(elemIdx)),
elemIdx,
axisCentroids),
permeability_[elemIdx]); permeability_[elemIdx]);
// normally there would be two half-transmissibilities that would be // normally there would be two half-transmissibilities that would be
@ -291,9 +270,7 @@ update(bool global, const TransUpdateQuantities update_quantities,
computeHalfDiffusivity_(transBoundaryEnergyIs, computeHalfDiffusivity_(transBoundaryEnergyIs,
faceAreaNormal, faceAreaNormal,
distanceVector_(faceCenterInside, distanceVector_(faceCenterInside,
intersection.indexInInside(), centroids_(elemIdx)),
elemIdx,
axisCentroids),
1.0); 1.0);
thermalHalfTransBoundary_[std::make_pair(elemIdx, boundaryIsIdx)] = thermalHalfTransBoundary_[std::make_pair(elemIdx, boundaryIsIdx)] =
transBoundaryEnergyIs; transBoundaryEnergyIs;
@ -373,17 +350,13 @@ update(bool global, const TransUpdateQuantities update_quantities,
faceAreaNormal, faceAreaNormal,
insideFaceIdx, insideFaceIdx,
distanceVector_(faceCenterInside, distanceVector_(faceCenterInside,
intersection.indexInInside(), centroids_(elemIdx)),
elemIdx,
axisCentroids),
permeability_[elemIdx]); permeability_[elemIdx]);
computeHalfTrans_(halfTrans2, computeHalfTrans_(halfTrans2,
faceAreaNormal, faceAreaNormal,
outsideFaceIdx, outsideFaceIdx,
distanceVector_(faceCenterOutside, distanceVector_(faceCenterOutside,
intersection.indexInOutside(), centroids_(outsideElemIdx)),
outsideElemIdx,
axisCentroids),
permeability_[outsideElemIdx]); permeability_[outsideElemIdx]);
applyNtg_(halfTrans1, insideFaceIdx, elemIdx, ntg); applyNtg_(halfTrans1, insideFaceIdx, elemIdx, ntg);
@ -470,16 +443,12 @@ update(bool global, const TransUpdateQuantities update_quantities,
computeHalfDiffusivity_(halfDiffusivity1, computeHalfDiffusivity_(halfDiffusivity1,
faceAreaNormal, faceAreaNormal,
distanceVector_(faceCenterInside, distanceVector_(faceCenterInside,
intersection.indexInInside(), centroids_(elemIdx)),
elemIdx,
axisCentroids),
1.0); 1.0);
computeHalfDiffusivity_(halfDiffusivity2, computeHalfDiffusivity_(halfDiffusivity2,
faceAreaNormal, faceAreaNormal,
distanceVector_(faceCenterOutside, distanceVector_(faceCenterOutside,
intersection.indexInOutside(), centroids_(outsideElemIdx)),
outsideElemIdx,
axisCentroids),
1.0); 1.0);
//TODO Add support for multipliers //TODO Add support for multipliers
thermalHalfTrans_[details::directionalIsId(elemIdx, outsideElemIdx)] = halfDiffusivity1; thermalHalfTrans_[details::directionalIsId(elemIdx, outsideElemIdx)] = halfDiffusivity1;
@ -495,16 +464,12 @@ update(bool global, const TransUpdateQuantities update_quantities,
computeHalfDiffusivity_(halfDiffusivity1, computeHalfDiffusivity_(halfDiffusivity1,
faceAreaNormal, faceAreaNormal,
distanceVector_(faceCenterInside, distanceVector_(faceCenterInside,
intersection.indexInInside(), centroids_(elemIdx)),
elemIdx,
axisCentroids),
porosity_[elemIdx]); porosity_[elemIdx]);
computeHalfDiffusivity_(halfDiffusivity2, computeHalfDiffusivity_(halfDiffusivity2,
faceAreaNormal, faceAreaNormal,
distanceVector_(faceCenterOutside, distanceVector_(faceCenterOutside,
intersection.indexInOutside(), centroids_(outsideElemIdx)),
outsideElemIdx,
axisCentroids),
porosity_[outsideElemIdx]); porosity_[outsideElemIdx]);
applyNtg_(halfDiffusivity1, insideFaceIdx, elemIdx, ntg); applyNtg_(halfDiffusivity1, insideFaceIdx, elemIdx, ntg);
@ -531,16 +496,12 @@ update(bool global, const TransUpdateQuantities update_quantities,
computeHalfDiffusivity_(halfDispersivity1, computeHalfDiffusivity_(halfDispersivity1,
faceAreaNormal, faceAreaNormal,
distanceVector_(faceCenterInside, distanceVector_(faceCenterInside,
intersection.indexInInside(), centroids_(elemIdx)),
elemIdx,
axisCentroids),
dispersion_[elemIdx]); dispersion_[elemIdx]);
computeHalfDiffusivity_(halfDispersivity2, computeHalfDiffusivity_(halfDispersivity2,
faceAreaNormal, faceAreaNormal,
distanceVector_(faceCenterOutside, distanceVector_(faceCenterOutside,
intersection.indexInOutside(), centroids_(outsideElemIdx)),
outsideElemIdx,
axisCentroids),
dispersion_[outsideElemIdx]); dispersion_[outsideElemIdx]);
applyNtg_(halfDispersivity1, insideFaceIdx, elemIdx, ntg); applyNtg_(halfDispersivity1, insideFaceIdx, elemIdx, ntg);
@ -1321,16 +1282,12 @@ computeHalfDiffusivity_(Scalar& halfDiff,
template<class Grid, class GridView, class ElementMapper, class CartesianIndexMapper, class Scalar> template<class Grid, class GridView, class ElementMapper, class CartesianIndexMapper, class Scalar>
typename Transmissibility<Grid,GridView,ElementMapper,CartesianIndexMapper,Scalar>::DimVector typename Transmissibility<Grid,GridView,ElementMapper,CartesianIndexMapper,Scalar>::DimVector
Transmissibility<Grid,GridView,ElementMapper,CartesianIndexMapper,Scalar>:: Transmissibility<Grid,GridView,ElementMapper,CartesianIndexMapper,Scalar>::
distanceVector_(const DimVector& center, distanceVector_(const DimVector& faceCenter,
int faceIdx, // in the reference element that contains the intersection const std::array<double,dimWorld>& cellCenter) const
unsigned elemIdx,
const std::array<std::vector<DimVector>, dimWorld>& axisCentroids) const
{ {
assert(faceIdx >= 0); DimVector x = faceCenter;
unsigned dimIdx = faceIdx/2; for (unsigned dimIdx = 0; dimIdx < dimWorld; ++dimIdx)
assert(dimIdx < dimWorld); x[dimIdx] -= cellCenter[dimIdx];
DimVector x = center;
x -= axisCentroids[dimIdx][elemIdx];
return x; return x;
} }