Changeset 19063

Timestamp:

02/02/15 14:46:37 (10 years ago)

Author:

Eric.Larour

Message:

CHG: added new capability to compute misfit when we are dealing with a global observation, such as Ice Volume or Mass.

File:

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

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

@@ -15 +15 @@
 #include "../modules/SurfaceAreax/SurfaceAreax.h"
 #include "../classes/Params/Parameters.h"
+#include "../classes/Inputs/Input.h"
+#include "../classes/gauss/Gauss.h"
 /*}}}*/
+IssmDouble OutputDefinitionsResponsex(FemModel* femmodel,int output_enum);
 
 class Misfit: public Object, public Definition{
@@ -27 +30 @@
                 int         weights_enum;
                 char*       timeinterpolation;
+                bool        local;     
                 
                 IssmDouble  misfit; //value carried over in time.
@@ -40 +44 @@
                         this->weights_enum = UNDEF;
                         this->timeinterpolation=NULL;
+                        this->local=true;
                         this->misfit=0;
                         this->lock=0;
@@ -45 +50 @@
                 }
                 /*}}}*/
-                Misfit(char* in_name, int in_definitionenum, int in_model_enum, int in_observation_enum, char* in_timeinterpolation, int in_weights_enum){/*{{{*/
+                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){/*{{{*/
 
                         this->definitionenum=in_definitionenum;
@@ -57 +62 @@
                         this->observation_enum=in_observation_enum;
                         this->weights_enum=in_weights_enum;
+                        this->local=in_local;
                         
                         this->misfit=0;
@@ -76 +82 @@
                         _printf_("    weights_enum: " << weights_enum << " " << EnumToStringx(weights_enum) << "\n");
                         _printf_("    timeinterpolation: " << timeinterpolation << "\n");
+                        _printf_("    local: " << local << "\n");
                 }
                 /*}}}*/
@@ -91 +98 @@
                 /*}}}*/
                 Object* copy() {/*{{{*/
-                        Misfit* mf = new Misfit(this->name,this->definitionenum, this->model_enum,this->observation_enum,this->timeinterpolation,this->weights_enum);
+                        Misfit* mf = new Misfit(this->name,this->definitionenum, this->model_enum,this->observation_enum,this->timeinterpolation,this->local,this->weights_enum);
                         mf->misfit=this->misfit;
                         mf->lock=this->lock;
@@ -112 +119 @@
                 /*}}}*/
                  IssmDouble Response(FemModel* femmodel){/*{{{*/
-
-                         int i;
-                         IssmDouble misfit_t=0.;
-                         IssmDouble all_misfit_t=0.;
+                                 
+                         /*diverse: */
+                         IssmDouble time,starttime,finaltime;
                          IssmDouble dt;
-                         IssmDouble area_t=0.;
-                         IssmDouble all_area_t;
-                         IssmDouble time,starttime,finaltime;
-
+                         
+                         /*recover time parameters: */
                          femmodel->parameters->FindParam(&starttime,TimesteppingStartTimeEnum);
                          femmodel->parameters->FindParam(&finaltime,TimesteppingFinalTimeEnum);
                          femmodel->parameters->FindParam(&time,TimeEnum);
-
-                         /*If we are locked, return time average: */
-                         if(this->lock)return misfit/(time-starttime);
-
-
-                         for(i=0;i<femmodel->elements->Size();i++){
-                                 Element* element=(Element*)femmodel->elements->GetObjectByOffset(i);
-                                 misfit_t+=element->Misfit(model_enum,observation_enum,weights_enum);
-                                 area_t+=element->MisfitArea(weights_enum);
-                         }
-
-                         ISSM_MPI_Allreduce ( (void*)&misfit_t,(void*)&all_misfit_t,1,ISSM_MPI_DOUBLE,ISSM_MPI_SUM,IssmComm::GetComm());
-                         ISSM_MPI_Allreduce ( (void*)&area_t,(void*)&all_area_t,1,ISSM_MPI_DOUBLE,ISSM_MPI_SUM,IssmComm::GetComm());
-                         area_t=all_area_t;
-                         misfit_t=all_misfit_t;
-                         
-                         /*Divide by surface area if not nill!: */
-                         if (area_t!=0) misfit_t=misfit_t/area_t;
-
-                         /*Recover delta_t: */
                          femmodel->parameters->FindParam(&dt,TimesteppingTimeStepEnum);
 
-                         /*Add this time's contribution to curent misfit: */
-                         misfit+=dt*misfit_t;
+                         if (this->local){ /*local computation: {{{*/
 
-                         /*Do we lock? i.e. are we at final_time? :*/
-                         if(time==finaltime)this->lock=1;
+                                 int i;
+                                 IssmDouble misfit_t=0.;
+                                 IssmDouble all_misfit_t=0.;
+                                 IssmDouble area_t=0.;
+                                 IssmDouble all_area_t;
 
-                         /*What we return is the value of misfit / time: */
-                         return misfit/(time-starttime);
+                        
+                                 /*If we are locked, return time average: */
+                                 if(this->lock)return misfit/(time-starttime);
+
+                                 for(i=0;i<femmodel->elements->Size();i++){
+                                         Element* element=(Element*)femmodel->elements->GetObjectByOffset(i);
+                                         misfit_t+=element->Misfit(model_enum,observation_enum,weights_enum);
+                                         area_t+=element->MisfitArea(weights_enum);
+                                 }
+
+                                 ISSM_MPI_Allreduce ( (void*)&misfit_t,(void*)&all_misfit_t,1,ISSM_MPI_DOUBLE,ISSM_MPI_SUM,IssmComm::GetComm());
+                                 ISSM_MPI_Allreduce ( (void*)&area_t,(void*)&all_area_t,1,ISSM_MPI_DOUBLE,ISSM_MPI_SUM,IssmComm::GetComm());
+                                 area_t=all_area_t;
+                                 misfit_t=all_misfit_t;
+                                 
+                                 /*Divide by surface area if not nill!: */
+                                 if (area_t!=0) misfit_t=misfit_t/area_t;
+
+                                 /*Add this time's contribution to curent misfit: */
+                                 misfit+=dt*misfit_t;
+
+                                 /*Do we lock? i.e. are we at final_time? :*/
+                                 if(time==finaltime)this->lock=1;
+
+                                 /*What we return is the value of misfit / time: */
+                                 return misfit/(time-starttime);
+                         } /*}}}*/
+                         else{ /*global computation: {{{ */
+                                 
+                                 IssmDouble model, observation;
+                                 
+                                 /*If we are locked, return time average: */
+                                 if(this->lock)return misfit/(time-starttime);
+
+                                 /*First, the global  model response: */
+                                 model=OutputDefinitionsResponsex(femmodel,this->model_enum);
+                                 /*Now, the observation is buried inside the elements, go fish it in the first element (cludgy, needs fixing): */
+                                 Element* element=(Element*)femmodel->elements->GetObjectByOffset(0); _assert_(element);
+                                 Input* input = element->GetInput(observation_enum); _assert_(input);
+                                 input->GetInputAverage(&observation);
+                                 
+                                 /*Add this time's contribution to curent misfit: */
+                                 misfit+=dt*(model-observation);
+                                 
+                                 /*Do we lock? i.e. are we at final_time? :*/
+                                 if(time==finaltime)this->lock=1;
+                                 
+                                 /*What we return is the value of misfit / time: */
+                                 return misfit/(time-starttime);
+                         } /*}}}*/
+
                  }
                         /*}}}*/