Changeset 26676


Ignore:
Timestamp:
11/30/21 05:29:36 (3 years ago)
Author:
vverjans
Message:

NEW: added FrictionWaterPressure to stochasticforcingfields, new functions Friction::IcePressure() and Friction::SubglacialWaterPressure(), clean-up Friction::EffectivePressure()

Location:
issm/trunk-jpl
Files:
3 added
11 edited

Legend:

Unmodified
Added
Removed

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

    r26601 r26676  
    545545
    546546        /*Get effective pressure and basal velocity*/
    547         IssmDouble Neff = EffectivePressure(gauss);
    548547        IssmDouble vmag = VelMag(gauss);
     548
     549        bool ispwStochastic;
     550        IssmDouble Neff;
     551        element->parameters->FindParam(&ispwStochastic,StochasticForcingIsWaterPressureEnum);
     552        if(ispwStochastic){
     553                /*Retrieve stochastic water pressure and compute ice pressure*/
     554                IssmDouble p_ice,p_water,Neff_limit;
     555                element->GetInputValue(&p_water,gauss,FrictionWaterPressureEnum);
     556                element->parameters->FindParam(&Neff_limit,FrictionEffectivePressureLimitEnum);
     557                p_ice = IcePressure(gauss);
     558                Neff  = max(Neff_limit*p_ice, p_ice - p_water);
     559        }       
     560        else{
     561                /*Compute effective pressure directly*/
     562                Neff = EffectivePressure(gauss);
     563        }
    549564
    550565        /*Check to prevent dividing by zero if vmag==0*/
    … …  
    845860        element->parameters->FindParam(&Neff_limit,FrictionEffectivePressureLimitEnum);
    846861
     862        /*Compute ice pressure*/
     863        p_ice = IcePressure(gauss);
     864
    847865        /*From base and thickness, compute effective pressure when drag is viscous, or get Neff from forcing:*/
    848866        switch(coupled_flag){
    849867                case 0:{
    850                         element->GetInputValue(&thickness, gauss,ThicknessEnum);
    851868                        element->GetInputValue(&base, gauss,BaseEnum);
    852869                        element->GetInputValue(&sealevel, gauss,SealevelEnum);
    853870                        IssmDouble rho_water = element->FindParam(MaterialsRhoSeawaterEnum);
    854                         IssmDouble rho_ice   = element->FindParam(MaterialsRhoIceEnum);
    855871                        IssmDouble gravity   = element->FindParam(ConstantsGEnum);
    856                         p_ice   = gravity*rho_ice*thickness;
    857872                        p_water = rho_water*gravity*(sealevel-base);
    858873                        Neff = p_ice - p_water;
    … …  
    860875                        break;
    861876                case 1:{
    862                         element->GetInputValue(&thickness, gauss,ThicknessEnum);
    863                         IssmDouble rho_ice   = element->FindParam(MaterialsRhoIceEnum);
    864                         IssmDouble gravity   = element->FindParam(ConstantsGEnum);
    865                         p_ice   = gravity*rho_ice*thickness;
    866877                        p_water = 0.;
    867878                        Neff = p_ice - p_water;
    … …  
    869880                        break;
    870881                case 2:{
    871                         element->GetInputValue(&thickness, gauss,ThicknessEnum);
    872882                        element->GetInputValue(&base, gauss,BaseEnum);
    873883                        element->GetInputValue(&sealevel, gauss,SealevelEnum);
    874884                        IssmDouble rho_water = element->FindParam(MaterialsRhoSeawaterEnum);
    875                         IssmDouble rho_ice   = element->FindParam(MaterialsRhoIceEnum);
    876885                        IssmDouble gravity   = element->FindParam(ConstantsGEnum);
    877                         p_ice   = gravity*rho_ice*thickness;
    878886                        p_water = max(0.,rho_water*gravity*(sealevel-base));
    879887                        Neff = p_ice - p_water;
    … …  
    882890                case 3:{
    883891                        element->GetInputValue(&Neff,gauss,FrictionEffectivePressureEnum);
    884                         element->GetInputValue(&thickness, gauss,ThicknessEnum);
    885                         IssmDouble rho_ice   = element->FindParam(MaterialsRhoIceEnum);
    886                         IssmDouble gravity   = element->FindParam(ConstantsGEnum);
    887                         p_ice   = gravity*rho_ice*thickness;
    888892                }
    889893                        break;
    890894                case 4:{
    891895                        element->GetInputValue(&Neff,gauss,EffectivePressureEnum);
    892                         element->GetInputValue(&thickness, gauss,ThicknessEnum);
    893                         IssmDouble rho_ice   = element->FindParam(MaterialsRhoIceEnum);
    894                         IssmDouble gravity   = element->FindParam(ConstantsGEnum);
    895                         p_ice   = gravity*rho_ice*thickness;
    896896                }
    897897                        break;
    … …  
    905905        /*Return effective pressure*/
    906906        return Neff;
     907
     908}/*}}}*/
     909IssmDouble Friction::IcePressure(Gauss* gauss){/*{{{*/
     910        /*Get ice pressure*/
     911
     912        IssmDouble  thickness,p_ice;
     913        /*Recover Inputs and Parameters*/
     914        element->GetInputValue(&thickness, gauss,ThicknessEnum);
     915        IssmDouble rho_ice = element->FindParam(MaterialsRhoIceEnum);
     916        IssmDouble gravity = element->FindParam(ConstantsGEnum);
     917
     918        /*Compute*/
     919        p_ice = gravity*rho_ice*thickness;
     920
     921        /*Return ice pressure*/
     922        return p_ice;
     923
     924}/*}}}*/
     925IssmDouble Friction::SubglacialWaterPressure(Gauss* gauss){/*{{{*/
     926        /*Get water pressure as a function of  flag */
     927
     928        int         coupled_flag;
     929        IssmDouble  base,sealevel,p_water;
     930
     931        /*Recover parameters: */
     932        element->parameters->FindParam(&coupled_flag,FrictionCouplingEnum);
     933
     934        switch(coupled_flag){
     935                case 0:{
     936                        element->GetInputValue(&base, gauss,BaseEnum);
     937                        element->GetInputValue(&sealevel, gauss,SealevelEnum);
     938                        IssmDouble rho_water = element->FindParam(MaterialsRhoSeawaterEnum);
     939                        IssmDouble gravity   = element->FindParam(ConstantsGEnum);
     940                        p_water = rho_water*gravity*(sealevel-base);
     941                }
     942                        break;
     943                case 1:{
     944                        p_water = 0.;
     945                }
     946                        break;
     947                case 2:{
     948                        element->GetInputValue(&base, gauss,BaseEnum);
     949                        element->GetInputValue(&sealevel, gauss,SealevelEnum);
     950                        IssmDouble rho_water = element->FindParam(MaterialsRhoSeawaterEnum);
     951                        IssmDouble gravity   = element->FindParam(ConstantsGEnum);
     952                        p_water = max(0.,rho_water*gravity*(sealevel-base));
     953                }
     954                        break;
     955                case 3:{
     956                        _error_("water pressure not computed for coupling==3 in friction law");
     957                }
     958                        break;
     959                case 4:{
     960                        _error_("water pressure not computed for coupling==4 in friction law");
     961                }
     962                        break;
     963                default:
     964                        _error_("not supported");
     965        }
     966
     967        /*Return water pressure*/
     968        return p_water;
    907969
    908970}/*}}}*/

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

    r26470 r26676  
    5252
    5353                IssmDouble EffectivePressure(Gauss* gauss);
     54                IssmDouble IcePressure(Gauss* gauss);
     55                IssmDouble SubglacialWaterPressure(Gauss* gauss);
    5456                IssmDouble VelMag(Gauss* gauss);
    5557                void GetBasalSlidingSpeeds(IssmDouble* pvx, Gauss* gauss);

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

    r26615 r26676  
    501501        bool isstochasticforcing;
    502502   parameters->FindParam(&isstochasticforcing,StochasticForcingIsStochasticForcingEnum);
     503        /*Stochastic Effective Pressure false by default*/
     504        parameters->AddObject(new BoolParam(StochasticForcingIsWaterPressureEnum,false));
    503505   if(isstochasticforcing){
    504506      int num_fields,stochastic_dim;

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

    r26640 r26676  
    44
    55#include "./StochasticForcingx.h"
     6#include "../../classes/Loads/Friction.h"
    67#include "../../shared/shared.h"
    78#include "../../toolkits/toolkits.h"
    … …  
    141142                                        }
    142143                                        break;
     144                                //VV working(29Nov2021)
     145                                case FrictionWaterPressureEnum:
     146                                        /*Specify that WaterPressure is stochastic*/
     147                                        femmodel->parameters->SetParam(true,StochasticForcingIsWaterPressureEnum);
     148                                        for(Object* &object:femmodel->elements->objects){
     149                  Element* element = xDynamicCast<Element*>(object);
     150                  int numvertices  = element->GetNumberOfVertices();
     151                  IssmDouble p_water_deterministic[numvertices];
     152                  IssmDouble p_water[numvertices];
     153                                                element->GetInputValue(&dimensionid,StochasticForcingDefaultIdEnum);
     154                                                Gauss* gauss=element->NewGauss();
     155                                                Friction* friction = new Friction(element);
     156                                                for(int i=0;i<numvertices;i++){
     157                                                        gauss->GaussVertex(i);
     158                                                        p_water_deterministic[i] = friction->SubglacialWaterPressure(gauss);
     159                                                        p_water[i]               = p_water_deterministic[i] + noisefield[dimensionid]; //make sure positive (29Nov2021)
     160                                                        p_water[i]               = max(0.0,p_water[i]);
     161                                                }
     162                                                element->AddInput(FrictionWaterPressureEnum,p_water,P1DGEnum);
     163                                                delete gauss;
     164                                                delete friction;
     165                                        }
     166                                        break;
     167
     168                                /*
     169                                case FrictionEffectivePressureEnum:
     170                                        femmodel->parameters->SetParam(true,StochasticForcingIsEffectivePressureEnum);
     171                                        for(Object* &object:femmodel->elements->objects){
     172                  Element* element = xDynamicCast<Element*>(object);
     173                  int numvertices  = element->GetNumberOfVertices();
     174                  IssmDouble Neff[numvertices];
     175                                                element->GetInputValue(&dimensionid,StochasticForcingDefaultIdEnum);
     176                                                Gauss* gauss=element->NewGauss();
     177                                                Friction* friction = new Friction(element);
     178                                                for(int i=0;i<numvertices;i++){
     179                                                        gauss->GaussVertex(i);
     180                                                        Neff[i] = friction->EffectivePressure(gauss);
     181                                                        Neff[i] = Neff[i]+noisefield[dimensionid];
     182                                                }
     183                                                element->AddInput(FrictionEffectivePressureEnum,Neff,P1DGEnum);
     184                                                delete gauss;
     185                                                delete friction;
     186                                        }
     187                                        break;
     188                                */
     189
    143190                                default:
    144191                                        _error_("Field "<<EnumToStringx(fields[j])<<" does not support stochasticity yet.");

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

    r26640 r26676  
    401401syn keyword cConstant StochasticForcingDimensionsEnum
    402402syn keyword cConstant StochasticForcingFieldsEnum
     403syn keyword cConstant StochasticForcingIsEffectivePressureEnum
    403404syn keyword cConstant StochasticForcingIsStochasticForcingEnum
     405syn keyword cConstant StochasticForcingIsWaterPressureEnum
    404406syn keyword cConstant StochasticForcingNumFieldsEnum
    405407syn keyword cConstant StochasticForcingRandomflagEnum
    … …  
    603605syn keyword cConstant BaselineBasalforcingsFloatingiceMeltingRateEnum
    604606syn keyword cConstant BaselineCalvingCalvingrateEnum
     607syn keyword cConstant BaselineFrictionEffectivePressureEnum
    605608syn keyword cConstant BedEnum
    606609syn keyword cConstant BedGRDEnum
    … …  
    10731076syn keyword cConstant WaterfractionEnum
    10741077syn keyword cConstant WaterheightEnum
     1078syn keyword cConstant FrictionWaterPressureEnum
     1079syn keyword cConstant FrictionWaterPressureNoiseEnum
    10751080syn keyword cConstant WeightsLevelsetObservationEnum
    10761081syn keyword cConstant WeightsSurfaceObservationEnum

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

    r26640 r26676  
    395395        StochasticForcingDimensionsEnum,
    396396        StochasticForcingFieldsEnum,
     397        StochasticForcingIsEffectivePressureEnum,
    397398        StochasticForcingIsStochasticForcingEnum,
     399        StochasticForcingIsWaterPressureEnum,
    398400        StochasticForcingNumFieldsEnum,
    399401        StochasticForcingRandomflagEnum,
    … …  
    599601        BaselineBasalforcingsFloatingiceMeltingRateEnum,
    600602        BaselineCalvingCalvingrateEnum,
     603        BaselineFrictionEffectivePressureEnum,
    601604        BedEnum,
    602605        BedGRDEnum,
    … …  
    10701073        WaterfractionEnum,
    10711074        WaterheightEnum,
     1075        FrictionWaterPressureEnum,
     1076        FrictionWaterPressureNoiseEnum,
    10721077        WeightsLevelsetObservationEnum,
    10731078        WeightsSurfaceObservationEnum,

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

    r26640 r26676  
    403403                case StochasticForcingDimensionsEnum : return "StochasticForcingDimensions";
    404404                case StochasticForcingFieldsEnum : return "StochasticForcingFields";
     405                case StochasticForcingIsEffectivePressureEnum : return "StochasticForcingIsEffectivePressure";
    405406                case StochasticForcingIsStochasticForcingEnum : return "StochasticForcingIsStochasticForcing";
     407                case StochasticForcingIsWaterPressureEnum : return "StochasticForcingIsWaterPressure";
    406408                case StochasticForcingNumFieldsEnum : return "StochasticForcingNumFields";
    407409                case StochasticForcingRandomflagEnum : return "StochasticForcingRandomflag";
    … …  
    605607                case BaselineBasalforcingsFloatingiceMeltingRateEnum : return "BaselineBasalforcingsFloatingiceMeltingRate";
    606608                case BaselineCalvingCalvingrateEnum : return "BaselineCalvingCalvingrate";
     609                case BaselineFrictionEffectivePressureEnum : return "BaselineFrictionEffectivePressure";
    607610                case BedEnum : return "Bed";
    608611                case BedGRDEnum : return "BedGRD";
    … …  
    10751078                case WaterfractionEnum : return "Waterfraction";
    10761079                case WaterheightEnum : return "Waterheight";
     1080                case FrictionWaterPressureEnum : return "FrictionWaterPressure";
     1081                case FrictionWaterPressureNoiseEnum : return "FrictionWaterPressureNoise";
    10771082                case WeightsLevelsetObservationEnum : return "WeightsLevelsetObservation";
    10781083                case WeightsSurfaceObservationEnum : return "WeightsSurfaceObservation";

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

    r26640 r26676  
    394394syn keyword juliaConstC StochasticForcingDimensionsEnum
    395395syn keyword juliaConstC StochasticForcingFieldsEnum
     396syn keyword juliaConstC StochasticForcingIsEffectivePressureEnum
    396397syn keyword juliaConstC StochasticForcingIsStochasticForcingEnum
     398syn keyword juliaConstC StochasticForcingIsWaterPressureEnum
    397399syn keyword juliaConstC StochasticForcingNumFieldsEnum
    398400syn keyword juliaConstC StochasticForcingRandomflagEnum
    … …  
    596598syn keyword juliaConstC BaselineBasalforcingsFloatingiceMeltingRateEnum
    597599syn keyword juliaConstC BaselineCalvingCalvingrateEnum
     600syn keyword juliaConstC BaselineFrictionEffectivePressureEnum
    598601syn keyword juliaConstC BedEnum
    599602syn keyword juliaConstC BedGRDEnum
    … …  
    10661069syn keyword juliaConstC WaterfractionEnum
    10671070syn keyword juliaConstC WaterheightEnum
     1071syn keyword juliaConstC FrictionWaterPressureEnum
     1072syn keyword juliaConstC FrictionWaterPressureNoiseEnum
    10681073syn keyword juliaConstC WeightsLevelsetObservationEnum
    10691074syn keyword juliaConstC WeightsSurfaceObservationEnum

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

    r26640 r26676  
    412412              else if (strcmp(name,"StochasticForcingDimensions")==0) return StochasticForcingDimensionsEnum;
    413413              else if (strcmp(name,"StochasticForcingFields")==0) return StochasticForcingFieldsEnum;
     414              else if (strcmp(name,"StochasticForcingIsEffectivePressure")==0) return StochasticForcingIsEffectivePressureEnum;
    414415              else if (strcmp(name,"StochasticForcingIsStochasticForcing")==0) return StochasticForcingIsStochasticForcingEnum;
     416              else if (strcmp(name,"StochasticForcingIsWaterPressure")==0) return StochasticForcingIsWaterPressureEnum;
    415417              else if (strcmp(name,"StochasticForcingNumFields")==0) return StochasticForcingNumFieldsEnum;
    416418              else if (strcmp(name,"StochasticForcingRandomflag")==0) return StochasticForcingRandomflagEnum;
    … …  
    504506              else if (strcmp(name,"Step")==0) return StepEnum;
    505507              else if (strcmp(name,"Steps")==0) return StepsEnum;
    506               else if (strcmp(name,"StressbalanceAbstol")==0) return StressbalanceAbstolEnum;
    507               else if (strcmp(name,"StressbalanceFSreconditioning")==0) return StressbalanceFSreconditioningEnum;
    508508         else stage=5;
    509509   }
    510510   if(stage==5){
    511               if (strcmp(name,"StressbalanceIsnewton")==0) return StressbalanceIsnewtonEnum;
     511              if (strcmp(name,"StressbalanceAbstol")==0) return StressbalanceAbstolEnum;
     512              else if (strcmp(name,"StressbalanceFSreconditioning")==0) return StressbalanceFSreconditioningEnum;
     513              else if (strcmp(name,"StressbalanceIsnewton")==0) return StressbalanceIsnewtonEnum;
    512514              else if (strcmp(name,"StressbalanceMaxiter")==0) return StressbalanceMaxiterEnum;
    513515              else if (strcmp(name,"StressbalanceNumRequestedOutputs")==0) return StressbalanceNumRequestedOutputsEnum;
    … …  
    617619              else if (strcmp(name,"BaselineBasalforcingsFloatingiceMeltingRate")==0) return BaselineBasalforcingsFloatingiceMeltingRateEnum;
    618620              else if (strcmp(name,"BaselineCalvingCalvingrate")==0) return BaselineCalvingCalvingrateEnum;
     621              else if (strcmp(name,"BaselineFrictionEffectivePressure")==0) return BaselineFrictionEffectivePressureEnum;
    619622              else if (strcmp(name,"Bed")==0) return BedEnum;
    620623              else if (strcmp(name,"BedGRD")==0) return BedGRDEnum;
    … …  
    626629              else if (strcmp(name,"BedSlopeY")==0) return BedSlopeYEnum;
    627630              else if (strcmp(name,"BottomPressure")==0) return BottomPressureEnum;
    628               else if (strcmp(name,"BottomPressureOld")==0) return BottomPressureOldEnum;
    629               else if (strcmp(name,"CalvingCalvingrate")==0) return CalvingCalvingrateEnum;
    630               else if (strcmp(name,"CalvingHabFraction")==0) return CalvingHabFractionEnum;
    631631         else stage=6;
    632632   }
    633633   if(stage==6){
    634               if (strcmp(name,"CalvingMeltingrate")==0) return CalvingMeltingrateEnum;
     634              if (strcmp(name,"BottomPressureOld")==0) return BottomPressureOldEnum;
     635              else if (strcmp(name,"CalvingCalvingrate")==0) return CalvingCalvingrateEnum;
     636              else if (strcmp(name,"CalvingHabFraction")==0) return CalvingHabFractionEnum;
     637              else if (strcmp(name,"CalvingMeltingrate")==0) return CalvingMeltingrateEnum;
    635638              else if (strcmp(name,"CalvingStressThresholdFloatingice")==0) return CalvingStressThresholdFloatingiceEnum;
    636639              else if (strcmp(name,"CalvingStressThresholdGroundedice")==0) return CalvingStressThresholdGroundediceEnum;
    … …  
    749752              else if (strcmp(name,"HydrologyHeadOld")==0) return HydrologyHeadOldEnum;
    750753              else if (strcmp(name,"HydrologyMoulinInput")==0) return HydrologyMoulinInputEnum;
    751               else if (strcmp(name,"HydrologyNeumannflux")==0) return HydrologyNeumannfluxEnum;
    752               else if (strcmp(name,"HydrologyReynolds")==0) return HydrologyReynoldsEnum;
    753               else if (strcmp(name,"HydrologySheetConductivity")==0) return HydrologySheetConductivityEnum;
    754754         else stage=7;
    755755   }
    756756   if(stage==7){
    757               if (strcmp(name,"HydrologySheetThickness")==0) return HydrologySheetThicknessEnum;
     757              if (strcmp(name,"HydrologyNeumannflux")==0) return HydrologyNeumannfluxEnum;
     758              else if (strcmp(name,"HydrologyReynolds")==0) return HydrologyReynoldsEnum;
     759              else if (strcmp(name,"HydrologySheetConductivity")==0) return HydrologySheetConductivityEnum;
     760              else if (strcmp(name,"HydrologySheetThickness")==0) return HydrologySheetThicknessEnum;
    758761              else if (strcmp(name,"HydrologySheetThicknessOld")==0) return HydrologySheetThicknessOldEnum;
    759762              else if (strcmp(name,"HydrologyTws")==0) return HydrologyTwsEnum;
    … …  
    872875              else if (strcmp(name,"SealevelUGrd")==0) return SealevelUGrdEnum;
    873876              else if (strcmp(name,"SealevelNGrd")==0) return SealevelNGrdEnum;
    874               else if (strcmp(name,"SealevelUEastEsa")==0) return SealevelUEastEsaEnum;
    875               else if (strcmp(name,"SealevelUNorthEsa")==0) return SealevelUNorthEsaEnum;
    876               else if (strcmp(name,"SealevelchangeIndices")==0) return SealevelchangeIndicesEnum;
    877877         else stage=8;
    878878   }
    879879   if(stage==8){
    880               if (strcmp(name,"SealevelchangeG")==0) return SealevelchangeGEnum;
     880              if (strcmp(name,"SealevelUEastEsa")==0) return SealevelUEastEsaEnum;
     881              else if (strcmp(name,"SealevelUNorthEsa")==0) return SealevelUNorthEsaEnum;
     882              else if (strcmp(name,"SealevelchangeIndices")==0) return SealevelchangeIndicesEnum;
     883              else if (strcmp(name,"SealevelchangeG")==0) return SealevelchangeGEnum;
    881884              else if (strcmp(name,"SealevelchangeGU")==0) return SealevelchangeGUEnum;
    882885              else if (strcmp(name,"SealevelchangeGE")==0) return SealevelchangeGEEnum;
    … …  
    995998              else if (strcmp(name,"SmbSizeini")==0) return SmbSizeiniEnum;
    996999              else if (strcmp(name,"SmbSmbCorr")==0) return SmbSmbCorrEnum;
    997               else if (strcmp(name,"SmbSmbref")==0) return SmbSmbrefEnum;
    998               else if (strcmp(name,"SmbSzaValue")==0) return SmbSzaValueEnum;
    999               else if (strcmp(name,"SmbT")==0) return SmbTEnum;
    10001000         else stage=9;
    10011001   }
    10021002   if(stage==9){
    1003               if (strcmp(name,"SmbTa")==0) return SmbTaEnum;
     1003              if (strcmp(name,"SmbSmbref")==0) return SmbSmbrefEnum;
     1004              else if (strcmp(name,"SmbSzaValue")==0) return SmbSzaValueEnum;
     1005              else if (strcmp(name,"SmbT")==0) return SmbTEnum;
     1006              else if (strcmp(name,"SmbTa")==0) return SmbTaEnum;
    10041007              else if (strcmp(name,"SmbTeValue")==0) return SmbTeValueEnum;
    10051008              else if (strcmp(name,"SmbTemperaturesAnomaly")==0) return SmbTemperaturesAnomalyEnum;
    … …  
    10991102              else if (strcmp(name,"Waterfraction")==0) return WaterfractionEnum;
    11001103              else if (strcmp(name,"Waterheight")==0) return WaterheightEnum;
     1104              else if (strcmp(name,"FrictionWaterPressure")==0) return FrictionWaterPressureEnum;
     1105              else if (strcmp(name,"FrictionWaterPressureNoise")==0) return FrictionWaterPressureNoiseEnum;
    11011106              else if (strcmp(name,"WeightsLevelsetObservation")==0) return WeightsLevelsetObservationEnum;
    11021107              else if (strcmp(name,"WeightsSurfaceObservation")==0) return WeightsSurfaceObservationEnum;
    … …  
    11161121              else if (strcmp(name,"Outputdefinition19")==0) return Outputdefinition19Enum;
    11171122              else if (strcmp(name,"Outputdefinition20")==0) return Outputdefinition20Enum;
    1118               else if (strcmp(name,"Outputdefinition21")==0) return Outputdefinition21Enum;
     1123         else stage=10;
     1124   }
     1125   if(stage==10){
     1126              if (strcmp(name,"Outputdefinition21")==0) return Outputdefinition21Enum;
    11191127              else if (strcmp(name,"Outputdefinition22")==0) return Outputdefinition22Enum;
    11201128              else if (strcmp(name,"Outputdefinition23")==0) return Outputdefinition23Enum;
    11211129              else if (strcmp(name,"Outputdefinition24")==0) return Outputdefinition24Enum;
    11221130              else if (strcmp(name,"Outputdefinition25")==0) return Outputdefinition25Enum;
    1123          else stage=10;
    1124    }
    1125    if(stage==10){
    1126               if (strcmp(name,"Outputdefinition26")==0) return Outputdefinition26Enum;
     1131              else if (strcmp(name,"Outputdefinition26")==0) return Outputdefinition26Enum;
    11271132              else if (strcmp(name,"Outputdefinition27")==0) return Outputdefinition27Enum;
    11281133              else if (strcmp(name,"Outputdefinition28")==0) return Outputdefinition28Enum;
    … …  
    12391244              else if (strcmp(name,"Boundary")==0) return BoundaryEnum;
    12401245              else if (strcmp(name,"BuddJacka")==0) return BuddJackaEnum;
    1241               else if (strcmp(name,"CalvingDev2")==0) return CalvingDev2Enum;
     1246         else stage=11;
     1247   }
     1248   if(stage==11){
     1249              if (strcmp(name,"CalvingDev2")==0) return CalvingDev2Enum;
    12421250              else if (strcmp(name,"CalvingHab")==0) return CalvingHabEnum;
    12431251              else if (strcmp(name,"CalvingLevermann")==0) return CalvingLevermannEnum;
    12441252              else if (strcmp(name,"CalvingVonmises")==0) return CalvingVonmisesEnum;
    12451253              else if (strcmp(name,"Cfdragcoeffabsgrad")==0) return CfdragcoeffabsgradEnum;
    1246          else stage=11;
    1247    }
    1248    if(stage==11){
    1249               if (strcmp(name,"Cfsurfacelogvel")==0) return CfsurfacelogvelEnum;
     1254              else if (strcmp(name,"Cfsurfacelogvel")==0) return CfsurfacelogvelEnum;
    12501255              else if (strcmp(name,"Cfsurfacesquare")==0) return CfsurfacesquareEnum;
    12511256              else if (strcmp(name,"Cflevelsetmisfit")==0) return CflevelsetmisfitEnum;
    … …  
    13621367              else if (strcmp(name,"Indexed")==0) return IndexedEnum;
    13631368              else if (strcmp(name,"IntExternalResult")==0) return IntExternalResultEnum;
    1364               else if (strcmp(name,"ElementInput")==0) return ElementInputEnum;
     1369         else stage=12;
     1370   }
     1371   if(stage==12){
     1372              if (strcmp(name,"ElementInput")==0) return ElementInputEnum;
    13651373              else if (strcmp(name,"IntMatExternalResult")==0) return IntMatExternalResultEnum;
    13661374              else if (strcmp(name,"IntMatParam")==0) return IntMatParamEnum;
    13671375              else if (strcmp(name,"IntParam")==0) return IntParamEnum;
    13681376              else if (strcmp(name,"IntVecParam")==0) return IntVecParamEnum;
    1369          else stage=12;
    1370    }
    1371    if(stage==12){
    1372               if (strcmp(name,"Inputs")==0) return InputsEnum;
     1377              else if (strcmp(name,"Inputs")==0) return InputsEnum;
    13731378              else if (strcmp(name,"Internal")==0) return InternalEnum;
    13741379              else if (strcmp(name,"Intersect")==0) return IntersectEnum;
    … …  
    14851490              else if (strcmp(name,"Riftfront")==0) return RiftfrontEnum;
    14861491              else if (strcmp(name,"SamplingAnalysis")==0) return SamplingAnalysisEnum;
    1487               else if (strcmp(name,"SamplingSolution")==0) return SamplingSolutionEnum;
     1492         else stage=13;
     1493   }
     1494   if(stage==13){
     1495              if (strcmp(name,"SamplingSolution")==0) return SamplingSolutionEnum;
    14881496              else if (strcmp(name,"SIAApproximation")==0) return SIAApproximationEnum;
    14891497              else if (strcmp(name,"SMBautoregression")==0) return SMBautoregressionEnum;
    14901498              else if (strcmp(name,"SMBcomponents")==0) return SMBcomponentsEnum;
    14911499              else if (strcmp(name,"SMBd18opdd")==0) return SMBd18opddEnum;
    1492          else stage=13;
    1493    }
    1494    if(stage==13){
    1495               if (strcmp(name,"SMBforcing")==0) return SMBforcingEnum;
     1500              else if (strcmp(name,"SMBforcing")==0) return SMBforcingEnum;
    14961501              else if (strcmp(name,"SMBgcm")==0) return SMBgcmEnum;
    14971502              else if (strcmp(name,"SMBgemb")==0) return SMBgembEnum;

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

    r26660 r26676  
    7272                                        end
    7373                                end
     74                                if(contains(field,'WaterPressure'))
     75               mdname = structstoch.mdnames(find(strcmp(structstoch.fields,char(field))));
     76               if~(isequal(class(md.friction),char(mdname)))
     77                  error('stochasticforcing field %s is only implemented for default friction', char(field));
     78               end
     79               if(md.friction.coupling~=0 && md.friction.coupling~=1 && md.friction.coupling~=2)
     80                  error('stochasticforcing field %s is only implemented for cases md.friction.coupling 0 or 1 or 2', char(field));
     81               end
     82               if(any(md.friction.q==0))
     83                  error('stochasticforcing field %s requires non-zero q exponent',char(field));
     84               end
     85            end
    7486                                %Checking for specific dimensions
    7587                                if ~(strcmp(field,'SMBautoregression') || strcmp(field,'FrontalForcingsRignotAutoregression'))
    … …  
    185197                'DefaultCalving',...
    186198                'FloatingMeltRate',...
     199                'FrictionWaterPressure',...
    187200                'FrontalForcingsRignotAutoregression',...
    188201                'SMBautoregression'
    … …  
    191204                'calving',...
    192205                'basalforcings',...
     206                'friction',...
    193207                'frontalforcingsrignotautoregression',...
    194208                'SMBautoregression'

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

    r26674 r26676  
    7070                mdname = structstoch[field]
    7171                if (type(md.smb).__name__ != mdname):
    72                     raise TypeError('md.smb does not agree with stochasticforcing field {}'.format(mdname))
     72                    raise TypeError('md.smb does not agree with stochasticforcing field {}'.format(field))
    7373            if 'FrontalForcings' in field:
    7474                mdname = structstoch[field]
    7575                if (type(md.frontalforcings).__name__ != mdname):
    76                     raise TypeError('md.frontalforcings does not agree with stochasticforcing field {}'.format(mdname))
     76                    raise TypeError('md.frontalforcings does not agree with stochasticforcing field {}'.format(field))
    7777            if 'Calving' in field:
    7878                mdname = structstoch[field]
    7979                if (type(md.calving).__name__ != mdname):
    80                     raise TypeError('md.calving does not agree with stochasticforcing field {}'.format(mdname))
     80                    raise TypeError('md.calving does not agree with stochasticforcing field {}'.format(field))
    8181            if 'BasalforcingsFloatingice' in field:
    8282                mdname = structstoch[field]
    8383                if (type(md.basalforcings).__name__ != mdname):
    84                     raise TypeError('md.basalforcings does not agree with stochasticforcing field {}'.format(mdname))
     84                    raise TypeError('md.basalforcings does not agree with stochasticforcing field {}'.format(field))
     85            if 'WaterPressure' in field:
     86                mdname = structstoch[field]
     87                if (type(md.friction).__name__ != mdname):
     88                    raise TypeError('stochasticforcing field {} is only implemented for default friction'.format(field))
     89                if (md.friction.coupling!=0 and md.friction.coupling!=1 and md.friction.coupling!=2):
     90                    raise TypeError('stochasticforcing field {} is only implemented for cases md.friction.coupling 0 or 1 or 2'.format(field))
     91                if (np.any(md.friction.q==0):
     92                        raise TypeError('stochasticforcing field {} requires non-zero q exponent'.format(field))
     93
    8594            # Checking for specific dimensions
    8695            if not (field == 'SMBautoregression' or field == 'FrontalForcingsRignotAutoregression'):
    … …  
    173182        structure = {'DefaultCalving': 'calving',
    174183                     'FloatingMeltRate': 'basalforcings',
     184                     'FrictionWaterPressure': 'friction',
    175185                     'FrontalForcingsRignotAutoregression': 'frontalforcingsrignotautoregression',
    176186                     'SMBautoregression': 'SMBautoregression'}
Note: See TracChangeset for help on using the changeset viewer.