From a376fbcfff69300737bbc7e2eb20a6bc8eeff17d Mon Sep 17 00:00:00 2001
From: Ned Coltman <edward.coltman@iws.uni-stuttgart.de>
Date: Thu, 19 Mar 2020 08:35:42 +0100
Subject: [PATCH] [test][containerstest] Create a unit test for the diffusion
 containers

---
 dumux/flux/fickiandiffusioncoefficients.hh |   7 +-
 test/CMakeLists.txt                        |   1 +
 test/flux/CMakeLists.txt                   |   2 +
 test/flux/fluidsystems.hh                  |  81 ++++++++++++
 test/flux/test_diffusioncontainers.cc      | 145 +++++++++++++++++++++
 5 files changed, 234 insertions(+), 2 deletions(-)
 create mode 100644 test/flux/CMakeLists.txt
 create mode 100644 test/flux/fluidsystems.hh
 create mode 100644 test/flux/test_diffusioncontainers.cc

diff --git a/dumux/flux/fickiandiffusioncoefficients.hh b/dumux/flux/fickiandiffusioncoefficients.hh
index 976f4bb669..11cb6386c2 100644
--- a/dumux/flux/fickiandiffusioncoefficients.hh
+++ b/dumux/flux/fickiandiffusioncoefficients.hh
@@ -69,11 +69,14 @@ public:
 
     const Scalar& operator()(int phaseIdx, int compIIdx, int compJIdx) const
     {
+        std::cout << phaseIdx << " " << compIIdx << " " << compJIdx << " | " << getIndex_(phaseIdx, compIIdx, compJIdx) << " | ";
         sortComponentIndices_(compIIdx, compJIdx);
-        assert(phaseIdx != compJIdx);
+        std::cout << phaseIdx << " " << compIIdx << " " << compJIdx << " | " << getIndex_(phaseIdx, compIIdx, compJIdx) << " | ";
+
         return diffCoeff_[getIndex_(phaseIdx, compIIdx, compJIdx)];
     }
 
