Changeset 16621


Ignore:
Timestamp:
11/06/13 10:02:04 (12 years ago)
Author:
jbondzio
Message:

CHG: Added harmonic mean for better convergence of steadystate and transient solver when handling a jump in enthalpy conductivity

Location:
issm/trunk-jpl
Files:
10 edited

Legend:

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

    r16614 r16621  
    37763776        IssmDouble epsvel=2.220446049250313e-16;
    37773777        IssmDouble heatcapacity,thermalconductivity,dt;
    3778         IssmDouble pressure[NUMVERTICES],enthalpy[NUMVERTICES];
     3778        IssmDouble enthalpy[NUMVERTICES], pressure[NUMVERTICES];
    37793779        IssmDouble latentheat,kappa;
    37803780        IssmDouble tau_parameter,diameter;
     
    38243824                /*Conduction: */ 
    38253825                GetBConduct(&B_conduct[0][0],&xyz_list[0][0],gauss);
    3826 
    3827                 //enthalpy_input->GetInputValue(&enthalpy, gauss);
    3828                 //pressure_input->GetInputValue(&pressure, gauss);
    3829                 GetInputListOnVertices(&enthalpy[0],EnthalpyPicardEnum);
     3826                GetInputListOnVertices(&enthalpy[0],EnthalpyEnum);
    38303827                GetInputListOnVertices(&pressure[0],PressureEnum);
    3831                 kappa=matpar->GetEnthalpyDiffusionParameterVolume(enthalpy,pressure); _assert_(kappa>0.);  // TODO: switch to pointers
     3828                kappa=matpar->GetEnthalpyDiffusionParameterVolume(enthalpy,pressure); _assert_(kappa>0.);
    38323829
    38333830                D_scalar_conduct=gauss->weight*Jdet*kappa/rho_ice;
     
    42324229        IssmDouble rho_ice,heatcapacity;
    42334230        IssmDouble thermalconductivity,kappa;
    4234         IssmDouble viscosity,pressure;
    4235         IssmDouble enthalpy;
    4236         IssmDouble enthalpy_v[NUMVERTICES], pressure_v[NUMVERTICES];
     4231        IssmDouble viscosity,pressure,enthalpy;
     4232        IssmDouble enthalpypicard[NUMVERTICES], pressure_p[NUMVERTICES];
    42374233        IssmDouble tau_parameter,diameter;
    42384234        IssmDouble u,v,w;
     
    43004296                        vy_input->GetInputValue(&v, gauss);
    43014297                        vz_input->GetInputValue(&w, gauss);
    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);
     4298                        GetInputListOnVertices(&enthalpypicard[0],EnthalpyPicardEnum);
     4299                        GetInputListOnVertices(&pressure_p[0],PressureEnum);
     4300                        kappa=matpar->GetEnthalpyDiffusionParameterVolume(enthalpypicard,pressure_p);
    43074301                        kappa/=rho_ice;
    43084302                        tau_parameter=GetStabilizationParameter(u,v,w,diameter,kappa);
  • issm/trunk-jpl/src/c/classes/Materials/Matpar.cpp

    r16614 r16621  
    397397
    398398        IssmDouble lambda; // fraction of cold ice
    399         IssmDouble kappa,kappa_c, kappa_t;
    400         IssmDouble Hc, Ht;
    401         IssmDouble PIE[6], dHpmp[6];
     399        IssmDouble kappa,kappa_c,kappa_t;  //enthalpy conductivities
     400        IssmDouble Hc,Ht;
     401        IssmDouble PIE[6],dHpmp[6];
    402402        int iv;
    403403
     
    419419        }
    420420        else {
    421                 /* return harmonic mean of thermal conductivities, weighted by fraction of cold/temperate ice*/
     421                /* return harmonic mean of thermal conductivities, weighted by fraction of cold/temperate ice,
     422                 cf Patankar 1980, pp44 */
    422423                kappa_c=GetEnthalpyDiffusionParameter(PureIceEnthalpy(0.)-1.,0.);
    423424                kappa_t=GetEnthalpyDiffusionParameter(PureIceEnthalpy(0.)+1.,0.);
Note: See TracChangeset for help on using the changeset viewer.