Commit 57a24c9d authored by Dennis Gläser's avatar Dennis Gläser Committed by Timo Koch
Browse files

[boxdfm] add output module for separate output of fracture solution

parent adf2f340
......@@ -235,6 +235,24 @@ public:
}
}
protected:
// some return functions for differing implementations to use
const Problem& problem() const { return problem_; }
const FVGridGeometry& fvGridGeometry() const { return gridGeom_; }
const GridVariables& gridVariables() const { return gridVariables_; }
const SolutionVector& sol() const { return sol_; }
const bool verbose() const { return verbose_; }
const std::string& name() const { return name_; }
const Dune::VTK::DataMode dataMode() const { return dm_; }
Dune::VTKWriter<GridView>& writer() { return *writer_; }
Dune::VTKSequenceWriter<GridView>& sequenceWriter() { return sequenceWriter_; }
const std::vector<VolVarScalarDataInfo>& volVarScalarDataInfo() const { return volVarScalarDataInfo_; }
const std::vector<VolVarVectorDataInfo>& volVarVectorDataInfo() const { return volVarVectorDataInfo_; }
const std::vector<Field>& fields() const { return fields_; }
private:
//! Assembles the fields and adds them to the writer (conforming output)
......
// -*- 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 <http://www.gnu.org/licenses/>. *
*****************************************************************************/
/*!
* \file
* \ingroup InputOutput
* \brief A VTK output module to simplify writing dumux simulation data to VTK format
*/
#ifndef POROUSMEDIUMFLOW_BOXDFM_VTK_OUTPUT_MODULE_HH
#define POROUSMEDIUMFLOW_BOXDFM_VTK_OUTPUT_MODULE_HH
#include <set>
#include <dune/grid/common/gridfactory.hh>
#include <dune/geometry/referenceelements.hh>
#include <dune/grid/common/mcmgmapper.hh>
#include <dumux/io/vtkoutputmodule.hh>
namespace Dumux {
/*!
* \ingroup InputOutput
* \brief A VTK output module to simplify writing dumux simulation data to VTK format.
* This output module is specialized for writing out data obtained by the box-dfm
* scheme. It writes out separate vtk files for the solution on the fracture. For
* this, a grid type to be used the fracture-conforming lower-dimensional grid has
* to be provided.
*
* \tparam TypeTag The TypeTag of the problem implementation
* \tparam FractureGrid The Type used for the lower-dimensional grid
* \tparam phaseIdxOffset Used for single-phase problems to retrieve the right phase name
*
* Handles the output of scalar and vector fields to VTK formatted file for multiple
* variables and timesteps. Certain predefined fields can be registered on
* 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 TypeTag, class FractureGrid, int phaseIdxOffset = 0>
class BoxDfmVtkOutputModule : public VtkOutputModule<TypeTag, phaseIdxOffset>
{
using ParentType = VtkOutputModule<TypeTag, phaseIdxOffset>;
using Problem = typename GET_PROP_TYPE(TypeTag, Problem);
using GridVariables = typename GET_PROP_TYPE(TypeTag, GridVariables);
using SolutionVector = typename GET_PROP_TYPE(TypeTag, SolutionVector);
using FVGridGeometry = typename GET_PROP_TYPE(TypeTag, FVGridGeometry);
using VelocityOutput = typename GET_PROP_TYPE(TypeTag, VelocityOutput);
using ModelTraits = typename GET_PROP_TYPE(TypeTag, ModelTraits);
static constexpr int numPhases = ModelTraits::numPhases();
using VV = typename GridVariables::VolumeVariables;
using FluidSystem = typename VV::FluidSystem;
using Scalar = typename GridVariables::Scalar;
using GridView = typename FVGridGeometry::GridView;
using FractureGridView = typename FractureGrid::LeafGridView;
using FractureVertexMapper = Dune::MultipleCodimMultipleGeomTypeMapper<FractureGridView>;
enum
{
dim = GridView::dimension,
dimWorld = GridView::dimensionworld
};
using IndexType = typename GridView::IndexSet::IndexType;
using Element = typename GridView::template Codim<0>::Entity;
using GlobalPosition = typename Element::Geometry::GlobalCoordinate;
using ReferenceElements = typename Dune::ReferenceElements<typename GridView::ctype, dim>;
using Field = Vtk::template Field<GridView>;
using FractureField = Vtk::template Field<FractureGridView>;
static_assert(dim > 1, "Box-Dfm output only works for dim > 1");
static_assert(FractureGrid::dimension == int(dim-1), "Fracture grid must be of codimension one!");
static_assert(FractureGrid::dimensionworld == int(dimWorld), "Fracture grid has to has the same coordinate dimension!");
static_assert(FVGridGeometry::discMethod == DiscretizationMethod::box, "Box-Dfm output module can only be used with the box scheme!");
public:
BoxDfmVtkOutputModule(const Problem& problem,
const FVGridGeometry& fvGridGeometry,
const GridVariables& gridVariables,
const SolutionVector& sol,
const std::string& name,
const std::string& paramGroup = "",
Dune::VTK::DataMode dm = Dune::VTK::conforming,
bool verbose = true)
: ParentType(problem, fvGridGeometry, gridVariables, sol, name, paramGroup, dm, verbose)
{
const auto& gridView = fvGridGeometry.gridView();
Dune::GridFactory<FractureGrid> gridFactory;
// insert fracture vertices
std::size_t fracVertexCount = 0;
vertexToFractureVertexIdx_.resize(gridView.size(dim));
for (const auto& v : vertices(gridView))
{
const auto vIdxGlobal = fvGridGeometry.vertexMapper().index(v);
if (fvGridGeometry.dofOnFracture(vIdxGlobal))
{
gridFactory.insertVertex(v.geometry().center());
vertexToFractureVertexIdx_[vIdxGlobal] = fracVertexCount++;
}
}
// insert fracture elements
std::set< std::pair<IndexType, unsigned int> > handledFacets;
auto vertIsOnFracture = [&fvGridGeometry] (auto idx) { return fvGridGeometry.dofOnFracture(idx); };
for (const auto& element : elements(gridView))
{
const auto referenceElement = ReferenceElements::general(element.geometry().type());
for (const auto& is : intersections(gridView, element))
{
// obtain all vertex indices on this intersection
const auto& isGeometry = is.geometry();
const auto numCorners = isGeometry.corners();
std::vector<IndexType> isVertexIndices(numCorners);
for (unsigned int i = 0; i < numCorners; ++i)
{
const auto vIdxLocal = referenceElement.subEntity(is.indexInInside(), 1, i, dim);
isVertexIndices[i] = fvGridGeometry.vertexMapper().subIndex(element, vIdxLocal, dim);
}
// determine if this is a fracture facet & if it has to be inserted
bool insertFacet = false;
if ( std::all_of(isVertexIndices.begin(), isVertexIndices.end(), vertIsOnFracture) )
{
insertFacet = true;
if (!is.boundary())
{
// only proceed if facet has not been handled yet
const auto pair = std::make_pair<IndexType, unsigned int>(fvGridGeometry.elementMapper().index(is.outside()),
is.indexInOutside());
if (handledFacets.count( pair ) != 0)
insertFacet = false;
}
}
if (insertFacet)
{
// transform intersection vertex indices to frac grid indices
std::for_each( isVertexIndices.begin(),
isVertexIndices.end(),
[&] (auto& idx) { idx = this->vertexToFractureVertexIdx_[idx]; } );
// insert the element
gridFactory.insertElement(isGeometry.type(), isVertexIndices);
// insert to set of handled facets
handledFacets.insert( std::make_pair<IndexType, unsigned int>(fvGridGeometry.elementMapper().index(element),
is.indexInInside()) );
}
}
}
auto grid = gridFactory.createGrid();
fractureGridView_ = std::make_unique<FractureGridView>(grid->leafGridView());
// update fracture vertex mapper
fractureVertexMapper_ = std::make_unique<FractureVertexMapper>(*fractureGridView_, Dune::mcmgVertexLayout());
// obtain map fracture vertex insertion index -> fracture grid vertex index
std::vector<IndexType> insToVertexIdx(fractureGridView_->size(FractureGridView::dimension));
for (const auto& v : vertices(*fractureGridView_))
insToVertexIdx[ gridFactory.insertionIndex(v) ] = fractureVertexMapper_->index(v);
// update the index map
for (IndexType dofIdx = 0; dofIdx < gridView.size(GridView::dimension); ++dofIdx)
if (fvGridGeometry.dofOnFracture(dofIdx))
vertexToFractureVertexIdx_[dofIdx] = insToVertexIdx[ vertexToFractureVertexIdx_[dofIdx] ];
// instantiate writers for the fracture
fractureWriter_ = std::make_shared< Dune::VTKWriter<FractureGridView> >(*fractureGridView_, dm);
fractureSequenceWriter_ = std::make_unique< Dune::VTKSequenceWriter<FractureGridView> >(fractureWriter_, name + "_fracture");
}
//////////////////////////////////////////////////////////////////////////////////////////////
//! Writing data
//////////////////////////////////////////////////////////////////////////////////////////////
//! Write the data for this timestep to file in four steps
//! (1) We assemble all registered variable fields
//! (2) We register them with the vtk writer
//! (3) The writer writes the output for us
//! (4) Clear the writer for the next time step
void write(double time, Dune::VTK::OutputType type = Dune::VTK::ascii)
{
Dune::Timer timer;
// write to file depending on data mode
const auto dm = this->dataMode();
if (dm == Dune::VTK::conforming)
writeConforming_(time, type);
else if (dm == Dune::VTK::nonconforming)
writeNonConforming_(time, type);
else
DUNE_THROW(Dune::NotImplemented, "Output for provided vtk data mode");
//! output
timer.stop();
if (this->verbose())
std::cout << "Writing output for problem \"" << this->name() << "\". Took " << timer.elapsed() << " seconds." << std::endl;
}
private:
//! Assembles the fields and adds them to the writer (conforming output)
void writeConforming_(double time, Dune::VTK::OutputType type)
{
//////////////////////////////////////////////////////////////
//! (1) Assemble all variable fields and add to writer
//////////////////////////////////////////////////////////////
// instatiate the velocity output
VelocityOutput velocityOutput(this->problem(), this->fvGridGeometry(), this->gridVariables(), this->sol());
std::array<std::vector<GlobalPosition>, numPhases> velocity;
// process rank
static bool addProcessRank = getParamFromGroup<bool>(this->paramGroup(), "Vtk.AddProcessRank");
std::vector<double> rank;
// volume variable data
std::vector<std::vector<Scalar>> volVarScalarData;
std::vector<std::vector<Scalar>> volVarScalarDataFracture;
std::vector<std::vector<GlobalPosition>> volVarVectorData;
std::vector<std::vector<GlobalPosition>> volVarVectorDataFracture;
// some references for convenience
const auto& gridView = this->fvGridGeometry().gridView();
const auto& volVarScalarDataInfo = this->volVarScalarDataInfo();
const auto& volVarVectorDataInfo = this->volVarVectorDataInfo();
//! Abort if no data was registered
if (!volVarScalarDataInfo.empty()
|| !volVarVectorDataInfo.empty()
|| !this->fields().empty()
|| velocityOutput.enableOutput()
|| addProcessRank)
{
const auto numCells = gridView.size(0);
const auto numDofs = gridView.size(dim);
const auto numFractureVert = fractureGridView_->size(FractureGridView::dimension);
// get fields for all volume variables
if (!this->volVarScalarDataInfo().empty())
{
volVarScalarData.resize(volVarScalarDataInfo.size(), std::vector<Scalar>(numDofs));
volVarScalarDataFracture.resize(volVarScalarDataInfo.size(), std::vector<Scalar>(numFractureVert));
}
if (!this->volVarVectorDataInfo().empty())
{
volVarVectorData.resize(volVarVectorDataInfo.size(), std::vector<GlobalPosition>(numDofs));
volVarVectorDataFracture.resize(volVarVectorDataInfo.size(), std::vector<GlobalPosition>(numFractureVert));
}
if (velocityOutput.enableOutput())
for (int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx)
velocity[phaseIdx].resize(numDofs);
// maybe allocate space for the process rank
if (addProcessRank) rank.resize(numCells);
for (const auto& element : elements(gridView, Dune::Partitions::interior))
{
const auto eIdxGlobal = this->fvGridGeometry().elementMapper().index(element);
auto fvGeometry = localView(this->fvGridGeometry());
auto elemVolVars = localView(this->gridVariables().curGridVolVars());
// If velocity output is enabled we need to bind to the whole stencil
// otherwise element-local data is sufficient
if (velocityOutput.enableOutput())
{
fvGeometry.bind(element);
elemVolVars.bind(element, fvGeometry, this->sol());
}
else
{
fvGeometry.bindElement(element);
elemVolVars.bindElement(element, fvGeometry, this->sol());
}
if (!volVarScalarDataInfo.empty() || !volVarVectorDataInfo.empty())
{
for (auto&& scv : scvs(fvGeometry))
{
const auto dofIdxGlobal = scv.dofIndex();
const auto& volVars = elemVolVars[scv];
// get the scalar-valued data
if (!scv.isOnFracture())
{
for (std::size_t i = 0; i < volVarScalarDataInfo.size(); ++i)
volVarScalarData[i][dofIdxGlobal] = volVarScalarDataInfo[i].get(volVars);
}
else
{
for (std::size_t i = 0; i < volVarScalarDataInfo.size(); ++i)
volVarScalarDataFracture[i][vertexToFractureVertexIdx_[dofIdxGlobal]] = volVarScalarDataInfo[i].get(volVars);
}
// get the vector-valued data
if (!scv.isOnFracture())
{
for (std::size_t i = 0; i < volVarVectorDataInfo.size(); ++i)
volVarVectorData[i][dofIdxGlobal] = volVarVectorDataInfo[i].get(volVars);
}
else
{
for (std::size_t i = 0; i < volVarVectorDataInfo.size(); ++i)
volVarVectorDataFracture[i][vertexToFractureVertexIdx_[dofIdxGlobal]] = volVarVectorDataInfo[i].get(volVars);
}
}
}
// velocity output
if (velocityOutput.enableOutput())
for (int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx)
velocityOutput.calculateVelocity(velocity[phaseIdx], elemVolVars, fvGeometry, element, phaseIdx);
//! the rank
if (addProcessRank)
rank[eIdxGlobal] = static_cast<double>(gridView.comm().rank());
}
//////////////////////////////////////////////////////////////
//! (2) Register data fields with the vtk writer
//////////////////////////////////////////////////////////////
// volume variables if any
for (std::size_t i = 0; i < volVarScalarDataInfo.size(); ++i)
{
this->sequenceWriter().addVertexData(volVarScalarData[i], volVarScalarDataInfo[i].name);
fractureSequenceWriter_->addVertexData(volVarScalarDataFracture[i], volVarScalarDataInfo[i].name);
}
for (std::size_t i = 0; i < volVarVectorDataInfo.size(); ++i)
{
this->sequenceWriter().addVertexData( Field(gridView, this->fvGridGeometry().vertexMapper(), volVarVectorData[i],
volVarVectorDataInfo[i].name, /*numComp*/dimWorld, /*codim*/dim).get() );
fractureSequenceWriter_->addVertexData( FractureField(*fractureGridView_, *fractureVertexMapper_, volVarVectorDataFracture[i],
volVarVectorDataInfo[i].name, /*numComp*/dimWorld, /*codim*/dim-1).get() );
}
// the velocity field
if (velocityOutput.enableOutput())
{
for (int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx)
this->sequenceWriter().addVertexData( Field(gridView, this->fvGridGeometry().vertexMapper(), velocity[phaseIdx],
"velocity_" + std::string(FluidSystem::phaseName(phaseIdx+phaseIdxOffset)) + " (m/s)",
/*numComp*/dimWorld, /*codim*/dim).get() );
}
// the process rank
if (addProcessRank)
this->sequenceWriter().addCellData( Field(gridView, this->fvGridGeometry().elementMapper(), rank,
"process rank", /*numComp*/1, /*codim*/0).get() );
// also register additional (non-standardized) user fields if any (only on matrix grid)
for (auto&& field : this->fields())
{
if (field.codim() == 0)
this->sequenceWriter().addCellData(field.get());
else if (field.codim() == dim)
this->sequenceWriter().addVertexData(field.get());
else
DUNE_THROW(Dune::RangeError, "Cannot add wrongly sized vtk scalar field!");
}
}
//////////////////////////////////////////////////////////////
//! (2) The writer writes the output for us
//////////////////////////////////////////////////////////////
this->sequenceWriter().write(time, type);
fractureSequenceWriter_->write(time, type);
//////////////////////////////////////////////////////////////
//! (3) Clear the writer
//////////////////////////////////////////////////////////////
this->writer().clear();
fractureWriter_->clear();
}
//! Assembles the fields and adds them to the writer (conforming output)
void writeNonConforming_(double time, Dune::VTK::OutputType type)
{
DUNE_THROW(Dune::NotImplemented, "Non-conforming output for Box-Dfm model");
}
std::unique_ptr<FractureVertexMapper> fractureVertexMapper_;
std::unique_ptr<FractureGridView> fractureGridView_;
std::shared_ptr<Dune::VTKWriter<FractureGridView>> fractureWriter_;
std::unique_ptr< Dune::VTKSequenceWriter<FractureGridView> > fractureSequenceWriter_;
std::vector<IndexType> vertexToFractureVertexIdx_;
};
} // end namespace Dumux
#endif
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment