Changeset 8487


Ignore:
Timestamp:
06/03/11 08:20:25 (14 years ago)
Author:
seroussi
Message:

trunk/src/c: work on enthalpy

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

Legend:

Unmodified
Added
Removed

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

    r8483 r8487  
    450450        EnthalpySolutionEnum,
    451451        EnthalpyAnalysisEnum,
    452         EnthalpyEnum
     452        EnthalpyEnum,
     453        WaterFractionEnum,
     454        ReferenceTemperatureEnum
    453455};
    454456

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

    r8483 r8487  
    394394                case EnthalpyAnalysisEnum : return "EnthalpyAnalysis";
    395395                case EnthalpyEnum : return "Enthalpy";
     396                case WaterFractionEnum : return "WaterFraction";
     397                case ReferenceTemperatureEnum : return "ReferenceTemperature";
    396398                default : return "unknown";
    397399

  • issm/trunk/src/c/modules/ModelProcessorx/CreateParameters.cpp

    r8392 r8487  
    6262        parameters->AddObject(new DoubleParam(BetaEnum,iomodel->beta));
    6363        parameters->AddObject(new DoubleParam(MeltingPointEnum,iomodel->meltingpoint));
     64        parameters->AddObject(new DoubleParam(ReferenceTemperatureEnum,iomodel->referencetemperature));
    6465        parameters->AddObject(new DoubleParam(LatentHeatEnum,iomodel->latentheat));
    6566        parameters->AddObject(new DoubleParam(HeatCapacityEnum,iomodel->heatcapacity));

  • issm/trunk/src/c/modules/ModelProcessorx/Enthalpy/UpdateElementsEnthalpy.cpp

    r8483 r8487  
    4040        IoModelFetchData(&iomodel->pressure,NULL,NULL,iomodel_handle,"pressure");
    4141        IoModelFetchData(&iomodel->temperature,NULL,NULL,iomodel_handle,"temperature");
     42        IoModelFetchData(&iomodel->waterfraction,NULL,NULL,iomodel_handle,"waterfraction");
    4243        IoModelFetchData(&iomodel->geothermalflux,NULL,NULL,iomodel_handle,"geothermalflux");
    4344        IoModelFetchData(&iomodel->vx,NULL,NULL,iomodel_handle,"vx");
    … …  
    7475        xfree((void**)&iomodel->pressure);
    7576        xfree((void**)&iomodel->temperature);
     77        xfree((void**)&iomodel->waterfraction);
    7678        xfree((void**)&iomodel->geothermalflux);
    7779        xfree((void**)&iomodel->vx);

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

    r8483 r8487  
    392392        else if (strcmp(name,"EnthalpyAnalysis")==0) return EnthalpyAnalysisEnum;
    393393        else if (strcmp(name,"Enthalpy")==0) return EnthalpyEnum;
     394        else if (strcmp(name,"WaterFraction")==0) return WaterFractionEnum;
     395        else if (strcmp(name,"ReferenceTemperature")==0) return ReferenceTemperatureEnum;
    394396        else _error_("Enum %s not found",name);
    395397

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

    r8483 r8487  
    48664866                for(i=0;i<6;i++)nodeinputs[i]=iomodel->temperature[penta_vertex_ids[i]-1];
    48674867                this->inputs->AddInput(new PentaVertexInput(TemperatureEnum,nodeinputs));
     4868        }
     4869        if (iomodel->waterfraction) {
     4870                for(i=0;i<6;i++)nodeinputs[i]=iomodel->waterfraction[penta_vertex_ids[i]-1];
     4871                this->inputs->AddInput(new PentaVertexInput(WaterFractionEnum,nodeinputs));
    48684872        }
    48694873        if (iomodel->dhdt) {
    … …  
    60106014                this->inputs->AddInput(new PentaVertexInput(TemperaturePicardEnum,values));
    60116015        }
     6016
     6017        /*Free ressources:*/
     6018        xfree((void**)&doflist);
     6019}
     6020/*}}}*/
     6021/*FUNCTION Penta::InputUpdateFromSolutionEnthalpy {{{1*/
     6022void  Penta::InputUpdateFromSolutionEnthalpy(double* solution){
     6023
     6024        const int    numdof=NDOF1*NUMVERTICES;
     6025
     6026        bool   converged;
     6027        int    i,rheology_law;
     6028        double xyz_list[NUMVERTICES][3];
     6029        double values[numdof];
     6030        double temperatures[numdof];
     6031        double waterfraction[numdof];
     6032        double B[numdof];
     6033        double B_average,s_average;
     6034        int*   doflist=NULL;
     6035
     6036        /*Get dof list: */
     6037        GetDofList(&doflist,NoneApproximationEnum,GsetEnum);
     6038
     6039        /*Use the dof list to index into the solution vector: */
     6040        for(i=0;i<numdof;i++){
     6041                values[i]=solution[doflist[i]];
     6042
     6043                /*Check solution*/
     6044                if(isnan(values[i])) _error_("NaN found in solution vector");
     6045        }
     6046
     6047        /*Get all inputs and parameters*/
     6048        GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES);
     6049        Input* surface_input=inputs->GetInput(SurfaceEnum); _assert_(surface_input);
     6050
     6051        this->inputs->GetParameterValue(&converged,ConvergedEnum);
     6052        if(converged){
     6053                /*Convert enthalpy into temperature and water fraction*/
     6054//              for(i=0;i<numdof;i++){
     6055//                      if values[i]<GetPureIceEnthalpy(){
     6056//                              temperatures[i]=values[i];
     6057//                              waterfraction[i]=0;
     6058//                      }
     6059//                      else{
     6060//                              temperatures[i]=values[i];
     6061//                              waterfraction[i]=values[i];
     6062//                      }
     6063//              }
     6064
     6065                this->inputs->AddInput(new PentaVertexInput(EnthalpyEnum,values));
     6066                this->inputs->AddInput(new PentaVertexInput(WaterFractionEnum,waterfraction));
     6067                this->inputs->AddInput(new PentaVertexInput(TemperatureEnum,temperatures));
     6068
     6069                /*Update Rheology only if converged (we must make sure that the temperature is below melting point
     6070                 * otherwise the rheology could be negative*/
     6071                this->parameters->FindParam(&rheology_law,RheologyLawEnum);
     6072                switch(rheology_law){
     6073                        case NoneEnum:
     6074                                /*Do nothing: B is not temperature dependent*/
     6075                                break;
     6076                        case PatersonEnum:
     6077                                B_average=Paterson((values[0]+values[1]+values[2]+values[3]+values[4]+values[5])/6.0);
     6078                                for(i=0;i<numdof;i++) B[i]=B_average;
     6079                                this->matice->inputs->AddInput(new PentaVertexInput(RheologyBEnum,B));
     6080                                break;
     6081                        case ArrheniusEnum:
     6082                                surface_input->GetParameterAverage(&s_average);
     6083                                B_average=Arrhenius((values[0]+values[1]+values[2]+values[3]+values[4]+values[5])/6.0,
     6084                                                        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),
     6085                                                        matice->GetN());
     6086                                for(i=0;i<numdof;i++) B[i]=B_average;
     6087                                this->matice->inputs->AddInput(new PentaVertexInput(RheologyBEnum,B));
     6088                                break;
     6089                        default:
     6090                                _error_("Rheology law %s not supported yet",EnumToStringx(rheology_law));
     6091
     6092                }
     6093        }
     6094//      else{
     6095//              this->inputs->AddInput(new PentaVertexInput(EnthalpyPicardEnum,values));
     6096//      }
    60126097
    60136098        /*Free ressources:*/
    … …  
    71117196                        this->inputs->AddInput(new PentaVertexInput(VyMeshEnum,nodeinputs));
    71127197                        this->inputs->AddInput(new PentaVertexInput(VzMeshEnum,nodeinputs));
     7198                        if (iomodel->temperature && iomodel->waterfraction) {
     7199                                for(i=0;i<6;i++){
     7200                                        if(iomodel->temperature[penta_vertex_ids[i]-1] < iomodel->meltingpoint-iomodel->beta*iomodel->pressure[penta_vertex_ids[i]-1]){
     7201                                                nodeinputs[i]=iomodel->heatcapacity*(iomodel->temperature[penta_vertex_ids[i]-1]-iomodel->referencetemperature);
     7202                                        }
     7203                                        else nodeinputs[i]=iomodel->heatcapacity*
     7204                                         (iomodel->meltingpoint-iomodel->beta*iomodel->pressure[penta_vertex_ids[i]-1]-iomodel->referencetemperature)
     7205                                                +iomodel->latentheat*iomodel->waterfraction[penta_vertex_ids[i]-1];
     7206                                }
     7207                                this->inputs->AddInput(new PentaVertexInput(EnthalpyEnum,nodeinputs));
     7208                        }
     7209                        else _error_("temperature and waterfraction required for the enthalpy solution");
    71137210                        break;
    71147211

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

    r8483 r8487  
    245245                void    InputUpdateFromSolutionPrognostic(double* solutiong);
    246246                void    InputUpdateFromSolutionThermal( double* solutiong);
     247                void    InputUpdateFromSolutionEnthalpy( double* solutiong);
    247248                void    InputUpdateFromSolutionOneDof(double* solutiong,int enum_type);
    248249                void    InputUpdateFromSolutionOneDofCollapsed(double* solutiong,int enum_type);

  • issm/trunk/src/c/objects/IoModel.cpp

    r8399 r8487  
    6969        xfree((void**)&this->pressure);
    7070        xfree((void**)&this->temperature);
     71        xfree((void**)&this->waterfraction);
    7172        xfree((void**)&this->drag_coefficient);
    7273        xfree((void**)&this->drag_p);
    … …  
    213214        IoModelFetchData(&this->hydro_kn,iomodel_handle,"hydro_kn");
    214215        IoModelFetchData(&this->meltingpoint,iomodel_handle,"meltingpoint");
     216        IoModelFetchData(&this->referencetemperature,iomodel_handle,"referencetemperature");
    215217        IoModelFetchData(&this->latentheat,iomodel_handle,"latentheat");
    216218        IoModelFetchData(&this->heatcapacity,iomodel_handle,"heatcapacity");
    … …  
    286288        this->pressure=NULL;
    287289        this->temperature=NULL;
     290        this->waterfraction=NULL;
    288291        this->gl_melting_rate=0;
    289292        this->basal_melting_rate=NULL;
    … …  
    383386        this->beta=0;
    384387        this->meltingpoint=0;
     388        this->referencetemperature=0;
    385389        this->latentheat=0;
    386390        this->heatcapacity=0;

  • issm/trunk/src/c/objects/IoModel.h

    r8399 r8487  
    5656                double* pressure;
    5757                double* temperature;
     58                double* waterfraction;
    5859
    5960                /*i/o: */
    … …  
    172173                double hydro_kn;
    173174                double meltingpoint;
     175                double referencetemperature;
    174176                double latentheat;
    175177                double  heatcapacity,thermalconductivity;
Note: See TracChangeset for help on using the changeset viewer.