Changeset 10628


Ignore:
Timestamp:
11/14/11 10:59:49 (13 years ago)
Author:
Mathieu Morlighem
Message:

Better artificial diffusivity in thermal model

Location:
issm/trunk/src/c/objects/Elements
Files:
4 edited

Legend:

Unmodified
Added
Removed

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

    r10597 r10628  
    822822        /*return PentaRef field*/
    823823        return this->element_type;
     824}
     825/*}}}*/
     826/*FUNCTION Penta::GetElementSizes{{{1*/
     827void Penta::GetElementSizes(double* hx,double* hy,double* hz){
     828
     829        double xyz_list[NUMVERTICES][3];
     830        double xmin,ymin,zmin;
     831        double xmax,ymax,zmax;
     832
     833        /*Get xyz list: */
     834        GetVerticesCoordinates(&xyz_list[0][0],nodes,NUMVERTICES);
     835        xmin=xyz_list[0][0]; xmax=xyz_list[0][0];
     836        ymin=xyz_list[0][1]; ymax=xyz_list[0][1];
     837        zmin=xyz_list[0][2]; zmax=xyz_list[0][2];
     838
     839        for(int i=1;i<NUMVERTICES;i++){
     840                if(xyz_list[i][0]<xmin) xmin=xyz_list[i][0];
     841                if(xyz_list[i][0]>xmax) xmax=xyz_list[i][0];
     842                if(xyz_list[i][1]<ymin) ymin=xyz_list[i][1];
     843                if(xyz_list[i][1]>ymax) ymax=xyz_list[i][1];
     844                if(xyz_list[i][2]<zmin) zmin=xyz_list[i][2];
     845                if(xyz_list[i][2]>zmax) zmax=xyz_list[i][2];
     846        }
     847
     848        *hx=xmax-xmin;
     849        *hy=ymax-ymin;
     850        *hz=zmax-zmin;
    824851}
    825852/*}}}*/
    … …  
    31803207                        K[0][0]=D_scalar_stab*pow(u,2);       K[0][1]=D_scalar_stab*fabs(u)*fabs(v);
    31813208                        K[1][0]=D_scalar_stab*fabs(u)*fabs(v);K[1][1]=D_scalar_stab*pow(v,2);
    3182 
    3183                         GetBArtdiff(&B_stab[0][0],&xyz_list[0][0],gauss);
     3209                        _error_("TO BE RECODED");
     3210
     3211                        //GetBArtdiff(&B_stab[0][0],&xyz_list[0][0],gauss);
    31843212
    31853213                        TripleMultiply(&B_stab[0][0],2,numdof,1,
    … …  
    33043332        int        stabilization;
    33053333        int        i,j,ig,found=0;
    3306         double     Jdet,u,v,w,um,vm,wm;
    3307         double     epsvel=2.220446049250313e-16;
     3334        double     Jdet,u,v,w,um,vm,wm,vel;
     3335        double     h,hx,hy,hz,vx,vy,vz;
    33083336        double     gravity,rho_ice,rho_water;
    33093337        double     heatcapacity,thermalconductivity,dt;
    … …  
    33213349        double     D_scalar_trans,D_scalar_stab;
    33223350        double     D[3][3];
    3323         double     K[2][2]={0.0};
     3351        double     K[3][3]={0.0};
    33243352        Tria*      tria=NULL;
    33253353        GaussPenta *gauss=NULL;
    … …  
    33743402                GetBprimeAdvec(&Bprime_advec[0][0],&xyz_list[0][0],gauss);
    33753403
    3376                 vx_input->GetInputValue(&u, gauss);
    3377                 vy_input->GetInputValue(&v, gauss);
    3378                 vz_input->GetInputValue(&w, gauss);
    3379                 vxm_input->GetInputValue(&um,gauss);
    3380                 vym_input->GetInputValue(&vm,gauss);
    3381                 vzm_input->GetInputValue(&wm,gauss);
     3404                vx_input->GetInputValue(&u, gauss); vxm_input->GetInputValue(&um,gauss); vx=u-um;
     3405                vy_input->GetInputValue(&v, gauss); vym_input->GetInputValue(&vm,gauss); vy=v-vm;
     3406                vz_input->GetInputValue(&w, gauss); vzm_input->GetInputValue(&wm,gauss); vz=w-wm;
    33823407
    33833408                D_scalar_advec=gauss->weight*Jdet;
    33843409                if(dt) D_scalar_advec=D_scalar_advec*dt;
    33853410
    3386                 D[0][0]=D_scalar_advec*(u-um);D[0][1]=0;                    D[0][2]=0;
    3387                 D[1][0]=0;                    D[1][1]=D_scalar_advec*(v-vm);D[1][2]=0;
    3388                 D[2][0]=0;                    D[2][1]=0;                    D[2][2]=D_scalar_advec*(w-wm);
     3411                D[0][0]=D_scalar_advec*vx;    D[0][1]=0;                    D[0][2]=0;
     3412                D[1][0]=0;                    D[1][1]=D_scalar_advec*vy;    D[1][2]=0;
     3413                D[2][0]=0;                    D[2][1]=0;                    D[2][2]=D_scalar_advec*vz;
    33893414
    33903415                TripleMultiply(&B_advec[0][0],3,numdof,1,
    … …  
    34103435                if(stabilization==1){
    34113436                        /*Build K: */
    3412                         D_scalar_stab=gauss->weight*Jdet/(pow(u-um,2)+pow(v-vm,2)+epsvel);
     3437                        GetElementSizes(&hx,&hy,&hz);
     3438                        vel=sqrt(pow(vx,2.)+pow(vy,2.)+pow(vz,2.))+1.e-14;
     3439                        h=sqrt( pow(hx*vx/vel,2.) + pow(hy*vy/vel,2.) + pow(hz*vz/vel,2.));
     3440
     3441                        K[0][0]=h/(2*vel)*vx*vx;  K[0][1]=h/(2*vel)*vx*vy; K[0][2]=h/(2*vel)*vx*vz;
     3442                        K[1][0]=h/(2*vel)*vy*vx;  K[1][1]=h/(2*vel)*vy*vy; K[1][2]=h/(2*vel)*vy*vz;
     3443                        K[2][0]=h/(2*vel)*vz*vx;  K[2][1]=h/(2*vel)*vz*vy; K[2][2]=h/(2*vel)*vz*vz;
     3444
     3445                        D_scalar_stab=gauss->weight*Jdet;
    34133446                        if(dt) D_scalar_stab=D_scalar_stab*dt;
    3414                         K[0][0]=D_scalar_stab*pow(u,2);       K[0][1]=D_scalar_stab*fabs(u)*fabs(v);
    3415                         K[1][0]=D_scalar_stab*fabs(u)*fabs(v);K[1][1]=D_scalar_stab*pow(v,2);
    3416 
    3417                         GetBArtdiff(&B_stab[0][0],&xyz_list[0][0],gauss);
    3418 
    3419                         TripleMultiply(&B_stab[0][0],2,numdof,1,
    3420                                                 &K[0][0],2,2,0,
    3421                                                 &B_stab[0][0],2,numdof,0,
     3447                        for(i=0;i<3;i++) for(j=0;j<3;j++) K[i][j] = D_scalar_stab*K[i][j];
     3448
     3449                        GetBprimeAdvec(&Bprime_advec[0][0],&xyz_list[0][0],gauss);
     3450
     3451                        TripleMultiply(&Bprime_advec[0][0],3,numdof,1,
     3452                                                &K[0][0],3,3,0,
     3453                                                &Bprime_advec[0][0],3,numdof,0,
    34223454                                                &Ke->values[0],1);
    34233455                }
    34243456                else if(stabilization==2){
    3425 
    34263457                        GetNodalFunctionsP1Derivatives(&dbasis[0][0],&xyz_list[0][0], gauss);
    3427 
    34283458                        tau_parameter=GetStabilizationParameter(u-um,v-vm,w-wm,diameter,rho_ice,heatcapacity,thermalconductivity);
    34293459

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

    r10571 r10628  
    170170                void    GetSidList(int* sidlist);
    171171                int     GetElementType(void);
     172                void    GetElementSizes(double* hx,double* hy,double* hz);
    172173                void    GetInputListOnVertices(double* pvalue,int enumtype);
    173174                void    GetInputListOnVertices(double* pvalue,int enumtype,double defaultvalue);

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

    r10514 r10628  
    384384        }
    385385
    386 }
    387 /*}}}*/
    388 /*FUNCTION PentaRef::GetBArtdiff {{{1*/
    389 void PentaRef::GetBArtdiff(double* B_stab, double* xyz_list, GaussPenta* gauss){
    390         /*Compute B  matrix. B=[B1 B2 B3 B4 B5 B6] where Bi is of size 5*NDOF1.
    391          * For node i, Bi' can be expressed in the actual coordinate system
    392          * by:
    393          *       Bi_stab=[ dh/dx ]
    394          *                 [ dh/dy ]
    395          * where h is the interpolation function for node i.
    396          *
    397          * We assume B has been allocated already, of size: 2x(NDOF1*NUMNODESP1)
    398          */
    399 
    400         /*Same thing in the actual coordinate system: */
    401         double dbasis[3][NUMNODESP1];
    402 
    403         /*Get dbasis in actual coordinates system : */
    404         GetNodalFunctionsP1Derivatives(&dbasis[0][0],xyz_list,gauss);
    405 
    406         /*Build B': */
    407         for (int i=0;i<NUMNODESP1;i++){
    408                 *(B_stab+NDOF1*NUMNODESP1*0+NDOF1*i)=dbasis[0][i];
    409                 *(B_stab+NDOF1*NUMNODESP1*1+NDOF1*i)=dbasis[1][i];
    410         }
    411386}
    412387/*}}}*/

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

    r10514 r10628  
    4545                void GetBprimeVert(double* B, double* xyz_list, GaussPenta* gauss);
    4646                void GetBAdvec(double* B_advec, double* xyz_list, GaussPenta* gauss);
    47                 void GetBArtdiff(double* B_stab, double* xyz_list, GaussPenta* gauss);
    4847                void GetBConduct(double* B_conduct, double* xyz_list, GaussPenta* gauss);
    4948                void GetBVert(double* B, double* xyz_list, GaussPenta* gauss);
Note: See TracChangeset for help on using the changeset viewer.