Changeset 26696


Ignore:
Timestamp:
12/02/21 08:09:53 (3 years ago)
Author:
aleahsommers
Message:

CHG: cleaned up SHAKTI code

File:
1 edited

Legend:

Unmodified
Added
Removed

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

    r26650 r26696  
    202202        IssmDouble Jdet;
    203203        IssmDouble* xyz_list = NULL;
    204         IssmDouble  gap,bed,thickness,head,g,rho_ice,rho_water,A,B,n,head_old;/***/
     204        IssmDouble  gap,bed,thickness,head,g,rho_ice,rho_water,A,B,n;
    205205
    206206        /*Fetch number of nodes and dof for this finite element*/
    … …  
    223223        element->FindParam(&dt,TimesteppingTimeStepEnum);
    224224
    225         /***Get all inputs and parameters*/
     225        /*Get all inputs and parameters*/
    226226        element->FindParam(&rho_water,MaterialsRhoFreshwaterEnum);
    227227        element->FindParam(&rho_ice,MaterialsRhoIceEnum);
    … …  
    233233        Input* head_input = element->GetInput(HydrologyHeadEnum);              _assert_(head_input);
    234234        Input* base_input = element->GetInput(BaseEnum);                      _assert_(base_input);
    235         Input* head_old_input = element->GetInput(HydrologyHeadOldEnum);        _assert_(head_old_input);
    236235
    237236        /* Start  looping on the number of gaussian points: */
    … …  
    248247                gap_input->GetInputValue(&gap,gauss);
    249248                head_input->GetInputValue(&head,gauss);
    250                 head_old_input->GetInputValue(&head_old,gauss);
    251 
    252                 /***Get ice A parameter*/
     249
     250                /*Get ice A parameter*/
    253251                B_input->GetInputValue(&B,gauss);
    254252                n_input->GetInputValue(&n,gauss);
    255253                A=pow(B,-n);
    256254
    257                 /***Get water and ice pressures*/
     255                /*Get water and ice pressures*/
    258256                IssmDouble pressure_ice   = rho_ice*g*thickness;    _assert_(pressure_ice>0.);
    259257                IssmDouble pressure_water = rho_water*g*(head-bed);
    260258                if(pressure_water>pressure_ice) pressure_water = pressure_ice;
    261259
    262                 /***Get water pressure from previous time step to use in lagged creep term*/
    263                 IssmDouble pressure_water_old = rho_water*g*(head_old-bed);
    264                 if(pressure_water_old>pressure_ice) pressure_water_old = pressure_ice;
    265260
    266261
    … …  
    293288   IssmDouble  PMPheat,dpressure_water[2],dbed[2];     
    294289        IssmDouble* xyz_list = NULL;
    295         IssmDouble dgapxx; /***/
    296         IssmDouble dgapyy; /***/
     290//        IssmDouble dgapxx; /***/
     291//      IssmDouble dgapyy; /***/
    297292
    298293        /*Fetch number of nodes and dof for this finite element*/
    … …  
    320315        Input* br_input             = element->GetInput(HydrologyBumpHeightEnum);        _assert_(br_input);
    321316   Input* headold_input        = element->GetInput(HydrologyHeadOldEnum);           _assert_(headold_input);
    322         Input* dgapxx_input         = element->GetInput(HydrologyGapHeightXXEnum); /***/
    323         Input* dgapyy_input         = element->GetInput(HydrologyGapHeightYYEnum); /***/
     317//      Input* dgapxx_input         = element->GetInput(HydrologyGapHeightXXEnum); /***/
     318//      Input* dgapyy_input         = element->GetInput(HydrologyGapHeightYYEnum); /***/
    324319
    325320        /*Get conductivity from inputs*/
    … …  
    351346                br_input->GetInputValue(&br,gauss);
    352347                headold_input->GetInputValue(&head_old,gauss);
    353                 dgapxx_input->GetInputValue(&dgapxx,gauss); /***/
    354                 dgapyy_input->GetInputValue(&dgapyy,gauss); /***/
     348//              dgapxx_input->GetInputValue(&dgapxx,gauss); /***/
     349//              dgapyy_input->GetInputValue(&dgapyy,gauss); /***/
    355350
    356351                /*Get ice A parameter*/
    … …  
    365360                 beta = 0.;
    366361
    367                 /* TEST Compute lc term */
    368 //              if(gap<br)
    369 //               lc = gap*(1-(br-gap)/br);
    370 //              else
    371                  lc = gap;
    372 //                lc=gap*tanh(pow(gap,3)/pow(br,3));
    373 //              lc=pow(gap,2)/lr;
    374 
    375362                /*Compute frictional heat flux*/
    376363                friction->GetAlpha2(&alpha2,gauss);
    … …  
    378365                frictionheat=alpha2*(vx*vx+vy*vy);
    379366
    380                 /* *** TEST take out frictional heat for large gap height *** */
     367                /* Take out frictional heat for large gap height */
    381368                if(gap>br)
    382369                 frictionheat=0.;
    … …  
    399386                _assert_(meltrate>0.);
    400387
    401 
    402 //              for(int i=0;i<numnodes;i++) pe->values[i]+=Jdet*gauss->weight*
    403 //               (
    404 //                meltrate*(1/rho_water-1/rho_ice)
    405 //                +A*pow(fabs(pressure_ice - pressure_water),n-1)*(pressure_ice - pressure_water)*gap
    406 //                -beta*sqrt(vx*vx+vy*vy)
    407 //                +ieb
    408 //                +storage*head_old/dt
    409 //                )*basis[i];     
    410 
    411                  /* With weighted creep term*/
    412 //                 for(int i=0;i<numnodes;i++) pe->values[i]+=Jdet*gauss->weight*
    413  //                 (
    414   //                 meltrate*(1/rho_water-1/rho_ice)
    415    //                +A*pow(fabs(pressure_ice - pressure_water),n-1)*(pressure_ice + rho_water*g*bed)*gap
    416     //               +(n-1)*A*pow(fabs(pressure_ice - pressure_water),n-1)*(rho_water*g*head)*gap
    417      //              -beta*sqrt(vx*vx+vy*vy)
    418       //             +ieb
    419        //            +storage*head_old/dt
    420         //           )*basis[i];
    421 
    422                  /* With weighted creep term */
    423388                  for(int i=0;i<numnodes;i++) pe->values[i]+=Jdet*gauss->weight*
    424389                   (
    … …  
    443408         //            )*basis[i];
    444409
    445 
    446410       
    447 //                 /*Test with linear creep rate*/
    448  //                for(int i=0;i<numnodes;i++) pe->values[i]+=Jdet*gauss->weight*
    449   //                (
    450    //                meltrate*(1/rho_water-1/rho_ice)
    451     //               +(pressure_ice - pressure_water)*gap/pow(10,13)
    452      //              -beta*sqrt(vx*vx+vy*vy)
    453       //             +ieb
    454        //            +storage*head_old/dt
    455         //           )*basis[i];
    456411        }
    457412        /*Clean up and return*/
    … …  
    597552        IssmDouble q = 0.;
    598553   IssmDouble channelization = 0.;
    599         IssmDouble dgapxx; /***/
    600         IssmDouble dgapyy; /***/
     554//      IssmDouble dgapxx; /***/
     555//      IssmDouble dgapyy; /***/
    601556
    602557        /*Retrieve all inputs and parameters*/
    … …  
    616571        Input* lr_input             = element->GetInput(HydrologyBumpSpacingEnum);       _assert_(lr_input);
    617572        Input* br_input             = element->GetInput(HydrologyBumpHeightEnum);        _assert_(br_input);
    618         Input* dgapxx_input         = element->GetInput(HydrologyGapHeightXXEnum); /***/
    619         Input* dgapyy_input         = element->GetInput(HydrologyGapHeightYYEnum); /***/
     573 //       Input* dgapxx_input         = element->GetInput(HydrologyGapHeightXXEnum); /***/
     574  //      Input* dgapyy_input         = element->GetInput(HydrologyGapHeightYYEnum); /***/
    620575
    621576        /*Get conductivity from inputs*/
    … …  
    643598                lr_input->GetInputValue(&lr,gauss);
    644599                br_input->GetInputValue(&br,gauss);
    645                 dgapxx_input->GetInputValue(&dgapxx,gauss); /***/
    646                 dgapyy_input->GetInputValue(&dgapyy,gauss); /***/
     600//              dgapxx_input->GetInputValue(&dgapxx,gauss); /***/
     601//              dgapyy_input->GetInputValue(&dgapyy,gauss); /***/
    647602
    648603                /*Get ice A parameter*/
    … …  
    657612                 beta = 0.;
    658613
    659                 /* TEST Compute lc term*/
    660 //              if(gap<br)
    661 //               lc = gap*(1-(br-gap)/br);
    662 //              else
    663                  lc = gap;
    664 //               lc=gap*tanh(pow(gap,3)/pow(br,3));
    665 //              lc = pow(gap,2)/lr;
    666 
    667614                /*Compute frictional heat flux*/
    668615                friction->GetAlpha2(&alpha2,gauss);
    … …  
    670617                frictionheat=alpha2*(vx*vx+vy*vy);
    671618
    672                 /* *** TEST take out frictional heat for large gap height *** */
     619                /* Take out frictional heat for large gap height */
    673620                if(gap>br)
    674621                 frictionheat=0.;
    … …  
    702649
    703650
    704                  /*TEST with linear creep rate*/
    705 //                 newgap += gauss->weight*Jdet*(gap+dt*(
    706  //                                        meltrate/rho_ice
    707   //                                       -(pressure_ice-pressure_water)*gap/pow(10,13)
    708    //                                      +beta*sqrt(vx*vx+vy*vy)
    709     //                                     ));
    710651                totalweights +=gauss->weight*Jdet;
    711652
Note: See TracChangeset for help on using the changeset viewer.