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/material/common/Valgrind.hpp>
35
36#include <dune/istl/bvector.hh>
37#include <dune/istl/matrix.hh>
38
39#include <dune/common/fvector.hh>
40#include <dune/common/fmatrix.hh>
41
42namespace Opm {
43// forward declaration
44template<class TypeTag>
45class FvBaseAdLocalLinearizer;
46}
47
48namespace Opm::Properties {
49
50// declare the property tags required for the finite differences local linearizer
51
52namespace TTag {
54} // namespace TTag
55
56// set the properties to be spliced in
57template<class TypeTag>
58struct LocalLinearizer<TypeTag, TTag::AutoDiffLocalLinearizer>
60
62template<class TypeTag>
63struct Evaluation<TypeTag, TTag::AutoDiffLocalLinearizer>
64{
65private:
66 static const unsigned numEq = getPropValue<TypeTag, Properties::NumEq>();
67
69
70public:
71 using type = DenseAd::Evaluation<Scalar, numEq>;
72};
73
74} // namespace Opm::Properties
75
76namespace Opm {
77
86template<class TypeTag>
88{
89private:
99 using Element = typename GridView::template Codim<0>::Entity;
100
101 enum { numEq = getPropValue<TypeTag, Properties::NumEq>() };
102
103 using ScalarVectorBlock = Dune::FieldVector<Scalar, numEq>;
104 // extract local matrices from jacobian matrix for consistency
106
107 using ScalarLocalBlockVector = Dune::BlockVector<ScalarVectorBlock>;
108 using ScalarLocalBlockMatrix = Dune::Matrix<ScalarMatrixBlock>;
109
110public:
113 { }
114
115 // copying local linearizer objects around is a very bad idea, so we explicitly
116 // prevent it...
118
120 { delete internalElemContext_; }
121
125 static void registerParameters()
126 { }
127
136 void init(Simulator& simulator)
137 {
138 simulatorPtr_ = &simulator;
140 internalElemContext_ = new ElementContext(simulator);
141 }
142
154 void linearize(const Element& element)
155 {
157 }
158
173 void linearize(ElementContext& elemCtx, const Element& elem)
174 {
175 elemCtx.updateStencil(elem);
176 elemCtx.updateAllIntensiveQuantities();
177
178 // update the weights of the primary variables for the context
179 model_().updatePVWeights(elemCtx);
180
181 // resize the internal arrays of the linearizer
182 resize_(elemCtx);
183
184 // compute the local residual and its Jacobian
185 unsigned numPrimaryDof = elemCtx.numPrimaryDof(/*timeIdx=*/0);
186 for (unsigned focusDofIdx = 0; focusDofIdx < numPrimaryDof; focusDofIdx++) {
187 elemCtx.setFocusDofIndex(focusDofIdx);
188 elemCtx.updateAllExtensiveQuantities();
189
190 // calculate the local residual
191 localResidual_.eval(elemCtx);
192
193 // convert the local Jacobian matrix and the right hand side from the data
194 // structures used by the automatic differentiation code to the conventional
195 // ones used by the linear solver.
196 updateLocalLinearization_(elemCtx, focusDofIdx);
197 }
198 }
199
203 LocalResidual& localResidual()
204 { return localResidual_; }
205
209 const LocalResidual& localResidual() const
210 { return localResidual_; }
211
220 const ScalarMatrixBlock& jacobian(unsigned domainScvIdx, unsigned rangeScvIdx) const
221 { return jacobian_[domainScvIdx][rangeScvIdx]; }
222
228 const ScalarVectorBlock& residual(unsigned dofIdx) const
229 { return residual_[dofIdx]; }
230
231protected:
232 Implementation& asImp_()
233 { return *static_cast<Implementation*>(this); }
234 const Implementation& asImp_() const
235 { return *static_cast<const Implementation*>(this); }
236
237 const Simulator& simulator_() const
238 { return *simulatorPtr_; }
239 const Problem& problem_() const
240 { return simulatorPtr_->problem(); }
241 const Model& model_() const
242 { return simulatorPtr_->model(); }
243
248 void resize_(const ElementContext& elemCtx)
249 {
250 size_t numDof = elemCtx.numDof(/*timeIdx=*/0);
251 size_t numPrimaryDof = elemCtx.numPrimaryDof(/*timeIdx=*/0);
252
253 residual_.resize(numDof);
254 if (jacobian_.N() != numDof || jacobian_.M() != numPrimaryDof)
255 jacobian_.setSize(numDof, numPrimaryDof);
256 }
257
261 void reset_(const ElementContext&)
262 {
263 residual_ = 0.0;
264 jacobian_ = 0.0;
265 }
266
271 void updateLocalLinearization_(const ElementContext& elemCtx,
272 unsigned focusDofIdx)
273 {
274 const auto& resid = localResidual_.residual();
275
276 for (unsigned eqIdx = 0; eqIdx < numEq; eqIdx++)
277 residual_[focusDofIdx][eqIdx] = resid[focusDofIdx][eqIdx].value();
278
279 size_t numDof = elemCtx.numDof(/*timeIdx=*/0);
280 for (unsigned dofIdx = 0; dofIdx < numDof; dofIdx++) {
281 for (unsigned eqIdx = 0; eqIdx < numEq; eqIdx++) {
282 for (unsigned pvIdx = 0; pvIdx < numEq; pvIdx++) {
283 // A[dofIdx][focusDofIdx][eqIdx][pvIdx] is the partial derivative of
284 // the residual function 'eqIdx' for the degree of freedom 'dofIdx'
285 // with regard to the focus variable 'pvIdx' of the degree of freedom
286 // 'focusDofIdx'
287 jacobian_[dofIdx][focusDofIdx][eqIdx][pvIdx] = resid[dofIdx][eqIdx].derivative(pvIdx);
288 Valgrind::CheckDefined(jacobian_[dofIdx][focusDofIdx][eqIdx][pvIdx]);
289 }
290 }
291 }
292 }
293
294 Simulator *simulatorPtr_;
295 Model *modelPtr_;
296
297 ElementContext *internalElemContext_;
298
299 LocalResidual localResidual_;
300
301 ScalarLocalBlockVector residual_;
302 ScalarLocalBlockMatrix jacobian_;
303};
304
305} // namespace Opm
306
307#endif
Calculates the local residual and its Jacobian for a single element of the grid.
Definition: fvbaseadlocallinearizer.hh:88
Simulator * simulatorPtr_
Definition: fvbaseadlocallinearizer.hh:294
const ScalarVectorBlock & residual(unsigned dofIdx) const
Returns the local residual of a sub-control volume.
Definition: fvbaseadlocallinearizer.hh:228
~FvBaseAdLocalLinearizer()
Definition: fvbaseadlocallinearizer.hh:119
void linearize(const Element &element)
Compute an element's local Jacobian matrix and evaluate its residual.
Definition: fvbaseadlocallinearizer.hh:154
Implementation & asImp_()
Definition: fvbaseadlocallinearizer.hh:232
const Problem & problem_() const
Definition: fvbaseadlocallinearizer.hh:239
ScalarLocalBlockVector residual_
Definition: fvbaseadlocallinearizer.hh:301
LocalResidual & localResidual()
Return reference to the local residual.
Definition: fvbaseadlocallinearizer.hh:203
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:271
void reset_(const ElementContext &)
Reset the all relevant internal attributes to 0.
Definition: fvbaseadlocallinearizer.hh:261
Model * modelPtr_
Definition: fvbaseadlocallinearizer.hh:295
const Implementation & asImp_() const
Definition: fvbaseadlocallinearizer.hh:234
FvBaseAdLocalLinearizer(const FvBaseAdLocalLinearizer &)=delete
ScalarLocalBlockMatrix jacobian_
Definition: fvbaseadlocallinearizer.hh:302
void resize_(const ElementContext &elemCtx)
Resize all internal attributes to the size of the element.
Definition: fvbaseadlocallinearizer.hh:248
const Simulator & simulator_() const
Definition: fvbaseadlocallinearizer.hh:237
const Model & model_() const
Definition: fvbaseadlocallinearizer.hh:241
const LocalResidual & localResidual() const
Return reference to the local residual.
Definition: fvbaseadlocallinearizer.hh:209
void init(Simulator &simulator)
Initialize the local Jacobian object.
Definition: fvbaseadlocallinearizer.hh:136
static void registerParameters()
Register all run-time parameters for the local jacobian.
Definition: fvbaseadlocallinearizer.hh:125
FvBaseAdLocalLinearizer()
Definition: fvbaseadlocallinearizer.hh:111
void linearize(ElementContext &elemCtx, const Element &elem)
Compute an element's local Jacobian matrix and evaluate its residual.
Definition: fvbaseadlocallinearizer.hh:173
const ScalarMatrixBlock & jacobian(unsigned domainScvIdx, unsigned rangeScvIdx) const
Returns the local Jacobian matrix of the residual of a sub-control volume.
Definition: fvbaseadlocallinearizer.hh:220
LocalResidual localResidual_
Definition: fvbaseadlocallinearizer.hh:299
ElementContext * internalElemContext_
Definition: fvbaseadlocallinearizer.hh:297
Declare the properties used by the infrastructure code of the finite volume discretizations.
Definition: blackoilmodel.hh:72
Definition: blackoilboundaryratevector.hh:37
typename Properties::Detail::GetPropImpl< TypeTag, Property >::type::type GetPropType
get the type alias defined in the property (equivalent to old macro GET_PROP_TYPE(....
Definition: propertysystem.hh:235
DenseAd::Evaluation< Scalar, numEq > type
Definition: fvbaseadlocallinearizer.hh:71
Representation of a function evaluation and all necessary derivatives with regard to the intensive qu...
Definition: fvbaseproperties.hh:66
The type of the local linearizer.
Definition: fvbaseproperties.hh:97
Definition: fvbaseadlocallinearizer.hh:53