Commit fba977d6 authored by Timo Koch's avatar Timo Koch
Browse files

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

Fix/tracer velocity

Closes #782

See merge request !1786
parents 08ade75c 35a86158
......@@ -18,5 +18,6 @@ maxwellstefanslaw.hh
referencesystemformulation.hh
shallowwaterflux.hh
stationaryvelocityfield.hh
traits.hh
upwindscheme.hh
DESTINATION ${CMAKE_INSTALL_INCLUDEDIR}/dumux/flux)
......@@ -27,6 +27,9 @@
#ifndef DUMUX_DISCRETIZATION_STATIONARY_VELOCITY_FIELD_HH
#define DUMUX_DISCRETIZATION_STATIONARY_VELOCITY_FIELD_HH
#include <type_traits>
#include <dumux/flux/traits.hh>
#include <dumux/discretization/method.hh>
#include <dumux/flux/fluxvariablescaching.hh>
......@@ -64,6 +67,11 @@ public:
}
};
//! Set stationary velocity field to true in the FluxTraits
template<class Scalar>
struct HasStationaryVelocityField<StationaryVelocityField<Scalar>>
: public std::true_type {};
} // end namespace Dumux
#endif
// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
// vi: set et ts=4 sw=4 sts=4:
/*****************************************************************************
* See the file COPYING for full copying permissions. *
* *
* This program is free software: you can redistribute it and/or modify *
* it under the terms of the GNU General Public License as published by *
* the Free Software Foundation, either version 3 of the License, or *
* (at your option) any later version. *
* *
* This program is distributed in the hope that it will be useful, *
* but WITHOUT ANY WARRANTY; without even the implied warranty of *
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
* GNU General Public License for more details. *
* *
* You should have received a copy of the GNU General Public License *
* along with this program. If not, see <http://www.gnu.org/licenses/>. *
*****************************************************************************/
/*!
* \file
* \ingroup Flux
* \brief Defines the flux traits.
*/
#ifndef DUMUX_FLUX_TRAITS_HH
#define DUMUX_FLUX_TRAITS_HH
#include <type_traits>
namespace Dumux {
/*!
* \ingroup Flux
* \brief Trait of an advection type stating whether it implements a stationary velocity field
*/
template<class AdvectionType>
struct HasStationaryVelocityField : public std::false_type {};
/*!
* \ingroup Flux
* \brief Traits of a flux variables type
*/
template<class FluxVariables>
struct FluxTraits
{
static constexpr bool hasStationaryVelocityField()
{
return HasStationaryVelocityField<typename FluxVariables::AdvectionType>::value;
}
};
} // namespace Dumux
#endif // DUMUX_FLUX_TRAITS_HH
......@@ -26,6 +26,7 @@
#define DUMUX_POROUSMEDIUMFLOW_VELOCITY_HH
#include <vector>
#include <dumux/flux/traits.hh>
#include <dune/common/fvector.hh>
#include <dune/common/float_cmp.hh>
......@@ -57,11 +58,14 @@ class PorousMediumFlowVelocity
using ElementVolumeVariables = typename GridVolumeVariables::LocalView;
using FluidSystem = typename VolumeVariables::FluidSystem;
using Scalar = typename GridVariables::Scalar;
using FluxTraits = typename Dumux::FluxTraits<FluxVariables>;
using AdvectionType = typename FluxVariables::AdvectionType;
static constexpr int dim = GridView::dimension;
static constexpr int dimWorld = GridView::dimensionworld;
static constexpr bool isBox = GridGeometry::discMethod == DiscretizationMethod::box;
static constexpr int dofCodim = isBox ? dim : 0;
static constexpr bool stationaryVelocityField = FluxTraits::hasStationaryVelocityField();
using GlobalPosition = typename Element::Geometry::GlobalCoordinate;
using ReferenceElements = Dune::ReferenceElements<typename GridView::ctype, dim>;
......@@ -282,21 +286,42 @@ public:
auto bcTypes = problemBoundaryTypes_(element, scvf);
if (bcTypes.hasNeumann())
{
// check if we have Neumann no flow, we can just use 0
const auto neumannFlux = Deprecated::neumann(problem_, element, fvGeometry, elemVolVars, elemFluxVarsCache, scvf);
using NumEqVector = std::decay_t<decltype(neumannFlux)>;
if (Dune::FloatCmp::eq<NumEqVector, Dune::FloatCmp::CmpStyle::absolute>(neumannFlux, NumEqVector(0.0), 1e-30))
scvfFluxes[scvfIndexInInside[localScvfIdx]] = 0;
// cubes
else if (dim == 1 || geomType.isCube())
// for stationary velocity fields we can easily compute the correct velocity
// this special treatment makes sure that the velocity field is also correct on Neumann boundaries
// of tracer problems where the velocity field is given.
// (For Dirichlet boundaries no special treatment is necessary.)
if (stationaryVelocityField)
{
const auto fIdx = scvfIndexInInside[localScvfIdx];
const auto fIdxOpposite = fIdx%2 ? fIdx-1 : fIdx+1;
scvfFluxes[fIdx] = -scvfFluxes[fIdxOpposite];
const auto flux = AdvectionType::flux(problem_, element, fvGeometry, elemVolVars,
scvf, phaseIdx, elemFluxVarsCache);
const auto& insideVolVars = elemVolVars[scvf.insideScvIdx()];
scvfFluxes[scvfIndexInInside[localScvfIdx]] += flux / insideVolVars.extrusionFactor();
}
else
{
// check if we have Neumann no flow, we can just use 0
const auto neumannFlux = Deprecated::neumann(problem_, element, fvGeometry, elemVolVars, elemFluxVarsCache, scvf);
using NumEqVector = std::decay_t<decltype(neumannFlux)>;
if (Dune::FloatCmp::eq<NumEqVector, Dune::FloatCmp::CmpStyle::absolute>(neumannFlux, NumEqVector(0.0), 1e-30))
scvfFluxes[scvfIndexInInside[localScvfIdx]] = 0;
// otherwise, we try some reconstruction (TODO: Can this be improved?)
// for cubes
else if (dim == 1 || geomType.isCube())
{
const auto fIdx = scvfIndexInInside[localScvfIdx];
const auto fIdxOpposite = fIdx%2 ? fIdx-1 : fIdx+1;
scvfFluxes[fIdx] = -scvfFluxes[fIdxOpposite];
}
// for simplices
else if (geomType.isSimplex())
scvfFluxes[scvfIndexInInside[localScvfIdx]] = 0;
else
DUNE_THROW(Dune::NotImplemented, "Velocity computation at Neumann boundaries for cell-centered and prism/pyramid");
}
// simplices
else if (geomType.isSimplex())
scvfFluxes[scvfIndexInInside[localScvfIdx]] = 0;
}
}
......@@ -304,6 +329,7 @@ public:
localScvfIdx++;
}
Velocity refVelocity;
// cubes: On the reference element simply average over opposite fluxes
// note that this is equal to a corner velocity interpolation method
......@@ -324,7 +350,7 @@ public:
}
// 3D prism and pyramids
else
DUNE_THROW(Dune::NotImplemented, "velocity output for cell-centered and prism/pyramid");
DUNE_THROW(Dune::NotImplemented, "Velocity computation for cell-centered and prism/pyramid");
Velocity scvVelocity(0);
jacobianT2.mtv(refVelocity, scvVelocity);
......
......@@ -637,19 +637,19 @@
1000 1000 1000 1000
</DataArray>
<DataArray type="Float32" Name="velocity_Groundwater (m/s)" NumberOfComponents="3" format="ascii">
-1.102022907e-06 2.590747044e-06 0 -1.472232611e-06 7.210132935e-07 0 6.135285275e-07 1.226460569e-07 0 -2.206547833e-06 7.605763585e-06 0
-3.04372179e-06 1.027269263e-06 0 7.624727516e-08 3.685960905e-07 0 3.742941544e-06 1.319326543e-06 0 9.045208458e-07 1.679732486e-05 0
-7.539802482e-06 6.671227311e-06 0 -3.873328296e-06 2.305055204e-06 0 8.939652503e-07 1.558288773e-06 0 2.047548548e-07 2.416057612e-07 0
-1.52349422e-08 6.436293916e-07 0 -1.746525129e-07 1.23548607e-05 0 -2.983554452e-07 3.483635567e-07 0 1.603786473e-07 8.051803889e-06 0
3.839931537e-07 8.730427794e-07 0 1.358515611e-08 1.950072829e-06 0 1.661969229e-08 3.158886841e-08 0 2.64002491e-08 2.536741249e-06 0
7.600337248e-09 3.272032245e-08 0 3.134172388e-08 1.546748877e-06 0 4.832586598e-08 7.22737866e-07 0 3.008161542e-08 9.104686427e-08 0
-7.760051268e-09 1.527252209e-07 0 -1.162628791e-06 3.392388862e-06 0 -1.843889891e-06 6.451628565e-07 0 -5.605528486e-07 1.533103841e-06 0
1.350256014e-07 3.578106202e-07 0 2.529583298e-08 4.790475941e-06 0 -2.5670181e-07 2.579756256e-06 0 9.307453475e-09 6.960911492e-07 0
8.341055491e-08 6.241996289e-06 0 -8.93788183e-07 2.25847225e-06 0 -7.029688618e-07 5.940153969e-08 0 2.833088047e-06 1.293294758e-06 0
3.096192358e-06 5.905750186e-06 0 2.584979768e-07 7.980457326e-07 0 -4.573736021e-07 1.249281581e-06 0 -7.240656146e-07 3.706823861e-07 0
3.782706415e-07 5.187374427e-06 0 7.419180292e-07 1.346798854e-06 0 -1.184608678e-07 6.10841937e-07 0 -1.537402568e-06 3.708992381e-06 0
-7.679751093e-07 6.854446838e-06 0 3.425169837e-07 1.287100986e-06 0 -3.235867894e-07 2.512623269e-06 0 -1.030137895e-07 7.294159786e-06 0
2.185198582e-06 2.858712378e-07 0 2.157496056e-06 7.005677617e-06 0 -3.726464683e-07 5.554140444e-06 0 -4.169659462e-07 1.113699568e-06 0
-1.10202e-06 4.07947e-06 0 -1.47223e-06 2.17384e-06 0 6.13529e-07 1.59924e-06 0 -2.20655e-06 1.10375e-05 0
-3.04372e-06 5.39139e-06 0 7.62473e-08 5.20311e-07 0 3.74294e-06 6.52223e-06 0 9.04521e-07 2.68727e-05 0
-7.5398e-06 1.16201e-05 0 -3.87333e-06 9.99891e-06 0 8.93965e-07 2.49507e-06 0 2.04755e-07 4.15509e-07 0
-1.52349e-08 1.13497e-06 0 -1.74653e-07 2.47026e-05 0 -2.98355e-07 5.80154e-07 0 1.60379e-07 1.66789e-05 0
3.83993e-07 1.39439e-06 0 1.35852e-08 3.88143e-06 0 1.66197e-08 8.49279e-08 0 2.64002e-08 5.06151e-06 0
7.60034e-09 5.86103e-08 0 3.13417e-08 3.12407e-06 0 4.83259e-08 1.43189e-06 0 3.00816e-08 1.77437e-07 0
-7.76005e-09 2.72265e-07 0 -1.16263e-06 5.66309e-06 0 -1.84389e-06 1.73075e-06 0 -5.60553e-07 3.90912e-06 0
1.35026e-07 5.68285e-07 0 2.52958e-08 9.61856e-06 0 -2.56702e-07 4.83991e-06 0 9.30745e-09 1.9778e-06 0
8.34106e-08 1.19725e-05 0 -8.93788e-07 4.05126e-06 0 -7.02969e-07 7.75311e-07 0 2.83309e-06 5.46614e-06 0
3.09619e-06 9.19506e-06 0 2.58498e-07 1.37484e-06 0 -4.57374e-07 2.00394e-06 0 -7.24066e-07 9.69294e-07 0
3.78271e-07 1.12492e-05 0 7.41918e-07 2.18284e-06 0 -1.18461e-07 8.72064e-07 0 -1.5374e-06 6.34866e-06 0
-7.67975e-07 1.55476e-05 0 3.42517e-07 1.84594e-06 0 -3.23587e-07 5.0874e-06 0 -1.03014e-07 1.47467e-05 0
2.1852e-06 2.70154e-06 0 2.1575e-06 1.18539e-05 0 -3.72646e-07 5.55414e-06 0 -4.16966e-07 1.1137e-06 0
-9.794098332e-08 2.545940845e-07 0 1.278802984e-07 1.497640369e-05 0 8.586937383e-08 2.331672704e-06 0 2.644526909e-08 5.194820574e-07 0
4.977602543e-07 2.385048219e-06 0 -3.412924116e-06 3.77589422e-05 0 -6.650605428e-06 1.241584687e-05 0 -3.491896223e-06 2.378009185e-06 0
7.142941172e-08 1.78535322e-06 0 9.691161722e-07 9.167491157e-07 0 -1.819464501e-06 3.642301863e-06 0 -5.231321666e-06 2.576653606e-05 0
......@@ -1249,20 +1249,20 @@
2.172896529e-08 3.099795265e-07 0 4.532921594e-09 2.586529568e-08 0 -2.13380531e-08 1.738166247e-05 0 -6.42608029e-08 2.74189432e-07 0
6.42349562e-07 4.586747764e-06 0 6.236807621e-07 3.6518461e-06 0 2.401273377e-06 1.525748667e-05 0 -4.558905403e-06 2.715373193e-05 0
-1.264782622e-05 9.306331776e-06 0 -5.868747394e-06 7.013909453e-06 0 -1.059488341e-06 2.054737024e-06 0 2.183707579e-07 4.280394478e-06 0
2.509133083e-06 6.195174365e-07 0 1.470599273e-06 6.810921604e-06 0 -2.127700469e-07 3.978108907e-06 0 7.233115866e-07 2.10731082e-06 0
-5.37226974e-09 8.649515621e-06 0 -1.036154003e-06 1.910766514e-06 0 -1.02696356e-06 6.17689409e-07 0 -3.375485448e-06 1.059637998e-05 0
-2.640285629e-06 8.958862963e-06 0 -2.203608176e-07 2.118862085e-06 0 1.135733214e-07 8.319661902e-07 0 4.274000105e-07 4.348061793e-06 0
1.349285526e-06 2.035795887e-06 0 -6.815530469e-07 1.057773957e-08 0 -2.550322506e-06 9.428581507e-06 0 -9.7740417e-07 1.858776614e-06 0
5.944171377e-08 1.410272603e-06 0 2.759835752e-06 4.99375119e-06 0 2.45686897e-06 2.34261006e-06 0 -7.284698e-08 1.584969311e-07 0
1.041478797e-07 5.023975973e-07 0 2.511210084e-07 3.517138794e-07 0 7.142982668e-07 1.043898123e-06 0 5.434893637e-07 2.297776064e-06 0
-2.160872192e-08 4.20829025e-08 0 1.517537029e-08 5.063321851e-08 0 -6.553140253e-08 5.963049716e-07 0 2.163441621e-07 8.636174016e-07 0
2.751651778e-07 3.418507504e-06 0 -1.090806109e-06 2.190308805e-07 0 -7.491701126e-07 2.730883352e-06 0 -1.694137381e-07 2.784076742e-06 0
-5.116960438e-07 4.857181921e-06 0 4.126272657e-08 1.639657086e-07 0 -8.242377589e-07 4.244644174e-07 0 -9.83650807e-07 2.0824084e-06 0
-1.201609336e-07 5.968308869e-07 0 4.871202819e-08 1.744449776e-07 0 -1.405182047e-06 1.646878474e-07 0 -4.7090316e-06 1.183258469e-08 0
-1.855365213e-06 8.704866559e-06 0 4.971587941e-07 1.445205413e-07 0 -9.198865314e-07 1.932642817e-06 0 9.967801873e-08 2.1959886e-06 0
3.178618044e-06 6.369881703e-06 0 1.097572067e-06 1.831581721e-05 0 -2.04117373e-06 3.958675279e-06 0 2.496433638e-07 8.119062613e-07 0
-2.406181636e-07 1.317787451e-06 0 -5.717812996e-07 1.210848836e-06 0 7.856112916e-08 9.372601539e-08 0 5.739914144e-08 4.140760211e-06 0
</DataArray>
2.50913e-06 6.19517e-07 0 1.4706e-06 6.81092e-06 0 -2.1277e-07 8.16899e-06 0 7.23312e-07 3.06577e-06 0
-5.37227e-09 1.91766e-05 0 -1.03615e-06 2.97478e-06 0 -1.02696e-06 2.07294e-06 0 -3.37549e-06 2.27037e-05 0
-2.64029e-06 1.56716e-05 0 -2.20361e-07 4.06396e-06 0 1.13573e-07 1.50376e-06 0 4.274e-07 8.54247e-06 0
1.34929e-06 3.30337e-06 0 -6.81553e-07 2.82022e-06 0 -2.55032e-06 1.79269e-05 0 -9.77404e-07 3.07493e-06 0
5.94417e-08 2.42632e-06 0 2.75984e-06 7.68133e-06 0 2.45687e-06 7.29436e-06 0 -7.2847e-08 2.37572e-07 0
1.04148e-07 9.07222e-07 0 2.51121e-07 6.54028e-07 0 7.14298e-07 1.67402e-06 0 5.43489e-07 5.18014e-06 0
-2.16087e-08 6.46774e-08 0 1.51754e-08 8.39707e-08 0 -6.55314e-08 1.29061e-06 0 2.16344e-07 1.34736e-06 0
2.75165e-07 7.15807e-06 0 -1.09081e-06 1.48298e-06 0 -7.4917e-07 4.07522e-06 0 -1.69414e-07 6.37495e-06 0
-5.11696e-07 9.24985e-06 0 4.12627e-08 2.39484e-07 0 -8.24238e-07 1.80288e-06 0 -9.83651e-07 3.37028e-06 0
-1.20161e-07 1.12471e-06 0 4.8712e-08 2.48972e-07 0 -1.40518e-06 1.88319e-06 0 -4.70903e-06 1.7737e-06 0
-1.85537e-06 1.2806e-05 0 4.97159e-07 2.54022e-06 0 -9.19887e-07 3.03115e-06 0 9.9678e-08 4.20655e-06 0
3.17862e-06 9.84625e-06 0 1.09757e-06 4.16062e-05 0 -2.04117e-06 6.08154e-06 0 2.49643e-07 1.16881e-06 0
-2.40618e-07 3.58084e-06 0 -5.71781e-07 1.80759e-06 0 7.85611e-08 1.51215e-07 0 5.73991e-08 8.33892e-06 0
</DataArray>
<DataArray type="Float32" Name="process rank" NumberOfComponents="1" format="ascii">
0 0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0 0 0
......
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