Changeset 25536

Timestamp:

09/04/20 20:15:03 (5 years ago)

Author:

Mathieu Morlighem

Message:

CHG: working on AD transient checkpointing

File:

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

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

@@ -288 +288 @@
         if(isslr) sealevelrise_core_geometry(femmodel);
 
+        std::vector<IssmDouble> time_all;
+   int                     Ysize = 0;
+        CountDoublesFunctor   *hdl_countdoubles = NULL;
+        RegisterInputFunctor  *hdl_regin        = NULL;
+        RegisterOutputFunctor *hdl_regout       = NULL;
+        SetAdjointFunctor     *hdl_setadjoint   = NULL;
+
         while(time < finaltime - (yts*DBL_EPSILON)){ //make sure we run up to finaltime.
 
@@ -298 +305 @@
                 femmodel->parameters->SetParam(step,StepEnum);
                 femmodel->parameters->SetParam(false,SaveResultsEnum);
+                time_all.push_back(time);
 
                 /*Store Model State at the beginning of the step*/
                 if(VerboseSolution()) _printf0_("   checkpointing model \n");
                 femmodel->CheckPointAD(step);
+
+                if(VerboseSolution()) _printf0_("   counting number of active variables\n");
+                hdl_countdoubles = new CountDoublesFunctor();
+                femmodel->Marshall(hdl_countdoubles);
+      if(hdl_countdoubles->DoubleCount()>Ysize) Ysize= hdl_countdoubles->DoubleCount();
+                delete hdl_countdoubles;
 
                 /*Run transient step!*/
@@ -314 +328 @@
         }
 
+        if(VerboseSolution()) _printf0_("   done with initial complete transient\n");
+
         /*__________________________________________________________________________________*/
 
@@ -320 +336 @@
         GetVectorFromControlInputsx(&X,&Xsize,femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,"value");
 
-        /*Get Y (model state) size*/
-        CountDoublesFunctor* marshallhandle1 = new CountDoublesFunctor();
-        femmodel->Marshall(marshallhandle1);
-        int Ysize = marshallhandle1->DoubleCount();
-        delete marshallhandle1;
-
-        /*Initialize Xb, Yb and Yin*/
-        double *Xb  = xNewZeroInit<double>(Xsize);
-        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);
+
 
         /*Start tracing*/
@@ -336 +346 @@
 
         /*Reverse dependent (f)*/
-        RegisterInputFunctor* marshallhandle2 = new RegisterInputFunctor(Yin);
-        femmodel->Marshall(marshallhandle2);
-        delete marshallhandle2;
+   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.;
+   IssmDouble J = 0.;
         for(int i=0;i<dependent_objects->Size();i++){
                 DependentObject* dep=(DependentObject*)dependent_objects->GetObjectByOffset(i);
                 J += dep->GetValue();
         }
-
-        if(IssmComm::GetRank()==0) {
-                tape_codi.registerOutput(J);
-        }
+        if(IssmComm::GetRank()==0) tape_codi.registerOutput(J);
+
         tape_codi.setPassive();
-
         J.gradient() = 1.0;
         tape_codi.evaluate();
 
+   /*Initialize Xb and Yb*/
+   double *Xb  = xNewZeroInit<double>(Xsize);
         for(int i=0;i<Xsize;i++) Xb[i] += X[i].gradient();
-        SetAdjointFunctor* marshallhandle3 = new SetAdjointFunctor(Yb);
-        femmodel->Marshall(marshallhandle3);
-        delete marshallhandle3;
+
+   hdl_setadjoint = new SetAdjointFunctor(Yb,Ysize);
+   femmodel->Marshall(hdl_setadjoint);
+   delete hdl_setadjoint;
 
         /*reverse loop for transient step (G)*/
@@ -370 +380 @@
 
                 /*We need to store the CoDiPack identifier here, since y is overwritten.*/
-                RegisterInputFunctor* marshallhandle4 = new RegisterInputFunctor(Yin);
-                femmodel->Marshall(marshallhandle4);
-                delete marshallhandle4;
+                hdl_regin = new RegisterInputFunctor(Yin,Ysize);
+                femmodel->Marshall(hdl_regin);
+                delete hdl_regin;
 
                 /*Tell codipack that X is the independent*/
@@ -382 +392 @@
 
                 /*Register output*/
-                RegisterOutputFunctor* marshallhandle5 = new RegisterOutputFunctor();
-                femmodel->Marshall(marshallhandle5);
-                delete marshallhandle5;
+                hdl_regout = new RegisterOutputFunctor();
+                femmodel->Marshall(hdl_regout);
+                delete hdl_regout;
 
                 /*stop tracing*/
@@ -391 +401 @@
                 /*Reverse transient step (G)*/
                 /* Using y_b here to seed the next reverse iteration there y_b is always overwritten*/
-                SetAdjointFunctor* marshallhandle6 = new SetAdjointFunctor(Yb);
-                femmodel->Marshall(marshallhandle6);
-                delete marshallhandle6;
+                hdl_setadjoint = new SetAdjointFunctor(Yb,Ysize);
+                femmodel->Marshall(hdl_setadjoint);
+                delete hdl_setadjoint;
 
                 tape_codi.evaluate();