My Project
Loading...
Searching...
No Matches
TransTpfa_impl.hpp
1#include <opm/common/utility/numeric/blas_lapack.h>
3#include <opm/grid/GridHelpers.hpp>
4
5#include <cmath>
6
7namespace Dune
8{
9class CpGrid;
10}
11
12namespace Opm
13{
14namespace UgGridHelpers
15{
16int dimensions(const Dune::CpGrid&);
17
18double faceArea(const Dune::CpGrid&, int);
19}
20}
21
22namespace
23{
24inline const double* multiplyFaceNormalWithArea(const Dune::CpGrid& grid, int face_index, const double* in)
25{
26 int d=Opm::UgGridHelpers::dimensions(grid);
27 double* out=new double[d];
28 double area=Opm::UgGridHelpers::faceArea(grid, face_index);
29
30 for(int i=0;i<d;++i)
31 out[i]=in[i]*area;
32 return out;
33}
34
35inline void maybeFreeFaceNormal(const Dune::CpGrid&, const double* array)
36{
37 delete[] array;
38}
39
40inline const double* multiplyFaceNormalWithArea(const UnstructuredGrid&, int, const double* in)
41{
42 return in;
43}
44
45inline void maybeFreeFaceNormal(const UnstructuredGrid&, const double*)
46{}
47}
48
49/* ---------------------------------------------------------------------- */
50/* htrans <- sum(C(:,i) .* K(cellNo,:) .* N(:,j), 2) ./ sum(C.*C, 2) */
51/* ---------------------------------------------------------------------- */
52template<class Grid>
53void
54tpfa_htrans_compute(const Grid* G, const double *perm, double *htrans)
55/* ---------------------------------------------------------------------- */
56{
57 using namespace Opm::UgGridHelpers;
58 int d, j;
59 double s, dist, denom;
60
61 double Kn[3];
62 typename CellCentroidTraits<Grid>::IteratorType cc = beginCellCentroids(*G);
63 typename Cell2FacesTraits<Grid>::Type c2f = cell2Faces(*G);
64 typename FaceCellTraits<Grid>::Type face_cells = faceCells(*G);
65
66 const double *n;
67 const double *K;
68
69 MAT_SIZE_T nrows, ncols, ldA, incx, incy;
70 double a1, a2;
71
72 d = dimensions(*G);
73
74 nrows = ncols = ldA = d;
75 incx = incy = 1 ;
76 a1 = 1.0; a2 = 0.0 ;
77
78 for (int c =0, i = 0; c < numCells(*G); c++) {
79 K = perm + (c * d * d);
80
81 typedef typename Cell2FacesTraits<Grid>::Type::row_type FaceRow;
82 FaceRow faces = c2f[c];
83
84 for(typename FaceRow::const_iterator f=faces.begin(), end=faces.end();
85 f!=end; ++f, ++i)
86 {
87 s = 2.0*(face_cells(*f, 0) == c) - 1.0;
88 n = faceNormal(*G, *f);
89 const double* nn=multiplyFaceNormalWithArea(*G, *f, n);
90 const double* fc = &(faceCentroid(*G, *f)[0]);
91 dgemv_("No Transpose", &nrows, &ncols,
92 &a1, K, &ldA, nn, &incx, &a2, &Kn[0], &incy);
93 maybeFreeFaceNormal(*G, nn);
94
95 htrans[i] = denom = 0.0;
96 for (j = 0; j < d; j++) {
97 dist = fc[j] - getCoordinate(cc, j);
98
99 htrans[i] += s * dist * Kn[j];
100 denom += dist * dist;
101 }
102
103 assert (denom > 0);
104 htrans[i] /= denom;
105 htrans[i] = std::abs(htrans[i]);
106 }
107 // Move to next cell centroid.
108 cc = increment(cc, 1, d);
109 }
110}
111
112
113/* ---------------------------------------------------------------------- */
114template<class Grid>
115void
116tpfa_trans_compute(const Grid* G, const double *htrans, double *trans)
117/* ---------------------------------------------------------------------- */
118{
119 using namespace Opm::UgGridHelpers;
120
121 for (int f = 0; f < numFaces(*G); f++) {
122 trans[f] = 0.0;
123 }
124
125 typename Cell2FacesTraits<Grid>::Type c2f = cell2Faces(*G);
126
127 for (int c = 0, i = 0; c < numCells(*G); c++) {
128 typedef typename Cell2FacesTraits<Grid>::Type::row_type FaceRow;
129 FaceRow faces = c2f[c];
130
131 for(typename FaceRow::const_iterator f=faces.begin(), end=faces.end();
132 f!=end; ++f, ++i)
133 {
134 trans[*f] += 1.0 / htrans[i];
135 }
136 }
137
138 for (int f = 0; f < numFaces(*G); f++) {
139 trans[f] = 1.0 / trans[f];
140 }
141}
142
143
144/* ---------------------------------------------------------------------- */
145 template<class Grid>
146void
147tpfa_eff_trans_compute(const Grid* G,
148 const double *totmob,
149 const double *htrans,
150 double *trans)
151/* ---------------------------------------------------------------------- */
152{
153 using namespace Opm::UgGridHelpers;
154
155 for (int f = 0; f < numFaces(*G); f++) {
156 trans[f] = 0.0;
157 }
158
159 typename Cell2FacesTraits<Grid>::Type c2f = cell2Faces(*G);
160
161 for (int c = 0, i = 0; c < numCells(*G); c++) {
162 typedef typename Cell2FacesTraits<Grid>::Type::row_type FaceRow;
163 FaceRow faces = c2f[c];
164
165 for(typename FaceRow::const_iterator f=faces.begin(), end=faces.end();
166 f!=end; ++f, ++i)
167 {
168 trans[*f] += 1.0 / (totmob[c] * htrans[i]);
169 }
170 }
171
172 for (int f = 0; f < numFaces(*G); f++) {
173 trans[f] = 1.0 / trans[f];
174 }
175}
[ provides Dune::Grid ]
Definition CpGrid.hpp:201
Copyright 2019 Equinor AS.
Definition CartesianIndexMapper.hpp:10
T area(const Point< T, 2 > *c)
Computes the area of a 2-dimensional triangle.
Definition Volumes.hpp:118
Holds the implementation of the CpGrid as a pimple.
Definition CellQuadrature.cpp:68
Maps the grid type to the associated type of the cell to faces mapping.
Definition GridHelpers.hpp:280
Traits of the cell centroids of a grid.
Definition GridHelpers.hpp:129
Traits of the face to attached cell mappping of a grid.
Definition GridHelpers.hpp:337
Data structure for an unstructured grid, unstructured meaning that any cell may have an arbitrary num...
Definition UnstructuredGrid.h:99
Routines to assist in the calculation of two-point transmissibilities.