Debug fv3

docs/html/mnt_regrid_edges_8h.html

@@ -223,6 +223,7 @@ <h2 class="memtitle"><span class="permalink"><a href="#af6d9017b55316e2236513b7e
   </table>
   </dd>
 </dl>
+<dl class="section note"><dt>Note</dt><dd>the orientation of edges is left to right, bottom to top when using CELL_BY_CELL 3---&gt;&mdash;2 | | ^ ^ | | 0---&gt;&mdash;1 The line integrals should be compatible with the orientation of the edges. </dd></dl>
 <dl class="section return"><dt>Returns</dt><dd>error code (0 is OK) </dd></dl>
 
 </div>

docs/html/mnt_regrid_edges_8h_source.html

@@ -237,29 +237,29 @@
 <div class="line"><a id="l00359" name="l00359"></a><span class="lineno">  359</span>LIBRARY_API</div>
 <div class="line"><a id="l00360" name="l00360"></a><span class="lineno"><a class="line" href="mnt_regrid_edges_8h.html#a81adced614297c90ede7bc25c341eadb">  360</a></span><span class="keywordtype">int</span> <a class="code hl_function" href="mnt_regrid_edges_8h.html#a81adced614297c90ede7bc25c341eadb">mnt_regridedges_getNumDstEdges</a>(<a class="code hl_struct" href="struct_regrid_edges__t.html">RegridEdges_t</a>** self, std::size_t* n);</div>
 <div class="line"><a id="l00361" name="l00361"></a><span class="lineno">  361</span></div>
-<div class="line"><a id="l00370" name="l00370"></a><span class="lineno">  370</span>LIBRARY_API</div>
-<div class="line"><a id="l00371" name="l00371"></a><span class="lineno"><a class="line" href="mnt_regrid_edges_8h.html#af6d9017b55316e2236513b7e81811a62">  371</a></span><span class="keywordtype">int</span> <a class="code hl_function" href="mnt_regrid_edges_8h.html#af6d9017b55316e2236513b7e81811a62">mnt_regridedges_apply</a>(<a class="code hl_struct" href="struct_regrid_edges__t.html">RegridEdges_t</a>** self, </div>
-<div class="line"><a id="l00372" name="l00372"></a><span class="lineno">  372</span>                          <span class="keyword">const</span> <span class="keywordtype">double</span> src_data[], <span class="keywordtype">double</span> dst_data[], <span class="keywordtype">int</span> placement);</div>
-<div class="line"><a id="l00373" name="l00373"></a><span class="lineno">  373</span></div>
-<div class="line"><a id="l00384" name="l00384"></a><span class="lineno">  384</span>LIBRARY_API</div>
-<div class="line"><a id="l00385" name="l00385"></a><span class="lineno"><a class="line" href="mnt_regrid_edges_8h.html#afe60e4a2e763c876c5e4310acae7c3bf">  385</a></span><span class="keywordtype">int</span> <a class="code hl_function" href="mnt_regrid_edges_8h.html#afe60e4a2e763c876c5e4310acae7c3bf">mnt_regridedges_vectorApply</a>(<a class="code hl_struct" href="struct_regrid_edges__t.html">RegridEdges_t</a>** self, </div>
-<div class="line"><a id="l00386" name="l00386"></a><span class="lineno">  386</span>                          <span class="keyword">const</span> <span class="keywordtype">double</span> src_u[], <span class="keyword">const</span> <span class="keywordtype">double</span> src_v[],</div>
-<div class="line"><a id="l00387" name="l00387"></a><span class="lineno">  387</span>                          <span class="keywordtype">double</span> dst_u[], <span class="keywordtype">double</span> dst_v[],</div>
-<div class="line"><a id="l00388" name="l00388"></a><span class="lineno">  388</span>                          <span class="keywordtype">int</span> fs);</div>
-<div class="line"><a id="l00389" name="l00389"></a><span class="lineno">  389</span></div>
-<div class="line"><a id="l00399" name="l00399"></a><span class="lineno">  399</span>LIBRARY_API</div>
-<div class="line"><a id="l00400" name="l00400"></a><span class="lineno"><a class="line" href="mnt_regrid_edges_8h.html#ac2900dcefc3194b04449df18956bb235">  400</a></span><span class="keywordtype">int</span> <a class="code hl_function" href="mnt_regrid_edges_8h.html#ac2900dcefc3194b04449df18956bb235">mnt_regridedges_loadWeights</a>(<a class="code hl_struct" href="struct_regrid_edges__t.html">RegridEdges_t</a>** self, </div>
-<div class="line"><a id="l00401" name="l00401"></a><span class="lineno">  401</span>                                <span class="keyword">const</span> <span class="keywordtype">char</span>* fort_filename, <span class="keywordtype">int</span> n);</div>
-<div class="line"><a id="l00402" name="l00402"></a><span class="lineno">  402</span></div>
-<div class="line"><a id="l00411" name="l00411"></a><span class="lineno">  411</span>LIBRARY_API</div>
-<div class="line"><a id="l00412" name="l00412"></a><span class="lineno"><a class="line" href="mnt_regrid_edges_8h.html#a601f00a4c33b27c4121603f63beea44c">  412</a></span><span class="keywordtype">int</span> <a class="code hl_function" href="mnt_regrid_edges_8h.html#a601f00a4c33b27c4121603f63beea44c">mnt_regridedges_dumpWeights</a>(<a class="code hl_struct" href="struct_regrid_edges__t.html">RegridEdges_t</a>** self, </div>
-<div class="line"><a id="l00413" name="l00413"></a><span class="lineno">  413</span>                                <span class="keyword">const</span> <span class="keywordtype">char</span>* fort_filename, <span class="keywordtype">int</span> n);</div>
-<div class="line"><a id="l00414" name="l00414"></a><span class="lineno">  414</span></div>
-<div class="line"><a id="l00420" name="l00420"></a><span class="lineno">  420</span>LIBRARY_API</div>
-<div class="line"><a id="l00421" name="l00421"></a><span class="lineno"><a class="line" href="mnt_regrid_edges_8h.html#a3cdc2e4c1523d2149b0f67d3bd926403">  421</a></span><span class="keywordtype">int</span> <a class="code hl_function" href="mnt_regrid_edges_8h.html#a3cdc2e4c1523d2149b0f67d3bd926403">mnt_regridedges_print</a>(<a class="code hl_struct" href="struct_regrid_edges__t.html">RegridEdges_t</a>** self);</div>
-<div class="line"><a id="l00422" name="l00422"></a><span class="lineno">  422</span> </div>
-<div class="line"><a id="l00423" name="l00423"></a><span class="lineno">  423</span> </div>
-<div class="line"><a id="l00424" name="l00424"></a><span class="lineno">  424</span><span class="preprocessor">#endif </span><span class="comment">// MNT_REGRID_EDGES</span></div>
+<div class="line"><a id="l00377" name="l00377"></a><span class="lineno">  377</span>LIBRARY_API</div>
+<div class="line"><a id="l00378" name="l00378"></a><span class="lineno"><a class="line" href="mnt_regrid_edges_8h.html#af6d9017b55316e2236513b7e81811a62">  378</a></span><span class="keywordtype">int</span> <a class="code hl_function" href="mnt_regrid_edges_8h.html#af6d9017b55316e2236513b7e81811a62">mnt_regridedges_apply</a>(<a class="code hl_struct" href="struct_regrid_edges__t.html">RegridEdges_t</a>** self, </div>
+<div class="line"><a id="l00379" name="l00379"></a><span class="lineno">  379</span>                          <span class="keyword">const</span> <span class="keywordtype">double</span> src_data[], <span class="keywordtype">double</span> dst_data[], <span class="keywordtype">int</span> placement);</div>
+<div class="line"><a id="l00380" name="l00380"></a><span class="lineno">  380</span></div>
+<div class="line"><a id="l00391" name="l00391"></a><span class="lineno">  391</span>LIBRARY_API</div>
+<div class="line"><a id="l00392" name="l00392"></a><span class="lineno"><a class="line" href="mnt_regrid_edges_8h.html#afe60e4a2e763c876c5e4310acae7c3bf">  392</a></span><span class="keywordtype">int</span> <a class="code hl_function" href="mnt_regrid_edges_8h.html#afe60e4a2e763c876c5e4310acae7c3bf">mnt_regridedges_vectorApply</a>(<a class="code hl_struct" href="struct_regrid_edges__t.html">RegridEdges_t</a>** self, </div>
+<div class="line"><a id="l00393" name="l00393"></a><span class="lineno">  393</span>                          <span class="keyword">const</span> <span class="keywordtype">double</span> src_u[], <span class="keyword">const</span> <span class="keywordtype">double</span> src_v[],</div>
+<div class="line"><a id="l00394" name="l00394"></a><span class="lineno">  394</span>                          <span class="keywordtype">double</span> dst_u[], <span class="keywordtype">double</span> dst_v[],</div>
+<div class="line"><a id="l00395" name="l00395"></a><span class="lineno">  395</span>                          <span class="keywordtype">int</span> fs);</div>
+<div class="line"><a id="l00396" name="l00396"></a><span class="lineno">  396</span></div>
+<div class="line"><a id="l00406" name="l00406"></a><span class="lineno">  406</span>LIBRARY_API</div>
+<div class="line"><a id="l00407" name="l00407"></a><span class="lineno"><a class="line" href="mnt_regrid_edges_8h.html#ac2900dcefc3194b04449df18956bb235">  407</a></span><span class="keywordtype">int</span> <a class="code hl_function" href="mnt_regrid_edges_8h.html#ac2900dcefc3194b04449df18956bb235">mnt_regridedges_loadWeights</a>(<a class="code hl_struct" href="struct_regrid_edges__t.html">RegridEdges_t</a>** self, </div>
+<div class="line"><a id="l00408" name="l00408"></a><span class="lineno">  408</span>                                <span class="keyword">const</span> <span class="keywordtype">char</span>* fort_filename, <span class="keywordtype">int</span> n);</div>
+<div class="line"><a id="l00409" name="l00409"></a><span class="lineno">  409</span></div>
+<div class="line"><a id="l00418" name="l00418"></a><span class="lineno">  418</span>LIBRARY_API</div>
+<div class="line"><a id="l00419" name="l00419"></a><span class="lineno"><a class="line" href="mnt_regrid_edges_8h.html#a601f00a4c33b27c4121603f63beea44c">  419</a></span><span class="keywordtype">int</span> <a class="code hl_function" href="mnt_regrid_edges_8h.html#a601f00a4c33b27c4121603f63beea44c">mnt_regridedges_dumpWeights</a>(<a class="code hl_struct" href="struct_regrid_edges__t.html">RegridEdges_t</a>** self, </div>
+<div class="line"><a id="l00420" name="l00420"></a><span class="lineno">  420</span>                                <span class="keyword">const</span> <span class="keywordtype">char</span>* fort_filename, <span class="keywordtype">int</span> n);</div>
+<div class="line"><a id="l00421" name="l00421"></a><span class="lineno">  421</span></div>
+<div class="line"><a id="l00427" name="l00427"></a><span class="lineno">  427</span>LIBRARY_API</div>
+<div class="line"><a id="l00428" name="l00428"></a><span class="lineno"><a class="line" href="mnt_regrid_edges_8h.html#a3cdc2e4c1523d2149b0f67d3bd926403">  428</a></span><span class="keywordtype">int</span> <a class="code hl_function" href="mnt_regrid_edges_8h.html#a3cdc2e4c1523d2149b0f67d3bd926403">mnt_regridedges_print</a>(<a class="code hl_struct" href="struct_regrid_edges__t.html">RegridEdges_t</a>** self);</div>
+<div class="line"><a id="l00429" name="l00429"></a><span class="lineno">  429</span> </div>
+<div class="line"><a id="l00430" name="l00430"></a><span class="lineno">  430</span> </div>
+<div class="line"><a id="l00431" name="l00431"></a><span class="lineno">  431</span><span class="preprocessor">#endif </span><span class="comment">// MNT_REGRID_EDGES</span></div>
 <div class="ttc" id="amnt_grid_8h_html"><div class="ttname"><a href="mnt_grid_8h.html">mntGrid.h</a></div></div>
 <div class="ttc" id="amnt_regrid_edges_8h_html_a25b0237e7a1851454eea8cba1108f6db"><div class="ttname"><a href="mnt_regrid_edges_8h.html#a25b0237e7a1851454eea8cba1108f6db">mnt_regridedges_del</a></div><div class="ttdeci">LIBRARY_API int mnt_regridedges_del(RegridEdges_t **self)</div><div class="ttdoc">Destructor.</div></div>
 <div class="ttc" id="amnt_regrid_edges_8h_html_a2d5bffef774b78622d203fbe0c09220d"><div class="ttname"><a href="mnt_regrid_edges_8h.html#a2d5bffef774b78622d203fbe0c09220d">mnt_regridedges_setDstGridFlags</a></div><div class="ttdeci">LIBRARY_API int mnt_regridedges_setDstGridFlags(RegridEdges_t **self, int fixLonAcrossDateline, int averageLonAtPole)</div><div class="ttdoc">Set the destination grid flags.</div></div>

