Changeset 25476

Timestamp:

08/25/20 21:01:51 (5 years ago)

Author:

Mathieu Morlighem

Message:

NEW: working on check pointing, not working just yet...

Location:

issm/trunk-jpl/src/c

Files:

• classes/FemModel.cpp (modified) (2 diffs)
• classes/FemModel.h (modified) (2 diffs)
• cores/transient_core.cpp (modified) (5 diffs)

issm/trunk-jpl/src/c/classes/FemModel.cpp

@@ -277 +277 @@
 }
 /*}}}*/
+void FemModel::CheckPointAD(int step){/*{{{*/
+
+        /*Get rank*/
+        int my_rank = IssmComm::GetRank();
+
+        /*Get string sizes*/
+        int rank_length = (my_rank == 0 ? 1 : int(log10(static_cast<double>(my_rank))+1));
+        int step_length = (step    == 0 ? 1 : int(log10(static_cast<double>(step))   +1));
+
+        /*Create restart file*/
+        char* restartfilename  = xNew<char>(strlen("AD_step_")+step_length+strlen("_rank_")+rank_length+strlen(".ckpt")+1);
+        sprintf(restartfilename,"%s%i%s%i%s","AD_step_",step,"_rank_",my_rank,".ckpt");
+        this->parameters->AddObject(new StringParam(RestartFileNameEnum,restartfilename));
+
+        /*Write files*/
+        this->CheckPoint();
+
+        /*Clean up and return*/
+        xDelete<char>(restartfilename);
+
+}/*}}}*/
 void FemModel::CleanUp(void){/*{{{*/
 
@@ -600 +621 @@
         xDelete<char>(restartfilename);
         xDelete<char>(femmodel_buffer);
+}/*}}}*/
+void FemModel::RestartAD(int step){ /*{{{*/
+
+        /*Get rank*/
+        int my_rank = IssmComm::GetRank();
+
+        /*Get string sizes*/
+        int rank_length = (my_rank == 0 ? 1 : int(log10(static_cast<double>(my_rank))+1));
+        int step_length = (step    == 0 ? 1 : int(log10(static_cast<double>(step))   +1));
+
+        /*Create restart file*/
+        char* restartfilename  = xNew<char>(strlen("AD_step_")+step_length+strlen("_rank_")+rank_length+strlen(".ckpt")+1);
+        sprintf(restartfilename,"%s%i%s%i%s","AD_step_",step,"_rank_",my_rank,".ckpt");
+        this->parameters->AddObject(new StringParam(RestartFileNameEnum,restartfilename));
+
+        /*Read files*/
+        this->Restart();
+
+        /*Clean up and return*/
+        xDelete<char>(restartfilename);
 }/*}}}*/
 void FemModel::SetCurrentConfiguration(int configuration_type,int analysis_type){/*{{{*/

issm/trunk-jpl/src/c/classes/FemModel.h

@@ -74 +74 @@
                 int  AnalysisIndex(int);
                 void CheckPoint(void);
+                void CheckPointAD(int step);
                 void CleanUp(void);
                 FemModel* copy();
@@ -82 +83 @@
                 void Marshall(char** pmarshalled_data, int* pmarshalled_data_size, int marshall_direction);
                 void Restart(void);
+                void RestartAD(int step);
                 void SetCurrentConfiguration(int configuration_type);
                 void SetCurrentConfiguration(int configuration_type,int analysis_type);

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

@@ -272 +272 @@
 
         /*parameters: */
-        IssmDouble finaltime,dt,yts;
-        bool       isoceancoupling,iscontrol,isautodiff,isslr;
-        int        timestepping;
-        int        output_frequency;
-        int        recording_frequency;
-        int        amr_frequency,amr_restart;
-
-        /*intermediary: */
-        int        step;
-        IssmDouble time;
-
-        /*first, figure out if there was a check point, if so, do a reset of the FemModel* femmodel structure. */
-        femmodel->parameters->FindParam(&recording_frequency,SettingsRecordingFrequencyEnum);
-        if(recording_frequency) femmodel->Restart();
+        IssmDouble output_value;
+        IssmDouble finaltime,dt,yts,time;
+        bool       isoceancoupling,isslr;
+        int        step,timestepping;
+        DataSet*   dependent_objects=NULL;
 
         /*then recover parameters common to all solutions*/
@@ -292 +283 @@
         femmodel->parameters->FindParam(&finaltime,TimesteppingFinalTimeEnum);
         femmodel->parameters->FindParam(&yts,ConstantsYtsEnum);
-        femmodel->parameters->FindParam(&output_frequency,SettingsOutputFrequencyEnum);
         femmodel->parameters->FindParam(&timestepping,TimesteppingTypeEnum);
         femmodel->parameters->FindParam(&isslr,TransientIsslrEnum);
-        femmodel->parameters->FindParam(&amr_frequency,TransientAmrFrequencyEnum);
-        femmodel->parameters->FindParam(&iscontrol,InversionIscontrolEnum);
-        femmodel->parameters->FindParam(&isautodiff,AutodiffIsautodiffEnum);
-
-        DataSet* dependent_objects=NULL;
-        if(iscontrol){
-                femmodel->parameters->FindParam(&dependent_objects,AutodiffDependentObjectsEnum);
-        }
-
+        femmodel->parameters->FindParam(&dependent_objects,AutodiffDependentObjectsEnum);
         if(isslr) sealevelrise_core_geometry(femmodel);
 
@@ -309 +291 @@
 
                 /*Time Increment*/
-                switch(timestepping){
-                        case AdaptiveTimesteppingEnum:
-                                femmodel->TimeAdaptx(&dt);
-                                if(time+dt>finaltime) dt=finaltime-time;
-                                femmodel->parameters->SetParam(dt,TimesteppingTimeStepEnum);
-                                break;
-                        case FixedTimesteppingEnum:
-                                femmodel->parameters->FindParam(&dt,TimesteppingTimeStepEnum);
-                                break;
-                        default:
-                                _error_("Time stepping \""<<EnumToStringx(timestepping)<<"\" not supported yet");
-                }
+                if(timestepping!=FixedTimesteppingEnum) _error_("not supported yet, but easy to handle...");
+                femmodel->parameters->FindParam(&dt,TimesteppingTimeStepEnum);
                 step+=1;
                 time+=dt;
@@ -327 +299 @@
                 femmodel->parameters->SetParam(false,SaveResultsEnum);
 
+                /*Store Model State at the beginning of the step*/
+                if(VerboseSolution()) _printf0_("   checkpointing model \n");
+                femmodel->CheckPointAD(step);
+
                 /*Run transient step!*/
                 transient_step(femmodel);
-
-                if(recording_frequency && (step%recording_frequency==0)){
-                        if(VerboseSolution()) _printf0_("   checkpointing model \n");
-                        femmodel->CheckPoint();
-                }
 
                 /*Go through our dependent variables, and compute the response:*/
                 for(int i=0;i<dependent_objects->Size();i++){
                         DependentObject* dep=(DependentObject*)dependent_objects->GetObjectByOffset(i);
-                        IssmDouble  output_value;
                         dep->Responsex(&output_value,femmodel);
                         dep->AddValue(output_value);
@@ -344 +314 @@
         }
 
-        if(!iscontrol) femmodel->RequestedDependentsx();
-        if(iscontrol) femmodel->parameters->SetParam(dependent_objects,AutodiffDependentObjectsEnum);
+        /*__________________________________________________________________________________*/
+
+        /*Get X (control) and Y (model state)*/
+        IssmDouble *X = NULL; int Xsize;
+        GetVectorFromControlInputsx(&X,&Xsize,femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,"value");
+        IssmDouble *Y = NULL; int Ysize;
+        _error_("don't know how to get Y and Ysize here...");
+
+        /*Initialize Xb, Yb and Yin*/
+        double *Xb = xNewZeroInit<double>(Xsize);
+        double *Yb  = xNewZeroInit<double>(Ysize);
+        double *Yin = xNewZeroInit<double>(Ysize);
+
+        /*Start tracing*/
+        auto& tape_codi = IssmDouble::getGlobalTape();
+        tape_codi.setActive();
+
+        /*Reverse dependent (f)*/
+        for(int i=0; i < Ysize; i++) tape_codi.registerInput(Y[i]);
+        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();
+        }
+
+        if(IssmComm::GetRank()==0) {
+                tape_codi.registerOutput(J);
+        }
+        tape_codi.setPassive();
+
+        J.gradient() = 1.0;
+        tape_codi.evaluate();
+
+        for(int i=0;i<Xsize;i++) Xb[i] += X[i].gradient();
+        for(int i=0;i<Ysize;i++) Yb[i]  = Y[i].gradient();
+
+        /*reverse loop for transient step (G)*/
+        for(int reverse_step = step;reverse_step>=0; reverse_step--){
+
+                /*Restore model from this step*/
+                tape_codi.reset();
+                femmodel->RestartAD(reverse_step);
+                tape_codi.setActive();
+
+                /*We need to store the CoDiPack identifier here, since y is overwritten.*/
+                for(int i=0; i<Ysize; i++) {
+                        tape_codi.registerInput(Y[i]);
+                        Yin[i] = Y[i].getGradientData();
+                }
+
+                /*Tell codipack that X is the independent*/
+                for(int i=0; i<Xsize; i++) tape_codi.registerInput(X[i]);
+                SetControlInputsFromVectorx(femmodel,X);
+
+                /*Get New state*/
+                transient_step(femmodel);
+                _error_("don't know how to get new Y again here...");
+                for(int i=0; i<Ysize; i++) tape_codi.registerOutput(Y[i]);
+
+                /*stop tracing*/
+                tape_codi.setPassive();
+
+                /*Reverse transient step (G)*/
+                /* Using y_b here to seed the next reverse iteration there y_b is always overwritten*/
+                for(int i=0; i<Ysize; i++) Y[i].gradient() = Yb[i]; 
+                tape_codi.evaluate();
+                /* here we access the gradient data via the stored identifiers.*/
+                for(int i=0; i<Ysize; i++) Xb[i]  = tape_codi.gradient(Yin[i]);
+                for(int i=0; i<Xsize; i++) Yb[i] += X[i].gradient();
+        }
+
+        /*Now set final gradient*/
+        IssmDouble* aG=xNew<IssmDouble>(Xsize);
+        for(int i=0;i<Xsize;i++){
+                aG[i] = reCast<IssmDouble>(Xb[i]);
+        }
+        ControlInputSetGradientx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,aG);
+        xDelete<IssmDouble>(aG);
+
+        xDelete<IssmDouble>(X);
+        xDelete<double>(Xb);
+        xDelete<IssmDouble>(Y);
+        xDelete<double>(Yb);
+        xDelete<double>(Yin);
 }/*}}}*/
 #endif