Commit cf0fb664 authored by Bernd Flemisch's avatar Bernd Flemisch
Browse files

Reverting previous commits and going back to revision 6808.

There is no need to rush things. We agreed on first having the meeting
and then porting the change to the non-MpNc-models. Moreover, the change
has not been carried out according to what has been discussed at the
meeting. 
E.g. there was no completeState functionality of the models, and the simple
models also introduced parameter caches, although it was thought to not
necessarily having to do so.


git-svn-id: svn://svn.iws.uni-stuttgart.de/DUMUX/dumux/trunk@6818 2fb0f335-1f38-0410-981e-8018bf24f1b0
parent abb68f4a
......@@ -45,16 +45,16 @@ namespace Dumux
struct Source
{
Dune::FieldVector<double,2> globalPos;
double q;
int index;
Dune::FieldVector<double,2> globalPos;
double q;
int index;
};
struct BoundarySegment
{
double from, to;
bool neumann;
double value;
double from, to;
bool neumann;
double value;
};
template<class TypeTag>
......@@ -140,74 +140,74 @@ public:
{
// this->spatialParameters().setDelta(delta_);
// Write input parameters into private variables
//Dune::FieldVector<int,2> resolution = Params::tree().template get<Dune::FieldVector<int,2> >("Geometry.numberOfCells");
// Write input parameters into private variables
Dune::FieldVector<int,2> resolution = Params::tree().template get<Dune::FieldVector<int,2> >("Geometry.numberOfCells");
domainSize_ = Params::tree().template get<GlobalPosition>("Geometry.domainSize");
geometryDepth_ = Params::tree().template get<double>("Geometry.depth");
resolution_ = Params::tree().template get<Dune::FieldVector<int,2>>("Geometry.numberOfCells");
// Read sources
std::vector<double> sources = Params::tree().template get<std::vector<double> >("Source.sources");
int NumberOfSources = std::trunc(sources.size()/3);
for (int sourceCount=0; sourceCount<NumberOfSources ; sourceCount++)
{
Source tempSource;
tempSource.globalPos[0]=sources[sourceCount*3];
tempSource.globalPos[1]=sources[sourceCount*3+1];
tempSource.q=sources[sourceCount*3+2];
tempSource.index = std::floor(tempSource.globalPos[0]
* resolution_[0]/domainSize_[0])
+ std::floor(tempSource.globalPos[1]*resolution_[1]/domainSize_[1])
* resolution_[0];
sources_.push_back(tempSource);
}
// Read Boundary Conditions
std::vector<double> BC = Params::tree().template get<std::vector<double> >("BoundaryConditions.left");
int NumberOfSegments = std::trunc(BC.size()/4);
for (int segmentCount=0; segmentCount<NumberOfSegments ; segmentCount++)
{
BoundarySegment tempSegment;
tempSegment.from = BC[segmentCount*4];
tempSegment.to = BC[segmentCount*4+1];
tempSegment.neumann = (BC[segmentCount*4+2]!=0);
tempSegment.value = BC[segmentCount*4+3];
boundaryConditions_[2].push_back(tempSegment);
}
BC = Params::tree().template get<std::vector<double> >("BoundaryConditions.right");
NumberOfSegments = std::trunc(BC.size()/4);
for (int segmentCount=0; segmentCount<NumberOfSegments ; segmentCount++)
{
BoundarySegment tempSegment;
tempSegment.from = BC[segmentCount*4];
tempSegment.to = BC[segmentCount*4+1];
tempSegment.neumann = (BC[segmentCount*4+2]!=0);
tempSegment.value = BC[segmentCount*4+3];
boundaryConditions_[3].push_back(tempSegment);
}
BC = Params::tree().template get<std::vector<double> >("BoundaryConditions.bottom");
NumberOfSegments = std::trunc(BC.size()/4);
for (int segmentCount=0; segmentCount<NumberOfSegments ; segmentCount++)
{
BoundarySegment tempSegment;
tempSegment.from = BC[segmentCount*4];
tempSegment.to = BC[segmentCount*4+1];
tempSegment.neumann = (BC[segmentCount*4+2]!=0);
tempSegment.value = BC[segmentCount*4+3];
boundaryConditions_[1].push_back(tempSegment);
}
BC = Params::tree().template get<std::vector<double> >("BoundaryConditions.top");
NumberOfSegments = std::trunc(BC.size()/4);
for (int segmentCount=0; segmentCount<NumberOfSegments ; segmentCount++)
{
BoundarySegment tempSegment;
tempSegment.from = BC[segmentCount*4];
tempSegment.to = BC[segmentCount*4+1];
tempSegment.neumann = (BC[segmentCount*4+2]!=0);
tempSegment.value = BC[segmentCount*4+3];
boundaryConditions_[0].push_back(tempSegment);
}
std::vector<double> sources = Params::tree().template get<std::vector<double> >("Source.sources");
int NumberOfSources = std::trunc(sources.size()/3);
for (int sourceCount=0; sourceCount<NumberOfSources ; sourceCount++)
{
Source tempSource;
tempSource.globalPos[0]=sources[sourceCount*3];
tempSource.globalPos[1]=sources[sourceCount*3+1];
tempSource.q=sources[sourceCount*3+2];
tempSource.index = std::floor(tempSource.globalPos[0]
* resolution_[0]/domainSize_[0])
+ std::floor(tempSource.globalPos[1]*resolution_[1]/domainSize_[1])
* resolution_[0];
sources_.push_back(tempSource);
}
// Read Boundary Conditions
std::vector<double> BC = Params::tree().template get<std::vector<double> >("BoundaryConditions.left");
int NumberOfSegments = std::trunc(BC.size()/4);
for (int segmentCount=0; segmentCount<NumberOfSegments ; segmentCount++)
{
BoundarySegment tempSegment;
tempSegment.from = BC[segmentCount*4];
tempSegment.to = BC[segmentCount*4+1];
tempSegment.neumann = (BC[segmentCount*4+2]!=0);
tempSegment.value = BC[segmentCount*4+3];
boundaryConditions_[2].push_back(tempSegment);
}
BC = Params::tree().template get<std::vector<double> >("BoundaryConditions.right");
NumberOfSegments = std::trunc(BC.size()/4);
for (int segmentCount=0; segmentCount<NumberOfSegments ; segmentCount++)
{
BoundarySegment tempSegment;
tempSegment.from = BC[segmentCount*4];
tempSegment.to = BC[segmentCount*4+1];
tempSegment.neumann = (BC[segmentCount*4+2]!=0);
tempSegment.value = BC[segmentCount*4+3];
boundaryConditions_[3].push_back(tempSegment);
}
BC = Params::tree().template get<std::vector<double> >("BoundaryConditions.bottom");
NumberOfSegments = std::trunc(BC.size()/4);
for (int segmentCount=0; segmentCount<NumberOfSegments ; segmentCount++)
{
BoundarySegment tempSegment;
tempSegment.from = BC[segmentCount*4];
tempSegment.to = BC[segmentCount*4+1];
tempSegment.neumann = (BC[segmentCount*4+2]!=0);
tempSegment.value = BC[segmentCount*4+3];
boundaryConditions_[1].push_back(tempSegment);
}
BC = Params::tree().template get<std::vector<double> >("BoundaryConditions.top");
NumberOfSegments = std::trunc(BC.size()/4);
for (int segmentCount=0; segmentCount<NumberOfSegments ; segmentCount++)
{
BoundarySegment tempSegment;
tempSegment.from = BC[segmentCount*4];
tempSegment.to = BC[segmentCount*4+1];
tempSegment.neumann = (BC[segmentCount*4+2]!=0);
tempSegment.value = BC[segmentCount*4+3];
boundaryConditions_[0].push_back(tempSegment);
}
}
/*!
......@@ -251,10 +251,10 @@ public:
{
values = 0;
Scalar density=Fluid::density(0,0);
for (int sourceCount = 0; sourceCount != sources_.size(); sourceCount++)
for (int sourceCount = 0; sourceCount != sources_.size(); sourceCount++)
{
if (this->variables().index(element) == sources_[sourceCount].index)
values+=sources_[sourceCount].q*density/element.geometry().volume()/geometryDepth_;
if (this->variables().index(element) == sources_[sourceCount].index)
values+=sources_[sourceCount].q*density/element.geometry().volume()/geometryDepth_;
}
}
......@@ -293,12 +293,12 @@ public:
boundaryIndex=0;
}
for (int segmentCount=0; segmentCount<boundaryConditions_[boundaryIndex].size();segmentCount++)
for (int segmentCount=0; segmentCount<boundaryConditions_[boundaryIndex].size();segmentCount++)
{
if ((boundaryConditions_[boundaryIndex][segmentCount].from < coordinate) &&
(coordinate < boundaryConditions_[boundaryIndex][segmentCount].to))
{
if (boundaryConditions_[boundaryIndex][segmentCount].neumann)
if ((boundaryConditions_[boundaryIndex][segmentCount].from < coordinate) &&
(coordinate < boundaryConditions_[boundaryIndex][segmentCount].to))
{
if (boundaryConditions_[boundaryIndex][segmentCount].neumann)
bcType.setAllNeumann();
else
bcType.setAllDirichlet();
......@@ -341,11 +341,11 @@ public:
boundaryIndex=0;
}
for (int segmentCount=0; segmentCount<boundaryConditions_[boundaryIndex].size();segmentCount++)
{
if ((boundaryConditions_[boundaryIndex][segmentCount].from < coordinate) &&
(coordinate < boundaryConditions_[boundaryIndex][segmentCount].to))
{
for (int segmentCount=0; segmentCount<boundaryConditions_[boundaryIndex].size();segmentCount++)
{
if ((boundaryConditions_[boundaryIndex][segmentCount].from < coordinate) &&
(coordinate < boundaryConditions_[boundaryIndex][segmentCount].to))
{
values = boundaryConditions_[boundaryIndex][segmentCount].value*Fluid::density(0,0)*9.81;
return;
}
......@@ -382,12 +382,12 @@ public:
coordinate=globalPos[0];
boundaryIndex=0;
}
for (int segmentCount=0; segmentCount<boundaryConditions_[boundaryIndex].size();segmentCount++)
{
if ((boundaryConditions_[boundaryIndex][segmentCount].from < coordinate) &&
(coordinate < boundaryConditions_[boundaryIndex][segmentCount].to))
{
for (int segmentCount=0; segmentCount<boundaryConditions_[boundaryIndex].size();segmentCount++)
{
if ((boundaryConditions_[boundaryIndex][segmentCount].from < coordinate) &&
(coordinate < boundaryConditions_[boundaryIndex][segmentCount].to))
{
values = boundaryConditions_[boundaryIndex][segmentCount].value*Fluid::density(0,0)*(-1);
return;
}
......@@ -398,73 +398,73 @@ public:
void writeOutput()
{
Dune::FieldVector<int,2> resolution = Params::tree().template get<Dune::FieldVector<int,2> >("Geometry.numberOfCells");
Scalar zmax, zmin;
zmax=this->variables().pressure()[0]/(Fluid::density(0,0)*9.81);
zmin=this->variables().pressure()[0]/(Fluid::density(0,0)*9.81);
for (int i=0; i< resolution[0]*resolution[1]; i++)
{
Scalar currentHead= this->variables().pressure()[i]/(Fluid::density(0,0)*9.81);
zmax = std::max(currentHead,zmax);
zmin = std::min(currentHead,zmin);
}
std::ofstream dataFile;
dataFile.open("dumux-out.vgfc");
dataFile << "Gridplot" << std::endl;
dataFile << "## This is a DuMuX output for the ViPLab Graphics driver. \n";
dataFile << "## This output file was generated at " << __TIME__ <<", "<< __DATE__<< "\n";
dataFile << "# x-range 0 "<< domainSize_[0] << "\n" ;
dataFile << "# y-range 0 "<< domainSize_[1] << "\n" ;
dataFile << "# x-count " << resolution_[0] << "\n" ;
dataFile << "# y-count " << resolution_[1] << "\n" ;
if ((zmax-zmin)/zmax>0.01)
dataFile << "# scale 1 1 "<< sqrt(domainSize_[0]*domainSize_[1])/(zmax-zmin) << "\n";
else
dataFile << "# scale 1 1 1\n";
dataFile << "# min-color 255 0 0\n";
dataFile << "# max-color 0 0 255\n";
dataFile << "# time 0 \n" ;
dataFile << "# label piezometric head \n";
for (int i=0; i< resolution_[1]; i++)
{
for (int j=0; j<resolution_[0]; j++)
{
int currentIdx = i*resolution_[0]+j;
dataFile << this->variables().pressure()[currentIdx]/(Fluid::density(0,0)*9.81);
if(j != resolution_[0]-1) // all but last entry
dataFile << " ";
else // write the last entry
dataFile << "\n";
}
}
dataFile.close();
//Textoutput:
std::cout << " x y h v_x v_y"<<std::endl;
std::cout << "------------------------------------------------------------"<<std::endl;
ElementIterator eItEnd = this->gridView().template end<0> ();
for (ElementIterator eIt = this->gridView().template begin<0> (); eIt != eItEnd; ++eIt)
{
int cellIndex = this->variables().index(*eIt);
double v_x,v_y,piezo,x,y;
v_x= (this->variables().velocity()[cellIndex][0][0]+this->variables().velocity()[cellIndex][1][0])/2;
v_y= (this->variables().velocity()[cellIndex][2][1]+this->variables().velocity()[cellIndex][3][1])/2;
if (std::abs(v_x)<1e-17)
v_x=0;
if (std::abs(v_y)<1e-17)
v_y=0;
piezo=this->variables().pressure()[cellIndex]/(Fluid::density(0,0)*9.81);
x = eIt->geometry().center()[0];
y = eIt->geometry().center()[1];
printf("%10.4g %10.4g %10.4g %13.4g %13.4g\n",x,y,piezo,v_x,v_y);
}
Scalar zmax, zmin;
zmax=this->variables().pressure()[0]/(Fluid::density(0,0)*9.81);
zmin=this->variables().pressure()[0]/(Fluid::density(0,0)*9.81);
for (int i=0; i< resolution[0]*resolution[1]; i++)
{
Scalar currentHead= this->variables().pressure()[i]/(Fluid::density(0,0)*9.81);
zmax = std::max(currentHead,zmax);
zmin = std::min(currentHead,zmin);
}
std::ofstream dataFile;
dataFile.open("dumux-out.vgfc");
dataFile << "Gridplot" << std::endl;
dataFile << "## This is a DuMuX output for the ViPLab Graphics driver. \n";
dataFile << "## This output file was generated at " << __TIME__ <<", "<< __DATE__<< "\n";
dataFile << "# x-range 0 "<< domainSize_[0] << "\n" ;
dataFile << "# y-range 0 "<< domainSize_[1] << "\n" ;
dataFile << "# x-count " << resolution_[0] << "\n" ;
dataFile << "# y-count " << resolution_[1] << "\n" ;
if ((zmax-zmin)/zmax>0.01)
dataFile << "# scale 1 1 "<< sqrt(domainSize_[0]*domainSize_[1])/(zmax-zmin) << "\n";
else
dataFile << "# scale 1 1 1\n";
dataFile << "# min-color 255 0 0\n";
dataFile << "# max-color 0 0 255\n";
dataFile << "# time 0 \n" ;
dataFile << "# label piezometric head \n";
for (int i=0; i< resolution_[1]; i++)
{
for (int j=0; j<resolution_[0]; j++)
{
int currentIdx = i*resolution_[0]+j;
dataFile << this->variables().pressure()[currentIdx]/(Fluid::density(0,0)*9.81);
if(j != resolution_[0]-1) // all but last entry
dataFile << " ";
else // write the last entry
dataFile << "\n";
}
}
dataFile.close();
//Textoutput:
std::cout << " x y h v_x v_y"<<std::endl;
std::cout << "------------------------------------------------------------"<<std::endl;
ElementIterator eItEnd = this->gridView().template end<0> ();
for (ElementIterator eIt = this->gridView().template begin<0> (); eIt != eItEnd; ++eIt)
{
int cellIndex = this->variables().index(*eIt);
double v_x,v_y,piezo,x,y;
v_x= (this->variables().velocity()[cellIndex][0][0]+this->variables().velocity()[cellIndex][1][0])/2;
v_y= (this->variables().velocity()[cellIndex][2][1]+this->variables().velocity()[cellIndex][3][1])/2;
if (std::abs(v_x)<1e-17)
v_x=0;
if (std::abs(v_y)<1e-17)
v_y=0;
piezo=this->variables().pressure()[cellIndex]/(Fluid::density(0,0)*9.81);
x = eIt->geometry().center()[0];
y = eIt->geometry().center()[1];
printf("%10.4g %10.4g %10.4g %13.4g %13.4g\n",x,y,piezo,v_x,v_y);
}
}
private:
......
......@@ -32,10 +32,10 @@ namespace Dumux
struct Lens
{
Dune::FieldMatrix<double,2,2> permeability;
//double permeability;
Dune::FieldVector<double,2> lowerLeft;
Dune::FieldVector<double,2> upperRight;
Dune::FieldMatrix<double,2,2> permeability;
//double permeability;
Dune::FieldVector<double,2> lowerLeft;
Dune::FieldVector<double,2> upperRight;
};
......@@ -91,33 +91,33 @@ public:
GroundwaterSpatialParams(const GridView& gridView)
: FVSpatialParametersOneP<TypeTag>(gridView), permeability_(0)
{
delta_=1e-3;
porosity_ = 0.2;
Scalar permFactor = 0.001/(1000*9.81);
permeability_[0][0] = Params::tree().template get<double>("SpatialParameters.permeability")*permFactor;
permeability_[1][1] = permeability_[0][0];
permeability_[0][1] = 0;
permeability_[1][0] = 0;
//Lenses:
std::vector<double> lenses = Params::tree().template get<std::vector<double>>("SpatialParameters.lenses");
int NumberOfLenses = std::trunc(lenses.size()/5);
for (int lensCount=0; lensCount<NumberOfLenses ; lensCount++)
{
Lens tempLens;
tempLens.lowerLeft[0]=lenses[lensCount*5];
tempLens.upperRight[0]=lenses[lensCount*5+1];
tempLens.lowerLeft[1]=lenses[lensCount*5+2];
tempLens.upperRight[1]=lenses[lensCount*5+3];
tempLens.permeability[0][0]=lenses[lensCount*5+4]*permFactor;
tempLens.permeability[1][1] = tempLens.permeability[0][0];
tempLens.permeability[0][1] = 0;
tempLens.permeability[1][0] = 0;
lenses_.push_back(tempLens);
}
delta_=1e-3;
porosity_ = 0.2;
Scalar permFactor = 0.001/(1000*9.81);
permeability_[0][0] = Params::tree().template get<double>("SpatialParameters.permeability")*permFactor;
permeability_[1][1] = permeability_[0][0];
permeability_[0][1] = 0;
permeability_[1][0] = 0;
//Lenses:
std::vector<double> lenses = Params::tree().template get<std::vector<double>>("SpatialParameters.lenses");
int NumberOfLenses = std::trunc(lenses.size()/5);
for (int lensCount=0; lensCount<NumberOfLenses ; lensCount++)
{
Lens tempLens;
tempLens.lowerLeft[0]=lenses[lensCount*5];
tempLens.upperRight[0]=lenses[lensCount*5+1];
tempLens.lowerLeft[1]=lenses[lensCount*5+2];
tempLens.upperRight[1]=lenses[lensCount*5+3];
tempLens.permeability[0][0]=lenses[lensCount*5+4]*permFactor;
tempLens.permeability[1][1] = tempLens.permeability[0][0];
tempLens.permeability[0][1] = 0;
tempLens.permeability[1][0] = 0;
lenses_.push_back(tempLens);
}
}
private:
......
......@@ -181,17 +181,17 @@ public:
// Write header for output file
std::ofstream dataFile;
dataFile.open("dumux-out.vgfc");
// dataFile << "Gridplot" << std::endl;
dataFile << "## This is a DuMuX output for the ViPLab graphics driver. \n";
dataFile << "## This output file was generated at " << __TIME__ <<", "<< __DATE__<< "\n";
dataFile << "# x-range " << this->bboxMin()[0] << " " << this->bboxMax()[0] << "\n" ;
dataFile << "# y-range " << this->bboxMin()[1] << " " << this->bboxMax()[1] << "\n" ;
dataFile << "# x-count " << resX_+1 << "\n" ;
dataFile << "# y-count " << resY_+1 << "\n" ;
// dataFile << "# min-color 0 0 0\n";
// dataFile << "# max-color 255 255 255\n";
dataFile.close();
dataFile.open("dumux-out.vgfc");
// dataFile << "Gridplot" << std::endl;
dataFile << "## This is a DuMuX output for the ViPLab graphics driver. \n";
dataFile << "## This output file was generated at " << __TIME__ <<", "<< __DATE__<< "\n";
dataFile << "# x-range " << this->bboxMin()[0] << " " << this->bboxMax()[0] << "\n" ;
dataFile << "# y-range " << this->bboxMin()[1] << " " << this->bboxMax()[1] << "\n" ;
dataFile << "# x-count " << resX_+1 << "\n" ;
dataFile << "# y-count " << resY_+1 << "\n" ;
// dataFile << "# min-color 0 0 0\n";
// dataFile << "# max-color 255 255 255\n";
dataFile.close();
}
/*!
......@@ -260,13 +260,13 @@ public:
*/
void dirichletAtPos(PrimaryVariables &values, const GlobalPosition& globalPos) const
{
// if (onInlet_(globalPos) && this->timeManager().time() <= infiltrationEndTime_)
// {
// values[pressureIdx] = upperPressure_;
// values[x1Idx] =0.8;
// }
// else
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;
......@@ -340,36 +340,36 @@ public:
//
void writeOutput()
{
//todo: Das hier ist noch nicht fertig.
//todo: Das hier ist noch nicht fertig.
Scalar time = this->timeManager().time();
if (time<0)
return;
SolutionVector &sol = this->model().curSol();
Scalar time = this->timeManager().time();
if (time<0)
return;
SolutionVector &sol = this->model().curSol();
std::ofstream dataFile;
dataFile.open("dumux-out.vgfc", std::fstream::app);
dataFile << "# time "<< time <<"\n" ;
dataFile << "# label Concentration \n";
dataFile << "# min-color 0 0 0\n";
dataFile << "# max-color 255 255 255\n";
for (int j=0; j < resY_+1; j++)
{
for (int i=0; i < resX_+1; i++)
{
int currentIdx = i*(resY_+1)+j;
dataFile << sol[currentIdx][x1Idx];
if(i != resX_) // all but last entry
dataFile << " ";
else // write the last entry
{
dataFile << "\n";
}
}
}
dataFile.close();
dataFile.open("dumux-out.vgfc", std::fstream::app);
dataFile << "# time "<< time <<"\n" ;
dataFile << "# label Concentration \n";
dataFile << "# min-color 0 0 0\n";
dataFile << "# max-color 255 255 255\n";
for (int j=0; j < resY_+1; j++)
{
for (int i=0; i < resX_+1; i++)
{
int currentIdx = i*(resY_+1)+j;
dataFile << sol[currentIdx][x1Idx];
if(i != resX_) // all but last entry
dataFile << " ";
else // write the last entry
{