From cf2ddada5a785c17d81cd97e4dc0e52bb3d59f23 Mon Sep 17 00:00:00 2001 From: hommel <johannes.hommel@iws.uni-stuttgart.de> Date: Wed, 25 Oct 2017 12:12:10 +0200 Subject: [PATCH] [2p] test with new .cc file structure --- .../2p/implicit/lensproblem.hh | 7 +- .../2p/implicit/test_box2p.cc | 194 +++++++++++++++++- .../porousmediumflow/2p/implicit/test_cc2p.cc | 193 ++++++++++++++++- 3 files changed, 381 insertions(+), 13 deletions(-) diff --git a/test/porousmediumflow/2p/implicit/lensproblem.hh b/test/porousmediumflow/2p/implicit/lensproblem.hh index 4897635d0d..8d7a63e0da 100644 --- a/test/porousmediumflow/2p/implicit/lensproblem.hh +++ b/test/porousmediumflow/2p/implicit/lensproblem.hh @@ -29,7 +29,7 @@ #include <dumux/material/components/simpleh2o.hh> #include <dumux/material/components/dnapl.hh> #include <dumux/porousmediumflow/2p/implicit/model.hh> -#include <dumux/porousmediumflow/implicit/problem.hh> +#include <dumux/porousmediumflow/problem.hh> #include <dumux/implicit/cellcentered/tpfa/properties.hh> #include <dumux/porousmediumflow/2p/implicit/gridadaptindicator.hh> #include <dumux/porousmediumflow/2p/implicit/adaptionhelper.hh> @@ -215,10 +215,7 @@ public: { temperature_ = 273.15 + 20; // -> 20°C - name_ = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, - std::string, - Problem, - Name); + name_ = getParam<std::string>("Problem.Name"); } /*! diff --git a/test/porousmediumflow/2p/implicit/test_box2p.cc b/test/porousmediumflow/2p/implicit/test_box2p.cc index 2d6beebc94..e0666e462e 100644 --- a/test/porousmediumflow/2p/implicit/test_box2p.cc +++ b/test/porousmediumflow/2p/implicit/test_box2p.cc @@ -23,7 +23,29 @@ */ #include <config.h> #include "lensproblem.hh" -#include <dumux/common/start.hh> + +#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 <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/common/parameterparser.hh> + +#include <dumux/linear/seqsolverbackend.hh> +#include <dumux/nonlinear/newtonmethod.hh> + +#include <dumux/assembly/fvassembler.hh> + +#include <dumux/io/vtkoutputmodule.hh> /*! * \brief Provides an interface for customizing error messages associated with @@ -57,8 +79,172 @@ void usage(const char *progName, const std::string &errorMsg) //////////////////////// // the main function //////////////////////// -int main(int argc, char** argv) +int main(int argc, char** argv) try { - typedef TTAG(LensBoxProblem) TypeTag; - return Dumux::start<TypeTag>(argc, argv, usage); + using namespace Dumux; + + // define the type tag for this problem + using TypeTag = TTAG(LensBoxProblem); + + // 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(Parameters::getTree()); + 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); + SolutionVector x(leafGridView.size(GridView::dimension)); + 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.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 = FVAssembler<TypeTag, DiffMethod::numeric>; + auto assembler = std::make_shared<Assembler>(problem, fvGridGeometry, gridVariables, timeLoop); + + // the linear solver + using LinearSolver = typename GET_PROP_TYPE(TypeTag, LinearSolver); + auto linearSolver = std::make_shared<LinearSolver>(leafGridView, fvGridGeometry->vertexMapper()); + + // the non-linear solver + using NewtonController = Dumux::NewtonController<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(); + + // 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; +} + +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 (...) +{ + std::cerr << "Unknown exception thrown! ---> Abort!" << std::endl; + return 4; +} \ No newline at end of file diff --git a/test/porousmediumflow/2p/implicit/test_cc2p.cc b/test/porousmediumflow/2p/implicit/test_cc2p.cc index 2cb7906ea8..8f3b8c5633 100644 --- a/test/porousmediumflow/2p/implicit/test_cc2p.cc +++ b/test/porousmediumflow/2p/implicit/test_cc2p.cc @@ -23,7 +23,29 @@ */ #include <config.h> #include "lensproblem.hh" -#include <dumux/common/start.hh> + +#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 <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/common/parameterparser.hh> + +#include <dumux/linear/seqsolverbackend.hh> +#include <dumux/nonlinear/newtonmethod.hh> + +#include <dumux/assembly/fvassembler.hh> + +#include <dumux/io/vtkoutputmodule.hh> /*! * \brief Provides an interface for customizing error messages associated with @@ -57,8 +79,171 @@ void usage(const char *progName, const std::string &errorMsg) //////////////////////// // the main function //////////////////////// -int main(int argc, char** argv) +int main(int argc, char** argv) try { - typedef TTAG(LensCCProblem) TypeTag; - return Dumux::start<TypeTag>(argc, argv, usage); + using namespace Dumux; + + // define the type tag for this problem + using TypeTag = TTAG(LensCCProblem); + + // 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(Parameters::getTree()); + 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 SolutionVector = typename GET_PROP_TYPE(TypeTag, SolutionVector); + SolutionVector x(leafGridView.size(0)); + 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.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 = FVAssembler<TypeTag, DiffMethod::numeric>; + auto assembler = std::make_shared<Assembler>(problem, fvGridGeometry, gridVariables, timeLoop); + + // the linear solver + using LinearSolver = typename GET_PROP_TYPE(TypeTag, LinearSolver); + auto linearSolver = std::make_shared<LinearSolver>(leafGridView, fvGridGeometry->elementMapper()); + + // the non-linear solver + using NewtonController = typename GET_PROP_TYPE(TypeTag, NewtonController); + 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(); + + // 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; +} + +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 (...) +{ + std::cerr << "Unknown exception thrown! ---> Abort!" << std::endl; + return 4; +} \ No newline at end of file -- GitLab