fvbasefdlocallinearizer.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_FD_LOCAL_LINEARIZER_HH
29 #define EWOMS_FV_BASE_FD_LOCAL_LINEARIZER_HH
30 
33 
34 #include <opm/material/common/MathToolbox.hpp>
35 #include <opm/common/Valgrind.hpp>
36 
37 #include <dune/istl/bvector.hh>
38 #include <dune/istl/matrix.hh>
39 
40 #include <dune/common/fvector.hh>
41 #include <dune/common/fmatrix.hh>
42 
43 #include <limits>
44 
45 namespace Ewoms {
46 // forward declaration
47 template<class TypeTag>
49 
50 namespace Properties {
51 // declare the property tags required for the finite differences local linearizer
52 NEW_TYPE_TAG(FiniteDifferenceLocalLinearizer);
53 
54 NEW_PROP_TAG(LocalLinearizer);
55 NEW_PROP_TAG(Evaluation);
56 NEW_PROP_TAG(NumericDifferenceMethod);
57 NEW_PROP_TAG(BaseEpsilon);
58 
59 NEW_PROP_TAG(LocalResidual);
61 NEW_PROP_TAG(Problem);
62 NEW_PROP_TAG(Model);
63 NEW_PROP_TAG(PrimaryVariables);
64 NEW_PROP_TAG(ElementContext);
65 NEW_PROP_TAG(Scalar);
66 NEW_PROP_TAG(Evaluation);
67 NEW_PROP_TAG(GridView);
68 NEW_PROP_TAG(NumEq);
69 
70 // set the properties to be spliced in
71 SET_TYPE_PROP(FiniteDifferenceLocalLinearizer, LocalLinearizer,
73 
74 SET_TYPE_PROP(FiniteDifferenceLocalLinearizer, Evaluation,
75  typename GET_PROP_TYPE(TypeTag, Scalar));
76 
84 SET_INT_PROP(FiniteDifferenceLocalLinearizer, NumericDifferenceMethod, +1);
85 
87 SET_SCALAR_PROP(FiniteDifferenceLocalLinearizer,
88  BaseEpsilon,
89  std::max<Scalar>(0.9123e-10, std::numeric_limits<Scalar>::epsilon()*1.23e3));
90 }
91 
128 template<class TypeTag>
130 {
131 private:
132  typedef typename GET_PROP_TYPE(TypeTag, LocalLinearizer) Implementation;
133  typedef typename GET_PROP_TYPE(TypeTag, LocalResidual) LocalResidual;
134  typedef typename GET_PROP_TYPE(TypeTag, Simulator) Simulator;
135  typedef typename GET_PROP_TYPE(TypeTag, Problem) Problem;
136  typedef typename GET_PROP_TYPE(TypeTag, Model) Model;
137  typedef typename GET_PROP_TYPE(TypeTag, PrimaryVariables) PrimaryVariables;
138  typedef typename GET_PROP_TYPE(TypeTag, ElementContext) ElementContext;
139  typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
140  typedef typename GET_PROP_TYPE(TypeTag, GridView) GridView;
141  typedef typename GridView::template Codim<0>::Entity Element;
142 
143  enum { numEq = GET_PROP_VALUE(TypeTag, NumEq) };
144  typedef Dune::FieldMatrix<Scalar, numEq, numEq> ScalarMatrixBlock;
145  typedef Dune::FieldVector<Scalar, numEq> ScalarVectorBlock;
146 
147  typedef Dune::BlockVector<ScalarVectorBlock> ScalarLocalBlockVector;
148  typedef Dune::Matrix<ScalarMatrixBlock> ScalarLocalBlockMatrix;
149 
150  typedef typename LocalResidual::LocalEvalBlockVector LocalEvalBlockVector;
151 
152 #if __GNUC__ == 4 && __GNUC_MINOR__ <= 6
153 public:
154  // make older GCCs happy by providing a public copy constructor (this is necessary
155  // for their implementation of std::vector, although the method is never called...)
156  FvBaseFdLocalLinearizer(const FvBaseFdLocalLinearizer&)
157  : internalElemContext_(0)
158  {}
159 
160 #else
161  // copying local residual objects around is a very bad idea, so we explicitly prevent
162  // it...
163  FvBaseFdLocalLinearizer(const FvBaseFdLocalLinearizer&) = delete;
164 #endif
165 public:
166  FvBaseFdLocalLinearizer()
167  : internalElemContext_(0)
168  { }
169 
170  ~FvBaseFdLocalLinearizer()
171  { delete internalElemContext_; }
172 
176  static void registerParameters()
177  {
178  EWOMS_REGISTER_PARAM(TypeTag, int, NumericDifferenceMethod,
179  "The method used for numeric differentiation (-1: backward "
180  "differences, 0: central differences, 1: forward differences)");
181  }
182 
191  void init(Simulator& simulator)
192  {
193  simulatorPtr_ = &simulator;
194  delete internalElemContext_;
195  internalElemContext_ = new ElementContext(simulator);
196  }
197 
209  void linearize(const Element& element)
210  {
211  linearize(*internalElemContext_, element);
212  }
213 
228  void linearize(ElementContext& elemCtx, const Element& elem)
229  {
230  elemCtx.updateAll(elem);
231 
232  // update the weights of the primary variables for the context
233  model_().updatePVWeights(elemCtx);
234 
235  resize_(elemCtx);
236  reset_(elemCtx);
237 
238  // calculate the local residual
239  localResidual_.eval(residual_, elemCtx);
240 
241  // calculate the local jacobian matrix
242  size_t numPrimaryDof = elemCtx.numPrimaryDof(/*timeIdx=*/0);
243  for (unsigned dofIdx = 0; dofIdx < numPrimaryDof; dofIdx++) {
244  for (unsigned pvIdx = 0; pvIdx < numEq; pvIdx++) {
245  asImp_().evalPartialDerivative_(elemCtx, dofIdx, pvIdx);
246 
247  // incorporate the partial derivatives into the local Jacobian matrix
248  updateLocalJacobian_(elemCtx, dofIdx, pvIdx);
249  }
250  }
251  }
252 
257  static Scalar baseEpsilon()
258  { return GET_PROP_VALUE(TypeTag, BaseEpsilon); }
259 
271  Scalar numericEpsilon(const ElementContext& elemCtx,
272  unsigned dofIdx,
273  unsigned pvIdx) const
274  {
275  unsigned globalIdx = elemCtx.globalSpaceIndex(dofIdx, /*timeIdx=*/0);
276  Scalar pvWeight = elemCtx.model().primaryVarWeight(globalIdx, pvIdx);
277  assert(pvWeight > 0 && std::isfinite(pvWeight));
278  Opm::Valgrind::CheckDefined(pvWeight);
279 
280  return baseEpsilon()/pvWeight;
281  }
282 
286  LocalResidual& localResidual()
287  { return localResidual_; }
288 
292  const LocalResidual& localResidual() const
293  { return localResidual_; }
294 
303  const ScalarMatrixBlock& jacobian(unsigned domainScvIdx, unsigned rangeScvIdx) const
304  { return jacobian_[domainScvIdx][rangeScvIdx]; }
305 
311  const ScalarVectorBlock& residual(unsigned dofIdx) const
312  { return residual_[dofIdx]; }
313 
314 protected:
315  Implementation& asImp_()
316  { return *static_cast<Implementation*>(this); }
317  const Implementation& asImp_() const
318  { return *static_cast<const Implementation*>(this); }
319 
320  const Simulator& simulator_() const
321  { return *simulatorPtr_; }
322  const Problem& problem_() const
323  { return simulatorPtr_->problem(); }
324  const Model& model_() const
325  { return simulatorPtr_->model(); }
326 
331  { return EWOMS_GET_PARAM(TypeTag, int, NumericDifferenceMethod); }
332 
337  void resize_(const ElementContext& elemCtx)
338  {
339  size_t numDof = elemCtx.numDof(/*timeIdx=*/0);
340  size_t numPrimaryDof = elemCtx.numPrimaryDof(/*timeIdx=*/0);
341 
342  residual_.resize(numDof);
343  jacobian_.setSize(numDof, numPrimaryDof);
344 
345  derivResidual_.resize(numDof);
346  }
347 
351  void reset_(const ElementContext& elemCtx)
352  {
353  size_t numDof = elemCtx.numDof(/*timeIdx=*/0);
354  size_t numPrimaryDof = elemCtx.numPrimaryDof(/*timeIdx=*/0);
355  for (unsigned primaryDofIdx = 0; primaryDofIdx < numPrimaryDof; ++ primaryDofIdx)
356  for (unsigned dof2Idx = 0; dof2Idx < numDof; ++ dof2Idx)
357  jacobian_[dof2Idx][primaryDofIdx] = 0.0;
358 
359  for (unsigned primaryDofIdx = 0; primaryDofIdx < numDof; ++ primaryDofIdx)
360  residual_[primaryDofIdx] = 0.0;
361  }
362 
405  void evalPartialDerivative_(ElementContext& elemCtx,
406  unsigned dofIdx,
407  unsigned pvIdx)
408  {
409  // save all quantities which depend on the specified primary
410  // variable at the given sub control volume
411  elemCtx.stashIntensiveQuantities(dofIdx);
412 
413  PrimaryVariables priVars(elemCtx.primaryVars(dofIdx, /*timeIdx=*/0));
414  Scalar eps = asImp_().numericEpsilon(elemCtx, dofIdx, pvIdx);
415  Scalar delta = 0.0;
416 
417  if (numericDifferenceMethod_() >= 0) {
418  // we are not using backward differences, i.e. we need to
419  // calculate f(x + \epsilon)
420 
421  // deflect primary variables
422  priVars[pvIdx] += eps;
423  delta += eps;
424 
425  // calculate the deflected residual
426  elemCtx.updateIntensiveQuantities(priVars, dofIdx, /*timeIdx=*/0);
427  elemCtx.updateAllExtensiveQuantities();
428  localResidual_.eval(derivResidual_, elemCtx);
429  }
430  else {
431  // we are using backward differences, i.e. we don't need
432  // to calculate f(x + \epsilon) and we can recycle the
433  // (already calculated) residual f(x)
434  derivResidual_ = residual_;
435  }
436 
437  if (numericDifferenceMethod_() <= 0) {
438  // we are not using forward differences, i.e. we don't
439  // need to calculate f(x - \epsilon)
440 
441  // deflect the primary variables
442  priVars[pvIdx] -= delta + eps;
443  delta += eps;
444 
445  // calculate the deflected residual again, this time we use the local
446  // residual's internal storage.
447  elemCtx.updateIntensiveQuantities(priVars, dofIdx, /*timeIdx=*/0);
448  elemCtx.updateAllExtensiveQuantities();
449  localResidual_.eval(elemCtx);
450 
451  derivResidual_ -= localResidual_.residual();
452  }
453  else {
454  // we are using forward differences, i.e. we don't need to
455  // calculate f(x - \epsilon) and we can recycle the
456  // (already calculated) residual f(x)
457  derivResidual_ -= residual_;
458  }
459 
460  assert(delta > 0);
461 
462  // divide difference in residuals by the magnitude of the
463  // deflections between the two function evaluation
464  derivResidual_ /= delta;
465 
466  // restore the original state of the element's volume
467  // variables
468  elemCtx.restoreIntensiveQuantities(dofIdx);
469 
470 #ifndef NDEBUG
471  for (unsigned i = 0; i < derivResidual_.size(); ++i)
472  Opm::Valgrind::CheckDefined(derivResidual_[i]);
473 #endif
474  }
475 
481  void updateLocalJacobian_(const ElementContext& elemCtx,
482  unsigned focusDofIdx,
483  unsigned pvIdx)
484  {
485  size_t numDof = elemCtx.numDof(/*timeIdx=*/0);
486  for (unsigned dofIdx = 0; dofIdx < numDof; dofIdx++) {
487  for (unsigned eqIdx = 0; eqIdx < numEq; eqIdx++) {
488  // A[dofIdx][focusDofIdx][eqIdx][pvIdx] is the partial derivative of the
489  // residual function 'eqIdx' for the degree of freedom 'dofIdx' with
490  // regard to the primary variable 'pvIdx' of the degree of freedom
491  // 'focusDofIdx'
492  jacobian_[dofIdx][focusDofIdx][eqIdx][pvIdx] = derivResidual_[dofIdx][eqIdx];
493  Opm::Valgrind::CheckDefined(jacobian_[dofIdx][focusDofIdx][eqIdx][pvIdx]);
494  }
495  }
496  }
497 
498  Simulator *simulatorPtr_;
499  Model *modelPtr_;
500 
501  ElementContext *internalElemContext_;
502 
503  LocalEvalBlockVector residual_;
504  LocalEvalBlockVector derivResidual_;
505  ScalarLocalBlockMatrix jacobian_;
506 
507  LocalResidual localResidual_;
508 };
509 
510 } // namespace Ewoms
511 
512 #endif
Scalar numericEpsilon(const ElementContext &elemCtx, unsigned dofIdx, unsigned pvIdx) const
Returns the epsilon value which is added and removed from the current solution.
Definition: fvbasefdlocallinearizer.hh:271
static void registerParameters()
Register all run-time parameters for the local jacobian.
Definition: fvbasefdlocallinearizer.hh:176
void linearize(ElementContext &elemCtx, const Element &elem)
Compute an element&#39;s local Jacobian matrix and evaluate its residual.
Definition: fvbasefdlocallinearizer.hh:228
void evalPartialDerivative_(ElementContext &elemCtx, unsigned dofIdx, unsigned pvIdx)
Compute the partial derivatives of a context&#39;s residual functions.
Definition: fvbasefdlocallinearizer.hh:405
Definition: baseauxiliarymodule.hh:37
Model & model()
Return the physical model used in the simulation.
Definition: simulator.hh:189
#define GET_PROP_VALUE(TypeTag, PropTagName)
Access the value attribute of a property for a type tag.
Definition: propertysystem.hh:469
static Scalar baseEpsilon()
Returns the unweighted epsilon value used to calculate the local derivatives.
Definition: fvbasefdlocallinearizer.hh:257
const LocalResidual & localResidual() const
Return reference to the local residual.
Definition: fvbasefdlocallinearizer.hh:292
Problem & problem()
Return the object which specifies the pysical setup of the simulation.
Definition: simulator.hh:202
#define NEW_TYPE_TAG(...)
Define a new type tag.
Definition: propertysystem.hh:169
#define GET_PROP_TYPE(TypeTag, PropTagName)
Access the type attribute of a property for a type tag.
Definition: propertysystem.hh:486
Calculates the Jacobian of the local residual for finite volume spatial discretizations using a finit...
Definition: fvbasefdlocallinearizer.hh:48
#define SET_INT_PROP(EffTypeTagName, PropTagName,...)
Set a property to a simple constant integer value.
Definition: propertysystem.hh:345
void updateLocalJacobian_(const ElementContext &elemCtx, unsigned focusDofIdx, unsigned pvIdx)
Updates the current local Jacobian matrix with the partial derivatives of all equations for primary v...
Definition: fvbasefdlocallinearizer.hh:481
#define EWOMS_REGISTER_PARAM(TypeTag, ParamType, ParamName, Description)
Register a run-time parameter.
Definition: parametersystem.hh:68
This file provides the infrastructure to retrieve run-time parameters.
static int numericDifferenceMethod_()
Returns the numeric difference method which is applied.
Definition: fvbasefdlocallinearizer.hh:330
void resize_(const ElementContext &elemCtx)
Resize all internal attributes to the size of the element.
Definition: fvbasefdlocallinearizer.hh:337
#define EWOMS_GET_PARAM(TypeTag, ParamType, ParamName)
Retrieve a runtime parameter.
Definition: parametersystem.hh:99
void linearize(const Element &element)
Compute an element&#39;s local Jacobian matrix and evaluate its residual.
Definition: fvbasefdlocallinearizer.hh:209
const ScalarVectorBlock & residual(unsigned dofIdx) const
Returns the local residual of a sub-control volume.
Definition: fvbasefdlocallinearizer.hh:311
LocalResidual & localResidual()
Return reference to the local residual.
Definition: fvbasefdlocallinearizer.hh:286
void init(Simulator &simulator)
Initialize the local Jacobian object.
Definition: fvbasefdlocallinearizer.hh:191
void reset_(const ElementContext &elemCtx)
Reset the all relevant internal attributes to 0.
Definition: fvbasefdlocallinearizer.hh:351
Provides the magic behind the eWoms property system.
Manages the initializing and running of time dependent problems.
Definition: simulator.hh:75
#define NEW_PROP_TAG(PTagName)
Define a property tag.
Definition: propertysystem.hh:247
#define SET_TYPE_PROP(EffTypeTagName, PropTagName,...)
Set a property which defines a type.
Definition: propertysystem.hh:377
const ScalarMatrixBlock & jacobian(unsigned domainScvIdx, unsigned rangeScvIdx) const
Returns the local Jacobian matrix of the residual of a sub-control volume.
Definition: fvbasefdlocallinearizer.hh:303
#define SET_SCALAR_PROP(EffTypeTagName, PropTagName,...)
Set a property to a simple constant scalar value.
Definition: propertysystem.hh:394