[tracer][reference] implements an outflow boundary for tracer problem

The tracer problem had neumann-noflow boundaries so far. We now
implement an outflow BC on the top so that the tracer is transported out
of the domain and  does not accumulate anymore.
The reference solution is adapted.

examples/1ptracer/README.md

@@ -665,10 +665,21 @@ We use convenient declarations that we derive from the property system.
     using SpatialParams = GetPropType<TypeTag, Properties::SpatialParams>;
     using Element = typename GridGeometry::GridView::template Codim<0>::Entity;
     using GlobalPosition = typename Element::Geometry::GlobalCoordinate;
+    using NumEqVector = GetPropType<TypeTag, Properties::NumEqVector>;
+    using GridVariables = GetPropType<TypeTag, Properties::GridVariables>;
+    using ElementVolumeVariables = typename GridVariables::GridVolumeVariables::LocalView;
+    using ElementFluxVariablesCache = typename GridVariables::GridFluxVariablesCache::LocalView;
+    using FVElementGeometry = typename GetPropType<TypeTag, Properties::GridGeometry>::LocalView;
+    using SubControlVolumeFace = typename FVElementGeometry::SubControlVolumeFace;
 ```
 We create a bool saying whether mole or mass fractions are used
 ```cpp
     static constexpr bool useMoles = getPropValue<TypeTag, Properties::UseMoles>();
+```
+We create additional int to make dimWorld and numComponents available in the problem
+```cpp
+    static constexpr int dimWorld = GridView::dimensionworld;
+    static const int numComponents = FluidSystem::numComponents;
 
 public:
 ```
@@ -716,8 +727,37 @@ We assign different values, depending on wether mole concentrations or mass conc
         }
         return initialValues;
     }
+```
+We implement an outflow boundary on the top of the domain and prescribe zero-flux Neumann boundary conditions on all other boundaries.
+```cpp
+        NumEqVector neumann(const Element& element,
+                            const FVElementGeometry& fvGeometry,
+                            const ElementVolumeVariables& elemVolVars,
+                            const ElementFluxVariablesCache& elemFluxVarsCache,
+                            const SubControlVolumeFace& scvf) const
+        {
+            NumEqVector values(0.0);
+            const auto& volVars = elemVolVars[scvf.insideScvIdx()];
+            const auto& globalPos = scvf.center();
+```
+This is the outflow boundary, where tracer is transported by advection
+with the given flux field and diffusive flux is enforced to be zero
+```cpp
+            if (globalPos[dimWorld-1] > this->gridGeometry().bBoxMax()[dimWorld-1] - eps_)
+            {
+                values = this->spatialParams().volumeFlux(element, fvGeometry, elemVolVars, scvf)
+                         * volVars.massFraction(0, 0) * volVars.density(0)
+                         / scvf.area();
+                assert(values>=0.0 && "Volume flux at outflow boundary is expected to have a positive sign");
+            }
+```
+prescribe zero-flux Neumann boundary conditions elsewhere
+```cpp
+            else
+                values = 0.0;
 
-
+            return values;
+        }
 
 private:
 ```

examples/1ptracer/problem_tracer.hh

@@ -158,8 +158,17 @@ class TracerTestProblem : public PorousMediumFlowProblem<TypeTag>
     using SpatialParams = GetPropType<TypeTag, Properties::SpatialParams>;
     using Element = typename GridGeometry::GridView::template Codim<0>::Entity;
     using GlobalPosition = typename Element::Geometry::GlobalCoordinate;
+    using NumEqVector = GetPropType<TypeTag, Properties::NumEqVector>;
+    using GridVariables = GetPropType<TypeTag, Properties::GridVariables>;
+    using ElementVolumeVariables = typename GridVariables::GridVolumeVariables::LocalView;
+    using ElementFluxVariablesCache = typename GridVariables::GridFluxVariablesCache::LocalView;
+    using FVElementGeometry = typename GetPropType<TypeTag, Properties::GridGeometry>::LocalView;
+    using SubControlVolumeFace = typename FVElementGeometry::SubControlVolumeFace;
     //We create a bool saying whether mole or mass fractions are used
     static constexpr bool useMoles = getPropValue<TypeTag, Properties::UseMoles>();
+    //We create additional int to make dimWorld and numComponents available in the problem
+    static constexpr int dimWorld = GridView::dimensionworld;
+    static const int numComponents = FluidSystem::numComponents;
 
 public:
     // This is the constructor of our problem class:
@@ -198,7 +207,32 @@ public:
         return initialValues;
     }
 
