Changeset 15161

Timestamp:

05/30/13 10:08:07 (12 years ago)

Author:

bdef

Message:

BUG: fixing the convergence issue in the hydrological solver

Location:

issm/trunk-jpl/src/c

Files:

• analyses/hydrology_core.cpp (modified) (1 diff)
• classes/Elements/Tria.cpp (modified) (5 diffs)
• modules/ModelProcessorx/HydrologyDCInefficient/UpdateElementsHydrologyDCInefficient.cpp (modified) (1 diff)
• shared/Enum/EnumDefinitions.h (modified) (1 diff)
• shared/Enum/EnumToStringx.cpp (modified) (1 diff)
• shared/Enum/StringToEnumx.cpp (modified) (5 diffs)
• solutionsequences/solutionsequence_hydro_nonlinear.cpp (modified) (2 diffs)

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

@@ -81 +81 @@
                         InputDuplicatex(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,SedimentHeadEnum,SedimentHeadOldEnum);
                         femmodel->parameters->FindParam(&isefficientlayer,HydrologydcIsefficientlayerEnum);
-
-                        if(VerboseSolution()) _printf0_("   computing water transfer\n");
+                        if (isefficientlayer){
+                                InputDuplicatex(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,EplHeadEnum,EplHeadOldEnum);
+                        }
 
                         if(VerboseSolution()) _printf0_("   computing water head\n");

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

@@ -2135 +2135 @@
                                 name==SedimentHeadOldEnum ||
                                 name==SedimentHeadEnum ||
+                                name==EplHeadOldEnum ||
+                                name==EplHeadEnum ||
                                 name==WaterTransferEnum ||
                                 name==BasisIntegralEnum ||
@@ -6218 +6220 @@
 
         if(reCast<bool,IssmDouble>(dt)){
-                old_wh_input=inputs->GetInput(EplHeadEnum); _assert_(old_wh_input);
+                old_wh_input=inputs->GetInput(EplHeadOldEnum); _assert_(old_wh_input);
         }
 
@@ -6407 +6409 @@
                 this->parameters->FindParam(&penalty_factor,HydrologydcPenaltyFactorEnum);
                 kappa=kmax*pow(10.,penalty_factor);
-
+                
                 for(int i=0;i<NUMVERTICES;i++){
                         this->GetHydrologyDCInefficientHmax(&h_max,nodes[i]);
@@ -6417 +6419 @@
                                 residual[i]=0.0;
                 }
+                if(this->id==1)printf("residual, %g \n", residual[1]);
+                if(this->id==1)printf("hSed, %g \n", values[1]);
         }
         /*Add input to the element: */
@@ -6422 +6426 @@
         this->inputs->AddInput(new TriaP1Input(SedimentHeadEnum,values));
         this->inputs->AddInput(new TriaP1Input(SedimentHeadResidualEnum,residual));
-        if(converged)this->inputs->AddInput(new TriaP1Input(SedimentHeadOldEnum,values));
 
         /*Free ressources:*/

issm/trunk-jpl/src/c/modules/ModelProcessorx/HydrologyDCInefficient/UpdateElementsHydrologyDCInefficient.cpp

@@ -43 +43 @@
         iomodel->FetchDataToInput(elements,BasalforcingsMeltingRateEnum);
         iomodel->FetchDataToInput(elements,SedimentHeadEnum);
-        /*      iomodel->FetchDataToInput(elements,SedimentHeadOldEnum);*/
 
         /*Free data: */

issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h

@@ -90 +90 @@
         SedimentHeadResidualEnum,
         EplHeadEnum,
+        EplHeadOldEnum,
   HydrologydcRelTolEnum,
         HydrologydcSpcsedimentHeadEnum,

issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp

@@ -98 +98 @@
                 case SedimentHeadResidualEnum : return "SedimentHeadResidual";
                 case EplHeadEnum : return "EplHead";
+                case EplHeadOldEnum : return "EplHeadOld";
                 case HydrologydcRelTolEnum : return "HydrologydcRelTol";
                 case HydrologydcSpcsedimentHeadEnum : return "HydrologydcSpcsedimentHead";

issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp

