fvbaseadlocallinearizer.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_AD_LOCAL_LINEARIZER_HH
29 #define EWOMS_FV_BASE_AD_LOCAL_LINEARIZER_HH
30 
31 #include "fvbaseproperties.hh"
32 
33 #include <opm/material/densead/Math.hpp>
34 #include <opm/common/Valgrind.hpp>
35 #include <opm/common/Unused.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 namespace Ewoms {
44 // forward declaration
45 template<class TypeTag>
47 
48 namespace Properties {
49 // declare the property tags required for the finite differences local linearizer
50 NEW_TYPE_TAG(AutoDiffLocalLinearizer);
51 
52 NEW_PROP_TAG(LocalLinearizer);
53 NEW_PROP_TAG(Evaluation);
54 
55 NEW_PROP_TAG(LocalResidual);
57 NEW_PROP_TAG(Problem);
58 NEW_PROP_TAG(Model);
59 NEW_PROP_TAG(PrimaryVariables);
60 NEW_PROP_TAG(ElementContext);
61 NEW_PROP_TAG(Scalar);
62 NEW_PROP_TAG(Evaluation);
63 NEW_PROP_TAG(GridView);
64 
65 // set the properties to be spliced in
66 SET_TYPE_PROP(AutoDiffLocalLinearizer, LocalLinearizer,
68 
70 SET_PROP(AutoDiffLocalLinearizer, Evaluation)
71 {
72 private:
73  static const unsigned numEq = GET_PROP_VALUE(TypeTag, NumEq);
74 
75  typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
76 
77 public:
78  typedef Opm::DenseAd::Evaluation<Scalar, numEq> type;
79 };
80 
81 } // namespace Properties
82 
91 template<class TypeTag>
92 class FvBaseAdLocalLinearizer
93 {
94 private:
95  typedef typename GET_PROP_TYPE(TypeTag, LocalLinearizer) Implementation;
96  typedef typename GET_PROP_TYPE(TypeTag, LocalResidual) LocalResidual;
97  typedef typename GET_PROP_TYPE(TypeTag, Simulator) Simulator;
98  typedef typename GET_PROP_TYPE(TypeTag, Problem) Problem;
99  typedef typename GET_PROP_TYPE(TypeTag, Model) Model;
100  typedef typename GET_PROP_TYPE(TypeTag, PrimaryVariables) PrimaryVariables;
101  typedef typename GET_PROP_TYPE(TypeTag, ElementContext) ElementContext;
102  typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
103  typedef typename GET_PROP_TYPE(TypeTag, GridView) GridView;
104  typedef typename GridView::template Codim<0>::Entity Element;
105 
106  enum { numEq = GET_PROP_VALUE(TypeTag, NumEq) };
107 
108  typedef Dune::FieldVector<Scalar, numEq> ScalarVectorBlock;
109  typedef Dune::FieldMatrix<Scalar, numEq, numEq> ScalarMatrixBlock;
110 
111  typedef Dune::BlockVector<ScalarVectorBlock> ScalarLocalBlockVector;
112  typedef Dune::Matrix<ScalarMatrixBlock> ScalarLocalBlockMatrix;
113 
114 public:
115  FvBaseAdLocalLinearizer()
116  : internalElemContext_(0)
117  { }
118 
119  // copying local linearizer objects around is a very bad idea, so we explicitly
120  // prevent it...
121  FvBaseAdLocalLinearizer(const FvBaseAdLocalLinearizer&) = delete;
122 
123  ~FvBaseAdLocalLinearizer()
124  { delete internalElemContext_; }
125 
129  static void registerParameters()
130  { }
131 
140  void init(Simulator& simulator)
141  {
142  simulatorPtr_ = &simulator;
143  delete internalElemContext_;
144  internalElemContext_ = new ElementContext(simulator);
145  }
146 
158  void linearize(const Element& element)
159  {
160  linearize(*internalElemContext_, element);
161  }
162 
177  void linearize(ElementContext& elemCtx, const Element& elem)
178  {
179  elemCtx.updateStencil(elem);
180  elemCtx.updateAllIntensiveQuantities();
181 
182  // update the weights of the primary variables for the context
183  model_().updatePVWeights(elemCtx);
184 
185  // resize the internal arrays of the linearizer
186  resize_(elemCtx);
187  reset_(elemCtx);
188 
189  // compute the local residual and its Jacobian
190  unsigned numPrimaryDof = elemCtx.numPrimaryDof(/*timeIdx=*/0);
191  for (unsigned focusDofIdx = 0; focusDofIdx < numPrimaryDof; focusDofIdx++) {
192  elemCtx.setFocusDofIndex(focusDofIdx);
193  elemCtx.updateAllExtensiveQuantities();
194 
195  // calculate the local residual
196  localResidual_.eval(elemCtx);
197 
198  // convert the local Jacobian matrix and the right hand side from the data
199  // structures used by the automatic differentiation code to the conventional
200  // ones used by the linear solver.
201  updateLocalLinearization_(elemCtx, focusDofIdx);
202  }
203  }
204 
208  LocalResidual& localResidual()
209  { return localResidual_; }
210 
214  const LocalResidual& localResidual() const
215  { return localResidual_; }
216 
225  const ScalarMatrixBlock& jacobian(unsigned domainScvIdx, unsigned rangeScvIdx) const
226  { return jacobian_[domainScvIdx][rangeScvIdx]; }
227 
233  const ScalarVectorBlock& residual(unsigned dofIdx) const
234  { return residual_[dofIdx]; }
235 
236 protected:
237  Implementation& asImp_()
238  { return *static_cast<Implementation*>(this); }
239  const Implementation& asImp_() const
240  { return *static_cast<const Implementation*>(this); }
241 
242  const Simulator& simulator_() const
243  { return *simulatorPtr_; }
244  const Problem& problem_() const
245  { return simulatorPtr_->problem(); }
246  const Model& model_() const
247  { return simulatorPtr_->model(); }
248 
253  void resize_(const ElementContext& elemCtx)
254  {
255  size_t numDof = elemCtx.numDof(/*timeIdx=*/0);
256  size_t numPrimaryDof = elemCtx.numPrimaryDof(/*timeIdx=*/0);
257 
258  residual_.resize(numDof);
259  jacobian_.setSize(numDof, numPrimaryDof);
260  }
261 
265  void reset_(const ElementContext& elemCtx OPM_UNUSED)
266  {
267  residual_ = 0.0;
268  jacobian_ = 0.0;
269  }
270 
275  void updateLocalLinearization_(const ElementContext& elemCtx,
276  unsigned focusDofIdx)
277  {
278  const auto& resid = localResidual_.residual();
279 
280  for (unsigned eqIdx = 0; eqIdx < numEq; eqIdx++)
281  residual_[focusDofIdx][eqIdx] = resid[focusDofIdx][eqIdx].value();
282 
283  size_t numDof = elemCtx.numDof(/*timeIdx=*/0);
284  for (unsigned dofIdx = 0; dofIdx < numDof; dofIdx++) {
285  for (unsigned eqIdx = 0; eqIdx < numEq; eqIdx++) {
286  for (unsigned pvIdx = 0; pvIdx < numEq; pvIdx++) {
287  // A[dofIdx][focusDofIdx][eqIdx][pvIdx] is the partial derivative of
288  // the residual function 'eqIdx' for the degree of freedom 'dofIdx'
289  // with regard to the focus variable 'pvIdx' of the degree of freedom
290  // 'focusDofIdx'
291  jacobian_[dofIdx][focusDofIdx][eqIdx][pvIdx] = resid[dofIdx][eqIdx].derivative(pvIdx);
292  Opm::Valgrind::CheckDefined(jacobian_[dofIdx][focusDofIdx][eqIdx][pvIdx]);
293  }
294  }
295  }
296  }
297 
298  Simulator *simulatorPtr_;
299  Model *modelPtr_;
300 
301  ElementContext *internalElemContext_;
302 
303  LocalResidual localResidual_;
304 
305  ScalarLocalBlockVector residual_;
306  ScalarLocalBlockMatrix jacobian_;
307 };
308 
309 } // namespace Ewoms
310 
311 #endif
void linearize(const Element &element)
Compute an element&#39;s local Jacobian matrix and evaluate its residual.
Definition: fvbaseadlocallinearizer.hh:158
Definition: baseauxiliarymodule.hh:37
void linearize(ElementContext &elemCtx, const Element &elem)
Compute an element&#39;s local Jacobian matrix and evaluate its residual.
Definition: fvbaseadlocallinearizer.hh:177
Model & model()
Return the physical model used in the simulation.
Definition: simulator.hh:189
void reset_(const ElementContext &elemCtx OPM_UNUSED)
Reset the all relevant internal attributes to 0.
Definition: fvbaseadlocallinearizer.hh:265
#define GET_PROP_VALUE(TypeTag, PropTagName)
Access the value attribute of a property for a type tag.
Definition: propertysystem.hh:469
const ScalarVectorBlock & residual(unsigned dofIdx) const
Returns the local residual of a sub-control volume.
Definition: fvbaseadlocallinearizer.hh:233
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
const LocalResidual & localResidual() const
Return reference to the local residual.
Definition: fvbaseadlocallinearizer.hh:214
const ScalarMatrixBlock & jacobian(unsigned domainScvIdx, unsigned rangeScvIdx) const
Returns the local Jacobian matrix of the residual of a sub-control volume.
Definition: fvbaseadlocallinearizer.hh:225
void resize_(const ElementContext &elemCtx)
Resize all internal attributes to the size of the element.
Definition: fvbaseadlocallinearizer.hh:253
Declare the properties used by the infrastructure code of the finite volume discretizations.
#define SET_PROP(EffTypeTagName, PropTagName)
Set a property for a specific type tag.
Definition: propertysystem.hh:297
Calculates the local residual and its Jacobian for a single element of the grid.
Definition: fvbaseadlocallinearizer.hh:46
void updateLocalLinearization_(const ElementContext &elemCtx, unsigned focusDofIdx)
Updates the current local Jacobian matrix with the partial derivatives of all equations for the degre...
Definition: fvbaseadlocallinearizer.hh:275
static void registerParameters()
Register all run-time parameters for the local jacobian.
Definition: fvbaseadlocallinearizer.hh:129
void init(Simulator &simulator)
Initialize the local Jacobian object.
Definition: fvbaseadlocallinearizer.hh:140
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
LocalResidual & localResidual()
Return reference to the local residual.
Definition: fvbaseadlocallinearizer.hh:208