Skip to content

Commit b92b347

Browse files
danshaperodanshapero
committed
Use continuation in Rainier notebook
1 parent 1398311 commit b92b347

File tree

1 file changed

+105
-57
lines changed

1 file changed

+105
-57
lines changed

‎notebooks/rainier.ipynb‎

Lines changed: 105 additions & 57 deletions
Original file line numberDiff line numberDiff line change
@@ -12,7 +12,7 @@
1212
"from numpy import pi as π\n",
1313
"import matplotlib.pyplot as plt\n",
1414
"import firedrake\n",
15-
"from firedrake import inner, grad, dx, ds, exp, min_value, max_value, Constant, derivative\n",
15+
"from firedrake import inner, grad, dx, exp, min_value, max_value, Constant, derivative\n",
1616
"import irksome\n",
1717
"import icepack\n",
1818
"from icepack2 import model"
@@ -152,7 +152,7 @@
152152
"source": [
153153
"r_h = Constant(5e3)\n",
154154
"H = Constant(100.0)\n",
155-
"expr = H * firedrake.max_value(0, 1 - inner(x, x) / r_h**2)\n",
155+
"expr = H * max_value(0, 1 - inner(x, x) / r_h**2)\n",
156156
"h = firedrake.Function(Q).interpolate(expr)\n",
157157
"h_0 = h.copy(deepcopy=True)"
158158
]
@@ -217,7 +217,7 @@
217217
"metadata": {},
218218
"outputs": [],
219219
"source": [
220-
"from icepack2.constants import gravity as g, ice_density as ρ_I, glen_flow_law as n\n",
220+
"from icepack2.constants import gravity, ice_density, glen_flow_law\n",
221221
"A = icepack.rate_factor(Constant(273.15))"
222222
]
223223
},
@@ -241,6 +241,10 @@
241241
"u = firedrake.Function(V)\n",
242242
"v = firedrake.TestFunction(V)\n",
243243
"\n",
244+
"ρ_I = Constant(ice_density)\n",
245+
"g = Constant(gravity)\n",
246+
"n = Constant(glen_flow_law)\n",
247+
"\n",
244248
"P = ρ_I * g * h\n",
245249
"S_n = inner(grad(s), grad(s))**((n - 1) / 2)\n",
246250
"u_shear = -2 * A * P ** n / (n + 2) * h * S_n * grad(s)\n",
@@ -403,66 +407,19 @@
403407
"metadata": {},
404408
"outputs": [],
405409
"source": [
406-
"h_min = Constant(1e-3)\n",
410+
"h_min = Constant(1e-5)\n",
407411
"\n",
408412
"rfields = {\n",
409413
" \"velocity\": u,\n",
410414
" \"membrane_stress\": M,\n",
411415
" \"basal_stress\": τ,\n",
412-
" \"thickness\": firedrake.max_value(h_min, h),\n",
416+
" \"thickness\": max_value(h_min, h),\n",
413417
" \"surface\": s,\n",
414418
"}\n",
415419
"\n",
416420
"L_r = sum(fn(**rfields, **rheology) for fn in fns)\n",
417-
"F_r = firedrake.derivative(L_r, z)\n",
418-
"J_r = firedrake.derivative(F_r, z)"
419-
]
420-
},
421-
{
422-
"cell_type": "markdown",
423-
"id": "468d2245-b155-4ead-8ca1-c37c474d128f",
424-
"metadata": {},
425-
"source": [
426-
"We'll also want to add a term to the linearization assuming a linear viscous flow law."
427-
]
428-
},
429-
{
430-
"cell_type": "code",
431-
"execution_count": null,
432-
"id": "37dba801-9394-4719-99d1-fdec7522f7a5",
433-
"metadata": {},
434-
"outputs": [],
435-
"source": [
436-
"v_c = firedrake.replace(u_c, {h: firedrake.max_value(h, h_min)})\n",
437-
"linear_rheology = {\n",
438-
" \"flow_law_exponent\": 1,\n",
439-
" \"flow_law_coefficient\": ε_c / τ_c,\n",
440-
" \"sliding_exponent\": 1,\n",
441-
" \"sliding_coefficient\": v_c / τ_c,\n",
442-
"}\n",
443-
"\n",
444-
"L_1 = sum(fn(**rfields, **linear_rheology) for fn in fns)\n",
445-
"F_1 = firedrake.derivative(L_1, z)\n",
446-
"J_1 = firedrake.derivative(F_1, z)"
447-
]
448-
},
449-
{
450-
"cell_type": "markdown",
451-
"id": "39fbac26-12e1-408a-8aab-6c32e321edd5",
452-
"metadata": {},
453-
"source": [
454-
"Finally, the linearization that we'll use for the nonlinear solver is a combination of the two (viscous and Glen)."
455-
]
456-
},
457-
{
458-
"cell_type": "code",
459-
"execution_count": null,
460-
"id": "3ad45f70-18a5-44f5-9d6f-2efd2a3c3e46",
461-
"metadata": {},
462-
"outputs": [],
463-
"source": [
464-
"λ = Constant(1e-3)\n",
465-
"J = J_r + λ * J_1"
421+
"F_r = derivative(L_r, z)\n",
422+
"J_r = derivative(F_r, z)"
466423
]
467424
},
468425
{
@@ -481,9 +438,9 @@
481438
"outputs": [],
482439
"source": [
483440
"degree = 1\n",
484-
"qdegree = max(8, degree ** n)\n",
441+
"qdegree = max(8, degree ** glen_flow_law)\n",
485442
"pparams = {\"form_compiler_parameters\": {\"quadrature_degree\": qdegree}}\n",
486-
"momentum_problem = firedrake.NonlinearVariationalProblem(F, z, J=J, **pparams)"
443+
"momentum_problem = firedrake.NonlinearVariationalProblem(F, z, J=J_r, **pparams)"
487444
]
488445
},
489446
{
@@ -503,6 +460,7 @@
503460
"outputs": [],
504461
"source": [
505462
"sparams = {\n",
463+
" \"snes_monitor\": \":rainier-output.txt\",\n",
506464
" \"snes_type\": \"newtonls\",\n",
507465
" \"snes_max_it\": 200,\n",
508466
" \"snes_linesearch_type\": \"nleqerr\",\n",
@@ -520,7 +478,10 @@
520478
"metadata": {},
521479
"outputs": [],
522480
"source": [
523-
"momentum_solver.solve()"
481+
"num_continuation_steps = 5\n",
482+
"for exponent in np.linspace(1.0, 3.0, num_continuation_steps):\n",
483+
" n.assign(exponent)\n",
484+
" momentum_solver.solve()"
524485
]
525486
},
526487
{
@@ -705,6 +666,93 @@
705666
"from IPython.display import HTML\n",
706667
"HTML(animation.to_html5_video())"
707668
]
669+
},
670+
{
671+
"cell_type": "markdown",
672+
"id": "a1519ecf-2c80-4de1-9f54-33b861435322",
673+
"metadata": {},
674+
"source": [
675+
"Finally, let's check that the solver was working right."
676+
]
677+
},
678+
{
679+
"cell_type": "code",
680+
"execution_count": null,
681+
"id": "6a0ec125-564b-4457-bd47-71eb911e7a8b",
682+
"metadata": {},
683+
"outputs": [],
684+
"source": [
685+
"junk = \"SNES Function norm\"\n",
686+
"with open(\"rainier-output.txt\", \"r\") as log_file:\n",
687+
" raw_data = [\n",
688+
" line.replace(junk, \"\").strip().split()\n",
689+
" for line in log_file.readlines()\n",
690+
" ]\n",
691+
"\n",
692+
"convergence_log = [\n",
693+
" (int(iter_text), float(norm_text)) for iter_text, norm_text in raw_data\n",
694+
"]\n",
695+
"\n",
696+
"reduced_log = []\n",
697+
"current_entry = []\n",
698+
"for count, error in convergence_log:\n",
699+
" if count == 0:\n",
700+
" reduced_log.append(current_entry)\n",
701+
" current_entry = []\n",
702+
"\n",
703+
" current_entry.append(error)\n",
704+
"\n",
705+
"reduced_log.pop(0);"
706+
]
707+
},
708+
{
709+
"cell_type": "markdown",
710+
"id": "a4d62643-ae8f-4a5f-9a22-d45166505a45",
711+
"metadata": {},
712+
"source": [
713+
"First, we want to see that the number of Newton iterations was bounded.\n",
714+
"If it starts reaching above 20 then you're probably in trouble."
715+
]
716+
},
717+
{
718+
"cell_type": "code",
719+
"execution_count": null,
720+
"id": "1768b1cc-a969-4baf-9b09-4d0e83cdafbb",
721+
"metadata": {},
722+
"outputs": [],
723+
"source": [
724+
"fig, ax = plt.subplots()\n",
725+
"ax.set_ylim((0, 10))\n",
726+
"ax.set_xlabel(\"Time step\")\n",
727+
"ax.set_ylabel(\"#Newton steps\")\n",
728+
"ax.plot([len(entry) for entry in reduced_log]);"
729+
]
730+
},
731+
{
732+
"cell_type": "markdown",
733+
"id": "a0448417-72b3-4e5c-abba-1ad4a20f319c",
734+
"metadata": {},
735+
"source": [
736+
"Next, we'll have a look at the initial and final residual values.\n",
737+
"We used a convergence criterion where the iteration stops if the residual is reduced by more than something like $10^4$.\n",
738+
"It would be cause for concern if the final value of the residual were edging upward, but this isn't the case."
739+
]
740+
},
741+
{
742+
"cell_type": "code",
743+
"execution_count": null,
744+
"id": "eedac523-9bd0-4b48-8ba4-81d2622df72c",
745+
"metadata": {},
746+
"outputs": [],
747+
"source": [
748+
"fig, ax = plt.subplots()\n",
749+
"ax.set_yscale(\"log\")\n",
750+
"ax.set_xlabel(\"Time step\")\n",
751+
"ax.set_ylabel(\"Residual\")\n",
752+
"ax.plot([entry[0] for entry in reduced_log], label=\"Initial\")\n",
753+
"ax.plot([entry[-1] for entry in reduced_log], label=\"Final\")\n",
754+
"ax.legend(loc=\"upper right\");"
755+
]
708756
}
709757
],
710758
"metadata": {

0 commit comments

Comments
 (0)