source: issm/trunk/src/c/classes/Misfit.cpp@ 22758

Last change on this file since 22758 was 22758, checked in by Mathieu Morlighem, 7 years ago

merged trunk-jpl and trunk for revision 22757
File size: 7.5 KB
Line 
1/*!\file Misfit.cpp
2 * \brief: Misfit object
3 */
4
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"
14#include "./ExternalResults/ExternalResult.h"
15#include "./ExternalResults/Results.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/Inputs/Input.h"
23#include "../classes/gauss/Gauss.h"
24/*}}}*/
25
26/*Misfit constructors, destructors :*/
27Misfit::Misfit(){/*{{{*/
28
29 this->definitionenum = -1;
30 this->name = NULL;
31 this->model_enum = UNDEF;
32 this->observation_enum = UNDEF;
33 this->weights_enum = UNDEF;
34 this->timeinterpolation=NULL;
35 this->local=1;
36 this->misfit=0;
37 this->lock=0;
38
39}
40/*}}}*/
41Misfit::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){/*{{{*/
42
43 this->definitionenum=in_definitionenum;
44
45 this->name = xNew<char>(strlen(in_name)+1);
46 xMemCpy<char>(this->name,in_name,strlen(in_name)+1);
47
48 this->timeinterpolation = xNew<char>(strlen(in_timeinterpolation)+1);
49 xMemCpy<char>(this->timeinterpolation,in_timeinterpolation,strlen(in_timeinterpolation)+1);
50
51 this->model_enum=in_model_enum;
52 this->observation_enum=in_observation_enum;
53 this->weights_enum=in_weights_enum;
54 this->local=in_local;
55
56 this->misfit=0;
57 this->lock=0;
58}
59/*}}}*/
60Misfit::~Misfit(){/*{{{*/
61 if(this->name)xDelete(this->name);
62 if(this->timeinterpolation)xDelete(this->timeinterpolation);
63 this->misfit=0;
64 this->lock=0;
65}
66/*}}}*/
67/*Object virtual function resolutoin: */
68Object* Misfit::copy() {/*{{{*/
69 Misfit* mf = new Misfit(this->name,this->definitionenum, this->model_enum,this->observation_enum,this->timeinterpolation,this->local,this->weights_enum);
70 mf->misfit=this->misfit;
71 mf->lock=this->lock;
72 return (Object*) mf;
73}
74/*}}}*/
75void Misfit::DeepEcho(void){/*{{{*/
76 this->Echo();
77}
78/*}}}*/
79void Misfit::Echo(void){/*{{{*/
80 _printf_(" Misfit: " << name << " " << this->definitionenum << "\n");
81 _printf_(" model_enum: " << model_enum << " " << EnumToStringx(model_enum) << "\n");
82 _printf_(" observation_enum: " << observation_enum << " " << EnumToStringx(observation_enum) << "\n");
83 _printf_(" weights_enum: " << weights_enum << " " << EnumToStringx(weights_enum) << "\n");
84 _printf_(" timeinterpolation: " << timeinterpolation << "\n");
85 _printf_(" local: " << local << "\n");
86}
87/*}}}*/
88int Misfit::Id(void){/*{{{*/
89 return -1;
90}
91/*}}}*/
92void Misfit::Marshall(char** pmarshalled_data,int* pmarshalled_data_size, int marshall_direction){/*{{{*/
93 _error_("not implemented yet!");
94}
95/*}}}*/
96int Misfit::ObjectEnum(void){/*{{{*/
97 return MisfitEnum;
98}
99/*}}}*/
100/*Definition virtual function resolutoin: */
101int Misfit::DefinitionEnum(){/*{{{*/
102 return this->definitionenum;
103}
104/*}}}*/
105char* Misfit::Name(){/*{{{*/
106 char* name2=xNew<char>(strlen(this->name)+1);
107 xMemCpy(name2,this->name,strlen(this->name)+1);
108
109 return name2;
110}
111/*}}}*/
112IssmDouble Misfit::Response(FemModel* femmodel){/*{{{*/
113
114 /*diverse: */
115 IssmDouble time,starttime,finaltime;
116 IssmDouble dt;
117
118 /*recover time parameters: */
119 femmodel->parameters->FindParam(&starttime,TimesteppingStartTimeEnum);
120 femmodel->parameters->FindParam(&finaltime,TimesteppingFinalTimeEnum);
121 femmodel->parameters->FindParam(&time,TimeEnum);
122 femmodel->parameters->FindParam(&dt,TimesteppingTimeStepEnum);
123
124 if (this->local==1){ /*area integration using elements: {{{*/
125
126 int i;
127 IssmDouble misfit_t=0.;
128 IssmDouble all_misfit_t=0.;
129 IssmDouble area_t=0.;
130 IssmDouble all_area_t;
131
132
133 /*If we are locked, return time average: */
134 if(this->lock)return misfit/(time-starttime);
135
136 for(i=0;i<femmodel->elements->Size();i++){
137 Element* element=(Element*)femmodel->elements->GetObjectByOffset(i);
138 misfit_t+=element->Misfit(model_enum,observation_enum,weights_enum);
139 area_t+=element->MisfitArea(weights_enum);
140 }
141
142 ISSM_MPI_Allreduce ( (void*)&misfit_t,(void*)&all_misfit_t,1,ISSM_MPI_DOUBLE,ISSM_MPI_SUM,IssmComm::GetComm());
143 ISSM_MPI_Allreduce ( (void*)&area_t,(void*)&all_area_t,1,ISSM_MPI_DOUBLE,ISSM_MPI_SUM,IssmComm::GetComm());
144 area_t=all_area_t;
145 misfit_t=all_misfit_t;
146
147 /*Divide by surface area if not nill!: */
148 if (area_t!=0) misfit_t=misfit_t/area_t;
149
150 /*Add this time's contribution to curent misfit: */
151 misfit+=dt*misfit_t;
152
153 /*Do we lock? i.e. are we at final_time? :*/
154 if(time==finaltime)this->lock=1;
155
156 /*What we return is the value of misfit / time if transient*/
157 if(time!=0.) return misfit/(time-starttime);
158 return misfit;
159 } /*}}}*/
160 else if (this->local==2){ /*vertex by vertex computation: {{{*/
161
162 IssmDouble* model = NULL;
163 IssmDouble* observation= NULL;
164 IssmDouble* weights= NULL;
165 int msize,osize,wsize;
166
167 /*Are we transient?:*/
168 if (time==0){
169 IssmDouble misfit_t=0.;
170
171 /*get global vectors: */
172 GetVectorFromInputsx(&model,&msize,femmodel,model_enum);
173 GetVectorFromInputsx(&observation,&osize,femmodel,observation_enum);_assert_(msize==osize);
174 GetVectorFromInputsx(&weights,&wsize,femmodel,weights_enum); _assert_(wsize==msize);
175
176 int count=0;
177 for (int i=0;i<msize;i++){
178 misfit_t += pow(model[i]-observation[i],2)*weights[i];
179 if (weights[i]!=0)count++;
180 }
181 misfit=sqrt(misfit_t/count);
182
183 /*Free ressources:*/
184 xDelete<IssmDouble>(model);
185 xDelete<IssmDouble>(observation);
186 xDelete<IssmDouble>(weights);
187
188 /*return value: */
189 return misfit;
190 }
191 else{
192
193 IssmDouble misfit_t=0.;
194 IssmDouble all_misfit_t=0.;
195
196 /*If we are locked, return time average: */
197 if(this->lock)return misfit/(time-starttime);
198
199 /*get global vectors: */
200 GetVectorFromInputsx(&model,&msize,femmodel,model_enum);
201 GetVectorFromInputsx(&observation,&osize,femmodel,observation_enum);_assert_(msize==osize);
202 GetVectorFromInputsx(&weights,&wsize,femmodel,weights_enum); _assert_(wsize==msize);
203
204 int count=0;
205 for (int i=0;i<msize;i++){
206 misfit_t += pow(model[i]-observation[i],2)*weights[i];
207 if (weights[i]!=0)count++;
208 }
209 misfit=sqrt(misfit_t/count);
210
211 /*Add this time's contribution to curent misfit: */
212 misfit=sqrt(misfit_t)/count;
213 misfit+=dt*misfit_t;
214
215 /*Do we lock? i.e. are we at final_time? :*/
216 if(time==finaltime)this->lock=1;
217
218 /*Free ressources:*/
219 xDelete<IssmDouble>(model);
220 xDelete<IssmDouble>(observation);
221 xDelete<IssmDouble>(weights);
222
223 /*What we return is the value of misfit / time: */
224 return misfit/(time-starttime);
225 }
226
227 } /*}}}*/
228 else{ /*global computation: {{{ */
229
230 IssmDouble model, observation;
231
232 /*If we are locked, return time average: */
233 if(this->lock)return misfit/(time-starttime);
234
235 /*First, the global model response: */
236 model=OutputDefinitionsResponsex(femmodel,this->model_enum);
237 /*Now, the observation is buried inside the elements, go fish it in the first element (cludgy, needs fixing): */
238 Element* element=(Element*)femmodel->elements->GetObjectByOffset(0); _assert_(element);
239 Input* input = element->GetInput(observation_enum); _assert_(input);
240 input->GetInputAverage(&observation);
241
242 /*Add this time's contribution to curent misfit: */
243 misfit+=dt*(model-observation);
244
245 /*Do we lock? i.e. are we at final_time? :*/
246 if(time==finaltime)this->lock=1;
247
248 /*What we return is the value of misfit / time: */
249 return misfit/(time-starttime);
250 } /*}}}*/
251
252 }
253 /*}}}*/
Note: See TracBrowser for help on using the repository browser.