Changeset 26420

Timestamp:

09/01/21 18:23:05 (4 years ago)

Author:

Mathieu Morlighem

Message:

CHG: added new stabilization schemes for levelset transport

File:

• issm/trunk-jpl/src/c/analyses/LevelsetAnalysis.cpp (modified) (7 diffs)

issm/trunk-jpl/src/c/analyses/LevelsetAnalysis.cpp

@@ -215 +215 @@
         IssmDouble*    basis    = xNew<IssmDouble>(numnodes);
         IssmDouble*    dbasis   = xNew<IssmDouble>(2*numnodes);
-        IssmDouble*    Bprime = NULL;
-        if(stabilization==2){
-                Bprime   = xNew<IssmDouble>(dim*numnodes);
-        }
 
         /*Retrieve all inputs and parameters*/
@@ -224 +220 @@
         basalelement->FindParam(&dt,TimesteppingTimeStepEnum);
         basalelement->FindParam(&migrationmax,MigrationMaxEnum);
+
+        h = element->CharacteristicLength();
 
         Input* mf_vx_input        = NULL;
@@ -291 +289 @@
                 switch(stabilization){
                         case 0:
-                                // no stabilization, do nothing
+                                /*Nothing to be done*/
                                 break;
                         case 1:
@@ -309 +307 @@
                                   {
                                         /* Streamline Upwinding */
-                                        basalelement->ElementSizes(&hx,&hy,&hz);
-                                        h=sqrt( pow(hx*w[0]/vel,2) + pow(hy*w[1]/vel,2) );
+                                        mf_vx_input->GetInputAverage(&w[0]);
+                                        mf_vy_input->GetInputAverage(&w[1]);
+                                        vel=sqrt(w[0]*w[0]+w[1]*w[1])+1.e-8;
+                                        IssmDouble tau=h/(2*vel);
                                         for(int i=0;i<numnodes;i++){
                                                 for(int j=0;j<numnodes;j++){
-                                                        Ke->values[i*numnodes+j] += D_scalar*h/(2.*vel)*(
-                                                                                dbasis[0*numnodes+i] *(w[0]*w[0]*dbasis[0*numnodes+j] + w[0]*w[1]*dbasis[1*numnodes+j]) +
-                                                                                dbasis[1*numnodes+i] *(w[1]*w[0]*dbasis[0*numnodes+j] + w[1]*w[1]*dbasis[1*numnodes+j]) 
-                                                                                );
+                                                        Ke->values[i*numnodes+j]+=dt*gauss->weight*Jdet*tau*(
+                                                                                w[0]*dbasis[0*numnodes+i]+w[1]*dbasis[1*numnodes+i])*(w[0]*dbasis[0*numnodes+j]+w[1]*dbasis[1*numnodes+j]);
                                                 }
                                         }
                                   }
                                 break;
+                        case 5:{
+                                /*SUPG*/
+                                mf_vx_input->GetInputAverage(&w[0]);
+                                mf_vy_input->GetInputAverage(&w[1]);
+                                vel=sqrt(w[0]*w[0]+w[1]*w[1])+1.e-8;
+                                IssmPDouble xi=1.;
+                                IssmDouble  tau=xi*h/(2*vel);
+
+                                IssmDouble vx,vy;
+                                vx = w[0]; vy = w[1];
+                                /*Mass matrix - part 2*/
+                                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]);
+                                        }
+                                }
+
+                                /*Advection matrix - part 2, 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])*(vx*dbasis[0*numnodes+i]+vy*dbasis[1*numnodes+i]);
+                                        }
+                                }
+
+                                break;
+                        }
                         default:
                                 _error_("unknown type of stabilization in LevelsetAnalysis.cpp");
@@ -330 +354 @@
         xDelete<IssmDouble>(basis);
         xDelete<IssmDouble>(dbasis);
-        xDelete<IssmDouble>(Bprime);
         delete gauss;
         if(basalelement->IsSpawnedElement()){basalelement->DeleteMaterials(); delete basalelement;};
@@ -341 +364 @@
 
         /*Intermediaries */
-        int domaintype;
+        int         domaintype,stabilization;
         IssmDouble  Jdet,dt;
         IssmDouble  lsf;
@@ -348 +371 @@
         /*Fetch number of nodes and dof for this finite element*/
         int numnodes = basalelement->GetNumberOfNodes();
