diff --git a/dumux/freeflow/staggered/model.hh b/dumux/freeflow/staggered/model.hh index 2856c75b03d06a2853f5e16ff367bb6535e87c19..16dbdaef447df3f0fbcb38c5d3e5462bcc894dd2 100644 --- a/dumux/freeflow/staggered/model.hh +++ b/dumux/freeflow/staggered/model.hh @@ -93,116 +93,6 @@ public: // NonIsothermalModel::maybeAddTemperature(vtkOutputModule); } - - /*! - * \brief \copybrief Dumux::ImplicitModel::addOutputVtkFields - * - * Specialization for the NavierStokesModel, adding the pressure and - * the process rank to the VTK writer. - */ - template - void addOutputVtkFields(const SolutionVector &sol, - MultiWriter &writer) - { - // TODO: implement vtk output properly, account for 3d - using VectorField = Dune::BlockVector >; - - // create the required scalar fields - const auto numElements = this->gridView_().size(0); - auto *p = writer.allocateManagedBuffer(numElements); - auto *delP = writer.allocateManagedBuffer(numElements); - - auto *v_x_pos = writer.allocateManagedBuffer(numElements); - auto *v_x_neg = writer.allocateManagedBuffer(numElements); - auto *v_y_pos = writer.allocateManagedBuffer(numElements); - auto *v_y_neg = writer.allocateManagedBuffer(numElements); - - VectorField *velocity = writer.template allocateManagedBuffer(numElements); - - // initialize velocity field - for (unsigned int i = 0; i < numElements; ++i) - { - (*velocity)[i] = Scalar(0); - } - - auto *rank = writer.allocateManagedBuffer(numElements); - - for (const auto& element : elements(this->gridView_(), Dune::Partitions::interior)) - { - auto eIdx = this->elementMapper().index(element); - (*rank)[eIdx] = this->gridView_().comm().rank(); - - // get the local fv geometry - auto fvGeometry = localView(this->globalFvGeometry()); - fvGeometry.bindElement(element); - - auto elemVolVars = localView(this->curGlobalVolVars()); - elemVolVars.bindElement(element, fvGeometry, this->curSol()); - - for (auto&& scv : scvs(fvGeometry)) - { - const auto& volVars = elemVolVars[scv]; - auto dofIdxGlobal = scv.dofIndex(); - - (*p)[dofIdxGlobal] = volVars.pressure(); - (*delP)[dofIdxGlobal] = volVars.pressure() - 1.1e5; - - GlobalPosition velocityVector(0.0); - for (auto&& scvf : scvfs(fvGeometry)) - { - auto& origFaceVars = this->curGlobalFaceVars().faceVars(scvf.dofIndex()); - auto dirIdx = scvf.directionIndex(); - - velocityVector[dirIdx] += 0.5*origFaceVars.velocity(); - - if(scvf.unitOuterNormal()[dirIdx] > 0.0) - { - if(dirIdx == 0) - (*v_x_pos)[dofIdxGlobal] = origFaceVars.velocity(); - if(dirIdx == 1) - (*v_y_pos)[dofIdxGlobal] = origFaceVars.velocity(); - } - else - { - if(dirIdx == 0) - (*v_x_neg)[dofIdxGlobal] = origFaceVars.velocity(); - if(dirIdx == 1) - (*v_y_neg)[dofIdxGlobal] = origFaceVars.velocity(); - } - } - (*velocity)[dofIdxGlobal] = velocityVector; - } - } - writer.attachDofData(*p, "p", isBox); - writer.attachDofData(*delP, "delP", isBox); - writer.attachDofData(*v_x_pos, "v_x_pos", isBox); - writer.attachDofData(*v_x_neg, "v_x_neg", isBox); - writer.attachDofData(*v_y_pos, "v_y_pos", isBox); - writer.attachDofData(*v_y_neg, "v_y_neg", isBox); - writer.attachCellData(*rank, "process rank"); - writer.attachDofData(*velocity, "velocity", isBox, dim); - } - - auto velocity(const Element& element) const - { - GlobalPosition velocityVector(0.0); - - // get the local fv geometry - auto fvGeometry = localView(this->globalFvGeometry()); - fvGeometry.bindElement(element); - for (auto&& scv : scvs(fvGeometry)) - { - for (auto&& scvf : scvfs(fvGeometry)) - { - auto& origFaceVars = this->curGlobalFaceVars().faceVars(scvf.dofIndex()); - auto dirIdx = scvf.directionIndex(); - velocityVector[dirIdx] += 0.5*origFaceVars.velocity(); - } - } - return velocityVector; - } - - }; } diff --git a/dumux/freeflow/staggered/vtkoutputmodule.hh b/dumux/freeflow/staggered/vtkoutputmodule.hh index f3fe5f22d386bd8e268e1ee26902bda61882ce9b..3a55f7c8281ebed939b0c0e6e933131c7fdddf43 100644 --- a/dumux/freeflow/staggered/vtkoutputmodule.hh +++ b/dumux/freeflow/staggered/vtkoutputmodule.hh @@ -71,13 +71,13 @@ private: * \param face The face */ template - void getVectorData_(Data& priVarVectorData, const Face& face) + void getPrivarVectorData_(Data& priVarVectorData, const Face& face) { const int dofIdxGlobal = face.dofIndex(); const int dirIdx = directionIndex(face.unitOuterNormal()); const Scalar velocity = this->problem().model().curSol()[faceIdx][dofIdxGlobal][0]; - for (int i = 0; i < this->faceData().priVarVectorDataInfo.size(); ++i) - priVarVectorData[i][dofIdxGlobal * this->faceData().priVarVectorDataInfo[i].pvIdx.size() + dirIdx] = velocity; + for (int i = 0; i < this->priVarVectorDataInfo_.size(); ++i) + priVarVectorData[i][dofIdxGlobal * this->priVarVectorDataInfo_[i].pvIdx.size() + dirIdx] = velocity; } }; diff --git a/dumux/io/pointcloudvtkwriter.hh b/dumux/io/pointcloudvtkwriter.hh new file mode 100644 index 0000000000000000000000000000000000000000..767ca96ebd78b76d2efc50dd88b0d0425e97df0c --- /dev/null +++ b/dumux/io/pointcloudvtkwriter.hh @@ -0,0 +1,402 @@ +// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- +// vi: set et ts=4 sw=4 sts=4: +/***************************************************************************** + * See the file COPYING for full copying permissions. * + * * + * This program 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. * + * * + * This program 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 this program. If not, see . * + *****************************************************************************/ +/*! + * \file + * \brief A VTK writer specialized for staggered grid implementations with dofs on the faces + */ +#ifndef STAGGERED_VTK_WRITER_HH +#define STAGGERED_VTK_WRITER_HH + +#include + +#include +#include +#include +#include + + +namespace Dumux +{ + +/*! + * \ingroup InputOutput + * \brief A VTK output module to simplify writing dumux simulation data to VTK format + * + * Handles the output of scalar and vector fields to VTK formatted file for multiple + * variables and timesteps. Certain predefined fields can be registered on problem / model + * initialization and/or be turned on/off using the designated properties. Additionally + * non-standardized scalar and vector fields can be added to the writer manually. + */ +template +class PointCloudVtkWriter +{ + using GlobalPosition = Dune::FieldVector; + + static constexpr unsigned int precision = 6; + static constexpr unsigned int numBeforeLineBreak = 15; + + /*! + * \brief A class holding a data container and additional information + */ + template + class VTKFunction + { + public: + /*! + * \brief A class holding a data container and additional information + * + * \param data The data container + * \param name The name of the data set + * \param numComponents The number of components of the data set + */ + VTKFunction(const ContainerType& data, const std::string& name, const int numComponents) : data_(data), name_(name), numComponents_(numComponents) + {} + + /*! + * \brief Returns the name of the data set + */ + const std::string& name() const + { + return name_; + } + + /*! + * \brief Returns the number of components of the data set + */ + int numComponents() const + { + return numComponents_; + } + + /*! + * \brief Allows random acces to data + * + * \param dofIdx The index + */ + auto& operator() (const int idx) const { return data_[idx]; } + + decltype(auto) begin() const + { + return data_.begin(); + } + + decltype(auto) end() const + { + return data_.end(); + } + + /*! + * \brief Returns the size of the data container + */ + int size() const + { + return data_.size(); + } + + private: + const ContainerType& data_; + const std::string name_; + const int numComponents_; + }; + + +public: + using ScalarFunction = VTKFunction>; + using VectorFunction = VTKFunction>; + + + PointCloudVtkWriter(const std::vector& coordinates) : coordinates_(coordinates) + {} + + /*! + * \brief Create an output file + * + * \param name The base name + * \param type The output type + */ + void write(const std::string& name, Dune::VTK::OutputType type = Dune::VTK::ascii) + { + auto filename = getSerialPieceName(name, ""); + + std::ofstream file; + file.open(filename); + writeHeader_(file); + writeCoordinates_(file, coordinates_); + writeDataInfo_(file); + + for(auto&& data : scalarPointData_) + { + writeData_(file, data); + } + for(auto&& data :vectorPointData_) + { + writeData_(file, data); + } + + file << "\n"; + file << "\n"; + file << "\n"; + file << ""; + + clear(); + + file.close(); + } + + /*! + * \brief Create an output file in parallel + * + * \param name The base name + * \param type The output type + */ + void pwrite(const std::string & name, const std::string & path, const std::string & extendpath, + Dune::VTK::OutputType type = Dune::VTK::ascii) + { + DUNE_THROW(Dune::NotImplemented, "parallel point cloud vtk output not supported yet"); + } + + /*! + * \brief Add a vector of scalar data that live on arbitrary points to the visualization. + * + * \param v The vector containing the data + * \param name The name of the data set + * \param ncomps The number of components of the data set + */ + void addPointData(const std::vector& v, const std::string &name, int ncomps = 1) + { + assert(v.size() == ncomps * coordinates_.size()); + scalarPointData_.push_back(ScalarFunction(v, name, ncomps)); + } + + /*! + * \brief Add a vector of vector data that live on arbitrary points to the visualization. + * + * \param v The vector containing the data + * \param name The name of the data set + * \param ncomps The number of components of the data set + */ + void addPointData(const std::vector& v, const std::string &name, int ncomps = 1) + { + assert(v.size() == coordinates_.size()); + vectorPointData_.push_back(VectorFunction(v, name, ncomps)); + } + + /*! + * \brief Clears the data lists + */ + void clear() + { + scalarPointData_.clear(); + vectorPointData_.clear(); + } + + /*! + * \brief Return name of a serial header file + * + * \param name Base name of the VTK output. This should be without + * any directory parts and without a filename extension. + * \param path Directory part of the resulting header name. May be + * empty, in which case the resulting name will not have a + * directory part. If non-empty, may or may not have a + * trailing '/'. If a trailing slash is missing, one is + * appended implicitly. + */ + std::string getSerialPieceName(const std::string& name, + const std::string& path) const + { + static const std::string extension = ".vtp"; + + return Dune::concatPaths(path, name+extension); + } + + /*! + * \brief Return name of a parallel header file + * + * \param name Base name of the VTK output. This should be without + * any directory parts and without a filename extension. + * \param path Directory part of the resulting header name. May be + * empty, in which case the resulting name will not have a + * directory part. If non-empty, may or may not have a + * trailing '/'. If a trailing slash is missing, one is + * appended implicitly. + * \param commSize Number of processes writing a parallel vtk output. + */ + std::string getParallelHeaderName(const std::string& name, + const std::string& path, + int commSize) const + { + std::ostringstream s; + if(path.size() > 0) { + s << path; + if(path[path.size()-1] != '/') + s << '/'; + } + s << 's' << std::setw(4) << std::setfill('0') << commSize << '-'; + s << name; + s << ".pvtp"; + return s.str(); + } + + +private: + /*! + * \brief Writes the header to the file + */ + void writeHeader_(std::ostream& file) + { + std::string header = "\n"; + header += "\n"; + header += "\n"; + header += "\n"; + file << header; + } + + /*! + * \brief Writes information about the data to the file + */ + void writeDataInfo_(std::ostream& file) + { + std::string scalarName; + std::string vectorName; + bool foundScalar = false; + bool foundVector = false; + + for(auto&& data : scalarPointData_) + { + if(data.numComponents() == 1 && !foundScalar) + { + scalarName = data.name(); + foundScalar = true; + continue; + } + + if(data.numComponents() > 1 && !foundVector) + { + vectorName = data.name(); + foundVector = true; + } + } + + for(auto&& data : vectorPointData_) + { + if(data.numComponents() > 1 && !foundVector) + { + vectorName = data.name(); + foundVector = true; + } + } + + if(foundScalar) + if(foundVector) + file << "\n"; + else + file << "\n"; + else if(foundVector) + file << "\n"; + else + return; + } + + /*! + * \brief Writes the coordinates to the file + * + * \param positions Container to store the positions + */ + void writeCoordinates_(std::ostream& file, const std::vector& positions) + { + // write the positions to the file + file << "\n"; + file << "\n"; + int counter = 0; + for(auto&& x : positions) + { + file << x ; + + if(x.size() == 1) + file << " 0 0 "; + if(x.size() == 2) + file << " 0 "; + + // introduce a line break after a certain time + if((++counter) > numBeforeLineBreak) + { + file << std::endl; + counter = 0; + } + } + file << "\n\n"; + file << "\n"; + } + + /*! + * \brief Writes data to the file + * + * \param data The data container which hold the data itself, as well as the name of the data set and the number of its components + */ + template + void writeData_(std::ostream& file, const T& data) + { + file << "\n"; + int counter = 0; + for(auto&& value : data) + { + // forward to specialized function + writeToFile_(file, value); + + // introduce a line break after a certain time + if((++counter) > numBeforeLineBreak) + { + file << std::endl; + counter = 0; + } + } + file << "\n\n"; + } + + /*! + * \brief Writes a scalar to the file + * + * \param s The scalar + */ + void writeToFile_(std::ostream& file, const Scalar& s) + { + file << s << " "; + } + + /*! + * \brief Writes a vector to the file + * + * \param g The vector + */ + void writeToFile_(std::ostream& file, const GlobalPosition& g) + { + assert(g.size() > 1 && g.size() < 4); + if(g.size() < 3) + file << g << " 0 "; + else + file << g; + } + + const std::vector& coordinates_; + std::list scalarPointData_; + std::list vectorPointData_; +}; +} // end namespace Dumux + +#endif diff --git a/dumux/io/staggeredvtkoutputmodule.hh b/dumux/io/staggeredvtkoutputmodule.hh index 4a6ab7c82ca44d4dc6a2675a44ea7a703e22bcf2..194701d5cfc6a74f41a4450e3f5bd6412386d107 100644 --- a/dumux/io/staggeredvtkoutputmodule.hh +++ b/dumux/io/staggeredvtkoutputmodule.hh @@ -26,7 +26,8 @@ #include #include -#include +#include +#include namespace Properties { @@ -55,6 +56,8 @@ class StaggeredVtkOutputModule : public VtkOutputModuleBase enum { dim = GridView::dimension }; enum { dimWorld = GridView::dimensionworld }; + using GlobalPosition = Dune::FieldVector; + using DofTypeIndices = typename GET_PROP(TypeTag, DofTypeIndices); typename DofTypeIndices::CellCenterIdx cellCenterIdx; typename DofTypeIndices::FaceIdx faceIdx; @@ -65,31 +68,18 @@ class StaggeredVtkOutputModule : public VtkOutputModuleBase using Positions = std::vector; using Data = std::vector>; - // a collection of types and variables to be passed to the staggered vtk writer - struct WriterData - { - using Positions = StaggeredVtkOutputModule::Positions; - using PriVarScalarDataInfo = StaggeredVtkOutputModule::PriVarScalarDataInfo; - using PriVarVectorDataInfo = StaggeredVtkOutputModule::PriVarVectorDataInfo; - using Data = StaggeredVtkOutputModule::Data; - Data priVarScalarData; - Data priVarVectorData; - Data secondVarScalarData; - Data secondVarVectorData; - Positions positions; - std::vector priVarScalarDataInfo; - std::vector priVarVectorDataInfo; - std::vector secondVarScalarDataInfo; - std::vector secondVarVectorDataInfo; - }; - public: StaggeredVtkOutputModule(const Problem& problem, - Dune::VTK::DataMode dm = Dune::VTK::conforming) : ParentType(problem, dm), faceWriter_(problem) + Dune::VTK::DataMode dm = Dune::VTK::conforming) : ParentType(problem, dm), + faceWriter_(std::make_shared>(coordinates_)), + sequenceWriter_(faceWriter_, problem.name() + "-face", "","", + problem.gridView().comm().rank(), + problem.gridView().comm().size() ) { writeFaceVars_ = GET_PARAM_FROM_GROUP(TypeTag, bool, Vtk, WriteFaceData); + coordinatesInitialized_ = false; } ////////////////////////////////////////////////////////////////////////////////////////////// @@ -97,12 +87,31 @@ public: //! Do not call these methods after initialization ////////////////////////////////////////////////////////////////////////////////////////////// + + //! Output a scalar field + //! \param name The name of the vtk field + //! \returns A reference to the resized scalar field to be filled with the actual data + std::vector& createFaceScalarField(const std::string& name) + { + faceScalarFields_.emplace_back(std::make_pair(std::vector(this->problem().model().numFaceDofs()), name)); + return faceScalarFields_.back().first; + } + + //! Output a vector field + //! \param name The name of the vtk field + //! \returns A reference to the resized vector field to be filled with the actual data + std::vector& createFaceVectorField(const std::string& name) + { + faceVectorFields_.emplace_back(std::make_pair(std::vector(this->problem().model().numFaceDofs()), name)); + return faceVectorFields_.back().first; + } + //! Output a scalar primary variable //! \param name The name of the vtk field //! \param pvIdx The index in the primary variables vector void addFacePrimaryVariable(const std::string& name, unsigned int pvIdx) { - faceData_.priVarScalarDataInfo.push_back(PriVarScalarDataInfo{pvIdx, name}); + priVarScalarDataInfo_.push_back(PriVarScalarDataInfo{pvIdx, name}); } //! Output a vector primary variable @@ -111,7 +120,7 @@ public: void addFacePrimaryVariable(const std::string& name, std::vector pvIndices) { assert(pvIndices.size() < 4 && "Vtk doesn't support vector dimensions greater than 3!"); - faceData_.priVarVectorDataInfo.push_back(PriVarVectorDataInfo{pvIndices, name}); + priVarVectorDataInfo_.push_back(PriVarVectorDataInfo{pvIndices, name}); } void write(double time, Dune::VTK::OutputType type = Dune::VTK::ascii) @@ -142,16 +151,25 @@ protected: return this->problem().model().curSol()[cellCenterIdx][dofIdxGlobal][pvIdx]; } - /*! - * \brief Returns a reference to the face data - */ - const WriterData& faceData() const - { - return faceData_; - } + std::vector priVarScalarDataInfo_; + std::vector priVarVectorDataInfo_; + std::vector secondVarScalarDataInfo_; + std::vector secondVarVectorDataInfo_; private: + void updateCoordinates_() + { + std::cout << "updating coordinates" << std::endl; + coordinates_.resize(this->problem().model().numFaceDofs()); + for(auto&& facet : facets(this->problem().gridView())) + { + const int dofIdxGlobal = this->problem().gridView().indexSet().index(facet); + coordinates_[dofIdxGlobal] = facet.geometry().center(); + } + coordinatesInitialized_ = true; + } + /*! * \brief Gathers all face-related data and invokes the face vtk-writer using these data. */ @@ -163,13 +181,14 @@ private: std::vector dofVisited(numPoints, false); // get fields for all primary coordinates and variables - Positions positions(numPoints*3); + if(!coordinatesInitialized_) + updateCoordinates_(); - Data priVarScalarData(faceData_.priVarScalarDataInfo.size(), std::vector(numPoints)); + Data priVarScalarData(priVarScalarDataInfo_.size(), std::vector(numPoints)); - Data priVarVectorData(faceData_.priVarVectorDataInfo.size()); - for (std::size_t i = 0; i < faceData_.priVarVectorDataInfo.size(); ++i) - priVarVectorData[i].resize(numPoints*faceData_.priVarVectorDataInfo[i].pvIdx.size()); + Data priVarVectorData(priVarVectorDataInfo_.size()); + for (std::size_t i = 0; i < priVarVectorDataInfo_.size(); ++i) + priVarVectorData[i].resize(numPoints*priVarVectorDataInfo_[i].pvIdx.size()); for(auto&& element : elements(this->problem().gridView())) { @@ -180,54 +199,35 @@ private: if(dofVisited[scvf.dofIndex()]) continue; - asImp_().getPositions_(positions, scvf); - asImp_().getScalarData_(priVarScalarData, scvf); - asImp_().getVectorData_(priVarVectorData, scvf); + asImp_().getPrivarScalarData_(priVarScalarData, scvf); + asImp_().getPrivarVectorData_(priVarVectorData, scvf); + dofVisited[scvf.dofIndex()] = true; } } - faceData_.priVarScalarData = std::move(priVarScalarData); - faceData_.priVarVectorData = std::move(priVarVectorData); - faceData_.positions = std::move(positions); -// results.secondVarScalarData = xxx; //TODO: implemented secondVarData -// results.secondVarScalarData = xxx; + // transfer priVar scalar data to writer + for(int i = 0; i < priVarScalarDataInfo_.size(); ++i) + faceWriter_->addPointData(priVarScalarData[i], priVarScalarDataInfo_[i].name); - faceWriter_.write(faceData_); - } + // transfer priVar vector data to writer + for(int i = 0; i < priVarVectorDataInfo_.size(); ++i) + faceWriter_->addPointData(priVarVectorData[i], priVarVectorDataInfo_[i].name, priVarVectorDataInfo_[i].pvIdx.size()); - /*! - * \brief Retrives the position (center) of the face. - * ParaView expects three-dimensional coordinates, therefore we fill empty entries with zero for dim < 3. - * - * \param positions Container to store the positions - * \param face The face - */ - template - void getPositions_(Positions& positions, const Face& face) - { - const int dofIdxGlobal = face.dofIndex(); - auto pos = face.center(); - if(dim == 1) - { - positions[dofIdxGlobal*3] = pos[0]; - positions[dofIdxGlobal*3 + 1] = 0.0; - positions[dofIdxGlobal*3 + 2] = 0.0; - } - else if(dim == 2) - { - positions[dofIdxGlobal*3] = pos[0]; - positions[dofIdxGlobal*3 + 1] = pos[1]; - positions[dofIdxGlobal*3 + 2] = 0.0; - } - else - { - positions[dofIdxGlobal*3] = pos[0]; - positions[dofIdxGlobal*3 + 1] = pos[1] ; - positions[dofIdxGlobal*3 + 2] = pos[2] ; - } + // transfer custom scalar data to writer + for(auto&& scalarField : faceScalarFields_) + faceWriter_->addPointData(scalarField.first, scalarField.second); + + // transfer custom vector data to writer + for(auto&& vectorField : faceVectorFields_) + faceWriter_->addPointData(vectorField.first, vectorField.second, 3); + + sequenceWriter_.write(this->problem().timeManager().time() + 1.0); + faceScalarFields_.clear(); + faceVectorFields_.clear(); } + /*! * \brief Retrives scalar-valued data from the face. * @@ -235,13 +235,11 @@ private: * \param face The face */ template - void getScalarData_(Data& priVarScalarData, const Face& face) + void getPrivarScalarData_(Data& priVarScalarData, const Face& face) { const int dofIdxGlobal = face.dofIndex(); - for(int pvIdx = 0; pvIdx < faceData_.priVarScalarDataInfo.size(); ++pvIdx) - { + for(int pvIdx = 0; pvIdx < priVarScalarDataInfo_.size(); ++pvIdx) priVarScalarData[pvIdx][dofIdxGlobal] = this->problem().model().curSol()[faceIdx][dofIdxGlobal][pvIdx]; - } } /*! @@ -251,19 +249,27 @@ private: * \param face The face */ template - void getVectorData_(Data& priVarVectorData, const Face& face) + void getPrivarVectorData_(Data& priVarVectorData, const Face& face) { const int dofIdxGlobal = face.dofIndex(); - for (int i = 0; i < faceData_.priVarVectorDataInfo.size(); ++i) - for (int j = 0; j < faceData_.priVarVectorDataInfo[i].pvIdx.size(); ++j) - priVarVectorData[i][dofIdxGlobal*faceData_.priVarVectorDataInfo[i].pvIdx.size() + j] - = this->problem().model().curSol()[faceIdx][dofIdxGlobal][faceData_.priVarVectorDataInfo[i].pvIdx[j]]; + for (int i = 0; i < priVarVectorDataInfo_.size(); ++i) + for (int j = 0; j < priVarVectorDataInfo_[i].pvIdx.size(); ++j) + priVarVectorData[i][dofIdxGlobal*priVarVectorDataInfo_[i].pvIdx.size() + j] + = this->problem().model().curSol()[faceIdx][dofIdxGlobal][priVarVectorDataInfo_[i].pvIdx[j]]; } - StaggeredVtkWriter faceWriter_; - WriterData faceData_; + std::shared_ptr> faceWriter_; + + VTKSequenceWriter> sequenceWriter_; + bool writeFaceVars_; + std::vector coordinates_; + bool coordinatesInitialized_; + + std::list, std::string>> faceScalarFields_; + std::list, std::string>> faceVectorFields_; + //! Returns the implementation of the problem (i.e. static polymorphism) Implementation &asImp_() { return *static_cast(this); } diff --git a/dumux/io/staggeredvtkwriter.hh b/dumux/io/staggeredvtkwriter.hh deleted file mode 100644 index 8db5bbbc59d3fe512c14c8899923f2d4db66b87b..0000000000000000000000000000000000000000 --- a/dumux/io/staggeredvtkwriter.hh +++ /dev/null @@ -1,267 +0,0 @@ -// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- -// vi: set et ts=4 sw=4 sts=4: -/***************************************************************************** - * See the file COPYING for full copying permissions. * - * * - * This program 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. * - * * - * This program 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 this program. If not, see . * - *****************************************************************************/ -/*! - * \file - * \brief A VTK writer specialized for staggered grid implementations with dofs on the faces - */ -#ifndef STAGGERED_VTK_WRITER_HH -#define STAGGERED_VTK_WRITER_HH - -#include - -#include -#include - -namespace Properties -{ -NEW_PROP_TAG(VtkAddVelocity); -NEW_PROP_TAG(VtkAddProcessRank); -} - -namespace Dumux -{ - -/*! - * \ingroup InputOutput - * \brief A VTK output module to simplify writing dumux simulation data to VTK format - * - * Handles the output of scalar and vector fields to VTK formatted file for multiple - * variables and timesteps. Certain predefined fields can be registered on problem / model - * initialization and/or be turned on/off using the designated properties. Additionally - * non-standardized scalar and vector fields can be added to the writer manually. - */ -template -class StaggeredVtkWriter -{ - using Problem = typename GET_PROP_TYPE(TypeTag, Problem); - using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar); - using GridView = typename GET_PROP_TYPE(TypeTag, GridView); - using DofTypeIndices = typename GET_PROP(TypeTag, DofTypeIndices); - typename DofTypeIndices::CellCenterIdx cellCenterIdx; - typename DofTypeIndices::FaceIdx faceIdx; - - enum { dim = GridView::dimension }; - enum { dimWorld = GridView::dimensionworld }; - using GlobalPosition = Dune::FieldVector; - - static constexpr unsigned int precision = 6; - static constexpr unsigned int numBeforeLineBreak = 15; - - using ScalarInfo = std::vector; - using VectorInfo = std::vector; - using Positions = typename WriterData::Positions; - using Data = typename WriterData::Data; - -public: - StaggeredVtkWriter(const Problem& problem) : problem_(problem) - {} - - void write(const WriterData& data) - { - if(data.priVarScalarDataInfo.empty() && - data.priVarVectorDataInfo.empty() && - data.secondVarScalarDataInfo.empty() && - data.secondVarVectorDataInfo.empty()) - return; - - file_.open(fileName_()); - write_(data); - file_.close(); - ++curWriterNum_; - } - -private: - void writeHeader_(const int numPoints) - { - std::string header = "\n"; - header += "\n"; - header += "\n"; - header += "\n"; - file_ << header; - } - - /*! - * \brief Writes the coordinates and the simulation data to the file - * - * \param priVarScalarData Container to store the scalar-valued data - * \param priVarVectorData Container to store the vector-valued data - * \param scalarInfo Information concerning the data (e.g. names) - * \param vectorInfo Information concerning the data (e.g. names) - * \param positions Container to store the positions - */ - void write_(const WriterData& data) - { - const int numPoints = problem_.model().numFaceDofs(); - writeHeader_(numPoints); - writeCoordinates_(data.positions); - - std::string scalarName; - std::string vectorName; - - bool scalarValuesPresent = false; - bool vectorValuesPresent = false; - - - if(!(data.priVarScalarDataInfo.empty() && data.secondVarScalarDataInfo.empty())) - { - scalarValuesPresent = true; - scalarName = data.priVarScalarDataInfo.empty() ? data.secondVarScalarDataInfo[0].name : - data.priVarScalarDataInfo[0].name; - } - if(!(data.priVarVectorDataInfo.empty() && data.secondVarVectorDataInfo.empty())) - { - vectorValuesPresent = true; - vectorName = data.priVarVectorDataInfo.empty() ? data.secondVarVectorDataInfo[0].name : - data.priVarVectorDataInfo[0].name; - } - - if(scalarValuesPresent) - if(vectorValuesPresent) - file_ << ""; - else - file_ << ""; - else if(vectorValuesPresent) - file_ << ""; - else - return; - - file_ << std::endl; - writeAllData_(data); - - file_ << "\n"; - file_ << "\n"; - file_ << "\n"; - file_ << ""; - } - - /*! - * \brief Writes the coordinates to the file - * - * \param positions Container to store the positions - */ - void writeCoordinates_(const Positions& positions) - { - // write the positions to the file - file_ << "\n"; - file_ << "\n"; - int counter = 0; - for(auto&& x : positions) - { - file_ << x << " "; - // introduce a line break after a certain time - if((++counter) > numBeforeLineBreak) - { - file_ << std::endl; - counter = 0; - } - } - file_ << "\n\n"; - file_ << "\n"; - } - - void writeAllData_(const WriterData& data) - { - // write scalar-valued primary variable data - if(!data.priVarScalarDataInfo.empty()) - writeScalarData_(data.priVarScalarData, data.priVarScalarDataInfo); - - // write vector-valued primary variable data - if(!data.priVarVectorDataInfo.empty()) - writeVectorData_(data.priVarVectorData, data.priVarVectorDataInfo); - - // write scalar-valued secondary variable data - if(!data.secondVarScalarDataInfo.empty()) - writeScalarData_(data.secondVarScalarData, data.secondVarScalarDataInfo); - - // write vector-valued primary variable data - if(!data.secondVarVectorDataInfo.empty()) - writeVectorData_(data.secondVarVectorData, data.secondVarVectorDataInfo); - } - - /*! - * \brief Writes scalar-valued data to the file - * - * \param priVarScalarData Container to store the data - * \param info Information concerning the data (e.g. names) - */ - void writeScalarData_(const Data& priVarScalarData, const ScalarInfo& info) - { - // write the priVars to the file - for(int pvIdx = 0; pvIdx < info.size(); ++pvIdx) - { - file_ << "\n"; - int counter = 0; - for(auto&& x : priVarScalarData[pvIdx]) - { - file_ << x << " " ; - // introduce a line break after a certain time - if((++counter) > numBeforeLineBreak) - { - file_ << std::endl; - counter = 0; - } - } - file_ << "\n\n"; - } - } - - /*! - * \brief Writes vector-valued data to the file - * - * \param priVarVectorData Container to store the data - * \param info Information concerning the data (e.g. names) - */ - void writeVectorData_(const Data& priVarVectorData, const VectorInfo& info) - { - // write the priVars to the file - for(int i = 0; i < info.size(); ++i) - { - file_ << "\n"; - int counter = 0; - for(auto&& x : priVarVectorData[i]) - { - file_ << x << " " ; - // introduce a line break after a certain time - if((++counter) > numBeforeLineBreak) - { - file_ << std::endl; - counter = 0; - } - } - file_ << "\n\n"; - } - } - - /*! - * \brief Returns the file name for each timestep - */ - std::string fileName_() const - { - std::ostringstream oss; - oss << problem_.name() << "-face-" << std::setw(5) << std::setfill('0') << curWriterNum_ << ".vtp"; - return oss.str(); - } - - const Problem& problem_; - std::ofstream file_; - int curWriterNum_{0}; -}; -} // end namespace Dumux - -#endif diff --git a/dumux/io/vtksequencewriter.hh b/dumux/io/vtksequencewriter.hh new file mode 100644 index 0000000000000000000000000000000000000000..3ee08aaec44e7b6dd22ad64d5a1f76ec68b2e2f1 --- /dev/null +++ b/dumux/io/vtksequencewriter.hh @@ -0,0 +1,125 @@ +// -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*- +// vi: set et ts=4 sw=2 sts=2: + +#ifndef DUMUX_VTKSEQUENCEWRITER_HH +#define DUMUX_VTKSEQUENCEWRITER_HH + +#include +#include +#include +#include +#include +#include + +#include +#include + + +namespace Dumux { + + /** \brief Base class to write pvd-files which contains a list of all collected vtk-files. + * This is a modified version of DUNE's pvd writer which takes a VTKWriter as template + * argument making it more general. + * + * Derive from this class to write pvd-file suitable for easy visualization with + * The Visualization Toolkit (VTK). + * + * \tparam VTKWriter The VTKWriter class + * + */ + + template + class VTKSequenceWriter + { + std::shared_ptr vtkWriter_; + std::vector timesteps_; + std::string name_,path_,extendpath_; + int rank_; + int size_; + public: + /** \brief Set up the VTKSequenceWriter class + * + * \param vtkWriter Writer object used to write the individual time step data files + * \param rank Process number in a multi-process setting + * \param size Total number of processes + */ + explicit VTKSequenceWriter( std::shared_ptr vtkWriter, + const std::string& name, + const std::string& path, + const std::string& extendpath, + int rank, + int size) + : vtkWriter_(vtkWriter), + name_(name), path_(path), + extendpath_(extendpath), + rank_(rank), + size_(size) + {} + + ~VTKSequenceWriter() {} + + + /** + * \brief Writes VTK data for the given time, + * \param time The time(step) for the data to be written. + * \param type VTK output type. + */ + void write (double time, Dune::VTK::OutputType type = Dune::VTK::ascii) + { + /* remember current time step */ + unsigned int count = timesteps_.size(); + timesteps_.push_back(time); + + /* write VTK file */ + if(size_==1) + vtkWriter_->write(Dune::concatPaths(path_,seqName(count)),type); + else + vtkWriter_->pwrite(seqName(count), path_,extendpath_,type); + + /* write pvd file ... only on rank 0 */ + if (rank_==0) { + std::ofstream pvdFile; + pvdFile.exceptions(std::ios_base::badbit | std::ios_base::failbit | + std::ios_base::eofbit); + std::string pvdname = name_ + ".pvd"; + pvdFile.open(pvdname.c_str()); + pvdFile << " \n" + << " \n" + << " \n"; + for (unsigned int i=0; i<=count; i++) + { + // filename + std::string piecepath; + std::string fullname; + if(size_==1) { + piecepath = path_; + fullname = vtkWriter_->getSerialPieceName(seqName(i), piecepath); + } + else { + piecepath = Dune::concatPaths(path_, extendpath_); + fullname = vtkWriter_->getParallelHeaderName(seqName(i), piecepath, size_); + } + pvdFile << " \n"; + } + pvdFile << " \n" + << " \n" << std::flush; + pvdFile.close(); + } + } + private: + + // create sequence name + std::string seqName(unsigned int count) const + { + std::stringstream n; + n.fill('0'); + n << name_ << "-" << std::setw(5) << count; + return n.str(); + } + }; + +} // end namespace Dumux + +#endif diff --git a/test/freeflow/staggered/doneatestproblem.hh b/test/freeflow/staggered/doneatestproblem.hh index 47ecd10d034afb742d426c10262fce182d0bd148..cfdc785122b75ed02ce7753249de7157d9487698 100644 --- a/test/freeflow/staggered/doneatestproblem.hh +++ b/test/freeflow/staggered/doneatestproblem.hh @@ -34,9 +34,6 @@ #include -#include - - namespace Dumux { template @@ -291,51 +288,44 @@ public: } /*! - * \brief Append all quantities of interest which can be derived - * from the solution of the current time step to the VTK - * writer. + * \brief Adds additional VTK output data to the VTKWriter. Function is called by the output module on every write. */ - void addOutputVtkFields() + template + void addVtkOutputFields(VtkOutputModule& outputModule) const { - //Here we calculate the analytical solution - const auto numElements = this->gridView().size(0); - auto& pExact = *(this->resultWriter().allocateManagedBuffer(numElements)); - auto& velocityExact = *(this->resultWriter()).template allocateManagedBuffer(numElements); - - using DofTypeIndices = typename GET_PROP(TypeTag, DofTypeIndices); - typename DofTypeIndices::FaceIdx faceIdx; - - const auto someElement = *(elements(this->gridView()).begin()); + auto& pressureExact = outputModule.createScalarField("pressureExact", 0); + auto& velocityExact = outputModule.createVectorField("velocityExact", 0); - auto someFvGeometry = localView(this->model().globalFvGeometry()); - someFvGeometry.bindElement(someElement); - const auto someScv = *(scvs(someFvGeometry).begin()); - - Scalar time = std::max(this->timeManager().time() + this->timeManager().timeStepSize(), 1e-10); + auto& scalarFaceVelocityExact = outputModule.createFaceScalarField("scalarFaceVelocityExact"); + auto& vectorFaceVelocityExact = outputModule.createFaceVectorField("vectorFaceVelocityExact"); for (const auto& element : elements(this->gridView())) { auto fvGeometry = localView(this->model().globalFvGeometry()); fvGeometry.bindElement(element); - for (auto&& scv : scvs(fvGeometry)) { - auto dofIdxGlobal = scv.dofIndex(); - const auto& globalPos = scv.dofPosition(); - - pExact[dofIdxGlobal] = dirichletAtPos(globalPos)[pressureIdx]; + auto ccDofIdx = scv.dofIndex(); + auto ccDofPosition = scv.dofPosition(); + auto analyticalSolutionAtCc = dirichletAtPos(ccDofPosition); GlobalPosition velocityVector(0.0); for (auto&& scvf : scvfs(fvGeometry)) { + auto faceDofIdx = scvf.dofIndex(); + auto faceDofPosition = scvf.center(); auto dirIdx = scvf.directionIndex(); - velocityVector[dirIdx] += 0.5*dirichletAtPos(globalPos)[faceIdx][dirIdx]; + auto analyticalSolutionAtFace = dirichletAtPos(faceDofPosition); + scalarFaceVelocityExact[faceDofIdx] = analyticalSolutionAtFace[faceIdx][dirIdx]; + + GlobalPosition tmp(0.0); + tmp[dirIdx] = analyticalSolutionAtFace[faceIdx][dirIdx]; + vectorFaceVelocityExact[faceDofIdx] = std::move(tmp); } - velocityExact[dofIdxGlobal] = velocityVector; + pressureExact[ccDofIdx] = analyticalSolutionAtCc[pressureIdx]; + velocityExact[ccDofIdx] = analyticalSolutionAtCc[faceIdx]; } } - this->resultWriter().attachDofData(pExact, "p_exact", false); - this->resultWriter().attachDofData(velocityExact, "velocity_exact", false, dim); } /*! diff --git a/test/freeflow/staggered/kovasznaytestproblem.hh b/test/freeflow/staggered/kovasznaytestproblem.hh index c274b3a1c2f63b05e6259397a1d5bb415cae594d..9fe0c3f2ccc8e97b13a9eacce82fd38c9b3bc8ff 100644 --- a/test/freeflow/staggered/kovasznaytestproblem.hh +++ b/test/freeflow/staggered/kovasznaytestproblem.hh @@ -31,9 +31,6 @@ #include #include - -#include - // solve Navier-Stokes equations #define ENABLE_NAVIERSTOKES 1 @@ -300,51 +297,44 @@ public: } /*! - * \brief Append all quantities of interest which can be derived - * from the solution of the current time step to the VTK - * writer. + * \brief Adds additional VTK output data to the VTKWriter. Function is called by the output module on every write. */ - void addOutputVtkFields() + template + void addVtkOutputFields(VtkOutputModule& outputModule) const { - //Here we calculate the analytical solution - const auto numElements = this->gridView().size(0); - auto& pExact = *(this->resultWriter().allocateManagedBuffer(numElements)); - auto& velocityExact = *(this->resultWriter()).template allocateManagedBuffer(numElements); - - using DofTypeIndices = typename GET_PROP(TypeTag, DofTypeIndices); - typename DofTypeIndices::FaceIdx faceIdx; - - const auto someElement = *(elements(this->gridView()).begin()); + auto& pressureExact = outputModule.createScalarField("pressureExact", 0); + auto& velocityExact = outputModule.createVectorField("velocityExact", 0); - auto someFvGeometry = localView(this->model().globalFvGeometry()); - someFvGeometry.bindElement(someElement); - const auto someScv = *(scvs(someFvGeometry).begin()); - - Scalar time = std::max(this->timeManager().time() + this->timeManager().timeStepSize(), 1e-10); + auto& scalarFaceVelocityExact = outputModule.createFaceScalarField("scalarFaceVelocityExact"); + auto& vectorFaceVelocityExact = outputModule.createFaceVectorField("vectorFaceVelocityExact"); for (const auto& element : elements(this->gridView())) { auto fvGeometry = localView(this->model().globalFvGeometry()); fvGeometry.bindElement(element); - for (auto&& scv : scvs(fvGeometry)) { - auto dofIdxGlobal = scv.dofIndex(); - const auto& globalPos = scv.dofPosition(); - - pExact[dofIdxGlobal] = dirichletAtPos(globalPos)[pressureIdx]; + auto ccDofIdx = scv.dofIndex(); + auto ccDofPosition = scv.dofPosition(); + auto analyticalSolutionAtCc = dirichletAtPos(ccDofPosition); GlobalPosition velocityVector(0.0); for (auto&& scvf : scvfs(fvGeometry)) { + auto faceDofIdx = scvf.dofIndex(); + auto faceDofPosition = scvf.center(); auto dirIdx = scvf.directionIndex(); - velocityVector[dirIdx] += 0.5*dirichletAtPos(globalPos)[faceIdx][dirIdx]; + auto analyticalSolutionAtFace = dirichletAtPos(faceDofPosition); + scalarFaceVelocityExact[faceDofIdx] = analyticalSolutionAtFace[faceIdx][dirIdx]; + + GlobalPosition tmp(0.0); + tmp[dirIdx] = analyticalSolutionAtFace[faceIdx][dirIdx]; + vectorFaceVelocityExact[faceDofIdx] = std::move(tmp); } - velocityExact[dofIdxGlobal] = velocityVector; + pressureExact[ccDofIdx] = analyticalSolutionAtCc[pressureIdx]; + velocityExact[ccDofIdx] = analyticalSolutionAtCc[faceIdx]; } } - this->resultWriter().attachDofData(pExact, "p_exact", false); - this->resultWriter().attachDofData(velocityExact, "velocity_exact", false, dim); } /*!