Changeset 17006

Timestamp:

12/04/13 16:59:05 (11 years ago)

Author:

bdef

Message:

BUG: modification in the way the transfer deals with its mask

Location:

issm/trunk-jpl/src/c

Files:

• classes/Elements/Tria.cpp (modified) (5 diffs)
• classes/Materials/Matpar.cpp (modified) (1 diff)
• solutionsequences/solutionsequence_hydro_nonlinear.cpp (modified) (15 diffs)

issm/trunk-jpl/src/c/classes/Elements/Tria.cpp

@@ -4659 +4659 @@
         bool       isefficientlayer;
         bool       active_element;
-        int        transfermethod;
+        int        transfermethod,step;
         IssmDouble sed_trans,sed_thick;
         IssmDouble leakage,h_max;
@@ -4666 +4666 @@
         IssmDouble epl_specificstoring[numdof],sedstoring[numdof];
         IssmDouble epl_head[numdof],sed_head[numdof];
+        IssmDouble old_transfer[numdof];
 
         Input* active_element_input=NULL;
 
+                
         GetDofList(&doflist,NoneApproximationEnum,GsetEnum);
 
         /*Get the flag to know if the efficient layer is present*/
         this->parameters->FindParam(&isefficientlayer,HydrologydcIsefficientlayerEnum);
+
+        this->parameters->FindParam(&step,StepEnum);
 
         if(isefficientlayer){
@@ -4685 +4689 @@
                 case 1:
 
+        
                         active_element_input=inputs->GetInput(HydrologydcMaskEplactiveEltEnum); _assert_(active_element_input);
                         active_element_input->GetInputValue(&active_element);
@@ -4698 +4703 @@
 
                         if(!active_element){
-
                                 /*No transfer if the EPL is not active*/
-                                for(int i=0;i<numdof;i++){
-                                        wh_trans=0.0;
-                                        /*Assign output pointer*/
-                                        transfer->SetValue(doflist[i],wh_trans,INS_VAL);
-                                }
+                                /* for(int i=0;i<numdof;i++){ */
+                                /*      wh_trans=0.0; */
+                                /*      /\*Assign output pointer*\/ */
+                                /*      transfer->SetValue(doflist[i],wh_trans,INS_VAL); */
+                                /* } */
                         }
                         else{
-
+                        
+                                //GetInputListOnVertices(&old_transfer[0],WaterTransferEnum);
+                        
                                 for(int i=0;i<numdof;i++){
                                         epl_specificstoring[i]=matpar->GetEplSpecificStoring();         
                                         sedstoring[i]=matpar->GetSedimentStoring();
-                                        
+                                        this->GetHydrologyDCInefficientHmax(&h_max,i);
+                                                
                                         /*EPL head higher than sediment head, transfer from the epl to the sediment*/
                                         if(epl_head[i]>sed_head[i]){
                                                 wh_trans=epl_specificstoring[i]*epl_thickness[i]*sed_trans*(epl_head[i]-sed_head[i])/(leakage*sed_thick);                               
-                                                
+
                                                 /*No transfer if the sediment head is allready at the maximum*/
-                                                this->GetHydrologyDCInefficientHmax(&h_max,i);
-                                                if(sed_head[i]>=h_max)wh_trans=0.0;
+                                                if(sed_head[i]>=h_max){
+                                                        wh_trans=0.0;
+                                                }
                                         }
                                         /*EPL head lower than sediment head, transfer from the sediment to the epl*/
                                         else if(epl_head[i]<=sed_head[i]){
-                                                wh_trans=sedstoring[i]*sed_trans*(epl_head[i]-sed_head[i])/(leakage*sed_thick);                         
+                                                wh_trans=sedstoring[i]*sed_trans*(epl_head[i]-sed_head[i])/(leakage*sed_thick);
                                         }
+
+                                        /*Introduce relaxation*/
+                                        //wh_trans=old_transfer[i]+0.9*(wh_trans-old_transfer[i]);
+                                
                                         /*Assign output pointer*/
                                         transfer->SetValue(doflist[i],wh_trans,INS_VAL);
-                                        /* if(nodes[i]->id>=54){ */
-                                        /*      printf("%i %e %e %e \n",nodes[i]->id-54,wh_trans,sed_head[i],epl_head[i]); */
-                                        /* } */
-                                        /* else{*/
-                                        /*      printf("%i %e %e %e \n",nodes[i]->id,wh_trans,sed_head[i],epl_head[i]); */
                                 }
                         }
@@ -4799 +4806 @@
                 }
                 /*If epl thickness gets under 0, close the layer*/
