Changeset 18491

Timestamp:

09/10/14 09:12:28 (11 years ago)

Author:

abuzzi

Message:

CNG: code more constrained

File:

• issm/trunk-jpl/src/m/contrib/gravity/vfsa_mpi.cpp (modified) (3 diffs)

issm/trunk-jpl/src/m/contrib/gravity/vfsa_mpi.cpp ¶

@@ -169 +169 @@
         // Pp->Echo();
 
-        double  rhoi = 890;           /* ice density     */
+        double  rhoi = 917;           /* ice density     */
         double  rhow = 1030;          /* water density   */
         // double  rhos = 2013;               /* sediment density */
@@ -186 +186 @@
 
         double ctr=1;            /* parameter for filtering */
-        double sd=0.1;
+        double sd=0.5;
 
         Matrix *xobs= new Matrix(ny,nx);
@@ -669 +669 @@
         Matrix* m1grsm=new Matrix(my,mx);
         Matrix* m1col=new Matrix(mx*my,1);
+        Matrix* m1gr2=new Matrix(my,mx);
+        Matrix* m1grsm2=new Matrix(my,mx);
+        Matrix* m1col2=new Matrix(mx*my,1);
+        Matrix* nptflag= new Matrix(mx*my,1);
         double u=0;
         double y=0;
         m1->MatrixEqual(m0);
+        nptflag->MatrixSum(icethick,bathy);
+        /* first layer: ice */
         for (int i=0;i<mx*my;i++){
-                if(icethick->GetValue(i,0)==0){
-                        u=double (rand())/ double(RAND_MAX);
+                if(nptflag->GetValue(i,0)==0){
+                        u=double(rand())/double(RAND_MAX);
                         y=signe(u-0.5)*T*(pow(1+1/T,fabs(2*u-1))-1);
                         m1->SetValue(i,1,m0->GetValue(i,1)+y*ptval);
+                        m1->SetValue(i,2,m1->GetValue(i,1)+1e-10);
                         if(m1->GetValue(i,1)<=m1->GetValue(i,0)){
                                 m1->SetValue(i,1,m1->GetValue(i,0)+1e-10);
+                                m1->SetValue(i,2,m1->GetValue(i,1)+1e-10);
                         }
                         if(m1->GetValue(i,1)>=m1->GetValue(i,0)+mmax){
                                 m1->SetValue(i,1,m1->GetValue(i,0)+mmax);
-                        }
-                }
-        }
+                                m1->SetValue(i,2,m1->GetValue(i,1)+1e-10);
+                        }
+                }
+        }
+
         m1->ExtractColumn(m1col,1);
         vec2gridsimple(m1col,m1gr,mx,my);
         filtergrav(m1grsm,m1gr,ctr,sd,mx,my);
         reshape(m1grsm,m1col,mx,my);
+        m1->ExtractColumn(m1col2,2);
+        vec2gridsimple(m1col2,m1gr2,mx,my);
+        filtergrav(m1grsm2,m1gr2,ctr,sd,mx,my);
+        reshape(m1grsm2,m1col2,mx,my);
+
         for (int i=0;i<mx*my;i++){
-                if(icethick->GetValue(i,0)==0){
+                if(nptflag->GetValue(i,0)==0){
                         m1->SetValue(i,1,m1col->GetValue(i,0));
+                        m1->SetValue(i,2,m1col2->GetValue(i,0));
+                        if(m1->GetValue(i,1)<m1->GetValue(i,0)){
+                                m1->SetValue(i,1,m1->GetValue(i,0)+1e-10);
+                                m1->SetValue(i,2,m1->GetValue(i,1)+1e-10);
+                        }
                 }
                 else{
                         m1->SetValue(i,1,m0->GetValue(i,1));
-                }
-                if(m1->GetValue(i,1)<=m1->GetValue(i,0)){
-                        m1->SetValue(i,1,m1->GetValue(i,0)+1e-10);
-                }
-        }
-
-        for (int i=0;i<mx*my;i++){
-                if(bathy->GetValue(i,0)==0){
-                        u=double (rand())/ double(RAND_MAX);
-                        y=signe(u-0.5)*T*(pow(1+1/T,fabs(2*u-1))-1);
-                        m1->SetValue(i,2,m0->GetValue(i,2)+y*ptval);
-                        if(m1->GetValue(i,2)<=m1->GetValue(i,1)){
-                                m1->SetValue(i,2,m1->GetValue(i,1)+1e-10);
-                        }
-                        if(m1->GetValue(i,2)>=m1->GetValue(i,1)+mmax2){
-                                m1->SetValue(i,2,m1->GetValue(i,1)+mmax2);
-                        }
-                }
-        }
-        m1->ExtractColumn(m1col,2);
-        vec2gridsimple(m1col,m1gr,mx,my);
-        filtergrav(m1grsm,m1gr,ctr,sd,mx,my);
-        reshape(m1grsm,m1col,mx,my);
-        for (int i=0;i<mx*my;i++){
-                if(bathy->GetValue(i,0)==0){
-                        m1->SetValue(i,2,m1col->GetValue(i,0));
-                }
-                else{
                         m1->SetValue(i,2,m0->GetValue(i,2));
                 }
-                if(m1->GetValue(i,2)<=m1->GetValue(i,1)){
-                        m1->SetValue(i,2,m1->GetValue(i,1)+1e-10);
-                }
-        }
+        }
+
+                                /* second layer: water */
+//      for (int i=0;i<mx*my;i++){
+//              if(bathy->GetValue(i,0)==0){
+//                      u=double (rand())/ double(RAND_MAX);
+//                      y=signe(u-0.5)*T*(pow(1+1/T,fabs(2*u-1))-1);
+//                      m1->SetValue(i,2,m0->GetValue(i,2)+y*ptval);
+//                      if(m1->GetValue(i,2)<=m1->GetValue(i,1)){
+//                              m1->SetValue(i,2,m1->GetValue(i,1)+1e-10);
+//                      }
+//                      if(m1->GetValue(i,2)>=m1->GetValue(i,1)+mmax2){
+//                              m1->SetValue(i,2,m1->GetValue(i,1)+mmax2);
+//                      }
+//              }
+//      }
+//      m1->ExtractColumn(m1col,2);
+//      vec2gridsimple(m1col,m1gr,mx,my);
+//      filtergrav(m1grsm,m1gr,ctr,sd,mx,my);
+//      reshape(m1grsm,m1col,mx,my);
+//      for (int i=0;i<mx*my;i++){
+//              if(bathy->GetValue(i,0)==0){
+//                      m1->SetValue(i,2,m1col->GetValue(i,0));
+//              }
+//              else{
+//                      m1->SetValue(i,2,m0->GetValue(i,2));
+//              }
+//              if(m1->GetValue(i,2)<=m1->GetValue(i,1)){
+//                      m1->SetValue(i,2,m1->GetValue(i,1)+1e-10);
+//              }
+//      }
         delete m1gr;
         delete m1grsm;
         delete m1col;
+        delete m1gr2;
+        delete m1grsm2;
+        delete m1col2;
+        delete nptflag;
 }/*}}}*/
 double signe(double a){/*{{{*/