Changeset 27161


Ignore:
Timestamp:
07/18/22 12:42:48 (3 years ago)
Author:
vverjans
Message:

CHG: allowing stochastic water pressure for frictionschoof

Location:
issm/trunk-jpl/src
Files:
9 edited

Legend:

Unmodified
Added
Removed

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

    r26676 r27161  
    749749         *               C |u_b|^(m-1)
    750750         * alpha2= __________________________
    751          *          (1+(C/(Cmax N))^1/m |u_b| )^m
     751         *          (1+(C/(Cmax Neff))^1/m |u_b| )^m
    752752         *
    753753         * */
    … …  
    764764        C = coeff*coeff;
    765765
    766         /*Get effective pressure and velocity magnitude*/
    767         IssmDouble N  = EffectivePressure(gauss);
     766        /*Get effective pressure*/
     767        bool ispwStochastic;
     768        IssmDouble Neff;
     769        element->parameters->FindParam(&ispwStochastic,StochasticForcingIsWaterPressureEnum);
     770        if(ispwStochastic){
     771                /*Retrieve stochastic water pressure and compute ice pressure*/
     772                IssmDouble p_ice,p_water,Neff_limit;
     773                element->GetInputValue(&p_water,gauss,FrictionSchoofWaterPressureEnum);
     774                element->parameters->FindParam(&Neff_limit,FrictionEffectivePressureLimitEnum);
     775                p_ice = IcePressure(gauss);
     776                Neff  = max(Neff_limit*p_ice, p_ice - p_water);
     777        }       
     778        else{
     779                /*Compute effective pressure directly*/
     780                Neff = EffectivePressure(gauss);
     781        }
     782
     783        /*Get velocity magnitude*/
    768784        IssmDouble ub = VelMag(gauss);
    769785
    770786        /*Compute alpha^2*/
    771         if((ub<1e-10) ||(N==0.0)){
     787        if((ub<1e-10) ||(Neff==0.0)){
    772788                alpha2 = 0.;
    773789        }
    774790        else{
    775                 alpha2= (C*pow(ub,m-1.)) / pow(1.+  pow(C/(Cmax*N),1./m)*ub,m);
     791                alpha2= (C*pow(ub,m-1.)) / pow(1.+  pow(C/(Cmax*Neff),1./m)*ub,m);
    776792        }
    777793

  • issm/trunk-jpl/src/c/modules/StochasticForcingx/StochasticForcingx.cpp

    r27160 r27161  
    214214                                        }
    215215                                        break;
    216                                 default:
     216                                case FrictionSchoofWaterPressureEnum:
     217                                        /*Specify that WaterPressure is stochastic*/
     218                                        femmodel->parameters->SetParam(true,StochasticForcingIsWaterPressureEnum);
     219                                        for(Object* &object:femmodel->elements->objects){
     220                  Element* element = xDynamicCast<Element*>(object);
     221                  int numvertices  = element->GetNumberOfVertices();
     222                  IssmDouble p_water_deterministic[numvertices];
     223                  IssmDouble p_water[numvertices];
     224                                                element->GetInputValue(&dimensionid,StochasticForcingDefaultIdEnum);
     225                                                Gauss* gauss=element->NewGauss();
     226                                                Friction* friction = new Friction(element);
     227                                                for(int i=0;i<numvertices;i++){
     228                                                        gauss->GaussVertex(i);
     229                                                        p_water_deterministic[i] = friction->SubglacialWaterPressure(gauss);
     230                                                        p_water[i]               = p_water_deterministic[i] + noisefield[dimensionid]; //make sure positive (29Nov2021)
     231                                                }
     232                                                element->AddInput(FrictionSchoofWaterPressureEnum,p_water,P1DGEnum);
     233                                                delete gauss;
     234                                                delete friction;
     235                                        }
     236                                        break;                          default:
    217237                                        _error_("Field "<<EnumToStringx(fields[j])<<" does not support stochasticity yet.");
    218238                        }

  • issm/trunk-jpl/src/c/shared/Enum/Enum.vim

    r27154 r27161  
    754754syn keyword cConstant FrictionSedimentCompressibilityCoefficientEnum
    755755syn keyword cConstant FrictionTillFrictionAngleEnum
     756syn keyword cConstant FrictionSchoofWaterPressureEnum
    756757syn keyword cConstant FrictionWaterLayerEnum
     758syn keyword cConstant FrictionWaterPressureEnum
    757759syn keyword cConstant FrictionfEnum
    758760syn keyword cConstant FrontalForcingsBasinIdEnum
    … …  
    11231125syn keyword cConstant WaterfractionEnum
    11241126syn keyword cConstant WaterheightEnum
    1125 syn keyword cConstant FrictionWaterPressureEnum
    1126 syn keyword cConstant FrictionWaterPressureNoiseEnum
    11271127syn keyword cConstant WeightsLevelsetObservationEnum
    11281128syn keyword cConstant WeightsSurfaceObservationEnum
    … …  
    16391639syn keyword cType Cfsurfacesquare
    16401640syn keyword cType Channel
     1641syn keyword cType classes
    16411642syn keyword cType Constraint
    16421643syn keyword cType Constraints
    … …  
    16451646syn keyword cType ControlInput
    16461647syn keyword cType Covertree
     1648syn keyword cType DatasetInput
    16471649syn keyword cType DataSetParam
    1648 syn keyword cType DatasetInput
    16491650syn keyword cType Definition
    16501651syn keyword cType DependentObject
    … …  
    16591660syn keyword cType ElementInput
    16601661syn keyword cType ElementMatrix
     1662syn keyword cType Elements
    16611663syn keyword cType ElementVector
    1662 syn keyword cType Elements
    16631664syn keyword cType ExponentialVariogram
    16641665syn keyword cType ExternalResult
    … …  
    16671668syn keyword cType Friction
    16681669syn keyword cType Gauss
     1670syn keyword cType GaussianVariogram
     1671syn keyword cType gaussobjects
    16691672syn keyword cType GaussPenta
    16701673syn keyword cType GaussSeg
    16711674syn keyword cType GaussTetra
    16721675syn keyword cType GaussTria
    1673 syn keyword cType GaussianVariogram
    16741676syn keyword cType GenericExternalResult
    16751677syn keyword cType GenericOption
    … …  
    16881690syn keyword cType IssmDirectApplicInterface
    16891691syn keyword cType IssmParallelDirectApplicInterface
     1692syn keyword cType krigingobjects
    16901693syn keyword cType Load
    16911694syn keyword cType Loads
    … …  
    16981701syn keyword cType Matice
    16991702syn keyword cType Matlitho
     1703syn keyword cType matrixobjects
    17001704syn keyword cType MatrixParam
    17011705syn keyword cType Misfit
    … …  
    17101714syn keyword cType Observations
    17111715syn keyword cType Option
     1716syn keyword cType Options
    17121717syn keyword cType OptionUtilities
    1713 syn keyword cType Options
    17141718syn keyword cType Param
    17151719syn keyword cType Parameters
    … …  
    17251729syn keyword cType Regionaloutput
    17261730syn keyword cType Results
     1731syn keyword cType Riftfront
    17271732syn keyword cType RiftStruct
    1728 syn keyword cType Riftfront
    17291733syn keyword cType SealevelGeometry
    17301734syn keyword cType Seg
    17311735syn keyword cType SegInput
     1736syn keyword cType Segment
    17321737syn keyword cType SegRef
    1733 syn keyword cType Segment
    17341738syn keyword cType SpcDynamic
    17351739syn keyword cType SpcStatic
    … …  
    17501754syn keyword cType Vertex
    17511755syn keyword cType Vertices
    1752 syn keyword cType classes
    1753 syn keyword cType gaussobjects
    1754 syn keyword cType krigingobjects
    1755 syn keyword cType matrixobjects
    17561756syn keyword cType AdjointBalancethickness2Analysis
    17571757syn keyword cType AdjointBalancethicknessAnalysis

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

    r27154 r27161  
    750750        FrictionSedimentCompressibilityCoefficientEnum,
    751751        FrictionTillFrictionAngleEnum,
     752        FrictionSchoofWaterPressureEnum,
    752753        FrictionWaterLayerEnum,
     754        FrictionWaterPressureEnum,
    753755        FrictionfEnum,
    754756        FrontalForcingsBasinIdEnum,
    … …  
    11201122        WaterfractionEnum,
    11211123        WaterheightEnum,
    1122         FrictionWaterPressureEnum,
    1123         FrictionWaterPressureNoiseEnum,
    11241124        WeightsLevelsetObservationEnum,
    11251125        WeightsSurfaceObservationEnum,

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

    r27154 r27161  
    756756                case FrictionSedimentCompressibilityCoefficientEnum : return "FrictionSedimentCompressibilityCoefficient";
    757757                case FrictionTillFrictionAngleEnum : return "FrictionTillFrictionAngle";
     758                case FrictionSchoofWaterPressureEnum : return "FrictionSchoofWaterPressure";
    758759                case FrictionWaterLayerEnum : return "FrictionWaterLayer";
     760                case FrictionWaterPressureEnum : return "FrictionWaterPressure";
    759761                case FrictionfEnum : return "Frictionf";
    760762                case FrontalForcingsBasinIdEnum : return "FrontalForcingsBasinId";
    … …  
    11251127                case WaterfractionEnum : return "Waterfraction";
    11261128                case WaterheightEnum : return "Waterheight";
    1127                 case FrictionWaterPressureEnum : return "FrictionWaterPressure";
    1128                 case FrictionWaterPressureNoiseEnum : return "FrictionWaterPressureNoise";
    11291129                case WeightsLevelsetObservationEnum : return "WeightsLevelsetObservation";
    11301130                case WeightsSurfaceObservationEnum : return "WeightsSurfaceObservation";

  • issm/trunk-jpl/src/c/shared/Enum/Enumjl.vim

    r27154 r27161  
    747747syn keyword juliaConstC FrictionSedimentCompressibilityCoefficientEnum
    748748syn keyword juliaConstC FrictionTillFrictionAngleEnum
     749syn keyword juliaConstC FrictionSchoofWaterPressureEnum
    749750syn keyword juliaConstC FrictionWaterLayerEnum
     751syn keyword juliaConstC FrictionWaterPressureEnum
    750752syn keyword juliaConstC FrictionfEnum
    751753syn keyword juliaConstC FrontalForcingsBasinIdEnum
    … …  
    11161118syn keyword juliaConstC WaterfractionEnum
    11171119syn keyword juliaConstC WaterheightEnum
    1118 syn keyword juliaConstC FrictionWaterPressureEnum
    1119 syn keyword juliaConstC FrictionWaterPressureNoiseEnum
    11201120syn keyword juliaConstC WeightsLevelsetObservationEnum
    11211121syn keyword juliaConstC WeightsSurfaceObservationEnum

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

    r27154 r27161  
    774774              else if (strcmp(name,"FrictionSedimentCompressibilityCoefficient")==0) return FrictionSedimentCompressibilityCoefficientEnum;
    775775              else if (strcmp(name,"FrictionTillFrictionAngle")==0) return FrictionTillFrictionAngleEnum;
     776              else if (strcmp(name,"FrictionSchoofWaterPressure")==0) return FrictionSchoofWaterPressureEnum;
    776777              else if (strcmp(name,"FrictionWaterLayer")==0) return FrictionWaterLayerEnum;
     778              else if (strcmp(name,"FrictionWaterPressure")==0) return FrictionWaterPressureEnum;
    777779              else if (strcmp(name,"Frictionf")==0) return FrictionfEnum;
    778780              else if (strcmp(name,"FrontalForcingsBasinId")==0) return FrontalForcingsBasinIdEnum;
    … …  
    873875              else if (strcmp(name,"RadarIcePeriod")==0) return RadarIcePeriodEnum;
    874876              else if (strcmp(name,"RadarPowerMacGregor")==0) return RadarPowerMacGregorEnum;
    875               else if (strcmp(name,"RadarPowerWolff")==0) return RadarPowerWolffEnum;
    876               else if (strcmp(name,"RheologyBAbsGradient")==0) return RheologyBAbsGradientEnum;
    877877         else stage=8;
    878878   }
    879879   if(stage==8){
    880               if (strcmp(name,"RheologyBInitialguess")==0) return RheologyBInitialguessEnum;
     880              if (strcmp(name,"RadarPowerWolff")==0) return RadarPowerWolffEnum;
     881              else if (strcmp(name,"RheologyBAbsGradient")==0) return RheologyBAbsGradientEnum;
     882              else if (strcmp(name,"RheologyBInitialguess")==0) return RheologyBInitialguessEnum;
    881883              else if (strcmp(name,"RheologyBInitialguessMisfit")==0) return RheologyBInitialguessMisfitEnum;
    882884              else if (strcmp(name,"RheologyBbarAbsGradient")==0) return RheologyBbarAbsGradientEnum;
    … …  
    996998              else if (strcmp(name,"SmbDzAdd")==0) return SmbDzAddEnum;
    997999              else if (strcmp(name,"SmbDz")==0) return SmbDzEnum;
    998               else if (strcmp(name,"SmbDzMin")==0) return SmbDzMinEnum;
    999               else if (strcmp(name,"SmbDzTop")==0) return SmbDzTopEnum;
    10001000         else stage=9;
    10011001   }
    10021002   if(stage==9){
    1003               if (strcmp(name,"SmbDzini")==0) return SmbDziniEnum;
     1003              if (strcmp(name,"SmbDzMin")==0) return SmbDzMinEnum;
     1004              else if (strcmp(name,"SmbDzTop")==0) return SmbDzTopEnum;
     1005              else if (strcmp(name,"SmbDzini")==0) return SmbDziniEnum;
    10041006              else if (strcmp(name,"SmbEAir")==0) return SmbEAirEnum;
    10051007              else if (strcmp(name,"SmbEC")==0) return SmbECEnum;
    … …  
    11191121              else if (strcmp(name,"ThicknessAlongGradient")==0) return ThicknessAlongGradientEnum;
    11201122              else if (strcmp(name,"Thickness")==0) return ThicknessEnum;
    1121               else if (strcmp(name,"ThicknessOld")==0) return ThicknessOldEnum;
    1122               else if (strcmp(name,"ThicknessPositive")==0) return ThicknessPositiveEnum;
    11231123         else stage=10;
    11241124   }
    11251125   if(stage==10){
    1126               if (strcmp(name,"ThicknessResidual")==0) return ThicknessResidualEnum;
     1126              if (strcmp(name,"ThicknessOld")==0) return ThicknessOldEnum;
     1127              else if (strcmp(name,"ThicknessPositive")==0) return ThicknessPositiveEnum;
     1128              else if (strcmp(name,"ThicknessResidual")==0) return ThicknessResidualEnum;
    11271129              else if (strcmp(name,"TransientAccumulatedDeltaIceThickness")==0) return TransientAccumulatedDeltaIceThicknessEnum;
    11281130              else if (strcmp(name,"Vel")==0) return VelEnum;
    … …  
    11521154              else if (strcmp(name,"Waterfraction")==0) return WaterfractionEnum;
    11531155              else if (strcmp(name,"Waterheight")==0) return WaterheightEnum;
    1154               else if (strcmp(name,"FrictionWaterPressure")==0) return FrictionWaterPressureEnum;
    1155               else if (strcmp(name,"FrictionWaterPressureNoise")==0) return FrictionWaterPressureNoiseEnum;
    11561156              else if (strcmp(name,"WeightsLevelsetObservation")==0) return WeightsLevelsetObservationEnum;
    11571157              else if (strcmp(name,"WeightsSurfaceObservation")==0) return WeightsSurfaceObservationEnum;

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

    r26859 r27161  
    9595                                end
    9696                                if(contains(field,'WaterPressure'))
    97                mdname = structstoch.mdnames(find(strcmp(structstoch.fields,char(field))));
    98                if~(isequal(class(md.friction),char(mdname)))
    99                   error('stochasticforcing field %s is only implemented for default friction', char(field));
     97                                        mdname = structstoch.mdnames(find(strcmp(structstoch.fields,char(field))));
     98                                        if~(isequal(class(md.friction),char(mdname)))
     99                  error('md.friction does not agree with stochasticforcing field %s', char(field));
    100100               end
    101                if(md.friction.coupling~=0 && md.friction.coupling~=1 && md.friction.coupling~=2)
    102                   error('stochasticforcing field %s is only implemented for cases md.friction.coupling 0 or 1 or 2', char(field));
    103                end
    104                if(any(md.friction.q==0))
    105                   error('stochasticforcing field %s requires non-zero q exponent',char(field));
    106                end
     101                                        if(strcmp(class(md.friction),'friction') || strcmp(class(md.friction),'frictionschoof'))
     102                  if(md.friction.coupling~=0 && md.friction.coupling~=1 && md.friction.coupling~=2)
     103                     error('stochasticforcing field %s is only implemented for cases md.friction.coupling 0 or 1 or 2', char(field));
     104                  end
     105                                        end
     106                                        if(strcmp(class(md.friction),'friction'))
     107                  if(any(md.friction.q==0))
     108                     error('stochasticforcing field %s requires non-zero q exponent',char(field));
     109                  end
     110                                        end
    107111            end
    108112                                %Checking for specific dimensions
    … …  
    267271                'FloatingMeltRate',...
    268272                'FrictionWaterPressure',...
     273                'FrictionSchoofWaterPressure',...
    269274                'FrontalForcingsRignotAutoregression',...
    270275                'SMBautoregression',...
    … …  
    277282                'basalforcings',...
    278283                'friction',...
     284                'frictionschoof',...
    279285                'frontalforcingsrignotautoregression',...
    280286                'SMBautoregression',...

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

    r26905 r27161  
    106106                mdname = structstoch[field]
    107107                if (type(md.friction).__name__ != mdname):
    108                     raise TypeError('stochasticforcing field {} is only implemented for default friction'.format(field))
    109                 if md.friction.coupling not in[0, 1, 2]:
    110                     raise TypeError('stochasticforcing field {} is only implemented for cases md.friction.coupling 0 or 1 or 2'.format(field))
    111                 if (np.any(md.friction.q == 0)):
    112                     raise TypeError('stochasticforcing field {} requires non-zero q exponent'.format(field))
     108                    raise TypeError('md.friction does not agree with stochasticforcing field {}'.format(field))
     109                if (type(md.friction).__name__=='friction' or type(md.friction).__name__=='frictionschoof'):
     110                    if md.friction.coupling not in[0, 1, 2]:
     111                        raise TypeError('stochasticforcing field {} is only implemented for cases md.friction.coupling 0 or 1 or 2'.format(field))
     112                if (type(md.friction).__name__=='friction'):
     113                    if (np.any(md.friction.q == 0)):
     114                        raise TypeError('stochasticforcing field {} requires non-zero q exponent'.format(field))
    113115
    114116            # Checking for specific dimensions
    … …  
    232234                     'FloatingMeltRate': 'basalforcings',
    233235                     'FrictionWaterPressure': 'friction',
     236                     'FrictionSchoofWaterPressure': 'frictionschoof',
    234237                     'FrontalForcingsRignotAutoregression': 'frontalforcingsrignotautoregression',
    235238                     'SMBautoregression': 'SMBautoregression',
Note: See TracChangeset for help on using the changeset viewer.