Changeset 21526


Ignore:
Timestamp:
02/07/17 13:22:30 (8 years ago)
Author:
aleahsommers
Message:

NEW:Added englacial storage term to subglacial hydrology model

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

Legend:

Unmodified
Added
Removed

  • TabularUnified issm/trunk-jpl/src/c/analyses/HydrologySommersAnalysis.cpp

    r21506 r21526  
    148148        parameters->AddObject(iomodel->CopyConstantObject("md.friction.law",FrictionLawEnum));
    149149   parameters->AddObject(iomodel->CopyConstantObject("md.hydrology.relaxation",HydrologyRelaxationEnum));
     150        parameters->AddObject(iomodel->CopyConstantObject("md.hydrology.storage",HydrologyStorageEnum));
    150151}/*}}}*/
    151152
    … …  
    173174        ElementMatrix* Ke     = element->NewElementMatrix();
    174175        IssmDouble*    dbasis = xNew<IssmDouble>(2*numnodes);
     176        IssmDouble*    basis  = xNew<IssmDouble>(numnodes);
    175177
    176178        /*Retrieve all inputs and parameters*/
    … …  
    179181        /*Get conductivity from inputs*/
    180182        IssmDouble conductivity = GetConductivity(element);
     183
     184        /*Get englacial storage coefficient*/
     185        IssmDouble storage,dt;
     186        element->FindParam(&storage,HydrologyStorageEnum);
     187        element->FindParam(&dt,TimesteppingTimeStepEnum);
    181188
    182189        /* Start  looping on the number of gaussian points: */
    … …  
    187194                element->JacobianDeterminant(&Jdet,xyz_list,gauss);
    188195                element->NodalFunctionsDerivatives(dbasis,xyz_list,gauss);
     196                element->NodalFunctions(basis,gauss);
    189197
    190198                for(int i=0;i<numnodes;i++){
    191199                        for(int j=0;j<numnodes;j++){
    192                                 Ke->values[i*numnodes+j] += conductivity*gauss->weight*Jdet*(dbasis[0*numnodes+i]*dbasis[0*numnodes+j] + dbasis[1*numnodes+i]*dbasis[1*numnodes+j]);
     200                                Ke->values[i*numnodes+j] += conductivity*gauss->weight*Jdet*(dbasis[0*numnodes+i]*dbasis[0*numnodes+j] + dbasis[1*numnodes+i]*dbasis[1*numnodes+j])
     201                                  + gauss->weight*Jdet*storage/dt*basis[i]*basis[j];
    193202                        }
    194203                }
    … …  
    244253        IssmDouble conductivity = GetConductivity(element);
    245254
     255        /*Get englacial storage coefficient*/
     256        IssmDouble storage,dt;
     257   element->FindParam(&storage,HydrologyStorageEnum);
     258   element->FindParam(&dt,TimesteppingTimeStepEnum);
     259
    246260        /*Build friction element, needed later: */
    247261        Friction* friction=new Friction(element,2);
    … …  
    308322                  -beta*sqrt(vx*vx+vy*vy)
    309323                  +ieb
     324                  +storage*head_old/dt
    310325                  )*basis[i];           
    311326        }

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

    r21518 r21526  
    176176   HydrologyRelaxationEnum,
    177177        HydrologyBasalFluxEnum,
     178        HydrologyStorageEnum,
    178179        InversionControlParametersEnum,
    179180        InversionControlScalingFactorsEnum,

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

    r21518 r21526  
    182182                case HydrologyRelaxationEnum : return "HydrologyRelaxation";
    183183                case HydrologyBasalFluxEnum : return "HydrologyBasalFlux";
     184                case HydrologyStorageEnum : return "HydrologyStorage";
    184185                case InversionControlParametersEnum : return "InversionControlParameters";
    185186                case InversionControlScalingFactorsEnum : return "InversionControlScalingFactors";

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

    r21518 r21526  
    185185              else if (strcmp(name,"HydrologyRelaxation")==0) return HydrologyRelaxationEnum;
    186186              else if (strcmp(name,"HydrologyBasalFlux")==0) return HydrologyBasalFluxEnum;
     187              else if (strcmp(name,"HydrologyStorage")==0) return HydrologyStorageEnum;
    187188              else if (strcmp(name,"InversionControlParameters")==0) return InversionControlParametersEnum;
    188189              else if (strcmp(name,"InversionControlScalingFactors")==0) return InversionControlScalingFactorsEnum;
    … …  
    259260              else if (strcmp(name,"CalvingMinthickness")==0) return CalvingMinthicknessEnum;
    260261              else if (strcmp(name,"DefaultCalving")==0) return DefaultCalvingEnum;
    261               else if (strcmp(name,"CalvinglevermannCoeff")==0) return CalvinglevermannCoeffEnum;
    262262         else stage=3;
    263263   }
    264264   if(stage==3){
    265               if (strcmp(name,"CalvinglevermannMeltingrate")==0) return CalvinglevermannMeltingrateEnum;
     265              if (strcmp(name,"CalvinglevermannCoeff")==0) return CalvinglevermannCoeffEnum;
     266              else if (strcmp(name,"CalvinglevermannMeltingrate")==0) return CalvinglevermannMeltingrateEnum;
    266267              else if (strcmp(name,"CalvingdevCoeff")==0) return CalvingdevCoeffEnum;
    267268              else if (strcmp(name,"Calvingratex")==0) return CalvingratexEnum;
    … …  
    382383              else if (strcmp(name,"SmbECini")==0) return SmbECiniEnum;
    383384              else if (strcmp(name,"SmbWini")==0) return SmbWiniEnum;
    384               else if (strcmp(name,"SmbAini")==0) return SmbAiniEnum;
    385385         else stage=4;
    386386   }
    387387   if(stage==4){
    388               if (strcmp(name,"SmbTini")==0) return SmbTiniEnum;
     388              if (strcmp(name,"SmbAini")==0) return SmbAiniEnum;
     389              else if (strcmp(name,"SmbTini")==0) return SmbTiniEnum;
    389390              else if (strcmp(name,"SmbSizeini")==0) return SmbSizeiniEnum;
    390391              else if (strcmp(name,"SMBforcing")==0) return SMBforcingEnum;
    … …  
    505506              else if (strcmp(name,"SurfaceSlopeY")==0) return SurfaceSlopeYEnum;
    506507              else if (strcmp(name,"Temperature")==0) return TemperatureEnum;
    507               else if (strcmp(name,"TemperaturePicard")==0) return TemperaturePicardEnum;
    508508         else stage=5;
    509509   }
    510510   if(stage==5){
    511               if (strcmp(name,"TemperaturePDD")==0) return TemperaturePDDEnum;
     511              if (strcmp(name,"TemperaturePicard")==0) return TemperaturePicardEnum;
     512              else if (strcmp(name,"TemperaturePDD")==0) return TemperaturePDDEnum;
    512513              else if (strcmp(name,"ThicknessAbsMisfit")==0) return ThicknessAbsMisfitEnum;
    513514              else if (strcmp(name,"SurfaceAbsMisfit")==0) return SurfaceAbsMisfitEnum;
    … …  
    628629              else if (strcmp(name,"Outputdefinition38")==0) return Outputdefinition38Enum;
    629630              else if (strcmp(name,"Outputdefinition39")==0) return Outputdefinition39Enum;
    630               else if (strcmp(name,"Outputdefinition40")==0) return Outputdefinition40Enum;
    631631         else stage=6;
    632632   }
    633633   if(stage==6){
    634               if (strcmp(name,"Outputdefinition41")==0) return Outputdefinition41Enum;
     634              if (strcmp(name,"Outputdefinition40")==0) return Outputdefinition40Enum;
     635              else if (strcmp(name,"Outputdefinition41")==0) return Outputdefinition41Enum;
    635636              else if (strcmp(name,"Outputdefinition42")==0) return Outputdefinition42Enum;
    636637              else if (strcmp(name,"Outputdefinition43")==0) return Outputdefinition43Enum;
    … …  
    751752              else if (strcmp(name,"Mpi")==0) return MpiEnum;
    752753              else if (strcmp(name,"Mumps")==0) return MumpsEnum;
    753               else if (strcmp(name,"Gsl")==0) return GslEnum;
    754754         else stage=7;
    755755   }
    756756   if(stage==7){
    757               if (strcmp(name,"Cuffey")==0) return CuffeyEnum;
     757              if (strcmp(name,"Gsl")==0) return GslEnum;
     758              else if (strcmp(name,"Cuffey")==0) return CuffeyEnum;
    758759              else if (strcmp(name,"BuddJacka")==0) return BuddJackaEnum;
    759760              else if (strcmp(name,"CuffeyTemperate")==0) return CuffeyTemperateEnum;
    … …  
    874875              else if (strcmp(name,"AdjointHorizAnalysis")==0) return AdjointHorizAnalysisEnum;
    875876              else if (strcmp(name,"DefaultAnalysis")==0) return DefaultAnalysisEnum;
    876               else if (strcmp(name,"BalancethicknessAnalysis")==0) return BalancethicknessAnalysisEnum;
    877877         else stage=8;
    878878   }
    879879   if(stage==8){
    880               if (strcmp(name,"BalancethicknessSolution")==0) return BalancethicknessSolutionEnum;
     880              if (strcmp(name,"BalancethicknessAnalysis")==0) return BalancethicknessAnalysisEnum;
     881              else if (strcmp(name,"BalancethicknessSolution")==0) return BalancethicknessSolutionEnum;
    881882              else if (strcmp(name,"Balancethickness2Analysis")==0) return Balancethickness2AnalysisEnum;
    882883              else if (strcmp(name,"Balancethickness2Solution")==0) return Balancethickness2SolutionEnum;
    … …  
    997998              else if (strcmp(name,"Materials")==0) return MaterialsEnum;
    998999              else if (strcmp(name,"Nodes")==0) return NodesEnum;
    999               else if (strcmp(name,"Contours")==0) return ContoursEnum;
    10001000         else stage=9;
    10011001   }
    10021002   if(stage==9){
    1003               if (strcmp(name,"Parameters")==0) return ParametersEnum;
     1003              if (strcmp(name,"Contours")==0) return ContoursEnum;
     1004              else if (strcmp(name,"Parameters")==0) return ParametersEnum;
    10041005              else if (strcmp(name,"Vertices")==0) return VerticesEnum;
    10051006              else if (strcmp(name,"Results")==0) return ResultsEnum;

  • TabularUnified issm/trunk-jpl/src/m/classes/hydrologysommers.m

    r21463 r21526  
    1616                neumannflux     = NaN;
    1717                relaxation      = 0;
     18                storage         = 0;
    1819        end
    1920        methods
    … …  
    3334              % Set under-relaxation parameter to be 1 (no under-relaxation of nonlinear iteration)     
    3435                        self.relaxation=1;
     36                        self.storage=0;
    3537                end % }}}
    3638                function md = checkconsistency(self,md,solution,analyses) % {{{
    … …  
    5153                        md = checkfield(md,'fieldname','hydrology.spchead','size',[md.mesh.numberofvertices 1]);       
    5254         md = checkfield(md,'fieldname','hydrology.relaxation','>=',0);
     55                        md = checkfield(md,'fieldname','hydrology.storage','>=',0);
    5356                end % }}}
    5457                function disp(self) % {{{
    … …  
    6467                        fielddisplay(self,'spchead','water head constraints (NaN means no constraint) (m)');
    6568                        fielddisplay(self,'relaxation','under-relaxation coefficient for nonlinear iteration');
     69                        fielddisplay(self,'storage','englacial storage coefficient (void ratio)');
    6670                end % }}}
    6771                function marshall(self,prefix,md,fid) % {{{
    … …  
    8084                        WriteData(fid,prefix,'object',self,'class','hydrology','fieldname','spchead','format','DoubleMat','mattype',1);
    8185         WriteData(fid,prefix,'object',self,'class','hydrology','fieldname','relaxation','format','Double');
     86                        WriteData(fid,prefix,'object',self,'class','hydrology','fieldname','storage','format','Double');
    8287                end % }}}
    8388        end
Note: See TracChangeset for help on using the changeset viewer.