Merge branch 'feature/napl-3p-test' into 'master'

[naplinfiltration] Separate 3p and 3p3c for pedagogical reasons

* Make the ctests run only about 30 seconds each to save computer power
  when testing automated

* [ ] This can be only merged if the corresponding merge request of dumux
branch fix/3p-parker-vangenuchten is merged.

See merge request !10

lecture/mm/columnxylene/columnxylenespatialparams.hh

@@ -152,6 +152,10 @@ public:
         coarseMaterialParams_.setRhoBulk(1500.);
         fineMaterialParams_.setKdNAPL(0.);
         fineMaterialParams_.setRhoBulk(1500.);
+
+        // use capping as regularization
+        coarseMaterialParams_.useConstRegularization(true);
+        fineMaterialParams_.useConstRegularization(true);
     }
 
     /*!

lecture/mm/naplinfiltration/3p/CMakeLists.txt

+add_input_file_links()
+
+# for the test only simulate one week
+add_dumux_test(naplinfiltration3p naplinfiltration3p naplinfiltration3p.cc
+               python ${dumux_INCLUDE_DIRS}/bin/runtest.py
+                 --script fuzzy
+                 --files ${CMAKE_SOURCE_DIR}/lecture/references/infiltration3p-reference.vtu
+                         ${CMAKE_CURRENT_BINARY_DIR}/infiltration3p-00004.vtu
+                 --command "${CMAKE_CURRENT_BINARY_DIR}/naplinfiltration3p -TimeManager.TEnd 604800 -TimeManager.EpisodeLength 302400")
+
+# headers for installation and headercheck
+install(FILES
+        problem.hh
+        naplinfiltration3p.cc
+        naplinfiltration3p.input
+        DESTINATION ${CMAKE_INSTALL_INCLUDEDIR}/lecture/mm/naplinfiltration/3p)

lecture/mm/naplinfiltration/naplinfiltrationexercise.cc → lecture/mm/naplinfiltration/3p/naplinfiltration3p.cc

@@ -19,15 +19,10 @@
 /*!
  * \file
  *
- * \brief application, using the 3p3c or 3p box model
+ * \brief application, using the 3p box model
  */
 #include "config.h"