+
 private:
     /*
      * \brief For the mpnc case, we have a rectangular diffCoeff matrix (numPhases x numComponents).
@@ -116,13 +119,13 @@ public:
 
     const Scalar& operator()(int phaseIdx, int compIIdx, int compJIdx) const
     {
-        sortComponentIndices_(compIIdx, compJIdx);
         assert(phaseIdx == 0);
         assert(phaseIdx != compJIdx);
         assert(compIIdx != compJIdx);
         assert(compIIdx < numComponents);
         assert(compJIdx < numComponents);
         assert(compIIdx <= compJIdx);
+        sortComponentIndices_(compIIdx, compJIdx);
         return diffCoeff_[compJIdx-1];
     }
 
diff --git a/test/CMakeLists.txt b/test/CMakeLists.txt
index a543ada950..9cec5e3ff8 100644
--- a/test/CMakeLists.txt
+++ b/test/CMakeLists.txt
@@ -1,5 +1,6 @@
 add_subdirectory(common)
 add_subdirectory(geomechanics)
+add_subdirectory(flux)
 add_subdirectory(freeflow)
 add_subdirectory(io)
 add_subdirectory(linear)
diff --git a/test/flux/CMakeLists.txt b/test/flux/CMakeLists.txt
new file mode 100644
index 0000000000..0d26745087
--- /dev/null
+++ b/test/flux/CMakeLists.txt
@@ -0,0 +1,2 @@
+dumux_add_test(SOURCES test_diffusioncontainers.cc
+               LABELS unit flux)
diff --git a/test/flux/fluidsystems.hh b/test/flux/fluidsystems.hh
new file mode 100644
index 0000000000..95be775fc7
--- /dev/null
+++ b/test/flux/fluidsystems.hh
@@ -0,0 +1,81 @@
+// -*- 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
+ * \brief Fluid systems for testing the diffusion containers
+ */
+#ifndef DUMUX_FLUIDSYSTEMS_HH
+#define DUMUX_FLUIDSYSTEMS_HH
+
+#include <string>
+#include <vector>
+
+namespace Dumux {
+
+/*!
+ * \brief A Fluid System for testing the MpNc container
+ */
+class MpNcFluidSystem
+{
+public:
+    //! Constructor
+    MpNcFluidSystem(const int numPhases, const int numComponents)
+    : numPhases_(numPhases)
+    , numComponents_(numComponents)
+    {
+        fillDiffCoeffMatrix_();
+    }
+
+    int binaryDiffusionCoefficient(const int phaseIdx, const int compIIdx, const int compJIdx) const
+    { return binaryDiffusionCoefficients_[phaseIdx][compJIdx]; }
+
+    int numPhases() const
+    { return numPhases_; }
+
+    int numComponents() const
+    { return numComponents_; }
+
+private:
+    void fillDiffCoeffMatrix_()
+    {
+        int count = 0;
+        binaryDiffusionCoefficients_.resize(numPhases_, std::vector<int>(numComponents_, -1));
+        for (unsigned int phaseIdx = 0; phaseIdx < numPhases_; ++ phaseIdx)
+        {
+            for (unsigned int compJIdx = 0; compJIdx < numComponents_; ++compJIdx)
+            {
+                if (compJIdx == phaseIdx)
+                    continue;
+                else
+                {
+                    binaryDiffusionCoefficients_[phaseIdx][compJIdx] = count;
+                    count++;
+                }
+            }
+        }
+    }
+
+    int numPhases_;
+    int numComponents_;
+    std::vector<std::vector<int>> binaryDiffusionCoefficients_;
+};
+
+} // end namespace Dumux
+
+#endif
diff --git a/test/flux/test_diffusioncontainers.cc b/test/flux/test_diffusioncontainers.cc
new file mode 100644
index 0000000000..44e96ebb45
--- /dev/null
+++ b/test/flux/test_diffusioncontainers.cc
@@ -0,0 +1,145 @@
+// -*- 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 FluxTests
+ * \brief This test makes sure that the diffusion coefficient containers
+ *        store the proper information, and are accessed correctly.
+ */
+
+#include <config.h>
+#include <dune/common/exceptions.hh>
+#include <dumux/common/parameters.hh>
+
+// include all diffusioncoefficient containers
+#include <dumux/flux/fickiandiffusioncoefficients.hh>
+#include <dumux/flux/maxwellstefandiffusioncoefficients.hh>
+
+template<int nP, int nC>
+struct TestFluidSystem
+{
+     static constexpr int numPhases = nP;
+     static constexpr int numComponents = nC;
+
+     int binaryDiffusionCoefficient(const int phaseIdx, const int compIIdx, const int compJIdx) const
+     { return phaseIdx*10000 + compIIdx*100 + compJIdx; }
+};
+
+int main()
+{
+    using namespace Dumux;
+    using Scalar = int;
+
+    const int numPhases = 3;
+    const int numComponents = 8;
+
+    TestFluidSystem<numPhases,numComponents> fluidSystem;
+
+    // temporary for debugging
+    // {
+    for (int phaseIdx = 0; phaseIdx < fluidSystem.numPhases; ++ phaseIdx)
+    {
+        int compIIdx = phaseIdx;
+        for (int compJIdx = 0; compJIdx < fluidSystem.numComponents; ++compJIdx)
+        {
+            std::cout << fluidSystem.binaryDiffusionCoefficient(phaseIdx,compIIdx,compJIdx) << "  ";
+        }
+        std::cout << "\n";
+    }
+    // }
+
+    const bool onlyTracers = false;
+    using DiffusionCoefficientsContainer = FickianDiffusionCoefficients<Scalar,
+                                                                        numPhases,
+                                                                        numComponents,
+                                                                        onlyTracers>;
+
+    auto getDiffusionCoefficient = [&](int phaseIdx, int compIIdx, int compJIdx)
+    {
+        return fluidSystem.binaryDiffusionCoefficient(phaseIdx,
+                                                      compIIdx,
+                                                      compJIdx);
+    };
+
+    DiffusionCoefficientsContainer mpncDiffCoeffs;
+    mpncDiffCoeffs.update(getDiffusionCoefficient);
+
+    for (int phaseIdx = 0; phaseIdx < fluidSystem.numPhases; ++ phaseIdx)
+    {
+        int compIIdx = phaseIdx;
+        std::cout << phaseIdx << " : \n";
+        for (int compJIdx = 0; compJIdx < fluidSystem.numComponents; ++compJIdx)
+        {
+            if (compIIdx == compJIdx)
+                continue;
+            int binaryCoeff = fluidSystem.binaryDiffusionCoefficient(phaseIdx, compIIdx, compJIdx);
+            int containerCoeff = mpncDiffCoeffs(phaseIdx, compIIdx, compJIdx);
+            if (containerCoeff != binaryCoeff)
+                DUNE_THROW(Dune::InvalidStateException, "Coefficients differ! " << "fluidSystem: " << binaryCoeff << " and "
+                                                                                << "container: " << containerCoeff <<
+                                                        ", for indicies (p,cI,cJ): (" << phaseIdx << ","
+                                                                                      << compIIdx << ","
+                                                                                      << compJIdx<<")" );
+            std::cout << " \n";
+        }
+    }
+
+    for (int phaseIdx = 0; phaseIdx < fluidSystem.numPhases; ++ phaseIdx)
+    {
+        int compJIdx = phaseIdx;
+        std::cout << phaseIdx << " : \n";
+        for (int compIIdx = 0; compIIdx < fluidSystem.numComponents; ++compIIdx)
+        {
+            if (compJIdx == compIIdx)
+                continue;
+            int binaryCoeff = fluidSystem.binaryDiffusionCoefficient(phaseIdx, compIIdx, compJIdx);
+            int containerCoeff = mpncDiffCoeffs(phaseIdx, compIIdx, compJIdx);
+            if (containerCoeff != binaryCoeff)
+                DUNE_THROW(Dune::InvalidStateException, "Coefficients differ! " << "fluidSystem: " << binaryCoeff << " and "
+                                                                                << "container: " << containerCoeff <<
+                                                        ", for indicies (p,cI,cJ): (" << phaseIdx << ","
+                                                                                      << compIIdx << ","
+                                                                                      << compJIdx<<")" );
+            std::cout << " \n";
+        }
+    }
+
+
+//     // test 1pnc
+//     template <class Scalar, int numComponents>
+//     class FickianDiffusionCoefficients<Scalar, 1, numComponents>
+
+//     // test 2p2c
+//     template <class Scalar>
+//     class FickianDiffusionCoefficients<Scalar, 2, 2>
+
+//     // test 3p2c
+//     template <class Scalar>
+//     class FickianDiffusionCoefficients<Scalar, 3, 2>
+
+//     // test tracer mpnc
+//     template <class Scalar, int numPhases, int numComponents>
+//     class FickianDiffusionCoefficients<Scalar, numPhases, numComponents, true>
+
+//     // test max stef
+//     template <class Scalar, int numPhases, int numComponents>
+//     class MaxwellStefanDiffusionCoefficients
+
+    return 0.0;
+}
-- 
GitLab