Changeset 1676


Ignore:
Timestamp:
08/14/09 09:13:47 (16 years ago)
Author:
Mathieu Morlighem
Message:

Added artifficial diffusivity in thermal model

Location:
issm/trunk/src
Files:
4 edited

Legend:

Unmodified
Added
Removed

  • issm/trunk/src/c/ModelProcessorx/Thermal/CreateElementsNodesAndMaterialsThermal.cpp

    r1651 r1676  
    198198                penta_epsvel=model->epsvel;
    199199                penta_onwater=(bool)*(model->elementonwater+i);
     200                penta_artdiff=model->artificial_diffusivity;
    200201
    201202                /*We need the field collapse for transient, so that we can use compute B with the average temperature*/

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

    r1480 r1676  
    33893389        double  gauss_weight,area_gauss_weight,vert_gauss_weight;
    33903390        double  gauss_coord[4];
     3391        double  gauss_l1l2l3[3];
    33913392
    33923393        int area_order=5;
    … …  
    33953396        int     dofs[3]={0,1,2};
    33963397        double  dt;
     3398        double  K[2][2]={0.0};
    33973399
    33983400        double  vxvyvz_list[numgrids][3];
    … …  
    34113413        double     Ke_gaussian_conduct[numdof][numdof];
    34123414        double     Ke_gaussian_advec[numdof][numdof];
     3415        double     Ke_gaussian_artdiff[numdof][numdof];
    34133416        double     Ke_gaussian_transient[numdof][numdof];
    34143417        double     B[3][numdof];
    … …  
    34163419        double     B_conduct[3][numdof];
    34173420        double     B_advec[3][numdof];
     3421        double     B_artdiff[2][numdof];
    34183422        double     Bprime_advec[3][numdof];
    34193423        double     L[numdof];
    … …  
    34253429        double     tBD_conduct[3][numdof];
    34263430        double     tBD_advec[3][numdof];
     3431        double     tBD_artdiff[3][numdof];
    34273432        double     tLD[numdof];
    34283433
    … …  
    34663471                vz_list[i]=vxvyvz_list[i][2];
    34673472        }
     3473
    34683474
    34693475        /* Get gaussian points and weights. Penta is an extrusion of a Tria, we therefore
    … …  
    35563562                                }
    35573563                        }
     3564
     3565                        /*Artifficial diffusivity*/
     3566                        if(artdiff){
     3567                                /*Build K: */
     3568                                D_scalar=gauss_weight*Jdet/(pow(u,2)+pow(v,2)+epsvel);
     3569                                if(sub_analysis_type!=SteadyAnalysisEnum()){
     3570                                        D_scalar=D_scalar*dt;
     3571                                }
     3572                                K[0][0]=D_scalar*pow(u,2);       K[0][1]=D_scalar*fabs(u)*fabs(v);
     3573                                K[1][0]=D_scalar*fabs(u)*fabs(v);K[1][1]=D_scalar*pow(v,2);
     3574
     3575                                /*Get B_artdiff: */
     3576                                GetB_artdiff(&B_artdiff[0][0],&xyz_list[0][0],gauss_coord);
     3577
     3578                                /*  Do the triple product B'*K*B: */
     3579                                MatrixMultiply(&B_artdiff[0][0],2,numdof,1,&K[0][0],2,2,0,&tBD_artdiff[0][0],0);
     3580                                MatrixMultiply(&tBD_artdiff[0][0],numdof,2,0,&B_artdiff[0][0],2,numdof,0,&Ke_gaussian_artdiff[0][0],0);
     3581                        }
     3582                        else{
     3583                                for(i=0;i<numdof;i++){
     3584                                        for(j=0;j<numdof;j++){
     3585                                                Ke_gaussian_artdiff[i][j]=0;
     3586                                        }
     3587                                }
     3588                        }
    35583589                                       
    35593590                        /*Add Ke_gaussian to pKe: */
    35603591                        for(i=0;i<numdof;i++){
    35613592                                for(j=0;j<numdof;j++){
    3562                                         K_terms[i][j]+=Ke_gaussian_conduct[i][j]+Ke_gaussian_advec[i][j]+Ke_gaussian_transient[i][j];
     3593                                        K_terms[i][j]+=Ke_gaussian_conduct[i][j]+Ke_gaussian_advec[i][j]+Ke_gaussian_transient[i][j]+Ke_gaussian_artdiff[i][j];
    35633594                                }
    35643595                        }
    … …  
    36223653
    36233654#undef __FUNCT__
    3624 #define __FUNCT__ "Penta::GetB_advec"
    3625 void Penta::GetB_advec(double* B_advec, double* xyz_list, double* gauss_coord){
     3655#define __FUNCT__ "Penta::GetB_artdiff"
     3656void Penta::GetB_artdiff(double* B_artdiff, double* xyz_list, double* gauss_coord){
    36263657
    36273658        /*Compute B  matrix. B=[B1 B2 B3 B4 B5 B6] where Bi is of size 5*DOFPERGRID.
    36283659         * For grid i, Bi' can be expressed in the basic coordinate system
    36293660         * by:
    3630          *       Bi_conduct_basic=[ h ]
     3661         *       Bi_artdiff_basic=[ dh/dx ]
     3662         *                       [ dh/dy ]
     3663         * where h is the interpolation function for grid i.
     3664         *
     3665         * We assume B has been allocated already, of size: 2x(DOFPERGRID*numgrids)
     3666         */
     3667
     3668        int i;
     3669        const int calculationdof=3;
     3670        const int numgrids=6;
     3671        int DOFPERGRID=1;
     3672
     3673        /*Same thing in the basic coordinate system: */
     3674        double dh1dh6_basic[calculationdof][numgrids];
     3675
     3676        /*Get dh1dh2dh3 in basic coordinates system : */
     3677        GetNodalFunctionsDerivativesBasic(&dh1dh6_basic[0][0],xyz_list,gauss_coord);
     3678
     3679        /*Build B': */
     3680        for (i=0;i<numgrids;i++){
     3681                *(B_artdiff+DOFPERGRID*numgrids*0+DOFPERGRID*i)=dh1dh6_basic[0][i];
     3682                *(B_artdiff+DOFPERGRID*numgrids*1+DOFPERGRID*i)=dh1dh6_basic[1][i];
     3683        }
     3684}
     3685
     3686#undef __FUNCT__
     3687#define __FUNCT__ "Penta::GetB_advec"
     3688void Penta::GetB_advec(double* B_advec, double* xyz_list, double* gauss_coord){
     3689
     3690        /*Compute B  matrix. B=[B1 B2 B3 B4 B5 B6] where Bi is of size 5*DOFPERGRID.
     3691         * For grid i, Bi' can be expressed in the basic coordinate system
     3692         * by:
     3693         *       Bi_advec_basic =[ h ]
    36313694         *                       [ h ]
    36323695         *                       [ h ]

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

    r1188 r1676  
    137137                void  GetB_advec(double* B_advec, double* xyz_list, double* gauss_coord);
    138138                void  GetBprime_advec(double* Bprime_advec, double* xyz_list, double* gauss_coord);
     139                void  GetB_artdiff(double* B_artdiff, double* xyz_list, double* gauss_coord);
    139140                void  CreateKMatrixMelting(Mat Kgg,void* inputs,int analysis_type,int sub_analysis_type);
    140141                void  CreatePVectorThermal( Vec pg, void* vinputs,int analysis_type,int sub_analysis_type);

  • issm/trunk/src/m/classes/public/solveparallel.m

    r1268 r1676  
    3030        waitonlock([executionpath '/' md.name '.lock']);
    3131        %load results
     32        displaystring(md.debug,'loading results from cluster');
    3233        md=loadresultsfromcluster(md);
    3334end
Note: See TracChangeset for help on using the changeset viewer.