-                /* else if(epl_thickness[i]<0.0){ */
-                /*      vec_mask->SetValue(nodes[i]->Sid(),0.,INS_VAL); */
-                /* } */
+                else if(epl_thickness[i]<0.0){
+                        vec_mask->SetValue(nodes[i]->Sid(),0.,INS_VAL);
+                }
                 /*Increase of the efficient system is needed if the epl head reach the maximum value (sediment max value for now)*/
                 this->GetHydrologyDCInefficientHmax(&h_max,i);

issm/trunk-jpl/src/c/classes/Materials/Matpar.cpp

@@ -382 +382 @@
 /*FUNCTION Matpar::GetSedimentStoring {{{*/
 IssmDouble Matpar::GetSedimentStoring(){
-        return this->rho_freshwater* this->g* this->sediment_porosity* this->sediment_thickness*
-    ( this->water_compressibility+( this->sediment_compressibility/ this->sediment_porosity));           
+        return this->rho_freshwater * this->g * this->sediment_porosity * this->sediment_thickness *
+    (this->water_compressibility + (this->sediment_compressibility / this->sediment_porosity));          
 }                
 /*}}}*/ 
 /*FUNCTION Matpar::GetEplSpecificStoring {{{*/
 IssmDouble Matpar::GetEplSpecificStoring(){
-        return this->rho_freshwater* this->g* this->epl_porosity* 
-    ( this->water_compressibility+( this->epl_compressibility/ this->epl_porosity));             
+        return this->rho_freshwater * this->g * this->epl_porosity * 
+    (this->water_compressibility + (this->epl_compressibility / this->epl_porosity));            
 }                
 /*}}}*/ 

issm/trunk-jpl/src/c/solutionsequences/solutionsequence_hydro_nonlinear.cpp

@@ -80 +80 @@
                 ug_sed->Copy(ug_sed_main_iter);
                 
-                //test
-                /* uf_sed_sub_iter=uf_sed->Duplicate(); */
-                /* uf_sed->Copy(uf_sed_sub_iter); */
-
                 if(isefficientlayer){
                         ug_epl_main_iter=ug_epl->Duplicate();
                         ug_epl->Copy(ug_epl_main_iter);
-                        //test
-                        ug_epl_sub_iter=ug_epl->Duplicate();
-                        ug_epl->Copy(ug_epl_sub_iter);
                 }
 
