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

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

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