+        //We implement an outflow boundary on the top of the domain and prescribe zero-flux Neumann boundary conditions on all other boundaries.
+        NumEqVector neumann(const Element& element,
+                            const FVElementGeometry& fvGeometry,
+                            const ElementVolumeVariables& elemVolVars,
+                            const ElementFluxVariablesCache& elemFluxVarsCache,
+                            const SubControlVolumeFace& scvf) const
+        {
+            NumEqVector values(0.0);
+            const auto& volVars = elemVolVars[scvf.insideScvIdx()];
+            const auto& globalPos = scvf.center();
+
+            // This is the outflow boundary, where tracer is transported by advection
+            // with the given flux field and diffusive flux is enforced to be zero
+            if (globalPos[dimWorld-1] > this->gridGeometry().bBoxMax()[dimWorld-1] - eps_)
+            {
+                values = this->spatialParams().volumeFlux(element, fvGeometry, elemVolVars, scvf)
+                         * volVars.massFraction(0, 0) * volVars.density(0)
+                         / scvf.area();
+                assert(values>=0.0 && "Volume flux at outflow boundary is expected to have a positive sign");
+            }
+            //prescribe zero-flux Neumann boundary conditions elsewhere
+            else
+                values = 0.0;
 
+            return values;
+        }
 
 private:
 // We assign a private global variable for the epsilon:

test/porousmediumflow/tracer/1ptracer/problem_tracer.hh

@@ -158,9 +158,17 @@ class TracerTestProblem : public PorousMediumFlowProblem<TypeTag>
     using PrimaryVariables = GetPropType<TypeTag, Properties::PrimaryVariables>;
     using FluidSystem = GetPropType<TypeTag, Properties::FluidSystem>;
     using SpatialParams = GetPropType<TypeTag, Properties::SpatialParams>;
+    using NumEqVector = GetPropType<TypeTag, Properties::NumEqVector>;
+    using GridVariables = GetPropType<TypeTag, Properties::GridVariables>;
+    using ElementVolumeVariables = typename GridVariables::GridVolumeVariables::LocalView;
+    using ElementFluxVariablesCache = typename GridVariables::GridFluxVariablesCache::LocalView;
+    using FVElementGeometry = typename GetPropType<TypeTag, Properties::GridGeometry>::LocalView;
+    using SubControlVolumeFace = typename FVElementGeometry::SubControlVolumeFace;
 
     //! property that defines whether mole or mass fractions are used
     static constexpr bool useMoles = getPropValue<TypeTag, Properties::UseMoles>();
+    static constexpr int dimWorld = GridView::dimensionworld;
+    static const int numComponents = FluidSystem::numComponents;
 
     using Element = typename GridGeometry::GridView::template Codim<0>::Entity;
     using GlobalPosition = typename Element::Geometry::GlobalCoordinate;
@@ -193,6 +201,43 @@ public:
         values.setAllNeumann();
         return values;
     }
