Adjoint new try

ewoms/common/parametersystem.hh

@@ -961,7 +961,6 @@ public:
         check_(Dune::className<ParamType>(), propTagName, paramName);
 #endif
 
-        typedef typename GET_PROP(TypeTag, ParameterMetaData) ParamsMeta;
         if (errorIfNotRegistered) {
             if (ParamsMeta::registrationOpen())
                 throw std::runtime_error("Parameters can only checked after _all_ of them have "
@@ -1038,7 +1037,6 @@ private:
         check_(Dune::className<ParamType>(), propTagName, paramName);
 #endif
 
-        typedef typename GET_PROP(TypeTag, ParameterMetaData) ParamsMeta;
         if (errorIfNotRegistered) {
             if (ParamsMeta::registrationOpen())
                 throw std::runtime_error("Parameters can only retieved after _all_ of them have "

ewoms/common/simulator.hh

@@ -56,6 +56,7 @@ NEW_PROP_TAG(EndTime);
 NEW_PROP_TAG(RestartTime);
 NEW_PROP_TAG(InitialTimeStepSize);
 NEW_PROP_TAG(PredeterminedTimeStepsFile);
+NEW_PROP_TAG(SimulatorManageTimeStep);
 
 END_PROPERTIES
 
@@ -107,6 +108,8 @@ class Simulator
     typedef typename GET_PROP_TYPE(TypeTag, Model) Model;
     typedef typename GET_PROP_TYPE(TypeTag, Problem) Problem;
 
+    enum { simulatorManageTimeStep = GET_PROP_VALUE(TypeTag, SimulatorManageTimeStep) };
+
 public:
     // do not allow to copy simulators around
     Simulator(const Simulator& ) = delete;
@@ -866,11 +869,11 @@ public:
      * name and uses the extension <tt>.ers</tt>. (Ewoms ReStart
      * file.)  See Ewoms::Restart for details.
      */
-    void serialize()
+    void serialize(const bool atEndOfStep = true)
     {
         typedef Ewoms::Restart Restarter;
         Restarter res;
-        res.serializeBegin(*this);
+        res.serializeBegin(*this, atEndOfStep);
         if (gridView().comm().rank() == 0)
             std::cout << "Serialize to file '" << res.fileName() << "'"
                       << ", next time step size: " << timeStepSize()
@@ -882,6 +885,23 @@ public:
         res.serializeEnd();
     }
 
+    void deserializeAll(Scalar t,bool only_reservoir)
+    {
+        typedef Ewoms::Restart Restarter;
+        Restarter res;
+        res.deserializeBegin(*this, t);
+        if (gridView().comm().rank() == 0)
+            std::cout << "Deserialize file '" << res.fileName() << "'"
+                      << ", next time step size: " << timeStepSize()
+                      << "\n" << std::flush;
+        this->deserialize(res);
+        if( not(only_reservoir) ){
+            problem_->deserialize(res, false);
+        }
+         model_->deserialize(res);
+
+        res.deserializeEnd();
+    }
     /*!
      * \brief Write the time manager's state to a restart file.
      *
@@ -899,6 +919,7 @@ public:
             << episodeLength_ << " "
             << startTime_ << " "
             << time_ << " "
+            << timeStepSize_ << " "
             << timeStepIdx_ << " ";
         restarter.serializeSectionEnd();
     }
@@ -920,6 +941,7 @@ public:
             >> episodeLength_
             >> startTime_
             >> time_
+            >> timeStepSize_
             >> timeStepIdx_;
         restarter.deserializeSectionEnd();
     }

ewoms/disc/common/fvbasediscretization.hh

@@ -268,6 +268,7 @@ SET_BOOL_PROP(FvBaseDiscretization, UseVolumetricResidual, true);
 
 //! eWoms is mainly targeted at research, so experimental features are enabled by
 //! default.
+SET_BOOL_PROP(FvBaseDiscretization, SimulatorManageTimeStep, true);
 SET_BOOL_PROP(FvBaseDiscretization, EnableExperiments, true);
 
 END_PROPERTIES
@@ -759,6 +760,15 @@ public:
     bool enableStorageCache() const
     { return enableStorageCache_; }
 
+    /*!
+     * \brief Set the value of enable storage cache  
+     *
+     * Be aware that calling the *CachedStorage() methods if the storage cache is
+     * disabled will crash the program.
+     */
+    void setEnableStorageCache(bool enableStorageCache) 
+    { enableStorageCache_= enableStorageCache; }
+
     /*!
      * \brief Retrieve an entry of the cache for the storage term.
      *

ewoms/disc/common/fvbaseelementcontext.hh

@@ -97,9 +97,10 @@ public:
     {
         // remember the simulator object
         simulatorPtr_ = &simulator;
-        enableStorageCache_ = EWOMS_GET_PARAM(TypeTag, bool, EnableStorageCache);
+        enableStorageCache_ = simulator.model().enableStorageCache();//EWOMS_GET_PARAM(TypeTag, bool, EnableStorageCache);
         stashedDofIdx_ = -1;
         focusDofIdx_ = -1;
+        focusTimeIdx_ = 0;
     }
 
     static void *operator new(size_t size)
@@ -260,6 +261,17 @@ public:
     void setFocusDofIndex(unsigned dofIdx)
     { focusDofIdx_ = dofIdx; }
 
+    /*!
+     * \brief Sets the time index on which the simulator is currently "focused" on
+     *
+     * I.e., in the case of automatic differentiation, all derivatives are with regard to
+     * the primary variables of that time index. Only "primary" DOFs can be
+     * focused on.
+     */
+
+    void setFocusTimeIndex(unsigned timeIdx)
+    { focusTimeIdx_ = timeIdx; }
+
     /*!
      * \brief Returns the degree of freedom on which the simulator is currently "focused" on
      *
@@ -268,6 +280,14 @@ public:
     unsigned focusDofIndex() const
     { return focusDofIdx_; }
 
+    /*!
+     * \brief Returns the time index on which the simulator is currently "focused" on
+     *
+     * \copydetails setFocusDof()
+     */ 
+    unsigned focusTimeIndex() const
+    { return focusTimeIdx_; }
+
     /*!
      * \brief Return a reference to the simulator.
      */
@@ -581,7 +601,7 @@ protected:
 #endif
 
         dofVars_[dofIdx].priVars[timeIdx] = priVars;
-        dofVars_[dofIdx].intensiveQuantities[timeIdx].update(/*context=*/asImp_(), dofIdx, timeIdx);
+        dofVars_[dofIdx].intensiveQuantities[timeIdx].update(/*context=*/asImp_(), dofIdx, timeIdx, focusTimeIdx_);
     }
 
     IntensiveQuantities intensiveQuantitiesStashed_;
@@ -599,6 +619,7 @@ protected:
 
     int stashedDofIdx_;
     int focusDofIdx_;
+    int focusTimeIdx_;
     bool enableStorageCache_;
 };
 

ewoms/disc/common/fvbaseintensivequantities.hh

@@ -67,7 +67,8 @@ public:
      */
     void update(const ElementContext& elemCtx,
                 unsigned dofIdx,
-                unsigned timeIdx)
+                unsigned timeIdx,
+                unsigned focusTimeIdx OPM_UNUSED)
     { extrusionFactor_ = elemCtx.problem().extrusionFactor(elemCtx, dofIdx, timeIdx); }
 
     /*!

ewoms/disc/common/fvbaselinearizer.hh

@@ -143,6 +143,12 @@ public:
     {
         simulatorPtr_ = &simulator;
         eraseMatrix();
+        auto it = elementCtx_.begin();
+        const auto& endIt = elementCtx_.end();
+        for (; it != endIt; ++it){
+            delete *it;
+        }
+        elementCtx_.resize(0);
     }
 
     /*!
@@ -158,14 +164,15 @@ public:
     }
 
     /*!
-     * \brief Linearize the full system of non-linear equations.
+     * \brief Linearize the full system of non-linear equations about focusTimeIdx
+     * defualt is 0 i.e. current time normally time to be found.
      *
      * This means the spatial domain plus all auxiliary equations.
      */
-    void linearize()
+    void linearize(unsigned focusTimeIdx = 0)
     {
-        linearizeDomain();
-        linearizeAuxiliaryEquations();
+        linearizeDomain(focusTimeIdx);
+        linearizeAuxiliaryEquations(focusTimeIdx);
     }
 
     /*!
@@ -178,7 +185,7 @@ public:
      * The current state of affairs (esp. the previous and the current solutions) is
      * represented by the model object.
      */
-    void linearizeDomain()
+    void linearizeDomain(unsigned focustimeIndex)
     {
         // we defer the initialization of the Jacobian matrix until here because the
         // auxiliary modules usually assume the problem, model and grid to be fully
@@ -188,7 +195,7 @@ public:
 
         int succeeded;
         try {
-            linearize_();
+            linearize_(focustimeIndex);
             succeeded = 1;
         }
 #if ! DUNE_VERSION_NEWER(DUNE_COMMON, 2,5)
@@ -227,7 +234,7 @@ public:
      * \brief Linearize the part of the non-linear system of equations that is associated
      *        with the spatial domain.
      */
-    void linearizeAuxiliaryEquations()
+    void linearizeAuxiliaryEquations(unsigned focusTimeIdx OPM_UNUSED = 0)
     {
         // flush possible local caches into matrix structure
         jacobian_->commit();
@@ -237,7 +244,7 @@ public:
         for (unsigned auxModIdx = 0; auxModIdx < model.numAuxiliaryModules(); ++auxModIdx) {
             bool succeeded = true;
             try {
-                model.auxiliaryModule(auxModIdx)->linearize(*jacobian_, residual_);
+                model.auxiliaryModule(auxModIdx)->linearize(*jacobian_, residual_);//, focusTimeIdx);
             }
             catch (const std::exception& e) {
                 succeeded = false;
@@ -330,6 +337,7 @@ private:
             elementCtx_[threadId] = new ElementContext(simulator_());
     }
 
+
     // Construct the BCRS matrix for the Jacobian of the residual function
     void createMatrix_()
     {
@@ -425,7 +433,7 @@ private:
     }
 
     // linearize the whole system
-    void linearize_()
+    void linearize_(unsigned focustimeindex)
     {
         resetSystem_();
 
@@ -472,8 +480,7 @@ private:
                     const Element& elem = *elemIt;
                     if (!linearizeNonLocalElements && elem.partitionType() != Dune::InteriorEntity)
                         continue;
-
-                    linearizeElement_(elem);
+                    linearizeElement_(elem, focustimeindex);
                 }
             }
             // If an exception occurs in the parallel block, it won't escape the
@@ -503,14 +510,15 @@ private:
     }
 
     // linearize an element in the interior of the process' grid partition
-    void linearizeElement_(const Element& elem)
+    void linearizeElement_(const Element& elem, unsigned focustimeindex)
     {
         unsigned threadId = ThreadManager::threadId();
 
         ElementContext *elementCtx = elementCtx_[threadId];
         auto& localLinearizer = model_().localLinearizer(threadId);
 
         // the actual work of linearization is done by the local linearizer class
+        elementCtx->setFocusTimeIndex(focustimeindex);
         localLinearizer.linearize(*elementCtx, elem);
 
         // update the right hand side and the Jacobian matrix

ewoms/disc/common/fvbaselocalresidual.hh

@@ -163,16 +163,17 @@ public:
         assert(residual.size() == elemCtx.numDof(/*timeIdx=*/0));
 
         residual = 0.0;
-
+        if(elemCtx.focusTimeIndex()==0){//NB: This is only valid if fluxeterm sis purely implicite, This was nesseary to not get false derivatives form wells in adjiont.
         // evaluate the flux terms
-        asImp_().evalFluxes(residual, elemCtx, /*timeIdx=*/0);
-
+            asImp_().evalFluxes(residual, elemCtx, /*timeIdx=*/0);
+        }
         // evaluate the storage and the source terms
-        asImp_().evalVolumeTerms_(residual, elemCtx);
 
+        asImp_().evalVolumeTerms_(residual, elemCtx);
+        if(elemCtx.focusTimeIndex()==0){//NB: This is only valid if fluxeterm sis purely implicite, This was nesseary to not get false derivatives form wells in adjiont.
         // evaluate the boundary conditions
-        asImp_().evalBoundary_(residual, elemCtx, /*timeIdx=*/0);
-
+            asImp_().evalBoundary_(residual, elemCtx, /*timeIdx=*/0);
+        }
         if (useVolumetricResidual) {
             // make the residual volume specific (i.e., make it incorrect mass per cubic
             // meter instead of total mass)
@@ -492,6 +493,7 @@ protected:
     {
         EvalVector tmp;
         EqVector tmp2;
+        EvalVector tmp2Der;// this is added to support calculation of storage term derivative focusTime!=0 for need for adjoint: always used if cachedStororage term si false
         RateVector sourceRate;
 
         tmp = 0.0;
@@ -572,24 +574,41 @@ protected:
             else {
                 // if the mass storage at the beginning of the time step is not cached,
                 // we re-calculate it from scratch.
-                tmp2 = 0.0;
-                asImp_().computeStorage(tmp2, elemCtx,  dofIdx, /*timeIdx=*/1);
-                Opm::Valgrind::CheckDefined(tmp2);
+                // if storage cache not enables we always use derivatives
+                tmp2Der = 0.0;
+                asImp_().computeStorage(tmp2Der, elemCtx,  dofIdx, /*timeIdx=*/1);
+                Opm::Valgrind::CheckDefined(tmp2Der);
             }
 
-            // Use the implicit Euler time discretization
-            for (unsigned eqIdx = 0; eqIdx < numEq; ++eqIdx) {
-                tmp[eqIdx] -= tmp2[eqIdx];
-                tmp[eqIdx] *= scvVolume / elemCtx.simulator().timeStepSize();
-
-                residual[dofIdx][eqIdx] += tmp[eqIdx];
+            if( !elemCtx.enableStorageCache() or elemCtx.focusTimeIndex()>0){
+                // Use the implicit Euler time discretization
+                for (unsigned eqIdx = 0; eqIdx < numEq; ++eqIdx) {
+                    if (elemCtx.focusTimeIndex() == 0) {
+                        tmp2Der[eqIdx] = Toolbox::value(tmp2Der[eqIdx]);
+                    } else {
+                        tmp[eqIdx] = Toolbox::value(tmp[eqIdx]);
+                    }
+                    tmp[eqIdx] -= tmp2Der[eqIdx];
+                    tmp[eqIdx] *= scvVolume / elemCtx.simulator().timeStepSize();
+                    residual[dofIdx][eqIdx] += tmp[eqIdx];
+                }
+            }else{
+                for (unsigned eqIdx = 0; eqIdx < numEq; ++eqIdx) {
+                    // cached storage should only be used in forward mode
+                    assert(elemCtx.focusTimeIndex()==0);
+                    tmp[eqIdx] -= tmp2[eqIdx];
+                    tmp[eqIdx] *= scvVolume / elemCtx.simulator().timeStepSize();
+                    residual[dofIdx][eqIdx] += tmp[eqIdx];
+                }
             }
 
             Opm::Valgrind::CheckDefined(residual[dofIdx]);
 
             // deal with the source term
-            asImp_().computeSource(sourceRate, elemCtx, dofIdx, /*timeIdx=*/0);
-
+            //assumes sourceterm is independent of prevois state
+            if(elemCtx.focusTimeIndex()==0){
+                asImp_().computeSource(sourceRate, elemCtx, dofIdx, /*timeIdx=*/0);
+            }
             // if the model uses extensive quantities in its storage term, and we use
             // automatic differention and current DOF is also not the one we currently
             // focus on, the storage term does not need any derivatives!