Changeset 8496


Ignore:
Timestamp:
06/03/11 10:27:30 (14 years ago)
Author:
seroussi
Message:

src/c enthalpy

Location:
issm/trunk/src/c
Files:
7 edited

Legend:

Unmodified
Added
Removed

  • issm/trunk/src/c/EnumDefinitions/EnumDefinitions.h

    r8487 r8496  
    451451        EnthalpyAnalysisEnum,
    452452        EnthalpyEnum,
     453        EnthalpyPicardEnum,
    453454        WaterFractionEnum,
    454455        ReferenceTemperatureEnum

  • issm/trunk/src/c/modules/EnumToStringx/EnumToStringx.cpp

    r8487 r8496  
    394394                case EnthalpyAnalysisEnum : return "EnthalpyAnalysis";
    395395                case EnthalpyEnum : return "Enthalpy";
     396                case EnthalpyPicardEnum : return "EnthalpyPicard";
    396397                case WaterFractionEnum : return "WaterFraction";
    397398                case ReferenceTemperatureEnum : return "ReferenceTemperature";

  • issm/trunk/src/c/modules/StringToEnumx/StringToEnumx.cpp

    r8487 r8496  
    392392        else if (strcmp(name,"EnthalpyAnalysis")==0) return EnthalpyAnalysisEnum;
    393393        else if (strcmp(name,"Enthalpy")==0) return EnthalpyEnum;
     394        else if (strcmp(name,"EnthalpyPicard")==0) return EnthalpyPicardEnum;
    394395        else if (strcmp(name,"WaterFraction")==0) return WaterFractionEnum;
    395396        else if (strcmp(name,"ReferenceTemperature")==0) return ReferenceTemperatureEnum;

  • issm/trunk/src/c/objects/Elements/Penta.cpp

    r8490 r8496  
    612612                        Ke=CreateKMatrixThermal();
    613613                        break;
     614                case EnthalpyAnalysisEnum:
     615                        Ke=CreateKMatrixEnthalpy();
     616                        break;
    614617                case MeltingAnalysisEnum:
    615618                        Ke=CreateKMatrixMelting();
    … …  
    18341837        double     gravity,rho_ice,rho_water;
    18351838        double     heatcapacity,thermalconductivity,dt;
    1836         double     latentheat,kappa=1;
     1839        double     pressure,enthalpy;
     1840        double     latentheat,kappa;
    18371841        double     tau_parameter,diameter;
    18381842        double     xyz_list[NUMVERTICES][3];
    … …  
    18701874        this->parameters->FindParam(&artdiff,ArtDiffEnum);
    18711875        this->parameters->FindParam(&epsvel,EpsVelEnum);
    1872         Input* vx_input=inputs->GetInput(VxEnum);      _assert_(vx_input);
    1873         Input* vy_input=inputs->GetInput(VyEnum);      _assert_(vy_input);
    1874         Input* vz_input=inputs->GetInput(VzEnum);      _assert_(vz_input);
    1875         Input* vxm_input=inputs->GetInput(VxMeshEnum); _assert_(vxm_input);
    1876         Input* vym_input=inputs->GetInput(VyMeshEnum); _assert_(vym_input);
    1877         Input* vzm_input=inputs->GetInput(VzMeshEnum); _assert_(vzm_input);
     1876        Input* pressure_input=inputs->GetInput(PressureEnum);      _assert_(pressure_input);
     1877        Input* enthalpy_input=inputs->GetInput(EnthalpyEnum);      _assert_(enthalpy_input);
     1878        Input* vx_input=inputs->GetInput(VxEnum);                  _assert_(vx_input);
     1879        Input* vy_input=inputs->GetInput(VyEnum);                  _assert_(vy_input);
     1880        Input* vz_input=inputs->GetInput(VzEnum);                  _assert_(vz_input);
     1881        Input* vxm_input=inputs->GetInput(VxMeshEnum);             _assert_(vxm_input);
     1882        Input* vym_input=inputs->GetInput(VyMeshEnum);             _assert_(vym_input);
     1883        Input* vzm_input=inputs->GetInput(VzMeshEnum);             _assert_(vzm_input);
    18781884        if (artdiff==2) diameter=MinEdgeLength(xyz_list);
    18791885
    … …  
    18911897                GetBConduct(&B_conduct[0][0],&xyz_list[0][0],gauss);
    18921898
    1893                 //kappa=GetEnthalpyDiffusivityParameter(rho_ice,heatcapacity,thermalconductivity,latentheat,enthalpy);
     1899                enthalpy_input->GetParameterValue(&enthalpy, gauss);
     1900                pressure_input->GetParameterValue(&pressure, gauss);
     1901                kappa=matpar->GetEnthalpyDiffusionParameter(enthalpy,pressure);
    18941902                D_scalar_conduct=gauss->weight*Jdet*kappa;
    18951903                if(dt) D_scalar_conduct=D_scalar_conduct*dt;
    … …  
    23772385                        pe=CreatePVectorThermal();
    23782386                        break;
     2387                case EnthalpyAnalysisEnum:
     2388                        pe=CreatePVectorEnthalpy();
     2389                        break;
    23792390                case MeltingAnalysisEnum:
    23802391                        pe=CreatePVectorMelting();
    … …  
    34123423        int        i,j,ig;
    34133424        double     Jdet2d;
     3425        double     heatcapacity,h_pmp;
    34143426        double     mixed_layer_capacity,thermal_exchange_velocity;
    34153427        double     rho_ice,rho_water,pressure,dt,scalar_ocean;
    3416         double     meltingpoint,beta,heatcapacity,h_pmp;
    34173428        double     xyz_list[NUMVERTICES][3];
    34183429        double     xyz_list_tria[NUMVERTICES2D][3];
    … …  
    34343445        rho_ice=matpar->GetRhoIce();
    34353446        heatcapacity=matpar->GetHeatCapacity();
    3436         beta=matpar->GetBeta();
    3437         meltingpoint=matpar->GetMeltingPoint();
    34383447        this->parameters->FindParam(&dt,DtEnum);
    34393448        Input* pressure_input=inputs->GetInput(PressureEnum); _assert_(pressure_input);
    … …  
    34493458
    34503459                pressure_input->GetParameterValue(&pressure,gauss);
    3451                 h_pmp=meltingpoint-beta*pressure;
     3460                h_pmp=matpar->PureIceEnthalpy(pressure);
    34523461
    34533462                scalar_ocean=gauss->weight*Jdet2d*rho_water*mixed_layer_capacity*thermal_exchange_velocity*(h_pmp)/(rho_ice*heatcapacity);
    … …  
    41774186                GetSolutionFromInputsThermal(solution);
    41784187        }
     4188        else if(analysis_type==EnthalpyAnalysisEnum){
     4189                GetSolutionFromInputsEnthalpy(solution);
     4190        }
    41794191        else{
    41804192                _error_("analysis: %i (%s) not supported yet",analysis_type,EnumToStringx(analysis_type));
    … …  
    43594371                t_input->GetParameterValue(&temp,gauss);
    43604372                values[i]=temp;
     4373        }
     4374
     4375        /*Add value to global vector*/
     4376        VecSetValues(solution,numdof,doflist,(const double*)values,INSERT_VALUES);
     4377
     4378        /*Free ressources:*/
     4379        delete gauss;
     4380        xfree((void**)&doflist);
     4381}
     4382/*}}}*/
     4383/*FUNCTION Penta::GetSolutionFromInputsEnthalpy{{{1*/
     4384void  Penta::GetSolutionFromInputsEnthalpy(Vec solution){
     4385
     4386        const int    numdof=NDOF1*NUMVERTICES;
     4387
     4388        int          i;
     4389        int*         doflist=NULL;
     4390        double       values[numdof];
     4391        double       enthalpy;
     4392        GaussPenta   *gauss=NULL;
     4393
     4394        /*Get dof list: */
     4395        GetDofList(&doflist,NoneApproximationEnum,GsetEnum);
     4396        Input* h_input=inputs->GetInput(EnthalpyEnum); _assert_(h_input);
     4397
     4398        gauss=new GaussPenta();
     4399        for(i=0;i<NUMVERTICES;i++){
     4400                /*Recover temperature*/
     4401                gauss->GaussVertex(i);
     4402                h_input->GetParameterValue(&enthalpy,gauss);
     4403                values[i]=enthalpy;
    43614404        }
    43624405
    … …  
    50425085                InputUpdateFromSolutionThermal( solution);
    50435086        }
     5087        else if (analysis_type==EnthalpyAnalysisEnum){
     5088                InputUpdateFromSolutionEnthalpy( solution);
     5089        }
    50445090        else if (analysis_type==MeltingAnalysisEnum){
    50455091                InputUpdateFromSolutionOneDof(solution,BasalMeltingRateEnum);
    … …  
    60266072        double xyz_list[NUMVERTICES][3];
    60276073        double values[numdof];
     6074        double pressure[NUMVERTICES];
    60286075        double temperatures[numdof];
    60296076        double waterfraction[numdof];
    … …  
    60456092        /*Get all inputs and parameters*/
    60466093        GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES);
     6094        GetParameterListOnVertices(&pressure[0],PressureEnum);
    60476095        Input* surface_input=inputs->GetInput(SurfaceEnum); _assert_(surface_input);
     6096       
    60486097
    60496098        this->inputs->GetParameterValue(&converged,ConvergedEnum);
    60506099        if(converged){
    60516100                /*Convert enthalpy into temperature and water fraction*/
    6052 //              for(i=0;i<numdof;i++){
    6053 //                      if values[i]<GetPureIceEnthalpy(){
    6054 //                              temperatures[i]=values[i];
    6055 //                              waterfraction[i]=0;
    6056 //                      }
    6057 //                      else{
    6058 //                              temperatures[i]=values[i];
    6059 //                              waterfraction[i]=values[i];
    6060 //                      }
    6061 //              }
    6062 
     6101                for(i=0;i<numdof;i++) matpar->EnthalpyToThermal(&temperatures[i],&waterfraction[i],values[i],pressure[i]);
     6102                       
    60636103                this->inputs->AddInput(new PentaVertexInput(EnthalpyEnum,values));
    60646104                this->inputs->AddInput(new PentaVertexInput(WaterFractionEnum,waterfraction));
    … …  
    60736113                                break;
    60746114                        case PatersonEnum:
    6075                                 B_average=Paterson((values[0]+values[1]+values[2]+values[3]+values[4]+values[5])/6.0);
     6115                                B_average=Paterson((temperatures[0]+temperatures[1]+temperatures[2]+temperatures[3]+temperatures[4]+temperatures[5])/6.0);
    60766116                                for(i=0;i<numdof;i++) B[i]=B_average;
    60776117                                this->matice->inputs->AddInput(new PentaVertexInput(RheologyBEnum,B));
    … …  
    60796119                        case ArrheniusEnum:
    60806120                                surface_input->GetParameterAverage(&s_average);
    6081                                 B_average=Arrhenius((values[0]+values[1]+values[2]+values[3]+values[4]+values[5])/6.0,
     6121                                B_average=Arrhenius((temperatures[0]+temperatures[1]+temperatures[2]+temperatures[3]+temperatures[4]+temperatures[5])/6.0,
    60826122                                                        s_average-((xyz_list[0][2]+xyz_list[1][2]+xyz_list[2][2]+xyz_list[3][2]+xyz_list[4][2]+xyz_list[5][2])/6.0),
    60836123                                                        matice->GetN());
    … …  
    60906130                }
    60916131        }
    6092 //      else{
    6093 //              this->inputs->AddInput(new PentaVertexInput(EnthalpyPicardEnum,values));
    6094 //      }
     6132        else{
     6133                this->inputs->AddInput(new PentaVertexInput(EnthalpyPicardEnum,values));
     6134        }
    60956135
    60966136        /*Free ressources:*/

  • issm/trunk/src/c/objects/Elements/Penta.h

    r8487 r8496  
    225225                void      GetSolutionFromInputsDiagnosticVert(Vec solutiong);
    226226                void      GetSolutionFromInputsThermal(Vec solutiong);
     227                void      GetSolutionFromInputsEnthalpy(Vec solutiong);
    227228                double GetStabilizationParameter(double u, double v, double w, double diameter, double rho_ice, double heatcapacity, double thermalconductivity);
    228229                void    GetStrainRate3dPattyn(double* epsilon,double* xyz_list, GaussPenta* gauss, Input* vx_input, Input* vy_input);

  • issm/trunk/src/c/objects/Materials/Matpar.cpp

    r8490 r8496  
    3232        double  matpar_beta;
    3333        double  matpar_meltingpoint;
     34        double  matpar_referencetemperature;
    3435        double  matpar_mixed_layer_capacity;
    3536        double  matpar_thermal_exchange_velocity;
    … …  
    4647        matpar_beta=iomodel->beta;
    4748        matpar_meltingpoint=iomodel->meltingpoint;
     49        matpar_referencetemperature=iomodel->referencetemperature;
    4850        matpar_mixed_layer_capacity=iomodel->mixed_layer_capacity;
    4951        matpar_thermal_exchange_velocity=iomodel->thermal_exchange_velocity;
    … …  
    5961        this->beta=matpar_beta;
    6062        this->meltingpoint=matpar_meltingpoint;
     63        this->referencetemperature=matpar_referencetemperature;
    6164        this->mixed_layer_capacity=matpar_mixed_layer_capacity;
    6265        this->thermal_exchange_velocity=matpar_thermal_exchange_velocity;
    … …  
    8689        printf("   beta: %g\n",beta);
    8790        printf("   meltingpoint: %g\n",meltingpoint);
     91        printf("   referencetemperature: %g\n",referencetemperature);
    8892        printf("   mixed_layer_capacity: %g\n",mixed_layer_capacity);
    8993        printf("   thermal_exchange_velocity: %g\n",thermal_exchange_velocity);
    … …  
    106110        printf("   beta: %g\n",beta);
    107111        printf("   meltingpoint: %g\n",meltingpoint);
     112        printf("   referencetemperature: %g\n",referencetemperature);
    108113        printf("   mixed_layer_capacity: %g\n",mixed_layer_capacity);
    109114        printf("   thermal_exchange_velocity: %g\n",thermal_exchange_velocity);
    … …  
    147152        memcpy(marshalled_dataset,&beta,sizeof(beta));marshalled_dataset+=sizeof(beta);
    148153        memcpy(marshalled_dataset,&meltingpoint,sizeof(meltingpoint));marshalled_dataset+=sizeof(meltingpoint);
     154        memcpy(marshalled_dataset,&referencetemperature,sizeof(referencetemperature));marshalled_dataset+=sizeof(referencetemperature);
    149155        memcpy(marshalled_dataset,&mixed_layer_capacity,sizeof(mixed_layer_capacity));marshalled_dataset+=sizeof(mixed_layer_capacity);
    150156        memcpy(marshalled_dataset,&thermal_exchange_velocity,sizeof(thermal_exchange_velocity));marshalled_dataset+=sizeof(thermal_exchange_velocity);
    … …  
    168174                sizeof(beta)+
    169175                sizeof(meltingpoint)+
     176                sizeof(referencetemperature)+
    170177                sizeof(mixed_layer_capacity)+
    171178                sizeof(thermal_exchange_velocity)+
    … …  
    195202        memcpy(&beta,marshalled_dataset,sizeof(beta));marshalled_dataset+=sizeof(beta);
    196203        memcpy(&meltingpoint,marshalled_dataset,sizeof(meltingpoint));marshalled_dataset+=sizeof(meltingpoint);
     204        memcpy(&referencetemperature,marshalled_dataset,sizeof(referencetemperature));marshalled_dataset+=sizeof(referencetemperature);
    197205        memcpy(&mixed_layer_capacity,marshalled_dataset,sizeof(mixed_layer_capacity));marshalled_dataset+=sizeof(mixed_layer_capacity);
    198206        memcpy(&thermal_exchange_velocity,marshalled_dataset,sizeof(thermal_exchange_velocity));marshalled_dataset+=sizeof(thermal_exchange_velocity);
    … …  
    282290                        break;
    283291
     292                case  ReferenceTemperatureEnum:
     293                        this->referencetemperature=constant;
     294                        break;
     295
    284296                case  MixedLayerCapacityEnum:
    285297                        this->mixed_layer_capacity=constant;
    … …  
    347359double Matpar::GetMeltingPoint(){
    348360        return meltingpoint;
     361}
     362/*}}}1*/
     363/*FUNCTION Matpar::GetReferenceTemperature {{{1*/
     364double Matpar::GetReferenceTemperature(){
     365        return referencetemperature;
    349366}
    350367/*}}}1*/
    … …  
    390407}
    391408/*}}}1*/
     409/*FUNCTION Matpar::PureIceEnthalpy{{{1*/
     410double Matpar::PureIceEnthalpy(double pressure){
     411        return heatcapacity*(TMeltingPoint(pressure)-referencetemperature);
     412}
     413/*}}}1*/
     414/*FUNCTION Matpar::GetEnthalpyDiffusionParameter{{{1*/
     415double Matpar::GetEnthalpyDiffusionParameter(double enthalpy,double pressure){
     416        return thermalconductivity/(rho_ice*heatcapacity);
     417}
     418/*}}}1*/
     419/*FUNCTION Matpar::EnthalpyToThermal {{{1*/
     420void Matpar::EnthalpyToThermal(double* ptemperature,double* pwaterfraction,double enthalpy,double pressure){
     421
     422        /*Ouput*/
     423        double temperature,waterfraction;
     424       
     425        if(enthalpy<PureIceEnthalpy(pressure)){
     426                temperature=referencetemperature+enthalpy/heatcapacity;
     427                waterfraction=0;
     428        }
     429        else{
     430                temperature=TMeltingPoint(pressure);
     431                waterfraction=(enthalpy-PureIceEnthalpy(pressure))/latentheat;
     432        }
     433
     434        /*Assign output pointers:*/
     435        *pwaterfraction=waterfraction;
     436        *ptemperature=temperature;
     437}
     438/*}}}1*/
     439/*FUNCTION Matpar::ThermalToEnthalpy {{{1*/
     440void Matpar::ThermalToEnthalpy(double * penthalpy,double temperature,double waterfraction,double pressure){
     441
     442        /*Ouput*/
     443        double enthalpy;
     444       
     445        if(temperature<TMeltingPoint(pressure)){
     446                enthalpy=heatcapacity*(temperature-referencetemperature);
     447        }
     448        else{
     449                enthalpy=PureIceEnthalpy(pressure)+latentheat*waterfraction;
     450        }
     451
     452        /*Assign output pointers:*/
     453        *penthalpy=enthalpy;
     454}
     455/*}}}1*/

  • issm/trunk/src/c/objects/Materials/Matpar.h

    r8490 r8496  
    2323                double  beta;
    2424                double  meltingpoint;
     25                double  referencetemperature;
    2526                double  mixed_layer_capacity;
    2627                double  thermal_exchange_velocity;
    … …  
    7778                double GetBeta();
    7879                double GetMeltingPoint();
     80                double GetReferenceTemperature();
    7981                double GetGamma();
    8082                double GetKn();
    8183                double TMeltingPoint(double pressure);
     84                double PureIceEnthalpy(double pressure);
     85                double GetEnthalpyDiffusionParameter(double enthalpy,double pressure);
     86                void   EnthalpyToThermal(double* ptemperature,double* pwaterfraction,double enthalpy,double pressure);
     87                void   ThermalToEnthalpy(double* penthalpy,double temperature,double waterfraction,double pressure);
    8288                /*}}}*/
    8389
Note: See TracChangeset for help on using the changeset viewer.