diff --git a/test/freeflow/navierstokes/angeli/main.cc b/test/freeflow/navierstokes/angeli/main.cc index d0ce122aa0af211456d1d4d31e82cde329ef5be1..918903063d46f5fb1996e845d109949e7463d40a 100644 --- a/test/freeflow/navierstokes/angeli/main.cc +++ b/test/freeflow/navierstokes/angeli/main.cc @@ -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()); diff --git a/test/freeflow/navierstokes/angeli/problem.hh b/test/freeflow/navierstokes/angeli/problem.hh index c793ab1e1205e63d81babf3e21be6804d9219497..f1478831cce51aa230f9e4e2d7e5dbd3f1d866d2 100644 --- a/test/freeflow/navierstokes/angeli/problem.hh +++ b/test/freeflow/navierstokes/angeli/problem.hh @@ -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