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

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

NEW: marhsall strings instead of enums
File size: 13.7 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: */
[20690]36 iomodel->FindConstant(&smb_model,"md.smb.model");
[19554]37
[19528]38 switch(smb_model){
39 case SMBforcingEnum:
[20690]40 iomodel->FetchDataToInput(elements,"md.smb.mass_balance",SmbMassBalanceEnum,0.);
[19528]41 break;
[19554]42 case SMBgembEnum:
[20690]43 iomodel->FetchDataToInput(elements,"md.smb.Ta",SmbTaEnum);
44 iomodel->FetchDataToInput(elements,"md.smb.V",SmbVEnum);
45 iomodel->FetchDataToInput(elements,"md.smb.dswrf",SmbDswrfEnum);
46 iomodel->FetchDataToInput(elements,"md.smb.dlwrf",SmbDlwrfEnum);
47 iomodel->FetchDataToInput(elements,"md.smb.P",SmbPEnum);
48 iomodel->FetchDataToInput(elements,"md.smb.eAir",SmbEAirEnum);
49 iomodel->FetchDataToInput(elements,"md.smb.pAir",SmbPAirEnum);
50 iomodel->FetchDataToInput(elements,"md.smb.zTop",SmbZTopEnum);
51 iomodel->FetchDataToInput(elements,"md.smb.dzTop",SmbDzTopEnum);
52 iomodel->FetchDataToInput(elements,"md.smb.dzMin",SmbDzMinEnum);
53 iomodel->FetchDataToInput(elements,"md.smb.zY",SmbZYEnum);
54 iomodel->FetchDataToInput(elements,"md.smb.zMax",SmbZMaxEnum);
55 iomodel->FetchDataToInput(elements,"md.smb.zMin",SmbZMinEnum);
56 iomodel->FetchDataToInput(elements,"md.smb.Tmean",SmbTmeanEnum);
57 iomodel->FetchDataToInput(elements,"md.smb.C",SmbCEnum);
58 iomodel->FetchDataToInput(elements,"md.smb.Tz",SmbTzEnum);
59 iomodel->FetchDataToInput(elements,"md.smb.Vz",SmbVzEnum);
[19554]60 break;
[19528]61 case SMBpddEnum:
[20690]62 iomodel->FindConstant(&isdelta18o,"md.smb.isdelta18o");
63 iomodel->FindConstant(&ismungsm,"md.smb.ismungsm");
64 iomodel->FetchDataToInput(elements,"md.thermal.spctemperature",ThermalSpctemperatureEnum);
65 iomodel->FetchDataToInput(elements,"md.smb.s0p",SmbS0pEnum);
66 iomodel->FetchDataToInput(elements,"md.smb.s0t",SmbS0tEnum);
[19528]67 if(isdelta18o || ismungsm){
[20690]68 iomodel->FetchDataToInput(elements,"md.smb.temperatures_lgm",SmbTemperaturesLgmEnum);
69 iomodel->FetchDataToInput(elements,"md.smb.temperatures_presentday",SmbTemperaturesPresentdayEnum);
70 iomodel->FetchDataToInput(elements,"md.smb.precipitations_presentday",SmbPrecipitationsPresentdayEnum);
71 iomodel->FetchDataToInput(elements,"md.smb.precipitations_lgm",SmbPrecipitationsLgmEnum);
[19528]72 }
73 else{
[20690]74 iomodel->FetchDataToInput(elements,"md.smb.precipitation",SmbPrecipitationEnum);
75 iomodel->FetchDataToInput(elements,"md.smb.monthlytemperatures",SmbMonthlytemperaturesEnum);
[19528]76 }
77 break;
78 case SMBd18opddEnum:
[20690]79 iomodel->FindConstant(&ismungsm,"md.smb.ismungsm");
80 iomodel->FindConstant(&isd18opd,"md.smb.isd18opd");
81 iomodel->FetchDataToInput(elements,"md.thermal.spctemperature",ThermalSpctemperatureEnum);
82 iomodel->FetchDataToInput(elements,"md.smb.s0p",SmbS0pEnum);
83 iomodel->FetchDataToInput(elements,"md.smb.s0t",SmbS0tEnum);
84 if(isd18opd){
85 iomodel->FetchDataToInput(elements,"md.smb.temperatures_presentday",SmbTemperaturesPresentdayEnum);
86 iomodel->FetchDataToInput(elements,"md.smb.precipitations_presentday",SmbPrecipitationsPresentdayEnum);
[19528]87 }
88
89 break;
90 case SMBgradientsEnum:
[20690]91 iomodel->FetchDataToInput(elements,"md.smb.href",SmbHrefEnum);
92 iomodel->FetchDataToInput(elements,"md.smb.smbref",SmbSmbrefEnum);
93 iomodel->FetchDataToInput(elements,"md.smb.b_pos",SmbBPosEnum);
94 iomodel->FetchDataToInput(elements,"md.smb.b_neg",SmbBNegEnum);
[19528]95 break;
96 case SMBhenningEnum:
[20690]97 iomodel->FetchDataToInput(elements,"md.smb.smbref",SmbSmbrefEnum,0.);
[19528]98 break;
99 case SMBcomponentsEnum:
[20690]100 iomodel->FetchDataToInput(elements,"md.smb.accumulation",SmbAccumulationEnum,0.);
101 iomodel->FetchDataToInput(elements,"md.smb.evaporation",SmbEvaporationEnum,0.);
102 iomodel->FetchDataToInput(elements,"md.smb.runoff",SmbRunoffEnum,0.);
[19528]103 break;
104 case SMBmeltcomponentsEnum:
[20690]105 iomodel->FetchDataToInput(elements,"md.smb.accumulation",SmbAccumulationEnum,0.);
106 iomodel->FetchDataToInput(elements,"md.smb.evaporation",SmbEvaporationEnum,0.);
107 iomodel->FetchDataToInput(elements,"md.smb.melt",SmbMeltEnum,0.);
108 iomodel->FetchDataToInput(elements,"md.smb.refreeze",SmbRefreezeEnum,0.);
[19528]109 break;
110 default:
111 _error_("Surface mass balance model "<<EnumToStringx(smb_model)<<" not supported yet");
112 }
[19554]113
114
[19528]115
116}/*}}}*/
117void SmbAnalysis::UpdateParameters(Parameters* parameters,IoModel* iomodel,int solution_enum,int analysis_enum){/*{{{*/
118
119 int numoutputs;
120 char** requestedoutputs = NULL;
121 bool isdelta18o,ismungsm,isd18opd,interp;
122 int smb_model;
123 IssmDouble *temp = NULL;
124 int N,M;
125
[20690]126 parameters->AddObject(iomodel->CopyConstantObject("md.smb.model",SmbEnum));
[19528]127
[20690]128 iomodel->FindConstant(&smb_model,"md.smb.model");
129 iomodel->FindConstant(&interp,"md.timestepping.interp_forcings");
[19528]130
131 switch(smb_model){
132 case SMBforcingEnum:
133 /*Nothing to add to parameters*/
134 break;
[19554]135 case SMBgembEnum:
[20690]136 parameters->AddObject(iomodel->CopyConstantObject("md.smb.aIdx",SmbAIdxEnum));
137 parameters->AddObject(iomodel->CopyConstantObject("md.smb.swIdx",SmbSwIdxEnum));
138 parameters->AddObject(iomodel->CopyConstantObject("md.smb.denIdx",SmbDenIdxEnum));
139 parameters->AddObject(iomodel->CopyConstantObject("md.smb.outputFreq",SmbOutputFreqEnum));
140 parameters->AddObject(iomodel->CopyConstantObject("md.smb.cldFrac",SmbCldFracEnum));
141 parameters->AddObject(iomodel->CopyConstantObject("md.smb.t0wet",SmbT0wetEnum));
142 parameters->AddObject(iomodel->CopyConstantObject("md.smb.t0dry",SmbT0dryEnum));
143 parameters->AddObject(iomodel->CopyConstantObject("md.smb.K",SmbKEnum));
144 parameters->AddObject(iomodel->CopyConstantObject("md.smb.aSnow",SmbASnowEnum));
145 parameters->AddObject(iomodel->CopyConstantObject("md.smb.aIce",SmbAIceEnum));
146 parameters->AddObject(iomodel->CopyConstantObject("md.smb.dt",SmbDtEnum));
147 parameters->AddObject(iomodel->CopyConstantObject("md.smb.isgraingrowth",SmbIsgraingrowthEnum));
148 parameters->AddObject(iomodel->CopyConstantObject("md.smb.isalbedo",SmbIsalbedoEnum));
149 parameters->AddObject(iomodel->CopyConstantObject("md.smb.isshortwave",SmbIsshortwaveEnum));
150 parameters->AddObject(iomodel->CopyConstantObject("md.smb.isthermal",SmbIsthermalEnum));
151 parameters->AddObject(iomodel->CopyConstantObject("md.smb.isaccumulation",SmbIsaccumulationEnum));
152 parameters->AddObject(iomodel->CopyConstantObject("md.smb.ismelt",SmbIsmeltEnum));
153 parameters->AddObject(iomodel->CopyConstantObject("md.smb.isdensification",SmbIsdensificationEnum));
154 parameters->AddObject(iomodel->CopyConstantObject("md.smb.isturbulentflux",SmbIsturbulentfluxEnum));
155 parameters->AddObject(iomodel->CopyConstantObject("md.smb.InitDensityScaling",SmbInitDensityScalingEnum));
[19554]156 break;
[19528]157 case SMBpddEnum:
[20690]158 parameters->AddObject(iomodel->CopyConstantObject("md.smb.isdelta18o",SmbIsdelta18oEnum));
159 parameters->AddObject(iomodel->CopyConstantObject("md.smb.ismungsm",SmbIsmungsmEnum));
160 parameters->AddObject(iomodel->CopyConstantObject("md.smb.desfac",SmbDesfacEnum));
161 parameters->AddObject(iomodel->CopyConstantObject("md.smb.rlaps",SmbRlapsEnum));
162 parameters->AddObject(iomodel->CopyConstantObject("md.smb.rlapslgm",SmbRlapslgmEnum));
163 iomodel->FindConstant(&isdelta18o,"md.smb.isdelta18o");
164 iomodel->FindConstant(&ismungsm,"md.smb.ismungsm");
[19528]165
166 if(ismungsm){
[20690]167 iomodel->FetchData(&temp,&N,&M,"md.smb.Pfac"); _assert_(N==2);
[19528]168 parameters->AddObject(new TransientParam(SmbPfacEnum,&temp[0],&temp[M],interp,M));
[20690]169 iomodel->DeleteData(temp,"md.smb.Pfac");
[19528]170
[20690]171 iomodel->FetchData(&temp,&N,&M,"md.smb.Tdiff"); _assert_(N==2);
[19528]172 parameters->AddObject(new TransientParam(SmbTdiffEnum,&temp[0],&temp[M],interp,M));
[20690]173 iomodel->DeleteData(temp,"md.smb.Tdiff");
[19528]174
[20690]175 iomodel->FetchData(&temp,&N,&M,"md.smb.sealev"); _assert_(N==2);
[19528]176 parameters->AddObject(new TransientParam(SmbSealevEnum,&temp[0],&temp[M],interp,M));
[20690]177 iomodel->DeleteData(temp,"md.smb.sealev");
[19528]178 }
179 if(isdelta18o){
[20690]180 iomodel->FetchData(&temp,&N,&M,"md.smb.delta18o"); _assert_(N==2);
[19528]181 parameters->AddObject(new TransientParam(SmbDelta18oEnum,&temp[0],&temp[M],interp,M));
[20690]182 iomodel->DeleteData(temp,"md.smb.delta18o");
[19528]183
[20690]184 iomodel->FetchData(&temp,&N,&M,"md.smb.delta18o_surface"); _assert_(N==2);
[19528]185 parameters->AddObject(new TransientParam(SmbDelta18oSurfaceEnum,&temp[0],&temp[M],interp,M));
[20690]186 iomodel->DeleteData(temp,"md.smb.delta18o_surface");
[19528]187 }
188 break;
189 case SMBd18opddEnum:
[20690]190 parameters->AddObject(iomodel->CopyConstantObject("md.smb.ismungsm",SmbIsmungsmEnum));
191 parameters->AddObject(iomodel->CopyConstantObject("md.smb.isd18opd",SmbIsd18opdEnum));
192 parameters->AddObject(iomodel->CopyConstantObject("md.smb.desfac",SmbDesfacEnum));
193 parameters->AddObject(iomodel->CopyConstantObject("md.smb.rlaps",SmbRlapsEnum));
194 parameters->AddObject(iomodel->CopyConstantObject("md.smb.rlapslgm",SmbRlapslgmEnum));
195 iomodel->FindConstant(&ismungsm,"md.smb.ismungsm");
196 iomodel->FindConstant(&isd18opd,"md.smb.isd18opd");
[19528]197 if(isd18opd){
[20690]198 iomodel->FetchData(&temp,&N,&M,"md.smb.delta18o"); _assert_(N==2);
[19528]199 parameters->AddObject(new TransientParam(SmbDelta18oEnum,&temp[0],&temp[M],interp,M));
[20690]200 iomodel->DeleteData(temp,"md.smb.delta18o");
[19528]201
[20690]202 parameters->AddObject(iomodel->CopyConstantObject("md.smb.dpermil",SmbDpermilEnum));
[19528]203 }
204 break;
205 case SMBgradientsEnum:
206 /*Nothing to add to parameters*/
207 break;
208 case SMBhenningEnum:
209 /*Nothing to add to parameters*/
210 break;
211 case SMBcomponentsEnum:
212 /*Nothing to add to parameters*/
213 break;
214 case SMBmeltcomponentsEnum:
215 /*Nothing to add to parameters*/
216 break;
217 default:
218 _error_("Surface mass balance model "<<EnumToStringx(smb_model)<<" not supported yet");
219 }
220
[20690]221 iomodel->FindConstant(&requestedoutputs,&numoutputs,"md.smb.requested_outputs");
[19528]222 parameters->AddObject(new IntParam(SmbNumRequestedOutputsEnum,numoutputs));
223 if(numoutputs)parameters->AddObject(new StringArrayParam(SmbRequestedOutputsEnum,requestedoutputs,numoutputs));
[20690]224 iomodel->DeleteData(&requestedoutputs,numoutputs,"md.smb.requested_outputs");
[19528]225
226}/*}}}*/
227
228/*Finite Element Analysis*/
229void SmbAnalysis::Core(FemModel* femmodel){/*{{{*/
230
231 int smb_model;
232
233 /*Figure out smb model: */
234 femmodel->parameters->FindParam(&smb_model,SmbEnum);
235
236 /*branch to correct module*/
237 switch(smb_model){
238 case SMBforcingEnum:
239 /*Nothing to be done*/
240 break;
[19554]241 case SMBgembEnum:
242 Gembx(femmodel);
243 break;
[19528]244 case SMBpddEnum:
245 bool isdelta18o,ismungsm;
246 femmodel->parameters->FindParam(&isdelta18o,SmbIsdelta18oEnum);
247 femmodel->parameters->FindParam(&ismungsm,SmbIsmungsmEnum);
248 if(isdelta18o){
249 if(VerboseSolution()) _printf0_(" call Delta18oParameterization module\n");
250 Delta18oParameterizationx(femmodel);
251 }
252 if(ismungsm){
253 if(VerboseSolution()) _printf0_(" call MungsmtpParameterization module\n");
254 MungsmtpParameterizationx(femmodel);
255 }
256 if(VerboseSolution()) _printf0_(" call positive degree day module\n");
257 PositiveDegreeDayx(femmodel);
258 break;
259 case SMBd18opddEnum:
260 bool isd18opd;
261 femmodel->parameters->FindParam(&isd18opd,SmbIsd18opdEnum);
262 if(isd18opd){
263 if(VerboseSolution()) _printf0_(" call Delta18opdParameterization module\n");
264 Delta18opdParameterizationx(femmodel);
265 if(VerboseSolution()) _printf0_(" call positive degree day module\n");
266 PositiveDegreeDayx(femmodel);
267 }
268 break;
269 case SMBgradientsEnum:
270 if(VerboseSolution())_printf0_(" call smb gradients module\n");
271 SmbGradientsx(femmodel);
272 break;
273 case SMBhenningEnum:
274 if(VerboseSolution())_printf0_(" call smb Henning module\n");
275 SmbHenningx(femmodel);
276 break;
277 case SMBcomponentsEnum:
278 if(VerboseSolution())_printf0_(" call smb Components module\n");
279 SmbComponentsx(femmodel);
280 break;
281 case SMBmeltcomponentsEnum:
282 if(VerboseSolution())_printf0_(" call smb Melt Components module\n");
283 SmbMeltComponentsx(femmodel);
284 break;
285 case SMBgcmEnum:
286 /*Nothing to be done*/
287 break;
288 default:
289 _error_("Surface mass balance model "<<EnumToStringx(smb_model)<<" not supported yet");
290 }
291
292}/*}}}*/
293ElementVector* SmbAnalysis::CreateDVector(Element* element){/*{{{*/
294 _error_("not implemented");
295}/*}}}*/
296ElementMatrix* SmbAnalysis::CreateJacobianMatrix(Element* element){/*{{{*/
297_error_("Not implemented");
298}/*}}}*/
299ElementMatrix* SmbAnalysis::CreateKMatrix(Element* element){/*{{{*/
300 _error_("not implemented yet");
301}/*}}}*/
302ElementVector* SmbAnalysis::CreatePVector(Element* element){/*{{{*/
303_error_("not implemented yet");
304}/*}}}*/
305void SmbAnalysis::GetSolutionFromInputs(Vector<IssmDouble>* solution,Element* element){/*{{{*/
306 _error_("not implemented yet");
307}/*}}}*/
308void SmbAnalysis::GradientJ(Vector<IssmDouble>* gradient,Element* element,int control_type,int control_index){/*{{{*/
309 _error_("Not implemented yet");
310}/*}}}*/
311void SmbAnalysis::InputUpdateFromSolution(IssmDouble* solution,Element* element){/*{{{*/
312 _error_("not implemented yet");
313}/*}}}*/
314void SmbAnalysis::UpdateConstraints(FemModel* femmodel){/*{{{*/
315 /*Default, do nothing*/
316 return;
317}/*}}}*/
Note: See TracBrowser for help on using the repository browser.