[singlefracture] use unit fluid and updated parameters

test/porousmediumflow/tracer/1ptracer/singlefracture/1ptestproblem.hh

@@ -29,7 +29,7 @@
 
 #include <dumux/porousmediumflow/1p/model.hh>
 #include <dumux/porousmediumflow/problem.hh>
-#include <dumux/material/components/simpleh2o.hh>
+#include <dumux/material/components/constant.hh>
 #include <dumux/material/fluidsystems/liquidphase.hh>
 #include <dumux/porousmediumflow/1p/incompressiblelocalresidual.hh>
 
@@ -61,7 +61,7 @@ SET_TYPE_PROP(IncompressibleTestProblem, LocalResidual, OnePIncompressibleLocalR
 SET_PROP(IncompressibleTestProblem, FluidSystem)
 {
     using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar);
-    using type = FluidSystems::LiquidPhase<Scalar, SimpleH2O<Scalar> >;
+    using type = FluidSystems::LiquidPhase<Scalar, Dumux::Components::Constant<0, Scalar> >;
 };
 
 // Enable caching
@@ -128,7 +128,7 @@ class OnePTestProblem : public PorousMediumFlowProblem<TypeTag>
 public:
     OnePTestProblem(std::shared_ptr<const FVGridGeometry> fvGridGeometry) : ParentType(fvGridGeometry)
     {
-        overPressure_ = getParam<Scalar>("Problem.OverPressure");
+        headAtInlet_ = getParam<Scalar>("Problem.HeadAtInlet");
     }
 
     //! Specifies which kind of boundary condition should be used
@@ -145,8 +145,9 @@ public:
     PrimaryVariables dirichletAtPos(const GlobalPosition& globalPos) const
     {
         if (isOnInflowBoundary(globalPos))
-            return PrimaryVariables(1.0e+5+overPressure_);
-        return PrimaryVariables(1.0e+5);
+            return PrimaryVariables(headAtInlet_);
+        else
+            return PrimaryVariables(1.0);
     }
 
     //! Returns the temperature in the domain
@@ -180,7 +181,7 @@ public:
     }
 
 private:
-    Scalar overPressure_;
+    Scalar headAtInlet_;
     std::vector<Scalar> kField_;
     std::vector<Scalar> poro_;
 };

test/porousmediumflow/tracer/1ptracer/singlefracture/singlefracture.input

@@ -3,9 +3,9 @@ File = ./single_equi.msh
 DomainMarkers = true
 
 [SpatialParams]
-BottomLayerPermeability = 1e-12
-UpperLayerPermeability = 1e-13
-FracturePermeability = 1e-10
+BottomLayerPermeability = 1e-5
+UpperLayerPermeability = 1e-6
+FracturePermeability = 1e-1
 BottomLayerPorosity = 0.25
 UpperLayerPorosity = 0.2
 FracturePorosity = 0.4
@@ -16,12 +16,12 @@ Name = single_reference.txt
 [TimeLoop]
 TEnd = 1e8 # [s]
 NumOutputFiles = 1
-MaxTimeStepSize = 1e8
+MaxTimeStepSize = 1e6
 
 [Problem]
 Name = singlefracture_reference # name passed to the output routines
-OverPressure = 3e5
-BoundaryMoleFrac = 1e-2
+HeadAtInlet = 4
+ConcentrationAtInlet = 1e-2
 EnableGravity = false # disable gravity
 
 [LinearSolver]

test/porousmediumflow/tracer/1ptracer/singlefracture/test_singlefracture.cc

@@ -171,8 +171,7 @@ int main(int argc, char** argv) try
         volumeFlux.resize(fvGridGeometry->numScvf(), 0.0);
 
         using FluxVariables =  typename GET_PROP_TYPE(OnePTypeTag, FluxVariables);
-        auto upwindTerm = [](const auto& volVars) { return volVars.mobility(0); };
-        auto massUpwindTerm = [](const auto& volVars) { return volVars.density(0)*volVars.mobility(0); };
+        auto upwindTerm = [](const auto& volVars) { return 1.0; };
         for (const auto& element : elements(leafGridView))
         {
             auto fvGeometry = localView(*fvGridGeometry);
@@ -225,9 +224,9 @@ int main(int argc, char** argv) try
                             sumInflux += volumeFlux[idx];
 
                         if (problemOneP->isOnInflowBoundary(scvf.ipGlobal()))
-                            massInflux += fluxVars.advectiveFlux(0, massUpwindTerm);
+                            massInflux += volumeFlux[idx];
                         else if (problemOneP->isOnOutflowBoundary(scvf.ipGlobal()))
-                            massOutflux += fluxVars.advectiveFlux(0, massUpwindTerm);
+                            massOutflux += volumeFlux[idx];
                         else
                             DUNE_THROW(Dune::InvalidStateException, "Wrong Dirichlet BCs set");
                     }
