Merge pull request #6 from gdmcbain/milestones

milestones for natural continuation

pacopy/natural.py

@@ -9,13 +9,14 @@ def natural(
     lambda0,
     callback,
     lambda_stepsize0=1.0e-1,
-    lambda_stepsize_max=1.0e0,
+    lambda_stepsize_max=float("inf"),
     lambda_stepsize_aggressiveness=2,
     max_newton_steps=5,
     newton_tol=1.0e-12,
     max_steps=float("inf"),
     verbose=True,
     use_first_order_predictor=True,
+    milestones=None,
 ):
     """Natural parameter continuation.
 
@@ -42,8 +43,11 @@ def natural(
             to bootstrap the Newton process for the next iteration (order 0). Another
             possibility is to use :math:`u - s J^{-1}(u, \\lambda)
             \\frac{df}{d\\lambda}`, a first-order approximation.
+        milestones (Optional[Iterable[float]]): Don't step over these values.
     """
     lmbda = lambda0
+    if milestones is not None:
+        milestones = iter(milestones)
 
     k = 0
     try:
@@ -63,6 +67,8 @@ def natural(
     k += 1
 
     lambda_stepsize = lambda_stepsize0
+    if milestones is not None:
+        milestone = next(milestones)
 
     while True:
         if k > max_steps:
@@ -77,6 +83,8 @@ def natural(
 
         # Predictor
         lmbda += lambda_stepsize
+        if milestones:
+            lmbda = min(lmbda, milestone)
         if use_first_order_predictor:
             du_dlmbda = problem.jacobian_solver(u, lmbda, -problem.df_dlmbda(u, lmbda))
             u0 = u + du_dlmbda * lambda_stepsize
@@ -100,14 +108,19 @@ def natural(
             lambda_stepsize /= 2
             continue
 
-        lambda_stepsize *= (
-            1
-            + lambda_stepsize_aggressiveness
-            * ((max_newton_steps - newton_steps) / (max_newton_steps - 1)) ** 2
-        )
-        lambda_stepsize = min(lambda_stepsize, lambda_stepsize_max)
-
         callback(k, lmbda, u)
         k += 1
+        if milestones is not None and lmbda == milestone:
+            try:
+                milestone = next(milestones)
+            except StopIteration:
+                break
+        else:
+            lambda_stepsize *= (
+                1
+                + lambda_stepsize_aggressiveness
+                * ((max_newton_steps - newton_steps) / (max_newton_steps - 1)) ** 2
+            )
+            lambda_stepsize = min(lambda_stepsize, lambda_stepsize_max)
 
     return

test/test_bratu1d.py

@@ -80,6 +80,10 @@ def test_bratu(max_steps=10, update_plot=False):
     # line3, = ax2.plot([], [], "-", color="red")
     # line3.set_xdata(numpy.linspace(0.0, 1.0, problem.n))
 
+    milestones = numpy.arange(0.5, 3.2, 0.5)
+    if update_plot:
+        profile_fig, profile_ax = plt.subplots()
+
     def callback(k, lmbda, sol):
         if update_plot:
             line1.set_xdata(numpy.append(line1.get_xdata(), lmbda))
@@ -99,9 +103,20 @@ def callback(k, lmbda, sol):
             fig.canvas.draw()
             fig.canvas.flush_events()
             # plt.savefig('bratu1d.png'.format(k), transparent=True, bbox_inches="tight")
+            if lmbda in milestones:
+                profile_ax.plot(numpy.linspace(0.0, 1.0, problem.n), sol, label=lmbda)
+
         return
 
-    pacopy.natural(problem, u0, lmbda0, callback, max_steps=max_steps)
+    pacopy.natural(
+        problem,
+        u0,
+        lmbda0,
+        callback,
+        max_steps=max_steps,
+        newton_tol=1e-10,
+        milestones=milestones,
+    )
 
     # The condition number of the Jacobian is about 10^4, so we can only expect Newton
     # to converge up to about this factor above machine precision.