source: issm/trunk/src/c/analyses/HydrologySommersAnalysis.cpp@ 22758

Last change on this file since 22758 was 22758, checked in by Mathieu Morlighem, 7 years ago

merged trunk-jpl and trunk for revision 22757
File size: 23.6 KB
Line 
1#include "./HydrologySommersAnalysis.h"
2#include "../toolkits/toolkits.h"
3#include "../classes/classes.h"
4#include "../shared/shared.h"
5#include "../modules/modules.h"
6
7/*Define 2 hardcoded parameters*/
8#define OMEGA 0.001 // parameter controlling transition to nonlinear resistance in basal system (dimensionless)
9#define NU 1.787e-6 //kinematic water viscosity m^2/s
10#define CT 7.5e-8 // Clapeyron slope (K/Pa)
11#define CW 4.22e3 // specific heat capacity of water (J/kg/K)
12
13/*Model processing*/
14void HydrologySommersAnalysis::CreateConstraints(Constraints* constraints,IoModel* iomodel){/*{{{*/
15
16 /*retrieve some parameters: */
17 int hydrology_model;
18 iomodel->FindConstant(&hydrology_model,"md.hydrology.model");
19
20 if(hydrology_model!=HydrologysommersEnum) return;
21
22 IoModelToConstraintsx(constraints,iomodel,"md.hydrology.spchead",HydrologySommersAnalysisEnum,P1Enum);
23
24}/*}}}*/
25void HydrologySommersAnalysis::CreateLoads(Loads* loads, IoModel* iomodel){/*{{{*/
26
27 /*Fetch parameters: */
28 int hydrology_model;
29 iomodel->FindConstant(&hydrology_model,"md.hydrology.model");
30
31 /*Now, do we really want Sommers?*/
32 if(hydrology_model!=HydrologysommersEnum) return;
33
34 /*Create discrete loads for Moulins*/
35 CreateSingleNodeToElementConnectivity(iomodel);
36 for(int i=0;i<iomodel->numberofvertices;i++){
37 if (iomodel->domaintype!=Domain3DEnum){
38 /*keep only this partition's nodes:*/
39 if(iomodel->my_vertices[i]){
40 loads->AddObject(new Moulin(iomodel->loadcounter+i+1,i,iomodel,HydrologySommersAnalysisEnum));
41 }
42 }
43 else if(reCast<int>(iomodel->Data("md.mesh.vertexonbase")[i])){
44 if(iomodel->my_vertices[i]){
45 loads->AddObject(new Moulin(iomodel->loadcounter+i+1,i,iomodel,HydrologySommersAnalysisEnum));
46 }
47 }
48 }
49 iomodel->DeleteData(1,"md.mesh.vertexonbase");
50
51 /*Deal with Neumann BC*/
52 int M,N;
53 int *segments = NULL;
54 iomodel->FetchData(&segments,&M,&N,"md.mesh.segments");
55
56 /*Check that the size seem right*/
57 _assert_(N==3); _assert_(M>=3);
58 for(int i=0;i<M;i++){
59 if(iomodel->my_elements[segments[i*3+2]-1]){
60 loads->AddObject(new Neumannflux(iomodel->loadcounter+i+1,i,iomodel,segments,HydrologySommersAnalysisEnum));
61 }
62 }
63
64 xDelete<int>(segments);
65
66}/*}}}*/
67void HydrologySommersAnalysis::CreateNodes(Nodes* nodes,IoModel* iomodel){/*{{{*/
68
69 /*Fetch parameters: */
70 int hydrology_model;
71 iomodel->FindConstant(&hydrology_model,"md.hydrology.model");
72
73 /*Now, do we really want Sommers?*/
74 if(hydrology_model!=HydrologysommersEnum) return;
75
76 if(iomodel->domaintype==Domain3DEnum) iomodel->FetchData(2,"md.mesh.vertexonbase","md.mesh.vertexonsurface");
77 ::CreateNodes(nodes,iomodel,HydrologySommersAnalysisEnum,P1Enum);
78 iomodel->DeleteData(2,"md.mesh.vertexonbase","md.mesh.vertexonsurface");
79}/*}}}*/
80int HydrologySommersAnalysis::DofsPerNode(int** doflist,int domaintype,int approximation){/*{{{*/
81 return 1;
82}/*}}}*/
83void HydrologySommersAnalysis::UpdateElements(Elements* elements,IoModel* iomodel,int analysis_counter,int analysis_type){/*{{{*/
84
85 /*Fetch data needed: */
86 int hydrology_model,frictionlaw;
87 iomodel->FindConstant(&hydrology_model,"md.hydrology.model");
88
89 /*Now, do we really want Sommers?*/
90 if(hydrology_model!=HydrologysommersEnum) return;
91
92 /*Update elements: */
93 int counter=0;
94 for(int i=0;i<iomodel->numberofelements;i++){
95 if(iomodel->my_elements[i]){
96 Element* element=(Element*)elements->GetObjectByOffset(counter);
97 element->Update(i,iomodel,analysis_counter,analysis_type,P1Enum);
98 counter++;
99 }
100 }
101
102 iomodel->FetchDataToInput(elements,"md.geometry.thickness",ThicknessEnum);
103 iomodel->FetchDataToInput(elements,"md.geometry.base",BaseEnum);
104 if(iomodel->domaintype!=Domain2DhorizontalEnum){
105 iomodel->FetchDataToInput(elements,"md.mesh.vertexonbase",MeshVertexonbaseEnum);
106 iomodel->FetchDataToInput(elements,"md.mesh.vertexonsurface",MeshVertexonsurfaceEnum);
107 }
108 iomodel->FetchDataToInput(elements,"md.mask.ice_levelset",MaskIceLevelsetEnum);
109 iomodel->FetchDataToInput(elements,"md.mask.groundedice_levelset",MaskGroundediceLevelsetEnum);
110 iomodel->FetchDataToInput(elements,"md.basalforcings.groundedice_melting_rate",BasalforcingsGroundediceMeltingRateEnum);
111 iomodel->FetchDataToInput(elements,"md.basalforcings.geothermalflux",BasalforcingsGeothermalfluxEnum);
112 iomodel->FetchDataToInput(elements,"md.hydrology.head",HydrologyHeadEnum);
113 iomodel->FetchDataToInput(elements,"md.hydrology.gap_height",HydrologyGapHeightEnum);
114 iomodel->FetchDataToInput(elements,"md.hydrology.englacial_input",HydrologyEnglacialInputEnum);
115 iomodel->FetchDataToInput(elements,"md.hydrology.moulin_input",HydrologyMoulinInputEnum);
116 iomodel->FetchDataToInput(elements,"md.hydrology.bump_spacing",HydrologyBumpSpacingEnum);
117 iomodel->FetchDataToInput(elements,"md.hydrology.bump_height",HydrologyBumpHeightEnum);
118 iomodel->FetchDataToInput(elements,"md.hydrology.reynolds",HydrologyReynoldsEnum);
119 iomodel->FetchDataToInput(elements,"md.hydrology.neumannflux",HydrologyNeumannfluxEnum);
120 iomodel->FetchDataToInput(elements,"md.initialization.vx",VxEnum);
121 iomodel->FetchDataToInput(elements,"md.initialization.vy",VyEnum);
122 iomodel->FindConstant(&frictionlaw,"md.friction.law");
123
124 /*Friction law variables*/
125 switch(frictionlaw){
126 case 1:
127 iomodel->FetchDataToInput(elements,"md.friction.coefficient",FrictionCoefficientEnum);
128 iomodel->FetchDataToInput(elements,"md.friction.p",FrictionPEnum);
129 iomodel->FetchDataToInput(elements,"md.friction.q",FrictionQEnum);
130 break;
131 case 8:
132 iomodel->FetchDataToInput(elements,"md.friction.coefficient",FrictionCoefficientEnum);
133 break;
134 default:
135 _error_("Friction law "<< frictionlaw <<" not supported");
136 }
137}/*}}}*/
138void HydrologySommersAnalysis::UpdateParameters(Parameters* parameters,IoModel* iomodel,int solution_enum,int analysis_enum){/*{{{*/
139
140 /*retrieve some parameters: */
141 int hydrology_model;
142 int numoutputs;
143 char** requestedoutputs = NULL;
144 iomodel->FindConstant(&hydrology_model,"md.hydrology.model");
145
146 /*Now, do we really want Sommers?*/
147 if(hydrology_model!=HydrologysommersEnum) return;
148
149 parameters->AddObject(new IntParam(HydrologyModelEnum,hydrology_model));
150 parameters->AddObject(iomodel->CopyConstantObject("md.friction.law",FrictionLawEnum));
151 parameters->AddObject(iomodel->CopyConstantObject("md.hydrology.relaxation",HydrologyRelaxationEnum));
152 parameters->AddObject(iomodel->CopyConstantObject("md.hydrology.storage",HydrologyStorageEnum));
153
154 /*Requested outputs*/
155 iomodel->FindConstant(&requestedoutputs,&numoutputs,"md.hydrology.requested_outputs");
156 parameters->AddObject(new IntParam(HydrologyNumRequestedOutputsEnum,numoutputs));
157 if(numoutputs)parameters->AddObject(new StringArrayParam(HydrologyRequestedOutputsEnum,requestedoutputs,numoutputs));
158 iomodel->DeleteData(&requestedoutputs,numoutputs,"md.hydrology.requested_outputs");
159}/*}}}*/
160
161/*Finite Element Analysis*/
162void HydrologySommersAnalysis::Core(FemModel* femmodel){/*{{{*/
163 _error_("not implemented");
164}/*}}}*/
165ElementVector* HydrologySommersAnalysis::CreateDVector(Element* element){/*{{{*/
166 /*Default, return NULL*/
167 return NULL;
168}/*}}}*/
169ElementMatrix* HydrologySommersAnalysis::CreateJacobianMatrix(Element* element){/*{{{*/
170_error_("Not implemented");
171}/*}}}*/
172ElementMatrix* HydrologySommersAnalysis::CreateKMatrix(Element* element){/*{{{*/
173
174 /*Intermediaries */
175 IssmDouble Jdet;
176 IssmDouble* xyz_list = NULL;
177
178 /*Fetch number of nodes and dof for this finite element*/
179 int numnodes = element->GetNumberOfNodes();
180
181 /*Initialize Element vector and other vectors*/
182 ElementMatrix* Ke = element->NewElementMatrix();
183 IssmDouble* dbasis = xNew<IssmDouble>(2*numnodes);
184 IssmDouble* basis = xNew<IssmDouble>(numnodes);
185
186 /*Retrieve all inputs and parameters*/
187 element->GetVerticesCoordinates(&xyz_list);
188
189 /*Get conductivity from inputs*/
190 IssmDouble conductivity = GetConductivity(element);
191
192 /*Get englacial storage coefficient*/
193 IssmDouble storage,dt;
194 element->FindParam(&storage,HydrologyStorageEnum);
195 element->FindParam(&dt,TimesteppingTimeStepEnum);
196
197 /* Start looping on the number of gaussian points: */
198 Gauss* gauss=element->NewGauss(1);
199 for(int ig=gauss->begin();ig<gauss->end();ig++){
200 gauss->GaussPoint(ig);
201
202 element->JacobianDeterminant(&Jdet,xyz_list,gauss);
203 element->NodalFunctionsDerivatives(dbasis,xyz_list,gauss);
204 element->NodalFunctions(basis,gauss);
205
206 for(int i=0;i<numnodes;i++){
207 for(int j=0;j<numnodes;j++){
208 Ke->values[i*numnodes+j] += conductivity*gauss->weight*Jdet*(dbasis[0*numnodes+i]*dbasis[0*numnodes+j] + dbasis[1*numnodes+i]*dbasis[1*numnodes+j])
209 + gauss->weight*Jdet*storage/dt*basis[i]*basis[j];
210 }
211 }
212 }
213
214 /*Clean up and return*/
215 xDelete<IssmDouble>(xyz_list);
216 xDelete<IssmDouble>(basis);
217 xDelete<IssmDouble>(dbasis);
218 delete gauss;
219 return Ke;
220}/*}}}*/
221ElementVector* HydrologySommersAnalysis::CreatePVector(Element* element){/*{{{*/
222
223 /*Skip if water or ice shelf element*/
224 if(element->IsFloating()) return NULL;
225
226 /*Intermediaries */
227 IssmDouble Jdet,meltrate,G,dh[2],B,A,n;
228 IssmDouble gap,bed,thickness,head,ieb,head_old;
229 IssmDouble lr,br,vx,vy,beta;
230 IssmDouble alpha2,frictionheat;
231 IssmDouble PMPheat,dpressure_water[2],dbed[2];
232 IssmDouble* xyz_list = NULL;
233
234 /*Fetch number of nodes and dof for this finite element*/
235 int numnodes = element->GetNumberOfNodes();
236
237 /*Initialize Element vector and other vectors*/
238 ElementVector* pe = element->NewElementVector();
239 IssmDouble* basis = xNew<IssmDouble>(numnodes);
240
241 /*Retrieve all inputs and parameters*/
242 element->GetVerticesCoordinates(&xyz_list);
243 IssmDouble latentheat = element->GetMaterialParameter(MaterialsLatentheatEnum);
244 IssmDouble g = element->GetMaterialParameter(ConstantsGEnum);
245 IssmDouble rho_ice = element->GetMaterialParameter(MaterialsRhoIceEnum);
246 IssmDouble rho_water = element->GetMaterialParameter(MaterialsRhoFreshwaterEnum);
247 Input* geothermalflux_input = element->GetInput(BasalforcingsGeothermalfluxEnum);_assert_(geothermalflux_input);
248 Input* head_input = element->GetInput(HydrologyHeadEnum); _assert_(head_input);
249 Input* gap_input = element->GetInput(HydrologyGapHeightEnum); _assert_(gap_input);
250 Input* thickness_input = element->GetInput(ThicknessEnum); _assert_(thickness_input);
251 Input* base_input = element->GetInput(BaseEnum); _assert_(base_input);
252 Input* B_input = element->GetInput(MaterialsRheologyBEnum); _assert_(B_input);
253 Input* n_input = element->GetInput(MaterialsRheologyNEnum); _assert_(n_input);
254 Input* englacial_input = element->GetInput(HydrologyEnglacialInputEnum); _assert_(englacial_input);
255 Input* vx_input = element->GetInput(VxEnum); _assert_(vx_input);
256 Input* vy_input = element->GetInput(VyEnum); _assert_(vy_input);
257 Input* lr_input = element->GetInput(HydrologyBumpSpacingEnum); _assert_(lr_input);
258 Input* br_input = element->GetInput(HydrologyBumpHeightEnum); _assert_(br_input);
259 Input* headold_input = element->GetInput(HydrologyHeadOldEnum); _assert_(headold_input);
260
261 /*Get conductivity from inputs*/
262 IssmDouble conductivity = GetConductivity(element);
263
264 /*Get englacial storage coefficient*/
265 IssmDouble storage,dt;
266 element->FindParam(&storage,HydrologyStorageEnum);
267 element->FindParam(&dt,TimesteppingTimeStepEnum);
268
269 /*Build friction element, needed later: */
270 Friction* friction=new Friction(element,2);
271
272 /* Start looping on the number of gaussian points: */
273 Gauss* gauss=element->NewGauss(2);
274 for(int ig=gauss->begin();ig<gauss->end();ig++){
275 gauss->GaussPoint(ig);
276
277 element->JacobianDeterminant(&Jdet,xyz_list,gauss);
278 element->NodalFunctions(basis,gauss);
279 geothermalflux_input->GetInputValue(&G,gauss);
280 base_input->GetInputValue(&bed,gauss);
281 base_input->GetInputDerivativeValue(&dbed[0],xyz_list,gauss);
282 thickness_input->GetInputValue(&thickness,gauss);
283 gap_input->GetInputValue(&gap,gauss);
284 head_input->GetInputValue(&head,gauss);
285 head_input->GetInputDerivativeValue(&dh[0],xyz_list,gauss);
286 englacial_input->GetInputValue(&ieb,gauss);
287 lr_input->GetInputValue(&lr,gauss);
288 br_input->GetInputValue(&br,gauss);
289 vx_input->GetInputValue(&vx,gauss);
290 vy_input->GetInputValue(&vy,gauss);
291 headold_input->GetInputValue(&head_old,gauss);
292
293 /*Get ice A parameter*/
294 B_input->GetInputValue(&B,gauss);
295 n_input->GetInputValue(&n,gauss);
296 A=pow(B,-n);
297
298 /*Compute beta term*/
299 if(gap<br)
300 beta = (br-gap)/lr;
301 else
302 beta = 0.;
303
304 /*Compute frictional heat flux*/
305 friction->GetAlpha2(&alpha2,gauss);
306 vx_input->GetInputValue(&vx,gauss);
307 vy_input->GetInputValue(&vy,gauss);
308 frictionheat=alpha2*(vx*vx+vy*vy);
309
310 /*Get water and ice pressures*/
311 IssmDouble pressure_ice = rho_ice*g*thickness; _assert_(pressure_ice>0.);
312 IssmDouble pressure_water = rho_water*g*(head-bed);
313 if(pressure_water>pressure_ice) pressure_water = pressure_ice;
314
315 /*Get water pressure from previous time step to use in lagged creep term*/
316 IssmDouble pressure_water_old = rho_water*g*(head_old-bed);
317 if(pressure_water_old>pressure_ice) pressure_water_old = pressure_ice;
318
319 /*Compute change in sensible heat due to changes in pressure melting point*/
320 dpressure_water[0] = rho_water*g*(dh[0] - dbed[0]);
321 dpressure_water[1] = rho_water*g*(dh[1] - dbed[1]);
322 PMPheat=-CT*CW*conductivity*(dh[0]*dpressure_water[0]+dh[1]*dpressure_water[1]);
323
324 meltrate = 1/latentheat*(G+frictionheat+rho_water*g*conductivity*(dh[0]*dh[0]+dh[1]*dh[1])-PMPheat);
325 _assert_(meltrate>0.);
326
327 for(int i=0;i<numnodes;i++) pe->values[i]+=Jdet*gauss->weight*
328 (
329 meltrate*(1/rho_water-1/rho_ice)
330 +A*pow(fabs(pressure_ice - pressure_water),n-1)*(pressure_ice - pressure_water)*gap
331 -beta*sqrt(vx*vx+vy*vy)
332 +ieb
333 +storage*head_old/dt
334 )*basis[i];
335 }
336 /*Clean up and return*/
337 xDelete<IssmDouble>(xyz_list);
338 xDelete<IssmDouble>(basis);
339 delete friction;
340 delete gauss;
341 return pe;
342}/*}}}*/
343void HydrologySommersAnalysis::GetSolutionFromInputs(Vector<IssmDouble>* solution,Element* element){/*{{{*/
344 element->GetSolutionFromInputsOneDof(solution,HydrologyHeadEnum);
345}/*}}}*/
346void HydrologySommersAnalysis::GradientJ(Vector<IssmDouble>* gradient,Element* element,int control_type,int control_index){/*{{{*/
347 _error_("Not implemented yet");
348}/*}}}*/
349void HydrologySommersAnalysis::InputUpdateFromSolution(IssmDouble* solution,Element* element){/*{{{*/
350
351 /*Intermediary*/
352 IssmDouble dh[3];
353 int* doflist = NULL;
354 IssmDouble* xyz_list = NULL;
355
356 /*Get gravity from parameters*/
357 IssmDouble g = element->GetMaterialParameter(ConstantsGEnum);
358
359 /*Fetch number of nodes for this finite element*/
360 int numnodes = element->GetNumberOfNodes();
361
362 /*Fetch dof list and allocate solution vector*/
363 element->GetDofList(&doflist,NoneApproximationEnum,GsetEnum);
364 IssmDouble* values = xNew<IssmDouble>(numnodes);
365
366 /*Get thickness and base on nodes to apply cap on water head*/
367 IssmDouble* eff_pressure = xNew<IssmDouble>(numnodes);
368 IssmDouble* thickness = xNew<IssmDouble>(numnodes);
369 IssmDouble* bed = xNew<IssmDouble>(numnodes);
370 IssmDouble rho_ice = element->GetMaterialParameter(MaterialsRhoIceEnum);
371 IssmDouble rho_water = element->GetMaterialParameter(MaterialsRhoFreshwaterEnum);
372 element->GetInputListOnNodes(&thickness[0],ThicknessEnum);
373 element->GetInputListOnNodes(&bed[0],BaseEnum);
374
375 /*Get head from previous time-step and under-relaxation coefficient to use in under-relaxation for nonlinear convergence*/
376 IssmDouble* head_old = xNew<IssmDouble>(numnodes);
377 element->GetInputListOnNodes(&head_old[0],HydrologyHeadEnum);
378 IssmDouble relaxation;
379 element->FindParam(&relaxation,HydrologyRelaxationEnum);
380
381 /*Use the dof list to index into the solution vector: */
382 for(int i=0;i<numnodes;i++){
383 values[i]=solution[doflist[i]];
384
385 /*make sure that p_water<p_ice -> h<rho_i H/rho_w + zb*/
386 if(values[i]>rho_ice*thickness[i]/rho_water+bed[i]){
387 values[i] = rho_ice*thickness[i]/rho_water+bed[i];
388 }
389
390 /*Make sure that negative pressure is not allowed*/
391 // if(values[i]<bed[i]){
392 // values[i] = bed[i];
393 // }
394
395 /*Under-relaxation*/
396 values[i] = head_old[i] - relaxation*(head_old[i]-values[i]);
397
398 /*Calculate effective pressure*/
399 eff_pressure[i] = rho_ice*g*thickness[i] - rho_water*g*(values[i]-bed[i]);
400
401 if(xIsNan<IssmDouble>(values[i])) _error_("NaN found in solution vector");
402 if(xIsInf<IssmDouble>(values[i])) _error_("Inf found in solution vector");
403 }
404
405 /*Add input to the element: */
406 element->AddInput(HydrologyHeadEnum,values,element->GetElementType());
407 element->AddInput(EffectivePressureEnum,eff_pressure,P1Enum);
408
409 /*Update reynolds number according to new solution*/
410 element->GetVerticesCoordinates(&xyz_list);
411 Input* head_input = element->GetInput(HydrologyHeadEnum);_assert_(head_input);
412 head_input->GetInputDerivativeAverageValue(&dh[0],xyz_list);
413 IssmDouble conductivity = GetConductivity(element);
414
415 IssmDouble reynolds = conductivity*sqrt(dh[0]*dh[0]+dh[1]*dh[1])/NU;
416 element->AddInput(HydrologyReynoldsEnum,&reynolds,P0Enum);
417
418 /*Free resources:*/
419 xDelete<IssmDouble>(values);
420 xDelete<IssmDouble>(thickness);
421 xDelete<IssmDouble>(bed);
422 xDelete<IssmDouble>(xyz_list);
423 xDelete<int>(doflist);
424 xDelete<IssmDouble>(eff_pressure);
425 xDelete<IssmDouble>(head_old);
426}/*}}}*/
427void HydrologySommersAnalysis::UpdateConstraints(FemModel* femmodel){/*{{{*/
428 /*Default, do nothing*/
429 return;
430}/*}}}*/
431
432/*Additional methods*/
433IssmDouble HydrologySommersAnalysis::GetConductivity(Element* element){/*{{{*/
434
435 /*Intermediaries */
436 IssmDouble gap,reynolds;
437
438 /*Get gravity from parameters*/
439 IssmDouble g = element->GetMaterialParameter(ConstantsGEnum);
440
441 /*Get Reynolds and gap average values*/
442 Input* reynolds_input = element->GetInput(HydrologyReynoldsEnum); _assert_(reynolds_input);
443 Input* gap_input = element->GetInput(HydrologyGapHeightEnum); _assert_(gap_input);
444 reynolds_input->GetInputAverage(&reynolds);
445 gap_input->GetInputAverage(&gap);
446
447 /*Compute conductivity*/
448 IssmDouble conductivity = pow(gap,3)*g/(12.*NU*(1+OMEGA*reynolds));
449 _assert_(conductivity>0);
450
451 /*Clean up and return*/
452 return conductivity;
453}/*}}}*/
454void HydrologySommersAnalysis::UpdateGapHeight(FemModel* femmodel){/*{{{*/
455
456
457 for(int j=0;j<femmodel->elements->Size();j++){
458 Element* element=(Element*)femmodel->elements->GetObjectByOffset(j);
459 UpdateGapHeight(element);
460 }
461
462}/*}}}*/
463void HydrologySommersAnalysis::UpdateGapHeight(Element* element){/*{{{*/
464
465 /*Skip if water or ice shelf element*/
466 if(element->IsFloating()) return;
467
468 /*Intermediaries */
469 IssmDouble newgap = 0.;
470 IssmDouble Jdet,meltrate,G,dh[2],B,A,n,dt;
471 IssmDouble gap,bed,thickness,head,ieb;
472 IssmDouble lr,br,vx,vy,beta;
473 IssmDouble alpha2,frictionheat;
474 IssmDouble* xyz_list = NULL;
475 IssmDouble dpressure_water[2],dbed[2],PMPheat;
476 IssmDouble q = 0.;
477 IssmDouble channelization = 0.;
478
479 /*Retrieve all inputs and parameters*/
480 element->GetVerticesCoordinates(&xyz_list);
481 element->FindParam(&dt,TimesteppingTimeStepEnum);
482 IssmDouble latentheat = element->GetMaterialParameter(MaterialsLatentheatEnum);
483 IssmDouble g = element->GetMaterialParameter(ConstantsGEnum);
484 IssmDouble rho_ice = element->GetMaterialParameter(MaterialsRhoIceEnum);
485 IssmDouble rho_water = element->GetMaterialParameter(MaterialsRhoFreshwaterEnum);
486 Input* geothermalflux_input = element->GetInput(BasalforcingsGeothermalfluxEnum);_assert_(geothermalflux_input);
487 Input* head_input = element->GetInput(HydrologyHeadEnum); _assert_(head_input);
488 Input* gap_input = element->GetInput(HydrologyGapHeightEnum); _assert_(gap_input);
489 Input* thickness_input = element->GetInput(ThicknessEnum); _assert_(thickness_input);
490 Input* base_input = element->GetInput(BaseEnum); _assert_(base_input);
491 Input* B_input = element->GetInput(MaterialsRheologyBEnum); _assert_(B_input);
492 Input* n_input = element->GetInput(MaterialsRheologyNEnum); _assert_(n_input);
493 Input* englacial_input = element->GetInput(HydrologyEnglacialInputEnum); _assert_(englacial_input);
494 Input* vx_input = element->GetInput(VxEnum); _assert_(vx_input);
495 Input* vy_input = element->GetInput(VyEnum); _assert_(vy_input);
496 Input* lr_input = element->GetInput(HydrologyBumpSpacingEnum); _assert_(lr_input);
497 Input* br_input = element->GetInput(HydrologyBumpHeightEnum); _assert_(br_input);
498
499 /*Get conductivity from inputs*/
500 IssmDouble conductivity = GetConductivity(element);
501
502 /*Build friction element, needed later: */
503 Friction* friction=new Friction(element,2);
504
505 /*Keep track of weights*/
506 IssmDouble totalweights=0.;
507
508 /* Start looping on the number of gaussian points: */
509 Gauss* gauss=element->NewGauss(2);
510 for(int ig=gauss->begin();ig<gauss->end();ig++){
511 gauss->GaussPoint(ig);
512
513 element->JacobianDeterminant(&Jdet,xyz_list,gauss);
514
515 geothermalflux_input->GetInputValue(&G,gauss);
516 base_input->GetInputValue(&bed,gauss);
517 base_input->GetInputDerivativeValue(&dbed[0],xyz_list,gauss);
518 thickness_input->GetInputValue(&thickness,gauss);
519 gap_input->GetInputValue(&gap,gauss);
520 head_input->GetInputValue(&head,gauss);
521 head_input->GetInputDerivativeValue(&dh[0],xyz_list,gauss);
522 englacial_input->GetInputValue(&ieb,gauss);
523 lr_input->GetInputValue(&lr,gauss);
524 br_input->GetInputValue(&br,gauss);
525 vx_input->GetInputValue(&vx,gauss);
526 vy_input->GetInputValue(&vy,gauss);
527
528 /*Get ice A parameter*/
529 B_input->GetInputValue(&B,gauss);
530 n_input->GetInputValue(&n,gauss);
531 A=pow(B,-n);
532
533 /*Compute beta term*/
534 if(gap<br)
535 beta = (br-gap)/lr;
536 else
537 beta = 0.;
538
539 /*Compute frictional heat flux*/
540 friction->GetAlpha2(&alpha2,gauss);
541 vx_input->GetInputValue(&vx,gauss);
542 vy_input->GetInputValue(&vy,gauss);
543 frictionheat=alpha2*(vx*vx+vy*vy);
544
545 /*Get water and ice pressures*/
546 IssmDouble pressure_ice = rho_ice*g*thickness; _assert_(pressure_ice>0.);
547 IssmDouble pressure_water = rho_water*g*(head-bed);
548 if(pressure_water>pressure_ice) pressure_water = pressure_ice;
549
550 /* Compute change in sensible heat due to changes in pressure melting point*/
551 dpressure_water[0] = rho_water*g*(dh[0] - dbed[0]);
552 dpressure_water[1] = rho_water*g*(dh[1] - dbed[1]);
553 PMPheat=-CT*CW*conductivity*(dh[0]*dpressure_water[0]+dh[1]*dpressure_water[1]);
554
555 meltrate = 1/latentheat*(G+frictionheat+rho_water*g*conductivity*(dh[0]*dh[0]+dh[1]*dh[1])-PMPheat);
556 _assert_(meltrate>0.);
557
558 newgap += gauss->weight*Jdet*(gap+dt*(
559 meltrate/rho_ice
560 -A*pow(fabs(pressure_ice-pressure_water),n-1)*(pressure_ice-pressure_water)*gap
561 +beta*sqrt(vx*vx+vy*vy)
562 ));
563 totalweights +=gauss->weight*Jdet;
564
565 /* Compute basal water flux */
566 q += gauss->weight*Jdet*(conductivity*sqrt(dh[0]*dh[0]+dh[1]*dh[1]));
567
568 /* Compute "degree of channelization" (ratio of melt opening to opening by sliding) */
569 channelization += gauss->weight*Jdet*(meltrate/rho_ice/(meltrate/rho_ice+beta*sqrt(vx*vx+vy*vy)));
570 }
571
572 /*Divide by connectivity*/
573 newgap = newgap/totalweights;
574 IssmDouble mingap = 1e-6;
575 if(newgap<mingap) newgap=mingap;
576
577 /*Limit gap height to grow to surface*/
578 if(newgap>thickness)
579 newgap = thickness;
580
581 /*Add new gap as an input*/
582 element->AddInput(HydrologyGapHeightEnum,&newgap,P0Enum);
583
584 /*Divide by connectivity, add basal flux as an input*/
585 q = q/totalweights;
586 element->AddInput(HydrologyBasalFluxEnum,&q,P0Enum);
587
588 /* Divide by connectivity, add degree of channelization as an input */
589 channelization = channelization/totalweights;
590 element->AddInput(DegreeOfChannelizationEnum,&channelization,P0Enum);
591
592 /*Clean up and return*/
593 xDelete<IssmDouble>(xyz_list);
594 delete friction;
595 delete gauss;
596}/*}}}*/
Note: See TracBrowser for help on using the repository browser.