@@ -104 +97 @@
                 sedconverged=false;
                 for(;;){
+
                         uf_sed_sub_iter=uf_sed->Duplicate();_assert_(uf_sed_sub_iter);
                         uf_sed->Copy(uf_sed_sub_iter);
@@ -121 +115 @@
                                 if(num_unstable_constraints==0) sedconverged = true;
                                 if (sedcount>=hydro_maxiter){
+                                        //sedconverged=true;
                                         _error_("   maximum number of Sediment iterations (" << hydro_maxiter << ") exceeded");
                                 }
@@ -128 +123 @@
 
                         if(sedconverged){
-                                sedconverged=false;
+                                if(isefficientlayer){ 
+                                        /*Updating Elemental Mask*/
+                                        HydrologyDCInefficientAnalysis* analysis = new HydrologyDCInefficientAnalysis();
+                                        analysis->ElementizeEplMask(femmodel);
+                                        delete analysis;
+                                        femmodel->HydrologyTransferx();
+                                }
+                                sedconverged=false;     
+                                /*Checking convegence on the value of the sediment head*/
                                 duf=uf_sed_sub_iter->Duplicate();_assert_(duf);
                                 uf_sed_sub_iter->Copy(duf);
@@ -137 +140 @@
                                 if (xIsNan<IssmDouble>(ndu_sed) || xIsNan<IssmDouble>(nu_sed)) _error_("convergence criterion is NaN!");
                                 if (ndu_sed==0.0 && nu_sed==0.0) nu_sed=1.0e-6; /*Hacking the case where the EPL is used but empty*/
+                                if(VerboseConvergence()) _printf0_(setw(50) << left << "   Inner Sediment Convergence criterion:" << ndu_sed/nu_sed*100 << " aiming lower than " << eps_hyd*100 << " %\n");
                                 if((ndu_sed/nu_sed)<eps_hyd){
                                         if(VerboseConvergence()) _printf0_("   # Inner sediment convergence achieve \n");
@@ -145 +149 @@
 
                         if(sedconverged){
-                                if(isefficientlayer){
-                                        /*Updating Nodal Mask*/
-                                        HydrologyDCInefficientAnalysis* analysis = new HydrologyDCInefficientAnalysis();
-                                        analysis->ElementizeEplMask(femmodel);
-                                        delete analysis;
-                                        femmodel->HydrologyTransferx();
-                                }
                                 femmodel->parameters->SetParam(sediment_kmax,HydrologySedimentKmaxEnum);
                                 InputUpdateFromConstantx(femmodel,sedconverged,ConvergedEnum);
@@ -166 +163 @@
                         InputUpdateFromConstantx(femmodel,true,ResetPenaltiesEnum);
                         InputUpdateFromConstantx(femmodel,false,ConvergedEnum);
-                        femmodel->HydrologyEPLupdateDomainx();
                         femmodel->parameters->SetParam(HydrologyEfficientEnum,HydrologyLayerEnum);
 
@@ -188 +184 @@
                                 //updating mask after the computation of the epl thickness (Allow to close too thin EPL)
                                 femmodel->HydrologyEPLupdateDomainx();
+                                        
                                 ug_epl_sub_iter=ug_epl->Duplicate();_assert_(ug_epl_sub_iter);
-                                ug_epl->Copy(ug_epl_sub_iter);_assert_(ug_epl);
-                                
+                                ug_epl->Copy(ug_epl_sub_iter);
                                 SystemMatricesx(&Kff,&Kfs,&pf,&df,NULL,femmodel);
                                 CreateNodalConstraintsx(&ys,femmodel->nodes,HydrologyDCEfficientAnalysisEnum);
@@ -204 +200 @@
                                 InputUpdateFromSolutionx(femmodel,ug_epl);
                                 ConstraintsStatex(&constraints_converged,&num_unstable_constraints,femmodel);
-                                femmodel->HydrologyEPLupdateDomainx();                  
+
 
                                 if (!eplconverged){
@@ -210 +206 @@
                                         if(num_unstable_constraints==0) eplconverged = true;
                                         if (eplcount>=hydro_maxiter){
-                                        /*Hacking to get the results of non converged runs*/
-                                        //eplconverged = true;
+                                                //eplconverged =true;
                                                 _error_("   maximum number of EPL iterations (" << hydro_maxiter << ") exceeded");
                                         }
@@ -220 +215 @@
                                 if(eplconverged){
                                         eplconverged=false;
+                                
                                         dug=ug_epl_sub_iter->Duplicate();_assert_(dug);
                                         ug_epl_sub_iter->Copy(dug);
@@ -229 +225 @@
                                         if (xIsNan<IssmDouble>(ndu_epl) || xIsNan<IssmDouble>(nu_epl)) _error_("convergence criterion is NaN!");
                                         if (ndu_epl==0.0 && nu_epl==0.0) nu_epl=1.0e-6; /*Hacking the case where the EPL is used but empty*/
+                                        if(VerboseConvergence()) _printf0_(setw(50) << left << "   Inner EPL Convergence criterion:" << ndu_epl/nu_epl*100 << " aiming lower than " << eps_hyd*100 << " %\n");
                                         if((ndu_epl/nu_epl)<eps_hyd)eplconverged=true;
                                 }
@@ -234 +231 @@
 
                                 if(eplconverged){
-
-                                        /*Updating Nodal Mask*/
+                                        /*Updating Elemental Mask*/
                                         HydrologyDCInefficientAnalysis* analysis = new HydrologyDCInefficientAnalysis();
                                         analysis->ElementizeEplMask(femmodel);
                                         delete analysis;
                                         femmodel->HydrologyTransferx();
-
+                                        
                                         InputUpdateFromConstantx(femmodel,eplconverged,ConvergedEnum);
                                         InputUpdateFromConstantx(femmodel,sediment_kmax,MeltingOffsetEnum);
@@ -248 +244 @@
                         }
                 }
-                
                 
                 /*System convergence check*/
@@ -303 +298 @@
                 hydrocount++;
                 if(hydroconverged)break;
-        }
+}
         
         InputUpdateFromSolutionx(femmodel,ug_sed);