Changeset 18772


Ignore:
Timestamp:
11/10/14 15:10:13 (10 years ago)
Author:
seroussi
Message:

NEW: added water layer friction law

Location:
issm/trunk-jpl
Files:
4 added
8 edited

Legend:

Unmodified
Added
Removed

  • issm/trunk-jpl/src/c/analyses/EnthalpyAnalysis.cpp

    r18732 r18772  
    107107                        iomodel->FetchDataToInput(elements,PressureEnum);
    108108                        iomodel->FetchDataToInput(elements,TemperatureEnum);
     109                        break;
     110                case 5:
     111                        iomodel->FetchDataToInput(elements,FrictionCoefficientEnum);
     112                        iomodel->FetchDataToInput(elements,FrictionPEnum);
     113                        iomodel->FetchDataToInput(elements,FrictionQEnum);
     114                        iomodel->FetchDataToInput(elements,FrictionWaterLayerEnum);
    109115                        break;
    110116                default:

  • issm/trunk-jpl/src/c/analyses/StressbalanceAnalysis.cpp

    r18732 r18772  
    274274                        iomodel->FetchDataToInput(elements,TemperatureEnum);
    275275                        break;
     276                case 5:
     277                        iomodel->FetchDataToInput(elements,FrictionCoefficientEnum);
     278                        iomodel->FetchDataToInput(elements,FrictionPEnum);
     279                        iomodel->FetchDataToInput(elements,FrictionQEnum);
     280                        iomodel->FetchDataToInput(elements,FrictionWaterLayerEnum);
     281                        break;
    276282                default:
    277283                        _error_("not supported");

  • issm/trunk-jpl/src/c/analyses/ThermalAnalysis.cpp

    r18311 r18772  
    8888                        iomodel->FetchDataToInput(elements,FrictionCEnum);
    8989                        iomodel->FetchDataToInput(elements,FrictionMEnum);
     90                        break;
     91                case 5:
     92                        iomodel->FetchDataToInput(elements,FrictionCoefficientEnum);
     93                        iomodel->FetchDataToInput(elements,FrictionPEnum);
     94                        iomodel->FetchDataToInput(elements,FrictionQEnum);
     95                        iomodel->FetchDataToInput(elements,FrictionWaterLayerEnum);
    9096                        break;
    9197                default:

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

    r18770 r18772  
    264264        IssmDouble  thickness,bed;
    265265        IssmDouble  vx,vy,vz,vmag;
     266        IssmDouble  drag_coefficient,water_layer;
     267        IssmDouble  alpha2;
     268
     269        /*Recover parameters: */
     270        element->GetInputValue(&drag_p,FrictionPEnum);
     271        element->GetInputValue(&drag_q,FrictionQEnum);
     272        element->GetInputValue(&thickness, gauss,ThicknessEnum);
     273        element->GetInputValue(&bed, gauss,BaseEnum);
     274        element->GetInputValue(&drag_coefficient, gauss,FrictionCoefficientEnum);
     275        element->GetInputValue(&water_layer, gauss,FrictionWaterLayerEnum);
     276        IssmDouble rho_water   = element->GetMaterialParameter(MaterialsRhoSeawaterEnum);
     277        IssmDouble rho_ice     = element->GetMaterialParameter(MaterialsRhoIceEnum);
     278        IssmDouble gravity     = element->GetMaterialParameter(ConstantsGEnum);
     279
     280        //compute r and q coefficients: */
     281        r=drag_q/drag_p;
     282        s=1./drag_p;
     283
     284        //From bed and thickness, compute effective pressure when drag is viscous:
     285        Neff=gravity*(rho_ice*thickness+rho_water*(bed+water_layer));
     286        if(Neff<0)Neff=0;
     287
     288        switch(dim){
     289                case 1:
     290                        element->GetInputValue(&vx,gauss,VxEnum);
     291                        vmag=sqrt(vx*vx);
     292                        break;
     293                case 2:
     294                        element->GetInputValue(&vx,gauss,VxEnum);
     295                        element->GetInputValue(&vy,gauss,VyEnum);
     296                        vmag=sqrt(vx*vx+vy*vy);
     297                        break;
     298                case 3:
     299                        element->GetInputValue(&vx,gauss,VxEnum);
     300                        element->GetInputValue(&vy,gauss,VyEnum);
     301                        element->GetInputValue(&vz,gauss,VzEnum);
     302                        vmag=sqrt(vx*vx+vy*vy+vz*vz);
     303                        break;
     304                default:
     305                        _error_("not supported");
     306        }
     307
     308        /*Check to prevent dividing by zero if vmag==0*/
     309        if(vmag==0. && (s-1.)<0.) alpha2=0.;
     310        else alpha2=drag_coefficient*drag_coefficient*pow(Neff,r)*pow(vmag,(s-1.));
     311        _assert_(!xIsNan<IssmDouble>(alpha2));
     312
     313        /*Assign output pointers:*/
     314        *palpha2=alpha2;
     315}/*}}}*/
     316void Friction::GetAlphaComplement(IssmDouble* palpha_complement, Gauss* gauss){/*{{{*/
     317
     318        /* FrictionGetAlpha2 computes alpha2= drag^2 * Neff ^r * vel ^s, with Neff=rho_ice*g*thickness+rho_ice*g*bed, r=q/p and s=1/p.
     319         * FrictionGetAlphaComplement is used in control methods on drag, and it computes:
     320         * alpha_complement= Neff ^r * vel ^s*/
     321
     322        if(this->law!=1)_error_("not supported");
     323
     324        /*diverse: */
     325        IssmDouble  r,s;
     326        IssmDouble  vx,vy,vz,vmag;
     327        IssmDouble  drag_p,drag_q;
     328        IssmDouble  Neff;
    266329        IssmDouble  drag_coefficient;
    267         IssmDouble  alpha2;
     330        IssmDouble  bed,thickness;
     331        IssmDouble  alpha_complement;
    268332
    269333        /*Recover parameters: */
    … …  
    285349        if(Neff<0)Neff=0;
    286350
    287         switch(dim){
    288                 case 1:
    289                         element->GetInputValue(&vx,gauss,VxEnum);
    290                         vmag=sqrt(vx*vx);
    291                         break;
    292                 case 2:
    293                         element->GetInputValue(&vx,gauss,VxEnum);
    294                         element->GetInputValue(&vy,gauss,VyEnum);
    295                         vmag=sqrt(vx*vx+vy*vy);
    296                         break;
    297                 case 3:
    298                         element->GetInputValue(&vx,gauss,VxEnum);
    299                         element->GetInputValue(&vy,gauss,VyEnum);
    300                         element->GetInputValue(&vz,gauss,VzEnum);
    301                         vmag=sqrt(vx*vx+vy*vy+vz*vz);
    302                         break;
    303                 default:
    304                         _error_("not supported");
    305         }
    306 
    307         /*Check to prevent dividing by zero if vmag==0*/
    308         if(vmag==0. && (s-1.)<0.) alpha2=0.;
    309         else alpha2=drag_coefficient*drag_coefficient*pow(Neff,r)*pow(vmag,(s-1.));
    310         _assert_(!xIsNan<IssmDouble>(alpha2));
    311 
    312         /*Assign output pointers:*/
    313         *palpha2=alpha2;
    314 }/*}}}*/
    315 void Friction::GetAlphaComplement(IssmDouble* palpha_complement, Gauss* gauss){/*{{{*/
    316 
    317         /* FrictionGetAlpha2 computes alpha2= drag^2 * Neff ^r * vel ^s, with Neff=rho_ice*g*thickness+rho_ice*g*bed, r=q/p and s=1/p.
    318          * FrictionGetAlphaComplement is used in control methods on drag, and it computes:
    319          * alpha_complement= Neff ^r * vel ^s*/
    320 
    321         if(this->law!=1)_error_("not supported");
    322 
    323         /*diverse: */
    324         IssmDouble  r,s;
    325         IssmDouble  vx,vy,vz,vmag;
    326         IssmDouble  drag_p,drag_q;
    327         IssmDouble  Neff;
    328         IssmDouble  drag_coefficient;
    329         IssmDouble  bed,thickness;
    330         IssmDouble  alpha_complement;
    331 
    332         /*Recover parameters: */
    333         element->GetInputValue(&drag_p,FrictionPEnum);
    334         element->GetInputValue(&drag_q,FrictionQEnum);
    335         element->GetInputValue(&thickness, gauss,ThicknessEnum);
    336         element->GetInputValue(&bed, gauss,BaseEnum);
    337         element->GetInputValue(&drag_coefficient, gauss,FrictionCoefficientEnum);
    338         IssmDouble rho_water   = element->GetMaterialParameter(MaterialsRhoSeawaterEnum);
    339         IssmDouble rho_ice     = element->GetMaterialParameter(MaterialsRhoIceEnum);
    340         IssmDouble gravity     = element->GetMaterialParameter(ConstantsGEnum);
    341 
    342         //compute r and q coefficients: */
    343         r=drag_q/drag_p;
    344         s=1./drag_p;
    345 
    346         //From bed and thickness, compute effective pressure when drag is viscous:
    347         Neff=gravity*(rho_ice*thickness+rho_water*bed);
    348         if(Neff<0)Neff=0;
    349 
    350351        //We need the velocity magnitude to evaluate the basal stress:
    351352        switch(dim){

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

    r18758 r18772  
    9696        FrictionLawEnum,
    9797        FrictionGammaEnum,
     98        FrictionWaterLayerEnum,
    9899        GeometryHydrostaticRatioEnum,
    99100        HydrologyModelEnum,

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

    r18758 r18772  
    104104                case FrictionLawEnum : return "FrictionLaw";
    105105                case FrictionGammaEnum : return "FrictionGamma";
     106                case FrictionWaterLayerEnum : return "FrictionWaterLayer";
    106107                case GeometryHydrostaticRatioEnum : return "GeometryHydrostaticRatio";
    107108                case HydrologyModelEnum : return "HydrologyModel";

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

    r18758 r18772  
    104104              else if (strcmp(name,"FrictionLaw")==0) return FrictionLawEnum;
    105105              else if (strcmp(name,"FrictionGamma")==0) return FrictionGammaEnum;
     106              else if (strcmp(name,"FrictionWaterLayer")==0) return FrictionWaterLayerEnum;
    106107              else if (strcmp(name,"GeometryHydrostaticRatio")==0) return GeometryHydrostaticRatioEnum;
    107108              else if (strcmp(name,"HydrologyModel")==0) return HydrologyModelEnum;
    … …  
    136137              else if (strcmp(name,"HydrologydcEplMaxThickness")==0) return HydrologydcEplMaxThicknessEnum;
    137138              else if (strcmp(name,"HydrologydcEplThickness")==0) return HydrologydcEplThicknessEnum;
    138               else if (strcmp(name,"HydrologydcEplThicknessOld")==0) return HydrologydcEplThicknessOldEnum;
    139139         else stage=2;
    140140   }
    141141   if(stage==2){
    142               if (strcmp(name,"HydrologydcEplConductivity")==0) return HydrologydcEplConductivityEnum;
     142              if (strcmp(name,"HydrologydcEplThicknessOld")==0) return HydrologydcEplThicknessOldEnum;
     143              else if (strcmp(name,"HydrologydcEplConductivity")==0) return HydrologydcEplConductivityEnum;
    143144              else if (strcmp(name,"HydrologydcIsefficientlayer")==0) return HydrologydcIsefficientlayerEnum;
    144145              else if (strcmp(name,"HydrologydcSedimentlimitFlag")==0) return HydrologydcSedimentlimitFlagEnum;
    … …  
    259260              else if (strcmp(name,"MeshY")==0) return MeshYEnum;
    260261              else if (strcmp(name,"MeshZ")==0) return MeshZEnum;
    261               else if (strcmp(name,"MeshElementtype")==0) return MeshElementtypeEnum;
    262262         else stage=3;
    263263   }
    264264   if(stage==3){
    265               if (strcmp(name,"DomainType")==0) return DomainTypeEnum;
     265              if (strcmp(name,"MeshElementtype")==0) return MeshElementtypeEnum;
     266              else if (strcmp(name,"DomainType")==0) return DomainTypeEnum;
    266267              else if (strcmp(name,"DomainDimension")==0) return DomainDimensionEnum;
    267268              else if (strcmp(name,"Domain2Dhorizontal")==0) return Domain2DhorizontalEnum;
    … …  
    382383              else if (strcmp(name,"AdjointBalancethicknessAnalysis")==0) return AdjointBalancethicknessAnalysisEnum;
    383384              else if (strcmp(name,"AdjointBalancethickness2Analysis")==0) return AdjointBalancethickness2AnalysisEnum;
    384               else if (strcmp(name,"AdjointHorizAnalysis")==0) return AdjointHorizAnalysisEnum;
    385385         else stage=4;
    386386   }
    387387   if(stage==4){
    388               if (strcmp(name,"AnalysisCounter")==0) return AnalysisCounterEnum;
     388              if (strcmp(name,"AdjointHorizAnalysis")==0) return AdjointHorizAnalysisEnum;
     389              else if (strcmp(name,"AnalysisCounter")==0) return AnalysisCounterEnum;
    389390              else if (strcmp(name,"DefaultAnalysis")==0) return DefaultAnalysisEnum;
    390391              else if (strcmp(name,"BalancethicknessAnalysis")==0) return BalancethicknessAnalysisEnum;
    … …  
    505506              else if (strcmp(name,"SpcTransient")==0) return SpcTransientEnum;
    506507              else if (strcmp(name,"StringArrayParam")==0) return StringArrayParamEnum;
    507               else if (strcmp(name,"StringParam")==0) return StringParamEnum;
    508508         else stage=5;
    509509   }
    510510   if(stage==5){
    511               if (strcmp(name,"Seg")==0) return SegEnum;
     511              if (strcmp(name,"StringParam")==0) return StringParamEnum;
     512              else if (strcmp(name,"Seg")==0) return SegEnum;
    512513              else if (strcmp(name,"SegInput")==0) return SegInputEnum;
    513514              else if (strcmp(name,"Tria")==0) return TriaEnum;
    … …  
    628629              else if (strcmp(name,"P2")==0) return P2Enum;
    629630              else if (strcmp(name,"P2bubble")==0) return P2bubbleEnum;
    630               else if (strcmp(name,"P2bubblecondensed")==0) return P2bubblecondensedEnum;
    631631         else stage=6;
    632632   }
    633633   if(stage==6){
    634               if (strcmp(name,"P2xP1")==0) return P2xP1Enum;
     634              if (strcmp(name,"P2bubblecondensed")==0) return P2bubblecondensedEnum;
     635              else if (strcmp(name,"P2xP1")==0) return P2xP1Enum;
    635636              else if (strcmp(name,"P1xP2")==0) return P1xP2Enum;
    636637              else if (strcmp(name,"P1xP3")==0) return P1xP3Enum;
    … …  
    751752              else if (strcmp(name,"Arrhenius")==0) return ArrheniusEnum;
    752753              else if (strcmp(name,"LliboutryDuval")==0) return LliboutryDuvalEnum;
    753               else if (strcmp(name,"TransientIslevelset")==0) return TransientIslevelsetEnum;
    754754         else stage=7;
    755755   }
    756756   if(stage==7){
    757               if (strcmp(name,"ExtrapolationVariable")==0) return ExtrapolationVariableEnum;
     757              if (strcmp(name,"TransientIslevelset")==0) return TransientIslevelsetEnum;
     758              else if (strcmp(name,"ExtrapolationVariable")==0) return ExtrapolationVariableEnum;
    758759              else if (strcmp(name,"IceMaskNodeActivation")==0) return IceMaskNodeActivationEnum;
    759760              else if (strcmp(name,"LevelsetfunctionSlopeX")==0) return LevelsetfunctionSlopeXEnum;

  • issm/trunk-jpl/src/m/classes/frictionwaterlayer.m

    r18770 r18772  
    44%      frictionwaterlayer=frictionwaterlayer();
    55
    6 classdef frictiontemp
     6classdef frictionwaterlayer
    77        properties (SetAccess=public)
    88                coefficient = NaN;
    … …  
    1717                                        obj=setdefaultparameters(obj);
    1818                                case 1
    19                                         obj=structtoobj(frictiontemp(),varargin{1});
     19                                        obj=structtoobj(frictionwaterlayer(),varargin{1});
    2020                                otherwise
    2121                                        error('constructor not supported');
Note: See TracChangeset for help on using the changeset viewer.