Try using VIs for bounds constraints

Author: danshapero

notebooks/rainier-monolithic.ipynb

@@ -247,7 +247,11 @@
     "\n",
     "ρ_I = Constant(ice_density)\n",
     "g = Constant(gravity)\n",
-    "n = Constant(glen_flow_law)\n",
+    "\n",
+    "\n",
+    "#n = Constant(glen_flow_law)\n",
+    "n = Constant(1.0)\n",
+    "\n",
     "\n",
     "P = ρ_I * g * h\n",
     "S_n = inner(grad(s), grad(s))**((n - 1) / 2)\n",
@@ -458,10 +462,13 @@
    "metadata": {},
    "outputs": [],
    "source": [
+    "\"\"\"\n",
     "num_continuation_steps = 5\n",
     "for exponent in np.linspace(1.0, 3.0, num_continuation_steps):\n",
     "    n.assign(exponent)\n",
-    "    momentum_solver.solve()"
+    "    momentum_solver.solve()\n",
+    "\"\"\"\n",
+    "momentum_solver.solve()"
    ]
   },
   {
@@ -580,7 +587,30 @@
     "tableau = irksome.BackwardEuler()\n",
     "t = Constant(0.0)\n",
     "dt = Constant(1.0 / 12)\n",
-    "solver = irksome.TimeStepper(F, tableau, t, dt, w, **sparams)"
+    "\n",
+    "lower = firedrake.Function(W)\n",
+    "upper = firedrake.Function(W)\n",
+    "lower.assign(-np.inf)\n",
+    "upper.assign(+np.inf)\n",
+    "lower.subfunctions[3].assign(0.0)\n",
+    "bounds = (\"stage\", lower, upper)\n",
+    "\n",
+    "bparams = {\n",
+    "    \"solver_parameters\": {\n",
+    "        \"snes_vi_monitor\": None,\n",
+    "        \"snes_type\": \"vinewtonrsls\",\n",
+    "        \"snes_max_it\": 200,\n",
+    "        #\"snes_linesearch_type\": \"nleqerr\",\n",
+    "        \"ksp_type\": \"gmres\",\n",
+    "        \"pc_type\": \"lu\",\n",
+    "        \"pc_factor_mat_solver_type\": \"mumps\",\n",
+    "    },\n",
+    "    \"stage_type\": \"value\",\n",
+    "    \"basis_type\": \"Bernstein\",\n",
+    "    #\"bounds\": bounds,\n",
+    "}\n",
+    "\n",
+    "solver = irksome.TimeStepper(F, tableau, t, dt, w, **bparams)"
    ]
   },
   {
@@ -603,10 +633,11 @@
     "\n",
     "final_time = 500.0\n",
     "num_steps = int(final_time / float(dt))\n",
-    "for step in trange(num_steps):\n",
-    "    solver.advance()\n",
+    "progress_bar = trange(num_steps)\n",
+    "for step in progress_bar:\n",
+    "    solver.advance(bounds=bounds)\n",
     "    h = w.subfunctions[3]\n",
-    "    w.subfunctions[3].interpolate(max_value(0, h))\n",
+    "    progress_bar.set_description(str(h.dat.data_ro.min()))\n",
     "    s.interpolate(b + h)\n",
     "    a.interpolate(smb(s))\n",
     "\n",