@@ -98 +98 @@
               else if (strcmp(name,"SedimentHeadResidual")==0) return SedimentHeadResidualEnum;
               else if (strcmp(name,"EplHead")==0) return EplHeadEnum;
+              else if (strcmp(name,"EplHeadOld")==0) return EplHeadOldEnum;
               else if (strcmp(name,"HydrologydcRelTol")==0) return HydrologydcRelTolEnum;
               else if (strcmp(name,"HydrologydcSpcsedimentHead")==0) return HydrologydcSpcsedimentHeadEnum;
@@ -136 +137 @@
               else if (strcmp(name,"InversionMinParameters")==0) return InversionMinParametersEnum;
               else if (strcmp(name,"InversionNsteps")==0) return InversionNstepsEnum;
-              else if (strcmp(name,"InversionNumControlParameters")==0) return InversionNumControlParametersEnum;
          else stage=2;
    }
    if(stage==2){
-              if (strcmp(name,"InversionNumCostFunctions")==0) return InversionNumCostFunctionsEnum;
+              if (strcmp(name,"InversionNumControlParameters")==0) return InversionNumControlParametersEnum;
+              else if (strcmp(name,"InversionNumCostFunctions")==0) return InversionNumCostFunctionsEnum;
               else if (strcmp(name,"InversionStepThreshold")==0) return InversionStepThresholdEnum;
               else if (strcmp(name,"InversionThicknessObs")==0) return InversionThicknessObsEnum;
@@ -259 +260 @@
               else if (strcmp(name,"TimesteppingStartTime")==0) return TimesteppingStartTimeEnum;
               else if (strcmp(name,"TimesteppingFinalTime")==0) return TimesteppingFinalTimeEnum;
-              else if (strcmp(name,"TimesteppingCflCoefficient")==0) return TimesteppingCflCoefficientEnum;
          else stage=3;
    }
    if(stage==3){
-              if (strcmp(name,"TimesteppingTimeAdapt")==0) return TimesteppingTimeAdaptEnum;
+              if (strcmp(name,"TimesteppingCflCoefficient")==0) return TimesteppingCflCoefficientEnum;
+              else if (strcmp(name,"TimesteppingTimeAdapt")==0) return TimesteppingTimeAdaptEnum;
               else if (strcmp(name,"TimesteppingTimeStep")==0) return TimesteppingTimeStepEnum;
               else if (strcmp(name,"TransientIsdiagnostic")==0) return TransientIsdiagnosticEnum;
@@ -382 +383 @@
               else if (strcmp(name,"SpcDynamic")==0) return SpcDynamicEnum;
               else if (strcmp(name,"SpcStatic")==0) return SpcStaticEnum;
-              else if (strcmp(name,"SpcTransient")==0) return SpcTransientEnum;
          else stage=4;
    }
    if(stage==4){
-              if (strcmp(name,"StokesIceFront")==0) return StokesIceFrontEnum;
+              if (strcmp(name,"SpcTransient")==0) return SpcTransientEnum;
+              else if (strcmp(name,"StokesIceFront")==0) return StokesIceFrontEnum;
               else if (strcmp(name,"StringArrayParam")==0) return StringArrayParamEnum;
               else if (strcmp(name,"StringParam")==0) return StringParamEnum;
@@ -505 +506 @@
               else if (strcmp(name,"PentaP1ElementResult")==0) return PentaP1ElementResultEnum;
               else if (strcmp(name,"StringExternalResult")==0) return StringExternalResultEnum;
-              else if (strcmp(name,"Step")==0) return StepEnum;
          else stage=5;
    }
    if(stage==5){
-              if (strcmp(name,"Time")==0) return TimeEnum;
+              if (strcmp(name,"Step")==0) return StepEnum;
+              else if (strcmp(name,"Time")==0) return TimeEnum;
               else if (strcmp(name,"TriaP1ElementResult")==0) return TriaP1ElementResultEnum;
               else if (strcmp(name,"WaterColumnOld")==0) return WaterColumnOldEnum;

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

@@ -154 +154 @@
                         duf->AYPX(uf_sed,-1.0);
                         ndu_sed=duf->Norm(NORM_TWO); nu_sed=aged_uf_sed->Norm(NORM_TWO);
-                        if (xIsNan<IssmDouble>(ndu_sed) || xIsNan<IssmDouble>(nu_sed)) _error_("convergence criterion is NaN!");
-                        if(isefficientlayer){
-                                duf=aged_uf_epl->Duplicate(); aged_uf_epl->Copy(duf); duf->AYPX(uf_epl,-1.0);
-                                ndu_epl=duf->Norm(NORM_TWO); nu_epl=aged_uf_epl->Norm(NORM_TWO);
-                                if (xIsNan<IssmDouble>(ndu_epl) || xIsNan<IssmDouble>(nu_epl)) _error_("convergence criterion is NaN!");
-                        }
-                        //print
+                        if (xIsNan<IssmDouble>(ndu_sed) || xIsNan<IssmDouble>(nu_sed)) _error_("Sed convergence criterion is NaN!");
                         if (!xIsNan<IssmDouble>(eps_hyd)){
-                                if((ndu_sed/nu_sed)<eps_hyd){
-                                        if(VerboseConvergence()) _printf0_(setw(50) << left << "   Sediment Convergence criterion:" << ndu_sed/nu_sed*100 << " < " << eps_hyd*100 << " %\n");
-                                        hydroconverged=true;
-                                }
-                                else{ 
-                                        if(VerboseConvergence()) _printf0_(setw(50) << left << "   Sediment Convergence criterion:" << ndu_sed/nu_sed*100 << " > " << eps_hyd*100 << " %\n");
-                                        hydroconverged=false;
-                                }
-                                if(isefficientlayer){
-                                        if((ndu_epl/nu_epl)<eps_hyd){
-                                                if(VerboseConvergence()) _printf0_(setw(50) << left << "   EPL Convergence criterion:" << ndu_epl/nu_epl*100 << " < " << eps_hyd*100 << " %\n");
+                                if (!isefficientlayer){
+                                        if ((ndu_sed/nu_sed)<eps_hyd){
+                                                if(VerboseConvergence()) _printf0_(setw(50) << left << "   Converged");
+                                                hydroconverged=true;
                                         }
                                         else{ 
-                                                if(VerboseConvergence()) _printf0_(setw(50) << left << "   EPL Convergence criterion:" << ndu_epl/nu_epl*100 << " > " << eps_hyd*100 << " %\n");
+                                                if(VerboseConvergence()) _printf0_(setw(50) << left << "   Sediment Convergence criterion:" << ndu_sed/nu_sed*100 << " > " << eps_hyd*100 << " %\n");
+                                                hydroconverged=false;
+                                        }
+                                }
+                                else{
+                                        duf=aged_uf_epl->Duplicate(); aged_uf_epl->Copy(duf); duf->AYPX(uf_epl,-1.0);
+                                        ndu_epl=duf->Norm(NORM_TWO); nu_epl=aged_uf_epl->Norm(NORM_TWO);
+                                        if (xIsNan<IssmDouble>(ndu_epl) || xIsNan<IssmDouble>(nu_epl)) _error_("EPL 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 ((ndu_epl/nu_epl)<eps_hyd && (ndu_sed/nu_sed)<eps_hyd){
+                                                if (VerboseConvergence()) _printf0_(setw(50) << left << "   Converged");
+                                                hydroconverged=true;
+                                        }
+                                        else{ 
+                                                if(VerboseConvergence()) _printf0_(setw(50) << left << "   Sediment Convergence criterion:" << ndu_sed/nu_sed*100 << " aiming lower than " << eps_hyd*100 << " % ;"<<ndu_sed <<";"<<nu_sed <<"\n");
+                                                if(VerboseConvergence()) _printf0_(setw(50) << left << "   EPL Convergence criterion:" << ndu_epl/nu_epl*100 << " aiming lower than " << eps_hyd*100 << ";"<<ndu_epl <<";"<<nu_epl <<" %\n");
                                                 hydroconverged=false;
                                         }
@@ -186 +188 @@
                 }
                 hydrocount++;
-
                 if(hydroconverged)break;
         }
 
         InputUpdateFromSolutionx(femmodel->elements,femmodel->nodes, femmodel->vertices, femmodel->loads, femmodel->materials, femmodel->parameters,ug);
-
+        
         /*Free ressources: */
         delete ug;