Changeset 16031


Ignore:
Timestamp:
08/30/13 09:51:45 (12 years ago)
Author:
seroussi
Message:

CHG: cosmetics

File:
1 edited

Legend:

Unmodified
Added
Removed

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

    r16026 r16031  
    47004700        IssmDouble scalar,enthalpy,enthalpyup;
    47014701        IssmDouble pressure,pressureup;
    4702     IssmDouble watercolumn;
     4702        IssmDouble watercolumn;
    47034703        IssmDouble basis[NUMVERTICES];
    47044704        Friction*  friction=NULL;
    … …  
    47254725        Input* pressure_input=inputs->GetInput(PressureEnum);             _assert_(pressure_input);
    47264726        Input* geothermalflux_input=inputs->GetInput(BasalforcingsGeothermalfluxEnum); _assert_(geothermalflux_input);
    4727     Input* watercolumn_input=inputs->GetInput(WatercolumnEnum); _assert_(watercolumn_input);
     4727        Input* watercolumn_input=inputs->GetInput(WatercolumnEnum); _assert_(watercolumn_input);
    47284728
    47294729        /*Build friction element, needed later: */
    … …  
    47464746                        enthalpy_input->GetInputValue(&enthalpyup,gaussup);
    47474747                        pressure_input->GetInputValue(&pressureup,gaussup);
    4748             if(enthalpyup>=matpar->PureIceEnthalpy(pressureup)){
     4748                        if(enthalpyup>=matpar->PureIceEnthalpy(pressureup)){
    47494749                                // temperate ice has positive thickness: grad enthalpy*n=0.
    47504750                        }
    … …  
    47544754                }
    47554755                else{
    4756             watercolumn_input->GetInputValue(&watercolumn,gauss);
    4757             if(watercolumn==0.){
    4758                 /*add geothermal heat flux and basal friction*/
    4759                 geothermalflux_input->GetInputValue(&geothermalflux_value,gauss);
    4760                 friction->GetAlpha2(&alpha2,gauss,VxEnum,VyEnum,VzEnum);
    4761                 vx_input->GetInputValue(&vx,gauss);
    4762                 vy_input->GetInputValue(&vy,gauss);
    4763                 basalfriction=alpha2*(pow(vx,2.0)+pow(vy,2.0));
    4764 
    4765                 scalar=gauss->weight*Jdet2d*(basalfriction+geothermalflux_value)/(rho_ice);
    4766                 if(reCast<bool,IssmDouble>(dt)) scalar=dt*scalar;
    4767 
    4768                 for(i=0;i<numdof;i++) pe->values[i]+=scalar*basis[i];
    4769             }
    4770             else { /*do nothing (water layer acts as insulation)*/ }
     4756                        watercolumn_input->GetInputValue(&watercolumn,gauss);
     4757                        if(watercolumn==0.){
     4758                                /*add geothermal heat flux and basal friction*/
     4759                                geothermalflux_input->GetInputValue(&geothermalflux_value,gauss);
     4760                                friction->GetAlpha2(&alpha2,gauss,VxEnum,VyEnum,VzEnum);
     4761                                vx_input->GetInputValue(&vx,gauss);
     4762                                vy_input->GetInputValue(&vy,gauss);
     4763                                basalfriction=alpha2*(pow(vx,2.0)+pow(vy,2.0));
     4764
     4765                                scalar=gauss->weight*Jdet2d*(basalfriction+geothermalflux_value)/(rho_ice);
     4766                                if(reCast<bool,IssmDouble>(dt)) scalar=dt*scalar;
     4767
     4768                                for(i=0;i<numdof;i++) pe->values[i]+=scalar*basis[i];
     4769                        }
     4770                        else  /*do nothing (water layer acts as insulation)*/
    47714771                }
    47724772        }
    … …  
    51485148        const int    numdof=NDOF1*NUMVERTICES;
    51495149
    5150         bool   converged=false;
    5151         int    i,rheology_law;
     5150        bool       converged=false;
     5151        int        i,rheology_law;
    51525152        IssmDouble xyz_list[NUMVERTICES][3];
    51535153        IssmDouble values[numdof];
    … …  
    51575157        IssmDouble B[numdof];
    51585158        IssmDouble B_average,s_average;
    5159         int*   doflist=NULL;
     5159        int*       doflist=NULL;
    51605160
    51615161        /*Get dof list: */
    … …  
    52085208                                this->material->inputs->AddInput(new PentaInput(MaterialsRheologyBEnum,B,P1Enum));
    52095209                                break;
    5210             case LliboutryDuvalEnum:
    5211                 B_average=LliboutryDuval((values[0]+values[1]+values[2]+values[3]+values[4]+values[5])/6.0,                                                   (pressure[0]+pressure[1]+pressure[2]+pressure[3]+pressure[4]+pressure[5])/6.0,                                       material->GetN());
    5212                 for(i=0;i<numdof;i++) B[i]=B_average;
     5210                        case LliboutryDuvalEnum:
     5211                                B_average=LliboutryDuval((values[0]+values[1]+values[2]+values[3]+values[4]+values[5])/6.0,
     5212                                                        (pressure[0]+pressure[1]+pressure[2]+pressure[3]+pressure[4]+pressure[5])/6.0,
     5213                                                        material->GetN());
     5214                                for(i=0;i<numdof;i++) B[i]=B_average;
    52135215                                this->material->inputs->AddInput(new PentaInput(MaterialsRheologyBEnum,B,P1Enum));
    52145216                                break;
    52155217                        default:
    52165218                                _error_("Rheology law " << EnumToStringx(rheology_law) << " not supported yet");
    5217 
    52185219                }
    52195220        }
    … …  
    53055306        IssmDouble  xyz_list[NUMVERTICES][3];
    53065307        IssmDouble  xyz_list_tria[NUMVERTICES2D][3];
    5307         IssmDouble  heatflux, geothermalflux_value;
     5308        IssmDouble  heatflux, geothermalflux_value,kappa;
    53085309        IssmDouble  vec_heatflux[3];
    53095310        IssmDouble  normal_base[3], d1enthalpy[3];
    5310         IssmDouble  kappa;
    5311     IssmDouble  basalmeltingrate[NUMVERTICES], watercolumn[NUMVERTICES], meltingrate_enthalpy;
     5311        IssmDouble  basalmeltingrate[NUMVERTICES], watercolumn[NUMVERTICES], meltingrate_enthalpy;
    53125312        IssmDouble  enthalpy[NUMVERTICES], enthalpyup;
    53135313        IssmDouble  pressure, pressureup;
    … …  
    53155315        IssmDouble  latentheat, rho_ice;
    53165316        IssmDouble  basalfriction, alpha2;
    5317         IssmDouble  vx,vy,vz;
    5318     IssmDouble  connectivity;
    5319         IssmDouble  dt;
     5317        IssmDouble  vx,vy,vz,dt,connectivity;
    53205318        Friction   *friction  = NULL;
    53215319
    … …  
    53295327        /*Fetch parameters and inputs */
    53305328        latentheat=matpar->GetLatentHeat();
    5331     rho_ice=this->matpar->GetRhoIce();
    5332         Input* watercolumn_input=inputs->GetInput(WatercolumnEnum);       _assert_(watercolumn_input);
    5333     Input* basalmeltingrate_input=inputs->GetInput(BasalforcingsMeltingRateEnum);   _assert_(basalmeltingrate_input);
    5334         Input* enthalpy_input=inputs->GetInput(EnthalpyEnum);             _assert_(enthalpy_input);
    5335         Input* pressure_input=inputs->GetInput(PressureEnum);             _assert_(pressure_input);
    5336         Input* vx_input=inputs->GetInput(VxEnum);                         _assert_(vx_input);
    5337         Input* vy_input=inputs->GetInput(VyEnum);                         _assert_(vy_input);
    5338         Input* vz_input=inputs->GetInput(VzEnum);                         _assert_(vz_input);
     5329        rho_ice=this->matpar->GetRhoIce();
     5330        Input* watercolumn_input=inputs->GetInput(WatercolumnEnum);                    _assert_(watercolumn_input);
     5331        Input* basalmeltingrate_input=inputs->GetInput(BasalforcingsMeltingRateEnum);  _assert_(basalmeltingrate_input);
     5332        Input* enthalpy_input=inputs->GetInput(EnthalpyEnum);                          _assert_(enthalpy_input);
     5333        Input* pressure_input=inputs->GetInput(PressureEnum);                          _assert_(pressure_input);
     5334        Input* vx_input=inputs->GetInput(VxEnum);                                      _assert_(vx_input);
     5335        Input* vy_input=inputs->GetInput(VyEnum);                                      _assert_(vy_input);
     5336        Input* vz_input=inputs->GetInput(VzEnum);                                      _assert_(vz_input);
    53395337        Input* geothermalflux_input=inputs->GetInput(BasalforcingsGeothermalfluxEnum); _assert_(geothermalflux_input);
    53405338
    … …  
    53515349                gauss->GaussVertex(iv);
    53525350                gaussup->GaussVertex(iv+3);
    5353                
    5354         checkpositivethickness=true;
     5351
     5352                checkpositivethickness=true;
    53555353                watercolumn_input->GetInputValue(&watercolumn[iv], gauss);
    5356         basalmeltingrate_input->GetInputValue(&basalmeltingrate[iv],gauss);
     5354                basalmeltingrate_input->GetInputValue(&basalmeltingrate[iv],gauss);
    53575355                enthalpy_input->GetInputValue(&enthalpy[iv], gauss);
    53585356                pressure_input->GetInputValue(&pressure, gauss);
    53595357
    53605358                /*Calculate basal meltingrate after Fig.5 of A.Aschwanden 2012*/
    5361         meltingrate_enthalpy=0.;
     5359                meltingrate_enthalpy=0.;
    53625360                if((watercolumn[iv]>0.) && (enthalpy[iv]<matpar->PureIceEnthalpy(pressure))){
    5363             /*ensure that no ice is at T<Tm(p), if water layer present*/
    5364             enthalpy[iv]=matpar->PureIceEnthalpy(pressure);
    5365             //meltingrate_enthalpy=0.;
     5361                        /*ensure that no ice is at T<Tm(p), if water layer present*/
     5362                        enthalpy[iv]=matpar->PureIceEnthalpy(pressure);
     5363                        //meltingrate_enthalpy=0.;
    53665364                }
    53675365                else if(enthalpy[iv]<matpar->PureIceEnthalpy(pressure)){
    5368             /*cold base: set q*n=q_geo*n+frictionheating as Neumann BC in Penta::CreatePVectorEnthalpySheet*/
     5366                        /*cold base: set q*n=q_geo*n+frictionheating as Neumann BC in Penta::CreatePVectorEnthalpySheet*/
    53695367                        meltingrate_enthalpy=0.;
    5370             checkpositivethickness=false;
    5371                 }
    5372         else {/*do nothing, go to next check*/}
    5373 
    5374         if(checkpositivethickness){
    5375             /*ok, from here on all basal ice is temperate. Check for positive thickness of layer of temperate ice. */
    5376             istemperatelayer=false;
    5377             enthalpy_input->GetInputValue(&enthalpyup, gaussup);
    5378             pressure_input->GetInputValue(&pressureup, gaussup);   
    5379             if(enthalpyup>=matpar->PureIceEnthalpy(pressureup))
    5380                 istemperatelayer=true;
    5381             if(istemperatelayer)
    5382                 for(i=0;i<3;i++) vec_heatflux[i]=0.;
    5383             else{
    5384                 enthalpy_input->GetInputDerivativeValue(&d1enthalpy[0],&xyz_list[0][0],gauss);
    5385                 kappa=matpar->GetEnthalpyDiffusionParameter(enthalpy[iv],pressure);
    5386                 for(i=0;i<3;i++) vec_heatflux[i]=-kappa*d1enthalpy[i];
    5387             }
    5388 
    5389             /*geothermal heatflux*/
    5390             BedNormal(&normal_base[0],xyz_list_tria);
    5391             heatflux=0.;
    5392             for(i=0;i<3;i++) heatflux+=(vec_heatflux[i])*normal_base[i];
    5393             geothermalflux_input->GetInputValue(&geothermalflux_value,gauss);
    5394 
    5395             /*basal friction*/
    5396             friction->GetAlpha2(&alpha2,gauss,VxEnum,VyEnum,VzEnum);
    5397             vx_input->GetInputValue(&vx,gauss);
    5398             vy_input->GetInputValue(&vy,gauss);
    5399             vz_input->GetInputValue(&vz,gauss);
    5400             basalfriction=alpha2*(pow(vx,2.0)+pow(vy,2.0)+pow(vz,2.0));
    5401 
    5402             matpar->EnthalpyToThermal(&temperature, &waterfraction, enthalpy[iv],pressure);
    5403             meltingrate_enthalpy=(basalfriction-(heatflux-geothermalflux_value))/((1-waterfraction)*latentheat*rho_ice); // m/yr water equivalent
    5404         }
     5368                        checkpositivethickness=false;
     5369                }
     5370                else {/*do nothing, go to next check*/}
     5371
     5372                if(checkpositivethickness){
     5373                        /*ok, from here on all basal ice is temperate. Check for positive thickness of layer of temperate ice. */
     5374                        istemperatelayer=false;
     5375                        enthalpy_input->GetInputValue(&enthalpyup, gaussup);
     5376                        pressure_input->GetInputValue(&pressureup, gaussup);   
     5377                        if(enthalpyup>=matpar->PureIceEnthalpy(pressureup)) istemperatelayer=true;
     5378                        if(istemperatelayer) for(i=0;i<3;i++) vec_heatflux[i]=0.;
     5379                        else{
     5380                                enthalpy_input->GetInputDerivativeValue(&d1enthalpy[0],&xyz_list[0][0],gauss);
     5381                                kappa=matpar->GetEnthalpyDiffusionParameter(enthalpy[iv],pressure);
     5382                                for(i=0;i<3;i++) vec_heatflux[i]=-kappa*d1enthalpy[i];
     5383                        }
     5384
     5385                        /*geothermal heatflux*/
     5386                        BedNormal(&normal_base[0],xyz_list_tria);
     5387                        heatflux=0.;
     5388                        for(i=0;i<3;i++) heatflux+=(vec_heatflux[i])*normal_base[i];
     5389                        geothermalflux_input->GetInputValue(&geothermalflux_value,gauss);
     5390
     5391                        /*basal friction*/
     5392                        friction->GetAlpha2(&alpha2,gauss,VxEnum,VyEnum,VzEnum);
     5393                        vx_input->GetInputValue(&vx,gauss);
     5394                        vy_input->GetInputValue(&vy,gauss);
     5395                        vz_input->GetInputValue(&vz,gauss);
     5396                        basalfriction=alpha2*(pow(vx,2.0)+pow(vy,2.0)+pow(vz,2.0));
     5397
     5398                        matpar->EnthalpyToThermal(&temperature, &waterfraction, enthalpy[iv],pressure);
     5399                        meltingrate_enthalpy=(basalfriction-(heatflux-geothermalflux_value))/((1-waterfraction)*latentheat*rho_ice); // m/yr water equivalent
     5400                }
    54055401
    54065402                /*Update water column, basal meltingrate*/
    5407         connectivity=IssmDouble(vertices[iv]->Connectivity());
    5408         meltingrate_enthalpy/=connectivity; // divide meltingrate_enthalpy by connectivity to address multiple iterations over vertex
    5409         basalmeltingrate[iv]+=meltingrate_enthalpy;
     5403                connectivity=IssmDouble(vertices[iv]->Connectivity());
     5404                meltingrate_enthalpy/=connectivity; // divide meltingrate_enthalpy by connectivity to address multiple iterations over vertex
     5405                basalmeltingrate[iv]+=meltingrate_enthalpy;
    54105406                this->parameters->FindParam(&dt,TimesteppingTimeStepEnum);
    54115407                if(reCast<bool,IssmDouble>(dt))
    5412             watercolumn[iv]+=dt*meltingrate_enthalpy;
     5408                watercolumn[iv]+=dt*meltingrate_enthalpy;
    54135409                else
    5414             watercolumn[iv]+=meltingrate_enthalpy;
     5410                watercolumn[iv]+=meltingrate_enthalpy;
    54155411        } 
    5416     /*feed updated variables back into model*/
    5417     this->inputs->AddInput(new PentaInput(EnthalpyEnum, enthalpy, P1Enum));
    5418     this->inputs->AddInput(new PentaInput(WatercolumnEnum, watercolumn, P1Enum));
    5419     this->inputs->AddInput(new PentaInput(BasalforcingsMeltingRateEnum, basalmeltingrate, P1Enum));
     5412        /*feed updated variables back into model*/
     5413        this->inputs->AddInput(new PentaInput(EnthalpyEnum,enthalpy,P1Enum));
     5414        this->inputs->AddInput(new PentaInput(WatercolumnEnum,watercolumn,P1Enum));
     5415        this->inputs->AddInput(new PentaInput(BasalforcingsMeltingRateEnum,basalmeltingrate,P1Enum));
    54205416
    54215417        /*Clean up and return*/
    54225418        delete gauss;
    54235419        delete gaussup;
    5424     delete friction;
     5420        delete friction;
    54255421}
    54265422/*}}}*/
    … …  
    54295425   
    54305426    /*Intermediaries*/
    5431     bool isenthalpy;
    5432     IssmDouble waterfraction[NUMVERTICES], temperature[NUMVERTICES];
    5433     IssmDouble enthalpy[NUMVERTICES], pressure[NUMVERTICES];
    5434     IssmDouble latentheat, dt;
    5435     IssmDouble connectivity;
    5436     GaussPenta* gauss;
    5437    
    5438     Input* watercolumn_input=inputs->GetInput(WatercolumnEnum);       _assert_(watercolumn_input);
     5427        bool isenthalpy;
     5428        IssmDouble waterfraction[NUMVERTICES], temperature[NUMVERTICES];
     5429        IssmDouble enthalpy[NUMVERTICES], pressure[NUMVERTICES];
     5430        IssmDouble latentheat, dt;
     5431        IssmDouble connectivity;
     5432        GaussPenta* gauss;
     5433
     5434        Input* watercolumn_input=inputs->GetInput(WatercolumnEnum);       _assert_(watercolumn_input);
    54395435        Input* enthalpy_input=inputs->GetInput(EnthalpyEnum);             _assert_(enthalpy_input);
    5440     Input* pressure_input=inputs->GetInput(PressureEnum);             _assert_(pressure_input);
    5441    
    5442     /*Check wether enthalpy is activated*/
     5436        Input* pressure_input=inputs->GetInput(PressureEnum);             _assert_(pressure_input);
     5437
     5438        /*Check wether enthalpy is activated*/
    54435439        parameters->FindParam(&isenthalpy,ThermalIsenthalpyEnum);
    54445440        if(!isenthalpy) return;       
    5445    
    5446     this->parameters->FindParam(&dt,TimesteppingTimeStepEnum);
    5447     latentheat=matpar->GetLatentHeat();
    5448 
    5449     gauss=new GaussPenta();
    5450     for(int iv=0;iv<NUMVERTICES;iv++){
    5451         gauss->GaussVertex(iv);
    5452         enthalpy_input->GetInputValue(&enthalpy[iv], gauss);
    5453         pressure_input->GetInputValue(&pressure[iv], gauss);
    5454         matpar->EnthalpyToThermal(&temperature[iv],&waterfraction[iv], enthalpy[iv],pressure[iv]);
    5455    
    5456         /*drain water fraction & update enthalpy*/
    5457         // TODO: replace connectivity-wise draining by draining once per node per timestep
    5458         connectivity=IssmDouble(vertices[iv]->Connectivity());
    5459         waterfraction[iv]-=DrainageFunctionWaterfraction(waterfraction[iv], dt)/connectivity;
    5460         matpar->ThermalToEnthalpy(&enthalpy[iv], temperature[iv], waterfraction[iv], pressure[iv]);       
    5461     }
    5462     /*feed updated results back into model*/
    5463     this->inputs->AddInput(new PentaInput(EnthalpyEnum,enthalpy,P1Enum));
    5464     this->inputs->AddInput(new PentaInput(WaterfractionEnum,waterfraction,P1Enum));
    5465     // this->inputs->AddInput(new PentaInput(TemperatureEnum,temperature,P1Enum));    // temperature should not change during drainage...
     5441
     5442        this->parameters->FindParam(&dt,TimesteppingTimeStepEnum);
     5443        latentheat=matpar->GetLatentHeat();
     5444
     5445        gauss=new GaussPenta();
     5446        for(int iv=0;iv<NUMVERTICES;iv++){
     5447                gauss->GaussVertex(iv);
     5448                enthalpy_input->GetInputValue(&enthalpy[iv], gauss);
     5449                pressure_input->GetInputValue(&pressure[iv], gauss);
     5450                matpar->EnthalpyToThermal(&temperature[iv],&waterfraction[iv], enthalpy[iv],pressure[iv]);
     5451
     5452                /*drain water fraction & update enthalpy*/
     5453                // TODO: replace connectivity-wise draining by draining once per node per timestep
     5454                connectivity=IssmDouble(vertices[iv]->Connectivity());
     5455                waterfraction[iv]-=DrainageFunctionWaterfraction(waterfraction[iv], dt)/connectivity;
     5456                matpar->ThermalToEnthalpy(&enthalpy[iv], temperature[iv], waterfraction[iv], pressure[iv]);       
     5457        }
     5458        /*feed updated results back into model*/
     5459        this->inputs->AddInput(new PentaInput(EnthalpyEnum,enthalpy,P1Enum));
     5460        this->inputs->AddInput(new PentaInput(WaterfractionEnum,waterfraction,P1Enum));
     5461        // this->inputs->AddInput(new PentaInput(TemperatureEnum,temperature,P1Enum));    // temperature should not change during drainage...
    54665462}
    54675463/*}}}*/
Note: See TracChangeset for help on using the changeset viewer.