I joint to you the file .input, he tutorialspatialparams_cdecoupled and the file .hh to show all parameters that i use. Thak's for the help. Best regards.
2015-10-21 21:10 GMT+02:00 Martin <[email protected]>: > That probably means that the matrix is singular and therefore the pressure > equation can not be solved. > What grid are you using and what intrinsic permeability values? Which > material law? > > > Regards, > Martin > > > On 10/21/2015 08:12 PM, Ait Mahiout Latifa wrote: > > Hi Martin, > ok, so i corrected the condition, so now, my code is: > > > void boundaryTypesAtPos(BoundaryTypes &bcTypes, const GlobalPosition& > globalPos) const /*@\label{tutorial-decoupled:bctype}@*/ > { > > > if ((globalPos[0] > 580- eps_) && (globalPos[1] > 580- > eps_) ) > { > bcTypes.setDirichlet(pressEqIdx); > bcTypes.setDirichlet(satEqIdx); > //bcTypes.setAllDirichlet(); // alternative if the same BC > is used for both types of equations > } > // all other boundaries > else if ((globalPos[0] < 20- eps_) && (globalPos[1] < 20- > eps_) ) > { > //bcTypes.setNeumann(pressEqIdx); > //bcTypes.setDirichlet(satEqIdx); > bcTypes.setAllNeumann(); // alternative if the same BC is > used for both types of equations > } > else > bcTypes.setAllNeumann(); > } > //! Value for dirichlet boundary condition at position globalPos. > /*! In case of a dirichlet BC for the pressure equation the pressure > \f$ [Pa] \f$, and for > * the transport equation the saturation [-] have to be defined on > boundaries. > * > * \param values Values of primary variables at the boundary > * \param intersection The boundary intersection > * > * Alternatively, the function dirichletAtPos(PrimaryVariables > &values, const GlobalPosition& globalPos) > * could be defined, where globalPos is the vector including the > global coordinates of the finite volume. > */ > void dirichletAtPos(PrimaryVariables &values, const GlobalPosition& > globalPos) const > { > values=0; > > if ((globalPos[0] > 580- eps_) && (globalPos[1] > 580- eps_) ) > { > values[pwIdx] = 3.5e7; > values[swIdx] = 0.0; > } > else if ((globalPos[0] < 20- eps_) && (globalPos[1] < 20- eps_) > ) > values[swIdx] = 1.0; > } > //! Value for neumann boundary condition \f$ [\frac{kg}{m^3 \cdot s}] > \f$ at position globalPos. > /*! In case of a neumann boundary condition, the flux of matter > * is returned as a vector. > * > * \param values Boundary flux values for the different phases > * \param globalPos The position of the center of the finite volume > * > * Alternatively, the function neumann(PrimaryVariables &values, > const Intersection& intersection) could be defined, > * where intersection is the boundary intersection. > */ > void neumannAtPos(PrimaryVariables &values, const GlobalPosition& > globalPos) const /*@\label{tutorial-decoupled:neumann}@*/ > { > values = 0; > if ((globalPos[0] < 20- eps_) && (globalPos[1] < 20- eps_) ) > { > values[nPhaseIdx] = -1e-8; > //values[wPhaseIdx] = 0.0; > } > else > > { > values[nPhaseIdx] = 0.0; > values[wPhaseIdx] = 0.0; > } > } > //! Initial condition at position globalPos. > /*! Only initial values for saturation have to be given! > * > * \param values Values of primary variables > * \param element The finite volume element > * > * Alternatively, the function initialAtPos(PrimaryVariables &values, > const GlobalPosition& globalPos) > * could be defined, where globalPos is the vector including the > global coordinates of the finite volume. > */ > void initial(PrimaryVariables &values, > const Element &element) const > /*@\label{tutorial-decoupled:initial}@*/ > { > values = 0; > } > > private: > const Scalar eps_; > }; > } //end namespace > > > > > with eps_=1e-6 > > > > and now i have this error > Don't panic... ! > > Rank 0: No parameter file given. Defaulting to > './tutorial_decoupled.input' for input file. > Initializing problem 'tutorial_decoupled' > Dune reported error: ISTLError > [apply:/home/latifa/Dumux_2.6.0/dune-istl-2.3.1/dune/istl/solvers.hh:651]: > breakdown in BiCGSTAB - rho 0 <= EPSILON 1e-80 after 8.5 iterations > > i don't understand why i have this error, and what does it mean? > > 2015-10-21 20:06 GMT+02:00 Martin Schneider < > [email protected]>: > >> Hi, >> >> the first if-query is wrong, corresponding to your boundary conditions: >> >> if x < 20m and y < 20 m: q_w.n = -1e-8 and Sw=1 >> if x > 580 and y > 580: pw= 150 bar and Sn=1 >> >> it should be >> if ((globalPos[0] > 580- eps_) && (globalPos[1] > 580- >> eps_) ) >> { >> bcTypes.setDirichlet(pressEqIdx); >> bcTypes.setDirichlet(satEqIdx); >> //bcTypes.setAllDirichlet(); // alternative if the same BC >> is used for both types of equations >> } >> >> instead of: >> if ((globalPos[0] > 600- eps_) && (globalPos[1] > 600- >> eps_) ) >> { >> bcTypes.setDirichlet(pressEqIdx); >> bcTypes.setDirichlet(satEqIdx); >> //bcTypes.setAllDirichlet(); // alternative if the same BC >> is used for both types of equations >> } >> >> Regards, >> Martin >> >> >> On 10/21/2015 07:34 PM, Ait Mahiout Latifa wrote: >> >> >> Hi, >> i want to change the boundary conditions in tutorial_decoupled problem, >> so that in an domain 600*600, are imposed the boundary conditions: >> if x < 20m and y < 20 m: q_w.n = -1e-8 and Sw=1 >> if x > 580 and y > 580: pw= 150 bar and Sn=1 >> and no flux in the other parts of the domain. >> So, i change x and y in the .input file, and i have the folowing >> modifications in the tutorial_decoupled.hh: >> void boundaryTypesAtPos(BoundaryTypes &bcTypes, const GlobalPosition& >> globalPos) const /*@\label{tutorial-decoupled:bctype}@*/ >> { >> if ((globalPos[0] > 600- eps_) && (globalPos[1] > 600- >> eps_) ) >> { >> bcTypes.setDirichlet(pressEqIdx); >> bcTypes.setDirichlet(satEqIdx); >> //bcTypes.setAllDirichlet(); // alternative if the same BC >> is used for both types of equations >> } >> // all other boundaries >> else if ((globalPos[0] < 20- eps_) && (globalPos[1] < >> 20- eps_) ) >> { >> bcTypes.setNeumann(pressEqIdx); >> bcTypes.setDirichlet(satEqIdx); >> //bcTypes.setAllNeumann(); // alternative if the same BC >> is used for both types of equations >> } >> else >> bcTypes.setAllNeumann(); >> } >> //! Value for dirichlet boundary condition at position globalPos. >> /*! In case of a dirichlet BC for the pressure equation the pressure >> \f$ [Pa] \f$, and for >> * the transport equation the saturation [-] have to be defined on >> boundaries. >> * >> * \param values Values of primary variables at the boundary >> * \param intersection The boundary intersection >> * >> * Alternatively, the function dirichletAtPos(PrimaryVariables >> &values, const GlobalPosition& globalPos) >> * could be defined, where globalPos is the vector including the >> global coordinates of the finite volume. >> */ >> void dirichletAtPos(PrimaryVariables &values, const GlobalPosition& >> globalPos) const >> { >> if ((globalPos[0] > 600- eps_) && (globalPos[1] > 600- eps_) ) >> { >> values[pwIdx] = 3.5e7; >> values[swIdx] = 0.0; >> } >> else if ((globalPos[0] < 20- eps_) && (globalPos[1] < 20- >> eps_) ) >> values[swIdx] = 1.0; >> } >> //! Value for neumann boundary condition \f$ [\frac{kg}{m^3 \cdot s}] >> \f$ at position globalPos. >> /*! In case of a neumann boundary condition, the flux of matter >> * is returned as a vector. >> * >> * \param values Boundary flux values for the different phases >> * \param globalPos The position of the center of the finite volume >> * >> * Alternatively, the function neumann(PrimaryVariables &values, >> const Intersection& intersection) could be defined, >> * where intersection is the boundary intersection. >> */ >> void neumannAtPos(PrimaryVariables &values, const GlobalPosition& >> globalPos) const /*@\label{tutorial-decoupled:neumann}@*/ >> { >> values = 0; >> if ((globalPos[0] < 20- eps_) && (globalPos[1] < 20- eps_) ) >> { >> values[nPhaseIdx] = -1e-8; >> values[wPhaseIdx] = 0.0; >> } >> else >> >> { >> values[nPhaseIdx] = 0.0; >> values[wPhaseIdx] = 0.0; >> } >> } >> //! Initial condition at position globalPos. >> /*! Only initial values for saturation have to be given! >> * >> * \param values Values of primary variables >> * \param element The finite volume element >> * >> * Alternatively, the function initialAtPos(PrimaryVariables >> &values, const GlobalPosition& globalPos) >> * could be defined, where globalPos is the vector including the >> global coordinates of the finite volume. >> */ >> void initial(PrimaryVariables &values, >> const Element &element) const >> /*@\label{tutorial-decoupled:initial}@*/ >> { >> values = 0; >> } >> >> >> There isn't a problem in compilation, but in execution, io have the >> folowing error: >> >> ./tutorial_decoupled >> >> Wherever he saw a hole he always wanted to know the depth of it. To him >> this was important. >> - Jules Verne, A journey to the center of the earth >> >> Rank 0: No parameter file given. Defaulting to >> './tutorial_decoupled.input' for input file. >> Initializing problem 'tutorial_decoupled' >> Dune reported error: ISTLError >> [apply:/home/latifa/Dumux_2.6.0/dune-istl-2.3.1/dune/istl/solvers.hh:679]: >> h=0 in BiCGSTAB >> >> So please, where os the problem in my definition of the boundary >> conditions? An how i can arrange it? >> Best regards. >> >> >> >> _______________________________________________ >> Dumux mailing >> [email protected]https://listserv.uni-stuttgart.de/mailman/listinfo/dumux >> >> >> >> -- >> M.Sc. Martin Schneider >> University of Stuttgart >> Institute for Modelling Hydraulic and Environmental Systems >> Department of Hydromechanics and Modelling of Hydrosystems >> Pfaffenwaldring 61 >> D-70569 Stuttgart >> Tel: (+49) 0711/ 685-69159 >> Fax: (+49) 0711/ 685-60430 >> E-Mail: [email protected] >> >> >> _______________________________________________ >> Dumux mailing list >> [email protected] >> https://listserv.uni-stuttgart.de/mailman/listinfo/dumux >> >> > > > _______________________________________________ > Dumux mailing > [email protected]https://listserv.uni-stuttgart.de/mailman/listinfo/dumux > > > > _______________________________________________ > Dumux mailing list > [email protected] > https://listserv.uni-stuttgart.de/mailman/listinfo/dumux > >
tutorial_coupled.input
Description: Binary data
// -*- 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 spatial parameters for the sequential tutorial */ #ifndef DUMUX_TUTORIAL_SPATIAL_PARAMS_DECOUPLED_HH #define DUMUX_TUTORIAL_SPATIAL_PARAMS_DECOUPLED_HH #include <dumux/material/spatialparams/fvspatialparams.hh> #include <dumux/material/fluidmatrixinteractions/2p/linearmaterial.hh> #include <dumux/material/fluidmatrixinteractions/2p/regularizedbrookscorey.hh> #include <dumux/material/fluidmatrixinteractions/2p/efftoabslaw.hh> namespace Dumux { //forward declaration template<class TypeTag> class TutorialSpatialParamsDecoupled; namespace Properties { // The spatial parameters TypeTag NEW_TYPE_TAG(TutorialSpatialParamsDecoupled); // Set the spatial parameters SET_TYPE_PROP(TutorialSpatialParamsDecoupled, SpatialParams, Dumux::TutorialSpatialParamsDecoupled<TypeTag>); /*@\label{tutorial-decoupled:set-spatialparameters}@*/ // Set the material law SET_PROP(TutorialSpatialParamsDecoupled, MaterialLaw) { private: // material law typedefs typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar; typedef RegularizedBrooksCorey<Scalar> RawMaterialLaw; public: typedef EffToAbsLaw<RawMaterialLaw> type; }; } //! Definition of the spatial parameters for the decoupled tutorial template<class TypeTag> class TutorialSpatialParamsDecoupled: public FVSpatialParams<TypeTag> { typedef FVSpatialParams<TypeTag> ParentType; typedef typename GET_PROP_TYPE(TypeTag, Grid) Grid; typedef typename GET_PROP_TYPE(TypeTag, GridView) GridView; typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar; typedef typename Grid::ctype CoordScalar; enum {dim=Grid::dimension, dimWorld=Grid::dimensionworld}; typedef typename Grid::Traits::template Codim<0>::Entity Element; typedef Dune::FieldVector<CoordScalar, dimWorld> GlobalPosition; typedef Dune::FieldMatrix<Scalar,dim,dim> FieldMatrix; public: typedef typename GET_PROP_TYPE(TypeTag, MaterialLaw) MaterialLaw; typedef typename MaterialLaw::Params MaterialLawParams; //! Intrinsic permeability tensor K \f$[m^2]\f$ depending /*! on the position in the domain * * \param element The finite volume element * * Alternatively, the function intrinsicPermeabilityAtPos(const GlobalPosition& globalPos) could be * defined, where globalPos is the vector including the global coordinates of the finite volume. */ const FieldMatrix& intrinsicPermeability (const Element& element) const { return K_; } //! Define the porosity \f$[-]\f$ of the porous medium depending /*! on the position in the domain * * \param element The finite volume element * * Alternatively, the function porosityAtPos(const GlobalPosition& globalPos) could be * defined, where globalPos is the vector including the global coordinates of the finite volume. */ double porosity(const Element& element) const { return 0.2; } /*! Return the parameter object for the material law (i.e. Brooks-Corey) * depending on the position in the domain * * \param element The finite volume element * * Alternatively, the function materialLawParamsAtPos(const GlobalPosition& globalPos) * could be defined, where globalPos is the vector including the global coordinates of * the finite volume. */ const MaterialLawParams& materialLawParams(const Element &element) const { return materialLawParams_; } //! Constructor TutorialSpatialParamsDecoupled(const GridView& gridView) : ParentType(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(500); // Brooks-Corey shape parameters materialLawParams_.setLambda(2); } private: MaterialLawParams materialLawParams_; FieldMatrix K_; }; } // end namespace #endif
// -*- 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 problem for the sequential tutorial */ #ifndef DUMUX_TUTORIALPROBLEM_DECOUPLED_HH // guardian macro /*@\label{tutorial-decoupled:guardian1}@*/ #define DUMUX_TUTORIALPROBLEM_DECOUPLED_HH // guardian macro /*@\label{tutorial-decoupled:guardian2}@*/ // the grid includes #include <dune/grid/yaspgrid.hh> #include <dumux/io/cubegridcreator.hh> // dumux 2p-decoupled environment #include <dumux/decoupled/2p/diffusion/fv/fvpressureproperties2p.hh> #include <dumux/decoupled/2p/transport/fv/fvtransportproperties2p.hh> #include <dumux/decoupled/2p/impes/impesproblem2p.hh> /*@\label{tutorial-decoupled:parent-problem}@*/ // assign parameters dependent on space (e.g. spatial parameters) #include "tutorialspatialparams_decoupled.hh" /*@\label{tutorial-decoupled:spatialparameters}@*/ // include cfl-criterion after coats: more suitable if the problem is not advection dominated #include<dumux/decoupled/2p/transport/fv/evalcflfluxcoats.hh> // the components that are used #include <dumux/material/components/h2o.hh> #include <dumux/material/components/lnapl.hh> namespace Dumux { template<class TypeTag> class TutorialProblemDecoupled; ////////// // Specify the properties for the lens problem ////////// namespace Properties { // create a new type tag for the problem NEW_TYPE_TAG(TutorialProblemDecoupled, INHERITS_FROM(FVPressureTwoP, FVTransportTwoP, IMPESTwoP, TutorialSpatialParamsDecoupled)); /*@\label{tutorial-decoupled:create-type-tag}@*/ // Set the problem property SET_PROP(TutorialProblemDecoupled, Problem) /*@\label{tutorial-decoupled:set-problem}@*/ { typedef Dumux::TutorialProblemDecoupled<TypeTag> type; }; // Set the grid type SET_TYPE_PROP(TutorialProblemDecoupled, Grid, Dune::YaspGrid<2>); /*@\label{tutorial-decoupled:set-grid-type}@*/ //Set the grid creator SET_TYPE_PROP(TutorialProblemDecoupled, GridCreator, Dumux::CubeGridCreator<TypeTag>); /*@\label{tutorial-decoupled:set-gridcreator}@*/ // Set the wetting phase SET_PROP(TutorialProblemDecoupled, WettingPhase) /*@\label{tutorial-decoupled:2p-system-start}@*/ { private: typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar; public: typedef Dumux::LiquidPhase<Scalar, Dumux::H2O<Scalar> > type; /*@\label{tutorial-decoupled:wettingPhase}@*/ }; // Set the non-wetting phase SET_PROP(TutorialProblemDecoupled, NonwettingPhase) { private: typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar; public: typedef Dumux::LiquidPhase<Scalar, Dumux::LNAPL<Scalar> > type; /*@\label{tutorial-decoupled:nonwettingPhase}@*/ }; /*@\label{tutorial-decoupled:2p-system-end}@*/ SET_TYPE_PROP(TutorialProblemDecoupled, EvalCflFluxFunction, Dumux::EvalCflFluxCoats<TypeTag>); /*@\label{tutorial-decoupled:cflflux}@*/ SET_SCALAR_PROP(TutorialProblemDecoupled, ImpetCFLFactor, 0.95); /*@\label{tutorial-decoupled:cflfactor}@*/ // Disable gravity SET_BOOL_PROP(TutorialProblemDecoupled, ProblemEnableGravity, false); /*@\label{tutorial-decoupled:gravity}@*/ } /*@\label{tutorial-decoupled:propertysystem-end}@*/ /*! \ingroup DecoupledProblems * @brief Problem class for the decoupled tutorial */ template<class TypeTag> class TutorialProblemDecoupled: public IMPESProblem2P<TypeTag> /*@\label{tutorial-decoupled:def-problem}@*/ { typedef IMPESProblem2P<TypeTag> ParentType; typedef typename GET_PROP_TYPE(TypeTag, GridView) GridView; typedef typename GET_PROP_TYPE(TypeTag, TimeManager) TimeManager; typedef typename GET_PROP_TYPE(TypeTag, Indices) Indices; typedef typename GET_PROP_TYPE(TypeTag, BoundaryTypes) BoundaryTypes; typedef typename GET_PROP(TypeTag, SolutionTypes) SolutionTypes; typedef typename SolutionTypes::PrimaryVariables PrimaryVariables; enum { dimWorld = GridView::dimensionworld }; enum { wPhaseIdx = Indices::wPhaseIdx, nPhaseIdx = Indices::nPhaseIdx, pwIdx = Indices::pwIdx, swIdx = Indices::swIdx, pressEqIdx = Indices::pressureEqIdx, satEqIdx = Indices::satEqIdx }; typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar; typedef typename GridView::Traits::template Codim<0>::Entity Element; typedef typename GridView::Intersection Intersection; typedef Dune::FieldVector<Scalar, dimWorld> GlobalPosition; public: TutorialProblemDecoupled(TimeManager &timeManager, const GridView &gridView) : ParentType(timeManager, gridView), eps_(1e-6)/*@\label{tutorial-decoupled:constructor-problem}@*/ { //write only every 10th time step to output file this->setOutputInterval(10);/*@\label{tutorial-decoupled:outputinterval}@*/ } //! The problem name. /*! This is used as a prefix for files generated by the simulation. */ const char *name() const /*@\label{tutorial-decoupled:name}@*/ { return "tutorial_decoupled"; } //! Returns true if a restart file should be written. /* The default behaviour is to write no restart file. */ bool shouldWriteRestartFile() const /*@\label{tutorial-decoupled:restart}@*/ { return false; } //! Returns the temperature within the domain at position globalPos. /*! This problem assumes a temperature of 10 degrees Celsius. * * \param element The finite volume element * * Alternatively, the function temperatureAtPos(const GlobalPosition& globalPos) could be * defined, where globalPos is the vector including the global coordinates of the finite volume. */ Scalar temperature(const Element& element) const /*@\label{tutorial-decoupled:temperature}@*/ { return 273.15 + 10; // -> 10°C } //! Returns a constant pressure to enter material laws at position globalPos. /* For incrompressible simulations, a constant pressure is necessary * to enter the material laws to gain a constant density etc. In the compressible * case, the pressure is used for the initialization of material laws. * * \param element The finite volume element * * Alternatively, the function referencePressureAtPos(const GlobalPosition& globalPos) could be * defined, where globalPos is the vector including the global coordinates of the finite volume. */ Scalar referencePressure(const Element& element) const /*@\label{tutorial-decoupled:refPressure}@*/ { return 2e5; } //! Source of mass \f$ [\frac{kg}{m^3 \cdot s}] \f$ of a finite volume. /*! Evaluate the source term for all phases within a given * volume. * * \param values Includes sources for the two phases * \param element The finite volume element * * The method returns the mass generated (positive) or * annihilated (negative) per volume unit. * * Alternatively, the function sourceAtPos(PrimaryVariables &values, const GlobalPosition& globalPos) * could be defined, where globalPos is the vector including the global coordinates of the finite volume. */ void source(PrimaryVariables &values, const Element& element) const /*@\label{tutorial-decoupled:source}@*/ { values = 0; } //! Type of boundary conditions at position globalPos. /*! Defines the type the boundary condition for the pressure equation, * either pressure (dirichlet) or flux (neumann), * and for the transport equation, * either saturation (dirichlet) or flux (neumann). * * \param bcTypes Includes the types of boundary conditions * \param globalPos The position of the center of the finite volume * * Alternatively, the function boundaryTypes(PrimaryVariables &values, const Intersection& * intersection) could be defined, where intersection is the boundary intersection. */ void boundaryTypesAtPos(BoundaryTypes &bcTypes, const GlobalPosition& globalPos) const /*@\label{tutorial-decoupled:bctype}@*/ { if ((globalPos[0] > 580- eps_) && (globalPos[1] > 580- eps_) ) { bcTypes.setDirichlet(pressEqIdx); bcTypes.setDirichlet(satEqIdx); //bcTypes.setAllDirichlet(); // alternative if the same BC is used for both types of equations } // all other boundaries //else if ((globalPos[0] < 20- eps_) && (globalPos[1] < 20- eps_) ) //{ //bcTypes.setNeumann(pressEqIdx); //bcTypes.setDirichlet(satEqIdx); //bcTypes.setAllNeumann(); // alternative if the same BC is used for both types of equations //} else bcTypes.setAllNeumann(); } //! Value for dirichlet boundary condition at position globalPos. /*! In case of a dirichlet BC for the pressure equation the pressure \f$ [Pa] \f$, and for * the transport equation the saturation [-] have to be defined on boundaries. * * \param values Values of primary variables at the boundary * \param intersection The boundary intersection * * Alternatively, the function dirichletAtPos(PrimaryVariables &values, const GlobalPosition& globalPos) * could be defined, where globalPos is the vector including the global coordinates of the finite volume. */ void dirichletAtPos(PrimaryVariables &values, const GlobalPosition& globalPos) const { values=0; if ((globalPos[0] > 580- eps_) && (globalPos[1] > 580- eps_) ) { values[pwIdx] = 1.5e7; values[swIdx] = 0.0; } else if ((globalPos[0] < 20- eps_) && (globalPos[1] < 20- eps_) ) values[swIdx] = 1.0; } //! Value for neumann boundary condition \f$ [\frac{kg}{m^3 \cdot s}] \f$ at position globalPos. /*! In case of a neumann boundary condition, the flux of matter * is returned as a vector. * * \param values Boundary flux values for the different phases * \param globalPos The position of the center of the finite volume * * Alternatively, the function neumann(PrimaryVariables &values, const Intersection& intersection) could be defined, * where intersection is the boundary intersection. */ void neumannAtPos(PrimaryVariables &values, const GlobalPosition& globalPos) const /*@\label{tutorial-decoupled:neumann}@*/ { values = 0; if ((globalPos[0] < 20- eps_) && (globalPos[1] < 20- eps_) ) { values[nPhaseIdx] = -1e-8; //values[wPhaseIdx] = 0.0; } else { values[nPhaseIdx] = 0.0; values[wPhaseIdx] = 0.0; } } //! Initial condition at position globalPos. /*! Only initial values for saturation have to be given! * * \param values Values of primary variables * \param element The finite volume element * * Alternatively, the function initialAtPos(PrimaryVariables &values, const GlobalPosition& globalPos) * could be defined, where globalPos is the vector including the global coordinates of the finite volume. */ void initial(PrimaryVariables &values, const Element &element) const /*@\label{tutorial-decoupled:initial}@*/ { values = 0; } private: const Scalar eps_; }; } //end namespace #endif
_______________________________________________ Dumux mailing list [email protected] https://listserv.uni-stuttgart.de/mailman/listinfo/dumux
