Nonlinear Poisson for glacier channel flow

[loganemann]

Hello. I am attempting to solve a nonlinear Poisson equation for velocity $u(y,z)$ in a cross section of the half-width of a rectangular glacier channel.
$$\nabla\cdot\left[q(u) \nabla u \right] = -\frac{\tau_d}{H}$$
where
$$q(u) = \frac{1}{2}A^{-1/n}\left[ \frac{1}{4} \left(u_z^2 + u_y^2\right)\right]^{\frac{1-n}{2n}}$$
with a no-slip boundary condition at the base ($z=0$) and right sidewall ($y=W$), and a stress-free boundary condition ($\boldsymbol{\sigma}\cdot\mathbf{n} = 0$) at the ice surface ($z=H$) and centerline ($y=0$).
$u = 0, \ \ z=0$
$u= 0, \ \ y=W$
$q(u) u_z = 0, \ \ z=H$
$q(u) u_y = 0, \ \ y=0$

The weak form should be

$$\int_{\Omega}q(u) \nabla u \nabla v dV = \int_{\Omega} \frac{\tau_d}{H}v dV + \int_{\Gamma} v q(u) \nabla u \cdot \mathbf{n} dS$$

I'm relatively new to firedrake, so I suspect that my issue likely stems from some misunderstanding of ufl.

from firedrake import *
length = 1000
width = 1000
mesh = RectangleMesh(40, 40, length, width)

V = FunctionSpace(mesh, "Lagrange", 1)
y, z = SpatialCoordinate(mesh)

# Parameters
rheo = Constant(1.5e-25)
n=Constant(3)              
rho = Constant(917)       
g = Constant(9.81)
H = Constant(1000)
s_x = Constant(-0.01)
tau_d = -rho*g*H*s_x
f = tau_d/H

def q(u):
    return 0.5*rheo**(-1/n)*(0.25*(Dx(u,0))**2 + 0.25*(Dx(u,1))**2 + 1e-15)**((1-n)/(2*n))

u = TrialFunction(V)
v = TestFunction(V)

F = (q(u)*dot(grad(u),grad(v)) + tau_d/H*v)*dx

a = (0.5/rheo*inner(grad(u),grad(v)))*dx
L = (f*v)*dx

bc1 = DirichletBC(V, 0, 2)
bc2 = DirichletBC(V, 0, 3)

U = Function(V)

solve(F == 0, U, bcs=[bc1, bc2])

I'm receiving this error message:

---------------------------------------------------------------------------
ValueError                                Traceback (most recent call last)
Cell In[9], line 5
      1 U = Function(V)
      2 #solve(a == L, uh, bcs=[bc1, bc2])
      3 #solve(a == L, uh)
----> 5 solve(F == 0, U, bcs=[bc1, bc2], solver_parameters={"snes_monitor": None})

