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

Last change on this file since 18804 was 18804, checked in by Mathieu Morlighem, 10 years ago

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