Changeset 25560

Timestamp:

09/11/20 09:53:19 (5 years ago)

Author:

Mathieu Morlighem

Message:

CHG: almost done with checkpointing AD

File:

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

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

@@ -276 +276 @@
         bool       isoceancoupling,isslr;
         int        step,timestepping;
+        int        recording_frequency,finalcheckpoint;
         DataSet*   dependent_objects=NULL;
 
@@ -286 +287 @@
         femmodel->parameters->FindParam(&isslr,TransientIsslrEnum);
         femmodel->parameters->FindParam(&dependent_objects,AutodiffDependentObjectsEnum);
+        femmodel->parameters->FindParam(&recording_frequency,SettingsRecordingFrequencyEnum); _assert_(recording_frequency>0);
         if(isslr) sealevelrise_core_geometry(femmodel);
 
@@ -298 +300 @@
 
                 /*Time Increment*/
-                if(timestepping!=FixedTimesteppingEnum) _error_("not supported yet, but easy to handle...");
-                femmodel->parameters->FindParam(&dt,TimesteppingTimeStepEnum);
+                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");
+                }
                 step+=1;
                 time+=dt;
@@ -307 +319 @@
                 time_all.push_back(time);
 
+                if(VerboseSolution()){
+                        _printf0_("iteration " << step << "/" << ceil((finaltime-time)/dt)+step << \
+                                                "  time [yr]: " <<std::fixed<<setprecision(2)<< time/yts << " (time step: " << dt/yts << ")\n");
+                }
+
                 /*Store Model State at the beginning of the step*/
-                if(VerboseSolution()) _printf0_("   checkpointing model \n");
-                femmodel->CheckPointAD(step);
+                if(step%recording_frequency==0){
+                        if(VerboseSolution()) _printf0_("   checkpointing model (step: "<<step<<")\n");
+                        femmodel->CheckPointAD(step);
+                        finalcheckpoint = step;
+                }
 
                 if(VerboseSolution()) _printf0_("   counting number of active variables\n");
@@ -328 +348 @@
         }
 
+        int finalstep = step;
         if(VerboseSolution()) _printf0_("   done with initial complete transient\n");
 
@@ -366 +387 @@
         }
         else{
-                femmodel->IceVolumex(&J,false);
+                femmodel->IceVolumeAboveFloatationx(&J,false);
         }
         if(IssmComm::GetRank()==0) tape_codi.registerOutput(J);
@@ -380 +401 @@
 
         /*reverse loop for transient step (G)*/
-        for(int reverse_step = step;reverse_step>0; reverse_step--){
+        for(int reverse_step = finalcheckpoint;reverse_step>0; reverse_step-=recording_frequency){
 
                 /*Restore model from this step*/
@@ -403 +424 @@
 
                 /*Get New state*/
-                transient_step(femmodel);
+                for(int ii=0;ii<recording_frequency;ii++){
+                        int        thisstep = reverse_step+ii;
+                        IssmDouble thistime = time_all[reverse_step+ii-1];
+                        _printf0_("step "<<thisstep<<" ("<<ii+1<<"/"<<recording_frequency<<") time [yr]: "<<std::fixed<<std::setprecision(2)<<thistime/yts<<"\n");
+                        transient_step(femmodel);
+                        if(thisstep==finalstep) break;
+                }
 
                 /*Register output*/
@@ -422 +449 @@
 
                 /* here we access the gradient data via the stored identifiers.*/
-                //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();