From d85c949c6d5ee51484dc64e778974f4a7ccd9d9a Mon Sep 17 00:00:00 2001
From: Michael Sinsbeck <michael@sinsbeck.com>
Date: Tue, 17 May 2011 09:44:00 +0000
Subject: [PATCH] DuMuX-equivalent for the old GRUWA-application
git-svn-id: svn://svn.iws.uni-stuttgart.de/DUMUX/dumux/trunk@5840 2fb0f335-1f38-0410-981e-8018bf24f1b0
---
.../mhs/groundwater/external_interface.hh | 458 +++++++++++++++++
appl/lecture/mhs/groundwater/groundwater.cc | 111 ++++
.../mhs/groundwater/groundwater_problem.hh | 480 ++++++++++++++++++
.../groundwater/groundwater_spatialparams.hh | 135 +++++
appl/lecture/mhs/groundwater/pseudoh2o.hh | 85 ++++
5 files changed, 1269 insertions(+)
create mode 100644 appl/lecture/mhs/groundwater/external_interface.hh
create mode 100644 appl/lecture/mhs/groundwater/groundwater.cc
create mode 100644 appl/lecture/mhs/groundwater/groundwater_problem.hh
create mode 100644 appl/lecture/mhs/groundwater/groundwater_spatialparams.hh
create mode 100644 appl/lecture/mhs/groundwater/pseudoh2o.hh
diff --git a/appl/lecture/mhs/groundwater/external_interface.hh b/appl/lecture/mhs/groundwater/external_interface.hh
new file mode 100644
index 0000000000..cad3cff335
--- /dev/null
+++ b/appl/lecture/mhs/groundwater/external_interface.hh
@@ -0,0 +1,458 @@
+// $Id: pcm_parameters.hh 1329 $
+/*****************************************************************************
+ * 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/>. *
+ *****************************************************************************/
+#ifndef GW_PARAMETERS_HH
+#define GW_PARAMETERS_HH
+
+#include <stdlib.h>
+#include <iostream>
+#include <fstream>
+
+/** \todo Please doc me! */
+namespace Dumux
+{
+
+struct Lens
+{
+ Dune::FieldMatrix<double,2,2> permeability;
+ //double permeability;
+ Dune::FieldVector<double,2> lowerLeft;
+ Dune::FieldVector<double,2> upperRight;
+};
+
+struct Source
+{
+ Dune::FieldVector<double,2> globalPos;
+ double q;
+ int index;
+};
+
+struct BoundaryCondition
+{
+ std::vector<bool> neumann;
+ std::vector<double> value;
+ std::vector<double> endPoint;
+ int segmentCount;
+};
+
+class InterfaceSoilProperties
+//Class for Interface Soil Properties
+//Integrate Parameters for Probabilistic Collocation Method (iPCM).
+//Contained design and uncertainty parameters
+{
+public:
+ //Interface Soil Properties (ISP):
+ double permeability;
+ double porosity;
+ std::vector <Lens> lenses;
+
+ InterfaceSoilProperties(const char* isp_filename)
+ //Initialization of ISP Parameters
+ {
+ using namespace std;
+ double viscosity, density;
+
+ std::cout << std::endl
+ << "-----> ISP: Interface Soil Properties Initialization ...\n";
+ //ISP input file defenition
+ ifstream input;
+
+ //ISP file check
+ input.open(isp_filename);
+ if (!input)
+ {
+ cout << "\n";
+ cout << "-----> ISP: Fatal error! - Data read \n";
+ cout << "-----> ISP: Could not open the input data file: \""
+ << isp_filename << "\n";
+ }
+
+ //iPCM input file reading
+ char reader[100]; // variable for input value
+ while (!input.eof())
+ {
+ input >> reader;
+ //if (reader==string("<SoilProperties>"))
+ if (reader == string("<Viscosity>"))
+ {
+ input >> reader;
+ viscosity = atof(reader);
+ }
+ if (reader == string("<Density>"))
+ {
+ input >> reader;
+ density = atof(reader);
+ }
+
+
+
+ //cout << "-----> ISP: Soil Properties reading ... \n";
+ //ISP parameters initialization:
+ if (reader == string("<GlobalPermeability>"))
+ {
+ input >> reader;
+ permeability = atof(reader);
+ cout << "-----> Global Permeability: " << permeability << "\n";
+ }
+ if (reader == string("<Porosity>"))
+ {
+ input >> reader;
+ porosity = atof(reader);
+ cout << "-----> Porosity: " << porosity << "\n";
+ }
+
+ if (reader == string("<Lens>"))
+ {
+ Lens templens;
+ templens.permeability=0;
+
+ // templens füllen
+ while (reader != string("</Lens>"))
+ {
+ input >> reader;
+ if (reader == string("<Permeability>"))
+ {
+ input >> reader;
+ templens.permeability[0][0] = atof(reader);
+ templens.permeability[1][1] = templens.permeability[0][0];
+ }
+ if (reader == string("<Left>"))
+ {
+ input >> reader;
+ templens.lowerLeft[0] = atof(reader);
+ }
+ if (reader == string("<Right>"))
+ {
+ input >> reader;
+ templens.upperRight[0] = atof(reader);
+ }
+ if (reader == string("<Bottom>"))
+ {
+ input >> reader;
+ templens.lowerLeft[1] = atof(reader);
+ }
+ if (reader == string("<Top>"))
+ {
+ input >> reader;
+ templens.upperRight[1] = atof(reader);
+ }
+ }
+ //testen: ist templens gut?
+ //an lenses-vector anhängen.
+ lenses.push_back(templens);
+ cout << "-----> Lens: Left, Bottom, Right, Top: " << templens.lowerLeft
+ << " " << templens.upperRight << ", Permeability; "<<
+ templens.permeability[0][0] << "\n";
+ }
+
+ }
+ input.close();
+ permeability*=(viscosity/(density*9.81));
+ for(int i=0; i< lenses.size();i++)
+ {
+ lenses[i].permeability[0][0]*=(viscosity/(density*9.81));
+ lenses[i].permeability[1][1]*=(viscosity/(density*9.81));
+ }
+
+ }
+
+};
+
+class InterfaceFluidProperties
+//Class for Interface Fluid Properties
+//Integrate Parameters for Probabilistic Collocation Method (iPCM).
+//Contained design and uncertainty parameters
+{
+public:
+ //Interface Fluid Properties (IFP):
+ float viscosity;// Gas Diffusion Coefficient
+ float density;// Residual Saturation of CO2
+
+ InterfaceFluidProperties(const char* ifp_filename)
+ //Initialization of IFP Parameters
+ {
+ using namespace std;
+ std::cout << std::endl
+ << "-----> IFP: Interface Fluid Properties Initialization ...\n";
+ //IFP input file defenition
+ ifstream input;
+
+ //IFP file check
+ input.open(ifp_filename);
+ if (!input)
+ {
+ cout << "\n";
+ cout << "-----> IFP: Fatal error! - Data read \n";
+ cout << "-----> IFP: Could not open the input data file: \""
+ << ifp_filename << "\n";
+ }
+
+ //iPCM input file reading
+ char reader[100]; // variable for input value
+ //double K;
+ while (!input.eof())
+ {
+ input >> reader;
+ //if (reader==string("<FluidProperties>"))
+ //cout << "-----> IFP: Fluid Properties reading ... \n";
+ //IFP perameters initialization:
+ if (reader == string("<Viscosity>"))
+ {
+ input >> reader;
+ viscosity = atof(reader);
+ cout << "-----> Viscosity: "
+ << viscosity << "\n";
+ }
+ if (reader == string("<Density>"))
+ {
+ input >> reader;
+ density = atof(reader);
+ cout << "-----> Density: "
+ << density << "\n";
+ }
+ }
+ input.close();
+ }
+
+};
+
+class InterfaceProblemProperties
+//Class for Interface Problem Properties
+//Integrate Parameters for Probabilistic Collocation Method (iPCM).
+//Contained design and uncertainty parameters
+{
+public:
+ //Interface Problem Properties (IPP):
+ Dune::FieldVector<int,2> resolution;
+ double depth;
+ Dune::FieldVector<double,2> size;
+ std::vector <Source> sources;
+ BoundaryCondition BCondition[4];
+ int plotMode;
+
+ InterfaceProblemProperties(const char* ipp_filename)
+ //Initialization of IPP Parameters
+ {
+ using namespace std;
+ std::cout
+ << "-----> IPP: Interface Problem Properties Initialization ...\n";
+ //IPP input file definition
+ ifstream input;
+
+ //IPP file check
+ input.open(ipp_filename);
+ if (!input)
+ {
+ cout << "\n";
+ cout << "-----> IPP: Fatal error! - Data read \n";
+ cout << "-----> IPP: Could not open the input data file: \""
+ << ipp_filename << "\n";
+ }
+
+ //iPCM input file reading
+ char reader[100]; // variable for input value
+ //double K;
+ while (!input.eof())
+ {
+ input >> reader;
+ //if (reader==string("<BoundaryAndInitialConditions>"))
+ //cout << "-----> IPP: Boundary and Initial Conditions reading ... \n";
+ //IPP perameters initialization:
+ if (reader == string("<XLength>"))
+ {
+ input >> reader;
+ size[0] = atof(reader);
+ cout << "-----> X-Length: " << size[0] << "\n";
+ }
+ if (reader == string("<YLength>"))
+ {
+ input >> reader;
+ size[1] = atof(reader);
+ cout << "-----> Y-Length: " << size[1] << "\n";
+ }
+ if (reader == string("<ZLength>"))
+ {
+ input >> reader;
+ depth = atof(reader);
+ cout << "-----> Z-Length: " << depth << "\n";
+ }
+ if (reader == string("<XResolution>"))
+ {
+ input >> reader;
+ resolution[0] = atoi(reader);
+ cout << "-----> X-Resolution: " << resolution[0] << "\n";
+ }
+ if (reader == string("<YResolution>"))
+ {
+ input >> reader;
+ resolution[1] = atoi(reader);
+ cout << "-----> Y-Resolution: " << resolution[1] << "\n";
+ }
+ if (reader == string("<PlotMode>"))
+ {
+ input >> reader;
+ plotMode = atoi(reader);
+ cout << "-----> Plot-mode: " << plotMode << "\n";
+ }
+ if (reader == string("<Source>"))
+ {
+ Source tempSource;
+ tempSource.q=0;
+
+ while (reader != string("</Source>"))
+ {
+ input >> reader;
+ if (reader == string("<X>"))
+ {
+ input >> reader;
+ tempSource.globalPos[0] = atof(reader);
+ }
+ if (reader == string("<Y>"))
+ {
+ input >> reader;
+ tempSource.globalPos[1] = atof(reader);
+ }
+ if (reader == string("<Q>"))
+ {
+ input >> reader;
+ tempSource.q = atof(reader);
+ }
+ }
+ //an sources-vector anhängen.
+ sources.push_back(tempSource);
+ cout << "-----> Sink/Source: Position: " << tempSource.globalPos
+ << " " << ", Q = "<< tempSource.q << "\n";
+ }
+
+ if (reader == string("<Boundary>"))
+ {
+
+ int boundaryIndex = -1;
+ while (boundaryIndex == -1)
+ {
+ input >> reader;
+ if (reader == string("<Top/>"))
+ boundaryIndex = 0;
+ if (reader == string("<Bottom/>"))
+ boundaryIndex = 1;
+ if (reader == string("<Left/>"))
+ boundaryIndex = 2;
+ if (reader == string("<Right/>"))
+ boundaryIndex = 3;
+ }
+
+ while (reader != string("</Boundary>"))
+ {
+ input >> reader;
+
+ if (reader == string("<Top/>"))
+ boundaryIndex = 0;
+ if (reader == string("<Bottom/>"))
+ boundaryIndex = 1;
+ if (reader == string("<Left/>"))
+ boundaryIndex = 2;
+ if (reader == string("<Right/>"))
+ boundaryIndex = 3;
+
+ if (reader == string("<Type>"))
+ {
+ input >> reader;
+ while (reader != string("</Type>"))
+ {
+ if (reader == string("neumann"))
+ BCondition[boundaryIndex].neumann.push_back(true);
+ else
+ BCondition[boundaryIndex].neumann.push_back(false);
+ input >> reader;
+ }
+ }
+
+ if (reader == string("<Value>"))
+ {
+ input >> reader;
+ while (reader != string("</Value>"))
+ {
+ BCondition[boundaryIndex].value.push_back(atof(reader));
+ input >> reader;
+ }
+ }
+
+ if (reader == string("<EndPoint>"))
+ {
+ input >> reader;
+ while (reader != string("</EndPoint>"))
+ {
+ BCondition[boundaryIndex].endPoint.push_back(atof(reader));
+ input >> reader;
+ }
+ }
+ }
+
+ for (int i=0; i<4 ; i++)
+ {
+ BCondition[i].segmentCount=std::min(std::min(
+ BCondition[i].value.size(),
+ BCondition[i].neumann.size()),
+ BCondition[i].endPoint.size()+1);
+ if (!BCondition[i].segmentCount)
+ {
+ BCondition[i].value.resize(1);
+ BCondition[i].value[0] = 0;
+ BCondition[i].neumann.resize(1);
+ BCondition[i].neumann[0] = true;
+ BCondition[i].endPoint.resize(0);
+ BCondition[i].segmentCount=1;
+// std::cout << "Invalid Boundary condition in Boundary "<<i<<". Set to no-flow."<<std::endl;
+ }
+
+ std::cout<< "-----> Boundary Conditions for ";
+ switch(i)
+ {
+ case 0: std::cout<< "top"; break;
+ case 1: std::cout<< "bottom"; break;
+ case 2: std::cout<< "left"; break;
+ case 3: std::cout<< "right";
+ }
+ std::cout << " Boundary:"<<std::endl;
+ for (int j=0; j<BCondition[i].segmentCount ;j++)
+ {
+ if (BCondition[i].neumann[j])
+ std::cout << " " << "Neumann ";
+ else
+ std::cout << " " << "Dirichlet ";
+ std::cout << BCondition[i].value[j] <<" ";
+ if (j < BCondition[i].segmentCount-1)
+ std::cout << BCondition[i].endPoint[j];
+ std::cout << std::endl;
+ }
+ }
+ }
+ }
+ input.close();
+
+ // Calculate for each source the containing element
+ for (int sourceNumber = 0; sourceNumber != sources.size(); sourceNumber++)
+ {
+ sources[sourceNumber].index = std::floor(sources[sourceNumber].globalPos[0]
+ * resolution[0]/size[0])
+ + std::floor(sources[sourceNumber].globalPos[1]*resolution[1]/size[1])
+ * resolution[0];
+ }
+ }
+
+};
+
+} // end namespace
+#endif
diff --git a/appl/lecture/mhs/groundwater/groundwater.cc b/appl/lecture/mhs/groundwater/groundwater.cc
new file mode 100644
index 0000000000..6d46fa8e71
--- /dev/null
+++ b/appl/lecture/mhs/groundwater/groundwater.cc
@@ -0,0 +1,111 @@
+// $Id: test_diffusion.cc 4144 2010-08-24 10:10:47Z bernd $
+/*****************************************************************************
+ * 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: <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
+ *
+ * \ingroup IMPETtests
+ * \brief test for the decoupled one-phase model.
+ */
+#include "config.h"
+#include <iostream>
+#include <boost/format.hpp>
+
+#include <dune/common/exceptions.hh>
+#include <dune/common/mpihelper.hh>
+#include <dune/grid/common/gridinfo.hh>
+
+#include "external_interface.hh"
+#include "groundwater_problem.hh"
+
+////////////////////////
+// the main function
+////////////////////////
+void usage(const char *progname)
+{
+ std::cout << boost::format("usage: %s #refine [delta]\n")%progname;
+ exit(1);
+}
+
+int main(int argc, char** argv)
+{
+ try {
+ typedef TTAG(GroundwaterProblem) TypeTag;
+ typedef GET_PROP_TYPE(TypeTag, PTAG(Scalar)) Scalar;
+ typedef GET_PROP_TYPE(TypeTag, PTAG(Grid)) Grid;
+ static const int dim = Grid::dimension;
+ typedef Dune::FieldVector<Scalar, dim> GlobalPosition;
+
+ // initialize MPI, finalize is done automatically on exit
+ Dune::MPIHelper::instance(argc, argv);
+
+ ////////////////////////////////////////////////////////////
+ // parse the command line arguments
+ ////////////////////////////////////////////////////////////
+// if (argc != 2 && argc != 3)
+// usage(argv[0]);
+//
+// int numRefine;
+// std::istringstream(argv[1]) >> numRefine;
+
+
+
+ ////////////////////////////////////////////////////////////
+ // Read Input file
+ ////////////////////////////////////////////////////////////
+
+ Dumux::InterfaceProblemProperties interfaceProbProps("interface_groundwater.xml");
+ Dumux::InterfaceFluidProperties interfaceFluidProps("interface_groundwater.xml");
+
+ ////////////////////////////////////////////////////////////
+ // create the grid
+ ////////////////////////////////////////////////////////////
+ GlobalPosition L(0.0);
+ Grid grid(interfaceProbProps.resolution,L,interfaceProbProps.size);
+
+ ////////////////////////////////////////////////////////////
+ // adjust fluid properties
+ ////////////////////////////////////////////////////////////
+ typedef GET_PROP_TYPE(TypeTag, PTAG(Fluid)) Fluid;
+ Fluid::Component::setViscosity(interfaceFluidProps.viscosity);
+ Fluid::Component::setDensity(interfaceFluidProps.density);
+
+ ////////////////////////////////////////////////////////////
+ // instantiate and run the concrete problem
+ ////////////////////////////////////////////////////////////
+ typedef GET_PROP_TYPE(TypeTag, PTAG(Problem)) Problem;
+ Problem problem(grid.leafView(), interfaceProbProps);
+ problem.init();
+ problem.writeOutput();
+
+ return 0;
+ }
+ catch (Dune::Exception &e) {
+ std::cerr << "Dune reported error: " << e << std::endl;
+ }
+ catch (...) {
+ std::cerr << "Unknown exception thrown!\n";
+ throw;
+ }
+
+ return 3;
+}
diff --git a/appl/lecture/mhs/groundwater/groundwater_problem.hh b/appl/lecture/mhs/groundwater/groundwater_problem.hh
new file mode 100644
index 0000000000..48f8ef6625
--- /dev/null
+++ b/appl/lecture/mhs/groundwater/groundwater_problem.hh
@@ -0,0 +1,480 @@
+// $Id: test_diffusion_problem.hh 3655 2010-05-26 17:13:50Z bernd $
+/*****************************************************************************
+* Copyright (C) 2007-2008 by Markus Wolff *
+* 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 Dumux-equivalent for GRUWA (1p-stationary, finite volumes)
+ */
+#ifndef DUMUX_GROUNDWATER_PROBLEM_HH
+#define DUMUX_GROUNDWATER_PROBLEM_HH
+
+#if HAVE_UG
+#include <dune/grid/uggrid.hh>
+#endif
+
+#include <dune/grid/yaspgrid.hh>
+#include <dune/grid/sgrid.hh>
+
+#include <dumux/material/fluidsystems/liquidphase.hh>
+//#include <dumux/material/components/unit.hh>
+#include "pseudoh2o.hh"
+
+#include <dumux/decoupled/1p/diffusion/diffusionproblem1p.hh>
+#include <dumux/decoupled/1p/diffusion/fv/fvvelocity1p.hh>
+
+#include "groundwater_spatialparams.hh"
+#include "vipWriter.hh"
+
+namespace Dumux
+{
+
+template<class TypeTag>
+class GroundwaterProblem;
+
+//////////
+// Specify the properties
+//////////
+namespace Properties
+{
+NEW_TYPE_TAG(GroundwaterProblem, INHERITS_FROM(DecoupledOneP))
+ ;
+
+// Set the grid type
+SET_PROP(GroundwaterProblem, Grid)
+{// typedef Dune::YaspGrid<2> type;
+ typedef Dune::SGrid<2, 2> type;
+};
+
+// Set the wetting phase
+SET_PROP(GroundwaterProblem, Fluid)
+{
+private:
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(Scalar)) Scalar;
+public:
+ typedef Dumux::LiquidPhase<Scalar, Dumux::PseudoH2O<Scalar> > type;
+};
+
+// Set the spatial parameters
+SET_PROP(GroundwaterProblem, SpatialParameters)
+{
+private:
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(Grid)) Grid;
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(Scalar)) Scalar;
+
+public:
+ typedef Dumux::GroundwaterSpatialParams<TypeTag> type;
+};
+
+// Enable gravity
+SET_BOOL_PROP(GroundwaterProblem, EnableGravity, false);
+
+// Set the model
+SET_TYPE_PROP(GroundwaterProblem, Model, Dumux::FVVelocity1P<TypeTag>);
+
+//Set the problem
+SET_TYPE_PROP(GroundwaterProblem, Problem, Dumux::GroundwaterProblem<TTAG(GroundwaterProblem)>);
+
+}
+
+/*!
+ * \ingroup IMPETtests
+ *
+ * \brief test problem for the decoupled one-phase model.
+ */
+template<class TypeTag = TTAG(GroundwaterProblem)>
+class GroundwaterProblem: public DiffusionProblem1P<TypeTag, GroundwaterProblem<TypeTag> >
+{
+ typedef GroundwaterProblem<TypeTag> ThisType;
+ typedef DiffusionProblem1P<TypeTag, ThisType> ParentType;
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(GridView)) GridView;
+
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(Fluid)) Fluid;
+
+ enum
+ {
+ dim = GridView::dimension, dimWorld = GridView::dimensionworld
+ };
+
+ 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<Scalar, dimWorld> GlobalPosition;
+ typedef typename GridView::template Codim<0>::Iterator ElementIterator;
+
+public:
+ GroundwaterProblem(const GridView &gridView, Dumux::InterfaceProblemProperties interfaceProbProps) :
+ ParentType(gridView), problemProps_(interfaceProbProps)
+ {
+// this->spatialParameters().setDelta(delta_);
+
+
+ }
+
+ /*!
+ * \name Problem parameters
+ */
+ // \{
+
+ /*!
+ * \brief The problem name.
+ *
+ * This is used as a prefix for files generated by the simulation.
+ */
+ const char *name() const
+ {
+ return "groundwater";
+ }
+
+ 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
+ }
+
+ // \}
+
+ //! Returns the reference pressure for evaluation of constitutive relations
+ Scalar referencePressure(const GlobalPosition& globalPos, const Element& element) const
+ {
+ return 1e5; // -> 10°C
+ }
+
+ //!source term [kg/(m^3 s)] (or in 2D [kg/(m^2 s)]).
+ Scalar source(const GlobalPosition& globalPos, const Element& element)
+ {
+ Scalar q=0;
+ Scalar density=Fluid::density(0,0);
+ for (int sourceNumber = 0; sourceNumber != problemProps_.sources.size(); sourceNumber++)
+ {
+ Source source = problemProps_.sources[sourceNumber];
+ if (this->variables().index(element) == source.index)
+ q+=source.q*density/element.geometry().volume()/problemProps_.depth;
+ }
+ return q;
+ }
+
+ /*!
+ * \brief Returns the type of boundary condition.
+ *
+ * BC can be dirichlet (pressure) or neumann (flux).
+ */
+ typename BoundaryConditions::Flags bctype(const GlobalPosition& globalPos, const Intersection& intersection) const
+ {
+ double coordinate=0;
+ int boundaryIndex=0;
+ if (globalPos[0]<0.0001)
+ {
+ //Left boundary
+ coordinate=globalPos[1];
+ boundaryIndex=2;
+ }
+ if (globalPos[1]<0.0001)
+ {
+ //Bottom boundary
+ coordinate=globalPos[0];
+ boundaryIndex=1;
+ }
+ if (globalPos[0]> problemProps_.size[0] -0.0001)
+ {
+ //Right boundary
+ coordinate=globalPos[1];
+ boundaryIndex=3;
+ }
+ if (globalPos[1]> problemProps_.size[1] -0.0001)
+ {
+ //Top boundary
+ coordinate=globalPos[0];
+ boundaryIndex=0;
+ }
+ for (int segmentCount=0; segmentCount<problemProps_.BCondition[boundaryIndex].segmentCount-1;segmentCount++)
+ {
+ if (coordinate< problemProps_.BCondition[boundaryIndex].endPoint[segmentCount])
+ {
+ if (problemProps_.BCondition[boundaryIndex].neumann[segmentCount])
+ return BoundaryConditions::neumann;
+ else
+ return BoundaryConditions::dirichlet;
+ }
+ }
+ if (problemProps_.BCondition[boundaryIndex].neumann[problemProps_.BCondition[boundaryIndex].segmentCount-1])
+ return BoundaryConditions::neumann;
+ else
+ return BoundaryConditions::dirichlet;
+ }
+
+ //! return dirichlet condition (pressure, [Pa])
+ Scalar dirichlet(const GlobalPosition& globalPos, const Intersection& intersection) const
+ {
+ double coordinate=0;
+ int boundaryIndex=0;
+ if (globalPos[0]<0.0001)
+ {
+ //Left boundary
+ coordinate=globalPos[1];
+ boundaryIndex=2;
+ }
+ if (globalPos[1]<0.0001)
+ {
+ //Bottom boundary
+ coordinate=globalPos[0];
+ boundaryIndex=1;
+ }
+ if (globalPos[0]> problemProps_.size[0] -0.0001)
+ {
+ //Right boundary
+ coordinate=globalPos[1];
+ boundaryIndex=3;
+ }
+ if (globalPos[1]> problemProps_.size[1] -0.0001)
+ {
+ //Top boundary
+ coordinate=globalPos[0];
+ boundaryIndex=0;
+ }
+ for (int segmentCount=0; segmentCount<problemProps_.BCondition[boundaryIndex].segmentCount-1;segmentCount++)
+ {
+ if (coordinate< problemProps_.BCondition[boundaryIndex].endPoint[segmentCount])
+ {
+ return problemProps_.BCondition[boundaryIndex].value[segmentCount]*Fluid::density(0,0)*9.81;
+ }
+ }
+ return problemProps_.BCondition[boundaryIndex].value[problemProps_.BCondition[boundaryIndex].segmentCount-1]*Fluid::density(0,0)*9.81;
+
+
+// return 2e5;
+// return (exact(globalPos));
+ }
+
+ //! return neumann condition (flux, [kg/(m^2 s)])
+ Scalar neumann(const GlobalPosition& globalPos, const Intersection& intersection) const
+ {
+ double coordinate=0;
+ int boundaryIndex=0;
+ if (globalPos[0]<0.0001)
+ {
+ //Left boundary
+ coordinate=globalPos[1];
+ boundaryIndex=2;
+ }
+ if (globalPos[1]<0.0001)
+ {
+ //Bottom boundary
+ coordinate=globalPos[0];
+ boundaryIndex=1;
+ }
+ if (globalPos[0]> problemProps_.size[0] -0.0001)
+ {
+ //Right boundary
+ coordinate=globalPos[1];
+ boundaryIndex=3;
+ }
+ if (globalPos[1]> problemProps_.size[1] -0.0001)
+ {
+ //Top boundary
+ coordinate=globalPos[0];
+ boundaryIndex=0;
+ }
+ for (int segmentCount=0; segmentCount<problemProps_.BCondition[boundaryIndex].segmentCount-1;segmentCount++)
+ {
+ if (coordinate< problemProps_.BCondition[boundaryIndex].endPoint[segmentCount])
+ {
+ return problemProps_.BCondition[boundaryIndex].value[segmentCount]*Fluid::density(0,0)
+ *(-1);
+ }
+ }
+ return problemProps_.BCondition[boundaryIndex].value[problemProps_.BCondition[boundaryIndex].segmentCount-1]
+ *Fluid::density(0,0)*(-1);
+
+
+// if (globalPos[1]>0.999 && globalPos[0]>0.5)
+// return 1;
+// return 0.0;
+ }
+
+ void writeOutput()
+ { //Achtung: Quick und dirty:
+ // ViPWriter vipWriter(*this, problemProps_.resolution[0], problemProps_.resolution[1]);
+ switch(problemProps_.plotMode)
+ {
+ case 0: //grid-plot (colors only)
+ {
+ std::ofstream dataFile;
+ dataFile.open("output.vgf");
+ dataFile << "Gridplot" << std::endl;
+ dataFile << "## This is an DuMuX output for the NumLab Grafics driver. \n";
+ dataFile << "## This output file was generated at " << __TIME__ <<", "<< __DATE__<< "\n";
+ // the following specifies necessary information for the numLab grafical ouptut
+ dataFile << "# x-range 0 "<< problemProps_.size[0] << "\n" ;
+ dataFile << "# y-range 0 "<< problemProps_.size[1] << "\n" ;
+ dataFile << "# x-count " << problemProps_.resolution[0] << "\n" ;
+ dataFile << "# y-count " << problemProps_.resolution[1] << "\n" ;
+ dataFile << "# min-color 255 0 0\n";
+ dataFile << "# max-color 0 0 255\n";
+
+ dataFile << "# time 0 \n" ;
+
+ dataFile << "# label piezometric head \n";
+
+
+// for (int i=problemProps_.resolution[1]-1; i>=0; i--)
+ for (int i=0; i< problemProps_.resolution[1]; i++)
+ {
+ for (int j=problemProps_.resolution[0]-1; j>=0; j--)
+ {
+ int currentIdx = i*problemProps_.resolution[0]+j;
+ dataFile << this->variables().pressure()[currentIdx]/(Fluid::density(0,0)*9.81);
+ if(j != 0) // all but last entry
+ dataFile << " ";
+ else // write the last entry
+ dataFile << "\n";
+ }
+ }
+ dataFile.close();
+ }
+ break;
+ case 1: //Piecewise constant with 3d-plotter
+ {
+ std::ofstream dataFile;
+ dataFile.open("output.dat");
+
+ ElementIterator eItEnd = this->gridView().template end<0> ();
+ for (ElementIterator eIt = this->gridView().template begin<0> (); eIt != eItEnd; ++eIt)
+ {
+ //Aktuell gibts noch probleme mit der Skalierung.
+ int cellIndex = this->variables().index(*eIt);
+ dataFile << eIt->geometry().corner(0)[0] << " "
+ << -1*eIt->geometry().corner(0)[1] << " "
+ << this->variables().pressure()[cellIndex]/(Fluid::density(0,0)*9.81)<< " "
+ << eIt->geometry().corner(1)[0] << " "
+ << -1*eIt->geometry().corner(1)[1] << " "
+ << this->variables().pressure()[cellIndex]/(Fluid::density(0,0)*9.81)<< " "
+ << eIt->geometry().corner(2)[0] << " "
+ << -1*eIt->geometry().corner(2)[1] << " "
+ << this->variables().pressure()[cellIndex]/(Fluid::density(0,0)*9.81)<<std::endl;
+
+ dataFile << eIt->geometry().corner(3)[0] << " "
+ << -1*eIt->geometry().corner(3)[1] << " "
+ << this->variables().pressure()[cellIndex]/(Fluid::density(0,0)*9.81)<< " "
+ << eIt->geometry().corner(1)[0] << " "
+ << -1*eIt->geometry().corner(1)[1] << " "
+ << this->variables().pressure()[cellIndex]/(Fluid::density(0,0)*9.81)<< " "
+ << eIt->geometry().corner(2)[0] << " "
+ << -1*eIt->geometry().corner(2)[1] << " "
+ << this->variables().pressure()[cellIndex]/(Fluid::density(0,0)*9.81)<<std::endl;
+ }
+ }
+ break;
+ case 2: //Interpolate between cell centers
+ {
+ std::ofstream dataFile;
+ dataFile.open("interpolate.dat");
+
+ for (int i=0; i< problemProps_.resolution[1]-1; i++)
+ {
+ for (int j=0; j< problemProps_.resolution[0]-1; j++)
+ {
+ int currentIdx = i*problemProps_.resolution[0]+j;
+ double dx=problemProps_.size[0]/problemProps_.resolution[0];
+ double dy=problemProps_.size[1]/problemProps_.resolution[1];
+
+ dataFile << dx*(j+0.5) << " "
+ << -dy*(i+0.5) << " "
+ << this->variables().pressure()[currentIdx]/(Fluid::density(0,0)*9.81)<< " "
+ << dx*(j+1.5) << " "
+ << -dy*(i+0.5) << " "
+ << this->variables().pressure()[currentIdx+1]/(Fluid::density(0,0)*9.81)<< " "
+ << dx*(j+1.5) << " "
+ << -dy*(i+1.5) << " "
+ << this->variables().pressure()[currentIdx+1+problemProps_.resolution[0]]/(Fluid::density(0,0)*9.81)<<std::endl;
+
+ dataFile << dx*(j+0.5) << " "
+ << -dy*(i+0.5) << " "
+ << this->variables().pressure()[currentIdx]/(Fluid::density(0,0)*9.81)<< " "
+ << dx*(j+0.5) << " "
+ << -dy*(i+1.5) << " "
+ << this->variables().pressure()[currentIdx+problemProps_.resolution[0]]/(Fluid::density(0,0)*9.81)<< " "
+ << dx*(j+1.5) << " "
+ << -dy*(i+1.5) << " "
+ << this->variables().pressure()[currentIdx+1+problemProps_.resolution[0]]/(Fluid::density(0,0)*9.81)<<std::endl;
+ }
+ }
+ dataFile.close();
+ }
+ break;
+ default:
+ {
+ //Textoutput:
+ std::cout << "x, y, piezometric head, v_x, v_y"<<std::endl;
+ ElementIterator eItEnd = this->gridView().template end<0> ();
+ for (ElementIterator eIt = this->gridView().template begin<0> (); eIt != eItEnd; ++eIt)
+ {
+ int cellIndex = this->variables().index(*eIt);
+ double v_x,v_y,piezo;
+ v_x= (this->variables().velocity()[cellIndex][0][0]+this->variables().velocity()[cellIndex][1][0])/2;
+ v_y= (this->variables().velocity()[cellIndex][2][1]+this->variables().velocity()[cellIndex][3][1])/2;
+
+ if (std::abs(v_x)<1e-17)
+ v_x=0;
+ if (std::abs(v_y)<1e-17)
+ v_y=0;
+ piezo=this->variables().pressure()[cellIndex]/(Fluid::density(0,0)*9.81);
+
+ std::cout << std::setw(10)<< eIt->geometry().center()[0]<<" "
+ << std::setw(10)<< eIt->geometry().center()[1]<<" "
+ << std::setw(10) <<std::setprecision(6) << piezo<<" "
+ << std::setw(11)<< std::setprecision(6) << v_x<<" "
+ << std::setw(11)<< std::setprecision(6) << v_y <<std::endl;
+ }
+ }
+ }
+ }
+
+private:
+ Scalar exact (const GlobalPosition& globalPos) const
+ {
+ double pi = 4.0*atan(1.0);
+
+ return (sin(pi*globalPos[0])*sin(pi*globalPos[1]));
+ }
+
+ Dune::FieldVector<Scalar,dim> exactGrad (const GlobalPosition& globalPos) const
+ {
+ Dune::FieldVector<Scalar,dim> grad(0);
+ double pi = 4.0*atan(1.0);
+ grad[0] = pi*cos(pi*globalPos[0])*sin(pi*globalPos[1]);
+ grad[1] = pi*cos(pi*globalPos[1])*sin(pi*globalPos[0]);
+
+ return grad;
+ }
+
+ Dumux::InterfaceProblemProperties problemProps_;
+};
+} //end namespace
+
+#endif
diff --git a/appl/lecture/mhs/groundwater/groundwater_spatialparams.hh b/appl/lecture/mhs/groundwater/groundwater_spatialparams.hh
new file mode 100644
index 0000000000..bc346e0fc4
--- /dev/null
+++ b/appl/lecture/mhs/groundwater/groundwater_spatialparams.hh
@@ -0,0 +1,135 @@
+// $Id: test_2p_spatialparams.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: <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 spatial parameters for the GRUWA-equivalent exercise
+ */
+#ifndef GROUNDWATER_SPATIALPARAMETERS_HH
+#define GROUNDWATER_SPATIALPARAMETERS_HH
+
+
+#include <dumux/material/fluidmatrixinteractions/2p/linearmaterial.hh>
+#include <dumux/material/fluidmatrixinteractions/2p/efftoabslaw.hh>
+
+namespace Dumux
+{
+
+/*!
+ * \ingroup IMPETtests
+ * \brief spatial parameters for the test problem for diffusion models.
+ */
+template<class TypeTag>
+class GroundwaterSpatialParams
+{
+ 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<CoordScalar, dimWorld> GlobalPosition;
+ typedef Dune::FieldVector<CoordScalar, dim> LocalPosition;
+ typedef Dune::FieldMatrix<Scalar,dim,dim> FieldMatrix;
+
+ typedef LinearMaterial<Scalar> RawMaterialLaw;
+public:
+ typedef EffToAbsLaw<RawMaterialLaw> MaterialLaw;
+ typedef typename MaterialLaw::Params MaterialLawParams;
+
+ void update (Scalar saturationW, const Element& element)
+ {
+
+ }
+
+ const FieldMatrix& intrinsicPermeability (const GlobalPosition& globalPos, const Element& element) const
+ {
+ if (lenses_.size())
+ {
+ for (int lensnumber = 0; lensnumber != lenses_.size(); lensnumber++)
+ {
+ if ((lenses_[lensnumber].lowerLeft[0] < globalPos[0])
+ && (lenses_[lensnumber].lowerLeft[1] < globalPos[1])
+ && (lenses_[lensnumber].upperRight[0] > globalPos[0])
+ && (lenses_[lensnumber].upperRight[1] > globalPos[1]))
+ {
+ return lenses_[lensnumber].permeability;
+ }
+
+
+ }
+ }
+ // Schleife über alle Linsen
+ return permeability_;
+ }
+
+ double porosity(const GlobalPosition& globalPos, const Element& element) const
+ {
+ return porosity_;
+ }
+
+
+ // return the parameter object for the Brooks-Corey material law which depends on the position
+ const MaterialLawParams& materialLawParams(const GlobalPosition& globalPos, const Element &element) const
+ {
+ return materialLawParams_;
+ }
+
+ void setDelta(const double delta)
+ {
+ delta_ = delta;
+ }
+
+ GroundwaterSpatialParams(const GridView& gridView)
+ : permeability_(0)
+ {
+ delta_=1e-3;
+ Dumux::InterfaceSoilProperties interfaceSoilProps("interface_groundwater.xml");
+ porosity_ = interfaceSoilProps.porosity;
+ permeability_[0][0] = interfaceSoilProps.permeability;
+ permeability_[0][1] = 0;
+ permeability_[1][0] = 0;
+ permeability_[1][1] = interfaceSoilProps.permeability;
+
+ lenses_ = interfaceSoilProps.lenses;
+
+ // residual saturations
+ materialLawParams_.setSwr(0.0);
+ materialLawParams_.setSnr(0.0);
+
+ // parameters for the linear entry pressure function
+ materialLawParams_.setEntryPC(0);
+ materialLawParams_.setMaxPC(0);
+ }
+
+private:
+ MaterialLawParams materialLawParams_;
+ mutable FieldMatrix permeability_;
+ Scalar porosity_;
+ double delta_;
+ std::vector <Lens> lenses_;
+};
+
+} // end namespace
+#endif
diff --git a/appl/lecture/mhs/groundwater/pseudoh2o.hh b/appl/lecture/mhs/groundwater/pseudoh2o.hh
new file mode 100644
index 0000000000..0fe8754bb8
--- /dev/null
+++ b/appl/lecture/mhs/groundwater/pseudoh2o.hh
@@ -0,0 +1,85 @@
+/*****************************************************************************
+ * Copyright (C) 2009 by Andreas Lauser
+ * 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 Properties of pure water \f$H_2O\f$.
+ */
+#ifndef DUMUX_PSEUDOH2O_HH
+#define DUMUX_PSEUDOH2O_HH
+
+
+#include <dumux/material/components/component.hh>
+//#include "interface_BL.xml"
+
+namespace Dumux
+{
+/*!
+ * \brief Rough estimate for testing purposes of some water.
+ */
+template <class ScalarT>
+class PseudoH2O : public Component<ScalarT, PseudoH2O<ScalarT> >
+{
+ typedef Component<ScalarT, PseudoH2O<ScalarT> > ParentType;
+ typedef ScalarT Scalar;
+public:
+ /*!
+ * \brief A human readable name for the water.
+ */
+ static const char *name()
+ { return "H2O"; }
+
+ /*!
+ * \brief Rough estimate of the density of oil [kg/m^3].
+ */
+ static Scalar liquidDensity(Scalar temperature, Scalar pressure)
+ {
+ return density_;
+ }
+
+ /*!
+ * \brief Rough estimate of the viscosity of oil kg/(ms).
+ */
+ static Scalar liquidViscosity(Scalar temperature, Scalar pressure)
+ {
+ return viscosity_;
+ };
+
+ static void setViscosity(Scalar viscosity)
+ {
+ viscosity_ = viscosity;
+ }
+
+ static void setDensity(Scalar density)
+ {
+ density_ = density;
+ }
+
+private:
+ static Scalar viscosity_;
+ static Scalar density_;
+};
+template <class ScalarT>
+typename PseudoH2O<ScalarT>::Scalar PseudoH2O<ScalarT>::viscosity_ = 0.001;
+template <class ScalarT>
+typename PseudoH2O<ScalarT>::Scalar PseudoH2O<ScalarT>::density_ = 1000;
+} // end namepace
+
+#endif
--
GitLab