source: issm/oecreview/Archive/21337-21723/ISSM-21550-21551.diff@ 21726

Last change on this file since 21726 was 21726, checked in by Mathieu Morlighem, 8 years ago

CHG added Archive/21337-21723
File size: 18.8 KB

  • ../trunk-jpl/src/m/classes/frictionjosh.m

     
    66classdef frictionjosh
    77        properties (SetAccess=public)
    88                coefficient = NaN;
    9                 temperature = NaN;
     9                pressure_adjusted_temperature = NaN;
    1010        end
    1111        methods
    1212                function self = extrude(self,md) % {{{
    1313                        self.coefficient=project3d(md,'vector',self.coefficient,'type','node','layer',1);
    14                         self.temperature=project3d(md,'vector',self.temperature,'type','node','layer',1);
     14                        self.pressure_adjusted_temperature=project3d(md,'vector',self.pressure_adjusted_temperature,'type','node','layer',1);
    1515                end % }}}
    1616                function self = frictionjosh(varargin) % {{{
    1717                        switch nargin
    … …  
    3131                        %Early return
    3232                        if ~ismember('StressbalanceAnalysis',analyses) & ~ismember('ThermalAnalysis',analyses), return; end
    3333
    34                         md = checkfield(md,'fieldname','friction.coefficient','timeseries',1,'NaN',1,'Inf',1);
    35                         md = checkfield(md,'fieldname','friction.temperature','timeseries',1,'NaN',1,'Inf',1);
     34                        md = checkfield(md,'fieldname','friction.coefficient','NaN',1,'Inf',1);
     35                        md = checkfield(md,'fieldname','friction.pressure_adjusted_temperature','NaN',1,'Inf',1);
    3636
    3737                        %Check that temperature is provided
    3838                        md = checkfield(md,'fieldname','initialization.temperature','NaN',1,'Inf',1,'size',[md.mesh.numberofvertices 1]);
    … …  
    4040                function disp(self) % {{{
    4141                        disp(sprintf('Basal shear stress parameters: tau_b = coefficient^2 * Neff ^r * |u_b|^(s-1) * u_b * 1/f(T)\n(effective stress Neff=rho_ice*g*thickness+rho_water*g*bed, r=q/p and s=1/p)'));
    4242                        fielddisplay(self,'coefficient','friction coefficient [SI]');
    43                         fielddisplay(self,'temperature','friction temperature [K]');
     43                        fielddisplay(self,'pressure_adjusted_temperature','friction pressure_adjusted_temperature (T - Tpmp) [K]');
    4444                end % }}}
    4545                function marshall(self,prefix,md,fid) % {{{
    4646
    4747                        WriteData(fid,prefix,'name','md.friction.law','data',9,'format','Integer');
    4848                        WriteData(fid,prefix,'class','friction','object',self,'fieldname','coefficient','format','DoubleMat','mattype',1,'timeserieslength',md.mesh.numberofvertices+1,'yts',md.constants.yts);
    49                         WriteData(fid,prefix,'class','friction','object',self,'fieldname','temperature','format','DoubleMat','mattype',1,'timeserieslength',md.mesh.numberofvertices+1,'yts',md.constants.yts);
     49                        WriteData(fid,prefix,'class','friction','object',self,'fieldname','pressure_adjusted_temperature','format','DoubleMat','mattype',1,'timeserieslength',md.mesh.numberofvertices+1,'yts',md.constants.yts);
    5050                end % }}}
    5151        end
    5252end

  • ../trunk-jpl/src/c/classes/Elements/Element.cpp

     
    20322032                                yearlytemperatures[2] = s[2];
    20332033                                this->inputs->AddInput(new PentaInput(TemperatureEnum,&yearlytemperatures[0],P1Enum));
    20342034
    2035                                 /*Convert that to enthalpy for the enthalpy model*/
    2036                                 IssmDouble enthalpy[6];
    2037                                 GetInputListOnVertices(&enthalpy[0],EnthalpyEnum);
    2038                                 ThermalToEnthalpy(&enthalpy[3],yearlytemperatures[3],0.,0.);
    2039                                 ThermalToEnthalpy(&enthalpy[4],yearlytemperatures[3],0.,0.);
    2040                                 ThermalToEnthalpy(&enthalpy[5],yearlytemperatures[3],0.,0.);
    2041                                 this->inputs->AddInput(new PentaInput(EnthalpyEnum,&enthalpy[0],P1Enum));
     2035                                bool isenthalpy;
     2036                                this->parameters->FindParam(&isenthalpy,ThermalIsenthalpyEnum);
     2037                                if(isenthalpy){
     2038                                        /*Convert that to enthalpy for the enthalpy model*/
     2039                                        IssmDouble enthalpy[6];
     2040                                        GetInputListOnVertices(&enthalpy[0],EnthalpyEnum);
     2041                                        ThermalToEnthalpy(&enthalpy[3],yearlytemperatures[3],0.,0.);
     2042                                        ThermalToEnthalpy(&enthalpy[4],yearlytemperatures[3],0.,0.);
     2043                                        ThermalToEnthalpy(&enthalpy[5],yearlytemperatures[3],0.,0.);
     2044                                        this->inputs->AddInput(new PentaInput(EnthalpyEnum,&enthalpy[0],P1Enum));
     2045                                }
    20422046                        }
    20432047                        this->inputs->AddInput(new PentaInput(SmbMassBalanceEnum,&agd[0],P1Enum));
    20442048                        this->inputs->AddInput(new PentaInput(TemperaturePDDEnum,&yearlytemperatures[0],P1Enum));

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

     
    246246                case 8:
    247247                        GetAlpha2Sommers(palpha2,gauss);
    248248                        break;
     249                case 9:
     250                        GetAlpha2Josh(palpha2,gauss);
     251                        break;
    249252          default:
    250253                        _error_("Friction law "<< this->law <<" not supported");
    251254        }
    … …  
    462465        /*Assign output pointers:*/
    463466        *palpha2=alpha2;
    464467}/*}}}*/
     468void Friction::GetAlpha2Josh(IssmDouble* palpha2, Gauss* gauss){/*{{{*/
     469        /*Here, we want to parameterize the friction as a function of temperature
     470         *
     471         * alpha2 = alpha2_viscous * 1/f(T)
     472         *
     473         * where f(T) = exp((T-Tpmp)/gamma)
     474         */
     475
     476        /*Intermediaries: */
     477        IssmDouble  T,Tpmp,deltaT,deltaTref,pressure;
     478        IssmDouble  alpha2,time,gamma;
     479        const IssmDouble yts = 365*24*3600.;
     480
     481        /*Get viscous part*/
     482        this->GetAlpha2Viscous(&alpha2,gauss);
     483
     484        /*Get delta Refs*/
     485        element->GetInputValue(&deltaTref,gauss,FrictionPressureAdjustedTemperatureEnum);
     486
     487        /*Compute delta T*/
     488        element->GetInputValue(&T,gauss,TemperatureEnum);
     489        element->GetInputValue(&pressure,gauss,PressureEnum);
     490        Tpmp = element->TMeltingPoint(pressure);
     491        deltaT = T-Tpmp;
     492
     493        /*Compute gamma*/
     494        element->parameters->FindParam(&time,TimeEnum);
     495        if(time<25e3*yts){
     496                gamma = 10.;
     497        }
     498        else{
     499                gamma = 5.;
     500        }
     501        gamma = 5.;
     502
     503        /*Compute scaling parameter*/
     504        alpha2 = alpha2 * exp((deltaTref - deltaT)/(2*gamma));
     505
     506        /*Assign output pointers:*/
     507        *palpha2=alpha2;
     508}/*}}}*/
    465509void Friction::GetAlpha2Viscous(IssmDouble* palpha2, Gauss* gauss){/*{{{*/
    466510
    467511        /*This routine calculates the basal friction coefficient

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

     
    3535                void  GetAlpha2(IssmDouble* palpha2,Gauss* gauss);
    3636                void  GetAlpha2Coulomb(IssmDouble* palpha2,Gauss* gauss);
    3737                void  GetAlpha2Hydro(IssmDouble* palpha2,Gauss* gauss);
     38                void  GetAlpha2Josh(IssmDouble* palpha2,Gauss* gauss);
    3839                void  GetAlpha2Sommers(IssmDouble* palpha2,Gauss* gauss);
    3940                void  GetAlpha2Temp(IssmDouble* palpha2,Gauss* gauss);
    4041                void  GetAlpha2Viscous(IssmDouble* palpha2,Gauss* gauss);

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

     
    822822                        iomodel->FetchDataToInput(elements,"md.friction.p",FrictionPEnum);
    823823                        iomodel->FetchDataToInput(elements,"md.friction.q",FrictionQEnum);
    824824                        break;
     825                case 9:
     826                        iomodel->FetchDataToInput(elements,"md.friction.coefficient",FrictionCoefficientEnum);
     827                        iomodel->FetchDataToInput(elements,"md.friction.pressure_adjusted_temperature",FrictionPressureAdjustedTemperatureEnum);
     828                        break;
    825829                default:
    826830                        _error_("not supported");
    827831        }

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

     
    108108        FrictionAsEnum,
    109109        FrictionCoefficientEnum,
    110110        FrictionCoefficientcoulombEnum,
     111        FrictionPressureAdjustedTemperatureEnum,
    111112        FrictionPEnum,
    112113        FrictionQEnum,
    113114        FrictionMEnum,
    … …  
    273274        MaterialsMantleDensityEnum,
    274275        MaterialsEarthDensityEnum,
    275276        MeshAverageVertexConnectivityEnum,
     277        MeshXEnum,
     278        MeshYEnum,
     279        MeshZEnum,
    276280        MeshElementsEnum,
    277281        MeshNumberofelementsEnum,
    278282        MeshNumberoflayersEnum,

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

     
    114114                case FrictionAsEnum : return "FrictionAs";
    115115                case FrictionCoefficientEnum : return "FrictionCoefficient";
    116116                case FrictionCoefficientcoulombEnum : return "FrictionCoefficientcoulomb";
     117                case FrictionPressureAdjustedTemperatureEnum : return "FrictionPressureAdjustedTemperature";
    117118                case FrictionPEnum : return "FrictionP";
    118119                case FrictionQEnum : return "FrictionQ";
    119120                case FrictionMEnum : return "FrictionM";
    … …  
    279280                case MaterialsMantleDensityEnum : return "MaterialsMantleDensity";
    280281                case MaterialsEarthDensityEnum : return "MaterialsEarthDensity";
    281282                case MeshAverageVertexConnectivityEnum : return "MeshAverageVertexConnectivity";
     283                case MeshXEnum : return "MeshX";
     284                case MeshYEnum : return "MeshY";
     285                case MeshZEnum : return "MeshZ";
    282286                case MeshElementsEnum : return "MeshElements";
    283287                case MeshNumberofelementsEnum : return "MeshNumberofelements";
    284288                case MeshNumberoflayersEnum : return "MeshNumberoflayers";

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

     
    114114              else if (strcmp(name,"FrictionAs")==0) return FrictionAsEnum;
    115115              else if (strcmp(name,"FrictionCoefficient")==0) return FrictionCoefficientEnum;
    116116              else if (strcmp(name,"FrictionCoefficientcoulomb")==0) return FrictionCoefficientcoulombEnum;
     117              else if (strcmp(name,"FrictionPressureAdjustedTemperature")==0) return FrictionPressureAdjustedTemperatureEnum;
    117118              else if (strcmp(name,"FrictionP")==0) return FrictionPEnum;
    118119              else if (strcmp(name,"FrictionQ")==0) return FrictionQEnum;
    119120              else if (strcmp(name,"FrictionM")==0) return FrictionMEnum;
    … …  
    135136              else if (strcmp(name,"EffectivePressure")==0) return EffectivePressureEnum;
    136137              else if (strcmp(name,"EplHead")==0) return EplHeadEnum;
    137138              else if (strcmp(name,"EplHeadOld")==0) return EplHeadOldEnum;
    138               else if (strcmp(name,"EplHeadSlopeX")==0) return EplHeadSlopeXEnum;
    139139         else stage=2;
    140140   }
    141141   if(stage==2){
    142               if (strcmp(name,"EplHeadSlopeY")==0) return EplHeadSlopeYEnum;
     142              if (strcmp(name,"EplHeadSlopeX")==0) return EplHeadSlopeXEnum;
     143              else if (strcmp(name,"EplHeadSlopeY")==0) return EplHeadSlopeYEnum;
    143144              else if (strcmp(name,"EplZigZagCounter")==0) return EplZigZagCounterEnum;
    144145              else if (strcmp(name,"HydrologydcMaxIter")==0) return HydrologydcMaxIterEnum;
    145146              else if (strcmp(name,"HydrologydcRelTol")==0) return HydrologydcRelTolEnum;
    … …  
    258259              else if (strcmp(name,"CalvingLevermann")==0) return CalvingLevermannEnum;
    259260              else if (strcmp(name,"CalvingDev")==0) return CalvingDevEnum;
    260261              else if (strcmp(name,"CalvingMinthickness")==0) return CalvingMinthicknessEnum;
    261               else if (strcmp(name,"DefaultCalving")==0) return DefaultCalvingEnum;
    262262         else stage=3;
    263263   }
    264264   if(stage==3){
    265               if (strcmp(name,"CalvinglevermannCoeff")==0) return CalvinglevermannCoeffEnum;
     265              if (strcmp(name,"DefaultCalving")==0) return DefaultCalvingEnum;
     266              else if (strcmp(name,"CalvinglevermannCoeff")==0) return CalvinglevermannCoeffEnum;
    266267              else if (strcmp(name,"CalvinglevermannMeltingrate")==0) return CalvinglevermannMeltingrateEnum;
    267268              else if (strcmp(name,"CalvingdevCoeff")==0) return CalvingdevCoeffEnum;
    268269              else if (strcmp(name,"Calvingratex")==0) return CalvingratexEnum;
    … …  
    285286              else if (strcmp(name,"MaterialsMantleDensity")==0) return MaterialsMantleDensityEnum;
    286287              else if (strcmp(name,"MaterialsEarthDensity")==0) return MaterialsEarthDensityEnum;
    287288              else if (strcmp(name,"MeshAverageVertexConnectivity")==0) return MeshAverageVertexConnectivityEnum;
     289              else if (strcmp(name,"MeshX")==0) return MeshXEnum;
     290              else if (strcmp(name,"MeshY")==0) return MeshYEnum;
     291              else if (strcmp(name,"MeshZ")==0) return MeshZEnum;
    288292              else if (strcmp(name,"MeshElements")==0) return MeshElementsEnum;
    289293              else if (strcmp(name,"MeshNumberofelements")==0) return MeshNumberofelementsEnum;
    290294              else if (strcmp(name,"MeshNumberoflayers")==0) return MeshNumberoflayersEnum;
    … …  
    378382              else if (strcmp(name,"SmbDzini")==0) return SmbDziniEnum;
    379383              else if (strcmp(name,"SmbDini")==0) return SmbDiniEnum;
    380384              else if (strcmp(name,"SmbReini")==0) return SmbReiniEnum;
    381               else if (strcmp(name,"SmbGdnini")==0) return SmbGdniniEnum;
     385         else stage=4;
     386   }
     387   if(stage==4){
     388              if (strcmp(name,"SmbGdnini")==0) return SmbGdniniEnum;
    382389              else if (strcmp(name,"SmbGspini")==0) return SmbGspiniEnum;
    383390              else if (strcmp(name,"SmbECini")==0) return SmbECiniEnum;
    384391              else if (strcmp(name,"SmbWini")==0) return SmbWiniEnum;
    385          else stage=4;
    386    }
    387    if(stage==4){
    388               if (strcmp(name,"SmbAini")==0) return SmbAiniEnum;
     392              else if (strcmp(name,"SmbAini")==0) return SmbAiniEnum;
    389393              else if (strcmp(name,"SmbTini")==0) return SmbTiniEnum;
    390394              else if (strcmp(name,"SmbSizeini")==0) return SmbSizeiniEnum;
    391395              else if (strcmp(name,"SMBforcing")==0) return SMBforcingEnum;
    … …  
    501505              else if (strcmp(name,"SurfaceAverageVelMisfit")==0) return SurfaceAverageVelMisfitEnum;
    502506              else if (strcmp(name,"SurfaceLogVelMisfit")==0) return SurfaceLogVelMisfitEnum;
    503507              else if (strcmp(name,"SurfaceLogVxVyMisfit")==0) return SurfaceLogVxVyMisfitEnum;
    504               else if (strcmp(name,"SurfaceRelVelMisfit")==0) return SurfaceRelVelMisfitEnum;
     508         else stage=5;
     509   }
     510   if(stage==5){
     511              if (strcmp(name,"SurfaceRelVelMisfit")==0) return SurfaceRelVelMisfitEnum;
    505512              else if (strcmp(name,"SurfaceSlopeX")==0) return SurfaceSlopeXEnum;
    506513              else if (strcmp(name,"SurfaceSlopeY")==0) return SurfaceSlopeYEnum;
    507514              else if (strcmp(name,"Temperature")==0) return TemperatureEnum;
    508          else stage=5;
    509    }
    510    if(stage==5){
    511               if (strcmp(name,"TemperaturePicard")==0) return TemperaturePicardEnum;
     515              else if (strcmp(name,"TemperaturePicard")==0) return TemperaturePicardEnum;
    512516              else if (strcmp(name,"TemperaturePDD")==0) return TemperaturePDDEnum;
    513517              else if (strcmp(name,"ThicknessAbsMisfit")==0) return ThicknessAbsMisfitEnum;
    514518              else if (strcmp(name,"SurfaceAbsMisfit")==0) return SurfaceAbsMisfitEnum;
    … …  
    624628              else if (strcmp(name,"Outputdefinition33")==0) return Outputdefinition33Enum;
    625629              else if (strcmp(name,"Outputdefinition34")==0) return Outputdefinition34Enum;
    626630              else if (strcmp(name,"Outputdefinition35")==0) return Outputdefinition35Enum;
    627               else if (strcmp(name,"Outputdefinition36")==0) return Outputdefinition36Enum;
     631         else stage=6;
     632   }
     633   if(stage==6){
     634              if (strcmp(name,"Outputdefinition36")==0) return Outputdefinition36Enum;
    628635              else if (strcmp(name,"Outputdefinition37")==0) return Outputdefinition37Enum;
    629636              else if (strcmp(name,"Outputdefinition38")==0) return Outputdefinition38Enum;
    630637              else if (strcmp(name,"Outputdefinition39")==0) return Outputdefinition39Enum;
    631          else stage=6;
    632    }
    633    if(stage==6){
    634               if (strcmp(name,"Outputdefinition40")==0) return Outputdefinition40Enum;
     638              else if (strcmp(name,"Outputdefinition40")==0) return Outputdefinition40Enum;
    635639              else if (strcmp(name,"Outputdefinition41")==0) return Outputdefinition41Enum;
    636640              else if (strcmp(name,"Outputdefinition42")==0) return Outputdefinition42Enum;
    637641              else if (strcmp(name,"Outputdefinition43")==0) return Outputdefinition43Enum;
    … …  
    747751              else if (strcmp(name,"Sset")==0) return SsetEnum;
    748752              else if (strcmp(name,"Dense")==0) return DenseEnum;
    749753              else if (strcmp(name,"MpiDense")==0) return MpiDenseEnum;
    750               else if (strcmp(name,"MpiSparse")==0) return MpiSparseEnum;
     754         else stage=7;
     755   }
     756   if(stage==7){
     757              if (strcmp(name,"MpiSparse")==0) return MpiSparseEnum;
    751758              else if (strcmp(name,"Seq")==0) return SeqEnum;
    752759              else if (strcmp(name,"Mpi")==0) return MpiEnum;
    753760              else if (strcmp(name,"Mumps")==0) return MumpsEnum;
    754          else stage=7;
    755    }
    756    if(stage==7){
    757               if (strcmp(name,"Gsl")==0) return GslEnum;
     761              else if (strcmp(name,"Gsl")==0) return GslEnum;
    758762              else if (strcmp(name,"Cuffey")==0) return CuffeyEnum;
    759763              else if (strcmp(name,"BuddJacka")==0) return BuddJackaEnum;
    760764              else if (strcmp(name,"CuffeyTemperate")==0) return CuffeyTemperateEnum;
    … …  
    870874              else if (strcmp(name,"GenericOption")==0) return GenericOptionEnum;
    871875              else if (strcmp(name,"OptionCell")==0) return OptionCellEnum;
    872876              else if (strcmp(name,"OptionStruct")==0) return OptionStructEnum;
    873               else if (strcmp(name,"AdjointBalancethicknessAnalysis")==0) return AdjointBalancethicknessAnalysisEnum;
     877         else stage=8;
     878   }
     879   if(stage==8){
     880              if (strcmp(name,"AdjointBalancethicknessAnalysis")==0) return AdjointBalancethicknessAnalysisEnum;
    874881              else if (strcmp(name,"AdjointBalancethickness2Analysis")==0) return AdjointBalancethickness2AnalysisEnum;
    875882              else if (strcmp(name,"AdjointHorizAnalysis")==0) return AdjointHorizAnalysisEnum;
    876883              else if (strcmp(name,"DefaultAnalysis")==0) return DefaultAnalysisEnum;
    877          else stage=8;
    878    }
    879    if(stage==8){
    880               if (strcmp(name,"BalancethicknessAnalysis")==0) return BalancethicknessAnalysisEnum;
     884              else if (strcmp(name,"BalancethicknessAnalysis")==0) return BalancethicknessAnalysisEnum;
    881885              else if (strcmp(name,"BalancethicknessSolution")==0) return BalancethicknessSolutionEnum;
    882886              else if (strcmp(name,"Balancethickness2Analysis")==0) return Balancethickness2AnalysisEnum;
    883887              else if (strcmp(name,"Balancethickness2Solution")==0) return Balancethickness2SolutionEnum;
    … …  
    993997              else if (strcmp(name,"Melange")==0) return MelangeEnum;
    994998              else if (strcmp(name,"Water")==0) return WaterEnum;
    995999              else if (strcmp(name,"DataSet")==0) return DataSetEnum;
    996               else if (strcmp(name,"Constraints")==0) return ConstraintsEnum;
     1000         else stage=9;
     1001   }
     1002   if(stage==9){
     1003              if (strcmp(name,"Constraints")==0) return ConstraintsEnum;
    9971004              else if (strcmp(name,"Loads")==0) return LoadsEnum;
    9981005              else if (strcmp(name,"Materials")==0) return MaterialsEnum;
    9991006              else if (strcmp(name,"Nodes")==0) return NodesEnum;
    1000          else stage=9;
    1001    }
    1002    if(stage==9){
    1003               if (strcmp(name,"Contours")==0) return ContoursEnum;
     1007              else if (strcmp(name,"Contours")==0) return ContoursEnum;
    10041008              else if (strcmp(name,"Parameters")==0) return ParametersEnum;
    10051009              else if (strcmp(name,"Vertices")==0) return VerticesEnum;
    10061010              else if (strcmp(name,"Results")==0) return ResultsEnum;
Note: See TracBrowser for help on using the repository browser.