Changeset 18511

Timestamp:

09/12/14 09:21:30 (11 years ago)

Author:

bdef

Message:

NEW:Adding Hydrology driven friction law

Location:

issm/trunk-jpl/src/c

Files:

• classes/Loads/Friction.cpp (modified) (2 diffs)
• classes/Loads/Friction.h (modified) (1 diff)
• shared/Enum/EnumDefinitions.h (modified) (1 diff)
• shared/Enum/EnumToStringx.cpp (modified) (1 diff)
• shared/Enum/StringToEnumx.cpp (modified) (7 diffs)

issm/trunk-jpl/src/c/classes/Loads/Friction.cpp

@@ -48 +48 @@
                 case 2:
                         GetAlpha2Weertman(palpha2,gauss);
-                        break;
-                default:
+          case 3:
+                        GetAlpha2Hydro(palpha2,gauss);
+                        break;
+          default:
                         _error_("not supported");
         }
@@ -155 +157 @@
         *palpha2=alpha2;
 }/*}}}*/
+void Friction::GetAlpha2Hydro(IssmDouble* palpha2, Gauss* gauss){/*{{{*/
+
+        /*This routine calculates the basal friction coefficient 
+                Based on Gagliardini 2007, needs a good effective pressure computation
+          alpha2= C*Neff*[|vel|/{(C^n*Neff^n*As)*(1+alpha*Chi^q)}]^(1/n)  with
+                -Chi=|vel|/(C^n*Neff^n*As)
+                -alpha=(q-1)^(q-1)/q^q  **/
+
+        /*diverse: */
+        IssmDouble  q_exp;
+        IssmDouble  C_param;
+        IssmDouble  As;
+
+        IssmDouble  Neff;
+        IssmDouble  n;
+
+        IssmDouble  alpha;
+        IssmDouble  Chi;
+
+        IssmDouble  vx,vy,vz,vmag;
+        IssmDouble  alpha2;
+
+        /*Recover parameters: */
+        element->GetInputValue(&q_exp,FrictionQEnum);
+        element->GetInputValue(&C_param,FrictionCEnum);
+        element->GetInputValue(&As,FrictionAsEnum);
+
+        element->GetInputValue(&Neff,gauss,EffectivePressureEnum);
+        element->GetInputValue(&n,gauss,MaterialsRheologyNEnum);
+        
+        if(Neff<0)Neff=0;
+
+        switch(dim){
+                case 1:
+                        element->GetInputValue(&vx,gauss,VxEnum);
+                        vmag=sqrt(vx*vx);
+                        break;
+                case 2:
+                        element->GetInputValue(&vx,gauss,VxEnum);
+                        element->GetInputValue(&vy,gauss,VyEnum);
+                        vmag=sqrt(vx*vx+vy*vy);
+                        break;
+                case 3:
+                        element->GetInputValue(&vx,gauss,VxEnum);
+                        element->GetInputValue(&vy,gauss,VyEnum);
+                        element->GetInputValue(&vz,gauss,VzEnum);
+                        vmag=sqrt(vx*vx+vy*vy+vz*vz);
+                        break;
+                default:
+                        _error_("not supported");
+        }
+
+        //compute alpha and Chi coefficients: */
+        alpha=(pow(q_exp-1,q_exp-1))/pow(q_exp,q_exp);
+        Chi=vmag/(pow(C_param,n)*pow(Neff,n)*As);
+
+        alpha2=C_param*Neff*pow((Chi/(1+alpha*pow(Chi,q_exp))),1./n);
+        _assert_(!xIsNan<IssmDouble>(alpha2));
+        
+        /*Assign output pointers:*/
+        *palpha2=alpha2;
+}/*}}}*/
 void Friction::GetAlphaComplement(IssmDouble* palpha_complement, Gauss* gauss){/*{{{*/
 

issm/trunk-jpl/src/c/classes/Loads/Friction.h

@@ -32 +32 @@
                 void  GetAlpha2Viscous(IssmDouble* palpha2,Gauss* gauss);
                 void  GetAlpha2Weertman(IssmDouble* palpha2,Gauss* gauss);
+                void  GetAlpha2Hydro(IssmDouble* palpha2,Gauss* gauss);
                 void  GetAlphaComplement(IssmDouble* alpha_complement,Gauss* gauss);
 };

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

@@ -88 +88 @@
         FlowequationFeFSEnum,
         FlowequationVertexEquationEnum,
+        FrictionAsEnum,
         FrictionCoefficientEnum,
         FrictionPEnum,

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

@@ -96 +96 @@
                 case FlowequationFeFSEnum : return "FlowequationFeFS";
                 case FlowequationVertexEquationEnum : return "FlowequationVertexEquation";
+                case FrictionAsEnum : return "FrictionAs";
                 case FrictionCoefficientEnum : return "FrictionCoefficient";
                 case FrictionPEnum : return "FrictionP";

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

@@ -96 +96 @@
               else if (strcmp(name,"FlowequationFeFS")==0) return FlowequationFeFSEnum;
               else if (strcmp(name,"FlowequationVertexEquation")==0) return FlowequationVertexEquationEnum;
+              else if (strcmp(name,"FrictionAs")==0) return FrictionAsEnum;
               else if (strcmp(name,"FrictionCoefficient")==0) return FrictionCoefficientEnum;
               else if (strcmp(name,"FrictionP")==0) return FrictionPEnum;
@@ -136 +137 @@
               else if (strcmp(name,"HydrologydcIsefficientlayer")==0) return HydrologydcIsefficientlayerEnum;
               else if (strcmp(name,"HydrologydcSedimentlimitFlag")==0) return HydrologydcSedimentlimitFlagEnum;