-
-// Boolean that decides whether to use the model with three miscible components (3p3c, default)
-// or the model with three immiscible components (3p, change this to 1)
-#define USE_THREE_P_MODEL 0
-
-#include "naplinfiltrationproblem.hh"
+#include "problem.hh"
 #include <dumux/common/start.hh>
 
 /*!
@@ -58,11 +53,6 @@ void usage(const char *progName, const std::string &errorMsg)
 
 int main(int argc, char** argv)
 {
-#if USE_THREE_P_MODEL
-    typedef TTAG(InfiltrationProblemThreeP) ProblemTypeTag;
-#else
-    typedef TTAG(InfiltrationProblemThreePThreeC) ProblemTypeTag;
-#endif
+    typedef TTAG(InfiltrationProblem) ProblemTypeTag;
     return Dumux::start<ProblemTypeTag>(argc, argv, usage);
 }
-

lecture/mm/naplinfiltration/naplinfiltrationexercise.input → lecture/mm/naplinfiltration/3p/naplinfiltration3p.input

@@ -2,7 +2,7 @@
 DtInitial = 60          # [s]
 TEnd = 31536000         # [s]
 #TEnd = 315360000000     # [s]
-EpisodeLength = 3153600 # [s]
+EpisodeLength = 86400 # [s]
 
 [Grid]
 UpperRight = 500 10  # coordinates of upper right grid corner [m]

lecture/mm/naplinfiltration/3p/problem.hh

+// -*- 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
+ *
+ * \brief Isothermal NAPL infiltration problem: LNAPL contaminates
+ *        the unsaturated and the saturated groundwater zone.
+ */
+#ifndef DUMUX_NAPLINFILTRATIONPROBLEM_3P_HH
+#define DUMUX_NAPLINFILTRATIONPROBLEM_3P_HH
+
+#include <dumux/porousmediumflow/3p/implicit/propertydefaults.hh>
+#include <dumux/porousmediumflow/implicit/problem.hh>
+
+#include <lecture/mm/naplinfiltration/h2oairnaplfluidsystem.hh>
+#include <lecture/mm/naplinfiltration/spatialparams.hh>
+
+namespace Dumux
+{
+template <class TypeTag>
+class InfiltrationProblem;
+
+namespace Properties
+{
+NEW_TYPE_TAG(InfiltrationProblem, INHERITS_FROM(BoxThreeP, InfiltrationSpatialParams));
+
+// Set the grid type
+SET_TYPE_PROP(InfiltrationProblem, Grid, Dune::YaspGrid<2>);
+
+// Set the problem property
+SET_TYPE_PROP(InfiltrationProblem, Problem, Dumux::InfiltrationProblem<TypeTag>);
+
+// Set the fluid system
+SET_TYPE_PROP(InfiltrationProblem,
+              FluidSystem,
+              Dumux::FluidSystems::H2OAirNAPL<TypeTag, typename GET_PROP_TYPE(TypeTag, Scalar)>);
+
+// Enable gravity?
+SET_BOOL_PROP(InfiltrationProblem, ProblemEnableGravity, true);
+
+// Write newton convergence?
+SET_BOOL_PROP(InfiltrationProblem, NewtonWriteConvergence, false);
+
+// Maximum tolerated relative error in the Newton method
+SET_SCALAR_PROP(InfiltrationProblem, NewtonMaxRelativeShift, 1e-4);
+
+// -1 backward differences, 0: central differences, +1: forward differences
+SET_INT_PROP(InfiltrationProblem, ImplicitNumericDifferenceMethod, 1);
+}
+
+/*!
+ * \ingroup ThreePThreeCBoxModel
+ * \ingroup BoxTestProblems
+ * \brief Isothermal NAPL infiltration problem: LNAPL contaminates
+ *        the unsaturated and the saturated groundwater zone.
+ *
+ * The 2D domain of this test problem is 500 m long and 10 m deep, where
+ * the lower part represents a slightly inclined groundwater table, and the
+ * upper part is the vadose zone.
+ * A LNAPL (Non-Aqueous Phase Liquid which is lighter than water) infiltrates
+ * (modelled with a Neumann boundary condition) into the vadose zone. Upon
+ * reaching the water table, it spreads (since lighter than water) and migrates
+ * on top of the water table in the direction of the slope.
+ * On its way through the vadose zone, it leaves a trace of residually trapped
+ * immobile NAPL, which can in the following dissolve and evaporate slowly,
+ * and eventually be transported by advection and diffusion.
+ *
+ * Left and right boundaries are constant head boundaries (Dirichlet),
+ * Top and bottom are Neumann boundaries, all no-flow except for the small
+ * infiltration zone in the upper left part.
+ *
+ * This problem uses the \ref ThreePThreeCModel or the \ref ThreePModel.
+ *
+ * This problem should typically be simulated for 30 days.
+ * A good choice for the initial time step size is 60 s.
+ * To adjust the simulation time it is necessary to edit the file test_naplfiltrationexercise.input
+ *
+ * To run the simulation execute the following line in shell:
+ * <tt>./test_naplfiltrationexercise -ParameterFile test_naplfiltrationexercise.input</tt>
+ *  */
+template <class TypeTag >
+class InfiltrationProblem : public ImplicitPorousMediaProblem<TypeTag>
+{
+    typedef ImplicitPorousMediaProblem<TypeTag> ParentType;
+
+    typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
+    typedef typename GET_PROP_TYPE(TypeTag, GridView) GridView;
+
+    // copy some indices for convenience
+    typedef typename GET_PROP_TYPE(TypeTag, Indices) Indices;
+    enum {
+        pressureIdx = Indices::pressureIdx,
+        switch1Idx = Indices::swIdx,
+        switch2Idx = Indices::snIdx,
+
+        contiWEqIdx = Indices::contiWEqIdx,
+        contiNEqIdx = Indices::contiNEqIdx,
+        contiGEqIdx = Indices::contiGEqIdx,
+
+        // Grid and world dimension
+        dim = GridView::dimension,
+        dimWorld = GridView::dimensionworld
+    };
+
+
+    typedef typename GET_PROP_TYPE(TypeTag, PrimaryVariables) PrimaryVariables;
+    typedef typename GET_PROP_TYPE(TypeTag, BoundaryTypes) BoundaryTypes;
+    typedef typename GET_PROP_TYPE(TypeTag, TimeManager) TimeManager;
+
+    typedef typename GridView::template Codim<0>::Entity Element;
+    typedef typename GridView::template Codim<dim>::Entity Vertex;
+
+    typedef typename GET_PROP_TYPE(TypeTag, FVElementGeometry) FVElementGeometry;
+    typedef typename GET_PROP_TYPE(TypeTag, FluidSystem) FluidSystem;
+    typedef typename GET_PROP_TYPE(TypeTag, MaterialLaw) MaterialLaw;
+    typedef typename GET_PROP_TYPE(TypeTag, MaterialLawParams) MaterialLawParams;
+
+    typedef Dune::FieldVector<Scalar, dimWorld> GlobalPosition;
+
+public:
+    /*!
+     * \brief The constructor
+     *
+     * \param timeManager The time manager
+     * \param gridView The grid view
+     */
+    InfiltrationProblem(TimeManager &timeManager, const GridView &gridView)
+        : ParentType(timeManager, gridView), eps_(1e-5)
+    {
+        temperature_ = 273.15 + 10.0; // -> 10 degrees Celsius
+
+        episodeLength_ = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, Scalar, TimeManager, EpisodeLength);
+        this->timeManager().startNextEpisode(episodeLength_);
+
+        FluidSystem::init(/*tempMin=*/temperature_ - 1,
+                          /*tempMax=*/temperature_ + 1,
+                          /*nTemp=*/3,
+                          /*pressMin=*/0.8*1e5,
+                          /*pressMax=*/3*1e5,
+                          /*nPress=*/200);
+    }
+
+    /*!
+     * \name Problem parameters
+     */
+    // \{
+
+    /*!
+     * \brief The problem name.
+     *
+     * This is used as a prefix for files generated by the simulation.
+     */
+    std::string name() const
+    { return "infiltration3p"; }
+
+    /*!
+     * \brief Returns the temperature within the domain.
+     *
+     * This problem assumes a temperature of 10 degrees Celsius.
+     */
+    Scalar temperature() const
+    {
+        return temperature_;
+    };
+
+    void sourceAtPos(PrimaryVariables &values,
+                     const GlobalPosition &globalPos) const
+    {
+        values = 0;
+    }
+
+    // \}
+
+    /*!
+     * \name Boundary conditions
+     */
+    // \{
+
+    /*!
+     * \brief Specifies which kind of boundary condition should be
+     *        used for which equation on a given boundary segment.
+     *
+     * \param values The boundary types for the conservation equations
+     * \param globalPos The position of the finite volume in global coordinates
+     */
+    void boundaryTypesAtPos(BoundaryTypes &values,
+                            const GlobalPosition &globalPos) const
+    {
+        if(globalPos[0] > this->bBoxMax()[0] - eps_)
+            values.setAllDirichlet();
+        else if(globalPos[0] < this->bBoxMin()[0] + eps_)
+            values.setAllDirichlet();
+        else
+            values.setAllNeumann();
+    }
+
+    /*!
+     * \brief Evaluate the boundary conditions for a dirichlet
+     *        control volume.
+     *
+     * \param values The dirichlet values for the primary variables
+     * \param globalPos The position of the center of the finite volume
+     *            for which the dirichlet condition ought to be
+     *            set in global coordinates
+     *
+     * For this method, the \a values parameter stores primary variables.
+     */
+    void dirichletAtPos(PrimaryVariables &values,
+                        const GlobalPosition &globalPos) const
+    {
+        initial_(values, globalPos);
+    }
+
+    /*!
+     * \brief Evaluate the boundary conditions for a neumann
+     *        boundary segment.
+     *
+     * \param values The neumann values for the conservation equations in units of
+     *                 \f$ [ \textnormal{unit of conserved quantity} / (m^2 \cdot s )] \f$
+     * \param globalPos The position of the boundary face's integration point in global coordinates
+     *
+     * For this method, the \a values parameter stores the mass flux
+     * in normal direction of each phase. Negative values mean influx.
+     */
+    void neumannAtPos(PrimaryVariables &values,
+                      const GlobalPosition &globalPos) const
+    {
+        values = 0;
+
+        // negative values for injection
+        if (this->timeManager().time() < 2592000)
+        {
+            if (globalPos[0] < 175.0 + eps_
+                && globalPos[0] > 150.0 - eps_
+                && globalPos[1] > this->bBoxMax()[1] - eps_) // upper boundary
+            {
+                values[contiWEqIdx] = 0.0;
+                // mole flux conversion via M(mesit.) = 0,120 kg/mol --> 1.2e-4  kg/(s*m)
+                values[contiNEqIdx] = -0.001*0.12; // 3p needs a mass fraction
+                values[contiGEqIdx] = 0.0;
+            }
+        }
+    }
+
+    // \}
+
+    /*!
+     * \name Volume terms
+     */
+    // \{
+
+    /*!
+     * \brief Evaluate the initial value for a control volume.
+     *
+     * \param values The dirichlet values for the primary variables
+     * \param globalPos The position of the center of the finite volume
+     *            for which the initial values ought to be
+     *            set (in global coordinates)
+     *
+     * For this method, the \a values parameter stores primary variables.
+     */
+    void initialAtPos(PrimaryVariables &values,
+                      const GlobalPosition &globalPos) const
+    {
+        initial_(values, globalPos);
+    }
+
+    bool shouldWriteOutput() const
+    {
+    return this->timeManager().timeStepIndex() == 0 ||
+           this->timeManager().episodeWillBeOver() ||
+           this->timeManager().willBeFinished();
+    }
+
+    void episodeEnd()
+    {
+        this->timeManager().startNextEpisode(episodeLength_);
+    }
+
+private:
+    // internal method for the initial condition (reused for the
+    // dirichlet conditions!)
+    void initial_(PrimaryVariables &values,
+                  const GlobalPosition &globalPos) const
+    {
+        const auto& materialLawParams = this->spatialParams().materialLawParams();
+        const Scalar swr = materialLawParams.swr();
+        const Scalar sgr = materialLawParams.sgr();
+
+        if(globalPos[1] >= -1e-3*globalPos[0] + 5)
+        {
+            const Scalar pc = std::max(0.0, 9.81*1000.0*(globalPos[1] + 5e-4*globalPos[0] - 5));
+            const Scalar sw = std::min(1.0-sgr, std::max(swr, invertPcGW_(pc, materialLawParams)));
+
+            values[pressureIdx] = 1.0e5;
+            values[switch1Idx] = sw;
+            values[switch2Idx] = 0.0;
+        }
+        else
+        {
+            values[pressureIdx] = 1.0e5 + 9.81*1000.0*(5 - 5.0e-4*globalPos[0] - globalPos[1]);
+            values[switch1Idx] = 1.0 - sgr;
+            values[switch2Idx] = 0.0;
+        }
+    }
+
+    static Scalar invertPcGW_(const Scalar pcIn,
+                              const MaterialLawParams &pcParams)
+    {
+        Scalar lower(0.0);
+        Scalar upper(1.0);
+        const unsigned int maxIterations = 25;
+        const Scalar bisLimit = 1.0;
+
+        Scalar sw, pcGW;
+        for (unsigned int k = 1; k <= maxIterations; k++)
+        {
+            sw = 0.5*(upper + lower);
+            pcGW = MaterialLaw::pcgw(pcParams, sw);
+            const Scalar delta = std::abs(pcGW - pcIn);
+            if (delta < bisLimit)
+                return sw;
+
+            if (k == maxIterations)
+                return sw;
+
+            if (pcGW > pcIn)
+                lower = sw;
+            else
+                upper = sw;
+        }
+        return sw;
+    }
+
+    Scalar temperature_;
+    const Scalar eps_;
+    Scalar episodeLength_;
+};
+} //end namespace
+
+#endif

