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