Changeset 17171

Timestamp:

01/24/14 15:57:39 (11 years ago)

Author:

jbondzio

Message:

CHG: Added artificial diffusion to extrapolation routine to stabilize solution

File:

• issm/trunk-jpl/src/c/analyses/ExtrapolationAnalysis.cpp (modified) (6 diffs)

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

@@ -15 +15 @@
 /*}}}*/
 void ExtrapolationAnalysis::UpdateElements(Elements* elements,IoModel* iomodel,int analysis_counter,int analysis_type){/*{{{*/
-        int    stabilization,finiteelement;
+        int    finiteelement;
 
         /*Finite element type*/
@@ -83 +83 @@
         IssmDouble dlevelset[dim],normal[dim];
         IssmDouble norm_dlevelset;
+        IssmDouble hx,hy,hz,kappa;
         IssmDouble* xyz_list = NULL;
 
@@ -101 +102 @@
         /* Start  looping on the number of gaussian points: */
         Gauss* gauss=element->NewGauss(2);
-        for(int ig=gauss->begin();ig<gauss->end();ig++){
+        for(int ig=gauss->begin();ig<gauss->end();ig++){/*{{{*/
                 gauss->GaussPoint(ig);
 
@@ -128 +129 @@
                                         &Ke->values[0],1);
 
+                /* Stabilization *//*{{{*/
                 stabilization=1;
                 if (stabilization==0){/* no stabilization, do nothing*/}
                 else if(stabilization==1){
-                        /*Streamline upwinding*/
-                        for(row=0;row<dim;row++)
-                                for(col=0;col<dim;col++)
-                                        D[row][col]=h/(2.*1.)*normal[row]*normal[col];
-
+                        /* Artificial Diffusion */
+                        element->ElementSizes(&hx,&hy,&hz);
+                        h=sqrt( pow(hx*normal[0],2) + pow(hy*normal[1],2));
+                        kappa=h/2.; 
+                        D[0][0]=D_scalar*kappa;
+                        D[0][1]=0.;
+                        D[1][0]=0.;
+                        D[1][1]=D_scalar*kappa;
                         TripleMultiply(Bprime,dim,numnodes,1,
                                                 &D[0][0],dim,dim,0,
@@ -141 +146 @@
                                                 &Ke->values[0],1);
                 }
-        }
+                else if(stabilization==2){
+                        /*Streamline upwinding - do not use this for extrapolation: yields oscillating results due to smoothing along normal, not across */
+                        for(row=0;row<dim;row++)
+                                for(col=0;col<dim;col++)
+                                        D[row][col]=h/(2.*1.)*normal[row]*normal[col];
+
+                        TripleMultiply(Bprime,dim,numnodes,1,
+                                                &D[0][0],dim,dim,0,
+                                                Bprime,dim,numnodes,0,
+                                                &Ke->values[0],1);
+                }/*}}}*/
+        }/*}}}*/
 
         /*Clean up and return*/
@@ -150 +166 @@
         return Ke;
 
-        }/*}}}*/
+}/*}}}*/
 ElementVector* ExtrapolationAnalysis::CreatePVector(Element* element){/*{{{*/