Changeset 16793

Timestamp:

11/15/13 16:18:54 (11 years ago)

Author:

Eric.Larour

Message:

CHG: finished implementing an arbitrary misfit routine, which returns the temporally and spatially averaged misfit
of the model with resepct to observations.

File:

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

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

@@ -13 +13 @@
 #include "./Elements/Elements.h"
 #include "./FemModel.h"
+#include "../modules/SurfaceAreax/SurfaceAreax.h"
+#include "../classes/Params/Parameters.h"
 /*}}}*/
 
@@ -25 +27 @@
                 char*       timeinterpolation;
                 
+                IssmDouble  misfit; //value carried over in time.
+                int         lock; // if lock is on, we just return the value stored in "misfit".  this is used so we don't compute misfit past the final_time
+                
                 /*Misfit constructors, destructors :*/
                 /*FUNCTION Misfit() {{{*/
@@ -34 +39 @@
                         this->weights_enum = UNDEF;
                         this->timeinterpolation=NULL;
+                        this->misfit=0;
+                        this->lock=0;
 
                 }
@@ -49 +56 @@
                         this->observation_enum=in_observation_enum;
                         this->weights_enum=in_weights_enum;
+                        
+                        this->misfit=0;
+                        this->lock=0;
                 }
                 /*}}}*/
@@ -55 +65 @@
                         if(this->name)xDelete(this->name);
                         if(this->timeinterpolation)xDelete(this->timeinterpolation);
+                        this->misfit=0;
+                        this->lock=0;
                 }
                 /*}}}*/
@@ -101 +113 @@
 
                          int i;
-                         IssmDouble response=0;
-                         IssmDouble all_response=0;
+                         IssmDouble misfit_t=0;
+                         IssmDouble all_misfit_t=0;
+                         IssmDouble dt;
+                         IssmDouble area;
+                         IssmDouble time,starttime,finaltime;
+
+                         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);
-                                 response+=element->Misfit(model_enum,observation_enum,weights_enum);
+                                 misfit_t+=element->Misfit(model_enum,observation_enum,weights_enum);
                          }
+
+                         ISSM_MPI_Allreduce ( (void*)&misfit_t,(void*)&all_misfit_t,1,ISSM_MPI_DOUBLE,ISSM_MPI_SUM,IssmComm::GetComm());
+                         misfit_t=all_misfit_t;
                          
-                         ISSM_MPI_Allreduce ( (void*)&response,(void*)&all_response,1,ISSM_MPI_DOUBLE,ISSM_MPI_SUM,IssmComm::GetComm());
-                         return all_response;
+                         /*Divide by surface area: */
+                         SurfaceAreax(&area,femmodel);
+                         misfit_t=misfit_t/area;
+
+                         /*Recover delta_t: */
+                         femmodel->parameters->FindParam(&dt,TimesteppingTimeStepEnum);
+
+                         /*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);
                  }
                         /*}}}*/