# -*- coding=utf-8 -*- # The Monge-Ampère equation # ========================= # # .. rst-class:: emphasis # # This tutorial was contributed by Colin Cotter # __, based on code from Andrew # McRae __ and Lawrence Mitchell # __. # # The Monge-Ampère equation provides the solution to the optimal # transportation problem between two measures. Here, we consider the # case where the target measure is the usual Lebesgue measure, and the # template measure is :math:f(x)d^nx, both defined on the same # domain. Then, in two dimensions, the optimal transportation plan is # given by # # .. math:: # (x,y) \mapsto (x,y) + \nabla u, # # where :math:u satisfies the the Monge-Ampère equation # # .. math:: # # \det\left(I + D^2 u\right) = f, # # where :math:I is the identity matrix, and :math:D^2 is # the Hessian matrix of second derivatives, subject to the boundary # conditions :math:\frac{\partial u}{\partial n}=0. # # Here we follow the approach of :cite:lakkis2013finite, namely # to use the mixed formulation # # .. math:: # \sigma = D^2 u, # # \det(I + \sigma) = f, # # where :math:\sigma is a :math:2\times 2 tensor. # # Written in weak form, our problem is to find :math:(u, \sigma) \in # V\times \Sigma = W such that # # .. math:: # \int_\Omega \tau:(\sigma + D^2u)\,\mathrm{d}x # - \int_{\partial\Omega} \tau_{12}n_2u_x + \tau_{21}n_1u_y\,\mathrm{d}s # &=0, \quad \forall \tau \in \Sigma # # \int_\Omega v\det(I + \sigma)\,\mathrm{d}x = \int_\Omega # fv\,\mathrm{d}x \quad \forall v \in V. # # This is called a nonvariational discretisation since the PDE is not in # a divergence form. Note that we have dropped the boundary terms that # vanish due to the boundary condition. To proceed in the # discretisation, we simply choose :math:V to be a continuous # degree-k finite element space, and :math:\Sigma to be the :math:2 \times 2 # tensor continuous finite element space of the same degree. Since we have # Neumann boundary conditions, this variational problem has a null space # consisting of the constant functions in :math:V. # # For Dirichlet boundary conditions, :cite:awanou2014quadratic proved # that this algorithm converges when :math:k>1. Note that # the Jacobian system arising from Newton is only elliptic when # :math:I + \sigma is positive-definite; it is observed that # positive-definiteness is preserved by Newton iteration and hence we # must be careful to choose an appropriate initial guess. This is one of # the reasons why we have set things up with :math:I + \sigma here # instead of :math:\sigma as is more conventional for these equations, # since then :math:u=0 is an appropriate initial guess. This setup # also makes the application of the weak boundary conditions easier. # # We now proceed to set up the problem in Firedrake using a square # mesh of quadrilaterals. :: from firedrake import * n = 100 mesh = UnitSquareMesh(n, n, quadrilateral=True) # We construct the quadratic function space for :math:u, :: V = FunctionSpace(mesh, "CG", 2) # and the function space for :math:\sigma. :: Sigma = TensorFunctionSpace(mesh, "CG", 2) # We then combine them together in a mixed function space. :: W = V*Sigma # Next, we set up the source function, which must integrate to the area # of the domain. Note how in the integration of the :class:~.Constant # one, we must explicitly specify the domain we wish to integrate over. :: x, y = SpatialCoordinate(mesh) fexpr = exp(-(cos(x)**2 + cos(y)**2)) f = Function(V).interpolate(fexpr) scaling = assemble(Constant(1, domain=mesh)*dx)/assemble(f*dx) f *= scaling assert abs(assemble(f*dx)-assemble(Constant(1, domain=mesh)*dx)) < 1.0e-8 # Now we build the UFL expression for the variational form. We will use # the nonlinear solve, so the form needs to be a 1-form that depends on # a Function, w. :: v, tau = TestFunctions(W) w = Function(W) u, sigma = split(w) n = FacetNormal(mesh) I = Identity(mesh.geometric_dimension()) L = inner(sigma, tau)*dx L += (inner(div(tau), grad(u))*dx - (tau[0, 1]*n[1]*u.dx(0) + tau[1, 0]*n[0]*u.dx(1))*ds) L -= (det(I + sigma) - f)*v*dx # We must specify the nullspace for the operator. First we define a constant # nullspace, :: V_basis = VectorSpaceBasis(constant=True) # then we use it to build a nullspace of the mixed function space :math:W. :: nullspace = MixedVectorSpaceBasis(W, [V_basis, W.sub(1)]) # Then we set up the variational problem. :: u_prob = NonlinearVariationalProblem(L, w) # We need to set quite a few solver options, so we'll put them into a # dictionary. :: sp_it = { # We'll only use stationary preconditioners in the Schur complement, so # we can get away with GMRES applied to the whole mixed system :: # "ksp_type": "gmres", # We set up a Schur preconditioner, which is of type "fieldsplit". We also # need to tell the preconditioner that we want to eliminate :math:\sigma, # which is field "1", to get an equation for :math:u, which is field "0". :: # "pc_type": "fieldsplit", "pc_fieldsplit_type": "schur", "pc_fieldsplit_0_fields": "1", "pc_fieldsplit_1_fields": "0", # The "selfp" option selects a diagonal approximation of the A00 block. :: # "pc_fieldsplit_schur_precondition": "selfp", # We just use ILU to approximate the inverse of A00, without a KSP solver, :: # "fieldsplit_0_pc_type": "ilu", "fieldsplit_0_ksp_type": "preonly", # and use GAMG to approximate the inverse of the Schur complement matrix. :: # "fieldsplit_1_ksp_type": "preonly", "fieldsplit_1_pc_type": "gamg", # Finally, we'd like to see some output to check things are working, and # to limit the KSP solver to 20 iterations. :: # "ksp_monitor": None, "ksp_max_it": 20, "snes_monitor": None } # We then put all of these options into the iterative solver, :: u_solv = NonlinearVariationalSolver(u_prob, nullspace=nullspace, solver_parameters=sp_it) # and output the solution to a file. :: u, sigma = w.split() u_solv.solve() File("u.pvd").write(u) # An image of the solution is shown below. # # .. figure:: ma.png # :align: center # # A python script version of this demo can be found here # __. # # .. rubric:: References # # .. bibliography:: demo_references.bib # :filter: docname in docnames