From 764429db0c444553b10ac143211020ed05dc77ca Mon Sep 17 00:00:00 2001 From: Bernd Flemisch <bernd@iws.uni-stuttgart.de> Date: Thu, 6 Dec 2012 16:44:59 +0000 Subject: [PATCH] implicit branch: enhance the mpnc model with cell-centered discretization. Unify the obstacle test. git-svn-id: svn://svn.iws.uni-stuttgart.de/DUMUX/dumux/branches/implicit@9786 2fb0f335-1f38-0410-981e-8018bf24f1b0 --- .../mpnc/energy/mpncvtkwriterenergy.hh | 21 ++--- dumux/implicit/mpnc/mass/mpncvtkwritermass.hh | 7 +- dumux/implicit/mpnc/mpncvtkwritercommon.hh | 80 +++++++++++-------- test/implicit/mpnc/Makefile.am | 5 +- test/implicit/mpnc/obstacleproblem.hh | 57 +++++++------ test/implicit/mpnc/test_boxmpnc.cc | 59 ++++++++++++++ test/implicit/mpnc/test_boxmpnc.input | 34 ++++++++ .../mpnc/{test_mpnc.cc => test_ccmpnc.cc} | 2 +- .../{test_mpnc.input => test_ccmpnc.input} | 3 + 9 files changed, 193 insertions(+), 75 deletions(-) create mode 100644 test/implicit/mpnc/test_boxmpnc.cc create mode 100644 test/implicit/mpnc/test_boxmpnc.input rename test/implicit/mpnc/{test_mpnc.cc => test_ccmpnc.cc} (98%) rename test/implicit/mpnc/{test_mpnc.input => test_ccmpnc.input} (97%) diff --git a/dumux/implicit/mpnc/energy/mpncvtkwriterenergy.hh b/dumux/implicit/mpnc/energy/mpncvtkwriterenergy.hh index 558d881b8d..79be576f16 100644 --- a/dumux/implicit/mpnc/energy/mpncvtkwriterenergy.hh +++ b/dumux/implicit/mpnc/energy/mpncvtkwriterenergy.hh @@ -61,7 +61,8 @@ class MPNCVtkWriterEnergy : public MPNCVtkWriterModule<TypeTag> typedef typename ParentType::ScalarVector ScalarVector; typedef typename ParentType::PhaseVector PhaseVector; - + enum { isBox = GET_PROP_VALUE(TypeTag, ImplicitIsBox) }; + public: MPNCVtkWriterEnergy(const Problem &problem) : ParentType(problem) @@ -76,7 +77,7 @@ public: template <class MultiWriter> void allocBuffers(MultiWriter &writer) { - if (temperatureOutput_) this->resizeScalarBuffer_(temperature_); + if (temperatureOutput_) this->resizeScalarBuffer_(temperature_, isBox); } /*! @@ -105,7 +106,7 @@ public: void commitBuffers(MultiWriter &writer) { if (temperatureOutput_) - this->commitScalarBuffer_(writer, "T", temperature_); + this->commitScalarBuffer_(writer, "T", temperature_, isBox); } private: @@ -142,7 +143,7 @@ class MPNCVtkWriterEnergy<TypeTag, /* enableEnergy = */ true, /* enableKineticEn typedef typename ParentType::ScalarVector ScalarVector; typedef typename ParentType::PhaseVector PhaseVector; - + enum { isBox = GET_PROP_VALUE(TypeTag, ImplicitIsBox) }; public: MPNCVtkWriterEnergy(const Problem &problem) @@ -160,9 +161,9 @@ public: template <class MultiWriter> void allocBuffers(MultiWriter &writer) { - if (temperatureOutput_) this->resizeScalarBuffer_(temperature_); - if (enthalpyOutput_) this->resizePhaseBuffer_(enthalpy_); - if (internalEnergyOutput_) this->resizePhaseBuffer_(internalEnergy_); + if (temperatureOutput_) this->resizeScalarBuffer_(temperature_, isBox); + if (enthalpyOutput_) this->resizePhaseBuffer_(enthalpy_, isBox); + if (internalEnergyOutput_) this->resizePhaseBuffer_(internalEnergy_, isBox); } /*! @@ -196,11 +197,11 @@ public: void commitBuffers(MultiWriter &writer) { if (temperatureOutput_) - this->commitScalarBuffer_(writer, "T", temperature_); + this->commitScalarBuffer_(writer, "T", temperature_, isBox); if (enthalpyOutput_) - this->commitPhaseBuffer_(writer, "h_%s", enthalpy_); + this->commitPhaseBuffer_(writer, "h_%s", enthalpy_, isBox); if (internalEnergyOutput_) - this->commitPhaseBuffer_(writer, "u_%s", internalEnergy_); + this->commitPhaseBuffer_(writer, "u_%s", internalEnergy_, isBox); } private: diff --git a/dumux/implicit/mpnc/mass/mpncvtkwritermass.hh b/dumux/implicit/mpnc/mass/mpncvtkwritermass.hh index 23583d2482..3968f741df 100644 --- a/dumux/implicit/mpnc/mass/mpncvtkwritermass.hh +++ b/dumux/implicit/mpnc/mass/mpncvtkwritermass.hh @@ -56,7 +56,8 @@ class MPNCVtkWriterMass : public MPNCVtkWriterModule<TypeTag> enum { dim = GridView::dimension }; enum { numComponents = GET_PROP_VALUE(TypeTag, NumComponents) }; - + enum { isBox = GET_PROP_VALUE(TypeTag, ImplicitIsBox) }; + typedef typename ParentType::ComponentVector ComponentVector; bool fugacityOutput_; @@ -74,7 +75,7 @@ public: template <class MultiWriter> void allocBuffers(MultiWriter &writer) { - if (fugacityOutput_) this->resizeComponentBuffer_(fugacity_); + if (fugacityOutput_) this->resizeComponentBuffer_(fugacity_, isBox); } /*! @@ -106,7 +107,7 @@ public: void commitBuffers(MultiWriter &writer) { if (fugacityOutput_) - this->commitComponentBuffer_(writer, "f_%s", fugacity_); + this->commitComponentBuffer_(writer, "f_%s", fugacity_, isBox); } private: diff --git a/dumux/implicit/mpnc/mpncvtkwritercommon.hh b/dumux/implicit/mpnc/mpncvtkwritercommon.hh index c92536bf6d..c6e6f74e10 100644 --- a/dumux/implicit/mpnc/mpncvtkwritercommon.hh +++ b/dumux/implicit/mpnc/mpncvtkwritercommon.hh @@ -62,7 +62,8 @@ class MPNCVtkWriterCommon : public MPNCVtkWriterModule<TypeTag> typedef Dune::FieldVector<Scalar, dim> DimVector; typedef Dune::BlockVector<DimVector> DimField; typedef Dune::array<DimField, numPhases> PhaseDimField; - + enum { isBox = GET_PROP_VALUE(TypeTag, ImplicitIsBox) }; + public: MPNCVtkWriterCommon(const Problem &problem) : ParentType(problem) @@ -78,6 +79,10 @@ public: massFracOutput_ = GET_PARAM_FROM_GROUP(TypeTag, bool, Vtk, AddMassFractions); moleFracOutput_ = GET_PARAM_FROM_GROUP(TypeTag, bool, Vtk, AddMoleFractions); molarityOutput_ = GET_PARAM_FROM_GROUP(TypeTag, bool, Vtk, AddMolarities); + + // velocity output currently only works for box + if (!isBox) + velocityOutput_ = false; } /*! @@ -88,42 +93,47 @@ public: void allocBuffers(MultiWriter &writer) { if (porosityOutput_) - this->resizeScalarBuffer_(porosity_); + this->resizeScalarBuffer_(porosity_, isBox); if (boundaryTypesOutput_) - this->resizeScalarBuffer_(boundaryTypes_); + this->resizeScalarBuffer_(boundaryTypes_, isBox); if (velocityOutput_) { - Scalar nVerts = this->problem_.gridView().size(dim); + Scalar nDofs = this->problem_.model().numDofs(); for (int phaseIdx = 0; phaseIdx < numPhases; ++ phaseIdx) { - velocity_[phaseIdx].resize(nVerts); + velocity_[phaseIdx].resize(nDofs); velocity_[phaseIdx] = 0; } - this->resizeScalarBuffer_(boxSurface_); + this->resizeScalarBuffer_(boxSurface_, isBox); } - if (saturationOutput_) this->resizePhaseBuffer_(saturation_); - if (pressureOutput_) this->resizePhaseBuffer_(pressure_); - if (densityOutput_) this->resizePhaseBuffer_(density_); - if (mobilityOutput_) this->resizePhaseBuffer_(mobility_); - if (averageMolarMassOutput_) this->resizePhaseBuffer_(averageMolarMass_); - if (moleFracOutput_) this->resizePhaseComponentBuffer_(moleFrac_); - if (massFracOutput_) this->resizePhaseComponentBuffer_(massFrac_); - if (molarityOutput_) this->resizePhaseComponentBuffer_(molarity_); + if (saturationOutput_) this->resizePhaseBuffer_(saturation_, isBox); + if (pressureOutput_) this->resizePhaseBuffer_(pressure_, isBox); + if (densityOutput_) this->resizePhaseBuffer_(density_, isBox); + if (mobilityOutput_) this->resizePhaseBuffer_(mobility_, isBox); + if (averageMolarMassOutput_) this->resizePhaseBuffer_(averageMolarMass_, isBox); + if (moleFracOutput_) this->resizePhaseComponentBuffer_(moleFrac_, isBox); + if (massFracOutput_) this->resizePhaseComponentBuffer_(massFrac_, isBox); + if (molarityOutput_) this->resizePhaseComponentBuffer_(molarity_, isBox); } /*! * \brief Modify the internal buffers according to the volume * variables seen on an element */ - void processElement(const Element &elem, + void processElement(const Element &element, const FVElementGeometry &fvGeometry, const ElementVolumeVariables &elemVolVars, const ElementBoundaryTypes &elemBcTypes) { - int numLocalVertices = elem.geometry().corners(); - for (int localVertexIdx = 0; localVertexIdx < numLocalVertices; ++localVertexIdx) { - int globalIdx = this->problem_.vertexMapper().map(elem, localVertexIdx, dim); - const VolumeVariables &volVars = elemVolVars[localVertexIdx]; + for (int scvIdx = 0; scvIdx < fvGeometry.numSCV; ++scvIdx) + { + int globalIdx; + if (isBox) // vertex data + globalIdx = this->problem_.vertexMapper().map(element, scvIdx, dim); + else + globalIdx = this->problem_.elementMapper().map(element); + + const VolumeVariables &volVars = elemVolVars[scvIdx]; if (porosityOutput_) porosity_[globalIdx] = volVars.porosity(); @@ -132,7 +142,7 @@ public: // is used for a dirichlet condition int tmp = 0; for (int eqIdx = 0; eqIdx < numEq; ++eqIdx) { - if (elemBcTypes[localVertexIdx].isDirichlet(eqIdx)) + if (elemBcTypes[scvIdx].isDirichlet(eqIdx)) tmp += (1 << eqIdx); } if (boundaryTypesOutput_) boundaryTypes_[globalIdx] = tmp; @@ -155,13 +165,13 @@ public: if (velocityOutput_) { for (int faceIdx = 0; faceIdx < fvGeometry.numEdges; ++ faceIdx) { int i = fvGeometry.subContVolFace[faceIdx].i; - int I = this->problem_.vertexMapper().map(elem, i, dim); + int I = this->problem_.vertexMapper().map(element, i, dim); int j = fvGeometry.subContVolFace[faceIdx].j; - int J = this->problem_.vertexMapper().map(elem, j, dim); + int J = this->problem_.vertexMapper().map(element, j, dim); FluxVariables fluxVars(this->problem_, - elem, + element, fvGeometry, faceIdx, elemVolVars); @@ -191,41 +201,41 @@ public: void commitBuffers(MultiWriter &writer) { if (saturationOutput_) - this->commitPhaseBuffer_(writer, "S_%s", saturation_); + this->commitPhaseBuffer_(writer, "S_%s", saturation_, isBox); if (pressureOutput_) - this->commitPhaseBuffer_(writer, "p_%s", pressure_); + this->commitPhaseBuffer_(writer, "p_%s", pressure_, isBox); if (densityOutput_) - this->commitPhaseBuffer_(writer, "rho_%s", density_); + this->commitPhaseBuffer_(writer, "rho_%s", density_, isBox); if (averageMolarMassOutput_) - this->commitPhaseBuffer_(writer, "M_%s", averageMolarMass_); + this->commitPhaseBuffer_(writer, "M_%s", averageMolarMass_, isBox); if (mobilityOutput_) - this->commitPhaseBuffer_(writer, "lambda_%s", mobility_); + this->commitPhaseBuffer_(writer, "lambda_%s", mobility_, isBox); if (porosityOutput_) - this->commitScalarBuffer_(writer, "porosity", porosity_); + this->commitScalarBuffer_(writer, "porosity", porosity_, isBox); if (boundaryTypesOutput_) - this->commitScalarBuffer_(writer, "boundary types", boundaryTypes_); + this->commitScalarBuffer_(writer, "boundary types", boundaryTypes_, isBox); if (moleFracOutput_) - this->commitPhaseComponentBuffer_(writer, "x_%s^%s", moleFrac_); + this->commitPhaseComponentBuffer_(writer, "x_%s^%s", moleFrac_, isBox); if (massFracOutput_) - this->commitPhaseComponentBuffer_(writer, "X_%s^%s", massFrac_); + this->commitPhaseComponentBuffer_(writer, "X_%s^%s", massFrac_, isBox); if(molarityOutput_) - this->commitPhaseComponentBuffer_(writer, "c_%s^%s", molarity_); + this->commitPhaseComponentBuffer_(writer, "c_%s^%s", molarity_, isBox); if (velocityOutput_) { - int nVerts = this->problem_.gridView().size(dim); + int nDofs = this->problem_.model().numDofs(); for (int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) { // first, divide the velocity field by the // respective finite volume's surface area - for (int i = 0; i < nVerts; ++i) + for (int i = 0; i < nDofs; ++i) velocity_[phaseIdx][i] /= boxSurface_[i]; // commit the phase velocity std::ostringstream oss; diff --git a/test/implicit/mpnc/Makefile.am b/test/implicit/mpnc/Makefile.am index eacd0d0306..a48443cff5 100644 --- a/test/implicit/mpnc/Makefile.am +++ b/test/implicit/mpnc/Makefile.am @@ -1,11 +1,12 @@ -check_PROGRAMS = test_mpnc test_forchheimer2p test_forchheimer1p +check_PROGRAMS = test_boxmpnc test_ccmpnc test_forchheimer2p test_forchheimer1p noinst_HEADERS = *.hh EXTRA_DIST=*reference.vtu *reference.vtp *.input CMakeLists.txt grids/*.dgf gridsdir=$(datadir)/dumux/grids grids_DATA=grids/*.dgf -test_mpnc_SOURCES = test_mpnc.cc +test_boxmpnc_SOURCES = test_boxmpnc.cc +test_ccmpnc_SOURCES = test_ccmpnc.cc test_forchheimer2p_SOURCES = test_forchheimer2p.cc test_forchheimer1p_SOURCES = test_forchheimer1p.cc diff --git a/test/implicit/mpnc/obstacleproblem.hh b/test/implicit/mpnc/obstacleproblem.hh index ba8b58ee65..61a630e33e 100644 --- a/test/implicit/mpnc/obstacleproblem.hh +++ b/test/implicit/mpnc/obstacleproblem.hh @@ -48,8 +48,10 @@ class ObstacleProblem; namespace Properties { -NEW_TYPE_TAG(ObstacleProblem, INHERITS_FROM(BoxMPNC, ObstacleSpatialParams)); - +NEW_TYPE_TAG(ObstacleProblem, INHERITS_FROM(MPNC, ObstacleSpatialParams)); +NEW_TYPE_TAG(ObstacleBoxProblem, INHERITS_FROM(BoxModel, ObstacleProblem)); +NEW_TYPE_TAG(ObstacleCCProblem, INHERITS_FROM(CCModel, ObstacleProblem)); + // Set the grid type SET_TYPE_PROP(ObstacleProblem, Grid, Dune::YaspGrid<2>); @@ -158,7 +160,8 @@ class ObstacleProblem typedef Dune::FieldVector<typename GridView::Grid::ctype, dimWorld> GlobalPosition; typedef Dune::FieldVector<Scalar, numPhases> PhaseVector; typedef typename GET_PROP_TYPE(TypeTag, TimeManager) TimeManager; - + enum { isBox = GET_PROP_VALUE(TypeTag, ImplicitIsBox) }; + public: ObstacleProblem(TimeManager &timeManager, const GridView &gridView) : ParentType(timeManager, gridView) @@ -176,6 +179,10 @@ public: int np = 1000; FluidSystem::init(Tmin, Tmax, nT, pmin, pmax, np); + name_ = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, std::string, Problem, Name); + + const Element& element = *gridView.template begin<0>(); + elementWidth_ = element.geometry().corner(1)[0] - element.geometry().corner(0)[0]; } /*! @@ -221,8 +228,8 @@ public: * * This is used as a prefix for files generated by the simulation. */ - const char *name() const - { return "obstacle"; } + const std::string name() const + { return name_; } /*! * \brief Returns the temperature [K] within the domain. @@ -239,15 +246,14 @@ public: /*! * \brief Specifies which kind of boundary condition should be - * used for which equation on a given boundary segment. + * used for which equation on a given boundary control volume. * * \param values The boundary types for the conservation equations - * \param vertex The vertex for which the boundary type is set + * \param globalPos The position of the center of the finite volume */ - void boundaryTypes(BoundaryTypes &values, const Vertex &vertex) const + void boundaryTypesAtPos(BoundaryTypes &values, + const GlobalPosition &globalPos) const { - const GlobalPosition &globalPos = vertex.geometry().center(); - if (onInlet_(globalPos) || onOutlet_(globalPos)) values.setAllDirichlet(); else @@ -255,18 +261,17 @@ public: } /*! - * \brief Evaluate the boundary conditions for a Dirichlet - * boundary segment. + * \brief Evaluate the boundary conditions for a dirichlet + * control volume. * - * \param values The Dirichlet values for the primary variables - * \param vertex The vertex for which the boundary type is set + * \param values The dirichlet values for the primary variables + * \param globalPos The center of the finite volume which ought to be set. * * For this method, the \a values parameter stores primary variables. */ - void dirichlet(PrimaryVariables &values, const Vertex &vertex) const + void dirichletAtPos(PrimaryVariables &values, + const GlobalPosition &globalPos) const { - const GlobalPosition &globalPos = vertex.geometry().center(); - initial_(values, globalPos); } @@ -313,16 +318,15 @@ public: /*! * \brief Evaluate the initial value for a control volume. * + * \param values The initial values for the primary variables + * \param globalPos The center of the finite volume which ought to be set. + * * For this method, the \a values parameter stores primary * variables. */ - void initial(PrimaryVariables &values, - const Element &element, - const FVElementGeometry &fvGeometry, - const unsigned int scvIdx) const + void initialAtPos(PrimaryVariables &values, + const GlobalPosition &globalPos) const { - const GlobalPosition &globalPos = element.geometry().corner(scvIdx); - initial_(values, globalPos); Valgrind::CheckDefined(values); } @@ -350,7 +354,10 @@ private: // set the fluid temperatures fs.setTemperature(this->temperatureAtPos(globalPos)); - if (onInlet_(globalPos)) { + // the second condition is required for the cell-centered disc + if ((isBox && onInlet_(globalPos)) + || (globalPos[0] > 60 - elementWidth_/2.0 - eps_ && globalPos[1] <= 10)) + { // only liquid on inlet refPhaseIdx = wPhaseIdx; otherPhaseIdx = nPhaseIdx; @@ -436,6 +443,8 @@ private: Scalar temperature_; Scalar eps_; + std::string name_; + Scalar elementWidth_; }; } //end namespace diff --git a/test/implicit/mpnc/test_boxmpnc.cc b/test/implicit/mpnc/test_boxmpnc.cc new file mode 100644 index 0000000000..cfe511d11b --- /dev/null +++ b/test/implicit/mpnc/test_boxmpnc.cc @@ -0,0 +1,59 @@ +// -*- 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 Test for the MpNc box model. + */ +#include "config.h" +#include "obstacleproblem.hh" +#include <dumux/common/start.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) +{ + if (errorMsg.size() > 0) { + std::string errorMessageOut = "\nUsage: "; + errorMessageOut += progName; + errorMessageOut += " [options]\n"; + errorMessageOut += errorMsg; + errorMessageOut += "\n\nThe list of mandatory options for this program is:\n" + "\t-TimeManager.TEnd End of the simulation [s] \n" + "\t-TimeManager.DtInitial Initial timestep size [s] \n" + "\t-Grid.File Name of the file containing the grid \n" + "\t definition in DGF format\n"; + + std::cout << errorMessageOut + << "\n"; + } +} + +int main(int argc, char** argv) +{ + typedef TTAG(ObstacleBoxProblem) ProblemTypeTag; + return Dumux::start<ProblemTypeTag>(argc, argv, usage); +} diff --git a/test/implicit/mpnc/test_boxmpnc.input b/test/implicit/mpnc/test_boxmpnc.input new file mode 100644 index 0000000000..e973bddd89 --- /dev/null +++ b/test/implicit/mpnc/test_boxmpnc.input @@ -0,0 +1,34 @@ +############################################################### +# Parameter file for test_mpnc. +# Everything behind a '#' is a comment. +# Type "./test_mpnc --help" for more information. +############################################################### + +############################################################### +# Mandatory arguments +############################################################### + +[TimeManager] +DtInitial = 250 # [s] +TEnd = 1e4 # [s] + +[Grid] +File = ./grids/obstacle_24x16.dgf + +[LinearSolver] +ResidualReduction = 1e-12 + +[Problem] +Name = obstaclebox + +############################################################### +# Simulation restart +# +# DuMux simulations can be restarted from *.drs files +# Set Restart to the value of a specific file, +# e.g.: 'Restart = 27184.1' for the restart file +# name_time=27184.1_rank=0.drs +# Please comment in the two lines below, if restart is desired. +############################################################### +# [TimeManager] +# Restart = ... diff --git a/test/implicit/mpnc/test_mpnc.cc b/test/implicit/mpnc/test_ccmpnc.cc similarity index 98% rename from test/implicit/mpnc/test_mpnc.cc rename to test/implicit/mpnc/test_ccmpnc.cc index 20de9a1d6e..8c43e97d14 100644 --- a/test/implicit/mpnc/test_mpnc.cc +++ b/test/implicit/mpnc/test_ccmpnc.cc @@ -54,6 +54,6 @@ void usage(const char *progName, const std::string &errorMsg) int main(int argc, char** argv) { - typedef TTAG(ObstacleProblem) ProblemTypeTag; + typedef TTAG(ObstacleCCProblem) ProblemTypeTag; return Dumux::start<ProblemTypeTag>(argc, argv, usage); } diff --git a/test/implicit/mpnc/test_mpnc.input b/test/implicit/mpnc/test_ccmpnc.input similarity index 97% rename from test/implicit/mpnc/test_mpnc.input rename to test/implicit/mpnc/test_ccmpnc.input index a1bf8b9930..dc825e13a9 100644 --- a/test/implicit/mpnc/test_mpnc.input +++ b/test/implicit/mpnc/test_ccmpnc.input @@ -18,6 +18,9 @@ File = ./grids/obstacle_24x16.dgf [LinearSolver] ResidualReduction = 1e-12 +[Problem] +Name = obstaclecc + ############################################################### # Simulation restart # -- GitLab