fvbasegradientcalculator.hh
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  This file is part of the Open Porous Media project (OPM).
5 
6  OPM is free software: you can redistribute it and/or modify
7  it under the terms of the GNU General Public License as published by
8  the Free Software Foundation, either version 2 of the License, or
9  (at your option) any later version.
10 
11  OPM is distributed in the hope that it will be useful,
12  but WITHOUT ANY WARRANTY; without even the implied warranty of
13  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
14  GNU General Public License for more details.
15 
16  You should have received a copy of the GNU General Public License
17  along with OPM. If not, see <http://www.gnu.org/licenses/>.
18 
19  Consult the COPYING file in the top-level source directory of this
20  module for the precise wording of the license and the list of
21  copyright holders.
22 */
28 #ifndef EWOMS_FV_BASE_GRADIENT_CALCULATOR_HH
29 #define EWOMS_FV_BASE_GRADIENT_CALCULATOR_HH
30 
31 #include "fvbaseproperties.hh"
32 
33 #include <opm/common/Unused.hpp>
34 
35 #include <dune/common/fvector.hh>
36 
37 namespace Ewoms {
38 template<class TypeTag>
40 
47 template<class TypeTag>
49 {
50  typedef typename GET_PROP_TYPE(TypeTag, GridView) GridView;
51  typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
52  typedef typename GET_PROP_TYPE(TypeTag, Evaluation) Evaluation;
53  typedef typename GET_PROP_TYPE(TypeTag, Discretization) Discretization;
54  typedef typename GET_PROP_TYPE(TypeTag, ElementContext) ElementContext;
55 
56  enum { dim = GridView::dimension };
57  enum { dimWorld = GridView::dimensionworld };
58 
59  // maximum number of flux approximation points. to calculate this,
60  // we assume that the geometry with the most pointsq is a cube.
61  enum { maxFap = 2 << dim };
62 
63  typedef Opm::MathToolbox<Evaluation> Toolbox;
64  typedef Dune::FieldVector<Scalar, dim> DimVector;
65  typedef Dune::FieldVector<Evaluation, dim> EvalDimVector;
66 
67 public:
72  static void registerParameters()
73  { }
74 
82  template <bool prepareValues = true, bool prepareGradients = true>
83  void prepare(const ElementContext& elemCtx OPM_UNUSED, unsigned timeIdx OPM_UNUSED)
84  { /* noting to do */ }
85 
97  template <class QuantityCallback>
98  auto calculateScalarValue(const ElementContext& elemCtx,
99  unsigned fapIdx,
100  const QuantityCallback& quantityCallback) const
101  -> typename std::remove_reference<decltype(quantityCallback.operator()(0))>::type
102  {
103  typedef decltype(quantityCallback.operator()(0)) RawReturnType;
104  typedef typename std::remove_const<typename std::remove_reference<RawReturnType>::type>::type ReturnType;
105 
106  typedef Opm::MathToolbox<ReturnType> Toolbox;
107 
108  Scalar interiorDistance;
109  Scalar exteriorDistance;
110  computeDistances_(interiorDistance, exteriorDistance, elemCtx, fapIdx);
111 
112  const auto& face = elemCtx.stencil(/*timeIdx=*/0).interiorFace(fapIdx);
113  auto i = face.interiorIndex();
114  auto j = face.exteriorIndex();
115  auto focusDofIdx = elemCtx.focusDofIndex();
116 
117  // use the average weighted by distance...
118  ReturnType value;
119  if (i == focusDofIdx)
120  value = quantityCallback(i)*interiorDistance;
121  else
122  value = Toolbox::value(quantityCallback(i))*interiorDistance;
123 
124  if (j == focusDofIdx)
125  value += quantityCallback(j)*exteriorDistance;
126  else
127  value += Toolbox::value(quantityCallback(j))*exteriorDistance;
128 
129  value /= interiorDistance + exteriorDistance;
130 
131  return value;
132  }
133 
145  template <class QuantityCallback>
146  auto calculateVectorValue(const ElementContext& elemCtx,
147  unsigned fapIdx,
148  const QuantityCallback& quantityCallback) const
149  -> typename std::remove_reference<decltype(quantityCallback.operator()(0))>::type
150  {
151  typedef decltype(quantityCallback.operator()(0)) RawReturnType;
152  typedef typename std::remove_const<typename std::remove_reference<RawReturnType>::type>::type ReturnType;
153 
154  typedef decltype(std::declval<ReturnType>()[0]) RawFieldType;
155  typedef typename std::remove_const<typename std::remove_reference<RawFieldType>::type>::type FieldType;
156 
157  typedef Opm::MathToolbox<FieldType> Toolbox;
158 
159  Scalar interiorDistance;
160  Scalar exteriorDistance;
161  computeDistances_(interiorDistance, exteriorDistance, elemCtx, fapIdx);
162 
163  const auto& face = elemCtx.stencil(/*timeIdx=*/0).interiorFace(fapIdx);
164  auto i = face.interiorIndex();
165  auto j = face.exteriorIndex();
166  auto focusDofIdx = elemCtx.focusDofIndex();
167 
168  // use the average weighted by distance...
169  ReturnType value;
170  if (i == focusDofIdx) {
171  value = quantityCallback(i);
172  for (int k = 0; k < value.size(); ++k)
173  value[k] *= interiorDistance;
174  }
175  else {
176  const auto& dofVal = Toolbox::value(quantityCallback(i));
177  for (int k = 0; k < dofVal.size(); ++k)
178  value[k] = Toolbox::value(dofVal[k])*interiorDistance;
179  }
180 
181  if (j == focusDofIdx) {
182  const auto& dofVal = quantityCallback(j);
183  for (int k = 0; k < dofVal.size(); ++k)
184  value[k] += dofVal[k]*exteriorDistance;
185  }
186  else {
187  const auto& dofVal = quantityCallback(j);
188  for (int k = 0; k < dofVal.size(); ++k)
189  value[k] += Toolbox::value(dofVal[k])*exteriorDistance;
190  }
191 
192  Scalar totDistance = interiorDistance + exteriorDistance;
193  for (int k = 0; k < value.size(); ++k)
194  value[k] /= totDistance;
195 
196  return value;
197  }
198 
210  template <class QuantityCallback>
211  void calculateGradient(EvalDimVector& quantityGrad,
212  const ElementContext& elemCtx,
213  unsigned fapIdx,
214  const QuantityCallback& quantityCallback) const
215  {
216  const auto& stencil = elemCtx.stencil(/*timeIdx=*/0);
217  const auto& face = stencil.interiorFace(fapIdx);
218 
219  auto i = face.interiorIndex();
220  auto j = face.exteriorIndex();
221  auto focusIdx = elemCtx.focusDofIndex();
222 
223  const auto& interiorPos = stencil.subControlVolume(i).globalPos();
224  const auto& exteriorPos = stencil.subControlVolume(j).globalPos();
225 
226  Evaluation deltay;
227  if (i == focusIdx) {
228  deltay =
229  Toolbox::value(quantityCallback(j))
230  - quantityCallback(i);
231  }
232  else if (j == focusIdx) {
233  deltay =
234  quantityCallback(j)
235  - Toolbox::value(quantityCallback(i));
236  }
237  else
238  deltay =
239  Toolbox::value(quantityCallback(j))
240  - Toolbox::value(quantityCallback(i));
241 
242  Scalar distSquared = 0.0;
243  for (unsigned dimIdx = 0; dimIdx < dimWorld; ++dimIdx) {
244  Scalar tmp = exteriorPos[dimIdx] - interiorPos[dimIdx];
245  distSquared += tmp*tmp;
246  }
247 
248  // divide the gradient by the squared distance between the centers of the
249  // sub-control volumes: the gradient is the normalized directional vector between
250  // the two centers times the ratio of the difference of the values and their
251  // distance, i.e., d/abs(d) * delta y / abs(d) = d*delta y / abs(d)^2.
252  for (unsigned dimIdx = 0; dimIdx < dimWorld; ++dimIdx) {
253  Scalar tmp = exteriorPos[dimIdx] - interiorPos[dimIdx];
254  quantityGrad[dimIdx] = deltay*(tmp/distSquared);
255  }
256  }
257 
272  template <class QuantityCallback>
273  auto calculateBoundaryValue(const ElementContext& elemCtx OPM_UNUSED,
274  unsigned fapIdx OPM_UNUSED,
275  const QuantityCallback& quantityCallback)
276  -> decltype(quantityCallback.boundaryValue())
277  { return quantityCallback.boundaryValue(); }
278 
293  template <class QuantityCallback>
294  void calculateBoundaryGradient(EvalDimVector& quantityGrad,
295  const ElementContext& elemCtx,
296  unsigned faceIdx,
297  const QuantityCallback& quantityCallback) const
298  {
299  typedef typename std::decay<decltype(quantityCallback(0))>::type Evaluation;
300  typedef typename std::decay<decltype(quantityCallback.boundaryValue())>::type BoundaryEvaluation;
301 
302  typedef Opm::MathToolbox<Evaluation> Toolbox;
303  typedef Opm::MathToolbox<BoundaryEvaluation> BoundaryToolbox;
304 
305  const auto& stencil = elemCtx.stencil(/*timeIdx=*/0);
306  const auto& face = stencil.boundaryFace(faceIdx);
307 
308  Evaluation deltay;
309  if (face.interiorIndex() == elemCtx.focusDofIndex())
310  deltay = quantityCallback.boundaryValue() - quantityCallback(face.interiorIndex());
311  else
312  deltay =
313  BoundaryToolbox::value(quantityCallback.boundaryValue())
314  - Toolbox::value(quantityCallback(face.interiorIndex()));
315 
316  const auto& boundaryFacePos = face.integrationPos();
317  const auto& interiorPos = stencil.subControlVolume(face.interiorIndex()).center();
318 
319  Scalar distSquared = 0;
320  for (unsigned dimIdx = 0; dimIdx < dimWorld; ++dimIdx) {
321  Scalar tmp = boundaryFacePos[dimIdx] - interiorPos[dimIdx];
322  distSquared += tmp*tmp;
323  }
324 
325  // divide the gradient by the squared distance between the center of the
326  // sub-control and the center of the boundary face: the gradient is the
327  // normalized directional vector between the two centers times the ratio of the
328  // difference of the values and their distance, i.e., d/abs(d) * deltay / abs(d)
329  // = d*deltay / abs(d)^2.
330  for (unsigned dimIdx = 0; dimIdx < dimWorld; ++dimIdx) {
331  Scalar tmp = boundaryFacePos[dimIdx] - interiorPos[dimIdx];
332  quantityGrad[dimIdx] = deltay*(tmp/distSquared);
333  }
334  }
335 
336 private:
337  void computeDistances_(Scalar& interiorDistance,
338  Scalar& exteriorDistance,
339  const ElementContext& elemCtx,
340  unsigned fapIdx) const
341  {
342  const auto& stencil = elemCtx.stencil(/*timeIdx=*/0);
343  const auto& face = stencil.interiorFace(fapIdx);
344 
345  // calculate the distances of the position of the interior and of the exterior
346  // finite volume to the position of the integration point.
347  const auto& normal = face.normal();
348  auto i = face.interiorIndex();
349  auto j = face.exteriorIndex();
350  const auto& interiorPos = stencil.subControlVolume(i).globalPos();
351  const auto& exteriorPos = stencil.subControlVolume(j).globalPos();
352  const auto& integrationPos = face.integrationPos();
353 
354  interiorDistance = 0.0;
355  exteriorDistance = 0.0;
356  for (unsigned dimIdx = 0; dimIdx < dimWorld; ++dimIdx) {
357  interiorDistance +=
358  (interiorPos[dimIdx] - integrationPos[dimIdx])
359  * normal[dimIdx];
360 
361  exteriorDistance +=
362  (exteriorPos[dimIdx] - integrationPos[dimIdx])
363  * normal[dimIdx];
364  }
365 
366  interiorDistance = std::sqrt(std::abs(interiorDistance));
367  exteriorDistance = std::sqrt(std::abs(exteriorDistance));
368  }
369 };
370 } // namespace Ewoms
371 
372 #endif
Definition: baseauxiliarymodule.hh:37
void calculateBoundaryGradient(EvalDimVector &quantityGrad, const ElementContext &elemCtx, unsigned faceIdx, const QuantityCallback &quantityCallback) const
Calculates the gradient of an arbitrary quantity at any flux approximation point on the boundary...
Definition: fvbasegradientcalculator.hh:294
void calculateGradient(EvalDimVector &quantityGrad, const ElementContext &elemCtx, unsigned fapIdx, const QuantityCallback &quantityCallback) const
Calculates the gradient of an arbitrary quantity at any flux approximation point. ...
Definition: fvbasegradientcalculator.hh:211
static void registerParameters()
Register all run-time parameters for the gradient calculator of the base class of the discretization...
Definition: fvbasegradientcalculator.hh:72
The base class for the element-centered finite-volume discretization scheme.
Definition: fvbasegradientcalculator.hh:39
Declare the properties used by the infrastructure code of the finite volume discretizations.
This class calculates gradients of arbitrary quantities at flux integration points using the two-poin...
Definition: fvbasegradientcalculator.hh:48
void prepare(const ElementContext &elemCtx OPM_UNUSED, unsigned timeIdx OPM_UNUSED)
Precomputes the common values to calculate gradients and values of quantities at every interior flux ...
Definition: fvbasegradientcalculator.hh:83
auto calculateScalarValue(const ElementContext &elemCtx, unsigned fapIdx, const QuantityCallback &quantityCallback) const -> typename std::remove_reference< decltype(quantityCallback.operator()(0))>::type
Calculates the value of an arbitrary scalar quantity at any interior flux approximation point...
Definition: fvbasegradientcalculator.hh:98
auto calculateBoundaryValue(const ElementContext &elemCtx OPM_UNUSED, unsigned fapIdx OPM_UNUSED, const QuantityCallback &quantityCallback) -> decltype(quantityCallback.boundaryValue())
Calculates the value of an arbitrary quantity at any flux approximation point on the grid boundary...
Definition: fvbasegradientcalculator.hh:273
auto calculateVectorValue(const ElementContext &elemCtx, unsigned fapIdx, const QuantityCallback &quantityCallback) const -> typename std::remove_reference< decltype(quantityCallback.operator()(0))>::type
Calculates the value of an arbitrary vectorial quantity at any interior flux approximation point...
Definition: fvbasegradientcalculator.hh:146