source: issm/trunk-jpl/src/c/analyses/SmbAnalysis.cpp@ 19554

Last change on this file since 19554 was 19554, checked in by Eric.Larour, 10 years ago

CHG: finished first prototype implementation of the GEMB model inside the Gembx module of the SurfaceMassBalancex module,
and inside the SMBgemb.m class, as well as the SMBGemb method of the Element.cpp class.
Had to introduce the concept also of a DoubleArrayInput to keep track of the vertical gridded datasets within each element.

File size: 12.2 KB
RevLine 
[19528]1#include "./SmbAnalysis.h"
2#include "../toolkits/toolkits.h"
3#include "../classes/classes.h"
4#include "../shared/shared.h"
5#include "../modules/modules.h"
6
7/*Model processing*/
8void SmbAnalysis::CreateConstraints(Constraints* constraints,IoModel* iomodel){/*{{{*/
9 /*No constraints*/
10}/*}}}*/
11void SmbAnalysis::CreateLoads(Loads* loads, IoModel* iomodel){/*{{{*/
12 /*No loads*/
13}/*}}}*/
14void SmbAnalysis::CreateNodes(Nodes* nodes,IoModel* iomodel){/*{{{*/
15 ::CreateNodes(nodes,iomodel,SmbAnalysisEnum,P1Enum);
16}/*}}}*/
17int SmbAnalysis::DofsPerNode(int** doflist,int domaintype,int approximation){/*{{{*/
18 return 1;
19}/*}}}*/
20void SmbAnalysis::UpdateElements(Elements* elements,IoModel* iomodel,int analysis_counter,int analysis_type){/*{{{*/
21
22 int smb_model;
23 bool isdelta18o,ismungsm,isd18opd;
24
25 /*Update elements: */
26 int counter=0;
27 for(int i=0;i<iomodel->numberofelements;i++){
28 if(iomodel->my_elements[i]){
29 Element* element=(Element*)elements->GetObjectByOffset(counter);
30 element->Update(i,iomodel,analysis_counter,analysis_type,P1Enum);
31 counter++;
32 }
33 }
34
[19554]35 /*Figure out smb model: */
36 iomodel->Constant(&smb_model,SmbEnum);
37
38
39 Input* bof=NULL;
40 Element* element=NULL;
[19528]41 switch(smb_model){
42 case SMBforcingEnum:
43 iomodel->FetchDataToInput(elements,SmbMassBalanceEnum,0.);
44 break;
[19554]45 case SMBgembEnum:
46 iomodel->FetchDataToInput(elements,SmbTaEnum);
47 iomodel->FetchDataToInput(elements,SmbVEnum);
48 iomodel->FetchDataToInput(elements,SmbDswrfEnum);
49 iomodel->FetchDataToInput(elements,SmbDlwrfEnum);
50 iomodel->FetchDataToInput(elements,SmbPEnum);
51 iomodel->FetchDataToInput(elements,SmbEAirEnum);
52 iomodel->FetchDataToInput(elements,SmbPAirEnum);
53 iomodel->FetchDataToInput(elements,SmbZTopEnum);
54 iomodel->FetchDataToInput(elements,SmbDzTopEnum);
55 iomodel->FetchDataToInput(elements,SmbDzMinEnum);
56 iomodel->FetchDataToInput(elements,SmbZYEnum);
57 iomodel->FetchDataToInput(elements,SmbZMaxEnum);
58 iomodel->FetchDataToInput(elements,SmbZMinEnum);
59 iomodel->FetchDataToInput(elements,SmbTmeanEnum);
60 iomodel->FetchDataToInput(elements,SmbCEnum);
61 iomodel->FetchDataToInput(elements,SmbTzEnum);
62 iomodel->FetchDataToInput(elements,SmbVzEnum);
63 break;
[19528]64 case SMBpddEnum:
65 iomodel->Constant(&isdelta18o,SmbIsdelta18oEnum);
66 iomodel->Constant(&ismungsm,SmbIsmungsmEnum);
67 iomodel->FetchDataToInput(elements,ThermalSpctemperatureEnum);
68 iomodel->FetchDataToInput(elements,SmbS0pEnum);
69 iomodel->FetchDataToInput(elements,SmbS0tEnum);
70 if(isdelta18o || ismungsm){
71 iomodel->FetchDataToInput(elements,SmbTemperaturesLgmEnum);
72 iomodel->FetchDataToInput(elements,SmbTemperaturesPresentdayEnum);
73 iomodel->FetchDataToInput(elements,SmbPrecipitationsPresentdayEnum);
74 iomodel->FetchDataToInput(elements,SmbPrecipitationsLgmEnum);
75 }
76 else{
77 iomodel->FetchDataToInput(elements,SmbPrecipitationEnum);
78 iomodel->FetchDataToInput(elements,SmbMonthlytemperaturesEnum);
79 }
80 break;
81 case SMBd18opddEnum:
82 iomodel->Constant(&ismungsm,SmbIsmungsmEnum);
83 iomodel->Constant(&isd18opd,SmbIsd18opdEnum);
84 iomodel->FetchDataToInput(elements,ThermalSpctemperatureEnum);
85 iomodel->FetchDataToInput(elements,SmbS0pEnum);
86 iomodel->FetchDataToInput(elements,SmbS0tEnum);
87 if (isd18opd){
88 iomodel->FetchDataToInput(elements,SmbTemperaturesPresentdayEnum);
89 iomodel->FetchDataToInput(elements,SmbPrecipitationsPresentdayEnum);
90 }
91
92 break;
93 case SMBgradientsEnum:
94 iomodel->FetchDataToInput(elements,SmbHrefEnum);
95 iomodel->FetchDataToInput(elements,SmbSmbrefEnum);
96 iomodel->FetchDataToInput(elements,SmbBPosEnum);
97 iomodel->FetchDataToInput(elements,SmbBNegEnum);
98 break;
99 case SMBhenningEnum:
100 iomodel->FetchDataToInput(elements,SmbSmbrefEnum,0.);
101 break;
102 case SMBcomponentsEnum:
103 iomodel->FetchDataToInput(elements,SmbAccumulationEnum,0.);
104 iomodel->FetchDataToInput(elements,SmbEvaporationEnum,0.);
105 iomodel->FetchDataToInput(elements,SmbRunoffEnum,0.);
106 break;
107 case SMBmeltcomponentsEnum:
108 iomodel->FetchDataToInput(elements,SmbAccumulationEnum,0.);
109 iomodel->FetchDataToInput(elements,SmbEvaporationEnum,0.);
110 iomodel->FetchDataToInput(elements,SmbMeltEnum,0.);
111 iomodel->FetchDataToInput(elements,SmbRefreezeEnum,0.);
112 break;
113 default:
114 _error_("Surface mass balance model "<<EnumToStringx(smb_model)<<" not supported yet");
115 }
[19554]116
117
[19528]118
119}/*}}}*/
120void SmbAnalysis::UpdateParameters(Parameters* parameters,IoModel* iomodel,int solution_enum,int analysis_enum){/*{{{*/
121
122 int numoutputs;
123 char** requestedoutputs = NULL;
124 bool isdelta18o,ismungsm,isd18opd,interp;
125 int smb_model;
126 IssmDouble *temp = NULL;
127 int N,M;
128
129 parameters->AddObject(iomodel->CopyConstantObject(SmbEnum));
130
131 iomodel->Constant(&smb_model,SmbEnum);
132 iomodel->Constant(&interp,TimesteppingInterpForcingsEnum);
133
134 switch(smb_model){
135 case SMBforcingEnum:
136 /*Nothing to add to parameters*/
137 break;
[19554]138 case SMBgembEnum:
139 parameters->AddObject(iomodel->CopyConstantObject(SmbAIdxEnum));
140 parameters->AddObject(iomodel->CopyConstantObject(SmbSwIdxEnum));
141 parameters->AddObject(iomodel->CopyConstantObject(SmbDenIdxEnum));
142 parameters->AddObject(iomodel->CopyConstantObject(SmbOutputFreqEnum));
143 parameters->AddObject(iomodel->CopyConstantObject(SmbCldFracEnum));
144 parameters->AddObject(iomodel->CopyConstantObject(SmbT0wetEnum));
145 parameters->AddObject(iomodel->CopyConstantObject(SmbT0dryEnum));
146 parameters->AddObject(iomodel->CopyConstantObject(SmbKEnum));
147 parameters->AddObject(iomodel->CopyConstantObject(SmbASnowEnum));
148 parameters->AddObject(iomodel->CopyConstantObject(SmbAIceEnum));
149 parameters->AddObject(iomodel->CopyConstantObject(SmbDtEnum));
150 parameters->AddObject(iomodel->CopyConstantObject(SmbIsgraingrowthEnum));
151 parameters->AddObject(iomodel->CopyConstantObject(SmbIsalbedoEnum));
152 parameters->AddObject(iomodel->CopyConstantObject(SmbIsshortwaveEnum));
153 parameters->AddObject(iomodel->CopyConstantObject(SmbIsthermalEnum));
154 parameters->AddObject(iomodel->CopyConstantObject(SmbIsaccumulationEnum));
155 parameters->AddObject(iomodel->CopyConstantObject(SmbIsmeltEnum));
156 parameters->AddObject(iomodel->CopyConstantObject(SmbIsdensificationEnum));
157 parameters->AddObject(iomodel->CopyConstantObject(SmbIsturbulentfluxEnum));
158 break;
[19528]159 case SMBpddEnum:
160 parameters->AddObject(iomodel->CopyConstantObject(SmbIsdelta18oEnum));
161 parameters->AddObject(iomodel->CopyConstantObject(SmbIsmungsmEnum));
162 parameters->AddObject(iomodel->CopyConstantObject(SmbDesfacEnum));
163 parameters->AddObject(iomodel->CopyConstantObject(SmbRlapsEnum));
164 parameters->AddObject(iomodel->CopyConstantObject(SmbRlapslgmEnum));
165 iomodel->Constant(&isdelta18o,SmbIsdelta18oEnum);
166 iomodel->Constant(&ismungsm,SmbIsmungsmEnum);
167
168 if(ismungsm){
169 iomodel->FetchData(&temp,&N,&M,SmbPfacEnum); _assert_(N==2);
170 parameters->AddObject(new TransientParam(SmbPfacEnum,&temp[0],&temp[M],interp,M));
171 iomodel->DeleteData(temp,SmbPfacEnum);
172
173 iomodel->FetchData(&temp,&N,&M,SmbTdiffEnum); _assert_(N==2);
174 parameters->AddObject(new TransientParam(SmbTdiffEnum,&temp[0],&temp[M],interp,M));
175 iomodel->DeleteData(temp,SmbTdiffEnum);
176
177 iomodel->FetchData(&temp,&N,&M,SmbSealevEnum); _assert_(N==2);
178 parameters->AddObject(new TransientParam(SmbSealevEnum,&temp[0],&temp[M],interp,M));
179 iomodel->DeleteData(temp,SmbSealevEnum);
180 }
181 if(isdelta18o){
182 iomodel->FetchData(&temp,&N,&M,SmbDelta18oEnum); _assert_(N==2);
183 parameters->AddObject(new TransientParam(SmbDelta18oEnum,&temp[0],&temp[M],interp,M));
184 iomodel->DeleteData(temp,SmbDelta18oEnum);
185
186 iomodel->FetchData(&temp,&N,&M,SmbDelta18oSurfaceEnum); _assert_(N==2);
187 parameters->AddObject(new TransientParam(SmbDelta18oSurfaceEnum,&temp[0],&temp[M],interp,M));
188 iomodel->DeleteData(temp,SmbDelta18oSurfaceEnum);
189 }
190 break;
191 case SMBd18opddEnum:
192 parameters->AddObject(iomodel->CopyConstantObject(SmbIsmungsmEnum));
193 parameters->AddObject(iomodel->CopyConstantObject(SmbIsd18opdEnum));
194 parameters->AddObject(iomodel->CopyConstantObject(SmbDesfacEnum));
195 parameters->AddObject(iomodel->CopyConstantObject(SmbRlapsEnum));
196 parameters->AddObject(iomodel->CopyConstantObject(SmbRlapslgmEnum));
197 iomodel->Constant(&ismungsm,SmbIsmungsmEnum);
198 iomodel->Constant(&isd18opd,SmbIsd18opdEnum);
199 if(isd18opd){
200 iomodel->FetchData(&temp,&N,&M,SmbDelta18oEnum); _assert_(N==2);
201 parameters->AddObject(new TransientParam(SmbDelta18oEnum,&temp[0],&temp[M],interp,M));
202 iomodel->DeleteData(temp,SmbDelta18oEnum);
203
204 parameters->AddObject(iomodel->CopyConstantObject(SmbDpermilEnum));
205 }
206 break;
207 case SMBgradientsEnum:
208 /*Nothing to add to parameters*/
209 break;
210 case SMBhenningEnum:
211 /*Nothing to add to parameters*/
212 break;
213 case SMBcomponentsEnum:
214 /*Nothing to add to parameters*/
215 break;
216 case SMBmeltcomponentsEnum:
217 /*Nothing to add to parameters*/
218 break;
219 default:
220 _error_("Surface mass balance model "<<EnumToStringx(smb_model)<<" not supported yet");
221 }
222
223 iomodel->FetchData(&requestedoutputs,&numoutputs,SmbRequestedOutputsEnum);
224 parameters->AddObject(new IntParam(SmbNumRequestedOutputsEnum,numoutputs));
225 if(numoutputs)parameters->AddObject(new StringArrayParam(SmbRequestedOutputsEnum,requestedoutputs,numoutputs));
226 iomodel->DeleteData(&requestedoutputs,numoutputs,SmbRequestedOutputsEnum);
227
228}/*}}}*/
229
230/*Finite Element Analysis*/
231void SmbAnalysis::Core(FemModel* femmodel){/*{{{*/
232
233 int smb_model;
234
235 /*Figure out smb model: */
236 femmodel->parameters->FindParam(&smb_model,SmbEnum);
237
238 /*branch to correct module*/
239 switch(smb_model){
240 case SMBforcingEnum:
241 /*Nothing to be done*/
242 break;
[19554]243 case SMBgembEnum:
244 Gembx(femmodel);
245 break;
[19528]246 case SMBpddEnum:
247 bool isdelta18o,ismungsm;
248 femmodel->parameters->FindParam(&isdelta18o,SmbIsdelta18oEnum);
249 femmodel->parameters->FindParam(&ismungsm,SmbIsmungsmEnum);
250 if(isdelta18o){
251 if(VerboseSolution()) _printf0_(" call Delta18oParameterization module\n");
252 Delta18oParameterizationx(femmodel);
253 }
254 if(ismungsm){
255 if(VerboseSolution()) _printf0_(" call MungsmtpParameterization module\n");
256 MungsmtpParameterizationx(femmodel);
257 }
258 if(VerboseSolution()) _printf0_(" call positive degree day module\n");
259 PositiveDegreeDayx(femmodel);
260 break;
261 case SMBd18opddEnum:
262 bool isd18opd;
263 femmodel->parameters->FindParam(&isd18opd,SmbIsd18opdEnum);
264 if(isd18opd){
265 if(VerboseSolution()) _printf0_(" call Delta18opdParameterization module\n");
266 Delta18opdParameterizationx(femmodel);
267 if(VerboseSolution()) _printf0_(" call positive degree day module\n");
268 PositiveDegreeDayx(femmodel);
269 }
270 break;
271 case SMBgradientsEnum:
272 if(VerboseSolution())_printf0_(" call smb gradients module\n");
273 SmbGradientsx(femmodel);
274 break;
275 case SMBhenningEnum:
276 if(VerboseSolution())_printf0_(" call smb Henning module\n");
277 SmbHenningx(femmodel);
278 break;
279 case SMBcomponentsEnum:
280 if(VerboseSolution())_printf0_(" call smb Components module\n");
281 SmbComponentsx(femmodel);
282 break;
283 case SMBmeltcomponentsEnum:
284 if(VerboseSolution())_printf0_(" call smb Melt Components module\n");
285 SmbMeltComponentsx(femmodel);
286 break;
287 case SMBgcmEnum:
288 /*Nothing to be done*/
289 break;
290 default:
291 _error_("Surface mass balance model "<<EnumToStringx(smb_model)<<" not supported yet");
292 }
293
294}/*}}}*/
295ElementVector* SmbAnalysis::CreateDVector(Element* element){/*{{{*/
296 _error_("not implemented");
297}/*}}}*/
298ElementMatrix* SmbAnalysis::CreateJacobianMatrix(Element* element){/*{{{*/
299_error_("Not implemented");
300}/*}}}*/
301ElementMatrix* SmbAnalysis::CreateKMatrix(Element* element){/*{{{*/
302 _error_("not implemented yet");
303}/*}}}*/
304ElementVector* SmbAnalysis::CreatePVector(Element* element){/*{{{*/
305_error_("not implemented yet");
306}/*}}}*/
307void SmbAnalysis::GetSolutionFromInputs(Vector<IssmDouble>* solution,Element* element){/*{{{*/
308 _error_("not implemented yet");
309}/*}}}*/
310void SmbAnalysis::GradientJ(Vector<IssmDouble>* gradient,Element* element,int control_type,int control_index){/*{{{*/
311 _error_("Not implemented yet");
312}/*}}}*/
313void SmbAnalysis::InputUpdateFromSolution(IssmDouble* solution,Element* element){/*{{{*/
314 _error_("not implemented yet");
315}/*}}}*/
316void SmbAnalysis::UpdateConstraints(FemModel* femmodel){/*{{{*/
317 /*Default, do nothing*/
318 return;
319}/*}}}*/
Note: See TracBrowser for help on using the repository browser.