From a832b289c14300fafb4b5518b93eadcda505e577 Mon Sep 17 00:00:00 2001
From: Edward 'Ned' Coltman <nedc@flip.iws.uni-stuttgart.de>
Date: Wed, 25 Apr 2018 09:13:48 +0200
Subject: [PATCH] [freeflow][komega] First set of K-Omega classes

---
 dumux/freeflow/rans/twoeq/CMakeLists.txt      |   1 +
 .../freeflow/rans/twoeq/komega/CMakeLists.txt |  12 +
 .../rans/twoeq/komega/fluxvariables.hh        |  49 ++++
 dumux/freeflow/rans/twoeq/komega/indices.hh   |  54 ++++
 .../rans/twoeq/komega/localresidual.hh        |  51 ++++
 dumux/freeflow/rans/twoeq/komega/model.hh     | 135 ++++++++++
 dumux/freeflow/rans/twoeq/komega/models.hh    |  49 ++++
 dumux/freeflow/rans/twoeq/komega/problem.hh   | 171 ++++++++++++
 .../twoeq/komega/staggered/CMakeLists.txt     |   4 +
 .../twoeq/komega/staggered/fluxvariables.hh   | 176 +++++++++++++
 .../twoeq/komega/staggered/localresidual.hh   | 127 +++++++++
 .../rans/twoeq/komega/volumevariables.hh      | 246 ++++++++++++++++++
 .../rans/twoeq/komega/vtkoutputfields.hh      |  67 +++++
 13 files changed, 1142 insertions(+)
 create mode 100644 dumux/freeflow/rans/twoeq/komega/CMakeLists.txt
 create mode 100644 dumux/freeflow/rans/twoeq/komega/fluxvariables.hh
 create mode 100644 dumux/freeflow/rans/twoeq/komega/indices.hh
 create mode 100644 dumux/freeflow/rans/twoeq/komega/localresidual.hh
 create mode 100644 dumux/freeflow/rans/twoeq/komega/model.hh
 create mode 100644 dumux/freeflow/rans/twoeq/komega/models.hh
 create mode 100644 dumux/freeflow/rans/twoeq/komega/problem.hh
 create mode 100644 dumux/freeflow/rans/twoeq/komega/staggered/CMakeLists.txt
 create mode 100644 dumux/freeflow/rans/twoeq/komega/staggered/fluxvariables.hh
 create mode 100644 dumux/freeflow/rans/twoeq/komega/staggered/localresidual.hh
 create mode 100644 dumux/freeflow/rans/twoeq/komega/volumevariables.hh
 create mode 100644 dumux/freeflow/rans/twoeq/komega/vtkoutputfields.hh

diff --git a/dumux/freeflow/rans/twoeq/CMakeLists.txt b/dumux/freeflow/rans/twoeq/CMakeLists.txt
index 3778f272e5..cb96f788ce 100644
--- a/dumux/freeflow/rans/twoeq/CMakeLists.txt
+++ b/dumux/freeflow/rans/twoeq/CMakeLists.txt
@@ -1,2 +1,3 @@
 add_subdirectory(kepsilon)
+add_subdirectory(komega)
 add_subdirectory(lowrekepsilon)
