From 463b1c686cac9e6feb77cf49cae8baa5d02f599c Mon Sep 17 00:00:00 2001 From: DennisGlaeser <dennis.glaeser@iws.uni-stuttgart.de> Date: Tue, 21 Nov 2017 22:08:00 +0100 Subject: [PATCH] [2p][fracture] use unified main file for cc and box --- .../2p/implicit/fracture/CMakeLists.txt | 41 +-- .../2p/implicit/fracture/problem.hh | 12 +- .../implicit/fracture/test_2p_fracture_cc.cc | 14 +- .../implicit/fracture/test_2p_fracture_fv.cc | 234 ++++++++++++++++++ 4 files changed, 263 insertions(+), 38 deletions(-) create mode 100644 test/porousmediumflow/2p/implicit/fracture/test_2p_fracture_fv.cc diff --git a/test/porousmediumflow/2p/implicit/fracture/CMakeLists.txt b/test/porousmediumflow/2p/implicit/fracture/CMakeLists.txt index d3dd26a717..c269eb9c7a 100644 --- a/test/porousmediumflow/2p/implicit/fracture/CMakeLists.txt +++ b/test/porousmediumflow/2p/implicit/fracture/CMakeLists.txt @@ -1,30 +1,31 @@ dune_symlink_to_source_files(FILES "test_fracture.input" "grids") dune_add_test(NAME test_2p_fracture_box - SOURCES test_2p_fracture_box.cc - COMMAND python ${CMAKE_SOURCE_DIR}/bin/testing/runtest.py - --script fuzzy - --files ${CMAKE_SOURCE_DIR}/test/references/fracturebox2p-reference.vtu - ${CMAKE_CURRENT_BINARY_DIR}/fracturebox-00023.vtu - --command "${CMAKE_CURRENT_BINARY_DIR}/test_2p_fracture_box test_fracture.input -Problem.Name fracturebox") + SOURCES test_2p_fracture_fv.cc + COMPILE_DEFINITIONS TYPETAG=FractureBoxTypeTag + COMMAND ${CMAKE_SOURCE_DIR}/bin/testing/runtest.py + CMD_ARGS --script fuzzy + --files ${CMAKE_SOURCE_DIR}/test/references/fracturebox2p-reference.vtu + ${CMAKE_CURRENT_BINARY_DIR}/fracturebox-00023.vtu + --command "${CMAKE_CURRENT_BINARY_DIR}/test_2p_fracture_box test_fracture.input -Problem.Name fracturebox") dune_add_test(NAME test_2p_fracture_tpfa - SOURCES test_2p_fracture_cc.cc - COMPILE_DEFINITIONS TYPETAG=FractureCCTpfaProblem - COMMAND python ${CMAKE_SOURCE_DIR}/bin/testing/runtest.py - --script fuzzy - --files ${CMAKE_SOURCE_DIR}/test/references/fracturecc2p-reference.vtu - ${CMAKE_CURRENT_BINARY_DIR}/fracturetpfa-00029.vtu - --command "${CMAKE_CURRENT_BINARY_DIR}/test_2p_fracture_tpfa test_fracture.input -Problem.Name fracturetpfa") + SOURCES test_2p_fracture_fv.cc + COMPILE_DEFINITIONS TYPETAG=FractureCCTpfaTypeTag + COMMAND ${CMAKE_SOURCE_DIR}/bin/testing/runtest.py + CMD_ARGS --script fuzzy + --files ${CMAKE_SOURCE_DIR}/test/references/fracturecc2p-reference.vtu + ${CMAKE_CURRENT_BINARY_DIR}/fracturetpfa-00029.vtu + --command "${CMAKE_CURRENT_BINARY_DIR}/test_2p_fracture_tpfa test_fracture.input -Problem.Name fracturetpfa") dune_add_test(NAME test_2p_fracture_mpfa - SOURCES test_2p_fracture_cc.cc - COMPILE_DEFINITIONS TYPETAG=FractureCCMpfaProblem - COMMAND python ${CMAKE_SOURCE_DIR}/bin/testing/runtest.py - --script fuzzy - --files ${CMAKE_SOURCE_DIR}/test/references/fractureccmpfa2p-reference.vtu - ${CMAKE_CURRENT_BINARY_DIR}/fracturempfa-00031.vtu - --command "${CMAKE_CURRENT_BINARY_DIR}/test_2p_fracture_mpfa test_fracture.input -Problem.Name fracturempfa") + SOURCES test_2p_fracture_fv.cc + COMPILE_DEFINITIONS TYPETAG=FractureCCMpfaTypeTag + COMMAND ${CMAKE_SOURCE_DIR}/bin/testing/runtest.py + CMD_ARGS --script fuzzy + --files ${CMAKE_SOURCE_DIR}/test/references/fractureccmpfa2p-reference.vtu + ${CMAKE_CURRENT_BINARY_DIR}/fracturempfa-00031.vtu + --command "${CMAKE_CURRENT_BINARY_DIR}/test_2p_fracture_mpfa test_fracture.input -Problem.Name fracturempfa") set(CMAKE_BUILD_TYPE Release) diff --git a/test/porousmediumflow/2p/implicit/fracture/problem.hh b/test/porousmediumflow/2p/implicit/fracture/problem.hh index dca7497e1b..ecc0170e7d 100644 --- a/test/porousmediumflow/2p/implicit/fracture/problem.hh +++ b/test/porousmediumflow/2p/implicit/fracture/problem.hh @@ -46,16 +46,12 @@ class FractureProblem; namespace Properties { NEW_TYPE_TAG(FractureProblem, INHERITS_FROM(TwoP, FractureSpatialParams)); -NEW_TYPE_TAG(FractureBoxProblem, INHERITS_FROM(BoxModel, FractureProblem)); -NEW_TYPE_TAG(FractureCCTpfaProblem, INHERITS_FROM(CCTpfaModel, FractureProblem)); -NEW_TYPE_TAG(FractureCCMpfaProblem, INHERITS_FROM(CCMpfaModel, FractureProblem)); +NEW_TYPE_TAG(FractureBoxTypeTag, INHERITS_FROM(BoxModel, FractureProblem)); +NEW_TYPE_TAG(FractureCCTpfaTypeTag, INHERITS_FROM(CCTpfaModel, FractureProblem)); +NEW_TYPE_TAG(FractureCCMpfaTypeTag, INHERITS_FROM(CCMpfaModel, FractureProblem)); -#if HAVE_DUNE_FOAMGRID +// set the grid property SET_TYPE_PROP(FractureProblem, Grid, Dune::FoamGrid<2, 3>); -#else -#warning "FoamGrid is required for this test!" -SET_TYPE_PROP(FractureProblem, Grid, Dune::YaspGrid<3>); -#endif // Set the problem property SET_TYPE_PROP(FractureProblem, Problem, Dumux::FractureProblem<TypeTag>); diff --git a/test/porousmediumflow/2p/implicit/fracture/test_2p_fracture_cc.cc b/test/porousmediumflow/2p/implicit/fracture/test_2p_fracture_cc.cc index d5117ca6a9..3bb9db15fe 100644 --- a/test/porousmediumflow/2p/implicit/fracture/test_2p_fracture_cc.cc +++ b/test/porousmediumflow/2p/implicit/fracture/test_2p_fracture_cc.cc @@ -67,7 +67,6 @@ void usage(const char *progName, const std::string &errorMsg) //////////////////////// int main(int argc, char** argv) try { -#if HAVE_DUNE_FOAMGRID using namespace Dumux; // define the type tag for this problem @@ -106,7 +105,7 @@ int main(int argc, char** argv) try // the solution vector using SolutionVector = typename GET_PROP_TYPE(TypeTag, SolutionVector); - SolutionVector x(leafGridView.size(0)); + SolutionVector x(fvGridGeometry->numDofs()); problem->applyInitialSolution(x); auto xOld = x; @@ -143,11 +142,11 @@ int main(int argc, char** argv) try // the linear solver using LinearSolver = AMGBackend<TypeTag>; - auto linearSolver = std::make_shared<LinearSolver>(leafGridView, fvGridGeometry->elementMapper()); + auto linearSolver = std::make_shared<LinearSolver>(leafGridView, fvGridGeometry->dofMapper()); // the non-linear solver - using NewtonController = NewtonController<TypeTag>; - using NewtonMethod = NewtonMethod<TypeTag, NewtonController, Assembler, LinearSolver>; + using NewtonController = Dumux::NewtonController<TypeTag>; + using NewtonMethod = Dumux::NewtonMethod<NewtonController, Assembler, LinearSolver>; auto newtonController = std::make_shared<NewtonController>(leafGridView.comm(), timeLoop); NewtonMethod nonLinearSolver(newtonController, assembler, linearSolver); @@ -208,11 +207,6 @@ int main(int argc, char** argv) try } return 0; -#else -#warning External grid module dune-foamgrid needed to run this example. - std::cerr << "Test skipped, it needs dune-foamgrid!" << std::endl; - return 77; -#endif } // end main catch (Dumux::ParameterException &e) { diff --git a/test/porousmediumflow/2p/implicit/fracture/test_2p_fracture_fv.cc b/test/porousmediumflow/2p/implicit/fracture/test_2p_fracture_fv.cc new file mode 100644 index 0000000000..3bb9db15fe --- /dev/null +++ b/test/porousmediumflow/2p/implicit/fracture/test_2p_fracture_fv.cc @@ -0,0 +1,234 @@ +// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- +// vi: set et ts=4 sw=4 sts=4: +/***************************************************************************** + * See the file COPYING for full copying permissions. * + * * + * This program is free software: you can redistribute it and/or modify * + * it under the terms of the GNU General Public License as published by * + * the Free Software Foundation, either version 2 of the License, or * + * (at your option) any later version. * + * * + * This program is distributed in the hope that it will be useful, * + * but WITHOUT ANY WARRANTY; without even the implied warranty of * + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the * + * GNU General Public License for more details. * + * * + * You should have received a copy of the GNU General Public License * + * along with this program. If not, see <http://www.gnu.org/licenses/>. * + *****************************************************************************/ +/*! + * \file + * + * \brief Fracture 2d in 3d test for the two-phase box model + */ +#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 "problem.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/amgbackend.hh> +#include <dumux/nonlinear/newtonmethod.hh> +#include <dumux/nonlinear/newtoncontroller.hh> + +#include <dumux/assembly/fvassembler.hh> +#include <dumux/assembly/diffmethod.hh> + +#include <dumux/discretization/methods.hh> + +#include <dumux/io/vtkoutputmodule.hh> + +/*! + * \brief Provides an interface for customizing error messages associated with + * reading in parameters. + * + * \param progName The name of the program, that was tried to be started. + * \param errorMsg The error message that was issued by the start function. + * Comprises the thing that went wrong and a general help message. + */ +void usage(const char *progName, const std::string &errorMsg) +{} + +//////////////////////// +// the main function +//////////////////////// +int main(int argc, char** argv) try +{ + using namespace Dumux; + + // define the type tag for this problem + using TypeTag = TTAG(TYPETAG); + + // 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 SolutionVector = typename GET_PROP_TYPE(TypeTag, SolutionVector); + SolutionVector x(fvGridGeometry->numDofs()); + 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 = AMGBackend<TypeTag>; + auto linearSolver = std::make_shared<LinearSolver>(leafGridView, fvGridGeometry->dofMapper()); + + // the non-linear solver + using NewtonController = Dumux::NewtonController<TypeTag>; + using NewtonMethod = Dumux::NewtonMethod<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; +} // 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 (...) +{ + std::cerr << "Unknown exception thrown! ---> Abort!" << std::endl; + return 4; +} -- GitLab