source: issm/trunk/src/c/classes/Misfit.h@ 20500

Last change on this file since 20500 was 20500, checked in by Mathieu Morlighem, 9 years ago

merged trunk-jpl and trunk for revision 20497
File size: 6.1 KB
Line 
1/*!\file Misfit.h
2 * \brief: header file for Misfit object
3 */
4
5#ifndef _MISFIT_H_
6#define _MISFIT_H_
7
8/*Headers:*/
9/*{{{*/
10#include "./Definition.h"
11#include "../datastructures/datastructures.h"
12#include "./Elements/Element.h"
13#include "./Elements/Elements.h"
14#include "./FemModel.h"
15#include "../modules/SurfaceAreax/SurfaceAreax.h"
16#include "../classes/Params/Parameters.h"
17#include "../classes/Inputs/Input.h"
18#include "../classes/gauss/Gauss.h"
19/*}}}*/
20IssmDouble OutputDefinitionsResponsex(FemModel* femmodel,int output_enum);
21
22class Misfit: public Object, public Definition{
23
24 public:
25
26 int definitionenum;
27 char* name;
28 int model_enum;
29 int observation_enum;
30 int weights_enum;
31 char* timeinterpolation;
32 bool local;
33
34 IssmDouble misfit; //value carried over in time.
35 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
36
37 /*Misfit constructors, destructors :*/
38 Misfit(){/*{{{*/
39
40 this->definitionenum = -1;
41 this->name = NULL;
42 this->model_enum = UNDEF;
43 this->observation_enum = UNDEF;
44 this->weights_enum = UNDEF;
45 this->timeinterpolation=NULL;
46 this->local=true;
47 this->misfit=0;
48 this->lock=0;
49
50 }
51 /*}}}*/
52 Misfit(char* in_name, int in_definitionenum, int in_model_enum, int in_observation_enum, char* in_timeinterpolation, bool in_local, int in_weights_enum){/*{{{*/
53
54 this->definitionenum=in_definitionenum;
55 this->name = xNew<char>(strlen(in_name)+1);
56 xMemCpy<char>(this->name,in_name,strlen(in_name)+1);
57
58 this->timeinterpolation = xNew<char>(strlen(in_timeinterpolation)+1);
59 xMemCpy<char>(this->timeinterpolation,in_timeinterpolation,strlen(in_timeinterpolation)+1);
60
61 this->model_enum=in_model_enum;
62 this->observation_enum=in_observation_enum;
63 this->weights_enum=in_weights_enum;
64 this->local=in_local;
65
66 this->misfit=0;
67 this->lock=0;
68 }
69 /*}}}*/
70 ~Misfit(){/*{{{*/
71 if(this->name)xDelete(this->name);
72 if(this->timeinterpolation)xDelete(this->timeinterpolation);
73 this->misfit=0;
74 this->lock=0;
75 }
76 /*}}}*/
77 /*Object virtual function resolutoin: */
78 void 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 void DeepEcho(void){/*{{{*/
88 this->Echo();
89 }
90 /*}}}*/
91 int Id(void){/*{{{*/
92 return -1;
93 }
94 /*}}}*/
95 int ObjectEnum(void){/*{{{*/
96 return MisfitEnum;
97 }
98 /*}}}*/
99 Object* copy() {/*{{{*/
100 Misfit* mf = new Misfit(this->name,this->definitionenum, this->model_enum,this->observation_enum,this->timeinterpolation,this->local,this->weights_enum);
101 mf->misfit=this->misfit;
102 mf->lock=this->lock;
103 return (Object*) mf;
104 }
105 /*}}}*/
106 void Marshall(char** pmarshalled_data,int* pmarshalled_data_size, int marshall_direction){/*{{{*/
107 _error_("not implemented yet!");
108 }
109 /*}}}*/
110 /*Definition virtual function resolutoin: */
111 char* Name(){/*{{{*/
112
113 char* name2=xNew<char>(strlen(this->name)+1);
114 xMemCpy(name2,this->name,strlen(this->name)+1);
115
116 return name2;
117 }
118 /*}}}*/
119 int DefinitionEnum(){/*{{{*/
120
121 return this->definitionenum;
122 }
123 /*}}}*/
124 IssmDouble Response(FemModel* femmodel){/*{{{*/
125
126 /*diverse: */
127 IssmDouble time,starttime,finaltime;
128 IssmDouble dt;
129
130 /*recover time parameters: */
131 femmodel->parameters->FindParam(&starttime,TimesteppingStartTimeEnum);
132 femmodel->parameters->FindParam(&finaltime,TimesteppingFinalTimeEnum);
133 femmodel->parameters->FindParam(&time,TimeEnum);
134 femmodel->parameters->FindParam(&dt,TimesteppingTimeStepEnum);
135
136 if (this->local){ /*local computation: {{{*/
137
138 int i;
139 IssmDouble misfit_t=0.;
140 IssmDouble all_misfit_t=0.;
141 IssmDouble area_t=0.;
142 IssmDouble all_area_t;
143
144
145 /*If we are locked, return time average: */
146 if(this->lock)return misfit/(time-starttime);
147
148 for(i=0;i<femmodel->elements->Size();i++){
149 Element* element=(Element*)femmodel->elements->GetObjectByOffset(i);
150 misfit_t+=element->Misfit(model_enum,observation_enum,weights_enum);
151 area_t+=element->MisfitArea(weights_enum);
152 }
153
154 ISSM_MPI_Allreduce ( (void*)&misfit_t,(void*)&all_misfit_t,1,ISSM_MPI_DOUBLE,ISSM_MPI_SUM,IssmComm::GetComm());
155 ISSM_MPI_Allreduce ( (void*)&area_t,(void*)&all_area_t,1,ISSM_MPI_DOUBLE,ISSM_MPI_SUM,IssmComm::GetComm());
156 area_t=all_area_t;
157 misfit_t=all_misfit_t;
158
159 /*Divide by surface area if not nill!: */
160 if (area_t!=0) misfit_t=misfit_t/area_t;
161
162 /*Add this time's contribution to curent misfit: */
163 misfit+=dt*misfit_t;
164
165 /*Do we lock? i.e. are we at final_time? :*/
166 if(time==finaltime)this->lock=1;
167
168 /*What we return is the value of misfit / time: */
169 return misfit/(time-starttime);
170 } /*}}}*/
171 else{ /*global computation: {{{ */
172
173 IssmDouble model, observation;
174
175 /*If we are locked, return time average: */
176 if(this->lock)return misfit/(time-starttime);
177
178 /*First, the global model response: */
179 model=OutputDefinitionsResponsex(femmodel,this->model_enum);
180 /*Now, the observation is buried inside the elements, go fish it in the first element (cludgy, needs fixing): */
181 Element* element=(Element*)femmodel->elements->GetObjectByOffset(0); _assert_(element);
182 Input* input = element->GetInput(observation_enum); _assert_(input);
183 input->GetInputAverage(&observation);
184
185 /*Add this time's contribution to curent misfit: */
186 misfit+=dt*(model-observation);
187
188 /*Do we lock? i.e. are we at final_time? :*/
189 if(time==finaltime)this->lock=1;
190
191 /*What we return is the value of misfit / time: */
192 return misfit/(time-starttime);
193 } /*}}}*/
194
195 }
196 /*}}}*/
197};
198
199#endif /* _MISFIT_H_ */
Note: See TracBrowser for help on using the repository browser.