Changeset 16734


Ignore:
Timestamp:
11/13/13 11:02:01 (11 years ago)
Author:
seroussi
Message:

CHG: adding Enthalpy GetInputFromSolution

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

Legend:

Unmodified
Added
Removed

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

    r16684 r16734  
    188188}/*}}}*/
    189189void EnthalpyAnalysis::InputUpdateFromSolution(IssmDouble* solution,Element* element){/*{{{*/
    190         _error_("not implemented yet");
    191 }/*}}}*/
     190
     191        bool        converged;
     192        int         i,rheology_law;
     193        IssmDouble  B_average,s_average,T_average=0.,P_average=0.;
     194        int        *doflist   = NULL;
     195        IssmDouble *xyz_list  = NULL;
     196
     197        /*Fetch number of nodes and dof for this finite element*/
     198        int numnodes    = element->GetNumberOfNodes();
     199
     200        /*Fetch dof list and allocate solution vector*/
     201        element->GetDofList(&doflist,NoneApproximationEnum,GsetEnum);
     202        IssmDouble* values        = xNew<IssmDouble>(numnodes);
     203        IssmDouble* pressure      = xNew<IssmDouble>(numnodes);
     204        IssmDouble* temperature   = xNew<IssmDouble>(numnodes);
     205        IssmDouble* waterfraction = xNew<IssmDouble>(numnodes);
     206
     207        /*Use the dof list to index into the solution vector: */
     208        for(i=0;i<numnodes;i++){
     209                values[i]=solution[doflist[i]];
     210
     211                /*Check solution*/
     212                if(xIsNan<IssmDouble>(values[i])) _error_("NaN found in solution vector");
     213        }
     214
     215        /*Get all inputs and parameters*/
     216        element->GetInputValue(&converged,ConvergedEnum);
     217        element->GetInputListOnNodes(pressure,PressureEnum);
     218        if(converged){
     219                for(i=0;i<numnodes;i++){
     220                        element->EnthalpyToThermal(&temperature[i],&waterfraction[i],values[i],pressure[i]);
     221                        if(waterfraction[i]<0.) _error_("Negative water fraction found in solution vector");
     222                        //if(waterfraction[i]>1.) _error_("Water fraction >1 found in solution vector");
     223                }
     224                element->AddInput(EnthalpyEnum,values,P1Enum);
     225                element->AddInput(WaterfractionEnum,waterfraction,P1Enum);
     226                element->AddInput(TemperatureEnum,temperature,P1Enum);
     227
     228                /*Update Rheology only if converged (we must make sure that the temperature is below melting point
     229                 * otherwise the rheology could be negative*/
     230                element->FindParam(&rheology_law,MaterialsRheologyLawEnum);
     231                switch(rheology_law){
     232                        case NoneEnum:
     233                                /*Do nothing: B is not temperature dependent*/
     234                                break;
     235                        case PatersonEnum:
     236                                for(i=0;i<numnodes;i++) T_average+=values[i]/reCast<IssmDouble>(numnodes);
     237                                B_average=Paterson(T_average);
     238                                element->AddMaterialInput(MaterialsRheologyBEnum,&B_average,P0Enum);
     239                                break;
     240                        case ArrheniusEnum:{
     241                                Input* surface_input=element->GetInput(SurfaceEnum); _assert_(surface_input);
     242                                surface_input->GetInputAverage(&s_average);
     243                                element->GetVerticesCoordinates(&xyz_list);
     244                                for(i=0;i<numnodes;i++) T_average+=values[i]/reCast<IssmDouble>(numnodes);
     245                                //B_average=Arrhenius(T_average,
     246                                                        //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),
     247                                                        //element->GetMaticeParameter(MaterialsRheologyNEnum));
     248                                element->AddMaterialInput(MaterialsRheologyBEnum,&B_average,P0Enum);
     249                                break;
     250                                }
     251                        case LliboutryDuvalEnum:
     252                                for(i=0;i<numnodes;i++) T_average+=values[i]/reCast<IssmDouble>(numnodes);
     253                                for(i=0;i<numnodes;i++) P_average+=pressure[i]/reCast<IssmDouble>(numnodes);
     254                                B_average=LliboutryDuval(T_average,P_average,
     255                                                        element->GetMaterialParameter(MaterialsRheologyNEnum),
     256                                                        element->GetMaterialParameter(MaterialsBetaEnum),
     257                                                        element->GetMaterialParameter(ConstantsReferencetemperatureEnum),
     258                                                        element->GetMaterialParameter(MaterialsHeatcapacityEnum),
     259                                                        element->GetMaterialParameter(MaterialsLatentheatEnum));
     260                                element->AddMaterialInput(MaterialsRheologyBEnum,&B_average,P0Enum);
     261                                break;
     262                        default:
     263                                _error_("Rheology law " << EnumToStringx(rheology_law) << " not supported yet");
     264                }
     265        }
     266        else{
     267                element->AddInput(EnthalpyPicardEnum,values,P1Enum);
     268        }
     269
     270        /*Free ressources:*/
     271        xDelete<IssmDouble>(values);
     272        xDelete<IssmDouble>(pressure);
     273        xDelete<IssmDouble>(temperature);
     274        xDelete<IssmDouble>(waterfraction);
     275        xDelete<IssmDouble>(xyz_list);
     276        xDelete<int>(doflist);
     277}/*}}}*/

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

    r16727 r16734  
    4646                virtual void   CreatePVector(Vector<IssmDouble>* pf)=0;
    4747                virtual void   CreateJacobianMatrix(Matrix<IssmDouble>* Jff)=0;
     48                virtual void   EnthalpyToThermal(IssmDouble* ptemperature,IssmDouble* pwaterfraction,IssmDouble enthalpy,IssmDouble pressure)=0;
    4849                virtual void   FindParam(int* pvalue,int paramenum)=0;
    4950                virtual void   FindParam(IssmDouble* pvalue,int paramenum)=0;

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

    r16727 r16734  
    963963void Penta::Echo(void){
    964964        this->DeepEcho();
     965}
     966/*}}}*/
     967/*FUNCTION Penta::EnthalpyToThermal{{{*/
     968void Penta::EnthalpyToThermal(IssmDouble* ptemperature,IssmDouble* pwaterfraction,IssmDouble enthalpy,IssmDouble pressure){
     969        matpar->EnthalpyToThermal(ptemperature,pwaterfraction,enthalpy,pressure);
    965970}
    966971/*}}}*/

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

    r16727 r16734  
    8484                void   CreateJacobianMatrix(Matrix<IssmDouble>* Jff);
    8585                void   Delta18oParameterization(void);
     86                void   EnthalpyToThermal(IssmDouble* ptemperature,IssmDouble* pwaterfraction,IssmDouble enthalpy,IssmDouble pressure);
    8687                void     GetDofList(int** pdoflist,int approximation_enum,int setenum);
    8788                void     GetDofListVelocity(int** pdoflist,int setenum);

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

    r16727 r16734  
    8383                void        CreateJacobianMatrix(Matrix<IssmDouble>* Jff){_error_("not implemented yet");};
    8484                void        Delta18oParameterization(void){_error_("not implemented yet");};
     85                void        EnthalpyToThermal(IssmDouble* ptemperature,IssmDouble* pwaterfraction,IssmDouble enthalpy,IssmDouble pressure){_error_("not implemented yet");};
    8586                void        FindParam(int* pvalue,int paramenum);
    8687                void        FindParam(IssmDouble* pvalue,int paramenum);

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

    r16727 r16734  
    7979                void        CreateJacobianMatrix(Matrix<IssmDouble>* Jff);
    8080                void        Delta18oParameterization(void);
     81                void        EnthalpyToThermal(IssmDouble* ptemperature,IssmDouble* pwaterfraction,IssmDouble enthalpy,IssmDouble pressure){_error_("not implemented yet");};
    8182                void        FindParam(int* pvalue,int paramenum);
    8283                void        FindParam(IssmDouble* pvalue,int paramenum);
Note: See TracChangeset for help on using the changeset viewer.