From 8b0957067fd1b2a6adf33c6fce8ccdbc9f086fc2 Mon Sep 17 00:00:00 2001 From: Michael Sinsbeck <michael@sinsbeck.com> Date: Thu, 18 Aug 2011 14:40:53 +0000 Subject: [PATCH] Corrected errors that were made in the last commit. Now it finally works git-svn-id: svn://svn.iws.uni-stuttgart.de/DUMUX/dumux/trunk@6521 2fb0f335-1f38-0410-981e-8018bf24f1b0 --- appl/lecture/msm/1p2cvs2p/dumux-in1p2c | 10 +- appl/lecture/msm/1p2cvs2p/dumux-in2p | 4 +- .../msm/1p2cvs2p/external_interface.hh | 389 ------------------ appl/lecture/msm/1p2cvs2p/lensproblem1p2c.hh | 53 +-- appl/lecture/msm/1p2cvs2p/lensproblem2p.hh | 38 +- .../msm/1p2cvs2p/lensspatialparameters1p2c.hh | 4 +- 6 files changed, 38 insertions(+), 460 deletions(-) delete mode 100644 appl/lecture/msm/1p2cvs2p/external_interface.hh diff --git a/appl/lecture/msm/1p2cvs2p/dumux-in1p2c b/appl/lecture/msm/1p2cvs2p/dumux-in1p2c index cc60aecacb..997a4cca5e 100644 --- a/appl/lecture/msm/1p2cvs2p/dumux-in1p2c +++ b/appl/lecture/msm/1p2cvs2p/dumux-in1p2c @@ -1,11 +1,11 @@ -MaxTimeStepSize = 500 -tEnd = 50000 -dt = 100 +MaxTimeStepSize = 4000 +tEnd = 64000 +dt = 1000 [Soil] -FinePermeability = 3.1e-11 +FinePermeability = 0.3e-11 CoarsePermeability = 5.89912e-11 -FinePorosity = 0.4 +FinePorosity = 0.5 CoarsePorosity = 0.5 [Boundary] diff --git a/appl/lecture/msm/1p2cvs2p/dumux-in2p b/appl/lecture/msm/1p2cvs2p/dumux-in2p index 5b06b58037..e0c4c8bd7c 100644 --- a/appl/lecture/msm/1p2cvs2p/dumux-in2p +++ b/appl/lecture/msm/1p2cvs2p/dumux-in2p @@ -1,5 +1,5 @@ -MaxTimeStepSize = 500 -tEnd = 50000 +MaxTimeStepSize = 5000 +tEnd = 500000 dt = 100 [Soil] diff --git a/appl/lecture/msm/1p2cvs2p/external_interface.hh b/appl/lecture/msm/1p2cvs2p/external_interface.hh deleted file mode 100644 index df16c95111..0000000000 --- a/appl/lecture/msm/1p2cvs2p/external_interface.hh +++ /dev/null @@ -1,389 +0,0 @@ -/***************************************************************************** - * This program is free software: you can redistribute it and/or modify * - * it under the terms of the GNU General Public License as published by * - * the Free Software Foundation, either version 2 of the License, or * - * (at your option) any later version. * - * * - * This program is distributed in the hope that it will be useful, * - * but WITHOUT ANY WARRANTY; without even the implied warranty of * - * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the * - * GNU General Public License for more details. * - * * - * You should have received a copy of the GNU General Public License * - * along with this program. If not, see <http://www.gnu.org/licenses/>. * - *****************************************************************************/ -#ifndef 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/lensproblem1p2c.hh b/appl/lecture/msm/1p2cvs2p/lensproblem1p2c.hh index 574828a93a..4c8e0aab93 100644 --- a/appl/lecture/msm/1p2cvs2p/lensproblem1p2c.hh +++ b/appl/lecture/msm/1p2cvs2p/lensproblem1p2c.hh @@ -82,8 +82,11 @@ SET_PROP(LensProblem, SpatialParameters) // Define whether mole(true) or mass(false) fractions are used SET_BOOL_PROP(LensProblem, UseMoles, true); -// Enable gravity +// Disable gravity SET_BOOL_PROP(LensProblem, EnableGravity, false); + +// Disable Jacobian recycling +SET_BOOL_PROP(LensProblem, EnableJacobianRecycling, false); } /*! @@ -214,18 +217,17 @@ public: * \param vertex The vertex on the boundary for which the * conditions needs to be specified */ - void boundaryTypes(BoundaryTypes &values, const Vertex &vertex) const + void boundaryTypesAtPos(BoundaryTypes &values, const GlobalPosition& globalPos) const { - const GlobalPosition globalPos = vertex.geometry().center(); - if (onUpperBoundary_(globalPos) || onLowerBoundary_(globalPos)) values.setAllDirichlet(); else values.setAllNeumann(); + if (onInlet_(globalPos)) values.setNeumann(transEqIdx); - if (onLowerBoundary_(globalPos)) - values.setNeumann(transEqIdx); +// if (onLowerBoundary_(globalPos)) +// values.setNeumann(transEqIdx); } /*! @@ -237,11 +239,15 @@ public: * * For this method, the \a values parameter stores primary variables. */ - void dirichlet(PrimaryVariables &values, const Vertex &vertex) const + void dirichletAtPos(PrimaryVariables &values, const GlobalPosition& globalPos) const { - const GlobalPosition globalPos = vertex.geometry().center(); - - if (onUpperBoundary_(globalPos)) +// if (onInlet_(globalPos) && this->timeManager().time() <= infiltrationEndTime_) +// { +// values[pressureIdx] = upperPressure_; +// values[x1Idx] =0.8; +// } +// else + if (onUpperBoundary_(globalPos)) { values[pressureIdx] = upperPressure_; values[x1Idx] = 0.0; @@ -262,21 +268,11 @@ public: * 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 + void neumannAtPos(PrimaryVariables &values, const GlobalPosition& globalPos) const { - const GlobalPosition &globalPos - = element.geometry().corner(scvIdx); - values = 0.0; - const Scalar& time = this->timeManager().time(); - - if (time >= infiltrationStartTime_ && time <= infiltrationEndTime_) + if (time <= infiltrationEndTime_ && infiltrationStartTime_ <= time) { if (onInlet_(globalPos)) values[transEqIdx] = -infiltrationRate_; @@ -297,10 +293,7 @@ public: * 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 + void sourceAtPos(PrimaryVariables &values, const GlobalPosition& globalPos) const { values = Scalar(0.0); } @@ -311,14 +304,8 @@ public: * For this method, the \a values parameter stores primary * variables. */ - void initial(PrimaryVariables &values, - const Element &element, - const FVElementGeometry &fvElemGeom, - int scvIdx) const + void initialAtPos(PrimaryVariables &values, const GlobalPosition& globalPos) 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]; diff --git a/appl/lecture/msm/1p2cvs2p/lensproblem2p.hh b/appl/lecture/msm/1p2cvs2p/lensproblem2p.hh index e62a29d3a7..c54d50739c 100644 --- a/appl/lecture/msm/1p2cvs2p/lensproblem2p.hh +++ b/appl/lecture/msm/1p2cvs2p/lensproblem2p.hh @@ -188,6 +188,7 @@ public: infiltrationRate_ = Params::tree().template get<double>("Boundary.InfiltrationRate"); infiltrationStartTime_= 1.0e-9;//The infiltrations starts always after the first time step! infiltrationEndTime_= Params::tree().template get<double>("Boundary.InfiltrationEndTime"); + backgroundSwr_ = Params::tree().template get<double>("Soil.CoarseResidualSaturationWetting"); } /*! @@ -232,18 +233,15 @@ public: * \param vertex The vertex on the boundary for which the * conditions needs to be specified */ - void boundaryTypes(BoundaryTypes &values, const Vertex &vertex) const + void boundaryTypesAtPos(BoundaryTypes &values, const GlobalPosition& globalPos) const { - const GlobalPosition globalPos = vertex.geometry().center(); - - if (onUpperBoundary_(globalPos) || onLowerBoundary_(globalPos)) values.setAllDirichlet(); else values.setAllNeumann(); if (onInlet_(globalPos)) - values.setNeumann(contiNEqIdx); + values.setAllNeumann(); } /*! @@ -255,10 +253,8 @@ public: * * For this method, the \a values parameter stores primary variables. */ - void dirichlet(PrimaryVariables &values, const Vertex &vertex) const + void dirichletAtPos(PrimaryVariables &values, const GlobalPosition& globalPos) const { - const GlobalPosition globalPos = vertex.geometry().center(); - if (onUpperBoundary_(globalPos)) { values[pwIdx] = upperPressure_; @@ -301,21 +297,13 @@ public: * 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 + void neumannAtPos(PrimaryVariables &values, const GlobalPosition& globalPos) const { - const GlobalPosition &globalPos - = element.geometry().corner(scvIdx); - values = 0.0; const Scalar& time = this->timeManager().time(); - if (time >= infiltrationStartTime_ && time <= infiltrationEndTime_) + if (time <= infiltrationEndTime_) { if (onInlet_(globalPos)) values[contiNEqIdx] = -infiltrationRate_; // kg / (m * s) @@ -337,10 +325,7 @@ public: * 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 + void sourceAtPos(PrimaryVariables &values, const GlobalPosition& globalPos) const { values = Scalar(0.0); } @@ -351,14 +336,8 @@ public: * For this method, the \a values parameter stores primary * variables. */ - void initial(PrimaryVariables &values, - const Element &element, - const FVElementGeometry &fvElemGeom, - int scvIdx) const + void initialAtPos(PrimaryVariables &values, const GlobalPosition& globalPos) 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]; @@ -406,6 +385,7 @@ private: Scalar infiltrationRate_; Scalar infiltrationStartTime_; Scalar infiltrationEndTime_; + Scalar backgroundSwr_; }; } //end namespace diff --git a/appl/lecture/msm/1p2cvs2p/lensspatialparameters1p2c.hh b/appl/lecture/msm/1p2cvs2p/lensspatialparameters1p2c.hh index 0738722993..ded344020c 100644 --- a/appl/lecture/msm/1p2cvs2p/lensspatialparameters1p2c.hh +++ b/appl/lecture/msm/1p2cvs2p/lensspatialparameters1p2c.hh @@ -72,8 +72,8 @@ public: lensK_(0), outerK_(0) { - lensPorosity_ = Params::tree().template get<double>("Soil.FinePermeability"); - outerPorosity_ = Params::tree().template get<double>("Soil.CoarsePermeability"); + lensPorosity_ = Params::tree().template get<double>("Soil.FinePorosity"); + outerPorosity_ = Params::tree().template get<double>("Soil.CoarsePorosity"); longitudinalDispersivity_ = 1.0e-5; transverseDispersivity_ = 1.0e-6; -- GitLab