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

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

BUG: fixing enthalpy KMatrix
File size: 32.0 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*/
[16539]8int EnthalpyAnalysis::DofsPerNode(int** doflist,int meshtype,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));
19
20 iomodel->FetchData(&requestedoutputs,&numoutputs,ThermalRequestedOutputsEnum);
21 parameters->AddObject(new IntParam(ThermalNumRequestedOutputsEnum,numoutputs));
22 if(numoutputs)parameters->AddObject(new StringArrayParam(ThermalRequestedOutputsEnum,requestedoutputs,numoutputs));
23 iomodel->DeleteData(&requestedoutputs,numoutputs,ThermalRequestedOutputsEnum);
[16539]24}/*}}}*/
25void EnthalpyAnalysis::UpdateElements(Elements* elements,IoModel* iomodel,int analysis_counter,int analysis_type){/*{{{*/
26
27 bool dakota_analysis;
28 bool isenthalpy;
29
30 /*Now, is the model 3d? otherwise, do nothing: */
31 if(iomodel->meshtype==Mesh2DhorizontalEnum)return;
32
33 /*Is enthalpy requested?*/
34 iomodel->Constant(&isenthalpy,ThermalIsenthalpyEnum);
35 if(!isenthalpy) return;
36
37 /*Fetch data needed: */
38 iomodel->FetchData(3,TemperatureEnum,WaterfractionEnum,PressureEnum);
39
40 /*Update elements: */
41 int counter=0;
42 for(int i=0;i<iomodel->numberofelements;i++){
43 if(iomodel->my_elements[i]){
44 Element* element=(Element*)elements->GetObjectByOffset(counter);
45 element->Update(i,iomodel,analysis_counter,analysis_type,P1Enum);
46 counter++;
47 }
48 }
49
50 iomodel->Constant(&dakota_analysis,QmuIsdakotaEnum);
51
52 iomodel->FetchDataToInput(elements,ThicknessEnum);
53 iomodel->FetchDataToInput(elements,SurfaceEnum);
54 iomodel->FetchDataToInput(elements,BedEnum);
55 iomodel->FetchDataToInput(elements,FrictionCoefficientEnum);
56 iomodel->FetchDataToInput(elements,FrictionPEnum);
57 iomodel->FetchDataToInput(elements,FrictionQEnum);
58 iomodel->FetchDataToInput(elements,MaskIceLevelsetEnum);
59 iomodel->FetchDataToInput(elements,MaskGroundediceLevelsetEnum);
60 iomodel->FetchDataToInput(elements,MeshElementonbedEnum);
61 iomodel->FetchDataToInput(elements,MeshElementonsurfaceEnum);
62 iomodel->FetchDataToInput(elements,FlowequationElementEquationEnum);
63 iomodel->FetchDataToInput(elements,MaterialsRheologyBEnum);
64 iomodel->FetchDataToInput(elements,MaterialsRheologyNEnum);
65 iomodel->FetchDataToInput(elements,PressureEnum);
66 iomodel->FetchDataToInput(elements,TemperatureEnum);
67 iomodel->FetchDataToInput(elements,WaterfractionEnum);
68 iomodel->FetchDataToInput(elements,EnthalpyEnum);
69 iomodel->FetchDataToInput(elements,BasalforcingsGeothermalfluxEnum);
70 iomodel->FetchDataToInput(elements,WatercolumnEnum);
71 iomodel->FetchDataToInput(elements,BasalforcingsMeltingRateEnum);
72 iomodel->FetchDataToInput(elements,VxEnum);
73 iomodel->FetchDataToInput(elements,VyEnum);
74 iomodel->FetchDataToInput(elements,VzEnum);
75 InputUpdateFromConstantx(elements,0.,VxMeshEnum);
76 InputUpdateFromConstantx(elements,0.,VyMeshEnum);
77 InputUpdateFromConstantx(elements,0.,VzMeshEnum);
78 if(dakota_analysis){
79 elements->InputDuplicate(TemperatureEnum,QmuTemperatureEnum);
80 elements->InputDuplicate(BasalforcingsMeltingRateEnum,QmuMeltingEnum);
81 elements->InputDuplicate(VxMeshEnum,QmuVxMeshEnum);
82 elements->InputDuplicate(VxMeshEnum,QmuVyMeshEnum);
83 elements->InputDuplicate(VxMeshEnum,QmuVzMeshEnum);
84 }
85
86 /*Free data: */
87 iomodel->DeleteData(3,TemperatureEnum,WaterfractionEnum,PressureEnum);
88}/*}}}*/
[16542]89void EnthalpyAnalysis::CreateNodes(Nodes* nodes,IoModel* iomodel){/*{{{*/
[16539]90
91 if(iomodel->meshtype==Mesh3DEnum) iomodel->FetchData(2,MeshVertexonbedEnum,MeshVertexonsurfaceEnum);
[16542]92 ::CreateNodes(nodes,iomodel,EnthalpyAnalysisEnum,P1Enum);
[16539]93 iomodel->DeleteData(2,MeshVertexonbedEnum,MeshVertexonsurfaceEnum);
94}/*}}}*/
[16542]95void EnthalpyAnalysis::CreateConstraints(Constraints* constraints,IoModel* iomodel){/*{{{*/
[16539]96
97 /*Intermediary*/
98 int count;
99 int M,N;
100 bool spcpresent = false;
101 IssmDouble heatcapacity;
102 IssmDouble referencetemperature;
103
104 /*Output*/
105 IssmDouble *spcvector = NULL;
106 IssmDouble* times=NULL;
107 IssmDouble* values=NULL;
108
109 /*Fetch parameters: */
110 iomodel->Constant(&heatcapacity,MaterialsHeatcapacityEnum);
111 iomodel->Constant(&referencetemperature,ConstantsReferencetemperatureEnum);
112
113 /*return if 2d mesh*/
114 if(iomodel->meshtype==Mesh2DhorizontalEnum) return;
115
116 /*Fetch data: */
117 iomodel->FetchData(&spcvector,&M,&N,ThermalSpctemperatureEnum);
118
119 //FIX ME: SHOULD USE IOMODELCREATECONSTRAINTS
120 /*Transient or static?:*/
121 if(M==iomodel->numberofvertices){
122 /*static: just create Constraints objects*/
123 count=0;
124
125 for(int i=0;i<iomodel->numberofvertices;i++){
126 /*keep only this partition's nodes:*/
127 if((iomodel->my_vertices[i])){
128
129 if (!xIsNan<IssmDouble>(spcvector[i])){
130
131 constraints->AddObject(new SpcStatic(iomodel->constraintcounter+count+1,iomodel->nodecounter+i+1,1,heatcapacity*(spcvector[i]-referencetemperature),EnthalpyAnalysisEnum));
132 count++;
133
134 }
135 }
136 }
137 }
138 else if (M==(iomodel->numberofvertices+1)){
139 /*transient: create transient SpcTransient objects. Same logic, except we need to retrieve
140 * various times and values to initialize an SpcTransient object: */
141 count=0;
142
143 /*figure out times: */
144 times=xNew<IssmDouble>(N);
145 for(int j=0;j<N;j++){
146 times[j]=spcvector[(M-1)*N+j];
147 }
148
149 /*Create constraints from x,y,z: */
150 for(int i=0;i<iomodel->numberofvertices;i++){
151
152 /*keep only this partition's nodes:*/
153 if((iomodel->my_vertices[i])){
154
155 /*figure out times and values: */
156 values=xNew<IssmDouble>(N);
157 spcpresent=false;
158 for(int j=0;j<N;j++){
159 values[j]=heatcapacity*(spcvector[i*N+j]-referencetemperature);
160 if(!xIsNan<IssmDouble>(values[j]))spcpresent=true; //NaN means no spc by default
161 }
162
163 if(spcpresent){
164 constraints->AddObject(new SpcTransient(iomodel->constraintcounter+count+1,iomodel->nodecounter+i+1,1,N,times,values,EnthalpyAnalysisEnum));
165 count++;
166 }
167 xDelete<IssmDouble>(values);
168 }
169 }
170 }
171 else{
172 _error_("Size of field " << EnumToStringx(ThermalSpctemperatureEnum) << " not supported");
173 }
174
175 /*Free ressources:*/
176 iomodel->DeleteData(spcvector,ThermalSpctemperatureEnum);
177 xDelete<IssmDouble>(times);
178 xDelete<IssmDouble>(values);
179}/*}}}*/
[16542]180void EnthalpyAnalysis::CreateLoads(Loads* loads, IoModel* iomodel){/*{{{*/
[16539]181
182 /*No loads */
183}/*}}}*/
[16675]184
[16782]185/*Finite Element Analysis*/
186ElementMatrix* EnthalpyAnalysis::CreateKMatrix(Element* element){/*{{{*/
[16888]187
188 /*compute all stiffness matrices for this element*/
189 ElementMatrix* Ke1=CreateKMatrixVolume(element);
190 ElementMatrix* Ke2=CreateKMatrixShelf(element);
191 ElementMatrix* Ke =new ElementMatrix(Ke1,Ke2);
192
193 /*clean-up and return*/
194 delete Ke1;
195 delete Ke2;
196 return Ke;
[16782]197}/*}}}*/
[16888]198ElementMatrix* EnthalpyAnalysis::CreateKMatrixVolume(Element* element){/*{{{*/
199
200 /*Intermediaries */
201 int stabilization;
202 IssmDouble Jdet,dt,u,v,w,um,vm,wm,vel;
203 IssmDouble h,hx,hy,hz,vx,vy,vz;
204 IssmDouble tau_parameter,diameter;
205 IssmDouble D_scalar;
206 IssmDouble* xyz_list = NULL;
207
208 /*Fetch number of nodes and dof for this finite element*/
209 int numnodes = element->GetNumberOfNodes();
210
211 /*Initialize Element vector and other vectors*/
212 ElementMatrix* Ke = element->NewElementMatrix();
213 IssmDouble* basis = xNew<IssmDouble>(numnodes);
214 IssmDouble* dbasis = xNew<IssmDouble>(3*numnodes);
215 IssmDouble* B = xNew<IssmDouble>(3*numnodes);
216 IssmDouble* Bprime = xNew<IssmDouble>(3*numnodes);
217 IssmDouble D[3][3] = {0.};
218 IssmDouble K[3][3];
219
220 /*Retrieve all inputs and parameters*/
221 element->GetVerticesCoordinates(&xyz_list);
222 element->FindParam(&dt,TimesteppingTimeStepEnum);
223 element->FindParam(&stabilization,ThermalStabilizationEnum);
224 IssmDouble rho_water = element->GetMaterialParameter(MaterialsRhoWaterEnum);
225 IssmDouble rho_ice = element->GetMaterialParameter(MaterialsRhoIceEnum);
226 IssmDouble gravity = element->GetMaterialParameter(ConstantsGEnum);
227 IssmDouble heatcapacity = element->GetMaterialParameter(MaterialsHeatcapacityEnum);
228 IssmDouble thermalconductivity = element->GetMaterialParameter(MaterialsThermalconductivityEnum);
229 Input* vx_input = element->GetInput(VxEnum); _assert_(vx_input);
230 Input* vy_input = element->GetInput(VyEnum); _assert_(vy_input);
231 Input* vz_input = element->GetInput(VzEnum); _assert_(vz_input);
232 Input* vxm_input = element->GetInput(VxMeshEnum); _assert_(vxm_input);
233 Input* vym_input = element->GetInput(VyMeshEnum); _assert_(vym_input);
234 Input* vzm_input = element->GetInput(VzMeshEnum); _assert_(vzm_input);
235 if(stabilization==2) diameter=element->MinEdgeLength(xyz_list);
236
237 /*Enthalpy diffusion parameter*/
[16895]238 IssmDouble kappa=this->EnthalpyDiffusionParameterVolume(element,EnthalpyEnum); _assert_(kappa>0.);
[16888]239
240 /* Start looping on the number of gaussian points: */
241 Gauss* gauss=element->NewGauss(2);
242 for(int ig=gauss->begin();ig<gauss->end();ig++){
243 gauss->GaussPoint(ig);
244
245 element->JacobianDeterminant(&Jdet,xyz_list,gauss);
246 D_scalar=gauss->weight*Jdet;
247 if(dt!=0.) D_scalar=D_scalar*dt;
248
249 /*Conduction: */
250 GetBConduct(B,element,xyz_list,gauss);
251 D[0][0]=D_scalar*kappa/rho_ice;
252 D[1][1]=D_scalar*kappa/rho_ice;
253 D[2][2]=D_scalar*kappa/rho_ice;
254 TripleMultiply(B,3,numnodes,1,
255 &D[0][0],3,3,0,
256 B,3,numnodes,0,
257 &Ke->values[0],1);
258
259 /*Advection: */
260 GetBAdvec(B,element,xyz_list,gauss);
261 GetBAdvecprime(Bprime,element,xyz_list,gauss);
262 vx_input->GetInputValue(&u,gauss); vxm_input->GetInputValue(&um,gauss); vx=u-um;
263 vy_input->GetInputValue(&v,gauss); vym_input->GetInputValue(&vm,gauss); vy=v-vm;
264 vz_input->GetInputValue(&w,gauss); vzm_input->GetInputValue(&wm,gauss); vz=w-wm;
265 D[0][0]=D_scalar*vx;
266 D[1][1]=D_scalar*vy;
267 D[2][2]=D_scalar*vz;
268 TripleMultiply(B,3,numnodes,1,
269 &D[0][0],3,3,0,
270 Bprime,3,numnodes,0,
271 &Ke->values[0],1);
272
273 /*Transient: */
274 if(dt!=0.){
275 D_scalar=gauss->weight*Jdet;
276 element->NodalFunctions(basis,gauss);
277 TripleMultiply(basis,numnodes,1,0,
278 &D_scalar,1,1,0,
279 basis,1,numnodes,0,
280 &Ke->values[0],1);
281 D_scalar=D_scalar*dt;
282 }
283
284 /*Artifficial diffusivity*/
285 if(stabilization==1){
286 element->ElementSizes(&hx,&hy,&hz);
287 vel=sqrt(vx*vx + vy*vy + vz*vz)+1.e-14;
288 h=sqrt( pow(hx*vx/vel,2) + pow(hy*vy/vel,2) + pow(hz*vz/vel,2));
[16894]289 K[0][0]=h/(2.*vel)*vx*vx; K[0][1]=h/(2.*vel)*vx*vy; K[0][2]=h/(2.*vel)*vx*vz;
290 K[1][0]=h/(2.*vel)*vy*vx; K[1][1]=h/(2.*vel)*vy*vy; K[1][2]=h/(2.*vel)*vy*vz;
291 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]292 for(int i=0;i<3;i++) for(int j=0;j<3;j++) K[i][j] = D_scalar*K[i][j];
293
294 GetBAdvecprime(Bprime,element,xyz_list,gauss);
295 TripleMultiply(Bprime,3,numnodes,1,
296 &K[0][0],3,3,0,
297 Bprime,3,numnodes,0,
298 &Ke->values[0],1);
299 }
300 else if(stabilization==2){
301 element->NodalFunctionsDerivatives(dbasis,xyz_list,gauss);
302 tau_parameter=element->StabilizationParameter(u-um,v-vm,w-wm,diameter,kappa/rho_ice);
303 for(int i=0;i<numnodes;i++){
304 for(int j=0;j<numnodes;j++){
305 Ke->values[i*numnodes+j]+=tau_parameter*D_scalar*
[16895]306 ((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]307 }
308 }
309 if(dt!=0.){
[16896]310 D_scalar=gauss->weight*Jdet;
[16888]311 for(int i=0;i<numnodes;i++){
312 for(int j=0;j<numnodes;j++){
[16895]313 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]314 }
315 }
316 }
317 }
318 }
319
320 /*Clean up and return*/
321 xDelete<IssmDouble>(xyz_list);
322 xDelete<IssmDouble>(basis);
323 xDelete<IssmDouble>(dbasis);
324 xDelete<IssmDouble>(B);
325 xDelete<IssmDouble>(Bprime);
326 delete gauss;
327 return Ke;
328}/*}}}*/
329ElementMatrix* EnthalpyAnalysis::CreateKMatrixShelf(Element* element){/*{{{*/
330
331 /*Initialize Element matrix and return if necessary*/
332 if(!element->IsOnBed() || !element->IsFloating()) return NULL;
333
334 IssmDouble dt,Jdet,D;
335 IssmDouble *xyz_list_base = NULL;
336
337 /*Get basal element*/
338 if(!element->IsOnBed() || !element->IsFloating()) return NULL;
339
340 /*Fetch number of nodes for this finite element*/
341 int numnodes = element->GetNumberOfNodes();
342
343 /*Initialize vectors*/
344 ElementMatrix* Ke = element->NewElementMatrix();
345 IssmDouble* basis = xNew<IssmDouble>(numnodes);
346
347 /*Retrieve all inputs and parameters*/
348 element->GetVerticesCoordinatesBase(&xyz_list_base);
349 element->FindParam(&dt,TimesteppingTimeStepEnum);
350 IssmDouble gravity = element->GetMaterialParameter(ConstantsGEnum);
351 IssmDouble rho_water = element->GetMaterialParameter(MaterialsRhoWaterEnum);
352 IssmDouble rho_ice = element->GetMaterialParameter(MaterialsRhoIceEnum);
353 IssmDouble heatcapacity = element->GetMaterialParameter(MaterialsHeatcapacityEnum);
354 IssmDouble mixed_layer_capacity= element->GetMaterialParameter(MaterialsMixedLayerCapacityEnum);
355 IssmDouble thermal_exchange_vel= element->GetMaterialParameter(MaterialsThermalExchangeVelocityEnum);
356
357 /* Start looping on the number of gaussian points: */
358 Gauss* gauss=element->NewGaussBase(2);
359 for(int ig=gauss->begin();ig<gauss->end();ig++){
360 gauss->GaussPoint(ig);
361
362 element->JacobianDeterminantBase(&Jdet,xyz_list_base,gauss);
363 element->NodalFunctions(basis,gauss);
364
365 D=gauss->weight*Jdet*rho_water*mixed_layer_capacity*thermal_exchange_vel/(heatcapacity*rho_ice);
366 if(reCast<bool,IssmDouble>(dt)) D=dt*D;
367 TripleMultiply(basis,numnodes,1,0,
368 &D,1,1,0,
369 basis,1,numnodes,0,
370 &Ke->values[0],1);
371
372 }
373
374 /*Clean up and return*/
375 delete gauss;
376 xDelete<IssmDouble>(basis);
377 xDelete<IssmDouble>(xyz_list_base);
378 return Ke;
379}/*}}}*/
[16782]380ElementVector* EnthalpyAnalysis::CreatePVector(Element* element){/*{{{*/
[16812]381
382 /*compute all load vectors for this element*/
383 ElementVector* pe1=CreatePVectorVolume(element);
384 ElementVector* pe2=CreatePVectorSheet(element);
385 ElementVector* pe3=CreatePVectorShelf(element);
386 ElementVector* pe =new ElementVector(pe1,pe2,pe3);
387
388 /*clean-up and return*/
389 delete pe1;
390 delete pe2;
391 delete pe3;
392 return pe;
[16782]393}/*}}}*/
[16812]394ElementVector* EnthalpyAnalysis::CreatePVectorVolume(Element* element){/*{{{*/
395
396 /*Intermediaries*/
397 int stabilization;
398 IssmDouble Jdet,phi,dt;
399 IssmDouble enthalpy;
400 IssmDouble kappa,tau_parameter,diameter;
401 IssmDouble u,v,w;
402 IssmDouble scalar_def,scalar_transient;
403 IssmDouble* xyz_list = NULL;
404
405 /*Fetch number of nodes and dof for this finite element*/
406 int numnodes = element->GetNumberOfNodes();
407 int numvertices = element->GetNumberOfVertices();
408
409 /*Initialize Element vector*/
410 ElementVector* pe = element->NewElementVector();
411 IssmDouble* basis = xNew<IssmDouble>(numnodes);
412 IssmDouble* dbasis = xNew<IssmDouble>(3*numnodes);
413
414 /*Retrieve all inputs and parameters*/
415 element->GetVerticesCoordinates(&xyz_list);
416 IssmDouble rho_ice = element->GetMaterialParameter(MaterialsRhoIceEnum);
417 IssmDouble thermalconductivity = element->GetMaterialParameter(MaterialsThermalconductivityEnum);
418 element->FindParam(&dt,TimesteppingTimeStepEnum);
419 element->FindParam(&stabilization,ThermalStabilizationEnum);
420 Input* vx_input=element->GetInput(VxEnum); _assert_(vx_input);
421 Input* vy_input=element->GetInput(VyEnum); _assert_(vy_input);
422 Input* vz_input=element->GetInput(VzEnum); _assert_(vz_input);
423 Input* enthalpy_input = NULL;
424 if(reCast<bool,IssmDouble>(dt)){enthalpy_input = element->GetInput(EnthalpyEnum); _assert_(enthalpy_input);}
425 if(stabilization==2){
426 diameter=element->MinEdgeLength(xyz_list);
[16895]427 kappa=this->EnthalpyDiffusionParameterVolume(element,EnthalpyPicardEnum); _assert_(kappa>0.);
[16812]428 }
429
430 /* Start looping on the number of gaussian points: */
[16895]431 Gauss* gauss=element->NewGauss(2);
[16812]432 for(int ig=gauss->begin();ig<gauss->end();ig++){
433 gauss->GaussPoint(ig);
434
435 element->JacobianDeterminant(&Jdet,xyz_list,gauss);
436 element->NodalFunctions(basis,gauss);
437 element->ViscousHeating(&phi,xyz_list,gauss,vx_input,vy_input,vz_input);
438
439 scalar_def=phi/rho_ice*Jdet*gauss->weight;
[16895]440 if(dt!=0.) scalar_def=scalar_def*dt;
[16812]441
442 /*TODO: add -beta*laplace T_m(p)*/
443 for(int i=0;i<numnodes;i++) pe->values[i]+=scalar_def*basis[i];
444
445 /* Build transient now */
446 if(reCast<bool,IssmDouble>(dt)){
447 enthalpy_input->GetInputValue(&enthalpy, gauss);
448 scalar_transient=enthalpy*Jdet*gauss->weight;
449 for(int i=0;i<numnodes;i++) pe->values[i]+=scalar_transient*basis[i];
450 }
451
452 if(stabilization==2){
453 element->NodalFunctionsDerivatives(dbasis,xyz_list,gauss);
454
455 vx_input->GetInputValue(&u,gauss);
456 vy_input->GetInputValue(&v,gauss);
457 vz_input->GetInputValue(&w,gauss);
[16895]458 tau_parameter=element->StabilizationParameter(u,v,w,diameter,kappa/rho_ice);
[16812]459
[16895]460 for(int 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]);
461
462 if(dt!=0.){
463 for(int 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]464 }
465 }
466 }
467
468 /*Clean up and return*/
469 xDelete<IssmDouble>(basis);
470 xDelete<IssmDouble>(dbasis);
471 xDelete<IssmDouble>(xyz_list);
472 delete gauss;
473 return pe;
474
475}/*}}}*/
476ElementVector* EnthalpyAnalysis::CreatePVectorSheet(Element* element){/*{{{*/
[16888]477
478 /* Geothermal flux on ice sheet base and basal friction */
479 if(!element->IsOnBed() || element->IsFloating()) return NULL;
480
481 IssmDouble dt,Jdet,enthalpy,pressure,watercolumn,geothermalflux,vx,vy,vz;
482 IssmDouble enthalpyup,pressureup,alpha2,scalar,basalfriction,heatflux;
483 IssmDouble *xyz_list_base = NULL;
484
485 /*Fetch number of nodes for this finite element*/
486 int numnodes = element->GetNumberOfNodes();
487
488 /*Initialize vectors*/
489 ElementVector* pe = element->NewElementVector();
490 IssmDouble* basis = xNew<IssmDouble>(numnodes);
491
492 /*Retrieve all inputs and parameters*/
493 element->GetVerticesCoordinatesBase(&xyz_list_base);
494 element->FindParam(&dt,TimesteppingTimeStepEnum);
495 Input* vx_input = element->GetInput(VxEnum); _assert_(vx_input);
496 Input* vy_input = element->GetInput(VyEnum); _assert_(vy_input);
497 Input* vz_input = element->GetInput(VzEnum); _assert_(vz_input);
498 Input* enthalpy_input = element->GetInput(EnthalpyPicardEnum); _assert_(enthalpy_input);
499 Input* pressure_input = element->GetInput(PressureEnum); _assert_(pressure_input);
500 Input* geothermalflux_input = element->GetInput(BasalforcingsGeothermalfluxEnum); _assert_(geothermalflux_input);
501 Input* watercolumn_input = element->GetInput(WatercolumnEnum); _assert_(watercolumn_input);
502 IssmDouble rho_ice = element->GetMaterialParameter(MaterialsRhoIceEnum);
503 IssmDouble heatcapacity = element->GetMaterialParameter(MaterialsHeatcapacityEnum);
504
505 /*Build friction element, needed later: */
506 Friction* friction=new Friction(element,3);
507
508 /* Start looping on the number of gaussian points: */
509 Gauss* gauss = element->NewGaussBase(2);
510 Gauss* gaussup = element->NewGaussTop(2);
511 for(int ig=gauss->begin();ig<gauss->end();ig++){
512 gauss->GaussPoint(ig);
513
514 element->JacobianDeterminantBase(&Jdet,xyz_list_base,gauss);
515 element->NodalFunctions(basis,gauss);
516
517 enthalpy_input->GetInputValue(&enthalpy,gauss);
518 pressure_input->GetInputValue(&pressure,gauss);
519 watercolumn_input->GetInputValue(&watercolumn,gauss);
520
521 if((watercolumn<=0.) && (enthalpy < PureIceEnthalpy(element,pressure))){
522 /* the above check is equivalent to
523 NOT ((watercolumn>0.) AND (enthalpy<PIE)) AND (enthalpy<PIE)*/
524 geothermalflux_input->GetInputValue(&geothermalflux,gauss);
525
526 friction->GetAlpha2(&alpha2,gauss,vx_input,vy_input,vz_input);
527 vx_input->GetInputValue(&vx,gauss);
528 vy_input->GetInputValue(&vy,gauss);
529 vz_input->GetInputValue(&vz,gauss);
530 basalfriction = alpha2*(vx*vx + vy*vy + vz*vz);
531 heatflux = (basalfriction+geothermalflux)/(rho_ice);
532
533 scalar = gauss->weight*Jdet*heatflux;
534 if(dt!=0.) scalar=dt*scalar;
535
536 for(int i=0;i<numnodes;i++) pe->values[i]+=scalar*basis[i];
537 }
538 else if(enthalpy >= PureIceEnthalpy(element,pressure)){
539 /* check positive thickness of temperate basal ice layer */
540 enthalpy_input->GetInputValue(&enthalpyup,gaussup);
541 pressure_input->GetInputValue(&pressureup,gaussup);
542 if(enthalpyup >= PureIceEnthalpy(element,pressureup)){
543 // TODO: temperate ice has positive thickness: grad enthalpy*n=0.
544 }
545 else{
546 // only base temperate, set Dirichlet BCs in Penta::UpdateBasalConstraintsEnthalpy()
547 }
548 }
549 else{
550 // base cold, but watercolumn positive. Set base to h_pmp.
551 }
552 }
553
554 /*Clean up and return*/
555 delete gauss;
556 delete gaussup;
557 delete friction;
558 xDelete<IssmDouble>(basis);
559 xDelete<IssmDouble>(xyz_list_base);
560 return pe;
561
[16812]562}/*}}}*/
563ElementVector* EnthalpyAnalysis::CreatePVectorShelf(Element* element){/*{{{*/
564
[16888]565 /*Get basal element*/
566 if(!element->IsOnBed() || !element->IsFloating()) return NULL;
567
[16813]568 IssmDouble h_pmp,dt,Jdet,scalar_ocean,pressure;
[16812]569 IssmDouble *xyz_list_base = NULL;
570
571 /*Fetch number of nodes for this finite element*/
572 int numnodes = element->GetNumberOfNodes();
573
574 /*Initialize vectors*/
575 ElementVector* pe = element->NewElementVector();
576 IssmDouble* basis = xNew<IssmDouble>(numnodes);
577
578 /*Retrieve all inputs and parameters*/
579 element->GetVerticesCoordinatesBase(&xyz_list_base);
580 element->FindParam(&dt,TimesteppingTimeStepEnum);
581 Input* pressure_input=element->GetInput(PressureEnum); _assert_(pressure_input);
582 IssmDouble gravity = element->GetMaterialParameter(ConstantsGEnum);
583 IssmDouble rho_water = element->GetMaterialParameter(MaterialsRhoWaterEnum);
584 IssmDouble rho_ice = element->GetMaterialParameter(MaterialsRhoIceEnum);
585 IssmDouble heatcapacity = element->GetMaterialParameter(MaterialsHeatcapacityEnum);
586 IssmDouble mixed_layer_capacity= element->GetMaterialParameter(MaterialsMixedLayerCapacityEnum);
587 IssmDouble thermal_exchange_vel= element->GetMaterialParameter(MaterialsThermalExchangeVelocityEnum);
588
589 /* Start looping on the number of gaussian points: */
590 Gauss* gauss=element->NewGaussBase(2);
591 for(int ig=gauss->begin();ig<gauss->end();ig++){
592 gauss->GaussPoint(ig);
593
594 element->JacobianDeterminantBase(&Jdet,xyz_list_base,gauss);
595 element->NodalFunctions(basis,gauss);
596
597 pressure_input->GetInputValue(&pressure,gauss);
[16813]598 h_pmp=element->PureIceEnthalpy(pressure);
[16812]599
[16813]600 scalar_ocean=gauss->weight*Jdet*rho_water*mixed_layer_capacity*thermal_exchange_vel*h_pmp/(heatcapacity*rho_ice);
[16812]601 if(reCast<bool,IssmDouble>(dt)) scalar_ocean=dt*scalar_ocean;
602
603 for(int i=0;i<numnodes;i++) pe->values[i]+=scalar_ocean*basis[i];
604 }
605
606 /*Clean up and return*/
607 delete gauss;
608 xDelete<IssmDouble>(basis);
609 xDelete<IssmDouble>(xyz_list_base);
610 return pe;
611}/*}}}*/
[16888]612void EnthalpyAnalysis::GetBConduct(IssmDouble* B,Element* element,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/
613 /*Compute B matrix. B=[B1 B2 B3 B4 B5 B6] where Bi is of size 5*NDOF1.
614 * For node i, Bi' can be expressed in the actual coordinate system
615 * by:
616 * Bi_conduct=[ dh/dx ]
617 * [ dh/dy ]
618 * [ dh/dz ]
619 * where h is the interpolation function for node i.
620 *
621 * We assume B has been allocated already, of size: 3x(NDOF1*numnodes)
622 */
623
624 /*Fetch number of nodes for this finite element*/
625 int numnodes = element->GetNumberOfNodes();
626
627 /*Get nodal functions derivatives*/
628 IssmDouble* dbasis=xNew<IssmDouble>(3*numnodes);
629 element->NodalFunctionsDerivatives(dbasis,xyz_list,gauss);
630
631 /*Build B: */
632 for(int i=0;i<numnodes;i++){
633 B[numnodes*0+i] = dbasis[0*numnodes+i];
634 B[numnodes*1+i] = dbasis[1*numnodes+i];
635 B[numnodes*2+i] = dbasis[2*numnodes+i];
636 }
637
638 /*Clean-up*/
639 xDelete<IssmDouble>(dbasis);
640}/*}}}*/
641void EnthalpyAnalysis::GetBAdvec(IssmDouble* B,Element* element,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/
642 /*Compute B matrix. B=[B1 B2 B3 B4 B5 B6] where Bi is of size 5*NDOF1.
643 * For node i, Bi' can be expressed in the actual coordinate system
644 * by:
645 * Bi_advec =[ h ]
646 * [ h ]
647 * [ h ]
648 * where h is the interpolation function for node i.
649 *
650 * We assume B has been allocated already, of size: 3x(NDOF1*NUMNODESP1)
651 */
652
653 /*Fetch number of nodes for this finite element*/
654 int numnodes = element->GetNumberOfNodes();
655
656 /*Get nodal functions*/
657 IssmDouble* basis=xNew<IssmDouble>(numnodes);
658 element->NodalFunctions(basis,gauss);
659
660 /*Build B: */
661 for(int i=0;i<numnodes;i++){
662 B[numnodes*0+i] = basis[i];
663 B[numnodes*1+i] = basis[i];
664 B[numnodes*2+i] = basis[i];
665 }
666
667 /*Clean-up*/
668 xDelete<IssmDouble>(basis);
669}/*}}}*/
670void EnthalpyAnalysis::GetBAdvecprime(IssmDouble* B,Element* element,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/
671 /*Compute B matrix. B=[B1 B2 B3 B4 B5 B6] where Bi is of size 5*NDOF1.
672 * For node i, Bi' can be expressed in the actual coordinate system
673 * by:
674 * Biprime_advec=[ dh/dx ]
675 * [ dh/dy ]
676 * [ dh/dz ]
677 * where h is the interpolation function for node i.
678 *
679 * We assume B has been allocated already, of size: 3x(NDOF1*numnodes)
680 */
681
682 /*Fetch number of nodes for this finite element*/
683 int numnodes = element->GetNumberOfNodes();
684
685 /*Get nodal functions derivatives*/
686 IssmDouble* dbasis=xNew<IssmDouble>(3*numnodes);
687 element->NodalFunctionsDerivatives(dbasis,xyz_list,gauss);
688
689 /*Build B: */
690 for(int i=0;i<numnodes;i++){
691 B[numnodes*0+i] = dbasis[0*numnodes+i];
692 B[numnodes*1+i] = dbasis[1*numnodes+i];
693 B[numnodes*2+i] = dbasis[2*numnodes+i];
694 }
695
696 /*Clean-up*/
697 xDelete<IssmDouble>(dbasis);
698}/*}}}*/
[16675]699void EnthalpyAnalysis::GetSolutionFromInputs(Vector<IssmDouble>* solution,Element* element){/*{{{*/
700 element->GetSolutionFromInputsOneDof(solution,EnthalpyEnum);
701}/*}}}*/
[16684]702void EnthalpyAnalysis::InputUpdateFromSolution(IssmDouble* solution,Element* element){/*{{{*/
[16734]703
704 bool converged;
705 int i,rheology_law;
706 IssmDouble B_average,s_average,T_average=0.,P_average=0.;
707 int *doflist = NULL;
708 IssmDouble *xyz_list = NULL;
709
710 /*Fetch number of nodes and dof for this finite element*/
711 int numnodes = element->GetNumberOfNodes();
712
713 /*Fetch dof list and allocate solution vector*/
714 element->GetDofList(&doflist,NoneApproximationEnum,GsetEnum);
715 IssmDouble* values = xNew<IssmDouble>(numnodes);
716 IssmDouble* pressure = xNew<IssmDouble>(numnodes);
[16745]717 IssmDouble* surface = xNew<IssmDouble>(numnodes);
718 IssmDouble* B = xNew<IssmDouble>(numnodes);
[16734]719 IssmDouble* temperature = xNew<IssmDouble>(numnodes);
720 IssmDouble* waterfraction = xNew<IssmDouble>(numnodes);
721
722 /*Use the dof list to index into the solution vector: */
723 for(i=0;i<numnodes;i++){
724 values[i]=solution[doflist[i]];
725
726 /*Check solution*/
727 if(xIsNan<IssmDouble>(values[i])) _error_("NaN found in solution vector");
728 }
729
730 /*Get all inputs and parameters*/
731 element->GetInputValue(&converged,ConvergedEnum);
[16745]732 element->GetInputListOnNodes(&pressure[0],PressureEnum);
[16734]733 if(converged){
734 for(i=0;i<numnodes;i++){
735 element->EnthalpyToThermal(&temperature[i],&waterfraction[i],values[i],pressure[i]);
736 if(waterfraction[i]<0.) _error_("Negative water fraction found in solution vector");
[16837]737 if(waterfraction[i]>1.) _error_("Water fraction >1 found in solution vector");
[16734]738 }
739 element->AddInput(EnthalpyEnum,values,P1Enum);
740 element->AddInput(WaterfractionEnum,waterfraction,P1Enum);
741 element->AddInput(TemperatureEnum,temperature,P1Enum);
742
743 /*Update Rheology only if converged (we must make sure that the temperature is below melting point
744 * otherwise the rheology could be negative*/
745 element->FindParam(&rheology_law,MaterialsRheologyLawEnum);
[16745]746 element->GetInputListOnNodes(&surface[0],SurfaceEnum);
[16734]747 switch(rheology_law){
748 case NoneEnum:
749 /*Do nothing: B is not temperature dependent*/
750 break;
751 case PatersonEnum:
[16748]752 for(i=0;i<numnodes;i++) B[i]=Paterson(temperature[i]);
[16745]753 element->AddMaterialInput(MaterialsRheologyBEnum,&B[0],P1Enum);
[16734]754 break;
[16745]755 case ArrheniusEnum:
[16734]756 element->GetVerticesCoordinates(&xyz_list);
[16748]757 for(i=0;i<numnodes;i++) B[i]=Arrhenius(temperature[i],surface[i]-xyz_list[i*3+2],element->GetMaterialParameter(MaterialsRheologyNEnum));
[16745]758 element->AddMaterialInput(MaterialsRheologyBEnum,&B[0],P1Enum);
[16734]759 break;
760 case LliboutryDuvalEnum:
[16750]761 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));
[16745]762 element->AddMaterialInput(MaterialsRheologyBEnum,&B[0],P1Enum);
763 break;
764 default: _error_("Rheology law " << EnumToStringx(rheology_law) << " not supported yet");
[16734]765 }
766 }
767 else{
768 element->AddInput(EnthalpyPicardEnum,values,P1Enum);
769 }
770
771 /*Free ressources:*/
772 xDelete<IssmDouble>(values);
773 xDelete<IssmDouble>(pressure);
[16745]774 xDelete<IssmDouble>(surface);
775 xDelete<IssmDouble>(B);
[16734]776 xDelete<IssmDouble>(temperature);
777 xDelete<IssmDouble>(waterfraction);
778 xDelete<IssmDouble>(xyz_list);
779 xDelete<int>(doflist);
[16684]780}/*}}}*/
[16888]781
782/*Intermediaries*/
783IssmDouble EnthalpyAnalysis::EnthalpyDiffusionParameter(Element* element,IssmDouble enthalpy,IssmDouble pressure){/*{{{*/
784
785 IssmDouble heatcapacity = element->GetMaterialParameter(MaterialsHeatcapacityEnum);
786 IssmDouble temperateiceconductivity = element->GetMaterialParameter(MaterialsTemperateiceconductivityEnum);
787 IssmDouble thermalconductivity = element->GetMaterialParameter(MaterialsThermalconductivityEnum);
788
789 if(enthalpy < PureIceEnthalpy(element,pressure)){
790 return thermalconductivity/heatcapacity;
791 }
792 else{
793 return temperateiceconductivity/heatcapacity;
794 }
795}/*}}}*/
[16895]796IssmDouble EnthalpyAnalysis::EnthalpyDiffusionParameterVolume(Element* element,int enthalpy_enum){/*{{{*/
[16888]797
798 int iv;
799 IssmDouble lambda; /* fraction of cold ice */
800 IssmDouble kappa ,kappa_c,kappa_t; /* enthalpy conductivities */
801 IssmDouble Hc,Ht;
802
803
804 /*Get pressures and enthalpies on vertices*/
805 int numvertices = element->GetNumberOfVertices();
806 IssmDouble* pressures = xNew<IssmDouble>(numvertices);
807 IssmDouble* enthalpies = xNew<IssmDouble>(numvertices);
808 IssmDouble* PIE = xNew<IssmDouble>(numvertices);
809 IssmDouble* dHpmp = xNew<IssmDouble>(numvertices);
810 element->GetInputListOnVertices(pressures,PressureEnum);
[16895]811 element->GetInputListOnVertices(enthalpies,enthalpy_enum);
[16888]812 for(iv=0;iv<numvertices;iv++){
813 PIE[iv] = PureIceEnthalpy(element,pressures[iv]);
814 dHpmp[iv] = enthalpies[iv]-PIE[iv];
815 }
816
817 bool allequalsign = true;
818 if(dHpmp[0]<0.){
819 for(iv=1; iv<numvertices;iv++) allequalsign=(allequalsign && (dHpmp[iv]<0.));
820 }
821 else{
822 for(iv=1; iv<numvertices;iv++) allequalsign=(allequalsign && (dHpmp[iv]>=0.));
823 }
824
825 if(allequalsign){
826 kappa = EnthalpyDiffusionParameter(element,enthalpies[0],pressures[0]);
827 }
828 else{
829 /* return harmonic mean of thermal conductivities, weighted by fraction of cold/temperate ice,
830 cf Patankar 1980, pp44 */
831 kappa_c = EnthalpyDiffusionParameter(element,PureIceEnthalpy(element,0.)-1.,0.);
832 kappa_t = EnthalpyDiffusionParameter(element,PureIceEnthalpy(element,0.)+1.,0.);
833 Hc=0.; Ht=0.;
834 for(iv=0; iv<numvertices;iv++){
835 if(enthalpies[iv]<PIE[iv])
836 Hc+=(PIE[iv]-enthalpies[iv]);
837 else
838 Ht+=(enthalpies[iv]-PIE[iv]);
839 }
840 _assert_((Hc+Ht)>0.);
841 lambda = Hc/(Hc+Ht);
842 kappa = 1./(lambda/kappa_c + (1.-lambda)/kappa_t);
843 }
844
845 /*Clean up and return*/
846 xDelete<IssmDouble>(PIE);
847 xDelete<IssmDouble>(dHpmp);
848 xDelete<IssmDouble>(pressures);
849 xDelete<IssmDouble>(enthalpies);
850 return kappa;
851}
852/*}}}*/
853IssmDouble EnthalpyAnalysis::PureIceEnthalpy(Element* element,IssmDouble pressure){/*{{{*/
854
855 IssmDouble heatcapacity = element->GetMaterialParameter(MaterialsHeatcapacityEnum);
856 IssmDouble referencetemperature = element->GetMaterialParameter(ConstantsReferencetemperatureEnum);
857
858 return heatcapacity*(TMeltingPoint(element,pressure)-referencetemperature);
859}/*}}}*/
860IssmDouble EnthalpyAnalysis::TMeltingPoint(Element* element,IssmDouble pressure){/*{{{*/
861
862 IssmDouble meltingpoint = element->GetMaterialParameter(MaterialsMeltingpointEnum);
863 IssmDouble beta = element->GetMaterialParameter(MaterialsBetaEnum);
864
865 return meltingpoint-beta*pressure;
866}/*}}}*/
Note: See TracBrowser for help on using the repository browser.