Changeset 8483


Ignore:
Timestamp:
06/02/11 15:34:39 (14 years ago)
Author:
seroussi
Message:

c: start to create enthalpy solution

Location:
issm/trunk/src/c
Files:
6 added
12 edited

Legend:

Unmodified
Added
Removed

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

    r8429 r8483  
    447447        IsDiagnosticEnum,
    448448        IsPrognosticEnum,
    449         IsThermalEnum
     449        IsThermalEnum,
     450        EnthalpySolutionEnum,
     451        EnthalpyAnalysisEnum,
     452        EnthalpyEnum
    450453};
    451454

  • issm/trunk/src/c/Makefile.am

    r8468 r8483  
    463463                                        ./modules/ModelProcessorx/Thermal/CreateConstraintsThermal.cpp\
    464464                                        ./modules/ModelProcessorx/Thermal/CreateLoadsThermal.cpp\
     465                                        ./modules/ModelProcessorx/Enthalpy/UpdateElementsEnthalpy.cpp\
     466                                        ./modules/ModelProcessorx/Enthalpy/CreateNodesEnthalpy.cpp\
     467                                        ./modules/ModelProcessorx/Enthalpy/CreateConstraintsEnthalpy.cpp\
     468                                        ./modules/ModelProcessorx/Enthalpy/CreateLoadsEnthalpy.cpp\
    465469                                        ./modules/ModelProcessorx/Hydrology/UpdateElementsHydrology.cpp\
    466470                                        ./modules/ModelProcessorx/Hydrology/CreateNodesHydrology.cpp\
    … …  
    11211125                                        ./modules/ModelProcessorx/Thermal/CreateConstraintsThermal.cpp\
    11221126                                        ./modules/ModelProcessorx/Thermal/CreateLoadsThermal.cpp\
     1127                                        ./modules/ModelProcessorx/Enthalpy/UpdateElementsEnthalpy.cpp\
     1128                                        ./modules/ModelProcessorx/Enthalpy/CreateNodesEnthalpy.cpp\
     1129                                        ./modules/ModelProcessorx/Enthalpy/CreateConstraintsEnthalpy.cpp\
     1130                                        ./modules/ModelProcessorx/Enthalpy/CreateLoadsEnthalpy.cpp\
    11231131                                        ./modules/ModelProcessorx/Hydrology/UpdateElementsHydrology.cpp\
    11241132                                        ./modules/ModelProcessorx/Hydrology/CreateNodesHydrology.cpp\
    … …  
    13241332                                        ./solutions/thermal_core.cpp\
    13251333                                        ./solutions/thermal_core_step.cpp\
     1334                                        ./solutions/enthalpy_core.cpp\
    13261335                                        ./solutions/WriteLockFile.cpp\
    13271336                                        ./solutions/control_core.cpp\

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

    r8429 r8483  
    391391                case IsPrognosticEnum : return "IsPrognostic";
    392392                case IsThermalEnum : return "IsThermal";
     393                case EnthalpySolutionEnum : return "EnthalpySolution";
     394                case EnthalpyAnalysisEnum : return "EnthalpyAnalysis";
     395                case EnthalpyEnum : return "Enthalpy";
    393396                default : return "unknown";
    394397

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

    r8429 r8483  
    7474                        break;
    7575               
     76                case EnthalpyAnalysisEnum:
     77                        CreateNodesEnthalpy(pnodes, iomodel,iomodel_handle);
     78                        CreateConstraintsEnthalpy(pconstraints,iomodel,iomodel_handle);
     79                        CreateLoadsEnthalpy(ploads,iomodel,iomodel_handle);
     80                        UpdateElementsEnthalpy(elements,iomodel,iomodel_handle,analysis_counter,analysis_type);
     81                        break;
     82               
    7683                case HydrologyAnalysisEnum:
    7784                        CreateNodesHydrology(pnodes, iomodel,iomodel_handle);

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

    r8303 r8483  
    7575                numdofs=1;
    7676        }
     77        else if (analysis_type==EnthalpyAnalysisEnum){
     78                numdofs=1;
     79        }
    7780        else if (analysis_type==HydrologyAnalysisEnum){
    7881                numdofs=1;

  • issm/trunk/src/c/modules/ModelProcessorx/ModelProcessorx.h

    r8365 r8483  
    6161void    UpdateElementsThermal(Elements* elements,IoModel* iomodel,FILE* iomodel_handle,int analysis_counter,int analysis_type);
    6262
     63/*enthalpy:*/
     64void    CreateNodesEnthalpy(Nodes** pnodes,IoModel* iomodel,FILE* iomodel_handle);
     65void    CreateConstraintsEnthalpy(Constraints** pconstraints,IoModel* iomodel,FILE* iomodel_handle);
     66void  CreateLoadsEnthalpy(Loads** ploads, IoModel* iomodel, FILE* iomodel_handle);
     67void    UpdateElementsEnthalpy(Elements* elements,IoModel* iomodel,FILE* iomodel_handle,int analysis_counter,int analysis_type);
     68
    6369/*hydrology:*/
    6470void    CreateNodesHydrology(Nodes** pnodes,IoModel* iomodel,FILE* iomodel_handle);
    … …  
    6672void  CreateLoadsHydrology(Loads** ploads, IoModel* iomodel, FILE* iomodel_handle);
    6773void    UpdateElementsHydrology(Elements* elements,IoModel* iomodel,FILE* iomodel_handle,int analysis_counter,int analysis_type);
    68 
    6974
    7075/*melting:*/

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

    r8429 r8483  
    389389        else if (strcmp(name,"IsPrognostic")==0) return IsPrognosticEnum;
    390390        else if (strcmp(name,"IsThermal")==0) return IsThermalEnum;
     391        else if (strcmp(name,"EnthalpySolution")==0) return EnthalpySolutionEnum;
     392        else if (strcmp(name,"EnthalpyAnalysis")==0) return EnthalpyAnalysisEnum;
     393        else if (strcmp(name,"Enthalpy")==0) return EnthalpyEnum;
    391394        else _error_("Enum %s not found",name);
    392395

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

    r8432 r8483  
    18031803        }
    18041804
     1805        /*Clean up and return*/
     1806        delete gauss;
     1807        return Ke;
     1808}
     1809/*}}}*/
     1810/*FUNCTION Penta::CreateKMatrixEnthalpy {{{1*/
     1811ElementMatrix* Penta::CreateKMatrixEnthalpy(void){
     1812       
     1813        /*compute all stiffness matrices for this element*/
     1814        ElementMatrix* Ke1=CreateKMatrixEnthalpyVolume();
     1815        ElementMatrix* Ke2=CreateKMatrixEnthalpyShelf();
     1816        ElementMatrix* Ke =new ElementMatrix(Ke1,Ke2);
     1817       
     1818        /*clean-up and return*/
     1819        delete Ke1;
     1820        delete Ke2;
     1821        return Ke;
     1822}
     1823/*}}}*/
     1824/*FUNCTION Penta::CreateKMatrixEnthalpyVolume {{{1*/
     1825ElementMatrix* Penta::CreateKMatrixEnthalpyVolume(void){
     1826
     1827        /*Constants*/
     1828        const int    numdof=NDOF1*NUMVERTICES;
     1829
     1830        /*Intermediaries */
     1831        int        artdiff;
     1832        int        i,j,ig,found=0;
     1833        double     Jdet,u,v,w,um,vm,wm,epsvel;
     1834        double     gravity,rho_ice,rho_water;
     1835        double     heatcapacity,thermalconductivity,dt;
     1836        double     latentheat,kappa=1;
     1837        double     tau_parameter,diameter;
     1838        double     xyz_list[NUMVERTICES][3];
     1839        double     B[3][numdof];
     1840        double     Bprime[3][numdof];
     1841        double     B_conduct[3][numdof];
     1842        double     B_advec[3][numdof];
     1843        double     B_artdiff[2][numdof];
     1844        double     Bprime_advec[3][numdof];
     1845        double     L[numdof];
     1846        double     dh1dh6[3][6];
     1847        double     D_scalar_conduct,D_scalar_advec;
     1848        double     D_scalar_trans,D_scalar_artdiff;
     1849        double     D[3][3];
     1850        double     K[2][2]={0.0};
     1851        double     Ke_gaussian_conduct[numdof][numdof];
     1852        double     Ke_gaussian_advec[numdof][numdof];
     1853        double     Ke_gaussian_artdiff[numdof][numdof];
     1854        double     Ke_gaussian_transient[numdof][numdof];
     1855        Tria*      tria=NULL;
     1856        GaussPenta *gauss=NULL;
     1857
     1858        /*Initialize Element matrix*/
     1859        ElementMatrix* Ke=new ElementMatrix(nodes,NUMVERTICES,this->parameters,NoneApproximationEnum);
     1860
     1861        /*Retrieve all inputs and parameters*/
     1862        GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES);
     1863        rho_water=matpar->GetRhoWater();
     1864        rho_ice=matpar->GetRhoIce();
     1865        gravity=matpar->GetG();
     1866        heatcapacity=matpar->GetHeatCapacity();
     1867        latentheat=matpar->GetLatentHeat();
     1868        thermalconductivity=matpar->GetThermalConductivity();
     1869        this->parameters->FindParam(&dt,DtEnum);
     1870        this->parameters->FindParam(&artdiff,ArtDiffEnum);
     1871        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);
     1878        if (artdiff==2) diameter=MinEdgeLength(xyz_list);
     1879
     1880        /* Start  looping on the number of gaussian points: */
     1881        gauss=new GaussPenta(2,2);
     1882        for (ig=gauss->begin();ig<gauss->end();ig++){
     1883
     1884                gauss->GaussPoint(ig);
     1885
     1886                GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss);
     1887
     1888                /*Conduction: */ 
     1889                /*Need to change that depending on enthalpy value -> cold or temperate ice: */ 
     1890
     1891                GetBConduct(&B_conduct[0][0],&xyz_list[0][0],gauss);
     1892
     1893                //kappa=GetEnthalpyDiffusivityParameter(rho_ice,heatcapacity,thermalconductivity,latentheat,enthalpy);
     1894                D_scalar_conduct=gauss->weight*Jdet*kappa;
     1895                if(dt) D_scalar_conduct=D_scalar_conduct*dt;
     1896
     1897                D[0][0]=D_scalar_conduct; D[0][1]=0; D[0][2]=0;
     1898                D[1][0]=0; D[1][1]=D_scalar_conduct; D[1][2]=0;
     1899                D[2][0]=0; D[2][1]=0; D[2][2]=D_scalar_conduct;
     1900
     1901                TripleMultiply(&B_conduct[0][0],3,numdof,1,
     1902                                        &D[0][0],3,3,0,
     1903                                        &B_conduct[0][0],3,numdof,0,
     1904                                        &Ke_gaussian_conduct[0][0],0);
     1905
     1906                /*Advection: */
     1907
     1908                GetBAdvec(&B_advec[0][0],&xyz_list[0][0],gauss);
     1909                GetBprimeAdvec(&Bprime_advec[0][0],&xyz_list[0][0],gauss);
     1910
     1911                vx_input->GetParameterValue(&u, gauss);
     1912                vy_input->GetParameterValue(&v, gauss);
     1913                vz_input->GetParameterValue(&w, gauss);
     1914                vxm_input->GetParameterValue(&um,gauss);
     1915                vym_input->GetParameterValue(&vm,gauss);
     1916                vzm_input->GetParameterValue(&wm,gauss);
     1917
     1918                D_scalar_advec=gauss->weight*Jdet;
     1919                if(dt) D_scalar_advec=D_scalar_advec*dt;
     1920
     1921                D[0][0]=D_scalar_advec*(u-um);D[0][1]=0;                    D[0][2]=0;
     1922                D[1][0]=0;                    D[1][1]=D_scalar_advec*(v-vm);D[1][2]=0;
     1923                D[2][0]=0;                    D[2][1]=0;                    D[2][2]=D_scalar_advec*(w-wm);
     1924
     1925                TripleMultiply(&B_advec[0][0],3,numdof,1,
     1926                                        &D[0][0],3,3,0,
     1927                                        &Bprime_advec[0][0],3,numdof,0,
     1928                                        &Ke_gaussian_advec[0][0],0);
     1929
     1930                /*Transient: */
     1931
     1932                if(dt){
     1933                        GetNodalFunctionsP1(&L[0], gauss);
     1934                        D_scalar_trans=gauss->weight*Jdet;
     1935                        D_scalar_trans=D_scalar_trans;
     1936
     1937                        TripleMultiply(&L[0],numdof,1,0,
     1938                                                &D_scalar_trans,1,1,0,
     1939                                                &L[0],1,numdof,0,
     1940                                                &Ke_gaussian_transient[0][0],0);
     1941                }
     1942                else{
     1943                        for(i=0;i<numdof;i++) for(j=0;j<numdof;j++) Ke_gaussian_transient[i][j]=0;
     1944                }
     1945
     1946                /*Artifficial diffusivity*/
     1947
     1948                if(artdiff==1){
     1949                        /*Build K: */
     1950                        D_scalar_artdiff=gauss->weight*Jdet/(pow(u-um,2)+pow(v-vm,2)+epsvel);
     1951                        if(dt) D_scalar_artdiff=D_scalar_artdiff*dt;
     1952                        K[0][0]=D_scalar_artdiff*pow(u,2);       K[0][1]=D_scalar_artdiff*fabs(u)*fabs(v);
     1953                        K[1][0]=D_scalar_artdiff*fabs(u)*fabs(v);K[1][1]=D_scalar_artdiff*pow(v,2);
     1954
     1955                        GetBArtdiff(&B_artdiff[0][0],&xyz_list[0][0],gauss);
     1956
     1957                        TripleMultiply(&B_artdiff[0][0],2,numdof,1,
     1958                                                &K[0][0],2,2,0,
     1959                                                &B_artdiff[0][0],2,numdof,0,
     1960                                                &Ke_gaussian_artdiff[0][0],0);
     1961                }
     1962                else if(artdiff==2){
     1963
     1964                        GetNodalFunctionsP1Derivatives(&dh1dh6[0][0],&xyz_list[0][0], gauss);
     1965
     1966                        tau_parameter=GetStabilizationParameter(u-um,v-vm,w-wm,diameter,rho_ice,heatcapacity,thermalconductivity);
     1967
     1968                        for(i=0;i<numdof;i++){
     1969                                for(j=0;j<numdof;j++){
     1970                                        Ke_gaussian_artdiff[i][j]=tau_parameter*D_scalar_advec*
     1971                                          ((u-um)*dh1dh6[0][i]+(v-vm)*dh1dh6[1][i]+(w-wm)*dh1dh6[2][i])*((u-um)*dh1dh6[0][j]+(v-vm)*dh1dh6[1][j]+(w-wm)*dh1dh6[2][j]);
     1972                                }
     1973                        }
     1974                        if(dt){
     1975                                for(i=0;i<numdof;i++){
     1976                                        for(j=0;j<numdof;j++){
     1977                                                Ke_gaussian_artdiff[i][j]+=tau_parameter*D_scalar_trans*L[j]*((u-um)*dh1dh6[0][i]+(v-vm)*dh1dh6[1][i]+(w-wm)*dh1dh6[2][i]);
     1978                                        }
     1979                                }
     1980                        }
     1981                }
     1982                else{
     1983                        for(i=0;i<numdof;i++) for(j=0;j<numdof;j++) Ke_gaussian_artdiff[i][j]=0;
     1984                }
     1985
     1986                for(i=0;i<numdof;i++) for(j=0;j<numdof;j++) Ke->values[i*numdof+j]+=Ke_gaussian_conduct[i][j]+Ke_gaussian_advec[i][j]+Ke_gaussian_transient[i][j]+Ke_gaussian_artdiff[i][j];
     1987        }
     1988
     1989        /*Clean up and return*/
     1990        delete gauss;
     1991        return Ke;
     1992}
     1993/*}}}*/
     1994/*FUNCTION Penta::CreateKMatrixEnthalpyShelf {{{1*/
     1995ElementMatrix* Penta::CreateKMatrixEnthalpyShelf(void){
     1996
     1997        /*Constants*/
     1998        const int    numdof=NDOF1*NUMVERTICES;
     1999
     2000        /*Intermediaries */
     2001        int       i,j,ig;
     2002        double    mixed_layer_capacity,thermal_exchange_velocity;
     2003        double    rho_ice,rho_water,heatcapacity;
     2004        double    Jdet2d,dt;
     2005        double    xyz_list[NUMVERTICES][3];
     2006        double   xyz_list_tria[NUMVERTICES2D][3];
     2007        double    l1l6[NUMVERTICES];
     2008        double    D_scalar;
     2009        double    Ke_gaussian[numdof][numdof]={0.0};
     2010        GaussPenta *gauss=NULL;
     2011
     2012        /*Initialize Element matrix and return if necessary*/
     2013        if (!IsOnBed() || !IsOnShelf()) return NULL;
     2014        ElementMatrix* Ke=new ElementMatrix(nodes,NUMVERTICES,this->parameters,NoneApproximationEnum);
     2015
     2016        /*Retrieve all inputs and parameters*/
     2017        this->parameters->FindParam(&dt,DtEnum);
     2018        mixed_layer_capacity=matpar->GetMixedLayerCapacity();
     2019        thermal_exchange_velocity=matpar->GetThermalExchangeVelocity();
     2020        rho_water=matpar->GetRhoWater();
     2021        rho_ice=matpar->GetRhoIce();
     2022        heatcapacity=matpar->GetHeatCapacity();
     2023        GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES);
     2024        for(i=0;i<NUMVERTICES2D;i++) for(j=0;j<3;j++) xyz_list_tria[i][j]=xyz_list[i][j];
     2025
     2026        /* Start looping on the number of gauss (nodes on the bedrock) */
     2027        gauss=new GaussPenta(0,1,2,2);
     2028        for (ig=gauss->begin();ig<gauss->end();ig++){
     2029
     2030                gauss->GaussPoint(ig);
     2031               
     2032                GetTriaJacobianDeterminant(&Jdet2d, &xyz_list_tria[0][0], gauss);
     2033                GetNodalFunctionsP1(&l1l6[0], gauss);
     2034                               
     2035                D_scalar=gauss->weight*Jdet2d*rho_water*mixed_layer_capacity*thermal_exchange_velocity/(rho_ice*heatcapacity);
     2036                if(dt) D_scalar=dt*D_scalar;
     2037
     2038                TripleMultiply(&l1l6[0],numdof,1,0,
     2039                                        &D_scalar,1,1,0,
     2040                                        &l1l6[0],1,numdof,0,
     2041                                        &Ke_gaussian[0][0],0);
     2042
     2043                for(i=0;i<numdof;i++) for(j=0;j<numdof;j++) Ke->values[i*numdof+j]+=Ke_gaussian[i][j];
     2044        }
     2045       
    18052046        /*Clean up and return*/
    18062047        delete gauss;
    … …  
    30573298        /*Clean up and return*/
    30583299        delete gauss;
     3300        return pe;
     3301}
     3302/*}}}*/
     3303/*FUNCTION Penta::CreatePVectorEnthalpy {{{1*/
     3304ElementVector* Penta::CreatePVectorEnthalpy(void){
     3305
     3306        /*compute all load vectors for this element*/
     3307        ElementVector* pe1=CreatePVectorEnthalpyVolume();
     3308        ElementVector* pe2=CreatePVectorEnthalpySheet();
     3309        ElementVector* pe3=CreatePVectorEnthalpyShelf();
     3310        ElementVector* pe =new ElementVector(pe1,pe2,pe3);
     3311
     3312        /*clean-up and return*/
     3313        delete pe1;
     3314        delete pe2;
     3315        delete pe3;
     3316        return pe;
     3317}
     3318/*}}}*/
     3319/*FUNCTION Penta::CreatePVectorEnthalpyVolume {{{1*/
     3320ElementVector* Penta::CreatePVectorEnthalpyVolume(void){
     3321
     3322        /*Constants*/
     3323        const int  numdof=NUMVERTICES*NDOF1;
     3324
     3325        /*Intermediaries*/
     3326        int    i,j,ig,found=0;
     3327        int    friction_type,artdiff;
     3328        double Jdet,phi,dt;
     3329        double rho_ice,heatcapacity;
     3330        double thermalconductivity;
     3331        double viscosity,temperature;
     3332        double tau_parameter,diameter;
     3333        double u,v,w;
     3334        double scalar_def,scalar_transient;
     3335        double temperature_list[NUMVERTICES];
     3336        double xyz_list[NUMVERTICES][3];
     3337        double L[numdof];
     3338        double dh1dh6[3][6];
     3339        double epsilon[6];
     3340        GaussPenta *gauss=NULL;
     3341
     3342        /*Initialize Element vector*/
     3343        ElementVector* pe=new ElementVector(nodes,NUMVERTICES,this->parameters);
     3344
     3345        /*Retrieve all inputs and parameters*/
     3346        GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES);
     3347        rho_ice=matpar->GetRhoIce();
     3348        heatcapacity=matpar->GetHeatCapacity();
     3349        thermalconductivity=matpar->GetThermalConductivity();
     3350        this->parameters->FindParam(&dt,DtEnum);
     3351        this->parameters->FindParam(&artdiff,ArtDiffEnum);
     3352        Input* vx_input=inputs->GetInput(VxEnum); _assert_(vx_input);
     3353        Input* vy_input=inputs->GetInput(VyEnum); _assert_(vy_input);
     3354        Input* vz_input=inputs->GetInput(VzEnum); _assert_(vz_input);
     3355        Input* temperature_input=NULL;
     3356        if (dt) temperature_input=inputs->GetInput(TemperatureEnum); _assert_(inputs);
     3357        if (artdiff==2) diameter=MinEdgeLength(xyz_list);
     3358
     3359        /* Start  looping on the number of gaussian points: */
     3360        gauss=new GaussPenta(2,3);
     3361        for (ig=gauss->begin();ig<gauss->end();ig++){
     3362
     3363                gauss->GaussPoint(ig);
     3364
     3365                GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss);
     3366                GetNodalFunctionsP1(&L[0], gauss);
     3367
     3368                this->GetStrainRate3d(&epsilon[0],&xyz_list[0][0],gauss,vx_input,vy_input,vz_input);
     3369                matice->GetViscosity3dStokes(&viscosity,&epsilon[0]);
     3370                GetPhi(&phi, &epsilon[0], viscosity);
     3371
     3372                scalar_def=phi/(rho_ice)*Jdet*gauss->weight;
     3373                if(dt) scalar_def=scalar_def*dt;
     3374
     3375                for(i=0;i<NUMVERTICES;i++)  pe->values[i]+=scalar_def*L[i];
     3376
     3377                /* Build transient now */
     3378                if(dt){
     3379                        temperature_input->GetParameterValue(&temperature, gauss);
     3380                        scalar_transient=temperature*Jdet*gauss->weight;
     3381                        for(i=0;i<NUMVERTICES;i++)  pe->values[i]+=scalar_transient*L[i];
     3382                }
     3383
     3384                if(artdiff==2){
     3385                        GetNodalFunctionsP1Derivatives(&dh1dh6[0][0],&xyz_list[0][0], gauss);
     3386
     3387                        vx_input->GetParameterValue(&u, gauss);
     3388                        vy_input->GetParameterValue(&v, gauss);
     3389                        vz_input->GetParameterValue(&w, gauss);
     3390
     3391                        tau_parameter=GetStabilizationParameter(u,v,w,diameter,rho_ice,heatcapacity,thermalconductivity);
     3392
     3393                        for(i=0;i<NUMVERTICES;i++)  pe->values[i]+=tau_parameter*scalar_def*(u*dh1dh6[0][i]+v*dh1dh6[1][i]+w*dh1dh6[2][i]);
     3394                        if(dt){
     3395                                for(i=0;i<NUMVERTICES;i++)  pe->values[i]+=tau_parameter*scalar_transient*(u*dh1dh6[0][i]+v*dh1dh6[1][i]+w*dh1dh6[2][i]);
     3396                        }
     3397                }
     3398        }
     3399
     3400        /*Clean up and return*/
     3401        delete gauss;
     3402        return pe;
     3403}
     3404/*}}}*/
     3405/*FUNCTION Penta::CreatePVectorEnthalpyShelf {{{1*/
     3406ElementVector* Penta::CreatePVectorEnthalpyShelf(void){
     3407
     3408        /*Constants*/
     3409        const int  numdof=NUMVERTICES*NDOF1;
     3410
     3411        /*Intermediaries */
     3412        int        i,j,ig;
     3413        double     Jdet2d;
     3414        double     mixed_layer_capacity,thermal_exchange_velocity;
     3415        double     rho_ice,rho_water,pressure,dt,scalar_ocean;
     3416        double     meltingpoint,beta,heatcapacity,h_pmp;
     3417        double     xyz_list[NUMVERTICES][3];
     3418        double     xyz_list_tria[NUMVERTICES2D][3];
     3419        double     l1l6[NUMVERTICES];
     3420        GaussPenta* gauss=NULL;
     3421
     3422        /*Initialize Element vector*/
     3423        ElementVector* pe=new ElementVector(nodes,NUMVERTICES,this->parameters);
     3424
     3425        /* Ice/ocean heat exchange flux on ice shelf base */
     3426        if (!IsOnBed() || !IsOnShelf()) return NULL;
     3427
     3428        /*Retrieve all inputs and parameters*/
     3429        GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES);
     3430        for(i=0;i<NUMVERTICES2D;i++) for(j=0;j<3;j++) xyz_list_tria[i][j]=xyz_list[i][j];
     3431        mixed_layer_capacity=matpar->GetMixedLayerCapacity();
     3432        thermal_exchange_velocity=matpar->GetThermalExchangeVelocity();
     3433        rho_water=matpar->GetRhoWater();
     3434        rho_ice=matpar->GetRhoIce();
     3435        heatcapacity=matpar->GetHeatCapacity();
     3436        beta=matpar->GetBeta();
     3437        meltingpoint=matpar->GetMeltingPoint();
     3438        this->parameters->FindParam(&dt,DtEnum);
     3439        Input* pressure_input=inputs->GetInput(PressureEnum); _assert_(pressure_input);
     3440
     3441        /* Start looping on the number of gauss 2d (nodes on the bedrock) */
     3442        gauss=new GaussPenta(0,1,2,2);
     3443        for(ig=gauss->begin();ig<gauss->end();ig++){
     3444
     3445                gauss->GaussPoint(ig);
     3446
     3447                GetTriaJacobianDeterminant(&Jdet2d, &xyz_list_tria[0][0], gauss);
     3448                GetNodalFunctionsP1(&l1l6[0], gauss);
     3449
     3450                pressure_input->GetParameterValue(&pressure,gauss);
     3451                h_pmp=meltingpoint-beta*pressure;
     3452
     3453                scalar_ocean=gauss->weight*Jdet2d*rho_water*mixed_layer_capacity*thermal_exchange_velocity*(h_pmp)/(rho_ice*heatcapacity);
     3454                if(dt) scalar_ocean=dt*scalar_ocean;
     3455
     3456                for(i=0;i<numdof;i++) pe->values[i]+=scalar_ocean*l1l6[i];
     3457        }
     3458
     3459        /*Clean up and return*/
     3460        delete gauss;
     3461        return pe;
     3462}
     3463/*}}}*/
     3464/*FUNCTION Penta::CreatePVectorEnthalpySheet {{{1*/
     3465ElementVector* Penta::CreatePVectorEnthalpySheet(void){
     3466
     3467        /*Constants*/
     3468        const int  numdof=NUMVERTICES*NDOF1;
     3469
     3470        /*Intermediaries */
     3471        int        i,j,ig;
     3472        int        analysis_type,drag_type;
     3473        double     xyz_list[NUMVERTICES][3];
     3474        double     xyz_list_tria[NUMVERTICES2D][3]={0.0};
     3475        double     Jdet2d,dt;
     3476        double     rho_ice,heatcapacity,geothermalflux_value;
     3477        double     basalfriction,alpha2,vx,vy;
     3478        double     scalar;
     3479        double     l1l6[NUMVERTICES];
     3480        Friction*  friction=NULL;
     3481        GaussPenta* gauss=NULL;
     3482
     3483        /* Geothermal flux on ice sheet base and basal friction */
     3484        if (!IsOnBed() || IsOnShelf()) return NULL;
     3485
     3486        /*Initialize Element vector*/
     3487        ElementVector* pe=new ElementVector(nodes,NUMVERTICES,this->parameters);
     3488
     3489        /*Retrieve all inputs and parameters*/
     3490        GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES);
     3491        for(i=0;i<NUMVERTICES2D;i++) for(j=0;j<2;j++) xyz_list_tria[i][j]=xyz_list[i][j];
     3492        parameters->FindParam(&analysis_type,AnalysisTypeEnum);
     3493        rho_ice=matpar->GetRhoIce();
     3494        heatcapacity=matpar->GetHeatCapacity();
     3495        this->parameters->FindParam(&dt,DtEnum);
     3496        Input* vx_input=inputs->GetInput(VxEnum);                         _assert_(vx_input);
     3497        Input* vy_input=inputs->GetInput(VyEnum);                         _assert_(vy_input);
     3498        Input* vz_input=inputs->GetInput(VzEnum);                         _assert_(vz_input);
     3499        Input* geothermalflux_input=inputs->GetInput(GeothermalFluxEnum); _assert_(geothermalflux_input);
     3500
     3501        /*Build frictoin element, needed later: */
     3502        inputs->GetParameterValue(&drag_type,DragTypeEnum);
     3503        if (drag_type!=2)_error_(" non-viscous friction not supported yet!");
     3504        friction=new Friction("3d",inputs,matpar,analysis_type);
     3505
     3506        /* Start looping on the number of gauss 2d (nodes on the bedrock) */
     3507        gauss=new GaussPenta(0,1,2,2);
     3508        for(ig=gauss->begin();ig<gauss->end();ig++){
     3509
     3510                gauss->GaussPoint(ig);
     3511
     3512                GetTriaJacobianDeterminant(&Jdet2d, &xyz_list_tria[0][0], gauss);
     3513                GetNodalFunctionsP1(&l1l6[0], gauss);
     3514
     3515                geothermalflux_input->GetParameterValue(&geothermalflux_value,gauss);
     3516                friction->GetAlpha2(&alpha2,gauss,VxEnum,VyEnum,VzEnum);
     3517                vx_input->GetParameterValue(&vx,gauss);
     3518                vy_input->GetParameterValue(&vy,gauss);
     3519                basalfriction=alpha2*(pow(vx,2.0)+pow(vy,2.0));
     3520               
     3521                scalar=gauss->weight*Jdet2d*(basalfriction+geothermalflux_value)/(rho_ice);
     3522                if(dt) scalar=dt*scalar;
     3523
     3524                for(i=0;i<numdof;i++) pe->values[i]+=scalar*l1l6[i];
     3525        }
     3526
     3527        /*Clean up and return*/
     3528        delete gauss;
     3529        delete friction;
    30593530        return pe;
    30603531}
    … …  
    66347105                        break;
    66357106
     7107                case EnthalpyAnalysisEnum:
     7108                        /*Initialize mesh velocity*/
     7109                        for(i=0;i<6;i++)nodeinputs[i]=0;
     7110                        this->inputs->AddInput(new PentaVertexInput(VxMeshEnum,nodeinputs));
     7111                        this->inputs->AddInput(new PentaVertexInput(VyMeshEnum,nodeinputs));
     7112                        this->inputs->AddInput(new PentaVertexInput(VzMeshEnum,nodeinputs));
     7113                        break;
     7114
    66367115                default:
    66377116                        /*No update for other solution types*/

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

    r8376 r8483  
    167167                ElementMatrix* CreateKMatrixDiagnosticVertVolume(void);
    168168                ElementMatrix* CreateKMatrixDiagnosticVertSurface(void);
     169                ElementMatrix* CreateKMatrixEnthalpy(void);
     170                ElementMatrix* CreateKMatrixEnthalpyVolume(void);
     171                ElementMatrix* CreateKMatrixEnthalpyShelf(void);
    169172                ElementMatrix* CreateKMatrixMelting(void);
    170173                ElementMatrix* CreateKMatrixPrognostic(void);
    … …  
    194197                ElementVector* CreatePVectorDiagnosticStokesViscous(void);
    195198                ElementVector* CreatePVectorDiagnosticStokesShelf(void);
     199                ElementVector* CreatePVectorEnthalpy(void);
     200                ElementVector* CreatePVectorEnthalpyVolume(void);
     201                ElementVector* CreatePVectorEnthalpyShelf(void);
     202                ElementVector* CreatePVectorEnthalpySheet(void);
    196203                ElementVector* CreatePVectorAdjointStokes(void);
    197204                ElementVector* CreatePVectorDiagnosticVert(void);

  • issm/trunk/src/c/solutions/CorePointerFromSolutionEnum.cpp

    r8429 r8483  
    3434                        solutioncore=&thermal_core;
    3535                        break;
     36                case EnthalpySolutionEnum:
     37                        solutioncore=&enthalpy_core;
     38                        break;
    3639                case PrognosticSolutionEnum:
    3740                        solutioncore=&prognostic_core;

  • issm/trunk/src/c/solutions/SolutionConfiguration.cpp

    r8429 r8483  
    5858                        analyses[0]=ThermalAnalysisEnum;
    5959                        analyses[1]=MeltingAnalysisEnum;
     60                        break;
     61               
     62                case EnthalpySolutionEnum:
     63                        numanalyses=1;
     64                        analyses=(int*)xmalloc(numanalyses*sizeof(int));
     65                        analyses[0]=EnthalpyAnalysisEnum;
    6066                        break;
    6167               

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

    r8429 r8483  
    2121void thermal_core(FemModel* femmodel);
    2222void thermal_core_step(FemModel* femmodel,int step, double time);
     23void enthalpy_core(FemModel* femmodel);
    2324void surfaceslope_core(FemModel* femmodel);
    2425void bedslope_core(FemModel* femmodel);
Note: See TracChangeset for help on using the changeset viewer.