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

Last change on this file was 27232, checked in by Mathieu Morlighem, 3 years ago

merged trunk-jpl and trunk for revision 27230
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(MarshallHandle* marshallhandle){/*{{{*/
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 IssmDouble misfit_t=0.;
126 IssmDouble all_misfit_t=0.;
127 IssmDouble area_t=0.;
128 IssmDouble all_area_t;
129
130 /*If we are locked, return time average: */
131 if(this->lock)return misfit/(time-starttime);
132
133 for(Object* & object : femmodel->elements->objects){
134 Element* element=xDynamicCast<Element*>(object);
135 misfit_t+=element->Misfit(model_enum,observation_enum,weights_enum);
136 area_t+=element->MisfitArea(weights_enum);
137 }
138
139 ISSM_MPI_Allreduce ( (void*)&misfit_t,(void*)&all_misfit_t,1,ISSM_MPI_DOUBLE,ISSM_MPI_SUM,IssmComm::GetComm());
140 ISSM_MPI_Allreduce ( (void*)&area_t,(void*)&all_area_t,1,ISSM_MPI_DOUBLE,ISSM_MPI_SUM,IssmComm::GetComm());
141 area_t=all_area_t;
142 misfit_t=all_misfit_t;
143
144 /*Divide by surface area if not nill!: */
145 if (area_t!=0) misfit_t=misfit_t/area_t;
146
147 /*Add this time's contribution to curent misfit: */
148 misfit+=dt*misfit_t;
149
150 /*Do we lock? i.e. are we at final_time? :*/
151 if(time==finaltime)this->lock=1;
152
153 /*What we return is the value of misfit / time if transient*/
154 if(time!=0.) return misfit/(time-starttime);
155 return misfit;
156 } /*}}}*/
157 else if (this->local==2){ /*vertex by vertex computation: {{{*/
158
159 IssmDouble* model = NULL;
160 IssmDouble* observation= NULL;
161 IssmDouble* weights= NULL;
162 int msize,osize,wsize;
163
164 /*Are we transient?:*/
165 if (time==0){
166 IssmDouble misfit_t=0.;
167
168 /*get global vectors: */
169 GetVectorFromInputsx(&model,&msize,femmodel,model_enum);
170 GetVectorFromInputsx(&observation,&osize,femmodel,observation_enum);_assert_(msize==osize);
171 GetVectorFromInputsx(&weights,&wsize,femmodel,weights_enum); _assert_(wsize==msize);
172
173 int count=0;
174 for (int i=0;i<msize;i++){
175 misfit_t += pow(model[i]-observation[i],2)*weights[i];
176 if (weights[i]!=0)count++;
177 }
178 misfit=sqrt(misfit_t/count);
179
180 /*Free resources:*/
181 xDelete<IssmDouble>(model);
182 xDelete<IssmDouble>(observation);
183 xDelete<IssmDouble>(weights);
184
185 /*return value: */
186 return misfit;
187 }
188 else{
189
190 IssmDouble misfit_t=0.;
191 IssmDouble all_misfit_t=0.;
192
193 /*If we are locked, return time average: */
194 if(this->lock)return misfit/(time-starttime);
195
196 /*get global vectors: */
197 GetVectorFromInputsx(&model,&msize,femmodel,model_enum);
198 GetVectorFromInputsx(&observation,&osize,femmodel,observation_enum);_assert_(msize==osize);
199 GetVectorFromInputsx(&weights,&wsize,femmodel,weights_enum); _assert_(wsize==msize);
200
201 int count=0;
202 for (int i=0;i<msize;i++){
203 misfit_t += pow(model[i]-observation[i],2)*weights[i];
204 if (weights[i]!=0)count++;
205 }
206 misfit=sqrt(misfit_t/count);
207
208 /*Add this time's contribution to curent misfit: */
209 misfit=sqrt(misfit_t)/count;
210 misfit+=dt*misfit_t;
211
212 /*Do we lock? i.e. are we at final_time? :*/
213 if(time==finaltime)this->lock=1;
214
215 /*Free resources:*/
216 xDelete<IssmDouble>(model);
217 xDelete<IssmDouble>(observation);
218 xDelete<IssmDouble>(weights);
219
220 /*What we return is the value of misfit / time: */
221 return misfit/(time-starttime);
222 }
223
224 } /*}}}*/
225 else{ /*global computation: {{{ */
226
227 IssmDouble model, observation;
228
229 /*If we are locked, return time average: */
230 if(this->lock) return misfit/(time-starttime);
231
232 /*First, the global model response: */
233 model=OutputDefinitionsResponsex(femmodel,this->model_enum);
234 /*Now, the observation is buried inside the elements, go fish it in the first element (cludgy, needs fixing): */
235 Element* element = (Element*)femmodel->elements->GetObjectByOffset(0); _assert_(element);
236 Input* input = element->GetInput(observation_enum); _assert_(input);
237 input->GetInputAverage(&observation);
238
239 /*Add this time's contribution to curent misfit: */
240 misfit+=dt*(model-observation);
241
242 /*Do we lock? i.e. are we at final_time? :*/
243 if(time==finaltime)this->lock=1;
244
245 /*What we return is the value of misfit / time: */
246 return misfit/(time-starttime);
247 } /*}}}*/
248
249 }
250 /*}}}*/
Note: See TracBrowser for help on using the repository browser.