Changeset 1158 for XIOS/dev/dev_olga/extern/remap

Timestamp:

06/06/17 17:58:16 (7 years ago)

Author:

oabramkina

Message:

Two server levels: merging with trunk r1137.
There are bugs.

Location:

XIOS/dev/dev_olga/extern/remap/src

Files:

• elt.hpp (modified) (1 diff)
• intersect.cpp (modified) (2 diffs)
• intersect.hpp (modified) (1 diff)
• mapper.cpp (modified) (16 diffs)
• mapper.hpp (modified) (2 diffs)
• tree.cpp (modified) (1 diff)
• triple.cpp (modified) (1 diff)
• triple.hpp (modified) (1 diff)

XIOS/dev/dev_olga/extern/remap/src/elt.hpp

@@ -105 +105 @@
         }
 
+  void insert_vertex(int i, const Coord& v)
+  {
+    for(int j=n; j > i ; j--)
+    {
+      vertex[j]=vertex[j-1] ;
+      edge[j]=edge[j-1] ;
+      d[j]=d[j-1] ;
+      neighbour[j]=neighbour[j-1] ;
+    }
+    vertex[i+1]=v ;
+    n++ ;
+  }
+  
         int neighbour[NMAX];
         double d[NMAX]; /**< distance of centre of small circle to origin, zero if great circle */

XIOS/dev/dev_olga/extern/remap/src/intersect.cpp

@@ -37 +37 @@
 }
 
