From 07920f5877feaf6b6e501f0ce11e89caa64e573a Mon Sep 17 00:00:00 2001
From: Klaus Mosthaf <klmos@env.dtu.dk>
Date: Tue, 19 Oct 2010 09:01:58 +0000
Subject: [PATCH] migrated appl/lecture/msm to dumux (stable part)
git-svn-id: svn://svn.iws.uni-stuttgart.de/DUMUX/dumux/trunk@4458 2fb0f335-1f38-0410-981e-8018bf24f1b0
---
appl/lecture/Makefile.am | 3 +
appl/lecture/msm/1p2cvs2p/Makefile.am | 11 +
.../msm/1p2cvs2p/external_interface.hh | 377 +++++++++++++++
appl/lecture/msm/1p2cvs2p/interface1p2c.xml | 54 +++
appl/lecture/msm/1p2cvs2p/interface2p.xml | 63 +++
appl/lecture/msm/1p2cvs2p/lens_1p2c.cc | 196 ++++++++
.../msm/1p2cvs2p/lens_1p2cnewtoncontroller.hh | 120 +++++
appl/lecture/msm/1p2cvs2p/lens_2p.cc | 196 ++++++++
.../msm/1p2cvs2p/lens_2pnewtoncontroller.hh | 86 ++++
appl/lecture/msm/1p2cvs2p/lensproblem1p2c.hh | 384 +++++++++++++++
appl/lecture/msm/1p2cvs2p/lensproblem2p.hh | 432 +++++++++++++++++
.../msm/1p2cvs2p/lensspatialparameters1p2c.hh | 184 +++++++
.../msm/1p2cvs2p/lensspatialparameters2p.hh | 191 ++++++++
.../lecture/msm/1p2cvs2p/water_contaminant.hh | 142 ++++++
appl/lecture/msm/Makefile.am | 8 +
appl/lecture/msm/buckleyleverett/Makefile.am | 8 +
.../buckleyleverett_analytic.hh | 377 +++++++++++++++
.../msm/buckleyleverett/buckleyleverett_ff.cc | 152 ++++++
.../buckleyleverett_spatialparams.hh | 111 +++++
.../buckleyleverett/buckleyleverettproblem.hh | 317 ++++++++++++
.../msm/buckleyleverett/external_interface.hh | 451 ++++++++++++++++++
.../msm/buckleyleverett/interface_BL.xml | 48 ++
appl/lecture/msm/buckleyleverett/pseudoh2o.hh | 80 ++++
appl/lecture/msm/buckleyleverett/pseudooil.hh | 83 ++++
appl/lecture/msm/mcwhorter/Makefile.am | 7 +
.../msm/mcwhorter/external_interface.hh | 443 +++++++++++++++++
appl/lecture/msm/mcwhorter/interface_MW.xml | 51 ++
.../msm/mcwhorter/mcwhorter_analytic.hh | 390 +++++++++++++++
appl/lecture/msm/mcwhorter/mcwhorter_ff.cc | 81 ++++
.../msm/mcwhorter/mcwhorter_spatialparams.hh | 112 +++++
.../lecture/msm/mcwhorter/mcwhorterproblem.hh | 265 ++++++++++
31 files changed, 5423 insertions(+)
create mode 100755 appl/lecture/Makefile.am
create mode 100755 appl/lecture/msm/1p2cvs2p/Makefile.am
create mode 100644 appl/lecture/msm/1p2cvs2p/external_interface.hh
create mode 100755 appl/lecture/msm/1p2cvs2p/interface1p2c.xml
create mode 100755 appl/lecture/msm/1p2cvs2p/interface2p.xml
create mode 100755 appl/lecture/msm/1p2cvs2p/lens_1p2c.cc
create mode 100644 appl/lecture/msm/1p2cvs2p/lens_1p2cnewtoncontroller.hh
create mode 100755 appl/lecture/msm/1p2cvs2p/lens_2p.cc
create mode 100644 appl/lecture/msm/1p2cvs2p/lens_2pnewtoncontroller.hh
create mode 100644 appl/lecture/msm/1p2cvs2p/lensproblem1p2c.hh
create mode 100644 appl/lecture/msm/1p2cvs2p/lensproblem2p.hh
create mode 100644 appl/lecture/msm/1p2cvs2p/lensspatialparameters1p2c.hh
create mode 100644 appl/lecture/msm/1p2cvs2p/lensspatialparameters2p.hh
create mode 100644 appl/lecture/msm/1p2cvs2p/water_contaminant.hh
create mode 100755 appl/lecture/msm/Makefile.am
create mode 100644 appl/lecture/msm/buckleyleverett/Makefile.am
create mode 100644 appl/lecture/msm/buckleyleverett/buckleyleverett_analytic.hh
create mode 100644 appl/lecture/msm/buckleyleverett/buckleyleverett_ff.cc
create mode 100644 appl/lecture/msm/buckleyleverett/buckleyleverett_spatialparams.hh
create mode 100644 appl/lecture/msm/buckleyleverett/buckleyleverettproblem.hh
create mode 100644 appl/lecture/msm/buckleyleverett/external_interface.hh
create mode 100755 appl/lecture/msm/buckleyleverett/interface_BL.xml
create mode 100644 appl/lecture/msm/buckleyleverett/pseudoh2o.hh
create mode 100644 appl/lecture/msm/buckleyleverett/pseudooil.hh
create mode 100644 appl/lecture/msm/mcwhorter/Makefile.am
create mode 100644 appl/lecture/msm/mcwhorter/external_interface.hh
create mode 100755 appl/lecture/msm/mcwhorter/interface_MW.xml
create mode 100644 appl/lecture/msm/mcwhorter/mcwhorter_analytic.hh
create mode 100644 appl/lecture/msm/mcwhorter/mcwhorter_ff.cc
create mode 100644 appl/lecture/msm/mcwhorter/mcwhorter_spatialparams.hh
create mode 100644 appl/lecture/msm/mcwhorter/mcwhorterproblem.hh
diff --git a/appl/lecture/Makefile.am b/appl/lecture/Makefile.am
new file mode 100755
index 0000000000..571f03eca0
--- /dev/null
+++ b/appl/lecture/Makefile.am
@@ -0,0 +1,3 @@
+SUBDIRS = msm
+
+include $(top_srcdir)/am/global-rules
diff --git a/appl/lecture/msm/1p2cvs2p/Makefile.am b/appl/lecture/msm/1p2cvs2p/Makefile.am
new file mode 100755
index 0000000000..c0e2d6f9d8
--- /dev/null
+++ b/appl/lecture/msm/1p2cvs2p/Makefile.am
@@ -0,0 +1,11 @@
+check_PROGRAMS = lens_1p2c lens_2p
+
+lens_1p2c_SOURCES = lens_1p2c.cc
+lens_1p2c_CXXFLAGS = $(MPI_CPPFLAGS)
+lens_1p2c_LDADD = $(MPI_LDFLAGS)
+
+lens_2p_SOURCES = lens_2p.cc
+lens_2p_CXXFLAGS = $(MPI_CPPFLAGS)
+lens_2p_LDADD = $(MPI_LDFLAGS)
+
+include $(top_srcdir)/am/global-rules
diff --git a/appl/lecture/msm/1p2cvs2p/external_interface.hh b/appl/lecture/msm/1p2cvs2p/external_interface.hh
new file mode 100644
index 0000000000..65a806dca7
--- /dev/null
+++ b/appl/lecture/msm/1p2cvs2p/external_interface.hh
@@ -0,0 +1,377 @@
+// $Id: pcm_parameters.hh 1329 $
+#ifndef PCM_PARAMETERS_HH
+#define PCM_PARAMETERS_HH
+
+#include <stdlib.h>
+#include <iostream>
+#include <fstream>
+
+/** \todo Please doc me! */
+namespace Dumux
+{
+
+class InterfaceSoilProperties
+//Class for Interface Soil Properties
+//Integrate Parameters for Probabilistic Collocation Method (iPCM).
+//Contained design and uncertainty parameters
+{
+public:
+ //Interface Soil Properties (ISP):
+ float ISP_Permeability; // Permeability
+ float ISP_FinePermeability; // Fine Permeability
+ float ISP_CoarsePermeability; // Coarse Permeability
+ float ISP_Porosity; // Porosity
+ float ISP_FinePorosity; // Fine Porosity
+ float ISP_CoarsePorosity; // Coarse Porosity
+ float ISP_LeakageWellPermeability; // Leakage Well Permeability
+ float ISP_LongitudinalDispersivity; //Longitudinal dispersivity
+ float ISP_TransverseDispersivity; //Transverse dispersivity
+ float ISP_FineBrooksCoreyLambda;
+ float ISP_FineBrooksCoreyEntryPressure;
+ float ISP_CoarseBrooksCoreyLambda;
+ float ISP_CoarseBrooksCoreyEntryPressure;
+ float ISP_FineResidualSaturationWetting;
+ float ISP_FineResidualSaturationNonWetting;
+ float ISP_CoarseResidualSaturationWetting;
+ float ISP_CoarseResidualSaturationNonWetting;
+
+ InterfaceSoilProperties(const char* isp_filename)
+ //Initialization of ISP Parameters
+ {
+ using namespace std;
+ std::cout
+ << "-----> 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>"))
+ //cout << "-----> ISP: Soil Properties reading ... \n";
+ //ISP perameters initialization:
+ if (reader == string("<Permeability>"))
+ {
+ input >> reader;
+ ISP_Permeability = atof(reader);
+ cout << "-----> ISP: Permeability: " << ISP_Permeability
+ << "\n";
+ }
+ if (reader == string("<FinePermeability>"))
+ {
+ input >> reader;
+ ISP_FinePermeability = atof(reader);
+ cout << "-----> ISP: Fine permeability: "
+ << ISP_FinePermeability << "\n";
+ }
+ if (reader == string("<CoarsePermeability>"))
+ {
+ input >> reader;
+ ISP_CoarsePermeability = atof(reader);
+ cout << "-----> ISP: Coarse permeability: "
+ << ISP_CoarsePermeability << "\n";
+ }
+ if (reader == string("<Porosity>"))
+ {
+ input >> reader;
+ ISP_Porosity = atof(reader);
+ cout << "-----> ISP: Porosity: " << ISP_Porosity << "\n";
+ }
+ if (reader == string("<FinePorosity>"))
+ {
+ input >> reader;
+ ISP_FinePorosity = atof(reader);
+ cout << "-----> ISP: Fine porosity: " << ISP_FinePorosity
+ << "\n";
+ }
+ if (reader == string("<CoarsePorosity>"))
+ {
+ input >> reader;
+ ISP_CoarsePorosity = atof(reader);
+ cout << "-----> ISP: Coarse porosity: " << ISP_CoarsePorosity
+ << "\n";
+ }
+ if (reader == string("<LeakageWellPermeability>"))
+ {
+ input >> reader;
+ ISP_LeakageWellPermeability = atof(reader);
+ cout << "-----> ISP: Leakage Well Permeability: "
+ << ISP_LeakageWellPermeability << "\n";
+ }
+ if (reader == string("<LongitudinalDispersivity>"))
+ {
+ input >> reader;
+ ISP_LongitudinalDispersivity = atof(reader);
+ cout << "-----> ISP: Longitudinal dispersivity: "
+ << ISP_LongitudinalDispersivity << "\n";
+ }
+ if (reader == string("<TransverseDispersivity>"))
+ {
+ input >> reader;
+ ISP_TransverseDispersivity = atof(reader);
+ cout << "-----> ISP: Transverse dispersivity: "
+ << ISP_TransverseDispersivity << "\n";
+ }
+ if (reader == string("<FineBrooksCoreyLambda>"))
+ {
+ input >> reader;
+ ISP_FineBrooksCoreyLambda = atof(reader);
+ cout << "-----> ISP: Brooks-Corey lambda, fine: "
+ << ISP_FineBrooksCoreyLambda << "\n";
+ }
+ if (reader == string("<FineBrooksCoreyEntryPressure>"))
+ {
+ input >> reader;
+ ISP_FineBrooksCoreyEntryPressure = atof(reader);
+ cout << "-----> ISP: Brooks-Corey entry pressure, fine: "
+ << ISP_FineBrooksCoreyEntryPressure << "\n";
+ }
+ if (reader == string("<CoarseBrooksCoreyLambda>"))
+ {
+ input >> reader;
+ ISP_CoarseBrooksCoreyLambda = atof(reader);
+ cout << "-----> ISP: Brooks-Corey lambda, coarse: "
+ << ISP_CoarseBrooksCoreyLambda << "\n";
+ }
+ if (reader == string("<CoarseBrooksCoreyEntryPressure>"))
+ {
+ input >> reader;
+ ISP_CoarseBrooksCoreyEntryPressure = atof(reader);
+ cout << "-----> ISP: Brooks-Corey entry pressure, coarse: "
+ << ISP_CoarseBrooksCoreyEntryPressure << "\n";
+ }
+ if (reader == string("<FineResidualSaturationWetting>"))
+ {
+ input >> reader;
+ ISP_FineResidualSaturationWetting = atof(reader);
+ cout << "-----> ISP: Residual saturation wetting phase, fine: "
+ << ISP_FineResidualSaturationWetting << "\n";
+ }
+ if (reader == string("<FineResidualSaturationNonWetting>"))
+ {
+ input >> reader;
+ ISP_FineResidualSaturationNonWetting = atof(reader);
+ cout << "-----> ISP: Residual saturation nonwetting phase, fine: "
+ << ISP_FineResidualSaturationNonWetting << "\n";
+ }
+ if (reader == string("<CoarseResidualSaturationWetting>"))
+ {
+ input >> reader;
+ ISP_CoarseResidualSaturationWetting = atof(reader);
+ cout << "-----> ISP: Residual saturation wetting phase, coarse: "
+ << ISP_CoarseResidualSaturationWetting << "\n";
+ }
+ if (reader == string("<CoarseResidualSaturationNonWetting>"))
+ {
+ input >> reader;
+ ISP_CoarseResidualSaturationNonWetting = atof(reader);
+ cout << "-----> ISP: Residual saturation nonwetting phase, coarse: "
+ << ISP_CoarseResidualSaturationNonWetting << "\n";
+ }
+ }
+ input.close();
+ }
+
+};
+
+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 IFP_GasDiffCoeff;// Gas Diffusion Coefficient
+ float IFP_CO2ResidSat;// Residual Saturation of CO2
+ float IFP_MolecularDiffusionCoefficient;
+
+ InterfaceFluidProperties(const char* ifp_filename)
+ //Initialization of IFP Parameters
+ {
+ using namespace std;
+ std::cout
+ << "-----> 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("<GasDiffusionCoeff>"))
+ {
+ input >> reader;
+ IFP_GasDiffCoeff = atof(reader);
+ cout << "-----> IFP: Gas Diffusion Coefficient: "
+ << IFP_GasDiffCoeff << "\n";
+ }
+ if (reader == string("<CO2ResidualSaturation>"))
+ {
+ input >> reader;
+ IFP_CO2ResidSat = atof(reader);
+ cout << "-----> IFP: Residual Saturation of CO2: "
+ << IFP_CO2ResidSat << "\n";
+ }
+ if (reader == string("<MolecularDiffusionCoefficient>"))
+ {
+ input >> reader;
+ IFP_MolecularDiffusionCoefficient = atof(reader);
+ cout << "-----> IFP: Molecular diffusion coefficient: "
+ << IFP_MolecularDiffusionCoefficient << "\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):
+ float IPP_DepthBOR; // Depth BOR
+ float IPP_InjectionWellRate; // Injection Well Rate
+ float IPP_InjectionWindowSize; // Injection Well Window Size
+ float IPP_UpperPressure; // Pressure at a top boundary
+ float IPP_LowerPressure; // Pressure at a lower boundary
+ float IPP_InfiltrationRate; // A Infiltration rate
+ float IPP_MaxTimeStepSize; // Maximum time step size
+ float IPP_InfiltrationStartTime; // time to stop an infiltration
+ float IPP_InfiltrationEndTime; // time to stop an infiltration
+ float IPP_SimulationNumber;
+
+ InterfaceProblemProperties(const char* ipp_filename)
+ //Initialization of IPP Parameters
+ {
+ using namespace std;
+ std::cout
+ << "-----> IPP: Interface Soil Properties Initialization ...\n";
+ //IPP input file defenition
+ 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("<DepthBOR>"))
+ {
+ input >> reader;
+ IPP_DepthBOR = atof(reader);
+ cout << "-----> IPP: Depth BOR: " << IPP_DepthBOR << "\n";
+ }
+ if (reader == string("<InjectionWellRate>"))
+ {
+ input >> reader;
+ IPP_InjectionWellRate = atof(reader);
+ cout << "-----> IPP: Injection Well Rate: "
+ << IPP_InjectionWellRate << "\n";
+ }
+ if (reader == string("<InjectionWellWindowSize>"))
+ {
+ input >> reader;
+ IPP_InjectionWindowSize = atof(reader);
+ cout << "-----> IPP: Injection Well Window Size: "
+ << IPP_InjectionWindowSize << "\n";
+ }
+ if (reader == string("<UpperPressure>"))
+ {
+ input >> reader;
+ IPP_UpperPressure = atof(reader);
+ cout << "-----> IPP: Upper pressure: "
+ << IPP_UpperPressure << "\n";
+ }
+ if (reader == string("<LowerPressure>"))
+ {
+ input >> reader;
+ IPP_LowerPressure = atof(reader);
+ cout << "-----> IPP: Lower pressure: "
+ << IPP_LowerPressure << "\n";
+ }
+ if (reader == string("<InfiltrationRate>"))
+ {
+ input >> reader;
+ IPP_InfiltrationRate = atof(reader);
+ cout << "-----> IPP: Infiltration rate: "
+ << IPP_InfiltrationRate << "\n";
+ }
+ if (reader == string("<MaxTimeStepSize>"))
+ {
+ input >> reader;
+ IPP_MaxTimeStepSize = atof(reader);
+ cout << "-----> IPP: Maximum time step size: "
+ << IPP_MaxTimeStepSize << "\n";
+ }
+ if (reader == string("<InfiltrationStartTime>"))
+ {
+ input >> reader;
+ IPP_InfiltrationStartTime = atof(reader);
+ cout << "-----> IPP: Start time of infiltration: "
+ << IPP_InfiltrationStartTime << "\n";
+ }
+ if (reader == string("<InfiltrationEndTime>"))
+ {
+ input >> reader;
+ IPP_InfiltrationEndTime = atof(reader);
+ cout << "-----> IPP: End time of infiltration: "
+ << IPP_InfiltrationEndTime << "\n";
+ }
+ if (reader == string("<SimulationNumber>"))
+ {
+ input >> reader;
+ IPP_SimulationNumber = atof(reader);
+ cout << "-----> IPP: Output Name: "
+ << IPP_SimulationNumber << "\n";
+ }
+ }
+ input.close();
+ }
+
+};
+
+} // end namespace
+#endif
diff --git a/appl/lecture/msm/1p2cvs2p/interface1p2c.xml b/appl/lecture/msm/1p2cvs2p/interface1p2c.xml
new file mode 100755
index 0000000000..af962e5a3f
--- /dev/null
+++ b/appl/lecture/msm/1p2cvs2p/interface1p2c.xml
@@ -0,0 +1,54 @@
+<iPCM_problem>
+
+ <SoilProperties>
+ <FinePermeability>
+ 3.1e-12
+ </FinePermeability>
+ <CoarsePermeability>
+ 4.6e-11
+ </CoarsePermeability>
+ <FinePorosity>
+ 0.35
+ </FinePorosity>
+ <CoarsePorosity>
+ 0.39
+ </CoarsePorosity>
+ <LongitudinalDispersivity>
+ 0.2
+ </LongitudinalDispersivity>
+ <TransverseDispersivity>
+ 0.02
+ </TransverseDispersivity>
+ </SoilProperties>
+
+ <FluidProperties>
+ <MolecularDiffusionCoefficient>
+ 1e-9
+ </MolecularDiffusionCoefficient>
+ </FluidProperties>
+
+ <BoundaryAndInitialConditions>
+ <UpperPressure>
+ 1.6e5
+ </UpperPressure>
+ <LowerPressure>
+ 1.0e5
+ </LowerPressure>
+ <InfiltrationRate>
+ 0.1
+ </InfiltrationRate>
+ <MaxTimeStepSize>
+ 500
+ </MaxTimeStepSize>
+ <InfiltrationStartTime>
+ 0
+ </InfiltrationStartTime>
+ <InfiltrationEndTime>
+ 1000
+ </InfiltrationEndTime>
+ <SimulationNumber>
+ 1
+ </SimulationNumber>
+ </BoundaryAndInitialConditions>
+
+</iPCM_problem>
diff --git a/appl/lecture/msm/1p2cvs2p/interface2p.xml b/appl/lecture/msm/1p2cvs2p/interface2p.xml
new file mode 100755
index 0000000000..480e4e4d8d
--- /dev/null
+++ b/appl/lecture/msm/1p2cvs2p/interface2p.xml
@@ -0,0 +1,63 @@
+<iPCM_problem>
+
+ <SoilProperties>
+ <FinePermeability>
+ 9e-12
+ </FinePermeability>
+ <CoarsePermeability>
+ 4.6e-10
+ </CoarsePermeability>
+ <FinePorosity>
+ 0.39
+ </FinePorosity>
+ <CoarsePorosity>
+ 0.43
+ </CoarsePorosity>
+ <FineBrooksCoreyLambda>
+ 3.5
+ </FineBrooksCoreyLambda>
+ <FineBrooksCoreyEntryPressure>
+ 1800
+ </FineBrooksCoreyEntryPressure>
+ <CoarseBrooksCoreyLambda>
+ 2.0
+ </CoarseBrooksCoreyLambda>
+ <CoarseBrooksCoreyEntryPressure>
+ 200
+ </CoarseBrooksCoreyEntryPressure>
+ <FineResidualSaturationWetting>
+ 0.18
+ </FineResidualSaturationWetting>
+ <FineResidualSaturationNonWetting>
+ 0.0
+ </FineResidualSaturationNonWetting>
+ <CoarseResidualSaturationWetting>
+ 0.05
+ </CoarseResidualSaturationWetting>
+ <CoarseResidualSaturationNonWetting>
+ 0.0
+ </CoarseResidualSaturationNonWetting>
+ </SoilProperties>
+
+ <BoundaryAndInitialConditions>
+ <LowerPressure>
+ 1.0e5
+ </LowerPressure>
+ <UpperPressure>
+ 1.1e5
+ </UpperPressure>
+ <InfiltrationRate>
+ 0.04
+ </InfiltrationRate>
+ <InfiltrationStartTime>
+ 100
+ </InfiltrationStartTime>
+ <InfiltrationEndTime>
+ 10000
+ </InfiltrationEndTime>
+ <SimulationNumber>
+ 1
+ </SimulationNumber>
+ </BoundaryAndInitialConditions>
+
+</iPCM_problem>
diff --git a/appl/lecture/msm/1p2cvs2p/lens_1p2c.cc b/appl/lecture/msm/1p2cvs2p/lens_1p2c.cc
new file mode 100755
index 0000000000..f9543a6afb
--- /dev/null
+++ b/appl/lecture/msm/1p2cvs2p/lens_1p2c.cc
@@ -0,0 +1,196 @@
+// $Id: test_1p2c.cc 2425 2009-08-27 16:53:36Z bernd $
+/*****************************************************************************
+ * Copyright (C) 2009 by Karin Erbertseder, Klaus Mosthaf *
+ * Copyright (C) 2009 by Andreas Lauser *
+ * Copyright (C) 2008 by Bernd Flemisch *
+ * 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, as long as this copyright notice *
+ * is included in its original form. *
+ * *
+ * This program is distributed WITHOUT ANY WARRANTY. *
+ *****************************************************************************/
+#include "config.h"
+
+#include "lensproblem1p2c.hh"
+
+#include <dune/common/exceptions.hh>
+#include <dune/grid/common/gridinfo.hh>
+
+#include <dune/common/mpihelper.hh>
+
+#include <iostream>
+#include <boost/format.hpp>
+
+////////////////////////
+// helper class for grid instantiation
+////////////////////////
+template <class Grid, class Scalar>
+class CreateGrid
+{
+};
+
+#if HAVE_UG
+template <class Scalar>
+class CreateGrid<Dune::UGGrid<2>, Scalar>
+{
+public:
+ static Dune::UGGrid<2> *create(const Dune::FieldVector<Scalar, 2> &upperRight,
+ const Dune::FieldVector<int, 2> &cellRes)
+ {
+ Dune::UGGrid<2> *grid = new Dune::UGGrid<2>;
+ Dune::GridFactory<Dune::UGGrid<2> > factory(grid);
+ for (int i=0; i<=cellRes[0]; i++) {
+ for (int j=0; j<=cellRes[1]; j++) {
+ Dune::FieldVector<double,2> pos;
+ pos[0] = upperRight[0]*double(i)/cellRes[0];
+ pos[1] = upperRight[1]*double(j)/cellRes[1];
+ factory.insertVertex(pos);
+ }
+ }
+
+ for (int i=0; i<cellRes[0]; i++) {
+ for (int j=0; j<cellRes[1]; j++) {
+ std::vector<unsigned int> v(4);
+ v[0] = i*(cellRes[1]+1) + j;
+ v[1] = i*(cellRes[1]+1) + j+1;
+ v[2] = (i+1)*(cellRes[1]+1) + j;
+ v[3] = (i+1)*(cellRes[1]+1) + j+1;
+ factory.insertElement(Dune::GeometryType(Dune::GeometryType::cube,2), v);
+ }
+ }
+
+ factory.createGrid();
+ grid->loadBalance();
+ return grid;
+ }
+};
+#endif
+
+template <class Scalar>
+class CreateGrid<Dune::YaspGrid<2>, Scalar>
+{
+public:
+ static Dune::YaspGrid<2> *create(const Dune::FieldVector<Scalar, 2> &upperRight,
+ const Dune::FieldVector<int, 2> &cellRes)
+ {
+ return new Dune::YaspGrid<2>(
+#ifdef HAVE_MPI
+ Dune::MPIHelper::getCommunicator(),
+#endif
+ upperRight, // upper right
+ cellRes, // number of cells
+ Dune::FieldVector<bool,2>(false), // periodic
+ 4); // overlap
+ };
+};
+
+template <class Scalar>
+class CreateGrid<Dune::SGrid<2, 2>, Scalar>
+{
+public:
+ static Dune::SGrid<2, 2> *create(const Dune::FieldVector<Scalar, 2> &upperRight,
+ const Dune::FieldVector<int, 2> &cellRes)
+ {
+ return new Dune::SGrid<2,2>(cellRes, // number of cells
+ Dune::FieldVector<Scalar, 2>(0.0), // lower left
+ upperRight); // upper right
+
+ };
+};
+
+void usage(const char *progname)
+{
+ std::cout << boost::format("usage: %s tEnd dt\n")%progname;
+ exit(1);
+};
+
+int main(int argc, char** argv)
+{
+ try {
+ typedef TTAG(LensProblem) TypeTag;
+ typedef GET_PROP_TYPE(TypeTag, PTAG(Scalar)) Scalar;
+ typedef GET_PROP_TYPE(TypeTag, PTAG(Grid)) Grid;
+ typedef GET_PROP_TYPE(TypeTag, PTAG(Problem)) Problem;
+ typedef GET_PROP_TYPE(TypeTag, PTAG(TimeManager)) TimeManager;
+ typedef Dune::FieldVector<Scalar, Grid::dimensionworld> GlobalPosition;
+
+ static const int dim = Grid::dimension;
+
+ // initialize MPI, finalize is done automatically on exit
+ Dune::MPIHelper::instance(argc, argv);
+
+ ////////////////////////////////////////////////////////////
+ // parse the command line arguments
+ ////////////////////////////////////////////////////////////
+ if (argc < 3)
+ usage(argv[0]);
+
+ // deal with the restart stuff
+ int argPos = 1;
+ bool restart = false;
+ double restartTime = 0;
+ if (std::string("--restart") == argv[argPos]) {
+ restart = true;
+ ++argPos;
+
+ std::istringstream(argv[argPos++]) >> restartTime;
+ }
+
+ if (argc - argPos != 2) {
+ usage(argv[0]);
+ }
+
+ // read the initial time step and the end time
+ double tEnd, dt;
+ std::istringstream(argv[argPos++]) >> tEnd;
+ std::istringstream(argv[argPos++]) >> dt;
+
+ ////////////////////////////////////////////////////////////
+ // create the grid
+ ////////////////////////////////////////////////////////////
+ GlobalPosition lowerLeft(0.0);
+ GlobalPosition upperRight;
+ Dune::FieldVector<int,dim> res; // cell resolution
+ upperRight[0] = 5.0;
+ upperRight[1] = 4.0;
+ res[0] = 40;
+ res[1] = 32;
+
+ std::auto_ptr<Grid> grid(CreateGrid<Grid, Scalar>::create(upperRight, res));
+
+ ////////////////////////////////////////////////////////////
+ // instantiate and run the concrete problem
+ ////////////////////////////////////////////////////////////
+
+ // specify dimensions of the low-permeable lens
+ GlobalPosition lowerLeftLens, upperRightLens;
+ lowerLeftLens[0] = 1.0;
+ lowerLeftLens[1] = 2.0;
+ upperRightLens[0] = 4.0;
+ upperRightLens[1] = 3.0;
+
+ // instantiate and run the concrete problem
+ TimeManager timeManager;
+ Problem problem(timeManager, grid->leafView(), lowerLeftLens, upperRightLens);
+ timeManager.init(problem, 0, dt, tEnd, !restart);
+ if (restart)
+ problem.restart(restartTime);
+ timeManager.run();
+ 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/msm/1p2cvs2p/lens_1p2cnewtoncontroller.hh b/appl/lecture/msm/1p2cvs2p/lens_1p2cnewtoncontroller.hh
new file mode 100644
index 0000000000..ba4e850062
--- /dev/null
+++ b/appl/lecture/msm/1p2cvs2p/lens_1p2cnewtoncontroller.hh
@@ -0,0 +1,120 @@
+// $Id$
+/*****************************************************************************
+ * Copyright (C) 2008 by Bernd Flemisch, 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, as long as this copyright notice *
+ * is included in its original form. *
+ * *
+ * This program is distributed WITHOUT ANY WARRANTY. *
+ *****************************************************************************/
+/*!
+ * \file
+ * \brief A 1p2c specific controller for the newton solver.
+ *
+ * This controller 'knows' what a 'physically meaningful' solution is
+ * which allows the newton method to abort quicker if the solution is
+ * way out of bounds.
+ */
+#ifndef DUMUX_LENS_1P2C_NEWTON_CONTROLLER_HH
+#define DUMUX_LENS_1P2C_NEWTON_CONTROLLER_HH
+
+#include <dumux/nonlinear/newtoncontroller.hh>
+
+namespace Dumux {
+/*!
+ * \ingroup TwoPTwoCBoxModel
+ * \brief A 2p2c specific controller for the newton solver.
+ *
+ * This controller 'knows' what a 'physically meaningful' solution is
+ * which allows the newton method to abort quicker if the solution is
+ * way out of bounds.
+ */
+template <class TypeTag>
+class LensOnePTwoCNewtonController
+ : public NewtonController<TypeTag >
+{
+ typedef LensOnePTwoCNewtonController<TypeTag> ThisType;
+ typedef NewtonController<TypeTag> ParentType;
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(NewtonController)) Implementation;
+
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(Scalar)) Scalar;
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(Model)) Model;
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(NewtonMethod)) NewtonMethod;
+
+ typedef typename GET_PROP(TypeTag, PTAG(SolutionTypes)) SolutionTypes;
+ typedef typename SolutionTypes::SolutionFunction SolutionFunction;
+
+
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(OnePTwoCIndices)) Indices;
+
+ enum {
+ konti = Indices::konti,
+ transport = Indices::transport
+ };
+
+public:
+ LensOnePTwoCNewtonController()
+ {
+ this->setRelTolerance(1e-7);
+ this->setTargetSteps(9);
+ this->setMaxSteps(18);
+
+ relDefect_ = 1e100;
+ //load interface-file
+
+ Dumux::InterfaceProblemProperties interfaceProbProps("interface1p2c.xml");
+
+ maxTimeStepSize_ = interfaceProbProps.IPP_MaxTimeStepSize;
+ };
+
+ //! Suggest a new time stepsize based either on the number of newton
+ //! iterations required or on the variable switch
+ void newtonEndStep(SolutionFunction &u, SolutionFunction &uOld)
+ {
+ // call the method of the base class
+ ParentType::newtonEndStep(u, uOld);
+
+ relDefect_ = 0;
+ for (unsigned i = 0; i < (*u).size(); ++i) {
+ Scalar normP = std::max(1e3,std::abs((*u)[i][konti]));
+ normP = std::max(normP, std::abs((*uOld)[i][konti]));
+
+ Scalar normTrans = std::max(1e-3,std::abs((*u)[i][transport]));
+ normTrans = std::max(normTrans, std::abs((*uOld)[i][transport]));
+
+ relDefect_ = std::max(relDefect_,
+ std::abs((*u)[i][konti] - (*uOld)[i][konti])/normP);
+ relDefect_ = std::max(relDefect_,
+ std::abs((*u)[i][transport] - (*uOld)[i][transport])/normTrans);
+ }
+ }
+
+ //! Suggest a new time stepsize based either on the number of newton
+ //! iterations required or on the variable switch
+ Scalar suggestTimeStepSize(Scalar oldTimeStep) const
+ {
+ // use function of the newtoncontroller
+ return std::min(maxTimeStepSize_, ParentType::suggestTimeStepSize(oldTimeStep));
+ }
+
+ //! Returns true iff the current solution can be considered to
+ //! be acurate enough
+ bool newtonConverged()
+ {
+ return relDefect_ <= this->tolerance_;
+ };
+
+protected:
+ Scalar relDefect_;
+private:
+ Scalar maxTimeStepSize_;
+};
+}
+
+#endif
diff --git a/appl/lecture/msm/1p2cvs2p/lens_2p.cc b/appl/lecture/msm/1p2cvs2p/lens_2p.cc
new file mode 100755
index 0000000000..ac8eff58b6
--- /dev/null
+++ b/appl/lecture/msm/1p2cvs2p/lens_2p.cc
@@ -0,0 +1,196 @@
+// $Id: test_2p.cc 2539 2009-10-13 14:49:40Z bernd $
+/*****************************************************************************
+ * Copyright (C) 2010 by Klaus Mosthaf *
+ * Copyright (C) 2007-2008 by Bernd Flemisch *
+ * Copyright (C) 2008-2009 by Andreas Lauser *
+ * Institute of Hydraulic Engineering *
+ * University of Stuttgart, Germany *
+ * email: <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, as long as this copyright notice *
+ * is included in its original form. *
+ * *
+ * This program is distributed WITHOUT ANY WARRANTY. *
+ *****************************************************************************/
+#include "config.h"
+
+#include "lensproblem2p.hh"
+
+#include <dune/common/exceptions.hh>
+#include <dune/grid/common/gridinfo.hh>
+
+#include <dune/common/mpihelper.hh>
+
+#include <iostream>
+#include <boost/format.hpp>
+
+////////////////////////
+// helper class for grid instantiation
+////////////////////////
+template <class Grid, class Scalar>
+class CreateGrid
+{
+};
+
+#if HAVE_UG
+template <class Scalar>
+class CreateGrid<Dune::UGGrid<2>, Scalar>
+{
+public:
+ static Dune::UGGrid<2> *create(const Dune::FieldVector<Scalar, 2> &upperRight,
+ const Dune::FieldVector<int, 2> &cellRes)
+ {
+ Dune::UGGrid<2> *grid = new Dune::UGGrid<2>;
+ Dune::GridFactory<Dune::UGGrid<2> > factory(grid);
+ for (int i=0; i<=cellRes[0]; i++) {
+ for (int j=0; j<=cellRes[1]; j++) {
+ Dune::FieldVector<double,2> pos;
+ pos[0] = upperRight[0]*double(i)/cellRes[0];
+ pos[1] = upperRight[1]*double(j)/cellRes[1];
+ factory.insertVertex(pos);
+ }
+ }
+
+ for (int i=0; i<cellRes[0]; i++) {
+ for (int j=0; j<cellRes[1]; j++) {
+ std::vector<unsigned int> v(4);
+ v[0] = i*(cellRes[1]+1) + j;
+ v[1] = i*(cellRes[1]+1) + j+1;
+ v[2] = (i+1)*(cellRes[1]+1) + j;
+ v[3] = (i+1)*(cellRes[1]+1) + j+1;
+ factory.insertElement(Dune::GeometryType(Dune::GeometryType::cube,2), v);
+ }
+ }
+
+ factory.createGrid();
+ grid->loadBalance();
+ return grid;
+ }
+};
+#endif
+
+template <class Scalar>
+class CreateGrid<Dune::YaspGrid<2>, Scalar>
+{
+public:
+ static Dune::YaspGrid<2> *create(const Dune::FieldVector<Scalar, 2> &upperRight,
+ const Dune::FieldVector<int, 2> &cellRes)
+ {
+ return new Dune::YaspGrid<2>(
+#ifdef HAVE_MPI
+ Dune::MPIHelper::getCommunicator(),
+#endif
+ upperRight, // upper right
+ cellRes, // number of cells
+ Dune::FieldVector<bool,2>(false), // periodic
+ 4); // overlap
+ };
+};
+
+template <class Scalar>
+class CreateGrid<Dune::SGrid<2, 2>, Scalar>
+{
+public:
+ static Dune::SGrid<2, 2> *create(const Dune::FieldVector<Scalar, 2> &upperRight,
+ const Dune::FieldVector<int, 2> &cellRes)
+ {
+ return new Dune::SGrid<2,2>(cellRes, // number of cells
+ Dune::FieldVector<Scalar, 2>(0.0), // lower left
+ upperRight); // upper right
+
+ };
+};
+
+void usage(const char *progname)
+{
+ std::cout << boost::format("usage: %s [--restart restartTime] tEnd dt\n")%progname;
+ exit(1);
+};
+
+int main(int argc, char** argv)
+{
+ try {
+ typedef TTAG(LensProblem) TypeTag;
+ typedef GET_PROP_TYPE(TypeTag, PTAG(Scalar)) Scalar;
+ typedef GET_PROP_TYPE(TypeTag, PTAG(Grid)) Grid;
+ typedef GET_PROP_TYPE(TypeTag, PTAG(Problem)) Problem;
+ typedef GET_PROP_TYPE(TypeTag, PTAG(TimeManager)) TimeManager;
+ typedef Dune::FieldVector<Scalar, Grid::dimensionworld> GlobalPosition;
+
+ static const int dim = Grid::dimension;
+
+ // initialize MPI, finalize is done automatically on exit
+ Dune::MPIHelper::instance(argc, argv);
+
+ ////////////////////////////////////////////////////////////
+ // parse the command line arguments
+ ////////////////////////////////////////////////////////////
+ if (argc < 3)
+ usage(argv[0]);
+
+ // deal with the restart stuff
+ int argPos = 1;
+ bool restart = false;
+ double restartTime = 0;
+ if (std::string("--restart") == argv[argPos]) {
+ restart = true;
+ ++argPos;
+
+ std::istringstream(argv[argPos++]) >> restartTime;
+ }
+
+ if (argc - argPos != 2) {
+ usage(argv[0]);
+ }
+
+ // read the initial time step and the end time
+ double tEnd, dt;
+ std::istringstream(argv[argPos++]) >> tEnd;
+ std::istringstream(argv[argPos++]) >> dt;
+
+ ////////////////////////////////////////////////////////////
+ // create the grid
+ ////////////////////////////////////////////////////////////
+ GlobalPosition lowerLeft(0.0);
+ GlobalPosition upperRight;
+ Dune::FieldVector<int,dim> res; // cell resolution
+ upperRight[0] = 5.0;
+ upperRight[1] = 4.0;
+ res[0] = 40;
+ res[1] = 32;
+
+ std::auto_ptr<Grid> grid(CreateGrid<Grid, Scalar>::create(upperRight, res));
+
+ ////////////////////////////////////////////////////////////
+ // instantiate and run the concrete problem
+ ////////////////////////////////////////////////////////////
+
+ // specify dimensions of the low-permeable lens
+ GlobalPosition lowerLeftLens, upperRightLens;
+ lowerLeftLens[0] = 1.0;
+ lowerLeftLens[1] = 2.0;
+ upperRightLens[0] = 4.0;
+ upperRightLens[1] = 3.0;
+
+ // instantiate and run the concrete problem
+ TimeManager timeManager;
+ Problem problem(timeManager, grid->leafView(), lowerLeftLens, upperRightLens);
+ timeManager.init(problem, 0, dt, tEnd, !restart);
+ if (restart)
+ problem.restart(restartTime);
+ timeManager.run();
+ 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/msm/1p2cvs2p/lens_2pnewtoncontroller.hh b/appl/lecture/msm/1p2cvs2p/lens_2pnewtoncontroller.hh
new file mode 100644
index 0000000000..5d6a5fa8dd
--- /dev/null
+++ b/appl/lecture/msm/1p2cvs2p/lens_2pnewtoncontroller.hh
@@ -0,0 +1,86 @@
+// $Id$
+/*****************************************************************************
+ * Copyright (C) 2008 by Bernd Flemisch, 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, as long as this copyright notice *
+ * is included in its original form. *
+ * *
+ * This program is distributed WITHOUT ANY WARRANTY. *
+ *****************************************************************************/
+/*!
+ * \file
+ * \brief A 1p2c specific controller for the newton solver.
+ *
+ * This controller 'knows' what a 'physically meaningful' solution is
+ * which allows the newton method to abort quicker if the solution is
+ * way out of bounds.
+ */
+#ifndef DUMUX_LENS_2P_NEWTON_CONTROLLER_HH
+#define DUMUX_LENS_2P_NEWTON_CONTROLLER_HH
+
+#include <dumux/nonlinear/newtoncontroller.hh>
+
+namespace Dumux {
+/*!
+ * \ingroup TwoPTwoCBoxModel
+ * \brief A 2p2c specific controller for the newton solver.
+ *
+ * This controller 'knows' what a 'physically meaningful' solution is
+ * which allows the newton method to abort quicker if the solution is
+ * way out of bounds.
+ */
+template <class TypeTag>
+class LensTwoPNewtonController : public NewtonController<TypeTag>
+{
+ typedef LensTwoPNewtonController<TypeTag> ThisType;
+ typedef NewtonController<TypeTag> ParentType;
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(NewtonController)) Implementation;
+
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(Scalar)) Scalar;
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(Model)) Model;
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(NewtonMethod)) NewtonMethod;
+
+ typedef typename GET_PROP(TypeTag, PTAG(SolutionTypes)) SolutionTypes;
+ typedef typename SolutionTypes::SolutionFunction SolutionFunction;
+
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(TwoPIndices)) Indices;
+
+ enum {
+ pressureIdx = Indices::pressureIdx,
+ saturationIdx = Indices::saturationIdx,
+ pW = Indices::pW,
+ sN = Indices::sN,
+ };
+
+public:
+ LensTwoPNewtonController()
+ {
+ this->setRelTolerance(1e-7);
+ this->setTargetSteps(9);
+ this->setMaxSteps(18);
+
+ Dumux::InterfaceProblemProperties interfaceProbProps("interface2p.xml");
+
+ maxTimeStepSize_ = interfaceProbProps.IPP_MaxTimeStepSize;
+ };
+
+ //! Suggest a new time stepsize based either on the number of newton
+ //! iterations required or on the variable switch
+ Scalar suggestTimeStepSize(Scalar oldTimeStep) const
+ {
+ // use function of the newtoncontroller
+ return std::min(maxTimeStepSize_, ParentType::suggestTimeStepSize(oldTimeStep));
+ }
+
+private:
+ Scalar maxTimeStepSize_;
+};
+}
+
+#endif
diff --git a/appl/lecture/msm/1p2cvs2p/lensproblem1p2c.hh b/appl/lecture/msm/1p2cvs2p/lensproblem1p2c.hh
new file mode 100644
index 0000000000..d1956646db
--- /dev/null
+++ b/appl/lecture/msm/1p2cvs2p/lensproblem1p2c.hh
@@ -0,0 +1,384 @@
+// $Id: lensproblem.hh 2539 2009-10-13 14:49:40Z bernd $
+/*****************************************************************************
+ * Copyright (C) 2007-2008 by Klaus Mosthaf *
+ * Copyright (C) 2007-2008 by Bernd Flemisch *
+ * Copyright (C) 2008-2009 by Andreas Lauser *
+ * Institute of Hydraulic Engineering *
+ * University of Stuttgart, Germany *
+ * email: <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, as long as this copyright notice *
+ * is included in its original form. *
+ * *
+ * This program is distributed WITHOUT ANY WARRANTY. *
+ *****************************************************************************/
+#ifndef DUMUX_LENSPROBLEM1P2C_HH
+#define DUMUX_LENSPROBLEM1P2C_HH
+
+#if HAVE_UG
+#include <dune/grid/uggrid.hh>
+#endif
+
+#include <dune/grid/yaspgrid.hh>
+#include <dune/grid/sgrid.hh>
+#include "external_interface.hh"
+
+#include <dumux/material/components/simpleh2o.hh>
+#include <dumux/material/components/simplednapl.hh>
+#include <dumux/material/fluidsystems/liquidphase.hh>
+
+#include <dumux/boxmodels/1p2c/1p2cmodel.hh>
+
+#include <water_contaminant.hh>
+
+#include "lensspatialparameters1p2c.hh"
+
+namespace Dumux
+{
+
+ template <class TypeTag>
+ class LensProblem;
+
+ //////////
+ // Specify the properties for the lens problem
+ //////////
+ namespace Properties
+ {
+ NEW_TYPE_TAG(LensProblem, INHERITS_FROM(BoxOnePTwoC));
+
+ // Set the grid type
+ SET_PROP(LensProblem, Grid)
+ {
+#if HAVE_UG
+ typedef Dune::UGGrid<2> type;
+#else
+ typedef Dune::YaspGrid<2> type;
+ //typedef Dune::SGrid<2, 2> type;
+#endif
+ };
+
+ SET_PROP(LensProblem, 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
+ };
+
+ // Set the problem property
+ SET_PROP(LensProblem, Problem)
+ {
+ typedef Dumux::LensProblem<TypeTag> type;
+ };
+
+ // Set fluid configuration
+ SET_PROP(LensProblem, FluidSystem)
+ {
+ typedef Dumux::WaterContaminant<TypeTag> type;
+ };
+
+ // Set the spatial parameters
+ SET_PROP(LensProblem, SpatialParameters)
+ {
+ typedef Dumux::LensSpatialParameters1p2c<TypeTag> type;
+ };
+
+ // Enable gravity
+ SET_BOOL_PROP(LensProblem, EnableGravity, true);
+ }
+
+ /*!
+ * \ingroup TwoPBoxProblems
+ * \brief Soil decontamination problem where DNAPL infiltrates a fully
+ * water saturated medium.
+ *
+ * The domain is sized 6m times 4m and features a rectangular lens
+ * with low permeablility which spans from (1 m , 2 m) to (4 m, 3 m)
+ * and is surrounded by a medium with higher permability.
+ *
+ * On the top and the bottom of the domain neumann boundary conditions
+ * are used, while dirichlet conditions apply on the left and right
+ * boundaries.
+ *
+ * DNAPL is injected at the top boundary from 3m to 4m at a rate of
+ * 0.04 kg/(s m), the remaining neumann boundaries are no-flow
+ * boundaries.
+ *
+ * The dirichlet boundaries on the left boundary is the hydrostatic
+ * pressure scaled by a factor of 1.125, while on the right side it is
+ * just the hydrostatic pressure. The DNAPL saturation on both sides
+ * is zero.
+ *
+ * This problem uses the \ref TwoPBoxModel.
+ *
+ * This problem should typically simulated until \f$t_{\text{end}} =
+ * 50\,000\;s\f$ is reached. A good choice for the initial time step size
+ * is \f$t_{\text{inital}} = 1\,000\;s\f$.
+ *
+ * To run the simulation execute the following line in shell:
+ * <tt>./lens_1p2c 50000 100</tt>
+ */
+ template <class TypeTag >
+ class LensProblem : public OnePTwoCBoxProblem<TypeTag>
+ {
+ typedef LensProblem<TypeTag> ThisType;
+ typedef OnePTwoCBoxProblem<TypeTag> ParentType;
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(GridView)) GridView;
+
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(OnePTwoCIndices)) Indices;
+ enum
+ {
+ // Grid and world dimension
+ dim = GridView::dimension,
+ dimWorld = GridView::dimensionworld,
+
+ // indices of the primary variables
+ contiEqIdx = Indices::contiEqIdx,
+ transEqIdx = Indices::transEqIdx
+ };
+
+ 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(Scalar)) Scalar;
+ typedef Dune::FieldVector<Scalar, dim> LocalPosition;
+ typedef Dune::FieldVector<Scalar, dimWorld> GlobalPosition;
+
+ public:
+ LensProblem(TimeManager &timeManager,
+ const GridView &gridView,
+ const GlobalPosition &lensLowerLeft,
+ const GlobalPosition &lensUpperRight)
+ : ParentType(timeManager, gridView)
+ {
+ this->spatialParameters().setLensCoords(lensLowerLeft, lensUpperRight);
+
+ bboxMin_ = 0.0;
+ bboxMax_[0] = 5.0;
+ bboxMax_[1] = 4.0;
+
+ //load interface-file
+ Dumux::InterfaceProblemProperties interfaceProbProps("interface1p2c.xml");
+
+ upperPressure_ = interfaceProbProps.IPP_UpperPressure;
+ lowerPressure_ = interfaceProbProps.IPP_LowerPressure;
+ infiltrationRate_ = interfaceProbProps.IPP_InfiltrationRate;
+ infiltrationStartTime_= interfaceProbProps.IPP_InfiltrationStartTime;
+ infiltrationEndTime_= interfaceProbProps.IPP_InfiltrationEndTime;
+ }
+
+ /*!
+ * \name Problem parameters
+ */
+ // \{
+
+ /*!
+ * \brief The problem name.
+ *
+ * This is used as a prefix for files generated by the simulation.
+ */
+ const char *name() const
+ {
+ std::string simName = "lens-1p2c_run";
+ Dumux::InterfaceProblemProperties interfaceProbProps("interface1p2c.xml");
+ Scalar simNum = interfaceProbProps.IPP_SimulationNumber;
+
+ return (str(boost::format("%s-%02d")
+ %simName%simNum).c_str());
+ }
+
+ /*!
+ * \brief Returns the temperature within the domain.
+ *
+ * This problem assumes a temperature of 10 degrees Celsius.
+ */
+ Scalar temperature(const Element &element,
+ const FVElementGeometry &fvElemGeom,
+ int scvIdx) const
+ {
+ return 273.15 + 10; // -> 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 on the boundary for which the
+ * conditions needs to be specified
+ */
+ void boundaryTypes(BoundaryTypes &values, const Vertex &vertex) const
+ {
+ const GlobalPosition globalPos = vertex.geometry().center();
+
+ if (onUpperBoundary_(globalPos) || onLowerBoundary_(globalPos))
+ values.setAllDirichlet();
+ else
+ values.setAllNeumann();
+ if (onInlet_(globalPos))
+ values.setNeumann(transEqIdx);
+ }
+
+ /*!
+ * \brief Evaluate the boundary conditions for a dirichlet
+ * control volume.
+ *
+ * \param values The dirichlet values for the primary variables
+ * \param vertex The vertex representing the "half volume on the boundary"
+ *
+ * For this method, the \a values parameter stores primary variables.
+ */
+ void dirichlet(PrimaryVariables &values, const Vertex &vertex) const
+ {
+ const GlobalPosition globalPos = vertex.geometry().center();
+
+ if (onUpperBoundary_(globalPos))
+ {
+ values[contiEqIdx] = upperPressure_;
+ values[transEqIdx] = 0.0;
+ }
+ else if (onLowerBoundary_(globalPos))
+ {
+ values[contiEqIdx] = lowerPressure_;
+ values[transEqIdx] = 0.0;
+ }
+ else
+ values = 0.0;
+ }
+
+ /*!
+ * \brief Evaluate the boundary conditions for a neumann
+ * boundary segment.
+ *
+ * For this method, the \a values parameter stores the mass flux
+ * in normal direction of each phase. Negative values mean influx.
+ */
+ void neumann(PrimaryVariables &values,
+ const Element &element,
+ const FVElementGeometry &fvElemGeom,
+ const Intersection &isIt,
+ int scvIdx,
+ int boundaryFaceIdx) const
+ {
+ const GlobalPosition &globalPos
+ = element.geometry().corner(scvIdx);
+
+ values = 0.0;
+
+ const Scalar& time = this->timeManager().time();
+
+ if (time >= infiltrationStartTime_ && time <= infiltrationEndTime_)
+ {
+ if (onInlet_(globalPos))
+ values[transEqIdx] = -infiltrationRate_; //
+ }
+ }
+ // \}
+
+ /*!
+ * \name Volume terms
+ */
+ // \{
+
+ /*!
+ * \brief Evaluate the source term for all phases within a given
+ * sub-control-volume.
+ *
+ * For this method, the \a values parameter stores the rate mass
+ * generated or annihilate per volume unit. Positive values mean
+ * that mass is created, negative ones mean that it vanishes.
+ */
+ void source(PrimaryVariables &values,
+ const Element &element,
+ const FVElementGeometry &,
+ int subControlVolumeIdx) const
+ {
+ values = Scalar(0.0);
+ }
+
+ /*!
+ * \brief Evaluate the initial value for a control volume.
+ *
+ * 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);
+
+ // no contaminant, hydrostatic pressure
+ const Scalar depth = this->bboxMax()[1] - globalPos[1];
+ const Scalar height = this->bboxMax()[1] - this->bboxMin()[1];
+
+ values[contiEqIdx] = upperPressure_ - depth/height*(upperPressure_-lowerPressure_);
+ values[transEqIdx] = 0.0;
+ }
+ // \}
+
+ private:
+ bool onLeftBoundary_(const GlobalPosition &globalPos) const
+ {
+ return globalPos[0] < this->bboxMin()[0] + eps_;
+ }
+
+ bool onRightBoundary_(const GlobalPosition &globalPos) const
+ {
+ return globalPos[0] > this->bboxMax()[0] - eps_;
+ }
+
+ bool onLowerBoundary_(const GlobalPosition &globalPos) const
+ {
+ return globalPos[1] < this->bboxMin()[1] + eps_;
+ }
+
+ bool onUpperBoundary_(const GlobalPosition &globalPos) const
+ {
+ return globalPos[1] > this->bboxMax()[1] - eps_;
+ }
+
+ bool onInlet_(const GlobalPosition &globalPos) const
+ {
+ Scalar width = this->bboxMax()[0] - this->bboxMin()[0];
+ Scalar lambda = (this->bboxMax()[0] - globalPos[0])/width;
+ return onUpperBoundary_(globalPos)
+ && (bboxMax_[0] - 0.35*width)/width > lambda
+ && lambda > (bboxMax_[0] - 0.55*width)/width;
+ }
+
+ static const Scalar eps_ = 3e-6;
+ GlobalPosition bboxMin_;
+ GlobalPosition bboxMax_;
+
+ Scalar upperPressure_;
+ Scalar lowerPressure_;
+ Scalar infiltrationRate_;
+ Scalar infiltrationStartTime_;
+ Scalar infiltrationEndTime_;
+ };
+} //end namespace
+
+#endif
diff --git a/appl/lecture/msm/1p2cvs2p/lensproblem2p.hh b/appl/lecture/msm/1p2cvs2p/lensproblem2p.hh
new file mode 100644
index 0000000000..febc631631
--- /dev/null
+++ b/appl/lecture/msm/1p2cvs2p/lensproblem2p.hh
@@ -0,0 +1,432 @@
+// $Id: lensproblem.hh 2539 2009-10-13 14:49:40Z bernd $
+/*****************************************************************************
+ * Copyright (C) 2007-2008 by Klaus Mosthaf *
+ * Copyright (C) 2007-2008 by Bernd Flemisch *
+ * Copyright (C) 2008-2009 by Andreas Lauser *
+ * Institute of Hydraulic Engineering *
+ * University of Stuttgart, Germany *
+ * email: <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, as long as this copyright notice *
+ * is included in its original form. *
+ * *
+ * This program is distributed WITHOUT ANY WARRANTY. *
+ *****************************************************************************/
+#ifndef DUMUX_LENSPROBLEM2P_HH
+#define DUMUX_LENSPROBLEM2P_HH
+
+#if HAVE_UG
+#include <dune/grid/uggrid.hh>
+#endif
+
+#include <dune/grid/yaspgrid.hh>
+#include <dune/grid/sgrid.hh>
+#include "external_interface.hh"
+
+#include <dumux/material/components/simpleh2o.hh>
+#include <dumux/material/components/simplednapl.hh>
+#include <dumux/material/fluidsystems/liquidphase.hh>
+
+#include <dumux/boxmodels/2p/2pmodel.hh>
+
+#include "lensspatialparameters2p.hh"
+
+namespace Dumux
+{
+
+template <class TypeTag>
+class LensProblem;
+
+//////////
+// Specify the properties for the lens problem
+//////////
+namespace Properties
+{
+NEW_TYPE_TAG(LensProblem, INHERITS_FROM(BoxTwoP));
+
+// Set the grid type
+SET_PROP(LensProblem, Grid)
+{
+#if HAVE_UG
+ typedef Dune::UGGrid<2> type;
+#else
+ typedef Dune::YaspGrid<2> type;
+ //typedef Dune::SGrid<2, 2> type;
+#endif
+};
+
+SET_PROP(LensProblem, 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
+};
+
+// Set the problem property
+SET_PROP(LensProblem, Problem)
+{
+ typedef Dumux::LensProblem<TypeTag> type;
+};
+
+// Set the wetting phase
+SET_PROP(LensProblem, WettingPhase)
+{
+private:
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(Scalar)) Scalar;
+public:
+ typedef Dumux::LiquidPhase<Scalar, Dumux::SimpleH2O<Scalar> > type;
+};
+
+// Set the non-wetting phase
+SET_PROP(LensProblem, NonwettingPhase)
+{
+private:
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(Scalar)) Scalar;
+public:
+ typedef Dumux::LiquidPhase<Scalar, Dumux::SimpleDNAPL<Scalar> > type;
+};
+
+// Set the spatial parameters
+SET_PROP(LensProblem, SpatialParameters)
+{
+ typedef Dumux::LensSpatialParameters2p<TypeTag> type;
+};
+
+// Enable gravity
+SET_BOOL_PROP(LensProblem, EnableGravity, true);
+}
+
+/*!
+ * \ingroup TwoPBoxProblems
+ * \brief Soil decontamination problem where DNAPL infiltrates a fully
+ * water saturated medium.
+ *
+ * The domain is sized 6m times 4m and features a rectangular lens
+ * with low permeablility which spans from (1 m , 2 m) to (4 m, 3 m)
+ * and is surrounded by a medium with higher permability.
+ *
+ * On the top and the bottom of the domain neumann boundary conditions
+ * are used, while dirichlet conditions apply on the left and right
+ * boundaries.
+ *
+ * DNAPL is injected at the top boundary from 3m to 4m at a rate of
+ * 0.04 kg/(s m), the remaining neumann boundaries are no-flow
+ * boundaries.
+ *
+ * The dirichlet boundaries on the left boundary is the hydrostatic
+ * pressure scaled by a factor of 1.125, while on the right side it is
+ * just the hydrostatic pressure. The DNAPL saturation on both sides
+ * is zero.
+ *
+ * This problem uses the \ref TwoPBoxModel.
+ *
+ * This problem should typically simulated until \f$t_{\text{end}} =
+ * 50\,000\;s\f$ is reached. A good choice for the initial time step size
+ * is \f$t_{\text{inital}} = 1\,000\;s\f$.
+ *
+ * To run the simulation execute the following line in shell:
+ * <tt>./lens_2p 50000 100</tt>
+ */
+template <class TypeTag >
+class LensProblem : public TwoPProblem<TypeTag>
+{
+ typedef LensProblem<TypeTag> ThisType;
+ typedef TwoPProblem<TypeTag> ParentType;
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(GridView)) GridView;
+
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(TwoPIndices)) Indices;
+
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(FluidSystem)) FluidSystem;
+ typedef TwoPFluidState<TypeTag> FluidState;
+
+ enum
+ {
+ numEq = GET_PROP_VALUE(TypeTag, PTAG(NumEq)),
+
+ // primary variable indices
+ pressureIdx = Indices::pressureIdx,
+ saturationIdx = Indices::saturationIdx,
+ pwIdx = Indices::pwIdx,
+ SnIdx = Indices::SnIdx,
+
+ // equation indices
+ contiWEqIdx = Indices::contiWEqIdx,
+ contiNEqIdx = Indices::contiNEqIdx,
+
+ // phase indices
+ wPhaseIdx = Indices::wPhaseIdx,
+ nPhaseIdx = Indices::nPhaseIdx,
+
+ // 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(Scalar)) Scalar;
+ typedef Dune::FieldVector<Scalar, dim> LocalPosition;
+ typedef Dune::FieldVector<Scalar, dimWorld> GlobalPosition;
+
+public:
+ LensProblem(TimeManager &timeManager,
+ const GridView &gridView,
+ const GlobalPosition &lensLowerLeft,
+ const GlobalPosition &lensUpperRight)
+ : ParentType(timeManager, gridView)
+ {
+ this->spatialParameters().setLensCoords(lensLowerLeft, lensUpperRight);
+ bboxMin_ = 0.0;
+ bboxMax_[0] = 5.0;
+ bboxMax_[1] = 4.0;
+
+ //load interface-file
+ Dumux::InterfaceProblemProperties interfaceProbProps("interface2p.xml");
+
+ lowerPressure_ = interfaceProbProps.IPP_LowerPressure;
+ upperPressure_ = interfaceProbProps.IPP_UpperPressure;
+ infiltrationRate_ = interfaceProbProps.IPP_InfiltrationRate;
+ infiltrationStartTime_= interfaceProbProps.IPP_InfiltrationStartTime;
+ infiltrationEndTime_= interfaceProbProps.IPP_InfiltrationEndTime;
+ }
+
+ /*!
+ * \name Problem parameters
+ */
+ // \{
+
+ /*!
+ * \brief The problem name.
+ *
+ * This is used as a prefix for files generated by the simulation.
+ */
+ const char *name() const
+ {
+ std::string simName = "lens-2p_run";
+ Dumux::InterfaceProblemProperties interfaceProbProps("interface_BL.xml");
+ Scalar simNum = interfaceProbProps.IPP_SimulationNumber;
+
+ return (str(boost::format("%s-%02d")
+ %simName%simNum).c_str());
+ }
+
+ /*!
+ * \brief Returns the temperature within the domain.
+ *
+ * This problem assumes a temperature of 10 degrees Celsius.
+ */
+ Scalar temperature(const Element &element,
+ const FVElementGeometry &fvElemGeom,
+ int scvIdx) const
+ {
+ return 273.15 + 10; // -> 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 on the boundary for which the
+ * conditions needs to be specified
+ */
+ void boundaryTypes(BoundaryTypes &values, const Vertex &vertex) const
+ {
+ const GlobalPosition globalPos = vertex.geometry().center();
+
+
+ if (onUpperBoundary_(globalPos) || onLowerBoundary_(globalPos))
+ values.setAllDirichlet();
+ else
+ values.setAllNeumann();
+
+ if (onInlet_(globalPos))
+ values.setNeumann(contiNEqIdx);
+ }
+
+ /*!
+ * \brief Evaluate the boundary conditions for a dirichlet
+ * control volume.
+ *
+ * \param values The dirichlet values for the primary variables
+ * \param vertex The vertex representing the "half volume on the boundary"
+ *
+ * For this method, the \a values parameter stores primary variables.
+ */
+ void dirichlet(PrimaryVariables &values, const Vertex &vertex) const
+ {
+ const GlobalPosition globalPos = vertex.geometry().center();
+
+ if (onUpperBoundary_(globalPos))
+ {
+ values[pwIdx] = upperPressure_;
+ values[SnIdx] = 0.0;
+ }
+ else if (onLowerBoundary_(globalPos))
+ {
+ values[pwIdx] = lowerPressure_;
+ values[SnIdx] = 0.0;
+ }
+ else
+ values = 0.0;
+
+// Scalar densityW = this->wettingPhase().density();
+//
+// if (onLeftBoundary_(globalPos))
+// {
+// Scalar depth = bboxMax_[1] - globalPos[1];
+//
+// // hydrostatic pressure scaled by alpha
+// values[pW] = lowerPressure_ - densityW*this->gravity()[1]*depth;
+// values[sN] = 0.0;
+// }
+// else if (onRightBoundary_(globalPos))
+// {
+// Scalar depth = bboxMax_[1] - globalPos[1];
+//
+// // hydrostatic pressure
+// values[pW] = lowerPressure_ - densityW*this->gravity()[1]*depth;
+// values[sN] = 0.0;
+// }
+// else
+// values = 0.0;
+ }
+
+ /*!
+ * \brief Evaluate the boundary conditions for a neumann
+ * boundary segment.
+ *
+ * For this method, the \a values parameter stores the mass flux
+ * in normal direction of each phase. Negative values mean influx.
+ */
+ void neumann(PrimaryVariables &values,
+ const Element &element,
+ const FVElementGeometry &fvElemGeom,
+ const Intersection &isIt,
+ int scvIdx,
+ int boundaryFaceIdx) const
+ {
+ const GlobalPosition &globalPos
+ = element.geometry().corner(scvIdx);
+
+ values = 0.0;
+
+ const Scalar& time = this->timeManager().time();
+
+ if (time >= infiltrationStartTime_ && time <= infiltrationEndTime_)
+ {
+ if (onInlet_(globalPos))
+ values[contiNEqIdx] = -infiltrationRate_; // kg / (m * s)
+ }
+
+ }
+ // \}
+
+ /*!
+ * \name Volume terms
+ */
+ // \{
+
+ /*!
+ * \brief Evaluate the source term for all phases within a given
+ * sub-control-volume.
+ *
+ * For this method, the \a values parameter stores the rate mass
+ * generated or annihilate per volume unit. Positive values mean
+ * that mass is created, negative ones mean that it vanishes.
+ */
+ void source(PrimaryVariables &values,
+ const Element &element,
+ const FVElementGeometry &,
+ int subControlVolumeIdx) const
+ {
+ values = Scalar(0.0);
+ }
+
+ /*!
+ * \brief Evaluate the initial value for a control volume.
+ *
+ * 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);
+
+ // no DNAPL, hydrostatic pressure
+ const Scalar depth = this->bboxMax()[1] - globalPos[1];
+ const Scalar height = this->bboxMax()[1] - this->bboxMin()[1];
+
+ values[pwIdx] = upperPressure_ - depth/height*(upperPressure_-lowerPressure_);
+ values[SnIdx] = 0.0;
+ }
+ // \}
+
+private:
+ bool onLeftBoundary_(const GlobalPosition &globalPos) const
+ {
+ return globalPos[0] < this->bboxMin()[0] + eps_;
+ }
+
+ bool onRightBoundary_(const GlobalPosition &globalPos) const
+ {
+ return globalPos[0] > this->bboxMax()[0] - eps_;
+ }
+
+ bool onLowerBoundary_(const GlobalPosition &globalPos) const
+ {
+ return globalPos[1] < this->bboxMin()[1] + eps_;
+ }
+
+ bool onUpperBoundary_(const GlobalPosition &globalPos) const
+ {
+ return globalPos[1] > this->bboxMax()[1] - eps_;
+ }
+
+ bool onInlet_(const GlobalPosition &globalPos) const
+ {
+ Scalar width = this->bboxMax()[0] - this->bboxMin()[0];
+ Scalar lambda = (this->bboxMax()[0] - globalPos[0])/width;
+ return onUpperBoundary_(globalPos) && (bboxMax_[0]-0.35*width)/width > lambda && lambda > (bboxMax_[0]-0.55*width)/width;
+ }
+
+ static const Scalar eps_ = 3e-6;
+ GlobalPosition bboxMin_;
+ GlobalPosition bboxMax_;
+
+ Scalar upperPressure_;
+ Scalar lowerPressure_;
+ Scalar infiltrationRate_;
+ Scalar infiltrationStartTime_;
+ Scalar infiltrationEndTime_;
+};
+} //end namespace
+
+#endif
diff --git a/appl/lecture/msm/1p2cvs2p/lensspatialparameters1p2c.hh b/appl/lecture/msm/1p2cvs2p/lensspatialparameters1p2c.hh
new file mode 100644
index 0000000000..0051f5f5d7
--- /dev/null
+++ b/appl/lecture/msm/1p2cvs2p/lensspatialparameters1p2c.hh
@@ -0,0 +1,184 @@
+// $Id$
+/*****************************************************************************
+ * Copyright (C) 2010 by Markus Wolff *
+ * Copyright (C) 2010 by Klaus Mosthaf *
+ * Copyright (C) 2007-2008 by Bernd Flemisch *
+ * Copyright (C) 2008-2009 by Andreas Lauser *
+ * Institute of Hydraulic Engineering *
+ * University of Stuttgart, Germany *
+ * email: <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, as long as this copyright notice *
+ * is included in its original form. *
+ * *
+ * This program is distributed WITHOUT ANY WARRANTY. *
+ *****************************************************************************/
+#ifndef DUMUX_LENSSPATIALPARAMETERS_1P2C_HH
+#define DUMUX_LENSSPATIALPARAMETERS_1P2C_HH
+
+#include <dumux/material/spatialparameters/boxspatialparameters.hh>
+#include <dumux/boxmodels/1p2c/1p2cmodel.hh>
+
+/**
+ * @file
+ * @brief Class for defining spatial parameters
+ * @author Bernd Flemisch, Klaus Mosthaf, Markus Wolff
+ */
+
+namespace Dumux
+{
+
+/** \todo Please doc me! */
+
+template<class TypeTag>
+class LensSpatialParameters1p2c : public BoxSpatialParameters<TypeTag>
+{
+ typedef BoxSpatialParameters<TypeTag> ParentType;
+ 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;
+
+// typedef typename GET_PROP_TYPE(TypeTag, PTAG(OnePTwoCIndices)) Indices;
+
+ enum {
+ dim=GridView::dimension,
+ dimWorld=GridView::dimensionworld,
+//
+// // indices of the primary variables
+// contiEqIdx = Indices::contiEqIdx,
+// transEqIdx = Indices::transEqIdx
+ };
+
+ typedef Dune::FieldVector<CoordScalar,dimWorld> GlobalPosition;
+ typedef Dune::FieldMatrix<CoordScalar,dimWorld,dimWorld> FieldMatrix;
+
+ typedef typename GridView::template Codim<0>::Entity Element;
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(FVElementGeometry)) FVElementGeometry;
+
+
+public:
+
+ LensSpatialParameters1p2c(const GridView& gridView)
+ : ParentType(gridView),
+ lensK_(0),
+ outerK_(0)
+ {
+ Dumux::InterfaceSoilProperties interfaceSoilProps("interface1p2c.xml");
+
+ lensPorosity_ = interfaceSoilProps.ISP_FinePorosity;
+ outerPorosity_ = interfaceSoilProps.ISP_CoarsePorosity;
+
+ longitudinalDispersivity_ = interfaceSoilProps.ISP_LongitudinalDispersivity;
+ transverseDispersivity_ = interfaceSoilProps.ISP_TransverseDispersivity;
+
+ for(int i = 0; i < dim; i++){
+ lensK_[i][i] = interfaceSoilProps.ISP_FinePermeability;
+ outerK_[i][i] = interfaceSoilProps.ISP_CoarsePermeability;
+ }
+ }
+
+ /*!
+ * \brief Apply the intrinsic permeability tensor to a pressure
+ * potential gradient.
+ *
+ * \param element The current finite element
+ * \param fvElemGeom The current finite volume geometry of the element
+ * \param scvIdx The index sub-control volume face where the
+ * intrinsic velocity ought to be calculated.
+ */
+ const FieldMatrix& intrinsicPermeability(const Element &element,
+ const FVElementGeometry &fvElemGeom,
+ int scvIdx) const
+ {
+ const GlobalPosition &globalPos = fvElemGeom.subContVol[scvIdx].global;
+ if (isInLens_(globalPos))
+ return lensK_;
+ return outerK_;
+ }
+
+ Scalar porosity(const Element &element,
+ const FVElementGeometry &fvElemGeom,
+ int scvIdx) const
+ {
+ const GlobalPosition &globalPos = fvElemGeom.subContVol[scvIdx].global;
+ if (isInLens_(globalPos))
+ return lensPorosity_;
+ return outerPorosity_;
+ }
+
+ //! Set the bounding box of the fine-sand lens
+ void setLensCoords(const GlobalPosition& lensLowerLeft,
+ const GlobalPosition& lensUpperRight)
+ {
+ lensLowerLeft_ = lensLowerLeft;
+ lensUpperRight_ = lensUpperRight;
+ }
+
+ /*!
+ * \brief Define the tortuosity \f$[?]\f$.
+ *
+ * \param element The finite element
+ * \param fvElemGeom The finite volume geometry
+ * \param scvIdx The local index of the sub-control volume where
+ */
+ Scalar tortuosity(const Element &element,
+ const FVElementGeometry &fvElemGeom,
+ int scvIdx) const
+ {
+ return 0;
+ }
+
+ /*!
+ * \brief Define the dispersivity \f$[m]\f$.
+ *
+ * \param element The finite element
+ * \param fvElemGeom The finite volume geometry
+ * \param scvIdx The local index of the sub-control volume where
+ */
+ Dune::FieldVector<Scalar,2> dispersivity(const Element &element,
+ const FVElementGeometry &fvElemGeom,
+ int scvIdx) const
+ {
+ Dune::FieldVector<Scalar,2> values (0);
+ values[0] = longitudinalDispersivity_; // alpha_l
+ values[1] = transverseDispersivity_; // alpha_t
+
+ return values;
+ }
+
+ bool useTwoPointGradient(const Element &elem,
+ int vertexI,
+ int vertexJ) const
+ {
+ return false;
+ }
+
+
+private:
+ bool isInLens_(const GlobalPosition &pos) const
+ {
+ for (int i = 0; i < dim; ++i) {
+ if (pos[i] < lensLowerLeft_[i] || pos[i] > lensUpperRight_[i])
+ return false;
+ }
+ return true;
+ }
+
+ GlobalPosition lensLowerLeft_;
+ GlobalPosition lensUpperRight_;
+
+ FieldMatrix lensK_;
+ FieldMatrix outerK_;
+ Scalar lensPorosity_;
+ Scalar outerPorosity_;
+ Scalar longitudinalDispersivity_;
+ Scalar transverseDispersivity_;
+};
+
+} // end namespace
+#endif
+
diff --git a/appl/lecture/msm/1p2cvs2p/lensspatialparameters2p.hh b/appl/lecture/msm/1p2cvs2p/lensspatialparameters2p.hh
new file mode 100644
index 0000000000..50bb1459f8
--- /dev/null
+++ b/appl/lecture/msm/1p2cvs2p/lensspatialparameters2p.hh
@@ -0,0 +1,191 @@
+// $Id$
+/*****************************************************************************
+ * Copyright (C) 2010 by Markus Wolff *
+ * Copyright (C) 2010 by Klaus Mosthaf *
+ * Copyright (C) 2007-2008 by Bernd Flemisch *
+ * Copyright (C) 2008-2009 by Andreas Lauser *
+ * Institute of Hydraulic Engineering *
+ * University of Stuttgart, Germany *
+ * email: <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, as long as this copyright notice *
+ * is included in its original form. *
+ * *
+ * This program is distributed WITHOUT ANY WARRANTY. *
+ *****************************************************************************/
+#ifndef DUMUX_LENSSPATIALPARAMETERS_2P_HH
+#define DUMUX_LENSSPATIALPARAMETERS_2P_HH
+
+#include <dumux/material/spatialparameters/boxspatialparameters.hh>
+#include <dumux/material/fluidmatrixinteractions/2p/regularizedvangenuchten.hh>
+#include <dumux/material/fluidmatrixinteractions/2p/regularizedbrookscorey.hh>
+#include <dumux/material/fluidmatrixinteractions/2p/linearmaterial.hh>
+#include <dumux/material/fluidmatrixinteractions/2p/efftoabslaw.hh>
+
+#include <dumux/boxmodels/2p/2pmodel.hh>
+
+/**
+ * @file
+ * @brief Class for defining spatial parameters
+ * @author Bernd Flemisch, Klaus Mosthaf, Markus Wolff
+ */
+
+namespace Dumux
+{
+
+/** \todo Please doc me! */
+
+template<class TypeTag>
+class LensSpatialParameters2p : public BoxSpatialParameters<TypeTag>
+{
+ typedef BoxSpatialParameters<TypeTag> ParentType;
+ 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;
+
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(TwoPIndices)) Indices;
+
+ enum {
+ dim=GridView::dimension,
+ dimWorld=GridView::dimensionworld,
+
+ wPhaseIdx = Indices::wPhaseIdx,
+ nPhaseIdx = Indices::nPhaseIdx
+ };
+
+ typedef Dune::FieldVector<CoordScalar,dimWorld> GlobalPosition;
+ typedef Dune::FieldMatrix<CoordScalar,dimWorld,dimWorld> FieldMatrix;
+
+ typedef typename GridView::template Codim<0>::Entity Element;
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(FVElementGeometry)) FVElementGeometry;
+
+ // define the material law which is parameterized by effective
+ // saturations
+// typedef RegularizedVanGenuchten<Scalar> EffectiveLaw;
+ typedef RegularizedBrooksCorey<Scalar> EffectiveLaw;
+
+public:
+ // define the material law parameterized by absolute saturations
+ typedef EffToAbsLaw<EffectiveLaw> MaterialLaw;
+ typedef typename MaterialLaw::Params MaterialLawParams;
+
+ LensSpatialParameters2p(const GridView& gridView)
+ : ParentType(gridView),
+ lensK_(0),
+ outerK_(0)
+ {
+ Dumux::InterfaceSoilProperties interfaceSoilProps("interface2p.xml");
+
+ lensPorosity_ = interfaceSoilProps.ISP_FinePorosity;
+ outerPorosity_ = interfaceSoilProps.ISP_CoarsePorosity;
+
+ for(int i = 0; i < dim; i++){
+ lensK_[i][i] = interfaceSoilProps.ISP_FinePermeability;
+ outerK_[i][i] = interfaceSoilProps.ISP_CoarsePermeability;
+ }
+
+ // residual saturations
+ lensMaterialParams_.setSwr(interfaceSoilProps.ISP_FineResidualSaturationWetting);
+ lensMaterialParams_.setSnr(interfaceSoilProps.ISP_FineResidualSaturationNonWetting);
+ outerMaterialParams_.setSwr(interfaceSoilProps.ISP_CoarseResidualSaturationWetting);
+ outerMaterialParams_.setSnr(interfaceSoilProps.ISP_CoarseResidualSaturationNonWetting);
+
+ // parameters for the Van Genuchten law
+ // alpha and n
+// lensMaterialParams_.setVgAlpha(0.00045);
+// lensMaterialParams_.setVgN(7.3);
+// outerMaterialParams_.setVgAlpha(0.0037);
+// outerMaterialParams_.setVgN(4.7);
+
+ // parameters for the Brooks-Corey law
+ lensMaterialParams_.setPe(interfaceSoilProps.ISP_FineBrooksCoreyEntryPressure);
+ lensMaterialParams_.setAlpha(interfaceSoilProps.ISP_FineBrooksCoreyLambda);
+ outerMaterialParams_.setPe(interfaceSoilProps.ISP_CoarseBrooksCoreyEntryPressure);
+ outerMaterialParams_.setAlpha(interfaceSoilProps.ISP_CoarseBrooksCoreyLambda);
+
+ // parameters for the linear law
+ // minimum and maximum pressures
+ // lensMaterialParams_.setEntryPC(0);
+// outerMaterialParams_.setEntryPC(0);
+// lensMaterialParams_.setMaxPC(0);
+// outerMaterialParams_.setMaxPC(0);
+ }
+
+ /*!
+ * \brief Apply the intrinsic permeability tensor to a pressure
+ * potential gradient.
+ *
+ * \param element The current finite element
+ * \param fvElemGeom The current finite volume geometry of the element
+ * \param scvIdx The index sub-control volume face where the
+ * intrinsic velocity ought to be calculated.
+ */
+ const FieldMatrix& intrinsicPermeability(const Element &element,
+ const FVElementGeometry &fvElemGeom,
+ int scvIdx) const
+ {
+ const GlobalPosition &globalPos = fvElemGeom.subContVol[scvIdx].global;
+ if (isInLens_(globalPos))
+ return lensK_;
+ return outerK_;
+ }
+
+ Scalar porosity(const Element &element,
+ const FVElementGeometry &fvElemGeom,
+ int scvIdx) const
+ {
+ const GlobalPosition &globalPos = fvElemGeom.subContVol[scvIdx].global;
+ if (isInLens_(globalPos))
+ return lensPorosity_;
+ return outerPorosity_;
+ }
+
+ // return the brooks-corey context depending on the position
+ const MaterialLawParams& materialLawParams(const Element &element,
+ const FVElementGeometry &fvElemGeom,
+ int scvIdx) const
+ {
+ const GlobalPosition &globalPos = fvElemGeom.subContVol[scvIdx].global;
+
+ if (isInLens_(globalPos))
+ return lensMaterialParams_;
+ return outerMaterialParams_;
+ }
+
+
+ //! Set the bounding box of the fine-sand lens
+ void setLensCoords(const GlobalPosition& lensLowerLeft,
+ const GlobalPosition& lensUpperRight)
+ {
+ lensLowerLeft_ = lensLowerLeft;
+ lensUpperRight_ = lensUpperRight;
+ }
+
+private:
+ bool isInLens_(const GlobalPosition &pos) const
+ {
+ for (int i = 0; i < dim; ++i) {
+ if (pos[i] < lensLowerLeft_[i] || pos[i] > lensUpperRight_[i])
+ return false;
+ }
+ return true;
+ }
+
+ GlobalPosition lensLowerLeft_;
+ GlobalPosition lensUpperRight_;
+
+ FieldMatrix lensK_;
+ FieldMatrix outerK_;
+ Scalar lensPorosity_;
+ Scalar outerPorosity_;
+ MaterialLawParams lensMaterialParams_;
+ MaterialLawParams outerMaterialParams_;
+};
+
+} // end namespace
+#endif
+
diff --git a/appl/lecture/msm/1p2cvs2p/water_contaminant.hh b/appl/lecture/msm/1p2cvs2p/water_contaminant.hh
new file mode 100644
index 0000000000..54b8e4af71
--- /dev/null
+++ b/appl/lecture/msm/1p2cvs2p/water_contaminant.hh
@@ -0,0 +1,142 @@
+/*****************************************************************************
+ * Copyright (C) 2010 by Klaus Mosthaf *
+ * Copyright (C) 2010 by Bernd Flemisch *
+ * 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, as long as this copyright notice *
+ * is included in its original form. *
+ * *
+ * This program is distributed WITHOUT ANY WARRANTY. *
+ *****************************************************************************/
+/*!
+ * \file
+ *
+ * \brief A fluid system with one phase and an arbitrary number of components.
+ */
+#ifndef WATERCONTAMINANT_HH
+#define WATERCONTAMINANT_HH
+
+#include <dune/common/exceptions.hh>
+
+#include <dumux/common/propertysystem.hh>
+#include <dumux/boxmodels/1p2c/1p2cproperties.hh>
+
+namespace Dumux
+{
+
+namespace Properties
+{
+ NEW_PROP_TAG(Scalar);
+ NEW_PROP_TAG(OnePTwoCIndices);
+};
+
+/*!
+ * \brief A fluid system with one phase and an arbitrary number of components.
+ */
+template <class TypeTag, bool verbose=true>
+class WaterContaminant
+{
+ typedef WaterContaminant<TypeTag, verbose> ThisType;
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(Scalar)) Scalar;
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(OnePTwoCIndices)) Indices;
+
+public:
+ enum {
+ // component indices
+ waterIdx = 0,
+ contaminantIdx = 1
+ };
+
+ static void init()
+ {}
+
+ /*!
+ * \brief Return the human readable name of a component
+ */
+ static const char *componentName(int compIdx)
+ {
+ switch(compIdx)
+ {
+ case waterIdx:
+ return "Water";
+ case contaminantIdx:
+ return "Contaminant";
+ default:
+ DUNE_THROW(Dune::InvalidStateException, "Invalid component index " << compIdx);
+ }
+ }
+
+ /*!
+ * \brief Return the molar mass of a component [kg/mol].
+ */
+ static Scalar molarMass(int compIdx)
+ { return 18e-3; }
+
+ /*!
+ * \brief Given all mole fractions in a phase, return the phase
+ * density [kg/m^3].
+ */
+ template <class FluidState>
+ static Scalar phaseDensity(int phaseIdx,
+ Scalar temperature,
+ Scalar pressure,
+ const FluidState &fluidState)
+ {
+ if (phaseIdx == 0)
+ return 1000;//constDensity_; // in [kg /m^3]
+
+ DUNE_THROW(Dune::InvalidStateException, "Invalid phase index " << phaseIdx);
+ }
+
+ /*!
+ * \brief Return the dynamic viscosity of a phase.
+ */
+ template <class FluidState>
+ static Scalar phaseViscosity(int phaseIdx,
+ Scalar temperature,
+ Scalar pressure,
+ const FluidState &fluidState)
+ {
+ if (phaseIdx == 0)
+ return 1e-03;
+
+ DUNE_THROW(Dune::InvalidStateException, "Invalid phase index " << phaseIdx);
+ }
+
+ /*!
+ * \brief Given all mole fractions, return the diffusion
+ * coefficent of a component in a phase.
+ */
+ template <class FluidState>
+ static Scalar diffCoeff(int phaseIdx,
+ int compIIdx,
+ int compJIdx,
+ Scalar temperature,
+ Scalar pressure,
+ const FluidState &fluidState)
+ {
+ return 1e-9; // in [m^2/s]
+ }
+
+
+ WaterContaminant( )
+ {
+ //load interface-file
+ Dumux::InterfaceFluidProperties interfaceFluidProps("interface1p2c.xml");
+
+ diffCoefficient_ = interfaceFluidProps.IFP_MolecularDiffusionCoefficient;
+ }
+
+
+//private:
+ static Scalar diffCoefficient_;
+};
+
+} // end namepace
+
+#endif
diff --git a/appl/lecture/msm/Makefile.am b/appl/lecture/msm/Makefile.am
new file mode 100755
index 0000000000..1952b9b290
--- /dev/null
+++ b/appl/lecture/msm/Makefile.am
@@ -0,0 +1,8 @@
+SUBDIRS = . \
+ 1p2cvs2p \
+ buckleyleverett \
+ mcwhorter
+
+
+include $(top_srcdir)/am/global-rules
+
diff --git a/appl/lecture/msm/buckleyleverett/Makefile.am b/appl/lecture/msm/buckleyleverett/Makefile.am
new file mode 100644
index 0000000000..0c64c7b542
--- /dev/null
+++ b/appl/lecture/msm/buckleyleverett/Makefile.am
@@ -0,0 +1,8 @@
+dist_noinst_DATA = buckleyleverett_ff
+
+noinst_PROGRAMS = buckleyleverett_ff
+
+buckleyleverett_ff_SOURCES = buckleyleverett_ff.cc
+buckleyleverett_ff_LDFLAGGS = static
+
+include $(top_srcdir)/am/global-rules
diff --git a/appl/lecture/msm/buckleyleverett/buckleyleverett_analytic.hh b/appl/lecture/msm/buckleyleverett/buckleyleverett_analytic.hh
new file mode 100644
index 0000000000..10185e563c
--- /dev/null
+++ b/appl/lecture/msm/buckleyleverett/buckleyleverett_analytic.hh
@@ -0,0 +1,377 @@
+// $Id$
+
+#ifndef DUMUX_BUCKLEYLEVERETT_ANALYTICAL_HH
+#define DUMUX_BUCKLEYLEVERETT_ANALYTICAL_HH
+
+#include <dumux/material/fluidmatrixinteractions/2p/linearmaterial.hh>
+
+/**
+ * @file
+ * @brief Analytical solution of the buckley-leverett problem
+ * @author Markus Wolff
+ */
+
+namespace Dumux
+{
+/**
+ * \ingroup fracflow
+ * @brief IMplicit Pressure Explicit Saturation (IMPES) scheme for the solution of
+ * the Buckley-Leverett problem
+ */
+
+template<class Scalar, class Law> struct CheckMaterialLaw
+{
+ static bool isLinear()
+ {
+ return false;
+ }
+};
+
+template<class Scalar> struct CheckMaterialLaw<Scalar, LinearMaterial<Scalar> >
+{
+ static bool isLinear()
+ {
+ return true;
+ }
+};
+
+template<class TypeTag> class BuckleyLeverettAnalytic
+{
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(Problem)) Problem;
+ 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 GET_PROP_TYPE(TypeTag, PTAG(SpatialParameters)) SpatialParameters;
+ typedef typename SpatialParameters::MaterialLaw MaterialLaw;
+ typedef typename MaterialLaw::Params MaterialLawParams;
+
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(FluidSystem)) FluidSystem;
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(FluidState)) FluidState;
+
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(TwoPIndices)) Indices;
+
+ enum
+ {
+ dim = GridView::dimension, dimworld = GridView::dimensionworld
+ };
+ enum
+ {
+ wPhaseIdx = Indices::wPhaseIdx, nPhaseIdx = Indices::nPhaseIdx
+ };
+
+ typedef Dune::BlockVector<Dune::FieldVector<Scalar, 1> > BlockVector;
+ typedef typename GridView::Traits::template Codim<0>::Entity Element;
+ typedef typename GridView::template Codim<0>::Iterator ElementIterator;
+ typedef Dune::FieldVector<Scalar, dimworld> GlobalPosition;
+
+public:
+
+ // functions needed for analytical solution
+
+ void initializeAnalytic()
+ {
+ analyticSolution_.resize(size_);
+ analyticSolution_ = 0;
+ error_.resize(size_);
+ error_ = 0;
+ elementVolume_.resize(size_);
+ elementVolume_ = 0;
+
+ return;
+ }
+
+ void prepareAnalytic()
+ {
+ const MaterialLawParams& materialLawParams(problem_.spatialParameters().materialLawParams(dummyGlobal_, dummyElement_));
+
+ swr_ = materialLawParams.Swr();
+ snr_ = materialLawParams.Snr();
+ porosity_ = problem_.spatialParameters().porosity(dummyGlobal_, dummyElement_);
+
+ time_ = 0;
+
+ satVec_ = swr_;
+ for (int i = 1; i < pointNum_; i++)
+ {
+ satVec_[i] = satVec_[i - 1] + (1 - snr_ - swr_) / intervalNum_;
+ }
+// std::cout<<"satVec_ = "<<satVec_<<std::endl;
+
+ FluidState fluidState;
+ Scalar temp = problem_.temperature(dummyGlobal_, dummyElement_);
+ Scalar press = problem_.referencePressure(dummyGlobal_, dummyElement_);
+ Scalar viscosityW = FluidSystem::phaseViscosity(wPhaseIdx, temp, press, fluidState);
+ Scalar viscosityNW = FluidSystem::phaseViscosity(nPhaseIdx, temp, press, fluidState);
+
+ for (int i = 0; i < pointNum_; i++)
+ {
+ fractionalW_[i] = MaterialLaw::krw(materialLawParams, satVec_[i])/viscosityW;
+ fractionalW_[i] /= (fractionalW_[i] + MaterialLaw::krn(materialLawParams, satVec_[i])/viscosityNW);
+ }
+// std::cout<<"fractionalW_ = "<<fractionalW_<<std::endl;
+
+ dfwdsw_ = 0;
+ for (int i = 1; i < intervalNum_; i++)
+ {
+ dfwdsw_[i] = (fractionalW_[i + 1] - fractionalW_[i - 1]) / (satVec_[i + 1] - satVec_[i - 1]);
+ }
+// std::cout<<"dfwdsw_ = "<<dfwdsw_<<std::endl;
+
+ for (int i = 0; i < pointNum_; i++)
+ {
+ if (dfwdsw_[i] > dfwdsw_[i + 1])
+ {
+ dfwdswmax_ = i;
+ break;
+ }
+ }
+
+ return;
+ }
+
+ void calcSatError(BlockVector &Approx)
+ {
+ error_ = 0;
+ elementVolume_ = 0;
+
+ ElementIterator eItEnd = problem_.gridView().template end<0> ();
+ for (ElementIterator eIt = problem_.gridView().template begin<0> (); eIt != eItEnd; ++eIt)
+ {
+ // get entity
+ const Element& element = *eIt;
+ int index = problem_.variables().index(*eIt);
+ elementVolume_[index] = element.geometry().volume();
+ // std::cout<<"elementVolume_ = "<<elementVolume_[index]<<std::endl;
+ }
+
+ Scalar globalVolume = elementVolume_.one_norm();
+ // std::cout<<"globalVolume = "<<globalVolume<<std::endl;
+
+ for (int i = 0; i < size_; i++)
+ {
+ error_[i] = analyticSolution_[i] - Approx[i];
+ // std::cout<<"error_ = "<<error_[i]<<std::endl;
+ // std::cout<<"analyticSolution_ = "<<analyticSolution_[i]<<std::endl;
+ // std::cout<<"Approx = "<<Approx[i]<<std::endl;
+ }
+ // std::cout<<"error_ = "<<error_<<std::endl;
+
+ Scalar diffNorm = error_.two_norm();
+ // std::cout<<"diffNorm = "<<diffNorm<<std::endl;
+
+ for (int i = 0; i < size_; i++)
+ {
+ error_[i] = diffNorm * pow((elementVolume_[i] / globalVolume), 0.5);
+ }
+
+ return;
+ }
+
+ void updateExSol()
+ {
+ // position of the fluid front
+ xf_ = 0;
+ for (int i = 0; i < pointNum_; i++)
+ {
+ xf_[i] = vTot_ * time_ / porosity_ * dfwdsw_[i];
+ }
+
+// std::cout<<"xf_ = "<<xf_<<std::endl;
+ int xhelp = pointNum_ / 3;
+ int xhelpold = 0;
+ int xhelpoldold = 0;
+ int xfmax = 0;
+
+ // position of maximum xf_
+ for (int i = 0; i < pointNum_; i++)
+ {
+ if (xf_[i] > xf_[i + 1])
+ {
+ xfmax = i;
+ break;
+ }
+ }
+
+ // balancing of the areas ahead of the front and below the curve
+ bool a = true;
+ Scalar A1;
+ Scalar A2;
+ Scalar b;
+ int xhelp2 = 0;
+
+ while (a)
+ {
+ if (CheckMaterialLaw<Scalar, MaterialLaw>::isLinear())
+ break;
+
+ A1 = 0;
+
+ for (int i = 0; i <= xhelp - 1; i++)
+ {
+ A1 += (satVec_[i] - swr_ + satVec_[i + 1] - swr_) * 0.5 * (xf_[i + 1] - xf_[i]);
+ }
+
+ A2 = 0;
+
+ for (int i = xhelp; i <= xfmax - 1; i++)
+ {
+ A2 += (satVec_[xfmax] - satVec_[i] + satVec_[xfmax] - satVec_[i + 1]) * 0.5 * (xf_[i + 1] - xf_[i]);
+ }
+
+ b = xf_[xfmax];
+ xhelp2 = xfmax;
+
+ while (b > xf_[xhelp])
+ {
+ xhelp2 += 1;
+ b = xf_[xhelp2];
+ }
+
+ for (int i = xfmax; i <= xhelp2; i++)
+ {
+ A2 += (satVec_[i] - satVec_[xfmax] + satVec_[i + 1] - satVec_[xfmax]) * 0.5 * (xf_[i] - xf_[i + 1]);
+ }
+
+ xhelpoldold = xhelpold;
+ xhelpold = xhelp;
+
+ if (fabs(A1) > fabs(A2))
+ {
+ xhelp = xhelp - 1;
+ }
+
+ if (fabs(A1) < fabs(A2))
+ {
+ xhelp = xhelp + 1;
+ }
+
+ if (xhelp == xhelpoldold)
+ {
+ a = false;
+ }
+ }
+ // std::cout<<"xf_ = "<<xf_<<std::endl;
+
+ // iterate over vertices and get analytic saturation solution
+ ElementIterator eItEnd = problem_.gridView().template end<0> ();
+ for (ElementIterator eIt = problem_.gridView().template begin<0> (); eIt != eItEnd; ++eIt)
+ {
+ // get global coordinate of cell center
+ const GlobalPosition& globalPos = eIt->geometry().center();
+
+ // find index of current vertex
+ int index = problem_.variables().index(*eIt);
+
+ // account for linear material law
+ if (CheckMaterialLaw<Scalar, MaterialLaw>::isLinear())
+ {
+ if (globalPos[0] <= xf_[1])
+ {
+ analyticSolution_[index] = 1 - snr_;
+ }
+ if (globalPos[0] > xf_[1])
+ {
+ analyticSolution_[index] = swr_;
+ }
+ }
+
+ // non-linear material law
+
+ else
+ {
+ // find x_f next to global coordinate of the vertex
+ int xnext = 0;
+ for (int i = intervalNum_; i >= 0; i--)
+ {
+ if (globalPos[0] < xf_[i])
+ {
+ xnext = i;
+ break;
+ }
+ }
+
+ // account for the area not yet reached by the front
+ if (globalPos[0] > xf_[xhelp2])
+ {
+ analyticSolution_[index] = swr_;
+ continue;
+ }
+
+ if (globalPos[0] <= xf_[xhelp2])
+ {
+ analyticSolution_[index] = satVec_[xnext];
+ continue;
+ }
+ }
+ }
+
+ // call function to calculate the saturation error_
+ calcSatError(problem_.variables().saturation());
+
+ return;
+ }
+
+ void calculateAnalyticSolution()
+ {
+ time_ += problem_.timeManager().timeStepSize();
+
+ updateExSol();
+ }
+
+ //Write saturation and pressure into file
+ template<class MultiWriter>
+ void addOutputVtkFields(MultiWriter &writer)
+ {
+ BlockVector *analyticSolution = writer.template createField<Scalar, 1> (size_);
+ BlockVector *error = writer.template createField<Scalar, 1> (size_);
+
+ *analyticSolution = analyticSolution_;
+ *error = error_;
+
+ writer.addCellData(analyticSolution, "saturation (exact solution)");
+ writer.addCellData(error, "error_");
+
+ return;
+ }
+
+ //! Construct an IMPES object.
+ BuckleyLeverettAnalytic(Problem& problem, Scalar totalVelocity = 3e-7) :
+ problem_(problem), analyticSolution_(0), error_(0), elementVolume_(0), size_(problem.gridView().size(0)), vTot_(totalVelocity), dummyElement_(
+ *(problem_.gridView().template begin<0> ())), dummyGlobal_(GlobalPosition(1))
+ {
+ initializeAnalytic();
+ prepareAnalytic();
+ }
+
+protected:
+
+
+private:
+ Problem& problem_;
+
+ BlockVector analyticSolution_;
+ BlockVector error_;
+ BlockVector elementVolume_;
+
+ Scalar time_;
+ int size_;
+ Scalar swr_;
+ Scalar snr_;
+ Scalar porosity_;
+ Scalar vTot_;
+ enum
+ {
+ intervalNum_ = 1000, pointNum_ = intervalNum_ + 1
+ };
+ Dune::FieldVector<Scalar, pointNum_> satVec_;
+ Dune::FieldVector<Scalar, pointNum_> fractionalW_;
+ Dune::FieldVector<Scalar, pointNum_> dfwdsw_;
+ Dune::FieldVector<Scalar, pointNum_> xf_;
+ int dfwdswmax_;
+ const Element& dummyElement_;
+ const GlobalPosition& dummyGlobal_;
+
+};
+}
+#endif
diff --git a/appl/lecture/msm/buckleyleverett/buckleyleverett_ff.cc b/appl/lecture/msm/buckleyleverett/buckleyleverett_ff.cc
new file mode 100644
index 0000000000..5379aaa3ee
--- /dev/null
+++ b/appl/lecture/msm/buckleyleverett/buckleyleverett_ff.cc
@@ -0,0 +1,152 @@
+// $Id: test_impes.cc 4144 2010-08-24 10:10:47Z bernd $
+/*****************************************************************************
+ * Copyright (C) 2010 by Markus Wolff, Klaus Mosthaf *
+ * Copyright (C) 2007-2008 by Bernd Flemisch *
+ * Copyright (C) 2008-2009 by Andreas Lauser *
+ * Institute of Hydraulic Engineering *
+ * University of Stuttgart, Germany *
+ * email: <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, as long as this copyright notice *
+ * is included in its original form. *
+ * *
+ * This program is distributed WITHOUT ANY WARRANTY. *
+ *****************************************************************************/
+#include "config.h"
+
+#include "external_interface.hh"
+#include "buckleyleverettproblem.hh"
+
+#include <dune/grid/common/gridinfo.hh>
+
+#include <dune/common/exceptions.hh>
+#include <dune/common/mpihelper.hh>
+
+#include <iostream>
+#include <boost/format.hpp>
+
+#include <iomanip>
+
+//#include "dumux/decoupled/2p/old_fractionalflow/variableclass2p.hh"
+//#include "dumux/decoupled/2p/old_fractionalflow/define2pmodel.hh"
+//#include "buckleyleverettproblem.hh"
+//#include "impes_buckleyleverett_analytic.hh"
+//#include "dumux/timedisc/expliciteulerstep.hh"
+//#include "dumux/timedisc/timeloop.hh"
+
+void usage(const char *progname)
+{
+ std::cout << boost::format("usage: %s tEnd\n")%progname;
+ exit(1);
+}
+
+
+////////////////////////
+// the main function
+////////////////////////
+int main(int argc, char** argv)
+{
+ try
+ {
+ typedef TTAG(BuckleyLeverettProblem) TypeTag;
+ typedef GET_PROP_TYPE(TypeTag, PTAG(Scalar)) Scalar;
+ typedef GET_PROP_TYPE(TypeTag, PTAG(Grid)) Grid;
+ typedef GET_PROP_TYPE(TypeTag, PTAG(Problem)) Problem;
+ typedef Dune::FieldVector<Scalar, Grid::dimensionworld> GlobalPosition;
+
+ static const int dim = Grid::dimension;
+
+ // initialize MPI, finalize is done automatically on exit
+ Dune::MPIHelper::instance(argc, argv);
+
+
+ //load interface-file
+ Dumux::InterfaceProblemProperties interfaceProbProps("interface_BL.xml");
+ double discretizationLength = interfaceProbProps.IPP_DiscretizationLength;
+
+ // define the problem dimensions
+ Dune::FieldVector<Scalar, dim> lowerLeft(0);
+ Dune::FieldVector<Scalar, dim> upperRight(300);
+ upperRight[1] = 75;
+
+ int cellNumberX = 300/discretizationLength;
+ int cellNumberY = 75/discretizationLength;
+
+ Dune::FieldVector<int, dim> cellNumbers(cellNumberX);
+ cellNumbers[1] = cellNumberY;
+ if (argc != 2)
+ usage(argv[0]);
+
+ std::string arg1(argv[1]);
+ std::istringstream is1(arg1);
+ double tEnd;
+ is1 >> tEnd;
+
+ double dt = tEnd;
+
+ // create a grid object
+ Grid grid(cellNumbers, lowerLeft, upperRight);
+ Dune::gridinfo(grid);
+
+ ////////////////////////////////////////////////////////////
+ // instantiate and run the concrete problem
+ ////////////////////////////////////////////////////////////
+
+ Dumux::InterfaceFluidProperties interfaceFluidProps("interface_BL.xml");
+ typedef GET_PROP_TYPE(TypeTag, PTAG(WettingPhase)) WettingPhase;
+ typedef GET_PROP_TYPE(TypeTag, PTAG(NonwettingPhase)) NonwettingPhase;
+
+ WettingPhase::Component::setViscosity(interfaceFluidProps.IFP_ViscosityWettingFluid);
+ NonwettingPhase::Component::setViscosity(interfaceFluidProps.IFP_ViscosityNonWettingFluid);
+
+ WettingPhase::Component::setDensity(interfaceFluidProps.IFP_DensityWettingFluid);
+ NonwettingPhase::Component::setDensity(interfaceFluidProps.IFP_DensityNonWettingFluid);
+
+ Problem problem(grid.leafView(), lowerLeft, upperRight);
+
+ problem.timeManager().init(problem, 0, dt, tEnd);
+ problem.timeManager().run();
+ return 0;
+ }
+ catch (Dune::Exception &e) {
+ std::cerr << "Dune reported error: " << e << std::endl;
+ }
+ catch (...) {
+ std::cerr << "Unknown exception thrown!\n";
+ throw;
+ }
+
+// OLD STUFF -> PROBLEM
+// // IMPES parameters
+// int iterFlag = 0;
+// int nIter = 30;
+// double maxDefect = 1e-5;
+//
+// // plotting parameters
+// const char* fileName = "buckleyleverett";
+// int modulo = 1;
+//
+// typedef Dumux::FVVelocity2P<GridView, Scalar, VariableType,
+// ProblemType> DiffusionType;
+// DiffusionType diffusion(gridView, problem, modelDef);
+//
+// typedef Dumux::FVSaturation2P<GridView, Scalar, VariableType,
+// ProblemType> TransportType;
+// TransportType transport(gridView, problem, modelDef);
+//
+// typedef Dune::IMPESBLAnalytic<GridView, DiffusionType, TransportType,
+// VariableType> IMPESType;
+// IMPESType impes(diffusion, transport, iterFlag, nIter, maxDefect);
+// Dumux::ExplicitEulerStep<GridType, IMPESType> timestep;
+// Dumux::TimeLoop<GridView, IMPESType> timeloop(gridView, tStart, tEnd, fileName,
+// modulo, cFLFactor, maxDt, firstDt, timestep);
+// Dune::Timer timer;
+// timer.reset();
+// timeloop.execute(impes, false);
+// std::cout << "timeloop.execute took " << timer.elapsed() << " seconds"
+// << std::endl;
+// //printvector(std::cout, *fractionalflow, "saturation", "row", 200, 1);
+}
diff --git a/appl/lecture/msm/buckleyleverett/buckleyleverett_spatialparams.hh b/appl/lecture/msm/buckleyleverett/buckleyleverett_spatialparams.hh
new file mode 100644
index 0000000000..150d22ece9
--- /dev/null
+++ b/appl/lecture/msm/buckleyleverett/buckleyleverett_spatialparams.hh
@@ -0,0 +1,111 @@
+// $Id: buckleyleverett_spatialparams.hh 4144 2010-08-24 10:10:47Z bernd $
+/*****************************************************************************
+ * Copyright (C) 2008-2009 by Markus Wolff *
+ * Copyright (C) 2010 by Klaus Mosthaf *
+ * 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, as long as this copyright notice *
+ * is included in its original form. *
+ * *
+ * This program is distributed WITHOUT ANY WARRANTY. *
+ *****************************************************************************/
+#ifndef BUCKLEYLEVERETT_SPATIALPARAMETERS_HH
+#define BUCKLEYLEVERETT_SPATIALPARAMETERS_HH
+
+
+#include <dumux/material/fluidmatrixinteractions/2p/linearmaterial.hh>
+#include <dumux/material/fluidmatrixinteractions/2p/regularizedbrookscorey.hh>
+#include <dumux/material/fluidmatrixinteractions/2p/efftoabslaw.hh>
+
+namespace Dumux
+{
+
+/** \todo Please doc me! */
+
+template<class TypeTag>
+class BuckleyLeverettSpatialParams
+{
+ 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 RegularizedBrooksCorey<Scalar> RawMaterialLaw;
+// 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
+ {
+ return constPermeability_;
+ }
+
+ Scalar porosity(const GlobalPosition& globalPos, const Element& element) const
+ {
+ return porosity_;
+ }
+
+
+ // return the brooks-corey context depending on the position
+ const MaterialLawParams& materialLawParams(const GlobalPosition& globalPos, const Element &element) const
+ {
+ return materialLawParams_;
+ }
+
+
+ BuckleyLeverettSpatialParams(const GridView& gridView)
+ : constPermeability_(0)
+ {
+ Dumux::InterfaceSoilProperties interfaceSoilProps("interface_BL.xml");
+
+ // residual saturations
+ materialLawParams_.setSwr(interfaceSoilProps.ISP_ResidualSaturationWetting);
+ materialLawParams_.setSnr(interfaceSoilProps.ISP_ResidualSaturationNonWetting);
+
+ porosity_ = interfaceSoilProps.ISP_Porosity;
+
+ // parameters for the Brooks-Corey Law
+ // entry pressures
+ materialLawParams_.setPe(interfaceSoilProps.ISP_BrooksCoreyEntryPressure);
+ // Brooks-Corey shape parameters
+ materialLawParams_.setAlpha(interfaceSoilProps.ISP_BrooksCoreyLambda);
+
+ // parameters for the linear
+ // entry pressures function
+// materialLawParams_.setEntryPC(0);
+// materialLawParams_.setMaxPC(0);
+// materialLawParams_.setEntryPC(2);
+// materialLawParams_.setMaxPC(3);
+
+ for(int i = 0; i < dim; i++)
+ constPermeability_[i][i] = interfaceSoilProps.ISP_Permeability;
+ }
+
+private:
+ MaterialLawParams materialLawParams_;
+ FieldMatrix constPermeability_;
+ Scalar porosity_;
+
+};
+
+} // end namespace
+#endif
diff --git a/appl/lecture/msm/buckleyleverett/buckleyleverettproblem.hh b/appl/lecture/msm/buckleyleverett/buckleyleverettproblem.hh
new file mode 100644
index 0000000000..e1a8debf77
--- /dev/null
+++ b/appl/lecture/msm/buckleyleverett/buckleyleverettproblem.hh
@@ -0,0 +1,317 @@
+// $Id: test_impes_problem.hh 4209 2010-09-01 12:31:45Z mwolff $
+/*****************************************************************************
+ * Copyright (C) 2007-2008 by Markus Wolff, Klaus Mosthaf *
+ * Copyright (C) 2007-2008 by Bernd Flemisch *
+ * Copyright (C) 2008-2009 by Andreas Lauser *
+ * Institute of Hydraulic Engineering *
+ * University of Stuttgart, Germany *
+ * email: <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, as long as this copyright notice *
+ * is included in its original form. *
+ * *
+ * This program is distributed WITHOUT ANY WARRANTY. *
+ *****************************************************************************/
+#ifndef DUMUX_BUCKLEYLEVERETTPROBLEM_HH
+#define DUMUX_BUCKLEYLEVERETTPROBLEM_HH
+
+#include <dune/grid/sgrid.hh>
+
+#include <dumux/material/fluidsystems/liquidphase.hh>
+//#include <dumux/material/components/simpleh2o.hh>
+//#include <dumux/material/components/oil.hh>
+#include <pseudooil.hh>
+#include <pseudoh2o.hh>
+
+#include <dumux/decoupled/2p/impes/impesproblem2p.hh>
+#include <dumux/decoupled/2p/diffusion/fv/fvvelocity2p.hh>
+#include <dumux/decoupled/2p/diffusion/fvmpfa/fvmpfaovelocity2p.hh>
+#include <dumux/decoupled/2p/transport/fv/fvsaturation2p.hh>
+#include <dumux/decoupled/2p/transport/fv/capillarydiffusion.hh>
+#include <dumux/decoupled/2p/transport/fv/gravitypart.hh>
+
+#include "buckleyleverett_spatialparams.hh"
+#include "buckleyleverett_analytic.hh"
+
+namespace Dumux
+{
+template <class TypeTag>
+class BuckleyLeverettProblem;
+
+//////////
+// Specify the properties
+//////////
+namespace Properties
+{
+NEW_TYPE_TAG(BuckleyLeverettProblem, INHERITS_FROM(DecoupledTwoP, MPFAProperties, Transport));
+
+// Set the grid type
+SET_PROP(BuckleyLeverettProblem, Grid)
+{
+ // typedef Dune::YaspGrid<2> type;
+ typedef Dune::SGrid<2, 2> type;
+};
+
+// Set the problem property
+SET_PROP(BuckleyLeverettProblem, Problem)
+{
+public:
+ typedef Dumux::BuckleyLeverettProblem<TypeTag> type;
+};
+
+// Set the model properties
+SET_PROP(BuckleyLeverettProblem, TransportModel)
+{
+ typedef Dumux::FVSaturation2P<TTAG(BuckleyLeverettProblem)> type;
+};
+SET_TYPE_PROP(BuckleyLeverettProblem, DiffusivePart, Dumux::CapillaryDiffusion<TypeTag>);
+SET_TYPE_PROP(BuckleyLeverettProblem, ConvectivePart, Dumux::GravityPart<TypeTag>);
+
+SET_PROP(BuckleyLeverettProblem, PressureModel)
+{
+ typedef Dumux::FVVelocity2P<TTAG(BuckleyLeverettProblem)> type;
+// typedef Dumux::FVMPFAOVelocity2P<TTAG(BuckleyLeverettProblem)> type;
+};
+
+//SET_INT_PROP(BuckleyLeverettProblem, VelocityFormulation,
+// GET_PROP_TYPE(TypeTag, PTAG(TwoPIndices))::velocityW);
+
+//SET_INT_PROP(BuckleyLeverettProblem, PressureFormulation,
+// GET_PROP_TYPE(TypeTag, PTAG(TwoPIndices))::pressureGlobal);
+
+// Set the wetting phase
+SET_PROP(BuckleyLeverettProblem, WettingPhase)
+{
+private:
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(Scalar)) Scalar;
+public:
+ typedef Dumux::LiquidPhase<Scalar, Dumux::PseudoOil<Scalar> > type;
+};
+
+// Set the non-wetting phase
+SET_PROP(BuckleyLeverettProblem, NonwettingPhase)
+{
+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(BuckleyLeverettProblem, SpatialParameters)
+{
+private:
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(Grid)) Grid;
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(Scalar)) Scalar;
+
+public:
+ typedef Dumux::BuckleyLeverettSpatialParams<TypeTag> type;
+};
+
+// Enable gravity
+SET_BOOL_PROP(BuckleyLeverettProblem, EnableGravity, false);
+
+SET_SCALAR_PROP(BuckleyLeverettProblem, CFLFactor, 0.95); //0.95
+}
+
+/*!
+ * \ingroup DecoupledProblems
+ */
+template<class TypeTag = TTAG(BuckleyLeverettProblem)>
+class BuckleyLeverettProblem: public IMPESProblem2P<TypeTag, BuckleyLeverettProblem<TypeTag> >
+{
+ typedef BuckleyLeverettProblem<TypeTag> ThisType;
+ typedef IMPESProblem2P<TypeTag, ThisType> ParentType;
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(GridView)) GridView;
+
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(TwoPIndices)) Indices;
+
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(FluidSystem)) FluidSystem;
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(FluidState)) FluidState;
+
+ enum
+ {
+ dim = GridView::dimension, dimWorld = GridView::dimensionworld
+ };
+ enum
+ {
+ wetting = 0, nonwetting = 1
+ };
+ 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;
+
+public:
+ BuckleyLeverettProblem(const GridView &gridView,
+ const GlobalPosition lowerLeft = 0,
+ const GlobalPosition upperRight = 0,
+ const Scalar pleftbc = 2e5)
+ : ParentType(gridView),
+ lowerLeft_(lowerLeft),
+ upperRight_(upperRight),
+ eps_(1e-8),
+ pLeftBc_(pleftbc),
+ analyticSolution_(*this, 3e-7)
+ {
+ //load interface-file
+ Dumux::InterfaceFluidProperties interfaceFluidProps("interface_BL.xml");
+
+ densityNonWetting_ = interfaceFluidProps.IFP_DensityNonWettingFluid;
+ }
+
+
+
+ /*!
+ * \name Problem parameters
+ */
+ // \{
+
+ /*!
+ * \brief The problem name.
+ *
+ * This is used as a prefix for files generated by the simulation.
+ */
+ const char *name() const
+ {
+ std::string simName = "buckleyleverett_linear_run";
+ Dumux::InterfaceProblemProperties interfaceProbProps("interface_BL.xml");
+ Scalar simNum = interfaceProbProps.IPP_SimulationNumber;
+
+ return (str(boost::format("%s-%02d")
+ %simName%simNum).c_str());
+ }
+
+ bool shouldWriteRestartFile() const
+ {
+ return false;
+ }
+
+ void postTimeStep()
+ {
+ analyticSolution_.calculateAnalyticSolution();
+
+ ParentType::postTimeStep();
+ };
+
+ void addOutputVtkFields()
+ {
+ ParentType::addOutputVtkFields();
+ analyticSolution_.addOutputVtkFields(this->resultWriter());
+ }
+
+ //! Write the fields current solution into an VTK output file.
+ void writeOutput()
+ {
+// if (this->timeManager().time() > 43.1999e6)
+// {
+ if (this->gridView().comm().rank() == 0)
+ std::cout << "Writing result file for current time step\n";
+
+ this->resultWriter().beginTimestep(this->timeManager().time() + this->timeManager().timeStepSize(),
+ this->gridView());
+ this->asImp_().addOutputVtkFields();
+ this->resultWriter().endTimestep();
+// }
+ }
+
+ /*!
+ * \brief Returns the temperature within the domain.
+ *
+ * This problem assumes a temperature of 10 degrees Celsius.
+ */
+ Scalar temperature(const GlobalPosition& globalPos, const Element& element) const
+ {
+ return 273.15 + 10; // -> 10°C
+ }
+
+ // \}
+
+
+ Scalar referencePressure(const GlobalPosition& globalPos, const Element& element) const
+ {
+ return 1e5;
+ }
+
+ std::vector<Scalar> source(const GlobalPosition& globalPos, const Element& element)
+ {
+ return std::vector<Scalar>(2, 0.0);
+ }
+
+ typename BoundaryConditions::Flags bctypePress(const GlobalPosition& globalPos, const Intersection& intersection) const
+ {
+ if (globalPos[0] < eps_)
+ return BoundaryConditions::dirichlet;
+
+ // all other boundaries
+ else
+ return BoundaryConditions::neumann;
+ }
+
+ BoundaryConditions::Flags bctypeSat(const GlobalPosition& globalPos, const Intersection& intersection) const
+ {
+ if (globalPos[0]> (upperRight_[0] - eps_) || globalPos[0] < eps_)
+ return Dumux::BoundaryConditions::dirichlet;
+ else
+ return Dumux::BoundaryConditions::neumann;
+ }
+
+ Scalar dirichletPress(const GlobalPosition& globalPos, const Intersection& intersection) const
+ {
+ return pLeftBc_;
+ }
+
+ Scalar dirichletSat(const GlobalPosition& globalPos, const Intersection& intersection) const
+ {
+ if (globalPos[0] < eps_)
+ return 0.8;
+
+ // all other boundaries
+
+ else
+ return 0.2;
+ }
+
+ std::vector<Scalar> neumannPress(const GlobalPosition& globalPos, const Intersection& intersection) const
+ {
+ std::vector<Scalar> neumannFlux(2, 0.0);
+ if (globalPos[0]> upperRight_[0] - eps_)
+ {
+ // the volume flux should remain constant, when density is changed
+ // here, we multiply by the density of the NonWetting Phase
+ const Scalar referenceDensity = 1000.0;
+ neumannFlux[nonwetting] = 3e-4 * densityNonWetting_/referenceDensity;
+ }
+ return neumannFlux;
+ }
+
+ Scalar neumannSat(const GlobalPosition& globalPos, const Intersection& intersection, Scalar factor) const
+ {
+ if (globalPos[0] > upperRight_[0] - eps_)
+ return factor;
+ return 0;
+ }
+ Scalar initSat(const GlobalPosition& globalPos, const Element& element) const
+ {
+ if (globalPos[0] < eps_)
+ return 0.8;
+
+ else
+ return 0.2;
+ }
+
+private:
+ GlobalPosition lowerLeft_;
+ GlobalPosition upperRight_;
+ Scalar eps_;
+ Scalar pLeftBc_;
+ Scalar simulationNumber_;
+ Scalar densityNonWetting_;
+ BuckleyLeverettAnalytic<TypeTag> analyticSolution_;
+};
+}
+#endif
diff --git a/appl/lecture/msm/buckleyleverett/external_interface.hh b/appl/lecture/msm/buckleyleverett/external_interface.hh
new file mode 100644
index 0000000000..24f1ee48d5
--- /dev/null
+++ b/appl/lecture/msm/buckleyleverett/external_interface.hh
@@ -0,0 +1,451 @@
+// $Id: pcm_parameters.hh 1329 $
+#ifndef PCM_PARAMETERS_HH
+#define PCM_PARAMETERS_HH
+
+#include <stdlib.h>
+#include <iostream>
+#include <fstream>
+
+/** \todo Please doc me! */
+namespace Dumux
+{
+
+class InterfaceSoilProperties
+//Class for Interface Soil Properties
+//Integrate Parameters for Probabilistic Collocation Method (iPCM).
+//Contained design and uncertainty parameters
+{
+public:
+ //Interface Soil Properties (ISP):
+ float ISP_Permeability; // Permeability
+ float ISP_FinePermeability; // Fine Permeability
+ float ISP_CoarsePermeability; // Coarse Permeability
+ float ISP_Porosity; // Porosity
+ float ISP_FinePorosity; // Fine Porosity
+ float ISP_CoarsePorosity; // Coarse Porosity
+ float ISP_LeakageWellPermeability; // Leakage Well Permeability
+ float ISP_LongitudinalDispersivity; //Longitudinal dispersivity
+ float ISP_TransverseDispersivity; //Transverse dispersivity
+ float ISP_BrooksCoreyLambda;
+ float ISP_BrooksCoreyEntryPressure;
+ float ISP_FineBrooksCoreyLambda;
+ float ISP_FineBrooksCoreyEntryPressure;
+ float ISP_CoarseBrooksCoreyLambda;
+ float ISP_CoarseBrooksCoreyEntryPressure;
+ float ISP_ResidualSaturationWetting;
+ float ISP_ResidualSaturationNonWetting;
+ float ISP_FineResidualSaturationWetting;
+ float ISP_FineResidualSaturationNonWetting;
+ float ISP_CoarseResidualSaturationWetting;
+ float ISP_CoarseResidualSaturationNonWetting;
+
+ InterfaceSoilProperties(const char* isp_filename)
+ //Initialization of ISP Parameters
+ {
+ using namespace std;
+ std::cout
+ << "-----> 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>"))
+ //cout << "-----> ISP: Soil Properties reading ... \n";
+ //ISP parameters initialization:
+ if (reader == string("<Permeability>"))
+ {
+ input >> reader;
+ ISP_Permeability = atof(reader);
+ cout << "-----> ISP: Permeability: " << ISP_Permeability
+ << "\n";
+ }
+ if (reader == string("<FinePermeability>"))
+ {
+ input >> reader;
+ ISP_FinePermeability = atof(reader);
+ cout << "-----> ISP: Fine permeability: "
+ << ISP_FinePermeability << "\n";
+ }
+ if (reader == string("<CoarsePermeability>"))
+ {
+ input >> reader;
+ ISP_CoarsePermeability = atof(reader);
+ cout << "-----> ISP: Coarse permeability: "
+ << ISP_CoarsePermeability << "\n";
+ }
+ if (reader == string("<Porosity>"))
+ {
+ input >> reader;
+ ISP_Porosity = atof(reader);
+ cout << "-----> ISP: Porosity: " << ISP_Porosity << "\n";
+ }
+ if (reader == string("<FinePorosity>"))
+ {
+ input >> reader;
+ ISP_FinePorosity = atof(reader);
+ cout << "-----> ISP: Fine porosity: " << ISP_FinePorosity
+ << "\n";
+ }
+ if (reader == string("<CoarsePorosity>"))
+ {
+ input >> reader;
+ ISP_CoarsePorosity = atof(reader);
+ cout << "-----> ISP: Coarse porosity: " << ISP_CoarsePorosity
+ << "\n";
+ }
+ if (reader == string("<LeakageWellPermeability>"))
+ {
+ input >> reader;
+ ISP_LeakageWellPermeability = atof(reader);
+ cout << "-----> ISP: Leakage Well Permeability: "
+ << ISP_LeakageWellPermeability << "\n";
+ }
+ if (reader == string("<LongitudinalDispersivity>"))
+ {
+ input >> reader;
+ ISP_LongitudinalDispersivity = atof(reader);
+ cout << "-----> ISP: Longitudinal dispersivity: "
+ << ISP_LongitudinalDispersivity << "\n";
+ }
+ if (reader == string("<TransverseDispersivity>"))
+ {
+ input >> reader;
+ ISP_TransverseDispersivity = atof(reader);
+ cout << "-----> ISP: Transverse dispersivity: "
+ << ISP_TransverseDispersivity << "\n";
+ }
+ if (reader == string("<BrooksCoreyLambda>"))
+ {
+ input >> reader;
+ ISP_BrooksCoreyLambda = atof(reader);
+ cout << "-----> ISP: Brooks-Corey lambda: "
+ << ISP_BrooksCoreyLambda << "\n";
+ }
+
+ if (reader == string("<FineBrooksCoreyLambda>"))
+ {
+ input >> reader;
+ ISP_FineBrooksCoreyLambda = atof(reader);
+ cout << "-----> ISP: Brooks-Corey lambda, fine: "
+ << ISP_FineBrooksCoreyLambda << "\n";
+ }
+ if (reader == string("<FineBrooksCoreyEntryPressure>"))
+ {
+ input >> reader;
+ ISP_FineBrooksCoreyEntryPressure = atof(reader);
+ cout << "-----> ISP: Brooks-Corey entry pressure, fine: "
+ << ISP_FineBrooksCoreyEntryPressure << "\n";
+ }
+ if (reader == string("<CoarseBrooksCoreyLambda>"))
+ {
+ input >> reader;
+ ISP_CoarseBrooksCoreyLambda = atof(reader);
+ cout << "-----> ISP: Brooks-Corey lambda, coarse: "
+ << ISP_CoarseBrooksCoreyLambda << "\n";
+ }
+
+ if (reader == string("<BrooksCoreyEntryPressure>"))
+ {
+ input >> reader;
+ ISP_BrooksCoreyEntryPressure = atof(reader);
+ cout << "-----> ISP: Brooks-Corey entry pressure: "
+ << ISP_BrooksCoreyEntryPressure << "\n";
+ }
+ if (reader == string("<CoarseBrooksCoreyEntryPressure>"))
+ {
+ input >> reader;
+ ISP_CoarseBrooksCoreyEntryPressure = atof(reader);
+ cout << "-----> ISP: Brooks-Corey entry pressure, coarse: "
+ << ISP_CoarseBrooksCoreyEntryPressure << "\n";
+ }
+ if (reader == string("<ResidualSaturationWetting>"))
+ {
+ input >> reader;
+ ISP_ResidualSaturationWetting = atof(reader);
+ cout << "-----> ISP: Residual saturation wetting phase: "
+ << ISP_ResidualSaturationWetting << "\n";
+ }
+ if (reader == string("<ResidualSaturationNonWetting>"))
+ {
+ input >> reader;
+ ISP_ResidualSaturationNonWetting = atof(reader);
+ cout << "-----> ISP: Residual saturation nonwetting phase: "
+ << ISP_ResidualSaturationNonWetting << "\n";
+ }
+ if (reader == string("<FineResidualSaturationWetting>"))
+ {
+ input >> reader;
+ ISP_FineResidualSaturationWetting = atof(reader);
+ cout << "-----> ISP: Residual saturation wetting phase, fine: "
+ << ISP_FineResidualSaturationWetting << "\n";
+ }
+ if (reader == string("<FineResidualSaturationNonWetting>"))
+ {
+ input >> reader;
+ ISP_FineResidualSaturationNonWetting = atof(reader);
+ cout << "-----> ISP: Residual saturation nonwetting phase, fine: "
+ << ISP_FineResidualSaturationNonWetting << "\n";
+ }
+ if (reader == string("<CoarseResidualSaturationWetting>"))
+ {
+ input >> reader;
+ ISP_CoarseResidualSaturationWetting = atof(reader);
+ cout << "-----> ISP: Residual saturation wetting phase, coarse: "
+ << ISP_CoarseResidualSaturationWetting << "\n";
+ }
+ if (reader == string("<CoarseResidualSaturationNonWetting>"))
+ {
+ input >> reader;
+ ISP_CoarseResidualSaturationNonWetting = atof(reader);
+ cout << "-----> ISP: Residual saturation nonwetting phase, coarse: "
+ << ISP_CoarseResidualSaturationNonWetting << "\n";
+ }
+ }
+ input.close();
+ }
+
+};
+
+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 IFP_GasDiffCoeff;// Gas Diffusion Coefficient
+ float IFP_CO2ResidSat;// Residual Saturation of CO2
+ float IFP_MolecularDiffusionCoefficient;
+ float IFP_ViscosityWettingFluid;
+ float IFP_ViscosityNonWettingFluid;
+ float IFP_DensityWettingFluid;
+ float IFP_DensityNonWettingFluid;
+
+ InterfaceFluidProperties(const char* ifp_filename)
+ //Initialization of IFP Parameters
+ {
+ using namespace std;
+ std::cout
+ << "-----> 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("<GasDiffusionCoeff>"))
+ {
+ input >> reader;
+ IFP_GasDiffCoeff = atof(reader);
+ cout << "-----> IFP: Gas Diffusion Coefficient: "
+ << IFP_GasDiffCoeff << "\n";
+ }
+ if (reader == string("<CO2ResidualSaturation>"))
+ {
+ input >> reader;
+ IFP_CO2ResidSat = atof(reader);
+ cout << "-----> IFP: Residual Saturation of CO2: "
+ << IFP_CO2ResidSat << "\n";
+ }
+ if (reader == string("<MolecularDiffusionCoefficient>"))
+ {
+ input >> reader;
+ IFP_MolecularDiffusionCoefficient = atof(reader);
+ cout << "-----> IFP: Molecular diffusion coefficient: "
+ << IFP_MolecularDiffusionCoefficient << "\n";
+ }
+ if (reader == string("<ViscosityWettingFluid>"))
+ {
+ input >> reader;
+ IFP_ViscosityWettingFluid = atof(reader);
+ cout << "-----> IFP: Viscosity of the wetting phase fluid: "
+ << IFP_ViscosityWettingFluid << "\n";
+ }
+ if (reader == string("<ViscosityNonWettingFluid>"))
+ {
+ input >> reader;
+ IFP_ViscosityNonWettingFluid = atof(reader);
+ cout << "-----> IFP: Viscosity of the non wetting phase fluid: "
+ << IFP_ViscosityNonWettingFluid << "\n";
+ }
+ if (reader == string("<DensityWettingFluid>"))
+ {
+ input >> reader;
+ IFP_DensityWettingFluid = atof(reader);
+ cout << "-----> IFP: Density of the wetting phase fluid: "
+ << IFP_DensityWettingFluid << "\n";
+ }
+ if (reader == string("<DensityNonWettingFluid>"))
+ {
+ input >> reader;
+ IFP_DensityNonWettingFluid = atof(reader);
+ cout << "-----> IFP: Density of the non wetting phase fluid: "
+ << IFP_DensityNonWettingFluid << "\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):
+ float IPP_DepthBOR; // Depth BOR
+ float IPP_InjectionWellRate; // Injection Well Rate
+ float IPP_InjectionWindowSize; // Injection Well Window Size
+ float IPP_UpperPressure; // Pressure at a top boundary
+ float IPP_LowerPressure; // Pressure at a lower boundary
+ float IPP_InfiltrationRate; // A Infiltration rate
+ float IPP_MaxTimeStepSize; // Maximum time step size
+ float IPP_InfiltrationStartTime; // time to stop an infiltration
+ float IPP_InfiltrationEndTime; // time to stop an infiltration
+ float IPP_DiscretizationLength;
+ float IPP_SimulationNumber; // possibility to enumerate the model run outputs
+
+ InterfaceProblemProperties(const char* ipp_filename)
+ //Initialization of IPP Parameters
+ {
+ using namespace std;
+ std::cout
+ << "-----> IPP: Interface Problem Properties Initialization ...\n";
+ //IPP input file defenition
+ 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("<DepthBOR>"))
+ {
+ input >> reader;
+ IPP_DepthBOR = atof(reader);
+ cout << "-----> IPP: Depth BOR: " << IPP_DepthBOR << "\n";
+ }
+ if (reader == string("<InjectionWellRate>"))
+ {
+ input >> reader;
+ IPP_InjectionWellRate = atof(reader);
+ cout << "-----> IPP: Injection Well Rate: "
+ << IPP_InjectionWellRate << "\n";
+ }
+ if (reader == string("<InjectionWellWindowSize>"))
+ {
+ input >> reader;
+ IPP_InjectionWindowSize = atof(reader);
+ cout << "-----> IPP: Injection Well Window Size: "
+ << IPP_InjectionWindowSize << "\n";
+ }
+ if (reader == string("<UpperPressure>"))
+ {
+ input >> reader;
+ IPP_UpperPressure = atof(reader);
+ cout << "-----> IPP: Upper pressure: "
+ << IPP_UpperPressure << "\n";
+ }
+ if (reader == string("<LowerPressure>"))
+ {
+ input >> reader;
+ IPP_LowerPressure = atof(reader);
+ cout << "-----> IPP: Lower pressure: "
+ << IPP_LowerPressure << "\n";
+ }
+ if (reader == string("<InfiltrationRate>"))
+ {
+ input >> reader;
+ IPP_InfiltrationRate = atof(reader);
+ cout << "-----> IPP: Infiltration rate: "
+ << IPP_InfiltrationRate << "\n";
+ }
+ if (reader == string("<MaxTimeStepSize>"))
+ {
+ input >> reader;
+ IPP_MaxTimeStepSize = atof(reader);
+ cout << "-----> IPP: Maximum time step size: "
+ << IPP_MaxTimeStepSize << "\n";
+ }
+ if (reader == string("<InfiltrationStartTime>"))
+ {
+ input >> reader;
+ IPP_InfiltrationStartTime = atof(reader);
+ cout << "-----> IPP: Start time of infiltration: "
+ << IPP_InfiltrationStartTime << "\n";
+ }
+ if (reader == string("<InfiltrationEndTime>"))
+ {
+ input >> reader;
+ IPP_InfiltrationEndTime = atof(reader);
+ cout << "-----> IPP: End time of infiltration: "
+ << IPP_InfiltrationEndTime << "\n";
+ }
+ if (reader == string("<DiscretizationLength>"))
+ {
+ input >> reader;
+ IPP_DiscretizationLength = atof(reader);
+ cout << "-----> IPP: Discretization length: "
+ << IPP_DiscretizationLength << "\n";
+ }
+ if (reader == string("<SimulationNumber>"))
+ {
+ input >> reader;
+ IPP_SimulationNumber = atof(reader);
+ cout << "-----> IPP: Output Name: "
+ << IPP_SimulationNumber << "\n";
+ }
+ }
+ input.close();
+ }
+
+};
+
+} // end namespace
+#endif
diff --git a/appl/lecture/msm/buckleyleverett/interface_BL.xml b/appl/lecture/msm/buckleyleverett/interface_BL.xml
new file mode 100755
index 0000000000..1d02f69c2f
--- /dev/null
+++ b/appl/lecture/msm/buckleyleverett/interface_BL.xml
@@ -0,0 +1,48 @@
+<iPCM_problem>
+
+ <SoilProperties>
+ <Permeability>
+ 1e-7
+ </Permeability>
+ <Porosity>
+ 0.2
+ </Porosity>
+ <BrooksCoreyLambda>
+ 2.0
+ </BrooksCoreyLambda>
+ <BrooksCoreyEntryPressure>
+ 0.0
+ </BrooksCoreyEntryPressure>
+ <ResidualSaturationWetting>
+ 0.2
+ </ResidualSaturationWetting>
+ <ResidualSaturationNonWetting>
+ 0.2
+ </ResidualSaturationNonWetting>
+ </SoilProperties>
+
+ <FluidProperties>
+ <DensityWettingFluid>
+ 1e3
+ </DensityWettingFluid>
+ <DensityNonWettingFluid>
+ 1e3
+ </DensityNonWettingFluid>
+ <ViscosityWettingFluid>
+ 1e-3
+ </ViscosityWettingFluid>
+ <ViscosityNonWettingFluid>
+ 1e-3
+ </ViscosityNonWettingFluid>
+ </FluidProperties>
+
+ <BoundaryAndInitialConditions>
+ <DiscretizationLength>
+ 10
+ </DiscretizationLength>
+ <SimulationNumber>
+ 1
+ </SimulationNumber>
+ </BoundaryAndInitialConditions>
+
+</iPCM_problem>
diff --git a/appl/lecture/msm/buckleyleverett/pseudoh2o.hh b/appl/lecture/msm/buckleyleverett/pseudoh2o.hh
new file mode 100644
index 0000000000..8cbde7c6f4
--- /dev/null
+++ b/appl/lecture/msm/buckleyleverett/pseudoh2o.hh
@@ -0,0 +1,80 @@
+/*****************************************************************************
+ * 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, as long as this copyright notice *
+ * is included in its original form. *
+ * *
+ * This program is distributed WITHOUT ANY WARRANTY. *
+ *****************************************************************************/
+/*!
+ * \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 oil.
+ */
+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.01;
+template <class ScalarT>
+typename PseudoH2O<ScalarT>::Scalar PseudoH2O<ScalarT>::density_ = 1000;
+} // end namepace
+
+#endif
diff --git a/appl/lecture/msm/buckleyleverett/pseudooil.hh b/appl/lecture/msm/buckleyleverett/pseudooil.hh
new file mode 100644
index 0000000000..fa2e4b2652
--- /dev/null
+++ b/appl/lecture/msm/buckleyleverett/pseudooil.hh
@@ -0,0 +1,83 @@
+/*****************************************************************************
+ * 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, as long as this copyright notice *
+ * is included in its original form. *
+ * *
+ * This program is distributed WITHOUT ANY WARRANTY. *
+ *****************************************************************************/
+/*!
+ * \file
+ *
+ * \brief Properties of pure water \f$H_2O\f$.
+ */
+#ifndef DUMUX_PSEUDOOIL_HH
+#define DUMUX_PSEUDOOIL_HH
+
+
+#include <dumux/material/components/component.hh>
+//#include "interface_BL.xml"
+
+namespace Dumux
+{
+/*!
+ * \brief Rough estimate for testing purposes of some oil.
+ */
+template <class ScalarT>
+class PseudoOil : public Component<ScalarT, PseudoOil<ScalarT> >
+{
+ typedef Component<ScalarT, PseudoOil<ScalarT> > ParentType;
+ typedef ScalarT Scalar;
+
+public:
+ /*!
+ * \brief A human readable name for the water.
+ */
+ static const char *name()
+ { return "Oil"; }
+
+ /*!
+ * \brief Rough estimate of the density of oil [kg/m^3].
+ */
+ static Scalar liquidDensity(Scalar temperature, Scalar pressure)
+ {
+ return density_;
+ }
+ // double densityOil = 1000;
+ // double viscosityWater = interfaceFluidProps.IFP_ViscosityWettingFluid;
+ // double viscosityOil = interfaceFluidProps.IFP_ViscosityNonWettingFluid;
+
+ /*!
+ * \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 PseudoOil<ScalarT>::Scalar PseudoOil<ScalarT>::viscosity_ = 0.01;
+template <class ScalarT>
+typename PseudoOil<ScalarT>::Scalar PseudoOil<ScalarT>::density_ = 1000;
+} // end namepace
+
+#endif
diff --git a/appl/lecture/msm/mcwhorter/Makefile.am b/appl/lecture/msm/mcwhorter/Makefile.am
new file mode 100644
index 0000000000..d05a200316
--- /dev/null
+++ b/appl/lecture/msm/mcwhorter/Makefile.am
@@ -0,0 +1,7 @@
+dist_noinst_DATA = mcwhorter_ff
+
+noinst_PROGRAMS = mcwhorter_ff
+
+mcwhorter_ff_SOURCES = mcwhorter_ff.cc
+
+include $(top_srcdir)/am/global-rules
diff --git a/appl/lecture/msm/mcwhorter/external_interface.hh b/appl/lecture/msm/mcwhorter/external_interface.hh
new file mode 100644
index 0000000000..0fd3bd7be2
--- /dev/null
+++ b/appl/lecture/msm/mcwhorter/external_interface.hh
@@ -0,0 +1,443 @@
+// $Id: pcm_parameters.hh 1329 $
+#ifndef PCM_PARAMETERS_HH
+#define PCM_PARAMETERS_HH
+
+#include <stdlib.h>
+#include <iostream>
+#include <fstream>
+
+/** \todo Please doc me! */
+namespace Dumux
+{
+
+class InterfaceSoilProperties
+//Class for Interface Soil Properties
+//Integrate Parameters for Probabilistic Collocation Method (iPCM).
+//Contained design and uncertainty parameters
+{
+public:
+ //Interface Soil Properties (ISP):
+ float ISP_Permeability; // Permeability
+ float ISP_FinePermeability; // Fine Permeability
+ float ISP_CoarsePermeability; // Coarse Permeability
+ float ISP_Porosity; // Porosity
+ float ISP_FinePorosity; // Fine Porosity
+ float ISP_CoarsePorosity; // Coarse Porosity
+ float ISP_LeakageWellPermeability; // Leakage Well Permeability
+ float ISP_LongitudinalDispersivity; //Longitudinal dispersivity
+ float ISP_TransverseDispersivity; //Transverse dispersivity
+ float ISP_BrooksCoreyLambda;
+ float ISP_BrooksCoreyEntryPressure;
+ float ISP_FineBrooksCoreyLambda;
+ float ISP_FineBrooksCoreyEntryPressure;
+ float ISP_CoarseBrooksCoreyLambda;
+ float ISP_CoarseBrooksCoreyEntryPressure;
+ float ISP_MinCapillaryPressure;
+ float ISP_MaxCapillaryPressure;
+ float ISP_ResidualSaturationWetting;
+ float ISP_ResidualSaturationNonWetting;
+ float ISP_FineResidualSaturationWetting;
+ float ISP_FineResidualSaturationNonWetting;
+ float ISP_CoarseResidualSaturationWetting;
+ float ISP_CoarseResidualSaturationNonWetting;
+
+ InterfaceSoilProperties(const char* isp_filename)
+ //Initialization of ISP Parameters
+ {
+ using namespace std;
+ std::cout
+ << "-----> 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>"))
+ //cout << "-----> ISP: Soil Properties reading ... \n";
+ //ISP perameters initialization:
+ if (reader == string("<Permeability>"))
+ {
+ input >> reader;
+ ISP_Permeability = atof(reader);
+ cout << "-----> ISP: Permeability: " << ISP_Permeability
+ << "\n";
+ }
+ if (reader == string("<FinePermeability>"))
+ {
+ input >> reader;
+ ISP_FinePermeability = atof(reader);
+ cout << "-----> ISP: Fine permeability: "
+ << ISP_FinePermeability << "\n";
+ }
+ if (reader == string("<CoarsePermeability>"))
+ {
+ input >> reader;
+ ISP_CoarsePermeability = atof(reader);
+ cout << "-----> ISP: Coarse permeability: "
+ << ISP_CoarsePermeability << "\n";
+ }
+ if (reader == string("<Porosity>"))
+ {
+ input >> reader;
+ ISP_Porosity = atof(reader);
+ cout << "-----> ISP: Porosity: " << ISP_Porosity << "\n";
+ }
+ if (reader == string("<FinePorosity>"))
+ {
+ input >> reader;
+ ISP_FinePorosity = atof(reader);
+ cout << "-----> ISP: Fine porosity: " << ISP_FinePorosity
+ << "\n";
+ }
+ if (reader == string("<CoarsePorosity>"))
+ {
+ input >> reader;
+ ISP_CoarsePorosity = atof(reader);
+ cout << "-----> ISP: Coarse porosity: " << ISP_CoarsePorosity
+ << "\n";
+ }
+ if (reader == string("<LeakageWellPermeability>"))
+ {
+ input >> reader;
+ ISP_LeakageWellPermeability = atof(reader);
+ cout << "-----> ISP: Leakage Well Permeability: "
+ << ISP_LeakageWellPermeability << "\n";
+ }
+ if (reader == string("<LongitudinalDispersivity>"))
+ {
+ input >> reader;
+ ISP_LongitudinalDispersivity = atof(reader);
+ cout << "-----> ISP: Longitudinal dispersivity: "
+ << ISP_LongitudinalDispersivity << "\n";
+ }
+ if (reader == string("<TransverseDispersivity>"))
+ {
+ input >> reader;
+ ISP_TransverseDispersivity = atof(reader);
+ cout << "-----> ISP: Transverse dispersivity: "
+ << ISP_TransverseDispersivity << "\n";
+ }
+ if (reader == string("<BrooksCoreyLambda>"))
+ {
+ input >> reader;
+ ISP_BrooksCoreyLambda = atof(reader);
+ cout << "-----> ISP: Brooks-Corey lambda: "
+ << ISP_BrooksCoreyLambda << "\n";
+ }
+ if (reader == string("<BrooksCoreyEntryPressure>"))
+ {
+ input >> reader;
+ ISP_BrooksCoreyEntryPressure = atof(reader);
+ cout << "-----> ISP: Brooks-Corey entry pressure: "
+ << ISP_BrooksCoreyEntryPressure << "\n";
+ }
+
+ if (reader == string("<FineBrooksCoreyLambda>"))
+ {
+ input >> reader;
+ ISP_FineBrooksCoreyLambda = atof(reader);
+ cout << "-----> ISP: Brooks-Corey lambda, fine: "
+ << ISP_FineBrooksCoreyLambda << "\n";
+ }
+ if (reader == string("<FineBrooksCoreyEntryPressure>"))
+ {
+ input >> reader;
+ ISP_FineBrooksCoreyEntryPressure = atof(reader);
+ cout << "-----> ISP: Brooks-Corey entry pressure, fine: "
+ << ISP_FineBrooksCoreyEntryPressure << "\n";
+ }
+ if (reader == string("<CoarseBrooksCoreyLambda>"))
+ {
+ input >> reader;
+ ISP_CoarseBrooksCoreyLambda = atof(reader);
+ cout << "-----> ISP: Brooks-Corey lambda, coarse: "
+ << ISP_CoarseBrooksCoreyLambda << "\n";
+ }
+ if (reader == string("<CoarseBrooksCoreyEntryPressure>"))
+ {
+ input >> reader;
+ ISP_CoarseBrooksCoreyEntryPressure = atof(reader);
+ cout << "-----> ISP: Brooks-Corey entry pressure, coarse: "
+ << ISP_CoarseBrooksCoreyEntryPressure << "\n";
+ }
+ if (reader == string("<MinCapillaryPressure>"))
+ {
+ input >> reader;
+ ISP_MinCapillaryPressure = atof(reader);
+ cout << "-----> ISP: Minimal capillary pressure: "
+ << ISP_MinCapillaryPressure << "\n";
+ }
+ if (reader == string("<MaxCapillaryPressure>"))
+ {
+ input >> reader;
+ ISP_MaxCapillaryPressure = atof(reader);
+ cout << "-----> ISP: Maximal capillary pressure: "
+ << ISP_MaxCapillaryPressure << "\n";
+ }
+ if (reader == string("<ResidualSaturationWetting>"))
+ {
+ input >> reader;
+ ISP_ResidualSaturationWetting = atof(reader);
+ cout << "-----> ISP: Residual saturation wetting phase: "
+ << ISP_ResidualSaturationWetting << "\n";
+ }
+ if (reader == string("<ResidualSaturationNonWetting>"))
+ {
+ input >> reader;
+ ISP_ResidualSaturationNonWetting = atof(reader);
+ cout << "-----> ISP: Residual saturation nonwetting phase: "
+ << ISP_ResidualSaturationNonWetting << "\n";
+ }
+ if (reader == string("<FineResidualSaturationWetting>"))
+ {
+ input >> reader;
+ ISP_FineResidualSaturationWetting = atof(reader);
+ cout << "-----> ISP: Residual saturation wetting phase, fine: "
+ << ISP_FineResidualSaturationWetting << "\n";
+ }
+ if (reader == string("<FineResidualSaturationNonWetting>"))
+ {
+ input >> reader;
+ ISP_FineResidualSaturationNonWetting = atof(reader);
+ cout << "-----> ISP: Residual saturation nonwetting phase, fine: "
+ << ISP_FineResidualSaturationNonWetting << "\n";
+ }
+ if (reader == string("<CoarseResidualSaturationWetting>"))
+ {
+ input >> reader;
+ ISP_CoarseResidualSaturationWetting = atof(reader);
+ cout << "-----> ISP: Residual saturation wetting phase, coarse: "
+ << ISP_CoarseResidualSaturationWetting << "\n";
+ }
+ if (reader == string("<CoarseResidualSaturationNonWetting>"))
+ {
+ input >> reader;
+ ISP_CoarseResidualSaturationNonWetting = atof(reader);
+ cout << "-----> ISP: Residual saturation nonwetting phase, coarse: "
+ << ISP_CoarseResidualSaturationNonWetting << "\n";
+ }
+ }
+ input.close();
+ }
+
+};
+
+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 IFP_GasDiffCoeff;// Gas Diffusion Coefficient
+ float IFP_CO2ResidSat;// Residual Saturation of CO2
+ float IFP_MolecularDiffusionCoefficient;
+ float IFP_ViscosityWettingFluid;
+ float IFP_ViscosityNonWettingFluid;
+
+ InterfaceFluidProperties(const char* ifp_filename)
+ //Initialization of IFP Parameters
+ {
+ using namespace std;
+ std::cout
+ << "-----> 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("<GasDiffusionCoeff>"))
+ {
+ input >> reader;
+ IFP_GasDiffCoeff = atof(reader);
+ cout << "-----> IFP: Gas Diffusion Coefficient: "
+ << IFP_GasDiffCoeff << "\n";
+ }
+ if (reader == string("<CO2ResidualSaturation>"))
+ {
+ input >> reader;
+ IFP_CO2ResidSat = atof(reader);
+ cout << "-----> IFP: Residual Saturation of CO2: "
+ << IFP_CO2ResidSat << "\n";
+ }
+ if (reader == string("<MolecularDiffusionCoefficient>"))
+ {
+ input >> reader;
+ IFP_MolecularDiffusionCoefficient = atof(reader);
+ cout << "-----> IFP: Molecular diffusion coefficient: "
+ << IFP_MolecularDiffusionCoefficient << "\n";
+ }
+ if (reader == string("<ViscosityWettingFluid>"))
+ {
+ input >> reader;
+ IFP_ViscosityWettingFluid = atof(reader);
+ cout << "-----> IFP: Viscosity of the wetting phase fluid: "
+ << IFP_ViscosityWettingFluid << "\n";
+ }
+ if (reader == string("<ViscosityNonWettingFluid>"))
+ {
+ input >> reader;
+ IFP_ViscosityNonWettingFluid = atof(reader);
+ cout << "-----> IFP: Viscosity of the non wetting phase fluid: "
+ << IFP_ViscosityNonWettingFluid << "\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):
+ float IPP_DepthBOR; // Depth BOR
+ float IPP_InjectionWellRate; // Injection Well Rate
+ float IPP_InjectionWindowSize; // Injection Well Window Size
+ float IPP_UpperPressure; // Pressure at a top boundary
+ float IPP_LowerPressure; // Pressure at a lower boundary
+ float IPP_InfiltrationRate; // A Infiltration rate
+ float IPP_MaxTimeStepSize; // Maximum time step size
+ float IPP_InfiltrationStartTime; // time to stop an infiltration
+ float IPP_InfiltrationEndTime; // time to stop an infiltration
+ float IPP_DiscretizationLength;
+
+ InterfaceProblemProperties(const char* ipp_filename)
+ //Initialization of IPP Parameters
+ {
+ using namespace std;
+ std::cout
+ << "-----> IPP: Interface Soil Properties Initialization ...\n";
+ //IPP input file defenition
+ 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("<DepthBOR>"))
+ {
+ input >> reader;
+ IPP_DepthBOR = atof(reader);
+ cout << "-----> IPP: Depth BOR: " << IPP_DepthBOR << "\n";
+ }
+ if (reader == string("<InjectionWellRate>"))
+ {
+ input >> reader;
+ IPP_InjectionWellRate = atof(reader);
+ cout << "-----> IPP: Injection Well Rate: "
+ << IPP_InjectionWellRate << "\n";
+ }
+ if (reader == string("<InjectionWellWindowSize>"))
+ {
+ input >> reader;
+ IPP_InjectionWindowSize = atof(reader);
+ cout << "-----> IPP: Injection Well Window Size: "
+ << IPP_InjectionWindowSize << "\n";
+ }
+ if (reader == string("<UpperPressure>"))
+ {
+ input >> reader;
+ IPP_UpperPressure = atof(reader);
+ cout << "-----> IPP: Upper pressure: "
+ << IPP_UpperPressure << "\n";
+ }
+ if (reader == string("<LowerPressure>"))
+ {
+ input >> reader;
+ IPP_LowerPressure = atof(reader);
+ cout << "-----> IPP: Lower pressure: "
+ << IPP_LowerPressure << "\n";
+ }
+ if (reader == string("<InfiltrationRate>"))
+ {
+ input >> reader;
+ IPP_InfiltrationRate = atof(reader);
+ cout << "-----> IPP: Infiltration rate: "
+ << IPP_InfiltrationRate << "\n";
+ }
+ if (reader == string("<MaxTimeStepSize>"))
+ {
+ input >> reader;
+ IPP_MaxTimeStepSize = atof(reader);
+ cout << "-----> IPP: Maximum time step size: "
+ << IPP_MaxTimeStepSize << "\n";
+ }
+ if (reader == string("<InfiltrationStartTime>"))
+ {
+ input >> reader;
+ IPP_InfiltrationStartTime = atof(reader);
+ cout << "-----> IPP: Start time of infiltration: "
+ << IPP_InfiltrationStartTime << "\n";
+ }
+ if (reader == string("<InfiltrationEndTime>"))
+ {
+ input >> reader;
+ IPP_InfiltrationEndTime = atof(reader);
+ cout << "-----> IPP: End time of infiltration: "
+ << IPP_InfiltrationEndTime << "\n";
+ }
+ if (reader == string("<DiscretizationLength>"))
+ {
+ input >> reader;
+ IPP_DiscretizationLength = atof(reader);
+ cout << "-----> IPP: Discretization length: "
+ << IPP_DiscretizationLength << "\n";
+ }
+ }
+ input.close();
+ }
+
+};
+
+} // end namespace
+#endif
diff --git a/appl/lecture/msm/mcwhorter/interface_MW.xml b/appl/lecture/msm/mcwhorter/interface_MW.xml
new file mode 100755
index 0000000000..0cae89193d
--- /dev/null
+++ b/appl/lecture/msm/mcwhorter/interface_MW.xml
@@ -0,0 +1,51 @@
+<iPCM_problem>
+
+ <SoilProperties>
+ <Permeability>
+ 1e-10
+ </Permeability>
+ <Porosity>
+ 0.4
+ </Porosity>
+ <BrooksCoreyLambda>
+ 2.0
+ </BrooksCoreyLambda>
+ <BrooksCoreyEntryPressure>
+ 5000
+ </BrooksCoreyEntryPressure>
+ <MinCapillaryPressure>
+ 5000
+ </MinCapillaryPressure>
+ <MaxCapillaryPressure>
+ 8000
+ </MaxCapillaryPressure>
+ <ResidualSaturationWetting>
+ 0.
+ </ResidualSaturationWetting>
+ <ResidualSaturationNonWetting>
+ 0.
+ </ResidualSaturationNonWetting>
+ </SoilProperties>
+
+ <FluidProperties>
+ <FluidDensityWetting>
+ 1e3
+ </FluidDensityWetting>
+ <DensityNonWettingFluid>
+ 1e3
+ </DensityNonWettingFluid>
+ <ViscosityWettingFluid>
+ 1e-3
+ </ViscosityWettingFluid>
+ <ViscosityNonWettingFluid>
+ 1e-3
+ </ViscosityNonWettingFluid>
+ </FluidProperties>
+
+ <BoundaryAndInitialConditions>
+ <DiscretizationLength>
+ 0.05
+ </DiscretizationLength>
+ </BoundaryAndInitialConditions>
+
+</iPCM_problem>
diff --git a/appl/lecture/msm/mcwhorter/mcwhorter_analytic.hh b/appl/lecture/msm/mcwhorter/mcwhorter_analytic.hh
new file mode 100644
index 0000000000..876b9bd74f
--- /dev/null
+++ b/appl/lecture/msm/mcwhorter/mcwhorter_analytic.hh
@@ -0,0 +1,390 @@
+// $Id$
+
+#ifndef DUMUX_MCWHORTER_ANALYTIC_HH
+#define DUMUX_MCWHORTER_ANALYTIC_HH
+
+/**
+ * @file
+ * @brief Analytic solution of
+ * the McWhorter problem
+ * @author Markus Wolff, Anneli Schöniger
+ */
+
+namespace Dumux
+{
+/**
+ * @brief Analytic solution of
+ * the McWhorter problem
+ *
+ * for naming of variables see "An Improved Semi-Analytical Solution for Verification
+ * of Numerical Models of Two-Phase Flow in Porous Media"
+ * (R. Fucik, J. Mikyska, M. Benes, T. H. Illangasekare; 2007)
+ */
+
+template<class TypeTag>
+class McWhorterAnalytic
+{
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(Problem)) Problem;
+ 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 GET_PROP_TYPE(TypeTag, PTAG(SpatialParameters)) SpatialParameters;
+ typedef typename SpatialParameters::MaterialLaw MaterialLaw;
+ typedef typename MaterialLaw::Params MaterialLawParams;
+
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(FluidSystem)) FluidSystem;
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(FluidState)) FluidState;
+
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(TwoPIndices)) Indices;
+
+ enum
+ {
+ dim = GridView::dimension, dimworld = GridView::dimensionworld
+ };
+ enum
+ {
+ wPhaseIdx = Indices::wPhaseIdx, nPhaseIdx = Indices::nPhaseIdx
+ };
+
+ typedef Dune::BlockVector<Dune::FieldVector<Scalar, 1> > BlockVector;
+ typedef typename GridView::Traits::template Codim<0>::Entity Element;
+ typedef typename GridView::template Codim<0>::Iterator ElementIterator;
+ typedef Dune::FieldVector<Scalar, dimworld> GlobalPosition;
+
+public:
+
+ // functions needed for analytical solution
+
+ void initializeAnalytic()
+ {
+ analyticSolution_.resize(size_);
+ analyticSolution_=0;
+ error_.resize(size_);
+ error_=0;
+ elementVolume_.resize(size_);
+ elementVolume_=0;
+
+ return;
+ }
+
+ void calcSatError(BlockVector &Approx)
+ {
+ error_=0;
+ elementVolume_=0;
+
+ ElementIterator eItEnd = problem_.gridView().template end<0> ();
+
+ for (ElementIterator eIt = problem_.gridView().template end<0> (); eIt!= eItEnd; ++eIt)
+ {
+ // get entity
+ const Element& element = *eIt;
+ int index = problem_.variables().index(*eIt);
+ elementVolume_[index]= element.geometry().volume();
+ // std::cout<<"elementVolume_ = "<<elementVolume_[index]<<std::endl;
+ }
+
+ double globalVolume = elementVolume_.one_norm();
+ //std::cout<<"globalVolume = "<<globalVolume<<std::endl; TODO: globalVolume always zero!
+
+ for (int i=0; i<size_; i++)
+ {
+ error_[i]=analyticSolution_[i]-Approx[i];
+ //std::cout<<"error_ = "<<error_[i]<<std::endl;
+ //std::cout<<"analyticSolution_ = "<<analyticSolution_[i]<<std::endl;
+ //std::cout<<"Approx = "<<Approx[i]<<std::endl;
+ }
+
+ // std::cout<<"error_ = "<<error_<<std::endl;
+
+ double diffNorm = error_.two_norm();
+ // std::cout<<"diffNorm = "<<diffNorm<<std::endl;
+
+ for (int i=0; i<size_; i++)
+ {
+ error_[i] = (globalVolume) ? diffNorm * pow((elementVolume_[i]/globalVolume), 0.5) : 0;
+ }
+
+ return;
+ }
+
+ void prepareAnalytic()
+ {
+ const MaterialLawParams& materialLawParams(problem_.spatialParameters().materialLawParams(dummyGlobal_, dummyElement_));
+
+ swr_ = materialLawParams.Swr();
+ snr_ = materialLawParams.Snr();
+ porosity_ = problem_.spatialParameters().porosity(dummyGlobal_, dummyElement_);
+ permeability_ = problem_.spatialParameters().intrinsicPermeability(dummyGlobal_, dummyElement_)[0][0];
+ sInit_ = problem_.initSat(dummyGlobal_, dummyElement_);
+ Scalar s0 =(1 - snr_ - swr_);
+ time_=0;
+
+ h_= (s0 - sInit_)/intervalNum_;
+// std::cout<<"h_= "<<h_<<std::endl;
+
+ // define saturation range for analytic solution
+ satVec_= swr_;
+ for (int i=1; i<pointNum_; i++)
+ {
+ satVec_[i]=satVec_[i-1]+h_;
+ }
+ FluidState fluidState;
+ Scalar temp = problem_.temperature(dummyGlobal_, dummyElement_);
+ Scalar press = problem_.referencePressure(dummyGlobal_, dummyElement_);
+ Scalar viscosityW = FluidSystem::phaseViscosity(wPhaseIdx, temp, press, fluidState);
+ Scalar viscosityNW = FluidSystem::phaseViscosity(nPhaseIdx, temp, press, fluidState);
+
+
+ // get fractional flow function vector
+ for (int i=0; i<pointNum_; i++)
+ {
+ fractionalW_[i] = MaterialLaw::krw(materialLawParams, satVec_[i])/viscosityW;
+ fractionalW_[i] /= (fractionalW_[i] + MaterialLaw::krn(materialLawParams, satVec_[i])/viscosityNW);
+ }
+
+ // get capillary pressure derivatives
+ dpcdsw_=0;
+
+ for (int i=0; i<pointNum_; i++)
+ {
+ dpcdsw_[i] = MaterialLaw::dpC_dSw(materialLawParams, satVec_[i]);
+ }
+// std::cout<<"dpcdsw = "<<dpcdsw_<<std::endl;
+
+ // set initial fW
+ if (sInit_ == 0)
+ fInit_=0;
+ else
+ fInit_=fractionalW_[0];
+
+ fractionalW_[0]=0;
+
+ // normalize fW
+ // with r_ = qt/q0
+ // qt: total volume flux, q0: displacing phase flux at boundary
+ // --> r_ = 1 for unidirectional displacement; r_ = 0 for impermeable boundary
+ for (int i=0; i<pointNum_; i++)
+ {
+ fn_[i]= r_ * (fractionalW_[i] - fInit_)/ (1 - r_ * fInit_);
+ }
+
+ // std::cout<<"r_ = "<<r_<<std::endl;
+ // std::cout<<"fn_ = "<<fn_<<std::endl;
+
+ // diffusivity function
+ for (int i=0; i<pointNum_; i++)
+ {d_[i] = fractionalW_[i]*(-dpcdsw_[i])*(MaterialLaw::krn(materialLawParams, satVec_[i])/viscosityNW)*permeability_;
+ }
+
+ // std::cout<<"fractionalW_ = "<<fractionalW_<<std::endl;
+ // std::cout<<"permeability_ = "<<permeability_<<std::endl;
+ // std::cout<<"d_ = "<<d_<<std::endl;
+
+
+ // gk_: part of fractional flow function
+ // initial guess for gk_
+ for (int i=0; i<pointNum_; i++)
+ {
+ gk_[i] = d_[i]/(1-fn_[i]);
+ }
+
+ gk_[0] = 0;
+
+// std::cout<<"gk_ = "<<gk_<<std::endl;
+
+ return;
+ }
+
+ void updateExSol()
+ {
+ // with displacing phase flux at boundary q0 = A * 1/sqrt(t)
+ // Akm1, Ak: successive approximations of A
+ double Ak = 0;
+ double Akm1 = 0;
+ double diff = 1e100;
+
+ // partial numerical integrals a_, b_
+ a_=0, b_=0;
+ fp_=0;
+
+ // approximation of integral I
+ double I0 = 0;
+ double Ii = 0;
+
+ int k = 0;
+
+ while (diff> tolAnalytic_)
+ {
+ k++;
+// std::cout<<" k = "<<k<<std::endl;
+ if (k> 50000)
+ {
+ std::cout<<"Analytic solution: Too many iterations!"<<std::endl;
+ break;
+ }
+
+ Akm1=Ak;
+ I0=0;
+ for (int i=0; i<intervalNum_; i++)
+ {
+ a_[i] = 0.5 * h_ * sInit_ *(gk_[i] + gk_[i+1])+ pow(h_, 2) / 6* ((3* i + 1) * gk_[i]
+ + (3 * i + 2) * gk_[i+1]);
+ b_[i] = 0.5 * h_ * (gk_[i] + gk_[i+1]);
+ I0 += (a_[i] - sInit_ * b_[i]);
+ }
+ // std::cout<<" I0 = "<<I0<<std::endl;
+
+ gk_[0]=0;
+ for (int i=1; i<pointNum_; i++)
+ {
+ Ii=0;
+ for (int j = i; j<intervalNum_; j++)
+ Ii += (a_[j] - satVec_[i] * b_[j]);
+ //gk_[i] = d_[i] + gk_[i]*(fn_[i] + Ii/I0); // method A
+ gk_[i] = (d_[i] + gk_[i]*fn_[i])/(1 - Ii/I0); // method B
+ }
+
+ // with f(sInit) = 0: relationship between A and sInit
+ Ak = pow((0.5*porosity_/pow((1 - fInit_), 2)*I0), 0.5);
+ diff=fabs(Ak - Akm1);
+ // std::cout<<"diff = "<<diff<<std::endl;
+ }
+
+ // std::cout<<" b_ = "<<b_<<std::endl;
+ // std::cout<<" Ak = "<<Ak<<std::endl;
+
+
+ // fp_: first derivative of f
+ for (int i = 0; i<pointNum_; i++)
+ {
+ for (int j = i; j<intervalNum_; j++)
+ fp_[i] += b_[j]/I0;
+ }
+
+ // std::cout<<" fp_ = "<<fp_<<std::endl;
+ // std::cout<<fInit_<<std::endl;
+
+ for (int i = 0; i<pointNum_; i++)
+ {
+ xf_[i]= 2 * Ak * (1 - fInit_ * r_)/ porosity_ * fp_[i]* pow(time_, 0.5);
+ }
+
+// std::cout<<" xf_ = "<<xf_<<std::endl;
+
+ // iterate over vertices and get analytical saturation solution
+ ElementIterator eItEnd = problem_.gridView().template end<0> ();
+ for (ElementIterator eIt = problem_.gridView().template begin<0> (); eIt!= eItEnd; ++eIt)
+ {
+ // get global coordinate of cell center
+ const GlobalPosition& globalPos = eIt->geometry().center();
+
+// std::cout<<"globalPos = "<<globalPos[0]<<", x0 = "<<xf_[0]<<"\n";
+
+ // find index of current vertex
+ int index = problem_.variables().index(*eIt);
+
+ // find x_f next to global coordinate of the vertex
+ int xnext = 0;
+ for (int i=intervalNum_; i>0; i--)
+ {
+ if (globalPos[0] <= xf_[i])
+ {
+ xnext = i;
+ break;
+ }
+ }
+
+ // account for the area not yet reached by the front
+ if (globalPos[0]> xf_[0])
+ {
+ analyticSolution_[index] = sInit_;
+ continue;
+ }
+
+ if (globalPos[0] <= xf_[0])
+ {
+ analyticSolution_[index] = satVec_[xnext];
+ continue;
+ }
+
+// std::cout<<"Analytical = "<<satVec_[xnext]<<std::endl;
+ }
+
+ // call function to calculate the saturation error
+ calcSatError(problem_.variables().saturation());
+ return;
+ }
+
+ void calculateAnalyticSolution()
+ {
+ time_ += problem_.timeManager().timeStepSize();
+ updateExSol();
+ }
+
+ //Write saturation and pressure into file
+ template<class MultiWriter>
+ void addOutputVtkFields(MultiWriter &writer)
+ {
+ BlockVector *analyticSolution = writer.template createField<Scalar, 1> (size_);
+ BlockVector *error = writer.template createField<Scalar, 1> (size_);
+
+ *analyticSolution = analyticSolution_;
+ *error = error_;
+
+ writer.addCellData(analyticSolution, "saturation (exact solution)");
+ writer.addCellData(error, "error_");
+
+ return;
+ }
+
+ McWhorterAnalytic(Problem& problem, Scalar tolAnalytic = 1e-14) :
+ problem_(problem), analyticSolution_(0), error_(0), elementVolume_(0), size_(problem.gridView().size(0)), dummyElement_(
+ *(problem_.gridView().template begin<0> ())), dummyGlobal_(GlobalPosition(1)), tolAnalytic_(tolAnalytic)
+ {
+ initializeAnalytic();
+ prepareAnalytic();
+ }
+
+protected:
+
+
+private:
+ Problem& problem_;
+
+ BlockVector analyticSolution_;
+ BlockVector error_;
+ BlockVector elementVolume_;
+
+ int size_;
+
+ const Element& dummyElement_;
+ const GlobalPosition& dummyGlobal_;
+
+ Scalar tolAnalytic_;
+ Scalar r_;
+
+ Scalar swr_;
+ Scalar snr_;
+ Scalar porosity_;
+ Scalar sInit_;
+ Scalar time_;
+ Scalar permeability_;
+ enum
+ { intervalNum_ = 1000, pointNum_ = intervalNum_+1};
+ Dune::FieldVector<Scalar, pointNum_> satVec_;
+ Dune::FieldVector<Scalar,pointNum_> fractionalW_;
+ Dune::FieldVector<Scalar, pointNum_> dpcdsw_;
+ Dune::FieldVector<Scalar, pointNum_> fn_;
+ Dune::FieldVector<Scalar, pointNum_> d_;
+ Dune::FieldVector<Scalar, pointNum_> gk_;
+ Dune::FieldVector<Scalar, pointNum_> xf_;
+ Scalar fInit_;
+ Scalar h_;
+ Dune::FieldVector<Scalar,intervalNum_> a_;
+ Dune::FieldVector<Scalar,intervalNum_> b_;
+ Dune::FieldVector<Scalar,pointNum_> fp_;
+
+};
+}
+#endif
diff --git a/appl/lecture/msm/mcwhorter/mcwhorter_ff.cc b/appl/lecture/msm/mcwhorter/mcwhorter_ff.cc
new file mode 100644
index 0000000000..78a721c139
--- /dev/null
+++ b/appl/lecture/msm/mcwhorter/mcwhorter_ff.cc
@@ -0,0 +1,81 @@
+#include "config.h"
+#include <iostream>
+#include <boost/format.hpp>
+#include <iomanip>
+
+#include <dune/grid/common/gridinfo.hh>
+#include "external_interface.hh"
+#include "mcwhorterproblem.hh"
+
+
+int main(int argc, char** argv)
+{
+ try
+ {
+ typedef TTAG(McWhorterProblem) TypeTag;
+ typedef GET_PROP_TYPE(TypeTag, PTAG(Scalar)) Scalar;
+ typedef GET_PROP_TYPE(TypeTag, PTAG(Grid)) Grid;
+ typedef GET_PROP_TYPE(TypeTag, PTAG(Problem)) Problem;
+ typedef Dune::FieldVector<Scalar, Grid::dimensionworld> GlobalPosition;
+
+ static const int dim = Grid::dimension;
+
+ // initialize MPI, finalize is done automatically on exit
+ Dune::MPIHelper::instance(argc, argv);
+
+ //load interface-file
+ Dumux::InterfaceProblemProperties interfaceProbProps("interface_MW.xml");
+ double discretizationLength = interfaceProbProps.IPP_DiscretizationLength;
+ int cellNumber = 2/discretizationLength;
+
+ // define the problem dimensions
+ Dune::FieldVector<Scalar, dim> lowerLeft(0);
+ Dune::FieldVector<Scalar, dim> upperRight(2.0);
+// UpperRight[1] = 1;
+ Dune::FieldVector<int, dim> cellNumbers(cellNumber);
+// cellNumbers[0] = 26;
+ if (argc != 2)
+ {
+ std::cout << "usage: tEnd" << std::endl;
+ return 0;
+ }
+ std::string arg1(argv[1]);
+ std::istringstream is1(arg1);
+ double tEnd;
+ is1 >> tEnd;
+ double dt = tEnd;
+ // create a grid object
+ typedef Dune::SGrid<dim, dim> GridType;
+ typedef GridType::LeafGridView GridView;
+
+ // grid reference
+ GridType grid(cellNumbers,lowerLeft,upperRight);
+ Dune::gridinfo(grid);
+
+
+ ////////////////////////////////////////////////////////////
+ // instantiate and run the concrete problem
+ ////////////////////////////////////////////////////////////
+
+ Dumux::InterfaceFluidProperties interfaceFluidProps("interface_MW.xml");
+ typedef GET_PROP_TYPE(TypeTag, PTAG(WettingPhase)) WettingPhase;
+ typedef GET_PROP_TYPE(TypeTag, PTAG(NonwettingPhase)) NonwettingPhase;
+
+ WettingPhase::Component::setViscosity(interfaceFluidProps.IFP_ViscosityWettingFluid);
+ NonwettingPhase::Component::setViscosity(interfaceFluidProps.IFP_ViscosityNonWettingFluid);
+
+ Problem problem(grid.leafView(), lowerLeft, upperRight);
+
+ problem.timeManager().init(problem, 0, dt, tEnd);
+ problem.timeManager().run();
+
+ return 0;
+
+ } catch (Dune::Exception &e)
+ {
+ std::cerr << "Dune reported error: " << e << std::endl;
+ } catch (...)
+ {
+ std::cerr << "Unknown exception thrown!" << std::endl;
+ }
+}
diff --git a/appl/lecture/msm/mcwhorter/mcwhorter_spatialparams.hh b/appl/lecture/msm/mcwhorter/mcwhorter_spatialparams.hh
new file mode 100644
index 0000000000..868caa590d
--- /dev/null
+++ b/appl/lecture/msm/mcwhorter/mcwhorter_spatialparams.hh
@@ -0,0 +1,112 @@
+// $Id: buckleyleverett_spatialparams.hh 4144 2010-08-24 10:10:47Z bernd $
+/*****************************************************************************
+ * Copyright (C) 2008-2009 by Markus Wolff *
+ * Copyright (C) 2010 by Klaus Mosthaf *
+ * 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, as long as this copyright notice *
+ * is included in its original form. *
+ * *
+ * This program is distributed WITHOUT ANY WARRANTY. *
+ *****************************************************************************/
+#ifndef MCWHORTER_SPATIALPARAMETERS_HH
+#define MCWHORTER_SPATIALPARAMETERS_HH
+
+
+#include <dumux/material/fluidmatrixinteractions/2p/linearmaterial.hh>
+#include <dumux/material/fluidmatrixinteractions/2p/regularizedbrookscorey.hh>
+#include <dumux/material/fluidmatrixinteractions/2p/efftoabslaw.hh>
+
+namespace Dumux
+{
+
+/** \todo Please doc me! */
+
+template<class TypeTag>
+class McWhorterSpatialParams
+{
+ 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 RegularizedBrooksCorey<Scalar> RawMaterialLaw;
+// 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
+ {
+ return constPermeability_;
+ }
+
+ double porosity(const GlobalPosition& globalPos, const Element& element) const
+ {
+ return porosity_;
+ }
+
+
+ // return the brooks-corey context depending on the position
+ const MaterialLawParams& materialLawParams(const GlobalPosition& globalPos, const Element &element) const
+ {
+ return materialLawParams_;
+ }
+
+
+ McWhorterSpatialParams(const GridView& gridView)
+ : constPermeability_(0)
+ {
+ Dumux::InterfaceSoilProperties interfaceSoilProps("interface_MW.xml");
+
+ // residual saturations
+ materialLawParams_.setSwr(interfaceSoilProps.ISP_ResidualSaturationWetting);
+ materialLawParams_.setSnr(interfaceSoilProps.ISP_ResidualSaturationNonWetting);
+
+ porosity_ = interfaceSoilProps.ISP_Porosity;
+
+// // parameters for the Brooks-Corey Law
+// // entry pressures
+ materialLawParams_.setPe(interfaceSoilProps.ISP_BrooksCoreyEntryPressure);
+// // Brooks-Corey shape parameters
+ materialLawParams_.setAlpha(interfaceSoilProps.ISP_BrooksCoreyLambda);
+
+ // parameters for the linear
+ // entry pressures function
+// materialLawParams_.setEntryPC(0);
+// materialLawParams_.setMaxPC(0);
+
+ for(int i = 0; i < dim; i++)
+ {
+ constPermeability_[i][i] = interfaceSoilProps.ISP_Permeability;
+ }
+
+ }
+
+private:
+ MaterialLawParams materialLawParams_;
+ FieldMatrix constPermeability_;
+ Scalar porosity_;
+
+};
+
+} // end namespace
+#endif
diff --git a/appl/lecture/msm/mcwhorter/mcwhorterproblem.hh b/appl/lecture/msm/mcwhorter/mcwhorterproblem.hh
new file mode 100644
index 0000000000..76b99a0130
--- /dev/null
+++ b/appl/lecture/msm/mcwhorter/mcwhorterproblem.hh
@@ -0,0 +1,265 @@
+#ifndef DUMUX_MCWHORTERPROBLEM_HH
+#define DUMUX_MCWHORTERPROBLEM_HH
+
+#include <dune/grid/sgrid.hh>
+
+#include <dumux/material/fluidsystems/liquidphase.hh>
+//#include <dumux/material/components/simpleh2o.hh>
+//#include <dumux/material/components/oil.hh>
+
+#include "../buckleyleverett/pseudooil.hh"
+#include "../buckleyleverett/pseudoh2o.hh"
+
+#include <dumux/decoupled/2p/impes/impesproblem2p.hh>
+#include <dumux/decoupled/2p/diffusion/fv/fvvelocity2p.hh>
+#include <dumux/decoupled/2p/transport/fv/fvsaturation2p.hh>
+#include <dumux/decoupled/2p/transport/fv/capillarydiffusion.hh>
+#include <dumux/decoupled/2p/transport/fv/gravitypart.hh>
+
+
+#include "mcwhorter_spatialparams.hh"
+#include "mcwhorter_analytic.hh"
+
+namespace Dumux
+{
+template <class TypeTag>
+class McWhorterProblem;
+
+//////////
+// Specify the properties
+//////////
+namespace Properties
+{
+NEW_TYPE_TAG(McWhorterProblem, INHERITS_FROM(DecoupledTwoP, Transport));
+
+// Set the grid type
+SET_PROP(McWhorterProblem, Grid)
+{
+ // typedef Dune::YaspGrid<2> type;
+ typedef Dune::SGrid<1, 1> type;
+};
+
+// Set the problem property
+SET_PROP(McWhorterProblem, Problem)
+{
+public:
+ typedef Dumux::McWhorterProblem<TypeTag> type;
+};
+
+// Set the model properties
+SET_PROP(McWhorterProblem, TransportModel)
+{
+ typedef Dumux::FVSaturation2P<TTAG(McWhorterProblem)> type;
+};
+SET_TYPE_PROP(McWhorterProblem, DiffusivePart, Dumux::CapillaryDiffusion<TypeTag>);
+SET_TYPE_PROP(McWhorterProblem, ConvectivePart, Dumux::GravityPart<TypeTag>);
+
+SET_PROP(McWhorterProblem, PressureModel)
+{
+ typedef Dumux::FVVelocity2P<TTAG(McWhorterProblem)> type;
+};
+
+//SET_INT_PROP(McWhorterProblem, VelocityFormulation,
+// GET_PROP_TYPE(TypeTag, PTAG(TwoPIndices))::velocityW);
+
+SET_INT_PROP(McWhorterProblem, PressureFormulation,
+ GET_PROP_TYPE(TypeTag, PTAG(TwoPIndices))::pressureNW);
+
+// Set the wetting phase
+SET_PROP(McWhorterProblem, WettingPhase)
+{
+private:
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(Scalar)) Scalar;
+public:
+ typedef Dumux::LiquidPhase<Scalar, Dumux::PseudoOil<Scalar> > type;
+};
+
+// Set the non-wetting phase
+SET_PROP(McWhorterProblem, NonwettingPhase)
+{
+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(McWhorterProblem, SpatialParameters)
+{
+private:
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(Grid)) Grid;
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(Scalar)) Scalar;
+
+public:
+ typedef Dumux::McWhorterSpatialParams<TypeTag> type;
+};
+
+// Disable gravity
+SET_BOOL_PROP(McWhorterProblem, EnableGravity, false);
+
+SET_SCALAR_PROP(McWhorterProblem, CFLFactor, 0.1);
+}
+
+//! \ingroup transportProblems
+//! @brief McWhorter transport problem
+
+template<class TypeTag = TTAG(McWhorterProblem)>
+class McWhorterProblem: public IMPESProblem2P<TypeTag, McWhorterProblem<TypeTag> >
+{
+ typedef McWhorterProblem<TypeTag> ThisType;
+ typedef IMPESProblem2P<TypeTag, ThisType> ParentType;
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(GridView)) GridView;
+
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(TwoPIndices)) Indices;
+
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(FluidSystem)) FluidSystem;
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(FluidState)) FluidState;
+
+ enum
+ {
+ dim = GridView::dimension, dimWorld = GridView::dimensionworld
+ };
+ enum
+ {
+ wPhaseIdx = Indices::wPhaseIdx, nPhaseIdx = Indices::nPhaseIdx
+ };
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(Scalar)) Scalar;
+
+ typedef typename GridView::Traits::template Codim<0>::Entity Element;
+ typedef typename GridView::Intersection Intersection;
+ typedef Dune::FieldVector<Scalar, dimWorld> GlobalPosition;
+ typedef Dune::FieldVector<Scalar, dim> LocalPosition;
+
+public:
+
+ McWhorterProblem(const GridView &gridView,
+ const GlobalPosition lowerLeft = 0,
+ const GlobalPosition upperRight = 0,
+ const Scalar pleftbc = 2.0e5)//? initial oil pressure?
+ : ParentType(gridView),
+ lowerLeft_(lowerLeft),
+ upperRight_(upperRight),
+ eps_(1e-6),
+ pLeftBc_(pleftbc),
+ analyticSolution_(*this)
+ {}
+
+ /*!
+ * \brief The problem name.
+ *
+ * This is used as a prefix for files generated by the simulation.
+ */
+ const char *name() const
+ {
+ return "McWhorter";
+ }
+
+ bool shouldWriteRestartFile() const
+ {
+ return false;
+ }
+
+ void postTimeStep()
+ {
+ analyticSolution_.calculateAnalyticSolution();
+
+ ParentType::postTimeStep();
+ };
+
+ void addOutputVtkFields()
+ {
+ ParentType::addOutputVtkFields();
+ analyticSolution_.addOutputVtkFields(this->resultWriter());
+ }
+
+ bool shouldWriteOutput() const
+ {
+ if (this->timeManager().timeStepIndex() % 10 == 0)
+ {
+ return true;
+ }
+ return false;
+ }
+
+ /*!
+ * \brief Returns the temperature within the domain.
+ *
+ * This problem assumes a temperature of 10 degrees Celsius.
+ */
+
+ Scalar temperature(const GlobalPosition& globalPos, const Element& element) const
+ {return 273.15 + 10; // -> 10°C
+ }
+
+ Scalar referencePressure(const GlobalPosition& globalPos, const Element& element) const
+ {
+ return 1e5;
+ }
+ std::vector<Scalar> source (const GlobalPosition& globalPos, const Element& element)
+ {
+ return std::vector<Scalar>(2,0.0);//no source and sink terms
+ }
+
+ BoundaryConditions::Flags bctypePress(const GlobalPosition& globalPos, const Intersection& intersection) const
+ {
+ if (globalPos[0] < eps_)//west
+ return BoundaryConditions::dirichlet;
+
+ // all other boundaries
+ else
+ return BoundaryConditions::neumann;
+ }
+
+ BoundaryConditions::Flags bctypeSat(const GlobalPosition& globalPos, const Intersection& intersection) const
+ {
+ if (globalPos[0] < eps_)// west and east
+ return Dumux::BoundaryConditions::dirichlet;
+ else
+ return Dumux::BoundaryConditions::neumann;
+ }
+
+ Scalar dirichletPress(const GlobalPosition& globalPos, const Intersection& intersection) const
+ {
+ return pLeftBc_;
+ }
+
+ Scalar dirichletSat(const GlobalPosition& globalPos, const Intersection& intersection) const
+ {
+ if (globalPos[0] < eps_)//west
+ return 1.0;
+ else //east, north and south are Neumann
+ return 0.0;
+ }
+
+ std::vector<Scalar> neumannPress(const GlobalPosition& globalPos, const Intersection& intersection) const
+ {
+ return std::vector<Scalar>(2,0.0);
+ }
+
+ Scalar neumannSat(const GlobalPosition& globalPos, const Intersection& intersection, Scalar factor) const
+ {
+ return 0;
+ }
+
+ Scalar initSat (const GlobalPosition& globalPos, const Element& element) const
+ {
+ return 0.0;
+ }
+
+ /* McWhorterProblem(VC& variables, Fluid& wettingphase, Fluid& nonwettingphase,
+ Matrix2p<Grid, Scalar>& soil, TwoPhaseRelations<Grid, Scalar>& materialLaw = *(new TwoPhaseRelations<Grid,Scalar>&),
+ const GlobalPosition Right = 0)
+ : FractionalFlowProblem<GridView, Scalar, VC>(variables, wettingphase, nonwettingphase, soil, materialLaw),
+ right_(Right[0]), eps_(1e-6)
+ {}*/
+
+private:
+ GlobalPosition lowerLeft_;
+ GlobalPosition upperRight_;
+ Scalar eps_;
+ Scalar pLeftBc_;
+ Scalar right_;
+ McWhorterAnalytic<TypeTag> analyticSolution_;
+ };
+}
+#endif
--
GitLab