From b56e9025d99b3cac32b352525a3c740510f8591c Mon Sep 17 00:00:00 2001
From: Holger Class <holger.class@iws.uni-stuttgart.de>
Date: Fri, 8 Jul 2011 12:49:49 +0000
Subject: [PATCH] git-svn-id:
svn://svn.iws.uni-stuttgart.de/DUMUX/dumux/trunk@6164
2fb0f335-1f38-0410-981e-8018bf24f1b0
---
test/boxmodels/2p2cni/infiltrationproblem.hh | 405 +++++++++++++++++++
1 file changed, 405 insertions(+)
create mode 100644 test/boxmodels/2p2cni/infiltrationproblem.hh
diff --git a/test/boxmodels/2p2cni/infiltrationproblem.hh b/test/boxmodels/2p2cni/infiltrationproblem.hh
new file mode 100644
index 0000000000..53a8afad7f
--- /dev/null
+++ b/test/boxmodels/2p2cni/infiltrationproblem.hh
@@ -0,0 +1,405 @@
+/*****************************************************************************
+ * Copyright (C) 2011 by Holger Class *
+ * Institute of Hydraulic Engineering *
+ * University of Stuttgart, Germany *
+ * email: <givenname>.<name>@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. *
+ * *
+ * 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 Non-isothermal gas injection problem where a gas (e.g. air)
+ * is injected into a fully water saturated medium.
+ */
+#ifndef DUMUX_INFILTRATIONPROBLEM_HH
+#define DUMUX_INFILTRATIONPROBLEM_HH
+
+#include <dune/grid/io/file/dgfparser/dgfug.hh>
+#include <dune/grid/io/file/dgfparser/dgfs.hh>
+#include <dune/grid/io/file/dgfparser/dgfyasp.hh>
+
+#include <dumux/material/fluidsystems/holle_h20_air_mesitylene_system.hh>
+
+#include <dumux/boxmodels/holle3p3c/holle3p3cmodel.hh>
+
+#include "infiltrationspatialparameters.hh"
+
+#define ISOTHERMAL 0
+
+namespace Dumux
+{
+template <class TypeTag>
+class InfiltrationProblem;
+
+namespace Properties
+{
+NEW_TYPE_TAG(InfiltrationProblem, INHERITS_FROM(BoxHolleThreePThreeC));
+
+// Set the grid type
+SET_PROP(InfiltrationProblem, Grid)
+{
+ typedef Dune::YaspGrid<2> type;
+};
+
+#if HAVE_DUNE_PDELAB
+SET_PROP(InfiltrationProblem, LocalFEMSpace)
+{
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(Scalar)) Scalar;
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(GridView)) GridView;
+ enum{dim = GridView::dimension};
+
+public:
+ typedef Dune::PDELab::Q1LocalFiniteElementMap<Scalar,Scalar,dim> type; // for cubes
+// typedef Dune::PDELab::P1LocalFiniteElementMap<Scalar,Scalar,dim> type; // for simplices
+};
+#endif // HAVE_DUNE_PDELAB
+
+// Set the problem property
+SET_PROP(InfiltrationProblem, Problem)
+{
+ typedef Dumux::InfiltrationProblem<TypeTag> type;
+};
+
+// Set the wetting phase
+SET_TYPE_PROP(InfiltrationProblem, FluidSystem, Dumux::Holle_H2O_Air_Mesitylene_System<TypeTag>);
+
+// Set the spatial parameters
+SET_TYPE_PROP(InfiltrationProblem,
+ SpatialParameters,
+ Dumux::InfiltrationSpatialParameters<TypeTag>);
+
+// Enable gravity
+SET_BOOL_PROP(InfiltrationProblem, EnableGravity, true);
+
+// Use forward differences instead of central differences
+SET_INT_PROP(InfiltrationProblem, NumericDifferenceMethod, 0);
+
+// Write newton convergence
+SET_BOOL_PROP(InfiltrationProblem, NewtonWriteConvergence, true);
+
+//! Set the formulation
+SET_INT_PROP(InfiltrationProblem, Formulation, HolleThreePThreeCFormulation::pgSwSn);
+}
+
+
+/*!
+ * \ingroup HolleTHreePThreeCBoxModel
+ *
+ * */
+template <class TypeTag >
+class InfiltrationProblem : public HolleThreePThreeCProblem<TypeTag>
+{
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(Scalar)) Scalar;
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(GridView)) GridView;
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(Model)) Model;
+ typedef typename GridView::Grid Grid;
+
+ typedef InfiltrationProblem<TypeTag> ThisType;
+ typedef HolleThreePThreeCProblem<TypeTag> ParentType;
+
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(MaterialLaw)) MaterialLaw;
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(MaterialLawParams)) MaterialLawParams;
+
+
+ // copy some indices for convenience
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(HolleThreePThreeCIndices)) Indices;
+ enum {
+ numEq = GET_PROP_VALUE(TypeTag, PTAG(NumEq)),
+
+ pressureIdx = Indices::pressureIdx,
+ switch1Idx = Indices::switch1Idx,
+ switch2Idx = Indices::switch2Idx,
+// #if !ISOTHERMAL
+// temperatureIdx = Indices::temperatureIdx,
+// energyEqIdx = Indices::energyEqIdx,
+// #endif
+
+ // Phase State
+ threePhases = Indices::threePhases,
+ wPhaseOnly = Indices::wPhaseOnly,
+ gnPhaseOnly = Indices::gnPhaseOnly,
+ wnPhaseOnly = Indices::wnPhaseOnly,
+ gPhaseOnly = Indices::gPhaseOnly,
+ wgPhaseOnly = Indices::wgPhaseOnly,
+
+ // Grid and world dimension
+ dim = GridView::dimension,
+ dimWorld = GridView::dimensionworld,
+ };
+
+
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(PrimaryVariables)) PrimaryVariables;
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(BoundaryTypes)) BoundaryTypes;
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(TimeManager)) TimeManager;
+
+ typedef typename GridView::template Codim<0>::Entity Element;
+ typedef typename GridView::template Codim<dim>::Entity Vertex;
+ typedef typename GridView::Intersection Intersection;
+
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(FVElementGeometry)) FVElementGeometry;
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(FluidSystem)) FluidSystem;
+
+ typedef Dune::FieldVector<Scalar, dim> LocalPosition;
+ typedef Dune::FieldVector<Scalar, dimWorld> GlobalPosition;
+
+public:
+ /*!
+ * \brief The constructor
+ *
+ * \param timeManager The time manager
+ * \param gridView The grid view
+ */
+ InfiltrationProblem(TimeManager &timeManager, const GridView &gridView)
+ : ParentType(timeManager, gridView)
+ {
+ FluidSystem::init();
+ }
+
+ /*!
+ * \name Problem parameters
+ */
+ // \{
+
+ /*!
+ * \brief The problem name.
+ *
+ * This is used as a prefix for files generated by the simulation.
+ */
+ const char *name() const
+ { return "infiltration"; }
+
+ /*!
+ * \brief Returns the temperature within the domain.
+ *
+ * \param element The element
+ * \param fvElemGeom The finite-volume geometry in the box scheme
+ * \param scvIdx The local vertex index (SCV index)
+ *
+ * This problem assumes a temperature of 10 degrees Celsius.
+ */
+ Scalar temperature(const Element &element,
+ const FVElementGeometry &fvElemGeom,
+ int scvIdx) const
+ {
+ return (273.15 + 10.0); // -> Temperatur 10°C
+ };
+
+ // \}
+
+ /*!
+ * \name Boundary conditions
+ */
+ // \{
+
+ /*!
+ * \brief Specifies which kind of boundary condition should be
+ * used for which equation on a given boundary segment.
+ *
+ * \param values The boundary types for the conservation equations
+ * \param vertex The vertex for which the boundary type is set
+ */
+ void boundaryTypes(BoundaryTypes &values, const Vertex &vertex) const
+ {
+ const GlobalPosition globalPos = vertex.geometry().center();
+
+ if(globalPos[0] > 125. - eps_)
+ values.setAllDirichlet();
+ else if(globalPos[0] < eps_)
+ values.setAllDirichlet();
+ else
+ values.setAllNeumann();
+
+//#if !ISOTHERMAL
+// values.setDirichlet(temperatureIdx, energyEqIdx);
+//#endif
+ }
+
+ /*!
+ * \brief Evaluate the boundary conditions for a dirichlet
+ * boundary segment.
+ *
+ * \param values The dirichlet values for the primary variables
+ * \param vertex The vertex for which the boundary type is set
+ *
+ * For this method, the \a values parameter stores primary variables.
+ */
+ void dirichlet(PrimaryVariables &values, const Vertex &vertex) const
+ {
+ const GlobalPosition globalPos = vertex.geometry().center();
+
+ Scalar y = globalPos[1];
+ Scalar Sw, Swr=0.12, Sgr=0.03;
+
+ if(y >5. )
+ {
+ Scalar pc = 9.81 * 1000.0 * (y - 5); /*capillary Eq.*/
+ if (pc < 0.0) pc = 0.0;
+
+ Sw = invertPCGW_(pc,
+ this->spatialParameters().materialLawParams());
+ if (Sw < Swr) Sw = Swr;
+ if (Sw > 1.-Sgr) Sw = 1.-Sgr;
+
+ values[pressureIdx] = 1e5 ;
+ values[switch1Idx] = Sw;
+ values[switch2Idx] = 0.001;
+ }else {
+ values[pressureIdx] = 1e5 + 9.81 * 1000.0 * (5 - y);
+ values[switch1Idx] = 1.-Sgr;
+ values[switch2Idx] = 0.001;
+ }
+
+ //initial_(values, globalPos, element);
+ //const MaterialLawParams& materialParams = this->spatialParameters().materialLawParams();;
+ //MaterialLaw::pCGW(materialParams, 1.0);
+ }
+
+ /*!
+ * \brief Evaluate the boundary conditions for a neumann
+ * boundary segment.
+ *
+ * \param values The neumann values for the conservation equations
+ * \param element The finite element
+ * \param fvElemGeom The finite-volume geometry in the box scheme
+ * \param is The intersection between element and boundary
+ * \param scvIdx The local vertex index
+ * \param boundaryFaceIdx The index of the boundary face
+ *
+ * For this method, the \a values parameter stores the mass flux
+ * in normal direction of each phase. Negative values mean influx.
+ */
+ using ParentType::neumann;
+ void neumann(PrimaryVariables &values,
+ const Element &element,
+ const FVElementGeometry &fvElemGeom,
+ const Intersection &is,
+ int scvIdx,
+ int boundaryFaceIdx) const
+ {
+ const GlobalPosition &globalPos = element.geometry().corner(scvIdx);
+ values = 0;
+
+ // negative values for injection
+ if ((globalPos[0] <= 80.+eps_) && (globalPos[0] >= 75.+eps_) && (globalPos[1] >= 10.-eps_))
+ {
+ values[Indices::contiWEqIdx] = -0.0;
+ values[Indices::contiCEqIdx] = -1.2e-4;
+ values[Indices::contiAEqIdx] = -0.0;
+ }
+ }
+
+ // \}
+
+ /*!
+ * \name Volume terms
+ */
+ // \{
+
+ /*!
+ * \brief Evaluate the initial value for a control volume.
+ *
+ * \param values The initial values for the primary variables
+ * \param element The finite element
+ * \param fvElemGeom The finite-volume geometry in the box scheme
+ * \param scvIdx The local vertex index
+ *
+ * For this method, the \a values parameter stores primary
+ * variables.
+ */
+ void initial(PrimaryVariables &values,
+ const Element &element,
+ const FVElementGeometry &fvElemGeom,
+ int scvIdx) const
+ {
+ const GlobalPosition &globalPos = element.geometry().corner(scvIdx);
+
+ initial_(values, globalPos, element);
+ }
+
+ /*!
+ * \brief Return the initial phase state inside a control volume.
+ *
+ * \param vert The vertex
+ * \param globalIdx The index of the global vertex
+ * \param globalPos The global position
+ */
+ int initialPhasePresence(const Vertex &vert,
+ int &globalIdx,
+ const GlobalPosition &globalPos) const
+ {
+ return threePhases;
+ }
+
+private:
+ // internal method for the initial condition (reused for the
+ // dirichlet conditions!)
+ void initial_(PrimaryVariables &values,
+ const GlobalPosition &globalPos, const Element &element) const
+ {
+ Scalar y = globalPos[1];
+ Scalar Sw, Swr=0.12, Sgr=0.03;
+
+ if(y >5. )
+ {
+ Scalar pc = 9.81 * 1000.0 * (y - 5); /*capillary Eq.*/
+ if (pc < 0.0) pc = 0.0;
+
+ Sw = invertPCGW_(pc,
+ this->spatialParameters().materialLawParams());
+ if (Sw < Swr) Sw = Swr;
+ if (Sw > 1.-Sgr) Sw = 1.-Sgr;
+
+ values[pressureIdx] = 1e5 ;
+ values[switch1Idx] = Sw;
+ values[switch2Idx] = 0.001;
+ }else {
+ values[pressureIdx] = 1e5 + 9.81 * 1000.0 * (5 - y);
+ values[switch1Idx] = 1.-Sgr;
+ values[switch2Idx] = 0.001;
+ }
+ }
+
+ static Scalar invertPCGW_(Scalar pcIn, const MaterialLawParams &pcParams)
+ {
+ Scalar lower,upper;
+ int k;
+ int maxIt = 50;
+ Scalar bisLimit = 1.;
+ Scalar Sw, pcGW;
+ lower=0.0; upper=1.0;
+ for (k=1; k<=25; k++)
+ {
+ Sw = 0.5*(upper+lower);
+ pcGW = MaterialLaw::pCGW(pcParams, Sw);
+ Scalar delta = pcGW-pcIn;
+ if (delta<0.) delta*=-1.;
+ if (delta<bisLimit)
+ {
+ return(Sw);
+ }
+ if (k==maxIt) {
+ return(Sw);
+ }
+ if (pcGW>pcIn) lower=Sw;
+ else upper=Sw;
+ }
+ return(Sw);
+ }
+
+ static const Scalar eps_ = 1e-6;
+};
+} //end namespace
+
+#endif
--
GitLab