Changeset 23985

Timestamp:

06/04/19 16:51:10 (6 years ago)

Author:

tsantos

Message:

NEW: added Streamline Updwind Petrov-Galerkin stabilization for Ice Thickness advection

Location:

issm/trunk-jpl/src

Files:

• c/analyses/MasstransportAnalysis.cpp (modified) (10 diffs)
• m/classes/masstransport.m (modified) (2 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 dvx[2],dvy[2];
         IssmDouble* xyz_list = NULL;
@@ -389 +390 @@
                                 }
                                 break;
+                        case 5:
+                                /*SUPG*/
+                                if(dim!=2) _error_("Stabilization "<<stabilization<<" not supported yet for dim != 2");
+                                vel=sqrt(vx*vx+vy*vy)+1.e-8;
+                                tau=0.25;
+                                lambda_x=tau*h*vx/vel;
+                                lambda_y=tau*h*vy/vel;
+                                break;
                         default:
                                 _error_("Stabilization "<<stabilization<<" not supported yet");
@@ -405 +414 @@
                                                 Bprime,dim,numnodes,0,
                                                 &Ke->values[0],1);
+                }
+                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);
+                        
+                        /*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);
+                        
+                        /*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);
+                
                 }
         }
@@ -516 +553 @@
 
         /*Intermediaries */
+        int                     stabilization,dim,domaintype;
         int         melt_style,point1;
         bool        mainlyfloating;
@@ -521 +559 @@
         IssmDouble  Jdet,dt;
         IssmDouble  ms,mb,gmb,fmb,thickness;
+        IssmDouble  vx,vy,vel,lambda_x,lambda_y,h,tau;
         IssmDouble  gllevelset,phi=1.;
         IssmDouble* xyz_list = NULL;
         Gauss*      gauss     = NULL;
 
+        /*Get problem dimension*/
+        element->FindParam(&domaintype,DomainTypeEnum);
+        switch(domaintype){
+                case Domain2DverticalEnum:   dim = 1; break;
+                case Domain2DhorizontalEnum: dim = 2; break;
+                case Domain3DEnum:           dim = 2; break;
+                default: _error_("mesh "<<EnumToStringx(domaintype)<<" not supported yet");
+        }
+          
         /*Fetch number of nodes and dof for this finite element*/
         int numnodes = element->GetNumberOfNodes();
@@ -531 +579 @@
         ElementVector* pe    = element->NewElementVector();
         IssmDouble*    basis = xNew<IssmDouble>(numnodes);
+        IssmDouble*    dbasis= xNew<IssmDouble>(dim*numnodes);
 
         /*Retrieve all inputs and parameters*/
@@ -536 +585 @@
         element->FindParam(&melt_style,GroundinglineMeltInterpolationEnum);
         element->FindParam(&dt,TimesteppingTimeStepEnum);
+        element->FindParam(&stabilization,MasstransportStabilizationEnum);
         Input* gmb_input        = element->GetInput(BasalforcingsGroundediceMeltingRateEnum);  _assert_(gmb_input);
         Input* fmb_input        = element->GetInput(BasalforcingsFloatingiceMeltingRateEnum);  _assert_(fmb_input);
@@ -541 +591 @@
         Input* ms_input         = element->GetInput(SmbMassBalanceEnum);                       _assert_(ms_input);
         Input* thickness_input  = element->GetInput(ThicknessEnum);                            _assert_(thickness_input);
-
+        Input* vxaverage_input  = element->GetInput(VxAverageEnum);                                                                             _assert_(vxaverage_input);
+        Input* vyaverage_input  = element->GetInput(VyAverageEnum);                                                                             _assert_(vyaverage_input);
+
+        h=element->CharacteristicLength();
+        
         /*Recover portion of element that is grounded*/
         phi=element->GetGroundedPortion(xyz_list);
@@ -584 +638 @@
 
                 for(int i=0;i<numnodes;i++) pe->values[i]+=Jdet*gauss->weight*(thickness+dt*(ms-mb))*basis[i];
+        
+                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);
+                        
+                        /*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]);
+                        }
+                }
+        
         }
 
@@ -589 +662 @@
         xDelete<IssmDouble>(xyz_list);
         xDelete<IssmDouble>(basis);
+        xDelete<IssmDouble>(dbasis);
         delete gauss;
         return pe;

issm/trunk-jpl/src/m/classes/masstransport.m

@@ -92 +92 @@
                         md = checkfield(md,'fieldname','masstransport.isfreesurface','values',[0 1]);
                         md = checkfield(md,'fieldname','masstransport.hydrostatic_adjustment','values',{'Absolute' 'Incremental'});
-                        md = checkfield(md,'fieldname','masstransport.stabilization','values',[0 1 2 3 4]);
+                        md = checkfield(md,'fieldname','masstransport.stabilization','values',[0 1 2 3 4 5]);
                         md = checkfield(md,'fieldname','masstransport.min_thickness','>',0);
                         md = checkfield(md,'fieldname','masstransport.requested_outputs','stringrow',1);
@@ -103 +103 @@
                         fielddisplay(self,'min_thickness','minimum ice thickness allowed [m]');
                         fielddisplay(self,'hydrostatic_adjustment','adjustment of ice shelves surface and bed elevations: ''Incremental'' or ''Absolute'' ');
-                        fielddisplay(self,'stabilization','0: no, 1: artificial_diffusivity, 2: streamline upwinding, 3: discontinuous Galerkin, 4: Flux Correction Transport');
+                        fielddisplay(self,'stabilization','0: no, 1: artificial_diffusivity, 2: streamline upwinding, 3: discontinuous Galerkin, 4: Flux Correction Transport, 5: streamline upwind Petrov-Galerkin (SUPG)');
 
                         disp(sprintf('\n      %s','Penalty options:'));