Changeset 16614


Ignore:
Timestamp:
11/05/13 08:37:14 (12 years ago)
Author:
jbondzio
Message:

BUG: ThermalSolution and SteadystateSolution give same results now. ADD: thermal conductivity at CTS is calculated using harmonic mean.

Location:
issm/trunk-jpl/src/c
Files:
4 edited

Legend:

Unmodified
Added
Removed

  • issm/trunk-jpl/src/c/classes/Elements/Penta.cpp

    r16600 r16614  
    37763776        IssmDouble epsvel=2.220446049250313e-16;
    37773777        IssmDouble heatcapacity,thermalconductivity,dt;
    3778         IssmDouble pressure,enthalpy;
     3778        IssmDouble pressure[NUMVERTICES],enthalpy[NUMVERTICES];
    37793779        IssmDouble latentheat,kappa;
    37803780        IssmDouble tau_parameter,diameter;
    … …  
    38153815
    38163816        /* Start  looping on the number of gaussian points: */
    3817         gauss=new GaussPenta(3,3);
     3817        gauss=new GaussPenta(2,2);
    38183818        for(int ig=gauss->begin();ig<gauss->end();ig++){
    38193819
    … …  
    38253825                GetBConduct(&B_conduct[0][0],&xyz_list[0][0],gauss);
    38263826
    3827                 enthalpy_input->GetInputValue(&enthalpy, gauss);
    3828                 pressure_input->GetInputValue(&pressure, gauss);
    3829                 kappa=matpar->GetEnthalpyDiffusionParameter(enthalpy,pressure); _assert_(kappa>0.);
     3827                //enthalpy_input->GetInputValue(&enthalpy, gauss);
     3828                //pressure_input->GetInputValue(&pressure, gauss);
     3829                GetInputListOnVertices(&enthalpy[0],EnthalpyPicardEnum);
     3830                GetInputListOnVertices(&pressure[0],PressureEnum);
     3831                kappa=matpar->GetEnthalpyDiffusionParameterVolume(enthalpy,pressure); _assert_(kappa>0.);  // TODO: switch to pointers
    38303832
    38313833                D_scalar_conduct=gauss->weight*Jdet*kappa/rho_ice;
    … …  
    42314233        IssmDouble thermalconductivity,kappa;
    42324234        IssmDouble viscosity,pressure;
    4233         IssmDouble enthalpy,enthalpypicard;
     4235        IssmDouble enthalpy;
     4236        IssmDouble enthalpy_v[NUMVERTICES], pressure_v[NUMVERTICES];
    42344237        IssmDouble tau_parameter,diameter;
    42354238        IssmDouble u,v,w;
    … …  
    42974300                        vy_input->GetInputValue(&v, gauss);
    42984301                        vz_input->GetInputValue(&w, gauss);
    4299                         pressure_input->GetInputValue(&pressure, gauss);
    4300                         enthalpypicard_input->GetInputValue(&enthalpypicard, gauss);
    4301                         kappa=matpar->GetEnthalpyDiffusionParameter(enthalpy,pressure);
     4302                        //pressure_input->GetInputValue(&pressure, gauss);
     4303                        //enthalpypicard_input->GetInputValue(&enthalpypicard, gauss);
     4304                        GetInputListOnVertices(&enthalpy_v[0],EnthalpyPicardEnum);
     4305                        GetInputListOnVertices(&pressure_v[0],PressureEnum);
     4306                        kappa=matpar->GetEnthalpyDiffusionParameterVolume(enthalpy_v,pressure_v);
    43024307                        kappa/=rho_ice;
    43034308                        tau_parameter=GetStabilizationParameter(u,v,w,diameter,kappa);
    … …  
    50215026                        else{
    50225027                                enthalpy_input->GetInputDerivativeValue(&d1enthalpy[0],&xyz_list[0][0],gauss);
    5023                                 kappa=matpar->GetEnthalpyDiffusionParameterVolume(enthalpy[iv],enthalpy[iv+NUMVERTICES2D], pressure[iv],pressure[iv+NUMVERTICES2D]); _assert_(kappa>0.);
     5028                                kappa=matpar->GetEnthalpyDiffusionParameterVolume(enthalpy, pressure); _assert_(kappa>0.);
    50245029                                for(i=0;i<3;i++) vec_heatflux[i]=-kappa*d1enthalpy[i];
    50255030                        }

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

    r16590 r16614  
    387387/*FUNCTION Matpar::GetEnthalpyDiffusionParameter{{{*/
    388388IssmDouble Matpar::GetEnthalpyDiffusionParameter(IssmDouble enthalpy,IssmDouble pressure){
    389         /*if (enthalpy<=PureIceEnthalpy(pressure))
     389        if (enthalpy<PureIceEnthalpy(pressure))
    390390                return thermalconductivity/heatcapacity;
    391391        else
    392                 return thermalconductivity/heatcapacity*pow(10.,-5);
    393         */
    394 
    395         IssmDouble eps=0.1*heatcapacity;
    396         IssmDouble hpmp=PureIceEnthalpy(pressure);
    397         IssmDouble kappa_c=thermalconductivity/heatcapacity;
    398         IssmDouble kappa_t=temperateiceconductivity/heatcapacity;
    399        
    400         if(enthalpy<=hpmp-eps)
    401                 return kappa_c;
    402         else if(enthalpy>=hpmp+eps)
    403                 return kappa_t;
     392                return temperateiceconductivity/heatcapacity;
     393}
     394/*}}}*/
     395/*FUNCTION Matpar::GetEnthalpyDiffusionParameterVolume{{{*/
     396IssmDouble Matpar::GetEnthalpyDiffusionParameterVolume(IssmDouble enthalpy[6], IssmDouble pressure[6]){
     397
     398        IssmDouble lambda; // fraction of cold ice
     399        IssmDouble kappa,kappa_c, kappa_t;
     400        IssmDouble Hc, Ht;
     401        IssmDouble PIE[6], dHpmp[6];
     402        int iv;
     403
     404        for (iv=0; iv<6; iv++){
     405                PIE[iv]=PureIceEnthalpy(pressure[iv]);
     406                dHpmp[iv]=enthalpy[iv]-PIE[iv];
     407        }
     408
     409        bool allequalsign=true;
     410        if(dHpmp[0]<0)
     411                for(iv=1; iv<6;iv++)
     412                        allequalsign=(allequalsign && (dHpmp[iv]<0));
     413        else
     414                for(iv=1; iv<6;iv++)
     415                        allequalsign=(allequalsign && (dHpmp[iv]>=0));
     416
     417        if(allequalsign){
     418                kappa=GetEnthalpyDiffusionParameter(enthalpy[0], pressure[0]);
     419        }
    404420        else {
    405                 IssmDouble xi=enthalpy-hpmp;
    406                 IssmDouble pi=3.141592653589793238462643;
    407                 return kappa_c + (kappa_t-kappa_c)*((xi+eps)/(2*eps) + sin(pi*xi/eps)/(2*pi));
    408         }
    409 
    410         //return 1./(1.+exp(-(enthalpy-(hpmp))/eps))*(kappa_t-kappa_c) + kappa_c;
    411 }
    412 /*}}}*/
    413 /*FUNCTION Matpar::GetEnthalpyDiffusionParameterVolume{{{*/
    414 IssmDouble Matpar::GetEnthalpyDiffusionParameterVolume(IssmDouble enthalpy0, IssmDouble enthalpy1, IssmDouble pressure0, IssmDouble pressure1){
    415         /*returns kappa depending on distribution of enthalpy over edge of element
    416                 lambda is the barycentric coordinate that solves H0+(H1-H0)*lambda=H_pureice.
    417                 it represents the fraction of the ice column which is temperate/cold like H0.
    418                 if lambda<0 or lambda>1, then the whole ice column is cold or temperate, respectively.
    419         */
    420         IssmDouble kappa, kappa0, kappa1;
    421         if (enthalpy0!=enthalpy1){
    422                 IssmDouble lambda=(PureIceEnthalpy(pressure0)-enthalpy0)/(enthalpy1-enthalpy0);
    423                 if ((lambda>=0.) && (lambda<=1.)){
    424                         kappa0=GetEnthalpyDiffusionParameter(enthalpy0,pressure0);
    425                         kappa1=GetEnthalpyDiffusionParameter(enthalpy1,pressure1);
    426                         kappa=lambda*kappa0+(1.-lambda)*kappa1;
     421                /* return harmonic mean of thermal conductivities, weighted by fraction of cold/temperate ice*/
     422                kappa_c=GetEnthalpyDiffusionParameter(PureIceEnthalpy(0.)-1.,0.);
     423                kappa_t=GetEnthalpyDiffusionParameter(PureIceEnthalpy(0.)+1.,0.);
     424                Hc=0.; Ht=0.;
     425                for(iv=0; iv<6;iv++){
     426                        if(enthalpy[iv]<PIE[iv])
     427                                Hc+=(PIE[iv]-enthalpy[iv]);
     428                        else
     429                                Ht+=(enthalpy[iv]-PIE[iv]);
    427430                }
    428         }
    429         else
    430                 kappa=GetEnthalpyDiffusionParameter(enthalpy0, pressure0);
     431                _assert_((Hc+Ht)>0.);
     432                lambda=Hc/(Hc+Ht);
     433                kappa=1./(lambda/kappa_c + (1.-lambda)/kappa_t);
     434        }
    431435        return kappa;
    432436}

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

    r16590 r16614  
    128128                IssmDouble PureIceEnthalpy(IssmDouble pressure);
    129129                IssmDouble GetEnthalpyDiffusionParameter(IssmDouble enthalpy,IssmDouble pressure);
    130                 IssmDouble GetEnthalpyDiffusionParameterVolume(IssmDouble enthalpy0, IssmDouble enthalpy1, IssmDouble pressure0, IssmDouble pressure1);
     130                IssmDouble GetEnthalpyDiffusionParameterVolume(IssmDouble enthalpy[6], IssmDouble pressure[6]);
    131131                IssmDouble GetLithosphereShearModulus();
    132132                IssmDouble GetLithosphereDensity();

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

    r16612 r16614  
    181181        MaterialsThermalExchangeVelocityEnum,
    182182        MaterialsThermalconductivityEnum,
    183   MaterialsTemperateiceconductivityEnum,
     183        MaterialsTemperateiceconductivityEnum,
    184184        MaterialsLithosphereShearModulusEnum,
    185185        MaterialsLithosphereDensityEnum,
Note: See TracChangeset for help on using the changeset viewer.