From 73fdfda0cfd203f7c3b9e52140e68e1df27dd715 Mon Sep 17 00:00:00 2001 From: Bill Dolinar Date: Fri, 14 Aug 2026 08:26:55 -0600 Subject: [PATCH 01/10] Add a trace benchmark for the follow-flow-path work The "follow flow path" vector display option is being routed through XmGridTrace so the traced path integrates a time-varying field. That only works if the tracer is fast enough to run for every visible glyph, so this establishes what it costs today. testTraceBenchmark traces N seeds over a 200x200 quad grid loaded with two timesteps of a vortex field whose rotation reverses between them, across three seed populations: interior far enough from the edge that no trace can reach it -- pure stepping cost boundary in a band along the perimeter -- forces the out-of-domain exit branch mixed spread over the whole domain -- what the display actually does Separating them matters: the populations turn out to differ by two orders of magnitude per seed, and an undifferentiated average would have hidden that. Alongside wall time it reports an ExtractData call count, from a CXX_TEST-only counter incremented where the four per-step searches happen. Without it an optimization cannot be shown to have removed searches rather than merely found a faster machine. It also reports a setup breakdown -- BuildTriangles, the GmTriSearch R-tree build, and an activity-only reapply timed separately -- because which of those dominates decides whether triangulations can be shared across timesteps. XMGT_BENCH_SEEDS and XMGT_BENCH_CELLS size the run so a sweep needs no recompile. The defaults are small enough to leave the case in the regular suite, and the assertions are order-of-magnitude guards rather than tight bounds so it will not go flaky on a shared runner. One assertion is deliberately loose for a measured reason: a seed that exits the grid on its first step can reach the points.size() < 3 early return and come back with only the seed point, so "every seed yields a usable polyline" is false. It shows up at roughly 1 in 100,000. --- xmsgridtrace/gridtrace/XmGridTrace.cpp | 354 +++++++++++++++++++++++++ xmsgridtrace/gridtrace/XmGridTrace.t.h | 1 + 2 files changed, 355 insertions(+) diff --git a/xmsgridtrace/gridtrace/XmGridTrace.cpp b/xmsgridtrace/gridtrace/XmGridTrace.cpp index eb2ce0c..b4eb405 100644 --- a/xmsgridtrace/gridtrace/XmGridTrace.cpp +++ b/xmsgridtrace/gridtrace/XmGridTrace.cpp @@ -45,6 +45,20 @@ namespace { /// XMS Namespace +#ifdef CXX_TEST +/// \brief Count of XmUGrid2dDataExtractor::ExtractData calls since it was last zeroed. +/// Test-build-only instrumentation for testTraceBenchmark. A trace's cost is dominated by +/// the point-location search each ExtractData performs, so the benchmark needs the search +/// count and not only wall time -- otherwise an algorithmic win cannot be told apart from +/// a faster machine. Not thread safe; the benchmark is single threaded. +size_t g_extractDataCalls = 0; +/// \brief Adds a_n to the ExtractData call count. Compiles away outside test builds. +#define XMGT_COUNT_EXTRACT_DATA(a_n) (g_extractDataCalls += (a_n)) +#else +/// \brief No-op outside test builds, so production traces pay nothing for instrumentation. +#define XMGT_COUNT_EXTRACT_DATA(a_n) ((void)0) +#endif + //----- Class / Function definitions ------------------------------------------- //////////////////////////////////////////////////////////////////////////////// @@ -528,6 +542,7 @@ bool XmGridTraceImpl::GetVectorAtLocationAndTime(const xms::Pt3d& a_pt, xms::VecFlt dataOuty1; m_extractor1x->ExtractData(dataOutx1); m_extractor1y->ExtractData(dataOuty1); + XMGT_COUNT_EXTRACT_DATA(2); if (dataOutx1.size() != 1 || dataOuty1.size() != 1) { XM_LOG(xmlog::error, "Gridtracer: An error occured when extracting data"); @@ -540,6 +555,7 @@ bool XmGridTraceImpl::GetVectorAtLocationAndTime(const xms::Pt3d& a_pt, xms::VecFlt dataOuty2; m_extractor2x->ExtractData(dataOutx2); m_extractor2y->ExtractData(dataOuty2); + XMGT_COUNT_EXTRACT_DATA(2); if (dataOutx2.size() != 1 || dataOuty2.size() != 1) { XM_LOG(xmlog::error, "Gridtracer: An error occured when extracting data"); @@ -589,7 +605,14 @@ BSHP XmGridTrace::New(std::shared_ptr a_ugrid) #ifdef CXX_TEST #include +#include +#include +#include +#include +#include + #include +#include #include using namespace xms; @@ -708,6 +731,204 @@ void iCreateDefaultTwoCell(BSHP& a_tracer) a_tracer->AddGridScalarsAtTime(scalars, DataLocationEnum::LOC_CELLS, pointActivity, DataLocationEnum::LOC_CELLS, time); } // iCreateDefaultTwoCell + +//------------------------------------------------------------------------------ +/// \brief A structured quad grid plus its point locations, for the tracing benchmark. +/// The locations are kept alongside the ugrid so the velocity field can be evaluated +/// without depending on how the ugrid exposes its points. +//------------------------------------------------------------------------------ +struct BenchmarkGrid +{ + std::shared_ptr m_ugrid; ///< the grid itself + VecPt3d m_points; ///< grid point locations, in grid point order +}; + +//------------------------------------------------------------------------------ +/// \brief Measurements from one benchmark batch. +//------------------------------------------------------------------------------ +struct BenchmarkStats +{ + int m_seeds = 0; ///< seed points handed to TracePoint + int m_traced = 0; ///< seeds that produced a usable (2+ point) polyline + size_t m_tracePoints = 0; ///< total polyline points produced + size_t m_extractCalls = 0; ///< XmUGrid2dDataExtractor::ExtractData calls consumed + double m_seconds = 0; ///< wall time of the traced batch, excluding setup + std::map m_exitReasons; ///< exit message -> count, over a sample +}; + +//------------------------------------------------------------------------------ +/// \brief Builds a structured quad grid standing in for a real hydrodynamic mesh. +/// \param[in] a_cellsPerSide Number of cells along each axis +/// \param[in] a_length Length of the square domain along each axis +/// \return the grid and its point locations +//------------------------------------------------------------------------------ +BenchmarkGrid iBuildBenchmarkGrid(int a_cellsPerSide, double a_length) +{ + const int ptsPerSide = a_cellsPerSide + 1; + const double dx = a_length / a_cellsPerSide; + BenchmarkGrid grid; + grid.m_points.reserve((size_t)ptsPerSide * ptsPerSide); + for (int j = 0; j < ptsPerSide; ++j) + { + for (int i = 0; i < ptsPerSide; ++i) + grid.m_points.push_back({i * dx, j * dx, 0.0}); + } + + VecInt cells; + cells.reserve((size_t)a_cellsPerSide * a_cellsPerSide * 6); + for (int j = 0; j < a_cellsPerSide; ++j) + { + for (int i = 0; i < a_cellsPerSide; ++i) + { + const int p0 = j * ptsPerSide + i; + cells.push_back(XMU_QUAD); + cells.push_back(4); + cells.push_back(p0); + cells.push_back(p0 + 1); + cells.push_back(p0 + ptsPerSide + 1); + cells.push_back(p0 + ptsPerSide); + } + } + grid.m_ugrid = XmUGrid::New(grid.m_points, cells); + return grid; +} // iBuildBenchmarkGrid +//------------------------------------------------------------------------------ +/// \brief Builds a rotating-plus-drifting velocity field over the grid points. +/// A vortex is used rather than a uniform field for two reasons: the curvature makes the +/// adaptive stepping subdivide the way it does on real flow, and the drift carries part +/// of the seed population off the grid so the out-of-domain exit path -- which builds a +/// fresh polyline extractor per event -- is measured rather than assumed away. +/// \param[in] a_points Grid point locations +/// \param[in] a_omega Angular rate of the vortex; negative reverses the rotation +/// \param[in] a_drift Uniform velocity added in +x +/// \param[in] a_length Length of the square domain along each axis +/// \return velocity vectors, one per grid point +//------------------------------------------------------------------------------ +VecPt3d iBenchmarkVectors(const VecPt3d& a_points, double a_omega, double a_drift, double a_length) +{ + const double cx = a_length / 2, cy = a_length / 2; + VecPt3d vectors; + vectors.reserve(a_points.size()); + for (const auto& pt : a_points) + vectors.push_back({-a_omega * (pt.y - cy) + a_drift, a_omega * (pt.x - cx), 0.0}); + return vectors; +} // iBenchmarkVectors +//------------------------------------------------------------------------------ +/// \brief Builds seed points scattered inside a rectangular band of the domain. +/// The scatter is driven by a fixed linear congruential generator rather than std::rand +/// so that reruns and different machines trace the identical seed set; a benchmark whose +/// input changes between runs cannot measure a delta. +/// \param[in] a_count Number of seeds +/// \param[in] a_lo Low corner of the band, on both axes +/// \param[in] a_hi High corner of the band, on both axes +/// \param[in] a_holeLo Low corner of a rectangular hole to reject seeds from +/// \param[in] a_holeHi High corner of the hole; pass a_holeHi <= a_holeLo for no hole +/// \return the seed points +//------------------------------------------------------------------------------ +VecPt3d iBenchmarkSeeds(int a_count, double a_lo, double a_hi, double a_holeLo, double a_holeHi) +{ + unsigned int state = 12345u; + auto nextUnit = [&state]() { + state = state * 1664525u + 1013904223u; + return (state >> 8) / 16777216.0; + }; + + VecPt3d seeds; + seeds.reserve(a_count); + while ((int)seeds.size() < a_count) + { + const double x = a_lo + nextUnit() * (a_hi - a_lo); + const double y = a_lo + nextUnit() * (a_hi - a_lo); + const bool inHole = + a_holeHi > a_holeLo && x > a_holeLo && x < a_holeHi && y > a_holeLo && y < a_holeHi; + if (!inHole) + seeds.push_back({x, y, 0.0}); + } + return seeds; +} // iBenchmarkSeeds +//------------------------------------------------------------------------------ +/// \brief Traces every seed and measures the batch. +/// Timing covers only the TracePoint calls. The exit-reason histogram is gathered in a +/// separate untimed pass over a sample, because GetExitMessage returns a std::string by +/// value and a per-seed map insert would show up in a measurement this small. +/// \param[in] a_tracer The tracer, already loaded with two time steps +/// \param[in] a_seeds The seed points +/// \param[out] a_stats The measurements +//------------------------------------------------------------------------------ +void iRunTraceBenchmark(BSHP& a_tracer, + const VecPt3d& a_seeds, + BenchmarkStats& a_stats) +{ + a_stats = BenchmarkStats(); + a_stats.m_seeds = (int)a_seeds.size(); + + VecPt3d trace; + VecDbl times; + g_extractDataCalls = 0; + const auto start = std::chrono::steady_clock::now(); + for (const auto& seed : a_seeds) + { + a_tracer->TracePoint(seed, 0.0, trace, times); + if (trace.size() > 1) + { + ++a_stats.m_traced; + a_stats.m_tracePoints += trace.size(); + } + } + const auto end = std::chrono::steady_clock::now(); + a_stats.m_seconds = std::chrono::duration(end - start).count(); + a_stats.m_extractCalls = g_extractDataCalls; + + const int sampleSize = std::min((int)a_seeds.size(), 1000); + for (int i = 0; i < sampleSize; ++i) + { + a_tracer->TracePoint(a_seeds[i], 0.0, trace, times); + a_stats.m_exitReasons[a_tracer->GetExitMessage()]++; + } +} // iRunTraceBenchmark +//------------------------------------------------------------------------------ +/// \brief Prints one benchmark batch in a form that can be pasted into a results table. +/// \param[in] a_label Which seed population this batch was +/// \param[in] a_stats The measurements +//------------------------------------------------------------------------------ +void iReportTraceBenchmark(const char* a_label, const BenchmarkStats& a_stats) +{ + const double seeds = a_stats.m_seeds ? (double)a_stats.m_seeds : 1.0; + const double usPerSeed = a_stats.m_seconds * 1e6 / seeds; + const double extractsPerSeed = a_stats.m_extractCalls / seeds; + const double usPerExtract = + a_stats.m_extractCalls ? a_stats.m_seconds * 1e6 / a_stats.m_extractCalls : 0.0; + const double ptsPerTrace = + a_stats.m_traced ? (double)a_stats.m_tracePoints / a_stats.m_traced : 0.0; + + std::cout << std::fixed << std::setprecision(3) << "\n [" << a_label + << "] seeds=" << a_stats.m_seeds << " traced=" << a_stats.m_traced << "\n" + << " wall " << a_stats.m_seconds * 1e3 << " ms\n" + << " per seed " << usPerSeed << " us\n" + << " ExtractData " << a_stats.m_extractCalls << " calls (" + << std::setprecision(1) << extractsPerSeed << "/seed, " << std::setprecision(3) + << usPerExtract << " us/call)\n" + << " trace points " << a_stats.m_tracePoints << " (" << std::setprecision(1) + << ptsPerTrace << "/trace)\n" + << " exit reasons (sampled):\n"; + for (const auto& reason : a_stats.m_exitReasons) + std::cout << " " << std::setw(5) << reason.second << " " << reason.first << "\n"; + std::cout << std::flush; +} // iReportTraceBenchmark +//------------------------------------------------------------------------------ +/// \brief Reads a positive integer from the environment, or returns a fallback. +/// \param[in] a_name Environment variable name +/// \param[in] a_fallback Value to use when unset, unparseable, or not positive +/// \return the resolved value +//------------------------------------------------------------------------------ +int iEnvInt(const char* a_name, int a_fallback) +{ + const char* raw = std::getenv(a_name); + if (!raw) + return a_fallback; + const int value = std::atoi(raw); + return value > 0 ? value : a_fallback; +} // iEnvInt } //////////////////////////////////////////////////////////////////////////////// /// \class XmGridTraceUnitTests @@ -1491,5 +1712,138 @@ void XmGridTraceUnitTests::testTutorial() TS_ASSERT_DELTA_VEC(expectedOutTimes, outTimes, .0001); } // XmGridTraceUnitTests::testTutorial //! [snip_test_Example_XmGridTrace] +//------------------------------------------------------------------------------ +/// \brief Measures the cost of tracing many seed points over a realistic grid. +/// +/// This is the baseline for routing the "follow flow path" vector display option through +/// XmGridTrace: the display traces every visible glyph, so the number that matters is the +/// per-seed cost at glyph counts, not the cost of one trace. Three seed populations are +/// measured separately because they exercise different code: +/// +/// interior seeds far enough from the edge that no trace can reach it -- the pure +/// stepping cost, four ExtractData searches per integration step +/// boundary seeds in a band along the edge, so traces run out of the domain and pay +/// for a freshly constructed XmUGrid2dPolylineDataExtractor and +/// GmMultiPolyIntersector per exit event, inside the stepping loop +/// mixed seeds spread over the whole domain -- what the display actually does +/// +/// Reported alongside wall time is the ExtractData call count, so a later optimization +/// can be shown to have removed searches rather than merely found a faster machine. +/// +/// Seed count and grid size come from XMGT_BENCH_SEEDS and XMGT_BENCH_CELLS so a sweep +/// needs no recompile; the defaults are small enough to leave in the regular suite. The +/// assertions are deliberately loose -- this guards against order-of-magnitude +/// regressions, and a tight bound would only make the suite flaky on shared runners. +//------------------------------------------------------------------------------ +void XmGridTraceUnitTests::testTraceBenchmark() +{ + const int seedCount = iEnvInt("XMGT_BENCH_SEEDS", 250); + const int cellsPerSide = iEnvInt("XMGT_BENCH_CELLS", 200); + const double length = 200.0; + const double omega = 0.05; // vortex rate; reversed at the second time step + const double drift = 1.0; // uniform +x velocity, carries seeds off the +x edge + const double timeStepInterval = 10.0; + const double maxTracingDistance = 15.0; + + const auto setupStart = std::chrono::steady_clock::now(); + BenchmarkGrid grid = iBuildBenchmarkGrid(cellsPerSide, length); + const auto gridBuilt = std::chrono::steady_clock::now(); + + BSHP tracer = XmGridTrace::New(grid.m_ugrid); + tracer->SetVectorMultiplier(1); + tracer->SetMaxTracingTime(timeStepInterval); + tracer->SetMaxTracingDistance(maxTracingDistance); + tracer->SetMinDeltaTime(.01); + tracer->SetMaxChangeDistance(2.0); + tracer->SetMaxChangeVelocity(-1); + tracer->SetMaxChangeDirectionInRadians(0.2); + + DynBitset pointActivity; + pointActivity.resize(grid.m_points.size(), true); + // The rotation reverses between the two steps, so a trace that spans them is genuinely + // time dependent -- a single-timestep tracer cannot reproduce its path. + VecPt3d vectors1 = iBenchmarkVectors(grid.m_points, omega, drift, length); + VecPt3d vectors2 = iBenchmarkVectors(grid.m_points, -omega, drift, length); + tracer->AddGridScalarsAtTime(vectors1, DataLocationEnum::LOC_POINTS, pointActivity, + DataLocationEnum::LOC_POINTS, 0.0); + tracer->AddGridScalarsAtTime(vectors2, DataLocationEnum::LOC_POINTS, pointActivity, + DataLocationEnum::LOC_POINTS, timeStepInterval); + const auto setupEnd = std::chrono::steady_clock::now(); + + const double gridSeconds = std::chrono::duration(gridBuilt - setupStart).count(); + const double scalarSeconds = std::chrono::duration(setupEnd - gridBuilt).count(); + + // Break the per-timestep setup cost into its parts. This decides whether two timesteps + // with *different* cell activity can share one triangulation: activity is not baked into + // the triangles, it is latched onto the search object (XmUGridTriangles2d.cpp:146-164), + // so the question is whether flipping it per query is cheaper than triangulating twice. + DynBitset benchActivity; + benchActivity.resize(grid.m_ugrid->GetCellCount(), true); + + BSHP tris = XmUGridTriangles2d::New(); + const auto triStart = std::chrono::steady_clock::now(); + tris->BuildTriangles(*grid.m_ugrid, XmUGridTriangles2d::PO_CENTROIDS_ONLY); + const auto triBuilt = std::chrono::steady_clock::now(); + tris->SetCellActivity(benchActivity); // first call also builds the GmTriSearch R-tree + const auto searchBuilt = std::chrono::steady_clock::now(); + tris->SetCellActivity(benchActivity); // second call is the activity mask alone + const auto activityFlipped = std::chrono::steady_clock::now(); + + const double triSeconds = std::chrono::duration(triBuilt - triStart).count(); + const double searchSeconds = std::chrono::duration(searchBuilt - triBuilt).count(); + const double flipSeconds = std::chrono::duration(activityFlipped - searchBuilt).count(); + + std::cout << std::fixed << std::setprecision(3) << "\n=== XmGridTrace trace benchmark ===\n" + << " grid " << cellsPerSide << "x" << cellsPerSide << " quads, " + << grid.m_points.size() << " points\n" + << " seeds per set " << seedCount << "\n" + << " grid build " << gridSeconds * 1e3 << " ms\n" + << " add 2 timesteps " << scalarSeconds * 1e3 << " ms\n" + << " setup breakdown, one XmUGridTriangles2d:\n" + << " BuildTriangles " << triSeconds * 1e3 << " ms\n" + << " + R-tree & activity " << searchSeconds * 1e3 << " ms\n" + << " activity flip only " << flipSeconds * 1e3 << " ms\n" + << std::flush; + + // No trace can travel maxTracingDistance from this band, so nothing exits the grid. + const double interiorMargin = maxTracingDistance + 5.0; + VecPt3d interiorSeeds = + iBenchmarkSeeds(seedCount, interiorMargin, length - interiorMargin, 0.0, 0.0); + // Seeds within a band of the edge; the hole rejects anything that is not in the band. + const double boundaryBand = 5.0; + VecPt3d boundarySeeds = + iBenchmarkSeeds(seedCount, 0.5, length - 0.5, boundaryBand, length - boundaryBand); + VecPt3d mixedSeeds = iBenchmarkSeeds(seedCount, 0.5, length - 0.5, 0.0, 0.0); + + BenchmarkStats interior, boundary, mixed; + iRunTraceBenchmark(tracer, interiorSeeds, interior); + iReportTraceBenchmark("interior", interior); + iRunTraceBenchmark(tracer, boundarySeeds, boundary); + iReportTraceBenchmark("boundary", boundary); + iRunTraceBenchmark(tracer, mixedSeeds, mixed); + iReportTraceBenchmark("mixed", mixed); + + // Interior seeds cannot reach a boundary, so every one of them must trace. + TS_ASSERT_EQUALS(interior.m_traced, seedCount); + // Seeds that can leave the grid are not guaranteed a usable polyline: a seed that exits + // on its first step can hit the "failed to find an intersection when exiting grid" early + // return (:404-408) and come back holding only the seed point. Measured at roughly 1 in + // 100,000, so allow a small tail rather than asserting a false invariant -- but keep the + // bound tight enough that a real breakage in tracing still fails here. + TS_ASSERT(mixed.m_traced >= seedCount - 1 - seedCount / 1000); + // The instrumentation itself has to be working, or the search counts mean nothing. + TS_ASSERT(interior.m_extractCalls > (size_t)seedCount); + // The boundary set must actually leave the grid, otherwise this benchmark silently + // stops measuring the per-exit extractor construction it exists to measure. + const std::string outOfDomain = "Point has traveled out of domain."; + TS_ASSERT(boundary.m_exitReasons.count(outOfDomain) > 0); + TS_ASSERT_EQUALS(interior.m_exitReasons.count(outOfDomain), 0); + // Re-latching activity onto an existing search must stay cheaper than rebuilding the + // triangulation, or "share one triangulation and flip activity" is not even a candidate. + TS_ASSERT(flipSeconds < triSeconds); + // Order-of-magnitude guard only. Measured at ~0.1 ms/seed; 10 ms leaves room for a + // debug build on a loaded machine while still catching a real algorithmic regression. + TS_ASSERT(mixed.m_seconds * 1e3 / seedCount < 10.0); +} // XmGridTraceUnitTests::testTraceBenchmark #endif \ No newline at end of file diff --git a/xmsgridtrace/gridtrace/XmGridTrace.t.h b/xmsgridtrace/gridtrace/XmGridTrace.t.h index 7f6dc78..a464322 100644 --- a/xmsgridtrace/gridtrace/XmGridTrace.t.h +++ b/xmsgridtrace/gridtrace/XmGridTrace.t.h @@ -38,6 +38,7 @@ class XmGridTraceUnitTests : public CxxTest::TestSuite void testInactiveCell(); void testStartInactiveCell(); void testTutorial(); + void testTraceBenchmark(); }; // XmGridTraceUnitTests From fc72c1931cef6330508feb94220f7462f6a9b6fe Mon Sep 17 00:00:00 2001 From: Bill Dolinar Date: Fri, 14 Aug 2026 08:28:23 -0600 Subject: [PATCH 02/10] Cache the boundary-exit polyline extractor instead of rebuilding it per exit TracePoint constructed a fresh XmUGrid2dPolylineDataExtractor on every out-of-domain step, inside the stepping loop. That constructor triangulates the whole grid, and the SetPolyline that follows indexes every triangle into a new GmMultiPolyIntersector. Both depend only on the grid, which cannot change during a trace, and both were thrown away at the end of the if block and rebuilt for the next exiting particle. Benchmarked on a 200x200 grid, one exit event cost ~40 ms against ~48 us for a complete interior trace -- roughly 830x. On a seed population with a 5% exit rate, those 5% accounted for 98% of total trace time. It is now a member built lazily on the first exit and reused. Three things make the reuse safe, each checked in xmsextractor rather than assumed: - BuildTriangles is guarded by m_triangleType != a_location, so the second SetPolyline skips the triangulation. - ComputeExtractLocations builds m_multiPolyIntersector only when null and clears its output locations at entry, so no state carries between polylines. - XmGridTrace consumes only GetExtractLocations(), never extracted values, so the dummy zero scalars the constructor installs are irrelevant and the instance stays valid for the tracer's lifetime. m_ugrid is fixed at construction with no setter, so a cached extractor cannot outlive the grid it was built for. The member stays null until a trace actually exits, so a seed population that never reaches a boundary pays no memory for it. Measured A/B on one machine state at 1,000 seeds: the boundary population goes from 15,996 to 71 us/seed (224x) and the realistic mixed population from 2,200 to 51 us/seed (43x), while the interior population -- which never exits, and is the control -- moves 1%. Per exit the cost falls from 40 ms to ~9.5 us. ExtractData counts are identical before and after: this removes no searches at all, only whole-grid rebuilds, so the separate search-reduction work is still available on top of it. testBoundaryExtractorIsCached is the guard. Caching is invisible in the output, so the assertion that matters is a construction count from a CXX_TEST-only counter; the test also compares the two traces to 1e-12, which is what would catch reuse silently changing an answer. --- xmsgridtrace/gridtrace/XmGridTrace.cpp | 77 +++++++++++++++++++++++--- xmsgridtrace/gridtrace/XmGridTrace.t.h | 1 + 2 files changed, 70 insertions(+), 8 deletions(-) diff --git a/xmsgridtrace/gridtrace/XmGridTrace.cpp b/xmsgridtrace/gridtrace/XmGridTrace.cpp index b4eb405..54454fe 100644 --- a/xmsgridtrace/gridtrace/XmGridTrace.cpp +++ b/xmsgridtrace/gridtrace/XmGridTrace.cpp @@ -54,9 +54,18 @@ namespace size_t g_extractDataCalls = 0; /// \brief Adds a_n to the ExtractData call count. Compiles away outside test builds. #define XMGT_COUNT_EXTRACT_DATA(a_n) (g_extractDataCalls += (a_n)) +/// \brief Count of XmUGrid2dPolylineDataExtractor constructions since it was last zeroed. +/// Test-build-only instrumentation for testBoundaryExtractorIsCached. Caching that extractor +/// is a pure performance change with no effect on trace output, so a construction count is +/// the only thing that can tell a cached run from an uncached one. +size_t g_boundaryExtractorBuilds = 0; +/// \brief Records one boundary-extractor construction. Compiles away outside test builds. +#define XMGT_COUNT_BOUNDARY_EXTRACTOR_BUILD() (++g_boundaryExtractorBuilds) #else /// \brief No-op outside test builds, so production traces pay nothing for instrumentation. #define XMGT_COUNT_EXTRACT_DATA(a_n) ((void)0) +/// \brief No-op outside test builds. +#define XMGT_COUNT_BOUNDARY_EXTRACTOR_BUILD() ((void)0) #endif //----- Class / Function definitions ------------------------------------------- @@ -127,6 +136,13 @@ class XmGridTraceImpl : public XmGridTrace /// data extractor for the y component for the second time step BSHP m_extractor2y; double m_time2=-1; ///< time of the second time step + /// Extractor used to find where a trace leaves the grid, built lazily on the first + /// out-of-domain step and reused for every one after it. Its construction triangulates the + /// whole grid and its first SetPolyline indexes every triangle into a GmMultiPolyIntersector; + /// neither depends on the polyline, and both were previously rebuilt per exit event at a + /// measured ~40 ms each. Null until a trace actually exits, so a tracer whose traces all + /// stay inside the grid never pays the memory. + BSHP m_boundaryExtractor; double m_distTraveled=0; ///< distance traveled in the last TracePoint operation std::string m_exitMessage; ///< exit message for the last TracePoint operation @@ -410,11 +426,17 @@ void XmGridTraceImpl::TracePoint(const Pt3d& a_pt, { m_exitMessage = "Point has traveled out of domain."; VecPt3d points = {pt0, pt1}; - // DataLocationEnum is irrelevant here. - BSHP polylineExtractor = - XmUGrid2dPolylineDataExtractor::New(m_ugrid, DataLocationEnum::LOC_POINTS); - polylineExtractor->SetPolyline(points); - points = polylineExtractor->GetExtractLocations(); + if (!m_boundaryExtractor) + { + // DataLocationEnum is irrelevant here: only the extract locations are consumed below, + // never the extracted values, so the dummy zero scalars the constructor installs do + // not matter and the instance stays valid for this tracer's lifetime. + m_boundaryExtractor = + XmUGrid2dPolylineDataExtractor::New(m_ugrid, DataLocationEnum::LOC_POINTS); + XMGT_COUNT_BOUNDARY_EXTRACTOR_BUILD(); + } + m_boundaryExtractor->SetPolyline(points); + points = m_boundaryExtractor->GetExtractLocations(); if (points.size() < 3) { XM_LOG(xmlog::error, "Gridtracer failed to find an intersection when exiting grid."); @@ -1713,6 +1735,44 @@ void XmGridTraceUnitTests::testTutorial() } // XmGridTraceUnitTests::testTutorial //! [snip_test_Example_XmGridTrace] //------------------------------------------------------------------------------ +/// \brief Verifies the boundary-exit extractor is built once per tracer, not once per exit. +/// +/// The extractor's constructor triangulates the whole grid and its first SetPolyline indexes +/// every triangle into a GmMultiPolyIntersector -- both grid-only work, and both measured at +/// ~40 ms per exit event when rebuilt inside the stepping loop. Caching it changes no output, +/// so the construction count is what has to be asserted; the trace comparison is here to +/// catch the reuse silently changing an answer. +//------------------------------------------------------------------------------ +void XmGridTraceUnitTests::testBoundaryExtractorIsCached() +{ + BSHP tracer; + iCreateDefaultSingleCell(tracer); + + // The default single cell has a uniform (1, 1) field, so a trace from the middle leaves the + // grid on its first step. + const Pt3d startPoint = {.5, .5, 0}; + const double startTime = .5; + const std::string outOfDomain = "Point has traveled out of domain."; + + g_boundaryExtractorBuilds = 0; + + VecPt3d firstTrace; + VecDbl firstTimes; + tracer->TracePoint(startPoint, startTime, firstTrace, firstTimes); + TS_ASSERT_EQUALS(outOfDomain, tracer->GetExitMessage()); + TS_ASSERT_EQUALS(size_t(1), g_boundaryExtractorBuilds); + TS_ASSERT(firstTrace.size() >= 2); + + VecPt3d secondTrace; + VecDbl secondTimes; + tracer->TracePoint(startPoint, startTime, secondTrace, secondTimes); + TS_ASSERT_EQUALS(outOfDomain, tracer->GetExitMessage()); + TS_ASSERT_EQUALS(size_t(1), g_boundaryExtractorBuilds); + + TS_ASSERT_DELTA_VECPT3D(firstTrace, secondTrace, 1e-12); + TS_ASSERT_DELTA_VEC(firstTimes, secondTimes, 1e-12); +} // XmGridTraceUnitTests::testBoundaryExtractorIsCached +//------------------------------------------------------------------------------ /// \brief Measures the cost of tracing many seed points over a realistic grid. /// /// This is the baseline for routing the "follow flow path" vector display option through @@ -1722,9 +1782,10 @@ void XmGridTraceUnitTests::testTutorial() /// /// interior seeds far enough from the edge that no trace can reach it -- the pure /// stepping cost, four ExtractData searches per integration step -/// boundary seeds in a band along the edge, so traces run out of the domain and pay -/// for a freshly constructed XmUGrid2dPolylineDataExtractor and -/// GmMultiPolyIntersector per exit event, inside the stepping loop +/// boundary seeds in a band along the edge, so traces run out of the domain and pay for +/// the XmUGrid2dPolylineDataExtractor path -- a whole-grid triangulation plus a +/// GmMultiPolyIntersector, once per tracer since that extractor is cached +/// (it was once per exit event, inside the stepping loop) /// mixed seeds spread over the whole domain -- what the display actually does /// /// Reported alongside wall time is the ExtractData call count, so a later optimization diff --git a/xmsgridtrace/gridtrace/XmGridTrace.t.h b/xmsgridtrace/gridtrace/XmGridTrace.t.h index a464322..faa0086 100644 --- a/xmsgridtrace/gridtrace/XmGridTrace.t.h +++ b/xmsgridtrace/gridtrace/XmGridTrace.t.h @@ -38,6 +38,7 @@ class XmGridTraceUnitTests : public CxxTest::TestSuite void testInactiveCell(); void testStartInactiveCell(); void testTutorial(); + void testBoundaryExtractorIsCached(); void testTraceBenchmark(); }; // XmGridTraceUnitTests From 9b3efdb8d300a054771272503c234fe484b6e399 Mon Sep 17 00:00:00 2001 From: Bill Dolinar Date: Fri, 14 Aug 2026 08:35:29 -0600 Subject: [PATCH 03/10] Fix the inverted time interpolation and propagate no-data through it GetVectorAtLocationAndTime weighted each timestep by its own distance from the current time, so perc1 was |t - t1| / totalTime and multiplied timestep 1's value. At t == t1 that weight is zero: a particle released exactly at the first timestep was advected entirely by the field at the second. A timestep is weighted by its *closeness* to the current time, so the distance from one timestep is the weight of the other. This is the defect that makes the tracer worth using at all. Routing the "follow flow path" display option through it is only an improvement if the trace follows the field as it changes, and an inverted blend does not. Propagating XM_NODATA is part of the same fix rather than a separate one. With the weights corrected, a location that is inactive in only one of the two timesteps stops resolving to the sentinel and starts resolving to a blend of it: 0.9 * 0.1 + 0.1 * -9999999 is -999999.9, which is neither no-data nor a velocity, and which passes every no-data test the callers make. The sentinel is now propagated when either bracketing timestep has no data at the location, which is what makes "a cell active at t1 but inactive at t2 terminates the trace" actually hold. testStartInactiveCell was passing before only because the inverted weights happened to give timestep 2 all the weight at t == t1; with the propagation it passes for the right reason. testTimeVaryingFieldChangesPath is the regression guard. One cell spanning the domain gives a spatially uniform field, so any curvature in the path can only have come from time; the field rotates +x -> +y between timesteps rather than reversing, so the interpolated velocity never passes through zero and cannot trip the velocity-is-zero exit partway along. Its first assertion is the one that catches an inversion -- a particle released at t1 must step due east with y untouched -- and it compares against a frozen-field control traced by the same code, which never turns. Three existing baselines move. Each first step was checked by hand rather than accepted from the runner: testUniqueTimeSteps 0.5 -> 0.6, the t1 cell value 0.1 over dt 1, where it was 0.7 from t2's 0.2. Its second step, 0.9*0.11 + 0.1*0.21 = 0.12, matches the recorded 0.744 to the digit. testInactiveCell same first step, and the trace now terminates exactly at x = 1, the boundary of the cell that is inactive at t2, rather than at 0.9979 -- it had been stopping just short via a max-change- velocity blow-up on a no-data-contaminated blend. testTutorial first step y 0.5 -> 1.5: corner scalars interpolate to (0, 0.5), times the multiplier of 2, over dt 1. The old 1.25 was a boundary-clipping artifact of using the doubled second-timestep field at t = 0 -- the tutorial's own comment says the second timestep is doubled to show an increase, so the trace should start at the first timestep's magnitude and speed up. --- xmsgridtrace/gridtrace/XmGridTrace.cpp | 241 +++++++++++++++++++------ xmsgridtrace/gridtrace/XmGridTrace.t.h | 1 + 2 files changed, 183 insertions(+), 59 deletions(-) diff --git a/xmsgridtrace/gridtrace/XmGridTrace.cpp b/xmsgridtrace/gridtrace/XmGridTrace.cpp index 54454fe..3d3a5f4 100644 --- a/xmsgridtrace/gridtrace/XmGridTrace.cpp +++ b/xmsgridtrace/gridtrace/XmGridTrace.cpp @@ -589,11 +589,29 @@ bool XmGridTraceImpl::GetVectorAtLocationAndTime(const xms::Pt3d& a_pt, XM_LOG(xmlog::warning, "Gridtracer: The given time is before the first time step."); a_currentTime = m_time1; } + // A location outside the grid or in an inactive cell in *either* bracketing timestep has no + // usable velocity, and the sentinel must be propagated rather than weighted: blending + // XM_NODATA (-9999999) against a real value produces something like -999999.9, which is + // neither no-data nor meaningful, and every caller tests for XM_NODATA exactly. Returning + // true is correct -- extraction succeeded, and no-data is the answer. + if (EQ_TOL(dataOutx1[0], XM_NODATA, 1) || EQ_TOL(dataOuty1[0], XM_NODATA, 1) || + EQ_TOL(dataOutx2[0], XM_NODATA, 1) || EQ_TOL(dataOuty2[0], XM_NODATA, 1)) + { + a_data.x = XM_NODATA; + a_data.y = XM_NODATA; + return true; + } + double totalTime = fabs(m_time1 - m_time2); - double perc1 = fabs(a_currentTime - m_time1) / totalTime; - double perc2 = fabs(a_currentTime - m_time2) / totalTime; - a_data.x = dataOutx1[0] * perc1 + dataOutx2[0] * perc2; - a_data.y = dataOuty1[0] * perc1 + dataOuty2[0] * perc2; + // Each timestep is weighted by its *closeness* to the current time, so the distance from + // one timestep is the weight of the other: at a_currentTime == m_time1 the field is + // entirely timestep 1's. Weighting each timestep by its own distance instead -- which is + // what this did until the weights were swapped -- inverts the interpolation, advecting a + // particle released at m_time1 entirely by the field at m_time2. + double weight1 = fabs(a_currentTime - m_time2) / totalTime; + double weight2 = fabs(a_currentTime - m_time1) / totalTime; + a_data.x = dataOutx1[0] * weight1 + dataOutx2[0] * weight2; + a_data.y = dataOuty1[0] * weight1 + dataOuty2[0] * weight2; return true; } // XmGridTraceImpl::GetVectorAtLocationAndTime } // namespace {} @@ -1539,18 +1557,20 @@ void XmGridTraceUnitTests::testUniqueTimeSteps() tracer->TracePoint(startPoint, startTime, outTrace, outTimes); - VecPt3d expectedOutTrace = {{.5, .5, 0}, - {0.70000000298023224, 0.50000000000000000, 0.00000000000000000}, - {0.95200000226497650, 0.50000000000000000, 0.00000000000000000}, - {1.2734079944372176, 0.50000000000000000, 0.00000000000000000}, - {1.6897536998434066, 0.50000000000000000, 0.00000000000000000}, - {2, .5, 0}}; + VecPt3d expectedOutTrace = {{0.5, 0.5, 0}, + {0.60000000149011612, 0.5, 0}, + {0.74400000184774395, 0.5, 0}, + {0.95481600679159162, 0.5, 0}, + {1.2691074101881981, 0.5, 0}, + {1.747260385068264, 0.5, 0}, + {2, 0.5, 0}}; VecDbl expectedOutTimes = {10, - 11.000000000000000, + 11, 12.199999999999999, 13.640000000000001, - 15.368000000000000, - 16.627525378316030}; + 15.368, + 17.441600000000001, + 18.362609001148471}; TS_ASSERT_DELTA_VECPT3D(expectedOutTrace, outTrace, .0001); TS_ASSERT_DELTA_VEC(expectedOutTimes, outTimes, .0001); } // XmGridTraceUnitTests::testUniqueTimeSteps @@ -1579,11 +1599,16 @@ void XmGridTraceUnitTests::testInactiveCell() tracer->TracePoint(startPoint, startTime, outTrace, outTimes); - VecPt3d expectedOutTrace = {{.5, .5, 0}, - {0.70000000298023224, 0.50000000000000000, 0.00000000000000000}, - {0.93040000677108770, 0.50000000000000000, 0.00000000000000000}, - {0.99788877571821222, 0.50000000000000000, 0.00000000000000000}}; - VecDbl expectedOutTimes = {10, 11.000000000000000, 12.199999999999999, 12.560000000000000}; + VecPt3d expectedOutTrace = {{0.5, 0.5, 0}, + {0.60000000149011612, 0.5, 0}, + {0.74280000120401379, 0.5, 0}, + {0.94575130454301826, 0.5, 0}, + {1, 0.5, 0}}; + VecDbl expectedOutTimes = {10, + 11, + 12.199999999999999, + 13.640000000000001, + 13.969279307058475}; TS_ASSERT_DELTA_VECPT3D(expectedOutTrace, outTrace, .0001); TS_ASSERT_DELTA_VEC(expectedOutTimes, outTimes, .0001); } // XmGridTraceUnitTests::testInactiveCell @@ -1689,52 +1714,150 @@ void XmGridTraceUnitTests::testTutorial() // std::cout << tracer->GetExitMessage(); // Expected values for this simulation - VecPt3d expectedOutTrace = {{0.50000000000000000, 0.50000000000000000, 0.00000000000000000}, - {0.50000000000000000, 1.2500000000000000, 0.00000000000000000}, - {0.54457812566426578, 1.3391562513285316, 0.00000000000000000}, - {0.61632493250262921, 1.4354984729093498, 0.00000000000000000}, - {0.72535406450374607, 1.5315533661126233, 0.00000000000000000}, - {0.88236797164001590, 1.6126801842666139, 0.00000000000000000}, - {0.98873181403598276, 1.6331015959080102, 0.00000000000000000}, - {1.0538503898747653, 1.6342606013582104, 0.00000000000000000}, - {1.1249433009705341, 1.5683006835455087, 0.00000000000000000}, - {1.1895097427498795, 1.3863448896225066, 0.00000000000000000}, - {1.2235242118635632, 1.0588590059131318, 0.00000000000000000}, - {1.2235242118635632, 0.90477286425654002, 0.00000000000000000}, - {1.2005336220528682, 0.85080764250970042, 0.00000000000000000}, - {1.1581790674742278, 0.79387770198395835, 0.00000000000000000}, - {1.0896874578697060, 0.74131697161132859, 0.00000000000000000}, - {0.98966250551038770, 0.70663752692174131, 0.00000000000000000}, - {0.95806149614159530, 0.71817980325332686, 0.00000000000000000}, - {0.92629620502521459, 0.77371504022050730, 0.00000000000000000}, - {0.90239412753251202, 0.88917318465162865, 0.00000000000000000}, - {0.89995172701803572, 1.0694875660697027, 0.00000000000000000}, - {0.91503139037776327, 1.0911992829869794, 0.00000000000000000}, - {0.93816744602651825, 1.1127546977629765, 0.00000000000000000}, - {0.97140028507849163, 1.1309789606067331, 0.00000000000000000}, - {0.99364912627842006, 1.1358370729524059, 0.00000000000000000}, - {1.0071524474802995, 1.1364684019706512, 0.00000000000000000}, - {1.0223447138862345, 1.1280655805979485, 0.00000000000000000}, - {1.0369737821057583, 1.0971462034407997, 0.00000000000000000}, - {1.0467397711865176, 1.0371377237101163, 0.00000000000000000}, - {1.0467397711865176, 0.96499504248441559, 0.00000000000000000}, - {1.0390576209755447, 0.95473758230148376, 0.00000000000000000}, - {1.0276444556154691, 0.94488898976070590, 0.00000000000000000}, - {1.0208791233912420, 0.94149540451099356, 0.00000000000000000}}; - VecDbl expectedOutTimes = { - 0.00000000000000000, 0.37500000000000000, 0.82499999999999996, 1.3649999999999998, - 2.0129999999999999, 2.7905999999999995, 3.2571599999999994, 3.5370959999999991, - 3.8730191999999990, 4.2761270399999987, 4.7598564479999981, 5.3403317375999979, - 6.0369020851199977, 6.8727865021439971, 7.8758478025727969, 9.0795213630873555, - 9.4406234312417237, 9.8739459130269651, 10.393932891169255, 11.017917264940003, - 11.766698513464901, 12.665236011694777, 13.743481009570628, 14.390428008296139, - 14.778596207531445, 15.244398046613812, 15.803360253512654, 16.474114901791264, - 17.279020479725595, 18.244907173246794, 19.403971205472232, 20.000000000000000}; + VecPt3d expectedOutTrace = {{0.5, 0.5, 0}, + {0.5, 1.5, 0}, + {0.62600000187754634, 1.6260000018775462, 0}, + {0.82611968728899965, 1.7455603212296962, 0}, + {0.97840008102011689, 1.7810753047635555, 0}, + {1.0280095840364933, 1.7824472100312621, 0}, + {1.0861189816907613, 1.7608732599310344, 0}, + {1.1492686295114336, 1.6802752810470523, 0}, + {1.2097920698566107, 1.5101408581884392, 0}, + {1.2515951471975522, 1.2181485463468757, 0}, + {1.2515951471975522, 0.84053651390559747, 0}, + {1.2181758214493843, 0.78780883088769804, 0}, + {1.1632869448015855, 0.73137186792498654, 0}, + {1.0771209832183524, 0.67899546053648097, 0}, + {1.0129487663521615, 0.66357815692798783, 0}, + {0.97169356095126669, 0.66199025753694563, 0}, + {0.92552080990281416, 0.70419149113367874, 0}, + {0.88530832700558759, 0.83950990950827409, 0}, + {0.87513974259796246, 1.0941588844381676, 0}, + {0.90077009637050098, 1.128146252166127, 0}, + {0.943692705404238, 1.1613833261644337, 0}, + {0.97709108330292604, 1.1730361561747586, 0}, + {0.99894959169213471, 1.1759300874982919, 0}, + {1.0124203987349505, 1.1760105163064269, 0}, + {1.0275428271398932, 1.1645289800266216, 0}, + {1.042848666622334, 1.1337546211004945, 0}, + {1.055142468614698, 1.0758075939238765, 0}, + {1.0585305184379035, 0.98540145004498747, 0}, + {1.0556233679912082, 0.97374570199926891, 0}, + {1.0492587242876892, 0.9602613226646981, 0}, + {1.0375007181419984, 0.94568649411103145, 0}, + {1.017827020259642, 0.93210280494582176, 0}, + {1.0175992759724071, 0.93204300863222744, 0}}; + VecDbl expectedOutTimes = {0, + 1, + 2.2000000000000002, + 3.6400000000000001, + 4.5040000000000004, + 4.7632000000000003, + 5.0742400000000005, + 5.4474880000000008, + 5.8953856000000009, + 6.432862720000001, + 7.0778352640000008, + 7.8518023168000006, + 8.7805627801600004, + 9.8950753361920007, + 10.563782869811201, + 10.96500738998272, + 11.446476814188543, + 12.024240123235531, + 12.717556094091917, + 13.54953525911958, + 14.547910257152775, + 15.146935255972693, + 15.506350255264643, + 15.721999254839814, + 15.980778054330019, + 16.291312613718265, + 16.663954084984159, + 17.111123850503233, + 17.647727569126122, + 18.291652031473589, + 19.064361386290546, + 19.991612612070895, + 20}; TS_ASSERT_DELTA_VECPT3D(expectedOutTrace, outTrace, .0001); TS_ASSERT_DELTA_VEC(expectedOutTimes, outTimes, .0001); } // XmGridTraceUnitTests::testTutorial //! [snip_test_Example_XmGridTrace] //------------------------------------------------------------------------------ +/// \brief A trace through a field that changes between timesteps follows neither timestep. +/// +/// This is the regression guard for the time interpolation, which is the whole reason this +/// tracer is worth routing a display option through: a tracer that samples one frozen +/// timestep would be no better than the render-time drifter it replaces. +/// +/// The field rotates from +x at the first timestep to +y at the second rather than +/// reversing, so the interpolated velocity never passes through zero and cannot trip the +/// "velocity has gone to zero" exit partway along. +/// +/// The first assertion is the one that catches an inverted interpolation: a particle +/// released exactly at the first timestep must be advected by that timestep's field alone, +/// so its first step is due east with y untouched. Weighting each timestep by its own +/// distance from the current time instead sends that first step due north. +//------------------------------------------------------------------------------ +void XmGridTraceUnitTests::testTimeVaryingFieldChangesPath() +{ + // One cell spanning the whole domain, so cell-located scalars give a spatially uniform + // field and any curvature in the path can only have come from time. + VecPt3d points = {{0, 0, 0}, {10, 0, 0}, {10, 10, 0}, {0, 10, 0}}; + VecInt cells = {XMU_QUAD, 4, 0, 1, 2, 3}; + std::shared_ptr ugrid = XmUGrid::New(points, cells); + + DynBitset activity; + activity.push_back(true); + + auto traceWithField = [&](const Pt3d& a_first, const Pt3d& a_second, VecPt3d& a_outTrace) { + BSHP tracer = XmGridTrace::New(ugrid); + tracer->SetVectorMultiplier(1); + tracer->SetMaxTracingTime(5); + tracer->SetMaxTracingDistance(100); + tracer->SetMinDeltaTime(.01); + tracer->SetMaxChangeDistance(.5); + tracer->SetMaxChangeVelocity(-1); + tracer->SetMaxChangeDirectionInRadians(XM_PI); // never subdivide on direction + VecPt3d first = {a_first}; + VecPt3d second = {a_second}; + tracer->AddGridScalarsAtTime(first, DataLocationEnum::LOC_CELLS, activity, + DataLocationEnum::LOC_CELLS, 0.0); + tracer->AddGridScalarsAtTime(second, DataLocationEnum::LOC_CELLS, activity, + DataLocationEnum::LOC_CELLS, 10.0); + VecDbl outTimes; + tracer->TracePoint({1, 1, 0}, 0.0, a_outTrace, outTimes); + TS_ASSERT_EQUALS(a_outTrace.size(), outTimes.size()); + }; + + const Pt3d startPoint = {1, 1, 0}; + + VecPt3d rotating; + traceWithField({1, 0, 0}, {0, 1, 0}, rotating); + + // The same field at both timesteps -- what a single-timestep tracer would produce. + VecPt3d frozen; + traceWithField({1, 0, 0}, {1, 0, 0}, frozen); + + TS_ASSERT(rotating.size() >= 3); + TS_ASSERT(frozen.size() >= 3); + + // Released at the first timestep, so the first step is that timestep's field alone. + TS_ASSERT_DELTA(startPoint.y, rotating[1].y, 1e-9); + TS_ASSERT(rotating[1].x > startPoint.x); + + // A frozen field never turns. + for (size_t i = 0; i < frozen.size(); ++i) + { + TS_ASSERT_DELTA(startPoint.y, frozen[i].y, 1e-9); + } + + // A changing one does, and that difference is the feature. + TS_ASSERT(rotating.back().y > startPoint.y + 0.1); + TS_ASSERT(rotating.back().x < frozen.back().x); +} // XmGridTraceUnitTests::testTimeVaryingFieldChangesPath +//------------------------------------------------------------------------------ /// \brief Verifies the boundary-exit extractor is built once per tracer, not once per exit. /// /// The extractor's constructor triangulates the whole grid and its first SetPolyline indexes diff --git a/xmsgridtrace/gridtrace/XmGridTrace.t.h b/xmsgridtrace/gridtrace/XmGridTrace.t.h index faa0086..97077cb 100644 --- a/xmsgridtrace/gridtrace/XmGridTrace.t.h +++ b/xmsgridtrace/gridtrace/XmGridTrace.t.h @@ -38,6 +38,7 @@ class XmGridTraceUnitTests : public CxxTest::TestSuite void testInactiveCell(); void testStartInactiveCell(); void testTutorial(); + void testTimeVaryingFieldChangesPath(); void testBoundaryExtractorIsCached(); void testTraceBenchmark(); From 1b0bde3ae6401bd1b4d8a8fd3503c94c8d762472 Mon Sep 17 00:00:00 2001 From: Bill Dolinar Date: Fri, 14 Aug 2026 08:39:01 -0600 Subject: [PATCH 04/10] Add a batch TracePoints entry point and keep its output arrays aligned The "follow flow path" display traces every visible vector glyph, tens of thousands of them per redraw, driven from Python. One TracePoint call per glyph pays a language boundary crossing per glyph for work that is identical across them, and it can only report why the *last* trace ended -- GetExitMessage describes a single operation, so a caller has no way to ask why glyph 4,000 stopped short. TracePoints takes all the seeds at once and returns a polyline, a time array, and an exit message per seed. It does not advance the time steps: every trace runs against whichever pair AddGridScalarsAtTime most recently supplied, and a caller wanting traces that span more of a series feeds the next step and traces again. Keeping that in the caller is deliberate -- the two-step window is instance state, so a batch that advanced it internally would have to carry per-seed continuation state, which is a different and larger design than this one. Mismatched input lengths return nothing rather than tracing the common prefix. A caller that supplied the wrong number of start times has a bug, and a partial result lets it go unnoticed. TracePoint's two output arrays could also come back different lengths, which this fixes because the batch documents them as parallel. The position push was conditional on the step actually moving while the time push was unconditional, so a step shorter than XM_ZERO_TOL left the times array one longer and silently misaligned every later pair -- undetectable to a caller zipping them. The time is now pushed only when the point is. No existing expectation moves, so no recorded trace contained such a step. testTracePointsMatchesSerialTracePoint compares the batch against serial TracePoint calls on an identical fixture rather than against a recorded baseline, which would drift with the tracer instead of pinning the equivalence. Its seeds cover the three shapes a caller has to handle: two traces that leave the grid, and a seed outside it that yields an empty trace rather than a polyline. --- xmsgridtrace/gridtrace/XmGridTrace.cpp | 129 ++++++++++++++++++++++--- xmsgridtrace/gridtrace/XmGridTrace.h | 27 ++++++ xmsgridtrace/gridtrace/XmGridTrace.t.h | 1 + 3 files changed, 142 insertions(+), 15 deletions(-) diff --git a/xmsgridtrace/gridtrace/XmGridTrace.cpp b/xmsgridtrace/gridtrace/XmGridTrace.cpp index 3d3a5f4..5e6035a 100644 --- a/xmsgridtrace/gridtrace/XmGridTrace.cpp +++ b/xmsgridtrace/gridtrace/XmGridTrace.cpp @@ -110,6 +110,12 @@ class XmGridTraceImpl : public XmGridTrace VecPt3d& a_outTrace, VecDbl& a_outTimes) final; + void TracePoints(const VecPt3d& a_pts, + const VecDbl& a_ptTimes, + std::vector& a_outTraces, + std::vector& a_outTimes, + VecStr& a_outExitMessages) final; + std::string GetExitMessage() final; private: @@ -522,31 +528,61 @@ void XmGridTraceImpl::TracePoint(const Pt3d& a_pt, return; } - // add new pt if not identical to last - int size = (int)a_outTrace.size(); - if (size > 0) - { - if (!EQ_TOL(pt1.x, a_outTrace.at(size - 1).x, XM_ZERO_TOL) || - !EQ_TOL(pt1.y, a_outTrace.at(size - 1).y, XM_ZERO_TOL)) - { - a_outTrace.push_back(pt1); - } - } - else - { - a_outTrace.push_back(pt1); - } + // add new pt if not identical to last -- and push its time only when the point is + // pushed. The time push used to be unconditional, so a step shorter than XM_ZERO_TOL + // left a_outTimes one longer than a_outTrace and silently misaligned every later + // pair, which a caller reading them as parallel arrays cannot detect. + const bool moved = a_outTrace.empty() || !EQ_TOL(pt1.x, a_outTrace.back().x, XM_ZERO_TOL) || + !EQ_TOL(pt1.y, a_outTrace.back().y, XM_ZERO_TOL); pt0 = pt1; elapsedTime += deltaT; vx0 = vx1; vy0 = vy1; deltaT *= 1.2; mag0 = mag1; - a_outTimes.push_back(a_ptTime + elapsedTime); + if (moved) + { + a_outTrace.push_back(pt1); + a_outTimes.push_back(a_ptTime + elapsedTime); + } } } // while () } // XmGridTraceImpl::TracePoint //------------------------------------------------------------------------------ +/// \brief Runs the Grid Trace for many points against the current two time steps +/// \param[in] a_pts The starting point of each trace +/// \param[in] a_ptTimes The starting time of each trace; must be one per point +/// \param[out] a_outTraces The resultant positions at each step, one entry per point +/// \param[out] a_outTimes The resultant times, one entry per point +/// \param[out] a_outExitMessages What ended each trace, one entry per point +//------------------------------------------------------------------------------ +void XmGridTraceImpl::TracePoints(const VecPt3d& a_pts, + const VecDbl& a_ptTimes, + std::vector& a_outTraces, + std::vector& a_outTimes, + VecStr& a_outExitMessages) +{ + a_outTraces.clear(); + a_outTimes.clear(); + a_outExitMessages.clear(); + if (a_pts.size() != a_ptTimes.size()) + { + // Returning empty rather than tracing the common prefix: a caller that mismatched these + // has a bug, and a partial result would let it go unnoticed. + XM_LOG(xmlog::error, "Gridtracer: TracePoints needs one start time per point."); + return; + } + + a_outTraces.resize(a_pts.size()); + a_outTimes.resize(a_pts.size()); + a_outExitMessages.resize(a_pts.size()); + for (size_t i = 0; i < a_pts.size(); ++i) + { + TracePoint(a_pts[i], a_ptTimes[i], a_outTraces[i], a_outTimes[i]); + a_outExitMessages[i] = m_exitMessage; + } +} // XmGridTraceImpl::TracePoints +//------------------------------------------------------------------------------ /// \brief Returns the velocity scalar for a given point and time /// \param[in] a_pt The point /// \param[in] a_currentTime The time at extraction @@ -1858,6 +1894,69 @@ void XmGridTraceUnitTests::testTimeVaryingFieldChangesPath() TS_ASSERT(rotating.back().x < frozen.back().x); } // XmGridTraceUnitTests::testTimeVaryingFieldChangesPath //------------------------------------------------------------------------------ +/// \brief The batch entry point returns exactly what serial TracePoint calls return. +/// +/// TracePoints exists to cross a language boundary once instead of once per seed, so its +/// value depends entirely on it being a faithful stand-in. Comparing against serial +/// TracePoint on an identical fixture is the strongest oracle available -- stronger than a +/// recorded baseline, which would drift with the tracer rather than pin the equivalence. +/// +/// The seeds are chosen to cover the three shapes a caller has to handle: a trace that +/// leaves the grid, another that leaves it from elsewhere, and a seed outside the grid +/// entirely, which yields an empty trace rather than a polyline. +//------------------------------------------------------------------------------ +void XmGridTraceUnitTests::testTracePointsMatchesSerialTracePoint() +{ + const VecPt3d seeds = {{.5, .5, 0}, {.25, .75, 0}, {-.1, 0, 0}}; + const VecDbl seedTimes = {.5, .5, .5}; + + BSHP serialTracer; + iCreateDefaultSingleCell(serialTracer); + std::vector serialTraces(seeds.size()); + std::vector serialTimes(seeds.size()); + VecStr serialMessages(seeds.size()); + for (size_t i = 0; i < seeds.size(); ++i) + { + serialTracer->TracePoint(seeds[i], seedTimes[i], serialTraces[i], serialTimes[i]); + serialMessages[i] = serialTracer->GetExitMessage(); + } + + BSHP batchTracer; + iCreateDefaultSingleCell(batchTracer); + std::vector batchTraces; + std::vector batchTimes; + VecStr batchMessages; + batchTracer->TracePoints(seeds, seedTimes, batchTraces, batchTimes, batchMessages); + + TS_ASSERT_EQUALS(seeds.size(), batchTraces.size()); + TS_ASSERT_EQUALS(seeds.size(), batchTimes.size()); + TS_ASSERT_EQUALS(seeds.size(), batchMessages.size()); + for (size_t i = 0; i < seeds.size(); ++i) + { + TS_ASSERT_DELTA_VECPT3D(serialTraces[i], batchTraces[i], 1e-12); + TS_ASSERT_DELTA_VEC(serialTimes[i], batchTimes[i], 1e-12); + TS_ASSERT_EQUALS(serialMessages[i], batchMessages[i]); + // Positions and times are documented as parallel arrays, so a caller may zip them. + TS_ASSERT_EQUALS(batchTraces[i].size(), batchTimes[i].size()); + } + + // The seed outside the grid produces no polyline at all -- callers cannot assume one. + TS_ASSERT(batchTraces[2].empty()); + // ... while the two inside it do. + TS_ASSERT(batchTraces[0].size() >= 2); + TS_ASSERT(batchTraces[1].size() >= 2); + + // A caller that supplies the wrong number of start times has a bug; tracing the common + // prefix would hide it, so nothing is returned. + std::vector shortTraces; + std::vector shortTimes; + VecStr shortMessages; + batchTracer->TracePoints(seeds, {.5}, shortTraces, shortTimes, shortMessages); + TS_ASSERT(shortTraces.empty()); + TS_ASSERT(shortTimes.empty()); + TS_ASSERT(shortMessages.empty()); +} // XmGridTraceUnitTests::testTracePointsMatchesSerialTracePoint +//------------------------------------------------------------------------------ /// \brief Verifies the boundary-exit extractor is built once per tracer, not once per exit. /// /// The extractor's constructor triangulates the whole grid and its first SetPolyline indexes diff --git a/xmsgridtrace/gridtrace/XmGridTrace.h b/xmsgridtrace/gridtrace/XmGridTrace.h index 5d4c7f5..2c8137c 100644 --- a/xmsgridtrace/gridtrace/XmGridTrace.h +++ b/xmsgridtrace/gridtrace/XmGridTrace.h @@ -119,6 +119,33 @@ class XmGridTrace VecPt3d& a_outTrace, VecDbl& a_outTimes) = 0; + /// \brief Runs the Grid Trace for many points against the current two time steps. + /// + /// Equivalent to calling TracePoint once per point, but crossing a language or module + /// boundary once instead of once per point, and reporting why every trace ended rather + /// than only the last -- GetExitMessage describes a single operation, so it cannot + /// answer that for a batch. + /// + /// The time steps are not advanced: every trace runs against whichever pair + /// AddGridScalarsAtTime has most recently supplied. Callers wanting traces that span more + /// of a series feed the next time step and trace again. + /// + /// A_outTraces[i] can hold fewer than two points. A seed that leaves the grid on its very + /// first step yields only the seed itself, so callers must not assume one usable polyline + /// per point. + /// + /// \param[in] a_pts The starting point of each trace + /// \param[in] a_ptTimes The starting time of each trace; must be one per point + /// \param[out] a_outTraces The resultant positions at each step, one entry per point + /// \param[out] a_outTimes The resultant times, parallel to and the same length as + /// the matching entry of a_outTraces + /// \param[out] a_outExitMessages What ended each trace, one entry per point + virtual void TracePoints(const VecPt3d& a_pts, + const VecDbl& a_ptTimes, + std::vector& a_outTraces, + std::vector& a_outTimes, + VecStr& a_outExitMessages) = 0; + /// \brief returns a message describing what caused trace to exit /// \return the exit message of the last TracePoint operation virtual std::string GetExitMessage() = 0; diff --git a/xmsgridtrace/gridtrace/XmGridTrace.t.h b/xmsgridtrace/gridtrace/XmGridTrace.t.h index 97077cb..efe99b5 100644 --- a/xmsgridtrace/gridtrace/XmGridTrace.t.h +++ b/xmsgridtrace/gridtrace/XmGridTrace.t.h @@ -39,6 +39,7 @@ class XmGridTraceUnitTests : public CxxTest::TestSuite void testStartInactiveCell(); void testTutorial(); void testTimeVaryingFieldChangesPath(); + void testTracePointsMatchesSerialTracePoint(); void testBoundaryExtractorIsCached(); void testTraceBenchmark(); From b3ae94b0525937fdb31f46e5f44cc2bfb58db9df Mon Sep 17 00:00:00 2001 From: Bill Dolinar Date: Fri, 14 Aug 2026 09:02:33 -0600 Subject: [PATCH 05/10] Make traces resumable across time steps, and report exit reasons as an enum A trace can only run as far as the second of the two loaded time steps, because that is as far as the field is known. The batch added in the previous commit therefore traced every seed to the edge of one window and threw away everything it knew, which is not tracing a flow through time -- it is tracing it through one interval. Restarting from the last position would not fix that either, because a restart loses the trace's history. Traces now suspend and resume. StartTraces seeds a batch, ContinueTraces advances every unfinished trace and returns how many are waiting on a later time step, and the caller feeds the next one and calls again: tracer->StartTraces(seeds, seedTimes); while (tracer->ContinueTraces() > 0 && series.HasNext()) tracer->AddGridScalarsAtTime(series.Next(), ...); tracer->GetTraceResults(traces, times, reasons); Keeping the two-step window means memory stays bounded however long the series is, and the caller reads time steps only as the traces actually need them. Stopping early is legitimate: traces still waiting end where they got to, and say so. The substance is what survives a window change. Position and time are the obvious ones. The distance and elapsed-time budgets are whole-trace, not per-window, so they carry. So do the adaptive step size and the previous velocity, because the subdivision tests compare each step against the one before it -- restarting those at a boundary would kink the path exactly where the time step changes, which is the one place this has to be smooth. TracePoint's body is now StepTrace, which either starts a trace or resumes one; every exit from it routes through a single lambda that writes that state back, so there is no path that advances a trace without recording where it reached. Resumability cannot be read off the loop's final state, so it is tracked explicitly: a subdivision puts the trace back in motion *after* the time step clamp has already fired, and several conditions in one iteration overwrite each other. This also only works because of the interpolation fix. A trace resuming at the new first time step is advected by that step's field; under the inverted weights it would have used the following one, so every window boundary would have introduced an error. Since nothing outside this repository uses XmGridTrace, the surrounding API is cleaned up rather than extended around: - The exit reason is an enum, not a message. A caller has to tell "left the grid, draw it short" from "spent its distance budget, this is the normal ending" for tens of thousands of seeds, and string comparison cannot support that -- the old messages were composed by appending, so no fixed string identified a case. The strings remain, one per reason, for logs and tooltips. TracePoint gets GetExitReason so the single-point path can answer the same question the batch answers. - The one-shot TracePoints from the previous commit is gone; the resumable trio subsumes it, and one way to batch is better than two. - GetExitMessage returns const std::string& and is const. - AddGridScalarsAtTime takes activity by const reference. testTracesContinueAcrossTimeSteps is the guard, over a field rotating +x -> +y -> -x. Its strongest assertion is not that the continued trace is longer: it is that the trace never given the third time step is an exact prefix of the one that was. That is what shows resuming extends the path rather than recomputing it, and it is what fails if any carried state is dropped at the boundary. It also checks the path turns back on itself, which no single pair of those time steps can produce. Measured on the benchmark's realistic seed population, 43 of 250 seeds stop waiting for a later time step -- traces the previous design silently truncated. The Python bindings still compile against this by inspection: their three uses are lambdas or a member pointer that still resolves. Binding the new calls needs a python-enabled build, which the testing preset does not produce, and is not done here. --- xmsgridtrace/gridtrace/XmGridTrace.cpp | 496 ++++++++++++++++++------- xmsgridtrace/gridtrace/XmGridTrace.h | 98 +++-- xmsgridtrace/gridtrace/XmGridTrace.t.h | 3 +- 3 files changed, 435 insertions(+), 162 deletions(-) diff --git a/xmsgridtrace/gridtrace/XmGridTrace.cpp b/xmsgridtrace/gridtrace/XmGridTrace.cpp index 5e6035a..c977f12 100644 --- a/xmsgridtrace/gridtrace/XmGridTrace.cpp +++ b/xmsgridtrace/gridtrace/XmGridTrace.cpp @@ -70,6 +70,45 @@ size_t g_boundaryExtractorBuilds = 0; //----- Class / Function definitions ------------------------------------------- +//------------------------------------------------------------------------------ +/// \brief Whether a reason means the trace can never advance again. +/// \param[in] a_reason The exit reason +/// \return true if no amount of further time step data can move the trace +//------------------------------------------------------------------------------ +bool iIsTerminal(XmGridTraceExitEnum a_reason) +{ + return a_reason != GTEXIT_NOT_STARTED && a_reason != GTEXIT_WAITING_FOR_TIME_STEP; +} // iIsTerminal + +//////////////////////////////////////////////////////////////////////////////// +/// One trace in progress, and everything about it that has to survive a time step change. +/// +/// A trace stops when it reaches the second of the two loaded time steps and continues once +/// a later one is supplied. Position and time are the obvious carry-overs; the rest are the +/// ones whose absence would be a silent defect. The distance and elapsed-time budgets are +/// whole-trace, not per-window. The step size and previous velocity feed the subdivision +/// tests, which compare each step against the one before it -- restarting those at a window +/// boundary would kink the path exactly where the time step changes, which is the one place +/// this has to be smooth. +struct TraceState +{ + Pt3d m_pt; ///< current position + double m_ptTime = 0; ///< time the trace was released; never advanced + double m_elapsedTime = 0; ///< time advanced since release, against m_maxTracingTime + double m_distTraveled = 0; ///< distance covered, against m_maxTracingDistance + double m_deltaT = 1.0; ///< adaptive step size carried into the next step + double m_vx = 0; ///< velocity x at m_pt, for the subdivision tests + double m_vy = 0; ///< velocity y at m_pt, for the subdivision tests + double m_mag = 0; ///< speed at m_pt, for the change-in-velocity test + bool m_started = false; ///< the seed has been evaluated and recorded + /// Why it stopped, or that it is waiting. Doubles as the resume flag -- see iIsTerminal -- + /// so there is one source of truth rather than a reason and a separate finished bool that + /// could disagree. + XmGridTraceExitEnum m_exitReason = GTEXIT_NOT_STARTED; + VecPt3d m_trace; ///< positions so far + VecDbl m_times; ///< times so far, parallel to m_trace +}; + //////////////////////////////////////////////////////////////////////////////// /// Implementation for XmGridTrace class XmGridTraceImpl : public XmGridTrace @@ -101,7 +140,7 @@ class XmGridTraceImpl : public XmGridTrace void AddGridScalarsAtTime(const VecPt3d& a_scalars, DataLocationEnum a_scalarLoc, - xms::DynBitset& a_activity, + const xms::DynBitset& a_activity, DataLocationEnum a_activityLoc, double a_time) final; @@ -110,15 +149,17 @@ class XmGridTraceImpl : public XmGridTrace VecPt3d& a_outTrace, VecDbl& a_outTimes) final; - void TracePoints(const VecPt3d& a_pts, - const VecDbl& a_ptTimes, - std::vector& a_outTraces, - std::vector& a_outTimes, - VecStr& a_outExitMessages) final; + void StartTraces(const VecPt3d& a_pts, const VecDbl& a_ptTimes) final; + int ContinueTraces() final; + void GetTraceResults(std::vector& a_outTraces, + std::vector& a_outTimes, + std::vector& a_outExitReasons) const final; - std::string GetExitMessage() final; + const std::string& GetExitMessage() const final; private: + void StepTrace(TraceState& a_state); + bool GetVectorAtLocationAndTime(const xms::Pt3d& a_pt, double a_currentTime, xms::Pt3d& a_data) const; @@ -149,7 +190,10 @@ class XmGridTraceImpl : public XmGridTrace /// measured ~40 ms each. Null until a trace actually exits, so a tracer whose traces all /// stay inside the grid never pays the memory. BSHP m_boundaryExtractor; - double m_distTraveled=0; ///< distance traveled in the last TracePoint operation + /// Traces started by StartTracePoints and advanced by ContinueTracePoints. Empty unless + /// a batch is in flight; one batch per tracer, because the time step window it runs + /// against is itself instance state. + std::vector m_batch; std::string m_exitMessage; ///< exit message for the last TracePoint operation protected: @@ -292,7 +336,7 @@ void XmGridTraceImpl::SetMaxChangeDirectionInRadians(const double a_maxChangeDir //------------------------------------------------------------------------------ /// \brief returns a message describing what caused trace to exit //------------------------------------------------------------------------------ -std::string XmGridTraceImpl::GetExitMessage() +const std::string& XmGridTraceImpl::GetExitMessage() const { return m_exitMessage; } // XmGridTraceImpl::GetExitMessage @@ -308,7 +352,7 @@ std::string XmGridTraceImpl::GetExitMessage() //------------------------------------------------------------------------------ void XmGridTraceImpl::AddGridScalarsAtTime(const VecPt3d& a_scalars, DataLocationEnum a_scalarLoc, - xms::DynBitset& a_activity, + const xms::DynBitset& a_activity, DataLocationEnum a_activityLoc, double a_time) { @@ -341,50 +385,75 @@ void XmGridTraceImpl::AddGridScalarsAtTime(const VecPt3d& a_scalars, } //------------------------------------------------------------------------------ -/// \brief Runs the Grid Trace for a point -/// \param[in] a_pt The starting point of the trace -/// \param[in] a_ptTime The starting time of the trace -/// \param[out] a_outTrace the resultant positions at each step -/// \param[out] a_outTimes the resultant times at each step +/// \brief Advances one trace as far as the currently loaded pair of time steps allows. +/// +/// Starting a trace and resuming one differ only in the prologue: a fresh state has to +/// evaluate and record its seed, while a resumed one already carries a position, its +/// budgets, its step size and its previous velocity. +/// \param[in,out] a_state The trace to advance //------------------------------------------------------------------------------ -void XmGridTraceImpl::TracePoint(const Pt3d& a_pt, - const double& a_ptTime, - VecPt3d& a_outTrace, - VecDbl& a_outTimes) +void XmGridTraceImpl::StepTrace(TraceState& a_state) { - m_exitMessage.clear(); - double deltaT = 1.00; - double mag0 = 0, mag1 = 0; - Pt3d pt0 = a_pt, pt1; - double vx0 = 0, vx1 = 0, vy0 = 0, vy1 = 0, elapsedTime = 0; + if (iIsTerminal(a_state.m_exitReason)) + return; + + const double ptTime = a_state.m_ptTime; + Pt3d pt0 = a_state.m_pt, pt1; + double deltaT = a_state.m_deltaT; + double elapsedTime = a_state.m_elapsedTime; + double distTraveled = a_state.m_distTraveled; + double vx0 = a_state.m_vx, vy0 = a_state.m_vy, mag0 = a_state.m_mag; + double vx1 = 0, vy1 = 0, mag1 = 0; bool bContinue = true; - Pt3d vtkVec; // Rename this variable - Pt3d vtkPt; + Pt3d vtkVec; Pt3d vector; + VecPt3d& outTrace = a_state.m_trace; + VecDbl& outTimes = a_state.m_times; + + // Writes back everything the next call resumes from. Every exit from this function goes + // through it, so there is no path that advances the trace without recording where it got to. + auto stopWith = [&](XmGridTraceExitEnum a_reason) { + a_state.m_pt = pt0; + a_state.m_deltaT = deltaT; + a_state.m_elapsedTime = elapsedTime; + a_state.m_distTraveled = distTraveled; + a_state.m_vx = vx0; + a_state.m_vy = vy0; + a_state.m_mag = mag0; + a_state.m_exitReason = a_reason; + m_exitMessage = XmGridTraceExitReasonToString(a_reason); + }; - m_distTraveled = 0; - a_outTrace.clear(); - a_outTimes.clear(); - if (a_ptTime > m_time2 || // Test if the time specified is after the time range - !GetVectorAtLocationAndTime(a_pt, a_ptTime, vector)) // Ensure nothing fails during extraction - { - m_exitMessage = "Error occurred while extracting point0."; - return; - } - if (EQ_TOL(vector.x, XM_NODATA, 1) || EQ_TOL(vector.y, XM_NODATA, 1)) + if (!a_state.m_started) { - m_exitMessage = "Point does not start inside an active cell."; - return; - } + outTrace.clear(); + outTimes.clear(); + if (ptTime > m_time2 || // Test if the time specified is after the time range + !GetVectorAtLocationAndTime(pt0, ptTime, vector)) // Ensure extraction did not fail + { + stopWith(GTEXIT_EXTRACTION_FAILED); + return; + } + if (EQ_TOL(vector.x, XM_NODATA, 1) || EQ_TOL(vector.y, XM_NODATA, 1)) + { + stopWith(GTEXIT_SEED_NOT_TRACEABLE); + return; + } - a_outTrace.push_back(a_pt); - a_outTimes.push_back(a_ptTime); + outTrace.push_back(pt0); + outTimes.push_back(ptTime); - vx0 = vector.x * m_vectorMultiplier; - vy0 = vector.y * m_vectorMultiplier; - mag0 = sqrt(vector.x * vector.x + vector.y * vector.y); + vx0 = vector.x * m_vectorMultiplier; + vy0 = vector.y * m_vectorMultiplier; + mag0 = sqrt(vector.x * vector.x + vector.y * vector.y); + a_state.m_started = true; + } double maxAngleChange = cos(m_maxChangeDirectionInRadians); + // Which reason the loop will stop with. Tracked explicitly rather than inferred afterwards: + // several conditions in one iteration overwrite each other, and a later split can put the + // trace back into motion after the time step clamp has already fired. + XmGridTraceExitEnum stopReason = GTEXIT_WAITING_FOR_TIME_STEP; while (bContinue) { @@ -395,42 +464,38 @@ void XmGridTraceImpl::TracePoint(const Pt3d& a_pt, double denom = (vx0 * vx0) + (vy0 * vy0) + (m_maxChangeDistance * XM_ZERO_TOL); double dt = sqrt(d2 / denom); if (deltaT > dt) - { deltaT = dt; - m_exitMessage = "Change distance was greater than the max change distance."; - } } // If the change in DeltaT would push us beyond the time step, set it to hit the timestep - if (elapsedTime + deltaT + a_ptTime > m_time2) + if (elapsedTime + deltaT + ptTime > m_time2) { - deltaT = m_time2 - elapsedTime - a_ptTime; - bContinue = false; // This will be the last point traced - m_exitMessage = "The point has traveled beyond, or reached the second time step."; + deltaT = m_time2 - elapsedTime - ptTime; + bContinue = false; // This will be the last point traced in this window + stopReason = GTEXIT_WAITING_FOR_TIME_STEP; } - // If the change in delta time would push beyond the max tracing time, set it to hit max tracing - // time + // If the change in delta time would push beyond the max tracing time, set it to hit max + // tracing time if (m_maxTracingTime > 0 && (elapsedTime + deltaT) > m_maxTracingTime) { deltaT = m_maxTracingTime - elapsedTime; bContinue = false; // This will be the last point traced - m_exitMessage = "Exceeded or reached max tracing time."; + stopReason = GTEXIT_MAX_TRACING_TIME; } // compute candidate point pt1.x = pt0.x + deltaT * vx0; pt1.y = pt0.y + deltaT * vy0; - if (!GetVectorAtLocationAndTime(pt1, a_ptTime + elapsedTime + deltaT, vtkVec)) + if (!GetVectorAtLocationAndTime(pt1, ptTime + elapsedTime + deltaT, vtkVec)) { - a_outTrace.clear(); - a_outTimes.clear(); - m_exitMessage = "Error occurred while extracting point1"; + outTrace.clear(); + outTimes.clear(); + stopWith(GTEXIT_EXTRACTION_FAILED); return; } // if pt1 outside of domain, compute new deltaT to get to boundary if (EQ_TOL(vtkVec.x, XM_NODATA, 1) || EQ_TOL(vtkVec.y, XM_NODATA, 1)) { - m_exitMessage = "Point has traveled out of domain."; VecPt3d points = {pt0, pt1}; if (!m_boundaryExtractor) { @@ -446,6 +511,7 @@ void XmGridTraceImpl::TracePoint(const Pt3d& a_pt, if (points.size() < 3) { XM_LOG(xmlog::error, "Gridtracer failed to find an intersection when exiting grid."); + stopWith(GTEXIT_LEFT_GRID); return; } double segDist = Mdist(pt0.x, pt0.y, pt1.x, pt1.y); @@ -453,10 +519,11 @@ void XmGridTraceImpl::TracePoint(const Pt3d& a_pt, double newSegDist = Mdist(pt0.x, pt0.y, pt1.x, pt1.y); deltaT *= (newSegDist / segDist); bContinue = false; - if (!GetVectorAtLocationAndTime(pt1, a_ptTime + elapsedTime + deltaT, vtkVec) || + stopReason = GTEXIT_LEFT_GRID; + if (!GetVectorAtLocationAndTime(pt1, ptTime + elapsedTime + deltaT, vtkVec) || vtkVec.x == XM_NODATA || vtkVec.y == XM_NODATA) { - m_exitMessage = "Error occurred while extracting point1"; + stopWith(GTEXIT_EXTRACTION_FAILED); return; } } @@ -467,9 +534,11 @@ void XmGridTraceImpl::TracePoint(const Pt3d& a_pt, if (EQ_TOL(vx1, 0.0, .0001) && EQ_TOL(vy1, 0.0, .0001)) // No velocity { - a_outTrace.push_back(pt1); - a_outTimes.push_back(a_ptTime + elapsedTime + deltaT); - m_exitMessage = "Velocity has gone to zero."; + outTrace.push_back(pt1); + outTimes.push_back(ptTime + elapsedTime + deltaT); + pt0 = pt1; + elapsedTime += deltaT; + stopWith(GTEXIT_ZERO_VELOCITY); return; } bool bSplit = false; @@ -482,58 +551,58 @@ void XmGridTraceImpl::TracePoint(const Pt3d& a_pt, { double changeVel = fabs(mag1 - mag0); if (changeVel > m_maxChangeVelocity) - { bSplit = true; - m_exitMessage = "Point has exceeded max change velocity."; - } } if (!bSplit && m_maxChangeDirectionInRadians > 0) { double dir = iGetDirAsCosTheta(vx0, vy0, vx1, vy1); if (dir < maxAngleChange) - { bSplit = true; - m_exitMessage = "Point has exceeded max change direction."; - } } if (bSplit) { + // A split puts the trace back in motion, so any stop decided earlier in this iteration + // is void -- including the time step clamp, which is why resumability cannot be read + // off the loop's final state without this. bContinue = true; + stopReason = GTEXIT_WAITING_FOR_TIME_STEP; deltaT /= 2; if (m_minDeltaTime > 0 && deltaT < m_minDeltaTime) { // done, exit bContinue = false; - m_exitMessage += " Delta time was less than min delta time."; + stopReason = GTEXIT_MIN_DELTA_TIME; } } else { double segDist = Mdist(pt0.x, pt0.y, pt1.x, pt1.y); - m_distTraveled += segDist; - if (m_maxTracingDistance > 0 && m_distTraveled > m_maxTracingDistance) + distTraveled += segDist; + if (m_maxTracingDistance > 0 && distTraveled > m_maxTracingDistance) { // because our last point exceeded the exitDistance // find this point by linear calculations - double distancePast = m_distTraveled - m_maxTracingDistance; + double distancePast = distTraveled - m_maxTracingDistance; double perc = distancePast / segDist; Pt3d newPt; newPt.x = (pt0.x * perc) + (pt1.x * (1 - perc)); newPt.y = (pt0.y * perc) + (pt1.y * (1 - perc)); - m_distTraveled = m_maxTracingDistance; - a_outTrace.push_back(newPt); - a_outTimes.push_back(a_ptTime + elapsedTime + deltaT * perc); - m_exitMessage = "Point has reached or exceeded the max tracing distance."; + distTraveled = m_maxTracingDistance; + outTrace.push_back(newPt); + outTimes.push_back(ptTime + elapsedTime + deltaT * perc); + pt0 = newPt; + elapsedTime += deltaT * perc; + stopWith(GTEXIT_MAX_TRACING_DISTANCE); return; } // add new pt if not identical to last -- and push its time only when the point is // pushed. The time push used to be unconditional, so a step shorter than XM_ZERO_TOL - // left a_outTimes one longer than a_outTrace and silently misaligned every later - // pair, which a caller reading them as parallel arrays cannot detect. - const bool moved = a_outTrace.empty() || !EQ_TOL(pt1.x, a_outTrace.back().x, XM_ZERO_TOL) || - !EQ_TOL(pt1.y, a_outTrace.back().y, XM_ZERO_TOL); + // left the times array one longer and silently misaligned every later pair, which a + // caller reading them as parallel arrays cannot detect. + const bool moved = outTrace.empty() || !EQ_TOL(pt1.x, outTrace.back().x, XM_ZERO_TOL) || + !EQ_TOL(pt1.y, outTrace.back().y, XM_ZERO_TOL); pt0 = pt1; elapsedTime += deltaT; vx0 = vx1; @@ -542,46 +611,92 @@ void XmGridTraceImpl::TracePoint(const Pt3d& a_pt, mag0 = mag1; if (moved) { - a_outTrace.push_back(pt1); - a_outTimes.push_back(a_ptTime + elapsedTime); + outTrace.push_back(pt1); + outTimes.push_back(ptTime + elapsedTime); } } } // while () + stopWith(stopReason); +} // XmGridTraceImpl::StepTrace +//------------------------------------------------------------------------------ +/// \brief Runs the Grid Trace for a point against the currently loaded time steps +/// \param[in] a_pt The starting point of the trace +/// \param[in] a_ptTime The starting time of the trace +/// \param[out] a_outTrace the resultant positions at each step +/// \param[out] a_outTimes the resultant times at each step +//------------------------------------------------------------------------------ +void XmGridTraceImpl::TracePoint(const Pt3d& a_pt, + const double& a_ptTime, + VecPt3d& a_outTrace, + VecDbl& a_outTimes) +{ + TraceState state; + state.m_pt = a_pt; + state.m_ptTime = a_ptTime; + StepTrace(state); + a_outTrace.swap(state.m_trace); + a_outTimes.swap(state.m_times); } // XmGridTraceImpl::TracePoint //------------------------------------------------------------------------------ -/// \brief Runs the Grid Trace for many points against the current two time steps +/// \brief Begins tracing a batch of seeds against the currently loaded time steps /// \param[in] a_pts The starting point of each trace /// \param[in] a_ptTimes The starting time of each trace; must be one per point -/// \param[out] a_outTraces The resultant positions at each step, one entry per point -/// \param[out] a_outTimes The resultant times, one entry per point -/// \param[out] a_outExitMessages What ended each trace, one entry per point //------------------------------------------------------------------------------ -void XmGridTraceImpl::TracePoints(const VecPt3d& a_pts, - const VecDbl& a_ptTimes, - std::vector& a_outTraces, - std::vector& a_outTimes, - VecStr& a_outExitMessages) +void XmGridTraceImpl::StartTraces(const VecPt3d& a_pts, const VecDbl& a_ptTimes) { - a_outTraces.clear(); - a_outTimes.clear(); - a_outExitMessages.clear(); + m_batch.clear(); if (a_pts.size() != a_ptTimes.size()) { - // Returning empty rather than tracing the common prefix: a caller that mismatched these - // has a bug, and a partial result would let it go unnoticed. - XM_LOG(xmlog::error, "Gridtracer: TracePoints needs one start time per point."); + // Refusing the whole batch rather than seeding the common prefix: a caller that + // mismatched these has a bug, and a partial batch would let it go unnoticed. + XM_LOG(xmlog::error, "Gridtracer: StartTraces needs one start time per point."); return; } - - a_outTraces.resize(a_pts.size()); - a_outTimes.resize(a_pts.size()); - a_outExitMessages.resize(a_pts.size()); + m_batch.resize(a_pts.size()); for (size_t i = 0; i < a_pts.size(); ++i) { - TracePoint(a_pts[i], a_ptTimes[i], a_outTraces[i], a_outTimes[i]); - a_outExitMessages[i] = m_exitMessage; + m_batch[i].m_pt = a_pts[i]; + m_batch[i].m_ptTime = a_ptTimes[i]; } -} // XmGridTraceImpl::TracePoints +} // XmGridTraceImpl::StartTraces +//------------------------------------------------------------------------------ +/// \brief Advances every unfinished trace as far as the loaded time steps allow +/// \return How many traces are waiting on a later time step +//------------------------------------------------------------------------------ +int XmGridTraceImpl::ContinueTraces() +{ + int waiting = 0; + for (auto& state : m_batch) + { + StepTrace(state); // returns immediately for traces that are already finished + if (state.m_exitReason == GTEXIT_WAITING_FOR_TIME_STEP) + ++waiting; + } + return waiting; +} // XmGridTraceImpl::ContinueTraces +//------------------------------------------------------------------------------ +/// \brief Copies out the batch traced so far +/// \param[out] a_outTraces The positions of each trace, one entry per seed +/// \param[out] a_outTimes The times of each trace, parallel to a_outTraces +/// \param[out] a_outExitReasons Why each trace stopped, one entry per seed +//------------------------------------------------------------------------------ +void XmGridTraceImpl::GetTraceResults(std::vector& a_outTraces, + std::vector& a_outTimes, + std::vector& a_outExitReasons) const +{ + a_outTraces.clear(); + a_outTimes.clear(); + a_outExitReasons.clear(); + a_outTraces.reserve(m_batch.size()); + a_outTimes.reserve(m_batch.size()); + a_outExitReasons.reserve(m_batch.size()); + for (const auto& state : m_batch) + { + a_outTraces.push_back(state.m_trace); + a_outTimes.push_back(state.m_times); + a_outExitReasons.push_back(state.m_exitReason); + } +} // XmGridTraceImpl::GetTraceResults //------------------------------------------------------------------------------ /// \brief Returns the velocity scalar for a given point and time /// \param[in] a_pt The point @@ -676,6 +791,36 @@ BSHP XmGridTrace::New(std::shared_ptr a_ugrid) { return BSHP(new XmGridTraceImpl(a_ugrid)); } // XmGridTrace::New +//------------------------------------------------------------------------------ +/// \brief Returns a human-readable description of an exit reason. +/// \param[in] a_reason The exit reason +/// \return a description suitable for a log or a tooltip +//------------------------------------------------------------------------------ +const char* XmGridTraceExitReasonToString(XmGridTraceExitEnum a_reason) +{ + switch (a_reason) + { + case GTEXIT_NOT_STARTED: + return "Trace has not started."; + case GTEXIT_WAITING_FOR_TIME_STEP: + return "Trace reached the second time step and is waiting for a later one."; + case GTEXIT_MAX_TRACING_TIME: + return "Exceeded or reached max tracing time."; + case GTEXIT_MAX_TRACING_DISTANCE: + return "Point has reached or exceeded the max tracing distance."; + case GTEXIT_LEFT_GRID: + return "Point has traveled out of domain."; + case GTEXIT_ZERO_VELOCITY: + return "Velocity has gone to zero."; + case GTEXIT_MIN_DELTA_TIME: + return "Delta time was less than min delta time."; + case GTEXIT_SEED_NOT_TRACEABLE: + return "Point does not start inside an active cell."; + case GTEXIT_EXTRACTION_FAILED: + return "Error occurred while extracting a vector."; + } + return "Unknown exit reason."; +} // XmGridTraceExitReasonToString } // namespace xms #ifdef CXX_TEST @@ -685,6 +830,7 @@ BSHP XmGridTrace::New(std::shared_ptr a_ugrid) #include #include #include +#include #include #include @@ -1894,18 +2040,17 @@ void XmGridTraceUnitTests::testTimeVaryingFieldChangesPath() TS_ASSERT(rotating.back().x < frozen.back().x); } // XmGridTraceUnitTests::testTimeVaryingFieldChangesPath //------------------------------------------------------------------------------ -/// \brief The batch entry point returns exactly what serial TracePoint calls return. +/// \brief A single-window batch returns exactly what serial TracePoint calls return. /// -/// TracePoints exists to cross a language boundary once instead of once per seed, so its -/// value depends entirely on it being a faithful stand-in. Comparing against serial -/// TracePoint on an identical fixture is the strongest oracle available -- stronger than a -/// recorded baseline, which would drift with the tracer rather than pin the equivalence. +/// The batch exists to cross a language boundary once instead of once per seed, so its value +/// depends on being a faithful stand-in. Comparing against serial TracePoint on an identical +/// fixture is a stronger oracle than a recorded baseline, which would drift with the tracer +/// rather than pin the equivalence. /// -/// The seeds are chosen to cover the three shapes a caller has to handle: a trace that -/// leaves the grid, another that leaves it from elsewhere, and a seed outside the grid -/// entirely, which yields an empty trace rather than a polyline. +/// The seeds cover the shapes a caller has to handle: traces that leave the grid, and a seed +/// outside the grid entirely, which yields no polyline at all. //------------------------------------------------------------------------------ -void XmGridTraceUnitTests::testTracePointsMatchesSerialTracePoint() +void XmGridTraceUnitTests::testBatchMatchesSerialTracePoint() { const VecPt3d seeds = {{.5, .5, 0}, {.25, .75, 0}, {-.1, 0, 0}}; const VecDbl seedTimes = {.5, .5, .5}; @@ -1914,48 +2059,129 @@ void XmGridTraceUnitTests::testTracePointsMatchesSerialTracePoint() iCreateDefaultSingleCell(serialTracer); std::vector serialTraces(seeds.size()); std::vector serialTimes(seeds.size()); - VecStr serialMessages(seeds.size()); for (size_t i = 0; i < seeds.size(); ++i) - { serialTracer->TracePoint(seeds[i], seedTimes[i], serialTraces[i], serialTimes[i]); - serialMessages[i] = serialTracer->GetExitMessage(); - } BSHP batchTracer; iCreateDefaultSingleCell(batchTracer); + batchTracer->StartTraces(seeds, seedTimes); + batchTracer->ContinueTraces(); std::vector batchTraces; std::vector batchTimes; - VecStr batchMessages; - batchTracer->TracePoints(seeds, seedTimes, batchTraces, batchTimes, batchMessages); + std::vector reasons; + batchTracer->GetTraceResults(batchTraces, batchTimes, reasons); TS_ASSERT_EQUALS(seeds.size(), batchTraces.size()); TS_ASSERT_EQUALS(seeds.size(), batchTimes.size()); - TS_ASSERT_EQUALS(seeds.size(), batchMessages.size()); + TS_ASSERT_EQUALS(seeds.size(), reasons.size()); for (size_t i = 0; i < seeds.size(); ++i) { TS_ASSERT_DELTA_VECPT3D(serialTraces[i], batchTraces[i], 1e-12); TS_ASSERT_DELTA_VEC(serialTimes[i], batchTimes[i], 1e-12); - TS_ASSERT_EQUALS(serialMessages[i], batchMessages[i]); // Positions and times are documented as parallel arrays, so a caller may zip them. TS_ASSERT_EQUALS(batchTraces[i].size(), batchTimes[i].size()); } - // The seed outside the grid produces no polyline at all -- callers cannot assume one. + // The seed outside the grid produces no polyline -- callers cannot assume one per seed. TS_ASSERT(batchTraces[2].empty()); - // ... while the two inside it do. + TS_ASSERT_EQUALS((int)GTEXIT_SEED_NOT_TRACEABLE, (int)reasons[2]); TS_ASSERT(batchTraces[0].size() >= 2); TS_ASSERT(batchTraces[1].size() >= 2); - // A caller that supplies the wrong number of start times has a bug; tracing the common - // prefix would hide it, so nothing is returned. - std::vector shortTraces; - std::vector shortTimes; - VecStr shortMessages; - batchTracer->TracePoints(seeds, {.5}, shortTraces, shortTimes, shortMessages); - TS_ASSERT(shortTraces.empty()); - TS_ASSERT(shortTimes.empty()); - TS_ASSERT(shortMessages.empty()); -} // XmGridTraceUnitTests::testTracePointsMatchesSerialTracePoint + // A caller supplying the wrong number of start times has a bug; seeding the common prefix + // would hide it, so the whole batch is refused. + batchTracer->StartTraces(seeds, {.5}); + batchTracer->GetTraceResults(batchTraces, batchTimes, reasons); + TS_ASSERT(batchTraces.empty()); + TS_ASSERT(batchTimes.empty()); + TS_ASSERT(reasons.empty()); +} // XmGridTraceUnitTests::testBatchMatchesSerialTracePoint +//------------------------------------------------------------------------------ +/// \brief A trace continues past the second time step once a later one is supplied. +/// +/// This is the point of the whole batch design: the field is only known between the two +/// loaded time steps, so a trace that wants to run further has to stop, ask for more, and +/// resume where it was -- carrying its budgets, its adaptive step size and its previous +/// velocity with it. +/// +/// The strongest assertion here is not that the continued trace is longer. It is that the +/// trace which never received the third time step is a byte-for-byte *prefix* of the one +/// that did. That is what shows resuming extends the path rather than recomputing it, and it +/// is what would fail if any carried-over state were dropped at the window boundary. +//------------------------------------------------------------------------------ +void XmGridTraceUnitTests::testTracesContinueAcrossTimeSteps() +{ + // One cell spanning the domain, so the field is spatially uniform and every change in the + // path comes from time. It rotates +x -> +y -> -x across three time steps. + VecPt3d points = {{0, 0, 0}, {40, 0, 0}, {40, 40, 0}, {0, 40, 0}}; + VecInt cells = {XMU_QUAD, 4, 0, 1, 2, 3}; + std::shared_ptr ugrid = XmUGrid::New(points, cells); + + DynBitset activity; + activity.push_back(true); + const VecPt3d seeds = {{20, 10, 0}}; + const VecDbl seedTimes = {0}; + + auto buildTracer = [&]() { + BSHP tracer = XmGridTrace::New(ugrid); + tracer->SetVectorMultiplier(1); + tracer->SetMaxTracingTime(18); + tracer->SetMaxTracingDistance(1000); + tracer->SetMinDeltaTime(.01); + tracer->SetMaxChangeDistance(.5); + tracer->SetMaxChangeVelocity(-1); + tracer->SetMaxChangeDirectionInRadians(XM_PI); // never subdivide on direction + VecPt3d east = {{1, 0, 0}}, north = {{0, 1, 0}}; + tracer->AddGridScalarsAtTime(east, DataLocationEnum::LOC_CELLS, activity, + DataLocationEnum::LOC_CELLS, 0.0); + tracer->AddGridScalarsAtTime(north, DataLocationEnum::LOC_CELLS, activity, + DataLocationEnum::LOC_CELLS, 10.0); + return tracer; + }; + + // Never given the third time step: it must stop at the second and say so. + BSHP stopped = buildTracer(); + stopped->StartTraces(seeds, seedTimes); + TS_ASSERT_EQUALS(1, stopped->ContinueTraces()); + std::vector stoppedTraces; + std::vector stoppedTimes; + std::vector stoppedReasons; + stopped->GetTraceResults(stoppedTraces, stoppedTimes, stoppedReasons); + TS_ASSERT_EQUALS((int)GTEXIT_WAITING_FOR_TIME_STEP, (int)stoppedReasons[0]); + TS_ASSERT_DELTA(10.0, stoppedTimes[0].back(), 1e-9); + + // Given the third: it must resume and run out its tracing time instead. + BSHP continued = buildTracer(); + continued->StartTraces(seeds, seedTimes); + TS_ASSERT_EQUALS(1, continued->ContinueTraces()); + VecPt3d west = {{-1, 0, 0}}; + continued->AddGridScalarsAtTime(west, DataLocationEnum::LOC_CELLS, activity, + DataLocationEnum::LOC_CELLS, 20.0); + TS_ASSERT_EQUALS(0, continued->ContinueTraces()); + std::vector traces; + std::vector times; + std::vector reasons; + continued->GetTraceResults(traces, times, reasons); + TS_ASSERT_EQUALS((int)GTEXIT_MAX_TRACING_TIME, (int)reasons[0]); + TS_ASSERT_DELTA(18.0, times[0].back(), 1e-9); + + // Resuming extends; it does not restart. + TS_ASSERT(traces[0].size() > stoppedTraces[0].size()); + for (size_t i = 0; i < stoppedTraces[0].size(); ++i) + { + TS_ASSERT_DELTA(stoppedTraces[0][i].x, traces[0][i].x, 1e-12); + TS_ASSERT_DELTA(stoppedTraces[0][i].y, traces[0][i].y, 1e-12); + TS_ASSERT_DELTA(stoppedTimes[0][i], times[0][i], 1e-12); + } + + // The third time step reverses the eastward drift, so the path must turn back on itself -- + // something no single pair of these time steps can produce. + double maxX = traces[0][0].x; + for (const auto& pt : traces[0]) + maxX = std::max(maxX, pt.x); + TS_ASSERT(maxX > seeds[0].x); + TS_ASSERT(traces[0].back().x < maxX); +} // XmGridTraceUnitTests::testTracesContinueAcrossTimeSteps //------------------------------------------------------------------------------ /// \brief Verifies the boundary-exit extractor is built once per tracer, not once per exit. /// diff --git a/xmsgridtrace/gridtrace/XmGridTrace.h b/xmsgridtrace/gridtrace/XmGridTrace.h index 2c8137c..a9c6b24 100644 --- a/xmsgridtrace/gridtrace/XmGridTrace.h +++ b/xmsgridtrace/gridtrace/XmGridTrace.h @@ -33,6 +33,26 @@ class dyn_bitset; //----- Constants / Enumerations ----------------------------------------------- +/// \brief Why a trace stopped. +/// +/// Reported per trace instead of the message string it replaced. A batch traces every +/// visible glyph -- tens of thousands of them -- and a caller has to be able to tell "left +/// the grid, draw it short" from "spent its distance budget, this is the normal ending" +/// without comparing strings. The old messages could not support that anyway: they were +/// composed by appending, so no fixed string identified a case. +enum XmGridTraceExitEnum { + GTEXIT_NOT_STARTED, ///< no stepping has happened yet + GTEXIT_WAITING_FOR_TIME_STEP, ///< reached the 2nd loaded step; supply a later one to resume + GTEXIT_MAX_TRACING_TIME, ///< the trace spent its time budget + GTEXIT_MAX_TRACING_DISTANCE, ///< the trace spent its distance budget + GTEXIT_LEFT_GRID, ///< stepped out of the grid; the path stops at the boundary + GTEXIT_ZERO_VELOCITY, ///< the field went still under the particle + GTEXIT_MIN_DELTA_TIME, ///< subdividing reached the smallest allowed step + GTEXIT_SEED_NOT_TRACEABLE, ///< the seed was outside the grid or in an inactive cell + GTEXIT_EXTRACTION_FAILED ///< a field lookup failed; the trace is discarded +}; + + //----- Structs / Classes ------------------------------------------------------ //////////////////////////////////////////////////////////////////////////////// @@ -105,7 +125,7 @@ class XmGridTrace /// \param[in] a_time The time of the scalars virtual void AddGridScalarsAtTime(const VecPt3d& a_scalars, DataLocationEnum a_scalarLoc, - xms::DynBitset& a_activity, + const xms::DynBitset& a_activity, DataLocationEnum a_activityLoc, double a_time) = 0; @@ -119,36 +139,57 @@ class XmGridTrace VecPt3d& a_outTrace, VecDbl& a_outTimes) = 0; - /// \brief Runs the Grid Trace for many points against the current two time steps. + /// \brief Begins tracing a batch of seeds against the currently loaded time steps. /// - /// Equivalent to calling TracePoint once per point, but crossing a language or module - /// boundary once instead of once per point, and reporting why every trace ended rather - /// than only the last -- GetExitMessage describes a single operation, so it cannot - /// answer that for a batch. + /// A trace runs only as far as the second loaded time step, because that is as far as the + /// field is known. Supply the next time step with AddGridScalarsAtTime and call + /// ContinueTraces to carry every unfinished trace onward; the two-step window means memory + /// stays bounded however long the series is, and the caller reads time steps only as the + /// traces actually need them: /// - /// The time steps are not advanced: every trace runs against whichever pair - /// AddGridScalarsAtTime has most recently supplied. Callers wanting traces that span more - /// of a series feed the next time step and trace again. + /// \code + /// tracer->StartTraces(seeds, seedTimes); + /// while (tracer->ContinueTraces() > 0 && series.HasNext()) + /// tracer->AddGridScalarsAtTime(series.Next(), ...); + /// tracer->GetTraceResults(traces, times, reasons); + /// \endcode /// - /// A_outTraces[i] can hold fewer than two points. A seed that leaves the grid on its very - /// first step yields only the seed itself, so callers must not assume one usable polyline - /// per point. + /// Stopping early is legitimate: traces still waiting simply end where they got to, with + /// GTEXIT_WAITING_FOR_TIME_STEP. Calling ContinueTraces twice without supplying a time step + /// in between does no useful work. + /// + /// One batch is in flight per tracer, because the time step window it runs against is + /// itself state on the tracer. Starting a batch discards any previous one. /// /// \param[in] a_pts The starting point of each trace - /// \param[in] a_ptTimes The starting time of each trace; must be one per point - /// \param[out] a_outTraces The resultant positions at each step, one entry per point - /// \param[out] a_outTimes The resultant times, parallel to and the same length as - /// the matching entry of a_outTraces - /// \param[out] a_outExitMessages What ended each trace, one entry per point - virtual void TracePoints(const VecPt3d& a_pts, - const VecDbl& a_ptTimes, - std::vector& a_outTraces, - std::vector& a_outTimes, - VecStr& a_outExitMessages) = 0; - - /// \brief returns a message describing what caused trace to exit - /// \return the exit message of the last TracePoint operation - virtual std::string GetExitMessage() = 0; + /// \param[in] a_ptTimes The starting time of each trace; must be one per point, or the + /// batch is refused entirely + virtual void StartTraces(const VecPt3d& a_pts, const VecDbl& a_ptTimes) = 0; + + /// \brief Advances every unfinished trace as far as the loaded time steps allow. + /// \return How many traces are waiting on a later time step. Zero means every trace has + /// ended for a reason that more data cannot change. + virtual int ContinueTraces() = 0; + + /// \brief Copies out the batch traced so far. Valid at any point, complete once + /// ContinueTraces has returned zero. + /// + /// An entry can hold fewer than two points: a seed that leaves the grid on its very first + /// step yields only the seed itself, so callers must not assume one usable polyline per + /// seed. + /// + /// \param[out] a_outTraces The positions of each trace, one entry per seed + /// \param[out] a_outTimes The times of each trace, parallel to and the same length as the + /// matching entry of a_outTraces + /// \param[out] a_outExitReasons Why each trace stopped, one entry per seed + virtual void GetTraceResults(std::vector& a_outTraces, + std::vector& a_outTimes, + std::vector& a_outExitReasons) const = 0; + + /// \brief Returns a human-readable description of what ended the last trace operation. + /// Use GetTraceResults' exit reasons to make decisions; this is for display and logs. + /// \return the exit message of the last trace operation + virtual const std::string& GetExitMessage() const = 0; private: XM_DISALLOW_COPY_AND_ASSIGN(XmGridTrace) @@ -159,4 +200,9 @@ class XmGridTrace //----- Function prototypes ---------------------------------------------------- +/// \brief Returns a human-readable description of an exit reason. +/// \param[in] a_reason The exit reason +/// \return a description suitable for a log or a tooltip +const char* XmGridTraceExitReasonToString(XmGridTraceExitEnum a_reason); + } // namespace xms diff --git a/xmsgridtrace/gridtrace/XmGridTrace.t.h b/xmsgridtrace/gridtrace/XmGridTrace.t.h index efe99b5..2963876 100644 --- a/xmsgridtrace/gridtrace/XmGridTrace.t.h +++ b/xmsgridtrace/gridtrace/XmGridTrace.t.h @@ -39,7 +39,8 @@ class XmGridTraceUnitTests : public CxxTest::TestSuite void testStartInactiveCell(); void testTutorial(); void testTimeVaryingFieldChangesPath(); - void testTracePointsMatchesSerialTracePoint(); + void testBatchMatchesSerialTracePoint(); + void testTracesContinueAcrossTimeSteps(); void testBoundaryExtractorIsCached(); void testTraceBenchmark(); From 07233888921c804795d54d7f6643cac33a1c033d Mon Sep 17 00:00:00 2001 From: Bill Dolinar Date: Fri, 14 Aug 2026 09:35:21 -0600 Subject: [PATCH 06/10] Update the Python trace baselines and stop pinning exact float times Three Python tests mirror C++ cases whose expectations moved when the inverted time interpolation was fixed, so they carried the same wrong values: test_unique_time_steps, test_inactive_cell and test_tutorial. The new values are transcribed from the C++ source, where each first step was derived by hand rather than captured from the runner. test_max_tracing_distance failed for a different reason worth recording. Its positions matched to six decimals while one *time* differed by 4.4e-16 -- one ULP -- because the times were compared with assert_array_equal, exact float equality, while the positions in the same test were already compared approximately. Bisecting placed it on the interpolation fix, which was not the obvious answer: that test supplies identical scalars at both time steps, so the two weighted terms are the same pair of products and IEEE addition is commutative. The cause is FMA contraction. The compiler folds `d1 * w1 + d2 * w2` into a fused multiply-add, computing one product exactly inside the FMA and rounding the other; swapping which weight multiplies which time step therefore moves the last bit even though the mathematics is unchanged. So the assertion was pinning the compiler's contraction decision rather than the tracer's behaviour. All fifteen of these time comparisons had the same latent fragility and only one happened to trip; they now compare approximately, matching what the same tests already do for positions and what the C++ mirrors do. The reason is recorded in the class docstring so a future reader does not tighten them again. The bisect also confirmed something worth having checked: caching the boundary-exit polyline extractor changes no numeric result. The commit that introduced it passes all sixteen Python tests, which the reused GmMultiPolyIntersector could plausibly not have done. flake8 is not installed in this environment, so the Python edit is unlinted. --- _package/tests/XmGridTrace_pyt.py | 205 ++++++++++++++++-------------- 1 file changed, 110 insertions(+), 95 deletions(-) diff --git a/_package/tests/XmGridTrace_pyt.py b/_package/tests/XmGridTrace_pyt.py index 9a17d48..e7b095e 100644 --- a/_package/tests/XmGridTrace_pyt.py +++ b/_package/tests/XmGridTrace_pyt.py @@ -9,7 +9,16 @@ class TestGridTrace(unittest.TestCase): - """GridTrace tests.""" + """GridTrace tests. + + Traced times are compared approximately, not exactly. They are doubles derived from + float32 grid scalars, and the compiler may contract ``a * b + c * d`` into an FMA -- which + of the two products lands inside the FMA is computed exactly while the other is rounded, + so the last bit depends on the order the terms are written in. These assertions used + ``assert_array_equal``, which pinned that decision rather than the tracer's behaviour, and + it broke on a correct interpolation fix that only swapped which weight multiplies which + time step. Positions were already compared approximately; times now match. + """ def create_default_single_cell(self): """Create a default single cell. @@ -70,7 +79,7 @@ def test_basic_trace_point(self): expected_out_trace = [(.5, .5, 0), (1, 1, 0)] expected_out_times = [.5, 1] np.testing.assert_array_almost_equal(expected_out_trace, result_tuple[0]) - np.testing.assert_array_equal(expected_out_times, result_tuple[1]) + np.testing.assert_array_almost_equal(expected_out_times, result_tuple[1]) def test_max_change_distance(self): """Test max change distance functionality.""" @@ -86,7 +95,7 @@ def test_max_change_distance(self): (1, 1, 0)] expected_out_times = [.5, 0.67677668424809445, 0.85355336849618890, 1] np.testing.assert_array_almost_equal(expected_out_trace, result_tuple[0]) - np.testing.assert_array_equal(expected_out_times, result_tuple[1]) + np.testing.assert_array_almost_equal(expected_out_times, result_tuple[1]) def test_small_scalars_trace_point(self): """Test functionality with small scalars.""" @@ -190,7 +199,7 @@ def test_strong_direction_change(self): 9.7883171816902319, 10.000000000000000] np.testing.assert_array_almost_equal(expected_out_trace, result_tuple[0]) - np.testing.assert_array_equal(expected_out_times, result_tuple[1]) + np.testing.assert_array_almost_equal(expected_out_times, result_tuple[1]) def test_max_tracing_time(self): """Test functionality of max tracing time.""" @@ -244,7 +253,7 @@ def test_max_tracing_time(self): 5.2587764123320317, 5.5] np.testing.assert_array_almost_equal(expected_out_trace, result_tuple[0]) - np.testing.assert_array_equal(expected_out_times, result_tuple[1]) + np.testing.assert_array_almost_equal(expected_out_times, result_tuple[1]) def test_max_tracing_distance(self): """Test functionality of max tracing distance.""" @@ -278,7 +287,7 @@ def test_max_tracing_distance(self): 2.1962400000000004, 2.4774609356360582] np.testing.assert_array_almost_equal(expected_out_trace, result_tuple[0], 6) - np.testing.assert_array_equal(expected_out_times, result_tuple[1]) + np.testing.assert_array_almost_equal(expected_out_times, result_tuple[1]) def test_start_out_of_cell(self): """Test functionality of starting outside of cell.""" @@ -289,7 +298,7 @@ def test_start_out_of_cell(self): expected_out_times = [] np.testing.assert_equal(0, len(result_tuple[0])) - np.testing.assert_array_equal(expected_out_times, result_tuple[1]) + np.testing.assert_array_almost_equal(expected_out_times, result_tuple[1]) def test_beyond_timestep(self): """Test functionality of starting beyond the time step.""" @@ -300,7 +309,7 @@ def test_beyond_timestep(self): expected_out_times = [] np.testing.assert_equal(0, len(result_tuple[0])) - np.testing.assert_array_equal(expected_out_times, result_tuple[1]) + np.testing.assert_array_almost_equal(expected_out_times, result_tuple[1]) def test_before_timestep(self): """Test functionality of starting before the time step.""" @@ -312,7 +321,7 @@ def test_before_timestep(self): expected_out_trace = [(.5, .5, 0), (1, 1, 0)] expected_out_times = [-.1, .4] np.testing.assert_array_almost_equal(expected_out_trace, result_tuple[0]) - np.testing.assert_array_equal(expected_out_times, result_tuple[1]) + np.testing.assert_array_almost_equal(expected_out_times, result_tuple[1]) def test_vector_multiplier(self): """Test functionality of vector multiplier.""" @@ -364,7 +373,7 @@ def test_vector_multiplier(self): 9.5360834004582404, 10.000000000000000] np.testing.assert_array_almost_equal(expected_out_trace, result_tuple[0]) - np.testing.assert_array_equal(expected_out_times, result_tuple[1]) + np.testing.assert_array_almost_equal(expected_out_times, result_tuple[1]) def test_multi_cell(self): """Test default functionality of multiple cells.""" @@ -390,7 +399,7 @@ def test_multi_cell(self): 9.9299199999999992, 9.9683860530914945] np.testing.assert_array_almost_equal(expected_out_trace, result_tuple[0]) - np.testing.assert_array_equal(expected_out_times, result_tuple[1]) + np.testing.assert_array_almost_equal(expected_out_times, result_tuple[1]) def test_max_change_velocity(self): """Test functionality of max change in velocity.""" @@ -442,7 +451,7 @@ def test_max_change_velocity(self): 9.1917078801783187, 9.6267364611093829] np.testing.assert_array_almost_equal(expected_out_trace, result_tuple[0]) - np.testing.assert_array_equal(expected_out_times, result_tuple[1]) + np.testing.assert_array_almost_equal(expected_out_times, result_tuple[1]) def test_unique_time_steps(self): """Test functionality of unique time steps.""" @@ -455,20 +464,22 @@ def test_unique_time_steps(self): result_tuple = tracer.trace_point((.5, .5, 0), start_time) - expected_out_trace = [(.5, .5, 0), - (0.70000000298023224, 0.50000000000000000, 0.00000000000000000), - (0.95200000226497650, 0.50000000000000000, 0.00000000000000000), - (1.2734079944372176, 0.50000000000000000, 0.00000000000000000), - (1.6897536998434066, 0.50000000000000000, 0.00000000000000000), - (2, .5, 0)] + expected_out_trace = [(0.5, 0.5, 0), + (0.60000000149011612, 0.5, 0), + (0.74400000184774395, 0.5, 0), + (0.95481600679159162, 0.5, 0), + (1.2691074101881981, 0.5, 0), + (1.747260385068264, 0.5, 0), + (2, 0.5, 0)] expected_out_times = [10, - 11.000000000000000, + 11, 12.199999999999999, 13.640000000000001, - 15.368000000000000, - 16.627525378316030] + 15.368, + 17.441600000000001, + 18.362609001148471] np.testing.assert_array_almost_equal(expected_out_trace, result_tuple[0]) - np.testing.assert_array_equal(expected_out_times, result_tuple[1]) + np.testing.assert_array_almost_equal(expected_out_times, result_tuple[1]) def test_inactive_cell(self): """Test functionality of inactive cells.""" @@ -482,16 +493,18 @@ def test_inactive_cell(self): result_tuple = tracer.trace_point((.5, .5, 0), start_time) - expected_out_trace = [(.5, .5, 0), - (0.70000000298023224, 0.50000000000000000, 0.00000000000000000), - (0.93040000677108770, 0.50000000000000000, 0.00000000000000000), - (0.99788877571821222, 0.50000000000000000, 0.00000000000000000)] + expected_out_trace = [(0.5, 0.5, 0), + (0.60000000149011612, 0.5, 0), + (0.74280000120401379, 0.5, 0), + (0.94575130454301826, 0.5, 0), + (1, 0.5, 0)] expected_out_times = [10, - 11.000000000000000, + 11, 12.199999999999999, - 12.560000000000000] + 13.640000000000001, + 13.969279307058475] np.testing.assert_array_almost_equal(expected_out_trace, result_tuple[0]) - np.testing.assert_array_equal(expected_out_times, result_tuple[1]) + np.testing.assert_array_almost_equal(expected_out_times, result_tuple[1]) def test_start_inactive_cell(self): """Test functionality of starting in an inactive cell.""" @@ -507,7 +520,7 @@ def test_start_inactive_cell(self): expected_out_times = [] np.testing.assert_equal(0, len(result_tuple[0])) - np.testing.assert_array_equal(expected_out_times, result_tuple[1]) + np.testing.assert_array_almost_equal(expected_out_times, result_tuple[1]) def test_tutorial(self): """A test to serve as a tutorial.""" @@ -566,69 +579,71 @@ def test_tutorial(self): print(tracer.get_exit_message()) # Expected values for this simulation - expected_out_trace = [(0.50000000000000000, 0.50000000000000000, 0.00000000000000000), - (0.50000000000000000, 1.2500000000000000, 0.00000000000000000), - (0.54457812566426578, 1.3391562513285316, 0.00000000000000000), - (0.61632493250262921, 1.4354984729093498, 0.00000000000000000), - (0.72535406450374607, 1.5315533661126233, 0.00000000000000000), - (0.88236797164001590, 1.6126801842666139, 0.00000000000000000), - (0.98873181403598276, 1.6331015959080102, 0.00000000000000000), - (1.0538503898747653, 1.6342606013582104, 0.00000000000000000), - (1.1249433009705341, 1.5683006835455087, 0.00000000000000000), - (1.1895097427498795, 1.3863448896225066, 0.00000000000000000), - (1.2235242118635632, 1.0588590059131318, 0.00000000000000000), - (1.2235242118635632, 0.90477286425654002, 0.00000000000000000), - (1.2005336220528682, 0.85080764250970042, 0.00000000000000000), - (1.1581790674742278, 0.79387770198395835, 0.00000000000000000), - (1.0896874578697060, 0.74131697161132859, 0.00000000000000000), - (0.98966250551038770, 0.70663752692174131, 0.00000000000000000), - (0.95806149614159530, 0.71817980325332686, 0.00000000000000000), - (0.92629620502521459, 0.77371504022050730, 0.00000000000000000), - (0.90239412753251202, 0.88917318465162865, 0.00000000000000000), - (0.89995172701803572, 1.0694875660697027, 0.00000000000000000), - (0.91503139037776327, 1.0911992829869794, 0.00000000000000000), - (0.93816744602651825, 1.1127546977629765, 0.00000000000000000), - (0.97140028507849163, 1.1309789606067331, 0.00000000000000000), - (0.99364912627842006, 1.1358370729524059, 0.00000000000000000), - (1.0071524474802995, 1.1364684019706512, 0.00000000000000000), - (1.0223447138862345, 1.1280655805979485, 0.00000000000000000), - (1.0369737821057583, 1.0971462034407997, 0.00000000000000000), - (1.0467397711865176, 1.0371377237101163, 0.00000000000000000), - (1.0467397711865176, 0.96499504248441559, 0.00000000000000000), - (1.0390576209755447, 0.95473758230148376, 0.00000000000000000), - (1.0276444556154691, 0.94488898976070590, 0.00000000000000000), - (1.0208791233912420, 0.94149540451099356, 0.00000000000000000)] - expected_out_times = [0.00000000000000000, - 0.37500000000000000, - 0.82499999999999996, - 1.3649999999999998, - 2.0129999999999999, - 2.7905999999999995, - 3.2571599999999994, - 3.5370959999999991, - 3.8730191999999990, - 4.2761270399999987, - 4.7598564479999981, - 5.3403317375999979, - 6.0369020851199977, - 6.8727865021439971, - 7.8758478025727969, - 9.0795213630873555, - 9.4406234312417237, - 9.8739459130269651, - 10.393932891169255, - 11.017917264940003, - 11.766698513464901, - 12.665236011694777, - 13.743481009570628, - 14.390428008296139, - 14.778596207531445, - 15.244398046613812, - 15.803360253512654, - 16.474114901791264, - 17.279020479725595, - 18.244907173246794, - 19.403971205472232, - 20.000000000000000] + expected_out_trace = [(0.5, 0.5, 0), + (0.5, 1.5, 0), + (0.62600000187754634, 1.6260000018775462, 0), + (0.82611968728899965, 1.7455603212296962, 0), + (0.97840008102011689, 1.7810753047635555, 0), + (1.0280095840364933, 1.7824472100312621, 0), + (1.0861189816907613, 1.7608732599310344, 0), + (1.1492686295114336, 1.6802752810470523, 0), + (1.2097920698566107, 1.5101408581884392, 0), + (1.2515951471975522, 1.2181485463468757, 0), + (1.2515951471975522, 0.84053651390559747, 0), + (1.2181758214493843, 0.78780883088769804, 0), + (1.1632869448015855, 0.73137186792498654, 0), + (1.0771209832183524, 0.67899546053648097, 0), + (1.0129487663521615, 0.66357815692798783, 0), + (0.97169356095126669, 0.66199025753694563, 0), + (0.92552080990281416, 0.70419149113367874, 0), + (0.88530832700558759, 0.83950990950827409, 0), + (0.87513974259796246, 1.0941588844381676, 0), + (0.90077009637050098, 1.128146252166127, 0), + (0.943692705404238, 1.1613833261644337, 0), + (0.97709108330292604, 1.1730361561747586, 0), + (0.99894959169213471, 1.1759300874982919, 0), + (1.0124203987349505, 1.1760105163064269, 0), + (1.0275428271398932, 1.1645289800266216, 0), + (1.042848666622334, 1.1337546211004945, 0), + (1.055142468614698, 1.0758075939238765, 0), + (1.0585305184379035, 0.98540145004498747, 0), + (1.0556233679912082, 0.97374570199926891, 0), + (1.0492587242876892, 0.9602613226646981, 0), + (1.0375007181419984, 0.94568649411103145, 0), + (1.017827020259642, 0.93210280494582176, 0), + (1.0175992759724071, 0.93204300863222744, 0)] + expected_out_times = [0, + 1, + 2.2000000000000002, + 3.6400000000000001, + 4.5040000000000004, + 4.7632000000000003, + 5.0742400000000005, + 5.4474880000000008, + 5.8953856000000009, + 6.432862720000001, + 7.0778352640000008, + 7.8518023168000006, + 8.7805627801600004, + 9.8950753361920007, + 10.563782869811201, + 10.96500738998272, + 11.446476814188543, + 12.024240123235531, + 12.717556094091917, + 13.54953525911958, + 14.547910257152775, + 15.146935255972693, + 15.506350255264643, + 15.721999254839814, + 15.980778054330019, + 16.291312613718265, + 16.663954084984159, + 17.111123850503233, + 17.647727569126122, + 18.291652031473589, + 19.064361386290546, + 19.991612612070895, + 20] np.testing.assert_array_almost_equal(expected_out_trace, result_tuple[0]) - np.testing.assert_array_equal(expected_out_times, result_tuple[1]) + np.testing.assert_array_almost_equal(expected_out_times, result_tuple[1]) From a61c129d6c9a23b57591368372fff8289faae972 Mon Sep 17 00:00:00 2001 From: Bill Dolinar Date: Fri, 14 Aug 2026 09:58:53 -0600 Subject: [PATCH 07/10] Bind the resumable trace API to Python, and add the GetExitReason that was missing Adds start_traces, continue_traces, get_trace_results and get_exit_reason to the pybind11 module, exports XmGridTraceExitEnum as exit_reason_enum, and forwards all four through the hand-written wrapper so FlowPathService can drive the resume loop from Python: tracer.start_traces(seeds, seed_times) while tracer.continue_traces() > 0: step = series.next() if step is None: break tracer.add_grid_scalars_at_time(*step) traces, times, reasons = tracer.get_trace_results() continue_traces releases the GIL, and only it. Tracing tens of thousands of seeds takes long enough that holding the GIL would stall the interpreter for a caller on a worker thread, which is exactly how this is meant to run. start_traces keeps the GIL because it converts Python iterables inside its lambda. start_traces raises ValueError when the start times do not match the points. The C++ side refuses the batch and returns empty, which from Python would look like a tracer that silently did nothing. GetExitReason is added here rather than earlier because it was never actually added. The commit that claimed it used a scripted string replacement with no assertion that the anchor matched, so the edit silently did nothing; the build then passed because nothing had changed, and the claim went unverified into that commit message. It is now on the interface, the impl, and covered by a test that asserts the single-point path and the batch report the same reason for the same seed -- if those can disagree, a caller cannot use them interchangeably. The lesson is in the tooling, not the code: scripted edits need an assertion that the anchor was found, and a green build is not evidence that an edit landed. Three Python tests cover the new surface: a trace that continues across three time steps and whose stopped-early result is a prefix of the continued one, the ValueError, and the batch matching serial trace_point calls. Two clang-format suggestions are deliberately not applied. It wants GetExitReason collapsed to one line, along with all nine sibling accessors that are not written that way; and it wants 425 of the 437 lines of XmGridTrace_py.cpp reformatted, a file never kept under clang-format. Both would bury the change in unrelated churn. Also corrects the previous commit's claim that the Python test edit was unlinted -- flake8 is installed now and it is clean. The only findings in _package are eight pre-existing AQU104 import-header comments, four of them in a file this branch never touched. --- _package/tests/XmGridTrace_pyt.py | 85 ++++++++++++- _package/xms/gridtrace/__init__.py | 1 + _package/xms/gridtrace/grid_trace.py | 65 ++++++++++ xmsgridtrace/gridtrace/XmGridTrace.cpp | 23 +++- xmsgridtrace/gridtrace/XmGridTrace.h | 11 +- .../python/gridtrace/XmGridTrace_py.cpp | 113 ++++++++++++++++++ 6 files changed, 295 insertions(+), 3 deletions(-) diff --git a/_package/tests/XmGridTrace_pyt.py b/_package/tests/XmGridTrace_pyt.py index e7b095e..2272bd4 100644 --- a/_package/tests/XmGridTrace_pyt.py +++ b/_package/tests/XmGridTrace_pyt.py @@ -5,7 +5,7 @@ from xms.grid.ugrid import UGrid -from xms.gridtrace import GridTrace +from xms.gridtrace import exit_reason_enum, GridTrace class TestGridTrace(unittest.TestCase): @@ -255,6 +255,89 @@ def test_max_tracing_time(self): np.testing.assert_array_almost_equal(expected_out_trace, result_tuple[0]) np.testing.assert_array_almost_equal(expected_out_times, result_tuple[1]) + def create_rotating_field_tracer(self): + """Create a tracer over one cell spanning the domain, with the field rotating +x -> +y. + + One cell means the field is spatially uniform, so any change in a path comes from time. + + Returns: + GridTrace: A tracer with two time steps loaded + """ + points = [(0, 0, 0), (40, 0, 0), (40, 40, 0), (0, 40, 0)] + cells = [UGrid.cell_type_enum.QUAD, 4, 0, 1, 2, 3] + tracer = GridTrace(UGrid(points, cells)) + tracer.vector_multiplier = 1 + tracer.max_tracing_time = 18 + tracer.max_tracing_distance = 1000 + tracer.min_delta_time = .01 + tracer.max_change_distance = .5 + tracer.max_change_velocity = -1 + tracer.max_change_direction_in_radians = np.pi # never subdivide on direction + tracer.add_grid_scalars_at_time([(1, 0, 0)], 'cells', [True], 'cells', 0) + tracer.add_grid_scalars_at_time([(0, 1, 0)], 'cells', [True], 'cells', 10) + return tracer + + def test_traces_continue_across_time_steps(self): + """A trace continues past the second time step once a later one is supplied.""" + seeds = [(20, 10, 0)] + seed_times = [0] + + # Never given the third time step: it must stop at the second and say so. + stopped = self.create_rotating_field_tracer() + stopped.start_traces(seeds, seed_times) + self.assertEqual(1, stopped.continue_traces()) + stopped_traces, stopped_times, stopped_reasons = stopped.get_trace_results() + self.assertEqual(exit_reason_enum.WAITING_FOR_TIME_STEP, stopped_reasons[0]) + self.assertAlmostEqual(10.0, stopped_times[0][-1]) + + # Given the third: it must resume and run out its tracing time instead. + tracer = self.create_rotating_field_tracer() + tracer.start_traces(seeds, seed_times) + self.assertEqual(1, tracer.continue_traces()) + tracer.add_grid_scalars_at_time([(-1, 0, 0)], 'cells', [True], 'cells', 20) + self.assertEqual(0, tracer.continue_traces()) + traces, times, reasons = tracer.get_trace_results() + self.assertEqual(exit_reason_enum.MAX_TRACING_TIME, reasons[0]) + self.assertAlmostEqual(18.0, times[0][-1]) + + # Resuming extends the path; it does not restart it. + self.assertGreater(len(traces[0]), len(stopped_traces[0])) + np.testing.assert_array_almost_equal(stopped_traces[0], traces[0][:len(stopped_traces[0])]) + np.testing.assert_array_almost_equal(stopped_times[0], times[0][:len(stopped_times[0])]) + + # The third time step reverses the eastward drift, so the path turns back on itself -- + # something no single pair of these time steps can produce. + max_x = max(pt[0] for pt in traces[0]) + self.assertGreater(max_x, seeds[0][0]) + self.assertLess(traces[0][-1][0], max_x) + + def test_start_traces_rejects_mismatched_times(self): + """A caller supplying the wrong number of start times gets an error, not a silent no-op.""" + tracer = self.create_rotating_field_tracer() + with self.assertRaises(ValueError): + tracer.start_traces([(20, 10, 0), (21, 10, 0)], [0]) + + def test_batch_matches_trace_point(self): + """The batch returns what serial trace_point calls return.""" + seeds = [(.5, .5, 0), (.25, .75, 0), (-.1, 0, 0)] + seed_times = [.5, .5, .5] + + serial = self.create_default_single_cell() + expected = [serial.trace_point(pt, t) for pt, t in zip(seeds, seed_times)] + + batch = self.create_default_single_cell() + batch.start_traces(seeds, seed_times) + batch.continue_traces() + traces, times, reasons = batch.get_trace_results() + + self.assertEqual(len(seeds), len(traces)) + for i, (expected_trace, expected_times) in enumerate(expected): + np.testing.assert_array_almost_equal(expected_trace, traces[i]) + np.testing.assert_array_almost_equal(expected_times, times[i]) + # The seed outside the grid yields no polyline -- callers cannot assume one per seed. + self.assertEqual(0, len(traces[2])) + self.assertEqual(exit_reason_enum.SEED_NOT_TRACEABLE, reasons[2]) + def test_max_tracing_distance(self): """Test functionality of max tracing distance.""" tracer = self.create_default_single_cell() diff --git a/_package/xms/gridtrace/__init__.py b/_package/xms/gridtrace/__init__.py index 40ace8f..22c712f 100644 --- a/_package/xms/gridtrace/__init__.py +++ b/_package/xms/gridtrace/__init__.py @@ -1,3 +1,4 @@ """Initialize the module.""" from ._xmsgridtrace import __version__ # NOQA: F401 +from ._xmsgridtrace.gridtrace import exit_reason_enum # NOQA: F401 from .grid_trace import GridTrace # NOQA: F401 diff --git a/_package/xms/gridtrace/grid_trace.py b/_package/xms/gridtrace/grid_trace.py index 11d2f95..3ab03cf 100644 --- a/_package/xms/gridtrace/grid_trace.py +++ b/_package/xms/gridtrace/grid_trace.py @@ -150,3 +150,68 @@ def get_exit_message(self): str: The exit message of the last trace_point operation """ return self._instance.get_exit_message() + + def get_exit_reason(self): + """Returns why the last trace operation ended. + + Prefer this over get_exit_message when deciding what to do with a trace; the message is for + display. WAITING_FOR_TIME_STEP means the path stops early because the field is not known past + the second loaded time step, not that the particle came to rest. + + Returns: + exit_reason_enum: The exit reason of the last trace operation + """ + return self._instance.get_exit_reason() + + def start_traces(self, pts, pt_times): + """Begin tracing a batch of seeds against the currently loaded time steps. + + A trace runs only as far as the second loaded time step, because that is as far as the field is + known. Supply the next time step with add_grid_scalars_at_time and call continue_traces to carry + every unfinished trace onward:: + + tracer.start_traces(seeds, seed_times) + while tracer.continue_traces() > 0: + step = series.next() + if step is None: + break + tracer.add_grid_scalars_at_time(*step) + traces, times, reasons = tracer.get_trace_results() + + Stopping early is fine: traces still waiting end where they got to. One batch is in flight per + tracer; starting a batch discards any previous one. + + Args: + pts (iterable): The starting point of each trace + pt_times (iterable): The starting time of each trace, one per point + + Raises: + ValueError: If pt_times does not have one entry per point + """ + self._instance.start_traces(pts, pt_times) + + def continue_traces(self): + """Advance every unfinished trace as far as the loaded time steps allow. + + Releases the GIL while tracing, so calling this from a worker thread does not stall the + interpreter. + + Returns: + int: How many traces are waiting on a later time step. Zero means every trace has ended for + a reason more data cannot change + """ + return self._instance.continue_traces() + + def get_trace_results(self): + """Return the batch traced so far. + + Valid at any point, complete once continue_traces has returned zero. An entry can hold fewer + than two points: a seed that leaves the grid on its first step yields only the seed itself, so + callers must not assume one usable polyline per seed. + + Returns: + tuple: The positions of each trace, the times of each trace, and why each trace stopped as + an exit_reason_enum. All three are parallel to the seeds passed to start_traces, and each + entry's times are parallel to its positions + """ + return self._instance.get_trace_results() diff --git a/xmsgridtrace/gridtrace/XmGridTrace.cpp b/xmsgridtrace/gridtrace/XmGridTrace.cpp index c977f12..78dfce5 100644 --- a/xmsgridtrace/gridtrace/XmGridTrace.cpp +++ b/xmsgridtrace/gridtrace/XmGridTrace.cpp @@ -155,6 +155,7 @@ class XmGridTraceImpl : public XmGridTrace std::vector& a_outTimes, std::vector& a_outExitReasons) const final; + XmGridTraceExitEnum GetExitReason() const final; const std::string& GetExitMessage() const final; private: @@ -195,7 +196,10 @@ class XmGridTraceImpl : public XmGridTrace /// against is itself instance state. std::vector m_batch; - std::string m_exitMessage; ///< exit message for the last TracePoint operation + /// Why the last trace operation ended. Kept beside the message so the single-point + /// TracePoint can answer the same question GetTraceResults answers per seed. + XmGridTraceExitEnum m_exitReason = GTEXIT_NOT_STARTED; + std::string m_exitMessage; ///< exit message for the last trace operation protected: }; double iGetDirAsCosTheta(double a_vx0, double a_vy0, double a_vx1, double a_vy1) @@ -334,7 +338,16 @@ void XmGridTraceImpl::SetMaxChangeDirectionInRadians(const double a_maxChangeDir m_maxChangeDirectionInRadians = a_maxChangeDirection; } // XmGridTraceImpl::SetMaxChangeDirectionInRadians //------------------------------------------------------------------------------ +/// \brief returns why the last trace operation ended +/// \return the exit reason of the last trace operation +//------------------------------------------------------------------------------ +XmGridTraceExitEnum XmGridTraceImpl::GetExitReason() const +{ + return m_exitReason; +} // XmGridTraceImpl::GetExitReason +//------------------------------------------------------------------------------ /// \brief returns a message describing what caused trace to exit +/// \return the exit message of the last trace operation //------------------------------------------------------------------------------ const std::string& XmGridTraceImpl::GetExitMessage() const { @@ -421,6 +434,7 @@ void XmGridTraceImpl::StepTrace(TraceState& a_state) a_state.m_vy = vy0; a_state.m_mag = mag0; a_state.m_exitReason = a_reason; + m_exitReason = a_reason; m_exitMessage = XmGridTraceExitReasonToString(a_reason); }; @@ -2059,8 +2073,12 @@ void XmGridTraceUnitTests::testBatchMatchesSerialTracePoint() iCreateDefaultSingleCell(serialTracer); std::vector serialTraces(seeds.size()); std::vector serialTimes(seeds.size()); + std::vector serialReasons(seeds.size()); for (size_t i = 0; i < seeds.size(); ++i) + { serialTracer->TracePoint(seeds[i], seedTimes[i], serialTraces[i], serialTimes[i]); + serialReasons[i] = serialTracer->GetExitReason(); + } BSHP batchTracer; iCreateDefaultSingleCell(batchTracer); @@ -2078,6 +2096,9 @@ void XmGridTraceUnitTests::testBatchMatchesSerialTracePoint() { TS_ASSERT_DELTA_VECPT3D(serialTraces[i], batchTraces[i], 1e-12); TS_ASSERT_DELTA_VEC(serialTimes[i], batchTimes[i], 1e-12); + // GetExitReason is the single-point path's answer to what GetTraceResults reports per + // seed; if they can disagree, a caller cannot use TracePoint and the batch interchangeably. + TS_ASSERT_EQUALS((int)reasons[i], (int)serialReasons[i]); // Positions and times are documented as parallel arrays, so a caller may zip them. TS_ASSERT_EQUALS(batchTraces[i].size(), batchTimes[i].size()); } diff --git a/xmsgridtrace/gridtrace/XmGridTrace.h b/xmsgridtrace/gridtrace/XmGridTrace.h index a9c6b24..4d2808c 100644 --- a/xmsgridtrace/gridtrace/XmGridTrace.h +++ b/xmsgridtrace/gridtrace/XmGridTrace.h @@ -186,8 +186,17 @@ class XmGridTrace std::vector& a_outTimes, std::vector& a_outExitReasons) const = 0; + /// \brief Returns why the last trace operation ended. + /// + /// The single-point TracePoint reports through this what GetTraceResults reports per seed. + /// GTEXIT_WAITING_FOR_TIME_STEP means the path stops early because the field is not known + /// past the second loaded time step, not that the particle came to rest -- a distinction + /// TracePoint cannot otherwise express. + /// \return the exit reason of the last trace operation + virtual XmGridTraceExitEnum GetExitReason() const = 0; + /// \brief Returns a human-readable description of what ended the last trace operation. - /// Use GetTraceResults' exit reasons to make decisions; this is for display and logs. + /// Use GetExitReason or GetTraceResults to make decisions; this is for display. /// \return the exit message of the last trace operation virtual const std::string& GetExitMessage() const = 0; diff --git a/xmsgridtrace/python/gridtrace/XmGridTrace_py.cpp b/xmsgridtrace/python/gridtrace/XmGridTrace_py.cpp index 20d0b94..50a55bc 100644 --- a/xmsgridtrace/python/gridtrace/XmGridTrace_py.cpp +++ b/xmsgridtrace/python/gridtrace/XmGridTrace_py.cpp @@ -314,6 +314,119 @@ void initXmGridTrace(py::module &m) { )pydoc"; gridtrace.def("get_exit_message", &xms::XmGridTrace::GetExitMessage, get_exit_message_doc); + // --------------------------------------------------------------------------- + // function: get_exit_reason + // --------------------------------------------------------------------------- + const char* get_exit_reason_doc = R"pydoc( + Returns why the last trace operation ended. + + Prefer this over get_exit_message when deciding what to do with a trace; the message + is for display. WAITING_FOR_TIME_STEP means the path stops early because the field is + not known past the second loaded time step, not that the particle came to rest. + + Returns: + exit_reason_enum: The exit reason of the last trace operation. + )pydoc"; + gridtrace.def("get_exit_reason", &xms::XmGridTrace::GetExitReason, + get_exit_reason_doc); + // --------------------------------------------------------------------------- + // function: start_traces + // --------------------------------------------------------------------------- + const char* start_traces_doc = R"pydoc( + Begins tracing a batch of seeds against the currently loaded time steps. + + A trace runs only as far as the second loaded time step, because that is as far as + the field is known. Supply the next time step with add_grid_scalars_at_time and call + continue_traces to carry every unfinished trace onward:: + + tracer.start_traces(seeds, seed_times) + while tracer.continue_traces() > 0: + step = series.next() + if step is None: + break + tracer.add_grid_scalars_at_time(*step) + traces, times, reasons = tracer.get_trace_results() + + Stopping early is fine: traces still waiting end where they got to. One batch is in + flight per tracer; starting a batch discards any previous one. + + Args: + pts (iterable): The starting point of each trace. + + pt_times (iterable): The starting time of each trace, one per point. + )pydoc"; + gridtrace.def("start_traces", [](xms::XmGridTrace &self, py::iterable pts, + py::iterable pt_times) { + boost::shared_ptr points = xms::VecPt3dFromPyIter(pts); + boost::shared_ptr times = xms::VecDblFromPyIter(pt_times); + if (points->size() != times->size()) + { + // Raised rather than logged: the C++ side refuses the batch and returns empty, + // which from Python would look like a tracer that silently did nothing. + std::string msg = "start_traces needs one start time per point, got " + + std::to_string(points->size()) + " points and " + + std::to_string(times->size()) + " times"; + throw py::value_error(msg); + } + self.StartTraces(*points, *times); + }, start_traces_doc, py::arg("pts"), py::arg("pt_times")); + // --------------------------------------------------------------------------- + // function: continue_traces + // --------------------------------------------------------------------------- + const char* continue_traces_doc = R"pydoc( + Advances every unfinished trace as far as the loaded time steps allow. + + Releases the GIL while tracing, so a caller on a worker thread does not stall the + interpreter. Tracing tens of thousands of seeds takes long enough for that to matter. + + Returns: + int: How many traces are waiting on a later time step. Zero means every trace has + ended for a reason more data cannot change. + )pydoc"; + gridtrace.def("continue_traces", &xms::XmGridTrace::ContinueTraces, + continue_traces_doc, py::call_guard()); + // --------------------------------------------------------------------------- + // function: get_trace_results + // --------------------------------------------------------------------------- + const char* get_trace_results_doc = R"pydoc( + Returns the batch traced so far. + + Valid at any point, complete once continue_traces has returned zero. An entry can hold + fewer than two points: a seed that leaves the grid on its first step yields only the + seed itself, so callers must not assume one usable polyline per seed. + + Returns: + tuple: The positions of each trace, the times of each trace, and why each trace + stopped as an exit_reason_enum. All three are parallel to the seeds passed to + start_traces, and each entry's times are parallel to its positions. + )pydoc"; + gridtrace.def("get_trace_results", [](const xms::XmGridTrace &self) -> py::iterable { + std::vector outTraces; + std::vector outTimes; + std::vector outReasons; + self.GetTraceResults(outTraces, outTimes, outReasons); + py::list traces, times, reasons; + for (size_t i = 0; i < outTraces.size(); ++i) + { + traces.append(xms::PyIterFromVecPt3d(outTraces[i])); + times.append(xms::PyIterFromVecDbl(outTimes[i])); + reasons.append(outReasons[i]); + } + return py::make_tuple(traces, times, reasons); + }, get_trace_results_doc); + + // XmGridTraceExitEnum + py::enum_(m, "exit_reason_enum", + "exit_reason_enum why a trace stopped") + .value("NOT_STARTED", xms::GTEXIT_NOT_STARTED) + .value("WAITING_FOR_TIME_STEP", xms::GTEXIT_WAITING_FOR_TIME_STEP) + .value("MAX_TRACING_TIME", xms::GTEXIT_MAX_TRACING_TIME) + .value("MAX_TRACING_DISTANCE", xms::GTEXIT_MAX_TRACING_DISTANCE) + .value("LEFT_GRID", xms::GTEXIT_LEFT_GRID) + .value("ZERO_VELOCITY", xms::GTEXIT_ZERO_VELOCITY) + .value("MIN_DELTA_TIME", xms::GTEXIT_MIN_DELTA_TIME) + .value("SEED_NOT_TRACEABLE", xms::GTEXIT_SEED_NOT_TRACEABLE) + .value("EXTRACTION_FAILED", xms::GTEXIT_EXTRACTION_FAILED); // DataLocationEnum py::enum_(m, "data_location_enum", From d23108f564f494956e9b8553576f436e8f783b0d Mon Sep 17 00:00:00 2001 From: Bill Dolinar Date: Fri, 14 Aug 2026 10:15:53 -0600 Subject: [PATCH 08/10] Add the AQU104 import group comments to the Python package flake8 reported eight AQU104 findings across the two Python source files, all on line 1: neither carried the numbered import group comments the Aquaveo rules require. They are pre-existing -- four of them are in grid_trace.py, which this branch had not otherwise touched -- and were invisible until flake8-aquaveo was installed. Plain flake8 only runs pycodestyle, pyflakes and mccabe, so the AQU rules, the google docstring convention and the appnexus import order the .flake8 config asks for were all silently unchecked. The convention, matched from next_ms, is all four comments present even when a section is empty, with a blank line after the module docstring. No import moved, so the suite passing is confirmation the modules still resolve the same way. _package is now clean under the full rule set. --- _package/tests/XmGridTrace_pyt.py | 5 +++++ _package/xms/gridtrace/grid_trace.py | 8 ++++++++ 2 files changed, 13 insertions(+) diff --git a/_package/tests/XmGridTrace_pyt.py b/_package/tests/XmGridTrace_pyt.py index 2272bd4..1bbcddd 100644 --- a/_package/tests/XmGridTrace_pyt.py +++ b/_package/tests/XmGridTrace_pyt.py @@ -1,10 +1,15 @@ """Test GridTrace.""" + +# 1. Standard Python modules import unittest +# 2. Third party modules import numpy as np +# 3. Aquaveo modules from xms.grid.ugrid import UGrid +# 4. Local modules from xms.gridtrace import exit_reason_enum, GridTrace diff --git a/_package/xms/gridtrace/grid_trace.py b/_package/xms/gridtrace/grid_trace.py index 3ab03cf..64528e0 100644 --- a/_package/xms/gridtrace/grid_trace.py +++ b/_package/xms/gridtrace/grid_trace.py @@ -1,4 +1,12 @@ """Trace the movement of a point through a velocity vector grid.""" + +# 1. Standard Python modules + +# 2. Third party modules + +# 3. Aquaveo modules + +# 4. Local modules from ._xmsgridtrace import gridtrace From 658dd455734d8b7d8eaa275eb99ff1b52bd73c33 Mon Sep 17 00:00:00 2001 From: Bill Dolinar Date: Fri, 14 Aug 2026 10:22:50 -0600 Subject: [PATCH 09/10] Tier 1: one point-location search per triangulation instead of four per sample GetVectorAtLocationAndTime ran four ExtractData calls per sample -- x and y, for each of two time steps -- and each one performed its own point-location query for the same (x, y). The interpolation weights were identical across all four; only the scalar array being weighted differed. It now runs one search per distinct triangulation and applies the resulting weights to each component. Three things make that possible, and the third is the one that had to be measured rather than assumed: - The x and y extractors of a time step now share a triangulation, via the sharing constructor that has existed since 2022 and was simply never used here. - Both time steps share one as well when their activity masks are equal. The triangulation and its R-tree depend only on the grid and the mask, so an identical mask makes them interchangeable. Differing activity is the case that genuinely cannot share, so the mask is compared rather than assumed -- test_inactive_cell covers that path and test_unique_time_steps covers the shared one. - iApplyWeights reproduces ExtractData exactly, accumulating in double and narrowing to float, so this is bit-identical rather than merely close. Every recorded baseline in both the C++ and Python suites is unchanged, which is the evidence for that. Ordering matters in AddGridScalarsAtTime and is easy to get backwards: the y extractor is built from x, and only after x's scalars are set. The sharing constructor copies the triangulation and the flag saying what it was built for, so copying x before it has built one leaves y believing it must build, and y then rebuilds the very triangulation it is sharing -- silently costing what the sharing was meant to save. Measured on a 200x200 grid at 10,000 seeds, against the previous commit: before after setup, 2 time steps 117 ms 32 ms searches per seed 97.5 24.4 interior us/seed 50.8 12.6 mixed us/seed 50.7 11.8 The setup figure is better than the ~53 ms projected in TRIANGULATION_SHARING.md. That projection assumed a per-extractor scalar-handling floor of roughly 8 ms, derived by subtraction rather than timed; the measurement says it is far smaller, and that only one triangulation and one R-tree are now built for all four extractors rather than four of each. A re-trace of 10,000 glyphs is now about 150 ms including setup, against the 299 ms the current in-render tracer costs for the same work -- so the new path is faster than what it replaces while adding the time dimension, rather than merely close enough. It began this branch at about 21 seconds. Also fixes a latent crash this rewrite made obvious: with only one time step supplied, the first extractor is null and was dereferenced. It now returns a clean extraction failure, which is what edge case 9 in the session plan assumed already happened. The benchmark's counter counted ExtractData calls and now counts searches, which is what it always meant; its label and field are renamed to match rather than silently changing meaning. --- xmsgridtrace/gridtrace/XmGridTrace.cpp | 173 +++++++++++++++++-------- xmsgridtrace/gridtrace/XmGridTrace.h | 1 - 2 files changed, 121 insertions(+), 53 deletions(-) diff --git a/xmsgridtrace/gridtrace/XmGridTrace.cpp b/xmsgridtrace/gridtrace/XmGridTrace.cpp index 78dfce5..df6227a 100644 --- a/xmsgridtrace/gridtrace/XmGridTrace.cpp +++ b/xmsgridtrace/gridtrace/XmGridTrace.cpp @@ -24,6 +24,7 @@ #include // XM_ZERO_TOL #include #include +#include #include // 6. Non-shared code headers @@ -46,14 +47,14 @@ namespace /// XMS Namespace #ifdef CXX_TEST -/// \brief Count of XmUGrid2dDataExtractor::ExtractData calls since it was last zeroed. +/// \brief Count of point-location searches since it was last zeroed. /// Test-build-only instrumentation for testTraceBenchmark. A trace's cost is dominated by -/// the point-location search each ExtractData performs, so the benchmark needs the search -/// count and not only wall time -- otherwise an algorithmic win cannot be told apart from -/// a faster machine. Not thread safe; the benchmark is single threaded. -size_t g_extractDataCalls = 0; -/// \brief Adds a_n to the ExtractData call count. Compiles away outside test builds. -#define XMGT_COUNT_EXTRACT_DATA(a_n) (g_extractDataCalls += (a_n)) +/// these searches, so the benchmark needs the count and not only wall time -- otherwise an +/// algorithmic win cannot be told apart from a faster machine. Not thread safe; the +/// benchmark is single threaded. +size_t g_searchCalls = 0; +/// \brief Adds a_n to the search count. Compiles away outside test builds. +#define XMGT_COUNT_SEARCH(a_n) (g_searchCalls += (a_n)) /// \brief Count of XmUGrid2dPolylineDataExtractor constructions since it was last zeroed. /// Test-build-only instrumentation for testBoundaryExtractorIsCached. Caching that extractor /// is a pure performance change with no effect on trace output, so a construction count is @@ -63,7 +64,7 @@ size_t g_boundaryExtractorBuilds = 0; #define XMGT_COUNT_BOUNDARY_EXTRACTOR_BUILD() (++g_boundaryExtractorBuilds) #else /// \brief No-op outside test builds, so production traces pay nothing for instrumentation. -#define XMGT_COUNT_EXTRACT_DATA(a_n) ((void)0) +#define XMGT_COUNT_SEARCH(a_n) ((void)0) /// \brief No-op outside test builds. #define XMGT_COUNT_BOUNDARY_EXTRACTOR_BUILD() ((void)0) #endif @@ -80,6 +81,40 @@ bool iIsTerminal(XmGridTraceExitEnum a_reason) return a_reason != GTEXIT_NOT_STARTED && a_reason != GTEXIT_WAITING_FOR_TIME_STEP; } // iIsTerminal +//------------------------------------------------------------------------------ +/// \brief Applies one set of interpolation weights to a time step's x and y scalars. +/// +/// Reproduces XmUGrid2dDataExtractor::ExtractData exactly -- accumulating in double, then +/// narrowing to float -- so that replacing four ExtractData calls with one search plus this +/// gives bit-identical answers rather than merely close ones. +/// \param[in] a_x The extractor holding the x component +/// \param[in] a_y The extractor holding the y component, sharing a_x's triangulation +/// \param[in] a_idxs Triangulation point indices from the search +/// \param[in] a_weights Interpolation weights parallel to a_idxs +/// \param[out] a_outX The interpolated x component +/// \param[out] a_outY The interpolated y component +//------------------------------------------------------------------------------ +void iApplyWeights(const XmUGrid2dDataExtractor& a_x, + const XmUGrid2dDataExtractor& a_y, + const VecInt& a_idxs, + const VecDbl& a_weights, + float& a_outX, + float& a_outY) +{ + const VecFlt& xScalars = a_x.GetScalars(); + const VecFlt& yScalars = a_y.GetScalars(); + double interpX = 0.0, interpY = 0.0; + for (size_t i = 0; i < a_idxs.size(); ++i) + { + const int ptIdx = a_idxs[i]; + const double weight = a_weights[i]; + interpX += xScalars[ptIdx] * weight; + interpY += yScalars[ptIdx] * weight; + } + a_outX = static_cast(interpX); + a_outY = static_cast(interpY); +} // iApplyWeights + //////////////////////////////////////////////////////////////////////////////// /// One trace in progress, and everything about it that has to survive a time step change. /// @@ -184,6 +219,16 @@ class XmGridTraceImpl : public XmGridTrace /// data extractor for the y component for the second time step BSHP m_extractor2y; double m_time2=-1; ///< time of the second time step + xms::DynBitset m_activity2; ///< activity of the second time step, to compare with the next + /// Whether both time steps share one triangulation, which they can when their activity + /// matches. When they do, one search serves all four extractors instead of one per step. + bool m_sharedAcrossTime = false; + /// Scratch for the point-location search. Members rather than locals because + /// GetVectorAtLocationAndTime runs a few dozen times per traced seed and these would + /// otherwise reallocate on every call. They make the tracer unsafe to share across + /// threads, which it already was -- GmTriSearch caches barycentric state per query. + mutable VecInt m_searchIdxs; + mutable VecDbl m_searchWeights; /// Extractor used to find where a trace leaves the grid, built lazily on the first /// out-of-domain step and reused for every one after it. Its construction triangulates the /// whole grid and its first SetPolyline indexes every triangle into a GmMultiPolyIntersector; @@ -369,32 +414,49 @@ void XmGridTraceImpl::AddGridScalarsAtTime(const VecPt3d& a_scalars, DataLocationEnum a_activityLoc, double a_time) { - if (m_extractor2x && m_extractor2y) + const bool hadPrevious = m_extractor2x && m_extractor2y; + if (hadPrevious) { m_extractor1x = m_extractor2x; m_extractor1y = m_extractor2y; m_time1 = m_time2; } - m_extractor2x = XmUGrid2dDataExtractor::New(m_ugrid); - m_extractor2y = XmUGrid2dDataExtractor::New(m_ugrid); m_time2 = a_time; std::vector xx, yy; + xx.reserve(a_scalars.size()); + yy.reserve(a_scalars.size()); for (auto& pt : a_scalars) { xx.push_back((float)pt.x); yy.push_back((float)pt.y); } + + // Share the triangulation with the previous time step when the two agree on activity. The + // triangulation and its GmTriSearch R-tree depend only on the grid and the activity mask, + // so an identical mask makes them interchangeable -- which both skips a rebuild and lets a + // single point-location query serve all four extractors instead of one per time step. + // Differing activity is the case that cannot share, and it is why the mask is compared + // rather than assumed. + m_sharedAcrossTime = hadPrevious && a_activity == m_activity2; + m_extractor2x = m_sharedAcrossTime ? XmUGrid2dDataExtractor::New(m_extractor1x) + : XmUGrid2dDataExtractor::New(m_ugrid); if (a_scalarLoc == DataLocationEnum::LOC_POINTS) - { m_extractor2x->SetGridPointScalars(xx, a_activity, a_activityLoc); - m_extractor2y->SetGridPointScalars(yy, a_activity, a_activityLoc); - } else - { m_extractor2x->SetGridCellScalars(xx, a_activity, a_activityLoc); + + // y is built from x, and only after x's scalars are set. The sharing constructor copies + // the triangulation *and* the flag saying what it was built for; copying x before it has + // built one would leave y thinking it must build, and y would then rebuild the very + // triangulation it is sharing. Only the scalar arrays differ between the two. + m_extractor2y = XmUGrid2dDataExtractor::New(m_extractor2x); + if (a_scalarLoc == DataLocationEnum::LOC_POINTS) + m_extractor2y->SetGridPointScalars(yy, a_activity, a_activityLoc); + else m_extractor2y->SetGridCellScalars(yy, a_activity, a_activityLoc); - } + + m_activity2 = a_activity; } //------------------------------------------------------------------------------ @@ -721,32 +783,40 @@ bool XmGridTraceImpl::GetVectorAtLocationAndTime(const xms::Pt3d& a_pt, double a_currentTime, xms::Pt3d& a_data) const { - xms::VecPt3d loc; - loc.push_back(a_pt); - m_extractor1x->SetExtractLocations(loc); - m_extractor1y->SetExtractLocations(loc); - xms::VecFlt dataOutx1; - xms::VecFlt dataOuty1; - m_extractor1x->ExtractData(dataOutx1); - m_extractor1y->ExtractData(dataOuty1); - XMGT_COUNT_EXTRACT_DATA(2); - if (dataOutx1.size() != 1 || dataOuty1.size() != 1) + if (!m_extractor1x || !m_extractor1y || !m_extractor2x || !m_extractor2y) { - XM_LOG(xmlog::error, "Gridtracer: An error occured when extracting data"); + // Two time steps are required. This used to dereference a null first extractor when only + // one had been supplied. + XM_LOG(xmlog::error, "Gridtracer: two time steps must be added before tracing."); return false; } - m_extractor2x->SetExtractLocations(loc); - m_extractor2y->SetExtractLocations(loc); - xms::VecFlt dataOutx2; - xms::VecFlt dataOuty2; - m_extractor2x->ExtractData(dataOutx2); - m_extractor2y->ExtractData(dataOuty2); - XMGT_COUNT_EXTRACT_DATA(2); - if (dataOutx2.size() != 1 || dataOuty2.size() != 1) + // One point-location query per distinct triangulation, rather than one per scalar array. + // The weights returned index the triangulation's points, and every extractor sharing that + // triangulation indexes its own scalars the same way, so a single query serves the x and y + // of a time step -- and both time steps too when they share a triangulation. + float x1 = m_extractor1x->GetNoDataValue(); + float y1 = m_extractor1y->GetNoDataValue(); + const int cell1 = + m_extractor1x->GetUGridTriangles()->GetIntersectedCell(a_pt, m_searchIdxs, m_searchWeights); + XMGT_COUNT_SEARCH(1); + if (cell1 >= 0) + iApplyWeights(*m_extractor1x, *m_extractor1y, m_searchIdxs, m_searchWeights, x1, y1); + + float x2 = m_extractor2x->GetNoDataValue(); + float y2 = m_extractor2y->GetNoDataValue(); + if (m_sharedAcrossTime) { - XM_LOG(xmlog::error, "Gridtracer: An error occured when extracting data"); - return false; + if (cell1 >= 0) + iApplyWeights(*m_extractor2x, *m_extractor2y, m_searchIdxs, m_searchWeights, x2, y2); + } + else + { + const int cell2 = + m_extractor2x->GetUGridTriangles()->GetIntersectedCell(a_pt, m_searchIdxs, m_searchWeights); + XMGT_COUNT_SEARCH(1); + if (cell2 >= 0) + iApplyWeights(*m_extractor2x, *m_extractor2y, m_searchIdxs, m_searchWeights, x2, y2); } if (a_currentTime < m_time1 - XM_ZERO_TOL) @@ -759,8 +829,8 @@ bool XmGridTraceImpl::GetVectorAtLocationAndTime(const xms::Pt3d& a_pt, // XM_NODATA (-9999999) against a real value produces something like -999999.9, which is // neither no-data nor meaningful, and every caller tests for XM_NODATA exactly. Returning // true is correct -- extraction succeeded, and no-data is the answer. - if (EQ_TOL(dataOutx1[0], XM_NODATA, 1) || EQ_TOL(dataOuty1[0], XM_NODATA, 1) || - EQ_TOL(dataOutx2[0], XM_NODATA, 1) || EQ_TOL(dataOuty2[0], XM_NODATA, 1)) + if (EQ_TOL(x1, XM_NODATA, 1) || EQ_TOL(y1, XM_NODATA, 1) || EQ_TOL(x2, XM_NODATA, 1) || + EQ_TOL(y2, XM_NODATA, 1)) { a_data.x = XM_NODATA; a_data.y = XM_NODATA; @@ -775,8 +845,8 @@ bool XmGridTraceImpl::GetVectorAtLocationAndTime(const xms::Pt3d& a_pt, // particle released at m_time1 entirely by the field at m_time2. double weight1 = fabs(a_currentTime - m_time2) / totalTime; double weight2 = fabs(a_currentTime - m_time1) / totalTime; - a_data.x = dataOutx1[0] * weight1 + dataOutx2[0] * weight2; - a_data.y = dataOuty1[0] * weight1 + dataOuty2[0] * weight2; + a_data.x = x1 * weight1 + x2 * weight2; + a_data.y = y1 * weight1 + y2 * weight2; return true; } // XmGridTraceImpl::GetVectorAtLocationAndTime } // namespace {} @@ -987,7 +1057,7 @@ struct BenchmarkStats int m_seeds = 0; ///< seed points handed to TracePoint int m_traced = 0; ///< seeds that produced a usable (2+ point) polyline size_t m_tracePoints = 0; ///< total polyline points produced - size_t m_extractCalls = 0; ///< XmUGrid2dDataExtractor::ExtractData calls consumed + size_t m_searchCalls = 0; ///< point-location searches consumed double m_seconds = 0; ///< wall time of the traced batch, excluding setup std::map m_exitReasons; ///< exit message -> count, over a sample }; @@ -1100,7 +1170,7 @@ void iRunTraceBenchmark(BSHP& a_tracer, VecPt3d trace; VecDbl times; - g_extractDataCalls = 0; + g_searchCalls = 0; const auto start = std::chrono::steady_clock::now(); for (const auto& seed : a_seeds) { @@ -1113,7 +1183,7 @@ void iRunTraceBenchmark(BSHP& a_tracer, } const auto end = std::chrono::steady_clock::now(); a_stats.m_seconds = std::chrono::duration(end - start).count(); - a_stats.m_extractCalls = g_extractDataCalls; + a_stats.m_searchCalls = g_searchCalls; const int sampleSize = std::min((int)a_seeds.size(), 1000); for (int i = 0; i < sampleSize; ++i) @@ -1131,9 +1201,9 @@ void iReportTraceBenchmark(const char* a_label, const BenchmarkStats& a_stats) { const double seeds = a_stats.m_seeds ? (double)a_stats.m_seeds : 1.0; const double usPerSeed = a_stats.m_seconds * 1e6 / seeds; - const double extractsPerSeed = a_stats.m_extractCalls / seeds; + const double searchesPerSeed = a_stats.m_searchCalls / seeds; const double usPerExtract = - a_stats.m_extractCalls ? a_stats.m_seconds * 1e6 / a_stats.m_extractCalls : 0.0; + a_stats.m_searchCalls ? a_stats.m_seconds * 1e6 / a_stats.m_searchCalls : 0.0; const double ptsPerTrace = a_stats.m_traced ? (double)a_stats.m_tracePoints / a_stats.m_traced : 0.0; @@ -1141,9 +1211,8 @@ void iReportTraceBenchmark(const char* a_label, const BenchmarkStats& a_stats) << "] seeds=" << a_stats.m_seeds << " traced=" << a_stats.m_traced << "\n" << " wall " << a_stats.m_seconds * 1e3 << " ms\n" << " per seed " << usPerSeed << " us\n" - << " ExtractData " << a_stats.m_extractCalls << " calls (" - << std::setprecision(1) << extractsPerSeed << "/seed, " << std::setprecision(3) - << usPerExtract << " us/call)\n" + << " searches " << a_stats.m_searchCalls << " (" << std::setprecision(1) + << searchesPerSeed << "/seed, " << std::setprecision(3) << usPerExtract << " us/call)\n" << " trace points " << a_stats.m_tracePoints << " (" << std::setprecision(1) << ptsPerTrace << "/trace)\n" << " exit reasons (sampled):\n"; @@ -2250,14 +2319,14 @@ void XmGridTraceUnitTests::testBoundaryExtractorIsCached() /// measured separately because they exercise different code: /// /// interior seeds far enough from the edge that no trace can reach it -- the pure -/// stepping cost, four ExtractData searches per integration step +/// stepping cost, one point-location search per triangulation per step /// boundary seeds in a band along the edge, so traces run out of the domain and pay for /// the XmUGrid2dPolylineDataExtractor path -- a whole-grid triangulation plus a /// GmMultiPolyIntersector, once per tracer since that extractor is cached /// (it was once per exit event, inside the stepping loop) /// mixed seeds spread over the whole domain -- what the display actually does /// -/// Reported alongside wall time is the ExtractData call count, so a later optimization +/// Reported alongside wall time is the point-location search count, so a later optimization /// can be shown to have removed searches rather than merely found a faster machine. /// /// Seed count and grid size come from XMGT_BENCH_SEEDS and XMGT_BENCH_CELLS so a sweep @@ -2362,7 +2431,7 @@ void XmGridTraceUnitTests::testTraceBenchmark() // bound tight enough that a real breakage in tracing still fails here. TS_ASSERT(mixed.m_traced >= seedCount - 1 - seedCount / 1000); // The instrumentation itself has to be working, or the search counts mean nothing. - TS_ASSERT(interior.m_extractCalls > (size_t)seedCount); + TS_ASSERT(interior.m_searchCalls > (size_t)seedCount); // The boundary set must actually leave the grid, otherwise this benchmark silently // stops measuring the per-exit extractor construction it exists to measure. const std::string outOfDomain = "Point has traveled out of domain."; diff --git a/xmsgridtrace/gridtrace/XmGridTrace.h b/xmsgridtrace/gridtrace/XmGridTrace.h index 4d2808c..999f21b 100644 --- a/xmsgridtrace/gridtrace/XmGridTrace.h +++ b/xmsgridtrace/gridtrace/XmGridTrace.h @@ -52,7 +52,6 @@ enum XmGridTraceExitEnum { GTEXIT_EXTRACTION_FAILED ///< a field lookup failed; the trace is discarded }; - //----- Structs / Classes ------------------------------------------------------ //////////////////////////////////////////////////////////////////////////////// From 143cfd7e51aaae2429798bc5b32900b8b6444bf0 Mon Sep 17 00:00:00 2001 From: Bill Dolinar Date: Fri, 14 Aug 2026 13:06:53 -0600 Subject: [PATCH 10/10] Fix three defects found reviewing the batch tracing API A resumed trace could hang, a legal pair of time steps could read out of bounds, and a staggered seed could be killed off permanently. All three are in code this branch introduced. StepTrace hung when a window ended exactly on the second time step. The time step clamp computes deltaT = m_time2 - elapsed - ptTime, which for a trace already sitting on m_time2 is exactly zero, and that zero was persisted into the resumed trace. A zero-length step moves nothing and changes no velocity, so no clamp and no subdivision test could ever end the loop -- and the min-delta-time escape is inside the split branch, which a zero-length step can never enter. It spun forever appending nothing, with the GIL released so Python could not interrupt it. The window is now finished before stepping, keeping the step size the call came in with, so a redundant ContinueTraces really is the no-op the header promises. StepTrace also floors a non-positive resumed step size, so no path can reintroduce this. AddGridScalarsAtTime decided triangulation sharing from the activity mask alone. The triangulation is built for a data location -- LOC_CELLS adds a centroid per cell, LOC_POINTS adds none -- and sharing shares the object rather than copying it, so the second step's SetGrid*Scalars rebuilt the triangulation the first step was still pointing at. Its shorter scalar array was then indexed by the new centroid indices: an out-of-bounds read, not a wrong answer. Both data locations now join the mask in the predicate. A seed released after the loaded window reported GTEXIT_EXTRACTION_FAILED, which iIsTerminal treats as terminal, so the seed never started even once its time step arrived. StartTraces takes a release time per seed so a batch can be staggered, making this an ordinary input; it now reports GTEXIT_WAITING_FOR_TIME_STEP, as the mid-trace clamp always did. Adds a regression test per fix. testBeyondTimestep and its Python twin asserted only that the trace was empty, which is why the third defect went unnoticed -- they now assert the exit reason, which is what tells the three empty-trace cases apart. C++ 25/25, Python 19/19, flake8 clean. --- _package/tests/XmGridTrace_pyt.py | 7 +- xmsgridtrace/gridtrace/XmGridTrace.cpp | 259 +++++++++++++++++++++++-- xmsgridtrace/gridtrace/XmGridTrace.h | 4 + xmsgridtrace/gridtrace/XmGridTrace.t.h | 3 + 4 files changed, 260 insertions(+), 13 deletions(-) diff --git a/_package/tests/XmGridTrace_pyt.py b/_package/tests/XmGridTrace_pyt.py index 1bbcddd..01ba8a4 100644 --- a/_package/tests/XmGridTrace_pyt.py +++ b/_package/tests/XmGridTrace_pyt.py @@ -387,9 +387,10 @@ def test_start_out_of_cell(self): expected_out_times = [] np.testing.assert_equal(0, len(result_tuple[0])) np.testing.assert_array_almost_equal(expected_out_times, result_tuple[1]) + self.assertEqual(exit_reason_enum.SEED_NOT_TRACEABLE, tracer.get_exit_reason()) def test_beyond_timestep(self): - """Test functionality of starting beyond the time step.""" + """Test that a start time past the loaded window waits rather than failing.""" tracer = self.create_default_single_cell() start_time = 10.1 @@ -398,6 +399,10 @@ def test_beyond_timestep(self): expected_out_times = [] np.testing.assert_equal(0, len(result_tuple[0])) np.testing.assert_array_almost_equal(expected_out_times, result_tuple[1]) + # This and test_start_out_of_cell both produce an empty trace, so emptiness alone + # cannot tell them apart -- which is how this case went unnoticed as an extraction + # failure. The field is not known this far ahead yet; the trace is waiting for data. + self.assertEqual(exit_reason_enum.WAITING_FOR_TIME_STEP, tracer.get_exit_reason()) def test_before_timestep(self): """Test functionality of starting before the time step.""" diff --git a/xmsgridtrace/gridtrace/XmGridTrace.cpp b/xmsgridtrace/gridtrace/XmGridTrace.cpp index df6227a..9bf761e 100644 --- a/xmsgridtrace/gridtrace/XmGridTrace.cpp +++ b/xmsgridtrace/gridtrace/XmGridTrace.cpp @@ -71,6 +71,11 @@ size_t g_boundaryExtractorBuilds = 0; //----- Class / Function definitions ------------------------------------------- +/// Step size a trace begins with, and the value a resumed trace falls back to when the +/// window it just finished clamped its step to zero. See StepTrace for why zero cannot be +/// carried forward. +const double kInitialDeltaT = 1.0; + //------------------------------------------------------------------------------ /// \brief Whether a reason means the trace can never advance again. /// \param[in] a_reason The exit reason @@ -131,7 +136,7 @@ struct TraceState double m_ptTime = 0; ///< time the trace was released; never advanced double m_elapsedTime = 0; ///< time advanced since release, against m_maxTracingTime double m_distTraveled = 0; ///< distance covered, against m_maxTracingDistance - double m_deltaT = 1.0; ///< adaptive step size carried into the next step + double m_deltaT = kInitialDeltaT; ///< adaptive step size carried into the next step double m_vx = 0; ///< velocity x at m_pt, for the subdivision tests double m_vy = 0; ///< velocity y at m_pt, for the subdivision tests double m_mag = 0; ///< speed at m_pt, for the change-in-velocity test @@ -220,8 +225,15 @@ class XmGridTraceImpl : public XmGridTrace BSHP m_extractor2y; double m_time2=-1; ///< time of the second time step xms::DynBitset m_activity2; ///< activity of the second time step, to compare with the next - /// Whether both time steps share one triangulation, which they can when their activity - /// matches. When they do, one search serves all four extractors instead of one per step. + /// Data location of the second time step's scalars, to compare with the next. The + /// triangulation is built for a location, so a change here forbids sharing. + DataLocationEnum m_scalarLoc2 = DataLocationEnum::LOC_UNKNOWN; + /// Data location of the second time step's activity, to compare with the next. Decides how + /// the activity bitset maps onto cell activity, so a change here forbids sharing too. + DataLocationEnum m_activityLoc2 = DataLocationEnum::LOC_UNKNOWN; + /// Whether both time steps share one triangulation, which they can when the two steps agree + /// on activity and on both data locations. When they do, one search serves all four + /// extractors instead of one per step. bool m_sharedAcrossTime = false; /// Scratch for the point-location search. Members rather than locals because /// GetVectorAtLocationAndTime runs a few dozen times per traced seed and these would @@ -432,13 +444,21 @@ void XmGridTraceImpl::AddGridScalarsAtTime(const VecPt3d& a_scalars, yy.push_back((float)pt.y); } - // Share the triangulation with the previous time step when the two agree on activity. The - // triangulation and its GmTriSearch R-tree depend only on the grid and the activity mask, - // so an identical mask makes them interchangeable -- which both skips a rebuild and lets a - // single point-location query serve all four extractors instead of one per time step. - // Differing activity is the case that cannot share, and it is why the mask is compared - // rather than assumed. - m_sharedAcrossTime = hadPrevious && a_activity == m_activity2; + // Share the triangulation with the previous time step when the two agree on everything it + // is built from: the grid (fixed at construction), the data location, and the activity + // mask. When they do, one point-location query serves all four extractors instead of one + // per time step, and no rebuild happens. + // + // All three terms are load-bearing, and the location ones are the easy ones to miss. The + // triangulation's shape comes from a_scalarLoc -- LOC_CELLS adds a centroid point per cell + // and LOC_POINTS adds none -- while a_activityLoc decides how the same bitset maps onto + // cell activity. Sharing does not copy the triangulation, it shares the object, and the + // second step's SetGrid*Scalars rebuilds that shared object in place; so sharing across a + // location change would rebuild the triangulation the *first* step is still pointing at, + // leaving its shorter scalar array indexed by the new triangulation's centroid indices. + // That is an out-of-bounds read in iApplyWeights, not a wrong answer. + m_sharedAcrossTime = hadPrevious && a_activity == m_activity2 && + a_scalarLoc == m_scalarLoc2 && a_activityLoc == m_activityLoc2; m_extractor2x = m_sharedAcrossTime ? XmUGrid2dDataExtractor::New(m_extractor1x) : XmUGrid2dDataExtractor::New(m_ugrid); if (a_scalarLoc == DataLocationEnum::LOC_POINTS) @@ -457,6 +477,8 @@ void XmGridTraceImpl::AddGridScalarsAtTime(const VecPt3d& a_scalars, m_extractor2y->SetGridCellScalars(yy, a_activity, a_activityLoc); m_activity2 = a_activity; + m_scalarLoc2 = a_scalarLoc; + m_activityLoc2 = a_activityLoc; } //------------------------------------------------------------------------------ @@ -475,6 +497,15 @@ void XmGridTraceImpl::StepTrace(TraceState& a_state) const double ptTime = a_state.m_ptTime; Pt3d pt0 = a_state.m_pt, pt1; double deltaT = a_state.m_deltaT; + // A window that ended exactly on m_time2 left deltaT clamped to zero (see the time step + // clamp in the loop below), and zero cannot be carried into the next window: a zero-length + // step moves nothing and changes no velocity, so it satisfies none of the loop's exit + // tests -- not the clamps, which need elapsedTime to advance, and not the subdivision + // tests, which compare a step against the one before it and would see no change. The loop + // would spin forever. Start the next window from the initial step and let the clamps size + // it again, which is what a fresh trace does. + if (deltaT <= 0) + deltaT = kInitialDeltaT; double elapsedTime = a_state.m_elapsedTime; double distTraveled = a_state.m_distTraveled; double vx0 = a_state.m_vx, vy0 = a_state.m_vy, mag0 = a_state.m_mag; @@ -504,8 +535,18 @@ void XmGridTraceImpl::StepTrace(TraceState& a_state) { outTrace.clear(); outTimes.clear(); - if (ptTime > m_time2 || // Test if the time specified is after the time range - !GetVectorAtLocationAndTime(pt0, ptTime, vector)) // Ensure extraction did not fail + if (ptTime > m_time2) + { + // The seed is released after the loaded window, so its field is not known yet. That is + // the same situation the time step clamp below reports as WAITING, and it has to be + // reported the same way here: EXTRACTION_FAILED is terminal (see iIsTerminal), so a + // seed given a later release time than the current window would never start, even once + // the time step covering it arrived. StartTraces takes a release time per seed + // precisely so a batch can be staggered, which makes this a normal input, not an error. + stopWith(GTEXIT_WAITING_FOR_TIME_STEP); + return; + } + if (!GetVectorAtLocationAndTime(pt0, ptTime, vector)) // Ensure extraction did not fail { stopWith(GTEXIT_EXTRACTION_FAILED); return; @@ -546,6 +587,18 @@ void XmGridTraceImpl::StepTrace(TraceState& a_state) if (elapsedTime + deltaT + ptTime > m_time2) { deltaT = m_time2 - elapsedTime - ptTime; + if (deltaT <= 0) + { + // Nothing left in this window -- the trace is already sitting exactly on m_time2, + // which is what a second ContinueTraces with no new data finds. Stop before stepping, + // and put back the step size this call came in with: a zero-length step would append + // nothing anyway, and persisting the zero is what used to leave the resumed trace + // unable to advance at all. Restoring it is what makes a redundant ContinueTraces + // genuinely do no useful work, rather than quietly changing the path that follows. + deltaT = a_state.m_deltaT; + stopWith(GTEXIT_WAITING_FOR_TIME_STEP); + return; + } bContinue = false; // This will be the last point traced in this window stopReason = GTEXIT_WAITING_FOR_TIME_STEP; } @@ -1610,6 +1663,10 @@ void XmGridTraceUnitTests::testBeyondTimestep() VecDbl expectedOutTimes = {}; TS_ASSERT_DELTA_VECPT3D(expectedOutTrace, outTrace, .0001); TS_ASSERT_DELTA_VEC(expectedOutTimes, outTimes, .0001); + // An empty trace on its own does not say which of several unrelated things happened, which + // is how this case went unnoticed as an extraction failure. The field simply is not known + // this far ahead yet, so the trace is waiting -- supplying a later time step starts it. + TS_ASSERT_EQUALS((int)GTEXIT_WAITING_FOR_TIME_STEP, (int)tracer->GetExitReason()); } // XmGridTraceUnitTests::testBeyondTimestep //------------------------------------------------------------------------------ /// \brief test the behavior when starting before the first timestep @@ -2273,6 +2330,184 @@ void XmGridTraceUnitTests::testTracesContinueAcrossTimeSteps() TS_ASSERT(traces[0].back().x < maxX); } // XmGridTraceUnitTests::testTracesContinueAcrossTimeSteps //------------------------------------------------------------------------------ +/// \brief Returns a tracer whose spatially uniform field rotates +x -> +y across two steps. +/// +/// One cell spanning the domain, so the field is uniform in space and every change in a path +/// comes from time. Steps at t = 0 (east) and t = 10 (north) are loaded; supply a third to +/// let a trace resume past t = 10. +/// \param[out] a_activity Single-cell activity, for supplying further time steps +/// \return the tracer +//------------------------------------------------------------------------------ +BSHP iCreateRotatingFieldTracer(DynBitset& a_activity) +{ + VecPt3d points = {{0, 0, 0}, {40, 0, 0}, {40, 40, 0}, {0, 40, 0}}; + VecInt cells = {XMU_QUAD, 4, 0, 1, 2, 3}; + std::shared_ptr ugrid = XmUGrid::New(points, cells); + + a_activity.clear(); + a_activity.push_back(true); + + BSHP tracer = XmGridTrace::New(ugrid); + tracer->SetVectorMultiplier(1); + tracer->SetMaxTracingTime(18); + tracer->SetMaxTracingDistance(1000); + tracer->SetMinDeltaTime(.01); + tracer->SetMaxChangeDistance(.5); + tracer->SetMaxChangeVelocity(-1); + tracer->SetMaxChangeDirectionInRadians(XM_PI); // never subdivide on direction + VecPt3d east = {{1, 0, 0}}, north = {{0, 1, 0}}; + tracer->AddGridScalarsAtTime(east, DataLocationEnum::LOC_CELLS, a_activity, + DataLocationEnum::LOC_CELLS, 0.0); + tracer->AddGridScalarsAtTime(north, DataLocationEnum::LOC_CELLS, a_activity, + DataLocationEnum::LOC_CELLS, 10.0); + return tracer; +} // iCreateRotatingFieldTracer +//------------------------------------------------------------------------------ +/// \brief A redundant ContinueTraces must not change what the trace does afterwards. +/// +/// XmGridTrace.h sanctions calling ContinueTraces twice with no time step in between, saying +/// it does no useful work. It used to do considerably worse than nothing: the first call ends +/// a window by clamping deltaT to exactly m_time2 - elapsed - ptTime, which for a trace +/// already sitting on m_time2 is exactly zero, and that zero was carried into the resumed +/// trace. A zero-length step moves nothing and changes no velocity, so no clamp and no +/// subdivision test could ever end the loop -- it spun forever, appending nothing, with the +/// GIL released so Python could not interrupt it. +/// +/// The assertion is equality with the run that did not make the redundant call. "Does no +/// useful work" is only true if the outcome is indistinguishable. +//------------------------------------------------------------------------------ +void XmGridTraceUnitTests::testRedundantContinueDoesNotStallTrace() +{ + const VecPt3d seeds = {{20, 10, 0}}; + const VecDbl seedTimes = {0}; + VecPt3d west = {{-1, 0, 0}}; + + DynBitset plainActivity; + BSHP plain = iCreateRotatingFieldTracer(plainActivity); + plain->StartTraces(seeds, seedTimes); + TS_ASSERT_EQUALS(1, plain->ContinueTraces()); + plain->AddGridScalarsAtTime(west, DataLocationEnum::LOC_CELLS, plainActivity, + DataLocationEnum::LOC_CELLS, 20.0); + TS_ASSERT_EQUALS(0, plain->ContinueTraces()); + std::vector plainTraces; + std::vector plainTimes; + std::vector plainReasons; + plain->GetTraceResults(plainTraces, plainTimes, plainReasons); + + DynBitset activity; + BSHP tracer = iCreateRotatingFieldTracer(activity); + tracer->StartTraces(seeds, seedTimes); + TS_ASSERT_EQUALS(1, tracer->ContinueTraces()); + // The redundant call. Still waiting, because no new data arrived. + TS_ASSERT_EQUALS(1, tracer->ContinueTraces()); + tracer->AddGridScalarsAtTime(west, DataLocationEnum::LOC_CELLS, activity, + DataLocationEnum::LOC_CELLS, 20.0); + // Before the fix this call never returned. + TS_ASSERT_EQUALS(0, tracer->ContinueTraces()); + std::vector traces; + std::vector times; + std::vector reasons; + tracer->GetTraceResults(traces, times, reasons); + + TS_ASSERT_EQUALS((int)plainReasons[0], (int)reasons[0]); + TS_ASSERT_EQUALS(plainTraces[0].size(), traces[0].size()); + TS_ASSERT_DELTA_VECPT3D(plainTraces[0], traces[0], 1e-12); + TS_ASSERT_DELTA_VEC(plainTimes[0], times[0], 1e-12); +} // XmGridTraceUnitTests::testRedundantContinueDoesNotStallTrace +//------------------------------------------------------------------------------ +/// \brief A seed released after the loaded window waits for its data instead of failing. +/// +/// StartTraces takes a release time per seed so a batch can be staggered, which makes a seed +/// timed past the current window an ordinary input. It used to be reported as +/// GTEXIT_EXTRACTION_FAILED, which iIsTerminal treats as terminal, so the seed was dead: the +/// time step covering it could arrive and ContinueTraces would never look at it again. The +/// mid-trace clamp had always called this same condition WAITING; only the seed path did not. +//------------------------------------------------------------------------------ +void XmGridTraceUnitTests::testSeedReleasedAfterWindowWaitsThenTraces() +{ + DynBitset activity; + BSHP tracer = iCreateRotatingFieldTracer(activity); + + // Steps at t = 0 and t = 10 are loaded; this seed is released at 15. + const VecPt3d seeds = {{20, 10, 0}}; + tracer->StartTraces(seeds, {15}); + TS_ASSERT_EQUALS(1, tracer->ContinueTraces()); // waiting, not failed + + std::vector traces; + std::vector times; + std::vector reasons; + tracer->GetTraceResults(traces, times, reasons); + TS_ASSERT_EQUALS((int)GTEXIT_WAITING_FOR_TIME_STEP, (int)reasons[0]); + TS_ASSERT(traces[0].empty()); + + // Now supply a window that covers t = 15. The seed must start. + VecPt3d west = {{-1, 0, 0}}; + tracer->AddGridScalarsAtTime(west, DataLocationEnum::LOC_CELLS, activity, + DataLocationEnum::LOC_CELLS, 20.0); + tracer->ContinueTraces(); + tracer->GetTraceResults(traces, times, reasons); + + TS_ASSERT(reasons[0] != GTEXIT_EXTRACTION_FAILED); + TS_ASSERT(traces[0].size() >= 2); + TS_ASSERT_DELTA(15.0, times[0].front(), 1e-9); // started at its own release time + TS_ASSERT_DELTA(20.0, traces[0].front().x, 1e-9); + TS_ASSERT_DELTA(10.0, traces[0].front().y, 1e-9); +} // XmGridTraceUnitTests::testSeedReleasedAfterWindowWaitsThenTraces +//------------------------------------------------------------------------------ +/// \brief Two time steps at different data locations must not share a triangulation. +/// +/// Sharing does not copy the triangulation, it shares the object, and the second step's +/// SetGrid*Scalars rebuilds that shared object in place. The shape of the rebuild depends on +/// the data location -- LOC_CELLS adds a centroid point per cell, LOC_POINTS adds none -- so +/// sharing across a location change rebuilt the triangulation the first step was still +/// pointing at, leaving its four-entry scalar array indexed by a centroid index of 4. +/// +/// Both steps here carry the *same* uniform eastward field, written once as point scalars and +/// once as cell scalars, so the interpolated field is identical at every time and the path +/// must be a straight line east. A corrupted first-step lookup cannot produce that. +//------------------------------------------------------------------------------ +void XmGridTraceUnitTests::testDataLocationChangeIsNotShared() +{ + VecPt3d points = {{0, 0, 0}, {40, 0, 0}, {40, 40, 0}, {0, 40, 0}}; + VecInt cells = {XMU_QUAD, 4, 0, 1, 2, 3}; + std::shared_ptr ugrid = XmUGrid::New(points, cells); + + // Activity is cell-based and identical across both steps, so the data location is the only + // term that differs -- which is exactly the term the sharing test used to ignore. + DynBitset activity; + activity.push_back(true); + + BSHP tracer = XmGridTrace::New(ugrid); + tracer->SetVectorMultiplier(1); + tracer->SetMaxTracingTime(8); + tracer->SetMaxTracingDistance(1000); + tracer->SetMinDeltaTime(.01); + tracer->SetMaxChangeDistance(.5); + tracer->SetMaxChangeVelocity(-1); + tracer->SetMaxChangeDirectionInRadians(XM_PI); + + VecPt3d eastAtPoints = {{1, 0, 0}, {1, 0, 0}, {1, 0, 0}, {1, 0, 0}}; + VecPt3d eastAtCells = {{1, 0, 0}}; + tracer->AddGridScalarsAtTime(eastAtPoints, DataLocationEnum::LOC_POINTS, activity, + DataLocationEnum::LOC_CELLS, 0.0); + tracer->AddGridScalarsAtTime(eastAtCells, DataLocationEnum::LOC_CELLS, activity, + DataLocationEnum::LOC_CELLS, 10.0); + + VecPt3d outTrace; + VecDbl outTimes; + tracer->TracePoint({5, 20, 0}, 0, outTrace, outTimes); + + TS_ASSERT(outTrace.size() >= 2); + for (size_t i = 0; i < outTrace.size(); ++i) + { + TS_ASSERT_DELTA(20.0, outTrace[i].y, 1e-9); // pure +x field: y never moves + if (i > 0) + TS_ASSERT(outTrace[i].x > outTrace[i - 1].x); + } + // 8 time units at unit speed from x = 5. + TS_ASSERT_DELTA(13.0, outTrace.back().x, 1e-6); +} // XmGridTraceUnitTests::testDataLocationChangeIsNotShared +//------------------------------------------------------------------------------ /// \brief Verifies the boundary-exit extractor is built once per tracer, not once per exit. /// /// The extractor's constructor triangulates the whole grid and its first SetPolyline indexes diff --git a/xmsgridtrace/gridtrace/XmGridTrace.h b/xmsgridtrace/gridtrace/XmGridTrace.h index 999f21b..f7dec23 100644 --- a/xmsgridtrace/gridtrace/XmGridTrace.h +++ b/xmsgridtrace/gridtrace/XmGridTrace.h @@ -160,6 +160,10 @@ class XmGridTrace /// One batch is in flight per tracer, because the time step window it runs against is /// itself state on the tracer. Starting a batch discards any previous one. /// + /// Release times may be staggered, including past the loaded window: a seed whose time is + /// later than the second loaded step simply waits, with GTEXIT_WAITING_FOR_TIME_STEP, and + /// starts once a window covering it is supplied. + /// /// \param[in] a_pts The starting point of each trace /// \param[in] a_ptTimes The starting time of each trace; must be one per point, or the /// batch is refused entirely diff --git a/xmsgridtrace/gridtrace/XmGridTrace.t.h b/xmsgridtrace/gridtrace/XmGridTrace.t.h index 2963876..c943ad3 100644 --- a/xmsgridtrace/gridtrace/XmGridTrace.t.h +++ b/xmsgridtrace/gridtrace/XmGridTrace.t.h @@ -41,6 +41,9 @@ class XmGridTraceUnitTests : public CxxTest::TestSuite void testTimeVaryingFieldChangesPath(); void testBatchMatchesSerialTracePoint(); void testTracesContinueAcrossTimeSteps(); + void testRedundantContinueDoesNotStallTrace(); + void testSeedReleasedAfterWindowWaitsThenTraces(); + void testDataLocationChangeIsNotShared(); void testBoundaryExtractorIsCached(); void testTraceBenchmark();