Changeset 23387

Timestamp:

10/05/18 07:12:10 (6 years ago)

Author:

bdef

Message:

CHG: streamlining in EPL activation and dealing with frozen nodes EffectivPressure

Location:

issm/trunk-jpl/src/c

Files:

• analyses/HydrologyDCEfficientAnalysis.cpp (modified) (2 diffs)
• analyses/HydrologyDCInefficientAnalysis.cpp (modified) (5 diffs)
• solutionsequences/solutionsequence_hydro_nonlinear.cpp (modified) (2 diffs)

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

@@ -607 +607 @@
                                 /*centered*/
                                 //thickness[i] = old_thickness[i]*(1.0+opening-closing)/(1.0-opening+closing);
-
-                                //if(element->nodes[i]->Sid()==2299){
-                                //      printf("for node %i \n",element->nodes[i]->Sid());
-                                //      printf(" old thickness is  %g \n",old_thickness[i]);
-                                //      printf(" new thickness is  %g \n",thickness[i]);
-                                //}
                                 /*Take care of otherthikening*/
                                 if(thickness[i]>max_thick){
@@ -709 +703 @@
                 if(old_active[i]>0.){
                         vec_mask->SetValue(basalelement->nodes[i]->Sid(),1.,INS_VAL);
-                        /* If epl thickness gets under colapse thickness, close the layer */
-                        if(epl_thickness[i]<colapse_thick){
-                                //if(element->nodes[i]->Sid()==2299){
-                                        //printf("masking node %i with thickness %f\n",element->nodes[i]->Sid(),epl_thickness[i]);}
+                        /* If epl thickness gets under colapse thickness, close the layer if there is no residual*/
+                        if(epl_thickness[i]<colapse_thick && residual[i]<=0.){
                                 vec_mask->SetValue(basalelement->nodes[i]->Sid(),0.,INS_VAL);
                                 recurence->SetValue(basalelement->nodes[i]->Sid(),1.,INS_VAL);
                         }
                 }
-                /*If node is now closed bring its thickness back to initial*/
                 if (old_active[i]==0.){
-                        epl_thickness[i]=init_thick;
                         /*Activate if we have a residual from sediment*/
                         if(residual[i]>0.){
                                 vec_mask->SetValue(basalelement->nodes[i]->Sid(),1.,INS_VAL);
-                                if(old_active[i]==0.){
-                                        recurence->SetValue(basalelement->nodes[i]->Sid(),1.,INS_VAL);
-                                }
+                                recurence->SetValue(basalelement->nodes[i]->Sid(),1.,INS_VAL);
+                        }
+                        else{
+                                /*If node is now closed bring its thickness back to initial*/
+                                epl_thickness[i]=init_thick;
+                                vec_mask->SetValue(basalelement->nodes[i]->Sid(),0.,INS_VAL);
                         }
                 }

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

@@ -28 +28 @@
 
         /*retrieve some parameters: */
-         bool   issmb;
+        bool   issmb;
         iomodel->FindConstant(&issmb,"md.transient.issmb");
         iomodel->FindConstant(&hydrology_model,"md.hydrology.model");
@@ -374 +374 @@
         }
         if(smb_model==SMBgradientscomponentsEnum){
-                surface_runoff_input = basalelement->GetInput(SmbRunoffEnum); _assert_(surface_runoff_input); 
+                surface_runoff_input = basalelement->GetInput(SmbRunoffEnum); _assert_(surface_runoff_input);
         }
 
@@ -518 +518 @@
         /*Use the dof list to index into the solution vector: */
 
-        /*need to introduce thawed test herre*/
         for(int i=0;i<numnodes;i++){
                 values[i] =solution[doflist[i]];
@@ -546 +545 @@
                 kappa=kmax*pow(10.,penalty_factor);
 
+
+                Input* thawed_element_input = basalelement->GetInput(HydrologydcMaskThawedEltEnum); _assert_(thawed_element_input);
+                thawed_element_input->GetInputValue(&thawed_element);
+
                 for(int i=0;i<numnodes;i++){
+                        /*frozen elements heads are set to base elevation*/
+                        if(!thawed_element){
+                                values[i]=base[i];
+                        }
                         GetHydrologyDCInefficientHmax(&h_max,basalelement,basalelement->GetNode(i));
                         if(values[i]>h_max) {
@@ -555 +562 @@
                                 residual[i] = 0.;
                         }
-                        //adding base in min to take into account heads under bed wich don't change N
-                        //pressure[i]=(rho_ice*g*thickness[i])-(rho_freshwater*g*(max((min(h_max,values[i])-base[i]),0.0)));
                         pressure[i]=(rho_ice*g*thickness[i])-(rho_freshwater*g*(values[i]-base[i]));
+
                 }
                 xDelete<IssmDouble>(thickness);

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

@@ -190 +190 @@
                                 Reduceloadx(pf,Kfs,ys); delete Kfs;
                                 delete uf_epl;
+
                                 femmodel->profiler->Start(SOLVER);
                                 Solverx(&uf_epl,Kff,pf,uf_epl_sub_iter,df,femmodel->parameters);
                                 femmodel->profiler->Stop(SOLVER);
+
                                 delete Kff; delete pf; delete df;
                                 delete uf_epl_sub_iter;
-                                uf_epl_sub_iter=uf_epl->Duplicate();
+                                uf_epl_sub_iter=uf_epl->Duplicate();_assert_(uf_epl_sub_iter);
                                 uf_epl->Copy(uf_epl_sub_iter);
                                 delete ug_epl;
@@ -215 +217 @@
                                         _error_("   maximum number of EPL iterations (" << hydro_maxiter << ") exceeded");
                                 }
-                                //If there is some colapse go through sediment again
-                                /* if(ThickCount<L2Count)eplconverged=true; */
                                 eplcount++;