Ice Sheet System Model  4.18
Code documentation
Misfit.cpp
Go to the documentation of this file.
1 
5 /*Headers:*/
6 /*{{{*/
7 #ifdef HAVE_CONFIG_H
8  #include <config.h>
9 #else
10 #error "Cannot compile with HAVE_CONFIG_H symbol! run configure first!"
11 #endif
12 
13 #include "./classes.h"
16 #include "../datastructures/datastructures.h"
17 #include "./Elements/Element.h"
18 #include "./Elements/Elements.h"
19 #include "./FemModel.h"
20 #include "../modules/SurfaceAreax/SurfaceAreax.h"
21 #include "../classes/Params/Parameters.h"
22 #include "../classes/gauss/Gauss.h"
23 /*}}}*/
24 
25 /*Misfit constructors, destructors :*/
26 Misfit::Misfit(){/*{{{*/
27 
28  this->definitionenum = -1;
29  this->name = NULL;
30  this->model_enum = UNDEF;
31  this->observation_enum = UNDEF;
32  this->weights_enum = UNDEF;
33  this->timeinterpolation=NULL;
34  this->local=1;
35  this->misfit=0;
36  this->lock=0;
37 
38 }
39 /*}}}*/
40 Misfit::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){/*{{{*/
41 
42  this->definitionenum=in_definitionenum;
43 
44  this->name = xNew<char>(strlen(in_name)+1);
45  xMemCpy<char>(this->name,in_name,strlen(in_name)+1);
46 
47  this->timeinterpolation = xNew<char>(strlen(in_timeinterpolation)+1);
48  xMemCpy<char>(this->timeinterpolation,in_timeinterpolation,strlen(in_timeinterpolation)+1);
49 
50  this->model_enum=in_model_enum;
51  this->observation_enum=in_observation_enum;
52  this->weights_enum=in_weights_enum;
53  this->local=in_local;
54 
55  this->misfit=0;
56  this->lock=0;
57 }
58 /*}}}*/
59 Misfit::~Misfit(){/*{{{*/
60  if(this->name)xDelete(this->name);
62  this->misfit=0;
63  this->lock=0;
64 }
65 /*}}}*/
66 /*Object virtual function resolutoin: */
67 Object* Misfit::copy() {/*{{{*/
68  Misfit* mf = new Misfit(this->name,this->definitionenum, this->model_enum,this->observation_enum,this->timeinterpolation,this->local,this->weights_enum);
69  mf->misfit=this->misfit;
70  mf->lock=this->lock;
71  return (Object*) mf;
72 }
73 /*}}}*/
74 void Misfit::DeepEcho(void){/*{{{*/
75  this->Echo();
76 }
77 /*}}}*/
78 void Misfit::Echo(void){/*{{{*/
79  _printf_(" Misfit: " << name << " " << this->definitionenum << "\n");
80  _printf_(" model_enum: " << model_enum << " " << EnumToStringx(model_enum) << "\n");
81  _printf_(" observation_enum: " << observation_enum << " " << EnumToStringx(observation_enum) << "\n");
82  _printf_(" weights_enum: " << weights_enum << " " << EnumToStringx(weights_enum) << "\n");
83  _printf_(" timeinterpolation: " << timeinterpolation << "\n");
84  _printf_(" local: " << local << "\n");
85 }
86 /*}}}*/
87 int Misfit::Id(void){/*{{{*/
88  return -1;
89 }
90 /*}}}*/
91 void Misfit::Marshall(char** pmarshalled_data,int* pmarshalled_data_size, int marshall_direction){/*{{{*/
92  _error_("not implemented yet!");
93 }
94 /*}}}*/
95 int Misfit::ObjectEnum(void){/*{{{*/
96  return MisfitEnum;
97 }
98 /*}}}*/
99 /*Definition virtual function resolutoin: */
101  return this->definitionenum;
102 }
103 /*}}}*/
104 char* Misfit::Name(){/*{{{*/
105  char* name2=xNew<char>(strlen(this->name)+1);
106  xMemCpy(name2,this->name,strlen(this->name)+1);
107 
108  return name2;
109 }
110 /*}}}*/
112 
113  /*diverse: */
114  IssmDouble time,starttime,finaltime;
115  IssmDouble dt;
116 
117  /*recover time parameters: */
122 
123  if (this->local==1){ /*area integration using elements: {{{*/
124 
125  int i;
126  IssmDouble misfit_t=0.;
127  IssmDouble all_misfit_t=0.;
128  IssmDouble area_t=0.;
129  IssmDouble all_area_t;
130 
131  /*If we are locked, return time average: */
132  if(this->lock)return misfit/(time-starttime);
133 
134  for(i=0;i<femmodel->elements->Size();i++){
136  misfit_t+=element->Misfit(model_enum,observation_enum,weights_enum);
137  area_t+=element->MisfitArea(weights_enum);
138  }
139 
140  ISSM_MPI_Allreduce ( (void*)&misfit_t,(void*)&all_misfit_t,1,ISSM_MPI_DOUBLE,ISSM_MPI_SUM,IssmComm::GetComm());
141  ISSM_MPI_Allreduce ( (void*)&area_t,(void*)&all_area_t,1,ISSM_MPI_DOUBLE,ISSM_MPI_SUM,IssmComm::GetComm());
142  area_t=all_area_t;
143  misfit_t=all_misfit_t;
144 
145  /*Divide by surface area if not nill!: */
146  if (area_t!=0) misfit_t=misfit_t/area_t;
147 
148  /*Add this time's contribution to curent misfit: */
149  misfit+=dt*misfit_t;
150 
151  /*Do we lock? i.e. are we at final_time? :*/
152  if(time==finaltime)this->lock=1;
153 
154  /*What we return is the value of misfit / time if transient*/
155  if(time!=0.) return misfit/(time-starttime);
156  return misfit;
157  } /*}}}*/
158  else if (this->local==2){ /*vertex by vertex computation: {{{*/
159 
160  IssmDouble* model = NULL;
161  IssmDouble* observation= NULL;
162  IssmDouble* weights= NULL;
163  int msize,osize,wsize;
164 
165  /*Are we transient?:*/
166  if (time==0){
167  IssmDouble misfit_t=0.;
168 
169  /*get global vectors: */
171  GetVectorFromInputsx(&observation,&osize,femmodel,observation_enum);_assert_(msize==osize);
172  GetVectorFromInputsx(&weights,&wsize,femmodel,weights_enum); _assert_(wsize==msize);
173 
174  int count=0;
175  for (int i=0;i<msize;i++){
176  misfit_t += pow(model[i]-observation[i],2)*weights[i];
177  if (weights[i]!=0)count++;
178  }
179  misfit=sqrt(misfit_t/count);
180 
181  /*Free ressources:*/
182  xDelete<IssmDouble>(model);
183  xDelete<IssmDouble>(observation);
184  xDelete<IssmDouble>(weights);
185 
186  /*return value: */
187  return misfit;
188  }
189  else{
190 
191  IssmDouble misfit_t=0.;
192  IssmDouble all_misfit_t=0.;
193 
194  /*If we are locked, return time average: */
195  if(this->lock)return misfit/(time-starttime);
196 
197  /*get global vectors: */
199  GetVectorFromInputsx(&observation,&osize,femmodel,observation_enum);_assert_(msize==osize);
200  GetVectorFromInputsx(&weights,&wsize,femmodel,weights_enum); _assert_(wsize==msize);
201 
202  int count=0;
203  for (int i=0;i<msize;i++){
204  misfit_t += pow(model[i]-observation[i],2)*weights[i];
205  if (weights[i]!=0)count++;
206  }
207  misfit=sqrt(misfit_t/count);
208 
209  /*Add this time's contribution to curent misfit: */
210  misfit=sqrt(misfit_t)/count;
211  misfit+=dt*misfit_t;
212 
213  /*Do we lock? i.e. are we at final_time? :*/
214  if(time==finaltime)this->lock=1;
215 
216  /*Free ressources:*/
217  xDelete<IssmDouble>(model);
218  xDelete<IssmDouble>(observation);
219  xDelete<IssmDouble>(weights);
220 
221  /*What we return is the value of misfit / time: */
222  return misfit/(time-starttime);
223  }
224 
225  } /*}}}*/
226  else{ /*global computation: {{{ */
227 
228  IssmDouble model, observation;
229 
230  /*If we are locked, return time average: */
231  if(this->lock)return misfit/(time-starttime);
232 
233  /*First, the global model response: */
235  /*Now, the observation is buried inside the elements, go fish it in the first element (cludgy, needs fixing): */
236  Element* element = (Element*)femmodel->elements->GetObjectByOffset(0); _assert_(element);
237  Input2* input = element->GetInput2(observation_enum); _assert_(input);
238  input->GetInputAverage(&observation);
239 
240  /*Add this time's contribution to curent misfit: */
241  misfit+=dt*(model-observation);
242 
243  /*Do we lock? i.e. are we at final_time? :*/
244  if(time==finaltime)this->lock=1;
245 
246  /*What we return is the value of misfit / time: */
247  return misfit/(time-starttime);
248  } /*}}}*/
249 
250  }
251  /*}}}*/
DataSet::Size
int Size()
Definition: DataSet.cpp:399
TimesteppingFinalTimeEnum
@ TimesteppingFinalTimeEnum
Definition: EnumDefinitions.h:430
Misfit::Id
int Id(void)
Definition: Misfit.cpp:87
_assert_
#define _assert_(ignore)
Definition: exceptions.h:37
IssmDouble
double IssmDouble
Definition: types.h:37
Misfit::name
char * name
Definition: Misfit.h:22
MisfitEnum
@ MisfitEnum
Definition: EnumDefinitions.h:656
Misfit::lock
int lock
Definition: Misfit.h:27
ISSM_MPI_Allreduce
int ISSM_MPI_Allreduce(void *sendbuf, void *recvbuf, int count, ISSM_MPI_Datatype datatype, ISSM_MPI_Op op, ISSM_MPI_Comm comm)
Definition: issmmpi.cpp:116
_printf_
#define _printf_(StreamArgs)
Definition: Print.h:22
Misfit::Name
char * Name()
Definition: Misfit.cpp:104
ISSM_MPI_SUM
#define ISSM_MPI_SUM
Definition: issmmpi.h:134
FemModel::parameters
Parameters * parameters
Definition: FemModel.h:46
TimeEnum
@ TimeEnum
Definition: EnumDefinitions.h:427
Misfit::DefinitionEnum
int DefinitionEnum()
Definition: Misfit.cpp:100
IssmComm::GetComm
static ISSM_MPI_Comm GetComm(void)
Definition: IssmComm.cpp:30
TimesteppingTimeStepEnum
@ TimesteppingTimeStepEnum
Definition: EnumDefinitions.h:433
GetVectorFromInputsx
void GetVectorFromInputsx(IssmDouble **pvector, int *pvector_size, FemModel *femmodel, int name)
Definition: GetVectorFromInputsx.cpp:81
Misfit::observation_enum
int observation_enum
Definition: Misfit.h:23
Misfit::model_enum
int model_enum
Definition: Misfit.h:21
Misfit::Echo
void Echo(void)
Definition: Misfit.cpp:78
Misfit::timeinterpolation
char * timeinterpolation
Definition: Misfit.h:24
Misfit::weights_enum
int weights_enum
Definition: Misfit.h:25
Element::GetInput2
virtual Input2 * GetInput2(int inputenum)=0
Element
Definition: Element.h:41
Elements.h
ISSM_MPI_DOUBLE
#define ISSM_MPI_DOUBLE
Definition: issmmpi.h:125
Element.h
abstract class for Element object This class is a place holder for the Tria and the Penta elements....
Object
Definition: Object.h:13
Misfit::misfit
IssmDouble misfit
Definition: Misfit.h:28
xDelete
void xDelete(T **&aT_pp)
Definition: MemOps.h:97
EnumToStringx
const char * EnumToStringx(int enum_in)
Definition: EnumToStringx.cpp:15
ExternalResult.h
abstract class for ExternalResult object
UNDEF
#define UNDEF
Definition: constants.h:8
Element::Misfit
virtual IssmDouble Misfit(int modelenum, int observationenum, int weightsenum)=0
Misfit::Misfit
Misfit()
Definition: Misfit.cpp:26
Misfit::Marshall
void Marshall(char **pmarshalled_data, int *pmarshalled_data_size, int marshall_direction)
Definition: Misfit.cpp:91
FemModel::elements
Elements * elements
Definition: FemModel.h:44
TimesteppingStartTimeEnum
@ TimesteppingStartTimeEnum
Definition: EnumDefinitions.h:432
Input2
Definition: Input2.h:18
FemModel
Definition: FemModel.h:31
Misfit::ObjectEnum
int ObjectEnum(void)
Definition: Misfit.cpp:95
Misfit::~Misfit
~Misfit()
Definition: Misfit.cpp:59
Misfit::definitionenum
int definitionenum
Definition: Misfit.h:19
_error_
#define _error_(StreamArgs)
Definition: exceptions.h:49
DataSet::GetObjectByOffset
Object * GetObjectByOffset(int offset)
Definition: DataSet.cpp:334
FemModel.h
OutputDefinitionsResponsex
IssmDouble OutputDefinitionsResponsex(FemModel *femmodel, int output_enum)
Definition: OutputDefinitionsResponsex.cpp:38
Parameters::FindParam
void FindParam(bool *pinteger, int enum_type)
Definition: Parameters.cpp:262
Misfit::local
int local
Definition: Misfit.h:20
xMemCpy
T * xMemCpy(T *dest, const T *src, unsigned int size)
Definition: MemOps.h:152
Element::MisfitArea
virtual IssmDouble MisfitArea(int weightsenum)=0
Misfit::copy
Object * copy()
Definition: Misfit.cpp:67
Results.h
Input2::GetInputAverage
virtual void GetInputAverage(IssmDouble *pvalue)
Definition: Input2.h:33
Misfit::DeepEcho
void DeepEcho(void)
Definition: Misfit.cpp:74
Misfit::Response
IssmDouble Response(FemModel *femmodel)
Definition: Misfit.cpp:111
classes.h
Misfit
Definition: Misfit.h:15
femmodel
FemModel * femmodel
Definition: esmfbinders.cpp:16