Changeset 22236

Timestamp:

11/05/17 08:11:00 (7 years ago)

Author:

bdef

Message:

BUG:Change in the initialization of epl head

File:

• issm/trunk-jpl/src/c/analyses/HydrologyDCEfficientAnalysis.cpp (modified) (16 diffs)

issm/trunk-jpl/src/c/analyses/HydrologyDCEfficientAnalysis.cpp

@@ -240 +240 @@
                 /*Diffusivity*/
                 D_scalar=epl_transmitivity*gauss->weight*Jdet;
+                //D_scalar=gauss->weight*Jdet;
                 if(dt!=0.) D_scalar=D_scalar*dt;
                 D[0][0]=D_scalar;
@@ -253 +254 @@
                         basalelement->NodalFunctions(&basis[0],gauss);
                         D_scalar=epl_storing*gauss->weight*Jdet;
+                        //D_scalar=(epl_storing/epl_transmitivity)*gauss->weight*Jdet;
                         TripleMultiply(basis,numnodes,1,0,
                                                 &D_scalar,1,1,0,
@@ -262 +264 @@
                         transfer=GetHydrologyKMatrixTransfer(basalelement);
                         D_scalar=dt*transfer*gauss->weight*Jdet;
+                        //D_scalar=dt*(transfer/epl_transmitivity)*gauss->weight*Jdet;
                         TripleMultiply(basis,numnodes,1,0,
                                                                                  &D_scalar,1,1,0,
@@ -312 +315 @@
         IssmDouble dt,scalar,water_head;
         IssmDouble water_load,transfer;
-        IssmDouble epl_storing;
+        IssmDouble epl_storing,epl_transmitivity;
         IssmDouble Jdet;
         IssmDouble residual,connectivity;
@@ -348 +351 @@
                 basalelement ->NodalFunctions(basis,gauss);
                 epl_storing     = EplStoring(basalelement,gauss,epl_thick_input,epl_head_input,base_input);
+                epl_transmitivity = EplTransmitivity(basalelement,gauss,epl_thick_input,epl_head_input,base_input);
+
                 /*Loading term*/
                 water_input->GetInputValue(&water_load,gauss);
                 scalar = Jdet*gauss->weight*(water_load);
+                //scalar = Jdet*gauss->weight*(water_load)/epl_transmitivity;
                 if(dt!=0.) scalar = scalar*dt;
                 for(int i=0;i<numnodes;i++)pe->values[i]+=scalar*basis[i];
@@ -360 +366 @@
                         transfer=GetHydrologyPVectorTransfer(basalelement,gauss,sed_head_input);
                         scalar = Jdet*gauss->weight*((water_head*epl_storing)+(dt*transfer));
+                        //scalar = Jdet*gauss->weight*((water_head*epl_storing)+(dt*transfer))/epl_transmitivity;
                         for(int i=0;i<numnodes;i++)pe->values[i]+=scalar*basis[i];
                 }
@@ -369 +376 @@
         for(int iv=0;iv<numvertices;iv++){
                 gauss->GaussVertex(iv);
+                epl_transmitivity = EplTransmitivity(basalelement,gauss,epl_thick_input,epl_head_input,base_input);
                 connectivity = IssmDouble(basalelement->VertexConnectivity(iv));
                 residual_input->GetInputValue(&residual,gauss);
                 pe->values[iv]+=residual/connectivity;
+                //pe->values[iv]+=residual/(epl_transmitivity*connectivity);
         }
         /*Clean up and return*/
@@ -415 +424 @@
         /*Fetch dof list and allocate solution vector*/
         IssmDouble* eplHeads    = xNew<IssmDouble>(numnodes);
+        IssmDouble* basevalue    = xNew<IssmDouble>(numnodes);
+        IssmDouble* initvalue    = xNew<IssmDouble>(numnodes);
         basalelement->GetDofList(&doflist,NoneApproximationEnum,GsetEnum);
 
@@ -422 +433 @@
         /*Use the dof list to index into the solution vector: */
         /*If the EPL is not active we revert to the bedrock elevation*/
+
         if(active_element){
                 for(int i=0;i<numnodes;i++){
@@ -430 +442 @@
         }
         else{
-                basalelement->GetInputListOnVertices(&eplHeads[0],BaseEnum);
+                //tradeof between keeping initial condition and not getting to far from head at deactivation
+                basalelement->GetInputListOnVertices(&basevalue[0],BaseEnum);
+                basalelement->GetInputListOnVertices(&initvalue[0],EplHeadEnum);
                 for(int i=0;i<numnodes;i++){
+                        eplHeads[i]=max(basevalue[i],initvalue[i]);
                         if(xIsNan<IssmDouble>(eplHeads[i])) _error_("NaN found in solution vector");
                         if(xIsInf<IssmDouble>(eplHeads[i])) _error_("Inf found in solution vector");
@@ -440 +455 @@
         /*Free ressources:*/
         xDelete<IssmDouble>(eplHeads);
+        xDelete<IssmDouble>(basevalue);
+        xDelete<IssmDouble>(initvalue);
         xDelete<int>(doflist);
         if(domaintype!=Domain2DhorizontalEnum){basalelement->DeleteMaterials(); delete basalelement;};
@@ -464 +481 @@
         case 1:
                 element->FindParam(&leakage,HydrologydcLeakageFactorEnum);
-                transfer=leakage;
+                transfer=+leakage;
                 break;
         default:
@@ -489 +506 @@
                 sed_head_input->GetInputValue(&sediment_head,gauss);
                 element->FindParam(&leakage,HydrologydcLeakageFactorEnum);
-                transfer=sediment_head*leakage;
+                transfer=+sediment_head*leakage;
                 break;
         default:
@@ -558 +575 @@
                 element->GetInputListOnVertices(&ice_thickness[0],ThicknessEnum);
                 element->GetInputListOnVertices(&bed[0],BaseEnum);
-                
+
                 if(!active_element){
                         /*Keeping thickness to initial value if EPL is not active*/
@@ -731 +748 @@
         base_input->GetInputValue(&base_elev,gauss);
         water_sheet=max(0.0,(prestep_head-base_elev));
-
         storing=rho_freshwater*g*epl_porosity*epl_thickness*(water_compressibility+(epl_compressibility/epl_porosity));
 
@@ -760 +776 @@
 
         water_sheet=max(0.0,(prestep_head-base_elev));
-        
         epl_transmitivity=epl_conductivity*epl_thickness;
         //epl_transmitivity=max(1.0e-6,(epl_conductivity*min(water_sheet,epl_thickness)));