Changeset 22977

Timestamp:

07/19/18 14:59:24 (7 years ago)

Author:

tpelle

Message:

CHG: advection now working

File:

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

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

@@ -62 +62 @@
         int        domaintype,dim;
         IssmDouble Jdet,D_scalar;
-        IssmDouble vel,vx,vy,dvxdx,dvydy;
-        IssmDouble dvx[2],dvy[2];
+        IssmDouble vel,vx,vy;
         IssmDouble* xyz_list = NULL;
 
@@ -80 +79 @@
         ElementMatrix* Ke     = element->NewElementMatrix();
         IssmDouble*    basis  = xNew<IssmDouble>(numnodes);
-        IssmDouble*    B      = xNew<IssmDouble>(dim*numnodes);
+        IssmDouble*    dbasis = xNew<IssmDouble>(dim*numnodes);
         IssmDouble*    Bprime = xNew<IssmDouble>(dim*numnodes);
-        IssmDouble    D[2][2] = {0.};
+        IssmDouble     D[2][2] = {0.};
 
         /*Retrieve all inputs and parameters*/
         element->GetVerticesCoordinates(&xyz_list);
-        Input* vx_input=element->GetInput(BaseSlopeXEnum); _assert_(vx_input);
-        Input* vy_input=element->GetInput(BaseSlopeYEnum); _assert_(vy_input);
+        //Input* vx_input=element->GetInput(BaseSlopeXEnum); _assert_(vx_input);
+        //Input* vy_input=element->GetInput(BaseSlopeYEnum); _assert_(vy_input);
+        Input* vx_input=element->GetInput(VxEnum); _assert_(vx_input);
+        Input* vy_input=element->GetInput(VyEnum); _assert_(vy_input);
 
         IssmDouble h = element->CharacteristicLength();
 
         /* Start  looping on the number of gaussian points: */
-        Gauss* gauss=element->NewGauss(2);
+        Gauss* gauss=element->NewGauss(4);
         for(int ig=gauss->begin();ig<gauss->end();ig++){
                 gauss->GaussPoint(ig);
 
                 element->JacobianDeterminant(&Jdet,xyz_list,gauss);
-                GetB(B,element,dim,xyz_list,gauss);
+                element->NodalFunctions(basis,gauss);
+                element->NodalFunctionsDerivatives(dbasis,xyz_list,gauss);
                 GetBprime(Bprime,element,dim,xyz_list,gauss);
 
+                D_scalar=gauss->weight*Jdet;
+
+                /*Diffusion */
+                if(sqrt(vx*vx+vy*vy)<10./31536000.){
+                        IssmPDouble kappa = 1.e-4;
+                        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: */
                 vx_input->GetInputValue(&vx,gauss);
                 vy_input->GetInputValue(&vy,gauss);
-                vx_input->GetInputDerivativeValue(&dvx[0],xyz_list,gauss);
-                vy_input->GetInputDerivativeValue(&dvy[0],xyz_list,gauss);
-                dvxdx=dvx[0];
-                dvydy=dvy[1];
-                D_scalar=gauss->weight*Jdet;
-
-                /*Advection*/
-                //D[0][0]=D_scalar*dvxdx;
-                //D[0][1]=0.;
-                //D[1][0]=0.;
-                //D[1][1]=D_scalar*dvydy;
-                //TripleMultiply(B,2,numnodes,1,
-                //                      &D[0][0],2,2,0,
-                //                      B,2,numnodes,0,
-                //                      &Ke->values[0],1);
-
-                //vx = 0.;
-                //vy = -.1;
-                if(sqrt(vx*vx+vy*vy)>0.28){
-                        /*Advection*/
-                        D[0][0]=D_scalar*vx;
-                        D[1][1]=D_scalar*vy;
-                        TripleMultiply(B,2,numnodes,1,
-                                                &D[0][0],2,2,0,
-                                                Bprime,2,numnodes,0,
-                                                &Ke->values[0],1);
-                }
-                else{
-                        /*Diffusion*/
-                        D[0][0]=1e5*h/2.0;
-                        D[1][1]=1e5*h/2.0;
-                        TripleMultiply(Bprime,2,numnodes,1,
-                                                &D[0][0],2,2,0,
-                                                Bprime,2,numnodes,0,
-                                                &Ke->values[0],1);
-                }
-
-                /*Stabilization*/
-                D[0][0]=h/2.0*fabs(vx);
-                D[0][1]=0.;
-                D[1][0]=0.;
-                D[1][1]=h/2.0*fabs(vy);
-                TripleMultiply(Bprime,2,numnodes,1,
+                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]);
+                        }
+                }
+
+                /*Artifficial 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,2,numnodes,0,
+                                        Bprime,dim,numnodes,0,
                                         &Ke->values[0],1);
         }
@@ -152 +135 @@
         /*Clean up and return*/
         xDelete<IssmDouble>(xyz_list);
-        xDelete<IssmDouble>(B);
+        xDelete<IssmDouble>(basis);
+        xDelete<IssmDouble>(dbasis);
         xDelete<IssmDouble>(Bprime);
         delete gauss;