From ec5e8e82f07b2901aa966b0afe4f57bfdd36e646 Mon Sep 17 00:00:00 2001
From: Timo Koch <timo.koch@iws.uni-stuttgart.de>
Date: Tue, 14 Nov 2017 19:18:31 +0100
Subject: [PATCH] [1p][test] make 1d3dnetwork test pass again
---
.../1p/implicit/test_1pfv_network1d3d.cc | 3 ++
.../1p/implicit/tubesconvergencetest.py | 5 +-
.../1p/implicit/tubesproblem.hh | 54 ++++++++-----------
3 files changed, 27 insertions(+), 35 deletions(-)
diff --git a/test/porousmediumflow/1p/implicit/test_1pfv_network1d3d.cc b/test/porousmediumflow/1p/implicit/test_1pfv_network1d3d.cc
index 9c22775178..087fce42d1 100644
--- a/test/porousmediumflow/1p/implicit/test_1pfv_network1d3d.cc
+++ b/test/porousmediumflow/1p/implicit/test_1pfv_network1d3d.cc
@@ -186,6 +186,9 @@ int main(int argc, char** argv) try
<< "have been saved to restart files.");
}
+ // output l2 norm for convergence analysis
+ problem->outputL2Norm(x);
+
// make the new solution the old solution
xOld = x;
gridVariables->advanceTimeStep();
diff --git a/test/porousmediumflow/1p/implicit/tubesconvergencetest.py b/test/porousmediumflow/1p/implicit/tubesconvergencetest.py
index 066651907e..bcb21d0864 100755
--- a/test/porousmediumflow/1p/implicit/tubesconvergencetest.py
+++ b/test/porousmediumflow/1p/implicit/tubesconvergencetest.py
@@ -4,11 +4,12 @@ from math import *
import subprocess
import sys
-if len(sys.argv) != 2:
+if len(sys.argv) < 2:
sys.stderr.write('Please provide a single argument <testname> to the script\n')
sys.exit(1)
testname = str(sys.argv[1])
+testargs = [str(i) for i in sys.argv][2:]
# remove the old log file
subprocess.call(['rm', testname + '.log'])
@@ -16,7 +17,7 @@ print("Removed old log file ({})!".format(testname + '.log'))
# do the runs with different refinement
for i in [0, 1, 2, 3, 4, 5]:
- subprocess.call(['./' + testname, '-Grid.Refinement', str(i),
+ subprocess.call(['./' + testname] + testargs + ['-Grid.Refinement', str(i),
'-Problem.Name', testname])
# check the rates and append them to the log file
diff --git a/test/porousmediumflow/1p/implicit/tubesproblem.hh b/test/porousmediumflow/1p/implicit/tubesproblem.hh
index e261f79fb9..6b601eb815 100644
--- a/test/porousmediumflow/1p/implicit/tubesproblem.hh
+++ b/test/porousmediumflow/1p/implicit/tubesproblem.hh
@@ -29,6 +29,8 @@
#include <dune/geometry/quadraturerules.hh>
#include <dumux/discretization/cellcentered/tpfa/properties.hh>
+#include <dumux/discretization/box/properties.hh>
+#include <dumux/discretization/methods.hh>
#include <dumux/porousmediumflow/1p/implicit/model.hh>
#include <dumux/porousmediumflow/problem.hh>
#include <dumux/material/components/constant.hh>
@@ -94,23 +96,24 @@ class TubesTestProblem : public PorousMediumFlowProblem<TypeTag>
using PrimaryVariables = typename GET_PROP_TYPE(TypeTag, PrimaryVariables);
using BoundaryTypes = typename GET_PROP_TYPE(TypeTag, BoundaryTypes);
- using TimeManager = typename GET_PROP_TYPE(TypeTag, TimeManager);
using Element = typename GridView::template Codim<0>::Entity;
+ using FVGridGeometry = typename GET_PROP_TYPE(TypeTag, FVGridGeometry);
+ using SolutionVector = typename GET_PROP_TYPE(TypeTag, SolutionVector);
using FVElementGeometry = typename GET_PROP_TYPE(TypeTag, FVElementGeometry);
using SubControlVolume = typename GET_PROP_TYPE(TypeTag, SubControlVolume);
using GlobalPosition = Dune::FieldVector<Scalar, dimWorld>;
- enum { isBox = GET_PROP_VALUE(TypeTag, ImplicitIsBox) };
+ enum { isBox = GET_PROP_VALUE(TypeTag, DiscretizationMethod) == DiscretizationMethods::Box };
public:
TubesTestProblem(std::shared_ptr<const FVGridGeometry> fvGridGeometry)
: ParentType(fvGridGeometry)
{
- name_ = getParam<std::string>( "Problem.Name");
+ name_ = getParam<std::string>("Problem.Name");
//get hMax_ of the grid
hMax_ = 0.0;
- for (const auto& element : elements(this->gridView()))
+ for (const auto& element : elements(fvGridGeometry->gridView()))
hMax_ = std::max(element.geometry().volume(), hMax_);
}
@@ -185,19 +188,6 @@ public:
return PrimaryVariables(0.0);
}
- /*!
- * \brief Evaluate the boundary conditions for a neumann
- * boundary segment.
- *
- * For this method, the \a priVars parameter stores the mass flux
- * in normal direction of each component. Negative values mean
- * influx.
- */
- PrimaryVariables neumannAtPos(const GlobalPosition& globalPos) const
- {
- return PrimaryVariables(0.0);
- }
-
// \}
/*!
@@ -232,8 +222,7 @@ public:
const auto& globalPos = scv.center();
const auto& volVars = elemVolVars[scv];
- const auto K = this->spatialParams().permeability(element, scv,
- this->model().elementSolution(element, this->model().curSol()));
+ const auto K = this->spatialParams().permeability(element, scv, ElementSolutionVector());
using std::sin;
if (globalPos[2] > 0.5 - eps_)
@@ -260,19 +249,18 @@ public:
// \}
- void postTimeStep()
+ void outputL2Norm(const SolutionVector solution) const
{
- const auto& solution = this->model().curSol();
-
// calculate the discrete L2-Norm
- Scalar LTwoNorm = 0.0;
+ Scalar lTwoNorm = 0.0;
// get the Gaussian quadrature rule for intervals
const auto& quad = Dune::QuadratureRules<Scalar, dim>::rule(Dune::GeometryType(1), 1);
- for (const auto& element : elements(this->gridView()))
+ const auto& gg = this->fvGridGeometry();
+ for (const auto& element : elements(gg.gridView()))
{
- const unsigned int eIdx = this->elementMapper().index(element);
+ const auto eIdx = gg.elementMapper().index(element);
const auto geometry = element.geometry();
for(auto&& qp : quad)
{
@@ -289,22 +277,23 @@ public:
const auto& localFiniteElement = feCache_.get(geometry.type());
localFiniteElement.localBasis().evaluateFunction(qp.position(), shapeValues);
for (unsigned int i = 0; i < shapeValues.size(); ++i)
- p += shapeValues[i]*solution[this->model().dofMapper().subIndex(element, i, dim)][pressureIdx];
+ p += shapeValues[i]*solution[gg.dofMapper().subIndex(element, i, dim)][pressureIdx];
}
- LTwoNorm += (p - pe)*(p - pe)*qp.weight()*integrationElement;
+ lTwoNorm += (p - pe)*(p - pe)*qp.weight()*integrationElement;
}
}
- LTwoNorm = std::sqrt(LTwoNorm);
+ lTwoNorm = std::sqrt(lTwoNorm);
// write the norm into a log file
- logFile_.open(this->name() + ".log", std::ios::app);
- logFile_ << "[ConvergenceTest] L2-norm(pressure) = " << LTwoNorm << " hMax = " << hMax_ << std::endl;
- logFile_.close();
+ std::ofstream logFile;
+ logFile.open(this->name() + ".log", std::ios::app);
+ logFile << "[ConvergenceTest] L2-norm(pressure) = " << lTwoNorm << " hMax = " << hMax_ << std::endl;
+ logFile.close();
}
private:
- Scalar exactPressure_(const GlobalPosition& globalPos)
+ Scalar exactPressure_(const GlobalPosition& globalPos) const
{
using std::sin;
return sin(4.0*M_PI*globalPos[2]);
@@ -315,7 +304,6 @@ private:
Scalar hMax_;
- std::ofstream logFile_;
typename Dune::PQkLocalFiniteElementCache<Scalar, Scalar, dim, 1> feCache_;
};
--
GitLab