mint/tests/test_vector_interp_cubedsphere.py

@@ -74,7 +74,10 @@
 import numpy
 import pytest
 
-from mint import VectorInterp, CELL_BY_CELL_DATA
+from mint import Grid, VectorInterp, CELL_BY_CELL_DATA
+from tempfile import TemporaryDirectory
+
+DATA_DIR = Path(__file__).absolute().parent.parent.parent / Path('data')
 
 # scripts/generate_cubedsphere_grid.py is not part of the mint package;
 # reach it via sys.path, as in test_cubedsphere.py.
@@ -236,6 +239,33 @@ def test_vector_interp_convergence():
             f'({e_prev:.3e} -> {e_next:.3e}): the reconstruction is not converging')
 
 
+def test_lonlat_vector_interp_fv3_grid():
+    # Read in the FV3 cubed-sphere grid.
+    src_grid = Grid()
+    src_grid.setFlags(fixLonAcrossDateline=1, averageLonAtPole=1, degrees=True)
+    filename = str(DATA_DIR / Path('C24_SCRIP_desc.181018.nc'))
+    src_grid.loadFromFV32DFile(filename)
+
+    # Construct target points that are regularly spaced
+    # in lon-lat coordinates
+    lon_reg_A = numpy.linspace(1, 359, 180)
+    lat_reg_A = numpy.linspace(-89, 89, 90)
+
+    LAT_REG_A, LON_REG_A = numpy.meshgrid(lat_reg_A, lon_reg_A)
+    ntarget = LON_REG_A.shape[0] * LAT_REG_A.shape[1]
+    target_points_A = numpy.zeros((ntarget, 3), numpy.float64)
+    target_points_A[:, 0] = LON_REG_A.flat
+    target_points_A[:, 1] = LAT_REG_A.flat  
+
+    # Check the points are valid
+    vi = VectorInterp()
+    vi.setGrid(src_grid)
+    vi.buildLocator(numCellsPerBucket=128, periodX=360., enableFolding=True)
+    numbad = vi.findPoints(target_points_A, tol2=1.e-12)    
+
+    assert numbad == 0, f'{numbad} bad cells found'
+
+
 if __name__ == '__main__':
     test_exact_gradient_matches_finite_difference()
     for M in (4, 8, 16, 32):

