Skip to content
Snippets Groups Projects
Commit d7bda623 authored by Dennis Gläser's avatar Dennis Gläser Committed by Timo Koch
Browse files

[1p] port ni tests to new interface

The tests fail, but they do so because the output has changed. This has
to be checked again.[1p] port ni tests to new interface

The tests fail, but they do so because the output has changed. This has
to be checked again.[1p] port ni tests to new interface

The tests fail, but they do so because the output has changed. This has
to be checked again.[1p] port ni tests to new interface

The tests fail, but they do so because the output has changed. This has
to be checked again.[1p] port ni tests to new interface

The tests fail, but they do so because the output has changed. This has
to be checked again.[1p] port ni tests to new interface

The tests fail, but they do so because the output has changed. This has
to be checked again.[1p] port ni tests to new interface

The tests fail, but they do so because the output has changed. This has
to be checked again.[1p] port ni tests to new interface

The tests fail, but they do so because the output has changed. This has
to be checked again.
parent 2995ce66
No related branches found
No related tags found
Loading
......@@ -141,7 +141,7 @@ private:
const auto& insideVolVars = elemVolVars[insideScvIdx];
auto insideLambda = ThermalConductivityModel::effectiveThermalConductivity(insideVolVars, problem.spatialParams(), element, fvGeometry, insideScv);
Scalar ti = calculateOmega_(problem, element, scvf, insideLambda, insideScv);
Scalar ti = calculateOmega_(scvf, insideLambda, insideScv, insideVolVars.extrusionFactor());
// for the boundary (dirichlet) or at branching points we only need ti
if (scvf.boundary() || scvf.numOutsideScvs() > 1)
......@@ -164,9 +164,9 @@ private:
Scalar tj;
if (dim == dimWorld)
// assume the normal vector from outside is anti parallel so we save flipping a vector
tj = -1.0*calculateOmega_(problem, outsideElement, scvf, outsideLambda, outsideScv);
tj = -1.0*calculateOmega_(scvf, outsideLambda, outsideScv, outsideVolVars.extrusionFactor());
else
tj = calculateOmega_(problem, outsideElement, fvGeometry.flipScvf(scvf.index()), outsideLambda, outsideScv);
tj = calculateOmega_(fvGeometry.flipScvf(scvf.index()), outsideLambda, outsideScv, outsideVolVars.extrusionFactor());
// check for division by zero!
if (ti*tj <= 0.0)
......@@ -178,11 +178,10 @@ private:
return tij;
}
static Scalar calculateOmega_(const Problem& problem,
const Element& element,
const SubControlVolumeFace& scvf,
const DimWorldMatrix &lambda,
const SubControlVolume &scv)
static Scalar calculateOmega_(const SubControlVolumeFace& scvf,
const DimWorldMatrix& lambda,
const SubControlVolume& scv,
Scalar extrusionFactor)
{
GlobalPosition lambdaNormal;
lambda.mv(scvf.unitOuterNormal(), lambdaNormal);
......@@ -192,25 +191,20 @@ private:
distanceVector /= distanceVector.two_norm2();
Scalar omega = lambdaNormal * distanceVector;
omega *= problem.boxExtrusionFactor(element, scv);
return omega;
return omega*extrusionFactor;
}
static Scalar calculateOmega_(const Problem& problem,
const Element& element,
const SubControlVolumeFace& scvf,
static Scalar calculateOmega_(const SubControlVolumeFace& scvf,
Scalar lambda,
const SubControlVolume &scv)
const SubControlVolume &scv,
Scalar extrusionFactor)
{
auto distanceVector = scvf.ipGlobal();
distanceVector -= scv.center();
distanceVector /= distanceVector.two_norm2();
Scalar omega = lambda * (distanceVector * scvf.unitOuterNormal());
omega *= problem.boxExtrusionFactor(element, scv);
return omega;
return omega*extrusionFactor;
}
};
......
......@@ -103,6 +103,7 @@ class OnePNIConductionProblem : public ImplicitPorousMediaProblem<TypeTag>
using TimeManager = typename GET_PROP_TYPE(TypeTag, TimeManager);
using ThermalConductivityModel = typename GET_PROP_TYPE(TypeTag, ThermalConductivityModel);
using ElementVolumeVariables = typename GET_PROP_TYPE(TypeTag, ElementVolumeVariables);
using ElementSolution = typename GET_PROP_TYPE(TypeTag, ElementSolutionVector);
using VolumeVariables = typename GET_PROP_TYPE(TypeTag, VolumeVariables);
// copy some indices for convenience
......@@ -168,20 +169,20 @@ public:
auto& temperature = *(this->resultWriter().allocateManagedBuffer(numDofs));
const auto someElement = *(elements(this->gridView()).begin());
const auto initialPriVars = initial_(GlobalPosition(0.0));
const auto someElemSol = this->model().elementSolution(someElement, this->model().curSol());
auto someFvGeometry = localView(this->model().globalFvGeometry());
someFvGeometry.bindElement(someElement);
const auto someScv = *(scvs(someFvGeometry).begin());
VolumeVariables volVars;
volVars.update(initialPriVars, *this, someElement, someScv);
volVars.update(someElemSol, *this, someElement, someScv);
const auto porosity = this->spatialParams().porosity(someScv);
const auto porosity = this->spatialParams().porosity(someElement, someScv, someElemSol);
const auto densityW = volVars.density();
const auto heatCapacityW = FluidSystem::heatCapacity(volVars.fluidState(), 0);
const auto densityS = this->spatialParams().solidDensity(someElement, someScv);
const auto heatCapacityS = this->spatialParams().solidHeatCapacity(someElement, someScv);
const auto densityS = this->spatialParams().solidDensity(someElement, someScv, someElemSol);
const auto heatCapacityS = this->spatialParams().solidHeatCapacity(someElement, someScv, someElemSol);
const auto storage = densityW*heatCapacityW*porosity + densityS*heatCapacityS*(1 - porosity);
const auto effectiveThermalConductivity = ThermalConductivityModel::effectiveThermalConductivity(volVars, this->spatialParams(),
someElement, someFvGeometry, someScv);
......@@ -197,7 +198,7 @@ public:
auto globalIdx = scv.dofIndex();
const auto& globalPos = scv.dofPosition();
temperatureExact[globalIdx] = temperatureHigh_ + (initialPriVars[temperatureIdx] - temperatureHigh_)
temperatureExact[globalIdx] = temperatureHigh_ + (someElemSol[0][temperatureIdx] - temperatureHigh_)
*std::erf(0.5*std::sqrt(globalPos[0]*globalPos[0]*storage/time/effectiveThermalConductivity));
temperature[globalIdx] = this->model().curSol()[globalIdx][temperatureIdx];
}
......
......@@ -176,21 +176,21 @@ public:
auto& temperature = *(this->resultWriter().allocateManagedBuffer(numDofs));
const auto someElement = *(elements(this->gridView()).begin());
const auto initialPriVars = initial_(GlobalPosition(0.0));
const auto someElemSol = this->model().elementSolution(someElement, this->model().curSol());
auto someFvGeometry = localView(this->model().globalFvGeometry());
someFvGeometry.bindElement(someElement);
const auto someScv = *(scvs(someFvGeometry).begin());
VolumeVariables volVars;
volVars.update(initialPriVars, *this, someElement, someScv);
volVars.update(someElemSol, *this, someElement, someScv);
const auto porosity = this->spatialParams().porosity(someScv);
const auto porosity = this->spatialParams().porosity(someElement, someScv, someElemSol);
const auto densityW = volVars.density();
const auto heatCapacityW = FluidSystem::heatCapacity(volVars.fluidState(), 0);
const auto storageW = densityW*heatCapacityW*porosity;
const auto densityS = this->spatialParams().solidDensity(someElement, someScv);
const auto heatCapacityS = this->spatialParams().solidHeatCapacity(someElement, someScv);
const auto densityS = this->spatialParams().solidDensity(someElement, someScv, someElemSol);
const auto heatCapacityS = this->spatialParams().solidHeatCapacity(someElement, someScv, someElemSol);
const auto storageTotal = storageW + densityS*heatCapacityS*(1 - porosity);
std::cout << "storage: " << storageTotal << std::endl;
......
......@@ -49,96 +49,56 @@ class OnePNISpatialParams : public ImplicitSpatialParamsOneP<TypeTag>
using VolumeVariables = typename GET_PROP_TYPE(TypeTag, VolumeVariables);
using Element = typename GridView::template Codim<0>::Entity;
static const int dimWorld = GridView::dimensionworld;
using GlobalPosition = Dune::FieldVector<Scalar, dimWorld>;
public:
OnePNISpatialParams(const Problem& problem, const GridView &gridView)
: ParentType(problem, gridView)
{
permeability_ = 1e-10;
porosity_ = 0.4;
// heat conductivity of granite
lambdaSolid_ = 2.8;
}
: ParentType(problem, gridView) {}
/*!
* \brief Define the intrinsic permeability \f$\mathrm{[m^2]}\f$.
*
* \param element The finite element
* \param scv The sub control volume
* \param globalPos The global position
*/
Scalar intrinsicPermeability(const SubControlVolume& scv,
const VolumeVariables& volVars) const
{
return permeability_;
}
Scalar permeabilityAtPos(const GlobalPosition& globalPos) const
{ return 1e-10; }
/*!
* \brief Define the porosity \f$\mathrm{[-]}\f$.
*
* \param element The finite element
* \param scv The sub control volume
* \param globalPos The global position
*/
Scalar porosity(const SubControlVolume& scv) const
{
return porosity_;
}
/*!
* \brief Define the dispersivity.
*
* \param element The finite element
* \param scv The sub control volume
*/
Scalar dispersivity(const Element &element,
const SubControlVolume& scv) const
{
return 0;
}
Scalar porosityAtPos(const GlobalPosition& globalPos) const
{ return 0.4; }
/*!
* \brief Returns the heat capacity \f$[J / (kg K)]\f$ of the rock matrix.
*
* This is only required for non-isothermal models.
*
* \param element The finite element
* \param scv The sub control volume
* \param globalPos The global position
*/
Scalar solidHeatCapacity(const Element &element,
const SubControlVolume& scv) const
{
return 790; // specific heat capacity of granite [J / (kg K)]
}
Scalar solidHeatCapacityAtPos(const GlobalPosition& globalPos) const
{ return 790; /*specific heat capacity of granite [J / (kg K)]*/ }
/*!
* \brief Returns the mass density \f$[kg / m^3]\f$ of the rock matrix.
*
* This is only required for non-isothermal models.
*
* \param element The finite element
* \param scv The sub control volume
* \param globalPos The global position
*/
Scalar solidDensity(const Element &element,
const SubControlVolume& scv) const
{
return 2700; // density of granite [kg/m^3]
}
Scalar solidDensityAtPos(const GlobalPosition& globalPos) const
{ return 2700; /*density of granite [kg/m^3]*/ }
/*!
* \brief Returns the thermal conductivity \f$\mathrm{[W/(m K)]}\f$ of the porous material.
*
* \param element The finite element
* \param scv The sub control volume
* \param globalPos The global position
*/
Scalar solidThermalConductivity(const Element &element,
const SubControlVolume& scv) const
{
return lambdaSolid_;
}
private:
Scalar permeability_;
Scalar porosity_;
Scalar lambdaSolid_;
Scalar solidThermalConductivityAtPos(const GlobalPosition& globalPos) const
{ return 2.8; }
};
} // end namespace Dumux
......
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment