On Wed, Aug 28, 2019 at 8:48 AM Sander Rhebergen <srhebergen@uwaterloo.ca> wrote:
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.
If sec.getDof(p) == 0 for your point p, there are no dofs associate with that point for the space V. Should this be true in your problem? If not, are you passing the correct V? You could always use sec.view() to see what is going on. Thanks, Matt
------------------------------ *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> 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 < 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