wip in updating the derivatives after the two phase solution

it will be a very draft way to do it first, hopefully the results are
okay.

Then looking for better refactorization.

opm/material/constraintsolvers/ChiFlash.hpp

@@ -178,6 +178,12 @@ public:
             std::cout << "********" << std::endl;
         }
 
+        // all the solution should be processed in scalar form
+        // now we should update the derivatives
+        // TODO: should be able the handle the derivatives directly, which will affect the final organization
+        // of the code
+        updateDerivatives_(fluid_state_scalar, fluid_state);
+
         // Update phases
         /* typename FluidSystem::template ParameterCache<InputEval> paramCache;
         paramCache.updatePhase(fluid_state, oilPhaseIdx);
@@ -956,6 +962,26 @@ protected:
         }
     }
 
+    template <typename FlashFluidStateScalar, typename FlashFluidState>
+    static void updateDerivatives_(const FlashFluidStateScalar& fluid_state_scalar,
+                                         FlashFluidState& fluid_state)
+    {
+        // getting the secondary Jocobian matrix
+        constexpr size_t num_equations = numMisciblePhases * numMiscibleComponents + 1;
+        constexpr size_t secondary_num_pv = numComponents + 1;
+        using SecondaryEval = Opm::DenseAd::Evaluation<double, secondary_num_pv>; // three z and one pressure
+        using SecondaryComponentVector = Dune::FieldVector<SecondaryEval, numComponents>;
+        using SecondaryFlashFluidState = Opm::CompositionalFluidState<SecondaryEval, FluidSystem>;
+
+        SecondaryFlashFluidState secondary_fluid_state;
+        SecondaryComponentVector secondary_z;
+
+        // TODO: trying to get the secondary matrix here,
+        // TODO: then getting the primary matrix,
+        // TODO: then beginning from that point
+
+    }
+
     /* template <class Vector, class Matrix, class Eval, class ComponentVector>
     static void evalJacobian(const ComponentVector& globalComposition,
                              const Vector& x,

tests/test_chiflash.cpp

@@ -104,6 +104,14 @@ void testChiFlash()
             }
         }
         zInit /= sumMoles;
+        // initialize the derivatives
+        // TODO: the derivative eventually should be from the reservoir flow equations
+        Evaluation z_last = 1.;
+        for (unsigned compIdx = 0; compIdx < numComponents - 1; ++compIdx) {
+            zInit[compIdx] = Evaluation::createVariable(Opm::getValue(zInit[compIdx]), compIdx + 1);
+            z_last -= zInit[compIdx];
+        }
+        zInit[numComponents - 1] = z_last;
     }
 
     // TODO: only, p, z need the derivatives.