-              else if (strcmp(name,"HydrologydcSedimentlimit")==0) return HydrologydcSedimentlimitEnum;
          else stage=2;
    }
    if(stage==2){
-              if (strcmp(name,"HydrologydcTransferFlag")==0) return HydrologydcTransferFlagEnum;
+              if (strcmp(name,"HydrologydcSedimentlimit")==0) return HydrologydcSedimentlimitEnum;
+              else if (strcmp(name,"HydrologydcTransferFlag")==0) return HydrologydcTransferFlagEnum;
               else if (strcmp(name,"HydrologydcLeakageFactor")==0) return HydrologydcLeakageFactorEnum;
               else if (strcmp(name,"HydrologydcPenaltyFactor")==0) return HydrologydcPenaltyFactorEnum;
@@ -259 +260 @@
               else if (strcmp(name,"MassFluxSegments")==0) return MassFluxSegmentsEnum;
               else if (strcmp(name,"MassFluxSegmentsPresent")==0) return MassFluxSegmentsPresentEnum;
-              else if (strcmp(name,"QmuMassFluxSegmentsPresent")==0) return QmuMassFluxSegmentsPresentEnum;
          else stage=3;
    }
    if(stage==3){
-              if (strcmp(name,"QmuNumberofpartitions")==0) return QmuNumberofpartitionsEnum;
+              if (strcmp(name,"QmuMassFluxSegmentsPresent")==0) return QmuMassFluxSegmentsPresentEnum;
+              else if (strcmp(name,"QmuNumberofpartitions")==0) return QmuNumberofpartitionsEnum;
               else if (strcmp(name,"QmuNumberofresponses")==0) return QmuNumberofresponsesEnum;
               else if (strcmp(name,"QmuPartition")==0) return QmuPartitionEnum;
@@ -382 +383 @@
               else if (strcmp(name,"MasstransportAnalysis")==0) return MasstransportAnalysisEnum;
               else if (strcmp(name,"MasstransportSolution")==0) return MasstransportSolutionEnum;
-              else if (strcmp(name,"FreeSurfaceBaseAnalysis")==0) return FreeSurfaceBaseAnalysisEnum;
          else stage=4;
    }
    if(stage==4){
-              if (strcmp(name,"FreeSurfaceTopAnalysis")==0) return FreeSurfaceTopAnalysisEnum;
+              if (strcmp(name,"FreeSurfaceBaseAnalysis")==0) return FreeSurfaceBaseAnalysisEnum;
+              else if (strcmp(name,"FreeSurfaceTopAnalysis")==0) return FreeSurfaceTopAnalysisEnum;
               else if (strcmp(name,"SurfaceNormalVelocity")==0) return SurfaceNormalVelocityEnum;
               else if (strcmp(name,"ExtrudeFromBaseAnalysis")==0) return ExtrudeFromBaseAnalysisEnum;
@@ -505 +506 @@
               else if (strcmp(name,"Friction")==0) return FrictionEnum;
               else if (strcmp(name,"Internal")==0) return InternalEnum;
-              else if (strcmp(name,"MassFlux")==0) return MassFluxEnum;
          else stage=5;
    }
    if(stage==5){
-              if (strcmp(name,"MeltingOffset")==0) return MeltingOffsetEnum;
+              if (strcmp(name,"MassFlux")==0) return MassFluxEnum;
+              else if (strcmp(name,"MeltingOffset")==0) return MeltingOffsetEnum;
               else if (strcmp(name,"Misfit")==0) return MisfitEnum;
               else if (strcmp(name,"Pressure")==0) return PressureEnum;
@@ -628 +629 @@
               else if (strcmp(name,"MisfitTimeinterpolation")==0) return MisfitTimeinterpolationEnum;
               else if (strcmp(name,"MisfitWeights")==0) return MisfitWeightsEnum;
-              else if (strcmp(name,"MisfitWeightsEnum")==0) return MisfitWeightsEnumEnum;
          else stage=6;
    }
    if(stage==6){
-              if (strcmp(name,"SurfaceObservation")==0) return SurfaceObservationEnum;
+              if (strcmp(name,"MisfitWeightsEnum")==0) return MisfitWeightsEnumEnum;
+              else if (strcmp(name,"SurfaceObservation")==0) return SurfaceObservationEnum;
               else if (strcmp(name,"WeightsSurfaceObservation")==0) return WeightsSurfaceObservationEnum;
               else if (strcmp(name,"MinVel")==0) return MinVelEnum;
@@ -751 +752 @@
               else if (strcmp(name,"MaterialsTimeRelaxationDamage")==0) return MaterialsTimeRelaxationDamageEnum;
               else if (strcmp(name,"MaterialsRidgingExponent")==0) return MaterialsRidgingExponentEnum;
-              else if (strcmp(name,"MaterialsCohesion")==0) return MaterialsCohesionEnum;
          else stage=7;
    }
    if(stage==7){
-              if (strcmp(name,"MaterialsInternalFrictionCoef")==0) return MaterialsInternalFrictionCoefEnum;
+              if (strcmp(name,"MaterialsCohesion")==0) return MaterialsCohesionEnum;
+              else if (strcmp(name,"MaterialsInternalFrictionCoef")==0) return MaterialsInternalFrictionCoefEnum;
               else if (strcmp(name,"MaterialsCompressionCoef")==0) return MaterialsCompressionCoefEnum;
               else if (strcmp(name,"MaterialsTractionCoef")==0) return MaterialsTractionCoefEnum;