Ignore:
Timestamp:
12/04/14 21:24:24 (10 years ago)
Author:
seroussi
Message:

CHG: removed useless functions, call GetMaterialParameter instead

File:
1 edited

Legend:

Unmodified
Added
Removed

  • issm/trunk-jpl/src/c/classes/Materials/Matpar.cpp

    r18717 r18946  
    304304
    305305/*Matpar management: */
    306 void  Matpar::Configure(Elements* elementsin){/*{{{*/
     306void       Matpar::Configure(Elements* elementsin){/*{{{*/
    307307
    308308        /*nothing done yet!*/
    309309
     310}
     311/*}}}*/
     312void       Matpar::EnthalpyToThermal(IssmDouble* ptemperature,IssmDouble* pwaterfraction,IssmDouble enthalpy,IssmDouble pressure){/*{{{*/
     313
     314        /*Ouput*/
     315        IssmDouble temperature,waterfraction;
     316
     317        if(enthalpy<PureIceEnthalpy(pressure)){
     318                temperature=referencetemperature+enthalpy/heatcapacity;
     319                waterfraction=0.;
     320        }
     321        else{
     322                temperature=TMeltingPoint(pressure);
     323                waterfraction=(enthalpy-PureIceEnthalpy(pressure))/latentheat;
     324        }
     325
     326        /*Assign output pointers:*/
     327        *pwaterfraction=waterfraction;
     328        *ptemperature=temperature;
     329}
     330/*}}}*/
     331IssmDouble Matpar::GetEnthalpyDiffusionParameter(IssmDouble enthalpy,IssmDouble pressure){/*{{{*/
     332        if (enthalpy<PureIceEnthalpy(pressure))
     333                return thermalconductivity/heatcapacity;
     334        else
     335                return temperateiceconductivity/heatcapacity;
     336}
     337/*}}}*/
     338IssmDouble Matpar::GetEnthalpyDiffusionParameterVolume(int numvertices,IssmDouble* enthalpy,IssmDouble* pressure){/*{{{*/
     339
     340        int         iv;
     341        IssmDouble  lambda;                 // fraction of cold ice
     342        IssmDouble  kappa,kappa_c,kappa_t;  //enthalpy conductivities
     343        IssmDouble  Hc,Ht;
     344        IssmDouble* PIE   = xNew<IssmDouble>(numvertices);
     345        IssmDouble* dHpmp = xNew<IssmDouble>(numvertices);
     346
     347        for(iv=0; iv<numvertices; iv++){
     348                PIE[iv]=PureIceEnthalpy(pressure[iv]);
     349                dHpmp[iv]=enthalpy[iv]-PIE[iv];
     350        }
     351
     352        bool allequalsign=true;
     353        if(dHpmp[0]<0)
     354                for(iv=1; iv<numvertices;iv++) allequalsign=(allequalsign && (dHpmp[iv]<0));
     355        else
     356                for(iv=1; iv<numvertices;iv++) allequalsign=(allequalsign && (dHpmp[iv]>=0));
     357
     358        if(allequalsign){
     359                kappa=GetEnthalpyDiffusionParameter(enthalpy[0], pressure[0]);
     360        }
     361        else {
     362                /* return harmonic mean of thermal conductivities, weighted by fraction of cold/temperate ice,
     363                 cf Patankar 1980, pp44 */
     364                kappa_c=GetEnthalpyDiffusionParameter(PureIceEnthalpy(0.)-1.,0.);
     365                kappa_t=GetEnthalpyDiffusionParameter(PureIceEnthalpy(0.)+1.,0.);
     366                Hc=0.; Ht=0.;
     367                for(iv=0; iv<numvertices;iv++){
     368                        if(enthalpy[iv]<PIE[iv])
     369                         Hc+=(PIE[iv]-enthalpy[iv]);
     370                        else
     371                         Ht+=(enthalpy[iv]-PIE[iv]);
     372                }
     373                _assert_((Hc+Ht)>0.);
     374                lambda = Hc/(Hc+Ht);
     375                kappa  = 1./(lambda/kappa_c + (1.-lambda)/kappa_t);
     376        }
     377
     378        /*Clean up and return*/
     379        xDelete<IssmDouble>(PIE);
     380        xDelete<IssmDouble>(dHpmp);
     381        return kappa;
    310382}
    311383/*}}}*/
    … …  
    314386        switch(enum_in){
    315387                case MaterialsRhoIceEnum:                    return this->rho_ice;
    316                 case MaterialsRhoSeawaterEnum:                  return this->rho_water;
     388                case MaterialsRhoSeawaterEnum:               return this->rho_water;
    317389                case MaterialsRhoFreshwaterEnum:             return this->rho_freshwater;
    318390                case MaterialsMuWaterEnum:                   return this->mu_water;
    … …  
    344416                case MaterialsCompressionCoefEnum:                              return this->compression_coef;
    345417                case MaterialsTractionCoefEnum:                                 return this->traction_coef;
     418                case SurfaceforcingsDesfacEnum:              return this->desfac;
     419                case SurfaceforcingsS0pEnum:                 return this->s0p;
     420                case SurfaceforcingsS0tEnum:                 return this->s0t;
     421                case SurfaceforcingsRlapsEnum:               return this->rlaps;
     422                case SurfaceforcingsRlapslgmEnum:            return this->rlapslgm;
     423                case MaterialsLithosphereShearModulusEnum:   return this->lithosphere_shear_modulus;
     424                case MaterialsLithosphereDensityEnum:        return this->lithosphere_density;
     425                case MaterialsMantleDensityEnum:             return this->mantle_density;
     426                case MaterialsMantleShearModulusEnum:        return this->mantle_shear_modulus;
    346427                case ConstantsOmegaEnum:                                                        return this->omega;
    347428                case MaterialsTimeRelaxationStressEnum:         return this->time_relaxation_stress;
    … …  
    352433}
    353434/*}}}*/
    354 IssmDouble Matpar::GetBeta(){/*{{{*/
    355         return beta;
    356 }
    357 /*}}}*/
    358 IssmDouble Matpar::GetG(){/*{{{*/
    359         return g;
    360 }
    361 /*}}}*/
    362 IssmDouble Matpar::GetMeltingPoint(){/*{{{*/
    363         return meltingpoint;
    364 }
    365 /*}}}*/
    366 IssmDouble Matpar::GetRhoIce(){/*{{{*/
    367 
    368         return rho_ice;
    369 }
    370 /*}}}*/
    371 IssmDouble Matpar::GetRhoWater(){/*{{{*/
    372         return rho_water;
    373 }
    374 /*}}}*/
    375 IssmDouble Matpar::GetRhoFreshwater(){/*{{{*/
    376         return rho_freshwater;
    377 }
    378 /*}}}*/
    379 IssmDouble Matpar::GetDesFac(){/*{{{*/
    380         return desfac;
    381 }
    382 /*}}}*/
    383 IssmDouble Matpar::GetS0p(){/*{{{*/
    384         return s0p;
    385 }
    386 /*}}}*/
    387 IssmDouble Matpar::GetS0t(){/*{{{*/
    388         return s0t;
    389 }
    390 /*}}}*/
    391 IssmDouble Matpar::GetRlaps(){/*{{{*/
    392         return rlaps;
    393 }
    394 /*}}}*/
    395 IssmDouble Matpar::GetRlapslgm(){/*{{{*/
    396         return rlapslgm;
     435IssmDouble Matpar::PureIceEnthalpy(IssmDouble pressure){/*{{{*/
     436        return heatcapacity*(TMeltingPoint(pressure)-referencetemperature);
     437}
     438/*}}}*/
     439void       Matpar::ResetHooks(){/*{{{*/
     440
     441        //Nothing to be done
     442        return;
     443}
     444/*}}}*/
     445void       Matpar::ThermalToEnthalpy(IssmDouble * penthalpy,IssmDouble temperature,IssmDouble waterfraction,IssmDouble pressure){/*{{{*/
     446
     447        /*Ouput*/
     448        IssmDouble enthalpy;
     449
     450        if(temperature<TMeltingPoint(pressure)){
     451                enthalpy=heatcapacity*(temperature-referencetemperature);
     452        }
     453        else{
     454                enthalpy=PureIceEnthalpy(pressure)+latentheat*waterfraction;
     455        }
     456
     457        /*Assign output pointers:*/
     458        *penthalpy=enthalpy;
    397459}
    398460/*}}}*/
    … …  
    401463}
    402464/*}}}*/
    403 IssmDouble Matpar::PureIceEnthalpy(IssmDouble pressure){/*{{{*/
    404         return heatcapacity*(TMeltingPoint(pressure)-referencetemperature);
    405 }
    406 /*}}}*/
    407 IssmDouble Matpar::GetEnthalpyDiffusionParameter(IssmDouble enthalpy,IssmDouble pressure){/*{{{*/
    408         if (enthalpy<PureIceEnthalpy(pressure))
    409                 return thermalconductivity/heatcapacity;
    410         else
    411                 return temperateiceconductivity/heatcapacity;
    412 }
    413 /*}}}*/
    414 IssmDouble Matpar::GetEnthalpyDiffusionParameterVolume(int numvertices,IssmDouble* enthalpy,IssmDouble* pressure){/*{{{*/
    415 
    416         int         iv;
    417         IssmDouble  lambda;                 // fraction of cold ice
    418         IssmDouble  kappa,kappa_c,kappa_t;  //enthalpy conductivities
    419         IssmDouble  Hc,Ht;
    420         IssmDouble* PIE   = xNew<IssmDouble>(numvertices);
    421         IssmDouble* dHpmp = xNew<IssmDouble>(numvertices);
    422 
    423         for(iv=0; iv<numvertices; iv++){
    424                 PIE[iv]=PureIceEnthalpy(pressure[iv]);
    425                 dHpmp[iv]=enthalpy[iv]-PIE[iv];
    426         }
    427 
    428         bool allequalsign=true;
    429         if(dHpmp[0]<0)
    430                 for(iv=1; iv<numvertices;iv++) allequalsign=(allequalsign && (dHpmp[iv]<0));
    431         else
    432                 for(iv=1; iv<numvertices;iv++) allequalsign=(allequalsign && (dHpmp[iv]>=0));
    433 
    434         if(allequalsign){
    435                 kappa=GetEnthalpyDiffusionParameter(enthalpy[0], pressure[0]);
    436         }
    437         else {
    438                 /* return harmonic mean of thermal conductivities, weighted by fraction of cold/temperate ice,
    439                  cf Patankar 1980, pp44 */
    440                 kappa_c=GetEnthalpyDiffusionParameter(PureIceEnthalpy(0.)-1.,0.);
    441                 kappa_t=GetEnthalpyDiffusionParameter(PureIceEnthalpy(0.)+1.,0.);
    442                 Hc=0.; Ht=0.;
    443                 for(iv=0; iv<numvertices;iv++){
    444                         if(enthalpy[iv]<PIE[iv])
    445                          Hc+=(PIE[iv]-enthalpy[iv]);
    446                         else
    447                          Ht+=(enthalpy[iv]-PIE[iv]);
    448                 }
    449                 _assert_((Hc+Ht)>0.);
    450                 lambda = Hc/(Hc+Ht);
    451                 kappa  = 1./(lambda/kappa_c + (1.-lambda)/kappa_t);
    452         }
    453 
    454         /*Clean up and return*/
    455         xDelete<IssmDouble>(PIE);
    456         xDelete<IssmDouble>(dHpmp);
    457         return kappa;
    458 }
    459 /*}}}*/
    460 void Matpar::EnthalpyToThermal(IssmDouble* ptemperature,IssmDouble* pwaterfraction,IssmDouble enthalpy,IssmDouble pressure){/*{{{*/
    461 
    462         /*Ouput*/
    463         IssmDouble temperature,waterfraction;
    464 
    465         if(enthalpy<PureIceEnthalpy(pressure)){
    466                 temperature=referencetemperature+enthalpy/heatcapacity;
    467                 waterfraction=0.;
    468         }
    469         else{
    470                 temperature=TMeltingPoint(pressure);
    471                 waterfraction=(enthalpy-PureIceEnthalpy(pressure))/latentheat;
    472         }
    473 
    474         /*Assign output pointers:*/
    475         *pwaterfraction=waterfraction;
    476         *ptemperature=temperature;
    477 }
    478 /*}}}*/
    479 void Matpar::ThermalToEnthalpy(IssmDouble * penthalpy,IssmDouble temperature,IssmDouble waterfraction,IssmDouble pressure){/*{{{*/
    480 
    481         /*Ouput*/
    482         IssmDouble enthalpy;
    483 
    484         if(temperature<TMeltingPoint(pressure)){
    485                 enthalpy=heatcapacity*(temperature-referencetemperature);
    486         }
    487         else{
    488                 enthalpy=PureIceEnthalpy(pressure)+latentheat*waterfraction;
    489         }
    490 
    491         /*Assign output pointers:*/
    492         *penthalpy=enthalpy;
    493 }
    494 /*}}}*/
    495 IssmDouble Matpar::GetLithosphereShearModulus(){                 /*{{{*/
    496         return lithosphere_shear_modulus;                       
    497 }               
    498 /*}}}*/
    499 IssmDouble Matpar::GetLithosphereDensity(){              /*{{{*/
    500         return lithosphere_density;                     
    501 }               
    502 /*}}}*/
    503 IssmDouble Matpar::GetMantleDensity(){           /*{{{*/
    504         return mantle_density;                   
    505 }               
    506 /*}}}*/
    507 IssmDouble Matpar::GetMantleShearModulus(){              /*{{{*/
    508         return mantle_shear_modulus;                     
    509 }               
    510 /*}}}*/
    511 void  Matpar::ResetHooks(){/*{{{*/
    512 
    513         //Nothing to be done
    514         return;
    515 }
    516 /*}}}*/
Note: See TracChangeset for help on using the changeset viewer.