Changeset 17952


Ignore:
Timestamp:
05/06/14 17:00:23 (11 years ago)
Author:
cborstad
Message:

CHG: added frictionweertman class for Weertman sliding law

Location:
issm/trunk-jpl/src
Files:
5 added
13 edited

Legend:

Unmodified
Added
Removed

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

    r17946 r17952  
    1717        parameters->AddObject(iomodel->CopyConstantObject(ThermalIsenthalpyEnum));
    1818        parameters->AddObject(iomodel->CopyConstantObject(ThermalIsdynamicbasalspcEnum));
     19        parameters->AddObject(iomodel->CopyConstantObject(FrictionLawEnum));
    1920
    2021        iomodel->FetchData(&requestedoutputs,&numoutputs,ThermalRequestedOutputsEnum);
    … …  
    2526void EnthalpyAnalysis::UpdateElements(Elements* elements,IoModel* iomodel,int analysis_counter,int analysis_type){/*{{{*/
    2627
    27         bool dakota_analysis, islevelset;
    28         bool isenthalpy;
     28        bool dakota_analysis,islevelset,isenthalpy;
     29        int frictionlaw;
    2930
    3031        /*Now, is the model 3d? otherwise, do nothing: */
    … …  
    5455        iomodel->FetchDataToInput(elements,SurfaceEnum);
    5556        iomodel->FetchDataToInput(elements,BaseEnum);
    56         iomodel->FetchDataToInput(elements,FrictionCoefficientEnum);
    57         iomodel->FetchDataToInput(elements,FrictionPEnum);
    58         iomodel->FetchDataToInput(elements,FrictionQEnum);
    5957        iomodel->FetchDataToInput(elements,MaskIceLevelsetEnum);
    6058        iomodel->FetchDataToInput(elements,MaskGroundediceLevelsetEnum);
    … …  
    9088                iomodel->FetchDataToInput(elements,MeshVertexonbaseEnum); // required for updating active nodes
    9189        }
    92 
     90       
     91        /*Friction law variables*/
     92        switch(frictionlaw){
     93                case 1:
     94                        iomodel->FetchDataToInput(elements,FrictionCoefficientEnum);
     95                        iomodel->FetchDataToInput(elements,FrictionPEnum);
     96                        iomodel->FetchDataToInput(elements,FrictionQEnum);
     97                        break;
     98                case 2:
     99                        iomodel->FetchDataToInput(elements,FrictionCEnum);
     100                        iomodel->FetchDataToInput(elements,FrictionMEnum);
     101                        break;
     102                default:
     103                        _error_("not supported");
     104        }
    93105        /*Free data: */
    94106        iomodel->DeleteData(3,TemperatureEnum,WaterfractionEnum,PressureEnum);

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

    r17886 r17952  
    1010}/*}}}*/
    1111void MeltingAnalysis::UpdateParameters(Parameters* parameters,IoModel* iomodel,int solution_enum,int analysis_enum){/*{{{*/
     12        parameters->AddObject(iomodel->CopyConstantObject(FrictionLawEnum));
    1213}/*}}}*/
    1314void MeltingAnalysis::UpdateElements(Elements* elements,IoModel* iomodel,int analysis_counter,int analysis_type){/*{{{*/
     15
     16        int frictionlaw;
    1417
    1518        /*Now, is the model 3d? otherwise, do nothing: */
    … …  
    3033        iomodel->FetchDataToInput(elements,SurfaceEnum);
    3134        iomodel->FetchDataToInput(elements,BaseEnum);
    32         iomodel->FetchDataToInput(elements,FrictionCoefficientEnum);
    33         iomodel->FetchDataToInput(elements,FrictionPEnum);
    34         iomodel->FetchDataToInput(elements,FrictionQEnum);
    3535        iomodel->FetchDataToInput(elements,MaskIceLevelsetEnum);
    3636        if(iomodel->domaintype!=Domain2DhorizontalEnum){
    … …  
    4343        iomodel->FetchDataToInput(elements,BasalforcingsMeltingRateEnum);
    4444        iomodel->FetchDataToInput(elements,PressureEnum);
     45       
     46        /*Friction law variables*/
     47        switch(frictionlaw){
     48                case 1:
     49                        iomodel->FetchDataToInput(elements,FrictionCoefficientEnum);
     50                        iomodel->FetchDataToInput(elements,FrictionPEnum);
     51                        iomodel->FetchDataToInput(elements,FrictionQEnum);
     52                        break;
     53                case 2:
     54                        iomodel->FetchDataToInput(elements,FrictionCEnum);
     55                        iomodel->FetchDataToInput(elements,FrictionMEnum);
     56                        break;
     57                default:
     58                        _error_("not supported");
     59        }
    4560}/*}}}*/
    4661void MeltingAnalysis::CreateNodes(Nodes* nodes,IoModel* iomodel){/*{{{*/

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

    r17949 r17952  
    106106        parameters->AddObject(iomodel->CopyConstantObject(StressbalanceShelfDampeningEnum));
    107107        parameters->AddObject(iomodel->CopyConstantObject(StressbalanceViscosityOvershootEnum));
     108        parameters->AddObject(iomodel->CopyConstantObject(FrictionLawEnum));
    108109
    109110        /*XTH parameters*/
    … …  
    125126        /*Intermediaries*/
    126127        int    materials_type,finiteelement;
    127         int    approximation;
     128        int    approximation,frictionlaw;
    128129        int*   finiteelement_list=NULL;
    129130        bool   isSSA,isL1L2,isHO,isFS,iscoupling;
    … …  
    142143        iomodel->Constant(&materials_type,MaterialsEnum);
    143144        iomodel->Constant(&islevelset,TransientIslevelsetEnum);
     145        iomodel->Constant(&frictionlaw,FrictionLawEnum);
    144146
    145147        /*return if no processing required*/
    … …  
    205207        iomodel->FetchDataToInput(elements,SurfaceEnum);
    206208        iomodel->FetchDataToInput(elements,BaseEnum);
    207         iomodel->FetchDataToInput(elements,FrictionCoefficientEnum);
    208         iomodel->FetchDataToInput(elements,FrictionPEnum);
    209         iomodel->FetchDataToInput(elements,FrictionQEnum);
    210209        iomodel->FetchDataToInput(elements,MaskIceLevelsetEnum);
    211210        iomodel->FetchDataToInput(elements,MaskGroundediceLevelsetEnum);
    … …  
    238237        if(islevelset){
    239238                iomodel->FetchDataToInput(elements,IceMaskNodeActivationEnum);
     239        }
     240
     241        /*Friction law variables*/
     242        switch(frictionlaw){
     243                case 1:
     244                        iomodel->FetchDataToInput(elements,FrictionCoefficientEnum);
     245                        iomodel->FetchDataToInput(elements,FrictionPEnum);
     246                        iomodel->FetchDataToInput(elements,FrictionQEnum);
     247                        break;
     248                case 2:
     249                        iomodel->FetchDataToInput(elements,FrictionCEnum);
     250                        iomodel->FetchDataToInput(elements,FrictionMEnum);
     251                        break;
     252                default:
     253                        _error_("not supported");
    240254        }
    241255

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

    r17946 r17952  
    2121        parameters->AddObject(iomodel->CopyConstantObject(ThermalIsenthalpyEnum));
    2222        parameters->AddObject(iomodel->CopyConstantObject(ThermalIsdynamicbasalspcEnum));
     23        parameters->AddObject(iomodel->CopyConstantObject(FrictionLawEnum));
    2324
    2425        iomodel->FetchData(&requestedoutputs,&numoutputs,ThermalRequestedOutputsEnum);
    … …  
    2930}/*}}}*/
    3031void ThermalAnalysis::UpdateElements(Elements* elements,IoModel* iomodel,int analysis_counter,int analysis_type){/*{{{*/
     32
     33        int frictionlaw;
    3134
    3235        /*Now, is the model 3d? otherwise, do nothing: */
    … …  
    5154        iomodel->FetchDataToInput(elements,SurfaceEnum);
    5255        iomodel->FetchDataToInput(elements,BaseEnum);
    53         iomodel->FetchDataToInput(elements,FrictionCoefficientEnum);
    54         iomodel->FetchDataToInput(elements,FrictionPEnum);
    55         iomodel->FetchDataToInput(elements,FrictionQEnum);
    5656        iomodel->FetchDataToInput(elements,MaskIceLevelsetEnum);
    5757        iomodel->FetchDataToInput(elements,MaskGroundediceLevelsetEnum);
    … …  
    8484                iomodel->FetchDataToInput(elements,IceMaskNodeActivationEnum);
    8585                iomodel->FetchDataToInput(elements,MeshVertexonbaseEnum); // required for updating active nodes
     86        }
     87        /*Friction law variables*/
     88        switch(frictionlaw){
     89                case 1:
     90                        iomodel->FetchDataToInput(elements,FrictionCoefficientEnum);
     91                        iomodel->FetchDataToInput(elements,FrictionPEnum);
     92                        iomodel->FetchDataToInput(elements,FrictionQEnum);
     93                        break;
     94                case 2:
     95                        iomodel->FetchDataToInput(elements,FrictionCEnum);
     96                        iomodel->FetchDataToInput(elements,FrictionMEnum);
     97                        break;
     98                default:
     99                        _error_("not supported");
    86100        }
    87101}/*}}}*/

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

    r17946 r17952  
    2020        this->element=NULL;
    2121        this->dim=0;
     22        this->law=0;
     23
    2224}
    2325/*}}}*/
    … …  
    2729        this->element=element_in;
    2830        this->dim=dim_in;
     31        element_in->FindParam(&this->law,FrictionLawEnum);
    2932}
    3033/*}}}*/
    … …  
    4346/*FUNCTION Friction::GetAlpha2{{{*/
    4447void Friction::GetAlpha2(IssmDouble* palpha2, Gauss* gauss){
     48
     49        switch(this->law){
     50                case 1:
     51                        GetAlpha2Viscous(palpha2,gauss);
     52                        break;
     53                case 2:
     54                        GetAlpha2Weertman(palpha2,gauss);
     55                        break;
     56                default:
     57                        _error_("not supported");
     58        }
     59
     60}/*}}}*/
     61/*FUNCTION Friction::GetAlpha2{{{*/
     62void Friction::GetAlpha2Viscous(IssmDouble* palpha2, Gauss* gauss){
    4563
    4664        /*This routine calculates the basal friction coefficient
    … …  
    102120        *palpha2=alpha2;
    103121}/*}}}*/
     122/*FUNCTION Friction::GetAlpha2Weertman{{{*/
     123void Friction::GetAlpha2Weertman(IssmDouble* palpha2, Gauss* gauss){
     124
     125        /*This routine calculates the basal friction coefficient alpha2= C^-1/m |v|^(1/m-1) */
     126
     127        /*diverse: */
     128        IssmDouble  C,m;
     129        IssmDouble  vx,vy,vz,vmag;
     130        IssmDouble  alpha2;
     131
     132        /*Recover parameters: */
     133        element->GetInputValue(&C,gauss,FrictionCEnum);
     134        element->GetInputValue(&m,FrictionMEnum);
     135
     136        switch(dim){
     137                case 1:
     138                        element->GetInputValue(&vx,gauss,VxEnum);
     139                        vmag=sqrt(vx*vx);
     140                        break;
     141                case 2:
     142                        element->GetInputValue(&vx,gauss,VxEnum);
     143                        element->GetInputValue(&vy,gauss,VyEnum);
     144                        vmag=sqrt(vx*vx+vy*vy);
     145                        break;
     146                case 3:
     147                        element->GetInputValue(&vx,gauss,VxEnum);
     148                        element->GetInputValue(&vy,gauss,VyEnum);
     149                        element->GetInputValue(&vz,gauss,VzEnum);
     150                        vmag=sqrt(vx*vx+vy*vy+vz*vz);
     151                        break;
     152                default:
     153                        _error_("not supported");
     154        }
     155
     156        /*Check to prevent dividing by zero if vmag==0*/
     157        if(vmag==0. && (1./m-1.)<0.) alpha2=0.;
     158        else alpha2=pow(C,-1./m)*pow(vmag,(1./m-1.));
     159        _assert_(!xIsNan<IssmDouble>(alpha2));
     160
     161        /*Assign output pointers:*/
     162        *palpha2=alpha2;
     163}/*}}}*/
    104164/*FUNCTION Friction::GetAlphaComplement(IssmDouble* palpha_complement, Gauss* gauss,int vxenum,int vyenum,int vzenum) {{{*/
    105165void Friction::GetAlphaComplement(IssmDouble* palpha_complement, Gauss* gauss){
    … …  
    108168         * FrictionGetAlphaComplement is used in control methods on drag, and it computes:
    109169         * alpha_complement= Neff ^r * vel ^s*/
     170
     171        if(this->law!=1)_error_("not supported");
    110172
    111173        /*diverse: */

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

    r17943 r17952  
    2121                Element* element;
    2222                int      dim;
     23                int      law;
    2324
    2425                /*methods: */
    … …  
    2930                void  Echo(void);
    3031                void  GetAlpha2(IssmDouble* palpha2,Gauss* gauss);
     32                void  GetAlpha2Viscous(IssmDouble* palpha2,Gauss* gauss);
     33                void  GetAlpha2Weertman(IssmDouble* palpha2,Gauss* gauss);
    3134                void  GetAlphaComplement(IssmDouble* alpha_complement,Gauss* gauss);
    3235};

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

    r17946 r17952  
    8484        FrictionPEnum,
    8585        FrictionQEnum,
     86        FrictionMEnum,
     87        FrictionCEnum,
     88        FrictionLawEnum,
    8689        GeometryHydrostaticRatioEnum,
    8790        HydrologyModelEnum,

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

    r17946 r17952  
    9292                case FrictionPEnum : return "FrictionP";
    9393                case FrictionQEnum : return "FrictionQ";
     94                case FrictionMEnum : return "FrictionM";
     95                case FrictionCEnum : return "FrictionC";
     96                case FrictionLawEnum : return "FrictionLaw";
    9497                case GeometryHydrostaticRatioEnum : return "GeometryHydrostaticRatio";
    9598                case HydrologyModelEnum : return "HydrologyModel";

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

    r17946 r17952  
    9292              else if (strcmp(name,"FrictionP")==0) return FrictionPEnum;
    9393              else if (strcmp(name,"FrictionQ")==0) return FrictionQEnum;
     94              else if (strcmp(name,"FrictionM")==0) return FrictionMEnum;
     95              else if (strcmp(name,"FrictionC")==0) return FrictionCEnum;
     96              else if (strcmp(name,"FrictionLaw")==0) return FrictionLawEnum;
    9497              else if (strcmp(name,"GeometryHydrostaticRatio")==0) return GeometryHydrostaticRatioEnum;
    9598              else if (strcmp(name,"HydrologyModel")==0) return HydrologyModelEnum;
    … …  
    134137              else if (strcmp(name,"HydrologyLayer")==0) return HydrologyLayerEnum;
    135138              else if (strcmp(name,"HydrologySediment")==0) return HydrologySedimentEnum;
    136               else if (strcmp(name,"HydrologyEfficient")==0) return HydrologyEfficientEnum;
    137               else if (strcmp(name,"HydrologySedimentKmax")==0) return HydrologySedimentKmaxEnum;
    138               else if (strcmp(name,"IndependentObject")==0) return IndependentObjectEnum;
    139139         else stage=2;
    140140   }
    141141   if(stage==2){
    142               if (strcmp(name,"InversionControlParameters")==0) return InversionControlParametersEnum;
     142              if (strcmp(name,"HydrologyEfficient")==0) return HydrologyEfficientEnum;
     143              else if (strcmp(name,"HydrologySedimentKmax")==0) return HydrologySedimentKmaxEnum;
     144              else if (strcmp(name,"IndependentObject")==0) return IndependentObjectEnum;
     145              else if (strcmp(name,"InversionControlParameters")==0) return InversionControlParametersEnum;
    143146              else if (strcmp(name,"InversionCostFunctionThreshold")==0) return InversionCostFunctionThresholdEnum;
    144147              else if (strcmp(name,"InversionCostFunctionsCoefficients")==0) return InversionCostFunctionsCoefficientsEnum;
    … …  
    257260              else if (strcmp(name,"QmuMaterialsRheologyB")==0) return QmuMaterialsRheologyBEnum;
    258261              else if (strcmp(name,"RiftsNumrifts")==0) return RiftsNumriftsEnum;
    259               else if (strcmp(name,"RiftsRiftstruct")==0) return RiftsRiftstructEnum;
    260               else if (strcmp(name,"SettingsResultsOnNodes")==0) return SettingsResultsOnNodesEnum;
    261               else if (strcmp(name,"SettingsIoGather")==0) return SettingsIoGatherEnum;
    262262         else stage=3;
    263263   }
    264264   if(stage==3){
    265               if (strcmp(name,"SettingsLowmem")==0) return SettingsLowmemEnum;
     265              if (strcmp(name,"RiftsRiftstruct")==0) return RiftsRiftstructEnum;
     266              else if (strcmp(name,"SettingsResultsOnNodes")==0) return SettingsResultsOnNodesEnum;
     267              else if (strcmp(name,"SettingsIoGather")==0) return SettingsIoGatherEnum;
     268              else if (strcmp(name,"SettingsLowmem")==0) return SettingsLowmemEnum;
    266269              else if (strcmp(name,"SettingsOutputFrequency")==0) return SettingsOutputFrequencyEnum;
    267270              else if (strcmp(name,"SettingsWaitonlock")==0) return SettingsWaitonlockEnum;
    … …  
    380383              else if (strcmp(name,"SmoothedSurfaceSlopeYAnalysis")==0) return SmoothedSurfaceSlopeYAnalysisEnum;
    381384              else if (strcmp(name,"ThermalAnalysis")==0) return ThermalAnalysisEnum;
    382               else if (strcmp(name,"ThermalSolution")==0) return ThermalSolutionEnum;
    383               else if (strcmp(name,"TransientSolution")==0) return TransientSolutionEnum;
    384               else if (strcmp(name,"GiaSolution")==0) return GiaSolutionEnum;
    385385         else stage=4;
    386386   }
    387387   if(stage==4){
    388               if (strcmp(name,"GiaAnalysis")==0) return GiaAnalysisEnum;
     388              if (strcmp(name,"ThermalSolution")==0) return ThermalSolutionEnum;
     389              else if (strcmp(name,"TransientSolution")==0) return TransientSolutionEnum;
     390              else if (strcmp(name,"GiaSolution")==0) return GiaSolutionEnum;
     391              else if (strcmp(name,"GiaAnalysis")==0) return GiaAnalysisEnum;
    389392              else if (strcmp(name,"MeshdeformationSolution")==0) return MeshdeformationSolutionEnum;
    390393              else if (strcmp(name,"MeshdeformationAnalysis")==0) return MeshdeformationAnalysisEnum;
    … …  
    503506              else if (strcmp(name,"QmuSurface")==0) return QmuSurfaceEnum;
    504507              else if (strcmp(name,"QmuMelting")==0) return QmuMeltingEnum;
    505               else if (strcmp(name,"QmuVxMesh")==0) return QmuVxMeshEnum;
    506               else if (strcmp(name,"QmuVyMesh")==0) return QmuVyMeshEnum;
    507               else if (strcmp(name,"QmuVzMesh")==0) return QmuVzMeshEnum;
    508508         else stage=5;
    509509   }
    510510   if(stage==5){
    511               if (strcmp(name,"AndroidFrictionCoefficient")==0) return AndroidFrictionCoefficientEnum;
     511              if (strcmp(name,"QmuVxMesh")==0) return QmuVxMeshEnum;
     512              else if (strcmp(name,"QmuVyMesh")==0) return QmuVyMeshEnum;
     513              else if (strcmp(name,"QmuVzMesh")==0) return QmuVzMeshEnum;
     514              else if (strcmp(name,"AndroidFrictionCoefficient")==0) return AndroidFrictionCoefficientEnum;
    512515              else if (strcmp(name,"ResetPenalties")==0) return ResetPenaltiesEnum;
    513516              else if (strcmp(name,"SegmentOnIceShelf")==0) return SegmentOnIceShelfEnum;
    … …  
    626629              else if (strcmp(name,"MinVx")==0) return MinVxEnum;
    627630              else if (strcmp(name,"MaxVx")==0) return MaxVxEnum;
    628               else if (strcmp(name,"MaxAbsVx")==0) return MaxAbsVxEnum;
    629               else if (strcmp(name,"MinVy")==0) return MinVyEnum;
    630               else if (strcmp(name,"MaxVy")==0) return MaxVyEnum;
    631631         else stage=6;
    632632   }
    633633   if(stage==6){
    634               if (strcmp(name,"MaxAbsVy")==0) return MaxAbsVyEnum;
     634              if (strcmp(name,"MaxAbsVx")==0) return MaxAbsVxEnum;
     635              else if (strcmp(name,"MinVy")==0) return MinVyEnum;
     636              else if (strcmp(name,"MaxVy")==0) return MaxVyEnum;
     637              else if (strcmp(name,"MaxAbsVy")==0) return MaxAbsVyEnum;
    635638              else if (strcmp(name,"MinVz")==0) return MinVzEnum;
    636639              else if (strcmp(name,"MaxVz")==0) return MaxVzEnum;

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

    r17928 r17952  
    4848                        yts=365.0*24.0*3600.0;
    4949
     50                        WriteData(fid,'enum',FrictionLawEnum,'data',1,'format','Integer');
    5051                        WriteData(fid,'object',obj,'fieldname','coefficient','format','DoubleMat','mattype',1,'forcinglength',md.mesh.numberofvertices+1);
    5152                        %WriteData(fid,'object',obj,'fieldname','coefficient','format','DoubleMat','mattype',1);

  • issm/trunk-jpl/src/m/classes/friction.py

    r17928 r17952  
    4545        # }}}
    4646        def marshall(self,md,fid):    # {{{
     47                WriteData(fid,'enum',FrictionLawEnum(),'data',1,'format','Integer')
    4748                WriteData(fid,'object',self,'fieldname','coefficient','format','DoubleMat','mattype',1)
    4849                WriteData(fid,'object',self,'fieldname','p','format','DoubleMat','mattype',2)

  • issm/trunk-jpl/src/m/enum/EnumDefinitions.py

    r17946 r17952  
    8484def FrictionPEnum(): return StringToEnum("FrictionP")[0]
    8585def FrictionQEnum(): return StringToEnum("FrictionQ")[0]
     86def FrictionMEnum(): return StringToEnum("FrictionM")[0]
     87def FrictionCEnum(): return StringToEnum("FrictionC")[0]
     88def FrictionLawEnum(): return StringToEnum("FrictionLaw")[0]
    8689def GeometryHydrostaticRatioEnum(): return StringToEnum("GeometryHydrostaticRatio")[0]
    8790def HydrologyModelEnum(): return StringToEnum("HydrologyModel")[0]

  • issm/trunk-jpl/src/m/plot/plot_overlay.py

    r17838 r17952  
    22from processmesh import processmesh
    33from processdata import processdata
    4 from osgeo import gdal
    54import matplotlib.pyplot as plt
    65import matplotlib as mpl
    76import os
     7try:
     8        from osgeo import gdal
     9except ImportError:
     10        print 'osgeo/gdal for python not installed, plot_overlay is disabled'
    811
    912def plot_overlay(md,data,options,ax):
Note: See TracChangeset for help on using the changeset viewer.