Changeset 26889


Ignore:
Timestamp:
02/16/22 10:39:29 (3 years ago)
Author:
Mathieu Morlighem
Message:

CHG: making checkpointing work with multiple cost functions

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

Legend:

Unmodified
Added
Removed

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

    r26623 r26889  
    1515#include <sstream> // for output of the CoDiPack tape
    1616#include <fenv.h>
    17 void transient_ad(FemModel* femmodel);
     17double transient_ad(FemModel* femmodel, double* G,double* Jlist);
    1818#endif
    1919
    … …  
    222222        int    *N = NULL;
    223223        int    *control_enum    = NULL;
     224        int     checkpoint_frequency;
    224225        femmodel->parameters->FindParam(&num_responses,InversionNumCostFunctionsEnum);
    225226        femmodel->parameters->FindParam(&num_controls,InversionNumControlParametersEnum);
    … …  
    228229        femmodel->parameters->FindParam(&M,NULL,ControlInputSizeMEnum);
    229230        femmodel->parameters->FindParam(&control_enum,NULL,InversionControlParametersEnum);
     231        femmodel->parameters->FindParam(&checkpoint_frequency,SettingsCheckpointFrequencyEnum);
    230232
    231233        /*Constrain input vector and update controls*/
    … …  
    236238
    237239        int offset = 0;
    238         for (int c=0;c<num_controls;c++){
     240        for(int c=0;c<num_controls;c++){
    239241                for(int i=0;i<M[c]*N[c];i++){
    240242                        int index = offset+i;
    … …  
    246248        }
    247249
    248         /*Start Tracing*/
    249         simul_starttrace(femmodel);
    250         /*Set X as our new control input and as INDEPENDENT*/
     250        /*Special case: do we need to run AD with checkpointing?*/
     251        #ifdef _HAVE_CODIPACK_
     252        if(checkpoint_frequency && solution_type == TransientSolutionEnum){
     253                SetControlInputsFromVectorx(femmodel,X);
     254                *pf = transient_ad(femmodel, G, &Jlist[(*Jlisti)*JlistN]);
     255        }
     256        else
     257        #endif
     258          {
     259
     260                /*Start Tracing*/
     261                simul_starttrace(femmodel);
     262                /*Set X as our new control input and as INDEPENDENT*/
    251263#ifdef _HAVE_AD_
    252         IssmDouble* aX=xNew<IssmDouble>(intn,"t");
     264                IssmDouble* aX=xNew<IssmDouble>(intn,"t");
    253265#else
    254         IssmDouble* aX=xNew<IssmDouble>(intn);
     266                IssmDouble* aX=xNew<IssmDouble>(intn);
    255267#endif
    256268
    257         #if defined(_HAVE_ADOLC_)
    258         if(my_rank==0){
    259                 for(int i=0;i<intn;i++){
    260                         aX[i]<<=X[i];
    261                 }
    262         }
    263         #elif defined(_HAVE_CODIPACK_)
    264 
    265         /*Get tape*/
    266         #if _CODIPACK_MAJOR_==2
    267         auto& tape_codi = IssmDouble::getTape();
    268         #elif _CODIPACK_MAJOR_==1
    269         auto& tape_codi = IssmDouble::getGlobalTape();
    270         #else
    271         #error "_CODIPACK_MAJOR_ not supported"
    272         #endif
    273 
    274         codi_global.input_indices.clear();
    275         if(my_rank==0){
    276                 for (int i=0;i<intn;i++) {
    277                         aX[i]=X[i];
    278                         tape_codi.registerInput(aX[i]);
    279                         #if _CODIPACK_MAJOR_==2
    280                         codi_global.input_indices.push_back(aX[i].getIdentifier());
    281                         #elif _CODIPACK_MAJOR_==1
    282                         codi_global.input_indices.push_back(aX[i].getGradientData());
    283                         #else
    284                         #error "_CODIPACK_MAJOR_ not supported"
    285                         #endif
    286 
    287                 }
    288         }
    289         #else
    290         _error_("not suppoted");
    291         #endif
    292 
    293         ISSM_MPI_Bcast(aX,intn,ISSM_MPI_DOUBLE,0,IssmComm::GetComm());
    294         SetControlInputsFromVectorx(femmodel,aX);
    295         xDelete<IssmDouble>(aX);
    296 
    297         /*Compute solution (forward)*/
    298         void (*solutioncore)(FemModel*)=NULL;
    299         CorePointerFromSolutionEnum(&solutioncore,femmodel->parameters,solution_type);
    300         solutioncore(femmodel);
    301 
    302         /*Get Dependents*/
    303         IssmDouble   output_value;
    304         int          num_dependents;
    305         IssmPDouble *dependents;
    306         IssmDouble   J = 0.;
    307         DataSet     *dependent_objects = ((DataSetParam*)femmodel->parameters->FindParamObject(AutodiffDependentObjectsEnum))->value;
    308         femmodel->parameters->FindParam(&num_dependents,AutodiffNumDependentsEnum);
    309 
    310         /*Go through our dependent variables, and compute the response:*/
    311         dependents=xNew<IssmPDouble>(num_dependents);
    312         #if defined(_HAVE_CODIPACK_)
    313         codi_global.output_indices.clear();
    314         #endif
    315         int i=-1;
    316         for(Object* & object:dependent_objects->objects){
    317                 i++;
    318                 DependentObject* dep=xDynamicCast<DependentObject*>(object);
    319                 if(solution_type==TransientSolutionEnum) output_value = dep->GetValue();
    320                 if(solution_type!=TransientSolutionEnum) dep->Responsex(&output_value,femmodel);
    321                 if(my_rank==0) {
     269                #if defined(_HAVE_ADOLC_)
     270                if(my_rank==0){
     271                        for(int i=0;i<intn;i++){
     272                                aX[i]<<=X[i];
     273                        }
     274                }
     275                #elif defined(_HAVE_CODIPACK_)
     276
     277                /*Get tape*/
     278                #if _CODIPACK_MAJOR_==2
     279                auto& tape_codi = IssmDouble::getTape();
     280                #elif _CODIPACK_MAJOR_==1
     281                auto& tape_codi = IssmDouble::getGlobalTape();
     282                #else
     283                #error "_CODIPACK_MAJOR_ not supported"
     284                #endif
     285
     286                codi_global.input_indices.clear();
     287                if(my_rank==0){
     288                        for (int i=0;i<intn;i++) {
     289                                aX[i]=X[i];
     290                                tape_codi.registerInput(aX[i]);
     291                                #if _CODIPACK_MAJOR_==2
     292                                codi_global.input_indices.push_back(aX[i].getIdentifier());
     293                                #elif _CODIPACK_MAJOR_==1
     294                                codi_global.input_indices.push_back(aX[i].getGradientData());
     295                                #else
     296                                #error "_CODIPACK_MAJOR_ not supported"
     297                                #endif
     298
     299                        }
     300                }
     301                #else
     302                _error_("not suppoted");
     303                #endif
     304
     305                ISSM_MPI_Bcast(aX,intn,ISSM_MPI_DOUBLE,0,IssmComm::GetComm());
     306                SetControlInputsFromVectorx(femmodel,aX);
     307                xDelete<IssmDouble>(aX);
     308
     309                /*Compute solution (forward)*/
     310                void (*solutioncore)(FemModel*)=NULL;
     311                CorePointerFromSolutionEnum(&solutioncore,femmodel->parameters,solution_type);
     312                solutioncore(femmodel);
     313
     314                /*Get Dependents*/
     315                IssmDouble   output_value;
     316                int          num_dependents;
     317                IssmPDouble *dependents;
     318                IssmDouble   J = 0.;
     319                DataSet     *dependent_objects = ((DataSetParam*)femmodel->parameters->FindParamObject(AutodiffDependentObjectsEnum))->value;
     320                femmodel->parameters->FindParam(&num_dependents,AutodiffNumDependentsEnum);
     321
     322                /*Go through our dependent variables, and compute the response:*/
     323                dependents=xNew<IssmPDouble>(num_dependents);
     324                #if defined(_HAVE_CODIPACK_)
     325                codi_global.output_indices.clear();
     326                #endif
     327                int i=-1;
     328                IssmDouble TEMP = 0.;
     329                for(Object* & object:dependent_objects->objects){
     330                        i++;
     331                        DependentObject* dep=xDynamicCast<DependentObject*>(object);
     332                        if(solution_type==TransientSolutionEnum) output_value = dep->GetValue();
     333                        if(solution_type!=TransientSolutionEnum) dep->Responsex(&output_value,femmodel);
    322334
    323335                        #if defined(_HAVE_CODIPACK_)
    … …  
    340352                        J+=output_value;
    341353                }
    342         }
    343 
    344         /*Turning off trace tape*/
    345         simul_stoptrace();
    346 
    347         /*intermediary: */
    348         int          num_independents=intn;
    349         IssmPDouble *aWeightVector=NULL;
    350         IssmPDouble *weightVectorTimesJac=NULL;
    351         IssmPDouble *totalgradient=xNewZeroInit<IssmPDouble>(num_independents);
    352 
    353         /*if no dependents, no point in running a driver: */
    354         if(!(num_dependents*num_independents)) _error_("this is not allowed");
    355 
    356         /*for adolc to run in parallel, we 0 out on rank~=0. But we still keep track of num_dependents:*/
    357         int num_dependents_old   = num_dependents;
    358         int num_independents_old = num_independents;
    359 
    360         #if defined(_HAVE_ADOLC_)
    361         /*Get gradient for ADOLC {{{*/
    362         if(my_rank!=0){
    363                 num_dependents   = 0;
    364                 num_independents = 0;
    365         }
    366 
    367         /*get the EDF pointer:*/
    368         ext_diff_fct *anEDF_for_solverx_p=xDynamicCast<GenericParam<Adolc_edf> * >(femmodel->parameters->FindParamObject(AdolcParamEnum))->GetParameterValue().myEDF_for_solverx_p;
    369 
    370         /* these are always needed regardless of the interpreter */
    371         anEDF_for_solverx_p->dp_x=xNew<double>(anEDF_for_solverx_p->max_n);
    372         anEDF_for_solverx_p->dp_y=xNew<double>(anEDF_for_solverx_p->max_m);
    373 
    374         /* Ok, now we are going to call the fos_reverse in a loop on the index, from 0 to num_dependents, so
    375          * as to generate num_dependents gradients: */
    376         for(int aDepIndex=0;aDepIndex<num_dependents_old;aDepIndex++){
    377 
    378                 /*initialize direction index in the weights vector: */
    379                 aWeightVector=xNewZeroInit<IssmPDouble>(num_dependents);
    380                 if (my_rank==0) aWeightVector[aDepIndex]=1.;
    381 
    382                 /*initialize output gradient: */
    383                 weightVectorTimesJac=xNew<IssmPDouble>(num_independents);
    384 
    385                 /*set the forward method function pointer: */
    386                 #ifdef _HAVE_GSL_
    387                 anEDF_for_solverx_p->fos_reverse=EDF_fos_reverse_for_solverx;
    388                 #endif
    389                 #ifdef _HAVE_MUMPS_
    390                 anEDF_for_solverx_p->fos_reverse_iArr=fos_reverse_mumpsSolveEDF;
    391                 #endif
    392 
    393                 anEDF_for_solverx_p->dp_U=xNew<IssmPDouble>(anEDF_for_solverx_p->max_m);
    394                 anEDF_for_solverx_p->dp_Z=xNew<IssmPDouble>(anEDF_for_solverx_p->max_n);
    395 
    396                 /*call driver: */
    397                 fos_reverse(my_rank,num_dependents,num_independents, aWeightVector, weightVectorTimesJac );
    398 
    399                 /*Add to totalgradient: */
    400                 if(my_rank==0) for(int i=0;i<num_independents;i++) {
    401                         totalgradient[i]+=weightVectorTimesJac[i];
    402                 }
    403 
    404                 /*free resources :*/
    405                 xDelete(weightVectorTimesJac);
    406                 xDelete(aWeightVector);
    407         }
    408         /*}}}*/
    409         #elif defined(_HAVE_CODIPACK_)
    410         /*Get gradient for CoDiPack{{{*/
    411         if(VerboseAutodiff())_printf0_("   CoDiPack fos_reverse\n");
    412 
    413         /* call the fos_reverse in a loop on the index, from 0 to num_dependents, so
    414          * as to generate num_dependents gradients: */
    415         for(int dep_index=0;dep_index<num_dependents_old;dep_index++){
    416 
    417                 /*initialize direction index in the weights vector: */
    418                 if(my_rank==0){
    419                         if(dep_index<0 || dep_index>=num_dependents || codi_global.output_indices.size() <= dep_index){
    420                                 _error_("index value for dependent index should be in [0,num_dependents-1]");
     354
     355                /*Turning off trace tape*/
     356                simul_stoptrace();
     357
     358                /*intermediary: */
     359                int          num_independents=intn;
     360                IssmPDouble *aWeightVector=NULL;
     361                IssmPDouble *weightVectorTimesJac=NULL;
     362                IssmPDouble *totalgradient=xNewZeroInit<IssmPDouble>(num_independents);
     363
     364                /*if no dependents, no point in running a driver: */
     365                if(!(num_dependents*num_independents)) _error_("this is not allowed");
     366
     367                /*for adolc to run in parallel, we 0 out on rank~=0. But we still keep track of num_dependents:*/
     368                int num_dependents_old   = num_dependents;
     369                int num_independents_old = num_independents;
     370
     371                #if defined(_HAVE_ADOLC_)
     372                /*Get gradient for ADOLC {{{*/
     373                if(my_rank!=0){
     374                        num_dependents   = 0;
     375                        num_independents = 0;
     376                }
     377
     378                /*get the EDF pointer:*/
     379                ext_diff_fct *anEDF_for_solverx_p=xDynamicCast<GenericParam<Adolc_edf> * >(femmodel->parameters->FindParamObject(AdolcParamEnum))->GetParameterValue().myEDF_for_solverx_p;
     380
     381                /* these are always needed regardless of the interpreter */
     382                anEDF_for_solverx_p->dp_x=xNew<double>(anEDF_for_solverx_p->max_n);
     383                anEDF_for_solverx_p->dp_y=xNew<double>(anEDF_for_solverx_p->max_m);
     384
     385                /* Ok, now we are going to call the fos_reverse in a loop on the index, from 0 to num_dependents, so
     386                 * as to generate num_dependents gradients: */
     387                for(int aDepIndex=0;aDepIndex<num_dependents_old;aDepIndex++){
     388
     389                        /*initialize direction index in the weights vector: */
     390                        aWeightVector=xNewZeroInit<IssmPDouble>(num_dependents);
     391                        if (my_rank==0) aWeightVector[aDepIndex]=1.;
     392
     393                        /*initialize output gradient: */
     394                        weightVectorTimesJac=xNew<IssmPDouble>(num_independents);
     395
     396                        /*set the forward method function pointer: */
     397                        #ifdef _HAVE_GSL_
     398                        anEDF_for_solverx_p->fos_reverse=EDF_fos_reverse_for_solverx;
     399                        #endif
     400                        #ifdef _HAVE_MUMPS_
     401                        anEDF_for_solverx_p->fos_reverse_iArr=fos_reverse_mumpsSolveEDF;
     402                        #endif
     403
     404                        anEDF_for_solverx_p->dp_U=xNew<IssmPDouble>(anEDF_for_solverx_p->max_m);
     405                        anEDF_for_solverx_p->dp_Z=xNew<IssmPDouble>(anEDF_for_solverx_p->max_n);
     406
     407                        /*call driver: */
     408                        fos_reverse(my_rank,num_dependents,num_independents, aWeightVector, weightVectorTimesJac );
     409
     410                        /*Add to totalgradient: */
     411                        if(my_rank==0) for(int i=0;i<num_independents;i++) {
     412                                totalgradient[i]+=weightVectorTimesJac[i];
    421413                        }
    422                         tape_codi.setGradient(codi_global.output_indices[dep_index],1.0);
    423                 }
    424                 //feclearexcept(FE_ALL_EXCEPT);
    425                 //feenableexcept(FE_DIVBYZERO | FE_INVALID | FE_OVERFLOW);
    426                 tape_codi.evaluate();
    427 
    428                 /*Get gradient for this dependent */
    429                 weightVectorTimesJac=xNew<IssmPDouble>(num_independents);
    430                 auto in_size = codi_global.input_indices.size();
    431                 for(size_t i = 0; i < in_size; ++i){
    432                         _assert_(i<num_independents);
    433                         weightVectorTimesJac[i] = tape_codi.getGradient(codi_global.input_indices[i]);
    434                 }
    435                 if(my_rank==0) for(int i=0;i<num_independents;i++){
    436                         totalgradient[i]+=weightVectorTimesJac[i];
    437                 }
    438                 xDelete(weightVectorTimesJac);
    439         }
    440 
    441         /*Clear tape*/
    442         tape_codi.reset();
    443         /*}}}*/
    444         #else
    445         _error_("not suppoted");
    446         #endif
    447 
    448         /*Broadcast gradient to other ranks*/
    449         ISSM_MPI_Bcast(totalgradient,num_independents_old,ISSM_MPI_PDOUBLE,0,IssmComm::GetComm());
    450         /*Check size of Jlist to avoid crashes*/
    451         _assert_((*Jlisti)<JlistM);
    452         _assert_(JlistN==num_responses+1);
    453 
    454         /*Compute objective function and broadcast it to other cpus*/
    455         *pf = reCast<double>(J);
    456         ISSM_MPI_Bcast(pf,1,ISSM_MPI_PDOUBLE,0,IssmComm::GetComm());
    457         _printf0_("f(x) = "<<setw(12)<<setprecision(7)<<*pf<<"  |  ");
    458 
    459         /*Record cost function values and delete Jtemp*/
    460         for(int i=0;i<num_responses;i++) Jlist[(*Jlisti)*JlistN+i] = dependents[i];
    461         Jlist[(*Jlisti)*JlistN+num_responses] = reCast<IssmPDouble>(J);
    462 
    463         if(*indic==0){
    464                 /*dry run, no gradient required*/
    465 
    466                 /*Retrieve objective functions independently*/
    467                 _printf0_("            N/A |\n");
    468                 for(int i=0;i<num_responses;i++) _printf0_(" "<<setw(12)<<setprecision(7)<<Jlist[(*Jlisti)*JlistN+i]);
    469                 _printf0_("\n");
    470 
    471                 *Jlisti = (*Jlisti) +1;
    472                 xDelete<double>(XU);
    473                 xDelete<double>(XL);
    474                 return;
    475         }
    476 
    477         /*Compute gradient*/
    478         for(long i=0;i<num_independents_old;i++) G[i] = totalgradient[i];
     414
     415                        /*free resources :*/
     416                        xDelete(weightVectorTimesJac);
     417                        xDelete(aWeightVector);
     418                }
     419                /*}}}*/
     420                #elif defined(_HAVE_CODIPACK_)
     421                /*Get gradient for CoDiPack{{{*/
     422                if(VerboseAutodiff())_printf0_("   CoDiPack fos_reverse\n");
     423
     424                /* call the fos_reverse in a loop on the index, from 0 to num_dependents, so
     425                 * as to generate num_dependents gradients: */
     426                for(int dep_index=0;dep_index<num_dependents_old;dep_index++){
     427
     428                        /*initialize direction index in the weights vector: */
     429                        if(my_rank==0){
     430                                if(dep_index<0 || dep_index>=num_dependents || codi_global.output_indices.size() <= dep_index){
     431                                        _error_("index value for dependent index should be in [0,num_dependents-1]");
     432                                }
     433                                tape_codi.setGradient(codi_global.output_indices[dep_index],1.0);
     434                        }
     435                        //feclearexcept(FE_ALL_EXCEPT);
     436                        //feenableexcept(FE_DIVBYZERO | FE_INVALID | FE_OVERFLOW);
     437                        tape_codi.evaluate();
     438
     439                        /*Get gradient for this dependent */
     440                        weightVectorTimesJac=xNew<IssmPDouble>(num_independents);
     441                        auto in_size = codi_global.input_indices.size();
     442                        for(size_t i = 0; i < in_size; ++i){
     443                                _assert_(i<num_independents);
     444                                weightVectorTimesJac[i] = tape_codi.getGradient(codi_global.input_indices[i]);
     445                        }
     446                        if(my_rank==0) for(int i=0;i<num_independents;i++){
     447                                totalgradient[i]+=weightVectorTimesJac[i];
     448                        }
     449                        xDelete(weightVectorTimesJac);
     450                }
     451
     452                /*Clear tape*/
     453                tape_codi.reset();
     454                /*}}}*/
     455                #else
     456                _error_("not suppoted");
     457                #endif
     458
     459                /*Broadcast gradient to other ranks*/
     460                ISSM_MPI_Bcast(totalgradient,num_independents_old,ISSM_MPI_PDOUBLE,0,IssmComm::GetComm());
     461
     462                /*Check size of Jlist to avoid crashes*/
     463                _assert_((*Jlisti)<JlistM);
     464                _assert_(JlistN==num_responses+1);
     465
     466                /*Compute objective function and broadcast it to other cpus*/
     467                *pf = reCast<double>(J);
     468                ISSM_MPI_Bcast(pf,1,ISSM_MPI_PDOUBLE,0,IssmComm::GetComm());
     469
     470                /*Record cost function values and delete Jtemp*/
     471                for(int i=0;i<num_responses;i++) Jlist[(*Jlisti)*JlistN+i] = dependents[i];
     472                Jlist[(*Jlisti)*JlistN+num_responses] = reCast<IssmPDouble>(J);
     473
     474                if(*indic==0){
     475                        /*dry run, no gradient required*/
     476
     477                        /*Retrieve objective functions independently*/
     478                        _printf0_("f(x) = "<<setw(12)<<setprecision(7)<<*pf<<"  |  ");
     479                        _printf0_("            N/A |\n");
     480                        for(int i=0;i<num_responses;i++) _printf0_(" "<<setw(12)<<setprecision(7)<<Jlist[(*Jlisti)*JlistN+i]);
     481                        _printf0_("\n");
     482
     483                        *Jlisti = (*Jlisti) +1;
     484                        xDelete<double>(XU);
     485                        xDelete<double>(XL);
     486                        return;
     487                }
     488
     489                /*Compute gradient*/
     490                for(long i=0;i<num_independents_old;i++) G[i] = totalgradient[i];
     491
     492                xDelete<IssmPDouble>(dependents);
     493                xDelete<IssmPDouble>(totalgradient);
     494          } /*====????*/
    479495
    480496        /*Constrain Gradient*/
    … …  
    496512
    497513        /*Print info*/
     514        _printf0_("f(x) = "<<setw(12)<<setprecision(7)<<*pf<<"  |  ");
    498515        _printf0_("       "<<setw(12)<<setprecision(7)<<Gnorm<<" |");
    499516        for(int i=0;i<num_responses;i++) _printf0_(" "<<setw(12)<<setprecision(7)<<Jlist[(*Jlisti)*JlistN+i]);
    … …  
    509526        xDelete<int>(N);
    510527        xDelete<double>(scaling_factors);
    511         xDelete<IssmPDouble>(dependents);
    512         xDelete<IssmPDouble>(totalgradient);
    513528}/*}}}*/
    514529void controladm1qn3_core(FemModel* femmodel){/*{{{*/
    515530
    516         int checkpoint_frequency;
    517         femmodel->parameters->FindParam(&checkpoint_frequency,SettingsCheckpointFrequencyEnum);
    518         if(checkpoint_frequency){
    519                 #ifdef _HAVE_CODIPACK_
    520                 transient_ad(femmodel);
    521                 femmodel->OutputControlsx(&femmodel->results);
    522                 printf("No optimization implemented yet, skipping\n");
    523                 return;
    524                 #else
    525                 _error_("checkpointing not implemented for ADOLC");
    526                 #endif
    527         }
    528531
    529532        /*Intermediaries*/

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

    r26747 r26889  
    1616#include "../modules/modules.h"
    1717#include "../solutionsequences/solutionsequences.h"
     18
     19#ifdef _HAVE_CODIPACK_
     20extern CoDi_global codi_global;
     21#endif
    1822
    1923/*Prototypes*/
    … …  
    275279
    276280#ifdef _HAVE_CODIPACK_
    277 void transient_ad(FemModel* femmodel){/*{{{*/
     281double transient_ad(FemModel* femmodel, double* G, double* Jlist){/*{{{*/
    278282
    279283        /*parameters: */
    … …  
    282286        bool       isoceancoupling;
    283287        int        step,timestepping;
    284         int        checkpoint_frequency;
     288        int        checkpoint_frequency,num_responses;
    285289
    286290        /*Get rank*/
    … …  
    293297        femmodel->parameters->FindParam(&yts,ConstantsYtsEnum);
    294298        femmodel->parameters->FindParam(&timestepping,TimesteppingTypeEnum);
     299        femmodel->parameters->FindParam(&num_responses,InversionNumCostFunctionsEnum);
    295300        femmodel->parameters->FindParam(&checkpoint_frequency,SettingsCheckpointFrequencyEnum); _assert_(checkpoint_frequency>0);
    296301
    … …  
    343348
    344349                /*Go through our dependent variables, and compute the response:*/
    345                 if(finaltime==time){
    346                         DataSet* dependent_objects=((DataSetParam*)femmodel->parameters->FindParamObject(AutodiffDependentObjectsEnum))->value;
    347                         for(Object* & object:dependent_objects->objects){
    348                                 DependentObject* dep=(DependentObject*)object;
    349                                 dep->Responsex(&output_value,femmodel);
    350                                 dep->AddValue(output_value);
    351                         }
     350                DataSet* dependent_objects=((DataSetParam*)femmodel->parameters->FindParamObject(AutodiffDependentObjectsEnum))->value;
     351                for(Object* & object:dependent_objects->objects){
     352                        DependentObject* dep=(DependentObject*)object;
     353                        dep->Responsex(&output_value,femmodel);
     354                        dep->AddValue(output_value);
    352355                }
    353356
    … …  
    395398        SetControlInputsFromVectorx(femmodel,X);
    396399
    397         IssmDouble J = 0.;
     400        IssmDouble J     = 0.;
     401        int        count = 0;
    398402        DataSet* dependent_objects=((DataSetParam*)femmodel->parameters->FindParamObject(AutodiffDependentObjectsEnum))->value;
    399403        for(Object* & object:dependent_objects->objects){
    400404                DependentObject* dep=(DependentObject*)object;
    401                 J += dep->GetValue();
    402         }
    403         tape_codi.registerOutput(J);
     405                IssmDouble       output_value = dep->GetValue();
     406
     407
     408                J += output_value;
     409
     410                tape_codi.registerOutput(J);
     411                #if _CODIPACK_MAJOR_==2
     412                codi_global.output_indices.push_back(J.getIdentifier());
     413                #elif _CODIPACK_MAJOR_==1
     414                codi_global.output_indices.push_back(J.getGradientData());
     415                #else
     416                #error "_CODIPACK_MAJOR_ not supported"
     417                #endif
     418
     419                /*Keep track of output for printing*/
     420                Jlist[count] = output_value.getValue();
     421      count++;
     422        }
     423        Jlist[count] = J.getValue();
     424        _assert_(count == num_responses);
    404425
    405426        tape_codi.setPassive();
    406         if(IssmComm::GetRank()==0) J.gradient() = 1.0;
    407         tape_codi.evaluate();
     427
     428        if(VerboseAutodiff())_printf0_("   CoDiPack fos_reverse\n");
     429        for(int i=0;i<num_responses;i++){
     430                if(my_rank==0) tape_codi.setGradient(codi_global.output_indices[i],1.0);
     431                tape_codi.evaluate();
     432        }
    408433
    409434        /*Initialize Xb and Yb*/
    … …  
    449474
    450475                        /*Go through our dependent variables, and compute the response:*/
    451                         if(thistime==finaltime){
    452                                 DataSet* dependent_objects=((DataSetParam*)femmodel->parameters->FindParamObject(AutodiffDependentObjectsEnum))->value;
    453                                 for(Object* & object:dependent_objects->objects){
    454                                         DependentObject* dep=(DependentObject*)object;
    455                                         dep->Responsex(&output_value,femmodel);
    456                                         dep->AddValue(output_value);
    457                                 }
     476                        DataSet* dependent_objects=((DataSetParam*)femmodel->parameters->FindParamObject(AutodiffDependentObjectsEnum))->value;
     477                        for(Object* & object:dependent_objects->objects){
     478                                DependentObject* dep=(DependentObject*)object;
     479                                dep->Responsex(&output_value,femmodel);
     480                                dep->AddValue(output_value);
    458481                        }
    459482
    … …  
    484507                for(int i=0; i<Xsize;   i++) Xb[i] += X[i].gradient();
    485508        }
     509
    486510        /*Clear tape*/
    487511        tape_codi.reset();
    488512
    489513        /*Broadcast gradient to other ranks (make sure to sum all gradients)*/
    490         double* gradient=xNew<double>(Xsize);
    491         ISSM_MPI_Allreduce(Xb,gradient,Xsize,ISSM_MPI_PDOUBLE,ISSM_MPI_SUM,IssmComm::GetComm());
    492 
    493         /*Now set final gradient*/
    494         IssmDouble* aG=xNew<IssmDouble>(Xsize);
    495         for(int i=0;i<Xsize;i++) aG[i] = reCast<IssmDouble>(gradient[i]);
    496         xDelete<double>(gradient);
    497         ControlInputSetGradientx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,aG);
    498         xDelete<IssmDouble>(aG);
    499 
     514        ISSM_MPI_Allreduce(Xb,G,Xsize,ISSM_MPI_PDOUBLE,ISSM_MPI_SUM,IssmComm::GetComm());
     515
     516        /*Cleanup and return misfit*/
    500517        xDelete<IssmDouble>(X);
    501518        xDelete<double>(Xb);
    502519        xDelete<double>(Yb);
    503520        xDelete<int>(Yin);
     521   return J.getValue();
    504522}/*}}}*/
    505523#endif

  • issm/trunk-jpl/src/c/modules/ModelProcessorx/CreateOutputDefinitions.cpp

    r26468 r26889  
    311311                                /*cfsurfacelogvel variables: */
    312312                                int          num_cfsurfacelogvels;
    313                                 char**       cfsurfacelogvel_name                                               = NULL;   
    314                                 char**           cfsurfacelogvel_definitionstring               = NULL;   
    315                                 IssmDouble** cfsurfacelogvel_vxobs                      = NULL;
    316                                 IssmDouble** cfsurfacelogvel_vyobs                      = NULL;
    317                                 char**           cfsurfacelogvel_vxobs_string   = NULL;
    318                                 char**           cfsurfacelogvel_vyobs_string   = NULL;
    319                                 int*         cfsurfacelogvel_observation_M                      = NULL;
    320                                 int*         cfsurfacelogvel_observation_N                      = NULL;
    321                                 IssmDouble** cfsurfacelogvel_weights                                    = NULL;
    322                                 int*         cfsurfacelogvel_weights_M                          = NULL;
    323                                 int*         cfsurfacelogvel_weights_N                          = NULL;
    324                                 char**       cfsurfacelogvel_weightstring               = NULL;
    325                                 IssmDouble*  cfsurfacelogvel_datatime                           = NULL;
    326 
    327                                 /*Fetch name, modeltring, observation, observationtring, etc ... (see src/m/classes/cfsurfacelogvel.m): */
    328                                 iomodel->FetchMultipleData(&cfsurfacelogvel_name,&num_cfsurfacelogvels,                                                        "md.cfsurfacelogvel.name");
    329                                 iomodel->FetchMultipleData(&cfsurfacelogvel_definitionstring,&num_cfsurfacelogvels,                                            "md.cfsurfacelogvel.definitionstring");
    330                                 iomodel->FetchMultipleData(&cfsurfacelogvel_vxobs,&cfsurfacelogvel_observation_M,&cfsurfacelogvel_observation_N,&num_cfsurfacelogvels, "md.cfsurfacelogvel.vxobs");
    331                                 iomodel->FetchMultipleData(&cfsurfacelogvel_vxobs_string,&num_cfsurfacelogvels,                                          "md.cfsurfacelogvel.vxobs_string");
    332                                 iomodel->FetchMultipleData(&cfsurfacelogvel_vyobs,&cfsurfacelogvel_observation_M,&cfsurfacelogvel_observation_N,&num_cfsurfacelogvels, "md.cfsurfacelogvel.vyobs");
    333                                 iomodel->FetchMultipleData(&cfsurfacelogvel_vyobs_string,&num_cfsurfacelogvels,                                          "md.cfsurfacelogvel.vyobs_string");                    iomodel->FetchMultipleData(&cfsurfacelogvel_weights,&cfsurfacelogvel_weights_M,&cfsurfacelogvel_weights_N,&num_cfsurfacelogvels,             "md.cfsurfacelogvel.weights");
    334                                 iomodel->FetchMultipleData(&cfsurfacelogvel_weightstring,&num_cfsurfacelogvels,                                              "md.cfsurfacelogvel.weights_string");
    335                                 iomodel->FetchMultipleData(&cfsurfacelogvel_datatime,&num_cfsurfacelogvels,                                                                                                                                      "md.cfsurfacelogvel.datatime");
     313                                char       **cfsurfacelogvel_name             = NULL;
     314                                char       **cfsurfacelogvel_definitionstring = NULL;
     315                                IssmDouble **cfsurfacelogvel_vxobs            = NULL;
     316                                IssmDouble **cfsurfacelogvel_vyobs            = NULL;
     317                                char       **cfsurfacelogvel_vxobs_string     = NULL;
     318                                char       **cfsurfacelogvel_vyobs_string     = NULL;
     319                                int         *cfsurfacelogvel_observation_M    = NULL;
     320                                int         *cfsurfacelogvel_observation_N    = NULL;
     321                                IssmDouble **cfsurfacelogvel_weights          = NULL;
     322                                int         *cfsurfacelogvel_weights_M        = NULL;
     323                                int         *cfsurfacelogvel_weights_N        = NULL;
     324                                char       **cfsurfacelogvel_weightstring     = NULL;
     325                                IssmDouble  *cfsurfacelogvel_datatime         = NULL;
     326
     327            /*Fetch name, modeltring, observation, observationtring, etc ... (see src/m/classes/cfsurfacelogvel.m): */
     328            iomodel->FetchMultipleData(&cfsurfacelogvel_name,&num_cfsurfacelogvels,"md.cfsurfacelogvel.name");
     329            iomodel->FetchMultipleData(&cfsurfacelogvel_definitionstring,&num_cfsurfacelogvels,"md.cfsurfacelogvel.definitionstring");
     330            iomodel->FetchMultipleData(&cfsurfacelogvel_vxobs,&cfsurfacelogvel_observation_M,&cfsurfacelogvel_observation_N,&num_cfsurfacelogvels,"md.cfsurfacelogvel.vxobs");
     331            iomodel->FetchMultipleData(&cfsurfacelogvel_vxobs_string,&num_cfsurfacelogvels,"md.cfsurfacelogvel.vxobs_string");
     332            iomodel->FetchMultipleData(&cfsurfacelogvel_vyobs,NULL,NULL,&num_cfsurfacelogvels,"md.cfsurfacelogvel.vyobs");
     333            iomodel->FetchMultipleData(&cfsurfacelogvel_vyobs_string,&num_cfsurfacelogvels,"md.cfsurfacelogvel.vyobs_string");
     334            iomodel->FetchMultipleData(&cfsurfacelogvel_weights,&cfsurfacelogvel_weights_M,&cfsurfacelogvel_weights_N,&num_cfsurfacelogvels,"md.cfsurfacelogvel.weights");
     335            iomodel->FetchMultipleData(&cfsurfacelogvel_weightstring,&num_cfsurfacelogvels,"md.cfsurfacelogvel.weights_string");
     336            iomodel->FetchMultipleData(&cfsurfacelogvel_datatime,&num_cfsurfacelogvels,"md.cfsurfacelogvel.datatime");
    336337
    337338                                for(j=0;j<num_cfsurfacelogvels;j++){
Note: See TracChangeset for help on using the changeset viewer.