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

Last change on this file since 24686 was 24686, checked in by Mathieu Morlighem, 5 years ago

merged trunk-jpl and trunk for revision 24684
File size: 7.4 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/gauss/Gauss.h"
23/*}}}*/
24
25/*Misfit constructors, destructors :*/
26Misfit::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/*}}}*/
40Misfit::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/*}}}*/
59Misfit::~Misfit(){/*{{{*/
60 if(this->name)xDelete(this->name);
61 if(this->timeinterpolation)xDelete(this->timeinterpolation);
62 this->misfit=0;
63 this->lock=0;
64}
65/*}}}*/
66/*Object virtual function resolutoin: */
67Object* 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/*}}}*/
74void Misfit::DeepEcho(void){/*{{{*/
75 this->Echo();
76}
77/*}}}*/
78void 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/*}}}*/
87int Misfit::Id(void){/*{{{*/
88 return -1;
89}
90/*}}}*/
91void Misfit::Marshall(char** pmarshalled_data,int* pmarshalled_data_size, int marshall_direction){/*{{{*/
92 _error_("not implemented yet!");
93}
94/*}}}*/
95int Misfit::ObjectEnum(void){/*{{{*/
96 return MisfitEnum;
97}
98/*}}}*/
99/*Definition virtual function resolutoin: */
100int Misfit::DefinitionEnum(){/*{{{*/
101 return this->definitionenum;
102}
103/*}}}*/
104char* 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/*}}}*/
111IssmDouble Misfit::Response(FemModel* femmodel){/*{{{*/
112
113 /*diverse: */
114 IssmDouble time,starttime,finaltime;
115 IssmDouble dt;
116
117 /*recover time parameters: */
118 femmodel->parameters->FindParam(&starttime,TimesteppingStartTimeEnum);
119 femmodel->parameters->FindParam(&finaltime,TimesteppingFinalTimeEnum);
120 femmodel->parameters->FindParam(&time,TimeEnum);
121 femmodel->parameters->FindParam(&dt,TimesteppingTimeStepEnum);
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++){
135 Element* element=(Element*)femmodel->elements->GetObjectByOffset(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: */
170 GetVectorFromInputsx(&model,&msize,femmodel,model_enum);
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: */
198 GetVectorFromInputsx(&model,&msize,femmodel,model_enum);
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: */
234 model=OutputDefinitionsResponsex(femmodel,this->model_enum);
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 /*}}}*/
Note: See TracBrowser for help on using the repository browser.