Changeset 16888


Ignore:
Timestamp:
11/22/13 10:24:39 (11 years ago)
Author:
Mathieu Morlighem
Message:

CHG: added Enthalpy but not working

Location:
issm/trunk-jpl/src/c/analyses
Files:
2 edited

Legend:

Unmodified
Added
Removed

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

    r16837 r16888  
    185185/*Finite Element Analysis*/
    186186ElementMatrix* EnthalpyAnalysis::CreateKMatrix(Element* element){/*{{{*/
    187         _error_("not implemented yet");
     187
     188        /*compute all stiffness matrices for this element*/
     189        ElementMatrix* Ke1=CreateKMatrixVolume(element);
     190        ElementMatrix* Ke2=CreateKMatrixShelf(element);
     191        ElementMatrix* Ke =new ElementMatrix(Ke1,Ke2);
     192
     193        /*clean-up and return*/
     194        delete Ke1;
     195        delete Ke2;
     196        return Ke;
     197}/*}}}*/
     198ElementMatrix* EnthalpyAnalysis::CreateKMatrixVolume(Element* element){/*{{{*/
     199
     200        /*Intermediaries */
     201        int         stabilization;
     202        IssmDouble  Jdet,dt,u,v,w,um,vm,wm,vel;
     203        IssmDouble  h,hx,hy,hz,vx,vy,vz;
     204        IssmDouble  tau_parameter,diameter;
     205        IssmDouble  D_scalar;
     206        IssmDouble* xyz_list = NULL;
     207
     208        /*Fetch number of nodes and dof for this finite element*/
     209        int numnodes    = element->GetNumberOfNodes();
     210
     211        /*Initialize Element vector and other vectors*/
     212        ElementMatrix* Ke       = element->NewElementMatrix();
     213        IssmDouble*    basis    = xNew<IssmDouble>(numnodes);
     214        IssmDouble*    dbasis   = xNew<IssmDouble>(3*numnodes);
     215        IssmDouble*    B        = xNew<IssmDouble>(3*numnodes);
     216        IssmDouble*    Bprime   = xNew<IssmDouble>(3*numnodes);
     217        IssmDouble     D[3][3]  = {0.};
     218        IssmDouble     K[3][3];
     219
     220        /*Retrieve all inputs and parameters*/
     221        element->GetVerticesCoordinates(&xyz_list);
     222        element->FindParam(&dt,TimesteppingTimeStepEnum);
     223        element->FindParam(&stabilization,ThermalStabilizationEnum);
     224        IssmDouble  rho_water           = element->GetMaterialParameter(MaterialsRhoWaterEnum);
     225        IssmDouble  rho_ice             = element->GetMaterialParameter(MaterialsRhoIceEnum);
     226        IssmDouble  gravity             = element->GetMaterialParameter(ConstantsGEnum);
     227        IssmDouble  heatcapacity        = element->GetMaterialParameter(MaterialsHeatcapacityEnum);
     228        IssmDouble  thermalconductivity = element->GetMaterialParameter(MaterialsThermalconductivityEnum);
     229        Input* vx_input  = element->GetInput(VxEnum);     _assert_(vx_input);
     230        Input* vy_input  = element->GetInput(VyEnum);     _assert_(vy_input);
     231        Input* vz_input  = element->GetInput(VzEnum);     _assert_(vz_input);
     232        Input* vxm_input = element->GetInput(VxMeshEnum); _assert_(vxm_input);
     233        Input* vym_input = element->GetInput(VyMeshEnum); _assert_(vym_input);
     234        Input* vzm_input = element->GetInput(VzMeshEnum); _assert_(vzm_input);
     235        if(stabilization==2) diameter=element->MinEdgeLength(xyz_list);
     236
     237        /*Enthalpy diffusion parameter*/
     238        IssmDouble kappa=this->EnthalpyDiffusionParameterVolume(element); _assert_(kappa>0.);
     239
     240        /* Start  looping on the number of gaussian points: */
     241        Gauss* gauss=element->NewGauss(2);
     242        for(int ig=gauss->begin();ig<gauss->end();ig++){
     243                gauss->GaussPoint(ig);
     244
     245                element->JacobianDeterminant(&Jdet,xyz_list,gauss);
     246                D_scalar=gauss->weight*Jdet;
     247                if(dt!=0.) D_scalar=D_scalar*dt;
     248
     249                /*Conduction: */
     250                GetBConduct(B,element,xyz_list,gauss);
     251                D[0][0]=D_scalar*kappa/rho_ice;
     252                D[1][1]=D_scalar*kappa/rho_ice;
     253                D[2][2]=D_scalar*kappa/rho_ice;
     254                TripleMultiply(B,3,numnodes,1,
     255                                        &D[0][0],3,3,0,
     256                                        B,3,numnodes,0,
     257                                        &Ke->values[0],1);
     258
     259                /*Advection: */
     260                GetBAdvec(B,element,xyz_list,gauss);
     261                GetBAdvecprime(Bprime,element,xyz_list,gauss);
     262                vx_input->GetInputValue(&u,gauss); vxm_input->GetInputValue(&um,gauss); vx=u-um;
     263                vy_input->GetInputValue(&v,gauss); vym_input->GetInputValue(&vm,gauss); vy=v-vm;
     264                vz_input->GetInputValue(&w,gauss); vzm_input->GetInputValue(&wm,gauss); vz=w-wm;
     265                D[0][0]=D_scalar*vx;
     266                D[1][1]=D_scalar*vy;
     267                D[2][2]=D_scalar*vz;
     268                TripleMultiply(B,3,numnodes,1,
     269                                        &D[0][0],3,3,0,
     270                                        Bprime,3,numnodes,0,
     271                                        &Ke->values[0],1);
     272
     273                /*Transient: */
     274                if(dt!=0.){
     275                        D_scalar=gauss->weight*Jdet;
     276                        element->NodalFunctions(basis,gauss);
     277                        TripleMultiply(basis,numnodes,1,0,
     278                                                &D_scalar,1,1,0,
     279                                                basis,1,numnodes,0,
     280                                                &Ke->values[0],1);
     281                        D_scalar=D_scalar*dt;
     282                }
     283
     284                /*Artifficial diffusivity*/
     285                if(stabilization==1){
     286                        element->ElementSizes(&hx,&hy,&hz);
     287                        vel=sqrt(vx*vx + vy*vy + vz*vz)+1.e-14;
     288                        h=sqrt( pow(hx*vx/vel,2) + pow(hy*vy/vel,2) + pow(hz*vz/vel,2));
     289                        K[0][0]=h/(2.*vel)*fabs(vx*vx);  K[0][1]=h/(2.*vel)*fabs(vx*vy); K[0][2]=h/(2.*vel)*fabs(vx*vz);
     290                        K[1][0]=h/(2.*vel)*fabs(vy*vx);  K[1][1]=h/(2.*vel)*fabs(vy*vy); K[1][2]=h/(2.*vel)*fabs(vy*vz);
     291                        K[2][0]=h/(2.*vel)*fabs(vz*vx);  K[2][1]=h/(2.*vel)*fabs(vz*vy); K[2][2]=h/(2.*vel)*fabs(vz*vz);
     292                        for(int i=0;i<3;i++) for(int j=0;j<3;j++) K[i][j] = D_scalar*K[i][j];
     293
     294                        GetBAdvecprime(Bprime,element,xyz_list,gauss);
     295
     296                        TripleMultiply(Bprime,3,numnodes,1,
     297                                                &K[0][0],3,3,0,
     298                                                Bprime,3,numnodes,0,
     299                                                &Ke->values[0],1);
     300                }
     301                else if(stabilization==2){
     302                        element->NodalFunctionsDerivatives(dbasis,xyz_list,gauss);
     303                        tau_parameter=element->StabilizationParameter(u-um,v-vm,w-wm,diameter,kappa/rho_ice);
     304                        for(int i=0;i<numnodes;i++){
     305                                for(int j=0;j<numnodes;j++){
     306                                        Ke->values[i*numnodes+j]+=tau_parameter*D_scalar*
     307                                          ((u-um)*dbasis[0*3+i]+(v-vm)*dbasis[1*3+i]+(w-wm)*dbasis[2*3+i])*((u-um)*dbasis[0*3+j]+(v-vm)*dbasis[1*3+j]+(w-wm)*dbasis[2*3+j]);
     308                                }
     309                        }
     310                        if(dt!=0.){
     311                                for(int i=0;i<numnodes;i++){
     312                                        for(int j=0;j<numnodes;j++){
     313                                                Ke->values[i*numnodes+j]+=tau_parameter*D_scalar*basis[j]*((u-um)*dbasis[0*3+i]+(v-vm)*dbasis[1*3+i]+(w-wm)*dbasis[2*3+i]);
     314                                        }
     315                                }
     316                        }
     317                }
     318        }
     319
     320        /*Clean up and return*/
     321        xDelete<IssmDouble>(xyz_list);
     322        xDelete<IssmDouble>(basis);
     323        xDelete<IssmDouble>(dbasis);
     324        xDelete<IssmDouble>(B);
     325        xDelete<IssmDouble>(Bprime);
     326        delete gauss;
     327        return Ke;
     328}/*}}}*/
     329ElementMatrix* EnthalpyAnalysis::CreateKMatrixShelf(Element* element){/*{{{*/
     330
     331        /*Initialize Element matrix and return if necessary*/
     332        if(!element->IsOnBed() || !element->IsFloating()) return NULL;
     333
     334        IssmDouble  dt,Jdet,D;
     335        IssmDouble *xyz_list_base = NULL;
     336
     337        /*Get basal element*/
     338        if(!element->IsOnBed() || !element->IsFloating()) return NULL;
     339
     340        /*Fetch number of nodes for this finite element*/
     341        int numnodes = element->GetNumberOfNodes();
     342
     343        /*Initialize vectors*/
     344        ElementMatrix* Ke    = element->NewElementMatrix();
     345        IssmDouble*    basis = xNew<IssmDouble>(numnodes);
     346
     347        /*Retrieve all inputs and parameters*/
     348        element->GetVerticesCoordinatesBase(&xyz_list_base);
     349        element->FindParam(&dt,TimesteppingTimeStepEnum);
     350        IssmDouble  gravity             = element->GetMaterialParameter(ConstantsGEnum);
     351        IssmDouble  rho_water           = element->GetMaterialParameter(MaterialsRhoWaterEnum);
     352        IssmDouble  rho_ice             = element->GetMaterialParameter(MaterialsRhoIceEnum);
     353        IssmDouble  heatcapacity        = element->GetMaterialParameter(MaterialsHeatcapacityEnum);
     354        IssmDouble  mixed_layer_capacity= element->GetMaterialParameter(MaterialsMixedLayerCapacityEnum);
     355        IssmDouble  thermal_exchange_vel= element->GetMaterialParameter(MaterialsThermalExchangeVelocityEnum);
     356
     357        /* Start  looping on the number of gaussian points: */
     358        Gauss* gauss=element->NewGaussBase(2);
     359        for(int ig=gauss->begin();ig<gauss->end();ig++){
     360                gauss->GaussPoint(ig);
     361
     362                element->JacobianDeterminantBase(&Jdet,xyz_list_base,gauss);
     363                element->NodalFunctions(basis,gauss);
     364
     365                D=gauss->weight*Jdet*rho_water*mixed_layer_capacity*thermal_exchange_vel/(heatcapacity*rho_ice);
     366                if(reCast<bool,IssmDouble>(dt)) D=dt*D;
     367                TripleMultiply(basis,numnodes,1,0,
     368                                        &D,1,1,0,
     369                                        basis,1,numnodes,0,
     370                                        &Ke->values[0],1);
     371
     372        }
     373
     374        /*Clean up and return*/
     375        delete gauss;
     376        xDelete<IssmDouble>(basis);
     377        xDelete<IssmDouble>(xyz_list_base);
     378        return Ke;
    188379}/*}}}*/
    189380ElementVector* EnthalpyAnalysis::CreatePVector(Element* element){/*{{{*/
    … …  
    289480}/*}}}*/
    290481ElementVector* EnthalpyAnalysis::CreatePVectorSheet(Element* element){/*{{{*/
    291         return NULL;
     482
     483        /* Geothermal flux on ice sheet base and basal friction */
     484        if(!element->IsOnBed() || element->IsFloating()) return NULL;
     485
     486        IssmDouble  dt,Jdet,enthalpy,pressure,watercolumn,geothermalflux,vx,vy,vz;
     487        IssmDouble  enthalpyup,pressureup,alpha2,scalar,basalfriction,heatflux;
     488        IssmDouble *xyz_list_base = NULL;
     489
     490        /*Fetch number of nodes for this finite element*/
     491        int numnodes = element->GetNumberOfNodes();
     492
     493        /*Initialize vectors*/
     494        ElementVector* pe    = element->NewElementVector();
     495        IssmDouble*    basis = xNew<IssmDouble>(numnodes);
     496
     497        /*Retrieve all inputs and parameters*/
     498        element->GetVerticesCoordinatesBase(&xyz_list_base);
     499        element->FindParam(&dt,TimesteppingTimeStepEnum);
     500        Input* vx_input             = element->GetInput(VxEnum);                          _assert_(vx_input);
     501        Input* vy_input             = element->GetInput(VyEnum);                          _assert_(vy_input);
     502        Input* vz_input             = element->GetInput(VzEnum);                          _assert_(vz_input);
     503        Input* enthalpy_input       = element->GetInput(EnthalpyPicardEnum);              _assert_(enthalpy_input);
     504        Input* pressure_input       = element->GetInput(PressureEnum);                    _assert_(pressure_input);
     505        Input* geothermalflux_input = element->GetInput(BasalforcingsGeothermalfluxEnum); _assert_(geothermalflux_input);
     506        Input* watercolumn_input    = element->GetInput(WatercolumnEnum);                 _assert_(watercolumn_input);
     507        IssmDouble  rho_ice             = element->GetMaterialParameter(MaterialsRhoIceEnum);
     508        IssmDouble  heatcapacity        = element->GetMaterialParameter(MaterialsHeatcapacityEnum);
     509
     510        /*Build friction element, needed later: */
     511        Friction* friction=new Friction(element,3);
     512
     513        /* Start  looping on the number of gaussian points: */
     514        Gauss* gauss   = element->NewGaussBase(2);
     515        Gauss* gaussup = element->NewGaussTop(2);
     516        for(int ig=gauss->begin();ig<gauss->end();ig++){
     517                gauss->GaussPoint(ig);
     518
     519                element->JacobianDeterminantBase(&Jdet,xyz_list_base,gauss);
     520                element->NodalFunctions(basis,gauss);
     521
     522                enthalpy_input->GetInputValue(&enthalpy,gauss);
     523                pressure_input->GetInputValue(&pressure,gauss);
     524                watercolumn_input->GetInputValue(&watercolumn,gauss);
     525
     526                if((watercolumn<=0.) && (enthalpy < PureIceEnthalpy(element,pressure))){
     527                        /* the above check is equivalent to
     528                         NOT ((watercolumn>0.) AND (enthalpy<PIE)) AND (enthalpy<PIE)*/
     529                        geothermalflux_input->GetInputValue(&geothermalflux,gauss);
     530
     531                        friction->GetAlpha2(&alpha2,gauss,vx_input,vy_input,vz_input);
     532                        vx_input->GetInputValue(&vx,gauss);
     533                        vy_input->GetInputValue(&vy,gauss);
     534                        vz_input->GetInputValue(&vz,gauss);
     535                        basalfriction = alpha2*(vx*vx + vy*vy + vz*vz);
     536                        heatflux      = (basalfriction+geothermalflux)/(rho_ice);
     537
     538                        scalar = gauss->weight*Jdet*heatflux;
     539                        if(dt!=0.) scalar=dt*scalar;
     540
     541                        for(int i=0;i<numnodes;i++) pe->values[i]+=scalar*basis[i];
     542                }
     543                else if(enthalpy >= PureIceEnthalpy(element,pressure)){
     544                        /* check positive thickness of temperate basal ice layer */
     545                        enthalpy_input->GetInputValue(&enthalpyup,gaussup);
     546                        pressure_input->GetInputValue(&pressureup,gaussup);
     547                        if(enthalpyup >= PureIceEnthalpy(element,pressureup)){
     548                                // TODO: temperate ice has positive thickness: grad enthalpy*n=0.
     549                        }
     550                        else{
     551                                // only base temperate, set Dirichlet BCs in Penta::UpdateBasalConstraintsEnthalpy()
     552                        }
     553                }
     554                else{
     555                        // base cold, but watercolumn positive. Set base to h_pmp.
     556                }
     557        }
     558
     559        /*Clean up and return*/
     560        delete gauss;
     561        delete gaussup;
     562        delete friction;
     563        xDelete<IssmDouble>(basis);
     564        xDelete<IssmDouble>(xyz_list_base);
     565        return pe;
     566
    292567}/*}}}*/
    293568ElementVector* EnthalpyAnalysis::CreatePVectorShelf(Element* element){/*{{{*/
     569
     570        /*Get basal element*/
     571        if(!element->IsOnBed() || !element->IsFloating()) return NULL;
    294572
    295573        IssmDouble  h_pmp,dt,Jdet,scalar_ocean,pressure;
    296574        IssmDouble *xyz_list_base = NULL;
    297 
    298         /*Get basal element*/
    299         if(!element->IsOnBed() || !element->IsFloating()) return NULL;
    300575
    301576        /*Fetch number of nodes for this finite element*/
    … …  
    339614        xDelete<IssmDouble>(xyz_list_base);
    340615        return pe;
     616}/*}}}*/
     617void EnthalpyAnalysis::GetBConduct(IssmDouble* B,Element* element,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/
     618        /*Compute B  matrix. B=[B1 B2 B3 B4 B5 B6] where Bi is of size 5*NDOF1.
     619         * For node i, Bi' can be expressed in the actual coordinate system
     620         * by:
     621         *       Bi_conduct=[ dh/dx ]
     622         *                  [ dh/dy ]
     623         *                  [ dh/dz ]
     624         * where h is the interpolation function for node i.
     625         *
     626         * We assume B has been allocated already, of size: 3x(NDOF1*numnodes)
     627         */
     628
     629        /*Fetch number of nodes for this finite element*/
     630        int numnodes = element->GetNumberOfNodes();
     631
     632        /*Get nodal functions derivatives*/
     633        IssmDouble* dbasis=xNew<IssmDouble>(3*numnodes);
     634        element->NodalFunctionsDerivatives(dbasis,xyz_list,gauss);
     635
     636        /*Build B: */
     637        for(int i=0;i<numnodes;i++){
     638                B[numnodes*0+i] = dbasis[0*numnodes+i];
     639                B[numnodes*1+i] = dbasis[1*numnodes+i];
     640                B[numnodes*2+i] = dbasis[2*numnodes+i];
     641        }
     642
     643        /*Clean-up*/
     644        xDelete<IssmDouble>(dbasis);
     645}/*}}}*/
     646void EnthalpyAnalysis::GetBAdvec(IssmDouble* B,Element* element,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/
     647        /*Compute B  matrix. B=[B1 B2 B3 B4 B5 B6] where Bi is of size 5*NDOF1.
     648         * For node i, Bi' can be expressed in the actual coordinate system
     649         * by:
     650         *       Bi_advec =[ h ]
     651         *                 [ h ]
     652         *                 [ h ]
     653         * where h is the interpolation function for node i.
     654         *
     655         * We assume B has been allocated already, of size: 3x(NDOF1*NUMNODESP1)
     656         */
     657
     658        /*Fetch number of nodes for this finite element*/
     659        int numnodes = element->GetNumberOfNodes();
     660
     661        /*Get nodal functions*/
     662        IssmDouble* basis=xNew<IssmDouble>(numnodes);
     663        element->NodalFunctions(basis,gauss);
     664
     665        /*Build B: */
     666        for(int i=0;i<numnodes;i++){
     667                B[numnodes*0+i] = basis[i];
     668                B[numnodes*1+i] = basis[i];
     669                B[numnodes*2+i] = basis[i];
     670        }
     671
     672        /*Clean-up*/
     673        xDelete<IssmDouble>(basis);
     674}/*}}}*/
     675void EnthalpyAnalysis::GetBAdvecprime(IssmDouble* B,Element* element,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/
     676        /*Compute B  matrix. B=[B1 B2 B3 B4 B5 B6] where Bi is of size 5*NDOF1.
     677         * For node i, Bi' can be expressed in the actual coordinate system
     678         * by:
     679         *       Biprime_advec=[ dh/dx ]
     680         *                     [ dh/dy ]
     681         *                     [ dh/dz ]
     682         * where h is the interpolation function for node i.
     683         *
     684         * We assume B has been allocated already, of size: 3x(NDOF1*numnodes)
     685         */
     686
     687        /*Fetch number of nodes for this finite element*/
     688        int numnodes = element->GetNumberOfNodes();
     689
     690        /*Get nodal functions derivatives*/
     691        IssmDouble* dbasis=xNew<IssmDouble>(3*numnodes);
     692        element->NodalFunctionsDerivatives(dbasis,xyz_list,gauss);
     693
     694        /*Build B: */
     695        for(int i=0;i<numnodes;i++){
     696                B[numnodes*0+i] = dbasis[0*numnodes+i];
     697                B[numnodes*1+i] = dbasis[1*numnodes+i];
     698                B[numnodes*2+i] = dbasis[2*numnodes+i];
     699        }
     700
     701        /*Clean-up*/
     702        xDelete<IssmDouble>(dbasis);
    341703}/*}}}*/
    342704void EnthalpyAnalysis::GetSolutionFromInputs(Vector<IssmDouble>* solution,Element* element){/*{{{*/
    … …  
    422784        xDelete<int>(doflist);
    423785}/*}}}*/
     786
     787/*Intermediaries*/
     788IssmDouble EnthalpyAnalysis::EnthalpyDiffusionParameter(Element* element,IssmDouble enthalpy,IssmDouble pressure){/*{{{*/
     789
     790        IssmDouble heatcapacity             = element->GetMaterialParameter(MaterialsHeatcapacityEnum);
     791        IssmDouble temperateiceconductivity = element->GetMaterialParameter(MaterialsTemperateiceconductivityEnum);
     792        IssmDouble thermalconductivity      = element->GetMaterialParameter(MaterialsThermalconductivityEnum);
     793
     794        if(enthalpy < PureIceEnthalpy(element,pressure)){
     795                return thermalconductivity/heatcapacity;
     796        }
     797        else{
     798                return temperateiceconductivity/heatcapacity;
     799        }
     800}/*}}}*/
     801IssmDouble EnthalpyAnalysis::EnthalpyDiffusionParameterVolume(Element* element){/*{{{*/
     802
     803        int         iv;
     804        IssmDouble  lambda;                   /* fraction of cold ice    */
     805        IssmDouble  kappa   ,kappa_c,kappa_t; /* enthalpy conductivities */
     806        IssmDouble  Hc,Ht;
     807
     808
     809        /*Get pressures and enthalpies on vertices*/
     810        int         numvertices = element->GetNumberOfVertices();
     811        IssmDouble* pressures   = xNew<IssmDouble>(numvertices);
     812        IssmDouble* enthalpies  = xNew<IssmDouble>(numvertices);
     813        IssmDouble* PIE         = xNew<IssmDouble>(numvertices);
     814        IssmDouble* dHpmp       = xNew<IssmDouble>(numvertices);
     815        element->GetInputListOnVertices(pressures,PressureEnum);
     816        element->GetInputListOnVertices(enthalpies,EnthalpyEnum);
     817        for(iv=0;iv<numvertices;iv++){
     818                PIE[iv]   = PureIceEnthalpy(element,pressures[iv]);
     819                dHpmp[iv] = enthalpies[iv]-PIE[iv];
     820        }
     821
     822        bool allequalsign = true;
     823        if(dHpmp[0]<0.){
     824                for(iv=1; iv<numvertices;iv++) allequalsign=(allequalsign && (dHpmp[iv]<0.));
     825        }
     826        else{
     827                for(iv=1; iv<numvertices;iv++) allequalsign=(allequalsign && (dHpmp[iv]>=0.));
     828        }
     829
     830        if(allequalsign){
     831                kappa = EnthalpyDiffusionParameter(element,enthalpies[0],pressures[0]);
     832        }
     833        else{
     834                /* return harmonic mean of thermal conductivities, weighted by fraction of cold/temperate ice,
     835                        cf Patankar 1980, pp44 */
     836                kappa_c = EnthalpyDiffusionParameter(element,PureIceEnthalpy(element,0.)-1.,0.);
     837                kappa_t = EnthalpyDiffusionParameter(element,PureIceEnthalpy(element,0.)+1.,0.);
     838                Hc=0.; Ht=0.;
     839                for(iv=0; iv<numvertices;iv++){
     840                        if(enthalpies[iv]<PIE[iv])
     841                         Hc+=(PIE[iv]-enthalpies[iv]);
     842                        else
     843                         Ht+=(enthalpies[iv]-PIE[iv]);
     844                }
     845                _assert_((Hc+Ht)>0.);
     846                lambda = Hc/(Hc+Ht);
     847                kappa  = 1./(lambda/kappa_c + (1.-lambda)/kappa_t);
     848        }
     849
     850        /*Clean up and return*/
     851        xDelete<IssmDouble>(PIE);
     852        xDelete<IssmDouble>(dHpmp);
     853        xDelete<IssmDouble>(pressures);
     854        xDelete<IssmDouble>(enthalpies);
     855        return kappa;
     856}
     857/*}}}*/
     858IssmDouble EnthalpyAnalysis::PureIceEnthalpy(Element* element,IssmDouble pressure){/*{{{*/
     859
     860        IssmDouble heatcapacity         = element->GetMaterialParameter(MaterialsHeatcapacityEnum);
     861        IssmDouble referencetemperature = element->GetMaterialParameter(ConstantsReferencetemperatureEnum);
     862
     863        return heatcapacity*(TMeltingPoint(element,pressure)-referencetemperature);
     864}/*}}}*/
     865IssmDouble EnthalpyAnalysis::TMeltingPoint(Element* element,IssmDouble pressure){/*{{{*/
     866
     867        IssmDouble meltingpoint = element->GetMaterialParameter(MaterialsMeltingpointEnum);
     868        IssmDouble beta         = element->GetMaterialParameter(MaterialsBetaEnum);
     869
     870        return meltingpoint-beta*pressure;
     871}/*}}}*/

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

    r16812 r16888  
    2222                /*Finite element Analysis*/
    2323                ElementMatrix* CreateKMatrix(Element* element);
     24                ElementMatrix* CreateKMatrixVolume(Element* element);
     25                ElementMatrix* CreateKMatrixShelf(Element* element);
    2426                ElementVector* CreatePVector(Element* element);
    2527                ElementVector* CreatePVectorVolume(Element* element);
    2628                ElementVector* CreatePVectorSheet(Element* element);
    2729                ElementVector* CreatePVectorShelf(Element* element);
     30                void GetBConduct(IssmDouble* B,Element* element,IssmDouble* xyz_list,Gauss* gauss);
     31                void GetBAdvec(IssmDouble* B,Element* element,IssmDouble* xyz_list,Gauss* gauss);
     32                void GetBAdvecprime(IssmDouble* B,Element* element,IssmDouble* xyz_list,Gauss* gauss);
    2833                void GetSolutionFromInputs(Vector<IssmDouble>* solution,Element* element);
    2934                void InputUpdateFromSolution(IssmDouble* solution,Element* element);
     35
     36                /*Intermediaries*/
     37                IssmDouble EnthalpyDiffusionParameter(Element* element,IssmDouble enthalpy,IssmDouble pressure);
     38                IssmDouble EnthalpyDiffusionParameterVolume(Element* element);
     39                IssmDouble PureIceEnthalpy(Element* element,IssmDouble pressure);
     40                IssmDouble TMeltingPoint(Element* element,IssmDouble pressure);
    3041};
    3142#endif
Note: See TracChangeset for help on using the changeset viewer.