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
>