Changeset 23462

Timestamp:

11/09/18 09:05:58 (6 years ago)

Author:

Mathieu Morlighem

Message:

BUG: make sure vx and vy are initialized before they are used and force them to be 0 on boundary

File:

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

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

@@ -38 +38 @@
                 iomodel->FetchDataToInput(elements,"md.mesh.vertexonsurface",MeshVertexonsurfaceEnum);
         }
+        iomodel->FetchDataToInput(elements,"md.mesh.vertexonboundary",MeshVertexonboundaryEnum);
 }/*}}}*/
 void GLheightadvectionAnalysis::UpdateParameters(Parameters* parameters,IoModel* iomodel,int solution_enum,int analysis_enum){/*{{{*/
@@ -61 +62 @@
         /*Intermediaries */
         int        domaintype,dim;
-        IssmDouble Jdet,D_scalar;
+        IssmDouble Jdet,D_scalar,onboundary;
         IssmDouble vel,vx,vy;
         IssmDouble* xyz_list = NULL;
@@ -89 +90 @@
    Input* vx_input=element->GetInput(VxEnum); _assert_(vx_input);
         Input* vy_input=element->GetInput(VyEnum); _assert_(vy_input);
+        Input* bc_input=element->GetInput(MeshVertexonboundaryEnum); _assert_(bc_input);
 
         IssmDouble h = element->CharacteristicLength();
@@ -104 +106 @@
                 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;
+                        IssmPDouble kappa = -10.;
                         for(int i=0;i<numnodes;i++){
                                 for(int j=0;j<numnodes;j++){
@@ -113 +124 @@
                         }
                 }
+
                 /*Advection: */
-                vx_input->GetInputValue(&vx,gauss);
-                vy_input->GetInputValue(&vy,gauss);
                 for(int i=0;i<numnodes;i++){
                         for(int j=0;j<numnodes;j++){
@@ -122 +132 @@
                 }
 
-                /*Artifficial diffusivity*/
+                /*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);