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

Last change on this file since 27470 was 27470, checked in by Mathieu Morlighem, 3 years ago

NEW: moving friction io to Friction.cpp and preparing linearization as an option
File size: 62.4 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#include "../solutionsequences/solutionsequences.h"
7#include "../cores/cores.h"
8
9/*Model processing*/
10void EnthalpyAnalysis::CreateConstraints(Constraints* constraints,IoModel* iomodel){/*{{{*/
11
12 /*Intermediary*/
13 int count;
14 int M,N;
15 bool spcpresent = false;
16 int finiteelement;
17 IssmDouble heatcapacity;
18 IssmDouble referencetemperature;
19
20 /*Output*/
21 IssmDouble *spcvector = NULL;
22 IssmDouble *spcvectorstatic = NULL;
23 IssmDouble* times=NULL;
24 IssmDouble* values=NULL;
25 IssmDouble* issurface = NULL;
26
27 /*Fetch parameters: */
28 iomodel->FindConstant(&heatcapacity,"md.materials.heatcapacity");
29 iomodel->FindConstant(&referencetemperature,"md.constants.referencetemperature");
30 iomodel->FindConstant(&finiteelement,"md.thermal.fe");
31
32 /*return if 2d mesh*/
33 if(iomodel->domaintype==Domain2DhorizontalEnum) return;
34
35 /*Fetch data: */
36 iomodel->FetchData(&issurface,&M,&N,"md.mesh.vertexonsurface"); _assert_(N>0); _assert_(M==iomodel->numberofvertices);
37 iomodel->FetchData(&spcvector,&M,&N,"md.thermal.spctemperature");
38 iomodel->FetchData(&spcvectorstatic,&M,&N,"md.thermal.spctemperature");
39
40 /*Specific case for PDD, we want the constaints to be updated by the PDD scheme itself*/
41 bool isdynamic = false;
42 if (iomodel->solution_enum==TransientSolutionEnum){
43 int smb_model;
44 iomodel->FindConstant(&smb_model,"md.smb.model");
45 if(smb_model==SMBpddEnum) isdynamic=true;
46 if(smb_model==SMBd18opddEnum) isdynamic=true;
47 if(smb_model==SMBpddSicopolisEnum) isdynamic=true;
48 }
49
50 /*Convert spcs from temperatures to enthalpy*/
51 _assert_(N>0); _assert_(M>=iomodel->numberofvertices);
52 for(int i=0;i<iomodel->numberofvertices;i++){
53 for(int j=0;j<N;j++){
54 if (isdynamic){
55 if (issurface[i]==1){
56 spcvector[i*N+j] = heatcapacity*(spcvector[i*N+j]-referencetemperature);
57 spcvectorstatic[i*N+j] = NAN;
58 }
59 else{
60 spcvector[i*N+j] = NAN;
61 spcvectorstatic[i*N+j] = heatcapacity*(spcvectorstatic[i*N+j]-referencetemperature);
62 }
63 }
64 else{
65 spcvector[i*N+j] = heatcapacity*(spcvector[i*N+j]-referencetemperature);
66 }
67 }
68 }
69
70 if(isdynamic){
71 IoModelToDynamicConstraintsx(constraints,iomodel,spcvector,iomodel->numberofvertices,1,EnthalpyAnalysisEnum,finiteelement);
72 IoModelToConstraintsx(constraints,iomodel,spcvectorstatic,M,N,EnthalpyAnalysisEnum,finiteelement);
73 }
74 else{
75 IoModelToConstraintsx(constraints,iomodel,spcvector,M,N,EnthalpyAnalysisEnum,finiteelement);
76 }
77
78 /*Free resources:*/
79 iomodel->DeleteData(spcvector,"md.thermal.spctemperature");
80 iomodel->DeleteData(spcvectorstatic,"md.thermal.spctemperature");
81 iomodel->DeleteData(issurface,"md.mesh.vertexonsurface");
82 xDelete<IssmDouble>(times);
83 xDelete<IssmDouble>(values);
84}/*}}}*/
85void EnthalpyAnalysis::CreateLoads(Loads* loads, IoModel* iomodel){/*{{{*/
86
87 /*No loads */
88}/*}}}*/
89void EnthalpyAnalysis::CreateNodes(Nodes* nodes,IoModel* iomodel,bool isamr){/*{{{*/
90
91 int finiteelement;
92 iomodel->FindConstant(&finiteelement,"md.thermal.fe");
93
94 if(iomodel->domaintype==Domain3DEnum) iomodel->FetchData(2,"md.mesh.vertexonbase","md.mesh.vertexonsurface");
95 ::CreateNodes(nodes,iomodel,EnthalpyAnalysisEnum,finiteelement);
96 iomodel->DeleteData(2,"md.mesh.vertexonbase","md.mesh.vertexonsurface");
97}/*}}}*/
98int EnthalpyAnalysis::DofsPerNode(int** doflist,int domaintype,int approximation){/*{{{*/
99 return 1;
100}/*}}}*/
101void EnthalpyAnalysis::UpdateElements(Elements* elements,Inputs* inputs,IoModel* iomodel,int analysis_counter,int analysis_type){/*{{{*/
102
103 bool dakota_analysis,ismovingfront,isenthalpy;
104 int basalforcing_model,materialstype;
105
106 /*Now, is the model 3d? otherwise, do nothing: */
107 if(iomodel->domaintype==Domain2DhorizontalEnum)return;
108
109 /*Is enthalpy requested?*/
110 iomodel->FindConstant(&isenthalpy,"md.thermal.isenthalpy");
111 if(!isenthalpy) return;
112
113 /*Fetch data needed: */
114 iomodel->FetchData(3,"md.initialization.temperature","md.initialization.waterfraction","md.initialization.pressure");
115
116 /*Finite element type*/
117 int finiteelement;
118 iomodel->FindConstant(&finiteelement,"md.thermal.fe");
119
120 /*Update elements: */
121 int counter=0;
122 for(int i=0;i<iomodel->numberofelements;i++){
123 if(iomodel->my_elements[i]){
124 Element* element=(Element*)elements->GetObjectByOffset(counter);
125 element->Update(inputs,i,iomodel,analysis_counter,analysis_type,finiteelement);
126 counter++;
127 }
128 }
129
130 iomodel->FindConstant(&dakota_analysis,"md.qmu.isdakota");
131 iomodel->FindConstant(&ismovingfront,"md.transient.ismovingfront");
132 iomodel->FindConstant(&materialstype,"md.materials.type");
133
134 iomodel->FetchDataToInput(inputs,elements,"md.geometry.thickness",ThicknessEnum);
135 iomodel->FetchDataToInput(inputs,elements,"md.geometry.surface",SurfaceEnum);
136 iomodel->FetchDataToInput(inputs,elements,"md.initialization.sealevel",SealevelEnum,0);
137 iomodel->FetchDataToInput(inputs,elements,"md.geometry.base",BaseEnum);
138 iomodel->FetchDataToInput(inputs,elements,"md.mask.ice_levelset",MaskIceLevelsetEnum);
139 iomodel->FetchDataToInput(inputs,elements,"md.mask.ocean_levelset",MaskOceanLevelsetEnum);
140 if(iomodel->domaintype!=Domain2DhorizontalEnum){
141 iomodel->FetchDataToInput(inputs,elements,"md.mesh.vertexonbase",MeshVertexonbaseEnum);
142 iomodel->FetchDataToInput(inputs,elements,"md.mesh.vertexonsurface",MeshVertexonsurfaceEnum);
143 }
144 iomodel->FetchDataToInput(inputs,elements,"md.initialization.pressure",PressureEnum);
145 iomodel->FetchDataToInput(inputs,elements,"md.initialization.temperature",TemperatureEnum);
146 iomodel->FetchDataToInput(inputs,elements,"md.initialization.waterfraction",WaterfractionEnum);
147 iomodel->FetchDataToInput(inputs,elements,"md.initialization.enthalpy",EnthalpyEnum);
148 iomodel->FetchDataToInput(inputs,elements,"md.initialization.watercolumn",WatercolumnEnum);
149 iomodel->FetchDataToInput(inputs,elements,"md.basalforcings.groundedice_melting_rate",BasalforcingsGroundediceMeltingRateEnum);
150 iomodel->FetchDataToInput(inputs,elements,"md.initialization.vx",VxEnum);
151 iomodel->FetchDataToInput(inputs,elements,"md.initialization.vy",VyEnum);
152 iomodel->FetchDataToInput(inputs,elements,"md.initialization.vz",VzEnum);
153 InputUpdateFromConstantx(inputs,elements,0.,VxMeshEnum);
154 InputUpdateFromConstantx(inputs,elements,0.,VyMeshEnum);
155 InputUpdateFromConstantx(inputs,elements,0.,VzMeshEnum);
156 if(ismovingfront){
157 iomodel->FetchDataToInput(inputs,elements,"md.mesh.vertexonbase",MeshVertexonbaseEnum); // required for updating active nodes
158 }
159
160 /*Basal forcings variables*/
161 iomodel->FindConstant(&basalforcing_model,"md.basalforcings.model");
162 switch(basalforcing_model){
163 case MantlePlumeGeothermalFluxEnum:
164 break;
165 default:
166 iomodel->FetchDataToInput(inputs,elements,"md.basalforcings.geothermalflux",BasalforcingsGeothermalfluxEnum);
167 break;
168 }
169
170 /*Rheology type*/
171 iomodel->FetchDataToInput(inputs,elements,"md.materials.rheology_B",MaterialsRheologyBEnum);
172 switch(materialstype){
173 case MatenhancediceEnum:
174 iomodel->FetchDataToInput(inputs,elements,"md.materials.rheology_n",MaterialsRheologyNEnum);
175 iomodel->FetchDataToInput(inputs,elements,"md.materials.rheology_E",MaterialsRheologyEEnum);
176 break;
177 case MatdamageiceEnum:
178 iomodel->FetchDataToInput(inputs,elements,"md.materials.rheology_n",MaterialsRheologyNEnum);
179 break;
180 case MatestarEnum:
181 iomodel->FetchDataToInput(inputs,elements,"md.materials.rheology_Ec",MaterialsRheologyEcEnum);
182 iomodel->FetchDataToInput(inputs,elements,"md.materials.rheology_Es",MaterialsRheologyEsEnum);
183 break;
184 case MaticeEnum:
185 iomodel->FetchDataToInput(inputs,elements,"md.materials.rheology_n",MaterialsRheologyNEnum);
186 break;
187 default:
188 _error_("not supported");
189 }
190
191 /*Friction*/
192 FrictionUpdateInputs(elements, inputs, iomodel);
193
194 /*Free data: */
195 iomodel->DeleteData(3,"md.initialization.temperature","md.initialization.waterfraction","md.initialization.pressure");
196
197}/*}}}*/
198void EnthalpyAnalysis::UpdateParameters(Parameters* parameters,IoModel* iomodel,int solution_enum,int analysis_enum){/*{{{*/
199
200 int numoutputs;
201 char** requestedoutputs = NULL;
202
203 parameters->AddObject(iomodel->CopyConstantObject("md.thermal.stabilization",ThermalStabilizationEnum));
204 parameters->AddObject(iomodel->CopyConstantObject("md.thermal.maxiter",ThermalMaxiterEnum));
205 parameters->AddObject(iomodel->CopyConstantObject("md.thermal.reltol",ThermalReltolEnum));
206 parameters->AddObject(iomodel->CopyConstantObject("md.thermal.isenthalpy",ThermalIsenthalpyEnum));
207 parameters->AddObject(iomodel->CopyConstantObject("md.thermal.isdynamicbasalspc",ThermalIsdynamicbasalspcEnum));
208 parameters->AddObject(iomodel->CopyConstantObject("md.thermal.isdrainicecolumn",ThermalIsdrainicecolumnEnum));
209 parameters->AddObject(iomodel->CopyConstantObject("md.thermal.watercolumn_upperlimit",ThermalWatercolumnUpperlimitEnum));
210
211 iomodel->FindConstant(&requestedoutputs,&numoutputs,"md.thermal.requested_outputs");
212 parameters->AddObject(new IntParam(ThermalNumRequestedOutputsEnum,numoutputs));
213 if(numoutputs)parameters->AddObject(new StringArrayParam(ThermalRequestedOutputsEnum,requestedoutputs,numoutputs));
214 iomodel->DeleteData(&requestedoutputs,numoutputs,"md.thermal.requested_outputs");
215
216 /*Friction*/
217 FrictionUpdateParameters(parameters, iomodel);
218}/*}}}*/
219
220/*Finite Element Analysis*/
221void EnthalpyAnalysis::ApplyBasalConstraints(IssmDouble* local_spc,Element* element){/*{{{*/
222
223 /* Do not check if ice in element, this may lead to inconsistencies between cpu partitions */
224 /* Only update constraints at the base. */
225 if(!(element->IsOnBase())) return;
226
227 /*Intermediary*/
228 bool isdynamicbasalspc;
229 int numindices;
230 int *indices = NULL;
231 IssmDouble pressure;
232
233 /*Check wether dynamic basal boundary conditions are activated */
234 element->FindParam(&isdynamicbasalspc,ThermalIsdynamicbasalspcEnum);
235 if(!isdynamicbasalspc) return;
236
237 /*Get parameters and inputs: */
238 Input* pressure_input = element->GetInput(PressureEnum); _assert_(pressure_input);
239
240 /*Fetch indices of basal & surface nodes for this finite element*/
241 Penta *penta = (Penta *) element; // TODO: add Basal-/SurfaceNodeIndices to element.h, and change this to Element*
242 penta->BasalNodeIndices(&numindices,&indices,element->GetElementType());
243
244 GaussPenta* gauss=new GaussPenta();
245 for(int i=0;i<numindices;i++){
246 gauss->GaussNode(element->GetElementType(),indices[i]);
247
248 pressure_input->GetInputValue(&pressure,gauss);
249
250 /*apply or release spc*/
251 Node* node=element->GetNode(indices[i]);
252 if(!node->IsActive()) continue;
253 if(local_spc[node->Lid()]==1.){
254 pressure_input->GetInputValue(&pressure, gauss);
255 node->ApplyConstraint(0,PureIceEnthalpy(element,pressure));
256 }
257 else {
258 node->DofInFSet(0);
259 }
260 }
261
262 /*Free resources:*/
263 xDelete<int>(indices);
264 delete gauss;
265}/*}}}*/
266void EnthalpyAnalysis::ComputeBasalMeltingrate(FemModel* femmodel){/*{{{*/
267 /*Compute basal melting rates: */
268 for(Object* & object : femmodel->elements->objects){
269 Element* element=xDynamicCast<Element*>(object);
270 ComputeBasalMeltingrate(element);
271 }
272
273 /*extrude inputs*/
274 femmodel->parameters->SetParam(BasalforcingsGroundediceMeltingRateEnum,InputToExtrudeEnum);
275 extrudefrombase_core(femmodel);
276}/*}}}*/
277void EnthalpyAnalysis::ComputeBasalMeltingrate(Element* element){/*{{{*/
278 /*Calculate the basal melt rates of the enthalpy model after Aschwanden 2012*/
279 /* melting rate is positive when melting, negative when refreezing*/
280
281 /* Check if ice in element */
282 if(!element->IsIceInElement()) return;
283
284 /* Only compute melt rates at the base of grounded ice*/
285 if(!element->IsOnBase() || element->IsAllFloating()) return;
286
287 /* Intermediaries */
288 bool converged;
289 const int dim=3;
290 int i,is,state;
291 int nodedown,nodeup,numnodes,numsegments;
292 int enthalpy_enum;
293 IssmDouble vec_heatflux[dim],normal_base[dim],d1enthalpy[dim],d1pressure[dim];
294 IssmDouble basalfriction,alpha2,geothermalflux,heatflux;
295 IssmDouble dt,yts;
296 IssmDouble melting_overshoot,lambda;
297 IssmDouble vx,vy,vz;
298 IssmDouble *xyz_list = NULL;
299 IssmDouble *xyz_list_base = NULL;
300 int *pairindices = NULL;
301
302 /*Fetch parameters*/
303 element->GetVerticesCoordinates(&xyz_list);
304 element->GetVerticesCoordinatesBase(&xyz_list_base);
305 element->GetInputValue(&converged,ConvergedEnum);
306 element->FindParam(&dt,TimesteppingTimeStepEnum);
307 element->FindParam(&yts, ConstantsYtsEnum);
308
309 if(dt==0. && !converged) enthalpy_enum=EnthalpyPicardEnum;
310 else enthalpy_enum=EnthalpyEnum;
311
312 IssmDouble latentheat = element->FindParam(MaterialsLatentheatEnum);
313 IssmDouble rho_ice = element->FindParam(MaterialsRhoIceEnum);
314 IssmDouble rho_water = element->FindParam(MaterialsRhoFreshwaterEnum);
315 IssmDouble beta = element->FindParam(MaterialsBetaEnum);
316 IssmDouble kappa = EnthalpyDiffusionParameterVolume(element,enthalpy_enum); _assert_(kappa>=0.);
317 IssmDouble kappa_mix;
318
319 /*retrieve inputs*/
320 Input* enthalpy_input = element->GetInput(enthalpy_enum); _assert_(enthalpy_input);
321 Input* pressure_input = element->GetInput(PressureEnum); _assert_(pressure_input);
322 Input* geothermalflux_input = element->GetInput(BasalforcingsGeothermalfluxEnum); _assert_(geothermalflux_input);
323
324 /*Build friction element, needed later: */
325 Friction* friction=new Friction(element,3);
326
327 /******** MELTING RATES ************************************//*{{{*/
328 element->NormalBase(&normal_base[0],xyz_list_base);
329 element->VerticalSegmentIndicesBase(&pairindices,&numsegments);
330 IssmDouble* meltingrate_enthalpy = xNew<IssmDouble>(numsegments);
331 IssmDouble* heating = xNew<IssmDouble>(numsegments);
332
333 numnodes=element->GetNumberOfNodes();
334 IssmDouble* enthalpies = xNew<IssmDouble>(numnodes);
335 IssmDouble* pressures = xNew<IssmDouble>(numnodes);
336 IssmDouble* watercolumns = xNew<IssmDouble>(numnodes);
337 IssmDouble* basalmeltingrates = xNew<IssmDouble>(numnodes);
338 element->GetInputListOnNodes(enthalpies,enthalpy_enum);
339 element->GetInputListOnNodes(pressures,PressureEnum);
340 element->GetInputListOnNodes(watercolumns,WatercolumnEnum);
341 element->GetInputListOnNodes(basalmeltingrates,BasalforcingsGroundediceMeltingRateEnum);
342
343 IssmDouble watercolumnupperlimit = element->FindParam(ThermalWatercolumnUpperlimitEnum);
344
345 Gauss* gauss=element->NewGauss();
346 for(is=0;is<numsegments;is++){
347 nodedown = pairindices[is*2+0];
348 nodeup = pairindices[is*2+1];
349 gauss->GaussNode(element->GetElementType(),nodedown);
350
351 state=GetThermalBasalCondition(element, enthalpies[nodedown], enthalpies[nodeup], pressures[nodedown], pressures[nodeup], watercolumns[nodedown], basalmeltingrates[nodedown]);
352 switch (state) {
353 case 0:
354 // cold, dry base: apply basal surface forcing
355 for(i=0;i<3;i++) vec_heatflux[i]=0.;
356 break;
357 case 1: case 2: case 3:
358 // case 1 : cold, wet base: keep at pressure melting point
359 // case 2: temperate, thin refreezing base: release spc
360 // case 3: temperate, thin melting base: set spc
361 enthalpy_input->GetInputDerivativeValue(&d1enthalpy[0],xyz_list,gauss);
362 for(i=0;i<3;i++) vec_heatflux[i]=-kappa*d1enthalpy[i];
363 break;
364 case 4:
365 // temperate, thick melting base: set grad H*n=0
366 kappa_mix=GetWetIceConductivity(element, enthalpies[nodedown], pressures[nodedown]);
367 pressure_input->GetInputDerivativeValue(&d1pressure[0],xyz_list,gauss);
368 for(i=0;i<3;i++) vec_heatflux[i]=kappa_mix*beta*d1pressure[i];
369 break;
370 default:
371 _printf0_(" unknown thermal basal state found!");
372 }
373 if(state==0) meltingrate_enthalpy[is]=0.;
374 else{
375 /*heat flux along normal*/
376 heatflux=0.;
377 for(i=0;i<3;i++) heatflux+=(vec_heatflux[i])*normal_base[i];
378
379 /*basal friction*/
380 friction->GetAlpha2(&alpha2,gauss);
381 friction->GetBasalSlidingSpeeds(&vx, &vy, &vz, gauss);
382 basalfriction=alpha2*(vx*vx + vy*vy + vz*vz);
383 geothermalflux_input->GetInputValue(&geothermalflux,gauss);
384 /* -Mb= Fb-(q-q_geo)/((1-w)*L*rho), and (1-w)*rho=rho_ice, cf Aschwanden 2012, eqs.1, 2, 66*/
385 heating[is]=(heatflux+basalfriction+geothermalflux);
386 meltingrate_enthalpy[is]=heating[is]/(latentheat*rho_ice); // m/s water equivalent
387 }
388 }/*}}}*/
389
390 /******** UPDATE MELTINGRATES AND WATERCOLUMN **************//*{{{*/
391 for(is=0;is<numsegments;is++){
392 nodedown = pairindices[is*2+0];
393 nodeup = pairindices[is*2+1];
394 if(dt!=0.){
395 if(watercolumns[nodedown]+meltingrate_enthalpy[is]*dt<0.){ // prevent too much freeze on
396 lambda = -watercolumns[nodedown]/(dt*meltingrate_enthalpy[is]); _assert_(lambda>=0.); _assert_(lambda<1.);
397 watercolumns[nodedown]=0.;
398 basalmeltingrates[nodedown]=lambda*meltingrate_enthalpy[is]; // restrict freeze on only to size of watercolumn
399 enthalpies[nodedown]+=(1.-lambda)*dt/yts*meltingrate_enthalpy[is]*latentheat*rho_ice; // use rest of energy to cool down base: dE=L*m, m=(1-lambda)*meltingrate*rho_ice
400 }
401 else{
402 basalmeltingrates[nodedown]=meltingrate_enthalpy[is];
403 watercolumns[nodedown]+=dt*meltingrate_enthalpy[is];
404 }
405 if(watercolumns[nodedown]>watercolumnupperlimit) watercolumns[nodedown]=watercolumnupperlimit;
406 }
407 else{
408 basalmeltingrates[nodedown]=meltingrate_enthalpy[is];
409 if(watercolumns[nodedown]+meltingrate_enthalpy[is]<0.)
410 watercolumns[nodedown]=0.;
411 else
412 watercolumns[nodedown]+=meltingrate_enthalpy[is];
413 }
414 basalmeltingrates[nodedown]*=rho_water/rho_ice; // convert meltingrate from water to ice equivalent
415 _assert_(watercolumns[nodedown]>=0.);
416 }/*}}}*/
417
418 /*feed updated variables back into model*/
419 int finite_element = element->GetElementType(); if(finite_element==P1Enum) finite_element = P1DGEnum;
420 if(dt!=0.){
421 element->AddInput(enthalpy_enum,enthalpies,finite_element);
422 element->AddInput(WatercolumnEnum,watercolumns,finite_element);
423 }
424 element->AddInput(BasalforcingsGroundediceMeltingRateEnum,basalmeltingrates,P1DGEnum);
425
426 /*Clean up and return*/
427 delete gauss;
428 delete friction;
429 xDelete<int>(pairindices);
430 xDelete<IssmDouble>(enthalpies);
431 xDelete<IssmDouble>(pressures);
432 xDelete<IssmDouble>(watercolumns);
433 xDelete<IssmDouble>(basalmeltingrates);
434 xDelete<IssmDouble>(meltingrate_enthalpy);
435 xDelete<IssmDouble>(heating);
436 xDelete<IssmDouble>(xyz_list);
437 xDelete<IssmDouble>(xyz_list_base);
438}/*}}}*/
439void EnthalpyAnalysis::Core(FemModel* femmodel){/*{{{*/
440
441 IssmDouble dt;
442 bool isdynamicbasalspc;
443
444 femmodel->parameters->FindParam(&dt,TimesteppingTimeStepEnum);
445 femmodel->parameters->FindParam(&isdynamicbasalspc,ThermalIsdynamicbasalspcEnum);
446
447 if(VerboseSolution()) _printf0_(" computing enthalpy\n");
448 femmodel->SetCurrentConfiguration(EnthalpyAnalysisEnum);
449 if((dt>0.) && isdynamicbasalspc) UpdateBasalConstraints(femmodel);
450 solutionsequence_thermal_nonlinear(femmodel);
451
452 /*transfer enthalpy to enthalpy picard for the next step: */
453 InputDuplicatex(femmodel,EnthalpyEnum,EnthalpyPicardEnum);
454
455 PostProcessing(femmodel);
456
457}/*}}}*/
458void EnthalpyAnalysis::PreCore(FemModel* femmodel){/*{{{*/
459 _error_("not implemented");
460}/*}}}*/
461ElementVector* EnthalpyAnalysis::CreateDVector(Element* element){/*{{{*/
462 /*Default, return NULL*/
463 return NULL;
464}/*}}}*/
465ElementMatrix* EnthalpyAnalysis::CreateJacobianMatrix(Element* element){/*{{{*/
466_error_("Not implemented");
467}/*}}}*/
468ElementMatrix* EnthalpyAnalysis::CreateKMatrix(Element* element){/*{{{*/
469
470 /* Check if ice in element */
471 if(!element->IsIceInElement()) return NULL;
472
473 /*compute all stiffness matrices for this element*/
474 ElementMatrix* Ke1=CreateKMatrixVolume(element);
475 ElementMatrix* Ke2=CreateKMatrixShelf(element);
476 ElementMatrix* Ke =new ElementMatrix(Ke1,Ke2);
477
478 /*clean-up and return*/
479 delete Ke1;
480 delete Ke2;
481 return Ke;
482}/*}}}*/
483ElementMatrix* EnthalpyAnalysis::CreateKMatrixVolume(Element* element){/*{{{*/
484
485 /* Check if ice in element */
486 if(!element->IsIceInElement()) return NULL;
487
488 /*Intermediaries */
489 int stabilization;
490 IssmDouble Jdet,dt,u,v,w,um,vm,wm,vel;
491 IssmDouble h,hx,hy,hz,vx,vy,vz;
492 IssmDouble tau_parameter,diameter;
493 IssmDouble tau_parameter_anisotropic[2],tau_parameter_hor,tau_parameter_ver;
494 IssmDouble D_scalar;
495 IssmDouble* xyz_list = NULL;
496
497 /*Fetch number of nodes and dof for this finite element*/
498 int numnodes = element->GetNumberOfNodes();
499
500 /*Initialize Element vector and other vectors*/
501 ElementMatrix* Ke = element->NewElementMatrix();
502 IssmDouble* basis = xNew<IssmDouble>(numnodes);
503 IssmDouble* dbasis = xNew<IssmDouble>(3*numnodes);
504 IssmDouble K[3][3];
505
506 /*Retrieve all inputs and parameters*/
507 element->GetVerticesCoordinates(&xyz_list);
508 element->FindParam(&dt,TimesteppingTimeStepEnum);
509 element->FindParam(&stabilization,ThermalStabilizationEnum);
510 IssmDouble rho_water = element->FindParam(MaterialsRhoSeawaterEnum);
511 IssmDouble rho_ice = element->FindParam(MaterialsRhoIceEnum);
512 IssmDouble gravity = element->FindParam(ConstantsGEnum);
513 IssmDouble heatcapacity = element->FindParam(MaterialsHeatcapacityEnum);
514 IssmDouble thermalconductivity = element->FindParam(MaterialsThermalconductivityEnum);
515 Input* vx_input = element->GetInput(VxEnum); _assert_(vx_input);
516 Input* vy_input = element->GetInput(VyEnum); _assert_(vy_input);
517 Input* vz_input = element->GetInput(VzEnum); _assert_(vz_input);
518 Input* vxm_input = element->GetInput(VxMeshEnum); _assert_(vxm_input);
519 Input* vym_input = element->GetInput(VyMeshEnum); _assert_(vym_input);
520 Input* vzm_input = element->GetInput(VzMeshEnum); _assert_(vzm_input);
521
522 /*Enthalpy diffusion parameter*/
523 IssmDouble kappa=this->EnthalpyDiffusionParameterVolume(element,EnthalpyPicardEnum); _assert_(kappa>=0.);
524
525 /* Start looping on the number of gaussian points: */
526 Gauss* gauss=element->NewGauss(4);
527 while(gauss->next()){
528
529 element->JacobianDeterminant(&Jdet,xyz_list,gauss);
530 element->NodalFunctions(basis,gauss);
531 element->NodalFunctionsDerivatives(dbasis,xyz_list,gauss);
532
533 D_scalar=gauss->weight*Jdet;
534 if(dt!=0.) D_scalar=D_scalar*dt;
535
536 /*Conduction: */
537 for(int i=0;i<numnodes;i++){
538 for(int j=0;j<numnodes;j++){
539 Ke->values[i*numnodes+j] += D_scalar*kappa/rho_ice*(
540 dbasis[0*numnodes+j]*dbasis[0*numnodes+i] + dbasis[1*numnodes+j]*dbasis[1*numnodes+i] + dbasis[2*numnodes+j]*dbasis[2*numnodes+i]
541 );
542 }
543 }
544
545 /*Advection: */
546 vx_input->GetInputValue(&u,gauss); vxm_input->GetInputValue(&um,gauss); vx=u-um;
547 vy_input->GetInputValue(&v,gauss); vym_input->GetInputValue(&vm,gauss); vy=v-vm;
548 vz_input->GetInputValue(&w,gauss); vzm_input->GetInputValue(&wm,gauss); vz=w-wm;
549 for(int i=0;i<numnodes;i++){
550 for(int j=0;j<numnodes;j++){
551 Ke->values[i*numnodes+j] += D_scalar*(
552 vx*dbasis[0*numnodes+j]*basis[i] + vy*dbasis[1*numnodes+j]*basis[i] +vz*dbasis[2*numnodes+j]*basis[i]
553 );
554 }
555 }
556
557 /*Transient: */
558 if(dt!=0.){
559 D_scalar=gauss->weight*Jdet;
560 for(int i=0;i<numnodes;i++){
561 for(int j=0;j<numnodes;j++){
562 Ke->values[i*numnodes+j] += D_scalar*basis[j]*basis[i];
563 }
564 }
565 D_scalar=D_scalar*dt;
566 }
567
568 /*Artificial diffusivity*/
569 if(stabilization==1){
570 element->ElementSizes(&hx,&hy,&hz);
571 vel=sqrt(vx*vx + vy*vy + vz*vz)+1.e-14;
572 h=sqrt( pow(hx*vx/vel,2) + pow(hy*vy/vel,2) + pow(hz*vz/vel,2));
573 K[0][0]=h/(2.*vel)*fabs(vx*vx); K[0][1]=h/(2.*vel)*fabs(vx*vy); K[0][2]=h/(2.*vel)*fabs(vx*vz);
574 K[1][0]=h/(2.*vel)*fabs(vy*vx); K[1][1]=h/(2.*vel)*fabs(vy*vy); K[1][2]=h/(2.*vel)*fabs(vy*vz);
575 K[2][0]=h/(2.*vel)*fabs(vz*vx); K[2][1]=h/(2.*vel)*fabs(vz*vy); K[2][2]=h/(2.*vel)*fabs(vz*vz);
576 for(int i=0;i<3;i++) for(int j=0;j<3;j++) K[i][j] = D_scalar*K[i][j];
577
578 for(int i=0;i<numnodes;i++){
579 for(int j=0;j<numnodes;j++){
580 Ke->values[i*numnodes+j] += (
581 dbasis[0*numnodes+i] *(K[0][0]*dbasis[0*numnodes+j] + K[0][1]*dbasis[1*numnodes+j]+ K[0][2]*dbasis[2*numnodes+j]) +
582 dbasis[1*numnodes+i] *(K[1][0]*dbasis[0*numnodes+j] + K[1][1]*dbasis[1*numnodes+j]+ K[1][2]*dbasis[2*numnodes+j]) +
583 dbasis[2*numnodes+i] *(K[2][0]*dbasis[0*numnodes+j] + K[2][1]*dbasis[1*numnodes+j]+ K[2][2]*dbasis[2*numnodes+j])
584 );
585 }
586 }
587 }
588 /*SUPG*/
589 else if(stabilization==2){
590 element->NodalFunctionsDerivatives(dbasis,xyz_list,gauss);
591 diameter=element->MinEdgeLength(xyz_list);
592 tau_parameter=element->StabilizationParameter(u-um,v-vm,w-wm,diameter,kappa/rho_ice);
593 for(int i=0;i<numnodes;i++){
594 for(int j=0;j<numnodes;j++){
595 Ke->values[i*numnodes+j]+=tau_parameter*D_scalar*
596 ((u-um)*dbasis[0*numnodes+i]+(v-vm)*dbasis[1*numnodes+i]+(w-wm)*dbasis[2*numnodes+i])*
597 ((u-um)*dbasis[0*numnodes+j]+(v-vm)*dbasis[1*numnodes+j]+(w-wm)*dbasis[2*numnodes+j]);
598 }
599 }
600 if(dt!=0.){
601 D_scalar=gauss->weight*Jdet;
602 for(int i=0;i<numnodes;i++){
603 for(int j=0;j<numnodes;j++){
604 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]);
605 }
606 }
607 }
608 }
609 /*anisotropic SUPG*/
610 else if(stabilization==3){
611 element->NodalFunctionsDerivatives(dbasis,xyz_list,gauss);
612 element->ElementSizes(&hx,&hy,&hz);
613 element->StabilizationParameterAnisotropic(&tau_parameter_anisotropic[0],u-um,v-vm,w-wm,hx,hy,hz,kappa/rho_ice);
614 tau_parameter_hor=tau_parameter_anisotropic[0];
615 tau_parameter_ver=tau_parameter_anisotropic[1];
616 for(int i=0;i<numnodes;i++){
617 for(int j=0;j<numnodes;j++){
618 Ke->values[i*numnodes+j]+=D_scalar*
619 (sqrt(tau_parameter_hor)*(u-um)*dbasis[0*numnodes+i]+sqrt(tau_parameter_hor)*(v-vm)*dbasis[1*numnodes+i]+sqrt(tau_parameter_ver)*(w-wm)*dbasis[2*numnodes+i])*
620 (sqrt(tau_parameter_hor)*(u-um)*dbasis[0*numnodes+j]+sqrt(tau_parameter_hor)*(v-vm)*dbasis[1*numnodes+j]+sqrt(tau_parameter_ver)*(w-wm)*dbasis[2*numnodes+j]);
621 }
622 }
623 }
624 }
625
626 /*Clean up and return*/
627 xDelete<IssmDouble>(xyz_list);
628 xDelete<IssmDouble>(basis);
629 xDelete<IssmDouble>(dbasis);
630 delete gauss;
631 return Ke;
632}/*}}}*/
633ElementMatrix* EnthalpyAnalysis::CreateKMatrixShelf(Element* element){/*{{{*/
634
635 /* Check if ice in element */
636 if(!element->IsIceInElement()) return NULL;
637
638 /*Initialize Element matrix and return if necessary*/
639 if(!element->IsOnBase() || !element->IsAllFloating()) return NULL;
640
641 /*Intermediaries*/
642 IssmDouble dt,Jdet,D;
643 IssmDouble *xyz_list_base = NULL;
644
645 /*Fetch number of nodes for this finite element*/
646 int numnodes = element->GetNumberOfNodes();
647
648 /*Initialize vectors*/
649 ElementMatrix* Ke = element->NewElementMatrix();
650 IssmDouble* basis = xNew<IssmDouble>(numnodes);
651
652 /*Retrieve all inputs and parameters*/
653 element->GetVerticesCoordinatesBase(&xyz_list_base);
654 element->FindParam(&dt,TimesteppingTimeStepEnum);
655 IssmDouble gravity = element->FindParam(ConstantsGEnum);
656 IssmDouble rho_water = element->FindParam(MaterialsRhoSeawaterEnum);
657 IssmDouble rho_ice = element->FindParam(MaterialsRhoIceEnum);
658 IssmDouble heatcapacity = element->FindParam(MaterialsHeatcapacityEnum);
659 IssmDouble mixed_layer_capacity= element->FindParam(MaterialsMixedLayerCapacityEnum);
660 IssmDouble thermal_exchange_vel= element->FindParam(MaterialsThermalExchangeVelocityEnum);
661
662 /* Start looping on the number of gaussian points: */
663 Gauss* gauss=element->NewGaussBase(4);
664 while(gauss->next()){
665
666 element->JacobianDeterminantBase(&Jdet,xyz_list_base,gauss);
667 element->NodalFunctions(basis,gauss);
668
669 D=gauss->weight*Jdet*rho_water*mixed_layer_capacity*thermal_exchange_vel/(heatcapacity*rho_ice);
670 if(reCast<bool,IssmDouble>(dt)) D=dt*D;
671 for(int i=0;i<numnodes;i++) for(int j=0;j<numnodes;j++) Ke->values[i*numnodes+j] += D*basis[i]*basis[j];
672 }
673
674 /*Clean up and return*/
675 delete gauss;
676 xDelete<IssmDouble>(basis);
677 xDelete<IssmDouble>(xyz_list_base);
678 return Ke;
679}/*}}}*/
680ElementVector* EnthalpyAnalysis::CreatePVector(Element* element){/*{{{*/
681
682 /* Check if ice in element */
683 if(!element->IsIceInElement()) return NULL;
684
685 /*compute all load vectors for this element*/
686 ElementVector* pe1=CreatePVectorVolume(element);
687 ElementVector* pe2=CreatePVectorSheet(element);
688 ElementVector* pe3=CreatePVectorShelf(element);
689 ElementVector* pe =new ElementVector(pe1,pe2,pe3);
690
691 /*clean-up and return*/
692 delete pe1;
693 delete pe2;
694 delete pe3;
695 return pe;
696}/*}}}*/
697ElementVector* EnthalpyAnalysis::CreatePVectorVolume(Element* element){/*{{{*/
698
699 /* Check if ice in element */
700 if(!element->IsIceInElement()) return NULL;
701
702 /*Intermediaries*/
703 int i, stabilization;
704 IssmDouble Jdet,phi,dt;
705 IssmDouble enthalpy, Hpmp;
706 IssmDouble enthalpypicard, d1enthalpypicard[3];
707 IssmDouble pressure, d1pressure[3], d2pressure;
708 IssmDouble waterfractionpicard;
709 IssmDouble kappa,tau_parameter,diameter,hx,hy,hz,kappa_w;
710 IssmDouble tau_parameter_anisotropic[2],tau_parameter_hor,tau_parameter_ver;
711 IssmDouble u,v,w;
712 IssmDouble scalar_def, scalar_sens ,scalar_transient;
713 IssmDouble* xyz_list = NULL;
714 IssmDouble d1H_d1P, d1P2;
715
716 /*Fetch number of nodes and dof for this finite element*/
717 int numnodes = element->GetNumberOfNodes();
718 int numvertices = element->GetNumberOfVertices();
719
720 /*Initialize Element vector*/
721 ElementVector* pe = element->NewElementVector();
722 IssmDouble* basis = xNew<IssmDouble>(numnodes);
723 IssmDouble* dbasis = xNew<IssmDouble>(3*numnodes);
724
725 /*Retrieve all inputs and parameters*/
726 element->GetVerticesCoordinates(&xyz_list);
727 IssmDouble rho_ice = element->FindParam(MaterialsRhoIceEnum);
728 IssmDouble heatcapacity = element->FindParam(MaterialsHeatcapacityEnum);
729 IssmDouble thermalconductivity = element->FindParam(MaterialsThermalconductivityEnum);
730 IssmDouble temperateiceconductivity = element->FindParam(MaterialsTemperateiceconductivityEnum);
731 IssmDouble beta = element->FindParam(MaterialsBetaEnum);
732 IssmDouble latentheat = element->FindParam(MaterialsLatentheatEnum);
733 element->FindParam(&dt,TimesteppingTimeStepEnum);
734 element->FindParam(&stabilization,ThermalStabilizationEnum);
735 Input* vx_input=element->GetInput(VxEnum); _assert_(vx_input);
736 Input* vy_input=element->GetInput(VyEnum); _assert_(vy_input);
737 Input* vz_input=element->GetInput(VzEnum); _assert_(vz_input);
738 Input* enthalpypicard_input=element->GetInput(EnthalpyPicardEnum); _assert_(enthalpypicard_input);
739 Input* pressure_input=element->GetInput(PressureEnum); _assert_(pressure_input);
740 Input* enthalpy_input=NULL;
741 if(dt>0.){
742 enthalpy_input = element->GetInput(EnthalpyEnum); _assert_(enthalpy_input);
743 }
744
745 /* Start looping on the number of gaussian points: */
746 Gauss* gauss=element->NewGauss(4);
747 while(gauss->next()){
748
749 element->JacobianDeterminant(&Jdet,xyz_list,gauss);
750 element->NodalFunctions(basis,gauss);
751
752 /*viscous dissipation*/
753 element->ViscousHeating(&phi,xyz_list,gauss,vx_input,vy_input,vz_input);
754
755 scalar_def=phi/rho_ice*Jdet*gauss->weight;
756 if(dt!=0.) scalar_def=scalar_def*dt;
757
758 for(i=0;i<numnodes;i++) pe->values[i]+=scalar_def*basis[i];
759
760 /*sensible heat flux in temperate ice*/
761 enthalpypicard_input->GetInputValue(&enthalpypicard,gauss);
762 pressure_input->GetInputValue(&pressure,gauss);
763 Hpmp=this->PureIceEnthalpy(element, pressure);
764
765 if(enthalpypicard>=Hpmp){
766 enthalpypicard_input->GetInputDerivativeValue(&d1enthalpypicard[0],xyz_list,gauss);
767 pressure_input->GetInputDerivativeValue(&d1pressure[0],xyz_list,gauss);
768 d2pressure=0.; // for linear elements, 2nd derivative is zero
769
770 d1H_d1P=0.;
771 for(i=0;i<3;i++) d1H_d1P+=d1enthalpypicard[i]*d1pressure[i];
772 d1P2=0.;
773 for(i=0;i<3;i++) d1P2+=pow(d1pressure[i],2.);
774
775 scalar_sens=-beta*((temperateiceconductivity - thermalconductivity)/latentheat*(d1H_d1P + beta*heatcapacity*d1P2))/rho_ice;
776 if(dt!=0.) scalar_sens=scalar_sens*dt;
777 for(i=0;i<numnodes;i++) pe->values[i]+=scalar_sens*basis[i];
778 }
779
780 /* Build transient now */
781 if(dt>0.){
782 enthalpy_input->GetInputValue(&enthalpy, gauss);
783 scalar_transient=enthalpy*Jdet*gauss->weight;
784 for(i=0;i<numnodes;i++) pe->values[i]+=scalar_transient*basis[i];
785 }
786
787 /* SUPG */
788 if(stabilization==2){
789 element->NodalFunctionsDerivatives(dbasis,xyz_list,gauss);
790 diameter=element->MinEdgeLength(xyz_list);
791 kappa=this->EnthalpyDiffusionParameterVolume(element,EnthalpyPicardEnum); _assert_(kappa>=0.);
792 vx_input->GetInputValue(&u,gauss);
793 vy_input->GetInputValue(&v,gauss);
794 vz_input->GetInputValue(&w,gauss);
795 tau_parameter=element->StabilizationParameter(u,v,w,diameter,kappa/rho_ice);
796
797 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]);
798
799 if(dt!=0.){
800 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]);
801 }
802 }
803 /* anisotropic SUPG */
804 else if(stabilization==3){
805 element->NodalFunctionsDerivatives(dbasis,xyz_list,gauss);
806 element->ElementSizes(&hx,&hy,&hz);
807 kappa=this->EnthalpyDiffusionParameterVolume(element,EnthalpyPicardEnum); _assert_(kappa>=0.);
808 vx_input->GetInputValue(&u,gauss);
809 vy_input->GetInputValue(&v,gauss);
810 vz_input->GetInputValue(&w,gauss);
811 element->StabilizationParameterAnisotropic(&tau_parameter_anisotropic[0],u,v,w,hx,hy,hz,kappa/rho_ice);
812 tau_parameter_hor=tau_parameter_anisotropic[0];
813 tau_parameter_ver=tau_parameter_anisotropic[1];
814
815 for(i=0;i<numnodes;i++) pe->values[i]+=scalar_def*(tau_parameter_hor*u*dbasis[0*numnodes+i]+tau_parameter_hor*v*dbasis[1*numnodes+i]+tau_parameter_ver*w*dbasis[2*numnodes+i]);
816 }
817 }
818
819 /*Clean up and return*/
820 xDelete<IssmDouble>(basis);
821 xDelete<IssmDouble>(dbasis);
822 xDelete<IssmDouble>(xyz_list);
823 delete gauss;
824 return pe;
825
826}/*}}}*/
827ElementVector* EnthalpyAnalysis::CreatePVectorSheet(Element* element){/*{{{*/
828
829 /* Check if ice in element */
830 if(!element->IsIceInElement()) return NULL;
831
832 /* implementation of the basal condition decision chart of Aschwanden 2012, Fig.5 */
833 if(!element->IsOnBase() || element->IsAllFloating()) return NULL;
834
835 bool converged, isdynamicbasalspc;
836 int i, state;
837 int enthalpy_enum;
838 IssmDouble dt,Jdet,scalar;
839 IssmDouble enthalpy, enthalpyup, pressure, pressureup, watercolumn, meltingrate;
840 IssmDouble vx,vy,vz;
841 IssmDouble alpha2,basalfriction,geothermalflux,heatflux;
842 IssmDouble *xyz_list_base = NULL;
843
844 /*Fetch number of nodes for this finite element*/
845 int numnodes = element->GetNumberOfNodes();
846
847 /*Initialize vectors*/
848 ElementVector* pe = element->NewElementVector();
849 IssmDouble* basis = xNew<IssmDouble>(numnodes);
850
851 /*Retrieve all inputs and parameters*/
852 element->GetVerticesCoordinatesBase(&xyz_list_base);
853 element->FindParam(&dt,TimesteppingTimeStepEnum);
854 element->FindParam(&isdynamicbasalspc,ThermalIsdynamicbasalspcEnum);
855 element->GetInputValue(&converged,ConvergedEnum);
856 if(dt==0. && !converged) enthalpy_enum=EnthalpyPicardEnum; // use enthalpy from last iteration
857 else enthalpy_enum=EnthalpyEnum; // use enthalpy from last time step
858 Input* enthalpy_input = element->GetInput(enthalpy_enum); _assert_(enthalpy_input);
859 Input* pressure_input = element->GetInput(PressureEnum); _assert_(pressure_input);
860 Input* watercolumn_input = element->GetInput(WatercolumnEnum); _assert_(watercolumn_input);
861 Input* meltingrate_input = element->GetInput(BasalforcingsGroundediceMeltingRateEnum); _assert_(meltingrate_input);
862 Input* geothermalflux_input = element->GetInput(BasalforcingsGeothermalfluxEnum); _assert_(geothermalflux_input);
863 IssmDouble rho_ice = element->FindParam(MaterialsRhoIceEnum);
864
865 /*Build friction element, needed later: */
866 Friction* friction=new Friction(element,3);
867
868 /* Start looping on the number of gaussian points: */
869 Gauss* gauss=element->NewGaussBase(4);
870 Gauss* gaussup=element->NewGaussTop(4);
871 while(gauss->next() && gaussup->next()){
872
873 element->JacobianDeterminantBase(&Jdet,xyz_list_base,gauss);
874 element->NodalFunctions(basis,gauss);
875
876 if(isdynamicbasalspc){
877 enthalpy_input->GetInputValue(&enthalpy,gauss);
878 enthalpy_input->GetInputValue(&enthalpyup,gaussup);
879 pressure_input->GetInputValue(&pressure,gauss);
880 pressure_input->GetInputValue(&pressureup,gaussup);
881 watercolumn_input->GetInputValue(&watercolumn,gauss);
882 meltingrate_input->GetInputValue(&meltingrate,gauss);
883 state=GetThermalBasalCondition(element, enthalpy, enthalpyup, pressure, pressureup, watercolumn, meltingrate);
884 }
885 else
886 state=0;
887
888 switch (state) {
889 case 0: case 1: case 2: case 3:
890 // cold, dry base; cold, wet base; refreezing temperate base; thin temperate base:
891 // Apply basal surface forcing.
892 // Interpolated values of enthalpy on gauss nodes may indicate cold base,
893 // although one node might have become temperate. So keep heat flux switched on.
894 geothermalflux_input->GetInputValue(&geothermalflux,gauss);
895 friction->GetAlpha2(&alpha2,gauss);
896 friction->GetBasalSlidingSpeeds(&vx, &vy, &vz, gauss);
897 basalfriction=alpha2*(vx*vx+vy*vy+vz*vz);
898 heatflux=(basalfriction+geothermalflux)/(rho_ice);
899 scalar=gauss->weight*Jdet*heatflux;
900 if(dt!=0.) scalar=dt*scalar;
901 for(i=0;i<numnodes;i++)
902 pe->values[i]+=scalar*basis[i];
903 break;
904 case 4:
905 // temperate, thick melting base: set grad H*n=0
906 for(i=0;i<numnodes;i++)
907 pe->values[i]+=0.;
908 break;
909 default:
910 _printf0_(" unknown thermal basal state found!");
911 }
912 }
913
914 /*Clean up and return*/
915 delete gauss;
916 delete gaussup;
917 delete friction;
918 xDelete<IssmDouble>(basis);
919 xDelete<IssmDouble>(xyz_list_base);
920 return pe;
921
922}/*}}}*/
923ElementVector* EnthalpyAnalysis::CreatePVectorShelf(Element* element){/*{{{*/
924
925 /* Check if ice in element */
926 if(!element->IsIceInElement()) return NULL;
927
928 /*Get basal element*/
929 if(!element->IsOnBase() || !element->IsAllFloating()) return NULL;
930
931 IssmDouble Hpmp,dt,Jdet,scalar_ocean,pressure;
932 IssmDouble *xyz_list_base = NULL;
933
934 /*Fetch number of nodes for this finite element*/
935 int numnodes = element->GetNumberOfNodes();
936
937 /*Initialize vectors*/
938 ElementVector* pe = element->NewElementVector();
939 IssmDouble* basis = xNew<IssmDouble>(numnodes);
940
941 /*Retrieve all inputs and parameters*/
942 element->GetVerticesCoordinatesBase(&xyz_list_base);
943 element->FindParam(&dt,TimesteppingTimeStepEnum);
944 Input* pressure_input=element->GetInput(PressureEnum); _assert_(pressure_input);
945 IssmDouble gravity = element->FindParam(ConstantsGEnum);
946 IssmDouble rho_water = element->FindParam(MaterialsRhoSeawaterEnum);
947 IssmDouble rho_ice = element->FindParam(MaterialsRhoIceEnum);
948 IssmDouble heatcapacity = element->FindParam(MaterialsHeatcapacityEnum);
949 IssmDouble mixed_layer_capacity= element->FindParam(MaterialsMixedLayerCapacityEnum);
950 IssmDouble thermal_exchange_vel= element->FindParam(MaterialsThermalExchangeVelocityEnum);
951
952 /* Start looping on the number of gaussian points: */
953 Gauss* gauss=element->NewGaussBase(4);
954 while(gauss->next()){
955
956 element->JacobianDeterminantBase(&Jdet,xyz_list_base,gauss);
957 element->NodalFunctions(basis,gauss);
958
959 pressure_input->GetInputValue(&pressure,gauss);
960 Hpmp=element->PureIceEnthalpy(pressure);
961
962 scalar_ocean=gauss->weight*Jdet*rho_water*mixed_layer_capacity*thermal_exchange_vel*Hpmp/(heatcapacity*rho_ice);
963 if(reCast<bool,IssmDouble>(dt)) scalar_ocean=dt*scalar_ocean;
964
965 for(int i=0;i<numnodes;i++) pe->values[i]+=scalar_ocean*basis[i];
966 }
967
968 /*Clean up and return*/
969 delete gauss;
970 xDelete<IssmDouble>(basis);
971 xDelete<IssmDouble>(xyz_list_base);
972 return pe;
973}/*}}}*/
974void EnthalpyAnalysis::DrainWaterfraction(FemModel* femmodel){/*{{{*/
975 /*Drain excess water fraction in ice column: */
976 ComputeWaterfractionDrainage(femmodel);
977 DrainageUpdateWatercolumn(femmodel);
978 DrainageUpdateEnthalpy(femmodel);
979}/*}}}*/
980void EnthalpyAnalysis::ComputeWaterfractionDrainage(FemModel* femmodel){/*{{{*/
981
982 int k,numnodes;
983 IssmDouble dt;
984 Element* element= NULL;
985
986 femmodel->parameters->FindParam(&dt,TimesteppingTimeStepEnum);
987
988 for(Object* & object : femmodel->elements->objects){
989 element=xDynamicCast<Element*>(object);
990 numnodes=element->GetNumberOfNodes();
991 IssmDouble* waterfractions= xNew<IssmDouble>(numnodes);
992 IssmDouble* drainage= xNew<IssmDouble>(numnodes);
993
994 element->GetInputListOnNodes(waterfractions,WaterfractionEnum);
995 for(k=0; k<numnodes;k++){
996 drainage[k]=DrainageFunctionWaterfraction(waterfractions[k], dt);
997 }
998 int finite_element = element->GetElementType(); if(finite_element==P1Enum) finite_element = P1DGEnum;
999 element->AddInput(WaterfractionDrainageEnum,drainage,finite_element);
1000
1001 xDelete<IssmDouble>(waterfractions);
1002 xDelete<IssmDouble>(drainage);
1003 }
1004}/*}}}*/
1005void EnthalpyAnalysis::DrainageUpdateWatercolumn(FemModel* femmodel){/*{{{*/
1006
1007 int k,numnodes, numbasalnodes;
1008 IssmDouble dt;
1009 int* basalnodeindices=NULL;
1010 Element* element= NULL;
1011
1012 femmodel->parameters->FindParam(&dt,TimesteppingTimeStepEnum);
1013
1014 /*depth-integrate the drained water fraction */
1015 femmodel->parameters->SetParam(WaterfractionDrainageEnum,InputToDepthaverageInEnum);
1016 femmodel->parameters->SetParam(WaterfractionDrainageIntegratedEnum,InputToDepthaverageOutEnum);
1017 depthaverage_core(femmodel);
1018 femmodel->parameters->SetParam(WaterfractionDrainageIntegratedEnum,InputToExtrudeEnum);
1019 extrudefrombase_core(femmodel);
1020 /*multiply depth-average by ice thickness*/
1021 for(Object* & object : femmodel->elements->objects){
1022 element=xDynamicCast<Element*>(object);
1023 numnodes=element->GetNumberOfNodes();
1024 IssmDouble* drainage_int= xNew<IssmDouble>(numnodes);
1025 IssmDouble* thicknesses= xNew<IssmDouble>(numnodes);
1026
1027 element->GetInputListOnNodes(drainage_int,WaterfractionDrainageIntegratedEnum);
1028 element->GetInputListOnNodes(thicknesses,ThicknessEnum);
1029 for(k=0;k<numnodes;k++){
1030 drainage_int[k]*=thicknesses[k];
1031 }
1032 int finite_element = element->GetElementType(); if(finite_element==P1Enum) finite_element = P1DGEnum;
1033 element->AddInput(WaterfractionDrainageIntegratedEnum, drainage_int,finite_element);
1034
1035 xDelete<IssmDouble>(drainage_int);
1036 xDelete<IssmDouble>(thicknesses);
1037 }
1038
1039 /*update water column*/
1040 for(Object* & object : femmodel->elements->objects){
1041 element=xDynamicCast<Element*>(object);
1042 /* Check if ice in element */
1043 if(!element->IsIceInElement()) continue;
1044 if(!element->IsOnBase()) continue;
1045
1046 numnodes=element->GetNumberOfNodes();
1047 IssmDouble* watercolumn= xNew<IssmDouble>(numnodes);
1048 IssmDouble* drainage_int= xNew<IssmDouble>(numnodes);
1049 element->GetInputListOnNodes(watercolumn,WatercolumnEnum);
1050 element->GetInputListOnNodes(drainage_int,WaterfractionDrainageIntegratedEnum);
1051
1052 element->BasalNodeIndices(&numbasalnodes,&basalnodeindices,element->GetElementType());
1053 for(k=0;k<numbasalnodes;k++){
1054 watercolumn[basalnodeindices[k]]+=dt*drainage_int[basalnodeindices[k]];
1055 }
1056 int finite_element = element->GetElementType(); if(finite_element==P1Enum) finite_element = P1DGEnum;
1057 element->AddInput(WatercolumnEnum, watercolumn,finite_element);
1058
1059 xDelete<IssmDouble>(watercolumn);
1060 xDelete<IssmDouble>(drainage_int);
1061 xDelete<int>(basalnodeindices);
1062 }
1063}/*}}}*/
1064void EnthalpyAnalysis::DrainageUpdateEnthalpy(FemModel* femmodel){/*{{{*/
1065
1066 int k,numnodes;
1067 IssmDouble dt;
1068 Element* element= NULL;
1069 femmodel->parameters->FindParam(&dt,TimesteppingTimeStepEnum);
1070
1071 for(Object* & object : femmodel->elements->objects){
1072 element=xDynamicCast<Element*>(object);
1073 numnodes=element->GetNumberOfNodes();
1074 IssmDouble* enthalpies= xNew<IssmDouble>(numnodes);
1075 IssmDouble* pressures= xNew<IssmDouble>(numnodes);
1076 IssmDouble* temperatures= xNew<IssmDouble>(numnodes);
1077 IssmDouble* waterfractions= xNew<IssmDouble>(numnodes);
1078 IssmDouble* drainage= xNew<IssmDouble>(numnodes);
1079
1080 element->GetInputListOnNodes(pressures,PressureEnum);
1081 element->GetInputListOnNodes(temperatures,TemperatureEnum);
1082 element->GetInputListOnNodes(waterfractions,WaterfractionEnum);
1083 element->GetInputListOnNodes(drainage,WaterfractionDrainageEnum);
1084
1085 for(k=0;k<numnodes;k++){
1086 if(dt==0.)
1087 waterfractions[k]-=drainage[k];
1088 else
1089 waterfractions[k]-=dt*drainage[k];
1090
1091 element->ThermalToEnthalpy(&enthalpies[k], temperatures[k], waterfractions[k], pressures[k]);
1092 }
1093 int finite_element = element->GetElementType(); if(finite_element==P1Enum) finite_element = P1DGEnum;
1094 element->AddInput(WaterfractionEnum,waterfractions,finite_element);
1095 element->AddInput(EnthalpyEnum,enthalpies,finite_element);
1096
1097 xDelete<IssmDouble>(enthalpies);
1098 xDelete<IssmDouble>(pressures);
1099 xDelete<IssmDouble>(temperatures);
1100 xDelete<IssmDouble>(waterfractions);
1101 xDelete<IssmDouble>(drainage);
1102 }
1103}/*}}}*/
1104IssmDouble EnthalpyAnalysis::EnthalpyDiffusionParameter(Element* element,IssmDouble enthalpy,IssmDouble pressure){/*{{{*/
1105
1106 IssmDouble heatcapacity = element->FindParam(MaterialsHeatcapacityEnum);
1107 IssmDouble temperateiceconductivity = element->FindParam(MaterialsTemperateiceconductivityEnum);
1108 IssmDouble thermalconductivity = element->FindParam(MaterialsThermalconductivityEnum);
1109
1110 if(enthalpy < PureIceEnthalpy(element,pressure))
1111 return thermalconductivity/heatcapacity;
1112 else
1113 return temperateiceconductivity/heatcapacity;
1114}/*}}}*/
1115IssmDouble EnthalpyAnalysis::EnthalpyDiffusionParameterVolume(Element* element,int enthalpy_enum){/*{{{*/
1116
1117 int iv;
1118 IssmDouble lambda; /* fraction of cold ice */
1119 IssmDouble kappa,kappa_c,kappa_t; /* enthalpy conductivities */
1120 IssmDouble Hc,Ht;
1121
1122 /*Get pressures and enthalpies on vertices*/
1123 int numvertices = element->GetNumberOfVertices();
1124 int effectiveconductivity_averaging;
1125 IssmDouble* pressures = xNew<IssmDouble>(numvertices);
1126 IssmDouble* enthalpies = xNew<IssmDouble>(numvertices);
1127 IssmDouble* PIE = xNew<IssmDouble>(numvertices);
1128 IssmDouble* dHpmp = xNew<IssmDouble>(numvertices);
1129 element->GetInputListOnVertices(pressures,PressureEnum);
1130 element->GetInputListOnVertices(enthalpies,enthalpy_enum);
1131 element->FindParam(&effectiveconductivity_averaging,MaterialsEffectiveconductivityAveragingEnum);
1132
1133 for(iv=0;iv<numvertices;iv++){
1134 PIE[iv] = PureIceEnthalpy(element,pressures[iv]);
1135 dHpmp[iv] = enthalpies[iv]-PIE[iv];
1136 }
1137
1138 bool allequalsign = true;
1139 if(dHpmp[0]<0.){
1140 for(iv=1; iv<numvertices;iv++) allequalsign=(allequalsign && (dHpmp[iv]<0.));
1141 }
1142 else{
1143 for(iv=1; iv<numvertices;iv++) allequalsign=(allequalsign && (dHpmp[iv]>=0.));
1144 }
1145
1146 if(allequalsign){
1147 kappa = EnthalpyDiffusionParameter(element,enthalpies[0],pressures[0]);
1148 }
1149 else{
1150 kappa_c = EnthalpyDiffusionParameter(element,PureIceEnthalpy(element,0.)-1.,0.);
1151 kappa_t = EnthalpyDiffusionParameter(element,PureIceEnthalpy(element,0.)+1.,0.);
1152
1153 Hc=0.; Ht=0.;
1154 for(iv=0; iv<numvertices;iv++){
1155 if(enthalpies[iv]<PIE[iv])
1156 Hc+=(PIE[iv]-enthalpies[iv]);
1157 else
1158 Ht+=(enthalpies[iv]-PIE[iv]);
1159 }
1160 _assert_((Hc+Ht)>0.);
1161 lambda = Hc/(Hc+Ht);
1162 _assert_(lambda>=0.);
1163 _assert_(lambda<=1.);
1164
1165 if(effectiveconductivity_averaging==0){
1166 /* return arithmetic mean (volume average) of thermal conductivities, weighted by fraction of cold/temperate ice */
1167 kappa=kappa_c*lambda+(1.-lambda)*kappa_t;
1168 }
1169 else if(effectiveconductivity_averaging==1){
1170 /* return harmonic mean (reciprocal avarage) of thermal conductivities, weighted by fraction of cold/temperate ice, cf Patankar 1980, pp44 */
1171 kappa=kappa_c*kappa_t/(lambda*kappa_t+(1.-lambda)*kappa_c);
1172 }
1173 else if(effectiveconductivity_averaging==2){
1174 /* return geometric mean (power law) of thermal conductivities, weighted by fraction of cold/temperate ice */
1175 kappa=pow(kappa_c,lambda)*pow(kappa_t,1.-lambda);
1176 }
1177 else{
1178 _error_("effectiveconductivity_averaging not supported yet");
1179 }
1180 }
1181
1182 /*Clean up and return*/
1183 xDelete<IssmDouble>(PIE);
1184 xDelete<IssmDouble>(dHpmp);
1185 xDelete<IssmDouble>(pressures);
1186 xDelete<IssmDouble>(enthalpies);
1187 return kappa;
1188}/*}}}*/
1189void EnthalpyAnalysis::GetBasalConstraints(Vector<IssmDouble>* vec_spc,Element* element){/*{{{*/
1190
1191 /*Intermediary*/
1192 bool isdynamicbasalspc;
1193 IssmDouble dt;
1194
1195 /*Check wether dynamic basal boundary conditions are activated */
1196 element->FindParam(&isdynamicbasalspc,ThermalIsdynamicbasalspcEnum);
1197 if(!isdynamicbasalspc) return;
1198
1199 element->FindParam(&dt,TimesteppingTimeStepEnum);
1200 if(dt==0.){
1201 GetBasalConstraintsSteadystate(vec_spc,element);
1202 }
1203 else{
1204 GetBasalConstraintsTransient(vec_spc,element);
1205 }
1206}/*}}}*/
1207void EnthalpyAnalysis::GetBasalConstraintsSteadystate(Vector<IssmDouble>* vec_spc,Element* element){/*{{{*/
1208
1209 /* Check if ice in element */
1210 if(!element->IsIceInElement()) return;
1211
1212 /* Only update constraints at the base.
1213 * Floating ice is not affected by basal BC decision chart. */
1214 if(!(element->IsOnBase()) || element->IsAllFloating()) return;
1215
1216 /*Intermediary*/
1217 int numindices, numindicesup, state;
1218 int *indices = NULL, *indicesup = NULL;
1219 IssmDouble enthalpy, enthalpyup, pressure, pressureup, watercolumn, meltingrate;
1220
1221 /*Get parameters and inputs: */
1222 Input* enthalpy_input = element->GetInput(EnthalpyPicardEnum); _assert_(enthalpy_input);
1223 Input* pressure_input = element->GetInput(PressureEnum); _assert_(pressure_input);
1224 Input* watercolumn_input = element->GetInput(WatercolumnEnum); _assert_(watercolumn_input);
1225 Input* meltingrate_input = element->GetInput(BasalforcingsGroundediceMeltingRateEnum); _assert_(meltingrate_input);
1226
1227 /*Fetch indices of basal & surface nodes for this finite element*/
1228 Penta *penta = (Penta *) element; // TODO: add Basal-/SurfaceNodeIndices to element.h, and change this to Element*
1229 penta->BasalNodeIndices(&numindices,&indices,element->GetElementType());
1230 penta->SurfaceNodeIndices(&numindicesup,&indicesup,element->GetElementType()); _assert_(numindices==numindicesup);
1231
1232 GaussPenta* gauss=new GaussPenta();
1233 GaussPenta* gaussup=new GaussPenta();
1234 for(int i=0;i<numindices;i++){
1235 gauss->GaussNode(element->GetElementType(),indices[i]);
1236 gaussup->GaussNode(element->GetElementType(),indicesup[i]);
1237
1238 enthalpy_input->GetInputValue(&enthalpy,gauss);
1239 enthalpy_input->GetInputValue(&enthalpyup,gaussup);
1240 pressure_input->GetInputValue(&pressure,gauss);
1241 pressure_input->GetInputValue(&pressureup,gaussup);
1242 watercolumn_input->GetInputValue(&watercolumn,gauss);
1243 meltingrate_input->GetInputValue(&meltingrate,gauss);
1244
1245 state=GetThermalBasalCondition(element, enthalpy, enthalpyup, pressure, pressureup, watercolumn, meltingrate);
1246 switch (state) {
1247 case 0:
1248 // cold, dry base: apply basal surface forcing
1249 vec_spc->SetValue(element->nodes[i]->Pid(),0.,INS_VAL);
1250 break;
1251 case 1:
1252 // cold, wet base: keep at pressure melting point
1253 vec_spc->SetValue(element->nodes[i]->Pid(),1.,INS_VAL);
1254 break;
1255 case 2:
1256 // temperate, thin refreezing base:
1257 vec_spc->SetValue(element->nodes[i]->Pid(),1.,INS_VAL);
1258 break;
1259 case 3:
1260 // temperate, thin melting base: set spc
1261 vec_spc->SetValue(element->nodes[i]->Pid(),1.,INS_VAL);
1262 break;
1263 case 4:
1264 // temperate, thick melting base:
1265 vec_spc->SetValue(element->nodes[i]->Pid(),1.,INS_VAL);
1266 break;
1267 default:
1268 _printf0_(" unknown thermal basal state found!");
1269 }
1270 }
1271
1272 /*Free resources:*/
1273 xDelete<int>(indices);
1274 xDelete<int>(indicesup);
1275 delete gauss;
1276 delete gaussup;
1277}/*}}}*/
1278void EnthalpyAnalysis::GetBasalConstraintsTransient(Vector<IssmDouble>* vec_spc,Element* element){/*{{{*/
1279
1280 /* Check if ice in element */
1281 if(!element->IsIceInElement()) return;
1282
1283 /* Only update constraints at the base.
1284 * Floating ice is not affected by basal BC decision chart.*/
1285 if(!(element->IsOnBase()) || element->IsAllFloating()) return;
1286
1287 /*Intermediary*/
1288 int numindices, numindicesup, state;
1289 int *indices = NULL, *indicesup = NULL;
1290 IssmDouble enthalpy, enthalpyup, pressure, pressureup, watercolumn, meltingrate;
1291
1292 /*Get parameters and inputs: */
1293 Input* enthalpy_input = element->GetInput(EnthalpyEnum); _assert_(enthalpy_input); //TODO: check EnthalpyPicard?
1294 Input* pressure_input = element->GetInput(PressureEnum); _assert_(pressure_input);
1295 Input* watercolumn_input = element->GetInput(WatercolumnEnum); _assert_(watercolumn_input);
1296 Input* meltingrate_input = element->GetInput(BasalforcingsGroundediceMeltingRateEnum); _assert_(meltingrate_input);
1297
1298 /*Fetch indices of basal & surface nodes for this finite element*/
1299 Penta *penta = (Penta *) element; // TODO: add Basal-/SurfaceNodeIndices to element.h, and change this to Element*
1300 penta->BasalNodeIndices(&numindices,&indices,element->GetElementType());
1301 penta->SurfaceNodeIndices(&numindicesup,&indicesup,element->GetElementType()); _assert_(numindices==numindicesup);
1302
1303 GaussPenta* gauss=new GaussPenta();
1304 GaussPenta* gaussup=new GaussPenta();
1305
1306 for(int i=0;i<numindices;i++){
1307 gauss->GaussNode(element->GetElementType(),indices[i]);
1308 gaussup->GaussNode(element->GetElementType(),indicesup[i]);
1309
1310 enthalpy_input->GetInputValue(&enthalpy,gauss);
1311 enthalpy_input->GetInputValue(&enthalpyup,gaussup);
1312 pressure_input->GetInputValue(&pressure,gauss);
1313 pressure_input->GetInputValue(&pressureup,gaussup);
1314 watercolumn_input->GetInputValue(&watercolumn,gauss);
1315 meltingrate_input->GetInputValue(&meltingrate,gauss);
1316
1317 state=GetThermalBasalCondition(element, enthalpy, enthalpyup, pressure, pressureup, watercolumn, meltingrate);
1318
1319 switch (state) {
1320 case 0:
1321 // cold, dry base: apply basal surface forcing
1322 vec_spc->SetValue(element->nodes[i]->Pid(),0.,INS_VAL);
1323 break;
1324 case 1:
1325 // cold, wet base: keep at pressure melting point
1326 vec_spc->SetValue(element->nodes[i]->Pid(),1.,INS_VAL);
1327 break;
1328 case 2:
1329 // temperate, thin refreezing base: release spc
1330 vec_spc->SetValue(element->nodes[i]->Pid(),0.,INS_VAL);
1331 break;
1332 case 3:
1333 // temperate, thin melting base: set spc
1334 vec_spc->SetValue(element->nodes[i]->Pid(),1.,INS_VAL);
1335 break;
1336 case 4:
1337 // temperate, thick melting base: set grad H*n=0
1338 vec_spc->SetValue(element->nodes[i]->Pid(),0.,INS_VAL);
1339 break;
1340 default:
1341 _printf0_(" unknown thermal basal state found!");
1342 }
1343
1344 }
1345
1346 /*Free resources:*/
1347 xDelete<int>(indices);
1348 xDelete<int>(indicesup);
1349 delete gauss;
1350 delete gaussup;
1351}/*}}}*/
1352void EnthalpyAnalysis::GetBConduct(IssmDouble* B,Element* element,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/
1353 /*Compute B matrix. B=[B1 B2 B3 B4 B5 B6] where Bi is of size 5*1.
1354 * For node i, Bi' can be expressed in the actual coordinate system
1355 * by:
1356 * Bi_conduct=[ dh/dx ]
1357 * [ dh/dy ]
1358 * [ dh/dz ]
1359 * where h is the interpolation function for node i.
1360 *
1361 * We assume B has been allocated already, of size: 3x(1*numnodes)
1362 */
1363
1364 /*Fetch number of nodes for this finite element*/
1365 int numnodes = element->GetNumberOfNodes();
1366
1367 /*Get nodal functions derivatives*/
1368 IssmDouble* dbasis=xNew<IssmDouble>(3*numnodes);
1369 element->NodalFunctionsDerivatives(dbasis,xyz_list,gauss);
1370
1371 /*Build B: */
1372 for(int i=0;i<numnodes;i++){
1373 B[numnodes*0+i] = dbasis[0*numnodes+i];
1374 B[numnodes*1+i] = dbasis[1*numnodes+i];
1375 B[numnodes*2+i] = dbasis[2*numnodes+i];
1376 }
1377
1378 /*Clean-up*/
1379 xDelete<IssmDouble>(dbasis);
1380}/*}}}*/
1381void EnthalpyAnalysis::GetSolutionFromInputs(Vector<IssmDouble>* solution,Element* element){/*{{{*/
1382 element->GetSolutionFromInputsOneDof(solution,EnthalpyEnum);
1383}/*}}}*/
1384int EnthalpyAnalysis::GetThermalBasalCondition(Element* element, IssmDouble enthalpy, IssmDouble enthalpyup, IssmDouble pressure, IssmDouble pressureup, IssmDouble watercolumn, IssmDouble meltingrate){/*{{{*/
1385
1386 /* Check if ice in element */
1387 if(!element->IsIceInElement()) return -1;
1388
1389 /* Only update Constraints at the base of grounded ice*/
1390 if(!(element->IsOnBase())) return -1;
1391
1392 /*Intermediary*/
1393 int state=-1;
1394 IssmDouble dt;
1395
1396 /*Get parameters and inputs: */
1397 element->FindParam(&dt,TimesteppingTimeStepEnum);
1398
1399 if(enthalpy<PureIceEnthalpy(element,pressure)){
1400 if(watercolumn<=0.) state=0; // cold, dry base
1401 else state=1; // cold, wet base (refreezing)
1402 }
1403 else{
1404 if(enthalpyup<PureIceEnthalpy(element,pressureup)){
1405 if((dt==0.) && (meltingrate<0.)) state=2; // refreezing temperate base (non-physical, only for steadystate solver)
1406 else state=3; // temperate base, but no temperate layer
1407 }
1408 else state=4; // temperate layer with positive thickness
1409 }
1410
1411 _assert_(state>=0);
1412 return state;
1413}/*}}}*/
1414IssmDouble EnthalpyAnalysis::GetWetIceConductivity(Element* element, IssmDouble enthalpy, IssmDouble pressure){/*{{{*/
1415
1416 IssmDouble temperature, waterfraction;
1417 IssmDouble kappa_w = 0.6; // thermal conductivity of water (in W/m/K)
1418 IssmDouble kappa_i = element->FindParam(MaterialsThermalconductivityEnum);
1419 element->EnthalpyToThermal(&temperature, &waterfraction, enthalpy, pressure);
1420
1421 return (1.-waterfraction)*kappa_i + waterfraction*kappa_w;
1422}/*}}}*/
1423void EnthalpyAnalysis::GradientJ(Vector<IssmDouble>* gradient,Element* element,int control_type,int control_interp,int control_index){/*{{{*/
1424 _error_("Not implemented yet");
1425}/*}}}*/
1426void EnthalpyAnalysis::InputUpdateFromSolution(IssmDouble* solution,Element* element){/*{{{*/
1427
1428 bool converged;
1429 int i,rheology_law;
1430 IssmDouble B_average,s_average,T_average=0.,P_average=0.;
1431 int *doflist = NULL;
1432 IssmDouble *xyz_list = NULL;
1433
1434 /*Fetch number of nodes and dof for this finite element*/
1435 int numnodes = element->GetNumberOfNodes();
1436
1437 /*Fetch dof list and allocate solution vector*/
1438 element->GetDofListLocal(&doflist,NoneApproximationEnum,GsetEnum);
1439 IssmDouble* values = xNew<IssmDouble>(numnodes);
1440 IssmDouble* pressure = xNew<IssmDouble>(numnodes);
1441 IssmDouble* surface = xNew<IssmDouble>(numnodes);
1442 IssmDouble* B = xNew<IssmDouble>(numnodes);
1443 IssmDouble* temperature = xNew<IssmDouble>(numnodes);
1444 IssmDouble* waterfraction = xNew<IssmDouble>(numnodes);
1445
1446 /*Use the dof list to index into the solution vector: */
1447 for(i=0;i<numnodes;i++){
1448 values[i]=solution[doflist[i]];
1449
1450 /*Check solution*/
1451 if(xIsNan<IssmDouble>(values[i])) _error_("NaN found in solution vector");
1452 if(xIsInf<IssmDouble>(values[i])) _error_("Inf found in solution vector");
1453 }
1454
1455 /*Get all inputs and parameters*/
1456 element->GetInputValue(&converged,ConvergedEnum);
1457 element->GetInputListOnNodes(&pressure[0],PressureEnum);
1458 int finite_element = element->GetElementType(); if(finite_element==P1Enum) finite_element = P1DGEnum;
1459 if(converged){
1460 for(i=0;i<numnodes;i++){
1461 element->EnthalpyToThermal(&temperature[i],&waterfraction[i],values[i],pressure[i]);
1462 if(waterfraction[i]<0.) _error_("Negative water fraction found in solution vector");
1463 //if(waterfraction[i]>1.) _error_("Water fraction >1 found in solution vector");
1464 }
1465 element->AddInput(EnthalpyEnum,values,finite_element);
1466 element->AddInput(WaterfractionEnum,waterfraction,finite_element);
1467 element->AddInput(TemperatureEnum,temperature,finite_element);
1468
1469 IssmDouble* n = xNew<IssmDouble>(numnodes);
1470 if(element->material->ObjectEnum()==MatestarEnum){
1471 for(i=0;i<numnodes;i++) n[i]=3.;
1472 }
1473 else{
1474 element->GetInputListOnNodes(&n[0],MaterialsRheologyNEnum);
1475 }
1476
1477 /*Update Rheology only if converged (we must make sure that the temperature is below melting point
1478 * otherwise the rheology could be negative*/
1479 element->FindParam(&rheology_law,MaterialsRheologyLawEnum);
1480 element->GetInputListOnNodes(&surface[0],SurfaceEnum);
1481 switch(rheology_law){
1482 case NoneEnum:
1483 /*Do nothing: B is not temperature dependent*/
1484 break;
1485 case BuddJackaEnum:
1486 for(i=0;i<numnodes;i++) B[i]=BuddJacka(temperature[i]);
1487 element->AddInput(MaterialsRheologyBEnum,&B[0],finite_element);
1488 break;
1489 case CuffeyEnum:
1490 for(i=0;i<numnodes;i++) B[i]=Cuffey(temperature[i]);
1491 element->AddInput(MaterialsRheologyBEnum,&B[0],finite_element);
1492 break;
1493 case CuffeyTemperateEnum:
1494 for(i=0;i<numnodes;i++) B[i]=CuffeyTemperate(temperature[i], waterfraction[i],n[i]);
1495 element->AddInput(MaterialsRheologyBEnum,&B[0],finite_element);
1496 break;
1497 case PatersonEnum:
1498 for(i=0;i<numnodes;i++) B[i]=Paterson(temperature[i]);
1499 element->AddInput(MaterialsRheologyBEnum,&B[0],finite_element);
1500 break;
1501 case NyeH2OEnum:
1502 for(i=0;i<numnodes;i++) B[i]=NyeH2O(values[i]);
1503 element->AddInput(MaterialsRheologyBEnum,&B[0],finite_element);
1504 break;
1505 case NyeCO2Enum:
1506 for(i=0;i<numnodes;i++) B[i]=NyeCO2(values[i]);
1507 element->AddInput(MaterialsRheologyBEnum,&B[0],finite_element);
1508 break;
1509 case ArrheniusEnum:{
1510 element->GetVerticesCoordinates(&xyz_list);
1511 for(i=0;i<numnodes;i++) B[i]=Arrhenius(temperature[i],surface[i]-xyz_list[i*3+2],n[i]);
1512 element->AddInput(MaterialsRheologyBEnum,&B[0],finite_element);
1513 break;
1514 }
1515 case LliboutryDuvalEnum:{
1516 for(i=0;i<numnodes;i++) B[i]=LliboutryDuval(values[i],pressure[i],n[i],element->FindParam(MaterialsBetaEnum),element->FindParam(ConstantsReferencetemperatureEnum),element->FindParam(MaterialsHeatcapacityEnum),element->FindParam(MaterialsLatentheatEnum));
1517 element->AddInput(MaterialsRheologyBEnum,&B[0],finite_element);
1518 break;
1519 }
1520 default: _error_("Rheology law " << EnumToStringx(rheology_law) << " not supported yet");
1521 }
1522 xDelete<IssmDouble>(n);
1523 }
1524 else{
1525 element->AddInput(EnthalpyPicardEnum,values,finite_element);
1526 }
1527
1528 /*Free resources:*/
1529 xDelete<IssmDouble>(values);
1530 xDelete<IssmDouble>(pressure);
1531 xDelete<IssmDouble>(surface);
1532 xDelete<IssmDouble>(B);
1533 xDelete<IssmDouble>(temperature);
1534 xDelete<IssmDouble>(waterfraction);
1535 xDelete<IssmDouble>(xyz_list);
1536 xDelete<int>(doflist);
1537}/*}}}*/
1538void EnthalpyAnalysis::PostProcessing(FemModel* femmodel){/*{{{*/
1539
1540 /*Intermediaries*/
1541 bool computebasalmeltingrates=true;
1542 bool isdrainicecolumn;
1543 IssmDouble dt;
1544
1545 femmodel->parameters->FindParam(&dt,TimesteppingTimeStepEnum);
1546 femmodel->parameters->FindParam(&isdrainicecolumn,ThermalIsdrainicecolumnEnum);
1547
1548 if(isdrainicecolumn){
1549 DrainWaterfraction(femmodel);
1550 }
1551 if(computebasalmeltingrates){
1552 ComputeBasalMeltingrate(femmodel);
1553 }
1554
1555}/*}}}*/
1556IssmDouble EnthalpyAnalysis::PureIceEnthalpy(Element* element,IssmDouble pressure){/*{{{*/
1557
1558 IssmDouble heatcapacity = element->FindParam(MaterialsHeatcapacityEnum);
1559 IssmDouble referencetemperature = element->FindParam(ConstantsReferencetemperatureEnum);
1560
1561 return heatcapacity*(TMeltingPoint(element,pressure)-referencetemperature);
1562}/*}}}*/
1563IssmDouble EnthalpyAnalysis::TMeltingPoint(Element* element,IssmDouble pressure){/*{{{*/
1564
1565 IssmDouble meltingpoint = element->FindParam(MaterialsMeltingpointEnum);
1566 IssmDouble beta = element->FindParam(MaterialsBetaEnum);
1567
1568 return meltingpoint-beta*pressure;
1569}/*}}}*/
1570void EnthalpyAnalysis::UpdateBasalConstraints(FemModel* femmodel){/*{{{*/
1571
1572 /*Update basal dirichlet BCs for enthalpy: */
1573 int numnodes = femmodel->nodes->NumberOfNodes();
1574 int localmasters = femmodel->nodes->NumberOfNodesLocal();
1575 Vector<IssmDouble>* spc = new Vector<IssmDouble>(localmasters,numnodes);
1576
1577 /*First create a vector to figure out what elements should be constrained*/
1578 for(Object* & object : femmodel->elements->objects){
1579 Element* element=xDynamicCast<Element*>(object);
1580 GetBasalConstraints(spc,element);
1581 }
1582
1583 /*Assemble*/
1584 spc->Assemble();
1585
1586 /*Get local vector with both masters and slaves:*/
1587 IssmDouble *local_spc = NULL;
1588 femmodel->GetLocalVectorWithClonesNodes(&local_spc,spc);
1589 delete spc;
1590
1591 /*Then update basal constraints nodes accordingly*/
1592 for(Object* & object : femmodel->elements->objects){
1593 Element* element=xDynamicCast<Element*>(object);
1594 ApplyBasalConstraints(local_spc,element);
1595 }
1596
1597 femmodel->UpdateConstraintsx();
1598
1599 /*Delete*/
1600 xDelete<IssmDouble>(local_spc);
1601}/*}}}*/
1602void EnthalpyAnalysis::UpdateConstraints(FemModel* femmodel){/*{{{*/
1603 SetActiveNodesLSMx(femmodel);
1604}/*}}}*/
Note: See TracBrowser for help on using the repository browser.