Changeset 28050


Ignore:
Timestamp:
01/17/24 17:08:12 (14 months ago)
Author:
sehrenfe
Message:

CHG: GlaDS- new input parameter rheology_B_base and new optional output HydrologySheetDischarge added to hydrology model

Location:
issm/trunk-jpl
Files:
8 edited

Legend:

Unmodified
Added
Removed

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

    r27937 r28050  
    151151        iomodel->FetchDataToInput(inputs,elements,"md.initialization.watercolumn",HydrologySheetThicknessEnum);
    152152        iomodel->FetchDataToInput(inputs,elements,"md.initialization.hydraulic_potential",HydraulicPotentialEnum);
     153        iomodel->FetchDataToInput(inputs,elements,"md.hydrology.rheology_B_base",HydrologyRheologyBBaseEnum);
     154
    153155        if(iomodel->domaintype==Domain2DhorizontalEnum){
    154156                iomodel->FetchDataToInput(inputs,elements,"md.initialization.vx",VxEnum);
    … …  
    249251        IssmDouble g         = element->FindParam(ConstantsGEnum);
    250252        IssmDouble e_v       = element->FindParam(HydrologyEnglacialVoidRatioEnum);
    251         Input* hr_input      = element->GetInput(HydrologyBumpHeightEnum);_assert_(hr_input);
     253        Input* hr_input  = element->GetInput(HydrologyBumpHeightEnum);       _assert_(hr_input);
    252254        Input* k_input   = element->GetInput(HydrologySheetConductivityEnum);_assert_(k_input);
    253         Input* phi_input = element->GetInput(HydraulicPotentialEnum);      _assert_(phi_input);
    254         Input* h_input   = element->GetInput(HydrologySheetThicknessEnum); _assert_(h_input);
    255         Input* H_input      = element->GetInput(ThicknessEnum); _assert_(H_input);
    256         Input* b_input      = element->GetInput(BedEnum); _assert_(b_input);
    257         Input* B_input      = element->GetInput(MaterialsRheologyBEnum);         _assert_(B_input);
    258         Input* n_input      = element->GetInput(MaterialsRheologyNEnum);         _assert_(n_input);
     255        Input* phi_input = element->GetInput(HydraulicPotentialEnum);        _assert_(phi_input);
     256        Input* h_input   = element->GetInput(HydrologySheetThicknessEnum);   _assert_(h_input);
     257        Input* H_input   = element->GetInput(ThicknessEnum);                _assert_(H_input);
     258        Input* b_input   = element->GetInput(BedEnum);                      _assert_(b_input);
     259        Input* B_input   = element->GetInput(HydrologyRheologyBBaseEnum);    _assert_(B_input);
     260        Input* n_input   = element->GetInput(MaterialsRheologyNEnum);        _assert_(n_input);
    259261
    260262        /* Start  looping on the number of gaussian points: */
    … …  
    275277                b_input->GetInputValue(&b,gauss);
    276278                H_input->GetInputValue(&H,gauss);
    277 
    278                 /*Hard code B*/
    279                 B = Cuffey(273.15-2);
    280279               
    281280                /*Get norm of gradient of hydraulic potential and make sure it is >0*/
    … …  
    364363        IssmDouble g         = element->FindParam(ConstantsGEnum);
    365364        IssmDouble e_v       = element->FindParam(HydrologyEnglacialVoidRatioEnum);
    366         Input* hr_input     = element->GetInput(HydrologyBumpHeightEnum);_assert_(hr_input);
    367         Input* h_input      = element->GetInput(HydrologySheetThicknessEnum);_assert_(h_input);
    368         Input* H_input      = element->GetInput(ThicknessEnum); _assert_(H_input);
    369         Input* b_input      = element->GetInput(BedEnum); _assert_(b_input);
    370         Input* G_input      = element->GetInput(BasalforcingsGeothermalfluxEnum);_assert_(G_input);
     365        Input* hr_input     = element->GetInput(HydrologyBumpHeightEnum);                _assert_(hr_input);
     366        Input* h_input      = element->GetInput(HydrologySheetThicknessEnum);            _assert_(h_input);
     367        Input* H_input      = element->GetInput(ThicknessEnum);                          _assert_(H_input);
     368        Input* b_input      = element->GetInput(BedEnum);                                _assert_(b_input);
     369        Input* G_input      = element->GetInput(BasalforcingsGeothermalfluxEnum);        _assert_(G_input);
    371370        Input* melt_input   = element->GetInput(BasalforcingsGroundediceMeltingRateEnum);_assert_(melt_input);
    372371        Input* RO_input     = NULL;
    373         Input* B_input      = element->GetInput(MaterialsRheologyBEnum);         _assert_(B_input);
    374         Input* n_input      = element->GetInput(MaterialsRheologyNEnum);         _assert_(n_input);
    375         Input* phiold_input = element->GetInput(HydraulicPotentialOldEnum);      _assert_(phiold_input);
    376         Input* phi_input    = element->GetInput(HydraulicPotentialEnum);         _assert_(phi_input);
     372        Input* B_input      = element->GetInput(HydrologyRheologyBBaseEnum);             _assert_(B_input);
     373        Input* n_input      = element->GetInput(MaterialsRheologyNEnum);                 _assert_(n_input);
     374        Input* phiold_input = element->GetInput(HydraulicPotentialOldEnum);              _assert_(phiold_input);
     375        Input* phi_input    = element->GetInput(HydraulicPotentialEnum);                 _assert_(phi_input);
    377376
    378377        /*Build friction element, needed later: */
    … …  
    397396                melt_input->GetInputValue(&melt,gauss);
    398397
    399                 /*Hard code B*/
    400                 B = Cuffey(273.15-2);
    401 
    402398                /*Get basal velocity*/
    403399                friction->GetBasalSlidingSpeeds(&vx, &vy ,gauss);
    … …  
    457453        element->InputUpdateFromSolutionOneDof(solution,HydraulicPotentialEnum);
    458454
    459         /*Compute Hydrology Vx and Vy for time stepping purposes (These inputs do not affect GlaDS)*/
     455        /*Compute Hydrology Vx and Vy for time stepping purposes, and Sheet Discharge as an optional output (These inputs do not affect GlaDS)*/
    460456
    461457        /*Intermediaries*/
    462458   IssmDouble  dphi[3],h,k,phi;
    463         IssmDouble  h_r;
     459        IssmDouble  h_r;
    464460        IssmDouble  oceanLS,iceLS;
    465461        IssmDouble* xyz_list = NULL;
    … …  
    468464        int numvertices = element->GetNumberOfVertices();
    469465
    470         /*Initialize new sheet thickness*/
     466        /*Initialize water sheet velocity and discharge*/
    471467        IssmDouble* vx = xNew<IssmDouble>(numvertices);
    472468        IssmDouble* vy = xNew<IssmDouble>(numvertices);
     469        IssmDouble* d = xNew<IssmDouble>(numvertices);
    473470
    474471        /*Set to 0 if inactive element*/
    … …  
    476473                for(int iv=0;iv<numvertices;iv++) vx[iv] = 0.;
    477474                for(int iv=0;iv<numvertices;iv++) vy[iv] = 0.;
     475                for(int iv=0;iv<numvertices;iv++) d[iv] = 0.;
    478476                element->AddInput(HydrologyWaterVxEnum,vx,P1DGEnum);
    479477                element->AddInput(HydrologyWaterVyEnum,vy,P1DGEnum);
     478                element->AddInput(HydrologySheetDischargeEnum,d,P1DGEnum);
    480479                xDelete<IssmDouble>(vx);
    481480                xDelete<IssmDouble>(vy);
     481                xDelete<IssmDouble>(d);
    482482                return;
    483483        }
    … …  
    499499        Input *iceLS_input   = element->GetInput(MaskIceLevelsetEnum);            _assert_(iceLS_input);
    500500
    501         /* Start  looping on the number of gaussian points: */
     501        /* Start looping on the number of gaussian points: */
    502502        Gauss* gauss=element->NewGauss();
    503503        for(int iv=0;iv<numvertices;iv++){
    … …  
    508508      phi_input->GetInputValue(&phi,gauss);
    509509      h_input->GetInputValue(&h,gauss);
    510         hr_input->GetInputValue(&h_r,gauss);
     510      hr_input->GetInputValue(&h_r,gauss);
    511511      k_input->GetInputValue(&k,gauss);
    512512                oceanLS_input->GetInputValue(&oceanLS,gauss);
    513513                iceLS_input->GetInputValue(&iceLS,gauss);
    514514
    515                 /*Set sheet thickness to zero if floating or no ice*/
     515                /*Set to zero if floating or no ice*/
    516516                if(oceanLS<0. || iceLS>0.){
    517517                        vx[iv] = 0.;
    518518         vy[iv] = 0.;
     519                        d[iv] = 0.;
    519520                }
    520521                else{
    … …  
    526527         /*If omega is zero, use standard model, otherwise transition model*/
    527528         IssmDouble nu = mu_water/rho_water;
    528         IssmDouble coeff;
    529         if(istransition==1 && omega>=AEPS){
    530                 IssmDouble hratio = fabs(h/h_r);
    531                 IssmDouble coarg = 1. + 4.*pow(hratio,3-2*alpha)*omega*k*pow(h,3)*normgradphi/nu;
    532                 coeff = nu/2./omega*pow(hratio,2*alpha-3) * (-1 + pow(coarg, 0.5))/normgradphi/max(AEPS,h);  // divide by h to get speed instead of discharge
    533         }
    534         else {
    535                 coeff = k*pow(h,alpha)*pow(normgradphi,beta-2.)/max(AEPS,h);  // divide by h to get speed instead of discharge
    536         }
    537 
    538 
    539                         vx[iv] = -coeff*dphi[0];
    540                         vy[iv] = -coeff*dphi[1];
     529                        IssmDouble coeff;
     530                        if(istransition==1 && omega>=AEPS){
     531                                IssmDouble hratio = fabs(h/h_r);
     532                                IssmDouble coarg = 1. + 4.*pow(hratio,3-2*alpha)*omega*k*pow(h,3)*normgradphi/nu;
     533                                coeff = nu/2./omega*pow(hratio,2*alpha-3) * (-1 + pow(coarg, 0.5))/normgradphi;  // coeff gives discharge; divide by h to get speed instead of discharge
     534                        }
     535                        else {
     536                        coeff = k*pow(h,alpha)*pow(normgradphi,beta-2.);  // coeff gives discharge; divide by h to get speed instead of discharge
     537                        }
     538
     539                        vx[iv] = -coeff/max(AEPS,h)*dphi[0];
     540                        vy[iv] = -coeff/max(AEPS,h)*dphi[1];
     541                               
     542                        d[iv] = coeff*normgradphi;
    541543                }
    542544        }
    … …  
    544546        element->AddInput(HydrologyWaterVxEnum,vx,P1DGEnum);
    545547        element->AddInput(HydrologyWaterVyEnum,vy,P1DGEnum);
     548        element->AddInput(HydrologySheetDischargeEnum,d,P1DGEnum);
    546549
    547550        /*Clean up and return*/
    … …  
    549552        xDelete<IssmDouble>(vx);
    550553        xDelete<IssmDouble>(vy);
     554        xDelete<IssmDouble>(d);
    551555        delete gauss;
    552556}/*}}}*/
    … …  
    646650        IssmDouble rho_water = element->FindParam(MaterialsRhoFreshwaterEnum);
    647651        IssmDouble g         = element->FindParam(ConstantsGEnum);
    648         Input* hr_input = element->GetInput(HydrologyBumpHeightEnum);_assert_(hr_input);
    649         Input* vx_input = element->GetInput(VxBaseEnum);_assert_(vx_input);
    650         Input* vy_input = element->GetInput(VyBaseEnum);_assert_(vy_input);
    651         Input* H_input  = element->GetInput(ThicknessEnum); _assert_(H_input);
    652         Input* b_input  = element->GetInput(BedEnum); _assert_(b_input);
    653         Input* hold_input  = element->GetInput(HydrologySheetThicknessOldEnum);_assert_(hold_input);
    654         Input* B_input  = element->GetInput(MaterialsRheologyBEnum);         _assert_(B_input);
    655         Input* n_input  = element->GetInput(MaterialsRheologyNEnum);         _assert_(n_input);
     652        Input* hr_input = element->GetInput(HydrologyBumpHeightEnum);         _assert_(hr_input);
     653        Input* vx_input = element->GetInput(VxBaseEnum);                      _assert_(vx_input);
     654        Input* vy_input = element->GetInput(VyBaseEnum);                      _assert_(vy_input);
     655        Input* H_input = element->GetInput(ThicknessEnum);                    _assert_(H_input);
     656        Input* b_input = element->GetInput(BedEnum);                          _assert_(b_input);
     657        Input* hold_input = element->GetInput(HydrologySheetThicknessOldEnum);_assert_(hold_input);
     658        Input* B_input = element->GetInput(HydrologyRheologyBBaseEnum);       _assert_(B_input);
     659        Input* n_input = element->GetInput(MaterialsRheologyNEnum);           _assert_(n_input);
    656660        Input* phi_input = element->GetInput(HydraulicPotentialEnum);         _assert_(phi_input);
    657         Input* oceanLS_input = element->GetInput(MaskOceanLevelsetEnum); _assert_(oceanLS_input);
    658         Input* iceLS_input = element->GetInput(MaskIceLevelsetEnum); _assert_(iceLS_input);
    659 
     661        Input* oceanLS_input = element->GetInput(MaskOceanLevelsetEnum);      _assert_(oceanLS_input);
     662        Input* iceLS_input = element->GetInput(MaskIceLevelsetEnum);          _assert_(iceLS_input);
    660663
    661664        /* Start  looping on the number of gaussian points: */
    … …  
    677680                iceLS_input->GetInputValue(&iceLS,gauss);
    678681
    679                 /*Hard code B*/
    680                 B = Cuffey(273.15-2);
    681 
    682682                /*Set sheet thickness to zero if floating or no ice*/
    683683                if(oceanLS<0. || iceLS>0.){

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

    r27936 r28050  
    354354        /*Initialize Element matrix and return if necessary*/
    355355        Tria*  tria=(Tria*)element;
    356         if(!tria->IsIceInElement()) return NULL;
     356        if(!tria->IsIceOnlyInElement()) return NULL;
    357357        _assert_(tria->FiniteElement()==P1Enum);
    358358        int index1=tria->GetVertexIndex(vertices[0]);
    … …  
    396396        Input* H_input      = element->GetInput(ThicknessEnum);                    _assert_(H_input);
    397397        Input* b_input      = element->GetInput(BedEnum);                          _assert_(b_input);
    398         Input* B_input      = element->GetInput(MaterialsRheologyBEnum);           _assert_(B_input);
     398        Input* B_input      = element->GetInput(HydrologyRheologyBBaseEnum);       _assert_(B_input);
    399399        Input* n_input      = element->GetInput(MaterialsRheologyNEnum);           _assert_(n_input);
    400400        Input* ks_input     = element->GetInput(HydrologySheetConductivityEnum);   _assert_(ks_input);
    401401        Input* kc_input     = element->GetInput(HydrologyChannelConductivityEnum); _assert_(kc_input);
    402         Input* hr_input     = element->GetInput(HydrologyBumpHeightEnum);_assert_(hr_input);
     402        Input* hr_input     = element->GetInput(HydrologyBumpHeightEnum);          _assert_(hr_input);
    403403        Input* phi_input    = element->GetInput(HydraulicPotentialEnum);           _assert_(phi_input);
    404404
    … …  
    432432                b_input->GetInputValue(&b,gauss);
    433433                H_input->GetInputValue(&H,gauss);
    434 
    435                 /*Hard code B*/
    436                 B = Cuffey(273.15-2);
    437434
    438435                /*Get values for a few potentials*/
    … …  
    493490                }
    494491
    495                 /*Closing rate term, see Gagliardini and Werder 2018 eq. A2 (v = v1*phi_i + v2(phi_{i+1}))*/
    496                 A=pow(B,-n);
    497                 v1 = 2./pow(n,n)*A*S*(pow(fabs(phi_0 - phi),n-1.)*( - n));
     492                /*Closing rate term*/
     493                /*See Gagliardini and Werder 2018 eq. A2 (v = v1*phi_i + v2(phi_{i+1}))*/
     494                A = pow(B,-n);
     495                if(phi_0-phi<0){
     496                        v1 = 0.;
     497                }
     498                else{
     499                        v1 = 2./pow(n,n)*A*S*(pow(fabs(phi_0-phi),n-1.)*( - n));
     500               
     501                }
     502
    498503                for(int i=0;i<numnodes;i++){
    499504                        for(int j=0;j<numnodes;j++){
    … …  
    512517        /*Initialize Element matrix and return if necessary*/
    513518        Tria* tria=(Tria*)element;
    514         if(!tria->IsIceInElement()) return NULL;
     519        if(!tria->IsIceOnlyInElement()) return NULL;
    515520        _assert_(tria->FiniteElement()==P1Enum);
    516521        int index1=tria->GetVertexIndex(vertices[0]);
    … …  
    550555        Input* H_input      = element->GetInput(ThicknessEnum);                    _assert_(H_input);
    551556        Input* b_input      = element->GetInput(BedEnum);                          _assert_(b_input);
    552         Input* B_input      = element->GetInput(MaterialsRheologyBEnum);           _assert_(B_input);
     557        Input* B_input      = element->GetInput(HydrologyRheologyBBaseEnum);       _assert_(B_input);
    553558        Input* n_input      = element->GetInput(MaterialsRheologyNEnum);           _assert_(n_input);
    554559        Input* ks_input     = element->GetInput(HydrologySheetConductivityEnum);   _assert_(ks_input);
    555560        Input* kc_input     = element->GetInput(HydrologyChannelConductivityEnum); _assert_(kc_input);
    556561        Input* phi_input    = element->GetInput(HydraulicPotentialEnum);           _assert_(phi_input);
    557         Input* hr_input     = element->GetInput(HydrologyBumpHeightEnum);_assert_(hr_input);
     562        Input* hr_input     = element->GetInput(HydrologyBumpHeightEnum);          _assert_(hr_input);
    558563
    559564        /*Get tangent vector*/
    … …  
    584589                hr_input->GetInputValue(&h_r,gauss);
    585590
    586                 /*Hard code B*/
    587                 B = Cuffey(273.15-2);
    588591               
    589592                /*Get values for a few potentials*/
    … …  
    625628                /*Compute closing rate*/
    626629                /*See Gagliardini and Werder 2018 eq. A2 (v = v2(phi_i) + v1*phi_{i+1})*/
    627                 A=pow(B,-n);
    628                 v2 = 2./pow(n,n)*A*this->S*(pow(fabs(phi_0 - phi),n-1.)*(phi_0 +(n-1.)*phi));
     630                A = pow(B,-n);
     631                if(phi_0-phi<0){
     632                        v2 = 0.;
     633                }
     634                else{
     635                        v2 = 2./pow(n,n)*A*this->S*(pow(fabs(phi_0 - phi),n-1.)*(phi_0 +(n-1.)*phi));
     636                }
    629637
    630638                for(int i=0;i<numnodes;i++){
    … …  
    648656        /*Initialize Element matrix and return if necessary*/
    649657        Tria*  tria=(Tria*)element;
    650         if(this->boundary){
     658        if(this->boundary || !tria->IsIceOnlyInElement()){
    651659                this->S = 0.;
    652660                return;
    … …  
    700708        Input* H_input      = element->GetInput(ThicknessEnum);                    _assert_(H_input);
    701709        Input* b_input      = element->GetInput(BedEnum);                          _assert_(b_input);
    702         Input* B_input      = element->GetInput(MaterialsRheologyBEnum);           _assert_(B_input);
     710        Input* B_input      = element->GetInput(HydrologyRheologyBBaseEnum);       _assert_(B_input);
    703711        Input* n_input      = element->GetInput(MaterialsRheologyNEnum);           _assert_(n_input);
    704712        Input* ks_input     = element->GetInput(HydrologySheetConductivityEnum);   _assert_(ks_input);
    705713        Input* kc_input     = element->GetInput(HydrologyChannelConductivityEnum); _assert_(kc_input);
    706714        Input* phi_input    = element->GetInput(HydraulicPotentialEnum);           _assert_(phi_input);
    707         Input* hr_input     = element->GetInput(HydrologyBumpHeightEnum);_assert_(hr_input);
     715        Input* hr_input     = element->GetInput(HydrologyBumpHeightEnum);          _assert_(hr_input);
    708716
    709717        /*Get tangent vector*/
    … …  
    728736
    729737
    730         /*Hard code B*/
    731         B = Cuffey(273.15-2);
    732 
    733738        /*Get values for a few potentials*/
    734739        phi_0   = rho_water*g*b + rho_ice*g*H;
    … …  
    754759
    755760        /*Ice rate factor*/
    756         A=pow(B,-n);
     761        A = pow(B,-n);
    757762
    758763        IssmDouble C = C_W*c_t*rho_water;
    … …  
    777782                                        + C*Qprime*pow(Snew,alpha_c-1.)*dPw
    778783                                        ) - 2./pow(n,n)*A*pow(fabs(N),n-1.)*N;
     784                if(N<0){
     785                        alpha = 1./(rho_ice*L)*(
     786               fabs(Qprime*pow(Snew,alpha_c-1.)*dphids)
     787               + C*Qprime*pow(Snew,alpha_c-1.)*dPw
     788               );
     789                }
    779790
    780791                IssmDouble beta = 1./(rho_ice*L)*( fabs(lc*qc*dphids) + C*fFactor*dPw );
    … …  
    786797                /*Constrain the cross section to be between 0 and 500 m^2*/
    787798                if(this->S<0.)   this->S = 0.;
    788                 if(this->S>100.) this->S = 100.;
    789                
    790                 /*Do not allow channels to grow in areas with no sheet thickness*/
    791                 if(H<200.) this->S = 0.;
     799                if(this->S>500.) this->S = 500.;
    792800
    793801                count++;

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

    r28036 r28050  
    899899        HydrologyNeumannfluxEnum,
    900900        HydrologyReynoldsEnum,
     901        HydrologyRheologyBBaseEnum,
    901902        HydrologySheetConductivityEnum,
     903        HydrologySheetDischargeEnum,
    902904        HydrologySheetThicknessEnum,
    903905        HydrologySheetThicknessOldEnum,

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

    r27937 r28050  
    1414                sheet_alpha               = NaN;
    1515                sheet_beta                = NaN;
     16                rheology_B_base           = NaN;
    1617
    1718                %Channels
    … …  
    8384                        md = checkfield(md,'fieldname','hydrology.sheet_alpha', 'numel', [1], '>', 0);
    8485                        md = checkfield(md,'fieldname','hydrology.sheet_beta', 'numel', [1], '>', 0);
    85 
     86                        md = checkfield(md,'fieldname','hydrology.rheology_B_base','size',[md.mesh.numberofvertices 1],'>=',0,'NaN',1,'Inf',1);
    8687                        %Channels
    8788                        md = checkfield(md,'fieldname','hydrology.ischannels','numel',[1],'values',[0 1]);
    … …  
    113114                        fielddisplay(self,'bump_height','typical bump height (h_r) [m]');
    114115                        fielddisplay(self,'omega','transition parameter (omega) []');
     116                        fielddisplay(self,'rheology_B_base','Ice rheology factor B at base of ice (B) [Pa s^(-1/3)]');
    115117                        disp(sprintf('      CHANNELS'));
    116118                        fielddisplay(self,'ischannels','Do we allow for channels? 1: yes, 0: no');
    … …  
    143145                        WriteData(fid,prefix,'object',self,'class','hydrology','fieldname','sheet_alpha','format','Double');
    144146                        WriteData(fid,prefix,'object',self,'class','hydrology','fieldname','sheet_beta','format','Double');
     147                        WriteData(fid,prefix,'object',self,'class','hydrology','fieldname','rheology_B_base','format','DoubleMat','mattype',1);
    145148
    146149                        %Channels

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

    r27937 r28050  
    2222        self.sheet_alpha = np.nan;
    2323        self.sheet_beta = np.nan;
     24        self.rheology_B_base = np.nan;
    2425
    2526        # Channels
    … …  
    5960        s += '{}\n'.format(fielddisplay(self, 'bump_height', 'typical bump height (h_r) [m]'))
    6061        s += '{}\n'.format(fielddisplay(self, 'omega', 'transition parameter (omega) []')) #TH
     62        s += '{}\n'.format(fielddisplay(self, 'rheology_B_base', 'ice rheology factor B at base of ice (B) [Pa s^(-1/3)]')) #SE
    6163        s = '\t--CHANNELS\n'
    6264        s += '{}\n'.format(fielddisplay(self, 'ischannels', 'Do we allow for channels? 1: yes, 0: no'))
    … …  
    129131        md = checkfield(md,'fieldname','hydrology.sheet_alpha', 'numel', [1], '>', 0);
    130132        md = checkfield(md,'fieldname','hydrology.sheet_beta', 'numel', [1], '>', 0);
     133        md = checkfield(md,'fieldname','hydrology.rheology_B_base', 'size', [md.mesh.numberofvertices], '>=', 0, 'np.nan', 1, 'Inf', 1)
    131134
    132135        # Channels
    … …  
    162165        WriteData(fid,prefix,'object',self,'class','hydrology','fieldname','sheet_alpha','format','Double');
    163166        WriteData(fid,prefix,'object',self,'class','hydrology','fieldname','sheet_beta','format','Double');
     167        WriteData(fid,prefix,'object',self,'class','hydrology','fieldname','rheology_B_base','format','DoubleMat', 'mattype', 1);
    164168
    165169        # Channels

  • issm/trunk-jpl/test/NightlyRun/test355.m

    r27447 r28050  
    5959md.hydrology.sheet_conductivity= 1e-3 * ones(md.mesh.numberofvertices,1);
    6060md.hydrology.channel_conductivity= 5.e-2 * ones(md.mesh.numberofvertices,1);
     61md.hydrology.rheology_B_base = cuffey(273.15 - 2)*ones(md.mesh.numberofvertices,1);
    6162
    6263% BCs for hydrology
    … …  
    7879        2e-10,2e-08,4e-07,...
    7980        3e-10,2e-08,4e-07,...
    80         4e-10,1e-08,4e-07};
     81        4e-10,2e-08,4e-07};
    8182field_values={...
    8283        md.results.TransientSolution(1).HydrologySheetThickness, ...

  • issm/trunk-jpl/test/NightlyRun/test355.py

    r27447 r28050  
    6868md.hydrology.sheet_conductivity = 1.e-3 * np.ones((md.mesh.numberofvertices))
    6969md.hydrology.channel_conductivity = 5.e-2 * np.ones((md.mesh.numberofvertices))
     70md.hydrology.rheology_B_base = 8.3788e7 * np.ones((md.mesh.numberofvertices))
    7071
    7172# BCs for hydrology
Note: See TracChangeset for help on using the changeset viewer.