Changeset 26490


Ignore:
Timestamp:
10/19/21 13:20:40 (4 years ago)
Author:
aleahsommers
Message:

CHG: changed the creep closure term in SHAKTI to reduce the nonlinearity and help convergence

File:
1 edited

Legend:

Unmodified
Added
Removed

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

    r26463 r26490  
    181181        IssmDouble Jdet;
    182182        IssmDouble* xyz_list = NULL;
     183        IssmDouble  gap,bed,thickness,head,g,rho_ice,rho_water,A,B,n,head_old;/***/
    183184
    184185        /*Fetch number of nodes and dof for this finite element*/
    … …  
    201202        element->FindParam(&dt,TimesteppingTimeStepEnum);
    202203
     204        /***Get all inputs and parameters*/
     205        element->FindParam(&rho_water,MaterialsRhoFreshwaterEnum);
     206        element->FindParam(&rho_ice,MaterialsRhoIceEnum);
     207        element->FindParam(&g,ConstantsGEnum);
     208        Input* B_input = element->GetInput(MaterialsRheologyBEnum);         _assert_(B_input);
     209        Input* n_input = element->GetInput(MaterialsRheologyNEnum);         _assert_(n_input);
     210        Input* gap_input = element->GetInput(HydrologyGapHeightEnum);         _assert_(gap_input);
     211        Input* thickness_input = element->GetInput(ThicknessEnum);                  _assert_(thickness_input);
     212        Input* head_input = element->GetInput(HydrologyHeadEnum);              _assert_(head_input);
     213        Input* base_input = element->GetInput(BaseEnum);                      _assert_(base_input);
     214        Input* head_old_input = element->GetInput(HydrologyHeadOldEnum);        _assert_(head_old_input);
     215
    203216        /* Start  looping on the number of gaussian points: */
    204217        Gauss* gauss=element->NewGauss(1);
    … …  
    209222                element->NodalFunctions(basis,gauss);
    210223
     224                /***/
     225                base_input->GetInputValue(&bed,gauss);
     226                thickness_input->GetInputValue(&thickness,gauss);
     227                gap_input->GetInputValue(&gap,gauss);
     228                head_input->GetInputValue(&head,gauss);
     229                head_old_input->GetInputValue(&head_old,gauss);
     230
     231                /***Get ice A parameter*/
     232                B_input->GetInputValue(&B,gauss);
     233                n_input->GetInputValue(&n,gauss);
     234                A=pow(B,-n);
     235
     236                /***Get water and ice pressures*/
     237                IssmDouble pressure_ice   = rho_ice*g*thickness;    _assert_(pressure_ice>0.);
     238                IssmDouble pressure_water = rho_water*g*(head-bed);
     239                if(pressure_water>pressure_ice) pressure_water = pressure_ice;
     240
     241                /***Get water pressure from previous time step to use in lagged creep term*/
     242                IssmDouble pressure_water_old = rho_water*g*(head_old-bed);
     243                if(pressure_water_old>pressure_ice) pressure_water_old = pressure_ice;
     244
     245
    211246                for(int i=0;i<numnodes;i++){
    212247                        for(int j=0;j<numnodes;j++){
    213248                                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])
    214                                   + gauss->weight*Jdet*storage/dt*basis[i]*basis[j];
     249                                  + gauss->weight*Jdet*storage/dt*basis[i]*basis[j]
     250                                        +gauss->weight*Jdet*A*(n)*(pow(fabs(pressure_ice-pressure_water),(n-1))*rho_water*g)*gap*basis[i]*basis[j];
    215251                        }
    216252                }
    … …  
    232268        IssmDouble  Jdet,meltrate,G,dh[2],B,A,n;
    233269        IssmDouble  gap,bed,thickness,head,ieb,head_old;
    234         IssmDouble  lr,br,vx,vy,beta;
     270        IssmDouble  lr,br,vx,vy,beta,lc;
    235271        IssmDouble  alpha2,frictionheat;
    236272   IssmDouble  PMPheat,dpressure_water[2],dbed[2];     
    237273        IssmDouble* xyz_list = NULL;
     274        IssmDouble dgap[2]; /***/
    238275
    239276        /*Fetch number of nodes and dof for this finite element*/
    … …  
    261298        Input* br_input             = element->GetInput(HydrologyBumpHeightEnum);        _assert_(br_input);
    262299   Input* headold_input        = element->GetInput(HydrologyHeadOldEnum);           _assert_(headold_input);
     300//      Input* dgapxx_input         = element->GetInput(HydrologyGapXXEnum); /***/
     301//      Input* dgapyy_input         = element->GetInput(HydrologyGapYYEnum); /***/
    263302
    264303        /*Get conductivity from inputs*/
    … …  
    289328                lr_input->GetInputValue(&lr,gauss);
    290329                br_input->GetInputValue(&br,gauss);
    291       headold_input->GetInputValue(&head_old,gauss);
     330                headold_input->GetInputValue(&head_old,gauss);
    292331
    293332                /*Get ice A parameter*/
    … …  
    302341                 beta = 0.;
    303342
     343                /* TEST Compute lc term */
     344//              if(gap<br)
     345//               lc = gap*(1-(br-gap)/br);
     346//              else
     347                 lc = gap;
     348//                lc=gap*tanh(pow(gap,3)/pow(br,3));
     349//              lc=pow(gap,2)/lr;
     350
    304351                /*Compute frictional heat flux*/
    305352                friction->GetAlpha2(&alpha2,gauss);
    … …  
    324371                _assert_(meltrate>0.);
    325372
    326                 for(int i=0;i<numnodes;i++) pe->values[i]+=Jdet*gauss->weight*
    327                  (
    328                   meltrate*(1/rho_water-1/rho_ice)
    329                   +A*pow(fabs(pressure_ice - pressure_water),n-1)*(pressure_ice - pressure_water)*gap
    330                   -beta*sqrt(vx*vx+vy*vy)
    331                   +ieb
    332                   +storage*head_old/dt
    333                   )*basis[i];           
     373
     374//              for(int i=0;i<numnodes;i++) pe->values[i]+=Jdet*gauss->weight*
     375//               (
     376//                meltrate*(1/rho_water-1/rho_ice)
     377//                +A*pow(fabs(pressure_ice - pressure_water),n-1)*(pressure_ice - pressure_water)*gap
     378//                -beta*sqrt(vx*vx+vy*vy)
     379//                +ieb
     380//                +storage*head_old/dt
     381//                )*basis[i];     
     382
     383                 /* With weighted creep term*/
     384                 for(int i=0;i<numnodes;i++) pe->values[i]+=Jdet*gauss->weight*
     385                  (
     386                   meltrate*(1/rho_water-1/rho_ice)
     387                   +A*pow(fabs(pressure_ice - pressure_water),n-1)*(pressure_ice + rho_water*g*bed)*gap
     388                   +(n-1)*A*pow(fabs(pressure_ice - pressure_water),n-1)*(rho_water*g*head)*gap
     389                   -beta*sqrt(vx*vx+vy*vy)
     390                   +ieb
     391                   +storage*head_old/dt
     392                   )*basis[i];
     393
     394       
     395//                 /*Test with linear creep rate*/
     396 //                for(int i=0;i<numnodes;i++) pe->values[i]+=Jdet*gauss->weight*
     397  //                (
     398   //                meltrate*(1/rho_water-1/rho_ice)
     399    //               +(pressure_ice - pressure_water)*gap/pow(10,13)
     400     //              -beta*sqrt(vx*vx+vy*vy)
     401      //             +ieb
     402       //            +storage*head_old/dt
     403        //           )*basis[i];
    334404        }
    335405        /*Clean up and return*/
    … …  
    382452                values[i]=solution[doflist[i]];
    383453
    384                 /*make sure that p_water<p_ice ->  h<rho_i H/rho_w + zb*/
    385                 if(values[i]>rho_ice*thickness[i]/rho_water+bed[i]){
    386                         values[i] = rho_ice*thickness[i]/rho_water+bed[i];
    387                 }
     454                /*overburden cap: make sure that p_water<p_ice ->  h<rho_i H/rho_w + zb*/
     455//              if(values[i]>rho_ice*thickness[i]/rho_water+bed[i]){
     456//                      values[i] = rho_ice*thickness[i]/rho_water+bed[i];
     457//              }
    388458
    389459                /*Make sure that negative pressure is not allowed*/
    390   //    if(values[i]<bed[i]){
    391         //              values[i] = bed[i];
    392         //      }
     460//      if(values[i]<bed[i]){
     461//                      values[i] = bed[i];
     462//              }
    393463
    394464                /*Under-relaxation*/
    … …  
    472542        IssmDouble  Jdet,meltrate,G,dh[2],B,A,n,dt;
    473543        IssmDouble  gap,bed,thickness,head;
    474         IssmDouble  lr,br,vx,vy,beta;
     544        IssmDouble  lr,br,vx,vy,beta,lc;
    475545        IssmDouble  alpha2,frictionheat;
    476546        IssmDouble* xyz_list = NULL;
    … …  
    532602                 beta = 0.;
    533603
     604                /* TEST Compute lc term*/
     605//              if(gap<br)
     606//               lc = gap*(1-(br-gap)/br);
     607//              else
     608                 lc = gap;
     609//               lc=gap*tanh(pow(gap,3)/pow(br,3));
     610//              lc = pow(gap,2)/lr;
     611
    534612                /*Compute frictional heat flux*/
    535613                friction->GetAlpha2(&alpha2,gauss);
    … …  
    555633                                        +beta*sqrt(vx*vx+vy*vy)
    556634                                        ));
     635
     636                /* TEST with gap height "diffusivity" for melting walls - not yet completed*/
     637//                newgap += gauss->weight*Jdet*(gap+dt*(
     638//                                         meltrate/rho_ice
     639//                                         -A*pow(fabs(pressure_ice-pressure_water),n-1)*(pressure_ice-pressure_water)*lc
     640//                                         +beta*sqrt(vx*vx+vy*vy)
     641//                                       +(meltrate*gap/rho_ice*(dgapxx+dgapyy))
     642//                                         ));
     643
     644
     645                 /*TEST with linear creep rate*/
     646//                 newgap += gauss->weight*Jdet*(gap+dt*(
     647 //                                        meltrate/rho_ice
     648  //                                       -(pressure_ice-pressure_water)*gap/pow(10,13)
     649   //                                      +beta*sqrt(vx*vx+vy*vy)
     650    //                                     ));
    557651                totalweights +=gauss->weight*Jdet;
    558652
Note: See TracChangeset for help on using the changeset viewer.