Changeset 25511


Ignore:
Timestamp:
09/01/20 16:33:34 (5 years ago)
Author:
Mathieu Morlighem
Message:

CHG: done with checkpointing... debugging time!

Location:
issm/trunk-jpl/src/c
Files:
2 edited

Legend:

Unmodified
Added
Removed

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

    r25497 r25511  
    316316        /*__________________________________________________________________________________*/
    317317
    318         /*Get X (control) and Y (model state)*/
     318        /*Get X (control)*/
    319319        IssmDouble *X = NULL; int Xsize;
    320320        GetVectorFromControlInputsx(&X,&Xsize,femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,"value");
    321         IssmDouble *Y = NULL; int Ysize;
    322         _error_("don't know how to get Y and Ysize here...");
     321
     322        /*Get Y (model state) size*/
     323        CountDoublesFunctor* marshallhandle1 = new CountDoublesFunctor();
     324        femmodel->Marshall(marshallhandle1);
     325        int Ysize = marshallhandle1->DoubleCount();
     326        delete marshallhandle1;
    323327
    324328        /*Initialize Xb, Yb and Yin*/
    325         double *Xb = xNewZeroInit<double>(Xsize);
     329        double *Xb  = xNewZeroInit<double>(Xsize);
    326330        double *Yb  = xNewZeroInit<double>(Ysize);
    327331        int    *Yin = xNewZeroInit<int>(Ysize);
    … …  
    332336
    333337        /*Reverse dependent (f)*/
    334         for(int i=0; i < Ysize; i++) tape_codi.registerInput(Y[i]);
     338        RegisterInputFunctor* marshallhandle2 = new RegisterInputFunctor(Yin);
     339        femmodel->Marshall(marshallhandle2);
     340        delete marshallhandle2;
    335341        for(int i=0; i < Xsize; i++) tape_codi.registerInput(X[i]);
    336342        SetControlInputsFromVectorx(femmodel,X);
    … …  
    351357
    352358        for(int i=0;i<Xsize;i++) Xb[i] += X[i].gradient();
    353         for(int i=0;i<Ysize;i++) Yb[i]  = Y[i].gradient();
     359        InitAdjointFunctor* marshallhandle3 = new InitAdjointFunctor(Yb);
     360        femmodel->Marshall(marshallhandle3);
     361        delete marshallhandle3;
    354362
    355363        /*reverse loop for transient step (G)*/
    … …  
    362370
    363371                /*We need to store the CoDiPack identifier here, since y is overwritten.*/
    364                 for(int i=0; i<Ysize; i++) {
    365                         tape_codi.registerInput(Y[i]);
    366                         Yin[i] = Y[i].getGradientData();
    367                 }
     372                RegisterInputFunctor* marshallhandle4 = new RegisterInputFunctor(Yin);
     373                femmodel->Marshall(marshallhandle4);
     374                delete marshallhandle4;
    368375
    369376                /*Tell codipack that X is the independent*/
    … …  
    373380                /*Get New state*/
    374381                transient_step(femmodel);
    375                 _error_("don't know how to get new Y again here...");
    376                 for(int i=0; i<Ysize; i++) tape_codi.registerOutput(Y[i]);
     382
     383                /*Register output*/
     384                RegisterOutputFunctor* marshallhandle5 = new RegisterOutputFunctor();
     385                femmodel->Marshall(marshallhandle5);
     386                delete marshallhandle5;
    377387
    378388                /*stop tracing*/
    … …  
    381391                /*Reverse transient step (G)*/
    382392                /* Using y_b here to seed the next reverse iteration there y_b is always overwritten*/
    383                 for(int i=0; i<Ysize; i++) Y[i].gradient() = Yb[i];
     393                SetAdjointFunctor* marshallhandle6 = new SetAdjointFunctor(Yb);
     394                femmodel->Marshall(marshallhandle6);
     395                delete marshallhandle6;
     396
    384397                tape_codi.evaluate();
     398
    385399                /* here we access the gradient data via the stored identifiers.*/
    386                 for(int i=0; i<Ysize; i++) Xb[i]  = tape_codi.gradient(Yin[i]);
    387                 for(int i=0; i<Xsize; i++) Yb[i] += X[i].gradient();
     400                for(int i=0; i<Ysize; i++) Yb[i]  = tape_codi.gradient(Yin[i]);
     401                for(int i=0; i<Xsize; i++) Xb[i] += X[i].gradient();
    388402        }
    389403
    … …  
    398412        xDelete<IssmDouble>(X);
    399413        xDelete<double>(Xb);
    400         xDelete<IssmDouble>(Y);
    401414        xDelete<double>(Yb);
    402415        xDelete<int>(Yin);

  • issm/trunk-jpl/src/c/shared/io/Marshalling/Marshalling.h

    r25510 r25511  
    2020        AD_REGISTERINPUT,
    2121        AD_REGISTEROUTPUT,
     22        AD_INITADJOINT,
    2223        AD_SETADJOINT,
    2324#endif
    … …  
    152153                        this->tape_codi    = &(IssmDouble::getGlobalTape());
    153154                }
    154                 int DoubleCount(void){return this->double_count;};
    155155                void Echo(void){
    156156                        printf("RegisterInputFunctor Echo:\n");
    … …  
    185185
    186186        public:
    187                 RegisterOutputFunctor(void) : MarshallHandle(AD_REGISTERINPUT){
     187                RegisterOutputFunctor(void) : MarshallHandle(AD_REGISTEROUTPUT){
    188188                        this->double_count = 0;
    189189                        this->tape_codi    = &(IssmDouble::getGlobalTape());
    190190                }
    191                 int DoubleCount(void){return this->double_count;};
    192191                void Echo(void){
    193192                        printf("RegisterOutputFunctor Echo:\n");
    … …  
    208207                                for(int i=0;i<size;i++){
    209208                                        this->tape_codi->registerOutput(value[i]);
     209                                        this->double_count++;
     210                                }
     211                        }
     212                }
     213}; /*}}}*/
     214class InitAdjointFunctor:     public MarshallHandle{ /*{{{*/
     215
     216        private:
     217                int                   double_count;
     218                IssmDouble::TapeType* tape_codi;
     219                double*               adjoint;
     220
     221        public:
     222                InitAdjointFunctor(double* adjoint_in) : MarshallHandle(AD_INITADJOINT){
     223                        this->double_count = 0;
     224                        this->tape_codi    = &(IssmDouble::getGlobalTape());
     225                        this->adjoint      = adjoint_in;
     226                }
     227                void Echo(void){
     228                        printf("InitAdjointFunctor Echo:\n");
     229                        printf("   double_count: %i\n",double_count);
     230                }
     231                template<typename T> void call(T & value){
     232                        /*General case: do nothing*/
     233                }
     234                template<typename T> void call(T* & value,int size){
     235                        /*General case: do nothing*/
     236                }
     237                void call(IssmDouble value){
     238                        value.gradient() = adjoint[this->double_count];
     239                        this->double_count++;
     240                }
     241                void call(IssmDouble* value,int size){
     242                        if(value){
     243                                for(int i=0;i<size;i++){
     244                                        value[i].gradient() = adjoint[this->double_count];
     245                                        this->double_count++;
     246                                }
     247                        }
     248                }
     249}; /*}}}*/
     250class SetAdjointFunctor:      public MarshallHandle{ /*{{{*/
     251
     252        private:
     253                int                   double_count;
     254                IssmDouble::TapeType* tape_codi;
     255                double*               adjoint;
     256
     257        public:
     258                SetAdjointFunctor(double* adjoint_in) : MarshallHandle(AD_SETADJOINT){
     259                        this->double_count = 0;
     260                        this->tape_codi    = &(IssmDouble::getGlobalTape());
     261                        this->adjoint      = adjoint_in;
     262                }
     263                void Echo(void){
     264                        printf("SetAdjointFunctor Echo:\n");
     265                        printf("   double_count: %i\n",double_count);
     266                }
     267                template<typename T> void call(T & value){
     268                        /*General case: do nothing*/
     269                }
     270                template<typename T> void call(T* & value,int size){
     271                        /*General case: do nothing*/
     272                }
     273                void call(IssmDouble value){
     274                        value.gradient() = this->adjoint[this->double_count];
     275                        this->double_count++;
     276                }
     277                void call(IssmDouble* value,int size){
     278                        if(value){
     279                                for(int i=0;i<size;i++){
     280                                        value[i].gradient() = this->adjoint[this->double_count];
    210281                                        this->double_count++;
    211282                                }
    … …  
    221292                case MARSHALLING_SIZE: {SizeCheckpointFunctor*  temp = xDynamicCast<SizeCheckpointFunctor*>(this);  temp->call(value); break;}
    222293#ifdef _HAVE_CODIPACK_
    223                 case AD_COUNTDOUBLES:  {CountDoublesFunctor*    temp = xDynamicCast<CountDoublesFunctor*>(this);    temp->call(value); break;}
    224                 case AD_REGISTERINPUT: {RegisterInputFunctor*   temp = xDynamicCast<RegisterInputFunctor*>(this);   temp->call(value); break;}
    225                 case AD_REGISTEROUTPUT:{RegisterOutputFunctor*  temp = xDynamicCast<RegisterOutputFunctor*>(this);  temp->call(value); break;}
    226                 //case AD_SETADJOINT:    {SetAdjointFunction*     temp = xDynamicCast<SetAdjointFunction*>(this);     temp->call(value); break;}
     294                case AD_COUNTDOUBLES:  {CountDoublesFunctor*   temp = xDynamicCast<CountDoublesFunctor*>(this);    temp->call(value); break;}
     295                case AD_REGISTERINPUT: {RegisterInputFunctor*  temp = xDynamicCast<RegisterInputFunctor*>(this);   temp->call(value); break;}
     296                case AD_REGISTEROUTPUT:{RegisterOutputFunctor* temp = xDynamicCast<RegisterOutputFunctor*>(this);  temp->call(value); break;}
     297                case AD_INITADJOINT:   {InitAdjointFunctor*    temp = xDynamicCast<InitAdjointFunctor*>(this);     temp->call(value); break;}
     298                case AD_SETADJOINT:    {SetAdjointFunctor*     temp = xDynamicCast<SetAdjointFunctor*>(this);      temp->call(value); break;}
    227299#endif
    228300                default: _error_("Operation "<<OperationNumber()<<" not supported yet");
    … …  
    235307                case MARSHALLING_SIZE: {SizeCheckpointFunctor*  temp = xDynamicCast<SizeCheckpointFunctor*>(this);  temp->call(value,size); break;}
    236308#ifdef _HAVE_CODIPACK_
    237                 case AD_COUNTDOUBLES:  {CountDoublesFunctor*    temp = xDynamicCast<CountDoublesFunctor*>(this);    temp->call(value,size); break;}
    238                 case AD_REGISTERINPUT: {RegisterInputFunctor*   temp = xDynamicCast<RegisterInputFunctor*>(this);   temp->call(value,size); break;}
    239                 case AD_REGISTEROUTPUT:{RegisterOutputFunctor*  temp = xDynamicCast<RegisterOutputFunctor*>(this);  temp->call(value,size); break;}
    240                 //case AD_SETADJOINT:    {SetAdjointFunction*     temp = xDynamicCast<SetAdjointFunction*>(this);     temp->call(value,size); break;}
     309                case AD_COUNTDOUBLES:  {CountDoublesFunctor*   temp = xDynamicCast<CountDoublesFunctor*>(this);    temp->call(value,size); break;}
     310                case AD_REGISTERINPUT: {RegisterInputFunctor*  temp = xDynamicCast<RegisterInputFunctor*>(this);   temp->call(value,size); break;}
     311                case AD_REGISTEROUTPUT:{RegisterOutputFunctor* temp = xDynamicCast<RegisterOutputFunctor*>(this);  temp->call(value,size); break;}
     312                case AD_INITADJOINT:   {InitAdjointFunctor*    temp = xDynamicCast<InitAdjointFunctor*>(this);     temp->call(value,size); break;}
     313                case AD_SETADJOINT:    {SetAdjointFunctor*     temp = xDynamicCast<SetAdjointFunctor*>(this);      temp->call(value,size); break;}
    241314#endif
    242315                default: _error_("Operation "<<OperationNumber() <<" not supported yet");
Note: See TracChangeset for help on using the changeset viewer.