+/** New methods to find an insert a neighbour in a cell of the source mesh.
+ *  return true/false if cell b is a neighbour of a. if "insert" is true, then b will be inserted as a neighbour
+ * in cell a . This is needed for 2 order interpolation that need neighboround for gradient computing.
+ * A cell is a neighbour if :
+ *  - it shares 2 countiguous vertex (ie an edge) with a
+ *  - A vertex of b is located on one of an edge of a.
+ **/
+bool insertNeighbour( Elt& a, const Elt& b, bool insert )
+{
+  // for now suppose pole -> Oz
+  Coord pole(0,0,1) ;
+  Coord O, Oa1, Oa2,Ob1,Ob2,V1,V2 ;
+  double da,db,alpha,alpha1,alpha2,delta ;
+  
+    
+  for (int i = 0; i < a.n; i++)
+  {
+    for (int j = 0; j < b.n; j++)
+    {
+// share a full edge ? be carefull at the orientation
+      assert(squaredist(a.vertex[ i       ], b.vertex[ j       ]) > 1e-10*1e-10 ||
+             squaredist(a.vertex[(i+1)%a.n], b.vertex[(j+1)%b.n]) > 1e-10*1e-10);
+      if (   squaredist(a.vertex[ i       ], b.vertex[ j           ]) < 1e-10*1e-10 &&
+             squaredist(a.vertex[(i+1)%a.n], b.vertex[(j+b.n-1)%b.n]) < 1e-10*1e-10)
+      {
+        if (insert) a.neighbour[i] = b.id.ind ;
+        return true;
+      }
+      
+
+// 1 or 2 vertices of b is located on an edge of a
+       da=a.d[i] ;
+       if (scalarprod(a.edge[i], pole) < 0) da=-da ;
+       db=b.d[(j+b.n-1)%b.n] ;
+       if (scalarprod(b.edge[(j+b.n-1)%b.n], pole) < 0) db=-db ;
+      
+      if ( fabs(da-db)<1e-10 ) 
+      {
+        O=pole*da ;
+        Oa1=a.vertex[i]-O ;
+        Oa2=a.vertex[(i+1)%a.n]-O ; 
+        Ob1=b.vertex[j]-O ;
+        Ob2=b.vertex[(j+b.n-1)%b.n]-O ;
+        V1=crossprod(Oa1,Oa2) ;
+        V2=crossprod(Ob1,Ob2) ;
+        if (norm(crossprod(V1,V2))/(norm(V1)*norm(V2)) < 1e-10)
+        {
+          alpha = vectAngle(Oa1,Oa2,V1) ;
+          alpha1= vectAngle(Oa1,Ob1,V1) ;
+          alpha2= vectAngle(Oa1,Ob2,V1) ;
+          delta= alpha2-alpha1 ;
+          if (delta >= M_PI) delta=2*M_PI-delta ;
+          else if (delta <= -M_PI) delta=2*M_PI+delta ;
+          
+          if (alpha >= 0)
+          {
+            if (alpha1 > 1e-10 && alpha1 < alpha-1e-10)
+            {
+              if (alpha2 > 1e-10 && alpha2 < alpha-1e-10)
+              {
+                assert(delta > 0) ;
+                if (insert)
+                {
+                // insert both
+                  a.insert_vertex(i,b.vertex[(j+b.n-1)%b.n]);
+                  a.insert_vertex(i,b.vertex[j]);
+                  a.neighbour[i+1] = b.id.ind ;
+                }
+                return true ;
+              }
+              else
+              {
+                assert( delta > 0 ) ;
+                if (insert)
+                {
+                //insert alpha1
+                  a.insert_vertex(i,b.vertex[j]);
+                  a.neighbour[i+1] = b.id.ind ;
+                }
+                return true ;
+              }
+            }
+            else if (alpha2 > 1e-10 && alpha2 < alpha-1e-10)
+            {
+              assert( delta > 0 ) ;
+              if (insert)
+              {
+              // insert alpha2
+                a.insert_vertex(i,b.vertex[(j+b.n-1)%b.n]);
+                a.neighbour[i] = b.id.ind ;
+              }
+              return true ;
+            }
+            else
+            {
+              // nothing to do
+            } 
+
+          }
+          else  // alpha < 0
+          {
+            if (alpha1 < -1e-10 && alpha1 > alpha+1e-10)
+            {
+              if (alpha2 < -1e-10 && alpha2 > alpha+1e-10)
+              {
+                assert(delta < 0) ;
+                if (insert)
+                {
+                // insert both
+                  a.insert_vertex(i,b.vertex[(j+b.n-1)%b.n]);
+                  a.insert_vertex(i,b.vertex[j]);
+                  a.neighbour[i+1] = b.id.ind ;
+                }
+                return true ;
+              }
+              else
+              {
+                assert(delta < 0) ;
+                if (insert)
+                {
+                //insert alpha1
+                  a.insert_vertex(i,b.vertex[j]);
+                  a.neighbour[i+1] = b.id.ind ;
+                }
+                return true ;
+              }
+            }
+            else if (alpha2 < -1e-10 && alpha2 > alpha+1e-10)
+            {
+              assert(delta < 0) ;
+              if (insert)
+              {
+              // insert alpha2
+                a.insert_vertex(i,b.vertex[(j+b.n-1)%b.n]);
+                a.neighbour[i] = b.id.ind ;
+              }
+              return true ;
+            }
+            else
+            {
+               // nothing to do
+            }
+          }
+        }
+      }
+    }
+  }
+  return false;
+}
 
 /**
@@ -44 +193 @@
 void set_neighbour(Elt& a, const Elt& b)
 {
-        if (b.id.ind == a.id.ind) return;
-        int idx = neighbour_idx(a, b);
-        if (idx != NOT_FOUND)
-                a.neighbour[idx] = b.id.ind;
+  if (b.id.ind == a.id.ind) return;
+/*
+  int idx = neighbour_idx(a, b);
+  if (idx != NOT_FOUND)
+  a.neighbour[idx] = b.id.ind;
+*/
+  insertNeighbour(a,b,true) ; 
 }
 
 /** return true if `a` and `b` share an edge */
-bool isNeighbour(const Elt& a, const Elt& b)
-{
-        return neighbour_idx(a, b) != NOT_FOUND;
+bool isNeighbour(Elt& a, const Elt& b)
+{
+        // return neighbour_idx(a, b) != NOT_FOUND;
+  return insertNeighbour(a,b,false) ;
 }
 

XIOS/dev/dev_olga/extern/remap/src/intersect.hpp

@@ -2 +2 @@
 
 void set_neighbour(Elt& elt, const Elt& elt2);
-bool isNeighbour(const Elt& elt, const Elt& elt2);
+bool isNeighbour(Elt& elt, const Elt& elt2);
 
 void intersect(Elt *a, Elt *b);

XIOS/dev/dev_olga/extern/remap/src/mapper.cpp

@@ -109 +109 @@
 }
 
