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

Last change on this file since 18589 was 18589, checked in by jbondzio, 10 years ago

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