Re: [firedrake] nullspace within nonlinearvariationalsolver?

27 Aug 2019

      Hi Patrick,

Thanks for the response. I used the piece of code you just sent and tried it. I get the error:

Traceback (most recent call last):
  File "meevc_nullspace.py", line 116, in <module>
    bc = PressureFixBC(W.sub(2), 0, 1)
  File "meevc_nullspace.py", line 45, in __init__
    x = dm.getVecClosure(coordsSection, coordsVec, pt).reshape(-1, dim).mean(axis=0)
ValueError: cannot reshape array of size 3 into shape (2)

I think this is related to me using a PeriodicSquareMesh, because when I use a unit square mesh there are no error messages and the solver seems to run. Could you tell me what is happening in

x = dm.getVecClosure(coordsSection, coordsVec, pt).reshape(-1, dim).mean(axis=0)

and whether it can be modified for the PeriodicSquareMesh?

Thanks,

Sander

________________________________
From: Patrick Farrell <patrick.farrell@maths.ox.ac.uk>
Sent: Tuesday, August 27, 2019 2:57:03 PM
To: Sander Rhebergen
Cc: Firedrake Project
Subject: Re: [firedrake] nullspace within nonlinearvariationalsolver?

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()
                 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 <firedrake-bounces@imperial.ac.uk> on behalf of firedrake-request@imperial.ac.uk <firedrake-request@imperial.ac.uk>
*Sent:* Tuesday, August 27, 2019 7:00 AM
*To:* firedrake@imperial.ac.uk
*Subject:* firedrake Digest, Vol 32, Issue 9
Send firedrake mailing list submissions to
         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
You can reach the person managing the list at
         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>
To: "firedrake@imperial.ac.uk" <firedrake@imperial.ac.uk>
Subject: [firedrake] nullspace within nonlinearvariationalsolver?
Message-ID: <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