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

Last change on this file since 24136 was 24136, checked in by rueckamp, 6 years ago

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