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

Last change on this file since 17952 was 17952, checked in by cborstad, 11 years ago

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