Changeset 21738


Ignore:
Timestamp:
05/22/17 21:11:02 (8 years ago)
Author:
Eric.Larour
Message:

CHG: introducing new way of computing misfit (local==2) in a discretized way.

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

Legend:

Unmodified
Added
Removed

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

    r21049 r21738  
    1919/*}}}*/
    2020IssmDouble OutputDefinitionsResponsex(FemModel* femmodel,int output_enum);
     21void    GetVectorFromInputsx( IssmDouble** pvector,FemModel* femmodel,int name,int type);
     22
    2123
    2224class Misfit: public Object, public Definition{
    … …  
    2527
    2628                int         definitionenum;
    27                 bool        local;     
     29                int         local;     
    2830                int         model_enum;
    2931                char*       name;
    … …  
    4446                        this->weights_enum = UNDEF;
    4547                        this->timeinterpolation=NULL;
    46                         this->local=true;
     48                        this->local=1;
    4749                        this->misfit=0;
    4850                        this->lock=0;
    … …  
    5052                }
    5153                /*}}}*/
    52                 Misfit(char* in_name, int in_definitionenum, int in_model_enum, int in_observation_enum, char* in_timeinterpolation, bool in_local, int in_weights_enum){/*{{{*/
     54                Misfit(char* in_name, int in_definitionenum, int in_model_enum, int in_observation_enum, char* in_timeinterpolation, int in_local, int in_weights_enum){/*{{{*/
    5355
    5456                        this->definitionenum=in_definitionenum;
    … …  
    133135                         femmodel->parameters->FindParam(&dt,TimesteppingTimeStepEnum);
    134136
    135                          if (this->local){ /*local computation: {{{*/
     137                         if (this->local==1){ /*area integration using elements: {{{*/
    136138
    137139                                 int i;
    … …  
    167169                                 /*What we return is the value of misfit / time: */
    168170                                 return misfit/(time-starttime);
     171                         } /*}}}*/
     172                         else if (this->local==2){ /*vertex by vertex computation: {{{*/
     173
     174                                 IssmDouble* model = NULL;
     175                                 IssmDouble* observation= NULL;
     176                                 IssmDouble* weights= NULL;
     177                                 
     178                                 /*Are we transient?:*/
     179                                 if (time==0){
     180                                         IssmDouble misfit_t=0.;
     181                                         
     182                                         /*get global vectors: */
     183                                         GetVectorFromInputsx(&model,femmodel,model_enum,VertexSIdEnum);
     184                                         GetVectorFromInputsx(&observation,femmodel,observation_enum,VertexSIdEnum);
     185                                         GetVectorFromInputsx(&weights,femmodel,weights_enum,VertexSIdEnum);
     186
     187                                         int count=0;
     188                                         for (int i=0;i<femmodel->vertices->NumberOfVertices();i++){
     189                                                 misfit_t += pow(model[i]-observation[i],2)*weights[i];
     190                                                 if (weights[i]!=0)count++;
     191                                         }
     192                                         misfit=sqrt(misfit_t)/count;
     193
     194                                         /*Free ressources:*/
     195                                         xDelete<IssmDouble>(model);
     196                                         xDelete<IssmDouble>(observation);
     197                                         xDelete<IssmDouble>(weights);
     198
     199                                         /*return value: */
     200                                         return misfit;
     201                                 }
     202                                 else{
     203                                         
     204                                         IssmDouble misfit_t=0.;
     205                                         IssmDouble all_misfit_t=0.;
     206
     207                                         /*If we are locked, return time average: */
     208                                         if(this->lock)return misfit/(time-starttime);
     209
     210                                         /*get global vectors: */
     211                                         GetVectorFromInputsx(&model,femmodel,model_enum,VertexSIdEnum);
     212                                         GetVectorFromInputsx(&observation,femmodel,observation_enum,VertexSIdEnum);
     213                                         GetVectorFromInputsx(&weights,femmodel,weights_enum,VertexSIdEnum);
     214
     215                                         int count=0;
     216                                         for (int i=0;i<femmodel->vertices->NumberOfVertices();i++){
     217                                                 misfit_t += pow(model[i]-observation[i],2)*weights[i];
     218                                                 if (weights[i]!=0)count++;
     219                                         }
     220
     221                                         /*Add this time's contribution to curent misfit: */
     222                                         misfit=sqrt(misfit_t)/count;
     223                                         misfit+=dt*misfit_t;
     224
     225                                         /*Do we lock? i.e. are we at final_time? :*/
     226                                         if(time==finaltime)this->lock=1;
     227                                         
     228                                         /*Free ressources:*/
     229                                         xDelete<IssmDouble>(model);
     230                                         xDelete<IssmDouble>(observation);
     231                                         xDelete<IssmDouble>(weights);
     232
     233                                         /*What we return is the value of misfit / time: */
     234                                         return misfit/(time-starttime);
     235                                 }
     236
    169237                         } /*}}}*/
    170238                         else{ /*global computation: {{{ */

  • issm/trunk-jpl/src/c/modules/ModelProcessorx/CreateOutputDefinitions.cpp

    r21690 r21738  
    9696
    9797                                        /*First create a misfit object for that specific string (misfit_model_string_s[j]):*/
    98                                         output_definitions->AddObject(new Misfit(misfit_name_s[j],StringToEnumx(misfit_definitionstring_s[j]),StringToEnumx(misfit_model_string_s[j]),StringToEnumx(misfit_observation_string_s[j]),misfit_timeinterpolation_s[j],(bool)misfit_local_s[j],StringToEnumx(misfit_weights_string_s[j])));
     98                                        output_definitions->AddObject(new Misfit(misfit_name_s[j],StringToEnumx(misfit_definitionstring_s[j]),StringToEnumx(misfit_model_string_s[j]),StringToEnumx(misfit_observation_string_s[j]),misfit_timeinterpolation_s[j],misfit_local_s[j],StringToEnumx(misfit_weights_string_s[j])));
    9999
    100100                                        /*Now, for this particular misfit object, make sure we plug into the elements: the observation, and the weights.*/

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

    r21722 r21738  
    787787        SealevelAbsoluteEnum,
    788788        SealevelEustaticEnum,
     789        SealevelObsEnum,
     790        SealevelWeightsEnum,
    789791        SealevelriseDeltathicknessEnum,
    790792        SealevelriseMaxiterEnum,

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

    r21722 r21738  
    767767                case SealevelAbsoluteEnum : return "SealevelAbsolute";
    768768                case SealevelEustaticEnum : return "SealevelEustatic";
     769                case SealevelObsEnum : return "SealevelObs";
     770                case SealevelWeightsEnum : return "SealevelWeights";
    769771                case SealevelriseDeltathicknessEnum : return "SealevelriseDeltathickness";
    770772                case SealevelriseMaxiterEnum : return "SealevelriseMaxiter";

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

    r21722 r21738  
    785785              else if (strcmp(name,"SealevelAbsolute")==0) return SealevelAbsoluteEnum;
    786786              else if (strcmp(name,"SealevelEustatic")==0) return SealevelEustaticEnum;
     787              else if (strcmp(name,"SealevelObs")==0) return SealevelObsEnum;
     788              else if (strcmp(name,"SealevelWeights")==0) return SealevelWeightsEnum;
    787789              else if (strcmp(name,"SealevelriseDeltathickness")==0) return SealevelriseDeltathicknessEnum;
    788790              else if (strcmp(name,"SealevelriseMaxiter")==0) return SealevelriseMaxiterEnum;
    … …  
    873875              else if (strcmp(name,"SpcStatic")==0) return SpcStaticEnum;
    874876              else if (strcmp(name,"SpcTransient")==0) return SpcTransientEnum;
    875               else if (strcmp(name,"StringArrayParam")==0) return StringArrayParamEnum;
    876               else if (strcmp(name,"StringParam")==0) return StringParamEnum;
    877877         else stage=8;
    878878   }
    879879   if(stage==8){
    880               if (strcmp(name,"Seg")==0) return SegEnum;
     880              if (strcmp(name,"StringArrayParam")==0) return StringArrayParamEnum;
     881              else if (strcmp(name,"StringParam")==0) return StringParamEnum;
     882              else if (strcmp(name,"Seg")==0) return SegEnum;
    881883              else if (strcmp(name,"SegInput")==0) return SegInputEnum;
    882884              else if (strcmp(name,"Tria")==0) return TriaEnum;
    … …  
    996998              else if (strcmp(name,"P2xP4")==0) return P2xP4Enum;
    997999              else if (strcmp(name,"P1P1")==0) return P1P1Enum;
    998               else if (strcmp(name,"P1P1GLS")==0) return P1P1GLSEnum;
    999               else if (strcmp(name,"MINI")==0) return MINIEnum;
    10001000         else stage=9;
    10011001   }
    10021002   if(stage==9){
    1003               if (strcmp(name,"MINIcondensed")==0) return MINIcondensedEnum;
     1003              if (strcmp(name,"P1P1GLS")==0) return P1P1GLSEnum;
     1004              else if (strcmp(name,"MINI")==0) return MINIEnum;
     1005              else if (strcmp(name,"MINIcondensed")==0) return MINIcondensedEnum;
    10041006              else if (strcmp(name,"TaylorHood")==0) return TaylorHoodEnum;
    10051007              else if (strcmp(name,"LATaylorHood")==0) return LATaylorHoodEnum;

  • issm/trunk-jpl/src/m/classes/misfit.m

    r21049 r21738  
    1818        properties (SetAccess=public)
    1919                %misfit
    20                 name                                            = '';
     20                name                                    = '';
    2121                definitionstring                = ''; %string that identifies this output definition uniquely, from 'Outputdefinition[1-10]'
    2222                model_string                    = ''; %string for field that is modeled
    … …  
    2424                observation_string      = ''; %string for observed field.
    2525                timeinterpolation               = '';
    26                 local                                           = 1;
     26                local                                   = 1;
    2727                weights                                 = NaN; %weight coefficients for every vertex
    2828                weights_string                  = ''; %string to identify this particular set of weights
Note: See TracChangeset for help on using the changeset viewer.