Changeset 24000

Timestamp:

06/07/19 13:46:04 (6 years ago)

Author:

tsantos

Message:

CHG: working on SUPG for ice thickness. Fixed a bug, added another part of the formulation

File:

• issm/trunk-jpl/src/c/analyses/MasstransportAnalysis.cpp (modified) (8 diffs)

issm/trunk-jpl/src/c/analyses/MasstransportAnalysis.cpp ¶

@@ -285 +285 @@
         IssmDouble Jdet,D_scalar,dt,h;
         IssmDouble vel,vx,vy,dvxdx,dvydy;
-        IssmDouble lambda_x,lambda_y,tau;
+        IssmDouble xi,tau;
         IssmDouble dvx[2],dvy[2];
         IssmDouble* xyz_list = NULL;
@@ -304 +304 @@
         ElementMatrix* Ke     = element->NewElementMatrix();
         IssmDouble*    basis  = xNew<IssmDouble>(numnodes);
+        IssmDouble*             dbasis = xNew<IssmDouble>(dim*numnodes);
         IssmDouble*    B      = xNew<IssmDouble>(dim*numnodes);
         IssmDouble*    Bprime = xNew<IssmDouble>(dim*numnodes);
@@ -367 +368 @@
                 switch(stabilization){
                         case 0:
-                                /*Nothing to be onde*/
+                                /*Nothing to be done*/
                                 break;
                         case 1:
@@ -393 +394 @@
                                 /*SUPG*/
                                 if(dim!=2) _error_("Stabilization "<<stabilization<<" not supported yet for dim != 2");
+                                element->NodalFunctionsDerivatives(dbasis,xyz_list,gauss);
+                                vxaverage_input->GetInputAverage(&vx);
+                                vyaverage_input->GetInputAverage(&vy);
                                 vel=sqrt(vx*vx+vy*vy)+1.e-8;
-                                tau=0.25;
-                                lambda_x=tau*h*vx/vel;
-                                lambda_y=tau*h*vy/vel;
+                                //xi=0.3130;
+                                xi=1;
+                                tau=xi*h/(2*vel);
                                 break;
                         default:
@@ -415 +419 @@
                                                 &Ke->values[0],1);
                 }
-                if(stabilization==5){
+                if(stabilization==5){/*{{{*/
                          /*Mass matrix - part 2*/
-                        D_scalar=gauss->weight*Jdet;
-                        D[0*dim+0]=D_scalar*lambda_x;
-                        D[1*dim+1]=D_scalar*lambda_y;
-                        TripleMultiply(Bprime,dim,numnodes,1,
-                                                                D,dim,dim,0,
-                                                                B,dim,numnodes,0,
-                                                                &Ke->values[0],1);
+                        for(int i=0;i<numnodes;i++){
+                                for(int j=0;j<numnodes;j++){
+                                        Ke->values[i*numnodes+j]+=gauss->weight*Jdet*tau*basis[j]*(vx*dbasis[0*numnodes+i]+vy*dbasis[1*numnodes+i]);    
+                                }
+                        }
+                        /*Mass matrix - part 3*/
+                        for(int i=0;i<numnodes;i++){
+                                for(int j=0;j<numnodes;j++){
+                                        Ke->values[i*numnodes+j]+=gauss->weight*Jdet*tau*basis[j]*(basis[i]*dvxdx+basis[i]*dvydy);      
+                                }
+                        }
                         
                         /*Advection matrix - part 2, A*/
-                        D_scalar=dt*gauss->weight*Jdet;
-                        D[0*dim+0]=D_scalar*dvxdx*lambda_x;
-                        D[1*dim+1]=D_scalar*dvydy*lambda_y;
-                        TripleMultiply(Bprime,dim,numnodes,1,
-                                                D,dim,dim,0,
-                                                B,dim,numnodes,0,
-                                                &Ke->values[0],1);
-                        
+                        for(int i=0;i<numnodes;i++){
+            for(int j=0;j<numnodes;j++){
+               Ke->values[i*numnodes+j]+=dt*gauss->weight*Jdet*tau*(vx*dbasis[0*numnodes+j]+vy*dbasis[1*numnodes+j])*(vx*dbasis[0*numnodes+i]+vy*dbasis[1*numnodes+i]);
+            }
+         }
+                        /*Advection matrix - part 3, A*/
+                        for(int i=0;i<numnodes;i++){
+            for(int j=0;j<numnodes;j++){
+                                        Ke->values[i*numnodes+j]+=dt*gauss->weight*Jdet*tau*(vx*dbasis[0*numnodes+j]+vy*dbasis[1*numnodes+j])*(basis[i]*dvxdx+basis[i]*dvydy);
+                                }
+         }
+
                         /*Advection matrix - part 2, B*/
-                        D[0*dim+0]=D_scalar*vx*lambda_x;
-                        D[1*dim+1]=D_scalar*vy*lambda_y;
-                        TripleMultiply(Bprime,dim,numnodes,1,
-                                                D,dim,dim,0,
-                                                Bprime,dim,numnodes,0,
-                                                &Ke->values[0],1);
-                
-                }
+                        for(int i=0;i<numnodes;i++){
+            for(int j=0;j<numnodes;j++){
+                                        Ke->values[i*numnodes+j]+=dt*gauss->weight*Jdet*tau*(basis[j]*dvxdx+basis[j]*dvydy)*(vx*dbasis[0*numnodes+i]+vy*dbasis[1*numnodes+i]);
+                                }
+         }
+                        /*Advection matrix - part 3, B*/
+                        for(int i=0;i<numnodes;i++){
+            for(int j=0;j<numnodes;j++){
+                                        Ke->values[i*numnodes+j]+=dt*gauss->weight*Jdet*tau*(basis[j]*dvxdx+basis[j]*dvydy)*(basis[i]*dvxdx+basis[i]*dvydy);    
+                                }
+                        }
+                }/*}}}*/
         }
 
