Skip to content

Commit 5fa57d9

Browse files
danshaperodanshapero
committed
Get rid of old numerical scheme
1 parent 699ecc5 commit 5fa57d9

File tree

2 files changed

+105
-55
lines changed

2 files changed

+105
-55
lines changed

‎notebooks/rainier.ipynb‎

Lines changed: 105 additions & 35 deletions
Original file line numberDiff line numberDiff line change
@@ -13,7 +13,9 @@
1313
"import matplotlib.pyplot as plt\n",
1414
"import firedrake\n",
1515
"from firedrake import inner, grad, dx, ds, exp, min_value, max_value, Constant, derivative\n",
16-
"import irksome"
16+
"import irksome\n",
17+
"import icepack\n",
18+
"from icepack2 import model"
1719
]
1820
},
1921
{
@@ -86,7 +88,7 @@
8688
"B = Constant(4e3)\n",
8789
"r_b = Constant(150e3 / (2 * π))\n",
8890
"expr = B * exp(-inner(x, x) / r_b**2)\n",
89-
"b = firedrake.interpolate(expr, S)"
91+
"b = firedrake.Function(S).interpolate(expr)"
9092
]
9193
},
9294
{
@@ -151,7 +153,7 @@
151153
"r_h = Constant(5e3)\n",
152154
"H = Constant(100.0)\n",
153155
"expr = H * firedrake.max_value(0, 1 - inner(x, x) / r_h**2)\n",
154-
"h = firedrake.interpolate(expr, Q)\n",
156+
"h = firedrake.Function(Q).interpolate(expr)\n",
155157
"h_0 = h.copy(deepcopy=True)"
156158
]
157159
},
@@ -183,8 +185,8 @@
183185
"metadata": {},
184186
"outputs": [],
185187
"source": [
186-
"s = firedrake.interpolate(b + h, Q)\n",
187-
"a = firedrake.interpolate(smb(s), Q)"
188+
"s = firedrake.Function(Q).interpolate(b + h)\n",
189+
"a = firedrake.Function(Q).interpolate(smb(s))"
188190
]
189191
},
190192
{
@@ -215,11 +217,8 @@
215217
"metadata": {},
216218
"outputs": [],
217219
"source": [
218-
"from icepack import rate_factor\n",
219220
"from icepack2.constants import gravity as g, ice_density as ρ_I, glen_flow_law as n\n",
220-
"\n",
221-
"temp = Constant(273.15)\n",
222-
"A = rate_factor(temp)"
221+
"A = icepack.rate_factor(Constant(273.15))"
223222
]
224223
},
225224
{
@@ -300,7 +299,6 @@
300299
"metadata": {},
301300
"outputs": [],
302301
"source": [
303-
"from icepack2 import model, solvers\n",
304302
"fns = [\n",
305303
" model.viscous_power,\n",
306304
" model.friction_power,\n",
@@ -390,13 +388,42 @@
390388
"F = derivative(L, z)"
391389
]
392390
},
391+
{
392+
"cell_type": "markdown",
393+
"id": "9b8cdb6d-1b52-4e65-abd6-968b5eaf545a",
394+
"metadata": {},
395+
"source": [
396+
"For the solution procedure, we'll want to use a linearization of the problem where the thickness is clamped from below."
397+
]
398+
},
399+
{
400+
"cell_type": "code",
401+
"execution_count": null,
402+
"id": "46a2aff8-bee5-4d9d-8f16-b4424484bf14",
403+
"metadata": {},
404+
"outputs": [],
405+
"source": [
406+
"h_min = Constant(1e-3)\n",
407+
"\n",
408+
"rfields = {\n",
409+
" \"velocity\": u,\n",
410+
" \"membrane_stress\": M,\n",
411+
" \"basal_stress\": τ,\n",
412+
" \"thickness\": firedrake.max_value(h_min, h),\n",
413+
" \"surface\": s,\n",
414+
"}\n",
415+
"\n",
416+
"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+
},
393421
{
394422
"cell_type": "markdown",
395423
"id": "468d2245-b155-4ead-8ca1-c37c474d128f",
396424
"metadata": {},
397425
"source": [
398-
"Now we'll create a regularized second derivative of the Lagrangian.\n",
399-
"Here we assume that all the exponents are equal to 1 and the thickness is bounded below."
426+
"We'll also want to add a term to the linearization assuming a linear viscous flow law."
400427
]
401428
},
402429
{
@@ -406,28 +433,36 @@
406433
"metadata": {},
407434
"outputs": [],
408435
"source": [
409-
"h_min = Constant(2.0)\n",
410-
"λ = Constant(1e-3)\n",
411436
"v_c = firedrake.replace(u_c, {h: firedrake.max_value(h, h_min)})\n",
412-
"regularized_rheology = {\n",
437+
"linear_rheology = {\n",
413438
" \"flow_law_exponent\": 1,\n",
414-
" \"flow_law_coefficient\": λ * ε_c / τ_c,\n",
439+
" \"flow_law_coefficient\": ε_c / τ_c,\n",
415440
" \"sliding_exponent\": 1,\n",
416-
" \"sliding_coefficient\": λ * v_c / τ_c,\n",
441+
" \"sliding_coefficient\": v_c / τ_c,\n",
417442
"}\n",
418443
"\n",
419-
"regularized_fields = {\n",
420-
" \"velocity\": u,\n",
421-
" \"membrane_stress\": M,\n",
422-
" \"basal_stress\": τ,\n",
423-
" \"thickness\": firedrake.max_value(h, h_min),\n",
424-
"}\n",
425-
"\n",
426-
"K = sum(\n",
427-
" fn(**regularized_fields, **regularized_rheology)\n",
428-
" for fn in [model.viscous_power, model.friction_power]\n",
429-
")\n",
430-
"J = derivative(derivative(L + K, z), z)"
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"
431466
]
432467
},
433468
{
@@ -456,7 +491,6 @@
456491
"id": "35f45cb5-f7b2-40d7-aced-514de7a1f8dc",
457492
"metadata": {},
458493
"source": [
459-
"The function `try_regularized_solve` will automatically choose a sensible value of how much to regularize $dF$ in order to make things go.\n",
460494
"Now we can see our initial value of the velocity.\n",
461495
"Note how it's a different magnitude but similar shape to the SIA velocity."
462496
]
@@ -468,8 +502,25 @@
468502
"metadata": {},
469503
"outputs": [],
470504
"source": [
471-
"momentum_solver = firedrake.NonlinearVariationalSolver(momentum_problem)\n",
472-
"solvers.try_regularized_solve(momentum_solver, z, λ, verbose=True)"
505+
"sparams = {\n",
506+
" \"snes_type\": \"newtonls\",\n",
507+
" \"snes_max_it\": 200,\n",
508+
" \"snes_linesearch_type\": \"nleqerr\",\n",
509+
" \"ksp_type\": \"gmres\",\n",
510+
" \"pc_type\": \"lu\",\n",
511+
" \"pc_factor_mat_solver_type\": \"umfpack\",\n",
512+
"}\n",
513+
"momentum_solver = firedrake.NonlinearVariationalSolver(momentum_problem, solver_parameters=sparams)"
514+
]
515+
},
516+
{
517+
"cell_type": "code",
518+
"execution_count": null,
519+
"id": "e618ec87-b03f-44a8-8f0d-84b7ba2f58d9",
520+
"metadata": {},
521+
"outputs": [],
522+
"source": [
523+
"momentum_solver.solve()"
473524
]
474525
},
475526
{
@@ -568,7 +619,7 @@
568619
" s.interpolate(b + h)\n",
569620
" a.interpolate(smb(s))\n",
570621
"\n",
571-
" solvers.try_regularized_solve(momentum_solver, z, λ)\n",
622+
" momentum_solver.solve()\n",
572623
"\n",
573624
" hs.append(h.copy(deepcopy=True))\n",
574625
" us.append(z.subfunctions[0].copy(deepcopy=True))"
@@ -594,7 +645,26 @@
594645
"axes.set_aspect(\"equal\")\n",
595646
"axes.set_xlim((0, 10e3))\n",
596647
"axes.set_ylim((0, 10e3))\n",
597-
"colors = firedrake.tripcolor(h, num_sample_points=1, shading=\"flat\", axes=axes)\n",
648+
"colors = firedrake.tripcolor(hs[-1], num_sample_points=1, shading=\"flat\", axes=axes)\n",
649+
"firedrake.triplot(mesh, axes=axes)\n",
650+
"firedrake.tricontour(a, levels=[0, 1], colors=\"tab:orange\", axes=axes)\n",
651+
"fig.colorbar(colors);"
652+
]
653+
},
654+
{
655+
"cell_type": "code",
656+
"execution_count": null,
657+
"id": "f01a49bf-f4c5-41b4-b5fc-3d76deaa55fe",
658+
"metadata": {},
659+
"outputs": [],
660+
"source": [
661+
"u, M, τ = z.subfunctions\n",
662+
"\n",
663+
"fig, axes = plt.subplots()\n",
664+
"axes.set_aspect(\"equal\")\n",
665+
"axes.set_xlim((0, 10e3))\n",
666+
"axes.set_ylim((0, 10e3))\n",
667+
"colors = firedrake.tripcolor(u, num_sample_points=1, shading=\"flat\", axes=axes)\n",
598668
"firedrake.triplot(mesh, axes=axes)\n",
599669
"firedrake.tricontour(a, levels=[0, 1], colors=\"tab:orange\", axes=axes)\n",
600670
"fig.colorbar(colors);"
@@ -653,7 +723,7 @@
653723
"name": "python",
654724
"nbconvert_exporter": "python",
655725
"pygments_lexer": "ipython3",
656-
"version": "3.11.5"
726+
"version": "3.11.9"
657727
}
658728
},
659729
"nbformat": 4,

‎src/icepack2/solvers.py‎

Lines changed: 0 additions & 20 deletions
Original file line numberDiff line numberDiff line change
@@ -3,30 +3,10 @@
33
from firedrake import inner, dx, derivative, action
44

55
__all__ = [
6-
"try_regularized_solve",
76
"ConstrainedOptimizationProblem",
87
"NewtonSolver",
98
]
109

11-
def try_regularized_solve(solver, state, parameter, verbose=False):
12-
state_initial = state.copy(deepcopy=True)
13-
converged = False
14-
while not converged:
15-
try:
16-
if verbose:
17-
print(float(parameter))
18-
solver.solve()
19-
converged = True
20-
except firedrake.ConvergenceError as error:
21-
message = str(error)
22-
if "DIVERGED_MAX_IT" in message:
23-
state_initial.assign(state)
24-
parameter.assign(parameter / 2)
25-
elif "DIVERGED_DTOL" in message:
26-
state.assign(state_initial)
27-
parameter.assign(3 * parameter / 2)
28-
29-
3010
class ConstrainedOptimizationProblem:
3111
def __init__(self, L, z, H=None, bcs=None, form_compiler_parameters=None):
3212
self.objective = L

0 commit comments

Comments
 (0)