Start using icepack2 library in gibbous demo

demos/gibbous-ice-shelf/gibbous.py

@@ -10,11 +10,12 @@
     grad,
     dx,
     ds,
+    derivative,
     NonlinearVariationalProblem,
     NonlinearVariationalSolver,
 )
-from icepack.constants import glen_flow_law
-from dualform import ice_shelf
+from icepack2 import model
+from icepack2.constants import glen_flow_law as n
 
 parser = argparse.ArgumentParser()
 parser.add_argument("--input")
@@ -138,31 +139,33 @@
     "thickness": h,
 }
 
+fns = [model.viscous_power, model.ice_shelf_momentum_balance]
+
+rheology = {
+    "flow_law_exponent": n,
+    "flow_law_coefficient": ε_c / τ_c ** n,
+}
+
+L = sum(fn(**fields, **rheology,) for fn in fns)
+F = derivative(L, z)
+
 h_min = firedrake.Constant(1.0)
 rfields = {
     "velocity": u,
     "membrane_stress": M,
     "thickness": firedrake.max_value(h_min, h),
 }
 
-params = {
-    "viscous_yield_strain": ε_c,
-    "viscous_yield_stress": τ_c,
-    "outflow_ids": (2,),
+λ = firedrake.Constant(1e-1)
+rrheology = {
+    "flow_law_exponent": 1,
+    "flow_law_coefficient": λ * ε_c / τ_c,
 }
 
-fns = [ice_shelf.viscous_power, ice_shelf.boundary, ice_shelf.constraint]
-L_1 = sum(fn(**fields, **params, flow_law_exponent=1) for fn in fns)
-F_1 = firedrake.derivative(L_1, z)
-
-L = sum(fn(**fields, **params, flow_law_exponent=glen_flow_law) for fn in fns)
-F = firedrake.derivative(L, z)
-
-L_r = sum(fn(**rfields, **params, flow_law_exponent=glen_flow_law) for fn in fns)
-F_r = firedrake.derivative(L_r, z)
-J_r = firedrake.derivative(F_r, z)
+K = model.viscous_power(**fields, **rrheology)
+J = derivative(derivative(L + K, z), z)
 
-qdegree = int(glen_flow_law) + 2
+qdegree = n + 2
 bc = firedrake.DirichletBC(Z.sub(0), u0, (1,))
 problem_params = {
     #"form_compiler_parameters": {"quadrature_degree": qdegree},
@@ -177,13 +180,12 @@
         "snes_rtol": 1e-2,
     },
 }
-firedrake.solve(F_1 == 0, z, **problem_params, **solver_params)
+firedrake.solve(F == 0, z, J=J, **problem_params, **solver_params)
 
-# Create a regularized Lagrangian
 
 # Set up the diagnostic problem and solver
 # TODO: possibly use a regularized Jacobian
-diagnostic_problem = NonlinearVariationalProblem(F, z, J=J_r, **problem_params)
+diagnostic_problem = NonlinearVariationalProblem(F, z, J=J, **problem_params)
 diagnostic_solver = NonlinearVariationalSolver(diagnostic_problem, **solver_params)
 diagnostic_solver.solve()