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

Last change on this file since 20500 was 20500, checked in by Mathieu Morlighem, 9 years ago

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