diff --git a/dumux/freeflow/rans/twoeq/komega/CMakeLists.txt b/dumux/freeflow/rans/twoeq/komega/CMakeLists.txt
new file mode 100644
index 0000000000..e6ef83c51b
--- /dev/null
+++ b/dumux/freeflow/rans/twoeq/komega/CMakeLists.txt
@@ -0,0 +1,12 @@
+add_subdirectory(staggered)
+
+install(FILES
+fluxvariables.hh
+indices.hh
+localresidual.hh
+model.hh
+models.hh
+problem.hh
+volumevariables.hh
+vtkoutputfields.hh
+DESTINATION ${CMAKE_INSTALL_INCLUDEDIR}/dumux/freeflow/rans/twoeq/komega)
diff --git a/dumux/freeflow/rans/twoeq/komega/fluxvariables.hh b/dumux/freeflow/rans/twoeq/komega/fluxvariables.hh
new file mode 100644
index 0000000000..70c3d6af60
--- /dev/null
+++ b/dumux/freeflow/rans/twoeq/komega/fluxvariables.hh
@@ -0,0 +1,49 @@
+// -*- 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 2 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 KOmegaModel
+  * \copydoc Dumux::KOmegaFluxVariables
+  */
+#ifndef DUMUX_KOMEGA_FLUXVARIABLES_HH
+#define DUMUX_KOMEGA_FLUXVARIABLES_HH
+
+#include <dumux/common/properties.hh>
+#include <dumux/freeflow/rans/twoeq/komega/staggered/fluxvariables.hh>
+
+namespace Dumux
+{
+
+// forward declaration
+template<class TypeTag, DiscretizationMethod discMethod>
+class KOmegaFluxVariablesImpl;
+
+/*!
+ * \ingroup KOmegaModel
+ * \brief The flux variables class for the k-omega model.
+          This is a convenience alias for that actual,
+          discretization-specific flux variables.
+ * \note  Not all specializations are currently implemented
+ */
+template<class TypeTag>
+using KOmegaFluxVariables = KOmegaFluxVariablesImpl<TypeTag, GET_PROP_TYPE(TypeTag, FVGridGeometry)::discMethod>;
+
+} // end namespace
+
+#endif
diff --git a/dumux/freeflow/rans/twoeq/komega/indices.hh b/dumux/freeflow/rans/twoeq/komega/indices.hh
new file mode 100644
index 0000000000..9d7284fde8
--- /dev/null
+++ b/dumux/freeflow/rans/twoeq/komega/indices.hh
@@ -0,0 +1,54 @@
+// -*- 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 2 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 KOmegaModel
+ * \copydoc Dumux::KOmegaIndices
+ */
+#ifndef DUMUX_KOMEGA_INDICES_HH
+#define DUMUX_KOMEGA_INDICES_HH
+
+#include <dumux/freeflow/navierstokes/indices.hh>
+
+namespace Dumux {
+
+// \{
+/*!
+ * \ingroup KOmegaModel
+ * \brief The common indices for the isothermal KOmega model.
+ *
+ * \tparam dimension The dimension of the problem
+ */
+template <int dimension>
+struct KOmegaIndices : public NavierStokesIndices<dimension>
+{
+private:
+    using ParentType = NavierStokesIndices<dimension>;
+
+public:
+    static constexpr auto turbulentKineticEnergyEqIdx = dimension + numComponents;
+    static constexpr auto turbulentKineticEnergyIdx = turbulentKineticEnergyEqIdx;
+    static constexpr auto dissipationEqIdx = turbulentKineticEnergyEqIdx + 1;
+    static constexpr auto dissipationIdx = dissipationEqIdx;
+};
+
+// \}
+} // end namespace
+
+#endif
diff --git a/dumux/freeflow/rans/twoeq/komega/localresidual.hh b/dumux/freeflow/rans/twoeq/komega/localresidual.hh
new file mode 100644
index 0000000000..61764cc545
--- /dev/null
+++ b/dumux/freeflow/rans/twoeq/komega/localresidual.hh
@@ -0,0 +1,51 @@
+// -*- 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 2 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 KOmega TwoEq Model
+ * \copydoc Dumux::KOmegaKOmegaResidual
+ */
+#ifndef DUMUX_KOMEGA_LOCAL_RESIDUAL_HH
+#define DUMUX_KOMEGA_LOCAL_RESIDUAL_HH
+
+#include <dumux/common/properties.hh>
+#include <dumux/discretization/methods.hh>
+#include <dumux/freeflow/navierstokes/localresidual.hh>
+#include <dumux/freeflow/rans/twoeq/komega/staggered/localresidual.hh>
+
+namespace Dumux
+{
+
+// forward declaration
+template<class TypeTag, DiscretizationMethod discMethod>
+class KOmegaResidualImpl;
+
+/*!
+ * \ingroup KOmega TwoEq Model
+ * \brief The local residual class for the k-omega model.
+          This is a convenience alias for the actual,
+          discretization-specific local residual.
+ * \note  Not all specializations are currently implemented
+ */
+template<class TypeTag>
+using KOmegaResidual = KOmegaResidualImpl<TypeTag, GET_PROP_TYPE(TypeTag, FVGridGeometry)::discMethod>;
+
+}
+
+#endif
diff --git a/dumux/freeflow/rans/twoeq/komega/model.hh b/dumux/freeflow/rans/twoeq/komega/model.hh
new file mode 100644
index 0000000000..c84ccf8ec2
--- /dev/null
+++ b/dumux/freeflow/rans/twoeq/komega/model.hh
@@ -0,0 +1,135 @@
+// -*- 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 2 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 KOmegaModel
+ *
+ * \brief A single-phase, isothermal k-omega 2-Eq. model
+ *
+ * \copydoc RANSModel
+ *
+ * These models calculate the eddy viscosity with two additional PDEs,
+ * one for the turbulentKineticEnergy (k) and a second for the dissipation (omega).
+ */
+
+#ifndef DUMUX_KOMEGA_MODEL_HH
+#define DUMUX_KOMEGA_MODEL_HH
+
+#include <dumux/common/properties.hh>
+#include <dumux/freeflow/properties.hh>
+#include <dumux/freeflow/rans/model.hh>
+
+#include "fluxvariables.hh"
+#include "indices.hh"
+#include "localresidual.hh"
+#include "volumevariables.hh"
+#include "vtkoutputfields.hh"
+
+namespace Dumux
+{
+namespace Properties {
+
+/*!
+ *\ingroup KOmegaModel
+ * \brief Traits for the k-omega model
+ */
+template<int dimension>
+struct KOmegaModelTraits
+{
+    //! The dimension of the model
+    static constexpr int dim() { return dimension; }
+
+    //! There are as many momentum balance equations as dimensions,
+    //! one mass balance equation and two turbulent transport equations
+    static constexpr int numEq() { return dim()+1+2; }
+
+    //! The number of phases is always 1
+    static constexpr int numPhases() { return 1; }
+
+    //! The number of components
+    static constexpr int numComponents() { return 1; }
+
+    //! Enable advection
+    static constexpr bool enableAdvection() { return true; }
+
+    //! The one-phase model has no molecular diffusion
+    static constexpr bool enableMolecularDiffusion() { return true; }
+
+    //! The model is isothermal
+    static constexpr bool enableEnergyBalance() { return false; }
+
+    //! The indices
+    using Indices = KOmegaIndices<dim()>;
+
+    //! The model includes a limiter to the production term
+    static constexpr bool enableKOmegaProductionLimiter() { return true; }
+};
+
+///////////////////////////////////////////////////////////////////////////
+// default property values for the isothermal k-omega single phase model
+///////////////////////////////////////////////////////////////////////////
+
+//! The type tag for the single-phase, isothermal k-omega model
+NEW_TYPE_TAG(KOmega, INHERITS_FROM(RANS));
+
+//!< states some specifics of the isothermal k-omega model
+SET_PROP(KOmega, ModelTraits)
+{
+private:
+    using GridView = typename GET_PROP_TYPE(TypeTag, FVGridGeometry)::GridView;
+    static constexpr int dim = GridView::dimension;
+public:
+    using type = KOmegaModelTraits<dim>;
+};
+
+//! The flux variables
+SET_TYPE_PROP(KOmega, FluxVariables, KOmegaFluxVariables<TypeTag>);
+
+//! The local residual
+SET_TYPE_PROP(KOmega, LocalResidual, KOmegaResidual<TypeTag>);
+
+//! The volume variables
+SET_TYPE_PROP(KOmega, VolumeVariables, KOmegaVolumeVariables<TypeTag>);
+
+//! The specific vtk output fields
+SET_PROP(KOmega, VtkOutputFields)
+{
+private:
+    using FVGridGeometry = typename GET_PROP_TYPE(TypeTag, FVGridGeometry);
+public:
+    using type = KOmegaVtkOutputFields<FVGridGeometry>;
+};
+
+///////////////////////////////////////////////////////////////////////////
+// default property values for the non-isothermal k-omega single phase model
+///////////////////////////////////////////////////////////////////////////
+
+
+//! The type tag for the single-phase, non-isothermal k-omega 2-Eq. model
+NEW_TYPE_TAG(KOmegaNI, INHERITS_FROM(RANSNI));
+
+//! The volume variables
+SET_TYPE_PROP(KOmegaNI, VolumeVariables, KOmegaVolumeVariables<TypeTag>);
+
+// \}
+}
+
+} // end namespace
+
+#endif // DUMUX_KOMEGA_MODEL_HH
diff --git a/dumux/freeflow/rans/twoeq/komega/models.hh b/dumux/freeflow/rans/twoeq/komega/models.hh
new file mode 100644
index 0000000000..c2e04ec8dc
--- /dev/null
+++ b/dumux/freeflow/rans/twoeq/komega/models.hh
@@ -0,0 +1,49 @@
+// -*- 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 2 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 KOmegaModel
+ * \copydoc Dumux::KOmegaIndices
+ */
+#ifndef DUMUX_KOMEGA_MODELS_HH
+#define DUMUX_KOMEGA_MODELS_HH
+
+namespace Dumux {
+
+/*!
+ * \ingroup KOmegaModel
+ * \brief The available eddy viscosity models
+ *
+ * \todo update
+ * The following models are available:
+ *  -# original model developed in \cite Wilcox08
+ *  -# original model developed in \cite Wilcox88
+ *
+ * A good overview and additional models are given in \cite Patel1985a and online at NASA Turbulence Reference (https://turbmodels.larc.nasa.gov/)
+ */
+class KOmegaModels
+{
+public:
+    static constexpr int wilcox08 = 0;
+    static constexpr int wilcox88 = 1;
+};
+
+} // end namespace Dumux
+
+#endif
diff --git a/dumux/freeflow/rans/twoeq/komega/problem.hh b/dumux/freeflow/rans/twoeq/komega/problem.hh
new file mode 100644
index 0000000000..f31bf89487
--- /dev/null
+++ b/dumux/freeflow/rans/twoeq/komega/problem.hh
@@ -0,0 +1,171 @@
+// -*- 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 2 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 KOmega TwoEq Model
+ * \copydoc Dumux::KomegaProblem
+ */
+#ifndef DUMUX_KOMEGA_PROBLEM_HH
+#define DUMUX_KOMEGA_PROBLEM_HH
+
+#include <dumux/common/properties.hh>
+#include <dumux/common/staggeredfvproblem.hh>
+#include <dumux/discretization/localview.hh>
+#include <dumux/discretization/staggered/elementsolution.hh>
+#include <dumux/discretization/methods.hh>
+#include <dumux/freeflow/rans/problem.hh>
+
+#include "model.hh"
+#include "models.hh"
+
+namespace Dumux
+{
+
+/*!
+ * \ingroup KOmega TwoEq Model
+ * \brief K-Omega turbulence model problem base class.
+ *
+ * This implements the 2-equation k-omega turbulence model developed in Wilcox08 and Wilcox88
+ */
+template<class TypeTag>
+class KOmegaProblem : public RANSProblem<TypeTag>
+{
+    using ParentType = RANSProblem<TypeTag>;
+    using Implementation = typename GET_PROP_TYPE(TypeTag, Problem);
+
+    using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar);
+
+    using FVGridGeometry = typename GET_PROP_TYPE(TypeTag, FVGridGeometry);
+    using FVElementGeometry = typename GET_PROP_TYPE(TypeTag, FVGridGeometry)::LocalView;
+    using SubControlVolumeFace = typename FVElementGeometry::SubControlVolumeFace;
+
+    using VolumeVariables = typename GET_PROP_TYPE(TypeTag, VolumeVariables);
+    using SolutionVector = typename GET_PROP_TYPE(TypeTag, SolutionVector);
+    using PrimaryVariables = typename GET_PROP_TYPE(TypeTag, PrimaryVariables);
+    using CellCenterPrimaryVariables = typename GET_PROP_TYPE(TypeTag, CellCenterPrimaryVariables);
+    using Indices = typename GET_PROP_TYPE(TypeTag, ModelTraits)::Indices;
+
+    using GridView = typename GET_PROP_TYPE(TypeTag, GridView);
+    using Grid = typename GridView::Grid;
+    enum {
+        dim = Grid::dimension,
+      };
+    using GlobalPosition = Dune::FieldVector<Scalar, dim>;
+    using DimVector = Dune::FieldVector<Scalar, dim>;
+    using DimMatrix = Dune::FieldMatrix<Scalar, dim, dim>;
+
+public:
+    KOmegaProblem(std::shared_ptr<const FVGridGeometry> fvGridGeometry) : ParentType(fvGridGeometry)
+    {
+        kOmegaModel_ = getParamFromGroup<int>(GET_PROP_VALUE(TypeTag, ModelParameterGroup), "KOmega.Model", 0);
+        if (kOmegaModel_ != KOmegaModels::wilcox08
+            && kOmegaModel_ != KOmegaModels::wilcox88)
+        {
+            DUNE_THROW(Dune::NotImplemented, "This k-omega model is not implemented: " << kOmegaModel_);
+        }
+        useStoredEddyViscosity_ = getParamFromGroup<bool>(GET_PROP_VALUE(TypeTag, ModelParameterGroup), "RANS.UseStoredEddyViscosity", true);
+    }
+
+    /*!
+     * \brief Correct size of the static (solution independent) wall variables
+     */
+    void updateStaticWallProperties()
+    {
+        ParentType::updateStaticWallProperties();
+
+        // update size and initial values of the global vectors
+        storedDissipation_.resize(this->fvGridGeometry().elementMapper().size(), 0.0);
+        storedDissipationGradient_.resize(this->fvGridGeometry().elementMapper().size(), DimMatrix(0.0));
+        storedDynamicEddyViscosity_.resize(this->fvGridGeometry().elementMapper().size(), 0.0);
+        storedTurbulentKineticEnergy_.resize(this->fvGridGeometry().elementMapper().size(), 0.0);
+        storedTurbulentKineticEnergyGradient_.resize(this->fvGridGeometry().elementMapper().size(), 0.0);
+    }
+
+    /*!
+     * \brief Update the dynamic (solution dependent) relations to the walls
+     *
+     * \param curSol The solution vector.
+     */
+    void updateDynamicWallProperties(const SolutionVector& curSol)
+    {
+        ParentType::updateDynamicWallProperties(curSol);
+
+        for (const auto& element : elements(this->fvGridGeometry().gridView()))
+        {
+            unsigned int elementID = this->fvGridGeometry().elementMapper().index(element);
+
+            auto fvGeometry = localView(this->fvGridGeometry());
+            fvGeometry.bindElement(element);
+            for (auto&& scv : scvs(fvGeometry))
+            {
+                const int dofIdx = scv.dofIndex();
+                const auto& cellCenterPriVars = curSol[FVGridGeometry::cellCenterIdx()][dofIdx];
+                PrimaryVariables priVars = makePriVarsFromCellCenterPriVars<PrimaryVariables>(cellCenterPriVars);
+                auto elemSol = elementSolution<typename FVGridGeometry::LocalView>(std::move(priVars));
+                // NOTE: first update the turbulence quantities
+                storedDissipation_[elementID] = elemSol[0][Indices::dissipationEqIdx];
+                storedTurbulentKineticEnergy_[elementID] = elemSol[0][Indices::turbulentKineticEnergyEqIdx];
+                // NOTE: then update the volVars
+                VolumeVariables volVars;
+                volVars.update(elemSol, asImp_(), element, scv);
+                volVars.calculateEddyViscosity();
+                storedKinematicEddyViscosity_[elementID] = volVars.kinematicEddyViscosity();
+            }
+        }
+    }
+
+    //! \brief Returns the \$f \beta_{\omega} \$f constant
+    const Scalar betaOmega() const
+    {
+          if(kOmegaModel_ == KOmegaModels::wilcox88)
+              return 0.0750;
+          else if (kOmegaModel_ == KOmegaModels::wilcox08)
+              return 0.0708;
+          else
+              DUNE_THROW(Dune::NotImplemented, "This Model is not implemented.");
+    }
+
+    /*!
+     * \brief Returns the index of the k-omega model applied
+     */
+    const int kOmegaModel() const
+    { return kOmegaModel_; }
+
+public:
+    std::vector<Scalar> storedDissipation_;
+    std::vector<Scalar> storedDynamicEddyViscosity_;
+    std::vector<Scalar> storedTurbulentKineticEnergy_;
+    std::vector<DimMatrix> storedDissipationGradient_;
+    std::vector<DimMatrix> storedTurbulentKineticEnergyGradient_;
+    int kOmegaModel_;
+    bool useStoredEddyViscosity_;
+
+private:
+    //! Returns the implementation of the problem (i.e. static polymorphism)
+    Implementation &asImp_()
+    { return *static_cast<Implementation *>(this); }
+
+    //! \copydoc asImp_()
+    const Implementation &asImp_() const
+    { return *static_cast<const Implementation *>(this); }
+};
+
+}
+
+#endif
diff --git a/dumux/freeflow/rans/twoeq/komega/staggered/CMakeLists.txt b/dumux/freeflow/rans/twoeq/komega/staggered/CMakeLists.txt
new file mode 100644
index 0000000000..e700e78b28
--- /dev/null
+++ b/dumux/freeflow/rans/twoeq/komega/staggered/CMakeLists.txt
@@ -0,0 +1,4 @@
+install(FILES
+fluxvariables.hh
+localresidual.hh
+DESTINATION ${CMAKE_INSTALL_INCLUDEDIR}/dumux/freeflow/rans/twoeq/komega/staggered)
diff --git a/dumux/freeflow/rans/twoeq/komega/staggered/fluxvariables.hh b/dumux/freeflow/rans/twoeq/komega/staggered/fluxvariables.hh
new file mode 100644
index 0000000000..49f8f581c9
--- /dev/null
+++ b/dumux/freeflow/rans/twoeq/komega/staggered/fluxvariables.hh
@@ -0,0 +1,176 @@
+// -*- 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 2 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 KOmegaModel
+  * \copydoc Dumux::KOmegaFluxVariablesImpl
+  */
+#ifndef DUMUX_KOMEGA_STAGGERED_FLUXVARIABLES_HH
+#define DUMUX_KOMEGA_STAGGERED_FLUXVARIABLES_HH
+
+#include <numeric>
+#include <dumux/common/properties.hh>
+#include <dumux/discretization/fluxvariablesbase.hh>
+#include <dumux/discretization/methods.hh>
+#include <dumux/freeflow/navierstokes/fluxvariables.hh>
+#include <dumux/freeflow/rans/twoeq/komega/fluxvariables.hh>
+#include <dumux/freeflow/rans/twoeq/komega/models.hh>
+
+namespace Dumux
+{
+
+// forward declaration
+template<class TypeTag, DiscretizationMethod discMethod>
+class KOmegaFluxVariablesImpl;
+
+/*!
+  * \ingroup KOmegaModel
+  * \brief The flux variables class for the k-omega model using the staggered grid discretization.
+ */
+template<class TypeTag>
+class KOmegaFluxVariablesImpl<TypeTag, DiscretizationMethod::staggered>
+: public NavierStokesFluxVariables<TypeTag>
+{
+    using ParentType = NavierStokesFluxVariables<TypeTag>;
+    using GridVariables = typename GET_PROP_TYPE(TypeTag, GridVariables);
+    using GridVolumeVariables = typename GridVariables::GridVolumeVariables;
+    using ElementVolumeVariables = typename GridVolumeVariables::LocalView;
+    using VolumeVariables = typename GridVolumeVariables::VolumeVariables;
+    using GridFluxVariablesCache = typename GridVariables::GridFluxVariablesCache;
+    using FluxVariablesCache = typename GridFluxVariablesCache::FluxVariablesCache;
+    using GridFaceVariables = typename GridVariables::GridFaceVariables;
+    using ElementFaceVariables = typename GridFaceVariables::LocalView;
+    using FaceVariables = typename GridFaceVariables::FaceVariables;
+
+    using ModelTraits = typename GET_PROP_TYPE(TypeTag, ModelTraits);
+    using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar);
+    using GridView = typename GET_PROP_TYPE(TypeTag, GridView);
+    using Problem = typename GET_PROP_TYPE(TypeTag, Problem);
+    using Element = typename GridView::template Codim<0>::Entity;
+    using FVElementGeometry = typename GET_PROP_TYPE(TypeTag, FVGridGeometry)::LocalView;
+    using Indices = typename GET_PROP_TYPE(TypeTag, ModelTraits)::Indices;
+    using SubControlVolumeFace = typename FVElementGeometry::SubControlVolumeFace;
+    using CellCenterPrimaryVariables = typename GET_PROP_TYPE(TypeTag, CellCenterPrimaryVariables);
+    enum {
+        turbulentKineticEnergyEqIdx = Indices::turbulentKineticEnergyEqIdx,
+        dissipationEqIdx = Indices::dissipationEqIdx,
+    };
+
+public:
+
+    /*!
+    * \brief Computes the flux for the cell center residual.
+    */
+    CellCenterPrimaryVariables computeFluxForCellCenter(const Problem& problem,
+                                                        const Element &element,
+                                                        const FVElementGeometry& fvGeometry,
+                                                        const ElementVolumeVariables& elemVolVars,
+                                                        const ElementFaceVariables& elemFaceVars,
+                                                        const SubControlVolumeFace &scvf,
+                                                        const FluxVariablesCache& fluxVarsCache)
+    {
+        CellCenterPrimaryVariables flux = ParentType::computeFluxForCellCenter(problem, element, fvGeometry,
+                                                                               elemVolVars, elemFaceVars, scvf, fluxVarsCache);
+
+        // calculate advective flux
+        const auto bcTypes = problem.boundaryTypesAtPos(scvf.center());
+        const bool isOutflowK = scvf.boundary() && bcTypes.isOutflow(turbulentKineticEnergyEqIdx);
+        const bool isOutflowOmega = scvf.boundary() && bcTypes.isOutflow(dissipationEqIdx);
+        auto upwindTermK = [](const auto& volVars)
+        {
+            return volVars.turbulentKineticEnergy();
+        };
+        auto upwindTermOmega = [](const auto& volVars)
+        {
+            return volVars.dissipation();
+        };
+
+        flux[turbulentKineticEnergyEqIdx - ModelTraits::dim()]
+            = ParentType::advectiveFluxForCellCenter(elemVolVars, elemFaceVars, scvf, upwindTermK, isOutflowK);
+        flux[dissipationEqIdx - ModelTraits::dim()]
+            = ParentType::advectiveFluxForCellCenter(elemVolVars, elemFaceVars, scvf, upwindTermOmega, isOutflowOmega);
+        Dune::dverb << " k_adv " << ParentType::advectiveFluxForCellCenter(elemVolVars, elemFaceVars, scvf, upwindTermK, isOutflowK);
+        Dune::dverb << " w_adv " << ParentType::advectiveFluxForCellCenter(elemVolVars, elemFaceVars, scvf, upwindTermOmega, isOutflowOmega);
+
+        // calculate diffusive flux
+        const auto& insideScv = fvGeometry.scv(scvf.insideScvIdx());
+        const auto& outsideScv = fvGeometry.scv(scvf.outsideScvIdx());
+        const auto& insideVolVars = elemVolVars[scvf.insideScvIdx()];
+        const auto& outsideVolVars = elemVolVars[scvf.outsideScvIdx()];
+
+        // effective diffusion coefficients
+        Scalar insideCoeff_k = insideVolVars.kinematicViscosity()
+                                + ( insideVolVars.sigmaK() * insideVolVars.turbulentKineticEnergy() / insideVolVars.dissipation() );
+        Scalar outsideCoeff_k = outsideVolVars.kinematicViscosity()
+                                + ( outsideVolVars.sigmaK() * outsideVolVars.turbulentKineticEnergy() / outsideVolVars.dissipation() );
+        Scalar insideCoeff_w = insideVolVars.kinematicViscosity()
+                                + ( insideVolVars.sigmaOmega() * insideVolVars.turbulentKineticEnergy() / insideVolVars.dissipation() );
+        Scalar outsideCoeff_w = outsideVolVars.kinematicViscosity()
+                                + ( outsideVolVars.sigmaOmega() * outsideVolVars.turbulentKineticEnergy() / outsideVolVars.dissipation() );
+
+        // scale by extrusion factor
+        insideCoeff_k *= insideVolVars.extrusionFactor();
+        outsideCoeff_k *= outsideVolVars.extrusionFactor();
+        insideCoeff_w *= insideVolVars.extrusionFactor();
+        outsideCoeff_w *= outsideVolVars.extrusionFactor();
+
+        // average and distance
+        Scalar coeff_k = (insideCoeff_k - outsideCoeff_k) * 0.5;
+        Scalar coeff_w = (insideCoeff_w - outsideCoeff_w) * 0.5;
+        Scalar distance = 0.0;
+
+        // adapt boundary handling
+        if (scvf.boundary())
+        {
+            distance = (insideScv.dofPosition() - scvf.ipGlobal()).two_norm();
+            coeff_k = insideCoeff_k;
+            coeff_w = insideCoeff_w;
+        }
+        else
+        {
+            distance = (outsideScv.dofPosition() - insideScv.dofPosition()).two_norm();
+        }
+
+        if (!isOutflowK)
+        {
+            flux[turbulentKineticEnergyEqIdx - ModelTraits::dim()]
+                += coeff_k / distance
+                   * (insideVolVars.turbulentKineticEnergy() - outsideVolVars.turbulentKineticEnergy())
+                   * scvf.area();
+            Dune::dverb << " k_diff " << coeff_k / distance
+                                                 * (insideVolVars.turbulentKineticEnergy() - outsideVolVars.turbulentKineticEnergy())
+                                                 * scvf.area();
+        }
+        if (!isOutflowOmega)
+        {
+            flux[dissipationEqIdx - ModelTraits::dim()]
+                += coeff_w / distance
+                   * (insideVolVars.dissipation() - outsideVolVars.dissipation())
+                   * scvf.area();
+            Dune::dverb << " w_diff " << coeff_w / distance
+                                                 * (insideVolVars.dissipation() - outsideVolVars.dissipation())
+                                                 * scvf.area();
+        }
+        return flux;
+    }
+};
+
+} // end namespace
+
+#endif
diff --git a/dumux/freeflow/rans/twoeq/komega/staggered/localresidual.hh b/dumux/freeflow/rans/twoeq/komega/staggered/localresidual.hh
new file mode 100644
index 0000000000..1f9aee5950
--- /dev/null
+++ b/dumux/freeflow/rans/twoeq/komega/staggered/localresidual.hh
@@ -0,0 +1,127 @@
+// -*- 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 2 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 KOmegaModel
+  * \copydoc Dumux::KOmegaResidualImpl
+  */
+#ifndef DUMUX_STAGGERED_KOMEGA_LOCAL_RESIDUAL_HH
+#define DUMUX_STAGGERED_KOMEGA_LOCAL_RESIDUAL_HH
+
+#include <dune/common/hybridutilities.hh>
+#include <dumux/common/properties.hh>
+#include <dumux/discretization/methods.hh>
+#include <dumux/freeflow/navierstokes/localresidual.hh>
+
+namespace Dumux {
+
+// forward declaration
+template<class TypeTag,  DiscretizationMethod discMethod>
+class KOmegaResidualImpl;
+
+/*!
+  * \ingroup KOmegaModel
+  * \brief Element-wise calculation of the residual for k-omega models using the staggered discretization
+ */
+template<class TypeTag>
+class KOmegaResidualImpl<TypeTag, DiscretizationMethod::staggered>
+: public NavierStokesResidual<TypeTag>
+{
+    using ParentType = NavierStokesResidual<TypeTag>;
+    friend class StaggeredLocalResidual<TypeTag>;
+    friend ParentType;
+
+    using GridVariables = typename GET_PROP_TYPE(TypeTag, GridVariables);
+
+    using GridVolumeVariables = typename GridVariables::GridVolumeVariables;
+    using ElementVolumeVariables = typename GridVolumeVariables::LocalView;
+    using VolumeVariables = typename GridVolumeVariables::VolumeVariables;
+
+    using GridFluxVariablesCache = typename GridVariables::GridFluxVariablesCache;
+    using ElementFluxVariablesCache = typename GridFluxVariablesCache::LocalView;
+    using FluxVariables = typename GET_PROP_TYPE(TypeTag, FluxVariables);
+
+    using GridFaceVariables = typename GridVariables::GridFaceVariables;
+    using ElementFaceVariables = typename GridFaceVariables::LocalView;
+
+    using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar);
+    using Problem = typename GET_PROP_TYPE(TypeTag, Problem);
+    using GridView = typename GET_PROP_TYPE(TypeTag, GridView);
+    using Element = typename GridView::template Codim<0>::Entity;
+    using FVElementGeometry = typename GET_PROP_TYPE(TypeTag, FVGridGeometry)::LocalView;
+    using SubControlVolume = typename FVElementGeometry::SubControlVolume;
+    using CellCenterPrimaryVariables = typename GET_PROP_TYPE(TypeTag, CellCenterPrimaryVariables);
+    using Indices = typename GET_PROP_TYPE(TypeTag, ModelTraits)::Indices;
+    using BoundaryTypes = typename GET_PROP_TYPE(TypeTag, BoundaryTypes);
+    using ModelTraits = typename GET_PROP_TYPE(TypeTag, ModelTraits);
+    using CellCenterResidual = typename GET_PROP_TYPE(TypeTag, CellCenterPrimaryVariables);
+
+public:
+    using ParentType::ParentType;
+
+    //! Evaluate fluxes entering or leaving the cell center control volume.
+    CellCenterPrimaryVariables computeStorageForCellCenter(const Problem& problem,
+                                                           const SubControlVolume& scv,
+                                                           const VolumeVariables& volVars) const
+    {
+        CellCenterPrimaryVariables storage = ParentType::computeStorageForCellCenter(problem, scv, volVars);
+
+        static constexpr int turbulentKineticEnergyEqIdx = Indices::turbulentKineticEnergyEqIdx - ModelTraits::dim();
+        static constexpr int dissipationEqIdx = Indices::dissipationEqIdx - ModelTraits::dim();
+        storage[turbulentKineticEnergyEqIdx] = volVars.turbulentKineticEnergy();
+        storage[dissipationEqIdx] = volVars.dissipation();
+
+        return storage;
+    }
+
+    CellCenterPrimaryVariables computeSourceForCellCenter(const Problem& problem,
+                                                          const Element &element,
+                                                          const FVElementGeometry& fvGeometry,
+                                                          const ElementVolumeVariables& elemVolVars,
+                                                          const ElementFaceVariables& elemFaceVars,
+                                                          const SubControlVolume &scv) const
+    {
+        CellCenterPrimaryVariables source = ParentType::computeSourceForCellCenter(problem, element, fvGeometry,
+                                                                                   elemVolVars, elemFaceVars, scv);
+
+        const auto& volVars = elemVolVars[scv];
+
+        static constexpr int turbulentKineticEnergyEqIdx = Indices::turbulentKineticEnergyEqIdx - ModelTraits::dim();
+        static constexpr int dissipationEqIdx = Indices::dissipationEqIdx - ModelTraits::dim();
+
+    // production
+        Scalar productionTerm = 2.0 * volVars.kinematicEddyViscosity()* volVars.stressTensorScalarProduct();
+        if (ModelTraits::enableKOmegaProductionLimiter())
+        {
+          Scalar productionAlternative = 20.0 * volVars.betaK() * volVars.turbulentKineticEnergy() * volVars.dissipation();
+          productionTerm = std::min(productionTerm, productionAlternative);
+        }
+        source[turbulentKineticEnergyEqIdx] += productionTerm;
+        source[dissipationEqIdx] += volVars.alpha() * ( volVars.dissipation() / volVars.turbulentKineticEnergy() ) * productionTerm;
+
+    // destruction
+        source[turbulentKineticEnergyEqIdx] -= volVars.betaK() * volVars.turbulentKineticEnergy() * volVars.dissipation();
+        source[dissipationEqIdx] -=  problem.betaOmega() * volVars.dissipation() * volVars.dissipation();
+
+        return source;
+    }
+};
+} // end namespace Dumux
+
+#endif
diff --git a/dumux/freeflow/rans/twoeq/komega/volumevariables.hh b/dumux/freeflow/rans/twoeq/komega/volumevariables.hh
new file mode 100644
index 0000000000..c0a04542ef
--- /dev/null
+++ b/dumux/freeflow/rans/twoeq/komega/volumevariables.hh
@@ -0,0 +1,246 @@
+// -*- 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 2 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 KOmegaEqModel
+ *
+ * \copydoc Dumux::KOmegaVolumeVariables
+ */
+#ifndef DUMUX_KOMEGA_VOLUME_VARIABLES_HH
+#define DUMUX_KOMEGA_VOLUME_VARIABLES_HH
+
+#include <dumux/common/properties.hh>
+#include <dumux/common/parameters.hh>
+#include <dumux/material/fluidstates/immiscible.hh>
+#include <dumux/freeflow/rans/volumevariables.hh>
+
+#include "models.hh"
+
+namespace Dumux
+{
+<<<<<<< HEAD
+
+// forward declaration
+template <class TypeTag, bool enableEnergyBalance>
+class KOmegaVolumeVariablesImplementation;
+
+/*!
+ * \ingroup KOmegaEqModel
+ * \brief Volume variables for the single-phase k-omega 2-Eq model
+ *        The class is specialized for isothermal and non-isothermal models.
+ */
+template <class TypeTag>
+using KOmegaVolumeVariables = KOmegaVolumeVariablesImplementation<TypeTag, GET_PROP_TYPE(TypeTag, ModelTraits)::enableEnergyBalance()>;
+=======
+>>>>>>> f6e896abb... First set of K-omega classes
+
+/*!
+ * \ingroup KOmegaEqModel
+ * \brief Volume variables for the isothermal single-phase k-omega 2-Eq model.
+ */
+template <class TypeTag>
+class KOmegaVolumeVariablesImplementation<TypeTag, false>
+: virtual public RANSVolumeVariablesImplementation<TypeTag, false>
+{
+    using ParentType = RANSVolumeVariablesImplementation<TypeTag, false>;
+    using Implementation = typename GET_PROP_TYPE(TypeTag, VolumeVariables);
+    using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar);
+    using Problem = typename GET_PROP_TYPE(TypeTag, Problem);
+    using FVElementGeometry = typename GET_PROP_TYPE(TypeTag, FVGridGeometry)::LocalView;
+    using SubControlVolume = typename FVElementGeometry::SubControlVolume;
+    using Indices = typename GET_PROP_TYPE(TypeTag, ModelTraits)::Indices;
+    using GridView = typename GET_PROP_TYPE(TypeTag, GridView);
+    using Element = typename GridView::template Codim<0>::Entity;
+
+public:
+    using FluidSystem = typename GET_PROP_TYPE(TypeTag, FluidSystem);
+
+    /*!
+     * \brief Update all quantities for a given control volume
+     *
+     * \param elemSol A vector containing all primary variables connected to the element
+     * \param problem The object specifying the problem which ought to be simulated
+     * \param element An element which contains part of the control volume
+     * \param scv The sub-control volume
+     */
+    template<class ElementSolution>
+    void update(const ElementSolution &elemSol,
+                const Problem &problem,
+                const Element &element,
+                const SubControlVolume& scv)
+    {
+        ParentType::update(elemSol, problem, element, scv);
+        updateRANSProperties(elemSol, problem, element, scv);
+    }
+
+    /*!
+     * \brief Update all turbulent quantities for a given control volume
+     *
+     * Wall and roughness related quantities are stored. Eddy viscosity is set.
+     *
+     * \param elemSol A vector containing all primary variables connected to the element
+     * \param problem The object specifying the problem which ought to be simulated
+     * \param element An element which contains part of the control volume
+     * \param scv The sub-control volume
+     */
+    template<class ElementSolution>
+    void updateRANSProperties(const ElementSolution &elemSol,
+                              const Problem &problem,
+                              const Element &element,
+                              const SubControlVolume& scv)
+    {
+        ParentType::updateRANSProperties(elemSol, problem, element, scv);
+        kOmegaModel_ = problem.kOmegaModel();
+        betaOmega_ = problem.betaOmega();
+        turbulentKineticEnergy_ = elemSol[0][Indices::turbulentKineticEnergyIdx];
+        dissipation_ = elemSol[0][Indices::dissipationIdx];
+        storedDissipation_ = problem.storedDissipation_[ParentType::elementID()];
+        storedTurbulentKineticEnergy_ = problem.storedTurbulentKineticEnergy_[ParentType::elementID()];
+        stressTensorScalarProduct_ = problem.stressTensorScalarProduct_[ParentType::elementID()];
+        dofPosition_ = scv.dofPosition();
+        if (problem.useStoredEddyViscosity_)
+            ParentType::setKinematicEddyViscosity(problem.storedKinematicEddyViscosity_[ParentType::elementID()]);
+        else
+            calculateEddyViscosity();
+    }
+
+    /*!
+     * \brief Calculate and set the dynamic eddy viscosity.
+     */
+    void calculateEddyViscosity()
+    {
+        using std::sqrt;
+        using std::max;
+        Scalar kinematicEddyViscosity = 0.0;
+        if(kOmegaModel_ == KOmegaModels::wilcox88)
+            kinematicEddyViscosity = turbulentKineticEnergy() / dissipation();
+        else if (kOmegaModel_ == KOmegaModels::wilcox08)
+          {
+            Scalar limitiedDissipation = (7.0 / 8.0) * std::sqrt( 2.0 * stressTensorScalarProduct() * stressTensorScalarProduct() / betaK() );
+            Scalar Wbar = std::max(dissipation(), limitiedDissipation);
+            kinematicEddyViscosity = turbulentKineticEnergy() / Wbar;
+          }
+        else
+              DUNE_THROW(Dune::NotImplemented, "This model is not implemented.");
+        ParentType::setKinematicEddyViscosity(kinematicEddyViscosity);
+    }
+
+    /*!
+     * \brief Returns the turbulent kinetic energy \f$ m^2/s^2 \f$
+     */
+    Scalar turbulentKineticEnergy() const
+    {
+        return turbulentKineticEnergy_;
+    }
+
+    /*!
+     * \brief Returns an effective dissipation \f$ m^2/s^3 \f$
+     */
+    Scalar dissipation() const
+    {
+        return dissipation_;
+    }
+
+    /*!
+     * \brief Returns the turbulent kinetic energy \f$ m^2/s^2 \f$
+     */
+    Scalar storedTurbulentKineticEnergy() const
+    {
+        return storedTurbulentKineticEnergy_;
+    }
+
+    /*!
+     * \brief Returns an effective dissipation \f$ m^2/s^3 \f$
+     */
+    Scalar storedDissipation() const
+    {
+        return storedDissipation_;
+    }
+
+    /*!
+     * \brief Returns the scalar product of the stress tensor
+     */
+    Scalar stressTensorScalarProduct() const
+    {
+        return stressTensorScalarProduct_;
+    }
+
+    /**
+     *  \brief Here is a list of the background coefficients that are called in the Komega Model
+     */
+
+    //! \brief Returns the \$f \alpha \$f value
+    const Scalar alpha() const
+    {
+        if(kOmegaModel_ == KOmegaModels::wilcox88)
+            return 0.55555555556;
+        else if (kOmegaModel_ == KOmegaModels::wilcox08)
+            return 0.520;
+        else
+            DUNE_THROW(Dune::NotImplemented, "This Model is not implemented.");
+    }
+
+    //! \brief Returns the \$f \sigma_k \$f constant
+    const Scalar sigmaK() const
+    {
+        if(kOmegaModel_ == KOmegaModels::wilcox88)
+            return 0.50;
+        else if (kOmegaModel_ == KOmegaModels::wilcox08)
+            return 0.60;
+        else
+            DUNE_THROW(Dune::NotImplemented, "This Model is not implemented.");
+    }
+
+    //! \brief Returns the \$f \sigma_{\omega} \$f constant
+    const Scalar sigmaOmega() const
+    {
+        if(kOmegaModel_ == KOmegaModels::wilcox88)
+            return 0.50;
+        else if (kOmegaModel_ == KOmegaModels::wilcox08)
+            return 0.50;
+        else
+            DUNE_THROW(Dune::NotImplemented, "This Model is not implemented.");
+    }
+
+    //! \brief Returns the \$f \beta_k \$f constant
+    const Scalar betaK() const
+    {
+        if(kOmegaModel_ == KOmegaModels::wilcox88)
+            return 0.09;
+        else if (kOmegaModel_ == KOmegaModels::wilcox08)
+            return 0.09;
+        else
+            DUNE_THROW(Dune::NotImplemented, "This Model is not implemented.");
+    }
+public:
+    Scalar betaOmega_;
+
+protected:
+    Dune::FieldVector<Scalar,2> dofPosition_;
+    int kOmegaModel_;
+    Scalar turbulentKineticEnergy_;
+    Scalar dissipation_;
+    Scalar storedTurbulentKineticEnergy_;
+    Scalar storedDissipation_;
+    Scalar stressTensorScalarProduct_;
+};
+
+}
+
+#endif
diff --git a/dumux/freeflow/rans/twoeq/komega/vtkoutputfields.hh b/dumux/freeflow/rans/twoeq/komega/vtkoutputfields.hh
new file mode 100644
index 0000000000..5cc68feba1
--- /dev/null
+++ b/dumux/freeflow/rans/twoeq/komega/vtkoutputfields.hh
@@ -0,0 +1,67 @@
+// -*- 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 2 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 KOmegaModel
+ * \copydoc Dumux::KOmegaVtkOutputFields
+ */
+#ifndef DUMUX_KOMEGA_VTK_OUTPUT_FIELDS_HH
+#define DUMUX_KOMEGA_VTK_OUTPUT_FIELDS_HH
+
+#include <dumux/freeflow/rans/vtkoutputfields.hh>
+
+namespace Dumux
+{
+
+/*!
+ * \ingroup KOmegaModel
+ * \brief Adds vtk output fields for the Reynolds-Averaged Navier-Stokes model
+ */
+template<class FVGridGeometry>
+class KOmegaVtkOutputFields : public RANSVtkOutputFields<FVGridGeometry>
+{
+    enum { dim = FVGridGeometry::GridView::dimension };
+
+public:
+    //! Initialize the Navier-Stokes specific vtk output fields.
+    template <class VtkOutputModule>
+    static void init(VtkOutputModule& vtk)
+    {
+        RANSVtkOutputFields<FVGridGeometry>::init(vtk);
+        add(vtk);
+    }
+
+    //! Add the KOmegaModel specific vtk output fields.
+    template <class VtkOutputModule>
+    static void add(VtkOutputModule& vtk)
+    {
+        vtk.addVolumeVariable([](const auto& v){ return v.turbulentKineticEnergy(); }, "turbulentKineticEnergy");
+        vtk.addVolumeVariable([](const auto& v){ return v.dissipation(); }, "dissipation");
+        vtk.addVolumeVariable([](const auto& v){ return v.stressTensorScalarProduct(); }, "stressTensorScalarProduct");
+        vtk.addVolumeVariable([](const auto& v){ return v.kinematicEddyViscosity(); }, "eddyViscosity");
+        vtk.addVolumeVariable([](const auto& v){ return 2.0 * v.kinematicEddyViscosity() * v.stressTensorScalarProduct();}, "production_k");
+        vtk.addVolumeVariable([](const auto& v){ return 2.0 * v.kinematicEddyViscosity() * v.stressTensorScalarProduct() * v.alpha() * ( v.dissipation() / v.turbulentKineticEnergy() ) ; }, "production_omega");
+        vtk.addVolumeVariable([](const auto& v){ return v.betaK() * v.turbulentKineticEnergy() * v.dissipation() ;}, "destruction_k");
+        vtk.addVolumeVariable([](const auto& v){ return v.betaOmega_ * v.dissipation() * v.dissipation() ;}, "destruction_omega");
+    }
+};
+
+} // end namespace Dumux
+
+#endif
-- 
GitLab