diff --git a/appl/lecture/mhs/groundwater/groundwater_problem.hh b/appl/lecture/mhs/groundwater/groundwater_problem.hh
index 9b4af5f783c79f404116c64fe2db1e8912251cf2..eb8f756b8863fbd4885976126253498ff63df225 100644
--- a/appl/lecture/mhs/groundwater/groundwater_problem.hh
+++ b/appl/lecture/mhs/groundwater/groundwater_problem.hh
@@ -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
+ // 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;
}
@@ -383,11 +383,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)*(-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:
diff --git a/appl/lecture/mhs/groundwater/groundwater_spatialparams.hh b/appl/lecture/mhs/groundwater/groundwater_spatialparams.hh
index fb545d42eeacd5cac5df89b209c1b51d1e250f0d..c8fefa30a0567e84fb77d9c575f47872ec3120de 100644
--- a/appl/lecture/mhs/groundwater/groundwater_spatialparams.hh
+++ b/appl/lecture/mhs/groundwater/groundwater_spatialparams.hh
@@ -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:
diff --git a/appl/lecture/msm/1p2cvs2p/lensproblem1p2c.hh b/appl/lecture/msm/1p2cvs2p/lensproblem1p2c.hh
index 958c3d44738701ec68b25a0babced8e40e9ddac5..337a3b2daf0c3cabcd8aebb4a87c80acd71176ee 100644
--- a/appl/lecture/msm/1p2cvs2p/lensproblem1p2c.hh
+++ b/appl/lecture/msm/1p2cvs2p/lensproblem1p2c.hh
@@ -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
+ {
+ dataFile << "\n";
+ }
+ }
+ }
+ dataFile.close();
}
diff --git a/appl/lecture/msm/1p2cvs2p/lensproblem2p.hh b/appl/lecture/msm/1p2cvs2p/lensproblem2p.hh
index 16a381d9e4e458729f7f40f8a54cd59d38e1ebc1..a5df406e877b848b53452bd717a8aa30aa186931 100644
--- a/appl/lecture/msm/1p2cvs2p/lensproblem2p.hh
+++ b/appl/lecture/msm/1p2cvs2p/lensproblem2p.hh
@@ -241,7 +241,7 @@ public:
values.setAllNeumann();
if (onInlet_(globalPos))
- values.setAllNeumann();
+ values.setAllNeumann();
}
/*!
diff --git a/appl/lecture/msm/buckleyleverett/buckleyleverett_analytic.hh b/appl/lecture/msm/buckleyleverett/buckleyleverett_analytic.hh
index 0c4131168b852f33c831df81eaf4bbc9e805fde7..f5cea99cf8a1c78ee2256edb00746fcf4b234dc7 100644
--- a/appl/lecture/msm/buckleyleverett/buckleyleverett_analytic.hh
+++ b/appl/lecture/msm/buckleyleverett/buckleyleverett_analytic.hh
@@ -333,7 +333,7 @@ public:
BlockVector AnalyticSolution() const
{
- return analyticSolution_;
+ return analyticSolution_;
}
//Write saturation and pressure into file
@@ -357,8 +357,8 @@ public:
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();
+ initializeAnalytic();
+ prepareAnalytic();
}
diff --git a/appl/lecture/msm/buckleyleverett/buckleyleverett_spatialparams.hh b/appl/lecture/msm/buckleyleverett/buckleyleverett_spatialparams.hh
index 2196716478edef836305fc16dfddd9804fa2c7c1..6e7dc14147832d032d73f6b8487d10eba4f45f1e 100644
--- a/appl/lecture/msm/buckleyleverett/buckleyleverett_spatialparams.hh
+++ b/appl/lecture/msm/buckleyleverett/buckleyleverett_spatialparams.hh
@@ -96,16 +96,16 @@ public:
BuckleyLeverettSpatialParams(const GridView& gridView)
:ParentType(gridView)
{
- Scalar permFactor = 0.001/(1000*9.81);
+ Scalar permFactor = 0.001/(1000*9.81);
- constPermeability_ = Params::tree().template get<double>("SpatialParameters.permeability")*permFactor;
- //Lenses:
- materialLawParams_.setSwr(Params::tree().template get<double>("SpatialParameters.swr"));
- materialLawParams_.setSnr(Params::tree().template get<double>("SpatialParameters.snr"));
- materialLawParams_.setPe(Params::tree().template get<double>("SpatialParameters.pe"));
- materialLawParams_.setLambda(Params::tree().template get<double>("SpatialParameters.lambda"));
+ constPermeability_ = Params::tree().template get<double>("SpatialParameters.permeability")*permFactor;
+ //Lenses:
+ materialLawParams_.setSwr(Params::tree().template get<double>("SpatialParameters.swr"));
+ materialLawParams_.setSnr(Params::tree().template get<double>("SpatialParameters.snr"));
+ materialLawParams_.setPe(Params::tree().template get<double>("SpatialParameters.pe"));
+ materialLawParams_.setLambda(Params::tree().template get<double>("SpatialParameters.lambda"));
- porosity_ = Params::tree().template get<double>("SpatialParameters.porosity");
+ porosity_ = Params::tree().template get<double>("SpatialParameters.porosity");
}
private:
diff --git a/appl/lecture/msm/buckleyleverett/buckleyleverettproblem.hh b/appl/lecture/msm/buckleyleverett/buckleyleverettproblem.hh
index cd5a69edaeb51d15451880b6dacb1f3edbe87ac5..b3bd2504d523d098bf058d7076a0c9659c051b26 100644
--- a/appl/lecture/msm/buckleyleverett/buckleyleverettproblem.hh
+++ b/appl/lecture/msm/buckleyleverett/buckleyleverettproblem.hh
@@ -154,9 +154,9 @@ public:
densityNonWetting_ = Params::tree().template get<double>("Fluid.densityNW");
//Write header for ViPLab-Outputfile
- std::ofstream dataFile;
- dataFile.open("vipplot.vgf");
- dataFile.close();
+ std::ofstream dataFile;
+ dataFile.open("vipplot.vgf");
+ dataFile.close();
}
@@ -275,85 +275,85 @@ public:
//Override outputfunction for ViPLab-Output
void writeOutput()
{
- double time = this->timeManager().time();
- if (time < 0)
- return;
+ double time = this->timeManager().time();
+ if (time < 0)
+ return;
double discretizationLength = Params::tree().template get<double>("Geometry.discretizationLength");
- int cellNumberX = static_cast<int>(300/discretizationLength);
- discretizationLength = 300/cellNumberX; //Might be slightly different from the input parameter
-
- std::ofstream dataFile;
- dataFile.open("vipplot.vgf", std::fstream::app);
-
- std::cout<<"Writing output for time step.\n";
-
- if (time > 0)
- dataFile << "# newframe\n";
-
- dataFile << "# title Buckley-Leverett-Problem\n";
- dataFile << "# x-label x\n";
- dataFile << "# y-label Saturation\n";
- dataFile << "# x-range -10 310\n";
- dataFile << "# y-range 0 1\n";
-
- dataFile << "# color 0 0 255\n";
-
- Scalar curSat = this->variables().saturation()[0];
-
- // Draw first piece separately, as no vertical line is needed
- dataFile << 0 << " " << curSat << " "
- << discretizationLength << " " << curSat <<std::endl;
- for (int i=1; i < cellNumberX; i++)
- {
- // Vertical Line
- dataFile << i*discretizationLength << " "
- << curSat << " ";
- curSat = this->variables().saturation()[i];
- dataFile << i*discretizationLength << " " << curSat << std::endl;
- //Horizontal Line
- dataFile << i*discretizationLength << " " << curSat << " "
- << (i+1)*discretizationLength << " " << curSat << std::endl;
- }
- dataFile << "# color 255 0 0\n";
-
- curSat = analyticSolution_.AnalyticSolution()[0];
-
- dataFile << 0 << " " << curSat << " "
- << discretizationLength << " " << curSat <<std::endl;
- for (int i=1; i < cellNumberX; i++)
- {
- // Vertical Line
- dataFile << i*discretizationLength << " "
- << curSat << " ";
- curSat = analyticSolution_.AnalyticSolution()[i];
- dataFile << i*discretizationLength << " " << curSat << std::endl;
- //Horizontal Line
- dataFile << i*discretizationLength << " " << curSat << " "
- << (i+1)*discretizationLength << " " << curSat << std::endl;
- }
-// ElementIterator eItEnd = this->gridView().template end<0>();
-// for (ElementIterator eIt = this->gridView().template begin<0>(); eIt != eItEnd; ++eIt)
-// {
-// int currentIdx = this->variables().index(*eIt);
-// Scalar sat = this->variables().saturation()[currentIdx];
-// Scalar leftPos = eIt->geometry().corner(0)[0];
-// Scalar rightPos = eIt->geometry().corner(1)[0];
-// dataFile << leftPos <<" "<<sat<<" "<<rightPos<<" "<<sat<<"\n";
-// }
-
-
-// for (ElementIterator eIt = this->gridView().template begin<0>(); eIt != eItEnd; ++eIt)
-// {
-// int currentIdx = this->variables().index(*eIt);
-// Scalar sat = analyticSolution_.AnalyticSolution()[currentIdx];
-// Scalar leftPos = eIt->geometry().corner(0)[0];
-// Scalar rightPos = eIt->geometry().corner(1)[0];
-// dataFile << leftPos <<" "<<sat<<" "<<rightPos<<" "<<sat<<"\n";
-// }
-
- dataFile << "# legend Numerical Solution,Analytic Solution\n";
- dataFile.close();
+ int cellNumberX = static_cast<int>(300/discretizationLength);
+ discretizationLength = 300/cellNumberX; //Might be slightly different from the input parameter
+
+ std::ofstream dataFile;
+ dataFile.open("vipplot.vgf", std::fstream::app);
+
+ std::cout<<"Writing output for time step.\n";
+
+ if (time > 0)
+ dataFile << "# newframe\n";
+
+ dataFile << "# title Buckley-Leverett-Problem\n";
+ dataFile << "# x-label x\n";
+ dataFile << "# y-label Saturation\n";
+ dataFile << "# x-range -10 310\n";
+ dataFile << "# y-range 0 1\n";
+
+ dataFile << "# color 0 0 255\n";
+
+ Scalar curSat = this->variables().saturation()[0];
+
+ // Draw first piece separately, as no vertical line is needed
+ dataFile << 0 << " " << curSat << " "
+ << discretizationLength << " " << curSat <<std::endl;
+ for (int i=1; i < cellNumberX; i++)
+ {
+ // Vertical Line
+ dataFile << i*discretizationLength << " "
+ << curSat << " ";
+ curSat = this->variables().saturation()[i];
+ dataFile << i*discretizationLength << " " << curSat << std::endl;
+ //Horizontal Line
+ dataFile << i*discretizationLength << " " << curSat << " "
+ << (i+1)*discretizationLength << " " << curSat << std::endl;
+ }
+ dataFile << "# color 255 0 0\n";
+
+ curSat = analyticSolution_.AnalyticSolution()[0];
+
+ dataFile << 0 << " " << curSat << " "
+ << discretizationLength << " " << curSat <<std::endl;
+ for (int i=1; i < cellNumberX; i++)
+ {
+ // Vertical Line
+ dataFile << i*discretizationLength << " "
+ << curSat << " ";
+ curSat = analyticSolution_.AnalyticSolution()[i];
+ dataFile << i*discretizationLength << " " << curSat << std::endl;
+ //Horizontal Line
+ dataFile << i*discretizationLength << " " << curSat << " "
+ << (i+1)*discretizationLength << " " << curSat << std::endl;
+ }
+// ElementIterator eItEnd = this->gridView().template end<0>();
+// for (ElementIterator eIt = this->gridView().template begin<0>(); eIt != eItEnd; ++eIt)
+// {
+// int currentIdx = this->variables().index(*eIt);
+// Scalar sat = this->variables().saturation()[currentIdx];
+// Scalar leftPos = eIt->geometry().corner(0)[0];
+// Scalar rightPos = eIt->geometry().corner(1)[0];
+// dataFile << leftPos <<" "<<sat<<" "<<rightPos<<" "<<sat<<"\n";
+// }
+
+
+// for (ElementIterator eIt = this->gridView().template begin<0>(); eIt != eItEnd; ++eIt)
+// {
+// int currentIdx = this->variables().index(*eIt);
+// Scalar sat = analyticSolution_.AnalyticSolution()[currentIdx];
+// Scalar leftPos = eIt->geometry().corner(0)[0];
+// Scalar rightPos = eIt->geometry().corner(1)[0];
+// dataFile << leftPos <<" "<<sat<<" "<<rightPos<<" "<<sat<<"\n";
+// }
+
+ dataFile << "# legend Numerical Solution,Analytic Solution\n";
+ dataFile.close();
}
private:
diff --git a/appl/lecture/msm/mcwhorter/mcwhorter_analytic.hh b/appl/lecture/msm/mcwhorter/mcwhorter_analytic.hh
index 9a264c8b1ba21e339740334d3dd5b800a107d075..5dd6d6d148a47a486d656481600b01094226aae4 100644
--- a/appl/lecture/msm/mcwhorter/mcwhorter_analytic.hh
+++ b/appl/lecture/msm/mcwhorter/mcwhorter_analytic.hh
@@ -345,7 +345,7 @@ public:
BlockVector AnalyticSolution() const
{
- return analyticSolution_;
+ return analyticSolution_;
}
//Write saturation and pressure into file
diff --git a/appl/lecture/msm/mcwhorter/mcwhorterproblem.hh b/appl/lecture/msm/mcwhorter/mcwhorterproblem.hh
index de6af96249e9b5addc62a75b99e3dbed2691ffd0..81f1925a2a3d31ffb703505e1dcf601ad7587ec1 100644
--- a/appl/lecture/msm/mcwhorter/mcwhorterproblem.hh
+++ b/appl/lecture/msm/mcwhorter/mcwhorterproblem.hh
@@ -155,9 +155,9 @@ public:
this->setOutputInterval(10);
//Write header for ViPLab-Outputfile
- std::ofstream dataFile;
- dataFile.open("vipplot.vgf");
- dataFile.close();
+ std::ofstream dataFile;
+ dataFile.open("vipplot.vgf");
+ dataFile.close();
}
/*!
@@ -249,64 +249,64 @@ public:
//Override outputfunction for ViPLab-Output
void writeOutput()
{
- double time = this->timeManager().time();
- if (time < 0)
- return;
+ double time = this->timeManager().time();
+ if (time < 0)
+ return;
double discretizationLength = Params::tree().template get<double>("problem.DiscretizationLength");
- int cellNumberX = static_cast<int>(2/discretizationLength);
- discretizationLength = 2.0/cellNumberX; //Might be slightly different from the input parameter
-
- std::ofstream dataFile;
- dataFile.open("vipplot.vgf", std::fstream::app);
-
- std::cout<<"Writing output for time step.\n";
-
- if (time > 0)
- dataFile << "# newframe\n";
-
- dataFile << "# title Mc Whorther Problem\n";
- dataFile << "# x-label x\n";
- dataFile << "# y-label Saturation\n";
- dataFile << "# x-range -0.1 2.1\n";
- dataFile << "# y-range 0 1\n";
-
- dataFile << "# color 0 0 255\n";
-
- Scalar curSat = this->variables().saturation()[0];
-
- // Draw first piece separately, as no vertical line is needed
- dataFile << 0 << " " << curSat << " "
- << discretizationLength << " " << curSat <<std::endl;
- for (int i=1; i < cellNumberX; i++)
- {
- // Vertical Line
- dataFile << i*discretizationLength << " "
- << curSat << " ";
- curSat = this->variables().saturation()[i];
- dataFile << i*discretizationLength << " " << curSat << std::endl;
- //Horizontal Line
- dataFile << i*discretizationLength << " " << curSat << " "
- << (i+1)*discretizationLength << " " << curSat << std::endl;
- }
- dataFile << "# color 255 0 0\n";
-
- curSat = analyticSolution_.AnalyticSolution()[0];
- dataFile << 0 << " " << curSat << " "
- << discretizationLength << " " << curSat <<std::endl;
- for (int i=1; i < cellNumberX; i++)
- {
- // Vertical Line
- dataFile << i*discretizationLength << " "
- << curSat << " ";
- curSat = analyticSolution_.AnalyticSolution()[i];
- dataFile << i*discretizationLength << " " << curSat << std::endl;
- //Horizontal Line
- dataFile << i*discretizationLength << " " << curSat << " "
- << (i+1)*discretizationLength << " " << curSat << std::endl;
- }
- dataFile << "# legend Numerical Solution,Analytic Solution\n";
- dataFile.close();
+ int cellNumberX = static_cast<int>(2/discretizationLength);
+ discretizationLength = 2.0/cellNumberX; //Might be slightly different from the input parameter
+
+ std::ofstream dataFile;
+ dataFile.open("vipplot.vgf", std::fstream::app);
+
+ std::cout<<"Writing output for time step.\n";
+
+ if (time > 0)
+ dataFile << "# newframe\n";
+
+ dataFile << "# title Mc Whorther Problem\n";
+ dataFile << "# x-label x\n";
+ dataFile << "# y-label Saturation\n";
+ dataFile << "# x-range -0.1 2.1\n";
+ dataFile << "# y-range 0 1\n";
+
+ dataFile << "# color 0 0 255\n";
+
+ Scalar curSat = this->variables().saturation()[0];
+
+ // Draw first piece separately, as no vertical line is needed
+ dataFile << 0 << " " << curSat << " "
+ << discretizationLength << " " << curSat <<std::endl;
+ for (int i=1; i < cellNumberX; i++)
+ {
+ // Vertical Line
+ dataFile << i*discretizationLength << " "
+ << curSat << " ";
+ curSat = this->variables().saturation()[i];
+ dataFile << i*discretizationLength << " " << curSat << std::endl;
+ //Horizontal Line
+ dataFile << i*discretizationLength << " " << curSat << " "
+ << (i+1)*discretizationLength << " " << curSat << std::endl;
+ }
+ dataFile << "# color 255 0 0\n";
+
+ curSat = analyticSolution_.AnalyticSolution()[0];
+ dataFile << 0 << " " << curSat << " "
+ << discretizationLength << " " << curSat <<std::endl;
+ for (int i=1; i < cellNumberX; i++)
+ {
+ // Vertical Line
+ dataFile << i*discretizationLength << " "
+ << curSat << " ";
+ curSat = analyticSolution_.AnalyticSolution()[i];
+ dataFile << i*discretizationLength << " " << curSat << std::endl;
+ //Horizontal Line
+ dataFile << i*discretizationLength << " " << curSat << " "
+ << (i+1)*discretizationLength << " " << curSat << std::endl;
+ }
+ dataFile << "# legend Numerical Solution,Analytic Solution\n";
+ dataFile.close();
}
diff --git a/dumux/boxmodels/2p/2pmodel.hh b/dumux/boxmodels/2p/2pmodel.hh
index 79ed8bce58346b6926bbcf6c4794c1af3601924d..fb2676b09cd4794b11b7a9686a58163f32e0a10d 100644
--- a/dumux/boxmodels/2p/2pmodel.hh
+++ b/dumux/boxmodels/2p/2pmodel.hh
@@ -155,9 +155,9 @@ public:
// initialize velocity fields
for (int i = 0; i < numVertices; ++i)
{
- (*velocityN)[i] = Scalar(0);
- (*velocityW)[i] = Scalar(0);
- (*cellNum)[i] = Scalar(0.0);
+ (*velocityN)[i] = Scalar(0);
+ (*velocityW)[i] = Scalar(0);
+ (*cellNum)[i] = Scalar(0.0);
}
}
unsigned numElements = this->gridView_().size(0);
@@ -172,9 +172,9 @@ public:
for (; elemIt != elemEndIt; ++elemIt)
{
#warning "currently, velocity output only works for cubes and is set to false for simplices"
- if(elemIt->geometry().type().isCube() == false){
- velocityOutput = false;
- }
+ if(elemIt->geometry().type().isCube() == false){
+ velocityOutput = false;
+ }
int idx = this->elementMapper().map(*elemIt);
(*rank)[idx] = this->gridView_().comm().rank();
@@ -224,13 +224,13 @@ public:
scvVelocityW = 0;
scvVelocityN = 0;
- ElementVolumeVariables elemVolVars;
-
- elemVolVars.update(this->problem_(),
- *elemIt,
- fvElemGeom,
- false /* oldSol? */);
-
+ ElementVolumeVariables elemVolVars;
+
+ elemVolVars.update(this->problem_(),
+ *elemIt,
+ fvElemGeom,
+ false /* oldSol? */);
+
for (int faceIdx = 0; faceIdx< fvElemGeom.numEdges; faceIdx++)
{
@@ -243,7 +243,7 @@ public:
for (int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx)
{
- // data attached to upstream and the downstream vertices
+ // data attached to upstream and the downstream vertices
// of the current phase
const VolumeVariables up =
elemVolVars[fluxDat.upstreamIdx(phaseIdx)];
@@ -266,7 +266,7 @@ public:
GlobalPosition localNormal(0);
jacobianT1.mv(globalNormal, localNormal);
- // note only works for cubes
+ // note only works for cubes
const Scalar localArea = pow(2,-(dim-1));
localNormal /= localNormal.two_norm();
@@ -276,10 +276,10 @@ public:
const Scalar massUpwindWeight = GET_PARAM(TypeTag, Scalar, MassUpwindWeight);
PhasesVector q;
q[phaseIdx] = normalFlux
- * (massUpwindWeight
- * up.mobility(phaseIdx)
- + (1- massUpwindWeight)
- * dn.mobility(phaseIdx)) / localArea;
+ * (massUpwindWeight
+ * up.mobility(phaseIdx)
+ + (1- massUpwindWeight)
+ * dn.mobility(phaseIdx)) / localArea;
// transform the normal Darcy velocity into a vector
tmpVelocity[phaseIdx] = localNormal;
@@ -299,21 +299,21 @@ public:
const Dune::FieldVector<Scalar, dim> &localPos
= ReferenceElements::general(elemIt->geometry().type()).position(0, 0);
- // get the transposed Jacobian of the element mapping
+ // get the transposed Jacobian of the element mapping
const Dune::FieldMatrix<CoordScalar, dim, dim> &jacobianT2
= elemIt->geometry().jacobianTransposed(localPos);
// transform vertex velocities from local to global coordinates
- for (int i = 0; i < numVerts; ++i)
+ for (int i = 0; i < numVerts; ++i)
{
- int globalIdx = this->vertexMapper().map(*elemIt, i, dim);
+ int globalIdx = this->vertexMapper().map(*elemIt, i, dim);
// calculate the subcontrolvolume velocity by the Piola transformation
Dune::FieldVector<CoordScalar, dim> scvVelocity(0);
jacobianT2.mtv(scvVelocityW[i], scvVelocity);
scvVelocity /= elemIt->geometry().integrationElement(localPos);
// add up the wetting phase subcontrolvolume velocities for each vertex
- (*velocityW)[globalIdx] += scvVelocity;
+ (*velocityW)[globalIdx] += scvVelocity;
jacobianT2.mtv(scvVelocityN[i], scvVelocity);
scvVelocity /= elemIt->geometry().integrationElement(localPos);
@@ -324,11 +324,11 @@ public:
}
if(velocityOutput)
{
- // divide the vertex velocities by the number of adjacent scvs i.e. cells
- for(int globalIdx = 0; globalIdx<numVertices; ++globalIdx){
- (*velocityW)[globalIdx] /= (*cellNum)[globalIdx];
- (*velocityN)[globalIdx] /= (*cellNum)[globalIdx];
- }
+ // divide the vertex velocities by the number of adjacent scvs i.e. cells
+ for(int globalIdx = 0; globalIdx<numVertices; ++globalIdx){
+ (*velocityW)[globalIdx] /= (*cellNum)[globalIdx];
+ (*velocityN)[globalIdx] /= (*cellNum)[globalIdx];
+ }
}
writer.attachVertexData(*Sn, "Sn");
@@ -344,8 +344,8 @@ public:
writer.attachVertexData(*Te, "temperature");
if(velocityOutput) // check if velocity output is demanded
{
- writer.attachVertexData(*velocityW, "velocityW", dim);
- writer.attachVertexData(*velocityN, "velocityN", dim);
+ writer.attachVertexData(*velocityW, "velocityW", dim);
+ writer.attachVertexData(*velocityN, "velocityN", dim);
}
writer.attachCellData(*rank, "process rank");
}
diff --git a/dumux/boxmodels/2p2c/2p2cmodel.hh b/dumux/boxmodels/2p2c/2p2cmodel.hh
index 61da2b84f24d024abd2e7d6605d04527c5ec8366..55d2e600780872351eb70ba0f13f0ceeb8870278 100644
--- a/dumux/boxmodels/2p2c/2p2cmodel.hh
+++ b/dumux/boxmodels/2p2c/2p2cmodel.hh
@@ -324,9 +324,9 @@ public:
// initialize velocity fields
for (int i = 0; i < numVertices; ++i)
{
- (*velocityG)[i] = Scalar(0);
- (*velocityL)[i] = Scalar(0);
- (*cellNum)[i] = Scalar(0.0);
+ (*velocityG)[i] = Scalar(0);
+ (*velocityL)[i] = Scalar(0);
+ (*cellNum)[i] = Scalar(0.0);
}
}
@@ -400,13 +400,13 @@ public:
scvVelocityL = 0;
scvVelocityG = 0;
- ElementVolumeVariables elemVolVars;
-
- elemVolVars.update(this->problem_(),
- *elemIt,
- fvElemGeom,
- false /* oldSol? */);
-
+ ElementVolumeVariables elemVolVars;
+
+ elemVolVars.update(this->problem_(),
+ *elemIt,
+ fvElemGeom,
+ false /* oldSol? */);
+
for (int faceIdx = 0; faceIdx< fvElemGeom.numEdges; faceIdx++)
{
@@ -419,7 +419,7 @@ public:
for (int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx)
{
- // data attached to upstream and the downstream vertices
+ // data attached to upstream and the downstream vertices
// of the current phase
const VolumeVariables up =
elemVolVars[fluxDat.upstreamIdx(phaseIdx)];
@@ -435,7 +435,7 @@ public:
GlobalPosition localNormal(0);
jacobianT1.mv(globalNormal, localNormal);
- // note only works for cubes
+ // note only works for cubes
const Scalar localArea = pow(2,-(dim-1));
localNormal /= localNormal.two_norm();
@@ -445,10 +445,10 @@ public:
massUpwindWeight_ = GET_PARAM(TypeTag, Scalar, MassUpwindWeight);
PhasesVector q;
q[phaseIdx] = fluxDat.KmvpNormal(phaseIdx)
- * (massUpwindWeight_
- * up.mobility(phaseIdx)
- + (1- massUpwindWeight_)
- * dn.mobility(phaseIdx)) / localArea;
+ * (massUpwindWeight_
+ * up.mobility(phaseIdx)
+ + (1- massUpwindWeight_)
+ * dn.mobility(phaseIdx)) / localArea;
// transform the normal Darcy velocity into a vector
tmpVelocity[phaseIdx] = localNormal;
@@ -469,35 +469,35 @@ public:
const Dune::FieldVector<Scalar, dim>& localPos =
ReferenceElements::general(elemIt->geometry().type()).position(0, 0);
- // get the transposed Jacobian of the element mapping
+ // get the transposed Jacobian of the element mapping
const Dune::FieldMatrix<CoordScalar, dim, dim> jacobianT2 = elemIt->geometry().jacobianTransposed(localPos);
// transform vertex velocities from local to global coordinates
- for (int i = 0; i < numVerts; ++i)
+ for (int i = 0; i < numVerts; ++i)
{
- int globalIdx = this->vertexMapper().map(*elemIt, i, dim);
+ int globalIdx = this->vertexMapper().map(*elemIt, i, dim);
// calculate the subcontrolvolume velocity by the Piola transformation
Dune::FieldVector<CoordScalar, dim> scvVelocity(0);
jacobianT2.mtv(scvVelocityL[i], scvVelocity);
scvVelocity /= elemIt->geometry().integrationElement(localPos);
- // add up the wetting phase subcontrolvolume velocities for each vertex
- (*velocityL)[globalIdx] += scvVelocity;
+ // add up the wetting phase subcontrolvolume velocities for each vertex
+ (*velocityL)[globalIdx] += scvVelocity;
jacobianT2.mtv(scvVelocityG[i], scvVelocity);
scvVelocity /= elemIt->geometry().integrationElement(localPos);
- // add up the nonwetting phase subcontrolvolume velocities for each vertex
- (*velocityG)[globalIdx] += scvVelocity;
+ // add up the nonwetting phase subcontrolvolume velocities for each vertex
+ (*velocityG)[globalIdx] += scvVelocity;
}
}
}
if(velocityOutput_)
{
// divide the vertex velocities by the number of adjacent scvs i.e. cells
- for(int globalIdx = 0; globalIdx<numVertices; ++globalIdx){
- (*velocityL)[globalIdx] /= (*cellNum)[globalIdx];
- (*velocityG)[globalIdx] /= (*cellNum)[globalIdx];
- }
+ for(int globalIdx = 0; globalIdx<numVertices; ++globalIdx){
+ (*velocityL)[globalIdx] /= (*cellNum)[globalIdx];
+ (*velocityG)[globalIdx] /= (*cellNum)[globalIdx];
+ }
}
@@ -526,8 +526,8 @@ public:
if(velocityOutput_) // check if velocity output is demanded
{
- writer.attachVertexData(*velocityL, "velocityL", dim);
- writer.attachVertexData(*velocityG, "velocityG", dim);
+ writer.attachVertexData(*velocityL, "velocityL", dim);
+ writer.attachVertexData(*velocityG, "velocityG", dim);
}
writer.attachCellData(*rank, "process rank");
}
diff --git a/dumux/decoupled/2p/diffusion/fv/fvpressure2padaptive.hh b/dumux/decoupled/2p/diffusion/fv/fvpressure2padaptive.hh
index 70ef1bbb48a10f4f55accd8f0453220290370f88..8f667d3e3fa4858824c54f7c63698b2539303ea7 100644
--- a/dumux/decoupled/2p/diffusion/fv/fvpressure2padaptive.hh
+++ b/dumux/decoupled/2p/diffusion/fv/fvpressure2padaptive.hh
@@ -342,12 +342,12 @@ void FVPressure2Padaptive<TypeTag>::initializeMatrix()
template<class TypeTag>
void FVPressure2Padaptive<TypeTag>::assemble(bool first)
{
- updateMaterialLaws();
- // Adapt size of matrix A_ and Vector f_
- {int size = problem_.variables().gridSize();
- A_.setSize(size,size);
- f_.resize(size);}
- initializeMatrix(); // define non-zero entries in matrix
+ updateMaterialLaws();
+ // Adapt size of matrix A_ and Vector f_
+ {int size = problem_.variables().gridSize();
+ A_.setSize(size,size);
+ f_.resize(size);}
+ initializeMatrix(); // define non-zero entries in matrix
// initialization: set matrix A_ to zero
A_ = 0;
@@ -414,408 +414,408 @@ void FVPressure2Padaptive<TypeTag>::assemble(bool first)
// "normal" case: neighbor has same level
if (neighborPointer->level()==eIt.level())
{
- // neighbor cell center in global coordinates
- const GlobalPosition& globalPosNeighbor = neighborPointer->geometry().center();
-
- // distance vector between barycenters
- GlobalPosition distVec = globalPosNeighbor - globalPos;
-
- // compute distance between cell centers
-// Scalar dist = distVec.two_norm();
- Scalar dist = distVec*unitOuterNormal;
-
- // compute vectorized permeabilities
- FieldMatrix meanPermeability(0);
-
- problem_.spatialParameters().meanK(meanPermeability,
- problem_.spatialParameters().intrinsicPermeability(*eIt),
- problem_.spatialParameters().intrinsicPermeability(*neighborPointer));
-
- Dune::FieldVector<Scalar, dim> permeability(0);
- meanPermeability.mv(unitOuterNormal, permeability);
-
- // get mobilities and fractional flow factors
- Scalar lambdaWJ = problem_.variables().mobilityWetting(globalIdxJ);
- Scalar lambdaNWJ = problem_.variables().mobilityNonwetting(globalIdxJ);
- Scalar fractionalWJ = problem_.variables().fracFlowFuncWetting(globalIdxJ);
- Scalar fractionalNWJ = problem_.variables().fracFlowFuncNonwetting(globalIdxJ);
- Scalar densityWJ = problem_.variables().densityWetting(globalIdxJ);
- Scalar densityNWJ = problem_.variables().densityNonwetting(globalIdxJ);
-
- Scalar pcJ = problem_.variables().capillaryPressure(globalIdxJ);
-
- Scalar rhoMeanW = 0.5 * (densityWI + densityWJ);
- Scalar rhoMeanNW = 0.5 * (densityNWI + densityNWJ);
- Scalar fMeanW = 0.5 * (fractionalWI + fractionalWJ);
- Scalar fMeanNW = 0.5 * (fractionalNWI + fractionalNWJ);
-
- // update diagonal entry
- Scalar entry;
-
- //calculate potential gradients
- Scalar potentialW = 0;
- Scalar potentialNW = 0;
-
- Scalar densityW = 0;
- Scalar densityNW = 0;
-
- //if we are at the very first iteration we can't calculate phase potentials
- if (!first)
- {
- // After grid adaption, the potentials from the previous time step
- // can not be used for upwinding. The mean values is assumed instead.
- densityW = rhoMeanW;
- densityNW = rhoMeanNW;
-
- switch (pressureType)
- {
- case pw:
- {
- potentialW = (problem_.variables().pressure()[globalIdxI] - problem_.variables().pressure()[globalIdxJ]);
- potentialNW = (problem_.variables().pressure()[globalIdxI] - problem_.variables().pressure()[globalIdxJ] + pcI
- - pcJ);
- break;
- }
- case pn:
- {
- potentialW
- = (problem_.variables().pressure()[globalIdxI] - problem_.variables().pressure()[globalIdxJ] - pcI + pcJ);
- potentialNW = (problem_.variables().pressure()[globalIdxI] - problem_.variables().pressure()[globalIdxJ]);
- break;
- }
- case pglobal:
- {
- potentialW = (problem_.variables().pressure()[globalIdxI] - problem_.variables().pressure()[globalIdxJ] - fMeanNW
- * (pcI - pcJ));
- potentialNW = (problem_.variables().pressure()[globalIdxI] - problem_.variables().pressure()[globalIdxJ] + fMeanW
- * (pcI - pcJ));
- break;
- }
- }
-
- potentialW += densityW * (distVec * gravity);
- potentialNW += densityNW * (distVec * gravity);
-
- //store potentials for further calculations (velocity, saturation, ...)
- problem_.variables().potentialWetting(globalIdxI, isIndex) = potentialW;
- problem_.variables().potentialNonwetting(globalIdxI, isIndex) = potentialNW;
- }
-
- //do the upwinding of the mobility depending on the phase potentials
- Scalar lambdaW = (potentialW > 0.) ? lambdaWI : lambdaWJ;
- lambdaW = (potentialW == 0) ? 0.5 * (lambdaWI + lambdaWJ) : lambdaW;
- Scalar lambdaNW = (potentialNW > 0) ? lambdaNWI : lambdaNWJ;
- lambdaNW = (potentialNW == 0) ? 0.5 * (lambdaNWI + lambdaNWJ) : lambdaNW;
-
- densityW = (potentialW > 0.) ? densityWI : densityWJ;
- densityNW = (potentialNW > 0.) ? densityNWI : densityNWJ;
-
- densityW = (potentialW == 0) ? rhoMeanW : densityW;
- densityNW = (potentialNW == 0) ? rhoMeanNW : densityNW;
-
- //calculate current matrix entry
- entry = (lambdaW + lambdaNW) * ((permeability * unitOuterNormal) / dist) * faceArea;
-
- //calculate right hand side
- Scalar rightEntry = (lambdaW * densityW + lambdaNW * densityNW) * (permeability * gravity) * faceArea;
-
- switch (pressureType)
- {
- case pw:
- {
- // calculate capillary pressure gradient
- Dune::FieldVector<Scalar, dim> pCGradient = unitOuterNormal;
- pCGradient *= (pcI - pcJ) / dist;
-
- //add capillary pressure term to right hand side
- rightEntry += 0.5 * (lambdaNWI + lambdaNWJ) * (permeability * pCGradient) * faceArea;
- break;
- }
- case pn:
- {
- // calculate capillary pressure gradient
- Dune::FieldVector<Scalar, dim> pCGradient = unitOuterNormal;
- pCGradient *= (pcI - pcJ) / dist;
-
- //add capillary pressure term to right hand side
- rightEntry -= 0.5 * (lambdaWI + lambdaWJ) * (permeability * pCGradient) * faceArea;
- break;
- }
- }
-
- //set right hand side
- f_[globalIdxI] -= rightEntry;
-
- // set diagonal entry
- A_[globalIdxI][globalIdxI] += entry;
-
- // set off-diagonal entry
- A_[globalIdxI][globalIdxJ] -= entry;
+ // neighbor cell center in global coordinates
+ const GlobalPosition& globalPosNeighbor = neighborPointer->geometry().center();
+
+ // distance vector between barycenters
+ GlobalPosition distVec = globalPosNeighbor - globalPos;
+
+ // compute distance between cell centers
+// Scalar dist = distVec.two_norm();
+ Scalar dist = distVec*unitOuterNormal;
+
+ // compute vectorized permeabilities
+ FieldMatrix meanPermeability(0);
+
+ problem_.spatialParameters().meanK(meanPermeability,
+ problem_.spatialParameters().intrinsicPermeability(*eIt),
+ problem_.spatialParameters().intrinsicPermeability(*neighborPointer));
+
+ Dune::FieldVector<Scalar, dim> permeability(0);
+ meanPermeability.mv(unitOuterNormal, permeability);
+
+ // get mobilities and fractional flow factors
+ Scalar lambdaWJ = problem_.variables().mobilityWetting(globalIdxJ);
+ Scalar lambdaNWJ = problem_.variables().mobilityNonwetting(globalIdxJ);
+ Scalar fractionalWJ = problem_.variables().fracFlowFuncWetting(globalIdxJ);
+ Scalar fractionalNWJ = problem_.variables().fracFlowFuncNonwetting(globalIdxJ);
+ Scalar densityWJ = problem_.variables().densityWetting(globalIdxJ);
+ Scalar densityNWJ = problem_.variables().densityNonwetting(globalIdxJ);
+
+ Scalar pcJ = problem_.variables().capillaryPressure(globalIdxJ);
+
+ Scalar rhoMeanW = 0.5 * (densityWI + densityWJ);
+ Scalar rhoMeanNW = 0.5 * (densityNWI + densityNWJ);
+ Scalar fMeanW = 0.5 * (fractionalWI + fractionalWJ);
+ Scalar fMeanNW = 0.5 * (fractionalNWI + fractionalNWJ);
+
+ // update diagonal entry
+ Scalar entry;
+
+ //calculate potential gradients
+ Scalar potentialW = 0;
+ Scalar potentialNW = 0;
+
+ Scalar densityW = 0;
+ Scalar densityNW = 0;
+
+ //if we are at the very first iteration we can't calculate phase potentials
+ if (!first)
+ {
+ // After grid adaption, the potentials from the previous time step
+ // can not be used for upwinding. The mean values is assumed instead.
+ densityW = rhoMeanW;
+ densityNW = rhoMeanNW;
+
+ switch (pressureType)
+ {
+ case pw:
+ {
+ potentialW = (problem_.variables().pressure()[globalIdxI] - problem_.variables().pressure()[globalIdxJ]);
+ potentialNW = (problem_.variables().pressure()[globalIdxI] - problem_.variables().pressure()[globalIdxJ] + pcI
+ - pcJ);
+ break;
+ }
+ case pn:
+ {
+ potentialW
+ = (problem_.variables().pressure()[globalIdxI] - problem_.variables().pressure()[globalIdxJ] - pcI + pcJ);
+ potentialNW = (problem_.variables().pressure()[globalIdxI] - problem_.variables().pressure()[globalIdxJ]);
+ break;
+ }
+ case pglobal:
+ {
+ potentialW = (problem_.variables().pressure()[globalIdxI] - problem_.variables().pressure()[globalIdxJ] - fMeanNW
+ * (pcI - pcJ));
+ potentialNW = (problem_.variables().pressure()[globalIdxI] - problem_.variables().pressure()[globalIdxJ] + fMeanW
+ * (pcI - pcJ));
+ break;
+ }
+ }
+
+ potentialW += densityW * (distVec * gravity);
+ potentialNW += densityNW * (distVec * gravity);
+
+ //store potentials for further calculations (velocity, saturation, ...)
+ problem_.variables().potentialWetting(globalIdxI, isIndex) = potentialW;
+ problem_.variables().potentialNonwetting(globalIdxI, isIndex) = potentialNW;
+ }
+
+ //do the upwinding of the mobility depending on the phase potentials
+ Scalar lambdaW = (potentialW > 0.) ? lambdaWI : lambdaWJ;
+ lambdaW = (potentialW == 0) ? 0.5 * (lambdaWI + lambdaWJ) : lambdaW;
+ Scalar lambdaNW = (potentialNW > 0) ? lambdaNWI : lambdaNWJ;
+ lambdaNW = (potentialNW == 0) ? 0.5 * (lambdaNWI + lambdaNWJ) : lambdaNW;
+
+ densityW = (potentialW > 0.) ? densityWI : densityWJ;
+ densityNW = (potentialNW > 0.) ? densityNWI : densityNWJ;
+
+ densityW = (potentialW == 0) ? rhoMeanW : densityW;
+ densityNW = (potentialNW == 0) ? rhoMeanNW : densityNW;
+
+ //calculate current matrix entry
+ entry = (lambdaW + lambdaNW) * ((permeability * unitOuterNormal) / dist) * faceArea;
+
+ //calculate right hand side
+ Scalar rightEntry = (lambdaW * densityW + lambdaNW * densityNW) * (permeability * gravity) * faceArea;
+
+ switch (pressureType)
+ {
+ case pw:
+ {
+ // calculate capillary pressure gradient
+ Dune::FieldVector<Scalar, dim> pCGradient = unitOuterNormal;
+ pCGradient *= (pcI - pcJ) / dist;
+
+ //add capillary pressure term to right hand side
+ rightEntry += 0.5 * (lambdaNWI + lambdaNWJ) * (permeability * pCGradient) * faceArea;
+ break;
+ }
+ case pn:
+ {
+ // calculate capillary pressure gradient
+ Dune::FieldVector<Scalar, dim> pCGradient = unitOuterNormal;
+ pCGradient *= (pcI - pcJ) / dist;
+
+ //add capillary pressure term to right hand side
+ rightEntry -= 0.5 * (lambdaWI + lambdaWJ) * (permeability * pCGradient) * faceArea;
+ break;
+ }
+ }
+
+ //set right hand side
+ f_[globalIdxI] -= rightEntry;
+
+ // set diagonal entry
+ A_[globalIdxI][globalIdxI] += entry;
+
+ // set off-diagonal entry
+ A_[globalIdxI][globalIdxJ] -= entry;
} // end of case "neighbor has same level"
// hanging node situation: neighbor has higher level
if (neighborPointer->level()==eIt.level()+1)
{
- // Count number of hanging nodes
- // nodecount counts each hanging node twice, we have to divide by 2 in the end
- // not really necessary
- nodecount++;
-
- int globalIdxK = 0;
- ElementPointer thirdCellPointer = isIt->outside();
- bool foundK=false;
- bool foundIJ=false;
- // We are looking for two things:
- // IsIndexJ, the index of the interface from the neighbor-cell point of view
- // GlobalIdxK, the index of the third cell
- // for efficienty this is done in one IntersectionIterator-Loop
-
- // Intersectioniterator around cell J
- int isIndexJ = -1;
- IntersectionIterator isItEndJ = this->problem().gridView().iend(*neighborPointer);
-
- for (IntersectionIterator isItJ = this->problem().gridView().ibegin(*neighborPointer); isItJ != isItEndJ; ++isItJ)
- {
- // increase isIndexJ, if it is not found yet
- if (!foundIJ)
- isIndexJ++;
- if (isItJ->neighbor())
- {
- ElementPointer neighborPointer2 = isItJ->outside();
-
- // Neighbor of neighbor is Cell I?
- if (this->problem().variables().index(*neighborPointer2)==globalIdxI)
- {
- foundIJ=true;
- }
- else
- {
- if (neighborPointer2->level()==eIt.level()+1)
- {
- // To verify, we found the correct Cell K, we check
- // - is level(K)=level(J)?
- // - is distance(IJ)=distance(IK)?
- // - K is neighbor of J.
- const GlobalPosition& globalPosNeighbor = neighborPointer->geometry().center();
- const GlobalPosition& globalPosThirdCell = neighborPointer2->geometry().center();
- Dune::FieldVector<Scalar, dimWorld> distVecIJ = globalPosNeighbor - globalPos;
- Dune::FieldVector<Scalar, dimWorld> distVecIK = globalPosThirdCell - globalPos;
- if (((distVecIK.two_norm()-distVecIJ.two_norm()))/distVecIJ.two_norm() < 0.001)
- {
- globalIdxK= this->problem().variables().index(*neighborPointer2);
- thirdCellPointer = neighborPointer2;
- foundK=true;
- }
-
- }
- }
- }
- if (foundIJ && foundK) break;
- }
-
- // neighbor cell center in global coordinates
- const GlobalPosition& globalPosNeighbor = neighborPointer->geometry().center();
-// const GlobalPosition& globalPosThirdCell = thirdCellPointer->geometry().center();
- const GlobalPosition& globalPosInterface = isIt->geometry().center();
-
- Dune::FieldVector<Scalar,dimWorld> distVec = globalPosInterface - globalPos;
- Scalar lI= distVec*unitOuterNormal;
- distVec = globalPosNeighbor - globalPosInterface;
- Scalar lJ= distVec*unitOuterNormal;
- Scalar l=lI+lJ;
-
- FieldMatrix permeabilityI(0);
- FieldMatrix permeabilityJ(0);
- FieldMatrix permeabilityK(0);
-
- problem_.spatialParameters().meanK(permeabilityI,
- problem_.spatialParameters().intrinsicPermeability(*eIt));
- problem_.spatialParameters().meanK(permeabilityJ,
- problem_.spatialParameters().intrinsicPermeability(*neighborPointer));
- problem_.spatialParameters().meanK(permeabilityK,
- problem_.spatialParameters().intrinsicPermeability(*thirdCellPointer));
-
- // Calculate permeablity component normal to interface
- Scalar kI, kJ, kK, kMean, ng;
- Dune::FieldVector<Scalar, dim> permI(0);
- Dune::FieldVector<Scalar, dim> permJ(0);
- Dune::FieldVector<Scalar, dim> permK(0);
-
- permeabilityI.mv(unitOuterNormal, permI);
- permeabilityJ.mv(unitOuterNormal, permJ);
- permeabilityK.mv(unitOuterNormal, permK);
-
- // kI,kJ,kK=(n^T)Kn
- kI=unitOuterNormal*permI;
- kJ=unitOuterNormal*permJ;
- kK=unitOuterNormal*permK;
-
- // See Diplomarbeit Michael Sinsbeck
- kMean=kI*kJ*kK*l/(kJ*kK*lI+kI*(kJ+kK)/2*lJ);
-
- ng=this->gravity*unitOuterNormal;
-
- // get mobilities
- // fractional flow function is not evaluated, since we do not use p_global
- Scalar lambdaWJ = problem_.variables().mobilityWetting(globalIdxJ);
- Scalar lambdaNWJ = problem_.variables().mobilityNonwetting(globalIdxJ);
- Scalar densityWJ = problem_.variables().densityWetting(globalIdxJ);
- Scalar densityNWJ = problem_.variables().densityNonwetting(globalIdxJ);
- Scalar fractionalWJ = problem_.variables().fracFlowFuncWetting(globalIdxJ);
- Scalar fractionalNWJ = problem_.variables().fracFlowFuncNonwetting(globalIdxJ);
-
-// Scalar lambdaWK = problem_.variables().mobilityWetting(globalIdxK);
-// Scalar lambdaNWK = problem_.variables().mobilityNonwetting(globalIdxK);
- Scalar densityWK = problem_.variables().densityWetting(globalIdxK);
- Scalar densityNWK = problem_.variables().densityNonwetting(globalIdxK);
- Scalar fractionalWK = problem_.variables().fracFlowFuncWetting(globalIdxK);
- Scalar fractionalNWK = problem_.variables().fracFlowFuncNonwetting(globalIdxK);
-
- Scalar pcJ = problem_.variables().capillaryPressure(globalIdxJ);
- Scalar pcK = problem_.variables().capillaryPressure(globalIdxK);
- Scalar pcJK=(pcJ+pcK)/2;
-
- // Potentials from previous time step are not available.
- // Instead, calculate mean density, then find potentials,
- // then upwind density.
- // pressure from previous time step might also be incorrect.
-
- Scalar rhoMeanWIJ = (lJ*densityWI + lI*densityWJ)/l;
- Scalar rhoMeanNWIJ = (lJ*densityNWI + lI*densityNWJ)/l;
- Scalar rhoMeanWIK = (lJ*densityWI + lI*densityWK)/l;
- Scalar rhoMeanNWIK = (lJ*densityNWI + lI*densityNWK)/l;
-
- Scalar densityWIJ = 0;
- Scalar densityNWIJ = 0;
- Scalar densityWIK = 0;
- Scalar densityNWIK = 0;
-
- Scalar fMeanWIJ = (lJ*fractionalWI+lI*fractionalWJ)/l;
- Scalar fMeanNWIJ = (lJ*fractionalNWI+lI*fractionalNWJ)/l;
- Scalar fMeanWIK = (lJ*fractionalWI+lI*fractionalWK)/l;
- Scalar fMeanNWIK = (lJ*fractionalNWI+lI*fractionalNWK)/l;
-
- Scalar potentialWIJ = 0;
- Scalar potentialNWIJ = 0;
- Scalar potentialWIK = 0;
- Scalar potentialNWIK = 0;
-
- //if we are at the very first iteration we can't calculate phase potentials
- if (!first)
- {
- // potentials from previous time step no available.
- densityWIJ = rhoMeanWIJ;
- densityNWIJ = rhoMeanNWIJ;
- densityWIK = rhoMeanWIK;
- densityNWIK = rhoMeanNWIK;
-
- // Old pressure field
- Scalar pressI=problem_.variables().pressure()[globalIdxI];
- Scalar pressJ=problem_.variables().pressure()[globalIdxJ];
- Scalar pressK=problem_.variables().pressure()[globalIdxK];
- Scalar pressJK=(pressJ+pressK)/2;
- Scalar pcJK=(pcJ+pcK)/2;
-
- switch (pressureType)
- {
- case pw:
- {
- potentialWIJ = (pressI-pressJK)/l+
- (densityWIJ-lJ/l*(kI+kK)/kI*(densityWIK-densityWIJ)/2)*ng;
- potentialNWIJ = (pressI+pcI-(pressJK+pcJK))/l+
- (densityNWIJ-lJ/l*(kI+kK)/kI*(densityNWIK-densityNWIJ)/2)*ng;
- potentialWIK = (pressI-pressJK)/l+
- (densityWIK-lJ/l*(kI+kK)/kI*(densityWIJ-densityWIK)/2)*ng;
- potentialNWIK = (pressI+pcI-(pressJK+pcJK))/l+
- (densityNWIK-lJ/l*(kI+kK)/kI*(densityNWIJ-densityNWIK)/2)*ng;
- break;
- }
- case pn:
- {
- potentialWIJ = (pressI-pcI-(pressJK-pcJK))/l+
- (densityWIJ-lJ/l*(kI+kK)/kI*(densityWIK-densityWIJ)/2)*ng;
- potentialNWIJ = (pressI-pressJK)/l+
- (densityNWIJ-lJ/l*(kI+kK)/kI*(densityNWIK-densityNWIJ)/2)*ng;
- potentialWIK = (pressI-pcI-(pressJK-pcJK))/l+
- (densityWIK-lJ/l*(kI+kK)/kI*(densityWIJ-densityWIK)/2)*ng;
- potentialNWIK = (pressI-pressJK)/l+
- (densityNWIK-lJ/l*(kI+kK)/kI*(densityNWIJ-densityNWIK)/2)*ng;
- break;
- }
- case pglobal:
- {
- potentialWIJ = (pressI-fMeanNWIJ*pcI-(pressJK-fMeanNWIJ*pcJK))/l+
- (densityWIJ-lJ/l*(kI+kK)/kI*(densityWIK-densityWIJ)/2)*ng;
- potentialNWIJ = (pressI+fMeanWIJ*pcI-(pressJK+fMeanWIJ*pcJK))/l+
- (densityNWIJ-lJ/l*(kI+kK)/kI*(densityNWIK-densityNWIJ)/2)*ng;
- potentialWIK = (pressI-fMeanNWIK*pcI-(pressJK-fMeanNWIK*pcJK))/l+
- (densityWIK-lJ/l*(kI+kK)/kI*(densityWIJ-densityWIK)/2)*ng;
- potentialNWIK = (pressI+fMeanWIK*pcI-(pressJK+fMeanWIK*pcJK))/l+
- (densityNWIK-lJ/l*(kI+kK)/kI*(densityNWIJ-densityNWIK)/2)*ng;
- break;
- }
- }
-
- //store potentials for further calculations (velocity, saturation, ...)
-
- problem_.variables().potentialWetting(globalIdxI, isIndex) = potentialWIJ;
- problem_.variables().potentialNonwetting(globalIdxI, isIndex) = potentialNWIJ;
- problem_.variables().potentialWetting(globalIdxJ, isIndexJ) = -potentialWIJ;
- problem_.variables().potentialNonwetting(globalIdxJ, isIndexJ) = -potentialNWIJ;
- }
-
- //do the upwinding of the mobility depending on the phase potentials
- Scalar lambdaWIJ = (potentialWIJ > 0.) ? lambdaWI : lambdaWJ;
- lambdaWIJ = (potentialWIJ == 0) ? 0.5 * (lambdaWI + lambdaWJ) : lambdaWIJ;
- Scalar lambdaNWIJ = (potentialNWIJ > 0) ? lambdaNWI : lambdaNWJ;
- lambdaNWIJ = (potentialNWIJ == 0) ? 0.5 * (lambdaNWI + lambdaNWJ) : lambdaNWIJ;
- densityWIJ = (potentialWIJ > 0.) ? densityWI : densityWJ;
- densityNWIJ = (potentialNWIJ > 0.) ? densityNWI : densityNWJ;
- densityWIJ = (potentialWIJ == 0) ? rhoMeanWIJ : densityWIJ;
- densityNWIJ = (potentialNWIJ == 0) ? rhoMeanNWIJ : densityNWIJ;
-
- densityWIK = (potentialWIK > 0.) ? densityWI : densityWK;
- densityNWIK = (potentialNWIK > 0.) ? densityNWI : densityNWK;
- densityWIK = (potentialWIK == 0) ? rhoMeanWIK : densityWIK;
- densityNWIK = (potentialNWIK == 0) ? rhoMeanNWIK : densityNWIK;
-
-
- // update diagonal entry and right hand side
-
- Scalar entry;
- Scalar rightEntry;
-
- entry = (lambdaWIJ + lambdaNWIJ)*kMean/l*faceArea;
- rightEntry = faceArea*lambdaWIJ*kMean*ng*((densityWIJ)-(lJ/l)*(kI+kK)/kI*(densityWIK-densityWIJ)/2);
- rightEntry += faceArea*lambdaNWIJ*kMean*ng*((densityNWIJ)-(lJ/l)*(kI+kK)/kI*(densityNWIK-densityNWIJ)/2);
-
- switch (pressureType)
- {
- case pw:
- {
- rightEntry += faceArea*lambdaNWIJ*kMean*(pcJK-pcI)/l;
- break;
- }
- case pn:
- {
- rightEntry -= faceArea*lambdaWIJ*kMean*(pcJK-pcI)/l;
- break;
- }
- }
-
- f_[globalIdxI] -= rightEntry;
- f_[globalIdxJ] += rightEntry;
-
- // set diagonal entry
- A_[globalIdxI][globalIdxI] += entry;
- A_[globalIdxI][globalIdxJ] -= 0.5*entry;
- A_[globalIdxI][globalIdxK] -= 0.5*entry;
-
- // set off-diagonal entry
- A_[globalIdxJ][globalIdxI] -= entry;
- A_[globalIdxJ][globalIdxJ] += 0.5*entry;
- A_[globalIdxJ][globalIdxK] += 0.5*entry;
+ // Count number of hanging nodes
+ // nodecount counts each hanging node twice, we have to divide by 2 in the end
+ // not really necessary
+ nodecount++;
+
+ int globalIdxK = 0;
+ ElementPointer thirdCellPointer = isIt->outside();
+ bool foundK=false;
+ bool foundIJ=false;
+ // We are looking for two things:
+ // IsIndexJ, the index of the interface from the neighbor-cell point of view
+ // GlobalIdxK, the index of the third cell
+ // for efficienty this is done in one IntersectionIterator-Loop
+
+ // Intersectioniterator around cell J
+ int isIndexJ = -1;
+ IntersectionIterator isItEndJ = this->problem().gridView().iend(*neighborPointer);
+
+ for (IntersectionIterator isItJ = this->problem().gridView().ibegin(*neighborPointer); isItJ != isItEndJ; ++isItJ)
+ {
+ // increase isIndexJ, if it is not found yet
+ if (!foundIJ)
+ isIndexJ++;
+ if (isItJ->neighbor())
+ {
+ ElementPointer neighborPointer2 = isItJ->outside();
+
+ // Neighbor of neighbor is Cell I?
+ if (this->problem().variables().index(*neighborPointer2)==globalIdxI)
+ {
+ foundIJ=true;
+ }
+ else
+ {
+ if (neighborPointer2->level()==eIt.level()+1)
+ {
+ // To verify, we found the correct Cell K, we check
+ // - is level(K)=level(J)?
+ // - is distance(IJ)=distance(IK)?
+ // - K is neighbor of J.
+ const GlobalPosition& globalPosNeighbor = neighborPointer->geometry().center();
+ const GlobalPosition& globalPosThirdCell = neighborPointer2->geometry().center();
+ Dune::FieldVector<Scalar, dimWorld> distVecIJ = globalPosNeighbor - globalPos;
+ Dune::FieldVector<Scalar, dimWorld> distVecIK = globalPosThirdCell - globalPos;
+ if (((distVecIK.two_norm()-distVecIJ.two_norm()))/distVecIJ.two_norm() < 0.001)
+ {
+ globalIdxK= this->problem().variables().index(*neighborPointer2);
+ thirdCellPointer = neighborPointer2;
+ foundK=true;
+ }
+
+ }
+ }
+ }
+ if (foundIJ && foundK) break;
+ }
+
+ // neighbor cell center in global coordinates
+ const GlobalPosition& globalPosNeighbor = neighborPointer->geometry().center();
+// const GlobalPosition& globalPosThirdCell = thirdCellPointer->geometry().center();
+ const GlobalPosition& globalPosInterface = isIt->geometry().center();
+
+ Dune::FieldVector<Scalar,dimWorld> distVec = globalPosInterface - globalPos;
+ Scalar lI= distVec*unitOuterNormal;
+ distVec = globalPosNeighbor - globalPosInterface;
+ Scalar lJ= distVec*unitOuterNormal;
+ Scalar l=lI+lJ;
+
+ FieldMatrix permeabilityI(0);
+ FieldMatrix permeabilityJ(0);
+ FieldMatrix permeabilityK(0);
+
+ problem_.spatialParameters().meanK(permeabilityI,
+ problem_.spatialParameters().intrinsicPermeability(*eIt));
+ problem_.spatialParameters().meanK(permeabilityJ,
+ problem_.spatialParameters().intrinsicPermeability(*neighborPointer));
+ problem_.spatialParameters().meanK(permeabilityK,
+ problem_.spatialParameters().intrinsicPermeability(*thirdCellPointer));
+
+ // Calculate permeablity component normal to interface
+ Scalar kI, kJ, kK, kMean, ng;
+ Dune::FieldVector<Scalar, dim> permI(0);
+ Dune::FieldVector<Scalar, dim> permJ(0);
+ Dune::FieldVector<Scalar, dim> permK(0);
+
+ permeabilityI.mv(unitOuterNormal, permI);
+ permeabilityJ.mv(unitOuterNormal, permJ);
+ permeabilityK.mv(unitOuterNormal, permK);
+
+ // kI,kJ,kK=(n^T)Kn
+ kI=unitOuterNormal*permI;
+ kJ=unitOuterNormal*permJ;
+ kK=unitOuterNormal*permK;
+
+ // See Diplomarbeit Michael Sinsbeck
+ kMean=kI*kJ*kK*l/(kJ*kK*lI+kI*(kJ+kK)/2*lJ);
+
+ ng=this->gravity*unitOuterNormal;
+
+ // get mobilities
+ // fractional flow function is not evaluated, since we do not use p_global
+ Scalar lambdaWJ = problem_.variables().mobilityWetting(globalIdxJ);
+ Scalar lambdaNWJ = problem_.variables().mobilityNonwetting(globalIdxJ);
+ Scalar densityWJ = problem_.variables().densityWetting(globalIdxJ);
+ Scalar densityNWJ = problem_.variables().densityNonwetting(globalIdxJ);
+ Scalar fractionalWJ = problem_.variables().fracFlowFuncWetting(globalIdxJ);
+ Scalar fractionalNWJ = problem_.variables().fracFlowFuncNonwetting(globalIdxJ);
+
+// Scalar lambdaWK = problem_.variables().mobilityWetting(globalIdxK);
+// Scalar lambdaNWK = problem_.variables().mobilityNonwetting(globalIdxK);
+ Scalar densityWK = problem_.variables().densityWetting(globalIdxK);
+ Scalar densityNWK = problem_.variables().densityNonwetting(globalIdxK);
+ Scalar fractionalWK = problem_.variables().fracFlowFuncWetting(globalIdxK);
+ Scalar fractionalNWK = problem_.variables().fracFlowFuncNonwetting(globalIdxK);
+
+ Scalar pcJ = problem_.variables().capillaryPressure(globalIdxJ);
+ Scalar pcK = problem_.variables().capillaryPressure(globalIdxK);
+ Scalar pcJK=(pcJ+pcK)/2;
+
+ // Potentials from previous time step are not available.
+ // Instead, calculate mean density, then find potentials,
+ // then upwind density.
+ // pressure from previous time step might also be incorrect.
+
+ Scalar rhoMeanWIJ = (lJ*densityWI + lI*densityWJ)/l;
+ Scalar rhoMeanNWIJ = (lJ*densityNWI + lI*densityNWJ)/l;
+ Scalar rhoMeanWIK = (lJ*densityWI + lI*densityWK)/l;
+ Scalar rhoMeanNWIK = (lJ*densityNWI + lI*densityNWK)/l;
+
+ Scalar densityWIJ = 0;
+ Scalar densityNWIJ = 0;
+ Scalar densityWIK = 0;
+ Scalar densityNWIK = 0;
+
+ Scalar fMeanWIJ = (lJ*fractionalWI+lI*fractionalWJ)/l;
+ Scalar fMeanNWIJ = (lJ*fractionalNWI+lI*fractionalNWJ)/l;
+ Scalar fMeanWIK = (lJ*fractionalWI+lI*fractionalWK)/l;
+ Scalar fMeanNWIK = (lJ*fractionalNWI+lI*fractionalNWK)/l;
+
+ Scalar potentialWIJ = 0;
+ Scalar potentialNWIJ = 0;
+ Scalar potentialWIK = 0;
+ Scalar potentialNWIK = 0;
+
+ //if we are at the very first iteration we can't calculate phase potentials
+ if (!first)
+ {
+ // potentials from previous time step no available.
+ densityWIJ = rhoMeanWIJ;
+ densityNWIJ = rhoMeanNWIJ;
+ densityWIK = rhoMeanWIK;
+ densityNWIK = rhoMeanNWIK;
+
+ // Old pressure field
+ Scalar pressI=problem_.variables().pressure()[globalIdxI];
+ Scalar pressJ=problem_.variables().pressure()[globalIdxJ];
+ Scalar pressK=problem_.variables().pressure()[globalIdxK];
+ Scalar pressJK=(pressJ+pressK)/2;
+ Scalar pcJK=(pcJ+pcK)/2;
+
+ switch (pressureType)
+ {
+ case pw:
+ {
+ potentialWIJ = (pressI-pressJK)/l+
+ (densityWIJ-lJ/l*(kI+kK)/kI*(densityWIK-densityWIJ)/2)*ng;
+ potentialNWIJ = (pressI+pcI-(pressJK+pcJK))/l+
+ (densityNWIJ-lJ/l*(kI+kK)/kI*(densityNWIK-densityNWIJ)/2)*ng;
+ potentialWIK = (pressI-pressJK)/l+
+ (densityWIK-lJ/l*(kI+kK)/kI*(densityWIJ-densityWIK)/2)*ng;
+ potentialNWIK = (pressI+pcI-(pressJK+pcJK))/l+
+ (densityNWIK-lJ/l*(kI+kK)/kI*(densityNWIJ-densityNWIK)/2)*ng;
+ break;
+ }
+ case pn:
+ {
+ potentialWIJ = (pressI-pcI-(pressJK-pcJK))/l+
+ (densityWIJ-lJ/l*(kI+kK)/kI*(densityWIK-densityWIJ)/2)*ng;
+ potentialNWIJ = (pressI-pressJK)/l+
+ (densityNWIJ-lJ/l*(kI+kK)/kI*(densityNWIK-densityNWIJ)/2)*ng;
+ potentialWIK = (pressI-pcI-(pressJK-pcJK))/l+
+ (densityWIK-lJ/l*(kI+kK)/kI*(densityWIJ-densityWIK)/2)*ng;
+ potentialNWIK = (pressI-pressJK)/l+
+ (densityNWIK-lJ/l*(kI+kK)/kI*(densityNWIJ-densityNWIK)/2)*ng;
+ break;
+ }
+ case pglobal:
+ {
+ potentialWIJ = (pressI-fMeanNWIJ*pcI-(pressJK-fMeanNWIJ*pcJK))/l+
+ (densityWIJ-lJ/l*(kI+kK)/kI*(densityWIK-densityWIJ)/2)*ng;
+ potentialNWIJ = (pressI+fMeanWIJ*pcI-(pressJK+fMeanWIJ*pcJK))/l+
+ (densityNWIJ-lJ/l*(kI+kK)/kI*(densityNWIK-densityNWIJ)/2)*ng;
+ potentialWIK = (pressI-fMeanNWIK*pcI-(pressJK-fMeanNWIK*pcJK))/l+
+ (densityWIK-lJ/l*(kI+kK)/kI*(densityWIJ-densityWIK)/2)*ng;
+ potentialNWIK = (pressI+fMeanWIK*pcI-(pressJK+fMeanWIK*pcJK))/l+
+ (densityNWIK-lJ/l*(kI+kK)/kI*(densityNWIJ-densityNWIK)/2)*ng;
+ break;
+ }
+ }
+
+ //store potentials for further calculations (velocity, saturation, ...)
+
+ problem_.variables().potentialWetting(globalIdxI, isIndex) = potentialWIJ;
+ problem_.variables().potentialNonwetting(globalIdxI, isIndex) = potentialNWIJ;
+ problem_.variables().potentialWetting(globalIdxJ, isIndexJ) = -potentialWIJ;
+ problem_.variables().potentialNonwetting(globalIdxJ, isIndexJ) = -potentialNWIJ;
+ }
+
+ //do the upwinding of the mobility depending on the phase potentials
+ Scalar lambdaWIJ = (potentialWIJ > 0.) ? lambdaWI : lambdaWJ;
+ lambdaWIJ = (potentialWIJ == 0) ? 0.5 * (lambdaWI + lambdaWJ) : lambdaWIJ;
+ Scalar lambdaNWIJ = (potentialNWIJ > 0) ? lambdaNWI : lambdaNWJ;
+ lambdaNWIJ = (potentialNWIJ == 0) ? 0.5 * (lambdaNWI + lambdaNWJ) : lambdaNWIJ;
+ densityWIJ = (potentialWIJ > 0.) ? densityWI : densityWJ;
+ densityNWIJ = (potentialNWIJ > 0.) ? densityNWI : densityNWJ;
+ densityWIJ = (potentialWIJ == 0) ? rhoMeanWIJ : densityWIJ;
+ densityNWIJ = (potentialNWIJ == 0) ? rhoMeanNWIJ : densityNWIJ;
+
+ densityWIK = (potentialWIK > 0.) ? densityWI : densityWK;
+ densityNWIK = (potentialNWIK > 0.) ? densityNWI : densityNWK;
+ densityWIK = (potentialWIK == 0) ? rhoMeanWIK : densityWIK;
+ densityNWIK = (potentialNWIK == 0) ? rhoMeanNWIK : densityNWIK;
+
+
+ // update diagonal entry and right hand side
+
+ Scalar entry;
+ Scalar rightEntry;
+
+ entry = (lambdaWIJ + lambdaNWIJ)*kMean/l*faceArea;
+ rightEntry = faceArea*lambdaWIJ*kMean*ng*((densityWIJ)-(lJ/l)*(kI+kK)/kI*(densityWIK-densityWIJ)/2);
+ rightEntry += faceArea*lambdaNWIJ*kMean*ng*((densityNWIJ)-(lJ/l)*(kI+kK)/kI*(densityNWIK-densityNWIJ)/2);
+
+ switch (pressureType)
+ {
+ case pw:
+ {
+ rightEntry += faceArea*lambdaNWIJ*kMean*(pcJK-pcI)/l;
+ break;
+ }
+ case pn:
+ {
+ rightEntry -= faceArea*lambdaWIJ*kMean*(pcJK-pcI)/l;
+ break;
+ }
+ }
+
+ f_[globalIdxI] -= rightEntry;
+ f_[globalIdxJ] += rightEntry;
+
+ // set diagonal entry
+ A_[globalIdxI][globalIdxI] += entry;
+ A_[globalIdxI][globalIdxJ] -= 0.5*entry;
+ A_[globalIdxI][globalIdxK] -= 0.5*entry;
+
+ // set off-diagonal entry
+ A_[globalIdxJ][globalIdxI] -= entry;
+ A_[globalIdxJ][globalIdxJ] += 0.5*entry;
+ A_[globalIdxJ][globalIdxK] += 0.5*entry;
}
}
@@ -954,11 +954,11 @@ void FVPressure2Padaptive<TypeTag>::assemble(bool first)
if (!first)
{
- // Comment: In the non-adaptive case, one can upwind the density using the
- // old potential. Here, we use the mean density instead.
- // For the incompressible case, it does not make a difference, anyway.
- densityW = rhoMeanW;
- densityNW = rhoMeanNW;
+ // Comment: In the non-adaptive case, one can upwind the density using the
+ // old potential. Here, we use the mean density instead.
+ // For the incompressible case, it does not make a difference, anyway.
+ densityW = rhoMeanW;
+ densityNW = rhoMeanNW;
//calculate potential gradient
switch (pressureType)
diff --git a/dumux/decoupled/2p/diffusion/fv/fvvelocity2padaptive.hh b/dumux/decoupled/2p/diffusion/fv/fvvelocity2padaptive.hh
index 5a5e5daf63e2c814cfb8678896715d0c56fdf9ef..a6a007e515ab059cfb06725c2b3b0b8f0754fb57 100644
--- a/dumux/decoupled/2p/diffusion/fv/fvvelocity2padaptive.hh
+++ b/dumux/decoupled/2p/diffusion/fv/fvvelocity2padaptive.hh
@@ -116,7 +116,7 @@ public:
FVVelocity2Padaptive(Problem& problem)
: FVPressure2Padaptive<TypeTag>(problem)
{
- // todo: kompatibilität prüfen
+ // todo: kompatibilität prüfen
if (GET_PROP_VALUE(TypeTag, PTAG(EnableCompressibility)) && velocityType_ == vt)
{
DUNE_THROW(Dune::NotImplemented, "Total velocity - global pressure - model cannot be used with compressible fluids!");
@@ -282,7 +282,7 @@ void FVVelocity2Padaptive<TypeTag>::calculateVelocity()
!=isItEnd; ++isIt)
{
// local number of facet
-// isIndex++;
+// isIndex++;
int isIndex = isIt->indexInInside();
@@ -297,507 +297,507 @@ void FVVelocity2Padaptive<TypeTag>::calculateVelocity()
if (neighborPointer->level()==eIt.level())
{
- // neighbor cell center in global coordinates
- const GlobalPosition& globalPosNeighbor = neighborPointer->geometry().center();
-
- // distance vector between barycenters
- Dune::FieldVector<Scalar,dimWorld> distVec = globalPosNeighbor - globalPos;
-
- // compute distance between cell centers
- Scalar dist = distVec * unitOuterNormal;
-// Scalar dist = distVec.two_norm();
-
- // compute vectorized permeabilities
- FieldMatrix meanPermeability(0);
-
- this->problem().spatialParameters().meanK(meanPermeability,
- this->problem().spatialParameters().intrinsicPermeability(*eIt),
- this->problem().spatialParameters().intrinsicPermeability(*neighborPointer));
-
- Dune::FieldVector<Scalar, dim> permeability(0);
- meanPermeability.mv(unitOuterNormal, permeability);
-
- Scalar pressJ = this->problem().variables().pressure()[globalIdxJ];
- Scalar pcJ = this->problem().variables().capillaryPressure(globalIdxJ);
- Scalar lambdaWJ = this->problem().variables().mobilityWetting(globalIdxJ);
- Scalar lambdaNWJ = this->problem().variables().mobilityNonwetting(globalIdxJ);
- Scalar fractionalWJ = this->problem().variables().fracFlowFuncWetting(globalIdxJ);
- Scalar fractionalNWJ = this->problem().variables().fracFlowFuncNonwetting(globalIdxJ);
- Scalar densityWJ = this->problem().variables().densityWetting(globalIdxJ);
- Scalar densityNWJ = this->problem().variables().densityNonwetting(globalIdxJ);
-
- //calculate potential gradients
- Scalar potentialW = 0;
- Scalar potentialNW = 0;
-
- potentialW = this->problem().variables().potentialWetting(globalIdxI, isIndex);
- potentialNW = this->problem().variables().potentialNonwetting(globalIdxI, isIndex);
-
- Scalar densityW = (potentialW> 0.) ? densityWI : densityWJ;
- Scalar densityNW = (potentialNW> 0.) ? densityNWI : densityNWJ;
-
- densityW = (potentialW == 0.) ? 0.5 * (densityWI + densityWJ) : densityW;
- densityNW = (potentialNW == 0.) ? 0.5 * (densityNWI + densityNWJ) : densityNW;
-
- switch (this->pressureType)
- {
- case pw:
- {
- potentialW = (pressI - pressJ);
- potentialNW = (pressI - pressJ+ pcI - pcJ);
- break;
- }
- case pn:
- {
- potentialW = (pressI - pressJ - pcI + pcJ);
- potentialNW = (pressI - pressJ);
- break;
- }
- case pglobal:
- {
- potentialW = (pressI - pressJ - 0.5 * (fractionalNWI+fractionalNWJ)*(pcI - pcJ));
- potentialNW = (pressI - pressJ + 0.5 * (fractionalWI+fractionalWJ)*(pcI - pcJ));
- break;
- }
- }
-
- potentialW += densityW * (distVec * this->gravity);//delta z/delta x in unitOuterNormal[z]
- potentialNW += densityNW * (distVec * this->gravity);
-
- //store potentials for further calculations (velocity, saturation, ...)
- this->problem().variables().potentialWetting(globalIdxI, isIndex) = potentialW;
- this->problem().variables().potentialNonwetting(globalIdxI, isIndex) = potentialNW;
-
- //do the upwinding of the mobility depending on the phase potentials
- Scalar lambdaW = (potentialW > 0.) ? lambdaWI : lambdaWJ;
- lambdaW = (potentialW == 0) ? 0.5 * (lambdaWI + lambdaWJ) : lambdaW;
- Scalar lambdaNW = (potentialNW > 0.) ? lambdaNWI : lambdaNWJ;
- lambdaNW = (potentialNW == 0) ? 0.5 * (lambdaNWI + lambdaNWJ) : lambdaNW;
- densityW = (potentialW> 0.) ? densityWI : densityWJ;
- densityW = (potentialW == 0.) ? 0.5 * (densityWI + densityWJ) : densityW;
- densityNW = (potentialNW> 0.) ? densityNWI : densityNWJ;
- densityNW = (potentialNW == 0.) ? 0.5 * (densityNWI + densityNWJ) : densityNW;
-
- //calculate the gravity term
- Dune::FieldVector<Scalar,dimWorld> velocityW(permeability);
- Dune::FieldVector<Scalar,dimWorld> velocityNW(permeability);
- Dune::FieldVector<Scalar,dimWorld> gravityTermW(unitOuterNormal);
- Dune::FieldVector<Scalar,dimWorld> gravityTermNW(unitOuterNormal);
-
- gravityTermW *= (this->gravity*permeability)*(lambdaW * densityW);
- gravityTermNW *= (this->gravity*permeability)*(lambdaNW * densityNW);
-
- //calculate velocity depending on the pressure used -> use pc = pn - pw
- switch (this->pressureType)
- {
- case pw:
- {
- velocityW *= lambdaW * (pressI - pressJ)/dist;
- velocityNW *= lambdaNW * (pressI - pressJ)/dist + 0.5 * (lambdaNWI + lambdaNWJ) * (pcI - pcJ) / dist;
- velocityW += gravityTermW;
- velocityNW += gravityTermNW;
- break;
- }
- case pn:
- {
- velocityW *= lambdaW * (pressI - pressJ)/dist - 0.5 * (lambdaWI + lambdaWJ) * (pcI - pcJ) / dist;
- velocityNW *= lambdaNW * (pressI - pressJ) / dist;
- velocityW += gravityTermW;
- velocityNW += gravityTermNW;
- break;
- }
- case pglobal:
- {
- this->problem().variables().velocity()[globalIdxI][isIndex] = permeability;
- this->problem().variables().velocity()[globalIdxI][isIndex]*= (lambdaW + lambdaNW)* (pressI - pressJ )/ dist;
- this->problem().variables().velocity()[globalIdxI][isIndex] += gravityTermW;
- this->problem().variables().velocity()[globalIdxI][isIndex] += gravityTermNW;
- break;
- }
- }
-
- //store velocities
- switch (velocityType_)
- {
- case vw:
- {
- this->problem().variables().velocity()[globalIdxI][isIndex] = velocityW;
- this->problem().variables().velocitySecondPhase()[globalIdxI][isIndex] = velocityNW;
- break;
- }
- case vn:
- {
- this->problem().variables().velocity()[globalIdxI][isIndex] = velocityNW;
- this->problem().variables().velocitySecondPhase()[globalIdxI][isIndex] = velocityW;
- break;
- }
- case vt:
- {
- switch (this->pressureType)
- {
- case pw:
- {
- this->problem().variables().velocity()[globalIdxI][isIndex] = velocityW + velocityNW;
- this->problem().variables().velocitySecondPhase()[globalIdxI][isIndex] = velocityNW;
- break;
- }
- case pn:
- {
- this->problem().variables().velocity()[globalIdxI][isIndex] = velocityW + velocityNW;
- this->problem().variables().velocitySecondPhase()[globalIdxI][isIndex] = velocityW;
- break;
- }
- }
- break;
- }
- } //end of switch (velocityType_)
+ // neighbor cell center in global coordinates
+ const GlobalPosition& globalPosNeighbor = neighborPointer->geometry().center();
+
+ // distance vector between barycenters
+ Dune::FieldVector<Scalar,dimWorld> distVec = globalPosNeighbor - globalPos;
+
+ // compute distance between cell centers
+ Scalar dist = distVec * unitOuterNormal;
+// Scalar dist = distVec.two_norm();
+
+ // compute vectorized permeabilities
+ FieldMatrix meanPermeability(0);
+
+ this->problem().spatialParameters().meanK(meanPermeability,
+ this->problem().spatialParameters().intrinsicPermeability(*eIt),
+ this->problem().spatialParameters().intrinsicPermeability(*neighborPointer));
+
+ Dune::FieldVector<Scalar, dim> permeability(0);
+ meanPermeability.mv(unitOuterNormal, permeability);
+
+ Scalar pressJ = this->problem().variables().pressure()[globalIdxJ];
+ Scalar pcJ = this->problem().variables().capillaryPressure(globalIdxJ);
+ Scalar lambdaWJ = this->problem().variables().mobilityWetting(globalIdxJ);
+ Scalar lambdaNWJ = this->problem().variables().mobilityNonwetting(globalIdxJ);
+ Scalar fractionalWJ = this->problem().variables().fracFlowFuncWetting(globalIdxJ);
+ Scalar fractionalNWJ = this->problem().variables().fracFlowFuncNonwetting(globalIdxJ);
+ Scalar densityWJ = this->problem().variables().densityWetting(globalIdxJ);
+ Scalar densityNWJ = this->problem().variables().densityNonwetting(globalIdxJ);
+
+ //calculate potential gradients
+ Scalar potentialW = 0;
+ Scalar potentialNW = 0;
+
+ potentialW = this->problem().variables().potentialWetting(globalIdxI, isIndex);
+ potentialNW = this->problem().variables().potentialNonwetting(globalIdxI, isIndex);
+
+ Scalar densityW = (potentialW> 0.) ? densityWI : densityWJ;
+ Scalar densityNW = (potentialNW> 0.) ? densityNWI : densityNWJ;
+
+ densityW = (potentialW == 0.) ? 0.5 * (densityWI + densityWJ) : densityW;
+ densityNW = (potentialNW == 0.) ? 0.5 * (densityNWI + densityNWJ) : densityNW;
+
+ switch (this->pressureType)
+ {
+ case pw:
+ {
+ potentialW = (pressI - pressJ);
+ potentialNW = (pressI - pressJ+ pcI - pcJ);
+ break;
+ }
+ case pn:
+ {
+ potentialW = (pressI - pressJ - pcI + pcJ);
+ potentialNW = (pressI - pressJ);
+ break;
+ }
+ case pglobal:
+ {
+ potentialW = (pressI - pressJ - 0.5 * (fractionalNWI+fractionalNWJ)*(pcI - pcJ));
+ potentialNW = (pressI - pressJ + 0.5 * (fractionalWI+fractionalWJ)*(pcI - pcJ));
+ break;
+ }
+ }
+
+ potentialW += densityW * (distVec * this->gravity);//delta z/delta x in unitOuterNormal[z]
+ potentialNW += densityNW * (distVec * this->gravity);
+
+ //store potentials for further calculations (velocity, saturation, ...)
+ this->problem().variables().potentialWetting(globalIdxI, isIndex) = potentialW;
+ this->problem().variables().potentialNonwetting(globalIdxI, isIndex) = potentialNW;
+
+ //do the upwinding of the mobility depending on the phase potentials
+ Scalar lambdaW = (potentialW > 0.) ? lambdaWI : lambdaWJ;
+ lambdaW = (potentialW == 0) ? 0.5 * (lambdaWI + lambdaWJ) : lambdaW;
+ Scalar lambdaNW = (potentialNW > 0.) ? lambdaNWI : lambdaNWJ;
+ lambdaNW = (potentialNW == 0) ? 0.5 * (lambdaNWI + lambdaNWJ) : lambdaNW;
+ densityW = (potentialW> 0.) ? densityWI : densityWJ;
+ densityW = (potentialW == 0.) ? 0.5 * (densityWI + densityWJ) : densityW;
+ densityNW = (potentialNW> 0.) ? densityNWI : densityNWJ;
+ densityNW = (potentialNW == 0.) ? 0.5 * (densityNWI + densityNWJ) : densityNW;
+
+ //calculate the gravity term
+ Dune::FieldVector<Scalar,dimWorld> velocityW(permeability);
+ Dune::FieldVector<Scalar,dimWorld> velocityNW(permeability);
+ Dune::FieldVector<Scalar,dimWorld> gravityTermW(unitOuterNormal);
+ Dune::FieldVector<Scalar,dimWorld> gravityTermNW(unitOuterNormal);
+
+ gravityTermW *= (this->gravity*permeability)*(lambdaW * densityW);
+ gravityTermNW *= (this->gravity*permeability)*(lambdaNW * densityNW);
+
+ //calculate velocity depending on the pressure used -> use pc = pn - pw
+ switch (this->pressureType)
+ {
+ case pw:
+ {
+ velocityW *= lambdaW * (pressI - pressJ)/dist;
+ velocityNW *= lambdaNW * (pressI - pressJ)/dist + 0.5 * (lambdaNWI + lambdaNWJ) * (pcI - pcJ) / dist;
+ velocityW += gravityTermW;
+ velocityNW += gravityTermNW;
+ break;
+ }
+ case pn:
+ {
+ velocityW *= lambdaW * (pressI - pressJ)/dist - 0.5 * (lambdaWI + lambdaWJ) * (pcI - pcJ) / dist;
+ velocityNW *= lambdaNW * (pressI - pressJ) / dist;
+ velocityW += gravityTermW;
+ velocityNW += gravityTermNW;
+ break;
+ }
+ case pglobal:
+ {
+ this->problem().variables().velocity()[globalIdxI][isIndex] = permeability;
+ this->problem().variables().velocity()[globalIdxI][isIndex]*= (lambdaW + lambdaNW)* (pressI - pressJ )/ dist;
+ this->problem().variables().velocity()[globalIdxI][isIndex] += gravityTermW;
+ this->problem().variables().velocity()[globalIdxI][isIndex] += gravityTermNW;
+ break;
+ }
+ }
+
+ //store velocities
+ switch (velocityType_)
+ {
+ case vw:
+ {
+ this->problem().variables().velocity()[globalIdxI][isIndex] = velocityW;
+ this->problem().variables().velocitySecondPhase()[globalIdxI][isIndex] = velocityNW;
+ break;
+ }
+ case vn:
+ {
+ this->problem().variables().velocity()[globalIdxI][isIndex] = velocityNW;
+ this->problem().variables().velocitySecondPhase()[globalIdxI][isIndex] = velocityW;
+ break;
+ }
+ case vt:
+ {
+ switch (this->pressureType)
+ {
+ case pw:
+ {
+ this->problem().variables().velocity()[globalIdxI][isIndex] = velocityW + velocityNW;
+ this->problem().variables().velocitySecondPhase()[globalIdxI][isIndex] = velocityNW;
+ break;
+ }
+ case pn:
+ {
+ this->problem().variables().velocity()[globalIdxI][isIndex] = velocityW + velocityNW;
+ this->problem().variables().velocitySecondPhase()[globalIdxI][isIndex] = velocityW;
+ break;
+ }
+ }
+ break;
+ }
+ } //end of switch (velocityType_)
} //End of case "same level"
if (neighborPointer->level()==eIt.level()+1)
{
- int globalIdxK = 0;
- ElementPointer thirdCellPointer = isIt->outside();
- bool foundK=false;
- bool foundIJ=false;
- // We are looking for two things:
- // IsIndexJ, the index of the interface from the neighbor-cell point of view
- // GlobalIdxK, the index of the third cell
- // for efficienty this is done in one IntersectionIterator-Loop
-
- Scalar areaWeight = 0;
-
- // Intersectioniterator around cell J
-// int isIndexJ = -1;
- IntersectionIterator isItEndJ = this->problem().gridView().iend(*neighborPointer);
-
- for (IntersectionIterator isItJ = this->problem().gridView().ibegin(*neighborPointer); isItJ != isItEndJ; ++isItJ)
- {
- // increase isIndexJJ, if it is not found yet
- if (!foundIJ)
-// isIndexJ++;
- if (isItJ->neighbor())
- {
- ElementPointer neighborPointer2 = isItJ->outside();
-
- // Neighbor of neighbor is Cell I?
- if (this->problem().variables().index(*neighborPointer2)==globalIdxI)
- {
- foundIJ=true;
- areaWeight = isItJ->geometry().volume();
- }
- else
- {
- if (neighborPointer2->level()==eIt.level()+1)
- {
- // To verify, we found the correct Cell K, we check
- // - is level(K)=level(J)?
- // - is distance(IJ)=distance(IK)?
- // - K is neighbor of J.
- const GlobalPosition& globalPosNeighbor = neighborPointer->geometry().center();
- const GlobalPosition& globalPosThirdCell = neighborPointer2->geometry().center();
- Dune::FieldVector<Scalar, dimWorld> distVecIJ = globalPosNeighbor - globalPos;
- Dune::FieldVector<Scalar, dimWorld> distVecIK = globalPosThirdCell - globalPos;
- if (((distVecIK.two_norm()-distVecIJ.two_norm()))/distVecIJ.two_norm() < 0.001)
- {
- globalIdxK= this->problem().variables().index(*neighborPointer2);
- thirdCellPointer = neighborPointer2;
- foundK=true;
- }
-
- }
- }
- }
- if (foundIJ && foundK) break;
- }
-
- int isIndexJ = isIt->indexInOutside();
- areaWeight /= isIt->geometry().volume();
-// areaWeight = 1.0;
-
- // neighbor cell center in global coordinates
- const GlobalPosition& globalPosNeighbor = neighborPointer->geometry().center();
- const GlobalPosition& globalPosInterface = isIt->geometry().center();
-
- Dune::FieldVector<Scalar,dimWorld> distVec = globalPosInterface - globalPos;
- Scalar lI= distVec*unitOuterNormal;
- distVec = globalPosNeighbor - globalPosInterface;
- Scalar lJ= distVec*unitOuterNormal;
- Scalar l=lI+lJ;
-
- FieldMatrix permeabilityI(0);
- FieldMatrix permeabilityJ(0);
- FieldMatrix permeabilityK(0);
-
- this->problem().spatialParameters().meanK(permeabilityI,
- this->problem().spatialParameters().intrinsicPermeability(*eIt));
- this->problem().spatialParameters().meanK(permeabilityJ,
- this->problem().spatialParameters().intrinsicPermeability(*neighborPointer));
- this->problem().spatialParameters().meanK(permeabilityK,
- this->problem().spatialParameters().intrinsicPermeability(*thirdCellPointer));
-
- // Calculate permeablity component normal to interface
- Scalar kI, kJ, kK, ng, kMean;//, gI, gJ, gK;
- Dune::FieldVector<Scalar, dim> permI(0);
- Dune::FieldVector<Scalar, dim> permJ(0);
- Dune::FieldVector<Scalar, dim> permK(0);
-
- permeabilityI.mv(unitOuterNormal, permI);
- permeabilityJ.mv(unitOuterNormal, permJ);
- permeabilityK.mv(unitOuterNormal, permK);
-
- // kI,kJ,kK = (n^T)Kn
- kI=unitOuterNormal*permI;
- kJ=unitOuterNormal*permJ;
- kK=unitOuterNormal*permK;
- kMean=kI*kJ*kK*l/(kJ*kK*lI+kI*(kJ+kK)/2*lJ);
-
- ng=this->gravity*unitOuterNormal;
-
- // find secondary variables in cell J and K
- Scalar pressJ = this->problem().variables().pressure()[globalIdxJ];
- Scalar pcJ = this->problem().variables().capillaryPressure(globalIdxJ);
- Scalar lambdaWJ = this->problem().variables().mobilityWetting(globalIdxJ);
- Scalar lambdaNWJ = this->problem().variables().mobilityNonwetting(globalIdxJ);
- Scalar densityWJ = this->problem().variables().densityWetting(globalIdxJ);
- Scalar densityNWJ = this->problem().variables().densityNonwetting(globalIdxJ);
- Scalar fractionalWJ = this->problem().variables().fracFlowFuncWetting(globalIdxJ);
- Scalar fractionalNWJ = this->problem().variables().fracFlowFuncNonwetting(globalIdxJ);
-
- Scalar pressK = this->problem().variables().pressure()[globalIdxK];
- Scalar pcK = this->problem().variables().capillaryPressure(globalIdxK);
- Scalar densityWK = this->problem().variables().densityWetting(globalIdxK);
- Scalar densityNWK = this->problem().variables().densityNonwetting(globalIdxK);
- Scalar fractionalWK = this->problem().variables().fracFlowFuncWetting(globalIdxK);
- Scalar fractionalNWK = this->problem().variables().fracFlowFuncNonwetting(globalIdxK);
-
- Scalar pressJK=(pressJ+pressK)/2;
- Scalar pcJK=(pcJ+pcK)/2;
-
-
- // calculate potential gradients
- // reuse potentials from fvpressure2padaptive
-
- Scalar potentialWIJ = this->problem().variables().potentialWetting(globalIdxI, isIndex);
- Scalar potentialNWIJ = this->problem().variables().potentialNonwetting(globalIdxI, isIndex);
- // We dont know the insideIndex of interface IK
- Scalar potentialWIK = potentialWIJ;
- Scalar potentialNWIK = potentialNWIJ;
- // preliminary potential. The "real" ones are found below
-
- // Comment: reversed weighting is plausible, too (swap lJ and lI)
- Scalar rhoMeanWIJ = (lJ*densityWI + lI*densityWJ)/l;
- Scalar rhoMeanNWIJ = (lJ*densityNWI + lI*densityNWJ)/l;
- Scalar rhoMeanWIK = (lJ*densityWI + lI*densityWK)/l;
- Scalar rhoMeanNWIK = (lJ*densityNWI + lI*densityNWK)/l;
-
- Scalar fMeanWIJ = (lJ*fractionalWI+lI*fractionalWJ)/l;
- Scalar fMeanNWIJ = (lJ*fractionalNWI+lI*fractionalNWJ)/l;
- Scalar fMeanWIK = (lJ*fractionalWI+lI*fractionalWK)/l;
- Scalar fMeanNWIK = (lJ*fractionalNWI+lI*fractionalNWK)/l;
-
- // Upwinding for finding the upwinding direction
- Scalar densityWIJ = (potentialWIJ> 0.) ? densityWI : densityWJ;
- Scalar densityNWIJ = (potentialNWIJ> 0.) ? densityNWI : densityNWJ;
- Scalar densityWIK = (potentialWIK> 0.) ? densityWI : densityWK;
- Scalar densityNWIK = (potentialNWIK> 0.) ? densityNWI : densityNWK;
-
- densityWIJ = (potentialWIJ == 0.) ? rhoMeanWIJ : densityWIJ;
- densityNWIJ = (potentialNWIJ == 0.) ? rhoMeanNWIJ : densityNWIJ;
- densityWIK = (potentialWIK == 0.) ? rhoMeanWIK : densityWIK;
- densityNWIK = (potentialNWIK == 0.) ? rhoMeanNWIK : densityNWIK;
-
- Scalar fractionalWIJ = (potentialWIJ> 0.) ? fractionalWI : fractionalWJ;
- Scalar fractionalNWIJ = (potentialNWIJ> 0.) ? fractionalNWI : fractionalNWJ;
- Scalar fractionalWIK = (potentialWIK> 0.) ? fractionalWI : fractionalWK;
- Scalar fractionalNWIK = (potentialNWIK> 0.) ? fractionalNWI : fractionalNWK;
-
- fractionalWIJ = (potentialWIJ == 0.) ? fMeanWIJ : fractionalWIJ;
- fractionalNWIJ = (potentialNWIJ == 0.) ? fMeanNWIJ : fractionalNWIJ;
- fractionalWIK = (potentialWIK == 0.) ? fMeanWIK : fractionalWIK;
- fractionalNWIK = (potentialNWIK == 0.) ? fMeanNWIK : fractionalNWIK;
-
- switch (this->pressureType)
- {
- case pw:
- {
- potentialWIJ = (pressI-pressJK)/l+
- (densityWIJ-lJ/l*(kI+kK)/kI*(densityWIK-densityWIJ)/2)*ng;
- potentialNWIJ = (pressI+pcI-(pressJK+pcJK))/l+
- (densityNWIJ-lJ/l*(kI+kK)/kI*(densityNWIK-densityNWIJ)/2)*ng;
- potentialWIK = (pressI-pressJK)/l+
- (densityWIK-lJ/l*(kI+kK)/kI*(densityWIJ-densityWIK)/2)*ng;
- potentialNWIK = (pressI+pcI-(pressJK+pcJK))/l+
- (densityNWIK-lJ/l*(kI+kK)/kI*(densityNWIJ-densityNWIK)/2)*ng;
- break;
- }
- case pn:
- {
- potentialWIJ = (pressI-pcI-(pressJK-pcJK))/l+
- (densityWIJ-lJ/l*(kI+kK)/kI*(densityWIK-densityWIJ)/2)*ng;
- potentialNWIJ = (pressI-pressJK)/l+
- (densityNWIJ-lJ/l*(kI+kK)/kI*(densityNWIK-densityNWIJ)/2)*ng;
- potentialWIK = (pressI-pcI-(pressJK-pcJK))/l+
- (densityWIK-lJ/l*(kI+kK)/kI*(densityWIJ-densityWIK)/2)*ng;
- potentialNWIK = (pressI-pressJK)/l+
- (densityNWIK-lJ/l*(kI+kK)/kI*(densityNWIJ-densityNWIK)/2)*ng;
- break;
- }
- case pglobal:
- {
- potentialWIJ = (pressI-fractionalNWIJ*pcI-(pressJK-fractionalNWIJ*pcJK))/l+
- (densityWIJ-lJ/l*(kI+kK)/kI*(densityWIK-densityWIJ)/2)*ng;
- potentialNWIJ = (pressI+fractionalWIJ*pcI-(pressJK+fractionalWIJ*pcJK))/l+
- (densityNWIJ-lJ/l*(kI+kK)/kI*(densityNWIK-densityNWIJ)/2)*ng;
- potentialWIK = (pressI-fractionalNWIK*pcI-(pressJK-fractionalNWIK*pcJK))/l+
- (densityWIK-lJ/l*(kI+kK)/kI*(densityWIJ-densityWIK)/2)*ng;
- potentialNWIK = (pressI+fractionalWIK*pcI-(pressJK+fractionalWIK*pcJK))/l+
- (densityNWIK-lJ/l*(kI+kK)/kI*(densityNWIJ-densityNWIK)/2)*ng;
- break;
- }
- }
-
- //store potentials for further calculations (velocity, saturation, ...)
- // these quantities only have correct sign (needed for upwinding)
- // potentials are defined slightly different for adaptive scheme
- this->problem().variables().potentialWetting(globalIdxI, isIndex) = potentialWIJ;
- this->problem().variables().potentialNonwetting(globalIdxI, isIndex) = potentialNWIJ;
- this->problem().variables().potentialWetting(globalIdxJ, isIndexJ) = -potentialWIJ;
- this->problem().variables().potentialNonwetting(globalIdxJ, isIndexJ) = -potentialNWIJ;
-
-
- //do the upwinding of the mobility depending on the phase potentials
- Scalar lambdaWIJ = (potentialWIJ > 0.) ? lambdaWI : lambdaWJ;
- lambdaWIJ = (potentialWIJ == 0) ? 0.5 * (lambdaWI + lambdaWJ) : lambdaWIJ;
- Scalar lambdaNWIJ = (potentialNWIJ > 0.) ? lambdaNWI : lambdaNWJ;
- lambdaNWIJ = (potentialNWIJ == 0) ? 0.5 * (lambdaNWI + lambdaNWJ) : lambdaNWIJ;
-
- densityWIJ = (potentialWIJ> 0.) ? densityWI : densityWJ;
- densityNWIJ = (potentialNWIJ> 0.) ? densityNWI : densityNWJ;
- densityWIK = (potentialWIK> 0.) ? densityWI : densityWK;
- densityNWIK = (potentialNWIK> 0.) ? densityNWI : densityNWK;
-
- densityWIJ = (potentialWIJ == 0.) ? rhoMeanWIJ : densityWIJ;
- densityNWIJ = (potentialNWIJ == 0.) ? rhoMeanNWIJ : densityNWIJ;
- densityWIK = (potentialWIK == 0.) ? rhoMeanWIK : densityWIK;
- densityNWIK = (potentialNWIK == 0.) ? rhoMeanNWIK : densityNWIK;
-
- fractionalWIJ = (potentialWIJ> 0.) ? fractionalWI : fractionalWJ;
- fractionalNWIJ = (potentialNWIJ> 0.) ? fractionalNWI : fractionalNWJ;
- fractionalWIK = (potentialWIK> 0.) ? fractionalWI : fractionalWK;
- fractionalNWIK = (potentialNWIK> 0.) ? fractionalNWI : fractionalNWK;
-
- fractionalWIJ = (potentialWIJ == 0.) ? fMeanWIJ : fractionalWIJ;
- fractionalNWIJ = (potentialNWIJ == 0.) ? fMeanNWIJ : fractionalNWIJ;
- fractionalWIK = (potentialWIK == 0.) ? fMeanWIK : fractionalWIK;
- fractionalNWIK = (potentialNWIK == 0.) ? fMeanNWIK : fractionalNWIK;
-
- //calculate velocities and the gravity term
- Dune::FieldVector<Scalar,dimWorld> velocityW(unitOuterNormal);
- Dune::FieldVector<Scalar,dimWorld> velocityNW(unitOuterNormal);
- Dune::FieldVector<Scalar,dimWorld> gravityTermW(unitOuterNormal);
- Dune::FieldVector<Scalar,dimWorld> gravityTermNW(unitOuterNormal);
-
- gravityTermW *= lambdaWIJ*kMean*ng;
- gravityTermW *= densityWIJ-(lJ/l)*(kI+kK)/kI*(densityWIK-densityWIJ)/2;
- gravityTermNW *= lambdaNWIJ*kMean*ng;
- gravityTermNW *= densityNWIJ-(lJ/l)*(kI+kK)/kI*(densityNWIK-densityNWIJ)/2;
-
- switch (this->pressureType)
- {
- case pw:
- {
- velocityW *= lambdaWIJ*kMean*(pressI-pressJK)/l;
- velocityNW *= lambdaNWIJ*kMean*(pressI+pcI-(pressJK+pcJK))/l;
-
- velocityW += gravityTermW;
- velocityNW += gravityTermNW;
- break;
- }
- case pn:
- {
- velocityNW *= lambdaNWIJ*kMean*(pressI-pressJK)/l;
- velocityW *= lambdaWIJ*kMean*(pressI-pcI-(pressJK-pcJK))/l;
-
- velocityW += gravityTermW;
- velocityNW += gravityTermNW;
- break;
- }
- case pglobal:
- {
- velocityW *= lambdaWIJ*kMean*(pressI-fractionalNWIJ*pcI-(pressJK-fractionalNWIJ*pcJK))/l;
- velocityNW *= lambdaNWIJ*kMean*(pressI+fractionalWIJ*pcI-(pressJK+fractionalWIJ*pcJK))/l;
-
- velocityW += gravityTermW;
- velocityNW += gravityTermNW;
- break;
- }
- }
-
- switch (velocityType_)
- {
- case vw:
- {
- Dune::FieldVector<Scalar,dimWorld> vel(velocityW);
- vel*=areaWeight;
- this->problem().variables().velocity()[globalIdxI][isIndex] += vel;
- vel = velocityNW;
- vel*=areaWeight;
- this->problem().variables().velocitySecondPhase()[globalIdxI][isIndex] += vel;
+ int globalIdxK = 0;
+ ElementPointer thirdCellPointer = isIt->outside();
+ bool foundK=false;
+ bool foundIJ=false;
+ // We are looking for two things:
+ // IsIndexJ, the index of the interface from the neighbor-cell point of view
+ // GlobalIdxK, the index of the third cell
+ // for efficienty this is done in one IntersectionIterator-Loop
+
+ Scalar areaWeight = 0;
+
+ // Intersectioniterator around cell J
+// int isIndexJ = -1;
+ IntersectionIterator isItEndJ = this->problem().gridView().iend(*neighborPointer);
+
+ for (IntersectionIterator isItJ = this->problem().gridView().ibegin(*neighborPointer); isItJ != isItEndJ; ++isItJ)
+ {
+ // increase isIndexJJ, if it is not found yet
+ if (!foundIJ)
+// isIndexJ++;
+ if (isItJ->neighbor())
+ {
+ ElementPointer neighborPointer2 = isItJ->outside();
+
+ // Neighbor of neighbor is Cell I?
+ if (this->problem().variables().index(*neighborPointer2)==globalIdxI)
+ {
+ foundIJ=true;
+ areaWeight = isItJ->geometry().volume();
+ }
+ else
+ {
+ if (neighborPointer2->level()==eIt.level()+1)
+ {
+ // To verify, we found the correct Cell K, we check
+ // - is level(K)=level(J)?
+ // - is distance(IJ)=distance(IK)?
+ // - K is neighbor of J.
+ const GlobalPosition& globalPosNeighbor = neighborPointer->geometry().center();
+ const GlobalPosition& globalPosThirdCell = neighborPointer2->geometry().center();
+ Dune::FieldVector<Scalar, dimWorld> distVecIJ = globalPosNeighbor - globalPos;
+ Dune::FieldVector<Scalar, dimWorld> distVecIK = globalPosThirdCell - globalPos;
+ if (((distVecIK.two_norm()-distVecIJ.two_norm()))/distVecIJ.two_norm() < 0.001)
+ {
+ globalIdxK= this->problem().variables().index(*neighborPointer2);
+ thirdCellPointer = neighborPointer2;
+ foundK=true;
+ }
+
+ }
+ }
+ }
+ if (foundIJ && foundK) break;
+ }
+
+ int isIndexJ = isIt->indexInOutside();
+ areaWeight /= isIt->geometry().volume();
+// areaWeight = 1.0;
+
+ // neighbor cell center in global coordinates
+ const GlobalPosition& globalPosNeighbor = neighborPointer->geometry().center();
+ const GlobalPosition& globalPosInterface = isIt->geometry().center();
+
+ Dune::FieldVector<Scalar,dimWorld> distVec = globalPosInterface - globalPos;
+ Scalar lI= distVec*unitOuterNormal;
+ distVec = globalPosNeighbor - globalPosInterface;
+ Scalar lJ= distVec*unitOuterNormal;
+ Scalar l=lI+lJ;
+
+ FieldMatrix permeabilityI(0);
+ FieldMatrix permeabilityJ(0);
+ FieldMatrix permeabilityK(0);
+
+ this->problem().spatialParameters().meanK(permeabilityI,
+ this->problem().spatialParameters().intrinsicPermeability(*eIt));
+ this->problem().spatialParameters().meanK(permeabilityJ,
+ this->problem().spatialParameters().intrinsicPermeability(*neighborPointer));
+ this->problem().spatialParameters().meanK(permeabilityK,
+ this->problem().spatialParameters().intrinsicPermeability(*thirdCellPointer));
+
+ // Calculate permeablity component normal to interface
+ Scalar kI, kJ, kK, ng, kMean;//, gI, gJ, gK;
+ Dune::FieldVector<Scalar, dim> permI(0);
+ Dune::FieldVector<Scalar, dim> permJ(0);
+ Dune::FieldVector<Scalar, dim> permK(0);
+
+ permeabilityI.mv(unitOuterNormal, permI);
+ permeabilityJ.mv(unitOuterNormal, permJ);
+ permeabilityK.mv(unitOuterNormal, permK);
+
+ // kI,kJ,kK = (n^T)Kn
+ kI=unitOuterNormal*permI;
+ kJ=unitOuterNormal*permJ;
+ kK=unitOuterNormal*permK;
+ kMean=kI*kJ*kK*l/(kJ*kK*lI+kI*(kJ+kK)/2*lJ);
+
+ ng=this->gravity*unitOuterNormal;
+
+ // find secondary variables in cell J and K
+ Scalar pressJ = this->problem().variables().pressure()[globalIdxJ];
+ Scalar pcJ = this->problem().variables().capillaryPressure(globalIdxJ);
+ Scalar lambdaWJ = this->problem().variables().mobilityWetting(globalIdxJ);
+ Scalar lambdaNWJ = this->problem().variables().mobilityNonwetting(globalIdxJ);
+ Scalar densityWJ = this->problem().variables().densityWetting(globalIdxJ);
+ Scalar densityNWJ = this->problem().variables().densityNonwetting(globalIdxJ);
+ Scalar fractionalWJ = this->problem().variables().fracFlowFuncWetting(globalIdxJ);
+ Scalar fractionalNWJ = this->problem().variables().fracFlowFuncNonwetting(globalIdxJ);
+
+ Scalar pressK = this->problem().variables().pressure()[globalIdxK];
+ Scalar pcK = this->problem().variables().capillaryPressure(globalIdxK);
+ Scalar densityWK = this->problem().variables().densityWetting(globalIdxK);
+ Scalar densityNWK = this->problem().variables().densityNonwetting(globalIdxK);
+ Scalar fractionalWK = this->problem().variables().fracFlowFuncWetting(globalIdxK);
+ Scalar fractionalNWK = this->problem().variables().fracFlowFuncNonwetting(globalIdxK);
+
+ Scalar pressJK=(pressJ+pressK)/2;
+ Scalar pcJK=(pcJ+pcK)/2;
+
+
+ // calculate potential gradients
+ // reuse potentials from fvpressure2padaptive
+
+ Scalar potentialWIJ = this->problem().variables().potentialWetting(globalIdxI, isIndex);
+ Scalar potentialNWIJ = this->problem().variables().potentialNonwetting(globalIdxI, isIndex);
+ // We dont know the insideIndex of interface IK
+ Scalar potentialWIK = potentialWIJ;
+ Scalar potentialNWIK = potentialNWIJ;
+ // preliminary potential. The "real" ones are found below
+
+ // Comment: reversed weighting is plausible, too (swap lJ and lI)
+ Scalar rhoMeanWIJ = (lJ*densityWI + lI*densityWJ)/l;
+ Scalar rhoMeanNWIJ = (lJ*densityNWI + lI*densityNWJ)/l;
+ Scalar rhoMeanWIK = (lJ*densityWI + lI*densityWK)/l;
+ Scalar rhoMeanNWIK = (lJ*densityNWI + lI*densityNWK)/l;
+
+ Scalar fMeanWIJ = (lJ*fractionalWI+lI*fractionalWJ)/l;
+ Scalar fMeanNWIJ = (lJ*fractionalNWI+lI*fractionalNWJ)/l;
+ Scalar fMeanWIK = (lJ*fractionalWI+lI*fractionalWK)/l;
+ Scalar fMeanNWIK = (lJ*fractionalNWI+lI*fractionalNWK)/l;
+
+ // Upwinding for finding the upwinding direction
+ Scalar densityWIJ = (potentialWIJ> 0.) ? densityWI : densityWJ;
+ Scalar densityNWIJ = (potentialNWIJ> 0.) ? densityNWI : densityNWJ;
+ Scalar densityWIK = (potentialWIK> 0.) ? densityWI : densityWK;
+ Scalar densityNWIK = (potentialNWIK> 0.) ? densityNWI : densityNWK;
+
+ densityWIJ = (potentialWIJ == 0.) ? rhoMeanWIJ : densityWIJ;
+ densityNWIJ = (potentialNWIJ == 0.) ? rhoMeanNWIJ : densityNWIJ;
+ densityWIK = (potentialWIK == 0.) ? rhoMeanWIK : densityWIK;
+ densityNWIK = (potentialNWIK == 0.) ? rhoMeanNWIK : densityNWIK;
+
+ Scalar fractionalWIJ = (potentialWIJ> 0.) ? fractionalWI : fractionalWJ;
+ Scalar fractionalNWIJ = (potentialNWIJ> 0.) ? fractionalNWI : fractionalNWJ;
+ Scalar fractionalWIK = (potentialWIK> 0.) ? fractionalWI : fractionalWK;
+ Scalar fractionalNWIK = (potentialNWIK> 0.) ? fractionalNWI : fractionalNWK;
+
+ fractionalWIJ = (potentialWIJ == 0.) ? fMeanWIJ : fractionalWIJ;
+ fractionalNWIJ = (potentialNWIJ == 0.) ? fMeanNWIJ : fractionalNWIJ;
+ fractionalWIK = (potentialWIK == 0.) ? fMeanWIK : fractionalWIK;
+ fractionalNWIK = (potentialNWIK == 0.) ? fMeanNWIK : fractionalNWIK;
+
+ switch (this->pressureType)
+ {
+ case pw:
+ {
+ potentialWIJ = (pressI-pressJK)/l+
+ (densityWIJ-lJ/l*(kI+kK)/kI*(densityWIK-densityWIJ)/2)*ng;
+ potentialNWIJ = (pressI+pcI-(pressJK+pcJK))/l+
+ (densityNWIJ-lJ/l*(kI+kK)/kI*(densityNWIK-densityNWIJ)/2)*ng;
+ potentialWIK = (pressI-pressJK)/l+
+ (densityWIK-lJ/l*(kI+kK)/kI*(densityWIJ-densityWIK)/2)*ng;
+ potentialNWIK = (pressI+pcI-(pressJK+pcJK))/l+
+ (densityNWIK-lJ/l*(kI+kK)/kI*(densityNWIJ-densityNWIK)/2)*ng;
+ break;
+ }
+ case pn:
+ {
+ potentialWIJ = (pressI-pcI-(pressJK-pcJK))/l+
+ (densityWIJ-lJ/l*(kI+kK)/kI*(densityWIK-densityWIJ)/2)*ng;
+ potentialNWIJ = (pressI-pressJK)/l+
+ (densityNWIJ-lJ/l*(kI+kK)/kI*(densityNWIK-densityNWIJ)/2)*ng;
+ potentialWIK = (pressI-pcI-(pressJK-pcJK))/l+
+ (densityWIK-lJ/l*(kI+kK)/kI*(densityWIJ-densityWIK)/2)*ng;
+ potentialNWIK = (pressI-pressJK)/l+
+ (densityNWIK-lJ/l*(kI+kK)/kI*(densityNWIJ-densityNWIK)/2)*ng;
+ break;
+ }
+ case pglobal:
+ {
+ potentialWIJ = (pressI-fractionalNWIJ*pcI-(pressJK-fractionalNWIJ*pcJK))/l+
+ (densityWIJ-lJ/l*(kI+kK)/kI*(densityWIK-densityWIJ)/2)*ng;
+ potentialNWIJ = (pressI+fractionalWIJ*pcI-(pressJK+fractionalWIJ*pcJK))/l+
+ (densityNWIJ-lJ/l*(kI+kK)/kI*(densityNWIK-densityNWIJ)/2)*ng;
+ potentialWIK = (pressI-fractionalNWIK*pcI-(pressJK-fractionalNWIK*pcJK))/l+
+ (densityWIK-lJ/l*(kI+kK)/kI*(densityWIJ-densityWIK)/2)*ng;
+ potentialNWIK = (pressI+fractionalWIK*pcI-(pressJK+fractionalWIK*pcJK))/l+
+ (densityNWIK-lJ/l*(kI+kK)/kI*(densityNWIJ-densityNWIK)/2)*ng;
+ break;
+ }
+ }
+
+ //store potentials for further calculations (velocity, saturation, ...)
+ // these quantities only have correct sign (needed for upwinding)
+ // potentials are defined slightly different for adaptive scheme
+ this->problem().variables().potentialWetting(globalIdxI, isIndex) = potentialWIJ;
+ this->problem().variables().potentialNonwetting(globalIdxI, isIndex) = potentialNWIJ;
+ this->problem().variables().potentialWetting(globalIdxJ, isIndexJ) = -potentialWIJ;
+ this->problem().variables().potentialNonwetting(globalIdxJ, isIndexJ) = -potentialNWIJ;
+
+
+ //do the upwinding of the mobility depending on the phase potentials
+ Scalar lambdaWIJ = (potentialWIJ > 0.) ? lambdaWI : lambdaWJ;
+ lambdaWIJ = (potentialWIJ == 0) ? 0.5 * (lambdaWI + lambdaWJ) : lambdaWIJ;
+ Scalar lambdaNWIJ = (potentialNWIJ > 0.) ? lambdaNWI : lambdaNWJ;
+ lambdaNWIJ = (potentialNWIJ == 0) ? 0.5 * (lambdaNWI + lambdaNWJ) : lambdaNWIJ;
+
+ densityWIJ = (potentialWIJ> 0.) ? densityWI : densityWJ;
+ densityNWIJ = (potentialNWIJ> 0.) ? densityNWI : densityNWJ;
+ densityWIK = (potentialWIK> 0.) ? densityWI : densityWK;
+ densityNWIK = (potentialNWIK> 0.) ? densityNWI : densityNWK;
+
+ densityWIJ = (potentialWIJ == 0.) ? rhoMeanWIJ : densityWIJ;
+ densityNWIJ = (potentialNWIJ == 0.) ? rhoMeanNWIJ : densityNWIJ;
+ densityWIK = (potentialWIK == 0.) ? rhoMeanWIK : densityWIK;
+ densityNWIK = (potentialNWIK == 0.) ? rhoMeanNWIK : densityNWIK;
+
+ fractionalWIJ = (potentialWIJ> 0.) ? fractionalWI : fractionalWJ;
+ fractionalNWIJ = (potentialNWIJ> 0.) ? fractionalNWI : fractionalNWJ;
+ fractionalWIK = (potentialWIK> 0.) ? fractionalWI : fractionalWK;
+ fractionalNWIK = (potentialNWIK> 0.) ? fractionalNWI : fractionalNWK;
+
+ fractionalWIJ = (potentialWIJ == 0.) ? fMeanWIJ : fractionalWIJ;
+ fractionalNWIJ = (potentialNWIJ == 0.) ? fMeanNWIJ : fractionalNWIJ;
+ fractionalWIK = (potentialWIK == 0.) ? fMeanWIK : fractionalWIK;
+ fractionalNWIK = (potentialNWIK == 0.) ? fMeanNWIK : fractionalNWIK;
+
+ //calculate velocities and the gravity term
+ Dune::FieldVector<Scalar,dimWorld> velocityW(unitOuterNormal);
+ Dune::FieldVector<Scalar,dimWorld> velocityNW(unitOuterNormal);
+ Dune::FieldVector<Scalar,dimWorld> gravityTermW(unitOuterNormal);
+ Dune::FieldVector<Scalar,dimWorld> gravityTermNW(unitOuterNormal);
+
+ gravityTermW *= lambdaWIJ*kMean*ng;
+ gravityTermW *= densityWIJ-(lJ/l)*(kI+kK)/kI*(densityWIK-densityWIJ)/2;
+ gravityTermNW *= lambdaNWIJ*kMean*ng;
+ gravityTermNW *= densityNWIJ-(lJ/l)*(kI+kK)/kI*(densityNWIK-densityNWIJ)/2;
+
+ switch (this->pressureType)
+ {
+ case pw:
+ {
+ velocityW *= lambdaWIJ*kMean*(pressI-pressJK)/l;
+ velocityNW *= lambdaNWIJ*kMean*(pressI+pcI-(pressJK+pcJK))/l;
+
+ velocityW += gravityTermW;
+ velocityNW += gravityTermNW;
+ break;
+ }
+ case pn:
+ {
+ velocityNW *= lambdaNWIJ*kMean*(pressI-pressJK)/l;
+ velocityW *= lambdaWIJ*kMean*(pressI-pcI-(pressJK-pcJK))/l;
+
+ velocityW += gravityTermW;
+ velocityNW += gravityTermNW;
+ break;
+ }
+ case pglobal:
+ {
+ velocityW *= lambdaWIJ*kMean*(pressI-fractionalNWIJ*pcI-(pressJK-fractionalNWIJ*pcJK))/l;
+ velocityNW *= lambdaNWIJ*kMean*(pressI+fractionalWIJ*pcI-(pressJK+fractionalWIJ*pcJK))/l;
+
+ velocityW += gravityTermW;
+ velocityNW += gravityTermNW;
+ break;
+ }
+ }
+
+ switch (velocityType_)
+ {
+ case vw:
+ {
+ Dune::FieldVector<Scalar,dimWorld> vel(velocityW);
+ vel*=areaWeight;
+ this->problem().variables().velocity()[globalIdxI][isIndex] += vel;
+ vel = velocityNW;
+ vel*=areaWeight;
+ this->problem().variables().velocitySecondPhase()[globalIdxI][isIndex] += vel;
// this->problem().variables().velocity()[globalIdxI][isIndex] = velocityW;
// this->problem().variables().velocitySecondPhase()[globalIdxI][isIndex] = velocityNW;
- this->problem().variables().velocity()[globalIdxJ][isIndexJ] = velocityW;
- this->problem().variables().velocitySecondPhase()[globalIdxJ][isIndexJ] = velocityNW;
- break;
- }
- case vn:
- {
+ this->problem().variables().velocity()[globalIdxJ][isIndexJ] = velocityW;
+ this->problem().variables().velocitySecondPhase()[globalIdxJ][isIndexJ] = velocityNW;
+ break;
+ }
+ case vn:
+ {
Dune::FieldVector<Scalar,dimWorld> vel(velocityNW);
vel*=areaWeight;
- this->problem().variables().velocity()[globalIdxI][isIndex] += (vel);
+ this->problem().variables().velocity()[globalIdxI][isIndex] += (vel);
vel = velocityW;
vel*=areaWeight;
- this->problem().variables().velocitySecondPhase()[globalIdxI][isIndex] += (vel);
+ this->problem().variables().velocitySecondPhase()[globalIdxI][isIndex] += (vel);
// this->problem().variables().velocity()[globalIdxI][isIndex] = velocityNW;
// this->problem().variables().velocitySecondPhase()[globalIdxI][isIndex] = velocityW;
- this->problem().variables().velocity()[globalIdxJ][isIndexJ] = velocityNW;
- this->problem().variables().velocitySecondPhase()[globalIdxJ][isIndexJ] = velocityW;
- break;
- }
- case vt:
- {
- switch (this->pressureType)
- {
- case pw:
- {
- Dune::FieldVector<Scalar,dimWorld> vel(velocityW);
- vel+=velocityNW;
- this->problem().variables().velocity()[globalIdxI][isIndex] += (vel *=areaWeight);
- vel = velocityNW;
- this->problem().variables().velocitySecondPhase()[globalIdxI][isIndex] += (vel*=areaWeight);
+ this->problem().variables().velocity()[globalIdxJ][isIndexJ] = velocityNW;
+ this->problem().variables().velocitySecondPhase()[globalIdxJ][isIndexJ] = velocityW;
+ break;
+ }
+ case vt:
+ {
+ switch (this->pressureType)
+ {
+ case pw:
+ {
+ Dune::FieldVector<Scalar,dimWorld> vel(velocityW);
+ vel+=velocityNW;
+ this->problem().variables().velocity()[globalIdxI][isIndex] += (vel *=areaWeight);
+ vel = velocityNW;
+ this->problem().variables().velocitySecondPhase()[globalIdxI][isIndex] += (vel*=areaWeight);
// this->problem().variables().velocity()[globalIdxI][isIndex] = (velocityW+velocityNW);
// this->problem().variables().velocitySecondPhase()[globalIdxI][isIndex] = velocityNW;
- this->problem().variables().velocity()[globalIdxJ][isIndexJ] = velocityW + velocityNW;
- this->problem().variables().velocitySecondPhase()[globalIdxJ][isIndexJ] = velocityNW;
+ this->problem().variables().velocity()[globalIdxJ][isIndexJ] = velocityW + velocityNW;
+ this->problem().variables().velocitySecondPhase()[globalIdxJ][isIndexJ] = velocityNW;
- break;
- }
- case pn:
- {
+ break;
+ }
+ case pn:
+ {
Dune::FieldVector<Scalar,dimWorld> vel(velocityW);
vel+=velocityNW;
- this->problem().variables().velocity()[globalIdxI][isIndex] += (vel *=areaWeight);
- vel = velocityW;
- this->problem().variables().velocitySecondPhase()[globalIdxI][isIndex] += (vel*=areaWeight);
-// this->problem().variables().velocity()[globalIdxI][isIndex] = (velocityW+velocityNW);
-// this->problem().variables().velocitySecondPhase()[globalIdxI][isIndex] = velocityW;
- this->problem().variables().velocity()[globalIdxJ][isIndexJ] = velocityW + velocityNW;
- this->problem().variables().velocitySecondPhase()[globalIdxJ][isIndexJ] = velocityW;
- break;
- }
- }
- break;
- }
- }
+ this->problem().variables().velocity()[globalIdxI][isIndex] += (vel *=areaWeight);
+ vel = velocityW;
+ this->problem().variables().velocitySecondPhase()[globalIdxI][isIndex] += (vel*=areaWeight);
+// this->problem().variables().velocity()[globalIdxI][isIndex] = (velocityW+velocityNW);
+// this->problem().variables().velocitySecondPhase()[globalIdxI][isIndex] = velocityW;
+ this->problem().variables().velocity()[globalIdxJ][isIndexJ] = velocityW + velocityNW;
+ this->problem().variables().velocitySecondPhase()[globalIdxJ][isIndexJ] = velocityW;
+ break;
+ }
+ }
+ break;
+ }
+ }
}
}//end intersection with neighbor
diff --git a/dumux/decoupled/2p/variableclass2p.hh b/dumux/decoupled/2p/variableclass2p.hh
index d17669bc66ddb2217111c3a7f58a89e0d36a8adc..4d3f37c78d2e82825f018d1a4ca9d558712f0c9f 100644
--- a/dumux/decoupled/2p/variableclass2p.hh
+++ b/dumux/decoupled/2p/variableclass2p.hh
@@ -269,7 +269,7 @@ public:
void adaptVariableSize(int size)
{
//resize to grid size
- this->setGridSize(size);
+ this->setGridSize(size);
for (int i=0; i<2; i++) //for both phases
{
density_[i].resize(size);
diff --git a/dumux/decoupled/2p/variableclass2p_gridadapt.hh b/dumux/decoupled/2p/variableclass2p_gridadapt.hh
index 207f66863bcb5f4ca2e54695a40e7738a6af52d7..f019f86db57d78fb46c73438f6a1a516881c32b6 100644
--- a/dumux/decoupled/2p/variableclass2p_gridadapt.hh
+++ b/dumux/decoupled/2p/variableclass2p_gridadapt.hh
@@ -90,17 +90,17 @@ private:
// TODO: Doc me! Besseren Namen finden!
struct RestrictedValue
{
- Scalar saturation;
- Scalar press;
- Scalar volCorr;
- int count;
- RestrictedValue()
- {
- saturation = 0.;
- press = 0.;
- count = 0;
+ Scalar saturation;
+ Scalar press;
+ Scalar volCorr;
+ int count;
+ RestrictedValue()
+ {
+ saturation = 0.;
+ press = 0.;
+ count = 0;
volCorr = 0;
- }
+ }
};
const Grid& grid_;
@@ -114,7 +114,7 @@ public:
*/
VariableClass2PGridAdapt(const GridView& gridView) :
- VariableClass2P<TypeTag> (gridView), grid_(gridView.grid()), restrictionmap_(grid_,0) //codim 0 map
+ VariableClass2P<TypeTag> (gridView), grid_(gridView.grid()), restrictionmap_(grid_,0) //codim 0 map
{}
@@ -128,40 +128,40 @@ public:
*/
void storePrimVars(const Problem& problem)
{
- // loop over all levels of the grid
- for (int level=grid_.maxLevel(); level>=0; level--)
- {
- //get grid view on level grid
- LevelGridView levelView = grid_.levelView(level);
- for (LevelIterator it = levelView.template begin<0>();
- it!=levelView.template end<0>(); ++it)
- {
- //get your map entry
- RestrictedValue &rv = restrictionmap_[*it];
-
- // put your value in the map
- if (it->isLeaf())
- {
- // get index
- int indexI = this->index(*it);
-
- rv.saturation = this->saturation()[indexI][0];
- rv.press = this->pressure()[indexI][0];
- rv.volCorr = this->volumecorrection(indexI);
- rv.count = 1;
- }
- //Average in father
- if (it.level()>0)
- {
- ElementPointer epFather = it->father();
- RestrictedValue& rvf = restrictionmap_[*epFather];
- rvf.saturation += rv.saturation/rv.count;
- rvf.press += rv.press/rv.count;
- rvf.volCorr += rv.volCorr/rv.count;
- rvf.count += 1;
- }
- }
- }
+ // loop over all levels of the grid
+ for (int level=grid_.maxLevel(); level>=0; level--)
+ {
+ //get grid view on level grid
+ LevelGridView levelView = grid_.levelView(level);
+ for (LevelIterator it = levelView.template begin<0>();
+ it!=levelView.template end<0>(); ++it)
+ {
+ //get your map entry
+ RestrictedValue &rv = restrictionmap_[*it];
+
+ // put your value in the map
+ if (it->isLeaf())
+ {
+ // get index
+ int indexI = this->index(*it);
+
+ rv.saturation = this->saturation()[indexI][0];
+ rv.press = this->pressure()[indexI][0];
+ rv.volCorr = this->volumecorrection(indexI);
+ rv.count = 1;
+ }
+ //Average in father
+ if (it.level()>0)
+ {
+ ElementPointer epFather = it->father();
+ RestrictedValue& rvf = restrictionmap_[*epFather];
+ rvf.saturation += rv.saturation/rv.count;
+ rvf.press += rv.press/rv.count;
+ rvf.volCorr += rv.volCorr/rv.count;
+ rvf.count += 1;
+ }
+ }
+ }
}
/*!
@@ -176,53 +176,53 @@ public:
{
restrictionmap_.reserve();
- for (int level=0; level<=grid_.maxLevel(); level++)
- {
- LevelGridView levelView = grid_.levelView(level);
- for (LevelIterator it = levelView.template begin <0>();
- it!=levelView.template end <0>(); ++it)
- {
- if (!it->isNew())
- {
- //entry is in map, write in leaf
- if (it->isLeaf())
- {
- RestrictedValue &rv = restrictionmap_[*it];
- int newIdxI = this->index(*it);
- this->saturation()[newIdxI][0] = rv.saturation/rv.count;
- this->pressure()[newIdxI][0] = rv.press/rv.count;
- this->volumecorrection(newIdxI) = rv.volCorr/rv.count;
- }
- }
- else
- {
- // value is not in map, interpolate from father element
- if (it.level()>0)
- {
- ElementPointer ep = it->father();
- RestrictedValue& rvf = restrictionmap_[*ep];
- if (it->isLeaf())
- {
- int newIdxI = this->index(*it);
- this->saturation()[newIdxI][0] = rvf.saturation/rvf.count;
- this->pressure()[newIdxI][0] = rvf.press/rvf.count;
- this->volumecorrection(newIdxI) = rvf.volCorr/rvf.count;
- }
- else
- {
- //create new entry
- RestrictedValue& rv = restrictionmap_[*it];
- rv.saturation =rvf.saturation/rvf.count;
- rv.press =rvf.press/rvf.count;
- rv.volCorr =rvf.volCorr/rvf.count;
- rv.count = 1;
- }
- }
- }
- }
- }
- // reset entries in restrictionmap
- restrictionmap_.clear();
+ for (int level=0; level<=grid_.maxLevel(); level++)
+ {
+ LevelGridView levelView = grid_.levelView(level);
+ for (LevelIterator it = levelView.template begin <0>();
+ it!=levelView.template end <0>(); ++it)
+ {
+ if (!it->isNew())
+ {
+ //entry is in map, write in leaf
+ if (it->isLeaf())
+ {
+ RestrictedValue &rv = restrictionmap_[*it];
+ int newIdxI = this->index(*it);
+ this->saturation()[newIdxI][0] = rv.saturation/rv.count;
+ this->pressure()[newIdxI][0] = rv.press/rv.count;
+ this->volumecorrection(newIdxI) = rv.volCorr/rv.count;
+ }
+ }
+ else
+ {
+ // value is not in map, interpolate from father element
+ if (it.level()>0)
+ {
+ ElementPointer ep = it->father();
+ RestrictedValue& rvf = restrictionmap_[*ep];
+ if (it->isLeaf())
+ {
+ int newIdxI = this->index(*it);
+ this->saturation()[newIdxI][0] = rvf.saturation/rvf.count;
+ this->pressure()[newIdxI][0] = rvf.press/rvf.count;
+ this->volumecorrection(newIdxI) = rvf.volCorr/rvf.count;
+ }
+ else
+ {
+ //create new entry
+ RestrictedValue& rv = restrictionmap_[*it];
+ rv.saturation =rvf.saturation/rvf.count;
+ rv.press =rvf.press/rvf.count;
+ rv.volCorr =rvf.volCorr/rvf.count;
+ rv.count = 1;
+ }
+ }
+ }
+ }
+ }
+ // reset entries in restrictionmap
+ restrictionmap_.clear();
}
};
}
diff --git a/dumux/decoupled/2p2c/fvpressure2p2c.hh b/dumux/decoupled/2p2c/fvpressure2p2c.hh
index 748d6cd526df659a15643eca5f1b405a39bc2610..a5d031eb1fd67b49edac26707eef8cfcfd42380a 100644
--- a/dumux/decoupled/2p2c/fvpressure2p2c.hh
+++ b/dumux/decoupled/2p2c/fvpressure2p2c.hh
@@ -148,16 +148,16 @@ public:
//pressure solution routine: update estimate for secants, assemble, solve.
void pressure(bool solveTwice = true)
{
- //pre-transport to estimate update vector
- Scalar dt_estimate = 0.;
- Dune::dinfo << "secant guess"<< std::endl;
- problem().transportModel().update(-1, dt_estimate, problem().variables().updateEstimate(), false);
- //last argument false in update() makes shure that this is estimate and no "real" transport step
+ //pre-transport to estimate update vector
+ Scalar dt_estimate = 0.;
+ Dune::dinfo << "secant guess"<< std::endl;
+ problem().transportModel().update(-1, dt_estimate, problem().variables().updateEstimate(), false);
+ //last argument false in update() makes shure that this is estimate and no "real" transport step
problem().variables().updateEstimate() *= problem().timeManager().timeStepSize();
problem().variables().communicateUpdateEstimate();
- assemble(false); Dune::dinfo << "pressure calculation"<< std::endl;
+ assemble(false); Dune::dinfo << "pressure calculation"<< std::endl;
solve();
return;
@@ -508,20 +508,20 @@ void FVPressure2P2C<TypeTag>::assemble(bool first)
}
else
{
- // derivatives of the fluid volume with respect to concentration of components, or pressure
- if (problem().variables().dv(globalIdxI, wCompIdx) == 0)
- volumeDerivatives(globalPos, *eIt,
- problem().variables().dv(globalIdxI, wPhaseIdx),
- problem().variables().dv(globalIdxI, nPhaseIdx),
- problem().variables().dv_dp(globalIdxI));
-
- source[contiWEqIdx] *= problem().variables().dv(globalIdxI, wPhaseIdx); // note: dV_[i][1] = dv_dC1 = dV/dm1
- source[contiNEqIdx] *= problem().variables().dv(globalIdxI, nPhaseIdx);
+ // derivatives of the fluid volume with respect to concentration of components, or pressure
+ if (problem().variables().dv(globalIdxI, wCompIdx) == 0)
+ volumeDerivatives(globalPos, *eIt,
+ problem().variables().dv(globalIdxI, wPhaseIdx),
+ problem().variables().dv(globalIdxI, nPhaseIdx),
+ problem().variables().dv_dp(globalIdxI));
+
+ source[contiWEqIdx] *= problem().variables().dv(globalIdxI, wPhaseIdx); // note: dV_[i][1] = dv_dC1 = dV/dm1
+ source[contiNEqIdx] *= problem().variables().dv(globalIdxI, nPhaseIdx);
}
- f_[globalIdxI] = volume * (source[contiWEqIdx] + source[contiNEqIdx]);
+ f_[globalIdxI] = volume * (source[contiWEqIdx] + source[contiNEqIdx]);
/***********************************/
- // get absolute permeability
+ // get absolute permeability
FieldMatrix permeabilityI(problem().spatialParameters().intrinsicPermeability(globalPos, *eIt));
// get mobilities and fractional flow factors
@@ -601,11 +601,11 @@ void FVPressure2P2C<TypeTag>::assemble(bool first)
if (first) // if we are at the very first iteration we can't calculate phase potentials
{
- // get fractional flow factors in neigbor
+ // get fractional flow factors in neigbor
Scalar fractionalWJ = lambdaWJ / (lambdaWJ+ lambdaNWJ);
Scalar fractionalNWJ = lambdaNWJ / (lambdaWJ+ lambdaNWJ);
- // perform central weighting
+ // perform central weighting
Scalar lambda = (lambdaWI + lambdaWJ) * 0.5 + (lambdaNWI + lambdaNWJ) * 0.5;
entry = fabs(lambda*faceArea*(permeability*unitOuterNormal)/(dist));
@@ -615,21 +615,21 @@ void FVPressure2P2C<TypeTag>::assemble(bool first)
}
else
{
- // determine volume derivatives
+ // determine volume derivatives
if (problem().variables().dv(globalIdxJ, wCompIdx) == 0)
- volumeDerivatives(globalPosNeighbor, *neighborPointer,
- problem().variables().dv(globalIdxJ, wPhaseIdx),
- problem().variables().dv(globalIdxJ, nPhaseIdx),
- problem().variables().dv_dp(globalIdxJ));
+ volumeDerivatives(globalPosNeighbor, *neighborPointer,
+ problem().variables().dv(globalIdxJ, wPhaseIdx),
+ problem().variables().dv(globalIdxJ, nPhaseIdx),
+ problem().variables().dv_dp(globalIdxJ));
dv_dC1 = (problem().variables().dv(globalIdxJ, wPhaseIdx)
- + problem().variables().dv(globalIdxI, wPhaseIdx)) / 2; // dV/dm1= dv/dC^1
+ + problem().variables().dv(globalIdxI, wPhaseIdx)) / 2; // dV/dm1= dv/dC^1
dv_dC2 = (problem().variables().dv(globalIdxJ, nPhaseIdx)
- + problem().variables().dv(globalIdxI, nPhaseIdx)) / 2;
+ + problem().variables().dv(globalIdxI, nPhaseIdx)) / 2;
Scalar graddv_dC1 = (problem().variables().dv(globalIdxJ, wPhaseIdx)
- - problem().variables().dv(globalIdxI, wPhaseIdx)) / dist;
+ - problem().variables().dv(globalIdxI, wPhaseIdx)) / dist;
Scalar graddv_dC2 = (problem().variables().dv(globalIdxJ, nPhaseIdx)
- - problem().variables().dv(globalIdxI, nPhaseIdx)) / dist;
+ - problem().variables().dv(globalIdxI, nPhaseIdx)) / dist;
// potentialW = problem().variables().potentialWetting(globalIdxI, isIndex);
@@ -666,78 +666,78 @@ void FVPressure2P2C<TypeTag>::assemble(bool first)
potentialW += densityW * (unitDistVec * gravity_);
potentialNW += densityNW * (unitDistVec * gravity_);
- // initialize convenience shortcuts
- Scalar lambdaW, lambdaN;
- Scalar dV_w(0.), dV_n(0.); // dV_a = \sum_k \rho_a * dv/dC^k * X^k_a
- Scalar gV_w(0.), gV_n(0.); // multipaper eq(3.3) line 3 analogon dV_w
-
-
- //do the upwinding of the mobility depending on the phase potentials
- if (potentialW >= 0.)
- {
- dV_w = (dv_dC1 * problem().variables().wet_X1(globalIdxI) + dv_dC2 * (1. - problem().variables().wet_X1(globalIdxI)));
- lambdaW = problem().variables().mobilityWetting(globalIdxI);
- gV_w = (graddv_dC1 * problem().variables().wet_X1(globalIdxI) + graddv_dC2 * (1. - problem().variables().wet_X1(globalIdxI)));
- dV_w *= densityWI; gV_w *= densityWI;
- }
- else
- {
- dV_w = (dv_dC1 * problem().variables().wet_X1(globalIdxJ) + dv_dC2 * (1. - problem().variables().wet_X1(globalIdxJ)));
- lambdaW = problem().variables().mobilityWetting(globalIdxJ);
- gV_w = (graddv_dC1 * problem().variables().wet_X1(globalIdxJ) + graddv_dC2 * (1. - problem().variables().wet_X1(globalIdxJ)));
- dV_w *= densityWJ; gV_w *= densityWJ;
- }
- if (potentialNW >= 0.)
- {
- dV_n = (dv_dC1 * problem().variables().nonwet_X1(globalIdxI) + dv_dC2 * (1. - problem().variables().nonwet_X1(globalIdxI)));
- lambdaN = problem().variables().mobilityNonwetting(globalIdxI);
- gV_n = (graddv_dC1 * problem().variables().nonwet_X1(globalIdxI) + graddv_dC2 * (1. - problem().variables().nonwet_X1(globalIdxI)));
- dV_n *= densityNWI; gV_n *= densityNWI;
- }
- else
- {
- dV_n = (dv_dC1 * problem().variables().nonwet_X1(globalIdxJ) + dv_dC2 * (1. - problem().variables().nonwet_X1(globalIdxJ)));
- lambdaN = problem().variables().mobilityNonwetting(globalIdxJ);
- gV_n = (graddv_dC1 * problem().variables().nonwet_X1(globalIdxJ) + graddv_dC2 * (1. - problem().variables().nonwet_X1(globalIdxJ)));
- dV_n *= densityNWJ; gV_n *= densityNWJ;
- }
-
- //calculate current matrix entry
+ // initialize convenience shortcuts
+ Scalar lambdaW, lambdaN;
+ Scalar dV_w(0.), dV_n(0.); // dV_a = \sum_k \rho_a * dv/dC^k * X^k_a
+ Scalar gV_w(0.), gV_n(0.); // multipaper eq(3.3) line 3 analogon dV_w
+
+
+ //do the upwinding of the mobility depending on the phase potentials
+ if (potentialW >= 0.)
+ {
+ dV_w = (dv_dC1 * problem().variables().wet_X1(globalIdxI) + dv_dC2 * (1. - problem().variables().wet_X1(globalIdxI)));
+ lambdaW = problem().variables().mobilityWetting(globalIdxI);
+ gV_w = (graddv_dC1 * problem().variables().wet_X1(globalIdxI) + graddv_dC2 * (1. - problem().variables().wet_X1(globalIdxI)));
+ dV_w *= densityWI; gV_w *= densityWI;
+ }
+ else
+ {
+ dV_w = (dv_dC1 * problem().variables().wet_X1(globalIdxJ) + dv_dC2 * (1. - problem().variables().wet_X1(globalIdxJ)));
+ lambdaW = problem().variables().mobilityWetting(globalIdxJ);
+ gV_w = (graddv_dC1 * problem().variables().wet_X1(globalIdxJ) + graddv_dC2 * (1. - problem().variables().wet_X1(globalIdxJ)));
+ dV_w *= densityWJ; gV_w *= densityWJ;
+ }
+ if (potentialNW >= 0.)
+ {
+ dV_n = (dv_dC1 * problem().variables().nonwet_X1(globalIdxI) + dv_dC2 * (1. - problem().variables().nonwet_X1(globalIdxI)));
+ lambdaN = problem().variables().mobilityNonwetting(globalIdxI);
+ gV_n = (graddv_dC1 * problem().variables().nonwet_X1(globalIdxI) + graddv_dC2 * (1. - problem().variables().nonwet_X1(globalIdxI)));
+ dV_n *= densityNWI; gV_n *= densityNWI;
+ }
+ else
+ {
+ dV_n = (dv_dC1 * problem().variables().nonwet_X1(globalIdxJ) + dv_dC2 * (1. - problem().variables().nonwet_X1(globalIdxJ)));
+ lambdaN = problem().variables().mobilityNonwetting(globalIdxJ);
+ gV_n = (graddv_dC1 * problem().variables().nonwet_X1(globalIdxJ) + graddv_dC2 * (1. - problem().variables().nonwet_X1(globalIdxJ)));
+ dV_n *= densityNWJ; gV_n *= densityNWJ;
+ }
+
+ //calculate current matrix entry
entry = faceArea * (lambdaW * dV_w + lambdaN * dV_n)
- - volume * faceArea / perimeter * (lambdaW * gV_w + lambdaN * gV_n); // = boundary integral - area integral
+ - volume * faceArea / perimeter * (lambdaW * gV_w + lambdaN * gV_n); // = boundary integral - area integral
entry *= fabs((permeability*unitOuterNormal)/(dist));
//calculate right hand side
rightEntry = faceArea * (unitOuterNormal * unitDistVec) * (densityW * lambdaW * dV_w + densityNW * lambdaN * dV_n);
- rightEntry -= volume * faceArea / perimeter * (densityW * lambdaW * gV_w + densityNW * lambdaN * gV_n);
- rightEntry *= (permeability * gravity_); // = multipaper eq(3.3) line 2+3
+ rightEntry -= volume * faceArea / perimeter * (densityW * lambdaW * gV_w + densityNW * lambdaN * gV_n);
+ rightEntry *= (permeability * gravity_); // = multipaper eq(3.3) line 2+3
// include capillary pressure fluxes
- switch (pressureType)
- {
- case pw:
- {
- // calculate capillary pressure gradient
- Dune::FieldVector<Scalar, dim> pCGradient = unitDistVec;
- pCGradient *= (pcI - pcJ) / dist;
-
- //add capillary pressure term to right hand side
- rightEntry += lambdaN * dV_n * (permeability * pCGradient) * faceArea
+ switch (pressureType)
+ {
+ case pw:
+ {
+ // calculate capillary pressure gradient
+ Dune::FieldVector<Scalar, dim> pCGradient = unitDistVec;
+ pCGradient *= (pcI - pcJ) / dist;
+
+ //add capillary pressure term to right hand side
+ rightEntry += lambdaN * dV_n * (permeability * pCGradient) * faceArea
- lambdaN * gV_n * (permeability * pCGradient) * volume * faceArea / perimeter;
- break;
- }
- case pn:
- {
- // calculate capillary pressure gradient
- Dune::FieldVector<Scalar, dim> pCGradient = unitDistVec;
- pCGradient *= (pcI - pcJ) / dist;
-
- //add capillary pressure term to right hand side
- rightEntry -= lambdaW * dV_w * (permeability * pCGradient) * faceArea
+ break;
+ }
+ case pn:
+ {
+ // calculate capillary pressure gradient
+ Dune::FieldVector<Scalar, dim> pCGradient = unitDistVec;
+ pCGradient *= (pcI - pcJ) / dist;
+
+ //add capillary pressure term to right hand side
+ rightEntry -= lambdaW * dV_w * (permeability * pCGradient) * faceArea
- lambdaW * gV_w * (permeability * pCGradient) * volume * faceArea / perimeter;
- break;
- }
- }
+ break;
+ }
+ }
} // end !first
//set right hand side
@@ -753,7 +753,7 @@ void FVPressure2P2C<TypeTag>::assemble(bool first)
/************* boundary face ************************/
else
{
- // get volume derivatives inside the cell
+ // get volume derivatives inside the cell
dv_dC1 = problem().variables().dv(globalIdxI, wCompIdx);
dv_dC2 = problem().variables().dv(globalIdxI, nCompIdx);
@@ -781,8 +781,8 @@ void FVPressure2P2C<TypeTag>::assemble(bool first)
Dune::FieldVector<Scalar, dim> permeability(0);
permeabilityI.mv(unitDistVec, permeability);
- // create a fluid state for the boundary
- FluidState BCfluidState;
+ // create a fluid state for the boundary
+ FluidState BCfluidState;
Scalar temperatureBC = problem().temperatureAtPos(globalPosFace);
//read boundary values
@@ -791,7 +791,7 @@ void FVPressure2P2C<TypeTag>::assemble(bool first)
if (first)
{
- Scalar lambda = lambdaWI+lambdaNWI;
+ Scalar lambda = lambdaWI+lambdaNWI;
A_[globalIdxI][globalIdxI] += lambda * faceArea * (permeability * unitOuterNormal) / (dist);
Scalar pressBC = primaryVariablesOnBoundary[Indices::pressureEqIdx];
f_[globalIdxI] += lambda * faceArea * pressBC * (permeability * unitOuterNormal) / (dist);
@@ -890,7 +890,7 @@ void FVPressure2P2C<TypeTag>::assemble(bool first)
//do the upwinding of the mobility depending on the phase potentials
Scalar lambdaW, lambdaNW;
- Scalar dV_w, dV_n; // gV_a weglassen, da dV/dc am Rand ortsunabhängig angenommen -> am rand nicht bestimmbar -> nur Randintegral ohne Gebietsintegral
+ Scalar dV_w, dV_n; // gV_a weglassen, da dV/dc am Rand ortsunabhängig angenommen -> am rand nicht bestimmbar -> nur Randintegral ohne Gebietsintegral
if (potentialW >= 0.)
{
@@ -904,7 +904,7 @@ void FVPressure2P2C<TypeTag>::assemble(bool first)
{
densityW = densityWBound;
dV_w = (dv_dC1 * BCfluidState.massFrac(wPhaseIdx, wCompIdx)
- + dv_dC2 * BCfluidState.massFrac(wPhaseIdx, nCompIdx));
+ + dv_dC2 * BCfluidState.massFrac(wPhaseIdx, nCompIdx));
dV_w *= densityW;
lambdaW = lambdaWBound;
}
@@ -919,7 +919,7 @@ void FVPressure2P2C<TypeTag>::assemble(bool first)
{
densityNW = densityNWBound;
dV_n = (dv_dC1 * BCfluidState.massFrac(nPhaseIdx, wCompIdx)
- + dv_dC2 * BCfluidState.massFrac(nPhaseIdx, nCompIdx));
+ + dv_dC2 * BCfluidState.massFrac(nPhaseIdx, nCompIdx));
dV_n *= densityNW;
lambdaNW = lambdaNWBound;
}
@@ -963,7 +963,7 @@ void FVPressure2P2C<TypeTag>::assemble(bool first)
A_[globalIdxI][globalIdxI] += entry;
f_[globalIdxI] += entry * primaryVariablesOnBoundary[Indices::pressureEqIdx];
f_[globalIdxI] -= rightEntry * (unitOuterNormal * unitDistVec);
- } //end of if(first) ... else{...
+ } //end of if(first) ... else{...
} // end dirichlet
/**********************************
@@ -971,18 +971,18 @@ void FVPressure2P2C<TypeTag>::assemble(bool first)
**********************************/
else if(bcType.isNeumann(Indices::pressureEqIdx))
{
- PrimaryVariables J(NAN);
- problem().neumann(J, *isIt);
- if (first)
- {
- J[contiWEqIdx] /= densityWI;
- J[contiNEqIdx] /= densityNWI;
- }
- else
- {
- J[contiWEqIdx] *= dv_dC1;
- J[contiNEqIdx] *= dv_dC2;
- }
+ PrimaryVariables J(NAN);
+ problem().neumann(J, *isIt);
+ if (first)
+ {
+ J[contiWEqIdx] /= densityWI;
+ J[contiNEqIdx] /= densityNWI;
+ }
+ else
+ {
+ J[contiWEqIdx] *= dv_dC1;
+ J[contiNEqIdx] *= dv_dC2;
+ }
f_[globalIdxI] -= (J[contiWEqIdx] + J[contiNEqIdx]) * faceArea;
}
@@ -1023,13 +1023,13 @@ void FVPressure2P2C<TypeTag>::assemble(bool first)
{
if (erri <= x_mi * maxErr)
f_[globalIdxI] +=
- problem().variables().errorCorrection(globalIdxI) =
- fac* (1-x_mi*(lofac-1)/(x_lo-x_mi) + (lofac-1)/(x_lo-x_mi)*erri/maxErr)
+ problem().variables().errorCorrection(globalIdxI) =
+ fac* (1-x_mi*(lofac-1)/(x_lo-x_mi) + (lofac-1)/(x_lo-x_mi)*erri/maxErr)
* problem().variables().volErr()[globalIdxI] * volume;
else
f_[globalIdxI] +=
- problem().variables().errorCorrection(globalIdxI) =
- fac * (1 + x_mi - hifac*x_mi/(1-x_mi) + (hifac/(1-x_mi)-1)*erri/maxErr)
+ problem().variables().errorCorrection(globalIdxI) =
+ fac * (1 + x_mi - hifac*x_mi/(1-x_mi) + (hifac/(1-x_mi)-1)*erri/maxErr)
* problem().variables().volErr()[globalIdxI] * volume;
}
// printmatrix(std::cout, A_, "global stiffness matrix", "row", 11, 3);
@@ -1073,10 +1073,10 @@ void FVPressure2P2C<TypeTag>::initialMaterialLaws(bool compositional)
problem().variables().communicateTransportedQuantity();
problem().variables().communicatePressure();
- // initialize the fluid system
+ // initialize the fluid system
FluidState fluidState;
- // iterate through leaf grid an evaluate c0 at cell center
+ // iterate through leaf grid an evaluate c0 at cell center
ElementIterator eItEnd = problem().gridView().template end<0>();
ElementIterator eIt = problem().gridView().template begin<0>();
for (; eIt != eItEnd; ++eIt)
@@ -1091,45 +1091,45 @@ void FVPressure2P2C<TypeTag>::initialMaterialLaws(bool compositional)
Scalar temperature_ = problem().temperatureAtPos(globalPos);
// initial conditions
- PhaseVector pressure(0.);
- Scalar sat_0=0.;
+ PhaseVector pressure(0.);
+ Scalar sat_0=0.;
- typename Indices::BoundaryFormulation icFormulation;
- problem().initialFormulation(icFormulation, *eIt); // get type of initial condition
+ typename Indices::BoundaryFormulation icFormulation;
+ problem().initialFormulation(icFormulation, *eIt); // get type of initial condition
if(!compositional) //means that we do the first approximate guess without compositions
{
- // phase pressures are unknown, so start with an exemplary
- Scalar exemplaryPressure = problem().referencePressure(*eIt);
- pressure[wPhaseIdx] = pressure[nPhaseIdx] = problem().variables().pressure()[globalIdx] = exemplaryPressure;
- problem().variables().capillaryPressure(globalIdx) = 0.;
- if (icFormulation == Indices::BoundaryFormulation::saturation) // saturation initial condition
- {
- sat_0 = problem().initSat(*eIt);
- fluidState.satFlash(sat_0, pressure, problem().spatialParameters().porosity(globalPos, *eIt), temperature_);
- }
- else if (icFormulation == Indices::BoundaryFormulation::concentration) // concentration initial condition
- {
- Scalar Z1_0 = problem().initConcentration(*eIt);
- fluidState.update(Z1_0, pressure, problem().spatialParameters().porosity(globalPos, *eIt), temperature_);
- }
+ // phase pressures are unknown, so start with an exemplary
+ Scalar exemplaryPressure = problem().referencePressure(*eIt);
+ pressure[wPhaseIdx] = pressure[nPhaseIdx] = problem().variables().pressure()[globalIdx] = exemplaryPressure;
+ problem().variables().capillaryPressure(globalIdx) = 0.;
+ if (icFormulation == Indices::BoundaryFormulation::saturation) // saturation initial condition
+ {
+ sat_0 = problem().initSat(*eIt);
+ fluidState.satFlash(sat_0, pressure, problem().spatialParameters().porosity(globalPos, *eIt), temperature_);
+ }
+ else if (icFormulation == Indices::BoundaryFormulation::concentration) // concentration initial condition
+ {
+ Scalar Z1_0 = problem().initConcentration(*eIt);
+ fluidState.update(Z1_0, pressure, problem().spatialParameters().porosity(globalPos, *eIt), temperature_);
+ }
}
- else if(compositional) //means we regard compositional effects since we know an estimate pressure field
+ else if(compositional) //means we regard compositional effects since we know an estimate pressure field
{
- if (icFormulation == Indices::BoundaryFormulation::saturation) // saturation initial condition
- {
- //get saturation, determine pc
- sat_0 = problem().initSat(*eIt);
- if(GET_PROP_VALUE(TypeTag, PTAG(EnableCapillarity)))
- {
+ if (icFormulation == Indices::BoundaryFormulation::saturation) // saturation initial condition
+ {
+ //get saturation, determine pc
+ sat_0 = problem().initSat(*eIt);
+ if(GET_PROP_VALUE(TypeTag, PTAG(EnableCapillarity)))
+ {
problem().variables().capillaryPressure(globalIdx)
= MaterialLaw::pC(problem().spatialParameters().materialLawParams(globalPos, *eIt),
sat_0);
- }
- else
- problem().variables().capillaryPressure(globalIdx) = 0.;
+ }
+ else
+ problem().variables().capillaryPressure(globalIdx) = 0.;
- //determine phase pressures from primary pressure variable
+ //determine phase pressures from primary pressure variable
switch (pressureType)
{
case pw:
@@ -1146,66 +1146,66 @@ void FVPressure2P2C<TypeTag>::initialMaterialLaws(bool compositional)
}
}
- fluidState.satFlash(sat_0, pressure, problem().spatialParameters().porosity(globalPos, *eIt), temperature_);
- }
- else if (icFormulation == Indices::BoundaryFormulation::concentration) // concentration initial condition
- {
- Scalar Z1_0 = problem().initConcentration(*eIt);
- // If total concentrations are given at the boundary, saturation is unknown.
- // This may affect pc and hence p_alpha and hence again saturation -> iteration.
-
- // iterations in case of enabled capillary pressure
- if(GET_PROP_VALUE(TypeTag, PTAG(EnableCapillarity)))
- {
- //start with pc from last TS
- Scalar pc(problem().variables().capillaryPressure(globalIdx));
-
- int maxiter = 3;
- //start iteration loop
- for(int iter=0; iter < maxiter; iter++)
- {
- //determine phase pressures from primary pressure variable
- switch (pressureType)
- {
- case pw:
- {
- pressure[wPhaseIdx] = problem().variables().pressure()[globalIdx];
- pressure[nPhaseIdx] = problem().variables().pressure()[globalIdx] + pc;
- break;
- }
- case pn:
- {
- pressure[wPhaseIdx] = problem().variables().pressure()[globalIdx] - pc;
- pressure[nPhaseIdx] = problem().variables().pressure()[globalIdx];
- break;
- }
- }
-
- //store old pc
- Scalar oldPc = pc;
- //update with better pressures
- fluidState.update(Z1_0, pressure, problem().spatialParameters().porosity(globalPos, *eIt),
- problem().temperatureAtPos(globalPos));
- pc = MaterialLaw::pC(problem().spatialParameters().materialLawParams(globalPos, *eIt),
- fluidState.saturation(wPhaseIdx));
- // TODO: get right criterion, do output for evaluation
- //converge criterion
- if (abs(oldPc-pc)<10)
- iter = maxiter;
+ fluidState.satFlash(sat_0, pressure, problem().spatialParameters().porosity(globalPos, *eIt), temperature_);
+ }
+ else if (icFormulation == Indices::BoundaryFormulation::concentration) // concentration initial condition
+ {
+ Scalar Z1_0 = problem().initConcentration(*eIt);
+ // If total concentrations are given at the boundary, saturation is unknown.
+ // This may affect pc and hence p_alpha and hence again saturation -> iteration.
+
+ // iterations in case of enabled capillary pressure
+ if(GET_PROP_VALUE(TypeTag, PTAG(EnableCapillarity)))
+ {
+ //start with pc from last TS
+ Scalar pc(problem().variables().capillaryPressure(globalIdx));
+
+ int maxiter = 3;
+ //start iteration loop
+ for(int iter=0; iter < maxiter; iter++)
+ {
+ //determine phase pressures from primary pressure variable
+ switch (pressureType)
+ {
+ case pw:
+ {
+ pressure[wPhaseIdx] = problem().variables().pressure()[globalIdx];
+ pressure[nPhaseIdx] = problem().variables().pressure()[globalIdx] + pc;
+ break;
+ }
+ case pn:
+ {
+ pressure[wPhaseIdx] = problem().variables().pressure()[globalIdx] - pc;
+ pressure[nPhaseIdx] = problem().variables().pressure()[globalIdx];
+ break;
+ }
+ }
+
+ //store old pc
+ Scalar oldPc = pc;
+ //update with better pressures
+ fluidState.update(Z1_0, pressure, problem().spatialParameters().porosity(globalPos, *eIt),
+ problem().temperatureAtPos(globalPos));
+ pc = MaterialLaw::pC(problem().spatialParameters().materialLawParams(globalPos, *eIt),
+ fluidState.saturation(wPhaseIdx));
+ // TODO: get right criterion, do output for evaluation
+ //converge criterion
+ if (abs(oldPc-pc)<10)
+ iter = maxiter;
pc = MaterialLaw::pC(problem().spatialParameters().materialLawParams(globalPos, *eIt),
fluidState.saturation(wPhaseIdx));
- }
- problem().variables().capillaryPressure(globalIdx) = pc; //complete iteration procedure
- }
- else // capillary pressure neglected
- {
- problem().variables().capillaryPressure(globalIdx) = 0.;
- pressure[wPhaseIdx] = pressure[nPhaseIdx]
- = problem().variables().pressure()[globalIdx];
- fluidState.update(Z1_0, pressure, problem().spatialParameters().porosity(globalPos, *eIt), temperature_);
- }
- } //end conc initial condition
+ }
+ problem().variables().capillaryPressure(globalIdx) = pc; //complete iteration procedure
+ }
+ else // capillary pressure neglected
+ {
+ problem().variables().capillaryPressure(globalIdx) = 0.;
+ pressure[wPhaseIdx] = pressure[nPhaseIdx]
+ = problem().variables().pressure()[globalIdx];
+ fluidState.update(Z1_0, pressure, problem().spatialParameters().porosity(globalPos, *eIt), temperature_);
+ }
+ } //end conc initial condition
} //end compositional
// initialize densities
@@ -1247,7 +1247,7 @@ void FVPressure2P2C<TypeTag>::updateMaterialLaws()
problem().variables().communicateTransportedQuantity();
problem().variables().communicatePressure();
- // instantiate a brandnew fluid state object
+ // instantiate a brandnew fluid state object
FluidState fluidState;
//get timestep for error term
@@ -1368,32 +1368,32 @@ void FVPressure2P2C<TypeTag>::updateMaterialLaws()
/******** update variables in variableclass **********/
// initialize saturation, capillary pressure
- problem().variables().saturation(globalIdx) = fluidState.saturation(wPhaseIdx);
- problem().variables().capillaryPressure(globalIdx) = pc; // pc=0 if EnableCapillarity=false
+ problem().variables().saturation(globalIdx) = fluidState.saturation(wPhaseIdx);
+ problem().variables().capillaryPressure(globalIdx) = pc; // pc=0 if EnableCapillarity=false
// initialize viscosities
- problem().variables().viscosityWetting(globalIdx)
- = FluidSystem::phaseViscosity(wPhaseIdx, temperature_, pressure[wPhaseIdx], fluidState);
- problem().variables().viscosityNonwetting(globalIdx)
- = FluidSystem::phaseViscosity(nPhaseIdx, temperature_, pressure[nPhaseIdx], fluidState);
+ problem().variables().viscosityWetting(globalIdx)
+ = FluidSystem::phaseViscosity(wPhaseIdx, temperature_, pressure[wPhaseIdx], fluidState);
+ problem().variables().viscosityNonwetting(globalIdx)
+ = FluidSystem::phaseViscosity(nPhaseIdx, temperature_, pressure[nPhaseIdx], fluidState);
// initialize mobilities
- problem().variables().mobilityWetting(globalIdx) =
+ problem().variables().mobilityWetting(globalIdx) =
MaterialLaw::krw(problem().spatialParameters().materialLawParams(globalPos, *eIt), fluidState.saturation(wPhaseIdx))
/ problem().variables().viscosityWetting(globalIdx);
- problem().variables().mobilityNonwetting(globalIdx) =
+ problem().variables().mobilityNonwetting(globalIdx) =
MaterialLaw::krn(problem().spatialParameters().materialLawParams(globalPos, *eIt), fluidState.saturation(wPhaseIdx))
/ problem().variables().viscosityNonwetting(globalIdx);
// initialize mass fractions
- problem().variables().wet_X1(globalIdx) = fluidState.massFrac(wPhaseIdx, wCompIdx);
- problem().variables().nonwet_X1(globalIdx) = fluidState.massFrac(nPhaseIdx, wCompIdx);
+ problem().variables().wet_X1(globalIdx) = fluidState.massFrac(wPhaseIdx, wCompIdx);
+ problem().variables().nonwet_X1(globalIdx) = fluidState.massFrac(nPhaseIdx, wCompIdx);
// initialize densities
problem().variables().densityWetting(globalIdx)
- = FluidSystem::phaseDensity(wPhaseIdx, temperature_, pressure[wPhaseIdx], fluidState);
+ = FluidSystem::phaseDensity(wPhaseIdx, temperature_, pressure[wPhaseIdx], fluidState);
problem().variables().densityNonwetting(globalIdx)
- = FluidSystem::phaseDensity(nPhaseIdx, temperature_, pressure[nPhaseIdx], fluidState);
+ = FluidSystem::phaseDensity(nPhaseIdx, temperature_, pressure[nPhaseIdx], fluidState);
// determine volume mismatch between actual fluid volume and pore volume
Scalar sumConc = (problem().variables().totalConcentration(globalIdx, wCompIdx)
@@ -1402,7 +1402,7 @@ void FVPressure2P2C<TypeTag>::updateMaterialLaws()
Scalar massn = problem().variables().numericalDensity(globalIdx, nPhaseIdx) = sumConc * fluidState.phaseMassFraction(nPhaseIdx);
if ((problem().variables().densityWetting(globalIdx)*problem().variables().densityNonwetting(globalIdx)) == 0)
- DUNE_THROW(Dune::MathError, "Decoupled2p2c::postProcessUpdate: try to divide by 0 density");
+ DUNE_THROW(Dune::MathError, "Decoupled2p2c::postProcessUpdate: try to divide by 0 density");
Scalar vol = massw / problem().variables().densityWetting(globalIdx) + massn / problem().variables().densityNonwetting(globalIdx);
if (dt != 0)
{
@@ -1441,7 +1441,7 @@ void FVPressure2P2C<TypeTag>::updateMaterialLaws()
template<class TypeTag>
void FVPressure2P2C<TypeTag>::volumeDerivatives(GlobalPosition globalPos, ElementPointer ep, Scalar& dv_dC1, Scalar& dv_dC2, Scalar& dv_dp)
{
- // cell index
+ // cell index
int globalIdx = problem().variables().index(*ep);
// get cell temperature
@@ -1450,10 +1450,10 @@ void FVPressure2P2C<TypeTag>::volumeDerivatives(GlobalPosition globalPos, Elemen
// initialize an Fluid state for the update
FluidState updFluidState;
- /**********************************
- * a) get necessary variables
- **********************************/
- //determine phase pressures from primary pressure variable
+ /**********************************
+ * a) get necessary variables
+ **********************************/
+ //determine phase pressures from primary pressure variable
PhaseVector pressure(0.);
switch (pressureType)
{
@@ -1482,18 +1482,18 @@ void FVPressure2P2C<TypeTag>::volumeDerivatives(GlobalPosition globalPos, Elemen
// actual fluid volume
Scalar volalt = (m1+m2) * (nuw1 * v_w + (1-nuw1) * v_g);
- /**********************************
- * b) define increments
- **********************************/
+ /**********************************
+ * b) define increments
+ **********************************/
// increments for numerical derivatives
Scalar inc1 = (fabs(problem().variables().updateEstimate(globalIdx, wCompIdx)) > 1e-8 / v_w) ? problem().variables().updateEstimate(globalIdx,wCompIdx) : 1e-8/v_w;
Scalar inc2 =(fabs(problem().variables().updateEstimate(globalIdx, nCompIdx)) > 1e-8 / v_g) ? problem().variables().updateEstimate(globalIdx,nCompIdx) : 1e-8 / v_g;
Scalar incp = 1e-2;
- /**********************************
- * c) Secant method for derivatives
- **********************************/
+ /**********************************
+ * c) Secant method for derivatives
+ **********************************/
// numerical derivative of fluid volume with respect to pressure
PhaseVector p_(incp);
@@ -1525,17 +1525,17 @@ void FVPressure2P2C<TypeTag>::volumeDerivatives(GlobalPosition globalPos, Elemen
// numerical derivative of fluid volume with respect to mass of component 1
m1 += inc1;
Z1 = m1 / (m1 + m2);
- updFluidState.update(Z1, pressure, problem().spatialParameters().porosity(globalPos, *ep), temperature_);
- Scalar satt = updFluidState.saturation(wPhaseIdx);
- Scalar nuw = satt / v_w / (satt/v_w + (1-satt)/v_g);
+ updFluidState.update(Z1, pressure, problem().spatialParameters().porosity(globalPos, *ep), temperature_);
+ Scalar satt = updFluidState.saturation(wPhaseIdx);
+ Scalar nuw = satt / v_w / (satt/v_w + (1-satt)/v_g);
dv_dC1 = ((m1+m2) * (nuw * v_w + (1-nuw) * v_g) - volalt) /inc1;
m1 -= inc1;
// numerical derivative of fluid volume with respect to mass of component 2
m2 += inc2;
Z1 = m1 / (m1 + m2);
- updFluidState.update(Z1, pressure, problem().spatialParameters().porosity(globalPos, *ep), temperature_);
- satt = updFluidState.saturation(wPhaseIdx);
+ updFluidState.update(Z1, pressure, problem().spatialParameters().porosity(globalPos, *ep), temperature_);
+ satt = updFluidState.saturation(wPhaseIdx);
nuw = satt / v_w / (satt/v_w + (1-satt)/v_g);
dv_dC2 = ((m1+m2) * (nuw * v_w + (1-nuw) * v_g) - volalt)/ inc2;
m2 -= inc2;
diff --git a/dumux/decoupled/2p2c/fvpressure2p2cmultiphysics.hh b/dumux/decoupled/2p2c/fvpressure2p2cmultiphysics.hh
index cc4efa21e632187c4d0ce263e9434b5ce370d95d..45c48cef9121c34552fad26554fc01b28c26330d 100644
--- a/dumux/decoupled/2p2c/fvpressure2p2cmultiphysics.hh
+++ b/dumux/decoupled/2p2c/fvpressure2p2cmultiphysics.hh
@@ -256,15 +256,15 @@ void FVPressure2P2CMultiPhysics<TypeTag>::assemble(bool first)
}
else
{
- // derivatives of the fluid volume with respect to concentration of components, or pressure
- if (problem().variables().dv(globalIdxI, wPhaseIdx) == 0)
- this->volumeDerivatives(globalPos, *eIt,
- problem().variables().dv(globalIdxI, wPhaseIdx),
- problem().variables().dv(globalIdxI, nPhaseIdx),
- problem().variables().dv_dp(globalIdxI));
-
- source[Indices::contiWEqIdx] *= problem().variables().dv(globalIdxI, wPhaseIdx); // note: dV_[i][1] = dv_dC1 = dV/dm1
- source[Indices::contiNEqIdx] *= problem().variables().dv(globalIdxI, nPhaseIdx);
+ // derivatives of the fluid volume with respect to concentration of components, or pressure
+ if (problem().variables().dv(globalIdxI, wPhaseIdx) == 0)
+ this->volumeDerivatives(globalPos, *eIt,
+ problem().variables().dv(globalIdxI, wPhaseIdx),
+ problem().variables().dv(globalIdxI, nPhaseIdx),
+ problem().variables().dv_dp(globalIdxI));
+
+ source[Indices::contiWEqIdx] *= problem().variables().dv(globalIdxI, wPhaseIdx); // note: dV_[i][1] = dv_dC1 = dV/dm1
+ source[Indices::contiNEqIdx] *= problem().variables().dv(globalIdxI, nPhaseIdx);
}
this->f_[globalIdxI] = volume * (source[Indices::contiWEqIdx] + source[Indices::contiNEqIdx]);
/********************************************************************/
@@ -408,19 +408,19 @@ void FVPressure2P2CMultiPhysics<TypeTag>::assemble(bool first)
{
// determine volume derivatives
if (problem().variables().dv(globalIdxJ, wPhaseIdx) == 0)
- this->volumeDerivatives(globalPosNeighbor, *neighborPointer,
- problem().variables().dv(globalIdxJ, wPhaseIdx),
- problem().variables().dv(globalIdxJ, nPhaseIdx),
- problem().variables().dv_dp(globalIdxJ));
+ this->volumeDerivatives(globalPosNeighbor, *neighborPointer,
+ problem().variables().dv(globalIdxJ, wPhaseIdx),
+ problem().variables().dv(globalIdxJ, nPhaseIdx),
+ problem().variables().dv_dp(globalIdxJ));
dv_dC1 = (problem().variables().dv(globalIdxJ, wPhaseIdx)
- + problem().variables().dv(globalIdxI, wPhaseIdx)) / 2; // dV/dm1= dV/dC^1
+ + problem().variables().dv(globalIdxI, wPhaseIdx)) / 2; // dV/dm1= dV/dC^1
dv_dC2 = (problem().variables().dv(globalIdxJ, nPhaseIdx)
- + problem().variables().dv(globalIdxI, nPhaseIdx)) / 2;
+ + problem().variables().dv(globalIdxI, nPhaseIdx)) / 2;
graddv_dC1 = (problem().variables().dv(globalIdxJ, wPhaseIdx)
- + problem().variables().dv(globalIdxI, wPhaseIdx)) / dist;
+ + problem().variables().dv(globalIdxI, wPhaseIdx)) / dist;
graddv_dC2 = (problem().variables().dv(globalIdxJ, nPhaseIdx)
- + problem().variables().dv(globalIdxI, nPhaseIdx)) / dist;
+ + problem().variables().dv(globalIdxI, nPhaseIdx)) / dist;
}
else
{
@@ -570,7 +570,7 @@ void FVPressure2P2CMultiPhysics<TypeTag>::assemble(bool first)
/************* boundary face ************************/
else
{
- // get volume derivatives inside the cell
+ // get volume derivatives inside the cell
dv_dC1 = problem().variables().dv(globalIdxI, wCompIdx);
dv_dC2 = problem().variables().dv(globalIdxI, nCompIdx);
@@ -927,13 +927,13 @@ void FVPressure2P2CMultiPhysics<TypeTag>::assemble(bool first)
{
if (erri <= x_mi * maxErr)
this->f_[globalIdxI] +=
- problem().variables().errorCorrection(globalIdxI) =
- fac* (1-x_mi*(lofac-1)/(x_lo-x_mi) + (lofac-1)/(x_lo-x_mi)*erri/maxErr)
+ problem().variables().errorCorrection(globalIdxI) =
+ fac* (1-x_mi*(lofac-1)/(x_lo-x_mi) + (lofac-1)/(x_lo-x_mi)*erri/maxErr)
* problem().variables().volErr()[globalIdxI] * volume;
else
this->f_[globalIdxI] +=
- problem().variables().errorCorrection(globalIdxI) =
- fac * (1 + x_mi - hifac*x_mi/(1-x_mi) + (hifac/(1-x_mi)-1)*erri/maxErr)
+ problem().variables().errorCorrection(globalIdxI) =
+ fac * (1 + x_mi - hifac*x_mi/(1-x_mi) + (hifac/(1-x_mi)-1)*erri/maxErr)
* problem().variables().volErr()[globalIdxI] * volume;
}
} // end grid traversal
diff --git a/dumux/decoupled/2p2c/fvtransport2p2c.hh b/dumux/decoupled/2p2c/fvtransport2p2c.hh
index fcfd85be8f16bce2d6622299ae3a2e54dedafbf3..2ecf3a62d781516601898b11eee5fac1bffe9b4c 100644
--- a/dumux/decoupled/2p2c/fvtransport2p2c.hh
+++ b/dumux/decoupled/2p2c/fvtransport2p2c.hh
@@ -290,8 +290,8 @@ void FVTransport2P2C<TypeTag>::update(const Scalar t, Scalar& dt, TransportSolut
// compute mean permeability
Dune::FieldMatrix<Scalar,dim,dim> meanK_(0.);
Dumux::harmonicMeanMatrix(meanK_,
- K_I,
- problem().spatialParameters().intrinsicPermeability(globalPosNeighbor, *neighborPointer));
+ K_I,
+ problem().spatialParameters().intrinsicPermeability(globalPosNeighbor, *neighborPointer));
Dune::FieldVector<Scalar,dim> K(0);
meanK_.umv(unitDistVec,K);
@@ -355,7 +355,7 @@ void FVTransport2P2C<TypeTag>::update(const Scalar t, Scalar& dt, TransportSolut
}
/******************************************
- * Boundary Face
+ * Boundary Face
******************************************/
if (isIt->boundary())
{
@@ -398,15 +398,15 @@ void FVTransport2P2C<TypeTag>::update(const Scalar t, Scalar& dt, TransportSolut
Scalar Xw1Bound = BCfluidState.massFrac(wPhaseIdx, wCompIdx);
Scalar Xn1Bound = BCfluidState.massFrac(nPhaseIdx, wCompIdx);
Scalar densityWBound = BCfluidState.density(wPhaseIdx);
- Scalar densityNWBound = BCfluidState.density(nPhaseIdx);
- Scalar viscosityWBound = FluidSystem::phaseViscosity(wPhaseIdx,
- problem().temperatureAtPos(globalPosFace),
- pressBound[wPhaseIdx], BCfluidState);
- Scalar viscosityNWBound = FluidSystem::phaseViscosity(nPhaseIdx,
- problem().temperatureAtPos(globalPosFace),
- pressBound[wPhaseIdx], BCfluidState);
- if(GET_PROP_VALUE(TypeTag, PTAG(EnableCapillarity)))
- pcBound = BCfluidState.capillaryPressure();
+ Scalar densityNWBound = BCfluidState.density(nPhaseIdx);
+ Scalar viscosityWBound = FluidSystem::phaseViscosity(wPhaseIdx,
+ problem().temperatureAtPos(globalPosFace),
+ pressBound[wPhaseIdx], BCfluidState);
+ Scalar viscosityNWBound = FluidSystem::phaseViscosity(nPhaseIdx,
+ problem().temperatureAtPos(globalPosFace),
+ pressBound[wPhaseIdx], BCfluidState);
+ if(GET_PROP_VALUE(TypeTag, PTAG(EnableCapillarity)))
+ pcBound = BCfluidState.capillaryPressure();
// average
double densityW_mean = (densityWI + densityWBound) / 2;
double densityNW_mean = (densityNWI + densityNWBound) / 2;
@@ -499,23 +499,23 @@ void FVTransport2P2C<TypeTag>::update(const Scalar t, Scalar& dt, TransportSolut
updFactor[nCompIdx] = - J[Indices::contiNEqIdx] * faceArea / volume;
// for timestep control
- #define cflIgnoresNeumann
- #ifdef cflIgnoresNeumann
+ #define cflIgnoresNeumann
+ #ifdef cflIgnoresNeumann
factor[0] = 0;
factor[1] = 0;
- #else
+ #else
double inflow = updFactor[wCompIdx] / densityW + updFactor[nCompIdx] / densityNW;
if (inflow>0)
- {
- factor[0] = updFactor[wCompIdx] / densityW + updFactor[nCompIdx] / densityNW; // =factor in
- factor[1] = -(updFactor[wCompIdx] / densityW /SwmobI + updFactor[nCompIdx] / densityNW / SnmobI); // =factor out
- }
+ {
+ factor[0] = updFactor[wCompIdx] / densityW + updFactor[nCompIdx] / densityNW; // =factor in
+ factor[1] = -(updFactor[wCompIdx] / densityW /SwmobI + updFactor[nCompIdx] / densityNW / SnmobI); // =factor out
+ }
else
{
- factor[0] = -(updFactor[wCompIdx] / densityW + updFactor[nCompIdx] / densityNW); // =factor in
- factor[1] = updFactor[wCompIdx] / densityW /SwmobI + updFactor[nCompIdx] / densityNW / SnmobI; // =factor out
+ factor[0] = -(updFactor[wCompIdx] / densityW + updFactor[nCompIdx] / densityNW); // =factor in
+ factor[1] = updFactor[wCompIdx] / densityW /SwmobI + updFactor[nCompIdx] / densityNW / SnmobI; // =factor out
}
- #endif
+ #endif
}//end neumann boundary
}//end boundary
// correct update Factor by volume error
@@ -533,7 +533,7 @@ void FVTransport2P2C<TypeTag>::update(const Scalar t, Scalar& dt, TransportSolut
}// end all intersections
- /*********** Handle source term ***************/
+ /*********** Handle source term ***************/
PrimaryVariables q(NAN);
problem().source(q, *eIt);
updateVec[wCompIdx][globalIdxI] += q[Indices::contiWEqIdx];
@@ -541,7 +541,7 @@ void FVTransport2P2C<TypeTag>::update(const Scalar t, Scalar& dt, TransportSolut
// account for porosity
sumfactorin = std::max(sumfactorin,sumfactorout)
- / problem().spatialParameters().porosity(globalPos, *eIt);
+ / problem().spatialParameters().porosity(globalPos, *eIt);
if ( 1./sumfactorin < dt)
{
diff --git a/dumux/material/MpNcfluidsystems/h2on2fluidsystem.hh b/dumux/material/MpNcfluidsystems/h2on2fluidsystem.hh
index 328662b0a782017466f540231d6e790f7d2a94d1..a29638048487099c45f533a9840b246eea2f7164 100644
--- a/dumux/material/MpNcfluidsystems/h2on2fluidsystem.hh
+++ b/dumux/material/MpNcfluidsystems/h2on2fluidsystem.hh
@@ -471,7 +471,7 @@ public:
const ParameterCache ¶mCache,
int phaseIdx)
{
-// TODO thermal conductivity is a function of:
+// TODO thermal conductivity is a function of:
// Scalar p = fluidState.pressure(phaseIdx);
// Scalar T = fluidState.temperature(phaseIdx);
// Scalar x = fluidState.moleFrac(phaseIdx,compIdx);
diff --git a/test/decoupled/2padaptive/test_impes_adaptive_problem.hh b/test/decoupled/2padaptive/test_impes_adaptive_problem.hh
index dfc5fb508b7d0bd2edc3628b2d91926c6ecf49e7..354ec108c75dbd98bb9332495b3702f87c6332c4 100644
--- a/test/decoupled/2padaptive/test_impes_adaptive_problem.hh
+++ b/test/decoupled/2padaptive/test_impes_adaptive_problem.hh
@@ -66,7 +66,7 @@ NEW_TYPE_TAG(TestIMPESAdaptiveProblem, INHERITS_FROM(DecoupledTwoP, Transport, T
// Set the grid type
SET_PROP(TestIMPESAdaptiveProblem, Grid)
{
- typedef Dune::UGGrid<2> type;
+ typedef Dune::UGGrid<2> type;
};
// Set the problem property