Changeset 14928

Timestamp:

05/06/13 13:02:15 (12 years ago)

Author:

bdef

Message:

CHG: Computation of the residual for the inefficent hydrological drainage system

Location:

issm/trunk-jpl/src/c

Files:

• EnumDefinitions/EnumDefinitions.h (modified) (1 diff)
• classes/objects/Elements/Tria.cpp (modified) (6 diffs)
• modules/EnumToStringx/EnumToStringx.cpp (modified) (1 diff)
• modules/StringToEnumx/StringToEnumx.cpp (modified) (5 diffs)
• solutions/hydrology_core.cpp (modified) (1 diff)
• solvers/solver_hydro_nonlinear.cpp (modified) (1 diff)

issm/trunk-jpl/src/c/EnumDefinitions/EnumDefinitions.h ¶

@@ -90 +90 @@
         HydrologydcEnum,
         SedimentHeadEnum,
+        SedimentHeadResidualEnum,
         EplHeadEnum,
         HydrologydcSpcsedimentHeadEnum,

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

@@ -2052 +2052 @@
                                 name==GradientEnum ||
                                 name==OldGradientEnum ||
+        name==ConvergedEnum ||
                                 name==QmuVxEnum ||
                                 name==QmuVyEnum ||
@@ -6201 +6202 @@
         IssmDouble xyz_list[NUMVERTICES][3];
         IssmDouble dt,scalar,water_head;
-        IssmDouble water_load;
+        IssmDouble residual_load;
         IssmDouble epl_storing;
         IssmDouble basis[numdof];
@@ -6213 +6214 @@
         GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
         this->parameters->FindParam(&dt,TimesteppingTimeStepEnum);
-        Input* water_input=inputs->GetInput(BasalforcingsMeltingRateEnum);  _assert_(water_input);
+        Input* residual_input=inputs->GetInput(SedimentHeadResidualEnum);  _assert_(residual_input);
         Input* old_wh_input=NULL; 
 
@@ -6230 +6231 @@
 
                 /*Loading term*/
-                water_input->GetInputValue(&water_load,gauss);
-                scalar = Jdet*gauss->weight*water_load;
+                residual_input->GetInputValue(&residual_load,gauss);
+                scalar = Jdet*gauss->weight*residual_load;
                 if(reCast<bool,IssmDouble>(dt)) scalar = scalar*dt;
                 for(int i=0;i<numdof;i++) pe->values[i]+=scalar*basis[i];
@@ -6315 +6316 @@
         const int   numdof         = NDOF1 *NUMVERTICES;
         int        *doflist        = NULL;
+        bool  converged;
         IssmDouble  values[numdof];
+        IssmDouble  residual[numdof];
+        IssmDouble  sediment_storing;
+        IssmDouble  penalty_factor;
+        IssmDouble  kmax, kappa, h_max;
 
         /*Get dof list: */
         GetDofList(&doflist,NoneApproximationEnum,GsetEnum);
+
 
         /*Use the dof list to index into the solution vector: */
@@ -6326 +6333 @@
         }
 
+        /*If converged keep the residual in mind*/
+        this->inputs->GetInputValue(&converged,ConvergedEnum);
+        if(converged){
+                this->parameters->FindParam(&kmax,HydrologySedimentKmaxEnum);
+                this->parameters->FindParam(&penalty_factor,HydrologydcPenaltyFactorEnum);
+                sediment_storing=matpar->GetSedimentStoring();
+                kappa=kmax*pow(10.,penalty_factor);
+                
+                for(int i=0;i<NUMVERTICES;i++){
+                        this->GetHydrologyDCInefficientHmax(&h_max,nodes[i]);
+                        if(values[i]>h_max)
+                                residual[i]=kappa*(values[i]-h_max)*sediment_storing;
+                        else
+                                residual[i]=0.0;
+                }
+        }
+
         /*Add input to the element: */
         this->inputs->AddInput(new TriaP1Input(SedimentHeadEnum,values));
+        this->inputs->AddInput(new TriaP1Input(SedimentHeadResidualEnum,residual));
 
         /*Free ressources:*/

issm/trunk-jpl/src/c/modules/EnumToStringx/EnumToStringx.cpp ¶

@@ -94 +94 @@
                 case HydrologydcEnum : return "Hydrologydc";
                 case SedimentHeadEnum : return "SedimentHead";
+                case SedimentHeadResidualEnum : return "SedimentHeadResidual";
                 case EplHeadEnum : return "EplHead";
                 case HydrologydcSpcsedimentHeadEnum : return "HydrologydcSpcsedimentHead";

issm/trunk-jpl/src/c/modules/StringToEnumx/StringToEnumx.cpp ¶

