Changeset 17027


Ignore:
Timestamp:
12/19/13 04:51:54 (11 years ago)
Author:
jbondzio
Message:

CHG: moved ComputeBasalMeltingrate and UpdateBasalConstraints to EnthalpyAnalysis. Minor fixes to the LliboutryDuval rheology.

Location:
issm/trunk-jpl/src
Files:
6 edited

Legend:

Unmodified
Added
Removed

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

    r17015 r17027  
    246246
    247247        /*Enthalpy diffusion parameter*/
    248         IssmDouble kappa=this->EnthalpyDiffusionParameterVolume(element,EnthalpyPicardEnum); _assert_(kappa>0.);
     248        IssmDouble kappa=this->EnthalpyDiffusionParameterVolume(element,EnthalpyPicardEnum); _assert_(kappa>=0.);
    249249
    250250        /* Start  looping on the number of gaussian points: */
    … …  
    443443        if(stabilization==2){
    444444                diameter=element->MinEdgeLength(xyz_list);
    445                 kappa=this->EnthalpyDiffusionParameterVolume(element,EnthalpyPicardEnum); _assert_(kappa>0.);
     445                kappa=this->EnthalpyDiffusionParameterVolume(element,EnthalpyPicardEnum); _assert_(kappa>=0.);
    446446        }
    447447
    … …  
    829829        }
    830830
    831         /*drain excess water fraction: */
    832         //for(int i=0;i<femmodel->elements->Size();i++){
    833         //      element=dynamic_cast<Element*>(femmodel->elements->GetObjectByOffset(i));
    834         //      element->DrainWaterfraction();
    835         //}
    836 
    837831        /*Update basal dirichlet BCs for enthalpy: */
    838832        for(int i=0;i<femmodel->elements->Size();i++){
    … …  
    841835        }
    842836}/*}}}*/
    843 
    844 
    845 
    846 
    847 
    848837
    849838void EnthalpyAnalysis::ComputeBasalMeltingrate(Element* element){/*{{{*/
    … …  
    878867        Input* vy_input               = element->GetInput(VyEnum);                          _assert_(vy_input);
    879868        Input* vz_input               = element->GetInput(VzEnum);                          _assert_(vz_input);
    880         IssmDouble kappa=EnthalpyDiffusionParameterVolume(element,EnthalpyEnum);     _assert_(kappa>0.);
     869        IssmDouble kappa=EnthalpyDiffusionParameterVolume(element,EnthalpyEnum);     _assert_(kappa>=0.);
    881870        element->NormalBase(&normal_base[0],xyz_list_base);
    882871        element->VerticalSegmentIndices(&pairindices,&numsegments);
    … …  
    990979                        basalmeltingrate[vertexdown]=meltingrate_enthalpy[is];
    991980                        watercolumn[vertexdown]+=meltingrate_enthalpy[is];
    992                 }              
     981                }       
    993982                _assert_(watercolumn[vertexdown]>=0.);
    994983        }
    … …  
    1002991        delete gauss;
    1003992        delete friction;
    1004         xDelete<IssmDouble>(enthalpy);
    1005         xDelete<IssmDouble>(pressure);
    1006         xDelete<IssmDouble>(watercolumn);
    1007         xDelete<IssmDouble>(basalmeltingrate);
     993        xDelete<IssmDouble>(enthalpy);
     994        xDelete<IssmDouble>(pressure);
     995        xDelete<IssmDouble>(watercolumn);
     996        xDelete<IssmDouble>(basalmeltingrate);
    1008997        xDelete<IssmDouble>(meltingrate_enthalpy);
    1009998        xDelete<IssmDouble>(heating);
    1010999        xDelete<IssmDouble>(drainrate_column);
    10111000        xDelete<IssmDouble>(drainrate_element);
     1001        xDelete<IssmDouble>(xyz_list);
    10121002}/*}}}*/
    10131003
    … …  
    10701060}/*}}}*/
    10711061
    1072 
    1073 
    1074 
    1075 
    10761062void EnthalpyAnalysis::UpdateBasalConstraints(Element* element){/*{{{*/
    10771063
    1078         /* /\*Intermediary*\/ */
    1079         /* bool        isdynamicbasalspc,setspc; */
    1080         /* int         numindices, numindicesup; */
    1081         /* IssmDouble  pressure, pressureup; */
    1082         /* IssmDouble  h_pmp, enthalpy, enthalpyup; */
    1083         /* IssmDouble  watercolumn; */
    1084         /* int        *indices = NULL, *indicesup = NULL; */
     1064        /*Intermediary*/
     1065        bool        isdynamicbasalspc,setspc;
     1066        int         numindices, numindicesup;
     1067        IssmDouble  pressure, pressureup;
     1068        IssmDouble  h_pmp, enthalpy, enthalpyup;
     1069        IssmDouble  watercolumn;
     1070        int        *indices = NULL, *indicesup = NULL;
     1071        Node*       node = NULL;
    10851072       
    1086         /* /\* Only update Constraints at the base of grounded ice*\/ */
    1087         /* if(!IsOnBed() || IsFloating()) return; */
    1088 
    1089         /* /\*Check wether dynamic basal boundary conditions are activated *\/ */
    1090         /* element->FindParam(&isdynamicbasalspc,ThermalIsdynamicbasalspcEnum); */
    1091         /* if(!isdynamicbasalspc) return; */
    1092 
    1093         /* /\*Fetch indices of basal & surface nodes for this finite element*\/ */
    1094         /* // BasalNodeIndices(&numindices,&indices,this->element_type); */
    1095         /* // SurfaceNodeIndices(&numindicesup,&indicesup,this->element_type); */
    1096         /* //   _assert_(numindices==numindicesup); */
    1097 
    1098         /* /\*Get parameters and inputs: *\/ */
    1099         /* Input* pressure_input=inputs->GetInput(PressureEnum); _assert_(pressure_input); */
    1100         /* Input* enthalpy_input=inputs->GetInput(EnthalpyEnum); _assert_(enthalpy_input); */
    1101         /* Input* watercolumn_input=inputs->GetInput(WatercolumnEnum); _assert_(watercolumn_input); */
    1102 
    1103         /* /\*if there is a temperate layer of zero thickness, set spc enthalpy=h_pmp at that node*\/ */
    1104         /* GaussPenta* gauss=new GaussPenta(); */
    1105         /* GaussPenta* gaussup=new GaussPenta(); */
    1106         /* for(int i=0;i<numindices;i++){ */
    1107         /*      gauss->GaussNode(this->element_type,indices[i]); */
    1108         /*      gaussup->GaussNode(this->element_type,indicesup[i]); */
    1109 
    1110         /*      /\*Check wether there is a temperate layer at the base or not *\/ */
    1111         /*      /\*check if node is temperate, if not, continue*\/ */
    1112         /*      enthalpy_input->GetInputValue(&enthalpy, gauss); */
    1113         /*      pressure_input->GetInputValue(&pressure, gauss); */
    1114         /*      watercolumn_input->GetInputValue(&watercolumn,gauss); */
    1115         /*      setspc = false; */
    1116         /*      // TODO: add case H<Hpmp && watercolumn>0.; */
    1117         /*      if (enthalpy>=element->PureIceEnthalpy(pressure)){ */
    1118         /*              /\*check if upper node is temperate, too. */
    1119         /*                      if yes, then we have a temperate layer of positive thickness and reset the spc. */
    1120         /*                      if not, apply dirichlet BC.*\/ */
    1121         /*              enthalpy_input->GetInputValue(&enthalpyup, gaussup); */
    1122         /*              pressure_input->GetInputValue(&pressureup, gaussup); */
    1123         /*              setspc=((enthalpyup<element->PureIceEnthalpy(pressureup)) && (watercolumn>=0.))?true:false; */
    1124         /*      } */
    1125 
    1126         /*      if (setspc) { */
    1127         /*              /\*Calculate enthalpy at pressure melting point *\/ */
    1128         /*              h_pmp=matpar->PureIceEnthalpy(pressure); */
    1129         /*              /\*Apply Dirichlet condition (dof = 0 here, since there is only one degree of freedom per node)*\/ */
    1130         /*              //element->ApplyConstraint(indices[i],1,h_pmp); */
    1131         /*              //nodes[indices[i]]->ApplyConstraint(1,h_pmp); */
    1132         /*      } */
    1133         /*      else { */
    1134         /*              /\*remove spc*\/ */
    1135         /*              //element->RemoveConstraint(indices[i],0); */
    1136         /*              //nodes[indices[i]]->DofInFSet(0); */
    1137         /*      } */
    1138         /* } */
    1139 
    1140         /* /\*Free ressources:*\/ */
    1141         /* xDelete<int>(indices); */
    1142         /* xDelete<int>(indicesup); */
    1143         /* delete gauss; */
    1144         /* delete gaussup; */
    1145 }/*}}}*/
    1146 
    1147 
    1148 
    1149 
     1073        /* Only update Constraints at the base of grounded ice*/
     1074        if(!(element->IsOnBed()) || element->IsFloating()) return;
     1075
     1076        /*Check wether dynamic basal boundary conditions are activated */
     1077        element->FindParam(&isdynamicbasalspc,ThermalIsdynamicbasalspcEnum);
     1078        if(!isdynamicbasalspc) return;
     1079
     1080        /*Fetch indices of basal & surface nodes for this finite element*/
     1081        Penta *penta =  (Penta *) element; // TODO: add Basal-/SurfaceNodeIndices to element.h, and change this to Element*
     1082        penta->BasalNodeIndices(&numindices,&indices,element->GetElementType());
     1083        penta->SurfaceNodeIndices(&numindicesup,&indicesup,element->GetElementType());
     1084        _assert_(numindices==numindicesup);
     1085
     1086        /*Get parameters and inputs: */
     1087        Input* pressure_input=element->GetInput(PressureEnum); _assert_(pressure_input);
     1088        Input* enthalpy_input=element->GetInput(EnthalpyEnum); _assert_(enthalpy_input);
     1089        Input* watercolumn_input=element->GetInput(WatercolumnEnum); _assert_(watercolumn_input);
     1090
     1091        /*if there is a temperate layer of zero thickness, set spc enthalpy=h_pmp at that node*/
     1092        GaussPenta* gauss=new GaussPenta();
     1093        GaussPenta* gaussup=new GaussPenta();
     1094        for(int i=0;i<numindices;i++){
     1095                gauss->GaussNode(element->GetElementType(),indices[i]);
     1096                gaussup->GaussNode(element->GetElementType(),indicesup[i]);
     1097
     1098                /*Check wether there is a temperate layer at the base or not */
     1099                /*check if node is temperate, else continue*/
     1100                enthalpy_input->GetInputValue(&enthalpy, gauss);
     1101                pressure_input->GetInputValue(&pressure, gauss);
     1102                watercolumn_input->GetInputValue(&watercolumn,gauss);
     1103                h_pmp=PureIceEnthalpy(element,pressure);
     1104                if (enthalpy>=h_pmp){
     1105                        /*check if upper node is temperate, too.
     1106                                if yes, then we have a temperate layer of positive thickness and reset the spc.
     1107                                if not, apply dirichlet BC.*/
     1108                        enthalpy_input->GetInputValue(&enthalpyup, gaussup);
     1109                        pressure_input->GetInputValue(&pressureup, gaussup);
     1110                        setspc=((enthalpyup<PureIceEnthalpy(element,pressureup)) && (watercolumn>=0.))?true:false;
     1111                }
     1112                else if(watercolumn>0.) // case H<Hpmp && watercolumn>0.
     1113                        setspc=true;
     1114                else
     1115                        setspc = false;
     1116
     1117                node=element->GetNode(indices[i]);
     1118                if (setspc)
     1119                        node->ApplyConstraint(1,h_pmp); /*apply spc*/ //nodes[indices[i]]->ApplyConstraint(1,h_pmp);
     1120                else                   
     1121                        node->DofInFSet(0); /*remove spc*/ //nodes[indices[i]]->DofInFSet(0);
     1122        }
     1123
     1124        /*Free ressources:*/
     1125        xDelete<int>(indices);
     1126        xDelete<int>(indicesup);
     1127        delete gauss;
     1128        delete gaussup;
     1129}/*}}}*/
    11501130
    11511131/*Intermediaries*/
    … …  
    11671147        int         iv;
    11681148        IssmDouble  lambda;                   /* fraction of cold ice    */
    1169         IssmDouble  kappa   ,kappa_c,kappa_t; /* enthalpy conductivities */
     1149        IssmDouble  kappa,kappa_c,kappa_t; /* enthalpy conductivities */
    11701150        IssmDouble  Hc,Ht;
    1171 
    11721151
    11731152        /*Get pressures and enthalpies on vertices*/
    … …  
    12091188                _assert_((Hc+Ht)>0.);
    12101189                lambda = Hc/(Hc+Ht);
    1211                 kappa  = 1./(lambda/kappa_c + (1.-lambda)/kappa_t);
    1212         }
     1190                kappa  = kappa_c*kappa_t/(lambda*kappa_t+(1.-lambda)*kappa_c); // ==(lambda/kappa_c + (1.-lambda)/kappa_t)^-1
     1191        }       
    12131192
    12141193        /*Clean up and return*/
    … …  
    12181197        xDelete<IssmDouble>(enthalpies);
    12191198        return kappa;
    1220 }
    1221 /*}}}*/
     1199}/*}}}*/
    12221200IssmDouble EnthalpyAnalysis::PureIceEnthalpy(Element* element,IssmDouble pressure){/*{{{*/
    12231201

  • issm/trunk-jpl/src/c/classes/Elements/Element.h

    r17014 r17027  
    168168                virtual IssmDouble GetYcoord(Gauss* gauss)=0;
    169169                virtual void   GetMaterialInputValue(IssmDouble* pvalue,Node* node,int enumtype)=0;
     170                virtual int    GetElementType(void)=0;
    170171
    171172                virtual IssmDouble SurfaceArea(void)=0;

  • issm/trunk-jpl/src/c/classes/Elements/Seg.h

    r17014 r17027  
    141141                IssmDouble  GetYcoord(Gauss* gauss){_error_("Not implemented");};
    142142                IssmDouble  GetZcoord(Gauss* gauss){_error_("not implemented yet");};
     143                int         GetElementType(void){_error_("not implemented yet");};
    143144                Gauss*      NewGauss(void){_error_("not implemented yet");};
    144145                Gauss*      NewGauss(int order);

  • issm/trunk-jpl/src/c/cores/thermal_core.cpp

    r16742 r17027  
    66#include "../toolkits/toolkits.h"
    77#include "../classes/classes.h"
     8#include "../analyses/analyses.h"
    89#include "../shared/shared.h"
    910#include "../modules/modules.h"
    … …  
    1718        int    solution_type,numoutputs;
    1819        char** requested_outputs = NULL;
     20        EnthalpyAnalysis * enthalpy_analysis = NULL;
    1921
    2022        /*first recover parameters common to all solutions*/
    … …  
    4749
    4850                /*Post process*/
    49                 if(solution_type!=SteadystateSolutionEnum) PostprocessingEnthalpyx(femmodel);
     51                enthalpy_analysis = new EnthalpyAnalysis();
     52                enthalpy_analysis->PostProcessing(femmodel);
     53                delete enthalpy_analysis;
    5054        }
    5155        else{
    … …  
    6670        /*Free ressources:*/   
    6771        if(numoutputs){for(int i=0;i<numoutputs;i++){xDelete<char>(requested_outputs[i]);} xDelete<char*>(requested_outputs);}
     72
    6873}

  • issm/trunk-jpl/src/c/shared/Elements/LliboutryDuval.cpp

    r16816 r17027  
    5353  if (enthalpy < H_sp){
    5454    Tstar = referencetemperature + enthalpy/heatcapacity - betaCC*pressure;     
    55     waterfraction = 0;
     55    waterfraction = 0.;
    5656  }
    5757  else{
    … …  
    6363
    6464  /*Get A*/
    65   if(Tstar<263.15){
    66     A=3.61*pow(10.,-13.) * exp(  -6.*pow(10.,4.)/(R*Tstar));
     65  if(Tstar<=263.15){
     66    A=3.61e-13 * exp(  -6.e+4/(R*Tstar));
    6767  }
    6868  else{
    69     A=1.73*pow(10.,  3.) * exp(-13.9*pow(10.,4.)/(R*Tstar));
     69    A=1.73e3   * exp(-13.9e+4/(R*Tstar));
    7070  }
    71   A*=(1 + 181.25*waterfraction);
     71  A*=(1. + 181.25*waterfraction);
    7272
    7373  /*Convert to B*/
    74   _assert_(n>0);
    7574  B=pow(A,-1./n);
    7675

  • issm/trunk-jpl/src/m/materials/arrhenius.m

    r16271 r17027  
    3131end
    3232
    33 wf_max=.1;
    34 if(waterfraction>wf_max)
     33wf_max=1.;
     34if(any(waterfraction>wf_max))
    3535    error(['waterfraction exceeds permitted maximum of ' num2str(wf_max) '.']);
    3636end
     37
     38%limit waterfraction to 1%
     39pos1p=find(waterfraction>0.01);
     40waterfraction(pos1p)=0.01;
    3741
    3842pos=find((temperature<TMeltingPoint(pressure)) & (waterfraction>0)); % cold, wet ice
    … …  
    7175        Qa=GetQa(T);
    7276        A0=GetA0(T);
     77        if(w>0.01)
     78           w=0.01;
     79        end
    7380        A=A0.*exp(-Qa./(R*T)).*(1+181.25*w);       
    7481    end
Note: See TracChangeset for help on using the changeset viewer.