diff --git a/dumux/flux/fickiandiffusioncoefficients.hh b/dumux/flux/fickiandiffusioncoefficients.hh
index 976f4bb66970077a77b0a93e7abf9c4b0bf49bd8..11cb6386c2245fa80058bca69df6c44034c43f19 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 a543ada95089e2348ada2f5c32a08bfd0877e6af..9cec5e3ff866b7f0fa7cd8f49b96a4873354d95f 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 0000000000000000000000000000000000000000..0d26745087d0ea0223412d470141220580dc031d
--- /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 0000000000000000000000000000000000000000..95be775fc75994a2e949baaf91dc7ff2fc664a13
--- /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 0000000000000000000000000000000000000000..44e96ebb456a1a2a31063dd31dba071fac7714aa
--- /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;
+}