Re: [firedrake] nullspace within nonlinearvariationalsolver?

27 Aug 2019

      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
-------------- next part --------------
HTML attachment scrubbed and removed
------------------------------
Message: 2
Date: Mon, 26 Aug 2019 18:58:00 -0400
From: Matthew Knepley <knepley@gmail.com>
To: Sander Rhebergen <srhebergen@uwaterloo.ca>
Cc: "firedrake@imperial.ac.uk" <firedrake@imperial.ac.uk>
Subject: Re: [firedrake] nullspace within nonlinearvariationalsolver?
Message-ID:
         <CAMYG4Gm1rLmrZX8otNz+5YB8BpV6djXzsU2HPDUshAKjYHxNmQ@mail.gmail.com>
Content-Type: text/plain; charset="utf-8"
On Mon, Aug 26, 2019 at 5:37 PM Sander Rhebergen <srhebergen@uwaterloo.ca>
wrote:
...
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
HI Sander,
Even if you set a nullspace here, you will get the same error because the
null space is not projected out, we just know
about it, thus LU will still fail. Getting the projected matrix would be
prohibitively expensive (I assume). Just to make sure,
the nullspace is constant pressure, right? I suggest either
a) Using a Schur complement scheme, sticking the pressure mass matrix as
the PC matrix, and then using LU
or
b) Using monolithic multigrid since you can use SVD on the coarse grid
with no problem
Or am I misunderstanding something. I have never used the VVP formulation.
Thanks,
Matt
...
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
_______________________________________________
firedrake mailing list
firedrake@imperial.ac.uk
https://mailman.ic.ac.uk/mailman/listinfo/firedrake
-- 
What most experimenters take for granted before they begin their
experiments is infinitely more interesting than any results to which their
experiments lead.
-- Norbert Wiener
https://www.cse.buffalo.edu/~knepley/ <http://www.cse.buffalo.edu/~knepley/>
-------------- next part --------------
HTML attachment scrubbed and removed
------------------------------
_______________________________________________
firedrake mailing list
firedrake@imperial.ac.uk
https://mailman.ic.ac.uk/mailman/listinfo/firedrake
End of firedrake Digest, Vol 32, Issue 9
****************************************
_______________________________________________
firedrake mailing list
firedrake@imperial.ac.uk
https://mailman.ic.ac.uk/mailman/listinfo/firedrake

Re: [firedrake] nullspace within nonlinearvariationalsolver?

Patrick Farrell