Re: [firedrake] firedrake Digest, Vol 32, Issue 9
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
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
Just a small adjustment to this. If you fix the pressure at a node, then solve, you should then subtract off the mean value of the soln so it is mean value zero. Even if you fix it at a corner (like he suggests) there can be a slight degradation in the convergence (in L^2) by a log factor. I know you can compute the mean value fairly easily in FEniCS; I assume you can do it in FireDrake as well. -Shawn -----Original Message----- From: firedrake-bounces@imperial.ac.uk [mailto:firedrake-bounces@imperial.ac.uk] On Behalf Of Patrick Farrell Sent: Tuesday, August 27, 2019 1:57 PM To: Sander Rhebergen <srhebergen@uwaterloo.ca> Cc: Firedrake Project <firedrake@imperial.ac.uk> 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://nam04.safelinks.protection.outlook.com/?url=https%3A%2F%2Fmail man.ic.ac.uk%2Fmailman%2Flistinfo%2Ffiredrake&data=02%7C01%7Cwalke r%40lsu.edu%7Ced9d57c5b8f2498e2ade08d72b206055%7C2d4dad3f50ae47d983a09 ae2b1f466f8%7C0%7C1%7C637025290393221667&sdata=nnRf%2B50i5w7BNyfjx ct1HYU2upGirb1P3yiFNFmcbjI%3D&reserved=0 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://nam04.safelinks.protection.outlook.com/?url=https%3A%2F%2Farxiv.org%2Fpdf%2F1604.00257.pdf&data=02%7C01%7Cwalker%40lsu.edu%7Ced9d57c5b8f2498e2ade08d72b206055%7C2d4dad3f50ae47d983a09ae2b1f466f8%7C0%7C1%7C637025290393221667&sdata=8r4g6RlPG8rszcUvom2V0YTlPcuKfyCo0IV19i1UNqs%3D&reserved=0) 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://nam04.safelinks.protection.outlook.com/?url=https%3A%2F%2Farx iv.org%2Fpdf%2F1604.00257.pdf&data=02%7C01%7Cwalker%40lsu.edu%7Ce d9d57c5b8f2498e2ade08d72b206055%7C2d4dad3f50ae47d983a09ae2b1f466f8%7C0%7C1%7C637025290393226659&sdata=mOkeF9n822N7VTw2t%2FTcXztR6vjRo7Q2YXGujbUYMpo%3D&reserved=0) 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://nam04.safelinks.protection.outlook.com/?url=https%3A%2F%2Fmai lman.ic.ac.uk%2Fmailman%2Flistinfo%2Ffiredrake&data=02%7C01%7Cwal ker%40lsu.edu%7Ced9d57c5b8f2498e2ade08d72b206055%7C2d4dad3f50ae47d983 a09ae2b1f466f8%7C0%7C1%7C637025290393226659&sdata=gCmvdTbgiqWyFJJ jUeZnTm%2BGmse5v5GYxOcUXuDuNe4%3D&reserved=0
-- 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://nam04.safelinks.protection.outlook.com/?url=https:%2F%2Fwww.cs e.buffalo.edu%2F~knepley%2F&data=02%7C01%7Cwalker%40lsu.edu%7Ced9d 57c5b8f2498e2ade08d72b206055%7C2d4dad3f50ae47d983a09ae2b1f466f8%7C0%7C 1%7C637025290393226659&sdata=lqH%2FTxH3LKRJKrDXWsVsSRuifeJ35jjMa0P KyZkEsrQ%3D&reserved=0 <https://nam04.safelinks.protection.outlook.com/?url=http:%2F%2Fwww.cs e.buffalo.edu%2F~knepley%2F&data=02%7C01%7Cwalker%40lsu.edu%7Ced9d 57c5b8f2498e2ade08d72b206055%7C2d4dad3f50ae47d983a09ae2b1f466f8%7C0%7C 1%7C637025290393226659&sdata=sSxuFrs4nGAixpv0ZQTWF1rJ5HnicRYQ0%2Fm iHIGfGaY%3D&reserved=0> -------------- next part -------------- HTML attachment scrubbed and removed
------------------------------
_______________________________________________ firedrake mailing list firedrake@imperial.ac.uk https://nam04.safelinks.protection.outlook.com/?url=https%3A%2F%2Fmail man.ic.ac.uk%2Fmailman%2Flistinfo%2Ffiredrake&data=02%7C01%7Cwalke r%40lsu.edu%7Ced9d57c5b8f2498e2ade08d72b206055%7C2d4dad3f50ae47d983a09 ae2b1f466f8%7C0%7C1%7C637025290393226659&sdata=gCmvdTbgiqWyFJJjUeZ nTm%2BGmse5v5GYxOcUXuDuNe4%3D&reserved=0
End of firedrake Digest, Vol 32, Issue 9 ****************************************
_______________________________________________ firedrake mailing list firedrake@imperial.ac.uk https://nam04.safelinks.protection.outlook.com/?url=https%3A%2F%2Fmail man.ic.ac.uk%2Fmailman%2Flistinfo%2Ffiredrake&data=02%7C01%7Cwalke r%40lsu.edu%7Ced9d57c5b8f2498e2ade08d72b206055%7C2d4dad3f50ae47d983a09 ae2b1f466f8%7C0%7C1%7C637025290393231663&sdata=U0L4Y1PQWFm605UYKxb ot%2BkOHwIOKL8mHPiWV2%2BZq%2Fg%3D&reserved=0
_______________________________________________ firedrake mailing list firedrake@imperial.ac.uk https://nam04.safelinks.protection.outlook.com/?url=https%3A%2F%2Fmailman.ic.ac.uk%2Fmailman%2Flistinfo%2Ffiredrake&data=02%7C01%7Cwalker%40lsu.edu%7Ced9d57c5b8f2498e2ade08d72b206055%7C2d4dad3f50ae47d983a09ae2b1f466f8%7C0%7C1%7C637025290393231663&sdata=U0L4Y1PQWFm605UYKxbot%2BkOHwIOKL8mHPiWV2%2BZq%2Fg%3D&reserved=0
On 27/08/2019 16:45, Shawn W Walker wrote:
Just a small adjustment to this.
If you fix the pressure at a node, then solve, you should then subtract off the mean value of the soln so it is mean value zero. Even if you fix it at a corner (like he suggests) there can be a slight degradation in the convergence (in L^2) by a log factor. I know you can compute the mean value fairly easily in FEniCS; I assume you can do it in FireDrake as well.
Great point. You should also manually orthogonalise if you specify the nullspace to PETSc, as PETSc enforces orthogonality in the Euclidean inner product: https://github.com/firedrakeproject/firedrake/issues/1329 Patrick
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
Hello, I tried Lawrence's suggestion and that seems to do the trick here. Thanks! Sander ________________________________ From: Sander Rhebergen Sent: Tuesday, August 27, 2019 3:44:12 PM To: Patrick Farrell Cc: Firedrake Project Subject: Re: [firedrake] nullspace within nonlinearvariationalsolver? 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
On 27/08/2019 17:14, Sander Rhebergen wrote:
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?
I think this is because of the way firedrake constructs periodic meshes (by first constructing a 3D mesh and then warping it). It looks like PETSc's opinion of the coordinates are the unwarped 3D coordinates. I think if you set dim = 3 it should work. Maybe Matt can chime in with the right code to calculate dim. Patrick
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
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
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
On 28 Aug 2019, at 15:02, Matthew Knepley <knepley@gmail.com> wrote:
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.
For periodic meshes, the coordinate field we have is DG, i.e. with no dofs associated with vertices. Cheers, Lawrence
On Wed, Aug 28, 2019 at 9:08 AM Lawrence Mitchell <wencel@gmail.com> wrote:
On 28 Aug 2019, at 15:02, Matthew Knepley <knepley@gmail.com> wrote:
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.
For periodic meshes, the coordinate field we have is DG, i.e. with no dofs associated with vertices.
But we are checking the solution space, not coordinate space. Right? Thanks, Matt
Cheers, Lawrence
-- 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/>
Hi Sander,
On 27 Aug 2019, at 20:46, Sander Rhebergen <srhebergen@uwaterloo.ca> 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
[...] It is sometimes possible to hit mumps hard enough to handle rank deficient systems. Perhaps something like -pc_type lu -pc_factor_mat_solver_type mumps And then additionally -mat_mumps_icntl_24 1 # enable detection of null pivots -mat_mumps_icntl_25 0 # return a possible solution See Section 5.9 of the mumps user manual for more: http://mumps.enseeiht.fr/doc/userguide_5.2.1.pdf Cheers, Lawrence
participants (6)
-
Lawrence Mitchell
-
Lawrence Mitchell
-
Matthew Knepley
-
Patrick Farrell
-
Sander Rhebergen
-
Shawn W Walker