@@ -284,10 +283,10 @@ int main(int argc, char** argv) try
 
         //! add caption for tracer data
         file << std::endl << "time [s] \t | \t "
-                             "tracer influx [kg/s] \t | \t "
-                             "tracer Outflux [kg/s] \t | \t "
-                             "tracer mass fracture [kg] \t | \t "
-                             "tracer mass lower layer [kg]" << std::endl;
+                             "tracer influx [1/s] \t | \t "
+                             "tracer Outflux [1/s] \t | \t "
+                             "poro*tracer integral fracture [-] \t | \t "
+                             "poro*tracer integral lower matrix layer [-]" << std::endl;
     }
 
     ////////////////////////////////////////////////////////////

test/porousmediumflow/tracer/1ptracer/singlefracture/tracertestproblem.hh

@@ -64,10 +64,49 @@ SET_TYPE_PROP(TracerTestProblem, Problem, TracerTestProblem<TypeTag>);
 SET_TYPE_PROP(TracerTestProblem, SpatialParams, TracerTestSpatialParams<TypeTag>);
 
 // Define whether mole(true) or mass (false) fractions are used
-SET_BOOL_PROP(TracerTestProblem, UseMoles, false);
 SET_BOOL_PROP(TracerTestCCProblem, EnableMolecularDiffusion, false);
 SET_BOOL_PROP(TracerTestCCProblem, SolutionDependentAdvection, false);
 
+// use the static interaction volume type with sizes known at compile time
+SET_PROP(TracerTestProblem, PrimaryInteractionVolume)
+{
+private:
+    using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar);
+    using NodalIndexSet = typename GET_PROP_TYPE(TypeTag, DualGridNodalIndexSet);
+
+    // use the default traits
+    using Traits = CCMpfaODefaultStaticInteractionVolumeTraits< NodalIndexSet, Scalar, 8, 12 >;
+public:
+    using type = CCMpfaOStaticInteractionVolume< Traits >;
+};
+
+// sizes of the ivs are known, use Dune::ReservedVector in index sets
+SET_PROP(TracerTestProblem, DualGridNodalIndexSet)
+{
+    using GV = typename GET_PROP_TYPE(TypeTag, GridView);
+    static constexpr int dim = GV::dimension;
+    static constexpr int dimWorld = GV::dimensionworld;
+private:
+    struct Traits
+    {
+        using GridView = GV;
+        using GridIndexType = typename GV::IndexSet::IndexType;
+        using LocalIndexType = std::uint8_t;
+
+        //! per default, we use dynamic data containers (iv size unknown)
+        template< class T > using NodalScvDataStorage = Dune::ReservedVector< T, 8 >;
+        template< class T > using NodalScvfDataStorage = Dune::ReservedVector< T, 24 >;
+
+        //! store data on neighbors of scvfs in static containers if possible
+        template< class T >
+        using ScvfNeighborDataStorage = typename std::conditional_t< (dim<dimWorld),
+                                                                     std::vector< T >,
+                                                                     Dune::ReservedVector< T, 2 > >;
+    };
+public:
+    using type = CCMpfaDualGridNodalIndexSet< Traits >;
+};
+
 //! A simple fluid system with one tracer component
 template<class TypeTag>
 class TracerFluidSystem : public FluidSystems::BaseFluidSystem<typename GET_PROP_TYPE(TypeTag, Scalar),
@@ -88,9 +127,9 @@ public:
     //! The number of components
     static constexpr int numComponents = 1;
     //! Human readable component name (index compIdx) (for vtk output)
-    static std::string componentName(int compIdx) { return "tracer_" + std::to_string(compIdx); }
+    static std::string componentName(int compIdx) { return "c"; }
     //! Human readable phase name (index phaseIdx) (for velocity vtk output)
-    static std::string phaseName(int phaseIdx = 0) { return "Groundwater"; }
+    static std::string phaseName(int phaseIdx = 0) { return "h2o"; }
     //! Molar mass in kg/mol of the component with index compIdx
     static constexpr Scalar molarMass(unsigned int compIdx) { return 1.; }
     //! binary diffusion coefficient