+        basalelement->FindParam(&stabilization,MasstransportStabilizationEnum);
 
         /*Initialize Element vector*/
         ElementVector* pe = basalelement->NewElementVector();
-        basalelement->FindParam(&dt,TimesteppingTimeStepEnum);
-
-        if(dt!=0.){
-                /*Initialize basis vector*/
-                IssmDouble*    basis = xNew<IssmDouble>(numnodes);
-
-                /*Retrieve all inputs and parameters*/
-                basalelement->GetVerticesCoordinates(&xyz_list);
-                Input* levelset_input  = basalelement->GetInput(MaskIceLevelsetEnum); _assert_(levelset_input);
-
-                /* Start  looping on the number of gaussian points: */
-                Gauss* gauss=basalelement->NewGauss(2);
-                while(gauss->next()){
-                        basalelement->JacobianDeterminant(&Jdet,xyz_list,gauss);
-                        basalelement->NodalFunctions(basis,gauss);
-
-                        /* old function value */
-                        levelset_input->GetInputValue(&lsf,gauss);
-                        for(int i=0;i<numnodes;i++) pe->values[i]+=Jdet*gauss->weight*lsf*basis[i];
-                }
-
-                /*Clean up and return*/
-                xDelete<IssmDouble>(xyz_list);
-                xDelete<IssmDouble>(basis);
-                basalelement->FindParam(&domaintype,DomainTypeEnum);
-                if(basalelement->IsSpawnedElement()){basalelement->DeleteMaterials(); delete basalelement;};
-                delete gauss;
-        }
+        basalelement->FindParam(&dt,TimesteppingTimeStepEnum); _assert_(dt>0.);
+
+        /*Initialize basis vector*/
+        IssmDouble*    basis = xNew<IssmDouble>(numnodes);
+        IssmDouble*    dbasis = NULL;
+        if(stabilization==5) dbasis= xNew<IssmDouble>(2*numnodes);
+
+        /*Retrieve all inputs and parameters*/
+        basalelement->GetVerticesCoordinates(&xyz_list);
+        Input* levelset_input = basalelement->GetInput(MaskIceLevelsetEnum); _assert_(levelset_input);
+        Input* mf_vx_input    = basalelement->GetInput(MovingFrontalVxEnum); _assert_(mf_vx_input);
+        Input* mf_vy_input    = basalelement->GetInput(MovingFrontalVyEnum); _assert_(mf_vy_input);
+
+        IssmDouble h=element->CharacteristicLength();
+
+        /* Start  looping on the number of gaussian points: */
+        Gauss* gauss=basalelement->NewGauss(2);
+        while(gauss->next()){
+                basalelement->JacobianDeterminant(&Jdet,xyz_list,gauss);
+                basalelement->NodalFunctions(basis,gauss);
+
+                /* old function value */
+                levelset_input->GetInputValue(&lsf,gauss);
+                for(int i=0;i<numnodes;i++) pe->values[i]+=Jdet*gauss->weight*lsf*basis[i];
+
+                if(stabilization==5){ /*SUPG*/
+         IssmDouble vx,vy,vel;
+                        basalelement->NodalFunctionsDerivatives(dbasis,xyz_list,gauss);
+                        mf_vx_input->GetInputAverage(&vx);
+                        mf_vy_input->GetInputAverage(&vy);
+                        vel=sqrt(vx*vx+vy*vy)+1.e-8;
+                        IssmPDouble xi=1;
+                        IssmDouble  tau=xi*h/(2*vel);
+
+                        /*Force vector - part 2*/
+                        for(int i=0;i<numnodes;i++){
+                                pe->values[i]+=Jdet*gauss->weight*lsf*(tau*vx*dbasis[0*numnodes+i]+tau*vy*dbasis[1*numnodes+i]);
+                        }
+                }
+        }
+
+        /*Clean up and return*/
+        xDelete<IssmDouble>(xyz_list);
+        xDelete<IssmDouble>(basis);
+   xDelete<IssmDouble>(dbasis);
+        basalelement->FindParam(&domaintype,DomainTypeEnum);
+        if(basalelement->IsSpawnedElement()){basalelement->DeleteMaterials(); delete basalelement;};
+        delete gauss;
 
         return pe;