Changeset 1220 for XIOS/dev/branch_openmp/extern/remap

Timestamp:

07/20/17 09:18:34 (7 years ago)

Author:

yushan

Message:

test_remap_omp tested on ADA except two fields

Location:

XIOS/dev/branch_openmp/extern/remap/src

Files:

• gridRemap.hpp (modified) (1 diff)
• intersect.cpp (modified) (8 diffs)
• intersection_ym.cpp (modified) (1 diff)
• mapper.cpp (modified) (1 diff)
• parallel_tree.cpp (modified) (1 diff)
• polyg.cpp (modified) (4 diffs)

XIOS/dev/branch_openmp/extern/remap/src/gridRemap.hpp

@@ -14 +14 @@
 Coord readPole(std::istream&);
 
-extern CRemapGrid srcGrid;
-extern CRemapGrid tgtGrid;
+
 
 }

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

@@ -15 +15 @@
 namespace sphereRemap {
 
+extern CRemapGrid srcGrid;
+#pragma omp threadprivate(srcGrid)
+
+extern CRemapGrid tgtGrid;
+#pragma omp threadprivate(tgtGrid)
+
 using namespace std;
 
@@ -21 +27 @@
 int neighbour_idx(const Elt& a, const Elt& b)
 {
-        for (int i = 0; i < a.n; i++)
-        {
-                for (int j = 0; j < b.n; j++)
-                {
-                        assert(squaredist(a.vertex[ i       ], b.vertex[ j       ]) > EPS*EPS ||
-                               squaredist(a.vertex[(i+1)%a.n], b.vertex[(j+1)%b.n]) > EPS*EPS);
-                        if (   squaredist(a.vertex[ i       ], b.vertex[ j           ]) < 1e-13*1e-13 &&
-                               squaredist(a.vertex[(i+1)%a.n], b.vertex[(j+b.n-1)%b.n]) < 1e-13*1e-13)
-                        {
-                                return i;
-                        }
-                }
-        }
-        return NOT_FOUND;
+  for (int i = 0; i < a.n; i++)
+  {
+    for (int j = 0; j < b.n; j++)
+    {
+      assert(squaredist(a.vertex[ i       ], b.vertex[ j       ]) > EPS*EPS ||
+             squaredist(a.vertex[(i+1)%a.n], b.vertex[(j+1)%b.n]) > EPS*EPS);
+      if (   squaredist(a.vertex[ i       ], b.vertex[ j           ]) < 1e-13*1e-13 &&
+             squaredist(a.vertex[(i+1)%a.n], b.vertex[(j+b.n-1)%b.n]) < 1e-13*1e-13)
+      {
+        return i;
+      }
+    }
+  }
+  return NOT_FOUND;
 }
 
@@ -205 +211 @@
 bool isNeighbour(Elt& a, const Elt& b)
 {
-        // return neighbour_idx(a, b) != NOT_FOUND;
+  // return neighbour_idx(a, b) != NOT_FOUND;
   return insertNeighbour(a,b,false) ;
 }
@@ -213 +219 @@
 void intersect(Elt *a, Elt *b)
 {
-        int na = a->n; /* vertices of a */
-        int nb = b->n; /* vertices of b */
-        Coord *c   = new Coord[na+nb];
-        Coord *c2  = new Coord[na+nb];
-        Coord *xc  = new Coord[na+nb];
-        Coord *xc2 = new Coord[na+nb];
-        Coord gc, gc2;
-        double *d = new double[na+nb];
-        double *d2 = new double[na+nb];
-        double are, are2;
-        Ipt ipt[NMAX*NMAX];
-        Ipt ipt2[NMAX*NMAX];
-        ptsec(a, b, ipt);
-        /* make ipt2 transpose of ipt */
-        for (int ii = 0; ii < na; ii++)
-                for (int jj = 0; jj < nb; jj++)
-                        ipt2[jj*na+ii] = ipt[ii*nb+jj];
-        list<Sgm> iscot;
-        recense(a, b, ipt, iscot, 0);
-        recense(b, a, ipt2, iscot, 1);
-
-        int nseg = iscot.size();
-        int nc = 0;
-        int nc2 = 0;
-        while (iscot.size() && nc < 2)
-                nc = assemble(iscot, c, d, xc);
-        while (iscot.size() && nc2 < 2)
-                nc2 = assemble(iscot, c2, d2, xc2);
-//      assert(nseg == nc + nc2 || nseg == 1); // unused segment
+  int na = a->n; /* vertices of a */
+  int nb = b->n; /* vertices of b */
+  Coord *c   = new Coord[na+nb];
+  Coord *c2  = new Coord[na+nb];
+  Coord *xc  = new Coord[na+nb];
+  Coord *xc2 = new Coord[na+nb];
+  Coord gc, gc2;
+  double *d = new double[na+nb];
+  double *d2 = new double[na+nb];
+  double are, are2;
+  Ipt ipt[NMAX*NMAX];
+  Ipt ipt2[NMAX*NMAX];
+  ptsec(a, b, ipt);
+  /* make ipt2 transpose of ipt */
+  for (int ii = 0; ii < na; ii++)
+    for (int jj = 0; jj < nb; jj++)
+      ipt2[jj*na+ii] = ipt[ii*nb+jj];
+  list<Sgm> iscot;
+  recense(a, b, ipt, iscot, 0);
+  recense(b, a, ipt2, iscot, 1);
+
+  int nseg = iscot.size();
+  int nc = 0;
+  int nc2 = 0;
+  while (iscot.size() && nc < 2)
+    nc = assemble(iscot, c, d, xc);
+  while (iscot.size() && nc2 < 2)
+    nc2 = assemble(iscot, c2, d2, xc2);
+//  assert(nseg == nc + nc2 || nseg == 1); // unused segment
 
         if (!(nseg == nc + nc2 || nseg == 1))
@@ -260 +266 @@
 
 //    intersect_ym(a,b) ;
-        if (nc == 1) nc = 0;
-        if (nc2 == 1) nc2 = 0;
-        gc = barycentre(xc, nc);
-        gc2 = barycentre(xc2, nc2);
-        orient(nc, xc, c, d, gc);
-
-        Coord pole = srcGrid.pole;
-        if (pole == ORIGIN) pole = tgtGrid.pole;
-        const double MINBASE = 1e-11;
-        if (nc == 2) /* nc is the number of vertices of super mesh element */
-        {
-                double base = arcdist(xc[0], xc[1]);
+  if (nc == 1) nc = 0;
+  if (nc2 == 1) nc2 = 0;
+  gc = barycentre(xc, nc);
+  gc2 = barycentre(xc2, nc2);
+  orient(nc, xc, c, d, gc);
+
+  Coord pole = srcGrid.pole;
+  if (pole == ORIGIN) pole = tgtGrid.pole;
+  const double MINBASE = 1e-11;
+  if (nc == 2) /* nc is the number of vertices of super mesh element */
+  {
+    double base = arcdist(xc[0], xc[1]);
 cerr << "DID ARRIVE " << base << xc[0] << xc[1] << endl;
-                gc = midpoint(gc, midpointSC(xc[0], xc[1]));
-                /* intersection area `are` must be zero here unless we have one great and one small circle */
-                are = alun(base, fabs(scalarprod(xc[0], pole)));
-        }
-        else
-        {
-                are = airbar(nc, xc, c, d, pole, gc);
-        }
-        if (nc2 == 2)
-        {
-                double base = arcdist(xc2[0], xc2[1]);
+    gc = midpoint(gc, midpointSC(xc[0], xc[1]));
+    /* intersection area `are` must be zero here unless we have one great and one small circle */
+    are = alun(base, fabs(scalarprod(xc[0], pole)));
+  }
+  else
+  {
+    are = airbar(nc, xc, c, d, pole, gc);
+  }
+  if (nc2 == 2)
+  {
+    double base = arcdist(xc2[0], xc2[1]);
 cerr << "DID ARRIVE " << base << xc2[0] << xc2[1] << endl;
-                assert(base > MINBASE);
-                gc2 = midpoint(gc2, midpointSC(xc2[0], xc2[1]));
-                are2 = alun(base, fabs(scalarprod(xc2[0], pole))); // 0
-        }
-        else
-        {
-                are2 = airbar(nc2, xc2, c2, d2, pole, gc2);
-        }
+    assert(base > MINBASE);
+    gc2 = midpoint(gc2, midpointSC(xc2[0], xc2[1]));
+    are2 = alun(base, fabs(scalarprod(xc2[0], pole))); // 0
+  }
+  else
+  {
+    are2 = airbar(nc2, xc2, c2, d2, pole, gc2);
+  }
 
 //  double ym_area=intersect_ym(a,b) ;
 
   if (nc > 1)
-        {
-                /* create one super mesh polygon that src and dest point to */
-                Polyg *is = new Polyg;
-                is->x = gc;
-                is->area = are;
-                is->id = b->id;
-                is->src_id = b->src_id;
-                is->n = nc;
-                (a->is).push_back(is);
-                (b->is).push_back(is);
+  {
+    /* create one super mesh polygon that src and dest point to */
+    Polyg *is = new Polyg;
+    is->x = gc;
+    is->area = are;
+    is->id = b->id;
+    is->src_id = b->src_id;
+    is->n = nc;
+    (a->is).push_back(is);
+    (b->is).push_back(is);
 /*
-          if (  2*fabs(are-ym_area)/(are+ym_area) > 1.1 && ym_area>1e-8)
+    if (  2*fabs(are-ym_area)/(are+ym_area) > 1.1 && ym_area>1e-8)
     {
       cout<<"Big area difference : "<<are<<"  "<<ym_area<<endl ;
@@ -314 +320 @@
     }
 */
-//              cout<<"intersection : "<<are<<" "<< ym_area<<"  diff : "<<fabs(are-ym_area)<<"  ratio : "<<fabs(are-ym_area)/(0.5*(are+ym_area))<<endl ;
-        }
-        if (nc2 > 1)
-        {
-                Polyg *is = new Polyg;
-                is->x = gc2;
-                is->area = are2;
-                is->id = b->id; /* intersection holds id of corresponding source element (see Elt class definition for details about id) */
-                is->src_id = b->src_id;
-                is->n = nc2;
-                (a->is).push_back(is);
-                (b->is).push_back(is);
+//    cout<<"intersection : "<<are<<" "<< ym_area<<"  diff : "<<fabs(are-ym_area)<<"  ratio : "<<fabs(are-ym_area)/(0.5*(are+ym_area))<<endl ;
+  }
+  if (nc2 > 1)
+  {
+    Polyg *is = new Polyg;
+    is->x = gc2;
+    is->area = are2;
+    is->id = b->id; /* intersection holds id of corresponding source element (see Elt class definition for details about id) */
+    is->src_id = b->src_id;
+    is->n = nc2;
+    (a->is).push_back(is);
+    (b->is).push_back(is);
 /*
-    if (        2*fabs(are-ym_area)/(are+ym_area) > 1.1 && ym_area>1e-8 )
+    if (  2*fabs(are-ym_area)/(are+ym_area) > 1.1 && ym_area>1e-8 )
     {
       cout<<"Big area difference : "<<are<<"  "<<ym_area<<endl ;
@@ -333 +339 @@
     }
 */
-//              cout<<"intersection : "<<are2<<" "<< ym_area<<"  diff : "<<fabs(are-ym_area)<<"  ratio : "<<fabs(are-ym_area)/(0.5*(are+ym_area))<<endl ;
-        }
+//    cout<<"intersection : "<<are2<<" "<< ym_area<<"  diff : "<<fabs(are-ym_area)<<"  ratio : "<<fabs(are-ym_area)/(0.5*(are+ym_area))<<endl ;
+  }
 /*
   if (nc<=1 && nc2<=1)
@@ -344 +350 @@
   }
 */
-        delete [] c;
-        delete [] c2;
-        delete [] xc;
-        delete [] xc2;
-        delete [] d;
-        delete [] d2;
-}
-
-}
+  delete [] c;
+  delete [] c2;
+  delete [] xc;
+  delete [] xc2;
+  delete [] d;
+  delete [] d2;
+}
+
+}