-vector<double> Mapper::computeWeights(int interpOrder, bool renormalize)
+vector<double> Mapper::computeWeights(int interpOrder, bool renormalize, bool quantity)
 {
         vector<double> timings;
@@ -152 +152 @@
         if (mpiRank == 0 && verbose) cout << "Remapping..." << endl;
         tic = cputime();
-        nWeights = remap(&targetElements[0], targetElements.size(), interpOrder, renormalize);
+        nWeights = remap(&targetElements[0], targetElements.size(), interpOrder, renormalize, quantity);
         timings.push_back(cputime() - tic);
 
@@ -166 +166 @@
    @param order is the order of interpolaton (must be 1 or 2).
 */
-int Mapper::remap(Elt *elements, int nbElements, int order, bool renormalize)
+int Mapper::remap(Elt *elements, int nbElements, int order, bool renormalize, bool quantity)
 {
         int mpiSize, mpiRank;
@@ -184 +184 @@
         int **sendElement = new int*[mpiSize]; /* indices of elements required from other rank */
         double **recvValue = new double*[mpiSize];
+        double **recvArea = new double*[mpiSize];
         Coord **recvGrad = new Coord*[mpiSize];
         GloId **recvNeighIds = new GloId*[mpiSize]; /* ids of the of the source neighbours which also contribute through gradient */
@@ -205 +206 @@
                         sendElement[rank] = new int[nbSendElement[rank]];
                         recvValue[rank]   = new double[nbSendElement[rank]];
+                        recvArea[rank]    = new double[nbSendElement[rank]];
                         if (order == 2)
                         {
@@ -234 +236 @@
         int **recvElement = new int*[mpiSize];
         double **sendValue = new double*[mpiSize];
+        double **sendArea = new double*[mpiSize];
         Coord **sendGrad = new Coord*[mpiSize];
         GloId **sendNeighIds = new GloId*[mpiSize];
-        MPI_Request *sendRequest = new MPI_Request[3*mpiSize];
-        MPI_Request *recvRequest = new MPI_Request[3*mpiSize];
+        MPI_Request *sendRequest = new MPI_Request[4*mpiSize];
+        MPI_Request *recvRequest = new MPI_Request[4*mpiSize];
         for (int rank = 0; rank < mpiSize; rank++)
         {
@@ -250 +253 @@
                         recvElement[rank] = new int[nbRecvElement[rank]];
                         sendValue[rank]   = new double[nbRecvElement[rank]];
+                        sendArea[rank]   = new double[nbRecvElement[rank]];
                         if (order == 2)
                         {
@@ -263 +267 @@
                 }
         }
-        MPI_Status *status = new MPI_Status[3*mpiSize];
-        MPI_Waitall(nbRecvRequest, recvRequest, status);
+        MPI_Status *status = new MPI_Status[4*mpiSize];
+        
         MPI_Waitall(nbSendRequest, sendRequest, status);
+        MPI_Waitall(nbRecvRequest, recvRequest, status);
 
         /* for all indices that have been received from requesting ranks: pack values and gradients, then send */
@@ -278 +283 @@
                         {
                                 sendValue[rank][j] = sstree.localElements[recvElement[rank][j]].val;
+                                sendArea[rank][j] = sstree.localElements[recvElement[rank][j]].area;
                                 if (order == 2)
                                 {
@@ -297 +303 @@
                         MPI_Issend(sendValue[rank],  nbRecvElement[rank], MPI_DOUBLE, rank, 0, communicator, &sendRequest[nbSendRequest]);
                         nbSendRequest++;
+                        MPI_Issend(sendArea[rank],  nbRecvElement[rank], MPI_DOUBLE, rank, 0, communicator, &sendRequest[nbSendRequest]);
+                        nbSendRequest++;
                         if (order == 2)
                         {
@@ -317 +325 @@
                         MPI_Irecv(recvValue[rank],  nbSendElement[rank], MPI_DOUBLE, rank, 0, communicator, &recvRequest[nbRecvRequest]);
                         nbRecvRequest++;
+                        MPI_Irecv(recvArea[rank],  nbSendElement[rank], MPI_DOUBLE, rank, 0, communicator, &recvRequest[nbRecvRequest]);
+                        nbRecvRequest++;
                         if (order == 2)
                         {
@@ -334 +344 @@
                 }
         }
+        
+        MPI_Waitall(nbSendRequest, sendRequest, status);
         MPI_Waitall(nbRecvRequest, recvRequest, status);
-        MPI_Waitall(nbSendRequest, sendRequest, status);
+        
 
         /* now that all values and gradients are available use them to computed interpolated values on target
@@ -357 +369 @@
                         int rank = (*it)->id.rank;
                         double fk = recvValue[rank][n1];
+                        double srcArea = recvArea[rank][n1];
                         double w = (*it)->area;
+      if (quantity) w/=srcArea ;
 
                         /* first order: src value times weight (weight = supermesh area), later divide by target area */
                         int kk = (order == 2) ? n1 * (NMAX + 1) : n1;
                         GloId neighID = recvNeighIds[rank][kk];
-                        wgt_map[neighID] += (*it)->area;
+                        wgt_map[neighID] += w;
 
                         if (order == 2)
@@ -383 +397 @@
     for (map<GloId,double>::iterator it = wgt_map.begin(); it != wgt_map.end(); it++)
                 {
-                        this->remapMatrix[i] = (it->second / e.area) / renorm;
+      if (quantity)  this->remapMatrix[i] = (it->second ) / renorm;
+                        else this->remapMatrix[i] = (it->second / e.area) / renorm;
                         this->srcAddress[i] = it->first.ind;
                         this->srcRank[i] = it->first.rank;
@@ -400 +415 @@
                         delete[] sendElement[rank];
                         delete[] recvValue[rank];
+                        delete[] recvArea[rank];
                         if (order == 2)
                         {
@@ -410 +426 @@
                         delete[] recvElement[rank];
                         delete[] sendValue[rank];
+                        delete[] sendArea[rank];
                         if (order == 2)
                                 delete[] sendGrad[rank];

XIOS/dev/dev_olga/extern/remap/src/mapper.hpp

@@ -34 +34 @@
        /** @param trgElts are the elements of the unstructured target grid
            Returns the timings for substeps: */
-       vector<double> computeWeights(int interpOrder, bool renormalize=false);
+       vector<double> computeWeights(int interpOrder, bool renormalize=false, bool quantity=false);
        int getNbWeights(void) { return nWeights ; }
 /*
@@ -51 +51 @@
 private:
        /** @return number of weights (local to cpu) */
-       int remap(Elt* elements, int nbElements, int order, bool renormalize=false);
+       int remap(Elt* elements, int nbElements, int order, bool renormalize=false, bool quantity=false);
 
        void buildMeshTopology();

XIOS/dev/dev_olga/extern/remap/src/tree.cpp

@@ -142 +142 @@
         root->parent = 0;
         root->leafCount = 0;
-        root->centre = ORIGIN;
+// initialize root node on the sphere
+  root->centre.x=1 ; root->centre.y=0 ; root->centre.z=0 ; 
         root->radius = 0.;
         root->reinserted = false;

XIOS/dev/dev_olga/extern/remap/src/triple.cpp

@@ -124 +124 @@
 }
 
+// return oriented vector angle in range [-pi..pi], pole must be orthogonal to a and b
+double vectAngle(const Coord &a, const Coord &b, const Coord &pole)
+{ 
+  double nab = 1./(norm(a)*norm(b)) ;
+  
+  Coord a_cross_b=crossprod(a, b)*nab ;
+  double sinVect ;
+  if (scalarprod(a_cross_b, pole) >= 0) sinVect=norm(a_cross_b) ;
+  else sinVect=-norm(a_cross_b) ;
+  double cosVect=scalarprod(a,b)*nab ;
+
+  return atan2(sinVect,cosVect) ;
 }
+
+}

XIOS/dev/dev_olga/extern/remap/src/triple.hpp

@@ -121 +121 @@
 
 double angle(const Coord &a, const Coord &b, const Coord &pole);
+
+// return oriented vector angle in range [-pi..pi], pole must be orthogonal to a and b
+double vectAngle(const Coord &a, const Coord &b, const Coord &pole) ;
+
 void print_Coord(Coord &p);