source: issm/oecreview/Archive/21337-21723/ISSM-21429-21430.diff

Last change on this file was 21726, checked in by Mathieu Morlighem, 8 years ago

CHG added Archive/21337-21723
File size: 23.2 KB
RevLine 
[21726]1Index: ../trunk-jpl/src/c/analyses/HydrologyDCInefficientAnalysis.cpp
2===================================================================
3--- ../trunk-jpl/src/c/analyses/HydrologyDCInefficientAnalysis.cpp (revision 21429)
4+++ ../trunk-jpl/src/c/analyses/HydrologyDCInefficientAnalysis.cpp (revision 21430)
5@@ -18,10 +18,10 @@
6 int penalty_lock;
7 int hydro_maxiter;
8 bool isefficientlayer;
9- IssmDouble sedimentlimit;
10 IssmDouble penalty_factor;
11- IssmDouble leakagefactor;
12 IssmDouble rel_tol;
13+ IssmDouble leakagefactor;
14+ IssmDouble sedimentlimit;
15
16 /*retrieve some parameters: */
17 iomodel->FindConstant(&hydrology_model,"md.hydrology.model");
18@@ -29,32 +29,30 @@
19 /*Now, do we really want DC?*/
20 if(hydrology_model!=HydrologydcEnum) return;
21
22- iomodel->FetchData(&isefficientlayer, "md.hydrology.isefficientlayer");
23 iomodel->FetchData(&sedimentlimit_flag, "md.hydrology.sedimentlimit_flag" );
24 iomodel->FetchData(&transfer_flag, "md.hydrology.transfer_flag" );
25 iomodel->FetchData(&penalty_factor, "md.hydrology.penalty_factor" );
26- iomodel->FetchData(&rel_tol, "md.hydrology.rel_tol" );
27- iomodel->FetchData(&penalty_lock, "md.hydrology.penalty_lock" );
28+ iomodel->FetchData(&isefficientlayer, "md.hydrology.isefficientlayer");
29 iomodel->FetchData(&hydro_maxiter, "md.hydrology.max_iter" );
30+ iomodel->FetchData(&penalty_lock, "md.hydrology.penalty_lock" );
31+ iomodel->FetchData(&rel_tol, "md.hydrology.rel_tol" );
32
33- if(sedimentlimit_flag==1){
34- iomodel->FetchData(&sedimentlimit,"md.hydrology.sedimentlimit");
35- parameters->AddObject(new DoubleParam(HydrologydcSedimentlimitEnum,sedimentlimit));
36- }
37-
38- if(transfer_flag==1){
39- iomodel->FetchData(&leakagefactor,"md.hydrology.leakage_factor");
40- parameters->AddObject(new DoubleParam(HydrologydcLeakageFactorEnum,leakagefactor));
41- }
42-
43- parameters->AddObject(new DoubleParam(HydrologydcPenaltyFactorEnum,penalty_factor));
44 parameters->AddObject(new IntParam(HydrologyModelEnum,hydrology_model));
45- parameters->AddObject(new BoolParam(HydrologydcIsefficientlayerEnum,isefficientlayer));
46 parameters->AddObject(new IntParam(HydrologydcSedimentlimitFlagEnum,sedimentlimit_flag));
47 parameters->AddObject(new IntParam(HydrologydcTransferFlagEnum,transfer_flag));
48- parameters->AddObject(new DoubleParam(HydrologydcRelTolEnum,rel_tol));
49 parameters->AddObject(new IntParam(HydrologydcPenaltyLockEnum,penalty_lock));
50 parameters->AddObject(new IntParam(HydrologydcMaxIterEnum,hydro_maxiter));
51+ parameters->AddObject(new BoolParam(HydrologydcIsefficientlayerEnum,isefficientlayer));
52+ parameters->AddObject(new DoubleParam(HydrologydcPenaltyFactorEnum,penalty_factor));
53+ parameters->AddObject(new DoubleParam(HydrologydcRelTolEnum,rel_tol));
54+ if(transfer_flag==1){
55+ iomodel->FetchData(&leakagefactor,"md.hydrology.leakage_factor");
56+ parameters->AddObject(new DoubleParam(HydrologydcLeakageFactorEnum,leakagefactor));
57+ }
58+ if(sedimentlimit_flag==1){
59+ iomodel->FetchData(&sedimentlimit,"md.hydrology.sedimentlimit");
60+ parameters->AddObject(new DoubleParam(HydrologydcSedimentlimitEnum,sedimentlimit));
61+ }
62 }/*}}}*/
63
64 void HydrologyDCInefficientAnalysis::UpdateElements(Elements* elements,IoModel* iomodel,int analysis_counter,int analysis_type){/*{{{*/
65@@ -108,7 +106,9 @@
66 /*Now, do we really want DC?*/
67 if(hydrology_model!=HydrologydcEnum) return;
68
69- if(iomodel->domaintype!=Domain2DhorizontalEnum) iomodel->FetchData(2,"md.mesh.vertexonbase","md.mesh.vertexonsurface");
70+ if(iomodel->domaintype!=Domain2DhorizontalEnum){
71+ iomodel->FetchData(2,"md.mesh.vertexonbase","md.mesh.vertexonsurface");
72+ }
73 ::CreateNodes(nodes,iomodel,HydrologyDCInefficientAnalysisEnum,P1Enum);
74 iomodel->DeleteData(2,"md.mesh.vertexonbase","md.mesh.vertexonsurface");
75 }/*}}}*/
76@@ -130,8 +130,9 @@
77 iomodel->FindConstant(&hydrology_model,"md.hydrology.model");
78 if(hydrology_model!=HydrologydcEnum) return;
79
80- if(iomodel->domaintype==Domain3DEnum) iomodel->FetchData(1,"md.mesh.vertexonbase");
81-
82+ if(iomodel->domaintype==Domain3DEnum){
83+ iomodel->FetchData(1,"md.mesh.vertexonbase");
84+ }
85 //create penalties for nodes: no node can have water above the max
86 CreateSingleNodeToElementConnectivity(iomodel);
87 for(int i=0;i<iomodel->numberofvertices;i++){
88@@ -189,7 +190,7 @@
89 bool active_element,isefficientlayer;
90 IssmDouble D_scalar,Jdet,dt;
91 IssmDouble sediment_transmitivity;
92- IssmDouble prestep_head,base_elev;
93+ IssmDouble base_elev;
94 IssmDouble transfer,sediment_storing;
95 IssmDouble *xyz_list = NULL;
96
97@@ -300,7 +301,6 @@
98 IssmDouble Jdet;
99
100 IssmDouble *xyz_list = NULL;
101- Input* old_wh_input = NULL;
102 Input* active_element_input = NULL;
103
104 /*Fetch number of nodes and dof for this finite element*/
105@@ -320,19 +320,16 @@
106 Input* thick_input = basalelement->GetInput(ThicknessEnum);
107 Input* base_input = basalelement->GetInput(BaseEnum);
108 Input* water_input = basalelement->GetInput(BasalforcingsGroundediceMeltingRateEnum); _assert_(water_input);
109- if(dt!= 0.){old_wh_input = basalelement->GetInput(SedimentHeadOldEnum); _assert_(old_wh_input);}
110
111 /*Transfer related Inputs*/
112 if(isefficientlayer){
113 active_element_input = basalelement->GetInput(HydrologydcMaskEplactiveEltEnum); _assert_(active_element_input);
114 }
115
116-
117 /* Start looping on the number of gaussian points: */
118 Gauss* gauss=basalelement->NewGauss(2);
119 for(int ig=gauss->begin();ig<gauss->end();ig++){
120 gauss->GaussPoint(ig);
121-
122 basalelement->JacobianDeterminant(&Jdet,xyz_list,gauss);
123 basalelement->NodalFunctions(basis,gauss);
124
125@@ -358,11 +355,9 @@
126 }
127 }
128
129-
130 /*Transient and transfer terms*/
131 if(dt!=0.){
132- old_wh_input ->GetInputValue(&water_head,gauss);
133- sediment_storing = SedimentStoring(basalelement,gauss,sed_head_input,base_input);//sed_head_input
134+ sediment_storing = SedimentStoring(basalelement,gauss,sed_head_input,base_input);
135 if(isefficientlayer){
136 /*Dealing with the sediment part of the transfer term*/
137 active_element_input->GetInputValue(&active_element);
138@@ -427,20 +422,26 @@
139
140 void HydrologyDCInefficientAnalysis::InputUpdateFromSolution(IssmDouble* solution,Element* element){/*{{{*/
141
142- int domaintype;
143- bool converged;
144- int* doflist=NULL;
145- Element* basalelement=NULL;
146+ /*Intermediaries*/
147+ int domaintype;
148+ Element* basalelement;
149+ bool converged;
150+ int* doflist = NULL;
151
152+ /*Get basal element*/
153 element->FindParam(&domaintype,DomainTypeEnum);
154- if(domaintype!=Domain2DhorizontalEnum){
155- if(!element->IsOnBase()) return;
156- basalelement=element->SpawnBasalElement();
157+ switch(domaintype){
158+ case Domain2DhorizontalEnum:
159+ basalelement = element;
160+ break;
161+ case Domain3DEnum:
162+ if(!element->IsOnBase()) return NULL;
163+ basalelement = element->SpawnBasalElement();
164+ break;
165+ default: _error_("mesh "<<EnumToStringx(domaintype)<<" not supported yet");
166 }
167- else{
168- basalelement = element;
169- }
170
171+
172 /*Fetch number of nodes for this finite element*/
173 int numnodes = basalelement->GetNumberOfNodes();
174
175@@ -478,7 +479,6 @@
176 kappa=kmax*pow(10.,penalty_factor);
177
178 for(int i=0;i<numnodes;i++){
179-
180 GetHydrologyDCInefficientHmax(&h_max,basalelement,basalelement->GetNode(i));
181 if(values[i]>h_max) {
182 residual[i] = kappa*(values[i]-h_max);
183@@ -526,18 +526,22 @@
184 sed_head_input->GetInputValue(&prestep_head,gauss);
185 water_sheet=max(0.0,(prestep_head-(base_elev-sediment_thickness)));
186 specific_porosity=sediment_porosity-rho_freshwater*g*sediment_porosity*sediment_thickness*(water_compressibility+(sediment_compressibility/sediment_porosity));
187- if (water_sheet<=0.9*sediment_thickness){//porosity for unconfined region
188+ //porosity for unconfined region
189+ if (water_sheet<=0.9*sediment_thickness){
190 sediment_storing=sediment_porosity;
191 }
192- else if((water_sheet<sediment_thickness) && (water_sheet>0.9*sediment_thickness)){//continuity ramp
193+ //continuity ramp
194+ else if((water_sheet<sediment_thickness) && (water_sheet>0.9*sediment_thickness)){
195 sediment_storing=(sediment_thickness-water_sheet)*specific_porosity/(0.1*sediment_thickness)+
196 rho_freshwater*g*sediment_porosity*sediment_thickness*(water_compressibility+(sediment_compressibility/sediment_porosity));
197 }
198- else{//storing coefficient for confined
199+ //storing coefficient for confined
200+ else{
201 sediment_storing=rho_freshwater*g*sediment_porosity*sediment_thickness*(water_compressibility+(sediment_compressibility/sediment_porosity));
202 }
203 return sediment_storing;
204 }/*}}}*/
205+
206 IssmDouble HydrologyDCInefficientAnalysis::SedimentTransmitivity(Element* element,Gauss* gauss,Input* sed_head_input, Input* base_input,Input* SedTrans_input){/*{{{*/
207 IssmDouble sediment_transmitivity;
208 IssmDouble FullLayer_transmitivity;
209@@ -573,13 +577,10 @@
210 element->FindParam(&h_max,HydrologydcSedimentlimitEnum);
211 break;
212 case 2:
213-
214 rho_water = element->GetMaterialParameter(MaterialsRhoFreshwaterEnum);
215 rho_ice = element->GetMaterialParameter(MaterialsRhoIceEnum);
216-
217 _assert_(thick_input);
218 _assert_(base_input);
219-
220 /*Compute max*/
221 thick_input->GetInputValue(&thickness,gauss);
222 base_input->GetInputValue(&bed,gauss);
223@@ -633,11 +634,7 @@
224 IssmDouble HydrologyDCInefficientAnalysis::GetHydrologyKMatrixTransfer(Element* element, Gauss* gauss, Input* sed_head_input, Input* epl_head_input, Input* thick_input, Input* base_input){/*{{{*/
225
226 int transfermethod;
227- IssmDouble hmax;
228- IssmDouble epl_head,sediment_head;
229 IssmDouble leakage,transfer;
230- IssmDouble continuum, factor;
231-
232 element->FindParam(&transfermethod,HydrologydcTransferFlagEnum);
233
234 /*Switch between the different transfer methods cases*/
235@@ -647,23 +644,8 @@
236 transfer=0.0;
237 break;
238 case 1:
239- _assert_(sed_head_input);
240- _assert_(epl_head_input);
241-
242- sed_head_input->GetInputValue(&sediment_head,gauss);
243- epl_head_input->GetInputValue(&epl_head,gauss);
244 element->FindParam(&leakage,HydrologydcLeakageFactorEnum);
245-
246- hmax=GetHydrologyDCInefficientHmax(element, gauss, thick_input, base_input);
247-
248- //Computing continuum function to apply to transfer term, transfer is null only if
249- //epl_head>sediment_head AND sediment_head>h_max
250- //let's try without that for a while
251- /* continuum=((1.0/(1.0+exp(-20.0*(sediment_head-epl_head)))))+(1.0/(1.0+exp(-20.0*(hmax-sediment_head)))); */
252- /* factor=min(continuum,1.0); */
253- /* transfer=leakage*factor; */
254 transfer=leakage;
255-
256 break;
257 default:
258 _error_("no case higher than 1 for the Transfer method");
259@@ -674,11 +656,8 @@
260 IssmDouble HydrologyDCInefficientAnalysis::GetHydrologyPVectorTransfer(Element* element, Gauss* gauss, Input* sed_head_input, Input* epl_head_input, Input* thick_input, Input* base_input){/*{{{*/
261
262 int transfermethod;
263- IssmDouble hmax;
264- IssmDouble epl_head,sediment_head;
265+ IssmDouble epl_head;
266 IssmDouble leakage,transfer;
267- IssmDouble continuum, factor;
268-
269 element->FindParam(&transfermethod,HydrologydcTransferFlagEnum);
270
271 /*Switch between the different transfer methods cases*/
272@@ -688,24 +667,10 @@
273 transfer=0.0;
274 break;
275 case 1:
276- _assert_(sed_head_input);
277 _assert_(epl_head_input);
278-
279- sed_head_input->GetInputValue(&sediment_head,gauss);
280 epl_head_input->GetInputValue(&epl_head,gauss);
281 element->FindParam(&leakage,HydrologydcLeakageFactorEnum);
282-
283- hmax=GetHydrologyDCInefficientHmax(element, gauss, thick_input, base_input);
284-
285- //Computing continuum function to apply to transfer term, transfer is null only if
286- //epl_head>sediment_head AND sediment_head>h_max
287- //let's try without that for a while
288- /* continuum=((1.0/(1.0+exp(-20.0*(sediment_head-epl_head)))))+(1.0/(1.0+exp(-20.0*(hmax-sediment_head)))); */
289- /* factor=min(continuum,1.0); */
290- /* transfer=epl_head*leakage*factor; */
291-
292 transfer=epl_head*leakage;
293-
294 break;
295 default:
296 _error_("no case higher than 1 for the Transfer method");
297@@ -730,4 +695,5 @@
298 }
299 element->AddInput(new BoolInput(HydrologydcMaskEplactiveEltEnum,element_active));
300 }
301+ delete element;
302 }/*}}}*/
303Index: ../trunk-jpl/src/c/analyses/HydrologyDCEfficientAnalysis.cpp
304===================================================================
305--- ../trunk-jpl/src/c/analyses/HydrologyDCEfficientAnalysis.cpp (revision 21429)
306+++ ../trunk-jpl/src/c/analyses/HydrologyDCEfficientAnalysis.cpp (revision 21430)
307@@ -8,6 +8,7 @@
308 int HydrologyDCEfficientAnalysis::DofsPerNode(int** doflist,int domaintype,int approximation){/*{{{*/
309 return 1;
310 }/*}}}*/
311+
312 void HydrologyDCEfficientAnalysis::UpdateParameters(Parameters* parameters,IoModel* iomodel,int solution_enum,int analysis_enum){/*{{{*/
313
314 int hydrology_model;
315@@ -32,9 +33,8 @@
316
317 iomodel->FetchData(&eplthickcomp,"md.hydrology.epl_thick_comp");
318 parameters->AddObject(new IntParam(HydrologydcEplThickCompEnum,eplthickcomp));
319-
320-
321 }/*}}}*/
322+
323 void HydrologyDCEfficientAnalysis::UpdateElements(Elements* elements,IoModel* iomodel,int analysis_counter,int analysis_type){/*{{{*/
324
325 bool isefficientlayer;
326@@ -82,7 +82,9 @@
327 iomodel->FindConstant(&isefficientlayer,"md.hydrology.isefficientlayer");
328 if(!isefficientlayer) return;
329
330- if(iomodel->domaintype!=Domain2DhorizontalEnum) iomodel->FetchData(2,"md.mesh.vertexonbase","md.mesh.vertexonsurface");
331+ if(iomodel->domaintype!=Domain2DhorizontalEnum){
332+ iomodel->FetchData(2,"md.mesh.vertexonbase","md.mesh.vertexonsurface");
333+ }
334 ::CreateNodes(nodes,iomodel,HydrologyDCEfficientAnalysisEnum,P1Enum);
335 iomodel->DeleteData(2,"md.mesh.vertexonbase","md.mesh.vertexonsurface");
336 }/*}}}*/
337@@ -105,14 +107,22 @@
338 void HydrologyDCEfficientAnalysis::CreateLoads(Loads* loads, IoModel* iomodel){/*{{{*/
339 /*Nothing for now*/
340
341- /*Fetch parameters: */
342+ /*Do we really want DC?*/
343 int hydrology_model;
344 iomodel->FindConstant(&hydrology_model,"md.hydrology.model");
345 if(hydrology_model!=HydrologydcEnum) return;
346
347- if(iomodel->domaintype==Domain3DEnum) iomodel->FetchData(1,"md.mesh.vertexonbase");
348+ /*Do we want an efficient layer*/
349+ bool isefficientlayer;
350+ iomodel->FindConstant(&isefficientlayer,"md.hydrology.isefficientlayer");
351+ if(!isefficientlayer) return;
352
353- //create penalties for nodes: no node can have water above the max
354+ /*Fetch parameters: */
355+ if(iomodel->domaintype==Domain3DEnum){
356+ iomodel->FetchData(1,"md.mesh.vertexonbase");
357+ }
358+
359+ //Add moulin inputs as loads
360 CreateSingleNodeToElementConnectivity(iomodel);
361 for(int i=0;i<iomodel->numberofvertices;i++){
362 if (iomodel->domaintype!=Domain3DEnum){
363@@ -131,8 +141,8 @@
364 }/*}}}*/
365
366 void HydrologyDCEfficientAnalysis::InitZigZagCounter(FemModel* femmodel){/*{{{*/
367+
368 int* eplzigzag_counter =NULL;
369-
370 eplzigzag_counter=xNewZeroInit<int>(femmodel->nodes->Size());
371 femmodel->parameters->AddObject(new IntVecParam(EplZigZagCounterEnum,eplzigzag_counter,femmodel->nodes->Size()));
372 xDelete<int>(eplzigzag_counter);
373@@ -141,11 +151,8 @@
374 void HydrologyDCEfficientAnalysis::ResetCounter(FemModel* femmodel){/*{{{*/
375
376 int* eplzigzag_counter=NULL;
377- Element* element=NULL;
378-
379 femmodel->parameters->FindParam(&eplzigzag_counter,NULL,EplZigZagCounterEnum);
380 for(int i=0;i<femmodel->nodes->Size();i++){
381-
382 eplzigzag_counter[i]=0;
383 }
384 femmodel->parameters->SetParam(eplzigzag_counter,femmodel->nodes->Size(),EplZigZagCounterEnum);
385@@ -298,7 +305,7 @@
386
387 /*Check that all nodes are active, else return empty matrix*/
388 if(!active_element) {
389- if(domaintype!=Domain2DhorizontalEnum){
390+ if(domaintype!=Domain2DhorizontalEnum){
391 basalelement->DeleteMaterials();
392 delete basalelement;
393 }
394@@ -341,7 +348,6 @@
395 Gauss* gauss=basalelement->NewGauss(2);
396 for(int ig=gauss->begin();ig<gauss->end();ig++){
397 gauss->GaussPoint(ig);
398-
399 basalelement ->JacobianDeterminant(&Jdet,xyz_list,gauss);
400 basalelement ->NodalFunctions(basis,gauss);
401
402@@ -395,20 +401,24 @@
403
404 void HydrologyDCEfficientAnalysis::InputUpdateFromSolution(IssmDouble* solution,Element* element){/*{{{*/
405
406- bool active_element;
407- int domaintype,i;
408- Element* basalelement=NULL;
409+ /*Intermediaries*/
410+ bool active_element;
411+ int domaintype;
412+ Element* basalelement;
413
414+ /*Get basal element*/
415 element->FindParam(&domaintype,DomainTypeEnum);
416-
417- if(domaintype!=Domain2DhorizontalEnum){
418- if(!element->IsOnBase()) return;
419- basalelement=element->SpawnBasalElement();
420+ switch(domaintype){
421+ case Domain2DhorizontalEnum:
422+ basalelement = element;
423+ break;
424+ case Domain3DEnum:
425+ if(!element->IsOnBase()) return NULL;
426+ basalelement = element->SpawnBasalElement();
427+ break;
428+ default: _error_("mesh "<<EnumToStringx(domaintype)<<" not supported yet");
429 }
430- else{
431- basalelement = element;
432- }
433-
434+
435 /*Intermediary*/
436 int* doflist = NULL;
437
438@@ -469,11 +479,7 @@
439 IssmDouble HydrologyDCEfficientAnalysis::GetHydrologyKMatrixTransfer(Element* element, Gauss* gauss, Input* sed_head_input, Input* epl_head_input, Input* thick_input, Input* base_input){/*{{{*/
440
441 int transfermethod;
442- IssmDouble hmax;
443- IssmDouble epl_head,sediment_head;
444 IssmDouble leakage,transfer;
445- IssmDouble continuum, factor;
446- HydrologyDCInefficientAnalysis* inefanalysis = NULL;
447
448 element->FindParam(&transfermethod,HydrologydcTransferFlagEnum);
449
450@@ -484,23 +490,7 @@
451 transfer=0.0;
452 break;
453 case 1:
454- _assert_(sed_head_input);
455- _assert_(epl_head_input);
456-
457- inefanalysis = new HydrologyDCInefficientAnalysis();
458- hmax = inefanalysis->GetHydrologyDCInefficientHmax(element,gauss,thick_input,base_input);
459- delete inefanalysis;
460-
461- sed_head_input->GetInputValue(&sediment_head,gauss);
462- epl_head_input->GetInputValue(&epl_head,gauss);
463 element->FindParam(&leakage,HydrologydcLeakageFactorEnum);
464-
465- //Computing continuum function to apply to transfer term, transfer is null only if
466- // epl_head>sediment_head AND sediment_head>h_max
467- //let's try without that for a while
468- /* continuum=((1.0/(1.0+exp(-20.0*(sediment_head-epl_head)))))+(1.0/(1.0+exp(-20.0*(hmax-sediment_head)))); */
469- /* factor=min(continuum,1.0); */
470- /* transfer=leakage*factor; */
471 transfer=leakage;
472 break;
473 default:
474@@ -512,13 +502,9 @@
475 IssmDouble HydrologyDCEfficientAnalysis::GetHydrologyPVectorTransfer(Element* element, Gauss* gauss, Input* sed_head_input, Input* epl_head_input, Input* thick_input, Input* base_input){/*{{{*/
476
477 int transfermethod;
478- IssmDouble hmax;
479- IssmDouble epl_head,sediment_head;
480+ IssmDouble sediment_head;
481 IssmDouble leakage,transfer;
482- IssmDouble continuum, factor;
483
484- HydrologyDCInefficientAnalysis* inefanalysis = NULL;
485-
486 element->FindParam(&transfermethod,HydrologydcTransferFlagEnum);
487
488 /*Switch between the different transfer methods cases*/
489@@ -529,24 +515,9 @@
490 break;
491 case 1:
492 _assert_(sed_head_input);
493- _assert_(epl_head_input);
494-
495- inefanalysis = new HydrologyDCInefficientAnalysis();
496- hmax = inefanalysis->GetHydrologyDCInefficientHmax(element,gauss,thick_input,base_input);
497- delete inefanalysis;
498-
499 sed_head_input->GetInputValue(&sediment_head,gauss);
500- epl_head_input->GetInputValue(&epl_head,gauss);
501 element->FindParam(&leakage,HydrologydcLeakageFactorEnum);
502-
503- //Computing continuum function to apply to transfer term, transfer is null only if
504- // epl_head>sediment_head AND sediment_head>h_max
505- //let's try without that for a while
506- /* continuum=((1.0/(1.0+exp(-20.0*(sediment_head-epl_head)))))+(1.0/(1.0+exp(-20.0*(hmax-sediment_head)))); */
507- /* factor=min(continuum,1.0); */
508- /* transfer=sediment_head*leakage*factor; */
509 transfer=sediment_head*leakage;
510-
511 break;
512 default:
513 _error_("no case higher than 1 for the Transfer method");
514@@ -563,7 +534,6 @@
515 IssmDouble dt,A,B;
516 IssmDouble EPLgrad2;
517 IssmDouble EPL_N;
518-
519
520 femmodel->parameters->FindParam(&domaintype,DomainTypeEnum);
521
522@@ -618,7 +588,6 @@
523 element->GetInputListOnVertices(&bed[0],BaseEnum);
524
525 if(!active_element){
526-
527 /*Keeping thickness to initial value if EPL is not active*/
528 for(int i=0;i<numnodes;i++){
529 thickness[i]=init_thick;
530@@ -626,22 +595,15 @@
531 }
532 else{
533 for(int i=0;i<numnodes;i++){
534-
535 /*Compute first the effective pressure in the EPL*/
536- //EPL_N=gravity*((rho_ice*ice_thickness[i])-(rho_water*max(0.0,(eplhead[i]-bed[i]))));
537- //allowing EPLN larger than Pi
538 EPL_N=gravity*((rho_ice*ice_thickness[i])-(rho_water*(eplhead[i]-bed[i])));
539 if(EPL_N<0.0)EPL_N=0.0;
540 /*Get then the square of the gradient of EPL heads*/
541 EPLgrad2 = (epl_slopeX[i]*epl_slopeX[i])+(epl_slopeY[i]*epl_slopeY[i]);
542-
543 /*And proceed to the real thing*/
544- /* thickness[i] = old_thickness[i]*(1.0+((rho_water*gravity*dt)/(rho_ice*latentheat))*epl_conductivity*EPLgrad2- */
545- /* (2.0*(A*dt/(pow(n,n)))*(pow(EPL_N,n)))); */
546 thickness[i] = old_thickness[i]*
547 (2.0+((rho_water*gravity*dt)/(rho_ice*latentheat))*epl_conductivity*EPLgrad2-(2.0*(A*dt/(pow(n,n)))*(pow(EPL_N,n))))/
548 (2.0-((rho_water*gravity*dt)/(rho_ice*latentheat))*epl_conductivity*EPLgrad2+(2.0*(A*dt/(pow(n,n)))*(pow(EPL_N,n))));
549- //thickness[i] = 0.5*(thickness[i]+old_thickness[i]);
550 /*Take care of otherthikening*/
551 if(thickness[i]>max_thick){
552 thickness[i] = max_thick;
553@@ -694,7 +656,7 @@
554 int domaintype;
555 IssmDouble h_max;
556 IssmDouble sedheadmin;
557- Element* basalelement=NULL;
558+ Element* basalelement;
559
560 /*Get basal element*/
561 element->FindParam(&domaintype,DomainTypeEnum);
562@@ -703,7 +665,7 @@
563 basalelement = element;
564 break;
565 case Domain3DEnum:
566- if(!element->IsOnBase()) return;
567+ if(!element->IsOnBase()) return NULL;
568 basalelement = element->SpawnBasalElement();
569 break;
570 default: _error_("mesh "<<EnumToStringx(domaintype)<<" not supported yet");
571@@ -736,11 +698,12 @@
572 if (old_active[i]==0.){
573 epl_thickness[i]=init_thick;
574 }
575-
576 /*Now starting to look at the activations*/
577 if(residual[i]>0.){
578 vec_mask->SetValue(basalelement->nodes[i]->Sid(),1.,INS_VAL);
579- if(old_active[i]==0.) recurence->SetValue(basalelement->nodes[i]->Sid(),1.,INS_VAL);
580+ if(old_active[i]==0.){
581+ recurence->SetValue(basalelement->nodes[i]->Sid(),1.,INS_VAL);
582+ }
583 }
584 /*If mask was already one, keep one or colapse*/
585 else if(old_active[i]>0.){
586@@ -758,7 +721,9 @@
587 /*Increase of the domain is on the downstream node in term of sediment head*/
588 if(sedhead[j] == sedheadmin){
589 vec_mask->SetValue(basalelement->nodes[j]->Sid(),1.,INS_VAL);
590- if(old_active[i]==0.) recurence->SetValue(basalelement->nodes[i]->Sid(),1.,INS_VAL);
591+ if(old_active[i]==0.){
592+ recurence->SetValue(basalelement->nodes[i]->Sid(),1.,INS_VAL);
593+ }
594 }
595 }
596 }
597@@ -777,7 +742,7 @@
598 void HydrologyDCEfficientAnalysis::HydrologyEPLGetActive(Vector<IssmDouble>* active_vec, Element* element){
599 /*Constants*/
600 int domaintype;
601- Element* basalelement=NULL;
602+ Element* basalelement;
603
604 /*Get basal element*/
605 element->FindParam(&domaintype,DomainTypeEnum);
606@@ -786,7 +751,7 @@
607 basalelement = element;
608 break;
609 case Domain3DEnum:
610- if(!element->IsOnBase()) return;
611+ if(!element->IsOnBase()) return NULL;
612 basalelement = element->SpawnBasalElement();
613 break;
614 default: _error_("mesh "<<EnumToStringx(domaintype)<<" not supported yet");
615@@ -795,10 +760,10 @@
616 const int numnodes = basalelement->GetNumberOfNodes();
617 IssmDouble flag = 0.;
618 IssmDouble* active = xNew<IssmDouble>(numnodes);
619+ bool active_element;
620
621 /*Pass the activity mask from elements to nodes*/
622 basalelement->GetInputListOnVertices(&active[0],HydrologydcMaskEplactiveNodeEnum);
623- bool active_element;
624 Input* active_element_input=basalelement->GetInput(HydrologydcMaskEplactiveEltEnum); _assert_(active_element_input);
625 active_element_input->GetInputValue(&active_element);
626
Note: See TracBrowser for help on using the repository browser.