diff --git a/examples/problems/diffusionproblem.hh b/examples/problems/diffusionproblem.hh index 4105f287c..4717c79ee 100644 --- a/examples/problems/diffusionproblem.hh +++ b/examples/problems/diffusionproblem.hh @@ -30,7 +30,7 @@ #include -#include +#include #include #include diff --git a/examples/problems/fingerproblem.hh b/examples/problems/fingerproblem.hh index c6d405432..7e4db6bec 100644 --- a/examples/problems/fingerproblem.hh +++ b/examples/problems/fingerproblem.hh @@ -28,7 +28,7 @@ #ifndef EWOMS_FINGER_PROBLEM_HH #define EWOMS_FINGER_PROBLEM_HH -#include +#include #include #include diff --git a/examples/problems/fractureproblem.hh b/examples/problems/fractureproblem.hh index 8b18b0ad0..06c8c40b8 100644 --- a/examples/problems/fractureproblem.hh +++ b/examples/problems/fractureproblem.hh @@ -38,7 +38,7 @@ #endif #include -#include +#include #include #include diff --git a/examples/problems/lensproblem.hh b/examples/problems/lensproblem.hh index c0ca86bd3..094d042a4 100644 --- a/examples/problems/lensproblem.hh +++ b/examples/problems/lensproblem.hh @@ -28,7 +28,7 @@ #ifndef EWOMS_LENS_PROBLEM_HH #define EWOMS_LENS_PROBLEM_HH -#include +#include #include #include #include diff --git a/examples/problems/powerinjectionproblem.hh b/examples/problems/powerinjectionproblem.hh index b88c06eef..50fc96753 100644 --- a/examples/problems/powerinjectionproblem.hh +++ b/examples/problems/powerinjectionproblem.hh @@ -29,7 +29,7 @@ #define EWOMS_POWER_INJECTION_PROBLEM_HH #include -#include +#include #include #include diff --git a/examples/tutorial1problem.hh b/examples/tutorial1problem.hh index a9fa5a71e..1d01cb2e0 100644 --- a/examples/tutorial1problem.hh +++ b/examples/tutorial1problem.hh @@ -45,7 +45,7 @@ // For the DUNE grid #include /*@\label{tutorial1:include-grid-manager}@*/ -#include /*@\label{tutorial1:include-grid-manager}@*/ +#include /*@\label{tutorial1:include-grid-manager}@*/ // For Dune::FieldMatrix #include diff --git a/opm/models/blackoil/blackoilenergymodules.hh b/opm/models/blackoil/blackoilenergymodules.hh index e94ab7f48..6cae6c4f3 100644 --- a/opm/models/blackoil/blackoilenergymodules.hh +++ b/opm/models/blackoil/blackoilenergymodules.hh @@ -29,7 +29,7 @@ #define EWOMS_BLACK_OIL_ENERGY_MODULE_HH #include "blackoilproperties.hh" -#include +#include #include #include diff --git a/opm/models/blackoil/blackoilfoammodules.hh b/opm/models/blackoil/blackoilfoammodules.hh index dbdf7fdd6..df95f5347 100644 --- a/opm/models/blackoil/blackoilfoammodules.hh +++ b/opm/models/blackoil/blackoilfoammodules.hh @@ -29,7 +29,7 @@ #define EWOMS_BLACK_OIL_FOAM_MODULE_HH #include "blackoilproperties.hh" -//#include +//#include #include #include diff --git a/opm/models/blackoil/blackoilmodel.hh b/opm/models/blackoil/blackoilmodel.hh index 15dc46647..ec5d8ecca 100644 --- a/opm/models/blackoil/blackoilmodel.hh +++ b/opm/models/blackoil/blackoilmodel.hh @@ -47,8 +47,8 @@ #include "blackoildarcyfluxmodule.hh" #include -#include -#include +#include +#include #include #include diff --git a/opm/models/blackoil/blackoilpolymermodules.hh b/opm/models/blackoil/blackoilpolymermodules.hh index 49e8c204f..6c3f7e9d4 100644 --- a/opm/models/blackoil/blackoilpolymermodules.hh +++ b/opm/models/blackoil/blackoilpolymermodules.hh @@ -29,7 +29,7 @@ #define EWOMS_BLACK_OIL_POLYMER_MODULE_HH #include "blackoilproperties.hh" -#include +#include #include #include diff --git a/opm/models/blackoil/blackoilsolventmodules.hh b/opm/models/blackoil/blackoilsolventmodules.hh index a03632ced..573ed6b9d 100644 --- a/opm/models/blackoil/blackoilsolventmodules.hh +++ b/opm/models/blackoil/blackoilsolventmodules.hh @@ -29,7 +29,7 @@ #define EWOMS_BLACK_OIL_SOLVENT_MODULE_HH #include "blackoilproperties.hh" -#include +#include #include #include diff --git a/opm/models/common/multiphasebaseproperties.hh b/opm/models/common/multiphasebaseproperties.hh index 7697fabee..ec391e49e 100644 --- a/opm/models/common/multiphasebaseproperties.hh +++ b/opm/models/common/multiphasebaseproperties.hh @@ -31,8 +31,8 @@ #define EWOMS_MULTI_PHASE_BASE_PROPERTIES_HH #include -#include -#include +#include +#include BEGIN_PROPERTIES diff --git a/opm/models/discretization/common/fvbasenewtonconvergencewriter.hh b/opm/models/discretization/common/fvbasenewtonconvergencewriter.hh index 16557fd8c..b66b44bbd 100644 --- a/opm/models/discretization/common/fvbasenewtonconvergencewriter.hh +++ b/opm/models/discretization/common/fvbasenewtonconvergencewriter.hh @@ -28,7 +28,7 @@ #ifndef EWOMS_FV_BASE_NEWTON_CONVERGENCE_WRITER_HH #define EWOMS_FV_BASE_NEWTON_CONVERGENCE_WRITER_HH -#include +#include #include #include diff --git a/opm/models/discretization/common/fvbaseproblem.hh b/opm/models/discretization/common/fvbaseproblem.hh index edfdf8c76..be9d98c3d 100644 --- a/opm/models/discretization/common/fvbaseproblem.hh +++ b/opm/models/discretization/common/fvbaseproblem.hh @@ -30,8 +30,8 @@ #include "fvbaseproperties.hh" -#include -#include +#include +#include #include #include diff --git a/opm/models/discretization/common/fvbaseproperties.hh b/opm/models/discretization/common/fvbaseproperties.hh index 67248ffdb..bf643f100 100644 --- a/opm/models/discretization/common/fvbaseproperties.hh +++ b/opm/models/discretization/common/fvbaseproperties.hh @@ -35,7 +35,7 @@ #include "fvbasefdlocallinearizer.hh" #include -#include +#include #include BEGIN_PROPERTIES diff --git a/opm/models/discretization/ecfv/ecfvbaseoutputmodule.hh b/opm/models/discretization/ecfv/ecfvbaseoutputmodule.hh index 234016996..0db7c3e59 100644 --- a/opm/models/discretization/ecfv/ecfvbaseoutputmodule.hh +++ b/opm/models/discretization/ecfv/ecfvbaseoutputmodule.hh @@ -30,7 +30,7 @@ #include "ecfvproperties.hh" -#include +#include namespace Opm { /*! diff --git a/opm/models/discretization/vcfv/vcfvbaseoutputmodule.hh b/opm/models/discretization/vcfv/vcfvbaseoutputmodule.hh index bccebe1f8..bda48d9ff 100644 --- a/opm/models/discretization/vcfv/vcfvbaseoutputmodule.hh +++ b/opm/models/discretization/vcfv/vcfvbaseoutputmodule.hh @@ -30,7 +30,7 @@ #include "vcfvproperties.hh" -#include +#include #include #include diff --git a/opm/models/flash/flashmodel.hh b/opm/models/flash/flashmodel.hh index 14fa5bda1..de04ae438 100644 --- a/opm/models/flash/flashmodel.hh +++ b/opm/models/flash/flashmodel.hh @@ -41,9 +41,9 @@ #include #include -#include -#include -#include +#include +#include +#include #include #include #include diff --git a/opm/models/flash/flashproperties.hh b/opm/models/flash/flashproperties.hh index a679a3f3c..108b6d63b 100644 --- a/opm/models/flash/flashproperties.hh +++ b/opm/models/flash/flashproperties.hh @@ -31,9 +31,9 @@ #define EWOMS_FLASH_PROPERTIES_HH #include -#include -#include -#include +#include +#include +#include BEGIN_PROPERTIES diff --git a/opm/models/immiscible/immisciblemodel.hh b/opm/models/immiscible/immisciblemodel.hh index 4193f7ce9..278bbe996 100644 --- a/opm/models/immiscible/immisciblemodel.hh +++ b/opm/models/immiscible/immisciblemodel.hh @@ -40,7 +40,7 @@ #include #include -#include +#include #include #include #include diff --git a/opm/models/immiscible/immiscibleproperties.hh b/opm/models/immiscible/immiscibleproperties.hh index 273d1ca0e..77380a375 100644 --- a/opm/models/immiscible/immiscibleproperties.hh +++ b/opm/models/immiscible/immiscibleproperties.hh @@ -31,7 +31,7 @@ #define EWOMS_IMMISCIBLE_PROPERTIES_HH #include -#include +#include BEGIN_PROPERTIES diff --git a/opm/models/io/baseoutputmodule.hh b/opm/models/io/baseoutputmodule.hh new file mode 100644 index 000000000..fc18a5ba9 --- /dev/null +++ b/opm/models/io/baseoutputmodule.hh @@ -0,0 +1,518 @@ +// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- +// vi: set et ts=4 sw=4 sts=4: +/* + This file is part of the Open Porous Media project (OPM). + + OPM is free software: you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation, either version 2 of the License, or + (at your option) any later version. + + OPM is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU General Public License + along with OPM. If not, see . + + Consult the COPYING file in the top-level source directory of this + module for the precise wording of the license and the list of + copyright holders. +*/ +/*! + * \file + * \copydoc Opm::BaseOutputModule + */ +#ifndef EWOMS_BASE_OUTPUT_MODULE_HH +#define EWOMS_BASE_OUTPUT_MODULE_HH + +#include "baseoutputwriter.hh" + +#include +#include + +#include + +#include +#include + +#include +#include +#include +#include + +#include + +BEGIN_PROPERTIES + +// forward definition of property tags +NEW_PROP_TAG(NumPhases); +NEW_PROP_TAG(NumComponents); +NEW_PROP_TAG(NumEq); + +NEW_PROP_TAG(Model); +NEW_PROP_TAG(Simulator); +NEW_PROP_TAG(Scalar); +NEW_PROP_TAG(Evaluation); +NEW_PROP_TAG(GridView); +NEW_PROP_TAG(ElementContext); +NEW_PROP_TAG(FluidSystem); +NEW_PROP_TAG(DiscBaseOutputModule); + +END_PROPERTIES + +namespace Opm { + +#if __GNUC__ || __clang__ +#pragma GCC diagnostic push +#pragma GCC diagnostic ignored "-Wpragmas" +#pragma GCC diagnostic ignored "-Wformat-nonliteral" +#endif + +/*! + * \brief The base class for writer modules. + * + * This class also provides some convenience methods for buffer + * management and is the base class for all other output writer + * modules. + */ +template +class BaseOutputModule +{ + typedef typename GET_PROP_TYPE(TypeTag, Simulator) Simulator; + typedef typename GET_PROP_TYPE(TypeTag, Model) Model; + typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar; + typedef typename GET_PROP_TYPE(TypeTag, GridView) GridView; + typedef typename GET_PROP_TYPE(TypeTag, ElementContext) ElementContext; + typedef typename GET_PROP_TYPE(TypeTag, FluidSystem) FluidSystem; + typedef typename GET_PROP_TYPE(TypeTag, DiscBaseOutputModule) DiscBaseOutputModule; + + enum { numPhases = GET_PROP_VALUE(TypeTag, NumPhases) }; + enum { numComponents = GET_PROP_VALUE(TypeTag, NumComponents) }; + enum { numEq = GET_PROP_VALUE(TypeTag, NumEq) }; + enum { dim = GridView::dimension }; + enum { dimWorld = GridView::dimensionworld }; + + typedef BaseOutputWriter::Tensor Tensor; + +public: + typedef BaseOutputWriter::ScalarBuffer ScalarBuffer; + typedef BaseOutputWriter::VectorBuffer VectorBuffer; + typedef BaseOutputWriter::TensorBuffer TensorBuffer; + + typedef std::array EqBuffer; + typedef std::array PhaseBuffer; + typedef std::array ComponentBuffer; + typedef std::array, numPhases> PhaseComponentBuffer; + + typedef std::array PhaseVectorBuffer; + + BaseOutputModule(const Simulator& simulator) + : simulator_(simulator) + {} + + virtual ~BaseOutputModule() + {} + + /*! + * \brief Allocate memory for the scalar fields we would like to + * write to disk. + * + * The module can dynamically cast the writer to the desired + * concrete class. If the writer is incompatible with the module, + * this method should become a no-op. + */ + virtual void allocBuffers() = 0; + + /*! + * \brief Modify the internal buffers according to the intensive quanties relevant + * for an element + * + * The module can dynamically cast the writer to the desired + * concrete class. If the writer is incompatible with the module, + * this method should become a no-op. + */ + virtual void processElement(const ElementContext& elemCtx) = 0; + + /*! + * \brief Add all buffers to the VTK output writer. + */ + virtual void commitBuffers(BaseOutputWriter& writer) = 0; + + /*! + * \brief Returns true iff the module needs to access the extensive quantities of a + * context to do its job. + * + * For example, this happens if velocities or gradients should be written. + * + * Always returning true here does not do any harm from the correctness perspective, + * but it slows down writing the output fields. Since most output modules only write + * intensive quantities, this method returns 'false' by default. + */ + virtual bool needExtensiveQuantities() const + { return false; } + +protected: + enum BufferType { + //! Buffer contains data associated with the degrees of freedom + DofBuffer, + + //! Buffer contains data associated with the grid's vertices + VertexBuffer, + + //! Buffer contains data associated with the grid's elements + ElementBuffer + }; + + /*! + * \brief Allocate the space for a buffer storing a scalar quantity + */ + void resizeScalarBuffer_(ScalarBuffer& buffer, + BufferType bufferType = DofBuffer) + { + size_t n; + if (bufferType == VertexBuffer) + n = static_cast(simulator_.gridView().size(dim)); + else if (bufferType == ElementBuffer) + n = static_cast(simulator_.gridView().size(0)); + else if (bufferType == DofBuffer) + n = simulator_.model().numGridDof(); + else + throw std::logic_error("bufferType must be one of Dof, Vertex or Element"); + + buffer.resize(n); + std::fill(buffer.begin(), buffer.end(), 0.0); + } + + /*! + * \brief Allocate the space for a buffer storing a tensorial quantity + */ + void resizeTensorBuffer_(TensorBuffer& buffer, + BufferType bufferType = DofBuffer) + { + size_t n; + if (bufferType == VertexBuffer) + n = static_cast(simulator_.gridView().size(dim)); + else if (bufferType == ElementBuffer) + n = static_cast(simulator_.gridView().size(0)); + else if (bufferType == DofBuffer) + n = simulator_.model().numGridDof(); + else + throw std::logic_error("bufferType must be one of Dof, Vertex or Element"); + + buffer.resize(n); + Tensor nullMatrix(dimWorld, dimWorld, 0.0); + std::fill(buffer.begin(), buffer.end(), nullMatrix); + } + + /*! + * \brief Allocate the space for a buffer storing a equation specific + * quantity + */ + void resizeEqBuffer_(EqBuffer& buffer, + BufferType bufferType = DofBuffer) + { + size_t n; + if (bufferType == VertexBuffer) + n = static_cast(simulator_.gridView().size(dim)); + else if (bufferType == ElementBuffer) + n = static_cast(simulator_.gridView().size(0)); + else if (bufferType == DofBuffer) + n = simulator_.model().numGridDof(); + else + throw std::logic_error("bufferType must be one of Dof, Vertex or Element"); + + for (unsigned i = 0; i < numEq; ++i) { + buffer[i].resize(n); + std::fill(buffer[i].begin(), buffer[i].end(), 0.0); + } + } + + /*! + * \brief Allocate the space for a buffer storing a phase-specific + * quantity + */ + void resizePhaseBuffer_(PhaseBuffer& buffer, + BufferType bufferType = DofBuffer) + { + size_t n; + if (bufferType == VertexBuffer) + n = static_cast(simulator_.gridView().size(dim)); + else if (bufferType == ElementBuffer) + n = static_cast(simulator_.gridView().size(0)); + else if (bufferType == DofBuffer) + n = simulator_.model().numGridDof(); + else + throw std::logic_error("bufferType must be one of Dof, Vertex or Element"); + + for (unsigned i = 0; i < numPhases; ++i) { + buffer[i].resize(n); + std::fill(buffer[i].begin(), buffer[i].end(), 0.0); + } + } + + /*! + * \brief Allocate the space for a buffer storing a component + * specific quantity + */ + void resizeComponentBuffer_(ComponentBuffer& buffer, + BufferType bufferType = DofBuffer) + { + size_t n; + if (bufferType == VertexBuffer) + n = static_cast(simulator_.gridView().size(dim)); + else if (bufferType == ElementBuffer) + n = static_cast(simulator_.gridView().size(0)); + else if (bufferType == DofBuffer) + n = simulator_.model().numGridDof(); + else + throw std::logic_error("bufferType must be one of Dof, Vertex or Element"); + + for (unsigned i = 0; i < numComponents; ++i) { + buffer[i].resize(n); + std::fill(buffer[i].begin(), buffer[i].end(), 0.0); + } + } + + /*! + * \brief Allocate the space for a buffer storing a phase and + * component specific buffer + */ + void resizePhaseComponentBuffer_(PhaseComponentBuffer& buffer, + BufferType bufferType = DofBuffer) + { + size_t n; + if (bufferType == VertexBuffer) + n = static_cast(simulator_.gridView().size(dim)); + else if (bufferType == ElementBuffer) + n = static_cast(simulator_.gridView().size(0)); + else if (bufferType == DofBuffer) + n = simulator_.model().numGridDof(); + else + throw std::logic_error("bufferType must be one of Dof, Vertex or Element"); + + for (unsigned i = 0; i < numPhases; ++i) { + for (unsigned j = 0; j < numComponents; ++j) { + buffer[i][j].resize(n); + std::fill(buffer[i][j].begin(), buffer[i][j].end(), 0.0); + } + } + } + + /*! + * \brief Add a buffer containing scalar quantities to the result file. + */ + void commitScalarBuffer_(BaseOutputWriter& baseWriter, + const char *name, + ScalarBuffer& buffer, + BufferType bufferType = DofBuffer) + { + if (bufferType == DofBuffer) + DiscBaseOutputModule::attachScalarDofData_(baseWriter, buffer, name); + else if (bufferType == VertexBuffer) + attachScalarVertexData_(baseWriter, buffer, name); + else if (bufferType == ElementBuffer) + attachScalarElementData_(baseWriter, buffer, name); + else + throw std::logic_error("bufferType must be one of Dof, Vertex or Element"); + } + + /*! + * \brief Add a buffer containing vectorial quantities to the result file. + */ + void commitVectorBuffer_(BaseOutputWriter& baseWriter, + const char *name, + VectorBuffer& buffer, + BufferType bufferType = DofBuffer) + { + if (bufferType == DofBuffer) + DiscBaseOutputModule::attachVectorDofData_(baseWriter, buffer, name); + else if (bufferType == VertexBuffer) + attachVectorVertexData_(baseWriter, buffer, name); + else if (bufferType == ElementBuffer) + attachVectorElementData_(baseWriter, buffer, name); + else + throw std::logic_error("bufferType must be one of Dof, Vertex or Element"); + } + + /*! + * \brief Add a buffer containing tensorial quantities to the result file. + */ + void commitTensorBuffer_(BaseOutputWriter& baseWriter, + const char *name, + TensorBuffer& buffer, + BufferType bufferType = DofBuffer) + { + if (bufferType == DofBuffer) + DiscBaseOutputModule::attachTensorDofData_(baseWriter, buffer, name); + else if (bufferType == VertexBuffer) + attachTensorVertexData_(baseWriter, buffer, name); + else if (bufferType == ElementBuffer) + attachTensorElementData_(baseWriter, buffer, name); + else + throw std::logic_error("bufferType must be one of Dof, Vertex or Element"); + } + + /*! + * \brief Add a buffer with as many variables as PDEs to the result file. + */ + void commitPriVarsBuffer_(BaseOutputWriter& baseWriter, + const char *pattern, + EqBuffer& buffer, + BufferType bufferType = DofBuffer) + { + char name[512]; + for (unsigned i = 0; i < numEq; ++i) { + std::string eqName = simulator_.model().primaryVarName(i); + snprintf(name, 512, pattern, eqName.c_str()); + + if (bufferType == DofBuffer) + DiscBaseOutputModule::attachScalarDofData_(baseWriter, buffer[i], name); + else if (bufferType == VertexBuffer) + attachScalarVertexData_(baseWriter, buffer[i], name); + else if (bufferType == ElementBuffer) + attachScalarElementData_(baseWriter, buffer[i], name); + else + throw std::logic_error("bufferType must be one of Dof, Vertex or Element"); + } + } + + /*! + * \brief Add a buffer with as many variables as PDEs to the result file. + */ + void commitEqBuffer_(BaseOutputWriter& baseWriter, + const char *pattern, + EqBuffer& buffer, + BufferType bufferType = DofBuffer) + { + char name[512]; + for (unsigned i = 0; i < numEq; ++i) { + std::ostringstream oss; + oss << i; + snprintf(name, 512, pattern, oss.str().c_str()); + + if (bufferType == DofBuffer) + DiscBaseOutputModule::attachScalarDofData_(baseWriter, buffer[i], name); + else if (bufferType == VertexBuffer) + attachScalarVertexData_(baseWriter, buffer[i], name); + else if (bufferType == ElementBuffer) + attachScalarElementData_(baseWriter, buffer[i], name); + else + throw std::logic_error("bufferType must be one of Dof, Vertex or Element"); + } + } + + /*! + * \brief Add a phase-specific buffer to the result file. + */ + void commitPhaseBuffer_(BaseOutputWriter& baseWriter, + const char *pattern, + PhaseBuffer& buffer, + BufferType bufferType = DofBuffer) + { + char name[512]; + for (unsigned i = 0; i < numPhases; ++i) { + snprintf(name, 512, pattern, FluidSystem::phaseName(i)); + + if (bufferType == DofBuffer) + DiscBaseOutputModule::attachScalarDofData_(baseWriter, buffer[i], name); + else if (bufferType == VertexBuffer) + attachScalarVertexData_(baseWriter, buffer[i], name); + else if (bufferType == ElementBuffer) + attachScalarElementData_(baseWriter, buffer[i], name); + else + throw std::logic_error("bufferType must be one of Dof, Vertex or Element"); + } + } + + /*! + * \brief Add a component-specific buffer to the result file. + */ + void commitComponentBuffer_(BaseOutputWriter& baseWriter, + const char *pattern, + ComponentBuffer& buffer, + BufferType bufferType = DofBuffer) + { + char name[512]; + for (unsigned i = 0; i < numComponents; ++i) { + snprintf(name, 512, pattern, FluidSystem::componentName(i)); + + if (bufferType == DofBuffer) + DiscBaseOutputModule::attachScalarDofData_(baseWriter, buffer[i], name); + else if (bufferType == VertexBuffer) + attachScalarVertexData_(baseWriter, buffer[i], name); + else if (bufferType == ElementBuffer) + attachScalarElementData_(baseWriter, buffer[i], name); + else + throw std::logic_error("bufferType must be one of Dof, Vertex or Element"); + } + } + + /*! + * \brief Add a phase and component specific quantities to the output. + */ + void commitPhaseComponentBuffer_(BaseOutputWriter& baseWriter, + const char *pattern, + PhaseComponentBuffer& buffer, + BufferType bufferType = DofBuffer) + { + char name[512]; + for (unsigned i= 0; i < numPhases; ++i) { + for (unsigned j = 0; j < numComponents; ++j) { + snprintf(name, 512, pattern, + FluidSystem::phaseName(i), + FluidSystem::componentName(j)); + + if (bufferType == DofBuffer) + DiscBaseOutputModule::attachScalarDofData_(baseWriter, buffer[i][j], name); + else if (bufferType == VertexBuffer) + attachScalarVertexData_(baseWriter, buffer[i][j], name); + else if (bufferType == ElementBuffer) + attachScalarElementData_(baseWriter, buffer[i][j], name); + else + throw std::logic_error("bufferType must be one of Dof, Vertex or Element"); + } + } + } + + void attachScalarElementData_(BaseOutputWriter& baseWriter, + ScalarBuffer& buffer, + const char *name) + { baseWriter.attachScalarElementData(buffer, name); } + + void attachScalarVertexData_(BaseOutputWriter& baseWriter, + ScalarBuffer& buffer, + const char *name) + { baseWriter.attachScalarVertexData(buffer, name); } + + void attachVectorElementData_(BaseOutputWriter& baseWriter, + VectorBuffer& buffer, + const char *name) + { baseWriter.attachVectorElementData(buffer, name); } + + void attachVectorVertexData_(BaseOutputWriter& baseWriter, + VectorBuffer& buffer, + const char *name) + { baseWriter.attachVectorVertexData(buffer, name); } + + void attachTensorElementData_(BaseOutputWriter& baseWriter, + TensorBuffer& buffer, + const char *name) + { baseWriter.attachTensorElementData(buffer, name); } + + void attachTensorVertexData_(BaseOutputWriter& baseWriter, + TensorBuffer& buffer, + const char *name) + { baseWriter.attachTensorVertexData(buffer, name); } + + const Simulator& simulator_; +}; + +#if __GNUC__ || __clang__ +#pragma GCC diagnostic pop +#endif + +} // namespace Opm + +#endif diff --git a/opm/models/io/baseoutputwriter.hh b/opm/models/io/baseoutputwriter.hh new file mode 100644 index 000000000..075b69620 --- /dev/null +++ b/opm/models/io/baseoutputwriter.hh @@ -0,0 +1,106 @@ +// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- +// vi: set et ts=4 sw=4 sts=4: +/* + This file is part of the Open Porous Media project (OPM). + + OPM is free software: you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation, either version 2 of the License, or + (at your option) any later version. + + OPM is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU General Public License + along with OPM. If not, see . + + Consult the COPYING file in the top-level source directory of this + module for the precise wording of the license and the list of + copyright holders. +*/ +/*! + * \file + * + * \copydoc Opm::BaseOutputWriter + */ +#ifndef EWOMS_BASE_OUTPUT_WRITER_HH +#define EWOMS_BASE_OUTPUT_WRITER_HH + +#include +#include + +#include + +namespace Opm { +/*! + * \brief The base class for all output writers. + * + * The sole purpose of this class is to enable RTTI + * (i.e. dynamic_cast) on writer objects. + */ +class BaseOutputWriter +{ +public: + typedef double Scalar; + typedef Dune::DynamicVector Vector; + typedef Dune::DynamicMatrix Tensor; + typedef std::vector ScalarBuffer; + typedef std::vector VectorBuffer; + typedef std::vector TensorBuffer; + + BaseOutputWriter() + {} + + virtual ~BaseOutputWriter() + {} + + /*! + * \brief Called when ever a new time step or a new grid + * must be written. + */ + virtual void beginWrite(double t) = 0; + + /*! + * \brief Add a scalar vertex centered vector field to the output. + */ + virtual void attachScalarVertexData(ScalarBuffer& buf, std::string name) = 0; + + /*! + * \brief Add a scalar element centered quantity to the output. + */ + virtual void attachScalarElementData(ScalarBuffer& buf, std::string name) = 0; + + /*! + * \brief Add a vectorial vertex centered vector field to the output. + */ + virtual void attachVectorVertexData(VectorBuffer& buf, std::string name) = 0; + + /*! + * \brief Add a vectorial element centered quantity to the output. + */ + virtual void attachVectorElementData(VectorBuffer& buf, std::string name) = 0; + + /*! + * \brief Add a tensorial vertex centered tensor field to the output. + */ + virtual void attachTensorVertexData(TensorBuffer& buf, std::string name) = 0; + + /*! + * \brief Add a tensorial element centered quantity to the output. + */ + virtual void attachTensorElementData(TensorBuffer& buf, std::string name) = 0; + + /*! + * \brief Finalizes the current writer. + * + * This means that everything will be written to disk, except if + * the onlyDiscard argument is true. In this case only all managed + * buffers are deleted, but no output is written. + */ + virtual void endWrite(bool onlyDiscard = false) = 0; +}; +} // namespace Opm + +#endif diff --git a/opm/models/io/basevanguard.hh b/opm/models/io/basevanguard.hh new file mode 100644 index 000000000..061de0764 --- /dev/null +++ b/opm/models/io/basevanguard.hh @@ -0,0 +1,163 @@ +// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- +// vi: set et ts=4 sw=4 sts=4: +/* + This file is part of the Open Porous Media project (OPM). + + OPM is free software: you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation, either version 2 of the License, or + (at your option) any later version. + + OPM is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU General Public License + along with OPM. If not, see . + + Consult the COPYING file in the top-level source directory of this + module for the precise wording of the license and the list of + copyright holders. +*/ +/*! + * \file + * \copydoc Opm::BaseVanguard + */ +#ifndef EWOMS_BASE_VANGUARD_HH +#define EWOMS_BASE_VANGUARD_HH + +#include +#include + +#include + +#if HAVE_DUNE_FEM +#include +#endif + +#include +#include + +BEGIN_PROPERTIES + +NEW_PROP_TAG(Grid); +NEW_PROP_TAG(Vanguard); +NEW_PROP_TAG(GridView); +NEW_PROP_TAG(GridPart); +NEW_PROP_TAG(GridViewLevel); +NEW_PROP_TAG(GridFile); +NEW_PROP_TAG(GridGlobalRefinements); +NEW_PROP_TAG(Simulator); + +END_PROPERTIES + +namespace Opm { + +/*! + * \brief Provides the base class for most (all?) simulator vanguards. + */ +template +class BaseVanguard +{ + typedef typename GET_PROP_TYPE(TypeTag, Simulator) Simulator; + typedef typename GET_PROP_TYPE(TypeTag, Grid) Grid; + typedef typename GET_PROP_TYPE(TypeTag, GridView) GridView; + typedef typename GET_PROP_TYPE(TypeTag, Vanguard) Implementation; + +#if HAVE_DUNE_FEM + typedef typename GET_PROP_TYPE(TypeTag, GridPart) GridPart; +#endif + +public: + BaseVanguard(Simulator& simulator) + : simulator_(simulator) + {} + + BaseVanguard(const BaseVanguard&) = delete; + + /*! + * \brief Returns a reference to the grid view to be used. + */ + const GridView& gridView() const + { return *gridView_; } + +#if HAVE_DUNE_FEM + /*! + * \brief Returns a reference to the grid part to be used. + */ + const GridPart& gridPart() const + { return *gridPart_; } + + /*! + * \brief Returns a reference to the grid part to be used. + */ + GridPart& gridPart() + { return *gridPart_; } +#endif + + /*! + * \brief Returns the number of times the grid has been changed since its creation. + * + * This basically says how often the grid has been adapted in the current simulation + * run. + */ + int gridSequenceNumber () const + { +#if HAVE_DUNE_FEM + typedef Dune::Fem::DofManager< Grid > FemDofManager; + return FemDofManager::instance( asImp_().grid() ).sequence(); +#else + return 0; // return the same sequence number >= 0 means the grid never changes +#endif + } + + + /*! + * \brief Distribute the grid (and attached data) over all + * processes. + */ + void loadBalance() + { + asImp_().grid().loadBalance(); + updateGridView_(); + } + +protected: + // this method should be called after the grid has been allocated + void finalizeInit_() + { + updateGridView_(); + } + + void updateGridView_() + { +#if HAVE_DUNE_FEM + // first delete old grid part + // this is due to a bug in dune-fem (dangling reference) + gridPart_.reset(); + gridPart_.reset(new GridPart(asImp_().grid())); + gridView_.reset(new GridView(static_cast(*gridPart_))); + assert(gridView_->size(0) == asImp_().grid().leafGridView().size(0)); +#else + gridView_.reset(new GridView(asImp_().grid().leafGridView())); +#endif + } + +private: + Implementation& asImp_() + { return *static_cast(this); } + + const Implementation& asImp_() const + { return *static_cast(this); } + + Simulator& simulator_; +#if HAVE_DUNE_FEM + std::unique_ptr gridPart_; +#endif + std::unique_ptr gridView_; +}; + +} // namespace Opm + +#endif diff --git a/opm/models/io/cubegridvanguard.hh b/opm/models/io/cubegridvanguard.hh new file mode 100644 index 000000000..46ee7f710 --- /dev/null +++ b/opm/models/io/cubegridvanguard.hh @@ -0,0 +1,156 @@ +// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- +// vi: set et ts=4 sw=4 sts=4: +/* + This file is part of the Open Porous Media project (OPM). + + OPM is free software: you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation, either version 2 of the License, or + (at your option) any later version. + + OPM is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU General Public License + along with OPM. If not, see . + + Consult the COPYING file in the top-level source directory of this + module for the precise wording of the license and the list of + copyright holders. +*/ +/*! + * \file + * \copydoc Opm::CubeGridVanguard + */ +#ifndef EWOMS_CUBE_GRID_VANGUARD_HH +#define EWOMS_CUBE_GRID_VANGUARD_HH + +#include +#include +#include +#include + +#include + +#include + +#include + +BEGIN_PROPERTIES + +NEW_PROP_TAG(Scalar); +NEW_PROP_TAG(Grid); + +NEW_PROP_TAG(DomainSizeX); +NEW_PROP_TAG(DomainSizeY); +NEW_PROP_TAG(DomainSizeZ); + +NEW_PROP_TAG(CellsX); +NEW_PROP_TAG(CellsY); +NEW_PROP_TAG(CellsZ); + +NEW_PROP_TAG(GridGlobalRefinements); + +END_PROPERTIES + +namespace Opm { + +/*! + * \brief Provides a simulator vanguad which creates a regular grid made of + * quadrilaterals. + * + * A quadrilateral is a line segment in 1D, a rectangle in 2D and a + * cube in 3D. + */ +template +class CubeGridVanguard : public BaseVanguard +{ + typedef BaseVanguard ParentType; + typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar; + typedef typename GET_PROP_TYPE(TypeTag, Simulator) Simulator; + typedef typename GET_PROP_TYPE(TypeTag, Grid) Grid; + + typedef Dune::shared_ptr GridPointer; + typedef typename Grid::ctype CoordScalar; + enum { dimWorld = Grid::dimensionworld }; + typedef Dune::FieldVector GlobalPosition; + +public: + /*! + * \brief Register all run-time parameters for the simulator vanguad. + */ + static void registerParameters() + { + EWOMS_REGISTER_PARAM(TypeTag, unsigned, GridGlobalRefinements, + "The number of global refinements of the grid " + "executed after it was loaded"); + EWOMS_REGISTER_PARAM(TypeTag, Scalar, DomainSizeX, + "The size of the domain in x direction"); + EWOMS_REGISTER_PARAM(TypeTag, unsigned, CellsX, + "The number of intervalls in x direction"); + if (dimWorld > 1) { + EWOMS_REGISTER_PARAM(TypeTag, Scalar, DomainSizeY, + "The size of the domain in y direction"); + EWOMS_REGISTER_PARAM(TypeTag, unsigned, CellsY, + "The number of intervalls in y direction"); + } + if (dimWorld > 2) { + EWOMS_REGISTER_PARAM(TypeTag, Scalar, DomainSizeZ, + "The size of the domain in z direction"); + EWOMS_REGISTER_PARAM(TypeTag, unsigned, CellsZ, + "The number of intervalls in z direction"); + } + } + + /*! + * \brief Create the grid + */ + CubeGridVanguard(Simulator& simulator) + : ParentType(simulator) + { + std::array cellRes; + GlobalPosition upperRight(0.0); + GlobalPosition lowerLeft(0.0); + + for (unsigned i = 0; i < dimWorld; ++i) + cellRes[i] = 0; + + upperRight[0] = EWOMS_GET_PARAM(TypeTag, Scalar, DomainSizeX); + cellRes[0] = EWOMS_GET_PARAM(TypeTag, unsigned, CellsX); + if (dimWorld > 1) { + upperRight[1] = EWOMS_GET_PARAM(TypeTag, Scalar, DomainSizeY); + cellRes[1] = EWOMS_GET_PARAM(TypeTag, unsigned, CellsY); + } + if (dimWorld > 2) { + upperRight[2] = EWOMS_GET_PARAM(TypeTag, Scalar, DomainSizeZ); + cellRes[2] = EWOMS_GET_PARAM(TypeTag, unsigned, CellsZ); + } + + unsigned numRefinements = EWOMS_GET_PARAM(TypeTag, unsigned, GridGlobalRefinements); + cubeGrid_ = Dune::StructuredGridFactory::createCubeGrid(lowerLeft, upperRight, cellRes); + cubeGrid_->globalRefine(static_cast(numRefinements)); + + this->finalizeInit_(); + } + + /*! + * \brief Returns a reference to the grid. + */ + Grid& grid() + { return *cubeGrid_; } + + /*! + * \brief Returns a reference to the grid. + */ + const Grid& grid() const + { return *cubeGrid_; } + +protected: + GridPointer cubeGrid_; +}; + +} // namespace Opm + +#endif diff --git a/opm/models/io/dgfvanguard.hh b/opm/models/io/dgfvanguard.hh new file mode 100644 index 000000000..4965ab407 --- /dev/null +++ b/opm/models/io/dgfvanguard.hh @@ -0,0 +1,208 @@ +// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- +// vi: set et ts=4 sw=4 sts=4: +/* + This file is part of the Open Porous Media project (OPM). + + OPM is free software: you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation, either version 2 of the License, or + (at your option) any later version. + + OPM is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU General Public License + along with OPM. If not, see . + + Consult the COPYING file in the top-level source directory of this + module for the precise wording of the license and the list of + copyright holders. +*/ +/*! + * \file + * \copydoc Opm::DgfVanguard + */ +#ifndef EWOMS_DGF_GRID_VANGUARD_HH +#define EWOMS_DGF_GRID_VANGUARD_HH + +#include +#include +#include + +#include +#include +#include + + +#include +#include + +BEGIN_PROPERTIES + +NEW_PROP_TAG(Grid); +NEW_PROP_TAG(GridFile); +NEW_PROP_TAG(Vanguard); +NEW_PROP_TAG(GridGlobalRefinements); +NEW_PROP_TAG(Scalar); +NEW_PROP_TAG(Simulator); + +END_PROPERTIES + +namespace Opm { + +/*! + * \brief Provides a simulator vanguard which creates a grid by parsing a Dune Grid + * Format (DGF) file. + */ +template +class DgfVanguard : public BaseVanguard +{ + typedef BaseVanguard ParentType; + typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar; + typedef typename GET_PROP_TYPE(TypeTag, Simulator) Simulator; + typedef typename GET_PROP_TYPE(TypeTag, Grid) Grid; + typedef Opm::FractureMapper FractureMapper; + + typedef std::unique_ptr< Grid > GridPointer; + +public: + /*! + * \brief Register all run-time parameters for the DGF simulator vanguard. + */ + static void registerParameters() + { + EWOMS_REGISTER_PARAM(TypeTag, std::string, GridFile, + "The file name of the DGF file to load"); + EWOMS_REGISTER_PARAM(TypeTag, unsigned, GridGlobalRefinements, + "The number of global refinements of the grid " + "executed after it was loaded"); + } + + /*! + * \brief Load the grid from the file. + */ + DgfVanguard(Simulator& simulator) + : ParentType(simulator) + { + const std::string dgfFileName = EWOMS_GET_PARAM(TypeTag, std::string, GridFile); + unsigned numRefinments = EWOMS_GET_PARAM(TypeTag, unsigned, GridGlobalRefinements); + + { + // create DGF GridPtr from a dgf file + Dune::GridPtr< Grid > dgfPointer( dgfFileName ); + + // this is only implemented for 2d currently + addFractures_( dgfPointer ); + + // store pointer to dune grid + gridPtr_.reset( dgfPointer.release() ); + } + + if (numRefinments > 0) + gridPtr_->globalRefine(static_cast(numRefinments)); + + this->finalizeInit_(); + } + + /*! + * \brief Returns a reference to the grid. + */ + Grid& grid() + { return *gridPtr_; } + + /*! + * \brief Returns a reference to the grid. + */ + const Grid& grid() const + { return *gridPtr_; } + + /*! + * \brief Distributes the grid on all processes of a parallel + * computation. + * + * This grid manager plays nice and also distributes the data of + * the DGF... + */ + void loadBalance() + { gridPtr_->loadBalance(); } + + /*! + * \brief Returns the fracture mapper + * + * The fracture mapper determines the topology of the fractures. + */ + FractureMapper& fractureMapper() + { return fractureMapper_; } + + /*! + * \brief Returns the fracture mapper + * + * The fracture mapper determines the topology of the fractures. + */ + const FractureMapper& fractureMapper() const + { return fractureMapper_; } + +protected: + void addFractures_(Dune::GridPtr& dgfPointer) + { + typedef typename Grid::LevelGridView LevelGridView; + + // check if fractures are available (only 2d currently) + if (dgfPointer.nofParameters(static_cast(Grid::dimension)) == 0) + return; + + LevelGridView gridView = dgfPointer->levelGridView(/*level=*/0); + const unsigned edgeCodim = Grid::dimension - 1; + +#if DUNE_VERSION_NEWER(DUNE_GRID, 2,6) + typedef Dune::MultipleCodimMultipleGeomTypeMapper VertexMapper; + VertexMapper vertexMapper(gridView, Dune::mcmgVertexLayout()); +#else + typedef Dune::MultipleCodimMultipleGeomTypeMapper VertexMapper; + VertexMapper vertexMapper(gridView); +#endif + + // first create a map of the dune to ART vertex indices + auto eIt = gridView.template begin(); + const auto eEndIt = gridView.template end(); + for (; eIt != eEndIt; ++eIt) { + const auto& element = *eIt; + const auto& refElem = + Dune::ReferenceElements::general(element.type()); + + const int edges = refElem.size( edgeCodim ); + for (int edge = 0; edge < edges; ++edge) { + const int vertices = refElem.size(edge, edgeCodim, Grid::dimension); + std::vector vertexIndices; + vertexIndices.reserve(Grid::dimension); + for (int vx = 0; vx < vertices; ++vx) { + // get local vertex number from edge + const int localVx = refElem.subEntity(edge, edgeCodim, vx, Grid::dimension); + + // get vertex + const auto vertex = element.template subEntity(localVx); + + // if vertex has parameter 1 insert as a fracture vertex + if (dgfPointer.parameters( vertex )[ 0 ] > 0) + vertexIndices.push_back( + static_cast(vertexMapper.subIndex(element, + static_cast(localVx), + Grid::dimension))); + } + // if 2 vertices have been found with flag 1 insert a fracture edge + if (static_cast(vertexIndices.size()) == Grid::dimension) + fractureMapper_.addFractureEdge(vertexIndices[0], vertexIndices[1]); + } + } + } + +private: + GridPointer gridPtr_; + FractureMapper fractureMapper_; +}; + +} // namespace Opm + +#endif diff --git a/opm/models/io/restart.hh b/opm/models/io/restart.hh new file mode 100644 index 000000000..6d22d465c --- /dev/null +++ b/opm/models/io/restart.hh @@ -0,0 +1,275 @@ +// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- +// vi: set et ts=4 sw=4 sts=4: +/* + This file is part of the Open Porous Media project (OPM). + + OPM is free software: you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation, either version 2 of the License, or + (at your option) any later version. + + OPM is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU General Public License + along with OPM. If not, see . + + Consult the COPYING file in the top-level source directory of this + module for the precise wording of the license and the list of + copyright holders. +*/ +/*! + * \file + * \copydoc Opm::Restart + */ +#ifndef EWOMS_RESTART_HH +#define EWOMS_RESTART_HH + +#include +#include +#include +#include + +namespace Opm { + +/*! + * \brief Load or save a state of a problem to/from the harddisk. + */ +class Restart +{ + /*! + * \brief Create a magic cookie for restart files, so that it is + * unlikely to load a restart file for an incorrectly. + */ + template + static const std::string magicRestartCookie_(const GridView& gridView) + { + static const std::string gridName = "blubb"; // gridView.grid().name(); + static const int dim = GridView::dimension; + + int numVertices = gridView.size(dim); + int numElements = gridView.size(0); + int numEdges = gridView.size(dim - 1); + int numCPUs = gridView.comm().size(); + int rank = gridView.comm().rank(); + + std::ostringstream oss; + oss << "eWoms restart file: " + << "gridName='" << gridName << "' " + << "numCPUs=" << numCPUs << " " + << "myRank=" << rank << " " + << "numElements=" << numElements << " " + << "numEdges=" << numEdges << " " + << "numVertices=" << numVertices; + return oss.str(); + } + + /*! + * \brief Return the restart file name. + */ + template + static const std::string restartFileName_(const GridView& gridView, + const std::string& outputDir, + const std::string& simName, + Scalar t) + { + std::string dir = outputDir; + if (dir == ".") + dir = ""; + else if (!dir.empty() && dir.back() != '/') + dir += "/"; + + int rank = gridView.comm().rank(); + std::ostringstream oss; + oss << dir << simName << "_time=" << t << "_rank=" << rank << ".ers"; + return oss.str(); + } + +public: + /*! + * \brief Returns the name of the file which is (de-)serialized. + */ + const std::string& fileName() const + { return fileName_; } + + /*! + * \brief Write the current state of the model to disk. + */ + template + void serializeBegin(Simulator& simulator) + { + const std::string magicCookie = magicRestartCookie_(simulator.gridView()); + fileName_ = restartFileName_(simulator.gridView(), + simulator.problem().outputDir(), + simulator.problem().name(), + simulator.time()); + + // open output file and write magic cookie + outStream_.open(fileName_.c_str()); + outStream_.precision(20); + + serializeSectionBegin(magicCookie); + serializeSectionEnd(); + } + + /*! + * \brief The output stream to write the serialized data. + */ + std::ostream& serializeStream() + { return outStream_; } + + /*! + * \brief Start a new section in the serialized output. + */ + void serializeSectionBegin(const std::string& cookie) + { outStream_ << cookie << "\n"; } + + /*! + * \brief End of a section in the serialized output. + */ + void serializeSectionEnd() + { outStream_ << "\n"; } + + /*! + * \brief Serialize all leaf entities of a codim in a gridView. + * + * The actual work is done by Serializer::serialize(Entity) + */ + template + void serializeEntities(Serializer& serializer, const GridView& gridView) + { + std::ostringstream oss; + oss << "Entities: Codim " << codim; + std::string cookie = oss.str(); + serializeSectionBegin(cookie); + + // write element data + typedef typename GridView::template Codim::Iterator Iterator; + + Iterator it = gridView.template begin(); + const Iterator& endIt = gridView.template end(); + for (; it != endIt; ++it) { + serializer.serializeEntity(outStream_, *it); + outStream_ << "\n"; + } + + serializeSectionEnd(); + } + + /*! + * \brief Finish the restart file. + */ + void serializeEnd() + { outStream_.close(); } + + /*! + * \brief Start reading a restart file at a certain simulated + * time. + */ + template + void deserializeBegin(Simulator& simulator, Scalar t) + { + fileName_ = restartFileName_(simulator.gridView(), simulator.problem().outputDir(), simulator.problem().name(), t); + + // open input file and read magic cookie + inStream_.open(fileName_.c_str()); + if (!inStream_.good()) { + throw std::runtime_error("Restart file '"+fileName_+"' could not be opened properly"); + } + + // make sure that we don't open an empty file + inStream_.seekg(0, std::ios::end); + auto pos = inStream_.tellg(); + if (pos == 0) { + throw std::runtime_error("Restart file '"+fileName_+"' is empty"); + } + inStream_.seekg(0, std::ios::beg); + + const std::string magicCookie = magicRestartCookie_(simulator.gridView()); + + deserializeSectionBegin(magicCookie); + deserializeSectionEnd(); + } + + /*! + * \brief The input stream to read the data which ought to be + * deserialized. + */ + std::istream& deserializeStream() + { return inStream_; } + + /*! + * \brief Start reading a new section of the restart file. + */ + void deserializeSectionBegin(const std::string& cookie) + { + if (!inStream_.good()) + throw std::runtime_error("Encountered unexpected EOF in restart file."); + std::string buf; + std::getline(inStream_, buf); + if (buf != cookie) + throw std::runtime_error("Could not start section '"+cookie+"'"); + } + + /*! + * \brief End of a section in the serialized output. + */ + void deserializeSectionEnd() + { + std::string dummy; + std::getline(inStream_, dummy); + for (unsigned i = 0; i < dummy.length(); ++i) { + if (!std::isspace(dummy[i])) { + throw std::logic_error("Encountered unread values while deserializing"); + } + } + } + + /*! + * \brief Deserialize all leaf entities of a codim in a grid. + * + * The actual work is done by Deserializer::deserialize(Entity) + */ + template + void deserializeEntities(Deserializer& deserializer, const GridView& gridView) + { + std::ostringstream oss; + oss << "Entities: Codim " << codim; + std::string cookie = oss.str(); + deserializeSectionBegin(cookie); + + std::string curLine; + + // read entity data + typedef typename GridView::template Codim::Iterator Iterator; + Iterator it = gridView.template begin(); + const Iterator& endIt = gridView.template end(); + for (; it != endIt; ++it) { + if (!inStream_.good()) { + throw std::runtime_error("Restart file is corrupted"); + } + + std::getline(inStream_, curLine); + std::istringstream curLineStream(curLine); + deserializer.deserializeEntity(curLineStream, *it); + } + + deserializeSectionEnd(); + } + + /*! + * \brief Stop reading the restart file. + */ + void deserializeEnd() + { inStream_.close(); } + +private: + std::string fileName_; + std::ifstream inStream_; + std::ofstream outStream_; +}; +} // namespace Opm + +#endif diff --git a/opm/models/io/simplexvanguard.hh b/opm/models/io/simplexvanguard.hh new file mode 100644 index 000000000..33891cd81 --- /dev/null +++ b/opm/models/io/simplexvanguard.hh @@ -0,0 +1,151 @@ +// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- +// vi: set et ts=4 sw=4 sts=4: +/* + This file is part of the Open Porous Media project (OPM). + + OPM is free software: you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation, either version 2 of the License, or + (at your option) any later version. + + OPM is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU General Public License + along with OPM. If not, see . + + Consult the COPYING file in the top-level source directory of this + module for the precise wording of the license and the list of + copyright holders. +*/ +/*! + * \file + * \copydoc Opm::SimplexGridVanguard + */ +#ifndef EWOMS_SIMPLEX_GRID_VANGUARD_HH +#define EWOMS_SIMPLEX_GRID_VANGUARD_HH + +#include +#include +#include + +#include +#include + +#include + +BEGIN_PROPERTIES + +NEW_PROP_TAG(Scalar); +NEW_PROP_TAG(Grid); + +NEW_PROP_TAG(DomainSizeX); +NEW_PROP_TAG(DomainSizeY); +NEW_PROP_TAG(DomainSizeZ); + +NEW_PROP_TAG(CellsX); +NEW_PROP_TAG(CellsY); +NEW_PROP_TAG(CellsZ); + +NEW_PROP_TAG(GridGlobalRefinements); + +END_PROPERTIES + +namespace Opm { +/*! + * \brief Provides a simulator vanguard which a creates regular grid made of simplices. + */ +template +class SimplexGridVanguard +{ + typedef BaseVanguard ParentType; + typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar; + typedef typename GET_PROP_TYPE(TypeTag, Simulator) Simulator; + typedef typename GET_PROP_TYPE(TypeTag, Grid) Grid; + + typedef Dune::shared_ptr GridPointer; + typedef typename Grid::ctype CoordScalar; + enum { dimWorld = Grid::dimensionworld }; + typedef Dune::FieldVector GlobalPosition; + +public: + /*! + * \brief Register all run-time parameters for the grid manager. + */ + static void registerParameters() + { + EWOMS_REGISTER_PARAM(TypeTag, unsigned, GridGlobalRefinements, + "The number of global refinements of the grid " + "executed after it was loaded"); + EWOMS_REGISTER_PARAM(TypeTag, Scalar, DomainSizeX, + "The size of the domain in x direction"); + EWOMS_REGISTER_PARAM(TypeTag, unsigned, CellsX, + "The number of intervalls in x direction"); + if (dimWorld > 1) { + EWOMS_REGISTER_PARAM(TypeTag, Scalar, DomainSizeY, + "The size of the domain in y direction"); + EWOMS_REGISTER_PARAM(TypeTag, unsigned, CellsY, + "The number of intervalls in y direction"); + } + if (dimWorld > 2) { + EWOMS_REGISTER_PARAM(TypeTag, Scalar, DomainSizeZ, + "The size of the domain in z direction"); + EWOMS_REGISTER_PARAM(TypeTag, unsigned, CellsZ, + "The number of intervalls in z direction"); + } + } + + /*! + * \brief Create the Grid + */ + SimplexGridVanguard(Simulator& simulator) + : ParentType(simulator) + { + Dune::array cellRes; + GlobalPosition upperRight; + GlobalPosition lowerLeft; + + lowerLeft[0] = 0.0; + upperRight[0] = EWOMS_GET_PARAM(TypeTag, Scalar, DomainSizeX); + cellRes[0] = EWOMS_GET_PARAM(TypeTag, unsigned, CellsX); + if (dimWorld > 1) { + lowerLeft[1] = 0.0; + upperRight[1] = EWOMS_GET_PARAM(TypeTag, Scalar, DomainSizeY); + cellRes[1] = EWOMS_GET_PARAM(TypeTag, unsigned, CellsY); + } + if (dimWorld > 2) { + lowerLeft[2] = 0.0; + upperRight[2] = EWOMS_GET_PARAM(TypeTag, Scalar, DomainSizeZ); + cellRes[2] = EWOMS_GET_PARAM(TypeTag, unsigned, CellsZ); + } + + simplexGrid_ = Dune::StructuredGridFactory::createSimplexGrid(lowerLeft, + upperRight, + cellRes); + + unsigned numRefinments = EWOMS_GET_PARAM(TypeTag, unsigned, GridGlobalRefinements); + simplexGrid_->globalRefine(numRefinments); + + this->finalizeInit_(); + } + + /*! + * \brief Returns a reference to the grid. + */ + Grid& grid() + { return simplexGrid_; } + + /*! + * \brief Returns a reference to the grid. + */ + const Grid& grid() const + { return *simplexGrid_; } + +private: + GridPointer simplexGrid_; +}; +} // namespace Opm + +#endif diff --git a/opm/models/io/structuredgridvanguard.hh b/opm/models/io/structuredgridvanguard.hh new file mode 100644 index 000000000..3e143737c --- /dev/null +++ b/opm/models/io/structuredgridvanguard.hh @@ -0,0 +1,199 @@ +// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- +// vi: set et ts=4 sw=4 sts=4: +/* + This file is part of the Open Porous Media project (OPM). + + OPM is free software: you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation, either version 2 of the License, or + (at your option) any later version. + + OPM is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU General Public License + along with OPM. If not, see . + + Consult the COPYING file in the top-level source directory of this + module for the precise wording of the license and the list of + copyright holders. +*/ +/*! + * \file + * \copydoc Opm::StructuredGridVanguard + */ +#ifndef EWOMS_STRUCTURED_GRID_VANGUARD_HH +#define EWOMS_STRUCTURED_GRID_VANGUARD_HH + +#include +#include +#include + +#include +#include + +#if HAVE_DUNE_ALUGRID +#include +#include +#endif + +#include +#include + +#include +#include + +namespace Opm { + +template +class StructuredGridVanguard; + +} // namespace Opm + +BEGIN_PROPERTIES + +NEW_TYPE_TAG(StructuredGridVanguard); + +// declare the properties required by the for the structured grid simulator vanguard +NEW_PROP_TAG(Grid); +NEW_PROP_TAG(Scalar); + +NEW_PROP_TAG(DomainSizeX); +NEW_PROP_TAG(DomainSizeY); +NEW_PROP_TAG(DomainSizeZ); + +NEW_PROP_TAG(CellsX); +NEW_PROP_TAG(CellsY); +NEW_PROP_TAG(CellsZ); + +NEW_PROP_TAG(GridGlobalRefinements); + +// GRIDDIM is only set by the finger problem +#ifndef GRIDDIM +static const int dim = 2; +#else +static const int dim = GRIDDIM; +#endif + +// set the Grid and Vanguard properties +#if HAVE_DUNE_ALUGRID +SET_TYPE_PROP(StructuredGridVanguard, Grid, Dune::ALUGrid< dim, dim, Dune::cube, Dune::nonconforming >); +#else +SET_TYPE_PROP(StructuredGridVanguard, Grid, Dune::YaspGrid< dim >); +#endif + +SET_TYPE_PROP(StructuredGridVanguard, Vanguard, Opm::StructuredGridVanguard); + +END_PROPERTIES + +namespace Opm { + +/*! + * \ingroup TestProblems + * + * \brief Helper class for grid instantiation of the lens problem. + */ +template +class StructuredGridVanguard : public BaseVanguard +{ + typedef BaseVanguard ParentType; + typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar; + typedef typename GET_PROP_TYPE(TypeTag, Simulator) Simulator; + typedef typename GET_PROP_TYPE(TypeTag, Grid) Grid; + + typedef std::unique_ptr GridPointer; + + static const int dim = Grid::dimension; + +public: + /*! + * \brief Register all run-time parameters for the structured grid simulator vanguard. + */ + static void registerParameters() + { + EWOMS_REGISTER_PARAM(TypeTag, unsigned, GridGlobalRefinements, + "The number of global refinements of the grid " + "executed after it was loaded"); + EWOMS_REGISTER_PARAM(TypeTag, Scalar, DomainSizeX, + "The size of the domain in x direction"); + EWOMS_REGISTER_PARAM(TypeTag, unsigned, CellsX, + "The number of intervalls in x direction"); + if (dim > 1) { + EWOMS_REGISTER_PARAM(TypeTag, Scalar, DomainSizeY, + "The size of the domain in y direction"); + EWOMS_REGISTER_PARAM(TypeTag, unsigned, CellsY, + "The number of intervalls in y direction"); + } + if (dim > 2) { + EWOMS_REGISTER_PARAM(TypeTag, Scalar, DomainSizeZ, + "The size of the domain in z direction"); + EWOMS_REGISTER_PARAM(TypeTag, unsigned, CellsZ, + "The number of intervalls in z direction"); + } + } + + /*! + * \brief Create the grid for the lens problem + */ + StructuredGridVanguard(Simulator& simulator) + : ParentType(simulator) + { + Dune::FieldVector cellRes; + + typedef double GridScalar; + Dune::FieldVector upperRight; + Dune::FieldVector lowerLeft( 0 ); + + upperRight[0] = EWOMS_GET_PARAM(TypeTag, Scalar, DomainSizeX); + upperRight[1] = EWOMS_GET_PARAM(TypeTag, Scalar, DomainSizeY); + + cellRes[0] = EWOMS_GET_PARAM(TypeTag, unsigned, CellsX); + cellRes[1] = EWOMS_GET_PARAM(TypeTag, unsigned, CellsY); + if (dim == 3) { + upperRight[2] = EWOMS_GET_PARAM(TypeTag, Scalar, DomainSizeZ); + cellRes[2] = EWOMS_GET_PARAM(TypeTag, unsigned, CellsZ); + } + + std::stringstream dgffile; + dgffile << "DGF" << std::endl; + dgffile << "INTERVAL" << std::endl; + dgffile << lowerLeft << std::endl; + dgffile << upperRight << std::endl; + dgffile << cellRes << std::endl; + dgffile << "#" << std::endl; + dgffile << "GridParameter" << std::endl; + dgffile << "overlap 1" << std::endl; + dgffile << "#" << std::endl; + dgffile << "Simplex" << std::endl; + dgffile << "#" << std::endl; + + // use DGF parser to create a grid from interval block + gridPtr_.reset( Dune::GridPtr< Grid >( dgffile ).release() ); + + unsigned numRefinements = EWOMS_GET_PARAM(TypeTag, unsigned, GridGlobalRefinements); + gridPtr_->globalRefine(static_cast(numRefinements)); + + this->finalizeInit_(); + } + + /*! + * \brief Return a reference to the grid object. + */ + Grid& grid() + { return *gridPtr_; } + + /*! + * \brief Return a constant reference to the grid object. + */ + const Grid& grid() const + { return *gridPtr_; } + +private: + GridPointer gridPtr_; +}; + +} // namespace Opm + +#endif diff --git a/opm/models/io/vtkblackoilenergymodule.hh b/opm/models/io/vtkblackoilenergymodule.hh new file mode 100644 index 000000000..54fed458c --- /dev/null +++ b/opm/models/io/vtkblackoilenergymodule.hh @@ -0,0 +1,233 @@ +// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- +// vi: set et ts=4 sw=4 sts=4: +/* + This file is part of the Open Porous Media project (OPM). + + OPM is free software: you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation, either version 2 of the License, or + (at your option) any later version. + + OPM is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU General Public License + along with OPM. If not, see . + + Consult the COPYING file in the top-level source directory of this + module for the precise wording of the license and the list of + copyright holders. +*/ +/*! + * \file + * \copydoc Opm::VtkBlackOilEnergyModule + */ +#ifndef EWOMS_VTK_BLACK_OIL_ENERGY_MODULE_HH +#define EWOMS_VTK_BLACK_OIL_ENERGY_MODULE_HH + +#include + +#include "vtkmultiwriter.hh" +#include "baseoutputmodule.hh" + +#include +#include +#include + +#include + +#include + +BEGIN_PROPERTIES + +// create new type tag for the VTK multi-phase output +NEW_TYPE_TAG(VtkBlackOilEnergy); + +// create the property tags needed for the energy module +NEW_PROP_TAG(VtkWriteRockInternalEnergy); +NEW_PROP_TAG(VtkWriteTotalThermalConductivity); +NEW_PROP_TAG(VtkWriteFluidInternalEnergies); +NEW_PROP_TAG(VtkWriteFluidEnthalpies); +NEW_PROP_TAG(VtkOutputFormat); +NEW_PROP_TAG(EnableVtkOutput); + +// set default values for what quantities to output +SET_BOOL_PROP(VtkBlackOilEnergy, VtkWriteRockInternalEnergy, true); +SET_BOOL_PROP(VtkBlackOilEnergy, VtkWriteTotalThermalConductivity, true); +SET_BOOL_PROP(VtkBlackOilEnergy, VtkWriteFluidInternalEnergies, true); +SET_BOOL_PROP(VtkBlackOilEnergy, VtkWriteFluidEnthalpies, true); + +END_PROPERTIES + +namespace Opm { +/*! + * \ingroup Vtk + * + * \brief VTK output module for the black oil model's energy related quantities. + */ +template +class VtkBlackOilEnergyModule : public BaseOutputModule +{ + typedef BaseOutputModule ParentType; + + typedef typename GET_PROP_TYPE(TypeTag, Simulator) Simulator; + typedef typename GET_PROP_TYPE(TypeTag, GridView) GridView; + typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar; + typedef typename GET_PROP_TYPE(TypeTag, Evaluation) Evaluation; + typedef typename GET_PROP_TYPE(TypeTag, ElementContext) ElementContext; + + static const int vtkFormat = GET_PROP_VALUE(TypeTag, VtkOutputFormat); + typedef Opm::VtkMultiWriter VtkMultiWriter; + + enum { enableEnergy = GET_PROP_VALUE(TypeTag, EnableEnergy) }; + enum { numPhases = GET_PROP_VALUE(TypeTag, NumPhases) }; + + typedef typename ParentType::ScalarBuffer ScalarBuffer; + typedef typename ParentType::PhaseBuffer PhaseBuffer; + +public: + VtkBlackOilEnergyModule(const Simulator& simulator) + : ParentType(simulator) + { } + + /*! + * \brief Register all run-time parameters for the multi-phase VTK output + * module. + */ + static void registerParameters() + { + if (!enableEnergy) + return; + + EWOMS_REGISTER_PARAM(TypeTag, bool, VtkWriteRockInternalEnergy, + "Include the volumetric internal energy of rock " + "in the VTK output files"); + EWOMS_REGISTER_PARAM(TypeTag, bool, VtkWriteTotalThermalConductivity, + "Include the total thermal conductivity of the medium and the fluids " + "in the VTK output files"); + EWOMS_REGISTER_PARAM(TypeTag, bool, VtkWriteFluidInternalEnergies, + "Include the internal energies of the fluids " + "in the VTK output files"); + EWOMS_REGISTER_PARAM(TypeTag, bool, VtkWriteFluidEnthalpies, + "Include the enthalpies of the fluids " + "in the VTK output files"); + } + + /*! + * \brief Allocate memory for the scalar fields we would like to + * write to the VTK file. + */ + void allocBuffers() + { + if (!EWOMS_GET_PARAM(TypeTag, bool, EnableVtkOutput)) + return; + + if (!enableEnergy) + return; + + if (rockInternalEnergyOutput_()) + this->resizeScalarBuffer_(rockInternalEnergy_); + if (totalThermalConductivityOutput_()) + this->resizeScalarBuffer_(totalThermalConductivity_); + if (fluidInternalEnergiesOutput_()) + this->resizePhaseBuffer_(fluidInternalEnergies_); + if (fluidEnthalpiesOutput_()) + this->resizePhaseBuffer_(fluidEnthalpies_); + } + + /*! + * \brief Modify the internal buffers according to the intensive quantities relevant for + * an element + */ + void processElement(const ElementContext& elemCtx) + { + if (!EWOMS_GET_PARAM(TypeTag, bool, EnableVtkOutput)) + return; + + if (!enableEnergy) + return; + + for (unsigned dofIdx = 0; dofIdx < elemCtx.numPrimaryDof(/*timeIdx=*/0); ++dofIdx) { + const auto& intQuants = elemCtx.intensiveQuantities(dofIdx, /*timeIdx=*/0); + unsigned globalDofIdx = elemCtx.globalSpaceIndex(dofIdx, /*timeIdx=*/0); + + if (rockInternalEnergyOutput_()) + rockInternalEnergy_[globalDofIdx] = + Opm::scalarValue(intQuants.rockInternalEnergy()); + + if (totalThermalConductivityOutput_()) + totalThermalConductivity_[globalDofIdx] = + Opm::scalarValue(intQuants.totalThermalConductivity()); + + for (int phaseIdx = 0; phaseIdx < numPhases; ++ phaseIdx) { + if (fluidInternalEnergiesOutput_()) + fluidInternalEnergies_[phaseIdx][globalDofIdx] = + Opm::scalarValue(intQuants.fluidState().internalEnergy(phaseIdx)); + + if (fluidEnthalpiesOutput_()) + fluidEnthalpies_[phaseIdx][globalDofIdx] = + Opm::scalarValue(intQuants.fluidState().enthalpy(phaseIdx)); + } + } + } + + /*! + * \brief Add all buffers to the VTK output writer. + */ + void commitBuffers(BaseOutputWriter& baseWriter) + { + VtkMultiWriter *vtkWriter = dynamic_cast(&baseWriter); + if (!vtkWriter) + return; + + if (!enableEnergy) + return; + + if (rockInternalEnergyOutput_()) + this->commitScalarBuffer_(baseWriter, "volumetric internal energy rock", rockInternalEnergy_); + + if (totalThermalConductivityOutput_()) + this->commitScalarBuffer_(baseWriter, "total thermal conductivity", totalThermalConductivity_); + + if (fluidInternalEnergiesOutput_()) + this->commitPhaseBuffer_(baseWriter, "internal energy_%s", fluidInternalEnergies_); + + if (fluidEnthalpiesOutput_()) + this->commitPhaseBuffer_(baseWriter, "enthalpy_%s", fluidEnthalpies_); + } + +private: + static bool rockInternalEnergyOutput_() + { + static bool val = EWOMS_GET_PARAM(TypeTag, bool, VtkWriteRockInternalEnergy); + return val; + } + + static bool totalThermalConductivityOutput_() + { + static bool val = EWOMS_GET_PARAM(TypeTag, bool, VtkWriteTotalThermalConductivity); + return val; + } + + static bool fluidInternalEnergiesOutput_() + { + static bool val = EWOMS_GET_PARAM(TypeTag, bool, VtkWriteFluidInternalEnergies); + return val; + } + + static bool fluidEnthalpiesOutput_() + { + static bool val = EWOMS_GET_PARAM(TypeTag, bool, VtkWriteFluidEnthalpies); + return val; + } + + ScalarBuffer rockInternalEnergy_; + ScalarBuffer totalThermalConductivity_; + PhaseBuffer fluidInternalEnergies_; + PhaseBuffer fluidEnthalpies_; +}; +} // namespace Opm + +#endif diff --git a/opm/models/io/vtkblackoilmodule.hh b/opm/models/io/vtkblackoilmodule.hh new file mode 100644 index 000000000..ec0d251ce --- /dev/null +++ b/opm/models/io/vtkblackoilmodule.hh @@ -0,0 +1,390 @@ +// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- +// vi: set et ts=4 sw=4 sts=4: +/* + This file is part of the Open Porous Media project (OPM). + + OPM is free software: you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation, either version 2 of the License, or + (at your option) any later version. + + OPM is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU General Public License + along with OPM. If not, see . + + Consult the COPYING file in the top-level source directory of this + module for the precise wording of the license and the list of + copyright holders. +*/ +/*! + * \file + * \copydoc Opm::VtkBlackOilModule + */ +#ifndef EWOMS_VTK_BLACK_OIL_MODULE_HH +#define EWOMS_VTK_BLACK_OIL_MODULE_HH + +#include + +#include "vtkmultiwriter.hh" +#include "baseoutputmodule.hh" + +#include +#include +#include + +#include + +#include + +BEGIN_PROPERTIES + +// create new type tag for the VTK multi-phase output +NEW_TYPE_TAG(VtkBlackOil); + +// create the property tags needed for the multi phase module +NEW_PROP_TAG(EnableVtkOutput); +NEW_PROP_TAG(VtkOutputFormat); +NEW_PROP_TAG(VtkWriteGasDissolutionFactor); +NEW_PROP_TAG(VtkWriteOilVaporizationFactor); +NEW_PROP_TAG(VtkWriteOilFormationVolumeFactor); +NEW_PROP_TAG(VtkWriteGasFormationVolumeFactor); +NEW_PROP_TAG(VtkWriteWaterFormationVolumeFactor); +NEW_PROP_TAG(VtkWriteOilSaturationPressure); +NEW_PROP_TAG(VtkWriteGasSaturationPressure); +NEW_PROP_TAG(VtkWriteSaturationRatios); +NEW_PROP_TAG(VtkWriteSaturatedOilGasDissolutionFactor); +NEW_PROP_TAG(VtkWriteSaturatedGasOilVaporizationFactor); +NEW_PROP_TAG(VtkWritePrimaryVarsMeaning); + +// set default values for what quantities to output +SET_BOOL_PROP(VtkBlackOil, VtkWriteGasDissolutionFactor, false); +SET_BOOL_PROP(VtkBlackOil, VtkWriteOilVaporizationFactor, false); +SET_BOOL_PROP(VtkBlackOil, VtkWriteOilFormationVolumeFactor, false); +SET_BOOL_PROP(VtkBlackOil, VtkWriteGasFormationVolumeFactor, false); +SET_BOOL_PROP(VtkBlackOil, VtkWriteWaterFormationVolumeFactor, false); +SET_BOOL_PROP(VtkBlackOil, VtkWriteOilSaturationPressure, false); +SET_BOOL_PROP(VtkBlackOil, VtkWriteGasSaturationPressure, false); +SET_BOOL_PROP(VtkBlackOil, VtkWriteSaturationRatios, false); +SET_BOOL_PROP(VtkBlackOil, VtkWriteSaturatedOilGasDissolutionFactor, false); +SET_BOOL_PROP(VtkBlackOil, VtkWriteSaturatedGasOilVaporizationFactor, false); +SET_BOOL_PROP(VtkBlackOil, VtkWritePrimaryVarsMeaning, false); +END_PROPERTIES + +namespace Opm { +/*! + * \ingroup Vtk + * + * \brief VTK output module for the black oil model's parameters. + */ +template +class VtkBlackOilModule : public BaseOutputModule +{ + typedef BaseOutputModule ParentType; + + typedef typename GET_PROP_TYPE(TypeTag, Simulator) Simulator; + typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar; + typedef typename GET_PROP_TYPE(TypeTag, Evaluation) Evaluation; + typedef typename GET_PROP_TYPE(TypeTag, ElementContext) ElementContext; + + typedef typename GET_PROP_TYPE(TypeTag, GridView) GridView; + typedef typename GET_PROP_TYPE(TypeTag, FluidSystem) FluidSystem; + + static const int vtkFormat = GET_PROP_VALUE(TypeTag, VtkOutputFormat); + typedef Opm::VtkMultiWriter VtkMultiWriter; + + enum { oilPhaseIdx = FluidSystem::oilPhaseIdx }; + enum { gasPhaseIdx = FluidSystem::gasPhaseIdx }; + enum { waterPhaseIdx = FluidSystem::waterPhaseIdx }; + + enum { gasCompIdx = FluidSystem::gasCompIdx }; + enum { oilCompIdx = FluidSystem::oilCompIdx }; + enum { waterCompIdx = FluidSystem::waterCompIdx }; + + typedef typename ParentType::ScalarBuffer ScalarBuffer; + +public: + VtkBlackOilModule(const Simulator& simulator) + : ParentType(simulator) + { } + + /*! + * \brief Register all run-time parameters for the multi-phase VTK output + * module. + */ + static void registerParameters() + { + EWOMS_REGISTER_PARAM(TypeTag, bool, VtkWriteGasDissolutionFactor, + "Include the gas dissolution factor (R_s) of the observed oil " + "in the VTK output files"); + EWOMS_REGISTER_PARAM(TypeTag, bool, VtkWriteOilVaporizationFactor, + "Include the oil vaporization factor (R_v) of the observed gas " + "in the VTK output files"); + EWOMS_REGISTER_PARAM(TypeTag, bool, VtkWriteOilFormationVolumeFactor, + "Include the oil formation volume factor (B_o) in the " + "VTK output files"); + EWOMS_REGISTER_PARAM(TypeTag, bool, VtkWriteGasFormationVolumeFactor, + "Include the gas formation volume factor (B_g) in the " + "VTK output files"); + EWOMS_REGISTER_PARAM(TypeTag, bool, VtkWriteWaterFormationVolumeFactor, + "Include the water formation volume factor (B_w) in the " + "VTK output files"); + EWOMS_REGISTER_PARAM(TypeTag, bool, VtkWriteOilSaturationPressure, + "Include the saturation pressure of oil (p_o,sat) in the " + "VTK output files"); + EWOMS_REGISTER_PARAM(TypeTag, bool, VtkWriteGasSaturationPressure, + "Include the saturation pressure of gas (p_g,sat) in the " + "VTK output files"); + EWOMS_REGISTER_PARAM(TypeTag, bool, VtkWriteSaturatedOilGasDissolutionFactor, + "Include the gas dissolution factor (R_s,sat) of gas saturated " + "oil in the VTK output files"); + EWOMS_REGISTER_PARAM(TypeTag, bool, VtkWriteSaturatedGasOilVaporizationFactor, + "Include the oil vaporization factor (R_v,sat) of oil saturated " + "gas in the VTK output files"); + EWOMS_REGISTER_PARAM(TypeTag, bool, VtkWriteSaturationRatios, + "Write the ratio of the actually and maximum dissolved component of " + "the mixtures"); + EWOMS_REGISTER_PARAM(TypeTag, bool, VtkWritePrimaryVarsMeaning, + "Include how the primary variables should be interpreted to the " + "VTK output files"); + } + + /*! + * \brief Allocate memory for the scalar fields we would like to + * write to the VTK file. + */ + void allocBuffers() + { + if (gasDissolutionFactorOutput_()) + this->resizeScalarBuffer_(gasDissolutionFactor_); + if (oilVaporizationFactorOutput_()) + this->resizeScalarBuffer_(oilVaporizationFactor_); + if (oilFormationVolumeFactorOutput_()) + this->resizeScalarBuffer_(oilFormationVolumeFactor_); + if (gasFormationVolumeFactorOutput_()) + this->resizeScalarBuffer_(gasFormationVolumeFactor_); + if (waterFormationVolumeFactorOutput_()) + this->resizeScalarBuffer_(waterFormationVolumeFactor_); + if (oilSaturationPressureOutput_()) + this->resizeScalarBuffer_(oilSaturationPressure_); + if (gasSaturationPressureOutput_()) + this->resizeScalarBuffer_(gasSaturationPressure_); + if (saturatedOilGasDissolutionFactorOutput_()) + this->resizeScalarBuffer_(saturatedOilGasDissolutionFactor_); + if (saturatedGasOilVaporizationFactorOutput_()) + this->resizeScalarBuffer_(saturatedGasOilVaporizationFactor_); + if (saturationRatiosOutput_()) { + this->resizeScalarBuffer_(oilSaturationRatio_); + this->resizeScalarBuffer_(gasSaturationRatio_); + } + if (primaryVarsMeaningOutput_()) + this->resizeScalarBuffer_(primaryVarsMeaning_); + } + + /*! + * \brief Modify the internal buffers according to the intensive quantities relevant for + * an element + */ + void processElement(const ElementContext& elemCtx) + { + if (!EWOMS_GET_PARAM(TypeTag, bool, EnableVtkOutput)) + return; + + for (unsigned dofIdx = 0; dofIdx < elemCtx.numPrimaryDof(/*timeIdx=*/0); ++dofIdx) { + const auto& fs = elemCtx.intensiveQuantities(dofIdx, /*timeIdx=*/0).fluidState(); + typedef typename std::remove_const::type>::type FluidState; + unsigned globalDofIdx = elemCtx.globalSpaceIndex(dofIdx, /*timeIdx=*/0); + + const auto& primaryVars = elemCtx.primaryVars(dofIdx, /*timeIdx=*/0); + + unsigned pvtRegionIdx = elemCtx.primaryVars(dofIdx, /*timeIdx=*/0).pvtRegionIndex(); + Scalar SoMax = std::max(Opm::getValue(fs.saturation(oilPhaseIdx)), + elemCtx.problem().maxOilSaturation(globalDofIdx)); + Scalar x_oG = Opm::getValue(fs.moleFraction(oilPhaseIdx, gasCompIdx)); + Scalar x_gO = Opm::getValue(fs.moleFraction(gasPhaseIdx, oilCompIdx)); + Scalar X_oG = Opm::getValue(fs.massFraction(oilPhaseIdx, gasCompIdx)); + Scalar X_gO = Opm::getValue(fs.massFraction(gasPhaseIdx, oilCompIdx)); + Scalar Rs = FluidSystem::convertXoGToRs(X_oG, pvtRegionIdx); + Scalar Rv = FluidSystem::convertXgOToRv(X_gO, pvtRegionIdx); + + Scalar RsSat = + FluidSystem::template saturatedDissolutionFactor(fs, + oilPhaseIdx, + pvtRegionIdx, + SoMax); + Scalar X_oG_sat = FluidSystem::convertRsToXoG(RsSat, pvtRegionIdx); + Scalar x_oG_sat = FluidSystem::convertXoGToxoG(X_oG_sat, pvtRegionIdx); + + Scalar RvSat = + FluidSystem::template saturatedDissolutionFactor(fs, + gasPhaseIdx, + pvtRegionIdx, + SoMax); + Scalar X_gO_sat = FluidSystem::convertRvToXgO(RvSat, pvtRegionIdx); + Scalar x_gO_sat = FluidSystem::convertXgOToxgO(X_gO_sat, pvtRegionIdx); + + if (gasDissolutionFactorOutput_()) + gasDissolutionFactor_[globalDofIdx] = Rs; + if (oilVaporizationFactorOutput_()) + oilVaporizationFactor_[globalDofIdx] = Rv; + if (oilFormationVolumeFactorOutput_()) + oilFormationVolumeFactor_[globalDofIdx] = + 1.0/FluidSystem::template inverseFormationVolumeFactor(fs, oilPhaseIdx, pvtRegionIdx); + if (gasFormationVolumeFactorOutput_()) + gasFormationVolumeFactor_[globalDofIdx] = + 1.0/FluidSystem::template inverseFormationVolumeFactor(fs, gasPhaseIdx, pvtRegionIdx); + if (waterFormationVolumeFactorOutput_()) + waterFormationVolumeFactor_[globalDofIdx] = + 1.0/FluidSystem::template inverseFormationVolumeFactor(fs, waterPhaseIdx, pvtRegionIdx); + if (oilSaturationPressureOutput_()) + oilSaturationPressure_[globalDofIdx] = + FluidSystem::template saturationPressure(fs, oilPhaseIdx, pvtRegionIdx); + if (gasSaturationPressureOutput_()) + gasSaturationPressure_[globalDofIdx] = + FluidSystem::template saturationPressure(fs, gasPhaseIdx, pvtRegionIdx); + if (saturatedOilGasDissolutionFactorOutput_()) + saturatedOilGasDissolutionFactor_[globalDofIdx] = RsSat; + if (saturatedGasOilVaporizationFactorOutput_()) + saturatedGasOilVaporizationFactor_[globalDofIdx] = RvSat; + if (saturationRatiosOutput_()) { + if (x_oG_sat <= 0.0) + oilSaturationRatio_[globalDofIdx] = 1.0; + else + oilSaturationRatio_[globalDofIdx] = x_oG / x_oG_sat; + + if (x_gO_sat <= 0.0) + gasSaturationRatio_[globalDofIdx] = 1.0; + else + gasSaturationRatio_[globalDofIdx] = x_gO / x_gO_sat; + } + + if (primaryVarsMeaningOutput_()) + primaryVarsMeaning_[globalDofIdx] = + primaryVars.primaryVarsMeaning(); + } + } + + /*! + * \brief Add all buffers to the VTK output writer. + */ + void commitBuffers(BaseOutputWriter& baseWriter) + { + VtkMultiWriter *vtkWriter = dynamic_cast(&baseWriter); + if (!vtkWriter) + return; + + if (gasDissolutionFactorOutput_()) + this->commitScalarBuffer_(baseWriter, "R_s", gasDissolutionFactor_); + if (oilVaporizationFactorOutput_()) + this->commitScalarBuffer_(baseWriter, "R_v", oilVaporizationFactor_); + if (oilFormationVolumeFactorOutput_()) + this->commitScalarBuffer_(baseWriter, "B_o", oilFormationVolumeFactor_); + if (gasFormationVolumeFactorOutput_()) + this->commitScalarBuffer_(baseWriter, "B_g", gasFormationVolumeFactor_); + if (waterFormationVolumeFactorOutput_()) + this->commitScalarBuffer_(baseWriter, "B_w", waterFormationVolumeFactor_); + if (oilSaturationPressureOutput_()) + this->commitScalarBuffer_(baseWriter, "p_o,sat", oilSaturationPressure_); + if (gasSaturationPressureOutput_()) + this->commitScalarBuffer_(baseWriter, "p_g,sat", gasSaturationPressure_); + if (saturatedOilGasDissolutionFactorOutput_()) + this->commitScalarBuffer_(baseWriter, "R_s,sat", saturatedOilGasDissolutionFactor_); + if (saturatedGasOilVaporizationFactorOutput_()) + this->commitScalarBuffer_(baseWriter, "R_v,sat", saturatedGasOilVaporizationFactor_); + if (saturationRatiosOutput_()) { + this->commitScalarBuffer_(baseWriter, "saturation ratio_oil", oilSaturationRatio_); + this->commitScalarBuffer_(baseWriter, "saturation ratio_gas", gasSaturationRatio_); + } + + if (primaryVarsMeaningOutput_()) + this->commitScalarBuffer_(baseWriter, "primary vars meaning", primaryVarsMeaning_); + } + +private: + static bool gasDissolutionFactorOutput_() + { + static bool val = EWOMS_GET_PARAM(TypeTag, bool, VtkWriteGasDissolutionFactor); + return val; + } + + static bool oilVaporizationFactorOutput_() + { + static bool val = EWOMS_GET_PARAM(TypeTag, bool, VtkWriteOilVaporizationFactor); + return val; + } + + static bool oilFormationVolumeFactorOutput_() + { + static bool val = EWOMS_GET_PARAM(TypeTag, bool, VtkWriteOilFormationVolumeFactor); + return val; + } + + static bool gasFormationVolumeFactorOutput_() + { + static bool val = EWOMS_GET_PARAM(TypeTag, bool, VtkWriteGasFormationVolumeFactor); + return val; + } + + static bool waterFormationVolumeFactorOutput_() + { + static bool val = EWOMS_GET_PARAM(TypeTag, bool, VtkWriteWaterFormationVolumeFactor); + return val; + } + + static bool oilSaturationPressureOutput_() + { + static bool val = EWOMS_GET_PARAM(TypeTag, bool, VtkWriteOilSaturationPressure); + return val; + } + + static bool gasSaturationPressureOutput_() + { + static bool val = EWOMS_GET_PARAM(TypeTag, bool, VtkWriteGasSaturationPressure); + return val; + } + + static bool saturatedOilGasDissolutionFactorOutput_() + { + static bool val = EWOMS_GET_PARAM(TypeTag, bool, VtkWriteSaturatedOilGasDissolutionFactor); + return val; + } + + static bool saturatedGasOilVaporizationFactorOutput_() + { + static bool val = EWOMS_GET_PARAM(TypeTag, bool, VtkWriteSaturatedGasOilVaporizationFactor); + return val; + } + + static bool saturationRatiosOutput_() + { + static bool val = EWOMS_GET_PARAM(TypeTag, bool, VtkWriteSaturationRatios); + return val; + } + + static bool primaryVarsMeaningOutput_() + { + static bool val = EWOMS_GET_PARAM(TypeTag, bool, VtkWritePrimaryVarsMeaning); + return val; + } + + ScalarBuffer gasDissolutionFactor_; + ScalarBuffer oilVaporizationFactor_; + ScalarBuffer oilFormationVolumeFactor_; + ScalarBuffer gasFormationVolumeFactor_; + ScalarBuffer waterFormationVolumeFactor_; + ScalarBuffer oilSaturationPressure_; + ScalarBuffer gasSaturationPressure_; + + ScalarBuffer saturatedOilGasDissolutionFactor_; + ScalarBuffer saturatedGasOilVaporizationFactor_; + ScalarBuffer oilSaturationRatio_; + ScalarBuffer gasSaturationRatio_; + + ScalarBuffer primaryVarsMeaning_; +}; +} // namespace Opm + +#endif diff --git a/opm/models/io/vtkblackoilpolymermodule.hh b/opm/models/io/vtkblackoilpolymermodule.hh new file mode 100644 index 000000000..03a9ea78c --- /dev/null +++ b/opm/models/io/vtkblackoilpolymermodule.hh @@ -0,0 +1,271 @@ +// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- +// vi: set et ts=4 sw=4 sts=4: +/* + This file is part of the Open Porous Media project (OPM). + + OPM is free software: you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation, either version 2 of the License, or + (at your option) any later version. + + OPM is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU General Public License + along with OPM. If not, see . + + Consult the COPYING file in the top-level source directory of this + module for the precise wording of the license and the list of + copyright holders. +*/ +/*! + * \file + * \copydoc Opm::VtkBlackOilPolymerModule + */ +#ifndef EWOMS_VTK_BLACK_OIL_POLYMER_MODULE_HH +#define EWOMS_VTK_BLACK_OIL_POLYMER_MODULE_HH + +#include + +#include "vtkmultiwriter.hh" +#include "baseoutputmodule.hh" + +#include +#include +#include + +#include + +#include + +BEGIN_PROPERTIES + +// create new type tag for the VTK multi-phase output +NEW_TYPE_TAG(VtkBlackOilPolymer); + +// create the property tags needed for the polymer output module +NEW_PROP_TAG(EnablePolymer); +NEW_PROP_TAG(EnableVtkOutput); +NEW_PROP_TAG(VtkWritePolymerConcentration); +NEW_PROP_TAG(VtkWritePolymerDeadPoreVolume); +NEW_PROP_TAG(VtkWritePolymerAdsorption); +NEW_PROP_TAG(VtkWritePolymerRockDensity); +NEW_PROP_TAG(VtkWritePolymerViscosityCorrection); +NEW_PROP_TAG(VtkWriteWaterViscosityCorrection); + +// set default values for what quantities to output +SET_BOOL_PROP(VtkBlackOilPolymer, VtkWritePolymerConcentration, true); +SET_BOOL_PROP(VtkBlackOilPolymer, VtkWritePolymerDeadPoreVolume, true); +SET_BOOL_PROP(VtkBlackOilPolymer, VtkWritePolymerViscosityCorrection, true); +SET_BOOL_PROP(VtkBlackOilPolymer, VtkWriteWaterViscosityCorrection, true); +SET_BOOL_PROP(VtkBlackOilPolymer, VtkWritePolymerRockDensity, true); +SET_BOOL_PROP(VtkBlackOilPolymer, VtkWritePolymerAdsorption, true); + +END_PROPERTIES + +namespace Opm { +/*! + * \ingroup Vtk + * + * \brief VTK output module for the black oil model's polymer related quantities. + */ +template +class VtkBlackOilPolymerModule : public BaseOutputModule +{ + typedef BaseOutputModule ParentType; + + typedef typename GET_PROP_TYPE(TypeTag, Simulator) Simulator; + typedef typename GET_PROP_TYPE(TypeTag, GridView) GridView; + typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar; + typedef typename GET_PROP_TYPE(TypeTag, Evaluation) Evaluation; + typedef typename GET_PROP_TYPE(TypeTag, ElementContext) ElementContext; + + static const int vtkFormat = GET_PROP_VALUE(TypeTag, VtkOutputFormat); + typedef Opm::VtkMultiWriter VtkMultiWriter; + + enum { enablePolymer = GET_PROP_VALUE(TypeTag, EnablePolymer) }; + + typedef typename ParentType::ScalarBuffer ScalarBuffer; + +public: + VtkBlackOilPolymerModule(const Simulator& simulator) + : ParentType(simulator) + { } + + /*! + * \brief Register all run-time parameters for the multi-phase VTK output + * module. + */ + static void registerParameters() + { + if (!enablePolymer) + return; + + EWOMS_REGISTER_PARAM(TypeTag, bool, VtkWritePolymerConcentration, + "Include the concentration of the polymer component in the water phase " + "in the VTK output files"); + EWOMS_REGISTER_PARAM(TypeTag, bool, VtkWritePolymerDeadPoreVolume, + "Include the fraction of the \"dead\" pore volume " + "in the VTK output files"); + EWOMS_REGISTER_PARAM(TypeTag, bool, VtkWritePolymerRockDensity, + "Include the amount of already adsorbed polymer component" + "in the VTK output files"); + EWOMS_REGISTER_PARAM(TypeTag, bool, VtkWritePolymerAdsorption, + "Include the adsorption rate of the polymer component" + "in the VTK output files"); + EWOMS_REGISTER_PARAM(TypeTag, bool, VtkWritePolymerViscosityCorrection, + "Include the viscosity correction of the polymer component " + "in the VTK output files"); + EWOMS_REGISTER_PARAM(TypeTag, bool, VtkWriteWaterViscosityCorrection, + "Include the viscosity correction of the water component " + "due to polymers in the VTK output files"); + } + + /*! + * \brief Allocate memory for the scalar fields we would like to + * write to the VTK file. + */ + void allocBuffers() + { + if (!EWOMS_GET_PARAM(TypeTag, bool, EnableVtkOutput)) + return; + + if (!enablePolymer) + return; + + if (polymerConcentrationOutput_()) + this->resizeScalarBuffer_(polymerConcentration_); + if (polymerDeadPoreVolumeOutput_()) + this->resizeScalarBuffer_(polymerDeadPoreVolume_); + if (polymerRockDensityOutput_()) + this->resizeScalarBuffer_(polymerRockDensity_); + if (polymerAdsorptionOutput_()) + this->resizeScalarBuffer_(polymerAdsorption_); + if (polymerViscosityCorrectionOutput_()) + this->resizeScalarBuffer_(polymerViscosityCorrection_); + if (waterViscosityCorrectionOutput_()) + this->resizeScalarBuffer_(waterViscosityCorrection_); + } + + /*! + * \brief Modify the internal buffers according to the intensive quantities relevant for + * an element + */ + void processElement(const ElementContext& elemCtx) + { + if (!EWOMS_GET_PARAM(TypeTag, bool, EnableVtkOutput)) + return; + + if (!enablePolymer) + return; + + for (unsigned dofIdx = 0; dofIdx < elemCtx.numPrimaryDof(/*timeIdx=*/0); ++dofIdx) { + const auto& intQuants = elemCtx.intensiveQuantities(dofIdx, /*timeIdx=*/0); + unsigned globalDofIdx = elemCtx.globalSpaceIndex(dofIdx, /*timeIdx=*/0); + + if (polymerConcentrationOutput_()) + polymerConcentration_[globalDofIdx] = + Opm::scalarValue(intQuants.polymerConcentration()); + + if (polymerDeadPoreVolumeOutput_()) + polymerDeadPoreVolume_[globalDofIdx] = + Opm::scalarValue(intQuants.polymerDeadPoreVolume()); + + if (polymerRockDensityOutput_()) + polymerRockDensity_[globalDofIdx] = + Opm::scalarValue(intQuants.polymerRockDensity()); + + if (polymerAdsorptionOutput_()) + polymerAdsorption_[globalDofIdx] = + Opm::scalarValue(intQuants.polymerAdsorption()); + + if (polymerViscosityCorrectionOutput_()) + polymerViscosityCorrection_[globalDofIdx] = + Opm::scalarValue(intQuants.polymerViscosityCorrection()); + + if (waterViscosityCorrectionOutput_()) + waterViscosityCorrection_[globalDofIdx] = + Opm::scalarValue(intQuants.waterViscosityCorrection()); + } + } + + /*! + * \brief Add all buffers to the VTK output writer. + */ + void commitBuffers(BaseOutputWriter& baseWriter) + { + VtkMultiWriter *vtkWriter = dynamic_cast(&baseWriter); + if (!vtkWriter) + return; + + if (!enablePolymer) + return; + + if (polymerConcentrationOutput_()) + this->commitScalarBuffer_(baseWriter, "polymer concentration", polymerConcentration_); + + if (polymerDeadPoreVolumeOutput_()) + this->commitScalarBuffer_(baseWriter, "dead pore volume fraction", polymerDeadPoreVolume_); + + if (polymerRockDensityOutput_()) + this->commitScalarBuffer_(baseWriter, "polymer rock density", polymerRockDensity_); + + if (polymerAdsorptionOutput_()) + this->commitScalarBuffer_(baseWriter, "polymer adsorption", polymerAdsorption_); + + if (polymerViscosityCorrectionOutput_()) + this->commitScalarBuffer_(baseWriter, "polymer viscosity correction", polymerViscosityCorrection_); + + if (waterViscosityCorrectionOutput_()) + this->commitScalarBuffer_(baseWriter, "water viscosity correction", waterViscosityCorrection_); + } + +private: + static bool polymerConcentrationOutput_() + { + static bool val = EWOMS_GET_PARAM(TypeTag, bool, VtkWritePolymerConcentration); + return val; + } + + static bool polymerDeadPoreVolumeOutput_() + { + static bool val = EWOMS_GET_PARAM(TypeTag, bool, VtkWritePolymerDeadPoreVolume); + return val; + } + + static bool polymerRockDensityOutput_() + { + static bool val = EWOMS_GET_PARAM(TypeTag, bool, VtkWritePolymerRockDensity); + return val; + } + + static bool polymerAdsorptionOutput_() + { + static bool val = EWOMS_GET_PARAM(TypeTag, bool, VtkWritePolymerAdsorption); + return val; + } + + static bool polymerViscosityCorrectionOutput_() + { + static bool val = EWOMS_GET_PARAM(TypeTag, bool, VtkWritePolymerViscosityCorrection); + return val; + } + + static bool waterViscosityCorrectionOutput_() + { + static bool val = EWOMS_GET_PARAM(TypeTag, bool, VtkWritePolymerViscosityCorrection); + return val; + } + + ScalarBuffer polymerConcentration_; + ScalarBuffer polymerDeadPoreVolume_; + ScalarBuffer polymerRockDensity_; + ScalarBuffer polymerAdsorption_; + ScalarBuffer polymerViscosityCorrection_; + ScalarBuffer waterViscosityCorrection_; +}; +} // namespace Opm + +#endif diff --git a/opm/models/io/vtkblackoilsolventmodule.hh b/opm/models/io/vtkblackoilsolventmodule.hh new file mode 100644 index 000000000..cbdcd75fa --- /dev/null +++ b/opm/models/io/vtkblackoilsolventmodule.hh @@ -0,0 +1,230 @@ +// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- +// vi: set et ts=4 sw=4 sts=4: +/* + This file is part of the Open Porous Media project (OPM). + + OPM is free software: you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation, either version 2 of the License, or + (at your option) any later version. + + OPM is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU General Public License + along with OPM. If not, see . + + Consult the COPYING file in the top-level source directory of this + module for the precise wording of the license and the list of + copyright holders. +*/ +/*! + * \file + * \copydoc Opm::VtkBlackOilSolventModule + */ +#ifndef EWOMS_VTK_BLACK_OIL_SOLVENT_MODULE_HH +#define EWOMS_VTK_BLACK_OIL_SOLVENT_MODULE_HH + +#include + +#include "vtkmultiwriter.hh" +#include "baseoutputmodule.hh" + +#include +#include +#include + +#include + +#include + +BEGIN_PROPERTIES + +// create new type tag for the VTK multi-phase output +NEW_TYPE_TAG(VtkBlackOilSolvent); + +// create the property tags needed for the solvent output module +NEW_PROP_TAG(EnableSolvent); +NEW_PROP_TAG(EnableVtkOutput); +NEW_PROP_TAG(VtkWriteSolventSaturation); +NEW_PROP_TAG(VtkWriteSolventDensity); +NEW_PROP_TAG(VtkWriteSolventViscosity); +NEW_PROP_TAG(VtkWriteSolventMobility); + +// set default values for what quantities to output +SET_BOOL_PROP(VtkBlackOilSolvent, VtkWriteSolventSaturation, true); +SET_BOOL_PROP(VtkBlackOilSolvent, VtkWriteSolventDensity, true); +SET_BOOL_PROP(VtkBlackOilSolvent, VtkWriteSolventViscosity, true); +SET_BOOL_PROP(VtkBlackOilSolvent, VtkWriteSolventMobility, true); + +END_PROPERTIES + +namespace Opm { +/*! + * \ingroup Vtk + * + * \brief VTK output module for the black oil model's solvent related quantities. + */ +template +class VtkBlackOilSolventModule : public BaseOutputModule +{ + typedef BaseOutputModule ParentType; + + typedef typename GET_PROP_TYPE(TypeTag, Simulator) Simulator; + typedef typename GET_PROP_TYPE(TypeTag, GridView) GridView; + typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar; + typedef typename GET_PROP_TYPE(TypeTag, Evaluation) Evaluation; + typedef typename GET_PROP_TYPE(TypeTag, ElementContext) ElementContext; + + static const int vtkFormat = GET_PROP_VALUE(TypeTag, VtkOutputFormat); + typedef Opm::VtkMultiWriter VtkMultiWriter; + + enum { enableSolvent = GET_PROP_VALUE(TypeTag, EnableSolvent) }; + + typedef typename ParentType::ScalarBuffer ScalarBuffer; + +public: + VtkBlackOilSolventModule(const Simulator& simulator) + : ParentType(simulator) + { } + + /*! + * \brief Register all run-time parameters for the multi-phase VTK output + * module. + */ + static void registerParameters() + { + if (!enableSolvent) + return; + + EWOMS_REGISTER_PARAM(TypeTag, bool, VtkWriteSolventSaturation, + "Include the \"saturation\" of the solvent component " + "in the VTK output files"); + EWOMS_REGISTER_PARAM(TypeTag, bool, VtkWriteSolventDensity, + "Include the \"density\" of the solvent component " + "in the VTK output files"); + EWOMS_REGISTER_PARAM(TypeTag, bool, VtkWriteSolventViscosity, + "Include the \"viscosity\" of the solvent component " + "in the VTK output files"); + EWOMS_REGISTER_PARAM(TypeTag, bool, VtkWriteSolventMobility, + "Include the \"mobility\" of the solvent component " + "in the VTK output files"); + } + + /*! + * \brief Allocate memory for the scalar fields we would like to + * write to the VTK file. + */ + void allocBuffers() + { + if (!EWOMS_GET_PARAM(TypeTag, bool, EnableVtkOutput)) + return; + + if (!enableSolvent) + return; + + if (solventSaturationOutput_()) + this->resizeScalarBuffer_(solventSaturation_); + if (solventDensityOutput_()) + this->resizeScalarBuffer_(solventDensity_); + if (solventViscosityOutput_()) + this->resizeScalarBuffer_(solventViscosity_); + if (solventMobilityOutput_()) + this->resizeScalarBuffer_(solventMobility_); + } + + /*! + * \brief Modify the internal buffers according to the intensive quantities relevant for + * an element + */ + void processElement(const ElementContext& elemCtx) + { + if (!EWOMS_GET_PARAM(TypeTag, bool, EnableVtkOutput)) + return; + + if (!enableSolvent) + return; + + typedef Opm::MathToolbox Toolbox; + for (unsigned dofIdx = 0; dofIdx < elemCtx.numPrimaryDof(/*timeIdx=*/0); ++dofIdx) { + const auto& intQuants = elemCtx.intensiveQuantities(dofIdx, /*timeIdx=*/0); + unsigned globalDofIdx = elemCtx.globalSpaceIndex(dofIdx, /*timeIdx=*/0); + + if (solventSaturationOutput_()) + solventSaturation_[globalDofIdx] = + Toolbox::scalarValue(intQuants.solventSaturation()); + + if (solventDensityOutput_()) + solventDensity_[globalDofIdx] = + Toolbox::scalarValue(intQuants.solventDensity()); + + if (solventViscosityOutput_()) + solventViscosity_[globalDofIdx] = + Toolbox::scalarValue(intQuants.solventViscosity()); + + if (solventMobilityOutput_()) + solventMobility_[globalDofIdx] = + Toolbox::scalarValue(intQuants.solventMobility()); + } + } + + /*! + * \brief Add all buffers to the VTK output writer. + */ + void commitBuffers(BaseOutputWriter& baseWriter) + { + VtkMultiWriter *vtkWriter = dynamic_cast(&baseWriter); + if (!vtkWriter) + return; + + if (!enableSolvent) + return; + + if (solventSaturationOutput_()) + this->commitScalarBuffer_(baseWriter, "saturation_solvent", solventSaturation_); + + if (solventDensityOutput_()) + this->commitScalarBuffer_(baseWriter, "density_solvent", solventDensity_); + + if (solventViscosityOutput_()) + this->commitScalarBuffer_(baseWriter, "viscosity_solvent", solventViscosity_); + + if (solventMobilityOutput_()) + this->commitScalarBuffer_(baseWriter, "mobility_solvent", solventMobility_); + } + +private: + static bool solventSaturationOutput_() + { + static bool val = EWOMS_GET_PARAM(TypeTag, bool, VtkWriteSolventSaturation); + return val; + } + + static bool solventDensityOutput_() + { + static bool val = EWOMS_GET_PARAM(TypeTag, bool, VtkWriteSolventDensity); + return val; + } + + static bool solventViscosityOutput_() + { + static bool val = EWOMS_GET_PARAM(TypeTag, bool, VtkWriteSolventViscosity); + return val; + } + + static bool solventMobilityOutput_() + { + static bool val = EWOMS_GET_PARAM(TypeTag, bool, VtkWriteSolventMobility); + return val; + } + + ScalarBuffer solventSaturation_; + ScalarBuffer solventDensity_; + ScalarBuffer solventViscosity_; + ScalarBuffer solventMobility_; +}; +} // namespace Opm + +#endif diff --git a/opm/models/io/vtkcompositionmodule.hh b/opm/models/io/vtkcompositionmodule.hh new file mode 100644 index 000000000..172e53bdc --- /dev/null +++ b/opm/models/io/vtkcompositionmodule.hh @@ -0,0 +1,295 @@ +// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- +// vi: set et ts=4 sw=4 sts=4: +/* + This file is part of the Open Porous Media project (OPM). + + OPM is free software: you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation, either version 2 of the License, or + (at your option) any later version. + + OPM is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU General Public License + along with OPM. If not, see . + + Consult the COPYING file in the top-level source directory of this + module for the precise wording of the license and the list of + copyright holders. +*/ +/*! + * \file + * \copydoc Opm::VtkCompositionModule + */ +#ifndef EWOMS_VTK_COMPOSITION_MODULE_HH +#define EWOMS_VTK_COMPOSITION_MODULE_HH + +#include "vtkmultiwriter.hh" +#include "baseoutputmodule.hh" + +#include +#include + +#include + +BEGIN_PROPERTIES + +// create new type tag for the VTK composition output +NEW_TYPE_TAG(VtkComposition); + +// create the property tags needed for the composition module +NEW_PROP_TAG(VtkWriteMassFractions); +NEW_PROP_TAG(VtkWriteMoleFractions); +NEW_PROP_TAG(VtkWriteTotalMassFractions); +NEW_PROP_TAG(VtkWriteTotalMoleFractions); +NEW_PROP_TAG(VtkWriteMolarities); +NEW_PROP_TAG(VtkWriteFugacities); +NEW_PROP_TAG(VtkWriteFugacityCoeffs); +NEW_PROP_TAG(VtkOutputFormat); +NEW_PROP_TAG(EnableVtkOutput); + +// set default values for what quantities to output +SET_BOOL_PROP(VtkComposition, VtkWriteMassFractions, false); +SET_BOOL_PROP(VtkComposition, VtkWriteMoleFractions, true); +SET_BOOL_PROP(VtkComposition, VtkWriteTotalMassFractions, false); +SET_BOOL_PROP(VtkComposition, VtkWriteTotalMoleFractions, false); +SET_BOOL_PROP(VtkComposition, VtkWriteMolarities, false); +SET_BOOL_PROP(VtkComposition, VtkWriteFugacities, false); +SET_BOOL_PROP(VtkComposition, VtkWriteFugacityCoeffs, false); + +END_PROPERTIES + +namespace Opm { + +/*! + * \ingroup Vtk + * + * \brief VTK output module for the fluid composition + * + * This module deals with the following quantities: + * - Mole fraction of a component in a fluid phase + * - Mass fraction of a component in a fluid phase + * - Molarity (i.e. molar concentration) of a component in a fluid phase + * - Fugacity of all components + * - FugacityCoefficient of all components in all phases + */ +template +class VtkCompositionModule : public BaseOutputModule +{ + typedef BaseOutputModule ParentType; + + typedef typename GET_PROP_TYPE(TypeTag, Simulator) Simulator; + typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar; + typedef typename GET_PROP_TYPE(TypeTag, Evaluation) Evaluation; + typedef typename GET_PROP_TYPE(TypeTag, ElementContext) ElementContext; + + typedef typename GET_PROP_TYPE(TypeTag, GridView) GridView; + + enum { numPhases = GET_PROP_VALUE(TypeTag, NumPhases) }; + enum { numComponents = GET_PROP_VALUE(TypeTag, NumComponents) }; + + static const int vtkFormat = GET_PROP_VALUE(TypeTag, VtkOutputFormat); + typedef Opm::VtkMultiWriter VtkMultiWriter; + + typedef typename ParentType::ComponentBuffer ComponentBuffer; + typedef typename ParentType::PhaseComponentBuffer PhaseComponentBuffer; + +public: + VtkCompositionModule(const Simulator& simulator) + : ParentType(simulator) + { } + + /*! + * \brief Register all run-time parameters for the Vtk output module. + */ + static void registerParameters() + { + EWOMS_REGISTER_PARAM(TypeTag, bool, VtkWriteMassFractions, + "Include mass fractions in the VTK output files"); + EWOMS_REGISTER_PARAM(TypeTag, bool, VtkWriteMoleFractions, + "Include mole fractions in the VTK output files"); + EWOMS_REGISTER_PARAM(TypeTag, bool, VtkWriteTotalMassFractions, + "Include total mass fractions in the VTK output files"); + EWOMS_REGISTER_PARAM(TypeTag, bool, VtkWriteTotalMoleFractions, + "Include total mole fractions in the VTK output files"); + EWOMS_REGISTER_PARAM(TypeTag, bool, VtkWriteMolarities, + "Include component molarities in the VTK output files"); + EWOMS_REGISTER_PARAM(TypeTag, bool, VtkWriteFugacities, + "Include component fugacities in the VTK output files"); + EWOMS_REGISTER_PARAM(TypeTag, bool, VtkWriteFugacityCoeffs, + "Include component fugacity coefficients in the VTK output files"); + } + + /*! + * \brief Allocate memory for the scalar fields we would like to + * write to the VTK file. + */ + void allocBuffers() + { + if (moleFracOutput_()) + this->resizePhaseComponentBuffer_(moleFrac_); + if (massFracOutput_()) + this->resizePhaseComponentBuffer_(massFrac_); + if (totalMassFracOutput_()) + this->resizeComponentBuffer_(totalMassFrac_); + if (totalMoleFracOutput_()) + this->resizeComponentBuffer_(totalMoleFrac_); + if (molarityOutput_()) + this->resizePhaseComponentBuffer_(molarity_); + + if (fugacityOutput_()) + this->resizeComponentBuffer_(fugacity_); + if (fugacityCoeffOutput_()) + this->resizePhaseComponentBuffer_(fugacityCoeff_); + } + + /*! + * \brief Modify the internal buffers according to the intensive quantities relevant + * for an element + */ + void processElement(const ElementContext& elemCtx) + { + typedef Opm::MathToolbox Toolbox; + + if (!EWOMS_GET_PARAM(TypeTag, bool, EnableVtkOutput)) + return; + + for (unsigned i = 0; i < elemCtx.numPrimaryDof(/*timeIdx=*/0); ++i) { + unsigned I = elemCtx.globalSpaceIndex(i, /*timeIdx=*/0); + const auto& intQuants = elemCtx.intensiveQuantities(i, /*timeIdx=*/0); + const auto& fs = intQuants.fluidState(); + + for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) { + for (unsigned compIdx = 0; compIdx < numComponents; ++compIdx) { + if (moleFracOutput_()) + moleFrac_[phaseIdx][compIdx][I] = Toolbox::value(fs.moleFraction(phaseIdx, compIdx)); + if (massFracOutput_()) + massFrac_[phaseIdx][compIdx][I] = Toolbox::value(fs.massFraction(phaseIdx, compIdx)); + if (molarityOutput_()) + molarity_[phaseIdx][compIdx][I] = Toolbox::value(fs.molarity(phaseIdx, compIdx)); + + if (fugacityCoeffOutput_()) + fugacityCoeff_[phaseIdx][compIdx][I] = + Toolbox::value(fs.fugacityCoefficient(phaseIdx, compIdx)); + } + } + + for (unsigned compIdx = 0; compIdx < numComponents; ++compIdx) { + if (totalMassFracOutput_()) { + Scalar compMass = 0; + Scalar totalMass = 0; + for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) { + totalMass += Toolbox::value(fs.density(phaseIdx)) * Toolbox::value(fs.saturation(phaseIdx)); + compMass += + Toolbox::value(fs.density(phaseIdx)) + *Toolbox::value(fs.saturation(phaseIdx)) + *Toolbox::value(fs.massFraction(phaseIdx, compIdx)); + } + totalMassFrac_[compIdx][I] = compMass / totalMass; + } + if (totalMoleFracOutput_()) { + Scalar compMoles = 0; + Scalar totalMoles = 0; + for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) { + totalMoles += + Toolbox::value(fs.molarDensity(phaseIdx)) + *Toolbox::value(fs.saturation(phaseIdx)); + compMoles += + Toolbox::value(fs.molarDensity(phaseIdx)) + *Toolbox::value(fs.saturation(phaseIdx)) + *Toolbox::value(fs.moleFraction(phaseIdx, compIdx)); + } + totalMoleFrac_[compIdx][I] = compMoles / totalMoles; + } + if (fugacityOutput_()) + fugacity_[compIdx][I] = Toolbox::value(intQuants.fluidState().fugacity(/*phaseIdx=*/0, compIdx)); + } + } + } + + /*! + * \brief Add all buffers to the VTK output writer. + */ + void commitBuffers(BaseOutputWriter& baseWriter) + { + VtkMultiWriter *vtkWriter = dynamic_cast(&baseWriter); + if (!vtkWriter) { + return; + } + + if (moleFracOutput_()) + this->commitPhaseComponentBuffer_(baseWriter, "moleFrac_%s^%s", moleFrac_); + if (massFracOutput_()) + this->commitPhaseComponentBuffer_(baseWriter, "massFrac_%s^%s", massFrac_); + if (molarityOutput_()) + this->commitPhaseComponentBuffer_(baseWriter, "molarity_%s^%s", molarity_); + if (totalMassFracOutput_()) + this->commitComponentBuffer_(baseWriter, "totalMassFrac^%s", totalMassFrac_); + if (totalMoleFracOutput_()) + this->commitComponentBuffer_(baseWriter, "totalMoleFrac^%s", totalMoleFrac_); + + if (fugacityOutput_()) + this->commitComponentBuffer_(baseWriter, "fugacity^%s", fugacity_); + if (fugacityCoeffOutput_()) + this->commitPhaseComponentBuffer_(baseWriter, "fugacityCoeff_%s^%s", fugacityCoeff_); + } + +private: + static bool massFracOutput_() + { + static bool val = EWOMS_GET_PARAM(TypeTag, bool, VtkWriteMassFractions); + return val; + } + + static bool moleFracOutput_() + { + static bool val = EWOMS_GET_PARAM(TypeTag, bool, VtkWriteMoleFractions); + return val; + } + + static bool totalMassFracOutput_() + { + static bool val = EWOMS_GET_PARAM(TypeTag, bool, VtkWriteTotalMassFractions); + return val; + } + + static bool totalMoleFracOutput_() + { + static bool val = EWOMS_GET_PARAM(TypeTag, bool, VtkWriteTotalMoleFractions); + return val; + } + + static bool molarityOutput_() + { + static bool val = EWOMS_GET_PARAM(TypeTag, bool, VtkWriteMolarities); + return val; + } + + static bool fugacityOutput_() + { + static bool val = EWOMS_GET_PARAM(TypeTag, bool, VtkWriteFugacities); + return val; + } + + static bool fugacityCoeffOutput_() + { + static bool val = EWOMS_GET_PARAM(TypeTag, bool, VtkWriteFugacityCoeffs); + return val; + } + + PhaseComponentBuffer moleFrac_; + PhaseComponentBuffer massFrac_; + PhaseComponentBuffer molarity_; + ComponentBuffer totalMassFrac_; + ComponentBuffer totalMoleFrac_; + + ComponentBuffer fugacity_; + PhaseComponentBuffer fugacityCoeff_; +}; + +} // namespace Opm + +#endif diff --git a/opm/models/io/vtkdiffusionmodule.hh b/opm/models/io/vtkdiffusionmodule.hh new file mode 100644 index 000000000..ac9d75c7c --- /dev/null +++ b/opm/models/io/vtkdiffusionmodule.hh @@ -0,0 +1,204 @@ +// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- +// vi: set et ts=4 sw=4 sts=4: +/* + This file is part of the Open Porous Media project (OPM). + + OPM is free software: you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation, either version 2 of the License, or + (at your option) any later version. + + OPM is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU General Public License + along with OPM. If not, see . + + Consult the COPYING file in the top-level source directory of this + module for the precise wording of the license and the list of + copyright holders. +*/ +/*! + * \file + * + * \copydoc Opm::VtkDiffusionModule + */ +#ifndef EWOMS_VTK_DIFFUSION_MODULE_HH +#define EWOMS_VTK_DIFFUSION_MODULE_HH + +#include "vtkmultiwriter.hh" +#include "baseoutputmodule.hh" + +#include +#include + +#include +#include + +BEGIN_PROPERTIES + +// create new type tag for the VTK output of the quantities for molecular +// diffusion +NEW_TYPE_TAG(VtkDiffusion); + +// create the property tags needed for the diffusion module +NEW_PROP_TAG(VtkWriteTortuosities); +NEW_PROP_TAG(VtkWriteDiffusionCoefficients); +NEW_PROP_TAG(VtkWriteEffectiveDiffusionCoefficients); +NEW_PROP_TAG(VtkOutputFormat); +NEW_PROP_TAG(EnableVtkOutput); + +// set default values for what quantities to output +SET_BOOL_PROP(VtkDiffusion, VtkWriteTortuosities, false); +SET_BOOL_PROP(VtkDiffusion, VtkWriteDiffusionCoefficients, false); +SET_BOOL_PROP(VtkDiffusion, VtkWriteEffectiveDiffusionCoefficients, false); + +END_PROPERTIES + +namespace Opm { +/*! + * \ingroup Vtk + * + * \brief VTK output module for quantities which make sense for models which + * incorperate molecular diffusion. + * + * This module deals with the following quantities: + * - Molecular diffusion coefficients of all components in all fluid phases + * - Effective molecular diffusion coefficients of the porous medium of all + * components in all fluid phases + */ +template +class VtkDiffusionModule : public BaseOutputModule +{ + typedef BaseOutputModule ParentType; + + typedef typename GET_PROP_TYPE(TypeTag, Simulator) Simulator; + typedef typename GET_PROP_TYPE(TypeTag, ElementContext) ElementContext; + typedef typename GET_PROP_TYPE(TypeTag, GridView) GridView; + typedef typename GET_PROP_TYPE(TypeTag, Evaluation) Evaluation; + + typedef Opm::MathToolbox Toolbox; + + typedef typename ParentType::PhaseComponentBuffer PhaseComponentBuffer; + typedef typename ParentType::PhaseBuffer PhaseBuffer; + + static const int vtkFormat = GET_PROP_VALUE(TypeTag, VtkOutputFormat); + typedef Opm::VtkMultiWriter VtkMultiWriter; + + enum { numPhases = GET_PROP_VALUE(TypeTag, NumPhases) }; + enum { numComponents = GET_PROP_VALUE(TypeTag, NumComponents) }; + +public: + VtkDiffusionModule(const Simulator& simulator) + : ParentType(simulator) + { } + + /*! + * \brief Register all run-time parameters for the Vtk output module. + */ + static void registerParameters() + { + EWOMS_REGISTER_PARAM(TypeTag, bool, VtkWriteTortuosities, + "Include the tortuosity for each phase in the VTK " + "output files"); + EWOMS_REGISTER_PARAM(TypeTag, bool, VtkWriteDiffusionCoefficients, + "Include the molecular diffusion coefficients in " + "the VTK output files"); + EWOMS_REGISTER_PARAM(TypeTag, bool, + VtkWriteEffectiveDiffusionCoefficients, + "Include the effective molecular diffusion " + "coefficients the medium in the VTK output files"); + } + + /*! + * \brief Allocate memory for the scalar fields we would like to + * write to the VTK file. + */ + void allocBuffers() + { + if (tortuosityOutput_()) + this->resizePhaseBuffer_(tortuosity_); + if (diffusionCoefficientOutput_()) + this->resizePhaseComponentBuffer_(diffusionCoefficient_); + if (effectiveDiffusionCoefficientOutput_()) + this->resizePhaseComponentBuffer_(effectiveDiffusionCoefficient_); + } + + /*! + * \brief Modify the internal buffers according to the intensive quanties relevant + * for an element + */ + void processElement(const ElementContext& elemCtx) + { + if (!EWOMS_GET_PARAM(TypeTag, bool, EnableVtkOutput)) + return; + + for (unsigned i = 0; i < elemCtx.numPrimaryDof(/*timeIdx=*/0); ++i) { + unsigned I = elemCtx.globalSpaceIndex(i, /*timeIdx=*/0); + const auto& intQuants = elemCtx.intensiveQuantities(i, /*timeIdx=*/0); + + for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) { + if (tortuosityOutput_()) + tortuosity_[phaseIdx][I] = Toolbox::value(intQuants.tortuosity(phaseIdx)); + for (unsigned compIdx = 0; compIdx < numComponents; ++compIdx) { + if (diffusionCoefficientOutput_()) + diffusionCoefficient_[phaseIdx][compIdx][I] = + Toolbox::value(intQuants.diffusionCoefficient(phaseIdx, compIdx)); + if (effectiveDiffusionCoefficientOutput_()) + effectiveDiffusionCoefficient_[phaseIdx][compIdx][I] = + Toolbox::value(intQuants.effectiveDiffusionCoefficient(phaseIdx, compIdx)); + } + } + } + } + + /*! + * \brief Add all buffers to the VTK output writer. + */ + void commitBuffers(BaseOutputWriter& baseWriter) + { + VtkMultiWriter *vtkWriter = dynamic_cast(&baseWriter); + if (!vtkWriter) { + return; + } + + if (tortuosityOutput_()) + this->commitPhaseBuffer_(baseWriter, "tortuosity", tortuosity_); + if (diffusionCoefficientOutput_()) + this->commitPhaseComponentBuffer_(baseWriter, "diffusionCoefficient", + diffusionCoefficient_); + if (effectiveDiffusionCoefficientOutput_()) + this->commitPhaseComponentBuffer_(baseWriter, + "effectiveDiffusionCoefficient", + effectiveDiffusionCoefficient_); + } + +private: + static bool tortuosityOutput_() + { + static bool val = EWOMS_GET_PARAM(TypeTag, bool, VtkWriteTortuosities); + return val; + } + + static bool diffusionCoefficientOutput_() + { + static bool val = EWOMS_GET_PARAM(TypeTag, bool, VtkWriteDiffusionCoefficients); + return val; + } + + static bool effectiveDiffusionCoefficientOutput_() + { + static bool val = EWOMS_GET_PARAM(TypeTag, bool, VtkWriteEffectiveDiffusionCoefficients); + return val; + } + + PhaseBuffer tortuosity_; + PhaseComponentBuffer diffusionCoefficient_; + PhaseComponentBuffer effectiveDiffusionCoefficient_; +}; + +} // namespace Opm + +#endif diff --git a/opm/models/io/vtkdiscretefracturemodule.hh b/opm/models/io/vtkdiscretefracturemodule.hh new file mode 100644 index 000000000..bf4dc6188 --- /dev/null +++ b/opm/models/io/vtkdiscretefracturemodule.hh @@ -0,0 +1,364 @@ +// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- +// vi: set et ts=4 sw=4 sts=4: +/* + This file is part of the Open Porous Media project (OPM). + + OPM is free software: you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation, either version 2 of the License, or + (at your option) any later version. + + OPM is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU General Public License + along with OPM. If not, see . + + Consult the COPYING file in the top-level source directory of this + module for the precise wording of the license and the list of + copyright holders. +*/ +/*! + * \file + * \copydoc Opm::VtkDiscreteFractureModule + */ +#ifndef EWOMS_VTK_DISCRETE_FRACTURE_MODULE_HH +#define EWOMS_VTK_DISCRETE_FRACTURE_MODULE_HH + +#include "vtkmultiwriter.hh" +#include "baseoutputmodule.hh" + +#include +#include + +#include + +#include + +#include + +BEGIN_PROPERTIES + +// create new type tag for the VTK multi-phase output +NEW_TYPE_TAG(VtkDiscreteFracture); + +// create the property tags needed for the multi phase module +NEW_PROP_TAG(Vanguard); +NEW_PROP_TAG(VtkWriteFractureSaturations); +NEW_PROP_TAG(VtkWriteFractureMobilities); +NEW_PROP_TAG(VtkWriteFractureRelativePermeabilities); +NEW_PROP_TAG(VtkWriteFracturePorosity); +NEW_PROP_TAG(VtkWriteFractureIntrinsicPermeabilities); +NEW_PROP_TAG(VtkWriteFractureFilterVelocities); +NEW_PROP_TAG(VtkWriteFractureVolumeFraction); +NEW_PROP_TAG(VtkOutputFormat); +NEW_PROP_TAG(EnableVtkOutput); +NEW_PROP_TAG(DiscBaseOutputModule); + +// set default values for what quantities to output +SET_BOOL_PROP(VtkDiscreteFracture, VtkWriteFractureSaturations, true); +SET_BOOL_PROP(VtkDiscreteFracture, VtkWriteFractureMobilities, false); +SET_BOOL_PROP(VtkDiscreteFracture, VtkWriteFractureRelativePermeabilities, true); +SET_BOOL_PROP(VtkDiscreteFracture, VtkWriteFracturePorosity, true); +SET_BOOL_PROP(VtkDiscreteFracture, VtkWriteFractureIntrinsicPermeabilities, false); +SET_BOOL_PROP(VtkDiscreteFracture, VtkWriteFractureFilterVelocities, false); +SET_BOOL_PROP(VtkDiscreteFracture, VtkWriteFractureVolumeFraction, true); + +END_PROPERTIES + +namespace Opm { +/*! + * \ingroup Vtk + * + * \brief VTK output module for quantities which make sense for all + * models which deal with discrete fractures in porous media. + * + * This module deals with the following quantities: + * - Saturations of all fluid phases in the fracture + * - Mobilities of all fluid phases in the fracture + * - Relative permeabilities of all fluid phases in the fracture + * - Porosity of the medium in the fracture + * - Norm of the intrinsic permeability of the medium in the fracture + */ +template +class VtkDiscreteFractureModule : public BaseOutputModule +{ + typedef BaseOutputModule ParentType; + + typedef typename GET_PROP_TYPE(TypeTag, Simulator) Simulator; + typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar; + typedef typename GET_PROP_TYPE(TypeTag, ElementContext) ElementContext; + + typedef typename GET_PROP_TYPE(TypeTag, Vanguard) Vanguard; + typedef typename GET_PROP_TYPE(TypeTag, GridView) GridView; + typedef typename GET_PROP_TYPE(TypeTag, FluidSystem) FluidSystem; + + typedef typename GET_PROP_TYPE(TypeTag, DiscBaseOutputModule) DiscBaseOutputModule; + + static const int vtkFormat = GET_PROP_VALUE(TypeTag, VtkOutputFormat); + typedef Opm::VtkMultiWriter VtkMultiWriter; + + enum { dim = GridView::dimension }; + enum { dimWorld = GridView::dimensionworld }; + enum { numPhases = GET_PROP_VALUE(TypeTag, NumPhases) }; + + typedef typename ParentType::ScalarBuffer ScalarBuffer; + typedef typename ParentType::PhaseBuffer PhaseBuffer; + typedef typename ParentType::PhaseVectorBuffer PhaseVectorBuffer; + +public: + VtkDiscreteFractureModule(const Simulator& simulator) + : ParentType(simulator) + { } + + /*! + * \brief Register all run-time parameters for the multi-phase VTK output module. + */ + static void registerParameters() + { + EWOMS_REGISTER_PARAM(TypeTag, bool, VtkWriteFractureSaturations, + "Include the phase saturations in the VTK output files"); + EWOMS_REGISTER_PARAM(TypeTag, bool, VtkWriteFractureMobilities, + "Include the phase mobilities in the VTK output files"); + EWOMS_REGISTER_PARAM(TypeTag, bool, VtkWriteFractureRelativePermeabilities, + "Include the phase relative permeabilities in the " + "VTK output files"); + EWOMS_REGISTER_PARAM(TypeTag, bool, VtkWriteFracturePorosity, + "Include the porosity in the VTK output files"); + EWOMS_REGISTER_PARAM(TypeTag, bool, VtkWriteFractureIntrinsicPermeabilities, + "Include the intrinsic permeability in the VTK output files"); + EWOMS_REGISTER_PARAM(TypeTag, bool, VtkWriteFractureFilterVelocities, + "Include in the filter velocities of the phases " + "the VTK output files"); + EWOMS_REGISTER_PARAM(TypeTag, bool, VtkWriteFractureVolumeFraction, + "Add the fraction of the total volume which is " + "occupied by fractures in the VTK output"); + } + + /*! + * \brief Allocate memory for the scalar fields we would like to + * write to the VTK file. + */ + void allocBuffers() + { + if (saturationOutput_()) + this->resizePhaseBuffer_(fractureSaturation_); + if (mobilityOutput_()) + this->resizePhaseBuffer_(fractureMobility_); + if (relativePermeabilityOutput_()) + this->resizePhaseBuffer_(fractureRelativePermeability_); + + if (porosityOutput_()) + this->resizeScalarBuffer_(fracturePorosity_); + if (intrinsicPermeabilityOutput_()) + this->resizeScalarBuffer_(fractureIntrinsicPermeability_); + if (volumeFractionOutput_()) + this->resizeScalarBuffer_(fractureVolumeFraction_); + + if (velocityOutput_()) { + size_t nDof = this->simulator_.model().numGridDof(); + for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) { + fractureVelocity_[phaseIdx].resize(nDof); + for (unsigned dofIdx = 0; dofIdx < nDof; ++dofIdx) { + fractureVelocity_[phaseIdx][dofIdx].resize(dimWorld); + fractureVelocity_[phaseIdx][dofIdx] = 0.0; + } + } + this->resizePhaseBuffer_(fractureVelocityWeight_); + } + } + + /*! + * \brief Modify the internal buffers according to the intensive quantities relevant for + * an element + */ + void processElement(const ElementContext& elemCtx) + { + if (!EWOMS_GET_PARAM(TypeTag, bool, EnableVtkOutput)) + return; + + const auto& fractureMapper = elemCtx.simulator().vanguard().fractureMapper(); + + for (unsigned i = 0; i < elemCtx.numPrimaryDof(/*timeIdx=*/0); ++i) { + unsigned I = elemCtx.globalSpaceIndex(i, /*timeIdx=*/0); + if (!fractureMapper.isFractureVertex(I)) + continue; + + const auto& intQuants = elemCtx.intensiveQuantities(i, /*timeIdx=*/0); + const auto& fs = intQuants.fractureFluidState(); + + if (porosityOutput_()) { + Opm::Valgrind::CheckDefined(intQuants.fracturePorosity()); + fracturePorosity_[I] = intQuants.fracturePorosity(); + } + if (intrinsicPermeabilityOutput_()) { + const auto& K = intQuants.fractureIntrinsicPermeability(); + fractureIntrinsicPermeability_[I] = K[0][0]; + } + + for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) { + if (saturationOutput_()) { + Opm::Valgrind::CheckDefined(fs.saturation(phaseIdx)); + fractureSaturation_[phaseIdx][I] = fs.saturation(phaseIdx); + } + if (mobilityOutput_()) { + Opm::Valgrind::CheckDefined(intQuants.fractureMobility(phaseIdx)); + fractureMobility_[phaseIdx][I] = intQuants.fractureMobility(phaseIdx); + } + if (relativePermeabilityOutput_()) { + Opm::Valgrind::CheckDefined(intQuants.fractureRelativePermeability(phaseIdx)); + fractureRelativePermeability_[phaseIdx][I] = + intQuants.fractureRelativePermeability(phaseIdx); + } + if (volumeFractionOutput_()) { + Opm::Valgrind::CheckDefined(intQuants.fractureVolume()); + fractureVolumeFraction_[I] += intQuants.fractureVolume(); + } + } + } + + if (velocityOutput_()) { + // calculate velocities if requested by the simulator + for (unsigned scvfIdx = 0; scvfIdx < elemCtx.numInteriorFaces(/*timeIdx=*/0); ++ scvfIdx) { + const auto& extQuants = elemCtx.extensiveQuantities(scvfIdx, /*timeIdx=*/0); + + unsigned i = extQuants.interiorIndex(); + unsigned I = elemCtx.globalSpaceIndex(i, /*timeIdx=*/0); + + unsigned j = extQuants.exteriorIndex(); + unsigned J = elemCtx.globalSpaceIndex(j, /*timeIdx=*/0); + + if (!fractureMapper.isFractureEdge(I, J)) + continue; + + for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) { + Scalar weight = + std::max(1e-16, std::abs(extQuants.fractureVolumeFlux(phaseIdx))); + Opm::Valgrind::CheckDefined(extQuants.extrusionFactor()); + assert(extQuants.extrusionFactor() > 0); + weight *= extQuants.extrusionFactor(); + + Dune::FieldVector v(extQuants.fractureFilterVelocity(phaseIdx)); + v *= weight; + + for (unsigned dimIdx = 0; dimIdx < dimWorld; ++dimIdx) { + fractureVelocity_[phaseIdx][I][dimIdx] += v[dimIdx]; + fractureVelocity_[phaseIdx][J][dimIdx] += v[dimIdx]; + } + + fractureVelocityWeight_[phaseIdx][I] += weight; + fractureVelocityWeight_[phaseIdx][J] += weight; + } + } + } + } + + /*! + * \brief Add all buffers to the VTK output writer. + */ + void commitBuffers(BaseOutputWriter& baseWriter) + { + VtkMultiWriter *vtkWriter = dynamic_cast(&baseWriter); + if (!vtkWriter) { + return; + } + + if (saturationOutput_()) + this->commitPhaseBuffer_(baseWriter, "fractureSaturation_%s", fractureSaturation_); + if (mobilityOutput_()) + this->commitPhaseBuffer_(baseWriter, "fractureMobility_%s", fractureMobility_); + if (relativePermeabilityOutput_()) + this->commitPhaseBuffer_(baseWriter, "fractureRelativePerm_%s", fractureRelativePermeability_); + + if (porosityOutput_()) + this->commitScalarBuffer_(baseWriter, "fracturePorosity", fracturePorosity_); + if (intrinsicPermeabilityOutput_()) + this->commitScalarBuffer_(baseWriter, "fractureIntrinsicPerm", fractureIntrinsicPermeability_); + if (volumeFractionOutput_()) { + // divide the fracture volume by the total volume of the finite volumes + for (unsigned I = 0; I < fractureVolumeFraction_.size(); ++I) + fractureVolumeFraction_[I] /= this->simulator_.model().dofTotalVolume(I); + this->commitScalarBuffer_(baseWriter, "fractureVolumeFraction", fractureVolumeFraction_); + } + + if (velocityOutput_()) { + size_t nDof = this->simulator_.model().numGridDof(); + + for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) { + // first, divide the velocity field by the + // respective finite volume's surface area + for (unsigned dofIdx = 0; dofIdx < nDof; ++dofIdx) + fractureVelocity_[phaseIdx][dofIdx] /= + std::max(1e-20, fractureVelocityWeight_[phaseIdx][dofIdx]); + // commit the phase velocity + char name[512]; + snprintf(name, 512, "fractureFilterVelocity_%s", FluidSystem::phaseName(phaseIdx)); + + DiscBaseOutputModule::attachVectorDofData_(baseWriter, fractureVelocity_[phaseIdx], name); + } + } + } + +private: + static bool saturationOutput_() + { + static bool val = EWOMS_GET_PARAM(TypeTag, bool, VtkWriteFractureSaturations); + return val; + } + + static bool mobilityOutput_() + { + static bool val = EWOMS_GET_PARAM(TypeTag, bool, VtkWriteFractureMobilities); + return val; + } + + static bool relativePermeabilityOutput_() + { + static bool val = EWOMS_GET_PARAM(TypeTag, bool, VtkWriteFractureRelativePermeabilities); + return val; + } + + static bool porosityOutput_() + { + static bool val = EWOMS_GET_PARAM(TypeTag, bool, VtkWriteFracturePorosity); + return val; + } + + static bool intrinsicPermeabilityOutput_() + { + static bool val = EWOMS_GET_PARAM(TypeTag, bool, VtkWriteFractureIntrinsicPermeabilities); + return val; + } + + static bool volumeFractionOutput_() + { + static bool val = EWOMS_GET_PARAM(TypeTag, bool, VtkWriteFractureVolumeFraction); + return val; + } + + static bool velocityOutput_() + { + static bool val = EWOMS_GET_PARAM(TypeTag, bool, VtkWriteFractureFilterVelocities); + return val; + } + + PhaseBuffer fractureSaturation_; + PhaseBuffer fractureMobility_; + PhaseBuffer fractureRelativePermeability_; + + ScalarBuffer fracturePorosity_; + ScalarBuffer fractureVolumeFraction_; + ScalarBuffer fractureIntrinsicPermeability_; + + PhaseVectorBuffer fractureVelocity_; + PhaseBuffer fractureVelocityWeight_; + + PhaseVectorBuffer potentialGradient_; + PhaseBuffer potentialWeight_; +}; + +} // namespace Opm + +#endif diff --git a/opm/models/io/vtkenergymodule.hh b/opm/models/io/vtkenergymodule.hh new file mode 100644 index 000000000..f0e6323c0 --- /dev/null +++ b/opm/models/io/vtkenergymodule.hh @@ -0,0 +1,218 @@ +// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- +// vi: set et ts=4 sw=4 sts=4: +/* + This file is part of the Open Porous Media project (OPM). + + OPM is free software: you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation, either version 2 of the License, or + (at your option) any later version. + + OPM is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU General Public License + along with OPM. If not, see . + + Consult the COPYING file in the top-level source directory of this + module for the precise wording of the license and the list of + copyright holders. +*/ +/*! + * \file + * \copydoc Opm::VtkEnergyModule + */ +#ifndef EWOMS_VTK_ENERGY_MODULE_HH +#define EWOMS_VTK_ENERGY_MODULE_HH + +#include "vtkmultiwriter.hh" +#include "baseoutputmodule.hh" + +#include +#include + +#include + +BEGIN_PROPERTIES + +// create new type tag for the VTK energy output +NEW_TYPE_TAG(VtkEnergy); + +// create the property tags needed for the energy module +NEW_PROP_TAG(VtkWriteSolidInternalEnergy); +NEW_PROP_TAG(VtkWriteThermalConductivity); +NEW_PROP_TAG(VtkWriteInternalEnergies); +NEW_PROP_TAG(VtkWriteEnthalpies); +NEW_PROP_TAG(VtkOutputFormat); +NEW_PROP_TAG(EnableVtkOutput); + +// set default values for what quantities to output +SET_BOOL_PROP(VtkEnergy, VtkWriteSolidInternalEnergy, false); +SET_BOOL_PROP(VtkEnergy, VtkWriteThermalConductivity, false); +SET_BOOL_PROP(VtkEnergy, VtkWriteInternalEnergies, false); +SET_BOOL_PROP(VtkEnergy, VtkWriteEnthalpies, false); + +END_PROPERTIES + +namespace Opm { +/*! + * \ingroup Vtk + * + * \brief VTK output module for quantities which make sense for models which + * assume thermal equilibrium. + * + * This module deals with the following quantities: + * - Specific enthalpy of all fluid phases + * - Specific internal energy of all fluid phases + * - Volumetric internal energy of the solid phase + * - Total thermal conductivity, i.e. the conductivity of the solid and all fluid phases + * combined + */ +template +class VtkEnergyModule : public BaseOutputModule +{ + typedef BaseOutputModule ParentType; + + typedef typename GET_PROP_TYPE(TypeTag, Simulator) Simulator; + typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar; + typedef typename GET_PROP_TYPE(TypeTag, Evaluation) Evaluation; + typedef typename GET_PROP_TYPE(TypeTag, ElementContext) ElementContext; + typedef typename GET_PROP_TYPE(TypeTag, GridView) GridView; + + typedef typename ParentType::ScalarBuffer ScalarBuffer; + typedef typename ParentType::PhaseBuffer PhaseBuffer; + + static const int vtkFormat = GET_PROP_VALUE(TypeTag, VtkOutputFormat); + enum { numPhases = GET_PROP_VALUE(TypeTag, NumPhases) }; + + typedef typename Opm::MathToolbox Toolbox; + typedef Opm::VtkMultiWriter VtkMultiWriter; + +public: + VtkEnergyModule(const Simulator& simulator) + : ParentType(simulator) + { + } + + /*! + * \brief Register all run-time parameters for the Vtk output module. + */ + static void registerParameters() + { + EWOMS_REGISTER_PARAM(TypeTag, bool, VtkWriteSolidInternalEnergy, + "Include the volumetric internal energy of solid" + "matrix in the VTK output files"); + EWOMS_REGISTER_PARAM(TypeTag, bool, VtkWriteThermalConductivity, + "Include the total thermal conductivity of the" + "medium in the VTK output files"); + EWOMS_REGISTER_PARAM(TypeTag, bool, VtkWriteEnthalpies, + "Include the specific enthalpy of the phases in " + "the VTK output files"); + EWOMS_REGISTER_PARAM(TypeTag, bool, VtkWriteInternalEnergies, + "Include the specific internal energy of the " + "phases in the VTK output files"); + } + + /*! + * \brief Allocate memory for the scalar fields we would like to + * write to the VTK file. + */ + void allocBuffers() + { + if (enthalpyOutput_()) + this->resizePhaseBuffer_(enthalpy_); + if (internalEnergyOutput_()) + this->resizePhaseBuffer_(internalEnergy_); + + if (solidInternalEnergyOutput_()) + this->resizeScalarBuffer_(solidInternalEnergy_); + if (thermalConductivityOutput_()) + this->resizeScalarBuffer_(thermalConductivity_); + } + + /*! + * \brief Modify the internal buffers according to the intensive quanties relevant + * for an element + */ + void processElement(const ElementContext& elemCtx) + { + if (!EWOMS_GET_PARAM(TypeTag, bool, EnableVtkOutput)) + return; + + for (unsigned i = 0; i < elemCtx.numPrimaryDof(/*timeIdx=*/0); ++i) { + unsigned I = elemCtx.globalSpaceIndex(i, /*timeIdx=*/0); + const auto& intQuants = elemCtx.intensiveQuantities(i, /*timeIdx=*/0); + const auto& fs = intQuants.fluidState(); + + if (solidInternalEnergyOutput_()) + solidInternalEnergy_[I] = Toolbox::value(intQuants.solidInternalEnergy()); + if (thermalConductivityOutput_()) + thermalConductivity_[I] = Toolbox::value(intQuants.thermalConductivity()); + + for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) { + if (enthalpyOutput_()) + enthalpy_[phaseIdx][I] = Toolbox::value(fs.enthalpy(phaseIdx)); + if (internalEnergyOutput_()) + internalEnergy_[phaseIdx][I] = Toolbox::value(fs.internalEnergy(phaseIdx)); + } + } + } + + /*! + * \brief Add all buffers to the VTK output writer. + */ + void commitBuffers(BaseOutputWriter& baseWriter) + { + VtkMultiWriter *vtkWriter = dynamic_cast(&baseWriter); + if (!vtkWriter) { + return; + } + + if (solidInternalEnergyOutput_()) + this->commitScalarBuffer_(baseWriter, "internalEnergySolid", solidInternalEnergy_); + if (thermalConductivityOutput_()) + this->commitScalarBuffer_(baseWriter, "thermalConductivity", thermalConductivity_); + + if (enthalpyOutput_()) + this->commitPhaseBuffer_(baseWriter, "enthalpy_%s", enthalpy_); + if (internalEnergyOutput_()) + this->commitPhaseBuffer_(baseWriter, "internalEnergy_%s", internalEnergy_); + } + +private: + static bool solidInternalEnergyOutput_() + { + static bool val = EWOMS_GET_PARAM(TypeTag, bool, VtkWriteSolidInternalEnergy); + return val; + } + + static bool thermalConductivityOutput_() + { + static bool val = EWOMS_GET_PARAM(TypeTag, bool, VtkWriteThermalConductivity); + return val; + } + + static bool enthalpyOutput_() + { + static bool val = EWOMS_GET_PARAM(TypeTag, bool, VtkWriteEnthalpies); + return val; + } + + static bool internalEnergyOutput_() + { + static bool val = EWOMS_GET_PARAM(TypeTag, bool, VtkWriteInternalEnergies); + return val; + } + + PhaseBuffer enthalpy_; + PhaseBuffer internalEnergy_; + + ScalarBuffer thermalConductivity_; + ScalarBuffer solidInternalEnergy_; +}; + +} // namespace Opm + +#endif diff --git a/opm/models/io/vtkmultiphasemodule.hh b/opm/models/io/vtkmultiphasemodule.hh new file mode 100644 index 000000000..2cc4763b5 --- /dev/null +++ b/opm/models/io/vtkmultiphasemodule.hh @@ -0,0 +1,482 @@ +// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- +// vi: set et ts=4 sw=4 sts=4: +/* + This file is part of the Open Porous Media project (OPM). + + OPM is free software: you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation, either version 2 of the License, or + (at your option) any later version. + + OPM is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU General Public License + along with OPM. If not, see . + + Consult the COPYING file in the top-level source directory of this + module for the precise wording of the license and the list of + copyright holders. +*/ +/*! + * \file + * \copydoc Opm::VtkMultiPhaseModule + */ +#ifndef EWOMS_VTK_MULTI_PHASE_MODULE_HH +#define EWOMS_VTK_MULTI_PHASE_MODULE_HH + +#include "vtkmultiwriter.hh" +#include "baseoutputmodule.hh" + +#include +#include + +#include +#include + +#include + +#include + +BEGIN_PROPERTIES + +// create new type tag for the VTK multi-phase output +NEW_TYPE_TAG(VtkMultiPhase); + +// create the property tags needed for the multi phase module +NEW_PROP_TAG(VtkWriteExtrusionFactor); +NEW_PROP_TAG(VtkWritePressures); +NEW_PROP_TAG(VtkWriteDensities); +NEW_PROP_TAG(VtkWriteSaturations); +NEW_PROP_TAG(VtkWriteMobilities); +NEW_PROP_TAG(VtkWriteRelativePermeabilities); +NEW_PROP_TAG(VtkWriteViscosities); +NEW_PROP_TAG(VtkWriteAverageMolarMasses); +NEW_PROP_TAG(VtkWritePorosity); +NEW_PROP_TAG(VtkWriteIntrinsicPermeabilities); +NEW_PROP_TAG(VtkWritePotentialGradients); +NEW_PROP_TAG(VtkWriteFilterVelocities); +NEW_PROP_TAG(VtkOutputFormat); +NEW_PROP_TAG(EnableVtkOutput); + +// set default values for what quantities to output +SET_BOOL_PROP(VtkMultiPhase, VtkWriteExtrusionFactor, false); +SET_BOOL_PROP(VtkMultiPhase, VtkWritePressures, true); +SET_BOOL_PROP(VtkMultiPhase, VtkWriteDensities, true); +SET_BOOL_PROP(VtkMultiPhase, VtkWriteSaturations, true); +SET_BOOL_PROP(VtkMultiPhase, VtkWriteMobilities, false); +SET_BOOL_PROP(VtkMultiPhase, VtkWriteRelativePermeabilities, true); +SET_BOOL_PROP(VtkMultiPhase, VtkWriteViscosities, false); +SET_BOOL_PROP(VtkMultiPhase, VtkWriteAverageMolarMasses, false); +SET_BOOL_PROP(VtkMultiPhase, VtkWritePorosity, true); +SET_BOOL_PROP(VtkMultiPhase, VtkWriteIntrinsicPermeabilities, false); +SET_BOOL_PROP(VtkMultiPhase, VtkWritePotentialGradients, false); +SET_BOOL_PROP(VtkMultiPhase, VtkWriteFilterVelocities, false); + +END_PROPERTIES + +namespace Opm { + +/*! + * \ingroup Vtk + * + * \brief VTK output module for quantities which make sense for all + * models which deal with multiple fluid phases in porous media + * that don't use flashy concepts like interfacial area. + * + * This module deals with the following quantities: + * - Pressures of all fluid phases + * - Densities of all fluid phases + * - Saturations of all fluid phases + * - Mobilities of all fluid phases + * - Relative permeabilities of all fluid phases + * - Viscosities of all fluid phases + * - Average molar masses of all fluid phases + * - Porosity of the medium + * - Norm of the intrinsic permeability of the medium + */ +template +class VtkMultiPhaseModule : public BaseOutputModule +{ + typedef BaseOutputModule ParentType; + + typedef typename GET_PROP_TYPE(TypeTag, Simulator) Simulator; + typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar; + typedef typename GET_PROP_TYPE(TypeTag, ElementContext) ElementContext; + + typedef typename GET_PROP_TYPE(TypeTag, GridView) GridView; + typedef typename GET_PROP_TYPE(TypeTag, FluidSystem) FluidSystem; + typedef typename GET_PROP_TYPE(TypeTag, DiscBaseOutputModule) DiscBaseOutputModule; + + static const int vtkFormat = GET_PROP_VALUE(TypeTag, VtkOutputFormat); + typedef Opm::VtkMultiWriter VtkMultiWriter; + + enum { dimWorld = GridView::dimensionworld }; + enum { numPhases = GET_PROP_VALUE(TypeTag, NumPhases) }; + + typedef typename ParentType::ScalarBuffer ScalarBuffer; + typedef typename ParentType::VectorBuffer VectorBuffer; + typedef typename ParentType::TensorBuffer TensorBuffer; + typedef typename ParentType::PhaseBuffer PhaseBuffer; + + typedef Dune::FieldVector DimVector; + + typedef std::array PhaseVectorBuffer; + +public: + VtkMultiPhaseModule(const Simulator& simulator) + : ParentType(simulator) + {} + + /*! + * \brief Register all run-time parameters for the multi-phase VTK output module. + */ + static void registerParameters() + { + EWOMS_REGISTER_PARAM(TypeTag, bool, VtkWriteExtrusionFactor, + "Include the extrusion factor of the degrees of freedom into the VTK output files"); + EWOMS_REGISTER_PARAM(TypeTag, bool, VtkWritePressures, + "Include the phase pressures in the VTK output files"); + EWOMS_REGISTER_PARAM(TypeTag, bool, VtkWriteDensities, + "Include the phase densities in the VTK output files"); + EWOMS_REGISTER_PARAM(TypeTag, bool, VtkWriteSaturations, + "Include the phase saturations in the VTK output files"); + EWOMS_REGISTER_PARAM(TypeTag, bool, VtkWriteMobilities, + "Include the phase mobilities in the VTK output files"); + EWOMS_REGISTER_PARAM(TypeTag, bool, VtkWriteRelativePermeabilities, + "Include the phase relative permeabilities in the VTK output files"); + EWOMS_REGISTER_PARAM(TypeTag, bool, VtkWriteViscosities, + "Include component phase viscosities in the VTK output files"); + EWOMS_REGISTER_PARAM(TypeTag, bool, VtkWriteAverageMolarMasses, + "Include the average phase mass in the VTK output files"); + EWOMS_REGISTER_PARAM(TypeTag, bool, VtkWritePorosity, + "Include the porosity in the VTK output files"); + EWOMS_REGISTER_PARAM(TypeTag, bool, VtkWriteIntrinsicPermeabilities, + "Include the intrinsic permeability in the VTK output files"); + EWOMS_REGISTER_PARAM(TypeTag, bool, VtkWriteFilterVelocities, + "Include in the filter velocities of the phases the VTK output files"); + EWOMS_REGISTER_PARAM(TypeTag, bool, VtkWritePotentialGradients, + "Include the phase pressure potential gradients in the VTK output files"); + } + + /*! + * \brief Allocate memory for the scalar fields we would like to + * write to the VTK file. + */ + void allocBuffers() + { + if (extrusionFactorOutput_()) this->resizeScalarBuffer_(extrusionFactor_); + if (pressureOutput_()) this->resizePhaseBuffer_(pressure_); + if (densityOutput_()) this->resizePhaseBuffer_(density_); + if (saturationOutput_()) this->resizePhaseBuffer_(saturation_); + if (mobilityOutput_()) this->resizePhaseBuffer_(mobility_); + if (relativePermeabilityOutput_()) this->resizePhaseBuffer_(relativePermeability_); + if (viscosityOutput_()) this->resizePhaseBuffer_(viscosity_); + if (averageMolarMassOutput_()) this->resizePhaseBuffer_(averageMolarMass_); + + if (porosityOutput_()) this->resizeScalarBuffer_(porosity_); + if (intrinsicPermeabilityOutput_()) this->resizeTensorBuffer_(intrinsicPermeability_); + + if (velocityOutput_()) { + size_t nDof = this->simulator_.model().numGridDof(); + for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++ phaseIdx) { + velocity_[phaseIdx].resize(nDof); + for (unsigned dofIdx = 0; dofIdx < nDof; ++ dofIdx) { + velocity_[phaseIdx][dofIdx].resize(dimWorld); + velocity_[phaseIdx][dofIdx] = 0.0; + } + } + this->resizePhaseBuffer_(velocityWeight_); + } + + if (potentialGradientOutput_()) { + size_t nDof = this->simulator_.model().numGridDof(); + for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++ phaseIdx) { + potentialGradient_[phaseIdx].resize(nDof); + for (unsigned dofIdx = 0; dofIdx < nDof; ++ dofIdx) { + potentialGradient_[phaseIdx][dofIdx].resize(dimWorld); + potentialGradient_[phaseIdx][dofIdx] = 0.0; + } + } + + this->resizePhaseBuffer_(potentialWeight_); + } + } + + /*! + * \brief Modify the internal buffers according to the intensive quantities seen on + * an element + */ + void processElement(const ElementContext& elemCtx) + { + if (!EWOMS_GET_PARAM(TypeTag, bool, EnableVtkOutput)) + return; + + const auto& problem = elemCtx.problem(); + for (unsigned i = 0; i < elemCtx.numPrimaryDof(/*timeIdx=*/0); ++i) { + unsigned I = elemCtx.globalSpaceIndex(i, /*timeIdx=*/0); + const auto& intQuants = elemCtx.intensiveQuantities(i, /*timeIdx=*/0); + const auto& fs = intQuants.fluidState(); + + if (extrusionFactorOutput_()) extrusionFactor_[I] = intQuants.extrusionFactor(); + if (porosityOutput_()) porosity_[I] = Opm::getValue(intQuants.porosity()); + + if (intrinsicPermeabilityOutput_()) { + const auto& K = problem.intrinsicPermeability(elemCtx, i, /*timeIdx=*/0); + for (unsigned rowIdx = 0; rowIdx < K.rows; ++rowIdx) + for (unsigned colIdx = 0; colIdx < K.cols; ++colIdx) + intrinsicPermeability_[I][rowIdx][colIdx] = K[rowIdx][colIdx]; + } + + for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) { + if (pressureOutput_()) + pressure_[phaseIdx][I] = Opm::getValue(fs.pressure(phaseIdx)); + if (densityOutput_()) + density_[phaseIdx][I] = Opm::getValue(fs.density(phaseIdx)); + if (saturationOutput_()) + saturation_[phaseIdx][I] = Opm::getValue(fs.saturation(phaseIdx)); + if (mobilityOutput_()) + mobility_[phaseIdx][I] = Opm::getValue(intQuants.mobility(phaseIdx)); + if (relativePermeabilityOutput_()) + relativePermeability_[phaseIdx][I] = Opm::getValue(intQuants.relativePermeability(phaseIdx)); + if (viscosityOutput_()) + viscosity_[phaseIdx][I] = Opm::getValue(fs.viscosity(phaseIdx)); + if (averageMolarMassOutput_()) + averageMolarMass_[phaseIdx][I] = Opm::getValue(fs.averageMolarMass(phaseIdx)); + } + } + + if (potentialGradientOutput_()) { + // calculate velocities if requested + for (unsigned faceIdx = 0; faceIdx < elemCtx.numInteriorFaces(/*timeIdx=*/0); ++ faceIdx) { + const auto& extQuants = elemCtx.extensiveQuantities(faceIdx, /*timeIdx=*/0); + + unsigned i = extQuants.interiorIndex(); + unsigned I = elemCtx.globalSpaceIndex(i, /*timeIdx=*/0); + + for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) { + Scalar weight = extQuants.extrusionFactor(); + + potentialWeight_[phaseIdx][I] += weight; + + const auto& inputPGrad = extQuants.potentialGrad(phaseIdx); + DimVector pGrad; + for (unsigned dimIdx = 0; dimIdx < dimWorld; ++dimIdx) + pGrad[dimIdx] = Opm::getValue(inputPGrad[dimIdx])*weight; + potentialGradient_[phaseIdx][I] += pGrad; + } // end for all phases + } // end for all faces + } + + if (velocityOutput_()) { + // calculate velocities if requested + for (unsigned faceIdx = 0; faceIdx < elemCtx.numInteriorFaces(/*timeIdx=*/0); ++ faceIdx) { + const auto& extQuants = elemCtx.extensiveQuantities(faceIdx, /*timeIdx=*/0); + + unsigned i = extQuants.interiorIndex(); + unsigned I = elemCtx.globalSpaceIndex(i, /*timeIdx=*/0); + + unsigned j = extQuants.exteriorIndex(); + unsigned J = elemCtx.globalSpaceIndex(j, /*timeIdx=*/0); + + for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) { + Scalar weight = std::max(1e-16, + std::abs(Opm::getValue(extQuants.volumeFlux(phaseIdx)))); + Opm::Valgrind::CheckDefined(extQuants.extrusionFactor()); + assert(extQuants.extrusionFactor() > 0); + weight *= extQuants.extrusionFactor(); + + const auto& inputV = extQuants.filterVelocity(phaseIdx); + DimVector v; + for (unsigned k = 0; k < dimWorld; ++k) + v[k] = Opm::getValue(inputV[k]); + if (v.two_norm() > 1e-20) + weight /= v.two_norm(); + v *= weight; + + velocity_[phaseIdx][I] += v; + velocity_[phaseIdx][J] += v; + + velocityWeight_[phaseIdx][I] += weight; + velocityWeight_[phaseIdx][J] += weight; + } // end for all phases + } // end for all faces + } + } + + /*! + * \brief Add all buffers to the VTK output writer. + */ + void commitBuffers(BaseOutputWriter& baseWriter) + { + VtkMultiWriter *vtkWriter = dynamic_cast(&baseWriter); + if (!vtkWriter) + return; + + if (extrusionFactorOutput_()) + this->commitScalarBuffer_(baseWriter, "extrusionFactor", extrusionFactor_); + if (pressureOutput_()) + this->commitPhaseBuffer_(baseWriter, "pressure_%s", pressure_); + if (densityOutput_()) + this->commitPhaseBuffer_(baseWriter, "density_%s", density_); + if (saturationOutput_()) + this->commitPhaseBuffer_(baseWriter, "saturation_%s", saturation_); + if (mobilityOutput_()) + this->commitPhaseBuffer_(baseWriter, "mobility_%s", mobility_); + if (relativePermeabilityOutput_()) + this->commitPhaseBuffer_(baseWriter, "relativePerm_%s", relativePermeability_); + if (viscosityOutput_()) + this->commitPhaseBuffer_(baseWriter, "viscosity_%s", viscosity_); + if (averageMolarMassOutput_()) + this->commitPhaseBuffer_(baseWriter, "averageMolarMass_%s", averageMolarMass_); + + if (porosityOutput_()) + this->commitScalarBuffer_(baseWriter, "porosity", porosity_); + if (intrinsicPermeabilityOutput_()) + this->commitTensorBuffer_(baseWriter, "intrinsicPerm", intrinsicPermeability_); + + if (velocityOutput_()) { + size_t numDof = this->simulator_.model().numGridDof(); + + for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) { + // first, divide the velocity field by the + // respective finite volume's surface area + for (unsigned i = 0; i < numDof; ++i) + velocity_[phaseIdx][i] /= velocityWeight_[phaseIdx][i]; + // commit the phase velocity + char name[512]; + snprintf(name, 512, "filterVelocity_%s", FluidSystem::phaseName(phaseIdx)); + + DiscBaseOutputModule::attachVectorDofData_(baseWriter, velocity_[phaseIdx], name); + } + } + + if (potentialGradientOutput_()) { + size_t numDof = this->simulator_.model().numGridDof(); + + for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) { + // first, divide the velocity field by the + // respective finite volume's surface area + for (unsigned i = 0; i < numDof; ++i) + potentialGradient_[phaseIdx][i] /= potentialWeight_[phaseIdx][i]; + // commit the phase velocity + char name[512]; + snprintf(name, 512, "gradP_%s", FluidSystem::phaseName(phaseIdx)); + + DiscBaseOutputModule::attachVectorDofData_(baseWriter, + potentialGradient_[phaseIdx], + name); + } + } + } + + /*! + * \brief Returns true iff the module needs to access the extensive quantities of a + * context to do its job. + * + * For example, this happens if velocities or gradients should be written. Always + * returning true here does not do any harm from the correctness perspective, but it + * slows down writing the output fields. + */ + virtual bool needExtensiveQuantities() const final + { + return velocityOutput_() || potentialGradientOutput_(); + } + +private: + static bool extrusionFactorOutput_() + { + static bool val = EWOMS_GET_PARAM(TypeTag, bool, VtkWriteExtrusionFactor); + return val; + } + + static bool pressureOutput_() + { + static bool val = EWOMS_GET_PARAM(TypeTag, bool, VtkWritePressures); + return val; + } + + static bool densityOutput_() + { + static bool val = EWOMS_GET_PARAM(TypeTag, bool, VtkWriteDensities); + return val; + } + + static bool saturationOutput_() + { + static bool val = EWOMS_GET_PARAM(TypeTag, bool, VtkWriteSaturations); + return val; + } + + static bool mobilityOutput_() + { + static bool val = EWOMS_GET_PARAM(TypeTag, bool, VtkWriteMobilities); + return val; + } + + static bool relativePermeabilityOutput_() + { + static bool val = EWOMS_GET_PARAM(TypeTag, bool, VtkWriteRelativePermeabilities); + return val; + } + + static bool viscosityOutput_() + { + static bool val = EWOMS_GET_PARAM(TypeTag, bool, VtkWriteViscosities); + return val; + } + + static bool averageMolarMassOutput_() + { + static bool val = EWOMS_GET_PARAM(TypeTag, bool, VtkWriteAverageMolarMasses); + return val; + } + + static bool porosityOutput_() + { + static bool val = EWOMS_GET_PARAM(TypeTag, bool, VtkWritePorosity); + return val; + } + + static bool intrinsicPermeabilityOutput_() + { + static bool val = EWOMS_GET_PARAM(TypeTag, bool, VtkWriteIntrinsicPermeabilities); + return val; + } + + static bool velocityOutput_() + { + static bool val = EWOMS_GET_PARAM(TypeTag, bool, VtkWriteFilterVelocities); + return val; + } + + static bool potentialGradientOutput_() + { + static bool val = EWOMS_GET_PARAM(TypeTag, bool, VtkWritePotentialGradients); + return val; + } + + ScalarBuffer extrusionFactor_; + PhaseBuffer pressure_; + PhaseBuffer density_; + PhaseBuffer saturation_; + PhaseBuffer mobility_; + PhaseBuffer relativePermeability_; + PhaseBuffer viscosity_; + PhaseBuffer averageMolarMass_; + + ScalarBuffer porosity_; + TensorBuffer intrinsicPermeability_; + + PhaseVectorBuffer velocity_; + PhaseBuffer velocityWeight_; + + PhaseVectorBuffer potentialGradient_; + PhaseBuffer potentialWeight_; +}; + +} // namespace Opm + +#endif diff --git a/opm/models/io/vtkmultiwriter.hh b/opm/models/io/vtkmultiwriter.hh new file mode 100644 index 000000000..884efc7b7 --- /dev/null +++ b/opm/models/io/vtkmultiwriter.hh @@ -0,0 +1,575 @@ +// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- +// vi: set et ts=4 sw=4 sts=4: +/* + This file is part of the Open Porous Media project (OPM). + + OPM is free software: you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation, either version 2 of the License, or + (at your option) any later version. + + OPM is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU General Public License + along with OPM. If not, see . + + Consult the COPYING file in the top-level source directory of this + module for the precise wording of the license and the list of + copyright holders. +*/ +/*! + * \file + * + * \copydoc Opm::VtkMultiWriter + */ +#ifndef EWOMS_VTK_MULTI_WRITER_HH +#define EWOMS_VTK_MULTI_WRITER_HH + +#include "vtkscalarfunction.hh" +#include "vtkvectorfunction.hh" +#include "vtktensorfunction.hh" + +#include +#include + +#include +#include + +#include +#include +#include + +#if HAVE_MPI +#include +#endif + +#include +#include +#include +#include +#include + +namespace Opm { +/*! + * \brief Simplifies writing multi-file VTK datasets. + * + * This class automatically keeps the meta file up to date and + * simplifies writing datasets consisting of multiple files. (i.e. + * multiple time steps or grid refinements within a time step.) + */ +template +class VtkMultiWriter : public BaseOutputWriter +{ + class WriteDataTasklet : public TaskletInterface + { + public: + WriteDataTasklet(VtkMultiWriter& multiWriter) + : multiWriter_(multiWriter) + { } + + void run() final + { + std::string fileName; + // write the actual data as vtu or vtp (plus the pieces file in the parallel case) + if (multiWriter_.commSize_ > 1) + fileName = multiWriter_.curWriter_->pwrite(/*name=*/multiWriter_.curOutFileName_, + /*path=*/multiWriter_.outputDir_, + /*extendPath=*/"", + static_cast(vtkFormat)); + else + fileName = multiWriter_.curWriter_->write(/*name=*/multiWriter_.outputDir_ + "/" + multiWriter_.curOutFileName_, + static_cast(vtkFormat)); + + // determine name to write into the multi-file for the + // current time step + multiWriter_.multiFile_.precision(16); + multiWriter_.multiFile_ << " \n"; + } + + private: + VtkMultiWriter& multiWriter_; + }; + + enum { dim = GridView::dimension }; + +#if DUNE_VERSION_NEWER(DUNE_GRID, 2,6) + typedef Dune::MultipleCodimMultipleGeomTypeMapper VertexMapper; + typedef Dune::MultipleCodimMultipleGeomTypeMapper ElementMapper; +#else + typedef Dune::MultipleCodimMultipleGeomTypeMapper ElementMapper; + typedef Dune::MultipleCodimMultipleGeomTypeMapper VertexMapper; +#endif + +public: + typedef BaseOutputWriter::Scalar Scalar; + typedef BaseOutputWriter::Vector Vector; + typedef BaseOutputWriter::Tensor Tensor; + typedef BaseOutputWriter::ScalarBuffer ScalarBuffer; + typedef BaseOutputWriter::VectorBuffer VectorBuffer; + typedef BaseOutputWriter::TensorBuffer TensorBuffer; + + typedef Dune::VTKWriter VtkWriter; +#if DUNE_VERSION_NEWER(DUNE_GRID, 2,5) + typedef std::shared_ptr< Dune::VTKFunction< GridView > > FunctionPtr; +#else + typedef typename VtkWriter::VTKFunctionPtr FunctionPtr; +#endif + + VtkMultiWriter(bool asyncWriting, + const GridView& gridView, + const std::string& outputDir, + const std::string& simName = "", + std::string multiFileName = "") + : gridView_(gridView) +#if DUNE_VERSION_NEWER(DUNE_GRID, 2,6) + , elementMapper_(gridView, Dune::mcmgElementLayout()) + , vertexMapper_(gridView, Dune::mcmgVertexLayout()) +#else + , elementMapper_(gridView) + , vertexMapper_(gridView) +#endif + , curWriter_(nullptr) + , curWriterNum_(0) + , taskletRunner_(/*numThreads=*/asyncWriting?1:0) + { + outputDir_ = outputDir; + if (outputDir == "") + outputDir_ = "."; + + simName_ = (simName.empty()) ? "sim" : simName; + multiFileName_ = multiFileName; + if (multiFileName_.empty()) + multiFileName_ = outputDir_+"/"+simName_+".pvd"; + + commRank_ = gridView.comm().rank(); + commSize_ = gridView.comm().size(); + } + + ~VtkMultiWriter() + { + taskletRunner_.barrier(); + releaseBuffers_(); + finishMultiFile_(); + + if (commRank_ == 0) + multiFile_.close(); + } + + /*! + * \brief Returns the number of the current VTK file. + */ + int curWriterNum() const + { return curWriterNum_; } + + /*! + * \brief Updates the internal data structures after mesh + * refinement. + * + * If the grid changes between two calls of beginWrite(), this + * method _must_ be called before the second beginWrite()! + */ + void gridChanged() + { + elementMapper_.update(); + vertexMapper_.update(); + } + + /*! + * \brief Called whenever a new time step must be written. + */ + void beginWrite(double t) + { + if (!multiFile_.is_open()) { + startMultiFile_(multiFileName_); + } + + // make sure that all previous output has been written and no other thread + // accesses the memory used as the target for the extracted quantities + taskletRunner_.barrier(); + releaseBuffers_(); + + curTime_ = t; + curOutFileName_ = fileName_(); + + curWriter_ = new VtkWriter(gridView_, Dune::VTK::conforming); + ++curWriterNum_; + } + + /*! + * \brief Allocate a managed buffer for a scalar field + * + * The buffer will be deleted automatically after the data has + * been written by to disk. + */ + ScalarBuffer *allocateManagedScalarBuffer(size_t numEntities) + { + ScalarBuffer *buf = new ScalarBuffer(numEntities); + managedScalarBuffers_.push_back(buf); + return buf; + } + + /*! + * \brief Allocate a managed buffer for a vector field + * + * The buffer will be deleted automatically after the data has + * been written by to disk. + */ + VectorBuffer *allocateManagedVectorBuffer(size_t numOuter, size_t numInner) + { + VectorBuffer *buf = new VectorBuffer(numOuter); + for (size_t i = 0; i < numOuter; ++ i) + (*buf)[i].resize(numInner); + + managedVectorBuffers_.push_back(buf); + return buf; + } + + /*! + * \brief Add a finished vertex centered vector field to the + * output. + * + * If the buffer is managed by the VtkMultiWriter, it must have + * been created using allocateManagedBuffer() and may not be used + * anywhere after calling this method. After the data is written + * to disk, it will be deleted automatically. + * + * If the buffer is not managed by the MultiWriter, the buffer + * must exist at least until the call to endWrite() + * finishes. + * + * In both cases, modifying the buffer between the call to this + * method and endWrite() results in _undefined behavior_. + */ + void attachScalarVertexData(ScalarBuffer& buf, std::string name) + { + sanitizeScalarBuffer_(buf); + + typedef Opm::VtkScalarFunction VtkFn; + FunctionPtr fnPtr(new VtkFn(name, + gridView_, + vertexMapper_, + buf, + /*codim=*/dim)); + curWriter_->addVertexData(fnPtr); + } + + /*! + * \brief Add a element centered quantity to the output. + * + * If the buffer is managed by the VtkMultiWriter, it must have + * been created using createField() and may not be used by + * anywhere after calling this method. After the data is written + * to disk, it will be deleted automatically. + * + * If the buffer is not managed by the MultiWriter, the buffer + * must exist at least until the call to endWrite() + * finishes. + * + * In both cases, modifying the buffer between the call to this + * method and endWrite() results in _undefined behaviour_. + */ + void attachScalarElementData(ScalarBuffer& buf, std::string name) + { + sanitizeScalarBuffer_(buf); + + typedef Opm::VtkScalarFunction VtkFn; + FunctionPtr fnPtr(new VtkFn(name, + gridView_, + elementMapper_, + buf, + /*codim=*/0)); + curWriter_->addCellData(fnPtr); + } + + /*! + * \brief Add a finished vertex centered vector field to the + * output. + * + * If the buffer is managed by the VtkMultiWriter, it must have + * been created using allocateManagedBuffer() and may not be used + * anywhere after calling this method. After the data is written + * to disk, it will be deleted automatically. + * + * If the buffer is not managed by the MultiWriter, the buffer + * must exist at least until the call to endWrite() + * finishes. + * + * In both cases, modifying the buffer between the call to this + * method and endWrite() results in _undefined behavior_. + */ + void attachVectorVertexData(VectorBuffer& buf, std::string name) + { + sanitizeVectorBuffer_(buf); + + typedef Opm::VtkVectorFunction VtkFn; + FunctionPtr fnPtr(new VtkFn(name, + gridView_, + vertexMapper_, + buf, + /*codim=*/dim)); + curWriter_->addVertexData(fnPtr); + } + + /*! + * \brief Add a finished vertex-centered tensor field to the output. + */ + void attachTensorVertexData(TensorBuffer& buf, std::string name) + { + typedef Opm::VtkTensorFunction VtkFn; + + for (unsigned colIdx = 0; colIdx < buf[0].N(); ++colIdx) { + std::ostringstream oss; + oss << name << "[" << colIdx << "]"; + + FunctionPtr fnPtr(new VtkFn(oss.str(), + gridView_, + vertexMapper_, + buf, + /*codim=*/dim, + colIdx)); + curWriter_->addVertexData(fnPtr); + } + } + + /*! + * \brief Add a element centered quantity to the output. + * + * If the buffer is managed by the VtkMultiWriter, it must have + * been created using createField() and may not be used by + * anywhere after calling this method. After the data is written + * to disk, it will be deleted automatically. + * + * If the buffer is not managed by the MultiWriter, the buffer + * must exist at least until the call to endWrite() + * finishes. + * + * In both cases, modifying the buffer between the call to this + * method and endWrite() results in _undefined behaviour_. + */ + void attachVectorElementData(VectorBuffer& buf, std::string name) + { + sanitizeVectorBuffer_(buf); + + typedef Opm::VtkVectorFunction VtkFn; + FunctionPtr fnPtr(new VtkFn(name, + gridView_, + elementMapper_, + buf, + /*codim=*/0)); + curWriter_->addCellData(fnPtr); + } + + /*! + * \brief Add a finished element-centered tensor field to the output. + */ + void attachTensorElementData(TensorBuffer& buf, std::string name) + { + typedef Opm::VtkTensorFunction VtkFn; + + for (unsigned colIdx = 0; colIdx < buf[0].N(); ++colIdx) { + std::ostringstream oss; + oss << name << "[" << colIdx << "]"; + + FunctionPtr fnPtr(new VtkFn(oss.str(), + gridView_, + elementMapper_, + buf, + /*codim=*/0, + colIdx)); + curWriter_->addCellData(fnPtr); + } + } + + /*! + * \brief Finalizes the current writer. + * + * This means that everything will be written to disk, except if + * the onlyDiscard argument is true. In this case only all managed + * buffers are deleted, but no output is written. + */ + void endWrite(bool onlyDiscard = false) + { + if (!onlyDiscard) { + auto tasklet = std::make_shared(*this); + taskletRunner_.dispatch(tasklet); + } + else + --curWriterNum_; + + // temporarily write the closing XML mumbo-jumbo to the mashup + // file so that the data set can be loaded even if the + // simulation is aborted (or not yet finished) + finishMultiFile_(); + } + + /*! + * \brief Write the multi-writer's state to a restart file. + */ + template + void serialize(Restarter& res) + { + res.serializeSectionBegin("VTKMultiWriter"); + res.serializeStream() << curWriterNum_ << "\n"; + + if (commRank_ == 0) { + std::streamsize fileLen = 0; + std::streamoff filePos = 0; + if (multiFile_.is_open()) { + // write the meta file into the restart file + filePos = multiFile_.tellp(); + multiFile_.seekp(0, std::ios::end); + fileLen = multiFile_.tellp(); + multiFile_.seekp(filePos); + } + + res.serializeStream() << fileLen << " " << filePos << "\n"; + + if (fileLen > 0) { + std::ifstream multiFileIn(multiFileName_.c_str()); + char *tmp = new char[fileLen]; + multiFileIn.read(tmp, static_cast(fileLen)); + res.serializeStream().write(tmp, fileLen); + delete[] tmp; + } + } + + res.serializeSectionEnd(); + } + + /*! + * \brief Read the multi-writer's state from a restart file. + */ + template + void deserialize(Restarter& res) + { + res.deserializeSectionBegin("VTKMultiWriter"); + res.deserializeStream() >> curWriterNum_; + + if (commRank_ == 0) { + std::string dummy; + std::getline(res.deserializeStream(), dummy); + + // recreate the meta file from the restart file + std::streamoff filePos; + std::streamsize fileLen; + res.deserializeStream() >> fileLen >> filePos; + std::getline(res.deserializeStream(), dummy); + if (multiFile_.is_open()) + multiFile_.close(); + + if (fileLen > 0) { + multiFile_.open(multiFileName_.c_str()); + + char *tmp = new char[fileLen]; + res.deserializeStream().read(tmp, fileLen); + multiFile_.write(tmp, fileLen); + delete[] tmp; + } + + multiFile_.seekp(filePos); + } + else { + std::string tmp; + std::getline(res.deserializeStream(), tmp); + } + res.deserializeSectionEnd(); + } + +private: + std::string fileName_() + { + // use a new file name for each time step + std::ostringstream oss; + oss << simName_ << "-" + << std::setw(5) << std::setfill('0') << curWriterNum_; + return oss.str(); + } + + std::string fileSuffix_() + { return (GridView::dimension == 1) ? "vtp" : "vtu"; } + + void startMultiFile_(const std::string& multiFileName) + { + // only the first process writes to the multi-file + if (commRank_ == 0) { + // generate one meta vtk-file holding the individual time steps + multiFile_.open(multiFileName.c_str()); + multiFile_ << "\n" + "\n" + " \n"; + } + } + + void finishMultiFile_() + { + // only the first process writes to the multi-file + if (commRank_ == 0) { + // make sure that we always have a working meta file + std::ofstream::pos_type pos = multiFile_.tellp(); + multiFile_ << " \n" + "\n"; + multiFile_.seekp(pos); + multiFile_.flush(); + } + } + + // make sure the field is well defined if running under valgrind + // and make sure that all values can be displayed by paraview + void sanitizeScalarBuffer_(ScalarBuffer& b OPM_UNUSED) + { + // nothing to do: this is done by VtkScalarFunction + } + + void sanitizeVectorBuffer_(VectorBuffer& b OPM_UNUSED) + { + // nothing to do: this is done by VtkVectorFunction + } + + // release the memory occupied by all buffer objects managed by the multi-writer + void releaseBuffers_() + { + // discard managed objects and the current VTK writer + delete curWriter_; + curWriter_ = nullptr; + while (managedScalarBuffers_.begin() != managedScalarBuffers_.end()) { + delete managedScalarBuffers_.front(); + managedScalarBuffers_.pop_front(); + } + while (managedVectorBuffers_.begin() != managedVectorBuffers_.end()) { + delete managedVectorBuffers_.front(); + managedVectorBuffers_.pop_front(); + } + } + + const GridView gridView_; + ElementMapper elementMapper_; + VertexMapper vertexMapper_; + + std::string outputDir_; + std::string simName_; + std::ofstream multiFile_; + std::string multiFileName_; + + int commSize_; // number of processes in the communicator + int commRank_; // rank of the current process in the communicator + + VtkWriter *curWriter_; + double curTime_; + std::string curOutFileName_; + int curWriterNum_; + + std::list managedScalarBuffers_; + std::list managedVectorBuffers_; + + TaskletRunner taskletRunner_; +}; +} // namespace Opm + +#endif diff --git a/opm/models/io/vtkphasepresencemodule.hh b/opm/models/io/vtkphasepresencemodule.hh new file mode 100644 index 000000000..21bd832cc --- /dev/null +++ b/opm/models/io/vtkphasepresencemodule.hh @@ -0,0 +1,141 @@ +// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- +// vi: set et ts=4 sw=4 sts=4: +/* + This file is part of the Open Porous Media project (OPM). + + OPM is free software: you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation, either version 2 of the License, or + (at your option) any later version. + + OPM is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU General Public License + along with OPM. If not, see . + + Consult the COPYING file in the top-level source directory of this + module for the precise wording of the license and the list of + copyright holders. +*/ +/*! + * \file + * \copydoc Opm::VtkPhasePresenceModule + */ +#ifndef EWOMS_VTK_PHASE_PRESENCE_MODULE_HH +#define EWOMS_VTK_PHASE_PRESENCE_MODULE_HH + +#include "vtkmultiwriter.hh" +#include "baseoutputmodule.hh" + +#include +#include + +BEGIN_PROPERTIES + +// create new type tag for the VTK primary variables output +NEW_TYPE_TAG(VtkPhasePresence); + +// create the property tags needed for the primary variables module +NEW_PROP_TAG(VtkWritePhasePresence); +NEW_PROP_TAG(VtkOutputFormat); +NEW_PROP_TAG(EnableVtkOutput); + +SET_BOOL_PROP(VtkPhasePresence, VtkWritePhasePresence, false); + +END_PROPERTIES + +namespace Opm { +/*! + * \ingroup Vtk + * + * \brief VTK output module for the fluid composition + */ +template +class VtkPhasePresenceModule : public BaseOutputModule +{ + typedef BaseOutputModule ParentType; + + typedef typename GET_PROP_TYPE(TypeTag, Simulator) Simulator; + typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar; + typedef typename GET_PROP_TYPE(TypeTag, ElementContext) ElementContext; + typedef typename GET_PROP_TYPE(TypeTag, GridView) GridView; + + static const int vtkFormat = GET_PROP_VALUE(TypeTag, VtkOutputFormat); + typedef Opm::VtkMultiWriter VtkMultiWriter; + + typedef typename ParentType::ScalarBuffer ScalarBuffer; + + +public: + VtkPhasePresenceModule(const Simulator& simulator) + : ParentType(simulator) + { } + + /*! + * \brief Register all run-time parameters for the Vtk output module. + */ + static void registerParameters() + { + EWOMS_REGISTER_PARAM(TypeTag, bool, VtkWritePhasePresence, + "Include the phase presence pseudo primary " + "variable in the VTK output files"); + } + + /*! + * \brief Allocate memory for the scalar fields we would like to + * write to the VTK file. + */ + void allocBuffers() + { + if (phasePresenceOutput_()) this->resizeScalarBuffer_(phasePresence_); + } + + /*! + * \brief Modify the internal buffers according to the intensive quanties relevant + * for an element + */ + void processElement(const ElementContext& elemCtx) + { + if (!EWOMS_GET_PARAM(TypeTag, bool, EnableVtkOutput)) + return; + + for (unsigned i = 0; i < elemCtx.numPrimaryDof(/*timeIdx=*/0); ++i) { + // calculate the phase presence + int phasePresence = elemCtx.primaryVars(i, /*timeIdx=*/0).phasePresence(); + unsigned I = elemCtx.globalSpaceIndex(i, /*timeIdx=*/0); + + if (phasePresenceOutput_()) + phasePresence_[I] = phasePresence; + } + } + + /*! + * \brief Add all buffers to the output writer. + */ + void commitBuffers(BaseOutputWriter& baseWriter) + { + VtkMultiWriter *vtkWriter = dynamic_cast(&baseWriter); + if (!vtkWriter) { + return; + } + + if (phasePresenceOutput_()) + this->commitScalarBuffer_(baseWriter, "phase presence", phasePresence_); + } + +private: + static bool phasePresenceOutput_() + { + static bool val = EWOMS_GET_PARAM(TypeTag, bool, VtkWritePhasePresence); + return val; + } + + ScalarBuffer phasePresence_; +}; + +} // namespace Opm + +#endif diff --git a/opm/models/io/vtkprimaryvarsmodule.hh b/opm/models/io/vtkprimaryvarsmodule.hh new file mode 100644 index 000000000..bace701ba --- /dev/null +++ b/opm/models/io/vtkprimaryvarsmodule.hh @@ -0,0 +1,185 @@ +// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- +// vi: set et ts=4 sw=4 sts=4: +/* + This file is part of the Open Porous Media project (OPM). + + OPM is free software: you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation, either version 2 of the License, or + (at your option) any later version. + + OPM is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU General Public License + along with OPM. If not, see . + + Consult the COPYING file in the top-level source directory of this + module for the precise wording of the license and the list of + copyright holders. +*/ +/*! + * \file + * \copydoc Opm::VtkPrimaryVarsModule + */ +#ifndef EWOMS_VTK_PRIMARY_VARS_MODULE_HH +#define EWOMS_VTK_PRIMARY_VARS_MODULE_HH + +#include +#include + +#include +#include + +BEGIN_PROPERTIES + +// create new type tag for the VTK primary variables output +NEW_PROP_TAG(EnableVtkOutput); +NEW_TYPE_TAG(VtkPrimaryVars); + +// create the property tags needed for the primary variables module +NEW_PROP_TAG(VtkWritePrimaryVars); +NEW_PROP_TAG(VtkWriteProcessRank); +NEW_PROP_TAG(VtkWriteDofIndex); +NEW_PROP_TAG(VtkOutputFormat); +NEW_PROP_TAG(EnableVtkOutput); + +SET_BOOL_PROP(VtkPrimaryVars, VtkWritePrimaryVars, false); +SET_BOOL_PROP(VtkPrimaryVars, VtkWriteProcessRank, false); +SET_BOOL_PROP(VtkPrimaryVars, VtkWriteDofIndex, false); + +END_PROPERTIES + +namespace Opm { + +/*! + * \ingroup Vtk + * + * \brief VTK output module for the fluid composition + */ +template +class VtkPrimaryVarsModule : public BaseOutputModule +{ + typedef BaseOutputModule ParentType; + + typedef typename GET_PROP_TYPE(TypeTag, Simulator) Simulator; + typedef typename GET_PROP_TYPE(TypeTag, ElementContext) ElementContext; + typedef typename GET_PROP_TYPE(TypeTag, GridView) GridView; + + static const int vtkFormat = GET_PROP_VALUE(TypeTag, VtkOutputFormat); + typedef Opm::VtkMultiWriter VtkMultiWriter; + + typedef typename ParentType::ScalarBuffer ScalarBuffer; + typedef typename ParentType::EqBuffer EqBuffer; + + enum { numEq = GET_PROP_VALUE(TypeTag, NumEq) }; + +public: + VtkPrimaryVarsModule(const Simulator& simulator) + : ParentType(simulator) + { } + + /*! + * \brief Register all run-time parameters for the Vtk output module. + */ + static void registerParameters() + { + EWOMS_REGISTER_PARAM(TypeTag, bool, VtkWritePrimaryVars, + "Include the primary variables into the VTK output files"); + EWOMS_REGISTER_PARAM(TypeTag, bool, VtkWriteProcessRank, + "Include the MPI process rank into the VTK output files"); + EWOMS_REGISTER_PARAM(TypeTag, bool, VtkWriteDofIndex, + "Include the index of the degrees of freedom into the VTK output files"); + } + + /*! + * \brief Allocate memory for the scalar fields we would like to + * write to the VTK file. + */ + void allocBuffers() + { + if (primaryVarsOutput_()) + this->resizeEqBuffer_(primaryVars_); + if (processRankOutput_()) + this->resizeScalarBuffer_(processRank_, + /*bufferType=*/ParentType::ElementBuffer); + if (dofIndexOutput_()) + this->resizeScalarBuffer_(dofIndex_); + } + + /*! + * \brief Modify the internal buffers according to the intensive quantities relevant for + * an element + */ + void processElement(const ElementContext& elemCtx) + { + if (!EWOMS_GET_PARAM(TypeTag, bool, EnableVtkOutput)) + return; + + const auto& elementMapper = elemCtx.model().elementMapper(); + unsigned elemIdx = static_cast(elementMapper.index(elemCtx.element())); + if (processRankOutput_() && !processRank_.empty()) + processRank_[elemIdx] = static_cast(this->simulator_.gridView().comm().rank()); + + for (unsigned i = 0; i < elemCtx.numPrimaryDof(/*timeIdx=*/0); ++i) { + unsigned I = elemCtx.globalSpaceIndex(i, /*timeIdx=*/0); + const auto& priVars = elemCtx.primaryVars(i, /*timeIdx=*/0); + + if (dofIndexOutput_()) + dofIndex_[I] = I; + + for (unsigned eqIdx = 0; eqIdx < numEq; ++eqIdx) { + if (primaryVarsOutput_() && !primaryVars_[eqIdx].empty()) + primaryVars_[eqIdx][I] = priVars[eqIdx]; + } + } + } + + /*! + * \brief Add all buffers to the VTK output writer. + */ + void commitBuffers(BaseOutputWriter& baseWriter) + { + VtkMultiWriter *vtkWriter = dynamic_cast(&baseWriter); + if (!vtkWriter) { + return; + } + + if (primaryVarsOutput_()) + this->commitPriVarsBuffer_(baseWriter, "PV_%s", primaryVars_); + if (processRankOutput_()) + this->commitScalarBuffer_(baseWriter, + "process rank", + processRank_, + /*bufferType=*/ParentType::ElementBuffer); + if (dofIndexOutput_()) + this->commitScalarBuffer_(baseWriter, "DOF index", dofIndex_); + } + +private: + static bool primaryVarsOutput_() + { + static bool val = EWOMS_GET_PARAM(TypeTag, bool, VtkWritePrimaryVars); + return val; + } + static bool processRankOutput_() + { + static bool val = EWOMS_GET_PARAM(TypeTag, bool, VtkWriteProcessRank); + return val; + } + static bool dofIndexOutput_() + { + static bool val = EWOMS_GET_PARAM(TypeTag, bool, VtkWriteDofIndex); + return val; + } + + EqBuffer primaryVars_; + ScalarBuffer processRank_; + ScalarBuffer dofIndex_; +}; + +} // namespace Opm + +#endif diff --git a/opm/models/io/vtkscalarfunction.hh b/opm/models/io/vtkscalarfunction.hh new file mode 100644 index 000000000..6163a433e --- /dev/null +++ b/opm/models/io/vtkscalarfunction.hh @@ -0,0 +1,126 @@ +// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- +// vi: set et ts=4 sw=4 sts=4: +/* + This file is part of the Open Porous Media project (OPM). + + OPM is free software: you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation, either version 2 of the License, or + (at your option) any later version. + + OPM is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU General Public License + along with OPM. If not, see . + + Consult the COPYING file in the top-level source directory of this + module for the precise wording of the license and the list of + copyright holders. +*/ +/*! + * \file + * \copydoc Opm::VtkScalarFunction + */ +#ifndef VTK_SCALAR_FUNCTION_HH +#define VTK_SCALAR_FUNCTION_HH + +#include + +#include +#include +#include +#include + +#include + +#include + +#include +#include +#include + +namespace Opm { + +/*! + * \brief Provides a vector-valued function using Dune::FieldVectors + * as elements. + */ +template +class VtkScalarFunction : public Dune::VTKFunction +{ + enum { dim = GridView::dimension }; + typedef typename GridView::ctype ctype; + typedef typename GridView::template Codim<0>::Entity Element; + + typedef BaseOutputWriter::ScalarBuffer ScalarBuffer; + +public: + VtkScalarFunction(std::string name, + const GridView& gridView, + const Mapper& mapper, + const ScalarBuffer& buf, + unsigned codim) + : name_(name) + , gridView_(gridView) + , mapper_(mapper) + , buf_(buf) + , codim_(codim) + { assert(int(buf_.size()) == int(mapper_.size())); } + + virtual std::string name() const + { return name_; } + + virtual int ncomps() const + { return 1; } + + virtual double evaluate(int mycomp OPM_UNUSED, + const Element& e, + const Dune::FieldVector& xi) const + { + unsigned idx; + if (codim_ == 0) { + // cells. map element to the index + idx = static_cast(mapper_.index(e)); + } + else if (codim_ == dim) { + // find vertex which is closest to xi in local + // coordinates. This code is based on Dune::P1VTKFunction + double min = 1e100; + int imin = -1; + Dune::GeometryType gt = e.type(); + int n = static_cast(e.subEntities(dim)); + for (int i = 0; i < n; ++i) { + Dune::FieldVector local = + Dune::ReferenceElements::general(gt).position(i, dim); + + local -= xi; + if (local.infinity_norm() < min) { + min = local.infinity_norm(); + imin = static_cast(i); + } + } + + // map vertex to an index + idx = static_cast(mapper_.subIndex(e, imin, codim_)); + } + else + throw std::logic_error("Only element and vertex based vector fields are" + " supported so far."); + + return static_cast(static_cast(buf_[idx])); + } + +private: + const std::string name_; + const GridView gridView_; + const Mapper& mapper_; + const ScalarBuffer& buf_; + unsigned codim_; +}; + +} // namespace Opm + +#endif diff --git a/opm/models/io/vtktemperaturemodule.hh b/opm/models/io/vtktemperaturemodule.hh new file mode 100644 index 000000000..ed9d7e9b0 --- /dev/null +++ b/opm/models/io/vtktemperaturemodule.hh @@ -0,0 +1,147 @@ +// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- +// vi: set et ts=4 sw=4 sts=4: +/* + This file is part of the Open Porous Media project (OPM). + + OPM is free software: you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation, either version 2 of the License, or + (at your option) any later version. + + OPM is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU General Public License + along with OPM. If not, see . + + Consult the COPYING file in the top-level source directory of this + module for the precise wording of the license and the list of + copyright holders. +*/ +/*! + * \file + * \copydoc Opm::VtkTemperatureModule + */ +#ifndef EWOMS_VTK_TEMPERATURE_MODULE_HH +#define EWOMS_VTK_TEMPERATURE_MODULE_HH + +#include "vtkmultiwriter.hh" +#include "baseoutputmodule.hh" + +#include +#include + +#include + +BEGIN_PROPERTIES + +// create new type tag for the VTK temperature output +NEW_TYPE_TAG(VtkTemperature); + +// create the property tags needed for the temperature module +NEW_PROP_TAG(VtkWriteTemperature); +NEW_PROP_TAG(VtkOutputFormat); +NEW_PROP_TAG(EnableVtkOutput); + +// set default values for what quantities to output +SET_BOOL_PROP(VtkTemperature, VtkWriteTemperature, true); + +END_PROPERTIES + +namespace Opm { + +/*! + * \ingroup Vtk + * + * \brief VTK output module for the temperature in which assume + * thermal equilibrium + */ +template +class VtkTemperatureModule : public BaseOutputModule +{ + typedef BaseOutputModule ParentType; + + typedef typename GET_PROP_TYPE(TypeTag, Simulator) Simulator; + typedef typename GET_PROP_TYPE(TypeTag, ElementContext) ElementContext; + + typedef typename GET_PROP_TYPE(TypeTag, GridView) GridView; + typedef typename GET_PROP_TYPE(TypeTag, Evaluation) Evaluation; + + typedef typename ParentType::ScalarBuffer ScalarBuffer; + + static const int vtkFormat = GET_PROP_VALUE(TypeTag, VtkOutputFormat); + typedef Opm::VtkMultiWriter VtkMultiWriter; + +public: + VtkTemperatureModule(const Simulator& simulator) + : ParentType(simulator) + {} + + /*! + * \brief Register all run-time parameters for the Vtk output module. + */ + static void registerParameters() + { + EWOMS_REGISTER_PARAM(TypeTag, bool, VtkWriteTemperature, + "Include the temperature in the VTK output files"); + } + + /*! + * \brief Allocate memory for the scalar fields we would like to + * write to the VTK file. + */ + void allocBuffers() + { + if (temperatureOutput_()) this->resizeScalarBuffer_(temperature_); + } + + /*! + * \brief Modify the internal buffers according to the intensive quantities relevant + * for an element + */ + void processElement(const ElementContext& elemCtx) + { + typedef Opm::MathToolbox Toolbox; + + if (!EWOMS_GET_PARAM(TypeTag, bool, EnableVtkOutput)) + return; + + for (unsigned i = 0; i < elemCtx.numPrimaryDof(/*timeIdx=*/0); ++i) { + unsigned I = elemCtx.globalSpaceIndex(i, /*timeIdx=*/0); + const auto& intQuants = elemCtx.intensiveQuantities(i, /*timeIdx=*/0); + const auto& fs = intQuants.fluidState(); + + if (temperatureOutput_()) + temperature_[I] = Toolbox::value(fs.temperature(/*phaseIdx=*/0)); + } + } + + /*! + * \brief Add all buffers to the VTK output writer. + */ + void commitBuffers(BaseOutputWriter& baseWriter) + { + VtkMultiWriter *vtkWriter = dynamic_cast(&baseWriter); + if (!vtkWriter) { + return; + } + + if (temperatureOutput_()) + this->commitScalarBuffer_(baseWriter, "temperature", temperature_); + } + +private: + static bool temperatureOutput_() + { + static bool val = EWOMS_GET_PARAM(TypeTag, bool, VtkWriteTemperature); + return val; + } + + ScalarBuffer temperature_; +}; + +} // namespace Opm + +#endif diff --git a/opm/models/io/vtktensorfunction.hh b/opm/models/io/vtktensorfunction.hh new file mode 100644 index 000000000..c99060533 --- /dev/null +++ b/opm/models/io/vtktensorfunction.hh @@ -0,0 +1,127 @@ +// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- +// vi: set et ts=4 sw=4 sts=4: +/* + This file is part of the Open Porous Media project (OPM). + + OPM is free software: you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation, either version 2 of the License, or + (at your option) any later version. + + OPM is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU General Public License + along with OPM. If not, see . + + Consult the COPYING file in the top-level source directory of this + module for the precise wording of the license and the list of + copyright holders. +*/ +/*! + * \file + * \copydoc Opm::VtkTensorFunction + */ +#ifndef VTK_TENSOR_FUNCTION_HH +#define VTK_TENSOR_FUNCTION_HH + +#include + +#include +#include +#include + +#include + +#include +#include +#include + +namespace Opm { + +/*! + * \brief Provides a tensor-valued function using Dune::FieldMatrix objects as elements. + */ +template +class VtkTensorFunction : public Dune::VTKFunction +{ + enum { dim = GridView::dimension }; + typedef typename GridView::ctype ctype; + typedef typename GridView::template Codim<0>::Entity Element; + + typedef BaseOutputWriter::TensorBuffer TensorBuffer; + +public: + VtkTensorFunction(std::string name, + const GridView& gridView, + const Mapper& mapper, + const TensorBuffer& buf, + unsigned codim, + unsigned matrixColumnIdx) + : name_(name) + , gridView_(gridView) + , mapper_(mapper) + , buf_(buf) + , codim_(codim) + , matrixColumnIdx_(matrixColumnIdx) + { assert(int(buf_.size()) == int(mapper_.size())); } + + virtual std::string name() const + { return name_; } + + virtual int ncomps() const + { return static_cast(buf_[0].M()); } + + virtual double evaluate(int mycomp, + const Element& e, + const Dune::FieldVector& xi) const + { + size_t idx; + if (codim_ == 0) { + // cells. map element to the index + idx = static_cast(mapper_.index(e)); + } + else if (codim_ == dim) { + // find vertex which is closest to xi in local + // coordinates. This code is based on Dune::P1VTKFunction + double min = 1e100; + int imin = -1; + Dune::GeometryType gt = e.type(); + int n = static_cast(e.subEntities(dim)); + for (int i = 0; i < n; ++i) { + Dune::FieldVector local = + Dune::ReferenceElements::general(gt).position(i, dim); + + local -= xi; + if (local.infinity_norm() < min) { + min = local.infinity_norm(); + imin = i; + } + } + + // map vertex to an index + idx = static_cast(mapper_.subIndex(e, imin, codim_)); + } + else + throw std::logic_error("Only element and vertex based tensor fields are supported so far."); + + unsigned i = static_cast(mycomp); + unsigned j = static_cast(matrixColumnIdx_); + + return static_cast(static_cast(buf_[idx][i][j])); + } + +private: + const std::string name_; + const GridView gridView_; + const Mapper& mapper_; + const TensorBuffer& buf_; + unsigned codim_; + unsigned matrixColumnIdx_; +}; + +} // namespace Opm + +#endif diff --git a/opm/models/io/vtkvectorfunction.hh b/opm/models/io/vtkvectorfunction.hh new file mode 100644 index 000000000..4fa29448a --- /dev/null +++ b/opm/models/io/vtkvectorfunction.hh @@ -0,0 +1,124 @@ +// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- +// vi: set et ts=4 sw=4 sts=4: +/* + This file is part of the Open Porous Media project (OPM). + + OPM is free software: you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation, either version 2 of the License, or + (at your option) any later version. + + OPM is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU General Public License + along with OPM. If not, see . + + Consult the COPYING file in the top-level source directory of this + module for the precise wording of the license and the list of + copyright holders. +*/ +/*! + * \file + * \copydoc Opm::VtkVectorFunction + */ +#ifndef VTK_VECTOR_FUNCTION_HH +#define VTK_VECTOR_FUNCTION_HH + +#include + +#include +#include +#include +#include + +#include + +#include +#include +#include + +namespace Opm { + +/*! + * \brief Provides a vector-valued function using Dune::FieldVectors + * as elements. + */ +template +class VtkVectorFunction : public Dune::VTKFunction +{ + enum { dim = GridView::dimension }; + typedef typename GridView::ctype ctype; + typedef typename GridView::template Codim<0>::Entity Element; + + typedef BaseOutputWriter::VectorBuffer VectorBuffer; + +public: + VtkVectorFunction(std::string name, + const GridView& gridView, + const Mapper& mapper, + const VectorBuffer& buf, + unsigned codim) + : name_(name) + , gridView_(gridView) + , mapper_(mapper) + , buf_(buf) + , codim_(codim) + { assert(int(buf_.size()) == int(mapper_.size())); } + + virtual std::string name() const + { return name_; } + + virtual int ncomps() const + { return static_cast(buf_[0].size()); } + + virtual double evaluate(int mycomp, + const Element& e, + const Dune::FieldVector& xi) const + { + unsigned idx; + if (codim_ == 0) { + // cells. map element to the index + idx = static_cast(mapper_.index(e)); + } + else if (codim_ == dim) { + // find vertex which is closest to xi in local + // coordinates. This code is based on Dune::P1VTKFunction + double min = 1e100; + int imin = -1; + Dune::GeometryType gt = e.type(); + int n = static_cast(e.subEntities(dim)); + for (int i = 0; i < n; ++i) { + Dune::FieldVector local = + Dune::ReferenceElements::general(gt).position(i, dim); + + local -= xi; + if (local.infinity_norm() < min) { + min = local.infinity_norm(); + imin = i; + } + } + + // map vertex to an index + idx = static_cast(mapper_.subIndex(e, imin, codim_)); + } + else + throw std::logic_error("Only element and vertex based vector fields are " + "supported so far."); + + return static_cast(static_cast(buf_[idx][static_cast(mycomp)])); + } + +private: + const std::string name_; + const GridView gridView_; + const Mapper& mapper_; + const VectorBuffer& buf_; + unsigned codim_; +}; + +} // namespace Opm + +#endif diff --git a/opm/models/utils/basicproperties.hh b/opm/models/utils/basicproperties.hh index ab760162a..2eb64d5a9 100644 --- a/opm/models/utils/basicproperties.hh +++ b/opm/models/utils/basicproperties.hh @@ -32,7 +32,7 @@ #include #include -#include +#include #if HAVE_DUNE_FEM #include diff --git a/opm/models/utils/simulator.hh b/opm/models/utils/simulator.hh index a2ffe073c..edd5d8132 100644 --- a/opm/models/utils/simulator.hh +++ b/opm/models/utils/simulator.hh @@ -28,7 +28,7 @@ #ifndef EWOMS_SIMULATOR_HH #define EWOMS_SIMULATOR_HH -#include +#include #include #include