Back to old regularization scheme

demos/gibbous-ice-shelf/gibbous.py

@@ -132,8 +132,7 @@
 else:
     z.sub(0).assign(u0)
 
-# Set up the diagnostic problem and compute an initial guess by solving a
-# Picard linearization of the problem
+# Set up the diagnostic problem
 u, M = firedrake.split(z)
 fields = {
     "velocity": u,
@@ -148,24 +147,17 @@
     "flow_law_coefficient": ε_c / τ_c ** n,
 }
 
-L = sum(fn(**fields, **rheology,) for fn in fns)
+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),
-}
-
-λ = firedrake.Constant(1e-1)
-rrheology = {
+# Make an initial guess by solving a Picard linearization
+linear_rheology = {
     "flow_law_exponent": 1,
-    "flow_law_coefficient": λ * ε_c / τ_c,
+    "flow_law_coefficient": ε_c / τ_c,
 }
 
-K = model.viscous_power(**rfields, **rrheology)
-J = derivative(derivative(L + K, z), z)
+L_1 = sum(fn(**fields, **linear_rheology) for fn in fns)
+F_1 = derivative(L_1, z)
 
 qdegree = n + 2
 bc = firedrake.DirichletBC(Z.sub(0), u0, (1,))
@@ -182,12 +174,23 @@
         "snes_rtol": 1e-2,
     },
 }
-firedrake.solve(F == 0, z, J=J, **problem_params, **solver_params)
+firedrake.solve(F_1 == 0, z, **problem_params, **solver_params)
+
+# Create an approximate linearization
+h_min = firedrake.Constant(1.0)
+rfields = {
+    "velocity": u,
+    "membrane_stress": M,
+    "thickness": firedrake.max_value(h_min, h),
+}
 
+L_r = sum(fn(**rfields, **rheology) for fn in fns)
+F_r = derivative(L_r, z)
+J_r = derivative(F_r, z)
 
 # Set up the diagnostic problem and solver
 # TODO: possibly use a regularized Jacobian
-diagnostic_problem = NonlinearVariationalProblem(F, z, J=J, **problem_params)
+diagnostic_problem = NonlinearVariationalProblem(F, z, J=J_r, **problem_params)
 diagnostic_solver = NonlinearVariationalSolver(diagnostic_problem, **solver_params)
 diagnostic_solver.solve()