lecture/mm/naplinfiltration/3p3c/CMakeLists.txt

+add_input_file_links()
+
+# for the test only simulate one week
+add_dumux_test(naplinfiltration3p3c naplinfiltration3p3c naplinfiltration3p3c.cc
+               python ${dumux_INCLUDE_DIRS}/bin/runtest.py
+                 --script fuzzy
+                 --files ${CMAKE_SOURCE_DIR}/lecture/references/infiltration3p3c-reference.vtu
+                         ${CMAKE_CURRENT_BINARY_DIR}/infiltration3p3c-00004.vtu
+                 --command "${CMAKE_CURRENT_BINARY_DIR}/naplinfiltration3p3c -TimeManager.TEnd 604800 -TimeManager.EpisodeLength 302400")
+
+# headers for installation and headercheck
+install(FILES
+        problem.hh
+        naplinfiltration3p3c.cc
+        naplinfiltration3p3c.input
+        DESTINATION ${CMAKE_INSTALL_INCLUDEDIR}/lecture/mm/naplinfiltration/3p3c)

lecture/mm/naplinfiltration/3p3c/naplinfiltration3p3c.cc

+// -*- 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
+ *
+ * \brief application, using the 3p3c box model
+ */
+#include "config.h"
+#include "problem.hh"
+#include <dumux/common/start.hh>
+
+/*!
+ * \brief Provides an interface for customizing error messages associated with
+ *        reading in parameters.
+ *
+ * \param progName  The name of the program, that was tried to be started.
+ * \param errorMsg  The error message that was issued by the start function.
+ *                  Comprises the thing that went wrong and a general help message.
+ */
+void usage(const char *progName, const std::string &errorMsg)
+{
+    if (errorMsg.size() > 0) {
+        std::string errorMessageOut = "\nUsage: ";
+                    errorMessageOut += progName;
+                    errorMessageOut += " [options]\n";
+                    errorMessageOut += errorMsg;
+                    errorMessageOut += "\n\nThe list of mandatory options for this program is:\n"
+                                        "\t-TimeManager.TEnd              End of the simulation [s] \n"
+                                        "\t-TimeManager.DtInitial         Initial timestep size [s] \n"
+                                        "\t-Grid.UpperRight               Coordinates of upper right grid corner [m] \n"
+                                        "\t-Grid.Cells                    Number of cells in x, y direction [-] \n";
+
+        std::cout << errorMessageOut
+                  << "\n";
+    }
+}
+
+int main(int argc, char** argv)
+{
+    typedef TTAG(InfiltrationProblem) ProblemTypeTag;
+    return Dumux::start<ProblemTypeTag>(argc, argv, usage);
+}

lecture/mm/naplinfiltration/3p3c/naplinfiltration3p3c.input

+[TimeManager]
+DtInitial = 60          # [s]
+TEnd = 31536000         # [s]
+#TEnd = 315360000000     # [s]
+EpisodeLength = 86400 # [s]
+
+[Grid]
+UpperRight = 500 10  # coordinates of upper right grid corner [m]
+Cells = 250 10       # number of cells in (x, y) direction [-]
+
+[SpatialParams]
+Permeability = 1.e-11        # [m^2]
+Porosity = 0.40              # [-]
+VanGenuchtenAlpha = 0.0005   # [-]
+VanGenuchtenN = 4.0          # [-]
+MultiplierNaplDensity = 1.0  # 1.0 means Mesitylene (an LNAPL) [-]
+
+[Newton]
+MaxSteps = 6
\ No newline at end of file

lecture/mm/naplinfiltration/naplinfiltrationproblem.hh → lecture/mm/naplinfiltration/3p3c/problem.hh

@@ -22,15 +22,14 @@
  * \brief Isothermal NAPL infiltration problem: LNAPL contaminates
  *        the unsaturated and the saturated groundwater zone.
  */
-#ifndef DUMUX_NAPLINFILTRATIONPROBLEM_HH
-#define DUMUX_NAPLINFILTRATIONPROBLEM_HH
+#ifndef DUMUX_NAPLINFILTRATIONPROBLEM_3P_3C_HH
+#define DUMUX_NAPLINFILTRATIONPROBLEM_3P_3C_HH
 
-#include <dumux/porousmediumflow/3p/implicit/propertydefaults.hh>
 #include <dumux/porousmediumflow/3p3c/implicit/propertydefaults.hh>
 #include <dumux/porousmediumflow/implicit/problem.hh>
 
-#include "h2oairnaplfluidsystem.hh"
-#include "naplinfiltrationspatialparams.hh"
+#include <lecture/mm/naplinfiltration/h2oairnaplfluidsystem.hh>
+#include <lecture/mm/naplinfiltration/spatialparams.hh>
 
 namespace Dumux
 {
@@ -39,9 +38,7 @@ class InfiltrationProblem;
 
 namespace Properties
 {
-NEW_TYPE_TAG(InfiltrationProblem, INHERITS_FROM(InfiltrationSpatialParams));
-NEW_TYPE_TAG(InfiltrationProblemThreeP, INHERITS_FROM(BoxThreeP, InfiltrationProblem));
-NEW_TYPE_TAG(InfiltrationProblemThreePThreeC, INHERITS_FROM(BoxThreePThreeC, InfiltrationProblem));
+NEW_TYPE_TAG(InfiltrationProblem, INHERITS_FROM(BoxThreePThreeC, InfiltrationSpatialParams));
 
 // Set the grid type
 SET_TYPE_PROP(InfiltrationProblem, Grid, Dune::YaspGrid<2>);
@@ -109,14 +106,10 @@ class InfiltrationProblem : public ImplicitPorousMediaProblem<TypeTag>
     typedef typename GET_PROP_TYPE(TypeTag, Indices) Indices;
     enum {
         pressureIdx = Indices::pressureIdx,
-#if USE_THREE_P_MODEL
-        switch1Idx = Indices::swIdx,
-        switch2Idx = Indices::snIdx,
-#else
+
         switch1Idx = Indices::switch1Idx,
         switch2Idx = Indices::switch2Idx,
         wgPhaseOnly = Indices::wgPhaseOnly,
-#endif
 
         contiWEqIdx = Indices::contiWEqIdx,
         contiNEqIdx = Indices::contiNEqIdx,
@@ -175,8 +168,8 @@ public:
      *
      * This is used as a prefix for files generated by the simulation.
      */
-    const char *name() const
-    { return "infiltration"; }
+    std::string name() const
+    { return "infiltration3p3c"; }
 
     /*!
      * \brief Returns the temperature within the domain.
@@ -260,8 +253,8 @@ public:
                 && globalPos[1] > this->bBoxMax()[1] - eps_) // upper boundary
             {
                 values[contiWEqIdx] = 0.0;
-                // /* mole flux, conversion via M(mesit.) = 0,120 kg/mol --> 1.2e-4  kg/(s*m)
-                values[contiNEqIdx] = -0.001;
+                // mole flux conversion via M(mesit.) = 0,120 kg/mol --> 1.2e-4  kg/(s*m)
+                values[contiNEqIdx] = -0.001; // 3p3c needs a mole fraction
                 values[contiGEqIdx] = 0.0;
             }
         }
@@ -290,7 +283,6 @@ public:
         initial_(values, globalPos);
     }
 
-#if !USE_THREE_P_MODEL
     /*!
      * \brief Return the initial phase state inside a control volume.
      *
@@ -304,7 +296,6 @@ public:
     {
         return Indices::wgPhaseOnly;
     }
-#endif
 
     bool shouldWriteOutput() const
     {