Changeset 23267


Ignore:
Timestamp:
09/12/18 11:25:22 (7 years ago)
Author:
Mathieu Morlighem
Message:

NEW: integrating CoDiPack in m1qn3, still crashes after first it

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

Legend:

Unmodified
Added
Removed

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

    r23250 r23267  
    1111#include "../solutionsequences/solutionsequences.h"
    1212
    13 #if defined (_HAVE_M1QN3_) && defined(_HAVE_ADOLC_)
     13#ifdef _HAVE_CODIPACK_
     14extern CoDi_global codi_global;
     15#include <sstream> // for output of the CoDiPack tape
     16#endif
     17
     18#if defined (_HAVE_M1QN3_) && defined(_HAVE_AD_)
    1419/*m1qn3 prototypes*/
    1520extern "C" void *ctonbe_; // DIS mode : Conversion
    … …  
    3540void simul_starttrace(FemModel* femmodel){/*{{{*/
    3641
     42        #if defined(_HAVE_ADOLC_)
    3743        /*Retrive ADOLC parameters*/
    3844        IssmDouble gcTriggerRatio;
    … …  
    5763        int my_rank=IssmComm::GetRank();
    5864        trace_on(my_rank,keepTaylors,reCast<size_t>(obufsize),reCast<size_t>(lbufsize),reCast<size_t>(cbufsize),reCast<size_t>(tbufsize),skipFileDeletion);
     65
     66        #elif defined(_HAVE_CODIPACK_)
     67
     68                //fprintf(stderr, "*** Codipack IoModel::StartTrace\n");
     69                /*
     70                 * FIXME codi
     71                 * - ADOL-C variant uses fine grained tracing with various arguments
     72                 * - ADOL-C variant sets a garbage collection parameter for its tape
     73                 * -> These parameters are not read for the CoDiPack ISSM version!
     74                 */
     75                auto& tape_codi = IssmDouble::getGlobalTape();
     76                tape_codi.setActive();
     77                #if _AD_TAPE_ALLOC_
     78                //alloc_profiler.Tag(StartInit, true);
     79                IssmDouble x_t(1.0), y_t(1.0);
     80                tape_codi.registerInput(y_t);
     81                int codi_allocn = 0;
     82                femmodel->parameters->FindParam(&codi_allocn,AutodiffTapeAllocEnum);
     83                for(int i = 0;i < codi_allocn;++i) {
     84                        x_t = y_t * y_t;
     85                }
     86                /*
     87                std::stringstream out_s;
     88                IssmDouble::getGlobalTape().printStatistics(out_s);
     89                _printf0_("StartTrace::Tape Statistics     : TapeAlloc count=[" << codi_allocn << "]\n" << out_s.str());
     90                */
     91                tape_codi.reset();
     92                //alloc_profiler.Tag(FinishInit, true);
     93                #endif
     94
     95        #else
     96        _error_("not implemented");
     97        #endif
    5998}/*}}}*/
    60 void simul_ad(long* indic,long* n,double* X,double* pf,double* G,long izs[1],float rzs[1],void* dzs){/*{{{*/
    61 
    62         /*Get rank*/
    63         int my_rank=IssmComm::GetRank();
    64 
    65         /*Recover Arguments*/
    66         m1qn3_struct *input_struct = (m1qn3_struct*)dzs;
    67 
    68         FemModel* femmodel = input_struct->femmodel;
    69         int num_responses,num_controls,numberofvertices,solution_type;
    70         femmodel->parameters->FindParam(&solution_type,SolutionTypeEnum);
    71         int* N = NULL;
    72         int N_add = 0;
    73         int* control_enum = NULL;
    74 
    75         if (solution_type == TransientSolutionEnum){
    76                 femmodel = input_struct->femmodel->copy();
    77         }
    78 
    79         IssmPDouble*  Jlist        = input_struct->Jlist;
    80         int           JlistM       = input_struct->M;
    81         int           JlistN       = input_struct->N;
    82         int*          Jlisti       = input_struct->i;
    83         int           intn         = (int)*n;
    84 
    85         /*Recover some parameters*/
    86         IssmDouble* scaling_factors = NULL;
    87         femmodel->parameters->FindParam(&num_responses,InversionNumCostFunctionsEnum);
    88         femmodel->parameters->FindParam(&num_controls,InversionNumControlParametersEnum);
    89         femmodel->parameters->FindParam(&scaling_factors,NULL,InversionControlScalingFactorsEnum);
    90         femmodel->parameters->FindParam(&N,NULL,ControlInputSizeNEnum);
    91         femmodel->parameters->FindParam(&control_enum,NULL,InversionControlParametersEnum);
    92         numberofvertices=femmodel->vertices->NumberOfVertices();
    93 
    94         /*Constrain input vector and update controls*/
    95         double  *XL = NULL;
    96         double  *XU = NULL;
    97         GetPassiveVectorFromControlInputsx(&XL,femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,"lowerbound");
    98         GetPassiveVectorFromControlInputsx(&XU,femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,"upperbound");
    99 
    100         N_add = 0;
    101         for (int c=0;c<num_controls;c++){
    102                 for(int i=0;i<numberofvertices*N[c];i++){
    103                         int index = N_add*numberofvertices+i;
    104                         X[index] = X[index]*reCast<double>(scaling_factors[c]);
    105                         if(X[index]>XU[index]) X[index]=XU[index];
    106                         if(X[index]<XL[index]) X[index]=XL[index];
    107                 }
    108                 N_add+=N[c];
    109         }
    110 
    111         /*Start Tracing*/
    112         simul_starttrace(femmodel);
    113         /*Set X as our new control input and as INDEPENDENT*/
    114 #ifdef _HAVE_AD_
    115         IssmDouble* aX=xNew<IssmDouble>(intn,"t");
    116 #else
    117         IssmDouble* aX=xNew<IssmDouble>(intn);
    118 #endif
    119         if(my_rank==0){
    120                 for(int i=0;i<intn;i++){
    121                         aX[i]<<=X[i];
    122                 }
    123         }
    124 
    125         ISSM_MPI_Bcast(aX,intn,ISSM_MPI_DOUBLE,0,IssmComm::GetComm());
    126         SetControlInputsFromVectorx(femmodel,aX);
    127         xDelete<IssmDouble>(aX);
    128 
    129         /*Compute solution (forward)*/
    130         void (*solutioncore)(FemModel*)=NULL;
    131         CorePointerFromSolutionEnum(&solutioncore,femmodel->parameters,solution_type);
    132         solutioncore(femmodel);
    133 
    134         /*Reset the time to zero for next optimization*/
    135         if(solution_type==TransientSolutionEnum){
    136                 IssmDouble restart_time;
    137                 femmodel->parameters->FindParam(&restart_time,TimesteppingStartTimeEnum);
    138                 femmodel->parameters->SetParam(restart_time,TimeEnum);
    139 
    140         }
    141 
    142         /*Get Dependents*/
    143         IssmDouble  output_value;
    144         int         num_dependents;
    145         IssmPDouble *dependents;
    146         DataSet*    dependent_objects=NULL;
    147         IssmDouble      J=0.;
    148         femmodel->parameters->FindParam(&num_dependents,AutodiffNumDependentsEnum);
    149         femmodel->parameters->FindParam(&dependent_objects,AutodiffDependentObjectsEnum);
    150 
    151         /*Go through our dependent variables, and compute the response:*/
    152         dependents=xNew<IssmPDouble>(num_dependents);
    153         for(int i=0;i<dependent_objects->Size();i++){
    154                 DependentObject* dep=(DependentObject*)dependent_objects->GetObjectByOffset(i);
    155                 if(solution_type==TransientSolutionEnum) output_value = dep->GetValue();
    156                 if(solution_type!=TransientSolutionEnum) dep->Responsex(&output_value,femmodel);
    157                 if (my_rank==0) {
    158                         output_value>>=dependents[i];
    159                         J+=output_value;
    160                 }
    161         }
    162 
    163         /*Turning off trace tape*/
     99void simul_stoptrace(){/*{{{*/
     100
     101        #if defined(_HAVE_ADOLC_)
    164102        trace_off();
    165         //time_t now = time(NULL);
    166         //if(my_rank==0) _printf_("\nTIME: "<<now<<"\n");
    167 
    168         /*Print tape statistics so that user can kill this run if something is off already:*/
    169103        if(VerboseAutodiff()){ /*{{{*/
    170104
    171105                #ifdef _HAVE_ADOLC_
     106                int my_rank=IssmComm::GetRank();
    172107                size_t  tape_stats[15];
    173108                tapestats(my_rank,tape_stats); //reading of tape statistics
    … …  
    222157        } /*}}}*/
    223158
     159        #elif defined(_HAVE_CODIPACK_)
     160        auto& tape_codi = IssmDouble::getGlobalTape();
     161        tape_codi.setPassive();
     162        if(VerboseAutodiff()){
     163                int my_rank=IssmComm::GetRank();
     164                if(my_rank == 0) {
     165                        // FIXME codi "just because" for now
     166                        tape_codi.printStatistics(std::cout);
     167                        codi_global.print(std::cout);
     168                }
     169        }
     170        #else
     171        _error_("not implemented");
     172        #endif
     173}/*}}}*/
     174void simul_ad(long* indic,long* n,double* X,double* pf,double* G,long izs[1],float rzs[1],void* dzs){/*{{{*/
     175
     176        /*Get rank*/
     177        int my_rank=IssmComm::GetRank();
     178
     179        /*Recover Arguments*/
     180        m1qn3_struct *input_struct = (m1qn3_struct*)dzs;
     181
     182        FemModel* femmodel = input_struct->femmodel;
     183        int num_responses,num_controls,numberofvertices,solution_type;
     184        femmodel->parameters->FindParam(&solution_type,SolutionTypeEnum);
     185        int* N = NULL;
     186        int N_add = 0;
     187        int* control_enum = NULL;
     188
     189        /*In transient, we need to make sure we do not modify femmodel at each iteration, make a copy*/
     190        if(solution_type == TransientSolutionEnum) femmodel = input_struct->femmodel->copy();
     191
     192        IssmPDouble*  Jlist        = input_struct->Jlist;
     193        int           JlistM       = input_struct->M;
     194        int           JlistN       = input_struct->N;
     195        int*          Jlisti       = input_struct->i;
     196        int           intn         = (int)*n;
     197
     198        /*Recover some parameters*/
     199        IssmDouble* scaling_factors = NULL;
     200        femmodel->parameters->FindParam(&num_responses,InversionNumCostFunctionsEnum);
     201        femmodel->parameters->FindParam(&num_controls,InversionNumControlParametersEnum);
     202        femmodel->parameters->FindParam(&scaling_factors,NULL,InversionControlScalingFactorsEnum);
     203        femmodel->parameters->FindParam(&N,NULL,ControlInputSizeNEnum);
     204        femmodel->parameters->FindParam(&control_enum,NULL,InversionControlParametersEnum);
     205        numberofvertices=femmodel->vertices->NumberOfVertices();
     206
     207        /*Constrain input vector and update controls*/
     208        double  *XL = NULL;
     209        double  *XU = NULL;
     210        GetPassiveVectorFromControlInputsx(&XL,femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,"lowerbound");
     211        GetPassiveVectorFromControlInputsx(&XU,femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,"upperbound");
     212
     213        N_add = 0;
     214        for (int c=0;c<num_controls;c++){
     215                for(int i=0;i<numberofvertices*N[c];i++){
     216                        int index = N_add*numberofvertices+i;
     217                        X[index] = X[index]*reCast<double>(scaling_factors[c]);
     218                        if(X[index]>XU[index]) X[index]=XU[index];
     219                        if(X[index]<XL[index]) X[index]=XL[index];
     220                }
     221                N_add+=N[c];
     222        }
     223
     224        /*Start Tracing*/
     225        simul_starttrace(femmodel);
     226        /*Set X as our new control input and as INDEPENDENT*/
     227#ifdef _HAVE_AD_
     228        IssmDouble* aX=xNew<IssmDouble>(intn,"t");
     229#else
     230        IssmDouble* aX=xNew<IssmDouble>(intn);
     231#endif
     232
     233        #if defined(_HAVE_ADOLC_)
     234        if(my_rank==0){
     235                for(int i=0;i<intn;i++){
     236                        aX[i]<<=X[i];
     237                }
     238        }
     239        #elif defined(_HAVE_CODIPACK_)
     240        auto& tape_codi = IssmDouble::getGlobalTape();
     241        if(my_rank==0){
     242                for (int i=0;i<intn;i++) {
     243                        aX[i]=X[i];
     244                        tape_codi.registerInput(aX[i]);
     245                        codi_global.input_indices.push_back(aX[i].getGradientData());
     246                }
     247        }
     248        #else
     249        _error_("not suppoted");
     250        #endif
     251
     252        ISSM_MPI_Bcast(aX,intn,ISSM_MPI_DOUBLE,0,IssmComm::GetComm());
     253        SetControlInputsFromVectorx(femmodel,aX);
     254        xDelete<IssmDouble>(aX);
     255
     256        /*Compute solution (forward)*/
     257        void (*solutioncore)(FemModel*)=NULL;
     258        CorePointerFromSolutionEnum(&solutioncore,femmodel->parameters,solution_type);
     259        solutioncore(femmodel);
     260
     261        /*Reset the time to zero for next optimization*/
     262        if(solution_type==TransientSolutionEnum){
     263                IssmDouble restart_time;
     264                femmodel->parameters->FindParam(&restart_time,TimesteppingStartTimeEnum);
     265                femmodel->parameters->SetParam(restart_time,TimeEnum);
     266
     267        }
     268
     269        /*Get Dependents*/
     270        IssmDouble  output_value;
     271        int         num_dependents;
     272        IssmPDouble *dependents;
     273        DataSet*    dependent_objects=NULL;
     274        IssmDouble      J=0.;
     275        femmodel->parameters->FindParam(&num_dependents,AutodiffNumDependentsEnum);
     276        femmodel->parameters->FindParam(&dependent_objects,AutodiffDependentObjectsEnum);
     277
     278        /*Go through our dependent variables, and compute the response:*/
     279        dependents=xNew<IssmPDouble>(num_dependents);
     280        for(int i=0;i<dependent_objects->Size();i++){
     281                DependentObject* dep=(DependentObject*)dependent_objects->GetObjectByOffset(i);
     282                if(solution_type==TransientSolutionEnum) output_value = dep->GetValue();
     283                if(solution_type!=TransientSolutionEnum) dep->Responsex(&output_value,femmodel);
     284                if(my_rank==0) {
     285
     286                        #if defined(_HAVE_CODIPACK_)
     287                        tape_codi.registerOutput(output_value);
     288                        dependents[i] = output_value.getValue();
     289                        codi_global.output_indices.push_back(output_value.getGradientData());
     290
     291                        #elif defined(_HAVE_ADOLC_)
     292                        output_value>>=dependents[i];
     293
     294                        #else
     295                        _error_("not suppoted");
     296                        #endif
     297                        J+=output_value;
     298                }
     299        }
     300
     301        /*Turning off trace tape*/
     302        simul_stoptrace();
     303        //time_t now = time(NULL);
     304        //if(my_rank==0) _printf_("\nTIME: "<<now<<"\n");
     305
    224306        /*diverse: */
    225307        int  dummy;
    … …  
    246328                num_independents = 0;
    247329        }
     330
     331        #if defined(_HAVE_ADOLC_)
     332        /*Get gradient for ADOLC {{{*/
    248333
    249334        /*get the EDF pointer:*/
    … …  
    290375                xDelete(aWeightVector);
    291376        }
     377        /*}}}*/
     378        #elif defined(_HAVE_CODIPACK_)
     379        /*Get gradient for CoDiPack{{{*/
     380        if(VerboseAutodiff())_printf0_("   CoDiPack fos_reverse\n");
     381        int     aDepIndex=0; /*FIXME: do we really need this?*/
     382
     383        /*retrieve direction index: */
     384        femmodel->parameters->FindParam(&aDepIndex,AutodiffFosReverseIndexEnum);
     385        if (my_rank==0) {
     386                if (aDepIndex<0 || aDepIndex>=num_dependents || codi_global.output_indices.size() <= aDepIndex){
     387                        _error_("index value for AutodiffFosReverseIndexEnum should be in [0,num_dependents-1]");
     388                }
     389                tape_codi.setGradient(codi_global.output_indices[aDepIndex], 1.0);
     390        }
     391        tape_codi.evaluate();
     392
     393        weightVectorTimesJac=xNew<double>(num_independents);
     394        /*call driver: */
     395        auto in_size = codi_global.input_indices.size();
     396        for(size_t i = 0; i < in_size; ++i) {
     397                weightVectorTimesJac[i] = tape_codi.getGradient(codi_global.input_indices[i]);
     398        }
     399
     400        /*Add to totalgradient: */
     401        totalgradient=xNewZeroInit<IssmPDouble>(num_independents_old);
     402        if(my_rank==0) for(int i=0;i<num_independents;i++) {
     403                totalgradient[i]+=weightVectorTimesJac[i];
     404        }
     405
     406        /*free resources :*/
     407        xDelete(weightVectorTimesJac);
     408        /*}}}*/
     409        #else
     410        _error_("not suppoted");
     411        #endif
    292412
    293413        /*Broadcast gradient to other ranks*/
    … …  
    551671
    552672#else
    553 void controladm1qn3_core(FemModel* femmodel){_error_("M1QN3 or ADOLC not installed");}
     673void controladm1qn3_core(FemModel* femmodel){_error_("M1QN3 or ADOLC/CoDiPack not installed");}
    554674#endif //_HAVE_M1QN3_

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

    r23252 r23267  
    103103
    104104        if(isautodiff){
    105 #if _HAVE_ADOLC_
    106                 /*Copy some parameters from IoModel to parameters dataset: {{{*/
     105                #if defined(_HAVE_ADOLC_)
     106                /*Copy some parameters from IoModel to parameters dataset*/
    107107                parameters->AddObject(iomodel->CopyConstantObject("md.autodiff.obufsize",AutodiffObufsizeEnum));
    108108                parameters->AddObject(iomodel->CopyConstantObject("md.autodiff.cbufsize",AutodiffCbufsizeEnum));
    … …  
    111111                parameters->AddObject(iomodel->CopyConstantObject("md.autodiff.gcTriggerRatio",AutodiffGcTriggerRatioEnum));
    112112                parameters->AddObject(iomodel->CopyConstantObject("md.autodiff.gcTriggerMaxSize",AutodiffGcTriggerMaxSizeEnum));
    113                 /*}}}*/
    114 #endif
     113
     114                #elif defined(_HAVE_CODIPACK_)
     115                parameters->AddObject(iomodel->CopyConstantObject("md.autodiff.tapeAlloc",AutodiffTapeAllocEnum));
     116
     117                #else
     118                _error_("not supported yet");
     119                #endif
     120
    115121                /*retrieve driver: {{{*/
    116122                iomodel->FindConstant(&autodiff_driver,"md.autodiff.driver");
Note: See TracChangeset for help on using the changeset viewer.