FacePropertiesTPSA_impl.hpp
Go to the documentation of this file.
1// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
2// vi: set et ts=4 sw=4 sts=4:
3/*
4 Copyright 2025 NORCE AS
5
6 This file is part of the Open Porous Media project (OPM).
7
8 OPM is free software: you can redistribute it and/or modify
9 it under the terms of the GNU General Public License as published by
10 the Free Software Foundation, either version 2 of the License, or
11 (at your option) any later version.
12
13 OPM is distributed in the hope that it will be useful,
14 but WITHOUT ANY WARRANTY; without even the implied warranty of
15 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
16 GNU General Public License for more details.
17
18 You should have received a copy of the GNU General Public License
19 along with OPM. If not, see <http://www.gnu.org/licenses/>.
20
21 Consult the COPYING file in the top-level source directory of this
22 module for the precise wording of the license and the list of
23 copyright holders.
24*/
25#ifndef TPSA_FACE_PROPERTIES_IMPL_HPP
26#define TPSA_FACE_PROPERTIES_IMPL_HPP
27
28#ifndef TPSA_FACE_PROPERTIES_HPP
29#include <config.h>
31#endif
32
33#include <dune/grid/common/mcmgmapper.hh>
34
35#include <opm/common/utility/ThreadSafeMapBuilder.hpp>
36
37#include <opm/grid/utility/ElementChunks.hpp>
38
39#include <opm/input/eclipse/EclipseState/EclipseState.hpp>
40
42
43#include <algorithm>
44#include <stdexcept>
45#include <utility>
46
47
48namespace Opm {
49
50// Copied from Opm::Transmissibility class
51namespace details {
52 constexpr unsigned elemIdxShift = 32; // bits
53
54 std::uint64_t isIdTPSA(std::uint32_t elemIdx1, std::uint32_t elemIdx2)
55 {
56 const std::uint32_t elemAIdx = std::min(elemIdx1, elemIdx2);
57 const std::uint64_t elemBIdx = std::max(elemIdx1, elemIdx2);
58
59 return (elemBIdx << elemIdxShift) + elemAIdx;
60 }
61} // namespace Opm::details
62
63// /////
64// Public functions
65// ////
75template<class Grid, class GridView, class ElementMapper, class CartesianIndexMapper, class Scalar>
77FacePropertiesTPSA(const EclipseState& eclState,
78 const GridView& gridView,
79 const CartesianIndexMapper& cartMapper,
80 const Grid& grid,
81 std::function<std::array<double,dimWorld>(int)> centroids)
82 : eclState_(eclState)
83 , gridView_(gridView)
84 , cartMapper_(cartMapper)
85 , grid_(grid)
86 , centroids_(centroids)
87 , lookUpData_(gridView)
88 , lookUpCartesianData_(gridView, cartMapper)
89{ }
90
94template<class Grid, class GridView, class ElementMapper, class CartesianIndexMapper, class Scalar>
97{
98 update();
99}
100
104template<class Grid, class GridView, class ElementMapper, class CartesianIndexMapper, class Scalar>
106update()
107{
108 // Number of elements
109 ElementMapper elemMapper(gridView_, Dune::mcmgElementLayout());
110 unsigned numElements = elemMapper.size();
111
112 // Extract shear modulus (Lame's sec. params)
113 extractSModulus_();
114
115 // Init. containers
116 // Note (from Transmissibility::update): Reserving some space in the hashmap upfront saves quite a bit of time
117 // because resizes are costly for hashmaps and there would be quite a few of them if we would not have a rough idea
118 // of how large the final map will be (the rough idea is a conforming Cartesian grid).
119 const int num_threads = ThreadManager::maxThreads();
120 weightsAvg_.clear();
121 weightsProd_.clear();
122 distance_.clear();
123 faceNormal_.clear();
124 if (num_threads == 1) {
125 weightsAvg_.reserve(numElements * 3 * 1.05);
126 weightsProd_.reserve(numElements * 3 * 1.05);
127 distance_.reserve(numElements * 3 * 1.05);
128 faceNormal_.reserve(numElements * 3 * 1.05);
129 }
130 weightsAvgBoundary_.clear();
131 weightsProdBoundary_.clear();
132 distanceBoundary_.clear();
133 faceNormalBoundary_.clear();
134
135 // Fill the centroids cache to avoid repeated calculations in loops below
136 centroids_cache_.resize(gridView_.size(0));
137 for (const auto& elem : elements(gridView_)) {
138 const unsigned elemIdx = elemMapper.index(elem);
139 centroids_cache_[elemIdx] = centroids_(elemIdx);
140 }
141
142 // Initialize thread safe insert_or_assign for face properties in the grid and separate for boundaries
143 ThreadSafeMapBuilder weightsAvgMap(weightsAvg_, num_threads, MapBuilderInsertionMode::Insert_Or_Assign);
144 ThreadSafeMapBuilder weightsProdMap(weightsProd_, num_threads, MapBuilderInsertionMode::Insert_Or_Assign);
145 ThreadSafeMapBuilder distanceMap(distance_, num_threads, MapBuilderInsertionMode::Insert_Or_Assign);
146 ThreadSafeMapBuilder faceNormalMap(faceNormal_, num_threads, MapBuilderInsertionMode::Insert_Or_Assign);
147
148 ThreadSafeMapBuilder weightsAvgBoundaryMap(weightsAvgBoundary_, num_threads, MapBuilderInsertionMode::Insert_Or_Assign);
149 ThreadSafeMapBuilder weightsProdBoundaryMap(weightsProdBoundary_, num_threads, MapBuilderInsertionMode::Insert_Or_Assign);
150 ThreadSafeMapBuilder distanceBoundaryMap(distanceBoundary_, num_threads, MapBuilderInsertionMode::Insert_Or_Assign);
151 ThreadSafeMapBuilder faceNormalBoundaryMap(faceNormalBoundary_, num_threads, MapBuilderInsertionMode::Insert_Or_Assign);
152
153#ifdef _OPENMP
154#pragma omp parallel for
155#endif
156 // Loop over grid element an compute face properties
157 for (const auto& chunk : ElementChunks(gridView_, Dune::Partitions::all, num_threads)) {
158 for (const auto& elem : chunk) {
159 // Init. face info for inside/outside cells
160 FaceInfo inside;
161 FaceInfo outside;
162 DimVector faceNormal;
163
164 // Set inside info
165 inside.elemIdx = elemMapper.index(elem);
166 inside.cartElemIdx = lookUpCartesianData_.
167 template getFieldPropCartesianIdx<Grid>(inside.elemIdx);
168
169 // Loop over intersection to neighboring cells
170 unsigned boundaryIsIdx = 0;
171 for (const auto& intersection : intersections(gridView_, elem)) {
172 // Handle grid boundary
173 if (intersection.boundary()) {
174 // One-sided cell properties
175 const auto& geometry = intersection.geometry();
176 inside.faceCenter = geometry.center();
177 faceNormal = intersection.centerUnitOuterNormal();
178
179 // Face properties on boundary
180 const auto index_pair = std::make_pair(inside.elemIdx, boundaryIsIdx);
181 Scalar distBound = computeDistance_(distanceVector_(inside.faceCenter, inside.elemIdx), faceNormal);
182 distanceBoundaryMap.insert_or_assign(index_pair, distBound);
183
184 Scalar weightsAvgBound = 1.0; // w_j = 0 -> w_avg = wi / (wi + 0)
185 Scalar weightsProdBound = 0.0; // w_j = 0 -> w_prod = wi * 0
186 weightsAvgBoundaryMap.insert_or_assign(index_pair, weightsAvgBound);
187 weightsProdBoundaryMap.insert_or_assign(index_pair, weightsProdBound);
188
189 faceNormalBoundaryMap.insert_or_assign(index_pair, faceNormal);
190
191 ++boundaryIsIdx;
192 continue;
193 }
194
195 // Handle intersection on process boundary (i.e., neighbor on different rank)
196 if (!intersection.neighbor()) {
197 ++boundaryIsIdx;
198 continue;
199 }
200
201 // Set outside info
202 const auto& outsideElem = intersection.outside();
203 outside.elemIdx = elemMapper.index(outsideElem);
204 outside.cartElemIdx = lookUpCartesianData_.
205 template getFieldPropCartesianIdx<Grid>(outside.elemIdx);
206
207 // Skip intersections that have already been processed
208 if (std::tie(inside.cartElemIdx, inside.elemIdx) > std::tie(outside.cartElemIdx, outside.elemIdx)) {
209 continue;
210 }
211
212 // Face indices for this intersection
213 inside.faceIdx = intersection.indexInInside();
214 outside.faceIdx = intersection.indexInOutside();
215
216 // Set NNC face properties to zero
217 if (inside.faceIdx == -1) {
218 const auto id = details::isIdTPSA(inside.elemIdx, outside.elemIdx);
219 weightsAvgMap.insert_or_assign(id, 0.0);
220 distanceMap.insert_or_assign(id, 0.0);
221 faceNormalMap.insert_or_assign(id, DimVector{0.0, 0.0, 0.0});
222
223 continue;
224 }
225
226 // Compute cell properties from grid
227 typename std::is_same<Grid, Dune::CpGrid>::type isCpGrid;
228 computeCellProperties(intersection,
229 inside,
230 outside,
231 faceNormal,
232 isCpGrid);
233
234 // Compute face properties
235 const auto id = details::isIdTPSA(inside.elemIdx, outside.elemIdx);
236 Scalar dist_in = computeDistance_(distanceVector_(inside.faceCenter, inside.elemIdx), faceNormal);
237 Scalar dist_out = computeDistance_(distanceVector_(outside.faceCenter, outside.elemIdx), faceNormal);
238 distanceMap.insert_or_assign(id, dist_in + dist_out);
239
240 Scalar weight_in = computeWeight_(dist_in, sModulus_[inside.elemIdx]);
241 Scalar weight_out = computeWeight_(dist_out, sModulus_[outside.elemIdx]);
242 Scalar weightsAvg = weight_in / (weight_in + weight_out);
243 Scalar weightsProd = weight_in * weight_out;
244 weightsAvgMap.insert_or_assign(id, weightsAvg);
245 weightsProdMap.insert_or_assign(id, weightsProd);
246
247 faceNormalMap.insert_or_assign(id, faceNormal);
248 }
249 }
250 }
251
252 // Clear centroid cache
253 centroids_cache_.clear();
254
255#ifdef _OPENMP
256#pragma omp parallel sections
257#endif
258 {
259#ifdef _OPENMP
260#pragma omp section
261#endif
262 weightsAvgMap.finalize();
263#ifdef _OPENMP
264#pragma omp section
265#endif
266 weightsProdMap.finalize();
267#ifdef _OPENMP
268#pragma omp section
269#endif
270 distanceMap.finalize();
271#ifdef _OPENMP
272#pragma omp section
273#endif
274 faceNormalMap.finalize();
275#ifdef _OPENMP
276#pragma omp section
277#endif
278 weightsAvgBoundaryMap.finalize();
279#ifdef _OPENMP
280#pragma omp section
281#endif
282 weightsProdBoundaryMap.finalize();
283#ifdef _OPENMP
284#pragma omp section
285#endif
286 distanceBoundaryMap.finalize();
287#ifdef _OPENMP
288#pragma omp section
289#endif
290 faceNormalBoundaryMap.finalize();
291 }
292}
293
301template<class Grid, class GridView, class ElementMapper, class CartesianIndexMapper, class Scalar>
303weightAverage(unsigned elemIdx1, unsigned elemIdx2) const
304{
305 auto tmp_whgt = weightsAvg_.at(details::isIdTPSA(elemIdx1, elemIdx2));
306
307 auto cartIdx1 = lookUpCartesianData_.
308 template getFieldPropCartesianIdx<Grid>(elemIdx1);
309 auto cartIdx2 = lookUpCartesianData_.
310 template getFieldPropCartesianIdx<Grid>(elemIdx2);
311 if (cartIdx1 < cartIdx2) {
312 return tmp_whgt;
313 }
314 else {
315 return 1.0 - tmp_whgt;
316 }
317}
318
326template<class Grid, class GridView, class ElementMapper, class CartesianIndexMapper, class Scalar>
328weightAverageBoundary(unsigned elemIdx, unsigned boundaryFaceIdx) const
329{
330 return weightsAvgBoundary_.at(std::make_pair(elemIdx, boundaryFaceIdx));
331}
332
340template<class Grid, class GridView, class ElementMapper, class CartesianIndexMapper, class Scalar>
342weightProduct(unsigned elemIdx1, unsigned elemIdx2) const
343{
344 return weightsProd_.at(details::isIdTPSA(elemIdx1, elemIdx2));
345}
346
354template<class Grid, class GridView, class ElementMapper, class CartesianIndexMapper, class Scalar>
356weightProductBoundary(unsigned elemIdx, unsigned boundaryFaceIdx) const
357{
358 return weightsProdBoundary_.at(std::make_pair(elemIdx, boundaryFaceIdx));
359}
360
368template<class Grid, class GridView, class ElementMapper, class CartesianIndexMapper, class Scalar>
370normalDistance(unsigned elemIdx1, unsigned elemIdx2) const
371{
372 return distance_.at(details::isIdTPSA(elemIdx1, elemIdx2));
373}
374
382template<class Grid, class GridView, class ElementMapper, class CartesianIndexMapper, class Scalar>
384normalDistanceBoundary(unsigned elemIdx, unsigned boundaryFaceIdx) const
385{
386 return distanceBoundary_.at(std::make_pair(elemIdx, boundaryFaceIdx));
387}
388
398template<class Grid, class GridView, class ElementMapper, class CartesianIndexMapper, class Scalar>
401cellFaceNormal(unsigned elemIdx1, unsigned elemIdx2)
402{
403 auto cartIdx1 = lookUpCartesianData_.
404 template getFieldPropCartesianIdx<Grid>(elemIdx1);
405 auto cartIdx2 = lookUpCartesianData_.
406 template getFieldPropCartesianIdx<Grid>(elemIdx2);
407 int sign = (cartIdx1 < cartIdx2) ? 1 : -1;
408 return sign * faceNormal_.at(details::isIdTPSA(elemIdx1, elemIdx2));
409}
410
420template<class Grid, class GridView, class ElementMapper, class CartesianIndexMapper, class Scalar>
423cellFaceNormalBoundary(unsigned elemIdx, unsigned boundaryFaceIdx) const
424{
425 return faceNormalBoundary_.at(std::make_pair(elemIdx, boundaryFaceIdx));
426}
427
428// /////
429// Protected functions
430// ////
438template<class Grid, class GridView, class ElementMapper, class CartesianIndexMapper, class Scalar>
440computeDistance_(const DimVector& distVec, const DimVector& faceNormal)
441{
442 return std::abs(Dune::dot(faceNormal, distVec));
443}
444
455template<class Grid, class GridView, class ElementMapper, class CartesianIndexMapper, class Scalar>
456template<class Intersection>
458computeCellProperties([[maybe_unused]] const Intersection& intersection,
459 [[maybe_unused]] FaceInfo& inside,
460 [[maybe_unused]] FaceInfo& outside,
461 [[maybe_unused]] DimVector& faceNormal,
462 /*isCpGrid=*/std::false_type) const
463{
464 throw std::runtime_error("TPSA not implemented for DUNE grid types other than CpGrid!");
465}
466
477template<class Grid, class GridView, class ElementMapper, class CartesianIndexMapper, class Scalar>
478template<class Intersection>
480computeCellProperties(const Intersection& intersection,
481 FaceInfo& inside,
482 FaceInfo& outside,
483 DimVector& faceNormal,
484 /*isCpGrid=*/std::true_type) const
485{
486 int faceIdx = intersection.id();
487 if (grid_.maxLevel() == 0) {
488 // Face center coordinates
489 inside.faceCenter = grid_.faceCenterEcl(inside.elemIdx, inside.faceIdx, intersection);
490 outside.faceCenter = grid_.faceCenterEcl(outside.elemIdx, outside.faceIdx, intersection);
491
492 // Face normal, ensuring it points from cell with lower to higher (global grid) index
493 faceNormal = grid_.faceNormal(faceIdx);
494 auto cartFaceCell0 = lookUpCartesianData_.
495 template getFieldPropCartesianIdx<Grid>(grid_.faceCell(faceIdx, 0));
496 auto cartFaceCell1 = lookUpCartesianData_.
497 template getFieldPropCartesianIdx<Grid>(grid_.faceCell(faceIdx, 1));
498 if (cartFaceCell0 > cartFaceCell1) {
499 faceNormal *= -1;
500 }
501 }
502 else {
503 throw std::runtime_error("TPSA not implemented with LGR");
504 }
505
506}
507
515template<class Grid, class GridView, class ElementMapper, class CartesianIndexMapper, class Scalar>
517computeWeight_(const Scalar distance, const Scalar smod)
518{
519 return distance / smod;
520}
521
529template<class Grid, class GridView, class ElementMapper, class CartesianIndexMapper, class Scalar>
532distanceVector_(const DimVector& faceCenter, const unsigned& cellIdx) const
533{
534 const auto& cellCenter = centroids_cache_.empty() ? centroids_(cellIdx)
535 : centroids_cache_[cellIdx];
536 DimVector x = faceCenter;
537 for (unsigned dimIdx = 0; dimIdx < dimWorld; ++dimIdx) {
538 x[dimIdx] -= cellCenter[dimIdx];
539 }
540
541 return x;
542}
543
551template<class Grid, class GridView, class ElementMapper, class CartesianIndexMapper, class Scalar>
554{
555 unsigned numElem = gridView_.size(/*codim=*/0);
556 sModulus_.resize(numElem);
557
558 const auto& fp = eclState_.fieldProps();
559 std::vector<double> sModulusData;
560 if (fp.has_double("SMODULUS")) {
561 sModulusData = this->lookUpData_.assignFieldPropsDoubleOnLeaf(fp, "SMODULUS");
562 }
563 else if (fp.has_double("YMODULE") && fp.has_double("LAME")) {
564 // Convert from Young's modulus and Lame's first parameter
565 const std::vector<double>& ymodulus = this->lookUpData_.assignFieldPropsDoubleOnLeaf(fp, "YMODULE");
566 const std::vector<double>& lameParam = this->lookUpData_.assignFieldPropsDoubleOnLeaf(fp, "LAME");
567 sModulusData.resize(numElem);
568 for (std::size_t i = 0; i < sModulusData.size(); ++i) {
569 const double r = std::sqrt(ymodulus[i] * ymodulus[i] + 9 * lameParam[i] * lameParam[i]
570 + 2 * ymodulus[i] * lameParam[i]);
571 sModulusData[i] = (ymodulus[i] - 3 * lameParam[i] + r) / 4.0;
572 }
573 }
574 else if (fp.has_double("YMODULE") && fp.has_double("PRATIO")) {
575 const std::vector<double>& ymodulus = this->lookUpData_.assignFieldPropsDoubleOnLeaf(fp, "YMODULE");
576 const std::vector<double>& pratio = this->lookUpData_.assignFieldPropsDoubleOnLeaf(fp, "PRATIO");
577 sModulusData.resize(numElem);
578 for (std::size_t i = 0; i < sModulusData.size(); ++i) {
579 sModulusData[i] = ymodulus[i] / (2 * (1 + pratio[i]));
580 }
581 }
582 else if (fp.has_double("LAME") && fp.has_double("PRATIO")) {
583 const std::vector<double>& lameParam = this->lookUpData_.assignFieldPropsDoubleOnLeaf(fp, "LAME");
584 const std::vector<double>& pratio = this->lookUpData_.assignFieldPropsDoubleOnLeaf(fp, "PRATIO");
585 sModulusData.resize(numElem);
586 for (std::size_t i = 0; i < sModulusData.size(); ++i) {
587 sModulusData[i] = lameParam[i] * (1 - 2 * pratio[i]) / (2 * pratio[i]);
588 }
589 }
590 else {
591 throw std::logic_error("Cannot read shear modulus data from ecl state, SMODULUS keyword missing, "
592 "and one of the following keyword pairs are missing for conversion: "
593 "(YMODULE, LAME), (YMODULE, PRATIO) and (LAME, PRATIO)!");
594 }
595
596 // Assign shear modulus
597 for (std::size_t dofIdx = 0; dofIdx < numElem; ++ dofIdx) {
598 sModulus_[dofIdx] = sModulusData[dofIdx];
599 }
600}
601
602} // namespace Opm
603
604#endif
void extractSModulus_()
Extract shear modulus from eclState.
Definition: FacePropertiesTPSA_impl.hpp:553
void finishInit()
Compute TPSA face properties.
Definition: FacePropertiesTPSA_impl.hpp:96
void computeCellProperties(const Intersection &intersection, FaceInfo &inside, FaceInfo &outside, DimVector &faceNormal, std::false_type) const
Compute face properties from general DUNE grid.
Definition: FacePropertiesTPSA_impl.hpp:458
Scalar weightProduct(unsigned elemIdx1, unsigned elemIdx2) const
Product of weights at interface between two elements.
Definition: FacePropertiesTPSA_impl.hpp:342
DimVector cellFaceNormal(unsigned elemIdx1, unsigned elemIdx2)
Cell face normal at interface between two elements.
Definition: FacePropertiesTPSA_impl.hpp:401
Dune::FieldVector< Scalar, dimWorld > DimVector
Definition: FacePropertiesTPSA.hpp:53
Scalar weightAverage(unsigned elemIdx1, unsigned elemIdx2) const
Average (half-)weight at interface between two elements.
Definition: FacePropertiesTPSA_impl.hpp:303
Scalar normalDistanceBoundary(unsigned elemIdx1, unsigned boundaryFaceIdx) const
Distance to boundary interface.
Definition: FacePropertiesTPSA_impl.hpp:384
Scalar weightAverageBoundary(unsigned elemIdx1, unsigned boundaryFaceIdx) const
Average (half-)weight at boundary interface.
Definition: FacePropertiesTPSA_impl.hpp:328
const DimVector & cellFaceNormalBoundary(unsigned elemIdx1, unsigned boundaryFaceIdx) const
Cell face normal of boundary interface.
Definition: FacePropertiesTPSA_impl.hpp:423
Scalar computeWeight_(const Scalar distance, const Scalar smod)
Compute weight ratio between distance and shear modulus.
Definition: FacePropertiesTPSA_impl.hpp:517
void update()
Compute TPSA face properties.
Definition: FacePropertiesTPSA_impl.hpp:106
Scalar computeDistance_(const DimVector &distVec, const DimVector &faceNormal)
Compute normal distance from cell center to face center.
Definition: FacePropertiesTPSA_impl.hpp:440
DimVector distanceVector_(const DimVector &faceCenter, const unsigned &cellIdx) const
Distance vector from cell center to face center.
Definition: FacePropertiesTPSA_impl.hpp:532
Scalar normalDistance(unsigned elemIdx1, unsigned elemIdx2) const
Distance between two elements.
Definition: FacePropertiesTPSA_impl.hpp:370
FacePropertiesTPSA(const EclipseState &eclState, const GridView &gridView, const CartesianIndexMapper &cartMapper, const Grid &grid, std::function< std::array< double, dimWorld >(int)> centroids)
Constructor.
Definition: FacePropertiesTPSA_impl.hpp:77
Scalar weightProductBoundary(unsigned elemIdx1, unsigned boundaryFaceIdx) const
Product of weights at boundary interface.
Definition: FacePropertiesTPSA_impl.hpp:356
static unsigned maxThreads()
Return the maximum number of threads of the current process.
Definition: threadmanager.hpp:66
constexpr unsigned elemIdxShift
Definition: FacePropertiesTPSA_impl.hpp:52
std::uint64_t isIdTPSA(std::uint32_t elemIdx1, std::uint32_t elemIdx2)
Definition: FacePropertiesTPSA_impl.hpp:54
Definition: blackoilbioeffectsmodules.hh:45
Definition: FacePropertiesTPSA.hpp:83
int faceIdx
Definition: FacePropertiesTPSA.hpp:85
unsigned cartElemIdx
Definition: FacePropertiesTPSA.hpp:87
unsigned elemIdx
Definition: FacePropertiesTPSA.hpp:86
DimVector faceCenter
Definition: FacePropertiesTPSA.hpp:84