Changeset 21738 for issm/trunk-jpl/src/c/classes/Misfit.h

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.

File:

• issm/trunk-jpl/src/c/classes/Misfit.h (modified) (6 diffs)

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

@@ -19 +19 @@
 /*}}}*/
 IssmDouble OutputDefinitionsResponsex(FemModel* femmodel,int output_enum);
+void    GetVectorFromInputsx( IssmDouble** pvector,FemModel* femmodel,int name,int type);
+
 
 class Misfit: public Object, public Definition{
@@ -25 +27 @@
 
                 int         definitionenum;
-                bool        local;     
+                int         local;     
                 int         model_enum;
                 char*       name;
@@ -44 +46 @@
                         this->weights_enum = UNDEF;
                         this->timeinterpolation=NULL;
-                        this->local=true;
+                        this->local=1;
                         this->misfit=0;
                         this->lock=0;
@@ -50 +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){/*{{{*/
+                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){/*{{{*/
 
                         this->definitionenum=in_definitionenum;
@@ -133 +135 @@
                          femmodel->parameters->FindParam(&dt,TimesteppingTimeStepEnum);
 
-                         if (this->local){ /*local computation: {{{*/
+                         if (this->local==1){ /*area integration using elements: {{{*/
 
                                  int i;
@@ -167 +169 @@
                                  /*What we return is the value of misfit / time: */
                                  return misfit/(time-starttime);
+                         } /*}}}*/
+                         else if (this->local==2){ /*vertex by vertex computation: {{{*/
+
+                                 IssmDouble* model = NULL;
+                                 IssmDouble* observation= NULL;
+                                 IssmDouble* weights= NULL;
+                                 
+                                 /*Are we transient?:*/
+                                 if (time==0){
+                                         IssmDouble misfit_t=0.;
+                                         
+                                         /*get global vectors: */
+                                         GetVectorFromInputsx(&model,femmodel,model_enum,VertexSIdEnum);
+                                         GetVectorFromInputsx(&observation,femmodel,observation_enum,VertexSIdEnum);
+                                         GetVectorFromInputsx(&weights,femmodel,weights_enum,VertexSIdEnum);
+
+                                         int count=0;
+                                         for (int i=0;i<femmodel->vertices->NumberOfVertices();i++){
+                                                 misfit_t += pow(model[i]-observation[i],2)*weights[i];
+                                                 if (weights[i]!=0)count++;
+                                         }
+                                         misfit=sqrt(misfit_t)/count;
+
+                                         /*Free ressources:*/
+                                         xDelete<IssmDouble>(model);
+                                         xDelete<IssmDouble>(observation);
+                                         xDelete<IssmDouble>(weights);
+
+                                         /*return value: */
+                                         return misfit;
+                                 }
+                                 else{
+                                         
+                                         IssmDouble misfit_t=0.;
+                                         IssmDouble all_misfit_t=0.;
+
+                                         /*If we are locked, return time average: */
+                                         if(this->lock)return misfit/(time-starttime);
+
+                                         /*get global vectors: */
+                                         GetVectorFromInputsx(&model,femmodel,model_enum,VertexSIdEnum);
+                                         GetVectorFromInputsx(&observation,femmodel,observation_enum,VertexSIdEnum);
+                                         GetVectorFromInputsx(&weights,femmodel,weights_enum,VertexSIdEnum);
+
+                                         int count=0;
+                                         for (int i=0;i<femmodel->vertices->NumberOfVertices();i++){
+                                                 misfit_t += pow(model[i]-observation[i],2)*weights[i];
+                                                 if (weights[i]!=0)count++;
+                                         }
+
+                                         /*Add this time's contribution to curent misfit: */
+                                         misfit=sqrt(misfit_t)/count;
+                                         misfit+=dt*misfit_t;
+
+                                         /*Do we lock? i.e. are we at final_time? :*/
+                                         if(time==finaltime)this->lock=1;
+                                         
+                                         /*Free ressources:*/
+                                         xDelete<IssmDouble>(model);
+                                         xDelete<IssmDouble>(observation);
+                                         xDelete<IssmDouble>(weights);
+
+                                         /*What we return is the value of misfit / time: */
+                                         return misfit/(time-starttime);
+                                 }
+
                          } /*}}}*/
                          else{ /*global computation: {{{ */