XIOS/dev/branch_openmp/extern/remap/src/intersection_ym.cpp

@@ -13 +13 @@
 
 namespace sphereRemap {
+
+extern CRemapGrid srcGrid;
+#pragma omp threadprivate(srcGrid)
+
+extern CRemapGrid tgtGrid;
+#pragma omp threadprivate(tgtGrid)
+
 
 using namespace std;

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

@@ -14 +14 @@
 
 namespace sphereRemap {
+
+extern CRemapGrid srcGrid;
+#pragma omp threadprivate(srcGrid)
+
+extern CRemapGrid tgtGrid;
+#pragma omp threadprivate(tgtGrid)
+
 
 /* A subdivition of an array into N sub-arays

XIOS/dev/branch_openmp/extern/remap/src/parallel_tree.cpp

@@ -14 +14 @@
 
 namespace sphereRemap {
+
+extern CRemapGrid srcGrid;
+#pragma omp threadprivate(srcGrid)
+
+extern CRemapGrid tgtGrid;
+#pragma omp threadprivate(tgtGrid)
 
 static const int assignLevel = 2;

XIOS/dev/branch_openmp/extern/remap/src/polyg.cpp

@@ -7 +7 @@
 
 #include "polyg.hpp"
+#include <stdio.h>
 
 namespace sphereRemap {
@@ -45 +46 @@
 Coord barycentre(const Coord *x, int n)
 {
+
         if (n == 0) return ORIGIN;
         Coord bc = ORIGIN;
         for (int i = 0; i < n; i++)
+        {
                 bc = bc + x[i];
+        }       
         /* both distances can be equal down to roundoff when norm(bc) < mashineepsilon 
            which can occur when weighted with tiny area */
-
-        assert(squaredist(bc, proj(bc)) <= squaredist(bc, proj(bc * (-1.0))));
-        //if (squaredist(bc, proj(bc)) > squaredist(bc, proj(bc * (-1.0)))) return proj(bc * (-1.0));
+  
+        
+  //assert(squaredist(bc, proj(bc)) <= squaredist(bc, proj(bc * (-1.0))));      
 
         return proj(bc);
@@ -173 +177 @@
                 t[1] = x[i];
                 t[2] = x[ii];
-                double sc=scalarprod(crossprod(t[1] - t[0], t[2] - t[0]), t[0]) ;
+    double sc=scalarprod(crossprod(t[1] - t[0], t[2] - t[0]), t[0]) ;
                 assert(sc >= -1e-10); // Error: tri a l'env (wrong orientation)
                 double area_gc = triarea(t[0], t[1], t[2]);
+                if(area_gc<=0) printf("area_gc = %e\n", area_gc);
                 double area_sc_gc_moon = 0;
                 if (d[i]) /* handle small circle case */
@@ -183 +188 @@
                         char sgl = (mext > 0) ? -1 : 1;
                         area_sc_gc_moon = sgl * alun(arcdist(t[1], t[2]), fabs(scalarprod(t[1], pole)));
+                        //if(area_sc_gc_moon<=0) printf("area_sc_gc_moon = %e\n", area_sc_gc_moon);
                         gg_exact = gg_exact + sc_gc_moon_normalintegral(t[1], t[2], pole);
                 }
                 area += area_gc + area_sc_gc_moon; /* for "spherical circle segment" sum triangular part (at) and "small moon" and => account for small circle */
                 g[i] = barycentre(t, 3) * (area_gc + area_sc_gc_moon);
+                //printf("g[%d] = (%e,%e,%e) * (%e+%e) = (%e,%e,%e) norm = %e\n", i, barycentre(t, 3).x, barycentre(t, 3).y, barycentre(t, 3).z, area_gc,  area_sc_gc_moon, g[i].x, g[i].y, g[i].z, norm(g[i]));
         }
         gg = barycentre(g, N);