src/MvMat_double.cpp

@@ -29,41 +29,37 @@ ColMat<double> eye(std::size_t n){
 
 ***********************************************************************/
 
-ColMat<double> load(const std::string& cfile) { 
+ColMat<double> load(const std::string& cfile) {
 
-  // copy to temp1 
+  // strip lines starting with % (comments), in memory -- no temp files, no shelling
+  // out to cp/sed/rm (which also meant cfile was being pasted, unescaped, into a shell
+  // command)
+  ifstream file(cfile.c_str());
+  std::ostringstream stripped;
+  std::string line;
+  while (std::getline(file, line)) {
+    if (line.empty() || line[0] != '%') {
+      stripped << line << '\n';
+    }
+  }
 
-  std::string temp1 = "temp1"; 
-  std::string s1 = "cp " + cfile + " " + temp1 + '\n'; 
-  system(s1.c_str()); 
+  std::istringstream in(stripped.str());
 
-  // get rid of lines starting with % (comments) 
+  unsigned int nr, nc;
+  in >> nr;
+  in >> nc;
 
-  std::string temp2 = "temp2"; 
-  std::string s2 = "sed 's;^%.*$;;g' " + temp1 + " > " + temp2 + '\n'; 
-  system(s2.c_str()); 
-
-  const char *c_temp2 = temp2.c_str(); 
-  ifstream file(c_temp2); 
 
-  unsigned int nr, nc; 
-  file >> nr; 
-  file >> nc; 
+  ColMat<double> a(nr, nc, 0.);
 
+  for(unsigned int i = 0; i < nr; ++i){
+    for(unsigned int j = 0; j < nc; ++j){
+      in >> a(i,j);
+    }
+  }
 
-  ColMat<double> a(nr, nc, 0.); 
-
-  for(unsigned int i = 0; i < nr; ++i){ 
-    for(unsigned int j = 0; j < nc; ++j){ 
-      file >> a(i,j); 
-    } 
-  } 
-
-  std::string s3 = "rm " + temp1 + ' ' + temp2 + '\n'; 
-  system(s3.c_str()); 
-
-  return a; 
-} 
+  return a;
+}
 
 // real matrix dot complex vector
 Vec_cmplx dot(const Mat &a, const Vec_cmplx &b)

src/mntPolylineIntegral.cpp

@@ -61,6 +61,14 @@ int mnt_polylineintegral_buildLocator(PolylineIntegral_t** self,
     if (enableFolding == 1) {
         (*self)->loc->enableFolding();
     }
+    // NOTE: deliberately not calling setCubedSphere here. containsPoint's spherical
+    // treatment (see vmtCellLocator.h) is only consistent with FindCell/VectorInterp,
+    // whose pcoords/weights come from the very same spherical model. PolylineIntegral's
+    // line/cell-edge intersections (collectIntersectionPoints) use a separate, still
+    // flat-(lon,lat)-straight-line algorithm (LineLineIntersector) that setCubedSphere
+    // would put out of sync with containsPoint near a pole -- fixing that would mean
+    // reworking the intersection math itself to use great-circle arcs, not just this
+    // locator, so it's left as flat/unchanged here rather than partially fixed.
     (*self)->loc->BuildLocator();
 
     return ier;

src/mntRegridEdges.cpp

@@ -562,6 +562,14 @@ int mnt_regridedges_buildLocator(RegridEdges_t** self, int numCellsPerBucket,
     if (enableFolding == 1) {
         (*self)->srcLoc->enableFolding();
     }
+    // NOTE: deliberately not calling setCubedSphere here. containsPoint's spherical
+    // treatment (see vmtCellLocator.h) is only consistent with FindCell/VectorInterp,
+    // whose pcoords/weights come from the very same spherical model. RegridEdges'
+    // line/cell-edge intersections (collectIntersectionPoints) use a separate, still
+    // flat-(lon,lat)-straight-line algorithm (LineLineIntersector) that setCubedSphere
+    // would put out of sync with containsPoint near a pole -- fixing that would mean
+    // reworking the intersection math itself to use great-circle arcs, not just this
+    // locator, so it's left as flat/unchanged here rather than partially fixed.
     (*self)->srcLoc->BuildLocator();
 
     return 0;

src/mntVectorInterp.cpp

@@ -56,6 +56,11 @@ int mnt_vectorinterp_buildLocator(VectorInterp_t** self, int numCellsPerBucket,
     if (enableFolding == 1) {
         (*self)->locator->enableFolding();
     }
+    // fixLonAcrossDateline and averageLonAtPole are only ever both set for a gnomonic
+    // cubed-sphere grid (see Grid_t/mnt_grid_setFlags) -- a plain (possibly rotated)
+    // lon-lat grid sets neither, even when it reaches a pole
+    (*self)->locator->setCubedSphere((*self)->grid->fixLonAcrossDateline &&
+                                      (*self)->grid->averageLonAtPole);
     (*self)->locator->BuildLocator();
 
     return 0;