From cd375a872ef59a17cfd8c95bdc360031c680bab1 Mon Sep 17 00:00:00 2001
From: Mathis Kelm <mathis.kelm@iws.uni-stuttgart.de>
Date: Wed, 13 Dec 2023 14:26:52 +0100
Subject: [PATCH] [boxdfm] remove corner storage in scv/scvf
---
dumux/experimental/ode/odesolver.hh | 226 ++++++++++++++++++
.../boxdfm/fvelementgeometry.hh | 64 +++--
.../boxdfm/subcontrolvolume.hh | 46 +---
.../boxdfm/subcontrolvolumeface.hh | 60 ++---
4 files changed, 307 insertions(+), 89 deletions(-)
create mode 100644 dumux/experimental/ode/odesolver.hh
diff --git a/dumux/experimental/ode/odesolver.hh b/dumux/experimental/ode/odesolver.hh
new file mode 100644
index 0000000000..bd1fc227c1
--- /dev/null
+++ b/dumux/experimental/ode/odesolver.hh
@@ -0,0 +1,226 @@
+//
+// SPDX-FileCopyrightInfo: Copyright © DuMux Project contributors, see AUTHORS.md in root folder
+// SPDX-License-Identifier: GPL-3.0-or-later
+//
+#include <cmath>
+#include <cassert>
+#include <iostream>
+#include <algorithm>
+#include <utility>
+#include <memory>
+#include <array>
+
+#include <dune/common/float_cmp.hh>
+#include <dune/common/exceptions.hh>
+
+#include <dumux/io/format.hh>
+#include <dumux/common/initialize.hh>
+#include <dumux/experimental/common/variables.hh>
+#include <dumux/nonlinear/newtonsolver.hh>
+
+#include <dumux/experimental/timestepping/timelevel.hh>
+#include <dumux/experimental/timestepping/multistagemethods.hh>
+#include <dumux/experimental/timestepping/multistagetimestepper.hh>
+
+namespace Dumux::Experimental {
+
+class ScalarAssembler
+{
+
+public:
+ using Scalar = double;
+ using SolutionVector = Scalar;
+ using ResidualType = Scalar;
+ using JacobianMatrix = Scalar;
+ using Variables = Experimental::Variables<SolutionVector>;
+ using StageParams = Experimental::MultiStageParams<Scalar>;
+
+ void setLinearSystem() {}
+ JacobianMatrix& jacobian() { return jac_; }
+ ResidualType& residual() { return res_; }
+
+ void assembleResidual(const Variables& vars)
+ {
+ if (stageParams_->size() != spatialResiduals_.size())
+ DUNE_THROW(Dune::InvalidStateException, "Wrong number of residuals");
+
+ // assemble residual components depending on current solution
+ assembleResiduals_(vars,
+ temporalResiduals_.back(),
+ spatialResiduals_.back()
+ );
+
+ // assemble residual for current stage solver
+ res_ = 0.0;
+ for (std::size_t k = 0; k < stageParams_->size(); ++k)
+ {
+ if (!stageParams_->skipTemporal(k))
+ res_ += stageParams_->temporalWeight(k)*temporalResiduals_[k];
+ if (!stageParams_->skipSpatial(k))
+ res_ += stageParams_->spatialWeight(k)*spatialResiduals_[k];
+ }
+ }
+
+ void assembleJacobianAndResidual(const Variables& vars)
+ {
+ assembleResidual(vars);
+ jac_ = 1.0;
+ }
+
+ void prepareStage(Variables& variables,
+ std::shared_ptr<const StageParams> params)
+ {
+ const auto curStage = params->size() - 1;
+ std::cout << "Assembler is preparing stage " << curStage << std::endl;
+ const auto prevStage = stageParams_ ? stageParams_->size() - 1 : 0;
+ if (curStage != prevStage+1)
+ DUNE_THROW(Dune::InvalidStateException,
+ "Can only prepare stages consecutively (current stage: " << curStage
+ << ", previous stage: " << prevStage << ")");
+
+ stageParams_ = params;
+
+ // in the first stage, also assemble the old residual
+ if (curStage == 1)
+ {
+ // we update the time level of the given grid variables
+ const auto t = params->timeAtStage(0);
+ const auto prevT = params->timeAtStage(0);
+ const auto dtFraction = params->timeStepFraction(0);
+ variables.updateTime(Experimental::TimeLevel{t, prevT, dtFraction});
+
+ // allocate memory for residuals of stage 0
+ spatialResiduals_.emplace_back(0.0);
+ temporalResiduals_.emplace_back(0.0);
+
+ // assemble stage 0 residuals
+ assembleResiduals_(variables,
+ temporalResiduals_.back(),
+ spatialResiduals_.back()
+ );
+ }
+
+ // we update the time level of the given grid variables
+ const auto t = params->timeAtStage(curStage);
+ const auto prevT = params->timeAtStage(0);
+ const auto dtFraction = params->timeStepFraction(curStage);
+ variables.updateTime(Experimental::TimeLevel{t, prevT, dtFraction});
+
+ // allocate memory for residuals of the upcoming stage
+ spatialResiduals_.emplace_back(0.0);
+ temporalResiduals_.emplace_back(0.0);
+ }
+
+ void clearStages()
+ {
+ std::cout << "Assembler is clearing stage data " << std::endl;
+ spatialResiduals_.clear();
+ temporalResiduals_.clear();
+ stageParams_.reset();
+ }
+
+private:
+ void assembleResiduals_(const Variables& variables,
+ ResidualType& temporal,
+ ResidualType& spatial) const
+ {
+ using std::exp;
+ temporal = variables.dofs();
+ spatial = -exp(variables.timeLevel().current());
+ }
+
+ ResidualType res_;
+ JacobianMatrix jac_;
+ std::vector<ResidualType> spatialResiduals_;
+ std::vector<ResidualType> temporalResiduals_;
+ std::shared_ptr<const StageParams> stageParams_;
+};
+
+class ScalarLinearSolver
+{
+public:
+ void setResidualReduction(double residualReduction) {}
+
+ bool solve(const double& A, double& x, const double& b) const
+ { x = b/A; return true; }
+
+ double norm(const double residual) const
+ {
+ using std::abs;
+ return abs(residual);
+ }
+};
+
+} // end namespace Dumux
+
+int main(int argc, char* argv[])
+{
+ using namespace Dumux;
+
+ // maybe initialize MPI and/or multithreading backend
+ Dumux::initialize(argc, argv);
+
+ using Assembler = ScalarAssembler;
+ using LinearSolver = ScalarLinearSolver;
+ using NewtonSolver = NewtonSolver<Assembler, LinearSolver, DefaultPartialReassembler>;
+
+ auto assembler = std::make_shared<Assembler>();
+ auto linearSolver = std::make_shared<LinearSolver>();
+ auto newtonSolver = std::make_shared<NewtonSolver>(assembler, linearSolver);
+
+ // initial solution and variables
+ using Scalar = typename Assembler::Scalar;
+ using Variables = typename Assembler::Variables;
+ using SolutionVector = typename Variables::SolutionVector;
+
+ const auto exact = [] (const Scalar t)
+ {
+ using std::exp;
+ return exp(t) - 1.0;
+ };
+
+ const auto computeError = [&] (const Variables& vars)
+ {
+ using std::abs;
+ const auto time = vars.timeLevel().current();
+ const auto exactSol = exact(time);
+ const auto absErr = abs(vars.dofs()-exactSol);
+ const auto relErr = absErr/exactSol;
+ return std::make_pair(absErr, relErr);
+ };
+
+ const auto testIntegration = [&] (auto method, Scalar expectedRelError)
+ {
+ std::cout << "\n-- Integration with " << method->name() << ":\n\n";
+ SolutionVector x = 0.0;
+ Variables vars(x);
+
+ using TimeStepper = Experimental::MultiStageTimeStepper<NewtonSolver>;
+ TimeStepper timeStepper(newtonSolver, method);
+
+ const Scalar dt = 0.01;
+ timeStepper.step(vars, /*time*/0.0, dt);
+
+ const auto [abs, rel] = computeError(vars);
+ std::cout << "\n"
+ << "-- Summary\n"
+ << "-- ===========================\n"
+ << "-- Exact solution: " << Fmt::format("{:.4e}\n", exact(dt))
+ << "-- Discrete solution: " << Fmt::format("{:.4e}\n", vars.dofs())
+ << "-- Absolute error: " << Fmt::format("{:.4e}\n", abs)
+ << "-- Relative error: " << Fmt::format("{:.4e}", rel) << std::endl;
+
+ if (Dune::FloatCmp::ne(rel, expectedRelError, 0.01))
+ DUNE_THROW(Dune::InvalidStateException,
+ "Detected deviation from reference > 1% for " << method->name());
+ };
+
+ using namespace Experimental::MultiStage;
+ testIntegration(std::make_shared<ExplicitEuler<Scalar>>(), 4.9917e-03);
+ testIntegration(std::make_shared<ImplicitEuler<Scalar>>(), 5.0083e-03);
+ testIntegration(std::make_shared<Theta<Scalar>>(0.5), 8.3333e-06);
+ testIntegration(std::make_shared<DIRKThirdOrderAlexander<Scalar>>(), 7.9214e-09);
+ testIntegration(std::make_shared<RungeKuttaExplicitFourthOrder<Scalar>>(), 3.4829e-12);
+
+ return 0;
+}
diff --git a/dumux/porousmediumflow/boxdfm/fvelementgeometry.hh b/dumux/porousmediumflow/boxdfm/fvelementgeometry.hh
index e564d16526..39bac01b7b 100644
--- a/dumux/porousmediumflow/boxdfm/fvelementgeometry.hh
+++ b/dumux/porousmediumflow/boxdfm/fvelementgeometry.hh
@@ -171,20 +171,34 @@ public:
const Element& element() const
{ return *element_; }
- // suppress warnings due to current implementation
- // these interfaces should be used!
- #pragma GCC diagnostic push
- #pragma GCC diagnostic ignored "-Wdeprecated-declarations"
-
//! Create the geometry of a given sub control volume
typename SubControlVolume::Traits::Geometry geometry(const SubControlVolume& scv) const
- { return scv.geometry(); }
+ {
+ using ScvGeometry = typename SubControlVolume::Traits::Geometry;
+ if (scv.isOnFracture())
+ DUNE_THROW(Dune::InvalidStateException, "The geometry object cannot be defined for fracture scvs "
+ "because the number of known corners is insufficient. "
+ "You can do this manually by extract the corners from this scv "
+ "and extrude them by the corresponding aperture. ");
+
+ const typename GG::GeometryHelper geometryHelper(element().geometry());
+ const auto corners = geometryHelper.getScvCorners(scv.index());
+ return ScvGeometry(Dune::GeometryTypes::cube(ScvGeometry::mydimension), corners);
+ }
//! Create the geometry of a given sub control volume face
typename SubControlVolumeFace::Traits::Geometry geometry(const SubControlVolumeFace& scvf) const
- { return scvf.geometry(); }
-
- #pragma GCC diagnostic pop
+ {
+ using ScvfGeometry = typename SubControlVolumeFace::Traits::Geometry;
+ if (scvf.isOnFracture())
+ DUNE_THROW(Dune::InvalidStateException, "The geometry object cannot be defined for fracture scvs "
+ "because the number of known corners is insufficient. "
+ "You can do this manually by extracting the corners from this scv "
+ "and extruding them by the corresponding aperture. ");
+ const typename GG::GeometryHelper geometryHelper(element().geometry());
+ const auto corners = geometryHelper.getScvfCorners(scvf.indexInElement());
+ return ScvfGeometry(Dune::GeometryTypes::cube(ScvfGeometry::mydimension), corners);
+ }
private:
const GridGeometry* gridGeometryPtr_;
@@ -327,20 +341,34 @@ public:
const Element& element() const
{ return *element_; }
- // suppress warnings due to current implementation
- // these interfaces should be used!
- #pragma GCC diagnostic push
- #pragma GCC diagnostic ignored "-Wdeprecated-declarations"
-
//! Create the geometry of a given sub control volume
typename SubControlVolume::Traits::Geometry geometry(const SubControlVolume& scv) const
- { return scv.geometry(); }
+ {
+ using ScvGeometry = typename SubControlVolume::Traits::Geometry;
+ if (scv.isOnFracture())
+ DUNE_THROW(Dune::InvalidStateException, "The geometry object cannot be defined for fracture scvs "
+ "because the number of known corners is insufficient. "
+ "You can do this manually by extract the corners from this scv "
+ "and extrude them by the corresponding aperture. ");
+
+ const GeometryHelper geometryHelper(element().geometry());
+ const auto corners = geometryHelper.getScvCorners(scv.index());
+ return ScvGeometry(Dune::GeometryTypes::cube(ScvGeometry::mydimension), corners);
+ }
//! Create the geometry of a given sub control volume face
typename SubControlVolumeFace::Traits::Geometry geometry(const SubControlVolumeFace& scvf) const
- { return scvf.geometry(); }
-
- #pragma GCC diagnostic pop
+ {
+ using ScvfGeometry = typename SubControlVolumeFace::Traits::Geometry;
+ if (scvf.isOnFracture())
+ DUNE_THROW(Dune::InvalidStateException, "The geometry object cannot be defined for fracture scvs "
+ "because the number of known corners is insufficient. "
+ "You can do this manually by extracting the corners from this scv "
+ "and extruding them by the corresponding aperture. ");
+ const GeometryHelper geometryHelper(element().geometry());
+ const auto corners = geometryHelper.getScvfCorners(scvf.indexInElement());
+ return ScvfGeometry(Dune::GeometryTypes::cube(ScvfGeometry::mydimension), corners);
+ }
private:
diff --git a/dumux/porousmediumflow/boxdfm/subcontrolvolume.hh b/dumux/porousmediumflow/boxdfm/subcontrolvolume.hh
index cf8504ea5a..38a7c2d5d0 100644
--- a/dumux/porousmediumflow/boxdfm/subcontrolvolume.hh
+++ b/dumux/porousmediumflow/boxdfm/subcontrolvolume.hh
@@ -69,7 +69,6 @@ class BoxDfmSubControlVolume
using LocalIndexType = typename T::LocalIndexType;
using Scalar = typename T::Scalar;
using GlobalPosition = typename T::GlobalPosition;
- using CornerStorage = typename T::CornerStorage;
enum { dim = Geometry::mydimension };
static_assert(dim == 2 || dim == 3, "Box-Dfm sub-control volume only implemented in 2d or 3d");
@@ -88,22 +87,22 @@ public:
GridIndexType elementIndex,
GridIndexType dofIndex)
: isFractureScv_(false)
- , corners_(geometryHelper.getScvCorners(scvIdx))
, center_(0.0)
- , volume_(Dumux::convexPolytopeVolume<T::dim>(
- Dune::GeometryTypes::cube(T::dim),
- [&](unsigned int i){ return corners_[i]; })
- )
, elementIndex_(elementIndex)
, vIdxLocal_(scvIdx)
, elemLocalScvIdx_(scvIdx)
, dofIndex_(dofIndex)
, facetIdx_(0)
{
+ auto corners = geometryHelper.getScvCorners(scvIdx);
+ dofPosition_ = corners[0];
+ volume_ = Dumux::convexPolytopeVolume<T::dim>(
+ Dune::GeometryTypes::cube(T::dim),
+ [&](unsigned int i){ return corners[i]; });
// compute center point
- for (const auto& corner : corners_)
+ for (const auto& corner : corners)
center_ += corner;
- center_ /= corners_.size();
+ center_ /= corners.size();
}
/*!
@@ -126,7 +125,6 @@ public:
GridIndexType elementIndex,
GridIndexType dofIndex)
: isFractureScv_(true)
- , corners_()
, center_(0.0)
, volume_(0.0)
, elementIndex_(elementIndex)
@@ -135,21 +133,17 @@ public:
, dofIndex_(dofIndex)
, facetIdx_(elemLocalFacetIdx)
{
- // copy corners
auto corners = geometryHelper.getBoundaryScvfCorners(intersection.indexInInside(), indexInIntersection);
- const auto numCorners = corners.size();
- corners_.resize(numCorners);
- for (unsigned int i = 0; i < numCorners; ++i)
- corners_[i] = corners[i];
+ dofPosition_ = corners[0];
// compute volume and scv center
volume_ = Dumux::convexPolytopeVolume<T::dim-1>(
Dune::GeometryTypes::cube(T::dim-1),
- [&](unsigned int i){ return corners_[i]; }
+ [&](unsigned int i){ return corners[i]; }
);
- for (const auto& corner : corners_)
+ for (const auto& corner : corners)
center_ += corner;
- center_ /= corners_.size();
+ center_ /= corners.size();
}
//! The center of the sub control volume
@@ -178,7 +172,7 @@ public:
// The position of the dof this scv is embedded in (list is defined such that first entry is vertex itself)
const GlobalPosition& dofPosition() const
- { return corners_[0]; }
+ { return dofPosition_; }
//! The global index of the element this scv is embedded in
GridIndexType elementIndex() const
@@ -188,23 +182,9 @@ public:
bool isOnFracture() const
{ return isFractureScv_; }
- //! The geometry of the sub control volume
- // e.g. for integration
- [[deprecated("Will be removed after 3.7. Use fvGeometry.geometry(scv).")]]
- Geometry geometry() const
- {
- if (isFractureScv_)
- DUNE_THROW(Dune::InvalidStateException, "The geometry object cannot be defined for fracture scvs "
- "because the number of known corners is insufficient. "
- "You can do this manually by extract the corners from this scv "
- "and extrude them by the corresponding aperture. ");
-
- return Geometry(Dune::GeometryTypes::cube(dim), corners_);
- }
-
private:
bool isFractureScv_;
- CornerStorage corners_;
+ GlobalPosition dofPosition_;
GlobalPosition center_;
Scalar volume_;
GridIndexType elementIndex_;
diff --git a/dumux/porousmediumflow/boxdfm/subcontrolvolumeface.hh b/dumux/porousmediumflow/boxdfm/subcontrolvolumeface.hh
index 85fb5af05d..4fecd088bc 100644
--- a/dumux/porousmediumflow/boxdfm/subcontrolvolumeface.hh
+++ b/dumux/porousmediumflow/boxdfm/subcontrolvolumeface.hh
@@ -69,7 +69,6 @@ class BoxDfmSubControlVolumeFace
using LocalIndexType = typename T::LocalIndexType;
using Scalar = typename T::Scalar;
using GlobalPosition = typename T::GlobalPosition;
- using CornerStorage = typename T::CornerStorage;
using Geometry = typename T::Geometry;
using BoundaryFlag = typename T::BoundaryFlag;
@@ -89,13 +88,7 @@ public:
const typename Element::Geometry& elemGeometry,
GridIndexType scvfIndex,
std::vector<LocalIndexType>&& scvIndices)
- : corners_(geometryHelper.getScvfCorners(scvfIndex))
- , center_(0.0)
- , unitOuterNormal_(geometryHelper.normal(corners_, scvIndices))
- , area_(Dumux::convexPolytopeVolume<T::dim-1>(
- Dune::GeometryTypes::cube(T::dim-1),
- [&](unsigned int i){ return corners_[i]; })
- )
+ : center_(0.0)
, scvfIndex_(scvfIndex)
, scvIndices_(std::move(scvIndices))
, boundary_(false)
@@ -103,9 +96,15 @@ public:
, boundaryFlag_{}
, facetIdx_(0)
{
- for (const auto& corner : corners_)
+ auto corners = geometryHelper.getScvfCorners(scvfIndex);
+ unitOuterNormal_ = geometryHelper.normal(corners, scvIndices_);
+ area_ = Dumux::convexPolytopeVolume<T::dim-1>(
+ Dune::GeometryTypes::cube(T::dim-1),
+ [&](unsigned int i){ return corners[i]; });
+
+ for (const auto& corner : corners)
center_ += corner;
- center_ /= corners_.size();
+ center_ /= corners.size();
}
//! Constructor for boundary scvfs
@@ -116,13 +115,8 @@ public:
LocalIndexType indexInIntersection,
GridIndexType scvfIndex,
std::vector<LocalIndexType>&& scvIndices)
- : corners_(geometryHelper.getBoundaryScvfCorners(intersection.indexInInside(), indexInIntersection))
- , center_(0.0)
+ : center_(0.0)
, unitOuterNormal_(intersection.centerUnitOuterNormal())
- , area_(Dumux::convexPolytopeVolume<T::dim-1>(
- Dune::GeometryTypes::cube(T::dim-1),
- [&](unsigned int i){ return corners_[i]; })
- )
, scvfIndex_(scvfIndex)
, scvIndices_(std::move(scvIndices))
, boundary_(true)
@@ -130,9 +124,13 @@ public:
, boundaryFlag_{intersection}
, facetIdx_(0)
{
- for (const auto& corner : corners_)
+ auto corners = geometryHelper.getBoundaryScvfCorners(intersection.indexInInside(), indexInIntersection);
+ area_ = Dumux::convexPolytopeVolume<T::dim-1>(
+ Dune::GeometryTypes::cube(T::dim-1),
+ [&](unsigned int i){ return corners[i]; });
+ for (const auto& corner : corners)
center_ += corner;
- center_ /= corners_.size();
+ center_ /= corners.size();
}
//! Constructor for inner fracture scvfs
@@ -144,8 +142,7 @@ public:
GridIndexType scvfIndex,
std::vector<LocalIndexType>&& scvIndices,
bool boundary)
- : corners_(geometryHelper.getFractureScvfCorners(intersection.indexInInside(), indexInIntersection))
- , center_(0.0)
+ : center_(0.0)
, scvfIndex_(scvfIndex)
, scvIndices_(std::move(scvIndices))
, boundary_(boundary)
@@ -153,21 +150,22 @@ public:
, boundaryFlag_{intersection}
, facetIdx_(intersection.indexInInside())
{
+ auto corners = geometryHelper.getFractureScvfCorners(intersection.indexInInside(), indexInIntersection);
// The area here is given in meters. In order to
// get the right dimensions, the user has to provide
// the appropriate aperture in the problem (via an extrusion factor)
if (T::dim == 3)
- area_ = (corners_[1]-corners_[0]).two_norm();
+ area_ = (corners[1]-corners[0]).two_norm();
else if (T::dim == 2)
area_ = 1.0;
// obtain the unit normal vector
- unitOuterNormal_ = geometryHelper.fractureNormal(corners_, intersection, indexInIntersection);
+ unitOuterNormal_ = geometryHelper.fractureNormal(corners, intersection, indexInIntersection);
// compute the scvf center
- for (const auto& corner : corners_)
+ for (const auto& corner : corners)
center_ += corner;
- center_ /= corners_.size();
+ center_ /= corners.size();
}
//! The center of the sub control volume face
@@ -224,21 +222,7 @@ public:
return static_cast<std::size_t>(!boundary());
}
- //! The geometry of the sub control volume face
- [[deprecated("Will be removed after 3.7. Use fvGeometry.geometry(scvf).")]]
- Geometry geometry() const
- {
- if (isFractureScvf_)
- DUNE_THROW(Dune::InvalidStateException, "The geometry object cannot be defined for fracture scvs "
- "because the number of known corners is insufficient. "
- "You can do this manually by extract the corners from this scv "
- "and extrude them by the corresponding aperture. ");
-
- return Geometry(Dune::GeometryTypes::cube(Geometry::mydimension), corners_);
- }
-
private:
- CornerStorage corners_;
GlobalPosition center_;
GlobalPosition unitOuterNormal_;
Scalar area_;
--
GitLab