Changeset 24790


Ignore:
Timestamp:
05/05/20 06:20:58 (5 years ago)
Author:
bdef
Message:

NEW: adding several averaging methods to pass from substep to main step

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

Legend:

Unmodified
Added
Removed

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

    r24565 r24790  
    336336        int        smbsubstepping;
    337337        int        hydrologysubstepping;
     338        int        smb_averaging;
    338339        IssmDouble dt,scalar,sediment_storing;
    339340        IssmDouble water_head,sediment_transmitivity;
    … …  
    359360        basalelement->FindParam(&isefficientlayer,HydrologydcIsefficientlayerEnum);
    360361        basalelement->FindParam(&smb_model,SmbEnum);
     362        basalelement->FindParam(&smb_averaging,SmbAveragingEnum);
    361363
    362364        Input2* sed_head_input   = basalelement->GetInput2(SedimentHeadSubstepEnum);
    … …  
    384386                }
    385387                else{
    386                         dummy_input = basalelement->GetInput2(SmbRunoffTransientEnum,time-dt,time); _assert_(dummy_input);
     388                        dummy_input = basalelement->GetInput2(SmbRunoffTransientEnum,time-dt,time,smb_averaging); _assert_(dummy_input);
    387389                }
    388390                surface_runoff_input=xDynamicCast<Input2*>(dummy_input); _assert_(surface_runoff_input);

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

    r24716 r24790  
    234234                virtual void       ControlToVectors(Vector<IssmPDouble>* vector_control, Vector<IssmPDouble>* vector_gradient,int control_enum)=0;
    235235                virtual void       CreateDistanceInputFromSegmentlist(IssmDouble* distances,int distanceenum){_error_("not implemented yet");};
    236                 virtual void       CreateInputTimeAverage(int transientinput_enum,int averagedinput_enum,IssmDouble init_time,IssmDouble end_time){_error_("not implemented yet "<<this->ObjectEnum());};
     236                virtual void       CreateInputTimeAverage(int transientinput_enum,int averagedinput_enum,IssmDouble init_time,IssmDouble end_time,int averaging_method){_error_("not implemented yet "<<this->ObjectEnum());};
    237237                virtual void       ElementResponse(IssmDouble* presponse,int response_enum)=0;
    238238                virtual void       ElementSizes(IssmDouble* phx,IssmDouble* phy,IssmDouble* phz)=0;
    … …  
    250250                virtual Input2*    GetInput2(int inputenum)=0;
    251251                virtual Input2*    GetInput2(int inputenum,IssmDouble time)=0;
    252                 virtual Input2*    GetInput2(int inputenum,IssmDouble start_time,IssmDouble end_time)=0;
     252                virtual Input2*    GetInput2(int inputenum,IssmDouble start_time,IssmDouble end_time,int averaging_method)=0;
    253253                virtual void       GetInputValue(IssmDouble* pvalue,Vertex* vertex,int enumtype){_error_("not implemented yet");};
    254254                virtual void       GetInputValue(IssmDouble* pvalue,Node* node,int enumtype){_error_("not implemented yet");};

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

    r24711 r24790  
    10081008}
    10091009/*}}}*/
    1010 void       Penta::CreateInputTimeAverage(int transientinput_enum,int averagedinput_enum,IssmDouble start_time,IssmDouble end_time){/*{{{*/
     1010void       Penta::CreateInputTimeAverage(int transientinput_enum,int averagedinput_enum,IssmDouble start_time,IssmDouble end_time,int averaging_method){/*{{{*/
    10111011
    10121012        _assert_(end_time>start_time);
    … …  
    10171017        IssmDouble dt;
    10181018        int        found,start_offset,end_offset;
    1019         int        averaging_method=0;
    1020 
    10211019
    10221020        /*Get transient input time steps*/

  • issm/trunk-jpl/src/c/classes/Elements/Penta.h

    r24711 r24790  
    6868                void                            CreateDistanceInputFromSegmentlist(IssmDouble* distances,int distanceenum);
    6969                ElementMatrix* CreateBasalMassMatrix(void);
    70                 void           CreateInputTimeAverage(int transientinput_enum,int averagedinput_enum,IssmDouble start_time,IssmDouble end_time);
     70                void           CreateInputTimeAverage(int transientinput_enum,int averagedinput_enum,IssmDouble start_time,IssmDouble end_time,int averaging_method);
    7171                void           ElementResponse(IssmDouble* presponse,int response_enum);
    7272                void           ElementSizes(IssmDouble* hx,IssmDouble* hy,IssmDouble* hz);
    … …  
    8686                Input2*        GetInput2(int enumtype);
    8787                Input2*        GetInput2(int enumtype,IssmDouble time);
    88                 Input2*        GetInput2(int inputenum,IssmDouble start_time,IssmDouble end_time){_error_("not implemented yet!");};
     88                Input2*        GetInput2(int inputenum,IssmDouble start_time,IssmDouble end_time,int averaging_method){_error_("not implemented yet!");};
    8989                void        GetInputListOnVertices(IssmDouble* pvalue,Input2* input,IssmDouble default_value);
    9090                void        GetInputListOnNodes(IssmDouble* pvalue,Input2* input,IssmDouble default_value);

  • issm/trunk-jpl/src/c/classes/Elements/Seg.h

    r24565 r24790  
    6565                Input2*     GetInput2(int enumtype);
    6666                Input2*     GetInput2(int enumtype,IssmDouble time);
    67                 Input2*     GetInput2(int inputenum,IssmDouble start_time,IssmDouble end_time){_error_("not implemented yet!");};
     67                Input2*     GetInput2(int inputenum,IssmDouble start_time,IssmDouble end_time,int averaging_method){_error_("not implemented yet!");};
    6868                void        GetInputListOnVertices(IssmDouble* pvalue,Input2* input,IssmDouble default_value);
    6969                void        GetInputListOnNodes(IssmDouble* pvalue,Input2* input,IssmDouble default_value);

  • issm/trunk-jpl/src/c/classes/Elements/Tetra.h

    r24565 r24790  
    6969                Input2*     GetInput2(int enumtype);
    7070                Input2*     GetInput2(int enumtype,IssmDouble time);
    71                 Input2*     GetInput2(int inputenum,IssmDouble start_time,IssmDouble end_time){_error_("not implemented yet!");};
     71                Input2*     GetInput2(int inputenum,IssmDouble start_time,IssmDouble end_time,int averaging_method){_error_("not implemented yet!");};
    7272                void        GetInputListOnVertices(IssmDouble* pvalue,Input2* input,IssmDouble default_value);
    7373                void        GetInputListOnNodes(IssmDouble* pvalue,Input2* input,IssmDouble default_value);

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

    r24750 r24790  
    19131913        }
    19141914}/*}}}*/
    1915 Input2*    Tria::GetInput2(int inputenum,IssmDouble start_time, IssmDouble end_time){/*{{{*/
     1915Input2*    Tria::GetInput2(int inputenum,IssmDouble start_time, IssmDouble end_time, int averaging_method){/*{{{*/
    19161916
    19171917        /*Get Input from dataset*/
    … …  
    19201920        }
    19211921        else{
    1922                 TriaInput2* input = this->inputs2->GetTriaInput(inputenum,start_time, end_time);
     1922                TriaInput2* input = this->inputs2->GetTriaInput(inputenum,start_time, end_time, int averaging_method);
    19231923                if(!input) return input;
    19241924
    … …  
    21092109        return datasetinput;
    21102110}/*}}}*/
    2111 void       Tria::CreateInputTimeAverage(int transientinput_enum,int averagedinput_enum,IssmDouble start_time,IssmDouble end_time){/*{{{*/
     2111void       Tria::CreateInputTimeAverage(int transientinput_enum,int averagedinput_enum,IssmDouble start_time,IssmDouble end_time, int averaging_method){/*{{{*/
    21122112
    21132113        _assert_(end_time>start_time);
    … …  
    21182118        IssmDouble dt;
    21192119        int        found,start_offset,end_offset;
    2120         int        averaging_method=0;
    2121 
    21222120
    21232121        /*Get transient input time steps*/
    … …  
    22132211                        break;
    22142212                case 2: /*Harmonic mean*/
    2215                         for(int iv=0;iv<NUMVERTICES;iv++) averaged_values[iv] = 1./(averaged_values[iv]/(end_time - start_time));
     2213                        for(int iv=0;iv<NUMVERTICES;iv++) averaged_values[iv] = (end_time - start_time)/(averaged_values[iv];
    22162214                        break;
    22172215                default:

  • issm/trunk-jpl/src/c/classes/Elements/Tria.h

    r24565 r24790  
    178178                void           AddControlInput(int input_enum,Inputs2* inputs2,IoModel* iomodel,IssmDouble* values,IssmDouble* values_min,IssmDouble* values_max, int interpolation_enum,int id);
    179179                void           DatasetInputCreate(IssmDouble* array,int M,int N,int* individual_enums,int num_inputs,Inputs2* inputs2,IoModel* iomodel,int input_enum);
    180                 void           CreateInputTimeAverage(int transientinput_enum,int averagedinput_enum,IssmDouble init_time,IssmDouble end_time);
     180                void           CreateInputTimeAverage(int transientinput_enum,int averagedinput_enum,IssmDouble init_time,IssmDouble end_time,int averaging_method);
    181181                void           GetInputAveragesUpToCurrentTime(int input_enum,IssmDouble** pvalues, IssmDouble** ptimes, int* pnumtimes, IssmDouble currenttime);
    182182                IssmDouble     GetArea(void);
    … …  
    189189                Input2*        GetInput2(int enumtype);
    190190                Input2*        GetInput2(int enumtype,IssmDouble time);
    191                 Input2*        GetInput2(int inputenum,IssmDouble start_time, IssmDouble end_time);
     191                Input2*        GetInput2(int inputenum,IssmDouble start_time,IssmDouble end_time,int averaging_method);
    192192                DatasetInput2* GetDatasetInput2(int inputenum);
    193193                void           GetInputValue(IssmDouble* pvalue,Node* node,int enumtype);

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

    r24656 r24790  
    52755275}
    52765276/*}}}*/
    5277 void FemModel::AverageTransientInputx(int* transientinput_enum,int* averagedinput_enum,IssmDouble init_time,IssmDouble end_time,int numoutputs){ /*{{{*/
     5277void FemModel::AverageTransientInputx(int* transientinput_enum,int* averagedinput_enum,IssmDouble init_time,IssmDouble end_time,int numoutputs,int averaging_method){ /*{{{*/
    52785278
    52795279        for(int i=0;i<numoutputs;i++){
    52805280                for(int j=0;j<this->elements->Size();j++){
    52815281                        Element* element = xDynamicCast<Element*>(elements->GetObjectByOffset(j));
    5282                         element->CreateInputTimeAverage(transientinput_enum[i],averagedinput_enum[i],init_time,end_time);
     5282                        element->CreateInputTimeAverage(transientinput_enum[i],averagedinput_enum[i],init_time,end_time,averaging_method);
    52835283                }
    52845284        }

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

    r24565 r24790  
    176176                void InitTransientInputx(int* transientinput_enum,int numoutputs);
    177177                void StackTransientInputx(int* input_enum,int* transientinput_enum,IssmDouble hydrotime,int numoutputs);
    178                 void AverageTransientInputx(int* transientinput_enum,int* averagedinput_enum,IssmDouble init_time,IssmDouble end_time,int numoutputs);
     178                void AverageTransientInputx(int* transientinput_enum,int* averagedinput_enum,IssmDouble init_time,IssmDouble end_time,int numoutputs,int averaging_method);
    179179                void UpdateConstraintsx(void);
    180180                int  UpdateVertexPositionsx(void);

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

    r24565 r24790  
    344344        }
    345345}/*}}}*/
    346 TriaInput2* Inputs2::GetTriaInput(int enum_in,IssmDouble start_time,IssmDouble end_time){/*{{{*/
     346TriaInput2* Inputs2::GetTriaInput(int enum_in,IssmDouble start_time,IssmDouble end_time,int averaging_method){/*{{{*/
    347347
    348348        /*Get input id*/
    … …  
    354354
    355355        if(input->ObjectEnum()==TransientInput2Enum){
    356                 return xDynamicCast<TransientInput2*>(input)->GetTriaInput(start_time,end_time);
     356                return xDynamicCast<TransientInput2*>(input)->GetTriaInput(start_time,end_time,averaging_method);
    357357        }
    358358        else{

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

    r24565 r24790  
    5454                TriaInput2*      GetTriaInput(int enum_type);
    5555                TriaInput2*      GetTriaInput(int enum_type,IssmDouble time);
    56                 TriaInput2*      GetTriaInput(int enum_in,IssmDouble start_time,IssmDouble end_time);
     56                TriaInput2*      GetTriaInput(int enum_in,IssmDouble start_time,IssmDouble end_time,int averaging_method);
    5757                PentaInput2*     GetPentaInput(int enum_type);
    5858                PentaInput2*     GetPentaInput(int enum_type,IssmDouble time);

  • issm/trunk-jpl/src/c/classes/Inputs2/PentaInput2.cpp

    r24556 r24790  
    155155
    156156        _assert_(this);
    157         _assert_(row>=0); 
    158         _assert_(row<this->M); 
     157        _assert_(row>=0);
     158        _assert_(row<this->M);
    159159        _assert_(this->N==1);
    160160
    … …  
    170170                for(int i=0;i<numindices;i++){
    171171                        int row = indices[i];
    172                         _assert_(row>=0); 
    173                         _assert_(row<this->M); 
     172                        _assert_(row>=0);
     173                        _assert_(row<this->M);
    174174                        this->values[row] = values_in[i];
    175175                }
    … …  
    179179                for(int i=0;i<numindices;i++){
    180180                        int row = indices[i];
    181                         _assert_(row>=0); 
    182                         _assert_(row<this->M); 
     181                        _assert_(row>=0);
     182                        _assert_(row<this->M);
    183183                        this->values[row] = values_in[i];
    184184                }
    … …  
    213213        for(int i=0;i<numindices;i++){
    214214                int row = indices[i];
    215                 _assert_(row>=0); 
    216                 _assert_(row<this->M); 
     215                _assert_(row>=0);
     216                _assert_(row<this->M);
    217217                this->element_values[i] = this->values[row];
    218218        }
    … …  
    389389}
    390390/*}}}*/
     391void PentaInput2::Pow(IssmDouble alpha){/*{{{*/
     392
     393        for(int i=0;i<this->M*this->N;i++) this->values[i] = pow(this->values[i],alpha);
     394        for(int i=0;i<PentaRef::NumberofNodes(this->interpolation);i++) this->element_values[i] = pow(this->element_values[i],alpha);
     395}
     396/*}}}*/
    391397void PentaInput2::AXPY(Input2* xinput,IssmDouble alpha){/*{{{*/
    392398
    … …  
    401407}
    402408/*}}}*/
     409void PentaInput2::PointWiseMult(Input2* xinput){/*{{{*/
     410
     411        /*xinput is of the same type, so cast it: */
     412        if(xinput->ObjectEnum()!=PentaInput2Enum) _error_("Operation not permitted because xinput is of type " << EnumToStringx(xinput->ObjectEnum()));
     413        PentaInput2* xpentainput=xDynamicCast<PentaInput2*>(xinput);
     414        if(xpentainput->GetInterpolation()!=this->interpolation) _error_("Operation not permitted because xinput is of type " << EnumToStringx(xinput->ObjectEnum()));
     415
     416        /* we need to check that the vector sizes are identical*/
     417        if(xpentainput->M!=this->M||xpentainput->N!=this->N) _error_("Operation not permitted because the inputs have different sizes");
     418
     419        /*Carry out the pointwise operation depending on type:*/
     420        for(int i=0;i<this->M*this->N;i++) this->values[i] = xpentainput->values[i] * this->values[i];
     421        for(int i=0;i<PentaRef::NumberofNodes(this->interpolation);i++) this->element_values[i] = xpentainput->element_values[i] * this->element_values[i];
     422}
     423/*}}}*/
    403424
    404425/*Object functions*/

  • issm/trunk-jpl/src/c/classes/Inputs2/PentaInput2.h

    r24335 r24790  
    3737                void GetInputValue(IssmDouble* pvalue,Gauss* gauss);
    3838                void Scale(IssmDouble scalar);
     39                void Pow(IssmDouble scalar);
    3940                void AXPY(Input2* xinput,IssmDouble scalar);
     41                void PointWiseMult(Input2* xinput);
    4042                void Serve(int numindices,int* indices);
    4143                void Serve(int row,int numindices);

  • issm/trunk-jpl/src/c/classes/Inputs2/SegInput2.cpp

    r24335 r24790  
    137137
    138138        _assert_(this);
    139         _assert_(row>=0); 
    140         _assert_(row<this->M); 
     139        _assert_(row>=0);
     140        _assert_(row<this->M);
    141141        _assert_(this->N==1);
    142142
    … …  
    152152                for(int i=0;i<numindices;i++){
    153153                        int row = indices[i];
    154                         _assert_(row>=0); 
    155                         _assert_(row<this->M); 
     154                        _assert_(row>=0);
     155                        _assert_(row<this->M);
    156156                        this->values[row] = values_in[i];
    157157                }
    … …  
    161161                for(int i=0;i<numindices;i++){
    162162                        int row = indices[i];
    163                         _assert_(row>=0); 
    164                         _assert_(row<this->M); 
     163                        _assert_(row>=0);
     164                        _assert_(row<this->M);
    165165                        this->values[row] = values_in[i];
    166166                }
    … …  
    170170                for(int i=0;i<numindices;i++){
    171171                        int row = indices[i];
    172                         _assert_(row>=0); 
    173                         _assert_(row<this->M); 
     172                        _assert_(row>=0);
     173                        _assert_(row<this->M);
    174174                        this->values[row] = values_in[i];
    175175                }
    … …  
    202202        for(int i=0;i<numindices;i++){
    203203                int row = indices[i];
    204                 _assert_(row>=0); 
    205                 _assert_(row<this->M); 
     204                _assert_(row>=0);
     205                _assert_(row<this->M);
    206206                this->element_values[i] = this->values[row];
    207207        }
    … …  
    309309}
    310310/*}}}*/
     311void SegInput2::Pow(IssmDouble alpha){/*{{{*/
     312
     313        for(int i=0;i<this->M*this->N;i++) this->values[i] = pow(this->values[i],alpha);
     314        for(int i=0;i<SegRef::NumberofNodes(this->interpolation);i++) this->element_values[i] = pow(this->element_values[i],alpha);
     315}
     316/*}}}*/
    311317void SegInput2::AXPY(Input2* xinput,IssmDouble alpha){/*{{{*/
    312318
    313319        /*xinput is of the same type, so cast it: */
    314320        if(xinput->ObjectEnum()!=SegInput2Enum) _error_("Operation not permitted because xinput is of type " << EnumToStringx(xinput->ObjectEnum()));
    315         SegInput2* xtriainput=xDynamicCast<SegInput2*>(xinput);
    316         if(xtriainput->GetInterpolation()!=this->interpolation) _error_("Operation not permitted because xinput is of type " << EnumToStringx(xinput->ObjectEnum()));
     321        SegInput2* xseginput=xDynamicCast<SegInput2*>(xinput);
     322        if(xseginput->GetInterpolation()!=this->interpolation) _error_("Operation not permitted because xinput is of type " << EnumToStringx(xinput->ObjectEnum()));
    317323
    318324        /*Carry out the AXPY operation depending on type:*/
    319         for(int i=0;i<this->M*this->N;i++) this->values[i] = alpha*xtriainput->values[i] + this->values[i];
    320         for(int i=0;i<SegRef::NumberofNodes(this->interpolation);i++) this->element_values[i] = alpha*xtriainput->element_values[i] + this->element_values[i];
     325        for(int i=0;i<this->M*this->N;i++) this->values[i] = alpha*xseginput->values[i] + this->values[i];
     326        for(int i=0;i<SegRef::NumberofNodes(this->interpolation);i++) this->element_values[i] = alpha*xseginput->element_values[i] + this->element_values[i];
     327}
     328/*}}}*/
     329void SegInput2::PointWiseMult(Input2* xinput){/*{{{*/
     330
     331        /*xinput is of the same type, so cast it: */
     332        if(xinput->ObjectEnum()!=SegInput2Enum) _error_("Operation not permitted because xinput is of type " << EnumToStringx(xinput->ObjectEnum()));
     333        SegInput2* xseginput=xDynamicCast<SegInput2*>(xinput);
     334        if(xseginput->GetInterpolation()!=this->interpolation) _error_("Operation not permitted because xinput is of type " << EnumToStringx(xinput->ObjectEnum()));
     335
     336        /* we need to check that the vector sizes are identical*/
     337        if(xseginput->M!=this->M||xseginput->N!=this->N) _error_("Operation not permitted because the inputs have different sizes");
     338
     339        /*Carry out the AXPY operation depending on type:*/
     340        for(int i=0;i<this->M*this->N;i++) this->values[i] = xseginput->values[i] * this->values[i];
     341        for(int i=0;i<SegRef::NumberofNodes(this->interpolation);i++) this->element_values[i] = xseginput->element_values[i] * this->element_values[i];
    321342}
    322343/*}}}*/

  • issm/trunk-jpl/src/c/classes/Inputs2/SegInput2.h

    r24335 r24790  
    3535                void GetInputValue(IssmDouble* pvalue,Gauss* gauss);
    3636                void Scale(IssmDouble scalar);
     37                void Pow(IssmDouble scalar);
    3738                void AXPY(Input2* xinput,IssmDouble scalar);
     39                void PointWiseMult(Input2* xinput);
    3840                void Serve(int numindices,int* indices);
    3941                void Serve(int row,int numindices);

  • issm/trunk-jpl/src/c/classes/Inputs2/TransientInput2.cpp

    r24711 r24790  
    290290}
    291291/*}}}*/
    292 TriaInput2* TransientInput2::GetTriaInput(IssmDouble start_time, IssmDouble end_time){/*{{{*/
     292TriaInput2* TransientInput2::GetTriaInput(IssmDouble start_time, IssmDouble end_time, int averaging_method){/*{{{*/
    293293
    294294        /*Set current time input*/
    295         this->SetAverageAsCurrentTimeInput(start_time,end_time);
     295        this->SetAverageAsCurrentTimeInput(start_time,end_time,averaging_method);
    296296        _assert_(this->current_input);
    297297
    … …  
    421421
    422422}/*}}}*/
    423 void TransientInput2::SetAverageAsCurrentTimeInput(IssmDouble start_time,IssmDouble end_time){/*{{{*/
    424 
    425         IssmDouble  dt;
     423void TransientInput2::SetAverageAsCurrentTimeInput(IssmDouble start_time,IssmDouble end_time, int averaging_method){/*{{{*/
     424
     425        IssmDouble  dt,durinv;
    426426        IssmPDouble eps=1.0e-6;
    427427        IssmDouble  dtsum=0;
    428428        int         found,start_offset,end_offset;
    429         int         averaging_method = 0;
    430429
    431430        /*go through the timesteps, and grab offset for start and end*/
    … …  
    463462                                break;
    464463                        case 1: /*Geometric mean*/
    465                                 _error_("Geometric not implemented yet");
     464                                if(offset==start_offset){
     465                                        this->current_input = stepinput->copy();
     466                                        this->current_input->Scale(dt);
     467                                }
     468                                else{
     469                                        this->stepinput->Scale(dt);
     470                                        this->current_input->PointWiseMult(stepinput)
     471                                }
     472                                break;
    466473                        case 2: /*Harmonic mean*/
    467                                 _error_("Harmonic not implemented yet");
     474                                if(offset==start_offset){
     475                                        this->current_input = stepinput->copy();
     476                                        this->current_input->Pow(-1);
     477                                        this->current_input->Scale(dt);
     478                                }
     479                                else{
     480                                        this->stepinput->Pow(-1);
     481                                        this->current_input->AXPY(stepinput,dt);
     482                                }
    468483                        default:
    469484                                _error_("averaging method is not recognised");
    … …  
    472487                offset+=1;
    473488        }
    474                 /*Integration done, now normalize*/
     489        _assert_(dtsum>0);
     490        durinv=1./dtsum;
     491        /*Integration done, now normalize*/
    475492        switch(averaging_method){
    476493                case 0: //Arithmetic mean
    477                         this->current_input->Scale(1/(dtsum));
     494                        this->current_input->Scale(durinv);
    478495                        break;
    479496                case 1: /*Geometric mean*/
    480                         _error_("Geometric not implemented yet");
     497                        this->current_input->Pow(durinv);
     498                        break;
    481499                case 2: /*Harmonic mean*/
    482                         _error_("Harmonic not implemented yet");
     500                        this->current_input->Scale(durinv);
     501                        this->current_input->Pow(-1);
    483502                default:
    484503                        _error_("averaging method is not recognised");

  • issm/trunk-jpl/src/c/classes/Inputs2/TransientInput2.h

    r24565 r24790  
    5050                TriaInput2*  GetTriaInput();
    5151                TriaInput2*  GetTriaInput(IssmDouble time);
    52                 TriaInput2*  GetTriaInput(IssmDouble start_time,IssmDouble end_time);
     52                TriaInput2*  GetTriaInput(IssmDouble start_time,IssmDouble end_time,int averaging_method);
    5353                TriaInput2*  GetTriaInput(int offset);
    5454                PentaInput2* GetPentaInput();
    … …  
    5959                int          GetTimeInputOffset(IssmDouble time);
    6060                void         SetCurrentTimeInput(IssmDouble time);
    61                 void         SetAverageAsCurrentTimeInput(IssmDouble start_time,IssmDouble end_time);
     61                void         SetAverageAsCurrentTimeInput(IssmDouble start_time,IssmDouble end_time,int averaging_method);
    6262                /*numerics:*/
    6363

  • issm/trunk-jpl/src/c/classes/Inputs2/TriaInput2.cpp

    r24360 r24790  
    142142
    143143        _assert_(this);
    144         _assert_(row>=0); 
    145         _assert_(row<this->M); 
     144        _assert_(row>=0);
     145        _assert_(row<this->M);
    146146        _assert_(this->N==1);
    147147
    … …  
    157157                for(int i=0;i<numindices;i++){
    158158                        int row = indices[i];
    159                         _assert_(row>=0); 
    160                         _assert_(row<this->M); 
     159                        _assert_(row>=0);
     160                        _assert_(row<this->M);
    161161                        this->values[row] = values_in[i];
    162162                }
    … …  
    166166                for(int i=0;i<numindices;i++){
    167167                        int row = indices[i];
    168                         _assert_(row>=0); 
    169                         _assert_(row<this->M); 
     168                        _assert_(row>=0);
     169                        _assert_(row<this->M);
    170170                        this->values[row] = values_in[i];
    171171                }
    … …  
    175175                for(int i=0;i<numindices;i++){
    176176                        int row = indices[i];
    177                         _assert_(row>=0); 
    178                         _assert_(row<this->M); 
     177                        _assert_(row>=0);
     178                        _assert_(row<this->M);
    179179                        this->values[row] = values_in[i];
    180180                }
    … …  
    207207        for(int i=0;i<numindices;i++){
    208208                int row = indices[i];
    209                 _assert_(row>=0); 
    210                 _assert_(row<this->M); 
     209                _assert_(row>=0);
     210                _assert_(row<this->M);
    211211                this->element_values[i] = this->values[row];
    212212        }
    … …  
    364364}
    365365/*}}}*/
     366void TriaInput2::Pow(IssmDouble alpha){/*{{{*/
     367
     368        for(int i=0;i<this->M*this->N;i++) this->values[i] = pow(this->values[i],alpha);
     369        for(int i=0;i<TriaRef::NumberofNodes(this->interpolation);i++) this->element_values[i] = pow(this->element_values[i],alpha);
     370}
     371/*}}}*/
    366372void TriaInput2::AXPY(Input2* xinput,IssmDouble alpha){/*{{{*/
    367373
    … …  
    376382}
    377383/*}}}*/
     384void TriaInput2::PointWiseMult(Input2* xinput){/*{{{*/
     385
     386        /*xinput is of the same type, so cast it: */
     387        if(xinput->ObjectEnum()!=TriaInput2Enum) _error_("Operation not permitted because xinput is of type " << EnumToStringx(xinput->ObjectEnum()));
     388        TriaInput2* xtriainput=xDynamicCast<TriaInput2*>(xinput);
     389        if(xtriainput->GetInterpolation()!=this->interpolation) _error_("Operation not permitted because xinput is of type " << EnumToStringx(xinput->ObjectEnum()));
     390
     391        /* we need to check that the vector sizes are identical*/
     392        if(xtriainput->M!=this->M||xtriainput->N!=this->N) _error_("Operation not permitted because the inputs have different sizes");
     393
     394        /*Carry out the AXPY operation depending on type:*/
     395        for(int i=0;i<this->M*this->N;i++) this->values[i] = xtriainput->values[i] * this->values[i];
     396        for(int i=0;i<TriaRef::NumberofNodes(this->interpolation);i++) this->element_values[i] = xtriainput->element_values[i] * this->element_values[i];
     397}
     398/*}}}*/
    378399
    379400/*Object functions*/

  • issm/trunk-jpl/src/c/classes/Inputs2/TriaInput2.h

    r24335 r24790  
    3838                void GetInputValue(IssmDouble* pvalue,Gauss* gauss);
    3939                void Scale(IssmDouble scalar);
     40                void Pow(IssmDouble scalar);
    4041                void AXPY(Input2* xinput,IssmDouble scalar);
     42                void PointWiseMult(Input2* xinput);
    4143                void Serve(int numindices,int* indices);
    4244                void Serve(int row,int numindices);

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

    r24565 r24790  
    6363
    6464                        if(dtslices>1){
    65                                 int        substep, numaveragedinput;
     65                                int        substep, numaveragedinput, hydro_averaging;
    6666                                IssmDouble global_time, subtime, yts;
    6767                                IssmDouble dt, subdt;
    … …  
    7070            femmodel->parameters->FindParam(&dt,TimesteppingTimeStepEnum);
    7171            femmodel->parameters->FindParam(&yts,ConstantsYtsEnum);
     72                                femmodel->parameters->FindParam(&hydro_averaging,HydrologyAveragingEnum);
    7273
    7374                                subtime=global_time-dt; //getting the time back to the start of the timestep
    … …  
    117118                                }
    118119                                /*averaging the stack*/
    119                                 femmodel->AverageTransientInputx(&transientinput[0],&averagedinput[0],global_time-dt,subtime,numaveragedinput);
     120                                femmodel->AverageTransientInputx(&transientinput[0],&averagedinput[0],global_time-dt,subtime,numaveragedinput,hydro_averaging);
    120121
    121122                                /*And reseting to global time*/

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

    r24656 r24790  
    5555        /*if yes compute necessary intermediaries and start looping*/
    5656        if (dtslices>1){
    57                 int        substep;
     57                int        substep,smb_averaging;
    5858                IssmDouble global_time,subtime,yts;
    5959                IssmDouble dt,subdt;
    … …  
    6262                femmodel->parameters->FindParam(&dt,TimesteppingTimeStepEnum);
    6363                femmodel->parameters->FindParam(&yts,ConstantsYtsEnum);
     64                femmodel->parameters->FindParam(&smb_averaging,SmbAveragingEnum);
    6465
    6566                subtime=global_time-dt; //getting the time back to the start of the timestep
    … …  
    8384                delete analysis;
    8485      /*averaging the transient input*/
    85                 femmodel->AverageTransientInputx(&transientinput[0],&averagedinput[0],global_time-dt,subtime,numaveragedinput);
     86                femmodel->AverageTransientInputx(&transientinput[0],&averagedinput[0],global_time-dt,subtime,numaveragedinput,smb_averaging);
    8687                /*and reset timesteping variables to original*/
    8788                femmodel->parameters->SetParam(global_time,TimeEnum);
Note: See TracChangeset for help on using the changeset viewer.