Commit 2b68068f authored by Timo Koch's avatar Timo Koch
Browse files

Merge branch 'fix/velocity-output' into 'master'

Fix/velocity output

Closes #1161

See merge request !3134
parents 307b13d4 6c2c9fda
Pipeline #17362 passed with stages
in 0 seconds
......@@ -34,6 +34,7 @@
#include <dune/geometry/referenceelements.hh>
#include <dumux/common/parameters.hh>
#include <dumux/discretization/extrusion.hh>
#include <dumux/discretization/method.hh>
#include <dumux/discretization/elementsolution.hh>
#include <dumux/flux/traits.hh>
......@@ -69,6 +70,7 @@ template<class GridVariables, class FluxVariables>
class PorousMediumFlowVelocity
{
using GridGeometry = typename GridVariables::GridGeometry;
using Extrusion = Extrusion_t<GridGeometry>;
using FVElementGeometry = typename GridGeometry::LocalView;
using SubControlVolume = typename GridGeometry::SubControlVolume;
using SubControlVolumeFace = typename GridGeometry::SubControlVolumeFace;
......@@ -161,7 +163,7 @@ public:
Scalar localArea = scvfReferenceArea_(geomType, scvf.index());
Scalar flux = fluxVars.advectiveFlux(phaseIdx, upwindTerm) / localArea;
const auto& insideVolVars = elemVolVars[scvf.insideScvIdx()];
flux /= insideVolVars.extrusionFactor();
flux /= insideVolVars.extrusionFactor() * Extrusion::area(scvf) / scvf.area();
tmpVelocity *= flux;
const int eIdxGlobal = gridGeometry_.elementMapper().index(element);
......@@ -206,7 +208,7 @@ public:
Scalar localArea = scvfReferenceArea_(geomType, scvf.index());
Scalar flux = fluxVars.advectiveFlux(phaseIdx, upwindTerm) / localArea;
const auto& insideVolVars = elemVolVars[scvf.insideScvIdx()];
flux /= insideVolVars.extrusionFactor();
flux /= insideVolVars.extrusionFactor() * Extrusion::area(scvf) / scvf.area();
// transform the volume flux into a velocity vector
Velocity tmpVelocity = localNormal;
......
......@@ -30,6 +30,7 @@ and compute the convergence rates.
#include <config.h>
#include <iostream>
#include <vector>
#include <dumux/common/initialize.hh>
#include <dumux/common/properties.hh> // for GetPropType
......@@ -90,22 +91,30 @@ int main(int argc, char** argv) try
// We define a function to update the discrete analytical solution vector
// using the exactSolution() function in the problem
const auto updateAnalyticalSolution = [&](auto& pExact)
const auto updateAnalyticalSolution = [&](auto& pExact, auto& vExact)
{
pExact.resize(gridGeometry->numDofs());
vExact.resize(gridGeometry->elementMapper().size());
for (const auto& element : elements(gridGeometry->gridView()))
{
auto fvGeometry = localView(*gridGeometry);
fvGeometry.bindElement(element);
for (auto&& scv : scvs(fvGeometry))
{
pExact[scv.dofIndex()] = problem->exactSolution(scv.dofPosition());
}
const auto eIdx = gridGeometry->elementMapper().index(element);
vExact[eIdx] = problem->exactVelocity(element.geometry().center());
}
};
// instantiate and initialize the discrete and exact solution vectors
using SolutionVector = GetPropType<TypeTag, Properties::SolutionVector>;
SolutionVector p(gridGeometry->numDofs());
SolutionVector pExact; updateAnalyticalSolution(pExact);
SolutionVector pExact;
std::vector<double> vExact;
updateAnalyticalSolution(pExact, vExact);
// instantiate and initialize the grid variables
using GridVariables = GetPropType<TypeTag, Properties::GridVariables>;
......@@ -119,6 +128,10 @@ int main(int argc, char** argv) try
VtkOutputModule<GridVariables, SolutionVector> vtkWriter(*gridVariables, p, problem->name());
GetPropType<TypeTag, Properties::IOFields>::initOutputModule(vtkWriter);
vtkWriter.addField(pExact, "pExact"); // add the exact solution to the output fields
vtkWriter.addField(vExact, "vExact"); // add the exact velocity to the output fields
using VelocityOutput = GetPropType<TypeTag, Properties::VelocityOutput>;
vtkWriter.addVelocityOutput(std::make_shared<VelocityOutput>(*gridVariables));
```
### Instantiate the solver
......@@ -179,7 +192,7 @@ in the input file.
assembler->updateAfterGridAdaption();
p.resize(gridGeometry->numDofs());
updateAnalyticalSolution(pExact);
updateAnalyticalSolution(pExact, vExact);
// solve problem on refined grid
solver.solve(p);
......
......@@ -213,6 +213,7 @@ properties and boundary conditions, from the input file.
// fluid properties
k_ = getParam<Scalar>("SpatialParams.Permeability");
nu_ = getParam<Scalar>("Component.LiquidKinematicViscosity");
rho_ = getParam<Scalar>("Component.LiquidDensity");
// The inner radius r1 can be determined from the grid
r1_ = gridGeometry->bBoxMin()[0];
......@@ -220,7 +221,6 @@ properties and boundary conditions, from the input file.
// boundary conditions
q1_ = getParam<Scalar>("Problem.Q1"); // mass flux into the domain at r1 in kg/s/m
p1_ = getParam<Scalar>("Problem.P1"); // pressure at the inner boundary at r1
}
```
......@@ -271,9 +271,16 @@ different levels of grid refinement.
return p;
}
const Scalar exactVelocity(const GlobalPosition& globalPos) const
{
const auto r = globalPos[0];
const auto v = q1_/(2*M_PI)/rho_/r;
return v;
}
private:
// private data members required for the analytical solution
Scalar q1_, k_, nu_, r1_, p1_;
Scalar q1_, k_, nu_, r1_, p1_, rho_;
};
} // end namespace Dumux
......
......@@ -27,6 +27,7 @@
#include <config.h>
#include <iostream>
#include <vector>
#include <dumux/common/initialize.hh>
#include <dumux/common/properties.hh> // for GetPropType
......@@ -83,22 +84,30 @@ int main(int argc, char** argv) try
// We define a function to update the discrete analytical solution vector
// using the exactSolution() function in the problem
const auto updateAnalyticalSolution = [&](auto& pExact)
const auto updateAnalyticalSolution = [&](auto& pExact, auto& vExact)
{
pExact.resize(gridGeometry->numDofs());
vExact.resize(gridGeometry->elementMapper().size());
for (const auto& element : elements(gridGeometry->gridView()))
{
auto fvGeometry = localView(*gridGeometry);
fvGeometry.bindElement(element);
for (auto&& scv : scvs(fvGeometry))
{
pExact[scv.dofIndex()] = problem->exactSolution(scv.dofPosition());
}
const auto eIdx = gridGeometry->elementMapper().index(element);
vExact[eIdx] = problem->exactVelocity(element.geometry().center());
}
};
// instantiate and initialize the discrete and exact solution vectors
using SolutionVector = GetPropType<TypeTag, Properties::SolutionVector>;
SolutionVector p(gridGeometry->numDofs());
SolutionVector pExact; updateAnalyticalSolution(pExact);
SolutionVector pExact;
std::vector<double> vExact;
updateAnalyticalSolution(pExact, vExact);
// instantiate and initialize the grid variables
using GridVariables = GetPropType<TypeTag, Properties::GridVariables>;
......@@ -110,6 +119,10 @@ int main(int argc, char** argv) try
VtkOutputModule<GridVariables, SolutionVector> vtkWriter(*gridVariables, p, problem->name());
GetPropType<TypeTag, Properties::IOFields>::initOutputModule(vtkWriter);
vtkWriter.addField(pExact, "pExact"); // add the exact solution to the output fields
vtkWriter.addField(vExact, "vExact"); // add the exact velocity to the output fields
using VelocityOutput = GetPropType<TypeTag, Properties::VelocityOutput>;
vtkWriter.addVelocityOutput(std::make_shared<VelocityOutput>(*gridVariables));
// ### Instantiate the solver
// We use the `LinearPDESolver` class, which is instantiated on the basis
......@@ -166,7 +179,7 @@ int main(int argc, char** argv) try
assembler->updateAfterGridAdaption();
p.resize(gridGeometry->numDofs());
updateAnalyticalSolution(pExact);
updateAnalyticalSolution(pExact, vExact);
// solve problem on refined grid
solver.solve(p);
......
......@@ -16,3 +16,6 @@ Permeability = 1e-10 # [m^2]
[Component]
LiquidKinematicViscosity = 1e-6
LiquidDensity = 1e3
[Vtk]
AddVelocity = true
......@@ -57,6 +57,7 @@ public:
// fluid properties
k_ = getParam<Scalar>("SpatialParams.Permeability");
nu_ = getParam<Scalar>("Component.LiquidKinematicViscosity");
rho_ = getParam<Scalar>("Component.LiquidDensity");
// The inner radius r1 can be determined from the grid
r1_ = gridGeometry->bBoxMin()[0];
......@@ -64,7 +65,6 @@ public:
// boundary conditions
q1_ = getParam<Scalar>("Problem.Q1"); // mass flux into the domain at r1 in kg/s/m
p1_ = getParam<Scalar>("Problem.P1"); // pressure at the inner boundary at r1
}
// [[/codeblock]]
......@@ -108,9 +108,16 @@ public:
return p;
}
const Scalar exactVelocity(const GlobalPosition& globalPos) const
{
const auto r = globalPos[0];
const auto v = q1_/(2*M_PI)/rho_/r;
return v;
}
private:
// private data members required for the analytical solution
Scalar q1_, k_, nu_, r1_, p1_;
Scalar q1_, k_, nu_, r1_, p1_, rho_;
};
} // end namespace Dumux
......
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