From 03ab43dbd607b8270ec37baeee4a129fed4b4d1b Mon Sep 17 00:00:00 2001
From: Katharina Heck <katharina.heck@iws.uni-stuttgart.de>
Date: Fri, 3 Nov 2017 16:50:11 +0100
Subject: [PATCH] [3p][test] add exact temperature calculation for ni tests
---
.../3p/implicit/3pniconductionproblem.hh | 46 +++++++++++--------
.../3p/implicit/3pniconvectionproblem.hh | 41 ++++++++++++-----
.../3p/implicit/test_box3pniconduction.cc | 19 ++++----
.../3p/implicit/test_box3pniconvection.cc | 5 ++
.../3p/implicit/test_cc3pniconduction.cc | 13 +++++-
.../3p/implicit/test_cc3pniconvection.cc | 13 +++++-
.../3p/implicit/test_cc3pniconvection.input | 3 ++
7 files changed, 97 insertions(+), 43 deletions(-)
diff --git a/test/porousmediumflow/3p/implicit/3pniconductionproblem.hh b/test/porousmediumflow/3p/implicit/3pniconductionproblem.hh
index f1a20768da..b237e14d31 100644
--- a/test/porousmediumflow/3p/implicit/3pniconductionproblem.hh
+++ b/test/porousmediumflow/3p/implicit/3pniconductionproblem.hh
@@ -112,6 +112,8 @@ class ThreePNIConductionProblem : public PorousMediumFlowProblem<TypeTag>
using ThermalConductivityModel = typename GET_PROP_TYPE(TypeTag, ThermalConductivityModel);
using VolumeVariables = typename GET_PROP_TYPE(TypeTag, VolumeVariables);
using ElementVolumeVariables = typename GET_PROP_TYPE(TypeTag, ElementVolumeVariables);
+ using ElementSolutionVector = typename GET_PROP_TYPE(TypeTag, ElementSolutionVector);
+ using SolutionVector = typename GET_PROP_TYPE(TypeTag, SolutionVector);
using SubControlVolumeFace = typename GET_PROP_TYPE(TypeTag, SubControlVolumeFace);
using IapwsH2O = H2O<Scalar>;
using NeumannFluxes = typename GET_PROP_TYPE(TypeTag, NumEqVector);
@@ -151,7 +153,8 @@ public:
outputInterval_ = getParam<int>("Problem.OutputInterval");
temperatureHigh_ = 300.0;
time_ = 0.0;
- }
+ temperatureExact_.resize(this->fvGridGeometry().gridView().size(GridView::dimension));
+ }
// get time from the mainfile timeloop
void setTime(Scalar time)
@@ -160,15 +163,24 @@ public:
/*!
* \brief Adds additional VTK output data to the VTKWriter. Function is called by the output module on every write.
*/
- const std::vector<Scalar>& getExactTemperature() const
+ /*!
+ * \brief Adds additional VTK output data to the VTKWriter. Function is called by the output module on every write.
+ */
+ const std::vector<Scalar>& getExactTemperature()
+ {
+ return temperatureExact_;
+ }
+
+ void updateExactTemperature(const SolutionVector& curSol)
{
const auto someElement = *(elements(this->fvGridGeometry().gridView()).begin());
- const auto someElemSol = this->model().elementSolution(someElement, this->model().curSol());
+
+ ElementSolutionVector someElemSol(someElement, curSol, this->fvGridGeometry());
const auto someInitSol = initialAtPos(someElement.geometry().center());
- auto someFvGeometry = localView(this->model().fvGridGeometry());
- someFvGeometry.bindElement(someElement);
- const auto someScv = *(scvs(someFvGeometry).begin());
+ auto fvGeometry = localView(this->fvGridGeometry());
+ fvGeometry.bindElement(someElement);
+ const auto someScv = *(scvs(fvGeometry).begin());
VolumeVariables volVars;
volVars.update(someElemSol, *this, someElement, someScv);
@@ -180,25 +192,23 @@ public:
const auto heatCapacityS = this->spatialParams().solidHeatCapacity(someElement, someScv, someElemSol);
const auto storage = densityW*heatCapacityW*porosity + densityS*heatCapacityS*(1 - porosity);
const auto effectiveThermalConductivity = ThermalConductivityModel::effectiveThermalConductivity(volVars, this->spatialParams(),
- someElement, someFvGeometry, someScv);
- using std::max;
-
- for (const auto& element : elements(this->gridView()))
+ someElement, fvGeometry, someScv);
+ for (const auto& element : elements(this->fvGridGeometry().gridView()))
{
- auto fvGeometry = localView(this->model().fvGridGeometry());
+ auto fvGeometry = localView(this->fvGridGeometry());
fvGeometry.bindElement(element);
for (auto&& scv : scvs(fvGeometry))
{
- auto globalIdx = scv.dofIndex();
- const auto& globalPos = scv.dofPosition();
- using std::erf;
- using std::sqrt;
- temperatureExact_[globalIdx] = temperatureHigh_ + (someInitSol[temperatureIdx] - temperatureHigh_)
+ auto globalIdx = scv.dofIndex();
+ const auto& globalPos = scv.dofPosition();
+ using std::erf;
+ using std::sqrt;
+ temperatureExact_[globalIdx] = temperatureHigh_ + (someInitSol[temperatureIdx] - temperatureHigh_)
*erf(0.5*sqrt(globalPos[0]*globalPos[0]*storage/time_/effectiveThermalConductivity));
+
}
}
- return temperatureExact_;
}
/*!
* \name Problem parameters
@@ -319,7 +329,7 @@ private:
static constexpr Scalar eps_ = 1e-6;
std::string name_;
int outputInterval_;
- std::vector<Scalar> temperatureExact_;
+ std::vector<double> temperatureExact_;
Scalar time_;
};
diff --git a/test/porousmediumflow/3p/implicit/3pniconvectionproblem.hh b/test/porousmediumflow/3p/implicit/3pniconvectionproblem.hh
index ce2e4bdc70..2bc1f5f254 100644
--- a/test/porousmediumflow/3p/implicit/3pniconvectionproblem.hh
+++ b/test/porousmediumflow/3p/implicit/3pniconvectionproblem.hh
@@ -109,6 +109,8 @@ class ThreePNIConvectionProblem : public PorousMediumFlowProblem<TypeTag>
using ThermalConductivityModel = typename GET_PROP_TYPE(TypeTag, ThermalConductivityModel);
using VolumeVariables = typename GET_PROP_TYPE(TypeTag, VolumeVariables);
using ElementVolumeVariables = typename GET_PROP_TYPE(TypeTag, ElementVolumeVariables);
+ using ElementSolutionVector = typename GET_PROP_TYPE(TypeTag, ElementSolutionVector);
+ using SolutionVector = typename GET_PROP_TYPE(TypeTag, SolutionVector);
using SubControlVolumeFace = typename GET_PROP_TYPE(TypeTag, SubControlVolumeFace);
using IapwsH2O = H2O<Scalar>;
@@ -155,24 +157,37 @@ public:
pressureHigh_ = 2e5;
pressureLow_ = 1e5;
+ time_ = 0.0;
+ temperatureExact_.resize(this->fvGridGeometry().gridView().size(GridView::dimension));
+
}
+ // get time from the mainfile timeloop
+ void setTime(Scalar time)
+ { time_ = time; }
+
+ /*!
+ * \brief Adds additional VTK output data to the VTKWriter. Function is called by the output module on every write.
+ */
/*!
* \brief Adds additional VTK output data to the VTKWriter. Function is called by the output module on every write.
*/
- template <class VtkOutputModule>
- void addVtkOutputFields(VtkOutputModule& outputModule) const
+ const std::vector<Scalar>& getExactTemperature()
+ {
+ return temperatureExact_;
+ }
+
+ void updateExactTemperature(const SolutionVector& curSol)
{
- auto& temperatureExact = outputModule.createScalarField("temperatureExact", dofCodim);
+ const auto someElement = *(elements(this->fvGridGeometry().gridView()).begin());
- const auto someElement = *(elements(this->gridView()).begin());
- const auto someElemSol = this->model().elementSolution(someElement, this->model().curSol());
+ ElementSolutionVector someElemSol(someElement, curSol, this->fvGridGeometry());
const auto someInitSol = initialAtPos(someElement.geometry().center());
- auto someFvGeometry = localView(this->model().fvGridGeometry());
- someFvGeometry.bindElement(someElement);
- const auto someScv = *(scvs(someFvGeometry).begin());
+ auto fvGeometry = localView(this->fvGridGeometry());
+ fvGeometry.bindElement(someElement);
+ const auto someScv = *(scvs(fvGeometry).begin());
VolumeVariables volVars;
volVars.update(someElemSol, *this, someElement, someScv);
@@ -187,19 +202,19 @@ public:
std::cout << "storage: " << storageTotal << '\n';
using std::max;
- const Scalar time = max(this->timeManager().time() + this->timeManager().timeStepSize(), 1e-10);
+ const Scalar time = max(time_, 1e-10);
const Scalar retardedFrontVelocity = darcyVelocity_*storageW/storageTotal/porosity;
std::cout << "retarded velocity: " << retardedFrontVelocity << '\n';
- for (const auto& element : elements(this->gridView()))
+ for (const auto& element : elements(this->fvGridGeometry().gridView()))
{
- auto fvGeometry = localView(this->model().fvGridGeometry());
+ auto fvGeometry = localView(this->fvGridGeometry());
fvGeometry.bindElement(element);
for (auto&& scv : scvs(fvGeometry))
{
auto dofIdxGlobal = scv.dofIndex();
auto dofPosition = scv.dofPosition();
- temperatureExact[dofIdxGlobal] = (dofPosition[0] < retardedFrontVelocity*time) ? temperatureHigh_ : temperatureLow_;
+ temperatureExact_[dofIdxGlobal] = (dofPosition[0] < retardedFrontVelocity*time) ? temperatureHigh_ : temperatureLow_;
}
}
}
@@ -338,6 +353,8 @@ private:
static constexpr Scalar eps_ = 1e-6;
std::string name_;
int outputInterval_;
+ std::vector<double> temperatureExact_;
+ Scalar time_;
};
} //end namespace
diff --git a/test/porousmediumflow/3p/implicit/test_box3pniconduction.cc b/test/porousmediumflow/3p/implicit/test_box3pniconduction.cc
index 2f964ae900..1271393c21 100644
--- a/test/porousmediumflow/3p/implicit/test_box3pniconduction.cc
+++ b/test/porousmediumflow/3p/implicit/test_box3pniconduction.cc
@@ -146,8 +146,11 @@ int main(int argc, char** argv) try
using VtkOutputFields = typename GET_PROP_TYPE(TypeTag, VtkOutputFields);
VtkOutputModule<TypeTag> vtkWriter(*problem, *fvGridGeometry, *gridVariables, x, problem->name());
VtkOutputFields::init(vtkWriter); //! Add model specific output fields
+ //add specific output
+ vtkWriter.addField(problem->getExactTemperature(), "temperatureExact");
//! Add model specific output fields
vtkWriter.write(0.0);
+ //add specific output
const auto outputInterval = getParam<int>("Problem.OutputInterval");
// instantiate time loop
@@ -200,22 +203,20 @@ int main(int argc, char** argv) try
// advance to the time loop to the next step
timeLoop->advanceTimeStep();
- // write vtk output
- vtkWriter.write(
- timeLoop->timeStepIndex()==0 ||
- timeLoop->timeStepIndex() % outputInterval == 0 ||
- timeLoop->willBeFinished());
-
+ problem->setTime(timeLoop->time()+timeLoop->timeStepSize());
+ problem->updateExactTemperature(x);
// report statistics of this time step
timeLoop->reportTimeStep();
// set new dt as suggested by newton controller
timeLoop->setTimeStepSize(newtonController->suggestTimeStepSize(timeLoop->timeStepSize()));
- problem->setTime(timeLoop->time()+timeLoop->timeStepSize());
- //add specific output
-// vtkWriter.addField(problem->getExactTemperature(), "exactTemperature");
+ // write vtk output
+ vtkWriter.write(timeLoop->timeStepIndex()==0 ||
+ timeLoop->timeStepIndex() % outputInterval == 0 ||
+ timeLoop->willBeFinished());
+
} while (!timeLoop->finished());
diff --git a/test/porousmediumflow/3p/implicit/test_box3pniconvection.cc b/test/porousmediumflow/3p/implicit/test_box3pniconvection.cc
index d6c2206c46..a180deac45 100644
--- a/test/porousmediumflow/3p/implicit/test_box3pniconvection.cc
+++ b/test/porousmediumflow/3p/implicit/test_box3pniconvection.cc
@@ -146,6 +146,8 @@ int main(int argc, char** argv) try
using VtkOutputFields = typename GET_PROP_TYPE(TypeTag, VtkOutputFields);
VtkOutputModule<TypeTag> vtkWriter(*problem, *fvGridGeometry, *gridVariables, x, problem->name());
VtkOutputFields::init(vtkWriter); //! Add model specific output fields
+ //add specific output
+ vtkWriter.addField(problem->getExactTemperature(), "temperatureExact");
vtkWriter.write(0.0);
const auto outputInterval = getParam<int>("Problem.OutputInterval");
@@ -199,6 +201,9 @@ int main(int argc, char** argv) try
// advance to the time loop to the next step
timeLoop->advanceTimeStep();
+ problem->setTime(timeLoop->time()+timeLoop->timeStepSize());
+ problem->updateExactTemperature(x);
+
// write vtk output
vtkWriter.write(
timeLoop->timeStepIndex()==0 ||
diff --git a/test/porousmediumflow/3p/implicit/test_cc3pniconduction.cc b/test/porousmediumflow/3p/implicit/test_cc3pniconduction.cc
index ced7b6eb03..c7f1f7b690 100644
--- a/test/porousmediumflow/3p/implicit/test_cc3pniconduction.cc
+++ b/test/porousmediumflow/3p/implicit/test_cc3pniconduction.cc
@@ -147,7 +147,11 @@ int main(int argc, char** argv) try
using VtkOutputFields = typename GET_PROP_TYPE(TypeTag, VtkOutputFields);
VtkOutputModule<TypeTag> vtkWriter(*problem, *fvGridGeometry, *gridVariables, x, problem->name());
VtkOutputFields::init(vtkWriter); //! Add model specific output fields
+ vtkWriter.addField(problem->getExactTemperature(), "temperatureExact");
+ //! Add model specific output fields
vtkWriter.write(0.0);
+ //add specific output
+ const auto outputInterval = getParam<int>("Problem.OutputInterval");
// instantiate time loop
auto timeLoop = std::make_shared<TimeLoop<Scalar>>(restartTime, dt, tEnd);
@@ -199,8 +203,8 @@ int main(int argc, char** argv) try
// advance to the time loop to the next step
timeLoop->advanceTimeStep();
- // write vtk output
- vtkWriter.write(timeLoop->time());
+ problem->setTime(timeLoop->time()+timeLoop->timeStepSize());
+ problem->updateExactTemperature(x);
// report statistics of this time step
timeLoop->reportTimeStep();
@@ -208,6 +212,11 @@ int main(int argc, char** argv) try
// set new dt as suggested by newton controller
timeLoop->setTimeStepSize(newtonController->suggestTimeStepSize(timeLoop->timeStepSize()));
+ // write vtk output
+ vtkWriter.write(timeLoop->timeStepIndex()==0 ||
+ timeLoop->timeStepIndex() % outputInterval == 0 ||
+ timeLoop->willBeFinished());
+
} while (!timeLoop->finished());
timeLoop->finalize(leafGridView.comm());
diff --git a/test/porousmediumflow/3p/implicit/test_cc3pniconvection.cc b/test/porousmediumflow/3p/implicit/test_cc3pniconvection.cc
index 3ffb002a0b..5e905d2eef 100644
--- a/test/porousmediumflow/3p/implicit/test_cc3pniconvection.cc
+++ b/test/porousmediumflow/3p/implicit/test_cc3pniconvection.cc
@@ -147,7 +147,11 @@ int main(int argc, char** argv) try
using VtkOutputFields = typename GET_PROP_TYPE(TypeTag, VtkOutputFields);
VtkOutputModule<TypeTag> vtkWriter(*problem, *fvGridGeometry, *gridVariables, x, problem->name());
VtkOutputFields::init(vtkWriter); //! Add model specific output fields
+ vtkWriter.addField(problem->getExactTemperature(), "temperatureExact");
+ //! Add model specific output fields
vtkWriter.write(0.0);
+ //add specific output
+ const auto outputInterval = getParam<int>("Problem.OutputInterval");
// instantiate time loop
auto timeLoop = std::make_shared<TimeLoop<Scalar>>(restartTime, dt, tEnd);
@@ -199,8 +203,8 @@ int main(int argc, char** argv) try
// advance to the time loop to the next step
timeLoop->advanceTimeStep();
- // write vtk output
- vtkWriter.write(timeLoop->time());
+ problem->setTime(timeLoop->time()+timeLoop->timeStepSize());
+ problem->updateExactTemperature(x);
// report statistics of this time step
timeLoop->reportTimeStep();
@@ -208,6 +212,11 @@ int main(int argc, char** argv) try
// set new dt as suggested by newton controller
timeLoop->setTimeStepSize(newtonController->suggestTimeStepSize(timeLoop->timeStepSize()));
+ // write vtk output
+ vtkWriter.write(timeLoop->timeStepIndex()==0 ||
+ timeLoop->timeStepIndex() % outputInterval == 0 ||
+ timeLoop->willBeFinished());
+
} while (!timeLoop->finished());
timeLoop->finalize(leafGridView.comm());
diff --git a/test/porousmediumflow/3p/implicit/test_cc3pniconvection.input b/test/porousmediumflow/3p/implicit/test_cc3pniconvection.input
index 8f4c56c4a4..e189a30b91 100644
--- a/test/porousmediumflow/3p/implicit/test_cc3pniconvection.input
+++ b/test/porousmediumflow/3p/implicit/test_cc3pniconvection.input
@@ -15,3 +15,6 @@ EnableGravity = 0 # disable gravity
[Vtk]
AddVelocity = 1 #Enable velocity output
+
+[Newton]
+EnableResidualCriterion = true
--
GitLab