26 #ifndef EWOMS_FV_BASE_GRADIENT_CALCULATOR_HH
27 #define EWOMS_FV_BASE_GRADIENT_CALCULATOR_HH
31 #include <dune/common/fvector.hh>
34 template<
class TypeTag>
43 template<
class TypeTag>
47 typedef typename GET_PROP_TYPE(TypeTag, GridView) GridView;
48 typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
49 typedef typename GET_PROP_TYPE(TypeTag, Evaluation) Evaluation;
50 typedef typename GET_PROP_TYPE(TypeTag, Discretization) Discretization;
51 typedef typename GET_PROP_TYPE(TypeTag, ElementContext) ElementContext;
53 enum { dim = GridView::dimension };
54 enum { dimWorld = GridView::dimensionworld };
58 enum { maxFap = 2 << dim };
60 typedef Opm::MathToolbox<Evaluation> Toolbox;
61 typedef Dune::FieldVector<Scalar, dim> DimVector;
62 typedef Dune::FieldVector<Evaluation, dim> EvalDimVector;
64 static_assert(std::is_same<Evaluation, Scalar>::value ||
66 "So far, automatic differentiation is only available for the "
67 "element-centered finite volume discretization!");
74 static void registerParameters()
84 template <
bool prepareValues = true,
bool prepareGradients = true>
85 void prepare(
const ElementContext &elemCtx,
int timeIdx)
87 const auto &stencil = elemCtx.stencil(timeIdx);
88 for (
int fapIdx = 0; fapIdx < stencil.numInteriorFaces(); ++ fapIdx) {
89 const auto &scvf = stencil.interiorFace(fapIdx);
90 const auto &normal = scvf.normal();
91 const auto &interiorPos = stencil.subControlVolume(scvf.interiorIndex()).globalPos();
92 const auto &exteriorPos = stencil.subControlVolume(scvf.exteriorIndex()).globalPos();
94 interiorDistance_[fapIdx] = 0;
95 exteriorDistance_[fapIdx] = 0;
96 for (
int dimIdx = 0; dimIdx < dimWorld; ++dimIdx) {
97 interiorDistance_[fapIdx] +=
98 (interiorPos[dimIdx] - scvf.integrationPos()[dimIdx])
101 exteriorDistance_[fapIdx] +=
102 (exteriorPos[dimIdx] - scvf.integrationPos()[dimIdx])
106 interiorDistance_[fapIdx] = std::abs(interiorDistance_[fapIdx]);
107 exteriorDistance_[fapIdx] = std::abs(exteriorDistance_[fapIdx]);
123 template <
class QuantityCallback>
124 auto calculateValue(
const ElementContext &elemCtx,
126 const QuantityCallback &quantityCallback)
const
127 ->
typename std::remove_reference<decltype(quantityCallback.operator()(0))>::type
129 const auto &face = elemCtx.stencil(0).interiorFace(fapIdx);
132 auto value(quantityCallback(face.interiorIndex()));
133 value *= interiorDistance_[fapIdx];
134 auto tmp(quantityCallback(face.exteriorIndex()));
135 tmp *= exteriorDistance_[fapIdx];
137 value /= interiorDistance_[fapIdx] + exteriorDistance_[fapIdx];
153 template <
class QuantityCallback>
154 void calculateGradient(EvalDimVector &quantityGrad,
155 const ElementContext &elemCtx,
157 const QuantityCallback &quantityCallback)
const
159 const auto &stencil = elemCtx.stencil(0);
160 const auto &face = stencil.interiorFace(fapIdx);
162 const auto& exteriorPos = stencil.subControlVolume(face.exteriorIndex()).center();
163 const auto& interiorPos = stencil.subControlVolume(face.interiorIndex()).center();
173 Toolbox::value(quantityCallback(face.exteriorIndex()))
174 - quantityCallback(face.interiorIndex());
176 Scalar distSquared = 0;
177 for (
int dimIdx = 0; dimIdx < dimWorld; ++dimIdx) {
178 quantityGrad[dimIdx] = deltay;
180 Scalar tmp = exteriorPos[dimIdx] - interiorPos[dimIdx];
181 distSquared += tmp*tmp;
188 for (
int dimIdx = 0; dimIdx < dimWorld; ++dimIdx) {
189 Scalar tmp = exteriorPos[dimIdx] - interiorPos[dimIdx];
190 quantityGrad[dimIdx] *= tmp/distSquared;
208 template <
class QuantityCallback>
209 auto calculateBoundaryValue(
const ElementContext &elemCtx,
211 const QuantityCallback &quantityCallback)
212 -> decltype(quantityCallback.boundaryValue())
213 {
return quantityCallback.boundaryValue(); }
229 template <
class QuantityCallback>
230 void calculateBoundaryGradient(EvalDimVector &quantityGrad,
231 const ElementContext &elemCtx,
233 const QuantityCallback &quantityCallback)
const
235 const auto &stencil = elemCtx.stencil(0);
236 const auto &face = stencil.boundaryFace(faceIdx);
238 auto deltay = quantityCallback.boundaryValue() - quantityCallback(face.interiorIndex());
240 const auto& boundaryFacePos = face.integrationPos();
241 const auto& interiorPos = stencil.subControlVolume(face.interiorIndex()).center();
243 Scalar distSquared = 0;
244 for (
int dimIdx = 0; dimIdx < dimWorld; ++dimIdx) {
245 Scalar tmp = boundaryFacePos[dimIdx] - interiorPos[dimIdx];
246 distSquared += tmp*tmp;
248 quantityGrad[dimIdx] = deltay;
256 for (
int dimIdx = 0; dimIdx < dimWorld; ++dimIdx) {
257 Scalar tmp = boundaryFacePos[dimIdx] - interiorPos[dimIdx];
258 quantityGrad[dimIdx] *= tmp/distSquared;
265 Scalar interiorDistance_[maxFap];
269 Scalar exteriorDistance_[maxFap];
Declare the properties used by the infrastructure code of the finite volume discretizations.
The base class for the element-centered finite-volume discretization scheme.
Definition: fvbasegradientcalculator.hh:35
Definition: baseauxiliarymodule.hh:35
This class calculates gradients of arbitrary quantities at flux integration points using the two-poin...
Definition: fvbasegradientcalculator.hh:44