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

Last change on this file since 17981 was 17981, checked in by Mathieu Morlighem, 11 years ago

BUG: fixed some memory leaks
File size: 49.7 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 domaintype,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 parameters->AddObject(iomodel->CopyConstantObject(FrictionLawEnum));
20
21 iomodel->FetchData(&requestedoutputs,&numoutputs,ThermalRequestedOutputsEnum);
22 parameters->AddObject(new IntParam(ThermalNumRequestedOutputsEnum,numoutputs));
23 if(numoutputs)parameters->AddObject(new StringArrayParam(ThermalRequestedOutputsEnum,requestedoutputs,numoutputs));
24 iomodel->DeleteData(&requestedoutputs,numoutputs,ThermalRequestedOutputsEnum);
25}/*}}}*/
26void EnthalpyAnalysis::UpdateElements(Elements* elements,IoModel* iomodel,int analysis_counter,int analysis_type){/*{{{*/
27
28 bool dakota_analysis,islevelset,isenthalpy;
29 int frictionlaw;
30
31 /*Now, is the model 3d? otherwise, do nothing: */
32 if(iomodel->domaintype==Domain2DhorizontalEnum)return;
33
34 /*Is enthalpy requested?*/
35 iomodel->Constant(&isenthalpy,ThermalIsenthalpyEnum);
36 if(!isenthalpy) return;
37
38 /*Fetch data needed: */
39 iomodel->FetchData(3,TemperatureEnum,WaterfractionEnum,PressureEnum);
40
41 /*Update elements: */
42 int counter=0;
43 for(int i=0;i<iomodel->numberofelements;i++){
44 if(iomodel->my_elements[i]){
45 Element* element=(Element*)elements->GetObjectByOffset(counter);
46 element->Update(i,iomodel,analysis_counter,analysis_type,P1Enum);
47 counter++;
48 }
49 }
50
51 iomodel->Constant(&dakota_analysis,QmuIsdakotaEnum);
52 iomodel->Constant(&islevelset,TransientIslevelsetEnum);
53 iomodel->Constant(&frictionlaw,FrictionLawEnum);
54
55 iomodel->FetchDataToInput(elements,ThicknessEnum);
56 iomodel->FetchDataToInput(elements,SurfaceEnum);
57 iomodel->FetchDataToInput(elements,BaseEnum);
58 iomodel->FetchDataToInput(elements,MaskIceLevelsetEnum);
59 iomodel->FetchDataToInput(elements,MaskGroundediceLevelsetEnum);
60 if(iomodel->domaintype!=Domain2DhorizontalEnum){
61 iomodel->FetchDataToInput(elements,MeshVertexonbaseEnum);
62 iomodel->FetchDataToInput(elements,MeshVertexonsurfaceEnum);
63 }
64 iomodel->FetchDataToInput(elements,FlowequationElementEquationEnum);
65 iomodel->FetchDataToInput(elements,MaterialsRheologyBEnum);
66 iomodel->FetchDataToInput(elements,MaterialsRheologyNEnum);
67 iomodel->FetchDataToInput(elements,PressureEnum);
68 iomodel->FetchDataToInput(elements,TemperatureEnum);
69 iomodel->FetchDataToInput(elements,WaterfractionEnum);
70 iomodel->FetchDataToInput(elements,EnthalpyEnum);
71 iomodel->FetchDataToInput(elements,BasalforcingsGeothermalfluxEnum);
72 iomodel->FetchDataToInput(elements,WatercolumnEnum);
73 iomodel->FetchDataToInput(elements,BasalforcingsMeltingRateEnum);
74 iomodel->FetchDataToInput(elements,VxEnum);
75 iomodel->FetchDataToInput(elements,VyEnum);
76 iomodel->FetchDataToInput(elements,VzEnum);
77 InputUpdateFromConstantx(elements,0.,VxMeshEnum);
78 InputUpdateFromConstantx(elements,0.,VyMeshEnum);
79 InputUpdateFromConstantx(elements,0.,VzMeshEnum);
80 if(dakota_analysis){
81 elements->InputDuplicate(TemperatureEnum,QmuTemperatureEnum);
82 elements->InputDuplicate(BasalforcingsMeltingRateEnum,QmuMeltingEnum);
83 elements->InputDuplicate(VxMeshEnum,QmuVxMeshEnum);
84 elements->InputDuplicate(VxMeshEnum,QmuVyMeshEnum);
85 elements->InputDuplicate(VxMeshEnum,QmuVzMeshEnum);
86 }
87 if(islevelset){
88 iomodel->FetchDataToInput(elements,IceMaskNodeActivationEnum);
89 iomodel->FetchDataToInput(elements,MeshVertexonbaseEnum); // required for updating active nodes
90 }
91
92 /*Friction law variables*/
93 switch(frictionlaw){
94 case 1:
95 iomodel->FetchDataToInput(elements,FrictionCoefficientEnum);
96 iomodel->FetchDataToInput(elements,FrictionPEnum);
97 iomodel->FetchDataToInput(elements,FrictionQEnum);
98 break;
99 case 2:
100 iomodel->FetchDataToInput(elements,FrictionCEnum);
101 iomodel->FetchDataToInput(elements,FrictionMEnum);
102 break;
103 default:
104 _error_("not supported");
105 }
106 /*Free data: */
107 iomodel->DeleteData(3,TemperatureEnum,WaterfractionEnum,PressureEnum);
108}/*}}}*/
109void EnthalpyAnalysis::CreateNodes(Nodes* nodes,IoModel* iomodel){/*{{{*/
110
111 if(iomodel->domaintype==Domain3DEnum) iomodel->FetchData(2,MeshVertexonbaseEnum,MeshVertexonsurfaceEnum);
112 ::CreateNodes(nodes,iomodel,EnthalpyAnalysisEnum,P1Enum);
113 iomodel->DeleteData(2,MeshVertexonbaseEnum,MeshVertexonsurfaceEnum);
114}/*}}}*/
115void EnthalpyAnalysis::CreateConstraints(Constraints* constraints,IoModel* iomodel){/*{{{*/
116
117 /*Intermediary*/
118 int count;
119 int M,N;
120 bool spcpresent = false;
121 IssmDouble heatcapacity;
122 IssmDouble referencetemperature;
123
124 /*Output*/
125 IssmDouble *spcvector = NULL;
126 IssmDouble* times=NULL;
127 IssmDouble* values=NULL;
128
129 /*Fetch parameters: */
130 iomodel->Constant(&heatcapacity,MaterialsHeatcapacityEnum);
131 iomodel->Constant(&referencetemperature,ConstantsReferencetemperatureEnum);
132
133 /*return if 2d mesh*/
134 if(iomodel->domaintype==Domain2DhorizontalEnum) return;
135
136 /*Fetch data: */
137 iomodel->FetchData(&spcvector,&M,&N,ThermalSpctemperatureEnum);
138
139 //FIX ME: SHOULD USE IOMODELCREATECONSTRAINTS
140 /*Transient or static?:*/
141 if(M==iomodel->numberofvertices){
142 /*static: just create Constraints objects*/
143 count=0;
144
145 for(int i=0;i<iomodel->numberofvertices;i++){
146 /*keep only this partition's nodes:*/
147 if((iomodel->my_vertices[i])){
148
149 if (!xIsNan<IssmDouble>(spcvector[i])){
150
151 constraints->AddObject(new SpcStatic(iomodel->constraintcounter+count+1,iomodel->nodecounter+i+1,0,heatcapacity*(spcvector[i]-referencetemperature),EnthalpyAnalysisEnum));
152 count++;
153
154 }
155 }
156 }
157 }
158 else if (M==(iomodel->numberofvertices+1)){
159 /*transient: create transient SpcTransient objects. Same logic, except we need to retrieve
160 * various times and values to initialize an SpcTransient object: */
161 count=0;
162
163 /*figure out times: */
164 times=xNew<IssmDouble>(N);
165 for(int j=0;j<N;j++){
166 times[j]=spcvector[(M-1)*N+j];
167 }
168
169 /*Create constraints from x,y,z: */
170 for(int i=0;i<iomodel->numberofvertices;i++){
171
172 /*keep only this partition's nodes:*/
173 if((iomodel->my_vertices[i])){
174
175 /*figure out times and values: */
176 values=xNew<IssmDouble>(N);
177 spcpresent=false;
178 for(int j=0;j<N;j++){
179 values[j]=heatcapacity*(spcvector[i*N+j]-referencetemperature);
180 if(!xIsNan<IssmDouble>(values[j]))spcpresent=true; //NaN means no spc by default
181 }
182
183 if(spcpresent){
184 constraints->AddObject(new SpcTransient(iomodel->constraintcounter+count+1,iomodel->nodecounter+i+1,0,N,times,values,EnthalpyAnalysisEnum));
185 count++;
186 }
187 xDelete<IssmDouble>(values);
188 }
189 }
190 }
191 else{
192 _error_("Size of field " << EnumToStringx(ThermalSpctemperatureEnum) << " not supported");
193 }
194
195 /*Free ressources:*/
196 iomodel->DeleteData(spcvector,ThermalSpctemperatureEnum);
197 xDelete<IssmDouble>(times);
198 xDelete<IssmDouble>(values);
199}/*}}}*/
200void EnthalpyAnalysis::CreateLoads(Loads* loads, IoModel* iomodel){/*{{{*/
201
202 /*No loads */
203}/*}}}*/
204
205/*Finite Element Analysis*/
206void EnthalpyAnalysis::Core(FemModel* femmodel){/*{{{*/
207 _error_("not implemented");
208}/*}}}*/
209ElementVector* EnthalpyAnalysis::CreateDVector(Element* element){/*{{{*/
210 /*Default, return NULL*/
211 return NULL;
212}/*}}}*/
213ElementMatrix* EnthalpyAnalysis::CreateJacobianMatrix(Element* element){/*{{{*/
214_error_("Not implemented");
215}/*}}}*/
216ElementMatrix* EnthalpyAnalysis::CreateKMatrix(Element* element){/*{{{*/
217
218 /* Check if ice in element */
219 if(!element->IsIceInElement()) return NULL;
220
221 /*compute all stiffness matrices for this element*/
222 ElementMatrix* Ke1=CreateKMatrixVolume(element);
223 ElementMatrix* Ke2=CreateKMatrixShelf(element);
224 ElementMatrix* Ke =new ElementMatrix(Ke1,Ke2);
225
226 /*clean-up and return*/
227 delete Ke1;
228 delete Ke2;
229 return Ke;
230}/*}}}*/
231ElementMatrix* EnthalpyAnalysis::CreateKMatrixVolume(Element* element){/*{{{*/
232
233 /* Check if ice in element */
234 if(!element->IsIceInElement()) return NULL;
235
236 /*Intermediaries */
237 int stabilization;
238 IssmDouble Jdet,dt,u,v,w,um,vm,wm,vel;
239 IssmDouble h,hx,hy,hz,vx,vy,vz;
240 IssmDouble tau_parameter,diameter;
241 IssmDouble D_scalar;
242 IssmDouble* xyz_list = NULL;
243
244 /*Fetch number of nodes and dof for this finite element*/
245 int numnodes = element->GetNumberOfNodes();
246
247 /*Initialize Element vector and other vectors*/
248 ElementMatrix* Ke = element->NewElementMatrix();
249 IssmDouble* basis = xNew<IssmDouble>(numnodes);
250 IssmDouble* dbasis = xNew<IssmDouble>(3*numnodes);
251 IssmDouble* B = xNew<IssmDouble>(3*numnodes);
252 IssmDouble* Bprime = xNew<IssmDouble>(3*numnodes);
253 IssmDouble D[3][3] = {0.};
254 IssmDouble K[3][3];
255
256 /*Retrieve all inputs and parameters*/
257 element->GetVerticesCoordinates(&xyz_list);
258 element->FindParam(&dt,TimesteppingTimeStepEnum);
259 element->FindParam(&stabilization,ThermalStabilizationEnum);
260 IssmDouble rho_water = element->GetMaterialParameter(MaterialsRhoSeawaterEnum);
261 IssmDouble rho_ice = element->GetMaterialParameter(MaterialsRhoIceEnum);
262 IssmDouble gravity = element->GetMaterialParameter(ConstantsGEnum);
263 IssmDouble heatcapacity = element->GetMaterialParameter(MaterialsHeatcapacityEnum);
264 IssmDouble thermalconductivity = element->GetMaterialParameter(MaterialsThermalconductivityEnum);
265 Input* vx_input = element->GetInput(VxEnum); _assert_(vx_input);
266 Input* vy_input = element->GetInput(VyEnum); _assert_(vy_input);
267 Input* vz_input = element->GetInput(VzEnum); _assert_(vz_input);
268 Input* vxm_input = element->GetInput(VxMeshEnum); _assert_(vxm_input);
269 Input* vym_input = element->GetInput(VyMeshEnum); _assert_(vym_input);
270 Input* vzm_input = element->GetInput(VzMeshEnum); _assert_(vzm_input);
271 if(stabilization==2) diameter=element->MinEdgeLength(xyz_list);
272
273 /*Enthalpy diffusion parameter*/
274 IssmDouble kappa=this->EnthalpyDiffusionParameterVolume(element,EnthalpyPicardEnum); _assert_(kappa>=0.);
275
276 /* Start looping on the number of gaussian points: */
277 Gauss* gauss=element->NewGauss(2);
278 for(int ig=gauss->begin();ig<gauss->end();ig++){
279 gauss->GaussPoint(ig);
280
281 element->JacobianDeterminant(&Jdet,xyz_list,gauss);
282 D_scalar=gauss->weight*Jdet;
283 if(dt!=0.) D_scalar=D_scalar*dt;
284
285 /*Conduction: */
286 GetBConduct(B,element,xyz_list,gauss);
287 D[0][0]=D_scalar*kappa/rho_ice;
288 D[1][1]=D_scalar*kappa/rho_ice;
289 D[2][2]=D_scalar*kappa/rho_ice;
290 TripleMultiply(B,3,numnodes,1,
291 &D[0][0],3,3,0,
292 B,3,numnodes,0,
293 &Ke->values[0],1);
294
295 /*Advection: */
296 GetBAdvec(B,element,xyz_list,gauss);
297 GetBAdvecprime(Bprime,element,xyz_list,gauss);
298 vx_input->GetInputValue(&u,gauss); vxm_input->GetInputValue(&um,gauss); vx=u-um;
299 vy_input->GetInputValue(&v,gauss); vym_input->GetInputValue(&vm,gauss); vy=v-vm;
300 vz_input->GetInputValue(&w,gauss); vzm_input->GetInputValue(&wm,gauss); vz=w-wm;
301 D[0][0]=D_scalar*vx;
302 D[1][1]=D_scalar*vy;
303 D[2][2]=D_scalar*vz;
304 TripleMultiply(B,3,numnodes,1,
305 &D[0][0],3,3,0,
306 Bprime,3,numnodes,0,
307 &Ke->values[0],1);
308
309 /*Transient: */
310 if(dt!=0.){
311 D_scalar=gauss->weight*Jdet;
312 element->NodalFunctions(basis,gauss);
313 TripleMultiply(basis,numnodes,1,0,
314 &D_scalar,1,1,0,
315 basis,1,numnodes,0,
316 &Ke->values[0],1);
317 D_scalar=D_scalar*dt;
318 }
319
320 /*Artifficial diffusivity*/
321 if(stabilization==1){
322 element->ElementSizes(&hx,&hy,&hz);
323 vel=sqrt(vx*vx + vy*vy + vz*vz)+1.e-14;
324 h=sqrt( pow(hx*vx/vel,2) + pow(hy*vy/vel,2) + pow(hz*vz/vel,2));
325 K[0][0]=h/(2.*vel)*vx*vx; K[0][1]=h/(2.*vel)*vx*vy; K[0][2]=h/(2.*vel)*vx*vz;
326 K[1][0]=h/(2.*vel)*vy*vx; K[1][1]=h/(2.*vel)*vy*vy; K[1][2]=h/(2.*vel)*vy*vz;
327 K[2][0]=h/(2.*vel)*vz*vx; K[2][1]=h/(2.*vel)*vz*vy; K[2][2]=h/(2.*vel)*vz*vz;
328 for(int i=0;i<3;i++) for(int j=0;j<3;j++) K[i][j] = D_scalar*K[i][j];
329
330 GetBAdvecprime(Bprime,element,xyz_list,gauss);
331 TripleMultiply(Bprime,3,numnodes,1,
332 &K[0][0],3,3,0,
333 Bprime,3,numnodes,0,
334 &Ke->values[0],1);
335 }
336 else if(stabilization==2){
337 element->NodalFunctionsDerivatives(dbasis,xyz_list,gauss);
338 tau_parameter=element->StabilizationParameter(u-um,v-vm,w-wm,diameter,kappa/rho_ice);
339 for(int i=0;i<numnodes;i++){
340 for(int j=0;j<numnodes;j++){
341 Ke->values[i*numnodes+j]+=tau_parameter*D_scalar*
342 ((u-um)*dbasis[0*numnodes+i]+(v-vm)*dbasis[1*numnodes+i]+(w-wm)*dbasis[2*numnodes+i])*((u-um)*dbasis[0*numnodes+j]+(v-vm)*dbasis[1*numnodes+j]+(w-wm)*dbasis[2*numnodes+j]);
343 }
344 }
345 if(dt!=0.){
346 D_scalar=gauss->weight*Jdet;
347 for(int i=0;i<numnodes;i++){
348 for(int j=0;j<numnodes;j++){
349 Ke->values[i*numnodes+j]+=tau_parameter*D_scalar*basis[j]*((u-um)*dbasis[0*numnodes+i]+(v-vm)*dbasis[1*numnodes+i]+(w-wm)*dbasis[2*numnodes+i]);
350 }
351 }
352 }
353 }
354 }
355
356 /*Clean up and return*/
357 xDelete<IssmDouble>(xyz_list);
358 xDelete<IssmDouble>(basis);
359 xDelete<IssmDouble>(dbasis);
360 xDelete<IssmDouble>(B);
361 xDelete<IssmDouble>(Bprime);
362 delete gauss;
363 return Ke;
364}/*}}}*/
365ElementMatrix* EnthalpyAnalysis::CreateKMatrixShelf(Element* element){/*{{{*/
366
367 /* Check if ice in element */
368 if(!element->IsIceInElement()) return NULL;
369
370 /*Initialize Element matrix and return if necessary*/
371 if(!element->IsOnBase() || !element->IsFloating()) return NULL;
372
373 /*Intermediaries*/
374 IssmDouble dt,Jdet,D;
375 IssmDouble *xyz_list_base = NULL;
376
377 /*Fetch number of nodes for this finite element*/
378 int numnodes = element->GetNumberOfNodes();
379
380 /*Initialize vectors*/
381 ElementMatrix* Ke = element->NewElementMatrix();
382 IssmDouble* basis = xNew<IssmDouble>(numnodes);
383
384 /*Retrieve all inputs and parameters*/
385 element->GetVerticesCoordinatesBase(&xyz_list_base);
386 element->FindParam(&dt,TimesteppingTimeStepEnum);
387 IssmDouble gravity = element->GetMaterialParameter(ConstantsGEnum);
388 IssmDouble rho_water = element->GetMaterialParameter(MaterialsRhoSeawaterEnum);
389 IssmDouble rho_ice = element->GetMaterialParameter(MaterialsRhoIceEnum);
390 IssmDouble heatcapacity = element->GetMaterialParameter(MaterialsHeatcapacityEnum);
391 IssmDouble mixed_layer_capacity= element->GetMaterialParameter(MaterialsMixedLayerCapacityEnum);
392 IssmDouble thermal_exchange_vel= element->GetMaterialParameter(MaterialsThermalExchangeVelocityEnum);
393
394 /* Start looping on the number of gaussian points: */
395 Gauss* gauss=element->NewGaussBase(2);
396 for(int ig=gauss->begin();ig<gauss->end();ig++){
397 gauss->GaussPoint(ig);
398
399 element->JacobianDeterminantBase(&Jdet,xyz_list_base,gauss);
400 element->NodalFunctions(basis,gauss);
401
402 D=gauss->weight*Jdet*rho_water*mixed_layer_capacity*thermal_exchange_vel/(heatcapacity*rho_ice);
403 if(reCast<bool,IssmDouble>(dt)) D=dt*D;
404 TripleMultiply(basis,numnodes,1,0,
405 &D,1,1,0,
406 basis,1,numnodes,0,
407 &Ke->values[0],1);
408
409 }
410
411 /*Clean up and return*/
412 delete gauss;
413 xDelete<IssmDouble>(basis);
414 xDelete<IssmDouble>(xyz_list_base);
415 return Ke;
416}/*}}}*/
417ElementVector* EnthalpyAnalysis::CreatePVector(Element* element){/*{{{*/
418
419 /* Check if ice in element */
420 if(!element->IsIceInElement()) return NULL;
421
422 /*compute all load vectors for this element*/
423 ElementVector* pe1=CreatePVectorVolume(element);
424 ElementVector* pe2=CreatePVectorSheet(element);
425 ElementVector* pe3=CreatePVectorShelf(element);
426 ElementVector* pe =new ElementVector(pe1,pe2,pe3);
427
428 /*clean-up and return*/
429 delete pe1;
430 delete pe2;
431 delete pe3;
432 return pe;
433}/*}}}*/
434ElementVector* EnthalpyAnalysis::CreatePVectorVolume(Element* element){/*{{{*/
435
436 /* Check if ice in element */
437 if(!element->IsIceInElement()) return NULL;
438
439 /*Intermediaries*/
440 int i, stabilization;
441 IssmDouble Jdet,phi,dt;
442 IssmDouble enthalpy, Hpmp;
443 IssmDouble enthalpypicard, d1enthalpypicard[3];
444 IssmDouble pressure, d1pressure[3], d2pressure;
445 IssmDouble waterfractionpicard;
446 IssmDouble kappa,tau_parameter,diameter,kappa_w;
447 IssmDouble u,v,w;
448 IssmDouble scalar_def, scalar_sens ,scalar_transient;
449 IssmDouble* xyz_list = NULL;
450 IssmDouble d1H_d1P, d1P2;
451
452 /*Fetch number of nodes and dof for this finite element*/
453 int numnodes = element->GetNumberOfNodes();
454 int numvertices = element->GetNumberOfVertices();
455
456 /*Initialize Element vector*/
457 ElementVector* pe = element->NewElementVector();
458 IssmDouble* basis = xNew<IssmDouble>(numnodes);
459 IssmDouble* dbasis = xNew<IssmDouble>(3*numnodes);
460
461 /*Retrieve all inputs and parameters*/
462 element->GetVerticesCoordinates(&xyz_list);
463 IssmDouble rho_ice = element->GetMaterialParameter(MaterialsRhoIceEnum);
464 IssmDouble heatcapacity = element->GetMaterialParameter(MaterialsHeatcapacityEnum);
465 IssmDouble thermalconductivity = element->GetMaterialParameter(MaterialsThermalconductivityEnum);
466 IssmDouble temperateiceconductivity = element->GetMaterialParameter(MaterialsTemperateiceconductivityEnum);
467 IssmDouble beta = element->GetMaterialParameter(MaterialsBetaEnum);
468 IssmDouble latentheat = element->GetMaterialParameter(MaterialsLatentheatEnum);
469 element->FindParam(&dt,TimesteppingTimeStepEnum);
470 element->FindParam(&stabilization,ThermalStabilizationEnum);
471 Input* vx_input=element->GetInput(VxEnum); _assert_(vx_input);
472 Input* vy_input=element->GetInput(VyEnum); _assert_(vy_input);
473 Input* vz_input=element->GetInput(VzEnum); _assert_(vz_input);
474 Input* enthalpypicard_input=element->GetInput(EnthalpyPicardEnum); _assert_(enthalpypicard_input);
475 Input* pressure_input=element->GetInput(PressureEnum); _assert_(pressure_input);
476 Input* enthalpy_input=NULL;
477 if(reCast<bool,IssmDouble>(dt)){enthalpy_input = element->GetInput(EnthalpyEnum); _assert_(enthalpy_input);}
478 if(stabilization==2){
479 diameter=element->MinEdgeLength(xyz_list);
480 kappa=this->EnthalpyDiffusionParameterVolume(element,EnthalpyPicardEnum); _assert_(kappa>=0.);
481 }
482
483 /* Start looping on the number of gaussian points: */
484 Gauss* gauss=element->NewGauss(3);
485 for(int ig=gauss->begin();ig<gauss->end();ig++){
486 gauss->GaussPoint(ig);
487
488 element->JacobianDeterminant(&Jdet,xyz_list,gauss);
489 element->NodalFunctions(basis,gauss);
490
491 /*viscous dissipation*/
492 element->ViscousHeating(&phi,xyz_list,gauss,vx_input,vy_input,vz_input);
493
494 scalar_def=phi/rho_ice*Jdet*gauss->weight;
495 if(dt!=0.) scalar_def=scalar_def*dt;
496
497 for(i=0;i<numnodes;i++) pe->values[i]+=scalar_def*basis[i];
498
499 /*sensible heat flux in temperate ice*/
500 enthalpypicard_input->GetInputValue(&enthalpypicard,gauss);
501 pressure_input->GetInputValue(&pressure,gauss);
502 Hpmp=this->PureIceEnthalpy(element, pressure);
503
504 if(enthalpypicard>=Hpmp){
505 enthalpypicard_input->GetInputDerivativeValue(&d1enthalpypicard[0],xyz_list,gauss);
506 pressure_input->GetInputDerivativeValue(&d1pressure[0],xyz_list,gauss);
507 d2pressure=0.; // for linear elements, 2nd derivative is zero
508
509 d1H_d1P=0.;
510 for(i=0;i<3;i++) d1H_d1P+=d1enthalpypicard[i]*d1pressure[i];
511 d1P2=0.;
512 for(i=0;i<3;i++) d1P2+=pow(d1pressure[i],2.);
513
514 scalar_sens=-beta*((temperateiceconductivity - thermalconductivity)/latentheat*(d1H_d1P + beta*heatcapacity*d1P2))/rho_ice;
515 if(dt!=0.) scalar_sens=scalar_sens*dt;
516 for(i=0;i<numnodes;i++) pe->values[i]+=scalar_sens*basis[i];
517 }
518
519 /* Build transient now */
520 if(reCast<bool,IssmDouble>(dt)){
521 enthalpy_input->GetInputValue(&enthalpy, gauss);
522 scalar_transient=enthalpy*Jdet*gauss->weight;
523 for(i=0;i<numnodes;i++) pe->values[i]+=scalar_transient*basis[i];
524 }
525
526 if(stabilization==2){
527 element->NodalFunctionsDerivatives(dbasis,xyz_list,gauss);
528
529 vx_input->GetInputValue(&u,gauss);
530 vy_input->GetInputValue(&v,gauss);
531 vz_input->GetInputValue(&w,gauss);
532 tau_parameter=element->StabilizationParameter(u,v,w,diameter,kappa/rho_ice);
533
534 for(i=0;i<numnodes;i++) pe->values[i]+=tau_parameter*scalar_def*(u*dbasis[0*numnodes+i]+v*dbasis[1*numnodes+i]+w*dbasis[2*numnodes+i]);
535
536 if(dt!=0.){
537 for(i=0;i<numnodes;i++) pe->values[i]+=tau_parameter*scalar_transient*(u*dbasis[0*numnodes+i]+v*dbasis[1*numnodes+i]+w*dbasis[2*numnodes+i]);
538 }
539 }
540 }
541
542 /*Clean up and return*/
543 xDelete<IssmDouble>(basis);
544 xDelete<IssmDouble>(dbasis);
545 xDelete<IssmDouble>(xyz_list);
546 delete gauss;
547 return pe;
548
549}/*}}}*/
550ElementVector* EnthalpyAnalysis::CreatePVectorSheet(Element* element){/*{{{*/
551
552 /* Check if ice in element */
553 if(!element->IsIceInElement()) return NULL;
554
555 /* implementation of the basal condition decision chart of Aschwanden 2012, Fig.5 */
556 if(!element->IsOnBase() || element->IsFloating()) return NULL;
557
558 IssmDouble dt,Jdet,enthalpy,pressure,watercolumn,geothermalflux,vx,vy,vz;
559 IssmDouble enthalpyup,pressureup,alpha2,scalar,basalfriction,heatflux;
560 IssmDouble *xyz_list_base = NULL;
561
562 /*Fetch number of nodes for this finite element*/
563 int numnodes = element->GetNumberOfNodes();
564
565 /*Initialize vectors*/
566 ElementVector* pe = element->NewElementVector();
567 IssmDouble* basis = xNew<IssmDouble>(numnodes);
568
569 /*Retrieve all inputs and parameters*/
570 element->GetVerticesCoordinatesBase(&xyz_list_base);
571 element->FindParam(&dt,TimesteppingTimeStepEnum);
572 Input* vx_input = element->GetInput(VxEnum); _assert_(vx_input);
573 Input* vy_input = element->GetInput(VyEnum); _assert_(vy_input);
574 Input* vz_input = element->GetInput(VzEnum); _assert_(vz_input);
575 Input* enthalpy_input = element->GetInput(EnthalpyPicardEnum); _assert_(enthalpy_input);
576 Input* pressure_input = element->GetInput(PressureEnum); _assert_(pressure_input);
577 Input* geothermalflux_input = element->GetInput(BasalforcingsGeothermalfluxEnum); _assert_(geothermalflux_input);
578 Input* watercolumn_input = element->GetInput(WatercolumnEnum); _assert_(watercolumn_input);
579 IssmDouble rho_ice = element->GetMaterialParameter(MaterialsRhoIceEnum);
580 IssmDouble heatcapacity = element->GetMaterialParameter(MaterialsHeatcapacityEnum);
581
582 /*Build friction element, needed later: */
583 Friction* friction=new Friction(element,3);
584
585 /* Start looping on the number of gaussian points: */
586 Gauss* gauss = element->NewGaussBase(2);
587 Gauss* gaussup = element->NewGaussTop(2);
588 for(int ig=gauss->begin();ig<gauss->end();ig++){
589 gauss->GaussPoint(ig);
590
591 element->JacobianDeterminantBase(&Jdet,xyz_list_base,gauss);
592 element->NodalFunctions(basis,gauss);
593
594 enthalpy_input->GetInputValue(&enthalpy,gauss);
595 pressure_input->GetInputValue(&pressure,gauss);
596 watercolumn_input->GetInputValue(&watercolumn,gauss);
597
598 if((watercolumn<=0.) && (enthalpy<PureIceEnthalpy(element,pressure))){
599 /* the above check is equivalent to
600 NOT [(watercolumn>0.) AND (enthalpy<PIE)] AND (enthalpy<PIE)*/
601 geothermalflux_input->GetInputValue(&geothermalflux,gauss);
602
603 friction->GetAlpha2(&alpha2,gauss);
604 vx_input->GetInputValue(&vx,gauss);
605 vy_input->GetInputValue(&vy,gauss);
606 vz_input->GetInputValue(&vz,gauss);
607 basalfriction = alpha2*(vx*vx + vy*vy + vz*vz);
608 heatflux = (basalfriction+geothermalflux)/(rho_ice);
609
610 scalar = gauss->weight*Jdet*heatflux;
611 if(dt!=0.) scalar=dt*scalar;
612
613 for(int i=0;i<numnodes;i++) pe->values[i]+=scalar*basis[i];
614 }
615 else if(enthalpy >= PureIceEnthalpy(element,pressure)){
616 /* check positive thickness of temperate basal ice layer */
617 enthalpy_input->GetInputValue(&enthalpyup,gaussup);
618 pressure_input->GetInputValue(&pressureup,gaussup);
619 if(enthalpyup >= PureIceEnthalpy(element,pressureup)){
620 // do nothing, set grad enthalpy*n=0.
621 }
622 else{
623 // only base temperate, set Dirichlet BCs in Penta::UpdateBasalConstraintsEnthalpy()
624 }
625 }
626 else{
627 // base cold, but watercolumn positive. Set base to pressure melting point enthalpy
628 }
629 }
630
631 /*Clean up and return*/
632 delete gauss;
633 delete gaussup;
634 delete friction;
635 xDelete<IssmDouble>(basis);
636 xDelete<IssmDouble>(xyz_list_base);
637 return pe;
638
639}/*}}}*/
640ElementVector* EnthalpyAnalysis::CreatePVectorShelf(Element* element){/*{{{*/
641
642 /* Check if ice in element */
643 if(!element->IsIceInElement()) return NULL;
644
645 /*Get basal element*/
646 if(!element->IsOnBase() || !element->IsFloating()) return NULL;
647
648 IssmDouble h_pmp,dt,Jdet,scalar_ocean,pressure;
649 IssmDouble *xyz_list_base = NULL;
650
651 /*Fetch number of nodes for this finite element*/
652 int numnodes = element->GetNumberOfNodes();
653
654 /*Initialize vectors*/
655 ElementVector* pe = element->NewElementVector();
656 IssmDouble* basis = xNew<IssmDouble>(numnodes);
657
658 /*Retrieve all inputs and parameters*/
659 element->GetVerticesCoordinatesBase(&xyz_list_base);
660 element->FindParam(&dt,TimesteppingTimeStepEnum);
661 Input* pressure_input=element->GetInput(PressureEnum); _assert_(pressure_input);
662 IssmDouble gravity = element->GetMaterialParameter(ConstantsGEnum);
663 IssmDouble rho_water = element->GetMaterialParameter(MaterialsRhoSeawaterEnum);
664 IssmDouble rho_ice = element->GetMaterialParameter(MaterialsRhoIceEnum);
665 IssmDouble heatcapacity = element->GetMaterialParameter(MaterialsHeatcapacityEnum);
666 IssmDouble mixed_layer_capacity= element->GetMaterialParameter(MaterialsMixedLayerCapacityEnum);
667 IssmDouble thermal_exchange_vel= element->GetMaterialParameter(MaterialsThermalExchangeVelocityEnum);
668
669 /* Start looping on the number of gaussian points: */
670 Gauss* gauss=element->NewGaussBase(2);
671 for(int ig=gauss->begin();ig<gauss->end();ig++){
672 gauss->GaussPoint(ig);
673
674 element->JacobianDeterminantBase(&Jdet,xyz_list_base,gauss);
675 element->NodalFunctions(basis,gauss);
676
677 pressure_input->GetInputValue(&pressure,gauss);
678 h_pmp=element->PureIceEnthalpy(pressure);
679
680 scalar_ocean=gauss->weight*Jdet*rho_water*mixed_layer_capacity*thermal_exchange_vel*h_pmp/(heatcapacity*rho_ice);
681 if(reCast<bool,IssmDouble>(dt)) scalar_ocean=dt*scalar_ocean;
682
683 for(int i=0;i<numnodes;i++) pe->values[i]+=scalar_ocean*basis[i];
684 }
685
686 /*Clean up and return*/
687 delete gauss;
688 xDelete<IssmDouble>(basis);
689 xDelete<IssmDouble>(xyz_list_base);
690 return pe;
691}/*}}}*/
692void EnthalpyAnalysis::GetBConduct(IssmDouble* B,Element* element,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/
693 /*Compute B matrix. B=[B1 B2 B3 B4 B5 B6] where Bi is of size 5*NDOF1.
694 * For node i, Bi' can be expressed in the actual coordinate system
695 * by:
696 * Bi_conduct=[ dh/dx ]
697 * [ dh/dy ]
698 * [ dh/dz ]
699 * where h is the interpolation function for node i.
700 *
701 * We assume B has been allocated already, of size: 3x(NDOF1*numnodes)
702 */
703
704 /*Fetch number of nodes for this finite element*/
705 int numnodes = element->GetNumberOfNodes();
706
707 /*Get nodal functions derivatives*/
708 IssmDouble* dbasis=xNew<IssmDouble>(3*numnodes);
709 element->NodalFunctionsDerivatives(dbasis,xyz_list,gauss);
710
711 /*Build B: */
712 for(int i=0;i<numnodes;i++){
713 B[numnodes*0+i] = dbasis[0*numnodes+i];
714 B[numnodes*1+i] = dbasis[1*numnodes+i];
715 B[numnodes*2+i] = dbasis[2*numnodes+i];
716 }
717
718 /*Clean-up*/
719 xDelete<IssmDouble>(dbasis);
720}/*}}}*/
721void EnthalpyAnalysis::GetBAdvec(IssmDouble* B,Element* element,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/
722 /*Compute B matrix. B=[B1 B2 B3 B4 B5 B6] where Bi is of size 5*NDOF1.
723 * For node i, Bi' can be expressed in the actual coordinate system
724 * by:
725 * Bi_advec =[ h ]
726 * [ h ]
727 * [ h ]
728 * where h is the interpolation function for node i.
729 *
730 * We assume B has been allocated already, of size: 3x(NDOF1*NUMNODESP1)
731 */
732
733 /*Fetch number of nodes for this finite element*/
734 int numnodes = element->GetNumberOfNodes();
735
736 /*Get nodal functions*/
737 IssmDouble* basis=xNew<IssmDouble>(numnodes);
738 element->NodalFunctions(basis,gauss);
739
740 /*Build B: */
741 for(int i=0;i<numnodes;i++){
742 B[numnodes*0+i] = basis[i];
743 B[numnodes*1+i] = basis[i];
744 B[numnodes*2+i] = basis[i];
745 }
746
747 /*Clean-up*/
748 xDelete<IssmDouble>(basis);
749}/*}}}*/
750void EnthalpyAnalysis::GetBAdvecprime(IssmDouble* B,Element* element,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/
751 /*Compute B matrix. B=[B1 B2 B3 B4 B5 B6] where Bi is of size 5*NDOF1.
752 * For node i, Bi' can be expressed in the actual coordinate system
753 * by:
754 * Biprime_advec=[ dh/dx ]
755 * [ dh/dy ]
756 * [ dh/dz ]
757 * where h is the interpolation function for node i.
758 *
759 * We assume B has been allocated already, of size: 3x(NDOF1*numnodes)
760 */
761
762 /*Fetch number of nodes for this finite element*/
763 int numnodes = element->GetNumberOfNodes();
764
765 /*Get nodal functions derivatives*/
766 IssmDouble* dbasis=xNew<IssmDouble>(3*numnodes);
767 element->NodalFunctionsDerivatives(dbasis,xyz_list,gauss);
768
769 /*Build B: */
770 for(int i=0;i<numnodes;i++){
771 B[numnodes*0+i] = dbasis[0*numnodes+i];
772 B[numnodes*1+i] = dbasis[1*numnodes+i];
773 B[numnodes*2+i] = dbasis[2*numnodes+i];
774 }
775
776 /*Clean-up*/
777 xDelete<IssmDouble>(dbasis);
778}/*}}}*/
779void EnthalpyAnalysis::GetSolutionFromInputs(Vector<IssmDouble>* solution,Element* element){/*{{{*/
780 element->GetSolutionFromInputsOneDof(solution,EnthalpyEnum);
781}/*}}}*/
782void EnthalpyAnalysis::InputUpdateFromSolution(IssmDouble* solution,Element* element){/*{{{*/
783
784 bool converged;
785 int i,rheology_law;
786 IssmDouble B_average,s_average,T_average=0.,P_average=0.;
787 int *doflist = NULL;
788 IssmDouble *xyz_list = NULL;
789
790 /*Fetch number of nodes and dof for this finite element*/
791 int numnodes = element->GetNumberOfNodes();
792
793 /*Fetch dof list and allocate solution vector*/
794 element->GetDofList(&doflist,NoneApproximationEnum,GsetEnum);
795 IssmDouble* values = xNew<IssmDouble>(numnodes);
796 IssmDouble* pressure = xNew<IssmDouble>(numnodes);
797 IssmDouble* surface = xNew<IssmDouble>(numnodes);
798 IssmDouble* B = xNew<IssmDouble>(numnodes);
799 IssmDouble* temperature = xNew<IssmDouble>(numnodes);
800 IssmDouble* waterfraction = xNew<IssmDouble>(numnodes);
801
802 /*Use the dof list to index into the solution vector: */
803 for(i=0;i<numnodes;i++){
804 values[i]=solution[doflist[i]];
805
806 /*Check solution*/
807 if(xIsNan<IssmDouble>(values[i])) _error_("NaN found in solution vector");
808 }
809
810 /*Get all inputs and parameters*/
811 element->GetInputValue(&converged,ConvergedEnum);
812 element->GetInputListOnNodes(&pressure[0],PressureEnum);
813 if(converged){
814 for(i=0;i<numnodes;i++){
815 element->EnthalpyToThermal(&temperature[i],&waterfraction[i],values[i],pressure[i]);
816 if(waterfraction[i]<0.) _error_("Negative water fraction found in solution vector");
817 //if(waterfraction[i]>1.) _error_("Water fraction >1 found in solution vector");
818 }
819 element->AddInput(EnthalpyEnum,values,element->GetElementType());
820 element->AddInput(WaterfractionEnum,waterfraction,element->GetElementType());
821 element->AddInput(TemperatureEnum,temperature,element->GetElementType());
822
823 /*Update Rheology only if converged (we must make sure that the temperature is below melting point
824 * otherwise the rheology could be negative*/
825 element->FindParam(&rheology_law,MaterialsRheologyLawEnum);
826 element->GetInputListOnNodes(&surface[0],SurfaceEnum);
827 switch(rheology_law){
828 case NoneEnum:
829 /*Do nothing: B is not temperature dependent*/
830 break;
831 case CuffeyEnum:
832 for(i=0;i<numnodes;i++) B[i]=Cuffey(temperature[i]);
833 element->AddInput(MaterialsRheologyBEnum,&B[0],element->GetElementType());
834 break;
835 case PatersonEnum:
836 for(i=0;i<numnodes;i++) B[i]=Paterson(temperature[i]);
837 element->AddInput(MaterialsRheologyBEnum,&B[0],element->GetElementType());
838 break;
839 case ArrheniusEnum:
840 element->GetVerticesCoordinates(&xyz_list);
841 for(i=0;i<numnodes;i++) B[i]=Arrhenius(temperature[i],surface[i]-xyz_list[i*3+2],element->GetMaterialParameter(MaterialsRheologyNEnum));
842 element->AddInput(MaterialsRheologyBEnum,&B[0],element->GetElementType());
843 break;
844 case LliboutryDuvalEnum:
845 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));
846 element->AddInput(MaterialsRheologyBEnum,&B[0],element->GetElementType());
847 break;
848 default: _error_("Rheology law " << EnumToStringx(rheology_law) << " not supported yet");
849 }
850 }
851 else{
852 element->AddInput(EnthalpyPicardEnum,values,element->GetElementType());
853 }
854
855 /*Free ressources:*/
856 xDelete<IssmDouble>(values);
857 xDelete<IssmDouble>(pressure);
858 xDelete<IssmDouble>(surface);
859 xDelete<IssmDouble>(B);
860 xDelete<IssmDouble>(temperature);
861 xDelete<IssmDouble>(waterfraction);
862 xDelete<IssmDouble>(xyz_list);
863 xDelete<int>(doflist);
864}/*}}}*/
865void EnthalpyAnalysis::UpdateConstraints(FemModel* femmodel){/*{{{*/
866
867 bool islevelset;
868 femmodel->parameters->FindParam(&islevelset,TransientIslevelsetEnum);
869 if(islevelset){
870 SetActiveNodesLSMx(femmodel);
871 }
872 return;
873}/*}}}*/
874
875/*Modules*/
876void EnthalpyAnalysis::PostProcessing(FemModel* femmodel){/*{{{*/
877
878 /*Intermediaries*/
879 int solution_type, i;
880 bool computebasalmeltingrates=true;
881 bool isdrainage=true;
882 bool updatebasalconstraints=true;
883
884 if(isdrainage){
885 /*Drain excess water fraction in ice column: */
886 for(i=0;i<femmodel->elements->Size();i++){
887 Element* element=dynamic_cast<Element*>(femmodel->elements->GetObjectByOffset(i));
888 DrainWaterfractionIcecolumn(element);
889 }
890 }
891
892 if(computebasalmeltingrates){
893 /*Compute basal melting rates: */
894 for(i=0;i<femmodel->elements->Size();i++){
895 Element* element=dynamic_cast<Element*>(femmodel->elements->GetObjectByOffset(i));
896 ComputeBasalMeltingrate(element);
897 }
898 }
899
900 if(updatebasalconstraints){
901 /*Update basal dirichlet BCs for enthalpy in transient runs: */
902 femmodel->parameters->FindParam(&solution_type,SolutionTypeEnum);
903 if(solution_type==TransientSolutionEnum){
904 for(i=0;i<femmodel->elements->Size();i++){
905 Element* element=dynamic_cast<Element*>(femmodel->elements->GetObjectByOffset(i));
906 UpdateBasalConstraints(element);
907 }
908 }
909 }
910}/*}}}*/
911void EnthalpyAnalysis::ComputeBasalMeltingrate(Element* element){/*{{{*/
912 /*Calculate the basal melt rates of the enthalpy model after Aschwanden 2012*/
913 /* melting rate is positive when melting, negative when refreezing*/
914
915 /* Check if ice in element */
916 if(!element->IsIceInElement()) return;
917
918 /* Only compute melt rates at the base of grounded ice*/
919 if(!element->IsOnBase() || element->IsFloating()) return;
920
921 /* Intermediaries */
922 const int dim=3;
923 int i,is,vertexdown,vertexup,numvertices,numsegments;
924 IssmDouble heatflux;
925 IssmDouble vec_heatflux[dim],normal_base[dim],d1enthalpy[dim];
926 IssmDouble basalfriction,alpha2;
927 IssmDouble dt,yts;
928 IssmDouble melting_overshoot,lambda;
929 IssmDouble geothermalflux;
930 IssmDouble vx,vy,vz;
931 IssmDouble *xyz_list = NULL;
932 IssmDouble *xyz_list_base = NULL;
933 int *pairindices = NULL;
934
935 /*Fetch parameters and inputs */
936 element->GetVerticesCoordinates(&xyz_list);
937 element->GetVerticesCoordinatesBase(&xyz_list_base);
938 IssmDouble latentheat = element->GetMaterialParameter(MaterialsLatentheatEnum);
939 IssmDouble rho_ice = element->GetMaterialParameter(MaterialsRhoIceEnum);
940 IssmDouble rho_water = element->GetMaterialParameter(MaterialsRhoFreshwaterEnum);
941 Input* enthalpy_input = element->GetInput(EnthalpyEnum); _assert_(enthalpy_input);
942 Input* geothermalflux_input = element->GetInput(BasalforcingsGeothermalfluxEnum); _assert_(geothermalflux_input);
943 Input* vx_input = element->GetInput(VxEnum); _assert_(vx_input);
944 Input* vy_input = element->GetInput(VyEnum); _assert_(vy_input);
945 Input* vz_input = element->GetInput(VzEnum); _assert_(vz_input);
946 IssmDouble kappa=EnthalpyDiffusionParameterVolume(element,EnthalpyEnum); _assert_(kappa>=0.);
947 element->NormalBase(&normal_base[0],xyz_list_base);
948 element->VerticalSegmentIndices(&pairindices,&numsegments);
949 IssmDouble* meltingrate_enthalpy = xNew<IssmDouble>(numsegments);
950 IssmDouble* heating = xNew<IssmDouble>(numsegments);
951
952 /*Build friction element, needed later: */
953 Friction* friction=new Friction(element,dim);
954
955 /******** MELTING RATES ************************************/
956 numvertices=element->GetNumberOfVertices();
957 IssmDouble* enthalpy = xNew<IssmDouble>(numvertices);
958 IssmDouble* pressure = xNew<IssmDouble>(numvertices);
959 IssmDouble* watercolumn = xNew<IssmDouble>(numvertices);
960 IssmDouble* basalmeltingrate = xNew<IssmDouble>(numvertices);
961 element->GetInputListOnVertices(enthalpy,EnthalpyEnum);
962 element->GetInputListOnVertices(pressure,PressureEnum);
963 element->GetInputListOnVertices(watercolumn,WatercolumnEnum);
964 element->GetInputListOnVertices(basalmeltingrate,BasalforcingsMeltingRateEnum);
965
966 Gauss* gauss=element->NewGauss();
967
968 for(int is=0;is<numsegments;is++){
969 vertexdown = pairindices[is*2+0];
970 vertexup = pairindices[is*2+1];
971 gauss->GaussVertex(vertexdown);
972
973 bool checkpositivethickness=true;
974 _assert_(watercolumn[vertexdown]>=0.);
975
976 /*Calculate basal meltingrate after Fig.5 of A.Aschwanden 2012*/
977 meltingrate_enthalpy[is]=0.;
978 heating[is]=0.;
979 if((watercolumn[vertexdown]>0.) && (enthalpy[vertexdown]<PureIceEnthalpy(element,pressure[vertexdown]))){
980 /*ensure that no ice is at T<Tm(p), if water layer present*/
981 enthalpy[vertexdown]=element->PureIceEnthalpy(pressure[vertexdown]);
982 }
983 else if(enthalpy[vertexdown]<element->PureIceEnthalpy(pressure[vertexdown])){
984 /*cold base: set q*n=q_geo*n+frictionheating as Neumann BC in Penta::CreatePVectorEnthalpySheet*/
985 checkpositivethickness=false; // cold base, skip next test
986 }
987 else{/*we have a temperate base, go to next test*/}
988
989 if(checkpositivethickness){
990 /*From here on all basal ice is temperate. Check for positive thickness of layer of temperate ice. */
991 bool istemperatelayer=false;
992 if(enthalpy[vertexup]>=element->PureIceEnthalpy(pressure[vertexup])) istemperatelayer=true;
993 if(istemperatelayer) for(i=0;i<dim;i++) vec_heatflux[i]=0.; // TODO: add -k*nabla T_pmp
994 else{
995 enthalpy_input->GetInputDerivativeValue(&d1enthalpy[0],xyz_list,gauss);
996 for(i=0;i<3;i++) vec_heatflux[i]=-kappa*d1enthalpy[i];
997 }
998
999 /*heat flux along normal*/
1000 heatflux=0.;
1001 for(i=0;i<3;i++) heatflux+=(vec_heatflux[i])*normal_base[i];
1002
1003 /*basal friction*/
1004 friction->GetAlpha2(&alpha2,gauss);
1005 vx_input->GetInputValue(&vx,gauss);
1006 vy_input->GetInputValue(&vy,gauss);
1007 vz_input->GetInputValue(&vz,gauss);
1008 basalfriction=alpha2*(vx*vx + vy*vy + vz*vz);
1009
1010 geothermalflux_input->GetInputValue(&geothermalflux,gauss);
1011 /* -Mb= Fb-(q-q_geo)/((1-w)*L*rho), and (1-w)*rho=rho_ice, cf Aschwanden 2012, eqs.1, 2, 66*/
1012 heating[is]=(heatflux+basalfriction+geothermalflux);
1013 meltingrate_enthalpy[is]=heating[is]/(latentheat*rho_ice); // m/s water equivalent
1014 }
1015 }
1016 /******** UPDATE MELTINGRATES AND WATERCOLUMN **************/
1017 element->FindParam(&dt,TimesteppingTimeStepEnum);
1018 for(is=0;is<numsegments;is++){
1019 vertexdown = pairindices[is*2+0];
1020 vertexup = pairindices[is*2+1];
1021 if(dt!=0.){
1022 if(watercolumn[vertexdown]+meltingrate_enthalpy[is]*dt<0.){ // prevent too much freeze on
1023 lambda = -watercolumn[vertexdown]/(dt*meltingrate_enthalpy[is]); _assert_(lambda>=0.); _assert_(lambda<1.);
1024 watercolumn[vertexdown]=0.;
1025 basalmeltingrate[vertexdown]=lambda*meltingrate_enthalpy[is]; // restrict freeze on only to size of watercolumn
1026 enthalpy[vertexdown]+=(1.-lambda)*meltingrate_enthalpy[is]*dt*latentheat; // use rest of energy to cool down base
1027 }
1028 else{
1029 basalmeltingrate[vertexdown]=meltingrate_enthalpy[is];
1030 watercolumn[vertexdown]+=dt*meltingrate_enthalpy[is];
1031 }
1032 }
1033 else{
1034 basalmeltingrate[vertexdown]=meltingrate_enthalpy[is];
1035 watercolumn[vertexdown]+=meltingrate_enthalpy[is];
1036 }
1037 basalmeltingrate[vertexdown]*=rho_water/rho_ice; // convert meltingrate from water to ice equivalent
1038 _assert_(watercolumn[vertexdown]>=0.);
1039 }
1040
1041 /*feed updated variables back into model*/
1042 element->AddInput(EnthalpyEnum,enthalpy,P1Enum);
1043 element->AddInput(WatercolumnEnum,watercolumn,P1Enum);
1044 element->AddInput(BasalforcingsMeltingRateEnum,basalmeltingrate,P1Enum);
1045
1046 /*Clean up and return*/
1047 delete gauss;
1048 delete friction;
1049 xDelete<int>(pairindices);
1050 xDelete<IssmDouble>(enthalpy);
1051 xDelete<IssmDouble>(pressure);
1052 xDelete<IssmDouble>(watercolumn);
1053 xDelete<IssmDouble>(basalmeltingrate);
1054 xDelete<IssmDouble>(meltingrate_enthalpy);
1055 xDelete<IssmDouble>(heating);
1056 xDelete<IssmDouble>(xyz_list);
1057 xDelete<IssmDouble>(xyz_list_base);
1058}/*}}}*/
1059void EnthalpyAnalysis::DrainWaterfractionIcecolumn(Element* element){/*{{{*/
1060
1061 /* Check if ice in element */
1062 if(!element->IsIceInElement()) return;
1063
1064 /* Only drain waterfraction of ice column from element at base*/
1065 if(!element->IsOnBase()) return; //FIXME: allow freeze on for floating elements
1066
1067 /* Intermediaries*/
1068 int is, numvertices, numsegments;
1069 int *pairindices = NULL;
1070
1071 numvertices=element->GetNumberOfVertices();
1072 element->VerticalSegmentIndices(&pairindices,&numsegments);
1073
1074 IssmDouble* watercolumn = xNew<IssmDouble>(numvertices);
1075 IssmDouble* drainrate_column = xNew<IssmDouble>(numsegments);
1076 IssmDouble* drainrate_element = xNew<IssmDouble>(numsegments);
1077
1078 element->GetInputListOnVertices(watercolumn,WatercolumnEnum);
1079
1080 for(is=0;is<numsegments;is++) drainrate_column[is]=0.;
1081 Element* elementi = element;
1082 for(;;){
1083 for(is=0;is<numsegments;is++) drainrate_element[is]=0.;
1084 DrainWaterfraction(elementi,drainrate_element); // TODO: make sure every vertex is only drained once
1085 for(is=0;is<numsegments;is++) drainrate_column[is]+=drainrate_element[is];
1086
1087 if(elementi->IsOnSurface()) break;
1088 elementi=elementi->GetUpperElement();
1089 }
1090 /* add drained water to water column*/
1091 for(is=0;is<numsegments;is++) watercolumn[is]+=drainrate_column[is];
1092 /* Feed updated water column back into model */
1093 element->AddInput(WatercolumnEnum,watercolumn,P1Enum);
1094
1095 xDelete<int>(pairindices);
1096 xDelete<IssmDouble>(drainrate_column);
1097 xDelete<IssmDouble>(drainrate_element);
1098 xDelete<IssmDouble>(watercolumn);
1099}/*}}}*/
1100void EnthalpyAnalysis::DrainWaterfraction(Element* element, IssmDouble* pdrainrate_element){/*{{{*/
1101
1102 /* Check if ice in element */
1103 if(!element->IsIceInElement()) return;
1104
1105 /*Intermediaries*/
1106 int iv,is,vertexdown,vertexup,numsegments;
1107 IssmDouble dt, height_element;
1108 IssmDouble rho_water, rho_ice;
1109 int numvertices = element->GetNumberOfVertices();
1110
1111 IssmDouble* xyz_list = NULL;
1112 IssmDouble* enthalpies = xNew<IssmDouble>(numvertices);
1113 IssmDouble* pressures = xNew<IssmDouble>(numvertices);
1114 IssmDouble* temperatures = xNew<IssmDouble>(numvertices);
1115 IssmDouble* waterfractions = xNew<IssmDouble>(numvertices);
1116 IssmDouble* deltawaterfractions = xNew<IssmDouble>(numvertices);
1117 int *pairindices = NULL;
1118
1119 rho_ice=element->GetMaterialParameter(MaterialsRhoIceEnum);
1120 rho_water=element->GetMaterialParameter(MaterialsRhoSeawaterEnum);
1121
1122 element->GetVerticesCoordinates(&xyz_list);
1123 element->GetInputListOnVertices(enthalpies,EnthalpyEnum);
1124 element->GetInputListOnVertices(pressures,PressureEnum);
1125
1126 element->FindParam(&dt,TimesteppingTimeStepEnum);
1127 for(iv=0;iv<numvertices;iv++){
1128 element->EnthalpyToThermal(&temperatures[iv],&waterfractions[iv], enthalpies[iv],pressures[iv]);
1129 deltawaterfractions[iv]=DrainageFunctionWaterfraction(waterfractions[iv], dt);
1130 }
1131
1132 /*drain waterfraction, feed updated variables back into model*/
1133 for(iv=0;iv<numvertices;iv++){
1134 if(reCast<bool,IssmDouble>(dt))
1135 waterfractions[iv]-=deltawaterfractions[iv]*dt;
1136 else
1137 waterfractions[iv]-=deltawaterfractions[iv];
1138 element->ThermalToEnthalpy(&enthalpies[iv], temperatures[iv], waterfractions[iv], pressures[iv]);
1139 }
1140 element->AddInput(EnthalpyEnum,enthalpies,P1Enum);
1141 element->AddInput(WaterfractionEnum,waterfractions,P1Enum);
1142
1143 /*return meltwater column equivalent to drained water*/
1144 element->VerticalSegmentIndices(&pairindices,&numsegments);
1145 for(is=0;is<numsegments;is++){
1146 vertexdown = pairindices[is*2+0];
1147 vertexup = pairindices[is*2+1];
1148 height_element=fabs(xyz_list[vertexup*3+2]-xyz_list[vertexdown*3+2]);
1149 pdrainrate_element[is]=(deltawaterfractions[vertexdown]+deltawaterfractions[vertexup])/2.*height_element; // return water equivalent of drainage
1150 _assert_(pdrainrate_element[is]>=0.);
1151 }
1152
1153 /*Clean up and return*/
1154 xDelete<int>(pairindices);
1155 xDelete<IssmDouble>(xyz_list);
1156 xDelete<IssmDouble>(enthalpies);
1157 xDelete<IssmDouble>(pressures);
1158 xDelete<IssmDouble>(temperatures);
1159 xDelete<IssmDouble>(waterfractions);
1160 xDelete<IssmDouble>(deltawaterfractions);
1161}/*}}}*/
1162void EnthalpyAnalysis::UpdateBasalConstraints(Element* element){/*{{{*/
1163
1164 /* Check if ice in element */
1165 if(!element->IsIceInElement()) return;
1166
1167 /* Only update Constraints at the base of grounded ice*/
1168 if(!(element->IsOnBase()) || element->IsFloating()) return;
1169
1170 /*Intermediary*/
1171 bool isdynamicbasalspc,setspc;
1172 int numindices, numindicesup;
1173 IssmDouble pressure, pressureup;
1174 IssmDouble h_pmp, enthalpy, enthalpyup;
1175 IssmDouble watercolumn;
1176 int *indices = NULL, *indicesup = NULL;
1177 Node* node = NULL;
1178
1179 /*Check wether dynamic basal boundary conditions are activated */
1180 element->FindParam(&isdynamicbasalspc,ThermalIsdynamicbasalspcEnum);
1181 if(!isdynamicbasalspc) return;
1182
1183 /*Fetch indices of basal & surface nodes for this finite element*/
1184 Penta *penta = (Penta *) element; // TODO: add Basal-/SurfaceNodeIndices to element.h, and change this to Element*
1185 penta->BasalNodeIndices(&numindices,&indices,element->GetElementType());
1186 penta->SurfaceNodeIndices(&numindicesup,&indicesup,element->GetElementType());
1187 _assert_(numindices==numindicesup);
1188
1189 /*Get parameters and inputs: */
1190 Input* pressure_input=element->GetInput(PressureEnum); _assert_(pressure_input);
1191 Input* enthalpy_input=element->GetInput(EnthalpyEnum); _assert_(enthalpy_input);
1192 Input* watercolumn_input=element->GetInput(WatercolumnEnum); _assert_(watercolumn_input);
1193
1194 /*if there is a temperate layer of zero thickness, set spc enthalpy=h_pmp at that node*/
1195 GaussPenta* gauss=new GaussPenta();
1196 GaussPenta* gaussup=new GaussPenta();
1197 for(int i=0;i<numindices;i++){
1198 gauss->GaussNode(element->GetElementType(),indices[i]);
1199 gaussup->GaussNode(element->GetElementType(),indicesup[i]);
1200
1201 /*Check wether there is a temperate layer at the base or not */
1202 /*check if node is temperate, else continue*/
1203 enthalpy_input->GetInputValue(&enthalpy, gauss);
1204 pressure_input->GetInputValue(&pressure, gauss);
1205 watercolumn_input->GetInputValue(&watercolumn,gauss);
1206 h_pmp=PureIceEnthalpy(element,pressure);
1207 if (enthalpy>=h_pmp){
1208 /*check if upper node is temperate, too.
1209 if yes, then we have a temperate layer of positive thickness and reset the spc.
1210 if not, apply dirichlet BC.*/
1211 enthalpy_input->GetInputValue(&enthalpyup, gaussup);
1212 pressure_input->GetInputValue(&pressureup, gaussup);
1213 setspc=((enthalpyup<PureIceEnthalpy(element,pressureup)) && (watercolumn>=0.))?true:false;
1214 }
1215 else
1216 setspc = false;
1217
1218 node=element->GetNode(indices[i]);
1219 if(setspc)
1220 node->ApplyConstraint(0,h_pmp); /*apply spc*/
1221 else
1222 node->DofInFSet(0); /*remove spc*/
1223 }
1224
1225 /*Free ressources:*/
1226 xDelete<int>(indices);
1227 xDelete<int>(indicesup);
1228 delete gauss;
1229 delete gaussup;
1230}/*}}}*/
1231
1232/*Intermediaries*/
1233IssmDouble EnthalpyAnalysis::EnthalpyDiffusionParameter(Element* element,IssmDouble enthalpy,IssmDouble pressure){/*{{{*/
1234
1235 IssmDouble heatcapacity = element->GetMaterialParameter(MaterialsHeatcapacityEnum);
1236 IssmDouble temperateiceconductivity = element->GetMaterialParameter(MaterialsTemperateiceconductivityEnum);
1237 IssmDouble thermalconductivity = element->GetMaterialParameter(MaterialsThermalconductivityEnum);
1238
1239 if(enthalpy < PureIceEnthalpy(element,pressure)){
1240 return thermalconductivity/heatcapacity;
1241 }
1242 else{
1243 return temperateiceconductivity/heatcapacity;
1244 }
1245}/*}}}*/
1246IssmDouble EnthalpyAnalysis::EnthalpyDiffusionParameterVolume(Element* element,int enthalpy_enum){/*{{{*/
1247
1248 int iv;
1249 IssmDouble lambda; /* fraction of cold ice */
1250 IssmDouble kappa,kappa_c,kappa_t; /* enthalpy conductivities */
1251 IssmDouble Hc,Ht;
1252
1253 /*Get pressures and enthalpies on vertices*/
1254 int numvertices = element->GetNumberOfVertices();
1255 IssmDouble* pressures = xNew<IssmDouble>(numvertices);
1256 IssmDouble* enthalpies = xNew<IssmDouble>(numvertices);
1257 IssmDouble* PIE = xNew<IssmDouble>(numvertices);
1258 IssmDouble* dHpmp = xNew<IssmDouble>(numvertices);
1259 element->GetInputListOnVertices(pressures,PressureEnum);
1260 element->GetInputListOnVertices(enthalpies,enthalpy_enum);
1261 for(iv=0;iv<numvertices;iv++){
1262 PIE[iv] = PureIceEnthalpy(element,pressures[iv]);
1263 dHpmp[iv] = enthalpies[iv]-PIE[iv];
1264 }
1265
1266 bool allequalsign = true;
1267 if(dHpmp[0]<0.){
1268 for(iv=1; iv<numvertices;iv++) allequalsign=(allequalsign && (dHpmp[iv]<0.));
1269 }
1270 else{
1271 for(iv=1; iv<numvertices;iv++) allequalsign=(allequalsign && (dHpmp[iv]>=0.));
1272 }
1273
1274 if(allequalsign){
1275 kappa = EnthalpyDiffusionParameter(element,enthalpies[0],pressures[0]);
1276 }
1277 else{
1278 /* return harmonic mean of thermal conductivities, weighted by fraction of cold/temperate ice,
1279 cf Patankar 1980, pp44 */
1280 kappa_c = EnthalpyDiffusionParameter(element,PureIceEnthalpy(element,0.)-1.,0.);
1281 kappa_t = EnthalpyDiffusionParameter(element,PureIceEnthalpy(element,0.)+1.,0.);
1282 Hc=0.; Ht=0.;
1283 for(iv=0; iv<numvertices;iv++){
1284 if(enthalpies[iv]<PIE[iv])
1285 Hc+=(PIE[iv]-enthalpies[iv]);
1286 else
1287 Ht+=(enthalpies[iv]-PIE[iv]);
1288 }
1289 _assert_((Hc+Ht)>0.);
1290 lambda = Hc/(Hc+Ht);
1291 kappa = kappa_c*kappa_t/(lambda*kappa_t+(1.-lambda)*kappa_c); // ==(lambda/kappa_c + (1.-lambda)/kappa_t)^-1
1292 }
1293
1294 /*Clean up and return*/
1295 xDelete<IssmDouble>(PIE);
1296 xDelete<IssmDouble>(dHpmp);
1297 xDelete<IssmDouble>(pressures);
1298 xDelete<IssmDouble>(enthalpies);
1299 return kappa;
1300}/*}}}*/
1301IssmDouble EnthalpyAnalysis::PureIceEnthalpy(Element* element,IssmDouble pressure){/*{{{*/
1302
1303 IssmDouble heatcapacity = element->GetMaterialParameter(MaterialsHeatcapacityEnum);
1304 IssmDouble referencetemperature = element->GetMaterialParameter(ConstantsReferencetemperatureEnum);
1305
1306 return heatcapacity*(TMeltingPoint(element,pressure)-referencetemperature);
1307}/*}}}*/
1308IssmDouble EnthalpyAnalysis::TMeltingPoint(Element* element,IssmDouble pressure){/*{{{*/
1309
1310 IssmDouble meltingpoint = element->GetMaterialParameter(MaterialsMeltingpointEnum);
1311 IssmDouble beta = element->GetMaterialParameter(MaterialsBetaEnum);
1312
1313 return meltingpoint-beta*pressure;
1314}/*}}}*/
Note: See TracBrowser for help on using the repository browser.