Changeset 24626


Ignore:
Timestamp:
03/03/20 11:37:10 (5 years ago)
Author:
Mathieu Morlighem
Message:

CHG: speeding up SMB model with reconstructed T and P

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

Legend:

Unmodified
Added
Removed

  • issm/trunk-jpl/src/c/analyses/SmbAnalysis.cpp

    r24583 r24626  
    114114                                iomodel->FetchDataToDatasetInput(inputs2,elements,"md.smb.precipitations_presentday",SmbPrecipitationsPresentdayEnum);
    115115                                if(!istemperaturescaled){
    116                                         iomodel->FetchDataToInput(inputs2,elements,"md.smb.temperatures_reconstructed",SmbTemperaturesReconstructedEnum);
     116                                        /*Fetch array*/
     117                                        IssmDouble* doublearray = NULL;
     118                                        int         M,N;
     119                                        iomodel->FetchData(&doublearray,&M,&N,"md.smb.temperatures_reconstructed");
     120                                        if(M!=iomodel->numberofvertices+1) _error_("md.smb.temperatures_reconstructed should have nbv+1 rows");
     121                                        if(N%12!=0) _error_("md.smb.temperatures_reconstructed should have a multiple of 12 columns (since it is monthly)");
     122
     123                                        /*Build indices*/
     124                                        int* ids = xNew<int>(N); for(int i=0;i<N;i++) ids[i] = i;
     125
     126                                        for(int i=0;i<elements->Size();i++){
     127                                                Element* element=xDynamicCast<Element*>(elements->GetObjectByOffset(i));
     128                                                element->DatasetInputCreate(doublearray,M-1,N,ids,N,inputs2,iomodel,SmbTemperaturesReconstructedEnum);
     129                                        }
     130                                        xDelete<int>(ids);
     131                                        iomodel->DeleteData(doublearray,"md.smb.temperatures_reconstructed");
    117132                                }
    118133                                if(!isprecipscaled){
    119                                         iomodel->FetchDataToInput(inputs2,elements,"md.smb.precipitations_reconstructed",SmbPrecipitationsReconstructedEnum);
     134                                        /*Fetch array*/
     135                                        IssmDouble* doublearray = NULL;
     136                                        int         M,N;
     137                                        iomodel->FetchData(&doublearray,&M,&N,"md.smb.precipitations_reconstructed");
     138                                        if(M!=iomodel->numberofvertices+1) _error_("md.smb.precipitations_reconstructed should have nbv+1 rows");
     139                                        if(N%12!=0) _error_("md.smb.precipitations_reconstructed should have a multiple of 12 columns (since it is monthly)");
     140
     141                                        /*Build indices*/
     142                                        int* ids = xNew<int>(N); for(int i=0;i<N;i++) ids[i] = i;
     143
     144                                        for(int i=0;i<elements->Size();i++){
     145                                                Element* element=xDynamicCast<Element*>(elements->GetObjectByOffset(i));
     146                                                element->DatasetInputCreate(doublearray,M-1,N,ids,N,inputs2,iomodel,SmbPrecipitationsReconstructedEnum);
     147                                        }
     148                                        xDelete<int>(ids);
     149                                        iomodel->DeleteData(doublearray,"md.smb.precipitations_reconstructed");
    120150                                }
    121151                        }
    … …  
    248278                                iomodel->DeleteData(temp,"md.smb.delta18o_surface");
    249279                        }
     280
    250281                        break;
    251282                case SMBpddSicopolisEnum:
    … …  
    266297                                parameters->AddObject(new TransientParam(SmbDelta18oEnum,&temp[0],&temp[M],interp,M));
    267298                                iomodel->DeleteData(temp,"md.smb.delta18o");
    268                         }
     299
     300                                IssmDouble yts;
     301                                bool istemperaturescaled,isprecipscaled;
     302                                iomodel->FindConstant(&yts,"md.constants.yts");
     303                                iomodel->FindConstant(&istemperaturescaled,"md.smb.istemperaturescaled");
     304                                iomodel->FindConstant(&isprecipscaled,"md.smb.isprecipscaled");
     305                                if(!istemperaturescaled){
     306                                        /*Fetch array*/
     307                                        IssmDouble* doublearray = NULL;
     308                                        int         M,N;
     309                                        iomodel->FetchData(&doublearray,&M,&N,"md.smb.temperatures_reconstructed");
     310                                        if(M!=iomodel->numberofvertices+1) _error_("md.smb.temperatures_reconstructed should have nbv+1 rows");
     311                                        if(N%12!=0) _error_("md.smb.temperatures_reconstructed should have a multiple of 12 columns (since it is monthly)");
     312                                        int numyears = N/12; _assert_(numyears*12==N);
     313
     314                                        /*Check times*/
     315                                        #ifdef _ISSM_DEBUG_
     316                                        for(int i=0;i<numyears;i++){
     317                                                for(int j=1;j<12;j++){
     318                                                        //_assert_(floor(doublearray[(M-1)*N+i*12+j]/yts)==floor(doublearray[(M-1)*N+i*12]/yts));
     319                                                        _assert_(doublearray[(M-1)*N+i*12+j]>doublearray[(M-1)*N+i*12+j-1]);
     320                                                }
     321                                        }
     322                                        #endif
     323
     324                                        /*Build time*/
     325                                        IssmDouble* times = xNew<IssmDouble>(numyears); for(int i=0;i<numyears;i++) times[i] = doublearray[(M-1)*N+i*12];
     326                                        parameters->AddObject(new DoubleVecParam(SmbTemperaturesReconstructedYearsEnum,times,numyears));
     327                                        xDelete<IssmDouble>(times);
     328                                        iomodel->DeleteData(doublearray,"md.smb.temperatures_reconstructed");
     329                                }
     330                                if(!isprecipscaled){
     331                                        /*Fetch array*/
     332                                        IssmDouble* doublearray = NULL;
     333                                        int         M,N;
     334                                        iomodel->FetchData(&doublearray,&M,&N,"md.smb.precipitations_reconstructed");
     335                                        if(M!=iomodel->numberofvertices+1) _error_("md.smb.precipitations_reconstructed should have nbv+1 rows");
     336                                        if(N%12!=0) _error_("md.smb.precipitations_reconstructed should have a multiple of 12 columns (since it is monthly)");
     337                                        int numyears = N/12; _assert_(numyears*12==N);
     338
     339                                        /*Check times*/
     340                                        #ifdef _ISSM_DEBUG_
     341                                        for(int i=0;i<numyears;i++){
     342                                                for(int j=1;j<12;j++){
     343                                                        //_assert_(floor(doublearray[(M-1)*N+i*12+j]/yts)==floor(doublearray[(M-1)*N+i*12]/yts));
     344                                                        _assert_(doublearray[(M-1)*N+i*12+j]>doublearray[(M-1)*N+i*12+j-1]);
     345                                                }
     346                                        }
     347                                        #endif
     348
     349                                        /*Build time*/
     350                                        IssmDouble* times = xNew<IssmDouble>(numyears); for(int i=0;i<numyears;i++) times[i] = doublearray[(M-1)*N+i*12];
     351                                        parameters->AddObject(new DoubleVecParam(SmbPrecipitationsReconstructedYearsEnum,times,numyears));
     352                                        xDelete<IssmDouble>(times);
     353                                        iomodel->DeleteData(doublearray,"md.smb.precipitations_reconstructed");
     354                                }
     355                        }
     356
    269357                        break;
    270358                case SMBgradientsEnum:

  • issm/trunk-jpl/src/c/classes/Elements/Element.cpp

    r24583 r24626  
    563563
    564564                /*Recover present day temperature and precipitation*/
    565         TransientInput2 *tinput3 = NULL;
    566         TransientInput2 *tinput4 = NULL;
     565        DatasetInput2 *dinput3 = NULL;
     566        DatasetInput2 *dinput4 = NULL;
     567        int            offset_t,offset_p,N;
    567568        if(!isTemperatureScaled){
    568                 tinput3=this->inputs2->GetTransientInput(SmbTemperaturesReconstructedEnum); _assert_(tinput3);
    569                 offset     = tinput3->GetTimeInputOffset(time_climt);
    570                 time_climt = tinput3->GetTimeByOffset(offset-11-fmod(offset-11,12.));
     569                IssmDouble* time_temp_scaled = NULL;
     570                parameters->FindParam(&time_temp_scaled,&N,SmbTemperaturesReconstructedYearsEnum);
     571                if(!binary_search(&offset_t,time_climt,time_temp_scaled,N)) _error_("time not sorted?");
     572                xDelete<IssmDouble>(time_temp_scaled);
     573                dinput3=this->GetDatasetInput2(SmbTemperaturesReconstructedEnum); _assert_(dinput3);
    571574        }
    572575        if(!isPrecipScaled){
    573                 tinput4=this->inputs2->GetTransientInput(SmbPrecipitationsReconstructedEnum); _assert_(tinput4);
    574                 offset     = tinput4->GetTimeInputOffset(time_climp);
    575                 time_climp = tinput4->GetTimeByOffset(offset-11-fmod(offset-11,12.));
     576                IssmDouble* time_precip_scaled = NULL;
     577                parameters->FindParam(&time_precip_scaled,&N,SmbPrecipitationsReconstructedYearsEnum);
     578                if(!binary_search(&offset_t,time_climt,time_precip_scaled,N)) _error_("time not sorted?");
     579                xDelete<IssmDouble>(time_precip_scaled);
     580                dinput4=this->GetDatasetInput2(SmbPrecipitationsReconstructedEnum); _assert_(dinput4);
    576581        }
    577582
    … …  
    583588        Gauss* gauss=this->NewGauss();
    584589        for(int month=0;month<12;month++) {
    585 
    586                 Input2 *input3 = NULL;
    587                 Input2 *input4 = NULL;
    588                 if(!isTemperatureScaled){
    589                         input3=this->GetInput2(SmbTemperaturesReconstructedEnum,time_climt+month/12.*yts); _assert_(input3);
    590                 }
    591                 if(!isPrecipScaled){
    592                         input4=this->GetInput2(SmbPrecipitationsReconstructedEnum,time_climp+month/12.*yts);
    593                 }
    594590                for(int iv=0;iv<NUM_VERTICES;iv++) {
    595591                        gauss->GaussVertex(iv);
    … …  
    599595
    600596                        if(!isTemperatureScaled){
    601                                 input3->GetInputValue(&TemperaturesReconstructed[iv*12+month],gauss);
     597                                dinput3->GetInputValue(&TemperaturesReconstructed[iv*12+month],gauss,offset_t*12+month);
    602598                        }
    603599                        if(!isPrecipScaled){
    604                                 input4->GetInputValue(&PrecipitationsReconstructed[iv*12+month],gauss);
     600                                dinput4->GetInputValue(&PrecipitationsReconstructed[iv*12+month],gauss,offset_p*12+month);
    605601                                PrecipitationsReconstructed[iv*12+month]=PrecipitationsReconstructed[iv*12+month]*yts;
    606602                        }

  • issm/trunk-jpl/src/c/shared/Enum/Enum.vim

    r24548 r24626  
    383383syn keyword cConstant SmbTdiffEnum
    384384syn keyword cConstant SmbThermoDeltaTScalingEnum
     385syn keyword cConstant SmbTemperaturesReconstructedYearsEnum
     386syn keyword cConstant SmbPrecipitationsReconstructedYearsEnum
    385387syn keyword cConstant SmoothThicknessMultiplierEnum
    386388syn keyword cConstant SolutionTypeEnum

  • issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h

    r24548 r24626  
    377377        SmbTdiffEnum,
    378378        SmbThermoDeltaTScalingEnum,
     379        SmbTemperaturesReconstructedYearsEnum,
     380        SmbPrecipitationsReconstructedYearsEnum,
    379381        SmoothThicknessMultiplierEnum,
    380382        SolutionTypeEnum,

  • issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp

    r24548 r24626  
    385385                case SmbTdiffEnum : return "SmbTdiff";
    386386                case SmbThermoDeltaTScalingEnum : return "SmbThermoDeltaTScaling";
     387                case SmbTemperaturesReconstructedYearsEnum : return "SmbTemperaturesReconstructedYears";
     388                case SmbPrecipitationsReconstructedYearsEnum : return "SmbPrecipitationsReconstructedYears";
    387389                case SmoothThicknessMultiplierEnum : return "SmoothThicknessMultiplier";
    388390                case SolutionTypeEnum : return "SolutionType";

  • issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp

    r24548 r24626  
    394394              else if (strcmp(name,"SmbTdiff")==0) return SmbTdiffEnum;
    395395              else if (strcmp(name,"SmbThermoDeltaTScaling")==0) return SmbThermoDeltaTScalingEnum;
     396              else if (strcmp(name,"SmbTemperaturesReconstructedYears")==0) return SmbTemperaturesReconstructedYearsEnum;
     397              else if (strcmp(name,"SmbPrecipitationsReconstructedYears")==0) return SmbPrecipitationsReconstructedYearsEnum;
    396398              else if (strcmp(name,"SmoothThicknessMultiplier")==0) return SmoothThicknessMultiplierEnum;
    397399              else if (strcmp(name,"SolutionType")==0) return SolutionTypeEnum;
    … …  
    504506              else if (strcmp(name,"CalvinglevermannCoeff")==0) return CalvinglevermannCoeffEnum;
    505507              else if (strcmp(name,"CalvingratexAverage")==0) return CalvingratexAverageEnum;
    506               else if (strcmp(name,"Calvingratex")==0) return CalvingratexEnum;
    507               else if (strcmp(name,"CalvingrateyAverage")==0) return CalvingrateyAverageEnum;
    508508         else stage=5;
    509509   }
    510510   if(stage==5){
    511               if (strcmp(name,"Calvingratey")==0) return CalvingrateyEnum;
     511              if (strcmp(name,"Calvingratex")==0) return CalvingratexEnum;
     512              else if (strcmp(name,"CalvingrateyAverage")==0) return CalvingrateyAverageEnum;
     513              else if (strcmp(name,"Calvingratey")==0) return CalvingrateyEnum;
    512514              else if (strcmp(name,"CalvingFluxLevelset")==0) return CalvingFluxLevelsetEnum;
    513515              else if (strcmp(name,"CalvingMeltingFluxLevelset")==0) return CalvingMeltingFluxLevelsetEnum;
    … …  
    627629              else if (strcmp(name,"InversionVelObs")==0) return InversionVelObsEnum;
    628630              else if (strcmp(name,"InversionVxObs")==0) return InversionVxObsEnum;
    629               else if (strcmp(name,"InversionVyObs")==0) return InversionVyObsEnum;
    630               else if (strcmp(name,"LevelsetfunctionSlopeX")==0) return LevelsetfunctionSlopeXEnum;
    631631         else stage=6;
    632632   }
    633633   if(stage==6){
    634               if (strcmp(name,"LevelsetfunctionSlopeY")==0) return LevelsetfunctionSlopeYEnum;
     634              if (strcmp(name,"InversionVyObs")==0) return InversionVyObsEnum;
     635              else if (strcmp(name,"LevelsetfunctionSlopeX")==0) return LevelsetfunctionSlopeXEnum;
     636              else if (strcmp(name,"LevelsetfunctionSlopeY")==0) return LevelsetfunctionSlopeYEnum;
    635637              else if (strcmp(name,"LoadingforceX")==0) return LoadingforceXEnum;
    636638              else if (strcmp(name,"LoadingforceY")==0) return LoadingforceYEnum;
    … …  
    750752              else if (strcmp(name,"SmbMeanSHF")==0) return SmbMeanSHFEnum;
    751753              else if (strcmp(name,"SmbMeanULW")==0) return SmbMeanULWEnum;
    752               else if (strcmp(name,"SmbMelt")==0) return SmbMeltEnum;
    753               else if (strcmp(name,"SmbMonthlytemperatures")==0) return SmbMonthlytemperaturesEnum;
    754754         else stage=7;
    755755   }
    756756   if(stage==7){
    757               if (strcmp(name,"SmbNetLW")==0) return SmbNetLWEnum;
     757              if (strcmp(name,"SmbMelt")==0) return SmbMeltEnum;
     758              else if (strcmp(name,"SmbMonthlytemperatures")==0) return SmbMonthlytemperaturesEnum;
     759              else if (strcmp(name,"SmbNetLW")==0) return SmbNetLWEnum;
    758760              else if (strcmp(name,"SmbNetSW")==0) return SmbNetSWEnum;
    759761              else if (strcmp(name,"SmbPAir")==0) return SmbPAirEnum;
    … …  
    873875              else if (strcmp(name,"Outputdefinition20")==0) return Outputdefinition20Enum;
    874876              else if (strcmp(name,"Outputdefinition21")==0) return Outputdefinition21Enum;
    875               else if (strcmp(name,"Outputdefinition22")==0) return Outputdefinition22Enum;
    876               else if (strcmp(name,"Outputdefinition23")==0) return Outputdefinition23Enum;
    877877         else stage=8;
    878878   }
    879879   if(stage==8){
    880               if (strcmp(name,"Outputdefinition24")==0) return Outputdefinition24Enum;
     880              if (strcmp(name,"Outputdefinition22")==0) return Outputdefinition22Enum;
     881              else if (strcmp(name,"Outputdefinition23")==0) return Outputdefinition23Enum;
     882              else if (strcmp(name,"Outputdefinition24")==0) return Outputdefinition24Enum;
    881883              else if (strcmp(name,"Outputdefinition25")==0) return Outputdefinition25Enum;
    882884              else if (strcmp(name,"Outputdefinition26")==0) return Outputdefinition26Enum;
    … …  
    996998              else if (strcmp(name,"Boundary")==0) return BoundaryEnum;
    997999              else if (strcmp(name,"BuddJacka")==0) return BuddJackaEnum;
    998               else if (strcmp(name,"CalvingDev2")==0) return CalvingDev2Enum;
    999               else if (strcmp(name,"CalvingHab")==0) return CalvingHabEnum;
    10001000         else stage=9;
    10011001   }
    10021002   if(stage==9){
    1003               if (strcmp(name,"CalvingLevermann")==0) return CalvingLevermannEnum;
     1003              if (strcmp(name,"CalvingDev2")==0) return CalvingDev2Enum;
     1004              else if (strcmp(name,"CalvingHab")==0) return CalvingHabEnum;
     1005              else if (strcmp(name,"CalvingLevermann")==0) return CalvingLevermannEnum;
    10041006              else if (strcmp(name,"CalvingVonmises")==0) return CalvingVonmisesEnum;
    10051007              else if (strcmp(name,"Cfdragcoeffabsgrad")==0) return CfdragcoeffabsgradEnum;
    … …  
    11191121              else if (strcmp(name,"Indexed")==0) return IndexedEnum;
    11201122              else if (strcmp(name,"IntExternalResult")==0) return IntExternalResultEnum;
    1121               else if (strcmp(name,"IntInput")==0) return IntInputEnum;
    1122               else if (strcmp(name,"ElementInput2")==0) return ElementInput2Enum;
    11231123         else stage=10;
    11241124   }
    11251125   if(stage==10){
    1126               if (strcmp(name,"SegInput2")==0) return SegInput2Enum;
     1126              if (strcmp(name,"IntInput")==0) return IntInputEnum;
     1127              else if (strcmp(name,"ElementInput2")==0) return ElementInput2Enum;
     1128              else if (strcmp(name,"SegInput2")==0) return SegInput2Enum;
    11271129              else if (strcmp(name,"TriaInput2")==0) return TriaInput2Enum;
    11281130              else if (strcmp(name,"PentaInput2")==0) return PentaInput2Enum;
    … …  
    12421244              else if (strcmp(name,"Riftfront")==0) return RiftfrontEnum;
    12431245              else if (strcmp(name,"SIAApproximation")==0) return SIAApproximationEnum;
    1244               else if (strcmp(name,"SMBcomponents")==0) return SMBcomponentsEnum;
    1245               else if (strcmp(name,"SMBd18opdd")==0) return SMBd18opddEnum;
    12461246         else stage=11;
    12471247   }
    12481248   if(stage==11){
    1249               if (strcmp(name,"SMBforcing")==0) return SMBforcingEnum;
     1249              if (strcmp(name,"SMBcomponents")==0) return SMBcomponentsEnum;
     1250              else if (strcmp(name,"SMBd18opdd")==0) return SMBd18opddEnum;
     1251              else if (strcmp(name,"SMBforcing")==0) return SMBforcingEnum;
    12501252              else if (strcmp(name,"SMBgcm")==0) return SMBgcmEnum;
    12511253              else if (strcmp(name,"SMBgemb")==0) return SMBgembEnum;
Note: See TracChangeset for help on using the changeset viewer.