Ignore:
Timestamp:
05/10/18 10:24:27 (7 years ago)
Author:
Mathieu Morlighem
Message:

merged trunk-jpl and trunk for revision 22757

Location:
issm/trunk
Files:
4 edited

Legend:

Unmodified
Added
Removed

  • issm/trunk

  • issm/trunk/src

  • issm/trunk/src/c

    • Property svn:ignore

      •  

        old new  
        2020kriging
        2121issm_slr
         22issm_ocean
         23lnb_param.mod
         24lovenb_sub.mod
         25model.mod
         26util.mod

  • issm/trunk/src/c/classes/Inputs/TransientInput.cpp

    r21341 r22758  
    7070        output->numtimesteps=this->numtimesteps;
    7171        output->timesteps=xNew<IssmDouble>(this->numtimesteps);
    72         xMemCpy(output->timesteps,this->timesteps,this->numtimesteps);
    73         output->inputs=static_cast<Inputs*>(this->inputs->Copy());
     72        if(this->numtimesteps>0){
     73                xMemCpy(output->timesteps,this->timesteps,this->numtimesteps);
     74                output->inputs=static_cast<Inputs*>(this->inputs->Copy());
     75        }
    7476        output->parameters=this->parameters;
    7577
    … …  
    133135        outinput->enum_type=this->enum_type;
    134136        outinput->numtimesteps=this->numtimesteps;
    135         outinput->timesteps=xNew<IssmDouble>(this->numtimesteps);
    136         xMemCpy(outinput->timesteps,this->timesteps,this->numtimesteps);
    137         outinput->inputs=(Inputs*)this->inputs->SpawnSegInputs(index1,index2);
     137        if(this->numtimesteps>0){
     138                outinput->timesteps=xNew<IssmDouble>(this->numtimesteps);
     139                xMemCpy(outinput->timesteps,this->timesteps,this->numtimesteps);
     140                outinput->inputs=(Inputs*)this->inputs->SpawnSegInputs(index1,index2);
     141        }
    138142        outinput->parameters=this->parameters;
    139143
    … …  
    152156        outinput->enum_type=this->enum_type;
    153157        outinput->numtimesteps=this->numtimesteps;
    154         outinput->timesteps=xNew<IssmDouble>(this->numtimesteps);
    155         xMemCpy(outinput->timesteps,this->timesteps,this->numtimesteps);
    156         outinput->inputs=(Inputs*)this->inputs->SpawnTriaInputs(index1,index2,index3);
     158        if(this->numtimesteps>0){
     159                outinput->timesteps=xNew<IssmDouble>(this->numtimesteps);
     160                xMemCpy<IssmDouble>(outinput->timesteps,this->timesteps,this->numtimesteps);
     161                outinput->inputs=(Inputs*)this->inputs->SpawnTriaInputs(index1,index2,index3);
     162        }
    157163        outinput->parameters=this->parameters;
    158164
    … …  
    226232
    227233        /*First, recover current time from parameters: */
    228         this->parameters->FindParam(&time,TimeEnum);
     234        parameters->FindParam(&time,TimeEnum);
    229235
    230236        /*Retrieve interpolated values for this time step: */
    … …  
    248254}
    249255/*}}}*/
     256int  TransientInput::GetTimeInputOffset(IssmDouble time){/*{{{*/
     257
     258        int        found;
     259        int        offset;
     260
     261        /*go through the timesteps, and figure out which interval we
     262         *     *fall within. Then interpolate the values on this interval: */
     263        found=binary_search(&offset,time,this->timesteps,this->numtimesteps);
     264        if(!found) _error_("Input not found (is TransientInput sorted ?)");
     265
     266        return offset;
     267}
     268/*}}}*/
     269IssmDouble  TransientInput::GetTimeByOffset(int offset){/*{{{*/
     270        if (offset < 0) offset=0;
     271        _assert_(offset<(this->numtimesteps));
     272        return this->timesteps[offset];
     273}
     274/*}}}*/
    250275void TransientInput::GetInputUpToCurrentTimeAverages(IssmDouble** pvalues, IssmDouble** ptimes, int* pnumtimes, IssmDouble currenttime){/*{{{*/
    251276
    252277        int         i;
    253         IssmDouble *times  = NULL;
    254         IssmDouble *values = NULL;
     278        IssmDouble* times  = NULL;
     279        IssmDouble* values = NULL;
    255280        int         numsteps;
    256281        bool        iscurrenttime_included = false;
    … …  
    270295
    271296        for(i=0;i<numsteps;i++){
    272 
    273297                if((iscurrenttime_included==false) && (i==(numsteps-1))){
    274 
    275298                        /*Retrieve interpolated values for current time step: */
    276299                        Input* input=GetTimeInput(currenttime);
    … …  
    290313}
    291314/*}}}*/
     315void TransientInput::GetInputAverageOnTimes(IssmDouble** pvalues, IssmDouble init_time){/*{{{*/
     316
     317        int                                     numnodes;
     318        IssmDouble  subtimestep;
     319        IssmDouble* values= NULL;
     320        Input* input=NULL;
     321
     322        /*Initialize numnode from transientinput out of time loop*/
     323        for(int i=0;i<this->numtimesteps;i++){
     324                /*First compute the lengt of the sub-timestep*/
     325                if(i==0){
     326                        subtimestep = this->timesteps[i]-init_time;
     327                }
     328                else{
     329                        subtimestep = this->timesteps[i]-this->timesteps[i-1];
     330                }       
     331                /*Figure out type of input and get it*/
     332                input = xDynamicCast<Input*>(this->inputs->GetObjectByOffset(i)); _assert_(input);
     333                switch(input->ObjectEnum()){
     334                case TriaInputEnum:{
     335                        TriaInput* triainput = (TriaInput*)this->inputs->GetObjectByOffset(i); _assert_(triainput);
     336                        if(i==0){
     337                                numnodes= triainput->NumberofNodes(triainput->interpolation_type);
     338                                values=xNewZeroInit<IssmDouble>(numnodes);
     339                        }
     340                        /*Sum the values weighted by their respective sub-timestep*/
     341                        for(int j=0;j<numnodes;j++){
     342                                values[j]+=(triainput->values[j]*subtimestep);
     343                        }
     344                        break;
     345                }
     346                case PentaInputEnum:{
     347                        PentaInput*     pentainput = (PentaInput*)this->inputs->GetObjectByOffset(i); _assert_(pentainput);
     348                        if(i==0){
     349                                numnodes= pentainput->NumberofNodes(pentainput->interpolation_type);
     350                                values=xNewZeroInit<IssmDouble>(numnodes);
     351                        }
     352                        /*Sum the values weighted by their respective sub-timestep*/
     353                        for(int j=0;j<numnodes;j++){
     354                                values[j]+=(pentainput->values[j]*subtimestep);
     355                        }
     356                        break;
     357                }
     358                case TetraInputEnum:{
     359                        TetraInput*     tetrainput = (TetraInput*)this->inputs->GetObjectByOffset(i); _assert_(tetrainput);
     360                        if(i==0){
     361                                numnodes= tetrainput->NumberofNodes(tetrainput->interpolation_type);
     362                                values=xNewZeroInit<IssmDouble>(numnodes);
     363                        }
     364                        /*Sum the values weighted by their respective sub-timestep*/
     365                        for(int j=0;j<numnodes;j++){
     366                                values[j]+=(tetrainput->values[j]*subtimestep);
     367                        }
     368                        break;
     369                }
     370                default: _error_("not implemented yet");
     371                }
     372        }
     373        /*Compute the average based on the length of the full timestep*/
     374        for(int j=0;j<numnodes;j++){
     375                values[j]=values[j]/(this->timesteps[this->numtimesteps-1]-init_time);
     376        }
     377        *pvalues=values;
     378}
     379/*}}}*/
    292380
    293381/*Intermediary*/
    … …  
    379467        int        found;
    380468        int        offset;
    381         bool       interp;
     469        bool       interp=false;
    382470
    383471        /*First, recover interp bool: */
    … …  
    438526
    439527} /*}}}*/
    440 IssmDouble TransientInput::InfinityNorm(void){/*{{{*/
    441 
    442         IssmDouble time;
    443         IssmDouble infnorm;
     528IssmDouble TransientInput::Max(void){/*{{{*/
     529
     530        IssmDouble time;
     531        IssmDouble max;
    444532
    445533        /*First, recover current time from parameters: */
    … …  
    450538
    451539        /*Call input function*/
    452         infnorm=input->InfinityNorm();
     540        max=input->Max();
     541
     542        delete input;
     543
     544        return max;
     545}
     546/*}}}*/
     547IssmDouble TransientInput::MaxAbs(void){/*{{{*/
     548
     549        IssmDouble time;
     550        IssmDouble maxabs;
     551
     552        /*First, recover current time from parameters: */
     553        parameters->FindParam(&time,TimeEnum);
     554
     555        /*Retrieve interpolated values for this time step: */
     556        Input* input=GetTimeInput(time);
     557
     558        /*Call input function*/
     559        maxabs=input->MaxAbs();
    453560
    454561        /*Clean-up and return*/
    455562        delete input;
    456         return infnorm;
    457 }
    458 /*}}}*/
    459 IssmDouble TransientInput::Max(void){/*{{{*/
    460 
    461         IssmDouble time;
    462         IssmDouble max;
     563        return maxabs;
     564
     565}
     566/*}}}*/
     567IssmDouble TransientInput::Min(void){/*{{{*/
     568
     569        IssmDouble time;
     570        IssmDouble min;
    463571
    464572        /*First, recover current time from parameters: */
    … …  
    469577
    470578        /*Call input function*/
    471         max=input->Max();
    472 
    473         delete input;
    474 
    475         return max;
    476 }
    477 /*}}}*/
    478 IssmDouble TransientInput::MaxAbs(void){/*{{{*/
    479 
    480         IssmDouble time;
    481         IssmDouble maxabs;
     579        min=input->Min();
     580
     581        /*Clean-up and return*/
     582        delete input;
     583        return min;
     584
     585}
     586/*}}}*/
     587IssmDouble TransientInput::MinAbs(void){/*{{{*/
     588
     589        IssmDouble time;
     590        IssmDouble minabs;
    482591
    483592        /*First, recover current time from parameters: */
    … …  
    488597
    489598        /*Call input function*/
    490         maxabs=input->MaxAbs();
     599        minabs=input->MinAbs();
    491600
    492601        /*Clean-up and return*/
    493602        delete input;
    494         return maxabs;
    495 
    496 }
    497 /*}}}*/
    498 IssmDouble TransientInput::Min(void){/*{{{*/
    499 
    500         IssmDouble time;
    501         IssmDouble min;
    502 
    503         /*First, recover current time from parameters: */
    504         parameters->FindParam(&time,TimeEnum);
    505 
    506    /*Retrieve interpolated values for this time step: */
    507         Input* input=GetTimeInput(time);
    508 
    509         /*Call input function*/
    510         min=input->Min();
    511 
    512         /*Clean-up and return*/
    513         delete input;
    514         return min;
    515 
    516 }
    517 /*}}}*/
    518 IssmDouble TransientInput::MinAbs(void){/*{{{*/
    519 
    520         IssmDouble time;
    521         IssmDouble minabs;
    522 
    523         /*First, recover current time from parameters: */
    524         parameters->FindParam(&time,TimeEnum);
    525 
    526         /*Retrieve interpolated values for this time step: */
    527         Input* input=GetTimeInput(time);
    528 
    529         /*Call input function*/
    530         minabs=input->MinAbs();
    531 
    532         /*Clean-up and return*/
    533         delete input;
    534603        return minabs;
    535604}
    536605/*}}}*/
    537 void TransientInput::SquareMin(IssmDouble* psquaremin,Parameters* parameters){/*{{{*/
    538 
    539         IssmDouble time;
    540 
    541         /*First, recover current time from parameters: */
    542         parameters->FindParam(&time,TimeEnum);
    543 
    544    /*Retrieve interpolated values for this time step: */
    545         Input* input=GetTimeInput(time);
    546 
    547         /*Call input function*/
    548         input->SquareMin(psquaremin,parameters);
    549 
    550         delete input;
    551 
    552 }
    553 /*}}}*/
Note: See TracChangeset for help on using the changeset viewer.