Changeset 22534


Ignore:
Timestamp:
03/12/18 16:14:34 (7 years ago)
Author:
Mathieu Morlighem
Message:

CHG: got rid of B and Bprime for FS, in hope to increase model performance

File:
1 edited

Legend:

Unmodified
Added
Removed

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

    r22513 r22534  
    31353135
    31363136        /*Intermediaries*/
    3137         int         i,dim,epssize;
     3137        int         i,dim;
    31383138        IssmDouble  viscosity,FSreconditioning,Jdet;
    31393139        IssmDouble *xyz_list = NULL;
    … …  
    31413141        /*Get problem dimension*/
    31423142        element->FindParam(&dim,DomainDimensionEnum);
    3143         if(dim==2) epssize = 3;
    3144         else       epssize = 6;
    31453143
    31463144        /*Fetch number of nodes and dof for this finite element*/
    … …  
    31483146        int pnumnodes = element->NumberofNodesPressure();
    31493147        int numdof    = vnumnodes*dim + pnumnodes;
    3150         int bsize     = epssize + 2;
    31513148
    31523149        /*Prepare coordinate system list*/
    … …  
    31573154
    31583155        /*Initialize Element matrix and vectors*/
    3159         ElementMatrix* Ke     = element->NewElementMatrix(FSvelocityEnum);
    3160         IssmDouble*    B      = xNew<IssmDouble>(bsize*numdof);
    3161         IssmDouble*    Bprime = xNew<IssmDouble>(bsize*numdof);
    3162         IssmDouble*    D      = xNewZeroInit<IssmDouble>(bsize*bsize);
     3156        ElementMatrix* Ke   = element->NewElementMatrix(FSvelocityEnum);
     3157        IssmDouble* vdbasis = xNew<IssmDouble>(dim*vnumnodes);
     3158        IssmDouble* pbasis  = xNew<IssmDouble>(pnumnodes);
    31633159
    31643160        /*Retrieve all inputs and parameters*/
    … …  
    31673163        Input* vx_input=element->GetInput(VxEnum);     _assert_(vx_input);
    31683164        Input* vy_input=element->GetInput(VyEnum);     _assert_(vy_input);
    3169         Input* vz_input;
     3165        Input* vz_input = NULL;
    31703166        if(dim==3){vz_input=element->GetInput(VzEnum); _assert_(vz_input);}
    31713167
    … …  
    31763172
    31773173                element->JacobianDeterminant(&Jdet,xyz_list,gauss);
    3178                 this->GetBFS(B,element,dim,xyz_list,gauss);
    3179                 this->GetBFSprime(Bprime,element,dim,xyz_list,gauss);
    3180 
     3174                element->NodalFunctionsDerivativesVelocity(vdbasis,xyz_list,gauss);
     3175                element->NodalFunctionsPressure(pbasis,gauss);
    31813176                element->material->ViscosityFS(&viscosity,dim,xyz_list,gauss,vx_input,vy_input,vz_input);
    3182                 for(i=0;i<epssize;i++)     D[i*bsize+i] = + 2.*viscosity*gauss->weight*Jdet;
    3183                 for(i=epssize;i<bsize;i++) D[i*bsize+i] = - FSreconditioning*gauss->weight*Jdet;
    3184 
    3185                 TripleMultiply(B,bsize,numdof,1,
    3186                                         D,bsize,bsize,0,
    3187                                         Bprime,bsize,numdof,0,
    3188                                         &Ke->values[0],1);
     3177
     3178                if(dim==3){
     3179                        /*Stress balance*/
     3180                        for(int i=0;i<vnumnodes;i++){
     3181                                for(int j=0;j<vnumnodes;j++){
     3182                                        Ke->values[(3*i+0)*numdof+3*j+0] += gauss->weight*Jdet*viscosity*(2.*vdbasis[0*vnumnodes+j]*vdbasis[0*vnumnodes+i] + vdbasis[1*vnumnodes+j]*vdbasis[1*vnumnodes+i] + vdbasis[2*vnumnodes+j]*vdbasis[2*vnumnodes+i]);
     3183                                        Ke->values[(3*i+0)*numdof+3*j+1] += gauss->weight*Jdet*viscosity*(vdbasis[0*vnumnodes+j]*vdbasis[1*vnumnodes+i]);
     3184                                        Ke->values[(3*i+0)*numdof+3*j+2] += gauss->weight*Jdet*viscosity*(vdbasis[0*vnumnodes+j]*vdbasis[2*vnumnodes+i]);
     3185                                        Ke->values[(3*i+1)*numdof+3*j+0] += gauss->weight*Jdet*viscosity*(vdbasis[1*vnumnodes+j]*vdbasis[0*vnumnodes+i]);
     3186                                        Ke->values[(3*i+1)*numdof+3*j+1] += gauss->weight*Jdet*viscosity*(vdbasis[0*vnumnodes+j]*vdbasis[0*vnumnodes+i] + 2.*vdbasis[1*vnumnodes+j]*vdbasis[1*vnumnodes+i] + vdbasis[2*vnumnodes+j]*vdbasis[2*vnumnodes+i]);
     3187                                        Ke->values[(3*i+1)*numdof+3*j+2] += gauss->weight*Jdet*viscosity*(vdbasis[1*vnumnodes+j]*vdbasis[2*vnumnodes+i]);
     3188                                        Ke->values[(3*i+2)*numdof+3*j+0] += gauss->weight*Jdet*viscosity*(vdbasis[2*vnumnodes+j]*vdbasis[0*vnumnodes+i]);
     3189                                        Ke->values[(3*i+2)*numdof+3*j+1] += gauss->weight*Jdet*viscosity*(vdbasis[2*vnumnodes+j]*vdbasis[1*vnumnodes+i]);
     3190                                        Ke->values[(3*i+2)*numdof+3*j+2] += gauss->weight*Jdet*viscosity*(vdbasis[0*vnumnodes+j]*vdbasis[0*vnumnodes+i] + vdbasis[1*vnumnodes+j]*vdbasis[1*vnumnodes+i] + 2.*vdbasis[2*vnumnodes+j]*vdbasis[2*vnumnodes+i]);
     3191                                }
     3192                                for(int k=0;k<pnumnodes;k++){
     3193                                        Ke->values[(3*i+0)*numdof+3*vnumnodes+k] += gauss->weight*Jdet*FSreconditioning*(-pbasis[k]*vdbasis[0*vnumnodes+i]);
     3194                                        Ke->values[(3*i+1)*numdof+3*vnumnodes+k] += gauss->weight*Jdet*FSreconditioning*(-pbasis[k]*vdbasis[1*vnumnodes+i]);
     3195                                        Ke->values[(3*i+2)*numdof+3*vnumnodes+k] += gauss->weight*Jdet*FSreconditioning*(-pbasis[k]*vdbasis[2*vnumnodes+i]);
     3196                                }
     3197                        }
     3198                        /*Incompressibility*/
     3199                        for(int k=0;k<pnumnodes;k++){
     3200                                for(int j=0;j<vnumnodes;j++){
     3201                                        Ke->values[(3*vnumnodes+k)*numdof+3*j+0] += gauss->weight*Jdet*(-FSreconditioning*vdbasis[0*vnumnodes+j]*pbasis[k]);
     3202                                        Ke->values[(3*vnumnodes+k)*numdof+3*j+1] += gauss->weight*Jdet*(-FSreconditioning*vdbasis[1*vnumnodes+j]*pbasis[k]);
     3203                                        Ke->values[(3*vnumnodes+k)*numdof+3*j+2] += gauss->weight*Jdet*(-FSreconditioning*vdbasis[2*vnumnodes+j]*pbasis[k]);
     3204                                }
     3205                        }
     3206                }
     3207                else{
     3208                        /*Stress balance*/
     3209                        for(int i=0;i<vnumnodes;i++){
     3210                                for(int j=0;j<vnumnodes;j++){
     3211                                        Ke->values[(2*i+0)*numdof+2*j+0] += gauss->weight*Jdet*viscosity*(2.*vdbasis[0*vnumnodes+j]*vdbasis[0*vnumnodes+i] + vdbasis[1*vnumnodes+j]*vdbasis[1*vnumnodes+i]);
     3212                                        Ke->values[(2*i+0)*numdof+2*j+1] += gauss->weight*Jdet*viscosity*(vdbasis[0*vnumnodes+j]*vdbasis[1*vnumnodes+i]);
     3213                                        Ke->values[(2*i+1)*numdof+2*j+0] += gauss->weight*Jdet*viscosity*(vdbasis[1*vnumnodes+j]*vdbasis[0*vnumnodes+i]);
     3214                                        Ke->values[(2*i+1)*numdof+2*j+1] += gauss->weight*Jdet*viscosity*(vdbasis[0*vnumnodes+j]*vdbasis[0*vnumnodes+i] + 2.*vdbasis[1*vnumnodes+j]*vdbasis[1*vnumnodes+i]);
     3215                                }
     3216                                for(int k=0;k<pnumnodes;k++){
     3217                                        Ke->values[(2*i+0)*numdof+2*vnumnodes+k] += gauss->weight*Jdet*FSreconditioning*(-pbasis[k]*vdbasis[0*vnumnodes+i]);
     3218                                        Ke->values[(2*i+1)*numdof+2*vnumnodes+k] += gauss->weight*Jdet*FSreconditioning*(-pbasis[k]*vdbasis[1*vnumnodes+i]);
     3219                                }
     3220                        }
     3221                        /*Incompressibility*/
     3222                        for(int k=0;k<pnumnodes;k++){
     3223                                for(int j=0;j<vnumnodes;j++){
     3224                                        Ke->values[(2*vnumnodes+k)*numdof+2*j+0] += gauss->weight*Jdet*(-FSreconditioning*vdbasis[0*vnumnodes+j]*pbasis[k]);
     3225                                        Ke->values[(2*vnumnodes+k)*numdof+2*j+1] += gauss->weight*Jdet*(-FSreconditioning*vdbasis[1*vnumnodes+j]*pbasis[k]);
     3226                                }
     3227                        }
     3228                }
    31893229        }
    31903230
    … …  
    31953235        delete gauss;
    31963236        xDelete<IssmDouble>(xyz_list);
    3197         xDelete<IssmDouble>(D);
    3198         xDelete<IssmDouble>(Bprime);
    3199         xDelete<IssmDouble>(B);
     3237        xDelete<IssmDouble>(pbasis);
     3238        xDelete<IssmDouble>(vdbasis);
    32003239        xDelete<int>(cs_list);
    32013240        return Ke;
    … …  
    37063745        /*Initialize Element matrix and vectors*/
    37073746        ElementMatrix* Ke = element->NewElementMatrix(FSvelocityEnum);
    3708         IssmDouble*    B  = xNew<IssmDouble>(dim*numdof);
    3709         IssmDouble*    D  = xNewZeroInit<IssmDouble>(dim*dim);
     3747        IssmDouble* vbasis=xNew<IssmDouble>(vnumnodes);
    37103748
    37113749        /*Retrieve all inputs and parameters*/
    … …  
    37413779                }
    37423780
    3743                 this->GetBFSFriction(B,element,dim,xyz_list_base,gauss);
    37443781                element->JacobianDeterminantBase(&Jdet,xyz_list_base,gauss);
    3745                 for(int i=0;i<dim;i++) D[i*dim+i] = alpha2*gauss->weight*Jdet; //taub_x = -alpha2 v_x (same for y)
    3746 
    3747                 TripleMultiply(B,dim,numdof,1,
    3748                                         D,dim,dim,0,
    3749                                         B,dim,numdof,0,
    3750                                         &Ke->values[0],1);
     3782                element->NodalFunctionsVelocity(vbasis,gauss);
     3783
     3784                if(dim==3){
     3785                        /*Stress balance*/
     3786                        for(int i=0;i<vnumnodes;i++){
     3787                                for(int j=0;j<vnumnodes;j++){
     3788                                        Ke->values[(3*i+0)*numdof+3*j+0] += alpha2*gauss->weight*Jdet*vbasis[i]*vbasis[j];
     3789                                        Ke->values[(3*i+1)*numdof+3*j+1] += alpha2*gauss->weight*Jdet*vbasis[i]*vbasis[j];
     3790                                        Ke->values[(3*i+2)*numdof+3*j+2] += alpha2*gauss->weight*Jdet*vbasis[i]*vbasis[j];
     3791                                }
     3792                        }
     3793                }
     3794                else{
     3795                        /*Stress balance*/
     3796                        for(int i=0;i<vnumnodes;i++){
     3797                                for(int j=0;j<vnumnodes;j++){
     3798                                        Ke->values[(2*i+0)*numdof+2*j+0] += alpha2*gauss->weight*Jdet*vbasis[i]*vbasis[j];
     3799                                        Ke->values[(2*i+1)*numdof+2*j+1] += alpha2*gauss->weight*Jdet*vbasis[i]*vbasis[j];
     3800                                }
     3801                        }
     3802                }
    37513803        }
    37523804
    … …  
    37573809        delete friction;
    37583810        xDelete<IssmDouble>(xyz_list_base);
    3759         xDelete<IssmDouble>(B);
    3760         xDelete<IssmDouble>(D);
     3811        xDelete<IssmDouble>(vbasis);
    37613812        return Ke;
    37623813}/*}}}*/
Note: See TracChangeset for help on using the changeset viewer.