3.3.0
DUNE for Multi-{Phase, Component, Scale, Physics, ...} flow and transport in porous media
staggerednewtonconvergencewriter.hh
Go to the documentation of this file.
1// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
2// vi: set et ts=4 sw=4 sts=4:
3/*****************************************************************************
4 * See the file COPYING for full copying permissions. *
5 * *
6 * This program is free software: you can redistribute it and/or modify *
7 * it under the terms of the GNU General Public License as published by *
8 * the Free Software Foundation, either version 3 of the License, or *
9 * (at your option) any later version. *
10 * *
11 * This program is distributed in the hope that it will be useful, *
12 * but WITHOUT ANY WARRANTY; without even the implied warranty of *
13 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
14 * GNU General Public License for more details. *
15 * *
16 * You should have received a copy of the GNU General Public License *
17 * along with this program. If not, see <http://www.gnu.org/licenses/>. *
18 *****************************************************************************/
26#ifndef DUMUX_STAGGERED_NEWTON_CONVERGENCE_WRITER_HH
27#define DUMUX_STAGGERED_NEWTON_CONVERGENCE_WRITER_HH
28
29#include <string>
30
31#include <dune/grid/io/file/vtk/vtksequencewriter.hh>
35
36namespace Dumux {
37
46template <class GridGeometry, class SolutionVector>
48{
49 using GridView = typename GridGeometry::GridView;
50
51 using CellCenterSolutionVector = typename std::decay_t<decltype(std::declval<SolutionVector>()[GridGeometry::cellCenterIdx()])>;
52 using FaceSolutionVector = typename std::decay_t<decltype(std::declval<SolutionVector>()[GridGeometry::faceIdx()])>;
53
54 using Scalar = typename CellCenterSolutionVector::block_type::value_type;
55
56 static constexpr auto numEqCellCenter = CellCenterSolutionVector::block_type::dimension;
57 static constexpr auto numEqFace = FaceSolutionVector::block_type::dimension;
58
59 using Element = typename GridView::template Codim<0>::Entity;
60 using GlobalPosition = typename Element::Geometry::GlobalCoordinate;
61
62 static_assert(GridGeometry::discMethod == DiscretizationMethod::staggered,
63 "This convergence writer does only work for the staggered method, use the NewtonConvergenceWriter instead");
64public:
70 StaggeredNewtonConvergenceWriter(const GridGeometry& gridGeometry,
71 const std::string& name = "newton_convergence")
72 : gridGeometry_(gridGeometry)
73 , ccWriter_(gridGeometry.gridView(), name, "", "")
74 , faceWriter_(std::make_shared<PointCloudVtkWriter<Scalar, GlobalPosition>>(coordinates_))
75 , faceSequenceWriter_(faceWriter_, name + "-face", "","",
76 gridGeometry.gridView().comm().rank(),
77 gridGeometry.gridView().comm().size())
78 {
79 resize();
80
81 for (int eqIdx = 0; eqIdx < numEqCellCenter; ++eqIdx)
82 {
83 ccWriter_.addCellData(xCellCenter_[eqIdx], "x_" + std::to_string(eqIdx));
84 ccWriter_.addCellData(deltaCellCenter_[eqIdx], "delta_" + std::to_string(eqIdx));
85 ccWriter_.addCellData(defCellCenter_[eqIdx], "defect_" + std::to_string(eqIdx));
86 }
87 }
88
90 void resize()
91 {
92 const auto numCellCenterDofs = gridGeometry_.numCellCenterDofs();
93 const auto numFaceDofs = gridGeometry_.numFaceDofs();
94
95 // resize the cell center output fields
96 for (int eqIdx = 0; eqIdx < numEqCellCenter; ++eqIdx)
97 {
98 defCellCenter_[eqIdx].resize(numCellCenterDofs);
99 deltaCellCenter_[eqIdx].resize(numCellCenterDofs);
100 xCellCenter_[eqIdx].resize(numCellCenterDofs);
101 }
102
103 // resize the face output fields
104 for (int eqIdx = 0; eqIdx < numEqFace; ++eqIdx)
105 {
106 defFace_[eqIdx].resize(numFaceDofs);
107 deltaFace_[eqIdx].resize(numFaceDofs);
108 xFace_[eqIdx].resize(numFaceDofs);
109 }
110
111 coordinates_.resize(numFaceDofs);
112 for (auto&& facet : facets(gridGeometry_.gridView()))
113 {
114 const auto dofIdxGlobal = gridGeometry_.gridView().indexSet().index(facet);
115 coordinates_[dofIdxGlobal] = facet.geometry().center();
116 }
117 }
118
121 void reset(std::size_t newId = 0UL)
122 { id_ = newId; iteration_ = 0UL; }
123
124 void write(const SolutionVector& uLastIter,
125 const SolutionVector& deltaU,
126 const SolutionVector& residual) override
127 {
128 assert(uLastIter.size() == deltaU.size() && uLastIter.size() == residual.size());
129
130 for (std::size_t dofIdxGlobal = 0; dofIdxGlobal < deltaU[GridGeometry::cellCenterIdx()].size(); ++dofIdxGlobal)
131 {
132 for (int eqIdx = 0; eqIdx < numEqCellCenter; ++eqIdx)
133 {
134 xCellCenter_[eqIdx][dofIdxGlobal] = uLastIter[GridGeometry::cellCenterIdx()][dofIdxGlobal][eqIdx];
135 deltaCellCenter_[eqIdx][dofIdxGlobal] = - deltaU[GridGeometry::cellCenterIdx()][dofIdxGlobal][eqIdx];
136 defCellCenter_[eqIdx][dofIdxGlobal] = residual[GridGeometry::cellCenterIdx()][dofIdxGlobal][eqIdx];
137 }
138 }
139
140 for (int eqIdx = 0; eqIdx < numEqFace; ++eqIdx)
141 {
142 faceWriter_->addPointData(xFace_[eqIdx], "x_" + std::to_string(eqIdx));
143 faceWriter_->addPointData(deltaFace_[eqIdx], "delta_" + std::to_string(eqIdx));
144 faceWriter_->addPointData(defFace_[eqIdx], "defect_" + std::to_string(eqIdx));
145 }
146
147 for (std::size_t dofIdxGlobal = 0; dofIdxGlobal < deltaU[GridGeometry::faceIdx()].size(); ++dofIdxGlobal)
148 {
149 for (int eqIdx = 0; eqIdx < numEqFace; ++eqIdx)
150 {
151 xFace_[eqIdx][dofIdxGlobal] = uLastIter[GridGeometry::faceIdx()][dofIdxGlobal][eqIdx];
152 deltaFace_[eqIdx][dofIdxGlobal] = - deltaU[GridGeometry::faceIdx()][dofIdxGlobal][eqIdx];
153 defFace_[eqIdx][dofIdxGlobal] = residual[GridGeometry::faceIdx()][dofIdxGlobal][eqIdx];
154 }
155 }
156
157 ccWriter_.write(static_cast<double>(id_) + static_cast<double>(iteration_)/1000);
158 faceSequenceWriter_.write(static_cast<double>(id_) + static_cast<double>(iteration_)/1000);
159 ++iteration_;
160 }
161
162private:
163 std::size_t id_ = 0UL;
164 std::size_t iteration_ = 0UL;
165
166 const GridGeometry& gridGeometry_;
167
168 Dune::VTKSequenceWriter<GridView> ccWriter_;
169
170 std::vector<GlobalPosition> coordinates_;
171 std::shared_ptr<PointCloudVtkWriter<Scalar, GlobalPosition>> faceWriter_;
173
174 std::array<std::vector<Scalar>, numEqCellCenter> defCellCenter_;
175 std::array<std::vector<Scalar>, numEqCellCenter> deltaCellCenter_;
176 std::array<std::vector<Scalar>, numEqCellCenter> xCellCenter_;
177
178 std::array<std::vector<Scalar>, numEqFace> defFace_;
179 std::array<std::vector<Scalar>, numEqFace> deltaFace_;
180 std::array<std::vector<Scalar>, numEqFace> xFace_;
181};
182
183} // end namespace Dumux
184
185#endif
Base class to write pvd-files which contains a list of all collected vtk-files. This is a modified ve...
A VTK writer specialized for staggered grid implementations with dofs on the faces.
This class provides the infrastructure to write the convergence behaviour of the newton method into a...
Definition: adapt.hh:29
A VTK output module to simplify writing dumux simulation data to VTK format.
Definition: pointcloudvtkwriter.hh:51
Base class to write pvd-files which contains a list of all collected vtk-files. This is a modified ve...
Definition: vtksequencewriter.hh:40
Definition: newtonconvergencewriter.hh:39
Writes the intermediate solutions for every Newton iteration (for staggered grid scheme)
Definition: staggerednewtonconvergencewriter.hh:48
void reset(std::size_t newId=0UL)
Definition: staggerednewtonconvergencewriter.hh:121
void write(const SolutionVector &uLastIter, const SolutionVector &deltaU, const SolutionVector &residual) override
Definition: staggerednewtonconvergencewriter.hh:124
StaggeredNewtonConvergenceWriter(const GridGeometry &gridGeometry, const std::string &name="newton_convergence")
Constructor.
Definition: staggerednewtonconvergencewriter.hh:70
void resize()
Resizes the output fields. This has to be called whenever the grid changes.
Definition: staggerednewtonconvergencewriter.hh:90