limit adaptive refinement only to gauss-kronrod intervals

tkwant/manybody.py

@@ -1704,10 +1704,11 @@ class State:
             Must have the calling signature of `kwant.operator`.
             Default: ``error_op`` from initialization.
         intervals : sequence of `tkwant.manybody.Interval`, optional
-            Apply the refinement process only to the intervals
+            Apply the refinement process only to the (Gauss-Kronrod) intervals
             given in the sequence. Note that in this case, all intervals
             must be present in the solver already. By default,
-            the refinement is done on all intervals stored in the solver.
+            the refinement is done on all Gauss-Kronrod intervals
+            stored in the solver.
 
         Returns
         -------
@@ -1811,7 +1812,7 @@ class State:
 
         # refine only intervals with Gauss-Kronrod quadrature
         intervals[:] = [interval for interval in intervals
-                        if interval.quadrature=='kronrod']
+                        if interval.quadrature == 'kronrod']
 
         results, errors = observable_with_error(intervals)
 
@@ -2013,7 +2014,7 @@ class State:
         .. [2] Gonnet, P., A Review of Error Estimation in Adaptive Quadrature,
                ACM Computing Surveys, Vol. 44, No. 4, Article 22, (2012).
         """
-        if not interval.quadrature=='kronrod':
+        if not interval.quadrature == 'kronrod':
             raise ValueError('quadpack error estimate works only on '
                              'Gauss-Kronrod quadrature intervals')
 

tkwant/tests/test_manybody.py

@@ -1633,6 +1633,36 @@ def test_manybody_state_initial_refine():
     assert len(intervals) == 13
 
 
+@pytest.mark.integtest  # marker for integration tests
+def test_manybody_state_refine_non_kronrod():
+
+    # test that refinement acts only on the gauss-kronrod intervals
+    # and that additional other intervals (here gauss-legendre in the example)
+    # are not touched by refine_intervals()
+
+    syst = make_system().finalized()
+
+    tmax = 500
+    params = {'v': 0.0001, 't0': 0, 'A': 0.00157, 'sigma': 24}
+
+    occupation = manybody.lead_occupation(chemical_potential=3)
+
+    state = manybody.State(syst, tmax, occupation, params=params, refine=False)
+    intervals = state.get_intervals()
+    assert len(intervals) == 2
+
+    for interval in intervals:
+        interval.quadrature = 'gausslegendre'
+        state._add_interval(interval)
+
+    intervals = state.get_intervals()
+    assert len(intervals) == 4  # 2 + 2 intervals
+
+    state.refine_intervals()
+    intervals = state.get_intervals()
+    assert len(intervals) == 15  # 13 + 2 intervals
+
+
 @pytest.mark.integtest
 def test_state_estimate_error():
 
@@ -1742,6 +1772,26 @@ def test_state_estimate_error():
     assert errors.shape == (len(intervals),) + shape_current
 
 
+@pytest.mark.integtest  # marker for integration tests
+def test_manybody_state_estimate_error_non_kronrod():
+
+    syst = make_system().finalized()
+
+    tmax = 500
+    params = {'v': 0.0001, 't0': 0, 'A': 0.00157, 'sigma': 24}
+
+    occupation = manybody.lead_occupation(chemical_potential=3)
+
+    state = manybody.State(syst, tmax, occupation, params=params, refine=False)
+    intervals = state.get_intervals()
+    gauss_interval = intervals[0]
+    gauss_interval.quadrature = 'gausslegendre'
+
+    with pytest.raises(ValueError) as exc:
+        state.estimate_error(intervals=gauss_interval)
+    assert str(exc.value) in "quadpack error estimate works only on Gauss-Kronrod quadrature intervals"
+
+
 @pytest.mark.integtest  # marker for integration tests
 def test_time_name_and_start_default_change_manybody_state():