@@ -148,12 +187,8 @@ class TracerTestProblem : public PorousMediumFlowProblem<TypeTag>
 public:
     TracerTestProblem(std::shared_ptr<const FVGridGeometry> fvGridGeom) : ParentType(fvGridGeom)
     {
-        boundaryMoleFrac_ = getParam<Scalar>("Problem.BoundaryMoleFrac");
+        concAtInlet_ = getParam<Scalar>("Problem.ConcentrationAtInlet");
         outputFileName_ = getParam<std::string>("OutputFile.Name");
-
-        // state in the console whether mole or mass fractions are used
-        if(useMoles) std::cout<<"problem uses mole fractions" << '\n';
-        else std::cout<<"problem uses mass fractions" << '\n';
     }
 
     //! writes the mass distribution in the layers to output file
@@ -188,12 +223,7 @@ public:
             const auto& vv = elemVolVars[eIdx];
             for (const auto& scv : scvs(fvGeometry))
             {
-                const auto mass = this->spatialParams().porosity(element, scv, elemSol)
-                                  *FluidSystem::molarMass(0)
-                                  *vv.moleFraction(0, 0)
-                                  /this->spatialParams().fluidMolarMass()
-                                  *this->spatialParams().fluidDensity()
-                                  *scv.volume();
+                const auto mass = vv.porosity()*vv.moleFraction(0, 0)*scv.volume();
 
                 if (marker == SpatialParams::LayerMarkers::lowerLayer)
                     massLowerLayer += mass;
@@ -207,7 +237,7 @@ public:
 
             // evaluate influx/outflux of tracer
             const auto c = element.geometry().center();
-            if (c[2] > 90.0 || c[2] < 10.0)
+            if (c[2] > 80.0 || c[2] < 10.0) // discard those elements that are clearly not connected to inflow/outflow boundaries
             {
                 for (const auto& scvf : scvfs(fvGeometry))
                 {
@@ -245,23 +275,15 @@ public:
         //! influx
         if (isOnInflowBoundary(scvf.ipGlobal()))
         {
-            const auto influx = this->spatialParams().volumeFlux(scvf)
-                                *this->spatialParams().fluidDensity()
-                                /this->spatialParams().fluidMolarMass()
-                                *boundaryMoleFrac_
-                                *FluidSystem::molarMass(0)
-                                /scvf.area();
+            const auto influx = this->spatialParams().volumeFlux(scvf)*concAtInlet_/scvf.area();
             return NumEqVector(influx);
         }
         //! outflux
         else if (isOnOutflowBoundary(scvf.ipGlobal()))
         {
             const auto outFlux = this->spatialParams().volumeFlux(scvf)
-                                 *this->spatialParams().fluidDensity()
-                                 /this->spatialParams().fluidMolarMass()
                                  *elemVolVars[scvf.insideScvIdx()].moleFraction(0, 0)
-                                 *FluidSystem::molarMass(0)/
-                                 scvf.area();
+                                 /scvf.area();
             return NumEqVector(outFlux);
         }
         //! remaining boundaries
@@ -269,7 +291,7 @@ public:
     }
 
     //! Evaluate the initial value for a control volume.
-    PrimaryVariables initialAtPos(const GlobalPosition &globalPos) const { return PrimaryVariables(0.0); }
+    PrimaryVariables initialAtPos(const GlobalPosition& globalPos) const { return PrimaryVariables(0.0); }
 
     //! Returns if a given position is on inflow boundary
     bool isOnInflowBoundary(const GlobalPosition& globalPos) const { return globalPos[0] < 1.0e-6 && globalPos[2] > 90.0; }
@@ -278,7 +300,7 @@ public:
     bool isOnOutflowBoundary(const GlobalPosition& globalPos) const { return globalPos[1] < 1.0e-6 && globalPos[2] < 10.0; }
 
 private:
-    Scalar boundaryMoleFrac_;
+    Scalar concAtInlet_;
     std::string outputFileName_;
 };
 

test/porousmediumflow/tracer/1ptracer/singlefracture/tracertestspatialparams.hh

@@ -93,11 +93,11 @@ public:
     //! They can possible vary with space but are usually constants
 
     //! fluid density
-    static constexpr Scalar fluidDensity() { return 1000.0; }
+    static constexpr Scalar fluidDensity() { return 1.0; }
     Scalar fluidDensity(const Element &element, const SubControlVolume& scv) const { return fluidDensity(); }
 
     //! fluid molar mass
-    static constexpr Scalar fluidMolarMass() { return 18.0e-3; }
+    static constexpr Scalar fluidMolarMass() { return 1.; }
     Scalar fluidMolarMass(const Element &element, const SubControlVolume& scv) const { return fluidMolarMass(); }
     Scalar fluidMolarMass(const GlobalPosition &globalPos) const { return fluidMolarMass(); }