+
+    /*!
+     * \brief Evaluates the boundary conditions for a Neumann boundary segment.
+     *
+     * Implements an outflow boundary on the top, where tracer is transported by
+     * advection with the given flux field and diffusive flux is enforced to be zero.
+     *
+     * \param element The element
+     * \param fvGeometry The finite volume geometry
+     * \param elemVolVars The element volume variables
+     * \param elemFluxVarsCache Flux variables caches for all faces in stencil
+     * \param scvf The subcontrolvolume face
+     */
+    NumEqVector neumann(const Element& element,
+                        const FVElementGeometry& fvGeometry,
+                        const ElementVolumeVariables& elemVolVars,
+                        const ElementFluxVariablesCache& elemFluxVarsCache,
+                        const SubControlVolumeFace& scvf) const
+    {
+        NumEqVector values(0.0);
+        const auto& volVars = elemVolVars[scvf.insideScvIdx()];
+        const auto& globalPos = scvf.center();
+
+        // This is the outflow boundary, where tracer is transported by advection
+        // with the given flux field and diffusive flux is enforced to be zero
+        if (globalPos[dimWorld-1] > this->gridGeometry().bBoxMax()[dimWorld-1] - eps_)
+        {
+            values = this->spatialParams().volumeFlux(element, fvGeometry, elemVolVars, scvf)
+                     * volVars.massFraction(0, 0) * volVars.density(0)
+                     / scvf.area();
+            assert(values>=0.0 && "Volume flux at outflow boundary is expected to have a positive sign");
+        }
+        else
+            values = 0.0;
+
+        return values;
+    }
     // \}
 
     /*!

test/references/test_1ptracer_transport-reference.vtu

@@ -208,11 +208,11 @@
           2.90532e-17 5.45141e-17 8.81528e-18 1.38004e-18 5.16647e-20 3.37864e-23 3.37094e-24 2.74128e-25 2.9408e-26 1.47463e-27 1.34525e-43 0
           2.47473e-32 7.6169e-29 3.34277e-27 4.23311e-27 7.01946e-25 2.16138e-23 5.24325e-24 5.09109e-24 1.4758e-24 1.94278e-23 1.01674e-22 1.69944e-21
           2.67936e-22 1.55879e-20 5.47453e-18 1.48655e-19 2.15992e-18 1.60201e-18 1.04019e-16 2.11704e-15 3.03978e-15 5.03726e-16 2.53974e-16 1.47194e-15
-          1.54598e-16 3.90078e-17 6.88191e-18 1.82039e-17 2.33099e-15 1.74998e-15 1.26179e-15 1.39141e-14 2.89783e-15 1.14903e-15 7.69323e-18 3.1605e-16
-          1.44819e-17 3.67297e-18 1.92558e-17 7.58995e-18 9.11791e-19 4.36327e-19 8.43538e-20 3.46616e-25 1.04941e-25 6.3449e-27 1.76634e-27 2.43065e-28
-          0 0 1.25917e-32 3.1183e-30 5.76142e-28 5.26306e-27 9.52037e-26 3.35952e-24 1.34232e-24 4.75994e-26 1.57447e-25 2.34604e-24
-          3.8149e-24 1.93092e-23 5.38689e-20 5.442e-19 2.63823e-18 2.74714e-19 2.96745e-19 2.61456e-19 4.48102e-17 3.68616e-15 9.62048e-16 3.43464e-17
-          3.84133e-17 1.55783e-16 1.52386e-18 1.3212e-17
+          1.54598e-16 3.90078e-17 5.22201e-18 1.70196e-17 1.19091e-15 1.60222e-15 6.19484e-16 6.75743e-15 1.94209e-15 9.75636e-16 7.30102e-18 2.34301e-16
+          1.16363e-17 2.14341e-18 1.22764e-17 6.98632e-18 8.5158e-19 3.75666e-19 5.58315e-20 3.44932e-25 1.02482e-25 6.22383e-27 1.70753e-27 2.09015e-28
+          0 0 1.19647e-32 3.05053e-30 4.87241e-28 4.6631e-27 8.94165e-26 2.81382e-24 1.09607e-24 4.74032e-26 1.3892e-25 2.18913e-24
+          3.69809e-24 1.92201e-23 3.64594e-20 4.04697e-19 2.16077e-18 1.97036e-19 2.76122e-19 2.28842e-19 3.65312e-17 1.19957e-15 8.28343e-16 3.3405e-17
+          3.20417e-17 1.48511e-16 1.50453e-18 9.95476e-18
         </DataArray>
         <DataArray type="Float32" Name="X^tracer_0" NumberOfComponents="1" format="ascii">
           4.17447e-13 4.2631e-13 4.0993e-13 6.60598e-14 2.69839e-16 6.9659e-12 3.14842e-17 1.03942e-17 1.23503e-19 5.05161e-20 2.66604e-13 8.47736e-12
@@ -419,11 +419,11 @@
           4.8422e-19 9.08568e-19 1.46921e-19 2.30007e-20 8.61078e-22 5.63107e-25 5.61824e-26 4.5688e-27 4.90133e-28 2.45772e-29 2.8026e-45 0
           4.12455e-34 1.26948e-30 5.57129e-29 7.05518e-29 1.16991e-26 3.6023e-25 8.73875e-26 8.48516e-26 2.45966e-26 3.23797e-25 1.69456e-24 2.8324e-23
           4.46561e-24 2.59799e-22 9.12422e-20 2.47759e-21 3.59987e-20 2.67001e-20 1.73365e-18 3.5284e-17 5.0663e-17 8.39543e-18 4.2329e-18 2.45323e-17
-          2.57663e-18 6.50129e-19 1.14699e-19 3.03399e-19 3.88498e-17 2.91663e-17 2.10298e-17 2.31902e-16 4.82972e-17 1.91505e-17 1.2822e-19 5.26751e-18
-          2.41365e-19 6.12161e-20 3.2093e-19 1.26499e-19 1.51965e-20 7.27211e-21 1.4059e-21 5.77693e-27 1.74902e-27 1.05748e-28 2.94389e-29 4.05108e-30
-          0 0 2.09862e-34 5.19717e-32 9.60237e-30 8.77177e-29 1.58673e-27 5.59921e-26 2.23721e-26 7.93323e-28 2.62412e-27 3.91006e-26
-          6.35817e-26 3.2182e-25 8.97814e-22 9.07001e-21 4.39705e-20 4.57857e-21 4.94576e-21 4.35761e-21 7.46837e-19 6.14361e-17 1.60341e-17 5.7244e-19
-          6.40222e-19 2.59638e-18 2.53977e-20 2.202e-19
+          2.57663e-18 6.50129e-19 8.70335e-20 2.8366e-19 1.98485e-17 2.67037e-17 1.03247e-17 1.12624e-16 3.23682e-17 1.62606e-17 1.21684e-19 3.90501e-18
+          1.93938e-19 3.57236e-20 2.04606e-19 1.16439e-19 1.4193e-20 6.26109e-21 9.30525e-22 5.74887e-27 1.70804e-27 1.03731e-28 2.84589e-29 3.48358e-30
+          0 0 1.99412e-34 5.08422e-32 8.12068e-30 7.77184e-29 1.49028e-27 4.6897e-26 1.82678e-26 7.90053e-28 2.31534e-27 3.64855e-26
+          6.16349e-26 3.20335e-25 6.07657e-22 6.74495e-21 3.60129e-20 3.28394e-21 4.60203e-21 3.81403e-21 6.08853e-19 1.99929e-17 1.38057e-17 5.5675e-19
+          5.34029e-19 2.47519e-18 2.50755e-20 1.65913e-19
         </DataArray>
         <DataArray type="Float32" Name="rho" NumberOfComponents="1" format="ascii">
           1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000
@@ -1262,7 +1262,7 @@
           -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>
         <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