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

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

NEW: added CreateDVector in all analyses
File size: 32.2 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 meshtype,int approximation){/*{{{*/
9 return 1;
10}/*}}}*/
11void EnthalpyAnalysis::UpdateParameters(Parameters* parameters,IoModel* iomodel,int solution_enum,int analysis_enum){/*{{{*/
12
13 int numoutputs;
14 char** requestedoutputs = NULL;
15
16 parameters->AddObject(iomodel->CopyConstantObject(ThermalStabilizationEnum));
17 parameters->AddObject(iomodel->CopyConstantObject(ThermalIsenthalpyEnum));
18 parameters->AddObject(iomodel->CopyConstantObject(ThermalIsdynamicbasalspcEnum));
19
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);
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}/*}}}*/
89void EnthalpyAnalysis::CreateNodes(Nodes* nodes,IoModel* iomodel){/*{{{*/
90
91 if(iomodel->meshtype==Mesh3DEnum) iomodel->FetchData(2,MeshVertexonbedEnum,MeshVertexonsurfaceEnum);
92 ::CreateNodes(nodes,iomodel,EnthalpyAnalysisEnum,P1Enum);
93 iomodel->DeleteData(2,MeshVertexonbedEnum,MeshVertexonsurfaceEnum);
94}/*}}}*/
95void EnthalpyAnalysis::CreateConstraints(Constraints* constraints,IoModel* iomodel){/*{{{*/
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}/*}}}*/
180void EnthalpyAnalysis::CreateLoads(Loads* loads, IoModel* iomodel){/*{{{*/
181
182 /*No loads */
183}/*}}}*/
184
185/*Finite Element Analysis*/
186ElementVector* EnthalpyAnalysis::CreateDVector(Element* element){/*{{{*/
187 /*Default, return NULL*/
188 return NULL;
189}/*}}}*/
190ElementMatrix* EnthalpyAnalysis::CreateJacobianMatrix(Element* element){/*{{{*/
191_error_("Not implemented");
192}/*}}}*/
193ElementMatrix* EnthalpyAnalysis::CreateKMatrix(Element* element){/*{{{*/
194
195 /*compute all stiffness matrices for this element*/
196 ElementMatrix* Ke1=CreateKMatrixVolume(element);
197 ElementMatrix* Ke2=CreateKMatrixShelf(element);
198 ElementMatrix* Ke =new ElementMatrix(Ke1,Ke2);
199
200 /*clean-up and return*/
201 delete Ke1;
202 delete Ke2;
203 return Ke;
204}/*}}}*/
205ElementMatrix* EnthalpyAnalysis::CreateKMatrixVolume(Element* element){/*{{{*/
206
207 /*Intermediaries */
208 int stabilization;
209 IssmDouble Jdet,dt,u,v,w,um,vm,wm,vel;
210 IssmDouble h,hx,hy,hz,vx,vy,vz;
211 IssmDouble tau_parameter,diameter;
212 IssmDouble D_scalar;
213 IssmDouble* xyz_list = NULL;
214
215 /*Fetch number of nodes and dof for this finite element*/
216 int numnodes = element->GetNumberOfNodes();
217
218 /*Initialize Element vector and other vectors*/
219 ElementMatrix* Ke = element->NewElementMatrix();
220 IssmDouble* basis = xNew<IssmDouble>(numnodes);
221 IssmDouble* dbasis = xNew<IssmDouble>(3*numnodes);
222 IssmDouble* B = xNew<IssmDouble>(3*numnodes);
223 IssmDouble* Bprime = xNew<IssmDouble>(3*numnodes);
224 IssmDouble D[3][3] = {0.};
225 IssmDouble K[3][3];
226
227 /*Retrieve all inputs and parameters*/
228 element->GetVerticesCoordinates(&xyz_list);
229 element->FindParam(&dt,TimesteppingTimeStepEnum);
230 element->FindParam(&stabilization,ThermalStabilizationEnum);
231 IssmDouble rho_water = element->GetMaterialParameter(MaterialsRhoWaterEnum);
232 IssmDouble rho_ice = element->GetMaterialParameter(MaterialsRhoIceEnum);
233 IssmDouble gravity = element->GetMaterialParameter(ConstantsGEnum);
234 IssmDouble heatcapacity = element->GetMaterialParameter(MaterialsHeatcapacityEnum);
235 IssmDouble thermalconductivity = element->GetMaterialParameter(MaterialsThermalconductivityEnum);
236 Input* vx_input = element->GetInput(VxEnum); _assert_(vx_input);
237 Input* vy_input = element->GetInput(VyEnum); _assert_(vy_input);
238 Input* vz_input = element->GetInput(VzEnum); _assert_(vz_input);
239 Input* vxm_input = element->GetInput(VxMeshEnum); _assert_(vxm_input);
240 Input* vym_input = element->GetInput(VyMeshEnum); _assert_(vym_input);
241 Input* vzm_input = element->GetInput(VzMeshEnum); _assert_(vzm_input);
242 if(stabilization==2) diameter=element->MinEdgeLength(xyz_list);
243
244 /*Enthalpy diffusion parameter*/
245 IssmDouble kappa=this->EnthalpyDiffusionParameterVolume(element,EnthalpyPicardEnum); _assert_(kappa>0.);
246
247 /* Start looping on the number of gaussian points: */
248 Gauss* gauss=element->NewGauss(2);
249 for(int ig=gauss->begin();ig<gauss->end();ig++){
250 gauss->GaussPoint(ig);
251
252 element->JacobianDeterminant(&Jdet,xyz_list,gauss);
253 D_scalar=gauss->weight*Jdet;
254 if(dt!=0.) D_scalar=D_scalar*dt;
255
256 /*Conduction: */
257 GetBConduct(B,element,xyz_list,gauss);
258 D[0][0]=D_scalar*kappa/rho_ice;
259 D[1][1]=D_scalar*kappa/rho_ice;
260 D[2][2]=D_scalar*kappa/rho_ice;
261 TripleMultiply(B,3,numnodes,1,
262 &D[0][0],3,3,0,
263 B,3,numnodes,0,
264 &Ke->values[0],1);
265
266 /*Advection: */
267 GetBAdvec(B,element,xyz_list,gauss);
268 GetBAdvecprime(Bprime,element,xyz_list,gauss);
269 vx_input->GetInputValue(&u,gauss); vxm_input->GetInputValue(&um,gauss); vx=u-um;
270 vy_input->GetInputValue(&v,gauss); vym_input->GetInputValue(&vm,gauss); vy=v-vm;
271 vz_input->GetInputValue(&w,gauss); vzm_input->GetInputValue(&wm,gauss); vz=w-wm;
272 D[0][0]=D_scalar*vx;
273 D[1][1]=D_scalar*vy;
274 D[2][2]=D_scalar*vz;
275 TripleMultiply(B,3,numnodes,1,
276 &D[0][0],3,3,0,
277 Bprime,3,numnodes,0,
278 &Ke->values[0],1);
279
280 /*Transient: */
281 if(dt!=0.){
282 D_scalar=gauss->weight*Jdet;
283 element->NodalFunctions(basis,gauss);
284 TripleMultiply(basis,numnodes,1,0,
285 &D_scalar,1,1,0,
286 basis,1,numnodes,0,
287 &Ke->values[0],1);
288 D_scalar=D_scalar*dt;
289 }
290
291 /*Artifficial diffusivity*/
292 if(stabilization==1){
293 element->ElementSizes(&hx,&hy,&hz);
294 vel=sqrt(vx*vx + vy*vy + vz*vz)+1.e-14;
295 h=sqrt( pow(hx*vx/vel,2) + pow(hy*vy/vel,2) + pow(hz*vz/vel,2));
296 K[0][0]=h/(2.*vel)*vx*vx; K[0][1]=h/(2.*vel)*vx*vy; K[0][2]=h/(2.*vel)*vx*vz;
297 K[1][0]=h/(2.*vel)*vy*vx; K[1][1]=h/(2.*vel)*vy*vy; K[1][2]=h/(2.*vel)*vy*vz;
298 K[2][0]=h/(2.*vel)*vz*vx; K[2][1]=h/(2.*vel)*vz*vy; K[2][2]=h/(2.*vel)*vz*vz;
299 for(int i=0;i<3;i++) for(int j=0;j<3;j++) K[i][j] = D_scalar*K[i][j];
300
301 GetBAdvecprime(Bprime,element,xyz_list,gauss);
302 TripleMultiply(Bprime,3,numnodes,1,
303 &K[0][0],3,3,0,
304 Bprime,3,numnodes,0,
305 &Ke->values[0],1);
306 }
307 else if(stabilization==2){
308 element->NodalFunctionsDerivatives(dbasis,xyz_list,gauss);
309 tau_parameter=element->StabilizationParameter(u-um,v-vm,w-wm,diameter,kappa/rho_ice);
310 for(int i=0;i<numnodes;i++){
311 for(int j=0;j<numnodes;j++){
312 Ke->values[i*numnodes+j]+=tau_parameter*D_scalar*
313 ((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]);
314 }
315 }
316 if(dt!=0.){
317 D_scalar=gauss->weight*Jdet;
318 for(int i=0;i<numnodes;i++){
319 for(int j=0;j<numnodes;j++){
320 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]);
321 }
322 }
323 }
324 }
325 }
326
327 /*Clean up and return*/
328 xDelete<IssmDouble>(xyz_list);
329 xDelete<IssmDouble>(basis);
330 xDelete<IssmDouble>(dbasis);
331 xDelete<IssmDouble>(B);
332 xDelete<IssmDouble>(Bprime);
333 delete gauss;
334 return Ke;
335}/*}}}*/
336ElementMatrix* EnthalpyAnalysis::CreateKMatrixShelf(Element* element){/*{{{*/
337
338 /*Initialize Element matrix and return if necessary*/
339 if(!element->IsOnBed() || !element->IsFloating()) return NULL;
340
341 /*Intermediaries*/
342 IssmDouble dt,Jdet,D;
343 IssmDouble *xyz_list_base = NULL;
344
345 /*Fetch number of nodes for this finite element*/
346 int numnodes = element->GetNumberOfNodes();
347
348 /*Initialize vectors*/
349 ElementMatrix* Ke = element->NewElementMatrix();
350 IssmDouble* basis = xNew<IssmDouble>(numnodes);
351
352 /*Retrieve all inputs and parameters*/
353 element->GetVerticesCoordinatesBase(&xyz_list_base);
354 element->FindParam(&dt,TimesteppingTimeStepEnum);
355 IssmDouble gravity = element->GetMaterialParameter(ConstantsGEnum);
356 IssmDouble rho_water = element->GetMaterialParameter(MaterialsRhoWaterEnum);
357 IssmDouble rho_ice = element->GetMaterialParameter(MaterialsRhoIceEnum);
358 IssmDouble heatcapacity = element->GetMaterialParameter(MaterialsHeatcapacityEnum);
359 IssmDouble mixed_layer_capacity= element->GetMaterialParameter(MaterialsMixedLayerCapacityEnum);
360 IssmDouble thermal_exchange_vel= element->GetMaterialParameter(MaterialsThermalExchangeVelocityEnum);
361
362 /* Start looping on the number of gaussian points: */
363 Gauss* gauss=element->NewGaussBase(2);
364 for(int ig=gauss->begin();ig<gauss->end();ig++){
365 gauss->GaussPoint(ig);
366
367 element->JacobianDeterminantBase(&Jdet,xyz_list_base,gauss);
368 element->NodalFunctions(basis,gauss);
369
370 D=gauss->weight*Jdet*rho_water*mixed_layer_capacity*thermal_exchange_vel/(heatcapacity*rho_ice);
371 if(reCast<bool,IssmDouble>(dt)) D=dt*D;
372 TripleMultiply(basis,numnodes,1,0,
373 &D,1,1,0,
374 basis,1,numnodes,0,
375 &Ke->values[0],1);
376
377 }
378
379 /*Clean up and return*/
380 delete gauss;
381 xDelete<IssmDouble>(basis);
382 xDelete<IssmDouble>(xyz_list_base);
383 return Ke;
384}/*}}}*/
385ElementVector* EnthalpyAnalysis::CreatePVector(Element* element){/*{{{*/
386
387 /*compute all load vectors for this element*/
388 ElementVector* pe1=CreatePVectorVolume(element);
389 ElementVector* pe2=CreatePVectorSheet(element);
390 ElementVector* pe3=CreatePVectorShelf(element);
391 ElementVector* pe =new ElementVector(pe1,pe2,pe3);
392
393 /*clean-up and return*/
394 delete pe1;
395 delete pe2;
396 delete pe3;
397 return pe;
398}/*}}}*/
399ElementVector* EnthalpyAnalysis::CreatePVectorVolume(Element* element){/*{{{*/
400
401 /*Intermediaries*/
402 int stabilization;
403 IssmDouble Jdet,phi,dt;
404 IssmDouble enthalpy;
405 IssmDouble kappa,tau_parameter,diameter;
406 IssmDouble u,v,w;
407 IssmDouble scalar_def,scalar_transient;
408 IssmDouble* xyz_list = NULL;
409
410 /*Fetch number of nodes and dof for this finite element*/
411 int numnodes = element->GetNumberOfNodes();
412 int numvertices = element->GetNumberOfVertices();
413
414 /*Initialize Element vector*/
415 ElementVector* pe = element->NewElementVector();
416 IssmDouble* basis = xNew<IssmDouble>(numnodes);
417 IssmDouble* dbasis = xNew<IssmDouble>(3*numnodes);
418
419 /*Retrieve all inputs and parameters*/
420 element->GetVerticesCoordinates(&xyz_list);
421 IssmDouble rho_ice = element->GetMaterialParameter(MaterialsRhoIceEnum);
422 IssmDouble thermalconductivity = element->GetMaterialParameter(MaterialsThermalconductivityEnum);
423 element->FindParam(&dt,TimesteppingTimeStepEnum);
424 element->FindParam(&stabilization,ThermalStabilizationEnum);
425 Input* vx_input=element->GetInput(VxEnum); _assert_(vx_input);
426 Input* vy_input=element->GetInput(VyEnum); _assert_(vy_input);
427 Input* vz_input=element->GetInput(VzEnum); _assert_(vz_input);
428 Input* enthalpy_input = NULL;
429 if(reCast<bool,IssmDouble>(dt)){enthalpy_input = element->GetInput(EnthalpyEnum); _assert_(enthalpy_input);}
430 if(stabilization==2){
431 diameter=element->MinEdgeLength(xyz_list);
432 kappa=this->EnthalpyDiffusionParameterVolume(element,EnthalpyPicardEnum); _assert_(kappa>0.);
433 }
434
435 /* Start looping on the number of gaussian points: */
436 Gauss* gauss=element->NewGauss(3);
437 for(int ig=gauss->begin();ig<gauss->end();ig++){
438 gauss->GaussPoint(ig);
439
440 element->JacobianDeterminant(&Jdet,xyz_list,gauss);
441 element->NodalFunctions(basis,gauss);
442 element->ViscousHeating(&phi,xyz_list,gauss,vx_input,vy_input,vz_input);
443
444 scalar_def=phi/rho_ice*Jdet*gauss->weight;
445 if(dt!=0.) scalar_def=scalar_def*dt;
446
447 /*TODO: add -beta*laplace T_m(p)*/
448 for(int i=0;i<numnodes;i++) pe->values[i]+=scalar_def*basis[i];
449
450 /* Build transient now */
451 if(reCast<bool,IssmDouble>(dt)){
452 enthalpy_input->GetInputValue(&enthalpy, gauss);
453 scalar_transient=enthalpy*Jdet*gauss->weight;
454 for(int i=0;i<numnodes;i++) pe->values[i]+=scalar_transient*basis[i];
455 }
456
457 if(stabilization==2){
458 element->NodalFunctionsDerivatives(dbasis,xyz_list,gauss);
459
460 vx_input->GetInputValue(&u,gauss);
461 vy_input->GetInputValue(&v,gauss);
462 vz_input->GetInputValue(&w,gauss);
463 tau_parameter=element->StabilizationParameter(u,v,w,diameter,kappa/rho_ice);
464
465 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]);
466
467 if(dt!=0.){
468 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]);
469 }
470 }
471 }
472
473 /*Clean up and return*/
474 xDelete<IssmDouble>(basis);
475 xDelete<IssmDouble>(dbasis);
476 xDelete<IssmDouble>(xyz_list);
477 delete gauss;
478 return pe;
479
480}/*}}}*/
481ElementVector* EnthalpyAnalysis::CreatePVectorSheet(Element* element){/*{{{*/
482
483 /* Geothermal flux on ice sheet base and basal friction */
484 if(!element->IsOnBed() || element->IsFloating()) return NULL;
485
486 IssmDouble dt,Jdet,enthalpy,pressure,watercolumn,geothermalflux,vx,vy,vz;
487 IssmDouble enthalpyup,pressureup,alpha2,scalar,basalfriction,heatflux;
488 IssmDouble *xyz_list_base = NULL;
489
490 /*Fetch number of nodes for this finite element*/
491 int numnodes = element->GetNumberOfNodes();
492
493 /*Initialize vectors*/
494 ElementVector* pe = element->NewElementVector();
495 IssmDouble* basis = xNew<IssmDouble>(numnodes);
496
497 /*Retrieve all inputs and parameters*/
498 element->GetVerticesCoordinatesBase(&xyz_list_base);
499 element->FindParam(&dt,TimesteppingTimeStepEnum);
500 Input* vx_input = element->GetInput(VxEnum); _assert_(vx_input);
501 Input* vy_input = element->GetInput(VyEnum); _assert_(vy_input);
502 Input* vz_input = element->GetInput(VzEnum); _assert_(vz_input);
503 Input* enthalpy_input = element->GetInput(EnthalpyPicardEnum); _assert_(enthalpy_input);
504 Input* pressure_input = element->GetInput(PressureEnum); _assert_(pressure_input);
505 Input* geothermalflux_input = element->GetInput(BasalforcingsGeothermalfluxEnum); _assert_(geothermalflux_input);
506 Input* watercolumn_input = element->GetInput(WatercolumnEnum); _assert_(watercolumn_input);
507 IssmDouble rho_ice = element->GetMaterialParameter(MaterialsRhoIceEnum);
508 IssmDouble heatcapacity = element->GetMaterialParameter(MaterialsHeatcapacityEnum);
509
510 /*Build friction element, needed later: */
511 Friction* friction=new Friction(element,3);
512
513 /* Start looping on the number of gaussian points: */
514 Gauss* gauss = element->NewGaussBase(2);
515 Gauss* gaussup = element->NewGaussTop(2);
516 for(int ig=gauss->begin();ig<gauss->end();ig++){
517 gauss->GaussPoint(ig);
518
519 element->JacobianDeterminantBase(&Jdet,xyz_list_base,gauss);
520 element->NodalFunctions(basis,gauss);
521
522 enthalpy_input->GetInputValue(&enthalpy,gauss);
523 pressure_input->GetInputValue(&pressure,gauss);
524 watercolumn_input->GetInputValue(&watercolumn,gauss);
525
526 if((watercolumn<=0.) && (enthalpy < PureIceEnthalpy(element,pressure))){
527 /* the above check is equivalent to
528 NOT ((watercolumn>0.) AND (enthalpy<PIE)) AND (enthalpy<PIE)*/
529 geothermalflux_input->GetInputValue(&geothermalflux,gauss);
530
531 friction->GetAlpha2(&alpha2,gauss,vx_input,vy_input,vz_input);
532 vx_input->GetInputValue(&vx,gauss);
533 vy_input->GetInputValue(&vy,gauss);
534 vz_input->GetInputValue(&vz,gauss);
535 basalfriction = alpha2*(vx*vx + vy*vy + vz*vz);
536 heatflux = (basalfriction+geothermalflux)/(rho_ice);
537
538 scalar = gauss->weight*Jdet*heatflux;
539 if(dt!=0.) scalar=dt*scalar;
540
541 for(int i=0;i<numnodes;i++) pe->values[i]+=scalar*basis[i];
542 }
543 else if(enthalpy >= PureIceEnthalpy(element,pressure)){
544 /* check positive thickness of temperate basal ice layer */
545 enthalpy_input->GetInputValue(&enthalpyup,gaussup);
546 pressure_input->GetInputValue(&pressureup,gaussup);
547 if(enthalpyup >= PureIceEnthalpy(element,pressureup)){
548 // TODO: temperate ice has positive thickness: grad enthalpy*n=0.
549 }
550 else{
551 // only base temperate, set Dirichlet BCs in Penta::UpdateBasalConstraintsEnthalpy()
552 }
553 }
554 else{
555 // base cold, but watercolumn positive. Set base to h_pmp.
556 }
557 }
558
559 /*Clean up and return*/
560 delete gauss;
561 delete gaussup;
562 delete friction;
563 xDelete<IssmDouble>(basis);
564 xDelete<IssmDouble>(xyz_list_base);
565 return pe;
566
567}/*}}}*/
568ElementVector* EnthalpyAnalysis::CreatePVectorShelf(Element* element){/*{{{*/
569
570 /*Get basal element*/
571 if(!element->IsOnBed() || !element->IsFloating()) return NULL;
572
573 IssmDouble h_pmp,dt,Jdet,scalar_ocean,pressure;
574 IssmDouble *xyz_list_base = NULL;
575
576 /*Fetch number of nodes for this finite element*/
577 int numnodes = element->GetNumberOfNodes();
578
579 /*Initialize vectors*/
580 ElementVector* pe = element->NewElementVector();
581 IssmDouble* basis = xNew<IssmDouble>(numnodes);
582
583 /*Retrieve all inputs and parameters*/
584 element->GetVerticesCoordinatesBase(&xyz_list_base);
585 element->FindParam(&dt,TimesteppingTimeStepEnum);
586 Input* pressure_input=element->GetInput(PressureEnum); _assert_(pressure_input);
587 IssmDouble gravity = element->GetMaterialParameter(ConstantsGEnum);
588 IssmDouble rho_water = element->GetMaterialParameter(MaterialsRhoWaterEnum);
589 IssmDouble rho_ice = element->GetMaterialParameter(MaterialsRhoIceEnum);
590 IssmDouble heatcapacity = element->GetMaterialParameter(MaterialsHeatcapacityEnum);
591 IssmDouble mixed_layer_capacity= element->GetMaterialParameter(MaterialsMixedLayerCapacityEnum);
592 IssmDouble thermal_exchange_vel= element->GetMaterialParameter(MaterialsThermalExchangeVelocityEnum);
593
594 /* Start looping on the number of gaussian points: */
595 Gauss* gauss=element->NewGaussBase(2);
596 for(int ig=gauss->begin();ig<gauss->end();ig++){
597 gauss->GaussPoint(ig);
598
599 element->JacobianDeterminantBase(&Jdet,xyz_list_base,gauss);
600 element->NodalFunctions(basis,gauss);
601
602 pressure_input->GetInputValue(&pressure,gauss);
603 h_pmp=element->PureIceEnthalpy(pressure);
604
605 scalar_ocean=gauss->weight*Jdet*rho_water*mixed_layer_capacity*thermal_exchange_vel*h_pmp/(heatcapacity*rho_ice);
606 if(reCast<bool,IssmDouble>(dt)) scalar_ocean=dt*scalar_ocean;
607
608 for(int i=0;i<numnodes;i++) pe->values[i]+=scalar_ocean*basis[i];
609 }
610
611 /*Clean up and return*/
612 delete gauss;
613 xDelete<IssmDouble>(basis);
614 xDelete<IssmDouble>(xyz_list_base);
615 return pe;
616}/*}}}*/
617void EnthalpyAnalysis::GetBConduct(IssmDouble* B,Element* element,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/
618 /*Compute B matrix. B=[B1 B2 B3 B4 B5 B6] where Bi is of size 5*NDOF1.
619 * For node i, Bi' can be expressed in the actual coordinate system
620 * by:
621 * Bi_conduct=[ dh/dx ]
622 * [ dh/dy ]
623 * [ dh/dz ]
624 * where h is the interpolation function for node i.
625 *
626 * We assume B has been allocated already, of size: 3x(NDOF1*numnodes)
627 */
628
629 /*Fetch number of nodes for this finite element*/
630 int numnodes = element->GetNumberOfNodes();
631
632 /*Get nodal functions derivatives*/
633 IssmDouble* dbasis=xNew<IssmDouble>(3*numnodes);
634 element->NodalFunctionsDerivatives(dbasis,xyz_list,gauss);
635
636 /*Build B: */
637 for(int i=0;i<numnodes;i++){
638 B[numnodes*0+i] = dbasis[0*numnodes+i];
639 B[numnodes*1+i] = dbasis[1*numnodes+i];
640 B[numnodes*2+i] = dbasis[2*numnodes+i];
641 }
642
643 /*Clean-up*/
644 xDelete<IssmDouble>(dbasis);
645}/*}}}*/
646void EnthalpyAnalysis::GetBAdvec(IssmDouble* B,Element* element,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/
647 /*Compute B matrix. B=[B1 B2 B3 B4 B5 B6] where Bi is of size 5*NDOF1.
648 * For node i, Bi' can be expressed in the actual coordinate system
649 * by:
650 * Bi_advec =[ h ]
651 * [ h ]
652 * [ h ]
653 * where h is the interpolation function for node i.
654 *
655 * We assume B has been allocated already, of size: 3x(NDOF1*NUMNODESP1)
656 */
657
658 /*Fetch number of nodes for this finite element*/
659 int numnodes = element->GetNumberOfNodes();
660
661 /*Get nodal functions*/
662 IssmDouble* basis=xNew<IssmDouble>(numnodes);
663 element->NodalFunctions(basis,gauss);
664
665 /*Build B: */
666 for(int i=0;i<numnodes;i++){
667 B[numnodes*0+i] = basis[i];
668 B[numnodes*1+i] = basis[i];
669 B[numnodes*2+i] = basis[i];
670 }
671
672 /*Clean-up*/
673 xDelete<IssmDouble>(basis);
674}/*}}}*/
675void EnthalpyAnalysis::GetBAdvecprime(IssmDouble* B,Element* element,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/
676 /*Compute B matrix. B=[B1 B2 B3 B4 B5 B6] where Bi is of size 5*NDOF1.
677 * For node i, Bi' can be expressed in the actual coordinate system
678 * by:
679 * Biprime_advec=[ dh/dx ]
680 * [ dh/dy ]
681 * [ dh/dz ]
682 * where h is the interpolation function for node i.
683 *
684 * We assume B has been allocated already, of size: 3x(NDOF1*numnodes)
685 */
686
687 /*Fetch number of nodes for this finite element*/
688 int numnodes = element->GetNumberOfNodes();
689
690 /*Get nodal functions derivatives*/
691 IssmDouble* dbasis=xNew<IssmDouble>(3*numnodes);
692 element->NodalFunctionsDerivatives(dbasis,xyz_list,gauss);
693
694 /*Build B: */
695 for(int i=0;i<numnodes;i++){
696 B[numnodes*0+i] = dbasis[0*numnodes+i];
697 B[numnodes*1+i] = dbasis[1*numnodes+i];
698 B[numnodes*2+i] = dbasis[2*numnodes+i];
699 }
700
701 /*Clean-up*/
702 xDelete<IssmDouble>(dbasis);
703}/*}}}*/
704void EnthalpyAnalysis::GetSolutionFromInputs(Vector<IssmDouble>* solution,Element* element){/*{{{*/
705 element->GetSolutionFromInputsOneDof(solution,EnthalpyEnum);
706}/*}}}*/
707void EnthalpyAnalysis::InputUpdateFromSolution(IssmDouble* solution,Element* element){/*{{{*/
708
709 bool converged;
710 int i,rheology_law;
711 IssmDouble B_average,s_average,T_average=0.,P_average=0.;
712 int *doflist = NULL;
713 IssmDouble *xyz_list = NULL;
714
715 /*Fetch number of nodes and dof for this finite element*/
716 int numnodes = element->GetNumberOfNodes();
717
718 /*Fetch dof list and allocate solution vector*/
719 element->GetDofList(&doflist,NoneApproximationEnum,GsetEnum);
720 IssmDouble* values = xNew<IssmDouble>(numnodes);
721 IssmDouble* pressure = xNew<IssmDouble>(numnodes);
722 IssmDouble* surface = xNew<IssmDouble>(numnodes);
723 IssmDouble* B = xNew<IssmDouble>(numnodes);
724 IssmDouble* temperature = xNew<IssmDouble>(numnodes);
725 IssmDouble* waterfraction = xNew<IssmDouble>(numnodes);
726
727 /*Use the dof list to index into the solution vector: */
728 for(i=0;i<numnodes;i++){
729 values[i]=solution[doflist[i]];
730
731 /*Check solution*/
732 if(xIsNan<IssmDouble>(values[i])) _error_("NaN found in solution vector");
733 }
734
735 /*Get all inputs and parameters*/
736 element->GetInputValue(&converged,ConvergedEnum);
737 element->GetInputListOnNodes(&pressure[0],PressureEnum);
738 if(converged){
739 for(i=0;i<numnodes;i++){
740 element->EnthalpyToThermal(&temperature[i],&waterfraction[i],values[i],pressure[i]);
741 if(waterfraction[i]<0.) _error_("Negative water fraction found in solution vector");
742 if(waterfraction[i]>1.) _error_("Water fraction >1 found in solution vector");
743 }
744 element->AddInput(EnthalpyEnum,values,P1Enum);
745 element->AddInput(WaterfractionEnum,waterfraction,P1Enum);
746 element->AddInput(TemperatureEnum,temperature,P1Enum);
747
748 /*Update Rheology only if converged (we must make sure that the temperature is below melting point
749 * otherwise the rheology could be negative*/
750 element->FindParam(&rheology_law,MaterialsRheologyLawEnum);
751 element->GetInputListOnNodes(&surface[0],SurfaceEnum);
752 switch(rheology_law){
753 case NoneEnum:
754 /*Do nothing: B is not temperature dependent*/
755 break;
756 case PatersonEnum:
757 for(i=0;i<numnodes;i++) B[i]=Paterson(temperature[i]);
758 element->AddMaterialInput(MaterialsRheologyBEnum,&B[0],P1Enum);
759 break;
760 case ArrheniusEnum:
761 element->GetVerticesCoordinates(&xyz_list);
762 for(i=0;i<numnodes;i++) B[i]=Arrhenius(temperature[i],surface[i]-xyz_list[i*3+2],element->GetMaterialParameter(MaterialsRheologyNEnum));
763 element->AddMaterialInput(MaterialsRheologyBEnum,&B[0],P1Enum);
764 break;
765 case LliboutryDuvalEnum:
766 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));
767 element->AddMaterialInput(MaterialsRheologyBEnum,&B[0],P1Enum);
768 break;
769 default: _error_("Rheology law " << EnumToStringx(rheology_law) << " not supported yet");
770 }
771 }
772 else{
773 element->AddInput(EnthalpyPicardEnum,values,P1Enum);
774 }
775
776 /*Free ressources:*/
777 xDelete<IssmDouble>(values);
778 xDelete<IssmDouble>(pressure);
779 xDelete<IssmDouble>(surface);
780 xDelete<IssmDouble>(B);
781 xDelete<IssmDouble>(temperature);
782 xDelete<IssmDouble>(waterfraction);
783 xDelete<IssmDouble>(xyz_list);
784 xDelete<int>(doflist);
785}/*}}}*/
786
787/*Intermediaries*/
788IssmDouble EnthalpyAnalysis::EnthalpyDiffusionParameter(Element* element,IssmDouble enthalpy,IssmDouble pressure){/*{{{*/
789
790 IssmDouble heatcapacity = element->GetMaterialParameter(MaterialsHeatcapacityEnum);
791 IssmDouble temperateiceconductivity = element->GetMaterialParameter(MaterialsTemperateiceconductivityEnum);
792 IssmDouble thermalconductivity = element->GetMaterialParameter(MaterialsThermalconductivityEnum);
793
794 if(enthalpy < PureIceEnthalpy(element,pressure)){
795 return thermalconductivity/heatcapacity;
796 }
797 else{
798 return temperateiceconductivity/heatcapacity;
799 }
800}/*}}}*/
801IssmDouble EnthalpyAnalysis::EnthalpyDiffusionParameterVolume(Element* element,int enthalpy_enum){/*{{{*/
802
803 int iv;
804 IssmDouble lambda; /* fraction of cold ice */
805 IssmDouble kappa ,kappa_c,kappa_t; /* enthalpy conductivities */
806 IssmDouble Hc,Ht;
807
808
809 /*Get pressures and enthalpies on vertices*/
810 int numvertices = element->GetNumberOfVertices();
811 IssmDouble* pressures = xNew<IssmDouble>(numvertices);
812 IssmDouble* enthalpies = xNew<IssmDouble>(numvertices);
813 IssmDouble* PIE = xNew<IssmDouble>(numvertices);
814 IssmDouble* dHpmp = xNew<IssmDouble>(numvertices);
815 element->GetInputListOnVertices(pressures,PressureEnum);
816 element->GetInputListOnVertices(enthalpies,enthalpy_enum);
817 for(iv=0;iv<numvertices;iv++){
818 PIE[iv] = PureIceEnthalpy(element,pressures[iv]);
819 dHpmp[iv] = enthalpies[iv]-PIE[iv];
820 }
821
822 bool allequalsign = true;
823 if(dHpmp[0]<0.){
824 for(iv=1; iv<numvertices;iv++) allequalsign=(allequalsign && (dHpmp[iv]<0.));
825 }
826 else{
827 for(iv=1; iv<numvertices;iv++) allequalsign=(allequalsign && (dHpmp[iv]>=0.));
828 }
829
830 if(allequalsign){
831 kappa = EnthalpyDiffusionParameter(element,enthalpies[0],pressures[0]);
832 }
833 else{
834 /* return harmonic mean of thermal conductivities, weighted by fraction of cold/temperate ice,
835 cf Patankar 1980, pp44 */
836 kappa_c = EnthalpyDiffusionParameter(element,PureIceEnthalpy(element,0.)-1.,0.);
837 kappa_t = EnthalpyDiffusionParameter(element,PureIceEnthalpy(element,0.)+1.,0.);
838 Hc=0.; Ht=0.;
839 for(iv=0; iv<numvertices;iv++){
840 if(enthalpies[iv]<PIE[iv])
841 Hc+=(PIE[iv]-enthalpies[iv]);
842 else
843 Ht+=(enthalpies[iv]-PIE[iv]);
844 }
845 _assert_((Hc+Ht)>0.);
846 lambda = Hc/(Hc+Ht);
847 kappa = 1./(lambda/kappa_c + (1.-lambda)/kappa_t);
848 }
849
850 /*Clean up and return*/
851 xDelete<IssmDouble>(PIE);
852 xDelete<IssmDouble>(dHpmp);
853 xDelete<IssmDouble>(pressures);
854 xDelete<IssmDouble>(enthalpies);
855 return kappa;
856}
857/*}}}*/
858IssmDouble EnthalpyAnalysis::PureIceEnthalpy(Element* element,IssmDouble pressure){/*{{{*/
859
860 IssmDouble heatcapacity = element->GetMaterialParameter(MaterialsHeatcapacityEnum);
861 IssmDouble referencetemperature = element->GetMaterialParameter(ConstantsReferencetemperatureEnum);
862
863 return heatcapacity*(TMeltingPoint(element,pressure)-referencetemperature);
864}/*}}}*/
865IssmDouble EnthalpyAnalysis::TMeltingPoint(Element* element,IssmDouble pressure){/*{{{*/
866
867 IssmDouble meltingpoint = element->GetMaterialParameter(MaterialsMeltingpointEnum);
868 IssmDouble beta = element->GetMaterialParameter(MaterialsBetaEnum);
869
870 return meltingpoint-beta*pressure;
871}/*}}}*/
Note: See TracBrowser for help on using the repository browser.