Changeset 27004

Timestamp:

05/13/22 06:38:23 (3 years ago)

Author:

Cheng Gong

Message:

ADD: stablization=6 for levelset functions using FAB+SUPG

Location:

issm/trunk-jpl/src

Files:

• c/analyses/LevelsetAnalysis.cpp (modified) (6 diffs)
• m/classes/levelset.m (modified) (1 diff)

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

@@ -10 +10 @@
 #include "../modules/modules.h"
 #include "../solutionsequences/solutionsequences.h"
+#include <math.h>
 
 void LevelsetAnalysis::CreateConstraints(Constraints* constraints,IoModel* iomodel){/*{{{*/
@@ -423 +424 @@
                                 break;
                         }
+                        case 6:{
+                                /*SUPG*/
+                                IssmDouble vx,vy;
+                                mf_vx_input->GetInputAverage(&vx);
+                                mf_vy_input->GetInputAverage(&vy);
+                                vel=sqrt(vx*vx+vy*vy)+1.e-8;
+                                IssmDouble ECN, K;
+                                ECN = vel *dt /h;
+                                K = 1./tanh(ECN) - 1./ECN;
+//                              if (ECN<1e-6) K = ECN /3.0;
+
+                                /*According to Hilmar, xi=K is too large*/
+                                IssmPDouble xi=0.1*K;
+
+                                IssmDouble  tau=xi*h/(2*vel);
+                                Input* levelset_input = NULL;
+
+
+                                IssmDouble kappa;
+                                IssmDouble p=4, q=4;
+                                IssmDouble phi[3];
+
+                           levelset_input=basalelement->GetInput(MaskIceLevelsetEnum); _assert_(levelset_input);
+                                levelset_input->GetInputValue(&phi[0], gauss);
+
+                                IssmDouble dphidx=0., dphidy=0.;
+                                IssmDouble nphi;
+
+                                for(int i=0;i<numnodes;i++){
+                                        dphidx += phi[i]*dbasis[0*numnodes+i];
+                                        dphidy += phi[i]*dbasis[1*numnodes+i];
+                                }
+                                nphi = sqrt(dphidx*dphidx+dphidy*dphidy);
+                        
+                                if (nphi >= 1) {
+                                        kappa = 1 - 1.0/nphi;
+                                }
+                                else {
+                                        kappa = 0.5/M_PI *sin(2*M_PI*nphi)/nphi;
+                                }
+
+                                kappa = kappa * vel / h;
+
+                                /*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]);
+                                        }
+                                }
+                                /*Add the pertubation term \nabla\cdot(\kappa*\nabla\phi)*/
+                                for(int i=0;i<numnodes;i++){
+                                        for(int j=0;j<numnodes;j++){
+                                                for(int k=0;k<dim;k++){
+                                                                Ke->values[i*numnodes+j]+= dt*gauss->weight*Jdet*kappa*dbasis[k*numnodes+j]*dbasis[k*numnodes+i];
+                                                }
+                                        }
+                                }
+
+                                break;
+                        }
                         default:
                                 _error_("unknown type of stabilization in LevelsetAnalysis.cpp");
@@ -458 +526 @@
         IssmDouble*    basis = xNew<IssmDouble>(numnodes);
         IssmDouble*    dbasis = NULL;
-        if(stabilization==5) dbasis= xNew<IssmDouble>(2*numnodes);
+        if((stabilization==5) |(stabilization == 6)) dbasis= xNew<IssmDouble>(2*numnodes);
 
         /*Retrieve all inputs and parameters*/
@@ -479 +547 @@
 
                 if(stabilization==5){ /*SUPG*/
-         IssmDouble vx,vy,vel;
+                        IssmDouble vx,vy,vel;
                         basalelement->NodalFunctionsDerivatives(dbasis,xyz_list,gauss);
                         mf_vx_input->GetInputAverage(&vx);
@@ -492 +560 @@
                         }
                 }
+                else if (stabilization ==6) {
+                        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;
+
+                        IssmDouble ECN, K;
+                        ECN = vel *dt /h;
+                        K = 1./tanh(ECN) - 1./ECN;
+        //              if (ECN<1e-6) K = ECN /3.0;
+
+                        /*According to Hilmar, xi=K is too large*/
+                        IssmPDouble xi=0.1*K;
+
+                        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]+vy*dbasis[1*numnodes+i]);
+                        }
+                }
         }
 
@@ -497 +587 @@
         xDelete<IssmDouble>(xyz_list);
         xDelete<IssmDouble>(basis);
-   xDelete<IssmDouble>(dbasis);
+        xDelete<IssmDouble>(dbasis);
         basalelement->FindParam(&domaintype,DomainTypeEnum);
         if(basalelement->IsSpawnedElement()){basalelement->DeleteMaterials(); delete basalelement;};

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

@@ -62 +62 @@
 
                         md = checkfield(md,'fieldname','levelset.spclevelset','Inf',1,'timeseries',1);
-                        md = checkfield(md,'fieldname','levelset.stabilization','values',[0 1 2 5]);
+                        md = checkfield(md,'fieldname','levelset.stabilization','values',[0 1 2 5 6]);
                         md = checkfield(md,'fieldname','levelset.kill_icebergs','numel',1,'values',[0 1]);
                         md = checkfield(md,'fieldname','levelset.migration_max','numel',1,'NaN',1,'Inf',1,'>',0);