Changeset 21389


Ignore:
Timestamp:
11/17/16 16:29:42 (8 years ago)
Author:
felicity
Message:

NEW: matenhancedice class, enhancement in Glen flow relationviscosity

Location:
issm/trunk-jpl/src/c
Files:
10 edited

Legend:

Unmodified
Added
Removed

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

    r21382 r21389  
    175175        iomodel->FetchDataToInput(elements,"md.materials.rheology_B",MaterialsRheologyBEnum);
    176176        switch(materialstype){
     177                case MatenhancediceEnum:
     178                        iomodel->FetchDataToInput(elements,"md.materials.rheology_n",MaterialsRheologyNEnum);
     179                        iomodel->FetchDataToInput(elements,"md.materials.rheology_E",MaterialsRheologyEEnum);
     180                        break;
    177181                case MatdamageiceEnum:
    178182                        iomodel->FetchDataToInput(elements,"md.materials.rheology_n",MaterialsRheologyNEnum);

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

    r21382 r21389  
    119119        iomodel->FetchDataToInput(elements,"md.materials.rheology_B",MaterialsRheologyBEnum);
    120120        switch(materialstype){
     121                case MatenhancediceEnum:
     122                        iomodel->FetchDataToInput(elements,"md.materials.rheology_n",MaterialsRheologyNEnum);
     123                        iomodel->FetchDataToInput(elements,"md.materials.rheology_E",MaterialsRheologyEEnum);
     124                        break;
    121125                case MatdamageiceEnum:
    122126                        iomodel->FetchDataToInput(elements,"md.materials.rheology_n",MaterialsRheologyNEnum);

  • issm/trunk-jpl/src/c/classes/Materials/Matice.cpp

    r20966 r21389  
    3131        this->element=NULL;
    3232        this->isdamaged=false;
     33        this->isenhanced=false;
    3334        return;
    3435}
    … …  
    5051   int    materialtype;
    5152   iomodel->FindConstant(&materialtype,"md.materials.type");
    52    if(materialtype==MatdamageiceEnum) this->isdamaged = true;
    53    else if(materialtype==MaticeEnum) this->isdamaged = false;
     53   if(materialtype==MatdamageiceEnum){
     54                this->isdamaged = true;
     55                this->isenhanced = false;
     56        }
     57        else if(materialtype==MaticeEnum){
     58                this->isdamaged = false;
     59                this->isenhanced = false;
     60        }
     61        else if(materialtype==MatenhancediceEnum){
     62                this->isdamaged = false;
     63                this->isenhanced = true;
     64        }
    5465   else _error_("Material type not recognized");
    55 
    5666        return;
    5767
    … …  
    7888        matice->element =(Element*)this->helement->delivers();
    7989        matice->isdamaged = this->isdamaged;
     90        matice->isenhanced = this->isenhanced;
    8091
    8192        return matice;
    … …  
    95106        matice->element =element_in;
    96107        matice->isdamaged = this->isdamaged;
     108        matice->isenhanced = this->isenhanced;
    97109
    98110        return matice;
    … …  
    124136        MARSHALLING(mid);
    125137        MARSHALLING(isdamaged);
     138        MARSHALLING(isenhanced);
    126139        this->helement->Marshall(pmarshalled_data,pmarshalled_data_size,marshall_direction);
    127140        this->element=(Element*)this->helement->delivers();
    … …  
    207220}
    208221/*}}}*/
     222IssmDouble Matice::GetE(){/*{{{*/
     223
     224        _assert_(this->isenhanced);
     225        /*Output*/
     226        IssmDouble E;
     227        if(this->isenhanced)element->inputs->GetInputAverage(&E,MaterialsRheologyEEnum);
     228        return E;
     229}
     230/*}}}*/
     231IssmDouble Matice::GetEbar(){/*{{{*/
     232
     233        _assert_(this->isenhanced);
     234        /*Output*/
     235        IssmDouble Ebar;
     236        if(this->isenhanced)element->inputs->GetInputAverage(&Ebar,MaterialsRheologyEbarEnum);
     237        return Ebar;
     238}
     239/*}}}*/
    209240IssmDouble Matice::GetN(){/*{{{*/
    210241
    … …  
    219250
    220251        return this->isdamaged;
     252}
     253/*}}}*/
     254bool Matice::IsEnhanced(){/*{{{*/
     255
     256        return this->isenhanced;
    221257}
    222258/*}}}*/
    … …  
    225261                                                                (1-D) B
    226262          viscosity= -------------------------
    227                                                   2 eps_eff ^[(n-1)/n]
    228 
    229           where viscosity is the viscotiy, B the flow law parameter , eps_eff is the effective strain rate
    230           and n the flow law exponent.
     263                                                  2 E^[1/n] eps_eff ^[(n-1)/n]
     264
     265          where viscosity is the viscosity, B the flow law parameter , eps_eff is the effective strain rate,
     266          n the flow law exponent, and E is the enhancement factor.
    231267
    232268          If eps_eff = 0 , it means this is the first time SystemMatrices is being run, and we
    … …  
    238274
    239275        /*Intermediary: */
    240         IssmDouble B,D=0.,n;
     276        IssmDouble B,D=0.,E=1.,n;
    241277
    242278        /*Get B and n*/
    … …  
    247283                _assert_(D>=0. && D<1.);
    248284        }
     285        if(this->isenhanced){
     286                E=GetE();
     287                _assert_(E>0.);
     288        }
    249289
    250290        if (n==1.){
    251                 /*Linear Viscous behavior (Newtonian fluid) viscosity=B/2: */
    252                 viscosity=(1.-D)*B/2.;
     291                /*Linear Viscous behavior (Newtonian fluid) viscosity=B/2E: */
     292                viscosity=(1.-D)*B/(2.*E);
    253293        }
    254294        else{
    … …  
    262302
    263303                else{
    264                         viscosity=(1.-D)*B/(2.*pow(eps_eff,(n-1.)/n));
     304                        viscosity=(1.-D)*B/(2.*pow(E,1./n)*pow(eps_eff,(n-1.)/n));
    265305                }
    266306        }
    … …  
    277317                                                                (1-D) B
    278318          viscosity= -------------------------
    279                                                   2 eps_eff ^[(n-1)/n]
    280 
    281           where B the flow law parameter, eps_eff is the effective strain rate and n the flow law exponent.
     319                                                  2 E^[1/n] eps_eff ^[(n-1)/n]
     320
     321          where B the flow law parameter, eps_eff is the effective strain rate, n the flow law exponent,
     322          and E is the enhancement factor.
    282323
    283324          If eps_eff = 0 , it means this is the first time SystemMatrices is being run, and we
    … …  
    289330
    290331        /*Intermediary: */
    291         IssmDouble B,D=0.,n;
     332        IssmDouble B,D=0.,E=1.,n;
    292333
    293334        /*Get B and n*/
    … …  
    298339                _assert_(D>=0. && D<1.);
    299340        }
     341        if(this->isenhanced){
     342                E=GetEbar();
     343                _assert_(E>0.);
     344        }
    300345
    301346        if (n==1.){
    302                 /*Linear Viscous behavior (Newtonian fluid) viscosity=B/2: */
    303                 viscosity=(1.-D)*B/2.;
     347                /*Linear Viscous behavior (Newtonian fluid) viscosity=B/2E: */
     348                viscosity=(1.-D)*B/(2.*E);
    304349        }
    305350        else{
    … …  
    312357
    313358                else{
    314                         viscosity=(1.-D)*B/(2.*pow(eps_eff,(n-1.)/n));
     359                        viscosity=(1.-D)*B/(2.*pow(E,1./n)*pow(eps_eff,(n-1.)/n));
    315360                }
    316361        }
    … …  
    475520
    476521        /*Intermediary: */
    477         IssmDouble D=0.,n;
     522        IssmDouble D=0.,E=1.,n;
    478523
    479524        /*Get B and n*/
    … …  
    483528                _assert_(D>=0. && D<1.);
    484529        }
     530        if(this->isenhanced){
     531                E=GetE();
     532                _assert_(E>0.);
     533        }
    485534
    486535        if(n==1.){
    487                 /*Linear Viscous behavior (Newtonian fluid) dmudB=B/2: */
    488                 dmudB=(1.-D)/2.;
     536                /*Linear Viscous behavior (Newtonian fluid) dmudB=B/2E: */
     537                dmudB=(1.-D)/(2.*E);
    489538        }
    490539        else{
    491540                if(eps_eff==0.) dmudB = 0.;
    492                 else            dmudB = (1.-D)/(2.*pow(eps_eff,(n-1.)/n));
     541                else            dmudB = (1.-D)/(2.*pow(E,1./n)*pow(eps_eff,(n-1.)/n));
    493542        }
    494543
    … …  
    503552
    504553        /*Intermediary: */
    505         IssmDouble n,B;
     554        IssmDouble n,B,E=1.;
    506555
    507556        /*Get B and n*/
    … …  
    509558        B=GetBbar();
    510559        _assert_(this->isdamaged);
     560        if(this->isenhanced){
     561                E=GetE();
     562                _assert_(E>0.);
     563        }
    511564
    512565        if(n==1.){
    513                 /*Linear Viscous behavior (Newtonian fluid) dmudB=B/2: */
    514                 dmudD=-B/2.;
     566                /*Linear Viscous behavior (Newtonian fluid) dmudB=B/2E: */
     567                dmudD=-B/(2.*E);
    515568        }
    516569        else{
    517570                if(eps_eff==0.) dmudD = 0.;
    518                 else            dmudD = -B/(2.*pow(eps_eff,(n-1.)/n));
     571                else            dmudD = -B/(2.*pow(E,1./n)*pow(eps_eff,(n-1.)/n));
    519572        }
    520573
    … …  
    636689        IssmDouble epsilon3d[5];/* epsilon=[exx,eyy,exy,exz,eyz];*/
    637690        IssmDouble epsilon2d[2];/* epsilon=[exx,exy];            */
    638         IssmDouble eps_eff;
     691        IssmDouble eps_eff,E=1.0;
    639692
    640693        if(dim==3){

  • issm/trunk-jpl/src/c/classes/Materials/Matice.h

    r20827 r21389  
    2727                int      mid;
    2828                bool     isdamaged;
     29                bool     isenhanced;
    2930                Hook    *helement;
    3031                Element *element;
    … …  
    7071                IssmDouble GetD();
    7172                IssmDouble GetDbar();
     73                IssmDouble GetE();
     74                IssmDouble GetEbar();
    7275                IssmDouble GetN();
    7376                bool       IsDamage();
     77                bool       IsEnhanced();
    7478                void       ResetHooks();
    7579                void       SetCurrentConfiguration(Elements* elements,Loads* loads,Nodes* nodes,Vertices* vertices,Materials* materials,Parameters* parameters);

  • issm/trunk-jpl/src/c/classes/Materials/Matpar.cpp

    r20827 r21389  
    8282                case MaticeEnum:
    8383                case MatdamageiceEnum:
     84                case MatenhancediceEnum:
    8485                case MatestarEnum:
    8586                        iomodel->FindConstant(&this->rho_ice,"md.materials.rho_ice");

  • issm/trunk-jpl/src/c/modules/ModelProcessorx/CreateElementsVerticesAndMaterials.cpp

    r21381 r21389  
    5656                                case 2:
    5757                                        elements->InputDuplicate(MaterialsRheologyBEnum,MaterialsRheologyBbarEnum);
     58                                        break;
     59                                case 3:
     60                                        break;
     61                                default:
     62                                        _error_("Mesh not supported yet");
     63                        }
     64                        break;
     65                case MatenhancediceEnum:
     66                        iomodel->FetchDataToInput(elements,"md.materials.rheology_B",MaterialsRheologyBEnum);
     67                        iomodel->FetchDataToInput(elements,"md.materials.rheology_n",MaterialsRheologyNEnum);
     68                        iomodel->FetchDataToInput(elements,"md.materials.rheology_E",MaterialsRheologyEEnum);
     69                        for (i=0;i<iomodel->numberofelements;i++) if(iomodel->my_elements[i]) materials->AddObject(new Matice(i+1,i,iomodel));
     70                        switch(iomodel->domaindim){
     71                                case 2:
     72                                        elements->InputDuplicate(MaterialsRheologyBEnum,MaterialsRheologyBbarEnum);
     73                                        elements->InputDuplicate(MaterialsRheologyEEnum,MaterialsRheologyEbarEnum);
    5874                                        break;
    5975                                case 3:

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

    r21381 r21389  
    214214        MaterialsRheologyLawEnum,
    215215        MaterialsRheologyNEnum,
     216        MaterialsRheologyEEnum,
     217        MaterialsRheologyEbarEnum,
    216218        MaterialsRheologyEcEnum,
    217219        MaterialsRheologyEcbarEnum,
    … …  
    828830        MaticeEnum,
    829831        MatdamageiceEnum,
     832        MatenhancediceEnum,
    830833        MatestarEnum,
    831834        MatparEnum,

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

    r21385 r21389  
    220220                case MaterialsRheologyLawEnum : return "MaterialsRheologyLaw";
    221221                case MaterialsRheologyNEnum : return "MaterialsRheologyN";
     222                case MaterialsRheologyEEnum : return "MaterialsRheologyE";
     223                case MaterialsRheologyEbarEnum : return "MaterialsRheologyEbar";
    222224                case MaterialsRheologyEcEnum : return "MaterialsRheologyEc";
    223225                case MaterialsRheologyEcbarEnum : return "MaterialsRheologyEcbar";
    … …  
    803805                case MaticeEnum : return "Matice";
    804806                case MatdamageiceEnum : return "Matdamageice";
     807                case MatenhancediceEnum : return "Matenhancedice";
    805808                case MatestarEnum : return "Matestar";
    806809                case MatparEnum : return "Matpar";

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

    r21385 r21389  
    223223              else if (strcmp(name,"MaterialsRheologyLaw")==0) return MaterialsRheologyLawEnum;
    224224              else if (strcmp(name,"MaterialsRheologyN")==0) return MaterialsRheologyNEnum;
     225              else if (strcmp(name,"MaterialsRheologyE")==0) return MaterialsRheologyEEnum;
     226              else if (strcmp(name,"MaterialsRheologyEbar")==0) return MaterialsRheologyEbarEnum;
    225227              else if (strcmp(name,"MaterialsRheologyEc")==0) return MaterialsRheologyEcEnum;
    226228              else if (strcmp(name,"MaterialsRheologyEcbar")==0) return MaterialsRheologyEcbarEnum;
    … …  
    258260              else if (strcmp(name,"CalvinglevermannMeltingrate")==0) return CalvinglevermannMeltingrateEnum;
    259261              else if (strcmp(name,"CalvingdevCoeff")==0) return CalvingdevCoeffEnum;
    260               else if (strcmp(name,"Calvingratex")==0) return CalvingratexEnum;
    261               else if (strcmp(name,"Calvingratey")==0) return CalvingrateyEnum;
    262262         else stage=3;
    263263   }
    264264   if(stage==3){
    265               if (strcmp(name,"CalvingratexAverage")==0) return CalvingratexAverageEnum;
     265              if (strcmp(name,"Calvingratex")==0) return CalvingratexEnum;
     266              else if (strcmp(name,"Calvingratey")==0) return CalvingrateyEnum;
     267              else if (strcmp(name,"CalvingratexAverage")==0) return CalvingratexAverageEnum;
    266268              else if (strcmp(name,"CalvingrateyAverage")==0) return CalvingrateyAverageEnum;
    267269              else if (strcmp(name,"StrainRateparallel")==0) return StrainRateparallelEnum;
    … …  
    381383              else if (strcmp(name,"SmbSizeini")==0) return SmbSizeiniEnum;
    382384              else if (strcmp(name,"SMBforcing")==0) return SMBforcingEnum;
    383               else if (strcmp(name,"SmbMassBalance")==0) return SmbMassBalanceEnum;
    384               else if (strcmp(name,"SMBgemb")==0) return SMBgembEnum;
    385385         else stage=4;
    386386   }
    387387   if(stage==4){
    388               if (strcmp(name,"SmbInitDensityScaling")==0) return SmbInitDensityScalingEnum;
     388              if (strcmp(name,"SmbMassBalance")==0) return SmbMassBalanceEnum;
     389              else if (strcmp(name,"SMBgemb")==0) return SMBgembEnum;
     390              else if (strcmp(name,"SmbInitDensityScaling")==0) return SmbInitDensityScalingEnum;
    389391              else if (strcmp(name,"SmbTa")==0) return SmbTaEnum;
    390392              else if (strcmp(name,"SmbV")==0) return SmbVEnum;
    … …  
    504506              else if (strcmp(name,"VxAverage")==0) return VxAverageEnum;
    505507              else if (strcmp(name,"Vx")==0) return VxEnum;
    506               else if (strcmp(name,"VxPicard")==0) return VxPicardEnum;
    507               else if (strcmp(name,"VyAverage")==0) return VyAverageEnum;
    508508         else stage=5;
    509509   }
    510510   if(stage==5){
    511               if (strcmp(name,"Vy")==0) return VyEnum;
     511              if (strcmp(name,"VxPicard")==0) return VxPicardEnum;
     512              else if (strcmp(name,"VyAverage")==0) return VyAverageEnum;
     513              else if (strcmp(name,"Vy")==0) return VyEnum;
    512514              else if (strcmp(name,"VyPicard")==0) return VyPicardEnum;
    513515              else if (strcmp(name,"Vz")==0) return VzEnum;
    … …  
    627629              else if (strcmp(name,"Outputdefinition46")==0) return Outputdefinition46Enum;
    628630              else if (strcmp(name,"Outputdefinition47")==0) return Outputdefinition47Enum;
    629               else if (strcmp(name,"Outputdefinition48")==0) return Outputdefinition48Enum;
    630               else if (strcmp(name,"Outputdefinition49")==0) return Outputdefinition49Enum;
    631631         else stage=6;
    632632   }
    633633   if(stage==6){
    634               if (strcmp(name,"Outputdefinition50")==0) return Outputdefinition50Enum;
     634              if (strcmp(name,"Outputdefinition48")==0) return Outputdefinition48Enum;
     635              else if (strcmp(name,"Outputdefinition49")==0) return Outputdefinition49Enum;
     636              else if (strcmp(name,"Outputdefinition50")==0) return Outputdefinition50Enum;
    635637              else if (strcmp(name,"Outputdefinition51")==0) return Outputdefinition51Enum;
    636638              else if (strcmp(name,"Outputdefinition52")==0) return Outputdefinition52Enum;
    … …  
    750752              else if (strcmp(name,"LliboutryDuval")==0) return LliboutryDuvalEnum;
    751753              else if (strcmp(name,"ExtrapolationVariable")==0) return ExtrapolationVariableEnum;
    752               else if (strcmp(name,"IceMaskNodeActivation")==0) return IceMaskNodeActivationEnum;
    753               else if (strcmp(name,"LevelsetfunctionSlopeX")==0) return LevelsetfunctionSlopeXEnum;
    754754         else stage=7;
    755755   }
    756756   if(stage==7){
    757               if (strcmp(name,"LevelsetfunctionSlopeY")==0) return LevelsetfunctionSlopeYEnum;
     757              if (strcmp(name,"IceMaskNodeActivation")==0) return IceMaskNodeActivationEnum;
     758              else if (strcmp(name,"LevelsetfunctionSlopeX")==0) return LevelsetfunctionSlopeXEnum;
     759              else if (strcmp(name,"LevelsetfunctionSlopeY")==0) return LevelsetfunctionSlopeYEnum;
    758760              else if (strcmp(name,"LevelsetfunctionPicard")==0) return LevelsetfunctionPicardEnum;
    759761              else if (strcmp(name,"LevelsetReinitFrequency")==0) return LevelsetReinitFrequencyEnum;
    … …  
    821823              else if (strcmp(name,"Matice")==0) return MaticeEnum;
    822824              else if (strcmp(name,"Matdamageice")==0) return MatdamageiceEnum;
     825              else if (strcmp(name,"Matenhancedice")==0) return MatenhancediceEnum;
    823826              else if (strcmp(name,"Matestar")==0) return MatestarEnum;
    824827              else if (strcmp(name,"Matpar")==0) return MatparEnum;
    … …  
    872875              else if (strcmp(name,"BalancevelocitySolution")==0) return BalancevelocitySolutionEnum;
    873876              else if (strcmp(name,"L2ProjectionEPLAnalysis")==0) return L2ProjectionEPLAnalysisEnum;
    874               else if (strcmp(name,"L2ProjectionBaseAnalysis")==0) return L2ProjectionBaseAnalysisEnum;
    875               else if (strcmp(name,"BedSlopeSolution")==0) return BedSlopeSolutionEnum;
    876               else if (strcmp(name,"DamageEvolutionSolution")==0) return DamageEvolutionSolutionEnum;
    877877         else stage=8;
    878878   }
    879879   if(stage==8){
    880               if (strcmp(name,"DamageEvolutionAnalysis")==0) return DamageEvolutionAnalysisEnum;
     880              if (strcmp(name,"L2ProjectionBaseAnalysis")==0) return L2ProjectionBaseAnalysisEnum;
     881              else if (strcmp(name,"BedSlopeSolution")==0) return BedSlopeSolutionEnum;
     882              else if (strcmp(name,"DamageEvolutionSolution")==0) return DamageEvolutionSolutionEnum;
     883              else if (strcmp(name,"DamageEvolutionAnalysis")==0) return DamageEvolutionAnalysisEnum;
    881884              else if (strcmp(name,"StressbalanceAnalysis")==0) return StressbalanceAnalysisEnum;
    882885              else if (strcmp(name,"StressbalanceSIAAnalysis")==0) return StressbalanceSIAAnalysisEnum;

  • issm/trunk-jpl/src/c/shared/io/Marshalling/IoCodeConversions.cpp

    r21049 r21389  
    8686                case 2: return MatestarEnum;
    8787                case 3: return MaticeEnum;
     88                case 4: return MatenhancediceEnum;
    8889                default: _error_("Marshalled materials code \""<<enum_in<<"\" not supported yet");
    8990        }
Note: See TracChangeset for help on using the changeset viewer.