@@ -448 +464 @@
         xDelete<IssmDouble>(xyz_list);
         xDelete<IssmDouble>(basis);
+        xDelete<IssmDouble>(dbasis);
         xDelete<IssmDouble>(B);
         xDelete<IssmDouble>(Bprime);
@@ -559 +576 @@
         IssmDouble  Jdet,dt;
         IssmDouble  ms,mb,gmb,fmb,thickness;
-        IssmDouble  vx,vy,vel,lambda_x,lambda_y,h,tau;
+        IssmDouble  vx,vy,vel,dvxdx,dvydy,xi,h,tau;
+        IssmDouble  dvx[2],dvy[2];
         IssmDouble  gllevelset,phi=1.;
         IssmDouble* xyz_list = NULL;
@@ -640 +658 @@
         
                 if(stabilization==5){ //SUPG
-                        /*Prepare coefficients*/
-                        vxaverage_input->GetInputValue(&vx,gauss);
-                        vxaverage_input->GetInputValue(&vy,gauss);
-                        vel = sqrt(vx*vx+vy*vy)+1.e-8;
-                        tau=0.25;
-                        lambda_x = tau*h*vx/vel;
-                        lambda_y = tau*h*vy/vel;
-        
-                        /*Get nodal derivatives*/
                         element->NodalFunctionsDerivatives(dbasis,xyz_list,gauss);
+                        vxaverage_input->GetInputAverage(&vx);
+                        vyaverage_input->GetInputAverage(&vy);
+                        vxaverage_input->GetInputDerivativeValue(&dvx[0],xyz_list,gauss);
+                        vyaverage_input->GetInputDerivativeValue(&dvy[0],xyz_list,gauss);
+                        vel=sqrt(vx*vx+vy*vy)+1.e-8;
+                        dvxdx=dvx[0];
+                        dvydy=dvy[1];
+                        //xi=0.3130;
+                        xi=1;
+                        tau=xi*h/(2*vel);
                         
                         /*Force vector - part 2*/
                         for(int i=0;i<numnodes;i++){
-                                pe->values[i]+=Jdet*gauss->weight*(thickness+dt*(ms-mb))*(lambda_x*dbasis[numnodes*0+i]+lambda_y*dbasis[numnodes*1+i]);
+                                pe->values[i]+=Jdet*gauss->weight*(thickness+dt*(ms-mb))*(tau*vx*dbasis[0*numnodes+i]+tau*vy*dbasis[1*numnodes+i]);
+                        }
+                        /*Force vector - part 3*/
+                        for(int i=0;i<numnodes;i++){
+                                pe->values[i]+=Jdet*gauss->weight*(thickness+dt*(ms-mb))*(tau*basis[i]*dvxdx+tau*basis[i]*dvydy);
                         }
                 }