Custom measure to refer to a part of the boundary

dolfinx
1

Hi all.

I’m working on a system of non-linear PDEs in a sphere. Accordingly, I have defined everything in spherical coordinates. Since the problem has rotational symmetry, things depend only on r and \theta and the mesh is simply a [0,R]\times[0,\pi] rectangle.

Following the Setting multiple Dirichlet, Neumann, and Robin conditions chapter of the FEniCSx tutorial, I define a custom measure

domain = mesh.create_rectangle(MPI.COMM_WORLD, np.array([[0,0],[R,np.pi]]),[10,int(10*np.pi)])
x = ufl.SpatialCoordinate(domain)

boundaries = [(1, lambda x: np.isclose(x[0], 0)), # origin
              (2, lambda x: np.isclose(x[0], R)), # surface of the sphere
              (3, lambda x: np.isclose(x[1], 0)),
              (4, lambda x: np.isclose(x[1], np.pi))]
facet_indices, facet_markers = [], []
fdim = domain.topology.dim - 1
for (marker, locator) in boundaries:
    facets = mesh.locate_entities(domain, fdim, locator)
    facet_indices.append(facets)
    facet_markers.append(np.full_like(facets, marker))
facet_indices = np.hstack(facet_indices).astype(np.int32)
facet_markers = np.hstack(facet_markers).astype(np.int32)
sorted_facets = np.argsort(facet_indices)
facet_tag = mesh.meshtags(domain, fdim, facet_indices[sorted_facets], facet_markers[sorted_facets])

ds = ufl.Measure("ds", domain=domain, subdomain_data=facet_tag)

so I can refer only to the r=R part of the boundary when writing the surface contribution to the energy:

F = fB(m,x)*Jacob(x)*ufl.dx + fS(m,m0)*Jacob(x)*ufl.ds(2)

I then compute the variation of this functional, which gives rise to Robin BCs in that region of the boundary, and solve the non-linear problem:

dF = ufl.derivative(F, m, phi) # variation of the energy
JF = ufl.derivative(dF, m, dm) # Jacobian of the energy
problem = fem.petsc.NonlinearProblem(dF, m, J=JF)

This procedure, however, doesn’t work:

Number of interations: 2
2022-08-09 14:55:33.806 (   0.110s) [main            ]              petsc.cpp:677   INFO| PETSc Krylov solver starting to solve system.
2022-08-09 14:55:33.839 (   0.143s) [main            ]              petsc.cpp:677   INFO| PETSc Krylov solver starting to solve system.
2022-08-09 14:55:33.846 (   0.151s) [main            ]       NewtonSolver.cpp:36    INFO| Newton iteration 2: r (abs) = 0 (tol = 1e-10) r (rel) = -nan(tol = 1e-06)
2022-08-09 14:55:33.847 (   0.151s) [main            ]       NewtonSolver.cpp:255   INFO| Newton solver finished in 2 iterations and 2 linear solver iterations.

Weirdly enough, if I define the energy by integrating over the complete boundary, F = fB(m,x)*Jacob(x)*ufl.dx + fS(m,m0)*Jacob(x)*ufl.ds, the solver converges well and the result is perfectly compatible with what is expected.

Am I defining the custom measure wrongly? Adding one more layer of complexity, my next step is to define different conditions for r=R,\,\theta<\pi/2 and r=R,\,\theta\geq\pi/2, so I will definitely need to separate the boundary in subregions.

2

It’s difficult to debug code when you don’t provide a MWE. But perhaps this is the issue: you’ve created a measure ds = ... but you haven’t used it, you’re using ufl.ds(2) where ufl.ds has no associated subdomain data.

1 Like
3

Thanks @nate! It was that. I got mixed between ufl.ds and ds.