Changeset 23318


Ignore:
Timestamp:
09/19/18 12:43:18 (7 years ago)
Author:
Mathieu Morlighem
Message:

CHG: make validation core work for AD (youhou)

File:
1 edited

Legend:

Unmodified
Added
Removed

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

    r23242 r23318  
    99#include "../modules/modules.h"
    1010
     11#ifdef _HAVE_CODIPACK_
     12extern CoDi_global codi_global;
     13#include <sstream> // for output of the CoDiPack tape
     14#endif
     15
     16#ifdef _HAVE_AD_
     17void simul_starttrace2(FemModel* femmodel){/*{{{*/
     18
     19        #if defined(_HAVE_ADOLC_)
     20        /*Retrive ADOLC parameters*/
     21        IssmDouble gcTriggerRatio;
     22        IssmDouble gcTriggerMaxSize;
     23        IssmDouble obufsize;
     24        IssmDouble lbufsize;
     25        IssmDouble cbufsize;
     26        IssmDouble tbufsize;
     27        femmodel->parameters->FindParam(&gcTriggerRatio,AutodiffGcTriggerRatioEnum);
     28        femmodel->parameters->FindParam(&gcTriggerMaxSize,AutodiffGcTriggerMaxSizeEnum);
     29        femmodel->parameters->FindParam(&obufsize,AutodiffObufsizeEnum);
     30        femmodel->parameters->FindParam(&lbufsize,AutodiffLbufsizeEnum);
     31        femmodel->parameters->FindParam(&cbufsize,AutodiffCbufsizeEnum);
     32        femmodel->parameters->FindParam(&tbufsize,AutodiffTbufsizeEnum);
     33
     34        /*Set garbage collection parameters: */
     35        setStoreManagerControl(reCast<IssmPDouble>(gcTriggerRatio),reCast<size_t>(gcTriggerMaxSize));
     36
     37        /*Start trace: */
     38        int skipFileDeletion=1;
     39        int keepTaylors=1;
     40        int my_rank=IssmComm::GetRank();
     41        trace_on(my_rank,keepTaylors,reCast<size_t>(obufsize),reCast<size_t>(lbufsize),reCast<size_t>(cbufsize),reCast<size_t>(tbufsize),skipFileDeletion);
     42
     43        #elif defined(_HAVE_CODIPACK_)
     44
     45                //fprintf(stderr, "*** Codipack IoModel::StartTrace\n");
     46                /*
     47                 * FIXME codi
     48                 * - ADOL-C variant uses fine grained tracing with various arguments
     49                 * - ADOL-C variant sets a garbage collection parameter for its tape
     50                 * -> These parameters are not read for the CoDiPack ISSM version!
     51                 */
     52                auto& tape_codi = IssmDouble::getGlobalTape();
     53                tape_codi.setActive();
     54                #if _AD_TAPE_ALLOC_
     55                //alloc_profiler.Tag(StartInit, true);
     56                IssmDouble x_t(1.0), y_t(1.0);
     57                tape_codi.registerInput(y_t);
     58                int codi_allocn = 0;
     59                femmodel->parameters->FindParam(&codi_allocn,AutodiffTapeAllocEnum);
     60                for(int i = 0;i < codi_allocn;++i) {
     61                        x_t = y_t * y_t;
     62                }
     63                /*
     64                std::stringstream out_s;
     65                IssmDouble::getGlobalTape().printStatistics(out_s);
     66                _printf0_("StartTrace::Tape Statistics     : TapeAlloc count=[" << codi_allocn << "]\n" << out_s.str());
     67                */
     68                tape_codi.reset();
     69                //alloc_profiler.Tag(FinishInit, true);
     70                #else
     71                tape_codi.reset();
     72                #endif
     73
     74        #else
     75        _error_("not implemented");
     76        #endif
     77}/*}}}*/
     78void simul_stoptrace2(){/*{{{*/
     79
     80        #if defined(_HAVE_ADOLC_)
     81        trace_off();
     82        if(VerboseAutodiff()){ /*{{{*/
     83
     84                #ifdef _HAVE_ADOLC_
     85                int my_rank=IssmComm::GetRank();
     86                size_t  tape_stats[15];
     87                tapestats(my_rank,tape_stats); //reading of tape statistics
     88                int commSize=IssmComm::GetSize();
     89                int *sstats=new int[7];
     90                sstats[0]=tape_stats[NUM_OPERATIONS];
     91                sstats[1]=tape_stats[OP_FILE_ACCESS];
     92                sstats[2]=tape_stats[NUM_LOCATIONS];
     93                sstats[3]=tape_stats[LOC_FILE_ACCESS];
     94                sstats[4]=tape_stats[NUM_VALUES];
     95                sstats[5]=tape_stats[VAL_FILE_ACCESS];
     96                sstats[6]=tape_stats[TAY_STACK_SIZE];
     97                int *rstats=NULL;
     98                if (my_rank==0) rstats=new int[commSize*7];
     99                ISSM_MPI_Gather(sstats,7,ISSM_MPI_INT,rstats,7,ISSM_MPI_INT,0,IssmComm::GetComm());
     100                if (my_rank==0) {
     101                        int offset=50;
     102                        int rOffset=(commSize/10)+1;
     103                        _printf_("   ADOLC statistics: \n");
     104                        _printf_("     "<<setw(offset)<<left<<"#independents: " <<setw(12)<<right<<tape_stats[NUM_INDEPENDENTS] << "\n");
     105                        _printf_("     "<<setw(offset)<<left<<"#dependents: " <<setw(12)<<right<<tape_stats[NUM_DEPENDENTS] << "\n");
     106                        _printf_("     "<<setw(offset)<<left<<"max #live active variables: " <<setw(12)<<right<<tape_stats[NUM_MAX_LIVES] << "\n");
     107                        _printf_("     operations: entry size "<< sizeof(unsigned char) << " Bytes \n");
     108                        _printf_("     "<<setw(offset)<<left<<"  #entries in buffer (AutodiffObufsizeEnum) " <<setw(12)<<right<<tape_stats[OP_BUFFER_SIZE] << "\n");
     109                        for (int r=0;r<commSize;++r)
     110                         _printf_("       ["<<setw(rOffset)<<right<<r<<"]"<<setw(offset-rOffset-4)<<left<<" #entries total" <<setw(12)<<right<<rstats[r*7+0] << (rstats[r*7+1]?" ->file":"") << "\n");
     111                        _printf_("     locations: entry size " << sizeof(locint) << " Bytes\n");
     112                        _printf_("     "<<setw(offset)<<left<<"  #entries in buffer (AutodiffLbufsizeEnum) " <<setw(12)<<right<<tape_stats[LOC_BUFFER_SIZE] << "\n");
     113                        for (int r=0;r<commSize;++r)
     114                         _printf_("       ["<<setw(rOffset)<<right<<r<<"]"<<setw(offset-rOffset-4)<<left<<" #entries total" <<setw(12)<<right<<rstats[r*7+2] << (rstats[r*7+3]?" ->file":"") << "\n");
     115                        _printf_("     constant values: entry size " << sizeof(double) << " Bytes\n");
     116                        _printf_("     "<<setw(offset)<<left<<"  #entries in buffer (AutodiffCbufsizeEnum) " <<setw(12)<<right<<tape_stats[VAL_BUFFER_SIZE] << "\n");
     117                        for (int r=0;r<commSize;++r)
     118                         _printf_("       ["<<setw(rOffset)<<right<<r<<"]"<<setw(offset-rOffset-4)<<left<<" #entries total" <<setw(12)<<right<<rstats[r*7+4] << (rstats[r*7+5]?" ->file":"") << "\n");
     119                        _printf_("     Taylor stack: entry size " << sizeof(revreal) << " Bytes\n");
     120                        _printf_("     "<<setw(offset)<<left<<"  #entries in buffer (AutodiffTbufsizeEnum) " <<setw(12)<<right<<tape_stats[TAY_BUFFER_SIZE] << "\n");
     121                        for (int r=0;r<commSize;++r)
     122                         _printf_("       ["<<setw(rOffset)<<right<<r<<"]"<<setw(offset-rOffset-4)<<left<<" #entries total" <<setw(12)<<right<<rstats[r*7+6] << (rstats[r*7+6]>tape_stats[TAY_BUFFER_SIZE]?" ->file":"") << "\n");
     123                        delete []rstats;
     124                }
     125                delete [] sstats;
     126                #endif
     127
     128                #ifdef _HAVE_CODIPACK_
     129                #ifdef _AD_TAPE_ALLOC_
     130                //_printf_("Allocation time  P(" << my_rank << "): " << alloc_profiler.DeltaTime(StartInit, FinishInit) << "\n");
     131                #endif
     132                std::stringstream out_s;
     133                IssmDouble::getGlobalTape().printStatistics(out_s);
     134                _printf0_("CoDiPack Profiling::Tape Statistics :\n" << out_s.str());
     135                #endif
     136        } /*}}}*/
     137
     138        #elif defined(_HAVE_CODIPACK_)
     139        auto& tape_codi = IssmDouble::getGlobalTape();
     140        tape_codi.setPassive();
     141        if(VerboseAutodiff()){
     142                int my_rank=IssmComm::GetRank();
     143                if(my_rank == 0) {
     144                        // FIXME codi "just because" for now
     145                        tape_codi.printStatistics(std::cout);
     146                        codi_global.print(std::cout);
     147                }
     148        }
     149        #else
     150        _error_("not implemented");
     151        #endif
     152}/*}}}*/
     153#endif
     154
    11155void controlvalidation_core(FemModel* femmodel){
    12156
    13157        int         solution_type,n;
    14158        int         num_responses;
    15         IssmDouble  j0,j,yts;
     159        IssmDouble  j0,j;
    16160        IssmDouble  Ialpha,exponent,alpha;
    17161        IssmDouble* scaling_factors = NULL;
    18162        IssmDouble* jlist = NULL;
    19         IssmDouble *G = NULL;
    20         IssmDouble *X = NULL;
    21         IssmDouble *X0= NULL;
    22 
    23         /*Solution and Adjoint core pointer*/
    24         void (*solutioncore)(FemModel*) = NULL;
    25         void (*adjointcore)(FemModel*)  = NULL;
     163        int my_rank=IssmComm::GetRank();
    26164
    27165        /*Recover parameters used throughout the solution*/
    … …  
    29167        femmodel->parameters->SetParam(false,SaveResultsEnum);
    30168        femmodel->parameters->FindParam(&num_responses,InversionNumCostFunctionsEnum);
    31         femmodel->parameters->FindParam(&yts,ConstantsYtsEnum);
    32169        femmodel->parameters->FindParam(&scaling_factors,NULL,InversionControlScalingFactorsEnum);
    33170
    34171        /*Get initial guess*/
    35         Vector<IssmDouble> *Xpetsc = NULL;
    36         GetVectorFromControlInputsx(&Xpetsc,femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,"value");
    37         Xpetsc->GetSize(&n);
    38         X0 = Xpetsc->ToMPISerial();
    39         delete Xpetsc;
    40 
    41         /*Allocate current vector*/
    42         X = xNew<IssmDouble>(n);
     172        IssmPDouble* X0 = NULL;
     173        GetPassiveVectorFromControlInputsx(&X0,&n,femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,"value");
     174
     175        /*Allocate vectors*/
     176        IssmDouble*  X = xNew<IssmDouble>(n);
     177        IssmPDouble* G = xNew<IssmPDouble>(n);
    43178
    44179        /*out of solution_type, figure out solution core and adjoint function pointer*/
     180        void (*solutioncore)(FemModel*)=NULL;
    45181        CorePointerFromSolutionEnum(&solutioncore,femmodel->parameters,solution_type);
     182
     183        #if defined(_HAVE_ADOLC_)
     184        /*{{{*/
     185        IssmDouble* aX=xNew<IssmDouble>(n);
     186        if(my_rank==0){
     187                for(int i=0;i<intn;i++){
     188                        aX[i]<<=X0[i];
     189                }
     190        }
     191        _error_("not implemented yet...");
     192        /*}}}*/
     193        #elif defined(_HAVE_CODIPACK_)
     194        simul_starttrace2(femmodel);
     195        IssmDouble* aX=xNew<IssmDouble>(n);
     196        auto& tape_codi = IssmDouble::getGlobalTape();
     197        codi_global.input_indices.clear();
     198        if(my_rank==0){
     199                for (int i=0;i<n;i++) {
     200                        aX[i]=X0[i];
     201                        tape_codi.registerInput(aX[i]);
     202                        codi_global.input_indices.push_back(aX[i].getGradientData());
     203                }
     204        }
     205        SetControlInputsFromVectorx(femmodel,aX);
     206        xDelete(aX);
     207
     208        if(VerboseControl()) _printf0_("   Compute Initial cost function\n");
     209        solutioncore(femmodel);
     210
     211        /*Get Dependents*/
     212        IssmDouble  output_value;
     213        int         num_dependents;
     214        IssmPDouble *dependents;
     215        DataSet*    dependent_objects=NULL;
     216        IssmDouble      J=0.;
     217        femmodel->parameters->FindParam(&num_dependents,AutodiffNumDependentsEnum);
     218        femmodel->parameters->FindParam(&dependent_objects,AutodiffDependentObjectsEnum);
     219
     220        /*Go through our dependent variables, and compute the response:*/
     221        dependents=xNew<IssmPDouble>(num_dependents);
     222        codi_global.output_indices.clear();
     223        for(int i=0;i<dependent_objects->Size();i++){
     224                DependentObject* dep=(DependentObject*)dependent_objects->GetObjectByOffset(i);
     225                if(solution_type==TransientSolutionEnum){
     226                        output_value = dep->GetValue();
     227                }
     228                else{
     229                        dep->Responsex(&output_value,femmodel);
     230                }
     231                _printf0_("=== output ="<<output_value<<" \n");
     232                if(my_rank==0) {
     233                        tape_codi.registerOutput(output_value);
     234                        dependents[i] = output_value.getValue();
     235                        codi_global.output_indices.push_back(output_value.getGradientData());
     236                        J+=output_value;
     237                }
     238        }
     239        j0 = J;
     240        _printf0_("Initial cost function J(x) = "<<setw(12)<<setprecision(7)<<j0<<"\n");
     241        _assert_(j0>0.);
     242        simul_stoptrace2();
     243        /*initialize direction index in the weights vector: */
     244        if(my_rank==0){
     245                tape_codi.setGradient(codi_global.output_indices[0],1.0);
     246        }
     247        tape_codi.evaluate();
     248
     249        /*Get gradient for this dependent */
     250        auto in_size = codi_global.input_indices.size();
     251        for(size_t i = 0; i < in_size; ++i) {
     252                G[i] = tape_codi.getGradient(codi_global.input_indices[i]);
     253        }
     254        #else
     255        /*{{{*/
     256        void (*adjointcore)(FemModel*)  = NULL;
    46257        AdjointCorePointerFromSolutionEnum(&adjointcore,solution_type);
    47258
    … …  
    59270        Gradjx(&G,NULL,femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters);
    60271        for(int i=0;i<n;i++) G[i] = -G[i];
     272        /*}}}*/
     273        #endif
    61274
    62275        /*Allocate output*/
    … …  
    76289                SetControlInputsFromVectorx(femmodel,X);
    77290                solutioncore(femmodel);
     291
     292                #if defined(_HAVE_CODIPACK_)
     293                j=0.;
     294                for(int i=0;i<dependent_objects->Size();i++){
     295                        DependentObject* dep=(DependentObject*)dependent_objects->GetObjectByOffset(i);
     296                        if(solution_type==TransientSolutionEnum){
     297                                output_value = dep->GetValue();
     298                        }
     299                        else{
     300                                dep->Responsex(&output_value,femmodel);
     301                        }
     302                        j+=output_value;
     303                }
     304                #else
    78305                femmodel->CostFunctionx(&j,NULL,NULL);
     306                #endif
    79307
    80308                IssmDouble Den = 0.;
    81309                for(int i=0;i<n;i++) Den += alpha* G[i] * scaling_factors[0];
    82310                Ialpha = fabs((j - j0)/Den - 1.);
     311                _assert_(fabs(Den)>0.);
    83312
    84313                _printf0_(" " << setw(11) << setprecision (5)<<alpha<<" " << setw(11) << setprecision (5)<<Ialpha<<"\n");
    … …  
    93322        femmodel->results->AddObject(new GenericExternalResult<IssmPDouble*>(femmodel->results->Size()+1,JEnum,J_passive,num,2,0,0));
    94323        xDelete<IssmPDouble>(J_passive);
     324        IssmDouble* aG=xNew<IssmDouble>(n);
     325        for(int i=0;i<n;i++) aG[i] = G[i];
     326        ControlInputSetGradientx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,aG);
     327        xDelete<IssmDouble>(aG);
    95328        #else
    96329        femmodel->results->AddObject(new GenericExternalResult<IssmPDouble*>(femmodel->results->Size()+1,JEnum,output,num,2,0,0));
     330        ControlInputSetGradientx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,G);
    97331        #endif
    98         ControlInputSetGradientx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,G);
    99332        femmodel->OutputControlsx(&femmodel->results);
    100333
    101334        /*Clean up and return*/
    102335        xDelete<IssmDouble>(output);
    103         xDelete<IssmDouble>(G);
     336        xDelete<IssmPDouble>(G);
    104337        xDelete<IssmDouble>(X);
    105         xDelete<IssmDouble>(X0);
     338        xDelete<double>(X0);
    106339        xDelete<IssmDouble>(scaling_factors);
    107340}
Note: See TracChangeset for help on using the changeset viewer.