@@ -95 +95 @@
               else if (strcmp(name,"Hydrologydc")==0) return HydrologydcEnum;
               else if (strcmp(name,"SedimentHead")==0) return SedimentHeadEnum;
+              else if (strcmp(name,"SedimentHeadResidual")==0) return SedimentHeadResidualEnum;
               else if (strcmp(name,"EplHead")==0) return EplHeadEnum;
               else if (strcmp(name,"HydrologydcSpcsedimentHead")==0) return HydrologydcSpcsedimentHeadEnum;
@@ -136 +137 @@
               else if (strcmp(name,"InversionVxObs")==0) return InversionVxObsEnum;
               else if (strcmp(name,"InversionVyObs")==0) return InversionVyObsEnum;
-              else if (strcmp(name,"InversionVzObs")==0) return InversionVzObsEnum;
          else stage=2;
    }
    if(stage==2){
-              if (strcmp(name,"MaskElementonfloatingice")==0) return MaskElementonfloatingiceEnum;
+              if (strcmp(name,"InversionVzObs")==0) return InversionVzObsEnum;
+              else if (strcmp(name,"MaskElementonfloatingice")==0) return MaskElementonfloatingiceEnum;
               else if (strcmp(name,"MaskElementongroundedice")==0) return MaskElementongroundediceEnum;
               else if (strcmp(name,"MaskElementonwater")==0) return MaskElementonwaterEnum;
@@ -259 +260 @@
               else if (strcmp(name,"TransientIsprognostic")==0) return TransientIsprognosticEnum;
               else if (strcmp(name,"TransientIsthermal")==0) return TransientIsthermalEnum;
-              else if (strcmp(name,"TransientIsgia")==0) return TransientIsgiaEnum;
          else stage=3;
    }
    if(stage==3){
-              if (strcmp(name,"TransientNumRequestedOutputs")==0) return TransientNumRequestedOutputsEnum;
+              if (strcmp(name,"TransientIsgia")==0) return TransientIsgiaEnum;
+              else if (strcmp(name,"TransientNumRequestedOutputs")==0) return TransientNumRequestedOutputsEnum;
               else if (strcmp(name,"TransientRequestedOutputs")==0) return TransientRequestedOutputsEnum;
               else if (strcmp(name,"SolutionType")==0) return SolutionTypeEnum;
@@ -382 +383 @@
               else if (strcmp(name,"Vertex")==0) return VertexEnum;
               else if (strcmp(name,"Air")==0) return AirEnum;
-              else if (strcmp(name,"Ice")==0) return IceEnum;
          else stage=4;
    }
    if(stage==4){
-              if (strcmp(name,"Melange")==0) return MelangeEnum;
+              if (strcmp(name,"Ice")==0) return IceEnum;
+              else if (strcmp(name,"Melange")==0) return MelangeEnum;
               else if (strcmp(name,"Water")==0) return WaterEnum;
               else if (strcmp(name,"Closed")==0) return ClosedEnum;
@@ -505 +506 @@
               else if (strcmp(name,"MaxVel")==0) return MaxVelEnum;
               else if (strcmp(name,"MinVx")==0) return MinVxEnum;
-              else if (strcmp(name,"MaxVx")==0) return MaxVxEnum;
          else stage=5;
    }
    if(stage==5){
-              if (strcmp(name,"MaxAbsVx")==0) return MaxAbsVxEnum;
+              if (strcmp(name,"MaxVx")==0) return MaxVxEnum;
+              else if (strcmp(name,"MaxAbsVx")==0) return MaxAbsVxEnum;
               else if (strcmp(name,"MinVy")==0) return MinVyEnum;
               else if (strcmp(name,"MaxVy")==0) return MaxVyEnum;

issm/trunk-jpl/src/c/solutions/hydrology_core.cpp ¶

@@ -84 +84 @@
                                 if(VerboseSolution()) _pprintLine_("   saving results ");
                                 InputToResultx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,SedimentHeadEnum);
+                                InputToResultx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,SedimentHeadResidualEnum);
 
                                 /*unload results*/

issm/trunk-jpl/src/c/solvers/solver_hydro_nonlinear.cpp ¶

@@ -67 +67 @@
 
                         if(converged){
+                                femmodel->parameters->SetParam(sediment_kmax,HydrologySedimentKmaxEnum);
+
                                 InputUpdateFromConstantx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,converged,ConvergedEnum);
-                                InputUpdateFromConstantx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,sediment_kmax,HydrologySedimentKmaxEnum);
                                 InputUpdateFromSolutionx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,ug);
+                                /*InputUpdateFromConstantx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,sediment_kmax,HydrologySedimentKmaxEnum);*/
                                 break;
                         }