vtkmultiphasemodule.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 */
27 #ifndef EWOMS_VTK_MULTI_PHASE_MODULE_HH
28 #define EWOMS_VTK_MULTI_PHASE_MODULE_HH
29 
30 #include "vtkmultiwriter.hh"
31 #include "baseoutputmodule.hh"
32 
35 
36 #include <opm/material/common/MathToolbox.hpp>
37 #include <opm/common/Valgrind.hpp>
38 
39 #include <dune/common/fvector.hh>
40 
41 #include <cstdio>
42 
43 namespace Ewoms {
44 namespace Properties {
45 // create new type tag for the VTK multi-phase output
46 NEW_TYPE_TAG(VtkMultiPhase);
47 
48 // create the property tags needed for the multi phase module
49 NEW_PROP_TAG(VtkWriteExtrusionFactor);
50 NEW_PROP_TAG(VtkWritePressures);
51 NEW_PROP_TAG(VtkWriteDensities);
52 NEW_PROP_TAG(VtkWriteSaturations);
53 NEW_PROP_TAG(VtkWriteMobilities);
54 NEW_PROP_TAG(VtkWriteRelativePermeabilities);
55 NEW_PROP_TAG(VtkWriteViscosities);
56 NEW_PROP_TAG(VtkWriteAverageMolarMasses);
57 NEW_PROP_TAG(VtkWritePorosity);
58 NEW_PROP_TAG(VtkWriteIntrinsicPermeabilities);
59 NEW_PROP_TAG(VtkWritePotentialGradients);
60 NEW_PROP_TAG(VtkWriteFilterVelocities);
61 NEW_PROP_TAG(VtkOutputFormat);
62 NEW_PROP_TAG(EnableVtkOutput);
63 NEW_PROP_TAG(Evaluation);
64 
65 // set default values for what quantities to output
66 SET_BOOL_PROP(VtkMultiPhase, VtkWriteExtrusionFactor, false);
67 SET_BOOL_PROP(VtkMultiPhase, VtkWritePressures, true);
68 SET_BOOL_PROP(VtkMultiPhase, VtkWriteDensities, true);
69 SET_BOOL_PROP(VtkMultiPhase, VtkWriteSaturations, true);
70 SET_BOOL_PROP(VtkMultiPhase, VtkWriteMobilities, false);
71 SET_BOOL_PROP(VtkMultiPhase, VtkWriteRelativePermeabilities, true);
72 SET_BOOL_PROP(VtkMultiPhase, VtkWriteViscosities, false);
73 SET_BOOL_PROP(VtkMultiPhase, VtkWriteAverageMolarMasses, false);
74 SET_BOOL_PROP(VtkMultiPhase, VtkWritePorosity, true);
75 SET_BOOL_PROP(VtkMultiPhase, VtkWriteIntrinsicPermeabilities, false);
76 SET_BOOL_PROP(VtkMultiPhase, VtkWritePotentialGradients, false);
77 SET_BOOL_PROP(VtkMultiPhase, VtkWriteFilterVelocities, false);
78 } // namespace Properties
79 
98 template<class TypeTag>
99 class VtkMultiPhaseModule : public BaseOutputModule<TypeTag>
100 {
102 
103  typedef typename GET_PROP_TYPE(TypeTag, Simulator) Simulator;
104  typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
105  typedef typename GET_PROP_TYPE(TypeTag, Evaluation) Evaluation;
106  typedef typename GET_PROP_TYPE(TypeTag, ElementContext) ElementContext;
107 
108  typedef typename GET_PROP_TYPE(TypeTag, GridView) GridView;
109  typedef typename GET_PROP_TYPE(TypeTag, FluidSystem) FluidSystem;
110  typedef typename GET_PROP_TYPE(TypeTag, DiscBaseOutputModule) DiscBaseOutputModule;
111 
112  static const int vtkFormat = GET_PROP_VALUE(TypeTag, VtkOutputFormat);
114 
115  typedef Opm::MathToolbox<Evaluation> Toolbox;
116 
117  enum { dimWorld = GridView::dimensionworld };
118  enum { numPhases = GET_PROP_VALUE(TypeTag, NumPhases) };
119 
120  typedef typename ParentType::ScalarBuffer ScalarBuffer;
121  typedef typename ParentType::VectorBuffer VectorBuffer;
122  typedef typename ParentType::TensorBuffer TensorBuffer;
123  typedef typename ParentType::PhaseBuffer PhaseBuffer;
124 
125  typedef Dune::FieldVector<Scalar, dimWorld> DimVector;
126 
127  typedef std::array<VectorBuffer, numPhases> PhaseVectorBuffer;
128 
129 public:
130  VtkMultiPhaseModule(const Simulator& simulator)
131  : ParentType(simulator)
132  {}
133 
137  static void registerParameters()
138  {
139  EWOMS_REGISTER_PARAM(TypeTag, bool, VtkWriteExtrusionFactor,
140  "Include the extrusion factor of the degrees of freedom into the VTK output files");
141  EWOMS_REGISTER_PARAM(TypeTag, bool, VtkWritePressures,
142  "Include the phase pressures in the VTK output files");
143  EWOMS_REGISTER_PARAM(TypeTag, bool, VtkWriteDensities,
144  "Include the phase densities in the VTK output files");
145  EWOMS_REGISTER_PARAM(TypeTag, bool, VtkWriteSaturations,
146  "Include the phase saturations in the VTK output files");
147  EWOMS_REGISTER_PARAM(TypeTag, bool, VtkWriteMobilities,
148  "Include the phase mobilities in the VTK output files");
149  EWOMS_REGISTER_PARAM(TypeTag, bool, VtkWriteRelativePermeabilities,
150  "Include the phase relative permeabilities in the VTK output files");
151  EWOMS_REGISTER_PARAM(TypeTag, bool, VtkWriteViscosities,
152  "Include component phase viscosities in the VTK output files");
153  EWOMS_REGISTER_PARAM(TypeTag, bool, VtkWriteAverageMolarMasses,
154  "Include the average phase mass in the VTK output files");
155  EWOMS_REGISTER_PARAM(TypeTag, bool, VtkWritePorosity,
156  "Include the porosity in the VTK output files");
157  EWOMS_REGISTER_PARAM(TypeTag, bool, VtkWriteIntrinsicPermeabilities,
158  "Include the intrinsic permeability in the VTK output files");
159  EWOMS_REGISTER_PARAM(TypeTag, bool, VtkWriteFilterVelocities,
160  "Include in the filter velocities of the phases the VTK output files");
161  EWOMS_REGISTER_PARAM(TypeTag, bool, VtkWritePotentialGradients,
162  "Include the phase pressure potential gradients in the VTK output files");
163  }
164 
170  {
171  if (extrusionFactorOutput_()) this->resizeScalarBuffer_(extrusionFactor_);
172  if (pressureOutput_()) this->resizePhaseBuffer_(pressure_);
173  if (densityOutput_()) this->resizePhaseBuffer_(density_);
174  if (saturationOutput_()) this->resizePhaseBuffer_(saturation_);
175  if (mobilityOutput_()) this->resizePhaseBuffer_(mobility_);
176  if (relativePermeabilityOutput_()) this->resizePhaseBuffer_(relativePermeability_);
177  if (viscosityOutput_()) this->resizePhaseBuffer_(viscosity_);
178  if (averageMolarMassOutput_()) this->resizePhaseBuffer_(averageMolarMass_);
179 
180  if (porosityOutput_()) this->resizeScalarBuffer_(porosity_);
181  if (intrinsicPermeabilityOutput_()) this->resizeTensorBuffer_(intrinsicPermeability_);
182 
183  if (velocityOutput_()) {
184  size_t nDof = this->simulator_.model().numGridDof();
185  for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++ phaseIdx) {
186  velocity_[phaseIdx].resize(nDof);
187  for (unsigned dofIdx = 0; dofIdx < nDof; ++ dofIdx) {
188  velocity_[phaseIdx][dofIdx].resize(dimWorld);
189  velocity_[phaseIdx][dofIdx] = 0.0;
190  }
191  }
192  this->resizePhaseBuffer_(velocityWeight_);
193  }
194 
195  if (potentialGradientOutput_()) {
196  size_t nDof = this->simulator_.model().numGridDof();
197  for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++ phaseIdx) {
198  potentialGradient_[phaseIdx].resize(nDof);
199  for (unsigned dofIdx = 0; dofIdx < nDof; ++ dofIdx) {
200  potentialGradient_[phaseIdx][dofIdx].resize(dimWorld);
201  potentialGradient_[phaseIdx][dofIdx] = 0.0;
202  }
203  }
204 
205  this->resizePhaseBuffer_(potentialWeight_);
206  }
207  }
208 
213  void processElement(const ElementContext& elemCtx)
214  {
215  typedef Opm::MathToolbox<Evaluation> Toolbox;
216 
217  if (!EWOMS_GET_PARAM(TypeTag, bool, EnableVtkOutput))
218  return;
219 
220  const auto& problem = elemCtx.problem();
221  for (unsigned i = 0; i < elemCtx.numPrimaryDof(/*timeIdx=*/0); ++i) {
222  unsigned I = elemCtx.globalSpaceIndex(i, /*timeIdx=*/0);
223  const auto& intQuants = elemCtx.intensiveQuantities(i, /*timeIdx=*/0);
224  const auto& fs = intQuants.fluidState();
225 
226  if (extrusionFactorOutput_()) extrusionFactor_[I] = intQuants.extrusionFactor();
227  if (porosityOutput_()) porosity_[I] = Toolbox::value(intQuants.porosity());
228  if (intrinsicPermeabilityOutput_()) {
229  const auto& K = problem.intrinsicPermeability(elemCtx, i, /*timeIdx=*/0);
230  intrinsicPermeability_[I].resize(K.rows, K.cols);
231  for (unsigned rowIdx = 0; rowIdx < K.rows; ++rowIdx)
232  for (unsigned colIdx = 0; colIdx < K.cols; ++colIdx)
233  intrinsicPermeability_[I][rowIdx][colIdx] = K[rowIdx][colIdx];
234  }
235 
236  for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
237  if (pressureOutput_())
238  pressure_[phaseIdx][I] = Toolbox::value(fs.pressure(phaseIdx));
239  if (densityOutput_())
240  density_[phaseIdx][I] = Toolbox::value(fs.density(phaseIdx));
241  if (saturationOutput_())
242  saturation_[phaseIdx][I] = Toolbox::value(fs.saturation(phaseIdx));
243  if (mobilityOutput_())
244  mobility_[phaseIdx][I] = Toolbox::value(intQuants.mobility(phaseIdx));
245  if (relativePermeabilityOutput_())
246  relativePermeability_[phaseIdx][I] = Toolbox::value(intQuants.relativePermeability(phaseIdx));
247  if (viscosityOutput_())
248  viscosity_[phaseIdx][I] = Toolbox::value(fs.viscosity(phaseIdx));
249  if (averageMolarMassOutput_())
250  averageMolarMass_[phaseIdx][I] = Toolbox::value(fs.averageMolarMass(phaseIdx));
251  }
252  }
253 
254  if (potentialGradientOutput_()) {
255  // calculate velocities if requested
256  for (unsigned faceIdx = 0; faceIdx < elemCtx.numInteriorFaces(/*timeIdx=*/0); ++ faceIdx) {
257  const auto& extQuants = elemCtx.extensiveQuantities(faceIdx, /*timeIdx=*/0);
258 
259  unsigned i = extQuants.interiorIndex();
260  unsigned I = elemCtx.globalSpaceIndex(i, /*timeIdx=*/0);
261 
262  for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
263  Scalar weight = extQuants.extrusionFactor();
264 
265  potentialWeight_[phaseIdx][I] += weight;
266 
267  const auto& inputPGrad = extQuants.potentialGrad(phaseIdx);
268  DimVector pGrad;
269  for (unsigned dimIdx = 0; dimIdx < dimWorld; ++dimIdx)
270  pGrad[dimIdx] = Toolbox::value(inputPGrad[dimIdx])*weight;
271  potentialGradient_[phaseIdx][I] += pGrad;
272  } // end for all phases
273  } // end for all faces
274  }
275 
276  if (velocityOutput_()) {
277  // calculate velocities if requested
278  for (unsigned faceIdx = 0; faceIdx < elemCtx.numInteriorFaces(/*timeIdx=*/0); ++ faceIdx) {
279  const auto& extQuants = elemCtx.extensiveQuantities(faceIdx, /*timeIdx=*/0);
280 
281  unsigned i = extQuants.interiorIndex();
282  unsigned I = elemCtx.globalSpaceIndex(i, /*timeIdx=*/0);
283 
284  unsigned j = extQuants.exteriorIndex();
285  unsigned J = elemCtx.globalSpaceIndex(j, /*timeIdx=*/0);
286 
287  for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
288  Scalar weight = std::max<Scalar>(1e-16,
289  std::abs(Toolbox::value(extQuants.volumeFlux(phaseIdx))));
290  Opm::Valgrind::CheckDefined(extQuants.extrusionFactor());
291  assert(extQuants.extrusionFactor() > 0);
292  weight *= extQuants.extrusionFactor();
293 
294  const auto& inputV = extQuants.filterVelocity(phaseIdx);
295  DimVector v;
296  for (unsigned k = 0; k < dimWorld; ++k)
297  v[k] = Toolbox::value(inputV[k]);
298  if (v.two_norm() > 1e-20)
299  weight /= v.two_norm();
300  v *= weight;
301 
302  velocity_[phaseIdx][I] += v;
303  velocity_[phaseIdx][J] += v;
304 
305  velocityWeight_[phaseIdx][I] += weight;
306  velocityWeight_[phaseIdx][J] += weight;
307  } // end for all phases
308  } // end for all faces
309  }
310  }
311 
315  void commitBuffers(BaseOutputWriter& baseWriter)
316  {
317  VtkMultiWriter *vtkWriter = dynamic_cast<VtkMultiWriter*>(&baseWriter);
318  if (!vtkWriter)
319  return;
320 
321  if (extrusionFactorOutput_())
322  this->commitScalarBuffer_(baseWriter, "extrusionFactor", extrusionFactor_);
323  if (pressureOutput_())
324  this->commitPhaseBuffer_(baseWriter, "pressure_%s", pressure_);
325  if (densityOutput_())
326  this->commitPhaseBuffer_(baseWriter, "density_%s", density_);
327  if (saturationOutput_())
328  this->commitPhaseBuffer_(baseWriter, "saturation_%s", saturation_);
329  if (mobilityOutput_())
330  this->commitPhaseBuffer_(baseWriter, "mobility_%s", mobility_);
331  if (relativePermeabilityOutput_())
332  this->commitPhaseBuffer_(baseWriter, "relativePerm_%s", relativePermeability_);
333  if (viscosityOutput_())
334  this->commitPhaseBuffer_(baseWriter, "viscosity_%s", viscosity_);
335  if (averageMolarMassOutput_())
336  this->commitPhaseBuffer_(baseWriter, "averageMolarMass_%s", averageMolarMass_);
337 
338  if (porosityOutput_())
339  this->commitScalarBuffer_(baseWriter, "porosity", porosity_);
340  if (intrinsicPermeabilityOutput_())
341  this->commitTensorBuffer_(baseWriter, "intrinsicPerm", intrinsicPermeability_);
342 
343  if (velocityOutput_()) {
344  size_t numDof = this->simulator_.model().numGridDof();
345 
346  for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
347  // first, divide the velocity field by the
348  // respective finite volume's surface area
349  for (unsigned i = 0; i < numDof; ++i)
350  velocity_[phaseIdx][i] /= velocityWeight_[phaseIdx][i];
351  // commit the phase velocity
352  char name[512];
353  snprintf(name, 512, "filterVelocity_%s", FluidSystem::phaseName(phaseIdx));
354 
355  DiscBaseOutputModule::attachVectorDofData_(baseWriter, velocity_[phaseIdx], name);
356  }
357  }
358 
359  if (potentialGradientOutput_()) {
360  size_t numDof = this->simulator_.model().numGridDof();
361 
362  for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
363  // first, divide the velocity field by the
364  // respective finite volume's surface area
365  for (unsigned i = 0; i < numDof; ++i)
366  potentialGradient_[phaseIdx][i] /= potentialWeight_[phaseIdx][i];
367  // commit the phase velocity
368  char name[512];
369  snprintf(name, 512, "gradP_%s", FluidSystem::phaseName(phaseIdx));
370 
371  DiscBaseOutputModule::attachVectorDofData_(baseWriter,
372  potentialGradient_[phaseIdx],
373  name);
374  }
375  }
376  }
377 
386  virtual bool needExtensiveQuantities() const final
387  {
388  return velocityOutput_() || potentialGradientOutput_();
389  }
390 
391 private:
392  static bool extrusionFactorOutput_()
393  { return EWOMS_GET_PARAM(TypeTag, bool, VtkWriteExtrusionFactor); }
394 
395  static bool pressureOutput_()
396  { return EWOMS_GET_PARAM(TypeTag, bool, VtkWritePressures); }
397 
398  static bool densityOutput_()
399  { return EWOMS_GET_PARAM(TypeTag, bool, VtkWriteDensities); }
400 
401  static bool saturationOutput_()
402  { return EWOMS_GET_PARAM(TypeTag, bool, VtkWriteSaturations); }
403 
404  static bool mobilityOutput_()
405  { return EWOMS_GET_PARAM(TypeTag, bool, VtkWriteMobilities); }
406 
407  static bool relativePermeabilityOutput_()
408  { return EWOMS_GET_PARAM(TypeTag, bool, VtkWriteRelativePermeabilities); }
409 
410  static bool viscosityOutput_()
411  { return EWOMS_GET_PARAM(TypeTag, bool, VtkWriteViscosities); }
412 
413  static bool averageMolarMassOutput_()
414  { return EWOMS_GET_PARAM(TypeTag, bool, VtkWriteAverageMolarMasses); }
415 
416  static bool porosityOutput_()
417  { return EWOMS_GET_PARAM(TypeTag, bool, VtkWritePorosity); }
418 
419  static bool intrinsicPermeabilityOutput_()
420  { return EWOMS_GET_PARAM(TypeTag, bool, VtkWriteIntrinsicPermeabilities); }
421 
422  static bool velocityOutput_()
423  { return EWOMS_GET_PARAM(TypeTag, bool, VtkWriteFilterVelocities); }
424 
425  static bool potentialGradientOutput_()
426  { return EWOMS_GET_PARAM(TypeTag, bool, VtkWritePotentialGradients); }
427 
428  ScalarBuffer extrusionFactor_;
429  PhaseBuffer pressure_;
430  PhaseBuffer density_;
431  PhaseBuffer saturation_;
432  PhaseBuffer mobility_;
433  PhaseBuffer relativePermeability_;
434  PhaseBuffer viscosity_;
435  PhaseBuffer averageMolarMass_;
436 
437  ScalarBuffer porosity_;
438  TensorBuffer intrinsicPermeability_;
439 
440  PhaseVectorBuffer velocity_;
441  PhaseBuffer velocityWeight_;
442 
443  PhaseVectorBuffer potentialGradient_;
444  PhaseBuffer potentialWeight_;
445 };
446 
447 } // namespace Ewoms
448 
449 #endif
void processElement(const ElementContext &elemCtx)
Modify the internal buffers according to the intensive quantities seen on an element.
Definition: vtkmultiphasemodule.hh:213
The base class for all output writers.
Definition: baseoutputwriter.hh:43
#define SET_BOOL_PROP(EffTypeTagName, PropTagName,...)
Set a property to a simple constant boolean value.
Definition: propertysystem.hh:361
Definition: baseauxiliarymodule.hh:37
Model & model()
Return the physical model used in the simulation.
Definition: simulator.hh:189
void commitPhaseBuffer_(BaseOutputWriter &baseWriter, const char *pattern, PhaseBuffer &buffer, BufferType bufferType=DofBuffer)
Add a phase-specific buffer to the result file.
Definition: baseoutputmodule.hh:408
VTK output module for quantities which make sense for all models which deal with multiple fluid phase...
Definition: vtkmultiphasemodule.hh:99
Simplifies writing multi-file VTK datasets.
Definition: vtkmultiwriter.hh:63
void commitScalarBuffer_(BaseOutputWriter &baseWriter, const char *name, ScalarBuffer &buffer, BufferType bufferType=DofBuffer)
Add a buffer containing scalar quantities to the result file.
Definition: baseoutputmodule.hh:305
#define GET_PROP_VALUE(TypeTag, PropTagName)
Access the value attribute of a property for a type tag.
Definition: propertysystem.hh:469
The base class for writer modules.
void commitBuffers(BaseOutputWriter &baseWriter)
Add all buffers to the VTK output writer.
Definition: vtkmultiphasemodule.hh:315
#define NEW_TYPE_TAG(...)
Define a new type tag.
Definition: propertysystem.hh:169
void resizeTensorBuffer_(TensorBuffer &buffer, BufferType bufferType=DofBuffer)
Allocate the space for a buffer storing a tensorial quantity.
Definition: baseoutputmodule.hh:190
void commitTensorBuffer_(BaseOutputWriter &baseWriter, const char *name, TensorBuffer &buffer, BufferType bufferType=DofBuffer)
Add a buffer containing tensorial quantities to the result file.
Definition: baseoutputmodule.hh:341
#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.
The base class for writer modules.
Definition: baseoutputmodule.hh:80
void allocBuffers()
Allocate memory for the scalar fields we would like to write to the VTK file.
Definition: vtkmultiphasemodule.hh:169
#define EWOMS_GET_PARAM(TypeTag, ParamType, ParamName)
Retrieve a runtime parameter.
Definition: parametersystem.hh:99
virtual bool needExtensiveQuantities() const final
Returns true iff the module needs to access the extensive quantities of a context to do its job...
Definition: vtkmultiphasemodule.hh:386
void resizePhaseBuffer_(PhaseBuffer &buffer, BufferType bufferType=DofBuffer)
Allocate the space for a buffer storing a phase-specific quantity.
Definition: baseoutputmodule.hh:235
static void registerParameters()
Register all run-time parameters for the multi-phase VTK output module.
Definition: vtkmultiphasemodule.hh:137
Provides the magic behind the eWoms property system.
void resizeScalarBuffer_(ScalarBuffer &buffer, BufferType bufferType=DofBuffer)
Allocate the space for a buffer storing a scalar quantity.
Definition: baseoutputmodule.hh:170
Simplifies writing multi-file VTK datasets.
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