Changeset 25550

Timestamp:

09/09/20 20:01:09 (5 years ago)

Author:

Mathieu Morlighem

Message:

CHG: fixing AD checkpointing

File:

• issm/trunk-jpl/src/c/cores/transient_core.cpp (modified) (8 diffs)

issm/trunk-jpl/src/c/cores/transient_core.cpp

@@ -4 +4 @@
 
 #ifdef HAVE_CONFIG_H
-        #include <config.h>
+#include <config.h>
 #else
 #error "Cannot compile with HAVE_CONFIG_H symbol! run configure first!"
@@ -289 +289 @@
 
         std::vector<IssmDouble> time_all;
-   int                     Ysize = 0;
+        int                     Ysize = 0;
         CountDoublesFunctor   *hdl_countdoubles = NULL;
         RegisterInputFunctor  *hdl_regin        = NULL;
@@ -314 +314 @@
                 hdl_countdoubles = new CountDoublesFunctor();
                 femmodel->Marshall(hdl_countdoubles);
-      if(hdl_countdoubles->DoubleCount()>Ysize) Ysize= hdl_countdoubles->DoubleCount();
+                if(hdl_countdoubles->DoubleCount()>Ysize) Ysize= hdl_countdoubles->DoubleCount();
                 delete hdl_countdoubles;
 
@@ -336 +336 @@
         GetVectorFromControlInputsx(&X,&Xsize,femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,"value");
 
-   /*Initialize model state adjoint (Yb)*/
-   double *Yb  = xNewZeroInit<double>(Ysize);
-   int    *Yin = xNewZeroInit<int>(Ysize);
-
+        /*Initialize model state adjoint (Yb)*/
+        double *Yb  = xNewZeroInit<double>(Ysize);
+        int    *Yin = xNewZeroInit<int>(Ysize);
+
+        /*Get final Ysize*/
+        hdl_countdoubles = new CountDoublesFunctor();
+        femmodel->Marshall(hdl_countdoubles);
+        int Ysize_i= hdl_countdoubles->DoubleCount();
+        delete hdl_countdoubles;
 
         /*Start tracing*/
@@ -346 +351 @@
 
         /*Reverse dependent (f)*/
-   hdl_regin = new RegisterInputFunctor(Yin,Ysize);
-   femmodel->Marshall(hdl_regin);
-   delete hdl_regin;
+        hdl_regin = new RegisterInputFunctor(Yin,Ysize);
+        femmodel->Marshall(hdl_regin);
+        delete hdl_regin;
         for(int i=0; i < Xsize; i++) tape_codi.registerInput(X[i]);
 
         SetControlInputsFromVectorx(femmodel,X);
-   IssmDouble J = 0.;
-        for(int i=0;i<dependent_objects->Size();i++){
-                DependentObject* dep=(DependentObject*)dependent_objects->GetObjectByOffset(i);
-                J += dep->GetValue();
+        IssmDouble J = 0.;
+        if(0){
+                for(int i=0;i<dependent_objects->Size();i++){
+                        DependentObject* dep=(DependentObject*)dependent_objects->GetObjectByOffset(i);
+                        J += dep->GetValue();
+                }
+        }
+        else{
+                femmodel->IceVolumex(&J,false);
         }
         if(IssmComm::GetRank()==0) tape_codi.registerOutput(J);
@@ -363 +373 @@
         tape_codi.evaluate();
 
-   /*Initialize Xb and Yb*/
-   double *Xb  = xNewZeroInit<double>(Xsize);
-        for(int i=0;i<Xsize;i++) Xb[i] += X[i].gradient();
-
-   hdl_setadjoint = new SetAdjointFunctor(Yb,Ysize);
-   femmodel->Marshall(hdl_setadjoint);
-   delete hdl_setadjoint;
+        /*Initialize Xb and Yb*/
+        double *Xb  = xNewZeroInit<double>(Xsize);
+        for(int i=0;i<Xsize  ;i++) Xb[i] += X[i].gradient();
+        for(int i=0;i<Ysize_i;i++) Yb[i]  = tape_codi.gradient(Yin[i]);
 
         /*reverse loop for transient step (G)*/
-        for(int reverse_step = step;reverse_step>=0; reverse_step--){
+        for(int reverse_step = step;reverse_step>0; reverse_step--){
 
                 /*Restore model from this step*/
@@ -378 +385 @@
                 femmodel->RestartAD(reverse_step);
                 tape_codi.setActive();
+
+                /*Get new Ysize*/
+                hdl_countdoubles = new CountDoublesFunctor();
+                femmodel->Marshall(hdl_countdoubles);
+                int Ysize_i= hdl_countdoubles->DoubleCount();
+                delete hdl_countdoubles;
 
                 /*We need to store the CoDiPack identifier here, since y is overwritten.*/
@@ -408 +421 @@
 
                 /* here we access the gradient data via the stored identifiers.*/
-                for(int i=0; i<Ysize; i++) Yb[i]  = tape_codi.gradient(Yin[i]);
-                for(int i=0; i<Xsize; i++) Xb[i] += X[i].gradient();
+                //for(int i=0; i<Ysize_i; i++) if(tape_codi.gradient(Yin[i])!=0.) printf(" %g (%i)",tape_codi.gradient(Yin[i]),Yin[i]);
+                for(int i=0; i<Ysize_i; i++) Yb[i]  = tape_codi.gradient(Yin[i]);
+                for(int i=0; i<Xsize;   i++) Xb[i] += X[i].gradient();
         }