source: issm/trunk-jpl/src/c/analyses/EnthalpyAnalysis.cpp@ 16837

Last change on this file since 16837 was 16837, checked in by jbondzio, 11 years ago

CHG: added error message when water fraction exceeds 100%
File size: 15.9 KB
Line 
1#include "./EnthalpyAnalysis.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*/
8int EnthalpyAnalysis::DofsPerNode(int** doflist,int meshtype,int approximation){/*{{{*/
9 return 1;
10}/*}}}*/
11void EnthalpyAnalysis::UpdateParameters(Parameters* parameters,IoModel* iomodel,int solution_enum,int analysis_enum){/*{{{*/
12
13 int numoutputs;
14 char** requestedoutputs = NULL;
15
16 parameters->AddObject(iomodel->CopyConstantObject(ThermalStabilizationEnum));
17 parameters->AddObject(iomodel->CopyConstantObject(ThermalIsenthalpyEnum));
18 parameters->AddObject(iomodel->CopyConstantObject(ThermalIsdynamicbasalspcEnum));
19
20 iomodel->FetchData(&requestedoutputs,&numoutputs,ThermalRequestedOutputsEnum);
21 parameters->AddObject(new IntParam(ThermalNumRequestedOutputsEnum,numoutputs));
22 if(numoutputs)parameters->AddObject(new StringArrayParam(ThermalRequestedOutputsEnum,requestedoutputs,numoutputs));
23 iomodel->DeleteData(&requestedoutputs,numoutputs,ThermalRequestedOutputsEnum);
24}/*}}}*/
25void EnthalpyAnalysis::UpdateElements(Elements* elements,IoModel* iomodel,int analysis_counter,int analysis_type){/*{{{*/
26
27 bool dakota_analysis;
28 bool isenthalpy;
29
30 /*Now, is the model 3d? otherwise, do nothing: */
31 if(iomodel->meshtype==Mesh2DhorizontalEnum)return;
32
33 /*Is enthalpy requested?*/
34 iomodel->Constant(&isenthalpy,ThermalIsenthalpyEnum);
35 if(!isenthalpy) return;
36
37 /*Fetch data needed: */
38 iomodel->FetchData(3,TemperatureEnum,WaterfractionEnum,PressureEnum);
39
40 /*Update elements: */
41 int counter=0;
42 for(int i=0;i<iomodel->numberofelements;i++){
43 if(iomodel->my_elements[i]){
44 Element* element=(Element*)elements->GetObjectByOffset(counter);
45 element->Update(i,iomodel,analysis_counter,analysis_type,P1Enum);
46 counter++;
47 }
48 }
49
50 iomodel->Constant(&dakota_analysis,QmuIsdakotaEnum);
51
52 iomodel->FetchDataToInput(elements,ThicknessEnum);
53 iomodel->FetchDataToInput(elements,SurfaceEnum);
54 iomodel->FetchDataToInput(elements,BedEnum);
55 iomodel->FetchDataToInput(elements,FrictionCoefficientEnum);
56 iomodel->FetchDataToInput(elements,FrictionPEnum);
57 iomodel->FetchDataToInput(elements,FrictionQEnum);
58 iomodel->FetchDataToInput(elements,MaskIceLevelsetEnum);
59 iomodel->FetchDataToInput(elements,MaskGroundediceLevelsetEnum);
60 iomodel->FetchDataToInput(elements,MeshElementonbedEnum);
61 iomodel->FetchDataToInput(elements,MeshElementonsurfaceEnum);
62 iomodel->FetchDataToInput(elements,FlowequationElementEquationEnum);
63 iomodel->FetchDataToInput(elements,MaterialsRheologyBEnum);
64 iomodel->FetchDataToInput(elements,MaterialsRheologyNEnum);
65 iomodel->FetchDataToInput(elements,PressureEnum);
66 iomodel->FetchDataToInput(elements,TemperatureEnum);
67 iomodel->FetchDataToInput(elements,WaterfractionEnum);
68 iomodel->FetchDataToInput(elements,EnthalpyEnum);
69 iomodel->FetchDataToInput(elements,BasalforcingsGeothermalfluxEnum);
70 iomodel->FetchDataToInput(elements,WatercolumnEnum);
71 iomodel->FetchDataToInput(elements,BasalforcingsMeltingRateEnum);
72 iomodel->FetchDataToInput(elements,VxEnum);
73 iomodel->FetchDataToInput(elements,VyEnum);
74 iomodel->FetchDataToInput(elements,VzEnum);
75 InputUpdateFromConstantx(elements,0.,VxMeshEnum);
76 InputUpdateFromConstantx(elements,0.,VyMeshEnum);
77 InputUpdateFromConstantx(elements,0.,VzMeshEnum);
78 if(dakota_analysis){
79 elements->InputDuplicate(TemperatureEnum,QmuTemperatureEnum);
80 elements->InputDuplicate(BasalforcingsMeltingRateEnum,QmuMeltingEnum);
81 elements->InputDuplicate(VxMeshEnum,QmuVxMeshEnum);
82 elements->InputDuplicate(VxMeshEnum,QmuVyMeshEnum);
83 elements->InputDuplicate(VxMeshEnum,QmuVzMeshEnum);
84 }
85
86 /*Free data: */
87 iomodel->DeleteData(3,TemperatureEnum,WaterfractionEnum,PressureEnum);
88}/*}}}*/
89void EnthalpyAnalysis::CreateNodes(Nodes* nodes,IoModel* iomodel){/*{{{*/
90
91 if(iomodel->meshtype==Mesh3DEnum) iomodel->FetchData(2,MeshVertexonbedEnum,MeshVertexonsurfaceEnum);
92 ::CreateNodes(nodes,iomodel,EnthalpyAnalysisEnum,P1Enum);
93 iomodel->DeleteData(2,MeshVertexonbedEnum,MeshVertexonsurfaceEnum);
94}/*}}}*/
95void EnthalpyAnalysis::CreateConstraints(Constraints* constraints,IoModel* iomodel){/*{{{*/
96
97 /*Intermediary*/
98 int count;
99 int M,N;
100 bool spcpresent = false;
101 IssmDouble heatcapacity;
102 IssmDouble referencetemperature;
103
104 /*Output*/
105 IssmDouble *spcvector = NULL;
106 IssmDouble* times=NULL;
107 IssmDouble* values=NULL;
108
109 /*Fetch parameters: */
110 iomodel->Constant(&heatcapacity,MaterialsHeatcapacityEnum);
111 iomodel->Constant(&referencetemperature,ConstantsReferencetemperatureEnum);
112
113 /*return if 2d mesh*/
114 if(iomodel->meshtype==Mesh2DhorizontalEnum) return;
115
116 /*Fetch data: */
117 iomodel->FetchData(&spcvector,&M,&N,ThermalSpctemperatureEnum);
118
119 //FIX ME: SHOULD USE IOMODELCREATECONSTRAINTS
120 /*Transient or static?:*/
121 if(M==iomodel->numberofvertices){
122 /*static: just create Constraints objects*/
123 count=0;
124
125 for(int i=0;i<iomodel->numberofvertices;i++){
126 /*keep only this partition's nodes:*/
127 if((iomodel->my_vertices[i])){
128
129 if (!xIsNan<IssmDouble>(spcvector[i])){
130
131 constraints->AddObject(new SpcStatic(iomodel->constraintcounter+count+1,iomodel->nodecounter+i+1,1,heatcapacity*(spcvector[i]-referencetemperature),EnthalpyAnalysisEnum));
132 count++;
133
134 }
135 }
136 }
137 }
138 else if (M==(iomodel->numberofvertices+1)){
139 /*transient: create transient SpcTransient objects. Same logic, except we need to retrieve
140 * various times and values to initialize an SpcTransient object: */
141 count=0;
142
143 /*figure out times: */
144 times=xNew<IssmDouble>(N);
145 for(int j=0;j<N;j++){
146 times[j]=spcvector[(M-1)*N+j];
147 }
148
149 /*Create constraints from x,y,z: */
150 for(int i=0;i<iomodel->numberofvertices;i++){
151
152 /*keep only this partition's nodes:*/
153 if((iomodel->my_vertices[i])){
154
155 /*figure out times and values: */
156 values=xNew<IssmDouble>(N);
157 spcpresent=false;
158 for(int j=0;j<N;j++){
159 values[j]=heatcapacity*(spcvector[i*N+j]-referencetemperature);
160 if(!xIsNan<IssmDouble>(values[j]))spcpresent=true; //NaN means no spc by default
161 }
162
163 if(spcpresent){
164 constraints->AddObject(new SpcTransient(iomodel->constraintcounter+count+1,iomodel->nodecounter+i+1,1,N,times,values,EnthalpyAnalysisEnum));
165 count++;
166 }
167 xDelete<IssmDouble>(values);
168 }
169 }
170 }
171 else{
172 _error_("Size of field " << EnumToStringx(ThermalSpctemperatureEnum) << " not supported");
173 }
174
175 /*Free ressources:*/
176 iomodel->DeleteData(spcvector,ThermalSpctemperatureEnum);
177 xDelete<IssmDouble>(times);
178 xDelete<IssmDouble>(values);
179}/*}}}*/
180void EnthalpyAnalysis::CreateLoads(Loads* loads, IoModel* iomodel){/*{{{*/
181
182 /*No loads */
183}/*}}}*/
184
185/*Finite Element Analysis*/
186ElementMatrix* EnthalpyAnalysis::CreateKMatrix(Element* element){/*{{{*/
187 _error_("not implemented yet");
188}/*}}}*/
189ElementVector* EnthalpyAnalysis::CreatePVector(Element* element){/*{{{*/
190
191 /*compute all load vectors for this element*/
192 ElementVector* pe1=CreatePVectorVolume(element);
193 ElementVector* pe2=CreatePVectorSheet(element);
194 ElementVector* pe3=CreatePVectorShelf(element);
195 ElementVector* pe =new ElementVector(pe1,pe2,pe3);
196
197 /*clean-up and return*/
198 delete pe1;
199 delete pe2;
200 delete pe3;
201 return pe;
202}/*}}}*/
203ElementVector* EnthalpyAnalysis::CreatePVectorVolume(Element* element){/*{{{*/
204
205 /*Intermediaries*/
206 int stabilization;
207 IssmDouble Jdet,phi,dt;
208 IssmDouble enthalpy;
209 IssmDouble kappa,tau_parameter,diameter;
210 IssmDouble u,v,w;
211 IssmDouble scalar_def,scalar_transient;
212 IssmDouble* xyz_list = NULL;
213
214 /*Fetch number of nodes and dof for this finite element*/
215 int numnodes = element->GetNumberOfNodes();
216 int numvertices = element->GetNumberOfVertices();
217
218 /*Initialize Element vector*/
219 ElementVector* pe = element->NewElementVector();
220 IssmDouble* basis = xNew<IssmDouble>(numnodes);
221 IssmDouble* dbasis = xNew<IssmDouble>(3*numnodes);
222 IssmDouble* pressure = xNew<IssmDouble>(numvertices);
223 IssmDouble* enthalpypicard = xNew<IssmDouble>(numvertices);
224
225 /*Retrieve all inputs and parameters*/
226 element->GetVerticesCoordinates(&xyz_list);
227 IssmDouble rho_ice = element->GetMaterialParameter(MaterialsRhoIceEnum);
228 IssmDouble thermalconductivity = element->GetMaterialParameter(MaterialsThermalconductivityEnum);
229 element->FindParam(&dt,TimesteppingTimeStepEnum);
230 element->FindParam(&stabilization,ThermalStabilizationEnum);
231 Input* vx_input=element->GetInput(VxEnum); _assert_(vx_input);
232 Input* vy_input=element->GetInput(VyEnum); _assert_(vy_input);
233 Input* vz_input=element->GetInput(VzEnum); _assert_(vz_input);
234 Input* enthalpy_input = NULL;
235 if(reCast<bool,IssmDouble>(dt)){enthalpy_input = element->GetInput(EnthalpyEnum); _assert_(enthalpy_input);}
236 if(stabilization==2){
237 element->GetInputListOnVertices(enthalpypicard,EnthalpyPicardEnum);
238 element->GetInputListOnVertices(pressure,PressureEnum);
239 diameter=element->MinEdgeLength(xyz_list);
240 }
241
242 /* Start looping on the number of gaussian points: */
243 Gauss* gauss=element->NewGauss(3);
244 for(int ig=gauss->begin();ig<gauss->end();ig++){
245 gauss->GaussPoint(ig);
246
247 element->JacobianDeterminant(&Jdet,xyz_list,gauss);
248 element->NodalFunctions(basis,gauss);
249 element->ViscousHeating(&phi,xyz_list,gauss,vx_input,vy_input,vz_input);
250
251 scalar_def=phi/rho_ice*Jdet*gauss->weight;
252 if(reCast<bool,IssmDouble>(dt)) scalar_def=scalar_def*dt;
253
254 /*TODO: add -beta*laplace T_m(p)*/
255 for(int i=0;i<numnodes;i++) pe->values[i]+=scalar_def*basis[i];
256
257 /* Build transient now */
258 if(reCast<bool,IssmDouble>(dt)){
259 enthalpy_input->GetInputValue(&enthalpy, gauss);
260 scalar_transient=enthalpy*Jdet*gauss->weight;
261 for(int i=0;i<numnodes;i++) pe->values[i]+=scalar_transient*basis[i];
262 }
263
264 if(stabilization==2){
265 element->NodalFunctionsDerivatives(dbasis,xyz_list,gauss);
266
267 vx_input->GetInputValue(&u,gauss);
268 vy_input->GetInputValue(&v,gauss);
269 vz_input->GetInputValue(&w,gauss);
270 kappa = element->EnthalpyDiffusionParameterVolume(numvertices,enthalpypicard,pressure) / rho_ice;
271 tau_parameter = element->StabilizationParameter(u,v,w,diameter,kappa);
272
273 for(int i=0;i<numnodes;i++) pe->values[i]+=tau_parameter*scalar_def*(u*dbasis[0*3+i]+v*dbasis[1*3+i]+w*dbasis[2*3+i]);
274 if(reCast<bool,IssmDouble>(dt)){
275 for(int i=0;i<numnodes;i++) pe->values[i]+=tau_parameter*scalar_transient*(u*dbasis[0*3+i]+v*dbasis[1*3+i]+w*dbasis[2*3+i]);
276 }
277 }
278 }
279
280 /*Clean up and return*/
281 xDelete<IssmDouble>(basis);
282 xDelete<IssmDouble>(dbasis);
283 xDelete<IssmDouble>(pressure);
284 xDelete<IssmDouble>(enthalpypicard);
285 xDelete<IssmDouble>(xyz_list);
286 delete gauss;
287 return pe;
288
289}/*}}}*/
290ElementVector* EnthalpyAnalysis::CreatePVectorSheet(Element* element){/*{{{*/
291 return NULL;
292}/*}}}*/
293ElementVector* EnthalpyAnalysis::CreatePVectorShelf(Element* element){/*{{{*/
294
295 IssmDouble h_pmp,dt,Jdet,scalar_ocean,pressure;
296 IssmDouble *xyz_list_base = NULL;
297
298 /*Get basal element*/
299 if(!element->IsOnBed() || !element->IsFloating()) return NULL;
300
301 /*Fetch number of nodes for this finite element*/
302 int numnodes = element->GetNumberOfNodes();
303
304 /*Initialize vectors*/
305 ElementVector* pe = element->NewElementVector();
306 IssmDouble* basis = xNew<IssmDouble>(numnodes);
307
308 /*Retrieve all inputs and parameters*/
309 element->GetVerticesCoordinatesBase(&xyz_list_base);
310 element->FindParam(&dt,TimesteppingTimeStepEnum);
311 Input* pressure_input=element->GetInput(PressureEnum); _assert_(pressure_input);
312 IssmDouble gravity = element->GetMaterialParameter(ConstantsGEnum);
313 IssmDouble rho_water = element->GetMaterialParameter(MaterialsRhoWaterEnum);
314 IssmDouble rho_ice = element->GetMaterialParameter(MaterialsRhoIceEnum);
315 IssmDouble heatcapacity = element->GetMaterialParameter(MaterialsHeatcapacityEnum);
316 IssmDouble mixed_layer_capacity= element->GetMaterialParameter(MaterialsMixedLayerCapacityEnum);
317 IssmDouble thermal_exchange_vel= element->GetMaterialParameter(MaterialsThermalExchangeVelocityEnum);
318
319 /* Start looping on the number of gaussian points: */
320 Gauss* gauss=element->NewGaussBase(2);
321 for(int ig=gauss->begin();ig<gauss->end();ig++){
322 gauss->GaussPoint(ig);
323
324 element->JacobianDeterminantBase(&Jdet,xyz_list_base,gauss);
325 element->NodalFunctions(basis,gauss);
326
327 pressure_input->GetInputValue(&pressure,gauss);
328 h_pmp=element->PureIceEnthalpy(pressure);
329
330 scalar_ocean=gauss->weight*Jdet*rho_water*mixed_layer_capacity*thermal_exchange_vel*h_pmp/(heatcapacity*rho_ice);
331 if(reCast<bool,IssmDouble>(dt)) scalar_ocean=dt*scalar_ocean;
332
333 for(int i=0;i<numnodes;i++) pe->values[i]+=scalar_ocean*basis[i];
334 }
335
336 /*Clean up and return*/
337 delete gauss;
338 xDelete<IssmDouble>(basis);
339 xDelete<IssmDouble>(xyz_list_base);
340 return pe;
341}/*}}}*/
342void EnthalpyAnalysis::GetSolutionFromInputs(Vector<IssmDouble>* solution,Element* element){/*{{{*/
343 element->GetSolutionFromInputsOneDof(solution,EnthalpyEnum);
344}/*}}}*/
345void EnthalpyAnalysis::InputUpdateFromSolution(IssmDouble* solution,Element* element){/*{{{*/
346
347 bool converged;
348 int i,rheology_law;
349 IssmDouble B_average,s_average,T_average=0.,P_average=0.;
350 int *doflist = NULL;
351 IssmDouble *xyz_list = NULL;
352
353 /*Fetch number of nodes and dof for this finite element*/
354 int numnodes = element->GetNumberOfNodes();
355
356 /*Fetch dof list and allocate solution vector*/
357 element->GetDofList(&doflist,NoneApproximationEnum,GsetEnum);
358 IssmDouble* values = xNew<IssmDouble>(numnodes);
359 IssmDouble* pressure = xNew<IssmDouble>(numnodes);
360 IssmDouble* surface = xNew<IssmDouble>(numnodes);
361 IssmDouble* B = xNew<IssmDouble>(numnodes);
362 IssmDouble* temperature = xNew<IssmDouble>(numnodes);
363 IssmDouble* waterfraction = xNew<IssmDouble>(numnodes);
364
365 /*Use the dof list to index into the solution vector: */
366 for(i=0;i<numnodes;i++){
367 values[i]=solution[doflist[i]];
368
369 /*Check solution*/
370 if(xIsNan<IssmDouble>(values[i])) _error_("NaN found in solution vector");
371 }
372
373 /*Get all inputs and parameters*/
374 element->GetInputValue(&converged,ConvergedEnum);
375 element->GetInputListOnNodes(&pressure[0],PressureEnum);
376 if(converged){
377 for(i=0;i<numnodes;i++){
378 element->EnthalpyToThermal(&temperature[i],&waterfraction[i],values[i],pressure[i]);
379 if(waterfraction[i]<0.) _error_("Negative water fraction found in solution vector");
380 if(waterfraction[i]>1.) _error_("Water fraction >1 found in solution vector");
381 }
382 element->AddInput(EnthalpyEnum,values,P1Enum);
383 element->AddInput(WaterfractionEnum,waterfraction,P1Enum);
384 element->AddInput(TemperatureEnum,temperature,P1Enum);
385
386 /*Update Rheology only if converged (we must make sure that the temperature is below melting point
387 * otherwise the rheology could be negative*/
388 element->FindParam(&rheology_law,MaterialsRheologyLawEnum);
389 element->GetInputListOnNodes(&surface[0],SurfaceEnum);
390 switch(rheology_law){
391 case NoneEnum:
392 /*Do nothing: B is not temperature dependent*/
393 break;
394 case PatersonEnum:
395 for(i=0;i<numnodes;i++) B[i]=Paterson(temperature[i]);
396 element->AddMaterialInput(MaterialsRheologyBEnum,&B[0],P1Enum);
397 break;
398 case ArrheniusEnum:
399 element->GetVerticesCoordinates(&xyz_list);
400 for(i=0;i<numnodes;i++) B[i]=Arrhenius(temperature[i],surface[i]-xyz_list[i*3+2],element->GetMaterialParameter(MaterialsRheologyNEnum));
401 element->AddMaterialInput(MaterialsRheologyBEnum,&B[0],P1Enum);
402 break;
403 case LliboutryDuvalEnum:
404 for(i=0;i<numnodes;i++) B[i]=LliboutryDuval(values[i],pressure[i],element->GetMaterialParameter(MaterialsRheologyNEnum),element->GetMaterialParameter(MaterialsBetaEnum),element->GetMaterialParameter(ConstantsReferencetemperatureEnum),element->GetMaterialParameter(MaterialsHeatcapacityEnum),element->GetMaterialParameter(MaterialsLatentheatEnum));
405 element->AddMaterialInput(MaterialsRheologyBEnum,&B[0],P1Enum);
406 break;
407 default: _error_("Rheology law " << EnumToStringx(rheology_law) << " not supported yet");
408 }
409 }
410 else{
411 element->AddInput(EnthalpyPicardEnum,values,P1Enum);
412 }
413
414 /*Free ressources:*/
415 xDelete<IssmDouble>(values);
416 xDelete<IssmDouble>(pressure);
417 xDelete<IssmDouble>(surface);
418 xDelete<IssmDouble>(B);
419 xDelete<IssmDouble>(temperature);
420 xDelete<IssmDouble>(waterfraction);
421 xDelete<IssmDouble>(xyz_list);
422 xDelete<int>(doflist);
423}/*}}}*/
Note: See TracBrowser for help on using the repository browser.