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

Last change on this file since 26463 was 26463, checked in by Cheng Gong, 3 years ago

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