source: issm/oecreview/Archive/16554-17801/ISSM-16637-16638.diff@ 17802

Last change on this file since 17802 was 17802, checked in by Mathieu Morlighem, 11 years ago

Added archives
File size: 34.1 KB
RevLine 
[17802]1Index: ../trunk-jpl/src/c/analyses/HydrologyDCEfficientAnalysis.cpp
2===================================================================
3--- ../trunk-jpl/src/c/analyses/HydrologyDCEfficientAnalysis.cpp (revision 16637)
4+++ ../trunk-jpl/src/c/analyses/HydrologyDCEfficientAnalysis.cpp (revision 16638)
5@@ -56,6 +56,8 @@
6 iomodel->FetchDataToInput(elements,MeshElementonsurfaceEnum);
7 iomodel->FetchDataToInput(elements,BasalforcingsMeltingRateEnum);
8 iomodel->FetchDataToInput(elements,EplHeadEnum);
9+ iomodel->FetchDataToInput(elements,HydrologydcEplThicknessEnum);
10+
11 }/*}}}*/
12 void HydrologyDCEfficientAnalysis::CreateNodes(Nodes* nodes,IoModel* iomodel){/*{{{*/
13
14Index: ../trunk-jpl/src/c/shared/Enum/EnumDefinitions.h
15===================================================================
16--- ../trunk-jpl/src/c/shared/Enum/EnumDefinitions.h (revision 16637)
17+++ ../trunk-jpl/src/c/shared/Enum/EnumDefinitions.h (revision 16638)
18@@ -113,7 +113,8 @@
19 HydrologydcEplCompressibilityEnum,
20 HydrologydcEplPorosityEnum,
21 HydrologydcEplThicknessEnum,
22- HydrologydcEplTransmitivityEnum,
23+ HydrologydcEplThicknessOldEnum,
24+ HydrologydcEplConductivityEnum,
25 HydrologydcIsefficientlayerEnum,
26 HydrologydcSedimentlimitFlagEnum,
27 HydrologydcSedimentlimitEnum,
28Index: ../trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp
29===================================================================
30--- ../trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp (revision 16637)
31+++ ../trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp (revision 16638)
32@@ -121,7 +121,8 @@
33 case HydrologydcEplCompressibilityEnum : return "HydrologydcEplCompressibility";
34 case HydrologydcEplPorosityEnum : return "HydrologydcEplPorosity";
35 case HydrologydcEplThicknessEnum : return "HydrologydcEplThickness";
36- case HydrologydcEplTransmitivityEnum : return "HydrologydcEplTransmitivity";
37+ case HydrologydcEplThicknessOldEnum : return "HydrologydcEplThicknessOld";
38+ case HydrologydcEplConductivityEnum : return "HydrologydcEplConductivity";
39 case HydrologydcIsefficientlayerEnum : return "HydrologydcIsefficientlayer";
40 case HydrologydcSedimentlimitFlagEnum : return "HydrologydcSedimentlimitFlag";
41 case HydrologydcSedimentlimitEnum : return "HydrologydcSedimentlimit";
42Index: ../trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp
43===================================================================
44--- ../trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp (revision 16637)
45+++ ../trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp (revision 16638)
46@@ -121,7 +121,8 @@
47 else if (strcmp(name,"HydrologydcEplCompressibility")==0) return HydrologydcEplCompressibilityEnum;
48 else if (strcmp(name,"HydrologydcEplPorosity")==0) return HydrologydcEplPorosityEnum;
49 else if (strcmp(name,"HydrologydcEplThickness")==0) return HydrologydcEplThicknessEnum;
50- else if (strcmp(name,"HydrologydcEplTransmitivity")==0) return HydrologydcEplTransmitivityEnum;
51+ else if (strcmp(name,"HydrologydcEplThicknessOld")==0) return HydrologydcEplThicknessOldEnum;
52+ else if (strcmp(name,"HydrologydcEplConductivity")==0) return HydrologydcEplConductivityEnum;
53 else if (strcmp(name,"HydrologydcIsefficientlayer")==0) return HydrologydcIsefficientlayerEnum;
54 else if (strcmp(name,"HydrologydcSedimentlimitFlag")==0) return HydrologydcSedimentlimitFlagEnum;
55 else if (strcmp(name,"HydrologydcSedimentlimit")==0) return HydrologydcSedimentlimitEnum;
56@@ -135,11 +136,11 @@
57 else if (strcmp(name,"HydrologySedimentKmax")==0) return HydrologySedimentKmaxEnum;
58 else if (strcmp(name,"WaterTransfer")==0) return WaterTransferEnum;
59 else if (strcmp(name,"IndependentObject")==0) return IndependentObjectEnum;
60- else if (strcmp(name,"InversionControlParameters")==0) return InversionControlParametersEnum;
61 else stage=2;
62 }
63 if(stage==2){
64- if (strcmp(name,"InversionCostFunctionThreshold")==0) return InversionCostFunctionThresholdEnum;
65+ if (strcmp(name,"InversionControlParameters")==0) return InversionControlParametersEnum;
66+ else if (strcmp(name,"InversionCostFunctionThreshold")==0) return InversionCostFunctionThresholdEnum;
67 else if (strcmp(name,"InversionCostFunctionsCoefficients")==0) return InversionCostFunctionsCoefficientsEnum;
68 else if (strcmp(name,"InversionCostFunctions")==0) return InversionCostFunctionsEnum;
69 else if (strcmp(name,"InversionGradientScaling")==0) return InversionGradientScalingEnum;
70@@ -258,11 +259,11 @@
71 else if (strcmp(name,"SteadystateNumRequestedOutputs")==0) return SteadystateNumRequestedOutputsEnum;
72 else if (strcmp(name,"SteadystateReltol")==0) return SteadystateReltolEnum;
73 else if (strcmp(name,"SteadystateRequestedOutputs")==0) return SteadystateRequestedOutputsEnum;
74- else if (strcmp(name,"Surface")==0) return SurfaceEnum;
75 else stage=3;
76 }
77 if(stage==3){
78- if (strcmp(name,"SurfaceforcingsPrecipitation")==0) return SurfaceforcingsPrecipitationEnum;
79+ if (strcmp(name,"Surface")==0) return SurfaceEnum;
80+ else if (strcmp(name,"SurfaceforcingsPrecipitation")==0) return SurfaceforcingsPrecipitationEnum;
81 else if (strcmp(name,"SurfaceforcingsMassBalance")==0) return SurfaceforcingsMassBalanceEnum;
82 else if (strcmp(name,"SurfaceforcingsIspdd")==0) return SurfaceforcingsIspddEnum;
83 else if (strcmp(name,"SurfaceforcingsDesfac")==0) return SurfaceforcingsDesfacEnum;
84@@ -381,11 +382,11 @@
85 else if (strcmp(name,"DoubleParam")==0) return DoubleParamEnum;
86 else if (strcmp(name,"DoubleVecParam")==0) return DoubleVecParamEnum;
87 else if (strcmp(name,"Element")==0) return ElementEnum;
88- else if (strcmp(name,"ExternalResult")==0) return ExternalResultEnum;
89 else stage=4;
90 }
91 if(stage==4){
92- if (strcmp(name,"FileParam")==0) return FileParamEnum;
93+ if (strcmp(name,"ExternalResult")==0) return ExternalResultEnum;
94+ else if (strcmp(name,"FileParam")==0) return FileParamEnum;
95 else if (strcmp(name,"Input")==0) return InputEnum;
96 else if (strcmp(name,"IntInput")==0) return IntInputEnum;
97 else if (strcmp(name,"InputToExtrude")==0) return InputToExtrudeEnum;
98@@ -504,11 +505,11 @@
99 else if (strcmp(name,"BasalFriction")==0) return BasalFrictionEnum;
100 else if (strcmp(name,"ViscousHeating")==0) return ViscousHeatingEnum;
101 else if (strcmp(name,"QmuTemperature")==0) return QmuTemperatureEnum;
102- else if (strcmp(name,"HydrologyWaterVx")==0) return HydrologyWaterVxEnum;
103 else stage=5;
104 }
105 if(stage==5){
106- if (strcmp(name,"HydrologyWaterVy")==0) return HydrologyWaterVyEnum;
107+ if (strcmp(name,"HydrologyWaterVx")==0) return HydrologyWaterVxEnum;
108+ else if (strcmp(name,"HydrologyWaterVy")==0) return HydrologyWaterVyEnum;
109 else if (strcmp(name,"StressTensor")==0) return StressTensorEnum;
110 else if (strcmp(name,"StressTensorxx")==0) return StressTensorxxEnum;
111 else if (strcmp(name,"StressTensorxy")==0) return StressTensorxyEnum;
112Index: ../trunk-jpl/src/c/modules/InputDuplicatex/InputDuplicatex.cpp
113===================================================================
114--- ../trunk-jpl/src/c/modules/InputDuplicatex/InputDuplicatex.cpp (revision 16637)
115+++ ../trunk-jpl/src/c/modules/InputDuplicatex/InputDuplicatex.cpp (revision 16638)
116@@ -8,7 +8,6 @@
117 #include "../../toolkits/toolkits.h"
118
119 void InputDuplicatex(FemModel* femmodel,int original_enum, int new_enum){
120-
121 /*Go through elemnets, and ask to reinitialie the input: */
122 for(int i=0;i<femmodel->elements->Size();i++){
123 Element* element=dynamic_cast<Element*>(femmodel->elements->GetObjectByOffset(i));
124Index: ../trunk-jpl/src/c/solutionsequences/solutionsequence_hydro_nonlinear.cpp
125===================================================================
126--- ../trunk-jpl/src/c/solutionsequences/solutionsequence_hydro_nonlinear.cpp (revision 16637)
127+++ ../trunk-jpl/src/c/solutionsequences/solutionsequence_hydro_nonlinear.cpp (revision 16638)
128@@ -128,9 +128,12 @@
129 femmodel->HydrologyEPLupdateDomainx();
130 femmodel->parameters->SetParam(HydrologyEfficientEnum,HydrologyLayerEnum);
131
132+ /*Reset constraint on the ZigZag Lock*/
133+ ResetConstraintsx(femmodel);
134 /*Iteration on the EPL layer*/
135 eplconverged = false;
136 for(;;){
137+ femmodel->HydrologyEPLThicknessx();
138 femmodel->HydrologyTransferx();
139 SystemMatricesx(&Kff,&Kfs,&pf,&df,NULL,femmodel);
140 CreateNodalConstraintsx(&ys,femmodel->nodes,HydrologyDCEfficientAnalysisEnum);
141Index: ../trunk-jpl/src/c/cores/hydrology_core.cpp
142===================================================================
143--- ../trunk-jpl/src/c/cores/hydrology_core.cpp (revision 16637)
144+++ ../trunk-jpl/src/c/cores/hydrology_core.cpp (revision 16638)
145@@ -83,6 +83,7 @@
146 femmodel->parameters->FindParam(&isefficientlayer,HydrologydcIsefficientlayerEnum);
147 if (isefficientlayer){
148 InputDuplicatex(femmodel,EplHeadEnum,EplHeadOldEnum);
149+ InputDuplicatex(femmodel,HydrologydcEplThicknessEnum,HydrologydcEplThicknessOldEnum);
150 }
151
152 /*Proceed now to heads computations*/
153@@ -91,8 +92,8 @@
154 if(save_results && ((i+1)%output_frequency==0 || (i+1)==nsteps)){
155 if(VerboseSolution()) _printf0_(" saving results \n");
156 if(isefficientlayer){
157- int outputs[6] = {SedimentHeadEnum,SedimentHeadResidualEnum,EplHeadEnum,HydrologydcMaskEplactiveEnum,EplHeadSlopeXEnum,EplHeadSlopeYEnum};
158- femmodel->RequestedOutputsx(&femmodel->results,&outputs[0],6);
159+ int outputs[7] = {SedimentHeadEnum,SedimentHeadResidualEnum,EplHeadEnum,HydrologydcMaskEplactiveEnum,EplHeadSlopeXEnum,EplHeadSlopeYEnum,HydrologydcEplThicknessEnum};
160+ femmodel->RequestedOutputsx(&femmodel->results,&outputs[0],7);
161 }
162 else{
163 int outputs[2] = {SedimentHeadEnum,SedimentHeadResidualEnum};
164Index: ../trunk-jpl/src/c/classes/Materials/Matice.cpp
165===================================================================
166--- ../trunk-jpl/src/c/classes/Materials/Matice.cpp (revision 16637)
167+++ ../trunk-jpl/src/c/classes/Materials/Matice.cpp (revision 16638)
168@@ -144,6 +144,20 @@
169 return pow(B,-n);
170 }
171 /*}}}*/
172+/*FUNCTION Matice::GetAbar {{{*/
173+IssmDouble Matice::GetAbar(){
174+ /*
175+ * A = 1/B^n
176+ */
177+
178+ IssmDouble B,n;
179+
180+ inputs->GetInputAverage(&B,MaterialsRheologyBbarEnum);
181+ n=this->GetN();
182+
183+ return pow(B,-n);
184+}
185+/*}}}*/
186 /*FUNCTION Matice::GetB {{{*/
187 IssmDouble Matice::GetB(){
188
189Index: ../trunk-jpl/src/c/classes/Materials/Material.h
190===================================================================
191--- ../trunk-jpl/src/c/classes/Materials/Material.h (revision 16637)
192+++ ../trunk-jpl/src/c/classes/Materials/Material.h (revision 16638)
193@@ -35,6 +35,7 @@
194 virtual void GetViscosityDerivativeEpsSquareFS(IssmDouble* pmu_prime, IssmDouble* pepsilon)=0;
195 virtual void GetViscosity2dDerivativeEpsSquare(IssmDouble* pmu_prime, IssmDouble* pepsilon)=0;
196 virtual IssmDouble GetA()=0;
197+ virtual IssmDouble GetAbar()=0;
198 virtual IssmDouble GetB()=0;
199 virtual IssmDouble GetBbar()=0;
200 virtual IssmDouble GetN()=0;
201Index: ../trunk-jpl/src/c/classes/Materials/Matpar.cpp
202===================================================================
203--- ../trunk-jpl/src/c/classes/Materials/Matpar.cpp (revision 16637)
204+++ ../trunk-jpl/src/c/classes/Materials/Matpar.cpp (revision 16638)
205@@ -65,8 +65,7 @@
206 if(isefficientlayer){
207 iomodel->Constant(&this->epl_compressibility,HydrologydcEplCompressibilityEnum);
208 iomodel->Constant(&this->epl_porosity,HydrologydcEplPorosityEnum);
209- iomodel->Constant(&this->epl_thickness,HydrologydcEplThicknessEnum);
210- iomodel->Constant(&this->epl_transmitivity,HydrologydcEplTransmitivityEnum);
211+ iomodel->Constant(&this->epl_conductivity,HydrologydcEplConductivityEnum);
212 }
213 }
214 else{
215@@ -353,9 +352,9 @@
216 ( this->water_compressibility+( this->sediment_compressibility/ this->sediment_porosity));
217 }
218 /*}}}*/
219-/*FUNCTION Matpar::GetEplStoring {{{*/
220-IssmDouble Matpar::GetEplStoring(){
221- return this->rho_freshwater* this->g* this->epl_porosity* this->epl_thickness*
222+/*FUNCTION Matpar::GetEplSpecificStoring {{{*/
223+IssmDouble Matpar::GetEplSpecificStoring(){
224+ return this->rho_freshwater* this->g* this->epl_porosity*
225 ( this->water_compressibility+( this->epl_compressibility/ this->epl_porosity));
226 }
227 /*}}}*/
228@@ -368,10 +367,10 @@
229 IssmDouble Matpar::GetSedimentThickness(){
230 return sediment_thickness;
231 }
232-/*}}}*/
233-/*FUNCTION Matpar::GetEplTransitivity {{{*/
234-IssmDouble Matpar::GetEplTransmitivity(){
235- return epl_transmitivity;
236+/*}}}*/
237+/*FUNCTION Matpar::GetEplConductivity {{{*/
238+IssmDouble Matpar::GetEplConductivity(){
239+ return epl_conductivity;
240 }
241 /*}}}*/
242 /*FUNCTION Matpar::TMeltingPoint {{{*/
243Index: ../trunk-jpl/src/c/classes/Materials/Matice.h
244===================================================================
245--- ../trunk-jpl/src/c/classes/Materials/Matice.h (revision 16637)
246+++ ../trunk-jpl/src/c/classes/Materials/Matice.h (revision 16638)
247@@ -62,6 +62,7 @@
248 void GetViscosityDerivativeEpsSquareFS(IssmDouble* pmu_prime, IssmDouble* pepsilon);
249 void GetViscosity2dDerivativeEpsSquare(IssmDouble* pmu_prime, IssmDouble* pepsilon);
250 IssmDouble GetA();
251+ IssmDouble GetAbar();
252 IssmDouble GetB();
253 IssmDouble GetBbar();
254 IssmDouble GetD();
255Index: ../trunk-jpl/src/c/classes/Materials/Matpar.h
256===================================================================
257--- ../trunk-jpl/src/c/classes/Materials/Matpar.h (revision 16637)
258+++ ../trunk-jpl/src/c/classes/Materials/Matpar.h (revision 16638)
259@@ -48,8 +48,7 @@
260
261 IssmDouble epl_compressibility;
262 IssmDouble epl_porosity;
263- IssmDouble epl_thickness;
264- IssmDouble epl_transmitivity;
265+ IssmDouble epl_conductivity;
266
267 /*gia: */
268 IssmDouble lithosphere_shear_modulus;
269@@ -93,6 +92,7 @@
270 void GetViscosityDerivativeEpsSquareFS(IssmDouble* pmu_prime, IssmDouble* pepsilon){_error_("not supported");};
271 void GetViscosity2dDerivativeEpsSquare(IssmDouble* pmu_prime, IssmDouble* pepsilon){_error_("not supported");};
272 IssmDouble GetA(){_error_("not supported");};
273+ IssmDouble GetAbar(){_error_("not supported");};
274 IssmDouble GetB(){_error_("not supported");};
275 IssmDouble GetBbar(){_error_("not supported");};
276 IssmDouble GetN(){_error_("not supported");};
277@@ -120,10 +120,10 @@
278 IssmDouble GetHydrologyCR();
279 IssmDouble GetHydrologyN();
280 IssmDouble GetSedimentStoring();
281- IssmDouble GetEplStoring();
282+ IssmDouble GetEplSpecificStoring();
283 IssmDouble GetSedimentTransmitivity();
284 IssmDouble GetSedimentThickness();
285- IssmDouble GetEplTransmitivity();
286+ IssmDouble GetEplConductivity();
287 IssmDouble TMeltingPoint(IssmDouble pressure);
288 IssmDouble PureIceEnthalpy(IssmDouble pressure);
289 IssmDouble GetEnthalpyDiffusionParameter(IssmDouble enthalpy,IssmDouble pressure);
290Index: ../trunk-jpl/src/c/classes/Elements/Tria.cpp
291===================================================================
292--- ../trunk-jpl/src/c/classes/Elements/Tria.cpp (revision 16637)
293+++ ../trunk-jpl/src/c/classes/Elements/Tria.cpp (revision 16638)
294@@ -1577,8 +1577,9 @@
295 void Tria::InputDuplicate(int original_enum,int new_enum){
296
297 /*Call inputs method*/
298- if (IsInput(original_enum)) inputs->DuplicateInput(original_enum,new_enum);
299-
300+ if (IsInput(original_enum)) {
301+ inputs->DuplicateInput(original_enum,new_enum);
302+ }
303 }
304 /*}}}*/
305 /*FUNCTION Tria::InputScale{{{*/
306@@ -2096,6 +2097,8 @@
307 name==SedimentHeadEnum ||
308 name==EplHeadOldEnum ||
309 name==EplHeadEnum ||
310+ name==HydrologydcEplThicknessOldEnum ||
311+ name==HydrologydcEplThicknessEnum ||
312 name==HydrologydcMaskEplactiveEnum ||
313 name==MeshVertexonbedEnum ||
314 name==WaterTransferEnum ||
315@@ -2454,10 +2457,6 @@
316 this->ComputeStressTensor();
317 input=this->inputs->GetInput(output_enum);
318 break;
319- case SurfaceNormalVelocityEnum:
320- this->ComputeSurfaceNormalVelocity();
321- input=this->inputs->GetInput(output_enum);
322- break;
323 default:
324 _error_("input "<<EnumToStringx(output_enum)<<" not found in element");
325 }
326@@ -2489,9 +2488,6 @@
327 for(int i=0;i<NUMVERTICES;i++) values[i] = values[i]/reCast<IssmDouble>(connectivity[i]);
328
329 vector->SetValues(NUMVERTICES,&sidlist[0],&values[0],ADD_VAL);
330-
331- /*FIXME: Inputs needs to be updated first in Tria::ResultInterpolation, this is a hack*/
332- this->inputs->DeleteInput(SurfaceNormalVelocityEnum);
333 break;
334 }
335 default:
336@@ -2650,53 +2646,6 @@
337 *(surface_normal+2) = normal[2]/normal_norm;
338 }
339 /*}}}*/
340-/*FUNCTION Tria::ComputeSurfaceNormalVelocity{{{*/
341-void Tria::ComputeSurfaceNormalVelocity(){
342-
343- IssmDouble sum,tangential_vector[2],normal_vector[2],time,ux,uy;
344- IssmDouble normal_velocity[NUMVERTICES],xyz_list[NUMVERTICES][3];
345- IssmDouble value[NUMVERTICES],verticesonsurface[NUMVERTICES];
346-
347- for(int iv=0;iv<NUMVERTICES;iv++){
348- normal_velocity[iv]=0.;
349- value[iv]=0.;
350- }
351-
352- GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
353-
354- GaussTria* gauss=new GaussTria();
355- Input* slope_input=inputs->GetInput(SurfaceSlopeXEnum); _assert_(slope_input);
356- // Input* slope_input= inputs->GetInput(SurfaceEnum); _assert_(slope_input);
357- Input* vx_input=inputs->GetInput(VxEnum); _assert_(vx_input);
358- Input* vy_input=inputs->GetInput(VyEnum); _assert_(vy_input);
359-
360-
361- /*Get list of nodes on surface*/
362- GetInputListOnVertices(&verticesonsurface[0],MeshVertexonsurfaceEnum);
363- sum = verticesonsurface[0]+verticesonsurface[1]+verticesonsurface[2];
364- _assert_(sum==0. || sum==1. || sum==2.);
365-
366- /*Compute normal velocity for surface nodes from L2 projected slope*/
367- for(int iv=0;iv<NUMVERTICES;iv++){
368- if(verticesonsurface[iv] == 1){
369- gauss->GaussNode(P1Enum,iv);
370- slope_input->GetInputValue(&value[iv],gauss);
371- // slope_input->GetInputDerivativeValue(&value[iv],&xyz_list[0][0],gauss);
372- vx_input->GetInputValue(&ux,gauss);
373- vy_input->GetInputValue(&uy,gauss);
374- tangential_vector[0]=sqrt(1./(pow(value[iv],2.)+1.));
375- tangential_vector[1]=value[iv]*tangential_vector[0];
376- normal_vector[0]=-1.*tangential_vector[1];
377- normal_vector[1]=tangential_vector[0];
378- normal_velocity[iv]=ux*normal_vector[0]+uy*normal_vector[1];
379- }
380- }
381-
382- delete gauss;
383- this->inputs->AddInput(new TriaInput(SurfaceNormalVelocityEnum,&normal_velocity[0],P1Enum));
384-
385-}
386-/*}}}*/
387 /*FUNCTION Tria::TimeAdapt{{{*/
388 IssmDouble Tria::TimeAdapt(void){
389
390@@ -3277,7 +3226,7 @@
391 /*Now get the average SMB over the element*/
392 Input* smb_input = inputs->GetInput(SurfaceforcingsMassBalanceEnum); _assert_(smb_input);
393 smb_input->GetInputAverage(&smb); // average smb on element in m ice s-1
394- Total_Smb=rho_ice*base*smb; // smb on element in kg s-1
395+ Total_Smb=rho_ice*base*smb; // smb on element in kg s-1
396
397 /*Return: */
398 return Total_Smb;
399@@ -6733,9 +6682,10 @@
400 ElementMatrix* Tria::CreateKMatrixHydrologyDCEfficient(void){
401
402 /* Intermediaries */
403- IssmDouble D_scalar,Jdet;
404- IssmDouble epl_transmitivity,dt;
405- IssmDouble epl_storing;
406+ IssmDouble D_scalar,Jdet,dt;
407+ IssmDouble epl_thickness;
408+ IssmDouble epl_conductivity;
409+ IssmDouble epl_specificstoring;
410 IssmDouble xyz_list[NUMVERTICES][3];
411
412 /*Check that all nodes are active, else return empty matrix*/
413@@ -6754,19 +6704,23 @@
414
415 /*Retrieve all inputs and parameters*/
416 GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
417+ Input* thickness_input=inputs->GetInput(HydrologydcEplThicknessEnum); _assert_(thickness);
418 this->parameters->FindParam(&dt,TimesteppingTimeStepEnum);
419- epl_storing = matpar->GetEplStoring();
420- epl_transmitivity = matpar->GetEplTransmitivity();
421+ epl_specificstoring = matpar->GetEplSpecificStoring();
422+ epl_conductivity = matpar->GetEplConductivity();
423
424+
425 /* Start looping on the number of gaussian points: */
426 GaussTria* gauss=new GaussTria(2);
427 for(int ig=gauss->begin();ig<gauss->end();ig++){
428-
429+
430+
431 gauss->GaussPoint(ig);
432 GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss);
433+ thickness_input->GetInputValue(&epl_thickness,gauss);
434
435 /*Diffusivity*/
436- D_scalar=epl_transmitivity*gauss->weight*Jdet;
437+ D_scalar=epl_conductivity*epl_thickness*gauss->weight*Jdet;
438 if(reCast<bool,IssmDouble>(dt)) D_scalar=D_scalar*dt;
439 D[0][0]=D_scalar; D[0][1]=0.;
440 D[1][0]=0.; D[1][1]=D_scalar;
441@@ -6779,7 +6733,7 @@
442 /*Transient*/
443 if(reCast<bool,IssmDouble>(dt)){
444 GetNodalFunctions(basis,gauss);
445- D_scalar=epl_storing*gauss->weight*Jdet;
446+ D_scalar=epl_specificstoring*epl_thickness*gauss->weight*Jdet;
447
448 TripleMultiply(basis,numnodes,1,0,
449 &D_scalar,1,1,0,
450@@ -6915,7 +6869,8 @@
451 IssmDouble xyz_list[NUMVERTICES][3];
452 IssmDouble dt,scalar,water_head;
453 IssmDouble transfer,residual;
454- IssmDouble epl_storing;
455+ IssmDouble epl_thickness;
456+ IssmDouble epl_specificstoring;
457 GaussTria* gauss = NULL;
458
459 /*Check that all nodes are active, else return empty matrix*/
460@@ -6931,11 +6886,12 @@
461 IssmDouble* basis = xNew<IssmDouble>(numnodes);
462
463 /*Retrieve all inputs and parameters*/
464- epl_storing = matpar->GetEplStoring();
465+ epl_specificstoring = matpar->GetEplSpecificStoring();
466 GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
467 this->parameters->FindParam(&dt,TimesteppingTimeStepEnum);
468- Input* residual_input=inputs->GetInput(SedimentHeadResidualEnum); _assert_(residual_input);
469- Input* transfer_input=inputs->GetInput(WaterTransferEnum); _assert_(transfer_input);
470+ Input* residual_input=inputs->GetInput(SedimentHeadResidualEnum); _assert_(residual_input);
471+ Input* transfer_input=inputs->GetInput(WaterTransferEnum); _assert_(transfer_input);
472+ Input* thickness_input=inputs->GetInput(HydrologydcEplThicknessEnum); _assert_(thickness_input);
473 Input* old_wh_input=NULL;
474
475 if(reCast<bool,IssmDouble>(dt)){
476@@ -6958,8 +6914,9 @@
477
478 /*Transient term*/
479 if(reCast<bool,IssmDouble>(dt)){
480+ thickness_input->GetInputValue(&epl_thickness,gauss);
481 old_wh_input->GetInputValue(&water_head,gauss);
482- scalar = Jdet*gauss->weight*water_head*epl_storing;
483+ scalar = Jdet*gauss->weight*water_head*epl_specificstoring*epl_thickness;
484 for(int i=0;i<numnodes;i++) pe->values[i]+=scalar*basis[i];
485 }
486 }
487@@ -7141,8 +7098,8 @@
488 IssmDouble sed_trans,sed_thick;
489 IssmDouble leakage,h_max;
490 IssmDouble wh_trans;
491- IssmDouble activeEpl[numdof];
492- IssmDouble eplstoring[numdof],sedstoring[numdof];
493+ IssmDouble activeEpl[numdof],epl_thickness[numdof];
494+ IssmDouble epl_specificstoring[numdof],sedstoring[numdof];
495 IssmDouble epl_head[numdof],sed_head[numdof];
496
497 GetDofList(&doflist,NoneApproximationEnum,GsetEnum);
498@@ -7164,6 +7121,7 @@
499 GetInputListOnVertices(&activeEpl[0],HydrologydcMaskEplactiveEnum);
500 GetInputListOnVertices(&sed_head[0],SedimentHeadEnum);
501 GetInputListOnVertices(&epl_head[0],EplHeadEnum);
502+ GetInputListOnVertices(&epl_thickness[0],HydrologydcEplThicknessEnum);
503
504 this->parameters->FindParam(&leakage,HydrologydcLeakageFactorEnum);
505
506@@ -7176,12 +7134,12 @@
507 wh_trans=0.0;
508 }
509 else{
510- eplstoring[i]=matpar->GetEplStoring();
511+ epl_specificstoring[i]=matpar->GetEplSpecificStoring();
512 sedstoring[i]=matpar->GetSedimentStoring();
513
514 /*EPL head higher than sediment head, transfer from the epl to the sediment*/
515 if(epl_head[i]>sed_head[i]){
516- wh_trans=eplstoring[i]*sed_trans*(epl_head[i]-sed_head[i])/(leakage*sed_thick);
517+ wh_trans=epl_specificstoring[i]*epl_thickness[i]*sed_trans*(epl_head[i]-sed_head[i])/(leakage*sed_thick);
518
519 /*No transfer if the sediment head is allready at the maximum*/
520 this->GetHydrologyDCInefficientHmax(&h_max,nodes[i]);
521@@ -7245,7 +7203,7 @@
522 GetInputListOnVertices(&eplhead[0],EplHeadEnum);
523 GetInputListOnVertices(&residual[0],SedimentHeadResidualEnum);
524
525- /*Get minimum sediment head*/
526+ /*Get minimum sediment head of the element*/
527 sedheadmin=sedhead[0];
528 for(i=1;i<numdof;i++) if(sedhead[i]<=sedheadmin)sedheadmin=sedhead[i];
529
530@@ -7279,6 +7237,63 @@
531 /*}}}*/
532 /*FUNCTION Tria::ComputeEPLThickness{{{*/
533 void Tria::ComputeEPLThickness(void){
534+
535+ int i;
536+ const int numdof = NDOF1 *NUMVERTICES;
537+ bool isefficientlayer;
538+ IssmDouble n,A,dt;
539+ IssmDouble rho_water,rho_ice;
540+ IssmDouble gravity,latentheat,EPLgrad;
541+ IssmDouble EPL_N,epl_conductivity;
542+ IssmDouble activeEpl[numdof],thickness[numdof];
543+ IssmDouble eplhead[numdof], old_thickness[numdof];
544+ IssmDouble epl_slopeX[numdof],epl_slopeY[numdof];
545+ IssmDouble ice_thickness[numdof],bed[numdof];
546+
547+
548+ /*Get the flag to know if the efficient layer is present*/
549+ this->parameters->FindParam(&isefficientlayer,HydrologydcIsefficientlayerEnum);
550+ this->parameters->FindParam(&dt,TimesteppingTimeStepEnum);
551+
552+ if(isefficientlayer){
553+ /*For now, assuming just one way to compute EPL thickness*/
554+ rho_water = matpar->GetRhoWater();
555+ rho_ice = matpar->GetRhoIce();
556+ gravity = matpar->GetG();
557+ latentheat = matpar->GetLatentHeat();
558+ n = material->GetN();
559+ A = material->GetAbar();
560+
561+ GetInputListOnVertices(&activeEpl[0],HydrologydcMaskEplactiveEnum);
562+ GetInputListOnVertices(&eplhead[0],EplHeadEnum);
563+ GetInputListOnVertices(&epl_slopeX[0],EplHeadSlopeXEnum);
564+ GetInputListOnVertices(&epl_slopeY[0],EplHeadSlopeYEnum);
565+ GetInputListOnVertices(&old_thickness[0],HydrologydcEplThicknessOldEnum);
566+ GetInputListOnVertices(&ice_thickness[0],ThicknessEnum);
567+ GetInputListOnVertices(&bed[0],BedEnum);
568+
569+ epl_conductivity = matpar->GetEplConductivity();
570+
571+ for(int i=0;i<numdof;i++){
572+ /*Keeping thickness to 1 if EPL is not active*/
573+ if(activeEpl[i]==0.0){
574+ thickness[i]=1.0;
575+ }
576+ else{
577+
578+ /*Compute first the effective pressure in the EPL*/
579+ EPL_N=gravity*((rho_ice*ice_thickness[i])-(rho_water*(eplhead[i]-bed[i])));
580+
581+ /*Get then the gradient of EPL heads*/
582+ EPLgrad = epl_slopeX[i]+epl_slopeY[i];
583+
584+ /*And proceed to the real thing*/
585+ thickness[i] = old_thickness[i]*(1-((rho_water*gravity*dt)/(rho_ice*latentheat))*epl_conductivity*pow(EPLgrad,2.0)-2.0*(A*dt/(pow(n,n)))*EPL_N);
586+
587+ }
588+ }
589+ this->inputs->AddInput(new TriaInput(HydrologydcEplThicknessEnum,thickness,P1Enum));
590+ }
591 }
592 /*}}}*/
593
594Index: ../trunk-jpl/src/c/classes/Elements/Penta.cpp
595===================================================================
596--- ../trunk-jpl/src/c/classes/Elements/Penta.cpp (revision 16637)
597+++ ../trunk-jpl/src/c/classes/Elements/Penta.cpp (revision 16638)
598@@ -2644,6 +2644,8 @@
599 name==EplHeadEnum ||
600 name==SedimentHeadOldEnum ||
601 name==EplHeadOldEnum ||
602+ name==HydrologydcEplThicknessOldEnum ||
603+ name==HydrologydcEplThicknessEnum ||
604 name==HydrologydcMaskEplactiveEnum ||
605 name==WaterTransferEnum
606
607@@ -2870,7 +2872,7 @@
608
609 switch(input->GetResultInterpolation()){
610 case P0Enum:
611- _error_("not implemented...");
612+ _error_("P0 not implemented yet for input "<<EnumToStringx(output_enum));
613 break;
614 case P1Enum:{
615 IssmDouble values[NUMVERTICES];
616@@ -10841,12 +10843,77 @@
617 /*FUNCTION Penta::ComputeEPLThickness{{{*/
618 void Penta::ComputeEPLThickness(void){
619
620+ int i;
621+ const int numdof = NDOF1 *NUMVERTICES;
622+ const int numdof2d = NDOF1 *NUMVERTICES2D;
623+ bool isefficientlayer;
624+ IssmDouble n,A,dt;
625+ IssmDouble rho_water,rho_ice;
626+ IssmDouble gravity,latentheat,EPLgrad;
627+ IssmDouble EPL_N,epl_conductivity;
628+ IssmDouble activeEpl[numdof],thickness[numdof];
629+ IssmDouble eplhead[numdof], old_thickness[numdof];
630+ IssmDouble epl_slopeX[numdof],epl_slopeY[numdof];
631+ IssmDouble ice_thickness[numdof],bed[numdof];
632+ Penta *penta = NULL;
633+ /*If not on bed, return*/
634 if (!IsOnBed())return;
635
636- Tria* tria=(Tria*)SpawnTria(0); //lower face is 0, upper face is 1.
637- tria->ComputeEPLThickness();
638- delete tria->material; delete tria;
639+ /*Get the flag to know if the efficient layer is present*/
640+ this->parameters->FindParam(&isefficientlayer,HydrologydcIsefficientlayerEnum);
641+ this->parameters->FindParam(&dt,TimesteppingTimeStepEnum);
642
643+ if(isefficientlayer){
644+ /*For now, assuming just one way to compute EPL thickness*/
645+ rho_water = matpar->GetRhoWater();
646+ rho_ice = matpar->GetRhoIce();
647+ gravity = matpar->GetG();
648+ latentheat = matpar->GetLatentHeat();
649+ epl_conductivity = matpar->GetEplConductivity();
650+ n = material->GetN();
651+ A = material->GetA();
652+
653+
654+ GetInputListOnVertices(&activeEpl[0],HydrologydcMaskEplactiveEnum);
655+ GetInputListOnVertices(&eplhead[0],EplHeadEnum);
656+ GetInputListOnVertices(&epl_slopeX[0],EplHeadSlopeXEnum);
657+ GetInputListOnVertices(&epl_slopeY[0],EplHeadSlopeYEnum);
658+ GetInputListOnVertices(&old_thickness[0],HydrologydcEplThicknessOldEnum);
659+ GetInputListOnVertices(&ice_thickness[0],ThicknessEnum);
660+ GetInputListOnVertices(&bed[0],BedEnum);
661+
662+ for(int i=0;i<numdof2d;i++){
663+ /*Keeping thickness to 1 if EPL is not active*/
664+ if(activeEpl[i]==0.0){
665+ thickness[i]=1.0;
666+ thickness[i+numdof2d]=thickness[i];
667+ }
668+ else{
669+
670+ /*Compute first the effective pressure in the EPL*/
671+ EPL_N=gravity*((rho_ice*ice_thickness[i])-(rho_water*(eplhead[i]-bed[i])));
672+
673+ /*Get then the gradient of EPL heads*/
674+ EPLgrad = epl_slopeX[i]+epl_slopeY[i];
675+
676+ /*And proceed to the real thing*/
677+ thickness[i] = old_thickness[i]*(1-((rho_water*gravity*dt)/(rho_ice* latentheat))*epl_conductivity*pow(EPLgrad,2.0)-2.0*(A*dt/(pow(n,n)))*EPL_N);
678+ thickness[i+numdof2d]=thickness[i];
679+ }
680+ }
681+ penta=this;
682+ for(;;){
683+
684+ /*Add input to the element: */
685+ penta->inputs->AddInput(new PentaInput(HydrologydcEplThicknessEnum,thickness,P1Enum));
686+
687+ /*Stop if we have reached the surface*/
688+ if (penta->IsOnSurface()) break;
689+
690+ /* get upper Penta*/
691+ penta=penta->GetUpperElement(); _assert_(penta->Id()!=this->id);
692+ }
693+ }
694 }
695 /*}}}*/
696 /*FUNCTION Penta::InputUpdateFromSolutionHydrologyDCInefficient{{{*/
697Index: ../trunk-jpl/src/m/classes/hydrologydc.m
698===================================================================
699--- ../trunk-jpl/src/m/classes/hydrologydc.m (revision 16637)
700+++ ../trunk-jpl/src/m/classes/hydrologydc.m (revision 16638)
701@@ -26,7 +26,7 @@
702 epl_compressibility = 0;
703 epl_porosity = 0;
704 epl_thickness = 0;
705- epl_transmitivity = 0;
706+ epl_conductivity = 0;
707 end
708 methods
709 % {{{ function obj = hydrologydc(varargin)
710@@ -60,7 +60,7 @@
711 obj.epl_compressibility = 1.0e-08;
712 obj.epl_porosity = 0.4;
713 obj.epl_thickness = 1.0;
714- obj.epl_transmitivity = 8.0e-02;
715+ obj.epl_conductivity = 8.0e-02;
716
717 end
718 % }}}
719@@ -96,7 +96,7 @@
720 md = checkfield(md,'hydrology.epl_compressibility','>',0,'numel',1);
721 md = checkfield(md,'hydrology.epl_porosity','>',0,'numel',1);
722 md = checkfield(md,'hydrology.epl_thickness','>',0,'numel',1);
723- md = checkfield(md,'hydrology.epl_transmitivity','>',0,'numel',1);
724+ md = checkfield(md,'hydrology.epl_conductivity','>',0,'numel',1);
725 end
726 end
727 % }}}
728@@ -136,8 +136,8 @@
729 fielddisplay(obj,'mask_eplactive','active (1) or not (0) EPL');
730 fielddisplay(obj,'epl_compressibility','epl compressibility [Pa^-1]');
731 fielddisplay(obj,'epl_porosity','epl [dimensionless]');
732- fielddisplay(obj,'epl_thickness','epl thickness [m]');
733- fielddisplay(obj,'epl_transmitivity','epl transmitivity [m^2/s]');
734+ fielddisplay(obj,'epl_thickness','epl initial thickness [m]');
735+ fielddisplay(obj,'epl_conductivity','epl conductivity [m^2/s]');
736 end
737
738 end
739@@ -171,7 +171,8 @@
740 WriteData(fid,'object',obj,'fieldname','epl_compressibility','format','Double');
741 WriteData(fid,'object',obj,'fieldname','epl_porosity','format','Double');
742 WriteData(fid,'object',obj,'fieldname','epl_thickness','format','Double');
743- WriteData(fid,'object',obj,'fieldname','epl_transmitivity','format','Double');
744+ %WriteData(fid,'object',obj,'fieldname','epl_transmitivity','format','Double');
745+ WriteData(fid,'object',obj,'fieldname','epl_conductivity','format','Double');
746 end
747 end
748 % }}}
749Index: ../trunk-jpl/src/m/enum/HydrologydcEplConductivityEnum.m
750===================================================================
751--- ../trunk-jpl/src/m/enum/HydrologydcEplConductivityEnum.m (revision 0)
752+++ ../trunk-jpl/src/m/enum/HydrologydcEplConductivityEnum.m (revision 16638)
753@@ -0,0 +1,11 @@
754+function macro=HydrologydcEplConductivityEnum()
755+%HYDROLOGYDCEPLCONDUCTIVITYENUM - Enum of HydrologydcEplConductivity
756+%
757+% WARNING: DO NOT MODIFY THIS FILE
758+% this file has been automatically generated by src/c/shared/Enum/Synchronize.sh
759+% Please read src/c/shared/Enum/README for more information
760+%
761+% Usage:
762+% macro=HydrologydcEplConductivityEnum()
763+
764+macro=StringToEnum('HydrologydcEplConductivity');
765Index: ../trunk-jpl/src/m/enum/EnumDefinitions.py
766===================================================================
767--- ../trunk-jpl/src/m/enum/EnumDefinitions.py (revision 16637)
768+++ ../trunk-jpl/src/m/enum/EnumDefinitions.py (revision 16638)
769@@ -113,7 +113,8 @@
770 def HydrologydcEplCompressibilityEnum(): return StringToEnum("HydrologydcEplCompressibility")[0]
771 def HydrologydcEplPorosityEnum(): return StringToEnum("HydrologydcEplPorosity")[0]
772 def HydrologydcEplThicknessEnum(): return StringToEnum("HydrologydcEplThickness")[0]
773-def HydrologydcEplTransmitivityEnum(): return StringToEnum("HydrologydcEplTransmitivity")[0]
774+def HydrologydcEplThicknessOldEnum(): return StringToEnum("HydrologydcEplThicknessOld")[0]
775+def HydrologydcEplConductivityEnum(): return StringToEnum("HydrologydcEplConductivity")[0]
776 def HydrologydcIsefficientlayerEnum(): return StringToEnum("HydrologydcIsefficientlayer")[0]
777 def HydrologydcSedimentlimitFlagEnum(): return StringToEnum("HydrologydcSedimentlimitFlag")[0]
778 def HydrologydcSedimentlimitEnum(): return StringToEnum("HydrologydcSedimentlimit")[0]
779Index: ../trunk-jpl/src/m/enum/HydrologydcEplThicknessOldEnum.m
780===================================================================
781--- ../trunk-jpl/src/m/enum/HydrologydcEplThicknessOldEnum.m (revision 0)
782+++ ../trunk-jpl/src/m/enum/HydrologydcEplThicknessOldEnum.m (revision 16638)
783@@ -0,0 +1,11 @@
784+function macro=HydrologydcEplThicknessOldEnum()
785+%HYDROLOGYDCEPLTHICKNESSOLDENUM - Enum of HydrologydcEplThicknessOld
786+%
787+% WARNING: DO NOT MODIFY THIS FILE
788+% this file has been automatically generated by src/c/shared/Enum/Synchronize.sh
789+% Please read src/c/shared/Enum/README for more information
790+%
791+% Usage:
792+% macro=HydrologydcEplThicknessOldEnum()
793+
794+macro=StringToEnum('HydrologydcEplThicknessOld');
Note: See TracBrowser for help on using the repository browser.