From f669136297c0b8a7c8bae75ec59f6760834d4390 Mon Sep 17 00:00:00 2001 From: Bernd Flemisch Date: Tue, 25 May 2010 13:27:16 +0000 Subject: [PATCH] adapted decoupled tutorial to the new material system --- examples/tutorial_decoupled.cc | 166 +++++----- examples/tutorialproblem_coupled.hh | 2 +- examples/tutorialproblem_decoupled.hh | 302 +++++++++++++----- examples/tutorialspatialparameters_coupled.hh | 4 +- .../tutorialspatialparameters_decoupled.hh | 98 ++++++ 5 files changed, 398 insertions(+), 174 deletions(-) create mode 100644 examples/tutorialspatialparameters_decoupled.hh diff --git a/examples/tutorial_decoupled.cc b/examples/tutorial_decoupled.cc index 07197ec37..bfe7c6b35 100644 --- a/examples/tutorial_decoupled.cc +++ b/examples/tutorial_decoupled.cc @@ -1,6 +1,8 @@ // $Id$ /***************************************************************************** - * Copyright (C) 2008-2009 by Markus Wolff * + * Copyright (C) 20010 by Markus Wolff * + * Copyright (C) 2007-2008 by Bernd Flemisch * + * Copyright (C) 2008-2009 by Andreas Lauser * * Institute of Hydraulic Engineering * * University of Stuttgart, Germany * * email: .@iws.uni-stuttgart.de * @@ -14,104 +16,98 @@ * This program is distributed WITHOUT ANY WARRANTY. * *****************************************************************************/ #include "config.h" -#include -#include -#include /*@\label{tutorial-decoupled:include-begin}@*/ -#include -#include -#include -#include "dumux/fractionalflow/variableclass2p.hh" -#include "dumux/fractionalflow/define2pmodel.hh" -#include "dumux/material/fluids/water.hh" -#include "dumux/material/fluids/lowviscosityoil.hh" -#include "tutorial_soilproperties_decoupled.hh" -#include "dumux/material/twophaserelations.hh" -#include "tutorialproblem_decoupled.hh" -#include "dumux/diffusion/fv/fvvelocity2p.hh" -#include "dumux/transport/fv/fvsaturation2p.hh" -#include "dumux/fractionalflow/impes/impes.hh" -#include "dumux/timedisc/timeloop.hh" /*@\label{tutorial-decoupled:include-end}@*/ +#include "tutorialproblem_decoupled.hh" + +#include + +#include +#include + +#include +#include + + +//////////////////////// +// the main function +//////////////////////// +void usage(const char *progname) +{ + std::cout << boost::format("usage: %s [--restart restartTime] tEnd\n")%progname; + exit(1); +} int main(int argc, char** argv) { - try{ - // define the problem dimensions - const int dim=2; /*@\label{tutorial-decoupled:dim}@*/ + try { + typedef TTAG(TutorialProblemDecoupled) TypeTag; + typedef GET_PROP_TYPE(TypeTag, PTAG(Scalar)) Scalar; + typedef GET_PROP_TYPE(TypeTag, PTAG(Grid)) Grid; + typedef GET_PROP_TYPE(TypeTag, PTAG(Problem)) Problem; + typedef Dune::FieldVector GlobalPosition; - // create a grid object - typedef double Scalar; /*@\label{tutorial-decoupled:grid-begin}@*/ - typedef Dune::SGrid Grid; - typedef Grid::LevelGridView GridView; - typedef Dune::FieldVector FieldVector; - Dune::FieldVector N(10); N[0] = 30; - FieldVector L(0); - FieldVector H(60); H[0] = 300; + static const int dim = Grid::dimension; + + // initialize MPI, finalize is done automatically on exit + Dune::MPIHelper::instance(argc, argv); + + //////////////////////////////////////////////////////////// + // parse the command line arguments + //////////////////////////////////////////////////////////// + if (argc < 2) + usage(argv[0]); + + // deal with the restart stuff + int argPos = 1; + bool restart = false; + double restartTime = 0; + if (std::string("--restart") == argv[argPos]) { + restart = true; + ++argPos; + + std::istringstream(argv[argPos++]) >> restartTime; + } + + if (argc - argPos != 1) { + usage(argv[0]); + } + + // read the initial time step and the end time + double tEnd, dt; + std::istringstream(argv[argPos++]) >> tEnd; + dt = tEnd; + + //////////////////////////////////////////////////////////// + // create the grid + //////////////////////////////////////////////////////////// + Dune::FieldVector N(1); N[0] = 100; + Dune::FieldVector L(0); + Dune::FieldVector H(60); H[0] = 300; Grid grid(N,L,H); - GridView gridView(grid.levelView(0));/*@\label{tutorial-decoupled:grid-end}@*/ + //////////////////////////////////////////////////////////// + // instantiate and run the concrete problem + //////////////////////////////////////////////////////////// - // define fluid and solid properties and constitutive relationships - Dumux::Water wettingfluid; /*@\label{tutorial-decoupled:water}@*/ - Dumux::LowViscosityOil nonwettingfluid; /*@\label{tutorial-decoupled:oil}@*/ - Dumux::TutorialSoil soil; /*@\label{tutorial-decoupled:soil}@*/ - Dumux::TwoPhaseRelations materialLaw(soil, wettingfluid, nonwettingfluid);/*@\label{tutorial-decoupled:twophaserelations}@*/ + Problem problem(grid.leafView(), L, H); - // create object containing the variables - typedef Dumux::VariableClass VariableClass; - VariableClass variables(gridView); + // load restart file if necessarry + if (restart) + problem.deserialize(restartTime); - //choose kind of two-phase model. Default: pw, Sw, vtotal - struct Dumux::DefineModel modelDef; -// modelDef.pressureType = modelDef.pressureW; -// modelDef.saturationType = modelDef.saturationW; -// modelDef.velocityType = modelDef.velocityTotal; - - // create object including the problem definition - typedef Dumux::TutorialProblemDecoupled Problem; - Problem problem(variables, wettingfluid, nonwettingfluid, soil, materialLaw,L, H); /*@\label{tutorial-decoupled:problem}@*/ - - // create object including the discretisation of the pressure equation - typedef Dumux::FVVelocity2P Diffusion; - Diffusion diffusion(gridView, problem, modelDef); /*@\label{tutorial-decoupled:diffusion}@*/ - - // create object including the space discretisation of the saturation equation - typedef Dumux::FVSaturation2P Transport; - Transport transport(gridView, problem, modelDef); /*@\label{tutorial-decoupled:transport}@*/ - - // some parameters used in the IMPES-object - int iterFlag = 0; - int nIter = 2; - double maxDefect = 1e-5; - - // create object including the IMPES (IMplicit Pressure Explicit Saturation) algorithm - typedef Dune::IMPES IMPES; - IMPES impes(diffusion, transport, iterFlag, nIter, maxDefect); /*@\label{tutorial-decoupled:impes}@*/ - - // some parameters needed for the TimeLoop-object - double tStart = 0; // start simulation at t = tStart - double tEnd = 4e7; // stop simulation at t = tEnd - const char* fileName = "tutorial_decoupled"; // name of the output files - int modulo = 1; // define time step interval in which output files are generated - double cFLFactor = 0.9; // security factor for the Courant-Friedrichs-Lewy-Criterion - - // create TimeLoop-object - Dumux::TimeLoop timeloop(gridView, tStart, tEnd, fileName, modulo, cFLFactor); /*@\label{tutorial-decoupled:timeloop}@*/ - - Dune::Timer timer; - timer.reset(); - - // start simulation - timeloop.execute(impes); /*@\label{tutorial-decoupled:execute}@*/ + // run the simulation + if (!problem.simulate(dt, tEnd)) + return 2; return 0; } - catch (Dune::Exception &e){ + catch (Dune::Exception &e) { std::cerr << "Dune reported error: " << e << std::endl; - return 1; } - catch (...){ - std::cerr << "Unknown exception thrown!" << std::endl; - return 1; + catch (...) { + std::cerr << "Unknown exception thrown!\n"; + throw; } + + return 3; } diff --git a/examples/tutorialproblem_coupled.hh b/examples/tutorialproblem_coupled.hh index 395e1e539..ed73369d5 100644 --- a/examples/tutorialproblem_coupled.hh +++ b/examples/tutorialproblem_coupled.hh @@ -77,7 +77,7 @@ SET_PROP(TutorialProblemCoupled, FluidSystem) /*@\label{tutorial-coupled:set-f // Set the spatial parameters SET_PROP(TutorialProblemCoupled, SpatialParameters) /*@\label{tutorial-coupled:set-spatialparameters}@*/ { - typedef Dumux::TutorialSpatialParameters type; + typedef Dumux::TutorialSpatialParametersCoupled type; }; // Disable gravity diff --git a/examples/tutorialproblem_decoupled.hh b/examples/tutorialproblem_decoupled.hh index ac559a901..6fa32ffdb 100644 --- a/examples/tutorialproblem_decoupled.hh +++ b/examples/tutorialproblem_decoupled.hh @@ -1,130 +1,260 @@ // $Id$ /***************************************************************************** - * Copyright (C) 2008-2009 by Markus Wolff * - * Institute of Hydraulic Engineering * - * University of Stuttgart, Germany * - * email: .@iws.uni-stuttgart.de * - * * - * 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, as long as this copyright notice * - * is included in its original form. * - * * - * This program is distributed WITHOUT ANY WARRANTY. * - *****************************************************************************/ -#ifndef TUTORIALPROBLEM_DECOUPLED_HH -#define TUTORIALPROBLEM_DECOUPLED_HH +* Copyright (C) 2007-2008 by Klaus Mosthaf * +* Copyright (C) 2007-2008 by Bernd Flemisch * +* Copyright (C) 2008-2009 by Andreas Lauser * +* Institute of Hydraulic Engineering * +* University of Stuttgart, Germany * +* email: .@iws.uni-stuttgart.de * +* * +* 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, as long as this copyright notice * +* is included in its original form. * +* * +* This program is distributed WITHOUT ANY WARRANTY. * +*****************************************************************************/ +#ifndef DUMUX_TUTORIALPROBLEM_DECOUPLED_HH +#define DUMUX_TUTORIALPROBLEM_DECOUPLED_HH -#include "dumux/fractionalflow/fractionalflowproblem.hh" +#if HAVE_UG +#include +#endif + +#include +#include + +#include +#include + +#include +#include +#include +#include +#include + +#include "tutorialspatialparameters_decoupled.hh" namespace Dumux { -/** \todo Please doc me! */ +template +class TutorialProblemDecoupled; -template class TutorialProblemDecoupled /*@\label{tutorial-decoupled:tutorialproblem}@*/ - : public FractionalFlowProblem +////////// +// Specify the properties for the lens problem +////////// +namespace Properties { - enum - {dim=GridView::dimension, dimWorld = GridView::dimensionworld}; - enum{wetting = 0, nonwetting = 1}; - typedef typename GridView::Grid Grid; - typedef typename GridView::Intersection Intersection; - typedef typename GridView::Traits::template Codim<0>::Entity Element; - typedef Dune::FieldVector LocalPosition; - typedef Dune::FieldVector GlobalPosition; +NEW_TYPE_TAG(TutorialProblemDecoupled, INHERITS_FROM(DecoupledTwoP, Transport)); + +// Set the grid type +SET_PROP(TutorialProblemDecoupled, Grid) +{ + // typedef Dune::YaspGrid<2> type; + typedef Dune::SGrid<2, 2> type; +}; + +// Set the problem property +SET_PROP(TutorialProblemDecoupled, Problem) +{ +public: + typedef Dumux::TutorialProblemDecoupled type; +}; + +// Set the model properties +SET_PROP(TutorialProblemDecoupled, SaturationModel) +{ + typedef Dumux::FVSaturation2P type; +}; + +SET_PROP(TutorialProblemDecoupled, PressureModel) +{ + typedef Dumux::FVVelocity2P type; +}; + +SET_INT_PROP(TutorialProblemDecoupled, VelocityFormulation, + GET_PROP_TYPE(TypeTag, PTAG(TwoPIndices))::velocityW); + +//SET_INT_PROP(TutorialProblemDecoupled, PressureFormulation, +// GET_PROP_TYPE(TypeTag, PTAG(TwoPIndices))::pressureGlobal); + +// Set the wetting phase +SET_PROP(TutorialProblemDecoupled, WettingPhase) +{ +private: + typedef typename GET_PROP_TYPE(TypeTag, PTAG(Scalar)) Scalar; +public: + typedef Dumux::LiquidPhase > type; +}; + +// Set the non-wetting phase +SET_PROP(TutorialProblemDecoupled, NonwettingPhase) +{ +private: + typedef typename GET_PROP_TYPE(TypeTag, PTAG(Scalar)) Scalar; +public: + typedef Dumux::LiquidPhase > type; +}; + +// Set the soil properties +SET_PROP(TutorialProblemDecoupled, SpatialParameters) +{ +private: + typedef typename GET_PROP_TYPE(TypeTag, PTAG(Grid)) Grid; + typedef typename GET_PROP_TYPE(TypeTag, PTAG(Scalar)) Scalar; public: - TutorialProblemDecoupled(VariableClass& variables, Fluid& wettingphase, Fluid& nonwettingphase, Matrix2p& soil, - TwoPhaseRelations& materialLaw = *(new TwoPhaseRelations), - const Dune::FieldVector Left = 0, const Dune::FieldVector Right = 0) - : FractionalFlowProblem(variables, wettingphase, nonwettingphase, soil, materialLaw), - Left_(Left[0]), Right_(Right[0]), eps_(1e-8) - {} + typedef Dumux::TutorialSpatialParametersDecoupled type; +}; - // function returning source/sink terms for the pressure equation - // depending on the position within the domain - virtual std::vector source(const GlobalPosition& globalPos, - const Element& e, /*@\label{tutorial-decoupled:qpress}@*/ - const LocalPosition& localPos) +SET_TYPE_PROP(TutorialProblemDecoupled, DiffusivePart, Dumux::CapillaryDiffusion); + +// Disable gravity +SET_BOOL_PROP(TutorialProblemDecoupled, EnableGravity, false); + +SET_SCALAR_PROP(TutorialProblemDecoupled, CFLFactor, 0.3); +} + +/*! +* \ingroup DecoupledProblems +*/ +template +class TutorialProblemDecoupled: public IMPESProblem2P > +{ + typedef TutorialProblemDecoupled ThisType; + typedef IMPESProblem2P ParentType; + typedef typename GET_PROP_TYPE(TypeTag, PTAG(GridView)) GridView; + + typedef typename GET_PROP_TYPE(TypeTag, PTAG(TwoPIndices)) Indices; + + typedef typename GET_PROP_TYPE(TypeTag, PTAG(FluidSystem)) FluidSystem; + typedef typename GET_PROP_TYPE(TypeTag, PTAG(FluidState)) FluidState; + + enum + { + dim = GridView::dimension, dimWorld = GridView::dimensionworld + }; + + enum + { + wPhaseIdx = Indices::wPhaseIdx, nPhaseIdx = Indices::nPhaseIdx + }; + + typedef typename GET_PROP_TYPE(TypeTag, PTAG(Scalar)) Scalar; + + typedef typename GridView::Traits::template Codim<0>::Entity Element; + typedef typename GridView::Intersection Intersection; + typedef Dune::FieldVector GlobalPosition; + typedef Dune::FieldVector LocalPosition; + +public: + TutorialProblemDecoupled(const GridView &gridView, const GlobalPosition lowerLeft = 0, const GlobalPosition upperRight = 0) : + ParentType(gridView), lowerLeft_(lowerLeft), upperRight_(upperRight) { - return std::vector(2,0.0); } - using FractionalFlowProblem::bctypePress; + /*! + * \name Problem parameters + */ + // \{ - // function returning the boundary condition type for solution - // of the pressure equation depending on the position within the domain - typename BoundaryConditions::Flags bctypePress(const GlobalPosition& globalPos, const Intersection& intersection) const /*@\label{tutorial-decoupled:bctypepress}@*/ + /*! + * \brief The problem name. + * + * This is used as a prefix for files generated by the simulation. + */ + const char *name() const { - if (globalPos[0] < eps_) + return "tutorial_decoupled"; + } + + bool shouldWriteRestartFile() const + { + return false; + } + + /*! + * \brief Returns the temperature within the domain. + * + * This problem assumes a temperature of 10 degrees Celsius. + */ + Scalar temperature(const GlobalPosition& globalPos, const Element& element) const + { + return 273.15 + 10; // -> 10°C + } + + // \} + + Scalar referencePressure(const GlobalPosition& globalPos, const Element& element) const + { + return 1e5; // -> 10°C + } + + std::vector source(const GlobalPosition& globalPos, const Element& element) { - return BoundaryConditions::dirichlet; + return std::vector(2, 0.0); } + + typename BoundaryConditions::Flags bctypePress(const GlobalPosition& globalPos, const Intersection& intersection) const + { + if ((globalPos[0] < lowerLeft_[0] + eps_)) + return BoundaryConditions::dirichlet; // all other boundaries return BoundaryConditions::neumann; } - // function returning the boundary condition type for solution - // of the saturation equation depending on the position within the domain - BoundaryConditions::Flags bctypeSat (const GlobalPosition& globalPos, const Intersection& intersection) const /*@\label{tutorial-decoupled:bctypesat}@*/ + BoundaryConditions::Flags bctypeSat(const GlobalPosition& globalPos, const Intersection& intersection) const { - if (globalPos[0]> (Right_ - eps_) || globalPos[0] < eps_) - { + if (globalPos[0] < lowerLeft_[0] + eps_) return Dumux::BoundaryConditions::dirichlet; - } + else + return Dumux::BoundaryConditions::neumann; + } + + Scalar dirichletPress(const GlobalPosition& globalPos, const Intersection& intersection) const + { + if (globalPos[0] < lowerLeft_[0] + eps_) + return 2e5; // all other boundaries - return Dumux::BoundaryConditions::neumann; + return 0; } - // function returning the Dirichlet boundary condition for the solution - // of the pressure equation depending on the position within the domain - Scalar dirichletPress(const GlobalPosition& globalPos, const Intersection& intersection) const /*@\label{tutorial-decoupled:gpress}@*/ + Scalar dirichletSat(const GlobalPosition& globalPos, const Intersection& intersection) const { - return 1e6; + if (globalPos[0] < lowerLeft_[0] + eps_) + return 1; + // all other boundaries + return 0; } - // function returning the Dirichlet boundary condition for the solution - // of the saturation equation depending on the position within the domain - Scalar dirichletSat(const GlobalPosition& globalPos, const Intersection& intersection) const /*@\label{tutorial-decoupled:gsat}@*/ + std::vector neumannPress(const GlobalPosition& globalPos, const Intersection& intersection) const { - if (globalPos[0] < eps_) + std::vector neumannFlux(2,0.0); + if (globalPos[0] > upperRight_[0] - eps_) { - return 1.0; + neumannFlux[nPhaseIdx] = 3e-4; } - // all other boundaries - return 0.0; - } - - using FractionalFlowProblem::neumannPress; - - // function returning the Neumann boundary condition for the solution - // of the pressure equation depending on the position within the domain - std::vector neumannPress(const GlobalPosition& globalPos, const Intersection& intersection) const /*@\label{tutorial-decoupled:jpress}@*/ - { - std::vector neumannFlux(2, 0.0); - if (globalPos[0]> Right_ - eps_) - { - neumannFlux[nonwetting] = 3e-4; - } - // all other boundaries return neumannFlux; } - // function returning the initial saturation - // depending on the position within the domain - Scalar initSat (const GlobalPosition& globalPos, const Element& e, /*@\label{tutorial-decoupled:initsat}@*/ - const Dune::FieldVector& xi) const + Scalar neumannSat(const GlobalPosition& globalPos, const Intersection& intersection, Scalar factor) const { - return 0.0; + return 0; + } + + Scalar initSat(const GlobalPosition& globalPos, const Element& element) const + { + return 0; } private: - Scalar Left_; - Scalar Right_; + GlobalPosition lowerLeft_; + GlobalPosition upperRight_; - Scalar eps_; + static const Scalar eps_ = 1e-6; }; -} // end namespace +} //end namespace + #endif diff --git a/examples/tutorialspatialparameters_coupled.hh b/examples/tutorialspatialparameters_coupled.hh index 93e31edbb..898d5719c 100644 --- a/examples/tutorialspatialparameters_coupled.hh +++ b/examples/tutorialspatialparameters_coupled.hh @@ -27,7 +27,7 @@ namespace Dumux { template -class TutorialSpatialParameters: public BoxSpatialParameters /*@\label{tutorial-coupled:tutorialSpatialParameters}@*/ +class TutorialSpatialParametersCoupled: public BoxSpatialParameters /*@\label{tutorial-coupled:tutorialSpatialParameters}@*/ { // Get informations for current implementation via property system typedef typename GET_PROP_TYPE(TypeTag, PTAG(Grid)) Grid; @@ -81,7 +81,7 @@ public: } // constructor - TutorialSpatialParameters(const GridView& gridView) : + TutorialSpatialParametersCoupled(const GridView& gridView) : BoxSpatialParameters(gridView), K_(0) { for (int i = 0; i < dim; i++) diff --git a/examples/tutorialspatialparameters_decoupled.hh b/examples/tutorialspatialparameters_decoupled.hh new file mode 100644 index 000000000..5a73e4015 --- /dev/null +++ b/examples/tutorialspatialparameters_decoupled.hh @@ -0,0 +1,98 @@ +// $Id: test_2p_spatialparamsinjection.hh 3456 2010-04-09 12:11:51Z mwolff $ +/***************************************************************************** + * Copyright (C) 2008-2009 by Markus Wolff * + * Institute of Hydraulic Engineering * + * University of Stuttgart, Germany * + * email: .@iws.uni-stuttgart.de * + * * + * 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, as long as this copyright notice * + * is included in its original form. * + * * + * This program is distributed WITHOUT ANY WARRANTY. * + *****************************************************************************/ +#ifndef TUTORIALSPATIALPARAMETERS_DECOUPLED_HH +#define TUTORIALSPATIALPARAMETERS_DECOUPLED_HH + + +//#include +#include +#include + +namespace Dumux +{ + +/** \todo Please doc me! */ + +template +class TutorialSpatialParametersDecoupled +{ + typedef typename GET_PROP_TYPE(TypeTag, PTAG(Grid)) Grid; + typedef typename GET_PROP_TYPE(TypeTag, PTAG(GridView)) GridView; + typedef typename GET_PROP_TYPE(TypeTag, PTAG(Scalar)) Scalar; + typedef typename Grid::ctype CoordScalar; + + enum + {dim=Grid::dimension, dimWorld=Grid::dimensionworld, numEq=1}; + typedef typename Grid::Traits::template Codim<0>::Entity Element; + + typedef Dune::FieldVector GlobalPosition; + typedef Dune::FieldVector LocalPosition; + typedef Dune::FieldMatrix FieldMatrix; + + typedef RegularizedBrooksCorey RawMaterialLaw; +// typedef LinearMaterial RawMaterialLaw; +public: + typedef EffToAbsLaw MaterialLaw; + typedef typename MaterialLaw::Params MaterialLawParams; + + void update (Scalar saturationW, const Element& element) + { + + } + + const FieldMatrix& intrinsicPermeability (const GlobalPosition& globalPos, const Element& element) const + { + return K_; + } + + double porosity(const GlobalPosition& globalPos, const Element& element) const + { + return 0.2; + } + + + // return the brooks-corey context depending on the position + const MaterialLawParams& materialLawParams(const GlobalPosition& globalPos, const Element &element) const + { + return materialLawParams_; + } + + + TutorialSpatialParametersDecoupled(const GridView& gridView) + : K_(0) + { + for (int i = 0; i < dim; i++) + K_[i][i] = 1e-7; + + // residual saturations + materialLawParams_.setSwr(0); + materialLawParams_.setSnr(0); + + // parameters for the Brooks-Corey Law + // entry pressures + materialLawParams_.setPe(10000); + + // Brooks-Corey shape parameters + materialLawParams_.setAlpha(2); + } + +private: + MaterialLawParams materialLawParams_; + FieldMatrix K_; +}; + +} // end namespace +#endif