Re: [firedrake] firedrake Digest, Vol 53, Issue 4
Thank you for the swift response Prof. Farrell. Indeed, I sloppily read through the spaces and changing the internal bubble to degree 4 no longer gives an error. However the resulting system cannot be solved by PETSc. That is, changing the code to: P2 = FiniteElement("CG", "tetrahedron", 2) B4 = FiniteElement("B", "tetrahedron", 4) FB3 = FiniteElement("FacetBubble", "tetrahedron", 3) Vh = VectorFunctionSpace(mesh, P2 + B4 + FB3) Qh = FunctionSpace(mesh, "DG", 1) gives the following traceback: File "/home/firedrake/firedrake/src/firedrake/firedrake/adjoint/variational_solver.py", line 75, in wrapper out = solve(self, **kwargs) File "/home/firedrake/firedrake/src/firedrake/firedrake/variational_solver.py", line 281, in solve solving_utils.check_snes_convergence(self.snes) File "/home/firedrake/firedrake/src/firedrake/firedrake/solving_utils.py", line 140, in check_snes_convergence raise ConvergenceError(r"""Nonlinear solve failed to converge after %d nonlinear iterations. firedrake.exceptions.ConvergenceError: Nonlinear solve failed to converge after 3 nonlinear iterations. Reason: DIVERGED_DTOL I would expect that for larger problems this could be an issue, but for a small test mesh such as this (UnitCubeMesh(5,5,5)) it shouldn't be, unless of course I am making some grave error. I am attaching a complete log with the monitors in case if someone wants to take a look and maybe spot what could be going wrong. On Fri, 25 Jun 2021 at 06:00, <firedrake-request@imperial.ac.uk> wrote:
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. Mixed formulation for incompressible hyperelasticity (Bhavesh Shrimali) 2. Re: Mixed formulation for incompressible hyperelasticity (Bhavesh Shrimali) 3. Re: Mixed formulation for incompressible hyperelasticity (Patrick Farrell)
----------------------------------------------------------------------
Message: 1 Date: Thu, 24 Jun 2021 22:11:25 -0500 From: Bhavesh Shrimali <bhavesh.shrimali@gmail.com> To: firedrake@imperial.ac.uk Subject: [firedrake] Mixed formulation for incompressible hyperelasticity Message-ID: < CANw4EA2dOXvzvGtWvphFHBA6azcD9yunrhu_Zad96fdT8rrJGg@mail.gmail.com> Content-Type: text/plain; charset="utf-8"
Hi everyone,
I was trying to use a mixed-formulation for an incompressible hyperelasticity problem. See the following MWE which uses a Taylor-Hood approximation -- P2 for displacement (u) and P1 for pressure (p):
mesh = UnitCubeMesh(5,5,5) Vh = VectorFunctionSpace(mesh, "CG", 2) Qh = FunctionSpace(mesh, "CG", 1) Wh = Vho * Qh w = Function(Wh) dw = TrialFunction(Wh) v = TestFunction(Wh) u, p = split(w) mu, lmbda = Constant(1.), Constant(1.e3) Js = (lmbda + p + sqrt((lmbda + p)**2. + 4.*lmbda*mu ))/(2.*lmbda) F = Identity(mesh.geometric_dimension()) + grad(u) J = det(F) I1 = inner(F, F) Pi = p * (J - Js) + mu/2 * (I1 - 3) - mu * ln(Js) + lmbda/2 * (Js - 1)**2 Pi *= dx a = derivative(Pi, w, v) Jac = derivative(a, w, dw) bcs = [ DirichletBC(Wh.sub(0).sub(0), Constant(0.), 1), DirichletBC(Wh.sub(0).sub(1), Constant(0.), 3), DirichletBC(Wh.sub(0).sub(2), Constant(0.), 5), DirichletBC(Wh.sub(0).sub(0), Constant(1.), 2) ] problem = NonlinearVariationalProblem(a, w, bcs, J=Jac) solver = NonlinearVariationalSolver(problem) solver.solve() u, p = w.split() File("displacement.pvd").write(u)
Everything works fine as expected. Now, I wanted to try another pair of approximations that is inf-sup stable, namely the one in the celebrated paper by Crouzeix-Raviart (1973). See (4.36) in http://archive.numdam.org/article/M2AN_1973__7_3_33_0.pdf.
This amounts to approximating `u` by a "P2 + Bubble" and `p` by a "DG-1" space. This is also given in (3.2.3) in https://link.springer.com/content/pdf/10.1007%2F978-3-642-87047-7_5.pdf and (2.38) in https://link.springer.com/content/pdf/10.1007%2F978-3-642-61623-5_2.pdf
Now, in 3D as far as I can tell this amounds to approximating `u` with "P2 + an internal bubble + FacetBubble". Would anyone know how this could be achieved within firedrake? Because as far as I can see, firedrake supports `FacetBubble`. However, when I try something like:
B3 = FiniteElement("Bubble", "tetrahedron", 3) FB3 = FiniteElement("FacetBubble", "tetrahedron", 3) P2 = FiniteElement("CG", "tetrahedron", 2) Vh = FunctionSpace(mesh, P2 + B3 + FB3)
I get the following error
File "/usr/lib/python3.8/functools.py", line 875, in wrapper return dispatch(args[0].__class__)(*args, **kw) File "/home/firedrake/firedrake/src/tsfc/tsfc/finatinterface.py", line 189, in convert_enrichedelement elements, deps = zip(*[_create_element(elem, **kwargs) File "/home/firedrake/firedrake/src/tsfc/tsfc/finatinterface.py", line 189, in <listcomp> elements, deps = zip(*[_create_element(elem, **kwargs) File "/home/firedrake/firedrake/src/tsfc/tsfc/finatinterface.py", line 310, in _create_element finat_element, deps = convert(ufl_element, **kwargs) File "/usr/lib/python3.8/functools.py", line 875, in wrapper return dispatch(args[0].__class__)(*args, **kw) File "/home/firedrake/firedrake/src/tsfc/tsfc/finatinterface.py", line 177, in convert_finiteelement return lmbda(cell, element.degree()), set() File "/home/firedrake/firedrake/src/FInAT/finat/fiat_elements.py", line 318, in __init__ super(Bubble, self).__init__(FIAT.Bubble(cell, degree)) File "/home/firedrake/firedrake/src/fiat/FIAT/bubble.py", line 33, in __init__ super(Bubble, self).__init__(ref_el, degree, codim=0) File "/home/firedrake/firedrake/src/fiat/FIAT/bubble.py", line 24, in __init__ raise RuntimeError('Bubble element of degree %d and codimension %d has no dofs' % (degree, codim)) RuntimeError: Bubble element of degree 3 and codimension 0 has no dofs
Any pointers would be appreciated --
*Bhavesh Shrimali* *Ph.D. Student, **Civil and Environmental Engineering,*
*University of Illinois Urbana-Champaign*
*Champaign, IL, United States*
*Email : bshrima2@illinois.edu <bshrima2@illinois.edu>* *Web : **https://bhaveshshrimali.github.io/ <https://bhaveshshrimali.github.io/>*
Also, for definiteness, the complete code: from firedrake import * parameters["form_compiler"]["quadrature_degree"]=5 mesh = UnitCubeMesh(5,5,5) P2 = FiniteElement("CG", "tetrahedron", 2) B4 = FiniteElement("B", "tetrahedron", 4) FB3 = FiniteElement("FacetBubble", "tetrahedron", 3) Vh = VectorFunctionSpace(mesh, P2 + B4 + FB3) Qh = FunctionSpace(mesh, "DG", 1) Wh = Vh * Qh w = Function(Wh) dw = TrialFunction(Wh) v = TestFunction(Wh) u, p = split(w) mu, lmbda = Constant(1.), Constant(1.e3) Js = (lmbda + p + sqrt((lmbda + p)**2. + 4.*lmbda*mu ))/(2.*lmbda) F = Identity(mesh.geometric_dimension()) + grad(u) J = det(F) I1 = inner(F, F) Pi = p * (J - Js) + mu/2 * (I1 - 3) - mu * ln(Js) + lmbda/2 * (Js - 1)**2 Pi *= dx a = derivative(Pi, w, v) Jac = derivative(a, w, dw) bcs = [ DirichletBC(Wh.sub(0).sub(0), Constant(0.), 1), DirichletBC(Wh.sub(0).sub(1), Constant(0.), 3), DirichletBC(Wh.sub(0).sub(2), Constant(0.), 5), DirichletBC(Wh.sub(0).sub(0), Constant(0.5), 2) ] problem = NonlinearVariationalProblem(a, w, bcs, J=Jac) solver = NonlinearVariationalSolver(problem) solver_parameters={'snes_monitor': None, 'snes_view': None, 'ksp_monitor_true_residual': None, 'snes_converged_reason': None, 'ksp_converged_reason': None, 'snes_type': 'test', 'mat_type': 'aij'} solver.parameters.update(solver_parameters) solver.solve() u, p = w.split() up = project(u, VectorFunctionSpace(mesh, "CG", 2)) File("displacement.pvd").write(u) On Fri, 25 Jun 2021 at 08:04, Bhavesh Shrimali <bhavesh.shrimali@gmail.com> wrote:
Thank you for the swift response Prof. Farrell.
Indeed, I sloppily read through the spaces and changing the internal bubble to degree 4 no longer gives an error. However the resulting system cannot be solved by PETSc. That is, changing the code to:
P2 = FiniteElement("CG", "tetrahedron", 2) B4 = FiniteElement("B", "tetrahedron", 4) FB3 = FiniteElement("FacetBubble", "tetrahedron", 3) Vh = VectorFunctionSpace(mesh, P2 + B4 + FB3) Qh = FunctionSpace(mesh, "DG", 1)
gives the following traceback:
File "/home/firedrake/firedrake/src/firedrake/firedrake/adjoint/variational_solver.py", line 75, in wrapper out = solve(self, **kwargs) File "/home/firedrake/firedrake/src/firedrake/firedrake/variational_solver.py", line 281, in solve solving_utils.check_snes_convergence(self.snes) File "/home/firedrake/firedrake/src/firedrake/firedrake/solving_utils.py", line 140, in check_snes_convergence raise ConvergenceError(r"""Nonlinear solve failed to converge after %d nonlinear iterations. firedrake.exceptions.ConvergenceError: Nonlinear solve failed to converge after 3 nonlinear iterations. Reason: DIVERGED_DTOL
I would expect that for larger problems this could be an issue, but for a small test mesh such as this (UnitCubeMesh(5,5,5)) it shouldn't be, unless of course I am making some grave error.
I am attaching a complete log with the monitors in case if someone wants to take a look and maybe spot what could be going wrong.
On Fri, 25 Jun 2021 at 06:00, <firedrake-request@imperial.ac.uk> wrote:
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. Mixed formulation for incompressible hyperelasticity (Bhavesh Shrimali) 2. Re: Mixed formulation for incompressible hyperelasticity (Bhavesh Shrimali) 3. Re: Mixed formulation for incompressible hyperelasticity (Patrick Farrell)
----------------------------------------------------------------------
Message: 1 Date: Thu, 24 Jun 2021 22:11:25 -0500 From: Bhavesh Shrimali <bhavesh.shrimali@gmail.com> To: firedrake@imperial.ac.uk Subject: [firedrake] Mixed formulation for incompressible hyperelasticity Message-ID: < CANw4EA2dOXvzvGtWvphFHBA6azcD9yunrhu_Zad96fdT8rrJGg@mail.gmail.com> Content-Type: text/plain; charset="utf-8"
Hi everyone,
I was trying to use a mixed-formulation for an incompressible hyperelasticity problem. See the following MWE which uses a Taylor-Hood approximation -- P2 for displacement (u) and P1 for pressure (p):
mesh = UnitCubeMesh(5,5,5) Vh = VectorFunctionSpace(mesh, "CG", 2) Qh = FunctionSpace(mesh, "CG", 1) Wh = Vho * Qh w = Function(Wh) dw = TrialFunction(Wh) v = TestFunction(Wh) u, p = split(w) mu, lmbda = Constant(1.), Constant(1.e3) Js = (lmbda + p + sqrt((lmbda + p)**2. + 4.*lmbda*mu ))/(2.*lmbda) F = Identity(mesh.geometric_dimension()) + grad(u) J = det(F) I1 = inner(F, F) Pi = p * (J - Js) + mu/2 * (I1 - 3) - mu * ln(Js) + lmbda/2 * (Js - 1)**2 Pi *= dx a = derivative(Pi, w, v) Jac = derivative(a, w, dw) bcs = [ DirichletBC(Wh.sub(0).sub(0), Constant(0.), 1), DirichletBC(Wh.sub(0).sub(1), Constant(0.), 3), DirichletBC(Wh.sub(0).sub(2), Constant(0.), 5), DirichletBC(Wh.sub(0).sub(0), Constant(1.), 2) ] problem = NonlinearVariationalProblem(a, w, bcs, J=Jac) solver = NonlinearVariationalSolver(problem) solver.solve() u, p = w.split() File("displacement.pvd").write(u)
Everything works fine as expected. Now, I wanted to try another pair of approximations that is inf-sup stable, namely the one in the celebrated paper by Crouzeix-Raviart (1973). See (4.36) in http://archive.numdam.org/article/M2AN_1973__7_3_33_0.pdf.
This amounts to approximating `u` by a "P2 + Bubble" and `p` by a "DG-1" space. This is also given in (3.2.3) in https://link.springer.com/content/pdf/10.1007%2F978-3-642-87047-7_5.pdf and (2.38) in https://link.springer.com/content/pdf/10.1007%2F978-3-642-61623-5_2.pdf
Now, in 3D as far as I can tell this amounds to approximating `u` with "P2 + an internal bubble + FacetBubble". Would anyone know how this could be achieved within firedrake? Because as far as I can see, firedrake supports `FacetBubble`. However, when I try something like:
B3 = FiniteElement("Bubble", "tetrahedron", 3) FB3 = FiniteElement("FacetBubble", "tetrahedron", 3) P2 = FiniteElement("CG", "tetrahedron", 2) Vh = FunctionSpace(mesh, P2 + B3 + FB3)
I get the following error
File "/usr/lib/python3.8/functools.py", line 875, in wrapper return dispatch(args[0].__class__)(*args, **kw) File "/home/firedrake/firedrake/src/tsfc/tsfc/finatinterface.py", line 189, in convert_enrichedelement elements, deps = zip(*[_create_element(elem, **kwargs) File "/home/firedrake/firedrake/src/tsfc/tsfc/finatinterface.py", line 189, in <listcomp> elements, deps = zip(*[_create_element(elem, **kwargs) File "/home/firedrake/firedrake/src/tsfc/tsfc/finatinterface.py", line 310, in _create_element finat_element, deps = convert(ufl_element, **kwargs) File "/usr/lib/python3.8/functools.py", line 875, in wrapper return dispatch(args[0].__class__)(*args, **kw) File "/home/firedrake/firedrake/src/tsfc/tsfc/finatinterface.py", line 177, in convert_finiteelement return lmbda(cell, element.degree()), set() File "/home/firedrake/firedrake/src/FInAT/finat/fiat_elements.py", line 318, in __init__ super(Bubble, self).__init__(FIAT.Bubble(cell, degree)) File "/home/firedrake/firedrake/src/fiat/FIAT/bubble.py", line 33, in __init__ super(Bubble, self).__init__(ref_el, degree, codim=0) File "/home/firedrake/firedrake/src/fiat/FIAT/bubble.py", line 24, in __init__ raise RuntimeError('Bubble element of degree %d and codimension %d has no dofs' % (degree, codim)) RuntimeError: Bubble element of degree 3 and codimension 0 has no dofs
Any pointers would be appreciated --
*Bhavesh Shrimali* *Ph.D. Student, **Civil and Environmental Engineering,*
*University of Illinois Urbana-Champaign*
*Champaign, IL, United States*
*Email : bshrima2@illinois.edu <bshrima2@illinois.edu>* *Web : **https://bhaveshshrimali.github.io/ <https://bhaveshshrimali.github.io/>*
Hi Bhavesh, two things (see below).
On 25 Jun 2021, at 14:12, Bhavesh Shrimali <bhavesh.shrimali@gmail.com> wrote:
Also, for definiteness, the complete code:
from firedrake import * parameters["form_compiler"]["quadrature_degree"]=5 mesh = UnitCubeMesh(5,5,5) P2 = FiniteElement("CG", "tetrahedron", 2) B4 = FiniteElement("B", "tetrahedron", 4) FB3 = FiniteElement("FacetBubble", "tetrahedron", 3) Vh = VectorFunctionSpace(mesh, P2 + B4 + FB3) Qh = FunctionSpace(mesh, "DG", 1) Wh = Vh * Qh w = Function(Wh) dw = TrialFunction(Wh) v = TestFunction(Wh) u, p = split(w) mu, lmbda = Constant(1.), Constant(1.e3) Js = (lmbda + p + sqrt((lmbda + p)**2. + 4.*lmbda*mu ))/(2.*lmbda) F = Identity(mesh.geometric_dimension()) + grad(u) J = det(F) I1 = inner(F, F) Pi = p * (J - Js) + mu/2 * (I1 - 3) - mu * ln(Js) + lmbda/2 * (Js - 1)**2 Pi *= dx a = derivative(Pi, w, v) Jac = derivative(a, w, dw) bcs = [ DirichletBC(Wh.sub(0).sub(0), Constant(0.), 1), DirichletBC(Wh.sub(0).sub(1), Constant(0.), 3), DirichletBC(Wh.sub(0).sub(2), Constant(0.), 5), DirichletBC(Wh.sub(0).sub(0), Constant(0.5), 2) ] problem = NonlinearVariationalProblem(a, w, bcs, J=Jac) solver = NonlinearVariationalSolver(problem) solver_parameters={'snes_monitor': None, 'snes_view': None, 'ksp_monitor_true_residual': None, 'snes_converged_reason': None, 'ksp_converged_reason': None, 'snes_type': 'test', 'mat_type': 'aij'}
solver.parameters.update(solver_parameters)
You need to provide this dictionary of solver parameters when you construct the NonlinearVariationalSolver, rather than updating parameters after the fact. so: solver = NLVS(problem, solver_parameters=...) Second, it looks like the problem is that the nonlinear is not converging, even though the linear solve is fine. I suspect that you need to start with better initial guesses. Can you do continuation in a parameter (such as lambda) to ramp up the nonlinearity? You may also need to add a linesearch (e.g. 'snes_linesearch_type': 'bt') If I set backtracking linesearch and do: for val in [1, 10, 100, 1000]: lmbda.assign(val) solver.solve() I can converge a problem on a 4x4x4 mesh. Lawrence
Hi Bhavesh,
On 25 Jun 2021, at 14:29, Lawrence Mitchell <wence@gmx.li> wrote:
Hi Bhavesh,
two things (see below).
On 25 Jun 2021, at 14:12, Bhavesh Shrimali <bhavesh.shrimali@gmail.com> wrote:
Also, for definiteness, the complete code:
from firedrake import * parameters["form_compiler"]["quadrature_degree"]=5 mesh = UnitCubeMesh(5,5,5) P2 = FiniteElement("CG", "tetrahedron", 2) B4 = FiniteElement("B", "tetrahedron", 4) FB3 = FiniteElement("FacetBubble", "tetrahedron", 3) Vh = VectorFunctionSpace(mesh, P2 + B4 + FB3)
I realised one final subtlety here. You are gluing together elements like this, but this construction doesn't end up with an interpolatory element (so your non-zero boundary condition on the mesh face with marker 2 is wrong). You should make a NodalEnrichedElement(P2, B4, FB3) here. Lawrence
participants (2)
-
Bhavesh Shrimali
-
Lawrence Mitchell