Changeset 16812


Ignore:
Timestamp:
11/17/13 08:55:19 (11 years ago)
Author:
Mathieu Morlighem
Message:

NEW: working on wnthalpy

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

Legend:

Unmodified
Added
Removed

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

    r16782 r16812  
    188188}/*}}}*/
    189189ElementVector* EnthalpyAnalysis::CreatePVector(Element* element){/*{{{*/
    190 _error_("not implemented yet");
     190
     191        /*compute all load vectors for this element*/
     192        ElementVector* pe1=CreatePVectorVolume(element);
     193        ElementVector* pe2=CreatePVectorSheet(element);
     194        ElementVector* pe3=CreatePVectorShelf(element);
     195        ElementVector* pe =new ElementVector(pe1,pe2,pe3);
     196
     197        /*clean-up and return*/
     198        delete pe1;
     199        delete pe2;
     200        delete pe3;
     201        return pe;
     202}/*}}}*/
     203ElementVector* EnthalpyAnalysis::CreatePVectorVolume(Element* element){/*{{{*/
     204
     205        /*Intermediaries*/
     206        int         stabilization;
     207        IssmDouble  Jdet,phi,dt;
     208        IssmDouble  enthalpy;
     209        IssmDouble  kappa,tau_parameter,diameter;
     210        IssmDouble  u,v,w;
     211        IssmDouble  scalar_def,scalar_transient;
     212        IssmDouble* xyz_list = NULL;
     213
     214        /*Fetch number of nodes and dof for this finite element*/
     215        int numnodes    = element->GetNumberOfNodes();
     216        int numvertices = element->GetNumberOfVertices();
     217
     218        /*Initialize Element vector*/
     219        ElementVector* pe             = element->NewElementVector();
     220        IssmDouble*    basis          = xNew<IssmDouble>(numnodes);
     221        IssmDouble*    dbasis         = xNew<IssmDouble>(3*numnodes);
     222        IssmDouble*    pressure       = xNew<IssmDouble>(numvertices);
     223        IssmDouble*    enthalpypicard = xNew<IssmDouble>(numvertices);
     224
     225        /*Retrieve all inputs and parameters*/
     226        element->GetVerticesCoordinates(&xyz_list);
     227        IssmDouble  rho_ice             = element->GetMaterialParameter(MaterialsRhoIceEnum);
     228        IssmDouble  thermalconductivity = element->GetMaterialParameter(MaterialsThermalconductivityEnum);
     229        element->FindParam(&dt,TimesteppingTimeStepEnum);
     230        element->FindParam(&stabilization,ThermalStabilizationEnum);
     231        Input* vx_input=element->GetInput(VxEnum); _assert_(vx_input);
     232        Input* vy_input=element->GetInput(VyEnum); _assert_(vy_input);
     233        Input* vz_input=element->GetInput(VzEnum); _assert_(vz_input);
     234        Input* enthalpy_input = NULL;
     235        if(reCast<bool,IssmDouble>(dt)){enthalpy_input = element->GetInput(EnthalpyEnum); _assert_(enthalpy_input);}
     236        if(stabilization==2){
     237                element->GetInputListOnVertices(enthalpypicard,EnthalpyPicardEnum);
     238                element->GetInputListOnVertices(pressure,PressureEnum);
     239                diameter=element->MinEdgeLength(xyz_list);
     240        }
     241
     242        /* Start  looping on the number of gaussian points: */
     243        Gauss* gauss=element->NewGauss(3);
     244        for(int ig=gauss->begin();ig<gauss->end();ig++){
     245                gauss->GaussPoint(ig);
     246
     247                element->JacobianDeterminant(&Jdet,xyz_list,gauss);
     248                element->NodalFunctions(basis,gauss);
     249                element->ViscousHeating(&phi,xyz_list,gauss,vx_input,vy_input,vz_input);
     250
     251                scalar_def=phi/rho_ice*Jdet*gauss->weight;
     252                if(reCast<bool,IssmDouble>(dt)) scalar_def=scalar_def*dt;
     253
     254                /*TODO: add -beta*laplace T_m(p)*/
     255                for(int i=0;i<numnodes;i++) pe->values[i]+=scalar_def*basis[i];
     256
     257                /* Build transient now */
     258                if(reCast<bool,IssmDouble>(dt)){
     259                        enthalpy_input->GetInputValue(&enthalpy, gauss);
     260                        scalar_transient=enthalpy*Jdet*gauss->weight;
     261                        for(int i=0;i<numnodes;i++) pe->values[i]+=scalar_transient*basis[i];
     262                }
     263
     264                if(stabilization==2){
     265                        element->NodalFunctionsDerivatives(dbasis,xyz_list,gauss);
     266
     267                        vx_input->GetInputValue(&u,gauss);
     268                        vy_input->GetInputValue(&v,gauss);
     269                        vz_input->GetInputValue(&w,gauss);
     270                        kappa          = element->EnthalpyDiffusionParameterVolume(numvertices,enthalpypicard,pressure) / rho_ice;
     271                        tau_parameter  = element->StabilizationParameter(u,v,w,diameter,kappa);
     272
     273                        for(int i=0;i<numnodes;i++) pe->values[i]+=tau_parameter*scalar_def*(u*dbasis[0*3+i]+v*dbasis[1*3+i]+w*dbasis[2*3+i]);
     274                        if(reCast<bool,IssmDouble>(dt)){
     275                                for(int i=0;i<numnodes;i++) pe->values[i]+=tau_parameter*scalar_transient*(u*dbasis[0*3+i]+v*dbasis[1*3+i]+w*dbasis[2*3+i]);
     276                        }
     277                }
     278        }
     279
     280        /*Clean up and return*/
     281        xDelete<IssmDouble>(basis);
     282        xDelete<IssmDouble>(dbasis);
     283        xDelete<IssmDouble>(pressure);
     284        xDelete<IssmDouble>(enthalpypicard);
     285        xDelete<IssmDouble>(xyz_list);
     286        delete gauss;
     287        return pe;
     288
     289}/*}}}*/
     290ElementVector* EnthalpyAnalysis::CreatePVectorSheet(Element* element){/*{{{*/
     291        return NULL;
     292}/*}}}*/
     293ElementVector* EnthalpyAnalysis::CreatePVectorShelf(Element* element){/*{{{*/
     294
     295        IssmDouble  t_pmp,dt,Jdet,scalar_ocean,pressure;
     296        IssmDouble *xyz_list_base = NULL;
     297
     298        /*Get basal element*/
     299        if(!element->IsOnBed() || !element->IsFloating()) return NULL;
     300
     301        /*Fetch number of nodes for this finite element*/
     302        int numnodes = element->GetNumberOfNodes();
     303
     304        /*Initialize vectors*/
     305        ElementVector* pe    = element->NewElementVector();
     306        IssmDouble*    basis = xNew<IssmDouble>(numnodes);
     307
     308        /*Retrieve all inputs and parameters*/
     309        element->GetVerticesCoordinatesBase(&xyz_list_base);
     310        element->FindParam(&dt,TimesteppingTimeStepEnum);
     311        Input*      pressure_input=element->GetInput(PressureEnum); _assert_(pressure_input);
     312        IssmDouble  gravity             = element->GetMaterialParameter(ConstantsGEnum);
     313        IssmDouble  rho_water           = element->GetMaterialParameter(MaterialsRhoWaterEnum);
     314        IssmDouble  rho_ice             = element->GetMaterialParameter(MaterialsRhoIceEnum);
     315        IssmDouble  heatcapacity        = element->GetMaterialParameter(MaterialsHeatcapacityEnum);
     316        IssmDouble  mixed_layer_capacity= element->GetMaterialParameter(MaterialsMixedLayerCapacityEnum);
     317        IssmDouble  thermal_exchange_vel= element->GetMaterialParameter(MaterialsThermalExchangeVelocityEnum);
     318
     319        /* Start  looping on the number of gaussian points: */
     320        Gauss* gauss=element->NewGaussBase(2);
     321        for(int ig=gauss->begin();ig<gauss->end();ig++){
     322                gauss->GaussPoint(ig);
     323
     324                element->JacobianDeterminantBase(&Jdet,xyz_list_base,gauss);
     325                element->NodalFunctions(basis,gauss);
     326
     327                pressure_input->GetInputValue(&pressure,gauss);
     328                t_pmp=element->TMeltingPoint(pressure);
     329
     330                scalar_ocean=gauss->weight*Jdet*rho_water*mixed_layer_capacity*thermal_exchange_vel*(t_pmp)/(heatcapacity*rho_ice);
     331                if(reCast<bool,IssmDouble>(dt)) scalar_ocean=dt*scalar_ocean;
     332
     333                for(int i=0;i<numnodes;i++) pe->values[i]+=scalar_ocean*basis[i];
     334        }
     335
     336        /*Clean up and return*/
     337        delete gauss;
     338        xDelete<IssmDouble>(basis);
     339        xDelete<IssmDouble>(xyz_list_base);
     340        return pe;
    191341}/*}}}*/
    192342void EnthalpyAnalysis::GetSolutionFromInputs(Vector<IssmDouble>* solution,Element* element){/*{{{*/

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

    r16782 r16812  
    2323                ElementMatrix* CreateKMatrix(Element* element);
    2424                ElementVector* CreatePVector(Element* element);
     25                ElementVector* CreatePVectorVolume(Element* element);
     26                ElementVector* CreatePVectorSheet(Element* element);
     27                ElementVector* CreatePVectorShelf(Element* element);
    2528                void GetSolutionFromInputs(Vector<IssmDouble>* solution,Element* element);
    2629                void InputUpdateFromSolution(IssmDouble* solution,Element* element);

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

    r16810 r16812  
    4949                virtual void   DeleteMaterials(void)=0;
    5050                virtual void   EnthalpyToThermal(IssmDouble* ptemperature,IssmDouble* pwaterfraction,IssmDouble enthalpy,IssmDouble pressure)=0;
     51                virtual IssmDouble EnthalpyDiffusionParameter(IssmDouble enthalpy,IssmDouble pressure)=0;
     52                virtual IssmDouble EnthalpyDiffusionParameterVolume(int numvertices,IssmDouble* enthalpy,IssmDouble* pressure)=0;
    5153                virtual void   FindParam(int* pvalue,int paramenum)=0;
    5254                virtual void   FindParam(IssmDouble* pvalue,int paramenum)=0;

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

    r16810 r16812  
    873873void Penta::EnthalpyToThermal(IssmDouble* ptemperature,IssmDouble* pwaterfraction,IssmDouble enthalpy,IssmDouble pressure){
    874874        matpar->EnthalpyToThermal(ptemperature,pwaterfraction,enthalpy,pressure);
     875}
     876/*}}}*/
     877/*FUNCTION Penta::EnthalpyDiffusionParameter{{{*/
     878IssmDouble Penta::EnthalpyDiffusionParameter(IssmDouble enthalpy,IssmDouble pressure){
     879        return matpar->GetEnthalpyDiffusionParameter(enthalpy,pressure);
     880}
     881/*}}}*/
     882/*FUNCTION Penta::EnthalpyDiffusionParameterVolume{{{*/
     883IssmDouble Penta::EnthalpyDiffusionParameterVolume(int numvertices,IssmDouble* enthalpy,IssmDouble* pressure){
     884        return matpar->GetEnthalpyDiffusionParameterVolume(numvertices,enthalpy,pressure);
    875885}
    876886/*}}}*/
    … …  
    37813791                GetInputListOnVertices(&enthalpy[0],EnthalpyPicardEnum);
    37823792                GetInputListOnVertices(&pressure[0],PressureEnum);
    3783                 kappa=matpar->GetEnthalpyDiffusionParameterVolume(enthalpy,pressure); _assert_(kappa>0.);
     3793                kappa=matpar->GetEnthalpyDiffusionParameterVolume(NUMVERTICES,enthalpy,pressure); _assert_(kappa>0.);
    37843794
    37853795                D_scalar_conduct=gauss->weight*Jdet*kappa/rho_ice;
    … …  
    42504260                        GetInputListOnVertices(&enthalpypicard[0],EnthalpyPicardEnum);
    42514261                        GetInputListOnVertices(&pressure_p[0],PressureEnum);
    4252                         kappa=matpar->GetEnthalpyDiffusionParameterVolume(enthalpypicard,pressure_p);
     4262                        kappa=matpar->GetEnthalpyDiffusionParameterVolume(NUMVERTICES,enthalpypicard,pressure_p);
    42534263                        kappa/=rho_ice;
    42544264                        tau_parameter=StabilizationParameter(u,v,w,diameter,kappa);
    … …  
    48084818                        else{
    48094819                                enthalpy_input->GetInputDerivativeValue(&d1enthalpy[0],&xyz_list[0][0],gauss);
    4810                                 kappa=matpar->GetEnthalpyDiffusionParameterVolume(enthalpy, pressure); _assert_(kappa>0.);
     4820                                kappa=matpar->GetEnthalpyDiffusionParameterVolume(NUMVERTICES,&enthalpy[0],&pressure[0]); _assert_(kappa>0.);
    48114821                                for(i=0;i<3;i++) vec_heatflux[i]=-kappa*d1enthalpy[i];
    48124822                        }

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

    r16810 r16812  
    204204                ElementVector* CreatePVectorFreeSurfaceBase(void);
    205205                ElementVector* CreatePVectorL2ProjectionBase(void);
     206                IssmDouble     EnthalpyDiffusionParameter(IssmDouble enthalpy,IssmDouble pressure);
     207                IssmDouble     EnthalpyDiffusionParameterVolume(int numvertices,IssmDouble* enthalpy,IssmDouble* pressure);
    206208                void           GetAreaCoordinates(IssmDouble *area_coordinates,IssmDouble xyz_zero[3][3],IssmDouble xyz_list[6][3],int numpoints);
    207209

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

    r16810 r16812  
    8484                void        Delta18oParameterization(void){_error_("not implemented yet");};
    8585                void        EnthalpyToThermal(IssmDouble* ptemperature,IssmDouble* pwaterfraction,IssmDouble enthalpy,IssmDouble pressure){_error_("not implemented yet");};
     86                IssmDouble  EnthalpyDiffusionParameter(IssmDouble enthalpy,IssmDouble pressure){_error_("not implemented");};
     87                IssmDouble  EnthalpyDiffusionParameterVolume(int numvertices,IssmDouble* enthalpy,IssmDouble* pressure){_error_("not implemented");};
    8688                void        FindParam(int* pvalue,int paramenum);
    8789                void        FindParam(IssmDouble* pvalue,int paramenum);

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

    r16810 r16812  
    239239                ElementVector* CreatePVectorL2Projection(void);
    240240                ElementVector* CreatePVectorL2ProjectionBase(void);
     241                IssmDouble     EnthalpyDiffusionParameter(IssmDouble enthalpy,IssmDouble pressure){_error_("not implemented");};
     242                IssmDouble     EnthalpyDiffusionParameterVolume(int numvertices,IssmDouble* enthalpy,IssmDouble* pressure){_error_("not implemented");};
    241243                IssmDouble     GetArea(void);
    242244                void           GetAreaCoordinates(IssmDouble *area_coordinates,IssmDouble xyz_zero[3][3],IssmDouble xyz_list[3][3],int numpoints);

  • issm/trunk-jpl/src/c/classes/Materials/Matpar.cpp

    r16783 r16812  
    417417/*}}}*/
    418418/*FUNCTION Matpar::GetEnthalpyDiffusionParameterVolume{{{*/
    419 IssmDouble Matpar::GetEnthalpyDiffusionParameterVolume(IssmDouble enthalpy[6], IssmDouble pressure[6]){
    420 
    421         IssmDouble lambda; // fraction of cold ice
    422         IssmDouble kappa,kappa_c,kappa_t;  //enthalpy conductivities
    423         IssmDouble Hc,Ht;
    424         IssmDouble PIE[6],dHpmp[6];
    425         int iv;
    426 
    427         for (iv=0; iv<6; iv++){
     419IssmDouble Matpar::GetEnthalpyDiffusionParameterVolume(int numvertices,IssmDouble* enthalpy,IssmDouble* pressure){
     420
     421        int         iv;
     422        IssmDouble  lambda;                 // fraction of cold ice
     423        IssmDouble  kappa,kappa_c,kappa_t;  //enthalpy conductivities
     424        IssmDouble  Hc,Ht;
     425        IssmDouble* PIE   = xNew<IssmDouble>(numvertices);
     426        IssmDouble* dHpmp = xNew<IssmDouble>(numvertices);
     427
     428        for(iv=0; iv<numvertices; iv++){
    428429                PIE[iv]=PureIceEnthalpy(pressure[iv]);
    429430                dHpmp[iv]=enthalpy[iv]-PIE[iv];
    … …  
    432433        bool allequalsign=true;
    433434        if(dHpmp[0]<0)
    434                 for(iv=1; iv<6;iv++)
    435                         allequalsign=(allequalsign && (dHpmp[iv]<0));
     435                for(iv=1; iv<numvertices;iv++) allequalsign=(allequalsign && (dHpmp[iv]<0));
    436436        else
    437                 for(iv=1; iv<6;iv++)
    438                         allequalsign=(allequalsign && (dHpmp[iv]>=0));
     437                for(iv=1; iv<numvertices;iv++) allequalsign=(allequalsign && (dHpmp[iv]>=0));
    439438
    440439        if(allequalsign){
    … …  
    447446                kappa_t=GetEnthalpyDiffusionParameter(PureIceEnthalpy(0.)+1.,0.);
    448447                Hc=0.; Ht=0.;
    449                 for(iv=0; iv<6;iv++){
     448                for(iv=0; iv<numvertices;iv++){
    450449                        if(enthalpy[iv]<PIE[iv])
    451                                 Hc+=(PIE[iv]-enthalpy[iv]);
     450                         Hc+=(PIE[iv]-enthalpy[iv]);
    452451                        else
    453                                 Ht+=(enthalpy[iv]-PIE[iv]);
     452                         Ht+=(enthalpy[iv]-PIE[iv]);
    454453                }
    455454                _assert_((Hc+Ht)>0.);
    456                 lambda=Hc/(Hc+Ht);
    457                 kappa=1./(lambda/kappa_c + (1.-lambda)/kappa_t);
    458         }
     455                lambda = Hc/(Hc+Ht);
     456                kappa  = 1./(lambda/kappa_c + (1.-lambda)/kappa_t);
     457        }
     458
     459        /*Clean up and return*/
     460        xDelete<IssmDouble>(PIE);
     461        xDelete<IssmDouble>(dHpmp);
    459462        return kappa;
    460463}

  • issm/trunk-jpl/src/c/classes/Materials/Matpar.h

    r16783 r16812  
    129129                IssmDouble PureIceEnthalpy(IssmDouble pressure);
    130130                IssmDouble GetEnthalpyDiffusionParameter(IssmDouble enthalpy,IssmDouble pressure);
    131                 IssmDouble GetEnthalpyDiffusionParameterVolume(IssmDouble enthalpy[6], IssmDouble pressure[6]);
     131                IssmDouble GetEnthalpyDiffusionParameterVolume(int numvertices,IssmDouble* enthalpy,IssmDouble* pressure);
    132132                IssmDouble GetLithosphereShearModulus();
    133133                IssmDouble GetLithosphereDensity();
Note: See TracChangeset for help on using the changeset viewer.