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 weightsProdMap.insert_or_assign(id, 0.0);
221 distanceMap.insert_or_assign(id, 1.0);
222 faceNormalMap.insert_or_assign(id, DimVector{0.0, 0.0, 0.0});
223
224 continue;
225 }
226
227 // Compute cell properties from grid
228 typename std::is_same<Grid, Dune::CpGrid>::type isCpGrid;
229 computeCellProperties(intersection,
230 inside,
231 outside,
232 faceNormal,
233 isCpGrid);
234
235 // Compute face properties
236 const auto id = details::isIdTPSA(inside.elemIdx, outside.elemIdx);
237 Scalar dist_in = computeDistance_(distanceVector_(inside.faceCenter, inside.elemIdx), faceNormal);
238 Scalar dist_out = computeDistance_(distanceVector_(outside.faceCenter, outside.elemIdx), faceNormal);
239 distanceMap.insert_or_assign(id, dist_in + dist_out);
240
241 Scalar weight_in = computeWeight_(dist_in, sModulus_[inside.elemIdx]);
242 Scalar weight_out = computeWeight_(dist_out, sModulus_[outside.elemIdx]);
243 Scalar weightsAvg = weight_in / (weight_in + weight_out);
244 Scalar weightsProd = weight_in * weight_out;
245 weightsAvgMap.insert_or_assign(id, weightsAvg);
246 weightsProdMap.insert_or_assign(id, weightsProd);
247
248 faceNormalMap.insert_or_assign(id, faceNormal);
249 }
250 }
251 }
252
253 // Clear centroid cache
254 centroids_cache_.clear();
255
256#ifdef _OPENMP
257#pragma omp parallel sections
258#endif
259 {
260#ifdef _OPENMP
261#pragma omp section
262#endif
263 weightsAvgMap.finalize();
264#ifdef _OPENMP
265#pragma omp section
266#endif
267 weightsProdMap.finalize();
268#ifdef _OPENMP
269#pragma omp section
270#endif
271 distanceMap.finalize();
272#ifdef _OPENMP
273#pragma omp section
274#endif
275 faceNormalMap.finalize();
276#ifdef _OPENMP
277#pragma omp section
278#endif
279 weightsAvgBoundaryMap.finalize();
280#ifdef _OPENMP
281#pragma omp section
282#endif
283 weightsProdBoundaryMap.finalize();
284#ifdef _OPENMP
285#pragma omp section
286#endif
287 distanceBoundaryMap.finalize();
288#ifdef _OPENMP
289#pragma omp section
290#endif
291 faceNormalBoundaryMap.finalize();
292 }
293}
294
302template<class Grid, class GridView, class ElementMapper, class CartesianIndexMapper, class Scalar>
304weightAverage(unsigned elemIdx1, unsigned elemIdx2) const
305{
306 auto tmp_whgt = weightsAvg_.at(details::isIdTPSA(elemIdx1, elemIdx2));
307
308 auto cartIdx1 = lookUpCartesianData_.
309 template getFieldPropCartesianIdx<Grid>(elemIdx1);
310 auto cartIdx2 = lookUpCartesianData_.
311 template getFieldPropCartesianIdx<Grid>(elemIdx2);
312 if (cartIdx1 < cartIdx2) {
313 return tmp_whgt;
314 }
315 else {
316 return 1.0 - tmp_whgt;
317 }
318}
319
327template<class Grid, class GridView, class ElementMapper, class CartesianIndexMapper, class Scalar>
329weightAverageBoundary(unsigned elemIdx, unsigned boundaryFaceIdx) const
330{
331 return weightsAvgBoundary_.at(std::make_pair(elemIdx, boundaryFaceIdx));
332}
333
341template<class Grid, class GridView, class ElementMapper, class CartesianIndexMapper, class Scalar>
343weightProduct(unsigned elemIdx1, unsigned elemIdx2) const
344{
345 return weightsProd_.at(details::isIdTPSA(elemIdx1, elemIdx2));
346}
347
355template<class Grid, class GridView, class ElementMapper, class CartesianIndexMapper, class Scalar>
357weightProductBoundary(unsigned elemIdx, unsigned boundaryFaceIdx) const
358{
359 return weightsProdBoundary_.at(std::make_pair(elemIdx, boundaryFaceIdx));
360}
361
369template<class Grid, class GridView, class ElementMapper, class CartesianIndexMapper, class Scalar>
371normalDistance(unsigned elemIdx1, unsigned elemIdx2) const
372{
373 return distance_.at(details::isIdTPSA(elemIdx1, elemIdx2));
374}
375
383template<class Grid, class GridView, class ElementMapper, class CartesianIndexMapper, class Scalar>
385normalDistanceBoundary(unsigned elemIdx, unsigned boundaryFaceIdx) const
386{
387 return distanceBoundary_.at(std::make_pair(elemIdx, boundaryFaceIdx));
388}
389
399template<class Grid, class GridView, class ElementMapper, class CartesianIndexMapper, class Scalar>
402cellFaceNormal(unsigned elemIdx1, unsigned elemIdx2)
403{
404 auto cartIdx1 = lookUpCartesianData_.
405 template getFieldPropCartesianIdx<Grid>(elemIdx1);
406 auto cartIdx2 = lookUpCartesianData_.
407 template getFieldPropCartesianIdx<Grid>(elemIdx2);
408 int sign = (cartIdx1 < cartIdx2) ? 1 : -1;
409 return sign * faceNormal_.at(details::isIdTPSA(elemIdx1, elemIdx2));
410}
411
421template<class Grid, class GridView, class ElementMapper, class CartesianIndexMapper, class Scalar>
424cellFaceNormalBoundary(unsigned elemIdx, unsigned boundaryFaceIdx) const
425{
426 return faceNormalBoundary_.at(std::make_pair(elemIdx, boundaryFaceIdx));
427}
428
429// /////
430// Protected functions
431// ////
439template<class Grid, class GridView, class ElementMapper, class CartesianIndexMapper, class Scalar>
441computeDistance_(const DimVector& distVec, const DimVector& faceNormal)
442{
443 return std::abs(Dune::dot(faceNormal, distVec));
444}
445
456template<class Grid, class GridView, class ElementMapper, class CartesianIndexMapper, class Scalar>
457template<class Intersection>
459computeCellProperties([[maybe_unused]] const Intersection& intersection,
460 [[maybe_unused]] FaceInfo& inside,
461 [[maybe_unused]] FaceInfo& outside,
462 [[maybe_unused]] DimVector& faceNormal,
463 /*isCpGrid=*/std::false_type) const
464{
465 throw std::runtime_error("TPSA not implemented for DUNE grid types other than CpGrid!");
466}
467
478template<class Grid, class GridView, class ElementMapper, class CartesianIndexMapper, class Scalar>
479template<class Intersection>
481computeCellProperties(const Intersection& intersection,
482 FaceInfo& inside,
483 FaceInfo& outside,
484 DimVector& faceNormal,
485 /*isCpGrid=*/std::true_type) const
486{
487 int faceIdx = intersection.id();
488 if (grid_.maxLevel() == 0) {
489 // Face center coordinates
490 inside.faceCenter = grid_.faceCenterEcl(inside.elemIdx, inside.faceIdx, intersection);
491 outside.faceCenter = grid_.faceCenterEcl(outside.elemIdx, outside.faceIdx, intersection);
492
493 // Face normal, ensuring it points from cell with lower to higher (global grid) index
494 faceNormal = grid_.faceNormal(faceIdx);
495 auto cartFaceCell0 = lookUpCartesianData_.
496 template getFieldPropCartesianIdx<Grid>(grid_.faceCell(faceIdx, 0));
497 auto cartFaceCell1 = lookUpCartesianData_.
498 template getFieldPropCartesianIdx<Grid>(grid_.faceCell(faceIdx, 1));
499 if (cartFaceCell0 > cartFaceCell1) {
500 faceNormal *= -1;
501 }
502 }
503 else {
504 throw std::runtime_error("TPSA not implemented with LGR");
505 }
506
507}
508
516template<class Grid, class GridView, class ElementMapper, class CartesianIndexMapper, class Scalar>
518computeWeight_(const Scalar distance, const Scalar smod)
519{
520 return distance / smod;
521}
522
530template<class Grid, class GridView, class ElementMapper, class CartesianIndexMapper, class Scalar>
533distanceVector_(const DimVector& faceCenter, const unsigned& cellIdx) const
534{
535 const auto& cellCenter = centroids_cache_.empty() ? centroids_(cellIdx)
536 : centroids_cache_[cellIdx];
537 DimVector x = faceCenter;
538 for (unsigned dimIdx = 0; dimIdx < dimWorld; ++dimIdx) {
539 x[dimIdx] -= cellCenter[dimIdx];
540 }
541
542 return x;
543}
544
552template<class Grid, class GridView, class ElementMapper, class CartesianIndexMapper, class Scalar>
555{
556 unsigned numElem = gridView_.size(/*codim=*/0);
557 sModulus_.resize(numElem);
558
559 const auto& fp = eclState_.fieldProps();
560 std::vector<double> sModulusData;
561 if (fp.has_double("SMODULUS")) {
562 sModulusData = this->lookUpData_.assignFieldPropsDoubleOnLeaf(fp, "SMODULUS");
563 }
564 else if (fp.has_double("YMODULE") && fp.has_double("LAME")) {
565 // Convert from Young's modulus and Lame's first parameter
566 const std::vector<double>& ymodulus = this->lookUpData_.assignFieldPropsDoubleOnLeaf(fp, "YMODULE");
567 const std::vector<double>& lameParam = this->lookUpData_.assignFieldPropsDoubleOnLeaf(fp, "LAME");
568 sModulusData.resize(numElem);
569 for (std::size_t i = 0; i < sModulusData.size(); ++i) {
570 const double r = std::sqrt(ymodulus[i] * ymodulus[i] + 9 * lameParam[i] * lameParam[i]
571 + 2 * ymodulus[i] * lameParam[i]);
572 sModulusData[i] = (ymodulus[i] - 3 * lameParam[i] + r) / 4.0;
573 }
574 }
575 else if (fp.has_double("YMODULE") && fp.has_double("PRATIO")) {
576 const std::vector<double>& ymodulus = this->lookUpData_.assignFieldPropsDoubleOnLeaf(fp, "YMODULE");
577 const std::vector<double>& pratio = this->lookUpData_.assignFieldPropsDoubleOnLeaf(fp, "PRATIO");
578 sModulusData.resize(numElem);
579 for (std::size_t i = 0; i < sModulusData.size(); ++i) {
580 sModulusData[i] = ymodulus[i] / (2 * (1 + pratio[i]));
581 }
582 }
583 else if (fp.has_double("LAME") && fp.has_double("PRATIO")) {
584 const std::vector<double>& lameParam = this->lookUpData_.assignFieldPropsDoubleOnLeaf(fp, "LAME");
585 const std::vector<double>& pratio = this->lookUpData_.assignFieldPropsDoubleOnLeaf(fp, "PRATIO");
586 sModulusData.resize(numElem);
587 for (std::size_t i = 0; i < sModulusData.size(); ++i) {
588 sModulusData[i] = lameParam[i] * (1 - 2 * pratio[i]) / (2 * pratio[i]);
589 }
590 }
591 else {
592 throw std::logic_error("Cannot read shear modulus data from ecl state, SMODULUS keyword missing, "
593 "and one of the following keyword pairs are missing for conversion: "
594 "(YMODULE, LAME), (YMODULE, PRATIO) and (LAME, PRATIO)!");
595 }
596
597 // Assign shear modulus
598 for (std::size_t dofIdx = 0; dofIdx < numElem; ++ dofIdx) {
599 sModulus_[dofIdx] = sModulusData[dofIdx];
600 }
601}
602
603} // namespace Opm
604
605#endif
void extractSModulus_()
Extract shear modulus from eclState.
Definition: FacePropertiesTPSA_impl.hpp:554
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:459
Scalar weightProduct(unsigned elemIdx1, unsigned elemIdx2) const
Product of weights at interface between two elements.
Definition: FacePropertiesTPSA_impl.hpp:343
DimVector cellFaceNormal(unsigned elemIdx1, unsigned elemIdx2)
Cell face normal at interface between two elements.
Definition: FacePropertiesTPSA_impl.hpp:402
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:304
Scalar normalDistanceBoundary(unsigned elemIdx1, unsigned boundaryFaceIdx) const
Distance to boundary interface.
Definition: FacePropertiesTPSA_impl.hpp:385
Scalar weightAverageBoundary(unsigned elemIdx1, unsigned boundaryFaceIdx) const
Average (half-)weight at boundary interface.
Definition: FacePropertiesTPSA_impl.hpp:329
const DimVector & cellFaceNormalBoundary(unsigned elemIdx1, unsigned boundaryFaceIdx) const
Cell face normal of boundary interface.
Definition: FacePropertiesTPSA_impl.hpp:424
Scalar computeWeight_(const Scalar distance, const Scalar smod)
Compute weight ratio between distance and shear modulus.
Definition: FacePropertiesTPSA_impl.hpp:518
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:441
DimVector distanceVector_(const DimVector &faceCenter, const unsigned &cellIdx) const
Distance vector from cell center to face center.
Definition: FacePropertiesTPSA_impl.hpp:533
Scalar normalDistance(unsigned elemIdx1, unsigned elemIdx2) const
Distance between two elements.
Definition: FacePropertiesTPSA_impl.hpp:371
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:357
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