File petsc4py[/PETSc/Log.pyx:188](http://localhost:8888/PETSc/Log.pyx#line=187), in petsc4py.PETSc.Log.EventDecorator.decorator.wrapped_func()

File petsc4py[/PETSc/Log.pyx:189](http://localhost:8888/PETSc/Log.pyx#line=188), in petsc4py.PETSc.Log.EventDecorator.decorator.wrapped_func()

File [~/envs/firedrake/src/firedrake/firedrake/adjoint_utils/solving.py:57](http://localhost:8888/lab/tree/~/envs/firedrake/src/firedrake/firedrake/adjoint_utils/solving.py#line=56), in annotate_solve.<locals>.wrapper(*args, **kwargs)
     54     tape.add_block(block)
     56 with stop_annotating():
---> 57     output = solve(*args, **kwargs)
     59 if annotate:
     60     if hasattr(args[1], "create_block_variable"):

File [~/envs/firedrake/src/firedrake/firedrake/solving.py:132](http://localhost:8888/lab/tree/~/envs/firedrake/src/firedrake/firedrake/solving.py#line=131), in solve(*args, **kwargs)
    130 # Call variational problem solver if we get an equation
    131 if isinstance(args[0], ufl.classes.Equation):
--> 132     _solve_varproblem(*args, **kwargs)
    133 else:
    134     # Solve pre-assembled system
    135     return _la_solve(*args, **kwargs)

File [~/envs/firedrake/src/firedrake/firedrake/solving.py:176](http://localhost:8888/lab/tree/~/envs/firedrake/src/firedrake/firedrake/solving.py#line=175), in _solve_varproblem(*args, **kwargs)
    174     raise TypeError("Only '0' support on RHS of nonlinear Equation, not %r" % eq.rhs)
    175 # Create problem
--> 176 problem = vs.NonlinearVariationalProblem(eq.lhs, u, bcs, J, Jp,
    177                                          form_compiler_parameters=form_compiler_parameters,
    178                                          restrict=restrict)
    179 # Create solver and call solve
    180 solver = vs.NonlinearVariationalSolver(problem, solver_parameters=solver_parameters,
    181                                        nullspace=nullspace,
    182                                        transpose_nullspace=nullspace_T,
    183                                        near_nullspace=near_nullspace,
    184                                        options_prefix=options_prefix,
    185                                        appctx=appctx)

File petsc4py[/PETSc/Log.pyx:188](http://localhost:8888/PETSc/Log.pyx#line=187), in petsc4py.PETSc.Log.EventDecorator.decorator.wrapped_func()

File petsc4py[/PETSc/Log.pyx:189](http://localhost:8888/PETSc/Log.pyx#line=188), in petsc4py.PETSc.Log.EventDecorator.decorator.wrapped_func()

File [~/envs/firedrake/src/pyadjoint/pyadjoint/tape.py:110](http://localhost:8888/lab/tree/~/envs/firedrake/src/pyadjoint/pyadjoint/tape.py#line=109), in no_annotations.<locals>.wrapper(*args, **kwargs)
    107 @wraps(function)
    108 def wrapper(*args, **kwargs):
    109     with stop_annotating():
--> 110         return function(*args, **kwargs)

File [~/envs/firedrake/src/firedrake/firedrake/adjoint_utils/variational_solver.py:15](http://localhost:8888/lab/tree/~/envs/firedrake/src/firedrake/firedrake/adjoint_utils/variational_solver.py#line=14), in NonlinearVariationalProblemMixin._ad_annotate_init.<locals>.wrapper(self, *args, **kwargs)
     11 @no_annotations
     12 @wraps(init)
     13 def wrapper(self, *args, **kwargs):
     14     from firedrake import derivative, adjoint, TrialFunction
---> 15     init(self, *args, **kwargs)
     16     self._ad_F = self.F
     17     self._ad_u = self.u_restrict

File [~/envs/firedrake/src/firedrake/firedrake/variational_solver.py:111](http://localhost:8888/lab/tree/~/envs/firedrake/src/firedrake/firedrake/variational_solver.py#line=110), in NonlinearVariationalProblem.__init__(self, F, u, bcs, J, Jp, form_compiler_parameters, is_linear, restrict)
    108 self.Jp_eq_J = Jp is None
    110 # Argument checking
--> 111 check_pde_args(self.F, self.J, self.Jp)
    113 # Store form compiler parameters
    114 self.form_compiler_parameters = form_compiler_parameters

File [~/envs/firedrake/src/firedrake/firedrake/variational_solver.py:25](http://localhost:8888/lab/tree/~/envs/firedrake/src/firedrake/firedrake/variational_solver.py#line=24), in check_pde_args(F, J, Jp)
     23     raise TypeError("Provided residual is a '%s', not a BaseForm or Slate Tensor" % type(F).__name__)
     24 if len(F.arguments()) != 1:
---> 25     raise ValueError("Provided residual is not a linear form")
     26 if not isinstance(J, (ufl.BaseForm, slate.slate.TensorBase)):
     27     raise TypeError("Provided Jacobian is a '%s', not a BaseForm or Slate Tensor" % type(J).__name__)

ValueError: Provided residual is not a linear form

[ksagiyam]

Linear and nonlinear problems are specified differently. If you want to use F as the residual for a nonlinear problem to be used as solve(F == 0, u,...), u (the solution Function) must be used to define F (instead of TrialFunction(V) as you currently do).

Please see https://www.firedrakeproject.org/variational-problems.html#expressing-a-variational-problem for a full explanation.

[loganemann]

Looks like defining u = Function(V) resolved the linear form problem. Now i suspect there may be an issue in how I'm defining the residual F = (q(u)*dot(grad(u),grad(v)) + tau_d/H*v)*dx

---------------------------------------------------------------------------
ValueError                                Traceback (most recent call last)
Cell In[15], line 5
      1 U = Function(V)
      2 #solve(a == L, uh, bcs=[bc1, bc2])
      3 #solve(a == L, uh)
----> 5 solve(F == 0, U, bcs=[bc1, bc2])
      6 # solve(F == 0, U)

File petsc4py[/PETSc/Log.pyx:188](http://localhost:8888/PETSc/Log.pyx#line=187), in petsc4py.PETSc.Log.EventDecorator.decorator.wrapped_func()

File petsc4py[/PETSc/Log.pyx:189](http://localhost:8888/PETSc/Log.pyx#line=188), in petsc4py.PETSc.Log.EventDecorator.decorator.wrapped_func()

File [~/envs/firedrake/src/firedrake/firedrake/adjoint_utils/solving.py:57](http://localhost:8888/lab/tree/~/envs/firedrake/src/firedrake/firedrake/adjoint_utils/solving.py#line=56), in annotate_solve.<locals>.wrapper(*args, **kwargs)
     54     tape.add_block(block)
     56 with stop_annotating():
---> 57     output = solve(*args, **kwargs)
     59 if annotate:
     60     if hasattr(args[1], "create_block_variable"):

File [~/envs/firedrake/src/firedrake/firedrake/solving.py:132](http://localhost:8888/lab/tree/~/envs/firedrake/src/firedrake/firedrake/solving.py#line=131), in solve(*args, **kwargs)
    130 # Call variational problem solver if we get an equation
    131 if isinstance(args[0], ufl.classes.Equation):
--> 132     _solve_varproblem(*args, **kwargs)
    133 else:
    134     # Solve pre-assembled system
    135     return _la_solve(*args, **kwargs)

File [~/envs/firedrake/src/firedrake/firedrake/solving.py:176](http://localhost:8888/lab/tree/~/envs/firedrake/src/firedrake/firedrake/solving.py#line=175), in _solve_varproblem(*args, **kwargs)
    174     raise TypeError("Only '0' support on RHS of nonlinear Equation, not %r" % eq.rhs)
    175 # Create problem
--> 176 problem = vs.NonlinearVariationalProblem(eq.lhs, u, bcs, J, Jp,
    177                                          form_compiler_parameters=form_compiler_parameters,
    178                                          restrict=restrict)
    179 # Create solver and call solve
    180 solver = vs.NonlinearVariationalSolver(problem, solver_parameters=solver_parameters,
    181                                        nullspace=nullspace,
    182                                        transpose_nullspace=nullspace_T,
    183                                        near_nullspace=near_nullspace,
    184                                        options_prefix=options_prefix,
    185                                        appctx=appctx)

File petsc4py[/PETSc/Log.pyx:188](http://localhost:8888/PETSc/Log.pyx#line=187), in petsc4py.PETSc.Log.EventDecorator.decorator.wrapped_func()

File petsc4py[/PETSc/Log.pyx:189](http://localhost:8888/PETSc/Log.pyx#line=188), in petsc4py.PETSc.Log.EventDecorator.decorator.wrapped_func()

File [~/envs/firedrake/src/pyadjoint/pyadjoint/tape.py:110](http://localhost:8888/lab/tree/~/envs/firedrake/src/pyadjoint/pyadjoint/tape.py#line=109), in no_annotations.<locals>.wrapper(*args, **kwargs)
    107 @wraps(function)
    108 def wrapper(*args, **kwargs):
    109     with stop_annotating():
--> 110         return function(*args, **kwargs)

File [~/envs/firedrake/src/firedrake/firedrake/adjoint_utils/variational_solver.py:26](http://localhost:8888/lab/tree/~/envs/firedrake/src/firedrake/firedrake/adjoint_utils/variational_solver.py#line=25), in NonlinearVariationalProblemMixin._ad_annotate_init.<locals>.wrapper(self, *args, **kwargs)
     20 try:
     21     # Some forms (e.g. SLATE tensors) are not currently
     22     # differentiable.
     23     dFdu = derivative(self.F,
     24                       self.u_restrict,
     25                       TrialFunction(self.u_restrict.function_space()))
---> 26     self._ad_adj_F = adjoint(dFdu)
     27 except (TypeError, NotImplementedError):
     28     self._ad_adj_F = None

File petsc4py[/PETSc/Log.pyx:188](http://localhost:8888/PETSc/Log.pyx#line=187), in petsc4py.PETSc.Log.EventDecorator.decorator.wrapped_func()

File petsc4py[/PETSc/Log.pyx:189](http://localhost:8888/PETSc/Log.pyx#line=188), in petsc4py.PETSc.Log.EventDecorator.decorator.wrapped_func()

File [~/envs/firedrake/src/firedrake/firedrake/ufl_expr.py:346](http://localhost:8888/lab/tree/~/envs/firedrake/src/firedrake/firedrake/ufl_expr.py#line=345), in adjoint(form, reordered_arguments, derivatives_expanded)
    339     v, u = extract_arguments(form)
    340     reordered_arguments = (Argument(u.function_space(),
    341                                     number=v.number(),
    342                                     part=v.part()),
    343                            Argument(v.function_space(),
    344                                     number=u.number(),
    345                                     part=u.part()))
--> 346 return ufl.adjoint(form, reordered_arguments, derivatives_expanded=derivatives_expanded)

File [~/envs/firedrake/src/ufl/ufl/formoperators.py:160](http://localhost:8888/lab/tree/~/envs/firedrake/src/ufl/ufl/formoperators.py#line=159), in adjoint(form, reordered_arguments, derivatives_expanded)
    158         form = expand_derivatives(form)
    159     if isinstance(form, Form):
--> 160         return compute_form_adjoint(form, reordered_arguments)
    161 return Adjoint(form)

File [~/envs/firedrake/src/ufl/ufl/algorithms/formtransformations.py:458](http://localhost:8888/lab/tree/~/envs/firedrake/src/ufl/ufl/algorithms/formtransformations.py#line=457), in compute_form_adjoint(form, reordered_arguments)
    455     raise ValueError("compute_form_adjoint cannot handle parts.")
    457 if len(arguments) != 2:
--> 458     raise ValueError("Expecting bilinear form.")
    460 v, u = arguments
    461 if v.number() >= u.number():

ValueError: Expecting bilinear form.

[loganemann]

Never mind, figured it out. I need to call solve(F == 0, u,...), rather than solve(F == 0, U,...). Thanks!