Commit 5520638e authored by Kilian Weishaupt's avatar Kilian Weishaupt
Browse files

[pnm][grid] Use PoreInscribedRadius and ThroatInscribedRadius consistently

parent 5ff479c5
......@@ -140,9 +140,9 @@ public:
{
const int eIdx = gridView.indexSet().index(element);
const auto& params = gridData.parameters(element);
static const auto throatRadiusIdx = gridData.parameterIndex("ThroatRadius");
static const auto throatInscribedRadiusIdx = gridData.parameterIndex("ThroatInscribedRadius");
static const auto throatLengthIdx = gridData.parameterIndex("ThroatLength");
throatInscribedRadius_[eIdx] = params[throatRadiusIdx];
throatInscribedRadius_[eIdx] = params[throatInscribedRadiusIdx];
throatLength_[eIdx] = params[throatLengthIdx];
// use a default value if no throat label is given by the grid
......
......@@ -247,7 +247,7 @@ public:
else if (paramName.find("Pore") != std::string::npos)
msg = "Make sure to include it in the vector of parameter names VertexParameters = " + paramName + " ... ...";
DUNE_THROW(Dumux::ParameterException, paramName << " not set in the grid or input file. \n" << msg << "\n" << list.str());
DUNE_THROW(Dumux::ParameterException, paramName << " not set in the grid file. \n" << msg << "\n" << list.str());
}
}
......@@ -317,7 +317,7 @@ private:
if (!isDgfData_)
{
vertexParameterNames_ = StringVector{"PoreInscribedRadius", "PoreVolume", "PoreLabel"};
elementParameterNames_ = StringVector{"ThroatRadius", "ThroatLength", "ThroatLabel"};
elementParameterNames_ = StringVector{"ThroatInscribedRadius", "ThroatLength", "ThroatLabel"};
if (numSubregions_ > 0)
{
vertexParameterNames_.push_back("PoreRegionId");
......
......@@ -194,15 +194,15 @@ public:
}
// treat throat parameters
auto defaultThroatRadius = throatRadiusGenerator_(-1, getParameter);
auto defaultThroatInscribedRadius = throatInscribedRadiusGenerator_(-1, getParameter);
auto defaultThroatLength = throatLengthGenerator_(-1, getParameter);
// get helper functions for throat radii and lengths on subregions
std::vector<decltype(defaultThroatRadius)> subregionThroatRadius;
std::vector<decltype(defaultThroatInscribedRadius)> subregionThroatInscribedRadius;
std::vector<decltype(defaultThroatLength)> subregionThroatLength;
for (int i = 0; i < numSubregions; ++i)
{
subregionThroatRadius.emplace_back(throatRadiusGenerator_(i, getParameter));
subregionThroatInscribedRadius.emplace_back(throatInscribedRadiusGenerator_(i, getParameter));
subregionThroatLength.emplace_back(throatLengthGenerator_(i, getParameter));
}
......@@ -211,14 +211,14 @@ public:
{
if (numSubregions == 0) // assign values if no subregions are specified
{
setParameter(element, "ThroatRadius", defaultThroatRadius(element));
setParameter(element, "ThroatInscribedRadius", defaultThroatInscribedRadius(element));
setParameter(element, "ThroatLength", defaultThroatLength(element));
}
else // assign values to throats if they are within a subregion
{
// default value for elements not belonging to a subregion
setParameter(element, "ThroatRegionId", -1);
setParameter(element, "ThroatRadius", defaultThroatRadius(element));
setParameter(element, "ThroatInscribedRadius", defaultThroatInscribedRadius(element));
setParameter(element, "ThroatLength", defaultThroatLength(element));
for (int id = 0; id < numSubregions; ++id)
......@@ -227,7 +227,7 @@ public:
if (intersectsPointGeometry(element.geometry().center(), subregion))
{
setParameter(element, "ThroatRegionId", id);
setParameter(element, "ThroatRadius", subregionThroatRadius[id](element));
setParameter(element, "ThroatInscribedRadius", subregionThroatInscribedRadius[id](element));
setParameter(element, "ThroatLength", subregionThroatLength[id](element));
}
}
......@@ -435,7 +435,7 @@ private:
};
};
const Scalar fixedPoreRadius = getParamFromGroup<Scalar>(paramGroup_, prefix + "PoreRadius", -1.0);
const Scalar fixedPoreRadius = getParamFromGroup<Scalar>(paramGroup_, prefix + "PoreInscribedRadius", -1.0);
// return random radius according to a user-specified distribution
if (fixedPoreRadius <= 0.0)
{
......@@ -443,12 +443,13 @@ private:
const auto seed = getParamFromGroup<unsigned int>(paramGroup_, prefix + "ParameterRandomNumberSeed", std::random_device{}());
generator.seed(seed);
const auto type = getParamFromGroup<std::string>(paramGroup_, prefix + "ParameterType");
const auto type = getParamFromGroup<std::string>(paramGroup_, prefix + "ParameterType", "lognormal");
if (type == "lognormal")
{
// if we use a lognormal distribution, get the mean and standard deviation from input file
const Scalar meanPoreRadius = getParamFromGroup<Scalar>(paramGroup_, prefix + "MeanPoreRadius");
const Scalar stddevPoreRadius = getParamFromGroup<Scalar>(paramGroup_, prefix + "StandardDeviationPoreRadius");
const auto [meanPoreRadius, stddevPoreRadius] = getDistributionInputParams_("Lognormal", prefix,
"MeanPoreInscribedRadius",
"StandardDeviationPoreInscribedRadius");
const Scalar variance = stddevPoreRadius*stddevPoreRadius;
using std::log;
......@@ -462,12 +463,11 @@ private:
else if (type == "uniform")
{
// if we use a uniform distribution, get the min and max from input file
const Scalar minPoreRadius = getParamFromGroup<Scalar>(paramGroup_, prefix + "MinPoreRadius");
const Scalar maxPoreRadius = getParamFromGroup<Scalar>(paramGroup_, prefix + "MaxPoreRadius");
const auto [minPoreRadius, maxPoreRadius] = getDistributionInputParams_("Uniform", prefix,
"MinPoreInscribedRadius",
"MaxPoreInscribedRadius");
Dumux::SimpleUniformDistribution<> poreRadiusDist(minPoreRadius, maxPoreRadius);
return generateFunction(poreRadiusDist);
}
else
DUNE_THROW(Dune::InvalidStateException, "Unknown parameter type " << type);
......@@ -481,6 +481,26 @@ private:
}
}
// print helpful error message if params are not properly provided
std::array<Scalar, 2> getDistributionInputParams_(const std::string& distributionName,
const std::string& prefix,
const std::string& paramName0,
const std::string& paramName1) const
{
try
{
return std::array{getParamFromGroup<Scalar>(paramGroup_, prefix + paramName0),
getParamFromGroup<Scalar>(paramGroup_, prefix + paramName1)};
}
catch (const Dumux::ParameterException& e)
{
std::cout << "\n" << distributionName << " pore-size distribution needs input parameters "
<< prefix + paramName0 << " and " << prefix + paramName1 << ".\n"
<< "Alternatively, use " << prefix << "PoreInscribedRadius to set a fixed inscribed pore radius." << std::endl;
DUNE_THROW(Dumux::ParameterException, e.what());
}
}
template <class GetParameter>
auto poreVolumeGenerator_(const GetParameter& getParameter) const
{
......@@ -531,16 +551,16 @@ private:
// returns a lambda taking a element and returning a radius
template <class GetParameter>
auto throatRadiusGenerator_(const int subregionId, const GetParameter& getParameter) const
auto throatInscribedRadiusGenerator_(const int subregionId, const GetParameter& getParameter) const
{
// adapt the parameter name if there are subregions
const std::string prefix = subregionId < 0 ? "Grid." : "Grid.Subregion" + std::to_string(subregionId) + ".";
// check for a user-specified fixed throat radius
const Scalar inputThroatRadius = getParamFromGroup<Scalar>(paramGroup_, prefix + "ThroatRadius", -1.0);
const Scalar inputThroatInscribedRadius = getParamFromGroup<Scalar>(paramGroup_, prefix + "ThroatInscribedRadius", -1.0);
// shape parameter for calculation of throat radius
const Scalar throatN = getParamFromGroup<Scalar>(paramGroup_, prefix + "ThroatRadiusN", 0.1);
const Scalar throatN = getParamFromGroup<Scalar>(paramGroup_, prefix + "ThroatInscribedRadiusN", 0.1);
return [=](const Element& element)
{
......@@ -548,8 +568,8 @@ private:
const std::array<Vertex, 2> vertices = {element.template subEntity<dim>(0), element.template subEntity<dim>(1)};
// the element parameters (throat radius and length)
if (inputThroatRadius > 0.0)
return inputThroatRadius;
if (inputThroatInscribedRadius > 0.0)
return inputThroatInscribedRadius;
else
{
const Scalar poreRadius0 = getParameter(vertices[0], "PoreInscribedRadius");
......
......@@ -10,8 +10,8 @@ UpperRight = 4e-4 4e-4 4e-4
NumPores = 20 20 20
ParameterType = lognormal
MeanPoreRadius = 1e-5
StandardDeviationPoreRadius = 4e-6
MeanPoreInscribedRadius = 1e-5
StandardDeviationPoreInscribedRadius = 4e-6
PoreGeometry = Sphere
ThroatCrossSectionShape = Square
ParameterRandomNumberSeed = 1
......
......@@ -5,9 +5,9 @@ Name = constant-params
UpperRight = 1e-3 1e-3 1e-3
NumPores = 5 5 5
PoreRadius = 1e-5
PoreInscribedRadius = 1e-5
PoreGeometry = Sphere
ThroatRadius = 1e-6
ThroatInscribedRadius = 1e-6
ThroatLength = 1e-4
BoundaryFaceMarker = 2 3 1 1 1 1
......@@ -28,9 +28,9 @@ Sanitize = true
UpperRight = 1e-3 1e-3
NumPores = 5 5
PoreRadius = 1e-5
PoreInscribedRadius = 1e-5
PoreGeometry = Cube
ThroatRadius = 1e-6
ThroatInscribedRadius = 1e-6
ThroatLength = 1e-4
BoundaryFaceMarker = 2 3 1 1
......@@ -50,10 +50,10 @@ Sanitize = true
UpperRight = 1e-3
NumPores = 20
PoreRadius = 1e-5
PoreInscribedRadius = 1e-5
PoreGeometry = Cylinder
PoreHeight = 1e-5
ThroatRadius = 1e-6
ThroatInscribedRadius = 1e-6
ThroatLength = 1e-4
BoundaryFaceMarker = 1 1
......
......@@ -6,8 +6,8 @@ UpperRight = 1e-3 1e-3 1e-3
NumPores = 5 5 5
ParameterType = lognormal
MeanPoreRadius = 1e-5
StandardDeviationPoreRadius = 3e-6
MeanPoreInscribedRadius = 1e-5
StandardDeviationPoreInscribedRadius = 3e-6
PoreGeometry = Sphere
BoundaryFaceMarker = 2 3 1 1 1 1
......@@ -29,8 +29,8 @@ UpperRight = 1e-3 1e-3
NumPores = 5 5
ParameterType = lognormal
MeanPoreRadius = 1e-5
StandardDeviationPoreRadius = 3e-6
MeanPoreInscribedRadius = 1e-5
StandardDeviationPoreInscribedRadius = 3e-6
PoreGeometry = Cube
BoundaryFaceMarker = 2 3 1 1
......@@ -52,8 +52,8 @@ UpperRight = 1e-3
NumPores = 20
ParameterType = lognormal
MeanPoreRadius = 1e-5
StandardDeviationPoreRadius = 3e-6
MeanPoreInscribedRadius = 1e-5
StandardDeviationPoreInscribedRadius = 3e-6
PoreGeometry = Cylinder
PoreHeight = 1e-5
......
......@@ -3,10 +3,10 @@ NumPores = 10 10 10
UpperRight = 1e-3 1e-3 1e-3
NumSubregions = 1
PoreRadius = 1e-5
PoreInscribedRadius = 1e-5
PoreGeometry = Sphere
ThroatCrossSectionShape = Circle
Subregion0.LowerLeft = 0 0 0
Subregion0.UpperRight = 0.5e-3 0.5e-3 0.5e-3
Subregion0.PoreRadius = 2e-5
Subregion0.PoreInscribedRadius = 2e-5
......@@ -6,7 +6,7 @@ WriteVTK = false
UpperRight = 1e-3 1e-3
NumPores = 5 5
PoreRadius = 1e-4
PoreInscribedRadius = 1e-4
ThroatRadius = 1e-5
PoreGeometry = Sphere
......
DGF
% Vertex parameters: PoreInscribedRadius PoreLabel
% Element parameters: ThroatRadius ThroatLength ThroatLabel
% Element parameters: ThroatInscribedRadius ThroatLength ThroatLabel
Vertex % Coordinates and volumes of the pore bodies
parameters 2
1 1 4 0.0002263 2.0
......
......@@ -4,8 +4,8 @@ ThroatCrossSectionShape = Square
UpperRight = 1e-3 1e-3 1e-3
NumPores = 4 4 4
PoreRadius = 2e-5
ThroatRadius = 2e-6
PoreInscribedRadius = 2e-5
ThroatInscribedRadius = 2e-6
BoundaryFaceMarker = 1 1 1 1 2 1
PriorityList = 4 5 0 1 2 3
......
......@@ -4,8 +4,8 @@ NumPores = 10 10 10
PoreGeometry = Cube
ThroatCrossSectionShape = Square
ParameterType = lognormal
MeanPoreRadius = 4.5e-5
StandardDeviationPoreRadius = 3e-6
MeanPoreInscribedRadius = 4.5e-5
StandardDeviationPoreInscribedRadius = 3e-6
BoundaryFaceMarker = 2 3 1 1 1 1
MinThroatLength = 1e-9
DeletionProbability = 0.9 0.9 0.9 0.9 0.9 0.9 0.9 0.9 0.9 0.9 0.9 0.9 0.9
......
......@@ -16,8 +16,8 @@
-1 -1 -1 -1 -1 -1 -1 1
</DataArray>
</PointData>
<CellData Scalars="throatRadius">
<DataArray type="Float32" Name="throatRadius" NumberOfComponents="1" format="ascii">
<CellData Scalars="throatInscribedRadius">
<DataArray type="Float32" Name="throatInscribedRadius" NumberOfComponents="1" format="ascii">
1e-06 1e-06 1e-06 1e-06 1e-06 1e-06 1e-06 1e-06 1e-06 1e-06 1e-06 1e-06
1e-06 1e-06 1e-06 1e-06 1e-06 1e-06 1e-06
</DataArray>
......
......@@ -16,8 +16,8 @@
-1 -1 -1 -1 -1 -1 -1 1
</DataArray>
</PointData>
<CellData Scalars="throatRadius">
<DataArray type="Float32" Name="throatRadius" NumberOfComponents="1" format="ascii">
<CellData Scalars="throatInscribedRadius">
<DataArray type="Float32" Name="throatInscribedRadius" NumberOfComponents="1" format="ascii">
4.14377e-06 4.74707e-06 4.78717e-06 5.06117e-06 4.85118e-06 5.15151e-06 5.76601e-06 4.33998e-06 4.34305e-06 5.23155e-06 5.23184e-06 9.47504e-06
7.18404e-06 6.61131e-06 6.932e-06 7.54417e-06 4.15076e-06 4.13528e-06 5.6777e-06
</DataArray>
......
......@@ -16,8 +16,8 @@
-1 -1 3 -1 -1 -1 3 1 2 1 1 3
</DataArray>
</PointData>
<CellData Scalars="throatRadius">
<DataArray type="Float32" Name="throatRadius" NumberOfComponents="1" format="ascii">
<CellData Scalars="throatInscribedRadius">
<DataArray type="Float32" Name="throatInscribedRadius" NumberOfComponents="1" format="ascii">
1e-06 1e-06 1e-06 1e-06 1e-06 1e-06 1e-06 1e-06 1e-06 1e-06 1e-06 1e-06
1e-06 1e-06 1e-06 1e-06 1e-06 1e-06 1e-06 1e-06 1e-06 1e-06 1e-06 1e-06
1e-06 1e-06 1e-06 1e-06 1e-06 1e-06 1e-06 1e-06 1e-06 1e-06 1e-06
......
......@@ -16,8 +16,8 @@
-1 -1 3 -1 -1 -1 3 1 2 1 1 3
</DataArray>
</PointData>
<CellData Scalars="throatRadius">
<DataArray type="Float32" Name="throatRadius" NumberOfComponents="1" format="ascii">
<CellData Scalars="throatInscribedRadius">
<DataArray type="Float32" Name="throatInscribedRadius" NumberOfComponents="1" format="ascii">
4.11658e-06 4.70756e-06 4.12337e-06 4.81161e-06 4.09397e-06 5.71113e-06 5.03676e-06 5.10376e-06 4.28953e-06 6.9563e-06 5.18869e-06 6.48355e-06
9.31931e-06 5.8116e-06 6.58115e-06 7.10058e-06 5.25119e-06 6.98718e-06 6.5364e-06 5.18384e-06 8.64324e-06 6.8266e-06 4.12387e-06 7.58904e-06
4.12097e-06 5.64464e-06 5.61358e-06 7.12156e-06 7.21317e-06 4.12475e-06 8.66678e-06 4.12122e-06 7.35925e-06 4.12238e-06 7.42437e-06
......
......@@ -25,8 +25,8 @@
1 1 2 1 3 1 1 2 3 1 1
</DataArray>
</PointData>
<CellData Scalars="throatRadius">
<DataArray type="Float32" Name="throatRadius" NumberOfComponents="1" format="ascii">
<CellData Scalars="throatInscribedRadius">
<DataArray type="Float32" Name="throatInscribedRadius" NumberOfComponents="1" format="ascii">
1e-06 1e-06 1e-06 1e-06 1e-06 1e-06 1e-06 1e-06 1e-06 1e-06 1e-06 1e-06
1e-06 1e-06 1e-06 1e-06 1e-06 1e-06 1e-06 1e-06 1e-06 1e-06 1e-06 1e-06
1e-06 1e-06 1e-06 1e-06 1e-06 1e-06 1e-06 1e-06 1e-06 1e-06 1e-06 1e-06
......
......@@ -25,8 +25,8 @@
1 1 2 1 3 1 1 2 3 1 1
</DataArray>
</PointData>
<CellData Scalars="throatRadius">
<DataArray type="Float32" Name="throatRadius" NumberOfComponents="1" format="ascii">
<CellData Scalars="throatInscribedRadius">
<DataArray type="Float32" Name="throatInscribedRadius" NumberOfComponents="1" format="ascii">
4.11658e-06 4.70756e-06 4.31042e-06 6.9563e-06 5.16132e-06 5.18517e-06 4.31144e-06 5.7185e-06 4.31344e-06 7.44744e-06 4.12163e-06 5.61949e-06
7.21941e-06 7.20563e-06 8.31374e-06 7.59153e-06 7.33705e-06 6.99583e-06 5.75207e-06 7.54291e-06 4.31091e-06 4.10902e-06 5.81255e-06 8.10476e-06
8.10679e-06 5.18072e-06 1.22845e-05 4.10978e-06 5.08512e-06 7.11641e-06 7.20984e-06 7.21875e-06 6.25225e-06 4.9708e-06 6.30755e-06 6.10565e-06
......
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment