[WIP ]Change some things

appl/dispersion1/richardsmimtestproblem.hh

@@ -71,8 +71,8 @@ SET_TYPE_PROP(RichardsMimTestTypeTag, LocalResidual, MIMDispersionLocalResidual<
 
 //the spatial params
 SET_TYPE_PROP(RichardsMimTestTypeTag, SpatialParams, RichardsMimTestSpatialParams<TypeTag>);
-    
-    
+
+
 //the solid system
 SET_PROP(RichardsMimTestTypeTag, SolidSystem)
 {
@@ -132,8 +132,8 @@ class RichardsMimTestProblem : public PorousMediumFlowProblem<TypeTag>
         wPhaseIdx =FluidSystem::liquidPhaseIdx,
         conti0EqIdx = Indices::conti0EqIdx,
         //immobile index
-        immobilePhaseIdx = FluidSystem::numComponents,
-        immobileCompIndex = SolidSystem::comp0Idx,
+        immobilePhaseIdx = FluidSystem::numComponents, // Index of the primary component
+        immobileCompIndex = SolidSystem::comp0Idx,  // index of the solid phase
         dimWorld = GridView::dimensionworld
     };
     using Element = typename GridView::template Codim<0>::Entity;
@@ -286,6 +286,57 @@ public:
 
     // \}
 
+    /*!
+     * \brief Evaluate the source term for all phases within a given
+     *        sub-control-volume.
+     *
+     * This is the method for the case where the source term is
+     * potentially solution dependent and requires some quantities that
+     * are specific to the fully-implicit method.
+     *
+     * \param values The source and sink values for the conservation equations in units of
+     *                 \f$ [ \textnormal{unit of conserved quantity} / (m^3 \cdot s )] \f$
+     * \param element The finite element
+     * \param fvGeometry The finite-volume geometry
+     * \param elemVolVars All volume variables for the element
+     * \param scv The subcontrolvolume
+     *
+     * For this method, the \a values parameter stores the conserved quantity rate
+     * generated or annihilate per volume unit. Positive values mean
+     * that the conserved quantity is created, negative ones mean that it vanishes.
+     * E.g. for the mass balance that would be a mass rate in \f$ [ kg / (m^3 \cdot s)] \f$.
+     */
+    NumEqVector source(const Element &element,
+                   const FVElementGeometry& fvGeometry,
+                   const ElementVolumeVariables& elemVolVars,
+                   const SubControlVolume &scv) const
+    {
+        NumEqVector source(0.0);
+
+        const auto& volVars = elemVolVars[scv];
+
+        Scalar moleFrac_wPhase = volVars.moleFraction(wPhaseIdx, compIdx);
+        //Scalar massFrac_Max_wPhase = this->spatialParams().solubilityLimit();
+        //Scalar moleFrac_Max_wPhase = massToMoleFrac_(massFrac_Max_wPhase);
+
+        // liquid phase
+        using std::abs;
+        Scalar imPhase = volVars.porosity() * volVars.molarDensity(wPhaseIdx)
+                                               * volVars.saturation(wPhaseIdx);
+                        //                       * abs(moleFrac_wPhase - moleFrac_Max_wPhase);
+      //  if (moleFrac_wPhase < moleFrac_Max_wPhase)
+          //  imPhase *= -1;
+
+        // gas phase
+      //  imPhase += volVars.porosity() * volVars.molarDensity(gasPhaseIdx)
+          //                               * volVars.saturation(gasPhaseIdx)
+        //                                 * abs(moleFracNaCl_nPhase - moleFracNaCl_Max_nPhase);
+
+        source[immobilePhaseIdx] += imPhase;
+        return source;
+
+    }
+
 private:
     PrimaryVariables initial_(const GlobalPosition &globalPos) const
     {