Re: [firedrake] nullspace within nonlinearvariationalsolver?

28 Aug 2019

      Thanks Matt, Patrick. That indeed took care of that error. There are some other errors now - the array of nodes that are fixed is empty so that I still have a zero pivot. Instead of fixing [0, 0] I fixed a point which I know lies in the mesh. This gives me a non-empty array of indices now (just the one point), but set.getDof(i)=0 for that point and so the nodes array is empty. I'll play around with the code a bit more, but for now I'll use Lawrence's suggestion.

________________________________
From: Matthew Knepley <knepley@gmail.com>
Sent: Tuesday, August 27, 2019 7:28 PM
To: Patrick Farrell
Cc: Sander Rhebergen; Firedrake Project
Subject: Re: [firedrake] nullspace within nonlinearvariationalsolver?

On Tue, Aug 27, 2019 at 2:57 PM Patrick Farrell <patrick.farrell@maths.ox.ac.uk<mailto:patrick.farrell@maths.ox.ac.uk>> wrote:
Hi Sander,

I find the easiest and most robust thing to do with these systems is just to pin
the pressure at a single point. That way you can use LU as you like, and can
focus on debugging the discretisation instead of the discretisation + solver
together.

It's not so straightforward to do this in firedrake. Here's some code I have
lying around that does this (probably written by Lawrence?):

         class PressureFixBC(DirichletBC):
             def __init__(self, V, val, subdomain, method="topological"):
                 super().__init__(V, val, subdomain, method)
                 sec = V.dm.getDefaultSection()
                 dm = V.mesh()._plex

                 coordsSection = dm.getCoordinateSection()
                 coordsDM = dm.getCoordinateDM()
                 dim = coordsDM.getDimension()

dim = dm.getCoordinateDim()

  Thanks,

    Matt

                 coordsVec = dm.getCoordinatesLocal()

                 (vStart, vEnd) = dm.getDepthStratum(0)
                 indices = []
                 for pt in range(vStart, vEnd):
                     x = dm.getVecClosure(coordsSection, coordsVec, pt).reshape(-1, dim).mean(axis=0)
                     if x.dot(x) == 0.0: # fix [0, 0] in original mesh coordinates (bottom left corner)
                         if dm.getLabelValue("pyop2_ghost", pt) == -1:
                             indices = [pt]
                         break

                 nodes = []
                 for i in indices:
                     if sec.getDof(i) > 0:
                         nodes.append(sec.getOffset(i))

                 self.nodes = numpy.asarray(nodes, dtype=IntType)

                 print("Fixing nodes %s" % self.nodes)

Use like bc = PressureFixBC(Z.sub(2), 0, 1). (The 1 is just a dummy label, it's
needed for the constructor but not used in the code.) If you're not on a unit
square, then you'd need to change the conditional x.dot(x) == 0.0 to something
else. (The reason why I search by coordinate is to pin the same value on every
level of a multigrid hierarchy.)

Hope this helps!

Patrick

On 27/08/2019 16:16, Sander Rhebergen wrote:
...
Hi Matt,
Thanks for the explanation. I'm new to this formulation as well and not sure about the best way to address the solver part, but I will try out your suggestions.
Sander
------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
*From:* firedrake-bounces@imperial.ac.uk<mailto:firedrake-bounces@imperial.ac.uk> <firedrake-bounces@imperial.ac.uk<mailto:firedrake-bounces@imperial.ac.uk>> on behalf of firedrake-request@imperial.ac.uk<mailto:firedrake-request@imperial.ac.uk> <firedrake-request@imperial.ac.uk<mailto:firedrake-request@imperial.ac.uk>>
*Sent:* Tuesday, August 27, 2019 7:00 AM
*To:* firedrake@imperial.ac.uk<mailto:firedrake@imperial.ac.uk>
*Subject:* firedrake Digest, Vol 32, Issue 9
Send firedrake mailing list submissions to
         firedrake@imperial.ac.uk<mailto:firedrake@imperial.ac.uk>
To subscribe or unsubscribe via the World Wide Web, visit
https://mailman.ic.ac.uk/mailman/listinfo/firedrake
or, via email, send a message with subject or body 'help' to
         firedrake-request@imperial.ac.uk<mailto:firedrake-request@imperial.ac.uk>
You can reach the person managing the list at
         firedrake-owner@imperial.ac.uk<mailto:firedrake-owner@imperial.ac.uk>
When replying, please edit your Subject line so it is more specific
than "Re: Contents of firedrake digest..."
Today's Topics:
1.  nullspace within nonlinearvariationalsolver? (Sander Rhebergen)
    2. Re:  nullspace within nonlinearvariationalsolver?
       (Matthew Knepley)
----------------------------------------------------------------------
Message: 1
Date: Mon, 26 Aug 2019 21:37:11 +0000
From: Sander Rhebergen <srhebergen@uwaterloo.ca<mailto:srhebergen@uwaterloo.ca>>
To: "firedrake@imperial.ac.uk<mailto:firedrake@imperial.ac.uk>" <firedrake@imperial.ac.uk<mailto:firedrake@imperial.ac.uk>>
Subject: [firedrake] nullspace within nonlinearvariationalsolver?
Message-ID: <d728cc7c78d84497a4a45da641289f9d@uwaterloo.ca<mailto:d728cc7c78d84497a4a45da641289f9d@uwaterloo.ca>>
Content-Type: text/plain; charset="iso-8859-1"
Hello,
Is it possible to set a constant nullspace for a NonlinearVariationalSolver?
I tried the following:
order = 2
V3 = FunctionSpace(mesh, "CG", order)
V1 = FunctionSpace(mesh, "BDM", order)
V2 = FunctionSpace(mesh, "DG", order-1)
W  = MixedFunctionSpace([V3, V1, V2])
v_basis = VectorSpaceBasis(constant=True)
null_space = MixedVectorSpaceBasis(W, [W.sub(0), W.sub(1), v_basis])
# Build Nonlinear Solver
uprob = NonlinearVariationalProblem(L, w1)
usolver = NonlinearVariationalSolver(uprob, nullspace=null_space, solver_parameters = {
     'snes_max_it': '50',
     'snes_rtol': '1e-12',
     'mat_type': 'aij',
     'ksp_type': 'preonly',
     'pc_type': 'lu',
     'snes_monitor': True,
     'snes_view': False,
     'snes_converged_reason': True,
     'ksp_converged_reason': True})
but the linear solver did not converge:
0 SNES Function norm 7.964422825545e-02
   Linear firedrake_1_ solve did not converge due to DIVERGED_PCSETUP_FAILED iterations 0
                  PCSETUP_FAILED due to FACTOR_NUMERIC_ZEROPIVOT
Nonlinear firedrake_1_ solve did not converge due to DIVERGED_LINEAR_SOLVE iterations 0
I can send the full code if necessary. I'm trying to solve Navier-Stokes in velocity-vorticity-pressure formulation (https://arxiv.org/pdf/1604.00257.pdf) but the solver doesn't converge. If I add a small parameter times a pressure mass-matrix, the solver converges, but I was hoping to avoid doing this.
Thanks,
Sander