Commit 75858e5c authored by Timo Koch's avatar Timo Koch Committed by hanchuan
Browse files

[test][ff][angleli] Cleanup time-dependent functions

parent fffd17dd
......@@ -49,11 +49,10 @@
/*!
* \brief Creates analytical solution.
* Returns a tuple of the analytical solution for the pressure, the velocity and the velocity at the faces
* \param time the time at which to evaluate the analytical solution
* \param problem the problem for which to evaluate the analytical solution
*/
template<class Scalar, class Problem>
auto createAnalyticalSolution(const Scalar time, const Problem& problem)
template<class Problem, class Scalar = double>
auto createAnalyticalSolution(const Problem& problem)
{
const auto& gridGeometry = problem.gridGeometry();
using GridView = typename std::decay_t<decltype(gridGeometry)>::GridView;
......@@ -80,7 +79,7 @@ auto createAnalyticalSolution(const Scalar time, const Problem& problem)
{
auto ccDofIdx = scv.dofIndex();
auto ccDofPosition = scv.dofPosition();
auto analyticalSolutionAtCc = problem.analyticalSolution(ccDofPosition, time);
auto analyticalSolutionAtCc = problem.analyticalSolution(ccDofPosition);
// velocities on faces
for (auto&& scvf : scvfs(fvGeometry))
......@@ -88,7 +87,7 @@ auto createAnalyticalSolution(const Scalar time, const Problem& problem)
const auto faceDofIdx = scvf.dofIndex();
const auto faceDofPosition = scvf.center();
const auto dirIdx = scvf.directionIndex();
const auto analyticalSolutionAtFace = problem.analyticalSolution(faceDofPosition, time);
const auto analyticalSolutionAtFace = problem.analyticalSolution(faceDofPosition);
analyticalVelocityOnFace[faceDofIdx][dirIdx] = analyticalSolutionAtFace[Indices::velocity(dirIdx)];
}
......@@ -148,7 +147,6 @@ int main(int argc, char** argv)
// the problem (initial and boundary conditions)
using Problem = GetPropType<TypeTag, Properties::Problem>;
auto problem = std::make_shared<Problem>(gridGeometry);
problem->updateTimeStepSize(timeLoop->timeStepSize());
// the solution vector
using SolutionVector = GetPropType<TypeTag, Properties::SolutionVector>;
......@@ -168,7 +166,7 @@ int main(int argc, char** argv)
using IOFields = GetPropType<TypeTag, Properties::IOFields>;
IOFields::initOutputModule(vtkWriter); // Add model specific output fields
auto analyticalSolution = createAnalyticalSolution(timeLoop->time(), *problem);
auto analyticalSolution = createAnalyticalSolution(*problem);
vtkWriter.addField(std::get<0>(analyticalSolution), "pressureExact");
vtkWriter.addField(std::get<1>(analyticalSolution), "velocityExact");
vtkWriter.addFaceField(std::get<2>(analyticalSolution), "faceVelocityExact");
......@@ -190,6 +188,8 @@ int main(int argc, char** argv)
const bool printL2Error = getParam<bool>("Problem.PrintL2Error");
timeLoop->start(); do
{
problem->setTime(timeLoop->time() + timeLoop->timeStepSize());
// solve the non-linear system with time step control
nonLinearSolver.solve(x, *timeLoop);
......@@ -204,22 +204,23 @@ int main(int argc, char** argv)
using PrimaryVariables = GetPropType<TypeTag, Properties::PrimaryVariables>;
using L2Error = NavierStokesTestL2Error<Scalar, ModelTraits, PrimaryVariables>;
const auto l2error = L2Error::calculateL2Error(*problem, x);
const auto [l2errorAbs, l2errorRel] = L2Error::calculateL2Error(*problem, x);
const int numCellCenterDofs = gridGeometry->numCellCenterDofs();
const int numFaceDofs = gridGeometry->numFaceDofs();
std::cout << std::setprecision(8) << "** L2 error (abs/rel) for "
<< std::setw(6) << numCellCenterDofs << " cc dofs and " << numFaceDofs << " face dofs (total: " << numCellCenterDofs + numFaceDofs << "): "
<< std::scientific
<< "L2(p) = " << l2error.first[Indices::pressureIdx] << " / " << l2error.second[Indices::pressureIdx]
<< ", L2(vx) = " << l2error.first[Indices::velocityXIdx] << " / " << l2error.second[Indices::velocityXIdx]
<< ", L2(vy) = " << l2error.first[Indices::velocityYIdx] << " / " << l2error.second[Indices::velocityYIdx]
<< "L2(p) = " << l2errorAbs[Indices::pressureIdx] << " / " << l2errorRel[Indices::pressureIdx]
<< ", L2(vx) = " << l2errorAbs[Indices::velocityXIdx] << " / " << l2errorRel[Indices::velocityXIdx]
<< ", L2(vy) = " << l2errorAbs[Indices::velocityYIdx] << " / " << l2errorRel[Indices::velocityYIdx]
<< std::endl;
}
// update the analytical solution
analyticalSolution = createAnalyticalSolution(*problem);
// advance to the time loop to the next step
timeLoop->advanceTimeStep();
problem->updateTime(timeLoop->time());
analyticalSolution = createAnalyticalSolution(timeLoop->time(), *problem);
// write vtk output
vtkWriter.write(timeLoop->time());
......@@ -229,7 +230,6 @@ int main(int argc, char** argv)
// set new dt as suggested by newton solver
timeLoop->setTimeStepSize(nonLinearSolver.suggestTimeStepSize(timeLoop->timeStepSize()));
problem->updateTimeStepSize(timeLoop->timeStepSize());
} while (!timeLoop->finished());
......
......@@ -74,21 +74,19 @@ public:
rho_ = getParam<Scalar>("Component.LiquidDensity", 1.0);
}
/*!
/*!
* \brief Returns the temperature within the domain in [K].
*
* This problem assumes a temperature of 10 degrees Celsius.
* This problem assumes a temperature of 20 degrees Celsius (unused)
*/
Scalar temperature() const
{ return 298.0; }
{ return 293.15; }
// \}
/*!
/*!
* \name Boundary conditions
*/
// \{
/*!
/*!
* \brief Specifies which kind of boundary condition should be
* used for which equation on a given boundary control volume.
*
......@@ -126,29 +124,25 @@ public:
return false;
}
/*!
/*!
* \brief Returns Dirichlet boundary values at a given position.
*
* \param globalPos The global position
*/
PrimaryVariables dirichletAtPos(const GlobalPosition & globalPos) const
{
// use the values of the analytical solution
return analyticalSolution(globalPos, time_+timeStepSize_);
return analyticalSolution(globalPos);
}
/*!
* \brief Returns the analytical solution of the problem at a given time and position.
*
* \param globalPos The global position
* \param time The current simulation time
*/
PrimaryVariables analyticalSolution(const GlobalPosition& globalPos, const Scalar time) const
PrimaryVariables analyticalSolution(const GlobalPosition& globalPos) const
{
const Scalar x = globalPos[0];
const Scalar y = globalPos[1];
const Scalar t = time;
const Scalar t = time_;
PrimaryVariables values;
......@@ -159,54 +153,35 @@ public:
return values;
}
/*!
* \brief Returns the analytical solution of the problem at a given position.
*
* \param globalPos The global position
*/
PrimaryVariables analyticalSolution(const GlobalPosition& globalPos) const
{
return analyticalSolution(globalPos, time_+timeStepSize_);
}
// \}
/*!
/*!
* \name Volume terms
*/
// \{
/*!
/*!
* \brief Evaluates the initial value for a control volume.
*
* \param globalPos The global position
*/
PrimaryVariables initialAtPos(const GlobalPosition &globalPos) const
{
return analyticalSolution(globalPos, 0);
return analyticalSolution(globalPos);
}
/*!
* \brief Updates the time
*/
void updateTime(const Scalar time)
{
time_ = time;
}
// \}
/*!
* \brief Updates the time step size
* \brief Updates the time information
*/
void updateTimeStepSize(const Scalar timeStepSize)
void setTime(Scalar t)
{
timeStepSize_ = timeStepSize;
time_ = t;
}
private:
Scalar kinematicViscosity_;
Scalar rho_;
Scalar kinematicViscosity_, rho_;
Scalar time_ = 0;
Scalar timeStepSize_ = 0;
};
} // end namespace Dumux
......
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment