From 9a6338a31bcddc00fb07399d6a64cfbffcc0b721 Mon Sep 17 00:00:00 2001 From: Simon Scholz Date: Thu, 31 Oct 2019 10:25:28 +0100 Subject: [PATCH 1/3] [flux] Introduce flux traits --- dumux/flux/CMakeLists.txt | 1 + dumux/flux/stationaryvelocityfield.hh | 8 ++++ dumux/flux/traits.hh | 53 +++++++++++++++++++++++++++ 3 files changed, 62 insertions(+) create mode 100644 dumux/flux/traits.hh diff --git a/dumux/flux/CMakeLists.txt b/dumux/flux/CMakeLists.txt index 30051be02f..e32e550170 100644 --- a/dumux/flux/CMakeLists.txt +++ b/dumux/flux/CMakeLists.txt @@ -18,5 +18,6 @@ maxwellstefanslaw.hh referencesystemformulation.hh shallowwaterflux.hh stationaryvelocityfield.hh +traits.hh upwindscheme.hh DESTINATION ${CMAKE_INSTALL_INCLUDEDIR}/dumux/flux) diff --git a/dumux/flux/stationaryvelocityfield.hh b/dumux/flux/stationaryvelocityfield.hh index 5f766d0fbf..4f04257adb 100644 --- a/dumux/flux/stationaryvelocityfield.hh +++ b/dumux/flux/stationaryvelocityfield.hh @@ -27,6 +27,9 @@ #ifndef DUMUX_DISCRETIZATION_STATIONARY_VELOCITY_FIELD_HH #define DUMUX_DISCRETIZATION_STATIONARY_VELOCITY_FIELD_HH +#include +#include + #include #include @@ -64,6 +67,11 @@ public: } }; +//! Set stationary velocity field to true in the FluxTraits +template +struct HasStationaryVelocityField> +: public std::true_type {}; + } // end namespace Dumux #endif diff --git a/dumux/flux/traits.hh b/dumux/flux/traits.hh new file mode 100644 index 0000000000..831a0a2527 --- /dev/null +++ b/dumux/flux/traits.hh @@ -0,0 +1,53 @@ +// -*- 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 . * + *****************************************************************************/ +/*! + * \file + * \ingroup Flux + * \brief Defines the flux traits. + */ +#ifndef DUMUX_FLUX_TRAITS_HH +#define DUMUX_FLUX_TRAITS_HH + +#include + +namespace Dumux { + +/*! + * \ingroup Flux + * \brief Trait of an advection type stating whether it implements a stationary velocity field + */ +template +struct HasStationaryVelocityField : public std::false_type {}; + +/*! + * \ingroup Flux + * \brief Traits of a flux variables type + */ +template +struct FluxTraits +{ + static constexpr bool hasStationaryVelocityField() + { + return HasStationaryVelocityField::value; + } +}; + +} // namespace Dumux + +#endif // DUMUX_FLUX_TRAITS_HH -- GitLab From a45789d56ed134aca819b2facb241972eec20322 Mon Sep 17 00:00:00 2001 From: Simon Scholz Date: Thu, 31 Oct 2019 10:32:57 +0100 Subject: [PATCH 2/3] [1ptracer][reference] fix velocity calculation on boundary for stationary fields This uses the newly introduced fluxtraits to check if the velocity field is stationary. With these changes the velocity output is corrected in the reference solution. --- dumux/porousmediumflow/velocity.hh | 46 +++++++++++----- .../test_1ptracer_transport-reference.vtu | 54 +++++++++---------- 2 files changed, 60 insertions(+), 40 deletions(-) diff --git a/dumux/porousmediumflow/velocity.hh b/dumux/porousmediumflow/velocity.hh index 5531251ccb..8ef2c6dadf 100644 --- a/dumux/porousmediumflow/velocity.hh +++ b/dumux/porousmediumflow/velocity.hh @@ -26,6 +26,7 @@ #define DUMUX_POROUSMEDIUMFLOW_VELOCITY_HH #include +#include #include #include @@ -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; + 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; @@ -282,21 +286,36 @@ 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; - if (Dune::FloatCmp::eq(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; + if (Dune::FloatCmp::eq(neumannFlux, NumEqVector(0.0), 1e-30)) + scvfFluxes[scvfIndexInInside[localScvfIdx]] = 0; + // 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]; + } + // simplices + else if (geomType.isSimplex()) + scvfFluxes[scvfIndexInInside[localScvfIdx]] = 0; } - // simplices - else if (geomType.isSimplex()) - scvfFluxes[scvfIndexInInside[localScvfIdx]] = 0; } } @@ -304,6 +323,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 diff --git a/test/references/test_1ptracer_transport-reference.vtu b/test/references/test_1ptracer_transport-reference.vtu index 219faa70d7..043d35e1bc 100644 --- a/test/references/test_1ptracer_transport-reference.vtu +++ b/test/references/test_1ptracer_transport-reference.vtu @@ -637,19 +637,19 @@ 1000 1000 1000 1000 - -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 - + 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 + 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 -- GitLab From 35a86158663efa1e5d524485340520c396e6737c Mon Sep 17 00:00:00 2001 From: Timo Koch Date: Wed, 20 Nov 2019 10:48:59 +0100 Subject: [PATCH 3/3] [pmflow][velocity][doc] Improve comments and add not implemented error message --- dumux/porousmediumflow/velocity.hh | 12 +++++++++--- 1 file changed, 9 insertions(+), 3 deletions(-) diff --git a/dumux/porousmediumflow/velocity.hh b/dumux/porousmediumflow/velocity.hh index 8ef2c6dadf..771a9207a1 100644 --- a/dumux/porousmediumflow/velocity.hh +++ b/dumux/porousmediumflow/velocity.hh @@ -305,16 +305,22 @@ public: using NumEqVector = std::decay_t; if (Dune::FloatCmp::eq(neumannFlux, NumEqVector(0.0), 1e-30)) scvfFluxes[scvfIndexInInside[localScvfIdx]] = 0; - // cubes + + // 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]; } - // simplices + + // 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"); } } } @@ -344,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); -- GitLab