From 926e2a40f6ac1a9542ff388d101a8d7d4eab691f Mon Sep 17 00:00:00 2001 From: Kilian Weishaupt <kilian.weishaupt@iws.uni-stuttgart.de> Date: Wed, 8 Nov 2017 16:40:49 +0100 Subject: [PATCH] [freeflow] Port test_kovasznay to new structure of next --- .../staggered/kovasznaytestproblem.hh | 191 ++++++++-------- test/freeflow/staggered/test_kovasznay.cc | 212 +++++++++++++++++- test/freeflow/staggered/test_kovasznay.input | 11 +- 3 files changed, 310 insertions(+), 104 deletions(-) diff --git a/test/freeflow/staggered/kovasznaytestproblem.hh b/test/freeflow/staggered/kovasznaytestproblem.hh index 03f7b6290d..30371b3564 100644 --- a/test/freeflow/staggered/kovasznaytestproblem.hh +++ b/test/freeflow/staggered/kovasznaytestproblem.hh @@ -24,13 +24,15 @@ #ifndef DUMUX_KOVASZNAY_TEST_PROBLEM_HH #define DUMUX_KOVASZNAY_TEST_PROBLEM_HH -#include <dumux/implicit/staggered/properties.hh> -#include <dumux/freeflow/staggered/model.hh> -#include <dumux/implicit/problem.hh> +#include <dumux/freeflow/staggered/problem.hh> +#include <dumux/discretization/staggered/properties.hh> #include <dumux/material/components/simpleh2o.hh> #include <dumux/material/fluidsystems/liquidphase.hh> #include <dumux/material/components/constant.hh> +#include <dumux/discretization/staggered/properties.hh> +#include <dumux/freeflow/staggered/properties.hh> + // solve Navier-Stokes equations #define ENABLE_NAVIERSTOKES 1 @@ -70,9 +72,6 @@ SET_BOOL_PROP(KovasznayTestProblem, EnableFVGridGeometryCache, true); SET_BOOL_PROP(KovasznayTestProblem, EnableGlobalFluxVariablesCache, true); SET_BOOL_PROP(KovasznayTestProblem, EnableGlobalVolumeVariablesCache, true); -// Enable gravity -SET_BOOL_PROP(KovasznayTestProblem, ProblemEnableGravity, true); - #if ENABLE_NAVIERSTOKES SET_BOOL_PROP(KovasznayTestProblem, EnableInertiaTerms, true); #else @@ -88,13 +87,13 @@ SET_BOOL_PROP(KovasznayTestProblem, EnableInertiaTerms, false); template <class TypeTag> class KovasznayTestProblem : public NavierStokesProblem<TypeTag> { - typedef NavierStokesProblem<TypeTag> ParentType; + using ParentType = NavierStokesProblem<TypeTag>; - typedef typename GET_PROP_TYPE(TypeTag, GridView) GridView; - typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar; + using GridView = typename GET_PROP_TYPE(TypeTag, GridView); + using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar); // copy some indices for convenience - typedef typename GET_PROP_TYPE(TypeTag, Indices) Indices; + using Indices = typename GET_PROP_TYPE(TypeTag, Indices); enum { // Grid and world dimension dim = GridView::dimension, @@ -110,16 +109,16 @@ class KovasznayTestProblem : public NavierStokesProblem<TypeTag> velocityYIdx = Indices::velocityYIdx }; - typedef typename GET_PROP_TYPE(TypeTag, BoundaryTypes) BoundaryTypes; - typedef typename GET_PROP_TYPE(TypeTag, TimeManager) TimeManager; + using BoundaryTypes = typename GET_PROP_TYPE(TypeTag, BoundaryTypes); - typedef typename GridView::template Codim<0>::Entity Element; - typedef typename GridView::Intersection Intersection; + using Element = typename GridView::template Codim<0>::Entity; + using Intersection = typename GridView::Intersection; - typedef typename GET_PROP_TYPE(TypeTag, FVElementGeometry) FVElementGeometry; - typedef typename GET_PROP_TYPE(TypeTag, SubControlVolume) SubControlVolume; + using FVElementGeometry = typename GET_PROP_TYPE(TypeTag, FVElementGeometry); + using SubControlVolume = typename GET_PROP_TYPE(TypeTag, SubControlVolume); + using FVGridGeometry = typename GET_PROP_TYPE(TypeTag, FVGridGeometry); - typedef Dune::FieldVector<Scalar, dimWorld> GlobalPosition; + using GlobalPosition = Dune::FieldVector<Scalar, dimWorld>; using CellCenterPrimaryVariables = typename GET_PROP_TYPE(TypeTag, CellCenterPrimaryVariables); using FacePrimaryVariables = typename GET_PROP_TYPE(TypeTag, FacePrimaryVariables); @@ -127,36 +126,29 @@ class KovasznayTestProblem : public NavierStokesProblem<TypeTag> using BoundaryValues = typename GET_PROP_TYPE(TypeTag, BoundaryValues); using InitialValues = typename GET_PROP_TYPE(TypeTag, BoundaryValues); using SourceValues = typename GET_PROP_TYPE(TypeTag, BoundaryValues); + using SolutionVector = typename GET_PROP_TYPE(TypeTag, SolutionVector); using DofTypeIndices = typename GET_PROP(TypeTag, DofTypeIndices); typename DofTypeIndices::CellCenterIdx cellCenterIdx; typename DofTypeIndices::FaceIdx faceIdx; public: - KovasznayTestProblem(TimeManager &timeManager, const GridView &gridView) - : ParentType(timeManager, gridView), eps_(1e-6) + KovasznayTestProblem(std::shared_ptr<const FVGridGeometry> fvGridGeometry) + : ParentType(fvGridGeometry), eps_(1e-6) { - name_ = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, - std::string, - Problem, - Name); - - printL2Error_ = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, - bool, - Problem, - PrintL2Error); + printL2Error_ = getParam<bool>("Problem.PrintL2Error"); - kinematicViscosity_ = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, Scalar, Problem, LiquidKinematicViscosity); + kinematicViscosity_ = getParam<Scalar>("Component.LiquidKinematicViscosity", 1.0); Scalar reynoldsNumber = 1.0 / kinematicViscosity_; lambda_ = 0.5 * reynoldsNumber - std::sqrt(reynoldsNumber * reynoldsNumber * 0.25 + 4.0 * M_PI * M_PI); using CellArray = std::array<unsigned int, dimWorld>; - const CellArray numCells = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, - CellArray, - Grid, - Cells); - cellSizeX_ = this->bBoxMax()[0] / numCells[0]; + const auto numCells = getParam<CellArray>("Grid.Cells"); + + cellSizeX_ = this->fvGridGeometry().bBoxMax()[0] / numCells[0]; + + createAnalyticalSolution_(); } /*! @@ -164,28 +156,18 @@ public: */ // \{ - /*! - * \brief The problem name. - * - * This is used as a prefix for files generated by the simulation. - */ - std::string name() const - { - return name_; - } - bool shouldWriteRestartFile() const { return false; } - void postTimeStep() const + void postTimeStep(const SolutionVector& curSol) const { if(printL2Error_) { - const auto l2error = calculateL2Error(); - const int numCellCenterDofs = this->model().numCellCenterDofs(); - const int numFaceDofs = this->model().numFaceDofs(); + const auto l2error = calculateL2Error(curSol); + const int numCellCenterDofs = this->fvGridGeometry().gridView().size(0); + const int numFaceDofs = this->fvGridGeometry().gridView().size(1); 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 @@ -294,56 +276,16 @@ public: return values; } - /*! - * \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 - { - auto& pressureExact = outputModule.createScalarField("pressureExact", 0); - auto& velocityExact = outputModule.createVectorField("velocityExact", 0); - - auto& scalarFaceVelocityExact = outputModule.createFaceScalarField("scalarFaceVelocityExact"); - auto& vectorFaceVelocityExact = outputModule.createFaceVectorField("vectorFaceVelocityExact"); - - for (const auto& element : elements(this->gridView())) - { - auto fvGeometry = localView(this->model().fvGridGeometry()); - fvGeometry.bindElement(element); - for (auto&& scv : scvs(fvGeometry)) - { - auto ccDofIdx = scv.dofIndex(); - auto ccDofPosition = scv.dofPosition(); - auto analyticalSolutionAtCc = dirichletAtPos(ccDofPosition); - - GlobalPosition velocityVector(0.0); - for (auto&& scvf : scvfs(fvGeometry)) - { - auto faceDofIdx = scvf.dofIndex(); - auto faceDofPosition = scvf.center(); - auto dirIdx = scvf.directionIndex(); - auto analyticalSolutionAtFace = dirichletAtPos(faceDofPosition); - scalarFaceVelocityExact[faceDofIdx] = analyticalSolutionAtFace[faceIdx][dirIdx]; - - GlobalPosition tmp(0.0); - tmp[dirIdx] = analyticalSolutionAtFace[faceIdx][dirIdx]; - vectorFaceVelocityExact[faceDofIdx] = std::move(tmp); - } - pressureExact[ccDofIdx] = analyticalSolutionAtCc[pressureIdx]; - velocityExact[ccDofIdx] = analyticalSolutionAtCc[faceIdx]; - } - } - } /*! * \brief Calculate the L2 error between the analytical solution and the numerical approximation. * */ - auto calculateL2Error() const + auto calculateL2Error(const SolutionVector& curSol) const { BoundaryValues sumError(0.0), sumReference(0.0), l2NormAbs(0.0), l2NormRel(0.0); - const int numFaceDofs = this->model().numFaceDofs(); + const int numFaceDofs = this->fvGridGeometry().gridView().size(1); std::vector<Scalar> staggeredVolume(numFaceDofs); std::vector<Scalar> errorVelocity(numFaceDofs); @@ -352,9 +294,9 @@ public: Scalar totalVolume = 0.0; - 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)) @@ -363,7 +305,7 @@ public: const auto dofIdxCellCenter = scv.dofIndex(); const auto& posCellCenter = scv.dofPosition(); const auto analyticalSolutionCellCenter = dirichletAtPos(posCellCenter)[cellCenterIdx]; - const auto numericalSolutionCellCenter = this->model().curSol()[cellCenterIdx][dofIdxCellCenter]; + const auto numericalSolutionCellCenter = curSol[cellCenterIdx][dofIdxCellCenter]; sumError[cellCenterIdx] += squaredDiff_(analyticalSolutionCellCenter, numericalSolutionCellCenter) * scv.volume(); sumReference[cellCenterIdx] += analyticalSolutionCellCenter * analyticalSolutionCellCenter * scv.volume(); totalVolume += scv.volume(); @@ -374,7 +316,7 @@ public: const int dofIdxFace = scvf.dofIndex(); const int dirIdx = scvf.directionIndex(); const auto analyticalSolutionFace = dirichletAtPos(scvf.center())[faceIdx][dirIdx]; - const auto numericalSolutionFace = this->model().curSol()[faceIdx][dofIdxFace][momentumBalanceIdx]; + const auto numericalSolutionFace = curSol[faceIdx][dofIdxFace][momentumBalanceIdx]; directionIndex[dofIdxFace] = dirIdx; errorVelocity[dofIdxFace] = squaredDiff_(analyticalSolutionFace, numericalSolutionFace); velocityReference[dofIdxFace] = squaredDiff_(analyticalSolutionFace, 0.0); @@ -407,7 +349,63 @@ public: return std::make_pair(l2NormAbs, l2NormRel); } + /*! + * \brief Returns the analytical solution for the pressure + */ + auto& getAnalyticalPressureSolution() const + { + return analyticalPressure_; + } + + /*! + * \brief Returns the analytical solution for the velocity + */ + auto& getAnalyticalVelocitySolution() const + { + return analyticalVelocity_; + } + private: + + /*! + * \brief Adds additional VTK output data to the VTKWriter. Function is called by the output module on every write. + */ + void createAnalyticalSolution_() + { + analyticalPressure_.resize(this->fvGridGeometry().gridView().size(0)); + analyticalVelocity_.resize(this->fvGridGeometry().gridView().size(0)); + + + for (const auto& element : elements(this->fvGridGeometry().gridView())) + { + auto fvGeometry = localView(this->fvGridGeometry()); + fvGeometry.bindElement(element); + for (auto&& scv : scvs(fvGeometry)) + { + auto ccDofIdx = scv.dofIndex(); + auto ccDofPosition = scv.dofPosition(); + auto analyticalSolutionAtCc = analyticalSolution(ccDofPosition); + + // TODO: velocities on faces + // GlobalPosition velocityVector(0.0); + // for (auto&& scvf : scvfs(fvGeometry)) + // { + // auto faceDofIdx = scvf.dofIndex(); + // auto faceDofPosition = scvf.center(); + // auto dirIdx = scvf.directionIndex(); + // auto analyticalSolutionAtFace = analyticalSolution(faceDofPosition); + // // scalarFaceVelocityExact[faceDofIdx] = analyticalSolutionAtFace[faceIdx][dirIdx]; + // + // GlobalPosition tmp(0.0); + // tmp[dirIdx] = analyticalSolutionAtFace[faceIdx][dirIdx]; + // vectorFaceVelocityExact[faceDofIdx] = std::move(tmp); + // } + analyticalPressure_[ccDofIdx] = analyticalSolutionAtCc[pressureIdx]; + analyticalVelocity_[ccDofIdx] = analyticalSolutionAtCc[faceIdx]; + } + } + } + template<class T> T squaredDiff_(const T& a, const T& b) const { @@ -416,16 +414,17 @@ private: bool isLowerLeftCell_(const GlobalPosition& globalPos) const { - return globalPos[0] < (this->bBoxMin()[0] + 0.5*cellSizeX_ + eps_); + return globalPos[0] < (this->fvGridGeometry().bBoxMin()[0] + 0.5*cellSizeX_ + eps_); } Scalar eps_; Scalar cellSizeX_; - std::string name_; Scalar kinematicViscosity_; Scalar lambda_; bool printL2Error_; + std::vector<Scalar> analyticalPressure_; + std::vector<GlobalPosition> analyticalVelocity_; }; } //end namespace diff --git a/test/freeflow/staggered/test_kovasznay.cc b/test/freeflow/staggered/test_kovasznay.cc index 67d5014aec..994d63d603 100644 --- a/test/freeflow/staggered/test_kovasznay.cc +++ b/test/freeflow/staggered/test_kovasznay.cc @@ -21,9 +21,37 @@ * * \brief Test for the staggered grid Stokes model (Donea et al., 2003) */ -#include <config.h> + #include <config.h> + + #include <ctime> + #include <iostream> + + #include <dune/common/parallel/mpihelper.hh> + #include <dune/common/timer.hh> + #include <dune/grid/io/file/dgfparser/dgfexception.hh> + #include <dune/grid/io/file/vtk.hh> + #include <dune/istl/io.hh> + + #include "kovasznaytestproblem.hh" -#include <dumux/common/start.hh> + +#include <dumux/common/propertysystem.hh> +#include <dumux/common/parameters.hh> +#include <dumux/common/valgrind.hh> +#include <dumux/common/dumuxmessage.hh> +#include <dumux/common/defaultusagemessage.hh> + +#include <dumux/linear/seqsolverbackend.hh> +#include <dumux/nonlinear/newtonmethod.hh> +#include <dumux/nonlinear/newtoncontroller.hh> + +#include <dumux/assembly/staggeredfvassembler.hh> +#include <dumux/assembly/diffmethod.hh> +#include <dumux/implicit/staggered/newtoncontroller.hh> + +#include <dumux/discretization/methods.hh> + +#include <dumux/io/vtkoutputmodule.hh> /*! * \brief Provides an interface for customizing error messages associated with @@ -57,8 +85,182 @@ void usage(const char *progName, const std::string &errorMsg) } } -int main(int argc, char** argv) +int main(int argc, char** argv) try +{ + using namespace Dumux; + + // define the type tag for this problem + using TypeTag = TTAG(KovasznayTestProblem); + + // initialize MPI, finalize is done automatically on exit + const auto& mpiHelper = Dune::MPIHelper::instance(argc, argv); + + // print dumux start message + if (mpiHelper.rank() == 0) + DumuxMessage::print(/*firstCall=*/true); + + // parse command line arguments and input file + Parameters::init(argc, argv, usage); + + // try to create a grid (from the given grid file or the input file) + using GridCreator = typename GET_PROP_TYPE(TypeTag, GridCreator); + GridCreator::makeGrid(); + GridCreator::loadBalance(); + + //////////////////////////////////////////////////////////// + // run instationary non-linear problem on this grid + //////////////////////////////////////////////////////////// + + // we compute on the leaf grid view + const auto& leafGridView = GridCreator::grid().leafGridView(); + + // create the finite volume grid geometry + using FVGridGeometry = typename GET_PROP_TYPE(TypeTag, FVGridGeometry); + auto fvGridGeometry = std::make_shared<FVGridGeometry>(leafGridView); + fvGridGeometry->update(); + + // the problem (initial and boundary conditions) + using Problem = typename GET_PROP_TYPE(TypeTag, Problem); + auto problem = std::make_shared<Problem>(fvGridGeometry); + + // the solution vector + using GridView = typename GET_PROP_TYPE(TypeTag, GridView); + using SolutionVector = typename GET_PROP_TYPE(TypeTag, SolutionVector); + using DofTypeIndices = typename GET_PROP(TypeTag, DofTypeIndices); + typename DofTypeIndices::CellCenterIdx cellCenterIdx; + typename DofTypeIndices::FaceIdx faceIdx; + const auto numDofsCellCenter = leafGridView.size(0); + const auto numDofsFace = leafGridView.size(1); + SolutionVector x; + x[cellCenterIdx].resize(numDofsCellCenter); + x[faceIdx].resize(numDofsFace); + problem->applyInitialSolution(x); + auto xOld = x; + + // the grid variables + using GridVariables = typename GET_PROP_TYPE(TypeTag, GridVariables); + auto gridVariables = std::make_shared<GridVariables>(problem, fvGridGeometry); + gridVariables->init(x, xOld); + + // get some time loop parameters + using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar); + const auto tEnd = getParam<Scalar>("TimeLoop.TEnd"); + const auto maxDivisions = getParam<int>("TimeLoop.MaxTimeStepDivisions"); + const auto maxDt = getParam<Scalar>("TimeLoop.MaxTimeStepSize"); + auto dt = getParam<Scalar>("TimeLoop.DtInitial"); + + // check if we are about to restart a previously interrupted simulation + Scalar restartTime = 0; + if (Parameters::getTree().hasKey("Restart") || Parameters::getTree().hasKey("TimeLoop.Restart")) + restartTime = getParam<Scalar>("TimeLoop.Restart"); + + // intialize the vtk output module + 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->getAnalyticalPressureSolution(), "pressureExact", 1); + vtkWriter.addField(problem->getAnalyticalVelocitySolution(), "velocityExact", GridView::dimensionworld); + vtkWriter.write(0.0); + + // instantiate time loop + auto timeLoop = std::make_shared<TimeLoop<Scalar>>(restartTime, dt, tEnd); + timeLoop->setMaxTimeStepSize(maxDt); + + // the assembler with time loop for instationary problem + using Assembler = StaggeredFVAssembler<TypeTag, DiffMethod::numeric>; + auto assembler = std::make_shared<Assembler>(problem, fvGridGeometry, gridVariables, timeLoop); + + // the linear solver + using LinearSolver = Dumux::UMFPackBackend<TypeTag>; + auto linearSolver = std::make_shared<LinearSolver>(); + + // the non-linear solver + using NewtonController = StaggeredNewtonController<TypeTag>; + using NewtonMethod = Dumux::NewtonMethod<TypeTag, NewtonController, Assembler, LinearSolver>; + auto newtonController = std::make_shared<NewtonController>(leafGridView.comm(), timeLoop); + NewtonMethod nonLinearSolver(newtonController, assembler, linearSolver); + + // time loop + timeLoop->start(); do + { + // set previous solution for storage evaluations + assembler->setPreviousSolution(xOld); + + // try solving the non-linear system + for (int i = 0; i < maxDivisions; ++i) + { + // linearize & solve + auto converged = nonLinearSolver.solve(x); + + if (converged) + break; + + if (!converged && i == maxDivisions-1) + DUNE_THROW(Dune::MathError, + "Newton solver didn't converge after " + << maxDivisions + << " time-step divisions. dt=" + << timeLoop->timeStepSize() + << ".\nThe solutions of the current and the previous time steps " + << "have been saved to restart files."); + } + + // make the new solution the old solution + xOld = x; + gridVariables->advanceTimeStep(); + + // advance to the time loop to the next step + timeLoop->advanceTimeStep(); + + problem->postTimeStep(x); + + // write vtk output + vtkWriter.write(timeLoop->time()); + + // report statistics of this time step + timeLoop->reportTimeStep(); + + // set new dt as suggested by newton controller + timeLoop->setTimeStepSize(newtonController->suggestTimeStepSize(timeLoop->timeStepSize())); + + } while (!timeLoop->finished()); + + timeLoop->finalize(leafGridView.comm()); + + //////////////////////////////////////////////////////////// + // finalize, print dumux message to say goodbye + //////////////////////////////////////////////////////////// + + // print dumux end message + if (mpiHelper.rank() == 0) + { + Parameters::print(); + DumuxMessage::print(/*firstCall=*/false); + } + + return 0; +} // end main +catch (Dumux::ParameterException &e) +{ + std::cerr << std::endl << e << " ---> Abort!" << std::endl; + return 1; +} +catch (Dune::DGFException & e) +{ + std::cerr << "DGF exception thrown (" << e << + "). Most likely, the DGF file name is wrong " + "or the DGF file is corrupted, " + "e.g. missing hash at end of file or wrong number (dimensions) of entries." + << " ---> Abort!" << std::endl; + return 2; +} +catch (Dune::Exception &e) +{ + std::cerr << "Dune reported error: " << e << " ---> Abort!" << std::endl; + return 3; +} +catch (...) { - typedef TTAG(KovasznayTestProblem) ProblemTypeTag; - return Dumux::start<ProblemTypeTag>(argc, argv, usage); + std::cerr << "Unknown exception thrown! ---> Abort!" << std::endl; + return 4; } diff --git a/test/freeflow/staggered/test_kovasznay.input b/test/freeflow/staggered/test_kovasznay.input index 02e71ddbcf..4372ca4e49 100644 --- a/test/freeflow/staggered/test_kovasznay.input +++ b/test/freeflow/staggered/test_kovasznay.input @@ -1,4 +1,4 @@ -[TimeManager] +[TimeLoop] DtInitial = 1 # [s] TEnd = 2 # [s] @@ -9,11 +9,16 @@ Cells = 50 50 [Problem] Name = test_kovasznay # name passed to the output routines -LiquidDensity = 1 EnableGravity = false -LiquidKinematicViscosity = 0.025 PrintL2Error = false +[Component] +LiquidDensity = 1 +LiquidKinematicViscosity = 0.025 + [ Newton ] MaxSteps = 10 MaxRelativeShift = 1e-5 + +[Vtk] +AddVelocity = true -- GitLab