Changeset 23966

Timestamp:

05/31/19 08:58:09 (6 years ago)

Author:

Mathieu Morlighem

Message:

CHG: added alternative to GL advection: streamline diffusion

File:

• issm/trunk-jpl/src/c/analyses/GLheightadvectionAnalysis.cpp (modified) (2 diffs)

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

@@ -61 +61 @@
 
         /*Intermediaries */
+        const IssmPDouble yts = 365*24*3600.;
         int        domaintype,dim;
         IssmDouble Jdet,D_scalar,onboundary;
@@ -104 +105 @@
                 GetBprime(Bprime,element,dim,xyz_list,gauss);
 
-                D_scalar=gauss->weight*Jdet;
-
                 /*Get velocity*/
                 vx_input->GetInputValue(&vx,gauss);
                 vy_input->GetInputValue(&vy,gauss);
-                bc_input->GetInputValue(&onboundary,gauss);
-                if(onboundary>0.){
-                        /*We do not want to advect garbage, make sure only diffusion is applied on boundary*/
-                        vx = 0.; vy = 0.;
-                }
-
-                /*Diffusion */
-                if(sqrt(vx*vx+vy*vy)<1000./31536000.){
-                        IssmPDouble kappa = -10.;
+
+                if(false){
+                        /*Streamline diffusion*/
+                        vel = sqrt(vx*vx+vy*vy);
+                        if(vel<10./yts){
+                                vx = 0.; vy = 0.;
+                                vel = 30./yts*500000.;
+                        }
+
                         for(int i=0;i<numnodes;i++){
                                 for(int j=0;j<numnodes;j++){
-                                        Ke->values[i*numnodes+j] += D_scalar*kappa*(dbasis[0*numnodes+j]*dbasis[0*numnodes+i] + dbasis[1*numnodes+j]*dbasis[1*numnodes+i]);
+                                        Ke->values[i*numnodes+j] += gauss->weight*Jdet*(
+                                                                (vx*dbasis[0*numnodes+i] + vy*dbasis[1*numnodes+i])*(vx*dbasis[0*numnodes+j] + vy*dbasis[1*numnodes+j])
+                                                                + vel/500000.*(dbasis[0*numnodes+i]*dbasis[0*numnodes+j] + dbasis[1*numnodes+i]*dbasis[1*numnodes+j]));
                                 }
                         }
                 }
-
-                /*Advection: */
-                for(int i=0;i<numnodes;i++){
-                        for(int j=0;j<numnodes;j++){
-                                Ke->values[i*numnodes+j] += (D_scalar*(vx*dbasis[0*numnodes+j]*basis[i] + vy*dbasis[1*numnodes+j]*basis[i]))*1e-2;
-                        }
-                }
-
-                /*Artificial diffusivity*/
-                vel=sqrt(vx*vx + vy*vy)+1.e-14;
-                D[0][0]=D_scalar*h/(2.*vel)*fabs(vx*vx);  D[0][1]=D_scalar*h/(2.*vel)*fabs(vx*vy);
-                D[1][0]=D_scalar*h/(2.*vel)*fabs(vy*vx);  D[1][1]=D_scalar*h/(2.*vel)*fabs(vy*vy);
-                TripleMultiply(Bprime,dim,numnodes,1,
-                                        &D[0][0],2,2,0,
-                                        Bprime,dim,numnodes,0,
-                                        &Ke->values[0],1);
+                else{
+                        D_scalar=gauss->weight*Jdet;
+
+                        bc_input->GetInputValue(&onboundary,gauss);
+                        if(onboundary>0.){
+                                /*We do not want to advect garbage, make sure only diffusion is applied on boundary*/
+                                vx = 0.; vy = 0.;
+                        }
+
+                        /*Diffusion */
+                        if(sqrt(vx*vx+vy*vy)<1000./31536000.){
+                                IssmPDouble kappa = -10.;
+                                for(int i=0;i<numnodes;i++){
+                                        for(int j=0;j<numnodes;j++){
+                                                Ke->values[i*numnodes+j] += D_scalar*kappa*(dbasis[0*numnodes+j]*dbasis[0*numnodes+i] + dbasis[1*numnodes+j]*dbasis[1*numnodes+i]);
+                                        }
+                                }
+                        }
+
+                        /*Advection: */
+                        for(int i=0;i<numnodes;i++){
+                                for(int j=0;j<numnodes;j++){
+                                        Ke->values[i*numnodes+j] += (D_scalar*(vx*dbasis[0*numnodes+j]*basis[i] + vy*dbasis[1*numnodes+j]*basis[i]))*1e-2;
+                                }
+                        }
+
+                        /*Artificial diffusivity*/
+                        vel=sqrt(vx*vx + vy*vy)+1.e-14;
+                        D[0][0]=D_scalar*h/(2.*vel)*fabs(vx*vx);  D[0][1]=D_scalar*h/(2.*vel)*fabs(vx*vy);
+                        D[1][0]=D_scalar*h/(2.*vel)*fabs(vy*vx);  D[1][1]=D_scalar*h/(2.*vel)*fabs(vy*vy);
+                        TripleMultiply(Bprime,dim,numnodes,1,
+                                                &D[0][0],2,2,0,
+                                                Bprime,dim,numnodes,0,
+                                                &Ke->values[0],1);
+                }
         }