Changeset 17030

Timestamp:

12/19/13 08:25:02 (11 years ago)

Author:

bdef

Message:

CHG: uniformisation of the mask gestion in Hydrology

Location:

issm/trunk-jpl/src/c

Files:

• analyses/HydrologyDCEfficientAnalysis.cpp (modified) (7 diffs)
• analyses/HydrologyDCInefficientAnalysis.cpp (modified) (3 diffs)
• analyses/L2ProjectionEPLAnalysis.cpp (modified) (5 diffs)
• classes/Elements/Tria.cpp (modified) (7 diffs)
• solutionsequences/solutionsequence_hydro_nonlinear.cpp (modified) (4 diffs)

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

@@ -58 +58 @@
         iomodel->FetchDataToInput(elements,EplHeadEnum);
         iomodel->FetchDataToInput(elements,HydrologydcEplInitialThicknessEnum);
+
+        //      iomodel->FetchDataToInput(elements,HydrologydcMaskEplactiveNodeEnum);
         
         elements->InputDuplicate(HydrologydcEplInitialThicknessEnum,HydrologydcEplThicknessEnum);
@@ -98 +100 @@
 
 /*Finite Element Analysis*/
-void           HydrologyDCEfficientAnalysis::Core(FemModel* femmodel){/*{{{*/
+void HydrologyDCEfficientAnalysis::Core(FemModel* femmodel){/*{{{*/
         _error_("not implemented");
 }/*}}}*/
+
 ElementVector* HydrologyDCEfficientAnalysis::CreateDVector(Element* element){/*{{{*/
         /*Default, return NULL*/
         return NULL;
 }/*}}}*/
+
 ElementMatrix* HydrologyDCEfficientAnalysis::CreateJacobianMatrix(Element* element){/*{{{*/
 _error_("Not implemented");
 }/*}}}*/
+
 ElementMatrix* HydrologyDCEfficientAnalysis::CreateKMatrix(Element* element){/*{{{*/
 
         /*Intermediaries*/
         int      meshtype;
+        bool     active_element;
         Element* basalelement;
+        Input*   active_element_input=NULL;
 
         /*Get basal element*/
@@ -127 +134 @@
         }
 
+        /* active_element_input=basalelement->GetInput(HydrologydcMaskEplactiveEltEnum); _assert_(active_element_input); */
+        /* active_element_input->GetInputValue(&active_element); */
+
         /*Check that all nodes are active, else return empty matrix*/
         if(!basalelement->AllActive()) {
-                if(meshtype!=Mesh2DhorizontalEnum){
+                //if(!active_element){
+        if(meshtype!=Mesh2DhorizontalEnum){
                         basalelement->DeleteMaterials(); 
                         delete basalelement;
@@ -199 +210 @@
         /*Intermediaries*/
         int      meshtype;
+        bool       active_element;
         Element* basalelement;
+        Input*   active_element_input=NULL;
+
 
         /*Get basal element*/
@@ -213 +227 @@
                 default: _error_("mesh "<<EnumToStringx(meshtype)<<" not supported yet");
         }
+        /* active_element_input=basalelement->GetInput(HydrologydcMaskEplactiveEltEnum); _assert_(active_element_input); */
+        /* active_element_input->GetInputValue(&active_element); */
 
         /*Check that all nodes are active, else return empty matrix*/
-        if(!basalelement->AllActive()){
-                if(meshtype!=Mesh2DhorizontalEnum){
+        if(!basalelement->AllActive()) {
+                //if(!active_element){
+        if(meshtype!=Mesh2DhorizontalEnum){
                         basalelement->DeleteMaterials(); 
                         delete basalelement;
@@ -289 +306 @@
         return pe;
 }/*}}}*/
+
 void HydrologyDCEfficientAnalysis::GetSolutionFromInputs(Vector<IssmDouble>* solution,Element* element){/*{{{*/
         element->GetSolutionFromInputsOneDof(solution,EplHeadEnum);
 }/*}}}*/
+
 void HydrologyDCEfficientAnalysis::InputUpdateFromSolution(IssmDouble* solution,Element* element){/*{{{*/
         int meshtype;
@@ -315 +334 @@
         return rho_freshwater*g*epl_porosity*(water_compressibility+(epl_compressibility/epl_porosity));                 
 }/*}}}*/
+
 void HydrologyDCEfficientAnalysis::GetB(IssmDouble* B,Element* element,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/
         /*Compute B  matrix. B=[B1 B2 B3] where Bi is of size 3*NDOF2. 

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

@@ -296 +296 @@
         return pe;
 }/*}}}*/
+
 void HydrologyDCInefficientAnalysis::GetB(IssmDouble* B,Element* element,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/
         /*Compute B  matrix. B=[B1 B2 B3] where Bi is of size 3*NDOF2. 
@@ -323 +324 @@
         xDelete<IssmDouble>(dbasis);
 }/*}}}*/
+
 void HydrologyDCInefficientAnalysis::GetSolutionFromInputs(Vector<IssmDouble>* solution,Element* element){/*{{{*/
         element->GetSolutionFromInputsOneDof(solution,SedimentHeadEnum);
 }/*}}}*/
+
 void HydrologyDCInefficientAnalysis::InputUpdateFromSolution(IssmDouble* solution,Element* element){/*{{{*/
 
@@ -395 +398 @@
         return rho_freshwater*g*sediment_porosity*sediment_thickness*(water_compressibility+(sediment_compressibility/sediment_porosity));               
 }/*}}}*/
+
 void HydrologyDCInefficientAnalysis::ElementizeEplMask(FemModel* femmodel){/*{{{*/
 

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

@@ -73 +73 @@
         _error_("not implemented");
 }/*}}}*/
+
 ElementVector* L2ProjectionEPLAnalysis::CreateDVector(Element* element){/*{{{*/
         /*Default, return NULL*/
         return NULL;
 }/*}}}*/
+
 ElementMatrix* L2ProjectionEPLAnalysis::CreateJacobianMatrix(Element* element){/*{{{*/
 _error_("Not implemented");
 }/*}}}*/
+
 ElementMatrix* L2ProjectionEPLAnalysis::CreateKMatrix(Element* element){/*{{{*/
 
         /*Intermediaries*/
         int      meshtype;
+        bool     active_element;
         Element* basalelement;
+        Input*   active_element_input=NULL;
 
         /*Get basal element*/
@@ -101 +106 @@
                         break;
                 default: _error_("mesh "<<EnumToStringx(meshtype)<<" not supported yet");
+        }
+
+        /* active_element_input=basalelement->GetInput(HydrologydcMaskEplactiveEltEnum); _assert_(active_element_input); */
+        /* active_element_input->GetInputValue(&active_element); */
+
+        /* Check that all nodes are active, else return empty matrix */
+        if(!basalelement->AllActive()) {
+                //      if(!active_element){
+        if(meshtype!=Mesh2DhorizontalEnum){
+                        basalelement->DeleteMaterials();
+                        delete basalelement;
+                }
+                return NULL;
         }
 
@@ -143 +161 @@
         /*Intermediaries*/
         int      meshtype;
+        bool     active_element;
         Element* basalelement;
+        Input*   active_element_input=NULL;
 
         /*Get basal element*/
@@ -158 +178 @@
         }
 
+        /* active_element_input=basalelement->GetInput(HydrologydcMaskEplactiveEltEnum); _assert_(active_element_input); */
+        /* active_element_input->GetInputValue(&active_element); */
+
+        /*Check that all nodes are active, else return empty matrix*/
+        if(!basalelement->AllActive()) {
+                /* if(!active_element){ */
+                if(meshtype!=Mesh2DhorizontalEnum){
+                        basalelement->DeleteMaterials();
+                        delete basalelement;
+                }
+                return NULL;
+        }
+        
         /*Intermediaries */
         int         input_enum,index;
@@ -199 +232 @@
         return pe;
 }/*}}}*/
+
 void L2ProjectionEPLAnalysis::GetSolutionFromInputs(Vector<IssmDouble>* solution,Element* element){/*{{{*/
            _error_("not implemented yet");
 }/*}}}*/
+
 void L2ProjectionEPLAnalysis::InputUpdateFromSolution(IssmDouble* solution,Element* element){/*{{{*/
         int inputenum,meshtype;

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

@@ -4729 +4729 @@
 
                                         /*Introduce relaxation*/
-                                        wh_trans=old_transfer[i]+0.66*(wh_trans-old_transfer[i]);
+                                        //                                      wh_trans=old_transfer[i]+0.66*(wh_trans-old_transfer[i]);
+
+                                        if(this->Id()==27){
+                                                printf("old %e, transfer is %e ,eplhead %e, sedhead %e, factor %e \n",old_transfer[i],wh_trans,epl_head[i],sed_head[i],epl_specificstoring[i]*epl_thickness[i]*sed_trans[i]);
+                                        }
                                 
                                         /*Assign output pointer*/
@@ -4782 +4786 @@
         IssmDouble  residual[numdof];
 
+        bool       active_element;
+        Input* active_element_input=NULL;
+
+        active_element_input=inputs->GetInput(HydrologydcMaskEplactiveEltEnum); _assert_(active_element_input);
+        active_element_input->GetInputValue(&active_element);
+
         GetInputListOnVertices(&old_active[0],HydrologydcMaskEplactiveNodeEnum);        
         GetInputListOnVertices(&epl_thickness[0],HydrologydcEplThicknessEnum);  
@@ -4808 +4818 @@
                 /*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);
-                if(eplhead[i]>=h_max && this->AnyActive()){
+                //if(eplhead[i]>=h_max && this->AnyActive()){
+                if(eplhead[i]>=h_max && active_element){
                         for(j=0;j<numdof;j++){
                                 if(old_active[j]>0.){
@@ -4838 +4849 @@
         IssmDouble  ice_thickness[numdof],bed[numdof];
 
+        bool       active_element;
+        Input* active_element_input=NULL;
 
         /*Get the flag to know if the efficient layer is present*/
@@ -4854 +4867 @@
                 A                = material->GetAbar();
                 
-                GetInputListOnVertices(&activeEpl[0],HydrologydcMaskEplactiveNodeEnum);
+                //      GetInputListOnVertices(&activeEpl[0],HydrologydcMaskEplactiveNodeEnum);
+
+                active_element_input=inputs->GetInput(HydrologydcMaskEplactiveEltEnum); _assert_(active_element_input);
+                active_element_input->GetInputValue(&active_element);
+                        
                 GetInputListOnVertices(&eplhead[0],EplHeadEnum);
                 GetInputListOnVertices(&epl_slopeX[0],EplHeadSlopeXEnum); 
@@ -4862 +4879 @@
                 GetInputListOnVertices(&bed[0],BedEnum);
                 
-                if(!this->AnyActive()){
+                //if(!this->AnyActive()){
+                if(!active_element){
                         /*Keeping thickness to initial value if EPL is not active*/
                         for(int i=0;i<numdof;i++){ 
@@ -4878 +4896 @@
                                         /*And proceed to the real thing*/
                                         thickness[i] = old_thickness[i]*(1+((rho_water*gravity*dt)/(rho_ice*latentheat))*epl_conductivity*EPLgrad2-2.0*(A*dt/(pow(n,n)))*(pow(EPL_N,n)));
+                                        //                                      thickness[i] = old_thickness[i]+0.66*(thickness[i]-old_thickness[i]);
+                                        if(this->Id()==27){
+                                                printf("old %e, thickness is %e ,eplhead %e, slopeX %e, slopeY %e\n",old_thickness[i],thickness[i],eplhead[i],epl_slopeX[i],epl_slopeY[i]);
+                                        }
                         }
                 }

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

@@ -169 +169 @@
                         femmodel->parameters->SetParam(HydrologyEfficientEnum,HydrologyLayerEnum);
                 
-                        eplconverged = false;
-                        
+                        eplconverged=false;
                         for(;;){
                                 ug_epl_sub_iter=ug_epl->Duplicate();_assert_(ug_epl_sub_iter);
                                 ug_epl->Copy(ug_epl_sub_iter);
-                                        
+                                
+                                /*Loop on the epl layer to deal with the penalization, have been removed TO CHECK*/
+                                                                                                
+                                SystemMatricesx(&Kff,&Kfs,&pf,&df,NULL,femmodel);
+                                CreateNodalConstraintsx(&ys,femmodel->nodes,HydrologyDCEfficientAnalysisEnum);
+                                Reduceloadx(pf,Kfs,ys); delete Kfs;
+                                delete uf_epl;
+                                Solverx(&uf_epl,Kff,pf,uf_epl_sub_iter,df,femmodel->parameters);
+                                delete Kff; delete pf; delete df;
+                                delete uf_epl_sub_iter;
+                                uf_epl_sub_iter=uf_epl->Duplicate();
+                                uf_epl->Copy(uf_epl_sub_iter);
+                                delete ug_epl; 
+                                Mergesolutionfromftogx(&ug_epl,uf_epl,ys,femmodel->nodes,femmodel->parameters); delete ys;
+                                InputUpdateFromSolutionx(femmodel,ug_epl);
+                                ConstraintsStatex(&constraints_converged,&num_unstable_constraints,femmodel);
+                                
                                 /* {{{ *//*Start by retrieving the EPL head slopes and compute EPL Thickness*/
                                 if(VerboseSolution()) _printf0_("computing EPL Head slope...\n");
@@ -183 +198 @@
                                 femmodel->parameters->SetParam(EplHeadSlopeYEnum,InputToL2ProjectEnum);
                                 solutionsequence_linear(femmodel);
+                                femmodel->SetCurrentConfiguration(HydrologyDCEfficientAnalysisEnum);
                                 femmodel->HydrologyEPLThicknessx();
+                                
                                 //updating mask after the computation of the epl thickness (Allow to close too thin EPL)
                                 femmodel->HydrologyEPLupdateDomainx();
                                 /* }}} */
-                                femmodel->SetCurrentConfiguration(HydrologyDCEfficientAnalysisEnum);
-                                
-                                /* {{{ *//*Loop on the epl layer to deal with the penalization*/
-                                for(;;){
-                                        SystemMatricesx(&Kff,&Kfs,&pf,&df,NULL,femmodel);
-                                        CreateNodalConstraintsx(&ys,femmodel->nodes,HydrologyDCEfficientAnalysisEnum);
-                                        Reduceloadx(pf,Kfs,ys); delete Kfs;
-                                        delete uf_epl;
-                                        Solverx(&uf_epl,Kff,pf,uf_epl_sub_iter,df,femmodel->parameters);
-                                        delete Kff; delete pf; delete df;
-                                        delete uf_epl_sub_iter;
-                                        uf_epl_sub_iter=uf_epl->Duplicate();
-                                        uf_epl->Copy(uf_epl_sub_iter);
-                                        delete ug_epl; 
-                                        Mergesolutionfromftogx(&ug_epl,uf_epl,ys,femmodel->nodes,femmodel->parameters); delete ys;
-                                        InputUpdateFromSolutionx(femmodel,ug_epl);
-                                        ConstraintsStatex(&constraints_converged,&num_unstable_constraints,femmodel);
-                                
-                                        if (!eplconverged){
-                                                if(VerboseConvergence()) _printf0_("   # EPL unstable constraints = " << num_unstable_constraints << "\n");
-                                                if(num_unstable_constraints==0) eplconverged = true;
-                                                if (eplcount>=hydro_maxiter){
-                                                        _error_("   maximum number of EPL iterations (" << hydro_maxiter << ") exceeded");
-                                                }
-                                        }
-                                        /*Add an iteration and get out of the loop if the penalisation is converged*/
-                                        eplcount++;
-                                        if(eplconverged) break;
-                                }
-                                /* }}} */ /*End of the EPL penalization loop*/
-
-                                eplconverged=false;
+                                
                                 /*Update Elemental Mask and compute transfer*/
                                 HydrologyDCInefficientAnalysis* analysis = new HydrologyDCInefficientAnalysis();
@@ -224 +210 @@
                                 delete analysis;
                                 femmodel->HydrologyTransferx();
-
+                                
                                 dug=ug_epl_sub_iter->Duplicate();_assert_(dug);
                                 ug_epl_sub_iter->Copy(dug);
@@ -235 +221 @@
                                 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;
-                        
+                                if (eplcount>=hydro_maxiter){
+                                        _error_("   maximum number of EPL iterations (" << hydro_maxiter << ") exceeded");
+                                }
+                                eplcount++;
+                                
                                 delete ug_epl_sub_iter;
                                 if(eplconverged){