diff --git a/_package/tests/XmGridTrace_pyt.py b/_package/tests/XmGridTrace_pyt.py index 9a17d48..01ba8a4 100644 --- a/_package/tests/XmGridTrace_pyt.py +++ b/_package/tests/XmGridTrace_pyt.py @@ -1,15 +1,29 @@ """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 -from xms.gridtrace import GridTrace +# 4. Local modules +from xms.gridtrace import exit_reason_enum, GridTrace 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 +84,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 +100,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 +204,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 +258,90 @@ 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 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.""" @@ -278,7 +375,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,10 +386,11 @@ 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]) + 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 @@ -300,7 +398,11 @@ 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]) + # 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.""" @@ -312,7 +414,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 +466,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 +492,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 +544,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 +557,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 +586,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 +613,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 +672,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]) 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..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 @@ -150,3 +158,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 eb2ce0c..9bf761e 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 @@ -45,8 +46,109 @@ namespace { /// XMS Namespace +#ifdef CXX_TEST +/// \brief Count of point-location searches since it was last zeroed. +/// Test-build-only instrumentation for testTraceBenchmark. A trace's cost is dominated by +/// 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 +/// 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_SEARCH(a_n) ((void)0) +/// \brief No-op outside test builds. +#define XMGT_COUNT_BOUNDARY_EXTRACTOR_BUILD() ((void)0) +#endif + //----- 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 +/// \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 + +//------------------------------------------------------------------------------ +/// \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. +/// +/// 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 = 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 + 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 @@ -78,7 +180,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; @@ -87,9 +189,18 @@ class XmGridTraceImpl : public XmGridTrace VecPt3d& a_outTrace, VecDbl& a_outTimes) final; - std::string GetExitMessage() 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; + + XmGridTraceExitEnum GetExitReason() const 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; @@ -113,9 +224,39 @@ 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 - double m_distTraveled=0; ///< distance traveled in the last TracePoint operation - - std::string m_exitMessage; ///< exit message for the last TracePoint operation + xms::DynBitset m_activity2; ///< activity of the second time step, to compare with the next + /// 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 + /// 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; + /// 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; + /// 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; + + /// 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) @@ -254,9 +395,18 @@ 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 //------------------------------------------------------------------------------ -std::string XmGridTraceImpl::GetExitMessage() +const std::string& XmGridTraceImpl::GetExitMessage() const { return m_exitMessage; } // XmGridTraceImpl::GetExitMessage @@ -272,83 +422,155 @@ 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) { - 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 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) - { 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; + m_scalarLoc2 = a_scalarLoc; + m_activityLoc2 = a_activityLoc; } //------------------------------------------------------------------------------ -/// \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; + // 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; + double vx1 = 0, vy1 = 0, mag1 = 0; bool bContinue = true; - Pt3d vtkVec; // Rename this variable - Pt3d vtkPt; + Pt3d vtkVec; Pt3d vector; - - 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)) + 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_exitReason = a_reason; + m_exitMessage = XmGridTraceExitReasonToString(a_reason); + }; + + if (!a_state.m_started) { - m_exitMessage = "Point does not start inside an active cell."; - return; - } + outTrace.clear(); + outTimes.clear(); + 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; + } + 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) { @@ -359,51 +581,66 @@ 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; + 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; } - // 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}; - // 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."); + stopWith(GTEXIT_LEFT_GRID); return; } double segDist = Mdist(pt0.x, pt0.y, pt1.x, pt1.y); @@ -411,10 +648,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; } } @@ -425,9 +663,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; @@ -440,77 +680,153 @@ 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 - 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 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; vy0 = vy1; deltaT *= 1.2; mag0 = mag1; - a_outTimes.push_back(a_ptTime + elapsedTime); + if (moved) + { + 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 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 +//------------------------------------------------------------------------------ +void XmGridTraceImpl::StartTraces(const VecPt3d& a_pts, const VecDbl& a_ptTimes) +{ + m_batch.clear(); + if (a_pts.size() != a_ptTimes.size()) + { + // 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; + } + m_batch.resize(a_pts.size()); + for (size_t i = 0; i < a_pts.size(); ++i) + { + m_batch[i].m_pt = a_pts[i]; + m_batch[i].m_ptTime = a_ptTimes[i]; + } +} // 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 /// \param[in] a_currentTime The time at extraction @@ -520,30 +836,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); - 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); - 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) @@ -551,11 +877,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(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; + 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 = x1 * weight1 + x2 * weight2; + a_data.y = y1 * weight1 + y2 * weight2; return true; } // XmGridTraceImpl::GetVectorAtLocationAndTime } // namespace {} @@ -584,12 +928,50 @@ 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 #include +#include +#include +#include +#include +#include +#include + #include +#include #include using namespace xms; @@ -708,6 +1090,203 @@ 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_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 +}; + +//------------------------------------------------------------------------------ +/// \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_searchCalls = 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_searchCalls = g_searchCalls; + + 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 searchesPerSeed = a_stats.m_searchCalls / seeds; + const double usPerExtract = + 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; + + 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" + << " 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"; + 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 @@ -1084,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 @@ -1296,18 +1879,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 @@ -1336,11 +1921,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 @@ -1446,50 +2036,648 @@ 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 A single-window batch returns exactly what serial TracePoint calls return. +/// +/// 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 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::testBatchMatchesSerialTracePoint() +{ + 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()); + 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); + batchTracer->StartTraces(seeds, seedTimes); + batchTracer->ContinueTraces(); + std::vector batchTraces; + std::vector batchTimes; + 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(), 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); + // 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()); + } + + // The seed outside the grid produces no polyline -- callers cannot assume one per seed. + TS_ASSERT(batchTraces[2].empty()); + 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 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 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 +/// 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 +/// 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, 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 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 +/// 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_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."; + 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.h b/xmsgridtrace/gridtrace/XmGridTrace.h index 5d4c7f5..f7dec23 100644 --- a/xmsgridtrace/gridtrace/XmGridTrace.h +++ b/xmsgridtrace/gridtrace/XmGridTrace.h @@ -33,6 +33,25 @@ 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 +124,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,9 +138,70 @@ class XmGridTrace VecPt3d& a_outTrace, VecDbl& a_outTimes) = 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; + /// \brief 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 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: + /// + /// \code + /// tracer->StartTraces(seeds, seedTimes); + /// while (tracer->ContinueTraces() > 0 && series.HasNext()) + /// tracer->AddGridScalarsAtTime(series.Next(), ...); + /// tracer->GetTraceResults(traces, times, reasons); + /// \endcode + /// + /// 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. + /// + /// 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 + 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 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 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; private: XM_DISALLOW_COPY_AND_ASSIGN(XmGridTrace) @@ -132,4 +212,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 7f6dc78..c943ad3 100644 --- a/xmsgridtrace/gridtrace/XmGridTrace.t.h +++ b/xmsgridtrace/gridtrace/XmGridTrace.t.h @@ -38,6 +38,14 @@ class XmGridTraceUnitTests : public CxxTest::TestSuite void testInactiveCell(); void testStartInactiveCell(); void testTutorial(); + void testTimeVaryingFieldChangesPath(); + void testBatchMatchesSerialTracePoint(); + void testTracesContinueAcrossTimeSteps(); + void testRedundantContinueDoesNotStallTrace(); + void testSeedReleasedAfterWindowWaitsThenTraces(); + void testDataLocationChangeIsNotShared(); + void testBoundaryExtractorIsCached(); + void testTraceBenchmark(); }; // XmGridTraceUnitTests 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",