diff --git a/CHANGELOG.md b/CHANGELOG.md
index 16ed7dbf9b1dbb7eedd2e2c11dae6989ea743e03..693fd71034413418929ced54b807b40e8b6993f3 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -7,13 +7,26 @@ Differences Between DuMuX 3.1 and DuMuX 3.0
- Added new porous medium model for the energy balance of a porous solid (general heat equation)
- __Multidomain__: It is now possible to use the facet coupling module together with the mpfa-o scheme in the bulk domain.
+- Added a `StaggeredNewtonConvergenceWriter` for the staggered grid discretization scheme
### Immediate interface changes not allowing/requiring a deprecation period
+- `NewtonConvergenceWriter`'s first template argument has changed from `GridView`, the `FVGridGeometry`. This allows to call the `resize()` method after a grid change without any arguments.
+ Here is an example of how to instatiate the convergence writer:
+ ```
+ using NewtonConvergenceWriter = Dumux::NewtonConvergenceWriter<FVGridGeometry, SolutionVector>;
+ auto convergenceWriter = std::make_shared<NewtonConvergenceWriter>(*fvGridGeometry);
+ ```
+
+
### Deprecated classes/files, to be removed after 3.1:
### Deprecated member functions, to be removed after 3.1:
+- The convergence writer is no longer passed to `NewtonSolver`'s `solve()` method.
+ For outputting convergence data, please use `newtonSolver.attachConvergenceWriter(convWriter)` in `main.cc` (directly after instantiating the writer).
+ To stop the output, the writer can also be detached again using `newtonSolver.detachConvergenceWriter()`.
+
### Deleted classes/files, property names, constants/enums
diff --git a/dumux/nonlinear/newtonconvergencewriter.hh b/dumux/nonlinear/newtonconvergencewriter.hh
index 0ebfcb6577246a00e018c3e2781105cf0e769c58..6d06350ed547ac3f5c3b8c1e11fe6aa3f37dfb4e 100644
--- a/dumux/nonlinear/newtonconvergencewriter.hh
+++ b/dumux/nonlinear/newtonconvergencewriter.hh
@@ -25,9 +25,10 @@
#ifndef DUMUX_NEWTON_CONVERGENCE_WRITER_HH
#define DUMUX_NEWTON_CONVERGENCE_WRITER_HH
-#include <dune/common/exceptions.hh>
+#include <string>
+
#include <dune/grid/io/file/vtk/vtksequencewriter.hh>
-#include <dumux/common/parameters.hh>
+#include <dumux/discretization/method.hh>
namespace Dumux {
@@ -49,12 +50,15 @@ struct ConvergenceWriterInterface
* to write out multiple Newton solves with a unique id, if you don't call use all
* Newton iterations just come after each other in the pvd file.
*/
-// template <class Scalar, class GridView, int numEq>
-template <class GridView, class SolutionVector>
+template <class FVGridGeometry, class SolutionVector>
class NewtonConvergenceWriter : public ConvergenceWriterInterface<SolutionVector>
{
+ using GridView = typename FVGridGeometry::GridView;
static constexpr auto numEq = SolutionVector::block_type::dimension;
using Scalar = typename SolutionVector::block_type::value_type;
+
+ static_assert(FVGridGeometry::discMethod != DiscretizationMethod::staggered,
+ "This convergence writer does not work for the staggered method, use the StaggeredNewtonConvergenceWriter instead");
public:
/*!
* \brief Constructor
@@ -62,14 +66,14 @@ public:
* \param size the size of the solution data
* \param name base name of the vtk output
*/
- NewtonConvergenceWriter(const GridView& gridView,
- std::size_t size,
+ NewtonConvergenceWriter(const FVGridGeometry& fvGridGeometry,
const std::string& name = "newton_convergence")
- : writer_(gridView, name, "", "")
+ : fvGridGeometry_(fvGridGeometry)
+ , writer_(fvGridGeometry.gridView(), name, "", "")
{
- resize(size);
+ resize();
- if (size == gridView.size(GridView::dimension))
+ if (FVGridGeometry::discMethod == DiscretizationMethod::box)
{
for (int eqIdx = 0; eqIdx < numEq; ++eqIdx)
{
@@ -78,7 +82,7 @@ public:
writer_.addVertexData(def_[eqIdx], "defect_" + std::to_string(eqIdx));
}
}
- else if (size == gridView.size(0))
+ else
{
for (int eqIdx = 0; eqIdx < numEq; ++eqIdx)
{
@@ -87,21 +91,19 @@ public:
writer_.addCellData(def_[eqIdx], "defect_" + std::to_string(eqIdx));
}
}
- else
- {
- DUNE_THROW(Dune::InvalidStateException, "Wrong size of output fields in the Newton convergence writer!");
- }
}
//! Resizes the output fields. This has to be called whenever the grid changes
- void resize(std::size_t size)
+ void resize()
{
+ const auto numDofs = fvGridGeometry_.numDofs();
+
// resize the output fields
for (int eqIdx = 0; eqIdx < numEq; ++eqIdx)
{
- def_[eqIdx].resize(size);
- delta_[eqIdx].resize(size);
- x_[eqIdx].resize(size);
+ def_[eqIdx].resize(numDofs);
+ delta_[eqIdx].resize(numDofs);
+ x_[eqIdx].resize(numDofs);
}
}
@@ -110,9 +112,9 @@ public:
void reset(std::size_t newId = 0UL)
{ id_ = newId; iteration_ = 0UL; }
- void write(const SolutionVector &uLastIter,
- const SolutionVector &deltaU,
- const SolutionVector &residual) override
+ void write(const SolutionVector& uLastIter,
+ const SolutionVector& deltaU,
+ const SolutionVector& residual) override
{
assert(uLastIter.size() == deltaU.size() && uLastIter.size() == residual.size());
@@ -134,6 +136,8 @@ private:
std::size_t id_ = 0UL;
std::size_t iteration_ = 0UL;
+ const FVGridGeometry& fvGridGeometry_;
+
Dune::VTKSequenceWriter<GridView> writer_;
std::array<std::vector<Scalar>, numEq> def_;
diff --git a/dumux/nonlinear/newtonsolver.hh b/dumux/nonlinear/newtonsolver.hh
index 9739eebd94e8f03d89937934916ffc7ba6ec2b81..4eeb779abffd974f15ae0a6a86cfd34d44125c49 100644
--- a/dumux/nonlinear/newtonsolver.hh
+++ b/dumux/nonlinear/newtonsolver.hh
@@ -49,6 +49,8 @@
#include "newtonconvergencewriter.hh"
+#include <dune/common/deprecated.hh>
+
namespace Dumux {
namespace Detail {
@@ -193,8 +195,21 @@ public:
* Does time step control when the Newton fails to converge
*/
template<class TimeLoop>
+ DUNE_DEPRECATED_MSG("Use attachConvergenceWriter(convWriter) and solve(x, *timeLoop) instead")
void solve(SolutionVector& uCurrentIter, TimeLoop& timeLoop,
- std::shared_ptr<ConvergenceWriter> convWriter = nullptr)
+ std::shared_ptr<ConvergenceWriter> convWriter)
+ {
+ if (!convergenceWriter_)
+ attachConvergenceWriter(convWriter);
+
+ solve(uCurrentIter, timeLoop);
+ }
+ /*!
+ * \brief Run the Newton method to solve a non-linear system.
+ * Does time step control when the Newton fails to converge
+ */
+ template<class TimeLoop>
+ void solve(SolutionVector& uCurrentIter, TimeLoop& timeLoop)
{
if (assembler_->isStationaryProblem())
DUNE_THROW(Dune::InvalidStateException, "Using time step control with stationary problem makes no sense!");
@@ -203,7 +218,7 @@ public:
for (std::size_t i = 0; i <= maxTimeStepDivisions_; ++i)
{
// linearize & solve
- const bool converged = solve_(uCurrentIter, convWriter);
+ const bool converged = solve_(uCurrentIter);
if (converged)
return;
@@ -238,9 +253,9 @@ public:
* \brief Run the Newton method to solve a non-linear system.
* The solver is responsible for all the strategic decisions.
*/
- void solve(SolutionVector& uCurrentIter, std::shared_ptr<ConvergenceWriter> convWriter = nullptr)
+ void solve(SolutionVector& uCurrentIter)
{
- const bool converged = solve_(uCurrentIter, convWriter);
+ const bool converged = solve_(uCurrentIter);
if (!converged)
DUNE_THROW(NumericalProblem, "Newton solver didn't converge after "
<< numSteps_ << " iterations.\n");
@@ -256,6 +271,15 @@ public:
{
numSteps_ = 0;
initPriVarSwitch_(u, HasPriVarsSwitch{});
+
+ // write the initial residual if a convergence writer was set
+ if (convergenceWriter_)
+ {
+ assembler_->assembleResidual(u);
+ SolutionVector delta(u);
+ delta = 0; // dummy vector, there is no delta before solving the linear system
+ convergenceWriter_->write(u, delta, assembler_->residual());
+ }
}
/*!
@@ -621,6 +645,18 @@ public:
const std::string& paramGroup() const
{ return paramGroup_; }
+ /*!
+ * \brief Attach a convergence writer to write out intermediate results after each iteration
+ */
+ void attachConvergenceWriter(std::shared_ptr<ConvergenceWriter> convWriter)
+ { convergenceWriter_ = convWriter; }
+
+ /*!
+ * \brief Detach the convergence writer to stop the output
+ */
+ void detachConvergenceWriter()
+ { convergenceWriter_ = nullptr; }
+
protected:
void computeResidualReduction_(const SolutionVector &uCurrentIter)
@@ -724,7 +760,7 @@ private:
* \brief Run the Newton method to solve a non-linear system.
* The solver is responsible for all the strategic decisions.
*/
- bool solve_(SolutionVector& uCurrentIter, std::shared_ptr<ConvergenceWriter> convWriter = nullptr)
+ bool solve_(SolutionVector& uCurrentIter)
{
// the given solution is the initial guess
SolutionVector uLastIter(uCurrentIter);
@@ -814,10 +850,10 @@ private:
newtonEndStep(uCurrentIter, uLastIter);
// if a convergence writer was specified compute residual and write output
- if (convWriter)
+ if (convergenceWriter_)
{
assembler_->assembleResidual(uCurrentIter);
- convWriter->write(uLastIter, deltaU, assembler_->residual());
+ convergenceWriter_->write(uCurrentIter, deltaU, assembler_->residual());
}
// detect if the method has converged
@@ -1240,6 +1276,9 @@ private:
std::unique_ptr<PrimaryVariableSwitch> priVarSwitch_;
//! if we switched primary variables in the last iteration
bool priVarsSwitchedInLastIteration_ = false;
+
+ //! convergence writer
+ std::shared_ptr<ConvergenceWriter> convergenceWriter_ = nullptr;
};
} // end namespace Dumux
diff --git a/dumux/nonlinear/staggerednewtonconvergencewriter.hh b/dumux/nonlinear/staggerednewtonconvergencewriter.hh
new file mode 100644
index 0000000000000000000000000000000000000000..1cb85da200eb9be4485731a8dbe9d380b29e785c
--- /dev/null
+++ b/dumux/nonlinear/staggerednewtonconvergencewriter.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:
+/*****************************************************************************
+ * 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 3 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 Nonlinear
+ * \brief This class provides the infrastructure to write the
+ * convergence behaviour of the newton method for the
+ * staggered discretization scheme into a VTK file.
+ */
+#ifndef DUMUX_STAGGERED_NEWTON_CONVERGENCE_WRITER_HH
+#define DUMUX_STAGGERED_NEWTON_CONVERGENCE_WRITER_HH
+
+#include <string>
+
+#include <dune/grid/io/file/vtk/vtksequencewriter.hh>
+#include <dumux/io/vtksequencewriter.hh>
+#include <dumux/io/pointcloudvtkwriter.hh>
+#include "newtonconvergencewriter.hh"
+
+namespace Dumux {
+
+/*!
+ * \ingroup Nonlinear
+ * \brief Writes the intermediate solutions for every Newton iteration (for staggered grid scheme)
+ * \note To use this create a shared_ptr to an instance of this class in the main file
+ * and pass it to newton.solve(x, convergencewriter). You can use the reset method
+ * to write out multiple Newton solves with a unique id, if you don't call use all
+ * Newton iterations just come after each other in the pvd file.
+ */
+template <class FVGridGeometry, class SolutionVector>
+class StaggeredNewtonConvergenceWriter : public ConvergenceWriterInterface<SolutionVector>
+{
+ using GridView = typename FVGridGeometry::GridView;
+
+ using CellCenterSolutionVector = typename std::decay_t<decltype(std::declval<SolutionVector>()[FVGridGeometry::cellCenterIdx()])>;
+ using FaceSolutionVector = typename std::decay_t<decltype(std::declval<SolutionVector>()[FVGridGeometry::faceIdx()])>;
+
+ using Scalar = typename CellCenterSolutionVector::block_type::value_type;
+
+ static constexpr auto numEqCellCenter = CellCenterSolutionVector::block_type::dimension;
+ static constexpr auto numEqFace = FaceSolutionVector::block_type::dimension;
+
+ using Element = typename GridView::template Codim<0>::Entity;
+ using GlobalPosition = typename Element::Geometry::GlobalCoordinate;
+
+ static_assert(FVGridGeometry::discMethod == DiscretizationMethod::staggered,
+ "This convergence writer does only work for the staggered method, use the NewtonConvergenceWriter instead");
+public:
+ /*!
+ * \brief Constructor
+ * \param fvGridGeometry The finite volume geometry on the grid view
+ * \param name Base name of the vtk output
+ */
+ StaggeredNewtonConvergenceWriter(const FVGridGeometry& fvGridGeometry,
+ const std::string& name = "newton_convergence")
+ : fvGridGeometry_(fvGridGeometry)
+ , ccWriter_(fvGridGeometry.gridView(), name, "", "")
+ , faceWriter_(std::make_shared<PointCloudVtkWriter<Scalar, GlobalPosition>>(coordinates_))
+ , faceSequenceWriter_(faceWriter_, name + "-face", "","",
+ fvGridGeometry.gridView().comm().rank(),
+ fvGridGeometry.gridView().comm().size())
+ {
+ resize();
+
+ for (int eqIdx = 0; eqIdx < numEqCellCenter; ++eqIdx)
+ {
+ ccWriter_.addCellData(xCellCenter_[eqIdx], "x_" + std::to_string(eqIdx));
+ ccWriter_.addCellData(deltaCellCenter_[eqIdx], "delta_" + std::to_string(eqIdx));
+ ccWriter_.addCellData(defCellCenter_[eqIdx], "defect_" + std::to_string(eqIdx));
+ }
+ }
+
+ //! Resizes the output fields. This has to be called whenever the grid changes
+ void resize()
+ {
+ const auto numCellCenterDofs = fvGridGeometry_.numCellCenterDofs();
+ const auto numFaceDofs = fvGridGeometry_.numFaceDofs();
+
+ // resize the cell center output fields
+ for (int eqIdx = 0; eqIdx < numEqCellCenter; ++eqIdx)
+ {
+ defCellCenter_[eqIdx].resize(numCellCenterDofs);
+ deltaCellCenter_[eqIdx].resize(numCellCenterDofs);
+ xCellCenter_[eqIdx].resize(numCellCenterDofs);
+ }
+
+ // resize the face output fields
+ for (int eqIdx = 0; eqIdx < numEqFace; ++eqIdx)
+ {
+ defFace_[eqIdx].resize(numFaceDofs);
+ deltaFace_[eqIdx].resize(numFaceDofs);
+ xFace_[eqIdx].resize(numFaceDofs);
+ }
+
+ coordinates_.resize(numFaceDofs);
+ for (auto&& facet : facets(fvGridGeometry_.gridView()))
+ {
+ const auto dofIdxGlobal = fvGridGeometry_.gridView().indexSet().index(facet);
+ coordinates_[dofIdxGlobal] = facet.geometry().center();
+ }
+ }
+
+ //! Reset the convergence writer for a possible next Newton step
+ //! You may set a different id in case you don't want the output to be overwritten by the next step
+ void reset(std::size_t newId = 0UL)
+ { id_ = newId; iteration_ = 0UL; }
+
+ void write(const SolutionVector& uLastIter,
+ const SolutionVector& deltaU,
+ const SolutionVector& residual) override
+ {
+ assert(uLastIter.size() == deltaU.size() && uLastIter.size() == residual.size());
+
+ for (std::size_t dofIdxGlobal = 0; dofIdxGlobal < deltaU[FVGridGeometry::cellCenterIdx()].size(); ++dofIdxGlobal)
+ {
+ for (int eqIdx = 0; eqIdx < numEqCellCenter; ++eqIdx)
+ {
+ xCellCenter_[eqIdx][dofIdxGlobal] = uLastIter[FVGridGeometry::cellCenterIdx()][dofIdxGlobal][eqIdx];
+ deltaCellCenter_[eqIdx][dofIdxGlobal] = - deltaU[FVGridGeometry::cellCenterIdx()][dofIdxGlobal][eqIdx];
+ defCellCenter_[eqIdx][dofIdxGlobal] = residual[FVGridGeometry::cellCenterIdx()][dofIdxGlobal][eqIdx];
+ }
+ }
+
+ for (int eqIdx = 0; eqIdx < numEqFace; ++eqIdx)
+ {
+ faceWriter_->addPointData(xFace_[eqIdx], "x_" + std::to_string(eqIdx));
+ faceWriter_->addPointData(deltaFace_[eqIdx], "delta_" + std::to_string(eqIdx));
+ faceWriter_->addPointData(defFace_[eqIdx], "defect_" + std::to_string(eqIdx));
+ }
+
+ for (std::size_t dofIdxGlobal = 0; dofIdxGlobal < deltaU[FVGridGeometry::faceIdx()].size(); ++dofIdxGlobal)
+ {
+ for (int eqIdx = 0; eqIdx < numEqFace; ++eqIdx)
+ {
+ xFace_[eqIdx][dofIdxGlobal] = uLastIter[FVGridGeometry::faceIdx()][dofIdxGlobal][eqIdx];
+ deltaFace_[eqIdx][dofIdxGlobal] = - deltaU[FVGridGeometry::faceIdx()][dofIdxGlobal][eqIdx];
+ defFace_[eqIdx][dofIdxGlobal] = residual[FVGridGeometry::faceIdx()][dofIdxGlobal][eqIdx];
+ }
+ }
+
+ ccWriter_.write(static_cast<double>(id_) + static_cast<double>(iteration_)/1000);
+ faceSequenceWriter_.write(static_cast<double>(id_) + static_cast<double>(iteration_)/1000);
+ ++iteration_;
+ }
+
+private:
+ std::size_t id_ = 0UL;
+ std::size_t iteration_ = 0UL;
+
+ const FVGridGeometry& fvGridGeometry_;
+
+ Dune::VTKSequenceWriter<GridView> ccWriter_;
+
+ std::vector<GlobalPosition> coordinates_;
+ std::shared_ptr<PointCloudVtkWriter<Scalar, GlobalPosition>> faceWriter_;
+ VTKSequenceWriter<PointCloudVtkWriter<Scalar, GlobalPosition>> faceSequenceWriter_;
+
+ std::array<std::vector<Scalar>, numEqCellCenter> defCellCenter_;
+ std::array<std::vector<Scalar>, numEqCellCenter> deltaCellCenter_;
+ std::array<std::vector<Scalar>, numEqCellCenter> xCellCenter_;
+
+ std::array<std::vector<Scalar>, numEqFace> defFace_;
+ std::array<std::vector<Scalar>, numEqFace> deltaFace_;
+ std::array<std::vector<Scalar>, numEqFace> xFace_;
+};
+
+} // end namespace Dumux
+
+#endif
diff --git a/test/freeflow/navierstokes/donea/main.cc b/test/freeflow/navierstokes/donea/main.cc
index 79d9334fb7791d01f6731bb18bc5226db1314d60..80e92746d2ca71a9ee2f3f4ff28f37c0a3d55564 100644
--- a/test/freeflow/navierstokes/donea/main.cc
+++ b/test/freeflow/navierstokes/donea/main.cc
@@ -43,6 +43,7 @@
#include <dumux/io/staggeredvtkoutputmodule.hh>
#include <dumux/linear/seqsolverbackend.hh>
#include <dumux/nonlinear/newtonsolver.hh>
+#include <dumux/nonlinear/staggerednewtonconvergencewriter.hh>
#include "problem.hh"
@@ -116,6 +117,10 @@ int main(int argc, char** argv) try
using NewtonSolver = Dumux::NewtonSolver<Assembler, LinearSolver>;
NewtonSolver nonLinearSolver(assembler, linearSolver);
+ using NewtonConvergenceWriter = StaggeredNewtonConvergenceWriter<FVGridGeometry, SolutionVector>;
+ auto convergenceWriter = std::make_shared<NewtonConvergenceWriter>(*fvGridGeometry);
+ nonLinearSolver.attachConvergenceWriter(convergenceWriter);
+
// linearize & solve
Dune::Timer timer;
nonLinearSolver.solve(x);
diff --git a/test/porousmediumflow/2p/implicit/nonisothermal/main.cc b/test/porousmediumflow/2p/implicit/nonisothermal/main.cc
index 17dffc5ede1e99e79324b148bada652ca8cf4dfc..3c5aa3face8e33f9969c09fe9548039e215aa34d 100644
--- a/test/porousmediumflow/2p/implicit/nonisothermal/main.cc
+++ b/test/porousmediumflow/2p/implicit/nonisothermal/main.cc
@@ -155,9 +155,9 @@ int main(int argc, char** argv) try
NewtonSolver nonLinearSolver(assembler, linearSolver);
//the convergence writer
- using GridView = GetPropType<TypeTag, Properties::GridView>;
- using NewtonConvergenceWriter = Dumux::NewtonConvergenceWriter<GridView, SolutionVector>;
- auto convergenceWriter = std::make_shared<NewtonConvergenceWriter>(leafGridView, fvGridGeometry->numDofs());
+ using NewtonConvergenceWriter = Dumux::NewtonConvergenceWriter<FVGridGeometry, SolutionVector>;
+ auto convergenceWriter = std::make_shared<NewtonConvergenceWriter>(*fvGridGeometry);
+ nonLinearSolver.attachConvergenceWriter(convergenceWriter);
// time loop
timeLoop->start(); do
@@ -165,8 +165,11 @@ int main(int argc, char** argv) try
// set previous solution for storage evaluations
assembler->setPreviousSolution(xOld);
+ // reset the convergence writer so that each time can be easily identified in the pvd file
+ convergenceWriter->reset(timeLoop->timeStepIndex());
+
// solve the non-linear system with time step control
- nonLinearSolver.solve(x, *timeLoop, convergenceWriter);
+ nonLinearSolver.solve(x, *timeLoop);
// make the new solution the old solution
xOld = x;