Changeset 18872


Ignore:
Timestamp:
11/28/14 09:42:15 (10 years ago)
Author:
Eric.Larour
Message:

CHG: new constructor for the ad control core solution, which allows to feed in an X independent vector and
a choice on where to start tracing or not. This is useful to allow updates on the control update through
feeding this update during the I/O phase of the FeMmodel constructor. This way, we make sure that for each
ad run of the control core, the tape initializes with the updated X vector instead of the value from disk.

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

Legend:

Unmodified
Added
Removed

  • issm/trunk-jpl/src/c/classes/FemModel.cpp

    r18852 r18872  
    4343
    4444/*Object constructors and destructor*/
    45 FemModel::FemModel(int argc,char** argv,ISSM_MPI_Comm incomm){/*{{{*/
     45FemModel::FemModel(int argc,char** argv,ISSM_MPI_Comm incomm,bool trace){/*{{{*/
    4646
    4747        /*configuration: */
    48         int* analyses=NULL;
    49         int  numanalyses;
    5048        int  solution_type;
    5149        int  ierr;
    … …  
    7775        /*Create femmodel from input files: */
    7876        profiler->Tag(StartInit);
    79         this->InitFromFiles(rootpath,binfilename,outbinfilename,petscfilename,lockfilename,solution_type);
     77        this->InitFromFiles(rootpath,binfilename,outbinfilename,petscfilename,lockfilename,solution_type,trace,NULL);
    8078        profiler->Tag(FinishInit);
    8179
    8280        /*Free resources */
    83         xDelete<int>(analyses);
    8481        xDelete<char>(lockfilename);
    8582        xDelete<char>(binfilename);
    … …  
    8784        xDelete<char>(petscfilename);
    8885        xDelete<char>(rootpath);
     86
     87}
     88/*}}}*/
     89FemModel::FemModel(char* rootpath, char* inputfilename, char* outputfilename, char* toolkitsfilename, char* lockfilename, ISSM_MPI_Comm incomm, int solution_type,IssmPDouble* X){ /*{{{*/
     90
     91        bool traceon=true;
     92       
     93        /*Store the communicator, but do not set it as a global variable, as this has already
     94         * been done by the FemModel that called this copy constructor: */
     95        this->comm=incomm;
     96        this->SetStaticComm();
     97
     98        /*Create femmodel from input files, with trace activated: */
     99        this->InitFromFiles(rootpath,inputfilename,outputfilename,toolkitsfilename,lockfilename,solution_type,traceon,X);
    89100
    90101}
    … …  
    131142}
    132143/*}}}*/
    133 void FemModel::InitFromFiles(char* rootpath, char* inputfilename, char* outputfilename, char* toolkitsfilename, char* lockfilename, const int in_solution_type){/*{{{*/
     144void FemModel::InitFromFiles(char* rootpath, char* inputfilename, char* outputfilename, char* toolkitsfilename, char* lockfilename, const int in_solution_type,bool trace,IssmPDouble* X){/*{{{*/
    134145
    135146        /*intermediary*/
    … …  
    148159        this->analysis_counter = nummodels-1;   //point to last analysis_type carried out.
    149160        this->results          = new Results(); //not initialized by CreateDataSets
    150 
    151161        /*Open input file on cpu 0 and create IoModel */
    152162        if(my_rank==0) IOMODEL = pfopen0(inputfilename ,"rb");
    153         IoModel* iomodel = new IoModel(IOMODEL);
     163        IoModel* iomodel = new IoModel(IOMODEL,trace,X);
    154164
    155165        /*Figure out what analyses are activated for this solution*/
    … …  
    188198        pfclose(toolkitsoptionsfid,toolkitsfilename);
    189199
    190         /*Open output file once for all and add output file name and file descriptor to parameters*/
     200        /*Open output file once for all and add output file descriptor to parameters*/
    191201        output_fid=pfopen(outputfilename,"wb");
     202        this->parameters->SetParam(output_fid,OutputFilePointerEnum);
     203       
     204        /*Now save all of these file names into parameters, you never know when you might need them: */
     205        this->parameters->AddObject(new StringParam(ToolkitsFileNameEnum,toolkitsfilename));
     206        this->parameters->AddObject(new StringParam(RootPathEnum,rootpath));
     207        this->parameters->AddObject(new StringParam(InputFileNameEnum,inputfilename));
    192208        this->parameters->AddObject(new StringParam(OutputFileNameEnum,outputfilename));
    193         this->parameters->SetParam(output_fid,OutputFilePointerEnum);
    194 
    195         /*Save lock file name for later: */
    196209        this->parameters->AddObject(new StringParam(LockFileNameEnum,lockfilename));
    197210
    … …  
    741754                }
    742755                else{
    743                         switch(output_enum){
    744 
    745                                 /*Scalar output*/
    746                                 case DivergenceEnum:               this->Divergencex(&double_result);               break;
    747                                 case MaxDivergenceEnum:            this->MaxDivergencex(&double_result);            break;
    748                                 case IceVolumeEnum:                this->IceVolumex(&double_result);                break;
    749                                 case IceVolumeAboveFloatationEnum: this->IceVolumeAboveFloatationx(&double_result); break;
    750                                 case MinVelEnum:                   this->MinVelx(&double_result);                   break;
    751                                 case MaxVelEnum:                   this->MaxVelx(&double_result);                   break;
    752                                 case MinVxEnum:                    this->MinVxx(&double_result);                    break;
    753                                 case MaxVxEnum:                    this->MaxVxx(&double_result);                    break;
    754                                 case MaxAbsVxEnum:                 this->MaxAbsVxx(&double_result);                 break;
    755                                 case MinVyEnum:                    this->MinVyx(&double_result);                    break;
    756                                 case MaxVyEnum:                    this->MaxVyx(&double_result);                    break;
    757                                 case MaxAbsVyEnum:                 this->MaxAbsVyx(&double_result);                 break;
    758                                 case MinVzEnum:                    this->MinVzx(&double_result);                    break;
    759                                 case MaxVzEnum:                    this->MaxVzx(&double_result);                    break;
    760                                 case MaxAbsVzEnum:                 this->MaxAbsVzx(&double_result);                 break;
    761                                 case MassFluxEnum:                 this->MassFluxx(&double_result);                 break;
    762                                 case TotalSmbEnum:                 this->TotalSmbx(&double_result);                 break;
     756                        /*last chance for the output definition, if the enum is one of Outputdefinition[1-10]Enum:*/
     757                        if(output_enum>=Outputdefinition1Enum && output_enum <=Outputdefinition10Enum){
     758                                double_result = OutputDefinitionsResponsex(this,output_enum);
     759                                if(save_results){
     760                                        results->AddResult(new GenericExternalResult<IssmPDouble>(results->Size()+1,output_string,reCast<IssmPDouble>(double_result),step,time));
     761                                        continue;
     762                                }
     763                        }
     764                        else{
     765                                switch(output_enum){
     766
     767                                        /*Scalar output*/
     768                                        case DivergenceEnum:               this->Divergencex(&double_result);               break;
     769                                        case MaxDivergenceEnum:            this->MaxDivergencex(&double_result);            break;
     770                                        case IceVolumeEnum:                this->IceVolumex(&double_result);                break;
     771                                        case IceVolumeAboveFloatationEnum: this->IceVolumeAboveFloatationx(&double_result); break;
     772                                        case MinVelEnum:                   this->MinVelx(&double_result);                   break;
     773                                        case MaxVelEnum:                   this->MaxVelx(&double_result);                   break;
     774                                        case MinVxEnum:                    this->MinVxx(&double_result);                    break;
     775                                        case MaxVxEnum:                    this->MaxVxx(&double_result);                    break;
     776                                        case MaxAbsVxEnum:                 this->MaxAbsVxx(&double_result);                 break;
     777                                        case MinVyEnum:                    this->MinVyx(&double_result);                    break;
     778                                        case MaxVyEnum:                    this->MaxVyx(&double_result);                    break;
     779                                        case MaxAbsVyEnum:                 this->MaxAbsVyx(&double_result);                 break;
     780                                        case MinVzEnum:                    this->MinVzx(&double_result);                    break;
     781                                        case MaxVzEnum:                    this->MaxVzx(&double_result);                    break;
     782                                        case MaxAbsVzEnum:                 this->MaxAbsVzx(&double_result);                 break;
     783                                        case MassFluxEnum:                 this->MassFluxx(&double_result);                 break;
     784                                        case TotalSmbEnum:                 this->TotalSmbx(&double_result);                 break;
    763785
    764786                           /*Scalar control output*/
    … …  
    778800                                case SurfaceAbsMisfitEnum:          SurfaceAbsMisfitx(&double_result); break;
    779801
    780                            /*Vector */
    781                                 default:
    782 
    783                                         /*Vector layout*/
    784                                         int interpolation,nodesperelement,size;
    785                                         int rank_interpolation=-1,rank_nodesperelement=-1;
    786 
    787                                         /*Get interpolation (and compute input if necessary)*/
    788                                         for(int j=0;j<elements->Size();j++){
    789                                                 Element* element=xDynamicCast<Element*>(this->elements->GetObjectByOffset(j));
    790                                                 element->ResultInterpolation(&rank_interpolation,&rank_nodesperelement,output_enum);
    791                                         }
    792 
    793                                         /*Broadcast for cpus that do not have any elements*/
    794                                         ISSM_MPI_Reduce(&rank_interpolation,&interpolation,1,ISSM_MPI_INT,ISSM_MPI_MAX,0,IssmComm::GetComm());
    795                                         ISSM_MPI_Reduce(&rank_nodesperelement,&nodesperelement,1,ISSM_MPI_INT,ISSM_MPI_MAX,0,IssmComm::GetComm());
    796                                         ISSM_MPI_Bcast(&interpolation,1,ISSM_MPI_INT,0,IssmComm::GetComm());
    797                                         ISSM_MPI_Bcast(&nodesperelement,1,ISSM_MPI_INT,0,IssmComm::GetComm());
    798 
    799                                         if(results_on_nodes){
    800 
    801                                                 /*Allocate matrices*/
    802                                                 int         nbe       = this->elements->NumberOfElements();
    803                                                 IssmDouble* values    = xNewZeroInit<IssmDouble>(nbe*nodesperelement);
    804                                                 IssmDouble* allvalues = xNew<IssmDouble>(nbe*nodesperelement);
    805 
    806                                                 /*Fill-in matrix*/
     802                                   /*Vector */
     803                                        default:
     804
     805                                                /*Vector layout*/
     806                                                int interpolation,nodesperelement,size;
     807                                                int rank_interpolation=-1,rank_nodesperelement=-1;
     808
     809                                                /*Get interpolation (and compute input if necessary)*/
    807810                                                for(int j=0;j<elements->Size();j++){
    808811                                                        Element* element=xDynamicCast<Element*>(this->elements->GetObjectByOffset(j));
    809                                                         element->ResultToPatch(values,nodesperelement,output_enum);
     812                                                        element->ResultInterpolation(&rank_interpolation,&rank_nodesperelement,output_enum);
    810813                                                }
    811814
    812                                                 /*Gather from all cpus*/
    813                                                 ISSM_MPI_Allreduce((void*)values,(void*)allvalues,nbe*nodesperelement,ISSM_MPI_PDOUBLE,ISSM_MPI_SUM,IssmComm::GetComm());
    814                                                 xDelete<IssmDouble>(values);
    815 
    816                                                 if(save_results)results->AddResult(new GenericExternalResult<IssmDouble*>(results->Size()+1,output_enum,allvalues,nbe,nodesperelement,step,time));
    817                                                 xDelete<IssmDouble>(allvalues);
    818 
    819                                         }
    820                                         else{
    821 
    822                                                 /*Allocate vector depending on interpolation*/
    823                                                 switch(interpolation){
    824                                                         case P0Enum: size = this->elements->NumberOfElements(); break;
    825                                                         case P1Enum: size = this->vertices->NumberOfVertices(); break;
    826                                                         default:     _error_("Interpolation "<<EnumToStringx(interpolation)<<" not supported yet");
     815                                                /*Broadcast for cpus that do not have any elements*/
     816                                                ISSM_MPI_Reduce(&rank_interpolation,&interpolation,1,ISSM_MPI_INT,ISSM_MPI_MAX,0,IssmComm::GetComm());
     817                                                ISSM_MPI_Reduce(&rank_nodesperelement,&nodesperelement,1,ISSM_MPI_INT,ISSM_MPI_MAX,0,IssmComm::GetComm());
     818                                                ISSM_MPI_Bcast(&interpolation,1,ISSM_MPI_INT,0,IssmComm::GetComm());
     819                                                ISSM_MPI_Bcast(&nodesperelement,1,ISSM_MPI_INT,0,IssmComm::GetComm());
     820
     821                                                if(results_on_nodes){
     822
     823                                                        /*Allocate matrices*/
     824                                                        int         nbe       = this->elements->NumberOfElements();
     825                                                        IssmDouble* values    = xNewZeroInit<IssmDouble>(nbe*nodesperelement);
     826                                                        IssmDouble* allvalues = xNew<IssmDouble>(nbe*nodesperelement);
     827
     828                                                        /*Fill-in matrix*/
     829                                                        for(int j=0;j<elements->Size();j++){
     830                                                                Element* element=xDynamicCast<Element*>(this->elements->GetObjectByOffset(j));
     831                                                                element->ResultToPatch(values,nodesperelement,output_enum);
     832                                                        }
     833
     834                                                        /*Gather from all cpus*/
     835                                                        ISSM_MPI_Allreduce((void*)values,(void*)allvalues,nbe*nodesperelement,ISSM_MPI_PDOUBLE,ISSM_MPI_SUM,IssmComm::GetComm());
     836                                                        xDelete<IssmDouble>(values);
     837
     838                                                        if(save_results)results->AddResult(new GenericExternalResult<IssmDouble*>(results->Size()+1,output_enum,allvalues,nbe,nodesperelement,step,time));
     839                                                        xDelete<IssmDouble>(allvalues);
    827840
    828841                                                }
    829                                                 Vector<IssmDouble> *vector_result = new Vector<IssmDouble>(size);
    830 
    831                                                 /*Fill in vector*/
    832                                                 for(int j=0;j<elements->Size();j++){
    833                                                         Element* element=(Element*)elements->GetObjectByOffset(j);
    834                                                         element->ResultToVector(vector_result,output_enum);
     842                                                else{
     843
     844                                                        /*Allocate vector depending on interpolation*/
     845                                                        switch(interpolation){
     846                                                                case P0Enum: size = this->elements->NumberOfElements(); break;
     847                                                                case P1Enum: size = this->vertices->NumberOfVertices(); break;
     848                                                                default:     _error_("Interpolation "<<EnumToStringx(interpolation)<<" not supported yet");
     849
     850                                                        }
     851                                                        Vector<IssmDouble> *vector_result = new Vector<IssmDouble>(size);
     852
     853                                                        /*Fill in vector*/
     854                                                        for(int j=0;j<elements->Size();j++){
     855                                                                Element* element=(Element*)elements->GetObjectByOffset(j);
     856                                                                element->ResultToVector(vector_result,output_enum);
     857                                                        }
     858                                                        vector_result->Assemble();
     859
     860                                                        if(save_results)results->AddResult(new GenericExternalResult<Vector<IssmDouble>*>(results->Size()+1,output_enum,vector_result,step,time));
    835861                                                }
    836                                                 vector_result->Assemble();
    837 
    838                                                 if(save_results)results->AddResult(new GenericExternalResult<Vector<IssmDouble>*>(results->Size()+1,output_enum,vector_result,step,time));
    839                                         }
    840                                         isvec = true;
    841                                         break;
    842                         }
     862                                                isvec = true;
     863                                                break;
     864                                }
     865                        }
     866
    843867                }
    844868

  • issm/trunk-jpl/src/c/classes/FemModel.h

    r18736 r18872  
    4646
    4747                /*constructors, destructors: */
    48                 FemModel(int argc,char** argv,ISSM_MPI_Comm comm_init);
     48                FemModel(int argc,char** argv,ISSM_MPI_Comm comm_init,bool trace=false);
     49                FemModel(char* rootpath, char* inputfilename, char* outputfilename, char* toolkitsfilename, char* lockfilename, ISSM_MPI_Comm incomm, int solution_type,IssmPDouble* X);
    4950                ~FemModel();
    5051
    … …  
    5253                void Echo();
    5354                FemModel* copy();
    54                 void InitFromFiles(char* rootpath, char* inputfilename, char* outputfilename, char* petscfilename, char* lockfilename, const int solution_type);
     55                void InitFromFiles(char* rootpath, char* inputfilename, char* outputfilename, char* petscfilename, char* lockfilename, const int solution_type,bool trace,IssmPDouble* X=NULL);
    5556                void SolutionAnalysesList(int** panalyses,int* pnumanalyses,IoModel* iomodel,int solutiontype);
    5657                void CleanUp(void);
Note: See TracChangeset for help on using the changeset viewer.